膨胀扇的前缘和后缘分别对应上游和下游状态的 Mach 角。

上游 Mach 角：

\[ \mu_1=\sin^{-1}\frac{1}{M_1} \]

下游 Mach 角：

\[ \mu_2=\sin^{-1}\frac{1}{M_2} \]

因为膨胀后 \(M_2>M_1\)，所以 \(\mu_2<\mu_1\)。

对量热完全气体，Prandtl–Meyer 函数为：

\[ f(M)= \sqrt{\frac{\gamma+1}{\gamma-1}} \tan^{-1} \sqrt{\frac{\gamma-1}{\gamma+1}(M^2-1)} -\tan^{-1}\sqrt{M^2-1} \]

膨胀波前后满足：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ f_2=f_1+\theta } $} \]

计算流程：

\[ M_1\rightarrow f_1\rightarrow f_2=f_1+\theta\rightarrow M_2 \]

其中 \(M_2\) 需要由 Prandtl–Meyer 函数反求得到。

得到 \(M_2\) 后，可由等熵关系计算波后热力学状态。

压力：

\[ p_2=p_1 \left( \frac{ 1+\frac{\gamma-1}{2}M_1^2 }{ 1+\frac{\gamma-1}{2}M_2^2 } \right)^{\gamma/(\gamma-1)} \]

温度：

\[ T_2=T_1 \left( \frac{ 1+\frac{\gamma-1}{2}M_1^2 }{ 1+\frac{\gamma-1}{2}M_2^2 } \right) \]

密度由状态方程给出：

\[ \rho_2=\frac{p_2}{RT_2} \]

本章数值计算需正确捕捉：

• 膨胀波扇区
• 流动转向
• 马赫数增加
• 压强、温度、密度降低
• 空间推进过程中的边界影响

二维定常无粘 Euler 方程可写成强守恒形式：

\[ \frac{\partial F}{\partial x} =J-\frac{\partial G}{\partial y} \]

本问题为等熵、绝热、无体力流动，因此：

\[ J=0 \]

于是控制方程化为：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial F}{\partial x} =-\frac{\partial G}{\partial y} } $} \]

这里的关键是：
方程左端只含一个 \(x\) 方向导数，因此可以沿 \(x\) 方向进行下游推进。

\(x\) 方向通量向量为：

\[ F= \begin{Bmatrix} \rho u\\ \rho u^2+p\\ \rho uv\\ \rho u\left(e+\frac{V^2}{2}\right)+pu \end{Bmatrix} \]

其中：

\[ V^2=u^2+v^2 \]

将其元素记为：

\[ F_1=\rho u, \qquad F_2=\rho u^2+p \]

\[ F_3=\rho uv,\qquad F_4=\rho u\left(e+\frac{u^2+v^2}{2}\right)+pu \]

\(y\) 方向通量向量为：

\[ G= \begin{Bmatrix} \rho v\\ \rho uv\\ \rho v^2+p\\ \rho v\left(e+\frac{V^2}{2}\right)+pv \end{Bmatrix} \]

将其元素记为：

\[ G_1=\rho v \]

\[ G_2=\rho uv \]

\[ G_3=\rho v^2+p \]

\[ G_4=\rho v\left(e+\frac{u^2+v^2}{2}\right)+pv \]

对量热完全气体：

\[ e=c_vT=\frac{RT}{\gamma-1} =\frac{1}{\gamma-1}\frac{p}{\rho} \]

因此，\(F_4\) 可写为：

\[ F_4= \frac{\gamma}{\gamma-1}pu +\rho u\frac{u^2+v^2}{2} \]

类似地：

\[ G_4= \frac{\gamma}{\gamma-1}pv +\rho v\frac{u^2+v^2}{2} \]

这样处理的目的是什么？

数值推进直接得到的是：

\[ F_1,\quad F_2,\quad F_3,\quad F_4 \]

而真正需要分析的物理变量是：

\[ \rho,\quad u,\quad v,\quad p,\quad T \]

密度由一个二次方程得到：

\[ \rho= \frac{-B+\sqrt{B^2-4AC}}{2A} \]

其中：

\[ A=\frac{F_3^2}{2F_1^2}-F_4, \qquad B=\frac{\gamma}{\gamma-1}F_1F_2, \qquad C=-\frac{\gamma+1}{2(\gamma-1)}F_1^3 \]

密度 \(\rho\) 得到后，其余变量为：

\[ u=\frac{F_1}{\rho} \]

\[ v=\frac{F_3}{F_1} \]

\[ p=F_2-F_1u \]

\[ T=\frac{p}{\rho R} \]

为了计算：

\[ \frac{\partial F}{\partial x} =-\frac{\partial G}{\partial y} \]

需要在每个推进位置计算 \(G\)。

一种方法是：

\[ F\rightarrow \rho,u,v,p,T\rightarrow G \]

但可以采用守恒处理方式：

直接由 \(F_1,F_2,F_3,F_4\) 构造 \(G_1,G_2,G_3,G_4\)。

由定义可得：

\[ G_1=\rho v =\rho\frac{F_3}{F_1} \]

又因为：

\[ G_2=\rho uv=F_3 \]

对于 \(G_3\)：

\[ G_3=\rho v^2+p =\rho\left(\frac{F_3}{F_1}\right)^2+p \]

其中：

\[ p=F_2-\frac{F_1^2}{\rho}\quad \Longrightarrow \quad G_3= \rho\left(\frac{F_3}{F_1}\right)^2 +F_2-\frac{F_1^2}{\rho} \]

对于能量通量：

\[ G_4= \frac{\gamma}{\gamma-1}pv +\rho v\frac{u^2+v^2}{2} \]

代入：

\[ u=\frac{F_1}{\rho},\qquad v=\frac{F_3}{F_1},\qquad p=F_2-\frac{F_1^2}{\rho} \]

得到：

\[ G_4= \frac{\gamma}{\gamma-1} \left(F_2-\frac{F_1^2}{\rho}\right) \frac{F_3}{F_1} +\frac{\rho F_3}{2F_1} \left[ \left(\frac{F_1}{\rho}\right)^2 +\left(\frac{F_3}{F_1}\right)^2 \right] \]

这一组关系使得数值推进过程始终保持在守恒变量框架内。

采用如下简单贴体变换：

\[ \xi=x \]

\[ \eta=\frac{y-y_s(x)}{h(x)} \]

其中：

• \(y_s(x)\)：固壁下边界位置
• \(h(x)\)：从下壁面到上边界的局部高度

该变换保证：

\[ y=y_s(x)\Rightarrow \eta=0 \]

\[ y=y_s(x)+h(x)\Rightarrow \eta=1 \]

由链式法则：

\[ \frac{\partial}{\partial x}= \frac{\partial}{\partial \xi} \frac{\partial \xi}{\partial x}+ \frac{\partial}{\partial \eta} \frac{\partial \eta}{\partial x} \]

\[ \frac{\partial}{\partial y}= \frac{\partial}{\partial \xi} \frac{\partial \xi}{\partial y}+ \frac{\partial}{\partial \eta} \frac{\partial \eta}{\partial y} \]

由 \(\xi=x\) 可得：

\[ \frac{\partial \xi}{\partial x}=1,\qquad \frac{\partial \xi}{\partial y}=0 \]

由 \(\eta=(y-y_s)/h\) 可得：

\[ \frac{\partial \eta}{\partial y}=\frac{1}{h} \]

由：

\[ \eta=\frac{y-y_s(x)}{h(x)} \]

在固定 \(y\) 下对 \(x\) 求导：

\[ \frac{\partial\eta}{\partial x}= -\frac{1}{h}\frac{dy_s}{dx} -\frac{\eta}{h}\frac{dh}{dx} \]

膨胀角位置记为 \(x=E\)。

对 \(x\le E\)：

\[ y_s=0,\qquad h=\text{constant}\qquad \Longrightarrow \qquad \frac{\partial\eta}{\partial x}=0 \]

对 \(x\ge E\)：

\[ y_s=-(x-E)\tan\theta \]

\[ h=H+(x-E)\tan\theta \]

因此：

\[ \frac{dy_s}{dx}=-\tan\theta,\qquad \frac{dh}{dx}=\tan\theta \qquad \Longrightarrow \qquad \frac{\partial\eta}{\partial x}= \frac{(1-\eta)\tan\theta}{h} \]

所以：

\[ \frac{\partial\eta}{\partial x}= \begin{cases} 0,&x\le E\\[4pt] \dfrac{(1-\eta)\tan\theta}{h},&x\ge E \end{cases} \]

最终得到：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial}{\partial x}= \frac{\partial}{\partial \xi}+ \frac{\partial\eta}{\partial x} \frac{\partial}{\partial\eta} } $} \]

以及：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial}{\partial y}= \frac{1}{h} \frac{\partial}{\partial\eta} } $} \]

这一步的意义是：
把物理平面中的偏导数，转化为计算平面规则网格上的偏导数。

原方程为：

\[ \frac{\partial F}{\partial x} =-\frac{\partial G}{\partial y} \]

代入导数变换：

\[ \frac{\partial F}{\partial \xi}+ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F}{\partial\eta}= -\frac{1}{h} \frac{\partial G}{\partial\eta} \]

整理得到：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \frac{\partial F}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F}{\partial\eta}+ \frac{1}{h} \frac{\partial G}{\partial\eta} \right] } $} \]

这就是下游推进所用的计算平面控制方程。

写成分量形式：

\[ \frac{\partial F_1}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F_1}{\partial\eta} + \frac{1}{h}\frac{\partial G_1}{\partial\eta} \right] \]

\[ \frac{\partial F_2}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F_2}{\partial\eta} + \frac{1}{h}\frac{\partial G_2}{\partial\eta} \right] \]

\[ \frac{\partial F_3}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F_3}{\partial\eta} + \frac{1}{h}\frac{\partial G_3}{\partial\eta} \right] \]

\[ \frac{\partial F_4}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F_4}{\partial\eta} + \frac{1}{h}\frac{\partial G_4}{\partial\eta} \right] \]

我们在本算例中直接使用有量纲变量。

使用有量纲变量并不会影响 CFD 求解的正确性。
关键是所有物理量必须采用一致单位。

本章采用 SI 单位制：

• 长度：m
• 密度：kg/m\(^3\)
• 压强：N/m\(^2\)
• 温度：K
• 速度：m/s

这样可以直接获得具有工程量级感的计算结果。

对该几何，局部高度 \(h=h(x)\) 为：

\[ h= \begin{cases} 40\,\mathrm{m},&0\le x\le 10\,\mathrm{m}\\[4pt] 40+(x-10)\tan\theta,&10\le x\le 65\,\mathrm{m} \end{cases} \]

该函数用于计算：

\[ \frac{\partial\eta}{\partial x}= \begin{cases} 0,&x\le E\\[4pt] \dfrac{(1-\eta)\tan\theta}{h},&x\ge E \end{cases} \]

几何变化通过 \(h(x)\) 和 \(\partial\eta/\partial x\) 进入控制方程。

初始数据线取在：

\[ x=0 \]

在这条竖直线上，每个网格点都给定相同的均匀来流条件。

我们将 \(x=0\) 的初始数据线划分为：

\[ 41 \]

个网格点，即：

\[ j=1,2,\dots,41 \]

这意味着数值计算从一条已知的竖直初始线出发，
沿 \(x\) 方向一步一步向下游推进。

在 \(x=0\) 的所有 \(j\) 点上：

\[ u=0.678\times10^3\,\mathrm{m/s} \]

\[ v=0 \]

\[ \rho=0.123\times10^1\,\mathrm{kg/m^3} \]

\[ p=0.101\times10^6\,\mathrm{N/m^2} \]

\[ T=0.286\times10^3\,\mathrm{K} \]

\[ M=0.200\times10^1 \]

也即：

\[ u\approx678\,\mathrm{m/s},\quad \rho\approx1.23\,\mathrm{kg/m^3},\quad p\approx1.01\times10^5\,\mathrm{N/m^2},\quad T\approx286\,\mathrm{K},\quad M=2 \]

这里的推进变量不是时间，而是：

\[ \xi \]

因此，MacCormack 方法用于空间推进：

\[ \xi_i\rightarrow \xi_{i+1} \]

类比前面时间推进中的预测—校正格式：

预测步：沿 \(\eta\) 方向采用前向差分
校正步：沿 \(\eta\) 方向采用后向差分

当前先给出预测步。

由：

\[ \frac{\partial F_1}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F_1}{\partial\eta} + \frac{1}{h}\frac{\partial G_1}{\partial\eta} \right] \]

预测步采用前向差分：

\[ \left( \frac{\partial F_1}{\partial \xi} \right)_{i,j}= \left(\frac{\partial\eta}{\partial x}\right)_{i,j} \frac{(F_1)_{i,j}-(F_1)_{i,j+1}}{\Delta\eta} + \frac{1}{h} \frac{(G_1)_{i,j}-(G_1)_{i,j+1}}{\Delta\eta} \]

由：

\[ \frac{\partial F_2}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F_2}{\partial\eta} + \frac{1}{h}\frac{\partial G_2}{\partial\eta} \right] \]

预测步为：

\[ \left( \frac{\partial F_2}{\partial \xi} \right)_{i,j}= \left(\frac{\partial\eta}{\partial x}\right)_{i,j} \frac{(F_2)_{i,j}-(F_2)_{i,j+1}}{\Delta\eta} + \frac{1}{h} \frac{(G_2)_{i,j}-(G_2)_{i,j+1}}{\Delta\eta} \]

由：

\[ \frac{\partial F_3}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F_3}{\partial\eta} + \frac{1}{h}\frac{\partial G_3}{\partial\eta} \right] \]

预测步为：

\[ \left( \frac{\partial F_3}{\partial \xi} \right)_{i,j}= \left(\frac{\partial\eta}{\partial x}\right)_{i,j} \frac{(F_3)_{i,j}-(F_3)_{i,j+1}}{\Delta\eta} + \frac{1}{h} \frac{(G_3)_{i,j}-(G_3)_{i,j+1}}{\Delta\eta} \]

由：

\[ \frac{\partial F_4}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F_4}{\partial\eta} + \frac{1}{h}\frac{\partial G_4}{\partial\eta} \right] \]

预测步为：

\[ \left( \frac{\partial F_4}{\partial \xi} \right)_{i,j}= \left(\frac{\partial\eta}{\partial x}\right)_{i,j} \frac{(F_4)_{i,j}-(F_4)_{i,j+1}}{\Delta\eta} + \frac{1}{h} \frac{(G_4)_{i,j}-(G_4)_{i,j+1}}{\Delta\eta} \]

目前已经完成：

❶ 物理问题：Prandtl–Meyer 膨胀波
❷ 解析解：\(f_2=f_1+\theta\)
❸ 控制方程：二维定常 Euler 强守恒形式
❹ 推进变量：\(F_1,F_2,F_3,F_4\)
❺ 坐标变换：\((x,y)\rightarrow(\xi,\eta)\)
❻ 数值方法：MacCormack 空间推进预测步

后续需要继续处理：

• 校正步
• 边界条件
• 由 \(F\) 恢复物理变量
• 数值结果与解析解对比

由预测步得到的 \(\partial F_k/\partial \xi\)，可将通量变量从站位 \(i\) 推进到预测站位 \(i+1\)：

\[ (\overline{F}_1)_{i+1,j} =(F_1)_{i,j} +\left( \frac{\partial F_1}{\partial \xi} \right)_{i,j}\Delta\xi \]

\[ (\overline{F}_2)_{i+1,j} =(F_2)_{i,j} +\left( \frac{\partial F_2}{\partial \xi} \right)_{i,j}\Delta\xi \]

\[ (\overline{F}_3)_{i+1,j} =(F_3)_{i,j} +\left( \frac{\partial F_3}{\partial \xi} \right)_{i,j}\Delta\xi \]

\[ (\overline{F}_4)_{i+1,j} =(F_4)_{i,j} +\left( \frac{\partial F_4}{\partial \xi} \right)_{i,j}\Delta\xi \]

这里的上划线表示预测值。
它们还不是最终的 \(i+1\) 站位结果，需要经过校正步修正。

进入校正步之前，需要先由预测通量变量恢复预测状态量。

预测密度由二次关系得到：

\[ (\overline{\rho})_{i+1,j}= \frac{-B+\sqrt{B^2-4AC}}{2A} \]

其中：

\[ A=\frac{ (\overline{F}_3)^2_{i+1,j} }{ 2(\overline{F}_1)^2_{i+1,j} } -(\overline{F}_4)_{i+1,j} \]

\[ B=\frac{\gamma}{\gamma-1} (\overline{F}_1)_{i+1,j} (\overline{F}_2)_{i+1,j} \]

\[ C= -\frac{\gamma+1}{2(\gamma-1)} (\overline{F}_1)^3_{i+1,j} \]

得到 \((\overline{\rho})_{i+1,j}\) 后，可构造预测的 \(G\) 通量。

\[ (\overline{G}_1)_{i+1,j}= (\overline{\rho})_{i+1,j} \left( \frac{\overline{F}_3}{\overline{F}_1} \right)_{i+1,j} \]

\[ (\overline{G}_2)_{i+1,j}= (\overline{F}_3)_{i+1,j} \]

\[ (\overline{G}_3)_{i+1,j}= (\overline{\rho})_{i+1,j} \left( \frac{\overline{F}_3}{\overline{F}_1} \right)^2_{i+1,j}+ (\overline{F}_2)_{i+1,j}- \frac{ (\overline{F}_1)^2_{i+1,j} }{ (\overline{\rho})_{i+1,j} } \]

\[ (\overline{G}_4)_{i+1,j}= \frac{\gamma}{\gamma-1} \left[ (\overline{F}_2)_{i+1,j}- \frac{ (\overline{F}_1)^2_{i+1,j} }{ (\overline{\rho})_{i+1,j} } \right] \left( \frac{\overline{F}_3}{\overline{F}_1} \right)_{i+1,j}+ \frac{(\overline{\rho})_{i+1,j}}{2} \left( \frac{\overline{F}_3}{\overline{F}_1} \right)_{i+1,j} \left[ \left( \frac{\overline{F}_1}{\overline{\rho}} \right)^2_{i+1,j}+ \left( \frac{\overline{F}_3}{\overline{F}_1} \right)^2_{i+1,j} \right] \]

到此为止，校正步所需的预测通量 \(\overline{F}\) 和 \(\overline{G}\) 已经准备好。

校正步仍从变换后的控制方程出发：

\[ \frac{\partial F_1}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F_1}{\partial\eta}+ \frac{1}{h} \frac{\partial G_1}{\partial\eta} \right] \]

此时采用预测值，并在 \(\eta\) 方向使用后向差分：

\[ \left( \frac{\overline{\partial F_1}}{\partial \xi} \right)_{i+1,j}= \left(\frac{\partial\eta}{\partial x}\right) \frac{ (\overline{F}_1)_{i+1,j-1}- (\overline{F}_1)_{i+1,j} }{\Delta\eta}+ \frac{1}{h} \frac{ (\overline{G}_1)_{i+1,j-1}- (\overline{G}_1)_{i+1,j} }{\Delta\eta} \]

类似地，有：

\[ \left( \frac{\overline{\partial F_2}}{\partial \xi} \right)_{i+1,j}= \left(\frac{\partial\eta}{\partial x}\right) \frac{ (\overline{F}_2)_{i+1,j-1}- (\overline{F}_2)_{i+1,j} }{\Delta\eta}+ \frac{1}{h} \frac{ (\overline{G}_2)_{i+1,j-1}- (\overline{G}_2)_{i+1,j} }{\Delta\eta} \]

\[ \left( \frac{\overline{\partial F_3}}{\partial \xi} \right)_{i+1,j}= \left(\frac{\partial\eta}{\partial x}\right) \frac{ (\overline{F}_3)_{i+1,j-1}- (\overline{F}_3)_{i+1,j} }{\Delta\eta}+ \frac{1}{h} \frac{ (\overline{G}_3)_{i+1,j-1}- (\overline{G}_3)_{i+1,j} }{\Delta\eta} \]

\[ \left( \frac{\overline{\partial F_4}}{\partial \xi} \right)_{i+1,j}= \left(\frac{\partial\eta}{\partial x}\right) \frac{ (\overline{F}_4)_{i+1,j-1}- (\overline{F}_4)_{i+1,j} }{\Delta\eta}+ \frac{1}{h} \frac{ (\overline{G}_4)_{i+1,j-1}- (\overline{G}_4)_{i+1,j} }{\Delta\eta} \]

预测步和校正步得到两个空间导数。

对它们取平均：

\[ \left( \frac{\partial F_1}{\partial \xi} \right)_{\mathrm{av}}= \frac{1}{2} \left[ \left( \frac{\partial F_1}{\partial \xi} \right)_{i,j}+ \left( \frac{\overline{\partial F_1}}{\partial \xi} \right)_{i+1,j} \right] \]

\[ \left( \frac{\partial F_2}{\partial \xi} \right)_{\mathrm{av}}= \frac{1}{2} \left[ \left( \frac{\partial F_2}{\partial \xi} \right)_{i,j}+ \left( \frac{\overline{\partial F_2}}{\partial \xi} \right)_{i+1,j} \right] \]

\[ \left( \frac{\partial F_3}{\partial \xi} \right)_{\mathrm{av}}= \frac{1}{2} \left[ \left( \frac{\partial F_3}{\partial \xi} \right)_{i,j}+ \left( \frac{\overline{\partial F_3}}{\partial \xi} \right)_{i+1,j} \right] \]

\[ \left( \frac{\partial F_4}{\partial \xi} \right)_{\mathrm{av}}= \frac{1}{2} \left[ \left( \frac{\partial F_4}{\partial \xi} \right)_{i,j}+ \left( \frac{\overline{\partial F_4}}{\partial \xi} \right)_{i+1,j} \right] \]

得到平均空间导数后，更新到下一个站位：

\[ (F_1)_{i+1,j}= (F_1)_{i,j}+ \left( \frac{\partial F_1}{\partial \xi} \right)_{\mathrm{av}} \Delta\xi \]

\[ (F_2)_{i+1,j}= (F_2)_{i,j}+ \left( \frac{\partial F_2}{\partial \xi} \right)_{\mathrm{av}} \Delta\xi \]

\[ (F_3)_{i+1,j}= (F_3)_{i,j}+ \left( \frac{\partial F_3}{\partial \xi} \right)_{\mathrm{av}} \Delta\xi \]

\[ (F_4)_{i+1,j}= (F_4)_{i,j}+ \left( \frac{\partial F_4}{\partial \xi} \right)_{\mathrm{av}} \Delta\xi \]

本问题中，膨胀角位于：

\[ x=10\,\mathrm{m} \]

该点是一个几何奇异点，壁面方向突然改变。

在控制方程中，这种突变体现为度量项：

\[ \frac{\partial\eta}{\partial x} \]

在膨胀角前后不连续。

这种不连续会在数值解中诱发振荡。
因此需要加入人工粘性来平滑高频误差。

预测步中，对 \(F_1\) 和 \(F_2\) 的人工粘性项可写为：

\[ (SF_1)_{i,j}= \frac{ C_y|p_{i,j+1}-2p_{i,j}+p_{i,j-1}| }{ p_{i,j+1}+2p_{i,j}+p_{i,j-1} } \left[ (F_1)_{i,j+1} -2(F_1)_{i,j} +(F_1)_{i,j-1} \right] \]

\[ (SF_2)_{i,j}= \frac{ C_y|p_{i,j+1}-2p_{i,j}+p_{i,j-1}| }{ p_{i,j+1}+2p_{i,j}+p_{i,j-1} } \left[ (F_2)_{i,j+1} -2(F_2)_{i,j} +(F_2)_{i,j-1} \right] \]

\(SF_3\) 和 \(SF_4\) 具有相同形式。

加入人工粘性后，预测更新变为：

\[ (\overline{F}_1)_{i+1,j}= (F_1)_{i,j}+ \left( \frac{\partial F_1}{\partial \xi} \right)_{i,j} \Delta\xi+ (SF_1)_{i,j} \]

\[ (\overline{F}_2)_{i+1,j}= (F_2)_{i,j}+ \left( \frac{\partial F_2}{\partial \xi} \right)_{i,j} \Delta\xi+ (SF_2)_{i,j} \]

类似地：

\[ (\overline{F}_3)_{i+1,j}= (F_3)_{i,j}+ \left( \frac{\partial F_3}{\partial \xi} \right)_{i,j} \Delta\xi+ (SF_3)_{i,j} \]

\[ (\overline{F}_4)_{i+1,j}= (F_4)_{i,j}+ \left( \frac{\partial F_4}{\partial \xi} \right)_{i,j} \Delta\xi+ (SF_4)_{i,j} \]

校正步也要加入人工粘性。此时使用预测变量：

\[ (\overline{SF}_1)_{i+1,j}= \frac{ C_y|\overline{p}_{i+1,j+1}-2\overline{p}_{i+1,j}+\overline{p}_{i+1,j-1}| }{ \overline{p}_{i+1,j+1}+2\overline{p}_{i+1,j}+\overline{p}_{i+1,j-1} } \left[ (\overline{F}_1)_{i+1,j+1} -2(\overline{F}_1)_{i+1,j} +(\overline{F}_1)_{i+1,j-1} \right] \]

\[ (\overline{SF}_2)_{i+1,j}= \frac{ C_y|\overline{p}_{i+1,j+1}-2\overline{p}_{i+1,j}+\overline{p}_{i+1,j-1}| }{ \overline{p}_{i+1,j+1}+2\overline{p}_{i+1,j}+\overline{p}_{i+1,j-1} } \left[ (\overline{F}_2)_{i+1,j+1} -2(\overline{F}_2)_{i+1,j} +(\overline{F}_2)_{i+1,j-1} \right] \]

\(\overline{SF}_3\) 和 \(\overline{SF}_4\) 类似。

最终更新时加入校正步人工粘性：

\[ (F_1)_{i+1,j}= (F_1)_{i,j}+ \left( \frac{\partial F_1}{\partial \xi} \right)_{\mathrm{av}} \Delta\xi+ (\overline{SF}_1)_{i+1,j} \]

\[ (F_2)_{i+1,j}= (F_2)_{i,j}+ \left( \frac{\partial F_2}{\partial \xi} \right)_{\mathrm{av}} \Delta\xi + (\overline{SF}_2)_{i+1,j} \]

类似地可得到 \((F_3)_{i+1,j}\) 和 \((F_4)_{i+1,j}\)。

完成这一步后，内部网格点 \(j=2\) 到 \(j=40\) 的流场变量已经可由 \(F_1\) 到 \(F_4\) 解码得到。

对无粘流，壁面处正确的物理条件是：

\[ \mathbf{V}\cdot\mathbf{n}=0 \]

也即速度必须与壁面相切。

这是壁面处唯一必须施加的物理边界条件。
其他变量应由内部流场和相容关系决定。

但在数值计算中，壁面位于边界点，缺少壁面外侧网格点。
因此需要特殊处理。

壁面试算 Mach 数为：

\[ (M_1)_{\mathrm{cal}}= \frac{ \sqrt{(u_1)^2_{\mathrm{cal}}+(v_1)^2_{\mathrm{cal}}} }{ (a_1)_{\mathrm{cal}} } \]

由此可计算对应的 Prandtl–Meyer 函数：

\[ f_{\mathrm{cal}}=f\left((M_1)_{\mathrm{cal}}\right) \]

接着设想通过一个局部 Prandtl–Meyer 波，把试算速度旋转到壁面切向。

这一步不是说壁面真实存在一个膨胀波，
而是用一个局部等熵修正来消除数值计算产生的法向速度。

对于膨胀角上游壁面，修正后的 Prandtl–Meyer 函数为：

\[ f_{\mathrm{act}}=f_{\mathrm{cal}}+\phi_1 \]

由 \(f_{\mathrm{act}}\) 反求：

\[ (M_1)_{\mathrm{act}} \]

然后利用等熵关系修正壁面热力学量：

\[ p_{\mathrm{act}}=p_{\mathrm{cal}} \left[ \frac{ 1+\frac{\gamma-1}{2}M^2_{\mathrm{cal}} }{ 1+\frac{\gamma-1}{2}M^2_{\mathrm{act}} } \right]^{\gamma/(\gamma-1)} \]

\[ T_{\mathrm{act}}= T_{\mathrm{cal}} \frac{ 1+\frac{\gamma-1}{2}M^2_{\mathrm{cal}} }{ 1+\frac{\gamma-1}{2}M^2_{\mathrm{act}} } \]

\[ \rho_{\mathrm{act}}= \frac{p_{\mathrm{act}}}{RT_{\mathrm{act}}} \]

膨胀角下游，壁面本身有倾角 \(\theta\)。

设试算速度方向与水平线夹角为：

\[ \psi=\tan^{-1}\frac{|v_2|}{u_2} \]

为了使速度最终与壁面相切，需要旋转的角度为：

\[ \phi_2=\theta-\psi \]

相应地：

\[ f_{\mathrm{act}}=f_{\mathrm{cal}}+\phi_2 \]

若试算速度指向壁内，则 \(\phi_2\) 可能为负。
这相当于通过一个局部等熵压缩修正，使速度重新与壁面相切。

在时间推进中，CFL 条件限制的是：

\[ \Delta t \]

在下游推进中，限制的是：

\[ \Delta x \quad\text{或}\quad \Delta \xi \]

原因相同：数值推进不能越过物理信息传播的影响范围。

对超声定常流，信息沿 Mach 线传播。
因此下游推进步长必须与 Mach 线斜率相容。

为了保证所有网格点都满足稳定性要求，应取最严格的限制。

因此：

\[ \Delta x= \frac{\Delta y}{ |\tan(\theta\pm\mu)|_{\max} } \]

在计算平面中，由于：

\[ \xi=x \]

推进步长应满足：

\[ \Delta\xi\leq\Delta x \]

引入 Courant 数 \(C\)，得到：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \Delta\xi=C \frac{\Delta y}{ |\tan(\theta\pm\mu)|_{\max} } } $} \]

其中：

\[ C\leq1 \]

公式：

\[ \Delta\xi=C \frac{\Delta y}{ |\tan(\theta\pm\mu)|_{\max} } \]

说明：

• \(\Delta y\) 越小，允许的推进步长越小
• Mach 线越陡，允许的推进步长越小
• 局部流动转角和 Mach 角共同决定稳定性限制
• \(C\) 起到安全系数作用

它是时间推进 CFL 条件在定常超声空间推进中的对应形式。

到此，二维 Prandtl–Meyer 膨胀波的空间推进算法已经基本完整。

算法流程：

❶ 已知站位 \(i\) 的 \(F\)
❷ 构造 \(G\)
❸ 预测步得到 \(\overline{F}\)
❹ 解码 \(\overline{\rho}\) 并构造 \(\overline{G}\)
❺ 校正步得到平均导数
❻ 加入人工粘性
❼ 更新 \(F\) 到站位 \(i+1\)
❽ 解码得到 \(\rho,u,v,p,T\)
❾ 对壁面点施加无穿透边界条件
❿ 根据 CFL 条件确定下一步 \(\Delta\xi\)

请思考：

1. 为什么该问题可以沿 \(x\) 方向推进，而不是必须全场联立求解？
   
2. 为什么要用强守恒形式中的 \(F_1,F_2,F_3,F_4\) 作为推进变量？
   
3. 膨胀角处为什么容易产生数值振荡？
   
4. 壁面处为什么不是指定 \(u\) 和 \(v\)，而是只要求速度与壁面相切？
   
5. 空间推进中的 CFL 条件和时间推进中的 CFL 条件有什么共同点？

教材将 \(x\) 方向速度分量 \(u\) 画成不同流向站位上的竖向剖面。

典型站位包括：

\[ x=0,\quad 16.17,\quad 32.31,\quad 48.99,\quad 66.23\,\mathrm{m} \]

图中同时给出：

• 黑点：Prandtl–Meyer 解析解
• 实线：下游推进数值解
• 斜线区域：膨胀波前缘与后缘
• 下边界：膨胀角后的壁面形状

该图的目的，是直接比较数值方法能否正确捕捉膨胀波内部和波后的速度分布。

膨胀角附近存在两个困难：

❶ 几何奇异性

\[ \frac{\partial\eta}{\partial x} \]

在膨胀角前后不连续。

❷ 网格分辨率不足

刚过膨胀角时，膨胀波很窄，壁面与膨胀波后缘之间只有少数几个网格点。

结果是：数值格式在最困难的位置“正面撞上”奇异点，
而且可用于吸收误差的网格点又很少。

整体网格点数大约增加为原来的 4 倍。

如果加密后的结果与原计算没有显著差别，说明前面的结果基本达到了网格无关性。

本算例中，计算区域尺寸约为：

\[ 65\,\mathrm{m}\times40\,\mathrm{m} \]

这个长度尺度本身并不是问题的核心。

如果采用毫米作为单位，也可以把计算区域写成：

\[ 65\,\mathrm{mm}\times40\,\mathrm{mm} \]

只要所有量使用一致单位，
数值方法本身并不依赖某一个特定长度尺度。

本章的核心是：

\[ \colorbox{yellow}{$\displaystyle \color{black}{ \text{steady supersonic flow} \Rightarrow \text{space marching} } $} \]

对于定常超声流，信息沿下游传播，因此可以从初始数据线出发，沿 \(x\) 或 \(\xi\) 方向逐步推进。

这与前一章的时间推进不同。
这里的“推进变量”不是时间，而是空间坐标。

控制方程采用强守恒形式：

\[ \frac{\partial F}{\partial x}= -\frac{\partial G}{\partial y} \]

在贴体坐标中变为：

\[ \frac{\partial F}{\partial \xi}= -\left[ \left(\frac{\partial\eta}{\partial x}\right) \frac{\partial F}{\partial\eta} + \frac{1}{h} \frac{\partial G}{\partial\eta} \right] \]

强守恒形式使得数值方法直接推进通量变量 \(F_1,F_2,F_3,F_4\)，
这对于波捕捉问题尤其重要。

由于膨胀角后壁面发生倾斜，物理区域不是规则矩形。

因此引入变换：

\[ \xi=x,\qquad \eta=\frac{y-y_s(x)}{h(x)} \]

将物理区域映射到规则计算区域。

贴体坐标的作用是：
让复杂边界变成简单的坐标线。

本章捕捉的是膨胀波，而不是激波。

但仍然可以看到波捕捉方法的基本特点：

• 使用守恒形式
• 允许波系在网格中自然展开
• 在奇异点附近加入人工粘性
• 通过解析解检查结果

人工粘性主要用于膨胀角附近，
远离奇异点后其作用可以很小。

无粘壁面条件为：

\[ \mathbf{V}\cdot\mathbf{n}=0 \]

即速度与壁面相切。

教材采用 Abbott 的壁面处理方法，通过局部 Prandtl–Meyer 修正，使壁面速度满足切向条件。

边界条件看似只是局部处理，
但在下游推进问题中，边界误差会沿流向传播和积累。