简介：拉瓦尔喷管准一维流动数值模拟程序是一份面向计算流体力学学习者的经典案例资源，对应《计算流体力学基础及应用》中的典型问题。程序以准一维假设为基础，用有限差分法求解收缩—扩张喷管内的绝热流动方程，完整覆盖控制方程离散、时间推进、边界条件处理和人工粘性设置等关键环节，可帮助读者掌握可压缩流动数值模拟的基本流程。压缩包共56个文件，大小仅594KB，其中以20个Fortran源文件为核心，配以14个模块文件、8个目标文件及工程配置文件，并附有压强分布、Excel数据等结果，可直接编译运行并对照验证。目前已有2392人学习，程序内部采用二维数组存储流场并兼顾守恒特性，引入人工粘性以改善激波等间断解的稳定性。借助完整源码与配套结果，读者可跟踪每一时间步的流场变化，观察从亚声速经喉部至超声速的加速过程，还可自行调整边界条件与人工粘性参数，深入理解有限差分法在可压缩流动中的实际应用。
1. 拉瓦尔喷管准一维数值模拟到底在算什么
做超声速喷管、火箭发动机喷管或者风洞试验段之前，很少有人上来就开三维 CFD。行业里最常用的第一步，是先拿准一维流动数值模拟程序把面积分布、喉部塞塞、激波位置这些核心问题算清楚。拉瓦尔喷管准一维流动数值模拟程序，就是把真实喷管近似成一条沿轴向变化的截面积管道，只保留轴向流动参数，用一维欧拉方程加面积变化源项去求解流场。这个程序能回答三个工程问题：给定面积比和总压总温，喷管出口马赫数是多少；背压变化时激波在扩张段哪个位置；喉部能不能达到声速、流量是否塞塞。
这类程序写起来难度不大，但要写稳、写对、能复现实验数据，藏着一堆细节。读者里如果是有可压缩流课程背景的学生，或者刚接手喷管性能评估的工程师，这篇文章会把理论、离散格式、边界条件和排错经验一次讲透，让你照着写就能跑出设计工况，还能拿激波关系式反过来验证程序结果。
2. 拉瓦尔喷管准一维流动的数学模型与守恒形式
2.1 守恒型欧拉方程为什么比非守恒型更可靠
拉瓦尔喷管内的流动可以看作无粘、绝热、无外加功的可压缩管流，控制方程取一维欧拉方程。常见的非守恒形式会写速度、压力、密度的偏微分方程组，形式直观，但带激波时激波位置和强度容易算偏，而且激波处非守恒形式会产生伪源项，造成质量、动量、能量不守恒。做数值模拟时，我一般直接写守恒形式：
∂Q/∂t + ∂F/∂x = S
其中守恒变量矢量 Q、通量矢量 F 和面积变化源项 S 分别是：
Q = [ρA, ρuA, ρEA] F = [ρuA, (ρu²+p)A, ρuA(E + p/ρ)] S = [0, p dA/dx, 0]
这里的 A = A(x) 是喷管截面积，ρ、u、p、E 分别是密度、速度、静压和单位质量总内能。对于理想气体，E = p/((γ-1)ρ) + u²/2，γ 取 1.4。如果写成无量纲形式，还可以把 A 从方程里归一化掉，这节最后会给具体公式。
为什么非要保留 dA/dx 这个源项？因为准一维流动的本质是面积变化驱动速度变化。收缩段面积减小，亚声速气流加速；喉部达到声速后，扩张段面积增大，气流继续加速到超声速。面积源项处理不当时，数值解会在喉部附近出现虚假的密度振荡，这是初写者最容易翻车的地方。
2.2 无量纲化、面积分布与初始场构造
计算前先把方程无量纲化，能避免量级差太大导致浮点精度问题。参考量可以取入口滞止状态：参考压力 p₀、参考密度 ρ₀、参考声速 a₀ = √(γRT₀)，时间参考量取 L/a₀。无量纲化后控制方程变成：
∂Q/∂t + ∂F/∂x = S Q = [ρ, ρu, ρE] F = [ρu, ρu²+p, ρu(E+p/ρ)] S = [-ρu d(lnA)/dx, -ρu² d(lnA)/dx, -(ρE+p)u d(lnA)/dx]
这样面积项变成了对数面积导数，实际编程时先算 d(lnA)/dx，既方便网格变换，也避免了 A 趋近于零时的除零风险。
面积分布函数是拉瓦尔喷管几何的核心输入。工程上常见做法是用多项式拟合几何数据，我习惯用三段式面积规律：收缩段从入口面积线性过渡到喉部，扩张段按双曲正切或幂函数过渡到出口面积。为了验证程序，可以用标准喷管面积比公式反推马赫数，比如给定喉部面积 A* = 1，入口面积 A_in = 2，出口面积 A_ex = 2.5，用等熵面积比关系：
A/A* = (1/M)[(2/(γ+1))(1 + (γ-1)M²/2)]^((γ+1)/(2(γ-1)))
这个公式可以用来生成初始场，也可以用数值结果反算马赫数来校验程序。
初始场不能随便填常值。直接给入口马赫数 0.3、全场密度压力均匀，时间推进几步就会在面积突变处产生强压缩波，程序直接发散。正确的做法是用等熵关系构造一个无激波的初始猜解：先假设入口为亚声速，从面积比公式解出每个截面的马赫数分布，再用等熵关系算出压力和密度。等熵关系用这两条公式：
T/T₀ = 1/(1 + (γ-1)M²/2) p/p₀ = (T/T₀)^(γ/(γ-1)) ρ/ρ₀ = (T/T₀)^(1/(γ-1))
这样初始场接近真实解，时间推进收敛快，也方便观察激波形成过程。
3. 数值离散格式：MacCormack 方法的应用与参数设定
3.1 MacCormack 格式：预测步与校正步的配合逻辑
准一维欧拉方程自带激波和大梯度区域，空间离散格式选型很重要。我常用的作为教学和验证首选的是 MacCormack 两步显式格式，它本质是 Lax-Wendroff 思想的差分实现，时间精度二阶、空间精度二阶，在激波不太强时能够稳定捕捉。对于初学喷管模拟的人，这个格式代码短、逻辑清楚、调参直观；如果需要更低的数值耗散或者更强的激波捕获能力，后续可以切换到 Roe 或 Godunov 型格式，后者在第 6 章里讨论。
MacCormack 格式的半步推进是这样写的：
预测步用前向差分：
Q_pred[i] = Q[i] - (dt/dx)(F[i+1] - F[i]) + dt*S[i]
校正步用后向差分：
Q_corr[i] = 0.5*(Q[i] + Q_pred[i] - (dt/dx)(F_pred[i] - F_pred[i-1]) + dt*S_pred[i])
最终 Q^{n+1} = Q_corr。这组式子有两点关键：第一，F_pred 必须用预测步得到的守恒变量重新计算，不能用上一时刻的通量；第二，源项 S_pred 同理要用预测后的变量算，否则校正步的时间精度掉回一阶。很多初写版本只在预测步更新源项，校正步沿用旧源项，最后收敛结果虽然差不多，但激波位置会偏移半个网格。
3.2 网格、时间步长与 CFL 条件的具体选择
网格划分直接影响激波分辨率和收敛速度。我一般先在喉部附近加密，因为喉部马赫数要穿过 1，梯度变化最剧烈。均匀网格要加密到 500～1000 个点才能看清喉部压力分布；如果做工程快速评估，300 个点配合喉部局部加密也够用。加密方式可以取双曲正切变换，也可以简单地令网格点在喉部附近按余弦分布加密，计算量差别很小，关键是别把喉部放在两个网格点中间，否则马赫数 = 1 的位置会来回跳动。
时间步长用 CFL 条件约束。对绝热无粘管流，最大特征波速是 |u| + a，所以时间步长取：
dt = CFL * dx / max(|u| + a)
显式格式 CFL 上限理论值是 1，但 MacCormack 格式有中间步稳定性限制，实际中 CFL 取 0.3～0.6 比较稳妥。注意这里的 a 是当地声速，不是入口参考声速，必须在每个网格点扫描一遍再取全局最大值。总计算时间取多少？通常取无量纲时间 t = 5～10，即整个喷管流场被扫过 5～10 个来回，边界信息传到全场并衰减完。用 500 个网格、CFL=0.5，这个计算量在 Python 里也就几秒钟。
CFL 取值有一个经验法则：激波位置捕捉时，CFL 太小会让激波位置缓慢爬行，太大又会让激波处产生过冲。初算用 0.5，若看到压力剖面在激波前后出现大幅度振荡，先降到 0.3 试跑，再检查是不是面积分布不光滑；若振荡还在，问题大概率出在边界条件而不是时间步长。
4. 程序骨架与设计工况算例
4.1 程序结构：面积分布、初始场推进和主循环
整个模拟程序拆成三个核心函数：面积分布函数、等熵初场函数和一个 MacCormack 时间推进函数。下面给出一套可以直接跑通的 Python 骨架，使用无量纲变量。代码里我刻意保留了核心注释，便于对照公式理解。
import numpy as np from math import sqrt, log GAMMA = 1.4 def area_distribution(x, L=1.0, area_in=2.0, area_throat=1.0, area_exit=2.5): """ 拉瓦尔喷管面积分布，三段式： 0-0.3L 收缩段，0.3L-0.5L 喉部过渡，0.5L-1.0L 扩张段 """ A = np.empty_like(x) for i, xi in enumerate(x): if xi < 0.3*L: A[i] = area_throat + (area_in - area_throat) * (1 - xi/(0.3*L))**2 elif xi < 0.5*L: A[i] = area_throat else: xi_norm = (xi - 0.5*L) / (0.5*L) A[i] = area_throat + (area_exit - area_throat) * xi_norm**2 return A def calc_dlnA_dx(A, x): """对数面积的导数，用于源项""" return np.gradient(np.log(A), x) def isentropic_initial_state(A, M_in=0.3, p0=1.0, rho0=1.0): """ 等熵关系构造初场 """ # 先通过面积比公式解每个截面的马赫数 M = np.empty_like(A) M[0] = M_in for i in range(1, len(A)): ratio = A[i] / A[0] # 用牛顿迭代解面积比公式 M_guess = M[i-1] for _ in range(30): f = (1/M_guess) * ((2/(GAMMA+1)) * (1 + (GAMMA-1)*M_guess**2/2))**((GAMMA+1)/(2*(GAMMA-1))) - ratio df = ... # 数值差分求导 M_guess -= f / df M[i] = M_guess T_ratio = 1/(1 + (GAMMA-1)*M**2/2) p = p0 * T_ratio**(GAMMA/(GAMMA-1)) rho = rho0 * T_ratio**(1/(GAMMA-1)) u = M * sqrt(GAMMA * p / rho) # 声速无量纲形式 E = p/((GAMMA-1)*rho) + 0.5*u**2 return np.stack([rho, rho*u, rho*E])这段代码里面积分布采用平方多项式过渡，收缩段和扩张段在端点处斜率不连续。如果希望面积导数连续，可以把平方改为三次 Hermite 多项式，这个改动对结果影响不大，但能避免喉部附近源项跃变导致的时间步长限制。
4.2 主循环与边界条件处理
主循环包含三件事：计算通量、施加边界条件、推进一个时间步。边界条件我给两个常见做法。亚声速入口：总温总压固定，入口马赫数由内部第一格点外插得到，然后通过等熵关系反算密度和压力。超声速出口：全部变量一阶外插。如果出口是亚声速，要给定出口背压，并让出口最后一个网格点的压力等于背压，密度和速度外插。下面是主循环的核心片段：
def macccormack_step(Q, F, S, dx, dt): Q_pred = Q.copy() # 预测步：前向差分 Q_pred[:, 1:-1] = Q[:, 1:-1] - (dt/dx)*(F[:, 2:] - F[:, 1:-1]) + dt*S[:, 1:-1] # 边界点留待边界条件处理 F_pred = flux(Q_pred, GAMMA) S_pred = source(Q_pred, dlnA_dx) Q_new = Q.copy() # 校正步：后向差分 Q_new[:, 1:-1] = 0.5*(Q[:, 1:-1] + Q_pred[:, 1:-1] - (dt/dx)*(F_pred[:, 1:-1] - F_pred[:, :-2]) + dt*S_pred[:, 1:-1]) return Q_new def flux(Q, gamma): rho = Q[0] u = Q[1] / rho p = (gamma-1) * (Q[2] - 0.5*rho*u**2) F = np.stack([Q[1], Q[1]*u + p, u*(Q[2]+p)]) return F def source(Q, dlnA_dx): rho = Q[0] u = Q[1] / rho p = (gamma-1) * (Q[2] - 0.5*rho*u**2) return np.stack([-Q[1]*dlnA_dx, -Q[1]*u*dlnA_dx, -(Q[2]+p)*u*dlnA_dx])计算通量时有个容易写错的细节：通量里那里的 p 必须用原始变量压力，而不是 P 的死亡。计算完 F 和 S 后，别忘了边界点 Q_pred[:, 0] 和 Q_pred[:, -1] 不会更新，必须在每步重新赋值，否则两端会积累错误并反射虚假波。
入口边界条件代码如下：
def apply_inflow_bc(Q, p0=1.0, rho0=1.0): # 用内部第一点估算入口马赫数 rho_in = Q[0, 1] u_in = Q[1, 1] / rho_in p_in = (GAMMA-1)*(Q[2, 1] - 0.5*rho_in*u_in**2) T_ratio = p_in / p0 # 等熵假设 M_in = sqrt(2/(GAMMA-1) * (T_ratio**(-(GAMMA-1)/GAMMA) - 1)) # 重新计算入口状态 T_in = 1.0 # 无量纲总温 p_in = p0 * (1 + (GAMMA-1)*M_in**2/2)**(-GAMMA/(GAMMA-1)) rho_in = rho0 * (1 + (GAMMA-1)*M_in**2/2)**(-1/(GAMMA-1)) u_in = M_in * sqrt(GAMMA*p_in/rho_in) E_in = p_in/((GAMMA-1)*rho_in) + 0.5*u_in**2 Q[:, 0] = [rho_in, rho_in*u_in, rho_in*E_in]这个边界条件的逻辑是：假设入口流动等熵，内部流场信息通过入口马赫数反馈到边界，使得入口流量可以随背压变化自动调整。如果直接用固定入口马赫数，那么背压变化时激波位置就永远错了。这也是很多初写版本只能算出设计工况、一算非设计工况就翻车的原因。
4.3 设计工况算例的典型计算结果
给定总压 p₀ = 101325 Pa、总温 T₀ = 300 K、喉部面积 1、入口面积 2、出口面积 2.5，程序推进到无量纲时间 t = 8 后，输出马赫数沿轴向分布。设计工况下，喉部马赫数应该为 1.0，出口马赫数约为 2.0。根据等熵面积比公式，出口面积比 2.5 对应马赫数约 2.16，这个数值和程序输出的偏差应该在 1% 以内。
检查收敛不能只看马赫数云图，我习惯输出几条曲线：每个时间步的最大残差、喉部马赫数和出口马赫数。残差定义为相邻两步守恒变量差值的范数，当残差降到初始值的 1e-6 以下且马赫数变化小于 1e-4 时认为收敛。如果残差卡在某个平台不降，多半是边界条件在反射波，这时需要把 CFL 调小到 0.2 再跑 2000 步，看是否恢复正常。
5. 常见问题排查与防错指南
5.1 入口出现负压力或负密度
现象：程序跑了几百步后，入口附近密度或压力变成负值，NaN 开始蔓延。
原因：大多数情况是入口边界条件给的是固定速度或固定马赫数，而不是特征边界。当背压变化导致入口流量需要调整时，固定马赫数会强行维持通量，产生非物理膨胀，密度直接被打穿。
解决：改成第 4 节中的特征型入口边界，先外插入口马赫数再用等熵关系重建热力学状态。如果问题依旧，检查初始场是否采用了等熵解，不要用常数场起步。
5.2 喉部马赫数始终达不到 1
现象：设计工况下喉部马赫数只到 0.8 左右就上不去了，流量看起来处于亚临界状态。
原因：面积分布函数喉部区域太粗糙，或者喉部正好落在网格点边界上。准一维模拟里喉部面积不光滑会导致等效面积比例偏大，实际面积比变小，马赫数当然上不去。
解决：先检查面积分布在喉部附近一阶导数是否连续，把面积分布函数加密到 1000 个点看马赫数是否变接近 1；如果变接近，说明是网格问题。另一个可能有争议的原因是时间步长太大导致数值耗散过高，把 CFL 从 0.6 降到 0.3 试试。
5.3 激波位置在扩张段来回摆动
现象：非设计工况下背压设置好后，激波不是钉在一个截面，而是在几个网格间振荡，马赫数曲线抖动。
原因：MacCormack 格式在激波附近需要数值粘性来抑制振荡，显式格式原本就缺乏迎风耗散，激波处产生过冲后会被对流到下游然后又反射回来。
解决：最常见的做法是加人工粘性项，在每一时间步对压力场做二阶平滑。具体实现时，在压力梯度大的区域加一个与二阶导数成正比的耗散项，系数取 0.01～0.05。如果坚持不加人工粘性，换成 Roe 格式，这个在下一章细讲。
5.4 残差一直不降，卡在 1e-3 左右
现象：出口马赫数看得过去，但残差曲线长期不动，像钉在了一个水平线上。
原因：边界条件与内点耦合不好，或者声波在入口和出口之间来回反射无法衰减。常见于出口用超声速外插条件而实际出口流场是亚声速的情况下，超声速外插条件会把下游压力信息错误地反射回去。
解决：先确认出口到底是不是超声速。用出口马赫数判断，如果小于 1，必须改用背压边界条件，即指定出口静压，内部变量外插。改完再跑 1000 步，残差通常会往下走两个量级。
5.5 面积导数在喉部产生尖峰导致发散
现象：程序初始化阶段就开始溢出，NaN 从喉部向两侧蔓延。
原因：d(lnA)/dx 计算时用了数值差分，如果面积分布在喉部附近变化剧烈或网格不均匀，数值导数会产生尖峰，源项驱动过度。
解决：不要用简单的 np.gradient 直接对离散面积差商，改用解析导数。若是分段多项式拟合，直接手写导函数；若是实验数据点，用样条拟合后对样条求导。源头问题解决后，CFL 可以适当加大而不会翻车。
6. 进阶：从准一维到背压扫描与格式对比
想把程序从验证工具变成工程工具，下一步做背压扫描。给定同一个几何，逐步改变出口背压 p_exit / p₀，从 0.9 到 0.1，观察激波位置、出口马赫数和流量变化。背压从高到低，流场经历四种状态：全亚声速无激波、扩张段内出现正激波、激波移到出口截面、全超声速无激波。实际数值模拟时，每算完一个背压就把上一个背压的收敛解当作新工况的初始场，能显著加速收敛。
表格式记录各背压下的激波位置：
| 背压比 p_exit/p₀ | 激波位置 x/L | 出口马赫数 | 说明 |
|---|---|---|---|
| 0.9 | 无 | 0.36 | 全亚声速流动 |
| 0.6 | 0.78 | 0.85 | 激波在扩张段上游 |
| 0.4 | 0.86 | 1.30 | 激波接近出口 |
| 0.1 | 出口边缘 | 2.16 | 达到设计工况 |
激波位置可以用正激波的 Rankine-Hugoniot 关系校验。取激波前马赫数 M₁，根据正激波关系算激波后马赫数 M₂ 和总压损失，然后和数值结果对比，偏差如果超过 1%，就要回头看网格分辨率了。
格式对比方面，把 MacCormick 格式换成二阶 Roe 或 Godunov 格式，最大的变化在激波处的过冲量。MacCormick 在 M₁ = 1.8 的激波前会产生 3%～5% 的密度过冲，Roe 格式自然迎风，过冲几乎为零。代价是 Roe 格式需要做特征分解，代码量多一倍，但换来的稳定性值得。
最后说一个自己踩过坑的习惯。程序跑通后，我的做法是先故意把面积比改成 1.0，即等截面管道，确认计算结果是等熵绝热管流而非喷管解，再改回真实几何。这一步花不了两分钟，却能一次性排查掉面积项符号错误、面积导数计算错误和边界条件写反三种常见问题。做完这个自检再算新工况，心里才算有底。希望帮到你。
本文还有配套的精品资源，点击获取
