Ocellaris: a discontinuous Galerkin finite element solver for two-phase flows with high density differences

TL;DR

Ocellaris is an open-source discontinuous Galerkin finite element solver designed for stable and accurate simulation of two-phase flows with high density differences, extending previous 2D methods to 3D with pressure coupling.

Contribution

The paper introduces a 3D extension of the Ocellaris solver with slope limiting and pressure-correction, enabling stable simulations of high-density ratio flows.

Findings

Successfully simulated dam-breaking and wave-structure interactions.

Results closely match experimental data for free-surface and pressure.

Extended the 2D slope-limiting approach to 3D with pressure coupling.

Abstract

In free-surface flows, such as breaking ocean waves, the momentum field will have a discontinuity at the interface between the two immiscible fluids, air and water, but still be smooth in most of the domain. Using a higher-order numerical method is more efficient than increasing the number of low-order computational cells in areas where the solution is smooth, but higher-order approximations cause convective instabilities at discontinuities. In Ocellaris we use slope limiting of discontinuous Galerkin solutions to stabilise finite element simulations of flows with large density jumps, which would otherwise blow up due to Gibbs oscillations resulting from approximating a factor 1000 sharp jump (air to water) by higher-order shape functions. We have previously shown a slope-limiting procedure for velocity fields that is able to stabilise 2D free-surface simulations running on a single CPU. In this paper our solver is extended to 3D and coupled to an algebraic pressure-correction scheme that retains the exact incompressibility of the direct solution used in the 2D simulations. We have tested the method on a common 3D dam-breaking test case and compared the free-surface evolution and impact pressures to experimental results. We also show how a forcing-zone approach can be used to simulate a surface-piercing vertical cylinder in an infinite wave field. In both cases the free-surface elevation and the forces and pressures compare well with published experiments. The Ocellaris solver is available as an open-source and well-documented program along with the input files needed to replicate the included results (www.ocellaris.org).

Equations

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}\right)$$

$$\nabla \cdot \mathbf{u}$$

$$\frac{\partial \rho }{\partial t}+\mathbf{u}\cdot \nabla \rho$$

$$\rho =\alpha {\rho }_{\text{water}}+(1-\alpha ){\rho }_{\text{air}},$$

$$\nu =\alpha {\nu }_{\text{water}}+(1-\alpha ){\nu }_{\text{air}}.$$

$$\mathcal{A}(\mathbf{u},\mathbf{v};\mathbf{w})+\mathcal{B}(p,\mathbf{v})$$

$$\mathcal{C}(\mathbf{u},q)$$

$$\mathcal{F}(\alpha ,\beta ;\mathbf{w})$$

$$\mathbf{u}$$

$$\frac{\partial \mathbf{u}}{\partial \mathbf{n}}$$

$$\left\{\left\{u\right\}\right\}$$

$$\displaystyle\left\llbracket\,u\,\right\rrbracket$$

$$\displaystyle\left\llbracket\,{\bm{{u}}}\,\right\rrbracket_{\bm{{n}}}$$

$${\int }_{\mathcal{T}}\frac{\rho }{\Delta t}({\gamma }_1\mathbf{u}+{\gamma }_2{\mathbf{u}}^n+{\gamma }_3{\mathbf{u}}^{n-1})\mathbf{v}\mathrm{d}\mathbf{x}$$

$$\displaystyle\qquad-\ \int_{\mathcal{T}}{\bm{{u}}}\cdot\nabla\cdot(\rho{\bm{{v}}}\otimes{\bm{{w}}})\,\mathrm{d}{\bm{{x}}}\ +\ \int_{\mathcal{S}}{\bm{{w}}}\cdot{\bm{{n}}}^{+}\,\hat{\bm{{u}}}^{\bm{{w}}}\cdot\left\llbracket\,\rho{\bm{{v}}}\,\right\rrbracket\,\mathrm{d}s$$

$$\displaystyle\qquad+\ \int_{\mathcal{T}}\mu\left(\nabla{\bm{{u}}}+(\nabla{\bm{{u}}})^{T}\right):\nabla{\bm{{v}}}\,\mathrm{d}{\bm{{x}}}\ +\ \int_{\mathcal{S}_{I}}\kappa_{\mu}\left\llbracket\,{\bm{{u}}}\,\right\rrbracket\cdot\left\llbracket\,{\bm{{v}}}\,\right\rrbracket\,\mathrm{d}s$$

$$\displaystyle\qquad-\ \int_{\mathcal{S}}(\left\{\!\left\{\mu\left(\nabla{\bm{{u}}}+(\nabla{\bm{{u}}})^{T}\right)\right\}\!\right\}\cdot{\bm{{n}}}^{+})\cdot\left\llbracket\,{\bm{{v}}}\,\right\rrbracket\,\mathrm{d}s$$

$$\displaystyle\qquad-\ \int_{\mathcal{S}_{I}}(\left\{\!\left\{\mu\left(\nabla{\bm{{v}}}+(\nabla{\bm{{v}}})^{T}\right)\right\}\!\right\}\cdot{\bm{{n}}}^{+})\cdot\left\llbracket\,{\bm{{u}}}\,\right\rrbracket\,\mathrm{d}s$$

$$\displaystyle\qquad-\ \int_{\mathcal{T}}p\,\nabla\cdot{\bm{{v}}}\,\mathrm{d}{\bm{{x}}}\ +\ \int_{\mathcal{S}}\hat{p}\,{\bm{{n}}}^{+}\cdot\left\llbracket\,{\bm{{v}}}\,\right\rrbracket\,\mathrm{d}s\ =\ \int_{\mathcal{T}}\rho\,{\bm{{g}}}\,\mathrm{d}{\bm{{x}}}.$$

$$\displaystyle\int_{\mathcal{S}}\hat{{\bm{{u}}}}^{p}\cdot{\bm{{n}}}^{+}\left\llbracket\,q\,\right\rrbracket\,\mathrm{d}s\ -\ \int_{\mathcal{T}}{\bm{{u}}}\cdot\nabla q\,\mathrm{d}{\bm{{x}}}\ =\ 0.$$

$$\displaystyle\int_{\mathcal{T}}\frac{1}{\Delta t}(\gamma_{1}\alpha^{n+1}+\gamma_{2}\alpha^{n}+\gamma_{3}\alpha^{n-1})\beta\,\mathrm{d}{\bm{{x}}}+\int_{\mathcal{S}}\hat{\alpha}^{n+1}\,{\bm{{w}}}\cdot{\bm{{n}}}^{+}\left\llbracket\,\alpha\,\right\rrbracket\,\mathrm{d}s=0,$$

$${\mathbf{w}}^{\ast }={\mathbf{w}}^{\ast ,n+1}=2{\mathbf{w}}^n-{\mathbf{w}}^{n-1},$$

$$\left[\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ \mathbf{C}& 0\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \mathbf{p}\end{array}\right]=\left[\begin{array}{c}\mathbf{d}\\ \mathbf{e}\end{array}\right],$$

$$\mathbf{A}{\mathbf{u}}^{\ast }=\mathbf{d}-\mathbf{B}{\mathbf{p}}^{\ast }.$$

$$\mathbf{M}(\mathbf{u}-{\mathbf{u}}^{\ast })=-\mathbf{B}(\mathbf{p}-{\mathbf{p}}^{\ast }).$$

$$\mathbf{C}{\mathbf{M}}^{-1}\mathbf{B}\mathbf{p}=\mathbf{C}{\mathbf{M}}^{-1}\mathbf{B}{\mathbf{p}}^{\ast }-\mathbf{e}+\mathbf{C}{\mathbf{u}}^{\ast }.$$

$$\mathbf{u}=-{\mathbf{M}}^{-1}\mathbf{B}(\mathbf{p}-{\mathbf{p}}^{\ast }).$$

$$\mathbf{C}\mathbf{u}-\mathbf{e}$$

$$\displaystyle\overset{\mathclap{\text{\eqref{eq:IPCSA_spliterr}}}}{=}\ \ {\mathbf{{C}}}{\mathbf{{M}}}^{-1}{\mathbf{{B}}}({\mathbf{{p}}}-{\mathbf{{p}}}^{*})+{\mathbf{{C}}}{\mathbf{{u}}}^{*}-{\mathbf{{e}}}$$

$$\displaystyle\overset{\mathclap{\text{\eqref{eq:IPCSA_eqp}}}}{=}\ \ 0.$$

$${\int }_{\mathcal{T}}{\kappa }_{\text{fz}}{\beta }_{\text{fz}}(\mathbf{u}-{\mathbf{u}}_D)\cdot \mathbf{v}\mathrm{d}\mathbf{x},$$

$${\Psi }_w(x,z,t)=B_{0}z+\sum _{j=1}^N{B}_j\frac{\sinh jkz}{\cosh jkD}\cos jkx,$$

$$\Psi$$

$$Z$$

$${\mathbf{u}}_x$$

$${\mathbf{u}}_z$$

$$\mathrm{Co}$$

$${\mathrm{Co}}_f$$

Full text

We have previously shown a slope-limiting procedure for velocity fields that is able to stabilise 2D free-surface simulations running on a single CPU. In this paper our solver is extended to 3D and coupled to an algebraic pressure-correction scheme that retains the exact incompressibility of the direct solution used in the 2D simulations. We have tested the method on a common 3D dam-breaking test case and compared the free-surface evolution and impact pressures to experimental results. We also show how a forcing-zone approach can be used to simulate a surface-piercing vertical cylinder in an infinite wave field. In both cases the free-surface elevation and the forces and pressures compare well with published experiments. The Ocellaris solver is available as an open-source and well-documented program along with the input files needed to replicate the included results (www.ocellaris.org).

1 Introduction

Increasing the polynomial approximation order is more efficient than increasing the number of computational cells when solving partial differential equations (PDEs) where the solution is smooth (Babuška and Dorr, 1981). Performing more work locally in higher-order cells is also advantageous for today’s parallel computers (Kubatko et al., 2009; Kirby et al., 2012; Huerta et al., 2013). However, for equations with convective operators, a jump in the solution or the coefficients will cause nonlinear convective instabilities—spurious Gibbs oscillations at discontinuities—that will eventually destroy the solution if the ratio between the high and low sides of the jump is large. Linear numerical schemes, such as the discontinuous Galerkin finite element method (DG FEM), produces linear algebraic equations from discretisation of linear PDEs. Such linear schemes cannot guarantee convective stability if they are not monotonic, and a linear monotonic scheme must be first-order and is hence highly dissipative (Krivodonova et al., 2004; Harten, 1983). One way to get around the problem of having to chose between high-order approximations and convective stability is to make the scheme nonlinear by the inclusion of a nonlinear cell-wise projection operator, a slope limiter.

We have previously shown that component-wise slope limiters can be used on the convected velocity field to stabilise two-phase flow simulations with large density jumps (Landet et al., 2018). Most existing methods used for the simulation of air/water two-phase flows employ low-order finite volume methods for the discretisation of the governing PDEs. In these methods the solution is piecewise constant and the convective instabilities can be dealt with by using a flux limiter, a non-linear facet-local parameter that blends the upwind and downwind fluxes to obtain stability without excessive diffusion (Hirt and Nichols, 1981; Ubbink, 1997; Weller et al., 1998; Popinet, 2003; Kleefsman et al., 2005).

Ocellaris solves the variable-density Navier–Stokes equations for two-phase flow (section 2). The numerical method is based on a higher-order interior-penalty DG FEM (section 3) with a component-wise velocity slope limiter for the convected velocity field (section 3.4) and an algebraic pressure-correction method for the pressure–velocity coupling (section 4). A regular wave model based on a stream-function formulation is used for initial and boundary values for water-wave simulations. The details of this model are presented in section 5, including how the stream-function approach can be used to specify one velocity field, valid in both the air and water phases, which is divergence free, satisfies all boundary conditions, and conforms to the free surface. Section 5 also describes how a forcing-zone approach is used for damping free-surface disturbances near the boundaries in order to avoid unwanted reflections. The implementation of Ocellaris is briefly described in section 6, and the main steps required to set up an Ocellaris simulation are described in section 6.1.

The BlendedAlgebraicVOF multiphase model used in this paper computes a piecewise constant density distribution based on an algebraic VOF method (Muzaferija et al., 1998). Other multiphase models are available, and it is also possible to define a custom model in the Ocellaris simulation input file (YAML format, see Landet (2019b) for details). Using such a very simple model is not optimal, but the velocity slope limiter will lower the effective order of the obtained velocities at the density jump, so perhaps not much is lost in terms of obtainable accuracy from using a relatively simple free-surface capturing method. An interface-capturing method that can take full advantage of the high order of the convecting velocity field is an interesting research topic, and would most likely require fewer cells in the free-surface region than what is used here. Still, the results show good agreement with lab experiments, and away from the free surface the overall method retains the high-order approximation properties of the Navier–Stokes discretisation since the true density field is constant here.

Section 7.1 shows the performance of the Ocellaris solver on a 3D dam-breaking test case. The test is meant to simulate a “green-water” event, a large wave breaking over a ship deck and impacting the cargo, in this case a single container outfitted with pressure gauges. The results show that both the free-surface evolution and the impact pressures compare well with published experiments. The second test case, presented in section 7.2, is a vertical cylinder exposed to steep regular waves. This is meant to simulate wave loads on an offshore windmill or another type of slender surface-piercing marine structure. This test shows good agreement in regard to the total force on the cylinder when compared to laboratory experiments. The test also shows that the forcing-zone approach presented in section 5 works well in a DG FEM setting to dampen free-surface disturbances near the inlet and outlet of the numerical wave tank. Finally, the discussion in section 8 concludes that Ocellaris’ high-order multiphase flow solver can successfully simulate complex 3D air/water free-surface flows at high Reynolds numbers.

2 Mathematical model of free-surface flow

Ocellaris solves the variable-density Navier–Stokes equations 1, 2 and 3 with piecewise constant density and viscosity. Standard notation is used for the unknown functions; ${\bm{{u}}}$ is the velocity, $p$ is the pressure, and $\rho$ is the fluid density. The coefficients are the dynamic viscosity $\mu$ and the acceleration of gravity $g$,

[TABLE]

Ocellaris is currently designed to study air/water free-surface flows at high Reynolds numbers, and does not include a turbulence closure model or the effect of surface tension. The error made by neglecting these effects is small for the included benchmark tests (Kleefsman et al., 2005; Paulsen et al., 2014). The program design is made to be easily extendable and there should be no fundamental problem with including such effects at a later time. For the first public release, the sole focus of the Ocellaris project has been to show that using higher-order DG methods for free-surface flows is feasible and can be made stable in regard to convective instabilities without compromising mass conservation.

A VOF colour function approach is used for density transport and free-surface capturing (Hirt and Nichols, 1981),

[TABLE]

where $\alpha\in[0,1]$ is the colour/indicator function which is used as the unknown instead of the fluid density. The kinematic viscosity, $\nu=\nicefrac{{\mu}}{{\rho}}$ is computed similarly,

[TABLE]

3 Discontinuous Galerkin discretisation

The numerical method is an extension of Landet et al. (2018) which builds on Cockburn et al. (2005) and uses the symmetric interior-penalty (SIP) method for the viscous term (Arnold, 1982). The domain is discretised as an irregular mesh comprised of tetrahedral cells. Let $\mathcal{T}$ be the set of all cells and $\mathcal{S}$ the set of all facets in the mesh. Polynomial function spaces of degree $k$ on each cell $K$ are denoted $P_{k}(K)$. These have no continuity at cell boundaries and no inherent boundary conditions.

We use calligraphic typeface to denote operators and sets, bold italic for vectors functions and italic for scalar functions. For nabla the conventions $(\nabla{\bm{{u}}})_{ij}=\partial_{j}{\bm{{u}}}_{i}$ and $(\nabla\cdot{\bm{{\sigma}}})_{i}=\partial_{j}\sigma_{ij}$ are used. For time derivatives ${\bm{{u}}}={\bm{{u}}}^{n+1}$ is the unknown trial function and ${\bm{{u}}}^{n}$ and ${\bm{{u}}}^{n-1}$ are the known values at the previous two time steps. A second-order backwards-differencing formulation, BDF2, is used for time integration. The parameters are $\{\gamma_{1},\gamma_{2},\gamma_{3}\}=\{\nicefrac{{3}}{{2}},-2,\nicefrac{{1}}{{2}}\}$, see the first integral in equation 14. The convecting velocity ${\bm{{w}}}$ is considered known—which linearises the momentum equation—and ${\bm{{w}}}$ is Hdiv-conforming such that the flux is continuous, \left\llbracket\,{\bm{{w}}}\,\right\rrbracket_{\bm{{n}}}=0, see section 4 for details.

The governing equations, 1, 2 and 3, are cast into the following form: find ${\bm{{u}}}\in[P_{2}(K)]^{3}$, $p\in P_{1}(K)$, and $\alpha\in P_{0}(K)$ such that

[TABLE]

in the tessellated domain $\mathcal{T}$ subject to

[TABLE]

The discontinuous Galerkin (DG) method works by breaking integrals over the whole domain into a sum of integrals over each mesh cell $K\in\mathcal{T}$, and defining fluxes of the unknown functions between these cells. The average and jump operators across an internal facet between two cells $K^{+}$ and $K^{-}$ are defined as

[TABLE]

where $u^{+}$ is the value of $u$ along the internal facet when computed using the shape functions and degrees of freedom related to the $K^{+}$ cell and vice versa for $u^{-}$. For the exterior facets, which only have one connected cell due to being located on the domain boundary, let the connected cell be denoted $K^{+}$ such that {\bm{{n}}}^{+}\!\cdot\left\llbracket\,{\bm{{u}}}\,\right\rrbracket={\bm{{n}}}^{+}\!\cdot{\bm{{u}}}^{+}={\bm{{n}}}\cdot{\bm{{u}}}. Take $\left\{\!\left\{u\right\}\!\right\}=u$ and otherwise let all $K^{-}$ values related to the non existing element outside the domain be zero.

3.1 Momentum equation

The momentum equation 1 is discretised using the SIP method for the elliptic term (Arnold, 1982) and otherwise using the fluxes from Cockburn et al. (2005). The application of boundary conditions and the choice of the penalty parameter $\kappa_{\mu}$ is described in detail in Landet et al. (2018). The flux of pressure is $\hat{p}=\left\{\!\left\{p\right\}\!\right\}$ and the convective flux $\hat{\bm{{u}}}^{\bm{{w}}}$ is a pure upwind flux. After multiplication with ${\bm{{v}}}$ and integration over $\mathcal{T}$, followed by integration by parts and application of the SIP method to the viscosity, ${\mathcal{{A}}}$ can be found as the bilinear part containing ${\bm{{u}}}={\bm{{u}}}^{n+1}$, ${\mathcal{{B}}}$ as the bilinear part containing $p$, and ${\mathcal{{D}}}$ as the linear part of

[TABLE]

On Dirichlet boundaries $\hat{\bm{{u}}}^{\bm{{w}}}$ is the upwind value. Depending on flow direction this is either ${\bm{{u}}}$ or ${\bm{{u}}}_{D}$. On Neumann boundaries let $\hat{\bm{{u}}}^{\bm{{w}}}={\bm{{u}}}$. The boundary conditions (BCs) can be set separately for each velocity component, allowing slip, ${\bm{{u}}}\cdot{\bm{{n}}}=0$, or non-slip, ${\bm{{u}}}=0$, BCs. On all boundaries take $\hat{p}=p$.

3.2 Continuity equation

The equation used to ensure a divergence-free velocity field (2) is multiplied by $q$ and integrated over $\mathcal{T}$. The flux is $\hat{{\bm{{u}}}}^{p}=\left\{\!\left\{{\bm{{u}}}\right\}\!\right\}$. After integration by parts ${\mathcal{{C}}}({\bm{{u}}},q)$ and ${\mathcal{{E}}}(q)$ can be found from

[TABLE]

The non-zero ${\mathcal{{E}}}(q)$ results from using the boundary conditions in the flux, $\hat{\bm{{u}}}^{p}={\bm{{u}}}_{D}$ on Dirichlet boundaries. On Neumann boundaries the unknown function is used directly, $\hat{\bm{{u}}}^{p}={\bm{{u}}}$, but when assembling the contribution from each velocity component separately it is beneficial to set $\hat{\bm{{u}}}^{p}=0$ on free-slip surfaces also for components with Neumann BCs. Otherwise, it is possible to get ${\bm{{u}}}\cdot{\bm{{n}}}\neq 0$ due to small errors in the generated meshes when the boundaries are not perfectly approximated by simplices.

3.3 Density transport

The HRIC method (Muzaferija et al., 1998) is used to define the colour function flux $\hat{\alpha}$. The HRIC flux limiter is stable and avoids the excessive interface diffusion of a standard upwind flux. It would be ideal to construct a better interface capturing method that can take advantage of the fact that the velocity is in $[P_{2}(K)]^{3}$, but in this first public release of Ocellaris a standard algebraic VOF method is used to validate the stability and applicability of the overall free-surface flow solver.

To construct the DG operators for the density transport equation 8, the strong form in equation 3 is multiplied by the test function $\beta$, and the result is integrated over the domain. After integration by parts and discarding derivatives of the piecewise constant functions, ${\mathcal{{F}}}$ can be found as the bilinear part and ${\mathcal{{G}}}$ as the linear part of

[TABLE]

3.4 Velocity slope limiter

Using functions spaces that are higher order than piecewise constants to approximate the velocity will cause Gibbs oscillation instabilities when the density field has jumps, such as the order $1000$ jump in density in free-surface simulations of air and water. Slope-limiting techniques can be used to remove such instabilities. The velocity slope-limiting approach taken here is the component-wise hierarchical Taylor-polynomial based slope limiter described in Landet et al. (2018), where each velocity component of the convected velocity, ${\bm{{u}}}$, is slope limited by a vertex-based scalar slope limiter (Kuzmin, 2010) while the convecting velocity, ${\bm{{w}}}$, is left unchanged. This approach requires no tuning and the time spent in the velocity slope limiter is only around $0.2\text{\,}\mathrm{\char 37\relax}$ of the total running time.

The effect of using a slope limiter is dramatic. Figure 1 shows the evolution of the total, potential and kinetic energy in the 3D dam-breaking test case presented more thoroughly in section 7.1. Water rushes out of the broken dam towards a small box-shaped object which is impacted after about $0.4\text{\,}\mathrm{s}$. As can be seen, the slope-limited simulation preserves total energy well, trading potential for kinetic energy up to the time of impact. After the impact the conservation of energy is not perfect,but the kinetic energy is always controlled and does not blow up. The second plot shows that the kinetic energy quickly blows up when using the same simulation setup with the velocity slope limiter deactivated. The simulation is automatically stopped after less than $50$ time steps when the Courant number passes $1000$.

4 Solution algorithm

The time-stepping procedure is shown in algorithm 1. To decouple the density field solver from the Navier–Stokes solver, and to linearise the momentum equation, a second-order extrapolation ${\bm{{w}}}^{*}$ is used to estimate the convective velocity ${\bm{{w}}}$,

[TABLE]

which means $\rho^{n+1}$ can be computed before the velocity at time $t=(n+1)\Delta t$.

The discretised matrix version of equations 6 and 7,

[TABLE]

is a saddle-point problem, which is solved by the incremental pressure-correction scheme on algebraic form (IPCS-A). This is an incremental version of the classic Chorin-Temam methods (Chorin, 1968; Temam, 1969). ${\mathbf{{A}}}$, ${\mathbf{{B}}}$ and ${\mathbf{{C}}}$ are sparse matrix versions of the ${\mathcal{{A}}}$, ${\mathcal{{B}}}$ and ${\mathcal{{C}}}$ operators, discretised using the discontinuous Galerkin method described in section 3. The unknowns ${\bm{{u}}}$ and $p$ are now vectors of degrees of freedom, ${\mathbf{{u}}}$ and ${\mathbf{{p}}}$.

The first IPCS-A step is momentum prediction, performed by using an approximate pressure field ${\mathbf{{p}}}^{*}$ inserted into the first row of equation 18,

[TABLE]

A splitting error is now introduced. First let ${\mathbf{{A}}}={\mathbf{{M}}}+{\mathbf{{R}}}$ where ${\mathbf{{M}}}$ is a scaled mass matrix resulting from assembly of the time derivative and ${\mathbf{{R}}}$ contains the convective and diffusive operators. Then make the assumption that ${\mathbf{{R}}}({\mathbf{{u}}}-{\mathbf{{u}}}^{*})\approx 0$ and use this when subtracting equation 19 from the first row of equation 18,

[TABLE]

The matrix ${\mathbf{{M}}}$ is block diagonal, and is hence cheap to invert. Use this property and the divergence-free criterion, ${\mathbf{{C}}}{\mathbf{{u}}}={\mathbf{{e}}}$, to remove the unknown ${\mathbf{{u}}}$ from equation 20, and reorganise this into an equation for ${\mathbf{{p}}}$,

[TABLE]

The velocity ${\mathbf{{u}}}$ can be recovered without solving a linear system, simply by substituting the pressure ${\mathbf{{p}}}$ from the solution of equation 21 into equation 20,

[TABLE]

The momentum-prediction equation 19 is solved followed by the pressure-correction equation 21 and the velocity update equation 22 in an iterative manner until the required accuracy is reached. The resulting velocity field is projected into a BDM-type function space (Brezzi et al., 1987) where it becomes exactly divergence free—the velocity flux is continuous across internal facets and the velocity field is pointwise divergence free inside each cell, see Landet and Mortensen (2019) for details. The result from this projection is stored both as ${\bm{{u}}}^{n+1}$ and ${\bm{{w}}}^{n+1}$. The last step is slope limiting ${\bm{{u}}}^{n+1}$ such that any instabilities are prevented from growing through the time derivative.

Solving the coupled problem in equation 18 without some form of pressure and velocity splitting requires a direct solver, and such solvers scale badly in terms of parallel computing efficiency. The IPCS-A method can be solved using standard iterative Krylov methods which scale much better. Exact mass conservation, ${\mathbf{{C}}}{\mathbf{{u}}}-{\mathbf{{e}}}=0$, is still ensured on the algebraic level in each iteration,

[TABLE]

For the example simulations described in section 7, the maximum cell wise divergence error in the convecting velocity range from $10^{-5}$ to $10^{-3}$, while the convected velocity (which is slope limited) has a maximum cell wise divergence error ranging from $10^{2}$ to $10^{4}$. This is with $10^{-8}$ relative error and $10^{-10}$ absolute error convergence criteria in the pressure-correction Krylov solver. The result is that the simulation in section 7.1, a dam breaking in a closed box, has less than $3\text{\times}{10}^{-5}\text{\,}\mathrm{\char 37\relax}$ change in total mass from time step 3 to the final time step (3415 steps, two seconds simulated time).

5 Incoming waves and boundary reflections

In our second numerical example, shown in detail in section 7.2, we will study a surface-piercing vertical cylinder in an infinite wave field. A truncated computational domain is inevitably required to compute the solution in a finite amount of time. Our solution to the problem of reflected waves from the boundaries is to use a forcing-zone approach to wave damping (Perić and Abdel-Maksoud, 2016, 2018). We impose Dirichlet boundary conditions to the velocity and density fields at both inlet and outlet boundaries. These boundaries are then padded by forcing zones inside the domain which penalise deviations from the undisturbed wave field. This damps out any free-surface disturbances caused by the structure without having to damp out the incident wave field itself. See Perić and Abdel-Maksoud (2018) for estimates of the minimum forcing zone size and the penalty magnitude needed to obtain a given reduction in reflected wave amplitudes.

Figure 2 shows the forcing zones used in the test case described in section 7.2. The inlet zone is $0.75\text{\,}\mathrm{m}$ long and the outlet zone is $1.0\text{\,}\mathrm{m}$ long. The shape of the zone is a quadratic polynomial with maximum value $1.0$ at the boundary and zero first-derivative towards the inner domain. To force the solution towards the incident wave field, a penalty term is added to the momentum equations,

[TABLE]

where $\kappa_{\text{fz}}$ is the penalty parameter, in our calculations $10$, $\beta_{\text{fz}}$ is the zone shape shown in figure 2 and ${\bm{{u}}}_{D}(x,z,t)$ is the incident-wave velocity field used for initial and boundary conditions. The same forcing-zone approach is added to the density-transport equation for the colour function with the same forcing-zone shapes and the same penalty parameter, $\kappa_{\text{fz}}=10$.

The initial and boundary conditions for the velocity and density fields are computed from Fenton stream function wave theory. This is a high-order regular wave theory based on approximating a stream function by a truncated Fourier series. This method of constructing non-linear regular waves was pioneered by Dean (1965). Our implementation is based on Rienecker and Fenton (1981), which is often referred to as “Fenton” stream function wave theory to differentiate it from the original “Dean” stream function wave theory.

The Fenton method is based on collocation in $N+1$ points on the free surface along half the wave length, $\lambda$. The Fenton stream function,

[TABLE]

is non-linear in the wave height $\eta$ since $z=\eta$ in the collocation points. $\Psi_{w}$ does a-priori satisfy the bottom boundary condition at $z=0$ and also the Laplace equation $\nabla^{2}\Psi_{w}=0$. The following conditions are imposed in the Newton–Raphson iterations that are applied to compute the unknown coefficients in equation 25: (i) the free surface is a stream line, $\Psi_{w}=$ const., (ii) the pressure is constant at the free surface, (iii) the wave height is $H$, such that $\eta(0)-\eta(\lambda/2)=H$, (iv) the mean wave elevation is $D$, such that $\int_{0}^{\lambda/2}\eta\,\mathrm{d}x=D\lambda/2$. See Rienecker and Fenton (1981) for details.

We use the same approach to find a compatible stream function for the air phase, $\Psi_{a}$, when $\Psi_{w}$ has been determined. The stream function in equation 25 is now linear in the unknowns since $\eta$ is known, so the expansion coefficients can be found by a single linear solve. This gives two separate domains where the boundary conditions are satisfied and the divergence is zero, but the velocity parallel to the free surface has a discontinuity at the free surface since the method is based on potential theory and does not contain the viscosity that causes a velocity shear at the free surface. To approximate this, a blended stream function is constructed,

[TABLE]

where the blending function $f(Z)$ is zero in the water and unity in the air above the blending zone. All blending is performed in the air phase, and the blending zone has height $d$, starting from the free surface. A fifth-order polynomial smooth step function is used for $f(Z)$. This function has zero first and second derivatives at the top and bottom of the blending zone. The resulting velocity field,

[TABLE]

is continuous and satisfies the continuity equation everywhere, the dynamic and kinematic free-surface boundary conditions at the interface, and the free-slip boundary conditions at the top and bottom of the domain. These properties make the blended solution a good choice to use as initial and boundary conditions in an exactly divergence-free Navier–Stokes solver.

6 Implementation

The Ocellaris solver (Landet, 2019a) is built on FEniCS (Logg et al., 2012) and PETSc (Davis, 2004; Balay et al., 1997, 2018; Dalcin et al., 2011; LLNL, ). The overwhelming majority of the code, including the definition of the weak form and the time-stepping procedure, is implemented in Python 3. The FEniCS form compiler, FFC, is used to compile the weak form, defined in UFL Python format, to optimised C++ code Kirby and Logg (2006); Ølgaard et al. (2008); Alnæs et al. (2014). Wherever tight loops over the mesh cells or facets are needed in other parts of the program, the loop is written in C++. For the simulation examples shown in section 7, the additional C++ code is restricted to parts of the VOF implementation and the velocity slope limiter.

Ocellaris is developed using automated unit and MMS (method of manufactured solutions) testing. Unfortunately, MMS testing of two-phase flows with discontinuous density fields is not well developed, so testing of this functionality involves running full time simulations, which is not done automatically on each change of the code as this would be too expensive. The Ocellaris program design is made up of independent pluggable components, letting the user define which combination of pressure-splitting scheme, free-surface model, slope limiter and more to use, simply by referencing the relevant components in the input file. This also means that there are automated MMS test of the implemented pressure-correction solvers, but they are tested with a single-phase flow model where classical analytical solutions are available, such as the Taylor–Green vortex (Green and Taylor, 1937).

6.1 Example input file

The following file listings show excerpts from the input file used to simulate regular waves passing a vertical, surface-piercing cylinder. More information about this example can be found in section 7.2. Complete documentation of the input file format can be found in (Landet, 2019b), and the full input file can be found in (Landet, 2019c). All Ocellaris input files are written on YAML format and start with a mandatory header section followed by an optional metadata section. An example of these two input file sections is shown in LABEL:lst:cyl_meta.

Ocellaris supports defining constants to be used throughout the input file. Defining these once on top of the file makes parameter studies and input changes easier to perform and makes the file easier to read. The definition of convenience constants and physical constants can be seen in LABEL:lst:cyl_constants.

The mesh section is shown in LABEL:lst:cyl_msh. The mesh file, created in gmsh (Geuzaine and Remacle, 2009), is loaded using the meshio Python package which implements readers and writers for many unstructured mesh file formats.

Known fields functions are defined in LABEL:lst:cyl_fields_wave and LABEL:lst:cyl_fields_zone. Known fields are used in Ocellaris to define initial and boundary conditions and also to define the location of the free-surface wave damping zones described in section 5. The waves are defined using the raschii Python package which produces C++ code for the Fenton and Stokes regular wave models for use in FEniCS-based solvers. Note also that Ocellaris accepts Python expressions such as "py$ H - d" instead of scalars, booleans and strings. These expressions can be used to define parameters in terms of the user-defined constants.

The damping zone in LABEL:lst:cyl_fields_zone is implemented as a generic scalar field. The C++ code is compiled inside a namespace that includes all user-specified constants and an array, x, defining the coordinates where the field is to be evaluated. Using C++ lambdas allows using multi-line expressions to compute the field value. Standard C++ expressions, such as "x[0] + x[1]*x[2]" ($x+yz$), can also be given for fields which do not need multiple statements to compute the field value. The inlet damping zone is defined equivalently to the outlet damping zone, but is not shown here for brevity.

A definition of a forcing zone is shown in LABEL:lst:cyl_fz. There are four such zones in the simulation, damping the momentum and density fields at the inlet and outlet, see section 5. The zone in the example shows the damping of the momentum equations at the outlet boundary. It uses the previously defined known fields waves and outlet zone to specify the target value and field location.

The initial conditions are defined in LABEL:lst:cyl_ic. The naming scheme uses the postfix p to specify the value of a field at $t=0$, the previous time-step value, and 0 and 2 to specify the $x$ and $z$ directions respectively. To use higher-order time stepping from the start, the values at $t=-\Delta t$ can be specified by using the pp prefix, but that is normally only done in convergence tests where the analytical solution is known.

LABEL:lst:cyl_bc shows the definition of boundary conditions for the inlet. The inside code is used to select the facets on the inlet. If a boundary region is marked with an integer identifier in the mesh generator, then facet selection can be based on this identifier instead. For the cylinder-in-waves example the boundary facets are all marked with C++ code as shown in LABEL:lst:cyl_bc. The on_boundary boolean flag is true for external facing facets.

The Navier–Stokes solver section specifies which velocity–pressure splitting scheme to use, the iteration tolerances and the PETSc Krylov solver parameters. Any PETSc configuration variable can be changed for maximum flexibility. In the example shown in LABEL:lst:cyl_solver, the Ocellaris default options for the momentum and pressure PETSc solvers are used and only the convergence criteria and the number of inner iterations (the number of pressure corrections per time step) are changed. After ten time steps, only two pressure corrections are performed per time step. The Krylov solver tolerances are given as three numbers; the tolerance for the first three pressure corrections, the value for the mid range of pressure corrections, and finally the values for the last five pressure corrections. Since, after the first ten time steps, there are only two pressure corrections per time step, only the last value in the lists matter. This gradual decrease of tolerances is done to avoid spending a lot of time in the Krylov solver in the beginning of the simulation when the pressure corrections are not converged and exact answers are not needed.

The multiphase VOF input sections are shown in LABEL:lst:cyl_mps. The HRIC VOF scheme is selected and 5 sub-cycles—the number of advection steps of the VOF colour function per time step of the Navier–Stokes solver—are applied to maximise sharpness of the interface by lowering the effective Courant number in the VOF solver.

The final excerpt from the input file is shown in LABEL:lst:cyl_slopelim. Here the velocity slope limiter is configured to use the scalar HierarchicalTaylor slope-limiting method by Kuzmin (2010) for each component of the convected velocity field.

7 Numerical examples

The following example simulations have been run in Ocellaris. Source code and input files, which can be used to reproduce all the results, can be found in (Landet, 2019c). Visualisations and isosurfaces have been made in Paraview (Ahrens et al., 2005). The meshes, all available for download, have been generated and optimised in Gmsh (Geuzaine and Remacle, 2009).

7.1 3D dam breaking

Kleefsman et al. (2005) presents experimental and numerical results for a 3D dam-breaking test case where a small container is subjected to a fast-moving front of water that splashes over and around the container. The tank is $3\text{\times}1\text{\times}1\text{\,}\mathrm{m}$ and the container is $16\text{\times}40\text{\times}16\text{\,}\mathrm{cm}$. There are three surface-height probes that are initially dry and one that is in the fluid (H4). The container is fitted with eight pressure gauges close to the centre line. The locations of the gauges can be found in Issa and Violeau (2006) along with an exact description of the tank and container geometries.

The Ocellaris simulations have been run based on irregular tetrahedral meshes. The same mesh input has been used with three different mesh densities, resulting in a coarse mesh with a total of $18\,676$ tetrahedra, a medium-density mesh with $40\,450$ tetrahedra and a fine mesh with $123\,628$ tetrahedra. A cross-sectional view of the medium-density mesh can be seen in figure 4. The time step is adaptively controlled based on the maximum cell-based (28) and facet-based (29) Courant numbers. The time step is halved if one of the two is above $0.3$ and doubled if they are both below $0.05$. The Courant numbers are computed for each cell and facet as

[TABLE]

where $D_{c}$ is the cell diameter, $S_{f}$ is the facet area, and $V_{c}$ is the cell volume. The facet-based Courant number is used in the HRIC transport scheme for the colour function and this is sub-cycled with $N_{\text{sc}}$ sub-cycles per time step. We have used $N_{\text{sc}}=5$ which makes the cell- and facet-based Courant numbers similar in magnitude.

Figure 5 shows the evolution of the four free-surface probes for the medium and fine mesh. The fine mesh fits best with the experiments, but both mesh resolutions show good comparison for probes H2–H4. The H1 probe, where the free-surface height is multi valued, does not compare well. It is not clear exactly what heigh is measured in the experiments at this location, we have reported the topmost surface intersection. A visualisation of the colour function at $t=$0.4075\text{,}\mathrm{s}, when pressure probe P2 spikes in the experimental results, is shown in [figure 7](#S7.F7). The colour function isosurface from $t=$1.1075\text{\,}\mathrm{s}, when the free surface is multi-valued behind the container, can be seen in figure 7.

The pressure probe time series are shown in figure 8. The pressure probes on top of the container, P5–P8, all show very similar behaviour, we have included the results for probe P7 as an example, the other plots would have been roughly identical. On the side of the container facing the wave, the mesh convergence can be seen clearly in the P1 and P2 probes. Further up the container wall the pressure peak is less pronounced in the numerical results, possibly due to the interface being smeared over more cells than optimal, creating a more gradual rise in density at the pressure sensors than what happened in the experiment.

7.2 Cylinder in regular waves

Grue and Huseby (2002) gives experimental results for the inline force on a cylinder in regular waves. Figure 3d in that reference is used by Paulsen et al. (2014) as a benchmark problem, and in this numerical example we will do the same. Fenton stream-function waves are applied in a numerical wave tank which is $4\text{\,}\mathrm{m}$ long and $0.5\text{\,}\mathrm{m}$ wide. Our numerical domain height is $1.2\text{\,}\mathrm{m}$ meter with free-slip BCs at the top. The still-water depth is $h=$0.6\text{,}\mathrm{m}, the wave heigh is $H=$0.112\text{\,}\mathrm{m} and the wave length is $\lambda=$1.204\text{,}\mathrm{m}, leading to a wave number of $k=$5.216\text{\,}{\mathrm{m}}^{-1}, a wave period of $T=$0.843\text{,}\mathrm{s} and a phase speed of $c=$1.429\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}. The situation is shown in figure 2 with a surface-piercing vertical cylinder of radius $R=$3\text{,}\mathrm{cm}. Nondimensional parameters are $kR=$0.16, $kh=$0.16 and $kH=$0.58. The order of the Fenton stream function wave in equation 25 is $N=$10$$.

The mesh is shown in figure 9. The mesh cells in the wave propagation zone have characteristic lengths of about $2.5\text{\,}\mathrm{cm}$, giving approximately 4.5 elements per wave height. To test wave propagation through the domain a similar mesh is used with the same refinement around the cylinder location, but without the cylinder being present. Symmetry boundary conditions are applied to the centre longitudinal plane shown prominently in figure 9. The inlet and outlet boundary conditions with associated forcing zones are implemented as described in section 5. The cylinder, when it is included, is non-slip in the horizontal plane, but the vertical velocity component is allowed to slip to avoid the free surface sticking to the cylinder. This approximation is done since our mesh resolution is not fine enough to model the true wetting dynamic with appropriate slip lengths. The longitudinal wall away from the cylinder is modelled as free-slip boundary, the same is true for the top and bottom surfaces of the tank.

The difference between the target wave elevation from the Fenton stream function and the VOF free surface is computed based on the intersection of the $\alpha=0.5$ isosurface and a centred longitudinal plane through the numerical domain without the cylinder. The region from $\nicefrac{{2}}{{3}}\,\lambda$ in front of the cylinder position to $\nicefrac{{1}}{{3}}\,\lambda$ behind the cylinder position is used to compute the difference between the phase of the numerical free surface and the stream function’s free surface. The diffusive error is computed in the same region as $\frac{1}{H}\sqrt{\langle\eta^{\prime 2}\rangle}$ where $\eta^{\prime}$ is the difference between the VOF and the stream-function wave elevation after correcting for the phase error and $\langle\cdot\rangle$ denotes the mean over about $300$ points where $\eta^{\prime}$ is sampled. Figure 10 shows the results compared to Paulsen et al. (2014, figures 2 and 3), but note that our results are from a 3D domain with forcing zones as described in section 5, while their results are from 2D simulations in a periodic domain which will let errors develop more easily over time.

Our numerical results are compared to the experimental results by Grue and Huseby (2002) in figures 11 and 12. The time series from figures 3c and 3d in their work have been shifted in time by $-15.63\text{\,}\mathrm{s}$ to align the wave phases. Figure 11 shows the wave elevation and also includes the Fenton stream-function solution to show the phase and diffusive errors in a more intuitive way. This also confirms that the applied stream-function wave is similar to the wave obtained in the laboratory wave flume. The inline force shown in figure 12 is computed as the total force on the cylinder in the longitudinal direction, positive towards the outlet. Note the bumps in the force signal right before each trough, the “secondary load cycle”. The associated free-surface ridge behind the cylinder at this time is shown in figure 13.

8 Discussion

Both test cases presented in section 7 show that the Ocellaris DG FEM solver is able to simulate complex two-phase flows in 3D. The pressure probes in the dam breaking test show good agreement with the experiments. Pressure probe P1 matches very well and the convergence between the mesh resolutions is clear. The impact pressures higher up on the vertical front wall of the container are less peaked than the measurements. They are most likely softened by the diffusive free surface, which can be seen in figure 7. The surface height probes match well, except for the probe behind the container where the flow is very complex and the free surface is multi valued. For the second test case the final comparison of the inline force on the cylinder shows an excellent match and the expected free-surface ridge behind the cylinder is observed.

The comparatively primitive free-surface capturing method is the weak point of the current implementation. Using a piecewise constant density field on the same mesh as the piecewise quadratic convecting velocity means that a lot of the data is thrown away when it comes to the advection of the free surface. Quantities directly related to the jump in the density field will always be limited to first-order convergence, but a more sophisticated free-surface capturing method could provide a sharper profile, hopefully eliminating the softening of the impact pressures in the first test case and removing much of the diffusive error in the wave propagation test. The excellent match of the inline force in the second test case should be treated with some scepticism considering the observed lowering of the peak height of the incident wave. The diffusive interface may have counteracted the expected lowering of the peak force by increasing the air phase pressures through an artificially high fluid density near the free surface.

The reason we selected to use only 4.5 elements per wave height for the wave propagation test is that the velocity function space has $30$ unknowns per mesh cell, so it is necessary to limit the number of cells more than in lower-order methods. The simulation has about $120\,000$ cells and is run on $32$ CPUs, so significantly increasing the mesh density will come at a large cost. This highlights the need for more research into free-surface capturing schemes that are designed specifically for high-order methods.

The main advantage of the method is the potential for using fewer elements away from the free surface due to the quadratic approximating polynomials. The dependency on having a very high mesh resolution in the free-surface region must be removed if the method is to be a preferable solution for general free-surface simulations. Another approach would be to make the method $h$–$p$ adaptive. A relevant issue for further research, which we have not touched on here, is that we observe that the number of iterations in the Krylov solver for the momentum equation increases a lot when the density ratio increases. The momentum equation can most likely be preconditioned or stabilised to stop this tendency.

We have shown that it is viable to use Ocellaris’ exactly divergence-free high-order discontinuous Galerkin finite element method with velocity slope limiting to simulate realistic air/water free-surface physics, but some work remains to be fully competitive with established low-order methods.

Acknowledgements

The simulations were performed on resources provided by UNINETT Sigma2, the National Infrastructure for High Performance Computing and Data Storage in Norway.