# TVD analysis of a (pseudo-)staggered scheme for the isentropic Euler equations

Katia Ait-Ameur, Michael Ndjinga

► To cite this version:

Katia Ait-Ameur, Michael Ndjinga. TVD analysis of a (pseudo-)staggered scheme for the isentropic Euler equations. Springer Proceedings in Mathematics & Statistics, 433, pp.249-257, 2023, FVCA 2023: Finite Volumes for Complex Applications X-Volume 2, Hyperbolic and Related Problems, 10.1007/978-3-031-40860-1_26. hal-04128390

HAL Id: hal-04128390

https://hal.science/hal-04128390

Submitted on 14 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
TVD analysis of a (pseudo-)staggered scheme for
the isentropic Euler equations

Katia Ait-Ameur¹ and Michaël Ndjinga²

¹ CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, France,
katia.ait-ameur@polytechnique.edu,

² Université Paris-Saclay, CEA Saclay, DES,ISAS,DM2S,STMF, France,
michael.ndjinga@cea.fr

**Abstract.** In this paper, we build and analyze the stability of a collocated scheme, involving a specific numerical diffusion operator, for the isentropic Euler equations. This scheme is based on the numerical diffusion operator of a family of staggered finite volume schemes introduced in [1]. The properties of this operator allowed to understand the $L^2$-stability of staggered finite volume methods. Staggered schemes are popular in the thermal hydraulics community for their reported robustness and lack of spurious oscillations. The contributions of this paper are twofold: we firstly build a colocated scheme with a staggered-based numerical diffusion operator, hence the proposed terminology of "pseudo-staggered" scheme, and we secondly present a rigorous TVD analysis, in order to explain the non-oscillatory behaviour of staggered discretisations.

**Keywords:** Euler equations, compressible flows, finite volumes, staggered grids, TVD analysis

# 1 Introduction

Staggered finite volume discretisations were introduced in the 1960s for incompressible flows [5], then nearly incompressible flows [2]. Although non conservative, they are nowadays very popular in the thermal hydraulics community of two phase flows [4], and used in many nuclear safety simulation platforms. However, their stability analysis remained based on heuristics [3], and practical simulations. Starting from 2008, entropy-stable staggered discretisations have been proposed for the Euler system [6]. In [1], a framework is proposed for the rigourous analysis of the linearised $L^2$-stability of these staggered schemes [4] in view of designing robust conservative $L^2$-stable for multiphase flows.

Beyond stability, one of the purported advantages of the staggered discretisation is their apparent lack of spurious oscillations. This paper proposes a TVD analysis to explain the good behaviour of staggered discretisation for low Mach flows. We perform a rigourous TVD analysis in the case of the 1D pseudo-staggered scheme introduced in [1] for the isentropic Euler equations as a first step before studying other discretisations and more advanced models. In section
2 we recall the numerical method, then we compute in section 3 the equations giving the evolution of the discrete Riemann invariants. Finally in section 4 we give a partial result obtained in the case of low Mach number flows. Some first numerical results are given in Section 5.

## 2 A pseudo-staggered scheme

We consider the following 1D isentropic Euler equations in conservative form:

$$ \begin{cases} \partial_t \rho + \partial_x q = 0 \\ \partial_t q + \partial_x \frac{q^2}{\rho} + \partial_x p = 0. \end{cases} \qquad (1) $$

where $t$ stands for the time, $\rho$, $q$ and $p$ are the density, momentum and pressure respectively. The pressure $p$ is a function of the density $p = f(\rho)$ with $f$ a strictly increasing function. System (1) takes the conservative form

$$ \partial_t U + \partial_x F(U) = 0, \quad U = \begin{pmatrix} \rho \\ q \end{pmatrix}, \quad F(U) = \begin{pmatrix} q \\ \frac{q^2}{\rho} + p \end{pmatrix}. \qquad (2) $$

We propose a conservative finite volume discretisation of (1)

$$ \partial_t U_i + \frac{F_{i,i+1} - F_{i-1,i}}{\Delta x} = 0, \quad U_i = \begin{pmatrix} \rho_i \\ q_i \end{pmatrix}, \qquad (3) $$

where the numerical flux takes the form

$$ F_{i,i+1}(U_i, U_{i+1}) = \frac{F(U_i) + F(U_{i+1})}{2} + D_{\text{stag}}(U_i, U_{i+1}) \frac{U_i - U_{i+1}}{2}, \qquad (4) $$

where the matrix $D_{\text{stag}}$ is the diffusion operator of the scheme and is the same as the one of the staggered scheme [1, Corollary 1]

$$ D_{\text{stag}}(U_i, U_j) = \begin{pmatrix} |u_{ij}| - u_{ij} & 1 \\ -c^2 - u_{ij}^2 & 2u_{ij} \end{pmatrix}. \qquad (5) $$

where $c$ is the sound speed, and $u_{ij}$ is the Roe average of the velocity between the two states $U_i$ and $U_j$. This matrix has complex eigenvalues $|u_{ij}| \pm ic$ but still satisfies the stability condition introduced in [1] : $\tilde{P}(D_{\text{stag}} + {}^t D_{\text{stag}})\tilde{P}^{-1} = 2|u_{ij}|\mathbb{I}_2 \ge 0$, where $\tilde{P} = \begin{pmatrix} c & 0 \\ -u_{ij} & 1 \end{pmatrix}$ is a transfer matrix that symmetrises $\nabla F$.

Hence the final form of the semi-discrete finite volume scheme

$$ \partial_t U_i + \frac{A^{-}(U_i, U_{i+1})(U_{i+1} - U_i) + A^{+}(U_{i-1}, U_i)(U_i - U_{i-1})}{\Delta x} = 0. \qquad (6) $$

The matrices $A^+$ and $A^-$ are defined by:

$$ A^{+}(U_i, U_j) = \frac{A(U_i, U_j) + D_{\text{stag}}(U_i, U_j)}{2} = \begin{pmatrix} \frac{|u_{ij}| - u_{ij}}{2} & 1 \\ -u_{ij}^2 & u_{ij} + \frac{|u_{ij}| + u_{ij}}{2} \end{pmatrix}, \qquad (7) $$
$$A^{-}(U_i, U_j) = \frac{A(U_i, U_j) - D_{\text{Stag}}(U_i, U_j)}{2} = \begin{pmatrix} -\frac{|u_{ij}| - u_{ij}}{2} & 0 \\ c^2 & u_{ij} - \frac{|u_{ij}| + u_{ij}}{2} \end{pmatrix}, \quad (8)$$

where $A(U_i, U_j)$ is the Roe matrix between the states $U_i$ and $U_j$.

In order to analyse the total variation of schemes (6), we study the behaviour of Riemann invariants that are connected to the left eigenvectors of $\nabla F$ (see [7, chapter 1.5]). Indeed for 1D smooth flows, the Riemann invariants $l_1, l_2$ follow a transport equation $\partial_t l_k + \lambda_k \partial_x l_k$, and therefore should follow a TVD property : $\partial_t TV(l_k)(t) = 0$, $k = 1, 2$ (see Kruzkov's theorem [7, section 3.3]).

### 3 Evolution of the discrete Riemann invariants

We linearise the Euler system (2) around a constant state $U_0 = (\rho_0, q_0)$:

$$\partial_t U + \nabla F(U_0) \partial_x U = 0, \quad (9)$$

as well as the (pseudo-)staggered scheme (6) which becomes

$$\partial_t U_i + \frac{A^{-}(U_0, U_0)(U_{i+1} - U_i) + A^{+}(U_0, U_0)(U_i - U_{i-1})}{\Delta x} = 0, \quad (10)$$

$\nabla F(U_0)$ is diagonalisable with left eigenvectors $L_1(U_0) = (u_0 + c, -1)$, and $L_2(U_0) = (-u_0 + c, 1)$. The transition matrix $P(U_0)$ for diagonalising $A(U_0, U_0)$ and its inverse $P^{-1}(U_0)$ are given by

$$P(U_0) = \begin{pmatrix} u_0 + c & -1 \\ -u_0 + c & 1 \end{pmatrix}, \quad P^{-1}(U_0) = \frac{1}{2c} \begin{pmatrix} 1 & 1 \\ u_0 - c & u_0 + c \end{pmatrix}. \quad (11)$$

In the linear case, the Riemann invariants are:

$$l_1(U_0, U) = \rho u_0 + \rho c - q, \quad l_2(U_0, U) = -\rho u_0 + \rho c + q. \quad (12)$$

We compute the evolution of the discrete invariants by multiplying (10) by the transition matrix $P(U_0)$,

- $U_i$ becomes $P(U_0)U_i = (l_{1,i}, l_{2,i})$,

- $A(U_0, U_0)$ becomes $P(U_0)A(U_0, U_0)P^{-1}(U_0) = \begin{pmatrix} u_0 - c & 0 \\ 0 & u_0 + c \end{pmatrix}$,

- $D_{\text{Stag}}(U_0, U_0)$ becomes

$$\bar{D}_{\text{Stag}}(U_0, U_0) = P(U_0)D_{\text{Stag}}(U_0, U_0)P^{-1}(U_0) = \begin{pmatrix} |u_0| & c \\ -c & |u_0| \end{pmatrix}, \quad (13)$$

- $A^+(U_0, U_0)$ becomes $\frac{1}{2} \begin{pmatrix} u_0 + |u_0| - c & c \\ -c & u_0 + |u_0| + c \end{pmatrix}$,
$$ - A^{-}(U_{0}, U_{0}) \text{ becomes } \frac{1}{2} \begin{pmatrix} u_{0} - |u_{0}| - c & -c \\ c & u_{0} - |u_{0}| + c \end{pmatrix}. $$

Hence the evolution of the discrete Riemann invariants is given by

$$ \partial_t l_{1,i} + (u_0^- - \frac{c}{2}) \frac{[l_1]_i^{i+1}}{\Delta x} - \frac{c}{2} \frac{[l_2]_i^{i+1}}{\Delta x} + (u_0^+ - \frac{c}{2}) \frac{[l_1]_{i-1}^i}{\Delta x} + \frac{c}{2} \frac{[l_2]_{i-1}^i}{\Delta x} = 0 \quad (14) $$

$$ \partial_t l_{2,i} + \frac{c}{2} \frac{[l_1]_i^{i+1}}{\Delta x} + (u_0^- + \frac{c}{2}) \frac{[l_2]_i^{i+1}}{\Delta x} - \frac{c}{2} \frac{[l_1]_{i-1}^i}{\Delta x} + (u_0^+ + \frac{c}{2}) \frac{[l_2]_{i-1}^i}{\Delta x} = 0. \quad (15) $$

where $[X]_j^i$ is the jump between cells j and i of the X variable: $[X]_j^i = X_i - X_j$.

## 4 Low Mach TVD analysis

The following theorem on the low Mach TVD analysis follows from the fact that for large sound speeds ($c \gg |u|$) we have $|l_1 - l_2| \ll \min\{l_1, l_2\}$.

**Theorem 1 (Low Mach TVD analysis).** *The numerical scheme (3)-(4) based on the diffusion matrix (5) satisfies*

$$ \partial_t \sum_i |[l_1]_i^{i+1}| \le 0 $$

provided

- assumption 1: $c$ (resp. $M$) is large enough (resp. small enough)

- assumption 2: $[\rho]_i^{i+1} = 0$ implies $2|[q]_i^{i+1}| + u^+|[\rho]_{i-1}^i| - [q - \rho u_0]_{i-1}^i \ge 0$.

*Proof.* We focus on $TV(l_1)$, but the same approach can be used to obtain $TV(l_2)$. We first compute the evolutions of $[l_1]_i^{i+1}$ from (14) for two cells i and i + 1.

$$ \begin{aligned} \partial_t [l_1]_i^{i+1} + (u^- - \frac{c}{2}) \frac{[l_1]_{i+1}^{i+2} - [l_1]_i^{i+1}}{\Delta x} + (u^+ - \frac{c}{2}) \frac{[l_1]_i^{i+1} - [l_1]_{i-1}^i}{\Delta x} \\ - \frac{c}{2} \frac{[l_2]_{i+1}^{i+2} - [l_2]_i^{i+1}}{\Delta x} + \frac{c}{2} \frac{[l_2]_i^{i+1} - [l_1]_{i-1}^i}{\Delta x} = 0. \end{aligned} \quad (16) $$

We multiply (16) by $sign([l_1]_i^{i+1})$ to obtain the evolution of $|[l_1]_i^{i+1}|$ and sum over i to obtain the evolution of the variations:

$$ \begin{aligned} \partial_t \sum_i |[l_1]_i^{i+1}| + (u^- - \frac{c}{2}) \frac{\sum_i \operatorname{sign}([l_1]_i^{i+1}) [l_1]_{i+1}^{i+2} - \sum_i |[l_1]_i^{i+1}|}{\Delta x} \\ + (u^+ - \frac{c}{2}) \frac{\sum_i |[l_1]_i^{i+1}| - \sum_i \operatorname{sign}([l_1]_i^{i+1}) [l_1]_{i-1}^i}{\Delta x} \\ - \frac{c}{2} \sum_i \operatorname{sign}([l_1]_i^{i+1}) \frac{[l_2]_{i+1}^{i+2} - [l_2]_i^{i+1}}{\Delta x} + \frac{c}{2} \sum_i \operatorname{sign}([l_1]_i^{i+1}) \frac{[l_2]_i^{i+1} - [l_2]_{i-1}^i}{\Delta x} = 0. \end{aligned} $$
We perform the shift of indices $i+1 \to i$ and $i-1 \to i$, and after some factori-
sation we get:

$$
\begin{align*}
& \partial_t \sum_i |[l_1]_i^{i+1}| + (u^- - \frac{c}{2}) \sum_i |[l_1]_i^{i+1}| \frac{\operatorname{sign}([l_1]_{i-1}^i)\operatorname{sign}([l_1]_i^{i+1}) - 1}{\Delta x} \\
& \qquad + (u^+ - \frac{c}{2}) \sum_i |[l_1]_i^{i+1}| \frac{1 - \operatorname{sign}([l_1]_{i+1}^{i+2})\operatorname{sign}([l_1]_i^{i+1})}{\Delta x} \\
& - \frac{c}{2} \sum_i [l_2]_i^{i+1} \frac{\operatorname{sign}([l_1]_{i-1}^i) - \operatorname{sign}([l_1]_i^{i+1})}{\Delta x} + \frac{c}{2} \sum_i [l_2]_i^{i+1} \frac{\operatorname{sign}([l_1]_i^{i+1}) - \operatorname{sign}([l_1]_{i+1}^{i+2})}{\Delta x} = 0.
\end{align*}
$$

We now expand $u^{\pm} - \frac{c}{2}$ and introduce the difference $[l_1]_i^{i+1}$ in the last two lines:

$$
\begin{align*}
& \partial_t \sum_i |[l_1]_i^{i+1}| - \frac{c}{\Delta x} \sum_i |[l_1]_i^{i+1}| (\operatorname{sign}([l_1]_i^{i+1})\operatorname{sign}([l_1]_{i-1}^i) - 1) \\
& \quad + \sum_i (u^- |[l_1]_i^{i+1}| - \operatorname{sign}([l_1]_i^{i+1}) \frac{c}{2} |[l_1 - l_2]_i^{i+1}|) \frac{\operatorname{sign}([l_1]_{i-1}^i)\operatorname{sign}([l_1]_i^{i+1}) - 1}{\Delta x} \\
& \quad + \sum_i (u^+ |[l_1]_i^{i+1}| - \operatorname{sign}([l_1]_i^{i+1}) \frac{c}{2} |[l_1 - l_2]_i^{i+1}|) \frac{1 - \operatorname{sign}([l_1]_{i+1}^{i+2})\operatorname{sign}([l_1]_i^{i+1})}{\Delta x} = 0.
\end{align*}
$$

Performing a shift of indices $i \to i-1$ in the last line we find

$$
\begin{align*}
& \partial_t \sum_i |[l_1]_i^{i+1}| - \frac{c}{\Delta x} \sum_i |[l_1]_i^{i+1}| (\operatorname{sign}([l_1]_i^{i+1})\operatorname{sign}([l_1]_{i-1}^i) - 1) \\
& \quad + \sum_i (u^{-} |[l_1]_i^{i+1}| - \operatorname{sign}([l_1]_i^{i+1}) \frac{c}{2} |[l_1 - l_2]_i^{i+1}|) \frac{\operatorname{sign}([l_1]_{i-1}^i)\operatorname{sign}([l_1]_i^{i+1}) - 1}{\Delta x} \\
& \quad + \sum_i (u^{+} |[l_1]_{i-1}^i| - \operatorname{sign}([l_1]_{i-1}^i) \frac{c}{2} |[l_1 - l_2]_{i-1}^i|) \frac{1 - \operatorname{sign}([l_1]_i^{i+1})\operatorname{sign}([l_1]_{i-1}^i)}{\Delta x} = 0,
\end{align*}
$$

and therefore

$$
\begin{align}
& \partial_t \sum_i |[l_1]_i^{i+1}| - \frac{1}{\Delta x} \sum_i (\operatorname{sign}([l_1]_i^{i+1})\operatorname{sign}([l_1]_{i-1}^i) - 1) \nonumber \\
& \quad \left( c |[l_1]_i^{i+1}| - u^{-} |[l_1]_i^{i+1}| + u^{+} |[l_1]_{i-1}^i| + \operatorname{sign}([l_1]_i^{i+1}) \frac{c}{2} |[l_1 - l_2]_i^{i+1}| \right. \nonumber \\
& \quad \left. - \operatorname{sign}([l_1]_{i-1}^i) \frac{c}{2} |[l_1 - l_2]_{i-1}^i| \right) = 0. \tag{17}
\end{align}
$$

In order to complete the proof we need to prove that the second line in (17)
is positive. Note that the first three terms in (17) are positive. Thus there is
competition between the respective jumps of $l_1$ and $l_1 - l_2$ in absolute value.

Now consider the limit in (17) when $c \to \infty$. From the definition of the
Riemann invariants (12)

$$
[l_1]_i^{i+1} = c[\rho]_i^{i+1} + [q - \rho u_0]_i^{i+1} \quad (18)
$$

$$
[l_1 - l_2]_i^{i+1} = 2[q - \rho u_0]_i^{i+1}, \quad (19)
$$
there are two cases :

- if $[\rho]_i^{i+1} \neq 0$ then $[l_1 - l_2]_i^{i+1}$ and $[l_1 - l_2]_{i-1}^i$ are negligible (since order one in $c$) compared to $[l_1]_i^{i+1}$ (order two in $c$) as $c \to \infty$, hence the assumption 1

- if $[\rho]_i^{i+1} = 0$ then $[l_1]_i^{i+1} = [q]_i^{i+1}$, $[l_1 - l_2]_i^{i+1} = 2[q]_i^{i+1}$ and the second line in (17) reduces to

$$ \begin{aligned} & c|[q]_i^{i+1}| - u^{-}|[q]_i^{i+1}| + u^{+}|[l_1]_{i-1}^i| + \operatorname{sign}([q]_i^{i+1})c[q]_i^{i+1} \\ & - \operatorname{sign}([l_1]_{i-1}^i)c[q - \rho u_0]_{i-1}^i \\ & = 2c|[q]_i^{i+1}| - u^{-}|[q]_i^{i+1}| + u^{+}|[l_1]_{i-1}^i| - \operatorname{sign}([l_1]_{i-1}^i)c[q - \rho u_0]_{i-1}^i \end{aligned} $$

The leading term for large $c$ is $2|[q]_i^{i+1}| + u^{+}|[\rho]_{i-1}^i| - [q - \rho u_0]_{i-1}^i$, hence the assumption 2.

<figure><img src="image_1.png" /></figure>

*Remark 1.* The assumption 1 in theorem 1 is met for flows at low Mach number although we do not have a specific bound on the Mach number to be met in practice. The assumption 2 is more tricky as it has to be checked for every numerical method.

## 5 Numerical results

As an illustration of our main result, we present the numerical simulation of a Riemann problem for a perfect gas with isothermal gas law $p = \rho c^2$ with $c = 300m/s$.

The initial data is made of two constant states separated by a discontinuity : a left state : ($\rho_L = 1, q_L = 10\rho_L$), and a right state : ($\rho_R = 2, q_R = -10\rho_L$).

We compare the numerical result obtained with three fully implicit finite volume methods taking the form (3) where the numerical flux takes the form

$$ F_{i,i+1}(U_i, U_{i+1}) = \frac{F(U_i) + F(U_{i+1})}{2} + D(U_i, U_{i+1}) \frac{U_i - U_{i+1}}{2}, \quad (20) $$

where the matrix $D$ can be $D_{upwind} = |\nabla F|$ for the upwind scheme, $D_{centered} = 0$ for the centered scheme, or $D = D_{stag}$ (equation 5) for the staggered scheme.

Picture 1 shows the numerical results obtained with CFL number 1,

and picture 2a shows the evolution of the Riemann invariant $l_1 = \frac{u}{c} + \ln(\rho)$.

Picture 1 shows as expected that the three conservative schemes capture the three constant states separated by two waves. The result of the upwind and staggered schemes are very similar. However the centered scheme oscillates since it is not TVD.

On picture 2a we observe as expected that the total variation of the centered scheme grows to a large final value. The total variation of the upwind scheme remains almost constant. The total variation of the staggered scheme is very close to that of the upwind scheme, however it increases in the first time step.
<figure><img src="image_2.png" /><figcaption>Fig. 1: Solution of the Riemann problem at $t = \frac{nx}{4}$</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(a) Evolution of the total variation
of the Riemann invariant $l_1$</figcaption></figure>

<caption>(b) Evolution of the total variation
of the Riemann invariant $l_1$</caption>
With the initial data chosen, the second assumption of theorem 1 is satisfied. We therefore expected that the scheme would become TVD as one increases the sound speed. However looking at picture 2b where we increased the sound speed from 30 to 30000, the total variation of the pseudo-staggered scheme never appears to be diminishing. One possible explanation is the fact that we analysed the linearised Euler system but simulated and plotted the nonlinear invariant. Another possible reason is that our theoretical analysis is based on a continuous in time semi-discrete equation whereas the numerical results are obtained with a fully discrete scheme.

## 6 Conclusion

We performed a TVD analysis for a stable conservative scheme inspired from staggered discretisations used in the multiphase flow community. The Riemann invariants were analysed and their total variation decreases provided the two conditions of theorem 1 are met. These assumptions are not necessarily met in practice, and given the numerical results, we believe our TVD criteria is not optimal. We are therefore working on improving our TVD analysis in order to match the observed behaviour of the scheme.

## References

1. Ndjinga, M., Ait-Ameur, K.: A New Class of L2-Stable Schemes for the Isentropic Euler Equations on Staggered Grids, Proceedings of Finite Volumes for Complex Applications - Methods, Theoretical Aspects, Examples: FVCA 9, Bergen, Norway, June 2020

2. Harlow, F., Amsden, A.: Numerical calculation of almost incompressible flow. Journal of Computational Physics, 3:8093 (1968)

3. Hirt, C.W.: Heuristic stability theory for finite difference equations. J. Comp. Phys., 2, pp. 339-355 (1968)

4. Prosperetti, A., Tryggvason, G.: Computational methods for multiphase flow. Cambridge University press (2009)

5. Harlow, F.H., Welch, J.E.: Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. The Physics of Fluids 8(12), 2182-2189 (1965)

6. Gallouët, T., Gastaldo, L., Herbin, R., Latché, J.-C.: An unconditionally stable pressure correction scheme for the compressible barotropic Navier-Stokes equations, Math. Model. Numer. Anal., 42, 2, 303-331, 2008

7. Godlewski, E., Raviart, P.A.: Numerical approximation of hyperbolic systems of conservation laws, Applied Mathematical Sciences, vol. 118., Springer-Verlag, New York (1996)