基于BEMT的螺旋桨性能分析：Matlab实现不同前进比下推力扭矩效率计算

做螺旋桨性能分析的时候，我习惯拿出来当第一板斧的工具就是叶片单元动量理论。这个理论骨架不算新，上世纪二十年代末由Glauert那批人逐步奠定，但直到今天，无论是无人机螺旋桨预研、多旋翼动力匹配，还是轻型飞机定距桨选型，它依然是方案阶段百分之八十以上计算的地基。标题里这个组合很典型：给定几何、恒定转速、扫不同前进比。翻译成人话就是——手里已经有一副具体桨叶，转速锁定不变，让来流速度从零开始一点点加大，看这副桨的推力、扭矩和效率分别怎么变。这就是一套标准的螺旋桨稳态性能曲线生成流程，而Matlab是把它落地最顺手的工具。

这篇内容就是完整复现这套流程：从BEMT的两个祖传理论怎么融合讲起，把几何参数化、翼型数据查表、诱导因子迭代、推力扭矩积分、无量纲系数画曲线全部走一遍。适合的人很明确：正在做飞行器动力选型的学生、搞多旋翼或固定翼动力匹配的工程师、想用Matlab快速评估一副新桨几何的爱好者。你不需要提前会CFD，也不需要懂太深的空气动力学，只要能把公式落到循环里，跑出来的曲线就能给你相当靠谱的参考。

1. 叶片单元动量理论到底在算什么

1.1 两个祖传理论的各自短板

动量理论把螺旋桨想象成一个均匀吸气的致动盘，站在桨盘前后做动量守恒和能量守恒。它给的东西很干净：诱导速度、推力、功率之间有明确的解析关系，静推力状态能直接算出来，而且物理图像清晰——桨把气流加速了，加速的代价就是推力反作用在桨上。但它有个致命短板：完全不看桨叶的真实几何。你给它两副桨，一副宽叶大扭转，一副窄叶小平桨距，只要桨盘面积一样，动量理论给出的结果一模一样。这在工程里显然没法接受，因为换一副桨就换了一套性能。

叶素理论则相反，它把桨叶沿展向切成一片片二维翼型，每个截面根据当地合速度和攻角查出升阻力系数，再沿径向积分得到总推力和总扭矩。这个理论能看到弦长分布、扭转分布，也能带进翼型的升阻特性，但它默认来流就是均匀自由流，没有考虑一个关键事实：桨叶本身在扰动流场，每一片叶素感受到的有效来流已经被自己和其他叶片改变了。忽略诱导速度，攻角就会算错，推力可能差出几倍甚至一个量级。

两个理论单独拿出来都有硬伤：动量理论知道诱导流动但看不到几何，叶素理论看得到几何却猜不到诱导。叶片单元动量理论（BEMT）做的事情，就是强行把两者在每一段微小环带上接起来：动量理论说这个环带上应该有多少推力和扭矩，叶素理论也用当地气动力给出一份答案，两者不一致就调整诱导速度，调到两边相等为止。这就是整个BEMT的魂。

1.2 把动量理论和叶素焊在一起的推导

具体推导是这么走的。把桨盘沿径向切成很多细环带，半径 ( r ) 处取宽度 ( dr ) 的微元环带，面积是 ( 2\pi r dr )。设远处来流速度为 ( V_0 )，盘面处轴向诱导速度为 ( aV_0 )，于是盘面处实际轴向速度为 ( V_0(1+a) )。再定义 ( a' ) 为切向诱导因子，盘面处气流获得的切向速度是 ( a'\Omega r )，这里的 ( \Omega ) 是螺旋桨角速度。

动量理论先给出推力和扭矩的微分：

[ dT = 4 \pi r \rho V_0^2 a(1+a) , dr ]

[ dQ = 4 \pi r^3 \rho V_0 (1+a) a' \Omega , dr ]

这两个公式的物理来源不复杂：推力等于单位时间流过环带的气流动量增量，扭矩等于角动量增量。关键是它们都跟诱导因子挂钩，而诱导因子正是我们要解的未知数。

再看叶素一侧。在半径 ( r ) 处，叶素看到的合速度由两部分组成：轴向分量 ( V_0(1+a) )，切向分量 ( \Omega r(1-a') )。合速度大小和入流角满足：

[ V_{rel} = \sqrt{ [V_0(1+a)]^2 + [\Omega r(1-a')]^2 } ]

[ \varphi = \arctan \frac{V_0(1+a)}{\Omega r(1-a')} ]

当地几何桨距角记为 ( \theta )，攻角就是：

[ \alpha = \theta - \varphi ]

接下来查翼型升阻力系数 ( C_l(\alpha) )、( C_d(\alpha) )，然后沿来流方向和垂直来流方向分解，再投影到螺旋桨轴向和切向。这里我习惯把两个投影系数直接算好：

[ C_y = C_l \cos\varphi - C_d \sin\varphi ]

[ C_x = C_l \sin\varphi + C_d \cos\varphi ]

每个叶素贡献的推力和扭矩为：

[ dT = B \cdot \frac{1}{2} \rho V_{rel}^2 c , C_y , dr ]

[ dQ = B \cdot \frac{1}{2} \rho V_{rel}^2 c , C_x , r , dr ]

其中 ( B ) 是桨叶数，( c ) 是当地弦长。

现在关键一步：令动量理论的 ( dT ) 等于叶素理论的 ( dT )，再利用 ( V_{rel}\sin\varphi = V_0(1+a) ) 消去 ( V_{rel} )，整理后得到轴向诱导因子的更新关系：

[ \frac{a}{1+a} = \frac{B c C_y}{8 \pi r F \sin^2\varphi} ]

同理，联立扭矩方程，用 ( V_{rel}\cos\varphi = \Omega r(1-a') ) 消元，得到：

[ \frac{a'}{(1-a')^2} = \frac{B \Omega c C_x}{8 \pi V_0 (1+a) F \cos^2\varphi} ]

公式里的 ( F ) 是叶尖损失修正因子，后面单独讲。这个推导过程看起来全是代数，但它在编程里对应一个非常清晰的动作：每一轮迭代，用当前的 ( a )、( a' ) 算出攻角，查表得到 ( C_l )、( C_d )，再用上面两个式子更新 ( a )、( a' )，反复循环直到收敛。整个Matlab核心求解器，说白了就是上面这两行公式的循环。

1.3 为什么这个标题适合用BEMT做

有经验的工程师都知道，螺旋桨性能分析有很多条路：CFD当然最接近真实，但网格生成、湍流模型、旋转域设置，一套算下来动辄几小时甚至几天，做参数扫描根本扛不住；风洞实验精度高，但成本和时间都不是方案前期能承受的。BEMT计算量小到可以忽略不计，Matlab里一个循环几毫秒就能算完一个工况，扫几十个前进比也只是眨眼的事。同时它又能保留几何细节——弦长分布、扭转分布、翼型升阻特性全部能带进去，这正好匹配“给定螺旋桨几何形状”这个前提条件。

当然BEMT有自己的边界：它假设每个环带之间流动互不干扰，忽略径向流动，在强三维效应、大失速、低前进比湍流尾迹区精度会明显下降。但作为分析不同前进比下性能趋势的第一轮工具，它几乎是最优选择。我们后面所有代码和结论，都是在这个边界内讨论的。

2. 几何、翼型和工况：把螺旋桨翻译成数据

2.1 弦长分布与扭转分布的建模

让程序跑起来的第一步，是把几何变成数组。一副螺旋桨的几何，在BEMT里归根结底就是两件事：沿展向的弦长分布 ( c(r) )，以及沿展向的几何桨距角分布 ( \theta(r) )。这两个分布一般来自桨叶设计图纸、三坐标扫描数据，或者直接用现成桨叶实测。没有实测数据的时候，用简化模型做演示也很常见。

我在代码里用了一个典型小型双叶螺旋桨的参数：半径 ( R = 0.1,\text{m} )，桨叶数 ( B = 2 )，展向取41个站点。弦长按一个单调变化规律分布，内段宽外段窄，这符合大多数小桨的特点；扭转则从根部30度左右线性减到桨尖20度左右。下面这段就是几何初始化的Matlab写法：

% 几何与工况基础参数 rho = 1.225; % 空气密度 kg/m^3 R = 0.1; % 桨叶半径 m B = 2; % 桨叶数 n = 6000 / 60; % 恒定转速 rps（6000 RPM） omega = 2 * pi * n; % 角速度 rad/s D = 2 * R; % 桨盘直径 m % 展向离散：根切起点放0.03R，避免桨毂区域奇点 r = linspace(0.03*R, R, 41)'; % 示例弦长分布：根部略宽，向桨尖收窄 c = 0.02 * (0.9 - 0.5 * (r / R).^2); % 示例扭转分布：线性扭转，根30度，尖20度 theta = deg2rad(30) - deg2rad(10) * (r - 0.03*R) / (R - 0.03*R);

这里有一个细节：站点规划尽量在根部做一个截断，不要从 ( r=0 ) 开始。因为BEMT在接近桨毂时，叶素环带周长趋近于零，模型本身不再成立，而且桨毂处通常有整流罩或安装结构，气动上也不是有效区域。一般取 ( r_{\min} = 0.1R \sim 0.15R ) 起步，演示代码取0.03R是为了让小桨看起来更连续，实际工程建议放宽到0.1R以上。

2.2 翼型气动数据的准备与插值

几何有了，下一步是翼型。每个叶素截面上的翼型，需要用升力系数 ( C_l ) 和阻力系数 ( C_d ) 随攻角 ( \alpha ) 的变化曲线。这些数据哪里来？三个途径：翼型风洞实验数据、XFOIL/CFD计算结果、设计手册经验数据。小螺旋桨常用的是RAF-6、Clark-Y、E63这类低速翼型，我这里演示就用一组RAF-6简化数据。

典型的翼型数据表长这样：

真实情况下数据点会更密，失速后区间的数据尤其重要。把这张表存成两个列向量，之后在迭代里用interp1做线性插值：

alpha_tab = [-4 0 4 8 12 16]; % 攻角表 deg Cl_tab = [-0.32 0.31 0.64 0.93 1.18 1.10]; Cd_tab = [0.012 0.008 0.010 0.016 0.030 0.080]; % 在BEMT迭代中，alpha_deg是当前叶素攻角 Cl = interp1(alpha_tab, Cl_tab, alpha_deg, 'linear', 'extrap'); Cd = interp1(alpha_tab, Cd_tab, alpha_deg, 'linear', 'extrap');

这里我特别提醒一下外推的问题。interp1的'extrap'选项会把数据表范围外的攻角强行线性外推，看起来能用，但物理上一团糟。失速后升力不会无限增长，阻力也不会一直线性变大。所以在正式计算里，我一般会手动限制攻角范围，超出数据表的部分做截断处理，或者用专门考虑大攻角行为的气动模型。后面第五章会详细展开这个坑。

2.3 前进比、转速和来流速度的关系

所有几何都准备好了，最后定义工况。标题里“不同前进比下恒定转速”这个设定，对应的数学关系非常直接。前进比定义为：

[ J = \frac{V}{nD} ]

其中 ( V ) 是来流速度，( n ) 是转速（转/秒，rps），( D ) 是螺旋桨直径。它描述的是“前方的气流在桨转一圈期间向前走了多少倍直径”。这个无量纲量决定了螺旋桨的工作状态。恒定转速意味着 ( n ) 不动，扫不同前进比本质上就是扫不同来流速度：

[ V = J \cdot n \cdot D ]

拿刚才那个0.2m直径、6000RPM的桨来说，( n = 100,\text{rps} )，( D = 0.2,\text{m} )，所以每个前进比对应的来流速度是：

从悬停到高速巡航，这一条J轴覆盖了螺旋桨几乎全部典型工作状态。在恒转速条件下做这个扫描，得到的 ( C_T(J) )、( C_P(J) )、( \eta(J) ) 就是一副桨的动力特性名片。

3. Matlab代码实现：从求解器到性能曲线

3.1 程序骨架怎么搭

整个Matlab实现分四个层次。最外层是主脚本，负责定义几何、工况、调用求解器、画图；第二层是单工况求解函数，输入某个 ( J ) 对应的来流速度，输出总推力和总扭矩；第三层是核心BEMT迭代函数，处理所有叶素的诱导因子求解；最底层是翼型数据查询函数。这样分层的好处很直白：以后换一副桨只需要改几何数组，换翼型只改数据表，换工况只需要改循环变量，完全不需要动核心迭代逻辑。

主脚本的结构大概这样：

% 主脚本：恒定转速下扫描前进比 J_array = 0:0.05:1.0; CT = zeros(size(J_array)); CP = zeros(size(J_array)); eta = zeros(size(J_array)); for k = 1:length(J_array) V = J_array(k) * n * D; [T, Q, P] = BEMT_solver(prop, airdata, V, omega, rho); CT(k) = T / (rho * n^2 * D^4); CP(k) = P / (rho * n^3 * D^5); eta(k)= J_array(k) * CT(k) / CP(k); end

这个循环里的BEMT_solver就是单工况求解函数。这里先提一个重要概念：无量纲系数。为什么用 ( \rho n^2 D^4 ) 和 ( \rho n^3 D^5 ) 来做分母？因为在恒定转速下，推力和功率随空气密度、转速、直径的变化大致服从这些比例关系，除以它们之后，( C_T )、( C_P ) 就完全由螺旋桨几何和前进比决定，可以跨尺寸、跨转速比较。这是工程上处理螺旋桨数据最标准的做法。

3.2 核心BEMT迭代求解代码

单工况求解函数BEMT_solver是整个程序的心脏。输入是几何数组、翼型数据、来流速度和角速度，输出是总推力、总扭矩和总功率。完整代码如下：

function [T, Q, P] = BEMT_solver(prop, airdata, V, omega, rho) r = prop.r; c = prop.c; th = prop.theta; R = prop.R; B = prop.B; N = length(r); a = zeros(N, 1); % 轴向诱导因子 ap = zeros(N, 1); % 切向诱导因子 % 阻尼因子：新值占比，经验取值0.5~0.8 damp = 0.6; maxIter = 300; tol = 1e-5; for it = 1:maxIter a_old = a; ap_old = ap; for i = 1:N % 入流角 phi = atan( V * (1 + a(i)) / (omega * r(i) * (1 - ap(i))) ); % 攻角 alpha_deg = rad2deg( th(i) - phi ); % 翼型查表 [Cl, Cd] = airfoil_lookup(airdata, alpha_deg); % 叶素投影系数 Cy = Cl * cos(phi) - Cd * sin(phi); Cx = Cl * sin(phi) + Cd * cos(phi); % 叶尖损失修正 f_tip = (B / 2) * (R - r(i)) / (r(i) * sin(phi)); F = 2 / pi * acos(exp(-f_tip)); % 局部实度 sigma = B * c(i) / (2 * pi * r(i)); % 轴向诱导更新 nomA = sigma * Cy / (8 * F * sin(phi)^2); a_new = nomA / (1 + nomA); % 切向诱导更新 nomB = sigma * omega * r(i) * Cx / (4 * V * (1 + a(i)) * F * cos(phi)^2); ap_new = nomB / (1 + 2 * nomB); % 近似形式，小扰动下用 % 阻尼更新，缓解振荡 a(i) = (1-damp) * a(i) + damp * a_new; ap(i) = (1-damp) * ap(i) + damp * ap_new; end err = max(abs(a - a_old)) + max(abs(ap - ap_old)); if err < tol break; end if it == maxIter warning('BEMT迭代达到上限，未完全收敛'); end end % ... 推力扭矩积分（下节） end

这段代码里有几个值得琢磨的细节。阻尼因子damp是我实际调试时加的，不加它，低前进比工况下诱导因子常常在两次迭代之间来回跳，加了阻尼之后系统明显稳定，代价只是收敛步数多一点点。airfoil_lookup是一个独立的查表子函数，避免主循环里塞太多插值代码。

关于切向诱导更新式的写法，我再解释一下。前面理论推导给出的严格形式是：

[ \frac{a'}{(1-a')^2} = \frac{B \Omega c C_x}{8 \pi V_0 (1+a) F \cos^2\varphi} ]

直接解这个二次方程，得到：

[ a' = \frac{2m + 1 - \sqrt{4m + 1}}{2m} ]

其中 ( m = \frac{B \Omega c C_x}{8 \pi V_0 (1+a) F \cos^2\varphi} )。在代码里我用了小扰动近似 ( a'/(1-a') \approx a' )，得到更简单的ap_new = nomB / (1 + 2*nomB)。这个近似在 ( a' < 0.2 ) 时误差很小，但对于低前进比、大载荷工况，建议用严格的二次方程解，或者直接用牛顿迭代求根。严格版本虽然代码长几行，但边界更干净，不容易在重载状态下跑飞。

3.3 推力和功率积分与无量纲系数

迭代收敛之后，每个叶素都有了收敛的 ( a )、( a' )，也有了当地合速度、攻角、升阻力系数。最后一步是沿展向积分：

T = 0; Q = 0; for i = 1:N phi = atan( V * (1 + a(i)) / (omega * r(i) * (1 - ap(i))) ); alpha_deg = rad2deg( th(i) - phi ); [Cl, Cd] = airfoil_lookup(airdata, alpha_deg); Cx = Cl * sin(phi) + Cd * cos(phi); Cy = Cl * cos(phi) - Cd * sin(phi); Vrel = V * (1 + a(i)) / sin(phi); dT = B * 0.5 * rho * Vrel^2 * c(i) * Cy; dQ = B * 0.5 * rho * Vrel^2 * c(i) * Cx * r(i); T = T + dT * (R / N); % dr近似 Q = Q + dQ * (R / N); end P = Q * omega;

这个地方有个常见的积分细节：如果站点是均匀离散的，( dr ) 可以近似为R/N；如果像代码里用linspace(0.03*R, R, N)，实际每个站点间距不完全相等，更严谨的写法是直接用diff(r)构造积分权重。我倾向于直接保留径向数组然后用trapz(r, dT_array)做梯形积分，这样无论站点怎么分布都不会出问题。

主脚本里算无量纲系数时，功率用 ( \rho n^3 D^5 ) 做分母。有些教材会在分母里加0.5、甚至把 ( n ) 写成角速度，不同定义之间差一个常数倍。我的建议是：自己代码里固定用一套定义，和文献对比时先确认定义是否一致，不要让换算错误毁掉一天的调试成果。

3.4 曲线绘制与结果保存

画图本身不复杂，但布局有讲究。我一般用subplot三行一列或者一行三列，把三张曲线放在同一张图里，方便把横轴J对齐着看：

figure('Color','w'); subplot(3,1,1); plot(J_array, CT, 'b-', 'LineWidth', 1.6); ylabel('C_T'); grid on; title('恒定转速下螺旋桨性能随前进比变化'); subplot(3,1,2); plot(J_array, CP, 'r-', 'LineWidth', 1.6); ylabel('C_P'); grid on; subplot(3,1,3); plot(J_array, eta, 'k-', 'LineWidth', 1.6); xlabel('J'); ylabel('\eta'); grid on;

同时建议把每个展向站点在特定J下的攻角分布、诱导因子分布也画一下，因为只看总体性能曲线，你很难发现某个截面是否已经处于失速状态。比如我常画一张 ( \alpha(r) ) 分布图，能看到低速大拉力状态下内段攻角是否已经远超失速角，这对理解曲线形状突变非常有帮助。

4. 结果怎么看：曲线趋势与工程译读

4.1 效率峰值、静推力与低速状态的物理边界

代码跑通之后，首先检查曲线是否符合物理常识。推力系数 ( C_T ) 随 ( J ) 增大单调下降：J小的时候来流速度低，桨叶每个截面攻角大，拉力大；J增大了，来流把攻角压小，拉力自然降下来。到J足够大时，攻角可能变成负的，推力甚至会变成负值——这在螺旋桨制动或顺桨状态才会出现，正常的推进分析里我们只关心 ( C_T>0 ) 的区域。

效率 ( \eta ) 的典型形状是两头低中间高，在某个J处出现峰值。峰值对应的前进比就是这副桨的设计点附近。峰值左侧效率下降是因为载荷很大、诱导损失和叶型损失都快速增加；峰值右侧效率下降则是因为攻角变小，升阻比变差，桨叶在“吃”多余的阻力。如果算出来的效率峰值超过1，那几乎肯定有问题，要么是损失项被忽略得太厉害，要么是前面的动量方程推导有符号或系数错误——注意，在纯损失为零的理想情况下效率可以趋近1，但带入实际翼型阻力后不可能超过1。

J=0是一个特殊的边界点。那意味着 ( V=0 )，整个流动处于悬停状态。这时候动量理论里盘面诱导速度远比来流大，( a ) 值很容易超过0.5，进入所谓的湍流尾迹状态。理想动量方程在这个区域的推力和诱导速度关系已经失准，所以很多BEMT程序会在 ( a>0.4 ) 时切入Glauert经验修正。如果你的分析重点是悬停或极低速工况，这个修正不能省。

4.2 CT、CP、eta之外，还要关心哪些依变量

三条主曲线之外，我建议至少额外检查两个分布。第一个是每个J下沿展向的攻角分布 ( \alpha(r) )，它能直观看出一副桨是否在某个工况下内段失速或者外段接近失速。第二个是沿展向的载荷分布 ( dC_T/dr ) 或局部升力系数 ( C_l c / R ) 的分布，这个分布决定结构受力和噪声来源。BEMT虽然给不出精确的三维流场，但给这些工程判断已经完全够用。

举个例子，以前我算过一副偏航速优化的桨，扭转设计比较激进，根部攻角在小J时达到了16度以上，而翼型数据表里失速角只有14度。只看 ( C_T(J) ) 曲线，只是觉得低速段推力增长比预期的慢一点，但看 ( \alpha(r) ) 分布立刻就能定位到内段翼型已经失速，问题一目了然。这种诊断能力，是单看总体性能曲线得不到的。

5. 我实际跑数时踩过的坑

5.1 低前进比迭代发散的处理

这个坑几乎谁跑BEMT都会撞上。J=0或者J很小的时候，来流速度接近零，入流角 ( \varphi ) 非常小，( \sin\varphi ) 趋向于零，而诱导更新公式分母里有 ( \sin^2\varphi )，数值上极易出现爆炸性的巨大更新量。我早期版本在J=0时直接跑出NaN，最后一行行查才定位到这个问题。

解决办法可以组合使用。首先，给诱导因子加阻尼，新值不要直接替换旧值，按比例混合。我上面代码里的damp = 0.6就是这个作用，调试时曾经为了把J=0跑稳，直接把阻尼压到0.2，收敛慢但稳。其次，在迭代过程中给 ( a ) 设上限，比如 ( a_{\max}=0.7 )，超过就截断。物理上，大载荷状态动量方程本来就不够准，强行算更大的a只是安慰自己。第三，改善初始化，不要从 ( a=0 ) 开始。对于J=0工况，先给一个 ( a=0.3 ) 左右的初值，往往能让迭代直接落到正确的解附近，省掉头几十次震荡。

5.2 推力系数曲线抖动的排查

如果算出来 ( C_T(J) ) 不是光滑曲线，而是局部有小锯齿，问题大概率出在翼型数据插值或展向站点数量。我遇到过两种情况：一种是翼型表数据点太稀疏，相邻两个攻角的升力系数跳变过大，导致迭代时攻角微小的变化引起查表结果大幅波动；另一种是展向站点太少，只有20个点左右，叶尖区域的损失修正对网格特别敏感。把站点提到40~60个、把翼型表在关键区的数据加密到每2度一个点，抖动基本会消失。

还有一种隐蔽的情况：叶尖损失因子F在接近叶尖时 ( f_tip ) 趋近无穷大，按公式 ( F = \frac{2}{\pi}\arccos(e^{-f}) ) 算出来的F趋近0，很正常。但如果站点恰好取到 ( r=R ) 的叶尖点，exp(-f_tip)会下溢成0，F也是0，这没问题；麻烦的是如果某些实现里F做了分母，会出现除零错误。我的习惯是最后一站取到0.98R或0.99R，不给叶尖点留任何找茬的机会。

5.3 翼型数据外推的坑

前面提到过'extrap'的危险性。实际飞行中螺旋桨攻角范围比大多数翼型表宽得多，低速高载荷时内段攻角可以轻松超过20度，而XFOIL算出来的翼型数据通常只到15度左右就停了。如果强制线性外推，( C_l ) 会无限上涨，算出来的推力和功率大到离谱，而且整个曲线看起来“性能特别好”。我第一次做BEMT就栽在这里，当时还兴奋地以为设计了一副超高性能桨，后来把攻角分布一打印，发现内段攻角全是28度、30度，升力系数外推到了2.5，物理上完全不可能。

现在我的做法是：外推前先手动审查攻角范围，把超出数据表的部分按失速后数据截断处理。失速后的翼型行为大致是升力系数不再增长甚至缓慢下降，阻力系数显著上升。如果不追求绝对精度，最稳妥的做法就是把数据表外的 ( C_l ) 固定为失速点的值，( C_d ) 按每度0.015~0.02的速率线性增加。这样做出来的曲线虽然保守一点，但至少物理趋势是可信的。

最后再分享一个小习惯：每次跑完新的几何或新的翼型数据，我都会先把J=0.5的单点结果手动验算一遍，拿出计算器，挑一个中间站位的叶素，按公式手算攻角、升力系数和推力贡献，跟程序输出的那个站点数据对比。这个动作看起来笨，但能最快发现符号错误、数组错位、单位混乱这些BEMT里最常见的低级问题。程序跑得再快，也快不过你一次靠谱的手工校验。这个习惯我保留到了现在，也建议每个做螺旋桨性能分析的人都保持住。