平板上超声速流动求解示例程序实验报告

本文分析零攻角平板上二维超声速黏性层流的数值求解。讨论对象是本仓库的 C++ 串行程序，算例参数取自 Anderson《计算流体力学入门》中文版第 9 章的等温壁示例。项目仓库：hYdr0xY-Zeng/CFD。

求解的问题

自由来流以马赫数 $M_\infty=4$ 平行于平板进入计算域。平板从 $x=0$ 开始，攻角为零，壁面无滑移且保持来流温度 $T_w=T_\infty$ 。求解目标是得到稳定后的密度、速度、压力和温度场，观察前缘附近的压缩、黏性边界层及其相互作用。

程序求解二维、层流、可压缩、量热完全气体的完整 Navier-Stokes 方程；忽略体积力和体积热源，采用 Stokes 黏性假设与 Fourier 导热定律。守恒形式为

\frac{\partial\mathbf U}{\partial t} +\frac{\partial\mathbf E}{\partial x} +\frac{\partial\mathbf F}{\partial y}=0, \qquad \mathbf U=\begin{bmatrix}\rho\\\rho u\\\rho v\\E_t\end{bmatrix}, \quad E_t=\rho\left(c_vT+\frac{u^2+v^2}{2}\right).

其中 $\mathbf E$ 、 $\mathbf F$ 分别是 $x$ 、 $y$ 方向通量。以 $x$ 方向为例，程序的 fluxes 函数计算

\mathbf E=\begin{bmatrix} \rho u\\ \rho u^2+p-\tau_{xx}\\ \rho uv-\tau_{xy}\\ (E_t+p)u-u\tau_{xx}-v\tau_{xy}+q_x \end{bmatrix}, \qquad \mathbf F=\begin{bmatrix} \rho v\\ \rho uv-\tau_{xy}\\ \rho v^2+p-\tau_{yy}\\ (E_t+p)v-u\tau_{xy}-v\tau_{yy}+q_y \end{bmatrix}.

连续方程中的黏性应力及热通量为

\begin{aligned} \tau_{xx}&=2\mu u_x-\frac23\mu(u_x+v_y),\\ \tau_{yy}&=2\mu v_y-\frac23\mu(u_x+v_y),\\ \tau_{xy}&=\mu(u_y+v_x),\\ q_x&=-kT_x,\qquad q_y=-kT_y. \end{aligned}

这里 $u_x$ 表示 $\partial u/\partial x$ 。状态方程为 $p=\rho RT$ ，并取 $c_v=R/(\gamma-1)$ 、 $c_p=\gamma c_v$ 。时间推进的作用是使初始场逐步趋近稳态，并非研究该平板流动的真实瞬态过程。

前处理

参数与计算域

参数	C++ 默认值	含义
`mach_inf`	4	自由来流马赫数
`sound_inf`	$340.28\;\mathrm{m/s}$	输入的自由来流声速
`pressure_inf`	$101325\;\mathrm{Pa}$	自由来流静压
`temperature_inf`	$288.16\;\mathrm K$	自由来流静温
`gamma`, `gas_constant`, `prandtl`	$1.4$ ， $287\;\mathrm{J/(kg\,K)}$ ， $0.71$	空气物性常数
`viscosity_ref`, `sutherland`	$1.7894\times10^{-5}\;\mathrm{kg/(m\,s)}$ ， $110\;\mathrm K$	Sutherland 公式参数
`plate_length`, `wall_temperature_ratio`	$10^{-5}\;\mathrm m$ ， $1$	板长和壁温比
`nx`, `ny`, `cfl`	$70$ ， $70$ ， $0.6$	节点数和时间步系数

Physics 结构体由这些常数计算 $u_\infty=M_\infty a_\infty=1361.12\;\mathrm{m/s}$ 、 $\rho_\infty=p_\infty/(RT_\infty)=1.22518\;\mathrm{kg/m^3}$ 以及基于板长的雷诺数

Re_L=\frac{\rho_\infty u_\infty L}{\mu_\infty}=931.944.

书中将雷诺数概述为“约 1000”；程序用上述物性计算值确定网格高度，没有直接把 1000 当作精确输入。依照书中的边界层厚度估计， $\delta=5L/\sqrt{Re_L}=1.638\times10^{-6}\;\mathrm m$ ，计算域高度取 $H=5\delta=8.189\times10^{-6}\;\mathrm m$ 。Grid 使用均匀的 $70\times70$ 节点网格，因此 $\Delta x=L/69=1.4493\times10^{-7}\;\mathrm m$ 、 $\Delta y=H/69=1.1869\times10^{-7}\;\mathrm m$ 。

网格节点用 $(i,j)$ 表示，其中 $i$ 沿平板， $j$ 垂直平板。Grid::id(j,i)=j\,n_x+i 将二维节点映射为一维下标，Grid::sid(j,i,c)=4\,\mathrm{id}(j,i)+c 再定位四个守恒分量。因而 State 是按“节点、分量”排列的 std::vector<double>，保证了访问时的缓存友好，Primitive 则分别保存 $\rho,u,v,T,p,\mu,k$ 等标量场。这种布局使 advance 中每个节点的四个守恒量能够依次更新。

黏度由 Physics::viscosity 按 Sutherland 公式计算，导热系数由固定 Prandtl 数得到：

\mu(T)=\mu_\infty\left(\frac{T}{T_\infty}\right)^{3/2} \frac{T_\infty+S}{T+S}, \qquad k(T)=\frac{\mu(T)c_p}{Pr}.

初始条件与边界条件

initial_state 先把各节点设为自由来流，再调用 apply_boundaries。边界值用 put 转换回守恒变量，以保证边界上的密度、动量和总能量一致。

边界位置	`apply_boundaries` 中的处理
左侧入口，除前缘	$u=u_\infty$ 、 $v=0$ 、 $p=p_\infty$ 、 $T=T_\infty$
上边界	固定为自由来流状态
下方平板，除前缘	$u=v=0$ 、 $T=T_w$ ；壁压由 $p_{i,0}=2p_{i,1}-p_{i,2}$ 外插
前缘 $(0,0)$	$u=v=0$ 、 $p=p_\infty$ 、 $T=T_\infty$
右侧出口，除上下角点	对 $u,v,p,T$ 沿 $x$ 线性外插，即 $\phi_{n_x-1,j}=2\phi_{n_x-2,j}-\phi_{n_x-3,j}$

入口位于左侧、出口位于右侧，而前缘单独处理；出口上角点服从上边界，下角点服从壁面条件。所有给定 $p,T$ 的边界均通过 $\rho=p/(RT)$ 得到密度。

计算求解

原始变量、通量与差分

每个时间步先由 decode 从 $\mathbf U$ 恢复速度、内能、温度与压力，并检查密度、温度是否为正且数值有限。fluxes 调用 difference 计算速度和温度梯度，再构造包含对流、压力、黏性应力与热通量的 $\mathbf E,\mathbf F$ 。

为了配合 MacCormack 格式，预测步在通量外层使用向前差分，而黏性通量中同方向的一阶导数使用向后差分；校正步把这两组方向对换。交叉方向导数在内点使用中心差分。difference 对缺少外侧节点的边界采用单侧差分。连续方程中的 $\tau_{xy}$ 是同一个物理量，代码中 tau_xy_e 与 tau_xy_f 的写法不同，是因为两个通量方向选用的离散导数不同。

利用相反方向的差分抵消一阶截断误差项, 提高到二阶精度

时间步与 MacCormack 预测校正

time_step 在所有内点计算局部稳定性限制(CFL条件)，再取最小允许时间步。其实现可写为

\Delta t=\frac{K}{\displaystyle\max_{i,j}\left[ \frac{|u|}{\Delta x}+\frac{|v|}{\Delta y} +a\sqrt{\frac1{\Delta x^2}+\frac1{\Delta y^2}} +2\nu'\left(\frac1{\Delta x^2}+\frac1{\Delta y^2}\right) \right]_{i,j}}, \qquad \nu'=\frac{\mu}{\rho}\max\left(\frac43,\frac\gamma{Pr}\right), \quad K=0.6.

式中 $a=\sqrt{\gamma RT}$ 是当地声速， $\nu'$ 是用于估算黏性及热扩散限制的有效运动黏度。

advance 只对内点更新守恒量。预测步使用当前场的向前通量差分：

\mathbf U^*_{i,j}=\mathbf U^n_{i,j} -\frac{\Delta t}{\Delta x}(\mathbf E^n_{i+1,j}-\mathbf E^n_{i,j}) -\frac{\Delta t}{\Delta y}(\mathbf F^n_{i,j+1}-\mathbf F^n_{i,j}).

程序随即对预测场施加边界条件，重新计算预测场通量。校正步使用向后通量差分并与预测值取平均：

\mathbf U^{n+1}_{i,j}=\frac12\left[ \mathbf U^n_{i,j}+\mathbf U^*_{i,j} -\frac{\Delta t}{\Delta x}(\mathbf E^*_{i,j}-\mathbf E^*_{i-1,j}) -\frac{\Delta t}{\Delta y}(\mathbf F^*_{i,j}-\mathbf F^*_{i,j-1}) \right].

校正后再次执行 apply_boundaries。main 在每一步比较**全网格（含边界）**新旧密度，满足 $\max|\rho^{n+1}-\rho^n|<10^{-8}\;\mathrm{kg/m^3}$ 时停止(对应程序中的变量drho)。若使用默认关闭的 --safe-mode，一次推进出现非物理状态后会把该步 $\Delta t$ 减半重试，最多 8 次；本报告的默认结果没有触发重试。

编译与运行

在 WSL 中从算例目录运行 C++ 版本：

cd /home/hydr0xy/CFD/cases/anderson-ch9-flat-plate
mkdir -p build
g++ -std=c++17 -O2 -Wall -Wextra -pedantic cpp/main.cpp -o build/flat_plate
./build/flat_plate --nx 70 --ny 70 --cfl 0.6 --max-steps 10000 --out results/cpp

write_outputs 生成 results/cpp/field.csv 和 results/cpp/history.csv。前者按 $y$ 行、 $x$ 列列出坐标及 $\rho,u,v,p,T,M$ ；后者记录步数、时间步、最大密度变化、进出口质量流率、相对差与重试次数。

后处理

守恒与收敛检查

以下数值直接读取 C++ 版输出的 field.csv 与 history.csv。mass_flow 对入口和出口的 $\rho u$ 沿 $y$ 使用梯形积分，按单位展向深度计算质量流率，并用 $|\dot m_{out}-\dot m_{in}|/|\dot m_{in}|$ 表示相对差。由于初始条件的设定,上下边界无质量流入/出，因此计算程序中mass_flow函数计算的是一列网格的流入/出。

网格	收敛步数	最后一步最大密度变化（kg/m³）	入口 / 出口质量流率（kg/(m·s)）	相对差
70 × 70	3745	$9.9655\times10^{-9}$	0.0135576 / 0.0136522	0.6973%
40 × 50	1999	$9.9956\times10^{-9}$	0.0135172 / 0.0136443	0.9398%

两个网格均达到设定的逐步密度变化阈值，且质量流率相对差小于 1%。这两项检查说明当前离散解满足本算例选定的停止和整体守恒指标；它们本身不能证明局部激波位置或所有细节已经网格收敛。

流场曲线及物理解释

下列图片是仓库已有的自绘结果图，均显示等温壁算例；本文对图中量级的核算采用 C++ 输出。仓库绘图脚本也支持把 results/cpp/field.csv 作为输入重绘这些图。横向剖面中的无量纲高度取尾缘 $x=L$ 处的 $\bar y=(y/L)\sqrt{Re_L}$ 。

图 1. 壁面与尾缘压力。前缘后方出现明显压力峰，C++ 场的最大壁压约为 $2.922p_\infty$ ；沿板向下游降低，尾缘壁压约为 $1.365p_\infty$ 。尾缘剖面从壁面向外先升后降，最后趋向自由来流压力。这些特征反映前缘压缩与黏性层附近的相互作用；下游壁压曲线与书中照片并未逐点吻合。

图 2. 尾缘温度。壁面强制保持 $T_w/T_\infty=1$ ，壁面上方温度先升高，再向远场回落；C++ 场的尾缘最高温度约为 $1.605T_\infty$ 。温度峰与高速黏性流中的压缩和耗散有关，但单凭这条剖面不能分离两者各自的贡献。

图 3. 尾缘 $x$ 方向速度。 $u/u_\infty$ 从无滑移壁面的 0 向外增长，逐渐接近 1，体现边界层的速度梯度。

图 4. 尾缘当地马赫数。壁面速度为零，所以当地马赫数为零；远场趋近设定的马赫数 4。马赫数同时依赖速度与当地温度，不能仅由图 3 的速度曲线按常数比例换算。

40 × 50 网格的前缘附近仍有压力峰，尾缘温度、速度和马赫数的主要趋势也保留。但最大壁压约为 $2.554p_\infty$ ，明显低于 70 × 70 网格的 $2.922p_\infty$ ；相比之下，尾缘壁压分别约为 $1.370p_\infty$ 与 $1.365p_\infty$ 。前缘峰值对网格较敏感，两个网格不足以给出严格的网格收敛阶或外推误差。

总结

C++ 程序把守恒量存储、物性计算、黏性通量、MacCormack 预测校正和边界更新串成了可追踪的求解流程。默认 70 × 70 等温壁算例达到 $10^{-8}\;\mathrm{kg/m^3}$ 的全网格密度变化阈值，进出口质量流率相对差为 0.6973%；结果图展示了前缘压力峰、壁面温度约束、速度边界层和当地马赫数的变化。

书中称等温壁算例约在 4339 步收敛，本程序为 3745 步；壁面下游压力也存在可见差异。边界离散、黏性导数处理和停止条件细节都可能影响这些差异，目前没有证据把差异归因于单一原因。该结果适合作为完整 Navier-Stokes 教学求解例和进一步验证的起点，不应表述为对书中曲线的严格逐点复现。

资料与代码： Anderson《计算流体力学入门》中文版第 9 章（用户提供的书页照片）；算例说明与参考资料；本仓库 C++ 源码；GitHub 仓库。