9.5.2 Implementation 实现[cfd-9-5-2]

由于我们以积分形式求解控制方程，必须把式(9.50)写成en

Since we solve the governing equations in their integral form, we have to write Eq. (9.50) as

\[\frac{\partial}{\partial t}\int_{\Omega}\vec{W}\,d\Omega + \bar{P}\bar{\Gamma}^{-1}\oint_{\partial\Omega}\left(\vec{F}_c - \vec{F}_v\right)dS = \bar{P}\bar{\Gamma}^{-1}\int_{\Omega}\vec{Q}\,d\Omega\,. \tag{9.51}\]

需要特别注意的是，对于非定常流动，式(9.51)不是守恒的。因此，为了获得时间精确的解，必须采用双时间步进方法（见6.3节）。en

It is important to note that Eq. (9.51) is not conservative for unsteady flows. Therefore, the dual time-stepping approach (see Section 6.3) has to be employed in order to obtain a time-accurate solution.

采用显式多级格式（见6.1节所述）求解预条件控制方程(9.51)时，按以下步骤进行：

  • 1. 基于矩阵\(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\)的谱半径计算新的时间步长。
  • 2. 利用谱半径（如JST格式）或预条件通量雅可比矩阵的特征值和特征向量（如Roe上风格式）计算人工耗散。耗散项既可以在守恒变量下构造，也可以在原始变量下构造[90]。
  • 3. 计算对流通量（或者不作改动，或者基于新变量\(\vec{W}_p\)）。
  • 4. 把耗散通量和对流通量相加。根据人工耗散的构造方式，或者将整个残差，或者仅将对流项乘以守恒变量预条件矩阵，即乘以\(\bar{P}\bar{\Gamma}^{-1}\)。
  • 5. 将残差乘以\(\alpha_k\Delta t/V\)。
  • 6. （可选）执行隐式残差光顺。
  • 7. 从旧守恒变量\(\vec{W}^{(0)}\)中减去残差，以得到新守恒变量\(\vec{W}^{n+1}\)。
  • 8. 更新边界条件（注意：所有基于特征变量的边界条件——如入流、出流、远场——都需要修改）。
en

The solution of the preconditioned governing equations (9.51) using an explicit multistage scheme (described in Section 6.1) proceeds according to the following steps:

  • 1. Compute a new time step based on the spectral radii of the matrix \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\).
  • 2. Evaluate artificial dissipation using the spectral radii (e.g., JST scheme) or the eigenvalues and eigenvectors of the preconditioned flux Jacobian (e.g., Roe upwind scheme). The dissipation can be formulated either in the conservative or in the primitive variables [90].
  • 3. Compute the convective fluxes (either without change or based on the new variables \(\vec{W}_p\)).
  • 4. Sum up the dissipative and convective fluxes. Depending on the formulation of the artificial dissipation, either the whole residual or just the convective terms are multiplied by the conservative variable preconditioning matrix, i.e. by \(\bar{P}\bar{\Gamma}^{-1}\).
  • 5. Multiply the residual by \(\alpha_k\Delta t/V\).
  • 6. Carry out the implicit residual smoothing (optionally).
  • 7. Subtract the residuals from the old conservative variables \(\vec{W}^{(0)}\) in oder to obtain the new conservative variables \(\vec{W}^{n+1}\).
  • 8. Update the boundary conditions (note that all boundary conditions which are based on the characteristic variables - like inflow, outflow, farfield - need to be changed).

隐式格式遵循类似的步骤，只是省略第5步和第6步，并且守恒变量以不同的方式更新（参见6.2节）。还应注意，由于雅可比矩阵为\(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\)，预条件必须包含在隐式算子中。en

Similar procedure is followed for an implicit scheme, only the steps 5. and 6. are omitted and the conservative variables are updated in a different way (cf. Section 6.2). It should also be noted that the preconditioning has to be included in the implicit operator since the Jacobian is \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\).

预条件后的标量耗散格式（JST）取如下形式（见式(4.50)）en

The preconditioned scalar dissipation scheme (JST) takes the form (see Eq. (4.50))

\[\begin{aligned} \vec{D}_{I+1/2} = \hat{\Lambda}^S_{I+1/2}\left(\bar{\Gamma}\bar{P}^{-1}\right)_{I+1/2}\Big[&\epsilon^{(2)}_{I+1/2}\left(\vec{W}_{I+1} - \vec{W}_I\right)\\ &- \epsilon^{(4)}_{I+1/2}\left(\vec{W}_{I+2} - 3\vec{W}_{I+1} + 3\vec{W}_I - \vec{W}_{I-1}\right)\Big]\,. \end{aligned} \tag{9.52}\]

注意，如果之后要把整个残差乘以\(\bar{P}\bar{\Gamma}^{-1}\)，就需要项\(\bar{\Gamma}\bar{P}^{-1}\)，因为谱半径\(\hat{\Lambda}^S\)已经包含了预条件矩阵（因此它与式(4.53)不同）。在式(9.52)中，也可以把\(\bar{P}^{-1}\)与\(\vec{W}\)合并为\(W_p\)。这样，我们可以用原始变量表述预条件的标量耗散格式：en

Note that the term \(\bar{\Gamma}\bar{P}^{-1}\) is required if the complete residual is later multiplied by \(\bar{P}\bar{\Gamma}^{-1}\), since the spectral radius \(\hat{\Lambda}^S\) already contains the preconditioning matrix (thus it is different from Eq. (4.53)). It is also possible to combine \(\bar{P}^{-1}\) and \(\vec{W}\) into \(W_p\) in Eq. (9.52). Then, we can formulate the preconditioned scalar dissipation scheme in primitive variables as

\[\begin{aligned} \vec{D}_{I+1/2} = \hat{\Lambda}^S_{I+1/2}\Big[&\epsilon^{(2)}_{I+1/2}\left(\vec{W}_{p,I+1} - \vec{W}_{p,I}\right)\\ &- \epsilon^{(4)}_{I+1/2}\left(\vec{W}_{p,I+2} - 3\vec{W}_{p,I+1} + 3\vec{W}_{p,I} - \vec{W}_{p,I-1}\right)\Big]\,. \end{aligned} \tag{9.53}\]

必须认识到，此时只有对流通量乘以\(\Gamma^{-1}\)，然后整个残差再乘以\(\bar{P}\)。式(9.53)的非守恒形式比关系式(9.52)稍微简单一些，但对于内流或含激波的流动存在精度问题。en

It is important to realize that in this case only the convective fluxes are multiplied by \(\Gamma^{-1}\) and the complete residual then by \(\bar{P}\). The nonconservative form in Eq. (9.53) is somewhat simpler than the relation (9.52), however there are problems with the accuracy for internal flows or for flows containing shocks.

预条件后的Roe上风格式可以写成（参见式(4.91)）en

The preconditioned Roe upwind scheme can be written as (cf. Eq. (4.91))

\[\begin{aligned} \left(\vec{F}_c\right)_{I+1/2} = \frac{1}{2}\Big[&\vec{F}_c(\vec{W}_R) + \vec{F}_c(\vec{W}_L)\\ &- \left(\bar{\Gamma}\bar{P}^{-1}\right)_{I+1/2}\left|\bar{A}_{Roe}\right|_{I+1/2}\left(\vec{W}_R - \vec{W}_L\right)\Big]\,, \end{aligned} \tag{9.54}\]

其中\(|\bar{A}_{Roe}| = \bar{T}_{c,p}|\bar{\Lambda}_{c,p}|\bar{T}_{c,p}^{-1}\)。左特征向量（\(\bar{T}_{c,p}^{-1}\)）、右特征向量（\(\bar{T}_{c,p}\)）以及特征值（\(\bar{\Lambda}_{c,p}\)）都是矩阵\(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\)的。在\((I + 1/2)\)处的流动变量值由式(4.89)给出的Roe平均获得，并且整个残差乘以\(\bar{P}\bar{\Gamma}^{-1}\)。同样，式(9.54)的预条件Roe格式也可以用原始变量\(\vec{W}_p\)表述，这就得到en

where \(|\bar{A}_{Roe}| = \bar{T}_{c,p}|\bar{\Lambda}_{c,p}|\bar{T}_{c,p}^{-1}\). The left (\(\bar{T}_{c,p}^{-1}\)) and right (\(\bar{T}_{c,p}\)) eigenvectors, as well as the eigenvalues (\(\bar{\Lambda}_{c,p}\)) are those of the matrix \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\). Values of the flow variables at \((I + 1/2)\) are obtained by Roe's averaging given in Eq. (4.89), and the whole residual is multiplied by \(\bar{P}\bar{\Gamma}^{-1}\). Again, the preconditioned Roe scheme in Eq. (9.54) can be formulated in the primitive variables \(\vec{W}_p\) leading us to

\[\begin{aligned} \left(\vec{F}_c\right)_{I+1/2} = \frac{1}{2}\Big[&\vec{F}_c(\vec{W}_{p,R}) + \vec{F}_c(\vec{W}_{p,L})\\ &- \bar{\Gamma}_{I+1/2}\left|\bar{A}_{Roe,p}\right|_{I+1/2}\left(\vec{W}_{p,R} - \vec{W}_{p,L}\right)\Big]\,. \end{aligned} \tag{9.55}\]

此时，构成\(\bar{A}_{Roe,p}\)的特征值和特征向量由矩阵\(\Gamma^{-1}\bar{A}_{c,p}\)确定（参见式(9.47)）。\(\Gamma^{-1}\bar{A}_{c,p}\)的特征值与\(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\)的特征值相同，即等于\(\bar{\Lambda}_{c,p}\)，但特征向量不同。因此，式(9.55)中的Roe矩阵构成为\(|\bar{A}_{Roe,p}| = \bar{T}_p|\bar{\Lambda}_{c,p}|\bar{T}_p^{-1}\)。应当指出，在这种情形下，残差同样必须乘以\(\bar{P}\bar{\Gamma}^{-1}\)，以便回到守恒变量。en

The eigenvalues and eigenvectors which compose \(\bar{A}_{Roe,p}\) are now determined by the matrix \(\Gamma^{-1}\bar{A}_{c,p}\) (cf. Eq. (9.47)). The eigenvalues of \(\Gamma^{-1}\bar{A}_{c,p}\) are identical to the eigenvalues of \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\), i.e. to \(\bar{\Lambda}_{c,p}\), but the eigenvectors are different. Thus, the Roe matrix in Eq. (9.55) is composed as \(|\bar{A}_{Roe,p}| = \bar{T}_p|\bar{\Lambda}_{c,p}|\bar{T}_p^{-1}\). It should be mentioned that also in this case the residual has to be multiplied by \(\bar{P}\bar{\Gamma}^{-1}\) in order to obtain conservative variables.