做螺旋桨气动分析这些年我越来越觉得叶片单元动量理论BEMT是工程阶段最值得优先掌握的估算工具。它不需要网格不需要湍流模型只要有一组螺旋桨几何参数——弦长分布、扭转角分布、翼型升阻极曲线——再加一个转速和一个来流速度就能在Matlab里把所有关键性能指标拉出来。这次的工作是把一套固定几何的螺旋桨放在恒定转速下从接近悬停的小前进比一直扫到高速巡航的大前进比完整观察拉力、扭矩、功率和效率怎么变化。文章会把理论怎么落进代码、代码里哪些参数最容易出错、曲线异常时怎么排查都写清楚。正在入门螺旋桨性能估算的人或者做无人机动力选型想快速出一版性能包的工程师都可以照着复现。1. 项目拆解恒定转速下的前进比扫描到底在写什么1.1 几何输入与性能输出的边界先定下来标题里有两组词需要先拆清楚。一组是“给定螺旋桨几何形状”另一组是“恒定转速下不同前进比”。前者决定计算对象后者决定计算工况。给定螺旋桨几何形状通常是指这样几样东西桨叶半径R、桨叶数B、沿展向的弦长分布c(r)、沿展向的桨距角扭转分布θ(r)以及每个截面翼型的升力系数和阻力系数随攻角变化的极曲线。前四样决定了几何外形最后一样决定了气动响应。很多人把注意力放在前四样上却低估了翼型数据对结果精度的影响。实际上同一把螺旋桨如果翼型极曲线取的是低雷诺数薄翼型取的是带襟翼的厚翼型算出来的推力和效率可以差出两成以上。性能输出的指标通常是总拉力T、总扭矩Q、总功率P和效率η。这四个量在不同前进比下变化趋势完全不同拉力随前进比升高而显著下降扭矩也是这样功率因为和转速绑定下降速度又和扭矩不完全同步效率则会在某个前进比附近出现峰值。理解这些曲线的形状比单看某一个点的数值重要得多。为了方便对比工程上会把T、Q、P无量纲化转成推力系数CT、扭矩系数CQ、功率系数CP配合前进比J一起看。后文的代码就是按这一套无量纲定义输出的。符号含义常用单位或量纲R桨叶半径mB桨叶数无c(r)当地弦长mθ(r)当地桨距角radω或Ω角速度rad/sJ前进比 J V∞/(nD)无CT推力系数 T/(ρn²D⁴)无CP功率系数 P/(ρn³D⁵)无η效率 η J·CT/CP无这里的n是转每秒D是螺旋桨直径。这套无量纲表示在螺旋桨行业里是通用的方便不同尺寸、不同转速的桨放在同一张图里比较。1.2 为什么选“恒定转速前进比扫描”这个研究方式螺旋桨性能研究有几种常见边界。第一种是固定转速变化来流看推力和功率曲线第二种是固定拉力反推转速和功率第三种是固定油门扭矩上限看转速和拉力变化。标题选的是第一种它的好处是直接对应大多数直流电机驱动的无人机系统电调输出转速指令转速环把RPM锁得很死随着飞行速度变化螺旋桨的负荷随之改变。前进比J V∞/(nD)把转速和尺寸的影响归一化了因此不同尺寸的螺旋桨可以在同一张图上比较。J小表示来流相对螺旋桨旋转速度很小对应悬停或低速爬升J大表示来流明显对应巡航或高速飞行。恒定转速下扫描J实际上就是固定RPM让来流速度从接近0慢慢加到高速。这种扫描方式还有一个工程上的好处它对“飞机从静止加速到巡航”这个完整飞行剖面给出了一张连续的推进系统性能图。螺旋桨在不同速度下会不会负荷过大、能不能保持效率、电机在哪个区间最省电都可以从J-CT、J-CP、J-η曲线上直接读出来。所以这篇工作虽然听着像课程设计实际上就是无人机动力系统设计里最常用的第一版性能评估。2. 叶片单元动量理论拆解动量视角和叶素视角如何合体2.1 动量理论控制体滑流模型给出一条约束方程动量理论看的是整个桨盘和它身后的滑流。对半径r处一个宽度dr的环形控制体来流速度V∞桨盘平面处的轴向速度是V∞vi远后方滑流最终减速到V∞2vi其中vi就是轴向诱导速度。质量流量是ρ乘以环形面积再乘以桨盘处的速度也就是ρ·2πr·dr·(V∞vi)。滑流速度增量是2vi因此环形控制体获得的轴向动量增量是质量流量乘以2vi这就是这个环带产生的拉力增量dT 4πrρ(V∞vi)·vi·dr这个式子干净、物理图像清楚但它只描述了“拉力与诱导速度之间的关系”并没有说清楚诱导速度到底由桨叶几何决定。换句话说动量理论给的是一个约束而不是答案。从动量理论出发还能得到理想状态下消耗于诱导的那部分功率dP dT(V∞vi)。悬停特例中V∞0dT 4πrρvi²dr全盘积分后能推出悬停时的理想诱导功率P sqrt(T³/(2ρA))A是桨盘面积。这是很多教材里出现的著名关系它也是BEMT低阶估算的参考底线真实螺旋桨效率永远低于这个理想值因为还有型阻功率和旋转尾流的额外损失。2.2 叶素理论把桨叶拆成一条条小机翼叶素理论的做法是把桨叶沿展向切成N个薄片每一片近似成一个二维翼型。每个小翼型看到的来流是两个分量的合成轴向分量V∞vi旋转分量Ωr-vt其中vt是当地切向诱导速度。由此可以得到当地入流角φ atan((V∞vi)/(Ωr-vt))当地攻角α θ(r) - φ。θ(r)是桨叶在当地的几何桨距角也就是零攻角位置和旋转平面之间的夹角实用中直接用扭转分布带入。有了攻角就能从翼型极曲线插值得到升力系数Cl(α)和阻力系数Cd(α)再算当地相对合速度Vrel sqrt((V∞vi)² (Ωr-vt)²)于是每片叶素单位展长上的升力和阻力都有了解析式dL 0.5ρVrel²·c·Cl(α)drdD 0.5ρVrel²·c·Cd(α)dr把升力沿轴向和切向投影分别得到单个叶片的推力和扭矩贡献dT dL·cosφ - dD·sinφdQ (dL·sinφ dD·cosφ)·r注意这里是单根叶片的量实际计算整个环形带的推力时要乘上桨叶数B。翼型极曲线的质量直接决定这条链路上结果的可靠性。风洞试验数据、XFOIL算出来的结果、CFD二维扫角结果都可以用但必须确认雷诺数和实际飞行工况大致匹配。2.3 联立求解两个理论给同一个量两个表达式迭代闭合现在关键的问题出现了dT和dQ的表达式里都含有诱导速度vi和vt而诱导速度本身又不独立它也受推力、扭矩影响。动量理论和叶素理论分别给了dT的两个表达式一个是来自叶素的dT_b B(dL·cosφ - dD·sinφ)另一个是来自动量的dT_m 4πrρ(V∞vi)vi·F dr。把两者联立就能解出满足自洽条件的诱导速度。数值上最常用的做法是先给vi和vt一个初值用叶素理论算出dT_b再把dT_b代入动量方程反解出新的vi。具体地说由4πrρF·vi² 4πrρF·V∞·vi - dT_b 0取正根vi_new [-V∞ sqrt(V∞² dT_b/(πrρF))] / 2其中F是Prandtl桨尖损失因子F (2/π)·arccos(exp(-f))f (B/2)(R-r)/(r·sinφ)。这个修正的意义在于桨尖处上下表面压差无法维持越靠近桨尖环带实际能承载的载荷越低。如果忽略F推力会被系统性高估效率曲线在高速段的形状也会失真。得到vi_new后用松弛因子缓慢更新重复上面的流程直到诱导速度不再变化。这一步是整个BEMT实现的核心也是最容易写错的环节。计算顺序稍微颠倒或者初始值选得太离谱迭代就会发散尤其是在V∞接近0的时候。2.4 BEMT关键假设与适用范围BEMT能成为工程主力是因为它在“计算成本”和“物理合理度”之间找到了平衡点。但也正因如此它的几个假设必须时刻放在心上。第一假设每个叶素是二维翼型忽略了展向流动和三维失速延迟第二假设流动定常、不可压马赫数不能太高第三假设桨盘载荷沿环带均匀分布对高实度比或强锥度桨会有偏差第四没有考虑桨叶之间的非定常干扰和桨尖涡的卷起。这些假设决定了BEMT最适合做参数扫描、方案对比和早期性能包输出。一旦进入详细设计、气动噪声评估或者强失速工况就需要更高保真度的手段来复核关键工作点。这不是BEMT的缺点而是任何低阶工具的边界。明白边界在哪里才能放心地在边界内用它。3. Matlab实现细节从几何数据到性能曲线3.1 参数输入弦长、扭转和翼型极曲线怎么组织Matlab里做这件事最直接的办法是把几何量离散成若干个径向站位每个站位对应弦长、扭转角、叶素所在半径。下面这组参数是我常用的起始模板r linspace(0.15R, 0.95R, 40)chord 0.022 - 0.008(r/R)twist 28(1-r/R)再换算成弧度。网格起点放在15%R而不是从桨毂开始是因为桨根附近的流动受到毂体、旋转圆柱效应和根部分离涡影响很大BEMT的二维假设在那里基本失效。从15%或20%半径开始起步数值上更稳工程上也更合理。翼型极曲线通常以表格形式存成三列攻角、升力系数、阻力系数。扫描攻角范围至少要覆盖桨叶可能出现的全部攻角区间。实际计算中要特别留意interp1插值在攻角越界时的行为。我一般把攻角限制在表格范围内并两端饱和而不是直接使用外推。因为升力曲线在失速后是非线性的默认的线性外推会得到明显离谱的升力系数。3.2 核心迭代循环一次一个叶素地更新诱导速度下面这段代码是主循环的骨架把注释读进去就能复现。为了让第一次接触的人不被细节淹没这里先固定切向诱导速度vt只让轴向诱导速度vi参与迭代低载荷粗算阶段这个简化是可接受的。要纳入完整切向诱导可以在同一循环里再加一条切向动量方程思路完全一样。% 基础参数 rho 1.225; % 空气密度 kg/m^3 R 0.254; % 桨叶半径 m D 2*R; % 直径 m B 2; % 桨叶数 RPM 6000; % 恒定转速 rpm Omega RPM*2*pi/60; % 角速度 rad/s n RPM/60; % 转每秒 % 叶素网格15%R 到 95%R40个站位 N 40; r linspace(0.15*R, 0.95*R, N); dr r(2) - r(1); % 几何分布示例换成你自己的螺旋桨即可 chord 0.022 - 0.008*(r/R); % 弦长分布 m twist (28*(1 - r/R)) .* pi/180; % 桨距角分布 rad % 翼型极曲线示例线性升力二次阻力 alpha_tab linspace(-0.35, 0.55, 100); cl_tab 2*pi*sin(alpha_tab); % 薄翼线性段近似 cd_tab 0.02 0.25*(alpha_tab - 0.05).^2; % 前进比扫描 J_list 0.05:0.05:0.8; CT zeros(size(J_list)); CP zeros(size(J_list)); eta zeros(size(J_list)); for k 1:length(J_list) J J_list(k); Vinf J * n * D; % 来流速度 % 诱导速度初值不能从0开始 vi ones(N,1); % 轴向诱导初值 vt 0.02*Omega*r; % 切向诱导初值本骨架中固定 relax 0.25; for it 1:400 vi_old vi; dT zeros(N,1); dQ zeros(N,1); % 叶素理论计算每个站位的推力和扭矩 for i 1:N phi atan2(Vinf vi(i), Omega*r(i) - vt(i)); alpha twist(i) - phi; cl interp1(alpha_tab, cl_tab, alpha, linear, extrap); cd interp1(alpha_tab, cd_tab, alpha, linear, extrap); Vrel sqrt((Vinf vi(i))^2 (Omega*r(i) - vt(i))^2); dL 0.5*rho*Vrel^2*chord(i)*cl*dr; dD 0.5*rho*Vrel^2*chord(i)*cd*dr; dT(i) B*(dL*cos(phi) - dD*sin(phi)); dQ(i) B*(dL*sin(phi) dD*cos(phi))*r(i); end % 动量理论反解轴向诱导速度逐叶素更新 for i 1:N phi_i atan2(Vinf vi(i), Omega*r(i) - vt(i)); f_tip (B/2)*(R - r(i))/(r(i)*sin(phi_i)); F_tip (2/pi)*acos(exp(-f_tip)); vi_new 0.5*(-Vinf sqrt(Vinf^2 dT(i)/(pi*r(i)*rho*F_tip))); vi(i) vi(i) relax*(vi_new - vi(i)); end % 收敛判据 if max(abs(vi - vi_old)) 1e-4 break; end end % 积分得到总力和功率 T sum(dT); Q sum(dQ); P Q*Omega; % 无量纲系数 CT(k) T/(rho*n^2*D^4); CP(k) P/(rho*n^3*D^5); eta(k) J*CT(k)/CP(k); end % 画图 figure; plot(J_list, CT, o-); xlabel(J); ylabel(CT); figure; plot(J_list, eta, s-); xlabel(J); ylabel(Efficiency);这个骨架里升力系数用的是薄翼线性段近似阻力的二次型也是典型的估算值。实际应用时把alpha_tab、cl_tab、cd_tab换成自己的极曲线表即可。两个细节初始诱导速度不能给0尤其从低速工况开始计算时给0会让第一次迭代的攻角异常大容易直接飞掉松弛因子取0.25左右比较稳太大会在低前进比段震荡太小则收敛缓慢。3.3 无量纲化输出与结果验证三个无量纲量的定义按常规螺旋桨标准来CT T/(ρn²D⁴)CP P/(ρn³D⁵)η J·CT/CP这三个量把转速、直径和空气密度的影响都归一化了不同条件下的计算结果可以直接比较。举例说一个直径为0.508m的桨以6000RPM旋转n 100转每秒桨盘面积约0.203平方米在J0.3时对应的来流约15m/s属于低速巡航J0.7对应来流35m/s接近小飞机的高速飞行。写完代码后先别急着扫参数做几步验证。第一步是网格独立性检查把叶素数量从40翻到80看推力和效率的变化是否在1%以内如果变化明显就说明离散太粗。第二步是看悬停特例把J设成一个很小的值对比总推力与动量理论给出的悬停估计是否处于同一量级。第三步如果有条件拿一组公开的螺旋桨风洞数据对一下CT和CP这是最令人安心的验证。我经常看到有人直接拿着代码去算连网格都没检查最后曲线形状怪异却找不到原因问题往往就出在这类基础验证没做。4. 前进比扫描结果分析与常见问题排查4.1 曲线怎么读推力、功率、效率随J的变化规律把CT对J画出来典型趋势是单调下降。J很小的时候桨叶攻角大每个叶素都在高升力区工作推力系数很高J变大使来流分量增大入流角φ变大攻角α θ - φ变小升力下降阻力占比上升CT随之下降。CP的下降通常比CT缓和因为扭矩还包含与诱导速度耦合的部分。效率曲线则是先升后降低J段型阻损耗占比大高J段推力本身变小效率也不划算。这些曲线对设计判断很有用。如果某型多旋翼在悬停设计点J0.05附近工作效率低不要紧因为悬停本来就在高负荷低效率区但巡航机在前进比0.40.6范围里效率拉不起来就需要重新考虑桨距角设计或更换螺旋桨。下面给一个趋势示意表不是具体某把桨的真实数据但形状是一致的JCTCPη0.050.0950.0620.0770.150.0880.0550.2400.350.0610.0410.5210.550.0380.0320.6530.700.0210.0270.544我经常提醒自己读曲线时不要只看总推力在多少J变为负要先看哪个半径开始“拖后腿”。当J增加到很大时内侧叶素会先出现攻角极低甚至为负的情况局部进入负推力状态。这部分叶素不仅不贡献拉力还增加阻力导致整体效率急剧下滑。这也是为什么BEMT要逐叶素输出结果而不是只看总积分。4.2 常见问题速查表发散、越界、负拉力下面这些是我在不同项目里反复遇过的问题整理成一个速查表。遇到曲线异常时对照着查大多数情况几分钟就能定位。现象可能原因处理办法低前进比时迭代发散残差越来越大诱导速度初值为0第一次攻角过大或松弛因子太大初值按悬停估算给relax降到0.2左右J从0.05起步interp1返回NaN或出现荒谬升力系数攻角超出翼型极曲线表范围用两端饱和代替extrap或扩展极曲线攻角范围J增大后内侧单元出现负推力入流角过大导致局部负攻角物理真实存在保留负推力不要强制归零检查扭转是否过小推力系数对网格密度敏感叶素数量少于20或网格分布在桨尖太粗用40到60个叶素做一次网格加倍测试确认变化1%计算出的效率在大J段异常升高忽略了Prandtl桨尖损失高估了桨尖单元贡献加入F修正观察趋势是否恢复正常这些坑的处理思路本质都是“不要让数值方法在物理失效的边界上继续工作”。BEMT是工程模型不是万能求解器遇到边界工况宁可把工况范围缩小也不要硬算一个不可信的数字出来。4.3 BEMT的边界以及我现在的固定使用习惯BEMT能做的事和不能做的事要心里有数。它假定每个叶素是二维翼型忽略了展向流动、旋转导致的失速延迟、桨尖涡的卷起以及动态失速这类非定常效果。在低前进比、大攻角工况下这些效应往往同时出现单纯靠BEMT算出的悬停效率会偏乐观误差可能达到10%到20%。所以我的使用习惯是BEMT用于方案筛选、参数扫参、初步性能包输出一旦进入详细设计再用更高保真度的方法复核关键工作点。还有一点值得单独提J接近0的悬停点是BEMT的天然薄弱点。我经常遇到从J0直接算就发散、从J0.05算完再手动外推反而稳定的情况。原因在于V∞0时动量方程退化迭代方程对初值极其敏感。现在我一般把扫描序列写成J_list 0.05:0.05:0.8需要悬停数据时用低J段的拟合曲线外推而不是硬着头皮去算J0。换翼型极曲线时也有一个细节实际螺旋桨用的翼型在失速前升力斜率往往低于2π雷诺数变化也会让极曲线整体平移。最好根据不同半径处的雷诺数选择对应雷诺数的极曲线表而不是全桨共用一条。这在低雷诺数小桨上尤其重要因为桨根和桨尖的弦长、速度差得很远雷诺数能差出好几倍。把这步做好了BEMT结果的可靠性会明显上一个台阶。最后说一个我自己的固定操作跑完所有前进比之后我会额外保存一组“诱导速度沿半径分布”的中间结果。它虽然只是中间量却能告诉我哪个半径范围正在高速消耗功率。如果某段叶素始终处于高阻力低升力状态下一步就该修改弦长或扭转分布。这个数据比最终的CT和CP曲线更能指导几何优化。