Hamiltonian Finite Element Discretization for Nonlinear Free Surface Water Waves

A novel ﬁnite element discretization for nonlinear potential ﬂow water waves is presented. Starting from Luke’s Lagrangian formulation we prove that an appropriate ﬁnite element discretization preserves the Hamiltonian structure of the potential ﬂow water wave equations, even on general time-dependent, deforming and unstructured meshes. For the time-integration we use a modiﬁed Störmer–Verlet method, since the Hamiltonian system is non-autonomous due to boundary surfaces with a prescribed motion, such as a wave maker. This results in a stable and accurate numerical discretization, even for large amplitude nonlinear water waves. The numerical algorithm is tested on various wave problems, including a comparison with experiments containing wave interactions resulting in a large amplitude splash.


Introduction
The numerical simulation of nonlinear water waves is a challenging problem. These waves appear naturally in the ocean, rivers and lakes and greatly affect the motion of ships and induce significant forces on floating and fixed structures. Since in many cases the wave motion can be considered as nearly inviscid and irrotational, we model the water waves using the nonlinear potential flow water wave equations. Any amount of numerical dissipation, either added explicitly to stabilize the numerical scheme or implicitly present in the numerical discretization, will then significantly influence the accuracy of the wave computations, in particular for long time simulations.
The potential flow water wave equations, when expressed in terms of the free surface potential and wave height, have a Hamiltonian structure, but this structure is generally lost in the numerical discretization. The main topic of this article is to develop a finite element discretization for which we can explicitly prove that it preserves the Hamiltonian structure, even on time-dependent meshes that are needed to follow the free surface motion. This will directly result in an energy preserving numerical discretization that is stable for large amplitude nonlinear waves.
The starting point for the derivation of the Hamiltonian finite element discretization for nonlinear potential flow water waves is Luke's variational formulation [16]. After introducing time-dependent basis functions in the variational formulation, we can prove in several steps that the numerical discretization exactly preserves the Hamiltonian structure, even on time-dependently deforming unstructured meshes suitable for large amplitude waves. This Hamiltonian structure, when combined with a symplectic time integration scheme, prevents energy drift. A crucial subtlety in maintaining the Hamiltonian structure stems from the mesh movement. Any movement of the free surface needs to be accommodated by a mesh movement in the interior of the domain, which results in an intricate coupling between the free surface motion and the solution of the Laplace equation for the velocity potential.
Many numerical discretizations have been proposed for the solution of the potential flow water wave equations. The most popular approach for solving these equations is the boundary element method, starting with the work by Longuet-Higgins and Cokelet [15]. More recent works in this direction include [3,4,8,9,11,13], while the older works are covered by the review paper by Tsai and Yue [26]. These methods, however, typically require evaluating a singular integration kernel and tend to require the evaluation of dense matrix-vector products, which have to be solved with a fast multipole method to keep the computational complexity approximately linear. Moreover, these methods do not automatically preserve the Hamiltonian structure of the potential flow water wave equations.
An alternative approach is to use the finite element method, computing the solution in the entirety of the domain. This is not necessarily more expensive, since all interesting physical phenomena still happen at the free surface, allowing the use of a limited number of elements in the vertical direction. The finite element discretization gives rise to a sparse system of equations, meaning it is much easier to solve in linear time. However, all previous attempts to use a finite element discretization require additional numerical stabilization or specially constructed meshes to prevent numerical instabilities [17][18][19]21,24,25,[27][28][29][30]32].
The direct precursor to this work is the work by Gagarina et al. [10]. The main difference is that we prove that the discrete equations retain the Hamiltonian structure of the potential flow water wave equations, even on unstructured, time dependent meshes, and that we provide explicit equations for the dependence of the unknowns on general mesh movement, thereby generalizing the result in [10].
In the remainder of this article, we will introduce the potential flow equations and a suitable Lagrangian in Sect. 2. In Sect. 3 we will use this Lagrangian to construct a discrete Hamiltonian formulation. We will introduce a time stepping scheme for the resulting nonautonomous Hamiltonian system in Sect. 4 and discuss an efficient technique to solve the resulting algebraic equations. Finally, in Sect. 5 we present some results that numerically verify the stability and accuracy of the numerical scheme.

Governing Equations
The equations describing water wave motions are defined on a time-dependent domain Ω t ⊆ R 3 , t ∈ (0, T ). Its boundary ∂Ω t is split into two parts, ∂Ω t = ∂Ω t,s ∪ ∂Ω t,R , with ∂Ω t,s the free surface. The other boundary ∂Ω t,R may consist e.g. of a wave maker, a beach or a bottom surface. Each point R ∈ ∂Ω t,R has a prescribed velocity ∂ R ∂t . The position of the free surface ∂Ω t,s is unknown a priori and is to be determined as part of the initial-boundary-value problem describing wave motions.
We assume that the free surface is single-valued. This allows us to define with x, y, z the coordinates in a Cartesian system, with the undisturbed water surface equal to z = z 0 and η : (0, T ) × D t,s → R representing the free surface ∂Ω t,s . The dynamics of inviscid potential flow water waves is now governed by the following initial-boundary-value problem with initial conditions where the operators ∇ and Δ are, respectively, the gradient and Laplace operator, g = g x , g y , g z T ∈ R 3 the gravity vector, ν the unit outward normal vector at ∂Ω t , and φ is the velocity potential.

Lagrangian and Hamiltonian Approach
The variational formulation of (2) using dimensionless variables is based on the Lagrangian withx = x, y, z T , which was presented by Luke [16]. To compute the variations of this functional we use Reynolds transport theorem [7] in the form of the identity where v(t, x) = dx dt is the velocity of the domain boundary ∂Ω t . The continuum equations can readily be recovered by computing variations with respect to η and φ. That is, we require The free surface height function η only appears in the functional L 0 in the description of the free surface boundary (1), so if we compute the functional derivative using Leibniz' theorem [7], we obtain Since δη is arbitrary, this recovers (2c) when δ η L 0 = 0. Before computing variations with respect to φ we rewrite (3) using (4) where in the integrands over D t,s we have z = η. The variations of L 0 with respect to φ, after integration by parts, are equal to where the variations δ φ at t = 0 and t = T are taken to be zero and we used the relation ν| ∂Ω t,s = ∇(z−η) |∇(z−η)| . Considering the arbitrary variations in the interior and at the different sections of the domain boundary separately, this recovers the other three equations in (2) after setting δ φ L 0 = 0.
Moreover, the governing equations in (2) can be recognized as a Hamiltonian system with respect to the unknown free surface height and free surface potential [33]. Preserving this Hamiltonian structure in the finite element discretization significantly improves the accuracy of free surface wave computations, but currently there are only a few attempts to preserve the Hamiltonian structure in a discretization [5]. Another benefit of a Hamiltonian (Galerkin) semidiscretization is that we can make use of the well-developed geometric integration theory [12] to construct an energy-preserving numerical discretization.

Discretization of the Variational Principle
The finite element discretization is based on a tessellation T h of the domain Ω t . The tessellation T h is changing in time to accommodate for the free surface motion and other moving boundaries, such as a wave maker.
Using a nodal Lagrangian basis, the set of nodes in Ω t is denoted with N . Within this set the nodes on ∂Ω t,s are denoted with N s , those on ∂Ω t,R with N R , while the other nodes are denoted with N r . Note that N s ∩ N R is in general non-empty, these nodes correspond to grid points located at the interface of the free surface and the other boundaries.
Using the various sets of nodes, φ t j j∈N denotes the set of basis functions used to approximate the velocity potential φ and with a slight abuse of notation η t j j∈N s the basis functions for the free surface height η.
With these basis functions we approximate the free surface height and velocity potential, respectively, as The computational domain determined by the numerical approximations at time t is denoted with Ω t,h and the corresponding numerical approximation of the free surface and other boundary surfaces with ∂Ω t,s,h and ∂Ω t,R,h , respectively. It is of critical importance to note that the tessellation T h , and therefore also the basis functions, have an explicit dependency on both the free surface height η and the prescribed domain boundary movement of ∂Ω t,R . This also implies that the basis functions have a dependency on t and each of the expansion coefficients a of the free surface height approximation.
We require that the basis functions φ t j j∈N and η t j j∈N s satisfy the following compatibility condition at the free surface ∂Ω t,s,h At the domain boundaries ∂Ω t,s and ∂Ω t,R the discretization φ h can be given as With these approximations, the Lagrangian functional (3) can be (semi)-discretized as We will compute the variations separately for a j (t) and b j (t). To write the corresponding variational derivatives in a more compact form, we introduce the matrices Φ t where τ is the tangential vector along the interface ∂Ω t,s,h ∩∂Ω t,R,h , e z = (0, 0, 1) T , and the indices a and t denote explicit, but hidden dependencies on the free-surface parametrization and time. Furthermore, we use the following decompositions: where the submatrix where Φ 1,a,R is Φ a,R [i] with i ∈ N s . We use without further reference the following:

Proposition 1
The matrices and vectors introduced above satisfy the following statements.
For simplicity, we usually do not denote the dependence of Φ t on a and the timedependence of the submatrices Φ i j i, j = 1, 2. We will also use the notations a, b and b s for the vectors composed of the coefficients {a j } j∈N s , {b j } j∈N and {b j } j∈N s , respectively.

Lemma 1
The variational principle for the Lagrangian L h (6) with semidiscretized variables a and b, and t ∈ (0, T ), can be given in the following matrix-vector form Proof We compute variations of the discrete functional L h with respect to a j and b j , which are assumed to be zero for t = 0 and t = T . We consider the contributions in (6) one by one, starting with partial derivatives with respect to a j Applying Leibniz' theorem we obtain where we used the fact that in this integral only the free surface boundary depends on a j . The second term introduces an additional complication, since not only the boundary, but also the integrand depends on a j . Split the domain according to partition (1). We can now make the dependency of the boundary on η explicit by splitting the integral into two parts where we used (1) to represent the free surface with the height functions η l h . Next, we apply first Leibniz' theorem to the first term on the right hand side in (11) and then Reynolds transport theorem to the second term Since we assume that δφ(0, ·) = δφ(T, ·) = 0 the integrals over Ω T and Ω 0 vanish. We also use the relations It is beneficial to further evaluate the second term on the right hand side in (12) using (7.2) in Flanders [7], splitting this contribution into a line integral that will be used when applying (7.2) in [7] to the integrals over D t,s,h and an expression where ∂ ∂a j is the outermost operator, which is convenient when constructing the Hamiltonian. This results in with τ the unit tangential vector at the interface ∂Ω t,R,h ∩ ∂Ω t,s,h . Note that τ is orthogonal to ν R . The vector ∂ x ∂a j links the mesh velocity to the free surface velocity. Since the mesh is a tessellation of the domain, the mesh at ∂Ω t,R,h can only move parallel to ∂Ω t,R,h , hence the first integral on the right hand side in (13) is zero. At the solid wall-free surface intersection ∂Ω t,R,h ∩ ∂Ω t,s,h we need to enforce that the free surface moves tangentially to the solid wall, hence we have to apply the correction with x s,R = x, y, η h . An alternative interpretation of this is that an infinitesimally small sliver of D t,s,h nearest to the wave maker is rotated, such that it aligns correctly. By introducing t R = τ × ν R (14) can also be written as The second integral on the right hand side in (13) is then equal to After collecting all terms, the final result for (13) then becomes Using the compatibility condition at the free surface (5) we obtain then Since the basis functions η t k do not depend on a j the last integral in (15) is zero and ∂ a j M 2 can be expressed as For the third term we simply write Combining the three terms we obtain for the variations of L h with respect to a j which is equivalent to (9a). For Eqs. (9b) and (9c) consider the variations with respect to b j . The first term does not depend on b j , hence For the second term use (4) again to obtain Finally, the third term is just a straightforward differentiation. We have After applying the decompositions (7) and (8) the terms can be combined to give (9b) and (9c).
To express the discretized variational principle in (9) as a Hamiltonian system we introduce the variableb We will also use the notation S = Φ 11 − Φ 12 Φ −1 22 Φ 21 for the Schur complement, possibly with an upper index to indicate its dependence on time.
Lemma 2 Using the variableb s , the system in (9) can be recasted aṡ . Now focus on (9b) Expand the vector product which can be reordered to form (18a). The other equality requires a more elaborate derivation. Split (9a) into two parts. The first part of (9a) becomes using the relationĖ which can be derived with help of (7.2) in [7] in a manner similar to obtain (13). For details, see "Appendix 2". Now use the product rule and the definition (17) to rewrite (19) in terms ofb s to obtain Next, consider the second part of (9a) Here, in the second step the block structured matrix is expanded into its components. Expanding the products and adding up similar terms finally results in which can be combined with (20) to obtain the statement of the Lemma. (6) is equivalent with the forced Hamiltonian systemȧ

Theorem 1 The discrete variational form corresponding to the Lagrangian
where Proof Obviously, On the other hand

Remark 1
We note that the only explicit dependence on time in the discrete Hamiltonian is caused by the wave maker motion. Therefore this semi-discrete formulation is energy conservative, viz. dH dt = 0 without a wave-maker, even on unstructured meshes. This motivates to integrate (22)-(23) with a symplectic time integrator, since this will then result in a stable numerical discretization without the need for the addition of any stabilization terms.
In order to simplify the derivation of the time discretization, which will be discussed in the next section, we use the following relation where we have used the identity ∂ ∂t (24) greatly shortens expressions whenever E andb s are to be evaluated at the same time levels in the time integration method. This expansion appears to be redundant in view of (21). However, the interior component of b represents a solution of the Laplace equation. This could cause a dependency of b on the boundary shape, so we do need (24) to show that this dependency is fully factored in Φ.

Time Integration for the Discrete Variational Formulation
The time discretization of the Hamiltonian finite element discretization (22) is performed with the second order accurate Störmer-Verlet time integration method. The Hamiltonian system (22) is, however, non-autonomous. This requires a modification of the Störmer-Verlet scheme for which we follow the procedure outlined in [14]. Given the non-autonomous Hamiltonian system we introduce the new variables P = ( p, τ 0 ) and Q = (q, τ ) and the HamiltonianH with H (P, Q) = H (τ, p, q)−τ 0 . Here τ corresponds to time and the fictitious variable τ 0 ensures that P and Q are of the same dimension. The Störmer-Verlet scheme for a non-autonomous Hamiltonian system can now be expressed as which in the original variables gives The update for τ n ensures that τ indeed represents time and will be taken for granted in the following. The fictitious variable τ 0 is of no interest to us, so its update scheme is not presented.
The non-autonomous Störmer-Verlet scheme applied to the Hamiltonian finite element discretization (22) using p = a and q =b s in (25) is now equal to where relation (24) has been used to shorten the expressions. In terms of the original variables a, b and b s we obtain now the algebraic equations Since the time stepping in (26) is implicit, we first solve the equation for a n+ 1 2 with a Newton method, followed by the equation for b s,n+1 . Finally, a n+1 can be obtained in an explicit way. A full derivation that prepares (26) for numerical treatment can be found in Appendix 1.
The full numerical scheme can be summarized as follows: -Interpolate the initial surface data. For simulations using a wave maker, a still free water surface is used. A more detailed description of the mesh movement, which is done after the free surface or the wavemaker updates, is given at the end of Sect. 4.2.

Computing Derivatives ∂ a Φ and ∂ a Φ R
In the derivation of the discrete Hamiltonian the derivatives ∂ a Φ and ∂ a Φ R have been left untreated, since this was beneficial for arriving at Eq. (23). In this section we will discuss the computation of the derivatives with respect to the free surface coefficients a. Consider ∂ a Φ element-wise, as a summation on the finite element tessellation T h : where the shape of the element and the basis functions depend on the free surface η h , hence implicitly on the coefficients a. Introduce a reference elementK . We will denote the image of the basis functions onK asφ, the reference coordinates as (x,ŷ,ẑ) and the gradient operator with respect to reference coordinates as∇. Assume, for every element K ∈ T h , that there is an invertible mapping F K :K → K . Since we use nodal basis functions, we have x = lx lφl , wherex l are the coordinates of the nodal points of element K . The Jacobian of F K with respect to the reference element coordinates is given by J = lx l∇φl T .
Perform the coordinate transformation, Using the matrix identities ∂ ∂t where the summation convention is used on repeated indices. Transforming back to the elements K ∈ T h we obtain the relation The coupling between the node locations and the free surface parametrization ∂x l ∂a k has to be constructed depending on the choice of the mesh movement algorithm.
For the computation of ∂ a Φ R we can use (13). We have For the first term use the chain rule and the mapping F K and for the second term the compatibility condition (5) to find We would like to consider the mesh deformation and the rest of the time stepping scheme separately. To this end, we introduce the matrices and we obtain the relations The matrices C and B are separated into free surface and interior parts similar to Φ. The node velocity ∂x i ∂a j follows from the mesh movement algorithm. Assume that the mesh movement algorithm, with free surface node positions fixed, is either based on maintaining a force balance or based on solving an additional PDE, see Sect. 4.2. In both cases, the node displacements u are given by where F is the Jacobian with respect to the node displacements of the mesh movement algorithm. Inverting the Jacobian gives The node displacements and the node position are linked by a constant offset, hence we directly obtain

Mesh Deformation Algorithm
We base the mesh deformation algorithm on Masud and Hughes [20]. The idea is to compute a displacement field u ∈ R d and apply the computed displacements to the node coordinates. We use the still-water domain to provide an initial grid corresponding to the zero displacement field. The displacement field is the solution of the boundary value problem where τ is a bounded nonnegative function and the zero displacement domain Ω z is also used to compute the displacements. The free surface height η is the instantaneous wave height, hence (28) computes, contrary to [20], the displacements with respect to the original mesh for every update in the free surface height η. While this is given as a continuous system of equations, they are discretized using linear basis functions on T h in order to guarantee that the computed nodal displacements u can directly be used to deform the mesh. The parameter τ is typically large in areas where the elements are small, to prevent grid inversion. It is also large near the free surface to ensure that the gridlines closely follow the free surface and wave maker motion. To compute ∂x l ∂a k we take the derivative ∂ a of (28). Since (28) is solved on a fixed domain, there are no hidden dependencies and we can directly write with the understanding that these equations have to be discretized in the same way as the equations for the displacement. The derivative ∂x ∂t can be approximated in a similar manner. In our simulations, the small elements reside mostly near the free surface, so we choose τ = e 1+cz , where c ∼ = 1 can be tuned to prevent inversion for very shallow water simulations and simulations involving very steep waves or tuned to improve conditioning for very deep water simulations.
For more general problems the variable τ in element K can be computed as with Δ min , Δ max the area (or volume) of the smallest and largest elements in the mesh and Δ K the area (or volume) of element K . In [20] it is shown that this results in τ K -values that are essentially independent of the ratio Δ min /Δ max . A more detailed way to compute the τ -values is presented in [1], where the ratio of the inverse of the element Jacobian at the quadrature points to a reference quantity, e.g. the minimum of the inverse Jacobian in the mesh, is used to control the mesh deformation. This helps to ensure that the Jacobian remains positive inside the element, which prevents grid inversion. During the mesh updates we keep the background mesh fixed where we solve (28) with a conforming nodal finite element method, using a h and R(t) as inputs. Next, we reconstruct the mesh by displacing the actual nodes from the background nodes with the computed displacements.
This algorithm is a simplified form of the mesh deformation algorithms based on the elasticity equations. These algorithms are widely used for complex fluid-structure interaction problem and allow complex mesh deformations, see e.g. [1,2,22,23].

Fenton and Rienecker Wave
As a first model problem we consider the two-dimensional semi-analytical steady wave solution of (2), computed by Fenton and Rienecker [6] using a combination of Fourier expansions and numerical methods. This solution provides a correction to the Stokes wave [31]. The Fenton wave is a standard test case suitable to investigate the accuracy of numerical methods for nonlinear water waves. See Fig. 1 for an impression. We compute the steady wave  9.5 × 10 −4 1.2 × 10 −9 1.56 3.93 The time step is given as a fraction of the wave period solution for a water depth H = 1, gravity coefficient g = 1 and domain length X ∼ = 4.9636. The zero displacement grid is a regular grid of N x by N y rectangles, see Table 1. This grid was further split into triangles by subdividing the rectangles along their diagonals in an alternating manner. We performed a convergence test for linear basis functions by simulating a water wave for 10 periods, with 12N x 32 time steps per period, comparing the free surface height η h computed with the Hamiltonian finite element discretization with the free surface height computed with the semi-analytical method proposed by Fenton and Rienecker [6]. Since the focus of the Hamiltonian finite element method is on energy conservation, we also compute the difference between the minimum and the maximum Hamiltonian energy during the simulation. The results in Table 1 show that second order accuracy for the free surface height is obtained.
Next, we also performed a long simulation for 100 wave periods attempting to detect if there exists a systematic drift in the Hamiltonian energy. The results of the energy variation in this simulation can be found in Fig. 2.

Soliton
The second test case for code verification is provided by a traveling wave solution. The initial wave profile in a 2D domain of depth 0.5 is given by After moving away from the boundary, this initial wave profile will deform into a traveling wave closely resembling for some offset c. A close approximation of this solution is depicted in Fig. 4. This test case was also considered by Westhuis [30], who used a combined finite difference-finite element discretization of (2). The travelling wave will be simulated for 120 s, with Δt = 0.05. The domain has a reflecting wall at x = 150 m. In order to verify numerically the stability of our new scheme, we choose a sequence of mesh sizes Δx ∈ {2 m, 1 m, 0.5 m, 0.25 m}. In the vertical direction, we reproduce the choices made in [30]. That is, we use 6 elements in The coarsest of these meshes is unable to sufficiently resolve the traveling wave profile, while the choice Δx = 1 m can resolve the traveling wave, but not the high frequency modes that are required to keep the wave stable. In these cases, we cannot expect to solve the equations with any accuracy, but we still find reasonable bounds on the energy. See Fig. 3 for an overview of the behavior of the energy. Figure 4 shows snapshots of the wave profile for various meshes.
Since Δx is the only parameter changed in these computations we expect that any changes in the energy are caused by the nonlinear exchange of energy with under-resolved wave modes. This is confirmed by the dip in the energy when the wave interacts with the wall, where the high frequency modes play a larger role. Since we are only interested in the first part of the domain, we use a numerical basin of 90 m × 1 m, which is still large enough to ensure that no spurious reflections from the end wall interfere with the computed waves of interest. From the Fenton wave test case we know that rectangular elements offer greater accuracy, so we use rectangles near the surface. Further away from the surface we use triangles. Moreover, since we already know in advance that there will be localized phenomena near x = 50 m we refine the mesh in that area. We note that in practice 3 N iterations are usually enough. Snapshots of the mesh in the transition zone and in the splash zone are provided in Fig. 5. Following [10] we used Δx = 0.0027 near the splash zone and Δx = 0.016 away from the splash. Comparing the measured wave height to the computed wave height at various locations in the model basin, both in the time domain, see Fig. 6, and in the frequency domain, see Fig. 7, we conclude that they are in good agreement with experiments.
with Δb (k) s,n+1 and where a (0) n+ 1 2 = a n and b s,n+1 = b s,n . The non-linear algebraic Eqs. (29) and (30) are iterated until convergence is reached. In the Jacobian in (29) we use that ∂ a i and ∂ b s, j commute. Recall that S = Φ 11 − Φ 12 Φ −1 22 Φ 21 . For numerical efficiency reasons it is crucial to avoid explicitly forming Φ −1 22 . The following auxiliary equation is therefore introduced This introduces the inverse matrix F −1 22 for which we introduce the auxiliary equation We can now write (29) in a form that does not require explicit inverses of matrices where Δu We again introduce an auxiliary equation s,n+1 +

Appendix 2: Time Derivative of the Mass Matrix
The time derivative of the mass matrix can be constructed similar to the free surface derivative of the wave maker boundary integral. This matrix is introduced in (10), (11) and (16). A careful look at these equations reveals that the mass matrix is more accurately represented as a flux trough a surface Since η t j does not depend on z, the first term in (36) vanishes. In addition, the internal contributions from two adjacent patches D l t,s,h cancel, so the third term reduces to an integral over the boundary of the free surface. At the boundary we know the patch velocity v = ∂ R ∂t and we have applied the correction (14) to ensure that the free surface moves tangentially to the solid wall. These considerations reduce (36) From here we use t R = τ × ν R and expand all products to arrive at