6.2 Implicit Time-Stepping Schemes 隐式时间推进格式[cfd-6-2]

Referring to Eq. (6.27), we can write the linearisation of the residual \(\vec{R}^{n+1}\) in the form

with \(N_F\) being the number of faces of the control volume \(\Omega\) (cf. Eq. (4.2) or (5.2)). Thus, the flux Jacobian reads

It should be stressed that the flux Jacobian has to be conceived as an operator which acts on the update \(\Delta\vec{W}\). As stated above, the convective and viscous fluxes in Eq. (6.30) do not necessarily have to be identical to the fluxes in the explicit operator.

Because of the significant differences between the implicit operators on structured and unstructured grids, we shall treat each case separately.

In order to derive the form of the system matrix in Eq. (6.28), let us consider the 1-D grid in Fig. 6.1. Let us further assume that the cell-vertex scheme with dual control volumes (Subsection 4.2.3) and a simple average of fluxes are used for the spatial discretisation (cf. Fig. 4.8). In the absence of viscous fluxes, the residual is given by

Furthermore, the derivative of the convective fluxes at the face \(m = i + 1/2\) in Eq. (6.30) can be expressed as follows

where \(\bar{A}_c\) denotes the convective flux Jacobian (Section A.9). Hence, according to Eqs. (6.29), (6.31), and (6.32), the residual at \((t+\Delta t)\) is approximated as

Finally, for \(\beta = 1\) and a lumped mass matrix we can derive from Eq. (6.26) the implicit scheme

where \(\bar{I}\) stands for the identity matrix. It should be noted that the unit normal vector in each matrix \(\bar{A}_c\) is evaluated at the same side of the control volume as the associated area \(\Delta S\). As we can see, the implicit operator involves the same 3-point stencil (nodes \(i-1\), \(i\), and \(i+1\)) as the spatial discretisation.

In order to visualise the system matrix, we denote all terms in the implicit operator associated with the central node \(i\) as \(\mathbf{D}\), i.e.,

with the downwind node \((i+1)\) as \(\mathbf{U}\), i.e.,

and with the upwind node \((i-1)\) as \(\mathbf{L}\), i.e.,

respectively. Writing down Eq. (6.34) for all eight nodes of the grid in Fig. 6.1, we obtain the \(8\times 8\) block-tridiagonal matrix displayed on the right side of Fig. 6.1. Each of the blocks \(\mathbf{L}\), \(\mathbf{D}\), and \(\mathbf{U}\) represents a \(3\times 3\) matrix in 1D (because of the three conservation equations).

The same ideas carry over to multiple dimensions. For example, if the spatial discretisation in 2D would involve the 5-point stencil sketched in Fig. 6.2, we would obtain a block-pentadiagonal matrix. This is shown on the right side of Fig. 6.2. The nodes were ordered such that the \(i\)-index runs faster than the \(j\)-index (corresponds to \(mat(i,j)\) in FORTRAN). It should be noted that the second off-diagonal is at the distance of eight elements from the main diagonal (the total number of nodes in \(i\)-direction). Finally, the system matrix would become block-septadiagonal in 3D, if we would employ the 7-point stencil of Fig. 4.9b for the spatial discretisation.

In summary, we can state that the system matrix always possesses a regular, sparse and banded form for structured grids. It should be further mentioned that the term

in Eq. (6.28), which is always located on the main diagonal, can cause difficulties. Namely, if the time step becomes large, some iterative matrix inversion schemes (e.g., Gauss-Seidel) may fail due to reduced diagonal dominance of the implicit operator.

The appearance of the system matrix changes completely when we proceed from structured to unstructured grids. This can be demonstrated with the aid of a small unstructured grid sketched in Fig. 6.3. We want assume that the spatial discretisation is given by the cell-centred scheme presented in Subsection 5.2.1, which uses flow values only from the nearest neighbours (e.g., like a first-order upwind scheme - see Subsections 5.3.2 and 5.3.3). Thus, for example, the stencil for cell 2 includes the cells 18, 10, and 13. The resulting system matrix is displayed on the right side of Fig. 6.3. It is obvious that since the grid cells (nodes in the case of a median-dual scheme) are in general numbered in an arbitrary order, no regular pattern can be expected for the matrix. Only the main diagonal, which contains at least the expression (6.35) and possibly the derivative of the source term, is always present.

The quasi-random distribution of nonzero elements in the system matrix is undesirable. It slows down the convergence of iterative inversion methods like Gauss-Seidel. Furthermore, preconditioning techniques for Krylov-subspace methods like ILU (Incomplete Lower-Upper) factorisation scheme cannot be used efficiently. Therefore, strategies were developed where the cells (nodes) are renumbered such that the bandwidth of the system matrix is considerably reduced, i.e., the nonzero elements are clustered close to the main diagonal. The best-known renumbering strategy is the Reverse-Cuthill-McKee (RCM) algorithm [27], [28]. Figure 6.4 shows the resulting cell numbering and the system matrix when the RCM algorithm is applied to the example grid of Fig. 6.3. As we can see, the bandwidth of the matrix is significantly reduced and the matrix obtains a more regular structure.

Several other renumbering strategies were developed in order to minimise cache misses or to allow for vectorisation of the numerical scheme. An overview can be found, e.g., in Ref. [29].

Depending on the type of the underlying spatial discretisation scheme, an analytical evaluation of the flux Jacobian \(\partial\vec{R}/\partial\vec{W}\) in Eq. (6.28) may become very complex if not impossible. In order to make the concepts more clear, we shall derive the flux Jacobian for inviscid flows and then discuss the extension to the Navier-Stokes equations.

The flux Jacobian is most easily formulated in the case of the central spatial discretisation. As we already saw for the example in Fig. 6.1, \(\partial\vec{R}/\partial\vec{W}\) consists of the convective flux Jacobians (cf. Eq. (6.34)), which can be derived analytically (see Section A.9). Artificial viscosity is usually included in a simplified form, without the non-linear pressure sensor (Eq. (4.55)). We shall return to this point below in Subsection 6.2.3.

The evaluation of the flux Jacobian becomes more involved when one of the flux-vector splitting schemes (Subsection 4.3.2) is used as the basis for its derivation. Let us, for illustration, consider the scheme due to Steger and Warming [30]. Previous investigations [31] revealed that the Steger-Warming splitting is preferable over, e.g., the Van Leer's flux-vector splitting scheme (Eq. (4.60)) for various upwind discretisations of the implicit operator.

The basic idea of the Steger-Warming flux-vector splitting scheme is to divide the convective fluxes into a positive and a negative part, i.e.,

with the fluxes defined as

In Eq. (6.37), \(\bar{A}_{SW}^{\pm}\) denotes the positive/negative Steger-Warming flux-splitting Jacobian. Furthermore, \(\bar{T}\) represents the matrix of right eigenvectors, \(\bar{T}^{-1}\) the matrix of left eigenvectors, and \(\bar{\Lambda}^{\pm}\) stands for the diagonal matrix of positive/negative eigenvalues, respectively (cf. Section A.11). The eigenvalue matrices are defined as [30]

where \(\bar{\Lambda}_c\) is given by Eq. (A.84). Using the splitting defined in Eq. (6.36), we obtain for the product of the flux Jacobian with the update \(\Delta\vec{W}^n\) in Eq. (6.29)

In the above Eq. (6.39), \(\Delta\vec{W}_{L,m}^n\) and \(\Delta\vec{W}_{R,m}^n\) denote the updates of the left and right state at the face \(m\), respectively. On structured grids, the left and right state can be evaluated by the MUSCL approach (Eq. (4.46)). On unstructured grids, the reconstruction methods discussed in Subsection 5.3.3 can be applied. However, the stencil becomes wider with increasing accuracy, which leads to larger bandwidth of the system matrix. Therefore, and in order to reduce the numerical complexity, only first-order accurate approximation is usually employed in Eq. (6.39). As a compromise, we could reconstruct the left and right state with higher accuracy but retain the stencil of the first-order scheme for the evaluation of the derivatives [32].

To proceed with the discussion on the evaluation of the derivatives \(\partial\vec{F}_c^{\pm}/\partial\vec{W}\) in Eq. (6.39), let us consider, e.g., the positive flux at face \(m\)

which becomes with Eqs. (6.37), (6.38)

An expression similar to Eq. (6.40) can also be found for the negative flux. As we can see, the first term in Eq. (6.40) consists of convective flux Jacobians (see Section A.9) and thus presents no difficulty. However, the second term involves derivatives of matrix elements. Although it is possible to obtain the derivatives analytically either by hand calculation or by using a symbolic algebra package, this will produce a large, computationally inefficient code [33]. Alternatively, it is possible to assume the matrix \(\bar{A}_c\) is locally constant so that the second term in Eq. (6.40) can be neglected. However, depending on the type of the implicit scheme, this may severely restrict the CFL number [34].

Other approaches that we could use to compute the derivatives \(\partial\vec{F}_c^{\pm}/\partial\vec{W}\) in Eq. (6.39) would be the automatic differentiation of the source code (e.g., using ADIFOR [35]) or the finite-difference method (see, e.g., [26], [33]). Herewith, the derivative of the \(i\)-th component of a vector \(\vec{F}\) with respect to the \(j\)-th component of a dependent variable \(\vec{X}\) can be approximated as

where \(\vec{e}^{\,j}\) denotes the \(j\)-th standard basis vector. Dennis and Schnabel [36] suggested a stepsize \(h_j\) of the form

with \(\epsilon\) being the machine accuracy and \(\text{typ}\,x_j\) a typical size of \(x_j\). The reader is also referred to [37] and [38] for hints on efficient numerical evaluation of Jacobian matrices.

In the case of the flux-difference splitting scheme due to Roe (Subsection 4.3.3, Eq. (4.91)), we can write the product of the flux Jacobian with the update in Eq. (6.29) as

Similar to flux-vector splitting, the expression (6.43) contains convective flux Jacobians as well as derivatives of the Roe matrix \(\bar{A}_{Roe}\). The derivatives were presented in [34]. Since the corresponding formulae are very complex (see also Ref. [33]), it is a better idea to evaluate the term \(\partial\vec{F}_c/\partial\vec{W}\) numerically, as discussed above.

However, we can also simplify Eq. (6.43) by assuming locally constant Roe matrices [39]

In contrast to the Steger-Warming flux-vector splitting scheme, the above approximate linearisation (6.44) degrades the performance of the implicit scheme only slightly [34].

For the Navier-Stokes equations, we have to account also for the viscous fluxes in the implicit operator. The derivative \(\partial\vec{F}_v/\partial\vec{W}\), i.e., the viscous flux Jacobian in Eq. (6.30) is in general not straightforward to obtain. Additional complexity arises due to the fact that the viscous flux vector contains derivatives of flow variables. For this reason, we have either to evaluate the viscous flux Jacobian by finite differences (Eq. (6.41)), or we have to use a simplified formulation.

In the case of the TSL approximation of the Navier-Stokes equations (cf. Subsection 2.4.3 and Section A.6), it is possible to find the viscous flux Jacobian analytically by assuming locally constant dynamic viscosity and thermal conductivity coefficients. Then, according to the discussion in Appendix A.10, the term related to the viscous fluxes in Eq. (6.29) becomes

In above Eq. (6.45), \(\bar{A}^{*}_v\) stands for the viscous flux Jacobian given by Eq. (A.71) or Eq. (A.75) but without the spatial operators \(\partial_y(\cdot)\) (cf. Eq. (A.74) and (A.79)). The Jacobians are evaluated using either the left or the right state for all variables except for the dynamic viscosity, which is determined from an arithmetic average. It should be mentioned that Eq. (6.45) leads to a second-order central difference approximation in the implicit operator, supposed the left and right state are computed with first-order accuracy.