6.1.3 Treatment of the Source Term 源项的处理[cfd-6-1-3]

在某些情形下ï¼方程(4.2)或(5.2)中的源项 $\vec{Q}$ 会变得占主导地位。采用化学或湍流模型时经常会遇到这种情况。问题在于ï¼大的源项会使流动变量在空间和时间上迅速变化。强源项引起的变化发生在比流动方程的时间尺度小得多的时间尺度上。这会显著增大控制方程的刚性(stiffness)。刚性定义为雅可比矩阵 $\partial\vec{R}/\partial\vec{W}$ 的最大特征值与最小特征值之比。刚性也可以看作最大时间尺度与最小时间尺度之比。en

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.

把上述显式多步格式(或任何其他纯显式格式)应用于刚性方程组时ï¼为使时间积分稳定ï¼必须大幅减小时间步长ï¼于是向定态的收敛会变得非常慢。更严重的是ï¼显式格式甚至可能找不到正确的解[7]。Curtiss等人[8]提出的补救办法是以隐式方式处理源项。为演示这一方法ï¼把方程(6.4)的基本显式格式改写如下(参见方程(4.2)或(5.2))en

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))

\begin{aligned}\frac{\Omega_I}{\Delta t_I}\Delta\vec{W}_I^n &= -\left[\sum_{k=1}^{N_F}\left(\vec{F}_c^n - \vec{F}_v^n\right)_k\Delta S_k - \Omega_I\vec{Q}_I^{n+1}\right] \\&= -\left(\vec{R}_Q\right)_I^n + \Omega_I\vec{Q}_I^{n+1},\end{aligned} \tag{6.9}

其中源项现在在新时间层 $$(n+1)$$ 处求值。为简单起见ï¼方程(6.9)中略去了质量矩阵。由于时刻 $$(n+1)$$ 的源项值未知ï¼必须对它作近似。为此ï¼把源项在当前时间层 $$n$$ 处线性化ï¼得到en

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

\vec{Q}^{n+1} \approx \vec{Q}^n + \frac{\partial\vec{Q}}{\partial\vec{W}}\Delta\vec{W}^n. \tag{6.10}

把方程(6.10)代入方程(6.9)并整理各项ï¼得到如下关系[9]、[10]en

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

\left[\frac{1}{\Delta t_I}\bar{I} - \left(\frac{\partial\vec{Q}}{\partial\vec{W}}\right)_I\right]\Delta\vec{W}_I^n = -\frac{1}{\Omega_I}\left[\left(\vec{R}_Q\right)_I^n - \Omega_I\vec{Q}_I^n\right], \tag{6.11}

其中 $\bar{I}$ 表示单位矩阵。式(6.11)称为点隐式(point implicit)ï¼因为左端方括号内的项——即隐式算子——只依赖控制体 $\Omega_I$ 自身的值。与方程(6.9)比较可知ï¼标量时间步长 $\Delta t$ 现在变成了一个矩阵。这样ï¼每个流动方程都被一个单独的参数(对应于相应的特征值)所缩放。如此一来ï¼时间尺度之间的差异得到补偿ï¼源项造成的时间步长限制得以缓解。en

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.

把上述点隐式方法应用于方程(6.5)的多步格式ï¼对第 $$k$$ 步得到en

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

\vec{W}_I^{(k)} = \vec{W}_I^{(0)} - \left[\left(\vec{R}_Q\right)_I^{(k-1)} - \Omega_I\vec{Q}_I^{(0)}\right]\left[\frac{\Omega_I}{\alpha_k\Delta t_I}\bar{I} - \left(\frac{\partial\vec{Q}}{\partial\vec{W}}\right)_I\right]^{-1} \tag{6.12}

其中 $\vec{R}_Q$ 由方程(6.9)定义。对方程(6.7)的混合多步格式也有类似的表达式。关于源项对稳定性影响的详细研究ï¼有兴趣的读者可参阅文献[11]-[13]。en

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

更精细的做法是用半隐式Runge-Kutta格式来处理刚性源项ï¼该格式最早由Rosenbrock[14]提出。这种格式在数值上比方程(6.12)的方法更稳定ï¼但要求在每个阶段对一个大型线性方程组求逆。半隐式Runge-Kutta格式的细节及应用可参见文献[15]-[18]。en

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