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}\)。

注意：

\(\bfW\)表示欧拉方程的守恒变量，例如\(\bfW_0=\rho\)，\(\bfW_1=\rho\bfU\)，\(\bfW_2=\rho E\)

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降低为普通插值格式。

AUSM格式

AUSM格式依据面网格的马赫数方向，来对通量进行区分。在可压缩领域，存在大量的矢通量格式、通量差分裂等格式。AUSM格式将压力与对流速度进行区分，理解相对更加简单。首先介绍几个背景概念：

音速\(c\)：\(c=\sqrt{\gamma\frac{p}{\rho}}\)。马赫数标量\(|\bfM|\)：\(\frac{|\bfU|}{c}\)。

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

面马赫数\(M_f\)： \(\frac{\bfU_f\cdot\bfS_f}{c |\bfS_f|}\)，其为一个定义在网格面上的标量。其值的大小可以用来判断插值的方向。

首先将通量均通过马赫数来表示，其中面的速度，可以写成网格面马赫数的形式：

(12)\[ \bfU_f \cdot \bfS_f= c_f M_f |\bfS_f| \]

进一步的有质量密度、动量密度、能量密度通量的对流项：

(13)\[ \bfW_f (\bfU_f \cdot \bfS_f)=\bfW_f(c_f M_f |\bfS_f|) =\bfW_f c_f |\bfS_f| M_f \]

动量密度\(\bfW_1\)通量的压力项：

(14)\[ p_f\bfS_f=\frac{c_f^2}{\gamma_f}\rho_f\bfS_f \]

能量密度\(\bfW_2\)通量压力对流项：

(15)\[ p_f\bfU_f \cdot \bfS_f= \frac{c_f^2}{\gamma_f}\rho_f(c_f M_f |\bfS_f|) \]

AUSM分裂认为对于网格单元面，可以认为其影响区域取决于上下游。在超音速情况下，网格面信息全部来源于上游。亚音速情况下，网格面信息来自于左右两侧的网格单元。因此AUSM分裂将网格单元面的贡献依据特征值的大小，区分为左右网格单元（上下游）的贡献：

\(M_f\geqslant 1\)：网格面的值全部来自于上游。
\(M_f\leqslant -1\)：网格面的值全部来自于下游。
\(-1<M_f<1\)：网格面的值来自于上下游的混合。

AUSM分裂认为：

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

(17)\[ M_f=M_f^+ + M_f^- \]

注意：

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

方程(16)满足\(M\)为连续的、可导的、并且为低阶的。\(M_f\)用于判断信息的方向并进行插值。至此，AUSM中的关于速度通量的这部分的计算过程如下：

通过方程(17)来计算面马赫数；
依据面马赫数定义的方向，来更新面守恒变量的值\(\bfW_f\)，\(c\)为定值；至此，方程(13)中的\(\bfW_f c_f |\bfS_f|\)可以更新；
结合方程(17)来更新方程(13)的\(\bfW_f c_f |\bfS_f| M_f \)；

下面来看压力的贡献，AUSM格式认为方程变量中的对流与压力来自于不同的贡献，因此在处理通量的情况下，需要将对流贡献与压力贡献单独处理。AUSM格式将其写为：

(18)\[ p^+_f=\frac{1}{2}\frac{c_f^2}{\gamma_f}\rho_f^+(1+M^+),p^-_f=\frac{1}{2}\frac{c_f^2}{\gamma_f}\rho_f^-(1-M^-) \]

(19)\[ p_f = \frac{c_f^2}{\gamma_f}\rho_f = p^+_f + p^-_f \]

至此，AUSM中的关于压力通量的这部分的计算过程如下：

通过方程(19)来计算\(p_f\)；
依据\(p_f\)方向，来更新方程(14)以及(15)；

近似黎曼求解器、Roe方法

Roe格式与前述格式的思想均不相同。Roe格式的本质在于将原本非线性的双曲偏微分方程转化为线性偏微分方程组。这是因为针对线性双曲系统，其特征值与特征变量是不变的。因此可以通过解耦特征变量的方法来进行推进求解。对于实际情况下的非线性双曲系统，特征向量特征值等都是变化的量会导致求解非常复杂。Roe格式作为近似黎曼求解器，就是把非线性系统转变为线性系统后，依据上文讨论的方法进行求解。在Roe方法中，原本变化的\(\alpha_i,\mathbf{r}_i,\lambda_i\)都会转变成不取决于网格左右变量的固定值后，通过特征变量的推进来进行。

Roe方法的详细讨论请参考《无痛苦NS方程笔记》。