5.3 Discretisation of the Convective Fluxes 对流通量的离散[cfd-5-3]

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 according to Eq. (5.20). Since this would lead to odd-even decoupling of the solution (generation of two independent solutions of the discretised equations) and wiggles at shocks, artificial dissipation has to be added for stability. The artificial dissipation is based on a blend of second- and fourth-order differences. The scheme was first implemented for the Euler equations on structured grids by Jameson et al. [43]. Because of the names of the authors, it is also abbreviated as the JST scheme.

The implementation of the JST scheme on unstructured grids utilises the Laplacian operator for the second-order differences and the Laplacian of Laplacian for the fourth-order differences [44], [27]. In order to reduce the computational cost, pseudo-Laplacians are employed instead of true Laplacians. For this purpose, a 2-D formulation was proposed first in [45] and then improved in [46]. Later on, the scheme was extended to 3D in [2]. It makes use of a distance-weighting procedure. In this way, the scheme leads to a vanishing pseudo-Laplacian for a linearly varying function on any grid. Applied to a general scalar quantity \(U\) in cell \(I\), the pseudo-Laplacian assumes the form

where \(N_A\) stands for the number of adjacent control volumes. The cell indices have to be substituted by node indices \((i,j)\) in the case of the median-dual scheme. The sum in Eq. (5.24) is best evaluated using either a loop over faces (cell-centred scheme) or a loop over edges (median-dual scheme) similar to the flux computation. The geometrical weights \(\theta\) are defined as

and result from the solution of an optimisation problem [2]. The optimisation problem is solved by means of Lagrange multipliers. Herewith, the geometrical weights are obtained from the expression

where \(x,y,z\) are the Cartesian coordinates of the cell centroids (nodes in the case of the median-dual scheme). The Lagrange multipliers \(\lambda\) are computed for each cell (node) and follow from [2]

with the coefficients

Written for a cell \(I\), the first-order moments read

Furthermore, the second-order moments are given by

The geometrical weights (5.25) can lead to a non-positive approximation of the Laplacian and hence to a lost of stability on severely distorted grids. Therefore, clipping the weights to the range \((0,2)\) was suggested in [45]. However, this measure impairs the accuracy of the discretisation. See also the discussion in Ref. [12] for further details.

The fourth-order differences are evaluated as the Laplacian of the Laplacian, i.e., \(L(U)\) is substituted for \(U\) in Eq. (5.24). Hence, the final form of the artificial dissipation term for a cell \(I\) is

With the artificial dissipation term added, the system of equations in Eq. (5.2) becomes

where \(N_F\) denotes the number of the faces of the control volume (which may differ from the number of adjacent control volumes, e.g., if a quadrilateral face is divided into two triangles).

The second- and the fourth-order terms in Eq. (5.31) are scaled by the spectral radius of the convective flux Jacobian. According to Ref. [28], the spectral radius for cell \(I\) can be evaluated as

where \(V_m\) represents the contravariant velocity (2.22) and \(c_m\) the speed of sound, respectively. Both quantities are computed from flow variables averaged at the face. The spectral radius at the face of the control volume is obtained from

A pressure-based sensor is used to switch off the fourth-order differences at shocks and the second-order differences in smooth portions of the flow field. Herewith, the coefficients \(\epsilon_{IJ}^{(2)}\) and \(\epsilon_{IJ}^{(4)}\) in Eq. (5.31) 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\).

As we already discussed in Subsection 4.3.1, the accuracy of the above central scheme can be improved when we substitute a matrix [47] for the spectral radius \((\hat{\Lambda}_c)_{IJ}\) in Eq. (5.31). The implementation of this so-called matrix dissipation scheme on unstructured grids proceeds in the same way as on structured grids, with the scaling matrix defined as in Eq. (4.59). Application of the matrix dissipation scheme to 3-D mixed grids is discussed, e.g., in Ref. [48].

It is important to note that for elements other than triangles/tetrahedra, the popular explicit Runge-Kutta type of temporal discretisation experiences severe stability problems when it is coupled to the central scheme [49]. The reason is the representation of the fourth-order differences by the Laplacian of the Laplacian. A remedy is to employ a difference of the left and the right state (cf. Section 4.3) for the approximation of the fourth-order differences [49], i.e.,

This approach leads on quadrilateral/hexahedral grids to the same stencil as the corresponding structured scheme. The left and right state are computed using, e.g., the linear reconstruction described below in Subsection 5.3.3.

Upwind schemes seem to have gained, at least for the moment, much more popularity on unstructured grids than the above central scheme. In fact, the flux-difference splitting scheme of Roe [50] is the most widely employed approach on unstructured grids. It is the considerably more accurate resolution of boundary layers and the lower sensitivity to grid distortions in comparison to the central scheme, which explains the attractiveness of Roe's scheme. However, the price to be paid for the improved performance is the higher computational effort, which becomes quite significant if a limiter has to be used to suppress oscillations of the solution (Subsection 5.3.5).

Any of the upwind schemes presented in Section 4.3 for structured grids are applicable to unstructured grids without modifications to the basic methodology. Only the computation of the left and right state (Eq. (5.22)), which is denoted as solution reconstruction, as well as the evaluation of the limiting function require new formulations. For this reason, only the solution reconstruction and the limiters are discussed here. For details on the various upwind methods, the reader is referred to Subsections 4.3.2-4.3.4. An example for the implementation of Roe's scheme on unstructured grids using the median-dual approach can be found, e.g., in Ref. [33].

As we saw in Subsections 4.3.2-4.3.4, upwind schemes require flow states to be specified on the left and the right side of a control volume face. The same holds also for the modified artificial dissipation scheme in Eq. (5.37).

As a first approach, we can assume that the solution is constant inside each control volume. The left and right state are then simply the flow variables computed for the left and the right control volume. For example, in the case of the median-dual scheme (Fig. 5.9) we would set

with \(U\) representing some scalar flow variable. This leads to a spatial discretisation which is only first-order accurate. For viscous flows, a first-order accurate solution is too diffusive and leads to an excessive growth of shear layers. Therefore, methods with higher spatial accuracy are a must for the computation of viscous flows.

We can achieve second- and higher-order accuracy if we assume the solution changes over the control volume. For second-order accurate methods, which are the most commonly employed higher-order methods, the solution is assumed to vary in a linear fashion over the control volume. In order to compute the left and right state, a reconstruction of the assumed solution variation becomes necessary. In what follows, we shall discuss the most popular approaches for the reconstruction of linear and quadratic variations. The interested reader is referred to [38] for a comparison of various linear reconstruction techniques.

One possibility to achieve second-order accuracy consists of the extension of the MUSCL approach [51] to unstructured grids. When applied to the median-dual scheme, the method generates for each edge \(ij\) two "phantom" nodes \(i'\) and \(j'\) [52]-[56]. These phantom nodes are located at the endpoints of the line obtained by extending the edge \(ij\) by its length in both directions as sketched in Fig. 5.12. After the solution is interpolated from the surrounding elements (gray coloured in Fig. 5.12) to the phantom nodes, we can evaluate the left and right state using the MUSCL formulae Eq. (4.46). Hence,

with forward \((\Delta_{+})\) and the backward \((\Delta_{-})\) difference operators defined as

The MUSCL interpolation (5.39) has to be enhanced by a limiter function (according to Subsection 4.3.5) in the case of strong discontinuities. A disadvantage of this methodology is the necessity to store for each edge the elements which contain the phantom nodes. A further conceptual disadvantage is that no unique gradient is reconstructed for a control volume. Furthermore, difficulties can arise at boundaries, where one of the phantom points lies outside the physical domain.

Barth and Jespersen presented in [30] a reconstruction method, which is closely related to the finite element schemes. Here, it is assumed that the solution is piecewise linearly distributed over the control volume. Then, we can find the left and right state for a cell-centred scheme from the relations

where \(\nabla U_I\) is the gradient of \(U\) (\(= [\partial U/\partial x,\, \partial U/\partial y,\, \partial U/\partial z]^T\)) at the cell centre \(I\) and \(\Psi\) denotes a limiter function (cf. Subsection 5.3.5), respectively. The vectors \(\vec{r}_L\) and \(\vec{r}_R\) point from the cell-centroid to the face-midpoint, as indicated in Fig. 5.13a.

The same approach applies to the median-dual scheme [30], i.e.,

According to Fig. 5.9 or Fig. 5.13b,

represents the vector from node \(i\) to node \(j\).

It can be easily seen that the method of Barth and Jespersen corresponds to a Taylor-series expansion around the neighbouring centres/nodes of the face, where only the linear term is retained. The linear reconstruction is formally second-order accurate on regular grids [38]. The scheme reconstructs a linear function exactly on any grid, provided the gradient \(\nabla U\) is evaluated without an error. The linear reconstruction is likely the most popular one among the reconstruction methods.

The above scheme requires the computation of gradients at cell centres or at nodes, respectively. This can be accomplished either by the Green-Gauss or the least-squares approach, which are presented below in Subsection 5.3.4. Furthermore, the implementation of the limiter function on unstructured grids is described in detail in Subsection 5.3.5.

It was demonstrated by Frink [26] that for the cell-centred scheme the linear reconstruction (5.41) does not require an explicit evaluation of the gradient on purely triangular or tetrahedral grids. The reason are two invariant geometric features of these elements. First, a line from a node through the cell-centroid will always intersect the midpoint of the opposing face. Second, the distance from the cell-centroid to the face-midpoint is one-fourth (one-third for a triangle) of the distance from the face-midpoint to the opposite node. Thus, the gradient at the cell centre can be approximated by a simple finite difference [26]. For example, if we were to reconstruct the solution at the face-midpoint \(F_3\) in Fig. 5.7a, the formulae (5.41) would become

with \(U_C\) being the value at the cell-centroid, \(U_1\), \(U_2\), etc. denoting the nodal values, and finally \(\Psi\) standing for a limiter, respectively.

Two different ways were devised by Frink in order to determine the nodal values. The first approach is based on inverse distance weighting. Here, the contribution to a node from the surrounding cells is inverse proportional to the distance from the node to the cell-centroid [26], [57], i.e.,

with the weights \(\theta_{iJ} = 1/r_{iJ}\). The distance is computed from

The subscripts \(J\) and \(i\) refer to the cell-centroid and to the node, respectively. The above methodology leads to a reconstruction which is less than second-order accurate. However, Frink pointed out that no limiter is needed at least for inviscid flows [26], which reduces the computational effort significantly.

The second approach is based on work of Holmes et al. [45] and Rausch et al. [46] in 2D. It was later extended to 3D by Frink [2]. Here, the weights \(\theta_{iJ}\) in Eq. (5.45) are defined such that the nodal values are computed exactly if the variation is linear. This leads to the same constraints as for the computation of the pseudo Laplacian (5.24). Consequently, the weights are also the same and follow from the Equations (5.25)-(5.30). The coordinates \(x_I\), \(y_I\), \(z_I\) of the cell-centroids are just replaced by the node coordinates \(x_i\), \(y_i\), \(z_i\). The scheme is formally second-order accurate because the nodal values are computed exactly for a linear function. In order to assure positivity on distorted grids, the weights have to be restricted to the range \((0,2)\) [45]. Unfortunately, this reduces the accuracy of the reconstruction. Frink et al. [4] also recently reported some anomalous behaviour of the reconstruction for the Navier-Stokes equations.

In order to achieve higher than second-order accuracy with a polynomial reconstruction, we have to keep further terms in the truncated Taylor-series expansion around the neighbouring cell-centres/nodes of the face. Based on the work of Barth and Frederickson [58], Barth developed the concept of k-exact reconstruction scheme [59], i.e., a reconstruction exact for a polynomial of degree \(k\). The polynomial in Barth's method is defined in a way which guarantees the conservation of the mean, or in other words, the average of the reconstruction polynomial is equal to the mean solution in the control volume. This property assures the conservation of mass, momentum, and energy during the reconstruction. The method was implemented for \(k = 3\) in a median-dual scheme. The coefficients of the polynomial were computed using a least-squares approach. Similar ideas were followed for the cell-centred scheme by Mitchell and Walters [60], and by Mitchell [61]. However, these methods require a prohibitively high numerical effort and a complex data structure which prevented their widespread use.

Delanaye and Essers [62] and Delanaye [63] developed a particular form of the quadratic reconstruction for cell-centred schemes, which is computationally more efficient than the method of Barth. The left and right state are approximated using Taylor series truncated after the quadratic term [62], [63]

In the above Eq. (5.47), \(\bar{H}_I\) denotes the Hessian matrix, i.e.,

evaluated at the cell-centroid \(I\). The variables \(\Psi_{I,1}\) and \(\Psi_{I,2}\) represent two different limiter functions for the linear and the quadratic term [62], respectively. The quadratic reconstruction method is third-order accurate on regular grids and at least second-order accurate on arbitrary grids due to cancellation of error terms [63]. Necessary conditions for achieving these properties are, however, that the gradient \(\nabla U\) in Eq. (5.47) is evaluated at least with second-order and the Hessian with first-order accuracy. This is accomplished by combining Green-Gauss gradient evaluation with least-squares based approximation of the second derivatives [62], [63], which leads to a numerically efficient scheme. But the memory and the CPU-time overheads are still quite significant in comparison to the linear reconstruction. The method utilises a fixed stencil composed of face and node neighbours. The stencil is shown in Fig. 5.14 together with the integration path employed for the Green-Gauss gradient computation. In order to determine all coefficients of the quadratic polynomial, at least six (ten in 3D) values must be provided by the stencil. To maintain the accuracy provided by the quadratic reconstruction, it is necessary to consider a linear variation of the solution over the face instead of a constant value. This implies that the solution must be reconstructed at two points - the so-called Gauss quadrature points (cf. Fig. 5.14) - of a 2-D face (at three points of a triangular face) and that the fluxes have to be integrated in a piecewise manner over the face of the control volume [37].

An open point, which remains from the discussion of the piecewise linear and the quadratic reconstruction, is the determination of the gradient. Gradients of the velocity components and the temperature are also required for the evaluation of the viscous fluxes (Section 5.4). Two approaches will be presented in the following: the first is based on the Green-Gauss theorem, and the second utilises the least-squares method.

This method approximates the gradient of some scalar function \(U\) as the surface integral of the product of \(U\) with an outward-pointing unit normal vector over some control volume \(\Omega'\), i.e.,

Barth and Jespersen [30] derived a particular discretisation of the Green-Gauss approach from the Galerkin finite element method. Later on, the discretisation was extended to 3D by Barth [64]. Barth and Jespersen applied Eq. (5.49) to the region formed by the union of the elements meeting at a node. They proved that the approach can be formulated such that it becomes compatible with the edge-based data structure. However, this works only for the median-dual scheme on triangular/tetrahedral grids. The resulting formula reads

Here, \(\Omega'\) in Eq. (5.49) equals to the volume of the median-dual control volume \(\Omega\). The summation extends over all \(N_F\) edges incident to node \(i\). Furthermore, \(\vec{n}_{ij}\) denotes the average unit normal vector according to Eq. (5.23), and \(\Delta S_{ij}\) is the total face area, respectively. The same formula (5.50) is applicable in two or in three dimensions. It is important to mention that the summation has to be changed at boundaries in order to obtain a consistent approximation [65] (see also Section 8.10).

We can use the Green-Gauss method in the cell-centred scheme as well. Hence, the gradient at some cell-centroid \(I\) can be obtained from

where the summation extends over all faces of the cell with the volume \(\Omega\). In Eq. (5.51), \(\vec{n}_{IJ}\) denotes the unit normal vector and \(\Delta S_{IJ}\) the face area, respectively.

The main attractiveness of the Green-Gauss gradient evaluation by Eq. (5.50) or (5.51) is its similarity to the computation of the fluxes (e.g., Eq. (5.19)). This means that no additional data structures are needed for the reconstruction of gradients. The main disadvantage is that the approximation in Eq. (5.50) or Eq. (5.51), respectively, fails on mixed grids. It was demonstrated in [49] that the gradient can become highly inaccurate, particularly where different element types meet. We can solve the problem in the case of the median-dual scheme when we keep the volume \(\Omega'\) in Eq. (5.49) identical to the union of all cells incident to node \(i\). Referring to the situation sketched in Fig. 5.8, the gradient results then in 2D from

with \(i = 0\), the number of outer faces \(N_O = 6\) (given by the points \(P_1\)-\(P_6\)), and \(j+1 = 1\) for \(j = 6\). Furthermore, \(\vec{n}_j\) and \(\Delta S_j\) stand for the unit normal vector and the area of the outer cell faces. In 3D, we can use

when we assume all faces are triangles - either naturally or by decomposition. The same remedy can be also employed for the cell-centred scheme. The surface of the control volume \(\Omega'\) is then defined by the centroids of the distant-one and distant-two neighbouring cells [62], [63], as it is rendered in Fig. 5.14. The gradient is computed correspondingly to Eq. (5.52) or (5.53) with cell instead of node indices.

The clear disadvantage of such a cell-based approach is the necessity of an additional data structure, which provides a link between the central node/centroid and the outer faces of \(\Omega'\). Thus, the approach is no longer grid-transparent (i.e., independent of the cell information) and an efficient gather-scatter loop is no longer possible. This renders the least-squares technique, which is described below, more attractive on mixed grids.

Using the edge-/face-based implementation (5.50) or (5.51) on triangular or tetrahedral grids, and the cell-based methodology (5.52) or (5.53) on mixed grids, the Green-Gauss approach is at least first-order accurate [63]. It is also consistent, i.e., the gradient of a linear function is computed to roundoff error. First-order accuracy is sufficient for the linear reconstruction. Second-order accuracy on arbitrary grids, which is required for the quadratic reconstruction, can be achieved by subtracting an estimate of the truncation error from the first-order approximation of the gradient [63].

The evaluation of gradients by the least-squares approach was first introduced by Barth [64], [36]. In order to illustrate the method, let us consider the median-dual scheme. Herewith, the least-squares approach is based upon the use of a first-order Taylor series approximation for each edge which is incident to the central node \(i\). The change of the solution along an edge \(ij\) can be computed from

where \(\vec{r}_{ij}\) is given by Eq. (5.43) and represents the vector from node \(i\) to node \(j\) (see Fig. 5.9 or Fig. 5.13b). When we apply the relation (5.54) to all edges incident to node \(i\), we obtain the following over-constrained system of linear equations

with \(\Delta(\cdot)_{ij} = (\cdot)_j - (\cdot)_i\) and \(\partial_m(\cdot) = \partial(\cdot)/\partial m\). Further, \(N_A\) denotes the number of adjacent nodes \(j\) connected to \(i\) by an edge and \(\theta_j\) stands for some weighting coefficient. The weights can depend on the geometry and/or on the solution (see, e.g., [40]). However, in practice \(\theta_j\) is usually set to unity. For convenience, we abbreviate the above system (5.55) as

Solving Eq. (5.56) for the gradient vector \(\vec{x}\) requires the inversion of the matrix \(\bar{A}\). To prevent problems with ill-conditioning (particularly on stretched grids), Anderson and Bonhaus suggested to decompose \(\bar{A}\) into the product of an orthogonal matrix \(\bar{Q}\) and an upper triangular matrix \(\bar{R}\) using the Gram-Schmidt process [66]. Their approach was recently extended to 3D in [49]. Hence, the solution to Eq. (5.56) immediately follows from

Using a lower case letter with double subscripts to denote a matrix element, we may write the Gram-Schmidt orthogonalisation of the matrix \(\bar{A} = [\vec{a}_1,\, \vec{a}_2,\, \vec{a}_3]\) as \(\bar{Q} = [\vec{q}_1,\, \vec{q}_2,\, \vec{q}_3]\), where

The entries in the upper triangular matrix \(\bar{R}\) are obtained from

Using Eqs. (5.57)-(5.59), the gradient at node \(i\) follows from the weighted sum of the edge differences

with the vector of weights \(\vec{w}_{ij}\) defined as

The terms in the above Equation (5.61) are given by

where

The formulation of the least-squares approach remains formally the same for a cell-centred scheme, only the nodes have to be substituted by cell-centroids. An example may be found in Ref. [16].

The least-squares method is first-order accurate on general grids [63]. It is also consistent, i.e., the gradient of a linear function is computed to roundoff error, regardless of the type of the elements. Therefore, the method is particularly suited to mixed grids. The computational costs are comparable to those of the Green-Gauss approach, since only a vector-scalar multiplication (Eq. (5.60)) is needed within a single loop over faces/edges. However, we have to pre-compute and store the six entries (Eq. (5.59)) of the upper triangular matrix \(\bar{R}\) at each node. A detailed investigation in Ref. [67] revealed that the weighting coefficient \(\theta_j\) in Eqs. (5.55), (5.60) has to be set equal to the inverse of the distance between the nodes \(i\) and \(j\) (similar to \(\theta_{ij}\) in Eq. (5.45)), in order to obtain an accurate approximation for the gradient of a non-linear function on a highly stretched and additionally curved grid. However, this does not help in the case of the cell-centred scheme on triangular (tetrahedral) grids [67].

Experience also shows that the least-squares approach requires some attention in the case of the median-dual scheme, if prismatic or hexahedral cells are employed on a viscous wall. Consider Fig. 5.15 and assume that we want to compute the gradients of the velocity components at node \(i\). It may become obvious that only the contribution from the edge \(ij\) is useful, since at other nodes connected to \(i\) by an edge \(u = v = w = 0\). To increase the support of the stencil, we can insert the so-called virtual edges [49], as they are rendered by the dashed lines in Fig. 5.15. The virtual edges help to improve the accuracy and the robustness of the discretisation scheme considerably. It should be stressed that they are employed only for the gradient reconstruction but not for the flux computation.

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 want to achieve 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. We discussed this point in Subsection 4.3.5 for the case of structured upwind discretisation schemes.

On unstructured grids, the purpose of a limiter is to reduce the gradients used to reconstruct the left and right state at the face of the control volume. The limiter function must be zero at strong discontinuities, in order to obtain a first-order upwind scheme which guarantees monotonicity. Setting the limiter to zero leads to the constant reconstruction of Eq. (5.38). Of course, the original unlimited reconstruction has to be retained in smooth flow regions, in order to keep the amount of numerical dissipation as low as possible. In the following, we shall describe two widely used limiter functions - namely the limiters of Barth and Jespersen [30], and of Venkatakrishnan [68], [69].

The first implementation of a limiter function on unstructured grids was presented in Ref. [30]. In the case of the median-dual scheme, it is defined at the node \(i\) as

with the abbreviations

In Equations (5.64) and (5.65), \(\min_j\) or \(\max_j\) means the minimum or maximum value of all direct neighbours \(j\) of node \(i\) (i.e., all nodes connected to \(i\) by an edge). Furthermore, the edge vector \(\vec{r}_{ij}\), which is shown in Fig. 5.9 or in Fig. 5.13b, is defined according to Eq. (5.43). Finally, \(U_j\) denotes a scalar quantity at some neighbouring node \(j\). Similar formulae to those above hold for the cell-centred scheme with cell instead of node indices and with

where \(\vec{r}_L\) denotes the vector from the cell-centroid to the midpoint of the corresponding cell face. In order to avoid division by a very small value of \(\Delta_2\) in Eq. (5.64), it is better to modify \(\Delta_2\) as \(\mathrm{Sign}(\Delta_2)(\left|\Delta_2\right| + \omega)\), where \(\omega\) is approximately the machine accuracy [68].

Barth's limiter enforces a monotone solution. However, it is rather dissipative and it tends to smear discontinuities. A further problem presents the activation of the limiter due to numerical noise in smooth flow regions. This usually prevents the full convergence to steady state [68], [38]. Therefore, the limiter function due to Venkatakrishnan became more popular.

Venkatakrishnan's limiter [68], [69] is widely used because of its superior convergence properties. The limiter reduces the reconstructed gradient \(\nabla U\) at the vertex \(i\) by the factor

where

In the above Eq. (5.68), \(U_{max}\) and \(U_{min}\) stand for the minimum/maximum values of all surrounding nodes \(j\) and including the node \(i\) itself. Definitions of \(U_{max}\), \(U_{min}\) and \(\Delta_2\) are given in Eq. (5.65). The parameter \(\epsilon^{2}\) is intended to control the amount of limiting. Setting \(\epsilon^{2}\) to zero results in full limiting, but this may stall the convergence. Contrary to that, if \(\epsilon^{2}\) is set to a large value, the limiter function will return a value of about unity. Hence, there will be no limiting at all and wiggles could occur in the solution. In practice, it was found that \(\epsilon^{2}\) should be proportional to a local length scale, i.e.,

where \(K\) is a constant of \(\mathcal{O}(1)\) and \(\Delta h\) is for example the cube-root of the volume (square-root of the area in 2D) of the control volume. It is important to notice that the limiter function (5.67) must be defined with non-dimensional quantities. The influence of the coefficient \(K\) in Eq. (5.69) on the resolution of a shock is demonstrated in Fig. 5.16. It can be seen that the fully limited \((K = 0)\) and the solution for \(K = 5\) are identical. However, the explicit time-stepping scheme converged only about three orders of magnitude for \(K = 0\), whereas for \(K = 5\) it converged to machine zero (Fig. 5.17). Figure 5.16 also shows that the solution becomes gradually unlimited with increasing values of \(K\). This manifests itself as an increasing overshoot at the shock.

The computational effort for the evaluation of one of the above limiter functions is relatively high. Two loops over edges (faces in the case of the cell-centred scheme) and one loop over nodes (cells) are necessary in order to compute \(U_{max}\), \(U_{min}\) as well as the limiter \(\Psi\) itself. Furthermore, \(U_{max}\), \(U_{min}\) and \(\Psi\) have to be stored node-(cell-)wise, separately for each flow variable.