A Stabilised Nodal Spectral Element Method for
Fully Nonlinear Water Waves
Part 2: Wavebody interaction
Abstract
We present a new stabilised and efficient highorder nodal spectral element method based on the Mixed Eulerian Lagrangian (MEL) method for generalpurpose simulation of fully nonlinear water waves and wavebody interactions. In this MEL formulation a standard Laplace formulation is used to handle arbitrary body shapes using unstructured – possibly hybrid – meshes consisting of highorder curvilinear isoparametric quadrilateral/triangular elements to represent the body surfaces and for the evolving free surface. Importantly, our numerical analysis highlights that a single top layer of quadrilaterals elements resolves temporal instabilities in the numerical MEL scheme that are known to be associated with mesh topology containing asymmetric element orderings. The ’surface variable only’ free surface formulation based on introducing a particlefollowing (Lagrangian) reference frame contains quartic nonlinear terms that require proper treatment by numerical discretisation due to the possibility of strong aliasing effects. We demonstrate how to stabilise this nonlinear MEL scheme using an efficient combination of (i) global projection without quadrature errors, (ii) mild nonlinear spectral filtering and (iii) remeshing techniques. Numerical experiments revisiting known benchmarks are presented, and highlights that modelling using a highorder spectral element method provides excellent accuracy in prediction of nonlinear and dispersive wave propagation, and of nonlinear waveinduced loads on fixed submerged and surfacepiercing bodies.
keywords:
Nonlinear and dispersive free surface waves, Wavebody interaction, Marine Hydrodynamics, Offshore Engineering, Mixed Eulerian Lagrangian, HighOrder Spectral Element Method, Unstructured mesh, Highorder discretisation.color \FXRegisterAuthorakapekAK \FXRegisterAuthordbadbDB
1 Introduction
In the last decades, numerical simulation of hydrodynamics for free surface flows and the design of marine structures has become an indispensable tool in engineering analysis. The numerical tools need to accurately and efficiently account for nonlinear wavewave and wavebody interactions in large and representative marine regions in order to give reliable estimates of environmental loads on offshore structures such as floating production systems, shipwave hydrodynamics, offshore wind turbine installations and wave energy converters.
The development of numerical models for timedomain simulation of fully nonlinear and dispersive water wave propagation have been a research topic since the 1960s and with real engineering applications since the 1970s BeckReed2001 (). Models based on fully nonlinear potential flow (FNPF) have been considered relative mature for about two decades, see the review papers Yeung1982 (); TsaiYue1996 () and the references therein. More recent research focus on improved modelling by incorporating more physics and finally enabling largescale simulation of nonlinear waves CavaleriEtAl2007 (); Ma10 (). In many applications, computational speed is far more important than the cost of hardware and therefore algorithms that provide speed and scalability of work effort is of key interest.
The Boundary Element Method (BEM) is a widely used method in applications for free surface flow Yeung1975 () and flow around complex bodies due to the ease in handling complex geometry KimEtAll1999 (); HarrisEtAl2014 (). However, scalability has been shown to favour Finite Element Methods WuTaylor1995 (); ShaoEtAl2014 (). The use of Finite Element Methods (FEM) for fully nonlinear potential flow has also received significant attention starting with the original work of Wu & Eatock Taylor (1994) WuTaylor1994 (). Studies of free surface solvers based on secondorder FEM are, e.g. WuTaylor1995 (); GreavesBorthwickTaylor1997 (); MaWuTaylor2001 (); MaWuTaylor2001b (); TurnbullEtAl2003 (); WangWu2011 (); Zienkiewicz2014 (). For wavewave interaction and wave propagation with no bodies, the transformed methods can be used to map the physical domain into a fixed computational domain CaiEtAl1998 (); WesthuisAndonowati1998 (); ClaussSteinhagen:1999tr (); TurnbullEtAl2003 () is maybe the most relevant approach due to the numerical efficiency. However, most numerical models for free surface flows use the Mixed Eulerian Lagrangian method (MEL) LonguetHigginsCokelet1976 () for updating the free surface variables, especially for solving wavebody interaction problems. The MEL approach requires remeshing of and development of techniques that improves the efficiency, e.g. QALEFEM MaYan2006 (), have been put forward. It is wellknown that highorder discretisation methods can give significant reduction in computational effort compared to use of conventional loworder methods, especially for longtime simulations KO72 (). Higherorder schemes that has been proposed for FNPF models include highorder Finite Difference Methods (FDM) EBL08 (); DucrozetEtAl2010 (); ShaoEtAl2014 (), highorder BEM WangYaoTulin1995 (); Ferrant1996 (); CelebiEtAl1998 (); XueEtAl2001 (); FerrantEtAl2003 (); MaYan2009 (); YanLiu2011 (); HarrisEtAl2014b (), the HighOrder Spectral (HOS) Method DOMYOU1987 (); WestEtAl1987 (); GouinEtAl2017 () and other pseudospectral methods ChernEtA1999 (); ChristiansenEtAl2012 () based on global basis functions in a single element (domain). Indeed, it has been shown that efficient algorithms LiFleming1997 (); KorsmeyerEtal1999 (); EngsigKarup2014 () together with other means of acceleration (such as multidomain approaches and massively parallel computing EngsigKarupEtAl2011 (); NimmalaEtAl2013 () via software implementations on modern manycore hardware EGNL13 ()) render FNPF models practically feasible for analysis of wave propagation on standard work stations  even with a realtime perspective within reach for applications with appreciable domain sizes as discussed in EngsigKarupEtAl2011 ().
Even though highorder FDM have shown to be very efficient, secondorder FEM models remain popular because of the geometrical flexibility and the sparse matrix patterns of the discretisation. In contrast, a highorder finite element method, such as the Spectral Element Method (SEM) due to Patera (1984) PAT84 (), has historically received the least attention for FNPF equations. However, SEM is thought to be highly attractive as it combines the highorder accuracy of spectral methods for problems with sufficiently smooth solutions with the geometric flexibility via adaptive meshing capability of finite element methods. A previous study employing SEM for the FNPF model include RobertsonSherwin1999 () where an Arbitrary Lagrangian Eulerian (ALE) technique was used to track the freesurface motion identified that numerical instabilities were caused by asymmetry in the mesh. In order to stabilise the model, they added a diffusive term proportional to the mesh skewness to the kinematic free surface condition. Recently, a SEM model allowing for the computation of very steep waves has been proposed EngsigKarupEtAl2016 (). The use of transformed domains partitioned into a single vertical layer of elements is shown to avoid the fundamental problem of instabilities caused by mesh asymmetry. It was also found that the quartic nonlinear terms present in the Zakharov form ZAK1968 () of the free surface conditions could cause severe aliasing problems and consequently numerical instability for marginally resolved or very steep nonlinear waves. This problem was mitigated through overintegration of the free surface equations and application of a gentle spectral filtering using a cap of 1% of the highest modal coefficient. Finally, in the context of SEM for FNPF models we also mention the HOSE method ZhuEtAl2003 (), however, it is only the boundary domains that are discretised with SEM, not the interior domain.
Thus, the challenge we seek to address in this work is the development of a robust SEM based FNPF solver with support for geometric flexibility for handling wavebody interaction problems. This is achieved by extending the results of EngsigKarupEtAl2016 (), keeping the same dealiasing techniques. In order to include arbitrary shaped bodies, we discard the transform in favour for the MEL method. To avoid the temporal stability issues associated with asymmetric mesh configurations near the free surface boundary Westhuis2001 (); RobertsonSherwin1999 () we propose to use hybrid meshes consisting of a layer of quadrilaterals at the surface and unstructured triangles below. Eigenvalue analysis of the semidiscrete formulation for smallamplitude waves is used to illustrate that this mesh configuration is linearly stable and during the simulation the vertical interfaces of the quadrilaterals are kept vertical. Global remeshing is employed where the vertices of the mesh topology are repositioned similar to the technique used in the QALEFEM method MaYan2006 () to improve general stability of the model for nonlinear waves. This is combined with local remeshing of the node distributions of the free surface elements to keep the discrete operators wellconditioned for all times, and is shown to increase the temporal stability significantly. In the following, several test cases are used to show the robustness and accuracy of the proposed SEM model based on MEL.
1.1 On highorder methods for fully nonlinear potential flow models
Highorder methods allow for convergence rates faster than quadratic with mesh refinement, and is attractive for improving efficiency of numerical methods  especially for longtime integration KO72 (). Indeed, highorder discretisation methods are needed for realistic largescale applications. To understand why, we need to first to understand that the key criterion for deciding between practical numerical tools is based on identifying the tool which is the most numerically efficient one, e.g. measured in terms of CPU time. In general, we can express the minimum work effort for a numerical scheme in terms of the discretization parameter as
where is the spatial dimension of the problem and a characteristic mesh size. At the same time, for a numerical method the spatial error is bound by
where characterise the order of accuracy (rate of convergence). For SEM the optimal order is with the order of the local polynomial expansions. So, to be efficient, we desire that the error can be reduced much faster than the work increases, i.e. that faster than when resolution is increased (). Thus,
indicating that for larger dimensions, it is necessary that a higher order spatial discretisation is used or that the work is balanced or minimised, e.g. through adapting the mesh resolution. So in three spatial dimensions () it is not possible to be efficient for largescale applications with loworder methods since the convergence rate is for a loworder method. Another key requirement is the development of methods which leads to scalable work effort, i.e. linear scaling with the problem size when subject to mesh refinement pursuing better accuracy. We see that this is only possible if the order of accuracy at least matches the dimension of the problem, i.e. .
1.2 Paper contribution
The key challenge we seek to address in this work is the development of a FNPF solver with highorder convergence rate that also has support for geometric flexibility for handling arbitrary body shapes and other structures. To this end we exploit and extend the results of the recent work of EngsigKarup, Eskilsson & Bigoni (2016) EngsigKarupEtAl2016 () on a stabilised transformed spectral element method for efficient and accurate marine hydrodynamics applications. We consider in this work a Mixed Eulerian Lagrangian (MEL) method LonguetHigginsCokelet1976 () based on explicit timestepping, we discard the transformation to be able to include arbitrary shaped bodies in the fluid domain, and propose a remedy to the temporal stability issues associated with the use of asymmetric mesh configurations near the free surface boundary as identified Westhuis2001 () for a classical loworder Galerkin finite element method and RobertsonSherwin1999 () using a highorder spectral element method for a fully nonlinear potential flow solver. These new developments constitutes a new robust SEM methodology for nonlinear wavebody interactions.
2 Mixed EulerianLagrangian (MEL) formulation
The governing equations for Fully Nonlinear Potential Flow (FNPF) is given in the following. Let the fluid domain () be a bounded, connected domain with piecewise smooth boundary and introduce restrictions to the free surface and the bathymetry . Let be the time domain.
We seek a scalar velocity potential function satisfying the Laplace problem
(1a)  
(1b)  
(1c) 
where describes variation in the still water depth. The evolution of the free surface boundary is described by . The notations are illustrated in Figure 1.
The MEL time marching technique assumes a particlefollowing (Lagrangian) reference frame for the free surface particles with changes in position in time given by
(2) 
where having assumed only two space dimensions.
Thus, at the free surface boundary this implies kinematic conditions of the form
(3a)  
that must be satisfied together with the dynamic boundary condition stated in terms of Bernoulli’s equation  
(3b) 
where the material derivative connects the Eulerian and Lagrangian reference frame through the relation
(4) 
where is a velocity vector for the moving frame of reference that can be chosen arbitrarily. If we obtain an Eulerian frame of reference and if we have a Lagrangian frame of reference.
Spatial and temporal differentiation of free surface variables are given by the chain rules
(5) 
and is useful for expressing the free surface equations valid at in terms of free surface variables only, cf. EGNL13 (), in the form
(6a)  
(6b)  
(6c) 
having assumed a zero reference pressure at the free surface. The tilde ’’ is used to denote that a variable is evaluated at the free surface, e.g. . Note, this formulation contains quartic nonlinear terms as in the Zakharov form ZAK1968 () for the Eulerian formulation of FNPF equations. These nonlinear terms needs proper treatment to deal with aliasing effects.
2.1 Boundary conditions
For the solution of the Laplace problem, the following free surface boundary condition is specified
(7) 
while at fixed vertical boundaries, impermeable wall boundary conditions are assumed
(8) 
where denotes an outward pointing unit normal vector to . At rigid vertical wall boundaries the domain boundary conditions at free surface variables imposed are
(9) 
Wave generation and absorption zones are included using the embedded penalty forcing technique described in EGNL13 ().
3 Numerical discretisation
Following EngsigKarupEtAl2016 (), we present the discretization of the governing equations in a general computational framework based on the method of lines, where first a semidiscrete system of ordinary differential equations is formed by spatial discretisation in two space dimensions using a nodal SEM. The semidiscrete system is subject to temporal integration performed using an explicit fourthorder RungeKutta method.
3.1 Weak Galerkin formulation and discretisation
We form a partition of the domain to obtain a tessellation of consisting of nonoverlapping shaperegular elements such that with denoting the ’th element. We introduce for any tessellation the spectral element approximation space of continuous, piecewise polynomial functions of degree at most ,
which is used to form finitedimensional nodal spectral element approximations
(10) 
where is the global finite element basis functions with cardinal property at mesh nodes with the Kronecker Symbol.
3.1.1 Unsteady free surface equations
The weak formulation of the free surface equations (6) is derived in the following form. Find where such that
(11a)  
(11b)  
(11c) 
for all . Substitute the expressions in (14) into (11) and choose . The discretisation in one spatial dimension becomes
(12a)  
(12b)  
(12c) 
where the following global matrices have been introduced
(13) 
where is a vector containing the set of discrete nodal values.
Following EngsigKarupEtAl2016 (), the gradients of the free surface state variables are recovered via a global gradient recovery technique based on global Galerkin projections that work for arbitrary unstructured meshes in the SEM framework as described in EngsigKarupEtAl2016 (). Aliasing effects are effectively handled using exact quadrature for nonlinear terms combined with a mild spectral filtering technique HW08 () that gently removes highfrequency noise that may arise as a result of marginal resolution.
Remark: The free surface node positions are changing in time, and this implies that the mesh must change accordingly. Thus, the scheme needs to recompute the global spectral element matrices (13) at every time step which impacts the computational efficiency of the scheme.
3.1.2 Curvilinear isoparametric elements
To handle arbitrary body shapes we partition the fluid domain to obtain another tessellation consisting of nonoverlapping shaperegular elements such that the tessellation can be formed by combining quadrilateral and triangular curvilinear elements into an unstructured hybrid mesh such that . In two space dimensions, the nodal spectral element approximations takes the form
(14) 
where is the total degrees of freedom in the discretisation.
The curvilinear elements makes it possible to treat the deformations in the free surface and the body surfaces as illustrated in Figure 2 with two different hybrid unstructured meshes. In both cases a single quadrilateral layer is used just below the free surface.
Consider the ’th element . On this element, we form a local polynomial expansion expressed as
(15) 
where both a modal and a nodal expansion in the reference element is given in terms of nodes/modes. We have introduced a map to take nodes from the physical element to a reference element, where is a single computational reference element. This dual representation can be exploited to form the curvilinear representations, since the coefficient vectors are related through
(16) 
where , is the set of orthonormal basis functions and with nodes in the reference element that defines the Lagrangian basis. For quadrilaterals a tensor product grid formed by LegendreGaussLobatto nodes in 1D is used. For triangles, the node distribution is determined using the explicit warp & blend procedure War06 (). The onetoone mapping from a general curvilinear element to the reference element is highlighted in Figure 5. The edges of the physical quadrilateral elements are defined by the functions , , and by introducing isoparametric polynomial interpolants Zienkewitch1971 () of the same order as the spectral approximations of the form
(17) 
it is possible to represent curved boundaries.
Using transfinite interpolation with linear blending GH73 () the affine transformation from the square reference quadrilateral to the physical domain is defined in terms of these edge curves in the form
(18) 
Similarly, using transfinite interpolation with linear blending, the transformation for triangles is given as
(19) 
The use of curvilinear elements, implies that the Jacobian of the mapping is no longer constant as for straightsided triangular elements. Thus, to avoid quadrature errors higher order quadratures are employed in the discrete Galerkin projections and this increases the cost of the scheme proportionally for the elements in question. Note, the introduction of curvilinear elements, can be setup in a preprocessing step and therefore does not add any significant additional complexity to the numerical scheme.
3.1.3 Spatial discretisation of the Laplace problem
Consider the discretisation of the governing equations for the Laplace problem (1). We seek to construct a linear system of the form
(20) 
The starting point is a weak Galerkin formulation that can be expressed as: find such that
(21) 
for all where the boundary integrals vanish at domain boundaries where impermeable walls are assumed. The discrete system operator is defined as
(22) 
The elemental integrals are approximated through change of variables as
(23) 
where is the Jacobian of the affine mapping . The global assembly of this operator preserves the symmetry, and the resulting linear system is modified to impose the Dirichlet boundary conditions (7) at the free surface. The vertical free surface velocity is recovered from the potential via a global Galerkin projection that involves a global matrix for the vertical derivative.
4 Numerical properties
We start out by considering the numerical properties of the model related to the temporal stability and convergence of the numerical MEL scheme. Results of comparison with the stabilised Eulerian formulation EngsigKarupEtAl2016 () is included since the two schemes are complementary.
4.1 Temporal linear stability analysis of semidiscrete system
We revisit the analysis of temporal instability following the works of Westhuis2001 (); RobertsonSherwin1999 () by considering the semidiscrete free surface formulation that arise under the assumption of smallamplitude waves
(24a)  
(24b)  
(24c) 
The discretization of this semidiscrete systems leads to
(25) 
where and
(26) 
where
(27) 
and having introduced matrix decompositions of the global matrices of the form
(28) 
where the subscript indices ’b’ refers to the free surface nodes and the ’i’ refers to all interior nodes. The eigenspectrum of determines the temporal stability of this system.
In the context of SEM, we seek to first confirm the results given in RobertsonSherwin1999 (). By doing a similar eigenanalysis using a triangulated asymmetric mesh, we confirm that we have temporal instability, e.g. see representative results in Figure 11. To fix this problem, we need to avoid asymmetric meshes near the free surface layer. Following RobertsonSherwin1999 () we can also consider a triangulated symmetric mesh in Figure (b)b which turns out to be stable for all polynomial orders except 3 and 4 used in the analysis. This is in line with the results and conclusions presented by Robertson & Sherwin (1999), but reveals that numerical issues may arise for such triangulated meshes in specific configurations. Furthermore, since our objective is to introduce arbitrary shaped bodies inside the fluid domain, in the next experiment we use a hybrid mesh. This mesh is consisting of a layer of quadrilateral below the free surface of the fluid similar to the meshes used in EngsigKarupEtAl2016 () and combine this layer with an additional triangulated layer used to represent complex geometries such as a submerged cylinder as illustrated in Figure LABEL:fig:AC_18. The eigenanalysis shows that we then reach temporal stability for arbitrary polynomial expansion orders. These results are in line with other experiments we have carried out that are not presented here, confirming that by introducing a quadrilateral layer we can fix the temporal stability problem to obtain purely imaginary eigenspectra (to machine precision). Using instead a similar mesh but with slightly skewed quadrilaterals, reveals again that the mesh asymmetry leads to temporal instability as shown in Figure 17. So, the temporal instability is associated with the accuracy of the vertical gradient approximation that is used to compute at the free surface and that determines the dispersive properties of the model. So poor accuracy in destroys the general applicability of the model since the wave propagation cannot be resolved accurately. This makes it clear that a quadrilateral layer with vertical alignment of nodes close to the free surface provides accurate recovery of the vertical free surface velocities and fixes the temporal instability problem.




Remark: This new insight paves the road towards also stabilising nonlinear simulations using a highorder accurate unstructured method such as SEM following the recipe for stabilisations recently laid out in EngsigKarupEtAl2016 () for an Eulerian scheme based on quadrilaterals only.
In the following, we consider the strenuous benchmark of nonlinear wave propagation of steep nonlinear stream function waves by setting up a mesh with periodicity conditions imposed to connect the eastern and western boundaries. Then, we examine the conservation of energy and accuracy to assess the stability of the model. This analysis shows that the exact integration used in the discrete global projections has a big effect on stabilising the solution. Following EngsigKarupEtAl2016 (), a test is performed for different amounts of filtering in the last mode and for different maximum steepness (). For example, the analysis in Figure 22 demonstrates that a 1% filter works well for the MEL Formulation.
To maintain temporal stability when using explicit timestepping methods, it is necessary to employ remeshing DOMYOU1987 (). The objective is to either prevent strict CFL restrictions although this is not severe for FNPF methods where numerical stability is influenced only weakly by spatial resolution in the horizontal EGNL13 (); EngsigKarupEtAl2016 (), and to retain a proper mesh quality to avoid temporal instabilities induced by the mesh or illconditioning in the global operators (discussed in BjontegaardRonquist2009 ()). Steep nonlinear waves are associated with larger node movements producing deformations in intranode positions in the free surface elements. These deformations changes the conditioning of Vandermonde matrices and may impact the accuracy of the numerical scheme. This problem can be prevented by including a local remesh operation which changes the interior node distribution while not changing the initial mesh topology. For example, whenever an element goes below 75% or above 125% of its original size and with the initial free surface elements assumed of uniform size. The idea is to stop tracking the original material nodes at the free surface, and reposition the intraelement nodes through a local operation via interpolation to the original node distribution used in each element based on LegendreGaussLobatto nodes. A simple local remesh operation of this kind helps to stabilise the solver by improving accuracy in the numerical solutions and is a technique used in similar FEM solvers, e.g. see WuTaylor1994 (). The effects of local remeshing on temporal stability for a steep nonlinear stream function wave is presented in Figure 22, and shows that local remeshing is essential to maintain temporal stability for longer integration times for steep nonlinear waves. A snapshot of a stream function wave after propagating for wave periods the wave is still represented very accurately with essentially no amplitude or dispersion errors, cf. Figure 22 (a). Furthermore, both mass and energy is conserved with high accuracy, cf. Figure 22 (b).




Also, we compare the MEL method and the Eulerian method due to EngsigKarupEtAl2016 () in terms of different stabilisation techniques in Figure 25. We find that with our current strategies the stabilised MEL method is more robust than the stabilised Eulerian scheme.


4.2 Convergence tests
To validate the highorder spectral element method, we demonstrate for the MEL scheme that the highorder convergence rate is in line with EngsigKarupEtAl2016 () and the accuracy of the method using convergence tests as depicted in Figure 26. These tests, demonstrate that we can exert control over approximations errors in the scheme by adjusting the resolution in terms of choosing the points per wave length. For the most nonlinear waves of 90% of maximum steepness, the curves find a plateau at the size of truncation error before convergence to machine precision due to insufficient accuracy in our numerical solution of stream function waves used.
5 Numerical experiments
We examine different strenuous test cases that serve as validation of the numerical spectral element model proposed.
5.1 Reflection of highamplitude solitary waves from a vertical wall
We setup a solitary wave as initial condition using the highorder accurate spectral numerical scheme due to DutCla14 () and consider the propagation of solitary waves of different amplitudes above a flat bed that are reflected by a vertical solid wall. In this experiment, a solitary wave approaches the wall with constant speed and starts to accelerate forward when the crest is at a distance of approximately from the wall. The water level at the wall position grows leading to the formation of a thin jet shooting up along the wall surface. When the maximum of free surface elevation is reached at the wall position it is said that the wave is attached to the wall, at this moment and the height of the crest is . Thereafter, the jet forms and reaches its maximum run up at time . After this event, the jet collapses slightly faster than it developed. The detachment time corresponds to the wave crest leaving the wall. The height of the wave at that instant is , always smaller than the attachment height (i.e. ). Then a reflected wave propagates in the opposite direction with the characteristics of a solitary wave of reduced amplitude. This adjustment in the height produces a dispersive trail behind the wave. The depression becomes more abrupt for increasing wave steepness.
This test case requires the fully nonlinear free surface boundary conditions to capture the steepest solitary waves and their nonlinear interaction with the wall. The case was first studied using a numerical method by Cooker et al. Cooker97 (), where a BoundaryIntegral Method (BIM) was used to solve the Euler equations with fully nonlinear boundary conditions. The results obtained using BIM showed to be in excellent agreement with the experimental data given in Max76 (), however, showed difficulties in the down run phase for the steepest wave of indicating that for the most nonlinear cases, numerical modelling is difficult. Accurate simulation of the steepest waves requires sufficient spatial resolution and high accuracy in the kinematics to capture the finest details of the changing kinematics of the fluid near the wall. This problem was recently been addressed using other highorder numerical models, e.g. a highorder Boussinesq model with same fully nonlinear free surface boundary conditions based on a highorder FDM model MBL02 () and a nodal discontinuous Galerkin spectral element method EHBM06 () and both studies showed excellent results for the wave runup and depthintegrated force histories up to .
The numerical experiments are carried out using a domain size m, with the initial position of the center of the solitary wave at m. We employ a structured mesh consisting of quadrilaterals with elements in a single layer, where is varied proportional to the wave height to resolve the waves in the range . Each element is based on polynomial expansion orders . The time step size is chosen in the interval s. Stabilisation of the numerical scheme is achieved using exact integration and mild spectral filtering using a top mode spectral filter every time step using the MEL scheme.




The attachment analysis of Cooker97 () is reproduced using SEM with the initial profile of the solitary wave DutCla14 () up to with excellent agreement, cf. Figures (a)a and (b)b. This time history of the wave propagation is illustrated using the MEL scheme in Figure 31 for a highamplitude wave with relative amplitude . The forces are computed numerically as a postprocessing step by integration of the local pressure obtained from Bernoulli’s equation
(29) 
along the wet part of the wall. Approximation of the time derivative is approximated by using the acceleration potential method due to Tanizawa1995 (). This method estimates by solving an additional boundary value problem based on a Laplace problem defined in terms of in the form
(30) 
subject to the boundary conditions
(31a)  
(31b) 
The wall force vector is determined numerically from pressure using
(32) 


5.2 Rectangular obstacle piercing the free surface
In Lin (2006) Lin06 () a numerical method for NavierStokes equations is proposed based on transformation of the fluid domain using a multiplelayer coordinate model. Wavestructure interactions of a solitary wave of amplitude (mildly nonlinear) with a rectangular obstacle fixed at different positions (seated, midsubmergence, floating) are investigated in a two dimensional wave flume of depth m. Free surface elevation histories at three different gauge locations are found in excellent agreement with computed VOF solutions. We consider here the experiment corresponding to the obstacle piercing the free surface (width: m, height: m, draught: m) for validation of the present methodology dealing with semisubmerged bodies.
A wave flume of m is sketched in Figure (a)a with the cylinder located at its center ( m). In this scenario, positions of gauges 1 and 3 in Lin06 () correspond to m and m, respectively. This domain is discretised using hybrid meshes with a top layer of curvilinear quadrilaterals. Below, unstructured grid of triangles is employed in the central part ( m), fitting the body contour (Triangles size is chosen proportional to their distance to the lower corners of the body). The zones away from the structure m m are meshed with regular quadrilaterals of width m, small enough to guarantee accurate propagation of this wave according to our tests in Section 5.1. A detail of this mesh corresponding to m is shown in Figure (b)b.


Similar to previous experiments with solitary waves, the initial condition is generated using the accurate numerical solution due to DutCla14 () corresponding to a wave peaked at m. Using the MEL scheme, the time step size is s. The evolution of the free surface elevations at gauge positions are compared with the Lin (2006) results in Figure 40, and excellent agreement is observed.


5.3 Solitary Wave Propagation Over a Submerged Semicircular Cylinder


We consider the numerical experiment described in WangEtAl2013 () on the interaction between a solitary wave of amplitude (mildly nonlinear) and a submerged semicircular cylinder with radius . This experiment serves to validate our force estimation, where the acceleration method due to Tanizawa (1995) is used. The initial mesh is illustrated in Figure 43 and consists of 60 elements in the central part and only 164 elements in total. The computed time evolution of the free surface is given in Figure 48 (a). In Figure 48 (d) it is seen how the solitary wave propagates undisturbed from the initial condition until it starts the interaction with the semicircular cylinder resulting in variation in the horizontal dynamic load. During the interaction, a dispersive trail develops after which the solitary wave restore to it’s original form. The results are in good agreement with results of Wang et al. The differences in the force curves given in Figure 48 (c) is understood to be related to numerical dispersion in the Desingularized Boundary Integral Equation Method (DBIEM) and approximation errors associated with the loworder accurate initial solitary wave. The force curve produced using the SEM solver is for a wellresolved flow and therefore is considered converged to high accuracy. In same figure, the force curve for Morison Equation that takes into account also viscous effects is included and is computed following KlettnerEames2012 (); WangEtAl2013 (). In figure (b) it is confirmed that the SEM conserves with high accuracy mass and energy during all of the simulation.




5.4 Solitary wave propagation over a fixed submerged cylinder
Several studies with both experimental and numerical results can be found regarding the interactions of solitary waves with a submerged circular cylinder Sib_etal82 (); Cooker_etal90 (); ChianErtekin1992 (); CleMas95 (). In the following, we consider two previously reported experiments.
First, we consider the introduction of a fixed submerged cylinder inside the fluid domain with a flat bed following the experiments due to Clement & Mas (1995) CleMas95 (). Analysis of the hydrodynamic horizontal forces exerted on the cylinder and comparison with experimental results Sib_etal82 () serve as validation in the experiments reported here.
Figure 51 shows a typical initial mesh for the experiments. The mesh consist of a small quadrilateral layer near the free surface and contains zones of quadrilateral elements only in first and last part of domain. These zones flank a central part, where a hybrid meshing strategy is employed to adapt using an unstructured triangulation to the submerged cylinder surface in the interior. The initial conditions , for the solitary waves is produced using the method described in DutCla14 ().
In our first experiment, the solitary wave height is and corresponds to zone 1 described in CleMas95 () where interaction is said weak. The cylinder has dimensionless radius and is positioned with submergence . Figure (a)a shows the resulting time series of the free surface, where we observe a small deformation of the wave during passing the cylinder. This deformation becomes wider at the end and more abrupt in positions just above the cylinder (centered at ). The computed mass and energy conservation measures given in Figure (b)b confirms stability and accuracy of the simulation. The dimensionless horizontal component is compared in Figure LABEL:fig:HC_15a with experimental results of Sib_etal82 () given in (CleMas95, , Figure 2b). The results are in good agreement for times before the collision (when horizontal force is zero) and significant discrepancies afterwards where viscous effects are important due to vortex shedding in the experiments Sib_etal82 (). These effects cannot be captured by potential flow models like the present one.






In our second experiment, the wave height is and the cylinder has radius and is positioned at a submergence . Following CleMas95 () this experiment is located in the crestcrest exchange zone 2. In this zone incident and transmitted waves can be clearly distinguished. The fluid motion in this experiment is studied when the wave passes the cylinder. In Figure (a)a and Figure (b)b the horizontal and vertical velocity components in a portion of the fluid domain () are illustrated at simultaneous times highlighting the changes in the subsurface kinematics during the solitary wave interaction with the cylinder. While approaching the position of the cylinder the wave looses height and slows down, cf. Figure (b)b and Figure (c)c. Then the birth of a transmitted wave can be observed, cf. Figure (b)b. The transmitted wave is not shifted and leaves behind a dispersive trail that becomes even more abrupt than observed in compared to wave interaction in our first experiment, cf. Figure (a)a.







6 Conclusions
We have presented a new stabilised nodal spectral element model for simulation of fully nonlinear water wave propagation based on a Mixed Eulerian Lagrangian (MEL) formulation. The stability issues associated with mesh asymmetry as reported in RobertsonSherwin1999 () is resolved by using a hybrid mesh with a quadrilateral layer with interfaces aligned with the vertical direction to resolve the free surface and layer just below this level. Our linear stability analysis confirms that the temporal instability associated with triangulated meshes can be fixed, paving the way for considering remaining nonlinear stability issues. By combining a hybrid mesh strategy with the ideas described in EngsigKarupEtAl2016 () on stabilisation of free surface formulation with quartic nonlinear terms, the model is stabilised by using exact quadratures to effectively reduce aliasing errors and mild spectral filtering to add some artificial viscosity to secure robustness for marginally resolved flows. This is combined with a remeshing strategy to counter element deformations that may lead to numerical illconditioning and in the worst cases breakdown if not used. Our numerical analysis confirm this strategy to work well for the steepest nonlinear water waves when using this stabilised nonlinear MEL formulation. In the MEL formulation we avoid the use of a transform of the vertical coordinate, to introduce bodies of arbitrary geometry that is resolved using highorder curvilinear elements of an unstructured hybrid mesh. By using these elements, only few elements are needed to resolve the kinematics in the water column both with and without complex body surfaces. With this new methodology, we validate the model by revisiting known strenuous benchmarks for fully nonlinear wave models, e.g. solitary wave propagation and reflection, wavebody interaction with a submerged fixed cylinder and wavebody interaction with a fixed surfacepiercing structure in the form of a pontoon. The numerical results obtained are excellent compared with other published results and demonstrate the high accuracy that can be achieved with the highorder spectral element method.
In ongoing work, we are extending the new stabilised spectral element solvers towards advanced and realistic nonlinear hydrodynamics applications in three space dimensions (cf. the model based on an Eulerian formulation in three space dimensionsEngsigKarupEtAl2016ISOPE ()) by extending the current computational framework to also handle moving and floating objects of arbitrary body shape. In ongoing work, we consider freely moving bodies and surfacepiercing structures with nonvertical boundary at the bodysurface intersections and hybrid modelling approaches PaulsenEtAl2014 (); VerbruggheEtAl2016 ().
References
References
 [1] R. Beck and A. Reed. Modern computational methods for ships in a seaway. Trans. Soc. of Naval Architects and Marine Eng., 109:1–52, 2001.
 [2] T. Bjøntegaard and E. M. Rønquist. Accurate interfacetracking for Arbitrary LagrangianEulerian schemes. Journal of Computational Physics, 228(12):4379 – 4399, 2009.
 [3] X. Cai, H. P. Langtangen, B. F. Nielsen, and A. Tveito. A finite element method for fully nonlinear water waves. J. Comput. Phys., 143:544–568, July 1998.
 [4] L. Cavaleri, J.H.G.M. Alves, F. Ardhuin, A. Babanin, M. Banner, K. Belibassakis, M. Benoit, M. Donelan, J. Groeneweg, T.H.C. Herbers, P. Hwang, P.A.E.M. Janssen, T. Janssen, I.V. Lavrenov, R. Magne, J. Monbaliu, M. Onorato, V. Polnikov, D. Resio, W.E. Rogers, A. Sheremet, J. McKee Smith, H.L. Tolman, G. van Vledder, J. Wolf, and I. Young. Wave modelling – the state of the art. Progress in Oceanography, 75(4):603 – 674, 2007.
 [5] M. S. Celebi, M. H. Kim, and R. F. Beck. Fully Nonlinear 3D Numerical Wave Tank Simulation. J. Ship Res., 42(1):35–45, 1998.
 [6] M. J. Chern, A. G. L. Borthwick, and R. Eatock Taylor. A Pseudospectral transformation model of 2D Nonlinear Waves. Journal of Fluids and Structures, 13(5):607 – 630, 1999.
 [7] C. Chian and R. C. Ertekin. Diffraction of solitary waves by submerged horizontal cylinders. Wave Motion, 15(2):121–142, 1992.
 [8] T. Christiansen, A. P. EngsigKarup, and H. B. Bingham. Efficient pseudospectral model for nonlinear water waves. In Proceedings of The 27th International Workshop on Water Waves and Floating Bodies (IWWWFB), 2012.
 [9] G. Clauss and U. Steinhagen. Numerical simulation of nonlinear transient waves and its validation by laboratory data. International Journal of the International Offshore and Polar Engineering Conference (ISOPE), Brest, France, III:368–375, 1999.
 [10] A. Clement and S. Mas. Hydrodynamic forces induced by a solitary wave on a submerged circular cylinder. Proceedings of the International Offshore and Polar Engineering Conference (ISOPE), 1:339–347, 1995.
 [11] M. J. Cooker, D. H. Peregrine, C. Vidal, and J. W. Dold. The interaction between a solitary wave and a submerged cylinder. Journal of Fluid Mechanics, 215(1):1–22, 1990.
 [12] M. J. Cooker, P. D. Weidman, and D. S. Bale. Reflection of a highamplitude solitary wave at a vertical wall. Oceanographic Literature Review, 1998.
 [13] D. G. Dommermuth and D. K. P. Yue. A highorder spectral method for the study of nonlinear gravity waves. J. Fluid Mech., 184:267–288, 1987.
 [14] G. Ducrozet, H. B. Bingham, A. P. EngsigKarup, and P. Ferrant. Highorder finite difference solution for 3D nonlinear wavestructure interaction. J. Hydrodynamics, Ser. B, 22(5):225–230, 2010.
 [15] D. Dutykh and D. Clamond. Efficient computation of steady solitary gravity waves. Wave Motion, 51(1):86–99, 2014.
 [16] A. P. EngsigKarup. Analysis of efficient preconditioned defect correction methods for nonlinear water waves. Int. J. Num. Meth. Fluids, 74(10):749–773, 2014.
 [17] A. P. EngsigKarup, H. B. Bingham, and O. Lindberg. An efficient flexibleorder model for 3D nonlinear water waves. J. Comput. Phys., 228:2100–2118, 2009.
 [18] A. P. EngsigKarup, C. Eskilsson, and D. Bigoni. A stabilised nodal spectral element method for fully nonlinear water waves. Journal of Computational Physics, 318(6):1–21, 2016.
 [19] A. P. EngsigKarup, C. Eskilsson, and D. Bigoni. Unstructured spectral element model for dispersive and nonlinear wave propagation. In Proc. 26th International Ocean and Polar Engineering Conference (ISOPE), Greece, 2016.
 [20] A. P. EngsigKarup, S. L. Glimberg, A. S. Nielsen, and O. Lindberg. Fast hydrodynamics on heterogenous manycore hardware. In Raphäel Couturier, editor, Designing Scientific Applications on GPUs, Lecture notes in computational science and engineering. CRC Press / Taylor & Francis Group, 2013.
 [21] A. P. EngsigKarup, J. S. Hesthaven, H. B. Bingham, and P. A. Madsen. Nodal DGFEM solutions of highorder Boussinesqtype equations. J. Engineering Math., 56:351–370, 2006.
 [22] A. P. EngsigKarup, M. G. Madsen, and S. L. Glimberg. A massively parallel GPUaccelerated model for analysis of fully nonlinear free surface waves. Int. J. Num. Meth. Fluids, 70(1), 2011.
 [23] P. Ferrant. Simulation of strongly nonlinear wave generation and wavebody interactions using a 3D MEL model. In Proceedings of the 21st ONR Symposium on Naval Hydrodynamics, Trondheim, Norway, pages 93–109, 1996.
 [24] P. Ferrant, D. Le Touzé, and K. Pelletier. Nonlinear timedomain models for irregular wave diffraction about offshore structures. International Journal for Numerical Methods in Fluids, 43(1011):1257–1277, 2003.
 [25] W. J. Gordon and C. A. Hall. Construction of curvilinear coordinate systems and their applications to mesh generation. Int. J. Num. Meth. Engng., 7:461–477, 1973.
 [26] M. Gouin, G. Ducrozet, and P. Ferrant. Propagation of 3d nonlinear waves over an elliptical mound with a highorder spectral method. European Journal of Mechanics  B/Fluids, 63:9 – 24, 2017.
 [27] D. M. Greaves, G. X. Wu, A. G. L. Borthwick, and R. Eatock Taylor. A moving boundary finite element method for fully nonlinear wave simulations. J. Ship Res., 41(3):181–194, 1997.
 [28] J. C. Harris, E. Dombre, M. Benoit, and S. T. Grilli. A comparison of methods in fully nonlinear boundary element numerical wave tank development. In 14émes Journées de l’Hydrodynamique, 2014.
 [29] J. C. Harris, E. Dombre, M. Benoit, and S. T. Grilli. Fast integral equation methods for fully nonlinear water wave modeling. In Proceedings of the Twentyfourth (2014) International Ocean and Polar Engineering Conference, 2014.
 [30] J. S. Hesthaven and T. Warburton. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer, 2008.
 [31] C. H. Kim, A. H. Clément, and K. Tanizawa. Recent research and development of numerical wave tanks  a review. Int. J. Offshore and Polar Engng., 9(4):241–256, 1999.
 [32] C. Klettner and I. Eames. Momentum and energy of a solitary wave interacting with a submerged semicircular cylinder. J. Fluid Mech., 708:576–595, 2012.
 [33] T. Korsmeyer, T. Klemas, J. White, and J. Phillips. Fast hydrodynamic analysis of large offshore structures. In Proceedings of The Annual International Offshore and Polar Engineering Conference (ISOPE), 1999.
 [34] H.O. Kreiss and J. Oliger. Comparison of accurate methods for the integration of hyperbolic equations. Tellus, 24:199–215, 1972.
 [35] B. Li and C. A. Fleming. A three dimensional multigrid model for fully nonlinear water waves. Coastal Engineering, 30:235–258, 1997.
 [36] P. Lin. A multiplelayer coordinate model for simulation of wave–structure interaction. Computers & Fluids, 35(2):147 – 167, 2006.
 [37] M. S. LonguetHiggins and E. D. Cokelet. The deformation of steep surface waves on water. I. A numerical method of computation. Proc. Roy. Soc. London Ser. A, 350(1660):1–26, 1976.
 [38] Q. W. Ma, editor. Advances in numerical simulation of nonlinear water waves, volume 11. World Scientific Publishing Co. Pte. Ltd., 2010.
 [39] Q. W. Ma, G. X. Wu, and R. Eatock Taylor. Finite element simulation of fully nonlinear interaction between vertical cylinders and steep waves. part 1: methodology and numerical procedure. Int. J. Num. Meth. Fluids, 36(3):265–285, 2001.
 [40] Q. W. Ma, G. X. Wu, and R. Eatock Taylor. Finite element simulations of fully nonlinear interaction between vertical cylinders and steep waves. part 2: numerical results and validation. Int. J. Num. Meth. Fluids, 36(3):287–308, 2001.
 [41] Q. W. Ma and S. Yan. Quasi ALE finite element method for nonlinear water waves. J. Comput. Phys., 212(1):52 – 72, 2006.
 [42] Q. W. Ma and S. Yan. QALEFEM for numerical modelling of nonlinear interaction between 3d moored floating bodies and steep waves. Int. J. Num. Meth. Engng., 78(6):713–756, 2009.
 [43] P. A. Madsen, H. B. Bingham, and H. Liu. A new boussinesq method for fully nonlinear waves from shallow to deep water. J. Fluid Mech., 462:1–30, 2002.
 [44] T. Maxworthy. Experiments on the collision between solitary waves. Journal of Fluid Mechanics, 76(1):177–185, 1976.
 [45] S. B. Nimmala, S. C. Yim, and S. T. Grilli. An efficient parallelized 3D FNPF numerical wave tank for largescale wave basin experiment simulation. Journal of Offshore Mechanics and Arctic Engineering, 135(2), 2013.
 [46] A. T. Patera. A Spectral element method for fluid dynamics: Laminar flow in a channel expansion. J. Comput. Phys., 54:468–488, 1984.
 [47] B. T. Paulsen, H. Bredmose, and H. B. Bingham. An efficient domain decomposition strategy for wave loads on surface piercing circular cylinders. Coastal Engineering, 86:57–76, 2014.
 [48] I. Robertson and S. J. Sherwin. Freesurface flow simulation using hp/spectral elements. J. Comput. Phys., 155:26–53, 1999.
 [49] Y.L. Shao and O. M. Faltinsen. A harmonic polynomial cell (HPC) method for 3D Laplace equation with application in marine hydrodynamics. J. Comput. Phys., 274(0):312 – 332, 2014.
 [50] P. Sibley, L. E. Coates, and K. Arumugam. Solitary wave forces on horizontal cylinders. Applied Ocean Research, 4(2):113–117, 1982.
 [51] K. Tanizawa. A nonlinear simulation method of 3D body motions in waves (1st report). J. Soc. Nav. Arch. Japan, 178:179–191, 1995.
 [52] W.T. Tsai and D. K. P. Yue. Computation of nonlinear freesurface flows. In Annual review of fluid mechanics, Vol. 28, pages 249–278. Annual Reviews, Palo Alto, CA, 1996.
 [53] M. S. Turnbull, A. G. L. Borthwick, and R. Eatock Taylor. Numerical wave tank based on a transformed finite element inviscid flow solver. Int. J. Num. Meth. Fluids, 42(6):641–663, 2003.
 [54] T. Verbrugghe, P. Troch, A. Kortenhaus, V. Stratigaki, and A.P. EngsigKarup. Development of a numerical modelling tool for combined near field and far field wave transformations using a coupling of potential flow solvers. In Proc. of 2nd International Conference on Renewable energies Offshore (RENEW 2016), At Lisbon, Portugal, 2016.
 [55] C.Z. Wang and G.X. Wu. A brief summary of finite element method applications to nonlinear wavestructure interactions. J. Marine. Sci. Appl., 10:127–138, 2011.
 [56] L. Wang, H. Tang, and Y. Wu. Study on interaction between a solitary wave and a submerged semicircular cylinder using acceleration potential. In Proceedings of The Annual International Offshore and Polar Engineering Conference (ISOPE), 2013.
 [57] P. Wang, Y. Yao, and M. P. Tulin. An efficient numerical tank for nonlinear water waves, based on the multisubdomain approach with bem. International Journal for Numerical Methods in Fluids, 20(12):1315–1336, 1995.
 [58] T. Warburton. An explicit construction for interpolation nodes on the simplex. J. Engineering Math., 56(3):247–262, 2006.
 [59] B. J. West, K. A. Brueckner, R. S. Janda, D. M. Milder, and R. L. Milton. A new numerical method for surface hydrodynamics. Journal of Geophysical Research: Oceans, 92(C11):11803–11824, 1987.
 [60] J.H. Westhuis. The numerical simulation of nonlinear waves in a hydrodynamic model test basin. PhD thesis, Department of Mathematics, University of Twente, The Netherlands, 2001.
 [61] J.H. Westhuis and A. J. Andonowati. Applying the finite element method in numerically solving the two dimensional freesurface water wave equations. In Proceedings of The 13th International Workshop on Water Waves and Floating Bodies (IWWWFB), 1998.
 [62] G. X. Wu and R. Eatock Taylor. Finite element analysis of twodimensional nonlinear transient water waves. Applied Ocean Res., 16(6):363 – 372, 1994.
 [63] G. X. Wu and R. Eatock Taylor. Time stepping solutions of the twodimensional nonlinear wave radiation problem. Ocean Engineering, 22(8):785 – 798, 1995.
 [64] M. Xue, H. Xu, Y. Liu, and D. K. P. Yue. Computations of fully nonlinear threedimensional wave–wave and wave–body interactions. part 1. dynamics of steep threedimensional waves. Journal of Fluid Mechanics, 438:11–39, 07 2001.
 [65] H. Yan and Y. Liu. An efficient highorder boundary element method for nonlinear wave–wave and wavebody interactions. Journal of Computational Physics, 230(2):402 – 424, 2011.
 [66] E. W. Yeung. Numerical methods in freesurface flows. In Annual review of fluid mechanics, Vol. 14, pages 395–442. Annual Reviews, Palo Alto, Calif., 1982.
 [67] R.W. Yeung. A hybridintegralequation method for timeharmonic freesurface flow. In Proc. 1st Int Conf. Num. Ship Hydrodynamics, Gaithersburg, Maryland, pages 581–607, 1975.
 [68] V. E. Zakharov. Stability of periodic waves of finite amplitude on the surface of a deep fluid. J. Appl. Mech. Tech. Phys., 9:190–194, 1968.
 [69] Q. Zhu, Y. Liu, and D. K. P. Yue. Threedimensional instability of standing waves. Journal of Fluid Mechanics, 496:213–242, 12 2003.
 [70] O. C. Zienkewitch. The finite element method in engineering science, 2nd edition. McGrawHill, New York, 1971.
 [71] O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu. Chapter 6  free surface and buoyancy driven flows. In O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu, editors, The Finite Element Method for Fluid Dynamics (Seventh Edition), pages 195 – 224. ButterworthHeinemann, Oxford, seventh edition edition, 2014.