chapter. Unstructured Finite Volume Schemes 第5章 非结构网格有限体积格式[cfd-0006]

Generally, we think of the flow in a plane as being 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 2D from the product of its area with the depth \(b\). Since the depth \(b\) is arbitrary, we may set \(b = 1\) for convenience. In the following discussion, we restrict ourselves to triangular and quadrilateral elements. Even though the control volume of a median-dual scheme can have a rather complex shape, it can always be decomposed into triangles and/or quadrilaterals.

The area of a general triangle can be most conveniently and exactly calculated by the formula of Gauss. Thus, using a node numbering in accordance with Fig. 5.3a, the volume results from

The nodes have to be numbered in the anti-clockwise direction in order to obtain a positive value for the volume.

The area of a general quadrilateral can be exactly calculated by Gauss' formula, which leads, after some algebra, to the expression

where the nodes are numbered according to Fig. 5.3b in the anti-clockwise direction. In the above, we assumed that the control volume is located in the \(x-y\)-plane and that the \(z\)-coordinate represents the symmetry axis.

The edges of a control volume are given by straight lines in 2D and therefore the unit normal vector is constant along them. When we integrate the fluxes according to the approximation of Eq. (5.2), we have to evaluate the product of the area of a face \(\Delta S\) and the corresponding unit normal vector \(\vec{n}\) which is the face vector \(\vec{S}\). Considering Fig. 5.3, the outward pointing face vector, e.g., at the side 2-3 is given by

Because of the symmetry, the \(z\)-component of the face vectors (and of the unit normal vector) is zero. It is therefore omitted in Eq. (5.6). The unit normal vector can be obtained from Eq. (5.6) with

where \(S_x\), \(S_y\) denote the Cartesian components of the face vector.

The centre of the triangle from Fig. 5.3a is defined as

with \(\vec{r}_{1/2/3}\) representing the Cartesian coordinates of the nodes. The centre of a quadrilateral element can be computed by the formula given in Ref. [34]. For this purpose, the quadrilateral is decomposed into two triangles which share two points. With the node numbering according to Fig. 5.3b, and 1 and 3 being the common nodes, the relation reads

The volumes of the two triangles \(\Omega_{123}\) and \(\Omega_{134}\) are evaluated using Eq. (5.4), and their centroids \(\vec{r}_c\) are obtained from Eq. (5.8).

As opposed to the previous 2-D case, the computation of face vectors and volumes poses in 3D some problems for elements or control volumes with quadrilateral faces. The main reason for this is that, in general, the four vertices of a quadrilateral face of a control volume may not lie in a plane. Then, the normal vector is no longer constant on such face (see Fig. 4.2). In order to overcome this difficulty, we could decompose each quadrilateral face into two or even more triangles. However, the gain in accuracy is hardly noticeable for a second-order scheme on a smooth grid. The additional effort can only be justified - and in fact it becomes necessary - for a third- and higher order spatial discretisations. Therefore, we shall apply a simplified treatment of the quadrilateral faces in the following considerations, which is based on an averaged normal vector.

The face vector \(\vec{S}\) can be exactly computed for a triangular face using Gauss' formula. Defining the nodes according to Fig. 5.4a, we obtain for the edge differences of the triangle 1-2-3

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

The averaged face vector \(\vec{S}\) of a quadrilateral face, like that rendered in Fig. 5.4b, is most conveniently computed using the same Gauss' formula as employed in 2-D for the area of a quadrilateral. Thus, for the face given by the nodes 5, 6, 7 and 8 in Fig. 5.4b, we first define the differences

Then, we obtain the outward pointing face vector \(\vec{S} = \vec{n}\Delta S\) from the relation

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 in both cases from \(\vec{n} = \vec{S}/\Delta S\) with

where \(S_x\), \(S_y\) and \(S_z\) denote the Cartesian components of the face vector given by Eq. (5.11) or Eq. (5.13), respectively.

As we already stated in the case of 3-D structured finite volume schemes, a very convenient approach for the computation of volumes is based on the divergence theorem [35]. The discussion in Subsection 4.1.2 led finally to the expression

for the volume, where \(N_F\) denotes the number of the faces of the control volume, \((\vec{r}_c)_m\) the centre of the face \(m\) of the control volume, and \(\vec{S}_m\) the face vector (outward directed) of the face \(m\), respectively. The formula (5.15) is directly applicable on unstructured grids. It is exact for a volume with triangular faces, or a volume with planar quadrilateral faces.

The previously mentioned median-dual type of control volume requires the knowledge of the centroid of the grid cell. The centroid of a general volume is defined as

According to the derivation in Ref. [34], the relation in Eq. (5.16) can be discretised as

where the centre of a face \(m\), i.e., \((\vec{r}_c)_m\) is obtained from Eq. (5.8) for a triangular face, or from Eq. (5.9) for a quadrilateral face. As we can note, the denominator in Eq. (5.17) is identical to \(\Omega/3\) from Eq. (5.15).

We speak of a cell-centred scheme if the control volumes are identical with the grid cells and if the flow variables are associated with their centroids, as it is sketched in Fig. 5.6. When we evaluate the discretised flow equations (5.2), we have to supply the convective and the viscous fluxes at the midpoints of the faces of the control volume, which is sufficient for a second-order accurate discretisation on smooth grids (averaged normal vector is employed for quadrilateral faces). The fluxes can be approximated in one of three 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 unit normal vector (generally applied only to the convective fluxes);
• by using an average of variables associated with the centroids of the grid cells adjacent to the left and to the right side of the cell face;
• by computing the fluxes from flow quantities reconstructed separately on both sides of the cell face from values in the surrounding cells (employed only for the convective fluxes).

Thus, considering for example the cell face with the unit normal vector \(\vec{n}_{01}\) in Fig. 5.6, the first approach - average of fluxes - reads in two dimensions

with the face area \(\Delta S_{01}\) computed from Eqs. (5.6) and (5.7).

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

where the conservative/dependent variables at the face with the unit normal vector \(\vec{n}_{01}\) are defined as the arithmetic average of values at the two adjacent cells. Thus,

The flux vector \(\vec{F}\) in Eq. (5.19) represents either the convective or the viscous fluxes.

The third methodology starts with an interpolation of flow quantities (usually velocity components, pressure, density and total enthalpy) separately to both sides of the cell face. The reconstructed quantities - termed the left and the right state (see Subsection 5.3.3) - differ in general. The fluxes through the cell face are then evaluated from the difference of the left and right state using an appropriate non-linear function. Hence,

where

represent the reconstructed states.

Of course, similar relations hold for the other control volume faces as well. The above approximations can be employed in the same way in three dimensions. The face vector \(\vec{S}\) is then evaluated using the formulae (5.11) or (5.13), respectively.

As we already stated in the introduction to this section, the basic data structure which describes the elements has to be extended in an appropriate way to support the discretisation methodology. It is obvious from the previous discussion that numerical operations are carried out using mainly the faces of the elements (control volumes) together with values at the centres of the adjacent cells. It is therefore quite natural to employ a face-based data structure for the spatial discretisation. Such data structure stores for each particular face in the grid (see Fig. 5.6):

• pointers to the two cells which share the respective face - this allows us to access the flow variables associated with the two cells (\(C_0\), \(C_1\));
• the face vector (\(\vec{S}_{01} = \vec{n}_{01}\Delta S_{01}\)) - must point consistently either outwards or inwards;
• two vectors from the centroid each cell to the midpoint of the face \(M_{01}\), i.e., (\(C_0-M_{01}\)), (\(C_1-M_{01}\)) - required for an accurate interpolation of flow variables to the face. This is not necessary for purely tetrahedral grids, where a simple extrapolation formula can be used [26], [2] (cf. Eq. (5.44)).

Hence, the integration of the fluxes (e.g., according to Eq. (5.19)) would be implemented as a loop over all (i.e., internal and boundary) faces contained in the grid:

After the loop is completed and the source term \(\vec{Q}_I\Omega_I\) is added, we obtain the final residuals (\(\vec{R}\)) in all cells. A less efficient approach would be to loop over elements because the face vectors would have to be stored twice and the fluxes would be computed twice (with the exception of the boundaries). Furthermore, because we use exactly the same face vector \(\vec{S}_{IJ}\) in order to to evaluate the partial fluxes into the volumes \(\Omega_I\) and \(\Omega_J\), the conservation properties of the governing equations are automatically retained.

Within the cell-vertex scheme, the flow variables are associated with the grid nodes (vertices). Median-dual control volumes are formed by connecting the centroids, face- and edge-midpoints of all cells sharing the particular node. This is depicted in Fig. 5.7a for a tetrahedron and in Fig. 5.7b for a hexahedron. The definition of a median-dual control volume results in a polyhedral hull around each grid node, as it is sketched in Fig. 5.8 for a 2-D mixed grid. This polyhedra can be viewed as a dual grid - hence the name of the scheme. It is interesting to note that the median-dual finite volume discretisation is equivalent to the Galerkin finite element scheme with linear elements (see, e.g., [36]).

In order to evaluate the discretised flow equations (5.2), we have to integrate the convective and viscous fluxes over the surface of the control volume. Hence, we would have to compute the fluxes for each partial face (e.g., \(F_1-M_{13}-F_2-C\) in Fig. 5.7a) separately. However, this is only required for a third- or higher-order accurate discretisations [37], [38]. In the case of a second-order scheme, which is most frequently employed, we may assume the flow variables to be constant for all faces grouped around a particular edge. The fluxes are then evaluated at the midpoint of the edge using the variables and the gradients from both nodes. This approach allows us to define a mean unit normal vector and a total face area associated with each edge. Thus referring to Fig. 5.8, the mean normal vector, e.g., for the edge \(P_0-P_1\) becomes

and the total face area is given by: \(\Delta S_{01} = \Delta S_L + \Delta S_R\). The same applies also in 3D, where the mean normal vector results from a sum over all partial faces having the particular edge-midpoint in common, as it is rendered in Fig. 5.9. The face vector (\(\vec{S} = \vec{n}\Delta S\)) is computed in 2D from Eq. (5.6). In three dimensions, where the partial faces are always quadrilaterals, we can either divide them into triangles and use Eq. (5.11), or we can employ a simplified treatment due to Eq. (5.13), which is sufficient for smooth grids. The element centroids and the face centres are obtained by the formulae (5.17), (5.8), or (5.9), respectively.

The fluxes can then be evaluated according to one of the three following methodologies:

• by the average of fluxes computed from values at both nodes of an edge, but using the same mean unit normal vector (generally applied only to the convective fluxes);
• by using an average of variables stored at the two nodes of an edge;
• by computing the fluxes from flow quantities reconstructed separately on both sides of the face of the control volume from values at the surrounding nodes (employed only for the convective fluxes).

The computation of fluxes follows formally the same approaches as for the cell-centred scheme. Thus, the formulae (5.18)-(5.22) are applicable also in the case of the median-dual scheme. If we utilise the above approach which associates each edge with a mean unit normal, the most efficient methodology is to employ an edge-based data structure for the spatial discretisation. The edge-based data structure stores for each particular edge in the grid (cf. Fig. 5.9):

• pointers to the two nodes which define the edge - this allows us to access the flow variables associated with the two control volumes \(\Omega_i\) and \(\Omega_j\);
• the face vector (\(\vec{S}_{ij} = \vec{n}_{ij}\Delta S_{ij}\)) - must point consistently either outwards or inwards;
• the edge vector from node \(i\) to node \(j\) - required for the interpolation of flow variables to the face (solution reconstruction). Alternatively, the edge vector can be computed on the fly from coordinates of the nodes. This is not required for the standard central scheme with artificial dissipation (Subsection 5.3.1).

With this, the integration of the fluxes (e.g., according to Eq. (5.19)) would be implemented as a loop over all edges in the grid:

After the loop is completed and the source term \(\vec{Q}_i\Omega_i\) is added, we obtain the final residuals (\(\vec{R}\)) in all nodes. This approach is significantly more efficient than summing up the fluxes over each control volume separately, because we store each mean face vector only once and we also visit each edge only once instead of twice. Furthermore, since we use exactly the same mean face vector \(\vec{S}_{ij}\) in order to to evaluate the partial fluxes into the volumes \(\Omega_i\) and \(\Omega_j\), the mass, momentum and energy remain exactly conserved.

The relative advantages and disadvantages of the cell-centred and the median-dual scheme are the subject of controversial debates. The main reason is the lack of fair comparisons of the two methodologies with respect to accuracy, computational time and memory for realistic configurations. Our intention here is to collect the most important arguments for and against each of the approaches regarding:

• accuracy,
• computational work,
• memory requirements, and
• flexibility.

This should lead to a greater understanding of the problems inherent to each scheme and should be of help in selecting the most suitable scheme for the intended applications.

A cell-centred scheme on a triangular/tetrahedral grid leads to about twice/six times as many control volumes and hence degrees of freedom as a median-dual scheme [36]. On typical mixed grids, which consist of tetrahedra and prisms, a cell-centred scheme gives roughly three times more unknowns than a median-dual scheme. This suggests that cell-centred schemes are more accurate than cell-vertex discretisations on an identical grid. However, the residual of a cell-centred scheme results from a much smaller number of fluxes as compared to a median-dual scheme (three versus approximately seven on a tetrahedral grid), which may impair the accuracy. Thus, there is no clear evidence about which scheme might be superior.

The median-dual scheme suffers from a particular problem on stretched triangular and tetrahedral grids. Consider, for example, Fig. 5.10, which shows a tessellation composed of right triangles, as it is often employed near solid walls for viscous flows. We can see in Fig. 5.10a that the face \(\Delta S_{ij}\) becomes highly skewed with respect to the edge \(ij\). However, spatial discretisation schemes mostly assume fluxes to be orthogonal to a face (especially Riemann solvers). Thus, an error is introduced which is particularly significant for a first-order scheme [38]. The situation can be improved using the so-called containment-dual control volume [39]. As depicted in Fig. 5.11, the containment-dual approach employs the centres of the minimum spanning circles/spheres instead of the cell centroids to define the faces. This leads to control volumes identical to those on quadrilateral grids (Fig. 5.10b). Notice that there is no face area associated with diagonal edges like \(ij'\). An additional effort is required for pre-processing, but the solution accuracy can be improved noticeably [40]. Of course, another possibility is to employ directly quadrilateral or hexahedral

cells within the boundary layers. Further discussion of grid-induced errors can be found in Ref. [22] and [23].

Another problem inherent to the median-dual scheme is the discretisation at boundaries of the physical domain. What happens is that there is only about one half of the control volume left at the boundary (cf. Fig. 4.6). The integration of fluxes around the faces results in a residual located inside - ideally at the centre - of the control volume. However, the residual is associated with the node, residing directly on the boundary. This mismatch leads to increased discretisation error in comparison to the cell-centred scheme, which is particularly undesirable on solid walls. The definition of the dual control volume causes also problems at sharp corners (like trailing edges), which show up as unphysical peaks in pressure or density. Further complications arise at periodic boundaries (see Chapter 8.8), where the fluxes from both parts of the control volume have to be summed up correctly.

The mismatch between the centroid of the control volume and the node where the residual is stored has also a further implication for the median-dual scheme. It arises as the mass matrix in the case of unsteady flows. We discussed this point already at the beginning of Section 3.2. The advantage of the cell-centred scheme is that the mass matrix can be eliminated from the equations, without compromising the solution accuracy. By contrast, the median-dual scheme requires a special treatment of the mass matrix [41], [42].

In order to judge the computational effort required for both schemes, we have to consider primarily the integration of the fluxes. We know from the previous discussion that the cell-centred scheme uses a loop over cell faces whereas the median-dual scheme loops over the edges. Since the evaluation of the fluxes at an interface is quite similar for both schemes, the ratio of the number of cell faces to the number of edges gives the ratio of the computational work. Thus, on a tetrahedral grid, where the number cell faces (if counted only once for every two cells) is approximately two times larger than the number of edges, the cell-centred scheme is computationally twice as much expensive as the median-dual scheme on an identical grid [36]. The cell-centred approach becomes however more competitive on mixed grids containing prismatic elements. Apart from boundary treatment, both methods are computationally equivalent on hexahedral grids, where the number of faces equals the number of edges.

Considering the memory requirements, the cell-centred scheme has to store about six times more flow variables on tetrahedral and about three times more variables on usual mixed grids as compared to the median-dual scheme. Furthermore, as we saw, both schemes require to store two integers and three reals (pointers and face vector) per cell face or edge, respectively. Additionally, the cell-centred scheme has to keep two vectors to the face-midpoint - 6 reals - per cell face in memory. On the contrary, the median-dual scheme can work with the node coordinates only, which are considerably fewer values. Thus in summary, the cell-centred scheme needs, on average, more than twice as much computer memory as the median-dual method.

One significant advantage of the cell-centred scheme appears in the case of non-conforming cell interfaces, like those at the letter "F" in Fig. 3.4. In contrast to the median-dual methodology, no special and expensive procedure is required for the computation of the fluxes at the interface. This allows for an increased flexibility in the grid generation and also in the grid adaptation.

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.

A common feature of this type of gradient computation is the necessity to store either information about the grid elements or some coefficients related to the geometry of the elements. Hence, we have to extend the data structure beyond the face-/edge-based formulation presented earlier for the convective fluxes. Below, we discuss three well-established methods for the cell-centred and the median-dual discretisation.

One possible way of evaluating the gradients at a face of the control volume is to define an auxiliary control volume centred at the face and to employ the Green-Gauss theorem. We already discussed this approach in Section 4.4 in the framework of the structured finite volume discretisation. For example, in the case of the median-dual scheme we can compute the gradient at the edge-midpoint as the volume average of gradients for all elements which share the edge [70]. The element-based gradients are evaluated according to Eq. (5.49) by looping over all grid cells and accumulating the gradients at the edges. The values of \(U\) at the cell-faces are obtained by averaging the nodal values, in a manner similar to Eq. (5.53). This approach is relatively costly in terms of the memory and the number of operations. However, it can be implemented for any mix of grid elements.

Another methodology, which is applicable to the median-dual scheme, was derived from the Galerkin finite element method [31]. Basically speaking, the approach transforms the integration of gradients over the surface of the control volume into an evaluation of the Hessian matrix (second derivatives) at the central node. The viscous terms then follow the differential form of the Navier-Stokes equations in Cartesian coordinates (Eq. (A.4) with \(\xi = x\), \(\eta = y\), \(\zeta = z\) and \(J^{-1} = 1\)), which contains terms such as

with \(\partial_m(\cdot) = \partial(\cdot)/\partial m\). Hence, no further integration of the viscous fluxes over the faces of the control volume is required.

The original scheme was formulated for purely triangular/tetrahedral grids. It employs a union of all elements that contain the particular node. In order to simplify the implementation, the dynamic viscosity coefficient is averaged from the nodal values, which is a difference to the Galerkin method. Then, the second derivatives can be evaluated at node \(i\) as follows [71]

The volume \(\Omega'\) contains all tetrahedra which share the node \(i\). The coefficient matrix \(\alpha\) in Eq. (5.71) is symmetric about the diagonal [71], i.e., \(\alpha_{xy} = \alpha_{yx}\), \(\alpha_{xz} = \alpha_{zx}\), and \(\alpha_{yz} = \alpha_{zy}\), respectively. Thus, it is necessary to store only six coefficients for each edge. The coefficients are given by [71]

where \(n\), \(k\) denote the \(x\), \(y\), \(z\) subscripts, and \((\vec{S}_e^{i})_k\), \((\vec{S}_e^{j})_n\) represent components of the outer face vectors \(\vec{S}_e^{i}\), \(\vec{S}_e^{j}\) displayed in Fig. 5.18. The summation is carried out over all tetrahedra (with particular volumes \(\Omega_e\)) which share the edge \(ij\). If the grid is stationary, the coefficients can be computed in a pre-processing step. A desirable feature is that the same edge-based data structure can be employed as for the convective fluxes.

The disadvantage of this approach is however that the full viscous terms are retained only on triangular or tetrahedral grids. For elements like prisms or hexahedra, this technique simplifies to a TSL-like approximation of the Navier-Stokes equations in all three coordinate directions [8]. An extension to non-simplex elements which conserves the full viscous terms was presented in [72]. But the efficient edge-based data structure can then no longer be used since the extended stencil involves nodes not connected by an edge to point \(i\).

This scheme is intended for the cell-centred type of control volume and purely tetrahedral grids. It employs a modified version [4] of the stencil introduced in [61] to evaluate the gradients at the cell faces. The approach is based on an average of the values at the three nodes which define the cell face, combined with the quantities at the cell centroids. The first derivatives at a cell face result from the solution of the linear system of equations [4]

Referring to Fig. 5.19, the subscripts \(I\), \(J\) denote the cell-centroids, and the subscripts \(1\), \(2\), \(3\) stand for the nodes \(P_1\), \(P_2\) and \(P_3\), respectively. Flow variables at the nodes can be determined either from the inverse-distance weighting (Eq. (5.45)) or by the pseudo-Laplacian weighting [2] (similar to Eq. (5.24)).

Since we already computed the gradients inside each control volume (e.g., using the piecewise linear reconstruction, Eq. (5.41) or (5.42)), it would be tempting to evaluate the gradient at the face-midpoint by the simple average [73]

This approach is particularly attractive, because it requires only the basic face- or edge-based data structure and no additional storage. However, as it was pointed out, e.g., in Ref. [71], it leads to a wide stencil with an unfavourable distribution of the weights [49]. Furthermore, it was demonstrated in [49] that the stencil allows for the decoupling of the solution on quadrilateral or hexahedral grids.

The properties of the method can be improved, and particularly the decoupling can be prevented, by using the directional derivative along the connection between the cell-centroids (in the case of the cell-centred scheme), i.e.,

where \(\ell_{IJ}\) represents the distance between the both cell-centroids \(I\) and \(J\) (dashed line in Fig. 5.19). A similar expression holds also for the median-dual scheme with \(\vec{r}_{ij}\) according to Eq. (5.43). With the definition of the unit vector \(\vec{t}_{IJ}\) along the line connecting \(I\) and \(J\),

the modified average may be written as [74], [75]

where \(\overline{\nabla U}_{IJ}\) is given by Eq. (5.74). The modification leads to strongly coupled stencils on tetrahedral as well as on prismatic or hexahedral grids [49]. The modified approach is also still compatible with the face-/edge-based data structure and requires no additional storage. It is therefore more attractive than the element-based methodology, provided the gradients inside control volumes are utilised for the convective fluxes anyway.