2-20特征线法MOC+Euler CFD喷嘴设计：从燃烧室参数出发，自动设计一个超声速二维喷嘴

一句话介绍： 这套程序先根据燃烧室压力、温度和环境压力寻找最大静推力对应的理想喷嘴出口状态，再用特征线法自动生成二维超声速喷嘴轮廓，随后建立计算网格，并通过二维 Euler 方程与 MacCormack 有限体积格式计算喷嘴内部的马赫数和压力分布。

先看它解决什么问题

超声速喷嘴不是简单把一根管子做成“越来越宽”。

喷嘴出口面积太小，气体没有充分膨胀；出口面积太大，又可能出现过膨胀甚至激波。真正合理的喷嘴需要同时回答：

喉部以后应该膨胀到多大？出口马赫数应该是多少？

喷嘴壁面应该怎样弯曲，才能让超声速流动平滑转向？

设计出来的轮廓放进二维流场计算以后，马赫数和压力是否符合预期？

这套程序把整个过程串成：

燃烧室 / 环境条件 → 最大推力出口状态 → 特征线法设计喷嘴轮廓 → 自动生成 CFD 网格 → Euler 方程求解 → 检查马赫数和压力分布

所以它既包含“喷嘴设计”，也包含“设计后的流场验证”。

1. 这个工程真正包含什么？

理解整个工程主要看 4 个文件。

2. 当前设计条件是什么？

源码默认参数为：

因此喉部面积为：

A = 0.025 × 0.1 = 0.0025 m²*

程序把喉部设为临界状态，也就是：M = 1

随后从喉部开始逐渐增大出口面积，寻找能够获得较大推力的超声速出口状态。

3. 第一步：为什么先搜索出口面积？

对于给定喉部面积和燃烧室状态，不同出口面积会对应不同的超声速出口马赫数。

程序使用等熵喷管的面积—马赫数关系：

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

其中：

A：当前出口面积；

A*：喉部临界面积；

M：当前马赫数；

γ：比热比。

程序并不是直接知道 M，而是每改变一次面积，就利用 Newton-Raphson 迭代反求对应的超声速马赫数。

这一步实际上建立了：出口面积 ↔ 出口马赫数之间的对应关系。

4. 出口压力怎样随马赫数变化？

求出马赫数以后，程序利用等熵关系计算静压：

P/P_c = [1 + (γ-1)M²/2]^[-γ/(γ-1)]

随着喷嘴继续膨胀、M 增大，静压会不断下降。

理想情况下，希望喷嘴出口压力与环境压力接近：

P_e ≈ P_amb

这样气体在喷嘴内部已经比较充分地完成膨胀，不需要在喷嘴外继续通过明显的压缩或膨胀结构重新适应环境压力。

5. 程序怎样判断哪个出口状态推力更大？

对于每个候选出口面积，程序都会计算静推力：

F = ṁV_e + (P_e - P_amb)A_e

其中：

ṁV_e：高速喷气带来的动量推力；

(P_e - P_amb)A_e：出口压力与环境压力差产生的压力推力。

所以喷嘴不能只追求“马赫数越大越好”。

出口面积继续增加以后：

喷气速度可能继续变化；出口压力继续下降；压力推力也会发生变化。

程序把这些影响一起算进去，然后从所有候选出口状态中寻找最大静推力。

6. 当前参数下的设计出口状态

按照 nozzle.m 的实际计算过程核对，最大静推力附近的设计参数约为：

可以看到出口静压已经非常接近：

P_amb = 101 kPa

这说明程序搜索出的最大推力设计点同时接近理想膨胀状态。

由于喷嘴宽度为 0.1 m，对应的目标全出口高度约为：

0.006025 / 0.1 = 0.06025 m

而喉部全高度只有 0.025 m，因此气流需要在喉后逐渐扩张到约 2.41 倍面积。

7. 为什么知道出口面积以后还不能直接绘制曲线？

如果只把喉部和出口用一条直线连接，虽然几何上得到了扩张通道，但并没有保证超声速流动在整个区域中满足合适的膨胀波关系。

超声速流动具有很强的方向传播特征。

壁面发生转折以后，会产生膨胀波；这些波在流场内部传播并相互作用，最终决定：局部马赫数；流动方向；喷嘴壁面的合理形状。

所以程序下一步采用：Method of Characteristics，特征线法（MOC）

来设计真正的二维超声速喷嘴轮廓。

8. 特征线法为什么能设计喷嘴？

超声速二维无旋流动中，可以定义 Prandtl-Meyer 膨胀函数 ν(M)。

程序使用的关系为：

ν(M) = √[(γ+1)/(γ-1)] × atan√[(γ-1)(M²-1)/(γ+1)] - atan√(M²-1)

这里 ν 表示气流从声速状态膨胀到马赫数 M 时所需要的总转角能力。

对于最短长度喷嘴，程序使用：θ_max = ν(M_e) / 2

当前：M_e ≈ 2.263

对应最大壁面转角约：θ_max ≈ 19.24°

也就是说，喷嘴在喉后先通过膨胀使流动逐渐转向，然后再通过特征线之间的相互作用把出口流动重新整理得接近平行。

9. C⁺ 和 C⁻ 特征线在程序中做什么？

程序利用两个特征不变量：

K⁺ = θ - ν

K⁻ = θ + ν

这里：

θ：局部流动方向角；

ν：Prandtl-Meyer 角。

沿不同族的特征线，这些组合量具有确定的传播关系。

程序用它们计算：

每个特征点的 θ；每个特征点的 ν；对应马赫数 M；马赫角 μ；特征线之间的交点位置。

其中马赫角为：μ = asin(1/M)

所以特征线法真正做的是：先求流动应该怎样转，再根据特征线斜率反推出这些流动状态在二维空间中应该出现在哪里。

最终这些特征点共同确定喷嘴壁面。

10. 为什么程序使用 15 条特征线？

当前设置：num = 15

也就是用 15 条主要特征线离散喷嘴膨胀区。

特征线数量越多：流场离散更细；喷嘴壁面点更多；

对连续膨胀过程的近似更精细；

但计算点和交点数量也会增加。

程序会把计算出的上半喷嘴轮廓保存到 noz，同时利用对称性得到下半喷嘴。

因此最终得到的是一个关于中心线对称的二维收敛后扩张喷嘴模型，其中本程序重点设计的是喉部以后的超声速扩张段。

11. 设计完轮廓以后，为什么还要重新生成网格？

特征线法得到的喷嘴壁面点在 x 方向并不是均匀分布。

CFD 求解更适合使用结构化计算网格，因此 noz_mesh.m 会：

找到特征线喷嘴轮廓中的较小 x 间距；建立均匀 x 分布；对喷嘴壁面进行线性插值；在中心线和壁面之间生成 y 方向节点。

最终形成：中心线 → 喷嘴上壁之间的一半结构化二维网格。

由于喷嘴上下对称，计算和显示时可以利用对称性扩展到完整喷嘴。

12. CFD 部分真正求解什么方程？

noz_cfd.m 和 solver.m 求解的是二维可压缩 Euler 方程。

守恒形式可以概括成：

∂Q/∂t + ∂F/∂x + ∂G/∂y = 0

程序中的守恒变量为：

Q = [ρ, ρu, ρv, e]

分别代表：

质量；x 方向动量；y 方向动量；总能量。

所以 CFD 部分不是只根据面积公式计算一条中心线马赫数，而是在整个二维喷嘴网格中推进：

密度 + 两个方向速度 + 压力 + 温度 + 马赫数

对应的流场。

13. MacCormack 方法怎样推进流场？

程序采用 MacCormack predictor-corrector 思想。

每个时间步主要经历：

当前守恒量 → 正向通量预测 → 施加边界条件 → 反向通量修正 → 得到下一时刻状态

这种方法通过前向和后向差分组合，提高时间推进的精度。

程序还自动根据 CFL 条件计算时间步长：

Δt ∝ CFL / 最大局部波传播速度

当前：CFL = 0.8

也就是说每一步时间推进都会根据当前速度、声速和网格大小调整稳定时间尺度，而不是固定使用任意时间步长。

14. CFD 为什么是对特征线设计的进一步验证？

在 nozzle.m 中，喷嘴设计主要依赖：等熵一维关系；

Prandtl-Meyer 膨胀；二维特征线几何。

而 CFD 部分重新从二维 Euler 守恒方程出发计算整个流场。

因此两部分关注的角度不同：

MOC：从目标出口状态反推合理喷嘴几何

CFD：给定喷嘴几何以后，重新求解内部二维流动

如果 CFD 得到的出口马赫数、压力变化与设计目标比较一致，就说明设计轮廓在二维流动计算中也具有较好的合理性。

这形成了一个完整闭环：

目标性能 → 几何设计 → 数值流场验证

15. 为什么不只用一维面积—马赫数关系设计喷嘴？

一维公式非常适合快速确定：

出口面积；出口马赫数；压力；理论推力。

但它把每一个截面上的流动看成近似均匀的一维状态。

它并不能直接告诉我们：

喉部到出口之间的壁面到底应该怎样弯曲。

MOC 补上了二维超声速波系和流向信息，而 CFD 又进一步计算实际二维流场。

因此这套程序不是重复计算，而是三个层次：

一维性能设计 → 二维几何设计 → 二维数值验证

16. 这种方法理论上为什么适合超声速喷嘴？

优势一：出口状态和推力直接关联

程序先从推力公式出发搜索出口面积，使喷嘴几何设计不是凭经验选一个扩张比，而是先确定一个与当前燃烧室和环境压力匹配的目标状态。

优势二：MOC 利用了超声速流动的传播规律

在超声速流中，扰动沿特征方向传播。

通过 C⁺、C⁻ 两族特征线，可以把流动角和马赫数变化转化为空间中的几何交点，因此特别适合设计二维超声速膨胀流道。

优势三：最短长度喷嘴兼顾膨胀和几何尺寸

程序采用 θ_max = ν(M_e)/2 的最短长度喷嘴设计思想，在达到目标出口马赫数的同时控制喷嘴轴向长度。

优势四：CFD 不再局限于中心线关系

Euler CFD 在完整二维网格上求解质量、动量和能量守恒，可以观察喷嘴内部局部马赫数和压力变化，从而给解析/特征线设计增加一层数值验证。

17. 这套程序适合用在哪里？

超声速喷嘴设计教学

可以完整观察：

面积—马赫数关系 → Prandtl-Meyer 膨胀 → 特征线网格 → 喷嘴壁面

的形成过程。

火箭或高压气体喷管基础设计

给定燃烧室压力、温度、气体性质和环境压力，可以研究不同工况对应的理想出口面积和出口马赫数。

MOC 与 CFD 联合实验

先用理论方法生成喷嘴，再用二维 Euler 方程验证，是非常典型的“设计 + 数值验证”教学流程。

喷嘴参数敏感性分析

改变燃烧室压力、环境压力或气体 γ，可以观察设计面积比、出口马赫数和喷嘴轮廓如何随工况变化。

18. 怎么运行？

这个工程的运行顺序非常重要。

第一步：运行 nozzle.m

它会自动完成：

出口面积搜索 → 最大推力状态 → 特征线计算 → 喷嘴轮廓设计

第二步：直接运行 noz_mesh.m

不要清空工作区。这个脚本会使用 nozzle.m 中已经得到的 noz，生成 CFD 网格。

第三步：直接运行 noz_cfd.m

同样不要清空变量。程序会使用刚才生成的 x、y 网格，进行默认 500 个时间步的 CFD 推进。

如果只做喷嘴几何设计，运行前两步即可；如果还希望验证二维流场，则继续执行第三步。

最值得调整的参数只有燃烧室压力/温度、环境压力和特征线数量 num：前两类决定设计工况，num 主要影响特征线离散精细程度。

19. 一句话看懂这个项目

这是一个 MATLAB 二维超声速喷嘴“性能设计 + 特征线几何设计 + CFD 验证”程序：它先通过等熵面积—马赫数关系和推力公式搜索最大静推力对应的出口状态，当前默认工况得到 A_e≈0.006025 m²、M_e≈2.263、推力约 3900.87 N；随后利用 Prandtl-Meyer 函数和 15 条特征线生成最短长度喷嘴轮廓，再自动建立结构化网格，并用二维 Euler 方程和 MacCormack 有限体积格式计算喷嘴内部的马赫数与压力场。