6.2.1 Matrix Form of the Implicit Operator 隐式算子的矩阵形式[cfd-6-2-1]

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