chapter. Acceleration Techniques 第9章 加速技术[cfd-000A]

The standard formulation of the implicit residual smoothing reads in 3D

where \(\vec{R}^{*}_{I,J,K}\), \(\vec{R}^{**}_{I,J,K}\), and \(\vec{R}^{***}_{I,J,K}\) denote the smoothed residuals in \(I\)-, \(J\)-, and \(K\)-direction, respectively. The parameters \(\epsilon^I\), \(\epsilon^J\), and \(\epsilon^K\) stand for the smoothing coefficients in the three computational coordinates. The implicit operator in Eq. (9.5) resembles second-order central difference. The term Central Implicit Residual Smoothing (CIRS) is therefore used. The implicit system in Eq. (9.5) is solved by the Thomas algorithm for the inversion of tridiagonal matrices.

The smoothing coefficients are usually defined as functions of spectral radii of the convective flux Jacobians [11]. The purpose is to apply only as much smoothing in each coordinate direction as it is necessary for stability and good error damping. A suitable formula for 2D was suggested in [12]

Here, \(\sigma^{*}/\sigma\) denotes the ratio of the CFL numbers of the smoothed and unsmoothed scheme. The variable \(r\) stands for the ratio of the convective spectral radii (Eq. (4.53)), i.e., \(r = \hat{\Lambda}_c^{J}/\hat{\Lambda}_c^{I}\). The parameter \(\Psi \approx 0.125\) ensures linear stability of the smoothing operation.

In 3D, the smoothing coefficients can be evaluated using an expression similar to Eq. (9.6) [12]

where

etc. The typical value of the parameter \(\Psi\) is 0.0625.

The maximum of the ratio of the CFL numbers \(\sigma^{*}/\sigma\) depends on the value of the smoothing coefficient and on the type of the spatial discretisation scheme. In the case of the central scheme (Subsection 4.3.1), the ratio is given by

In practice, value of \(\sigma^{*}/\sigma \approx 2\) can be reached (\(\epsilon = 0.8\)). Higher ratios reduce the damping of the time-stepping scheme. There is no such simple condition like (9.8) for the upwind spatial discretisation since the maximum of \(\sigma^{*}/\sigma\) depends also on the stage coefficients. But normally, the CFL number (see Tables 6.1 and 6.2) can be raised by about factor two for inviscid and up to factor five for viscous flows (using the hybrid scheme from Table 6.2). Figures 9.2 and 9.3 demonstrate the effect of the central IRS on the convergence for an inviscid and viscous flow on structured grids. The comparison is done in terms of the CPU time, in order to account for the additional numerical effort due to the IRS. As we can see, the solution time can be significantly reduced especially for a viscous flow. Even larger savings can be however realised when the IRS is employed together with multigrid (Section 9.4).

The limitation of the time step due to the viscous spectral radius \(\hat{\Lambda}_v\) in Eq. (6.18) can be offset with the aid of the implicit residual smoothing. The maximum time step is calculated according to the formula (6.14) without \(\hat{\Lambda}_v\). In flow regions where the viscous spectral radius dominates, smoothing is carried out using higher coefficients \(\epsilon^I\), \(\epsilon^J\), \(\epsilon^K\). Smoothing coefficients based on the viscous spectral radii can be computed from [12], [13]

where the constant \(C \approx 5/4\). The maximum of the coefficients from Eqs. (9.7) and (9.9) is then taken as the resulting smoothing coefficient.

CIRS is implemented on unstructured grids by applying the Laplacian operator (see Eq. (5.24)) to the residual. Thus, the smoothed residual \(\vec{R}^{*}_{I}\) in a control volume \(I\) is obtained from the implicit relation [14]

The sum includes all \(N_A\) adjacent control volumes. The relation (9.10) is solved for \(\vec{R}^{*}_{I}\) using Jacobi iteration. Useful values of the smoothing coefficient are \(0.5 \le \epsilon \le 0.8\). With these values of \(\epsilon\) it is possible to increase the CFL number (and hence the time step) two to five times. Due to the diagonal dominance of the matrix, the Jacobi iteration converges in about two steps.

The previously discussed CIRS method works satisfactorily for subsonic and transonic flows. It is also helpful in viscous dominated regions. However, CIRS exhibits poor error-damping characteristics in conjunction with upwind spatial discretisation schemes. Furthermore, the robustness of a multistage scheme accelerated by CIRS suffers in the case of strong shocks. Therefore, the so-called Upwind Implicit Residual Smoothing (UIRS) method was developed [15], [16], which is particularly suited for high Mach-number flows.

By contrast to CIRS, the UIRS methodology takes the sign of the convective eigenvalues \(\bar{\Lambda}_c\) into account. The idea is to smooth the residuals only in the direction of the characteristic of the Euler equations (i.e., \(dx/dt = \mathrm{const.} = \Lambda_c\)). This approach prevents unphysical influences on the upstream residuals. The UIRS method requires transformation of the residuals into the characteristic variables (cf. Appendix A.11). In this way, each component of the residual vector can be smoothed independently according to the sign of the corresponding eigenvalue. For the \(l\)-th component of the residual, the implicit operator is defined in 1D as [15], [16]

where \(\vec{R}^{c}\) denotes the residual transformed into characteristic variables, i.e.,

The smoothed residuals \(\vec{R}^{*}\) are obtained by the solution of a tridiagonal (bidiagonal if \(\Lambda_c^l\) does not change its sign) equation system using the Thomas algorithm. Afterwards, the residuals \(\vec{R}^{*}\) are transformed back into the conservative variables and the solution can be updated.

The most desirable feature of UIRS is that it leads to very favourable damping properties of the multistage scheme, particularly in connection with an upwind spatial discretisation. The ability to damp solution errors remains or even improves for high smoothing coefficients. It was demonstrated for 1-D Euler equations that values like \(\epsilon = 500\) and \(\sigma^{*} = 1000\) result in a stable and very fast explicit multistage scheme [15], [16] (see also in analysis/mstage on the CD-ROM). However, the problem is the implementation of UIRS in multiple dimensions. Since the convective flux Jacobians cannot be diagonalised simultaneously in all coordinate directions (see Appendix A.11), the transformation Eq. (9.12) and the smoothing Eq. (9.11) have to be carried out separately for each computational coordinate. The effect of the coordinate splitting is a reduced maximum smoothing coefficient to \(2 \le \epsilon \le 6\). Despite this, the convergence to the steady state is strongly accelerated as compared to CIRS [17], [16]. The largest improvements in terms of convergence speed and robustness occur in combination with multigrid [18], [16].

The smoothing coefficients in multiple dimensions are scaled by the eigenvalues. For example, in 2D the following formula can be employed

The relation between the CFL number of the smoothed scheme and the coefficient \(\epsilon\) reads

The constant \(C\) depends on the kind of the spatial discretisation. In the case of the central scheme \(C = 1\). For the 1st- or 2nd-order upwind scheme the value is \(C = 2\).

In order to circumvent the numerical effort of the transformation to the characteristic variables, a simplified version of the UIRS method was suggested [17], [18], [16]. Written in the \(I\)-direction it reads

The Mach number \(M\) is based on velocity projected into the direction of the particular computational coordinate (here: \(I\)). Due to the low operation count, the simplified UIRS is especially suitable for 3-D flow problems. In Ref. [16] it was demonstrated that the CPU time needed to reach the steady state can be halved as compared to CIRS (hypersonic flow past a blunt cylinder, \(M_\infty = 8\)).

Before the geometric multigrid scheme can be applied, the coarser grids have to be generated. The standard way is to coarsen the grid evenly in all coordinate directions. However, Mulder [49] proposed an approach called semicoarsening. Here, the grid is coarsened only in one direction, which is changed from one coarse level to another. The semicoarsening methodology is especially suited for flow problems, where the governing equations are stiff in one spatial direction. An example is the direction normal to the wall in boundary layers. Applications of semicoarsening to the Navier-Stokes equations were reported, e.g., in Refs. [26], [50], [51].

The discretised governing equations on the fine grid read in accordance with the relationship (6.1)

In the following, the finest grid will be denoted by subscript \(h\) in reference to the spacing of the grid lines (characteristic dimension of the control volume on unstructured grids). Starting from a known solution \(\vec{W}^{n}_{h}\), a new solution \(\vec{W}^{n+1}_{h}\) is obtained after one time step with some suitable iterative scheme. A new residual \(\vec{R}^{n+1}_{h}\) is evaluated with this solution. In order to improve the solution \(\vec{W}^{n+1}_{h}\) using a coarse grid, the following three steps are carried out:

The solution is transferred to the coarse grid by means of the interpolation

where the subscript \(2h\) denotes the coarse grid¹ and \(\hat{I}^{2h}_{h}\) is the interpolation operator. The residuals have to be transferred to the coarse grid as well, so that their low-frequency error components can be smoothed. A conservative transfer operator is employed for this purpose. This means that when the control volume size increases, the value of the residual must increase by the same amount. The residuals of the fine grid are also required in order to retain the solution accuracy of the fine grid on the coarse grid. For this purpose, a source term, the so-called forcing function [19], [22], is formed as the difference between the residual transferred from the fine grid and the residual computed using the initial solution \(\vec{W}^{(0)}_{2h}\) (Eq. (9.17)) on the coarse grid, i.e,

Here, \(I^{2h}_{h}\) represents the restriction operator which transfers residuals from the fine to the coarse grid. This type of multigrid scheme is known as the Full Approximation Storage (FAS) method [19]. The FAS method is particularly suited for non-linear equations because the nonlinearities in the system are carried down to the coarse levels through the re-discretisation.

The solution is evaluated on the coarse grid in the same way as on the fine grid. The forcing function in Eq. (9.18) is added to the residual of the coarse grid, i.e.,

Hence, the time-stepping scheme can be written in the form

In the case of the explicit multistage scheme (Subsection 6.1.1), this results in

in accordance with Eq. (6.5). It has to be noted that during the first iteration (stage in Eq. (9.21)), \((\vec{R}_F)_{2h}\) is identical to the residual transferred from the fine grid (i.e., \((\vec{R}_F)_{2h} = I^{2h}_{h}\vec{R}^{n+1}_{h}\) from Eq. (9.18)). This guarantees that the solution on the coarse grid depends on the residual of the fine grid and thus retains the accuracy of the fine grid.

An important question is the accuracy of the spatial discretisation scheme on coarse grids. Since the coarse grids do not influence the accuracy of the fine-grid solution, first-order schemes are sufficient. The advantages of first-order accurate discretisation on coarse grids are the increased robustness, better damping properties, and lower numerical effort in comparison to higher-order schemes.

After one or several time steps (iterations) were carried out on the coarse grid, the correction with respect to the initial - interpolated - solution (Eq. (9.17)) is computed. This so-called coarse grid correction is given by

The coarse-grid correction is interpolated to the fine grid in order to improve the solution there. Hence, the new solution on the fine grid reads

where \(I^{h}_{2h}\) is denoted as the prolongation operator.

The basic multigrid scheme described above consists of one coarse grid only. If multiple coarse grids are present, steps 1 and 2 are repeated until the coarsest grid is reached. It is important to realize that the forcing function on the coarse grids is formed from the restricted corrected residual of Eq. (9.19)). For example, on the coarse grid \(4h\), the forcing function is obtained from

In this way, the residual of the finest grid controls the accuracy of the solution on all coarse grids. After a given number of time steps on the coarsest grid, step 3 can be successively repeated until the finest grid is reached again. This procedure is known as a saw-tooth or V-cycle (see Fig. 9.4a). However, it is also possible to conduct more cycles on the coarse grids. This strategy, termed the W-cycle, is displayed in Fig. 9.4b. It is employed particularly frequently for transonic flows. In the case of supersonic and hypersonic flows, the V-cycle proved to be more efficient.

The optimum number of time steps before the restriction and after the prolongation depends on the type of the time-stepping scheme. In the case of the explicit multistage scheme (Subsection 6.1.1 or 6.1.2), it is common to carry out only one time step before the restriction of residuals and no time step after the prolongation. However, the robustness of the multigrid scheme can be improved by smoothing the coarse grid corrections (Eq. (9.22)) before adding them to the fine grid solution \(\vec{W}^{n+1}_{h}\) (Eq. (9.23)). The same central implicit smoothing (with constant coefficients) as described in Section 9.3 is utilised.

The other popular time-stepping method, the implicit LU-SGS scheme (see Subsection 6.2.4), requires two iterations before the restriction for the best multigrid efficiency [34]. The number of time steps after the prolongation depends on the spatial discretisation. In the case of the central scheme (Subsection 4.3.1), no time step is necessary [34]-[36], but the solution correction can be smoothed. On the contrary, one time step should carried out after the prolongation if an upwind spatial discretisation is used. This (2,1)-strategy proved to be an optimum with respect to robustness and computing time for various flow conditions [35], [36].

It should be pointed out that in practice the multigrid scheme is not started directly from the finest grid. Instead, several multigrid cycles are executed from one of the coarse grids. The approximate solution is interpolated to the next finer grid (using the same operator as for the prolongation), few more cycles are performed, the solution is again interpolated to the next finer grid and so on, until the finest grid is reached. In this way, a good starting solution is obtained on the finest grid with only a moderate numerical effort. This very efficient procedure is termed the Full Multigrid (FMG) method [19].

The restriction (9.18) and the prolongation (9.23) operator must fulfil certain accuracy requirements, namely [52]

where \(m_R\) and \(m_P\) denote the degree plus 1 of the polynomial, which is exactly interpolated by the restriction and the prolongation operator, respectively. For example, \(m_R\) or \(m_P\) are equal to two in the case of linear interpolation. Furthermore, \(m_E\) represents the order of the governing equations. Thus, \(m_E = 1\) for the Euler equations, and \(m_E = 2\) in the case of the Navier-Stokes equations. If the condition (9.25) is violated, the additional errors introduced by the restriction and/or prolongation will disturb the fine-grid solution. Hence, such multigrid scheme will converge only slowly or it will even diverge.

The implementation of multigrid on structured grids is straightforward since the coarse grids can be easily generated by deleting every second grid line in the respective coordinate direction. The spacing of the grid lines is therefore \(2h\), \(4h\), etc. This guarantees that the numerical effort on the coarse grids stays low as compared to the finest grid. Several representative examples are provided in Refs. [11]-[13], [20]-[22], [26], [32]-[36].

As we can conclude from Fig. 9.5, the operators for the solution interpolation, the restriction of residuals and the prolongation of corrections have to be defined differently for cell-centred and node-centred (cell-vertex) schemes. For example, two successive grids have every second grid point in common. On the contrary, the cell centres are always at different locations. Therefore, we shall discuss the standard forms of the transfer operators separately for the cell-centred and node-centred (identical to cell-vertex) finite-volume schemes.

Apart from the symmetrical, purely geometrically defined restriction and prolongation operators, which will be presented next, upwind-biased forms were suggested in [16], Chapter 4. Upwind restriction and prolongation, which account both for the characteristics of the flow equations, improve the robustness of the multigrid scheme for hypersonic flows. In order to save space, we shall discuss the implementation of an upwind prolongation operator for the node-centred spatial discretisation only.

The solution is transfered from the fine to the coarse grid by using a volume weighted interpolation. In 2D, Eq. (9.17) becomes (see Fig. 9.6a)

Similar transfer operator is employed in 3D, where the summation is over the eight fine-grid control volumes, which form one coarse-grid cell.

The restriction operator is defined as a sum of the residuals from all cells which are contained in one coarse-grid control volume. Hence, in 2D we have (cf. Fig. 9.6a)

and likewise in 3D. Those residuals \(\vec{R}^{n+1}_{h}\) in Eq. (9.27), which are located outside the physical domain, have to be set to zero.

The prolongation of the coarse-grid correction Eq. (9.23) can be conducted in two different ways. First, a zeroth-order prolongation operator results, if \(\delta\vec{W}_{2h}\) is equally distributed to all surrounding cell centres as indicated in Fig. 9.6b. Thus, for instance

The second possibility, which leads to a faster convergence of the multigrid scheme, consists of two steps. In a first step, \(\delta\vec{W}_{2h}\) is interpolated to the grid nodes like for the cell-vertex scheme (see below). In a second step, the nodal values are averaged to obtain the value in the centre of the fine-grid cell. Referring to Fig. 9.6c, the following final relationship can be derived

The corresponding expression in 3D can be found by a similar procedure. It reads

Because of the common nodes between the fine and coarse grid, the solution can be transfered simply by injection, i.e.,

The standard central restriction operator represents a linear interpolation from the nodes of all four (eight in 3D) fine-grid cells, which resemble one coarse-grid cell. According to Fig. 9.7, the restricted residual is computed in 2D as

In 3D, the fine-grid residuals are collected as follows

with the factors

In the above Eq. (9.34), only those indices are shown which are different from \(i, j, k\). It should be noted that the residuals \(\vec{R}^{n+1}_{h}\) must be set to zero at all physical boundary and dummy points before the restriction.

The prolongation of the coarse-grid correction can be implemented as a loop over the points of the coarse grid. Within the loop, the values \((\delta\vec{W}_{2h})_{i,j,k}\) are distributed to the fine-grid points using the same weights as for the restriction (cf. Fig. 9.7b). The particular contributions are summed up in order to obtain the complete transfered correction at each point of the fine grid.

In principle, upwind prolongation can be formulated either in characteristic or in conservative variables. A particularly efficient implementation in conservative variables was proposed in [16]. The methodology employs upwind-biased interpolation of the corrections in Eq. (9.23) according to the Mach number and the velocity direction. The numerical effort is very low, nevertheless the robustness of the multigrid scheme can be significantly improved for high Mach-number flows [16].

The interpolation of the solution corrections \(\delta\vec{W}_{2h}\) to the finer grid is accomplished in two steps. In the first step, the corrections at the points \(A\), \(B\), \(C\), and \(D\) (see Fig. 9.8), which are common to both grids, are transferred directly to the finer grid. In the second step, the corrections are interpolated to the points \(e\), \(f\), \(e'\), \(f'\) and \(g\), which are contained only on the finer grid. The interpolation depends on the sign and the absolute value of the Mach number at the corresponding point. The Mach number in this case is calculated using

is employed. The normal vector \(\vec{n}\) in Eq. (9.35) can be obtained either by averaging the face vectors of the control volume (in the direction \(A-B\)), or by normalising the vector from point \(A\) to point \(B\). Then, the correction is transferred to point \(e\) according to the rule

The same procedure applies also to the points \(f\) to \(f'\). The upwinding helps to match the information exchange between the grids better to the real physics. The interpolation to the point \(g\) is more difficult. In Ref. [16], the values at the surrounding points \(A\) to \(D\) were simply averaged

However, some sort of upwind weighted interpolation would be more appropriate. The upwind prolongation can be implemented in similar way also in 3D. Despite the simplification, encouraging results were obtained in a number of test cases [16].

As compared to the structured grids, the construction of the coarse grids is much more involved in the case of unstructured grids. The problem is how to construct an uniformly coarsened grid from a set of elements (grid cells) which have no particular ordering. Additionally, the ratio of the cell volumes of the coarse to the fine grid has to stay within a certain margins (about 4 in 2D and 8 in 3D). One possibility how to solve this problem is to apply the AMG methodology [41]-[48], which we briefly discussed at the beginning of Section 9.4. However, the geometric multigrid is still more widely used. Therefore, we shall concentrate here on this approach.

Three main methods for the generation of coarse grids can be identified:

• nonnested-grids approach,
• topological methods, and
• agglomeration of control volumes.

Reviews of the above methods were presented in Refs. [53] and [54].

The standard restriction and prolongation operators are based on purely geometrically defined interpolation. Leclerq and Stoufflet [55] suggested upwind transfer operators, which are particularly promising for flows with strong shocks. Their upwind restriction/prolongation is based on the transformation of the residuals/corrections into the characteristic variables. After an upwind-biased interpolation, the restricted/prolongated values are transformed back into the physical variables. Numerically much less expensive upwind multigrid method was presented in Ref. [16] (see also Subsection 9.4.3, Eq. (9.36)).

The most obvious idea is to generate a sequence of completely independent, increasingly coarser grids [56]-[61]. It is not necessary that the grids contain any common nodes. Therefore, we speak of nonnested grids. However, it is important that the main geometrical features (leading and trailing edge, fuselage nose, etc.) are retained on all coarse grids. This is not easy to accomplish, particularly in the case of a geometrically complex configuration. Multigrid based on nonnested grids is hardly used today.

One particular approach applies graph-based algorithms in order to remove certain nodes from the fine grid. The remaining nodes are then re-triangulated [62], [63]. On the contrary to the nonnested-grids approach, the interpolation between the grids becomes easier, since the successive grids contain common nodes. However, the method inherits the drawback of the nonnested grids with respect to geometry conformance.

A further topological method employs grid refinement [64], [29], [65]. The technique starts from a coarse grid and generates finer grids by element division. The methodology can be applied either over the whole physical domain or only locally (e.g., at boundary layers). The disadvantage of this approach is that the quality of the finest grid strongly depends on the initial coarse grid and the refinement procedure. The problem can be partially cured by edge swapping [66], [67].

Another idea for the generation of coarse grids is based on edge collapsing [68]. It was initially developed for inviscid flows on tetrahedral grids. The edge-collapsing method was further extended to viscous flows on mixed-element grids in [69].

A very efficient methodology for unstructured grids is the so-called agglomeration multigrid. It was first presented by Lallemand [70], Lallemand et al. [71] and by Koobus et al. [72]. Later on, the agglomeration multigrid was adopted by various authors [73]-[77], [31]. The method generates a coarse grid by fusing the control volumes of the finer grid with their neighbours. The resulting coarse grids consist of successively larger, irregularly shaped polyhedral cells. This is depicted in Fig. 9.9. As we can see, the agglomeration technique retains the full discretisation of the boundary surfaces. This represents a significant advantage over all previously discussed unstructured multigrid methods. However, it should be mentioned that up to now, implementations of the agglomeration multigrid were based mostly on the median-dual cell-vertex scheme (Subsection 5.2.2). The application of the agglomeration multigrid to a cell-centred scheme (Subsection 5.2.1) was described in [73].

The volume agglomeration for a node-centred scheme proceeds in the following steps:

• 1. build a list of the so-called seed points. Seed points are grid points selected to agglomerate the surrounding control volumes. The list of seed points can contain either those points which form an approximate maximal independent set [75], or simply all points of the current grid level.
• 2. Loop over all seed points.
• 3. If the seed point is unagglomerated, agglomerate all its nearest neighbours (connected by an edge), which were not already agglomerated.
• 4. Check the coarsening ratio (i.e., how many fine-grid control volumes are contained within a coarse-grid volume). If the ratio is less than four (eight in 3D), the neighbours of the already agglomerated nearest neighbours are added (if not associated with another seed point), until the optimum coarsening ratio is achieved. In Ref. [27], it was proposed to agglomerate those distance-two neighbours first, which are connected to at least two (three in 3D) agglomerated nearest neighbours.
• 5. If there are still seed points in the list, goto step 2.
• 6. Eliminate singletons. These are single control volumes which could not be agglomerated, because there were no unagglomerated neighbours. A singleton can be eliminated by agglomeration with such neighbouring control volume, which has the smallest coarsening ratio. This leads to coarse-grid levels with a more regular distribution of control-volume areas [31].

The above procedure is repeated until all coarse grids are generated.

The volume agglomeration has to start from the boundary in order to preserve grid isotropy. In 3D, user intervention may be required to prescribe the agglomeration direction depending on the shape of the boundary. To overcome this difficulty, Okamoto et al. [77] proposed another algorithm denoted as global coarsening. The method employs a global partitioning scheme, which is based on edge colouring. The partitioning scheme is used to generate an independent set of edges. In a second step, all control volumes, which share an edge of the independent set, are agglomerated. The procedure is repeated until the prescribed coarsening ratio is achieved. The method does not require the specification of an initial seed point or agglomeration direction. It can also treat any type of grid cells.

The implementation of agglomeration multigrid presents little difficulty for inviscid flows. The Euler equations are in general discretised as fluxes over individual control-volume faces. In this respect, it does not matter how complex is the shape of the control volume (in fact, averaged face vectors are used). Furthermore, a first-order accurate spatial scheme requires the knowledge of flow quantities in the neighbouring control volumes only. This information is readily available on the coarse grids.

In the case of viscous flows, the discretisation of the diffusive fluxes on arbitrary shaped control volumes is no longer straightforward. The problem is the evaluation of gradients at face midpoints (see the discussion in [31]). Another, and even more serious, difficulty is related to the required accuracy of the prolongation operator. The inequality in Eq. (9.25) suggests \(m_P = 2\), since the common restriction operator (sum of residuals) leads to \(m_R = 1\) only. However, the construction of a linear interpolation is not easy on the coarse grids.

In order to circumvent the construction of a first-order accurate prolongation operator, Mavriplis [53] proposed to use constant prolongation (i.e., all points of the fine grid contained within an agglomerated coarse-grid volume get the same solution correction) and a scaling of the viscous fluxes. However, this approach does not lead to optimal multigrid efficiency.

Haselbacher [31] suggested to retain the fine-grid discretisation of the viscous fluxes also on the coarse grids and to enforce the boundary conditions like on the finest grid. Moreover, he proposed a piecewise linear prolongation operator. For a scalar quantity \(U\), it can be written as

The gradient \((\nabla\delta U_{2h})_i\) in Eq. (9.38) is calculated by using the linear least-squares reconstruction described in Subsection 5.3.4, Eq. (5.55). The values of the limiter function \(\Psi_i\) are evaluated according to the Barth-Jespersen limiter function presented in Subsection 5.3.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.