Convergent finite differences for 1D viscous isentropic flow in Eulerian coordinates

We construct a new finite difference method for the flow of ideal viscous isentropic gas in one spatial dimension. For the continuity equation, the method is a standard upwind discretization. For the momentum equation, the method is an uncommon upwind discretization, where the moment and the velocity are solved on dual grids. Our main result is convergence of the method as discretization parameters go to zero. Convergence is proved by adapting the mathematical existence theory of Lions and Feireisl to the numerical setting.


Introduction
In this paper, we will develop a convergent finite difference method for the flow of an ideal viscous isentropic gas in one spatial dimension. We will assume that the flow may be modeled by the Navier-Stokes system (cf. [16]): The unknowns in this system are the fluid density = (t, x) and the fluid velocity u = u(t, x). For an isentropic flow, the ideal pressure law takes the form p = a γ , a > 0, where the value of γ is determined by the specific gas in question. In this paper, we will be forced require that to prove convergence of the method. Note that this significantly limits the physical applicability of the convergence result. While monoatomic gases (γ ∼ 5 3 ), such as helium, are included, diatomic gases (γ ∼ 7 5 ), such as air, are not. The condition is necessary for convergence, but not stability, of the method. Existence of global classical solutions for the system (1.1)-(1.2) is wellknown and was first established by Kanel' [9] for smooth data (see also [14]). For non-smooth initial data, corresponding results have been obtained by David Hoff [1,7]. These works are all based on the same approach: find pointwise upper and lower bounds on the density, then use these bounds to estimate the second derivative of the velocity. With smooth initial density, the density remain H 1 , while with discontinuous initial density, the density is at most BV as initial discontinuities persists for all time (see [7]). The pointwise bounds on the density are obtained by tracking certain quantities along streamlines rendering the entire existence theory essentially Lagrangian.
In the literature, one can find a huge variety of numerical methods appropriate for (1.1)-(1.2). However, very few of these methods have been proven to converge. This lack of rigorous results are most likely a consequence of the Lagrangian nature of the existence theory. In fact, prior to this paper, all convergent methods have been discretized in Lagrangian coordinates and are due to David Hoff and collaborators [18,19,20]. That being said, there are also several existence results that utilize discrete Lagrangian approximation schemes (e.g. [1,7,8]) to construct solutions. For practical applications, results in Eulerian coordinates is often desirable and a mapping from discrete Lagrange to discrete Euler is costly. For this reason, most practitioners would employ an Eulerian discretization.
In this paper, we will discretize the equation in Eulerian coordinates. As a consequence, obtaining pointwise bounds on the density becomes highly involved and will not be pursued in this paper. Instead, we will develop a convergence theory in the spirit of the continuous existence theory [4,15,17] as sparked by P. L. Lions. In particular, we will not obtain any form of continuity of the discrete density and instead prove strong convergence of the density using renormalization and what is known as the effective viscous flux (cf. [17]).
In more than one spatial dimension, the literature is almost void of convergent numerical methods. For the full system, the only result is the paper [13], by the author, in which a convergent finite element discretization of (1.1)-(1.2) is developed. Over the last years, there have also been developed a series of convergent methods [2,3,5,10,11,12] for the Stokes version of (1.1)-(1.2). The method we shall develop and analyze in this paper can be considered as an archtype for all the methods mentioned above. That is, in one spatial dimension, all of these methods are in a sense equivalent and of the form we shall consider here. In addition, the method we present here is the first finite difference method for which convergence is proved.
The paper is organized as follows: In the next section, we will define the numerical method and state the main convergence result. Then, we will derive stability of the method and some other basic properties. In Section 4, we will derive an equation for the effective viscous flux and also provide a higher integrability estimate for the density. Section 5 concerns passing to the limit in the method. In particular, we prove that the limit almost satisfies (1.1)-(1.2). The remaining ingredient is to prove strong convergence of the density which is proved in the final section.

The numerical method and main result
Our method will be posed on a uniform grid in both space and time. In time, we shall approximate at discrete points t k = k∆t, where ∆t is assumed to be of the order ∆x. In space, we will divide the domain [0, L] into N intervals of length ∆x = L/N . It will be convenient to write this as where we have introduced the slightly confusing notation We shall also need the dual grid given by the midpoint nodes The numerical method will approximate the density on the dual grid and the velocity on the standard grid: Staggered grid of this kind is necessary for incompressible flows, and as a consequence also widely used to develop all speed compressible flow methods.

Discrete operators.
To discretize the convective terms, we shall utilize an upwind method. To facilitate this, we introduce the notation u + = max{u, 0}, u − = min{u, 0}.
We shall also need the average velocity over an interval For the continuity equation (1.1), we shall use the following upwind flux This upwind flux will also be used for the momentum equation, where we upwind the (averaged) momentum The remaining derivatives will be discretized using the operators which defines the following (standard) Laplace operator For time discretization, we will use implicit time-stepping The numerical method is defined as follows: Definition 2.1. Given ∆t, ∆x < 1, and initial data (1.3), define the numerical initial data Determine sequentially the numbers ( k i , u k i ), i = 0, . . . , N − 1, k = 0, . . . , M, solving the nonlinear system It is not completely trivial that our numerical method is well-defined. Since the system of equations (2.1)-(2.2) are both nonlinear and implicit, existence of a solutions needs to be established. Lemma 2.2. Let 0 < ∆t, ∆x < 1 be fixed. The nonlinear implicit system of equations (2.1)-(2.2) admits at least one solution.
Proof. The existence can be proved using a toplogical degree argument. Since the proof is very similar to the corresponding result in [6,11,13], we do not give the details here and instead refer the reader to any of the mentioned papers.
From [6] and [11], we also have that the method preserves strict positivity of the density: To analyze the finite difference method, it will be of great convenience to extend the numerical solution to all of [0, T )×(0, L). Since the density appears nonlinearly in the pressure, we will use piecewise constants to extend it: For the velocity, we shall utilize piecewise continuous linears in space and piecewise constants in time.
Note that an extension of the average velocity u k i is now given by the L 2 projection onto piecewise constant: Then, as h → 0, u h u in L 2 (0, T ; W 1,2 0 (0, L)), h → a.e on (0, T ) × (0, L), where ( , u) is a weak solution of (1.1)-(1.2): Theorem 2.4 will follow as a consequence of the various results stated and proved in the upcoming sections. The proof will be completed in Section 6.

Stability and energy estimates
In this section, we we will derive a numerical analog of the continuous energy estimate. This estimate yields in particular stability of the method, but it will also provides us with necessary L p bounds uniform in the discretization parameters.
3.1. Rernormalized continuity scheme. We shall need the following renormalized continuity scheme at several occasions. The term renormalized is motivated by the corresponding continuous equation and its role in the existence theory (cf. [15]).
satisfies the continuity scheme (2.1), then the following identity holds where P is given by respectively. Proof. We begin by multiplying (2.1) with B ( ) to obtain . By applying Taylor expansion, we find that where * is some number in [ k−1 i , k i ]. By adding and subtracting and applying another Taylor expansion, Here, † i and ‡ i are some numbers in [ i , i+1 ] and [ i , i−1 ], respectively. Next, we calculate We conclude the proof by combining (3.4) in (3.3) and (3.2).

3.2.
The convection operator. To prove stability of the method, we shall multiply the momentum scheme (2.2) by u k i+1/2 and sum over all i. For this purpose, we shall need the following identity for the convective discretization.
Lemma 3.2. The following identity holds where the numerical diffusion term N 1 is given by Proof. By applying summation by parts, we see that (3.5) Next, we apply the definition of Up(·) and add and subtract to deduce We then conclude the proof by setting this identity in (3.5).
3.3. Energy estimate. We are now ready to prove stability of the method.
where the numerical diffusion terms are Proof. We begin by multiplying (2.2) by u k i+1/2 ∆x and sum over all i to obtain where we have also applied summation by parts together with the boundary condition u k −1/2 = u k N −1/2 = 0. Next, we apply Lemma 3.2 and Lemma 3. (3.7) To proceed, we observe the following identities and We conclude by multiplying with ∆t, summing over k = 1, . . . , M , and recalling the definition of P[·].
The previous stability estimate also provides some uniform integrability estimates on various quantities. To state these, let us introduce the notation to denote the case when f h is bounded in L p uniformly with respect to both ∆t and ∆x.

As a consequence,
The remaining bounds follows directly from these using Sobolev embedding and interpolation estimates.

Control in time.
We end this section by establishing some weak control on the time-derivatives. For this purpose, let for all k = 1, . . . , M .
Proof. 1. We first prove (3.10). Let φ ∈ C ∞ (0, L) be arbitrary and define Now, multiply the continuity scheme (2.1) with ∆xφ i , sum over all i, and apply summation by parts, to deduce .
The last inequality is an application of the Hölder inequality. Now, since (3.12) holds for all φ, we can conclude that it continues to hold for all φ ∈ W 1, γ γ−1 (0,L) By multiplying with ∆t and summing over all k, we obtain By multiplying the momentum scheme (2.2) with ∆xv i+1/2 , summing over i, and applying summation by parts, we obtain . Next, we multiply with ∆t and sum over all k to conclude where the last inequality is Corollary 3.4. Finally, we calculate form which (3.11) follows.

The effective viscous flux
The purpose of this section is to derive an equation for the quantity This quantity is often termed the effective viscous flux and stands at the center of both the available existence results (in more than 1D) [15,4] and the available numerical results (cf. [13]). Specifically, both higher (than L 1 ) integrability on the pressure and strong convergence of the density is proved using the upcoming equation (see Proposition 4.2 for details).
To derive the desired equation, we shall need the following discrete Neumann Laplace operator and we observe that −∆ i is nothing but the standard 3-point Laplacian on the density grid. Moreover, the Neumann condition is realized by adding a shadow cell on both sides of the domain (0, L). Due to the Neumann condition, ∆ i is only well-defined for sources f h of zero mean.
In the upcoming analysis, we shall not need ∆ i directly, but its discrete derivative: where we observe that the Neumann condition renders We shall also need the inverse of the ∆ i+1/2 operator occurring in the momentum scheme (2.2): We will mostly be interested in the discrete derivative The following result follows from standard summation by parts.
Proof. By direct calculation, using the definition of ∆ i and ∆ i+1/2 , and the respective boundary conditions, we deduce which concludes the proof.
The following proposition gives the effective viscous flux equation we shall need in the convergence analysis. The two error terms appearing in (4.1) will be bounded below.
where the numerical error terms E m 1 and E m 2 are given by (4.6) and (4.8), respectively.

CONVERGENT FINITE DIFFERENCES FOR 1D VISCOUS ISENTROPIC FLOW 13
and summing over all i and k = 0, . . . , m, we obtain the starting point Here, we have introduced the quantities We will now rewrite the sums S 1 and S 2 using the definition of v k i+1/2 and the continuity scheme. We will treat the least involved term, namely S 2 , first. For this purpose, we apply summation by parts and the definition of v k i+1/2 to calculate For the S 1 term, we first apply summation by parts in time followed by an application of Lemma 4.1 to deduce Next, we multiply the continuity scheme (2.1) by ∆t∆xq k i , where and sum over all i and k = 0, . . . , m to obtain ∆t∆x m k=0 i where the last identity follows from summation by parts, the definition of ∆ −1 i+1/2 , and where we have introduced the error term To conclude the proof, we shall need to further rewrite (4.5). By adding and subtracting, we obtain the identity Consequently, by combining this identity (4.5) and (4.4), we conclude where Finally, we apply (4.3) and (4.7) to (4.2) to discover This and a final summation by parts, concludes the proof.

4.1.
Bound on the error terms in (4.1). In order for the previous proposition to be useful, we will need to prove suitable bounds on the error terms E h 1 and E h 2 . It is at this stage we will need to impose some requirements on the value of γ. No other results in the paper imposes any unphysical restrictions on γ.
We begin by deriving a bound on E h 2 . and let E h 2 be given by (4.8). There is a constant C > 0, independent of ∆t and ∆x, such that for any m = 1, . . . , M , Proof. By two applications of the Cauchy-Schwartz inequality and a standard inverse estimate, we deduce where the norms are bounded due to Corollary 3.4.
We are now ready to bound the other error term in (4.1). The following lemma is the sole reason for the requirement γ > 3 2 .
Lemma 4.4. Assume that the adiabatic coefficient satisfies and let E h 1 be given by (4.6). There exists a constant C > 0, independent of discretization parameters, such that for any m = 0, . . . , M , Proof. From (4.6), we have that Let us now examine one of the terms in E h 1 . By adding and subtracting, we write (4.10) To bound the S 1 term, we multiply and divide by √ p ( ‡ i ) and apply the Cauchy-Schwartz inequality (4.11) Next, we make several applications of the Hölder inequality together with standard inverse estimates, to deduce where we have used Corollary 3.4 to conclude the last inequality. Combining (4.12)-(4.11) and recalling that ∆t = ∆x, yields Next, we turn to the S 2 term in (4.10). An application of the Cauchy-Schwartz inequality yields  Next, we proceed as in (4.12) to discover where we have used Corollary 3.4 in the last inequality. By combining (4.14) and(4.15), and recalling that ∆t = ∆x, we conclude By setting (4.13) and (4.16) in (4.10) we obtain and hence we have the desired bound for the first term in (4.9). The second term in (4.9) can be bounded by the exact same arguments.

4.2.
Higher integrability on the density. From Corollary 3.4, we only know that p( h ) is uniformly (in ∆t and ∆x) bounded in L ∞ (0, T ; L 1 (0, L)). Hence, it is unclear whether p( h ) actually converges to an integrable function. In the following lemma, we prove that the pressure has more integrability than provided by the energy estimate. There is a constant C > 0, independent of discretization parameters, such that In other words, p( h ) ∈ L ∞ (0, T ; L γ+1 γ (0, L)).
Proof. First, we rewrite the equation (4.1) in the form (4.1) To bound the terms involving the discrete inverse Laplacian, we shall use the elementary bound The proof is completed by fixing sufficiently small.

Weak convergence
In this section, we will pass to the limit in the numerical method and prove that the limit is almost a weak solution to the compressible Navier-Stokes equations. Our starting point is that Corollary 3.4 allow us to assert the existence of functions u ∈ L 2 (0, T ; W 1,2 0 (0, L)), ∈ L ∞ (0, T ; L γ (0, L)),

(5.2)
Note that we cannot make the identification p( ) = p( ) as this would require the density to converge strongly. We will prove that this is indeed true in the next section.
To conclude convergence of the product terms appearing in the method, we shall need the following lemma from [11]: converges weakly to f and g in L p 1 (0, T ; L q 1 (0, L)) and L p 2 (0, T ; L q 2 (0, L)), respectively, where 1 < p 1 , q 1 < ∞ and 1 gf in the sense of distributions on (0, T ) × Ω.

(5.3)
Proof. From Corollary 3.4 and Lemma 3.5, we have that . We can then apply Lemma 5.1, with g h = h and f h = u h , to conclude h u h u in L 2 (0, T ; L γ (0, L)).
Next, we notice that where the last inequality is the Poincaré inequality. Hence, by writing (5.4) and passing to the limit, we conclude the second convergence of (5.3). From Corollary 3.4 and Lemma 3.5, we have that The remaining convergences can be obtained similarly (i.e (5.4)).

5.1.
Convergence of the density scheme. We now prove that the limit ( , u) is a weak solution of the continuity equation. Proof. 1. We first claim that ( h , u h ) satisfies the equation To prove the claim (5.5), let φ ∈ C ∞ 0 ([0, T ) × [0, L]) be arbitrary and define Now, multiply (2.1) with φ k i ∆x, and sum over all i, to discover To proceed, we add and subtract to derive the identities and By applying (5.7)-(5.8) in (5.6), we obtain We then apply Green's theorem to obtain where we have defined (5.10) To proceed, we observe the two identities and ∆t k i where the last equality follows by definition of φ k i By combining (5.9)-(5.12), we obtain (5.5), which was our claim.
2. Next, we prove that the P 1 (φ) converges to zero as h → 0. More precisely, we claim that (5.13) To prove this, we apply the Cauchy-Schwartz inequality to obtain (5.14) Note that the first term after the inequality is bounded by the energy estimate (3.6). Now, the Hölder and Poincaré inequalities provides the bound 3. Let us now send h → 0 in (5.5) and thereby conclude the proof. For this purpose, we shall need the following elementary identity where we have also used that φ(T, ·) = 0. With this identity, (5.5) tell us that From (5.13), we have that P 1 (φ) → 0 as h → 0. Hence, there is no problem with sending h → 0, using (5.2) and (5.3), to conclude that in the sense of distributions. This concludes the proof.

5.2.
Weak limit of the momentum scheme. In this subsection, we pass to the limit in the momentum scheme to conclude that ( , u) is almost a weak solution of the momentum equation. We begin by deriving an integral formulation of the momentum scheme (2.2).
where the numerical error term is given by Proof.
Let v ∈ C ∞ 0 ([0, T ) × 0, L) be arbitrary and introduce the notation Now, multiply (2.2) with v k i+1/2 ∆x, and sum over all i to obtain We will now write each integral term in (5.16) in the form (5.15). Let us begin with the right-hand side. 1. Using summation by parts, we calculate which is of the form we wanted.

2.
For the time derivative term, we perform summation by parts to deduce 3. For the last term, we begin applying summation by parts Now, by adding and subtracting, we develop the identity Similarly, we derive the identity By applying (5.20) and (5.21) in (5.19), we obtain (5.22) 4. By applying (5.17), (5.18), and (5.22), to (5.16), multiplying with ∆t, and summing over all k,we discover which is (5.15) with This concludes the proof.
Next, we prove that the error term in (5.15) converges to zero as h → 0.
Lemma 5.5. Let P 2 (v), be as in the previous lemma. There is a constant C > 0, independent of h, such that Proof. By definition, we have that Let us now bound the S 1 and S 2 term separately.
1. Bound on S 1 : By adding and subtracting, we write (5.24) To bound the K 1 term, we apply the Cauchy-Schwartz and Poincaré inequalities to obtain where we have used that the term involving ∂ k t u h is bounded by the energy estimate (3.6).
To bound the K 2 term, we apply yet another application of the Cauchy-Schwartz and Poincaré inequalities to discover Setting (5.25)-(5.26) in (5.24) yields 2. Bound on S 2 : By adding and subtracting, we calculate (5.28) To bound the K 3 term, we apply the Cauchy-Schwartz inequality, followed by the Hölder inequality, and obtain where we have also utilized a standard inverse estimate in time and the energy estimate (3.6).
5.3. The effective viscous flux limit. We are going to end this section by passing to the limit in the effective viscous flux equation (4.1). To achieve this, we shall need the following result.
• The discrete time derivate of f h satisfies Then, for any t ∈ (0, T ), where k is given by k = t/h and ∆ D denotes the Laplace operator with homogenous Dirichlet conditions.
Proof. Let q h be the extension of From the requirements on f h , it is clear that Then, we can apply Lemma 5.1, with f h = g h = q h to conclude that q 2 h q 2 in the sense of distributions. Thus, q h → q a.e as h → 0, where the convergence may take place along a subsequence. Observe that the limit must satisfy , where ∆ D is the Laplacian with Dirichlet conditions. Next, let q L h be the linear-in-time interpolation of q h : for k = 0, . . . , M − 1. Since we have that we have that q L h → q in C(0, T ; L p (0, L)) and sup t q L h − q h L p (0,L) → 0, for any p < ∞. Finally, we write and send h → 0 to discover (5.34).
Equipped with the previous lemma, we now pass to the limit in the effective viscous flux equation (4.1).
Proof. Let m = t/∆t , then (4.1) provides the identity (5.36) By virtue of the previous lemma, we know that the product terms involv- converges to the product of the limits. Furthermore, from Lemmas 4.4 and 4.3, we know that the E m 1 and E m 2 terms converges to zero as h → 0. Hence, the only term that might be problematic is the term involving Up( uu).
By adding and subtracting, we write where the last term is bounded as where the last inequality is a standard inverse estimate. The two previous equations tell us that Hence, there is no problems with passing to the limit h → 0 in (5.36) and conclude the proof.

Strong convergence (proof of Theorem 2.4)
In the previous section, we proved that the method almost converges to a weak solution of the compressible Navier-Stokes equations. To conclude convergence of the method, and thereby conclude the main theorem, it only remains to prove that p( ) = a.e. This is the topic of this section. To make this identification, we shall prove that h converges strongly a.e. Our strategy will be to adopt the continuous arguments of Lions [15] to the numerical setting. For this purpose, we shall need the following well-known lemma. The reader may consult [4] for a proof. From Lemma 5.3, we know that the limit ( , u) is a weak solution of the continuity equation and moreover that ∈ L 2 (0, T ; L 2 (0, L)) and u ∈ L 2 (0, T ; W 1,2 0 (0, L)). The previous lemma is then applicable. By setting B(z) = z log z and integrating over the domain, we obtain where we have also used convexity of the map z → B(z) to conclude a sign on the error terms. By subtracting (6.1) from (6.2) and passing to the limit h → 0, we deduce for any t ∈ (0, T ). Here, the overbar denotes weak L 1 -limit and the first inequality is a consequence of convexity. The remaining ingredient to prove strong convergence of the density is to prove that the right-hand side in (6.3) is negative: As a consequence, h → a.e as h → 0.
Proof. Let ∆ N denote the Laplace operator with homogenous Neumann conditions on (0, L). From (5.1) and Lemma 3.5, we have that ∈ L ∞ (0, T ; L γ (0, L)) and t ∈ L 2 (0, T ; W −1,γ (0, L)), respectively. As a consequence, satisfies v x , v t ∈ L 2 (0, T ; L γ (0, L)). Hence, we can set v as test-function in (5.32) to obtain where the last equality is integration by parts and the duality of the Neumann and Dirichlet Laplacians ∆ N and ∆ D , respectively. To proceed, we note that φ = ∂ x ∆ −1 D [ u] is a valid test-function for the continuity equation (1.1). Since ( , u) is a weak solution of the continuity equation (Lemma 5.3), we deduce the identity By setting (6.6) in (6.5), we discover the identity where the last equality is (5.35). Finally, we rewrite (6.7) in the form where the last inequality follows from the convexity of z → p(z). This concludes the proof of (6.4).
Proof of Theorem 2.4: From Lemma 5.6, we have that ( , u) solves in the sense of distributions. The previous lemma tell us that p( ) = p( ) a.e in (0, T ) × (0, L).