5.3.1 Central Schemes with Artificial Dissipation 带人工耗散的中心格式[cfd-5-3-1]

中心格式的基本思想ï¼是按照方程(5.20)ï¼用控制体面两侧守恒变量的算术平均来计算该面上的对流通量。这会导致解的奇偶失联(odd-even decouplingï¼即离散方程产生两个相互独立的解)以及激波处的波动ï¼因此为了保持稳定性必须加入人工耗散。人工耗散基于二阶与四阶差分的混合。该格式由Jameson等人[43]最先在结构网格上针对欧拉方程实现。以各位作者姓氏的首字母命名ï¼它也简称为JST格式。en

The basic idea of the central scheme is to compute the convective fluxes at a face of the control volume from the arithmetic average of the conservative variables on both sides of the face according to Eq. (5.20). Since this would lead to odd-even decoupling of the solution (generation of two independent solutions of the discretised equations) and wiggles at shocks, artificial dissipation has to be added for stability. The artificial dissipation is based on a blend of second- and fourth-order differences. The scheme was first implemented for the Euler equations on structured grids by Jameson et al. [43]. Because of the names of the authors, it is also abbreviated as the JST scheme.

在非结构网格上实现JST格式时ï¼二阶差分采用拉普拉斯(Laplacian)算子ï¼四阶差分采用拉普拉斯的拉普拉斯[44]、[27]。为了降低计算成本ï¼采用伪拉普拉斯算子(pseudo-Laplacian)代替真正的拉普拉斯算子。为此,[45]首先提出了二维格式ï¼随后[46]作了改进。后来,[2]把该格式推广到三维。它利用了一种距离加权过程。这样ï¼对于任何网格上线性变化的函数ï¼该格式的伪拉普拉斯算子都为零。应用于单元 $$I$$ 中的某一般标量 $$U$$ 时ï¼伪拉普拉斯算子的形式为en

The implementation of the JST scheme on unstructured grids utilises the Laplacian operator for the second-order differences and the Laplacian of Laplacian for the fourth-order differences [44], [27]. In order to reduce the computational cost, pseudo-Laplacians are employed instead of true Laplacians. For this purpose, a 2-D formulation was proposed first in [45] and then improved in [46]. Later on, the scheme was extended to 3D in [2]. It makes use of a distance-weighting procedure. In this way, the scheme leads to a vanishing pseudo-Laplacian for a linearly varying function on any grid. Applied to a general scalar quantity $$U$$ in cell $$I$$ , the pseudo-Laplacian assumes the form

L(U_I) = \sum_{J=1}^{N_A} \theta_{IJ}\left(U_J - U_I\right), \tag{5.24}

其中 $$N_A$$ 表示相邻控制体的数目。在采用中点对偶格式(median-dual scheme)时ï¼单元指标须换成节点指标 $$(i,j)$$ 。式(5.24)中的求和最好像通量计算那样ï¼用对面(单元中心格式)或对边(中点对偶格式)的循环来计算。几何权重 $\theta$ 定义为en

where $$N_A$$ stands for the number of adjacent control volumes. The cell indices have to be substituted by node indices $$(i,j)$$ in the case of the median-dual scheme. The sum in Eq. (5.24) is best evaluated using either a loop over faces (cell-centred scheme) or a loop over edges (median-dual scheme) similar to the flux computation. The geometrical weights $\theta$ are defined as

\theta_{IJ} = 1 + \Delta\theta_{IJ} \tag{5.25}

权重由一个优化问题的解得到[2]。该优化问题借助拉格朗日乘子求解。由此ï¼几何权重由下式给出en

and result from the solution of an optimisation problem [2]. The optimisation problem is solved by means of Lagrange multipliers. Herewith, the geometrical weights are obtained from the expression

\Delta\theta_{IJ} = \lambda_{x,I}\left(x_J - x_I\right) + \lambda_{y,I}\left(y_J - y_I\right) + \lambda_{z,I}\left(z_J - z_I\right), \tag{5.26}

其中 $$x$$ 、 $$y$$ 、 $$z$$ 为单元形心的笛卡尔坐标(中点对偶格式则为节点坐标)。拉格朗日乘子 $\lambda$ 对每个单元(节点)计算ï¼由[2]得en

where $$x,y,z$$ are the Cartesian coordinates of the cell centroids (nodes in the case of the median-dual scheme). The Lagrange multipliers $\lambda$ are computed for each cell (node) and follow from [2]

\begin{aligned} \lambda_x &= \frac{R_x\,a_{11} + R_y\,a_{12} + R_z\,a_{13}}{d}\\ \lambda_y &= \frac{R_x\,a_{21} + R_y\,a_{22} + R_z\,a_{23}}{d}\\ \lambda_z &= \frac{R_x\,a_{31} + R_y\,a_{32} + R_z\,a_{33}}{d} \end{aligned} \tag{5.27}

其中的系数为en

with the coefficients

\begin{aligned} a_{11} &= I_{yy}I_{zz} - I_{yz}^{2}\\ a_{12} &= I_{xz}I_{yz} - I_{xy}I_{zz}\\ a_{13} &= I_{xy}I_{yz} - I_{xz}I_{yy}\\ a_{21} &= I_{xz}I_{yz} - I_{xy}I_{zz}\\ a_{22} &= I_{xx}I_{zz} - I_{xz}^{2}\\ a_{23} &= I_{xy}I_{xz} - I_{xx}I_{yz}\\ a_{31} &= I_{xy}I_{yz} - I_{xz}I_{yy}\\ a_{32} &= I_{xz}I_{xy} - I_{xx}I_{yz}\\ a_{33} &= I_{xx}I_{yy} - I_{xy}^{2}\\ d &= I_{xx}I_{yy}I_{zz} - I_{xx}I_{yz}^{2} - I_{yy}I_{xz}^{2} - I_{zz}I_{xy}^{2} + 2I_{xy}I_{xz}I_{yz}\,. \end{aligned} \tag{5.28}

对单元 $$I$$ 写出ï¼一阶矩为en

Written for a cell $$I$$ , the first-order moments read

\begin{aligned} R_{x,I} &= \sum_{J=1}^{N_A}\left(x_J - x_I\right)\\ R_{y,I} &= \sum_{J=1}^{N_A}\left(y_J - y_I\right)\\ R_{z,I} &= \sum_{J=1}^{N_A}\left(z_J - z_I\right). \end{aligned} \tag{5.29}

此外ï¼二阶矩由下式给出en

Furthermore, the second-order moments are given by

\begin{aligned} I_{xx,I} &= \sum_{J=1}^{N_A}\left(x_J - x_I\right)^{2}\\ I_{yy,I} &= \sum_{J=1}^{N_A}\left(y_J - y_I\right)^{2}\\ I_{zz,I} &= \sum_{J=1}^{N_A}\left(z_J - z_I\right)^{2}\\ I_{xy,I} &= \sum_{J=1}^{N_A}\left(x_J - x_I\right)\left(y_J - y_I\right)\\ I_{xz,I} &= \sum_{J=1}^{N_A}\left(x_J - x_I\right)\left(z_J - z_I\right)\\ I_{yz,I} &= \sum_{J=1}^{N_A}\left(y_J - y_I\right)\left(z_J - z_I\right). \end{aligned} \tag{5.30}

几何权重(5.25)在严重扭曲的网格上可能导致拉普拉斯算子的非正近似ï¼从而使稳定性丧失。因此,[45]建议把权重限制在 $$(0,2)$$ 范围内。不过ï¼这一措施会损害离散的精度。更多细节另见文献[12]中的讨论。en

The geometrical weights (5.25) can lead to a non-positive approximation of the Laplacian and hence to a lost of stability on severely distorted grids. Therefore, clipping the weights to the range $$(0,2)$$ was suggested in [45]. However, this measure impairs the accuracy of the discretisation. See also the discussion in Ref. [12] for further details.

四阶差分用拉普拉斯的拉普拉斯来计算ï¼即在式(5.24)中用 $$L(U)$$ 代替 $$U$$ 。于是ï¼单元 $$I$$ 的人工耗散项的最终形式为en

The fourth-order differences are evaluated as the Laplacian of the Laplacian, i.e., $$L(U)$$ is substituted for $$U$$ in Eq. (5.24). Hence, the final form of the artificial dissipation term for a cell $$I$$ is

\begin{aligned} \vec{D}_I &= \sum_{J=1}^{N_A}(\hat{\Lambda}_c)_{IJ}\,\epsilon_{IJ}^{(2)}\,\theta_{IJ}\left(\vec{W}_J - \vec{W}_I\right)\\ &-\sum_{J=1}^{N_A}(\hat{\Lambda}_c)_{IJ}\,\epsilon_{IJ}^{(4)}\,\theta_{IJ}\left[L(\vec{W}_J) - L(\vec{W}_I)\right]\,. \end{aligned} \tag{5.31}

加入人工耗散项后ï¼方程(5.2)中的方程组变为en

With the artificial dissipation term added, the system of equations in Eq. (5.2) becomes

\Omega_I\,\frac{d\vec{W}_I}{dt} = -\left[\sum_{m=1}^{N_F}\left(\vec{F}_c - \vec{F}_v\right)_m\Delta S_m\right] + \vec{D}_I + \vec{Q}_I\,\Omega_I\,, \tag{5.32}

其中 $$N_F$$ 表示控制体的面数(它可能与相邻控制体的数目不同ï¼例如当一个四边形面被分成两个三角形时)。en

where $$N_F$$ denotes the number of the faces of the control volume (which may differ from the number of adjacent control volumes, e.g., if a quadrilateral face is divided into two triangles).

式(5.31)中的二阶项与四阶项用对流通量雅可比矩阵的谱半径来缩放。按照文献[28]ï¼单元 $$I$$ 的谱半径可按下式计算en

The second- and the fourth-order terms in Eq. (5.31) are scaled by the spectral radius of the convective flux Jacobian. According to Ref. [28], the spectral radius for cell $$I$$ can be evaluated as

(\hat{\Lambda}_c)_I = \sum_{m=1}^{N_F}\left(\left|V_m\right| + c_m\right)\Delta S_m\,, \tag{5.33}

其中 $$V_m$$ 表示逆变速度(2.22), $$c_m$$ 为声速。这两个量都用面上平均的流动变量计算。控制体面上的谱半径由下式得到en

where $$V_m$$ represents the contravariant velocity (2.22) and $$c_m$$ the speed of sound, respectively. Both quantities are computed from flow variables averaged at the face. The spectral radius at the face of the control volume is obtained from

(\hat{\Lambda}_c)_{IJ} = \frac{1}{2}\left[(\hat{\Lambda}_c)_I + (\hat{\Lambda}_c)_J\right]\,. \tag{5.34}

利用一个基于压力的传感器ï¼在激波处关闭四阶差分ï¼在流场的光滑区域关闭二阶差分。据此ï¼式(5.31)中的系数 $\epsilon_{IJ}^{(2)}$ 与 $\epsilon_{IJ}^{(4)}$ 定义为en

A pressure-based sensor is used to switch off the fourth-order differences at shocks and the second-order differences in smooth portions of the flow field. Herewith, the coefficients $\epsilon_{IJ}^{(2)}$ and $\epsilon_{IJ}^{(4)}$ in Eq. (5.31) are defined as

\begin{aligned} \epsilon_{IJ}^{(2)} &= k^{(2)}\max(\Upsilon_I, \Upsilon_J)\\ \epsilon_{IJ}^{(4)} &= \max\left[0,\,\left(k^{(4)} - \epsilon_{IJ}^{(2)}\right)\right] \end{aligned} \tag{5.35}

其中的压力传感器由下式给出en

with the pressure sensor given by

\Upsilon_I = \frac{\left|\displaystyle\sum_{J=1}^{N_A}\theta_{IJ}\left(p_J - p_I\right)\right|}{\displaystyle\sum_{J=1}^{N_A}\left(p_J + p_I\right)} \tag{5.36}

参数的典型取值为 $k^{(2)} = 1/2$ ï¼以及 $1/128 \le k^{(4)} \le 1/64$ 。en

Typical values of the parameters are $k^{(2)} = 1/2$ and $1/128 \le k^{(4)} \le 1/64$ .

正如我们在4.3.1小节中已经讨论过的ï¼若在式(5.31)中用一个矩阵[47]代替谱半径 $(\hat{\Lambda}_c)_{IJ}$ ï¼上述中心格式的精度可以得到改进。这种所谓的矩阵耗散(matrix dissipation)格式在非结构网格上的实现方式与结构网格相同ï¼缩放矩阵的定义如同式(4.59)。矩阵耗散格式在三维混合网格上的应用可参见例如文献[48]。en

As we already discussed in Subsection 4.3.1, the accuracy of the above central scheme can be improved when we substitute a matrix [47] for the spectral radius $(\hat{\Lambda}_c)_{IJ}$ in Eq. (5.31). The implementation of this so-called matrix dissipation scheme on unstructured grids proceeds in the same way as on structured grids, with the scaling matrix defined as in Eq. (4.59). Application of the matrix dissipation scheme to 3-D mixed grids is discussed, e.g., in Ref. [48].

需要特别注意的是ï¼对于三角形/四面体以外的单元ï¼常用的显式Runge-Kutta型时间离散在与中心格式耦合时会出现严重的稳定性问题[49]。原因在于四阶差分是用拉普拉斯的拉普拉斯来表示的。一种补救办法是用左、右状态之差(参见4.3节)来近似四阶差分[49]ï¼即en

It is important to note that for elements other than triangles/tetrahedra, the popular explicit Runge-Kutta type of temporal discretisation experiences severe stability problems when it is coupled to the central scheme [49]. The reason is the representation of the fourth-order differences by the Laplacian of the Laplacian. A remedy is to employ a difference of the left and the right state (cf. Section 4.3) for the approximation of the fourth-order differences [49], i.e.,

\vec{D}_I = \sum_{J=1}^{N_A}(\hat{\Lambda}_c)_{IJ}\,\epsilon_{IJ}^{(2)}\,\theta_{IJ}\left(\vec{W}_J - \vec{W}_I\right) - \sum_{J=1}^{N_A}4\,(\hat{\Lambda}_c)_{IJ}\,\epsilon_{IJ}^{(4)}\left(\vec{W}_L - \vec{W}_R\right)\,. \tag{5.37}

这种做法在四边形/六面体网格上给出与相应结构格式相同的模板。左、右状态可用例如下文5.3.3小节所述的线性重构来计算。en

This approach leads on quadrilateral/hexahedral grids to the same stencil as the corresponding structured scheme. The left and right state are computed using, e.g., the linear reconstruction described below in Subsection 5.3.3.