做螺旋桨选型或者气动评估的时候，我很少直接拿厂家给的那条效率曲线当唯一依据。原因很简单：厂家曲线通常是在某一个转速、某一个来流速度下测出来的，换一个使用环境，桨的拉力和效率完全就是另一回事。所以我习惯用BEMT（Blade 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变成盘面处的V+v_i，再到远后方变成V+2v_i，v_i就是轴向诱导速度。它告诉我们是这个诱导速度让气流带走了动量和能量，从而产生了推力和扭矩。按动量方程推导，在半径r处、宽度dr的圆环上，推力微元可以写成：

dT_mom = 4π ρ r V^2 F a (1+a) dr

这里的a是轴向诱导因子，v_i = aV，F是桨尖损失修正。旋转方向也有类似关系，引入切向诱导因子a'，扭矩微元为：

dQ_mom = 4π ρ r^3 V Ω F a' (1+a) dr

其中Ω是角速度。这一套公式非常“强”，它不管桨叶长什么样，只关心圆盘从气流里抽走了多少动量。但问题是，它给不出a和a'到底是多少，必须靠叶素理论从桨叶的几何和翼型性能反推。

2.2 叶素理论的“切片”受力：攻角从哪来

把桨叶沿着展向切成很多小段，每一段半径r处，看当地气流合速度。轴向分量是V(1+a)，切向分量是Ωr(1-a')。于是当地气流角φ满足：

φ = atan2( V(1+a), Ω 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φ) dr

W是合速度的大小。你马上会发现一个循环关系：要算攻角，需要知道a和a'；要算a和a'，又需要叶素算出来的力；力又依赖攻角。这就是BEMT要迭代解的核心原因。只有叶素理论，攻角是猜的；只有动量理论，载荷是空的；两边必须联立。

2.3 联立迭代：a和a'就是两个未知数的闭合方程组

工程上最常用的办法是把动量公式倒过来。叶素理论算出一个dT_bem，动量理论要求dT_mom等于这个值，代入1+a和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 (1+a))

这里每一个半径r处都有一组(a, a')，再加上径向的一整圈离散点，就变成了一个分布式的非线性系统。实际写代码时，我用松弛迭代：先给初值a=0、a'=0，算叶素力，更新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。如果是海平面，rho=1.225，标准情况下直接用。想换高度就换密度，而不是改速度。单位制检查我会在程序开头加assert，输入n>100时自动警告一次，能省掉大量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点之间的间隔作为该点的展向宽度。如果你用均匀r，dr就是常数。注意最靠近轮毂的那一段，宽度是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趋近0，f趋近0，exp(-f)趋近1，acos(1)=0，F趋近0。在动量方程里F在分母上，F趋近0会让a_new变成无穷大。我的处理是给F加一个下限0.05，相当于认为桨尖有效载荷衰减到5%就不要再细算。这个方法简单粗暴，但工程上足够稳定。还要注意sinφ在小半径处可能很小，f会变得很大，这时exp(-f)趋近0，acos(0) = π/2，F趋近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.4，a'用0.3，问题立刻解决。

6. 算例结果怎么读：从效率曲线判断螺旋桨的工作区间

6.1 一组典型结果：CT、CP、η随前进比的走向

下面这组数是某型两叶桨在恒定转速下的示意计算结果，桨叶几何是我随手构造的，不代表任何真实产品，但趋势非常典型。

J=0时的η为0，因为飞行速度为0，推力做的是零功，但电机依然要消耗功率维持螺旋桨转动。从0到0.4，推力系数一路往下掉，但效率在爬升，说明飞行速度带来的收益开始体现。到J=0.6到0.8，效率基本到顶，之后如果继续增大J，推力会越来越小，效率也会下降。所以这条曲线最有价值的点就是最高效率对应的J区间，也就是所谓的设计巡航点。做飞机设计时，我会根据巡航速度和转速反推需要的J，然后看这个J是不是落在效率峰值附近，不是的话就要改桨距或换桨。

6.2 和实验数据对比时要注意什么

BEMT计算和实际台架测试的差距，很多时候不完全是BEMT本身的问题，而是测试状态和计算状态对不上。台架测试通常有支架阻力、电机散热气流、桨毂遮挡，这些东西都会影响实际进流。另一个通常被忽略的是雷诺数。BEMT里用的翼型极线一般来自风洞数据，但真实桨叶叶素雷诺数可能只有几十万，和风洞数据状态不一定一致。C_d偏小，算出来的效率就偏高。所以对比实验时，我习惯在表格里同时记录转速、来流速度和环境密度，然后把C_T和C_P对不上先归因于翼型数据，而不是急着改几何参数。很多时候把翼型极线换成同型号更接近真实雷诺数的数据，误差能缩小一半。

6.3 静拉力点J=0的特殊处理：一个小速度的工程折中

严格意义上，V=0时动量方程里的V^2会变成0，a_new公式直接除零，程序会崩。但这个点在实际工程里又非常重要，很多多旋翼选电机的第一步就是看静拉力够不够。我的做法是用V=0.05 m/s代替0，这样前进比J其实是一个很小但不是零的数，计算出来的C_T可以近似当作静拉力系数。为什么不用更小的0.001？因为数值上会让诱导速度迭代变得非常困难，收敛极慢。0.05在大多数螺旋桨尺寸下对应的J只有0.001量级，对结果影响远小于数值误差。更精细的工程做法是单独写一个静拉力求解函数，用专门处理悬停状态的动量理论版本，但作为初版工具，小速度替代几乎零成本，值得保留。

最后再说一个我实际操作中的体会：BEMT这套东西最值钱的不是那几个方程，而是你肯花时间把几何、翼型极线、单位制、数值稳定这些边角料都理清楚。代码本身不到一百行，但每一次结果对不上试验，几乎都是数据准备或者边界处理出了问题。如果后面你还想做得更细，可以在此基础上加入雷诺数随半径变化的修正，把三维修正项叠进去，或者反过来用一维优化器扫描最佳弦长分布，这个骨架都是现成的。别急着追求花哨，先把恒定转速下的性能曲线稳定跑出来，你就已经赢过大多数拍脑袋看标签选桨的人了。