3.2.2. 无粘通量与流场间断
(1) 双曲型守恒方程中的间断
(1.a) 广义解
通常，我们认为偏微分方程的解是连续可微的。 但是，一般的非线性双曲型守恒方程组 (Eq.3.2.1) 的解中都有可能存在间断， 因此，必须拓展双曲型守恒律解的概念。
以一维双曲型守恒方程组为例， 引入广义解 (generalized solution) 或弱解 (weak solution) 概念来表示包含间断的解：
设分片连续可微的函数 \(\mathbf{U}(x,t)\), 若任意与 \(\mathbf{U}\) 只有有限个交点的 光滑封闭曲线 \(\Gamma\) 都满足 \(\oint (\mathbf{F} dt - \mathbf{U} dx) = 0\), 则称 \(\mathbf{U}\) 是偏微分方程组在初值 \(\mathbf{U}(x,0)=\mathbf{U}_0(x)\) 下的广义解。 可以看出，在光滑区弱解就是通常的连续可微解。而在间断处，需要额外的条件来确定间断两侧的关系。
广义解推广了偏微分方程初值问题连续可微解的概念，但其后果是导致了广义解不唯一。 但是对于有明确物理意义的守恒律，其中只有一个解是有物理意义的，我们称之为 物理解 。 因此需要添加额外的条件来收敛到物理解，由于这个条件的作用与热力学第二定律相同， 因此被称为 “熵条件” (entropy condition)。
对于 CFD 而言，激波的两侧通过 R-H 关系相联系。
(1.b) 守恒格式和相容性
守恒格式 (conservative numerical scheme) 指格式有统一且守恒形式的通量 (with consistent and conservative numerical fluxes)。 即，在相邻单元 L 和 R 的界面上，从 L 指向 R 的通量等于从 R 指向 L 通量的相反数。 注意，前者使用 L 的模板点计算，后者使用 R 的模板点计算。
Lax-Wendroff theorem : a conservative numerical scheme for a hyperbolic system of conservation laws converges, then it converges towards a weak solution.
数值通量的相容条件：假设数值通量的计算函数为 \(\hat{\mathbf{F}}(\mathbf{U}_k | k= \cdots)\) (即该单元某个面的通量由若干个相应模板点计算得到), 那么, 如果满足 \(\hat{\mathbf{F}}(\mathbf{U}, \cdots, \mathbf{U}) = \mathbf{F}(\mathbf{U})\), 则称该数值通量有相容性。
结论: 相容的守恒格式 具有自动捕捉激波和接触间断的能力。
Euler方程 (和其它基于物理定律的守恒方程一样) 具有旋转不变性。 因此，可以将基于一维 Riemann 问题的通量计算格式推广到高维问题中。
(2) 一维迎风型有限差分格式
双曲型方程的差分格式收敛 (在Lax等价性定理满足时亦即稳定) 的必要条件 是差分格式的依赖域包含微分方程的依赖域, 即满足 CFL (Courant-Friedrichs-Lewy) 条件。
非线性双曲型方程组的的迎风格式本质上是线性格式形式上的推广。 最为常用的是 ROE 方法。
(2.a) 一维 Euler 方程
一维 Euler 方程 (理想气体) 的守恒形式为：
写为拟线性形式:
注意到 Euler 方程的通量 \(\mathbf{F}\) 是守恒变量 \(\mathbf{U}\) 的一次齐次函数:
那么，有
系数矩阵 \(A\) 的特征值和特征向量是:
记, \(\Lambda_{ij} = \lambda_i \delta_{ij}\), \(R = [\mathbf{v_1}, \mathbf{v_2}, \mathbf{v_3}]\)。 其中, \(L = R^{-1}\)。
对式 Eq.3.2.27 左乘 \(L\), 记 \(\mathbf{W}=L\mathbf{U}\) 为特征变量, 得：
形式上类似于线性波动方程组，适用于迎风格式。
(2.b) 通量微分分裂方法: ROE 格式
式 Eq.3.2.27 的守恒型差分格式为：
其中, \(\tilde{\mathbf{F}}_{j+1/2}\) 由 \(\mathbf{U}_k (k=\cdots)\) 计算得到。
一阶 Roe 格式通量:
其中, \(|\tilde A|=\tilde R |\tilde \Lambda| \tilde L\), \(|\tilde \Lambda| = \text{sign}(\tilde \lambda_i)\delta_{ij}\)。
将 \(\mathbf{U}\) 分解为特征振幅 \(\mathbf{U}_j=\sum_p^3 \beta_p \mathbf{v}_p\)。 那么, 界面处的守恒变量变化 (jump)为:
其中,
界面处的守恒变量 \(\tilde{\mathbf{U}}_{j+1/2}\) 由左右状态 \(\mathbf{U}_L = \mathbf{U}_{j}, \mathbf{U}_R = \mathbf{U}_{j+1}\) 平均得到 (Roe 平均):
根据线性 Reimann 问题的信息传播方向, 界面处的解为:
那么, 根据 \(\tilde A \tilde{\mathbf{v}}_p = \tilde{\lambda}_p \tilde{\mathbf{v}}_p\), 式 Eq.3.2.39 得
熵修正 : 经验上, 速度 \(u \approx a\) 时 Roe solver 会给出错误答案, 这是因为线性化在此时违背了热力学第二定律。因此, 当 \(|\tilde{\lambda}_p| < \epsilon \; (p=1,3)\) 时，对特征值进行修正:
以下为 ROE 格式通量在一维 Euler 问题中的实现，为了更清晰，代码的实现不考虑效率。
class Roe():
def __init__(self) -> None:
super().__init__()
@staticmethod
def average(rhoL, rhoR, uL, uR, tHL, tHR):
'''
Calculate ROE average rho, u, H, a
'''
ratio = np.sqrt(rhoR/rhoL)
roe_rho = np.sqrt(rhoL*rhoR)
roe_u = (uL + ratio*uR )/(1+ratio)
roe_H = (tHL + ratio*tHR)/(1+ratio)
roe_a = np.sqrt((GAMMA-1)*(roe_H-0.5*roe_u**2))
return roe_rho, roe_u, roe_H, roe_a
@staticmethod
def eigenvalue(u: float, a: float) -> np.ndarray:
'''
Calculate ROE eigenvalue[3] (with entropy fix)
'''
eps = 0.3*(a+abs(u))
eig = np.abs(np.array([u-a, u, u+a]))
for i in range(3):
if eig[i]<eps:
eig[i] = 0.5*(eig[i]**2/eps+eps)
return eig
@staticmethod
def eigenvector(u: float, a: float, tH: float) -> List[np.ndarray]:
'''
Calculate the eigenvectors v[3][3]
'''
vec = [None, None, None]
vec[0] = np.array([1.0, u-a, tH-u*a ])
vec[1] = np.array([1.0, u, 0.5*u**2])
vec[2] = np.array([1.0, u+a, tH+u*a ])
return vec
@staticmethod
def jump_weight(rhoL, rhoR, uL, uR, pL, pR, roe_rho, roe_a) -> np.ndarray:
'''
Calculate the weight (alpha[3]) of eigenvectors for
the jump on the face between the Left and Right values.
'''
dr = rhoR - rhoL
dp = pR - pL
du = uR - uL
alpha = np.zeros(3)
alpha[0] = 0.5*(dp-roe_a*roe_rho*du)/roe_a**2
alpha[1] = dr - dp/roe_a**2
alpha[2] = 0.5*(dp+roe_a*roe_rho*du)/roe_a**2
return alpha
@staticmethod
def flux_face(UL: np.ndarray, UR: np.ndarray) -> np.ndarray:
'''
Calculate ROE scheme flux on the face between Left and Right values.
'''
rhoL = UL[0]; uL = UL[1]/rhoL; pL=(GAMMA-1)*(UL[2]-0.5*rhoL*uL**2)
rhoR = UR[0]; uR = UR[1]/rhoR; pR=(GAMMA-1)*(UR[2]-0.5*rhoR*uR**2)
tHL = (UL[2] + pL)/rhoL
tHR = (UR[2] + pR)/rhoR
fL = np.array([UL[1], UL[1]*uL+pL, uL*(UL[2]+pL)])
fR = np.array([UR[1], UR[1]*uR+pR, uR*(UR[2]+pR)])
roe_rho, roe_u, roe_H, roe_a \
= Roe.average(rhoL, rhoR, uL, uR, tHL, tHR)
eigenvalue = Roe.eigenvalue (roe_u, roe_a)
eigenvector = Roe.eigenvector(roe_u, roe_a, roe_H)
alpha = Roe.jump_weight(rhoL, rhoR, uL, uR, pL, pR, roe_rho, roe_a)
fUpwind = np.zeros(3)
for i in range(3):
fUpwind += alpha[i]*eigenvalue[i]*eigenvector[i]
fFace = 0.5*(fL+fR) - 0.5*fUpwind
return fFace
(2.c) 矢通量分裂方法
矢通量分裂方法 (flux vector splitting, FVS) 最早由 Steger-Warming 提出。 将通量按特征值分解为
其中，\(|\lambda|_{\text{max}}\) 为特征值绝对值的最大值 (可以在该项前面乘上一个系数如 1.01 来确保分解后的矩阵特征值非负/非正)。 \(\nu_a\) 是针对物面法向的人工稳定因子，是为了弥补不能对角化的物面法向上粘性项的影响， 以增加数值计算的稳定性，仅在壁面附近的物面法向添加 (\(\hat{l}_{\text{cell}}\) 是网格空间尺度，如物面法向方向的长度)。
所以，式 Eq.3.2.38 中的数值通量为：
也就是说，考虑到波的传播方向，界面 \(j+1/2\) 左侧的 \(\mathbf{F}^{+}_{j}\) 和右侧的 \(\mathbf{F}^{-}_{j+1}\) 对界面的数值通量都有贡献。
与 ROE 格式的区别: ROE 取 \(j+1/2\) 处的特征结构, FVS 通量分裂逐点进行。
注意到，如果通量 \(\mathbf{F}\) 不是守恒变量 \(\mathbf{U}\) 的一次齐次函数， 则不能构造相应的 FVS 方法。另外，通量分裂的方法不是唯一的。
(3) 一维迎风型有限体积格式
记网格单元的体积平均 \(\bar{\mathbf{U}}_j = \int_{x_{j-1/2}}^{x_{j+1/2}} {\mathbf{U}(x)} dx\) 仍为 \(\mathbf{U}_j\), 即 Eq.3.2.38 形式不变。
首先对 \(\mathbf{U}(x)\) 进行重构，从而获得界面上的值。 以一阶重构 (空间二阶精度) 为例, 界面 \(j+1/2\) 的左右端值为:
那么, \(\tilde{\mathbf{F}}_{j+1/2}\) 的计算按照式 Eq.3.2.45 进行。 当然，在存在间断的流场中，直接使用一阶重构经常会出现数值发散 (密度或压力在重构时变为负数) 的情况。 因此，通常需要引入限制器或高分辨率格式。