5.3.4 Evaluation of the Gradients 梯度的计算[cfd-5-3-4]

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.