4.3 Discretisation of the Convective Fluxes 对流通量的离散化[cfd-4-3]

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}\).