Isogeometric analysis of linear free-surface potential flow

TL;DR

This paper introduces a new isogeometric analysis formulation for linear free-surface potential flow, achieving higher accuracy and energy conservation compared to traditional methods, demonstrated through numerical simulations.

Contribution

It develops a monolithic variational formulation compatible with IGA, integrating interior and free-surface problems for improved simulation accuracy.

Findings

Exact energy conservation achieved.

Superior accuracy over traditional finite elements.

Virtually no dispersion error on coarse meshes.

Abstract

This paper presents a novel variational formulation to simulate linear free-surface flow. The variational formulation is suitable for higher-order finite elements and higher-order and higher-continuity shape functions as employed in Isogeometric Analysis (IGA). The novel formulation combines the interior and free-surface problems in one monolithic formulation. This leads to exact energy conservation and superior performance in terms of accuracy when compared to a traditional segregated formulation. This is confirmed by the numerical computation of traveling waves in a periodic domain and a three-dimensional sloshing problem. The isogeometric approach shows significant improved performance compared to traditional finite elements. Even on very coarse quadratic NURBS meshes the dispersion error is virtually absent.

Equations

$$\Delta \varphi =$$

$$\displaystyle\phi_{+}\mbox{$\frac{1}{2}$}\nabla\phi\cdot\nabla\phi+g\eta=$$

$${\eta }_t+\tilde{\nabla }\varphi \tilde{\nabla }\eta -{\varphi }_z=$$

$$\displaystyle\nabla\cdot\mbox{\boldmath{$u$}}=$$

$$\displaystyle\mbox{\boldmath{$u$}}=$$

$$\displaystyle\mbox{\boldmath{$u$}}\cdot\mbox{\boldmath{$n$}}=\mbox{\boldmath{$n$}}\cdot\nabla\phi=$$

$$\Delta \varphi =$$

$${\varphi }_t+g\eta =$$

$${\varphi }_z=$$

$$\displaystyle\mbox{\boldmath{$n$}}\cdot\nabla\phi=$$

$$\Delta \varphi =$$

$${\varphi }_{tt}+g{\varphi }_z=$$

$$\displaystyle\mbox{\boldmath{$n$}}\cdot\nabla\phi=$$

$$(w,v{)}_{\Omega }=$$

$$(w,v{)}_{{\Gamma }_{fs}}=$$

$$\Vert w{\Vert }_{\Omega }^2=$$

$$\Vert w{\Vert }_{{\Gamma }_{fs}}^2=$$

$$B_r(w,\varphi )=0$$

$$B_r(w,\varphi ):=(\nabla w,\nabla \varphi {)}_{\Omega }+\frac{1}{g}(w,{\varphi }_{tt}{)}_{{\Gamma }_{fs}}.$$

$$B_r({\varphi }_t,\varphi )=$$

$$=$$

$$E_{\mathrm{kin}}=$$

$$E_{\mathrm{pot}}=$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{E}_{\mathrm{kin}}+\frac{\mathrm{d}}{\mathrm{d}t}{E}_{\mathrm{pot}}=0.$$

$$\Delta \varphi =$$

$$\varphi =$$

$$n\cdot \nabla \varphi =$$

$${\hat{\varphi }}_t+g\eta =$$

$${\eta }_t=$$

$$B_{int}(w,\varphi )=0$$

$$B_{int}(w,\varphi ):=(\nabla w,\nabla \varphi ).$$

$$B_{fs}(\{w,v\};\{\varphi ,\eta \})=\frac{g^2}{{\alpha }^2}(v,{\varphi }_z{)}_{{\Gamma }_{fs}},$$

$$B_{fs}(\{w,v\};\{\varphi ,\eta \}):=(w,{\varphi }_t+g\eta {)}_{{\Gamma }_{fs}}+\frac{g^2}{{\alpha }^2}(v,{\eta }_t{)}_{{\Gamma }_{fs}}.$$

$$(\nabla w,\nabla \varphi {)}_{\Omega }+(w,\lambda {)}_{{\Gamma }_{fs}}+(q,\varphi {)}_{{\Gamma }_{fs}}=(q,\hat{\varphi }{)}_{{\Gamma }_{fs}}.$$

$$(\nabla {\varphi }_t,\nabla \varphi )+({\varphi }_t,\lambda )=0,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{E}_{\mathrm{kin}}=-({\varphi }_t,\lambda ).$$

$$0=$$

Full text

The novel formulation combines the interior and free-surface problems in one monolithic formulation. This leads to exact energy conservation and superior performance in terms of accuracy when compared to a traditional segregated formulation. This is confirmed by the numerical computation of traveling waves in a periodic domain and a three-dimensional sloshing problem. The isogeometric approach shows significant improved performance compared to traditional finite elements. Even on very coarse quadratic NURBS meshes the dispersion error is virtually absent.

1 Introduction

In maritime and offshore applications potential flow is often a very useful model for prescribing relevant design scenarios. This is particularly the case when waves, which are dominated by inertia effects, are of main concern.

Numerical methods for predicting potential flow are often based on boundary integral formulations, referred to as boundary element methods (BEM) or panel methods. Examples of successful application of the BEM are (i) wave drift forces [1], (ii) ship wave resistance [2] and (iii) ship motion and maneuvering [3]. These references indicate the starting point of the development of BEM and numerous works have been published to improve these methods.

All these methods exploit the Greens functions and only require the boundary to be discretized. This leads to very low degree-of-freedom counts, however, the resulting matrices are dense and expensive to assemble. Moreover, straightforward implementations have a very unfavorable scaling in terms of computational time versus degree-of-freedom.

One way to circumvent this scaling issue is the use of the fast multipole method [4]. An alternative approach, which is explored in this paper, is to forgo the boundary integral approach. Instead variational methods discretizing the entire volume are used, this encompasses the classical finite element method (FEM) as well as the novel Isogeometric Analysis (IGA) approach [5, 6]. IGA uses splines as shape functions as opposed to the more standard polynomials in FEM. These spline functions are also used in CAD. One of the main goals of IGA is to create a seamless integration of analysis into design processes by using the same geometry description. It turns out that the spline shape functions have several other benefits, in terms of efficiency and accuracy [7, 8, 9].

Linear finite elements are used for the simulation of water waves by for instance Wu and Eatock-Taylor [10], Kim and Bai et al [11], and Westhuis [12], just to name a few. In the first two publications it is the starting point of a large body of work. An important issue in these methods is the coupling of the interior Laplace problem with the free-surface evolution, in particular the horizontal velocity at the free-surface. The mentioned bodies of work deal with this issue in different manners. In Wu and Eatock-Taylor [10] the recovered velocity is a projection on a new finite element space. This is referred to as Galerkin projection and results in a mixed formulation. Kim and Bai [11] employ a similar potential-velocity mixed formulation but employ the Hamiltonian formalism. In Westhuis [12] the vertical velocity at the free-surface is reconstructed using finite differences. Zienkiewicz and Zhu [13] have used other reconstruction methods leading to adequate results.

The previous works do allow for nonlinear waves. In this work we take a step back and limit ourselves to linear waves. However, we perform a thorough investigation of the numerical behavior of the methods considered with a particular focus on the coupling between the interior and free-surface problem. We investigate several weak formulations that deal with the coupling in different manners. Moreover, we establish the well-posedness of these formulations. In the discrete case, we derive accuracy estimates and analyze the energy behavior. Using an appropriate time integrator for some weak formulations provides exact energy conservation. The spatial discretization employs higher order elements.

The outline of the paper is as follows. In section 2 the strong formulation of the problem is introduced. In section 3 several corresponding weak formulations are presented and their energy conservation properties are analyzed. Energy conservation is only guaranteed when an appropriate time integrator is used. In section 4 midpoint time integration is shown to achieve this. After the temporal discretisation, we will focus on the spatial discretisation in section 5. In this section particular attention is given to Isogeometric analysis. In section 6 the energy conservation and dispersion properties of the methods are investigated. This reveals the superiority of the monolithic approach over the segregated approach. In section 7 and 8, the convergence of monolithic formulation is analyzed and verified, respectively. Finally, the paper ends with a 3D show case problem in section 9, and conclusions in section 10.

2 Strong forms of the free-surface problem

In this section we introduce three different sets of governing equations describing potential waves. We start off with the nonlinear equations, which are then linearized and condensed to arrive at two alternative formulations.

2.1 Non-linear strong form

The governing equations read:

[TABLE]

where $\phi:\Omega\rightarrow\mathbb{R}$ is the potential, $\eta:{\Gamma_{\!\!f\!s}}\rightarrow\mathbb{R}$ is the water height, and $\widetilde{\nabla}$ is the gradient restricted to exclude the vertical direction. Here $\Gamma_{fs}$ denotes the free-surface.

The Laplace problem (1a) is a consequence of the conservation of mass for a constant density fluid,

[TABLE]

in conjunction with the assumption of an inviscid and irrotational fluid, for which we can write

[TABLE]

The dynamic boundary condition (1b) is a condition on the pressure, whereas the kinematic condition (1c) ensures that the surface moves with the water.

The problem needs to be augmented with boundary conditions on the remaining boundaries. In this work we will assume either periodic boundary conditions or no-penetration boundary conditions. In the first case no boundary term is present, as such no boundary condition needs to be enforced. In case of the no-penetration boundary conditions, we get a Neumann boundary condition

[TABLE]

The problem is completed by specifying appropriate initial conditions for $\phi$ and $\eta$ at the free-surface. Note that in the fully non-linear problem, the location of the free-surface, denoted with ${\Gamma_{\!\!f\!s}}$, is part of the solution as it is determined by $z=\eta$.

2.2 Linear strong form

The nonlinear equations (1) can be simplified by assuming small disturbances. By neglecting the quadratic terms we arrive at

[TABLE]

where ${\Gamma_{\!\!f\!s}}$ is now assumed frozen on the undisturbed location.

By combining the dynamic and kinematic boundary condition, (5b) and (5c), we can eliminate the water height $\eta$ from the problem. The problem reduces to:

[TABLE]

where the problem is now second order in time. This means that the initial condition for $\eta$ at the free-surface is replaced by an initial condition for $\phi_{t}$ at the free-surface.

3 Weak forms of the linear free-surface problem

In this section we present several weak formulations of the linear wave problem. We start with a weak form of (6), and subsequently propose several formulations for problem (5). These later formulations of mixed character are more amiable in the nonlinear case or in situations with currents present. Furthermore, we analyze the energy conservation of the solution for each of the formulations. An analysis of the existence and accuracy of the solution for each of the formulations, is postponed to section 7.

We introduce the notation

[TABLE]

to denote the innerproduct over the entire domain and free-surface, respectively. The corresponding norms are denoted with

[TABLE]

3.1 Weak formulation of the reduced problem

A weak form of (6) follows when multiplying (6a) with an arbitrary function and integrating over the domain. The order of the required derivatives can be reduced by using Green’s identities and the resulting boundary terms can be simplified by using (6b) and (6c). Let $\mathcal{W}=H^{1}(\Omega)$ denote the function-space. The variational formulation of the reduced problem reads:

Find $\phi\in\mathcal{W}$ such that for all $w\in\mathcal{W}$:

[TABLE]

where

[TABLE]

An energy conservation statement for the reduced form can be derived by choosing $w=\phi_{t}$ in (9). This selection gives

[TABLE]

Realizing that the definition for kinetic and potential are

[TABLE]

we arrive at the following statement

[TABLE]

Consequently, the reduced weak form is exactly energy conservative assuming appropriate time integration.

3.2 Weak formulations of the segregated problem

Here we present a segregated weak formulation serving as a reference method based on the work of Wu and Taylor [10] and Wu et al. [14], Westhuis [12], Kyoung et al. [15], Bai et al. [16] and Kim et al. [17]. It is based on the strong form (5) and decouples the interior and surface parts of the problem.

The interior problem in strong form reads:

[TABLE]

where $\hat{\phi}$ is input from the free–surface problem:

[TABLE]

Here is $\phi_{z}$ is given by the interior problem. As such the two problems are artificially decoupled.

Let $\mathcal{W}_{0}$ and $\mathcal{W}_{\hat{\phi}}$ denote the subspaces of $H^{1}(\Omega)$ satisfying homogeneous Dirichlet boundary condition, $\phi=0$ on ${\Gamma_{\!\!f\!s}}$, and the homogeneous Dirichlet boundary condition $\phi=\hat{\phi}$ on ${\Gamma_{\!\!f\!s}}$, respectively. Furthermore let $\mathcal{W}_{\Gamma_{\!\!f\!s}}$ denote the trace space of $\mathcal{W}$. The associated weak formulation for the interior problem reads:

Find $\phi\in\mathcal{W}_{\hat{\phi}}$ such that for all $w\in\mathcal{W}_{0}$:

[TABLE]

where

[TABLE]

For the free-surface the weak form reads:

Find $(\phi,\eta)\in\mathcal{W}_{\Gamma}\times\mathcal{W}_{\Gamma}$ such that for all $(w,v)\in\mathcal{W}_{\Gamma}\times\mathcal{W}_{\Gamma}$:

[TABLE]

where

[TABLE]

Here $\alpha$ is a parameter, that eventually will be chosen to depend on the time integrator.

3.2.1 Energy conservation of the segregated formulation

To arrive and an energy statement we would like to select $w=\phi_{t}$ in both the interior problem (16) and free-surface problem (18). However, due to the Dirichlet boundary condition on ${\Gamma_{\!\!f\!s}}$ this is not allowed for the interior problem. To remedy this, the weak form (16) is written in an equivalent Lagrange multiplier formulation. from which yields exactly the same solution if the Lagrange multiplier space is constructed appropriately. This Lagrange multiplier formulation reads as follows:

Find $(\phi,\lambda)\in\mathcal{W}\times\mathcal{W}_{\Gamma}$ such that for all $(w,q)\in\mathcal{W}\times\mathcal{W}_{\Gamma}$:

[TABLE]

The new formulation allows the selection of $w=\phi_{t}$ and additionally we set $q=0$. This yields:

[TABLE]

which can be rewritten as:

[TABLE]

For the free-surface problem we select $w=\phi_{z}$ and $v=\eta\alpha^{2}/g$ which results in:

[TABLE]

This can be rewritten as:

[TABLE]

Combining the interior and boundary problem we arrive at the following energy statement:

[TABLE]

which results in conservation of energy if $\lambda=-\phi_{z}$. Unfortunately this relation only holds on sufficiently smooth meshes, where the solution is converged. This means conservation of energy can not be guaranteed.

In the following two subsections two alternative approaches to remedy this energy error are discussed. The first approach is to reconstruct the Lagrange multiplier and use this as a forcing in the boundary formulation. The other approach is to solve the boundary and interior problem in one monolithic formulation.

3.3 Segregated formulation, with LM reconstruction

The Lagrange multiplier from equation (20) can be reconstructed by selecting $q=0$, which results in

[TABLE]

for all $w\in\mathcal{W}$.

Energy conservation can be recovered if this Lagrange multiplier is directly used in the weak form for the free-surface problem given in (18). The right-hand side would need to be modified to:

[TABLE]

where in the last integral the functions $v$ need to be arbitrarily extended into the domain.

For the modified free-surface problem we again select $w=\phi_{z}$ and $v=\eta\alpha^{2}/g$ which results in:

[TABLE]

This can be combined with already obtained kinetic energy statement (22) to yield:

[TABLE]

and as such recovering the energy conservation.

Remark 1

Note that in the work of Westhuis [12] a higher-order reconstruction of $\phi_{z}$ is used. This reconstruction does not provide exact energy conservation.

3.4 Monolithic weak formulation

Here we present a monolithic formulation that bypasses the need of the Lagrange multiplier construction. When the divergence theorem on the interior problem creates boundary terms the kinematic boundary condition can be substituted in, this eliminates the problematic derivative. The dynamic boundary condition is added to the weak form in a way that guarantees coercivity.

We propose the weak formulation:

Find $(\phi,\eta)\in\mathcal{W}\times\mathcal{W}$ such that for all $(w,v)\in\mathcal{W}\times\mathcal{W}$,

[TABLE]

with

[TABLE]

Again $\alpha$ is a parameter that eventually will be chosen based on the time integrator.

3.4.1 Conservation of energy of the monolithic form

To establish an energy statement we select $w=\phi_{t}$ and $v=2\eta_{t}-\frac{\alpha}{g}\phi_{t}$ and substitute this in (3.4) which gives:

[TABLE]

From this it follows directly that

[TABLE]

which indicates that the total energy is conserved.

4 Time integration

In sections 3.1 and 3.4.1 it was proven for the respective weak forms that the total energy is conserved. This is stated as

[TABLE]

When defining the energies as

[TABLE]

and assuming the time-integrator has the correct behavior the time continous statement can be translated to a time discrete conservation statement, namely

[TABLE]

This translation holds when the time integrator satisfies the following relations,

[TABLE]

Here the tilde denotes the value used by the time integrator to evaluate the weak formulation.

In this paper midpoint time integration is adopted for both the first and second-order problems. In both cases the required translation is valid, as will be shown, while having second-order convergence.

Note, that as shown in [18] generalized midpoint time integration almost satisfies (37). In addition to (37) it also features a diffusion that scales with $\mathcal{O}(\Delta t)$.

4.1 Midpoint time integration for a first-order problem

For problems with first-order time derivatives the midpoint time-integration is determined by the following relations.

The ordinary differential equation:

[TABLE]

• 2.

The interpolation relation:

[TABLE]

• 3.

The kinematic relation:

[TABLE]

This results in three relations for three unknowns that are in principle solvable.

In 7.3 and B it will be shown that using

[TABLE]

results in favourable properties of the formulations. For the midpoint time-integration this results in:

[TABLE]

4.1.1 Kinetic and potential energy behavior

To see whether the kinetic energy requirement (37a) holds, we rewrite the kinematic relation (40) as:

[TABLE]

Using this relation we find that:

[TABLE]

This demonstrates that the kinetic energy requirement is satisfied. The potential energy relation goes analogously.

4.2 Midpoint time-integration for a 2nd order problem

For problems with 2nd order time derivatives the midpoint time-integration is determined by the following relations.

The ordinary differential equation

[TABLE]

• 2.

Two interpolation relations,

[TABLE]

• 3.

Two kinematic relations,

[TABLE]

This results in five relations for five unknowns that are in principle solvable.

4.2.1 Kinetic energy behavior

Again, to see whether the kinetic energy requirement (37a) holds we rewrite the kinematic relation (47a) as

[TABLE]

Using this relation we again find that

[TABLE]

demonstrating that the kinetic energy requirement is satisfied.

4.2.2 Potential energy behavior

For the potential energy we need to make the following substitution,

[TABLE]

which is justified by the dynamics boundary condition (5c). The potential energy requirement (37b) holds, namely,

[TABLE]

5 Spatial discretization

In this section the spatial discretisation is discussed. For the spatial discretisation both finite elements (FE) and NURBS based Isogeometric analysis (IGA) are used. IGA can be seen as an extension of FE. Both methods will be explained separately. But first the commonalities are discussed. Both approaches approximate the unknown exact solution, denoted as \phi(\mbox{\boldmath{x}}), by a weighted sum of known shape functions N_{b}(\mbox{\boldmath{x}}):

[TABLE]

Here $\phi_{b}$ are the unknown parameters that need to be determined. To convert the weak form in a set of equations we select:

[TABLE]

Using these approximations the terms in the weak form can be determined. For instance, a term in the weak form such as $(w,\phi)$ results in:

[TABLE]

where we dropped the argument $x$ to simplify the notation. The matrix is defined as:

[TABLE]

Here $\Omega_{e}$ denotes the element domain,of which the union covers the entire domain, viz. $\Omega=\cup_{e}\Omega_{e}$. The element integral is approximated using Gauss quadrature. Both the domain of the integral and gradients incorporate the mapping between the physical and reference domain, following the isoparametric paradigm.

Note that in contrast to the boundary element method the global matrices in FE and IGA are usually very sparse due to the compact support of the shape functions. An entry in $M_{ab}$ is only non-zero if the associated shape functions $N_{a}$ and $N_{b}$ share at least one element where both functions are non-zero.

5.1 Finite elements

In the case of finite elements the shape functions are often simple piecewise polynomials that are continuous at element interfaces.

The simplest example being the 1D linear shape functions as depicted in figure 1.

Higher-order shape functions, such as the quadratic Lagrangian shape functions depicted in figure 2, are also available.

The expected convergence of the overall numerical method is often directly related to the polynomial order of the shapefunctions used. Note that the functions are still continuous at element interfaces, but that higher order derivatives are not.

5.2 Isogeometric analysis

Isogeometric analysis is an extension of finite elements attempting to bridge the gap between design (CAD) and analysis (FE). The idea is to use the same NURBS (Non-Uniform Rational B-Splines) of CAD for the analysis. This results in an exact representation of CAD geometries, such as circles and ellipses. These shape functions allow for higher-order continuity at element interfaces, see Figure 3.

The higher-order continuity leads to additional constraints, resulting in higher-order approximations with limited degrees-of-freedom. Additionally, the higher-order continuity allows a strict tailoring of approximations spaces. For the incompressible Navier-Stokes equations this can be exploited to create velocity-pressure approximations that result in exactly solenoidal solutions [19]. In [20, 21] this exact divergence is essential for getting correct energy behavior of the single and two-fluid Navier-Stokes problem, respectively. In the current context the improved approximation and spectral properties are of paramount importance.

In [7] it is demonstrated that the extra resolution provided by the element kinks are not beneficial. Evans et al. [9] show that these kinks lead to bad performance of higher-order finite elements in wave propagation problems. Cotrell et al. [8] formally prove that NURBS have optimal approximation properties. This superior behavior is noticed in numerous application areas.

Due to the success of IGA on the one hand and the specific requirements on the other hand, a multitude of alternative spline technologies have emerged. Examples are T-Splines [22], LR-Splines [23], U-Splines [24] and many others. These approaches all allow higher-order continuity at element interfaces but avoid the rigid tensorial construction required for NURBS. Assuring the shape functions are linearly independent and the resulting system matrices are solvable is one of the main issues of these alternative splines.

6 Numerical comparison of the formulations

In this section the different weak formulations are numerically investigated. To this purpose we present a traveling wave case. We compare the energy behavior and assess convergence under mesh refinement.

6.1 Two-dimensional traveling wave

The performance of the different formulations is assesses using a simple traveling wave in a two-dimensional periodic domain.

Periodic boundary conditions are enforced on the sides ($\Gamma_{l}$ and $\Gamma_{r}$), no penetration on the bottom ($\Gamma_{b}$) and a free-surface on the top (${\Gamma_{\!\!f\!s}}$), as depicted in figure 4.

The initial condition is specified to be an airy wave. For a given wave height $\xi$ the water elevation and flow potential are given as:

[TABLE]

with $x$ and $z$ representing the spatial coordinates in horizontal and vertical direction. The wavenumber $k=2\pi/\lambda$ with wavelength $\lambda$, and angular frequency $\omega$ are related through the dispersion relation:

[TABLE]

where $H$ is the water depth. Snapshots for the solution of an airy-wave with a wavelength equal to the size of the domain are given in Figure 5.

6.2 Mesh convergence of energy traces

The simulation are done on a unit size square. The time step is chosen in alignment with the mesh resolution:

[TABLE]

where $T_{ex}$ is the theoretically expected period of the wave.

A straightforward computation reveals that the relations for the wave energy in this domain are:

[TABLE]

Figures 6 and 7 show the convergence of the energy time trace for linear finite elements when employing the segregated and monolithic formulation.

The energy traces for the reduced formulation are virtually identical to those of the monolithic formulation and are therefore not plotted. The segregated formulation displays quite severe fluctuations in the energy components and the sum of both does not remain constant. These fluctuations disappear with mesh refinement indicating the method is in principle valid. The monolithic formulation gives much better results on the same mesh when compared with the segregated formulation. The fluctuations of the kinetic and potential energy are significantly less than for the segregated formulation. On the 12 x 12 mesh the fluctuations are barely visible and on the 24 x 24 mesh they are essentially gone. Moreover, the total energy stays perfectly conserved, even on the coarsest mesh. Note that on the coarser meshes the total energy is underestimated. This is due to the loss of energy in projecting the initial condition on to the discrete space. This mismatch disappears under mesh refinement.

Figures 8 and 9 show a similar comparison but for quadratic finite elements. To account for the additional degrees-of-freedom the timestep half compared with the linear cases. The results are much better owing to the improved resolution from the quadratic elements, even for similar degree-of-freedom counts. The fluctuations of kinetic and potential energy are significantly reduced, for the monolithic formulation they are virtually non-existent even on the coarsest mesh. As a consequence, the total energy mismatch for the monolithic formulation is also virtually absent. The monolithic formulation already has a converged energy behavior on the coarsest mesh of only 3 x 3 quadratic elements.

6.3 Verification: Mesh convergence of period

Here we focus on the prediction of the period of the wave. In practice this can be of major importance as it determines when wave groups or wave crests arrive at certain points. In the study of wave interference the exact timing of the arrival of each wave becomes critical.

Here we use the same setup as in the previous section. The simulation is performed for 10 periods. The first 2 zero crossings are discarded and the period is computed based on the following 18 zero crossings.

In figure 10 we present the convergence for the different formulations for linear finite elements (figure 10(a)) and quadratic finite elements (figure 10(b)). The reduced and monolithic formulations have virtually the same behavior. In both cases the segregated formulation performs worse. All formulations benefit from the use of quadratic over linear finite elements.

Figure 11 shows the period convergence of the monolithic formulation using different basis functions. Increasing the order and continuity of the shape functions clearly demonstrates improved behavior.

Figure 12 illustrates the benefit of the higher-order higher-order continuity NURBS basis functions over the standard linear finite elements. In both cases the horizontal discretization consists of only 6 degrees-of-freedom. For the linear finite element case this leads to a rough estimate and it results in a bad approximation of the period. In case of cubic NURBS, the spatial representation seems nearly perfect. The consequence is that the period is predicted very accurately (an error of less than $0.015\%$).

Given that the monolithic and reduced formulations are nearly indistinguishable and show superior performance, we will focus on the monolithic formulation in the following.

7 Analysis of the weak formulations

In this section we analyze the existence, uniqueness and accuracy of the solution for the monolithic weak formulation. For completeness the analysis for the reduced and staggered formulations are given in A and B, respectively. Before we perform the analysis we will shortly introduce the required mathematical tools.

7.1 Preliminaries

To facilitate the analysis of the existence, uniqueness and accuracy of the approximate solutions we first introduce some essential technical concepts in this section. Existence and uniqueness of the solutions of the weak formulations are guaranteed by the Lax-Milgram theorem. See for instance [25] or other standard works on finite element theory.

Theorem 1** (Lax-Milgram)**

Let $V$ be a Hilbert space, $V^{\prime}$ its dual space, and ${B}(\cdot,\cdot)$ a bilinear form on $V$ that is both bounded and coercive:

[TABLE]

where $|\!|\cdot|\!|$ is norm on $V$ and $C_{b},C_{c}$ are positive scalars. Then, for any $f\in V^{\prime}$, there is a unique solution $\psi\in V$ to the equation

[TABLE]

The solution $\psi\in V$ satisfies the a-priori estimate:

[TABLE]

The weak formulations can be converted into a semi-discrete formulation by straightforwardly approximating the infinite-dimensional function space by a conforming subspace. This is called the Galerkin method. Existence and uniqueness of solutions obtained with the Galerkin method are automatically inherited from the continuous formulations. The approximate error of the Galerkin solution is given by Cea’s lemma.

Lemma 2** (Cea)**

Given a bounded and coercive bilinear operator $B(\cdot,\cdot)$ and linear operator $f\in V^{\prime}$, and finite dimensional space $V^{h}$ approximating $V$. There is a unique solution $\psi^{h}\in V^{h}$ to the equation

[TABLE]

The error of the approximate solution is bounded:

[TABLE]

where $C_{b}$ and $C_{c}$ are the constants in (60).

Let us assume that we have a decomposition of the domain into a mesh ${M}=\left\{\Omega_{e}\right\}_{e=1}^{n_{el}}$ with $\Omega=\cup_{e}\Omega_{e}$. The element size is denote by $h_{e}$ with maximum $h=\max_{e}h_{e}$. The boundary $\Gamma$ is partitioned as $\Gamma=\cup_{f}\Gamma_{f}$. The boundary element size is denote by $h_{f}$ with maximum $h_{b}=\max_{f}h_{f}$. In the finite element case the decomposition represents the finite element mesh, whereas in the IGA case it represents the NURBS mesh.

Assuming sufficient regularity of the solution $\psi$, a classical convergence analysis provides the a-priori interpolation estimates:

[TABLE]

where $C_{\Omega}$, $C_{T}$ and $C_{\Gamma}$ are constants, $p$ is the minimum degree of the shape functions and $\|\cdot\|_{p,\Omega}$ and $\|\cdot\|_{p,\Gamma}$ are the norms of the $p^{th}$ derivative over $\Omega$ and $\Gamma$, respectively. Note that for simplicity we assume $h\geq h_{b}$, the validity of this assumption depends on the definition of the element and face sizes $h_{e}$ and $h_{f}$. For more details see for instance [26, 27].

7.2 Time-discrete monolithic weak formulation

Let the test function pair be \mbox{\boldmath{W}}:=(w,v)\in\mathcal{V} and the trial function pair denote $\boldsymbol{\varPhi}:=(\phi,\eta)\in\mathcal{V}$. The time-discrete weak formulation then becomes:

Find $\boldsymbol{\varPhi}^{n+1/2}\in\mathcal{V}$ such that for all \mbox{\boldmath{W}}\in\mathcal{V}:

[TABLE]

We combine the relations (46)-(47) to arrive at:

[TABLE]

Employing these relations we arrive at the time-discrete problem:

Given $\phi^{n},\eta^{n}$, find $\boldsymbol{\varPhi}^{n+1/2}\in\mathcal{V}$ such that for all \mbox{\boldmath{W}}\in\mathcal{V}:

[TABLE]

7.3 Existence and accuracy of the monolithic form

The coercivity estimate is

[TABLE]

For $\alpha=2/\Delta t$ the coercivity estimate is sharp:

[TABLE]

where the norm is defined as,

[TABLE]

Next, we consider boundedness. We write:

[TABLE]

By defining the following shorthand notation

[TABLE]

we can write,

[TABLE]

Here $A$ is a symmetric matrix:

[TABLE]

with maximum eigenvalue $\lambda_{\rm max}=2$. Using this $\lambda_{\rm max}$ and applying Cauchy-Schwarz on (7.3) we arrive at the following boundedness estimate:

[TABLE]

where we used the identities \|\mbox{\boldmath{x}}\|_{2}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mbox{\boldmath{W}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{m} and \|\mbox{\boldmath{y}}\|_{2}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\boldsymbol{\varPhi}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{m}.

Using Cea’s lemma we arrive at the following, accuracy estimate:

[TABLE]

Using the definition of the norm (70) and the interpolation estimates (65) this yields:

[TABLE]

This indicates that depending on the relative importance of the terms we could pick an additional full order of convergence order.

8 V&V of the monolithic formulation

8.1 Verification: Error convergence

Here we present the error convergence for the monolithic formulation. We measure the error in the norm (70). Note that this norm includes the timestep size $\Delta t$. Therefore we choose to perform the mesh convergence study with a fixed timestep. In this way the definition of the norm does not change when refining the mesh. The timestep is chosen as $\Delta t=T/1000$, with end time $T$, leading to negligible time stepping errors.

The convergence of the error is shown in figure 13. Note that the theoretical convergence rates are $\mathcal{O}(h^{p})$ and $\mathcal{O}(h^{p+1})$, with $p$ the order of the basis functions, for the volumetric and boundary term. This indicates that for all discretizations, except for linear finite elements, the potential boundary term is the dominant term. Remark that obtained convergence rates are larger than the theoretical values, which confirms the viability of the formulation.

8.2 Validation: Dispersion relation

To show the adequacy of the monolithic formulation we will perform a validation case which consists of predicting the phase velocity of the wave for different heights. The theoretical solution is:

[TABLE]

with phase velocity $c_{p}$. The numerical approximation is obtained on a mesh with 8 x 8 cubic NURBS elements resulting in a 10 x 11 degrees-of-freedom resolution. The time step is taken such that $80$ periods fit in the time domain.

The phase speed is computed from the wave period in the same manner as in the previous subsections. Figure 14 shows that the numerical method gives a near perfect prediction of the wave speed over a wide range of water depths. This demonstrates that the monolithic formulation produces the correct physical answer.

9 Showcase: Sloshing in closed 3D container

Here we consider a sloshing case in a three-dimensional container with a circular obstruction in the middle of the domain. The domain is represented with a 32x32x32-mesh with no penetration boundary conditions on all boundaries. Note that IGA offers the ability to represent circles exactly. A picture of the setup is presented in figure 15.

The water surface is given an initial perturbation in one quarter of the domain. Some snapshots of the resulting water surface using the monolithic formulation are shown in figure 16.

The resulting traces of kinetic, potential and total energy are depicted in figure 17. It shows that kinetic and potential energy are exchanged in an erratic fashion due to the bouncing waves between the cylinder and the outer boundaries. Just as in the 2D cases the total energy is perfectly conserved.

10 Conclusion

This paper presents a novel variational formulation to simulate linear free-surface flow. Similar to formulations found in literature the potential and water height are the primary variables. This should facilitate the extension to non-linear cases and situations with forward speed.

The novelty of the formulation is that the interior problem and kinetic and dynamic boundary conditions are combined in one formulation. This monolithic formulation results in provable energy conservation, if a proper time integrator is selected. Additional to the superior energy behavior, the monolithic formulation exhibits significantly reduced dispersion errors.

The dispersion properties of the numerical method are further improved by the adoption of Isogeometric analysis. Isogeometric analysis is a discretization approach that at the same time allows for higher-order finite elements and higher-order and higher-continuity basis functions. This later feature results in more efficient and improved behavior of the numerical solutions.

We demonstrate that the benefit of employing Isogeometric analysis over traditional finite elements is substantial. The error in the kinetic and potential energy is very small for quadratic NURBS on very coarse meshes. Moreover, on very coarse discretizations (4 to 6 degrees-of freedom per wavelength) the dispersion error is negligible.

The analysis of numerical results is performed for a two-dimensional linear wave problem. A simple sloshing case demonstrates that the results directly translate to three-dimensions. Possible further work entails the extension to nonlinear cases with forward speed.

Acknowledgements

The authors gratefully acknowledge support from Delft University of Technology. The numerical simulation are performed using the open-source library MFEM [28], while results are visualized using VisIt [29].

Appendix A Analysis of the reduced form

To analyze the properties of this variational form, we first discretize in time. Approximating $\phi$ by $\phi^{n+1/2}$ and $\phi_{tt}$ by $\phi_{tt}^{n+1/2}$ gives:

Given $\phi^{n},\phi_{t}^{n},\phi_{tt}^{n}$, find $\phi^{n+1/2}\in\mathcal{W}$ such that for all $w\in\mathcal{W}$:

[TABLE]

In order to prove existence, uniqueness and accuracy estimates, we use Lax-Milgram and Cea’s lemma. To that purpose we establish coercivity and boundedness of the weak form.

We combine the relations (46)-(47) and get:

[TABLE]

Using this relation we arrive at the following weak form:

Given $\phi_{n},\phi_{t}^{n},\phi_{tt}^{n}$, find $\phi^{n+1/2}\in\mathcal{W}$ such that for all $w\in\mathcal{W}$:

[TABLE]

By defining a problem-dependent norm as

[TABLE]

the coercivity and boundedness estimates are sharp in this norm:

[TABLE]

The coercivity follows directly from the definition of the norm, while the boundedness estimate requires the Cauchy-Schwarz inequality.

Using Cea’s lemma we arrive at the following, accuracy estimate:

[TABLE]

Using the definition of the norm (86) and the interpolation estimates (65) this yields:

[TABLE]

This indicates that depending on the relative importance of the terms we could pick an additional half order of convergence order.

Appendix B Analysis of the segregated form

B.1 Interior problem

The weak form (16) is trivially coercive:

[TABLE]

Boundedness of the weak form follows directly from Cauchy-Schwartz:

[TABLE]

in the standard $H^{1}(\Omega)$-norm.

Using Cea’s lemma the accuracy is estimate to be:

[TABLE]

Using the interpolation estimates (65) this yields:

[TABLE]

This is the standard optimal estimate for standard Poisson problems.

B.2 Boundary problem

To analyze the properties of this variational form, we first discretize in time. Substitution of the relations (67) into (18) gives:

Given $\phi^{n},\phi_{t}^{n},\eta^{n},\eta_{t}^{n}$ and $\phi_{z}$, find $\boldsymbol{\varPhi}^{n+1/2}\in\mathcal{V}_{\Gamma_{\!\!f\!s}}$ such that for all \mbox{\boldmath{W}}\in\mathcal{V}_{\Gamma_{\!\!f\!s}}:

[TABLE]

Coercivity of the boundary problem follows via:

[TABLE]

Selecting $\alpha=2/\Delta t$ yields:

[TABLE]

with the norm defined as:

[TABLE]

Next, we prove boundedness. Applying Cauchy-Schwartz we get:

[TABLE]

By defining the following shorthand notation:

[TABLE]

we can write,

[TABLE]

Here $A$ is a symmetric matrix:

[TABLE]

with maximum eigenvalue $\lambda_{\rm max}=2$. Using this $\lambda_{\rm max}$ and applying Cauchy-Schwarz on (100) we arrive at:

[TABLE]

where we used the identities \|\mbox{\boldmath{x}}\|_{2}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mbox{\boldmath{W}}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{s} and \|\mbox{\boldmath{y}}\|_{2}={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\boldsymbol{\varPhi}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{s}.

Using Cea’s lemma the accuracy is estimated as:

[TABLE]

Using the definition of the norm (96) and the interpolation estimates (65) this yields:

[TABLE]