1  Modeling

This chapter describes the physical models used in CERF to simulate flows (without erosion in the version V1.0), as well as how these models are implemented in the code. It is divided into three sections. The first section presents the finite-volume method in the general case. The second section describes the Saint-Venant shallow-water flow model. The third section presents the Eulerian two-fluid flow model used in CERF.

1.1 The finite-volume method

The finite volume method is a numerical method used to solve partial differential equations. This method is used in CERF to solve the Saint-Venant equations and the two-fluid flow equations.

1.1.1 Principle of the finite volume method

In this section, we briefly recall the finite volume approximation, in dimension \(n_{dim}\), of a system of conservation equations of the form: \[ \frac{\partial \overrightarrow{w}}{\partial t} + \nabla \cdot \overline{\overline{F}}(\overrightarrow{w}) = \overrightarrow{S}(\overrightarrow{w}) \tag{1.1}\] where \(\overrightarrow{w}(\overrightarrow{x},t)\) is the vector of conservative variables of size \(n_{dof}\) depending on space and time, \(\overline{\overline{F}}\) the flux matrix of size \(n_{dof} \times n_{dim}\), \(\overrightarrow{S}\) the vector of source terms of size \(n_{dof}\), and \(\nabla \cdot\) the divergence operator.

The system is solved in \(\mathbb{R}^{n_{dim}}\). The computational domain \(\Omega \subset \mathbb{R}^{n_{dim}}\) is divided into a set of control volumes, also called cells, \(\Omega = \cup_{k} C_k\). We begin by integrating Equation 1.1 over a cell \(C_k\) using Green’s formula:

\[\int _{C_k} \frac{\partial \overrightarrow{w}}{\partial t} d \Omega + \int _{C_k} \nabla \cdot \overline{\overline{F}} d \Omega = \int _{C_k} \overrightarrow{S}(\overrightarrow{w}) d \Omega\]

\[ \int _{C_k} \frac{\partial \overrightarrow{w}}{\partial t} d \Omega + \sum _{a} \int _{\partial C_{k/a}} \overline{\overline{F}} \cdot \overrightarrow{n}_{k/a} d \partial\Omega = \int _{C_k} \overrightarrow{S}(\overrightarrow{w}) d \Omega \tag{1.2}\]

where \(\overrightarrow{n}_{k/a}\) denotes the unit normal vector on the boundary \(\partial C_{k/a}\) between cell \(C_k\) and the neighboring cell \(C_a\). The summation is carried out over all the neighboring cells of \(C_k\).

The finite volume method then consists in approximating the solution by piecewise constant functions per cell, such that by setting \[\overrightarrow{w}_k(t) \approx \frac {1}{|C_k|}\int _{C_k} \overrightarrow{w}(\overrightarrow{x},t) d \Omega\] where \(|C_k|\) represents the volume of cell \(C_k\), we approximate Equation 1.2 by \[ |C_k| \frac{\partial \overrightarrow{w}_k(t)}{\partial t} + \sum _{a} |\partial C_{k/a}| \: \widehat{\overrightarrow{F}} \left(\overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a} \right) = |C_k|S_k(t) \tag{1.3}\]

where \(S_k(t)\approx \frac {1}{|C_k|}\int _{C_k}S(\overrightarrow{x},t) d \Omega\) represents the average value of the source term vector over the cell, \(|\partial C_{k/a}|\) the surface area of the face between cells \(C_k\) and \(C_a\), and \(\widehat{\overrightarrow{F}}\) an approximation of the flux.

Since the unknowns \(\overrightarrow{w}\) are not defined at the interface between cells, the flux cannot be computed directly. It is therefore approximated by a numerical flux depending on the interface and on the values of the variables in the adjacent cells. We denote:

\[\overrightarrow{F}_{k/a}(t) = \widehat{\overrightarrow{F}} \left(\overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a} \right) \approx \frac{1}{|\partial C_{k/a}|}\int _{\partial C_{k/a}} \overline{\overline{F}} \cdot \overrightarrow{n}_{k/a} ds\]

There are many possible approaches: Rusanov, HLLC, etc. In CERF, the numerical flux is computed using Godunov’s method, based on the exact solution of the Riemann problem at the interface \(k/a\) (for more details, see Section 1.2.5 and Section 1.3.5).

Second-order accuracy in space can be achieved using the MUSCL (Monotone Upstream Scheme for Conservation Laws) method to reconstruct the primitive variables at the interfaces. From the values of the variables on the neighboring cells, the gradient of the primitive variables is estimated. Let \({w_p}_i\) be the \(i\)-th component of the primitive variables and \({w_p}_i^k\) its value at the center of the cell; the gradient over cell \(C_k\) can be estimated as: \[ |C_k| \widetilde{\nabla} {w_p}_i = \int _{C_k} \nabla {w_p}_i d\Omega =\sum _{a} \int _{\partial C_{k/a}} w_{p_i} \overrightarrow{n}_{k/a} d\partial\Omega =\sum_{a} |\partial C_{k/a}| \frac{1}{2}\left( {w_p}_i^k + {w_p}_i^a \right) \overrightarrow{n}_{k/a}\] and the primitive variables of cell \(C_k\) with center \(K\) can then be reconstructed at the centers \(F\) of the cell faces such that: \[w_{p_i}^F = w_{p_i}^k + \widetilde\nabla w_{p_i} \cdot \overrightarrow{KF}\]

To achieve second-order accuracy, \(\widehat{\overrightarrow{F}} \left(\overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a} \right)\) is then replaced by \(\widehat{\overrightarrow{F}} \left(\overrightarrow{w}_k^F(t),\overrightarrow{w}_a^F(t),\overrightarrow{n}_{k/a} \right)\)

Note

In general, a limiter must be combined with the MUSCL method. In CERF, the Barth-Jespersen limiter is used: \(Min(w_{p_F},w_{p_F}) \leq w_{p_i}^F \leq Max(w_{p_F},w_{p_F})\).

1.1.2 Time discretization

To facilitate the parallelization of the code, CERF uses only explicit schemes: the first- or second-order Runge-Kutta scheme. The time interval \(\left[0,T\right]\) is discretized such that \(t_{n+1} = t_n + \delta t_n\), where \(\delta t_n\) represents the time step at time \(t_n\). The first-order Runge-Kutta method, also known as the first-order explicit Euler method, consists in rewriting Equation 1.3 at time \(t_n\) by approximating the time derivative as a finite difference, such that: \[ |C_k|\frac{ \overrightarrow{w}_k(t_{n+1}) - \overrightarrow{w}_k(t_{n})}{\delta t_n} + \sum _{a} |\partial C_{k/a}| \: \overrightarrow{F}_{k/a}(t_n) = |C_k|\overrightarrow{S}_k(t_n)\]

By denoting \(\overrightarrow{w}_k(t_{n})=\overrightarrow{w}_k^{n}\), the new value of the conservative variables at time \(t_{n+1}\) is thus computed explicitly: \[ \overrightarrow{w}_k^{n+1} = \overrightarrow{w}_k^{n}+\frac{\delta t_n}{|C_k|} \sum _{a} |\partial C_{k/a}| \: \overrightarrow{F}_{k/a}^n - \delta t_n \overrightarrow{S}_k^n\]

The second-order Runge-Kutta method consists in computing the solution at time \(t_{n+1}\) in two steps. The first step consists in computing the solution at time \(t_{n+\frac{1}{2}}\) using the first-order explicit Euler scheme. The second step consists in computing the solution at time \(t_{n+1}\) using the first-order explicit Euler scheme together with the solution obtained at time \(t_{n+\frac{1}{2}}\). \[ \overrightarrow{w}_k^{n+\frac{1}{2}} = \overrightarrow{w}_k^{n}+\frac{\delta t_n}{2 |C_k|} \sum _{a} |\partial C_{k/a}| \: \overrightarrow{F}_{k/a}^n - \frac{\delta t_n}{2} \overrightarrow{S}_k^n \] \[ \overrightarrow{w}_k^{n+1} = \overrightarrow{w}_k^{n}+\frac{\delta t_n}{|C_k|} \sum _{a} |\partial C_{k/a}| \: \overrightarrow{F}_{k/a}^{n+\frac{1}{2}} - \delta t_n\overrightarrow{S}_k^{n+\frac{1}{2}} \]

Important 1.1: Courant-Friedrichs-Lewy stability condition

Explicit methods are simple to implement, but their stability is constrained by the Courant-Friedrichs-Lewy (CFL) condition. The CFL condition is a necessary condition to guarantee the stability of the numerical scheme. This condition depends on the wave propagation speed in the medium and on the size of the cells. In practice, it is determined for each hyperbolic model considered. In CERF, the cell size \(d_k\) is determined by the ratio between the volume of the cell and the surface area of its boundary: \[ d_k = \frac{|C_k|}{\sum_{a}|\partial C_{k/a}|}\]
In the 1D or 2D case, symmetry faces are not taken into account when computing the boundary surface area. The CFL condition is then given by: \[ \delta t_n \leq \alpha_{CFL} \frac{\min_{k}d_k}{\max_{k} \lambda_k}\]
where \(\alpha_{CFL}\) is a coefficient to be defined by the user and \(\lambda_k\) is the wave propagation speed in cell \(C_k\), depending on the model used (Note 1.1 for Saint-Venant and Note 1.2 for two-fluid flow).

1.1.3 Spatial discretization

Important

In its V1.0 version, CERF uses exclusively hexahedral cells and therefore quadrangular faces! Generalization is currently under development.

Note

CERF allows the use of dynamic meshes, meaning that the mesh can be regularly refined or coarsened (AMR: Adaptive Mesh Refinement).

CERF is a natively three-dimensional code. In its initial version, CERF uses only hexahedral cells and therefore quadrangular faces. Cells and faces are not necessarily regular. Nevertheless, as illustrated in Figure 1.1, two-dimensional simulations can be performed by creating a mesh with a single cell in the \(z\) direction and imposing “mirror”-type conditions on the faces in the \(x,y\) plane. Similarly, one-dimensional simulations can be performed by creating a mesh with a single cell in the \(y\) and \(z\) directions and imposing “mirror”-type conditions on the faces in the \(x,z\) and \(x,y\) planes.

Figure 1.1: Three-, two-, and one-dimensional meshes with hexahedral cells

1.1.3.1 Meshing strategy

CERF provides mesh refinement and coarsening tools (AMR: Adaptive Mesh Refinement). To avoid penalizing computation time, the mesh is not modified at every time step, but at instants determined by the user. This approach makes it possible to maintain good mesh quality while limiting computational cost. This method is referred to as BB-AMR (Block Based Adaptive Mesh Refinement). The various steps of the meshing strategy used in CERF are described below.

Figure 1.2: Initial conforming mesh: the “Master Mesh”.

Step 1: For the sake of simplicity, the mesh shown in Figure 1.2 consists of regular quadrangles, which can represent a hexahedral mesh viewed along \(z\) with a single cell in the \(z\) direction. As illustrated in Figure 1.2, we start from an initial conforming mesh made up solely of hexahedra of arbitrary shape. In a simple case, the mesh can be created by CERF; otherwise it can be imported in GMSH V2.2 format. This mesh is referred to as the master mesh. It is essential that the mesh be conforming in order to easily define the faces connecting the cells. The cells of this initial mesh are referred to as blocks.

Figure 1.3: Numbering of the initial mesh according to Morton codes.

Step 2: The blocks of the initial mesh are ordered by numbering them according to the Z-order or Morton code (Figure 1.3). With this approach, known as a Space Filling Curve (Section 1.1.3.2), the numbering makes it possible to traverse the blocks, even in three dimensions, following a numbering such that successive blocks in the numbering are, most of the time, neighbors in space.

Figure 1.4: Definition of the block refinement level.

Step 3: The refinement level of the blocks is then defined, as illustrated in Figure 1.4. Initially, the refinement level is defined by the user. During the computation, the refinement level is defined by a criterion Section 1.1.3.3. As verified in Ersoy, Golay, and Yushchenko (2013) and Thomas Altazin and Yushchenko (2016), the difference in refinement level between two blocks must not exceed \(1\). During refinement, the block is divided by \(2\) in each direction, for each level required. Thus, denoting by \(Niv\) the required refinement level (where \(Niv=0\) represents a block that is not refined) and \(dim=1,2,3\) the dimension of the problem, the number of cells in a block will be: \(2^{dim \times Niv}\). In the example shown in Figure 1.4, block \(1\) is not refined, block \(2\) is refined by \(1\) level, and block \(4\) by \(2\) levels. This leads, as illustrated in Figure 1.5 for the case \(n_{dim}=2\), to \(1\) cell in block \(1\), \(4\) cells in block \(2\), and \(16\) cells in block \(4\).

Figure 1.5: Decomposition of the mesh into domains.

Step 4: Since the time integration schemes in CERF are explicit, it is straightforward to parallelize the solution algorithm. To do this, the mesh must be split into several domains, ideally with a balanced number of cells and minimal interfaces, in order to minimize data exchange between domains. Many algorithms exist, some of them quite complex, that perform this operation efficiently. In CERF, we chose a compromise between simplicity and efficiency. We opted for a decomposition driven by Morton codes. Knowing the number of cells per block, the blocks are traversed following the Morton numbering and assigned to a domain. The number of cells per domain is kept balanced. As illustrated in Figure 1.5, the initial mesh is split into 3 domains made up of 25, 21, and 18 cells. Each computational domain is then assigned to a computing process.

Figure 1.6: Mesh generation.

Step 5: The final step consists in generating the mesh. To exchange information between domains, CERF creates ghost cells. These cells are copies of the neighboring cells of the domain. They are used to compute the numerical fluxes at the interfaces between domains. As illustrated in Figure 1.6, ghost cells are shown in gray. Note that the mesh generated in this way is not conforming! For example, the first cell at the bottom left of domain \(1\) in Figure 1.6 is surrounded by 6 faces (in the \(x,y\) plane). The refinement of the blocks defining the mesh naturally depends on the dimension of the problem being treated (Figure 1.7). Finally, if the mesh is modified during the computation, it is necessary to project the conservative variables onto the new mesh. In a remeshing operation, the level can only increase or decrease by at most \(1\), so that the projection step is relatively simple.

Figure 1.7: Refinement depending on the dimension of the problem.

1.1.3.2 Morton code

Space-filling curves are curves that make it possible to traverse a higher-dimensional space. First introduced by Peano, then by Hilbert or Lebesgue, they are used for example to index data, and in particular in the construction of quadtrees or octrees. In CERF, numbering by Morton code, also known as the Z-order curve, is used. This numbering makes it possible to index a set of data (points, cells, faces, etc.) such that data that are neighbors in space are also neighbors in the indexing. This property is used in CERF to split the mesh into computational domains by indexing the blocks of the master mesh. Consider a point with coordinates \((x,y,z)\) in a three-dimensional space. To index this point, the coordinates are first converted into positive integers. To do this, the smallest coordinates of the study domain \((x_{min},y_{min},z_{min})\) are defined, along with a distance \(d_{morton}\) defining the grid of the domain, for example one tenth of the size of the smallest cell in the domain. The integer coordinates of the point \((ix,iy,iz)\) are then defined by: \[ix = \frac{x-x_{min}}{d_{morton}} ,\: iy = \frac{y-y_{min}}{d_{morton}} ,\: iz = \frac{z-z_{min}}{d_{morton}} \] Morton numbering then consists in converting these integer coordinates into base \(2\) \((ix_2,iy_2,iz_2)\), then interleaving the bits of each coordinate to form a single base-\(2\) integer \(M_2\). Finally, converting back to base \(10\), the Morton index \(M\) of the point \((x,y,z)\) is obtained.

Example 1.1 \(ix=2\), \(iy=5\), \(iz=0\).
We then have \(ix_2=010\), \(iy_2=101\), and \(iz_2=000\).
Interleaving the bits gives \(M_2=010100010\), and in base \(10\), \(M=162\).

Note

All of CERF’s indexing tools can be found in the module ~/sources/UTI/mod_zorder.f90.

1.1.3.3 Refinement criterion

In general, a refinement criterion is defined for each physical model considered. This criterion is often derived from physical considerations based, for example, on the gradient of one of the primitive variables, or, more relevantly and more elaborately, on an a posteriori error estimate or the construction of an error indicator. In the case of hyperbolic problems, there exists a more general criterion derived from the work of Croisille (1990). This criterion is based on the numerical production of entropy. It was used by Golay (2009), Ersoy, Golay, and Yushchenko (2013), and studied at length by Puppo (2004). The entropy inequality of Lax (1971) allows us to state that for any hyperbolic system of conservation laws of the form Equation 1.1, there exists a mathematical entropy \(s(\overrightarrow{w})\) and an entropy flux \(\overrightarrow{\psi}_s(\overrightarrow{w})\) such that: \[\frac{\partial s(\overrightarrow{w})}{\partial t} + \nabla \cdot \overrightarrow{\psi}_s(\overrightarrow{w}) \leq 0 \tag{1.4}\] where the entropy flux \(\overrightarrow{\psi}_s(\overrightarrow{w})\) must satisfy a compatibility equation \(\nabla_w \psi_s = \nabla_w \overline{\overline{F}} \cdot \nabla_w s\). For smooth solutions, Equation 1.4 is an equality. In the case of a shock, it is a strict inequality. However, in numerical simulations, the numerical production of entropy is observed not to be zero! It turns out that the numerical production of entropy is a very good error indicator, except at shocks, where it tends to infinity. We therefore define \(p^s_k\), the numerical entropy production over cell \(C_k\), as:

\[p^s_k = \frac{1}{|C_k|} \int_{C_k} \left( \frac{\partial s(\overrightarrow{w})}{\partial t} + \nabla \cdot \overrightarrow{\psi}_s(\overrightarrow{w}) \right) d\Omega\].

The refinement criterion \(C_b\) is defined per block as the ratio of the numerical entropy production in the block to the total numerical entropy production over the domain. \[ C_b = \frac{\int_{Bloc} \left( \frac{\partial s(\overrightarrow{w})}{\partial t} + \nabla \cdot \overrightarrow{\psi}_s(\overrightarrow{w}) \right) d\Omega}{\int_{\Omega} \left( \frac{\partial s(\overrightarrow{w})}{\partial t} + \nabla \cdot \overrightarrow{\psi}_s(\overrightarrow{w}) \right) d\Omega}\] In CERF, the refinement criterion is defined per block. If the criterion exceeds a threshold value \(\alpha_{max}\), the block is refined. If the criterion is below a threshold value \(\alpha_{min}\), the block is coarsened. The threshold values are defined by the user.

1.2 Saint-Venant model

1.2.1 Governing equations

The shallow-water flow model, or Saint-Venant model, is a simplified flow model. It is derived from the Navier-Stokes equations by neglecting the viscosity terms, assuming a hydrostatic pressure distribution, and assuming that the water depth is much smaller than the characteristic length of the flow. The original Barré de Saint-Venant model is a 1D model Saint-Venant (1871); the derivation of the model is explained, for example, by Marche (2007) or Poussel (2024). Development of the model in CERF was initiated by Pons (2018). In a frame of reference \((O,x,y,z)\), we consider a shallow-water flow over a domain where the bathymetry is defined by \(z_b(x,y)\), with a water depth \(h(x,y)\), so that the free-surface level is defined by \(\zeta=h+z_b\), as illustrated in Figure 1.8.

Figure 1.8: Notation of the shallow-water flow model

The Saint-Venant equations are as follows: \[ \frac{\partial h}{\partial t} + \nabla \cdot (h \overrightarrow{u}) = 0 \tag{1.5}\] \[ \frac{\partial h\overrightarrow{u}}{\partial t} + \nabla \cdot (h \overrightarrow{u} \otimes \overrightarrow{u} + g\frac{h^2}{2} \overline{\overline{I}} ) = -gh \nabla z_b -gh\overrightarrow{S_f} \tag{1.6}\]

where \(\overrightarrow{u}=<u,v>^T\) represents the flow velocity, \(g\) the acceleration due to gravity, \(S_f\) the friction term, and where \(\overline{\overline{I}}\) denotes the identity matrix, \(\nabla\) the gradient operator, \(\nabla \cdot\) the divergence operator, and \(\otimes\) the tensor product. We define the conservative variables \(\overrightarrow{w}=<h,hu,hv>^T\) and the primitive variables \(\overrightarrow{w_p}=<h,u,v>^T\). This corresponds to the general form Equation 1.1 with: \[\overline{\overline{F}}(\overrightarrow{w}) = \left[ \begin{array}{cc} hu & hv \\ hu^2 + g\frac{h^2}{2} & h u v \\ h u v & hv^2 + g\frac{h^2}{2} \end{array}\right]\]

and
\[\overrightarrow{S} = -gh \nabla z_b -gh\overrightarrow{S_f} \]

1.2.2 Handling bathymetry

Handling bathymetry in the finite-volume resolution scheme requires particular care. Indeed, the bathymetry jump inherent to the fact that bathymetry is considered constant per cell can generate spurious waves. The reconstruction method proposed by Audusse et al. (2004) is therefore implemented in CERF. Let \(\Delta z_{k/a}\) be the bathymetry jump at the interface between cells. At the interface, we then distinguish the flux \(\overrightarrow{F}_k\) contributing to cell \(C_k\) and the flux \(\overrightarrow{F}_a\) contributing to cell \(C_a\), such that:

\[\overrightarrow{F}_k \left( \overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a}, \Delta z_{k/a}\right)=\overrightarrow{F} \left(\overrightarrow{w}_k^*(t),\overrightarrow{w}_a^*(t),\overrightarrow{n}_{k/a} \right) + \left\{ \begin{array}{cc} 0 \\ \frac{g}{2} \left(h_k^2 - {h_k^*}^2 \right) \overrightarrow{n}_{k/a} \end{array} \right\} \] \[\overrightarrow{F}_a \left( \overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a}, \Delta z_{k/a}\right)=\overrightarrow{F} \left(\overrightarrow{w}_k^*(t),\overrightarrow{w}_a^*(t),\overrightarrow{n}_{k/a} \right) + \left\{ \begin{array}{cc} 0 \\ \frac{g}{2} \left(h_a^2 - {h_a^*}^2 \right) \overrightarrow{n}_{k/a} \end{array} \right\} \] where \(\overrightarrow{w}_k^*=<h_k^*,\overrightarrow{u}_k>^T\) and \(\overrightarrow{w}_a^*=<h_a^*,\overrightarrow{u}_a>^T\), with: \[h_k^* = max \left(0, h_k - max \left( 0, \Delta z_{k/a}\right) \right)\] \[h_a^* = max \left(0, h_a - max \left( 0, \Delta z_{k/a}\right) \right)\]

1.2.3 Refinement criterion

For the Saint-Venant model, to define the numerical entropy production Equation 1.4, we use the entropy defined by: \[ s(\overrightarrow{w}) = \frac{1}{2} \left( h ||\overrightarrow{u}||^2 + gh^2 + 2 g h z_b \right) \] and the entropy flux \[ \overrightarrow{\psi}_s(\overrightarrow{w}) = \left( s + \frac{g h^2}{2} \right) \overrightarrow{u} .\]

1.2.4 Handling friction

In the momentum conservation equation Equation 1.6, the friction term is taken into account. In CERF, the friction term is modeled using either the Manning-Strickler term or the Darcy-Weisbach term. The Manning-Strickler term is defined by: \[ \overrightarrow{S_f} = C_f \frac{||\overrightarrow{u}|| \overrightarrow{u}}{h^{4/3}}\] where \(C_f = n^2\) is the Manning-Strickler friction coefficient.

The Darcy-Weisbach term is defined by: \[ \overrightarrow{S_f} = C_f \frac{||\overrightarrow{u}|| \overrightarrow{u}}{h}\]

The friction term is handled using a “splitting” method. We first solve problem Equation 1.1 without the friction term, \[ \frac{\partial \overrightarrow{w}}{\partial t} + \nabla \cdot \overline{\overline{F}}(\overrightarrow{w}) = -gh \nabla z_b \tag{1.7}\]
then we solve \[ \frac{\partial \overrightarrow{w}}{\partial t} = -gh\overrightarrow{S_f} \tag{1.8}\]

We then assume that friction modifies neither the water depth nor the direction of the flow, but only the magnitude of the flow velocity. The friction problem for the velocity can then be solved using a first-order implicit scheme over a time step \(\delta t_n\) and the velocity “penalized” as in Varra et al. (2024).
Let \(\overrightarrow{u}^*\) be the solution of problem Equation 1.7; problem Equation 1.8 can then be solved and the velocity penalized by: \[ \overrightarrow{u} = \frac{2\overrightarrow{u}^*}{1+\sqrt{1 + 4\beta}} \]
where
\[\beta = \delta t_n g C_f h^{-q}\]

with \(q=4/3\) for the Manning-Strickler model and \(q=1\) for the Darcy-Weisbach model.

1.2.5 Riemann solver for the Saint-Venant model

1.2.5.1 Overview

Numerical fluxes in CERF are computed using a Riemann solver. There are many references concerning this solver, such as Toro (1997). This section describes the Riemann solver used for the Saint-Venant model, based on the course by Gloria Faccanoni.

We therefore seek to solve the Saint-Venant problem in one dimension. We consider a shallow-water flow with depth \(h(x,t)\) and horizontal velocity \(u(x,t)\) over a flat bottom, described by the following problem: \[ \left \{ \begin{array}{ll} \displaystyle\frac{\partial h}{\partial t} + \frac{\partial (hu)}{\partial x} = 0 \\ \displaystyle\frac{\partial (hu)}{\partial t} + \frac{\partial (hu^2 + \frac{1}{2}gh^2)}{\partial x} = 0 \end{array} \right. \tag{1.9}\] With initial conditions such that: \[ \overrightarrow{w_p}(x,0) =\left\{ \begin{array}{ll} \overrightarrow{w_p}_L &\mbox{ if } x<0 \\ \overrightarrow{w_p}_R &\mbox{ if } x>0 \end{array} \right. \] where \(\overrightarrow{w_p}_L=<h_L,u_L>^T\) and \(\overrightarrow{w_p}_R=<h_R,u_R>^T\) are two constant states.

By expanding Equation 1.9, the system can be rewritten in the following form: \[ \frac{\partial}{\partial t} \left\{ \begin{array}{c} h \\ u \end{array} \right\} + \left[ \begin{array}{cc} u & h \\ g & u \end{array} \right] \frac{\partial}{\partial t} \left\{ \begin{array}{c} h \\ u \end{array} \right\} = \left\{ \begin{array}{c} 0 \\ 0 \end{array} \right\} \]

We then seek the eigenvalues \(\lambda_1\) and \(\lambda_2\) as well as the eigenvectors \(\overrightarrow{r}_1\) and \(\overrightarrow{r}_2\) of the matrix thus obtained. A quick calculation gives: \[ \lambda_1=u-\sqrt{gh} , \lambda_2=u+\sqrt{gh} \] with \(\lambda_1 \leq \lambda_2\) and \[\overrightarrow{r}_1=\left\{ \begin{array}{c} -\sqrt{h} \\ \sqrt{g} \end{array} \right\} , \overrightarrow{r}_2=\left\{ \begin{array}{c} \sqrt{h} \\ \sqrt{g} \end{array} \right\}\].

Note 1.1

The eigenvalues are used to determine the wave propagation speed within the cell Important 1.1 as: \(\lambda_k = Max\left( |u-\sqrt{gh}| , |u+\sqrt{gh}| \right)\)

To determine the nature of the two characteristic fields, we compute \(\nabla_{w_p}\lambda_1 \cdot \overrightarrow{r}_1 = 3/2 \neq 0\) and \(\nabla_{w_p}\lambda_2 \cdot \overrightarrow{r}_2 = 3/2 \neq 0\). The characteristic fields are therefore genuinely nonlinear. There will therefore be no contact discontinuity, and the two waves are either shocks or rarefaction waves.

Since there are \(2\) waves, \(1\) Riemann invariant \(I_{1;1}\), \(I_{2;1}\) can be defined for each wave, such that \(\nabla_{w_p}I_{1;1} \cdot \overrightarrow{r}_1 = 0\) and \(\nabla_{w_p}I_{2;1} \cdot \overrightarrow{r}_2 = 0\). It is easy to verify that \(I_{1;1}=u+2\sqrt{gh}\) and \(I_{2;1}=u-2\sqrt{gh}\) are Riemann invariants.

The entropy weak solution of the Riemann problem consists of at most three constant states \(\overrightarrow{w_p}_L\), \(\overrightarrow{w_p}_*\), \(\overrightarrow{w_p}_R\), separated by two discontinuities. These discontinuities are either shocks or rarefaction waves. In what follows, we determine the state \(\overrightarrow{w_p}_*\).

1.2.5.2 Study of the rarefaction waves

Consider a state \(\overrightarrow{w_p}_g\) and a state \(\overrightarrow{w_p}_d\), connected by a k-rarefaction wave. The following conditions must then be satisfied:

  • \(\lambda_k(\overrightarrow{w_p}_g)<\lambda_k(\overrightarrow{w_p}_d)\)
  • The Riemann invariant is conserved along the rarefaction wave \(I_{k;1}(\overrightarrow{w_p}_g)=I_{k;1}(\overrightarrow{w_p}_d)\)

Let us apply these conditions to the Saint-Venant problem.

Consider first a 1-rarefaction connecting a state \(\overrightarrow{w_p}_L\) to a state \(\overrightarrow{w_p}_*\), separated by a rarefaction wave. The Riemann invariants give \[ I_{k;1}(\overrightarrow{w_p}_L)=I_{k;1}(\overrightarrow{w_p}_*) \rightarrow u_L+2\sqrt{g\ h_L}=u_*+2\sqrt{g\ h_*} \] that is \[ u_* = u_L + 2\sqrt{g} \left( \sqrt{h_L} - \sqrt{h_*}\right) \tag{1.10}\] Furthermore, \[\lambda_k(\overrightarrow{w_p}_L)<\lambda_k(\overrightarrow{w_p}_*) \rightarrow u_L-\sqrt{g\ h_L} < u_*-\sqrt{g\ h_*} \rightarrow h_* < h_L\] and \[ u_* > u_L\]

Consider finally a 2-rarefaction connecting a state \(\overrightarrow{w_p}_*\) to a state \(\overrightarrow{w_p}_R\), separated by a rarefaction wave. The Riemann invariants give \[ I_{k;1}(\overrightarrow{w_p}_*)=I_{k;1}(\overrightarrow{w_p}_R) \rightarrow u_*-2\sqrt{g\ h_*}=u_R-2\sqrt{g\ h_R} \] that is \[ u_* = u_R + 2\sqrt{g} \left(\sqrt{h_*} - \sqrt{h_R}\right) \tag{1.11}\] Furthermore, \[\lambda_k(\overrightarrow{w_p}_*)<\lambda_k(\overrightarrow{w_p}_R) \rightarrow u_*+\sqrt{g\ h_*} < u_R+\sqrt{g\ h_R} \rightarrow h_* < h_R\] and \[ u_* < u_R\]

1.2.5.3 Study of the shocks

Consider a state \(\overrightarrow{w_p}_g\) and a state \(\overrightarrow{w_p}_d\), connected by a \(k\) shock wave. The following conditions must then be satisfied:

  • The Rankine-Hugoniot conditions, where \(\dot\sigma_k\) denotes the shock speed: \(\displaystyle \dot\sigma_k=\frac{\left[\overrightarrow{F}(\overrightarrow{w})\right]}{\left[\overrightarrow{w}\right]}\)
  • The Lax entropy conditions for a genuinely nonlinear field \(k\): \[ \left\{ \begin{array}{l} \lambda_{k-1}(\overrightarrow{w_p}_g) < \dot\sigma_k \\ \lambda_{k}(\overrightarrow{w_p}_d) < \dot\sigma_k < \lambda_{k}(\overrightarrow{w_p}_g) \\ \dot\sigma_k < \lambda_{k+1}(\overrightarrow{w_p}_d) \end{array} \right. \]

Applied to the Saint-Venant problem, the Rankine-Hugoniot conditions give: \[ \dot\sigma_k = \frac{h_d u_d - h_g u_g}{h_d - h_g} = \frac{h_d u_d^2 + \frac{g}{2}h_d^2 - h_g u_g^2 + \frac{g}{2}h_g^2}{h_d u_d - h_g u_g} \] in particular, we can define \(j_k\) such that

\[ j_k = h_g u_g - \dot\sigma_k h_g = h_d u_d - \dot\sigma_k h_d .\]

By eliminating the shock speed, we obtain: \[\left( u_g - u_d \right)^2 = \frac{g}{2 h_d h_g} \left( h_g - h_d \right)^2 \left( h_g + h_d \right)\] and therefore \[ j_k = \frac{g}{2} \frac{h_g - h_d}{u_d - u_g} \]

Consider first a 1-shock connecting a state \(\overrightarrow{w_p}_L\) to a state \(\overrightarrow{w_p}_*\), separated by a shock wave. The Lax conditions give: \[ \displaystyle\left\{ \begin{array}{l} \lambda_1(\overrightarrow{w_p}_*) < \dot\sigma_1 < \lambda_1(\overrightarrow{w_p}_L) \\ \dot\sigma_1 < \lambda_2(\overrightarrow{w_p}_*) \end{array} \right. \rightarrow \left\{ \begin{array}{l} u_*-\sqrt{g h_*} < \dot\sigma_1 < u_L-\sqrt{g h_L} \\ \dot\sigma_1 < u_*+\sqrt{g h_*} \end{array} \right.\]

The Rankine-Hugoniot conditions give: \[ \displaystyle\left\{ \begin{array}{l} j_1 = h_*\left( u_*-\dot\sigma_1\right) = h_L\left( u_L-\dot\sigma_1\right) \\ \displaystyle\dot\sigma_1 = \frac{h_* u_*^2 + \frac{g}{2}h_*^2 - h_L u_L^2 - \frac{g}{2}h_L^2}{h_* u_* - h_L u_L} \end{array} \right.\]

So, from the Lax conditions and using the first Rankine-Hugoniot condition, we obtain: \[ \displaystyle\left\{ \begin{array}{l} u_*-\sqrt{g h_*} < \dot\sigma_1 \\ u_L-\sqrt{g h_L} > \dot\sigma_1 \\ u_*+\sqrt{g h_*} > \dot\sigma_1 \end{array} \right. \rightarrow \left\{ \begin{array}{l} h_*\sqrt{g h_*} > j_1 \\ h_L\sqrt{g h_L} < j_1 \\ - h_*\sqrt{g h_*} < j_1 \end{array} \right.\] That is \[ h_* > h_L \] and since \(j_1=\frac{g}{2} \frac{h_L - h_*}{u_* - u_L} > 0\), we have \[ u_* < u_L \] and using the expression for \(\dot\sigma_1\), we find the velocity \[ u_*=u_L+\left(h_L-h_*\right)\sqrt{\frac{g}{2}\frac{h_L + h_*}{h_L h_*}} \tag{1.12}\]

Consider finally a 2-shock connecting a state \(\overrightarrow{w_p}_*\) to a state \(\overrightarrow{w_p}_R\), separated by a shock wave. The Lax conditions give: \[ \displaystyle\left\{ \begin{array}{l} \lambda_1(\overrightarrow{w_p}_*) < \dot\sigma_2 \\ \lambda_2(\overrightarrow{w_p}_R) < \dot\sigma_2 < \lambda_2(\overrightarrow{w_p}_*) \end{array} \right. \rightarrow \left\{ \begin{array}{l} u_*-\sqrt{g h_*} < \dot\sigma_2 \\ < u_R+\sqrt{g h_R} < \dot\sigma_2 < u_*+\sqrt{g h_*} \end{array} \right.\]

The Rankine-Hugoniot conditions give: \[ j_2 = h_*\left( u_*-\dot\sigma_2 \right) = h_R \left( u_R-\dot\sigma_2 \right)= \frac{g}{2} \frac{h_* - h_R}{u_R - u_*}\]

So, from the Lax conditions and using the first Rankine-Hugoniot condition, we obtain: \[ \displaystyle\left\{ \begin{array}{l} u_*-\sqrt{g h_*} < \dot\sigma_2 \\ u_R+\sqrt{g h_R} < \dot\sigma_2 \\ u_*+\sqrt{g h_*} > \dot\sigma_2 \end{array} \right. \rightarrow \left\{ \begin{array}{l} h_*\sqrt{g h_*} > j_2 \\ -h_R\sqrt{g h_R} > j_2 \\ - h_*\sqrt{g h_*} > j_2 \end{array} \right.\]

That is \[ h_* > h_R \] and since \(j_2 < 0\), we have \[ u_* > u_L \] and using the expression for \(\dot\sigma_2\), we find the velocity \[ u_* = u_R + \left( h_* - h_R \right) \sqrt{\frac{g}{2} \frac{h_R + h_*}{h_R h_*}} \tag{1.13}\]

1.2.5.4 Constructing the solution

To determine the solution of the Riemann problem, the state \(\overrightarrow{w_p}_*\) must be determined as a function of the states \(\overrightarrow{w_p}_L\) and \(\overrightarrow{w_p}_R\). The intermediate state is connected to the left state by a 1-wave, using Equation 1.12 and Equation 1.10: \[ u_* = v^L(h_*) = \left\{ \begin{array}{l} u_L - \displaystyle\left( \sqrt{h_*} - \sqrt{h_L}\right) \frac{2\sqrt{g}}{\sqrt{h_L} + \sqrt{h_*}} \mbox{ if } h_* \leq h_L \mbox{ (1-rarefaction)}\\ u_L - \displaystyle\left( \sqrt{h_*} - \sqrt{h_L}\right) \sqrt{\frac{g}{2} \frac{h_* + h_L}{h_* h_L}} \mbox{ if } h_* > h_L \mbox{ (1-shock)}\end{array} \right. \] and the intermediate state is connected to the right state by a 2-wave, using Equation 1.11 and Equation 1.13: \[ u_* = v^R(h_*) = \left\{ \begin{array}{l} u_R + \displaystyle\left( \sqrt{h_*} - \sqrt{h_R}\right) \frac{2\sqrt{g}}{\sqrt{h_R} + \sqrt{h_*}} \mbox{ if } h_* \leq h_R \mbox{ (2-rarefaction)}\\ u_R + \left( \sqrt{h_*} - \sqrt{h_R}\right) \displaystyle\sqrt{\frac{g}{2} \frac{h_* + h_R}{h_* h_R}} \mbox{ if } h_* > h_R \mbox{ (2-shock)} \end{array} \right. \]

So to find \(h_*\), the equation \(v^L(h_*)-v^R(h_*)=0\) must be solved using a Newton-Raphson scheme.

By introducing the function \(\zeta(h,\psi)=\left\{ \begin{array}{l} \displaystyle\frac{2\sqrt{g}}{\sqrt{h} + \sqrt{\psi}} \mbox{ if } h \leq \psi \\ \displaystyle\sqrt{\frac{g}{2} \frac{h + \psi}{h \psi}} \mbox{ if } h > \psi \end{array} \right.\), we will simply write \[u_* = u_R + \left( \sqrt{h_*} - \sqrt{h_R}\right) \zeta(h_*,h_R)\]

The previous calculations allow us to define the speeds of the shock waves \(\dot\sigma_1 = u_L - h_* \displaystyle\sqrt{\frac{g}{2} \frac{h_L + h_*}{h_L h_*}}\) and \(\dot\sigma_2 = u_R + h_* \displaystyle\sqrt{\frac{g}{2} \frac{h_R + h_*}{h_R h_*}}\) .

and the speeds of the rarefaction waves \(\lambda_1 = u_L - \sqrt{g h_L}\) , \(\lambda_1^* = u_* - \sqrt{g h_*}\), \(\lambda_2^* = u_* + \sqrt{g h_*}\) and \(\lambda_2 = u_R + \sqrt{g h_R}\).

The solution in a rarefaction zone is given by: \[\overrightarrow{w_p}_{1det}(x,t) = \left\{ \begin{array}{l} \displaystyle \frac{1}{9g} \left( u_L+2 \sqrt{g h_L} - \frac{x}{t} \right)^2 \\ \displaystyle \frac{1}{3} \left( u_L + 2 \sqrt{g h_L} + 2 \frac{x}{t} \right) \end{array} \right\}\] and \[\overrightarrow{w_p}_{2det}(x,t) = \left\{ \begin{array}{l} \displaystyle \frac{1}{9g} \left( -u_R + 2 \sqrt{g h_R} + \frac{x}{t} \right)^2 \\ \displaystyle \frac{1}{3} \left( u_R - 2 \sqrt{g h_R} + 2 \frac{x}{t} \right) \end{array} \right\} .\]

We can now write the complete solution of the Riemann problem.

Case 1: \(h_* > h_L\) and \(h_* > h_R\), the solution consists of a 1-shock and a 2-shock:

\[\overrightarrow{w_p}(x,t) = \left\{ \begin{array}{ll} \overrightarrow{w_p}_L & \mbox{ if } x < \dot\sigma_1 t \\ \overrightarrow{w_p}_* & \mbox{ if } \dot\sigma_1 t < x < \dot\sigma_2 t \\ \overrightarrow{w_p}_R & \mbox{ if } x > \dot\sigma_2 t \end{array} \right.\]

Case 2: \(h_L < h_* < h_R\), the solution consists of a 1-shock and a 2-rarefaction: \[\overrightarrow{w_p}(x,t) = \left\{ \begin{array}{ll} \overrightarrow{w_p}_L & \mbox{ if } x < \dot\sigma_1 t \\ \overrightarrow{w_p}_* & \mbox{ if } \dot\sigma_1 t < x < \lambda_2^* t \\ \overrightarrow{w_p}_{2det}(x,t) & \mbox{ if } \lambda_2^* t < x < \lambda_2^R t \\ \overrightarrow{w_p}_R & \mbox{ if } x > \lambda_2^R t \end{array} \right.\]

Case 3: \(h_R < h_* < h_L\), the solution consists of a 1-rarefaction and a 2-shock: \[\overrightarrow{w_p}(x,t) = \left\{ \begin{array}{ll} \overrightarrow{w_p}_L & \mbox{ if } x < \lambda_1^L t \\ \overrightarrow{w_p}_{1det}(x,t) & \mbox{ if } \lambda_1^L t < x < \lambda_1^* t \\ \overrightarrow{w_p}_* & \mbox{ if } \lambda_1^* t < x < \dot\sigma_2 t \\ \overrightarrow{w_p}_R & \mbox{ if } x > \dot\sigma_2 t \end{array} \right.\]

Case 4: \(h_* < h_L\) and \(h_* < h_R\), the solution consists of a 1-rarefaction and a 2-rarefaction: \[\overrightarrow{w_p}(x,t) = \left\{ \begin{array}{ll} \overrightarrow{w_p}_L & \mbox{ if } x < \lambda_1^L t \\ \overrightarrow{w_p}_{1det}(x,t) & \mbox{ if } \lambda_1^L t < x < \lambda_1^* t \\ \overrightarrow{w_p}_* & \mbox{ if } \lambda_1^* t < x < \lambda_2^* t \\ \overrightarrow{w_p}_{2det}(x,t) & \mbox{ if } \lambda_2^* t < x < \lambda_2^R t\\ \overrightarrow{w_p}_R & \mbox{ if } x > \lambda_2^R t \end{array} \right.\]

Note

In the case where dry zones appear, the algorithm is adapted as explained, for example, in this course by Gloria Faccanoni or in the thesis of Pons (2018).

Note

This Riemann solver is used in CERF. It is implemented in Fortran90 in ~/sources/PHY/svriemann.F90.

Important 1.2

In practice, it is necessary to define a threshold value defining the “dry” zone. In CERF, we therefore define \(h_{crit}\) such that if \(h<h_{crit}\), the zone is considered dry and \(h =0\) and \(\overrightarrow{u} = \overrightarrow{0}\) are imposed. This value can be modified by the user.

1.3 Two-fluid flow model

1.3.1 Governing equations

This model was originally developed to simulate flows in marine hydrodynamics, wave propagation and wave breaking. For breaking, the most physically relevant approach is to consider the two fluids, air and water, with behavior governed by an incompressible Navier-Stokes model, with or without turbulence. Handling the interface is, moreover, a topic in its own right! Whatever discretization approach is used, these simulations are very costly in terms of computation time. Moreover, due to numerical dissipation, it is difficult to simulate flows at very large scales. This is why, to simulate wave propagation phenomena, models of the Saint-Venant, Serre-Green-Naghdi, or inviscid irrotational flow type are often preferred. These models, however, model breaking phenomena poorly, or not at all. The literature on these topics is very extensive, and we do not claim to provide a state of the art here. We will simply offer a partial comparative study in Helluy, Philippe et al. (2005).
The two-fluid model presented here stems from the desire to find a compromise between the relevance of Navier-Stokes-type models and the simplicity of Saint-Venant-type models. We therefore consider a mixture of two fluids, so that no interface “tracking” is required and no surface tension can be taken into account. We further assume, which holds true in most cases studied in hydrodynamics, that viscosity can be neglected, which, from a numerical standpoint, advantageously frees us from modeling the Laplacian operator. Finally, as in Chorin (1967), we consider that the two fluids are very weakly compressible. This frees us from the need to impose incompressibility, which is costly in computation time; however, it becomes necessary to add an equation of state for the mixture. After investigation, according to Golay and Helluy (2007), it appears that an isothermal equation of state is sufficient for most of the cases studied, so that we can also dispense with solving the energy equation. We thus obtain a Eulerian model of the form:

\[ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \overrightarrow{u}) = 0 \tag{1.14}\] \[ \frac{\partial \rho\overrightarrow{u}}{\partial t} + \nabla \cdot \left(\rho \overrightarrow{u} \otimes \overrightarrow{u} + p \overline{\overline{I}} \right) = \rho\overrightarrow{g} \tag{1.15}\]

where \(\rho\) is the density of the mixture, \(\overrightarrow{u}=<u,v,w>^T\) is the velocity of the mixture, \(\overrightarrow{g}\) is the gravity vector, and where \(\overline{\overline{I}}\) denotes the identity matrix, \(\nabla \cdot\) the divergence operator, and \(\otimes\) the tensor product. Equation 1.14 represents conservation of mass and Equation 1.15 represents conservation of momentum. \(p\) is the pressure of the mixture, determined by the equation of state:
\[ p = p_0 + c_0^2 \left(\rho - \left( \phi \rho_a + (1 - \phi) \rho_w \right) \right) \tag{1.16}\]
where \(\rho_a\) represents the density of air, \(\rho_w\) the density of water, \(p_0\) the reference pressure, and \(c_0\) the speed of sound in the mixture. \(\phi\) represents the volume fraction of air in the mixture. \(\phi\) satisfies a transport equation, which we take in non-conservative form (Golay and Helluy (2007)), as follows:
\[ \frac{\partial \phi}{\partial t} + \overrightarrow{u} \cdot \nabla \phi = 0 \tag{1.17}\]

Typically, we take \(p_0 = 10^5 Pa\), \(\rho_a = 1. kg/m^3\), \(\rho_w = 1000. kg/m^3\) and \(c_0 = 20. m/s\).
In hydrodynamics, flows are on the order of a few \(m/s\), so that a speed of sound on the order of \(20 m/s\) leads to a Mach number \(M = \displaystyle \frac{||\overrightarrow{u}||}{c_0} < 0.3\). The mixture can then be considered nearly incompressible. This speed of sound is entirely artificial and does not correspond to any physical reality. It is chosen for numerical reasons. Indeed, the choice of \(c_0\) influences the stability of the explicit numerical scheme used to solve the equations, through the CFL condition (Important 1.1).

However, as it happens, according to Wood’s equation (Figure 1.9), while the speed of sound in water is on the order of 1481 \(m/s\) and in air on the order of 343 \(m/s\), the speed of sound in an air-water mixture is on the order of 23 \(m/s\). This is in no way a justification for the model used here, but it is interesting to note that the speed of sound in an air-water mixture is of the same order as the artificial speed of sound used in the two-fluid model.

Figure 1.9: Speed of sound in an air-water mixture according to Wood’s equation

Equations Equation 1.14, Equation 1.15 and Equation 1.17, together with Equation 1.16, constitute the hyperbolic problem solved by CERF to model two-fluid flows. We define the conservative variables \(\overrightarrow{w}=<\rho,\rho u,\rho v,\rho w,\phi>^T\) and the primitive variables \(\overrightarrow{w_p}=<\rho, u, v, w, p>^T\). This corresponds to the general form Equation 1.1.

1.3.2 Refinement criterion

For the two-fluid model, to define the numerical entropy production Equation 1.4, we use the entropy defined by: \[ s(\overrightarrow{w}) = \frac{1}{2}\rho ||\overrightarrow{u}||^2 + c_0^2 \rho ln(\rho) - c_0^2 \left( \rho_w - \rho_a \right) \phi\] and the entropy flux \[ \overrightarrow{\psi}_s(\overrightarrow{w}) = \left( s + c_0^2 \rho + c_0^2 \left( \rho_w - \rho_a \right) \right) \overrightarrow{u} = \left( \frac{1}{2}\rho ||\overrightarrow{u}||^2 + c_0^2 \rho \left( ln(\rho) + 1 \right) \right) \overrightarrow{u} .\]

1.3.3 Free-surface sharpening

The finite-volume resolution scheme for the hyperbolic system Equation 1.14, Equation 1.15 and Equation 1.17, even at second order, is known to be dissipative. As a result, the air-water interface diffuses and widens. To counteract this phenomenon, an interface sharpening model can be used. The idea is to add a source term to the transport equation for \(\phi\) Equation 1.17, which penalizes mixing zones using a splitting method. Over a fictitious time step, we therefore compute: \[ \frac{\partial \phi}{\partial \tau} = \phi^2 (1 - \phi)^2(\phi - c)\] where \(c\) is determined so as to conserve the volume fraction, i.e.:
\[ c = \frac{\int _{\Omega} \phi^3 (1 - \phi)^2 d\Omega}{\int _{\Omega} \phi^2 (1 - \phi)^2 d\Omega}\]

Note that for \(\phi=0\) or \(\phi=1\), this source term has no effect. Finally, since \(\phi\) is modified, for the sake of conservation we also modify the density and momentum of the mixture:
\[ \frac{\partial \rho}{\partial \tau} = \phi^2 (1 - \phi)^2(\phi - c)(\rho_a - \rho_w)\] \[ \frac{\partial \left(\rho \overrightarrow{u} \right)}{\partial \tau} = \phi^2 (1 - \phi)^2(\phi - c) \overrightarrow{u}(\rho_a - \rho_w)\]

The fictitious time step is calibrated to represent a fraction of the flux. This fraction is determined by a sharpening coefficient \(\alpha_{sharp}\). In practice, this coefficient is on the order of \(10^{-3} - 10^{-2}\).

1.3.4 The “rigid fluid” concept

Following the fictitious domain principle, part of the cells can be penalized to simulate an obstacle by imposing a zero velocity. Following the work of Coquerelle and Cottet (2008), part of the cells can also be penalized to obtain “rigid solid”-type behavior. The extension of the work of Coquerelle and Cottet (2008) to a two-fluid model was carried out by Altazin (2017). Consider a “solid” occupying the domain \(\Omega_s\), of uniform density \(\rho_s\), center of gravity \(G\), mass \(M_s\), and inertia matrix \(\overline{\overline{J_s}}\). The velocity of any point \(P\) of the solid is given by:
\[ \overrightarrow{v}_s = \overrightarrow{v}_G + \overrightarrow{\omega}_G \wedge \overrightarrow{GP}\]
where \(\overrightarrow{v}_G\) is the velocity of the center of gravity, \(\overrightarrow{\omega}_G\) is the angular velocity of the solid, and \(\overrightarrow{GP}\) is the vector connecting the center of gravity to the point \(P\). By minimizing \(\int _{\Omega_s} || \rho\overrightarrow{u} - \rho_s\overrightarrow{v}_s||^2 d\Omega_s\), we obtain, for a uniform solid density: \[\overrightarrow{v}_G=\frac{\int _{\Omega_s} \rho\overrightarrow{u} d\Omega_s}{M_s}\]
and
\[\overrightarrow{\omega}_G = \overline{\overline{J_s}}^{-1} \int _{\Omega_s} \overrightarrow{GP} \wedge \overrightarrow{u} .\]
At the end of each time step, the position and velocity of the solid’s center of gravity and the solid’s rotational velocity are computed. \(\rho_s\) and \(\overrightarrow{v}_s\) are then assigned to all the cells of the “solid”. In CERF, the solid is determined by its envelope. This envelope is defined by a triangular mesh. The cells inside the envelope are determined using a “ray casting” technique. In CERF, a solid can be considered “fixed”, “with imposed velocity”, or “free”.

1.3.5 Riemann solver for the two-fluid model

1.3.5.1 Overview

Numerical fluxes in CERF are computed using a Riemann solver. There are many references concerning this solver, such as Toro (1997). This section describes the Riemann solver used for the two-fluid flow model, based on Golay and Helluy (2007).

We therefore seek to solve the isothermal two-fluid problem in one dimension. We consider a flow with density \(\rho(x,t)\), horizontal velocity \(u(x,t)\), pressure \(p(x,t)\), and volume fraction \(\phi(x,t)\), described by the following problem: \[ \left \{ \begin{array}{ll} \displaystyle\frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} = 0 \\ \displaystyle\frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u^2 + p)}{\partial x} = 0 \\ \displaystyle\frac{\partial \rho \phi}{\partial t} + \frac{\partial \rho \phi u}{\partial x} = 0 \\ \end{array} \right. \tag{1.18}\] With an equation of state of the form: \[ p = p_0 + c_0^2 \left(\rho - \left( \phi \rho_a + (1 - \phi) \rho_w \right) \right) \] and initial conditions such that: \[ \overrightarrow{w_p}(x,0) =\left\{ \begin{array}{ll} \overrightarrow{w_p}_L &\mbox{ if } x<0 \\ \overrightarrow{w_p}_R &\mbox{ if } x>0 \end{array} \right. \] where \(\overrightarrow{w_p}_L=<\rho_L, u_L, \phi_L>^T\) and \(\overrightarrow{w_p}_R=<\rho_R, u_R, \phi_R>^T\) are two constant states.

By expanding Equation 1.18, the system can be rewritten in the following form: \[ \frac{\partial}{\partial t} \left\{ \begin{array}{c} \rho \\ u \\ \phi \end{array} \right\} + \left[ \begin{array}{ccc} u & \rho & 0 \\ \displaystyle\frac{c_0^2}{\rho} & u & \displaystyle\frac{c_0^2\left( \rho_w - \rho_a \right)}{\rho} \\ 0 & 0 & u \end{array} \right] \frac{\partial}{\partial x} \left\{ \begin{array}{c} \rho \\ u \\ \phi \end{array} \right\} = \left\{ \begin{array}{c} 0 \\ 0 \\ 0\end{array} \right\} . \]

We then seek the eigenvalues \(\lambda_1\), \(\lambda_2\), and \(\lambda_3\) as well as the eigenvectors \(\overrightarrow{r}_1\), \(\overrightarrow{r}_2\), and \(\overrightarrow{r}_3\) of the matrix thus obtained. A quick calculation gives: \[ \lambda_1=u - c_0 , \lambda_2=u, \lambda_3=u + c_0 \] with \(\lambda_1 \leq \lambda_2 \leq \lambda_3\) and \[\overrightarrow{r}_1=\left\{ \begin{array}{c} -\rho/c_0 \\ 1 \\0 \end{array} \right\} , \overrightarrow{r}_2=\left\{ \begin{array}{c} \left( \rho_w - \rho_a \right) \\ 0 \\ -1 \end{array} \right\}, \overrightarrow{r}_3=\left\{ \begin{array}{c} \rho/c_0 \\ 1 \\0 \end{array} \right\} .\]

Note 1.2

The eigenvalues are used to determine the wave propagation speed within the cell Important 1.1 as: \(\lambda_k = Max\left( |u - c_0| , |u|, |u + c_0| \right)\)

To determine the nature of the three characteristic fields, we compute \(\nabla_{w_p}\lambda_1 \cdot \overrightarrow{r}_1 = 1 \neq 0\), \(\nabla_{w_p}\lambda_2 \cdot \overrightarrow{r}_2 = 0\) and \(\nabla_{w_p}\lambda_3 \cdot \overrightarrow{r}_3 = 1 \neq 0\). Characteristic fields 1 and 3 are therefore genuinely nonlinear, while field 2 is linearly degenerate. There will therefore be one contact discontinuity and two waves, which are either shocks or rarefaction waves.

Since there are \(3\) waves, \(2\) Riemann invariants can be defined for each wave. The Riemann invariants of wave \(k\) are such that \(\nabla_{w_p}I_{k;1} \cdot \overrightarrow{r}_k = 0\) and \(\nabla_{w_p}I_{k;2} \cdot \overrightarrow{r}_k = 0\). It is easy to verify that:
\(I_{1;1} = \phi\) , \(I_{1;2} = u + c_0 ln(\rho)\),
\(I_{2;1} = u\) , \(I_{2;2} = p\),
\(I_{3;1} = \phi\) , \(I_{3;2} = u - c_0 ln(\rho)\).

The entropy weak solution of the Riemann problem consists of at most four constant states \(\overrightarrow{w_p}_L\), \(\overrightarrow{w_p}_I\), \(\overrightarrow{w_p}_{II}\), \(\overrightarrow{w_p}_R\), separated by a shock or a rarefaction wave, a contact discontinuity, then a shock or a rarefaction wave (Figure 1.10). In what follows, we determine the states \(\overrightarrow{w_p}_I\) and \(\overrightarrow{w_p}_{II}.\)

Figure 1.10: Diagram of the solutions to the Riemann problem with the isothermal two-fluid model

1.3.5.2 Study of the rarefaction waves

Consider a state \(\overrightarrow{w_p}_g\) and a state \(\overrightarrow{w_p}_d\), connected by a k-rarefaction wave. The following conditions must then be satisfied:

  • \(\lambda_k(\overrightarrow{w_p}_g)<\lambda_k(\overrightarrow{w_p}_d)\)
  • The Riemann invariants are conserved along the rarefaction wave \(I_{k;i}(\overrightarrow{w_p}_g)=I_{k;i}(\overrightarrow{w_p}_d)\) for \(i=1,2\)

Let us apply these conditions to the isothermal two-fluid problem.

Consider first a 1-rarefaction connecting a state \(\overrightarrow{w_p}_L\) to a state \(\overrightarrow{w_p}_I\), separated by a rarefaction wave. The Riemann invariants give:
\[\left\{ \begin{array}{l}I_{1;1}(\overrightarrow{w_p}_L) = I_{1;1}(\overrightarrow{w_p}_I) \\ I_{1;2}(\overrightarrow{w_p}_L) = I_{1;2}(\overrightarrow{w_p}_I)\end{array} \right. \rightarrow \left\{ \begin{array}{l}\phi_L=\phi_I \\ u_L + c_0 ln(\rho_L) = u_I + c_0 ln(\rho_I)\end{array} \right.\] that is \[ u_I = u_L + c_0 ln\left( \frac{\rho_L}{\rho_I}\right) \tag{1.19}\] Furthermore, \[\lambda_1(\overrightarrow{w_p}_L)<\lambda_1(\overrightarrow{w_p}_I) \rightarrow u_L < u_I \] and consequently, from Equation 1.19, \(\rho_L > \rho_I\).

Consider finally a 3-rarefaction connecting a state \(\overrightarrow{w_p}_{II}\) to a state \(\overrightarrow{w_p}_R\), separated by a rarefaction wave. The condition on the velocities is written: \[\lambda_3(\overrightarrow{w_p}_{II})<\lambda_3(\overrightarrow{w_p}_R) \rightarrow u_{II} < u_R .\]

The Riemann invariants give \[\left\{ \begin{array}{l}I_{3;1}(\overrightarrow{w_p}_{II}) = I_{3;1}(\overrightarrow{w_p}_R) \\ I_{3;2}(\overrightarrow{w_p}_{II}) = I_{3;2}(\overrightarrow{w_p}_R)\end{array} \right. \rightarrow \left\{ \begin{array}{l}\phi_{II}=\phi_R \\ u_{II} - c_0 ln(\rho_{II}) = u_R - c_0 ln(\rho_R)\end{array} \right.\] that is \[ u_{II} = u_R - c_0 ln\left( \frac{\rho_R}{\rho_{II}}\right) \tag{1.20}\]
so since \(u_{II} < u_R\), we have \(\rho_{II} < \rho_R\).

1.3.5.3 Study of the shocks

Consider a state \(\overrightarrow{w_p}_g\) and a state \(\overrightarrow{w_p}_d\), connected by a \(k\) shock wave. The following conditions must then be satisfied:

  • The Rankine-Hugoniot conditions, where \(\dot\sigma_k\) denotes the shock speed: \(\displaystyle \dot\sigma_k=\frac{\left[\overrightarrow{F}(\overrightarrow{w})\right]}{\left[\overrightarrow{w}\right]}\)
  • The Lax entropy conditions for a genuinely nonlinear field \(k\): \[ \left\{ \begin{array}{l} \lambda_{k-1}(\overrightarrow{w_p}_g) < \dot\sigma_k < \lambda_{k}(\overrightarrow{w_p}_g) \\ \lambda_{k}(\overrightarrow{w_p}_d) < \dot\sigma_k < \lambda_{k+1}(\overrightarrow{w_p}_d) \end{array} \right. \]

Applied to the two-fluid problem, the Rankine-Hugoniot conditions give: \[ \dot\sigma_k = \frac{\rho_d u_d - \rho_g u_g}{\rho_d - \rho_g} = \frac{\rho_d u_d^2 + p_d - \rho_g u_g^2 - p_g}{\rho_d u_d - \rho_g u_g} = \frac{\rho_d u_d \phi_d - \rho_g u_g \phi_g}{\rho_d \phi_d - \rho_g \phi_g} \] \[ \left\{ \begin{array}{l} \rho_d \dot\sigma_k - \rho_g \dot\sigma_k = \rho_d u_d - \rho_g u_g \\ \rho_d u_d \dot\sigma_k - \rho_g u_g \dot\sigma_k = \rho_d u_d^2 + p_d - \rho_g u_g^2 - p_g \\ \rho_d \phi_d \dot\sigma_k - \rho_g \phi_g \dot\sigma_k = \rho_d u_d \phi_d - \rho_g u_g \phi_g\end{array} \right.\] \[ \left\{ \begin{array}{l} \rho_g (u_g - \dot\sigma_k) = \rho_d (u_d - \dot\sigma_k) \\ u_d \rho_d (u_d - \dot\sigma_k) -u_g \rho_g (u_g - \dot\sigma_k) + p_d - p_g = 0 \\ \phi_g \rho_g (u_g - \dot\sigma_k) = \phi_d \rho_d (u_d - \dot\sigma_k) \end{array} \right.\] So, denoting \(j_k = \rho_g (u_g - \dot\sigma_k) = \rho_d (u_d - \dot\sigma_k)\), we obtain: \[\phi_g = \phi_d\]
and
\[j_k = \frac{p_d - p_g}{u_g - u_d}\]

By eliminating the shock speed, we also have: \[\left( \rho_d u_d - \rho_g u_g \right)^2 = \left( \rho_d - \rho_g \right) \left( \rho_d u_d^2 - \rho_g u_g^2 + p_d - p_g \right)\] or equivalently
\[\rho_d \rho_g \left( u_d - u_g \right)^2 = \left( \rho_d - \rho_g \right) \left( p_d - p_g \right)\]

Consider first a 1-shock connecting a state \(\overrightarrow{w_p}_L\) to a state \(\overrightarrow{w_p}_I\), separated by a shock wave. The Lax admissibility conditions are written: \[ \displaystyle\left\{ \begin{array}{l} \dot\sigma_1 < \lambda_1(\overrightarrow{w_p}_L) \\ \lambda_1(\overrightarrow{w_p}_I) <\dot\sigma_1 < \lambda_2(\overrightarrow{w_p}_I) \end{array} \right. \rightarrow \left\{ \begin{array}{l} \dot\sigma_1 < u_L - c_0 \\ u_I - c_0 < \dot\sigma_1 \\ \dot\sigma_1 < u_I \end{array} \right. \rightarrow \left\{ \begin{array}{l} \rho_L c_0 < j_1 \\ j_1 < \rho_I c_0 \\ \dot\sigma_1 < u_I \end{array} \right.\]

So in particular \(u_I < u_L\) and \(\rho_L < \rho_I\). The Rankine-Hugoniot conditions give: \[ \displaystyle\left\{ \begin{array}{l} \phi_L = \phi_I \\ \left( u_I - u_L \right)^2 = \displaystyle\frac{\rho_I - \rho_L}{\rho_I \rho_L}\left( p_I - p_L \right) \end{array} \right.\] from which, since \(u_I < u_L\) and \(\rho_L < \rho_I\): \[ u_I=u_L -\sqrt{\frac{\rho_I - \rho_L}{\rho_I \rho_L}\left( p_I - p_L \right)} = u_L + c_0 \frac{\rho_L - \rho_I}{\sqrt{\rho_I \rho_L}} \tag{1.21}\]

Consider finally a 3-shock connecting a state \(\overrightarrow{w_p}_{II}\) to a state \(\overrightarrow{w_p}_R\), separated by a shock wave. The Lax conditions give: \[ \displaystyle\left\{ \begin{array}{l} \lambda_2(\overrightarrow{w_p}_{II}) < \dot\sigma_3 < \lambda_3(\overrightarrow{w_p}_{II}) \\ \lambda_3(\overrightarrow{w_p}_R) < \dot\sigma_3 \end{array} \right. \rightarrow \left\{ \begin{array}{l} u_{II} < \dot\sigma_3 \\ \dot\sigma_3 < u_{II} + c_0 \\ u_R + c_0 < \dot\sigma_3 \end{array} \right. \rightarrow \left\{ \begin{array}{l} u_{II} < \dot\sigma_3 \\ - c_0 \rho_{II} < j_3 \\ j_3 < - c_0 \rho_R \end{array} \right.\]

So in particular \(u_R < u_{II}\) and \(\rho_R < \rho_{II}\). The Rankine-Hugoniot conditions are written: \[ \displaystyle\left\{ \begin{array}{l} \phi_{II} = \phi_R \\ \left( u_R - u_{II} \right)^2 = \displaystyle\frac{\rho_R - \rho_{II}}{\rho_R \rho_{II}}\left( p_R - p_{II} \right) \end{array} \right.\] from which, since \(u_R < u_{II}\) and \(\rho_R < \rho_{II}\): \[ u_{II}=u_R + \sqrt{\frac{\rho_R - \rho_{II}}{\rho_R \rho_{II}}\left( p_R - p_{II} \right)} = u_R - c_0 \frac{\rho_R - \rho_{II}}{\sqrt{\rho_R \rho_{II}}} \tag{1.22}\]

1.3.5.4 Study of the contact discontinuity

Consider a state \(\overrightarrow{w_p}_g\) and a state \(\overrightarrow{w_p}_d\), connected by a contact discontinuity. Across a k-contact discontinuity, entropy and the Riemann invariants are conserved. Since the discontinuity moves at speed \(\dot\sigma_k\), we have: \[\displaystyle \dot\sigma_k = \lambda_{k}(\overrightarrow{w_p}_g) = \lambda_{k}(\overrightarrow{w_p}_d) = \frac{\left[\overrightarrow{F}(\overrightarrow{w})\right]}{\left[\overrightarrow{w}\right]}\]

Consider then a 2-contact discontinuity connecting the state \(\overrightarrow{w_p}_I\) to the state \(\overrightarrow{w_p}_{II}\). We then have \[ \dot\sigma_2 = u_I = u_{II} = \frac{\rho_{II} u_{II} - \rho_I u_I}{\rho_{II} - \rho_I} = \frac{\rho_{II} u_{II}^2 + p_{II} - \rho_I u_I^2 - p_I}{\rho_{II} u_{II} - \rho_I u_I} = \frac{\rho_{II} u_{II} \phi_{II} - \rho_I u_I \phi_I}{\rho_{II} \phi_{II} - \rho_I \phi_I} \]

From which we deduce: \[ u_I = u_{II} \mbox{ and } p_I = p_{II} \tag{1.23}\]

1.3.5.5 Constructing the solution

To simplify notation, let us define a function \(H(a,b)\) such that \(H(a,b) = \displaystyle\left\{ \begin{array}{ll} c_0 \displaystyle\frac{\rho_a - \rho_b}{\sqrt{\rho_a \rho_b}} & \mbox{ if } \rho_a < \rho_b \\ c_0 ln \left(\displaystyle\frac{\rho_a}{\rho_b} \right) & \mbox{ if } \rho_a > \rho_b \end{array} \right.\).
For the first wave, using Equation 1.19, Equation 1.21, we can then write \(u_I = u_L + H(\rho_L, \rho_I)\), and for the third wave, using Equation 1.20, Equation 1.22, \(u_{II} = u_R - H(\rho_R, \rho_{II})\).

We know from the study of the contact discontinuity that: \[u_I - u_{II} = 0\] so
\[u_L + H(\rho_L, \rho_I) - u_R + H(\rho_R, \rho_{II}) = 0\] Now, \(\rho_I=\rho_I(p_I,\phi_I)\) and \(\rho_{II}=\rho_{II}(p_{II},\phi_{II})\). We also know that \(p_I = p_{II}\); we will therefore denote the intermediate pressure \(p^*\) such that \(p^* = p_I = p_{II}\). Moreover, we know \(\phi_I = \phi_L\) and \(\phi_{II} = \phi_R\). We can then write: \[u_L - u_R + H(\rho_L, \rho(p^*,\phi_L)) + H(\rho_R, \rho(p^*,\phi_R)) = 0\] This is therefore a nonlinear equation in \(p^*\) that can be solved using a Newton-Raphson method. Once \(p^*\) is known, the states \(\overrightarrow{w_p}_I\) and \(\overrightarrow{w_p}_{II}\) can be determined.

Note

This Riemann solver is used in CERF. It is implemented in Fortran90 in ~/sources/PHY/riemisot.F90.