3.2.2 Implicit Schemes 隐式格式[cfd-3-2-2]

A family of implicit time integration schemes is obtained from Eq. (3.4) by setting \(\beta\neq 0\). Very popular for the simulation of unsteady flows is the 3-point implicit backward-difference scheme with \(\beta=1\) and \(\omega=1/2\), which is 2nd-order accurate in time. In this case, the scheme is mostly employed within the so-called dual time-stepping approach [148]-[150], [110], [111], where a steady-state problem is solved in pseudo-time at each physical time step.

For the solution of stationary flow problems, a scheme with \(\omega=0\) is more suitable, since it requires less computer memory. Herewith, if we linearise the residual \(\vec{R}^{n+1}\) in Eq. (3.4) about the current time level, we obtain the scheme

The term \(\partial\vec{R}/\partial\vec{W}\) is denoted as the flux Jacobian. It constitutes a large sparse matrix. The expression enclosed in parenthesis on the left-hand side of Eq. (3.8) is also referred to as the implicit operator. As already discussed above, the mass matrix \(\overline{M}\) can be replaced by the identity matrix, without influencing the steady state solution. The parameter \(\beta\) in Eq. (3.8) is generally set to 1, which results in a 1st-order accurate temporal discretisation. A 2nd-order time accurate scheme is obtained for \(\beta=1/2\). However, this is not advised since the scheme with \(\beta=1\) is much more robust, and the time accuracy plays no role for steady problems anyway.

The principal advantage of implicit schemes as compared to explicit ones is that significantly larger time steps can be used, without hampering the stability of the time integration process. In fact, for \(\Delta t\rightarrow\infty\) the scheme (3.8) transforms into standard Newton's method, which exhibits quadratic convergence. However, the condition for quadratic convergence is that the flux Jacobian contains the complete linearisation of the residual. Another important advantage of implicit schemes is their superior robustness and convergence speed in the case of stiff equation systems and/or source terms, which are often encountered in real gas simulations, turbulence modelling, or in the case of highly stretched grids (high Reynolds number flows). On the other hand, the faster (in terms of time steps or iterations) and the more robust an implicit scheme is, the higher is usually the computational effort per time step or iteration. Therefore, an explicit scheme accelerated by multigrid can be equally or even more efficient. Furthermore, implicit schemes are significantly more difficult to vectorise or to parallelise than their explicit counterparts.

Written down for each control volume, the implicit scheme in Eq. (3.8) represents a large system of linear equations, which has to be solved for the update \(\Delta\vec{W}^{n}\) at each time step \(\Delta t\). This task can be accomplished using either a direct or an iterative method.

The direct methods are based on the exact inversion of the left-hand side of Eq. (3.8) using either the Gaussian elimination or some direct sparse matrix method [151], [152]. Although quadratic convergence was demonstrated on structured [153]-[156] as well as on unstructured grids [157], direct methods are not an option for 3-D problems because they require an excessively high computational effort and a huge amount of computer memory.

Thus, the only practical method for larger grids or 3-D problems are iterative methods. Here, the linear system is solved for \(\Delta\vec{W}^{n}\) at each time step using some iterative matrix inversion methodology. In order to reduce the memory requirements and also to increase the diagonal dominance, the flux Jacobian \(\partial\vec{R}/\partial\vec{W}\) is mostly based on linearisation of a 1st-order accurate spatial discretisation of the right-hand side. The two main consequences of this approximation are that the quadratic convergence of Newton's scheme cannot be achieved and that the maximum time step becomes limited. On the other hand, the numerical effort of an iteration step is significantly reduced, which leads to a numerically highly efficient scheme.

In the case of structured grids, iterative methods like the Alternating Direction Implicit (ADI) scheme [158]-[161], the (line) Jacobi or the Gauss-Seidel relaxation scheme [162]-[166], and particularly the Lower-Upper Symmetric Gauss-Seidel (LU-SGS; also referenced to as LU-SSOR - Lower-Upper Symmetric Successive Overrelaxation) scheme [167]-[171] are mainly employed. All these methods are based on splitting of the implicit operator into a sum or product of parts, which can be each inverted more easily. Because of the associated factorisation error (the difference with respect to the original matrix) and also the simplification of the flux Jacobian, it does not pay off to solve the linear system very accurately. In fact, only one iteration is carried out at each time step of the ADI and the LU-SGS method.

Implicit iterative methods for unstructured grids are in the most cases based on the Gauss-Seidel relaxation scheme [172]-[175]. In order to improve the convergence, it is possible to use the red-black Gauss-Seidel methodology. Its extension to unstructured grids was demonstrated in Refs. [176]-[178]. A particularly interesting possibility is also offered by an implementation of the LU-SGS scheme on unstructured grids [18], [179], [180], because of its very low memory requirements and numerical effort.

Because of the success of the line-implicit methods on structured grids, a few attempts were made to adopt this methodology on unstructured grids [181], [182]. The approach was to construct continuous lines such that each grid point or each grid cell (in the case of a cell-centred scheme) is visited only once - the so-called Hamiltonian tour [183]. The lines were oriented primarily in coordinate directions, but they were folded at the boundaries and where necessary (therefore they were nicknamed "snakes"). A tri-diagonal solver was then employed to invert the left-hand side of Eq. (3.8). Later on, it was recognised that folding the lines can slow down the convergence. To overcome this, each line was broken up into multiple linelets [184]. However, the performance on a vector computer was rather poor. The idea of linelets was also employed to improve the convergence of an explicit scheme on highly stretched viscous unstructured grids using an implicit solver in the direction across the boundary layer [145].

More sophisticated iterative techniques, which treat the linear equation system in a more global way, are the so-called Krylov subspace methods. Their development was triggered by the introduction of an efficient iterative scheme for solving large, sparse linear systems - namely the conjugate gradient method by Hestenes and Stiefel [185]. The original conjugate gradient method is restricted to Hermitian positive definite matrices only, but for an \(n\times n\) matrix it converges in at most \(n\) iterations. Since then, a variety of Krylov subspace methods was proposed for the solution of arbitrary non-singular matrices, as they occur in CFD applications. For example, there are methods like the Conjugate Gradient Squared (CGS) [186], the Bi-Conjugate Gradient Stabilised (Bi-CGSTAB) [187], or the Transpose-Free Quasi-Minimum Residual (TFQMR) [188] scheme.

However, the most widely employed method is the Generalised Minimal Residual (GMRES) scheme developed by Saad and Schultz [189]. If we rewrite the implicit scheme in Eq. (3.8) as

then \(\bar{J}\) represents a large, sparse, and non-symmetric matrix (the left-hand side). Starting from an initial guess \(\Delta\vec{W}_0^{n}\), the GMRES(\(m\)) method seeks a solution \(\Delta\vec{W}^{n}\) in the form \(\Delta\vec{W}^{n}=\Delta\vec{W}_0^{n}+\vec{y}_m\), where \(\vec{y}_m\) belongs to the Krylov subspace

such that the residual \(\Vert\bar{J}\,\Delta\vec{W}^{n}+\vec{R}^{n}\Vert\) becomes a minimum. The parameter \(m\) specifies the dimension of the Krylov subspace, or in other words the number of search directions (\(\bar{J}^{i}\vec{r}_0\)). Since all directions have to be stored, \(m\) is usually chosen between 10 and 40, the higher number being necessary for poorly conditioned matrices (which arise in the simulation of turbulent flows, real gas, etc.). GMRES has to be restarted, if no convergence is achieved within \(m\) sub-iterations. The GMRES method requires significantly more memory than, e.g., Bi-CGSTAB or TFQMR, but it is more robust, smoothly converging and usually also faster. A very detailed comparison of the various methodologies can be found in Ref. [190].

Nevertheless, as with other conjugate gradient methods, preconditioning is absolutely essential for CFD problems. Here, we solve

instead of the system in Eq. (3.9). The matrices \(\bar{P}_L\) and \(\bar{P}_R\) denote left and right preconditioners, respectively. The preconditioner should approximate \(\bar{J}^{-1}\) as close as possible, in order to cluster the eigenvalues near unity. On the other hand, it should be of course easy to invert. One particularly efficient preconditioner is the Incomplete Lower Upper factorisation method [191] with zero fill-in (ILU(0)). For the discussion of different preconditioning techniques in connection with GMRES the reader is referred to [97], [192]-[195].

Since the GMRES method requires a considerable amount of computer memory for storing the search directions and possibly also the preconditioning matrix, it is a good idea to circumvent an explicit formation and storage of the flux Jacobian \(\partial\vec{R}/\partial\vec{W}\). This is offered by the so-called matrix-free approach. The idea is based on the observation that GMRES (and some other Krylov subspace methods) employs only matrix vector products of the form

which can be simply approximated by finite-differences as

thus requiring only residual evaluations. The parameter \(\epsilon\) has to be chosen with some care, in order to minimise the numerical error (see, e.g., [196] or [197]). Another, and even more important, advantage of the matrix-free approach is that (numerically) accurate linearisation of a high-order residual \(\vec{R}^{n}\) can be easily utilised in the implicit scheme. Hence, the quadratic convergence of Newton's scheme can be achieved at moderate costs. In this case we speak of Newton-Krylov approach [197]-[201], [110]. Practical experience indicates that from all Krylov subspace methods, GMRES is best suited for the matrix-free implementation [202]. An interesting possibility is to utilise the LU-SGS scheme as a preconditioner for the matrix-free GMRES method. Since the LU-SGS scheme also does not require an explicit storage of the flux Jacobian, the memory requirements can be even further reduced. The computational efficiency of this approach was recently demonstrated for 3-D inviscid and laminar flows on unstructured grids [203].

The convergence of an implicit scheme can also be enhanced by using multigrid. There are basically two possible ways. First, we can employ multigrid inside an implicit scheme - as a solver for the linear equation system (3.9) arising at each time step, or as a preconditioner for one of the conjugate gradient methods [204], [205]. Second, the implicit scheme itself can serve as a smoother within the FAS multigrid method, which is applied directly to the governing equations [206]-[209], [178]. Some investigations show that at least for purely aerodynamic problems, rather "simple" implicit schemes (like Gauss-Seidel) combined with multigrid result in computationally more efficient solvers (in terms of the CPU time) than, e.g., GMRES [178], [198].