黃 娜，關(guān)正偉，黨曉圓，李 潔

(重慶郵電大學(xué)移通學(xué)院，重慶 401520)

1 引言

在風(fēng)力發(fā)電機(jī)的設(shè)計中，首先需要考慮到風(fēng)翼的設(shè)計是否滿足應(yīng)用地區(qū)的條件和發(fā)電機(jī)的功率需求。這就需要在設(shè)計時擬定風(fēng)翼的形狀，葉數(shù)等因素，其中風(fēng)翼的翼型直接影響風(fēng)力發(fā)電機(jī)的風(fēng)力特性參數(shù)[1]。

在設(shè)計中需要對發(fā)電機(jī)翼型進(jìn)行模擬以保證風(fēng)翼可以滿足設(shè)計需求。而對翼型特性的分析中，早期人們采用模擬的方法，文獻(xiàn)[2]提出Gurney襟翼對風(fēng)力機(jī)翼型氣動噪聲影響的數(shù)值模擬，當(dāng)迎角為4°～20°時，原始翼型及帶不同高度襟翼的翼型氣動特性和流場分布通過Fluent軟件LES模型進(jìn)行計算，然后根據(jù)FW-H聲學(xué)模擬法，用Acoustics模塊精確求解遠(yuǎn)場氣動噪聲。文獻(xiàn)[3]提出風(fēng)力機(jī)翼型定常粘性分離流動分析與數(shù)值模擬，針對風(fēng)力機(jī)翼型表面流體分離的實際情況，利用普朗特分離準(zhǔn)則，給出了風(fēng)力機(jī)翼型表面流體分離的充要條件。針對N-S方程所描述的實際粘性流，對其分離點附近的實際粘性流的特性進(jìn)行了分析，得到了分離流線、零u線和零渦線的相對位置，并利用 Fluent流體分析軟件進(jìn)行了數(shù)值計算，得到了S809翼型的氣動特性及其周圍流場分布。文獻(xiàn)[4]提出小攻角時風(fēng)力機(jī)翼型邊界層特性的數(shù)值模擬，發(fā)現(xiàn)沿翼型弦方向，從前緣到后緣，邊界層名義厚度、位移厚度、動量損失和能量損失厚度都有增加的趨勢，沿翼型吸力平面外法線方向，四種厚度的最大值分別為翼型弦長的1.25%、0.36%、0.17%和0.29%；在邊界層內(nèi)，邊界層內(nèi)流體的黏度影響顯著，沿翼型吸力平面外法線方向，流體切向速度的最大值隨邊界層的名義厚度增加而逐漸減少，當(dāng)法向高度大于邊界層的名義厚度時，邊界層速度剪切特性基本保持不變。

雖然上述方法取得一定進(jìn)展，但反映的特性有限。后續(xù)計算機(jī)的應(yīng)用，人們開始利用算式進(jìn)行數(shù)值仿真，為了進(jìn)一步提高仿真結(jié)果的準(zhǔn)確性，研究者利用有限元分析軟件對模擬風(fēng)翼的運動過程和風(fēng)阻狀況來進(jìn)行數(shù)值仿真。但在實際研究中發(fā)現(xiàn)，由于有限元分析軟件無法計算因特定情況下產(chǎn)生的翼身抖振現(xiàn)象，導(dǎo)致數(shù)值模擬分析結(jié)果準(zhǔn)確性不足，難以滿足實際需求。因此對于風(fēng)力發(fā)電機(jī)的翼型特性數(shù)值仿真仍需要進(jìn)一步的改善和研究。

2 直驅(qū)風(fēng)力發(fā)電機(jī)翼型特性數(shù)值仿真方法設(shè)計

2.1 直驅(qū)風(fēng)力發(fā)電機(jī)組整機(jī)結(jié)構(gòu)

直驅(qū)風(fēng)力發(fā)電機(jī)由傳動器、風(fēng)輪、制動器、液壓機(jī)、塔架以及發(fā)電機(jī)等基本設(shè)備組成，其中，傳動器主要由主軸、齒輪箱和聯(lián)軸節(jié)組成；風(fēng)輪主要由葉片和輪轂組成；制動器主要由氣動制動和機(jī)械制動組成；液壓器主要由電機(jī)、油箱、液壓閥等組成；塔架采用筒形結(jié)構(gòu)，底座采用鋼筋混凝土結(jié)構(gòu)；發(fā)電機(jī)主要由異步發(fā)電機(jī)和低速永磁發(fā)電機(jī)組成[5-6]。直驅(qū)風(fēng)力發(fā)電機(jī)組整機(jī)結(jié)構(gòu)如圖1所示。

圖1 直驅(qū)風(fēng)力發(fā)電機(jī)整機(jī)結(jié)構(gòu)圖

風(fēng)力機(jī)負(fù)荷的種類有：靜態(tài)負(fù)荷、穩(wěn)態(tài)負(fù)荷、周期負(fù)荷、極限負(fù)荷和疲勞負(fù)荷。這些負(fù)荷中，靜態(tài)負(fù)荷為風(fēng)電機(jī)組不轉(zhuǎn)動時所產(chǎn)生的負(fù)荷，穩(wěn)態(tài)負(fù)荷為風(fēng)電機(jī)組轉(zhuǎn)動時所產(chǎn)生的負(fù)荷，周期負(fù)荷為機(jī)組運行時所產(chǎn)生的負(fù)荷，周期負(fù)荷為機(jī)組結(jié)構(gòu)在某一周期變化期間所產(chǎn)生的負(fù)荷，負(fù)荷變化周期與風(fēng)輪轉(zhuǎn)速成正比[7]；極限負(fù)荷需要綜合考慮實際情況，區(qū)分機(jī)組正常狀態(tài)和故障狀態(tài)；兆瓦級風(fēng)電機(jī)組整機(jī)每轉(zhuǎn)一周，產(chǎn)生低速軸和葉片重力完全相反的力，從而產(chǎn)生疲勞負(fù)荷。

2.2 風(fēng)力空氣動力模型

在進(jìn)行風(fēng)力發(fā)電機(jī)翼型特性的數(shù)值仿真中，首先要確定風(fēng)力空氣動力模型，并作為數(shù)值仿真的運算基礎(chǔ)[8]。在本文研究中，根據(jù)動力理論的理想翼型的氣流模型，如圖2所示。

圖2 氣流管樁模型圖

在動力理論中的翼型理想氣流模型中，空氣被設(shè)置為不可壓縮的狀態(tài)，同時在氣流流經(jīng)風(fēng)翼時，氣流質(zhì)量相等，并在單位時間內(nèi)通過截面，而此時的氣流質(zhì)量為ρAU，而氣流通過風(fēng)輪時則應(yīng)滿足

ρA∞U∞=ρAdUd=ρAwUw

(1)

在式(1)中，ρ代表當(dāng)前氣流中的空氣密度，A代表氣流的截面面積，U代表當(dāng)前流經(jīng)風(fēng)翼的氣流流速，∞代表風(fēng)翼的無窮遠(yuǎn)處，d代表風(fēng)力翼型的風(fēng)輪圓盤處，w代表風(fēng)翼的下游尾端。由于風(fēng)翼的阻擋，導(dǎo)致氣流在經(jīng)過風(fēng)翼時，氣流的流速出現(xiàn)變化，從而在運算中導(dǎo)入軸向誘導(dǎo)因子a，以及風(fēng)翼對氣流的軸線速度誘變aU∞，而氣流流經(jīng)風(fēng)翼的速度則為

Ud=(1-a)U∞

(2)

而根據(jù)氣體動量理論，氣流經(jīng)過風(fēng)翼后的動力改變則為

ΔP=ρAdUd(U∞-Uw)

(3)

根據(jù)空氣動量改變公式，可以推出風(fēng)翼經(jīng)過的壓力差值，如下所示

(p+-p-)Ad=ρAdUd(1-a)(U∞-Uw)

(4)

而式(4)的壓力差，也可以反映出當(dāng)前壓力和風(fēng)能之間的關(guān)系，從而確定當(dāng)前風(fēng)翼翼型的利用系數(shù)Cp定義為

(5)

而風(fēng)翼利用系數(shù)的最大值，則代表風(fēng)翼翼型的理想狀態(tài)。

2.3 數(shù)值仿真翼型參數(shù)確定

風(fēng)翼的葉片從轉(zhuǎn)動中心到葉尖半徑中的剖面翼型弦長為Li，通常條件下風(fēng)翼葉片可以將空氣動力的平均分配至整個葉片中，這往往是受到葉片的扭曲形態(tài)以及翼片的翼型弦長等因素的影響[9]。因此需要將風(fēng)翼葉片中參數(shù)更全面的輸入。而葉片在轉(zhuǎn)動狀態(tài)下，葉尖在不同半徑下的葉片弦長Li的計算為

(6)

在式(6)中，r代表風(fēng)翼葉片的轉(zhuǎn)動圓心和該葉片葉尖之間的距離，單位為m。Ce代表風(fēng)翼葉片的形狀參數(shù)，該數(shù)值可經(jīng)由風(fēng)葉葉片的尖速值來獲得。CL代表風(fēng)翼葉片的升力系數(shù)，該數(shù)值可由翼型的升阻比曲線，并在其中選取最佳比值來獲得[10]。進(jìn)而根據(jù)葉片的實際安裝效果，給出葉片安裝時的安裝角度，即：θi=φi-αm，其中θi代表風(fēng)翼轉(zhuǎn)動中心在葉尖的不同位置中，半徑ri對應(yīng)的安裝角度，φi代表在風(fēng)翼的ri處對應(yīng)的葉片迎風(fēng)角，αm代表風(fēng)翼葉片的平均迎角。其中αm值需要進(jìn)行另外的計算，計算方法為

(7)

2.4 翼身抖振非定常數(shù)計算

由于預(yù)設(shè)了風(fēng)力空氣動力模型，因此可以根據(jù)該模型來結(jié)合當(dāng)前翼身特定狀態(tài)下的氣動特性，來計算風(fēng)翼的翼身抖振非定常數(shù)，首先確定風(fēng)翼翼身的邊界時間步條件，如下所示

△t=(0.03×c)/Um

(8)

在(8)中，c代表風(fēng)翼翼身的最大弦長，Um代表翼身來流速度。同時采用UDF方法來改變風(fēng)翼自由流角方向，并定義其中流迎角α隨時間t的規(guī)律產(chǎn)生的變化。根據(jù)K-ε湍流能方程，得出

(9)

在式(9)中，k代表其中的湍流動能，ε代表風(fēng)翼的湍流消耗率，Gk代表在不同的風(fēng)速影響下湍流的生成狀態(tài)，Gb代表受到浮力的影響導(dǎo)致的湍流動能生成狀態(tài)，YM代表湍流在空氣氣流出現(xiàn)脈動現(xiàn)象時導(dǎo)致的膨脹的產(chǎn)生率，σk代表普朗特常數(shù)，μt代表湍流狀態(tài)下的粘性因素。在進(jìn)行數(shù)值仿真時，需要將相應(yīng)的計算區(qū)域劃分為旋轉(zhuǎn)域以及靜止域，在計算中，控制體內(nèi)發(fā)電風(fēng)翼的外部區(qū)域，而旋轉(zhuǎn)域則是考慮風(fēng)翼在旋轉(zhuǎn)坐標(biāo)系內(nèi)的求解計算，并在靜止域中，經(jīng)過非旋轉(zhuǎn)坐標(biāo)系對控制進(jìn)行求解計算，并在不同區(qū)域中對風(fēng)翼的抖振非定常數(shù)進(jìn)行求解[12]。求解完成后，再進(jìn)行數(shù)據(jù)過交界面進(jìn)行傳遞，而在靜止域中，風(fēng)力發(fā)電機(jī)的空氣流暢計算控制方程則需要采用翼身抖振非定常數(shù)進(jìn)行計算，如下所示

(10)

式中，φ代表氣流輸送變量，pa代表風(fēng)翼空氣密度，ua則在靜止域內(nèi)的空氣流場速度，Γφ代表空氣流場內(nèi)的擴(kuò)散系數(shù)，qφ代表靜止域流場內(nèi)的源項。而在旋轉(zhuǎn)域中則可以寫為

qφr=qφ-ρ[2ω×uar+ω×(ω×r)]

(11)

式中，uar代表空氣流場內(nèi)旋轉(zhuǎn)坐標(biāo)系的速度，qφr帶旋轉(zhuǎn)坐標(biāo)中源項，而當(dāng)式(9)中的φ=1、qφr=0時，風(fēng)翼則處于連續(xù)狀態(tài)，則流向則作為源項動量來計算。

2.5 數(shù)值仿真機(jī)翼生成

在數(shù)值仿真中，往往風(fēng)翼葉片的生成過程通過點線面的生成過程，而在使用該方法進(jìn)行數(shù)值仿真過程中，由于需要逐點進(jìn)行生成，因此往往工作量較大，且費時費力[7]。本文中選擇采用Pre/E方法，來提供數(shù)據(jù)文件生成曲線的造型方式，并將對應(yīng)的翼型數(shù)據(jù)根據(jù)相應(yīng)的格式編寫為IBL格式，并導(dǎo)入對應(yīng)的Pro/E的翼型邊界曲線，由于在Fluent的分析過程中，使用的單位為m，因此在Pro/E中進(jìn)行二次導(dǎo)入的數(shù)據(jù)也應(yīng)為m為單位，同時在單一組別中的發(fā)電機(jī)翼型需要選擇相同的旋轉(zhuǎn)坐標(biāo)系，并生成對應(yīng)的葉片截面邊界，如圖3所示。

圖3 葉片截面邊界載入方式

葉片外表面，則可以運用相同的載入方式來實現(xiàn)，并對邊界進(jìn)行混合，依次選取各個截面上存在的邊界曲線，并生成在葉片的下表面上，僅葉片上下表面分別進(jìn)行生成，由于存在翼型尾側(cè)尖角，葉片容易出現(xiàn)局部扭曲，在生成時可以將兩個端部進(jìn)行填充，并獲取端部截面。經(jīng)過該方法獲得的翼型，無法進(jìn)行布爾運算，因此后續(xù)需要在數(shù)值仿真中將葉片實現(xiàn)實體化，并生成對應(yīng)的風(fēng)翼細(xì)節(jié)。同時采用Gambit軟件來輸入相關(guān)參數(shù)，并定義翼型中存在的復(fù)雜扭曲面，將保存的STP文件導(dǎo)入Gambit中，并利用其中幾何修正方法，連接其中的點、線、面，并保證原始的集合精度，通過軟件功能將其中的小縫隙進(jìn)行縫合，實現(xiàn)對風(fēng)力發(fā)電機(jī)翼型的仿真生成，并導(dǎo)入至上述模型中，實現(xiàn)數(shù)值仿真。

3 仿真分析

為了驗證設(shè)計的數(shù)值仿真方法的可行性，本文以某地區(qū)的直驅(qū)風(fēng)力發(fā)電機(jī)作為實驗對象，并使用文獻(xiàn)[2]、文獻(xiàn)[3]、文獻(xiàn)[4]中的數(shù)值仿真方法以及本文設(shè)計的數(shù)值仿真方法對該發(fā)電機(jī)的翼型特性進(jìn)行數(shù)值仿真。

3.1 仿真數(shù)值來源

仿真中的數(shù)值由該直驅(qū)風(fēng)力發(fā)電機(jī)中的技術(shù)參數(shù)以及該風(fēng)力發(fā)電機(jī)運行時接收到的數(shù)據(jù)作為數(shù)值仿真的技術(shù)來源。

3.2 風(fēng)力發(fā)電機(jī)技術(shù)參數(shù)

該風(fēng)力發(fā)電機(jī)的額定功率為400kW，風(fēng)力發(fā)電機(jī)所在環(huán)境風(fēng)速在3～15m/s左右，且發(fā)電機(jī)的風(fēng)能利用系數(shù)為0.41Cp。風(fēng)能轉(zhuǎn)換效率為0.76。地區(qū)空氣密度為1.316kg/m3。本文實驗中的風(fēng)力發(fā)電機(jī)，葉片數(shù)為3片，葉片之間夾角為120o，葉片質(zhì)量為均勻?qū)ΨQ分布。風(fēng)輪實度為12%。該發(fā)電機(jī)的翼型為NACA63-2XX系列，雷諾數(shù)為106。

3.3 發(fā)電機(jī)翼型特性仿真

本文實驗就當(dāng)?shù)爻Ｒ姷娜N不同的風(fēng)場環(huán)境時，發(fā)電機(jī)翼型的特性進(jìn)行數(shù)值仿真。第一種風(fēng)場環(huán)境為北-西風(fēng)向，風(fēng)速在11～12m/s之間，數(shù)值仿真結(jié)果如表1所示。

表1 第一種風(fēng)場環(huán)境數(shù)值仿真結(jié)果

第二種風(fēng)場環(huán)境中，風(fēng)向為南-西風(fēng)向，風(fēng)速為6～7m/s左右，數(shù)值仿真結(jié)果如表2所示。

表2 第二種風(fēng)場環(huán)境數(shù)值仿真結(jié)果

第三種風(fēng)場環(huán)境中，風(fēng)向為北-東風(fēng)向，風(fēng)速為12～13m/s左右，數(shù)值仿真結(jié)果如表3所示：

表3 第三種風(fēng)場環(huán)境數(shù)值仿真結(jié)果

在上述實驗結(jié)果中數(shù)值仿真方法1代表本文設(shè)計的數(shù)值仿真方法，數(shù)值仿真方法2代表文獻(xiàn)[2]中的數(shù)值仿真方法，數(shù)值仿真方法3代表文獻(xiàn)[3]中的數(shù)值仿真方法，數(shù)值仿真方法4代表文獻(xiàn)[4]中的數(shù)值仿真方法。在上述實驗結(jié)果中可以發(fā)現(xiàn)，本文設(shè)計的數(shù)值仿真方法獲得的翼型特性數(shù)值更接近實際的特性系數(shù)，證明本文設(shè)計的數(shù)值仿真方法具有較高的準(zhǔn)確性。

4 結(jié)束語

利用通過改善對發(fā)電機(jī)翼的抖振非常數(shù)值計算，提出了數(shù)值仿真方法，同時實驗證明該方法提高了數(shù)值仿真方法的準(zhǔn)確性。但本文研究中，由于采用了風(fēng)力空氣動力模型作為仿真的基礎(chǔ)模型，因此設(shè)計的數(shù)值仿真方法，無法對機(jī)翼自身磨損以及軸承老化等因素進(jìn)行模擬，仍需要進(jìn)一步的改進(jìn)。