萬全喜，張明輝，吳家龍

(山東科技大學(xué)機(jī)械電子工程學(xué)院，山東青島 266590)

隨著環(huán)境惡化和能源缺乏等問題的日益突出，各國(guó)都把目光聚集在可再生能源上。風(fēng)能作為一種無污染的可再生綠色能源，越來越引起人們的關(guān)注。其中，風(fēng)力發(fā)電是風(fēng)能利用的主要形式之一［1］。風(fēng)力機(jī)的核心部分是葉片，其主要功能就是將風(fēng)能轉(zhuǎn)化成機(jī)械能，它的轉(zhuǎn)化效率的高低將直接影響并且在很大程度上決定著風(fēng)力機(jī)的性能。因此風(fēng)力機(jī)葉片的優(yōu)化設(shè)計(jì)非常重要。傳統(tǒng)的Wilson方法是選擇同一的最佳攻角α，這樣設(shè)計(jì)出的葉片外形并不是最優(yōu)，氣動(dòng)性能也有所降低。作者根據(jù)葉片的最佳設(shè)計(jì)攻角α、升力因數(shù)CL沿葉片展長(zhǎng)方向呈非線性分布，對(duì)傳統(tǒng)的Wilson方法進(jìn)行改進(jìn)，利用MATLAB軟件計(jì)算葉片外形各參數(shù)，優(yōu)化并修正了200 kW風(fēng)力機(jī)葉片的氣動(dòng)外形。

1 風(fēng)力機(jī)葉片氣動(dòng)外形設(shè)計(jì)

風(fēng)力機(jī)的葉片氣動(dòng)外形設(shè)計(jì)非常重要，設(shè)計(jì)氣動(dòng)性能良好的葉片是風(fēng)力機(jī)獲得最大風(fēng)能利用系數(shù)及優(yōu)良經(jīng)濟(jì)效益的關(guān)鍵。風(fēng)力機(jī)葉片的優(yōu)化設(shè)計(jì)主要包括對(duì)葉片數(shù)B、葉片直徑D、翼型、安裝角θ和弦長(zhǎng)C等參數(shù)的確定。

1.1 風(fēng)輪直徑D的確定

風(fēng)力發(fā)電機(jī)風(fēng)輪直徑，由下面公式確定:

式中:P=200 kW，空氣密度ρ=1.225 kg/m3，設(shè)計(jì)風(fēng)速v1=13 m/s，風(fēng)能利用系數(shù)CP=0.4，葉片數(shù)B=3，風(fēng)力機(jī)的機(jī)電效率 η=η1η2=0.85(η1為傳動(dòng)效率，η2為發(fā)電機(jī)效率)。風(fēng)輪的直徑D為23.6 m，除去中心轂直徑后，風(fēng)輪的半徑取R=11.5 m。

1.2 翼型的確定

翼型的選擇對(duì)風(fēng)力機(jī)的風(fēng)能利用效率十分重要。一種性能良好的翼型應(yīng)該是在某一攻角范圍內(nèi)所對(duì)應(yīng)的升力系數(shù)CL較高、阻力系數(shù)Cd較小;同時(shí)所選雷諾數(shù)應(yīng)與風(fēng)力機(jī)實(shí)際運(yùn)行情況的雷諾數(shù)接近;此外，還應(yīng)具有良好的制造加工工藝性。設(shè)計(jì)風(fēng)力機(jī)時(shí)，應(yīng)該根據(jù)不同的設(shè)計(jì)需要選擇不同的翼型。此次設(shè)計(jì)選擇NACA2412翼型。利用Profili軟件可以得出NACA2412翼型在不同雷諾數(shù)下升阻比隨攻角變化的曲線圖，如圖1所示。

圖1 不同雷諾數(shù)下升阻比隨攻角變化的曲線圖

1.3 Wilson優(yōu)化設(shè)計(jì)方法的數(shù)學(xué)模型

傳統(tǒng)的Wilson方法是選擇同一的最佳攻角α，繼而得到安裝角θ=φ－α，這樣設(shè)計(jì)出的葉片外形并不是最優(yōu)，氣動(dòng)性能也有所降低。此次設(shè)計(jì)是根據(jù)葉片的最佳設(shè)計(jì)攻角α、升力因數(shù)CL沿葉片展長(zhǎng)方向呈非線性分布，對(duì)傳統(tǒng)的Wilson方法進(jìn)行改進(jìn)。這里在考慮葉尖和輪轂的損失系數(shù)的同時(shí) (阻力對(duì)氣動(dòng)外形設(shè)計(jì)影響不大，在此不加考慮)，還研究了風(fēng)輪在非設(shè)計(jì)狀態(tài)下的氣動(dòng)性能［2］。從而推導(dǎo)出氣動(dòng)外形設(shè)計(jì)的基本數(shù)學(xué)模型:

式中:CP為風(fēng)能利用系數(shù);a為軸向誘導(dǎo)因子;b為周向誘導(dǎo)因子;Ftip為葉尖修正系數(shù);Fhub為葉根修正系數(shù);F為普朗特修正系數(shù);λ0為葉尖速比;B為葉片數(shù)目;φ為葉素入流角 (°);θ為葉素安裝角(°);α為葉素攻角 (°);v1為來流風(fēng)速 (m/s);ν為運(yùn)動(dòng)黏度;R為葉輪半徑 (m);r為葉素剖面到風(fēng)輪中心的距離 (m);C為葉素弦長(zhǎng) (m);CL為葉素升力系數(shù)。

1.4 Wilson 優(yōu)化設(shè)計(jì)步驟［3］

(1)根據(jù)葉素理論，把葉片沿展長(zhǎng)方向平均分成23等份，得到23個(gè)葉素;

(2)針對(duì)每個(gè)葉素，求解以式 (1)為目標(biāo)函數(shù)、以式 (2)— (5)為約束條件的最優(yōu)問題，得到每一葉素的a、b、F，進(jìn)而計(jì)算出式 (6)及S;

(3)根據(jù)式 (7)和翼型的性能實(shí)驗(yàn)數(shù)據(jù)得到各葉素截面雷諾數(shù)Re所對(duì)應(yīng)的最佳攻角α及升力系數(shù)CL，進(jìn)而根據(jù)式 (8)求得安裝角θ，根據(jù)式 (9)求得弦長(zhǎng)C;

(4)對(duì)求得的葉片弦長(zhǎng)C和安裝角θ進(jìn)行擬合修正，使其滿足結(jié)構(gòu)、加工工藝等方面的要求;

(5)利用Profili軟件可得到翼型坐標(biāo)原始數(shù)據(jù)(x0，y0)，最后轉(zhuǎn)變成三維空間坐標(biāo) (x，y，z)。

2 葉片優(yōu)化設(shè)計(jì)的MATLAB編程

2.1 程序流程

程序流程框圖如圖2所示。

圖2 程序流程框圖

2.2 利用MATLAB編程解決的關(guān)鍵問題

2.2.1 用MATLAB優(yōu)化工具箱求解干涉因子a、b

用MATLAB優(yōu)化工具箱中求解最優(yōu)化問題最為出色的fmincon函數(shù)來求解干涉因子a、b。主程序編寫為turbine.m文件，另外創(chuàng)建兩個(gè)m文件用來存放目標(biāo)函數(shù)和條件函數(shù)，分別命名為objfun.m和confun.m。

通過查相關(guān)文獻(xiàn)，找到以下經(jīng)驗(yàn)性初值公式［4］:

式中:λ0為葉尖速比;R為風(fēng)輪半徑 (m);r為某一截面葉根的距離 (m)。

2.2.2 求解雷諾數(shù)Re及升力系數(shù)CL

通過Profili軟件，可以求得任何雷諾數(shù)Re在各攻角下所對(duì)應(yīng)升力系數(shù)CL和阻力系數(shù)Cd，及其所對(duì)應(yīng)的升阻比L=CL/Cd。升阻比是葉片設(shè)計(jì)的一個(gè)重要性能指標(biāo)，為了得到最大風(fēng)能利用系數(shù)，就要使升阻比趨于最大值。利用MATLAB優(yōu)化工具箱中的max函數(shù)求解各給定雷諾數(shù)Re_s下接近最大升阻比Lmax所對(duì)應(yīng)的升力系數(shù)CL_g和攻角α_g。利用MATLAB優(yōu)化工具箱中擬合函數(shù)polyfit(Re_s，CL_g，4)擬合雷諾數(shù)與升力系數(shù)，擬合多項(xiàng)式如下:

根據(jù)公式 (7)和 (9)推算出:

將公式 (15)和 (16)聯(lián)立，利用MATLAB優(yōu)化工具箱中的非線性方程組求解函數(shù)fsolve求解方程組，得到各截面最大升阻比Lmax下所對(duì)應(yīng)的Re和CL。

2.2.3 求解最佳攻角α

利用上述方法計(jì)算求得的Re，利用MATLAB中線性插值函數(shù) interp1(Re_s，α_g，Re，’*spline’)［5］，解得最佳攻角 α。

2.2.4 葉片弦長(zhǎng)C和安裝角θ的修正

經(jīng)過主程序計(jì)算獲得了理想葉片弦長(zhǎng)C和葉片安裝角θ。可是，此時(shí)葉片的弦長(zhǎng)C和安裝角θ是沿葉片展長(zhǎng)方向呈非線性分布的，并且葉片根部弦長(zhǎng)太大。這種形狀的葉片加工困難，用料也不經(jīng)濟(jì)。為解決此問題，需要同時(shí)考慮葉片氣動(dòng)性能及制造兩因素對(duì)理想的葉片弦長(zhǎng)C和安裝角θ進(jìn)行修正。葉片從風(fēng)中獲得功率的75%是由葉片的前半部分汲取的，葉片根部主要起連接和支撐作用，對(duì)風(fēng)能利用系數(shù)CP影響較小。因此選取0.5R～0.9R的各截面，利用MATLAB 工具箱中擬合函數(shù) polyfit(r，C，4)［5］對(duì)葉片根部弦長(zhǎng)修正，使葉片弦長(zhǎng)C適當(dāng)減小。擬合多項(xiàng)式如下:

同樣利用擬合函數(shù)polyfit(r，θ，4)對(duì)葉片根部安裝角修正，使葉片安裝角θ適當(dāng)減小。擬合多項(xiàng)式如下:

將0.2R～R各截面的r值代入公式 (17)和 (18)，得到修正后的弦長(zhǎng)C和安裝角θ。經(jīng)過多次實(shí)驗(yàn)，采用4次擬合多項(xiàng)式保證了擬合精度，同時(shí)也使葉根處弦長(zhǎng)C和安裝角θ得到減小。

葉片弦長(zhǎng)C、安裝角θ修正前后對(duì)比如圖3、圖4所示，各參數(shù)計(jì)算結(jié)果如表1所示。

圖3 修正前后弦長(zhǎng)隨半徑變化圖

圖4 修正前后安裝角隨半徑變化圖

表1 部分葉素截面各參數(shù)計(jì)算結(jié)果數(shù)據(jù)表

續(xù)表1

3 翼型坐標(biāo)轉(zhuǎn)換［6］

葉片各截面的空間坐標(biāo)求解實(shí)質(zhì)上就是圖形變換，即對(duì)組成圖形的各坐標(biāo)點(diǎn)進(jìn)行變換?；舅悸啡缦?翼型上下弦的原始坐標(biāo)數(shù)據(jù)(x0，y0→—) 翼型以氣動(dòng)中心為原點(diǎn)的二維坐標(biāo)數(shù)據(jù)(x1，y1→—) 葉素截面各離散點(diǎn)的空間實(shí)際坐標(biāo) (x，y，z)。原理示意見圖5。

圖5 翼型坐標(biāo)變換示意圖

具體步驟如下:

(1)通過 Profili軟件獲取翼型原始數(shù)據(jù) (x0，y0)。

(2)實(shí)現(xiàn)翼型1到翼型2的變換。計(jì)算以翼型氣動(dòng)中心作為原點(diǎn)、翼弦軸線為x軸的二維坐標(biāo)(x1，y1)。氣動(dòng)中心通常在到翼型前緣的距離為0.25～0.35弦長(zhǎng)處［7］，此處氣動(dòng)中心的 x軸坐標(biāo)設(shè)為0.25C，y軸的坐標(biāo)設(shè)為0，設(shè)氣動(dòng)中心坐標(biāo)為 (x2，y2)。

(3)實(shí)現(xiàn)翼型2到翼型3的轉(zhuǎn)變。結(jié)合弦長(zhǎng)計(jì)算各葉素坐標(biāo):

(4)實(shí)現(xiàn)翼型3到翼型4的轉(zhuǎn)變。利用下述公式旋轉(zhuǎn)葉素二維坐標(biāo)得到三維空間實(shí)際坐標(biāo) (x，y，z):

z=r

通過Excel辦公軟件可完成各葉素上所有離散點(diǎn)空間實(shí)際坐標(biāo)的計(jì)算，并把坐標(biāo)保存為3列數(shù)字的txt文件的格式。

4 葉片三維建模［8］

根據(jù)上述翼型坐標(biāo)變換獲得的空間實(shí)際坐標(biāo)(x，y，z)，通過SolidWorks軟件繪制葉片實(shí)體模型。具體步驟:利用SolidWorks的“通過xyz曲線”命令繪制各葉素輪廓線 — →利用“平面區(qū)域”命令生成葉素平面 — →利用“曲面放樣”命令實(shí)現(xiàn)在各葉素面間放樣生成立體圖。

基于表1各參數(shù)數(shù)據(jù)，結(jié)合葉柄部分?jǐn)?shù)據(jù)，利用SolidWorks軟件繪出修正前后葉片實(shí)體模型如圖6、圖7所示。

圖6 修正前葉片實(shí)體模型

圖7 修正后葉片實(shí)體模型

5 總結(jié)

選用NACA2412翼型，基于MATLAB軟件，結(jié)合Wilson優(yōu)化設(shè)計(jì)方法對(duì)葉片各葉素的弦長(zhǎng)和安裝角進(jìn)行了求解，考慮到制造性和經(jīng)濟(jì)性兩方面的因素，又對(duì)弦長(zhǎng)C和安裝角θ進(jìn)行了修正;然后利用Excel辦公軟件采用齊次坐標(biāo)轉(zhuǎn)換的方法對(duì)翼型二維坐標(biāo)進(jìn)行轉(zhuǎn)換，得到空間三維坐標(biāo) (x，y，z);最后應(yīng)用SolidWorks軟件繪出葉片實(shí)體模型。此過程簡(jiǎn)單、通用性強(qiáng)、易于進(jìn)一步分析，為葉片和其他復(fù)雜實(shí)體建模提供了新思路。

【1】王曉蓉，王偉勝，戴慧珠.我國(guó)風(fēng)力發(fā)電現(xiàn)狀和展望［J］.中國(guó)電力，2004，37(1):81－84.

【2】林閩，張崇巍，張艷紅，等.小型風(fēng)力發(fā)電機(jī)葉輪設(shè)計(jì)［J］.風(fēng)機(jī)技術(shù)，2007(1):28－30.

【3】包飛.風(fēng)力機(jī)葉片幾何設(shè)計(jì)與空氣動(dòng)力學(xué)仿真［D］.大連:大連理工大學(xué)，2008:26.

【4】王凡.風(fēng)力發(fā)電機(jī)的葉片設(shè)計(jì)方法研究［D］.南京:南京理工大學(xué)，2007.

【5】王沫然.MATLAB6.0與科學(xué)計(jì)算［M］.北京:電子工業(yè)出版社，2002.

【6】王學(xué)永.風(fēng)力發(fā)電機(jī)葉片設(shè)計(jì)及三維建模［D］.北京:華北電力大學(xué)，2008:34－36.

【7】賀德馨.風(fēng)工程與工業(yè)空氣動(dòng)力學(xué)［M］.北京:國(guó)防工業(yè)出版社，2006.

【8】陳家權(quán)，楊新彥.風(fēng)力機(jī)葉片立體圖設(shè)計(jì)［J］.機(jī)電工程，2006，23(4):37－40.