5.3.3 Solution Reconstruction 解重构[cfd-5-3-3]

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].