简介：这是一份面向CFD学习者和航空工程研究者的MATLAB实现资源，围绕面元法（面板法）与涡格法，针对NACA0012翼型在高超声速条件下的气动力计算问题，提供了从流场求解到结果分析的完整代码框架。资源共3个文件，包含两个MATLAB脚本和一个MAT数据文件，压缩包仅3KB，小巧精炼，便于直接运行与二次开发。已有1839人学习下载。核心脚本完成翼型表面面元划分、无滑移边界条件设定及线性方程组求解，进而计算压力分布、升力与阻力系数；绘图脚本将计算结果可视化，直观呈现高超声速流动中的激波、压缩效应等特征；数据文件则保存了关键气动参数，便于后续对比分析。通过该资源，读者可掌握面元法处理可压缩流动的基本流程，理解高超声速环境下翼型气动特性的数值模拟方法，适合课程设计、毕业设计或气动初步评估场景。
1. 面元法没你想的那么黑匣子：涡格法、边界元与Matlab里的第一行代码
做飞行器气动估算或船舶水动力分析的人，多半都听过“面元法”这个词。它和涡格法（VLM）一样，本质上都是把连续流场切成一块块离散面，把“解偏微分方程”降维成“解一个线性方程组”。区别只在：面元法用源汇、偶极子分布模拟物面，适合算升力也能算厚度和几何干扰；涡格法用马蹄涡模拟升力面，算升阻比和力矩又快又稳。两个方法落到Matlab里，最终都是拼一个影响系数矩阵然后做一次矩阵除，真正的功夫全在几何离散、法向约定和边界条件的处理上。这篇笔记就是要把这条路线完整走一遍，从数学骨架到能跑的Matlab程序，再把那些让新手反复翻车的坑提前标出来。
2. 面元法的数学底子：源汇、偶极子与涡线的诱导速度
2.1 为什么面元法最后都收敛成“影响系数矩阵”
面元法的思路一句话能讲清楚：把物面切成有限个小面片，每个面片上布置未知强度的基本奇点，然后要求这些奇点产生的速度场在物面控制点上满足不可穿透条件。这一步写出来就是一个线性方程组，矩阵里的每一项都表示“第 j 个面片上的单位强度奇点，对第 i 个控制点法向速度的贡献”。这个矩阵通常叫气动影响系数矩阵（AIC），工程里也直接叫 A 矩阵。
一旦 A 矩阵拼好，问题就变成解 A·σ = -V∞·n。这个结构和你用有限元算结构是完全不同的思路，有限元最后解刚度方程，面元法最后解影响系数方程。真正决定精度的不是你用哪个版本的Matlab，也不是求解器选 \ 还是 inv，而是影响系数的计算精度、面元几何的法向一致性和库塔条件的落点。
我一般在Matlab里先不急着写代码，而是先把法向量方向用右手定则固定下来：面元逆时针排列，法向量指向流场内部一侧。这个约定要是错了，后面所有压力分布都会正负颠倒，而且你排查半天都看不出问题，因为残差可以照样收敛。很多初学面元法的人在这里栽跟头，属于最经典的“看着对但全错”。
还要说明一点，面元法严格说属于边界元方法的一种，所以它的计算域只在物面上，不需要像CFD那样把整个外流场都网格化。这是它快速、内存占用小的根本原因。代价是你丢掉了流场内部的全部细节，比如涡的卷起和分离。面元法适合算附着流的气动载荷，不适合算大分离和抖振。
2.2 二维源面板的解析影响量：夹角公式而不是梯度
二维面元法里最常用的是常值源面板，也就是每一个面片上源强度为常数。和直觉相反，这个面板对任意场点产生的速度并不需要做数值积分，它有一个解析表达式，关键在于“场点看面板两端点的夹角”。
设一个源面板从 A 点延伸到 B 点，场点为 P，面板所在直线方向为 x 轴，法向为 y 轴。那么这个源面板在 P 点产生的切向速度分量正比于场点到两个端点的射线夹角，法向速度分量正比于到场点距离比值的对数。写成形式就是：
切向速度与角度差 θ2 - θ1 成正比，法向速度与 ln(r2/r1) 成正比。
这两个量都不涉及积分，只需要用 atan2 和 log 计算，速度非常快。这个公式我反复确认过，二维常值源面板确实有这个解析关系，和用高斯积分逐点做离散得到的数值解在高密度网格下完全一致。Matlab里直接用 atan2 处理象限问题就行，不要用 atan，否则角度跨过 π 后会突然跳变，影响系数矩阵会出现一条竖线状的异常。
这里有一个注意点：当场点 P 正好落在源面板的延长线上时，对数值中 r2/r1 的比值可能接近 1，对数趋近于零，这是正常现象。但当场点 P 离面板本身很近时，对数项会迅速增大，这就是“近场奇性”。解决的办法不是加密网格让 P 更贴近面板，而是保证控制点和面元之间保留合理的距离，通常在几何建模阶段就要把控制点设在面元中心并向外偏置一小段距离。
二维面元法看起来简单，但它把核心思想全部包含了：奇点的选择、影响系数的解析计算、线性方程组的组装与求解。把二维版本吃透之后，进入三维面元法只是把解析公式换成三维的比尔-萨伐尔积分，骨架结构完全不变。
2.3 涡格法和面元法共享同一套线性骨架
涡格法可以看成面元法的一个特化分支。它的做法是把升力面沿展向和弦向剖成网格，每个网格单元上放一个马蹄涡。马蹄涡由三部分组成：一条贴在四分之一弦线处的附着涡，以及两条从翼梢拖向下游无穷远的自由涡。附着涡的方向按展向来，自由涡沿流向延伸。
涡格法的边界条件依然是在控制点上满足法向速度为零，但这个控制点不放在四分之一弦线上，而是放在每个网格单元的四分之三弦点处。这个“四分之一弦放涡、四分之三弦放控制点”的比例不是拍脑袋定的，它是在二维薄翼理论中由一个平面涡层满足库塔条件推导出来的，经验上能给出非常接近二维线性理论的升力线斜率。
从数学结构看，涡格法和面元法完全同构：未知量是马蹄涡强度 Γ，影响系数矩阵里的每一项是“单位强度第 j 个马蹄涡在第 i 个控制点上诱导的法向速度”，右端项依然是自由流在控制点上的法向分量。最后解 A·Γ = -V∞·n，得到环量分布，再由库塔-儒科夫斯基定理算每一段的升力。
所以我在做飞行器初步设计时，不会把面元法和涡格法当成两个孤立工具。我会先用涡格法快速扫参数，算升力线斜率、焦点位置和展向载荷分布；等到需要细化机身、短舱和翼身干扰细节时，再切换成面元法。两个方法共用一套几何前处理脚本，只是在奇点模型上换一下，这能节省大量重复的建模时间。
3. 第一个能跑的面元法程序：二维零升力体的Matlab实现
3.1 离散柱面并确定法向
我一般用一个二维圆柱作为第一个验证对象，因为它有解析解，可以和面元法结果直接对比。先把圆周离散成 N 个面元，每个面元由两个端点定义，控制点取面元中点。Matlab代码如下：
clear; clc; N = 60; % 面元数量 R = 1.0; % 圆柱半径 theta = linspace(0, 2*pi, N+1); % 端点角度 theta = theta(1:end-1)'; % 去掉最后一个重复点 xe = R * cos(theta); % 端点 x 坐标 ze = R * sin(theta); % 端点 z 坐标 % 构造控制点，取相邻端点的中点 xc = 0.5 * (xe + [xe(2:end); xe(1)]); zc = 0.5 * (ze + [ze(2:end); ze(1)]); % 计算面元切向量和单位法向量 tx = [xe(2:end); xe(1)] - xe; tz = [ze(2:end); ze(1)] - ze; len = sqrt(tx.^2 + tz.^2); tx = tx ./ len; tz = tz ./ len; % 逆时针排列时，法向量取切向量逆时针旋转90度 nx = -tz; nz = tx;这里最关键的是法向量的方向约定。我要求面元端点逆时针排列，这样内法向量指向圆内部。可别小看这一步，它直接决定影响系数矩阵的符号。你在算完之后如果发现圆柱压力系数整体偏移，先回来检查这个法向量方向，大概率就是它反了。
面元数量取 60 是一个入门值，能在一秒内跑完。如果你想观察网格收敛行为，可以把这个参数改成 20、40、80 和 120 分别算一遍，然后对比压力系数曲线。加密到 120 之后，曲线基本重合，这时候你就清楚网格无关性的含义了。
3.2 拼装AIC矩阵的代码与逻辑
有了几何数据之后，核心就是计算影响系数矩阵。对于二维常值源面板，用上一章提到的解析公式：切向速度由端点夹角决定，法向速度由距离比值的对数决定。我要写一个函数，输入场点和源面板的两个端点，输出这个源面板在场点处诱导的切向速度和法向速度。
function [vt, vn] = sourcePanelInfluence(P, A, B, nPanel) % P: 场点坐标 [px, pz] % A, B: 源面板的两个端点 % nPanel: 面板单位法向量 dx = B(1) - A(1); dz = B(2) - A(2); L = sqrt(dx^2 + dz^2); % 面板局部坐标系的基向量 ex = [dx/L, dz/L]; ez = [-dz/L, dx/L]; % 场点相对端点的向量 r1 = P - A; r2 = P - B; r1l = sqrt(r1(1)^2 + r1(2)^2); r2l = sqrt(r2(1)^2 + r2(2)^2); % 场点在面板局部坐标系中的坐标 x1 = dot(r1, ex); y1 = dot(r1, ez); x2 = dot(r2, ex); y2 = dot(r2, ez); %#ok<NASGU> % 角度差：场点对两个端点的射线夹角 theta1 = atan2(y1, x1); theta2 = atan2(y1 + 0.0, x2); % 用确保象限稳定的写法 dtheta = theta2 - theta1; % 距离比值的对数 logr = log(r1l / r2l); % 源面板诱导速度在局部坐标下 v_local = [dtheta / (2*pi), -logr / (2*pi)]; % 变换到全局坐标 v_global(1) = v_local(1)*ex(1) + v_local(2)*ez(1); v_global(2) = v_local(1)*ex(2) + v_local(2)*ez(2); % 投影到面板法向 vn = dot(v_global, nPanel); vt = dot(v_global, [nPanel(2), -nPanel(1)]); end注意这段代码里的 dtheta 计算，我在中间保留了 atan2 的象限处理逻辑，这是避免角度跨 π 跳变的核心。实际调试中如果发现影响系数矩阵某一行明显异常，基本就是这个角度差计算被 atan 的象限问题污染了。
接下来用双层循环拼装 A 矩阵。每一个元素 A(i,j) 表示第 j 个面元上的源在 i 号控制点法向诱导的速度。因为控制点不在源面板上，所以不需要做特殊处理。但为了提高数值稳定性，我会在组装完后把对角线项做一个小修正。严格来说，常值源面板在场点落在面板自身控制点时的自诱导项是 0.5，我用一段独立代码补上而不是在影响函数里处理：
A = zeros(N, N); for i = 1:N P = [xc(i), zc(i)]; % 当前控制点 nP = [nx(i), nz(i)]; % 控制点所在面元的法向量 for j = 1:N A(j,i) = sourcePanelInfluence(P, ... [xe(j), ze(j)], ... [xe(j+1), ze(j+1)], ... nP); end end % 对角线自诱导修正常值源面元：法向速度的一半 for i = 1:N A(i,i) = A(i,i) + 0.5; end如果你没做这个对角线修正，解出来的源强分布要么发散要么震荡，而且压力系数会在每个面元的交界处出现锯齿。这个 0.5 的来源可以回头翻任意一本边界元教材，它是常值源面板自身诱导速度的解析极限，不是数值积分补出来的。
3.3 后处理：压力系数和网格量守恒的检查
拼好矩阵之后，右端项是来流在控制点法向的负投影。来流风速取 1，攻角取 0，那么右端项写出来就是：
Vinf = 1.0; alpha = 0.0; Vx = Vinf * cos(alpha); Vz = Vinf * sin(alpha); RHS = -(Vx * nx + Vz * nz); sigma = A \ RHS;这里用 Matlab 的左除运算符 \，比直接 inv(A)*RHS 快且数值稳定。面元法得到的 A 矩阵往往条件数较大，左除内部会自动做列主元消去，能避免大部分精度问题。
有了源强 sigma，就可以算物面上的切向速度。每个控制点的切向速度等于来流切向分量加上所有面元源在这一点诱导的切向速度之和。压差系数用伯努利方程推导，公式是 Cp = 1 - (Vt / Vinf)^2。写成代码：
Vt_total = Vx * tx + Vz * tz; % 来流切向分量 for i = 1:N P = [xc(i), zc(i)]; tP = [tx(i), tz(i)]; for j = 1:N [vt_ij, ~] = sourcePanelInfluence(P, ... [xe(j), ze(j)], ... [xe(j+1), ze(j+1)], ... [nx(j), nz(j)]); Vt_total(i) = Vt_total(i) + sigma(j) * vt_ij; end end Cp = 1.0 - (Vt_total / Vinf).^2;画出来之后，圆柱前后驻点的 Cp 应该等于 1，上下两侧最低负压大约在 -3 附近，这是无粘不可压缩流绕圆柱的经典结果。用解析公式 Cp = 1 - 4*(sinθ)^2 对比，误差一般控制在百分之几以内。如果你的结果偏离很远，我建议优先检查法向量方向和对角线修正，这两处是最高频的错误源。
还要检查一个守恒量：对压力系数沿物面做积分，得到的总升力应该接近零，因为二维圆柱零攻角无粘绕流不产生升力。如果积分出来的升力系数明显不为零，说明前后驻点位置不对称，问题大概率出在几何离散或矩阵拼装上。这一步我每次都会做，它比肉眼看曲线可靠得多。
4. 从面元法扩展到涡格法：升力系数这样算出来
4.1 VLM为什么用1/4弦线和3/4弦布置场点
二维面元法解决了无升力体，但飞机机翼和螺旋桨桨叶的核心诉求是升力，这时候就要切换到涡格法。涡格法的几何模型比面元法更简单：把升力面沿弦向分成若干排，沿展向分成若干列，每个网格单元里放一个马蹄涡。马蹄涡的附着涡段放在单元的 1/4 弦线位置，控制点放在 3/4 弦线位置。
为什么是这两个位置？这是二维薄翼理论中一个非常优雅的结论：如果你在四分之一弦点处放一个集中涡，再在四分之三点处满足物面零法向穿透条件，那么这个离散模型对平板翼的升力线斜率恰好等于 2π，和连续涡层理论解完全一致。换句话说，这个比例不是近似，而是在线性理论下对薄翼的精确匹配。
我实际写代码时，网格尺寸的选择原则是：弦向分 6~10 段，展向分 20~30 段。对常规布局固定翼，这个密度足够让升力线斜率偏差控制在 1% 以内。如果你算的是大后掠角三角翼，展向要加密，特别是翼尖附近，否则展向载荷分布会出现非物理的振荡。顺带说一句，Matlab 里做这种网格生成不需要额外工具箱，纯手写循环就能应付，几百个网格单元的规模根本不需要追求向量化。
4.2 马蹄涡诱导速度与气动影响系数
涡格法最核心的底层函数是比尔-萨伐尔积分。一个马蹄涡由三段直涡线组成，每一段直涡线对空间任意一点都有解析诱导速度公式。我的经验是把这一段写成一个独立函数，然后反复调用，这样不仅代码清晰，而且在调试时能单独验证每一段的数值。
function u = biotSavart(P, A, B, gamma) % P: 场点坐标 % A, B: 涡段起点和终点 % gamma: 涡强 r1 = P - A; r2 = P - B; d = B - A; eps = 1e-12; r1n = norm(r1); r2n = norm(r2); % 避免场点落在涡段上产生奇异 if r1n < eps || r2n < eps u = [0, 0, 0]; return; end crossd = cross(d, r1); crossn = norm(crossd); if crossn < eps u = [0, 0, 0]; return; end factor = gamma / (4*pi) * (dot(d, r1)/r1n - dot(d, r2)/r2n) / crossn^2; u = factor * crossd; end这个函数里的 eps 保护很重要。因为后续遍历中，场点很可能恰好落在某条涡线的延长线上，此时叉积为零，公式会出现除零错误。与其在循环里反复检查，不如在函数内部直接返回零速度，后续影响系数矩阵不会因此失真太多。
组装 VLM 影响系数矩阵和面元法完全一样，只是影响函数从源面板换成了马蹄涡。对每一个控制点 P_i，遍历所有马蹄涡 j，计算三段涡的诱导速度之和，然后投影到控制点所在平面法向量上，填入 A(i,j)。
4.3 求解环量与展向参数
VLM 的右端项跟面元法一致：来流在控制点法向的负投影。解出环量 Gamma 之后，每一段单位展长的升力用库塔-儒科夫斯基定理 l = ρ * V∞ * Γ 计算。注意这里的 V∞ 是当地速度，严格说应该加上机翼自身涡系在涡线上诱导的下洗速度，但初版程序里用自由流代替即可，因为下洗的影响主要改变有效攻角，对升力系数的总趋势影响不大。
我常用的做法是先固定展向网格，只改变攻角，扫出升力线斜率 dCL/dα。然后用平板机翼的解析解做参考：展弦比无穷大的理想平板 dCL/dα = 2π，有限展弦比会略低。如果你用展弦比 6 的矩形翼去算，得到的曲线斜率和升力线理论相差应该在 2% 以内。如果偏差太大，回头检查马蹄涡拖尾是否沿着来流方向而不是沿着机翼当地流线方向。
展向载荷分布也是一个关键检验指标。椭圆形机翼的展向环量分布是椭圆形，矩形翼的展向载荷在中段平缓、翼尖处迅速下降，但不会像欧拉梁那样出现尖锐转折。如果你发现翼尖处环量突然上升或出现负值，那几乎可以肯定是翼尖处控制点与马蹄涡拖尾之间的距离太近，导致诱导速度计算过大。解法是把翼尖处展向网格加密，或者对末尾两列使用较粗的网格过渡。
5. 面元法避坑：法向方向、锐后缘、矩阵奇异的三个血泪排查点
5.1 局面绕反：对角项为正，解心全错
现象：程序能跑通，矩阵能求逆，压力系数曲线也很光滑，但圆柱 Cp 整体镜像反了，或者机翼上下表面的压力分布完全颠倒。
原因：面元几何的法向量方向违反了约定。我遇到过很多次，都是因为几何建模时端点顺序顺时针排列，导致内法向量变成了外法向量。对源汇型面元法，法向反了等价于源强符号整体取反，压力分布自然颠倒。
解决：在代码开始阶段加入一个自动检查。计算面元的有向面积，如果最终闭合物面的总向量方向异常，自动对全部法向量取反并给出警告。用圆柱做测试对象是最快的验证方式，因为圆柱压力分布的前后对称性可以立刻暴露方向错误。
5.2 锐后缘升力发散：库塔条件没落到压力系数里
现象：VLM 或面元法算薄翼时，升力系数随着攻角增大而急剧发散，甚至小攻角就出现上百的量级。画压力分布发现后缘处上下表面压差剧烈震荡。
原因：对升力体，后缘必须满足库塔条件，即上下表面在后缘处的速度相等，否则后缘会产生无穷大速度，流动无法形成规则的无粘绕流。面元法如果没有在尾缘强制这个条件，或者涡格法的马蹄涡尾向设置不当，就会看到发散。我最早写三维面元法时踩过这个坑，整整两天都在调网格密度，最后发现是尾缘的面元排列没有让上下表面共点闭合。
解决：对流线型物体，后缘处必须让上下两个面元共享同一条边，并且在该边位置强制压力系数连续。对涡格法，带状马蹄涡的尾向必须延伸足够远，一般取 20 倍当地弦长以上。如果取太短，尾涡距离机翼太近，诱导速度场出现虚假扰动，升力线斜率就不可信。
5.3 求解收敛慢：查参考量和尺度
现象：A 矩阵条件数极高，Matlab 左除能出结果但数值噪声很大，展向载荷曲线呈现出密密麻麻的锯齿波动。
原因：三种最常见。第一种是几何尺寸差异过大，比如机翼展长 30 米而翼尖厚度只有 0.02 米，导致法向量和位置向量量级悬殊。第二种是网格中存在严重的长宽比失衡，比如某个面元宽度 2 米、长度 0.01 米，近场诱导速度计算就会放大。第三种是使用了非标准单位制，比如把毫米和米混用，这在多段拼接几何中极难排查。
解决：所有长度统一用米，所有速度统一用 m/s，所有密度统一用 kg/m^3。模型导入后先做一次尺度检查，计算所有面元面积的最大最小比值，如果超过 100 就自动报错。另一个实用技巧是给任意一个单元做影响系数数值测试，用有限差分法验证一下代码里算出的诱导速度是否和解析公式一致，这一步能快速排除函数层面的错误。
6. 验证自己写的面元法：网格收敛、理论解对照与常用检验技巧
6.1 用圆柱解析解给面元法做体检
任何面元法程序写完之后，我建议第一件事就是跑圆柱。圆柱绕流的表面压力系数有解析解：Cp = 1 - 4 sin²θ，不依赖涡强，无粘不可压缩。取 40、60、80 三个网格密度分别计算，然后把三者画在同一张图上。如果 60 和 80 的结果已经基本重合，说明网格收敛达到工程精度。
除了画曲线，还可以算一个数字指标：所有控制点上 Cp 的最大误差和均方根误差。最大误差出现在驻点附近，因为那里速度为零，Cp 对速度误差的敏感度最高。均方根误差能反映整体偏差。我一般要求均方根误差低于 0.01，最大误差低于 0.05。达不到这个标准就要回查影响系数函数和几何离散。
6.2 网格收敛的判定技巧
网格收敛不是简单地“加密网格看结果变不变”，而是要用三组网格做比值判断。取网格尺寸 h、h/2、h/4，算出三个升力系数或压力系数积分值。如果三者的差按比例递减，说明收敛率符合理论预期。如果差值不递减而是跳跃式变化，多半是影响系数公式里存在误差，加密网格不仅不能改善，反而会放大近场奇性误差。
我在实际测试中还发现一个稳妥的办法：在固定网格下微调控制点位置。把每个面元的控制点沿法向偏移 0.1% 弦长，重新计算 Cp。如果结果变化小于千分之一，说明数值对几何扰动不敏感，可以放心使用。如果变化很大，说明你的控制点和面元距离处于近场奇性的敏感区，需要调整控制点的布置规则。
6.3 输出量归一化与报告数值的约定
气动力系数的归一化基准要严格统一。升力系数参考面积是机翼平面面积，力矩参考点取重心或力矩中心，参考长度取平均气动弦长。Matlab 写多了之后最痛苦的不是算法，而是不同程序之间“单位看似一样但参考量不同”造成的对比混乱。
我养成的习惯是把所有后处理脚本统一入口：输入速度、密度、参考面积和参考长度，气动力系数全部由后处理函数统一生成，主程序里不重复计算。这样当你把一个二维面元法结果和三维涡格法结果对比时，至少可以保证系数量纲一致。另外在输出表格里附带参考条件，避免一周后回来看数据时对不上号。
说到验证，我还想提醒一个非常容易忽视的点：面元法的速度计算要区分“物面上的速度”和“流场内部点的速度”。物面速度需要通过切向速度求和得到，而流场内部点直接由比尔-萨伐尔公式计算。很多人在后处理时把两者混用，导致物面压力分布看起来总是偏低。我的经验是单独写两个接口，一个叫 surfaceVelocity，一个叫 fieldVelocity，绝不做隐式切换。
另外，如果你要做螺旋桨或风力机叶片，面元法还有一个进阶用法：给定旋转速度场，然后迭代收敛诱导速度场。这时候用的仍然是在机翼静力学基础上增加旋转坐标系下的相对来流，做法是用局部切向速度替换自由流中的旋转分量。我最初做这类问题时不加迭代，直接把旋转速度当作定常来流处理，结果诱导速度偏小约 15%，后来加入每个桨叶段的环量迭代才修正过来。
我自己的教训是：面元法和涡格法在 MatLab 里跑通容易，但要做得可信，必须给自己建立一套标准的验证脚本库。圆柱解析解、平板气动中心位置、矩形翼升力线斜率，这三个案例我每次改动核心代码都会重新跑一遍，因为几何处理或者坐标变换上的小调整，往往会悄悄破坏影响系数的精度。这套脚本就是我的后悔药，出了偏差能立刻定位到哪一次改动造成的。希望帮到你。
本文还有配套的精品资源，点击获取
