6.1 Explicit Time-Stepping Schemes 显式时间推进格式[cfd-6-1]

The concept of explicit multistage schemes was first presented by Jameson et al. [1]. The multistage scheme advances the solution in a number of steps - so-called stages - which can be viewed as a sequence of updates according to Eq. (6.4). Applied to the discretised governing equations (6.1), where the mass matrix was lumped, an \(m\)-stage scheme reads

In the above expressions (6.5), \(\alpha_k\) represents the stage coefficients. Furthermore, the denotation \(\vec{R}_I^{(k)}\) means that the residual is evaluated with the solution \(\vec{W}_I^{(k)}\) from the \(k\)-th stage.

Unlike in the classical Runge-Kutta schemes, only the zeroth solution and the newest residual are stored here in order to reduce the memory requirements. The stage coefficients can be tuned to increase the maximum time step and to improve the stability for a particular spatial discretisation [2]-[4]. For consistency, it is only required that \(\alpha_m = 1\). A consequence of the modification to the Runge-Kutta scheme is that second-order time accuracy can be realised only if \(\alpha_{m-1} = 1/2\). Otherwise, the multistage scheme is first-order accurate in time.

The above multistage approach (6.5) is particularly suitable for upwind spatial discretisation on structured as well as unstructured grids. Central discretisation schemes perform more efficiently with the hybrid multistage methodology, which will be described next. Sets of optimised stage coefficients for first- and second-order upwind schemes are presented in Table 6.1 for three- to five-stage schemes [2]. Practical experience shows that the coefficients for the first-order scheme should be preferred in cases, where the flow field contains strong shocks, regardless of the order of the spatial discretisation. This can be explained by the fact that every higher-order scheme switches to first order at shocks to prevent oscillations of the solution. However, the residuals at strong shocks influence the convergence to steady state most significantly.

The main disadvantage of every explicit scheme is that the time step (\(\Delta t\)) is severely restricted by the characteristics of the governing equations as well as by the grid geometry. We shall discuss the computation of the maximum allowable time step in Subsection 6.1.4. Theoretical aspects of the determination of the time step and the so-called CFL number will be considered in Section 10.3 on stability analysis.

The computational work of an explicit multistage scheme (6.5), applied to the system (6.1), can be substantially reduced if the viscous fluxes and the dissipation are not re-evaluated at each stage. Additionally, the dissipation terms from different stages can be blended to increase the stability of the scheme. Methods of this type were devised by Martinelli [5] and by Mavriplis et al. [6]. They are known as hybrid multistage schemes. Provided the stage coefficients are carefully optimised, the hybrid schemes are as robust as the basic multistage schemes.

Let us for illustration consider a popular 5-stage hybrid scheme, where the dissipative terms are evaluated at odd stages - generally denoted as the (5,3)-scheme. First, we split the spatial discretisation into two parts, i.e.,

The first part, \(\vec{R}_c\), contains the central discretisation of the convective fluxes, which can be either the average of variables or the average of fluxes. It also includes the source term. The second part, \(\vec{R}_d\), is composed of the viscous fluxes and the numerical dissipation. For example, in the case of the central scheme with artificial dissipation (Subsections 4.3.1 or 5.3.1) we would set

where \(\vec{W}_{av}\) represents the arithmetic average of flow variables from the left and the right side of face \(k\).

With the residual split according to Eq. (6.6), the (5,3)-scheme can be formulated as

where

The stage coefficients \(\alpha_m\) and the blending coefficients \(\beta_m\) in the above relations (6.7), (6.8) are given in Table 6.2 for central and upwind schemes. Both sets of coefficients are particularly optimised for the multigrid method (Section 9.4). We shall discuss the properties of the above hybrid multistage scheme later in Section 10.3.

It should be mentioned that it is also popular to evaluate the dissipation term \(\vec{R}_d\) in the first two stages only, without any blending. A well-known (5,2)-scheme, which is often employed with the central spatial discretisation, uses the stage coefficients of Table 6.2. However, the (5,2)-scheme is less suitable for viscous flows and multigrid than the above (5,3)-scheme.

There are certain cases in which the source term \(\vec{Q}\) in Eq. (4.2) or (5.2) becomes dominant. Such situation is often encountered when chemistry or turbulence models are employed. The problem is that a large source term changes the flow variables rapidly in space and in time. The changes due to a strong source term happen at much smaller time scales than those of the flow equations. This increases the stiffness of the governing equations significantly. The stiffness is defined as the ratio of the largest to the smallest eigenvalue of the Jacobian matrix \(\partial\vec{R}/\partial\vec{W}\). The stiffness can also be viewed as the ratio of the largest to the smallest time scale.

When we apply one of the above explicit multistage schemes (or any other purely explicit scheme) to a stiff system of equations, we will have to reduce the time step considerably in order to stabilise the time integration. Hence, the convergence to the steady state will become very slow. More seriously, an explicit scheme can even fail to find the correct solution [7]. A remedy suggested by Curtiss et al. [8] is to treat the source term in an implicit way. In order to demonstrate the approach, we rewrite the basic explicit scheme in Eq. (6.4) as follows (cf. Eq. (4.2) or (5.2))

where the source term is now evaluated at the new time level \((n+1)\). For simplicity, the mass matrix was omitted from Eq. (6.9). Since the value of the source term at the time \((n+1)\) is unknown, we have to approximate it. For this purpose, we linearise the source term about the current time level \(n\), resulting in

If we insert Eq. (6.10) into Eq. (6.9) and rearrange the terms, we obtain the following relation [9], [10]

where \(\bar{I}\) represents the identity matrix. The formulation (6.11) is called point implicit because the term in square brackets on the left-hand side - the implicit operator - depends only on values in the control volume \(\Omega_I\) itself. A comparison with Eq. (6.9) reveals that the scalar time step \(\Delta t\) changed now to a matrix. Thus, each flow equation becomes scaled by an individual parameter, corresponding to the associated eigenvalue. In this way, the disparity between the time scales is offset and the time step restriction due to the source term is alleviated.

When we apply the above point-implicit approach to the multistage scheme in Eq. (6.5), we obtain for the \(k\)-th stage

with \(\vec{R}_Q\) defined in Eq. (6.9). A similar expression holds also for the hybrid multistage scheme in Eq. (6.7). The interested reader may find a detailed investigation of the influence of the source term on stability in Refs. [11]-[13].

A more elaborate approach is to treat the stiff source term by means of an semi-implicit Runge-Kutta scheme, originally suggested by Rosenbrock [14]. The scheme is numerically more stable than the approach in Eq. (6.12), however it requires the inversion of a large system of linear equations at each stage. Details of the semi-implicit Runge-Kutta scheme and applications can be found in Refs. [15]-[18].

Every explicit time-stepping scheme remains stable only up to a certain value of the time step \(\Delta t\). To be stable, a time-stepping scheme has to fulfil the so-called Courant-Friedrichs-Lewy (CFL) condition [19]. It states that the domain of dependence of the numerical method has to include the domain of dependence of the partial differential equation. The CFL condition means for the basic explicit scheme (6.4) that the time step should be equal to or smaller than the time required to transport information across the stencil of the spatial discretisation scheme. Hence, in 1D the condition for the time step would read for the linear convection equation

where \(\Delta x/|\Lambda_c|\) represents the time necessary to propagate information over the cell size \(\Delta x\) with the velocity \(\Lambda_c\). The velocity \(\Lambda_c\) corresponds to the maximum eigenvalue of the convective flux Jacobian. The positive coefficient \(\sigma\) denotes the CFL number. The magnitude of the CFL number depends on the type and the parameters of the time-stepping scheme, as well as on the form of the spatial discretisation scheme. We shall investigate the dependency of \(\sigma\) in Section 10.3 for two model problems. Tables 6.1 and 6.2 list the CFL numbers for various multistage schemes and discretisations.

The maximum time step can be determined for linear model equations with the aid of Von Neumann stability analysis (Section 10.3). However, the maximum time step can be calculated only approximately in multiple dimensions and for non-linear governing equations. In the following, we will present relations for the estimation of the time step on structured and unstructured grids, and for inviscid as well as for viscous flows.

On a structured grid, the time step \(\Delta t\) can be determined for a control volume \(\Omega_I\) from the approximate relation [20]-[22]

The CFL number \(\sigma\) is given for multistage schemes in Table 6.1 and for hybrid schemes in Table 6.2. The spectral radii of the convective flux Jacobians (A.9) read for the three grid directions

The normal vectors and face areas in Eq. (6.15) are obtained by averaging the corresponding values from the two opposite sides of the control volume in the respective direction. For example, if a dual control volume (Subsection 4.2.3) would be oriented as sketched in Fig. 4.1b, we would use in the \(I\)-direction

Similar expressions hold for the \(J\)- and \(K\)-direction. With Eq. (6.14), we obtain a local time step, which is valid for one control volume only. If we are interested in a steady state solution, we may use the local time step to accelerate the convergence (cf. Section 9.1). However, if time accuracy is important, we have to employ one global time step for all volumes, i.e.,

where the minimum over all control volumes is taken.

For viscous flows, the spectral radii of the viscous flux Jacobians (A.10) have to be included in the computation of \(\Delta t\). They can severely limit the maximum time step in boundary layers. The time step can be evaluated from [5], [21], [22]

The constant which multiplies the viscous spectral radii is usually set as \(C = 4\) for central spatial discretisations, \(C = 2\) for first-order upwind and \(C = 1\) for second-order upwind discretisations. If we assume that an eddy-viscosity turbulence model is employed, the viscous spectral radii are given by [21], [22]

and similarly for the other directions. In Equation (6.19), \(\mu_L\) denotes the laminar and \(\mu_T\) the turbulent dynamic viscosity coefficient, respectively. Furthermore, \(Pr_L\) and \(Pr_T\) are the laminar and the turbulent Prandtl numbers. The CFL numbers in Tables 6.1 and 6.2 apply also for viscous flows. Particularly efficient for viscous flows is the (5,3) hybrid scheme from Eq. (6.7).

Several approaches were suggested for the estimation of the maximum time step on unstructured grids. We shall present two different approaches below.

One proven method, which closely follows the implementation on structured grids, reads [6]

where \(\hat{\Lambda}_c\) and \(\hat{\Lambda}_v\) represent a sum of the convective and viscous spectral radii over all faces of the control volume. As on structured grids, \(1 \le C \le 4\) is usually used. In the case of the cell-centred scheme, the spectral radii are defined as [6]

The values of the flow variables at the faces of the control volume are obtained by arithmetic averaging.

The spectral radii predicted by Eq. (6.21) are too large, particularly on mixed element grids. This leads to a smaller time step than required for the stability. The implementation in Ref. [23] offers a more accurate estimation of the time step, namely

with the convective spectral radii

and with the viscous spectral radii (eddy-viscosity turbulence model assumed)

The variables \(\Delta\hat{S}^x\), \(\Delta\hat{S}^y\) and \(\Delta\hat{S}^z\), respectively, represent projections of the control volume on the \(y\)-\(z\)-, \(x\)-\(z\)- and the \(x\)-\(y\)-plane. They are given by the formulae

where \(S_x\), \(S_y\) and \(S_z\) denote the \(x\)-, \(y\)- and the \(z\)-component of the face vector \(\vec{S} = \vec{n}\cdot\Delta S\). The CFL numbers stay in general the same as on structured grids. Thus, the values collected in Tables 6.1 and 6.2 still apply. Furthermore, the convergence to steady state can also be accelerated by local time stepping, in the same way as on the structured grids. A global time step, necessary for simulating unsteady flows, can be obtained from Eq. (6.17) as before.