簡(jiǎn)介基于動(dòng)量葉素理論BEM的風(fēng)力機(jī)葉片設(shè)計(jì)程序面向風(fēng)能工程領(lǐng)域的工程師與研究者可用于計(jì)算軸向和周向誘導(dǎo)因子為風(fēng)力機(jī)功率曲線(xiàn)、剪切層交互及湍流特性分析奠定基礎(chǔ)。該程序以MATLAB腳本文件.m格式形式提供壓縮包內(nèi)僅含1個(gè)文件整體僅3KB通過(guò)輸入葉片幾何參數(shù)、翼型特性和流場(chǎng)條件即可調(diào)用動(dòng)量葉素理論算法快速求出關(guān)鍵誘導(dǎo)因子從而輔助完成葉片氣動(dòng)外形設(shè)計(jì)與性能分析。軸向誘導(dǎo)因子反映葉片對(duì)來(lái)流風(fēng)速的阻滯程度周向誘導(dǎo)因子則體現(xiàn)旋轉(zhuǎn)扭矩效應(yīng)二者對(duì)優(yōu)化發(fā)電效率具有重要意義。資源包輕量緊湊適合用于教學(xué)演示、課題研究及工程初算也可借助云計(jì)算平臺(tái)進(jìn)行大規(guī)模多工況模擬目前已有746人學(xué)習(xí)下載是風(fēng)能領(lǐng)域可快速上手的實(shí)用工具。1. 這個(gè) BEM.m 到底能幫你算出什么不止是誘導(dǎo)因子做風(fēng)力機(jī)葉片設(shè)計(jì)的人對(duì)動(dòng)量葉素理論Blade Element Momentum TheoryBEM這個(gè)名字不會(huì)陌生。但真正動(dòng)手寫(xiě)過(guò)一次 BEM 程序的人都知道這個(gè)理論看似簡(jiǎn)單實(shí)際迭代起來(lái)全是細(xì)節(jié)軸向誘導(dǎo)因子 a 和周向誘導(dǎo)因子 a 的耦合解算、葉尖損失修正、失速狀態(tài)的發(fā)散處理每一步都能讓計(jì)算結(jié)果偏離實(shí)驗(yàn)值一大截。這份壓縮包里的 BEM.m 是一個(gè)用 MATLAB 寫(xiě)的單文件程序核心目標(biāo)就是把葉片沿展向離散成若干葉素逐一求解每個(gè)截面上的 a 和 a最終拿到整臺(tái)風(fēng)力機(jī)的氣動(dòng)性能和載荷分布。適合正在做葉片氣動(dòng)設(shè)計(jì)、課程設(shè)計(jì)、或者想驗(yàn)證自己手寫(xiě) BEM 代碼正確性的工程師和學(xué)生。2. 動(dòng)量葉素理論的原理與程序結(jié)構(gòu)先把 a 和 a 的耦合關(guān)系講透2.1 什么是軸向誘導(dǎo)因子和周向誘導(dǎo)因子在風(fēng)力機(jī)運(yùn)行時(shí)來(lái)流經(jīng)過(guò)旋轉(zhuǎn)的葉輪后速度會(huì)發(fā)生兩個(gè)方向的變化沿軸向的速度分量被葉片吸收了一部分導(dǎo)致通過(guò)葉輪平面的風(fēng)速低于自由來(lái)流周向方向則因?yàn)槿~片對(duì)氣流施加切向力讓尾流帶上旋轉(zhuǎn)分量。直接看公式更清楚。取輪盤(pán)半徑 r 處的一個(gè)環(huán)形截面軸向誘導(dǎo)因子 a 的定義是a 1 - V_a / V_∞其中 V_a 是葉輪平面處的軸向速度V_∞ 是自由來(lái)流速度。a 的值越大說(shuō)明氣流經(jīng)過(guò)葉輪時(shí)減速越多。周向誘導(dǎo)因子 a 的定義則是a V_θ / (Ω · r)V_θ 是尾流在葉輪平面處的切向速度Ω·r 是當(dāng)?shù)厝~片的切向速度。a 描述的是氣流被葉片帶動(dòng)旋轉(zhuǎn)的程度。這兩個(gè)因子不是獨(dú)立存在的。葉片吸收能量產(chǎn)生推力推力的大小反過(guò)來(lái)影響軸向速度的衰減葉片產(chǎn)生扭矩扭矩又反過(guò)來(lái)誘導(dǎo)出周向流動(dòng)。所以 BEM 方法的本質(zhì)就是把動(dòng)量理論和葉素理論聯(lián)立成方程組解出每個(gè)截面上的 a 和 a再用它們計(jì)算力系數(shù)、功率系數(shù)和載荷。2.2 BEM.m 的整體程序框架從壓縮包內(nèi)的文件結(jié)構(gòu)來(lái)看BEM.m 是一個(gè)單文件腳本。這類(lèi)程序的典型結(jié)構(gòu)通常分為四段參數(shù)輸入段定義葉片幾何各截面的弦長(zhǎng)、扭轉(zhuǎn)角、安裝半徑、翼型升阻力特性表、來(lái)流條件和轉(zhuǎn)速。離散初始化段把葉片沿展向分成 N 個(gè)葉素通常 20 到 40 個(gè)截面就足夠。迭代求解段對(duì)每個(gè)葉素單獨(dú)計(jì)算 a 和 a采用某種迭代格式反復(fù)修正直到相鄰兩次結(jié)果誤差小于設(shè)定閾值。結(jié)果匯總段把所有截面的推力和扭矩累加求出整機(jī)功率系數(shù)、推力系數(shù)同時(shí)輸出每個(gè)截面的局部攻角、升力系數(shù)、阻力系數(shù)和載荷分布。理解這個(gè)框架后再打開(kāi)代碼你就能清楚知道每個(gè)變量對(duì)應(yīng)的物理含義。如果程序里沒(méi)做葉尖修正大概率用的是簡(jiǎn)化模型適合初學(xué)驗(yàn)證但不適合工程計(jì)算。2.3 BEM 方程組是怎么閉合的單個(gè)葉素的平衡方程組包括兩個(gè)核心方程。動(dòng)量理論給出的軸向力平衡關(guān)系是dT 4·π·ρ·V_∞2·a·(1-a)·r·dr葉素理論給出的同一截面上的推力表達(dá)式是dT 0.5·ρ·W2·(C_l·cos φ C_d·sin φ)·c·dr其中 W 是相對(duì)合速度φ 是入流角c 是弦長(zhǎng)C_l 和 C_d 是升力系數(shù)和阻力系數(shù)這兩個(gè)系數(shù)由翼型數(shù)據(jù)表和當(dāng)?shù)毓ソ遣榈谩Ｂ?lián)立兩個(gè)方程就得到了軸向誘導(dǎo)因子的迭代式。周向方向同理動(dòng)量理論給出的扭矩方程是dQ 4·π·ρ·V_∞·(1-a)·a·Ω·r3·dr葉素理論給出的扭矩方程是dQ 0.5·ρ·W2·(C_l·sin φ - C_d·cos φ)·c·r·dr聯(lián)立后得到一個(gè)包含 a 和 a 的非線(xiàn)性方程組。MATLAB 里常見(jiàn)做法是先假設(shè)初始值 a0、a0然后按公式反復(fù)更新直到兩次迭代的相對(duì)偏差小于 1e-5 或者 1e-6。3. BEM.m 代碼實(shí)操?gòu)膮?shù)輸入到結(jié)果輸出的完整流程3.1 翼型數(shù)據(jù)與葉片幾何參數(shù)的輸入拿到 BEM.m 后第一件事不是急著運(yùn)行而是把里面的輸入?yún)?shù)和自己手頭的葉片數(shù)據(jù)對(duì)齊。一個(gè)典型的葉素截面輸入段長(zhǎng)這樣% 葉片幾何參數(shù)r為葉素半徑chord為弦長(zhǎng)twist為扭角 r [2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0]; % 單位m chord [1.20 1.05 0.92 0.80 0.70 0.60 0.52 0.45 0.40]; % 弦長(zhǎng)遞減 twist [25.0 18.0 12.0 8.0 5.0 3.0 1.5 0.5 -1.0]; % 扭角單位度翼型特性表通常以矩陣形式存儲(chǔ)第一列是攻角第二列是升力系數(shù)第三列是阻力系數(shù)% 翼型數(shù)據(jù)表攻角 / Cl / Cd airfoil_data [ -6.0 -0.30 0.012 -4.0 -0.15 0.010 -2.0 0.05 0.008 0.0 0.25 0.007 2.0 0.45 0.008 4.0 0.65 0.009 6.0 0.85 0.011 8.0 1.00 0.015 10.0 1.10 0.022 12.0 1.15 0.035 14.0 1.10 0.060 16.0 1.00 0.100 ];這里的攻角輸入常用度數(shù)制但公式內(nèi)部計(jì)算三角函數(shù)時(shí)必須統(tǒng)一轉(zhuǎn)為弧度否則查表算入流角時(shí)會(huì)出現(xiàn)嚴(yán)重偏差。我建議在代碼開(kāi)頭就寫(xiě)一行twist deg2rad(twist);把所有角度在一開(kāi)始統(tǒng)一轉(zhuǎn)成弧度避免后續(xù)反復(fù)切換導(dǎo)致踩坑。拿到幾何參數(shù)后程序會(huì)計(jì)算每個(gè)葉素的葉尖速比分布。葉尖速比 λ 是設(shè)計(jì)工況下的核心無(wú)量綱參數(shù)定義為葉片尖端線(xiàn)速度除以來(lái)流速度% 運(yùn)行工況 V_inf 10; % 自由來(lái)流速度m/s R 10; % 風(fēng)輪半徑m Omega 2.5; % 轉(zhuǎn)速rad/s與葉尖速比匹配 lambda Omega * R / V_inf; % 葉尖速比常規(guī)水平軸風(fēng)機(jī)在6-9之間這個(gè) λ 值決定了整個(gè)葉片的攻角分布是否合理。如果 λ 太低葉片大部分截面的攻角會(huì)偏大容易進(jìn)入失速區(qū)λ 太高則攻角偏小升力沒(méi)做起來(lái)。通常設(shè)計(jì)時(shí)先定額定風(fēng)速和額定轉(zhuǎn)速反推出 λ再反推每段弦長(zhǎng)和扭角。3.2 誘導(dǎo)因子的迭代求解核心段BEM.m 的核心循環(huán)在代碼中段。每個(gè)葉素單獨(dú)計(jì)算內(nèi)層迭代 a 和 a% 對(duì)每個(gè)葉素截面進(jìn)行迭代求解 for i 1:length(r) a 0; % 軸向誘導(dǎo)因子初值 a_prime 0; % 周向誘導(dǎo)因子初值 sigma B * chord(i) / (2 * pi * r(i)); % 當(dāng)?shù)貙?shí)度 for iter 1:200 % 最多迭代兩百次 phi atan2(V_inf * (1 - a), Omega * r(i) * (1 a_prime)); % 入流角 alpha phi - twist(i); % 當(dāng)?shù)毓ソ?% 查表獲得升阻力系數(shù) Cl interp1(airfoil_data(:,1), airfoil_data(:,2), rad2deg(alpha), linear, extrap); Cd interp1(airfoil_data(:,1), airfoil_data(:,3), rad2deg(alpha), linear, extrap); % BEM 方程組 C_T Cl * cos(phi) Cd * sin(phi); C_N Cl * sin(phi) - Cd * cos(phi); % 軸向因子更新 a_new 1 / (4 * sin(phi)^2 / (sigma * C_T) 1); % 周向因子更新 a_prime_new 1 / (4 * sin(phi) * cos(phi) / (sigma * C_N) - 1); % 收斂判斷 if abs(a_new - a) 1e-5 abs(a_prime_new - a_prime) 1e-5 a a_new; a_prime a_prime_new; break; end a a_new; a_prime a_prime_new; end % 保存當(dāng)前截面的結(jié)果 result_a(i) a; result_ap(i) a_prime; result_alpha(i) rad2deg(alpha); end這段代碼包含了 BEM 方法最核心的收斂邏輯。有幾個(gè)參數(shù)需要特別注意sigma是當(dāng)?shù)貙?shí)度表示該半徑處葉片弦長(zhǎng)占圓周的比例實(shí)度越高葉片對(duì)氣流的阻礙越強(qiáng)誘導(dǎo)因子越大。phi的計(jì)算用atan2而不是atan因?yàn)槿肓鹘强赡茉谌我庀笙轪tan2能正確處理正負(fù)值。阻尼處理如果直接令a a_new在高實(shí)度或大攻角下容易震蕩。工程中常用欠松弛迭代即a 0.5 * (a a_new)犧牲一點(diǎn)收斂速度換取穩(wěn)定性。上面代碼沒(méi)有寫(xiě)這一步但你在實(shí)際跑數(shù)據(jù)時(shí)遇到震蕩可以自行加上。迭代收斂的核心機(jī)制很簡(jiǎn)單先用當(dāng)前 a 和 a 計(jì)算入流角和攻角查表得 Cl、Cd再代入 BEM 方程組得到一組新的 a 和 a循環(huán)往復(fù)。物理含義是葉素受力和動(dòng)量守恒最終達(dá)到自洽狀態(tài)即葉片受到的力和氣流動(dòng)量變化完全平衡。3.3 葉尖損失修正與高誘導(dǎo)因子修正如果沒(méi)有做任何修正直接運(yùn)行上面那段代碼在小葉尖速比工況下算出來(lái)的功率系數(shù)會(huì)明顯偏高因?yàn)槿~片尖端附近的渦脫落導(dǎo)致實(shí)際載荷低于理想值。工程上最常用的是 Prandtl 葉尖修正因子% Prandtl 葉尖修正 F_tip (2 / pi) * acos(exp(- (B / 2) * (R - r(i)) / (r(i) * sin(phi)))); F_hub (2 / pi) * acos(exp(- (B / 2) * (r(i) - R_hub) / (r(i) * sin(phi)))); F F_tip * F_hub;修正后的軸向力平衡式變成a_new 1 / (4 * sin(phi)^2 / (sigma * C_T) 1); % 未修正 a_new_F 1 / (4 * F * sin(phi)^2 / (sigma * C_T) 1); % 加入修正注意加入修正后只是分母多乘了一個(gè) F物理含義是葉尖附近有效損失增大需要更強(qiáng)的動(dòng)量差來(lái)維持平衡。如果不加修正誘導(dǎo)因子會(huì)偏小計(jì)算功率偏高加了修正后更接近真實(shí)葉片的表現(xiàn)。對(duì)于軸向誘導(dǎo)因子超過(guò) 0.4 的情況普通 BEM 方程已經(jīng)失效因?yàn)閯?dòng)量理論假設(shè)尾流是均勻的而實(shí)際已經(jīng)處于湍流尾流狀態(tài)。此時(shí)需要切換到 Glauert 修正經(jīng)驗(yàn)公式。判斷邏輯如下if a 0.4 % Glauert 修正 a_new 2 * a - a^2; % 近似處理具體系數(shù)查文獻(xiàn) a a_new; end這個(gè)修正如果不做程序在低葉尖速比下算出來(lái)的推力和功率都是錯(cuò)的而且錯(cuò)誤量不是百分之幾而是百分之幾十。3.4 整體功率和推力匯總每個(gè)截面算完后把所有葉素的貢獻(xiàn)沿展向積分。數(shù)值上就是累加每個(gè)環(huán)形截面上的推力和扭矩dT(i) 4 * pi * rho * V_inf^2 * a * (1 - a) * r(i) * dr(i); dQ(i) 4 * pi * rho * V_inf * (1 - a) * a_prime * Omega * r(i)^3 * dr(i); T_total sum(dT); % 總推力 Q_total sum(dQ); % 總扭矩 P_total Q_total * Omega; % 軸功率這里dr(i)是葉素的寬度如果葉素分割不均勻建議在初始化時(shí)計(jì)算每個(gè)截面的實(shí)際寬度而不是直接取平均。功率系數(shù)和推力系數(shù)分別為Cp P_total / (0.5 * rho * V_inf^3 * pi * R^2); Ct T_total / (0.5 * rho * V_inf^2 * pi * R^2);Cp 值能直接看出設(shè)計(jì)是否合理。理論極限是 16/27也就是 0.593實(shí)際風(fēng)力機(jī)最高跑到 0.5 左右。如果程序算出來(lái) Cp 超過(guò) 0.593不用懷疑你的修正沒(méi)做對(duì)或者插值出錯(cuò)了。4. 風(fēng)力機(jī)葉片設(shè)計(jì)里的常見(jiàn)避坑記錄BEM.m 運(yùn)行翻車(chē)與排除4.1 誘導(dǎo)因子迭代發(fā)散報(bào)錯(cuò)或者算成負(fù)數(shù)現(xiàn)象程序運(yùn)行后某個(gè)葉素的 a 值突然跳到負(fù)值或者直接 NaN檢查發(fā)現(xiàn)是迭代那一步求解a_new時(shí)分母變成了零或接近零。原因當(dāng)入流角 φ 非常小時(shí)sin(phi)幾乎為零此時(shí) BEM 方程中的分母項(xiàng)4 * sin(phi)^2 / (sigma * C_T) 1趨于無(wú)窮a 反而趨近于零。但更常見(jiàn)的是在接近輪轂的截面上葉片實(shí)度很高初始攻角猜得偏大Cl 查出來(lái)是失速區(qū)的負(fù)斜率段導(dǎo)致C_T出現(xiàn)負(fù)值a_new隨之變負(fù)。解決第一給迭代設(shè)置上下邊界強(qiáng)制 a 在 0 到 1 之間a 在 -1 到 1 之間第二在查表時(shí)對(duì)攻角做限幅處理超出翼型數(shù)據(jù)范圍時(shí)用端點(diǎn)值外插而不是直接報(bào)錯(cuò)。我一般會(huì)寫(xiě)成alpha_clamped max(min(alpha_deg, max_alpha), min_alpha);這樣可以保證 Cl 和 Cd 永遠(yuǎn)落在合理區(qū)間。第三如果某個(gè)截面始終不收斂檢查是不是弦長(zhǎng)和扭角在該處有跳變幾何參數(shù)不連續(xù)會(huì)導(dǎo)致數(shù)值振蕩。4.2 不收斂導(dǎo)致計(jì)算風(fēng)速下 Cp 出現(xiàn)異常峰值現(xiàn)象在某個(gè)來(lái)流速度下Cp 突然飆升到 0.6 以上已經(jīng)突破貝茨極限。原因這是最典型的 BEM 數(shù)值問(wèn)題通常出現(xiàn)在軸向誘導(dǎo)因子 a 大于 0.4 的工況也就是重載荷狀態(tài)下。普通動(dòng)量方程假設(shè)尾流速度為正但當(dāng)葉片吸能過(guò)多時(shí)實(shí)際尾流可能反向方程失真導(dǎo)致 a 被高估。解決啟用 Glauert 修正但要注意修正公式有很多版本。最簡(jiǎn)單實(shí)用的是 BEM 文獻(xiàn)里常見(jiàn)的 Glauert 修正模型if a 0.4 a 0.4 (a - 0.4) * 0.6; % 阻尼修正 end更精確的做法是使用 Spera 修正公式% 定義 Glauert 修正系數(shù) if a 1/3 k 4 * sin(phi)^2 / (sigma * C_T); a_new (1 k * (1 - 2 * a)) / (1 k * (1 - a)); end這個(gè)修正要在每次迭代更新 a 時(shí)即時(shí)生效不能隔一次迭代再修否則耦合不上仍然會(huì)發(fā)散。4.3 葉尖損失修正過(guò)度功率計(jì)算結(jié)果整體偏低現(xiàn)象加上 Prandtl 葉尖修正后中小風(fēng)速段功率計(jì)算比實(shí)驗(yàn)值低了 15% 以上。原因Prandtl 修正原本針對(duì)理想螺旋槳尾流推導(dǎo)直接用在風(fēng)力機(jī)上如果葉尖速比比較低修正因子 F 在葉尖附近趨近于零導(dǎo)致該段幾乎不貢獻(xiàn)載荷。但實(shí)際葉片葉尖雖然載荷降低并不是完全沒(méi)有貢獻(xiàn)。我常見(jiàn)到有人在所有截面上統(tǒng)一乘了 F導(dǎo)致葉片中段被過(guò)度削弱。解決只對(duì)靠近葉尖的最后 10% 到 15% 展向區(qū)域啟用完整的 Prandtl 修正葉片中段保留未修正的計(jì)算值。實(shí)際做法是設(shè)定一個(gè)修正起始半徑比如r_start 0.85 * R在這個(gè)半徑以?xún)?nèi)設(shè)定 F1以外才逐步過(guò)渡到 Prandtl 修正值。4.4 翼型數(shù)據(jù)插值精度不夠?qū)е鹿ソ峭蛔儠r(shí)計(jì)算值跳變現(xiàn)象Cp 曲線(xiàn)在小范圍內(nèi)劇烈波動(dòng)看起來(lái)不是連續(xù)曲線(xiàn)而是一節(jié)一節(jié)的臺(tái)階。原因翼型升阻力表攻角間隔太大或者插值方式用了nearest而不是linear。BEM 程序?qū)?Cl 隨攻角的變化很敏感特別是葉尖附近攻角變化幾度功率可能差出 5%。解決用interp1時(shí)顯式指定線(xiàn)性插值并且保證攻角表覆蓋從負(fù)失速到正失速的完整范圍。如果翼型原始數(shù)據(jù)點(diǎn)稀疏在輸入程序前自己先做一次加密插值把步長(zhǎng)細(xì)化到 0.5 度以?xún)?nèi)。alpha_fine -6:0.5:18; Cl_fine interp1(raw_alpha, raw_Cl, alpha_fine, pchip); Cd_fine interp1(raw_alpha, raw_Cd, alpha_fine, pchip);PCHIP 插值能避免線(xiàn)性插值帶來(lái)的尖角也會(huì)讓迭代過(guò)程更平穩(wěn)。4.5 轉(zhuǎn)速和風(fēng)速的單位不統(tǒng)一導(dǎo)致整體結(jié)果偏了一百倍現(xiàn)象算出來(lái)的功率數(shù)值顯然不合理不是太大就是太小而且所有工況下的誤差方向一致。原因MATLAB 默認(rèn)使用國(guó)際單位制但很多翼型數(shù)據(jù)表里攻角用度轉(zhuǎn)速可能從 r/min 直接抄進(jìn)來(lái)沒(méi)轉(zhuǎn)換為 rad/s。如果轉(zhuǎn)速單位寫(xiě)錯(cuò)a 的物理含義就全變了整個(gè)迭代方程崩塌。解決所有輸入?yún)?shù)統(tǒng)一轉(zhuǎn)換為 SI 單位。轉(zhuǎn)速的換算關(guān)系是Omega_rad Omega_rpm * 2 * pi / 60。我在程序開(kāi)頭會(huì)固定寫(xiě)好一行注釋標(biāo)準(zhǔn)塊檢查所有輸入量的單位再運(yùn)行。5. 從 BEM.m 到工程落地結(jié)果驗(yàn)證、參數(shù)調(diào)整與批量工況計(jì)算5.1 用已知算例驗(yàn)證程序正確性拿到 BEM.m 并跑通后第一件事不是改參數(shù)而是驗(yàn)證。找一組文獻(xiàn)中已有實(shí)驗(yàn)數(shù)據(jù)的葉片參數(shù)比如 NREL 5MW 或較小的 NREL Phase VI 葉片數(shù)據(jù)。Phase VI 葉片在 NASA 風(fēng)洞做過(guò)全尺度實(shí)驗(yàn)有不同風(fēng)速下的功率和載荷數(shù)據(jù)非常適合用來(lái)校驗(yàn) BEM 程序。驗(yàn)證方法是固定輸入風(fēng)速?gòu)?5 m/s 到 25 m/s 變化每隔 1 m/s 計(jì)算一個(gè)工況畫(huà)出功率和推力曲線(xiàn)。對(duì)照實(shí)驗(yàn)數(shù)據(jù)的趨勢(shì)額定風(fēng)速以下功率隨風(fēng)速的三次方上升額定風(fēng)速附近功率增速放緩逐漸趨平超過(guò)切出風(fēng)速后功率開(kāi)始下降。如果自算曲線(xiàn)在低風(fēng)速段斜率對(duì)高風(fēng)速段對(duì)不上優(yōu)先檢查失速區(qū)域的翼型數(shù)據(jù)是否準(zhǔn)確而不是懷疑程序。NREL Phase VI 在失速區(qū)的 Cl 系數(shù)受雷諾數(shù)影響很大低雷諾數(shù)下失速提前升力下降更快。5.2 參數(shù)掃描與多工況模擬單工況跑通只是第一步做葉片設(shè)計(jì)時(shí)要看的是整條功率曲線(xiàn)。BEM.m 的向量化程度如果不夠高跑 20 個(gè)風(fēng)速工況可能需要幾十秒?？梢愿膶?xiě)最外層循環(huán)把所有風(fēng)速作為一個(gè)數(shù)組傳入V_list 4:1:25; % 風(fēng)速掃描范圍 Cp_array zeros(size(V_list)); for j 1:length(V_list) V_inf V_list(j); % 在這里調(diào)用 BEM 核心計(jì)算函數(shù)或者直接用內(nèi)層循環(huán)重算 Cp_array(j) compute_BEM(V_inf, blade_geo, airfoil_data); end這種參數(shù)掃描的價(jià)值在于能快速看到葉片設(shè)計(jì)點(diǎn)外的表現(xiàn)。比如設(shè)計(jì)風(fēng)速 10 m/s 下 Cp 是 0.48但在 7 m/s 以下可能掉到 0.35 以下說(shuō)明葉片在低風(fēng)速下匹配不佳。這時(shí)調(diào)整扭角分布或者轉(zhuǎn)速策略比改弦長(zhǎng)來(lái)的更直接。5.3 用程序優(yōu)化葉片設(shè)計(jì)參數(shù)的使用模板BEM.m 不僅能算固定葉片的性能也能反過(guò)來(lái)用于參數(shù)尋優(yōu)。常見(jiàn)做法是保持翼型和弦長(zhǎng)不變掃描扭角分布整體偏移量。扭角整體加 2 度和減 2 度葉片在不同風(fēng)速下的功率表現(xiàn)差別很大是因?yàn)楣ソ钦w平移后葉片運(yùn)行點(diǎn)偏離了翼型最大升阻比對(duì)應(yīng)的攻角失速點(diǎn)也隨之移動(dòng)。這里有一個(gè)我自己的使用習(xí)慣每調(diào)整一輪設(shè)計(jì)參數(shù)保留一份算出來(lái)的攻角分布表。葉片設(shè)計(jì)最大的隱性約束是葉尖攻角不能低于 2 度否則在葉片污染或來(lái)流偏差時(shí)容易提前失速。根部截面攻角應(yīng)該控制在 8 到 12 度太低了葉根段幾乎不產(chǎn)生升力。5.4 程序驗(yàn)證的正確參考基準(zhǔn)如果你的目標(biāo)不是設(shè)計(jì)新葉片而是把 BEM.m 的結(jié)果和 CFD 結(jié)果對(duì)比注意兩者載荷分布的定義一致性。BEM 輸出的每個(gè)葉素的局部推力是沿軸向的而 CFD 常輸出法向力系數(shù) Cn需要換算Cn Cl * cos(phi) Cd * sin(phi); % 法向力系數(shù)這個(gè)系數(shù)歸一化后和 CFD 結(jié)果的對(duì)比才有意義。直接拿 Cl 或 Cd 去對(duì)比因?yàn)閮烧叨x的參考不同怎么比對(duì)都對(duì)不上。最后提醒一句BEM.m 這類(lèi)程序最大的隱患不在代碼邏輯而在輸入數(shù)據(jù)的質(zhì)量。翼型表、弦長(zhǎng)分布、扭角分布每一項(xiàng)都要有據(jù)可查。我跑 BEM 程序這幾年踩過(guò)最大的坑就是拿了一個(gè)沒(méi)有任何來(lái)源標(biāo)注的翼型數(shù)據(jù)表直接算結(jié)果在失速區(qū)查出的 Cl 完全不連續(xù)整個(gè) Cp 曲線(xiàn)都是錯(cuò)的。從那以后我每次跑 BEM 前都強(qiáng)制自己先把翼型數(shù)據(jù)來(lái)源整理成一份帶參考注釋的表格確認(rèn)每個(gè)攻角對(duì)應(yīng)的 Cl、Cd 來(lái)自風(fēng)洞實(shí)驗(yàn)還是 XFOIL 計(jì)算再放進(jìn)程序。這個(gè)習(xí)慣節(jié)省了無(wú)數(shù)排查時(shí)間。你有 BEM.m 在手也要記得自己的程序永遠(yuǎn)只對(duì)你輸入的數(shù)據(jù)負(fù)責(zé)希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取