旋轉(zhuǎn)風(fēng)力機(jī)葉片氣動(dòng)彈性偏移變形的數(shù)值模擬研究

王旭東 陳 進(jìn) 王立存 陸群峰

1．重慶工商大學(xué)廢油資源化技術(shù)與裝備教育部工程研究中心，重慶，400067 2．重慶大學(xué)機(jī)械傳動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，重慶，400030

0 引言

風(fēng)力機(jī)系統(tǒng)作為一個(gè)多自由度系統(tǒng)，在對(duì)其進(jìn)行動(dòng)力學(xué)特性研究時(shí)就必須研究系統(tǒng)的穩(wěn)定性。在進(jìn)行風(fēng)力機(jī)系統(tǒng)設(shè)計(jì)時(shí)，為了保證這個(gè)多自由度系統(tǒng)的穩(wěn)定性，首先需要對(duì)每個(gè)部件的載荷和穩(wěn)定性進(jìn)行研究［1－2］。

風(fēng)力機(jī)葉片作為風(fēng)力機(jī)的關(guān)鍵部件，其良好的設(shè)計(jì)、可靠的質(zhì)量和優(yōu)越的性能是保證機(jī)組正常穩(wěn)定運(yùn)行的決定因素。隨著風(fēng)力機(jī)向大功率發(fā)展，葉片的載荷及振動(dòng)偏移特性的研究就更突出其重要性，而當(dāng)風(fēng)力機(jī)塔架的剛度很大時(shí)，風(fēng)力機(jī)系統(tǒng)的動(dòng)力學(xué)特性主要取決于風(fēng)力機(jī)葉片動(dòng)力學(xué)特性，可以說(shuō)，葉片作為風(fēng)力機(jī)的主要部件，其振動(dòng)偏移和穩(wěn)定性研究直接決定了風(fēng)力機(jī)系統(tǒng)設(shè)計(jì)的可靠性［3－5］。因此，對(duì)風(fēng)力機(jī)葉片的振動(dòng)偏移和穩(wěn)定性進(jìn)行研究就顯得尤為必要。本文從結(jié)構(gòu)動(dòng)力學(xué)理論出發(fā)，基于風(fēng)力機(jī)旋轉(zhuǎn)風(fēng)輪的氣動(dòng)彈性模型，提出了旋轉(zhuǎn)葉片偏移變形特性的數(shù)值模擬計(jì)算理論和方法。

1 旋轉(zhuǎn)風(fēng)輪的載荷特性模型

作用在風(fēng)輪上的空氣動(dòng)力學(xué)載荷包括各個(gè)葉片拍打方向所承受的推力和擺振方向所承受的剪力，推力和剪力的計(jì)算采用動(dòng)量葉素理論，利用Prandtl葉尖損失修正因子和Shen葉尖損失修正模型，對(duì)計(jì)算過(guò)程中葉片每個(gè)葉素中的法向力和切向力系數(shù) （Cn，Ct）進(jìn)行修正［6］。那么，葉片上推力和剪力的計(jì)算公式就可以表示為

它所產(chǎn)生的力矩的計(jì)算公式為

式中，B為葉片個(gè)數(shù)；λ為葉尖速比；c為葉片弦長(zhǎng)。

當(dāng)風(fēng)輪旋轉(zhuǎn)時(shí)，將會(huì)對(duì)葉片產(chǎn)生離心力，并且方向沿葉片向外。同時(shí)，在風(fēng)輪運(yùn)轉(zhuǎn)的過(guò)程中，作用在葉片上的空氣動(dòng)力學(xué)載荷所產(chǎn)生的推力使得葉片偏離旋轉(zhuǎn)平面，進(jìn)而使得作用在葉片上的離心力方向發(fā)生了改變，如圖1所示。此時(shí)葉片所承受的離心力可以表示為

式中，R為風(fēng)輪半徑；Ω為風(fēng)輪轉(zhuǎn)速；mi為風(fēng)輪葉素質(zhì)量。

圖1 葉片的離心力示意圖

從圖1可以看出，葉片離心力在拍打方向的分量為

式中，δ為葉片偏離旋轉(zhuǎn)平面的角度；θ為風(fēng)輪的錐角。

可以看出，由于葉片離心力在拍打方向的分量Fcn與葉片所承受的推力Tx方向相反，因此，該離心力的存在能夠減小拍打彎矩帶來(lái)的葉片偏離距離。

2 風(fēng)輪葉片偏移變形的動(dòng)力學(xué)模型

根據(jù)虛功原理，基于結(jié)構(gòu)動(dòng)力學(xué)理論建立葉片的動(dòng)力學(xué)方程：

式中，M為質(zhì)量矩陣；C為阻尼矩陣；K為剛度矩陣；Fg為外載荷，即空氣動(dòng)力學(xué)載荷和離心力；x為位移向量。考慮葉片前三階振動(dòng)模態(tài)的影響，取11個(gè)自由度，即

式中，w為風(fēng)輪彈性位移；φ為葉片的角度位置；xji為葉片i在第j階模態(tài)時(shí)所產(chǎn)生的位移；i，j＝1，2；上標(biāo)f、e分別表示揮舞方向和擺動(dòng)方向。

通過(guò)對(duì)式（7）的迭代求解得到式（8）中葉片在各自由度方向的位移。對(duì)于葉片在外載荷作用下產(chǎn)生的位移、速度和加速度，采用葉片前三階振動(dòng)模態(tài)的線性關(guān)系表達(dá)為［7］

式中，d為葉片產(chǎn)生的位移；u1f為一階揮舞方向模態(tài)單位特征值；u1e為一階擺振方向模態(tài)單位特征值；u2f為二階揮舞方向模態(tài)單位特征值。

3 旋轉(zhuǎn)風(fēng)輪葉片載荷的數(shù)值模擬

根據(jù)上述建立的風(fēng)輪載荷計(jì)算模型，在MATLAB軟件中采用數(shù)值迭代法對(duì)某2MW風(fēng)力機(jī)風(fēng)輪的載荷進(jìn)行了計(jì)算分析。該風(fēng)輪由3個(gè)葉片組成，每個(gè)葉片長(zhǎng)30．56m，設(shè)計(jì)額定風(fēng)速為14m／s，葉尖速比λ＝5，風(fēng)輪的設(shè)計(jì)旋轉(zhuǎn)速度為2．3rad／s。以葉素理論為基礎(chǔ)，將葉片沿展向劃分成20段，通過(guò)對(duì)風(fēng)輪誘導(dǎo)因子的迭代求解，確定風(fēng)輪各個(gè)葉片的空氣動(dòng)力學(xué)載荷［8］。圖2所示為額定風(fēng)速14m／s，風(fēng)輪旋轉(zhuǎn)12s時(shí)，即旋轉(zhuǎn)了4．43r后，所求得的葉片軸向和周向誘導(dǎo)因子在葉片展向的分布。圖3所示為該時(shí)刻葉片的空氣動(dòng)力學(xué)載荷法向力和切向力分布，其中切向力主要表現(xiàn)為風(fēng)輪的輸出功率，法向力主要表現(xiàn)為風(fēng)輪所承受的推力。風(fēng)輪在該推力的作用下，在垂直旋轉(zhuǎn)平面的方向會(huì)發(fā)生拍打位移，該位移的產(chǎn)生使得葉片的離心力方向發(fā)生變化，進(jìn)而使得葉片的離心力在拍打和擺振方向產(chǎn)生了分量，由于拍打方向的離心力和推力方向相反，因此使得葉片的推力減小，圖4所示為相同時(shí)刻條件下計(jì)算得到的葉片離心力在拍打和擺振方向的分布。

圖2 葉片誘導(dǎo)因子沿展向的分布

圖3 葉片空氣動(dòng)力學(xué)載荷的分布

圖4 葉片離心力的分布

4 風(fēng)輪偏移變形特性的算例分析

為了很好地研究風(fēng)力機(jī)葉片的偏移變形特性和驗(yàn)證該模型的可靠性，對(duì)風(fēng)輪在不同風(fēng)速和不同位置時(shí)刻的葉片的偏移變形特性進(jìn)行了計(jì)算分析。同樣針對(duì)上述2MW風(fēng)力機(jī)，在已得的載荷基礎(chǔ)上，在MATLAB軟件中編寫(xiě)程序?qū)θ~片在該載荷下的偏移變形位移進(jìn)行了數(shù)值模擬。葉片在載荷作用下的偏移變形與其所處的運(yùn)轉(zhuǎn)位置有關(guān)，本文以圖5所示各葉片所處的實(shí)線位置為風(fēng)輪的起始旋轉(zhuǎn)位置，風(fēng)輪以圖中箭頭方向即沿逆時(shí)針?lè)较蛐D(zhuǎn)。圖6所示為風(fēng)輪在額定風(fēng)速14m／s下，運(yùn)轉(zhuǎn)8s時(shí)，風(fēng)輪三個(gè)葉片沿展向的不同位移偏移變形。因運(yùn)轉(zhuǎn)時(shí)間為8s，即剛好旋轉(zhuǎn)了3r，那么該時(shí)刻各葉片所處的位置和起始位置相同。圖7所示為相同風(fēng)速下，風(fēng)輪在旋轉(zhuǎn)12s以后，即旋轉(zhuǎn)了4．43r，各個(gè)葉片處于圖5中虛線位置時(shí)，葉片的展向偏移變形。和圖6對(duì)比，盡管來(lái)流風(fēng)速和風(fēng)輪旋轉(zhuǎn)速度相同，但是由于風(fēng)輪旋轉(zhuǎn)時(shí)間不同，各葉片計(jì)算時(shí)刻所處的位置不一樣，因此，三個(gè)葉片的偏移變形和圖6所示的偏移變形有著不同。圖8所示為該風(fēng)輪運(yùn)轉(zhuǎn)12s時(shí)，風(fēng)速為10m／s下，風(fēng)輪各個(gè)葉片的偏移變形分布?？梢钥闯觯瑘D8中每個(gè)葉片的偏移變形規(guī)律和圖7一致，這是由于兩種工況的運(yùn)轉(zhuǎn)時(shí)間相同，即計(jì)算時(shí)刻每個(gè)葉片所處的位置相同，但是由于圖8風(fēng)速為10m／s，比圖7中風(fēng)速14m／s要小，即圖8葉片所承受的載荷要小些，因此葉片展向所產(chǎn)生的偏移變形也相對(duì)較小。

圖5 風(fēng)輪運(yùn)轉(zhuǎn)時(shí)葉片所處的位置

圖6 額定風(fēng)速14m／s下風(fēng)輪旋轉(zhuǎn)3轉(zhuǎn)時(shí)各個(gè)葉片的偏移變形分布

圖7 額定風(fēng)速14m／s下風(fēng)輪旋轉(zhuǎn)4．43轉(zhuǎn)時(shí)各個(gè)葉片的偏移變形分布

圖8 額定風(fēng)速10m／s下風(fēng)輪旋轉(zhuǎn)4．43轉(zhuǎn)時(shí)各個(gè)葉片的偏移變形分布

5 結(jié)論

（1）應(yīng)用動(dòng)量葉素理論對(duì)風(fēng)力機(jī)風(fēng)輪旋轉(zhuǎn)條件下的空氣動(dòng)力學(xué)載荷和離心力進(jìn)行了理論模型的建立和數(shù)值模擬。在充分考慮該載荷特性條件下，基于風(fēng)輪的結(jié)構(gòu)動(dòng)力學(xué)方程，提出了葉片的偏移變形特性計(jì)算數(shù)學(xué)模型。

（2）針對(duì)建立的葉片偏移變形計(jì)算數(shù)學(xué)模型，采用迭代法對(duì)某2MW風(fēng)力機(jī)在不同工況下的葉片展向不同位置的偏移變形進(jìn)行了數(shù)值模擬，得到了該風(fēng)輪在不同風(fēng)速和不同運(yùn)轉(zhuǎn)時(shí)刻的各個(gè)葉片的偏移變形分布。研究結(jié)果對(duì)風(fēng)力機(jī)葉片的振動(dòng)特性和疲勞壽命預(yù)測(cè)具有理論指導(dǎo)意義。

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

［2］ Maheri A，Noroozi S，Vinney J．Combined Analytical／FEA－based Coupled Aero Structure Simulation of a Wind Turbine with Bend-twist Adaptive Blades［J］．Renewable Energy，2007，32：916－930．

［3］ Larsen J W，Nielsen S R K．Nonlinear Parametric Instability of Wind Turbine Wings［J］．Journal of Sound and Vibration，2007，299：64－82．

［4］ Sarkar S，Bijl H．Nonlinear Aeroelastic Behavior of an Oscillating Airfoil During Stall－induced Vibration［J］．Journal of Fluids and Structures，2008，24：757－777．

［5］ Turhan O，Bulut G．On Nonlinear Vibrations of a Rotating Beam［J］．Journal of Sound and Vibration，2009，322：314－335．

［6］ Shen W Z，Mikkelsen R，Sorensen J N．Tip Loss Correction for Wind Turbines Computations［J］．Wind Energy，2005，8：457－475．

［7］ Hansen M O L．Aerodynamics of Wind Turbines［M］．London：James ＆ James （Science Publishers）Ltd，2000．

［8］ 陳進(jìn)，王旭東，沈文忠，等．風(fēng)力機(jī)葉片的形狀優(yōu)化設(shè)計(jì)［J］．機(jī)械工程學(xué)報(bào)，2010，46（3）：131－134．