# Augmented Skew-Symetric System for Shallow-Water System with Surface Tension Allowing Large Gradient of Density

In this paper, we introduce a new extended version of the shallow water equations with surface tension which is skew-symmetric with respect to the L2 scalar product and allows for large gradients of fluid height. This result is a generalization of the results published by P. Noble and J.-P. Vila in [SIAM J. Num. Anal. (2016)] and by D. Bresch, F. Couderc, P. Noble and J.P. Vila in [C.R. Acad. Sciences Paris (2016)] which are restricted to quadratic forms of the capillary energy respectively in the one dimensional and two dimensional setting.This is also an improvement of the results by J. Lallement, P. Villedieu et al. published in [AIAA Aviation Forum 2018] where the augmented version is not skew-symetric with respect to the L2 scalar product. Based on this new formulation, we propose a new numerical scheme and perform a nonlinear stability analysis.Various numerical simulations of the shallow water equations are presented to show differences between quadratic (w.r.t the gradient of the height) and general surface tension energy when high gradients of the fluid height occur.

There are no comments yet.

## Authors

• 1 publication
• 2 publications
• 1 publication
• 1 publication
• 3 publications
• 22 publications
• 1 publication
• 2 publications
06/26/2019

### A decoupled staggered scheme for the shallow water equations

We present a first order scheme based on a staggered grid for the shallo...
12/11/2020

### Stochastic dynamics of storm surge with stable noise

The Advanced Circulation (ADCIRC) and Simulating Nearshore Waves (SWAN) ...
06/05/2020

### Ensuring 'well-balanced' shallow water flows via a discontinuous Galerkin finite element method: issues at lowest order

The discontinuous Galerkin finite element method (DGFEM) developed by Rh...
10/26/2021

### An Arbitrary High Order and Positivity Preserving Method for the Shallow Water Equations

In this paper, we develop and present an arbitrary high order well-balan...
09/17/2019

### Hamiltonian regularisation of shallow water equations with uneven bottom

The regularisation of nonlinear hyperbolic conservation laws has been a ...
03/28/2021

### Application of Graphics Processing Units for self-consistent modelling of shallow water dynamics and sediment transport

In this paper, we describe a numerical algorithm for the self-consistent...
05/09/2018

### Presburger Arithmetic with algebraic scalar multiplications

We study complexity of integer sentences in S_α = (R, <, +,Z, x α x), wh...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In this paper, we consider shallow-water type equations with a full surface tension term issued from Hamiltonian formulation of P. Casal and H. Gouin ([4]) (see also D. Serre ([13])):

 {∂th+div(hu)=0(i)∂t(hu)+div(hu⊗u)+∇P=−div(∇h⊗∇pE)+∇(hdiv(∇pE))(ii) (1)

with the fluid height,

the fluid velocity vector field. The internal energy

is defined by

 (2)

with and . Finally, the pressure is given by

 P(h,p):=h∂hE(h,p)−E(h,p)=π(h)−(κ(h)−hκ′(h))Ecap(∥p∥), (3)

where

The three given functions , and are assumed to be positive and invertible from to with . Moreover, we suppose so that is strictly convex as soon as . In this context, the system (1) admits an additional energy conservation law which reads

 ∂t(12h∥u∥2+E)+div((12h∥u∥2+E+P)u)−div(div(∇pE)hu)+div(div(hu)∇pE)=0. (4)

For specific choices of the capillary energy, we note that the system (1) reduces to classical model of the fluid mechanics literature like the Euler-Korteweg isothermal system when :

 E(h,∇h)=Φ(h)+12κ(h)∥∇h∥2

or the shallow water type system for thin film flows both in quadratic capillary case

 E(h,∇h)=Φ(h)+12σ∥∇h∥2

and the fully nonlinear capillary case :

 E(h,∇h)=Φ(h)+σ(√1+∥∇h∥2−1).

Note that the fully nonlinear case admits the two following asymptotics

 E(h,∇h)=Φ(h)+σ∥∇h∥22+o∥∇h∥→0(∥∇h∥2),E(h,∇h)=Φ(h)+σ∥∇h∥+o∥∇h∥→∞(∥∇h∥).

It is a hard problem to propose a discretization of the shallow water equations (1) that is compatible with the energy equation (4

). The main issue is that one cannot adapt the proof of the energy estimate (

4) derived from (1) at a discretized level due to the presence of high order derivatives. The readers interested to understand the mathematical and numerical difficulties are referred to [10] and important references cited therein. The strategy then consists in performing a reduction of order in spatial derivatives and introducing an alternative system which contains lower order derivatives. This strategy was applied successfully in the context of Euler-Korteweg isothermal system when the internal energy is quadratic with respect to : see [9] in the one dimensional case and [1] in the two dimensional case. In both cases, the augmented version is obtained by introducing an auxiliary velocity which is proportional to and it admits an additional skew-symmetric structure with respect to the scalar product which makes the proof of energy estimates and the design of compatible numerical scheme easier. However, this approach does not work in the context of general internal energy. In [7], the authors consider the following augmented version in order to deal with more general internal energies:

 (5)

where . However, in the -dimensional setting, the assumption has to be made to show the conservation of the total energy and therefore it has to be satisfied initially: The interested reader is referred page 166–168. In order to avoid such a constraint which is hardly guaranteed in the discrete case, one could use instead the following modified formulation

for which it is easy to prove the conservation of the total energy

 ∂t(Etot)+div(u(Etot+π))=(div(h(ut∇)(∇pEtot))−div(h(∇pEttot∇)u))−div(u(pt∇pEtot−(κ−hκ′)Ecap))

for any smooth solution of the above system. However, this formulation introduces nonconservative terms in the momentum equation and it is then hard to satisfy for conservation of momentum and energy.

In this paper, a new skew-symmetric augmented version which is a second order differential system is proved for (1)–(3). In the small gradient limit, this fomulation is equivalent to the one derived by D. Bresch, F. Couderc, P. Noble and J.–P. Vila in [1]. This formulation is valid for any internal energy in the form . When specified to , we see that in the high gradient limit, which is a capillary term found usually in two fluids systems. We thus expect our approach to be useful in the context of bi-fluid flows. Note also that our paper could be also of practical interest to deal with generalization of Euler-Korteweg system: see [6] and [5] for discussions on compressible Korteweg type systems. We rely on the new augmented system to propose a numerical scheme which is energetically stable and extends what was done in [1] and [9]. Note that skew-symmetric augmented versions of the capillary shallow water equations are also useful from a theoretical point of view: see e.g. [3] for the proof of existence of dissipative solutions to the Euler-Korteweg isothermal system. Our present work will be the starting point to improve the work by Lallement and Villedieu (see [7] and [8]) related to disjunction term for triple point simulations: see [2].

The paper is divided in three parts: The first part introduces the augmented version with full surface tension and discuss its connection with the system derived in [1]. In the second part, we propose a numerical scheme satisfying energy stability. Finally, we present numerical illustrations based on our numerical scheme showing the importance to consider the augmented system with the full surface tension.

Notations. Throughout this paper, we will write vectors in column forms and the transpose of a matrix or a vector is defined as: . The symbol will denote the gradient operator considered as a vector: . When applied respectively to a vector or a matrix, the divergence operator is defined as

 div(u)=∂1u1+∂2u2,div(A)i=∂1Ai1+∂2Ai2.

## 2 Augmented version

Extending ideas from [9] in the one dimensional case, an augmented formulation of the shallow water equations (1) with was proposed in [1]: it is a second order system of PDEs which is skew symmetric with respect to the scalar product. The additional quantity was given by with : it is thus colinear to . An alternative extended form was proposed in [7], [8] to deal with the general case by introducing the additional unknown : although it seems to be efficient in the one dimensional case, this approach hardly extends to the two dimensional setting due to a lack of energy consistency as soon as the condition is not satisfied.

We now introduce our new formulation of (1) which is valid in the fully decoupled case and provides a dual formulation of capillary terms which ensures a straightforward consistent energy balance. To this end we introduce an additional unknown, denoted , which is colinear to and satisfies

 12h∥v∥2=κ(h)Ecap(∥∇h∥)

where . To do so, we define v as

 v=α(q2)p√κ(h)h (6)

where is given by

 α(q2)=√2Ecap(q)q.

Note that using the definition , we have the following relations

 ∥v∥2=α2(∥p∥2)∥p∥2κh,12α2(q2)q2κ(h)=κ(h)Ecap(q).

Note that, in this context, has the dimension of velocity and transforms the capillary energy into some kinetic energy. This interpretation of the capillary energy in terms of kinetic energy in our augmented system defined below motivates surely the robustness of our results.

Let us now write a system related to the unknowns where is given by (6) with . This will give a first order hyperbolic system on together with a second order part which has a skew-symmetric structure (for the scalar product). More precisely, we have the following result.

###### Lemma 2.1

i) Let

 U=⎛⎜⎝hhuhv⎞⎟⎠,F(U)=⎛⎜⎝chuhu⊗u+π(h)Idhv⊗u⎞⎟⎠ (7)

and

 M=⎛⎜ ⎜⎝0div(h∇(f(h,v)v)t)−∇(g(h,v)tv)−f(h,v)div(h∇ut)−g(h,v)divu⎞⎟ ⎟⎠. (8)

where

is a symmetric tensor and

a vector field given by

 f(h,v)=√κ(h)√h(2α′(q2)hα(q2)2κ(h)v⊗v+α(q2)Id)
 g(h,v)=((κ′(h)h2κ(h)+12)+2α′(q2)α(q2)q2)hv

where

 α(q2)=√2Ecap(q)qwithq=E−1cap(h∥v∥22κ(h)).

The augmented system reads

 ∂tU+div(F(U))=M. (9)

ii) If is regular enough then it also satisfies the following energy balance

 ∂t(12h∥u∥2+E) +div(u(12h∥u∥2+E+π)) =(div(hut∇t(f(h,v)tv))−div(h∇uf(h,v)v))−div(u(g(h,v)tv))

where

iii) If is regular enough and the initial velocity satisfies

 v0=α(∥∇h0∥2)√κ(h0)h0∇h0

then satisfies also

 v=α(∥∇h∥2)√κ(h)h∇h

and solves the original Equations with the full surface tension term given by (1)–(3).

Proof of Lemma 2.1.

Part i) and iii) Equation satisfied by . Let us first recall that and therefore

 hv=α(q2)√κ(h)h3/2∇hh:=α(q2)F(h)a,

with and . In order to write an evolution equation on , the first step is to calculate evolution equations on , and . For that purpose, we consider the mass conservation law written as

 ∂th+ut∇h+hdiv(u)=0. (11)

By dividing (11) by and differentiating with respect to , one finds

 ∂ta+∇(uta)+∇(div(u))=0. (12)

By multiplying (11) by , one finds

 ∂tF(h)+ut∇F(h)=−hF′(h)div(u),F′(h)=κ′(h)h3/22√κ(h)+32√κ(h)h. (13)

By differentiating directly (11), we have

 ∂t∇h+(ut∇)∇h=−div(h∇ut)−∇hdivu. (14)

The derivatives of are given by

 q∂tq=(∇h)t∂t∇h,q∂iq=(∇h)t∂i∇h,i=1,2.

This allows to calculate the equation on . Indeed, we can write:

 ∂tα(q2)+ut∇α(q2)=α′(q2)(∂tq2+(ut∇)q2)=2α′(q2)(∇h)t(∂t∇h+(ut∇)∇h)

By substituting the value of given by (14) into the former equation, one finds

 ∂tα(q2)+ut∇α(q2)=−2α′(q2)((∇h)tdiv(h∇ut)+q2divu). (15)

Finally, by using the fact that , one finds that the advective term is given by

 div(α(q2)F(h)a⊗u)=α(q2)F(h)((ut∇)a+div(u)a)+((ut∇)(α(q2)F(h)))a.

We can now calculate the equation satisfied by using formula (12)–(15). More precisely we have

 ∂t(hv)+div(hv⊗u) = ∂t(α(q2)F(h)a)+div(α(q2)F(h)a⊗u) = α(q2)F(h)((∂t+ut∇)a+div(u)a)+((∂t+ut∇)(α(q2)F(h)))a = α(q2)F(h)((ut∇)a+div(u)a−∇(uta+div(u))) −(α(q2)hF′(h)div(u)+2F(h)α′(q2)((∇h)tdiv(h∇ut)+div(u)q2))a = α(q2)F(h)((ut∇)a−∇(uta+div(u))) −((hF′(h)F(h)−1)div(u)+2α′(q2)α(q2)((∇h)tdiv(h∇ut)+div(u)q2))hv.

Note that we have the relation

 α(q2)F(h)((ut∇)a−∇(uta+div(u)))=−α(q2)F(h)hdiv(h∇ut).

and therefore, by using the relation , one finds

 ∂t(hv)+div(hv⊗u) = −div(u)(hF′(h)F(h)−1+2α′(q2)q2α(q2))hv −2α′(q2)α(q2)(h2vα(q2)F(h))tdiv(h∇ut)hv−α(q2)F(h)hdiv(h∇ut). = −div(u)(hF′(h)F(h)−1+2α′(q2)q2α(q2))hv −(2α′(q2)α(q2)2h3F(h)v⊗v+α(q2)F(h)hId)div(h∇ut).

This yields the conclusion on the evolution of .

Equation satisfied by . Let us first note that

 p=vα(q2)√κ(h)h

and therefore

 ∇pE=κ(h)(α(q2)2+2α(q2)α′(q2)∥p∥2)p=√κ(h)√h(α(q2)+2α′(q2)q2)v

Next, we expand and . First, one has

 f(h,v)v=√κ(h)√h(2α′(q2)hα(q2)2κ(h)∥v∥2+α(q2))v=√κ(h)√h(2α′q2+α)v=∇pE.

Now we observe that

 pt∇pE=(2α′(q2)q2+α(q2))α(q2)q2κ(h)

This yields

 g(h,v)tv =((κ′(h)h2κ(h)+12)+2α′(q2)α(q2)q2)h∥v∥2 =((κ′(h)h2κ(h)+12)α(q2)+2α′(q2)q2)α(q2)q2κ(h) =(2α′(q2)q2+α(q2))α(q2)q2κ(h)−(1−κ′(h)hκ(h))12(α(q2))2q2κ(h)

and thus

 g(h,v)tv=pt∇pE−(κ(h)−hκ′(h))Ecap(q)

Note that the momentum conservation equation of (1) can be written as:

 ∂t(hu)+div(hu⊗u)+∇π=−div(∇h⊗∇pE)+∇(hdiv(∇pE))+∇((κ(h)−hκ′(h))Ecap(∥∇h∥))

We now remark that

Then, by taking , we obtain

 div(h∇(∇pE)t)−∇(pt∇pE)=div(h∇(∇pE)t)−∇((∇h)t∇pE)=−div(∇h⊗∇pE)+∇(hdiv(∇pE))

and consequently the right-hand side of the momentum equation in the augmented system is :

 div(h∇(f(h,v)v)t)−∇(g(h,v)tv)=−div(∇h⊗∇pE)+∇(hdiv(∇pE))+∇((κ(h)−hκ′(h))Ecap(∥∇h∥))

and the momentum equation in the original system is satisfied, which gives the conclusion on for i).

Note that if is regular enough and the initial velocity satisfies

 v0=α(∥∇h0∥2)√κ(h0)h0∇h0

then satisfies also (6) and solves the original system.

Part ii). Recall that

 Etot(U)=12h(∥hu∥2+∥hv∥2)+Φ(h)

where is given by (7) and

 (Φh)′=πh2.

Let us consider the augmented system written as

 ∂tU+div(F(U))=M (16)

with the first order part given by

 U=⎛⎜⎝hhuhv⎞⎟⎠,F(U)=⎛⎜⎝chuhu⊗u+π(h)Idhv⊗u⎞⎟⎠

whereas the capillary term on the right hand side of (16) is given by

 M=⎛⎜ ⎜⎝0div(h∇(f(h,v)v)t)−∇(g(h,v)tv)−f(h,v)div(h∇ut)−g(h,v)divu⎞⎟ ⎟⎠.

The entropy variables for the first order part of (16) are given by

 Vt=(∇UEtot)t=(−12(∥u∥2+∥v∥2)+Φ′(h),ut,vt)

The energy equation is thus

 ∂tEtot+div(u(Etot+π))= (∇UEtot)tM = utdiv(h∇(f(h,v)v)t)−ut∇(g(h,v)tv) −vtf(h,v)div(h∇ut)−vtg(h,v)div(u) = utdiv(h∇(f(h,v)v)t)−(f(h,v)v)tdiv(h∇ut) −div(ug(h,v)tv) = div(h(ut∇)(f(h,v)v))−div(h((f(h,v)v)t∇)u)−div(ug(h,v)tv)

Recall that and and therefore

 ∂t(Etot)+div(u(Etot+π))=(div(h(ut∇)(∇pEtot))−div(h(∇pEttot∇)u))−div(u(pt∇pEtot−(κ−hκ′)Ecap)).

By chosing , we easily verify that

 (div(h(ut∇)(∇pEtot))−div(h(∇pEttot∇)u))−div(upt∇pEtot)=(div(h(ut∇)(∇pEtot))−div(h(∇pEttot∇)u))−div(u(∇h)t∇pEtot)=div(hdiv(∇pEtot)u)−div(div(hu)∇pEtot).

Then we get

 ∂t(Etot)+div(u(Etot+π−(κ−hκ′)Ecap))=div(hdiv(∇pEtot)u)−div(div(hu)∇pEtot)

which is exactly the formulation (4) of the Energy balance of the original system.

## 3 Energetically stable numerical scheme

The augmented formulation in Lemma 2.1 reads

 ∂tU+div(F(U))=M

with definitions (7,8) of and . The first order part of the augmented formulation in Lemma 2.1 is the classical Euler barotropic model with additional transport and admits an additional conservation law related to the total energy:

 Etot=∥hu∥22h+Φ(h)+∥hv∥22h.

whereas the entropy variables are

 (∇UEtot)t=Vt=(−12(∥u∥2+∥v∥2)+Φ′(h),ut,vt).

This total energy is the total energy of the shallow water equations with surface tension whereas the potential energy associated to surface tension is transformed into kinetic energy associated to the artificial velocity . The full system admits also an energy equation:

 ∂tEtot+div(u(Etot+π(h))) = VtM = div(hut∇(∇pE))−div(h(∇pE)t∇u) −div(u(pt∇pE−(κ−hκ′)Ecap(q)))

with the right-hand side in conservation form. One of the aim of this paper is to design a numerical scheme that is free from a CFL condition associated to surface tension. For that purpose, we follow the strategy in [1] and introduce an IMplicit-EXplicit strategy where the hyperbolic step is explicit in time whereas the step associated to surface tension is implicit in time. The spatial discretization is based on an entropy dissipative scheme for the first order part whereas we mimic the skew symmetric structure found at the continuous level to discretize the right hand side . We prove that this strategy is energetically stable in the case of periodic boundary conditions.

### 3.1 IMplicit - EXplicit formulation

Following [1], we consider the following IMplicit-EXplicit time discretization: the hyperbolic step is explicit

 Un+1/2−UnΔt+div(F(Un))=0 (17)

and the capillary skew symmetric second order step

 Un+1−Un+1/2Δt=Mn+1 (18)

with

 Mn+1=⎛⎜ ⎜⎝0div(hn+1∇(f(hn+1,vn+1/2)vn+1)t)−∇(g(hn+1,vn+1/2)tvn+1)−f(hn+1,vn+1/2)div(hn+1∇(un+1)t)−g(hn+1,vn+1/2)divun+1⎞⎟ ⎟⎠.

The second step is not fully implicit: instead it is semi-implicit so that the problem to solve for is linear. This does not affect the order of the time discretization since the time splitting is already first order in time. Let us now consider the spatial discretization. We will use a generic Finite Volume context. We introduce a spatial discretization of and operators with finite volume methods. For that purpose, we denote a cell of the mesh , an edge of and a neighboring cell of : see figure 1 for an illustration.

We use a classical entropy satisfying scheme of numerical flux

 Gne,K=G(UnK,UnKe,ne,K)

where is the outward normal to the cell (of measure ) at the edge (of measure ). We denote the average of the vector on the cell K. The hyperbolic step then reads

 Un+1/2K=UnK−ΔtmK∑e∈∂KmeGne,K (19)

and we assume that it is entropy dissipative in the sense that it satisfies the following discrete Entropy inequality

 Etot(Un+1/2K)≤Etot(UnK)−ΔtmK∑e∈∂KmeHne,K (20)

where is the entropy numerical flux associated with In the particular case of Euler Barotropic equations such numerical schemes exist and satisfy this inequality provided a hyperbolic CFL condition of the type

 maxKΔtmKme∥∥∇UF(UnK)∥∥

is satisfied for some . Moreover, under a similar CFL condition, the positivity of the fluid is preserved and the total energy is strictly convex: this will be a useful property to prove entropy stability for numerical schemes. The second step is

 Un+1K=Un+1/2K+ΔtMn+1K (22)

with

 Mn+1K=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝0−∇3,Δ(g(hn+1K,vn+1/2K)tvn+1K)+div1,Δ(hn+1K∇1,Δ(f(hn+1K,vn+1/2K)vn+1K)T)−g(hn+1K,vn+1/2K)div3,Δ(un+1K)−f(hn+1K,vn+1/2K)div1,Δ(hn+1K∇1,Δ(un+1K)T)⎞⎟ ⎟ ⎟ ⎟ ⎟⎠ (23)

where , , , are linear discrete operators approximating the corresponding ones in the definition of the operator and that will be defined hereafter. In particular shall be chosen as the dual discrete operator of