做螺旋桨选型或者气动评估的时候我很少直接拿厂家给的那条效率曲线当唯一依据。原因很简单厂家曲线通常是在某一个转速、某一个来流速度下测出来的换一个使用环境桨的拉力和效率完全就是另一回事。所以我习惯用BEMTBlade Element Momentum Theory动量理论剖面刀片方法在Matlab里搭一套快速分析程序给定桨叶几何、恒定转速和一系列前进比把这条桨的性能曲线自己算出来。BEMT说白了就是动量理论和叶素理论各取所长拼出来的工程方法计算量小、收敛快特别适合做参数扫掠很多飞行器气动设计和飞控评估的前期工作都靠它。这篇文章就把我的实现思路、代码骨架和踩过的坑完整讲清楚想自己搭一套BEMT工具的朋友可以直接参照。1. 先搞清楚BEMT要回答什么问题前进比与恒定转速的关系1.1 前进比、推力系数和效率三个无量纲量先立住BEMT算出来的不是绝对推力而是无量纲系数这是整个工具的核心。螺旋桨气动分析里最常用的两个系数是推力系数C_T和功率系数C_P定义为C_T T / (ρ n^2 D^4) C_P P / (ρ n^3 D^5)其中T是推力P是轴功率ρ是空气密度n是螺旋桨转速转每秒D是螺旋桨直径。引入前进比J V / (n D)V是来流速度。这样螺旋桨效率就是η J * C_T / C_P。这个公式非常优雅它把“飞行速度带来的收益”和“消耗的轴功率”放在一起本质上是一个无因次的推进效率。前进比越大说明同样的转速和直径下飞行速度越快但这不意味着效率越高——后面你会看到J过了峰值以后效率反而往下掉。1.2 为什么是“恒定转速扫前进比”而不是扫油门实际飞行中电调或发动机控制器通常会维持螺旋桨转速在某个恒定值附近这是最常见的工作模式。比如多旋翼悬停、固定翼巡航飞控给的是转速指令电机转速基本稳定但飞行速度会随油门、气流、姿态变化。这时候转速不变、直径不变来流速度V一变前进比J就跟着变桨叶每个剖面的有效攻角也全部变了。所以工程上最自然的做法就是固定转速n扫一组V换算成J看C_T、C_P、η怎么变。这也正是标题里“不同前进比下以恒定转速”这句话的物理来源。单独给一个J点不算数一条曲线才有说服力。2. 动量理论加叶素理论为什么非得“拼”在一起2.1 动量理论的“圆盘”约束滑流把诱导速度定下来了动量理论把螺旋桨简化成一个圆盘气体穿过圆盘时速度从远前方来流V变成盘面处的Vv_i再到远后方变成V2v_iv_i就是轴向诱导速度。它告诉我们是这个诱导速度让气流带走了动量和能量从而产生了推力和扭矩。按动量方程推导在半径r处、宽度dr的圆环上推力微元可以写成dT_mom 4π ρ r V^2 F a (1a) dr这里的a是轴向诱导因子v_i aVF是桨尖损失修正。旋转方向也有类似关系引入切向诱导因子a扭矩微元为dQ_mom 4π ρ r^3 V Ω F a (1a) dr其中Ω是角速度。这一套公式非常“强”它不管桨叶长什么样只关心圆盘从气流里抽走了多少动量。但问题是它给不出a和a到底是多少必须靠叶素理论从桨叶的几何和翼型性能反推。2.2 叶素理论的“切片”受力攻角从哪来把桨叶沿着展向切成很多小段每一段半径r处看当地气流合速度。轴向分量是V(1a)切向分量是Ωr(1-a)。于是当地气流角φ满足φ atan2( V(1a), Ω r(1-a) )这个φ非常关键它等于来流速度方向与旋转平面的夹角。桨叶在r处的几何桨距角是θ当地攻角就是α θ - φ。然后查翼型升力系数C_l和阻力系数C_d得到单位展长的力和力矩dT_bem 0.5 ρ W^2 c (C_l cosφ - C_d sinφ) dr dQ_bem 0.5 ρ W^2 c r (C_l sinφ C_d cosφ) drW是合速度的大小。你马上会发现一个循环关系要算攻角需要知道a和a要算a和a又需要叶素算出来的力力又依赖攻角。这就是BEMT要迭代解的核心原因。只有叶素理论攻角是猜的只有动量理论载荷是空的两边必须联立。2.3 联立迭代a和a就是两个未知数的闭合方程组工程上最常用的办法是把动量公式倒过来。叶素理论算出一个dT_bem动量理论要求dT_mom等于这个值代入1a和a的关系式就能解出a。整理一下令 X dT_bem / (4π ρ r V^2 F)则 a_new 0.5 * (sqrt(1 4X) - 1)同理a_new dQ_bem / (4π ρ r^3 V Ω F (1a))这里每一个半径r处都有一组(a, a)再加上径向的一整圈离散点就变成了一个分布式的非线性系统。实际写代码时我用松弛迭代先给初值a0、a0算叶素力更新a和a再拿新的a和a重新算攻角如此反复直到变化量小于容差。这个过程很像双臂互搏迭代十几二十步通常就稳定了。3. Matlab实现前的数据准备几何离散和翼型极线3.1 桨叶几何参数化半径、弦长、扭转怎么进结构体在写主循环之前先把输入数据组织好。我用一个结构体geo保存桨叶几何hub_radius是轮毂半径R是桨尖半径B是桨叶数r是各叶素半径chord是各站弦长theta是各站几何桨距角。r通常从轮毂到桨尖取30到50个点。离散方式我个人喜欢直接用linspace简单直接但如果你想让桨尖和轮毂附近的分辨率更高可以用余弦分布。余弦分布本质是把有剧烈流动变化的区域加密桨尖处有桨尖损失和强三维效应轮毂处速度低、载荷占比小所以优先加密桨尖。代码里判断几何是否合理时我习惯把弦长和扭转画出来看一眼很多时候错误一眼就能暴露比如弦长在某站跳变那就是数据给错了。3.2 翼型极线查表插值方法和大攻角外推的边界翼型数据是一张表包含攻角α度、升力系数C_l、阻力系数C_d。Matlab里直接interp1线性插值就行。这里有两个坑一是攻角单位必须和几何扭转一致我全部统一用度进主循环前算好θ反正代码里φ用atan2得到的是弧度转成度以后才能和θ相减二是角度可能超出数据表范围。静拉力附近每个叶素的攻角非常接近失速区甚至更大高前进比时攻角又会很负线性外推很容易推出一个离谱的C_l。稳妥做法是先用‘linear’插值越界时把C_l按失速后的趋势压住或者干脆截断在边界值。我更推荐补一组大攻角数据哪怕粗糙一点也比外推靠谱。某次我在一个项目里偷懒用了外推静拉力计算结果比台架测试高了15%后来排查半天才发现是这个原因。3.3 单位制约定一个能毁掉整个程序的隐藏坑许多第一次写BEMT的人死在这步。我强调三件事转速必须用转每秒不是RPM长度必须是米密度必须和高度匹配。转速出现RPM而没除以60是最常见的错误会让前进比全部差60倍推力系数差3600倍。另一个坑是转速单位用“转每分钟”但直径用厘米然后算出来的系数更像随机数。所以我把入口参数设计成明确物理量n输入转每秒V输入米每秒R输入米rho输入kg/m^3。如果是海平面rho1.225标准情况下直接用。想换高度就换密度而不是改速度。单位制检查我会在程序开头加assert输入n100时自动警告一次能省掉大量debug时间。4. 主程序骨架恒定转速、扫前进比的完整循环4.1 外层循环用一组来流速度变成一整条性能曲线主函数核心框架是两层循环。外层循环遍历来流速度数组V_list每个V对应一个前进比J最终得到J_vec、CT、CP、eta四个数组。内层是每个J下的诱导速度迭代。先看代码骨架function [J_vec, CT, CP, eta] bem_scan(geo, airfoil, n, V_list, rho) % BEMT主扫描函数 % geo: 桨叶几何结构体包含hub_radius, R, B, r, chord, theta % airfoil: 翼型极线表包含alpha(deg), cl, cd % n: 恒定转速转/秒 % V_list: 来流速度数组m/s % rho: 空气密度 Omega 2 * pi * n; R geo.R; B geo.B; r geo.r(:); % 叶素半径列向量 chord geo.chord(:); theta geo.theta(:); % 几何桨距角单位度 dr diff([geo.hub_radius; r]); % 每个叶素对应的展向宽度 relax 0.4; tol 1e-6; J_vec zeros(size(V_list)); CT zeros(size(V_list)); CP zeros(size(V_list)); eta zeros(size(V_list)); for k 1:numel(V_list) V max(V_list(k), 0.05); % 避免静拉力奇点见正文说明 a zeros(size(r)); ap zeros(size(r)); for iter 1:800 a_old a; ap_old ap; T 0; Q 0; for i 1:numel(r) phi atan2(V * (1 a(i)), Omega * r(i) * (1 - ap(i))); alpha theta(i) - phi * 180 / pi; [Cl, Cd] get_airfoil(airfoil, alpha); W sqrt((V * (1 a(i)))^2 (Omega * r(i) * (1 - ap(i)))^2); dT 0.5 * rho * W^2 * chord(i) * (Cl * cos(phi) - Cd * sin(phi)); dQ 0.5 * rho * W^2 * chord(i) * r(i) * (Cl * sin(phi) Cd * cos(phi)); f B * (R - r(i)) / (2 * r(i) * sin(phi)); F (2 / pi) * acos(exp(-f)); F max(F, 0.05); % 防止桨尖处除零 X dT / (4 * pi * rho * r(i) * V^2 * F); a_new 0.5 * (sqrt(1 4 * X) - 1); ap_new dQ / (4 * pi * rho * r(i)^3 * V * Omega * F * (1 a(i))); a(i) a_old(i) relax * (a_new - a_old(i)); ap(i) ap_old(i) relax * (ap_new - ap_old(i)); T T dT * dr(i); Q Q dQ * dr(i); end if max(abs(a - a_old)) tol max(abs(ap - ap_old)) tol break; end end J_vec(k) V / (2 * n * R); CT(k) T / (rho * n^2 * (2 * R)^4); CP(k) Q * Omega / (rho * n^3 * (2 * R)^5); eta(k) J_vec(k) * CT(k) / CP(k); end end这段代码已经把主要结构全部包含了。你注意到外层循环里我加了V max(V_list(k), 0.05)就是用一个小速度代替严格的零速度避免静拉力点动量方程除零。后面第6节我们再细说这个处理。4.2 内层迭代诱导速度更新为什么不能一步到位拿到dT之后直接算a_new按理说把a_new直接赋给a不就行了吗不行。因为dT是用旧的a和a算出来的旧值本身含有误差一步到位会让迭代像过山车一样发散。实际中高载荷工况下尤其静拉力附近a可以到0.3甚至0.5直接替换会在几个迭代步内震荡到负数。所以我在两个更新式子里都加了松弛因子relax当前解只往目标方向走40%。这个做法不是最快的但它稳。等代码跑通以后你可以再研究自适应松弛、牛顿迭代加速没必要第一次就追求最高效率。翼型查表函数很直接function [Cl, Cd] get_airfoil(airfoil, alpha) Cl interp1(airfoil.alpha, airfoil.cl, alpha, linear, extrap); Cd interp1(airfoil.alpha, airfoil.cd, alpha, linear, extrap); end这里用的是线性外推只适合做演示。实际使用前我强烈建议把airfoil表补成大攻角范围比如从-20度到30度并且让失速后C_l随攻角增大而下降而不是线性外推继续涨。翼型极线表越完整BEMT才越可信。4.3 积分推力扭矩dr怎么算决定结果一致性内层循环里每个叶素都算出了单位展长的dT和dQ最后要乘以展向宽度dr再累加。这里有个细节我用的dr是diff([hub_radius; r])意思是r数组第一项是轮毂半径然后每个r点之间的间隔作为该点的展向宽度。如果你用均匀rdr就是常数。注意最靠近轮毂的那一段宽度是r(1)-hub_radius这样扭矩积分不会丢面积。有些写法直接取R/N如果离散不均匀就错了。另一个细节是dT和dQ必须对应“单位展长”而不是乘以dr之后再代入动量公式求a_new否则量纲就乱了。我习惯在代码注释里写清楚每个变量的量纲这个习惯救过我很多次。一轮迭代结束后用累加出来的T和Q计算无量纲系数。这里的P是轴功率不是有用功率所以用QOmega。η则是用JCT/CP这个式子和无量纲定义完全一致不需要再手算一遍。5. 最容易翻车的三个数值细节失速、桨尖损失和松弛5.1 大攻角查表越界的连锁反应静态或低速来流下螺旋桨剖面的攻角经常超过翼型失速角。很多翼型极线表只给到-5度到15度但BEMT迭代过程中攻角可能到25度甚至更高。线性外推在超过15度后还让升力继续上升这在物理上完全错误真实薄翼在失速后升力会下降。这个错误会直接导致静拉力和小J工况的推力被高估扭矩也可能偏小效率当然对不上台架试验。我的建议是两件事第一给翼型极线表增加大攻角数据只要C_l在失速后有明确回落哪怕粗糙也比外推强第二在get_airfoil里加一个越界标记如果alpha超过数据范围就在结果里加warning。不要悄悄截断至少要让使用者知道这个工况已经超出常规翼型数据边界。5.2 Prandtl桨尖修正的边界行为桨尖损失因子F的公式是f B(R-r)/(2r sinφ)F (2/π) acos(exp(-f))。物理上它描述桨尖附近因为上下表面压差沟通而损失的载荷。问题是当r逼近R时R-r趋近0f趋近0exp(-f)趋近1acos(1)0F趋近0。在动量方程里F在分母上F趋近0会让a_new变成无穷大。我的处理是给F加一个下限0.05相当于认为桨尖有效载荷衰减到5%就不要再细算。这个方法简单粗暴但工程上足够稳定。还要注意sinφ在小半径处可能很小f会变得很大这时exp(-f)趋近0acos(0) π/2F趋近1反而没问题。实际调试时如果发现靠近桨尖的a出现跳动优先怀疑F的下限设置和离散密度。5.3 松弛因子怎么选不只用0.4还要看残差趋势松弛因子选大了发散选小了收敛慢。我经验上的起点是0.3到0.5低载荷巡航工况可以放宽到0.7但高载荷静拉力附近建议压到0.2。怎么判断要不要调我习惯在迭代循环里同时打印max(abs(a-a_old))和推力残差如果残差一直震荡不下降就把relax减半。还有一种情况是a在桨尖附近震荡这说明F处理可能还有问题不单纯是松弛的事。此外内外两层循环不要都用同一个松弛因子。诱导速度a和a的收敛速度不同a通常更慢尤其在切向诱导速度贡献大的桨根区域。我曾经遇到迭代500步后a已经稳定但a还在漂移的情况后来把a和a的松弛因子分开设置a用0.4a用0.3问题立刻解决。6. 算例结果怎么读从效率曲线判断螺旋桨的工作区间6.1 一组典型结果CT、CP、η随前进比的走向下面这组数是某型两叶桨在恒定转速下的示意计算结果桨叶几何是我随手构造的不代表任何真实产品但趋势非常典型。JC_TC_Pη00.120.1800.10.1000.1300.080.20.0850.1100.150.40.0550.0750.290.60.0300.0500.360.80.0180.0400.36J0时的η为0因为飞行速度为0推力做的是零功但电机依然要消耗功率维持螺旋桨转动。从0到0.4推力系数一路往下掉但效率在爬升说明飞行速度带来的收益开始体现。到J0.6到0.8效率基本到顶之后如果继续增大J推力会越来越小效率也会下降。所以这条曲线最有价值的点就是最高效率对应的J区间也就是所谓的设计巡航点。做飞机设计时我会根据巡航速度和转速反推需要的J然后看这个J是不是落在效率峰值附近不是的话就要改桨距或换桨。6.2 和实验数据对比时要注意什么BEMT计算和实际台架测试的差距很多时候不完全是BEMT本身的问题而是测试状态和计算状态对不上。台架测试通常有支架阻力、电机散热气流、桨毂遮挡这些东西都会影响实际进流。另一个通常被忽略的是雷诺数。BEMT里用的翼型极线一般来自风洞数据但真实桨叶叶素雷诺数可能只有几十万和风洞数据状态不一定一致。C_d偏小算出来的效率就偏高。所以对比实验时我习惯在表格里同时记录转速、来流速度和环境密度然后把C_T和C_P对不上先归因于翼型数据而不是急着改几何参数。很多时候把翼型极线换成同型号更接近真实雷诺数的数据误差能缩小一半。6.3 静拉力点J0的特殊处理一个小速度的工程折中严格意义上V0时动量方程里的V^2会变成0a_new公式直接除零程序会崩。但这个点在实际工程里又非常重要很多多旋翼选电机的第一步就是看静拉力够不够。我的做法是用V0.05 m/s代替0这样前进比J其实是一个很小但不是零的数计算出来的C_T可以近似当作静拉力系数。为什么不用更小的0.001因为数值上会让诱导速度迭代变得非常困难收敛极慢。0.05在大多数螺旋桨尺寸下对应的J只有0.001量级对结果影响远小于数值误差。更精细的工程做法是单独写一个静拉力求解函数用专门处理悬停状态的动量理论版本但作为初版工具小速度替代几乎零成本值得保留。最后再说一个我实际操作中的体会BEMT这套东西最值钱的不是那几个方程而是你肯花时间把几何、翼型极线、单位制、数值稳定这些边角料都理清楚。代码本身不到一百行但每一次结果对不上试验几乎都是数据准备或者边界处理出了问题。如果后面你还想做得更细可以在此基础上加入雷诺数随半径变化的修正把三维修正项叠进去或者反过来用一维优化器扫描最佳弦长分布这个骨架都是现成的。别急着追求花哨先把恒定转速下的性能曲线稳定跑出来你就已经赢过大多数拍脑袋看标签选桨的人了。