简介：这是一份基于MATLAB的特征线法喷管流动CFD计算源码，面向流体力学、计算流体力学方向的科研人员、工程师及高年级学生。喷管内部高速气流涉及可压缩性与非定常效应，特征线法通过追踪流场特征信息传播，对连续方程与动量方程进行离散迭代，尤其适合非均匀网格下的一维、二维流动求解；这一数值方法被封装为单个MATLAB脚本文件，使用者只需设置初始条件与边界条件，即可自动计算喷管内的速度、压力等分布，并以图形化方式查看结果。压缩包内仅包含1个m文件，整体大小约1KB，代码体积虽小但结构清晰，便于在MATLAB中直接运行、断点调试与二次开发，可针对不同喷管构型与工况灵活修改参数。目前已有583人学习下载，适合希望从理论走向代码实现、快速上手CFD喷管数值模拟的读者。
1. 不需要求解器也能算喷管流场：这条技术路线很多人没见过
给定一个收缩-扩张喷管的壁面几何，想知道马赫数在轴线上怎么走、出口能不能达到设计马赫数，最常见做法是开 Fluent 或 OpenFOAM 跑一遍密度基求解器。但有一个更轻的算法：在做特征线法（Method of Characteristics, MOC）之前，先意识到超音速喷管流场是一个初值问题而非边界值问题。流场里每一点的信息只来自上游某个扇形区域，因此可以从一条已知初值线出发，沿马赫波一步步“走”出整个流场——不需要迭代、不需要解大型线性方程组，更不需要人工粘性。trysome.m 这套 MATLAB CFD 源码做的就是这件事：用特征线法求解二维轴对称喷管内的定常无粘超音速流场。它解决的是喷管气动设计中最常用的一类问题：已知壁面型面，求设计工况下的马赫数场、压力场和流动角。适合气体动力学课程设计、风洞喷管预研，或者想用最小代价验证自己写的 CFD 求解器的人。下面的内容会把理论、可运行的源码骨架、验证方法和最容易翻车的细节一次讲完。
2. 特征线法怎么把喷管流场变成“沿线积分”
2.1 从速度位方程到马赫波：超音速才有的“单行道”
二维定常无粘等熵流动可以用速度位方程描述。设 x 方向马赫数分量为 Mx、y 方向为 My，速度位方程写成：
(1 - Mx²)φxx - 2MxMyφxy + (1 - My²)φyy = 0
这是一个二阶拟线性偏微分方程。用判别式 B² - 4AC 判断类型：当 M < 1 时方程为椭圆型，扰动向全空间传播，必须给定整个边界上的条件；当 M > 1 时方程为双曲型，扰动只能在下游一个有限锥形区域内传播。这个锥形的边界就是马赫波，马赫角满足 μ = arcsin(1/M)。
物理含义很清楚：亚音速流动中，下游的扰动会逆流影响上游，所以必须全局联立求解；超音速流动中信息沿特征线单方向传播，流场任何一个下游点只依赖上游的初始数据和边界。这个区别决定了求解策略完全不同——超音速区可以做空间推进，而不是时间推进或全场迭代。特征线法就是利用这个“单行道”属性，把偏微分方程沿特征方向拆成常微分方程，然后逐点推进。
2.2 C+ 与 C-：黎曼不变量把两个方程拆成四个常微分方程
在二维无旋流动中，特征线有两个方向，分别记为 C+ 和 C-，其斜率为：
C+：dy/dx = tan(θ + μ) C-：dy/dx = tan(θ - μ)
其中 θ 是流动方向角。沿这两条特征线，物理量满足黎曼不变量条件：
沿 C+：θ + ν(M) = 常数 沿 C-：θ - ν(M) = 常数
这里 ν 是普朗特-迈耶角，是马赫数的单调函数：
ν(M) = sqrt((γ+1)/(γ-1)) · arctan(sqrt((γ-1)(M²-1)/(γ+1))) - arctan(sqrt(M²-1))
这个公式是特征线法源码里出现频率最高的表达式。它的作用是：在特征线上流动角 θ 和 ν 的和或差保持不变，于是只要知道特征线一端的状态，另一端的状态就通过代数关系直接确定——不需要解微分方程。这正是 MOC 计算效率高的原因。
| 符号 | 含义 | 典型初值 |
|---|---|---|
| M | 马赫数 | 初值线上 1.01 ~ 1.2 |
| θ | 流动方向角，单位弧度 | 轴线上 0，壁面处最大 |
| μ | 马赫角，asin(1/M) | 随 M 减小 |
| ν | 普朗特-迈耶角 | M=1 时为 0 |
| γ | 比热比 | 空气取 1.4 |
注意 ν 的反函数没有解析式，需要做数值反解。后面源码里的 inversePM 函数就是干这个的，而且它做了两件事：一是用牛顿法迭代，二是把 dν/dM 的解析导数写进去保证收敛稳定。
2.3 为什么这套 MATLAB CFD 源码用 MOC 而不是有限体积推进
喷管设计工况的流场本质是一簇膨胀波，不包含激波。膨胀波本身就是马赫波的包络，用特征线法天然贴合物理过程，壁面上的边界条件可以直接嵌入特征线交点计算中，不需要像有限体积法那样构造通量函数和人工粘性。MOC 的网格也不是传统意义上的结构化网格，而是由特征线交织成的曲线网格，网格点就是特征线交点，每个点的值由上游两个点直接算出，没有隐式耦合。
更关键的是精度可控。特征线法的误差来源主要是几何线性化和数值插值，而不是格式耗散，因此可以精确捕捉膨胀波扇和壁面拐角影响区。对工程预研来说，MOC 给出的是解的特征结构，而不是被数值耗散抹平的平均场。当然边界条件是有限制的：只能处理定常无粘流，设计工况下无激波。如果喷管过膨胀或欠膨胀，出现内激波，特征线法就不适用了，这时才需要上 Euler 求解器。理解了这一点，再看源码就不会困惑为什么它没有激波捕捉模块。
3. 用 MATLAB 复现 MOC 核心推进：Trysome.m 的最小骨架
3.1 初值线生成：从音速线到第一条 C- 特征线
特征线推进需要一个起点。喷管喉部附近流动从亚音速过渡到超音速，MOC 不能直接穿过音速线，因为 M=1 时 μ=90°，特征线与流动方向垂直，推进格式退化。常见做法是在喉部下游取一条靠近音速的初始数据线，线上每点的马赫数、流动角、位置都是已知的。
一个稳定的做法是把初值线选成一条 C- 特征线，这样 C- 族特征线在初始段保持平行，网格不会一开始就严重扭曲。代码如下：
function [x0, y0, theta0, nu0, M0] = initLine(gamma, thWall, MachAxis, nPts) % 生成初值线：取为一条C-特征线，从轴线到壁面 % gamma : 比热比 % thWall : 壁面处流动角 (rad) % MachAxis : 轴线上初始马赫数，通常1.01~1.05 % nPts : 初值线上离散点数 nuAxis = prandtlMeyer(gamma, MachAxis); % 轴线处PM角 Qline = -nuAxis; % C-特征线上的黎曼不变量 theta-nu = const nuWall = thWall - Qline; % 壁面处nu值 MWall = inversePM(gamma, nuWall); % 对应马赫数，后续要检查是否合理 s = linspace(0, 1, nPts)'; theta0 = s * thWall; % 流动角沿初值线线性分布 nu0 = theta0 - Qline; % 由C-不变量直接得到nu M0 = arrayfun(@(n) inversePM(gamma, n), nu0); % 逐点反解马赫数 % 初值线几何：取垂直x轴的直线，y从0到壁面高度，壁面高度由质量流量确定 yWall = 1.0; % 归一化喉道半高 y0 = s * yWall; x0 = zeros(nPts, 1); % 初值线放在x=0平面上 end这段代码的逻辑是：先在轴线上给定一个略大于 1 的初始马赫数，计算出该点的 ν 值，从而确定 C- 特征线上的不变量 Q = θ - ν。整条初值线共享同一个 Q 值，线上每个点的流动角一旦确定，ν 和马赫数就由不变量关系直接推出，不需要额外假设。这是一种自洽的初值构造方式，实际源码里通常会在这个基础上再加一个喉部几何的质量守恒修正，但上面的版本已经能跑通。
3.2 内部点推进：两条特征线求交点
流场内部的一个新点，是由上游一条 C+ 特征线和一条 C- 特征线相交确定的。设上游有两个点 p1、p2，p1 的 C+ 特征线向下游延伸，p2 的 C- 特征线向下游延伸，两条线的交点就是新点 p3。
先算黎曼不变量，再解几何交点，最后反解马赫数。核心代码如下：
function [x3, y3, th3, nu3, M3] = interiorPoint(gamma, p1, p2) % p1: 上游C+点，结构体字段 [x, y, th, nu, M] % p2: 上游C-点，结构体字段 [x, y, th, nu, M] % 返回新点p3的物理量 mu1 = asin(1/p1.M); % p1处马赫角 mu2 = asin(1/p2.M); m1 = tan(p1.th + mu1); % C+特征线斜率 dy/dx m2 = tan(p2.th - mu2); % C-特征线斜率 dy/dx % 两条直线求交点：y-y1=m1(x-x1), y-y2=m2(x-x2) x3 = (p2.y - p1.y + m1*p1.x - m2*p2.x) / (m1 - m2); y3 = p1.y + m1 * (x3 - p1.x); % 黎曼不变量：沿C+传递 P=th+nu，沿C-传递 Q=th-nu P = p1.th + p1.nu; Q = p2.th - p2.nu; th3 = 0.5 * (P + Q); % 新点流动角 nu3 = 0.5 * (P - Q); % 新点PM角 M3 = inversePM(gamma, nu3); % 反解马赫数 end几何上这里做了线性化：把特征线段近似成直线，斜率用上游点的值。步长越小，线性化误差越小。实际源码中如果要提高精度，可以用平均斜率迭代校正一次，但对大多数喷管设计问题，直接线性化已经足够。代码里的 m1、m2 是 dy/dx 形式，注意不要写反。
3.3 壁面点和轴线点：一边是几何约束，一边是对称约束
流场推进到壁面时，新点落在壁面上。壁面几何已知，壁面角 θW 是给定值，这比内部点多了一个约束条件。上游 C- 特征线上的不变量 Q 在壁面点仍然成立，所以新点的 ν 可以直接算出：νW = θW - Q。
壁面是直线段时 θW 不变，计算很简单；壁面是圆弧或任意曲线时，需要把壁面离散成多段折线，每个推进步用当前那一段的斜率作为 θW。实现如下：
function [x3, y3, th3, nu3, M3] = wallPoint(gamma, p2, xWall, yWall, thWall) % p2 : 上游C-特征线上的点 % xWall,yWall,thWall : 当前壁面折线段的起点和斜率 mu2 = asin(1/p2.M); m2 = tan(p2.th - mu2); % C-特征线斜率 % 直线交点：C-特征线与壁面直线 % 壁面直线方程: y-yWall = tan(thWall)*(x-xWall) mW = tan(thWall); x3 = (yWall - p2.y + m2*p2.x - mW*xWall) / (m2 - mW); y3 = p2.y + m2 * (x3 - p2.x); % 流场状态：theta 由壁面给定，nu 由C-不变量确定 Q = p2.th - p2.nu; th3 = thWall; nu3 = thWall - Q; M3 = inversePM(gamma, nu3); end轴线点的处理更简单：由于流动对称，轴线上 θ = 0。这本质上是一个“流动角被给定的壁面点”，把 wallPoint 里的 θW 换成 0，几何约束换成 y = 0 即可。这就是为什么大部分二维喷管 MOC 源码只算上半平面，下半面镜像即可。
3.4 主推进循环的推进秩序和参数调节
有了内部点、壁面点、轴线点三个函数，主程序就能按“行”推进。第一行是初值线，从某一行出发，相邻两点各向前发一条特征线，生成下一行的内部点；行首和行尾分别用壁面点和轴线点补齐。这样依次推进直到到达喷管出口。
% 主循环：逐行推进 % field(1, :) = 初值线 % 对每一行 i，用 field(i, k) 和 field(i, k+1) 生成 field(i+1, k) % 行首补壁面点，行尾补轴线点（或反过来，取决于几何朝向） for i = 1:nRow-1 for k = 1:nCol-1 field(i+1, k) = interiorPoint(gamma, field(i, k), field(i, k+1)); end % 用壁面折线段求上游壁面点 field(i+1, end+1) = wallPoint(gamma, field(i+1, end-1), xW(i), yW(i), thW(i)); % 用轴线条件求轴线点 field(i+1, 1) = axisPoint(field(i+1, 2)); end这个主循环里的关键参数有三个：初值线上的点数 nPts、推进行数 nRow、壁面折线段长度。nPts 太少，插值误差大，等值线出现锯齿；太多则计算量线性增长，但对 MATLAB 来说几千个点毫无压力。推进行数决定出口分辨率，一般取 nPts 的 2 到 3 倍即可。壁面折线段长度要与步长相匹配，每段对应的转角不宜超过 1°~2°，否则壁面点计算的 θW 与实际曲率偏差过大。
4. 跑通与验证：这套特征线源码怎么确认没算错
4.1 最小算例参数组：先跑通再谈精度
验证 MOC 源码最适合的算例是二维对称喷管：喉道半高 1.0，设计出口马赫数约 2.0，壁面扩张半角 10°，比热比 γ = 1.4。这个算例的收敛性好，壁面没有强压缩，出口马赫数通过面积比公式可以独立验证。
| 参数 | 取值 | 说明 |
|---|---|---|
| gamma | 1.4 | 空气 |
| 喉道半高 | 1.0 | 归一化长度 |
| 设计马赫数 | 2.0 | 面积比决定出口高度 |
| 壁面扩张半角 | 10° | 折线段角度 |
| 初值线马赫数 | 1.05 | 轴线上初始值 |
| 初值线点数 | 20 | 越多越光滑 |
| 推进行数 | 60 | 与点数配合 |
| 出口面积比 | 1.38 | 按等熵关系给定 |
需要强调一点：设计马赫数 2.0 对应的面积比大约 1.38，这个值决定了出口半高。代码里不直接给面积比，而是把面积比换算成出口高度作为壁面终点。换算公式就是等熵面积比公式，在验证阶段会用到。
4.2 守恒性验证：用质量流量而不是“看起来像”
特征线法没有显式守恒格式，所以算完不等于算对。最有效的验证指标是单位宽度质量流量沿流向守恒。对二维平面流动，任一横截面上的质量流量应当等于入口质量流量，误差小于 0.5% 说明推进过程没有不可接受的插值损耗。
% 沿某一列（x固定）积分质量流量 % 等熵关系: rho*u 与 M 的关系 % rho*u = rho0 * a0 * M / (1 + 0.5*(gamma-1)*M^2)^((gamma+1)/(2*(gamma-1))) % 对截面上每个点算 rho*u，再对 y 积分 mdotRef = 1.0; % 由初值线积分得到入口质量流量 for i = 1:nCol ycol = [field(:, i).y]; Mcol = [field(:, i).M]; thcol = [field(:, i).th]; rhoU = Mcol ./ (1 + 0.5*(gamma-1)*Mcol.^2).^((gamma+1)/(2*(gamma-1))); ux = cos(thcol); % 轴向速度分量 mdot(i) = trapz(ycol, rhoU .* ux); end residual = abs(mdot - mdotRef) / mdotRef;这段代码的核心是等熵关系：在绝热无粘流动中，ρu 与马赫数之间存在解析关系，不需要单独求密度。trapz 做梯形积分，注意 y 方向的网格不等距——特征线网格天然不等距，trapz 能正确处理。如果残差超过 1%，优先怀疑初值线上的马赫数分布不满足 C- 不变量约束，其次检查壁面折线离散过粗。
4.3 用等熵面积比公式核对壁面马赫数
一个更直接的验证是沿壁面取若干点，用局部流管面积比反推马赫数，与 MOC 算出的壁面马赫数对比。面积比公式是：
A/A* = (1/M) · [(2/(γ+1)) · (1 + (γ-1)/2 · M²)]^((γ+1)/(2(γ-1)))
这里的 A 是当地流管截面积，A* 是喉道面积。对轴对称喷管，A/A* = (y/y* )²；对二维喷管，A/A* = y/y*。写个反函数做对比：
function M = areaRatioToMach(gamma, AR) % 给定面积比，反解马赫数，用于与MOC结果对照 M = 1.1; % 超音速分支初值 for it = 1:100 f = (1/M) * (2/(gamma+1) * (1 + 0.5*(gamma-1)*M^2))^((gamma+1)/(2*(gamma-1))) - AR; % 数值导数 fp = (f(M*1.001) - f) / (0.001*M); M = M - f/fp; if abs(f) < 1e-10, break; end end end注意反解要从超音速分支初值开始。亚音速分支的初值会收敛到 M < 1 的解，那并不是喷管扩张段想要的。把 MOC 壁面点的马赫数按当地 y 坐标换算成面积比，再用上面的函数反算，两者偏差应在 1% 以内。这个验证比整体质量守恒更挑剔，能直接指出问题所在的行号——如果偏差只在某一行之后出现，说明问题在那一段特征线推进的几何处理上。
4.4 云图与等值线：特征线网格的插值坑
特征线网格是不规则四边形网格，不能用 contourf 直接画。需要先插值到规则网格：
% 不规则网格插值到规则网格 F = scatteredInterpolant(X(:), Y(:), M(:), 'linear', 'nearest'); xq = linspace(min(X(:)), max(X(:)), 100); yq = linspace(0, max(Y(:)), 50); [Xq, Yq] = meshgrid(xq, yq); Mq = F(Xq, Yq); contourf(Xq, Yq, Mq, 20);插值方法选 linear，外推用 nearest。喷管壁面附近网格点少，linear 插值会穿过壁面产生不真实的凹陷，这时把数值设为 NaN 再画图能避免误导。可视化只是辅助，真正确认代码正确还得靠 4.2 和 4.3 的定量对比。
5. 特征线源码最容易翻车的三个位置
5.1 初值线的马赫数和流动角必须自洽
很多人把初值线随便取成一条竖线，给每个点一个相同的马赫数，结果推进两三行就出现负 ν 值或马赫数振荡。原因是不变量关系被破坏了：初值线上每个点都是独立的，相邻点的 C- 特征线在物理上应该对应上游同一条马赫波，但人为给定的分布并不满足。解决方法是按 3.1 节的方式构造初值线：先确定一条 C- 特征线的 Q 值，再沿线分配 θ，最后反推 ν 和 M。这个约束保证了初始网格的自洽性。
5.2 壁面曲率与步长的匹配
喷管壁面在喉部下游转弯最急，特征线网格在这个区域也最密。如果壁面折线段每段转角超过 2°，C- 特征线与壁面的交点在相邻两行之间会大幅跳变，导致马赫数等值线出现“折痕”。我一般会要求壁面点间距不超过当地特征线间距的一半。另一种做法是在壁面角突变的点做一次扇形膨胀波处理，把连续转弯离散成若干微小折转角，每段对应一条马赫波。这样虽然计算量增加，但能显著改善出口马赫数均匀性。
5.3 比热比 γ 改动后的连锁反应
γ 不只出现在 ν 公式里，还影响面积比、温度密度关系和出口条件判定。很多人在源码里只改了 γ 的数值，却忘了逆函数 inversePM 里的迭代初值也需要调整：γ 变大会让 ν 的最大值变小，同样的迭代初值可能落在非物理区。验证方法是把初值线的马赫数改成 1.2，重跑一遍，看两条初值线算出的同一出口截面质量流量是否一致。这个“初值无关性”测试比任何画图都更能证明源码的推进是可靠的。
最后一个实用技巧：把 MOC 算出的出口截面马赫数、流动角分布导成数据文件，直接作为 Euler 求解器的入口边界条件。这样特征线源码就从“画图工具”升级成了 CFD 前置设计模块，这也是工程上把快速气动设计和高精度仿真串起来的常规做法。
本文还有配套的精品资源，点击获取
