简介拉瓦尔喷管准一维流动数值模拟程序是一份面向计算流体力学学习者的经典案例资源对应《计算流体力学基础及应用》中的典型问题。程序以准一维假设为基础用有限差分法求解收缩—扩张喷管内的绝热流动方程完整覆盖控制方程离散、时间推进、边界条件处理和人工粘性设置等关键环节可帮助读者掌握可压缩流动数值模拟的基本流程。压缩包共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(Ep/ρ)] S [-ρu d(lnA)/dx, -ρu² d(lnA)/dx, -(ρEp)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[i1] - 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^{n1} Q_corr。这组式子有两点关键第一F_pred 必须用预测步得到的守恒变量重新计算不能用上一时刻的通量第二源项 S_pred 同理要用预测后的变量算否则校正步的时间精度掉回一阶。很多初写版本只在预测步更新源项校正步沿用旧源项最后收敛结果虽然差不多但激波位置会偏移半个网格。3.2 网格、时间步长与 CFL 条件的具体选择网格划分直接影响激波分辨率和收敛速度。我一般先在喉部附近加密因为喉部马赫数要穿过 1梯度变化最剧烈。均匀网格要加密到 5001000 个点才能看清喉部压力分布如果做工程快速评估300 个点配合喉部局部加密也够用。加密方式可以取双曲正切变换也可以简单地令网格点在喉部附近按余弦分布加密计算量差别很小关键是别把喉部放在两个网格点中间否则马赫数 1 的位置会来回跳动。时间步长用 CFL 条件约束。对绝热无粘管流最大特征波速是 |u| a所以时间步长取dt CFL * dx / max(|u| a)显式格式 CFL 上限理论值是 1但 MacCormack 格式有中间步稳定性限制实际中 CFL 取 0.30.6 比较稳妥。注意这里的 a 是当地声速不是入口参考声速必须在每个网格点扫描一遍再取全局最大值。总计算时间取多少通常取无量纲时间 t 510即整个喷管流场被扫过 510 个来回边界信息传到全场并衰减完。用 500 个网格、CFL0.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, L1.0, area_in2.0, area_throat1.0, area_exit2.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_in0.3, p01.0, rho01.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/(GAMMA1)) * (1 (GAMMA-1)*M_guess**2/2))**((GAMMA1)/(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, p01.0, rho01.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.010.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.60.780.85激波在扩张段上游0.40.861.30激波接近出口0.1出口边缘2.16达到设计工况激波位置可以用正激波的 Rankine-Hugoniot 关系校验。取激波前马赫数 M₁根据正激波关系算激波后马赫数 M₂ 和总压损失然后和数值结果对比偏差如果超过 1%就要回头看网格分辨率了。格式对比方面把 MacCormick 格式换成二阶 Roe 或 Godunov 格式最大的变化在激波处的过冲量。MacCormick 在 M₁ 1.8 的激波前会产生 3%5% 的密度过冲Roe 格式自然迎风过冲几乎为零。代价是 Roe 格式需要做特征分解代码量多一倍但换来的稳定性值得。最后说一个自己踩过坑的习惯。程序跑通后我的做法是先故意把面积比改成 1.0即等截面管道确认计算结果是等熵绝热管流而非喷管解再改回真实几何。这一步花不了两分钟却能一次性排查掉面积项符号错误、面积导数计算错误和边界条件写反三种常见问题。做完这个自检再算新工况心里才算有底。希望帮到你。本文还有配套的精品资源点击获取