非守恒型方程+MacCormack格式：拟一维等熵喷管CFD从零实现

1. 从一根拉瓦尔喷管的物理图像说起

如果你接触过火箭发动机、超声速风洞或者喷气推进相关的课题，大概率绕不开一个经典模型：拉瓦尔喷管。它的几何特征很简单——先收缩、后扩张，中间有个面积最小的喉部。但就是这么一个看起来平平无奇的管道，却能把气体从亚声速一路加速到超声速，背后靠的正是可压缩流动里最核心的一条规律：面积-速度关系。

很多人在学这部分内容时，公式推导都能跟上，但一到自己动手写代码算，就发现结果要么发散、要么在喉部附近出现莫名其妙的振荡，要么算出来的马赫数分布跟理论对不上。问题往往不在物理理解，而在数值方法的选择和实现细节上。这篇内容就是围绕这个痛点展开的：用非守恒型控制方程配合MacCormack显式格式，把拟一维等熵喷管流动从零算出来。

先把场景说清楚。所谓"拟一维"，是指我们假设流动参数只沿轴向变化，每个截面上参数均匀，忽略横向梯度。这个假设在喷管扩张角不大的时候相当好用，能把三维问题降成一维问题，计算量骤降，同时保留了可压缩流动最本质的特征。所谓"等熵"，是指忽略黏性和激波，流动过程可逆绝热。这套假设组合起来，就是CFD入门里最经典的验证算例之一。

那为什么偏偏要用非守恒型方程？现在主流CFD不都强调守恒型格式吗？这个问题很关键，我在后面会专门用一节来讲清楚它的来龙去脉。简单说，非守恒型在拟一维等熵问题里有它独特的便利性，尤其是配合MacCormack这种预测-校正格式时，实现起来非常直接，而且能拿到解析解做对照，是理解数值格式行为的最佳练手对象。

这篇文章适合谁看？如果你已经学过流体力学基础，知道什么是马赫数、总温总压、等熵关系式，并且想亲手把CFD代码跑起来，那这篇就是为你写的。如果你只是想了解数值方法的大致思路，也能从里面的原理拆解和踩坑记录里拿到不少东西。我会把每一步的"为什么"讲透，而不是甩一堆公式让你自己悟。

2. 非守恒型方程到底解决的是什么问题

2.1 守恒型与非守恒型的本质区别

要理解非守恒型，得先明白守恒型长什么样。守恒型控制方程的一般形式是：

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

其中U是守恒变量向量，F是对流通量向量，S是源项。它的特点是：通量项写在微分号内部，整个方程在离散后能保证物理量在控制体界面上严格守恒。有限体积法天然就是守恒的，这也是它在工程CFD里占据主导的原因。

非守恒型则是把方程写成所谓的"拟线性"形式：

∂W/∂t + A ∂W/∂x = S

这里W是原始变量（比如密度、速度、压力），A是系数矩阵。通量不再写在微分号里，而是系数矩阵乘以原始变量的空间导数。数学上这两种形式对光滑解是等价的，但在存在间断（比如激波）时，非守恒型会给出错误的结果——这也是它被主流CFD"嫌弃"的根本原因。

那为什么拟一维等熵喷管还要用它？因为这个算例里根本没有激波。等熵假设意味着流动全程光滑，非守恒型和守恒型的解在数学上完全一致。既然没有间断，非守恒型的"缺陷"就暴露不出来，反而能利用它形式简单、便于推导原始变量方程的优势。

2.2 拟一维等熵流动的原始变量方程推导

我们直接给出拟一维、非守恒、等熵假设下的控制方程组。设管道截面积为A(x)，流动参数为密度ρ、速度u、压力p、温度T。连续方程、动量方程和能量方程经过整理后，可以写成关于原始变量的形式。

连续方程：

∂ρ/∂t + ρ ∂u/∂x + u ∂ρ/∂x + ρu (1/A) dA/dx = 0

动量方程（无黏、无体积力）：

∂u/∂t + u ∂u/∂x + (1/ρ) ∂p/∂x = 0

能量方程在等熵条件下可以替换为等熵关系。对于量热完全气体，等熵关系给出：

p / ρ^γ = const

或者用声速表示更方便。声速 a = sqrt(γ p / ρ)，等熵条件下有：

a^2 = γ R T

以及总温总压守恒的关系式。实际编程时，最常用的做法是：用连续方程和动量方程推进密度和速度，然后用等熵关系反算压力和温度。这样只需要解两个偏微分方程，能量方程自动满足，计算量小且逻辑清晰。

这里有个细节值得强调：等熵关系是代数关系，不是微分方程。这意味着每一步时间推进后，我们直接用当前密度和某个参考状态算出压力，而不需要再积分一个能量方程。这是等熵假设带来的最大便利，也是这个算例适合入门的原因。

2.3 为什么这个场景下非守恒型反而更顺手

我实际写代码时对比过两种形式。守恒型需要处理通量向量，还要做通量分裂或者Riemann求解器，代码量大、调试麻烦。而非守恒型直接把原始变量的导数写出来，MacCormack格式一套就上，几十行核心代码就能跑。

更重要的是，边界条件的施加方式。非守恒型里我们直接操作密度、速度这些物理量，入口给总温总压对应的状态，出口做外推，非常直观。守恒型里你得先把原始变量转成守恒变量，施加边界后再转回来，多一层转换就多一层出错的可能。

当然，这不是说非守恒型更好。它只是在这个特定场景下更合适。一旦流动出现激波，或者你要算跨声速翼型这种有间断的问题，就必须老老实实回到守恒型。选方法要看问题，不能迷信某一种。

3. MacCormack预测-校正格式的运作机理

3.1 显式格式的时间推进逻辑

MacCormack方法是一种显式、二阶精度的有限差分格式，1969年由MacCormack提出，在CFD早期是绝对的主力。它的核心思想是"预测-校正"两步走：先用前向差分预测一个中间时刻的值，再用后向差分做校正，两步平均后得到二阶精度的结果。

对于形如 ∂W/∂t = -R(W) 的方程，其中R(W)代表空间离散后的右端项，MacCormack格式可以写成：

预测步（用前向差分）：

W_pred = W^n - Δt * R_forward(W^n)

校正步（用后向差分）：

W^{n+1} = 0.5 * (W^n + W_pred - Δt * R_backward(W_pred))

注意校正步里那个0.5的系数，它是对初始值和预测值做平均。这个平均是二阶精度的来源，也是格式稳定性的关键。

为什么前向和后向要交替用？因为单独用前向或后向差分只有一阶精度，而且会引入数值耗散或色散。交替使用后，误差项里的偶数阶项相互抵消，剩下的是二阶精度。这个思路和Lax-Wendroff格式有异曲同工之妙，但MacCormack不需要计算雅可比矩阵，实现上简单太多。

3.2 空间导数的离散细节

在拟一维问题里，空间导数包括 ∂ρ/∂x、∂u/∂x、∂p/∂x 这几项。用均匀网格的话，前向差分就是 (W[i+1] - W[i]) / Δx，后向差分就是 (W[i] - W[i-1]) / Δx。

但这里有个坑：面积变化项 (1/A) dA/dx 怎么处理。A(x)是已知的几何函数，dA/dx可以解析求导，也可以用差分。我建议用解析求导，因为几何函数通常很简单（比如双曲型或多项式型），解析导数精度更高，而且不引入额外的离散误差。

另一个细节是源项的离散位置。面积变化项本质上是个源项，它应该和通量项在同一时间层上计算。在预测步用n时刻的A和dA/dx，在校正步用预测步的A和dA/dx。如果偷懒全用n时刻的，会降到一阶精度，这个坑我踩过。

3.3 稳定性条件与时间步长选取

显式格式的命门是CFL条件。对于MacCormack格式，稳定条件大致是：

Δt ≤ CFL * Δx / (|u| + a)

其中CFL数一般取0.5到0.9之间。取1.0理论上临界，但实际因为有非线性效应和源项，很容易发散，所以我一般从0.5开始试，稳定后再往上调。

这里有个经验：喷管喉部附近是稳定性最脆弱的地方。因为那里速度接近声速，|u|+a 最大，而且面积变化率也大。如果全局用统一时间步，喉部会限制整个计算。我通常用全局最小时间步，或者做局部时间步（定常问题加速收敛时用）。对于非定常问题，必须用全局统一时间步，否则时间精度就没了。

还有一点，初始条件不能太离谱。如果你直接给一个均匀流场然后突然启动，初始瞬间会产生强烈的数值振荡，可能直接把格式冲垮。稳妥的做法是先用一阶耗散大的格式跑几百步，等流场平滑了再切换到MacCormack。或者干脆从解析解附近启动，收敛会快很多。

4. 从几何到边界：把问题搭起来的完整流程

4.1 喷管几何与网格生成

先定义喷管形状。经典的算例是喉部在中间，两侧用光滑曲线过渡。常用的几何函数是：

A(x) = A_throat * (1 + 某个关于x的光滑函数)

比如用双曲型：A(x) = A_throat + (A_in - A_throat) * ((x - x_throat)/(x_in - x_throat))^2 之类的形式。具体选什么不重要，关键是喉部面积最小、两侧单调变化、一阶导数连续。

网格方面，均匀网格最简单，但喉部附近梯度大，均匀网格要么全局加密浪费算力，要么喉部精度不够。我推荐用代数拉伸网格，在喉部附近加密。比如用tanh函数做坐标变换，把物理空间映射到计算空间的均匀网格上。这样既保证喉部分辨率，又不至于全局网格数爆炸。

网格数取多少？我试过200到1000个点。200个点能跑出大致趋势，但喉部马赫数会有几个百分点的误差。500个点基本收敛，1000个点用于验证网格无关性。对于学习目的，500个点足够了。

4.2 入口出口边界条件的处理

边界条件是CFD里最容易出错的地方，没有之一。拟一维等熵喷管通常分两种情况：亚声速入口、超声速出口，或者亚声速入口、亚声速出口（后者需要背压匹配）。

对于亚声速入口，我们给定总温T0和总压p0，以及流动方向。入口速度由内点外推或者用等熵关系反算。具体做法是：入口密度和速度用内点线性外推，然后用总焓守恒和等熵关系算出入口的声速和压力。这样入口边界是"软"的，允许数值波传出。

对于超声速出口，所有信息都从内向外传，所以直接零阶外推即可：W[出口] = W[出口-1]。超声速出口不需要给定任何物理量，给了反而会反射数值波。

对于亚声速出口，情况复杂一些，需要给定背压。出口压力固定为背压，密度和速度用特征线关系或者简单外推。这个场景容易出问题，因为亚声速出口有信息从外向内传，处理不好会反射。

我踩过的一个坑：入口边界用固定值而不是外推。一开始我图省事，入口密度速度全固定，结果算出来的流场在入口附近有非物理的边界层。后来改成外推+等熵反算，流场立刻干净了。这个教训是：边界条件要反映物理上的信息传播方向，不能想当然。

4.3 初始流场的构造技巧

初始条件直接影响收敛速度和稳定性。最粗暴的是全场给静止状态，但这会激发强烈的启动波。稍微好一点的是给一个线性分布的近似解。

我的做法是：用等熵面积-马赫数关系算一个近似解析解作为初始场。虽然这个解析解需要迭代求解（因为面积-马赫数关系是隐式的），但迭代几轮就能得到相当好的近似。用这个作为初始条件，MacCormack格式几乎不会振荡，收敛非常平滑。

面积-马赫数关系是：

A/A* = (1/M) * [(2/(γ+1)) * (1 + (γ-1)/2 * M^2)]^((γ+1)/(2(γ-1)))

给定A/A*，用牛顿迭代解出M，然后反算密度、速度、压力。这个初始场虽然不是精确解，但已经非常接近，数值格式只需要做小幅修正。

5. 代码实现中的关键细节与实测记录

5.1 主循环的骨架结构

核心代码结构大致如下（用Python示意，其他语言同理）：

# 初始化 W = initialize_flow() for n in range(num_steps): # 预测步 R_fwd = compute_residual(W, direction='forward') W_pred = W - dt * R_fwd # 施加边界条件到W_pred W_pred = apply_bc(W_pred) # 校正步 R_bwd = compute_residual(W_pred, direction='backward') W_new = 0.5 * (W + W_pred - dt * R_bwd) # 施加边界条件 W_new = apply_bc(W_new) W = W_new # 可选：每若干步输出或检查残差

注意边界条件要在预测步和校正步之后都施加。只在最后施加会导致中间状态越界，格式不稳定。这个细节很多教程不讲，但实际写代码时非常关键。

5.2 残差计算与面积源项

残差R(W)的计算是核心。以连续方程为例：

def compute_residual(W, direction): rho, u, p = unpack(W) if direction == 'forward': drho_dx = (rho[2:] - rho[1:-1]) / dx du_dx = (u[2:] - u[1:-1]) / dx else: drho_dx = (rho[1:-1] - rho[:-2]) / dx du_dx = (u[1:-1] - u[:-2]) / dx # 面积源项 A_face = 0.5 * (A[1:-1] + A[2:]) # 或类似平均 dA_dx = (A[2:] - A[1:-1]) / dx source = rho[1:-1] * u[1:-1] * dA_dx / A_face R_continuity = -(rho[1:-1] * du_dx + u[1:-1] * drho_dx + source) # 动量方程类似 ... return R

这里有个容易忽略的点：面积源项里的A用哪个位置的值。理论上应该用控制体界面的面积，实践中用相邻网格点的平均值即可。如果用网格中心值，会引入一阶误差。我实测过，用平均值和用中心值，在500点网格上马赫数分布差异在0.1%量级，但收敛速度有差别。

5.3 数值振荡的识别与抑制

MacCormack格式在光滑区表现很好，但在大梯度区（比如喉部附近）可能出现数值振荡。振荡的表现是：马赫数分布上出现高频的锯齿状波动，振幅随时间不衰减。

识别振荡很简单：把马赫数沿x画出来，如果曲线不光滑、有周期性锯齿，就是振荡了。抑制方法有几种：

• 降低CFL数：最直接，但收敛变慢。
• 加入人工黏性：在残差里加一个二阶或四阶耗散项，能有效压制振荡，但会抹平真实梯度。
• 使用通量限制器：把MacCormack改造成TVD格式，但这就偏离了原始格式的简洁性。
• 改用守恒型：如果振荡严重到影响结果，说明非守恒型在这个网格下不够用，换守恒型是正解。

我的经验是：500点网格、CFL=0.7、不加人工黏性，喉部振荡在可接受范围内。如果网格降到200点，振荡就明显了。所以网格分辨率是抑制振荡的第一道防线。

5.4 收敛判据与结果验证

定常问题的收敛判据通常是残差范数降到某个阈值以下，比如1e-6。但残差范数有时候会"假收敛"——残差降下去了，但解还没稳定。更可靠的做法是监控某个物理量（比如喉部马赫数或总流量）的变化率，当它连续几百步不变时，认为收敛。

验证结果最直接的办法是和解析解对比。拟一维等熵喷管有解析解（通过面积-马赫数关系得到），把数值解和解析解画在一起，看最大误差和分布形状。我实测下来，500点网格、MacCormack格式，马赫数分布的最大误差在0.5%以内，喉部附近稍大，约1%。这个精度对于学习目的是完全够的。

还有一个验证手段是网格收敛性分析：用200、400、800点分别算，看误差是否按二阶收敛。如果误差下降速率接近2，说明格式的二阶精度得到了体现。如果只有一阶，那说明某个地方实现错了（比如源项离散位置不对）。

6. 那些文档里不会写的踩坑实录

6.1 时间步长算错导致的"神秘发散"

我最开始写的时候，时间步长用了全局统一的 Δt = CFL * Δx / (|u_max| + a_max)，其中u_max和a_max是初始场的最大值。结果跑了没几百步就发散了。排查了半天，发现问题是：初始场是静止的，u_max≈0，但启动后速度迅速增大，实际需要的Δt比初始估计的小得多。

解决办法是每步都重新计算最大特征速度，动态调整Δt。或者更稳妥，用全局最大声速（对应总温）加上一个估计的最大速度来算。我后来改成每步计算 max(|u|+a)，发散问题就消失了。这个坑的教训是：显式格式的稳定性条件是局部的，不能用初始状态一劳永逸。

6.2 边界条件反射引发的非物理波动

另一个坑是出口边界。我一开始用超声速出口，直接零阶外推，没问题。后来想试试亚声速出口，给定背压，结果出口附近出现了强烈的反射波，整个流场都在振荡。

原因是我用了硬性给定出口压力，而出口速度用内点外推。这在亚声速出口是不自洽的——压力波会从出口反射回来。正确的做法是用特征线边界条件：根据出口马赫数判断传入和传出的特征波，传入波用外部条件，传出波用内点外推。这个实现起来复杂一些，但能消除反射。

对于入门学习，我建议先用超声速出口把流程跑通，等熟悉了再挑战亚声速出口。超声速出口简单且稳定，是验证格式的好选择。

6.3 面积源项的符号错误

这个坑很隐蔽。面积源项在连续方程里是 +ρu(1/A)dA/dx，在动量方程里是 +u^2(1/A)dA/dx（对于非守恒型）。符号取决于dA/dx的正负和方程的具体形式。

我一开始把连续方程的源项符号搞反了，结果算出来的密度分布在扩张段反而增大，完全违反物理。排查时我用了极限检验：假设一个等截面管道（dA/dx=0），源项应该为零，方程退化为标准一维流动。如果源项符号错了，等截面情况也会出错。用这个检验能快速定位符号问题。

6.4 网格拉伸带来的度量项遗漏

用拉伸网格时，坐标变换会引入度量项（雅可比）。如果直接在物理空间做差分，网格不均匀会导致差分公式的精度下降。正确做法是在计算空间做差分，然后乘以度量项。

我一开始偷懒，直接在物理空间用非均匀网格的差分公式，结果精度只有一阶，而且喉部误差特别大。后来改成计算空间均匀网格+度量项，精度立刻恢复到二阶。这个坑的教训是：网格拉伸不是免费的，度量项必须正确处理。

7. 从能跑到跑好：精度与效率的调优经验

7.1 网格分辨率与计算成本的平衡

前面提到500点网格够用，但具体多少点取决于你要多高的精度。我做过一组测试：

可以看到，网格加密一倍，误差降为约1/4（二阶收敛），但计算成本增加不止一倍（因为时间步也变小了）。所以存在一个性价比最优点。对于学习目的，500点是甜点区。

7.2 CFL数的选取策略

CFL数不是越大越好。理论上MacCormack的稳定极限在CFL=1附近，但实际因为非线性效应，CFL=0.9就可能发散。我一般用0.7作为默认值，稳定且收敛速度可接受。

如果想加速收敛，可以用当地时间步长：每个网格点用自己局部的Δt，而不是全局最小Δt。这样在特征速度小的地方可以用大步长，整体收敛速度能提升2到3倍。但注意，当地时间步长只适用于定常问题，非定常问题必须用全局统一时间步。

7.3 残差监控与收敛加速

残差监控是判断收敛的核心手段。我通常监控三个量：连续方程残差的L2范数、动量方程残差的L2范数、以及总流量的变化率。当三者都降到阈值以下且不再下降时，认为收敛。

如果收敛太慢，可以尝试残差光滑或者多重网格。但对于这个简单算例，这些高级技术有点杀鸡用牛刀。更实用的加速手段是：用解析解作为初始场，能把收敛步数减少一半以上。

7.4 结果后处理与可视化要点

算完之后，可视化是验证结果的重要环节。我通常画三张图：马赫数沿轴分布、压力沿轴分布、以及数值解与解析解的对比。

画图时注意：横坐标用x/L归一化，纵坐标马赫数范围取0到可能的最大值。如果出口是超声速，马赫数可能到2以上。对比图里，解析解用实线，数值解用离散点，这样能清楚看到误差分布。

还有一个技巧：画残差收敛曲线（对数坐标）。如果残差曲线是直线下降，说明收敛健康；如果出现平台或振荡，说明有问题。这个曲线是调试时的第一手资料。

8. 这套方法能走多远：适用边界与扩展方向

MacCormack加非守恒型这套组合，在拟一维等熵喷管上表现很好，但它有明确的适用边界。一旦流动出现激波、分离、或者强非等熵效应，这套方法就不够用了。比如超声速进气道里的激波串、喷管过膨胀状态下的激波，都需要守恒型格式配合激波捕捉技术。

如果想从这个算例往更复杂的问题扩展，我建议的路径是：先把非守恒型换成守恒型，保持MacCormack格式不变，看看结果差异。然后加入能量方程，去掉等熵假设。再然后加入黏性项，变成Navier-Stokes。每一步都只改一个因素，这样能清楚看到每个假设的影响。

另一个扩展方向是非定常问题。拟一维非定常喷管流动可以用来研究启动过程、阀门调节等瞬态现象。这时候时间精度变得重要，MacCormack的二阶时间精度就有价值了。但要注意，非定常问题必须用全局统一时间步，而且边界条件的处理要更加小心。

还有一个有意思的方向是与实验数据对比。虽然我这里没有实验条件，但如果你能找到公开的喷管实验数据（比如某型风洞的壁面压力分布），可以把数值结果和实验对比，看看等熵假设在什么范围内成立。这种对比能让你对模型的适用性有更直观的认识。

我个人在实际操作中的体会是：这个算例最大的价值不是算出一个喷管流场，而是让你亲手体会数值格式的每一个细节如何影响结果。从时间步长到边界条件，从源项离散到网格拉伸，每一个选择都会在结果里留下痕迹。把这些问题都踩一遍，你对CFD的理解就不再停留在公式层面，而是真正有了"手感"。这种手感，是看多少篇论文都换不来的。