# A Stabilised Nodal Spectral Element Method for Fully Nonlinear Water Waves

## Abstract

We present an arbitrary-order spectral element method for general-purpose simulation of non-overturning water waves, described by fully nonlinear potential theory. The method can be viewed as a high-order extension of the classical finite element method proposed by Cai et al (1998) [5], although the numerical implementation differs greatly. Features of the proposed spectral element method include: nodal Lagrange basis functions, a general quadrature-free approach and gradient recovery using global projections. The quartic nonlinear terms present in the Zakharov form of the free surface conditions can cause severe aliasing problems and consequently numerical instability for marginally resolved or very steep waves. We show how the scheme can be stabilised through a combination of over-integration of the Galerkin projections and a mild spectral filtering on a per element basis. This effectively removes any aliasing driven instabilities while retaining the high-order accuracy of the numerical scheme. The additional computational cost of the over-integration is found insignificant compared to the cost of solving the Laplace problem. The model is applied to several benchmark cases in two dimensions. The results confirm the high order accuracy of the model (exponential convergence), and demonstrate the potential for accuracy and speedup. The results of numerical experiments are in excellent agreement with both analytical and experimental results for strongly nonlinear and irregular dispersive wave propagation. The benefit of using a high-order – possibly adapted – spatial discretisation for accurate water wave propagation over long times and distances is particularly attractive for marine hydrodynamics applications.

## 1Introduction

Robust and cost-efficient time-dependent simulation of the propagation and transformation of water waves in both shallow near-shore and deeper off-shore areas is a computationally challenging and longstanding scientific problem for ocean, coastal and naval engineering applications. For example, fully non-linear wave simulations have been subject to research for a long time, and have still not yet entered common coastal and ocean engineering practice. One remaining key challenge is to resolve accurately highly nonlinear and dispersive wave propagation in maritime areas while taking into account varying bathymetry, the geometry of complex structures and their nonlinear interaction with fixed and floating structures. Resolving this problem lead to improved opportunities for using simulations in realistic marine regions as well as enabling experiments in numerical wave tanks of increasing fidelity. Furthermore, it is attractive to design a flexible computational framework that entails use on modern commodity workstations as well as high-performance computing systems. These goals dictate stringent requirements on the design of engineering tools, and this suggests that a generic tool for wave propagation should be based on

a general modelling basis for broadly describing relevant

*physics*,a generalised computational framework for the

*numerics*, andsoftware design for efficient mapping to modern and emerging many-core

*architectures*.

Our main objective is to meet each of these requirements via a spectral element based computational framework. Such a framework has already been used for several different applications within marine hydrodynamics such as the fully nonlinear potential flow equations (present work), Boussinesq and shallow water equations [20]. In this work we focus exclusively on the first two requirements and we pave the way for the fulfilment of the last one, which is being addressed in ongoing work. The use of modern many-core hardware is attractive for acceleration of the high-order spectral element framework and for enabling practical computations of realistic engineering problems [49].

### 1.1Choice of modelling basis for description of physics

During the last decades computationally efficient depth-integrated Boussinesq-type models have been widely adopted as essential tools for water wave modelling in the near-shore region; see e.g. [48]. For shorter waves, such as the ones arising in offshore and naval engineering, Boussinesq-type models are not applicable due to the limited accuracy in terms of dispersive and nonlinear properties. For these cases we have to turn to Computational Fluid Dynamics (CFD) models based on the Navier-Stokes equations, or fully nonlinear potential flow (FNPF) models. The CFD models take viscous effects into account; effects that may be important for breaking waves, load computations, boundary layer effects, etc. Even though CFD is often prohibitively expensive in terms of computational resources when considering simulation of entire sea states [22], it is widely used to quantify breaking wave loads and wave run-up on structures. CFD models are typically too dissipative as a result of the low-order accuracy imposed by computational limitations for large-scale wave simulations. In contrast, already today FNPF models can be used for long-time and large-scale wave simulations [12]. FNPF solvers can be used for resolution of full sea states in large marine or coastal areas where nonlinear waves interact with fixed or floating structures. The cons of the FNPF models are that they cannot account for non-overturning breaking waves and viscous effects. For these reasons it can be attractive to combine FNPF models (far-field) with CFD (near-field) in hybrid modelling approaches for wave structure interaction, cf. [55]. This hybrid approach enables better simulation of strong nonlinear wave structure interactions in areas where the local wave climates can not be predicted accurately via a FNPF model.

### 1.2On the quest towards developing numerical strategies for real-world applications

A review of existing conventional discretisation methods and applications reveals that historically the main emphasis has been on Finite Difference Methods (FDM), boundary element methods (BEM) and finite element methods (FEM) [37]. These methods have been designed for the concept of Numerical Wave Tanks (NWT). The main computational bottleneck in all such numerical solvers is the solution of a large linear system. In FDM and FEM the discretisation procedure leads to sparse linear systems due to the local support of discrete operators, while in BEM it is only the domain boundary that needs representation. The discretisation procedure for BEM is based on a surface integral formulation together with Green’s identities. This leads to dense non-symmetric matrix operators that cannot be solved in a straightforward way with linear asymptotic scaling. There has been some recent progress in bringing the asymptotic cost (scaling rate) down for BEM [29] for both matrix-vector multiplications and storage requirements using both high-order basis functions and the Fast Multipole Method (FMM) [27]. While this strategy can asymptotically achieve linear complexity n ( number of computational nodes) in work effort for the spatial solver, it has a large constant in front of this asymptotic scaling term due to the need of solving a dense linear system of equations. This leaves BEM less efficient compared to volume-based discretisation methods such as FDM and FEM solvers as suggested in [72]. We note, that BEM is particularly attractive as a near-field solver for cases where waves interact with complex geometries [75] and may be combined with a far-field solver such as FEM [73]. The overall efficiency and scalability of BEM [28] can be compared to efficient and massively parallel free surface hydrodynamics solvers such as [19] which can achieve very high efficiency and scalability using multigrid-type methods [42] for arbitrary sized discrete problems, in particular when the (possibly curvilinear multi-block) meshes are logically structured, e.g. as in [23].

### 1.3State-of-the-art in finite element methods for fully nonlinear water waves

Reviews on the state-of-the art of numerical models for freely propagating water waves are given in [37]. Our scope in the present work is restricted to FNPF solvers and FEM.

The use of FEM for fully nonlinear water waves was pioneered by Wu & Eatock Taylor (1994) [71]. Since then, the majority of solvers for fully nonlinear potential flow equations have been limited to second (low) order FEM schemes [72] based on a Mixed Eulerian Lagrangian method [43] for updating the free surface variables. This approach requires expensive mesh update techniques and may suffer from stability problems for deformed meshes [56]. However, it is particularly well suited for dealing with the interaction between waves and fixed or freely floating bodies [64]. To overcome this expensive re-meshing problem a less expensive Quasi Arbitrary Lagrangian Eulerian finite element method (QALE-FEM) was developed [47] and a novel mesh update technique for ALE is described in [3].

Very little work has been done based on purely Eulerian formulations, e.g. the classical approach where a -transformed formulation is used. Previous studies based on using the -transformation [5] have been limited to two spatial dimensions and flat bathymetries for wave-structure interaction applications. These are essentially limited to the inclusion of bottom-mounted possibly surface-piercing structures, despite that the method may be used as an efficient base solver for very large domains and from shallow to deep waters. More details about free surface solvers based on FEM can be found in [76].

While there are very few studies on high-order finite element methods for water wave applications, it is well known that large efficiency gains for time-dependent wave problems can be achieved by the use of high-order accurate methods [40]. One attractive class of methods are the Spectral Element Methods (SEM) [21] which rely on the strong theoretical foundations of Spectral Methods [6]. The Spectral Element Method was first used by Patera (1984) [54] for fluid dynamics problems and has since then gained popularity. It combines the best properties of spectral methods and classical finite element methods, namely, high accuracy and flexibility in the spatial representation of domains. For smooth problems, the use of high-order discretisation is generally an efficient way to balance accuracy and cost, since a high-order discretisation allow for more coarse spatial representation compared to a low-order method. High-order methods may also reduce the robustness of explicit time-stepping methods due to global restriction on the stable time step sizes, although, in recent works [16] it is demonstrated that for certain wave models that require operator inversion in the time-stepping, high-order discretization combined with explicit time stepping methods need not reduce robustness - not even for unstructured mesh methods – and thereby provide a strong basis for efficient tools. Furthermore, with the geometric flexibility provided by adaptive meshes and the support for cost-efficient simulation through - and -refinement strategies, the SEM is suitable for large-scale computing due to accurate temporal integration over long times via high-order discretisation that lead to small dispersive and dissipative errors. However, a key challenge for unstructured solvers such as FEM and SEM is achieving work effort that is linearly proportional to the number of computational processing nodes (weak scalability) for unstructured meshes used to represent general domains. This requires the design of efficient and scalable preconditioned iterative methods [5] and of data structures that map efficiently to modern many-core architectures [25].

### 1.4Paper contributions

The main objective of this work is the proper design and validation of the SEM framework for – ultimately large-scale – dispersive and nonlinear water wave propagation. This effort is crucial to enable engineering analysis of water waves in marine regions. This paper proposes for the first time a stabilised high-order spectral element method for solving the fully nonlinear potential flow equations. A rigorous assessment of the stabilised numerical model is carried out using numerical experiments and known benchmarks in two space dimensions.

## 2Governing equations

For the description of inviscid and irrotational fluid flows we introduce the set of fully nonlinear and dispersive free surface equations described by a two-variable potential flow formulation. This formulation can be derived from the Navier-Stokes equations, cf. [16].

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. The mathematical problem is to find a scalar velocity potential function and to determine the evolution of the free surface elevation . The notations are illustrated in Figure 1.

The Eulerian description of the unsteady kinematic and dynamic free surface boundary conditions can be expressed in the Zakharov form [74]. Find such that

We have introduced the ’’ symbol to denote functionals defined only on the free surface plane. The vertical component of the velocity is determined by first solving a Laplace problem

where describes the still water depth.

We are interested in a set of governing equations that can be used as a basis for efficient simulations with support for representing structures accurately. A basis for efficient simulations is the classical -transformation of the vertical coordinate

where is the height of the water column above the bottom. This transforms the fluid domain to a time-independent computational domain at the expense of time-varying coefficients. The notations are illustrated in Figure 2.

The main drawback of using the -transformation is that it needs to be non-singular, and thus exclude the geometric modelling of breaking waves. For general wave-structure problems, this restriction can be removed by discretising and solving the Laplace problem directly. Following [5], we express the -transformed system in a form where variable depth is accounted for. Let be the time-independent computational domain . The Jacobian of the map is then

The artificial scalar velocity function contains all information about the flow kinematics in the entire fluid volume. The velocity field can be determined from using the relation . All metric coefficients in can be evaluated from the known two-dimensional free surface and bottom positions at given instants of time. It is possible to discretise after completing the differentiations, cf. [15], however, this increases the complexity of this formulation.

### 2.1Boundary conditions

For the solution of the Laplace problem at every time step, the following free surface boundary condition is specified

while at vertical boundaries, impermeable wall boundary conditions are assumed

where is an outward pointing unit normal vector at . Due to symmetry at wall boundaries the domain boundary conditions at free surface variables imposed are

### 2.2Wave generation and absorption zones

We employ a general-purpose embedded penalty forcing technique [16] similar to the technique described in [8]. This technique can be used for both wave generation of regular and irregular waves as well as absorption in a numerical wave tank setup. It can be derived from the relaxation method described in [41] to avoid the preprocessing step and turn it into an equivalent source term to the governing equations

where is the nonlinear function for the PDE, a free surface state variable (, in Eq. ), the relaxation functions can be defined to avoid minimal reflections [13], and the source function defines the analytical representation of the wave input signal to be generated. We choose equal to time step size .

## 3Numerical discretisation

The governing partial differential equations are discretised in a generic computational framework based on the method of lines, where first a semi-discrete system of ordinary differential equations is formed by spatial discretisation using a nodal Spectral Element Method. We present the 2D formulation before presenting results of 2D numerical experiments.

### 3.1Weak Galerkin formulation and discretisation

We form a partition of the domain to obtain a tessellation of consisting of non-overlapping shape-regular elements such that with denoting the ’th element. We introduce the finite element approximation space of continuous, piece-wise polynomial functions of degree at most , .

#### Unsteady free surface equations

The weak formulation of the free surface equations takes the following form. Find where such that

for all . We introduce the finite-dimensional approximations

where is the set of global finite element basis functions with cardinal property at mesh nodes with the Kronecker Symbol. Substitute these expressions into and choose . The discretisation in one spatial dimension becomes

where the following global matrices have been introduced

The gradients of the free surface state variables are recovered as described in Section ?. Temporal integration of is performed using an explicit fourth-order Runge-Kutta method.

**Remark**: The free surface equations contain strongly nonlinear terms, up to fourth order. The discretisation of these terms calls for proper treatment to avoid aliasing effects. This will be addressed in Sections Section 3.3 and Section 3.4.

#### Spatial discretisation of the Laplace problem

Consider the discretisation of the governing equations for the -transformed Laplace problem . We seek to construct a linear system of the form

where is the total degrees of freedom in the discretisation.

The starting point is the weak formulation of the symmetric formulation that can be expressed as: find such that

for all where the boundary integrals vanish at domain boundaries where impermeable walls are assumed. The discrete system operator is defined via domain decomposition as

The elemental integrals are approximated through the change of variable

where is the Jacobian of the affine mapping and is the computational reference element. The global assembly of this operator preserves symmetry, and the resulting linear system is then modified to impose the Dirichlet boundary conditions at the free surface. Typically, is a large and sparse operator with a narrow band structure of non-zero elements in two space dimensions obtained by a proper permutation of the rows and numbering of nodes to reduce fill-in in the factorisation of Gaussian Elimination procedures. By using a symmetric reverse Cuthill-Mackee permutation, the bandwidth of the sparse matrix is minimised and the system can be efficiently solved by a sparse direct Gaussian elimination procedure. In 2D this leads to optimal and scalable work effort and is used in this work. For sparse symmetric positive definite systems such as Eq. , the iterative preconditioned conjugate gradient solver is an attractive choice when system sizes become large [5] as would be expected for typical applications in three space dimensions, since convergence is guaranteed and memory footprint is minimal.

#### A generic technique for gradient recovery

The gradients of the globally piece-wise continuous basis functions will be discontinuous across element interfaces in the classical sense. To guarantee global continuity of derivatives a gradient recovery procedure can be used. Several gradient recovery techniques are reviewed in [34]. In this work, a global gradient recovery technique is used within our spectral element framework.

The global approximation of components of the horizontal first derivative as a function is expressed as

By a global Galerkin projection

we can generate the linear systems of equations

to recover the coefficients of the expansion . This procedure is similar to the one described by [30] and used in other models, cf. [71].

**Remark**: In the FEM model described in [69] the gradient recovery step is related to stability of the numerics. It was found that a global projection method may lead to instability when using FEM. As we shall see in the numerical experiments in Section 5 we do not reach the same empirical conclusion.

Using this gradient recovering technique and after the solution of , it is possible to estimate the vertical free surface velocity to obtain closure in the free surface problem . The weak formulation is

From this we can construct a linear system of the form

where the discrete operators takes the form

From the solution of we obtain the vertical free surface velocities

where denotes an index set consisting of global numbers for the free surface nodes.

#### Nodal prismatic Lagrange finite elements in two space dimensions

We consider first construction of elements in two space dimensions, needed for representing the solutions of the governing equations. The expansions can be based on quadrilateral elements in the time-constant computational domain. These elements are formed by a subdivision of the horizontal free surface plane with possible irregular sized non-overlapping elements. Each of these elements can then be extended in the vertical from the surface to the bottom to form the quadrilateral elements. A spectrally accurate multivariate hierarchical polynomial expansion can be constructed by a tensor product of one-dimensional basis functions. This results in a quadrilateral element . We introduce the element basis functions

where is the ’th order orthonormal Jacobi polynomial on the interval with orthogonality property

These polynomials can be evaluated efficiently using a simple recurrence relation [31]. The basis have arbitrary polynomial orders in the horizontal plane and in the vertical plane. This makes it possible to tune the orders of the approximations to balance accuracy and efficiency needs in simulations.

For interpolating polynomials, the Unisolvence Theorem guarantees a unique connection between the hierarchical polynomial (modal) expansion and the corresponding Lagrange polynomial (nodal) expansion. Thus, for all , , and , and for each element , we have

which defines a relationship between modal and nodal coefficients in the form [33]

where a non-singular generalised Vandermonde matrix has been introduced. This connection can be exploited together with the duality in polynomial representation to define the local Lagrange basis functions with property and their derivatives

and the matrices defined in terms of the first derivatives of the modal basis functions

For use with integration on the local elements, the local nodal mass matrix is introduced

where orthonormality of the basis functions is exploited to avoid the use of discrete quadrature rules in the constructions, leaving the implementations *quadrature-free*.

By evaluating expressions such as at the chosen inter element node distribution we obtain the generic arbitrary-order elemental operators

that can be used in connection with the chain rule to define derivatives in physical space

with discrete counterparts based on collocation expressed in the form

The choice of the nodal distribution on the simplexes in two horizontal dimensions can be based on explicit construction of the nodal distribution set using, e.g., a blend and warp procedure [65]. This is used to define the high-order nodal basis functions [67].

### 3.2Generalised local element matrices via quadrature-free matrix-based operations

All global operators can be assembled from generalised local elemental operators. Consider global integrals in the general form

Approximate each of the integrands using finite element basis functions of the form

The global integrals can be reduced to local integrals through domain decomposition

and by inserting the approximated integrands, we obtain local integrals of the form

Expand the approximate integrands using nodal expansions such that

In the special cases, where and are differential operators, e.g.

the local integrals in can be evaluated by exploiting that for nodal differentiation matrices , cf. [6]. An implementation of the local elemental operator following the derivation just outlined, enables generic implementation in a single framework that has support for -adaptivity to balance accuracy and cost in computations. Furthermore, the integration of the integrands approximated using polynomial basis functions are without quadrature errors when quadrature-free matrix-based operations are carried out using evaluations on sufficiently fine meshes as described next.

### 3.3Removal of quadrature errors via higher-order numerical integration

Numerical integration of nonlinear terms is handled via super-collocation [38] and is expressed in the general form

where we have introduced the interpolation operator that maps the representation of a function on coarse mesh (subscript ) to that of a fine mesh (subscript ). Exploiting the result in , the integration (projection) of the polynomial representation can be done without quadrature errors via mass matrix based operations for integration of the form

where and . The computational cost of this integration step has complexities and due to the involved matrix operation. The relation between size of the coarse and fine variables is for exact integration of the polynomial representation of the coarse space basis functions when quartic terms are present. Therefore, the increase in cost is about four times for these operations.

### 3.4Dealing with aliasing by spectral filtering

The solution of governing equations that contain strongly nonlinear terms may pose a challenge for maintaining both accuracy and numerical stability in time when simulations are marginally resolved. In such cases, numerical instability may be related to aliasing effects. Aliasing results from evaluating interpolated products of functions, which when represented with insufficient or marginal resolution introduces errors in the functional representation. Possible remedies [38] focus on increasing resolution or introducing a proper stabilisation strategy to artificially dissipate such errors. While the former is simple, it is generally not considered feasible for large-scale systems or for long-time integration.

In this work, we employ a spectral filtering strategy exploiting the duality in the local element representation, cf. . On the ’th element the filtered local solution can be expressed as (assuming an order expansion in one space dimension here)

The filtering is applied only for the time-dependent free surface variables and . An exponential filter [32] can be used with cut-off low-pass filter index , to only affect the highest modes , such that

For example, by choosing the parameters a mild damping is achieved which gently removes five percent of the energy from only the highest mode in the basis. The filtering is done on a per element basis using a matrix-vector product

which has a work complexity of with total work effort proportional to . This local filtering matrix is constructed and used on all elements and applied repeatedly as necessary. We use the model basis in one space dimension [35]

to avoid introducing interface jumps or affect the mean through only filtering the higher modes (). If the filter is applied gently, i.e. by removing very little ’energy’ from the highest modes only, spectral accuracy can be recovered. Excessive filtering may reduce the convergence rate of the method albeit improve stability.

## 4Numerical properties of the model

### 4.1Temporal stability

For general schemes, the connection between nonlinear instability and explicit time-integration is important to understand [11]. To solve the governing equations efficiently, an explicit time integration method is preferred. Explicit schemes come with conditional stability in the form of the global CFL condition that dictates an upper bound for the choice of stable time step sizes where the constant . It is well-known that a main challenge for many high-order Spectral Element Methods is an unattractive scaling of the modulus of the eigenvalues of the discrete operators of the form where is expansion order and is the highest order of differentiation operator in the evolution equations. Typically, the constant is dependent on the minimum mesh size for an element in the mesh. This property may pose a severe problem in the accurate local representation of geometric features with small elements.

As discussed and shown in [16] for a linear finite difference scheme of the same governing equations, the stable time step size for explicit schemes has an upper bound given by the CFL condition , where is here still water depth and the gravitational acceleration constant. This bound is only dependent on the scale of the physics (still water depth) and the resolution chosen in the vertical. This result is demonstrated and shown in Figure 4, following [15] via discretisation and numerical eigenvalue analysis of the linearised system. This property is similar to the property inherent in several other wave models as described in the works [56]. This property implies that the CFL condition is tamed [66] in the sense that the time step size is not dictated by the numerics but only the physics (depth). Along the same line of the experiments described in [16], we find that there are only small changes in the CFL properties for nonlinear problems and hence small elements in the mesh do not impose any severe time step size restriction.

### 4.2Linear accuracy and dispersion properties

The accuracy of the numerical model is compared to the theoretical solution to the system of equations that arise when the system is subject to the assumption of small-amplitude waves (assume where wave height relative to wave length ). The theoretical solution for linear progressive monochromatic waves in one space dimension is given by [61]

with linear dispersion relation from Stokes theory , where is the wave number, is the gravitational acceleration (assumed as 9.81 ) and is the angular velocity with being the wave period.

In Figure 8 (a) results of - and -refinement strategies are presented and highlight an important advantage of the SEM. When solutions are smooth, it is possible to use high-order basis functions to improve the cost-efficiency of the method by allowing fewer degrees of freedom to be used to achieve the same level of accuracy of lower order methods. Most previous works have focused exclusively on classical FEM methods based on piece-wise linear approximations, which has a convergence rate that matches the curve for linear () basis functions.

The decisive criterion for choosing between different numerical strategies is to understand what is the amount of work (cost) to achieve a given level of accuracy. This is illustrated via results obtained via a semi-optimal brute-force procedure in Figures Figure 8 (b) and (c) where the work effort has been minimised^{1}

Clustering mesh nodes more densely closer to the free surface can improve accuracy in linear dispersion without increasing the CPU time significantly [45]. This is confirmed by results presented in Figure 10. The results show that the use of a high-order method or clustered vertical distribution of low-order elements is a must for an accurate approximation of dispersion in applications where is large (short wave length relative to depth). This highlights that one can tune the accuracy by proper choice of discretisation parameters. An important implication of these results is that the vertical node distribution can be used to control the range of validity of the model in terms of dispersive properties, i.e. a numerical truncation counterpart to the analytic truncation used in Boussinesq-type models. Interestingly, flexible-order finite difference solvers with cosine-clustered vertical distributions of nodes appears to be superior in linear dispersion compared to single layer SEM for equal number of vertical nodes. However, multi-layer SEM with cosine clustered elements of vertical order almost match the FDM, cf. [16].

### 4.3On nonlinear accuracy, stability and kinematics properties

The accurate computation of kinematics is essential for load predictions in wave-structure applications, e.g. for offshore foundations of wind turbines. For nonlinear waves, exact stream function (SF) wave solutions of permanent form [9] based on assuming a flat sea bed can be used to assess the accuracy with respect to variation in dimensionless nonlinearity () and dispersion () parameters. This is done by solving the Laplace problem first. Then, from the scalar velocity potential solution we can calculate the vertical free surface velocity and compare with exact results. Numerical results for linear, weakly nonlinear and strongly nonlinear SF waves in combinations of shallow, intermediate and deep waters are presented in Figures Figure 13 and Figure 16. Here we use only one layer of elements in the vertical and a fixed number of elements in the horizontal direction. Convergence is achieved by the variation of the polynomial order to achieve fast -convergence. All tests show convergence with increasing resolution as expected. When depth or nonlinearity increases more resolution is required. Similar tests were carried out for a flexible-order finite difference model in [2]. An immediate conclusion is that the SEM has similar resolution requirements as the corresponding finite difference solver to match the order of accuracy for nonlinear applications. This highlights that from a purely algorithmic (implementation independent) viewpoint there seems to be no significant tradeoffs in introducing geometric flexibility through this choice of discretisation. In Figure 18, the accuracy of the kinematics computations is shown for intermediate depth and deep water for very steep nonlinear stream function waves. Excellent agreement is found between exact and computed results for both intermediate and deep waters, which is difficult to represent in conventional wave propagation models due to lack of resolution.

## 5Numerical experiments

We examine different test cases and benchmarks, that inspect different properties of the numerical model that serves as validation.

### 5.1Stabilised nonlinear wave propagation of stream function waves

To test the robustness and accuracy of the numerical method, we consider the strenuous test case of propagating nonlinear stream function waves near the theoretical limit of wave steepness and nonlinearity. Few numerical schemes can resolve such waves accurately. Results from representative numerical experiments are presented in Table 1. The table shows how a standard Galerkin scheme unsurprisingly fails to be stable for all temporal resolutions chosen. Improvement in stability is achieved by using higher order numerical integration (cf. Section 3.3) integrands in the Galerkin integrals. This is done by interpolating the integrands to a finer mesh on the reference element corresponding to a basis with polynomial basis expansion order of resulting in terms of less than quartic order nonlinear to be subject to *over-integration*. This implies that exact integration is used for all nonlinear terms to within the order of accuracy of the numerical scheme. This is compared to using no over-integration but instead trying to stabilise only using a *gentle* spectral filtering strategy that caps off the highest 5% in the highest modes of the modal expansion. Both of these strategies fail to be stable. Instead, if the two strategies are combined, where the over-integration effectively removes any aliasing in the evaluation of the strongly nonlinear terms, the mild spectral filter dissipates just enough energy for the model to stabilise completely. We find that over-integration is only needed in the free-surface equations, leading to a marginal increase in computational cost, which is anyway driven by the Laplace solver. This is highlighted in the table, where the time stepping cost is only increased by approximately 15% when using over-integration and spectral filtering in comparison with a standard Galerkin formulation based on expansion order . This additional cost of over-integration would be even less significant for larger simulations.

Accuracy |
Cost/ | ||||
---|---|---|---|---|---|

40 | 80 | 160 | |||

No filter + | NaN | NaN | NaN | 1.00 | |

No over-integration | NaN | NaN | NaN | ||

No filter + | 1.3608e-03 | 7.5021e-04 | 7.6731e-04 | 1.02 | |

Over-integration | NaN | NaN | NaN | ||

Filter + | NaN | 5.6019e-04 | 8.0508e-03 | 1.12 | |

No over-integration | NaN | NaN | 1.2705e-03 | ||

NaN | NaN | 7.8374e-03 | |||

NaN | NaN | NaN | |||

Filter + | 1.3943e-03 | 7.0651e-04 | 1.0102e-03 | 1.15 | |

Over-integration | 7.4032e-03 | 4.3313e-03 | 7.0332e-03 | ||

7.2826e-02 | 5.7642e-02 | 7.5093e-02 | |||

### 5.2Convergence tests and high-order accuracy

To confirm the high-order accuracy of the model and evaluate the influence of spectral filtering we have carried out convergence tests using the exact nonlinear stream function wave solutions for parameters (dispersion) and , and of maximum wave steepness (nonlinearity). The results for tests of the proposed Galerkin scheme with over-integration without spectral filtering and with spectral filtering using a cap of 1% of the highest modal mode are presented in Figure 21. The results confirm the high-order convergence for the spatial spectral element discretisation. Particularly, this is clear for the mildly nonlinear wave. With increasing nonlinearity more spatial resolution is required to accurately resolve all harmonic modes of the solution. The gentle filtering is found to reduce accuracy and has some detrimental effect on the convergence. For increasing order of the local basis functions, these effects become less significant.

### 5.3Harmonic generation over a submerged bar

We consider the classical benchmark for wave transformation due to a submerged bar to test the accuracy of the SEM model. The test is often used for validation of deterministic dispersive and nonlinear wave models since it can be compared to experimental laboratory data and other known results in the literature, e.g., see [24]. The numerical wave tank is illustrated in Figure 22.

The experiment was originally proposed in [1] and subsequently an equivalent scaled experiment was carried out and described in [44]. We consider the setup for Case A in the original experiment. The weakly nonlinear input wave is generated in the numerical model using a regular stream function solution at undisturbed depth m with a wave height of cm and a wave period of s. The input wave is generated and propagates towards the submerged bar on a constant bottom. During the propagation over the bar, the wave will undergo a transformation resulting in a steepening and shortening of wavelengths due to nonlinear shoaling effects. At the top of the bar, the bound harmonics will be released as free harmonics (harmonic generation) decomposing the wave into shorter waves that propagate freely. Thus, to attain high accuracy in the calculations we need to use a model that can handle nonlinear wave-wave interactions and have accurate dispersion properties to capture the correct wave speed of the free harmonics after the bar. Taking advantage of the unstructured SEM the elements are of similar sizes but adjusted to have the interfaces positioned where the depth function has kinks in the first order gradients. Thus the -transformation is also local to the elements and entirely valid throughout the domain. The results presented in Figure 30 are based on a deterministic simulation where we have used 103 elements in the horizontal and one vertical layer with a multivariate basis of polynomial order in both the horizontal and vertical directions. A CFL condition with Courant number is used for defining the time step size. All results are found to be in excellent agreement with the experimental data with some qualitatively minor differences in phase between experimental and computed results which compare well to other published results.

### 5.4Irregular waves shoaling on a slope

Mase and Kirby [51] conducted laboratory experiments of shoaling and breaking of irregular waves on a constant slope. The water depth offshore the slope is 0.47 m and the slope is 1/20. The incoming irregular waves were generated by a Pierson-Moskowitz (PM) spectrum with peak frequency of 1.0 Hz. This set of experiments has been used for testing Boussinesq models in e.g. [68].

As the SEM model presently is limited to non-breaking waves we follow the approach of [68], i.e. the constant slope part of the wave tank is truncated and exchanged with a flat bottom at a depth such that no breaking occurs. For the SEM model we limit the still water depth to be no less than 0.19 m. The numerical wave tank is illustrated in Figure 31. The wave gauges are located at and m, respectively. Wave generation and absorption is performed with 4 m long relaxation zones. Further, as the random phase angle of the spectrum is unknown, the incoming wave train is generated in the numerical model based on a FFT of the measured waves at wave gauge 1 located at the toe of the slope. The wave kinematics is given by superposition of the linear wave modes.

The computational domain is partitioned into elements of length 0.1 m with an expansion basis of polynomial order in the horizontal and vertical directions. Figure 35 shows the simulated and experimental free surface elevation at the different wave gauge locations. There is a reasonably good fit to the experimental data, with some minor differences in peak and through amplitudes. There are two main reasons for the discrepancy. First, the shorter waves of the PM spectra would require a very high polynomial order in order to be properly resolved and are thus diffused. Secondly, the assumption of linear superposition in the generation zone is not correct, giving that the incoming wave train does not exactly match the recorded.

## 6Conclusions

We have presented a spectral element model for simulation of fully nonlinear water wave propagation. The main advantages of using the spectral element method is the opportunity for balancing high accuracy with unstructured meshes which can be adapted to geometry of arbitrary shape (sharp corners, curvilinear features, etc.) or features of the solution such as large gradients. The spectral element method’s dual roads towards convergence, namely - and -adaptivity, allows for balancing accuracy and cost effectively.

A two-dimensional spectral element model for fully nonlinear potential flow is implemented in a general computational framework using quadrature-free construction of local element operators. We have used arbitrary-order multivariate Lagrange (nodal) basis functions in space and an explicit fourth order Runge-Kutta method in time. The explicit time stepping is efficient as the model has a bounded eigenspectrum, and the stable time step sizes are governed by still water depth and the vertical resolution. The model is stabilised by using over-integration to effectively reduce aliasing errors and mild spectral model filtering to add some artificial viscosity to secure robustness for marginally resolved flows.

The proposed model was shown to have a convergence rate of order , although it is difficult to keep the sharp convergence rate for the most nonlinear waves due to high resolution requirements. Numerical experiments demonstrate that the stabilised model is both robust and accurate. It was shown how the spectral accuracy can be used to substantially reduce the number of degrees of freedom per wavelength and compare well to properties that can also be achieved with finite difference time domain schemes, but with the additional benefit of geometric flexibility. Also, we illustrate how the order of the basis function in vertical dimension can be used to control the range of validity of the model in terms of dispersive properties, i.e. a numerical truncation counterpart to the analytic truncation used in standard Boussinesq-type models. While the methodology is efficient in two space dimensions, particular attention is to be given to further improve numerical efficiency via efficient preconditioning methods that maintain high efficiency for general unstructured grids.

In ongoing work, we aim at considering advanced nonlinear hydrodynamics problems by extending the current framework to also handle moving and floating objects. The present model needs to be further improved to handle run-up for calculations in the swash zone and a breaking wave model needs to be included for realistic applications. The proposed methodology in two space dimensions can be extended without conceptual modifications for three space dimensions, paving the road for marine hydrodynamics applications in 3D for offshore structures. Main challenges for enabling efficient computations in three space dimensions would be to (i) switch from a direct solver to an iterative solver strategy and identify an efficient and scalable preconditioning strategy, and (ii) via proper software design enable efficient mapping to modern and emerging many-core architectures for accelerated performance of the computational framework.

## References

### Footnotes

- Tests done using a sequential proof-of-concept code on a laptop with a 2,3 GHz Intel Core I7 processor and 16 GB 1600 MHz DDR3 RAM.

### References

**Numerical simulation of nonlinear-wave propagation over a bar.**

S. Beji and J. A. Battjes.*Coastal Engineering*, 23:1–16, 1994.**On the accuracy of finite-difference solutions for nonlinear water waves.**

H. B. Bingham and H. Zhang.*J. Engineering Math.*, 58:211–228, 2007.**Mesh update techniques for free-surface flow solvers using spectral element method.**

E. Bouffanais and M. O. Deville.*J. Sci. Comput.*, 27(1-3):137–149, June 2006.**A reasoned overview on boussinesq-type models: the interplay between physics, mathematics and numerics.**

M. Brocchini.*Proc. Roy. Soc. London Ser. A*, 469(2160):1–27, 2013.**A finite element method for fully nonlinear water waves.**

X. Cai, H. P. Langtangen, B. F. Nielsen, and A. Tveito.*J. Comput. Phys.*, 143:544–568, July 1998.*Spectral methods - Fundamentals in single domains*.

C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Springer, 2006.**Numerical simulation of nonlinear transient waves and its validation by laboratory data.**

G. Clauss and U. Steinhagen.*International Journal of the International Offshore and Polar Engineering Conference (ISOPE), Brest, France*, III:368–375, 1999.**Nonlinear simulation of transient free surface flows.**

R. Cointe. In*In proceedings of the 5th International Conference in Numerical Ship Hydrodynamics*, 1989.**Stream function representation of nonlinear ocean waves.**

R. G. Dean.*J. Geophys. Res.*, 70:4561–4572, 1965.*High Order Methods for Incompressible Fluid Flow*.

M. O. Deville, P. F. Fischer, and E. H. Mund. Cambridge University Press, 2002.**The numerical computation of freely propagating time-dependent irrotational water waves.**

F. Dias and T. J. Bridges.*Fluid Dynam. Res.*, 38(12):803–830, 2006.**Rogue waves in large-scale fully-non-linear high-order-spectral simulations.**

G. Ducrozet and P. Ferrant. In*Proc. 22nd International Workshop on Water Waves and Floating Bodies (IWWWFB), Croatia*, 2007.*Unstructured Nodal DG-FEM solution of high-order Boussinesq-type equations*.

A. P. Engsig-Karup. PhD thesis, Department of Mechanical Engineering, Technical University of Denmark, 2006.**Analysis of efficient preconditioned defect correction methods for nonlinear water waves.**

A. P. Engsig-Karup.*Int. J. Num. Meth. Fluids*, 74(10):749–773, 2014.**An efficient flexible-order model for 3D nonlinear water waves.**

A. P. Engsig-Karup, H. B. Bingham, and O. Lindberg.*J. Comput. Phys.*, 228:2100–2118, 2009.**Fast hydrodynamics on heterogenous many-core hardware.**

A. P. Engsig-Karup, L. S. Glimberg, A. S. Nielsen, and O. Lindberg. In Raphäel Couturier, editor,*Designing Scientific Applications on GPUs*, Lecture notes in computational science and engineering. CRC Press / Taylor & Francis Group, 2013.**Nodal DG-FEM solutions of high-order Boussinesq-type equations.**

A. P. Engsig-Karup, J. S. Hesthaven, H. B. Bingham, and P. Madsen.*J. Engineering Math.*, 56:351–370, 2006.**DG-FEM solution for nonlinear wave-structure interaction using boussinesq-type equations.**

A. P. Engsig-Karup, J. S. Hesthaven, H. B. Bingham, and T. Warburton.*Coastal Engineering*, 55:197–208, 2008.**A massively parallel GPU-accelerated model for analysis of fully nonlinear free surface waves.**

A. P. Engsig-Karup, M. G. Madsen, and S. L. Glimberg.*Int. J. Num. Meth. Fluids*, 70(1), 2011.**On devising boussinesq-type models with bounded eigenspectra: One horizontal dimension.**

C. Eskilsson and A. P. Engsig-Karup.*J. Comput. Phys.*, 271(0):261 – 280, 2014.**The next step in coastal numerical models: spectral/hp element methods?**

C. Eskilsson, A. P. Engsig-Karup, S. J. Sherwin, J. S. Hesthaven, and L. Bergdahl. In*Proceedings of the WAVES2005 Conference, Madrid*, 2005.**CFD study of the overtopping discharge of the wave dragon wave energy converter.**

C. Eskilsson, J. Palm, J. P. Kofoed, and E. Friis-Madsen. In*Proc. 1st International Conference of Renewable Energies Offshore*. ASCE, 2015.*Designing Scientific Software for Heterogenous Computing - With application to large-scale water wave simulations*.

S. L. Glimberg. PhD thesis, Department of Applied Mathematics and Computer Science, Technical University of Denmark, Kongens Lyngby, Denmark., 2013.**Wave evolution over submerged sills: tests of a high-order boussinesq model.**

M. F. Gobbi and J. T. Kirby.*Coastal Engineering*, 37:57–96, 1999.**Using GPUs to improve multigrid solver performance on a cluster.**

D. Göddeke, R. Strzodka, J. Mohd-Yusof, P. McCormick, H. Wobker, C. Becker, and S. Turek.*Int. J. Comput. Sci. Eng.*, 4(1):36–55, 2008.**A moving boundary finite element method for fully nonlinear wave simulations.**

D. M. Greaves, G. X. Wu, A. G. L. Borthwick, and R. Eatock Taylor.*J. Ship Res.*, 41(3):181–194, 1997.**A fast algorithm for particle simulations.**

L. Greengard and V. Rokhlin.*Journal of Computational Physics*, 135(2):280 – 292, 1997.**A comparison of methods in fully nonlinear boundary element numerical wave tank development.**

J. C. Harris, E. Dombre, M. Benoit, and S. T. Grilli. In*14émes Journées de l’Hydrodynamique*, 2014.**Fast integral equation methods for fully nonlinear water wave modeling.**

J. C. Harris, E. Dombre, M. Benoit, and S. T. Grilli. In*Proceedings of the Twenty-fourth (2014) International Ocean and Polar Engineering Conference*, 2014.**A comparison of gradient recovery methods in finite-element calculations.**

D. M. Hawken, P. Townsend, and M. F. Webster.*Communications in Applied Numerical Methods*, 7(3):195–204, 1991.*Spectral Methods for Time-Dependent Problems*.

J. S. Hesthaven, S. Gottlieb, and D. Gottlieb. Cambridge Monographs on Applied And Computational Mathematics 21. Cambridge University Press, Cambridge, UK, 2007.**Filtering in legendre spectral methods.**

J. S. Hesthaven and R. M. Kirby.*Math. Comp.*, 77(263):1425–1452, 2003.*Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications*.

J. S. Hesthaven and T. Warburton. Springer, 2008.**Local and global smoothing of discontinuous finite element functions using a least squares method.**

E. Hinton and J. S. Campbell.*Int. J. Num. Meth. Engng.*, 8(3):461–480, 1974.*Spectral/*.*hp*element methods for CFD

G. E. Karniadakis and S. J. Sherwin. Oxford University Press, 1999.*Spectral/hp element methods for computational fluid dynamics*.

G. E. Karniadakis and S. J. Sherwin. Oxford University Press, 2 edition, 2005.**Recent research and development of numerical wave tanks - a review.**

C. H. Kim, A. H. Clément, and K. Tanizawa.*Int. J. Offshore and Polar Engng.*, 9(4):241–256, 1999.**De-alising on non-uniform grids: algorithms and applications.**

R. M. Kirby and G. E. Karniadakis.*J. Comput. Phys.*, 191:249–264, 2003.*Implementing Spectral Methods for Partial Differential Equations - Algorithms for Scientists and Engineers*.

D. Kopriva. Springer, 2009.**Comparison of accurate methods for the integration of hyperbolic equations.**

H.-O. Kreiss and J. Oliger.*Tellus*, 24:199–215, 1972.**Open boundaries in short wave simulations - a new approach.**

J. Larsen and H. Dancy.*Coastal Engineering*, 7:285–297, 1983.**A three dimensional multigrid model for fully nonlinear water waves.**

B. Li and C. A. Fleming.*Coastal Engineering*, 30:235–258, 1997.**The deformation of steep surface waves on water. I. A numerical method of computation.**

M. S. Longuet-Higgins and E. D. Cokelet.*Proc. Roy. Soc. London Ser. A*, 350(1660):1–26, 1976.**Projects 13G: Kinematics of waves breaking partially on an offshore bar: LDV measurements for waves with and without a net onshore current.**

H. R. Luth, B. Klopman, and N. Kitou., 1994.*Technical report H1573*, Delft Hydraulics**Finite element simulation of fully non-linear interaction between vertical cylinders and steep waves. part 1: methodology and numerical procedure.**

Q. W. Ma, G. X. Wu, and R. Eatock Taylor.*Int. J. Num. Meth. Fluids*, 36(3):265–285, 2001.**Finite element simulations of fully non-linear interaction between vertical cylinders and steep waves. part 2: numerical results and validation.**

Q. W. Ma, G. X. Wu, and R. Eatock Taylor.*Int. J. Num. Meth. Fluids*, 36(3):287–308, 2001.**Quasi ALE finite element method for nonlinear water waves.**

Q. W. Ma and S. Yan.*J. Comput. Phys.*, 212(1):52 – 72, 2006.**A review of boussinesq-type equations for gravity waves.**

P. A. Madsen and H. A. Schäffer.*In Advances in Coastal and Ocean Engineering*, 5:1–95, 1999.**Finite element assembly strategies on multi-core and many-core architectures.**

G. R. Markall, A. Slemmer, D. A. Ham, P. H. J. Kelly, C. D. Cantwell, and S. J. Sherwin.*International Journal for Numerical Methods in Fluids*, 71(1):80–97, 2013.**Openacc acceleration of the nek5000 spectral element code.**

S. Markidis, J. Gong, M. Schliephake, E. Laure, A. Hart, D. Henty, K. Heisey, and P. Fischer.*International Journal of High Performance Computing Applications*, 29(3):311–319, 2015.**Hybrid frequency-domain KdV equation for random wave transformation.**

H. Mase and J. T. Kirby. In*Proc. 23rd International Conference on Coastal Engineering*, pages 474–487. ASCE, 1992.**Dealiasing techniques for high-order spectral element methods on regular and irregular grids.**

G. Mengaldo, D. De Grazia, D. Moxey, P. E. Vincent, and S. J. Sherwin.*J. Comput. Phys.*, 299:56 – 81, 2015.**An efficient 3d-fnpf numerical wave tank for virtual large-scale wave basin experiments.**

S. B. Nimmala, S. C. Yim, and S. T. Grilli. In*In Proc. 31st Intl. Conf. on Ocean, Offshore and Arctic Engineering*, 2012.**A Spectral element method for fluid dynamics: Laminar flow in a channel expansion.**

A. T. Patera.*J. Comput. Phys.*, 54:468–488, 1984.**An efficient domain decomposition strategy for wave loads on surface piercing circular cylinders.**

B. T. Paulsen, H. Bredmose, and H. B. Bingham.*Coastal Engineering*, 86:57–76, 2014.**Free-surface flow simulation using hp/spectral elements.**

I. Robertson and S. J. Sherwin.*J. Comput. Phys.*, 155:26–53, 1999.**A harmonic polynomial cell (HPC) method for 3D Laplace equation with application in marine hydrodynamics.**

Y.-L. Shao and O. M. Faltinsen.*J. Comput. Phys.*, 274(0):312 – 332, 2014.**A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation.**

F. Shi, J. T. Kirby, J. C. Harris, J. D. Geiman, and S. T. Grilli.*Ocean Modelling*, 43–44:36–51, 2012.**Assessment of an advanced finite element tool for the simulation of fully-nonlinear gravity water waves.**

J. Spinneken, V. Heller, S. Kramer, M. Pigott, and A. Viré. In*Proceedings of The Annual International Offshore and Polar Engineering Conference (ISOPE) 2012*, 2012.**Velocity calculation methods in finite element based mel formulation.**

V. Sriram, S. A. Sannasiraj, and V. Sundar. In Q. W. Ma, editor,*Advances in Numerical Simulation of Nonlinear Water Waves*, volume 11 of*Advances in Coastal and Ocean Engineering*, pages 203–244. World Scientific, 2010.*Hydrodynamics of coastal regions*.

I. A. Svendsen and I. G. Jonsson. Technical University of Denmark, 2001.**The state of the art on numerical wave tank.**

K. Tanizawa. In*Proceedings of 4th Osaka Colloquium on Seakeeping Performance of Ships*, pages 95–14, 2000.**Numerical wave tank based on a -transformed finite element inviscid flow solver.**

M. S. Turnbull, A. G. L. Borthwick, and R. Eatock Taylor.*Int. J. Num. Meth. Fluids*, 42(6):641–663, 2003.**A brief summary of finite element method applications to nonlinear wave-structure interactions.**

C.-Z. Wang and G.-X. Wu.*J. Marine. Sci. Appl.*, 10:127–138, 2011.**An explicit construction for interpolation nodes on the simplex.**

T. Warburton.*J. Engineering Math.*, 56(3):247–262, 2006.**Taming the CFL number for discontinuous Galerkin methods on structured meshes.**

T. Warburton and T. Hagstrom.*SIAM J. Numer. Anal.*, 46(6):3151–3180, September 2008.**Basis functions for triangular and quadrilateral high-order elements.**

T. C. Warburton, S. J. Sherwin, and G. E. Karniadakis.*SIAM J. Sci. Comput.*, 20(5):1671–1695, April 1999.**Time-dependent numerical code for extended boussinesq equations.**

G. Wei and J. T. Kirby.*J. Waterway, Port, Coastal and Ocean Eng., ASCE*, 121(5):251–261, 1995.*The numerical simulation of nonlinear waves in a hydrodynamic model test basin*.

J.-H. Westhuis. PhD thesis, Department of Mathematics, University of Twente, The Netherlands, 2001.**Applying the finite element method in numerically solving the two dimensional free-surface water wave equations.**

J.-H. Westhuis and A. J. Andonowati. In*Proceedings of The 13th International Workshop on Water Waves and Floating Bodies (IWWWFB)*, 1998.**Finite element analysis of two-dimensional non-linear transient water waves.**

G. X. Wu and R. Eatock Taylor.*Applied Ocean Res.*, 16(6):363 – 372, 1994.**Time stepping solutions of the two-dimensional nonlinear wave radiation problem.**

G. X. Wu and R. Eatock Taylor.*Ocean Engineering*, 22(8):785 – 798, 1995.**The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies.**

G. X. Wu and R. Eatock Taylor.*Ocean Engineering*, 30(3):387 – 400, 2003.**Stability of periodic waves of finite amplitude on the surface of a deep fluid.**

V. E. Zakharov.*J. Appl. Mech. Tech. Phys.*, 9:190–194, 1968.**Fully nonlinear wave interaction with freely floating non-wall-sided structures.**

B. Z. Zhou, G. X. Wu, and B. Teng.*Engineering Analysis with Boundary Elements*, 50:117 – 132, 2015.**Chapter 6 - free surface and buoyancy driven flows.**

O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu. In O. C. Zienkiewicz, R. L. Taylor, and P. Nithiarasu, editors,*The Finite Element Method for Fluid Dynamics (Seventh Edition)*, pages 195 – 224. Butterworth-Heinemann, Oxford, seventh edition edition, 2014.