简介这是一份基于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-其斜率为Cdy/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 减小ν普朗特-迈耶角M1 时为 0γ比热比空气取 1.4注意 ν 的反函数没有解析式需要做数值反解。后面源码里的 inversePM 函数就是干这个的而且它做了两件事一是用牛顿法迭代二是把 dν/dM 的解析导数写进去保证收敛稳定。2.3 为什么这套 MATLAB CFD 源码用 MOC 而不是有限体积推进喷管设计工况的流场本质是一簇膨胀波不包含激波。膨胀波本身就是马赫波的包络用特征线法天然贴合物理过程壁面上的边界条件可以直接嵌入特征线交点计算中不需要像有限体积法那样构造通量函数和人工粘性。MOC 的网格也不是传统意义上的结构化网格而是由特征线交织成的曲线网格网格点就是特征线交点每个点的值由上游两个点直接算出没有隐式耦合。更关键的是精度可控。特征线法的误差来源主要是几何线性化和数值插值而不是格式耗散因此可以精确捕捉膨胀波扇和壁面拐角影响区。对工程预研来说MOC 给出的是解的特征结构而不是被数值耗散抹平的平均场。当然边界条件是有限制的只能处理定常无粘流设计工况下无激波。如果喷管过膨胀或欠膨胀出现内激波特征线法就不适用了这时才需要上 Euler 求解器。理解了这一点再看源码就不会困惑为什么它没有激波捕捉模块。3. 用 MATLAB 复现 MOC 核心推进Trysome.m 的最小骨架3.1 初值线生成从音速线到第一条 C- 特征线特征线推进需要一个起点。喷管喉部附近流动从亚音速过渡到超音速MOC 不能直接穿过音速线因为 M1 时 μ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); % 初值线放在x0平面上 end这段代码的逻辑是先在轴线上给定一个略大于 1 的初始马赫数计算出该点的 ν 值从而确定 C- 特征线上的不变量 Q θ - ν。整条初值线共享同一个 Q 值线上每个点的流动角一旦确定ν 和马赫数就由不变量关系直接推出不需要额外假设。这是一种自洽的初值构造方式实际源码里通常会在这个基础上再加一个喉部几何的质量守恒修正但上面的版本已经能跑通。3.2 内部点推进两条特征线求交点流场内部的一个新点是由上游一条 C 特征线和一条 C- 特征线相交确定的。设上游有两个点 p1、p2p1 的 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-y1m1(x-x1), y-y2m2(x-x2) x3 (p2.y - p1.y m1*p1.x - m2*p2.x) / (m1 - m2); y3 p1.y m1 * (x3 - p1.x); % 黎曼不变量沿C传递 Pthnu沿C-传递 Qth-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, k1) 生成 field(i1, k) % 行首补壁面点行尾补轴线点或反过来取决于几何朝向 for i 1:nRow-1 for k 1:nCol-1 field(i1, k) interiorPoint(gamma, field(i, k), field(i, k1)); end % 用壁面折线段求上游壁面点 field(i1, end1) wallPoint(gamma, field(i1, end-1), xW(i), yW(i), thW(i)); % 用轴线条件求轴线点 field(i1, 1) axisPoint(field(i1, 2)); end这个主循环里的关键参数有三个初值线上的点数 nPts、推进行数 nRow、壁面折线段长度。nPts 太少插值误差大等值线出现锯齿太多则计算量线性增长但对 MATLAB 来说几千个点毫无压力。推进行数决定出口分辨率一般取 nPts 的 2 到 3 倍即可。壁面折线段长度要与步长相匹配每段对应的转角不宜超过 1°~2°否则壁面点计算的 θW 与实际曲率偏差过大。4. 跑通与验证这套特征线源码怎么确认没算错4.1 最小算例参数组先跑通再谈精度验证 MOC 源码最适合的算例是二维对称喷管喉道半高 1.0设计出口马赫数约 2.0壁面扩张半角 10°比热比 γ 1.4。这个算例的收敛性好壁面没有强压缩出口马赫数通过面积比公式可以独立验证。参数取值说明gamma1.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)^((gamma1)/(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).^((gamma1)/(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/(gamma1) * (1 0.5*(gamma-1)*M^2))^((gamma1)/(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 前置设计模块这也是工程上把快速气动设计和高精度仿真串起来的常规做法。本文还有配套的精品资源点击获取