做小型螺旋桨性能评估时我最常用的第一轮工具并不是CFD而是叶素动量理论BEMT。标题里那个“动量理论剖面刀片方法”其实就是这个理论的另一种直译国内工程圈更习惯叫叶素动量理论或者片条理论。它的任务很明确给定一套螺旋桨几何在转速恒定的前提下把来流速度从零逐步往上推算出不同前进比下的拉力、扭矩、功率和效率。这个工作用Matlab实现并不难难在诱导速度的迭代收敛处理以及把翼型气动数据、几何参数和动量约束接在一起不出错。我把实际跑通这套流程的经验整理在下面希望能帮你少踩几个坑。1. 动手之前先搞懂BEMT到底在算什么1.1 为什么第一轮评估用BEMT而不是CFD螺旋桨气动分析的方法其实不少从经验公式、动量理论、叶素动量理论到全三维RANS、Lattice Boltzmann甚至直接上风洞台架测力。放在同一个天平上比BEMT的定位是“方案阶段的快速筛选器”。我做预研时选它有三个理由。第一计算量小到可以忽略把一片桨叶沿展向切几十个叶素每个叶素解一个二维非线性方程Matlab跑一个工况就是毫秒级扫一遍前进比也就几秒钟。第二物理因果链非常清楚弦长变化、扭转变化、转速变化、来流变化都会直接映射到某个剖面的迎角、升阻力系数和整桨性能曲线上调试时一眼就能看出问题在哪。第三输入数据全是设计阶段已经确定的参数不需要网格、不需要湍流模型、不需要壁面处理。CFD的优势在细节比如分离流、桨尖涡、三维失速、机翼与滑流的强干扰这些BEMT天然给不了。但在“给定几何、恒定转速、扫前进比”这类任务里CFD的单工况耗时从小时到天不等几何稍有改动又要重新来一轮网格根本不划算。我的经验是先用BEMT把设计空间快速缩小选定少数几个候选几何再让CFD去做精确校核最后实验验证。这比一上来就全流程CFD高效得多。方法单工况耗时输入要求适用范围典型用途BEMT毫秒级弦长、扭角、翼型极曲线附着流为主、载荷中等选型、趋势分析、动力匹配CFD RANS小时到数天三维几何、网格、湍流模型分离流、复杂干扰、流场细节最终校核、诊断特殊现象台架试验天到周样机、台架、传感器真实完备定型验证、标定1.2 动量理论和叶素理论是怎么互补的“BEMT”这个名字里的两个部分不是在讲同一件事。动量理论把螺旋桨简化成一个无限薄的致动盘假设来流均匀、不可压、无粘并忽略尾流旋转的细节只研究桨盘前后的速度变化和压力突跳。它给出的是宏观结论桨盘从流场中吸收了动量气流通过桨盘后轴向速度增加同时获得与旋转方向相反的切向速度于是形成推力和反扭矩。动量理论的优点是不依赖桨叶长什么样缺点是它根本不认识桨叶几何。叶素理论则反过来。它把桨叶沿展向切成一个个二维剖面每个剖面看成一段独立翼型用当地相对速度、迎角和二维翼型升阻力数据算出剖面气动力再沿展向积分。叶素的优点是能够考虑真实的弦长、扭转和翼型缺点是它没有回答“来流进入桨盘后到底被加速了多少”。这个加速量恰恰决定了当地迎角。BEMT把两边接起来动量理论给出轴向和周向诱导速度与载荷的关系叶素理论给出桨叶几何产生的载荷关系两者必须相等于是诱导速度被唯一确定。这就是它能把几何参数变成性能数据的核心机制。标题里的“剖面刀片方法”强调的正是叶素那一侧的“切剖面”思想翻译上可以接受但习惯上还是叫叶素动量理论。1.3 输入、输出和适用边界输入部分非常明确。几何上需要螺旋桨半径、沿展向的弦长分布和安装角分布安装角实际上就是每个剖面的几何扭转角。气动上需要翼型极曲线也就是升力系数Cl和阻力系数Cd随迎角α的查表数据。工况上需要转速rpm、来流速度Vinf、来流密度ρ和叶片数N。输出部分首先是每个叶素的诱导因子a和a、当地迎角、受力然后沿展向积分得到整桨推力T、扭矩Q、功率P和推进效率η。为了方便跨工况对比我还会额外算无量纲系数推进系数C_T T/(ρ n²D⁴)、功率系数C_P P/(ρ n³D⁵)、前进比J Vinf/(nD)其中n是每秒转速D是直径。适用边界也必须要讲清楚。BEMT建立在均匀定常来流、附着流、低速亚声速这些假设上所以桨尖马赫数最好不要超过0.7翼型数据要覆盖实际遇到的雷诺数范围且不需要处理失速后的复杂分离。如果来流很不均匀比如螺旋桨装在机翼后方或者在大攻角状态下工作BEMT只能给出趋势级的结果。这个约束不是缺点而是工具本身的定位快速、可解释、覆盖大多数正常设计点。2. 把BEMT公式拆到能直接写代码的程度2.1 截面速度三角与微元气动力取半径r处的一个小段厚度为dr。当地轴向速度是来流速度加上诱导速度写成Vinf×(1a)这里的a叫轴向诱导因子。切向速度稍微绕一些桨叶以角速度Ω旋转当地旋转速度为Ωr但由于桨盘把动量给了气流气流获得与旋转方向相反的切向速度所以在叶素坐标系中切向来流速度表示为Ωr×(1-a)a是切向诱导因子。不同教材对a的符号定义会有差异但只要自己在代码里保持一致结果不会出问题。相对速度Vrel就是由轴向和切向两个分量合成的Vaxial Vinf * (1 a); Vtang Omega * r * (1 - ap); Vrel sqrt(Vaxial^2 Vtang^2); phi atan2(Vaxial, Vtang); alpha theta - phi;phi是流入角theta是几何安装角alpha就是当地迎角。有了alpha就能从翼型极曲线插值得到Cl和Cd。单个叶素的升力是0.5ρVrel²c dr×Cl阻力是0.5ρVrel²c dr×Cd。把升力阻力分别沿轴向和切向投影轴向分量形成拉力切向分量形成对桨毂的旋转阻力矩。乘上叶片数N就得到整个桨盘在这一小段上的气动贡献。这里有个地方特别容易写错升力方向定义为垂直于相对速度阻力方向定义为平行于相对速度。投影到轴向时需要用流入角phi的正弦余弦公式是dT_el N×0.5ρVrel²c dr(Cl cosφ − Cd sinφ)dQ_el N×0.5ρVrel²c dr(Cl sinφ Cd cosφ)×r。如果角度定义反了推力变成阻力整个迭代直接发散。2.2 动量方程那一侧的约束同一小段还可以从动量角度建立方程。单位时间流过半径r处环形微元的质量流量是ρ×2πr dr×Vinf(1a)。气流通过桨盘后轴向速度增量在下游逐渐达到2aVinf由动量定理轴向力的增量就是质量流量乘以速度增量整理后得到dT_mom 4πρr Vinf²(1a)a dr。这是BEMT中轴向诱导因子a的动量约束。切向方向同理。气流通过桨盘后获得旋转速度角动量变化产生扭矩。最常见的表达形式是dQ_mom 4πρr³ Ω Vinf(1a)a dr。这个式子里的倍率在不同教材中可能差一个2的因子因为对“盘面处的切向速度”和“尾流远场切向速度”的定义不同。我在代码里采用上面这组形式推导过程中只要保证质量流量和速度增量的定义一致最终结果就是自洽的。从物理上看动量方程说的是“流场给定诱导速度因此能容纳多大载荷”叶素方程说的是“桨叶几何在某个诱导速度下能产生多大载荷”。两者必须相等这正是迭代求解的出发点。2.3 为什么必须迭代而不是直接解公式如果把叶素方程和动量方程合起来看a和a出现在速度合成、流入角、迎角、Cl/Cd插值里而Cl/Cd通常是一张离散的查表和α的非线性函数。整个方程组是强非线性隐式方程不可能像线性方程组那样一步解出闭式解。实际操作中几乎都用数值迭代。最稳妥的做法是对每个叶素列两个残差方程F1(a,a) dT_el − dT_mom F2(a,a) dQ_el − dQ_mom然后用二维牛顿法或者Matlab的fsolve把a和a解出来。注意每个叶素的a和a只对这个半径位置有效所以整个桨叶在计算时是逐段独立求解的唯一的联系只是把力积分起来。这个“逐段独立”的好处是代码结构清晰坏处是一旦某一段解飞了整条曲线就会在对应位置出现一个尖刺。3. Matlab实现从几何数据到性能曲线3.1 输入数据怎么准备我习惯把几何参数放在一个表里每一行对应一个径向截面归一化半径rNorm、弦长c、安装角thetaDeg。翼型极曲线单独放两个文件或者一个矩阵第一列是alpha第二列Cl第三列Cd。Matlab里用readmatrix读入方便替换翼型。下面的示例描述一个半径0.254米、两片桨叶的几何采样30个截面R 0.254; % 桨半径单位米 r linspace(0.08*R, R, 30); c 0.05 * ones(size(r)) .* (1 - 0.5 * r / R); % 弦长从根到尖递减 thetaDeg 30 - 20 * r / R; % 安装角从根到尖递减 Nb 2; % 叶片数 rho 1.225; % 来流密度 rpm 6000; Omega rpm * 2 * pi / 60; D 2 * R; n_rps rpm / 60;实际的几何数据通常来自CAD或三维扫描弦长和扭角分布不一定是光滑函数但只要每个径向位置都有对应值BEMT完全可以直接处理。翼型极曲线需要注意迎角范围。低速小型螺旋桨常用的翼型在-5度到15度范围内能有完整的Cl/Cd表但如果扫前进比时当地迎角超出了这个范围插值函数会返回NaN迭代直接失败。这个问题我在第5节专门讲。3.2 核心求解器代码骨架核心求解函数接收半径数组、弦长数组、安装角数组和工况参数返回整桨推力、扭矩、诱导因子分布。残差函数就是上一节的F1和F2。这里给出一个能直接修改使用的骨架function [T, Q, P, a, ap] bemSolve(r, c, thetaDeg, Vinf, Omega, rho, Nb) n length(r); theta deg2rad(thetaDeg); a zeros(n,1); ap zeros(n,1); T 0; Q 0; x0 [0.05; 0.01]; % 第一个叶素用保守初值 for i 1:n if i 1 x0 [a(i-1); ap(i-1)]; % 用上一叶素的解作为初值 end opt optimoptions(fsolve, Display, off, ... MaxIterations, 600, ... FunctionTolerance, 1e-12, ... StepTolerance, 1e-12); x fsolve((x) bemResidual(x, r(i), R, c(i), theta(i), ... Vinf, Omega, rho, Nb), x0, opt); a(i) x(1); ap(i) x(2); va Vinf * (1 a(i)); vt Omega * r(i) * (1 - ap(i)); Vrel sqrt(va^2 vt^2); phi atan2(va, vt); alpha theta(i) - phi; [Cl, Cd] lookupAirfoil(alpha); dT Nb * 0.5 * rho * Vrel^2 * c(i) * ... (Cl * cos(phi) - Cd * sin(phi)); dQ Nb * 0.5 * rho * Vrel^2 * c(i) * ... (Cl * sin(phi) Cd * cos(phi)) * r(i); T T dT; Q Q dQ; end P Q * Omega; end function F bemResidual(x, ri, R, c, theta, Vinf, Omega, rho, Nb) a x(1); ap x(2); va Vinf * (1 a); vt Omega * ri * (1 - ap); Vrel sqrt(va^2 vt^2); phi atan2(va, vt); alpha theta - phi; [Cl, Cd] lookupAirfoil(alpha); dT_el Nb * 0.5 * rho * Vrel^2 * c * ... (Cl * cos(phi) - Cd * sin(phi)); dQ_el Nb * 0.5 * rho * Vrel^2 * c * ... (Cl * sin(phi) Cd * cos(phi)) * ri; dT_mom 4 * pi * rho * ri * Vinf^2 * (1 a) * a; dQ_mom 4 * pi * rho * ri^3 * Vinf * Omega * (1 a) * ap; F [dT_el - dT_mom; dQ_el - dQ_mom]; endfsolve需要优化工具箱如果没有工具箱可以自己写一个二维牛顿迭代核心就是数值差分算雅可比矩阵逻辑并不复杂。但对于大多数场景直接使用fsolve更省事把容差收紧到1e-12能明显减少积分累积误差。lookupAirfoil函数没有写在骨架里因为不同项目的翼型表格式不一样。我建议在函数里用interp1对alpha做线性插值同时加一句边界判断把超出范围的迎角做受限外插具体做法后面会说。3.3 前进比扫描主程序给定恒定转速前进比JVinf/(nD)所以扫描前进比就是在同一个转速下改变来流速度。主程序循环调用bemSolve即可Jlist 0.05:0.05:1.0; for k 1:length(Jlist) Vinf Jlist(k) * n_rps * D; [T(k), Q(k), P(k), a_last, ap_last] ... bemSolve(r, c, thetaDeg, Vinf, Omega, rho, Nb); if k 1 % 保存初值供下一工况 end eta(k) T(k) * Vinf / P(k); end这里有几个工程细节。第一J0不能直接算因为动量方程里Vinf0会让dT_mom变成0迭代失去依据。我的做法是从J0.05起步然后利用相邻工况的解作为初值一路向前推进。第二每个工况之间传递初值非常关键如果每次都从固定初值重新开始中高前进比区域很可能收敛到错误分支或直接不收敛。第三效率η的定义是T×Vinf/P在J很低时分子很小效率接近0此时算出来的η噪声会很大画图前可以先做一次平滑。4. 不同前进比下性能到底怎么变化4.1 悬停点J0J0对应来流速度为零的悬停状态这类工况在固定翼巡航评估里不常用但在多旋翼和垂直起降方案里是必须校核的设计点。此时桨盘完全靠诱导速度把静止空气加速往下推当地迎角最大拉力也最大功率达到峰值而推进效率按定义严格等于零因为有用功T×Vinf为零。工程上更关注的是另一个指标单位功率能产生多大拉力即悬停效率或FoM。J0对BEMT是个数值陷阱。动量方程在Vinf0时右端是0叶素方程仍然有载荷两边无法匹配解不再唯一。我的处理是给一个很小的来流速度比如0.2米/秒把它当作“准悬停”计算结果与真实悬停的偏差通常很小因为诱导速度远大于这一微小的来流扰动。如果你有实验悬停数据用这个点做校验更能暴露翼型数据和动量假设的系统误差。4.2 巡航点效率峰值在中间随着J增大来流对桨叶的冲刷作用增强当地迎角逐渐减小升力系数下降所以推力和功率都下降。但这两者下降的速度不一样推力大致线性下降功率在高J区间下降得更快结果是效率曲线先上升后下降在某个中间J处出现峰值。峰值出现在哪里直接取决于弦长分布、扭角分布和翼型特性。常规设计里如果桨叶扭转合理峰值效率落在J0.5到J0.8之间。做动力匹配时我会把峰值效率对应的J乘以转速和直径换算成巡航速度再和飞行器的设计巡航速度对比。如果峰值点偏快了就增加桨距角或降低转速如果偏慢了就减小桨距角。这个匹配过程用BEMT扫一遍J就完成比靠经验调桨高效得多。需要提醒的是效率曲线峰值往往比较平缓峰值附近±0.1的J区间性能差别并不大选择工作点还应兼顾噪声、电机转速限制和低推力的稳定性。4.3 风车点拉力归零以后继续增大J拉力会持续下降直到某一前进比下推力变为零这个点通常叫风车点。此时螺旋桨不消耗功率来流反而带动桨叶旋转产生负扭矩相当于一个阻力面或回收装置。风车点两侧的物理特性差异很大推力为正时桨是推进器功率从轴端输入推力为负时桨变成涡轮可能从气流中提取能量。BEMT在这个区间仍然能算出结果但翼型表往往需要包含负迎角数据否则Cl/Cd插值会失真。如果在设计任务里只需要评估正常飞行包线扫到推力和效率变为负值之前就可以停。但如果项目涉及自转下滑、能量回收、空中风车启动风车点附近的特性就很重要。我实际测过几个小直径螺旋桨风车点对应的C_T通常随转速变化不大用BEMT预测的零推力J和实验误差大约在百分之几足够支撑概念评估。下面对典型前进比区间的特征做个汇总方便直观对照前进比范围推力趋势功率趋势效率表现主要关注点J接近0最大最大接近0悬停拉力、电机功率上限中J区间下降下降较快有峰值巡航匹配、最高效率点接近风车点接近0接近0趋近0零推力点、阻力特性超过风车点变负变负定义失效风车回收、自转下滑5. 调试过程与常见问题速查5.1 迭代不收敛的根因与初值策略我在这套代码上遇到最多的报错就是fsolve提前退出提示找不到解。绝大多数原因不在方程本身而在初值给得不好。BEMT的每个叶素方程可能存在多解或非物理解初值落在错误的吸引域里就会收敛到一个看起来合理但完全错误的诱导因子甚至干脆发散。我的固定做法分两步。第一步始终用相邻叶素的解作为当前叶素的初值因为弦长和扭转沿展向变化是连续的相邻剖面的诱导速度也应当连续。第二步在扫描前进比时把上一工况的所有解作为下一工况的初值每步J增量控制在0.05以内。如果遇到跳跃性变化比如从悬停直接切到大前进比初值离解太远就会失败。遇到失败时我的策略是先把J步长减小回退几步重新算而不是硬调fsolve容差。5.2 翼型数据插值和外插的坑翼型极曲线是查表数据默认只覆盖某个迎角区间。小型螺旋桨常用翼型表可能只覆盖-5度到15度但在低前进比高载荷下根部叶素的迎角轻松超过20度在大前进比风车工况下迎角又可能跑到-10度以下。查询返回NaN后残差函数直接失效。我在lookupAirfoil里做了三层保护。第一层采用interp1的线性插值不放高阶样条因为样条容易在表两端震荡。第二层对超范围迎角做线性外插外插斜率取表尾两点的斜率而不是任意延长。第三层如果外插后的Cl或Cd出现在物理上不可能的范围比如Cl远大于翼型最大升力系数就返回一个带警告的截断值作为提示。调试模式下我会把每个叶素的alpha和Cl、Cd打印出来看到底是哪一段在超界再决定是补充翼型数据还是限制前进比范围。5.3 动量理论断裂与涡环状态BEMT在中等以上载荷或者大诱导速度时会遇到一个理论缺陷动量方程的解可能不唯一甚至不存在这叫动量理论断裂。最典型的情况是轴向导诱导因子a大于0.3左右对应桨盘处于“涡轮状态”或“涡环状态”。在这种状态下流管假设不再成立直接用原始动量方程迭代要么算不出解要么算出非物理的突变。工程上最常见的做法是引入Glauert修正在a超过某个阈值后不再使用原始动量方程而是用一个经验公式把推力系数与a平滑连接避免迭代跨越不连续区。很多开源BEM代码里都有这段“超阈值修正”我第一次自己写时漏掉了结果在低前进比大拉力点上方出现了一段完全错误的负推力。如果你在J较小时发现推力曲线出现尖刺首先就该怀疑这里。5.4 结果合理性校验BEMT代码本身不复杂复杂的是判断结果是否可信。我每次跑完都会做几件检查。首先功率必须满足PQ×Ω因为我是用扭矩乘角速度算功率这个等式是恒真的但可以检查积分时有没有漏乘叶片数或半径。其次效率ηT×Vinf/P不能超过1更不能出现远大于1的情况一旦出现必然是高J段推力或功率的计算出了问题。更系统的做法是算无量纲系数并与经验区间对比。典型的小型螺旋桨在峰值效率附近C_T大约在0.05到0.15之间C_P在0.03到0.1之间具体值随几何差异很大但数量级如果偏离太远就要回头查几何单位。我遇到过用厘米和米混用的版本所有结果差了100倍好在无量纲系数一下就暴露了。最后检查各叶素的诱导因子和当地迎角沿展向是否光滑。如果某个叶素和邻居差了很大一截几乎可以断定那个点收敛到了错误分支。6. 这个工具后续还能怎么扩展6.1 变螺距与动力匹配优化现在这个程序还只能评估“给定几何、恒定转速、扫前进比”的性能稍微改一下就可以做变螺距优化。把安装角thetaDeg设计成变量让它随转速或来流速度调整BEMT可以快速生成不同桨距角下的效率云图。结合电机的效率曲线和电池/发动机的功率限制我就能搜索出一条最优控制策略在低速段用大桨距、高速段用小桨距或者反过来。这一步做得好整机续航提升百分之几并不夸张。6.2 滑流干扰和旋翼场景迁移BEMT假设来流均匀但实际螺旋桨经常装在机翼前方。把均匀来流Vinf替换成某个展向位置的速度剖面比如机翼下洗产生的下洗速度、机身附面层的减速区、甚至另一组旋翼的尾流速度BEMT就能做一阶滑流干扰评估。这个扩展对做分布式电推进的布局论证特别有用虽然精度不如CFD但能在几秒钟内给出“哪个展向位置放桨更有利”的方向性结论。把固定翼巡航的“前进比”概念换一换BEMT同样可以用于直升机旋翼的悬停和前飞性能只是需要引入挥舞角和来流倾斜角。多旋翼场景里J0就是悬停大家最关心的拉力功率比可以直接输出。如果还想加地面效应可以在动量方程中引入镜像速度项这属于另一个独立的扩展模块。6.3 把代码改造成工程组件目前这个Matlab函数虽然能用但还称不上一个工具。我建议把它封装成独立的类或函数库输入输出统一用结构体让几何、翼型表、工况参数分开传入。这样后面可以批量扫描几十个候选桨叶自动记录每个桨的峰值效率、零推力前进比、悬停拉力输出一张选型对比表。如果团队里其他人用Python或者C同样的核心逻辑可以逐行迁移因为BEMT本身并不依赖Matlab的特殊语法本质就是一组循环和查表。我自己跑完这套Matlab代码最大的体会是BEMT的价值不在精确而在把设计变量和性能之间的因果链压缩到可以一眼看穿。每当拿到一套新螺旋桨几何我第一件事就是先扫一遍J看看效率和功率峰值落在哪里再决定要不要上CFD或实验。如果你照上面的思路写出来遇到不收敛别急着调求解器先检查翼型数据表覆盖范围和初值传递多数情况下问题都出在这两处。 叶素动量理论BEMT：Matlab实现螺旋桨性能快速评估
做小型螺旋桨性能评估时我最常用的第一轮工具并不是CFD而是叶素动量理论BEMT。标题里那个“动量理论剖面刀片方法”其实就是这个理论的另一种直译国内工程圈更习惯叫叶素动量理论或者片条理论。它的任务很明确给定一套螺旋桨几何在转速恒定的前提下把来流速度从零逐步往上推算出不同前进比下的拉力、扭矩、功率和效率。这个工作用Matlab实现并不难难在诱导速度的迭代收敛处理以及把翼型气动数据、几何参数和动量约束接在一起不出错。我把实际跑通这套流程的经验整理在下面希望能帮你少踩几个坑。1. 动手之前先搞懂BEMT到底在算什么1.1 为什么第一轮评估用BEMT而不是CFD螺旋桨气动分析的方法其实不少从经验公式、动量理论、叶素动量理论到全三维RANS、Lattice Boltzmann甚至直接上风洞台架测力。放在同一个天平上比BEMT的定位是“方案阶段的快速筛选器”。我做预研时选它有三个理由。第一计算量小到可以忽略把一片桨叶沿展向切几十个叶素每个叶素解一个二维非线性方程Matlab跑一个工况就是毫秒级扫一遍前进比也就几秒钟。第二物理因果链非常清楚弦长变化、扭转变化、转速变化、来流变化都会直接映射到某个剖面的迎角、升阻力系数和整桨性能曲线上调试时一眼就能看出问题在哪。第三输入数据全是设计阶段已经确定的参数不需要网格、不需要湍流模型、不需要壁面处理。CFD的优势在细节比如分离流、桨尖涡、三维失速、机翼与滑流的强干扰这些BEMT天然给不了。但在“给定几何、恒定转速、扫前进比”这类任务里CFD的单工况耗时从小时到天不等几何稍有改动又要重新来一轮网格根本不划算。我的经验是先用BEMT把设计空间快速缩小选定少数几个候选几何再让CFD去做精确校核最后实验验证。这比一上来就全流程CFD高效得多。方法单工况耗时输入要求适用范围典型用途BEMT毫秒级弦长、扭角、翼型极曲线附着流为主、载荷中等选型、趋势分析、动力匹配CFD RANS小时到数天三维几何、网格、湍流模型分离流、复杂干扰、流场细节最终校核、诊断特殊现象台架试验天到周样机、台架、传感器真实完备定型验证、标定1.2 动量理论和叶素理论是怎么互补的“BEMT”这个名字里的两个部分不是在讲同一件事。动量理论把螺旋桨简化成一个无限薄的致动盘假设来流均匀、不可压、无粘并忽略尾流旋转的细节只研究桨盘前后的速度变化和压力突跳。它给出的是宏观结论桨盘从流场中吸收了动量气流通过桨盘后轴向速度增加同时获得与旋转方向相反的切向速度于是形成推力和反扭矩。动量理论的优点是不依赖桨叶长什么样缺点是它根本不认识桨叶几何。叶素理论则反过来。它把桨叶沿展向切成一个个二维剖面每个剖面看成一段独立翼型用当地相对速度、迎角和二维翼型升阻力数据算出剖面气动力再沿展向积分。叶素的优点是能够考虑真实的弦长、扭转和翼型缺点是它没有回答“来流进入桨盘后到底被加速了多少”。这个加速量恰恰决定了当地迎角。BEMT把两边接起来动量理论给出轴向和周向诱导速度与载荷的关系叶素理论给出桨叶几何产生的载荷关系两者必须相等于是诱导速度被唯一确定。这就是它能把几何参数变成性能数据的核心机制。标题里的“剖面刀片方法”强调的正是叶素那一侧的“切剖面”思想翻译上可以接受但习惯上还是叫叶素动量理论。1.3 输入、输出和适用边界输入部分非常明确。几何上需要螺旋桨半径、沿展向的弦长分布和安装角分布安装角实际上就是每个剖面的几何扭转角。气动上需要翼型极曲线也就是升力系数Cl和阻力系数Cd随迎角α的查表数据。工况上需要转速rpm、来流速度Vinf、来流密度ρ和叶片数N。输出部分首先是每个叶素的诱导因子a和a、当地迎角、受力然后沿展向积分得到整桨推力T、扭矩Q、功率P和推进效率η。为了方便跨工况对比我还会额外算无量纲系数推进系数C_T T/(ρ n²D⁴)、功率系数C_P P/(ρ n³D⁵)、前进比J Vinf/(nD)其中n是每秒转速D是直径。适用边界也必须要讲清楚。BEMT建立在均匀定常来流、附着流、低速亚声速这些假设上所以桨尖马赫数最好不要超过0.7翼型数据要覆盖实际遇到的雷诺数范围且不需要处理失速后的复杂分离。如果来流很不均匀比如螺旋桨装在机翼后方或者在大攻角状态下工作BEMT只能给出趋势级的结果。这个约束不是缺点而是工具本身的定位快速、可解释、覆盖大多数正常设计点。2. 把BEMT公式拆到能直接写代码的程度2.1 截面速度三角与微元气动力取半径r处的一个小段厚度为dr。当地轴向速度是来流速度加上诱导速度写成Vinf×(1a)这里的a叫轴向诱导因子。切向速度稍微绕一些桨叶以角速度Ω旋转当地旋转速度为Ωr但由于桨盘把动量给了气流气流获得与旋转方向相反的切向速度所以在叶素坐标系中切向来流速度表示为Ωr×(1-a)a是切向诱导因子。不同教材对a的符号定义会有差异但只要自己在代码里保持一致结果不会出问题。相对速度Vrel就是由轴向和切向两个分量合成的Vaxial Vinf * (1 a); Vtang Omega * r * (1 - ap); Vrel sqrt(Vaxial^2 Vtang^2); phi atan2(Vaxial, Vtang); alpha theta - phi;phi是流入角theta是几何安装角alpha就是当地迎角。有了alpha就能从翼型极曲线插值得到Cl和Cd。单个叶素的升力是0.5ρVrel²c dr×Cl阻力是0.5ρVrel²c dr×Cd。把升力阻力分别沿轴向和切向投影轴向分量形成拉力切向分量形成对桨毂的旋转阻力矩。乘上叶片数N就得到整个桨盘在这一小段上的气动贡献。这里有个地方特别容易写错升力方向定义为垂直于相对速度阻力方向定义为平行于相对速度。投影到轴向时需要用流入角phi的正弦余弦公式是dT_el N×0.5ρVrel²c dr(Cl cosφ − Cd sinφ)dQ_el N×0.5ρVrel²c dr(Cl sinφ Cd cosφ)×r。如果角度定义反了推力变成阻力整个迭代直接发散。2.2 动量方程那一侧的约束同一小段还可以从动量角度建立方程。单位时间流过半径r处环形微元的质量流量是ρ×2πr dr×Vinf(1a)。气流通过桨盘后轴向速度增量在下游逐渐达到2aVinf由动量定理轴向力的增量就是质量流量乘以速度增量整理后得到dT_mom 4πρr Vinf²(1a)a dr。这是BEMT中轴向诱导因子a的动量约束。切向方向同理。气流通过桨盘后获得旋转速度角动量变化产生扭矩。最常见的表达形式是dQ_mom 4πρr³ Ω Vinf(1a)a dr。这个式子里的倍率在不同教材中可能差一个2的因子因为对“盘面处的切向速度”和“尾流远场切向速度”的定义不同。我在代码里采用上面这组形式推导过程中只要保证质量流量和速度增量的定义一致最终结果就是自洽的。从物理上看动量方程说的是“流场给定诱导速度因此能容纳多大载荷”叶素方程说的是“桨叶几何在某个诱导速度下能产生多大载荷”。两者必须相等这正是迭代求解的出发点。2.3 为什么必须迭代而不是直接解公式如果把叶素方程和动量方程合起来看a和a出现在速度合成、流入角、迎角、Cl/Cd插值里而Cl/Cd通常是一张离散的查表和α的非线性函数。整个方程组是强非线性隐式方程不可能像线性方程组那样一步解出闭式解。实际操作中几乎都用数值迭代。最稳妥的做法是对每个叶素列两个残差方程F1(a,a) dT_el − dT_mom F2(a,a) dQ_el − dQ_mom然后用二维牛顿法或者Matlab的fsolve把a和a解出来。注意每个叶素的a和a只对这个半径位置有效所以整个桨叶在计算时是逐段独立求解的唯一的联系只是把力积分起来。这个“逐段独立”的好处是代码结构清晰坏处是一旦某一段解飞了整条曲线就会在对应位置出现一个尖刺。3. Matlab实现从几何数据到性能曲线3.1 输入数据怎么准备我习惯把几何参数放在一个表里每一行对应一个径向截面归一化半径rNorm、弦长c、安装角thetaDeg。翼型极曲线单独放两个文件或者一个矩阵第一列是alpha第二列Cl第三列Cd。Matlab里用readmatrix读入方便替换翼型。下面的示例描述一个半径0.254米、两片桨叶的几何采样30个截面R 0.254; % 桨半径单位米 r linspace(0.08*R, R, 30); c 0.05 * ones(size(r)) .* (1 - 0.5 * r / R); % 弦长从根到尖递减 thetaDeg 30 - 20 * r / R; % 安装角从根到尖递减 Nb 2; % 叶片数 rho 1.225; % 来流密度 rpm 6000; Omega rpm * 2 * pi / 60; D 2 * R; n_rps rpm / 60;实际的几何数据通常来自CAD或三维扫描弦长和扭角分布不一定是光滑函数但只要每个径向位置都有对应值BEMT完全可以直接处理。翼型极曲线需要注意迎角范围。低速小型螺旋桨常用的翼型在-5度到15度范围内能有完整的Cl/Cd表但如果扫前进比时当地迎角超出了这个范围插值函数会返回NaN迭代直接失败。这个问题我在第5节专门讲。3.2 核心求解器代码骨架核心求解函数接收半径数组、弦长数组、安装角数组和工况参数返回整桨推力、扭矩、诱导因子分布。残差函数就是上一节的F1和F2。这里给出一个能直接修改使用的骨架function [T, Q, P, a, ap] bemSolve(r, c, thetaDeg, Vinf, Omega, rho, Nb) n length(r); theta deg2rad(thetaDeg); a zeros(n,1); ap zeros(n,1); T 0; Q 0; x0 [0.05; 0.01]; % 第一个叶素用保守初值 for i 1:n if i 1 x0 [a(i-1); ap(i-1)]; % 用上一叶素的解作为初值 end opt optimoptions(fsolve, Display, off, ... MaxIterations, 600, ... FunctionTolerance, 1e-12, ... StepTolerance, 1e-12); x fsolve((x) bemResidual(x, r(i), R, c(i), theta(i), ... Vinf, Omega, rho, Nb), x0, opt); a(i) x(1); ap(i) x(2); va Vinf * (1 a(i)); vt Omega * r(i) * (1 - ap(i)); Vrel sqrt(va^2 vt^2); phi atan2(va, vt); alpha theta(i) - phi; [Cl, Cd] lookupAirfoil(alpha); dT Nb * 0.5 * rho * Vrel^2 * c(i) * ... (Cl * cos(phi) - Cd * sin(phi)); dQ Nb * 0.5 * rho * Vrel^2 * c(i) * ... (Cl * sin(phi) Cd * cos(phi)) * r(i); T T dT; Q Q dQ; end P Q * Omega; end function F bemResidual(x, ri, R, c, theta, Vinf, Omega, rho, Nb) a x(1); ap x(2); va Vinf * (1 a); vt Omega * ri * (1 - ap); Vrel sqrt(va^2 vt^2); phi atan2(va, vt); alpha theta - phi; [Cl, Cd] lookupAirfoil(alpha); dT_el Nb * 0.5 * rho * Vrel^2 * c * ... (Cl * cos(phi) - Cd * sin(phi)); dQ_el Nb * 0.5 * rho * Vrel^2 * c * ... (Cl * sin(phi) Cd * cos(phi)) * ri; dT_mom 4 * pi * rho * ri * Vinf^2 * (1 a) * a; dQ_mom 4 * pi * rho * ri^3 * Vinf * Omega * (1 a) * ap; F [dT_el - dT_mom; dQ_el - dQ_mom]; endfsolve需要优化工具箱如果没有工具箱可以自己写一个二维牛顿迭代核心就是数值差分算雅可比矩阵逻辑并不复杂。但对于大多数场景直接使用fsolve更省事把容差收紧到1e-12能明显减少积分累积误差。lookupAirfoil函数没有写在骨架里因为不同项目的翼型表格式不一样。我建议在函数里用interp1对alpha做线性插值同时加一句边界判断把超出范围的迎角做受限外插具体做法后面会说。3.3 前进比扫描主程序给定恒定转速前进比JVinf/(nD)所以扫描前进比就是在同一个转速下改变来流速度。主程序循环调用bemSolve即可Jlist 0.05:0.05:1.0; for k 1:length(Jlist) Vinf Jlist(k) * n_rps * D; [T(k), Q(k), P(k), a_last, ap_last] ... bemSolve(r, c, thetaDeg, Vinf, Omega, rho, Nb); if k 1 % 保存初值供下一工况 end eta(k) T(k) * Vinf / P(k); end这里有几个工程细节。第一J0不能直接算因为动量方程里Vinf0会让dT_mom变成0迭代失去依据。我的做法是从J0.05起步然后利用相邻工况的解作为初值一路向前推进。第二每个工况之间传递初值非常关键如果每次都从固定初值重新开始中高前进比区域很可能收敛到错误分支或直接不收敛。第三效率η的定义是T×Vinf/P在J很低时分子很小效率接近0此时算出来的η噪声会很大画图前可以先做一次平滑。4. 不同前进比下性能到底怎么变化4.1 悬停点J0J0对应来流速度为零的悬停状态这类工况在固定翼巡航评估里不常用但在多旋翼和垂直起降方案里是必须校核的设计点。此时桨盘完全靠诱导速度把静止空气加速往下推当地迎角最大拉力也最大功率达到峰值而推进效率按定义严格等于零因为有用功T×Vinf为零。工程上更关注的是另一个指标单位功率能产生多大拉力即悬停效率或FoM。J0对BEMT是个数值陷阱。动量方程在Vinf0时右端是0叶素方程仍然有载荷两边无法匹配解不再唯一。我的处理是给一个很小的来流速度比如0.2米/秒把它当作“准悬停”计算结果与真实悬停的偏差通常很小因为诱导速度远大于这一微小的来流扰动。如果你有实验悬停数据用这个点做校验更能暴露翼型数据和动量假设的系统误差。4.2 巡航点效率峰值在中间随着J增大来流对桨叶的冲刷作用增强当地迎角逐渐减小升力系数下降所以推力和功率都下降。但这两者下降的速度不一样推力大致线性下降功率在高J区间下降得更快结果是效率曲线先上升后下降在某个中间J处出现峰值。峰值出现在哪里直接取决于弦长分布、扭角分布和翼型特性。常规设计里如果桨叶扭转合理峰值效率落在J0.5到J0.8之间。做动力匹配时我会把峰值效率对应的J乘以转速和直径换算成巡航速度再和飞行器的设计巡航速度对比。如果峰值点偏快了就增加桨距角或降低转速如果偏慢了就减小桨距角。这个匹配过程用BEMT扫一遍J就完成比靠经验调桨高效得多。需要提醒的是效率曲线峰值往往比较平缓峰值附近±0.1的J区间性能差别并不大选择工作点还应兼顾噪声、电机转速限制和低推力的稳定性。4.3 风车点拉力归零以后继续增大J拉力会持续下降直到某一前进比下推力变为零这个点通常叫风车点。此时螺旋桨不消耗功率来流反而带动桨叶旋转产生负扭矩相当于一个阻力面或回收装置。风车点两侧的物理特性差异很大推力为正时桨是推进器功率从轴端输入推力为负时桨变成涡轮可能从气流中提取能量。BEMT在这个区间仍然能算出结果但翼型表往往需要包含负迎角数据否则Cl/Cd插值会失真。如果在设计任务里只需要评估正常飞行包线扫到推力和效率变为负值之前就可以停。但如果项目涉及自转下滑、能量回收、空中风车启动风车点附近的特性就很重要。我实际测过几个小直径螺旋桨风车点对应的C_T通常随转速变化不大用BEMT预测的零推力J和实验误差大约在百分之几足够支撑概念评估。下面对典型前进比区间的特征做个汇总方便直观对照前进比范围推力趋势功率趋势效率表现主要关注点J接近0最大最大接近0悬停拉力、电机功率上限中J区间下降下降较快有峰值巡航匹配、最高效率点接近风车点接近0接近0趋近0零推力点、阻力特性超过风车点变负变负定义失效风车回收、自转下滑5. 调试过程与常见问题速查5.1 迭代不收敛的根因与初值策略我在这套代码上遇到最多的报错就是fsolve提前退出提示找不到解。绝大多数原因不在方程本身而在初值给得不好。BEMT的每个叶素方程可能存在多解或非物理解初值落在错误的吸引域里就会收敛到一个看起来合理但完全错误的诱导因子甚至干脆发散。我的固定做法分两步。第一步始终用相邻叶素的解作为当前叶素的初值因为弦长和扭转沿展向变化是连续的相邻剖面的诱导速度也应当连续。第二步在扫描前进比时把上一工况的所有解作为下一工况的初值每步J增量控制在0.05以内。如果遇到跳跃性变化比如从悬停直接切到大前进比初值离解太远就会失败。遇到失败时我的策略是先把J步长减小回退几步重新算而不是硬调fsolve容差。5.2 翼型数据插值和外插的坑翼型极曲线是查表数据默认只覆盖某个迎角区间。小型螺旋桨常用翼型表可能只覆盖-5度到15度但在低前进比高载荷下根部叶素的迎角轻松超过20度在大前进比风车工况下迎角又可能跑到-10度以下。查询返回NaN后残差函数直接失效。我在lookupAirfoil里做了三层保护。第一层采用interp1的线性插值不放高阶样条因为样条容易在表两端震荡。第二层对超范围迎角做线性外插外插斜率取表尾两点的斜率而不是任意延长。第三层如果外插后的Cl或Cd出现在物理上不可能的范围比如Cl远大于翼型最大升力系数就返回一个带警告的截断值作为提示。调试模式下我会把每个叶素的alpha和Cl、Cd打印出来看到底是哪一段在超界再决定是补充翼型数据还是限制前进比范围。5.3 动量理论断裂与涡环状态BEMT在中等以上载荷或者大诱导速度时会遇到一个理论缺陷动量方程的解可能不唯一甚至不存在这叫动量理论断裂。最典型的情况是轴向导诱导因子a大于0.3左右对应桨盘处于“涡轮状态”或“涡环状态”。在这种状态下流管假设不再成立直接用原始动量方程迭代要么算不出解要么算出非物理的突变。工程上最常见的做法是引入Glauert修正在a超过某个阈值后不再使用原始动量方程而是用一个经验公式把推力系数与a平滑连接避免迭代跨越不连续区。很多开源BEM代码里都有这段“超阈值修正”我第一次自己写时漏掉了结果在低前进比大拉力点上方出现了一段完全错误的负推力。如果你在J较小时发现推力曲线出现尖刺首先就该怀疑这里。5.4 结果合理性校验BEMT代码本身不复杂复杂的是判断结果是否可信。我每次跑完都会做几件检查。首先功率必须满足PQ×Ω因为我是用扭矩乘角速度算功率这个等式是恒真的但可以检查积分时有没有漏乘叶片数或半径。其次效率ηT×Vinf/P不能超过1更不能出现远大于1的情况一旦出现必然是高J段推力或功率的计算出了问题。更系统的做法是算无量纲系数并与经验区间对比。典型的小型螺旋桨在峰值效率附近C_T大约在0.05到0.15之间C_P在0.03到0.1之间具体值随几何差异很大但数量级如果偏离太远就要回头查几何单位。我遇到过用厘米和米混用的版本所有结果差了100倍好在无量纲系数一下就暴露了。最后检查各叶素的诱导因子和当地迎角沿展向是否光滑。如果某个叶素和邻居差了很大一截几乎可以断定那个点收敛到了错误分支。6. 这个工具后续还能怎么扩展6.1 变螺距与动力匹配优化现在这个程序还只能评估“给定几何、恒定转速、扫前进比”的性能稍微改一下就可以做变螺距优化。把安装角thetaDeg设计成变量让它随转速或来流速度调整BEMT可以快速生成不同桨距角下的效率云图。结合电机的效率曲线和电池/发动机的功率限制我就能搜索出一条最优控制策略在低速段用大桨距、高速段用小桨距或者反过来。这一步做得好整机续航提升百分之几并不夸张。6.2 滑流干扰和旋翼场景迁移BEMT假设来流均匀但实际螺旋桨经常装在机翼前方。把均匀来流Vinf替换成某个展向位置的速度剖面比如机翼下洗产生的下洗速度、机身附面层的减速区、甚至另一组旋翼的尾流速度BEMT就能做一阶滑流干扰评估。这个扩展对做分布式电推进的布局论证特别有用虽然精度不如CFD但能在几秒钟内给出“哪个展向位置放桨更有利”的方向性结论。把固定翼巡航的“前进比”概念换一换BEMT同样可以用于直升机旋翼的悬停和前飞性能只是需要引入挥舞角和来流倾斜角。多旋翼场景里J0就是悬停大家最关心的拉力功率比可以直接输出。如果还想加地面效应可以在动量方程中引入镜像速度项这属于另一个独立的扩展模块。6.3 把代码改造成工程组件目前这个Matlab函数虽然能用但还称不上一个工具。我建议把它封装成独立的类或函数库输入输出统一用结构体让几何、翼型表、工况参数分开传入。这样后面可以批量扫描几十个候选桨叶自动记录每个桨的峰值效率、零推力前进比、悬停拉力输出一张选型对比表。如果团队里其他人用Python或者C同样的核心逻辑可以逐行迁移因为BEMT本身并不依赖Matlab的特殊语法本质就是一组循环和查表。我自己跑完这套Matlab代码最大的体会是BEMT的价值不在精确而在把设计变量和性能之间的因果链压缩到可以一眼看穿。每当拿到一套新螺旋桨几何我第一件事就是先扫一遍J看看效率和功率峰值落在哪里再决定要不要上CFD或实验。如果你照上面的思路写出来遇到不收敛别急着调求解器先检查翼型数据表覆盖范围和初值传递多数情况下问题都出在这两处。