螺旋桨BEMT性能计算：基于Matlab的拉力、功率与效率分析

这件事我琢磨了好一阵子：螺旋桨这东西，从外面看就是个一转就飞的“大风扇”，但真要把它的拉力、扭矩、效率随飞行状态的变化规律算明白，牵涉到的流体力学细节足够让初学者挠头。我最初接触叶片单元动量理论（Blade Element Momentum Theory，BEMT）是在大学做飞行器设计课设的时候，当时只想着“能出个结果交差就行”，结果被诱导速度那个迭代环节卡了整整一个周末。后来工作里做小型无人机的动力选型，才发现BEMT这套东西在工程上的用处比我想象中大得多：不需要跑复杂的CFD，不需要昂贵的许可证，只要把桨叶几何形状和翼型数据喂给Matlab，就能在几秒钟内得到螺旋桨在任意前进比下的性能曲线。这篇文章就把我这些年用Matlab实现BEMT的经验完整梳理一遍，从理论框架到代码细节再到调试心得，能帮你避掉大多数新手会踩的坑。

这个项目的目标一句话说清楚：给一个具体的螺旋桨几何（半径、弦长沿径向分布、扭转角沿径向分布、桨叶数、翼型气动数据），假设转速恒定不变，然后扫描一串前进比（advance ratio，J），算出每个J对应的拉力系数（C_T）、功率系数（C_P）和效率（η），最后画出曲线。前进比J是螺旋桨无因次分析里最核心的参数，定义为J = V / (n·D)，其中V是来流速度，n是转速（转/秒），D是桨盘直径。它直观地反映了“螺旋桨每转一圈实际向前推进的距离”和“直径”的比值。这活儿最适合用Matlab干，因为涉及迭代求解、数组索引、绘图，都是Matlab的舒适区。适合谁来读呢？飞行器气动专业的学生、无人机动力系统设计工程师、模型爱好者自制螺旋桨验证，甚至做风力机的人也能把同一套代码改吧改吧直接用。

1. 先拆解计算框架：BEMT到底在算一个什么方程

1.1 两类理论的拼接：动量守恒与叶素受力

BEMT的名字已经剧透了，它是两个经典理论的杂交体。动量理论（Momentum Theory）把螺旋桨看作一个致动盘（actuator disc），只关注气流通过桨盘前后的速度变化和压力差，好处是能给出理想效率和诱导速度的宏观关系，缺点是它不关心叶片长什么样子——同样一个拉力，它可以由宽桨叶低速旋转产生，也可以由窄桨叶高速旋转产生，动量理论完全分不出来。叶素理论（Blade Element Theory）则反过来，把每片桨叶沿展向切成无数小薄片，每一片当作一个二维翼型，用当地速度三角形算出攻角，再查翼型极曲线（Cl-α、Cd-α）得到升力和阻力，沿径向积分就得到整个桨叶的贡献。它非常依赖几何输入，但如果没有动量理论提供“自洽的诱导速度”，叶素理论里攻角怎么取就成了无源之水——你总不能用自由来流速度算攻角，那会让计算出的拉力大到离谱。

所以BEMT的精髓就是把这两个视角缝合起来：动量理论说“因为气流获得了动量，所以桨盘处的诱导速度是某个值”，叶素理论说“如果我采用这个诱导速度，桨叶能产生的升力是另一个值”，让这两个值互相逼近，解出每个径向位置的诱导速度，然后再反算载荷。这个自洽迭代的过程，就是整段代码的心脏。

1.2 为什么在工程初算阶段不用CFD

我知道肯定有人问：现在CFD都这么成熟了，为什么还要用BEMT？原因很现实。我在实际项目里做动力选型时，经常要在一天之内对比十几款候选螺旋桨，“同一款桨在不同转速、不同飞行速度下性能如何”这种矩阵式的扫描，如果用CFD来做，一个工况就可能要跑上几百核时，一周都不一定收敛。而BEMT把问题降维到了准一维的径向分布上，单工况计算时间以毫秒计，几百个工况的扫描也就几秒钟的事。精度上，对于正常设计的桨叶（没有严重的失速分离、没有跨音速效应），BEMT给出的效率和拉力系数和风洞实验的误差通常能控制在5%到10%以内，这个精度对工程预研和概念设计阶段完全够用。CFD当然能做BEMT看不到的细节——比如桨尖涡的螺旋结构、叶片表面流动分离形态——但那是后置的精细化验证，不是前期设计的工具。我的原则是：先用BEMT把设计空间探索清楚，锁定几个候选方案，最后用CFD甚至风洞实验确认，这才符合工程节奏。

1.3 前进比和恒定转速的意义

标题特意强调了“恒定转速”，这其实暴露了一个典型工程场景：多数中小型无人机用的是无刷电机加电子调速器，巡航时转速由飞控闭环稳定在某一固定值，飞行速度则随风况在变化。所以恒定转速下扫前进比，本质上是在回答“转速固定时，飞得快和飞得慢对螺旋桨效率有什么影响”。这比“固定速度扫转速”更贴近实际操纵方式。另外，从无因次角度讲，BEMT的输入参数几乎全部可以无因次化，在给定螺旋桨几何下，无因次系数只取决于前进比J和桨尖马赫数相关的雷诺数效应，所以“恒定转速”给了一个很好的基准。需要提前说明的是：严格说同一个转速在不同前进比下，桨叶各径向位置的雷诺数会随着合速度大小变化，但雷诺数对翼型极曲线的影响通常不大（除了小雷诺数下的层流分离泡问题），在工程初算里可以忽略或者用固定雷诺数的极曲线近似，这一点到后面讲翼型数据准备时会再提。

2. 理论公式与数值迭代细节的逐行拆解

2.1 几何输入：弦长、扭转角和翼型数据如何影响计算

给定一个螺旋桨几何形状，BEMT需要的最基本参数是这些：桨叶数B、半径R、弦长分布c(r)、几何扭转角分布θ(r)、翼型的升阻特性、来流速度V、转速n。其中弦长和扭转角的径向分布可以用离散点表示，Matlab里用两个等长的向量就行。我通常的做法是，把桨叶从毂半径r_hub到叶尖R等分成几十段，每一段取中点的弦长和扭转角。注意这里扭转角定义的是桨叶截面弦线和旋转平面的夹角，也就是桨距角。几何扭转角加上来流迎角的变化，共同决定当地实际攻角α = θ - φ，其中φ是入流角，也就是合速度方向和旋转平面的夹角。合速度的轴向分量是V + v_i（来流速度加诱导速度），切向分量是Ωr - v_t（旋转速度减切向诱导速度，注意旋向）。从这个三维速度三角形可以看到，如果诱导速度算错了，攻角就错，升阻力就错，性能自然全错，所以整个BEMT的成败都压在诱导速度迭代上。

翼型数据方面，典型的小型螺旋桨用的是Clark-Y、Eppler E63或者NACA系列薄翼型。Eppler等专业软件可以提供不同雷诺数下的Cl、Cd随攻角变化曲线。如果手上只有一份特定雷诺数的数据，我会用XFOIL在几个相关雷诺数下各算一遍，把数据表准备好。需要注意极曲线要覆盖负攻角到失速攻角以上的范围，因为螺旋桨内侧截面在低前进比时可能遇到接近甚至超过失速攻角的情况，程序里要做外插，否则查表会越界报错。

2.2 核心迭代公式与收敛判据

BEMT的标准解法是把动量理论结果和叶素结果联立求解两个未知数：轴向诱导因子a（轴向诱导速度除以自由来流速度），切向诱导因子a'（切向诱导速度除以桨尖旋转速度）。经典的Glauert模型给出如下方程组（为了在代码里清晰展示，我写成Matlab可读的形式）：

% 在每个径向站r处，当地实度比 sigma = B*c/(2*pi*r) % 入流角 phi = atan( (V*(1+a)) / (Omega*r*(1-a')) ) % 当地攻角 alpha = theta - phi % 升力系数Cl、阻力系数Cd 由alpha查翼型表获得 % 轴向动量方程（含Glauert大诱导速度修正）： % dT = 4*pi*r*rho*V^2*a*(1+a)*F*dr （小诱导速度状态） % 叶素方程：dT = 0.5*rho*Vrel^2*c*dr*(Cl*cos(phi) - Cd*sin(phi)) % 切向动量方程： % dQ = 4*pi*r*rho*V*Omega*r^2*a'*(1+a)*F*dr % 叶素切向力：dQ = 0.5*rho*Vrel^2*c*dr*r*(Cl*sin(phi) + Cd*cos(phi))

其中Vrel是当地合速度，F是叶尖和毂部修正因子（Prandtl修正），rho是空气密度。符号上要注意Cd项贡献的是阻力，它在轴向方程里会减小拉力，在切向方程里会增加扭矩，所以不能像很多简化教程那样只保留Cl项，否则在小前进比大攻角工况下误差特别大。迭代方法用最直接的松弛迭代：先假设a和a'全为0，算出攻角，再根据叶素公式算出dT和dQ，代入动量方程反解出a和a'，把它和上一轮的值做加权平均更新，然后重复，直到两次迭代的a和a'变化小于1e-6。松弛因子我一般取0.3到0.5，取太大容易振荡，取太小收敛慢，实际代码里可以用自适应松弛加快收敛。

2.3 Prandtl修正和Glauert修正的工程意义

如果不加修正，BEMT算出来的拉力在叶尖附近会偏高，因为动量理论假设桨盘上诱导速度均匀，但真实的螺旋桨叶尖会有涡脱落，叶片载荷迅速降为零。Prandtl叶尖损失因子F_tip是一个解析近似，它在接近叶尖时快速趋向于0，有效削弱叶尖段的载荷贡献。类似地，靠近毂部也有一个根部修正因子，虽然影响不如叶尖大，但在低前进比大推力的工况下，靠根部的低効率段会产生明显的扭矩，若不修正会高估效率。两个因子相乘得到总修正F。公式写出来是：

% 叶尖修正 f_tip = (B/2) * ((R-r) ./ (r * sin(phi_rad))); F_tip = (2/pi) * acos(exp(-f_tip)); % 根部修正（把R换成hub半径r_hub） f_hub = (B/2) * ((r-r_hub) ./ (r_hub * sin(phi_rad))); F_hub = (2/pi) * acos(exp(-f_hub)); F = F_tip .* F_hub; % 注意对exp参数做限幅，避免acos取值范围越界

另一个必须处理的是Glauert修正：经典的动量理论在诱导速度因子a超过0.5左右时会失效，因为此时滑流速度太大，理想动量假设已经进入“湍流风车”甚至涡环状态。对小螺旋桨在静止或低前进比大拉力工况，轴向诱导因子a很容易超过0.5，如果不修正，动量方程求出的a会偏大甚至发散。工程上常用的Clauert经验修正是把动量方程里的dT改写为包含一个由a主导的修正项的形式。我在代码里实际用的一种稳健做法是：对轴向诱导因子a，优先用“叶素求出的拉力反推动量方程a”的公式，若a > 0.5，就改用经验修正曲线的渐近线。这个修正对悬停和低前进比下的计算结果影响非常明显，建议无论如何都要加。

2.4 性能系数的无因次化与效率定义

迭代收敛后，把每个径向站算出的拉力微元dT和扭矩微元dQ沿径向积分，得到总拉力T和总扭矩Q。然后按标准的螺旋桨无因次系数定义换算：

• 拉力系数 C_T = T / (ρ·n²·D⁴)
• 功率系数 C_P = P / (ρ·n³·D⁵) = 2π·n·Q / (ρ·n³·D⁵)
• 前进比 J = V / (n·D)
• 效率 η = J · C_T / C_P

需要留意这里转速n用的单位是转每秒，如果代码里输入的是转每分，记得先除以60。另外，效率的定义本质上是“有用功率”（T·V）除以“轴功率”（2πnQ），在静止（V=0）时效率分子为零，这是正常的，那对应的是悬停状态，只看C_T和C_P就够了，不要试图在一个图里把悬停效率也画成无限大/零，会误导读图的人。

3. Matlab代码实现：模块划分与完整可运行骨架

3.1 几何数据与翼型数据表该怎么组织

我建议用函数化的思路组织代码，而不是把所有步骤都堆在一个for循环里。第一步是几何输入。假设我们算一个典型的10英寸（D≈0.254m）小型多旋翼桨，桨叶数B=2，半径R=0.127m，毂半径r_hub=0.012m。弦长和扭转角沿径向的分布可以用多项式拟合实测点获得，比如弦长从根部到叶尖从0.02m线性降到0.008m，扭转角根部35度、叶尖10度左右按线性递减，就是一个很粗糙但可用的输入。翼型数据表我用一个N行3列的矩阵存储攻角、升力系数、阻力系数，查表用interp1，线性插值足够。要注意攻角向量必须单调递增，这是interp1的硬性要求。

% 几何参数 B = 2; % 桨叶数 R = 0.127; % 叶尖半径 m r_hub = 0.012; % 毂半径 m rho = 1.225; % 密度 kg/m3 n_rev = 6000/60; % 恒定转速 6000 RPM 转每秒 Omega = 2*pi*n_rev; % 角速度 rad/s % 径向离散：40个站 nr = 40; r = linspace(r_hub, R, nr)'; c = interp1([r_hub R], [0.020 0.008], r, 'linear')'; % 弦长分布 theta_deg = interp1([r_hub R], [35 10], r, 'linear')'; % 几何扭转角分布 theta = deg2rad(theta_deg); % 翼型表：alpha_deg, Cl, Cd airfoil = load('clarky_polars.txt'); alpha_tab = deg2rad(airfoil(:,1)); Cl_tab = airfoil(:,2); Cd_tab = airfoil(:,3);

这里提前用一个线性分布的弦长和扭转角只是为了演示代码结构，真正的项目里应当用实测或设计图纸的离散点。实测点往往在不同径向位置的间距不均匀，处理时可以用Matlab的interp1把原始测量点重采样到等间距径向网格上，这样后续积分可以用更简单的梯形法。

3.2 核心求解循环的Matlab实现

螺旋桨性能求解的主体可以封装在一个函数里：输入是几何分布、翼型表、转速、来流速度，输出是每个径向站的a、a'、dT、dQ、总拉力、总扭矩、系数和效率。函数返回结构体是我比较推荐的方式，方便批量扫描前进比时统一收集结果。

function perf = bem_solve(r, c, theta, alpha_tab, Cl_tab, Cd_tab, ...) % 参数初始化 a = zeros(size(r)); ap = zeros(size(r)); % 迭代主循环 for iter = 1:200 a_old = a; ap_old = ap; for i = 1:length(r) % 当前半径处实度比、合速度、入流角 sigma = B*c(i) / (2*pi*r(i)); phi = atan( (V*(1+a(i))) / (Omega*r(i)*(1-ap(i))) ); alpha = theta(i) - phi; % 查表获取Cl、Cd if alpha > max(alpha_tab) alpha = max(alpha_tab); % 限制，避免插值越界 elseif alpha < min(alpha_tab) alpha = min(alpha_tab); end Cl = interp1(alpha_tab, Cl_tab, alpha, 'linear', 'extrap'); Cd = interp1(alpha_tab, Cd_tab, alpha, 'linear', 'extrap'); Vrel = sqrt((V*(1+a(i)))^2 + (Omega*r(i)*(1-ap(i)))^2); dT = 0.5*rho*Vrel^2*c(i)*(Cl*cos(phi) - Cd*sin(phi)); dQ = 0.5*rho*Vrel^2*c(i)*r(i)*(Cl*sin(phi) + Cd*cos(phi)); % 动量方程反推 a, ap（含Prandtl修正） F = prandtl_loss(B, R, r(i), phi); lhs_a = dT / (4*pi*r(i)*rho*V^2*F); a(i) = 0.5 * ( -1 + sqrt(1 + 4*lhs_a) ); % 注意当V接近0时该式退化，需要特殊处理，见5.3节 ... end if max(abs(a-a_old)) < 1e-6 && max(abs(ap-ap_old)) < 1e-6 break; end a = 0.4*a + 0.6*a_old; % 松弛 ap = 0.4*ap + 0.6*ap_old; end % 积分，输出系数 end

上面这段代码做了相当大的简化，重点在于展示骨架。实际细节里至少还要补三件事：一是攻角限制后要加警告或记录，否则可能掩盖叶素数据不足的问题；二是当V接近零时，轴向动量方程里含有V²的项会整体趋于零，直接导致除法爆掉，需要切换成悬停专用的退化方程；三是对力矩的动量方程要同样做切向因子a'的求解，它的表达式是另一个代数方程，采用和a类似的松弛迭代即可。这些我在第5节里会展开讲。

3.3 低前进比工况的数值退化与特殊处理

前进比J趋近于0，代表螺旋桨在原地不动高负荷运转（比如无人机垂直起降阶段）。这时候V=0，动量方程里的V(a)项变成0，用上面代码会有0/0风险。我会单独写一个分支：当V < 1e-6时，用直接求解诱导速度的退化公式。此时动量理论给出的理想功率关系是dT = 4πrρ·(rΩ)²·a'·(1+a)·F·dr，注意轴向诱导因子a的动量方程变成关于a的隐式方程，我没有用通用数值求根，而是直接用叶素的dT反推一个“等效诱导速度”，再用迭代强制收敛到同一个a，这种工程处理虽然缺少严格的数学推导，但在工程实践中表现稳定，速度也快。再一个办法是给V设一个下限比如0.01 m/s，这个值远小于正常工作速度，又能保证公式不炸，扫描曲线时直接把最低前进比设为0.005左右即可。两种方法我都试过，推荐后者配合代码分支：代码里保留一个if V > threshold的分支走标准动量方程，否则走悬停退化分支，两条路径算出的C_T在V=0.02 m/s处能够平滑衔接，曲线不跳变，说明处理合理。

3.4 批量扫描前进比并绘制性能曲线

恒定转速、不同前进比的扫描，本质就是在外层套一个for循环。转速固定在比如6000RPM，前进比从0.05到预期的最大值（通常到达效率峰值后继续增加到J小于1的区域），对应每个J的来流速度V = J·n·D。计算完所有J后，把C_T、C_P、η汇总成向量，一次plot出来。

J_list = linspace(0.05, 0.9, 40); CT = zeros(size(J_list)); CP = zeros(size(J_list)); ETA = zeros(size(J_list)); for k = 1:length(J_list) V = J_list(k) * n_rev * (2*R); perf = bem_solve(..., V, ...); CT(k) = perf.CT; CP(k) = perf.CP; ETA(k) = perf.eta; end % 绘制效率曲线 figure; yyaxis left; plot(J_list, CT, '-o'); ylabel('C_T'); yyaxis right; plot(J_list, ETA, '-s'); ylabel('\eta'); xlabel('J = V/(nD)'); grid on; axis tight;

真正做出来的效率曲线形状是一个先上升后下降的拱形，最佳效率点对应的J一般在0.3到0.6之间，这和我的实际测试经验一致——一个10英寸航拍桨在6000RPM下，效率最高对应的飞行速度约在8到12m/s，再快效率反而因为桨尖马赫数和失速/阻力增加而掉下来。这些物理特征如果曲线里没有出现，那多半是代码哪里出错了，而不是螺旋桨特殊。

4. 结果分析：怎么判断计算结果合理还是离谱

4.1 三个关键曲线的物理趋势与量级校验

先说量级。一个两叶10英寸桨在6000RPM下悬停，拉力一般在0.35到0.5公斤力左右，也就是3.5到5牛顿。C_T的量级大约是0.08到0.12。如果算出来是0.5或者0.001，基本可以断定某处有bug。功率呢，这种尺寸的桨悬停功率大约在50到100瓦，对应C_P量级约在0.04到0.06。这些经验值对常规的小型无人机桨有相当强的参考价值，可以作为第一次跑通代码的“冒烟测试”标准。

效率曲线的形状同样有严格的物理约束：J=0时效率是0，随着J增大，效率先上升，在某个J处达到峰值，之后下降。峰值效率对于设计良好的螺旋桨应该在0.7到0.85之间。如果你看到效率在低J时就超过了1，那必然是能量不守恒——通常在积分时搞错了单位，或者转速n误用了转每分钟而没有换算到转每秒，导致C_P计算偏小。我看到过初学者把RPM直接当Hz代入，结果效率普遍翻了几十倍，这个坑太典型了。

4.2 径向载荷分布：除了总性能，分布也值得看

只画总性能曲线还不够验算。我会额外画一个“沿径向的拉力密度分布dT/dr”和“局部攻角分布α(r)”，这个图的信息量非常大。正确的趋势应该是：低前进比（比如J=0.05）时，桨叶内侧攻角大，外侧攻角适中；沿径向向外，攻角逐渐减小。如果算出来某个中间位置攻角超过了翼型失速攻角（比如12度以上）并且大面积存在，说明该工况下桨叶已经部分失速，效率会明显下降，这和实验观察一致。另一个特征是叶尖处的dT/dr应该衰减到接近0，因为Prandtl修正压掉了载荷。如果叶尖载荷不但不衰减反而翘起来，就是Prandtl因子或入流角φ计算出了错，最常见的是角度制/弧度制混用导致cos/sin的计算结果南辕北辙。

4.3 收敛性分析与网格敏感性

BEMT这类数值方法有一个好习惯：先做收敛性检查再谈结果。我会做两个层面的检查。第一是迭代收敛，程序里记录每次迭代的均方根误差，理想情况是误差平滑下降，在几十步内跌到1e-6以下。如果误差震荡不收敛，多半是松弛因子太大或者攻角限制导致查表值跳变。第二是网格收敛，把径向站数从20逐步加密到80，看总拉力变化量。正常的BEMT在40个站以上时结果变化小于1%，如果20到80站的曲线有明显变化，说明弦长或扭转角分布有剧烈梯度，加密网格是必须的。我在自己的代码里默认用60个站，兼顾速度和精度。

4.4 拿实验数据或公开桨效测试来“对答案”

做工程最忌讳闭门造车。我强烈建议第一次跑通代码后，用一组公开的螺旋桨实验数据来校准。比如UIUC的螺旋桨数据库里有很多桨的BEMT对比数据，Propeller Database上也有一堆RTF厂商公布的风洞测试结果。我自己当时挑了一个Master Airscrew的MR 10x4.5桨，找它的官方拉力测试曲线，把几何量简化为矩形弦长分布和恒定扭转角，跑出来的C_T曲线在中等前进比段偏差在8%以内，效率峰值位置也基本对上。这说明BEMT在输入几何比较粗糙时的鲁棒性其实比很多教程里说的要好。当然如果偏差始终稳定偏大或偏小，可以检查翼型极曲线是否准确——很多公开桨用的翼型并非教科书上那几条标准的Clark-Y，而是厂家自己修型的，这里会产生系统偏差。

5. 调试实录：我踩过的坑和快速排查方法

5.1 角度制与弧度制的混用是头号杀手

再强调一遍，这个错误我犯过，找我请教的人也犯过。Matlab的三角函数sin、cos、tan默认输入是弧度，而气动数据表里攻角习惯用度。我的处理是：所有查表前把度数转成弧度，所有输入几何扭转角彻底转成弧度后再参与运算，只在输出画图的时候转回度数。程序里多处用到角度时，命名上明确区分theta_deg和theta_rad，不要在同一个变量上来回覆盖。另外，atan2和atan的选择也值得注意，BEMT里入流角φ = atan(Va / Vt)永远不会为负（因为Va、Vt都是正值），用atan就够了，但如果代码以后要扩展到风车状态，来流可能反向，建议直接用atan2，提前防一手。

5.2 查表超出范围与翼型数据外插的隐性错误

螺旋桨内侧叶素在低前进比时的攻角经常超过翼型极曲线的数据范围。我最初的做法是简单把攻角限制在表的最大值，结果内侧叶素的Cl被强制封顶，算出来的总推力在悬停状态明显偏低。后来我改成两端外插：对攻角大于最大数据点的情况，按最后两点的斜率继续外推，但要设一个上限比如25度，超过上限就用失速后的平板理论近似Cl = 2·sin(α)·cos(α)，Cd则按平板阻力增大。这个处理虽然粗糙，但至少在物理上比“封顶”要合理得多。另外注意翼型表的攻角范围至少要覆盖-5度到20度，否则大多数工况都会被截断，那一版代码基本不能用来扫描。

5.3 悬停点（V=0）迭代发散

低前进比扫描时，V设到0.01以下，轴向动量方程反推出来的a会突然异常大，迭代发散。这个现象的根源是方程结构本身在V=0时退化了，不是代码“哪里写错”。我后来用的解决方案是：把J的下限设置成0.02~0.03，对应的V在0.5m/s以上，既回避了奇异点，又对工程曲线没影响。如果真的要算精确悬停，走专门的悬停BEMT版本，输入就是转速和拉力需求，输出需要的桨距角——那是另一个话题了。

5.4 Prandtl修正里acos越界

acos的参数要求必须在[-1,1]之间，但在迭代初期a、a'还没收敛时，f_tip可能取得很大的正值或者很小的负值，exp(-f_tip)有时候会非常接近1，导致acos的参数由于浮点误差超过1。我在代码里加了明确的限幅：

arg = exp(-f_tip); arg(arg<0) = 0; arg(arg>1) = 1; F_tip = (2/pi) * acos(arg);

加了这个保护之后，程序稳健了很多，再也没出现过NaN。

5.5 常见问题速查表

5.6 提高代码稳健性和可读性的几个习惯

最后提几个我后来一直沿用的习惯。所有算法代码写成一个函数，输入输出用结构体，主脚本只负责定义参数和调用，这样批量扫描不同几何或转速时不用复制代码。在函数里写断言（assert），比如转速必须为正、弦长必须大于零、翼型表攻角必须单调，这些断言在参数写错时能立刻报错，而不是让程序带着错误数据跑出个精美但荒唐的曲线。代码里注释不写“这是什么”而写“为什么这样算”，比如Glauert修正旁边注释一句“避免低前进比动量方程发散”，三个月后再看代码能秒懂当时的决策逻辑，这对一个长期演进的项目尤其重要。

回到开头说的这个项目本身，把给定螺旋桨几何在恒定转速下的前进比扫描跑通之后，这套代码的复用价值远比一次作业大得多。我后来基于这套代码做了好几个变体：把Clark-Y极曲线换成NACA系列来对比不同翼型的效果，把矩形弦长分布换成实测桨叶离散点去预测大疆某款桨的悬停电流，甚至把输出端的性能系数连接到电机效率模型上，估算整机续航时间。BEMT的潜力不在公式本身多高深，而在于它把复杂的螺旋桨气动问题拆解成了几个可以用数值方法稳健求解的方程，而Matlab恰好是把这个过程表述出来最顺手的工具。自己在跑代码的过程中养成的“先想物理、再看公式、最后动手写”的习惯，才是这个项目留给我的更长期的收获。