CFD中的跨音速通量格式

超音速问题相对于不可压缩流存在一些问题。如果考虑流场各处都是超音速的情况，也即双曲方程。流动信息全部从上游向下游进行传输，相对简单。但对于低音速的区域，流动特征既可以向上游也可以向下游传递。在二者混合的情况下，理想的通量格式，应该满足超音速与低音速俩种情况下的物理特征。即对于一个网格面，在超音速的情况下，从上游网格单元插值。对于低音速的情况下，从网格单元的两侧来获取插值。

HLL格式

在引入HLL格式之前。简单介绍几个背景问题：

own网格与nei网格：在有限体积法中，可以将网格面毗连的网格单元定义为这个面的own单元以及nei单元。每个内部面均存在own单元与nei单元。

网格体心值：定义在网格单元点上的值。考虑变量为\(\mathbf{W}\)，网格体心值可以表示为\(\frac{1}{V}\int_{V}\mathbf{W}\rd V\)。

变量重组：这里区分为一阶以及高阶重组。考虑下图，在使用一阶格式的时候，网格内的各个位置的变量为均一的。考虑二阶格式，单一网格内不同位置的变量值不同，存在一个梯度。若使用更高阶的格式，图中的直线将会变成曲线。不管是一阶格式还是二阶格式，对于某个网格面上的值，均存在从不同的变量\(\mathbf{W}_{own}\)与\(\mathbf{W}_{nei}\)。

1) 网格值重组示意图

守恒法则：在下面的网格示意图中，黑线表示网格单元面，中间的黑线中的圆点表示面心，其中会存在黎曼问题。蓝色虚线包围的区域（W, evolve）表示影响区域，其中的\(\bfW\)的值会随着时间的推进而变化。对于时间步长比较小的情况下，影响区域仅仅存在于网格单元面的外围小区域。其他区域（W，constant）表示未受影响区域，其中的\(\bfW\)的值不会随着时间变化。现考虑红色虚线表示的网格，积分形式的守恒法则可以表示为：

(1)\[ \int_{V}\bfW^{t+\dt}\rd V= \int_{V}\bfW^t\rd V +\int_t^{t+\dt}F(\bfW_L)\rd t - \int_t^{t+\dt}F(\bfW_R)\rd t \]

方程(2)是精准的。同时，由于红色虚线位于未受影响区域，因此红色虚线上的通量并无变化。因此方程(2)可以简化为：

(2)\[ \int_{V}\bfW^{t+\dt}\rd V= \bfW^t V +\bigl(F(\bfW_L^t) - F(\bfW_R^t)\bigr)\dt \]

其中

(3)\[ \bfW^t V\approx\bfW_{own}^t (V_{own}+ V_{own}')+\bfW_{nei}^t (V_{nei}+ V_{nei}') \]

即

(4)\[ \int_{V}\bfW^{t+\dt}\rd V= \bfW_{own}^t (V_{own}+ V_{own}')+\bfW_{nei}^t (V_{nei}+ V_{nei}') +\bigl(F(\bfW_L^t) - F(\bfW_R^t)\bigr)\dt \]

2) 网格示意图

从另外的角度来看，\(\int_{V}\bfW^{t+\dt}\rd V\)可以分为三个区域：

(5)\[ \int_{V}\bfW^{t+\dt}\rd V= \int_{V_{own}}\bfW^{t+\dt}\rd V +\int_{V'}\bfW^{t+\dt}\rd V +\int_{V_{nei}}\bfW^{t+\dt}\rd V \]

再一次的，由于\(V_{own}\)以及\(V_{nei}\)区域为非影响区域，因此可以简化为：

(6)\[ \int_{V}\bfW^{t+\dt}\rd V= V_{own}\bfW^{t}_{own} +\int_{V'}\bfW^{t+\dt}\rd V +V_{nei}\bfW^{t}_{nei} \]

结合方程(6)与(4)，有：

(7)\[ \int_{V'}\bfW^{t+\dt}\rd V = \bfW_{own}^t V_{own}'+\bfW_{nei}^t V_{nei}' +\bigl(F(\bfW_L^t) - F(\bfW_R^t)\bigr)\dt \]

左右除以体积\(V'\)有：

(8)\[ \frac{1}{V'}\int_{V'}\bfW^{t+\dt}\rd V = \Bigl(\bfW_{own}^t V_{own}'+\bfW_{nei}^t V_{nei}' +\bigl(F(\bfW_L^t) - F(\bfW_R^t)\bigr)\dt\Bigr)\frac{1}{V'} \]

现在来看影响区域的体积计算方法。对于欧拉方程，存在三个不同的特征值，对于一个网格面，分别向own以及nei网格进行传播。定义nei网格单元在面处向nei网格传播的最快波速为\(a_{nei}^f\)，own网格单元在面处向own网格传播的最快波速为\(a_{own}^f\)。因此，有：

(9)\[ V_{nei}'\approx\Delta t|\bfS_f| a_{nei}^f, V_{own}'\approx\Delta t|\bfS_f| a_{own}^f \]

下一步考虑通量\(F\)，\(F\)为被传输变量\(\bfW\)的通量函数，其可以写为：

(10)\[ F(\bfW_R^t)\approx\phi_{nei}\bfW_{nei}^t,F(\bfW_L^t)\approx\phi_{own}\bfW_{own}^t \]

\(\frac{1}{V'}\int_{V'}\bfW^{t+\dt}\rd V\)可以理解为\(V'\)控制体的平均的\(\bfW\)的值。在HLL方法中，\(\bfW\)被认为是\(\bfW_{hll}\)。

(11)\[ \bfW_{hll} = \frac{\bfW_{own}^t |\bfS_f| a_{own}^f+\bfW_{nei}^t |\bfS_f| a_{nei}^f +\bigl(\phi_{own}\bfW_{own}^t - \phi_{nei}\bfW_{nei}^t\bigr)}{|\bfS_f| a_{nei}^f+|\bfS_f| a_{own}^f} \]

\(V'\)控制体的任一点的值均为\(\bfW_{hll}\)，因此面上的值为\(\bfW_{hll}\)。需要注意的是，方程(11)定义的是黎曼问题部分向左或者部分向右传递的情况。如果特征值均向左或者向右传递，HLL降低为普通插值格式。

矢通量分裂（vanLeer、AUSM）

在可压缩领域，AUSM格式与vanLeer格式等，都属于矢通量分裂格式。在引入vanLeer以及AUSM格式之前。简单介绍几个背景问题：

矢通量分裂： 对于网格单元面，可以认为其影响区域取决于上下游。在超音速情况下，网格面信息全部来源于上游。亚音速情况下，网格面信息来自于左右两侧的网格单元。因此矢通量分离将网格单元面的贡献依据特征值的大小，区分为左右网格单元（上下游）的贡献。矢通量分裂包含一系列的格式，比如Steger-Warning，vanLeer，以及AUSM等。

马赫数\(\bfM\)： \(\bfU/c\)，\(c\)表示音速。其中马赫数的分量表示速度各个方向的马赫数。

矢通量分裂采用马赫数来表示上游（迎风）方向。考虑\(x\)方向，\(M\geqslant 1\)的情况下，网格面的值全部来自于上游。反之，\(M\leqslant -1\)的情况下，网格面的值全部来自于下游。在\(-1<M<1\)的情况下，网格面的值来自于上下游的混合。AUSM分裂与vanLeer分裂针对马赫数的判断相同，均认为：

(12)\[ M^+=\frac{1}{4}(M+1)^2,M^-=-\frac{1}{4}(M-1)^2 \]

注意：

其中\(^+,^-\)表示上下游。具体是own网格还是nei网格，可以通过通量的大小来判断。

方程(12)为vanleer提出。其满足\(M\)为连续的、可导的、并且为低阶的。因此，矢通量分裂类格式在设计的情况下：

\(M_f\geqslant 1: F(\bfW)_f=F(\bfW^{+})_f\)；
\(M_f\leqslant -1: F(\bfW)_f=F(\bfW^{-})_f\)；
\(-1<M_f<1: F(\bfW)_f=F(\bfW^{+})_f+F(\bfW^{-})_f\)；

首先看一维欧拉方程的通量形式：

质量通量：\(F(\bfW_1)_f=\phi_f\rho_f c_f M_f\);
动量通量：\(F(\bfW_2)_f=\phi_f\rho_f c_f^2 (M_f^2+\frac{1}{\gamma_f}) \);
能量通量：\(F(\bfW_3)=\phi_f\rho_f c_f^3 M_f (\frac{1}{2}M_f^2+\frac{1}{\gamma_f-1}) \);

在使用vanLeer分裂对方程进行离散的情况下，当\(-1<M_f<1\)时，需要进一步的改写，例如质量通量：

(13)\[ F(\bfW_1)_f=\phi_f\rho_f c_f M_f=\phi_f\rho_f c_f\frac{1}{4}(M_f+1)^2-\phi_f\rho_f c_f\frac{1}{4}(-M_f+1)^2 \]

方程(13)右侧展开后与原始形式是相同的。同理，有动量通量：

(14)\[\begin{split} F(\bfW_2)_f=\phi_f\rho_f c_f \Bigl(\frac{1}{4}(M_f+1)^2\Bigr)\left((\gamma_f-1)M_fc_f+2c_f\right)\frac{1}{\gamma_f} \\ +\phi_f\rho_f c_f \Bigl(\frac{1}{4}(-M_f+1)^2\Bigr)\left((\gamma_f-1)M_fc_f+2c_f\right)\frac{1}{\gamma_f} \end{split}\]

以及能量通量：

(15)\[\begin{split} F(\bfW_3)_f=\phi_f\rho_f c_f \Bigl(\frac{1}{4}(M_f+1)^2\Bigr)\left((\gamma_f-1)M_fc_f+2c_f\right)^2\frac{1}{2\gamma_f^2-2} \\ +\phi_f\rho_f c_f \Bigl(\frac{1}{4}(-M_f+1)^2\Bigr)\left((\gamma_f-1)M_fc_f+2c_f\right)^2\frac{1}{2\gamma_f^2-2} \end{split}\]

在考虑通量分裂方向的情况下，方程(13)可以写为

(16)\[\begin{split} F(\bfW_1)_f=\phi_f^+\rho_f^+c_f^+M_f^+&+\phi_f^-\rho_f^-c_f^-M_f^- = \\ &\Bigl(\phi_f^+\rho_f^{+}c_f^{+}\frac{1}{4}(M_f^{+}+1)^2\Bigr)+\Bigl(-\phi_f^-\rho_f^{-}c_f^{-}\frac{1}{4}(M_f^{-}-1)^2\Bigr) \end{split}\]

方程(16)即为连续性方程采用vanLeer通量矢量分裂的计算方法。其他的动量方程以及能量方程，如方程(14)、(15)类似，在此不详细写出。

AUSM格式与vanLeer格式的主要区别，在于AUSM格式认为方程变量中的对流与压力来自于不同的贡献，因此在处理通量的情况下，需要将对流贡献与压力贡献单独处理。将压力贡献提取出来之后，回看一维欧拉方程的通量形式：

质量通量对流贡献：\(F(\bfW_1)_f=\phi_f\rho_f c_f M_f\);
动量通量的对流贡献以及压力贡献之和：\(F(\bfW_2)_f=\phi_f\rho_f c_f^2 M_f^2+\phi_fp \);
能量通量对流贡献：\(F(\bfW_3)=\phi_f\rho_f c_f M_f H_f \);

其中的对流贡献与vanLeer格式相同。对于压力贡献，在\(-1<M<1\)的情况下，参考方程(12)，AUSM格式将其写为：

(17)\[ p^+=\frac{1}{2}p(1+M),p^-=\frac{1}{2}p(1-M) \]

方程(18)也可以拓展为二阶：

(18)\[ p^+=\frac{1}{4}p(M+1)^2(2-M),p^-=\frac{1}{4}p(M-1)^2(2+M) \]

因此，动量通量的压力贡献，当\(-1<M_f<1\)时：

(19)\[ \phi_fp=\phi_f^+p^+_f+\phi_f^-p^-_f= \frac{1}{2}\phi_f^+p^+_f(1+M^+_f) + \frac{1}{2}\phi_f^-p^-_f(1-M^-_f) \]

通量差分裂（Roe）

Roe格式与前述格式的思想均不相同。Roe格式的本质在于将原本非线性的双曲偏微分方程转化为线性偏微分方程组。这一套线性的偏微分方程组可以进一步的获得精确解。那么Roe方法的误差就主要来源于将原始偏微分方程线性化的这个过程中。在Roe方法的提出过程中，线性方程组的雅克比矩阵需要满足若干条件。在这里原始的Roe方法推导过程以及设计思想不会进行详细描述。仅仅介绍Roe算法的计算流程。

针对不同的非线性偏微分方程组，Roe方法需要设计不同的线性偏微分方程组。在CFD领域，Roe方法主要用于求解欧拉方程。