On the behaviour of upwind schemes in the low Mach number limit

Abstract

This paper presents an asymptotic analysis in power of the Mach number of the flux difference splitting approximation of the compressible Euler equations in the low Mach number limit. We prove that the solutions of the discrete system contain pressure fluctuations of order of the Mach number while the continuous pressure scales with the square of the Mach number. This explains in a rigorous manner why this approximation of the compressible equations fails to compute very subsonic flow. We then show that a preconditioning of the numerical dissipation tensor allows to recover a correct scaling of the pressure. These theoretical results are totally confirmed by numerical experiments.

Introduction

The computation of low Mach number flows is continuing to be a challenge for compressible flow solvers and a general algorithm valid for all flow speeds is still to be found. While many different flow regimes can be encountered for small Mach number[1], this paper considers exclusively the simple case where the solutions of the equations for compressible flows converge to solutions satisfying the equations for incompressible flows. From a theoretical point of view, the situation is now well understood in the inviscid case: if the initial pressure field scales with the square of the Mach number: p(x, 0)=P₀+M²_∗p₂(x), and if the velocity field at time t=0 is close to a divergence free field: u(x, 0)=u₀ (x)+M_∗u₁(x) with div(u₀)=0, then it is known that the solutions of the equations for compressible flows remain uniformly bounded as the Mach number tends to zero, and that the limit solutions satisfy the equations for incompressible flows. A rigorous proof of this result can be found in Ref.[2], while a general treatment of hyperbolic problems exhibiting different time scales can also be found in works by Kreiss, Gustafsson and their co-workers (see Refs3, 4). These last works point out that with proper initialization (see for instance Ref.[5], where it is shown how to construct initial data such that the solution and their derivatives remain bounded independently of the ‘small parameter’) the solutions of the complete equations converge to solutions satisfying a set of ‘reduced’ equations that for fluid dynamics problems are the incompressible Euler equations.

However, if on the continuous case, the situation is clear, in the discrete case, the problem of finding on a fixed grid an approximation of the equations for compressible flows whose accuracy will not deteriorate as the Mach number decreases, is still a problem of great practical and theoretical interest.

Actually, it is well known that in addition to convergence and round-off error difficulties, the approximations of the compressible fluid flow equations suffer from accuracy problems in the low Mach number limit. There are experimental evidences showing that on a fixed mesh, the discretized solution of the compressible fluid flow equations are not an accurate approximation of the incompressible equations (e.g. see Refs6, 7).

A common explanation of this failure of the compressible solvers is to put forward the role of the numerical dissipation. In this paper, we are interested in upwind schemes, where the numerical dissipation associated with the two grid nodes i and l is of the form:

$$\text{δ∣}\text{A}\text{∣Δ}{}_{\mathrm{<mtext>il</mtext>}}\text{q=δ}\text{∑}\text{i}\text{∣λ}{}_{\mathrm{i}}\text{∣α}{}_{\mathrm{i}}\text{(Δ}{}_{\mathrm{<mtext>il</mtext>}}\text{)}\text{R}{}_{\mathrm{i}}\text{,}$$

with δ the mesh size, λ_i and R_i respectively the eigenvalues and eigenvectors of the Jacobian matrix A, α_i(Δ_il) the coordinate of the jump Δ_ilq between the values of the function q at the grid nodes i and l in the basis formed by the eigenvectors. In the above expression, the sum runs over all waves. If we introduce a reference Mach number M_∗ we see that this implies that some of the wave speeds λ_i become proportional to 1/M_∗. Thus, in the low Mach number limit these terms can grow without bound, resulting in an excessive numerical dissipation that pollutes the discrete solution.

We believe that this explanation is only partially true and needs to be refined. Firstly, we remark that, if the dissipative terms are of order 1/M_∗, the centred terms in the approximation scale with 1/M²_∗ (see Section 2). Thus, although the amplitude of the dissipative terms may go to infinity, they become negligible with respect to the centred ones as M_∗ goes to zero. Therefore, the accuracy problem may actually result from a lack of numerical dissipation, instead of resulting from an excess of numerical viscosity. Note that this fact has been observed before in related work carried out by Turkel, Fitermann and Van Leer[8], where for scalar numerical viscosity it was observed that a density based viscosity gives a numerical dissipation term that is too small. Secondly, in upwind schemes, the numerical dissipation is a 4×4 tensor (in 2D) and although some of the λ_i goes to infinity, it may well be that the discrete equations force the associated jumps α_i (Δ_il) to go to zero.

Our aim in this work is thus to explain in a detailed manner the mechanism by which compressible upwind schemes fail to compute very subsonic flows. The tool we will use for this study is only a standard asymptotic analysis in power of the Mach number. However, this analysis will be performed on the discrete equations instead of the continuous ones. More specifically, we will look for the asymptotic form that takes the discrete system in the low Mach number limit. The analysis we performed is a fully non-linear one that does not require any linearization of the governing equations. Moreover, it applies to the steady-state equations as well as the time-dependent ones. Its main result can be stated as follows: in the low Mach number limit, the solution of the discrete equations contains pressure fluctuations of the order of the Mach number. This is in clear contrast with the continuous case, where pressure fluctuations scale with the square of the Mach number.

For the steady-state equations, a common strategy to overcome the convergence problems encountered by compressible solvers is the use of preconditioning8, 9, 10, 11, 12. Following a suggestion in Ref.[8], we tried in Ref.[13] to use preconditioning as a cure to the accuracy problem. The results were extremely convincing and it appears that preconditioning is a powerful remedy to cure the accuracy problem. More specifically, the dissipative term in Eq. (1)δ∣A∣Δq is changed into:

$$\text{δ}\text{P}{}^{\mathrm{-1}}\text{∣}\text{PA}\text{∣Δq}$$

where P is the preconditioning matrix. However, the temporal and centred terms of the approximation remain unchanged (i.e. with no preconditioning). Therefore, the scheme is always consistent with the time-dependent equations and only the numerical dissipation is altered. Although more complex than in the case with no preconditioning, an asymptotic analysis of the resulting discrete equations is tractable. It provides a detailed explanation of the mechanism by which preconditioning increases the accuracy of the schemes. In particular, it shows that preconditioning increases the numerical dissipation terms associated to the continuity and energy equations by a factor of 1/M_∗. Since the original non-preconditioned equations contain terms of order 1/M_∗, the resulting preconditioned system contains only terms of order 1/M²_∗ and 1, as in the continuous case. Solving these equations gives pressure that possesses the same M²_∗ scaling as for the continuous system.

The analysis presented in this work uses Roe’s approximate Riemann solver. This scheme possesses some nice algebraic features that simplify, to some extent, the discrete equations. However, the analysis is general and applies to all the schemes that can be cast under the form in Eq. (1), with a matrix A that is the Jacobian of the continuous fluxes evaluated at some average between i and l. The outline of this paper is as follows: in the next section, we recall how to derive the convergence of the compressible Euler equations toward the incompressible ones. Next, we apply the same asymptotic expansions as in the continuous case to the discrete Roe scheme, and finally, in the last section, we examine the case of a preconditioned dissipation.