# On the convergence of Godunov scheme with a centered discretization of the pressure gradient

Jonathan Jung, Ibtissem Lannabi, Vincent Perrier

► To cite this version:

Jonathan Jung, Ibtissem Lannabi, Vincent Perrier. On the convergence of Godunov scheme with a centered discretization of the pressure gradient. FVCA X 2023 - Finite Volumes for Complex Applications X, Oct 2023, Strasbourg, France. hal-04260535

HAL Id: hal-04260535

https://hal.science/hal-04260535v1

Submitted on 26 Oct 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.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On the convergence of Godunov scheme with a centered discretization of the pressure gradient

Jonathan Jung*†

Ibtissem Lannabi†‡

Vincent Perrier‡†

October 24, 2023

## Abstract

This paper deals with the numerical resolution of a linear wave system using the Godunov scheme with a centered discretization of the pressure gradient. The interest in such schemes is motivated by the low Mach number accuracy problem. We have shown that for both steady and unsteady flows, an oscillatory mode appears in the numerical solution. This can be explained by the loss of the Total Variation Diminishing property on the characteristic variables. Moreover, we have illustrated numerically that the long time numerical solution does not converge to the expected steady state. In addition, the existence of the oscillatory mode in the numerical solution jeopardizes the mesh convergence rate of the scheme.

## Introduction

It is well-known that explicit Godunov-type schemes for solving the Euler equations fail, in general, to be accurate at low Mach number [7], in particular, on Cartesian grids. The most common argument to explain this failure is based on an asymptotic analysis of both the continuous equations and the numerical scheme. Initially, the low Mach number accuracy problem was explained via a single-scale asymptotic analysis that takes into account only the convective waves [7]. Afterwards, both acoustic and convective waves were taken into account by performing a two-scale asymptotic analysis [10, 8, 2], which led to a coupling of the zeroth-order momentum and the first-order pressure through a wave system.

A natural question that arises when considering a discretization of a Cauchy-problem for a wave system which has a long-time limit, is whether a discretization of this system also has a long-time limit and whether this long-time limit is consistent with that of the continuous system [9]. In [8], a link was elaborated between the low Mach number accuracy problem and the long-time consistency problem of a discrete wave system: a Godunov-type scheme is low Mach number accurate if and only if the discrete linear wave system derived from the two-scale asymptotic analysis is long-time limit consistent.

To ensure accuracy at low Mach number for the numerical resolution of the Euler equations on Cartesian grids, several fixes have been proposed. In this paper, we are interested in fixes of type [3, 11]. These schemes are also known to be two-scale asymptotic consistent with a discretization of the wave system, in which a centered discretization of the pressure gradient is used [2]. These schemes will be called pressure-centered type schemes. Since these fixes may introduce other issues, for instance, the appearance of an oscillatory mode (or checkerboard mode) in the numerical solution, it is interesting to address the study of the long-time limit of the wave system. Consider the linear wave system

$$ \left\{ \begin{aligned} \partial_{\tau} \hat{p} + \frac{1}{\rho} \nabla_{\mathbf{x}} \cdot \hat{\mathbf{u}} &= 0, \\ \partial_{\tau} \hat{\mathbf{u}} + \hat{\kappa} \nabla_{\mathbf{x}} \hat{p} &= 0, \end{aligned} \right. \qquad (1) $$

solved on a domain $\Omega$, associated with the initial condition $\hat{p}(\tau = 0, \mathbf{x}) = \hat{p}^0(\mathbf{x})$ and $\hat{\mathbf{u}}(\tau = 0, \mathbf{x})^T = \hat{\mathbf{u}}^0(\mathbf{x})^T$, where $\hat{p}$ denotes the pressure, $\hat{\mathbf{u}}$ denotes the velocity field, $\hat{\kappa}$, $\hat{\rho} > 0$ and $\hat{c} = \sqrt{\hat{\kappa}/\hat{\rho}}$ is the sound velocity. Let $T_h = \bigcup_i \Omega_i$ denotes a mesh of the domain $\Omega$, on which a cell-centered finite volume discretization is applied. We denote by $|\Omega_i|$ the volume of the cell $\Omega_i$, by $|\Gamma_{ij}|$ the size of the side $\Gamma_{ij}$ linking the cell $\Omega_i$ to the cell $\Omega_j$, by $\mathbf{n}_{ij}$ the unit normal outgoing from the cell $\Omega_i$ and going to the cell $\Omega_j$, by $\mathcal{V}_{\text{int}}(i)$ the set of the cells sharing a face with the cell $\Omega_i$, by $\mathcal{V}_b(i)$ the set of the boundary sides of $\Omega_i$. We consider the numerical scheme

$$ \begin{aligned} |\Omega_i| \left( \frac{\hat{\mathbf{U}}_i^{n+1} - \hat{\mathbf{U}}_i^n}{\delta\tau} \right) + \sum_{j \in \mathcal{V}_{\text{int}}(i)} |\Gamma_{ij}| \mathbf{F}^\theta (\hat{\mathbf{U}}_i^n, \hat{\mathbf{U}}_j^n, \mathbf{n}_{ij}) \\ + \sum_{j \in \mathcal{V}_b(i)} |\Gamma_{ij}| \mathbf{F}^b (\hat{\mathbf{U}}_i^n, \hat{\mathbf{U}}_b, \mathbf{n}_{ij}) = \mathbf{0}, \end{aligned} \qquad (2) $$

<footnote>*CNRS/Univ Pau and Pays Adour/E2S UPPA, Laboratoire de Mathématiques et de leurs Applications de Pau - Fédération IPRA, UMR5142 64000, Pau, France; jonathan.jung@univ-pau.fr & ibtissem.lannabi@univ-pau.fr & vincent.perrier@inria.fr</footnote>

<footnote>†Cagire team, Inria Bordeaux Sud-Ouest, France</footnote>
where $\delta\tau$ denotes the time-step, $\hat{\mathbf{U}}_i^n = (\hat{p}_i^n, \hat{\mathbf{u}}_i^n)^T$ is the vector of the discretized variables, $\hat{\mathbf{U}}_b = (\hat{p}_b, \hat{\mathbf{u}}_b)^T$ is a state weakly imposed on the boundary of the domain $\Omega$, $\mathbf{F}^\theta(\hat{\mathbf{U}}_i^n, \hat{\mathbf{U}}_j^n, \mathbf{n}_{ij})$ is the numerical flux inside the domain and is given by

$$ \mathbf{F}^\theta(\hat{\mathbf{U}}_i^n, \hat{\mathbf{U}}_j^n, \mathbf{n}_{ij}) = \begin{pmatrix} \frac{1}{\hat{\rho}} \frac{\hat{\mathbf{u}}_i^n + \hat{\mathbf{u}}_j^n}{2} \cdot \mathbf{n}_{ij} + \frac{\hat{c}}{2} (\hat{p}_i^n - \hat{p}_j^n) \\ \hat{\kappa} \frac{\hat{p}_i^n + \hat{p}_j^n}{2} \mathbf{n}_{ij} + \theta \frac{\hat{c}}{2} ((\hat{\mathbf{u}}_i^n - \hat{\mathbf{u}}_j^n) \cdot \mathbf{n}_{ij}) \mathbf{n}_{ij} \end{pmatrix}, \quad (3) $$

and $\mathbf{F}^b(\hat{\mathbf{U}}_i^n, \hat{\mathbf{U}}_b, \mathbf{n}_{ij})$ is the numerical flux on the boundary of $\Omega$. The parameter $\theta \in \{0, 1\}$ is introduced in order to consider both the Godunov ($\theta = 1$) and the pressure-centered ($\theta = 0$) schemes. Note that (2) with (3) is $L^2$ stable for $\theta = 1$ with CFL 1 and for $\theta = 0$ with CFL 0.5, see [3]. In section 1, periodic boundary conditions will be considered, so $V_b(i) = \emptyset$ for all $i$. In section 2, wall and Steger-Warming boundary conditions will be considered. We denote by $V_{\text{wall}}(i)/V_{\text{SW}}(i)$ the set of the boundary sides of $\Omega_i$ on which wall/Steger-Warming boundary conditions are imposed: $V_b(i) = V_{\text{wall}}(i) \cup V_{\text{SW}}(i)$. The flux $\mathbf{F}^b(\hat{\mathbf{U}}_i, \hat{\mathbf{U}}_b, \mathbf{n}_{ij})$ is then given by

$$ \mathbf{F}^b(\hat{\mathbf{U}}_i, \hat{\mathbf{U}}_b, \mathbf{n}_{ij}) = \begin{cases} \begin{pmatrix} 0 \\ \hat{\kappa}\hat{p}_i\mathbf{n}_{ij} + \hat{c}(\hat{\mathbf{u}}_i \cdot \mathbf{n}_{ij})\mathbf{n}_{ij} \end{pmatrix}, & \text{if } j \in V_{\text{wall}}(i), \\ \mathbf{F}^{\theta=1}(\hat{\mathbf{U}}_i, \hat{\mathbf{U}}_b, \mathbf{n}_{ij}), & \text{if } j \in V_{\text{SW}}(i). \end{cases} \quad (4) $$

In section 1, we show that the Total Variation Diminishing (TVD) property is lost with the pressure-centered scheme. In section 2, a problem in which the continuous wave system has a long-time limit is addressed. We show that the long-time limit of the numerical solution contains an oscillatory mode, which we are able to isolate. Lastly, we show that the long-time limit consistency is lost.

# 1 Loss of the TVD property

In this section, we are interested in the one-dimensional first-order wave system (1), solved on an infinite domain $\Omega$, using a finite volume scheme of the form (2)-(3). Let us recall that the TV norm of $v = (v_i)_{i \in \mathbb{N}}$ is given by $\text{TV}(v) = \sum_j |v_{j+1} - v_j|$. In what follows, the characteristic variables $\frac{\hat{p}}{2} - \frac{\hat{u}}{2\hat{\rho}\hat{c}}$ and $\frac{\hat{p}}{2} + \frac{\hat{u}}{2\hat{\rho}\hat{c}}$ are denoted by $\hat{C}_-$ and $\hat{C}_+$, respectively. The following proposition is presented.

**Proposition 1.** (Total Variation Diminishing property)

* The explicit Godunov scheme (2)-(3) with $\theta = 1$ is TVD on the characteristic variables $\hat{C}_-$ and $\hat{C}_+$ under the CFL condition $0 \le \hat{c}\delta\tau/\delta x \le 1$, i.e

$$ \text{TV}(\hat{C}_{-}^{n+1}) \le \text{TV}(\hat{C}_{-}^{n}), \quad \forall \hat{C}_{-}^{n}. $$

* The explicit pressure-centered scheme (2)-(3) with $\theta = 0$ is not TVD, i.e, for all $\delta\tau \in ]0, \delta x/(2\hat{c})[$, there exists an initial condition $(\hat{p}^0, \hat{u}^0)$ for which

$$ \text{TV}(\hat{C}_{-}^{1}) > \text{TV}(\hat{C}_{-}^{0}). $$

*Proof.* Regarding the conservation of the TVD property for the explicit Godunov scheme, it follows immediately since the Godunov scheme written in terms of the characteristic variables is an upwind scheme. We refer to [4], for a detailed proof.

Regarding the loss of the TVD property for the explicit pressure-centered scheme, it is obtained by considering the following initial condition: $\hat{p}^0/2 = \hat{u}^0/(2\hat{\rho}\hat{c})$ where $\hat{u}^0$ will be taken such that it has a discontinuity at some point. For example, we take the discrete initial condition defined by $\hat{u}_j^0 = -1$ for $j > k$ and $\hat{u}_j^0 = 1$ for $j \le k$. We thus obtain

$$ \begin{aligned} \text{TV}(\hat{C}_{-}^{1}) &= \frac{\hat{c}\delta t}{2\delta x} \sum_{j} \left| \frac{\hat{u}_{j+2}^{0} - 3\hat{u}_{j+1}^{0} + 3\hat{u}_{j}^{0} - \hat{u}_{j-1}^{0}}{2\hat{\rho}\hat{c}} \right| = \frac{\delta t}{\hat{\rho}\delta x} \left| \hat{u}_{k}^{0} - \hat{u}_{k+1}^{0} \right| \\ &= \frac{2\delta t}{\hat{\rho}\delta x} > 0 = \text{TV}(\hat{C}_{-}^{0}), \end{aligned} $$

which ends the proof. $\square$
In order to numerically illustrate the loss of the TVD property we consider the one dimensional Riemann test case for unsteady flows, associated with the first-order wave system (1). The initial condition considered for the selected test case is the following

$$ \hat{p}_L^0 = \hat{u}_L^0 = 1, \quad \hat{p}_R^0 = \hat{u}_R^0 = -1, \qquad (5) $$

with final time equal to 0.1, $\hat{c} = \hat{\rho} = 1$ and a CFL number equal to 0.45.

With the chosen initial condition, we remark that the characteristic variable $\hat{C}_-$ is initialized as zero. In Figure 1, we observe that with the Godunov scheme, the uniformity of $\hat{C}_-$ is preserved over time since the Godunov scheme is TVD on the characteristic variables. However, with the pressure-centered scheme, the uniformity of this characteristic variable is not preserved. Moreover, oscillations appear on the characteristic variables, as well as on the pressure and velocity. Note that the same phenomenon is observed when considering an initial uniform characteristic variable $\hat{C}_+$. In Figure 2, the plot of the TV norms is presented. It highlights that the TV norms of the two characteristic variables are preserved for the Godunov scheme but increase for the pressure-centered scheme. Similar results were obtained by considering an implicit time integration for both the Godunov and the pressure-centered schemes, but these results are not presented in this paper.

<figure><img src="image_2.png" /><figcaption>Figure 1: One-dimensional Riemann problem: characteristic variables at $t = 0.1$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: One-dimensional Riemann problem: TV norm of the characteristic variables as a function of time iterations.</figcaption></figure>
## 2 Loss of mesh convergence of the numerical solution in the long time limit

In this section, we are interested in the linear wave system (1), solved on a bounded domain $\Omega$, with wall and Steger-Warming boundary conditions. In this framework, the long-time limit of the solution of both the continuous wave system (1) and the numerical schemes (2)-(3)-(4) with $\theta \in \{0, 1\}$ exists (see [9]). This limit can be identified thanks to a Hodge-Helmholtz Decomposition (HHD). Indeed, at the continuous level, by considering a vector field $\hat{u} \in L^2(\Omega)$ defined on a bounded domain $\Omega$ with a smooth boundary $\partial\Omega$, $\hat{u}$ can be uniquely decomposed into a sum of an irrotational component (which is a gradient of a scalar potential function), and a divergence free one (which is a curl of a vector potential function). The proof for this decomposition in the case where $\hat{u}$ satisfies wall/Steger-Warming boundary conditions can be found in [8] and is based on solving a Poisson equation with Neumann boundary conditions. At the discrete level, a similar exact decomposition is obtained on triangular grids in [8], see also [1, 5, 4, 6] for other boundary conditions. By considering a piecewise constant vector field $\hat{u}_h$, its discrete HHD is given by $\hat{u}_h = P^\psi[\hat{u}_h] + P^\varphi[\hat{u}_h]$, where $P^\psi[\hat{u}_h]$ is the curl of a continuous $P^1$ potential vector, $P^\varphi[\hat{u}_h]$ is the gradient of a Crouzeix-Raviart scalar potential and $P^\psi[\hat{u}_h] \cdot n = \hat{u}_h \cdot n$ on the boundary (see [8, Proposition 5], for details about the proof). Back to the identification of the long-time limit, at the continuous level, it is identified as the uniform pressure $\hat{p}_b$ and the divergence free component of the HHD of the initial velocity $\hat{u}^0$. On triangular meshes, with the Godunov scheme ((2)-(3)-(4) with $\theta = 1$), the long-time limit of $\hat{p}_h$ and $\hat{u}_h$ is consistent with the one obtained at the continuous level: it corresponds to $\hat{p}_b$ and $P^\psi[\hat{u}_h^0]$. Since the long-time limit of the solution of the pressure-centered scheme ((2) with $\theta = 0$) has not been identified yet, our objective in this section is to characterize it numerically. Unfortunately, the Cartesian grid case cannot be addressed due to the non existence of a discrete HHD on such a grid (at least adapted to such boundary conditions).

The test case selected for this study is the scattering of a first-order wave system by a cylinder. The computational domain is an annulus $[r_0, r_1] \times [0, 2\pi[$, where $r_0 = 0.5$ and $r_1 = 5.5$. The boundary conditions considered for this test are: a wall boundary condition on the internal circle of radius $r_0$ and a Steger-Warming boundary condition on the external circle of radius $r_1$ with $\hat{p}_b = 0$ and $\hat{u}_b = (1, 0)^T$. The initial condition is set to $\hat{p}^0 = 0$ and $\hat{u}^0 = (0, 0)^T$.

It is worth emphasizing that the exact steady solution is known [9] and is given by $\hat{p}^{\text{exact}}(r, \theta) = 0$ and $\hat{u}^{\text{exact}}(r, \theta) = \frac{r_1^2}{r_1^2 - r_0^2} \left( 1 - \frac{r_0^2}{r_1^2} \cos(2\theta), -\frac{r_0^2}{r_1^2} \sin(2\theta) \right)^T$.

The mesh type considered is an unstructured triangular mesh generated using GMSH, with a characteristic length of $2\pi r_0 / 4N$ on the internal circle and $2\pi r_1 / 4N$ on the external circle. The mesh selected for the next figures is produced with $N = 40$ and contains 1332 triangular cells. The test is run with a CFL number equal to 0.45 and $\hat{\rho} = \hat{\kappa} = 1$.

In what follows, numerical results presented correspond to a solution for which the pressure and velocity residuals have converged with a numerical precision of $10^{-14}$.

In Figure 3, the iso-contours of the $L^2$ norm of the long-time limit of the exact and numerical velocity fields are plotted. We observe that the numerical velocity field computed using the Godunov scheme is close to the exact velocity field, while the one obtained with the pressure-centered scheme is also close but polluted by spurious oscillations. Using the long-time consistency result proved in [8, 9] and recalled above, the long-time limit of the numerical velocity field $\hat{u}_h^\infty$ computed using the Godunov scheme is $P^\psi[\hat{u}_h^0]$, the divergence free component of the discrete HHD of the initial velocity<fnref n="1" /> $\hat{u}_h^0$. If $\hat{u}_h^\infty$ is the long-time limit of the numerical velocity field, we define

$$ \hat{u}_h^\infty, \text{Spurious} = \hat{u}_h^\infty - P^\psi[\hat{u}_h^0]. \qquad (6) $$

The difference (6) is plotted in Figure 4 for the Godunov and the pressure-centered schemes. With the Godunov scheme, the difference is zero, as expected, while with the pressure-centered scheme, (6) is a spurious oscillatory mode and its scale is of order $10^{-1}$, which is not negligible compared to the magnitude of the velocity field on the boundary which is of magnitude 1. In order to study the long-time consistency of the scheme, a grid convergence study is performed. To do so, unstructured triangular meshes obtained with $N \in \{12, 20, 40, 80, 160, 320\}$ containing 96, 332, 1332, 5154, 20546, 83810 cells respectively, are considered. In Figure 5, the $L^2$ norm of the error between the exact solution and the long-time numerical solution computed using the Godunov and the pressure-centered schemes is plotted. We observe that with both schemes, the long-time limit of the pressure matches a uniform pressure, which is consistent with the long-time limit of the continuous pressure, as explained in [9]. However, regarding the long-time velocity field, the order of convergence is 1 with the Godunov scheme, while the pressure-centered scheme does not converge. This provides numerical evidence that the presence of the spurious mode jeopardizes the long-time limit consistency for the velocity.

<footnote n="1">¹$\hat{u}_h^0$ is equal to 0 for this test case (note that $P^\psi[\hat{u}_h^0]$ is not equal to 0 due to non homogeneous boundary conditions)</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 3: Numerical velocity field computed using the Godunov and the pressure-centered schemes: Exact (left), Godunov (middle), Pressure-centered (right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Difference (6) computed using the Godunov and the pressure-centered schemes: Godunov (left), pressure-centered (right).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: $L^2$ norm of the error between the exact solution and the long-time numerical solution computed using the Godunov and the pressure-centered schemes as a function of the mesh step.</figcaption></figure>

Conclusion

In this paper, the study of the numerical resolution of a first-order wave system using the pressure-centered scheme and explicit time integration is addressed. We have shown that, with this type of scheme,

• the Total Variation Diminishing property is lost.

• Oscillations appear in the unsteady numerical solution.

• The long-time numerical solution contains a spurious mode in the velocity field, which is an oscillatory mode.

• The mesh convergence is compromised due to the existence of this spurious mode.

Further study needs to be conducted for the pressure-centered scheme with implicit time integration. Moreover, the study was done only for triangular cell geometries, thus it needs to be extended to the quadrangular case.
References

[1] Arnold, D. N., Falk, R. S.: A uniformly accurate finite element method for the Reissner-Mindlin plate. SIAM J. Numer. Anal., 26:1276-1290 (1989). https://doi.org/10.1137/0726074

[2] Bruel, P., Delmas, S., Jung, J., Perrier, V.: A low Mach correction able to deal with low Mach acoustics. Journal of Computational Physics, 378, 723-759 (2019). https://doi.org/10.1016/j.jcp.2018.11.020

[3] Dellacherie, S.: Analysis of Godunov type schemes applied to the compressible Euler system at low Mach number. Journal of Computational Physics, 229(4), 978-1016 (2010). https://doi.org/10.1016/j.jcp.2009.09.044

[4] Dellacherie, S., Jung, J., Omnes, P., Raviart, P. A.: Construction of modified Godunov-type schemes accurate at any Mach number for the compressible Euler system. Mathematical Models and Methods in Applied Sciences, 26(13), 2525-2615 (2016). https://hal.science/file/index/docid/776629/filename/article_bas_mach_dellacherie_omnes_et_par.pdf

[5] Dellacherie, S., Omnes, P., Rieper, F.: The influence of cell geometry on the Godunov scheme applied to the linear wave equation. Journal of Computational Physics, vol. 229, no 14, p. 5315-5338 (2010). https://doi.org/10.1016/j.jcp.2010.03.012

[6] Dellacherie S., Jung J., Omnes P.: Construction of a low Mach finite volume scheme for the isentropic Euler system with porosity. ESAIM: Mathematical Modelling and Numerical Analysis, 1;55(3):1199-237 (2021). https://doi.org/10.1051/m2an/2021016

[7] Guillard, H., Viozat, C.: On the behaviour of upwind schemes in the low Mach number limit. Computers & fluids, 28(1), 63-86 (1999). https://doi.org/10.1016/S0045-7930(98)00017-6

[8] Jung, J., Perrier, V.: Steady low Mach number flows: identification of the spurious mode and filtering method. Journal of Computational Physics, 468, 111462 (2022). https://doi.org/10.1016/j.jcp.2022.111462

[9] Jung, J., Perrier, V.: Long time behavior of finite volume discretization of symmetrizable linear hyperbolic systems. IMA Journal of Numerical Analysis, 43(1), 326-356 (2023). https://doi.org/10.1093/imanum/drab092

[10] Müller, B.: Low-Mach-number asymptotics of the Navier-Stokes equations. Floating, Flowing, Flying: Pieter J. Zandbergen's Life as Innovator, Inspirator and Instigator in Numerical Fluid Dynamics, 97-109 (1998). https://doi.org/10.1007/978-94-017-1564-5_6

[11] Rieper, F.: A low-Mach number fix for Roe's approximate Riemann solver. Journal of Computational Physics, 230(13), 5263-5287 (2011). https://doi.org/10.1016/j.jcp.2011.03.025