9.5 Preconditioning for Low Mach Numbers 低马赫数预条件[cfd-9-5]

We demonstrate preconditioning with the aid of 1-D Euler equations. They can be written in differential form as

where \(\vec{W} = [\rho,\ \rho u,\ \rho E]^T\) is the vector of conservative variables and \(\vec{F}_c\) denotes the convective fluxes. In order to see the effect of preconditioning on the spectral radii (important for the computation of the time step) and on the convective flux Jacobian (important for upwind dissipation), the Euler equations (9.41) are rewritten in the quasilinear form

with \(\bar{A}_c = \partial\vec{F}_c/\partial\vec{W}\) being the convective flux Jacobian (cf. Appendix A.2).

The idea behind low Mach-number preconditioning is to transform the governing equations (9.41) such that the new equations have more favourable properties at low Mach numbers (i.e. below \(M \approx 0.2\)). First of all, we want to equalise the convective eigenvalues in order to bound the condition number (see Eq. (9.40)) and thus to remove the stiffness at \(M \rightarrow 0\). We further want to change the scaling of the numerical dissipation in order to improve the accuracy. Finally, we want to couple the pressure and the velocity as it is done in the pressure-based schemes. The transformation replaces the conservative variables \(\vec{W}\) by a different set of flow variables. Various choices are possible [90], but the most often used are the pressure, the velocity components and the temperature. Transforming the quasilinear form in Eq. (9.42) into the new variables \(\vec{W}_p\), we obtain

In the above Eq. (9.43), \(\bar{P} = \partial\vec{W}/\partial\vec{W}_p\) represents the transformation matrix from the new variables \(\vec{W}_p\) into the conservative variables \(\vec{W}\). Introducing a new flux Jacobian

we can write Eq. (9.43) as

If we now replace \(\bar{P}\) in front of the time derivative by a suitable preconditioning matrix \(\bar{\Gamma}\) (to be specified later), the preconditioned Eq. (9.45) reads

or equivalently

It is now possible to solve Eq. (9.47) in the new variables \(\vec{W}_p\) using any standard spatial and temporal discretisation scheme. Another possibility is to transform the preconditioned equations (9.47) back into the conservative variables \(\vec{W}\). We have

where \(\bar{P}^{-1} = \partial\vec{W}_p/\partial\vec{W}\) is the transformation matrix from the conservative variables. Equation (9.48) can be for convenience written as

or

The term \(\bar{P}\bar{\Gamma}^{-1}\) in Eq. (9.49) and (9.50) is called the conservative variable preconditioning matrix.

We can see now that the preconditioned equations in the conservative variables (9.49), (9.50) have the following features:

• Spatial derivatives are multiplied by \(\bar{P}\bar{\Gamma}^{-1}\).
• Unsteady equations are different from the original form in Eq. (9.42) and hence the solution is no longer time accurate.
• Stationary solution (i.e. \(\partial\vec{W}/\partial t = 0\)) remains unchanged.
• Eigenvalues and eigenvectors of the preconditioned system correspond to those of the matrix \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\) and hence are different from those of the original system in Eq. (9.41).

The main task is now to find a suitable set of variables \(\vec{W}_p\) and the matrix \(\bar{\Gamma}\) such that the convective eigenvalues of the preconditioned system (9.50) are equalised as close as possible. But it is also important that the matrices remain defined for \(M \rightarrow 0\). Furthermore, it is desirable that the preconditioned system converts into the original system for higher Mach numbers.

Before we present the transformation and preconditioning matrices in Subsection 9.5.3, we shall discuss the implementation of the low Mach number preconditioning in a flow solver.

Since we solve the governing equations in their integral form, we have to write Eq. (9.50) as

It is important to note that Eq. (9.51) is not conservative for unsteady flows. Therefore, the dual time-stepping approach (see Section 6.3) has to be employed in order to obtain a time-accurate solution.

The solution of the preconditioned governing equations (9.51) using an explicit multistage scheme (described in Section 6.1) proceeds according to the following steps:

• 1. Compute a new time step based on the spectral radii of the matrix \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\).
• 2. Evaluate artificial dissipation using the spectral radii (e.g., JST scheme) or the eigenvalues and eigenvectors of the preconditioned flux Jacobian (e.g., Roe upwind scheme). The dissipation can be formulated either in the conservative or in the primitive variables [90].
• 3. Compute the convective fluxes (either without change or based on the new variables \(\vec{W}_p\)).
• 4. Sum up the dissipative and convective fluxes. Depending on the formulation of the artificial dissipation, either the whole residual or just the convective terms are multiplied by the conservative variable preconditioning matrix, i.e. by \(\bar{P}\bar{\Gamma}^{-1}\).
• 5. Multiply the residual by \(\alpha_k\Delta t/V\).
• 6. Carry out the implicit residual smoothing (optionally).
• 7. Subtract the residuals from the old conservative variables \(\vec{W}^{(0)}\) in oder to obtain the new conservative variables \(\vec{W}^{n+1}\).
• 8. Update the boundary conditions (note that all boundary conditions which are based on the characteristic variables - like inflow, outflow, farfield - need to be changed).

Similar procedure is followed for an implicit scheme, only the steps 5. and 6. are omitted and the conservative variables are updated in a different way (cf. Section 6.2). It should also be noted that the preconditioning has to be included in the implicit operator since the Jacobian is \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\).

The preconditioned scalar dissipation scheme (JST) takes the form (see Eq. (4.50))

Note that the term \(\bar{\Gamma}\bar{P}^{-1}\) is required if the complete residual is later multiplied by \(\bar{P}\bar{\Gamma}^{-1}\), since the spectral radius \(\hat{\Lambda}^S\) already contains the preconditioning matrix (thus it is different from Eq. (4.53)). It is also possible to combine \(\bar{P}^{-1}\) and \(\vec{W}\) into \(W_p\) in Eq. (9.52). Then, we can formulate the preconditioned scalar dissipation scheme in primitive variables as

It is important to realize that in this case only the convective fluxes are multiplied by \(\Gamma^{-1}\) and the complete residual then by \(\bar{P}\). The nonconservative form in Eq. (9.53) is somewhat simpler than the relation (9.52), however there are problems with the accuracy for internal flows or for flows containing shocks.

The preconditioned Roe upwind scheme can be written as (cf. Eq. (4.91))

where \(|\bar{A}_{Roe}| = \bar{T}_{c,p}|\bar{\Lambda}_{c,p}|\bar{T}_{c,p}^{-1}\). The left (\(\bar{T}_{c,p}^{-1}\)) and right (\(\bar{T}_{c,p}\)) eigenvectors, as well as the eigenvalues (\(\bar{\Lambda}_{c,p}\)) are those of the matrix \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\). Values of the flow variables at \((I + 1/2)\) are obtained by Roe's averaging given in Eq. (4.89), and the whole residual is multiplied by \(\bar{P}\bar{\Gamma}^{-1}\). Again, the preconditioned Roe scheme in Eq. (9.54) can be formulated in the primitive variables \(\vec{W}_p\) leading us to

The eigenvalues and eigenvectors which compose \(\bar{A}_{Roe,p}\) are now determined by the matrix \(\Gamma^{-1}\bar{A}_{c,p}\) (cf. Eq. (9.47)). The eigenvalues of \(\Gamma^{-1}\bar{A}_{c,p}\) are identical to the eigenvalues of \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\), i.e. to \(\bar{\Lambda}_{c,p}\), but the eigenvectors are different. Thus, the Roe matrix in Eq. (9.55) is composed as \(|\bar{A}_{Roe,p}| = \bar{T}_p|\bar{\Lambda}_{c,p}|\bar{T}_p^{-1}\). It should be mentioned that also in this case the residual has to be multiplied by \(\bar{P}\bar{\Gamma}^{-1}\) in order to obtain conservative variables.

Among the various choices for the primitive variables \(\vec{W}_p\), the form

appears most often. Therefore, we shall restrict the further discussion to this particular form of \(\vec{W}_p\). In the following, we will present the transformation and the preconditioning matrices together with the eigenvalues and the left and right eigenvectors for a general fluid, as well as a perfect gas.

For a general fluid, the transformation matrix from the conservative into the primitive variables \(\bar{P}^{-1} = \partial\vec{W}_p/\partial\vec{W}\) is given by [91]

with

The transformation matrix from the primitive into the conservative variables \(\bar{P} = \partial\vec{W}/\partial\vec{W}_p\) reads [91]

Derivatives of the density and of the enthalpy with respect to the pressure and the temperature in Eqs. (9.57)-(9.59) can be written in the form

where \(\alpha_p\) and \(\alpha_T\) are the compressibility coefficients at constant pressure and temperature, respectively. The speed of sound can be computed from

In the case of a perfect gas (see Subsection 2.4.1), the compressibility coefficients in Eq. (9.60) become \(\alpha_p = 1/p\) and \(\alpha_T = 1/T\). In this case, the transformation matrix from the conservative into the primitive variables \(\bar{P}^{-1} = \partial\vec{W}_p/\partial\vec{W}\) from Eq. (9.57) can be cast into

The transformation matrix from the primitive into the conservative variables \(\bar{P} = \partial\vec{W}/\partial\vec{W}_p\) reads for a perfect gas

with \(q^2\) as defined in Eq. (9.58).

The construction of a preconditioning matrix is relatively easy in the case of the Euler equations. The formulation proposed by van Leer at al. [92] achieves the lowest attainable condition number. The methodology was discussed in detail in [79]. However, the preconditioning of the Navier-Stokes equations is more involved. The reason is that the viscous terms lead to complex wave speeds, which makes the preconditioned system difficult to analyse. The most recognised preconditioners for viscous flows were proposed by Choi and Merkle [93], [94], Turkel [90], [95]-[97], Lee and van Leer [98], [99] and Lee [80], Jorgenson and Pletcher [100], and by Weiss and Smith [101], [45], respectively. Examples of applications can be found in Refs. [91], [102]-[109].

The preconditioning matrix \(\bar{\Gamma}\) due to Weiss and Smith [101], [45] has, in the case of a general fluid, a form identical to \(\bar{P}\) in Eq. (9.59) with \(\rho_p\) replaced by a suitable preconditioning parameter \(\theta\). The same holds also for the inverse of the preconditioning matrix, i.e., \(\bar{\Gamma}^{-1}\) which resembles \(\bar{P}^{-1}\) from Eq. (9.57). Consequently, the parameter \(a_1\) from Eq. (9.58) is changed into

In the case of a perfect gas, the preconditioning matrix \(\bar{\Gamma}\) can be written as

where \(\theta\) is again the preconditioning parameter, which will be defined later. The inverse of the preconditioning matrix is given by

with the abbreviations

Several choices exist for the preconditioning parameter \(\theta\) in Equations (9.64)-(9.67). One definition, which is based on a reference velocity \(u_r\) was provided in Ref. [45]. It reads for a general fluid

or, correspondingly for a perfect gas,

The reference velocity \(u_r\) in Eq. (9.68) or Eq. (9.69) is obtained from the relation

where \(\Delta h\) is a measure of the control volume size, the quantity \(\Delta p\) stands for the pressure difference between the adjacent control volumes and \(\epsilon\) is a small number (\(\approx 10^{-3}\)). As we can see from the definition in Eq. (9.70), the reference velocity is bounded by the local transport velocity. The terms \(\nu/\Delta h\) and \(\kappa/\Delta h\) become important in boundary layers with dominant diffusion or heat conduction. The pressure term is intended to prevent \(u_r\) from vanishing at stagnation points. In the case that \(u_r = c\), \(\bar{\Gamma}\) becomes identical to \(\bar{P}\) and \(\bar{\Gamma}^{-1}\) is converted into \(\bar{P}^{-1}\). Hence, the preconditioning is turned off for a supersonic flow as intended.

Another possibility, which was devised for a perfect gas, is to set [104]

The parameter \(\beta\) in Eq. (9.71) is defined as

with the reference Mach number given by

and \(M\) being the local Mach number (\(M^2 = \|\vec{v}\|_2^2/c^2\)). According to Ref. [104], parameter \(M^2_{min} = K M_\infty^2\) and \(K \approx 3\) (others choose \(K = 1\) or even \(K = 0.15\); the exact value seems to depend on the number of control volumes in the stagnation region or inside the boundary layer). As it can be easily verified, when \(M \ge 1\) the parameter \(\beta\) equals to \(1/\gamma\) and the matrix \(\bar{\Gamma}\) becomes identical to \(\bar{P}\). It should be mentioned that both definitions of \(\theta\) in Eq. (9.69) and in Eq. (9.71) are equivalent. Hence, the reference Mach number could be determined as \(M_r = u_r/c\) with the help of Eq. (9.70).

The matrix of the eigenvalues of the preconditioned system Eq. (9.49), i.e, \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\) is given by

where \(V = \vec{v}\cdot\vec{n}\) represents the contravariant velocity, and

denotes the modified speed of sound. The parameters \(a_4\) and \(a_5\) in Eq. (9.74) and Eq. (9.75) read in the general case

where \(a_1\) is defined in Eq. (9.58) and \(a_1^{\Gamma}\) in Eq. (9.64). In the case of a perfect gas and the definition of the preconditioning parameter \(\theta\) according to Eq. (9.69), the parameters become \(a_4 = \phi\) and \(a_5 = \phi c^2\), with \(\phi\) given by Eq. (9.67). When using the second definition of \(\theta\) from Eq. (9.71), the parameters in Eq. (9.76) simplify to \(a_4 = M_r^2\) and \(a_5 = M_r^2 c^2\), respectively. As we can see, \(c' \approx (V/2)\sqrt{5}\) for \(|M| \rightarrow 0\) and hence the eigenvalues become equalised as intended. In this way, the stiffness represented by the condition number Eq. (9.40) is reduced (the condition number is \(C_N \approx 2.6\)) and the convergence of the time-stepping or iterative solution process is dramatically enhanced. On the other hand, \(a_4 = 1\) and \(c' = c\) for \(|M| \ge 1\), and thus the eigenvalues of \(\bar{A}_c\) are recovered.

Based on Eq. (9.74), the spectral radius of the preconditioned system becomes

with \(\Delta S\) denoting the face area. This expression is employed to compute the time step in Eq. (6.14), (6.18), (6.20) or (6.22). It also replaces the spectral radius in Eq. (4.53) in the case of the central artificial dissipation (see also Eqs. (5.33), (9.52), (9.53)).

The left and right eigenvectors of the convective flux Jacobian (cf. Section A.11) will be changed by the preconditioning. This is of importance for spatial discretisations based on characteristic variables like Roe's upwind scheme (Subsection 4.3.3), or the upwind TVD scheme presented in the Subsection 4.3.4.

The matrix of the left eigenvectors of the preconditioned new flux Jacobian \(\Gamma^{-1}\bar{A}_{c,p}\) in Eq. (9.47) can be written as [91]

The matrix of the right eigenvectors of \(\Gamma^{-1}\bar{A}_{c,p}\) multiplied by \(\bar{\Gamma}\), as indicated in Eq. (9.55), is given by the formula [91]

where \(a_1^{\Gamma}\) is defined in Eq. (9.64). Further parameters in the above Eqs. (9.78) and (9.79) read

with \(a_4\), \(a_5\) defined in Eq. (9.76) and \(c'\) in Eq. (9.75), respectively.

The eigenvectors \(\bar{T}_{c,p}\) and \(\bar{T}_{c,p}^{-1}\) of the matrix \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\), which are required for the preconditioned Roe scheme formulated in the conservative variables (see Eq. (9.54)), can be obtained by a modification of the above eigenvectors Eqs. (9.78) and (9.79). It can be shown that the columns of the right-eigenvector matrix of \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\) are composed of \(\bar{P}\vec{x}\), where \(\vec{x}\) are the right eigenvectors of \(\Gamma^{-1}\bar{A}_{c,p}\). Hence, the matrices of the right and left eigenvectors of \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\) can be expressed as

Inserting the above relations (9.81) into the formula for the preconditioned Roe scheme in the conservative variables Eq. (9.54), we obtain

with \(\bar{\Gamma}\bar{T}_p\) given by Eq. (9.79) and \(\bar{T}_p^{-1}\) by Eq. (9.78), respectively.

An example of the application of the preconditioned finite-volume scheme from Eq. (9.51) with \(\bar{\Gamma}\) according to Eq. (9.65) to airfoil flow is presented in Figs. 9.10-9.12. As we can observe, both the central scheme as well as Roe's upwind scheme fail to deliver the correct solution for an inflow Mach number of 0.01. The schemes without preconditioning cannot predict the pressure distribution and hence the lift coefficient (\(C_L = 0.323\) and \(0.324\) versus the correct \(0.352\)). As demonstrated in Fig. 9.12, preconditioning helps to obtain the correct solution (\(C_L = 0.353\)), and it also significantly accelerates the convergence.