9.3 Residual Smoothing 残差光顺[cfd-9-3]

The standard formulation of the implicit residual smoothing reads in 3D

where \(\vec{R}^{*}_{I,J,K}\), \(\vec{R}^{**}_{I,J,K}\), and \(\vec{R}^{***}_{I,J,K}\) denote the smoothed residuals in \(I\)-, \(J\)-, and \(K\)-direction, respectively. The parameters \(\epsilon^I\), \(\epsilon^J\), and \(\epsilon^K\) stand for the smoothing coefficients in the three computational coordinates. The implicit operator in Eq. (9.5) resembles second-order central difference. The term Central Implicit Residual Smoothing (CIRS) is therefore used. The implicit system in Eq. (9.5) is solved by the Thomas algorithm for the inversion of tridiagonal matrices.

The smoothing coefficients are usually defined as functions of spectral radii of the convective flux Jacobians [11]. The purpose is to apply only as much smoothing in each coordinate direction as it is necessary for stability and good error damping. A suitable formula for 2D was suggested in [12]

Here, \(\sigma^{*}/\sigma\) denotes the ratio of the CFL numbers of the smoothed and unsmoothed scheme. The variable \(r\) stands for the ratio of the convective spectral radii (Eq. (4.53)), i.e., \(r = \hat{\Lambda}_c^{J}/\hat{\Lambda}_c^{I}\). The parameter \(\Psi \approx 0.125\) ensures linear stability of the smoothing operation.

In 3D, the smoothing coefficients can be evaluated using an expression similar to Eq. (9.6) [12]

where

etc. The typical value of the parameter \(\Psi\) is 0.0625.

The maximum of the ratio of the CFL numbers \(\sigma^{*}/\sigma\) depends on the value of the smoothing coefficient and on the type of the spatial discretisation scheme. In the case of the central scheme (Subsection 4.3.1), the ratio is given by

In practice, value of \(\sigma^{*}/\sigma \approx 2\) can be reached (\(\epsilon = 0.8\)). Higher ratios reduce the damping of the time-stepping scheme. There is no such simple condition like (9.8) for the upwind spatial discretisation since the maximum of \(\sigma^{*}/\sigma\) depends also on the stage coefficients. But normally, the CFL number (see Tables 6.1 and 6.2) can be raised by about factor two for inviscid and up to factor five for viscous flows (using the hybrid scheme from Table 6.2). Figures 9.2 and 9.3 demonstrate the effect of the central IRS on the convergence for an inviscid and viscous flow on structured grids. The comparison is done in terms of the CPU time, in order to account for the additional numerical effort due to the IRS. As we can see, the solution time can be significantly reduced especially for a viscous flow. Even larger savings can be however realised when the IRS is employed together with multigrid (Section 9.4).

The limitation of the time step due to the viscous spectral radius \(\hat{\Lambda}_v\) in Eq. (6.18) can be offset with the aid of the implicit residual smoothing. The maximum time step is calculated according to the formula (6.14) without \(\hat{\Lambda}_v\). In flow regions where the viscous spectral radius dominates, smoothing is carried out using higher coefficients \(\epsilon^I\), \(\epsilon^J\), \(\epsilon^K\). Smoothing coefficients based on the viscous spectral radii can be computed from [12], [13]

where the constant \(C \approx 5/4\). The maximum of the coefficients from Eqs. (9.7) and (9.9) is then taken as the resulting smoothing coefficient.

CIRS is implemented on unstructured grids by applying the Laplacian operator (see Eq. (5.24)) to the residual. Thus, the smoothed residual \(\vec{R}^{*}_{I}\) in a control volume \(I\) is obtained from the implicit relation [14]

The sum includes all \(N_A\) adjacent control volumes. The relation (9.10) is solved for \(\vec{R}^{*}_{I}\) using Jacobi iteration. Useful values of the smoothing coefficient are \(0.5 \le \epsilon \le 0.8\). With these values of \(\epsilon\) it is possible to increase the CFL number (and hence the time step) two to five times. Due to the diagonal dominance of the matrix, the Jacobi iteration converges in about two steps.

The previously discussed CIRS method works satisfactorily for subsonic and transonic flows. It is also helpful in viscous dominated regions. However, CIRS exhibits poor error-damping characteristics in conjunction with upwind spatial discretisation schemes. Furthermore, the robustness of a multistage scheme accelerated by CIRS suffers in the case of strong shocks. Therefore, the so-called Upwind Implicit Residual Smoothing (UIRS) method was developed [15], [16], which is particularly suited for high Mach-number flows.

By contrast to CIRS, the UIRS methodology takes the sign of the convective eigenvalues \(\bar{\Lambda}_c\) into account. The idea is to smooth the residuals only in the direction of the characteristic of the Euler equations (i.e., \(dx/dt = \mathrm{const.} = \Lambda_c\)). This approach prevents unphysical influences on the upstream residuals. The UIRS method requires transformation of the residuals into the characteristic variables (cf. Appendix A.11). In this way, each component of the residual vector can be smoothed independently according to the sign of the corresponding eigenvalue. For the \(l\)-th component of the residual, the implicit operator is defined in 1D as [15], [16]

where \(\vec{R}^{c}\) denotes the residual transformed into characteristic variables, i.e.,

The smoothed residuals \(\vec{R}^{*}\) are obtained by the solution of a tridiagonal (bidiagonal if \(\Lambda_c^l\) does not change its sign) equation system using the Thomas algorithm. Afterwards, the residuals \(\vec{R}^{*}\) are transformed back into the conservative variables and the solution can be updated.

The most desirable feature of UIRS is that it leads to very favourable damping properties of the multistage scheme, particularly in connection with an upwind spatial discretisation. The ability to damp solution errors remains or even improves for high smoothing coefficients. It was demonstrated for 1-D Euler equations that values like \(\epsilon = 500\) and \(\sigma^{*} = 1000\) result in a stable and very fast explicit multistage scheme [15], [16] (see also in analysis/mstage on the CD-ROM). However, the problem is the implementation of UIRS in multiple dimensions. Since the convective flux Jacobians cannot be diagonalised simultaneously in all coordinate directions (see Appendix A.11), the transformation Eq. (9.12) and the smoothing Eq. (9.11) have to be carried out separately for each computational coordinate. The effect of the coordinate splitting is a reduced maximum smoothing coefficient to \(2 \le \epsilon \le 6\). Despite this, the convergence to the steady state is strongly accelerated as compared to CIRS [17], [16]. The largest improvements in terms of convergence speed and robustness occur in combination with multigrid [18], [16].

The smoothing coefficients in multiple dimensions are scaled by the eigenvalues. For example, in 2D the following formula can be employed

The relation between the CFL number of the smoothed scheme and the coefficient \(\epsilon\) reads

The constant \(C\) depends on the kind of the spatial discretisation. In the case of the central scheme \(C = 1\). For the 1st- or 2nd-order upwind scheme the value is \(C = 2\).

In order to circumvent the numerical effort of the transformation to the characteristic variables, a simplified version of the UIRS method was suggested [17], [18], [16]. Written in the \(I\)-direction it reads

The Mach number \(M\) is based on velocity projected into the direction of the particular computational coordinate (here: \(I\)). Due to the low operation count, the simplified UIRS is especially suitable for 3-D flow problems. In Ref. [16] it was demonstrated that the CPU time needed to reach the steady state can be halved as compared to CIRS (hypersonic flow past a blunt cylinder, \(M_\infty = 8\)).