with the fluid height,
the fluid velocity vector field. The internal energyis defined by
with and . Finally, the pressure is given by
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
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 :
or the shallow water type system for thin film flows both in quadratic capillary case
and the fully nonlinear capillary case :
Note that the fully nonlinear case admits the two following asymptotics
). 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  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  in the one dimensional case and  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 , the authors consider the following augmented version in order to deal with more general internal energies:
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
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 . 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  and  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  and . Note that skew-symmetric augmented versions of the capillary shallow water equations are also useful from a theoretical point of view: see e.g.  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  and ) related to disjunction term for triple point simulations: see .
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 . 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
2 Augmented version
Extending ideas from  in the one dimensional case, an augmented formulation of the shallow water equations (1) with was proposed in : 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 ,  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
where . To do so, we define v as
where is given by
Note that using the definition , we have the following relations
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.
where is a symmetric tensor and
is a symmetric tensor anda vector field given by
The augmented system reads
ii) If is regular enough then it also satisfies the following energy balance
Proof of Lemma 2.1.
Part i) and iii) Equation satisfied by . Let us first recall that and therefore
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
By dividing (11) by and differentiating with respect to , one finds
By multiplying (11) by , one finds
By differentiating directly (11), we have
The derivatives of are given by
This allows to calculate the equation on . Indeed, we can write:
By substituting the value of given by (14) into the former equation, one finds
Finally, by using the fact that , one finds that the advective term is given by
Note that we have the relation
and therefore, by using the relation , one finds
This yields the conclusion on the evolution of .
Equation satisfied by . Let us first note that
Next, we expand and . First, one has
Now we observe that
Note that the momentum conservation equation of (1) can be written as:
We now remark that
Then, by taking , we obtain
and consequently the right-hand side of the momentum equation in the augmented system is :
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
then satisfies also (6) and solves the original system.
Part ii). Recall that
where is given by (7) and
Let us consider the augmented system written as
with the first order part given by
whereas the capillary term on the right hand side of (16) is given by
The entropy variables for the first order part of (16) are given by
The energy equation is thus
Recall that and and therefore
By chosing , we easily verify that
Then we get
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
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:
whereas the entropy variables are
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:
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  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 , we consider the following IMplicit-EXplicit time discretization: the hyperbolic step is explicit
and the capillary skew symmetric second order step
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
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
and we assume that it is entropy dissipative in the sense that it satisfies the following discrete Entropy inequality
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
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
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