NUMERICAL INVESTIGATIONS OF THE COMPRESSIBLE NAVIER-STOKES SYSTEM

Abstract. In this paper we write, analyze and experimentally compare three different numerical schemes dedicated to the one dimensional barotropic Navier-Stokes equations: • a staggered scheme based on the Rusanov one for the inviscid (Euler) system, • a staggered pseudo-Lagrangian scheme in which the mesh “follows” the fluid, • the Eulerian projection (on a fixed mesh) of the preceding scheme. All these schemes only involve the resolution of linear systems (all the nonlinear terms are solved in an explicit way). We propose numerical illustrations of their behaviors on particular solutions in which the density has discontinuities (hereafter called Hoff solutions). We show that the three schemes seem to converge to the same solutions, and we compare the evolution of the amplitude of the discontinuity of the numerical solution (with the pseudo-Lagrangian scheme) with the one predicted by Hoff and observe a good agreement.


Introduction
We are interested in the simulation of the following Compressible Navier-Stokes system ∂ t ρ + ∂ x (ρu) = 0, In (1), the unknowns (t, x) → ρ(t, x) and (t, x) → u(t, x) stand respectively for the density and the velocity of a fluid. The quantity µ > 0 is the viscosity of the fluid. We restrict to the isentropic case, where the pressure is 1. Numerical methods

A splitting scheme for compressible Navier-Stokes equations
In this section, we propose a very simple scheme that consists in using a classical hyperbolic type scheme for the Euler system and adding, in a splitting procedure, the effect of the viscosity. To approximate the solution of the Euler system, we use the finite volume framework. The computational domain is the slab (0, 1), endowed with periodic boundary conditions for the sake of simplicity. It is discretized into I cells with equal length C i = (x i−1/2 , x i+1/2 ), and x i being the center of the cell. This amounts to define, for I ∈ N * , ∆x = 1/I, x i+1/2 = i∆x, x i = (i + 1/2)∆x for i = 0, . . . , I. Note that we thus have x I+1/2 = 1, and we define also x 1/2 = 0. For the time discretization, we discuss the schemes with the generic notation ∆t > 0, but it has to be kept in mind that actually ∆t will have to satisfy a local in time Courant-Friedrichs-Lewy type constraint. Therefore, in fact ∆t may depend on the time step (and the local discrete time should be appropriately denoted t n = n k=1 ∆t k instead of n∆t). We denote by X n i = (ρ n i , q n i ) the approximate solution at time t n on the cell C i , where q n i = ρ n i u n i is the discrete momentum: X n i is intended to be an approximation of 1 ∆x´C i X(t n , x) dx. The unknowns are updated as follows Here and below, we consider only first-order three points schemes where the flux F n i+1/2 at the interface x i+1/2 is defined as a function of the unknowns in the neighboring cells: There are many relevant choices for the numerical fluxes. Here, we shall use the Rusanov flux, given by where are the eigenvalues of the Jacobian ∇ ρ,q F (ρ, q). We refer the reader for instance to the textbook [21] for further details about this classical scheme for conservation laws. What is crucial is the following consistency property F is continuous and satisfies F(X, X) = F (X).
In the convergence analysis made below, the details of the numerical fluxes are not important, except property (3). The "splitting" scheme for (1) reads where we have incorporated the diffusion flux at x i+1/2 , here defined by: Periodic boundary conditions are taken into account by setting with ρ n+1 Note that in the case where the viscosity depends on the density, (5)- (6) involves the density at time t n+1 : since it has been already updated by (4) it can be used without any extra cost to update the velocity.

An explicit scheme on staggered grids
We describe now a scheme on staggered grids: the velocity variables are stored on the cells C i+1/2 = [x i , x i+1 ], while discrete densities and pressures are stored on C i = [x i−1/2 , x i+1/2 ]. The discrete mass conservation equation reads where the mass flux F n i+1/2 is defined from the velocity u n i+1/2 known at the interface x i+1/2 . Several definitions of the mass fluxes have been introduced, for instance by using the UpWind flux based on the material velocity [14], an idea reminiscent of the AUSM scheme [23]. Here we use instead the mass fluxes proposed in [2] which are constructed by using the characteristic speeds of the hyperbolic system. Denoting c = p (ρ) the sound speed, we set We note that F ± ≷ 0; the definition relies on the upwinding principles, according to the sign of the characteristic speeds u ± c. For the momentum equation, we set In (10), a natural choice for the discretization of the pressure term would be Π n+1/2 i = p(ρ n i ). However, it turns out that the semi-implicit formula Π , is a better choice, motivated by the consistency analysis of the scheme [2]. The convection fluxes are defined by which, again, rely on upwinding principles. Finally, for the diffusion term, the definition is quite similar to (6) Relations (9) and (10) hold for i = 1, ..., I, imposing furthermore ρ n I+1 = ρ n In the following, this scheme will be referred to as Eulerian staggered scheme 1.

A staggered pseudo-Lagrangian scheme
In this section we propose a pseudo-Lagrangian scheme. We use the term "pseudo-Lagangrian" because although it approximates the solution in the Euler variables, it strongly uses the conservation of the quantities in material volumes. In this scheme, the density unknowns at time t n , ρ n j , are inside the cells (x n j−1/2 , x n j+1/2 ), which are time-dependent, and the velocity unknowns u n j+1/2 are associated with the boundaries of the cells, the x n j+1/2 's. The latter are advected at velocity u n j+1/2 . The discretization of the mass equation corresponds to the discretization of the relation d dtˆy with x(t) and y(t) moving at the fluid velocity. It expresses that the mass in every cell is constant in time. Namely, we obtain The second equation in (11) accounts for the periodicity of the problem. We turn to the momentum equation

It is approximated by
where the last four equations stand for the periodicity. In these formulae, we use the quantities ρ n+1 j+1/2 , ∆x n+1 j+1/2 and p n+1 j that were not defined previously: we choose for them the following rather natural definition: Complemented with initial conditions ρ 0 j , u 0 j+1/2 , x 0 j+1/2 , the scheme is defined. Note that we have chosen an implicit approximation for the diffusion, in order to avoid any parabolic stability condition (and indeed we observe the stability of the algorithm).

A last staggered scheme
In this section we propose another staggered scheme, that can be viewed as a Eulerian version of the preceding pseudo-Lagrangian one. Here the cells are fixed and, starting from the pseudo-Lagrangian scheme, we have to take the transport into account. With the same notations, but now the length ∆x of the cells is constant, the mass conservation is replaced with where ρ n j+1/2,U is upwinded, that is to say that it is ρ n j if u n j+1/2 ≥ 0 and ρ n j+1/2 if u n j+1/2 < 0. The momentum equation is approximated by where the last four equations stand for the periodicity. In these formulae, we use the quantities ρ n+1 j+1/2 and p n+1 j that were not defined previously: we choose for them the following rather natural definition: Complemented with initial conditions ρ 0 j , u 0 j+1/2 , x 0 j+1/2 , the scheme is defined. Note that we have chosen an implicit approximation for the diffusion, in order to avoid any parabolic stability condition (and indeed we observe the stability of the algorithm). The convection term in the momentum equation is centered, but for fixed (not too small) viscosity coefficient, this does not make instabilities appear.
In the following, this scheme will be referred to as Eulerian staggered scheme 2.

Convergence analysis à la Lax-Wendroff
In this section, we investigate the convergence of the Rusanov scheme. For the sake of simplicity, we restrict the analysis to the case of uniform meshes: all cells have the same length. We consider a sequence of meshes, parametrized by k ∈ N, such that ∆x k and ∆t k both tend to 0, with ∆x k ∆t k = λ for some λ > 0 to ensure some stability property, for instance the positivity of the density.
The discrete initial data are given by The discrete approximation of ρ and u are denoted by We bear in mind that ρ n j and u n j also depend on k but this dependency is omitted to simplify the notation. The following statement is in the same spirit as the standard Lax-Wendroff theorem for conservation laws, see [21, Section 12.10].
Remark 2.2. The statement remains true when µ(ρ) = Cρ at the price of assuming u ∈ L r (0, T ; W 1,r (0, 1)) and u It equally applies for a continuous function µ, with the additional assumptions that the convergence of (ρ k , u k ) to (ρ, u) holds almost everywhere.
Proof. Consider a test function ζ = (ϕ, χ) ∈ C ∞ 0 ([0, T ) × (0, 1); R 2 ). We assume that ∆t k and ∆x k are small enough so that supp(ζ) ⊂ [0, T − ∆t k ) × (2∆x k , 1 − 2∆x k ). For n ∈ {0, ..., N − 1} and j ∈ {1, ..., I}, let us denote ζ n j = ζ(t n , x j ), and, for all t ∈ (0, T ), x ∈ (0, 1), We multiply the discrete mass conservation equation by ∆x k ∆t k ϕ n j , and we sum over j and n. We obtain The first term on the left hand side recasts as while the second term becomes Hence we get Using (3), the estimates and convergence assumed for (ρ k , u k ), we obtain as k → ∞ We turn to the momentum equation (5). Proceeding similarly, we obtain The first two sums can be written as which can be treated as for the mass conservation equation, by using that the product ρ k u k still converges in L r ((0, T ) × (0, 1)) for 1 ≤ r < +∞. Let us study the viscous term The case where the viscosity µ is constant can be handled by using two summations by parts, which allow us to write (19) as This can be cast aŝ and using the uniform convergence of χ k and its derivatives, together with the convergence in L r of u k , we are able to pass to the limit: as k → ∞, we obtain which completes the proof in the case of constant viscosity.
When µ depends on the density, we cannot proceed this way, and we need further assumptions. We start by writing (19) as To pass to the limit in the previous term we use a weak-strong convergence argument observing that This combines to the assumption that ∇u k converges weakly in L r ((0, T ) × (0, 1)) so that (20) tends tô T 0ˆ1 0 µ(ρ(t, x))∂ x u(t, x) ∂ x χ(t, x) dx dt, which concludes the proof.

Numerical results
Here, we present some numerical results that illustrate the behavior of the proposed schemes.

Comparisons
In this subsection, we compare the results obtained with the 4 different schemes that are dealt with in the paper. The results here are obtained at time t = 0.1 with µ = 0.1 for the initial condition u 0 (x) = 0 and ρ 0 (x) = 0.125 + 1.875χ [1/4,3/4]   We see that the four schemes seem to converge to the same solution, and that the pseudo-Lagrange scheme is less diffusive (the scheme based on the Rusanov flux is the most diffusive). In particular, the pseudo-Lagrange scheme is able to maintain some discontinuities in the density and in the velocity derivative.
For the same test-case, we now more precisely compare these four schemes, in the velocity variable. In Table 1 we compare, for different numbers of cells, the non-diffusive pseudo-Lagrangian scheme from Section 1.3 and the Eulerian Rusanov scheme from Section 1.1. The difference between these schemes seems to be of order 1/2, both in the L 1 and in the L ∞ norms. Note that the schemes are subject to some hyperbolic-typed Courant-Friedrichs-Lewy conditions, so that the time step is of the same order as the space step. Note also that to compare the solutions, as they are not obtained on the same meshes (the pseudo-Lagrange scheme is running on a moving mesh), we interpolate the difference of the two solutions on a mesh that is the intersection of the two and consider both solutions as constant-by-cell. In Table 2 we compare the pseudo-Lagrange and the Eulerian staggered scheme 1 presented in Section 1.2. In Table 3 we compare the pseudo-Lagrange and the Eulerian staggered scheme 2 form Section 1.4. We observe the same order of convergence.  Table 3. Difference between the Pseudo-Lagrange scheme and the Euler staggered scheme 2.

Behavior of the pseudo-Lagrange scheme on Hoff-type discontinuous solutions
In this section, we would like to illustrate a very interesting property of discontinuous solutions à la Hoff. These solutions, among which one finds the one illustrated in the previous section, have been studied by Hoff (see [15]). They consist in solutions with discontinuities in density, and thus in pressure and in the derivative of the velocity (because the effective flux p − µ∂ x u is continuous in space for almost every time). In [15] it is shown that the amplitude of the jump in log(ρ) decreases exponentially in time, with a rate that is at least some constant over the viscosity coefficient µ. We do not precise the result here because the constant depends on the maximum value of the density in the flow, and because actually the context of [15] is a bit different: it considers the Navier-Stokes system in Lagrange variables with a space variable ranging over R, with a single discontinuity, while here we consider the Navier-Stokes system on the torus. We propose to analyse numerically the behavior of the jump of log(ρ). This is made possible thanks to the use of a pseudo-Lagrange scheme that does not smooth the discontinuity in ρ artificially. With the same initial condition as in the previous section, Figure 2 presents the amplitude of the jump in log(ρ) with respect to time, with 500 cells in space. We also present its best approximation as exp(a + bt) by computing a and b by the least square method. Note that it is not relevant to do this on a too large time interval (with a given cell size), because the amplitude of the jump decays exponentially fast and it rapidly becomes of the same order as the cell size so that it cannot be distinguished from the natural steps of the discrete solution in regular regions.  Now we want to evaluate the rate at which the jump decreases as a function of the viscosity coefficient.  Table 4. Rate for the (decreasing) amplitude of the jump of the logarithm of ρ.
We note in Table 4 that the jump [log(ρ)] (with usual notations for the jump of a quantity) seems to behave as exp(−Ct/µ), which is exactly the majoration proved in [15]. Note that this is a very interesting phenomenon: as µ tends to 0, this kind of discontinuities do not persist at all (they disappear at time t = 0 + ), what is a paradox, as for compressible inviscid gases discontinuities are expected (and they are called shocks). This result by Hoff thus proves that the shocks in inviscid gases are limits of smooth parts in the solutions of viscous gases, in which the amplitude of discontinuities tends to 0.