chapter. Structured Finite Volume Schemes 第4章 结构网格有限体积格式[cfd-0005]

In general, we consider flow in a plane as a special case of a 3-D problem, where the solution is symmetric with respect to one coordinate direction (e.g., to the \(z\)-direction). Because of the symmetry and in order to obtain correct physical units for volume, pressure, etc., we set the depth of all grid cells and control volumes equal to a constant value \(b\). The volume of a control volume results then in two dimensions from the product of its area with the depth \(b\). The area of a quadrilateral can be exactly calculated by the formula of Gauss. Hence, for a control volume like that displayed in Fig. 4.1a, we get after some algebra

In the above, we have assumed that the control volume is located in the \(x\)–\(y\)-plane and that the \(z\)-coordinate is the symmetry axis. Since the depth \(b\) is arbitrary, we may set \(b = 1\) for convenience. In two dimensions, the faces of a control volume are given by straight lines and therefore the unit normal vector is constant along them. When we integrate the fluxes according to the approximation of Eq. (4.2), we have to evaluate the product of the area of a face \(\Delta S\) and the corresponding unit normal vector \(\vec{n}\) - the face vector \(\vec{S}\)

Because of the symmetry, the \(z\)-component of the face vectors (and of the unit normal vector) is zero. It is therefore dropped from the expressions. The face vectors of the control volume from Fig. 4.1a are given by the relations

The unit normal vector at face \(m\) is then obtained from Eq. (4.5) as

with

In practice, only the face vectors \(\vec{S}_1\) and \(\vec{S}_4\) are computed and stored for each control volume \(\Omega_{I,J}\). The face vectors \(\vec{S}_2\) as well as \(\vec{S}_3\) are taken (with reversed signs to become outward facing) from the appropriate neighbouring control volumes in order to save memory and to reduce the number of point operations.

As opposed to the previous 2-D case, the calculation of face vectors and volumes poses some problems in three dimensions. The main reason for this is that, in general, the four vertices of the face of a control volume may not lie in a plane. Then, the normal vector is no longer constant on the face (Fig. 4.2). In

order to overcome this difficulty, we could decompose all six faces of the control volume into two or more triangles each. The volume itself could then be built of tetrahedra. Performing this subdivision in an appropriate manner would lead to a discretisation scheme which is at least first-order accurate on arbitrary grids [3]. Of course, the numerical effort would be increased substantially, because the fluxes would have to be integrated over each partial triangle separately. Hence, the number of point operations would be at least doubled. However, in [3], [4] it is shown that for reasonably smooth grids, where the control volume faces approach parallelograms, the decomposition into triangles does not noticeably improve the solution accuracy. Therefore, we shall employ a simplified treatment of the quadrilateral faces in the following considerations, which is based on an averaged normal vector.

A face vector \(\vec{S}\) of an hexahedral control volume, like that rendered in Fig. 4.1b, is most conveniently computed using the same Gauss's formula as employed in 2D for the area of a quadrilateral. Thus, e.g., for the face \(m = 1\) (points 1, 5, 8 and 4 in Fig. 4.1b) we first define the differences

The face vector \(\vec{S}_1 = \vec{n}_1\Delta S_1\) results then from

The five remaining face vectors are calculated in a similar manner. It is again very convenient to store only three of the six the face vectors (e.g., \(\vec{S}_1\), \(\vec{S}_3\), and \(\vec{S}_5\)) for each control volume \(\Omega_{I,J,K}\). The remaining face vectors \(\vec{S}_2\), \(\vec{S}_4\) as well as \(\vec{S}_6\) are obtained (with reversed signs to become outward facing) from the appropriate neighbouring control volumes. The above expressions in Eq. (4.8) and (4.9) deliver an average face vector. The approximation becomes exact when the face approaches a parallelogram, i.e., when the vertices of the face lie all in one plane. The unit normal vector is obtained from Eq. (4.7) with

Various, more or less accurate formulae are available for the calculation of the volume of a general hexahedron (see, e.g., [1]). One approach, which performed very well in various applications, is based on the divergence theorem [5]. This relates the volume integral of the divergence of some vector quantity to its surface integral. The key idea is to use the location in space of some point of the control volume \(\Omega\), let us call it \(\vec{r} = [r_x, r_y, r_z]^T\), as the vector quantity. Herewith, the divergence theorem reads

We can easily evaluate the left-hand side of Eq. (4.11) which gives us the volume of \(\Omega\) that we are looking for

If we assume now the unit normal vector is constant on all faces of the control volume, we can solve the surface integral on the right-hand side of Eq. (4.11) as follows

In Eq. (4.13), \(\vec{r}_{\mathrm{mid},m}\) denotes the midpoint of the control volume face \(m\). For example,

where the vectors \(\vec{r}_1\), \(\vec{r}_5\), \(\vec{r}_8\), and \(\vec{r}_4\) correspond to the vertices 1, 5, 8, and 4 of the face \(m = 1\) in Fig. 4.1b. Similar relations hold for the midpoints of the remaining faces. The area \(\Delta S_m\) of the face \(m\) in Eq. (4.13) is obtained from Eq. (4.10). Combining Equations (4.12) and (4.13) together, and inserting the face vector \(\vec{S}\) for the product \(\vec{n}_m\Delta S_m\), we have finally the relationship

for the volume of the control volume \(\Omega_{I,J,K}\).

The origin of the coordinate system can be in principle moved to any place without affecting the volume calculation in Eq. (4.14). This leads us to the advice to locate the origin in one vertex of the control volume (e.g., point 1 in Fig. 4.1b), in order to achieve a better scaling of the numerical values. Thus, we may replace \(\vec{r}\) in the above expressions (4.11)–(4.14) by a transformed vector \(\vec{r}^{*}\), which is defined as

It is important to note that the volume computed with aid of Eq. (4.14) is exact for a control volume with planar faces.

We speak of a cell-centred scheme if the control volumes are identical with the grid cells and if the flow variables are located at the centroids of the grid cells as indicated in Fig. 4.3. When we evaluate the discretised flow equations (4.2), we have to supply the convective and the viscous fluxes at the faces of a cell [6]. They can be approximated in one of the three following ways:

• by the average of fluxes computed from values at the centroids of the grid cells to the left and to the right of the cell face, but using the same face vector (generally applied only to the convective fluxes);
• by using an average of variables associated with the centroids of the grid cells to the left and to the right of the cell face;
• by computing the fluxes from flow quantities interpolated separately to the left and to the right side of the cell face (employed only for the convective fluxes).

Thus, taking the cell face \(\vec{n}_{I+1/2,J}\) in Fig. 4.3 as an example, the first approach - average of fluxes - reads in two dimensions

with \(\Delta S_{I+1/2,J}\) computed from Eqs. (4.6) and (4.7).

The second possible approach - average of variables - can be formulated as follows

where the conservative/dependent variables at the face \(\vec{n}_{I+1/2,J}\) of the control volume are defined as the arithmetic average of values at the two adjacent cells, i.e.,

The flux vector \(\vec{F}\) in Eq. (4.16) stands either for the convective or for the viscous fluxes.

The third methodology starts with an interpolation of flow quantities (being mostly velocity components, pressure, density and total enthalpy) separately to both sides of the cell face. The interpolated quantities - termed the left and the right state (see the begin of Section 4.3) - differ in general between both sides. The fluxes through the cell face are then evaluated from the difference of the left and right state using some non-linear function. Hence,

where

represent the interpolated states. Of course, similar relations like (4.15)–(4.19) hold also for the other cell faces.

The same approximations are employed in three dimensions. For example, at the cell face \(\vec{n}_{I+1/2,J,K}\) (e.g., identical to \(\vec{n}_2\) in Fig. 4.1b) the average of fluxes in Eq. (4.15) becomes

with \(\Delta S_{I+1/2,J,K}\) being defined correspondingly to Equations (4.8) and (4.9). The average of variables reads similarly to Eq. (4.16) as

with

The way over the interpolation of the flow variables results similarly to Eq. (4.18) in

where \(\vec{U}_L\) and \(\vec{U}_R\) are the interpolated values at the cell face.

The last term in the discretised flow equations (4.2) which remains to be evaluated is the source term \(\vec{Q}\). As we already stated in the introduction, the source term is usually supposed to be constant inside the control volume. For this reason, it is calculated using the flow variables from the corresponding cell centre. Hence, we may define

Using the above relations, the fluxes through the faces can be computed and the numerical integration over the boundary of \(\Omega_{I,J,K}\) may be performed according to (4.2). In other words, the complete residual \(\vec{R}_{I,J,K}\) is obtained. In Sections 4.3 and 4.4, we shall learn more about the details of the evaluation of the convective and viscous fluxes.

In a cell-vertex scheme, all flow variables are associated with the nodes of the computational grid. Within the approach based on overlapping control volumes, the grid cells still represent the control volumes, just as in the case of the cell-centred scheme. The difference is that now the residuals computed for the control volumes have to be distributed to the grid points [4], [7], [8]. The situation is sketched in Fig. 4.4.

Let us consider the control volume \(\Omega_{I,J}\) in Fig. 4.4, which is defined by the nodes

Note that the point \((i,j)\) is located at the lower left corner of \(\Omega_{I,J}\). The convective fluxes, e.g., for the face \(\Delta S_{I,J-1/2}\), which is given by the points \((i,j)\) and \((i+1,j)\), are approximated as

The variables at the midpoint of the face are evaluated using an arithmetic average of the variables at the nodes defining the face, i.e.,

The face area \(\Delta S_{I,J-1/2}\) is computed from the relations (4.6) and (4.7).

The approach remains the same in three dimensions. For example, for the face associated with the normal vector \(\vec{n}_3\) in Fig. 4.1b, the averaged variables read

If we assume the edge 1-2 being oriented in the \(i\)-direction, edge 1-5 in the \(j\)-direction, edge 1-4 in the \(k\)-direction, and if we finally associate the point \((i,j,k)\) with the corner 1 in Fig. 4.1b, the average in Eq. (4.27) can also be written as

The convective flux is then again obtained from

where the face area \(\Delta S_{I,J,K-1/2}\) results from the formulae (4.8) and (4.9). The viscous fluxes are normally computed employing the same approach as for the dual control-volume scheme [9], [10], which results in a more compact scheme (i.e., one which involves fewer nodal values).

Summing up all face contributions given by relations (4.25), (4.26) and (4.28), (4.29), respectively, we obtain intermediate residuals \(\vec{R}_{I,J,K}\) for all grid cells. In order to relate the cell-based to the node-based residuals, a further approximation is made using a residual distribution formula. It is basically a function, which evaluates the unknown node-based residual from a weighted sum of all cells having the particular grid node in common. The following distribution formulae were devised:

• volume weighted sum due to Ni [7];
• non-weighted sum due to Hall [8];
• characteristic (upwind) weighting procedure of Rossow [11], [12].

Theoretical investigations of the truncation error [4] suggest that Ni's scheme is more accurate than Hall's approach. However, in practice Ni's distribution formula leads to problems in places, where the grid is strongly distorted and stretched. For example, strong oscillations of the pressure field near the trailing edge of an airfoil were observed when using O-grids [13], [4]. Furthermore, convergence could only be achieved when the numerical viscosity was increased considerably. The underlying idea of the upwind weighting procedure [11], [12] is quite similar to that of the fluctuation-splitting schemes [14]-[17] (cf. Subsection 3.1.5), but the implementation is numerically much simpler. Basically, the residuals are sent only upstream in the characteristic direction.

Of the three approaches, Hall's distribution scheme proved to be the most robust. In Hall's scheme, the residual at a particular node results from a simple sum of all intermediate residuals \(\vec{R}_{I,J,K}\), which cells share the node. Thus, in the 2-D case rendered in Fig. 4.4, we get

In three dimensions, in total eight cell-based residuals have to be summed up in the same way. A close inspection of (4.30) reveals that \(\vec{R}_{i,j}\) is just the net flux through the boundary of the supercell

due to the fact that the fluxes across the inner faces cancel each other. The supercell also represents the "total" control volume, centred at the point \((i,j)\). In the 3-D case, the total volume consists of the cells

As it can be seen from Fig. 4.4, the control volumes overlap by at least one cell, which gave the scheme its name.

The source term is calculated using the flow variables from the corresponding grid node, i.e.,

With the above definitions, the time-stepping scheme in Eq. (4.2) becomes

This represents a system of ordinary differential equations, which has to be solved in each grid point \((i,j,k)\) in the same way as for the cell-centred scheme. Note that the total volume (4.31) or (4.32), respectively, is utilised in Eq. (4.34).

The central scheme with artificial dissipation is very simple compared to other discretisation methods. It is easy to implement with either the cell-centred scheme or with both cell-vertex schemes. For these reasons the scheme became very wide-spread.

The basic idea of the central scheme is to compute the convective fluxes at a face of the control volume from the arithmetic average of the conservative variables on both sides of the face. Since this would allow for odd-even decoupling of the solution (generation of two independent solutions of the discretised equations) and for overshoots at shocks, artificial dissipation (which is similar to the viscous fluxes) has to be added for stability. The scheme was first implemented for the Euler equations by Jameson et al. [6]. Because of the names of the authors, it is also abbreviated as the JST scheme.

The central scheme is generally less accurate in the resolution of discontinuities and boundary layers than, let say, the upwind schemes. However, it is computationally considerably cheaper. Therefore, attempts were made to improve the accuracy of the scheme, while still keeping the numerical effort low. For example, improvements were devised to reduce the amount of the artificial dissipation in boundary layers [31], [32], or to enhance the shock resolution [33]-[35]. Another idea, followed in [35], is to utilise the Jacobian matrix of the convective fluxes, in order to scale the dissipation independently for each conservation equation. This successful approach is known as the matrix dissipation scheme. It can be viewed as a compromise between the original scalar scheme and the upwind schemes. The basic scheme also employs a single, pressure-based sensor to switch from second- to first-order accuracy at discontinuities to prevent non-physical oscillations of the flow variables. In [36], Jameson developed a concept called the SLIP (Symmetric Limited Positive) scheme, where a limiter is applied separately for each conservation equation. A brief description of the previous approaches and comparisons for inviscid and viscous 2-D flows was presented in [37].

The convective fluxes (Eq. (2.21)) through a face of the control volume are approximated using the average of variables, according to the Equations (4.16), (4.21), (4.36), or (4.41), respectively. Artificial dissipation is then added to the central fluxes for stability [6], [38]. Thus, the total convective flux at face \((I+1/2,J,K)\) reads

where the flow variables are averaged as (see also Fig. 4.8)

In the case of the cell-vertex scheme with dual control volumes, node indices \((i,j,k)\) would be used instead. For simplicity, \((I+1/2,J,K)\) will be abbreviated as \((I+1/2)\) hereafter. The artificial dissipation flux consists of a blend of adaptive second- and fourth-order differences, which result from the sum of first- and third-order difference operators

From Eq. (4.50) we can see that the scheme possesses a compact 9-point stencil in two dimensions and a 13-point stencil in three dimensions. The dissipation is scaled by the sum of the spectral radii of the convective flux Jacobians in all coordinate directions

The spectral radius at the cell face \((I+1/2)\), e.g., in I-direction (represented by the superscript \(I\)), results from the average

It is evaluated using the formula

where \(V\) stands for the contravariant velocity (2.22) and \(c\) for the speed of sound, respectively. A pressure-based sensor is used to switch off the fourth-order differences at shocks, where they would lead to strong oscillation of the solution. The sensor also switches off the second-order differences in smooth parts of the flow field, in order to reduce the dissipation to the lowest possible level. Herewith, the coefficients \(\epsilon^{(2)}\) and \(\epsilon^{(4)}\) in Eq. (4.50) are defined as

with the pressure sensor given by

Typical values of the parameters are \(k^{(2)} = 1/2\) and \(1/128 \le k^{(4)} \le 1/64\). In order to reduce the amount of artificial dissipation across a viscous shear layer, we can re-define the scaling factors in Eq. (4.50) as follows [31], [32]

These are then employed instead of Eq. (4.51). The directionally dependent coefficients \(\phi\) are given by the relations

The parameter \(\sigma\) is usually set equal to 1/2 or 2/3. This formulation decreases the scaling of the dissipation terms in the direction along the shorter side of a control volume, whose longer side is aligned with the flow.

In order to improve the accuracy by reducing the numerical dissipation, the preceding JST scheme can be modified to become more like an upwind scheme. The idea is to use a matrix - the convective flux Jacobian - instead of the scalar value \(\hat{\Lambda}^{S}\) to scale the dissipation terms [35]. In this way, each equation is scaled properly by the corresponding eigenvalue. Hence, the Eq. (4.50) becomes

The scaling matrix corresponds to the convective flux Jacobian \((\bar{A}_c = \partial\vec{F}_c/\partial\vec{W})\) diagonalised with absolute values of the eigenvalues

The matrices of right \((\bar{T})\) and left \((\bar{T}^{-1})\) eigenvectors as well as the diagonal matrix of the eigenvalues \(\bar{\Lambda}_c\) can be found in the Appendix A.11. The eigenvalues must be limited at stagnation points and sonic lines to prevent the dissipation from becoming zero. An efficient way of computing the product of \(|\bar{A}_c|\) with \(\vec{W}\) is provided in [35]. It should be noted that by setting \(\epsilon^{(2)} = 1/2\) and \(\epsilon^{(4)} = 0\), we obtain a first-order accurate, fully upwind scheme.

As it was already stated before, the idea here was to develop a scheme which accuracy is close to that of upwind schemes, but which is still computationally only slightly more expensive (about 15-20%) than the scalar dissipation approach. Results of comparisons with flux-vector splitting schemes (CUSP and AUSM) were recently reported in [37].

The flux-vector splitting methods can be viewed as the first level of upwind schemes, since they account only for the direction of wave propagation. The flux-vector splitting schemes decompose the vector of the convective fluxes into two parts - either according to the sign of certain characteristic variables, or into a convective and a pressure part. The well-known Van Leer's flux-vector splitting scheme [39] belongs to the first category based on characteristic decomposition. The second approach is followed by more recent methods like the Advection Upstream Splitting Method (AUSM) of Liou et al. [40], [41], or the Convective Upwind Split Pressure (CUSP) scheme of Jameson [42], [43], respectively. Further similar approaches are the Low-Diffusion Flux-Splitting Scheme (LDFSS) introduced by Edwards [44], or the Mach number-based Advection Pressure Splitting (MAPS) scheme of Rossow [45], [46].

The flux-vector splitting schemes can be implemented only in the framework of the cell-centred scheme (Subsection 4.2.1), or the cell-vertex scheme with dual control volumes (Subsection 4.2.3). Their advantage can be seen in only a moderately increased numerical effort but a much better resolution of shocks and boundary layers, as compared to the central scheme with scalar artificial dissipation. However, the matrix dissipation scheme (Eq. (4.58)), can also produce results of comparable accuracy [37]. Because of certain numerical difficulties, many modifications to the basic schemes (particularly to AUSM) were devised by various researchers and the development still continues. In the following, we shall present the basics of the Van Leer, AUSM and CUSP schemes, give some hints with respect to the most important modifications, and provide references to the corresponding literature.

Van Leer's flux-vector splitting scheme [39] is based on characteristic decomposition of the convective fluxes. An extension of the approach to body-fitted grids was presented in [47], [48]. The convective flux is split into a positive and a negative part, i.e.,

according to the Mach number normal to the face of the control volume (e.g., at \((I+1/2)\) - see Fig. 4.8)

where \(V\) represents the contravariant velocity (2.22) and \(c\) the speed of sound, respectively. In the case of the cell-vertex scheme with dual control volumes, the cell indices have to be replaced by node indices. The values of the flow variables \(\rho\), \(u\), \(v\), \(w\), and \(p\), respectively, have to be interpolated first to the faces of the control volume correspondingly to Eq. (4.19) or Eq. (4.39). Then, the positive fluxes are computed with the left state and the negative fluxes with the right state. The advection Mach number \((M_n)_{I+1/2}\) is obtained from the relation [39]

where the split Mach numbers are defined as

and

The Mach numbers \(M_L\) and \(M_R\) are evaluated using the left and right state, respectively, i.e.,

In the case of \(|M_n| < 1\) (subsonic flow), the positive and the negative flux parts are given by

The mass and energy flux components are defined as

For supersonic flow, i.e., for \(|M_n| \ge 1\), the fluxes are evaluated from

The evaluation of the left and right state follows generally the MUSCL approach [29], which is given by Eqs. (4.46). The higher order schemes \(\hat{\kappa} = -1\), \(\hat{\kappa} = 0\) and \(\hat{\kappa} = 1/3\), respectively, require a limiter if the flow field contains discontinuities like shocks. More details are provided in Subsection 4.3.5.

The flux-vector splitting scheme of Van Leer performs very well in the case of the Euler equations. But several investigations [49], [50], carried out with the Navier-Stokes equations revealed that splitting errors in the momentum and the energy equations smear the boundary layers and also lead to inaccurate stagnation and wall temperatures. A modification to the momentum flux in the direction normal to the boundary layer was therefore suggested in Ref. [51]. A similar remedy for the energy flux was proposed in Ref. [52]. Both modifications together remove the splitting errors, and hence they improve the solution accuracy considerably [53].

The Advection Upstream Splitting Method (AUSM) was introduced by Liou and Steffen [40], and Liou [54]. It was subsequently modified by Wada and Liou [55] and renamed as AUSMD/V. Finally, an improved version termed AUSM+ was presented by Liou [41], [56].

The underlying idea of the approach is based on the observation that the vector of convective fluxes (2.21) consists of two physically distinct parts, namely the convective and the pressure part

The first term in Eq. (4.69) represents scalar quantities, which are convected by the contravariant velocity \(V\). By contrast, the pressure term is governed by the acoustic wave speed. The idea now is to discretise the convective term in purely upwind manner by taking either the left or the right state, depending on the sign of \(V\) (even for subsonic flow). On the other hand, the pressure term includes both states in the subsonic case. It becomes fully upwind only for a supersonic flow.

Following the basic AUSM from [40], we introduce an advection Mach number \((M_n)_{I+1/2}\) from Eq. (4.61). Herewith, we can recast the convective flux at the face \((I+1/2)\), or \((i+1/2)\) of the control volume, respectively, into

where

Similar to Van Leer's flux-vector splitting scheme, the advection Mach number is evaluated as a sum of the left and right split Mach numbers according to the relations (4.62) and (4.63)-(4.65). The computation of the left and right state (flow quantities: \(\rho\), \(u\), \(v\), \(w\), \(p\), \(H\)) is based again on a separate interpolation to the faces of the control volume according to Eq. (4.19) or Eq. (4.39). The interpolation follows the MUSCL methodology [29], as it is given in Eq. (4.46). All higher-order MUSCL schemes (\(\hat{\kappa} = -1\), \(\hat{\kappa} = 0\), and \(\hat{\kappa} = 1/3\)) require a limiter, if the flow field contains strong gradients like shocks. Please refer to Subsection 4.3.5 for more details.

The pressure at the face \((I+1/2)\) of the control volume is obtained from the splitting [40]

with the split pressures given by [39]

and

It is also possible to use the following lower-order expansion for \(|M_{L/R}| < 1\)

It should be noted that we can write AUSM also in the form

The first term on the right-hand side of the above Eq. (4.76) represents a Mach number-weighted average of the left and right state - similar to the average of fluxes Eq. (4.15), (4.20) or Eq. (4.35), (4.40), respectively. The second term has a dissipative character. It is scaled by the scalar value \(|(M_n)_{I+1/2}|\).

AUSM proved to deliver a crisp resolution of strong shocks and accurate results for boundary layers. However, the original AUSM [40], [54] was found to generate local pressure oscillations at shocks and in cases where the flow is aligned with the grid [57]. In [57], [58] it was therefore suggested to switch at shocks to Van Leer's scheme. When the advection Mach number \((M_n)_{I+1/2}\) tends to zero, the dissipation term in Eq. (4.76) will approach zero as well. Thus, any disturbances cannot be damped by the scheme. In order to solve the flow alignment problem, it was proposed in [57] to modify the scaling of the dissipation term as follows

where \(\delta\) is a small value \((0 < \delta \le 0.5)\). Hence, there will always be a sufficient amount of numerical dissipation. In order to retain the accuracy of AUSM for boundary layers, the parameter \(\delta\) could be reduced in the wall normal direction using the same idea as given for the central scheme by Eq. (4.57). Further improvements of the basic AUSM, with respect to better behaviour in the vicinity of shocks, was presented in [41], [56] as AUSM+. The modifications consist of new Mach and pressure splittings, which replace the relations (4.63), (4.64) and (4.73), (4.74), respectively.

The concept of the Convective Upwind Split Pressure (CUSP) scheme is quite similar to that of AUSM. However, the CUSP approach has the advantage to be formulated as an average of fluxes (but without weighting like within AUSM) minus a dissipation term. This feature is crucial for the implementation in an explicit, hybrid multistage scheme. Furthermore, because of the different scaling factors as compared to AUSM, the CUSP scheme behaves more favourably in the case of flow alignment. The CUSP scheme was introduced by Jameson [42], [59], [60], and subsequently modified by Tatsumi et al. [43], [61]. It can be implemented either within the cell-centred or the (cell-vertex) dual control-volume type of spatial discretisation.

The convective fluxes (Eq. (2.21)) through a face of the control volume are approximated using the arithmetic average of fluxes according to the Equations (4.15), (4.20), (4.35), or (4.40), respectively. The dissipation term is then subtracted from the central fluxes for stabilisation. Thus, the total convective fluxes at the face \((I+1/2)\) read

In the case of the dual control-volume discretisation, \((i+1/2)\) would apply instead. The dissipation term, which is composed of a linear combination of the differences of the state and the flux vector, can be expressed as

There are two choices for the term \(\phi\) in Eq. (4.79). Either, \(\phi\) is set equal to the total energy \(E\), or it is set to the total enthalpy \(H\). In the first case we speak of the E-CUSP scheme [62], the second choice is consequently called the H-CUSP scheme. The formulation with \(\phi = H\) preserves the total enthalpy [43] and is therefore suitable for inviscid flows. The left (\(L\)) and right (\(R\)) state is evaluated similarly to the MUSCL approach [29], using the limited interpolation (4.118)-(4.121). The two factors \(\alpha^{*}c\) and \(\beta\) in Eq. (4.79) are defined as

and

with \(M_n = V/c\). The contravariant velocity \(V\) in the above formulae (4.80) and (4.81) is computed from the arithmetic mean of the velocity components, i.e.,

The speed of sound \(c\) in the evaluation of \(M_n\) is also obtained from arithmetically averaged quantities. The positive and negative eigenvalues \(\Lambda^{+}\) and \(\Lambda^{-}\) in Eqs. (4.80), (4.81) are those of the so-called Roe matrix [63], which will be discussed in the Subsection 4.3.3. The eigenvalues are given by [60]

where \(\tilde{V}\) denotes the contravariant velocity, \(\gamma\) is the ratio of specific heat coefficients and \(\tilde{c}\) stands for the speed of sound. The flow variables in Eq. (4.83), which have to be evaluated at the faces of the control volumes, are obtained using the so-called Roe averages [63]

The factors \(\alpha^{*}c\) and \(\beta\) are defined such that full upwinding of the convective fluxes results for supersonic flow, i.e., \(\alpha^{*}c = 0\) and \(\beta = \mathrm{sign}(M_n)\). On the other hand, in subsonic flow (when \(\beta = 0\)) the dissipation is scaled by \(|V|\). This is a desirable property for the computation of viscous layers. In cases of large aspect ratio cells, explicit time-stepping schemes usually require increased numerical dissipation in the direction of the longer cell side in order to stay robust. This can be accomplished by employing ratios of the spectral radii similar to Eq. (4.57). More details can be found in Ref. [37], which also contains comparisons between the CUSP scheme and the scalar as well as the matrix artificial dissipation (Subsection 4.3.1) scheme.

The implementation of the CUSP scheme can be considerably simplified by replacing the Roe averages from Eq. (4.84) by arithmetic averages and the eigenvalues from Eq. (4.83) by \(\Lambda^{\pm} = V \pm c\). Introducing the definition

the factors of the modified scheme read [60], [43]

and

The contravariant velocity in \(M_n = V/c\) is still computed as indicated in Eq. (4.82). The parameter \(\delta\) in Eq. (4.86) is intended to prevent the dissipation from disappearing at stagnation points, but this does not seem to be always necessary. The above simplifications lead to a computationally very efficient scheme, however the shock resolution is slightly reduced as compared to the original formulation with Roe averages.

The flux-difference splitting schemes evaluate the convective fluxes at a face of the control volume from the (in general discontinuous) left and right state by solving the Riemann (shock tube) problem. The idea was first introduced by Godunov [64]. In contrast to the flux-vector splitting schemes, the flux-difference splitting considers not only the direction of wave (information) propagation, but also the waves themselves. In order to reduce the computational effort of Godunov's scheme for the exact solution of the Riemann problem, approximate Riemann solvers were developed, e.g., by Osher et al. [65] and by Roe [63]. In particular, Roe's method is applied quite often because of its high accuracy in boundary layer flows and good resolution of shocks. Therefore, the Roe solver shall be presented in more detail in the following subsection.

Roe's approximate Riemann solver can be implemented either in the framework of the cell-centred scheme or the dual control-volume scheme. It is based on the decomposition of the flux difference over a face of the control volume into a sum of wave contributions, while ensuring the conservation properties of the Euler equations. On the face \((I+1/2)\) or \((i+1/2)\), respectively, the difference is expressed as [63]

In the above Eq. (4.88), \(\bar{A}_{Roe}\) denotes the so-called Roe matrix, and \(L\) or \(R\) the left and right state (see Fig. 4.8), respectively. The Roe matrix is identical to the convective flux Jacobian \(\bar{A}_c\) (see Appendix A.9), where the flow variables are replaced by the so-called Roe-averaged variables. The flux difference in Eq. (4.88) is exact, if the Roe's averages are computed from the left and the right state by the following formulae [63], [66]

We can make the decomposition into waves in Roe's scheme clearer when we insert the diagonalisation of the Roe matrix, i.e., \(\bar{A}_{Roe} = \bar{T}\bar{\Lambda}_c\bar{T}^{-1}\), into the Eq. (4.88)

The matrix of left \((\bar{T}^{-1})\) and right \((\bar{T})\) eigenvectors, as well as the diagonal matrix of eigenvalues \((\bar{\Lambda}_c)\) are evaluated using Roe's averaging (4.89). In the above Eq. (4.90), the characteristic variables \(\vec{C}\) represent the wave amplitudes, the eigenvalues \(\bar{\Lambda}_c\) are the associated wave speeds of the approximate Riemann problem, and finally the right eigenvectors are the waves themselves.

Following from the previous discussion, the convective fluxes are evaluated at the faces of a control volume faces according to the formula [63]

The product of \(|\bar{A}_{Roe}|\) and the difference of the left and right state can be efficiently evaluated as follows

where

The jump condition is defined as \(\Delta(\bullet) = (\bullet)_R - (\bullet)_L\) and the Roe-averaged variables are given in Eq. (4.89), respectively.

The left and the right state are determined using the MUSCL scheme [29], which is given in Eq. (4.46). All higher-order schemes (\(\hat{\kappa} = -1\), \(\hat{\kappa} = 0\), and \(\hat{\kappa} = 1/3\)) have to be supplemented by limiters (Subsection 4.3.5), if the flow field contains any discontinuities.

Because of the formulation in Eq. (4.88), Roe's approximate Riemann solver will produce an unphysical expansion shock in the case of stationary expansion, for which \((\vec{F}_c)_L = (\vec{F}_c)_R\) but \(\vec{W}_L \ne \vec{W}_R\). Furthermore, the so-called carbuncle phenomenon may occur, where a perturbation grows ahead of a strong bow shock along the stagnation line [67], [68]. See also the discussion in Ref. [69]. The underlying difficulty is that the original scheme does not recognise the sonic point. In order to solve this problem, the modulus of the eigenvalues \(|\bar{\Lambda}_c| = |\tilde{V} \pm \tilde{c}|\) is modified using Harten's entropy correction [70], [71]

where \(\delta\) is a small value, which can be conveniently set equal to some fraction (e.g., 1/10) of the local speed of sound. In order to prevent the linear waves \(|\Delta\vec{F}_{2,3,4}|\) from disappearing for \(\tilde{V} \rightarrow 0\) (e.g., at stagnation points or for grid-aligned flow), the above modification can also be applied to \(|\tilde{V}|\). A clear disadvantage of the Roe solver as compared to the central scheme or to the flux-vector splitting schemes shows up for a real gas simulation. Namely, the Roe matrix and averaging have to be changed correspondingly, which may become quite complicated. The reader may find examples of formulations for equilibrium as well as non-equilibrium real gas flows in [72]-[75] and in the references cited therein. An implementation of Roe's scheme for arbitrary compressible and incompressible fluids was recently described in Ref. [76].

The idea of Total Variation Diminishing (TVD) schemes was first pursued by Harten [77]. The TVD schemes are based on a concept aimed at preventing the generation of new extrema in the flow solution. The principal condition for a TVD scheme is that the total variation of the solution, defined as

for a scalar conservation equation, decreases in time. This implies that maxima in the solution must be non-increasing and minima non-decreasing. Hence no new local extrema may be created during the time evolution. Thus, a discretisation methodology with TVD properties allows it to resolve strong shock waves accurately, without any spurious oscillations of the solution, as they are for example generated by the central scheme with scalar or matrix artificial dissipation (Subsection 4.3.1).

The TVD schemes are implemented as an average of the convective fluxes combined with an additional dissipation term (flux-limited dissipation), which complies with the TVD conditions [77], [78]. If the dissipation term depends on the sign of the characteristic speeds, we speak of a symmetric TVD scheme [79], [80], otherwise of an upwind TVD scheme [81]-[85]. Experience shows that the upwind TVD scheme offers higher accuracy than the symmetric TVD scheme [86]. The upwind TVD scheme is particularly suitable for the simulation of supersonic and hypersonic flow fields [87]. It is also capable of accurate resolution of boundary layers [53], especially if the modification described in Ref. [88] is applied.

In this framework, the convective fluxes through the face \((I+1/2)\) of the control volume (see Fig. 4.8) can be expressed as

In the case of the cell-vertex scheme with dual control volumes (Subsection 4.2.3), the indices would read \((i+1/2)\), \((i+1)\), etc. The matrix \(\bar{T}\) contains the right eigenvectors of the Jacobian \(\bar{A}_c = \partial\vec{F}_c/\partial\vec{W}\). The entries of the matrix can be found in the Appendix A.11. In Equation (4.98), the term \(\vec{\Theta}\) takes account of the direction of the characteristic speeds. It controls the upwind direction of the difference operator. The l-th component of the vector \(\vec{\Theta}\) is defined as (cf. [84])

where \(\Lambda^{l}\) represents the individual eigenvalues of the diagonal matrix \(\bar{\Lambda}_c\) (see Appendix A.11), and \(\Psi\) the limiter function (Eq. (4.122)), respectively. Furthermore,

and finally \(\Delta C^{l}\) are the elements of the difference of characteristic variables, i.e.,

with \(\bar{T}^{-1}\) being the matrix of left eigenvectors. The so-called entropy correction of Harten [70], [71], i.e.,

prevents the value \(\psi(z)\) from vanishing for \(|z| \rightarrow 0\). The parameter \(\delta_1\) is best formulated as function of the velocity components and the speed of sound [84]

where \(0.05 \le \delta \le 0.5\). Values of the primitive variables at the face \((I+1/2)\) are obtained either from Roe's (4.89) or from simple arithmetic averaging of the states at \(I\) and \((I+1)\). The limiter function \(\Psi\), which prevents the generation of spurious solutions near strong gradients, will be presented in the next subsection. It should be stressed that the above upwind TVD scheme does not employ the MUSCL approach to achieve higher order accuracy.

One can show that the upwind TVD method is precisely of first-order in space when the limiter function \(\Psi\) in Eqs. (4.99), (4.100) is set equal to zero [89], which happens at discontinuities. Otherwise, the upwind TVD scheme, as presented above, is second-order accurate in smooth flow regions.

Second- and higher-order upwind spatial discretisations require the use of so-called limiters or limiter functions in order to prevent the generation of oscillations and spurious solutions in regions with large gradients (e.g., at shocks). Hence, what we are looking for is at least a monotonicity preserving scheme. This means that maxima in the flow field must be non-increasing, minima non-decreasing, and no new local extrema may be created during the time evolution. Or in other words, if the initial data is monotone then the solution has to remain monotone. The rather stringent conditions for monotonicity preserving schemes (or the more rigorous ones for TVD schemes) are often given up in favour of the Local Extremum Diminishing (LED) conditions [60]. Here, a local extremum contained only within the stencil has to decrease.

However, due to Godunov's theorem there is no possibility for a higher-order linear scheme (such as the MUSCL approach) to be monotonicity preserving [90]. It is therefore necessary to employ non-linear limiter functions in order to construct a monotonicity preserving or a TVD discretisation. This is demonstrated in Fig. 4.10, where the upwind TVD scheme of Eq. (4.98) was used with and without a limiter to compute 2-D transonic flow past the NACA 0012 airfoil. It can be clearly seen that without limiter, the solution exhibits large oscillations in the neighbourhood of the shocks on the upper and the lower side of the airfoil. On the other hand, the limited and the unlimited solutions become nearly identical away from the shocks.

The purpose of a limiter is to reduce the slopes (i.e., \((U_{I+1} - U_I)/\Delta x\)) used to interpolate a flow variable to the face of a control volume in order to constrain the solution variations. At strong discontinuities, the limiter has to reduce the slopes to zero to prevent the generation of new extrema. This implies for the MUSCL approach as well as for the TVD schemes that the (monotone) first-order upwind scheme (\(\epsilon = 0\) in Eq. (4.46)) is recovered in the immediate vicinity of large gradients. The last requirement to be imposed on a limiter is quite obvious - the original unlimited discretisation has to be obtained in smooth flow regions, in order to keep the amount of numerical dissipation as low as possible. The effect of a limiter on the interpolation of the left and right states is sketched in Fig. 4.11. The example shows the slope reduction at the local minimum at \(I\) and the change of the slope at the cells \((I+1)\), \((I+2)\) to achieve a monotone solution. It is important to realise that a difference between the left and right state at a face may (and generally will) still be present.

In the following, we shall describe four different limiter functions, which are well-established and proven in practice. We shall consider limiters for the second-order MUSCL, for the CUSP and for the upwind TVD scheme.

Van Leer's MUSCL approach [29] is turned into a monotonicity preserving scheme by employing a limiter function to reduce the differences \(\Delta_{+}U_I\) and \(\Delta_{-}U_I\) in Eq. (4.47) when necessary. Introducing slope limiters \(\Phi^{\pm}\), the MUSCL interpolation formulae in Eq. (4.46) are modified as follows (see also Fig. 4.8)

The parameter \(\epsilon\) in Eq. (4.46) was set equal to unity. The slope limiters are functions of the ratios of the consecutive solution variations, i.e., \(\Phi^{\pm}_{I+1/2} = \Phi(r^{\pm}_{I+1/2})\), with the definitions [1]

If we substitute now \(r_L\) for \(r^{+}_{I-1/2}\) and \(r_R\) for \(r^{-}_{I+3/2}\), thus

we can write Eq. (4.104) in the form

The above relationships Eq. (4.107) can be simplified if we consider only slope limiters with the symmetry property

With this definition, the limited MUSCL interpolation Eq. (4.104) becomes [91]

with the limiter function defined as

Different formulations of the slope limiter \(\Phi\) in Eq. (4.110) are now possible, which can be tailored to specific values of \(\hat{\kappa}\) to give the most accurate but stable and monotonicity preserving MUSCL scheme.

One particularly suitable combination for the second-order, upwind-biased scheme with \(\hat{\kappa} = 0\) is [92]

In this case, the function \(\Psi(r)\) corresponds to the Van Albada limiter [93]

and we obtain with Eq. (4.109) the following expressions for the left and right state

The function \(\delta\) is formally identical for both states. It reads

The coefficients \(a\) and \(b\) are defined for the left and right state as

and the difference operators \(\Delta_{\pm}\) are given by Eq. (4.47). The additional parameter \(\epsilon\) in Eq. (4.114) prevents the activation of the limiter in smooth flow regions due to small-scale oscillations [92]. This is sometimes necessary in order to achieve a fully converged steady-state solution. The parameter \(\epsilon\) is conveniently set proportional to the local grid scale, in 3D for example to \(\Omega^{1/3}\) [92], [94]. Additional scaling of the parameter \(\epsilon\) is required if the particular state variable \(U\) is given in physical units. It can be shown that the relations in Eq. (4.113) are identical to the original (unlimited) MUSCL scheme (4.46) with \(\hat{\kappa} = 0\) for smoothly varying flow. Thus, the accuracy of the solution is not influenced. On the other hand, the function \(\delta\) becomes zero at local extrema, reducing the accuracy to first order as desired.

Another limiter function was devised for the three-point, second-order accurate upwind-biased MUSCL scheme with \(\hat{\kappa} = 1/3\). Here, the slope limiter is given by

In this case, the function \(\Psi(r)\) corresponds to the limiter of Hemker and Koren [95]. Following the same way as in the previous case, we obtain formulae for the left and right state at the face \((I+1/2)\) which are identical to Eq. (4.113), but now with [92]

The definitions of the coefficients \(a\), \(b\), and of the parameter \(\epsilon\) are retained.

In the framework of the CUSP scheme (Subsection 4.3.2), the left (\(L\)) and right (\(R\)) states are evaluated to second-order accuracy according to [43]

where

In the above Eqs. (4.118) and (4.119), \(U\) represents a dependent variable and \(L()\) the limited average

respectively. The limiter itself is defined as

where \(\sigma\) is a positive coefficient which is usually set equal to two. The constant \(\epsilon\) is required to prevent division by zero (e.g., \(\epsilon = 10^{-20}\)). If \(\Delta_1\) and \(\Delta_2\) happen to have opposite sign but the same magnitude, the limiter becomes \(\Psi = 0\). This means that we obtain only a first-order accurate approximation for the left and the right state.

It should be mentioned that it is also possible to employ the MUSCL scheme with \(\hat{\kappa} = 0\) and the Van Albada limiter from Eqs. (4.113)-(4.115) instead of the above relations.

In comparison to the previous cases, the limiter here acts not on the conservative or the primitive variables, but on the characteristic variables \(\vec{C}\). One particularly suitable limiter function is given by [84]

where the \(\Delta C^{l}_{I+1/2}\) represents the difference of the characteristic variables at face \((I+1/2)\) of the control volume (Eq. (4.101)). The positive constant \(\epsilon \approx 10^{-20}\) in the denominator prevents division by zero. In regions with high gradients, the limiter function becomes zero, which leads with Eq. (4.99) and Eq. (4.98) to first-order accurate upwind scheme. The upwind TVD scheme in Eq. (4.98) retains second-order accuracy in areas with smoothly varying flow variables, where \(\Psi^{l}_{I} = C^{l}_{I} - C^{l}_{I-1}\).

In order to apply Green's theorem, which relates the volume integral of the first derivative to the surface integral of \(U\), we have to define a suitable control volume first. Since we need the derivatives at the midpoints of the faces for the summation in Eq. (4.2), we construct an auxiliary control volume centred at the face by connecting the midpoints of the edges defining adjacent grid cells [31], [19], [98] as shown in Fig. 4.12a. In order to evaluate the first derivative at the face \((I+1/2)\) - marked by a diamond symbol in Fig. 4.12a - we have to integrate the corresponding flow variable \(U\) over the boundary of the auxiliary control volume (denoted by the superscript ' in the following). Thus, e.g., for the derivative in the x-direction

where \(N_F\) stands for the number of faces (\(N_F = 4\) in 2D and \(N_F = 6\) in 3D). The volume \(\Omega'\) and the components of the face vector \(\vec{S}'_m = [S'_{x,m}, S'_{y,m}, S'_{z,m}]^T\), respectively, are computed as already presented in Section 4.1. The face values \(U_m\) are obtained either directly as cell-centred values (i.e., \(U_{i,j}\) and \(U_{i+1,j}\) on the left and the right face), or by averaging like on the upper and the lower face, e.g., at \(J+1/2\)

We can apply the same approach in three dimensions, where again four cell-centred values can be utilised for the averaging. Hence,

The above scheme is quite compact, with the computational stencil extending over only nine cells in two dimensions and over 15 in three dimensions. It should be noted that this approach for computing the first derivatives cannot suppress the generation of two types of spurious modes (decoupled solutions at neighbouring cell centres) [99], [19]: the chequer-board mode, arising from the form of the integral around the control volume, and a pair of corrugated or washboard modes, arising from the averaging of values in neighbouring cells. However, there are generally no difficulties with this in practice. A more serious problem would occur, if the gradients would be first evaluated for each cell (similar to the convective fluxes) and then averaged at the cell faces. Although this approach may appear more attractive than the current methodology, it is not recommended since it leads to strong odd-even decoupling.

A disadvantage of above scheme is a loss of accuracy if the grid is not uniform [98], [99]. Namely, for arbitrarily stretched grids the approximation of the derivatives becomes inconsistent (the centroid of the auxiliary control volume does no longer correspond to the face centre). Thus, the viscous fluxes are discretised with second-order accuracy only for moderately and smoothly stretched grids. Finally, it should be noted that the TSL approximation of the Navier-Stokes equations (Subsection 2.4.3) can easily be realised by omitting the appropriate contributions when computing the gradients. For example, if the boundary layer would be oriented along the \(I\)-direction in Fig. 4.12a, contributions from the left \((I, J)\) and the right side \((I+1, J)\) of the auxiliary control volume would be dropped.

As already mentioned, both types of cell-vertex schemes resort to the dual control volume (Subsection 4.2.3) for the discretisation of the viscous fluxes. Hence, the question is how to evaluate the first derivatives at the faces of this control volume. Considering Fig. 4.12b, one possible alternative is to calculate the gradients at the cell centres first by integrating over the grid cells, which yields first-order accuracy on arbitrarily stretched grids. In a next step, the cell-based gradients are averaged at the faces of the control volume \(\Omega\) like in Refs. [31], [100]. However, this approach cannot prevent an odd-even decoupling of the solution.

Another possibility, similar to the cell-centred scheme, is to construct an auxiliary control volume around the face by connecting the midpoints of the edges defining adjacent grid cells [101], [102]. This is depicted in Fig. 4.12b. The evaluation of the first differences proceeds along the same lines as discussed for the cell-centred scheme, with averaged quantities where necessary. It should be noted that this approach is formally identical to the finite difference approximation [96]-[98]. This scheme leads to first-order accurate discretisation of the viscous fluxes on arbitrarily stretched grids and to second-order accuracy on smooth grids [96], [98]. Another positive feature is that the computational stencil is confined to only nine nodes in two dimensions and to 15 nodes in three dimensions.

Finally, one further approach should be mentioned, where a more complex integration path was chosen, with averaging incorporating all neighbouring nodes [20]. A serious disadvantage of this scheme is that it encompasses a 25-point stencil even in two dimensions, which adds in general more numerical diffusion than compact stencils. Furthermore, if an implicit scheme would be envisioned for the time integration, the bandwidth of the flux Jacobian would become prohibitively large. A detailed discussion of various methodologies for the gradient evaluation can also be found in Ref. [103].