Construction of a low Mach finite volume scheme for the isentropic Euler system with porosity

Classical ﬁnite volume schemes for the Euler system are not accurate at low Mach number and some ﬁxes have to be used and were developed in a vast literature over the last two decades. The question we are interested in in this article is: What about if the porosity is no longer uniform? We ﬁrst show that this problem may be understood on the linear wave equation taking into account porosity. We explain the inﬂuence of the cell geometry on the accuracy property at low Mach number. In the triangular case, the stationary space of the Godunov scheme approaches well enough the continuous space of constant pressure and divergence-free velocity, while this is not the case in the Cartesian case. On Cartesian meshes, a ﬁx is proposed and accuracy at low Mach number is proved to be recovered. Based on the linear study, a numerical scheme and a low Mach ﬁx for the non-linear system, with a non-conservative source term due to the porosity variations, is proposed and tested.


Introduction
In this paper, we are interested in low Mach compressible fluid flows in porous media. In industrial processes, porous media are used to simulate the flow in a nuclear reactor core. The porosity appears because there are section reductions in a nuclear reactor core. Moreover, if we want to simulate an accidental scenario, we sometimes need to take into account the compressibility effects. Another class of problems motivated by industrial consideration is the simulation of a gas flow across a grid. Since the grid is in general too small to be meshed, a homogenization process is used to model the interactions between the grid and the flow [41]. Then, we consider the barotropic Euler equation. Since the porosity is not constant, a non-conservative term appears in the equations during the homogenization process [4] and the equations write (1.1) In (1.1), ≥ 0 and x ∈ Ω are respectively time and space variables and (x) is the porosity. We suppose here that (x) is known and does not depend on time. Unknowns , u and ( ) are respectively the density, the velocity and the pressure of the fluid. The pressure law satisfies ′ ( ) > 0. System (1.1) is a non conservative hyperbolic system [27] with eigenvalues in direction n given by u · n − , u · n and u · n + . Studies of flows in a variable cross section duct consider the same model and variations of the cross section are modeled through (possibly discontinuous) changes in porosity.
The non-conservative term in (1.1) introduces mathematical and numerical difficulties. In [11], the authors give a mathematical sense to the non-conservative product and introduce some schemes, named well-balanced schemes, that solve correctly the non-conservative term [8,20]. The treatment of the non-conservative term is essential to preserve steady states solutions. In this paper, we propose a well-balanced scheme that exactly preserves steady solutions over time in one space dimension [21,26]. The proposed scheme is based on a VFRoe scheme, established in [18] for the shallow water equation with topography and derived for the Euler system with porosity in [39,40]. The VFRoe solver consists into a local linearization of a Riemann problem which is simpler to handle since it only deals with linear problems and avoids the complex exact resolution of the Riemann problem with porosity jump. This construction allows to easily build schemes that exactly preserve one dimensional steady states. To our knowledge, no generalization for purely multidimensional problems exists and the behavior of the numerical scheme with respect to multi-dimensional steady solutions must be studied on a case-by-case basis. In this paper, we focus on the behavior of the scheme in the low Mach limit and, as we will see, steady solutions of the numerical scheme will play a determining role.
Finite volume Godunov type schemes applied to the compressible Euler system with uniform porosity are known to be inaccurate at low Mach number [13,23]. Indeed, they do not allow to recover the incompressible limit as the Mach number tends to zero. Over the two last decades, a large amount of work has been dedicated to deriving fixes for the uniform porosity case: [5,7,12,13,16,23,24,28,29,33,37]. Some recent works have been done on low Mach fix for non-conservative systems, we refer to [2,3,45] for the Euler equation with gravity or to [34,35] for two-phases flows. In these last studies, flux preconditioning techniques, initially proposed by Turkel [46], are applied and quadrangular meshes are considered. Here, we propose to also study the behavior of the numerical scheme on triangular meshes. Indeed, it was shown that in the uniform porosity case, if the mesh is composed of triangles in 2D or tetrahedra in 3D, the accuracy at low Mach number with the Roe scheme is recovered [14,22,38]. To our knowledge, this is the first study on the behavior of classical schemes at low Mach number on triangular meshes for non-conservative systems.
In this article, we study the accuracy, at low Mach number, on triangular and Cartesian meshes, of a numerical scheme for the non conservative system (1.1). Since the accuracy problem appears also in the linear case, we base our study on the linear wave equation with porosity. The low Mach accuracy problem is then understood and fixed in the linear case for Cartesian meshes, and the reason for its correct behavior on triangular meshes is underlined. In particular, preliminary results obtained in [15] based on a modified equation approach are extended to the discrete Cartesian case. Based on the linear study, a well-balanced scheme accurate at low Mach number for the non-linear system (1.1) is proposed and numerical tests are performed. They confirm that both the non corrected and corrected schemes are able to recover the low Mach asymptotics on triangular meshes, while this is the case only for the corrected scheme on Cartesian meshes.

Low Mach limit
To study the behavior of system (1.1) at low Mach number, four characteristic scales are supposed to be known: a time scale 0 , a density scale 0 , a velocity scale 0 and a porosity scale 0 . Then, the following dimensionless variables are defined˜= It is natural to scale the length by 0 = 0 × 0 , the sound speed by 2 0 = ′ ( 0 ) and the pressure by 0 = 0 2 0 . If the corresponding dimensionless variables are used, system (1.1) reads

Formal asymptotic expansion when the Mach number goes to 0
We are interested in the solutions of (2.2) when → 0. We recall formally the theoretical results of [25] in order to take the porosity into account. All the variables The case˜( 0) = 0 is out of the scope of this paper.
By injecting these quantities in (2.2), the momentum equation at order −2 and −1 gives and then, since is a regular function of , this leads tõ Then, if the initial and boundary conditions are well prepared in the sense that and if on the domain boundary˜( 0) (resp.˜( 1) ) is uniformly and constantly equals to˜0 (resp. 0) and if ∫︀ Ω (˜ũ) (0) · n = 0, the solution of (2.2) satisfies Note that in order to obtain the second equation in (2.7), we have chosen˜0 = 1, which is always possible up to a change of density scale from 0 to 0˜0 . Equations (2.6) mean that at low Mach number, if the initial and boundary conditions are well prepared, the solution of the compressible Euler system with porosity (2.2) is close to the solution of the incompressible Euler equation with porosity (2.7). Results (2.6) are formally proven here. For a rigorous proof in the uniform porosity case, we refer to [30,42]. For classical finite volume schemes, relations (2.6) are not always satisfied at the discrete level: this is the so-called accuracy problem at low Mach number, which expresses that a spurious component˜( 1) ̸ = 0 could be introduced at the discrete level [23] due to numerical approximations. In the current contribution, we consider that a numerical scheme is accurate at low Mach number for system (1.1) if relations (2.6) are satisfied at the discrete level.

Wave equation with porosity
To study the low Mach behavior, we change the variables to symmetrize the problem.

Model
For this purpose, we set the reference sound speed to 1/ and we define (˜,x) such that where formally ≪ 1. By injecting (2.8) in (2.2), we obtain the system Linearizing around ( ,ũ) = (0, 0), taking into account that˜′ (˜0) = 1 when˜0 = 1 as explained above and simplifying the notation by removing all the·, we obtain the linear wave equation with porosity and ⋆ = 1.

Weighted incompressible space ℰ and acoustic space ℰ ⊥
We are interested in the properties of System (2.9) solved on a torus T ⊂ R ∈{1,2,3} with periodic boundary conditions. For this, we assume that is a periodic function on T and we define the weighted Hilbert space }︂ endowed with the scalar product Of course, the space 2 should not be mistaken for the acoustic operator . We also define the spaces 1 (T) and 2 (T) that are generalizations of 1 (T) and 2 (T) to weighted spaces. We note that since (x) ∈ [ min , 1] with min > 0, the functions and 1 are in ∞ (T), and we have 2 (T) = 2 (T), 1 (T) = 1 (T) and 2 (T) = 2 (T). Nevertheless, we keep the index to define these spaces to refer to the scalar product (2.10). At last, we define the space (2.11) When = 1, ℰ is named the incompressible space (see [13]). We have the following result: In other words, any = ( , u) ∈ 2 (T) 1+ can be decomposed into whereˆ= (ˆ,û) ∈ ℰ and ⊥ = ( ⊥ , u ⊥ ) ∈ ℰ ⊥ and this decomposition is unique and orthogonal with respect to the scalar product defined by (2.10).
We call ℰ ⊥ the acoustic space. This is a generalization of the Hodge decomposition. Decomposition (2.13) defines an orthogonal projection ↦ −→ P :=ˆ.

Properties of the linear wave equation with porosity
We now detail some properties of the linear wave equation with porosity. These properties will not be always satisfied in the discrete case. Lemma 2.3. Let ( , x) be the solution of (2.9) on T ⊂ R ∈{1,2,3} with initial condition 0 . Then: For all ∈ 2 (T) 1+ , we now define the energy := ⟨ , ⟩ . The following lemma is an extension of the energy conservation property of the classical linear wave equation: Lemma 2.4. Let ( , x) be the solution of (2.9) on T ⊂ R ∈{1,2,3} . Then, for all ≥ 0,

The low Mach asymptotics
With Lemma 2.3 and by linearity, we get that if ( , x) is the solution of (2.9) on T ⊂ R ∈{1,2,3} with initial condition 0 , then We note that since P 0 ∈ ℰ is a stationary solution of (2.9), then P = P 0 ; hence (2.15) can be written as In fact, (2.15) is a version of (2.6) for the linear case. Indeed, the left condition in (2.15) just means that the initial condition is well-prepared. In the non-linear case, the projection P in the incompressible space ℰ is replaced by the incompressible solution of (2.7).
In this article, we consider that a numerical scheme for the linear system (2.9) is accurate at low Mach number if (2.16) is satisfied at the discrete level. We will study this property on Cartesian and triangular meshes.

Godunov scheme for the linear wave equation with porosity and its kernels
In [14,16], we explained the satisfactory behavior of the Godunov scheme at low Mach number on triangular meshes and its wrong behavior on Cartesian rectangular meshes on the Euler system without porosity ( uniformly equal to 1) by studying the kernel of the discrete spatial operator associated to the Godunov scheme. We also remarked that the accuracy of the Godunov scheme at low Mach number on Cartesian meshes can be recovered by deleting the diffusion term on the velocity field in the Godunov scheme. In [15], we discussed the case with porosity with the help of the modified equation approach; the limitations of this approach is that it only gives hints (but does not provide with a complete proof) on what happens on Cartesian meshes, and does not apply to triangular meshes. Our aim here is to analyse the behavior of the schemes on triangular and rectangular Cartesian meshes by directly studying them rather than their modified equations.
We now recall the Godunov scheme for the linear wave equation with porosity, recall why the study of its kernel is so important to study its low Mach accuracy and compute explicitly its kernels on triangular and Cartesian rectangular meshes. In particular, we underline that the kernel is strongly linked to the numerical dissipation of the Godunov scheme.

Godunov scheme
Let us suppose that the domain T ⊂ R 2 is discretized by cells Ω . Let Γ be the common edge of the two neighboring cells Ω and Ω and n the unit vector normal to Γ pointing from Ω to Ω . We assume that the data , and the unknowns and u are defined on the cells Ω in the following way and then set ( ) = and ( u) = u . The semi-discrete Godunov scheme applied to the resolution of the linear wave equation is obtained by integrating (2.9) over each cell Ω and then solving a Riemann problem on each Γ to express interface fluxes as functions of cell-centered values. Details are provided in [15]. This results in with = 1 and where is a mean-value of on Γ which depends on ( , ) (e.g. arithmetic or harmonic mean). The numerical flux in (3.2) is non-conservative because of the term that multiplies the flux on the momentum equation. Moreover, it is easy to prove the following properties: Remark 3.1. The numerical scheme (3.1) is well-balanced in the sense that it preserves exactly the onedimensional steady states ( = cte, u = cte).
Remark 3.2. The numerical scheme (3.1) can also be viewed as the VFRoe scheme [6,17] obtained with the variables ( , , u) for system (2.9) where the linearized Riemann problem is solved considering that satisfies = 0.
Scheme (3.1) can be written in compact form where the subscript · ℎ recalls that (3.2) comes from a spatial discretization of (2.9).

The low Mach problem
We want to study whether the Godunov scheme is accurate at low Mach number in the sense that it satisfies a version of (2.16) at the discrete level. Then, discrete incompressible spaces ℰ ℎ and (︀ ℰ ℎ )︀ ⊥ and a discrete orthogonal projection P ℎ have to be defined on triangular or Cartesian meshes. Moreover, the key points to obtain (2.16) at the continuous level are that ℰ = Ker and that (2.9) conserves energy (see Lem. 2.4). Then, the relationship between the discrete incompressible space ℰ ℎ and the kernel of the Godunov scheme Ker L ℎ have to be studied. The following theorem explains why this study is so important: is a positive constant independent of the Mach number and suppose moreover that ℰ ℎ ⊆ Ker L ℎ . Then, we have For a proof, we refer to [13,16]. In Theorem 3.3, system (3.2) is assumed to be well-posed. In particular, stability will be studied in more details in Section 4. In the current section, we focus on the kernel of the Godunov scheme on Cartesian and triangular meshes.

Kernels of the Godunov scheme
We first study the discrete kernel of the Godunov scheme ( = 1 in (3.1)) on different types of meshes and of its low-Mach modification ( = 0) on Cartesian rectangular meshes. The kernel Ker L ℎ , of the discrete acoustic operator L ℎ , is defined by On any type of mesh we have the following result, whose proof is postponed to Appendix A: and

Kernel on a triangular mesh
We now study some particular properties of the behavior of the Godunov scheme on a triangular mesh. Especially, we study the relation between the kernel of the Godunov scheme on a triangular mesh and a discrete version of the space ℰ defined by (2.11).

Construction of ℰ ℎ,△ and
(︀ ℰ ℎ,△ )︀ ⊥ . We construct an accurate discrete version of the well-prepared subspace ℰ defined by (2.11). Let us suppose that all Ω are triangles arranged so that the computational domain is periodic. Moreover, let us denote by ℎ the standard 1 (first-order polynomial functions) Lagrange finite element space associated with this triangular mesh Let us also denote by ℎ the nonconforming Crouzeix-Raviart 1 finite element space associated with this triangular mesh ℎ := { ℎ ∈ 2 (T), ℎ periodic on T such that ∀Ω : ( ℎ ) |Ω ∈ 1 (Ω ) and ℎ is continuous at the edge midpoints}.
Note that since the functions in ℎ (resp. ℎ ) are 1 on each cell, their curls (resp. their gradients) are constant vectors on each cell. Let us also define the discrete vector subspace Then, we define the space of constant piecewise functions }︂ endowed with the scalar product (2.10) which may be written for ( ℎ ) 1 and ( ℎ ) 2 in 2 (T) 3 as Adapting the proof of Theorem 4.1 in [1] (see also [32]) to the case of periodic elements in ℎ and ℎ and weighted spaces, we may prove the following lemma: Lemma 3.5. Assume that (Ω ) =1,..., is a triangular periodic mesh of a rectangular domain with no internal holes. For any ( , u) ∈ R 3 , there exist unique ( , ) ∈ R 2 , a unique ℎ ∈ ℎ and a unique ℎ ∈ ℎ with . Moreover this decomposition is orthogonal for the scalar product (3.8).
Proof. We firstly prove the orthogonality of decomposition (3.9). The orthogonality between¯and −¯is obvious because, by definition of¯we have: Now, we prove the orthogonality for the decomposition of u. For any ( , ) ∈ R 2 and ℎ ∈ ℎ (then ∇ ℎ is a constant vector on each cell Ω ), we have: denotes the jump of ℎ through the edge Γ . To obtain the last equality, we used the fact that each interface Γ contributes twice in the sum over the cell boundaries. Since ℎ is a 1 function, its integral on the edge Γ is equal to the length |Γ | multiplied by the value of ℎ at its midpoint. Thus, since ℎ is continuous at the edge midpoints, we have ∫︀ Γ [ ℎ ] n d = 0 on any edge, which proves the orthogonality between the field 1 ( , ) and the gradient of any element in ℎ . Moreover, for any ℎ ∈ ℎ and ℎ ∈ ℎ (then ∇ × ℎ and ∇ ℎ are constant vectors on each cell Ω ), it holds that Since ∇ · (∇×) = 0, the second sum vanishes. Moreover, denoting by t a unit vector such that (n, t) is a direct orthonormal system, the equality (∇ × ℎ ) · n = (∇ ℎ ) · t and the fact that ∇ ℎ · t is continuous along any interface Γ (since ℎ ∈ ℎ is a 1 nodal Lagrange function) imply that But on Γ , the product ∇ ℎ · t [ ℎ ] is a 1 function, and its integral over Γ is equal to the length |Γ | multiplied by the value of this function at its midpoint. Thus, since ℎ is continuous at the midpoint, then ∫︀ Γ ∇ ℎ · t [ ℎ ] d = 0 on any edge, which proves orthogonality between 1 ∇ × ℎ and ∇ ℎ . Then, the orthogonality of the decomposition is proved.
We secondly prove the existence and the uniqueness of decomposition (3.9). For there is no difficulty. Thus, we only consider the decomposition for u. We have to prove that the function defined by . Firstly, we prove injectivity. As is a linear function, we just have to prove that Assume that for all ∈ 1, , 1 By the orthogonality that we proved above, this implies ∀ ∈ 1, : Since ℎ is continuous at the edge midpoints and since ℎ is continuous on T, ( ) =1,..., and ( ) =1,..., do not depend on . Then, we have Since ∫︀ T ℎ dx = 0, we obtain ℎ = 0. Since ℎ ( , ) = − + is periodic on T, we have = = 0 which implies that ℎ = . And since ∫︀ T ℎ dx = 0, we obtain ℎ = 0. The conclusion is that and the function in injective.
To prove surjectivity, we prove that dim Any function ℎ ∈ ℎ is completely and uniquely determined by its values at the independent nodes of the mesh, which implies that dim ℎ = . Moreover the vanishing mean-value of ℎ implies a constraint that links the values on the various nodes. Thus, we have dim 0 ℎ = − 1. On the other hand, any ℎ ∈ ℎ is completely and uniquely determined by its values at the independent edge midpoints of the mesh, then dim ℎ = . Moreover, the vanishing mean-value of ℎ implies a constraint that links the values on the various edges. Thus, we have dim 0 Now, in a triangular periodic mesh of a rectangular domain with no internal holes, it is well known that + = 2 (proof by recurrence on the number of cells using the Descartes-Euler formula for a periodic domain), which proves the bijectivity of the function .
Let us underline that Lemma 3.5 with Corollary 3.6 is the discrete version of Lemma 2.2 on a triangular mesh.
A first explanation of the satisfying behavior of the Godunov scheme on triangular meshes. Here, we prove that, on triangular meshes, the kernel of the Godunov scheme corresponds exactly to the discretized space ℰ △ . This property shows that the discrete stationary space discretizes well the continuous one. This gives a (partial) explanation of the satisfactory behavior of the Godunov scheme on a triangular mesh.
for some ( , , ℎ , ℎ ) ∈ R 2 × ℎ × ℎ , this decomposition being orthogonal. Thus, we just have to prove that ∇ ℎ = 0. By orthogonality, we have because u and ∇ ℎ are constant on each triangle Ω . Then, we can write Since ℎ ∈ Ker L ℎ,△ =1, , we have ( u) · n = ( u) · n and we denote by ( ) this common value. Thus where [ ℎ ] denotes the jump of ℎ through the edge Γ . As already explained in the proof of Lemma 3.5, the fact that ℎ is a 1 function, which is continuous at the edge midpoints implies that ∫︀ Γ [ ℎ ] d = 0 on any edge, which proves that is to say for all ∈ 1, , (∇ ℎ ) |Ω = 0. This proves that ℎ ∈ ℰ ℎ,△ .

Kernel on a Cartesian mesh
We now study some particular properties of the behavior of the Godunov scheme on a rectangular uniform Cartesian mesh. Especially, we study the relation between the kernels of the standard Godunov scheme ( = 1) and of its modification ( = 0) on a uniform Cartesian mesh and a discrete version of the space ℰ defined by (2.11).

Construction of ℰ ℎ, and
(︀ ℰ ℎ, )︀ ⊥ . We construct an accurate discrete version of the well-prepared subspace ℰ defined by (2.11). Suppose that the computational domain is a rectangle and that the mesh is made up of × rectangles of constant size ∆ ×∆ where and are the numbers of cells in the and directions. In what follows, we shall suppose that both and are odd. Indeed, if this is not the case, the situation is a little more involved due to even/odd decoupling which may produce checkerboard modes. We introduce the following operators, which are accurate approximations of their continuous counterparts: In these definitions, it is implicitly meant that ( , ) ∈ R is periodic, that is to say (3.11) Let us now define the following subspace, which is an accurate discrete version of ℰ defined by (2.11).
We shall also need the following weighted discrete scalar product: We introduce in the following lemma a discrete Hodge decomposition for a collocated Cartesian mesh with periodic boundary conditions. The orthogonality is to be understood with respect to the discrete scalar product defined by (3.13). Nicolaides [31] also proved some kind of similar result but did not consider periodic boundary conditions, weighted spaces and collocated meshes (he did the proof for a staggered mesh). The proof presented here does not use the same techniques as Nicolaides'.
. Moreover this decomposition is orthogonal for the scalar product (3.13).
Proof. Let us first prove orthogonality. The orthogonality between¯and −¯is obvious. Now we prove the orthogonality for the decomposition of u. We have, for any ( , ) ∈ R 2 and periodic sequence ( , ) , ∈ R in the sense of (3.11) because of (3.11). Moreover, for any , ∈ R and , ∈ R periodic in the sense of (3.11), because of the periodicity of ( , ) and ( ). Then, orthogonality of the decomposition is proved.
We shall now prove existence and uniqueness of decomposition (3.14). For there is no problem, so we only consider the equation in u. We have to prove that the function defined by }︁ . Firstly, we prove injectivity. As is a linear function, we just have to prove that Assume that for all ( , ) ∈ 1, × 1, , By the orthogonality property proved above, this implies Then, for all ∈ 1, , ( , ) 2 is an arithmetic sequence of step −2 ∆ . By periodicity, we deduce that = 0. We obtain that , +1 = , −1 for all ( , ) ∈ 1, × 1, . Then, because is odd, this implies that Note that if were not odd, there would be an even/odd decoupling here (there would exist constants odd and even such that ,2 = even and ,2 +1 = odd ).
Both equalities on , can happen simultaneously only if the values do not depend on and , and thus , is constant. Since ∑︀ , ∆ ∆ , = 0 we obtain Similarly, we obtain for all ( , ) ∈ 1, × 1, , , = 0 and the function is injective. Moreover, injectivity and the following space dimension equality ensure bijectivity: A first explanation of the wrong behavior of the Godunov scheme on a Cartesian mesh. In this section, we show that the kernel of the standard ( = 1) Godunov scheme is not an accurate approximation of the kernel of the continuous wave equation. On the other hand, we show that the kernel of the modified ( = 0) Godunov scheme does approximate correctly the kernel of the continuous wave equation.
On Cartesian meshes, this proves that deleting the diffusion term on the velocity field ( = 0) allows to recover a kernel that is an accurate approximation of its continuous counterpart.

Right or wrong behavior of the Godunov scheme in the linear discrete case
We now study the low Mach accuracy of the Godunov scheme in the sense that the numerical solution (3.2) satisfies a discrete version of (2.16). As explained in Section 3.2, the two key points to prove this kind of property is that the kernel of the scheme satisfies ℰ ℎ ⊂ Ker L ℎ and that system (3.2) is well-posed (see Thm. 3.3). The study of the kernel was performed in Section 3, > 0 and = 0 define discrete operators L , whose kernels are very different. We now prove the well-posed property ( 2 -stability) and then study the low Mach accuracy of the Godunov scheme on Cartesian rectangular and on triangular meshes.

2 -stability of the Godunov scheme
We now prove stability of the semi-discrete scheme both when = 0 and when > 0. This property is essential in the sequel. Let us define the energy Theorem 4.1. Let ( , u ) be the solution of the semi-discrete scheme (3.1). We have Then, for ≥ 0 the Godunov scheme is dissipative since Proof. We multiply the first equation of (3.1) with 2|Ω | and sum with respect to . Since does not depend on time, we obtain d d Taking the scalar product of the second equation of (3.1) with 2|Ω |u and summing with respect to , we obtain d d

The triangular mesh case
In Section 2.3.4, we saw that in the continuous setting and for low values of the Mach number, the solution of the continuous system remains close to an incompressible state for any > 0, if this was the case at the initial time = 0. This section will show that this is also the case at the discrete level for the numerical solution of the Godunov scheme applied to the linear wave equation with porosity on triangular meshes. This explains the satisfactory behavior of this scheme on this type of meshes.

The Cartesian mesh case
We saw that the Godunov scheme ( = 1) on Cartesian meshes does not preserve an incompressible state 0 ℎ ∈ ℰ ℎ, , but it preserves it if we delete the numerical diffusion on the velocity by setting = 0. From Lemma 3.8, we can define an orthogonal projection P ℎ, : We want to study the evolution over time of the initial condition when it consists in the sum of an element in the discrete incompressible space (︀ ℰ ℎ, )︀ and of a perturbation of order in (︀ ℰ ℎ, )︀ ⊥ . This will give an explanation of the wrong behavior of the standard ( = 1) Godunov scheme on a Cartesian mesh and of the satisfactory behavior of the modified ( = 0) scheme. Moreover, since completely deleting the numerical diffusion by setting = 0 was shown in [16] to present stability issues in the non-linear case, we shall also study the intermediate case = .

Explanation of the wrong behavior of the Godunov scheme on a Cartesian mesh
The next theorem shows that for the standard Godunov scheme ( = 1) on Cartesian meshes, starting from a perturbation of an incompressible field, the numerical solution will substantially deviate from the initial condition after a short time that scales like ( ), when the space discretization parameters (∆ , ∆ ) are larger than the Mach number. 2) with initial condition 0 ℎ on a Cartesian mesh with discretization parameters (∆ , ∆ ). Then, when = 1, there exists 2 > 0, depending only on , ⋆ and on T such that for almost all 0 ℎ ∈ 2 (T) 3 and for all 1 > 0, there exists 3 depending only on ( 1 , 0 ℎ ) such that for any Proof. By linearity of L ℎ, =1, , the solution ℎ of the Godunov scheme (3.2) with initial condition 0 ℎ can be written as  We have because scheme (3.2) is dissipative when ≥ 0 (see Theorem 4.1). If ⃦ ⃦ 0 ℎ − P ℎ, 0 ℎ ⃦ ⃦ 2 = 1 , then (4.7) shows that we need to find a lower bound for the function where ℎ,2 is the solution of (4.6). Before proceeding to the detailed proof of this proposition, let us briefly mention the ideas behind it: the initial condition of (4.6) will be diffused by the operator onto its orthogonal projection in the kernel Ker(L ℎ, =1, ) (this orthogonal projection is denoted by P ℎ, =1, in the sequel) and we shall prove that the solution of (4.6) will tend to P ℎ, =1, (P ℎ, 0 ℎ ) exponentially fast with a convergence rate that depends on min(Δ ,Δ ) . As a consequence, after a time that scales like ( ), the solution of (4.6) will be close enough to its projection, and thus far enough from the initial condition. To prove this in detail, we shall follow the lines below: (2) we verify thatˆℎ := ℎ,2 − P ℎ, =1, P ℎ, 0 ℎ is solution of (3.2) and thatˆℎ( ) ∈ Ker P ℎ, =1, , for all ≥ 0, (3) we use an energy estimate for solutions of (3.2) and a discrete Poincaré-Wirtinger inequality forˆℎ that is satisfied on Ker P ℎ, =1, , to estimate how fastˆℎ tends to 0, (4) we obtain (4.4) by considering times of order .
In order to obtain these results, we first prove a series of Lemmas. We start by some considerations on the orthogonal projection onto Ker L ℎ, =1, : Lemma 4.4. The operator is the orthogonal projection P ℎ, =1, on Ker L ℎ, =1, . Moreover, if ℎ is a solution of (3.2) on T with initial condition 0 ℎ , then: Proof. Recalling that Ker L ℎ,

=1,
is characterized by (3.16), it is first straightforward to see that P( ℎ ) ∈ Ker L ℎ, =1, . Next, lengthy but easy algebra leads to ⟨( ℎ − P ℎ ), ℎ ⟩ ,ℎ = 0 for all ℎ ∈ Ker L ℎ, =1, . These two properties exactly prove that P = P ℎ, =1, . Moreover, when 0 ℎ ∈ Ker P ℎ,  (4.10). As far as is concerned, this is a direct consequence of the conservativity of fluxes in the first equation of (3.1). As far as , is concerned, extracting the component of (3.1) and specializing to a Cartesian mesh, we get , −1, , where we recall that − 1 2 , is the value of the porosity at the interface between cells Ω −1, and Ω , . Multiplying the equality above with (∆ ), summing over and noting that We now write a discrete Poincaré-Wirtinger inequality for a function ℎ ∈ Ker P ℎ, =1, .

·
By multiplying by ∆ and by summing over , we have The same analysis holds for such that for all ∈ 1, ,ℓ ∆ , and we finally obtain With (4.12)-(4.14), the result follows.
To prove inequality (4.4), we first prove the following lemma which shows that ℎ,2 tends exponentially fast to the projection of its initial condition on Ker L ℎ, =1, (Items 2. and 3. above): Lemma 4.6. There exists a constant (T) > 0 depending on T and such that Proof. Let us defineˆℎ = ℎ,2 −P ℎ, =1, P ℎ, 0 ℎ := (ˆ,û) ℎ . The idea is to apply the energy estimate of Theorem 4.1 toˆℎ and then the Poincaré inequality of Lemma 4.5. For this, we first remark thatˆℎ satisfies (3.2). Indeed, ℎ,2 satisfies (4.6), and P ℎ, =1, P ℎ, 0 ℎ does not depend on time and is in the kernel of L ℎ, =1, . Then,ˆℎ is solution of (3.2) and with Theorem 4.1, we have Moreover, the initial condition ofˆℎ is P ℎ, 0 ℎ −P ℎ, =1, P ℎ, 0 ℎ , which belongs to KerP ℎ, =1, . Thus, applying (4.10) of Lemma 4.4, it follows thatˆℎ( ) belongs to Ker P ℎ, =1, for all ≥ 0 and we can apply Lemma 4.5 to estimate the right-hand side of (4.16). This leads to Then Applying Grönwall's lemma, we obtain (4.15) becauseˆ0 ℎ = (1 − P ℎ, =1, ) ∘ P ℎ, 0 ℎ . Now, we are able to prove Theorem 4.3 (Item 4. above). By applying Lemma 4.6, we have for all ≥ 0 Since the right-hand side of (4.17) is a growing function of time, we can obtain a lower bound by evaluating it at any time; we set and we obtain: Using that 1 − exp (− /2) ≥ /3 for ∈ [0; 1], Eq. (4.18) implies that for min(∆ , ∆ ) ≤ 1 In the sequel, we assume that is strictly positive, which is the case for almost all functions 0 ℎ ∈ 2 (T) 3 . Let us now suppose that then we obtain from (4.7) and (4.19) that for any ≤ 3 1 min(∆ , ∆ ). Theorem 4.3 tells us that the wrong behavior of the standard Godunov scheme is due at the same time to a wrong kernel (the image of (Id − P h, =1, ) ∘ P h, is "too large") and to a fast diffusion rate, at least proportional to min(Δ ,Δ ) . There are thus two options to propose a correction to this scheme, namely restoring a correct kernel by setting = 0 or drastically diminishing the diffusion rate by setting = . If none of these solutions is used, then a possible (but expensive) solution is to choose (∆ , ∆ ) of the size of . These three possibilities are studied in the next subsections.
Proof. By linearity of L ℎ, , , the solution ℎ of (2.9) given by scheme (3.2) with initial condition 0 ℎ can be written as  We have , then (4.24) shows that we need to find an upper bound for the function where ℎ,2 is the solution of (4.23). Assume that = 0. Since P ℎ, 0 ℎ ∈ ℰ ℎ, = Ker L ℎ, =0, we have L ℎ,
Remark 4.8. It is important to stress that the constant 3 in item 2 of Theorem 4.7 depends on a concept of discrete smoothness for 0 ℎ detailed in the next subsection and that, in the worst case, it may behave proportionally to 1 min(Δ ,Δ ) .

The case of a very fine mesh
We observe that if the right-hand side of (4.27) can be bounded by max(∆ , ∆ ) with not depending on ( , , ∆ , ∆ ), then we shall also have a bound of the type (4.21) if = 1 (uncorrected Godunov scheme) and max(∆ , ∆ ) ≤ 0 . For this, we introduce the definition of discrete regularity: Definition 4.9. Let ℎ := ( ℎ , ℎ, , ℎ, ) be a family of discrete fields parameterized by (∆ , ∆ ); then we define 2 (T) 3 to be the set of families of discrete fields such that with the following definitions for the centered and staggered finite differences in the horizontal and vertical directions This concept allows us to prove that with discrete regular initial conditions, refining the mesh is also a possibility to obtain acceptable results on an ( ) time scale when using the standard Godunov scheme. Indeed, the following theorem holds: Theorem 4.10. Let ℎ ( ) be a solution of scheme (3.2) with initial condition 0 ℎ . When = 1, for all 0 ℎ such that P ℎ, 0 ℎ ∈ 2 (T) 3 , and all 0 , 1 , 2 > 0, there exists 3 Proof. For any ℎ := ( ℎ , ℎ, , ℎ, ) , , a direct calculation shows that (4.30) Therefore, if P ℎ, 0 ℎ ∈ 2 (T) 3 , then (4.27) and the fact that ⋆ ℎ ( = 0) = 0 show that

Numerical results on the wave equation
In this section, we perform some numerical simulations on the linear wave equation with porosity (2.9) using the Godunov scheme (3.1). The aim is to illustrate all the theoretical results of the article. A 2D periodic domain T = [0, 1[×[0, 1[ is considered. All simulations were run with an Euler explicit time stepping with a CFL number of 0.4. The parameters ⋆ and are set to ⋆ = 1 and = 10 −4 . We consider a regular Cartesian mesh containing 1 600 cells (∆ = ∆ = 0.025) and an unstructured triangular mesh containing 2 326 cells generated by GMSH [19].

A stationary case
We firstly illustrate the influence of the mesh type (triangular or Cartesian) on the kernel of the Godunov scheme. The initial condition 0 = ( 0 , u 0 ) is chosen such that 0 ∈ ℰ . We take This expression of corresponds to the "vortex in a box" test case of [9]. We note that is very important from a numerical point of view because it allows to define 0 ℎ such that 0 ℎ ∈ ℰ ℎ, on Cartesian meshes and such that 0 ℎ ∈ ℰ ℎ,△ on triangular meshes. Since 0 ∈ ℰ , the field defined by is solution of the linear wave equation with porosity (2.9). We study if (5.4) is or is not satisfied at the discrete level when we solve system (2.9) with Godunov's scheme (3.1) on a Cartesian or a triangular mesh with = 0 or = 1.
In Figure 1, we plot the norm of u obtained after 1000 iterations on Cartesian and triangular meshes with = 1 and = 0. The solution is preserved over time on triangular meshes with = 1 and = 0 but is also preserved over time on Cartesian meshes with = 0. This result illustrates Propositions 3.7 and 3.9.

A well-prepared initial condition
We now consider a well-prepared initial condition. It means that the initial condition can be split into two components, a component in the kernel ℰ plus a component of order in the orthogonal set to the kernel, ℰ ⊥ . We illustrate the theoretical results Theorems 4.3 and 4.7 on the evolution with respect to time of the deviation ⃦ ⃦ ℎ − P ℎ,△ or 0 ℎ ⃦ ⃦ 2 with the different schemes on triangular and Cartesian meshes. The initial condition 0 ℎ is given by where 0 ℎ,2 ∈ ℰ ℎ, or △ is given by (5.1) and 0 ℎ,1 ∈ (︀ ℰ ℎ, or △ )︀ ⊥ satisfying ‖ 0 ℎ,1 ‖ 2 = 1. More precisely, we take  where ℎ ( , ) = (sin(2 ) cos(2 )) ℎ .
The discrete field ℎ,1 is defined at the cell centers, and so is ℎ on Cartesian meshes; but on a triangular mesh ℎ ∈ ℎ , then ℎ is defined at the edge midpoints.
In Figure 2, we plot the evolution with respect to time of the deviation ⃦ ⃦ 0 ℎ − P ℎ,△ or 0 ℎ ⃦ ⃦ 2 with scheme 6. The non linear case

Numerical schemes
Since is regular and does not depend on time, we can write system (1.1) as where W = ( , , u) and the flux f and the source term (W) are given by The numerical scheme for system (6.1) is given by where W = ( , , u) and − is the non conservative numerical flux. In this paper, we use two different fluxes, a VFRoe flux and the well-balanced Lax-Friedrich scheme of [26].

Well-balanced Lax-Friedrich scheme of [26]
The well-balanced Lax-Friedrich scheme of [26] allows to maintain equilibrium states and is easy to implement. The non-conservative numerical flux is given by where LF corresponds to the standard Lax-Friedrich numerical flux . For the existence and uniqueness of − , we refer to [26].

VFRoe scheme
We want to write a non-linear scheme that is consistent with the study we did in the linear case (see Sects. 3 and 4). We recall that in the linear case, Godunov's scheme (3.1) can be interpreted as a VFRoe scheme [6,17] in variables ( , , u) (see Rem. 3.2). The VFRoe solver consists in a local linearization of a Riemann problem which is simpler to handle since it only deals with linear problems and avoids the complex exact resolution of the Riemann problem with porosity jump. Then, in the non-linear case, we write a VFRoe scheme in variables Y = ( , , u) . Another advantage of this set of variables is that we get a scheme that is well-balanced in the sense that it exactly preserves the one-dimensional steady states. For the VFRoe scheme (with another set of variables) applied to system (6.1), we refer to [39,40]. The VFRoe numerical flux is given by where (0 − , Y , Y , n) corresponds to the solution in / = 0 − of the linearized Riemann problem that is detailed in Appendix C. The VFRoe solver considered does not allow to treat the resonant cases when eigenvalues 1 = u · n − or

All-Mach VFRoe scheme
The flux in the all-Mach VFRoe scheme is given by whereˆandˆcorrespond to VFRoe average states (see (C.1)) and = min(1, max( , )) = min(1, max(‖u ‖/ , ‖u ‖/ )). We remark that in (6.5), we recover the classical VFRoe scheme if = 1. This means that we correct the numerical flux only if both states W and W are subsonic.

Numerical results
We perform a one dimensional test to check the robustness of the low Mach corrected scheme but also the capability of the scheme to maintain equilibrium states across a discontinuous cross-section. Indeed, since the low Mach correction reduces the numerical diffusion of the scheme, stability of this scheme for unsteady low Mach flow has to be tested. Moreover, it is well-known that schemes which do not maintain the equilibrium states may give unsatisfactory results when refining the mesh [26], so that the well-balanced property also has to be tested. Then, we perform a two dimensional test to check the low Mach accuracy of the different schemes on triangular and Cartesian meshes.
For all simulations, we use the following pressure law ( ) = where = 1 and = 1.5 and CFL = 0.4.

A one dimensional unsteady subsonic flow
Let us denote U = ( , , ). The initial condition is a Riemann problem where the left state U and the right state U are given by The domain is [0, 1] and the discontinuity in the initial condition is set to = 0.5. The exact solution is 1-rarefaction followed by a stationary contact, then followed by a 3-shock. For an exact solution, we refer to [27]. The Mach number of the solution varies from 4 × 10 −4 to 0.85 and then allows to test the robustness of the all-Mach VFRoe scheme. Moreover, since is discontinuous between U and U , we also test the capability of the scheme to preserve the two invariants of the stationary contact and 2 /2 + ℎ( ) where ℎ( ) = −1 /( − 1). In Figure 3, we plot the porosity (or cross-section) , the density , the velocity , the Mach number, and 2 /2 + ℎ( ) at time = 0.25 obtained with the well-balanced Lax-Friedrich scheme [26], the VFRoe scheme and the all-Mach VFRoe scheme. The all-Mach VFRoe scheme is stable. In fact, as for the constant porosity case, numerical tests seem to show that the all-Mach scheme is stable under a degenerated CFL condition which is exactly the half of the classical one (see [5,13] for more details). This justifies why all numerical results are obtained with CFL = 0.4. As expected, the all-Mach VFRoe scheme is the least diffusive scheme and the wellbalanced Lax-Friedrich scheme is the most diffusive one. Looking at the stationary contact in = 0.5, we remark that the two invariants of the stationary contact and 2 /2 + ℎ( ) are preserved across the discontinuity of . Then, the VFRoe and all-Mach VFRoe schemes are also well-balanced, like the well-balanced Lax-Friedrich scheme.

Two-dimensional low Mach flow
We consider a two-dimensional low Mach vortex flow. Domain, meshes and boundary conditions are the same as for the wave equation (see Sect. 5). The initial condition is an exact, steady and regular solution of the incompressible system (2.7). Then, (2.6) tells us that the solution of (2.2) will remain close to the initial condition since the latter solves (2.7) for all times. Note that in order to build an exact solution of the incompressible system (2.7), we adapted the isentropic vortex solution of [43,44] to the case of variable porosity fields. The initial condition is given by where Ω = exp We firstly study from a numerical point of view if the background (order 0 in the asymptotic expansion) steady incompressible solution is preserved over time and secondly if the different schemes are accurate at low Mach number in the sense that the amplitude of the perturbation with respect to the background incompressible solution satisfies (2.6) at the discrete level. In Figures 4 and 5, we plot the norm of u obtained at time = 2 s with = 10 −4 on Cartesian and triangular meshes with the well-balanced Lax-Friedrich, the VFRoe and the all-Mach VFRoe schemes. The incompressible steady velocity seems to be preserved over time with the all-Mach VFRoe scheme on triangular and Cartesian meshes and with the VFRoe scheme on triangular meshes. With the other schemes, the solution is extremely diffused. Note that the accuracy problem of the Lax-Friedrich scheme at low Mach number on triangular mesh was already illustrated in [36] for the uniform porosity case.
In Figure 6, we study the low Mach accuracy of the different numerical schemes in the sense that we check whether (2.6) is or is not satisfied at the discrete level. For that, we study the amplitude of the deviation of the numerical solution from the incompressible solution (which is the initial condition) with respect to the Mach number. We plot the norm of the deviation for the dimensionless density˜and the dimensionless field Figure 6. Norm of the deviation from the incompressible solution for the dimensionless density and dimensionless field˜ũ at final time = 2s for Mach numbers ranging from 10 −1 to 10 −7 with the well-balanced Lax-Friedrich scheme (6.3), the VFRoe scheme (6.4) and the all-Mach VFRoe scheme (6.5) on Cartesian and triangular meshes.ũ for Mach numbers ranging from 10 −1 to 10 −7 . Recall that˜,˜andũ are defined by (2.1) where here 0 = 0 = 1 and 0 = ( 0 ) × . Note that it is very important to initialize the field u and compute the incompressible solution by using (3.7) for triangular meshes and (3.12) for Cartesian meshes, in which the discrete values of are interpolated from the analytical expression of , because otherwise (i.e. if the discrete values of u are initialized directly from their analytical expression) an error of the order of the space step will be introduced and will hide the deviation that scales like the Mach number. We observe that the VFRoe scheme is accurate at low Mach number on triangles while the all-Mach VFRoe scheme is the only one which is accurate on Cartesian and triangular meshes. Indeed, for these schemes, the density deviation is of order 2 and the velocity deviation is of order , as expected. The well-balanced Lax-Friedrich scheme and the VFRoe scheme on Cartesian meshes are not accurate at low Mach number because their velocity deviation is of order 0 , and, moreover, the density deviation of the VFRoe scheme scales like .

Conclusion
In this article, we proposed a well-balanced compressible scheme accurate at low Mach number for the Euler equations with porosity. The proposed scheme is based on the study that is performed on the linear wave equation with porosity. Indeed, the low Mach accuracy problem of the Godunov scheme can be understood and cured in the linear case. For this, we extended the discrete Hodge decomposition of [14] to a weighted 2 space in order to take into account the porosity, and we extended to the discrete level the properties that were proven by studying the modified equation related to the Godunov scheme in [15]. We enlightened the influence of the cell geometry on the accuracy of this scheme. In the triangular case, the stationary space of the Godunov scheme approaches well enough the continuous space of constant pressures and divergence-free velocity fields (up to the porosity factor), while this is not the case in the Cartesian case. On Cartesian meshes, we have to delete the usual numerical diffusion on the velocity field to preserve constant pressure fields and divergence-free velocity fields (up to the porosity factor). Moreover, as the aim was to design an all Mach regime scheme, the correction that is introduced varies continuously with respect to the Mach number. As a result, on Cartesian meshes, we propose to multiply the numerical diffusion on the velocity field by the Mach number when is smaller then 1. We check with numerical tests that this corrected scheme is accurate at low Mach number. Note that these conclusions are only valid when the boundary conditions are periodic: non-periodic boundary conditions may require additional analysis that was not performed in the present work.
The proposed non-linear scheme is based on a VFRoe solver and is a non linear extension of the Godunov scheme proposed for the linear case. The VFRoe solver avoids the exact resolution of a Rieamnn problem with variable porosity and is easy to implement. Like in the linear case, the VFRoe scheme for the non linear system is not accurate at low Mach number on Cartesian meshes but is accurate at low Mach number on triangular meshes. Based on the linear study, a fix is proposed for Cartesian meshes. This fix is easy to implement, requires only the modification of a few lines of code and allows to recover the accuracy at low Mach number on Cartesian meshes.
Further research could be driven by the following issues: First, if the porosity is discontinuous, then care must be taken in the interpretation and the numerical treatment of System (1.1). Such questions are dealt with for example in [27] and we note that, in the particular case of Section 6.5.1, the scheme proposed in the present work computes a relevant numerical approximation. A second topic that needs to be studied is the extension of the approach presented here to the full Euler system with energy balance. which allows to write (3.4). If = 0, we can only deduce (A.2) from (A.1). Nevertheless, by injecting = in the first relation of (3.3), we obtain ∀ ∈ 1, , ∀ ∈ {neighboring cell of }, ∑︁ Γ ⊂ Ω |Γ | (︀ ( u) + ( u) )︀ · n = 0, which allows to write (3.5). Let us prove that Ker L ℎ >0, Ker L ℎ =0, . Let ℎ ∈ Ker L ℎ >0, . We have for all ∈ 1, and then ℎ ∈ Ker L ℎ =0, .
The solution (0 − , Y , Y , n) is given by