Hamiltonian Finite Element Discretization for Nonlinear Free Surface Water Waves

1840 Accesses
7 Citations
Explore all metrics

Abstract

A novel finite element discretization for nonlinear potential flow water waves is presented. Starting from Luke’s Lagrangian formulation we prove that an appropriate finite element discretization preserves the Hamiltonian structure of the potential flow water wave equations, even on general time-dependent, deforming and unstructured meshes. For the time-integration we use a modified 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.

The Pressure Boundary Condition and the Pressure as Lagrangian for Water Waves

Article Open access 11 March 2019

Computation of Nonlinear Free-Surface Flows Using the Method of Fundamental Solutions

Free-Boundary Problems for Wave–Structure Interactions in Shallow-Water: DG-ALE Description and Local Subcell Correction

Article 17 January 2024

Use our pre-submission checklist

Avoid common mistakes on your manuscript.

1 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 non-autonomous 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.

2 Governing Equations

The equations describing water wave motions are defined on a time-dependent domain ${\varOmega }_t\subseteq {\mathbb {R}}^3, t\in (0,T)$. Its boundary $\partial {\varOmega }_t$ is split into two parts, $\partial {\varOmega }_t=\partial {\varOmega }_{t,s}\cup \partial {\varOmega }_{t,R}$, with $\partial {\varOmega }_{t,s}$ the free surface. The other boundary $\partial {\varOmega }_{t,R}$ may consist e.g. of a wave maker, a beach or a bottom surface. Each point $R\in \partial {\varOmega }_{t,R}$ has a prescribed velocity $\frac{\partial R}{\partial t}$. The position of the free surface $\partial {\varOmega }_{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

$$\begin{aligned} {\varOmega }_{t,\eta }=\left\{ (x, y, z)\in {\mathbb {R}}^3|(x, y)\in \mathcal {D}_{t,s}\subseteq {\mathbb {R}}^2, 0< z < \eta (t,x, y)\right\} , \end{aligned}$$

(1)

with x, y, z the coordinates in a Cartesian system, with the undisturbed water surface equal to $z=z_0$ and $\eta \,{:}\,(0,T)\times \mathcal {D}_{t,s}\rightarrow {\mathbb {R}}$ representing the free surface $\partial {\varOmega }_{t,s}$.

The dynamics of inviscid potential flow water waves is now governed by the following initial-boundary-value problem

where the operators $\nabla $ and ${\varDelta }$ are, respectively, the gradient and Laplace operator, $g{}=\begin{pmatrix}g_x,&g_y,&g_z\end{pmatrix}^T\in {\mathbb {R}}^3$ the gravity vector, $\nu $ the unit outward normal vector at $\partial {\varOmega }_t$, and $\phi $ is the velocity potential.

2.1 Lagrangian and Hamiltonian Approach

The variational formulation of (2) using dimensionless variables is based on the Lagrangian

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_0(\phi ,\eta )&= - \int _0^T \int _{{\varOmega }_t} g\cdot {\bar{x}} + \frac{\partial \phi }{\partial t} + \frac{1}{2}|\nabla \phi |^2\,\mathrm {d}{\varOmega }\,\mathrm {d}t, \end{aligned} \end{aligned}$$

(3)

with ${\bar{x}}=\begin{pmatrix}x,&y,&z\end{pmatrix}^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

$$\begin{aligned} \begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\int _{{\varOmega }_{t}}\phi \,\mathrm {d}{\varOmega }&=\int _{{\varOmega }_t}\frac{\partial \phi }{\partial t}\,\mathrm {d}{\varOmega }+ \int _{\partial {\varOmega }_t}\phi {v}\cdot \nu \,\mathrm {d}S \\&=\int _{{\varOmega }_{t}}\frac{\partial \phi }{\partial t}\,\mathrm {d}{\varOmega }+\int _{\partial {\varOmega }_{t,R}}\frac{\partial R}{\partial t}\cdot \nu \phi \,\mathrm {d}S+\int _{D_{t,s}}\frac{\partial \eta }{\partial t}\phi \,\mathrm {d}x\,\mathrm {d}y, \end{aligned} \end{aligned}$$

(4)

where ${v}(t,x)=\frac{\mathrm {d}{\bar{x}}}{\mathrm {d}t}$ is the velocity of the domain boundary $\partial {\varOmega }_t$. The continuum equations can readily be recovered by computing variations with respect to $\eta $ and $\phi $. That is, we require

$$\begin{aligned} 0=\delta _{\eta }{\mathcal {L}}_0=\lim _{\epsilon \downarrow 0}\frac{{\mathcal {L}}_0(\phi , \eta +\epsilon \delta \eta )-{\mathcal {L}}_0(\phi , \eta )}{\epsilon } \end{aligned}$$

and

$$\begin{aligned} 0=\delta _{\phi }{\mathcal {L}}_0=\lim _{\epsilon \downarrow 0}\frac{{\mathcal {L}}_0(\phi +\epsilon \delta \phi , \eta )-{\mathcal {L}}_0(\phi , \eta )}{\epsilon }. \end{aligned}$$

The free surface height function $\eta $ only appears in the functional ${\mathcal {L}}_0$ in the description of the free surface boundary (1), so if we compute the functional derivative using Leibniz’ theorem [7], we obtain

$$\begin{aligned} \delta _{\eta }{\mathcal {L}}_0=-\int _0^T\int _{\mathcal {D}_{t,s}}\big (g\cdot \begin{pmatrix} x,&y,&\eta \end{pmatrix} + \partial _t\phi +\frac{1}{2}|\nabla \phi |^2\big )\delta \eta \,\mathrm {d}x\,\mathrm {d}y\,\mathrm {d}t. \end{aligned}$$

Since $\delta {\eta }$ is arbitrary, this recovers (2c) when $\delta _{\eta }{\mathcal {L}}_0=0$. Before computing variations with respect to $\phi $ we rewrite (3) using (4)

$$\begin{aligned} {\mathcal {L}}_0(\phi ,\eta )&= -\int _0^T \int _{{\varOmega }_t} g\cdot {\bar{x}} + \frac{\partial \phi }{\partial t} + \frac{1}{2}|\nabla \phi |^2\,\mathrm {d}{\varOmega }\,\mathrm {d}t\\&=-\int _0^T \int _{{\varOmega }_t} g\cdot {\bar{x}} + \frac{1}{2}|\nabla \phi |^2\,\mathrm {d}{\varOmega }\,\mathrm {d}t + \int _{{\varOmega }_0}\phi (0,\cdot )\,\mathrm {d}{\varOmega }- \int _{{\varOmega }_T}\phi (T,\cdot )\,\mathrm {d}{\varOmega }\\&\quad +\,\int _0^T\int _{\partial {\varOmega }_{t,R}}\frac{\partial R}{\partial t}\cdot \nu \phi \,\mathrm {d}S\,\mathrm {d}t + \int _0^T\int _{\mathcal {D}_{t,s}}\frac{\partial \eta }{\partial t}\phi \,\mathrm {d}x\,\mathrm {d}y\,\mathrm {d}t, \end{aligned}$$

where in the integrands over $\mathcal {D}_{t,s}$ we have $ z=\eta $. The variations of ${\mathcal {L}}_0$ with respect to $\phi $, after integration by parts, are equal to

$$\begin{aligned} \delta _\phi {\mathcal {L}}_0&=\int _0^T \int _{{\varOmega }_t}{\varDelta }\phi \delta \phi \,\mathrm {d}{\varOmega }\,\mathrm {d}t\\&\quad +\int _0^T\int _{\partial {\varOmega }_{t,R}}\delta \phi \left( \frac{\partial R}{\partial t}\cdot \nu -\nabla \phi \cdot \nu \right) \,\mathrm {d}S\,\mathrm {d}t \\&\quad + \int _0^T\int _{\mathcal {D}_{t,s}}\delta \phi \left( \frac{\partial \eta }{\partial t}-\nabla \phi \cdot \begin{pmatrix}-\frac{\partial \eta }{\partial x},&-\frac{\partial \eta }{\partial y},&1\end{pmatrix}\right) \,\mathrm {d}x\,\mathrm {d}y\,\mathrm {d}t, \end{aligned}$$

where the variations $\delta _\phi $ at $t=0$ and $t=T$ are taken to be zero and we used the relation $\nu |_{\partial {\varOmega }_{t,s}}=\frac{\nabla (z-\eta )}{|\nabla (z-\eta )|}$. 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 $\delta _\phi {\mathcal {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.

3 Discretization of the Variational Principle

The finite element discretization is based on a tessellation $\mathcal {T}_h$ of the domain ${\varOmega }_t$. The tessellation ${\mathcal {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 ${\overline{{\varOmega }}}_t$ is denoted with ${\mathcal {N}}$. Within this set the nodes on $\partial {\varOmega }_{t,s}$ are denoted with $N_s$, those on $\partial {\varOmega }_{t,R}$ with $N_R$, while the other nodes are denoted with $N_r$. Note that $N_s\cap 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, $\left\{ \phi _j^t\right\} _{j\in {\mathcal {N}}}$ denotes the set of basis functions used to approximate the velocity potential $\phi $ and with a slight abuse of notation $\left\{ \eta _j^t\right\} _{j\in N_s}$ the basis functions for the free surface height $\eta $.

With these basis functions we approximate the free surface height and velocity potential, respectively, as

$$\begin{aligned} \eta (t,x,y) \cong \eta _h(t,x,y) = \sum _{j\in N_s} a_j(t)\eta _j^t(x,y), \end{aligned}$$

and

$$\begin{aligned} \phi (t,x,y,z) \cong \phi _h(t,x,y,z) = \sum _{j\in {\mathcal {N}}} b_j(t)\phi _j^t(x,y,z). \end{aligned}$$

The computational domain determined by the numerical approximations at time t is denoted with ${\varOmega }_{t,h}$ and the corresponding numerical approximation of the free surface and other boundary surfaces with $\partial {\varOmega }_{t,s,h}$ and $\partial {\varOmega }_{t,R,h}$, respectively.

It is of critical importance to note that the tessellation ${\mathcal {T}}_h$, and therefore also the basis functions, have an explicit dependency on both the free surface height $\eta $ and the prescribed domain boundary movement of $\partial {\varOmega }_{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 $\left\{ \phi _j^t\right\} _{j\in {\mathcal {N}}}$ and $\left\{ \eta _j^t\right\} _{j\in N_s}$ satisfy the following compatibility condition at the free surface $\partial {\varOmega }_{t,s,h}$

$$\begin{aligned} \eta _j^t(x,y) = \phi _j^t(x,y,z)\quad \text {for}\;(x,y,z)\in \partial {\varOmega }_{t,s,h}, j\in N_s. \end{aligned}$$

(5)

At the domain boundaries $\partial {\varOmega }_{t,s}$ and $\partial {\varOmega }_{t,R}$ the discretization $\phi _h$ can be given as

$$\begin{aligned} \phi _h(t,x,y,z) = \sum _{j\in N_s} b_j(t)\phi _ j^t(x,y,z) \quad \text {for}\;(x,y,z) \in \partial {\varOmega }_{t,s,h} \end{aligned}$$

and

$$\begin{aligned} \phi _h(t,x,y,z) = \sum _{j\in N_R} b_j(t)\phi _ j^t(x,y,z) \quad \text {for}\;(x,y,z) \in \partial {\varOmega }_{t,R,h}. \end{aligned}$$

With these approximations, the Lagrangian functional (3) can be (semi)-discretized as

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_h(\phi _{h},\eta _h)&=- \int _0^T \int _{{\varOmega }_{t,h}} g\cdot {\bar{x}}+\frac{\partial }{\partial t}\left( \sum _{j\in {\mathcal {N}}}b_j\phi _j^t\right) +\frac{1}{2} \left| \sum _{j\in {\mathcal {N}}} b_j(t)\nabla \phi _j^t\right| ^2\,\mathrm {d}{\varOmega }\,\mathrm {d}t\\&:= \int _0^T M_1(t) + M_2(t) + M_3(t) \,\mathrm {d}t. \end{aligned} \end{aligned}$$

(6)

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 ${\varPhi }^t_{a}\in {\mathbb {R}}^{|{\mathcal {N}}|\times |{\mathcal {N}}|}, \mathcal {E}^t\in {\mathbb {R}}^{|N_s|\times |N_s|}, \mathcal {E}_{R}^t\in {\mathbb {R}}^{|N_s|\times |N_s|}$ and $\mathcal {D}^t\in {\mathbb {R}}^{|N_s|\times |N_s|}$ and the vectors ${\varPhi }_{a,R}^t\in {\mathbb {R}}^{|{\mathcal {N}}|}$ and $G\in {\mathbb {R}}^{|N_s|}$, which are defined as

$$\begin{aligned} {\varPhi }^t_{{a}} [i,j]= & {} \int _{{\varOmega }_{t,h}} \nabla \phi _i^t\cdot \nabla \phi _j^t\,\mathrm {d}{\varOmega }, \quad \mathcal {E}^t[i,j] = \int _{\mathcal {D}_{t,s,h}} \eta _i^t \eta _j^t \,\mathrm {d}x\,\mathrm {d}y,\\ \mathcal {D}^t[i,j]= & {} \int _{\mathcal {D}_{t,s,h}}\eta _i^t \frac{\partial \eta _j^t}{\partial t} \,\mathrm {d}x\,\mathrm {d}y,\\ \mathcal {E}_{R}^t [i,j]= & {} \int _{\partial {\varOmega }_{t,s,h}\cap \partial {\varOmega }_{t,R,h}} \eta ^t_i \eta ^t_j \frac{\partial R}{\partial t}\cdot \nu _R((\tau \times \nu _R)\cdot e_z)\,\mathrm {d}l \end{aligned}$$

and

$$\begin{aligned} {\varPhi }_{{a},R}^t [i] = \int _{\partial {\varOmega }_{t,R,h}} \frac{\partial R}{\partial t}\cdot \nu \phi _i^t \,\mathrm {d}S,\quad G^t[i]=\int _{\mathcal {D}_{t,s,h}}g\cdot \begin{pmatrix} x,&y,&0\end{pmatrix}\eta _i^t\,\mathrm {d}x\,\mathrm {d}y, \end{aligned}$$

where $\tau $ is the tangential vector along the interface $\partial {\varOmega }_{t,s,h}\cap \partial {\varOmega }_{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:

$$\begin{aligned} {\varPhi }^t_{a} = \begin{pmatrix} {\varPhi }_{11}&{}{\varPhi }_{12}\\ {\varPhi }_{21}&{}{\varPhi }_{22} \end{pmatrix}, \end{aligned}$$

(7)

where the submatrix ${\varPhi }_{11}\in {\mathbb {R}}^{|N_s|\times |N_s|}$ is ${\varPhi }_a^t[i,j]$ with $i,j\in N_s$ and ${\varPhi }_{{a},R}$ is split into

$$\begin{aligned} {\varPhi }_{{a},R} = \begin{pmatrix} {\varPhi }_{1,{a},R}^T&{\varPhi }_{2,{a},R}^T\end{pmatrix}^T, \end{aligned}$$

(8)

where ${\varPhi }_{1,{a},R}$ is ${\varPhi }_{a,R}[i]$ with $i\in N_s$. We use without further reference the following:

Proposition 1

The matrices and vectors introduced above satisfy the following statements.

(i)
Both ${\varPhi }^t_{a}$ and $\mathcal {E}^t$ are symmetric, ${\varPhi }^t_{a}$ is positive semidefinite, with ${\varPhi }_{21} = {\varPhi }_{12}^T$, and $\mathcal {E}^t$ is positive definite.
(ii)
The matrix $\mathcal {E}_{R}^t$ has only non-zero elements if $i,j\in N_s\cap N_R$, hence this matrix has a non-zero block of size $|N_s \cap N_R|\times |N_s\cap N_R|$, the remaining elements are zero.
(iii)
The non-zero components of ${\varPhi }_{{a},R}^t [i]$ are those with $i\in N_R$.

For simplicity, we usually do not denote the dependence of ${\varPhi }^t$ on ${a}$ and the time-dependence of the submatrices ${\varPhi }_{ij} \; 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\in N_s}, \{b_j\}_{j\in {\mathcal {N}}}$ and $\{b_j\}_{j\in N_s}$, respectively.

Lemma 1

The variational principle for the Lagrangian ${\mathcal {L}}_h$ (6) with semidiscretized variables ${a}$ and ${b}$, and $t\in (0,T)$, can be given in the following matrix-vector form

$$\begin{aligned} 0&=\mathcal {E}{{\dot{{b}}}}_s(t)+\mathcal {D}{b}_s(t)+\mathcal {E}_R{b}_s(t)+\mathcal {E}{a}(t)g_z\nonumber \\&\quad +\,\partial _{a}\left( \frac{1}{2}{b}(t)^T{\varPhi }{b}(t)-{\varPhi }_R\cdot {b}(t)\right) +G, \end{aligned}$$

(9a)

$$\begin{aligned} 0&=\mathcal {E}{{\dot{{a}}}}(t)+\mathcal {D}{a}(t)-\left[ {\varPhi }_{11}\,{\varPhi }_{12}\right] {b}(t)+{\varPhi }_{1,{a},R}, \end{aligned}$$

(9b)

$$\begin{aligned} 0&=-\left[ {\varPhi }_{21}\,{\varPhi }_{22}\right] {b}(t)+{\varPhi }_{2,{a},R}. \end{aligned}$$

(9c)

Here, an overdot represents differentiation with respect to time.

Proof

We compute variations of the discrete functional ${\mathcal {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$

$$\begin{aligned} \partial _{a_j}M_1 = -\partial _{a_j}\int _{{\varOmega }_{t,h}} g\cdot {\bar{x}}\,\mathrm {d}{\varOmega }= -\partial _{a_j}\int _{\mathcal {D}_{t,s,h}}\int _{ z=0}^{\eta _h}g\cdot ( x, y, z)\,\mathrm {d}z\,\mathrm {d}x\,\mathrm {d}y. \end{aligned}$$

Applying Leibniz’ theorem we obtain

$$\begin{aligned} \partial _{a_j}M_1&= -\int _{\mathcal {D}_{t,s,h}}g\cdot \begin{pmatrix} x,&y,&\eta _h\end{pmatrix}\eta _j^t\,\mathrm {d}x\,\mathrm {d}y =-g_z(\mathcal {E}^t{a})[j] - G^t[j], \end{aligned}$$

(10)

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 $\eta $ explicit by splitting the integral into two parts

$$\begin{aligned} \begin{aligned} \int _0^T\partial _{a_j} M_2(t)\,\mathrm {d}t&=-\int _0^T \partial _{a_j}\int _{\mathcal {D}_{t,s,h}}\int _{ z=0}^{\eta _h}\frac{\partial \phi _h^t}{\partial t}( x, y, z)\,\mathrm {d}z\,\mathrm {d}x\,\mathrm {d}y\,\mathrm {d}t\\&\quad -\int _0^T\int _{{\varOmega }_{t,h}\setminus {\varOmega }_{t,\eta }}\frac{\partial }{\partial a_j}\left( \frac{\partial \phi _h^t}{\partial t}\right) \,\mathrm {d}{\varOmega }\,\mathrm {d}t, \end{aligned} \end{aligned}$$

(11)

where we used (1) to represent the free surface with the height functions $\eta _h^l$. 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

$$\begin{aligned} \int _0^T\partial _{a_j} M_2(t)\,\mathrm {d}t&=-\int _0^T\int _{\mathcal {D}_{t,s,h}}\frac{\partial \phi _h^t}{\partial t}( x, y, \eta _h)\eta _j^t( x, y)\,\mathrm {d}x\,\mathrm {d}y\,\mathrm {d}t\\&\quad -\int _{{\varOmega }_T}\frac{\partial \phi ^T}{\partial a_j}\,\mathrm {d}{\varOmega }+\int _{{\varOmega }_0}\frac{\partial \phi ^0}{\partial a_j}\,\mathrm {d}{\varOmega }+\int _0^T\int _{\partial {\varOmega }_{t,h}}\frac{\partial \phi _h^t}{\partial a_j}v\cdot \nu \,\mathrm {d}S\,\mathrm {d}t, \end{aligned}$$

with $v=\frac{\mathrm {d}x}{\mathrm {d}t}|_{\partial {\varOmega }_t}$. Since we assume that $\delta \phi (0,\cdot )=\delta \phi (T,\cdot )=0$ the integrals over ${\varOmega }_T$ and ${\varOmega }_0$ vanish. We also use the relations

$$\begin{aligned} v\cdot \nu |_{\partial {\varOmega }_{t,s,h}}\,\mathrm {d}S&=\frac{\mathrm {d}}{\mathrm {d}t}( x, y,\eta _h)\cdot \left( -\frac{\partial \eta _h}{\partial x}, -\frac{\partial \eta _h}{\partial y}, 1\right) \,\mathrm {d}x\,\mathrm {d}y=\frac{\partial \eta _h}{\partial t}\,\mathrm {d}x\,\mathrm {d}y \end{aligned}$$

and

$$\begin{aligned} v\cdot \nu |_{\partial {\varOmega }_{t,R,h}}&=\frac{\partial R}{\partial t}\cdot \nu _R, \end{aligned}$$

with $\nu _R=\nu |_{\partial {\varOmega }_{t,R,h}}$, resulting in

$$\begin{aligned} \int _0^T\partial _{a_j} M_2(t)\,\mathrm {d}t&=-\int _0^T\int _{\mathcal {D}_{t,s,h}}\frac{\partial \phi _h^t}{\partial t}( x, y, \eta _h)\eta _j^t( x, y)\,\mathrm {d}x\,\mathrm {d}y\,\mathrm {d}t\nonumber \\&\quad +\int _0^T\int _{\partial {\varOmega }_{t,R,h}}\frac{\partial \phi _h^t}{\partial a_j}\frac{\partial R}{\partial t}\cdot \nu _R\,\mathrm {d}S\,\mathrm {d}t\nonumber \\&\quad +\int _0^T\int _{\mathcal {D}_{t,s,h}}\frac{\partial \phi _h^t}{\partial a_j}( x, y, \eta _h)\frac{\partial \eta _h}{\partial t}(t, x, y)\,\mathrm {d}x\,\mathrm {d}y\,\mathrm {d}t. \end{aligned}$$

(12)

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 $\mathcal {D}_{t,s,h}$ and an expression where $\frac{\partial }{\partial a_j}$ is the outermost operator, which is convenient when constructing the Hamiltonian. This results in

$$\begin{aligned} \int _{\partial {\varOmega }_{t, R, h}}\frac{\partial \phi _h}{\partial a_j}\frac{\partial R}{\partial t}\cdot \nu _R\,\mathrm {d}S&=-\int _{\partial {\varOmega }_{t,R,h}}\nabla \cdot \left( \phi _h\frac{\partial R}{\partial t}\right) \frac{\partial x}{\partial a_j}\cdot \nu _R\,\mathrm {d}S\nonumber \\&\quad +\int _{\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}}\left( \frac{\partial x}{\partial a_j}\times \phi _h\frac{\partial R}{\partial t}\right) \cdot \tau \,\mathrm {d}l\\&\quad +\frac{\partial }{\partial a_j}\int _{\partial {\varOmega }_{t,R,h}}\phi _h\frac{\partial R}{\partial t}\cdot \nu _R\,\mathrm {d}S,\nonumber \end{aligned}$$

(13)

with $\tau $ the unit tangential vector at the interface $\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}$. Note that $\tau $ is orthogonal to $\nu _R$. The vector $\frac{\partial x}{\partial a_j}$ links the mesh velocity to the free surface velocity. Since the mesh is a tessellation of the domain, the mesh at $\partial {\varOmega }_{t,R,h}$ can only move parallel to $\partial {\varOmega }_{t,R,h}$, hence the first integral on the right hand side in (13) is zero. At the solid wall-free surface intersection $\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}$ we need to enforce that the free surface moves tangentially to the solid wall, hence we have to apply the correction

$$\begin{aligned} \frac{\partial x}{\partial a_j}=\frac{\partial x_{s,R}}{\partial a_j}-\left( \frac{\partial x_{s,R}}{\partial a_j}\cdot \nu _R\right) \nu _R, \end{aligned}$$

(14)

with $x_{s,R}=\begin{pmatrix}x,&y,&\eta _h\end{pmatrix}$. An alternative interpretation of this is that an infinitesimally small sliver of $\mathcal {D}_{t,s,h}$ nearest to the wave maker is rotated, such that it aligns correctly. By introducing $t_R=\tau \times \nu _R$ (14) can also be written as

$$\begin{aligned} \frac{\partial x}{\partial a_j}=\left( \frac{\partial x_{s,R}}{\partial a_j}\cdot \tau \right) \tau +\left( \frac{\partial x_{s,R}}{\partial a_j}\cdot t_R\right) t_R. \end{aligned}$$

The second integral on the right hand side in (13) is then equal to

$$\begin{aligned}&\int _{\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}}\left( \frac{\partial x}{\partial a_j}\times \phi _h\frac{\partial R}{\partial t}\right) \cdot \tau \,\mathrm {d}l\\&\quad = \int _{\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}} \left( \tau \times \left( \left( \frac{\partial x_{s,R}}{\partial a_j}\cdot \tau \right) \tau +\left( \frac{\partial x_{s,R}}{\partial a_j}\cdot t_R\right) t_R\right) \right) \cdot \phi _h\frac{\partial R}{\partial t}\,\mathrm {d}l\\&\quad =\int _{\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}}-\left( \frac{\partial x_{s,R}}{\partial a_j}\cdot t_R\right) \left( \nu _R\cdot \frac{\partial R}{\partial t}\right) \phi _h\,\mathrm {d}l. \end{aligned}$$

Since $\frac{\partial x_{s,R}}{\partial a_j}=\begin{pmatrix}0,&0,&\eta _j^t\end{pmatrix}$ we obtain $t_R\cdot \frac{\partial x_{s,R}}{\partial a_j}=(\tau \times \nu _R)\cdot e_z\eta _j^t$ with $e_z=(0,0,1)^T$. After collecting all terms, the final result for (13) then becomes

$$\begin{aligned} \int _{\partial {\varOmega }_{t, R, h}}\frac{\partial \phi _h}{\partial a_j}\frac{\partial R}{\partial t}\cdot \nu _R\,\mathrm {d}S&=-\int _{\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}}\left( \frac{\partial R}{\partial t}\cdot \nu _R\right) \left( \tau \times \nu _R\right) \cdot e_z\eta _j^t\phi _h^t\,\mathrm {d}l\\&\quad +\frac{\partial }{\partial a_j}\int _{\partial {\varOmega }_{t,R,h}}\phi _h\frac{\partial R}{\partial t}\cdot \nu _R\,\mathrm {d}S. \end{aligned}$$

Using the compatibility condition at the free surface (5) we obtain then

$$\begin{aligned} \begin{aligned} \partial _{a_j}M_2&= -\int _{\mathcal {D}_{t,s,h}}\partial _t\left( \sum _{k\in N_s}b_{k}\eta _k\right) \eta _j\,\mathrm {d}x\,\mathrm {d}y \\&\quad + \partial _{a_j}\int _{\partial {\varOmega }_{t,R,h}}\frac{\partial R}{\partial t}\cdot \nu _R\sum _{i\in {\mathcal {N}}}b_i\phi _i\,\mathrm {d}S\\&\quad -\int _{\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}}\left( \frac{\partial R}{\partial t}\cdot \nu _R\right) \left( \tau \times \nu _R\right) \cdot e_z\eta _j\sum _{k\in {\mathcal {N}}}b_k\eta _k\,\mathrm {d}l\\&\quad +\int _{\mathcal {D}_{t,s,h}}\frac{\partial }{\partial a_j}\left( \sum _{k\in N_s}b_k\eta _k^t\right) \frac{\partial \eta _h}{\partial t}\,\mathrm {d}x\,\mathrm {d}y. \end{aligned} \end{aligned}$$

(15)

Since the basis functions $\eta _k^t$ do not depend on $a_j$ the last integral in (15) is zero and $\partial _{a_j} M_2$ can be expressed as

$$\begin{aligned} \partial _{a_j}M_2=\left( -\mathcal {E}{\dot{{b}}}_s-\mathcal {D}{b}_s-\mathcal {E}_R{b}_s\right) [j]+\partial _{a_j}({\varPhi }_R\cdot {b}). \end{aligned}$$

For the third term we simply write

$$\begin{aligned} \partial _{a_j}M_3 = -\partial _{a_j}\int _{{\varOmega }_{t,h}}\frac{1}{2} \left| \sum _{i\in {\mathcal {N}}}b_i\nabla \phi _i^t\right| ^2 \,\mathrm {d}{\varOmega }= -\partial _{a_j}\left( \frac{1}{2}{b}^T{\varPhi }{b}\right) . \end{aligned}$$

Combining the three terms we obtain for the variations of ${\mathcal {L}}_h$ with respect to $a_j$

$$\begin{aligned} 0&=-\mathcal {E}{{\dot{{b}}}}_s(t)-\mathcal {D}{b}_s(t)-\mathcal {E}_R{b}_s(t)-\mathcal {E}{a}(t)g_z\\&\quad -\partial _{a}\left( \frac{1}{2}{b}(t)^T{\varPhi }{b}(t)-{\varPhi }_R\cdot {b}(t)\right) - G, \end{aligned}$$

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

$$\begin{aligned} \partial _{b_j}M_1=0. \end{aligned}$$

For the second term use (4) again to obtain

$$\begin{aligned} \partial _{b_j}M_2&=\partial _{b_j}\left( \int _{\partial {\varOmega }_{t,R,h}}\frac{\partial R}{\partial t}\cdot \nu \sum _{i\in N_R}b_i\phi _i\,\mathrm {d}S\nonumber \right. \\&\quad \left. +\,\int _{\mathcal {D}_{t,s,h}}\partial _t\left( \sum _{k\in N_s} a_k\eta _k\right) \sum _{i\in N_s}b_i\phi _i\,\mathrm {d}x\,\mathrm {d}y\right) \nonumber \\&=\int _{\partial {\varOmega }_{t,R,h}}\frac{\partial R}{\partial t}\cdot \nu \phi _j\,\mathrm {d}S+\int _{\mathcal {D}_{t,s,h}}\partial _t\left( \sum _{k\in N_s} a_k\eta _k\right) \eta _j\,\mathrm {d}x\,\mathrm {d}y\nonumber \\&={\left\{ \begin{array}{ll}({\varPhi }_{1,{a},R}+\mathcal {E}{{\dot{{a}}}}+\mathcal {D}{a})[j]&{}j\in N_s,\\ {\varPhi }_{2,{a},R}[j]&{}j\notin N_s. \end{array}\right. } \end{aligned}$$

(16)

Finally, the third term is just a straightforward differentiation. We have

$$\begin{aligned} \partial _{b_j}M_3=-\int _{{\varOmega }_{t,h}}\sum _{i\in {\mathcal {N}}}b_i\nabla \phi _i\cdot \nabla \phi _j\,\mathrm {d}{\varOmega }=-({\varPhi }{b})[j]. \end{aligned}$$

After applying the decompositions (7) and (8) the terms can be combined to give (9b) and (9c). $\square $

To express the discretized variational principle in (1) as a Hamiltonian system we introduce the variable

$$\begin{aligned} {\tilde{{b}}}_s = \mathcal {E}{b}_s. \end{aligned}$$

(17)

We will also use the notation ${S}= {\varPhi }_{11}-{\varPhi }_{12}{\varPhi }_{22}^{-1}{\varPhi }_{21}$ for the Schur complement, possibly with an upper index to indicate its dependence on time.

Lemma 2

Using the variable $\tilde{{b}}_s$, the system in (1) can be recasted as

$$\begin{aligned} \dot{{a}}(t)= & {} \mathcal {E}^{-1}\left( -\mathcal {D}{a}(t)+ {S}\mathcal {E}^{-1} \tilde{b}_s - {\varPhi }_{1,{a},R}+ {\varPhi }_{12} {\varPhi }_{22}^{-1} {\varPhi }_{2,{a},R}\right) ,\quad \end{aligned}$$

(18a)

$$\begin{aligned} \dot{\tilde{{b}}}_s(t)= & {} \mathcal {D}^T\mathcal {E}^{-1}\tilde{{b}}_s(t)-\mathcal {E}{a}(t)g_z-\partial _{{a}}\left( \frac{1}{2}\tilde{{b}}_s(t)^T\mathcal {E}^{-1}{S}\mathcal {E}^{-1}\tilde{{b}}_s(t)\right. \nonumber \\&\left. -\,{\varPhi }_{1,{a},R}^T\mathcal {E}^{-1}\tilde{{b}}_s+{\varPhi }_{2,{a},R}{\varPhi }_{22}^{-1}{\varPhi }_{21}\mathcal {E}^{-1}\tilde{{b}}_s-\frac{1}{2}{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{2,{a},R}\right) \nonumber \\&-\,G. \end{aligned}$$

(18b)

Proof

First, we note that (9c) immediately implies

$$\begin{aligned} {b}= \begin{pmatrix}{b}_s\\ -{\varPhi }_{22}^{-1}({\varPhi }_{21}{b}_s - {\varPhi }_{2,{a},R})\end{pmatrix} = \begin{pmatrix}\mathcal {E}^{-1} {\tilde{{b}}}_s\\ -{\varPhi }_{22}^{-1}({\varPhi }_{21}\mathcal {E}^{-1} {\tilde{{b}}}_s - {\varPhi }_{2,{a},R})\end{pmatrix}. \end{aligned}$$

Now focus on (9b)

$$\begin{aligned} \begin{aligned} 0&={\dot{{a}}}+\mathcal {E}^{-1}\mathcal {D}{a}\\&\quad -\mathcal {E}^{-1}[{\varPhi }_{11}\,{\varPhi }_{12}] \begin{pmatrix}\mathcal {E}^{-1} {\tilde{{b}}}_s\\ -{\varPhi }_{22}^{-1}({\varPhi }_{21}\mathcal {E}^{-1} {\tilde{{b}}}_s - {\varPhi }_{2,{a},R})\end{pmatrix} +\mathcal {E}^{-1}{\varPhi }_{1,{a},R}. \end{aligned} \end{aligned}$$

Expand the vector product

$$\begin{aligned} 0={\dot{{a}}}+\mathcal {E}^{-1}\mathcal {D}{a}-\mathcal {E}^{-1}{\varPhi }_{11}\mathcal {E}^{-1} {\tilde{{b}}}_s+\mathcal {E}^{-1}{\varPhi }_{12}{\varPhi }_{22}^{-1}({\varPhi }_{21}\mathcal {E}^{-1} {\tilde{{b}}}_s - {\varPhi }_{2,{a},R})+\mathcal {E}^{-1}{\varPhi }_{1,{a},R}, \end{aligned}$$

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

$$\begin{aligned} \mathcal {E}{{\dot{{b}}}}_s+\mathcal {D}{b}_s+\mathcal {E}_R{b}_s+\mathcal {E}{a}g_z + G&=\mathcal {E}\dot{{b}}_s+\dot{\mathcal {E}}{b}_s-\mathcal {D}^T{b}_s+\mathcal {E}{a}g_z + G, \end{aligned}$$

(19)

using the relation

$$\begin{aligned} {\dot{\mathcal {E}}}=\mathcal {D}+\mathcal {D}^T+\mathcal {E}_R, \end{aligned}$$

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 of $\tilde{{b}}_s$ to obtain

$$\begin{aligned} \mathcal {E}{\dot{{b}}}_s+\mathcal {D}{b}_s+\mathcal {E}_R{b}_s+\mathcal {E}{a}g_z + G=\dot{\tilde{{b}}}_s-\mathcal {D}^T\mathcal {E}^{-1}\tilde{{b}}_s+\mathcal {E}{a}g_z + G. \end{aligned}$$

(20)

Next, consider the second part of (9a)

$$\begin{aligned}&\partial _{{a}}\left( \frac{1}{2}{b}^T{\varPhi }{b}-{\varPhi }_R^T{b}\right) =\partial _{{a}}\left( -{\varPhi }_R^T\begin{pmatrix}\mathcal {E}^{-1} {\tilde{{b}}}_s\\ -{\varPhi }_{22}^{-1}({\varPhi }_{21}\mathcal {E}^{-1} {\tilde{{b}}}_s - {\varPhi }_{2,{a},R})\end{pmatrix}\right. \\&\qquad \left. +\,\frac{1}{2} \begin{pmatrix}\mathcal {E}^{-1} {\tilde{{b}}}_s\\ -{\varPhi }_{22}^{-1}({\varPhi }_{21}\mathcal {E}^{-1} {\tilde{{b}}}_s - {\varPhi }_{2,{a},R})\end{pmatrix}^T{\varPhi }\begin{pmatrix}\mathcal {E}^{-1} {\tilde{{b}}}_s\\ -{\varPhi }_{22}^{-1}({\varPhi }_{21}\mathcal {E}^{-1} {\tilde{{b}}}_s - {\varPhi }_{2,{a},R})\end{pmatrix} \right) \\&\quad =\partial _{{a}}\left( -{\varPhi }_{1,{a},R}^T\mathcal {E}^{-1}\tilde{{b}}_s+{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}\left( {\varPhi }_{21}\mathcal {E}^{-1}\tilde{{b}}_s-{\varPhi }_{2,{a},R}\right) \right. \\&\qquad +\frac{1}{2}\tilde{{b}}_s^T\mathcal {E}^{-1}{\varPhi }_{11}\mathcal {E}^{-1}\tilde{{b}}_s-\left( \tilde{{b}}_s^T\mathcal {E}^{-1}{\varPhi }_{12}-{\varPhi }_{2,{a},R}^T\right) {\varPhi }_{22}^{-1}{\varPhi }_{21}\mathcal {E}^{-1}\tilde{{b}}_s\\&\qquad \left. +\,\frac{1}{2}\left( \tilde{{b}}_s^T\mathcal {E}^{-1}{\varPhi }_{12}-{\varPhi }_{2,{a},R}^T\right) {\varPhi }_{22}^{-1} \left( {\varPhi }_{21}\mathcal {E}^{-1} \tilde{{b}}_s-{\varPhi }_{2,{a},R}\right) \right) . \end{aligned}$$

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

$$\begin{aligned} \begin{aligned} \partial _{{a}}\left( \frac{1}{2}{b}^T{\varPhi }{b}-{\varPhi }_R^T{b}\right)&=\partial _{{a}}\left( -{\varPhi }_{1,{a},R}^T\mathcal {E}^{-1}\tilde{{b}}_s+{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{21}\mathcal {E}^{-1}\tilde{{b}}_s\right. \\&\quad \left. -\,\frac{1}{2}{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{2,{a},R}+ \frac{1}{2}\tilde{{b}}_s^T\mathcal {E}^{-1}S\mathcal {E}^{-1}\tilde{{b}}_s\right) , \end{aligned} \end{aligned}$$

(21)

which can be combined with (20) to obtain the statement of the Lemma. $\square $

Theorem 1

The discrete variational form corresponding to the Lagrangian (6) is equivalent with the forced Hamiltonian system

$$\begin{aligned} \begin{aligned} \dot{{a}}(t)&= \partial _{{\tilde{{b}}}_s} \mathcal {H}(t,{a}, {\tilde{{b}}}_s)\\ \dot{\tilde{{b}}}_s(t)&=-\partial _{{a}} \mathcal {H}(t,{a}, {\tilde{{b}}}_s) , \end{aligned} \end{aligned}$$

(22)

where

$$\begin{aligned} \begin{aligned} \mathcal {H}(t,{a}, {\tilde{{b}}}_s)&=\frac{{a}^T\mathcal {E}{a}}{2}g_z-\tilde{{b}}_s^T\mathcal {E}^{-1}\mathcal {D}{a}-{\varPhi }_{1,{a},R}^T\mathcal {E}^{-1}\tilde{{b}}_s+{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{21}\mathcal {E}^{-1}\tilde{{b}}_s\\&\quad -\frac{1}{2}{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{2,{a},R}+\frac{1}{2}\tilde{{b}}_s^T\mathcal {E}^{-1}{S}\mathcal {E}^{-1}\tilde{{b}}_s+G\cdot {a}. \end{aligned} \end{aligned}$$

(23)

Proof

Obviously,

$$\begin{aligned} \partial _{{a}}\mathcal {H}&=\mathcal {E}{a}g_z-\mathcal {D}^T\mathcal {E}^{-1}\tilde{{b}}_s+\partial _{{a}}\left( -{\varPhi }_{1,{a},R}^T\mathcal {E}^{-1}{\tilde{{b}}}_s+{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{21}\mathcal {E}^{-1}\tilde{{b}}_s\right. \\&\quad \left. -\,\frac{1}{2}{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{2,{a},R}+\frac{1}{2}\tilde{{b}}_s^T\mathcal {E}^{-1}{S}\mathcal {E}^{-1}\tilde{{b}}_s\right) +G. \end{aligned}$$

On the other hand

$$\begin{aligned} \begin{aligned} \partial _{\tilde{{b}}_s}\mathcal {H}&=-\mathcal {E}^{-1}\mathcal {D}{{a}}-\mathcal {E}^{-1}{\varPhi }_{1,{a},R}+\mathcal {E}^{-1}{\varPhi }_{12}{\varPhi }_{22}^{-1}{\varPhi }_{2,{a},R}+\mathcal {E}^{-1}{S}\mathcal {E}^{-1}\tilde{{b}}_s. \end{aligned} \end{aligned}$$

$\square $

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. $\frac{\mathrm {d}\mathcal H}{\mathrm {d}t}=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

$$\begin{aligned}&\partial _{a_j}\left( -{\varPhi }_{1,{a},R}^T\mathcal {E}^{-1}{\tilde{{b}}}_s+ {\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{21}\mathcal {E}^{-1}{\tilde{{b}}}_s\right. \nonumber \\&\left. \qquad -\,\frac{1}{2}{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{2,{a},R}+ \frac{1}{2}{\tilde{{b}}}_{s}^T\mathcal {E}^{-1}{S}\mathcal {E}^{-1}{\tilde{{b}}}_{s}\right) \nonumber \\&\quad =\partial _{a_j}\left( -{\varPhi }_{1,{a},R}^T{b}_s+ {\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{21}{b}_s- \frac{1}{2}{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}{\varPhi }_{2,{a},R}+ \frac{1}{2}{b}_{s}^T{S}{b}_{s}\right) \nonumber \\&\quad =-\partial _{a_j}{\varPhi }_{1,{a},R}^T{b}_s+ \partial _{a_j}({\varPhi }_{2,{a},R}^T){\varPhi }_{22}^{-1}{\varPhi }_{21}{b}_s- {\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}\partial _{a_j}({\varPhi }_{22}){\varPhi }_{22}^{-1}{\varPhi }_{21}{b}_s\nonumber \\&\qquad + {\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}\partial _{a_j}({\varPhi }_{21}){b}_s- \frac{1}{2}\partial _{a_j}({\varPhi }_{2,{a},R}^T){\varPhi }_{22}^{-1}{\varPhi }_{2,{a},R}\nonumber \\&\qquad + \frac{1}{2}{\varPhi }_{2,{a},R}^T {\varPhi }_{22}^{-1}\partial _{a_j}({\varPhi }_{22}){\varPhi }_{22}^{-1}{\varPhi }_{2,{a},R}- \frac{1}{2}{\varPhi }_{2,{a},R}^T{\varPhi }_{22}^{-1}\partial _{a_j}({\varPhi }_{2,{a},R})\nonumber \\&\qquad + \frac{1}{2}{b}_s^T\partial _{a_j}{\varPhi }_{11}{b}_s- \frac{1}{2}{b}_s^T\partial _{a_j}({\varPhi }_{12}){\varPhi }_{22}^{-1}{\varPhi }_{21}{b}_s\nonumber \\&\qquad + \frac{1}{2}{b}_s^T{\varPhi }_{12}{\varPhi }_{22}^{-1}\partial _{a_j}({\varPhi }_{22}){\varPhi }_{22}^{-1}{\varPhi }_{21}{b}_s- \frac{1}{2}{b}_s^T{\varPhi }_{12}{\varPhi }_{22}^{-1}\partial _{a_j}({\varPhi }_{21}){b}_s\nonumber \\&\quad = -\partial _{a_j}{\varPhi }_R^T{b}+\frac{1}{2}{b}^T\partial _{a_j}{\varPhi }{b}, \end{aligned}$$

(24)

where we have used the identity $\frac{\partial }{\partial t}A^{-1}=A^{-1}\frac{\partial A}{\partial t}A^{-1}$. Equation (24) greatly shortens expressions whenever $\mathcal {E}$ and ${\tilde{{b}}}_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 ${\varPhi }$.

4 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

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{{p}} = \partial _{q}H(t,{p},{q})\\ \dot{{q}} = -\partial _{p}H(t,{p},{q}), \end{array}\right. } \end{aligned}$$

we introduce the new variables ${P}=({p},\tau _0)$ and ${Q}=({q},\tau )$ and the Hamiltonian ${\tilde{H}}$ with ${\tilde{H}}({P},{Q})= H(\tau ,{p},{q})-\tau _0$. Here $\tau $ corresponds to time and the fictitious variable $\tau _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

$$\begin{aligned} {\left\{ \begin{array}{ll} {P}_{n+\frac{1}{2}} = {P}_{n} + \frac{{\varDelta }t}{2} \partial _{Q}{\tilde{H}}({P}_{n+\frac{1}{2}}, {Q}_{n})\\ {Q}_{n+1} = {Q}_{n} - \frac{{\varDelta }t}{2} \left( \partial _{P}{\tilde{H}}({P}_{n+\frac{1}{2}}, {Q}_{n}) + \partial _{P}{\tilde{H}}({P}_{n+\frac{1}{2}}, {Q}_{n+1})\right) \\ {P}_{n+1} = {P}_{n+\frac{1}{2}} + \frac{{\varDelta }t}{2} \partial _{Q}{\tilde{H}}({P}_{n+\frac{1}{2}}, {Q}_{n+1}) , \end{array}\right. } \end{aligned}$$

which in the original variables gives

$$\begin{aligned} {\left\{ \begin{array}{ll} {p}_{n+\frac{1}{2}} = {p}_{n} + \frac{{\varDelta }t}{2} \partial _{q}H(t_n,{p}_{n+\frac{1}{2}}, {q}_{n}) \\ {q}_{n+1} = {q}_{n} - \frac{{\varDelta }t}{2} \left( \partial _{p}H(t_n,{p}_{n+\frac{1}{2}}, {q}_{n}) + \partial _{p}H(t_{n+1},{p}_{n+\frac{1}{2}}, {q}_{n+1})\right) \\ \tau _{n+1} = \tau _{n}+{\varDelta }t\\ {p}_{n+1} = {p}_{n+\frac{1}{2}} + \frac{{\varDelta }t}{2} \partial _{q}H(t_{n+1},{p}_{n+\frac{1}{2}}, {q}_{n+1}). \end{array}\right. } \end{aligned}$$

(25)

The update for $\tau _n$ ensures that $\tau $ indeed represents time and will be taken for granted in the following. The fictitious variable $\tau _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={\tilde{b}}_s$ in (25) is now equal to

$$\begin{aligned} {a}_{n+\frac{1}{2}}&={a}_n+\frac{{\varDelta }t}{2}\left( -\mathcal {E}_n^{-1}\mathcal {D}_n{a}_{n+\frac{1}{2}}-\mathcal {E}_n^{-1}{\varPhi }_{1,{a}_{n+\frac{1}{2}},R_n}\right. \\&\quad \left. +\,\mathcal {E}_n^{-1}{\varPhi }_{12,n,{a}_{n+\frac{1}{2}}}{\varPhi }_{22,n,{a}_{n+\frac{1}{2}}}^{-1}{\varPhi }_{2,a_{n+\frac{1}{2}},R_n}+ \mathcal {E}_n^{-1}{S}_{n,{a}_{n+\frac{1}{2}}}\mathcal {E}_n^{-1}\tilde{{b}}_{s,n}\right) \\ \tilde{{b}}_{s,n+1}&=\tilde{{b}}_{s,n}-\frac{{\varDelta }t}{2}\left( g_z(\mathcal {E}_n+\mathcal {E}_{n+1}){a}_{n+\frac{1}{2}}-\mathcal {D}_n^T\mathcal {E}_n^{-1}\tilde{{b}}_{s,n}-\mathcal {D}_{n+1}^T\mathcal {E}_{n+1}^{-1}\tilde{{b}}_{s,n+1}\right. \\&\quad +\frac{1}{2}{b}_{n,a_{n+\frac{1}{2}}}^T\partial _{a}{\varPhi }_{n,{a}_{n+\frac{1}{2}}}{b}_{n,a_{n+\frac{1}{2}}}+ \frac{1}{2}{b}_{n+1,a_{n+\frac{1}{2}}}^T\partial _{a}{\varPhi }_{n+1,{a}_{n+\frac{1}{2}}}{b}_{n+1,a_{n+\frac{1}{2}}}\\&\quad \left. -\,\partial _{a}{\varPhi }_{a_{n+\frac{1}{2}},R_n}\cdot {b}_{n,a_{n+\frac{1}{2}}}- \partial _{a}{\varPhi }_{a_{n+\frac{1}{2}},R_{n+1}}\cdot {b}_{n+1,a_{n+\frac{1}{2}}} +G^n+G^{n+1}\right) \\ {a}_{n+1}&={a}_{n+\frac{1}{2}}+\frac{{\varDelta }t}{2}\left( -\mathcal {E}_{n+1}^{-1}\mathcal {D}_{n+1}{a}_{n+\frac{1}{2}}-\mathcal {E}_{n+1}^{-1}{\varPhi }_{1,{a}_{n+\frac{1}{2}},R_{n+1}}\right. \\&\quad +\mathcal {E}_{n+1}^{-1}{\varPhi }_{12,n+1,{a}_{n+\frac{1}{2}}}{\varPhi }_{22,n+1,{a}_{n+\frac{1}{2}}}^{-1}{\varPhi }_{2,{a}_{n+\frac{1}{2}},R_{n+1}}\\&\quad \left. +\,\mathcal {E}_{n+1}^{-1}{S}_{n+1,{a}_{n+\frac{1}{2}}}\mathcal {E}_{n+1}^{-1}\tilde{{b}}_{s,n+1}\right) , \end{aligned}$$

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

$$\begin{aligned} \mathcal {E}_{n}{a}_{n+\frac{1}{2}}&=\mathcal {E}_n{a}_n+\frac{{\varDelta }t}{2}\left( -\mathcal {D}_n{a}_{n+\frac{1}{2}}-{\varPhi }_{1,{a}_{n+\frac{1}{2}},R_n}\right. \nonumber \\&\quad \left. +\,{\varPhi }_{12,n,{a}_{n+\frac{1}{2}}}{\varPhi }_{22,n,{a}_{n+\frac{1}{2}}}^{-1}{\varPhi }_{2,{a}_{n+\frac{1}{2}},R_n}+{S}_{n,{a}_{n+\frac{1}{2}}}{b}_{s,n}\right) \nonumber \\ \mathcal {E}_{n+1}{b}_{s,n+1}&=\mathcal {E}_n{b}_{s,n}-\frac{{\varDelta }t}{2}\left( g_z(\mathcal {E}_n+\mathcal {E}_{n+1}){a}_{n+\frac{1}{2}} -\mathcal {D}_n^T{b}_{s,n}-\mathcal {D}_{n+1}^T{b}_{s,n+1}\right. \nonumber \\&\quad + \frac{1}{2}{b}_{n,a_{n+\frac{1}{2}}}^T\partial _{a}{\varPhi }_{n,{a}_{n+\frac{1}{2}}}{b}_{n,a_{n+\frac{1}{2}}}+ \frac{1}{2}{b}_{n+1,a_{n+\frac{1}{2}}}^T\partial _{a}{\varPhi }_{n+1,{a}_{n+\frac{1}{2}}}{b}_{n+1,a_{n+\frac{1}{2}}}\nonumber \\&\quad \left. -\,\partial _{a}{\varPhi }_{a_{n+\frac{1}{2}},R_n}^T{b}_{n,a_{n+\frac{1}{2}}}- \partial _{a}{\varPhi }_{a_{n+\frac{1}{2}},R_{n+1}}^T{b}_{n+1,a_{n+\frac{1}{2}}} +G^n+G^{n+1}\right) \nonumber \\ \mathcal {E}_{n+1}{a}_{n+1}&=\mathcal {E}_{n+1}{a}_{n+\frac{1}{2}}+\frac{{\varDelta }t}{2}\left( -\mathcal {D}_{n+1}{a}_{n+\frac{1}{2}}-{\varPhi }_{1,{a}_{n+\frac{1}{2}},R_{n+1}}\right. \nonumber \\&\quad \left. +\,{\varPhi }_{12,n+1,{a}_{n+\frac{1}{2}}}{\varPhi }_{22,n+1,{a}_{n+\frac{1}{2}}}^{-1}{\varPhi }_{2,{a}_{n+\frac{1}{2}},R_{n+1}}+{S}_{n+1,{a}_{n+\frac{1}{2}}}{b}_{s,n+1}\right) ,\nonumber \\ \end{aligned}$$

(26)

Since the time stepping in (26) is implicit, we first solve the equation for ${a}_{n+\frac{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.
Evaluate the matrices $\mathcal {E}_t, \mathcal {D}_t, {\varPhi }_{t,a_h}, {\varPhi }_{R,t,a_h}$ and $G_t$ on the current mesh.
while $t<t_{{ end}}$
- Iterate the Newton algorithm (33) until it converges, while moving the mesh using the new free surface height $\eta _h$ in (28) and updating ${\varPhi }_{t,a_h}$ and ${\varPhi }_{R,t,a_h}$ to account for the new free surface position.
- Increase $t=t+dt$, update the mesh to satisfy the new position of the wave maker and reevaluate the matrices $\mathcal {E}_t, \mathcal {D}_t, {\varPhi }_{t,a_h}, {\varPhi }_{R,t,a_h}$ and $G_t$.
- Iterate the Newton algorithm (35) until it converges.
- Solve the third equation from (26).
- Move the mesh and update ${\varPhi }_{t,a_h}$ and ${\varPhi }_{R,t,a_h}$ to account for the new free surface position.

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.

4.1 Computing Derivatives $\partial _{a}{\varPhi }$ and $\partial _{a}{\varPhi }_R$

In the derivation of the discrete Hamiltonian the derivatives $\partial _a{\varPhi }$ and $\partial _a{\varPhi }_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 $\partial _a{\varPhi }$ element-wise, as a summation on the finite element tessellation ${\mathcal {T}}_h$:

$$\begin{aligned} \frac{\partial }{\partial a_k}\sum _{K\in \mathcal T_h}\int _K\nabla \phi _i\cdot \nabla \phi _j\,\mathrm {d}K, \end{aligned}$$

where the shape of the element and the basis functions depend on the free surface $\eta _h$, hence implicitly on the coefficients a. Introduce a reference element ${\hat{K}}$. We will denote the image of the basis functions on ${\hat{K}}$ as ${\hat{\phi }}$, the reference coordinates as $({\hat{x}}, {\hat{y}}, {\hat{z}})$ and the gradient operator with respect to reference coordinates as ${\hat{\nabla }}$. Assume, for every element $K\in \mathcal T_h$, that there is an invertible mapping $F_K\,{:}\,{\hat{K}}\rightarrow K$. Since we use nodal basis functions, we have $x=\sum _l {\hat{x}}_l\hat{\phi _l}$, where ${\hat{x}}_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=\sum _l {\hat{x}}_l{\hat{\nabla }}\hat{\phi _l}^T$.

Perform the coordinate transformation,

$$\begin{aligned}&\frac{\partial }{\partial a_k}\int _K\nabla \phi _i\cdot \nabla \phi _j\,\mathrm {d}K = \int _{\hat{K}}\frac{\partial }{\partial a_k}\left( J^{-T}\hat{\nabla }\hat{\phi }_i\right) ^T\left( J^{-T}\hat{\nabla }\hat{\phi }_j\right) \left| J\right| \,\mathrm {d}\hat{K}. \end{aligned}$$

Using the matrix identities

$$\begin{aligned} \frac{\partial }{\partial t}A^{-1}=-A^{-1}\frac{\partial A}{\partial t}A^{-1} \end{aligned}$$

and

$$\begin{aligned} \frac{\partial }{\partial t}|A| = \mathrm {Tr}\left( A^{-1}\frac{\partial A}{\partial t}\right) |A|, \end{aligned}$$

with $|A|=\det (A)$, we obtain

$$\begin{aligned}&\frac{\partial }{\partial a_k}\int _K\nabla \phi _i\cdot \nabla \phi _j\,\mathrm {d}K = \int _{\hat{K}}\left[ -\left( J^{-T}\hat{\nabla }\hat{\phi }_l\left( \frac{\partial \hat{x}_l}{\partial a_k}\right) ^TJ^{-T}\hat{\nabla }\hat{\phi }_i\right) ^TJ^{-T}\hat{\nabla }\hat{\phi }_j\right. \\&\quad -\left( J^{-T}\hat{\nabla }\hat{\phi }_i\right) ^TJ^{-T}\hat{\nabla }\hat{\phi }_l\left( \frac{\partial \hat{x}_l}{\partial a_k}\right) ^TJ^{-T}\hat{\nabla }\hat{\phi }_j\\&\quad \left. +\,\left( J^{-T}\hat{\nabla }\hat{\phi }_i\right) ^TJ^{-T}\hat{\nabla }\hat{\phi }_j\mathrm {Tr}{\left( \frac{\partial \hat{x}_l}{\partial a_k}\hat{\nabla }\hat{\phi }_l^TJ^{-1}\right) }\right] \left| J\right| \,\mathrm {d}\hat{K}, \end{aligned}$$

where the summation convention is used on repeated indices. Transforming back to the elements $K\in {\mathcal {T}}_h$ we obtain the relation

$$\begin{aligned} \frac{\partial }{\partial a_k}\int _K\nabla \phi _i\cdot \nabla \phi _j\,\mathrm {d}K&= \int _{K}\left( -\left( \nabla \phi _i\cdot \frac{\partial {\hat{x}}_l}{\partial a_k}\right) \left( \nabla \phi _l\cdot \nabla \phi _j\right) \right. \\&\quad -\left( \nabla \phi _j\cdot \frac{\partial {\hat{x}}_l}{\partial a_k}\right) \left( \nabla \phi _l\cdot \nabla \phi _i\right) \\&\quad \left. +\,\left( \nabla \phi _l\cdot \frac{\partial {\hat{x}}_l}{\partial a_k}\right) \left( \nabla \phi _i\cdot \nabla \phi _j\right) \right) \,\mathrm {d}K. \end{aligned}$$

The coupling between the node locations and the free surface parametrization $\frac{\partial {\hat{x}}_l}{\partial a_k}$ has to be constructed depending on the choice of the mesh movement algorithm.

For the computation of $\partial _{a}{\varPhi }_R$ we can use (13). We have

$$\begin{aligned} \partial _{a_k}{\varPhi }_R&=\frac{\partial }{\partial {{a}_k}}\int _{\partial {\varOmega }_{t,R,h}} \phi _i \frac{\partial R}{\partial t}\cdot \nu _R\,\mathrm {d}S\\&=\int _{\partial {\varOmega }_{t, R, h}}\frac{\partial \phi _i}{\partial a_k}\frac{\partial R}{\partial t}\cdot \nu _R\,\mathrm {d}S\\&\quad +\int _{\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}}\left( \frac{\partial R}{\partial t}\cdot \nu _R\right) (\tau \times \nu _R)\cdot e_z\eta _k\phi _i\,\mathrm {d}l. \end{aligned}$$

For the first term use the chain rule and the mapping $F_K$ and for the second term the compatibility condition (5) to find

$$\begin{aligned} \partial _{a_k}{\varPhi }_R=\int _{\partial {\varOmega }_{t,R,h}}\left( \frac{\partial R}{\partial t}\cdot \nu _R\right) \phi _l\frac{\partial {\hat{x}}_l}{\partial a_k}\cdot \nabla \phi _i \,\mathrm {d}S+\mathcal {E}_R. \end{aligned}$$

We would like to consider the mesh deformation and the rest of the time stepping scheme separately. To this end, we introduce the matrices

$$\begin{aligned}{}[{C}]_{j,l}&=b_i\sum _{K\in {\mathcal {T}}_h}\int _K\left( -\nabla \phi _i(\nabla \phi _l\cdot \nabla \phi _j)-\nabla \phi _j(\nabla \phi _l\cdot \nabla \phi _i)\right. \\&\quad \left. +\,\nabla \phi _l(\nabla \phi _i\cdot \nabla \phi _j)\right) \,\mathrm {d}K \end{aligned}$$

and

$$\begin{aligned}{}[{B}]_{i,l}=\int _{\partial {\varOmega }_{t,R,h}}\left( \frac{\partial R}{\partial t}\cdot \nu _R\right) \phi _l\nabla \phi _i\,\mathrm {d}S \end{aligned}$$

and we obtain the relations

$$\begin{aligned} b_i\frac{\partial }{\partial a_k}\sum _{K\in \mathcal T_h}\int _K\nabla \phi _i\cdot \nabla \phi _j\,\mathrm {d}K={C}_{j,l}\frac{\partial \hat{x}_l}{\partial a_k} \end{aligned}$$

and

$$\begin{aligned} \partial _{a_k}{\varPhi }_{R,j}=\mathcal {E}_R + {B}_{j,l}\frac{\partial \hat{x}_l}{\partial a_k}. \end{aligned}$$

The matrices ${C}$ and ${B}$ are separated into free surface and interior parts

$$\begin{aligned} {C}=\begin{pmatrix} {C}_{11}&{}{C}_{12}\\ {C}_{21}&{}{C}_{22} \end{pmatrix}\quad {B}=\begin{pmatrix} {B}_{11}&{}{B}_{12}\\ {B}_{21}&{}{B}_{22} \end{pmatrix}, \end{aligned}$$

similar to ${\varPhi }$. The node velocity $\frac{\partial {\hat{x}}_i}{\partial 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

$$\begin{aligned} \begin{pmatrix} {I}&{}0\\ {F}_{21}&{}{F}_{22} \end{pmatrix} \begin{pmatrix} {u}_s\\ {u}_i \end{pmatrix} = \begin{pmatrix} {a}\\ 0 \end{pmatrix}, \end{aligned}$$

where ${F}$ is the Jacobian with respect to the node displacements of the mesh movement algorithm. Inverting the Jacobian gives

$$\begin{aligned} \begin{pmatrix} {u}_s\\ {u}_i \end{pmatrix} = \begin{pmatrix} {I}&{}\,0\\ -{F}_{22}^{-1}{F}_{21}&{}\,{F}_{22}^{-1} \end{pmatrix} \begin{pmatrix} {a}\\ 0 \end{pmatrix}. \end{aligned}$$

The node displacements and the node position are linked by a constant offset, hence we directly obtain

$$\begin{aligned} \frac{\partial \hat{x}_i}{\partial a_j}=\begin{pmatrix} {I}\\ -{F}_{22}^{-1}{F}_{21} \end{pmatrix}_{i,j} \end{aligned}$$

(27)

4.2 Mesh Deformation Algorithm

We base the mesh deformation algorithm on Masud and Hughes [20]. The idea is to compute a displacement field ${u}\in {\mathbb {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

$$\begin{aligned} {\left\{ \begin{array}{ll} \nabla \cdot \left( \left( 1+\tau \right) \nabla {u}\right) = 0&{} \quad \text {on}\;{\varOmega }_z\\ n\cdot {u} = \eta &{}\quad \text {on}\;\partial {\varOmega }_{z,s,h} \\ n\cdot {u} = R&{}\quad \text {on}\;\partial {\varOmega }_{z,R,h}, \end{array}\right. } \end{aligned}$$

(28)

where $\tau $ is a bounded nonnegative function and the zero displacement domain ${\varOmega }_z$ is also used to compute the displacements. The free surface height $\eta $ 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 $\eta $. While this is given as a continuous system of equations, they are discretized using linear basis functions on ${\mathcal {T}}_h$ in order to guarantee that the computed nodal displacements u can directly be used to deform the mesh. The parameter $\tau $ 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 $\frac{\partial {\hat{x}}_l}{\partial a_k}$ we take the derivative $\partial _{a}$ of (28). Since (28) is solved on a fixed domain, there are no hidden dependencies and we can directly write

$$\begin{aligned} {\left\{ \begin{array}{ll} \nabla \cdot \left( \left( 1+\tau \right) \nabla \left( \frac{\partial \hat{x}}{\partial a_k}\right) \right) = 0 &{}\quad \text {on}\;{\varOmega }_z\\ n\cdot \frac{\partial \hat{x}}{\partial a_k} = \phi _k &{}\quad \text {on}\;\partial {\varOmega }_{z,s,h} \\ n\cdot \frac{\partial \hat{x}}{\partial a_k} = 0&{}\quad \text {on}\;\partial {\varOmega }_{z,R,h}, \end{array}\right. } \end{aligned}$$

with the understanding that these equations have to be discretized in the same way as the equations for the displacement. The derivative $\frac{\partial {\hat{x}}}{\partial t}$ can be approximated in a similar manner.

In our simulations, the small elements reside mostly near the free surface, so we choose $\tau =\mathrm {e}^{1+cz}$, where $c\cong 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 $\tau $ in element K can be computed as

$$\begin{aligned} \tau _K=\frac{1-{\varDelta }_{{ min}}/{\varDelta }_{{ max}}}{{\varDelta }_K/{\varDelta }_{{ max}}}, \end{aligned}$$

with ${\varDelta }_{{ min}}, {\varDelta }_{{ max}}$ the area (or volume) of the smallest and largest elements in the mesh and ${\varDelta }_K$ the area (or volume) of element K. In [20] it is shown that this results in $\tau _K$-values that are essentially independent of the ratio ${\varDelta }_{{ min}}/{\varDelta }_{{ max}}$. A more detailed way to compute the $\tau $-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].

5 Results

5.1 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 solution for a water depth $H=1$, gravity coefficient $g=1$ and domain length $X\cong 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 $\frac{12N_x}{32}$ time steps per period, comparing the free surface height $\eta _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.

Table 1 L2-norm of the error in the free surface height and the difference between the maximum and minimum energy measured for various numbers of elements

Full size table

5.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

$$\begin{aligned} \eta _0&=0.215\,\mathrm {sech}{(1.18x)},\\ \phi _0&=0. \end{aligned}$$

After moving away from the boundary, this initial wave profile will deform into a traveling wave closely resembling

$$\begin{aligned} \eta (x,t)=0.1\,\mathrm {sech}^2\left( \frac{x+c-\sqrt{0.6g}t}{\sqrt{2}}\right) , \end{aligned}$$

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\,\text {s}$, with ${\varDelta }t=0.05$. The domain has a reflecting wall at $x=150\,\text {m}$. In order to verify numerically the stability of our new scheme, we choose a sequence of mesh sizes ${\varDelta }x\in \{2\,\text {m}, 1\,\text {m}, 0.5\,\text {m}, 0.25\,\text {m}\}$. In the vertical direction, we reproduce the choices made in [30]. That is, we use 6 elements in the vertical direction, placing the mesh lines at $z=0.5\left( \frac{\mathrm {cosh}(-0.1\pi \{0:1/6:1\})-\mathrm {cosh}(-0.1\pi )}{1-\mathrm {cosh}(-0.1\pi )}-1\right) $. The coarsest of these meshes is unable to sufficiently resolve the traveling wave profile, while the choice ${\varDelta }x=1\,\text {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 ${\varDelta }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.

5.3 Comparison with Experiments

Finally, we made a comparison with experiments. For this purpose, we used the data set from Run 202002, which was provided by the Maritime Research Institute Netherlands (MARIN). In this experiment a piston wave maker generates a wave train of successively faster moving waves that focus into a splash near $x=50\,\text {m}$ in a model wave basin with dimension $195.4\,\text {m} \times 1\,\text {m}$. The wave maker motion in the computations is identical to the wave maker motion used in the experiments. At time $t=0$ there are no waves present in the basin. Since we are only interested in the first part of the domain, we use a numerical basin of $90\,\text {m} \times 1\,\text {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\,\text {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 ${\varDelta }x=0.0027$ near the splash zone and ${\varDelta }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.

References

Bar-Yoseph, P.Z., Mereu, S., Chippada, S., Kalro, V.J.: Automatic monitoring of element shape quality in 2-d and 3-d computational mesh dynamics. Comput. Mech. 27(5), 378–395 (2001)
Article MATH Google Scholar
Bazilevs, Y., Calo, V.M., Hughes, T.J., Zhang, Y.: Isogeometric fluid–structure interaction: theory, algorithms, and computations. Comput. Mech. 43(1), 3–37 (2008)
Article MathSciNet MATH Google Scholar
Beale, J.: A convergent boundary integral method for three-dimensional water waves. Math. Comput. 70(235), 977–1029 (2001)
Article MathSciNet MATH Google Scholar
Broeze, J., van Daalen, E.F., Zandbergen, P.J.: A three-dimensional panel method for nonlinear free surface waves on vector computers. Comput. Mech. 13(1–2), 12–28 (1993)
Article MATH Google Scholar
Craig, W., Guyenne, P., Sulem, C.: A Hamiltonian approach to nonlinear modulation of surface water waves. Wave Motion 47(8), 552–563 (2010)
Article MathSciNet MATH Google Scholar
Fenton, J., Rienecker, M.: A Fourier method for solving nonlinear water-wave problems: application to solitary–wave interactions. J. Fluid Mech. 118, 411–443 (1982)
Article MathSciNet MATH Google Scholar
Flanders, H.: Differentiation under the integral sign. Am. Math. Mon. 80(6), 615–627 (1973)
Article MathSciNet MATH Google Scholar
Fochesato, C., Dias, F.: A fast method for nonlinear three-dimensional free-surface waves. Proc. R. Soc. A Math. Phys. Eng. Sci. 462(2073), 2715–2735 (2006)
Article MathSciNet MATH Google Scholar
Fochesato, C., Grilli, S., Dias, F.: Numerical modeling of extreme rogue waves generated by directional energy focusing. Wave Motion 44(5), 395–416 (2007)
Article MathSciNet MATH Google Scholar
Gagarina, E., Ambati, V.R., van der Vegt, J.J.W., Bokhove, O.: Variational space–time (dis)continuous Galerkin method for nonlinear free surface water waves. J. Comput. Phys. 275, 459–483 (2014)
Article MathSciNet MATH Google Scholar
Guerber, E., Benoit, M., Grilli, S.T., Buvat, C.: A fully nonlinear implicit model for wave interactions with submerged structures in forced or free motion. Eng. Anal. Bound. Elem. 36(7), 1151–1163 (2012)
Article MathSciNet MATH Google Scholar
Hairer, E., Lubich, C., Wanner, G.: Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, vol. 31. Springer, Berlin (2006)
MATH Google Scholar
Hou, T.Y., Zhang, P.: Convergence of a boundary integral method for 3-d water waves. Discret. Contin. Dyn. Syst. Ser. B 2(1), 1–34 (2002)
Article MathSciNet MATH Google Scholar
Leimkuhler, B., Reich, S.: Simulating Hamiltonian Dynamics. Cambridge Monographs on Applied and Computational Mathematics, 1st edn. Cambridge University Press, Cambridge (2004)
MATH Google Scholar
Longuet-Higgins, M.S., Cokelet, E.: The deformation of steep surface waves on water. I. A numerical method of computation. Proc. R. Soc. Lond. A Math. Phys. Sci. 350(1660), 1–26 (1976)
Article MathSciNet MATH Google Scholar
Luke, J.: A variational principle for a fluid with a free surface. J. Fluid Mech. 27(02), 395–397 (1967)
Article MathSciNet MATH Google Scholar
Ma, Q., Wu, G., Eatock Taylor, R.: Finite element simulation of fully non-linear interaction between vertical cylinders and steep waves. Part 1: methodology and numerical procedure. Int. J. Numer. Methods Fluids 36(3), 265–285 (2001)
Article MATH Google Scholar
Ma, Q., Wu, G., Eatock, R.: Finite element simulations of fully non-linear interaction between vertical cylinders and steep waves. Part 2: numerical results and validation. Int. J. Numer. Methods Fluids 36(3), 287–308 (2001)
Article Google Scholar
Ma, Q., Yan, S.: Quasi ALE finite element method for nonlinear water waves. J. Comput. Phys. 212(1), 52–72 (2006)
Article MathSciNet MATH Google Scholar
Masud, A., Hughes, T.J.R.: A space–time Galerkin/least-squares finite element formulation of the Navier–Stokes equations for moving domain problems. Comput. Methods Appl. Mech. Eng. 146, 91–126 (1997)
Article MathSciNet MATH Google Scholar
Robertson, I., Sherwin, S.: Free-surface flow simulation using hp/spectral elements. J. Comput. Phys. 26(53), 26–53 (1999)
Article MathSciNet MATH Google Scholar
Stein, K., Tezduyar, T.E., Benney, R.: Automatic mesh update with the solid-extension mesh moving technique. Comput. Methods Appl. Mech. Eng. 193(21), 2019–2032 (2004)
Article MATH Google Scholar
Tezduyar, T., Benney, R.: Mesh moving techniques for fluid–structure interactions with large displacements. J. Appl. Mech. 70(1), 58–63 (2003)
Article MATH Google Scholar
Tomar, S., van der Vegt, J.J.W.: A Runge–Kutta discontinuous Galerkin method for linear free-surface gravity waves using high order velocity recovery. Comput. Methods Appl. Mech. Eng. 196(13), 1984–1996 (2007)
Article MathSciNet MATH Google Scholar
Toth, F., Kaltenbacher, M.: Fully coupled linear modelling of incompressible free-surface flow, compressible air and flexible structures. Int. J. Numer. Methods Eng. 107(11), 947–969 (2015)
Tsai, W.-T., Yue, D.K.: Computation of nonlinear free-surface flows. Annu. Rev. Fluid Mech. 28(1), 249–278 (1996)
Article MathSciNet Google Scholar
van der Vegt, J.J.W., Tomar, S.: Discontinuous Galerkin method for linear free-surface gravity waves. J. Sci. Comput. 22(1–3), 531–567 (2005)
Article MathSciNet MATH Google Scholar
Wang, C., Wu, G.: Interactions between fully nonlinear water waves and cylinder arrays in a wave tank. Ocean Eng. 37(4), 400–417 (2010)
Article MathSciNet Google Scholar
Wang, C., Wu, G.: A brief summary of finite element method applications to nonlinear wave–structure interactions. J. Mar. Sci. Appl. 10(2), 127–138 (2011)
Article Google Scholar
Westhuis J.: The numerical simulation of nonlinear waves in a hydrodynamic model test basin. Ph.D. thesis. University of Twente (2001)
Whitham, G.B.: Linear and Nonlinear Waves, vol. 42. Wiley, London (2011)
MATH Google Scholar
Wu, G., Hu, Z.: Simulation of nonlinear interactions between waves and floating bodies through a finite-element-based numerical tank. Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 460(2050), 2797–2817 (2004)
Article MathSciNet MATH Google Scholar
Zakharov, V.E.: Stability of periodic waves of finite amplitude on the surface of a deep fluid. J. Appl. Mech. Tech. Phys. 9(2), 190–194 (1968)
Article Google Scholar

Download references

Acknowledgements

We acknowledge the financial support of the Technology Foundation STW in the project “FastFEM: Behavior of Fast Ships in Waves” and we thank the Maritime Research Institute Netherlands (MARIN) for providing the experimental data. F. Izsák was supported by the János Bolyai Research Fellowship of the Hungarian Academy of Sciences.

Author information

Authors and Affiliations

Mathematics of Computational Science, University of Twente, Enschede, The Netherlands
Freekjan Brink & J. J. W. van der Vegt
Department of Applied Analysis and Computational Mathematics, Eötvös University, Pázmány Sétány 1/C, Budapest, 1117, Hungary
Ferenc Izsák

Authors

Freekjan Brink
View author publications
You can also search for this author in PubMed Google Scholar
Ferenc Izsák
View author publications
You can also search for this author in PubMed Google Scholar
J. J. W. van der Vegt
View author publications
You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Freekjan Brink.

Appendices

Appendix 1: Solution of Algebraic Equations

Solving the algebraic equations for the Hamiltonian finite element discretization requires some special care. The first two equations in (26) are nonlinear and are solved with a Newton method. The equations for the Newton updates ${\varDelta }a_{n+\frac{1}{2}}$ and ${\varDelta }b_{s,n+1}$ are, respectively,

$$\begin{aligned} \begin{aligned}&\left( \mathcal {E}_{n}- \frac{{\varDelta }t}{2} \partial _{{a}} \left[ {S}_{n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{s,n} - {\varPhi }_{1,{a}_{n+\frac{1}{2}},R_{n}}^{(k)}\right. \right. \\&\qquad \left. \left. +\,{\varPhi }_{12,n,{a}_{n+\frac{1}{2}}}^{(k)}\left( {\varPhi }_{22,n,{a}_{n+\frac{1}{2}}}^{(k)}\right) ^{-1}{\varPhi }_{2,{a}_{n+\frac{1}{2}},R_{n}}^{(k)}\right] +\frac{{\varDelta }t}{2}\mathcal {D}_n \right) {\varDelta }a_{n+\frac{1}{2}}^{(k)}\\&\quad =- \left( \mathcal {E}_{n}{a}_{n+\frac{1}{2}}^{(k)}-\mathcal {E}_n{a}_n-\frac{{\varDelta }t}{2}\left( -\mathcal {D}_n{a}_{n+\frac{1}{2}}^{(k)}-{\varPhi }_{1,{a}_{n+\frac{1}{2}},R_n}^{(k)}\right. \right. \\&\qquad \left. \left. +\,{\varPhi }_{12,n,{a}_{n+\frac{1}{2}}}^{(k)}\left( {\varPhi }_{22,n,{a}_{n+\frac{1}{2}}}^{(k)}\right) ^{-1}{\varPhi }_{2,{a}_{n+\frac{1}{2}},R_n}^{(k)}+{S}_{n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{s,n}\right) \right) , \end{aligned} \end{aligned}$$

(29)

with ${\varDelta }a_{n+\frac{1}{2}}^{(k)}=a_{n+\frac{1}{2}}^{(k+1)}-a_{n+\frac{1}{2}}^{(k)}$ and

$$\begin{aligned}&\left( \mathcal {E}_{n+1}+\frac{{\varDelta }t}{2}\left[ \partial _{{a}}{\varPhi }_{11,n+1,{a}_{n+\frac{1}{2}}}{b}_{s,n+1}^{(k)} +\partial _{{a}}{\varPhi }_{12,n+1,{a}_{n+\frac{1}{2}}}{b}_{i,n+1,a_{n+\frac{1}{2}}} \right. \right. \nonumber \\&\qquad \left. \left. -\,\partial _{{a}}{\varPhi }_{1,{a}_{n+\frac{1}{2}},R_{n+1}} -\mathcal {D}_{n+1}\right] ^T\right) {\varDelta }{b}_{s,n+1}^{(k)}\nonumber \\&\quad = -\left\{ \mathcal {E}_{n+1}{b}_{s,n+1}^{(k)}-\mathcal {E}_n{b}_{s,n}+\frac{{\varDelta }t}{2}\left( g_z(\mathcal {E}_{n}+\mathcal {E}_{n+1}){a}_{n+\frac{1}{2}}-\mathcal {D}_n^T{b}_{s,n}\right. \right. \nonumber \\&\qquad -\mathcal {D}_{n+1}^T{b}_{s,n+1}^{(k)}+ \frac{1}{2}{b}_{n, a_{n+\frac{1}{2}}}^T\partial _{a}{\varPhi }_{n,{a}_{n+\frac{1}{2}}}{b}_{n, a_{n+\frac{1}{2}}}+ \frac{1}{2}{b}_{n+1, a_{n+\frac{1}{2}}}^{(k)\;T}\partial _{a}{\varPhi }_{n+1,{a}_{n+\frac{1}{2}}}{b}_{n+1, a_{n+\frac{1}{2}}}^{(k)}\nonumber \\&\qquad \left. \left. -\,\partial _{a}{\varPhi }_{{a}_{n+\frac{1}{2}},R_{n}}\cdot {b}_{n, a_{n+\frac{1}{2}}}- \partial _{a}{\varPhi }_{{a}_{n+\frac{1}{2}},R_{n+1}}\cdot {b}_{n+1, a_{n+\frac{1}{2}}}^{(k)}+G^n+G^{n+1}\right) \right\} ,\nonumber \\ \end{aligned}$$

(30)

with ${\varDelta }b_{s,n+1}^{(k)}=b_{s,n+1}^{(k+1)}-b_{s,n+1}^{(k)}$ and where ${a}_{n+\frac{1}{2}}^{(0)}={a}_n$ and ${b}_{s,n+1}^{(0)}={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 $\partial _{{a}_i}$ and $\partial _{{b}_{s,j}}$ commute. Recall that ${S}={\varPhi }_{11}-{\varPhi }_{12}{\varPhi }_{22}^{-1}{\varPhi }_{21}$.

For numerical efficiency reasons it is crucial to avoid explicitly forming ${\varPhi }_{22}^{-1}$. The following auxiliary equation is therefore introduced

$$\begin{aligned} {\varPhi }_{22,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{i,n,{a}_{n+\frac{1}{2}}}^{(k)} =-\left( {\varPhi }_{21,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{s,n}-{\varPhi }_{2,{a}_{n+\frac{1}{2}}, R_n}^{(k)}\right) . \end{aligned}$$

Substituting this into (29), we find

$$\begin{aligned}&\left( \begin{array}{cc} {A}_{11}&{}{A}_{12}\\ {A}_{21}&{}{A}_{22} \end{array}\right) \begin{pmatrix} {\varDelta }{a}_{n+\frac{1}{2}}^{(k)}\\ {\varDelta }{b}_{i,n,{a}_{n+\frac{1}{2}}}^{(k)} \end{pmatrix}=\begin{pmatrix}{v}_1\\ {v}_2\end{pmatrix},\\ {A}_{11}&=\mathcal {E}_n-\frac{{\varDelta }t}{2} \partial _{{a}} \left[ {\varPhi }_{11,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{s,n}+{\varPhi }_{12,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{i,n,{a}_{n+\frac{1}{2}}}^{(k)} - {\varPhi }_{1,{a}_{n+\frac{1}{2},k},R_{n}}\right] \nonumber \\&\quad +\frac{{\varDelta }t}{2}\mathcal {D}_n,\nonumber \\ {A}_{12}&=-\frac{{\varDelta }t}{2}{\varPhi }_{12,n,{a}_{n+\frac{1}{2}}}^{(k)},\nonumber \\ {A}_{21}&=-\frac{{\varDelta }t}{2} \partial _{{a}} \left[ {\varPhi }_{21,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{s,n}+{\varPhi }_{22,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{i,n,{a}_{n+\frac{1}{2}}}^{(k)} - {\varPhi }_{2,{a}_{n+\frac{1}{2},k},R_{n}}\right] ,\nonumber \\ {A}_{22}&=-\frac{{\varDelta }t}{2}{\varPhi }_{22,n,{a}_{n+\frac{1}{2}}}^{(k)},\nonumber \\ {v}_1&=-\left( \mathcal {E}_{n}{a}_{n+\frac{1}{2}}^{(k)}-\mathcal {E}_n{a}_n-\frac{{\varDelta }t}{2}\left( -\mathcal {D}_n{a}_{n+\frac{1}{2}}^{(k)}-{\varPhi }_{1,{a}_{n+\frac{1}{2}},R_n}^{(k)}\right. \right. \nonumber \\&\quad \left. \left. +\,{\varPhi }_{11,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{s,n}+{\varPhi }_{12,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{i,n,{a}_{n+\frac{1}{2}}}^{(k)}\right) \right) ,\nonumber \\ {v}_2&=-\frac{{\varDelta }t}{2}\left( {\varPhi }_{2,{a}_{n+\frac{1}{2}},R_n}^{(k)}-{\varPhi }_{21,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{s,n}-{\varPhi }_{22,n,{a}_{n+\frac{1}{2}}}^{(k)}{b}_{i,n,{a}_{n+\frac{1}{2}}}^{(k)}\right) ,\nonumber \end{aligned}$$

(31)

where ${\varDelta }b_{i,n+1,a_{n+\frac{1}{2}}}^{(k)}=b_{i,n+1,a_{n+\frac{1}{2}}}^{(k+1)}-b_{i,n+1,a_{n+\frac{1}{2}}}^{(k)}$ and the auxiliary equations are scaled to be of the same magnitude as the original equations. In (30) the Schur complement has already been reverted, but the value of $b_{i,n+1,a_{n+\frac{1}{2}}}$ depends on $b_{s,n+1}^{(k)}$, so it is better to update both $b_{i,n+1,a_{n+\frac{1}{2}}}$ and $b_{s,n+1}^{(k)}$ every Newton step. The combined update looks like

$$\begin{aligned}&\left( \begin{array}{cc} {\tilde{A}}_{11} &{} {\tilde{A}}_{12}\\ {\tilde{A}}_{21}&{} {\tilde{A}}_{22} \end{array}\right) \begin{pmatrix} {\varDelta }{b}_{s,n+1}^{(k)} \\ {\varDelta }{b}_{i,n+1,{a}_{n+\frac{1}{2}}}^{(k)} \end{pmatrix}= \begin{pmatrix}\tilde{v}_1\\ \tilde{v}_2 \end{pmatrix},\\ {\tilde{A}}_{11}&=\mathcal {E}_{n+1}+\frac{{\varDelta }t}{2} \left[ \partial _{{a}}{\varPhi }_{11,n+1,{a}_{n+\frac{1}{2}}}{b}_{s,n+1}^{(k)}+\partial _{{a}}{\varPhi }_{12,n+1,{a}_{n+\frac{1}{2}}}{b}_{i,n+1,{a}_{n+\frac{1}{2}}}^{(k)}\right. \\&\quad \left. -\,\partial _{{a}}{\varPhi }_{1,{a}_{n+\frac{1}{2}},R_{n+1}}-\frac{{\varDelta }t}{2}\mathcal {D}_{n+1}\right] ^T,\\ {\tilde{A}}_{12}&=\frac{{\varDelta }t}{2} \left[ \partial _{{a}} {\varPhi }_{21,n+1,{a}_{n+\frac{1}{2}}}{b}_{s,n+1}^{(k)}+\partial _{{a}} {\varPhi }_{22,n+1,{a}_{n+\frac{1}{2}}}{b}_{i,n+1,{a}_{n+\frac{1}{2}}}^{(k)} - \partial _{{a}} {\varPhi }_{2,{a}_{n+\frac{1}{2}},R_{n+1}}\right] ^T, \\ {\tilde{A}}_{21}&=\frac{{\varDelta }t}{2}{\varPhi }_{21,n+1,{a}_{n+\frac{1}{2}}},\\ {\tilde{A}}_{22}&=\frac{{\varDelta }t}{2}{\varPhi }_{22,n+1,{a}_{n+\frac{1}{2}}},\\ \tilde{v}_1&=- \left\{ \mathcal {E}_{n+1}{b}_{s,n+1}^{(k)}-\mathcal {E}_n{b}_{s,n}+\frac{{\varDelta }t}{2}\left( g_z(\mathcal {E}_{n}+\mathcal {E}_{n+1}){a}_{n+\frac{1}{2}}\right. \right. \\&\quad -\mathcal {D}_n^T{b}_{s,n}-\mathcal {D}_{n+1}^T{b}_{s,n+1}^{(k)}\\&\quad + \frac{1}{2}{b}_{n,{a}_{n+\frac{1}{2}}}^T\partial _{a}{\varPhi }_{n,{a}_{n+\frac{1}{2}}}{b}_{n,{a}_{n+\frac{1}{2}}}+ \frac{1}{2}{b}_{n+1,{a}_{n+\frac{1}{2}}}^{(k)\;T}\partial _{a}{\varPhi }_{n+1,{a}_{n+\frac{1}{2}}}{b}_{n+1,{a}_{n+\frac{1}{2}}}^{(k)}\\&\quad \left. \left. -\,\partial _{a}{\varPhi }_{{a}_{n+\frac{1}{2}},R_{n}}\cdot {b}_{n,a_{n+\frac{1}{2}}}- \partial _{a}{\varPhi }_{{a}_{n+\frac{1}{2}},R_{n+1}}\cdot {b}_{n+1,a_{n+\frac{1}{2}}}^{(k)}+G^n+G^{n+1}\right) \right\} ,\\ \tilde{v}_2&=-\frac{{\varDelta }t}{2}\left( -{\varPhi }_{2,{a}_{n+\frac{1}{2}},R_{n+1}}+{\varPhi }_{21,n+1,{a}_{n+\frac{1}{2}}}{b}_{s,n+1}^{(k)}+{\varPhi }_{22,n+1,{a}_{n+\frac{1}{2}}}{b}_{i,n+1,{a}_{n+\frac{1}{2}}}^{(k)}\right) . \end{aligned}$$

Substituting (27) into (31) then results in

$$\begin{aligned}&\left( \begin{array}{cc} {\hat{A}}_{11}&{} A_{12} \\ {\hat{A}}_{21}&{} A_{22} \\ \end{array}\right) \begin{pmatrix} {\varDelta }{a}_{n+\frac{1}{2}}^{(k)} \\ {\varDelta }{b}_{i,n,{a}_{n+\frac{1}{2}}}^{(k)} \end{pmatrix}= \begin{pmatrix} {v}_1\\ {v}_2 \end{pmatrix}, \\ {\hat{A}}_{11}&= \mathcal {E}_{n}-\frac{{\varDelta }t}{2} \begin{pmatrix} {C}_{11,n,{a}_{n+\frac{1}{2}}}^{(k)\;T}-{B}_{11,n,{a}_{n+\frac{1}{2}}}^{(k)\;T} \\ {C}_{12,n,{a}_{n+\frac{1}{2}}}^{(k)\;T}-{B}_{12,n,{a}_{n+\frac{1}{2}}}^{(k)\;T} \end{pmatrix}^T \begin{pmatrix} {I} \\ -{F}_{22}^{-1}{F}_{21} \end{pmatrix} +\frac{{\varDelta }t}{2}\mathcal {D}_n,\\ {\hat{A}}_{21}&=-\frac{{\varDelta }t}{2} \begin{pmatrix} {C}_{21,n,{a}_{n+\frac{1}{2}}}^{(k)\;T}-{B}_{21,n,{a}_{n+\frac{1}{2}}}^{(k)\;T} \\ {C}_{22,n,{a}_{n+\frac{1}{2}}}^{(k)\;T}-{B}_{22,n,{a}_{n+\frac{1}{2}}}^{(k)\;T} \end{pmatrix}^T \begin{pmatrix} {I} \\ -{F}_{22}^{-1}{F}_{21}\end{pmatrix}.\\ \end{aligned}$$

This introduces the inverse matrix $F_{22}^{-1}$ for which we introduce the auxiliary equation

$$\begin{aligned} {F}_{22}{u}_{i}=-{F}_{21}{a}. \end{aligned}$$

(32)

We can now write (29) in a form that does not require explicit inverses of matrices

$$\begin{aligned}&\left( \begin{array}{ccc} {\bar{A}}_{11}&{}{{A}}_{12}&{}{\bar{A}}_{13}\\ {\bar{A}}_{21}&{}{{A}}_{22}&{} {\bar{A}}_{23}\\ {\bar{A}}_{31}&{}0&{}{\bar{A}}_{33} \end{array}\right) \begin{pmatrix} {\varDelta }{a}_{n+\frac{1}{2}}^{(k)}\nonumber \\ {\varDelta }{b}_{i,n,{a}_{n+\frac{1}{2}}}^{(k)}\nonumber \\ {\varDelta }{u}_{i,n+\frac{1}{2}}^{(k)} \end{pmatrix}= \begin{pmatrix}{v}_1\\ {v}_2\\ \bar{v}_3 \end{pmatrix}, \nonumber \\ {\bar{A}}_{11}&=\mathcal {E}_{n}-\frac{{\varDelta }t}{2}\left[ {C}_{11,n,{a}_{n+\frac{1}{2}}}^{(k)}-{B}_{11,n,{a}_{n+\frac{1}{2}}}^{(k)}-\mathcal {D}_n\right] ,\nonumber \\ {\bar{A}}_{13}&=-\frac{{\varDelta }t}{2}\left( {C}_{12,n,{a}_{n+\frac{1}{2}}}^{(k)}-{B}_{12,n,{a}_{n+\frac{1}{2}}}^{(k)}\right) ,\nonumber \\ {\bar{A}}_{21}&=-\frac{{\varDelta }t}{2}\left( {C}_{21,n,{a}_{n+\frac{1}{2}}}^{(k)}-{B}_{21,n,{a}_{n+\frac{1}{2}}}^{(k)}\right) ,\nonumber \\ {\bar{A}}_{23}&=-\frac{{\varDelta }t}{2}\left( {C}_{22,n,{a}_{n+\frac{1}{2}}}^{(k)}-{B}_{22,n,{a}_{n+\frac{1}{2}}}^{(k)}\right) ,\nonumber \\ {\bar{A}}_{31}&=\frac{{\varDelta }t}{4}{{F}}_{21},\nonumber \\ {\bar{A}}_{33}&=\frac{{\varDelta }t}{4}{{F}}_{22},\nonumber \\ \bar{v}_3&=-\frac{{\varDelta }t}{4}{{F}}_{21}{a}_{n+\frac{1}{2}}^{(k)}-\frac{{\varDelta }t}{4}{{F}}_{22}{u}_{i,n+\frac{1}{2}}^{(k)}, \end{aligned}$$

(33)

where ${\varDelta }u_{i,n+\frac{1}{2}}^{(k)}=u_{i,n+\frac{1}{2}}^{(k+1)}-u_{i,n+\frac{1}{2}}^{(k)}$. Following the same steps for (30) as for (29) we obtain

$$\begin{aligned}&\left( \begin{array}{cc} {\check{A}}_{11}&{}{\tilde{A}}_{21}^T\nonumber \\ {\check{A}}_{21}&{}{\tilde{A}}_{22} \end{array}\right) ^T\begin{pmatrix} {\varDelta }{b}_{s,n+1}^{(k)}\\ {\varDelta }{b}_{i,n+1,{a}_{n+\frac{1}{2}}}^{(k)} \end{pmatrix}= \begin{pmatrix}\check{v}_1\\ \tilde{v}_2 \end{pmatrix}\nonumber ,\\ {\check{A}}_{11}&=\mathcal {E}_{n+1}+\frac{{\varDelta }t}{2} \begin{pmatrix} {C}_{11,n+1,{a}_{n+\frac{1}{2}}}^{(k)\;T}-{B}_{11,n+1,{a}_{n+\frac{1}{2}}}^{T}\\ {C}_{12,n+1,{a}_{n+\frac{1}{2}}}^{(k)\;T}-{B}_{12,n+1,{a}_{n+\frac{1}{2}}}^{T}\end{pmatrix}^T \begin{pmatrix}{I}\nonumber \\ -{F}_{22}^{-1}{F}_{21}\end{pmatrix} -\frac{{\varDelta }t}{2}\mathcal {D}_{n+1}, \nonumber \\ \check{A}_{12}&=\frac{{\varDelta }t}{2} \begin{pmatrix} {C}_{21,n+1,{a}_{n+\frac{1}{2}}}^{(k)\;T}-{B}_{21,n+1,{a}_{n+\frac{1}{2}}}^{T}\\ {C}_{22,n+1,{a}_{n+\frac{1}{2}}}^{(k)\;T}-{B}_{22,n+1,{a}_{n+\frac{1}{2}}}^{T}\end{pmatrix}^T \begin{pmatrix}{I}\\ -{F}_{22}^{-1}{F}_{21}\end{pmatrix},\nonumber \\ \check{v}_1&=- \left\{ \mathcal {E}_{n+1}{b}_{s,n+1}^{(k)}-\mathcal {E}_n{b}_{s,n}+\frac{{\varDelta }t}{2}\left( g_z(\mathcal {E}_{n}+\mathcal {E}_{n+1}){a}_{n+\frac{1}{2}}-\mathcal {D}_n^T{b}_{s,n}\right. -\mathcal {D}_{n+1}^T{b}_{s,n+1}^{(k)}\right. \nonumber \\&\quad +\,\frac{1}{2}\begin{pmatrix}{I}\\ -{F}_{22}^{-1}{F}_{21}\end{pmatrix}^T{C}^T_{n,{a}_{n+\frac{1}{2}}}{b}_n+ \frac{1}{2}\begin{pmatrix}{I}\\ -{F}_{22}^{-1}{F}_{21}\end{pmatrix}^T{C}_{n+1,{a}_{n+\frac{1}{2}}}^{(k)\;T}{b}_{n+1,a_{n+\frac{1}{2}}}^{(k)}\nonumber \\&\quad \left. \left. +\,\begin{pmatrix} -I\\ {F}_{22}^{-1}{F}_{21}\end{pmatrix}^T{B}^T_{n,{a}_{n+\frac{1}{2}}}{b}_n+ \begin{pmatrix}-I\\ {F}_{22}^{-1}{F}_{21}\end{pmatrix}^T{B}^T_{n+1,{a}_{n+\frac{1}{2}}}{b}_{n+1,a_{n+\frac{1}{2}}}^{(k)}+G^n+G^{n+1}\right) \right\} . \end{aligned}$$

(34)

We again introduce an auxiliary equation

$$\begin{aligned} {F}_{22}{d}=-\left( ({C}_{12}-2{B}_{12})^Tb_s + ({C}_{22}-2{B}_{22})^Tb_{i}\right) . \end{aligned}$$

This allows us to remove the inverses from (34), resulting in

$$\begin{aligned}&\left( \begin{array}{ccc} \breve{A}_{11}&{}{\tilde{A}}_{21}^T&{}{{\breve{A}}}_{13}\\ {\breve{A}}_{21}&{}{\tilde{A}}_{22}&{}{{\breve{A}}}_{23}\\ {{\breve{A}}}_{31}&{}0&{}{{\breve{A}}}_{33} \end{array}\right) ^T\begin{pmatrix} {\varDelta }{b}_{s,n+1}^{(k)}\\ {\varDelta }{b}_{i,n+1,{a}_{n+\frac{1}{2}}}^{(k)}\\ {\varDelta }{d}_{n+1,{a}_{n+\frac{1}{2}}}^{(k)} \end{pmatrix}= \begin{pmatrix}\breve{v}_1\\ \tilde{v}_2\\ \breve{v}_3 \end{pmatrix},\nonumber \\ {{\breve{A}}}_{11}&=\mathcal {E}_{n+1}+\frac{{\varDelta }t}{2}\left[ {C}_{11,n+1,{a}_{n+\frac{1}{2}}}^{(k)}-{B}_{11,n+1,{a}_{n+\frac{1}{2}}} -\mathcal {D}_{n+1}\right] , \nonumber \\ {{\breve{A}}}_{13}&=\frac{{\varDelta }t}{2}\left( {C}_{12,n+1,{a}_{n+\frac{1}{2}}}^{(k)}-{B}_{12,n+1,{a}_{n+\frac{1}{2}}}\right) ,\nonumber \\ {{\breve{A}}}_{21}&=\frac{{\varDelta }t}{2}\left( {C}_{21,n+1,{a}_{n+\frac{1}{2}}}^{(k)}-{B}_{21,n+1,{a}_{n+\frac{1}{2}}}\right) ,\nonumber \\ {{\breve{A}}}_{23}&=\frac{{\varDelta }t}{2}\left( {C}_{22,n+1,{a}_{n+\frac{1}{2}}}^{(k)}-{B}_{22,n+1,{a}_{n+\frac{1}{2}}}\right) ,\nonumber \\ {{\breve{A}}}_{31}&=\frac{{\varDelta }t}{4}{F}_{21},\nonumber \\ {{\breve{A}}}_{33}&=\frac{{\varDelta }t}{4}{F}_{22},\nonumber \\ {\breve{v}}_1&=- \left\{ \mathcal {E}_{n+1}{b}_{s,n+1}^{(k)}-\mathcal {E}_n{b}_{s,n}+\frac{{\varDelta }t}{2}\left( g_z(\mathcal {E}_{n}+\mathcal {E}_{n+1}){a}_{n+\frac{1}{2}}-\mathcal {D}_n^T{b}_{s,n}\right. \right. \nonumber \\&\quad -\mathcal {D}_{n+1}^T{b}_{s,n+1}^{(k)}+ \frac{1}{2}{C}^T_{11,n,{a}_{n+\frac{1}{2}}}{b}_{s,n}+ \frac{1}{2}{C}^T_{21,n,{a}_{n+\frac{1}{2}}}{b}_{i,n,{a}_{n+\frac{1}{2}}}\nonumber \\&\quad + \frac{1}{2}{C}_{11,n+1,{a}_{n+\frac{1}{2}}}^{(k)\;T}{b}_{s,n+1}^{(k)}+ \frac{1}{2}{C}_{21,n+1,{a}_{n+\frac{1}{2}}}^{(k)\;T}{b}_{i,n+1,{a}_{n+\frac{1}{2}}}^{(k)}\nonumber \\&\quad + {F}_{21}^T{d}_{n,{a}_{n+\frac{1}{2}}}+{F}_{21}^T{d}_{n+1, {a}_{n+\frac{1}{2}}}^{(k)}\nonumber \\&\quad -{B}^T_{11,n,{a}_{n+\frac{1}{2}}}{b}_{s,n}- {B}^T_{21,n,{a}_{n+\frac{1}{2}}}{b}_{i,n,{a}_{n+\frac{1}{2}}}\nonumber \\&\quad \left. \left. -\,{B}^T_{11,n+1,{a}_{n+\frac{1}{2}}}{b}_{s,n+1}^{(k)}- {B}^T_{21,n+1,{a}_{n+\frac{1}{2}}}{b}_{i,n+1,{a}_{n+\frac{1}{2}}}^{(k)}+G^n+G^{n+1}\right) \right\} ,\nonumber \\ \breve{v}_3&=-\frac{{\varDelta }t}{4}\left( \left( {C}_{12,n+1,{a}_{n+\frac{1}{2}}}^{(k)}-2{B}_{12,n+1,{a}_{n+\frac{1}{2}}}\right) {b}_{s,n+1}^{(k)}\right. \nonumber \\&\quad \left. +\,\left( {C}_{22,n+1,{a}_{n+\frac{1}{2}}}^{(k)}-2{B}_{22,n+1,{a}_{n+\frac{1}{2}}}\right) {b}_{i,n+1,{a}_{n+\frac{1}{2}}}^{(k)}+ {F}_{22}{d}_{n+1,{a}_{n+\frac{1}{2}}}^{(k)}\right) , \end{aligned}$$

(35)

where ${\varDelta }d_{n+1,{a}_{n+\frac{1}{2}}}^{(k)}=d_{n+1,{a}_{n+\frac{1}{2}}}^{(k+1)}-d_{n+1,{a}_{n+\frac{1}{2}}}^{(k)}$.

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

$$\begin{aligned} \int _{\mathcal {D}_{t,s,h}^l}\eta _i^t\eta _j^t\,\mathrm {d}x\,\mathrm {d}y=\int _{\mathcal {D}_{t,s,h}^l}\begin{pmatrix}0\\ 0\\ \eta _i^t\eta _j^t\end{pmatrix}\cdot e_{ z}\,\mathrm {d}x\,\mathrm {d}y. \end{aligned}$$

We can now apply (7.2) from [7] to find

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\sum _{l=1}^m\int _{\mathcal {D}_{t,s,h}^l}\eta _i^t\eta _j^t\,\mathrm {d}x\,\mathrm {d}y&= \sum _{l=1}^m\left( \int _{\mathcal {D}_{t,s,h}^l}\nabla \cdot \begin{pmatrix}0\\ 0\\ \eta _i^t\eta _j^t\end{pmatrix}(v\cdot e_{ z})\,\mathrm {d}x\,\mathrm {d}y \right. \nonumber \\&\quad +\int _{\mathcal {D}_{t,s,h}^l}\frac{\partial }{\partial t}\left( \begin{pmatrix}0\\ 0\\ \eta _i^t\eta _j^t\end{pmatrix}\cdot e_{ z}\right) \,\mathrm {d}x\,\mathrm {d}y\\&\quad \left. +\,\int _{\partial \mathcal {D}_{t,s,h}^l}\left( \begin{pmatrix}0\\ 0\\ \eta _i^t\eta _j^t\end{pmatrix}\times v\right) \cdot \tau \,\mathrm {d}l\right) .\nonumber \end{aligned}$$

(36)

Since $\eta ^t_j$ does not depend on z, the first term in (36) vanishes. In addition, the internal contributions from two adjacent patches $\mathcal {D}_{t,s,h}^l$ 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=\frac{\partial R}{\partial t}$ and we have applied the correction (14) to ensure that the free surface moves tangentially to the solid wall. These considerations reduce (36) to

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\sum _{l=1}^m\int _{\mathcal {D}_{t,s,h}^l}\eta _i^t\eta _j^t\,\mathrm {d}x\,\mathrm {d}y&= \sum _{l=1}^m\int _{\mathcal {D}_{t,s,h}^l}\frac{\partial \eta _i^t}{\partial t}\eta _j+\frac{\partial \eta _j^t}{\partial t}\eta _i\,\mathrm {d}x\,\mathrm {d}y\\&\quad +\int _{\partial {\varOmega }_{t,R,h}\cap \partial {\varOmega }_{t,s,h}}\eta _i\eta _j\left[ \left( (e_{ z}\cdot \tau )\tau +(e_{ z}\cdot t_R)t_R\right) \times \frac{\partial R}{\partial t}\right] \cdot \tau \,\mathrm {d}l. \end{aligned}$$

From here we use $t_R=\tau \times \nu _R$ and expand all products to arrive at

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\mathcal {E}=\mathcal {D}+\mathcal {D}^T+\mathcal {E}_R. \end{aligned}$$

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions

About this article

Cite this article

Brink, F., Izsák, F. & van der Vegt, J.J.W. Hamiltonian Finite Element Discretization for Nonlinear Free Surface Water Waves. J Sci Comput 73, 366–394 (2017). https://doi.org/10.1007/s10915-017-0416-9

Download citation

Received: 11 August 2016
Revised: 06 March 2017
Accepted: 09 March 2017
Published: 23 March 2017
Issue Date: October 2017
DOI: https://doi.org/10.1007/s10915-017-0416-9

Hamiltonian Finite Element Discretization for Nonlinear Free Surface Water Waves

Abstract

Similar content being viewed by others

The Pressure Boundary Condition and the Pressure as Lagrangian for Water Waves

Computation of Nonlinear Free-Surface Flows Using the Method of Fundamental Solutions

Free-Boundary Problems for Wave–Structure Interactions in Shallow-Water: DG-ALE Description and Local Subcell Correction

1 Introduction

2 Governing Equations

2.1 Lagrangian and Hamiltonian Approach

3 Discretization of the Variational Principle

Proposition 1

Lemma 1

Proof

Lemma 2

Proof

Theorem 1

Proof

Remark 1

4 Time Integration for the Discrete Variational Formulation

4.1 Computing Derivatives \(\partial _{a}{\varPhi }\) and \(\partial _{a}{\varPhi }_R\)

4.2 Mesh Deformation Algorithm

5 Results

5.1 Fenton and Rienecker Wave

5.2 Soliton

5.3 Comparison with Experiments

References

Acknowledgements

Author information

Authors and Affiliations

Corresponding author

Appendices

Appendix 1: Solution of Algebraic Equations

Appendix 2: Time Derivative of the Mass Matrix

Rights and permissions

About this article

Cite this article

Keywords

Mathematics Subject Classification

Navigation

Hamiltonian Finite Element Discretization for Nonlinear Free Surface Water Waves

Abstract

Similar content being viewed by others

The Pressure Boundary Condition and the Pressure as Lagrangian for Water Waves

Computation of Nonlinear Free-Surface Flows Using the Method of Fundamental Solutions

Free-Boundary Problems for Wave–Structure Interactions in Shallow-Water: DG-ALE Description and Local Subcell Correction

1 Introduction

2 Governing Equations

2.1 Lagrangian and Hamiltonian Approach

3 Discretization of the Variational Principle

Proposition 1

Lemma 1

Proof

Lemma 2

Proof

Theorem 1

Proof

Remark 1

4 Time Integration for the Discrete Variational Formulation

4.1 Computing Derivatives \(\partial _{a}{\varPhi }\) and \(\partial _{a}{\varPhi }_R\)

4.2 Mesh Deformation Algorithm

5 Results

5.1 Fenton and Rienecker Wave

5.2 Soliton

5.3 Comparison with Experiments

References

Acknowledgements

Author information

Authors and Affiliations

Corresponding author

Appendices

Appendix 1: Solution of Algebraic Equations

Appendix 2: Time Derivative of the Mass Matrix

Rights and permissions

About this article

Cite this article

Share this article

Keywords

Mathematics Subject Classification

Search

Navigation