A CLOSER LOOK AT LOW-SPEED PRECONDITIONING TECHNIQUES FOR THE EULER EQUATIONS OF GAS DYNAMICS

. It is widely known that compressible ﬂow solvers usually behave badly in terms of both accuracy and convergence rate at the incompressible limit. The present paper focuses on presenting the fundamentals of low-speed preconditioning techniques in a detailed approach, with particular emphasis on the Euler equations in two-dimensional generalized coordinates. Issues related to eigenvalue scaling, local time-stepping correction, artiﬁcial dissipation adjustment and boundary conditions implementation are discussed. The computational tests are conducted using a ﬁnite difference code, which is constructed for structured grids in general curvilinear coordinates. Spatial discretization uses central differences plus added artiﬁcial dissipation. The time march can be performed with either explicit or implicit time discretization methods. Although the investigation is performed with a speciﬁc code, the conclusions reached in the process of the present research are quite general and, therefore, useful for CFD practitioners regardless of the type of method they work with. Several airfoil simulations are addressed in order to demonstrate the beneﬁts or even the necessity of using low-speed preconditioning techniques in the context of aerospace applications.


INTRODUCTION
Historically, once academia began to recognize difficulties in simulating low-speed flows with codes based on general compressible formulations, there was an active effort in devising solvers based directly on the incompressible flow equations.In return of such effort, successful techniques were developed for the simulation of incompressible flows, as the well known pressure-based schemes [6] or the pseudo-compressibility methods [3].On the other hand, low-speed preconditioning techniques were created not only because it is more practical to have a single code capable of dealing with both compressible and incompressible regimes, instead of separately maintaining different dedicated codes, but mainly because there is a variety of situations in which only a solver capable of simultaneously handling well both regimes can give satisfactory results, such as in Ref. [1]: • high-speed flows with large embedded regions of low velocity, e.g., subsonic flows upstream of a strongly converging nozzle; • low-speed flows that are compressible due to density changes induced by heat sources, e.g., surface heat transfer and combustion simulations; • flows where compressible and incompressible regimes occur side by side and interact freely, e.g., propulsion and high-lift configurations.
The simplest way to understand the need for preconditioning is perhaps by noticing the disparity between convection and acoustic speeds at the incompressible limit.Considering the one-dimensional Euler equations, one would have the set of eigenvalues {u, u+a, u−a}, where u and a are the particle and sound speeds, respectively.The so called condition number (C N ), defined as the ratio between the highest eigenvalue (in absolute value) and the lowest eigenvalue (in absolute value), may be understood as a measure of time step restriction, or stiffness, when considering the eigenvalues of the Jacobian matrix of a hyperbolic system of equations to be numerically simulated [1].For the subsonic regime, where u − a < 0 < u + a, the condition number can be written as a function of the local Mach number (M ) by Once lim M →0 C N = +∞, common flow solvers present serious time step restrictions at low-speed regions and, consequently, overall convergence problems.The approach in lowspeed preconditioning techniques consists, essentially, in altering the original partial differential equations in order to have more appropriate eigenvalues, therefore avoiding, or at least delaying, the condition number divergence at the incompressible limit.This way, codes based on general compressible formulations become capable of simulating nearly incompressible flows with relative easiness.
At the moment, there is a pacified understanding among researchers that preconditioning methods consist of successful and well established techniques for dealing with simulations involving low-speed flows.Despite of that, one can have great difficulty in finding an introductory-level contribution that discuss the physical essence common to all of the preconditioning approaches as well as the basic implementation issues in a detailed manner.The present study comes in to address such difficulty and basically intends to popularize and make lowspeed preconditioning techniques better understood among general CFD practitioners.To that end, the two-dimensional Euler equations of gas dynamics are solved using the classical finitedifference discretization in generalized curvilinear coordinates, although the discussed concepts and principles have much wider applicability.This paper is organized as follows.Section 2 reviews the two-dimensional Euler equations in Cartesian and curvilinear coordinates, introduces the preconditioned formulation and discusses how the eigenvalues are modified by the preconditioning.Section 3 addresses important details such as how to alter the CFL-based local time step in order to obtain maximum efficiency from the preconditioning technique, as well as issues related to artificial dissipation adjustment and boundary conditions implementation.In section 4, several airfoil simulations are analyzed in order to demonstrate the benefits or even the necessity of using low-speed preconditioning techniques in the context of aerospace applications.

The General Compressible Formulation
In conservative form, the two-dimensional Euler equations in Cartesian coordinates are given by ∂Q ∂t where Q = [ρ, ρu, ρv, E] T is the vector of conserved variables and F x = [ρu, ρu 2 +p, ρuv, (E + p)u] T and F y = [ρv, ρvu, ρv 2 + p, (E + p)v] T are the inviscid flux vectors.Moreover, ρ stands for density, u and v are the Cartesian components of the velocity vector, is the total energy per unit volume, and e i is the specific internal energy.The static pressure p is obtained using the equation of state for a perfect gas, namely, p = ρ(γ − 1)e i , where γ is the fluid ratio of specific heats, which assumes the value γ = 7/5 for the air.
For codes based on structured grids, the Euler equations are re-written in generalized curvilinear coordinates as where Q = J −1 Q, Fξ = J −1 (ξ x F x + ξ y F y ), Fη = J −1 (η x F x + η y F y ), and J −1 = x ξ y η − x η y ξ is the reciprocal value of the transformation Jacobian.The metric terms composing the flux vectors can be calculated through the identities [4] ξ x ξ y η x η y = ∂ξ/∂x ∂ξ/∂y ∂η/∂x ∂η/∂y As will be seen in the following sections, it is convenient, for preconditioning purposes, to work with the Euler equations written in terms of the vector of primitive variables Naturally, Z = J −1 Z, A = (∂F x /∂Z)ξ x + (∂F y /∂Z)ξ y and B = (∂F x /∂Z)η x + (∂F y /∂Z)η y .Furthermore, the Jacobian matrices ∂Q/∂Z, ∂F x /∂Z and ∂F y /∂Z are given explicitly by Therefore, the Euler equations in general curvilinear coordinates can be rewritten as where it should be noted that absolutely no preconditioning was applied up to this point.In order to obtain the eigenvalues associated with Eq. ( 11), it should be pre-multiplied by (∂Q/∂Z) −1 , which gives where whose sets of eigenvalues are the usual (curvilinear) ones, namely, where a = (γp/ρ) 1/2 is the speed of sound and U = uξ x + vξ y and V = uη x + vη y are the well-known contravariant velocity components [4].

Introducing the Preconditioning Matrix
Different preconditioning methods are obtained depending on how the eigenvalues are modified.The present work follows the technique introduced in Ref. [9], which was found to perform very well for a diversity of flow cases.The main idea is to multiply the first column of the matrix ∂Q/∂Z by a 1/ factor, which can be incorporated into the Γ matrix, given by Thus, Eq. ( 11), after preconditioning, becomes In order to obtain the eigenvalues associated with Eq. ( 18), it should be premultiplied by Γ −1 , which gives The new matrices are defined as whose sets of preconditioned eigenvalues are The preconditioned eigenvalues return the original ones as → 1.
It is now important to consider the effect of preconditioning in the steady-state solution of the modified equation.For this sake, Eq. ( 18) can be written as where R stands for the residue of the Euler equations which, recalling Eqs. ( 6) and (7), is given by This term is, naturally, equal to zero when the steady-state solution is reached by means of the preconditioned formulation, as indicated in Eq. ( 24).Therefore, in both formulations, i.e., original or preconditioned, the same result is obtained at the steady state.
Finally, Eq. ( 24) can also be solved for the vector of the conserved variables through Eq. ( 5), which gives where M pc = (∂Q/∂Z) Γ −1 is called the preconditioning matrix in the present study.This matrix is given explicitly by It is worth noting that Eq. ( 26) makes very clear that, even though the steady-state solution of the preconditioned equation is precisely the same of the original one, the intermediate states are not, since the preconditioning matrix does alter the coherence of the temporal evolution.Hence, preconditioning techniques alone are of no use when there is interest in transient states or when simulating unsteady flows.Fortunately, preconditioning techniques are easily extendable to unsteady flow applications in conjunction with dual time stepping algorithms.

Analyzing the Eigenvalue Scaling
Up to this point, no definition was given for the factor .In order to better visualize how this parameter affect the eigenvalues, it is useful to consider once again (as in section 1) the one-dimensional Euler equations in Cartesian components.The eigenvalues obtained by preconditioning such equations would be The main strategy is to define as a function of the Mach number in order to maintain the condition number "well behaved" in low-speed flow regions.This work follows the simple suggestion of Ref. [8], where being M the local Mach number (since the preconditioning matrix has to be evaluated at each grid node when performing a numerical simulation) and M ∞ is a reference free-stream Mach number value (for external airfoil flows of interest here, M ∞ is the prescribed Mach number at infinity).
Figures 1 to 4 illustrate how the set of eigenvalues in Eq. (28) behaves as function of the local Mach number M for different values of the free-stream Mach number M ∞ when approaching the incompressible limit.In such figures, each color is associated with one eigenvalue (u + a in blue, u in green and u − a in red).The condition number C N is given in black.Also, curves shown with "empty dots" refers to the original eigenvalues and curves with "full dots" refers to the preconditioned ones.As can be seen in figures 1 to 4, the original and preconditioned eigenvalues coincides if M > 1.For M = 1, a well-known C N divergence is observed since the lowest eigenvalue (in absolute value) approaches zero.Such divergence is not critical for most flow cases and is not addressed here.The reader is referred to Ref. [2] and references therein for further discussions.
For subsonic cases, figures 1 to 4 show effective difference between the original and preconditioned eigenvalues, but such difference becomes effective for the condition number only when M < M ∞ , precisely where a discontinuous change of inclination can be noted in the curves associated with the preconditioned eigenvalues.Clearly, the preconditioned condition number shows better behavior in terms of divergence when compared with the original condition number, effect caused by the clustering of eigenvalues around the "central" one, which is not modified by the preconditioning and goes to zero at the incompressible limit, resulting in a bounded ratio of eigenvalues.
It is also important to point out that the divergence of the preconditioned condition number is delayed appropriately with decreasing M ∞ , assuring that the stiffness problem is properly bypassed, as expected of a good low-speed preconditioning method.

Adjusting the CFL-Based Time Step Restriction
Numerical experiments show that if the same constant, or fixed, time step is used to simulate a given flow case, the preconditioned formulation will be less efficient in terms of convergence rate than the original formulation, and this difficulty becomes worse with decreasing Mach number.Such assertion sounds contradictory since the preconditioning main purpose is obviously the gain in terms of convergence rate, specifically when dealing with low-speed flows.
In order to understand such problem, it is useful to consider the classical shock-tube problem [7].In such problem, the temporal evolution happens through the traffic of information carried by shock and expansion waves (and also by the contact discontinuity) which are associated with the standard set of eigenvalues {u, u + a, u − a}.Essentially, the information speed is directly related with the magnitude of the problem's eigenvalues.The fact is that such concepts can be extended naturally for general numerical flow simulations since, beginning from a given initial condition, information waves transmit modifications on the physical proprieties during convergence until a steady-state is achieved.
As seen if figures 1 to 4, the magnitude of the preconditioned eigenvalues is very little (especially at the incompressible limit) when compared with the original ones, and that is the reason why, when using a same constant/fixed time step, one gets slower convergence for the preconditioned formulation.Such disadvantage, however, comes from the same root that makes low-speed preconditioning techniques so efficient for low-speed flows, namely, the clustering of eigenvalues.The explanation is that the preconditioned system of equations admits a much larger time step (at the limit of numerical divergence) than the original system of equations, since the original eigenvalues easily extrapolates the boundaries of the stability envelope (having in mind constant/fixed time-stepping algorithms) because of their magnitude disparity.Therefore, the effectiveness of preconditioning consists essentially in allowing larger numerical time steps.The simplest way to explore such capacity is by using a CFL-based preconditioned local time step.
The alteration in the CFL-based time step is straightforward and basically consists in using the preconditioned eigenvalues instead of the original ones.For the one-dimensional Euler equations with no preconditioning, a typical expression for the time step based on the CFL number would be ∆t For the preconditioned equations, the eigenvalues in (28) would be used instead, giving Now, for the bi-dimensional formulation (in curvilinear coordinates) with no preconditioning, a typical CFL-based time step would be where r ξ and r η are obtained from ( 15) and ( 16), being given by Analogously, for the preconditioned equations, the preconditioned eigenvalues should be used, giving where r pc ξ and r pc η are now obtained from ( 22) and ( 23), being given by

Dealing with Numerical Dissipation and Boundary Conditions
In this work, where artificial dissipation is added explicitly to the formulation, it was found convenient to leave the dissipation term unaltered (in its original form) and to accommodate it with the residue R given in (25), being the stabilized version of equation ( 26) given by d Q where R stab = R + D would be the stabilized residue, being D the original dissipation term.By doing this, the steady state solution obtained is precisely the same given by the original formulation, namely, the one corresponding to R stab = 0.Such implementation is obviously independent of the specific numerical viscosity model being used and is, therefore, quite general.
It has been shown that numerical dissipation generally scales badly as the Mach number is reduced [5], but it was found in the present study that such problem becomes significant only for extremely low free-stream Mach numbers, being of no practical importance for a wide variety of low-speed test flow cases simulated.Additionally, it is also known that better performance can in fact be extracted of the preconditioned formulation if the dissipation term is modified according to the preconditioned eigenvalues.In order to avoid the overhead associated with the evaluation of preconditioning matrices for the artificial viscosity in each grid node, the present work did not followed such modification ideas further, since the performance gain was found already very satisfactory.The reader interested in how to modify the dissipation term may consult Ref. [8].
Concerning the implementation of eigenvalue-based boundary conditions such as the usual non-reflecting far-field ones (very common in the context of external flow simulations), a precise approach should consider the modified characteristic variables (or Riemann invariants) of the preconditioned system.In the present study, however, it was found that no significant differences become apparent by using the original (unmodified) boundary treatment, at least for steady-state problems.Heuristically (yet having in mind the shock-tube problem), in order to minimize errors from such source, it is important that the symmetric form of the eigenvalues {u − f (a), u, u + f (a)} is maintained by the preconditioning method being used.The reader interested in how to modify the far-field boundary conditions may also consult Ref. [8].

NUMERICAL RESULTS
In this section, some airfoil low-speed simulations are addressed in order to investigate the gain in convergence rate and accuracy of the preconditioned formulation over the original one.All simulations were performed using the classical 2nd order-accurate Runge-Kutta (explicit) time-stepping scheme.Also, in all the pairs of compared convergence curves, the same CFL-based time step was used, as explicitly informed in each figure.The well-known NACA0012 and RAE2822 airfoils were chosen as examples of typical profiles used in the aerospace context.Their related computational meshes can be seen in figures 7 and 10, respectively.It should be said that both meshes are of circular type and have 81x40 (NACA) and 129x200 (RAE) grid points along surface and radial coordinates, in that order.For each case, the far-field was "safely" fixed at a considerable distance (typically tenths of chords) away from the airfoils.
Figure 5 shows how convergence speed is loosed with decreasing free-stream Mach number.This effect is true for both preconditioned and original formulations, but it is far more critical without preconditioning.All charts of figure 5 refer to simulations for the NACA0012 profile with zero angle of atack (AOA).It is important to note that better convergence speed is obtained with the preconditioned formulation when the such rate is based on iterations.If based on computational time, the convergence speed is not always greater because the evaluation of the preconditioning matrix for all grid nodes in each iteration is somewhat costly.
The gain in convergence speed, however, is not the only benefit of using preconditioning techniques.In fact, figure 6 illustrates how steady-state converged flow fields can be wrongly predicted by using the original formulation at low-speeds.The distributions of pressure coefficient (Cp) are compared to the "exact" ones, which were obtained with XFOIL R .For each airfoil referred in figure 6, the error in the Cp distribution may lead to significant errors in the force and (specially) moment coefficients, since greater differences appears generally at trailing and leading edges.
To better compare the flow-field differences for the cases referred in figure 6, Mach number contours obtained with and without preconditioning are put side-by-side in figures 8 and 9 for the NACA0012, and in figures 11 and 12 for the RAE2822.In these four figures, the colormaps are adjusted for best contrast "individually", but all the isolines displayed corresponds to the exact same values in order to highlight the clear flow-field differences between steady-state converged solutions with and without preconditioning.

CONCLUDING REMARKS
This work presented a simple way to implement a low-speed preconditioning method for the Euler equations of gas dynamics with particular emphasis on the finite difference formulation based on generalized curvilinear coordinates.Nevertheless, concepts and principles common to all of the preconditioning methods were discussed in a detailed and didactically manner.Several airfoil simulations were performed in order to investigate the benefits or even the necessity of using preconditioning algorithms within the context of aerospace applications.It is hoped that the present study may contribute in turning low-speed preconditioning techniques better understood among general CFD practitioners.