Variational discretization of non-isothermal compressible flow

# A mixed variational discretization for non-isothermal compressible flow in pipelines

## Abstract.

We consider the non-isothermal flow of a compressible fluid through pipes. Starting from the full set of Euler equations, we propose a variational characterization of solutions that encodes the conservation of mass, energy, and entropy in a very direct manner. This variational principle is suitable for a conforming Galerkin approximation in space which automatically inherits the basic physical conservation laws. Three different spaces are used for approximation of density, mass flux, and temperature, and we consider a mixed finite element method as one possible choice of suitable approximation spaces. We also investigate the subsequent discretization in time by a problem adapted implicit time stepping scheme for which exact conservation of mass as well as a slight dissipation of energy and increase of entropy are proven which are due to the numerical dissipation of the implicit time discretization. The main arguments of our analysis are rather general and allow us to extend the approach with minor modification to more general boundary conditions and flow models taking into account friction, viscosity, heat conduction, and heat exchange with the surrounding medium.

Keywords: compressible flow, gas transport, variational principle, Galerkin approximation, energy estimates, mixed finite elements, implicit time stepping

AMS-classification (2000): 35D30,35R02,37L65,76M10,76N99

## 1. Introduction

We consider the systematic numerical approximation of compressible flow in pipelines. Such problems arise for instance in the simulation and optimization of gas networks [4, 25]. Motivated by such applications we tacitly assume throughout the text that the flow is subsonic. For ease of presentation, we neglect for the moment the presence of friction, viscosity, and heat transfer, and therefore consider the Euler equations

 ∂tρ+∂xm =0, (1.1) ∂tm+∂x(m2ρ+p) =0, (1.2) ∂tE+∂x(mρ(E+p)) =0. (1.3)

More general flow models will be considered below. Here denotes the density, the mass flux, the pressure, and the total energy; further is the kinetic energy, and the internal energy. The above equations are assumed to hold on a bounded closed interval representing the pipeline and for all . We further assume for the moment that the pipe is closed, i.e.,

 m =0on the boundary. (1.4)

It is well-known that for smooth solutions of the Euler equations (1.1)–(1.3), one can also deduce a conservation law for the entropy [13, 24], namely

 ∂t(ρs)+∂x(ms) =0. (1.5)

One may replace one of the equations in (1.1)–(1.3) by the entropy equation (1.5) and thus obtain an equivalent system for the unknown density , mass flux , and temperature . This viewpoint will play an important role in our considerations below.

The system (1.1)–(1.5) is complemented by equations of state. Here we choose to characterize the pressure , the specific internal energy , and the specific entropy as functions of density and temperature , i.e.,

 p=p(ρ,θ),e=e(ρ,θ),ands=s(ρ,θ). (1.6)

In order to comply to the basic laws of thermodynamics, these functions have to satisfy certain compatibility conditions; see [2, 13] and Section 2 below.

Under appropriate assumptions on the initial data and the constitutive relations, the local existence of smooth solutions to the system (1.1)–(1.6) can be guaranteed [13, 24]. In the context of gas pipelines, solutions are expected to remain smooth for all time, which can be explained by the stabilizing effect of friction [21, 22]. Moreover, the flow takes place at low Mach number and therefore no shocks should be generated.

###### Notation 1.1.

A triple is called smooth positive solution of the system (1.1)–(1.6) above, if all equations hold in a point wise sense and .

Due to the many important applications, a vast amount of literature has been devoted to the study of numerical methods for compressible flow problems, and the one dimensional problem discussed above is typically used as a starting point. Despite that fact, the convergence analysis for the flow of inviscid fluids is not completely settled, not even in one space dimension. Finite volume methods are probably most widely used for the numerical approximation of compressible flow. While a rather complete convergence theory has been established for scalar conservation laws, only partial results are available concerning the analysis of finite volume methods for the Euler equations; we refer to [18, 19] and the references given there. Similar results hold for discontinuous Galerkin methods [5, 6] which can be understood as high-order generalization of finite volume schemes. Using the stabilizing effect of viscosity, some truly implementable numerical schemes for the isothermal compressible Navier-Stokes equations have been shown to be globally convergent to weak solutions [28, 29, 30]. These are however formulated in Lagrangian coordinates which is prohibitive for a possible extension to pipe networks. In [14, 17], a globally convergent non-conforming finite element method for the isothermal compressible Navier-Stokes equations in Eulerian coordinates has been proposed and analyzed. This method is based on the approximation of the velocity and density field and makes use of several stabilization terms which leads to a rather strong violation of the conservation laws. Moreover, the method degenerates in the inviscid limit; see [17, Sec. 3.1].

The main goal of the current manuscript is to construct a numerical scheme that can handle general flow models for viscous and inviscid fluids. Moreover, the method should have good stability properties and preserve the basic conservation laws that encode the underlying the physical principles as good as possible. For this purpose, we consider here an extension of our previous work [10], which dealt with the conforming Galerkin approximation of isentropic flow on networks. The extension to the full set of Euler equations and generalization thereof, which are the subject of this paper, will require some non-trivial extensions. In contrast to the isentropic setting considered in [10], we are not yet able to systematically handle pipe networks and such extension are therefore left as a topic for future research. One obstacle here is the formulation of appropriate coupling conditions at pipe junctions. In contrast to the isentropic case [10, 23, 23], this issue seems not completely settled in the general case; we refer to [1, 3, 7, 15, 16] for positive as well as negative examples.

The remainder of the manuscript is organized as follows: In the first part of the paper, we discuss in detail the Euler equations on a closed pipe. We first recall some basic relations of equilibrium thermodynamics, and then introduce a variational characterization of smooth solutions which rather directly encodes the conservation of mass, energy, and entropy. In the second part, we investigate the numerical approximation of this variational principle in space by a conforming Galerkin method and then consider the subsequent discretization in time by a problem adapted implicit time stepping scheme. We prove strict conservation of mass, energy, and entropy for the semi-discretization and a slight dissipation of energy and monotonic increase in entropy for the fully discrete scheme. The proposed method therefore perfectly complies to the basic principles of thermodynamics. In the third part of the paper, we investigate the extension of our approach to more general flow models. The last part of the paper is devoted to numerical tests which illustrate the theoretical results and demonstrate the stability and conservation properties of the fully discrete schemes. Our research is motivated mainly by gas transport in pipelines, which takes place at low Mach number and therefore avoids the generation of shocks. We however also consider a shock tube problem to demonstrate the correct handling of shocks, rarefaction waves, and contact discontinuities that may arise in more general applications.

## Part I: Analysis on the continuous level

In the following three sections, we first review some basic relations of equilibrium thermodynamics and then present and analyze a particular variational principle for the Euler equations which will serve as the basis for our further considerations.

## 2. Auxiliary results

Let us first consider in a bit more detail the relations of pressure, internal energy, and entropy. It is well known, see e.g. [2, 8, 13], that for thermodynamical consistency, the pressure has to be related to the specific internal energy by

 eρ=1ρ2(p−θpθ). (2.1)

Subscripts denote partial derivatives and functions may in general depend on and . By integration with respect to , we can then express the internal energy as

 e(ρ,θ)=P(ρ,θ)−θPθ(ρ,θ)+Q(θ), (2.2)

Here denotes a pressure potential defined by

 P(ρ,θ)=∫ρ1p(r,θ)r2dr, (2.3)

and accordingly we call the thermal potential which is independent of . The two potentials allow us to rewrite various other terms that arise in our analysis later on in a common form. The spatial derivative of the pressure, for instance, can be expressed as

 1ρ∂xp =1ρ∂x(ρ2Pρ)=∂x((ρP)ρ)−Pθ∂xθ. (2.4)

Another important quantity is the specific enthalpy , which is given by

 h=e+pρ =P−θPθ+Q+ρPρ=(ρP)ρ−θPθ+Q. (2.5)

The temporal change of the internal energy can then be expanded as

 ∂t(ρe) =h∂tρ+(e−h)∂tρ+ρ∂te=h∂tρ+ρ∂te−pρ∂tρ. (2.6)

The last two terms of this splitting turn out to be related to the specific entropy , which is again determined up to a constant by thermodynamic relations [2, 13], namely

 θsρ=eρ−pρ2andθsθ=eθ. (2.7)

A suitable entropy can then be found by integration and reads

 s(ρ,θ)=∫θ1Qθ(t)tdt−Pθ(ρ,θ). (2.8)

The time derivative of this entropy can again be expressed via the potentials as

 ρθ∂ts =−ρθPθρ∂tρ+ρ(Qθ−θPθθ)∂tθ (2.9) =ρ(Pρ−θPθρ)∂tρ+ρ(Pθ−Pθ−θPθθ+Qθ)∂tθ−ρPρ∂tρ=ρ∂te−pρ∂tρ.

Note that the last term already appeared in formula (2.6) for the time derivative of the internal energy. The spatial derivative of the entropy can finally be expressed as

 θ∂xs =θ(sρ∂xρ+sθ∂xθ)=θsρ∂xρ+θsθ∂xθ (2.10) =−θPθρ∂xρ+(Pθ−(θPθ)θ+Qθ)∂xθ=∂x(Q−θPθ)+Pθ∂xθ.

A combination of these two formulas and the Euler equations (1.1)–(1.3) shows, see Remark 3.2 for details, that the evolution of the entropy is governed by

 ρ∂ts+m∂xs =0. (2.11)

Together with (1.1) this yields the conservation law (1.5) for the entropy.

###### Remark 2.1.

To completely describe the constitutive relations for the fluid under investigation, it is thus sufficient to prescribe the two potentials and . Pressure , internal energy , entropy , and derivatives of these quantities can then be expressed in terms of these potentials. This will substantially simplify our analysis later on.

Following physical intuition and to guarantee the well-posedness of the problem under investigation, we make the following structural assumptions about the potentials.

###### Assumption 2.2.

Let and and assume that , , and for all and with some constant .

These assumptions will tacitly be utilized at several places in our analysis below and they are therefore assumed to hold for the rest of the manuscript. To show that they are reasonable, let us briefly discuss a particular example that has been discussed in literature.

###### Example 1 (Admissible state equations).

Let us define

 P(ρ,θ)=cγγ−1ργ−1+C(ρ)θ+c′andQ(θ)=∫θ1cv(t)dt+c′′.

with , , and such that , , and . The constants do not play a significant role in the following. Via (2.2) and (2.3), we can then express the pressure and internal energy as

 p(ρ,θ)=cγργ+c(ρ)θ% ande(ρ,θ)=cγγ−1ργ+Q(θ)+c′.

This yields an extension of the state equations for an ideal gas that has already been considered in [13]. From the definition of the pressure potential, we can further see that

 Pρ(ρ,θ)=cγργ−2+Cρ(ρ)θ≥0for all ρ,θ>0

and

 (ρPρ)ρ=cγ(γ−1)ργ−2+(ρCρ)ρ≥0for all ρ,θ>0.

 Qθ−θPθθ=Qθ=cv(θ)≥c–v>0for all θ>0.

Hence all assumptions about the potentials made above are valid. Note that this simple example already covers, as special cases, the state equations of polytropic gases [8, 19] and of barotropic flow [13, 24]. Our results are therefore directly applicable in these situations.

## 3. An equivalent formulation

We next present an equivalent formulation for the Euler equations which turns out to be particularly well suited for numerical approximation. Recall the definition of the total energy density . With the formulas of the previous section, we then obtain

 ∂tE =∂t(m22ρ2)+∂t(ρe) (3.1) =m(1ρ∂tm−m2ρ2∂tρ)+h(∂tρ)+(ρ∂te−pρ∂tρ).

Note that the last term could also be expressed in terms of the entropy by , which will in fact be utilized below. By the Euler equations (1.1)–(1.3) and the relations between the constitutive equations and the potentials derived in the previous section, the terms in parenthesis on the right hand side of (3.1) can be expressed as

 ∂tρ =−∂xm, (3.2) 1ρ∂tm−m2ρ2∂tρ =−∂x(m22ρ2+(ρP)ρ)−m2ρ2∂xm+Pθ∂xθ, (3.3) ρ∂te−pρ∂tρ =−m∂x(Q−θPθ)−mPθ∂xθ. (3.4)

These equations provide an equivalent formulation for the problem under investigation.

###### Lemma 3.1 (Equivalence).

Let the state equations be defined as in Section 2. Then any smooth positive solution of the Euler equations (1.1)–(1.3) also solves (3.2)–(3.4), and vice versa. The two systems are equivalent in this sense.

###### Proof.

The first equation (3.2) is the same as (1.1).
Next consider (3.3): From equations (1.1)–(1.2), we can deduce that

 1ρ∂tm−m2ρ2∂tρ =−1ρ∂x(m2ρ+p)+m2ρ2∂xm=(i)+(ii)+(iii).

The first term can be expanded as

 −1ρ∂x(m2ρ) =−mρ∂x(mρ)−mρ2∂xm=−∂x(m22ρ2)−mρ2∂xm.

The second term (ii) can be replaced by the expression for derived in (2.4). Adding up all three terms then directly yields the second equation (3.3).
Now consider the third identity (3.4): From equation (3.1), we immediately obtain

 ρ∂te−pρ∂tρ =∂tE−m(1ρ∂tm−m2ρ2∂tρ)−h∂tρ=(∗).

We can now use (1.3), (3.3), and (3.2) to replace the time derivative terms, which yields

 (∗) =−∂x(mρ(E+p))+m(∂x(m22ρ2+(ρP)ρ)+m2ρ2∂xm−Pθ∂xθ)+h∂xm.

Using that and , we can expand the first term as

 −∂x(mρ(E+p)) =−m∂x(m22ρ2)−∂xm(m22ρ2)−m∂x((ρP)ρ−θPθ+Q)−h∂xm.

A combination with the remaining terms then yields the third identity (3.4). Reversing the individual steps yields the other direction and completes the proof of the lemma. ∎

###### Remark 3.2.

Using the formulas (2.9) and (2.10) for the derivatives of the entropy, the third equation (3.4) could also be expressed as

 ρθ∂ts=−mθ∂xs.

Since for any smooth positive solution, this is equivalent to (2.11) and together with equation (1.1) leads to the conservation law (1.5) for the entropy. The above three equations (3.2)–(3.4) thus actually describe the conservation of mass and the evolution of kinetic energy and of entropy. The unknown fields here are the density , the mass flux , and the temperature . This is the basic framework for our further considerations.

## 4. A variational principle

The equivalent formulation (3.2)–(3.4) allows us to establish the following variational characterization of solutions to the Euler equations which will be the starting point and main ingredient for the construction and analysis of numerical approximations later on.

###### Lemma 4.1 (Variational characterization).

Let the state equations be given as in Section 2. Then any smooth positive solution of (1.1)–(1.4) also solves

 (∂tρ,q) +(∂xm,q)=0, (4.1) (1ρ∂tρ−m2ρ2∂tρ,v) −(m22ρ2+(ρP)ρ,∂xv)+(m2ρ2∂xm−Pθ∂xθ,v)=0, (4.2) (ρ∂te−pρ∂tρ,wθ) −(Q−θPθ,∂x(mwθ))+(mPθ∂xθ,wθ)=0, (4.3)

for all test functions , , , and all . Vice versa, any smooth positive solution of (4.1)–(4.3) and (1.4) also solves (1.1)–(1.3).

###### Notation 4.2.

Here denotes the standard scalar product of and is the usual Sobolev space. In addition, we denote by the space of smooth flux functions that additionally vanish at the boundary. Note that on networks [10] or in multiple dimensions, . We therefore use this specific notation already here.

###### Proof of Lemma 4.1.

The assertion follows directly from the equivalence of the Euler equations with (3.2)–(3.4), by multiplying these latter set of equations with appropriate test functions, integration over the domain, and integration-by-parts for some of the terms. Note that all boundary terms vanish due to the homogeneous boundary conditions for and . This shows that any solution of the Euler equations satisfies the above variational principle. The other direction follows by reversing the individual steps. ∎

From the previous two lemmas and the physical principles underlying the Euler equations, one can now immediately deduce global conservation laws for mass, energy, and entropy for all smooth positive solutions of the above variational principle.

###### Lemma 4.3 (Global conservation).

Let be a smooth positive solution of (4.1)–(4.3) and (1.4). Then

 ddt∫ωρdx=0,ddt∫ωEdx=0,andddt∫ωρsdx=0.
###### Proof.

We provide a detailed proof that is only based on the particular form of the variational principle and which carries over directly to the discrete setting considered later on. Conservation of mass follows by testing (4.1) with which yields

 ∫ω∂tρdx=(∂tρ,1)=−(∂xm,1)=−∫ω∂xmdx=0,

where we used the boundary conditions (1.4) in the last step. The second identity is obtained as follows: Using formula (3.1) for the derivative of the total energy, we get

 ddt∫ωEdx =(∂tE,1) Missing or unrecognized delimiter for \Big

Testing the variational principle with , , and leads to

 ddt∫ωEdx =(m22ρ2+(ρP)ρ,∂xm)−(m2ρ2∂xm,m)+(Pθ∂xθ,m) −(∂xm,(ρP)ρ−θPθ+Q)+(Q−θPθ,∂xm)−(Pθ∂xθ,m)=0.

Note that all terms vanish due to the particular structure of the variational principle which is thus responsible for the conservation of energy. In order to verify the conservation of entropy, we proceed as follows: We start by observing that

 ddt∫ωρsdx =(s∂tρ+ρ∂ts,1)=(∂tρ,s)+(θρ∂ts,1θ)=(i)+(ii).

Testing the first equation in the variational principle with leads to

 (i) =−(∂xm,s)=−(∂xm,∫θ1Qθ(t)tdt−Pθ)=(m,1θ∂xQ)+(∂xm,Pθ).

In the last step, we used integration-by-parts and the boundary conditions (1.4) for . The expression (2.9) for the derivative of the entropy and the third equation in the variational principle tested with then allows us to express the second term as

 (ii) =(ρ∂te−pρ∂tρ,1θ)=(Q,∂xmθ)−(θPθ,1θ∂xm−m1θ2∂xθ)−(mPθ∂xθ,1θ) =−(1θ∂xQ,m)−(Pθ,∂xm).

A summation of the two terms now directly yields the conservation of entropy. ∎

## Part II: Discretization

In the following two sections, we discuss the discretization of the variational principle for the Euler equations. We start by a Galerkin approximation in space and then study the subsequent discretization in time by an implicit time stepping scheme.

## 5. Galerkin approximation

For the numerical approximation in space, we consider a conforming Galerkin approximation of the variational principle (4.1)–(4.3) stated above. For this purpose, we choose finite dimensional subspaces , , and . The semi-discrete solution is then defined by the following discrete variational problem.

###### Problem 5.1 (Galerkin semi-discretization).

Let , , be given. Find such that , , , and such that

 (∂tρh,qh) +(∂xmh,qh)=0, (1ρh∂tmh−mh2ρ2h∂tρh,vh) −(m2h2ρ2h+(ρhPh)ρ,∂xvh)+(mh2ρ2h∂xmh−Pθ,h∂xθh,vh)=0, (ρh∂teh−phρh∂tρh,whθh) −(Qh−θhPθ,h,∂x(mhwhθh))+(mhPθ,h∂xθh,whθh)=0,

holds for all test functions , , , and for all .

Here , , , and denote the respective functions evaluated at the discrete solutions. Also note that the solution components , , and depend on , while the test functions , , and are independent of time. We next establish the local well-posedness of this problem.

###### Lemma 5.2.

Let , , and let be given as in Section 2. Then Problem 5.1 admits a unique local solution on for some .

###### Proof.

By definition of the discrete internal energy, we have

 eh=e(ρh,θh)=P(ρh,θh)−θhPθ(ρh,θh)+Q(θh).

Differentiation with respect to time thus yields

 ∂teh =(Pρ,h−θhPθρ,h)∂tρ+(Pθ,h−Pθ,h−θPθθ,h+Qθ,h)∂tθh =(Pρ,h−θhPθρ,h)∂tρ+(Qθ,h−θPθθ,h)∂tθh.

After choosing a basis for , , and , the semi-discrete problem thus leads to a system of ordinary differential equations of the form

 M(Y)Y′+F(Y)=0, (5.1)

where denotes the coordinate vector for the functions , , and . Due to our assumptions on and , the function and the matrix are continuously differentiable as long as and are strictly positive. Moreover, has block triangular form and therefore is regular, if the diagonal blocks are regular. The diagonal blocks, on the other hand, correspond to the mass matrices for the function spaces , , and with weight functions , , and . By our assumptions on the potentials, these functions are strictly positive as long as and are strictly positive. Therefore, the matrix is regular if and are strictly positive. Local existence of a unique solution then follows from the Picard-Lindelöf theorem. ∎

###### Remark 5.3.

By inspection of the proof, one can see that the solution can be extended uniquely in time as long as it remains bounded and and stay strictly positive.

Due to the conforming Galerkin approximation and the particular form of the variational principle, conservation of mass, energy, and entropy also hold on the discrete level.

###### Lemma 5.4.

Let be a solution of Problem 5.1 with and . Then

 ddt∫ωρhdx=0,ddt∫ωEhdx=0,andddt∫ωρhshdx=0.

Here and denote the discrete energy and entropy.

###### Proof.

The assertions can be proved with literally the same arguments as already employed in the proof of Lemma 4.3 on the continuous level. ∎

###### Remark 5.5.

The conservation of energy also allows to obtain appropriate bounds for the norms of , , and . As a consequence, one can expect global existence of the solution, as long as and remain uniformly bounded away from zero. This may be used as a first step towards a complete convergence analysis of the proposed method.

## 6. Time discretization

As a final step of our investigations for the flow on a single pipe, we now discuss an appropriate discretization in time which preserves the underlying conservation laws as good as possible. Given a time step , we define for , and we denote by

 ¯∂τdn:=dn−dn−1τ% for n≥0,

the backward differences which are taken as approximations for the time derivatives. As will become clear below, adaptive time steps could be considered as well.

### 6.1. Fully discrete scheme

For the time discretization of the Galerkin approximation stated in Problem 5.1, we now consider the following implicit time stepping scheme.

###### Problem 6.1 (Fully discrete method).

Set , , and with initial values as in Problem 5.1. For find , such that

 (¯∂τρnh,qh) +(∂xmnh,qh)=0, (1ρn−1h¯∂τmnh−mnh2(ρnh)2¯∂τρnh,vh) −((mnh)22(ρnh)2+(ρP)nρ,h,∂xvh) +(mnh2(ρnh)2∂′xmnh−Pnθ,h∂xθnh,vh)=0, (ρn−1h¯∂τenh−pnhρnh¯∂τρnh,whθnh) −(Qnh−θnhPnh,∂x(mnhwhθnh))+(mnhPnh∂xθnh,whθnh)=0,

holds for all test functions , , and .

As before and similar expressions denote the corresponding functions evaluated at the discrete solutions and . The local well-posedness of this fully discrete scheme can be obtained with similar arguments as that of the semi-discrete problem.

###### Lemma 6.2.

Let be given with and . Furthermore, let the state equations be defined as in Section 2. Then for sufficiently small, the system in Problem 6.1 has a locally unique solution with and .

###### Proof.

With similar arguments as in the proof of Lemma 5.2, one can see that the problem for the th time step can be formulated in algebraic form as

 M(Yn;Yn−1)¯∂τYn+F(Yn)=0.

For we may set with matrix as in the proof of Lemma 5.2. There it was shown that is regular under the assumptions of the lemma. Moreover, and are continuously differentiable with respect to their arguments. Existence of a locally unique solution then follows by the implicit function theorem. For sufficiently small , we can deduce the positivity of and from their continuous dependence on the time step size. ∎

We comment in more detail on the size of the admissible time step below. Before, let us summarize the basic stability and conservation properties of the fully discrete scheme.

###### Lemma 6.3 (Conservation of mass, dissipation of energy, and non-decrease of entropy).

Let denote a positive solution of Problem 6.1. Then

 ∫ωρnhdx=∫ωρkhdx,∫ωEnhdx≤∫ωEkhdx,and∫ωρnhsnhdx≥∫ωρkhskhdx

for all . As before and denote the functions evaluated at , , and .

We call the solution of the discrete problem positive, if and for all .

###### Proof.

It suffices to consider the case ; the results for follow by induction. Testing the first equation with yields the exact conservation of mass

 ∫ωρnhdx =(ρnh,1)=(ρn−1h,1)+τ(¯∂τρnh,1)=(ρn−1h,1)−τ(∂xmnh,1)=∫ωρn−1hdx.

In the last step we used here that on the boundary.
Next consider the energy balance: Using the definition of the total energy , we obtain

 1τ(Enh−En−1h) =¯∂τEnh=¯∂τ((mnh)22ρnh)+(enh+pnhρnh)¯∂τρnh+(ρn−1h¯∂τenh−pnhρnh¯∂τρnh).

The first term in this expression can be further expanded as

 1τ((mnh)22ρnh−(mnh)22ρn−1h)+1τ((mnh)22ρn−1h−(mn−1h)22ρn−1h)≤−(mnh)22(ρnh)2¯∂τρnh+mnhρn−1h¯∂τmnh.

In the last step, we used that for any convex differentiable function there holds

 f(yn)−f(yn−1)≤f′(yn)(yn−yn−1).

This was applied here to the convex functions and , respectively. Substituting in the second term, we obtain after integration over that

 ∫ωEnhdx−∫ωEn−1hdx=(Enh−En−1h,1) ≤τ(1ρn−1h¯∂τmnh−(mnh)22(ρnh)2¯∂τρnh,mnh)+τ(¯∂τρnh,hnh)+τ(ρn−1h¯∂τenh−pnhρnh¯∂τρnh,1).

The time differences can be further evaluated by testing the fully discrete problem with test functions , , and , similar as in the proof of Lemma 4.3. One can then see again that the right hand side above sums up to zero.

The change of the discrete entropy can finally be expressed as

 ∫ωρnhsnhdx−∫ωρn−1hsn−1hdx =(ρnhsnh,1)−(ρn−1hsn−1h,1) =(ρnh−ρn−1h,snh)+(ρn−1h,sn