企业官网建设流程全解析

简介：拟一维喷管流动的 MATLAB 数值仿真程序包，面向 CFD 初学者及流体力学课程学习者。程序利用 MacCormack 预测-校正格式求解简化后的喷管流动控制方程，可在低计算成本下模拟速度、压力、密度沿轴向的变化规律，适合用于理解可压缩流场时间推进与差分格式的构建。包内仅包含 1 个 .m 脚本文件，整体压缩后约 1KB，代码量精简、结构清晰，便于直接阅读和调试运行。目前已有 904 人学习下载。通过运行该脚本，读者可掌握从网格划分、初始与边界条件设置、预测步与校正步迭代到结果可视化的一整套 MATLAB 实现思路，也可在此基础上修改入口参数、扩展为含激波等更复杂的算例，是拟一维喷管数值模拟入门的实用参考。

1. 拟一维喷管与MacCormack方法：为什么一段变截面管道就能说清喷管流动的本质

很多刚接触可压缩流动的人，第一次接触“拟一维喷管”这个概念时都会愣一下：喷管明明是三维几何，怎么就“一维”了？这里的“拟一维”指的是流动参数只沿轴向变化，但截面积沿轴向可以缩、可以扩，通过面积变化把真实的二维收缩-扩张效应“投影”到一条线上。配合MacCormack方法，用几十行代码就能把亚音速气流加速到超音速、甚至在扩张段捕捉到激波。我见过不少做课程设计或入门CFD的人在这一步会卡住，不是因为公式看不明白，而是不知道边界条件怎么给、CFL取多大、为什么算着算着就NaN了。这篇笔记按我实际调试的顺序，把方程、代码、参数和踩坑一次说清，适合正在写第一个可压缩流求解器的人照着复现。

MacCormack方法在喷管问题里地位特殊：它是显式、二阶精度、预测-校正型格式，不需要解大型线性方程组，编程量极小，却能正确捕捉激波和膨胀波，所以几乎所有CFD教材都会拿它做入门范例。但“能跑通”和“能跑对”是两码事。接下来我会从控制方程开始，逐步落到一段可复现的Python代码，再把我在调试中遇到的反直觉问题列成排查清单。

2. 控制方程与预测-校正格式：先把守恒型方程和两步格式对齐

2.1 拟一维喷管流动的守恒型方程

拟一维喷管的控制方程是一组带面积源项的双曲型守恒方程。对无粘、无热传导、完全气体的可压缩流动，采用守恒型变量的形式是唯一稳妥的选择——尤其在后文要捕捉激波的时候，守恒型能保证激波位置正确。很常见的一个错误是直接对速度、压力写非守恒型方程，激波附近会出现非物理的振荡，而且质量流量沿管长守不住。我一般推荐直接从守恒型出发，五个字：密度、动量、总能量。

方程组的无量纲形式我直接给出：

∂Q/∂t + ∂F/∂x = S

其中：

Q = [ρ, ρu, e]ᵀ F = [ρu, ρu² + p, (e + p)u]ᵀ S = -(1/A)·(dA/dx)·[ρu, ρu², (e + p)u]ᵀ

这里ρ是密度，u是速度，p是压力，e是单位体积总能量，A(x)是喷管截面积。完全气体状态方程把压力和总能连起来：p = (γ - 1)(e - ½ρu²)，γ取1.4。面积源项是拟一维方程和普通一维欧拉方程的唯一区别，它代表管道面积变化对流动的“几何驱动”作用。在喉部附近dA/dx变化最快，源项贡献最大，这也是喉部流动最难算准的原因之一。

无量纲化的好处是让方程数值量级均匀，避免压力几十万帕而密度只有一点几带来的收敛困难。我用喉部声速和喉部密度做特征量进行归一化，这样初始场可以给得很粗糙，算几步后自动调整到合理范围。面积函数我用了教材里最经典的对称双曲形喷管：

A(x) = 1 + 2.2*(x - 1.5)² （x在0到3之间，喉部在x=1.5处）

喉部面积恒为1，入口面积约5.95，面积比足够让气流在出口达到马赫数2以上。读者可以拿这个几何做基准，之后换成自己设计的型面只需要改这一行。

2.2 MacCormack预测-校正步的离散骨架

MacCormack格式是“预测步用一侧差分，校正步用另一侧差分”的典型代表。对守恒方程做有限体积式的半离散近似，空间导数先向前差分得到预测值，再用向后差分做校正，最后把两个时间层的空间导数取平均。这样做时间精度是二阶，空间精度也是二阶，而实现成本只比一阶迎风多一次通量计算。

预测步的离散形式是：

Q_pred[i] = Q_n[i] - (dt/dx)·(F_n[i+1] - F_n[i]) + dt·S_n[i]

校正步是：

Q_corr[i] = Q_n[i] - (dt/(2dx))·(F_pred[i] - F_pred[i-1] + F_n[i+1] - F_n[i]) + (dt/2)·(S_pred[i] + S_n[i])

这里的下标i对应网格单元中心，F_n、S_n是用本时间层Q算出的通量和源项，F_pred、S_pred是用预测步Q_pred重新计算的通量和源项。注意校正步里空间导数是“预测步后差 + 本层前差”的平均，方向和时间层都做了交叉，这是MacCormack格式能拿到二阶精度的关键。如果你把校正步的F_pred与F_n写反位置，格式实际上退化成某种耗散型格式，精度降到一阶还可能出现奇怪的相位误差。

稳定性方面，MacCormack格式属于显式格式，受CFL条件约束。对一维流动，时间步长必须满足：

dt = CFL · dx / (|u| + a)_max

其中a是当地声速，CFL一般取0.5到0.8。这里有个容易忽略的坑：(|u| + a)_max必须扫描整个网格取最大值，而不能用入口值或喉部值代替。拟一维喷管从亚音速到超音速，特征速度变化很大，如果固定某一个参考值算dt，很可能在出口附近违反CFL条件导致发散。

3. 边界条件与网格：入口总压、出口静压、喉部加密怎么配

3.1 边界条件的物理设定与数值外推

拟一维喷管最常见的工况是：入口接高压气源（总压给定），出口排向环境（背压给定）。这种情况下入口用亚音速入流，需要给两个物理条件（总压、总温），出口用亚音速出流，需要给一个物理条件（静压），其余变量用特征线或外推获得。如果出口流动已经是超音速，则所有变量都从内部外推，出口静压信息向上游传播不过去。

我在离散代码里习惯把边界节点单独处理，用“虚拟单元法”而不是“单边差分法”。入口虚拟单元的值这样构造：驻点条件（总压p0、总温T0）给定，从内部第一格取马赫数，按等熵关系反推入口静压、静温、密度和速度，再把守恒变量填进虚拟单元。这样做的效果是入口像一个“声学吸收边界”，反射波很小。如果直接给固定密度、速度，入口会反复反射数值波，经常出现上游振荡污染全场。

出口虚拟单元分成两种情况：亚音速出流时给定出口静压p_exit，从内部最后一格外推动量、总能量，再结合状态方程求密度；超音速出流时全部一阶外推。判断依据是最后一个网格单元的马赫数是否超过1。在同一个代码里同时写这两个分支，就能复现从完全亚音速到带激波的超音速流动的完整物理图像。

我见过很多人试图在入口直接给声速或马赫数，这等于把本该由计算确定的量人为固定了，结果入口段的压力分布和解析解对不上。记住一点：亚音速边界只能“指定”沿特征线传入的物理量，其余必须靠内部外推，这是双曲型方程数学性质决定的，不是编程习惯问题。

3.2 网格划分与面积函数选择

网格密度对MacCormack方法的影响比很多人想象的大。这个格式的数值耗散比较小，激波捕捉时能压出比较锋利的间断，但代价是对网格敏感性高——网格太粗时激波前后会出现明显的过冲和下冲，网格太密又会放大计算代价。对拟一维喷管，我通常先试200个网格，如果激波附近振荡减不下来再加密到400。工程经验是：MacCormack格式的激波过渡区在3到5个网格宽度，如果振荡超过峰值10%就需要处理。

网格分布最简单用均匀网格就能工作，但喉部附近的面积变化梯度大，最好在喉部局部加密。我常用的办法是对实际坐标做一次余弦变换，让网格点向喉部集中。均匀网格在喉部只有2到3个点，压力、马赫数的峰值会偏低，而加密后喉部能有5到8个点，解就更接近解析值。

面积函数的选择决定了喷管的工况。我做代码调试时用的是上述对称双曲型面积，简单且解析解存在，适合验证程序正确性。实际工程要模拟真实型面时，一般把面积离散点做成表格，在求解器里用线性插值或三次样条计算每个网格中心的A和dA/dx。这里注意：dA/dx必须用解析导数或样条导数，如果用相邻点面积差除以步长，二阶精度会被破坏，激波位置可能偏半个网格。

4. 从零写一个可复现的MacCormack求解器：代码与参数对照

4.1 最小Python实现

下面给出一段我能直接跑通的最小实现，用NumPy做数组运算，不需要任何其他第三方库。它把上一章说到的守恒型方程、MacCormack预测-校正步、特征边界条件全部落在代码里，核心时间推进循环只有二十多行。

import numpy as np import matplotlib.pyplot as plt gamma = 1.4 nx = 201 # 网格数（含边界虚拟单元） x = np.linspace(0, 3, nx) dx = x[1] - x[0] # 面积函数：喉部在x=1.5，面积比约5.95 A = 1.0 + 2.2 * (x - 1.5)**2 dA = 4.4 * (x - 1.5) # 解析导数 # 初始场：全流道给一个低速亚音速均匀流 rho = 0.8 + 0.0 * x u = 0.2 + 0.0 * x p = 1.0 / gamma * np.ones_like(x) # 由原始变量转守恒变量 e = p / (gamma - 1.0) + 0.5 * rho * u**2 def primitive(Q): rho = Q[0] mom = Q[1] e = Q[2] u = mom / rho p = (gamma - 1.0) * (e - 0.5 * rho * u**2) return rho, u, p def flux(Q): rho, u, p = primitive(Q) mom = Q[1] e = Q[2] F = np.zeros_like(Q) F[0] = mom F[1] = mom * u + p F[2] = (e + p) * u return F def source(Q, idx): # idx是网格索引数组，返回对应位置的源项S rho, u, p = primitive(Q) mom = Q[1] e = Q[2] S = np.zeros_like(Q) factor = -dA[idx] / A[idx] S[0] = factor * mom S[1] = factor * mom * u S[2] = factor * (e + p) * u return S # 边界条件：入口总压、总温；出口静压 p0_in = 1.0 / gamma # 无量纲总压 T0_in = 1.0 # 无量纲总温 p_exit = 0.75 / gamma # 出口背压，可通过改变该值切换工况 def apply_bc(Q): # 入口：用等熵关系从内部马赫数反推 rho_int, u_int, p_int = primitive(Q[:, 1]) Ma_in = u_int / np.sqrt(gamma * p_int / rho_int) p_in = p0_in / (1 + 0.5*(gamma-1)*Ma_in**2)**(gamma/(gamma-1)) T_in = T0_in / (1 + 0.5*(gamma-1)*Ma_in**2) rho_in = p_in / T_in # 无量纲完全气体 u_in = Ma_in * np.sqrt(gamma * T_in) e_in = p_in/(gamma-1) + 0.5*rho_in*u_in**2 Q[:, 0] = [rho_in, rho_in*u_in, e_in] # 出口：亚音速给静压，超音速全部外推 rho_out, u_out, p_out = primitive(Q[:, -2]) Ma_out = u_out / np.sqrt(gamma * p_out / rho_out) if Ma_out < 1.0: Q[:, -1] = Q[:, -2].copy() e_out = Q[2, -1] rho_out = (gamma-1)*e_out / (p_exit*gamma/(rho_out*u_out**2) + 0.5*(gamma-1)*u_out**2) if False else rho_out # 更稳妥：固定压力，密度/动量外推，再修正能量 p_fixed = p_exit rho_new = p_fixed / (gamma - 1) / (Q[2,-2]/Q[0,-2] - 0.5*(Q[1,-2]/Q[0,-2])**2) Q[:, -1] = [rho_new, Q[1,-2]*rho_new/Q[0,-2], p_fixed/(gamma-1)+0.5*rho_new*(Q[1,-2]/Q[0,-2])**2] else: Q[:, -1] = Q[:, -2].copy() return Q CFL = 0.6 nt = 30000 Q = np.vstack([rho, rho*u, e]).copy() Q = apply_bc(Q) for n in range(nt): # 计算全局稳定时间步 rho_all, u_all, p_all = primitive(Q) a_all = np.sqrt(gamma * p_all / rho_all) dt = CFL * dx / np.max(np.abs(u_all) + a_all) Q_pred = Q.copy() # 预测步：前向差分 F = flux(Q) S = source(Q, np.arange(nx)) Q_pred[:, 1:-1] = Q[:, 1:-1] - dt/dx * (F[:, 2:] - F[:, 1:-1]) + dt * S[:, 1:-1] # 校正步：后向差分（混合差分） F_pred = flux(Q_pred) S_pred = source(Q_pred, np.arange(nx)) 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]) + 0.5*dt*(S[:, 1:-1] + S_pred[:, 1:-1])) Q = Q_new.copy() Q = apply_bc(Q) # 后处理：转回原始变量绘图 rho, u, p = primitive(Q) Ma = u / np.sqrt(gamma * p / rho)

这段代码里有几个设计决策值得单独说明。首先是预测步和校正步的差分方向故意相反，这是MacCormack格式的核心；其次时间步长dt每步都重新计算，扫描全场所有单元的特征速度最大值，保证任何位置都不违反CFL条件；最后边界条件在每步末尾调用，确保数值解一直满足物理约束。如果你把出口背压p_exit改到0.75/γ，得到的解是典型的带激波工况——流动在喉部达到音速后在扩张段加速到超音速，再通过一道激波跳回亚音速，出口压力与背压匹配。

4.2 关键参数（CFL、迭代步数、网格数）怎么取

CFL参数是第一个要调的。对MacCormack格式，我试过从0.9开始，预测步直接出现负密度；降到0.5左右稳定，0.6到0.7之间能得到最小的数值振荡。不要追求极限CFL，这个格式不像隐式格式那样能硬撑，0.7以上就开始出现高频噪声。你可以在代码里把CFL设成0.9跑一次，观察密度残差随迭代步数变化，大概率看到指数增长——这就是稳定性边界。

迭代步数要分两个阶段。第一阶段是瞬态传播阶段，入口和出口边界条件的信息要穿过整个流道，至少需要L/((u+a)max)数量级的时间，换算成步数通常是几千步。第二阶段是收敛阶段，残差会缓慢下降，可能出现周期性震荡，这是激波在网格间“跳跃”造成的，不属于发散。我会固定跑30000步再输出，因为MacCormack格式收敛速度就是慢，10000步时压力分布看着差不多，但质量流量误差可能还在1%量级。

网格数从201起步，加密到401时激波位置移动通常在0.5个网格以内，说明网格已经基本收敛。再往上加密收益不大，反而开始暴露MacCormack格式在高频段耗散不足的问题，表现为出口附近出现小尺度的数值噪声。这个现象在网格从201到401时基本消失，但在801时重新出现，所以日常调试用201网格足够。

5. 常见问题排查与避坑：发散、振荡、质量流量不守恒

5.1 现象：CFL超过0.9直接NaN

我在调这个求解器时最早遇到的就是这个问题，把CFL从0.6改到0.9，跑了不到500步，密度直接变成负数，随后所有变量全部变成NaN。第一反应以为是边界条件写错了，反复检查没有发现问题。后来把逐步的dt打印出来才发现，喉部附近dA/dx很大的区域，源项贡献特别大，而且出口马赫数超过1之后，特征速度一下从1.6跳到2.8，全局dt被出口限制得很小，但校正步里用了固定dt，失稳就从出口附近开始。

解决方法是全局统一dt，并且把CFL降到0.6。有人习惯在每个单元用当地dt，那样虽然也能跑，但会破坏激波的传播方向性，造成虚假的相移。另外，给代码加一层保护：每步检查rho的最小值，如果出现负密度就自动把步长减半重试本步。这个保护在调试期很值，能帮你快速定位是哪一步导致的发散。

5.2 现象：激波处出现锯齿振荡

完全亚音速工况下程序一切正常，一旦出口背压调到激波工况，激波前后各3个网格单元内出现明显的锯齿状振荡，压力峰值比解析解高出15%以上。最初怀疑是激波两侧差分方向选择不当，后来查资料发现这是MacCormack格式的固有属性——中心型预测-校正格式没有迎风耗散，在间断附近必然产生吉布斯式振荡。

解决的办法不是改格式，而是加人工粘性项。在预测步之前对守恒变量做一次二阶人工耗散：

# 人工粘性项，系数епсилон通常取0.05~0.2 eps = 0.1 Q[:, 1:-1] += eps * (Q[:, 2:] - 2*Q[:, 1:-1] + Q[:, :-2]) * 0.5

这个操作会略微抹平激波过渡区，但能让压力振荡从15%降到2%以内。ε系数是有讲究的：太小压不住振荡，太大激波被抹宽到8个网格以上。我一般从0.1起步，观察激波前后最大过冲量，在过冲小于2%的前提下选尽可能小的ε。这个技巧本质上是在MacCormack格式的截断误差里注入一点二阶耗散，不会降低总体的二阶精度，因为激波附近本来就不是光滑区域。

5.3 现象：出口背压设错，整个流场解漂移

我第一次试着复现不同工况时，把出口背压直接从0.75/γ改到0.4/γ，结果入口的总压等熵关系不成立了，喉部前出现大面积的亚音速减速区，流动甚至出现回流。这是因为背压比低于临界值时，喷管处于过膨胀状态，正常数值解会出现出口斜激波，但拟一维模型没有横向自由度，激波无法存在，解只能靠剧烈的数值耗散强行“适应”错误的背压条件。

这个问题的正确做法是先算临界背压比。等熵关系给出喉部达到声速时的出口压力参考值，再根据正激波关系算出激波后的压力，只有当背压落在“激波后压力”和“完全超音速出口压力”之间时，激波才会稳定存在于扩张段内部。我的经验是先把背压设在较高值（比如0.85/γ）确保全流道亚音速，再逐步降低，每次降5%，观察激波出现位置。如果你想直接算纯超音速出口，出口马赫数超过1后边界分支会自动切到全外推，但你得把背压设得非常低，低到接近真空，这个物理上等价于排气进入太空。

5.4 现象：质量流量沿程不守恒

稳态解出来后，我把ρuA沿x轴画出来，发现喉部前后质量流量差达到3%——按守恒方程理论应该是完全平坦的直线。排查了压力、马赫数分布，都和解析解吻合，唯独流量不守恒。问题出在时间推进没有完全收敛：MacCormack格式的收敛判据是残差降到机器精度，而我跑了20000步就停止了，此时瞬态扰动还没完全排干净。

后来我把迭代步数从20000加到50000，质量流量误差从3%降到0.3%，再加密网格到401点，误差降到0.1%以下。另外还发现dA/dx如果用了数值差分（A[i+1]-A[i-1]）/2dx，流量误差会额外增加0.5%左右，换用解析导数后立刻改善。这提醒我：拟一维问题里几何导数的影响比想象中敏感，面积函数给得不光滑，整条流量线都会波动。

5.5 现象：网格粗时入口总温漂移

某次调试时发现，入口总温在计算过程中从1.0漂移到了1.03，虽然变化不大，但影响激波位置。原因是我入口虚拟单元里用等熵关系反推静温时，引入了内部网格马赫数的误差，而马赫数误差来自于边界附近网格的高频振荡。解决措施是把入口虚拟单元推到离物理边界更远一些的位置，让等熵反推使用流道内部第二个网格的值而不是紧挨边界的第一格，相当于给边界条件做了一次二阶精度的空间滤波。

这个问题比较隐蔽，因为它不影响稳定性和激波捕捉，只体现在长时间的漂移上。如果你的代码需要跑很多个工况做对比，建议每5000步输出一次入口总温，如果漂移超过1%，就去检查边界附近网格的振荡幅值，而不是急着改格式。人工粘性对这个现象也有帮助，ε取0.1时入口总温漂移几乎消失。

6. 验证与进阶：用等熵关系做基准，再往二维喷管走

6.1 解析解校验：Ma数与压力分布对照

写完求解器第一件事不是调参，而是验证。完全亚音速工况下，拟一维喷管流动存在等熵关系闭式解：给定入口总压、总温和喉部面积，任意截面上的马赫数满足面积-马赫数关系式，再通过等熵关系求出压力、密度和温度。把代码跑出来的Ma数分布和这个解析解画在一起，误差应该在1%以内；如果超出，优先排查边界条件而不是网格。

激波工况的验证稍微复杂一点，但没有解析解可用时需要做网格收敛性检查：画同一工况在201、401、801三个网格数下的压力分布，看激波位置是否随网格加密单调趋于一个固定值。MacCormack格式的激波位置误差会随网格加密按二阶收敛，如果加密后位置跳来跳去，说明人工粘性系数没调好或者出口背压正好落在临界值附近。这一步做完，程序才算真正可信。

6.2 从一维到二维：面积项换成几何边界

把求解器扩展到二维轴对称喷管时，MacCormack方法依然可用，但有一个决定性的区别：一维的面积源项被替换成几何边界条件，网格变成二维结构化网格，时间步长受两个方向的CFL条件共同约束，实际能取的CFL会从0.6降到0.3左右。我在扩展时发现，一维代码里练熟的边界处理经验——入口给总压总温、出口按马赫数判断外推、激波附近加人工粘性——在二维里依然适用，只是每步要处理角点和壁面反射。

进阶方向还有两类值得尝试：一是把时间推进从MacCormack换成Runge-Kutta显式格式，稳定性边界几乎不变但代码更容易加自适应时间步；二是引入限制器把人工粘性改为CUSP或HLLC通量，彻底解决激波振荡问题，但代价是代码量翻倍。我的习惯是先用MacCormack理解可压缩流的数值行为，等一维调通后再开始碰二维网格和Riemann求解器。

最后说一个习惯：所有参数选项做成配置文件，无论是CFL、网格数、背压比还是人工粘性系数，都从外部读入，求解器本体只做纯计算。这样当你凑齐一组能跑出漂亮激波的参数后，不会因为一次误改代码而丢掉结果——我的后悔药就是Git提交没做好，某次重构时把整个调试记录冲掉了，之后所有实验都先commit再动手。希望这篇笔记能帮你在拟一维喷管上少走几段弯路，参数调不出来的时候，回头检查边界条件和dA/dx，大概率问题出在那里。

本文还有配套的精品资源，点击获取