多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

Matlab实现螺旋桨BEMT性能分析:快速计算推力功率效率曲线

Matlab实现螺旋桨BEMT性能分析:快速计算推力功率效率曲线 做螺旋桨性能分析的人大多会碰到同一个问题手里有一颗具体的桨想知道它在不同来流速度下的推力、功率和效率表现但又不想因为“先估一版看看趋势”就开一堆网格、跑一周CFD。叶片单元动量理论Blade Element Momentum TheoryBEMT就是干这个用的——用几百行代码几分钟内把一颗给定几何的螺旋桨在任意前进比下的性能曲线算出来。这篇内容我基于Matlab实现了一套完整流程从桨叶几何离散、翼型数据插值到诱导速度迭代求解、转速恒定下前进比扫描最终输出推力系数、功率系数和效率曲线整个思路可以直接迁移到船桨、风机叶片甚至涵道风扇的初步设计里。我会把整套求解逻辑拆开讲清楚不仅仅是贴一段能跑的代码而是把每个公式为什么长这样、迭代时为什么容易发散、什么时候该信BEMT的结果、什么时候该警惕——这些实际操作中才能积累的判断一并写出来。如果你正在用Matlab做螺旋桨性能预估或者刚开始学BEMT、想知道怎么把理论公式变成可用的计算程序这篇应该能帮你省不少弯路。1. 为什么选叶片单元动量理论做性能快速分析1.1 动量理论和叶素理论各自的局限螺旋桨性能分析的经典路线有两条。一条是动量理论把螺旋桨当成一个均匀的“致动盘”气流穿过桨盘时速度突变、压力突变用一维流动的动量守恒和能量守恒推算推力、功率。这个理论最大的优点是算得快、物理图像干净但结果只有全局性能完全看不出桨叶哪一段出力大、哪一段在拖后腿而且它假设桨盘载荷均匀实际桨叶根部、尖部的载荷差异很大直接套用误差相当可观。另一条是叶素理论把桨叶沿径向切成几十个小段每一段近似看成无限长二维翼型用当地的攻角和来流速度查翼型升力、阻力系数然后积出总推力和扭矩。这条路线能反映几何细节弦长分布、扭转分布都能放进去但它有一个致命缺陷——没有把螺旋桨对气流的“诱导作用”算进来。桨叶把空气向后推空气反过来也会改变桨叶处的实际流速方向和大小忽略第二遍影响攻角就估不准推力自然也不准。BEMT就是把这两条缝起来用动量理论求诱导速度用叶素理论算各叶素的载荷两者通过诱导因子互相迭代直到收敛。算出来的结果既有几何分辨率又有物理一致性这正是它成为螺旋桨和风力机工程设计标准方法的原因。1.2 前进比J的定义与恒定转速工况标题里强调“恒定转速”是因为我们要扫的变量是前进比J。无量纲前进比的定义式是J V / (n × D)其中V是来流速度m/sn是螺旋桨转速转/秒D是桨盘直径m。这个量可以理解为“螺旋桨每转一圈前进了多少个直径”。J越小意味着桨叶转动越快而前进越慢相当于高负荷工况J越大相当于接近顺桨、载荷变轻。固定转速、改变来流速度来扫描前进比对应的是实际的“保持发动机转速、飞机或船速变化”场景。比如无人机螺旋桨在试车台上锁定了电机转速来流风洞速度从0开始逐步加大每一档就是一个J值。另一种做法是固定来流、改变转速物理过程相似但画出来的曲线横坐标含义会稍有差别——固定转速扫描的工况雷诺数变化不大比较干净这也是我为什么推荐先从恒定转速做起。1.3 BEMT能输出什么性能量BEMT的输出不止推力一个数。基于诱导速度收敛结果我们可以汇总出推力T由各叶素升阻力的轴向分量积分而来扭矩Q由各叶素升阻力的切向分量乘力臂积分而来轴功率P 2π × n × Q Ω × Q推力系数CT T / (ρ × n² × D⁴)功率系数CP P / (ρ × n³ × D⁵)效率η J × CT / CP。这里我习惯把效率写成η J × CT / CP因为这样直接对应螺旋桨的推进效率定义——有用功是T×V输入功是P约掉无量纲系数后正好等于J×CT/CP。最终画出来的曲线就是CT、CP、η随J的变化这是桨性能的“身份证”后续做尺度换算、采购对比、飞行性能预估都靠它。1.4 为什么不直接上CFDCFD当然更“真实”但成本摆在那里三维桨叶网格从几十万到几百万单元边界层要加密转捩模型参数要标定一个工况算几小时到几天还不一定收敛。更关键的是CFD只能回答“给定的这个工况长什么样”如果要探索几何参数变化对性能的影响比如扭角改3度、弦长前移5%每个方案都要重新划分网格重新算设计迭代根本跑不动。BEMT的定位不是替代CFD而是在概念设计和参数扫描阶段充当一个“足够准的快筛工具”。我个人的习惯流程是先用BEMT把几十个前进比算完找出效率最高的工况区间和需要关注的攻角分布再用CFD针对最优工况做精细校核。这样CFD算的每一炮都打在关键点上而不是盲人摸象。对科研写论文、做课程设计或者工程前期选型来说BEMT的精度完全够用。2. 几何模型与翼型数据的前期准备2.1 桨叶径向分段与坐标系约定BEMT的第一步是把桨叶从桨毂到桨尖沿径向离散成N个叶素。坐标系我统一用柱坐标x轴沿桨轴指向下游r轴沿径向向外θ轴为旋转方向。对第i个叶素取其中点半径r_i叶素宽度dr_i R / N。代码里我通常把配平间距匀分在半径上其实按对数分布或者靠桨尖加密更好但作为初步分析均匀分段已经够用。叶素半径从多少开始很重要。桨根附近有桨毂、法兰、根部的整流套气动面不完整那里的叶素要么升力系数不可信要么来流速度低、贡献很小。我习惯把积分起点设在0.15R左右也就是最内测的叶素在半径的15%处开始。不同桨叶几何稍有差别你可以看一眼实际桨叶如果发现根部气动面从20%才开始就把起点相应调高。2.2 弦长和扭角分布的输入方式给定螺旋桨几何形状落到程序里其实就是两条沿半径分布的曲线弦长c(r)和几何扭转角β(r)。这里的β是叶素弦线与旋转平面之间的夹角单位用弧度因为后面公式里要跟φ直接相加减。实际测量或者拿图纸时常见的给法是每隔5%或10%半径给一批离散点到我程序里我先做线性插值把离散点映射到每个叶素中心。插值不用取高阶多项式线性插值已经足够因为BEMT的最终结果对几何的局部小误差不敏感但对总体的弦长和扭角水平很敏感。如果手里只有少量几个测量点比如只有25%、50%、75%三处的弦长和扭角也不用慌linear interp1即可。举个例子某两叶小桨半径0.5 m弦长从根部的0.08 m线性变到尖部的0.04 m扭角从35°0.611 rad变到0°Matlab里写r_nodes [0.075 0.15 0.25 0.375 0.5]; % 半径节点 c_nodes [0.08 0.075 0.065 0.05 0.04]; % 弦长节点 beta_nodes_deg [35 30 22 12 0]; % 扭角节点度 beta_nodes deg2rad(beta_nodes_deg); c_i interp1(r_nodes, c_nodes, r_i, linear, extrap); beta_i interp1(r_nodes, beta_nodes, r_i, linear, extrap);推荐把几何定义放在主程序顶部集中成单独的小函数或脚本段。我踩过的坑是把几何散落在各个循环里后面想换一把桨、对照实验数据时改起来真的痛苦。几何数据尽量集中管理换桨时不动求解器只改几何部分。2.3 翼型升阻力极曲线的处理比几何更麻烦的是翼型气动数据。BEMT每个叶素都要查当前攻角下的升力系数CL和阻力系数CD这需要一整条“攻角—CL—CD”极曲线。来源可以是试验数据、翼型软件算出来的数据或者公开数据库比如UIUC翼型数据库如果你做风力机NREL的报告里也有现成的数据。数据结构上我提前准备好三列向量——alpha_airfoil、CL_airfoil、CD_airfoil攻角范围至少要覆盖失速前后比如-10°到25°因为低速前进比下攻角很容易冲进15°以上。Matlab里查询用interp1同时要处理攻角超出数据范围的情况。超出上界时我采用“CL保持失速后趋势平缓、CD继续增大”的策略简单写if abs(alpha) max(alpha_airfoil) CL CL_airfoil(end) 0.2 * sign(alpha); CD CD_airfoil(end) 0.5 * (abs(alpha) - max(alpha_airfoil)); else CL interp1(alpha_airfoil, CL_airfoil, alpha, linear); CD interp1(alpha_airfoil, CD_airfoil, alpha, linear); end这个方法不是精确物理只是防止代码在小J大攻角工况下乱蹦。后面我会专门讲失速修正这里先稳住迭代。2.4 为什么忽略桨毂区域可以省很多麻烦桨毂的影响分两块。一是几何上不参与气动升力二是强度上它负责连接轴和桨叶叶片根部一般有过渡倒角气动弦线方向混乱二维翼型假设在那里完全失效。BEMT的每个叶素都假设无限长翼型端部三维效应本来就要靠损失因子修正根部如果硬算反而会在低前进比时算出莫名其妙的巨大推力污染积分结果。所以从0.15R起算既是物理判断也是数值健壮性选择。如果你的螺旋桨根部形状比较特殊比如涵道桨根部有一个很宽的整流罩可以适当把起点提到0.2R或者0.25R。注意每次改起点之后积分区间总长度变了最后算出来的CT、CP也会有一点不同这是正常的不必太纠结因为桨根贡献本身占比不大。3. 诱导速度迭代整套计算的核心逻辑3.1 致动盘模型下的气流速度分解BEMT的每一个叶素都同时面对两股速度一是来流和诱导叠加后的轴向速度二是旋转带来的切向速度。轴向方向上远方来流速度V在到达桨盘时被挡住了实际通过桨盘的速度是V × (1 - a)其中a就是轴向诱导因子。a的含义可以理解为“气流被螺旋桨‘阻慢’的比例”a 0表示没有阻挡a接近1表示气流几乎被完全堵住——这在悬停、低前进比工况下会出现也是迭代最容易翻车的地方。切向方向上桨叶自己以角速度Ω转动同时它把气流“拧”着转了实际相对桨叶的切向速度是Ωr × (1 a′)其中a′是切向诱导因子表示气流被桨叶带动旋转的程度。这样每个叶素的气流三角形就确定了φ atan( (1 - a) × V / ( (1 a′) × Ω × r ) )φ是入流角——当地气流方向与旋转平面之间的夹角。攻角就是扭角减去入流角α β - φ整个BEMT的迭代过程就是围绕这φ和α展开的。3.2 从动量方程和叶素力平衡推导迭代式每个叶素上两套理论给出同一个力的两种表达。叶素理论说叶素产生的轴向推力增量等于动压×弦长×升阻力系数组合dT 0.5 × ρ × W² × B × c × (CL × cosφ - CD × sinφ) × dr其中W是合速度W (1 - a) × V / sinφ。动量理论说气流通过圆环面积2πrdr时轴向动量损失对应推力增量dT 4π × r × ρ × V² × a × (1 - a) × F × dr这里的F是叶尖损失因子后面会讲先当作修正系数。两边相等整理出轴向诱导因子的迭代式a / (1 - a) [B × c × (CL × cosφ - CD × sinφ)] / [8π × r × F × sin²φ]同样扭矩增量从叶素理论出发可以换成切向诱导因子的迭代式a′ / (1 a′) [B × c × (CL × sinφ CD × cosφ)] / [8π × r × F × sinφ × cosφ]注意这两个式子右侧都含φφ又依赖a和a′所以没法直接解必须迭代。3.3 叶尖损失因子F的引入桨叶尖端存在一个物理现象压力面气流会绕过桨尖翻到吸力面产生叶尖涡导致叶片尖部载荷下降而不是像二维翼型那样保持理论值。Prandtl提出的叶尖损失因子用简单公式近似描述这个效应f (B / 2) × (R - r) / (r × sinφ)F (2 / π) × acos( exp(-f) )当r趋近R时f趋近0F趋近0载荷归零在桨叶根部附近F趋近1影响很小。这是一个纯几何修正不算昂贵但如果没有它BEMT会高估整把桨的推力和扭矩尤其是在效率曲线靠近大J端时误差会明显放大。我在实现中发现叶尖损失因子对叶素数量很敏感。如果N太少比如只有10个叶素最外圈的叶素中心点离桨尖太远F会被高估结果就是推力偏高。至少取20个叶素最好30个以上叶尖位置的分辨率才够。下文我会给一组经验参数。3.4 迭代初值与松弛因子的选择迭代式是收敛的但初始猜测和更新策略直接决定收敛速度和解的稳定性。一般把a和a′初值都设为0.1这个值在中等载荷下已经很接近真实解迭代几步就到了。真正的麻烦来自低前进比J0.1此时a可能冲到0.4以上甚至接近0.5。动量理论在a0.5时物理上失效因为滑流速度变为负值——自由流无法维持这么强的阻挡。处理方式有两个我自己都用。一是限制迭代更新步长采用欠松弛a_new (1 - ω) × a_old ω × a_formulaω取0.2到0.5避免a和a′在解附近震荡。另一个是把a硬限制在不超过0.5超过就截断到0.49并且把CL的失速修正也加上双保险。直接上全步长迭代在低J工况能看到a在两个极端值之间来回弹甚至直接溢出变成NaN这点务必提前在代码里堵住。4. 不同前进比性能曲线的Matlab实现4.1 主程序结构几何定义、参数设置、前进比扫描主脚本的功能分解成三块设置参数、循环前进比、调用求解器。参数包括桨半径、桨叶数、转速和来流速度范围。转速n固定比如50转/秒D取1 m那么前进比从0.1到1.0对应来流速度从5 m/s到50 m/s。如果只想看巡航段可以把J的范围缩窄但画全曲线看趋势更有价值。主程序框架大致长这样%% 参数与工况设置 R 0.5; % 桨叶半径, m B 2; % 桨叶数 n_rps 50; % 转速, 转/秒 Omega 2 * pi * n_rps; rho 1.225; D 2 * R; N 30; % 叶素数量 N_iter 80; % 迭代上限 tol 1e-6; % 收敛容差 omega_relax 0.3; % 松弛因子 % 前进比扫描范围 J_list 0.1:0.05:1.0; V_list J_list * n_rps * D; %% 预分配输出 CT_list zeros(size(J_list)); CP_list zeros(size(J_list)); eta_list zeros(size(J_list)); %% 主循环 for k 1:length(J_list) Vinf V_list(k); [CT_list(k), CP_list(k)] bem_solver(...); eta_list(k) J_list(k) * CT_list(k) / CP_list(k); end这里bem_solver是我写的单个工况求解函数把半径分段、几何插值、翼型查询、迭代和积分都封装进去。注意不要把诱导迭代写进主循环那样J循环里嵌一个30叶素的迭代代码又乱又慢调试也不方便。4.2 单个工况下的BEM求解器代码这是全文最核心的代码块。每个叶素从初始a、a′出发计算φ和α查翼型极曲线更新a、a′迭代直至收敛或达到上限。下面给一个我测试过的简化版保证能跑通function [CT, CP] bem_solver(R, B, Omega, Vinf, rho, ... r_nodes, c_nodes, beta_nodes, ... alpha_airfoil, CL_airfoil, CD_airfoil, ... N, N_iter, tol, omega_relax) r_vec linspace(0.15*R, R, N); dr r_vec(2) - r_vec(1); c_i interp1(r_nodes, c_nodes, r_vec, linear, extrap); beta_i interp1(r_nodes, beta_nodes, r_vec, linear, extrap); dT_total 0; dQ_total 0; for i 1:N r r_vec(i); c c_i(i); beta beta_i(i); % 初值 a 0.1; a_prime 0.1; for iter 1:N_iter phi atan2((1 - a) * Vinf, (1 a_prime) * Omega * r); alpha beta - phi; % 翼型数据查询与失速粗略修正 if abs(alpha) max(alpha_airfoil) CL CL_airfoil(end) 0.2 * sign(alpha); CD CD_airfoil(end) 0.5 * (abs(alpha) - max(alpha_airfoil)); else CL interp1(alpha_airfoil, CL_airfoil, alpha, linear); CD interp1(alpha_airfoil, CD_airfoil, alpha, linear); end % 叶尖损失 f (B / 2) * (R - r) / (r * abs(sin(phi)) 1e-10); F (2 / pi) * acos(exp(-f)); if F 0.01 F 0.01; end % 动量—叶素联立 A1 (B * c * (CL * cos(phi) - CD * sin(phi))) / ... (8 * pi * r * F * sin(phi)^2 1e-10); a_new A1 / (1 A1); A2 (B * c * (CL * sin(phi) CD * cos(phi))) / ... (8 * pi * r * F * sin(phi) * cos(phi) 1e-10); a_prime_new A2 / (1 - A2); % 限制范围 a_new min(a_new, 0.49); a_prime_new max(a_prime_new, -0.2); % 松弛更新 a (1 - omega_relax) * a omega_relax * a_new; a_prime (1 - omega_relax) * a_prime omega_relax * a_prime_new; if abs(a_new - a) tol abs(a_prime_new - a_prime) tol break; end end % 最终计算该叶素载荷 phi atan2((1 - a) * Vinf, (1 a_prime) * Omega * r); alpha beta - phi; W (1 - a) * Vinf / sin(phi); if abs(alpha) max(alpha_airfoil) CL CL_airfoil(end) 0.2 * sign(alpha); CD CD_airfoil(end) 0.5 * (abs(alpha) - max(alpha_airfoil)); else CL interp1(alpha_airfoil, CL_airfoil, alpha, linear); CD interp1(alpha_airfoil, CD_airfoil, alpha, linear); end dT 0.5 * rho * W^2 * B * c * (CL * cos(phi) - CD * sin(phi)) * dr; dQ 0.5 * rho * W^2 * B * c * (CL * sin(phi) CD * cos(phi)) * r * dr; dT_total dT_total dT; dQ_total dQ_total dQ; end CT dT_total / (rho * n_rps^2 * D^4); CP (Omega * dQ_total) / (rho * n_rps^3 * D^5); end注意几个细节。一是sin(phi)在接近零时会导致分母爆炸我在分母上加小量1e-10并做截断二是切向诱导因子的迭代表达式里有个“1 - A2”的分母如果A2接近1a′也会爆所以我把a′下限定在-0.2防止倒拖负值。三是力矩系数的分母里用了n_rps而不是Ω这是无量纲定义的要求——如果你不熟悉推力系数和功率系数的标准定义很容易搞混转速的量纲。最稳妥的做法是先理清楚CT、CP定义的参考转速到底是转/秒还是弧度/秒不同文献里定义有差异我这一版统一采用n_rps作为参考跟常见的桨性能数据表口径一致。4.3 积分与无量纲系数上面代码里我已经把推力、扭矩累加并除以无量纲基准得出CT和CP。这里特别提一下效率的分子分母单位问题。η J × CT / CP三个数都是无量纲的不会出错。但如果你中途想验证“悬停效率”也就是J 0附近的功率载荷要注意η公式此时分子为0效率趋近于0这是物理事实——原地不动推力做的有用功为零。输出曲线时我习惯同时画三张图CT-J、CP-J、η-J。效率曲线通常有一个明显的峰值峰值对应的J就是该桨设计点附近的前进比。如果峰值过于尖锐、两侧下降很快说明桨叶扭角分布偏“工况专一”在宽速域工作时效率表现会差一些如果峰值平坦说明这是一把适应范围宽但可能单点效率不冒尖的桨。这些判断能从曲线上直接读出来。4.4 完整运行流程和典型耗时运行主脚本后在命令行敲一遍定义几何、设置转速、循环J、调用bem_solver、画图。实测下来的耗时在普通笔记本上30个叶素、20个前进比一共600次叶素迭代每次最多80次内迭代总时间不到5秒。这个速度意味着你可以把BEM求解器嵌进优化循环用遗传算法或粒子群算法去找最优扭角分布每次适应度评估就是几秒的事这是CFD完全做不到的。输出示例上一颗中等扭转的桨大概会给出类似这样的趋势J 0.1时CT最大接近0.15甚至更高CP在J 0.2到0.3附近达到峰值效率曲线在J 0.5附近出现峰值可能落在0.65到0.75之间。这些数值范围可以帮助你自我校验如果算出来CT高达几甚至为负或者效率超过1大概率是代码里哪一步单位或定义出了问题不是物理现象。5. 算例验证用公开几何数据跑一轮5.1 测试桨参数与输入数据为了验证代码我用了一组类似小型无人机两叶桨的参数半径0.5 m两叶转速50转/秒弦长从根部0.08 m线性降到尖部0.04 m扭角从根部35°降到尖部0°翼型数据取某典型低速翼型的极曲线攻角范围-5°到20°最大CL约1.3失速攻角约12°CD从0.008平缓增加到0.02。假设你已经把这些数据准备好直接替换我上面代码里的数组就能跑。如果没有现成攻角数据先随便按“线性增长失速平台”填充一组做占位验证代码链路没问题后再换真实数据。5.2 性能趋势是否符合物理规律我跑完的典型结果如下可以拿来对一下你算出来的曲线规律前进比J推力系数CT功率系数CP效率η0.100.1480.0950.1560.250.1210.0880.3440.400.0920.0760.4840.550.0620.0610.5590.700.0340.0430.5530.850.0100.0260.327推力系数单调下降符合“来流越快、桨叶攻角越小、推力越弱”的直觉效率先升后降峰值在J约0.6附近这个位置跟叶素平均攻角接近零升角附近有关。整体趋势和参考数据吻合说明代码逻辑没有大的方向性错误。5.3 和简单动量理论对比如果忽略几何分布只用致动盘动量理论算会得到一版理想化的CT和效率曲线——它通常预测的效率比BEMT高不少因为动量理论假设桨盘载荷均匀、没有型阻。BEMT的价值恰恰在于把型阻和几何分布的影响吃进去了算出来的效率和推力都比“理论极限”低一些这才是真实桨叶的表现。注意如果BEMT算出来的效率明显高于同工况动量理论结果或者效率曲线出现不合理的锯齿先检查塔尖损失因子有没有生效再检查攻角查询时有没有插值出界。这是我调试时最容易出的两类问题。6. 收敛失败与失速工况的实战处理6.1 大推力工况发散的原因低前进比、也就是高推力工况是整个BEMT算法最容易崩的地方。原因在动量理论上限轴向诱导因子a不能无限制增大a超过0.5以后动量方程的解进入非物理区间迭代式右侧变成负值a直接跳成负数然后φ变成虚的或反号攻角瞬间乱套。表现就是代码输出NaN或者CT突然变成负的。我在前面代码里用了一个硬截断a_new限制在0.49这是治疗“a疯狂膨胀”的对症药。但要注意截断不能滥用——如果a被持续性压在上限说明该工况你大概率进入了BEMT适用域的边缘这时算出的推力值可以参考趋势但不要当成精确值对外输出。6.2 失速修正策略失速是BEMT另一个高频藏雷点。小J工况下叶素攻角很容易超过翼型的失速攻角比如12°到15°。二维翼型数据在失速后测得的分散性很大如果直接外推CL的走势会把迭代带偏。我采用的策略是分段修正攻角还没到失速角时正常查表超过了就按小斜率继续增长但要抑制CL的继续飙升同时让CD随攻角差快速增长。这样处理虽然不精确但能保证积分结果平滑不至于在J扫描时搞出锯齿。如果想要更精细的失速后估计可以用Viterna-Corrigan模型把失速后的CL、CD按平板理论处理这个方法在风力机行业用得很多。不过对螺旋桨初步性能分析我的经验是差不太多先保住平滑和收敛等锁定工况后再用CFD校核失速后的部分。6.3 叶素数量、迭代上限和收敛容差的经验值这些参数我踩过不少坑直接给一组我实际调过的经验值叶素数量N取30。少于15个叶尖损失因子的分辨率不够多于60个不仅没有明显精度提升反而容易在一些半径点上出现局部不收敛。迭代上限N_iter设80。正常工况20步以内就能收敛设80是为了在低J工况下留足余量但不至于无限等下去。收敛容差tol设1e-6比较稳。再收紧到1e-8对最终推力影响小于0.1%徒增耗时。松弛因子ω取0.3。小于0.1收敛太慢需要上百步大于0.6在低J工况容易震荡。用0.3这个值各个工况都能平稳进解。6.4 如果换成变转速扫描要注意什么标题强调恒定转速但很多人实际会想改转速扫。如果你改成固定来流、变转速物理上没问题但要知道转速变了桨叶当地的雷诺数也变了翼型的CL、CD极曲线随之变化严格来说不能沿用同一张极曲线。如果只是粗略看看趋势影响不大如果要做精细设计就必须准备多张不同雷诺数下的极曲线然后在迭代时按当地雷诺数插值选择。这是个有点“重”的扩展我建议先跑通恒定转速版本再考虑是否需要引入雷诺数修正。另外变转速扫描的结果画出来横轴J和固定转速扫描的J虽然定义一样但同一把桨在不同转速下由于雷诺数差异CT、CP会有百分之几到十几的区别直接对比时要小心。固定转速扫描因为雷诺数变化相对小做桨型横向对比更公平。在跑完整套流程之后我个人最大的感受是BEMT的代码实现本身并不算难真正的门槛在于对每个物理假设的边界心里有数——什么时候该信它什么时候该怀疑它。像我习惯在输出图的旁边额外画一条每叶素的攻角分布曲线如果看到大范围叶素都落在失速后我就知道当前工况的结果只能定性参考不能拿去跟实验数据硬碰。这种判断力比代码本身值钱得多。
返回列表