萬全喜,張明輝,吳家龍
(山東科技大學機械電子工程學院,山東青島 266590)
隨著環(huán)境惡化和能源缺乏等問題的日益突出,各國都把目光聚集在可再生能源上。風能作為一種無污染的可再生綠色能源,越來越引起人們的關注。其中,風力發(fā)電是風能利用的主要形式之一[1]。風力機的核心部分是葉片,其主要功能就是將風能轉化成機械能,它的轉化效率的高低將直接影響并且在很大程度上決定著風力機的性能。因此風力機葉片的優(yōu)化設計非常重要。傳統(tǒng)的Wilson方法是選擇同一的最佳攻角α,這樣設計出的葉片外形并不是最優(yōu),氣動性能也有所降低。作者根據葉片的最佳設計攻角α、升力因數CL沿葉片展長方向呈非線性分布,對傳統(tǒng)的Wilson方法進行改進,利用MATLAB軟件計算葉片外形各參數,優(yōu)化并修正了200 kW風力機葉片的氣動外形。
風力機的葉片氣動外形設計非常重要,設計氣動性能良好的葉片是風力機獲得最大風能利用系數及優(yōu)良經濟效益的關鍵。風力機葉片的優(yōu)化設計主要包括對葉片數B、葉片直徑D、翼型、安裝角θ和弦長C等參數的確定。
風力發(fā)電機風輪直徑,由下面公式確定:
式中:P=200 kW,空氣密度ρ=1.225 kg/m3,設計風速v1=13 m/s,風能利用系數CP=0.4,葉片數B=3,風力機的機電效率 η=η1η2=0.85(η1為傳動效率,η2為發(fā)電機效率)。風輪的直徑D為23.6 m,除去中心轂直徑后,風輪的半徑取R=11.5 m。
翼型的選擇對風力機的風能利用效率十分重要。一種性能良好的翼型應該是在某一攻角范圍內所對應的升力系數CL較高、阻力系數Cd較小;同時所選雷諾數應與風力機實際運行情況的雷諾數接近;此外,還應具有良好的制造加工工藝性。設計風力機時,應該根據不同的設計需要選擇不同的翼型。此次設計選擇NACA2412翼型。利用Profili軟件可以得出NACA2412翼型在不同雷諾數下升阻比隨攻角變化的曲線圖,如圖1所示。
圖1 不同雷諾數下升阻比隨攻角變化的曲線圖
傳統(tǒng)的Wilson方法是選擇同一的最佳攻角α,繼而得到安裝角θ=φ-α,這樣設計出的葉片外形并不是最優(yōu),氣動性能也有所降低。此次設計是根據葉片的最佳設計攻角α、升力因數CL沿葉片展長方向呈非線性分布,對傳統(tǒng)的Wilson方法進行改進。這里在考慮葉尖和輪轂的損失系數的同時 (阻力對氣動外形設計影響不大,在此不加考慮),還研究了風輪在非設計狀態(tài)下的氣動性能[2]。從而推導出氣動外形設計的基本數學模型:
式中:CP為風能利用系數;a為軸向誘導因子;b為周向誘導因子;Ftip為葉尖修正系數;Fhub為葉根修正系數;F為普朗特修正系數;λ0為葉尖速比;B為葉片數目;φ為葉素入流角 (°);θ為葉素安裝角(°);α為葉素攻角 (°);v1為來流風速 (m/s);ν為運動黏度;R為葉輪半徑 (m);r為葉素剖面到風輪中心的距離 (m);C為葉素弦長 (m);CL為葉素升力系數。
(1)根據葉素理論,把葉片沿展長方向平均分成23等份,得到23個葉素;
(2)針對每個葉素,求解以式 (1)為目標函數、以式 (2)— (5)為約束條件的最優(yōu)問題,得到每一葉素的a、b、F,進而計算出式 (6)及S;
(3)根據式 (7)和翼型的性能實驗數據得到各葉素截面雷諾數Re所對應的最佳攻角α及升力系數CL,進而根據式 (8)求得安裝角θ,根據式 (9)求得弦長C;
(4)對求得的葉片弦長C和安裝角θ進行擬合修正,使其滿足結構、加工工藝等方面的要求;
(5)利用Profili軟件可得到翼型坐標原始數據(x0,y0),最后轉變成三維空間坐標 (x,y,z)。
程序流程框圖如圖2所示。
圖2 程序流程框圖
2.2.1 用MATLAB優(yōu)化工具箱求解干涉因子a、b
用MATLAB優(yōu)化工具箱中求解最優(yōu)化問題最為出色的fmincon函數來求解干涉因子a、b。主程序編寫為turbine.m文件,另外創(chuàng)建兩個m文件用來存放目標函數和條件函數,分別命名為objfun.m和confun.m。
通過查相關文獻,找到以下經驗性初值公式[4]:
式中:λ0為葉尖速比;R為風輪半徑 (m);r為某一截面葉根的距離 (m)。
2.2.2 求解雷諾數Re及升力系數CL
通過Profili軟件,可以求得任何雷諾數Re在各攻角下所對應升力系數CL和阻力系數Cd,及其所對應的升阻比L=CL/Cd。升阻比是葉片設計的一個重要性能指標,為了得到最大風能利用系數,就要使升阻比趨于最大值。利用MATLAB優(yōu)化工具箱中的max函數求解各給定雷諾數Re_s下接近最大升阻比Lmax所對應的升力系數CL_g和攻角α_g。利用MATLAB優(yōu)化工具箱中擬合函數polyfit(Re_s,CL_g,4)擬合雷諾數與升力系數,擬合多項式如下:
根據公式 (7)和 (9)推算出:
將公式 (15)和 (16)聯立,利用MATLAB優(yōu)化工具箱中的非線性方程組求解函數fsolve求解方程組,得到各截面最大升阻比Lmax下所對應的Re和CL。
2.2.3 求解最佳攻角α
利用上述方法計算求得的Re,利用MATLAB中線性插值函數 interp1(Re_s,α_g,Re,’*spline’)[5],解得最佳攻角 α。
2.2.4 葉片弦長C和安裝角θ的修正
經過主程序計算獲得了理想葉片弦長C和葉片安裝角θ。可是,此時葉片的弦長C和安裝角θ是沿葉片展長方向呈非線性分布的,并且葉片根部弦長太大。這種形狀的葉片加工困難,用料也不經濟。為解決此問題,需要同時考慮葉片氣動性能及制造兩因素對理想的葉片弦長C和安裝角θ進行修正。葉片從風中獲得功率的75%是由葉片的前半部分汲取的,葉片根部主要起連接和支撐作用,對風能利用系數CP影響較小。因此選取0.5R~0.9R的各截面,利用MATLAB 工具箱中擬合函數 polyfit(r,C,4)[5]對葉片根部弦長修正,使葉片弦長C適當減小。擬合多項式如下:
同樣利用擬合函數polyfit(r,θ,4)對葉片根部安裝角修正,使葉片安裝角θ適當減小。擬合多項式如下:
將0.2R~R各截面的r值代入公式 (17)和 (18),得到修正后的弦長C和安裝角θ。經過多次實驗,采用4次擬合多項式保證了擬合精度,同時也使葉根處弦長C和安裝角θ得到減小。
葉片弦長C、安裝角θ修正前后對比如圖3、圖4所示,各參數計算結果如表1所示。
圖3 修正前后弦長隨半徑變化圖
圖4 修正前后安裝角隨半徑變化圖
表1 部分葉素截面各參數計算結果數據表
續(xù)表1
葉片各截面的空間坐標求解實質上就是圖形變換,即對組成圖形的各坐標點進行變換。基本思路如下:翼型上下弦的原始坐標數據(x0,y0→—) 翼型以氣動中心為原點的二維坐標數據(x1,y1→—) 葉素截面各離散點的空間實際坐標 (x,y,z)。原理示意見圖5。
圖5 翼型坐標變換示意圖
具體步驟如下:
(1)通過 Profili軟件獲取翼型原始數據 (x0,y0)。
(2)實現翼型1到翼型2的變換。計算以翼型氣動中心作為原點、翼弦軸線為x軸的二維坐標(x1,y1)。氣動中心通常在到翼型前緣的距離為0.25~0.35弦長處[7],此處氣動中心的 x軸坐標設為0.25C,y軸的坐標設為0,設氣動中心坐標為 (x2,y2)。
(3)實現翼型2到翼型3的轉變。結合弦長計算各葉素坐標:
(4)實現翼型3到翼型4的轉變。利用下述公式旋轉葉素二維坐標得到三維空間實際坐標 (x,y,z):
z=r
通過Excel辦公軟件可完成各葉素上所有離散點空間實際坐標的計算,并把坐標保存為3列數字的txt文件的格式。
根據上述翼型坐標變換獲得的空間實際坐標(x,y,z),通過SolidWorks軟件繪制葉片實體模型。具體步驟:利用SolidWorks的“通過xyz曲線”命令繪制各葉素輪廓線 — →利用“平面區(qū)域”命令生成葉素平面 — →利用“曲面放樣”命令實現在各葉素面間放樣生成立體圖。
基于表1各參數數據,結合葉柄部分數據,利用SolidWorks軟件繪出修正前后葉片實體模型如圖6、圖7所示。
圖6 修正前葉片實體模型
圖7 修正后葉片實體模型
選用NACA2412翼型,基于MATLAB軟件,結合Wilson優(yōu)化設計方法對葉片各葉素的弦長和安裝角進行了求解,考慮到制造性和經濟性兩方面的因素,又對弦長C和安裝角θ進行了修正;然后利用Excel辦公軟件采用齊次坐標轉換的方法對翼型二維坐標進行轉換,得到空間三維坐標 (x,y,z);最后應用SolidWorks軟件繪出葉片實體模型。此過程簡單、通用性強、易于進一步分析,為葉片和其他復雜實體建模提供了新思路。
【1】王曉蓉,王偉勝,戴慧珠.我國風力發(fā)電現狀和展望[J].中國電力,2004,37(1):81-84.
【2】林閩,張崇巍,張艷紅,等.小型風力發(fā)電機葉輪設計[J].風機技術,2007(1):28-30.
【3】包飛.風力機葉片幾何設計與空氣動力學仿真[D].大連:大連理工大學,2008:26.
【4】王凡.風力發(fā)電機的葉片設計方法研究[D].南京:南京理工大學,2007.
【5】王沫然.MATLAB6.0與科學計算[M].北京:電子工業(yè)出版社,2002.
【6】王學永.風力發(fā)電機葉片設計及三維建模[D].北京:華北電力大學,2008:34-36.
【7】賀德馨.風工程與工業(yè)空氣動力學[M].北京:國防工業(yè)出版社,2006.
【8】陳家權,楊新彥.風力機葉片立體圖設計[J].機電工程,2006,23(4):37-40.