chapter. Temporal Discretisation 第6章 时间离散[cfd-0007]

An explicit scheme starts from a known solution \(\vec{W}^n\) and employs the corresponding residual \(\vec{R}^n\) in order to obtain a new solution at time \((t+\Delta t)\). In other words, the new solution \(\vec{W}^{n+1}\) depends solely on values already known. This fact makes the explicit schemes very simple and easy to implement.

As we discussed it in Subsection 3.2.1, a basic explicit scheme can be derived from Eq. (6.2) by setting \(\beta = 0\) and \(\omega = 0\). This results in

which is termed the forward Euler approximation. The mass matrix \(\bar{M}\) can be lumped (i.e., substituted by the identity matrix) for steady problems or for the cell-centred discretisation.

The by far most popular and widespread explicit method is the multistage (Runge-Kutta) time-stepping scheme and its variant the hybrid multistage scheme. Therefore, we shall describe both methods below.

Various implicit time integration schemes can be obtained by setting \(\beta \ne 0\) in Eq. (6.2). An implicit scheme with \(\omega = 0\) was found to be best suited for the solution of stationary flow problems (for unsteady flows see Section 6.3). Herewith, Eq. (6.2) simplifies to

As we can see, the implicit formulation leads to a set of non-linear equations for the unknown flow variables at the time \((t+\Delta t)\). The solution of Eq. (6.26) requires the evaluation of the residual at the new time level, i.e., \(\vec{R}^{n+1}\). Since we do not know \(\vec{W}^{n+1}\), this cannot be done directly. However, we can linearise the residual \(\vec{R}^{n+1}\) in Eq. (6.26) about the current time level, i.e.,

where the term \(\partial\vec{R}/\partial\vec{W}\) is referred to as the flux Jacobian. We should mention that the flux Jacobian is often derived from a rather crude approximation to the spatial discretisation represented by \(\vec{R}^n\). For example, in the case of higher-order upwind discretisations, it is quite common to base the flux Jacobian solely on a first-order upwind scheme. However, for best efficiency and robustness, the flux Jacobian should still reflect the most important features of the spatial discretisation.

If we substitute now the linearisation in Eq. (6.27) for \(\vec{R}^{n+1}\) into Eq. (6.26), we obtain the following implicit scheme

The term in square brackets on the left-hand side of Eq. (6.28) is referred to as the implicit operator or the system matrix. Consequently, the right-hand side of Eq. (6.28) is called the explicit operator. It is only the explicit operator that determines the spatial accuracy of the solution.

The implicit operator constitutes a large, sparse, and non-symmetric block matrix with dimensions equal to the total number of cells (cell-centred scheme) or grid points (cell-vertex scheme). Below we will discuss further the form of the implicit operator for structured as well as for unstructured grids. As we already saw in Section 3.2, the mass matrix \(\bar{M}\) can be replaced by the identity matrix, without influencing the steady state solution. The parameter \(\beta\) in Eq. (6.28) 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 recommended since the scheme with \(\beta = 1\) is much more robust, and the time accuracy is of no importance for steady problems.

In the case of stiff governing equations, the source term has to be included in the implicit operator. This happens quite naturally with the linearisation of the residual in Eq. (6.27), which leads to a formulation identical to Eq. (6.10). As demonstrated in Ref. [12], the above implicit scheme (6.28) remains stable for any time step if the eigenvalues of \(\partial\vec{Q}/\partial\vec{W}\) are all negative or zero.

The solution of the linear equation system (6.28) requires the inversion of the implicit operator, i.e., the inversion of a very large matrix. In principle, this can be done in two ways. The first one consists of a direct matrix inversion, using either the Gaussian elimination or some direct sparse matrix method [24], [25]. However, because of the excessive amount of memory and the very high computational effort, this approach is not suited for practical problems [26].

The second possibility of inverting the implicit operator represent iterative methods. We mentioned the most widely used ones in Subsection 3.2.2. Iterative methods can be divided roughly into two groups. The first one consists of approaches which decompose the implicit operator into several parts - a process called factorisation. The factors are constructed in such a way that they can be more easily inverted than the original implicit operator. To the second group belong schemes, which employ a Krylov-subspace method for the inversion of the implicit operator. In this case, the implicit time-stepping scheme (6.28) is usually turned into Newton's method by setting \(\Delta t \rightarrow \infty\). The scheme is then named Newton-Krylov method.

In the following, we shall discuss first the matrix structure of the implicit operator. Then, we shall investigate the possibilities of computing the flux Jacobian \(\partial\vec{R}/\partial\vec{W}\) in Eq. (6.28). Finally, we shall present the three most popular iterative methods in detail.

The Alternating Direction Implicit (ADI) scheme was one of the first iterative implicit schemes [40]. The ADI scheme can be implemented on structured grids only. It is based on an approximate splitting (or, in other words, factorisation) of the implicit operator in Eq. (6.28) into two (in 2D) or three (in 3D) factors. Each factor contains the linearisation of the convective and viscous fluxes for one particular direction in the computational space. In 3D, this leads to the formulation [40], [41]

For clarity, the node indices \(i, j, k\) were omitted from Eq. (6.46) where not required. Furthermore, the superscripts \(I, J, K\) in Eq. (6.46) mark the flux vector or of the face area associated with certain coordinate in the computational space. For example, \(\Delta S^{I}_{i+1/2,j,k}\) corresponds to \(\Delta S_2\) in Fig. 4.1b. By replacing the node indices by cell indices, the scheme in Eq. (6.46) can be applied to a cell-centred discretisation.

The derivatives of the convective and viscous fluxes in Eq. (6.46) can be evaluated as discussed in the previous subsection. The ADI method is traditionally coupled to the central scheme with artificial dissipation. In such a case, a formulation similar to that in Eq. (6.34) holds for each factor (of course, the linearisation of the source term is included in one factor only). In order to obtain a robust and efficient scheme, it is necessary to include a linearisation of the artificial dissipation terms in the implicit operator [42]-[44]. The linearisation is in general simplified by treating the spectral radii and the dissipation coefficients as independent of \(\vec{W}\). Hence, according to Eq. (4.48) the factor, e.g., in the I-direction becomes

with the implicit artificial dissipation term JST scheme, (cf. Eq. (4.50))

The implicit dissipation terms in other directions are defined in similar way. It is also possible to retain only the second-order differences in the implicit operator [43], [44]. This reduces the matrix for each of the factors from block-pentadiagonal to block-tridiagonal form (as sketched in Fig. 6.1). However, as pointed out in Ref. [44], this restricts the stability of the scheme.

The inversion of the implicit operator in Eq. (6.46) proceeds in three steps (two steps in 2D), i.e.,

where \(\mathbf{D}\) represents the diagonal, \(\mathbf{L}\) the lower-diagonal and \(\mathbf{U}\) the upper-diagonal terms, respectively. Each step requires the inversion of a block-tridiagonal or a block-pentadiagonal matrix (if \(\epsilon^{(4)}_{IM} > 0\) in Eq. (6.47)). This is done by a direct solution method. In order to reduce the numerical effort, Pulliam and Chaussee [45] suggested a diagonalised form of the ADI scheme. Herewith, the block matrices (composed of convective and viscous flux Jacobians) are transformed into diagonal matrices. Hence, only non-block tri- or pentadiagonal matrices have to be inverted, which results in significant savings of computational work and memory [43]. Although the diagonalisation is strictly valid only for Euler equations, it can be employed for viscous flows as well [43]. The linearisation of the viscous fluxes (e.g., like in Eq. (6.45)) is then either omitted in the implicit operator, or it can be approximated by the viscous eigenvalue (Eq. (6.19)).

The splitting of the implicit operator introduces what is called the factorisation error. It is the difference between the implicit operator of the base scheme in Eq. (6.28) and the factorised operator. In the case of the ADI scheme, this error term is scaled by the factor \((\Delta t)^{N}\), where \(N\) denotes the number of space dimensions. This term causes the ADI scheme to loose its unconditional stability in 3D [46]. However, the stability is improved if the fourth-order differences of the artificial viscosity are included in the implicit operator [44]. In fact, the ADI scheme was successfully used for the solution of various 3-D problems [43], [47]. An interesting implementation of the ADI scheme on multiblock grids, which treats the block boundaries implicitly, was presented by Rosenfeld et al. [48].

The time step \(\Delta t\) can be computed in the same way as presented in Subsection 6.1.4, using Eq. (6.14). The optimal CFL number varies between 20 and 50 depending on the flow case.

The implicit Lower-Upper Symmetric Gauss-Seidel (LU-SGS) scheme, which is also called the Lower-Upper Symmetric Successive Overrelaxation (LU-SSOR) scheme, became widely-used because of its low numerical complexity and modest memory requirements, which are both comparable to an explicit multistage scheme. Furthermore, the LU-SGS scheme can be implemented easily on vector and parallel computers. It can also be used on structured as well as on unstructured grids.

The LU-SGS scheme has its origins in the work of Jameson and Turkel [49], who considered decompositions of the implicit operator into lower and upper diagonally dominant factors. The LU-SGS method itself was introduced by Yoon and Jameson [50]-[52] as a relaxation method for solving the unfactored implicit scheme in Eq. (6.28). It was further developed and applied to 3-D viscous flow fields by Rieger and Jameson [53]. Since then, various researchers applied the LU-SGS scheme to viscous flows on structured [54]-[61] and on unstructured grids [62]-[66]. The LU-SGS approach is often used for the simulation of chemically reacting flows [67]-[71].

The LU-SGS scheme employs a simplification of the first-order accurate flux-vector splitting approach due to Steger and Warming (see Subsection 6.2.2) for the linearisation of the convective fluxes in Eq. (6.29). The linearisation is always kept the same regardless of the discretisation of the explicit operator. The LU-SGS scheme is further based on the factorisation of the implicit operator in Eq. (6.28) into the following three parts

The factors are constructed such that \(\mathbf{L}\) consists only of terms in the strictly lower triangular matrix, \(\mathbf{U}\) of terms in the strictly upper triangular matrix and \(\mathbf{D}\) of diagonal terms. It is important to remark that the number of factors remains always the same independent of the number of space dimensions.

The system matrix of the LU-SGS scheme (Eq. (6.49))) can be inverted in two steps - a forward and a backward sweep, i.e.,

with \(\vec{W}^{n+1} = \vec{W}^{n} + \Delta\vec{W}^{n}\). The operators \(\mathbf{L}\), \(\mathbf{D}\), and \(\mathbf{U}\) and also the inversion procedure differ on structured and unstructured grids in some respects. Therefore, we shall discuss below each case separately.

On structured grids, the operators are defined as (see [50]-[52], and [68], [55])

For better readability, only those node indices (or cell indices in the case of a cell-centred scheme) are shown in Eq. (6.51), which differ from \(i, j, k\). The superscripts \(i, j, k\) at \(\Delta S\) indicate the direction in the computational space. The unit normal vectors in the positive/negative flux Jacobians \(\bar{A}^{\pm}\) and in the viscous flux Jacobians \(\bar{A}_{v}\) are evaluated at the same side of the control volume like the associated face areas \(\Delta S\). Note that the unit normal vectors are assumed to point outwards of the control volume. In contrast, in various references it is supposed that the unit normal vectors from opposite sides of the control volume point in the same direction.

The viscous flux Jacobians in Eq. (6.51) are either computed numerically, or are replaced by their TSL approximation, corresponding to Eq. (6.45). It is possible to apply the TSL approximation in all computational coordinates, regardless of the actual orientation of the boundary layer(s). A further simplification consists of substituting the viscous flux Jacobians by the viscous spectral radii (Eq. (6.19)), i.e., \(\bar{A}_{v}\Delta S \approx \hat{\Lambda}_{v}\), as suggested in [63].

The split convective flux Jacobians \(\bar{A}^{\pm}\) are constructed in such a way that the eigenvalues of the (+) matrices are all non-negative, and of the (-) matrices are all non-positive. In general, the matrices are defined as [49]

where \(\bar{A}_{c}\) stands for the convective flux Jacobian (Section A.9) and \(\hat{\Lambda}_{c}\) represents the spectral radius of the convective flux Jacobian (given by Eq. (4.53) or Eq. (6.15)), respectively. Note the similarity between the above approximation (6.52) and Eq. (6.40), when the derivatives of \(\bar{A}_{c}\) are neglected. The factor \(\omega\) in Eq. (6.52) represents an overrelaxation parameter. It also determines the amount of implicit dissipation and hence influences the convergence properties of the scheme. The factor can be chosen in the range \(1 < \omega \le 2\). Higher values of \(\omega\) increase the stability of the LU-SGS scheme, but may slow down the convergence to steady state. The definition of the Jacobians \(\bar{A}^{\pm}\) in Eq. (6.52) ensures a diagonally dominant system matrix, which is very important for the efficiency and robustness of the iterative inversion procedure (6.50).

The splitting according to Eq. (6.52) together with averaged face vectors allow a simplified evaluation of the diagonal operator \(\mathbf{D}\)

The spectral radii of the convective flux Jacobians \(\hat{\Lambda}_{c}\) are given in Eq. (6.15). The face areas and normal vectors are averaged in the respective I-, J-, or K-direction according to Eq. (6.16). As we shall see immediately, this approximation helps to reduce the operation count and the memory requirements significantly.

A distinguishing feature of the LU-SGS method is how the forward and the backward sweep in Eq. (6.50) are carried out. In 2D, the sweeps are accomplished along diagonal lines \((i+j) = \mathrm{const.}\) in computational space. This is depicted in Fig. 6.5 for the forward sweep (first line of Eq. (6.50))). In this way, the off-diagonal terms involved in the \(\mathbf{L}\) and the \(\mathbf{U}\) operator become known from the previous part of a sweep (denoted by crosses in Fig. 6.5). In 3D, the implicit operator is inverted on \(i+j+k = \mathrm{const.}\) planes, as sketched in Fig. 6.6. Hence, the LU-SGS scheme can be written as

As we can see from Eq. (6.54), the only term which needs to be inverted is the diagonal term \(\mathbf{D}\). Thus, the LU-SGS methodology transforms the inversion of a sparse banded matrix into the inversion of a block-diagonal matrix. Furthermore, if the viscous flux Jacobians in Eq. (6.53) are approximated by the viscous spectral radii, the operator \(\mathbf{D}\) becomes a diagonal matrix (except for the source term). Hence, the LU-SGS scheme requires a very small computational effort as compared to other implicit schemes (e.g., the ADI scheme discussed previously). Furthermore, the inversion of the diagonal operator can be carried out independently for each node (cell) of the diagonal plane, which makes the scheme easy to vectorise. The indices of the nodes/cells on the diagonal planes can be obtained with the following pseudo-code [72]:

The number of diagonal planes is: nplanes = imax + jmax + kmax - 2. Obviously, the above code can be optimised for higher computational efficiency.

In order to avoid explicit evaluation and storage of the convective flux Jacobians in \(\mathbf{L}\) and \(\mathbf{U}\), the products \(\bar{A}^{\pm}\Delta\vec{W}^{n}\) can be substituted by Taylor series expansion of the fluxes [53]. Using Eq. (6.52), we can write

with the update of the convective fluxes

The simplification given by Eq. (6.55) is possible due to the sweeping along diagonal planes, since \(\vec{F}^{n+1}_{c}\) is then known. This leads to a further significant decrease of the numerical effort of the LU-SGS scheme.

The time step \(\Delta t\) can be computed in the same way as presented in Subsection 6.1.4, using Eq. (6.14). However, it should be noted that the implicit LU-SGS scheme in Eq. (6.49) represents an approximate Newton iteration in the case of \(\Delta t \to \infty\) as stated by Rieger and Jameson [53]. Thus in general, CFL numbers of the order of \(10^{4}\) to \(10^{6}\) are used in practice for stationary flows. The convergence is then controlled by the overrelaxation parameter \(\omega\). For the simulation of unsteady flows, we may employ the formulation presented below in Section 6.3. Another possibility is to use the modified version of the LU-SGS scheme described in Ref. [73].

Here, the operators read for a median-dual scheme [63]-[65]

In Eq. (6.57), \(L(i)\), and \(U(i)\) denote the nearest neighbours of node \(i\) which belong to the lower (upper) matrix, \(\Delta S_{ij}\) represents the face area associated with the edge \(ij\) (see Fig. 5.9), and \(N_{F}\) stands for the number of faces of the control volume \(\Omega_{i}\), respectively. The spectral radius of the convective fluxes \((\hat{\Lambda}_{c})_{i}\) is computed by Eq. (6.21). The viscous flux Jacobian \(\bar{A}_{v}\) can be again approximated by its spectral radius [63]. In this case, the diagonal operator becomes

where \((\hat{\Lambda}_{v})_{i}\) is evaluated according to Eq. (6.21). Formulae similar to Eq. (6.57) and (6.58) can be obtained in the case of the cell-centred scheme. The major difference is that the summation in the \(\mathbf{L}\) and the \(\mathbf{U}\) operator is conducted over faces of the cell instead of incident edges.

The sets \(L(i)\) and \(U(i)\) in Eq. (6.57) should fulfil the same function as the diagonal planes on structured grids. For this reason, it is necessary to arrange the nodes (cells) into layers such that [63]:

• nodes \(i\) (cells \(I\)) of a current layer have connections to layers with previously updated flow variables - otherwise the LU-SGS scheme degenerates to a Jacobi iteration,
• nodes (cells) in a layer are not connected to each other - otherwise the scheme could not be vectorised.

The layers can be generated for the median-dual scheme with a procedure described in Ref. [63]. Approaches for the cell-centred scheme were suggested in [74], [75].

An appropriate definition of the sets \(L(i)\) and \(U(i)\) allows for the following two-step inversion procedure [63]-[65]

where the viscous Jacobians were approximated by their spectral radii and where the positive/negative Jacobians were linearised according to Eq. (6.55). Furthermore, the factor \((r^{*}_{A})_{j}\) is defined as

with \(\|\vec{r}_{j} - \vec{r}_{i}\|_{2}\) being the length of the edge \(ij\).

The time step can be computed in the same way as for the explicit scheme (Eq. (6.20))). However, the viscous eigenvalue should be omitted, since the viscous terms are already contained in the implicit operator and hence do not reduce the time step as in the case of an explicit scheme. The CFL number can be chosen in the range from \(10^{4}\) to \(10^{6}\) for steady flows. The convergence speed and the robustness of the LU-SGS scheme can then be tuned using the overrelaxation parameter \(\omega\).

First of all, let us rewrite the implicit scheme given by Eq. (6.28) as

where \(\bar{J}\) represents the implicit operator (system matrix). As we already saw, \(\bar{J}\) constitutes a large, sparse, and generally non-symmetric matrix. In the previous subsections, we discussed two methods that decompose \(\bar{J}\) into several factors which can be each more easily inverted than \(\bar{J}\) itself. However, due to the factorisation error (and approximate linearisation of \(\vec{R}^{n+1}\)), only a linear convergence to steady state can be achieved. In order to obtain the quadratic convergence of Newton's method for the solution of non-linear equations, four conditions must be fulfilled:

• the linearisation of the residual must be exact,
• \(\bar{J}\) must be accurately inverted,
• the time step has to be \(\Delta t \to \infty\),
• initial solution must be, in some sense, close to the final solution.

Obviously, the main obstacles that have to be overcome are the linearisation and the inversion of the full system matrix.

A particularly suitable class of iterative techniques for the solution of large linear equation systems are the so-called Krylov-subspace methods. Several were proposed for the inversion of matrices which arise in CFD. Examples are the Conjugate Gradient Squared (CGS) method [76], the Bi-Conjugate Gradient Stabilised (Bi-CGSTAB) scheme [77], or the Transpose-Free Quasi-Minimum Residual (TFQMR) approach [78]. However, the most successful Krylov subspace method became the Generalised Minimal Residual (GMRES) technique, which was originally suggested by Saad and Schulz [79], [80]. Since then, the GMRES method was improved and augmented by several researchers [81]-[84]. Because of its popularity, we shall focus on the GMRES approach in the following. Nevertheless, most of what we shall discuss also applies to the other Krylov subspace methods.

As we mentioned in Subsection 3.2.2 (see also Appendix A.12), the GMRES method minimises the norm of the global residual, i.e., \(\|\bar{J}\Delta\vec{W}^{n} + \vec{R}^{n}\|\) over a set of \(m\) orthonormal vectors (search directions), which span the Krylov subspace \(\mathcal{K}_{m}\) given by Eq. (3.10). The GMRES algorithm can be summarised as follows:

• guess a starting solution \(\Delta\vec{W}^{n}_{0}\) and evaluate the initial residual vector \(\vec{r}_{0} = \bar{J}\Delta\vec{W}^{n}_{0} + \vec{R}^{n}\),
• generate the \(m\) search directions (by Gram-Schmidt orthogonalisation),
• solve the minimisation problem,
• form an approximate solution of Eq. (6.61) as \(\Delta\vec{W}^{n} = \Delta\vec{W}^{n}_{0} + \vec{y}_{m}\).

Since the memory requirements increase linearly with the number of search directions, \(m\) is restricted to values between 10 and 40 in practice. This might not be sufficient for a converged solution \(\Delta\vec{W}^{n}\). Thus, the GMRES method has to be restarted, i.e., we set \(\Delta\vec{W}^{n}_{0} = \Delta\vec{W}^{n}\), compute \(\vec{r}_{0}\) and proceed with step 2. As pointed out in Ref. [85], instead of working with a constant number of search directions, \(m\) should be reduced if the norm of the global residual drops below a specified tolerance. In this way, a large number of operations can be saved in later stages of the Newton iteration.

GMRES and other Krylov subspace methods allow us to circumvent an explicit computation and storage of the flux Jacobian \(\partial\vec{R}/\partial\vec{W}\). The idea is based on the observation that the methods rely only on matrix-vector products of the form \(\bar{J}\Delta\vec{W}^{n}\) and do not need the matrix \(\bar{J}\) explicitly. The product of the flux Jacobian with the solution update can be approximated by a simple finite difference as

which requires only two evaluations of the residual. The stepsize \(h\) has to be chosen with some care, in order to minimise the numerical error [36]. One particularly suitable formulation reads [86]

where \(\epsilon\) denotes the machine accuracy, \(d\) the scalar product \(\vec{W}^{n}\cdot\Delta\vec{W}^{n}\), \(|\Delta\vec{W}^{n}|\) is the vector \(\Delta\vec{W}^{n}\) with all elements set to their absolute values, and finally \(\mathrm{typ}\,U\) represents a typical size of \(U\). Apart from saving memory and operations, there is an even more important advantage of the finite-difference approximation. Namely, numerically accurate linearisation of a high-order residual \(\vec{R}^{n}\) (including boundary conditions, limiters, source terms, etc.) can be easily achieved.

Thus, the quadratic convergence of Newton's scheme can be realised at moderate costs. For this reason, we speak of such a scheme as of Newton-Krylov approach [86]-[88].

The efficiency of Krylov-subspace methods depends strongly on a good preconditioner. Its purpose is to cluster the eigenvalues of the system matrix \(\bar{J}\) around unity. Thus, instead of Equation (6.61), the left- or right-preconditioned system according to Eq. (3.11) is solved. Using Eqs. (6.61) and (6.62) together with the condition \(\Delta t \to \infty\), the Newton-Krylov method becomes

with left preconditioning, and

in the case of right preconditioning, respectively. The main difference between the two preconditioning methodologies is that left preconditioning scales the residual \(\vec{R}^{n}\) whereas right preconditioning does not. This has to be kept in mind when the convergence of the Krylov method is monitored.

Obviously, the preconditioner should be as close as possible to the inverse of the system matrix (\(\bar{P}_{L,R} \approx \bar{J}^{-1}\)). But on the other hand, it should be invertible with low numerical effort. Therefore, we have to find an optimal tradeoff between the convergence speed of the Krylov method and the time spend for inverting the preconditioning matrix. One of the most successful preconditioners is the Incomplete Lower Upper factorisation method [89], [90] with varying level of fill-in (mostly with zero, designated as ILU(0)). The ILU preconditioner is especially efficient in the case of viscous, turbulent flows, i.e., for stiff equations [91]. In order to obtain a good performance on unstructured grids, it is necessary to reorder the elements of the system matrix such that the bandwidth is reduced. We already discussed the RCM renumbering strategy [27], [28] in Section 6.2.1.

A serious disadvantage of the ILU preconditioning scheme is that elements of the matrix \(\bar{J}\) have to be computed (see Subsection 6.2.2) and stored. Therefore, some authors suggested to employ the LU-SGS scheme as a preconditioner [92], [93]. Hence, when Eq. (6.62) is applied, the formation and storage of \(\bar{J}\) is completely avoided. However, it was demonstrated in Ref. [91] on behalf of several 2-D cases that the GMRES method with LU-SGS is inferior to GMRES combined with ILU(0) in terms of the CPU-time. Nevertheless, the LU-SGS scheme (best when coupled with multigrid) still represents an attractive alternative, particularly in 3D.

The time step of the implicit Newton-Krylov method is infinitely large. However, it is advisable to use small time steps at the beginning of the Newton iteration process. The reason is that the flow solution is in general far from the steady state at the beginning of the solution process, i.e., the root of the nonlinear equation

and this may cause a breakdown of the Newton iteration. One possible remedy is the so-called Switched Evolution Relaxation (SER) technique [94]. Here, the term \(\Omega/\Delta t\) is retained in Eq. (6.61). The time step is evaluated in the same way as presented for the explicit scheme (Eq. (6.14) or Eq. (6.20) without the viscous eigenvalue). The CFL number \(\sigma\) is increased starting from a small initial value correspondingly to the reduction of the 2-norm of the residual, i.e.,

Hence, the convergence of the iteration procedure (6.61) will be at first linear, but it approaches the quadratic convergence of Newton's method for large CFL numbers. A further effect of the time term is the increased diagonal dominance of \(\bar{J}\) (inverse proportional to \(\Omega/\Delta t\)), which will help to stabilise the iteration. In Ref. [88], it was suggested to clip \(\sigma^{n+1}\) in Eq. (6.66) such that it increases by maximum factor of two and decreases less than factor of ten.

A further approach which can be used to overcome the start-up problems of Newton's method consist of grid sequencing, where the initial solution is obtained on a sequence of coarser grids and interpolated onto finer grids. Another possibility is to use a numerically cheap but robust iteration scheme for the initial guess. For example, we could start with a multigrid scheme driven by the LU-SGS method and then switch to GMRES with LU-SGS as preconditioner. This might be particularly interesting for viscous turbulent flows, where the convergence of a multigrid scheme usually slows down after the initial phase. However, the global flow solution is then already close to the steady state.

The simulation of unsteady flow phenomena is becoming increasingly important in many engineering disciplines. Examples are the interaction between stationary and rotating parts in turbomachinery, piston engines, fluid-structure interaction, helicopter aerodynamics, aeroacoustics, DNS or LES of turbulent flows, detonations, etc. Clearly, the simulation has to be conducted efficiently and with accuracy adequate for the problem being solved.

Explicit schemes represent the best choice for certain unsteady applications when the time scales are comparable to the spatial scales over the eigenvalue, i.e., when the CFL number dictated by the physics is of the order of unity. This is for example the case in aeroacoustics, DNS, and LES. For such applications, explicit Runge-Kutta schemes are quite popular. Since the global physical phenomena evolve much slower than the solution changes locally in these applications, it is necessary to integrate over a long period of physical time. In order to do this accurately, the temporal resolution of the explicit scheme has to be of 3rd or higher order. This requires the use of Runge-Kutta methods different to that we presented in Section 6.1 (see, e.g., [95], [96]).

In other cases, where the physical time scales are large in comparison to the spatial scales divided by the eigenvalue (e.g., flutter, rotor-stator interaction, etc.), the CFL number can be chosen in the order of several hundreds or even thousands without impairing the accuracy of the simulation. Obviously, in such cases an implicit scheme is more appropriate. In the following, we shall discuss a particular technique known as the dual time-stepping approach, which is very often employed for unsteady flows.

The dual time-stepping approach is based on the second-order time accurate version of the basic non-linear scheme in Eq. (6.2). For this purpose we set \(\beta = 1\) and \(\omega = 1/2\) in Eq. (6.2). Hence, we obtain

where \(\Delta t\) denotes the global physical time step and \(\bar{M}\) the mass matrix, respectively. Equation (6.67) constitutes a 3-point backward-difference (backward Euler) approximation of the time derivative in Eq. (6.1). In order to solve the system of non-linear equations given by Eq. (6.67), we can use either Newton's method or a time-stepping methodology. The latter can be written as

where \(\vec{W}^{*}\) is the approximation to \(\vec{W}^{n+1}\) and \(t^{*}\) denotes a pseudo-time variable. Note that there is no mass matrix in the time derivative. The unsteady residual is defined as

All terms which are constant during the time-stepping in Eq. (6.68) are gathered in a source term, i.e.,

In the case of moving and/or deforming grids, the new size of the control volume, i.e., \(\Omega^{n+1}\) in Eq. (6.68) has to satisfy the Geometry Conservation Law (see Appendix A.5 for details and references).

The stationary solution of Eq. (6.68) corresponds to the flow variables at the new time level, i.e., \(\vec{W}^{*} = \vec{W}^{n+1}\). Since \(\vec{R}^{*}_{I} = 0\) at steady state in pseudo time, Equation (6.67) is fulfilled. Any of the previously presented explicit or implicit time-marching schemes can be employed for the solution of the system of equations (6.68) in the pseudo time. In the following, we shall discuss the implementation of the dual time-stepping approach for explicit multistage and implicit schemes.