简介面向空气动力学研究者与飞行器设计人员这份代码包以离散涡法为核心解决翼型表面压强系数分布数值计算问题。离散涡法将流动中的涡结构离散为涡元素通过迭代求解涡元素间的相互作用获得流场速度与压力分布进而计算翼型压强系数特别适合处理自由边界与复杂几何外形的气动分析。压缩包共4个文件包括MATLAB主程序脚本、NACA翼型坐标文本以及两个Excel数据表体量仅22KB轻量便携主程序支持自定义翼型参数、调节涡元素数量与分布并能改变攻角研究不同飞行状态下的升力与阻力特性Excel与文本文件则提供NACA4409、NACA4412等翼型的几何和计算结果数据便于对照验证。已有277人学习下载使用者具备基础流体力学与MATLAB编程能力即可直接运行获得直观数值结果适合教学演示、科研验证及初步翼型设计优化能有效降低对商业CFD软件的依赖。1. 离散涡法求翼型压强系数分布为什么你值得把这段代码跑通离散涡法求翼型压强系数分布并不是一个只能在论文里看到的算法名字。简单说它把翼型表面替换成一圈离散的“涡”通过求解每个涡的强度再根据当地速度算出压力系数。我第一次需要评估一款低速翼型的压力分布时手里没有商业软件纠结了两天要不要硬着头皮谢一大套面元法最后发现离散涡法用不到一百行核心代码就能给出可信的结果。这篇文章会把背后的势流假设、离散方式、矩阵求解、压力恢复和参数调试讲透适合气动课程设计、风机叶片初步设计和需要用自研代码交叉验证的从业者。2. 离散涡法求翼型表面速度势流模型与离散方式2.1 离散涡法为什么能算压强系数伯努利关系的使用前提在低速无粘不可压流场中离散涡法把流动看成是均匀来流和一组离散点涡的线性叠加。点涡本身满足拉普拉斯方程因此叠加后的速度场依然是无旋的可以用伯努利方程把速度与压力联系起来。只要边界条件给得正确涡强就能解出来然后表面速度已知压力系数就只剩下一个简单公式$$ C_p 1 - \left(\frac{V}{V_\infty}\right)^2 $$这里 $V$ 是物面当地切向速度$V_\infty$ 是远场来流速度。这个公式对不可压势流严格成立也因此成了离散涡法求压力分布的关键桥梁。需要注意这个关系里的 $V$ 必须是“物面外侧”的当地速度不是某个点涡单独诱导出来的速度。如果只是把离散涡当成一组孤立点涡堆在流场里而不考虑物面边界条件算出来的速度场既不满足不可穿透也不满足物面压力物理C_p 曲线自然不可信。2.2 翼型几何离散用余弦分布重采样物面离散涡法的第一步是确定涡元和控制点。常见做法是沿翼型表面周向布置 N 个涡元每个涡元一个待定强度。网格越密前缘曲面变化越能捕捉但矩阵规模也越大。为了用较少点达到较高精度我习惯用余弦角度分布重新采样翼型表面让前缘和后缘的节点加密中段相对稀疏。import numpy as np def resample_airfoil(XRaw, YRaw, n_panels): # 输入必须是闭合的翼型多边形首尾点重合 # 先按弧长参数化原始点 s np.cumsum(np.hypot(np.diff(XRaw), np.diff(YRaw))) s np.concatenate([[0.0], s]) s_total s[-1] # 用余弦角度分布生成目标弧长位置前缘后缘更密集 beta np.linspace(0, 2 * np.pi, n_panels 1) s_target s_total * 0.5 * (1.0 - np.cos(beta)) # 线性插值得到新的翼面点 X np.interp(s_target, s, XRaw) Y np.interp(s_target, s, YRaw) # 强制首尾闭合于后缘点 X[-1] X[0] Y[-1] Y[0] return X, Y这段代码里最关键的是s_target的构造。beta从 0 均匀变到 2π取余弦后近 0 和近 2π 处的弧长步长变小对应翼型前缘和后缘附近加密。np.interp相当于在原始翼型坐标上做线性插值所以原始坐标需要足够平滑否则插值后物面会有小的折角后续影响系数也会跟着抖。闭合处理是为了让涡元首尾相接这也是避免后缘泄漏的一个前提。2.3 影响系数矩阵单位点涡在控制点上的作用把翼面离散成 N 段后每个涡元可以近似看成一个放置在涡元中心的点涡。对第 j 个涡心它在第 i 个控制点处产生的诱导速度由二维点涡公式给出$$ u_{ij} -\frac{\Gamma_j}{2\pi}\frac{y_i-y_j}{r_{ij}^2}, \qquad v_{ij} \frac{\Gamma_j}{2\pi}\frac{x_i-x_j}{r_{ij}^2} $$这里的 $r_{ij}$ 是控制点到涡心的距离。由于控制点不能和涡心重合我在教学实现里把控制点沿外法向偏移一个极小量工程代码里则用涡段积分代替点涡避免人为偏移带来的误差。影响系数的核心函数可以写成def point_vortex_induced_velocity(px, py, vx, vy): # 单位强度点涡在 (px, py) 处诱导的速度分量 rx px - vx ry py - vy r2 rx * rx ry * ry r2 max(r2, 1e-12) # 防止除零 return -ry / (2 * np.pi * r2), rx / (2 * np.pi * r2)这个函数的输入是控制点坐标和涡心坐标输出单位强度点涡诱导的两个速度分量。注意分母里的1e-12只是兜底真正控制点已经和涡心错开所以正常路径上不会触发。把每个涡元对每个控制点的影响都算一遍按法向投影填入矩阵就得到了影响系数矩阵。这里省略了投影部分因为投影要和法向量一起做后面第 3 章会给完整组装方式。3. 解涡强与恢复压力从矩阵到C_p曲线的代码链路3.1 组装边界条件矩阵不可穿透加后缘条件边界条件很好理解来流速度加上所有离散涡诱导的速度在翼型表面控制点上的法向分量必须为零否则流体就会穿过物面。用公式表示就是$$ \sum_{j1}^{N} A_{ij}\Gamma_j -V_\infty \cdot n_i $$其中 $A_{ij}$ 是第 j 个涡元对第 i 个控制点法向诱导速度的系数$n_i$ 是控制点处外法向单位向量。方程组里 N 个未知涡强N 个物面边界方程还差一个条件。这个条件由后缘条件补上对于尖后缘翼型后缘上下两侧的速度必须连续。在离散涡法里我使用的弱后缘条件是令后缘处相邻两个涡元的涡强互为相反数。组装代码我一般写成这样def assemble_system(X, Y, alpha_deg): # X, Y 是闭合翼型坐标顺序为逆时针 n_panel len(X) - 1 dx np.diff(X) dy np.diff(Y) ds np.hypot(dx, dy) # 外法向逆时针多边形的右侧法向量 (dy,-dx) nx dy / ds ny -dx / ds # 涡心取在每个涡元中点 Xj (X[:-1] X[1:]) / 2 Yj (Y[:-1] Y[1:]) / 2 # 控制点取在涡心沿外法向偏移一个极小量处 eps 1e-6 Xc Xj eps * nx Yc Yj eps * ny A np.zeros((n_panel 1, n_panel 1)) rhs np.zeros(n_panel 1) for i in range(n_panel): for j in range(n_panel): u, v point_vortex_induced_velocity(Xc[i], Yc[i], Xj[j], Yj[j]) A[i, j] u * nx[i] v * ny[i] # 来流速度法向分量的负值作为右端项 alpha np.deg2rad(alpha_deg) rhs[:n_panel] - (np.cos(alpha) * nx np.sin(alpha) * ny) # 后缘条件第一个涡元与最后一个涡元涡强之和为 0 A[n_panel, 0] 1.0 A[n_panel, n_panel - 1] 1.0 return A, rhs, X, Y, nx, ny这里我故意把控制点外移eps是为了在点涡有限元实现里避开奇异点。eps取得越小越接近真实物面但矩阵条件数会变差。A[i, j]的含义是单位强度第 j 个涡元在第 i 个控制点产生的法向速度。右端项则代表来流在物面法向的负贡献。后缘条件行放到矩阵最后一行未知数仍然是 N 个涡强所以矩阵是(N1) × N但我这里初始化了(N1) × (N1)多了一个未知数。实际实现中后缘条件应该用于消去一个未知数或者用拉格朗日乘子。上面的代码为了让结构清晰把右端最后一项设成 0但矩阵多了一列最终求解需要调整。更标准的做法是去掉一个未知数或者把方程数增加到 N1、未知数也 N1但多出的未知数可以解释为一个参考环量。这个细节会在 4.2 节说明。3.2 求解的数值注意点后缘条件与法向量符号np.linalg.solve直接解方阵前要确认法向量方向正确。对于逆时针排列的翼型多边形(dx, dy)是沿翼面边界前进方向右侧法向量恰好是(dy, -dx)也就是代码里的nxdy/ds, ny-dx/ds。这个方向指向流场外侧。如果翼型坐标是顺时针那么这个法向量朝里方程组的解会整体反号压力分布上下表面互换。我在第一次写这个代码时就在这个符号上翻过车最后用零攻角对称翼型做测试才抓出来。后缘条件的实现也有两种常见写法。一种是强后缘条件要求后缘上下两个控制点的速度模值相等另一种是弱后缘条件只约束涡强关系。弱后缘条件简单但对带厚度后缘的翼型会有偏差。工程中我更倾向用强后缘条件把后缘速度差写成一行方程加入系统。这样虽然多花一点装配时间但后缘附近的 C_p 曲线更平滑。3.3 从涡强恢复速度并输出压强系数求解得到涡强数组gamma之后控制点处的速度可以重新计算一遍来流速度加上所有涡元的诱导速度。与组装边界条件不同的是这里不再投影到法向而是保留速度矢量再投影到该点切向得到表面切向速度。因为 C_p 只和速度大小有关取绝对值不影响结果。def compute_cp(X, Y, gamma, alpha_deg, V_inf1.0): # 沿用 3.1 的离散方式 n_panel len(X) - 1 dx np.diff(X); dy np.diff(Y) ds np.hypot(dx, dy) tx dx / ds; ty dy / ds # 逆时针切向 Xj (X[:-1] X[1:]) / 2 Yj (Y[:-1] Y[1:]) / 2 nx dy / ds; ny -dx / ds eps 1e-6 Xc Xj eps * nx Yc Yj eps * ny alpha np.deg2rad(alpha_deg) vx V_inf * np.cos(alpha) * np.ones(n_panel) vy V_inf * np.sin(alpha) * np.ones(n_panel) for i in range(n_panel): for j in range(n_panel): u, v point_vortex_induced_velocity(Xc[i], Yc[i], Xj[j], Yj[j]) vx[i] gamma[j] * u vy[i] gamma[j] * v vt np.abs(vx * tx vy * ty) cp 1.0 - (vt / V_inf) ** 2 return cp这段代码的逻辑和组装矩阵完全对称只是把法向投影换成了切向投影。gamma[j]是第 j 个涡元的环量point_vortex_induced_velocity返回的是单位强度的速度所以乘上gamma[j]再累加。表面切向速度的绝对值用于 C_p是因为我们只关心速度大小不再关心方向。输出cp后按控制点在翼型上的位置就可以画常见的上下表面 C_p 曲线。4. 离散涡法三个必调参数点数、后缘条件与涡核半径4.1 离散点数N先做收敛性扫描再干活离散点数 N 是第一个要确定的参数。N 太小前缘吸力峰会被抹掉后缘压力会出现没有物理依据的波浪N 太大矩阵规模变大而且接近奇异的概率增加。我实践下来中等翼型、攻角不超过 12 度时N 取 100 到 160 已经足够。真正重要的是每次换新翼型时要做一次收敛性扫描。做法很简单固定攻角分别计算 N40、80、160、320 下的前缘最小 C_p 和后缘附近 C_p 值。如果前缘 C_p 峰值在 N 翻倍后变化小于 0.01就认为 N 合格。如果峰值还在明显变化就继续加 N。这个扫描对初学阶段尤其重要因为很多看起来“不光滑”的 C_p 曲线其实就是 N 不够而不是物理上就该振荡。N前缘峰值偏差矩阵规模单次求解耗时40约 5%40×40毫秒级80约 1%80×80毫秒级160小于 0.2%160×160毫秒级320基本稳定320×320百毫秒级上面是我的经验值具体翼型会有差异。需要注意前缘越薄的翼型对 N 越敏感。如果是尖前缘甚至要把前缘附近的局部网格单独加密否则全局加密效率很低。4.2 后缘条件的选择强约束、弱约束与矩阵条件数后缘条件直接影响总环量也就直接影响升力和前后缘压力差。弱后缘条件只令后缘两个涡元强度之和为零代码简单但存在两个问题一是当后缘有厚度时后缘两侧速度连续并不等价于涡强反号二是它可能和物面边界条件行产生线性相关导致矩阵奇异。强后缘条件则要求后缘上下两个控制点的切向速度相等。我通常构造一行系数后缘上控制点速度减去下控制点速度令其等于零。这个条件多一个方程同时会给系统带来一个额外的自由度需要用一个拉格朗日乘子或者去掉一个未知数来平衡。直接把它加在矩阵最后一行效果往往也足够但要注意该行系数数量级与边界条件行相差很多最好乘以一个特征长度比如弦长让矩阵各行量纲一致。否则np.linalg.solve可能给出很大的伪解。4.3 涡核半径数值奇异的灭火器如果实现中用点涡直接表示涡元涡核半径是一个不能回避的参数。涡核半径 r_core 修正点涡诱导速度公式使速度在距离涡心很近时不再趋于无穷$$ v_\theta \frac{\Gamma}{2\pi}\frac{r}{r^2 r_{core}^2} $$当 r 远大于 r_core 时公式退化为普通点涡当 r 接近零时速度趋向有限值。r_core 一般取弦长的 0.001 到 0.01。太小收敛困难太大压力峰值被压平。我的习惯是先取 0.005再做一次峰值敏感性测试如果峰值随 r_core 变化剧烈就要减小它并观察是否矩阵求解变得不稳定。5. 离散涡法求压强系数常见的5个翻车现场与排查方法5.1 现象矩阵奇异或线性求解报错这是最常遇到的问题现象是np.linalg.solve抛出奇异矩阵错误或者返回一堆量级极大的数值。原因通常有两个一是后缘条件没加导致方程数少于未知数二是翼型坐标不闭合首尾点没有严格重合涡元出现缺口。解决办法是打印矩阵条件数检查X[0]X[-1]和Y[0]Y[-1]是否严格成立。浮点误差会导致判断不通过建议手动把首尾赋值成同一个点。若加了后缘条件仍奇异检查后缘条件行是否与某一行几乎线性相关必要时乘以弦长来调节量纲。5.2 现象零攻角对称翼型算出非零升力对称翼型零攻角时无论厚度分布如何升力都应为零上下表面 C_p 应完全对称。如果算出来有明显的升力首先怀疑法向量方向。法向量取反后来流在物面法向投影也反号涡强解会镜像翻转最终上下表面压力互换积分后的升力不为零。排查方法很简单打印第一个涡元中点的外法向量手动判断它是否指向翼型外侧。另外一个容易忽略的是翼型坐标绕向如果绕向是顺时针外法向量公式的符号也要相应调整。5.3 现象前缘吸力峰值偏低前缘吸力峰是翼型压力分布最重要的特征之一。如果 N 已经取得足够大但峰值仍然比理论值低问题多半出在网格分布上。均匀按弦向布点会在前缘留下过宽的间隙控制点无法捕捉前缘附近陡峭的压力梯度。解决方法是改用余弦角度分布并检查前缘附近的第一个控制点距离前缘是否在 0.1% 弦长以内。还有一个容易被忽略的点如果使用的是点涡近似且控制点外移偏移量eps过大前缘峰值会被系统性压低把eps从 1e-4 降到 1e-6 再比较一次。5.4 现象后缘附近的C_p曲线出现扭曲后缘扭曲通常表现为压力系数在靠近后缘时突然上冲或下跌曲线不光滑。原因有两类一类是后缘点没有闭合上下翼面最后一小段涡元之间存在缺口另一类是在后缘附近使用了过短的涡元导致相邻涡元夹角极小点涡影响系数计算时角度跨过 π 边界。解决方法是强制后缘点重合并把后缘附近网格适当放松避免涡元长度比平均长度小一个数量级。如果用的是涡段积分公式还要在角度差计算里加上跨 π 修正。5.5 现象攻角稍微增大C_p就振荡离散涡法本质上是无粘势流方法没有模拟流动分离。攻角超过一定范围后真实流动已经分离势流解自然不稳定。常见表现是前缘吸力峰越来越尖最终出现非物理的局部振荡。解决方法是限制攻角范围一般不超过 15 度。如果必须计算大攻角就需要升级到非定常离散涡法在后缘释放自由涡模拟尾迹让涡量随流场输运。否则再加密网格也解决不了物理模型缺失带来的振荡这时候不是代码 bug而是方法适用范围到了尽头。6. 验证你的压强系数分布对称性、升力定理与收敛习惯6.1 零攻角对称翼型上下表面曲线重合性拿到第一版 C_p 分布后先用一个对称翼型、零攻角做冒烟测试。理想情况下上下表面的 C_p 曲线应当完全重合升力系数为零。如果两条曲线打得厉害优先检查坐标绕向和法向量符号如果只是小范围不重合再检查后缘条件。这一步可以过滤掉一大半低级问题比直接和实验数据对比更快速。6.2 用总升力系数做自检两种算法互验离散涡法解出的涡强本身可以积分出总环量。根据库塔-儒可夫斯基定理二维升力系数与总环量有直接关系$$ C_l \frac{2\Gamma_{total}}{V_\infty c} $$另一方面把 C_p 沿翼型表面积分也能得到法向力系数再换算到升力系数。这两个数值应当高度一致。如果差异超过 1%说明速度恢复或者积分方向有问题。我一般把这条当作最有价值的数值自检因为它把压力场和涡量场联系到了一起能捕捉到矩阵装配中隐藏的错误。6.3 把收敛性检查变成习惯不管多么着急每次更换翼型或攻角我都先跑一组 N80、160、320 的收敛性测试再决定最终数据取哪一个网格。这看起来多花几分钟实际上比事后发现结果不可信再返工省得多。我自己就曾因为想省时间直接用 N60 算了一个后缘压力分布结果局部隆起差点被当成真实气动特征。从那以后收敛性扫描就成了固定流程。希望这套判断方法也能帮到你。本文还有配套的精品资源点击获取
