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）

待更新。