简介面向流体力学数值仿真学习者的MATLAB程序资源围绕拟一维喷管模型演示MacCormack预测-校正格式在可压缩流动求解中的完整实现。内容紧扣喷管理论与有限差分方法帮助读者理解从网格划分、边界条件设置到时间推进、收敛判断的代码化流程适合具备一定MATLAB基础并希望上手CFD编程的本科高年级学生或工程初学者。压缩包仅含1个m文件大小约1KB代码紧凑、无冗余依赖便于逐行研读和修改参数进行扩展实验。已有904人学习浏览程序可直接在MATLAB中运行通过绘图输出观察速度、压力、密度随空间的演化可作为入门一维可压缩流模拟的参考模板。1. 拟一维喷管碰上 MacCormack 方法一节课代码换来整个超声速流场拟一维喷管算例加 MacCormack 方法是 CFD 入门里最典型的「小模型、大物理」组合。它用一个带面积变化的管道把可压缩效应、激波捕捉、显式时间推进一次性全演示完代码量却只有百行上下。很多第一次接触的人以为把网格画好、套个格式就能跑结果频繁卡在 NaN 爆炸、喉道马赫数到不了 1、出口反射波把全场污染这些问题上——问题恰恰出在对格式的理解和对初始/边界的处理上。这套资源是一份完整可跑的拟一维喷管 MacCormack 求解器覆盖守恒型方程离散、预测-校正时间步进、亚声速入口/超声速出口边界处理、收敛监控与后处理。适合三类人做课程设计、准备复试上机、或者想用最短代码验证自己对可压缩流动理解的从业者。拿到手第一步是看面积分布函数怎么写第二步是把边界条件换回你自己算例的压比条件剩下的就是调 CFL 和耐心等残差下降。2. 守恒型方程与 MacCormack 格式先搞懂离散再写循环2.1 为什么要把面积变化塞进源项而不是改成几何网格拟一维喷管的控制方程是带源项的一维欧拉方程形式上比纯一维多了面积这一项。关键在于喷管的截面积变化不能靠「把网格画成收缩再扩张」来体现因为方程本身是沿轴向的一维坐标几何收缩必须显式地出现在方程里否则动量方程里不会有驱动流动加速的力。守恒形式写成∂U/∂t ∂F/∂x S其中U [ρA, ρuA, ρEA]ᵀ F [ρuA, (ρu²p)A, (ρEp)uA]ᵀ S [0, p·dA/dx, 0]ᵀ这里的 A 是当地截面积E 是单位质量总能即 E p/((γ-1)ρ) u²/2γ 1.4。源项第三行是零因为面积变化只对动量方程产生推力效应能量方程里压力做功已经包含在通量的 (ρEp)uA 里了。很多初学版本把面积塞进通量而不是源项算出来喉道两侧压力分布是错的。判断自己写没写对的简单办法取恒定面积段源项必须自动退化为零方程回到标准一维欧拉方程。非守恒形式用原始变量 ρ、u、p 写在这个问题里很不推荐。喷管扩张段一旦出现正激波非守恒形式在激波处的数值解会产生额外的伪源项导致激波位置偏移几个网格守恒形式把 Rankine-Hugoniot 关系天然包含在离散里激波捕捉是「顺带完成」的不需要特殊处理。提示写代码时把 ρA、ρuA、ρEA 当作三个独立变量存而不是存 ρ、u、p。激波处密度和马赫数跳变很大但守恒变量在通量层保持守恒这是 MacCormack 这类差分格式能捕捉激波的前提。2.2 预测-校正两步差分方向交替的物理含义MacCormack 方法是二阶精度显式格式核心是预测步和校正步各用一次单侧差分方向相反然后取平均。预测步前差U_pred_i U_iⁿ − (Δt/Δx)(F_{i1}ⁿ − F_iⁿ) Δt·S_i校正步后差U_iⁿ⁺¹ 0.5[U_iⁿ U_iⁿ⁺¹_corr] U_iⁿ⁺¹_corr U_iⁿ − (Δt/Δx)(F_pred_i − F_pred_{i−1}) Δt·S_pred预测步的差分公式含入误差的截断项是二阶校正步再反向截一次两项相加后一阶误差互相抵消整体达到空间二阶精度。这个「前差后差再平均」的思路等价于对时间做了一次改进欧拉Heun 法推进所以时间方向也是二阶的。很多学员自己实现的时候容易犯一个误区预测步用前差校正步也想当然用前差。这会导致空间误差变回一阶激波被抹平两三个网格喉道附近的马赫数峰值明显偏低。判断格式是否写对看预测步用 F[i1] 还是 F[i−1]预测必须是前差校正必须是后差顺序不能反。这个格式还有个特性它没有显式的人工耗散项因此它的色散误差会在激波两侧产生小幅振荡这是格式内禀行为不是代码 bug。后处理时不要对着这些锯齿做任何滤波去抹掉物理量——先确认振幅相对于激波跳跃是否在 5% 量级以内。2.3 CFL 条件与声速特征显式步长的「最短命」逻辑显式时间步进的稳定性由 CFL 条件约束物理上信号传播的最大速度是 |u| a流动速度叠加声速一个时间步内信息不能跑过一个网格间距所以Δt CFL · Δx / max(|u| a)CFL ≤ 1实际使用时 CFL 取 0.8 比较稳妥。取 1.0 时在喉道附近那里速度梯度最大容易触发振荡取 0.5 以下虽然安全但收敛慢喷管算例要跑几万步才到稳态没必要。这里有个刚入门容易忽略的细节Δt 要在整个网格上取最小值用全域最大波速算而不是每个网格各自的时间步。因为显式格式的推进是同步的每个网格用各自的 Δt 会造成相邻网格状态错位相当于引入额外的非物理相位差。声速公式 a √(γp/ρ) 在每一步都要从原始变量重新算。初始场里如果 p 或 ρ 出现非正数a 直接变 NaN整个程序瞬间崩。所以主循环里第一步必须做正系数检查任何网格出现 p ≤ 0 或 ρ ≤ 0立刻停止并输出当前步数而不是让它带着负压力继续跑。提示CFL 的「玄学」在于它不是越大越好。MacCormack 在 CFL 0.9 以上时激波前后的过冲会明显加大如果只追求稳态解不想盯瞬态0.60.8 是性价比最高的区间。3. 代码落地网格、变量换算、边界与主循环3.1 网格与面积分布喉道处导数要光滑网格用均匀网格就够了拟一维问题的主要矛盾在面积函数 A(x) 的光滑性不在网格分布。经典算例用抛物线型面积分布A(x) 1 2.2·(x − 1.5)²喉道位于 x 1.5喉道面积 A* 1管道总长 L 3出口面积约 5.95。这个函数的好处是导数连续且解析式简单dA/dx 4.4(x − 1.5)不会在喉道出现导数跳变。import numpy as np gamma 1.4 L 3.0 nx 201 x np.linspace(0.0, L, nx) dx x[1] - x[0] # 面积分布喉道处 A 1两端 A ≈ 5.95 A 1.0 2.2 * (x - 1.5) ** 2 # 面积导数直接用解析式不要用 np.gradient dA_dx 4.4 * (x - 1.5)面积导数为什么必须用解析式因为 np.gradient 这类数值差分在喉道附近的误差会被源项放大尤其是喉道两侧导数方向相反一个网格的差异就可能让喉道马赫数停在 0.98 上不去。解析导数一行搞定别偷懒。如果想把算例改成自己的几何注意两点一是 A(x) 必须单调收缩再单调扩张喉道处导数过零二是面积函数要保证 C¹ 连续否则源项里 p·dA/dx 会突变激波位置会偏移。3.2 三套变量怎么换算守恒、原始、通量程序里需要维护三套量守恒变量 U、原始变量ρ, u, p、通量 F。它们的换算贯穿整个循环写错任何一个索引都会让结果悄悄出错。def cons_to_prim(U): rho U[0] / A u U[1] / U[0] p (gamma - 1.0) * (U[2] - 0.5 * U[1] * u) / A return rho, u, p def flux(U): rho, u, p cons_to_prim(U) F np.zeros_like(U) F[0] U[1] F[1] (rho * u * u p) * A F[2] u * (U[2] p * A) return F def source(U): rho, u, p cons_to_prim(U) S np.zeros_like(U) S[1] p * dA_dx return Scons_to_prim 里最容易被抄错的是压力换算。守恒变量第三行是 ρEA即 ρA·E其中 E p/((γ−1)ρ) u²/2。所以 ρEA pA/(γ−1) 0.5·ρu²A而 0.5·ρu²A 0.5·U[1]·u。于是 p (γ−1)(U[2] − 0.5·U[1]·u)/A。通量第三行写成 u·(U[2] pA) 而不是重新算一遍 ρE省一次计算也避免舍入误差叠加。这几个函数每个时间步会被调用两次预测步一次、校正步一次所以性能敏感点就在这保持精简是对的。3.3 入口滞止与出口外推两种最省事的边界写法边界条件是喷管算例里最容易砸锅的地方但用对套路后其实很简单入口是亚声速出口是超声速两边的处理逻辑完全不一样。入口如果按「固定三变量」写死等于把入口速度也钉住喉道永远无法壅塞。正确做法是速度 u 从内部第一个网格外推出来密度和压力用滞止关系重建。这样入口可以根据内部流动状态自行调整马赫数。def apply_bc(U): # 入口亚声速给定滞止参数速度从内点外推 u_in U[1, 1] / U[0, 1] a2 1.0 - 0.5 * (gamma - 1.0) * u_in * u_in M_in2 u_in * u_in / a2 fac 1.0 0.5 * (gamma - 1.0) * M_in2 rho_in rho0 / fac ** (1.0 / (gamma - 1.0)) p_in p0 / fac ** (gamma / (gamma - 1.0)) e_in p_in / ((gamma - 1.0) * rho_in) 0.5 * u_in * u_in U[0, 0] rho_in * A[0] U[1, 0] rho_in * u_in * A[0] U[2, 0] rho_in * e_in * A[0] # 出口超声速流动零梯度外推 U[:, -1] U[:, -2] return U这里 rho0 和 p0 是无量纲化的滞止值取 rho0 1.0p0 1/gamma这样喉道声速恰好为 1方便核对马赫数。入口速度外推的写法在稳态时完全够用如果做瞬态精度要求高的研究可以换特征线边界但喷管稳态算例没必要。出口零梯度外推的前提是出口马赫数大于 1。设计压比下这个条件自动满足但如果背压调高出现正激波出口会变成亚声速这时必须换成固定静压边界。3.4 主时间循环预测、校正、残差监控主循环结构是固定的先做边界 → 算通量源项 → 算 Δt → 预测步 → 再做一次边界 → 校正步 → 收敛判断。def maccormack_step(U): U apply_bc(U) F flux(U) S source(U) rho, u, p cons_to_prim(U) a np.sqrt(gamma * p / rho) dt cfl * dx / np.max(np.abs(u) a) U_pred np.zeros_like(U) for i in range(1, nx - 1): U_pred[:, i] U[:, i] - dt / dx * (F[:, i 1] - F[:, i]) dt * S[:, i] U_pred apply_bc(U_pred) F_pred flux(U_pred) S_pred source(U_pred) U_new np.zeros_like(U) for i in range(1, nx - 1): U_corr U[:, i] - dt / dx * (F_pred[:, i] - F_pred[:, i - 1]) dt * S_pred[:, i] U_new[:, i] 0.5 * (U[:, i] U_corr) return apply_bc(U_new), dt循环体里每个细节都有讲究预测步用 F[i1] − F[i]校正步用 F_pred[i] − F_pred[i−1]方向必须相反校正步的源项 S_pred 要用预测步的守恒变量重算而不是沿用预测步前的源项预测步结束后必须先做边界再算通量否则边界网格的通量是上一步的旧值。初始场给全场滞止条件加一个微小正速度u0 0.01相当于让流动从静止启动U np.zeros((3, nx)) for i in range(nx): rho, u, p 1.0, 0.01, 1.0 / gamma e p / ((gamma - 1.0) * rho) 0.5 * u * u U[0, i] rho * A[i] U[1, i] rho * u * A[i] U[2, i] rho * e * A[i] res_history [] for step in range(20000): U_old U.copy() U, dt maccormack_step(U) res np.sqrt(np.mean((U - U_old) ** 2)) / np.mean(U) if step % 1000 0: res_history.append((step, res)) if res 1e-10: break print(converged at step, step, residual, res)初始速度给 0 会引发入口回流给 0.01 这个量级既触发流动又不至于让入口边界在第一步就把速度外推成负值。残差用 L2 范数除以全场平均守恒变量归一化后大概在 1e-2 量级起步降到 1e-10 需要几千步。nx 201 时通常 30006000 步收敛如果想看激波细节建议 nx 401。提示不要每步都打印残差IO 会拖慢循环每 5001000 步记录一次即可。收敛后立即停止比硬跑满 20000 步省时且便于复现。4. 两个标准算例全超声速喷管与管内正激波4.1 算例一设计压比下跑到稳态设计压比对应出口静压约 0.095·p0。用零梯度出口边界跑出来的就是这条等熵设计工况喉道马赫数精确到 1扩张段单调加速到出口马赫数约 2.19。跑完后处理代码rho, u, p cons_to_prim(U) a np.sqrt(gamma * p / rho) M u / a # 输出喉道和出口的关键值 throat np.argmin(A) print(throat M , M[throat]) print(exit M , M[-1]) print(exit p/p0 , p[-1] / p0)喉道马赫数 0.9991.001 是合格的说明壅塞发生正确。出口马赫数 2.182.21 都在解析值附近压力比接近 0.095。这三个数任何一个偏了都说明代码有系统性错误而不是调参问题。压力分布从入口到喉道单调下降喉道之后继续下降马赫数在喉道前单调上升、喉道后继续上升全程无激波、无回流。整个流场等熵熵增应该为零——可以把 p/ρ^γ 沿程算一遍理想情况下是全常数实际会有千分位的小波动那就是 MacCormack 格式的数值耗散造成的。4.2 算例二提高背压逼出正激波当出口背压高于设计值但低于「全亚声速」临界值时扩张段会出现正激波。激波位置由背压唯一确定背压越高激波越往喉道方向移。把出口边界从零梯度换成固定静压def apply_bc_with_backpressure(U, p_b): # 出口亚声速密度和速度外推压力给定 rho_b U[0, -2] / A[-2] u_b U[1, -2] / U[0, -2] e_b p_b / ((gamma - 1.0) * rho_b) 0.5 * u_b * u_b U[0, -1] rho_b * A[-1] U[1, -1] rho_b * u_b * A[-1] U[2, -1] rho_b * e_b * A[-1] return U取 p_b 0.3·p0 时激波大约出现在 x ≈ 2.0 附近p_b 0.2·p0 时激波后移到 x ≈ 2.4。激波前后马赫数满足正常激波关系 M₂ ≈ 1/M₁ 的超声速支换算。观察要点激波前压力单调下降激波处压力突跳激波后到出口压力平台缓慢回升因为面积还在扩张亚声速段面积增大压力上升。马赫数曲线在激波处从约 1.7 突跳回约 0.6然后缓慢下降。4.3 结果怎么读马赫数曲线、压力比与流量守恒三组算例放一起对比规律非常清晰工况出口背压 p_b/p0流动特征出口马赫数工况 A≈ 0.095全超声速扩张≈ 2.19工况 B0.30扩张段正激波≈ 0.45工况 C≥ 0.50全亚声速 1工况 C 全程亚声速喉道马赫数小于 1流动不会被壅塞喉道处马赫数跟随背压变化。判段一个算例到底跑在哪个工况先看喉道马赫数是 1 说明壅塞不是 1 说明背压不够低或不够高。质量流量守恒是最硬核的自检方法。喷管稳态流量 m_dot ρuA 沿程应该是常数。把每个网格的 ρuA 算出来画一条线如果两端差超过 2%要么没收敛要么边界写错。常见错误是入口边界速度外推方向写反导致入口流量和出口流量系统性差 5% 以上。激波位置的细微差别还可以用激波前后压力比来验证上游马赫数 M₁ 对应的理论压力比 p₂/p₁ [2γM₁² − (γ−1)]/(γ1)把数值结果代进去偏差应该在 1% 内。超过这个量级说明激波附近网格太粗加密 nx 到 401 再试。5. 常见问题与排查NaN、反射波、喉道卡住的实战记录5.1 NaN 爆炸先查 CFL 和面积导数现象跑了几百步残差正常某一步突然所有变量变 NaN。原因分三种CFL 超过 1显式格式直接失稳面积导数用了数值差分喉道处导数不连续产生伪源项初始场某网格压力或密度出现非正值声速开出根号负数。解决先把 cfl 降到 0.8 再跑一遍dA_dx 换成解析式初始化循环里加一句断言确保 rho 和 p 全为正。另外建议在每个时间步后面加一个最小压力监测一旦 p ≤ 0 立即打印错误步数退出rho, u, p cons_to_prim(U) if np.min(p) 0 or np.min(rho) 0: print(non-physical state at step, step) break5.2 喉道马赫数卡在 0.98不是格式问题现象马赫数到喉道附近增速变缓最大只有 0.98 左右扩张段全是亚声速。原因大概率是背压比选错出口压力还没低到让喉道壅塞的程度。另一个常见原因是最外层循环步数给得少几千步就停了流动根本没到稳态。解决先核对工况。把出口背压降到设计值约 0.095·p0再看如果是零梯度出口还卡 0.98说明出口压力还没通过膨胀波释放这需要更多步数。MacCormack 的收敛速度在这个问题上是线性的什么加速收敛技巧都不用耐心跑够 5000 步再说。5.3 出口反射波污染全场现象喉道下游出现周期性传播的压力脉动像波浪一样从出口往上游跑稳态迟迟到不了。原因出口实际马赫数小于 1 却用了零梯度外推。超声速出口的信息全部向外传播零梯度外推正确但亚声速出口有一条特征线指向流场内必须有外部信息背压约束零梯度相当于把这个约束变成了「无反射」结果产生伪反射波。解决判断出口马赫数M 1 用零梯度外推M 1 就要用带背压的边界详见 4.2。更省事的做法是直接把出口背压按工况给定不要依赖零梯度。若瞬态要求高可以在出口前加一段 sponge 区阻尼系数从内到外渐增但稳态算例我一般不这么干。5.4 激波前后锯齿振荡格式天性现象激波位置对了但激波前后出现约 5% 幅度的振荡像 Gibbs 现象。原因MacCormack 是色散型格式没有内置耗散在激波这种间断处必然产生过冲。这不是 bug。很多人看到锯齿就加整体平滑结果梯度信息被抹掉激波位置反而偏了。解决先确认振幅相对激波跳跃小于 5%如果大于 10%加密网格nx 从 201 加到 401或者显式加入工耗散项。我一般只在动量方程里加 ε·(u_{i1} − 2u_i u_{i−1})ε 取 0.010.05加多了会把激波抹成三四个网格宽的斜波。5.5 残差停在 1e-4 不再下降现象残差曲线下降很快到 1e-4 量级变成平台再也下不去。原因激波处稳态本来就有小幅周期性振荡数值解在极限环上或者收敛判据里 U_old 和 U 的差值被激波那几个网格主导全局残差被「污染」。解决改收敛判据用喉道和出口马赫数每 500 步的相对变化来判断或者对残差做局域平均跳过激波前后 5 个网格。正规做法是看喉道马赫数连续 200 步变化小于 1e-6。别硬追求 1e-12MacCormack 格式在激波算例里达不到那么低。6. 验收与进阶用面积-马赫数关系给算例背书再谈两步改动跑通只是第一步验证才算数。拟一维等熵流动有一个严格的解析判据面积-马赫数关系。对任意截面积 A 和喉道面积 A*等熵流场任意位置满足 A/A* f(M)。把数值解的马赫数代进去反推面积比再跟真实的 A(x)/A* 对比误差就是整个求解器的「体检报告」。实现方法很简单对每个网格位置根据它处于喉道上游还是下游用二分法求面积-马赫数关系的对应根。def area_mach_inv(M): gg (gamma 1.0) / (2.0 * (gamma - 1.0)) return (1.0 / M) * ((2.0 / (gamma 1.0)) * (1.0 0.5 * (gamma - 1.0) * M * M)) ** gg def mach_from_area(ar, supersonic): lo, hi (1.001, 5.0) if supersonic else (0.05, 0.999) for _ in range(100): mid 0.5 * (lo hi) if (area_mach_inv(lo) - ar) * (area_mach_inv(mid) - ar) 0: hi mid else: lo mid return 0.5 * (lo hi) # 喉道前取亚声速支喉道后取超声速支 throat np.argmin(A) M_ana np.zeros_like(M) M_ana[:throat1] [mach_from_area(A[i], False) for i in range(throat1)] M_ana[throat1:] [mach_from_area(A[i], True) for i in range(throat1, nx)] max_err np.max(np.abs(M - M_ana)) print(max Mach error , max_err)nx 201 时 max_err 通常在 0.5%1% 量级加密到 401 后应该减半左右这就是网格收敛的直观证据。如果误差反而变大检查是不是喉道位置找错了——A 的最小值索引必须和 dA_dx 过零的位置严格一致。带激波的算例不能用这条关系验证激波下游那里已经不等熵只用它验证激波上游和喉道位置。往后走两步这个算例就能从「课程作业」升级成「能写进论文的方法验证」一是把固定几何的喷管改成喉道面积随时间变化的变几何算例加入移动边界二是把 MacCormack 替换成守恒更强的 TVD 型格式比如 Harten-Yee 的 TVD 修正对比同一背压下激波位置的差异你会直观看到高阶格式在激波分辨率上的优势。保守一点也能做把面积函数换成带喉道圆弧的真实喷管截面用这套代码跑一遍设计点输出的压力分布可以当作实验数据的验证基准。从那以后我每跑一个 MacCormack 算例都强制先做两件事画面积导数曲线确认几何光滑算完先对面积-马赫数解析解再谈流场细节。这套流程帮我挡掉过至少三次把边界写成「固定三变量」的翻车现场。希望帮到你——代码拿到手先把第 3 章的主循环跑通再回头改面积函数和背压你的第一个超声速流场就落地了。本文还有配套的精品资源点击获取