劉一鳴，李 春，王淵博，李 根，閆陽天

（上海理工大學(xué) 能源與動力工程學(xué)院/上海市動力工程多相流與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093）

化石燃料的過度使用及能源結(jié)構(gòu)的不合理導(dǎo)致溫室效應(yīng)和環(huán)境污染等問題，且在不久的將來還可能導(dǎo)致資源枯竭之境況。為此，國家“十四五”規(guī)劃綱要對發(fā)展風(fēng)電等清潔能源提出了明確要求。不僅我國，世界各國均遭遇同樣問題并都積極開發(fā)風(fēng)能資源。此外，現(xiàn)代風(fēng)力機(jī)通常均安裝葉片變槳系統(tǒng)，然而變槳系統(tǒng)為風(fēng)力機(jī)所有部件中故障率最高，且一旦發(fā)生故障將導(dǎo)致風(fēng)力機(jī)無法正常運(yùn)行時間最長的關(guān)鍵子系統(tǒng)之一，但鮮見與此相關(guān)氣動及結(jié)構(gòu)方面的研究。因此，在同時兼顧流體域和固體域基礎(chǔ)上，對風(fēng)力機(jī)特定狀態(tài)下變槳故障葉片進(jìn)行數(shù)值計算以探究其氣動特性及結(jié)構(gòu)響應(yīng)顯得尤為重要。

對于風(fēng)力機(jī)變槳故障，國內(nèi)外學(xué)者進(jìn)行了一些研究。Ronsten通過實(shí)測相關(guān)數(shù)據(jù)提出風(fēng)力機(jī)停機(jī)后葉片狀態(tài)將影響其力矩大小及壓力分布。據(jù)此不難推斷，停機(jī)后變槳故障與變槳成功葉片由于狀態(tài)差別極大，受載亦將存在較大差異。Tian等基于數(shù)值模擬方法并采用獨(dú)立變槳技術(shù)控制大型三葉片水平軸風(fēng)機(jī)運(yùn)行，結(jié)果證明該策略可有效減小葉片載荷變化范圍及風(fēng)輪轉(zhuǎn)矩波動幅度。因此，該獨(dú)立變槳技術(shù)能保證葉片時刻處于最佳攻角，旋轉(zhuǎn)過程中葉片狀態(tài)始終較為穩(wěn)定，而風(fēng)力機(jī)發(fā)生變槳故障時則與此情況相反。Chen等通過研究臺風(fēng)“天兔”襲擊汕尾紅海灣風(fēng)電場事件，指出變槳系統(tǒng)故障與風(fēng)力機(jī)大面積嚴(yán)重?fù)p毀情況存在很大聯(lián)系。閆學(xué)勤等采取將葉片方位角權(quán)系數(shù)分配與葉根氣動載荷反饋相結(jié)合的變槳新技術(shù)進(jìn)行葉根揮舞力矩及輪轂傾覆力矩控制，結(jié)果顯示新型獨(dú)立變槳技術(shù)可保證葉片具有低變槳故障率和小變槳誤差。任年鑫等研究了臺風(fēng)達(dá)到兆瓦級水平軸風(fēng)力機(jī)切出風(fēng)速時，風(fēng)力機(jī)順槳停機(jī)后變槳故障與變槳成功葉片的受載情況，結(jié)果證明，風(fēng)力機(jī)葉輪方位角及臺風(fēng)切入方向會對葉片受力產(chǎn)生影響，且變槳故障葉片載荷總大于變槳成功葉片。

由上述研究可知，風(fēng)力機(jī)變槳故障對葉片受載影響極大，故此時葉片結(jié)構(gòu)亦最危險。然而，在絕大多數(shù)研究中僅考慮了氣動側(cè)載荷卻未涉及結(jié)構(gòu)側(cè)響應(yīng)。此外，切出風(fēng)速下變槳故障葉片受載及入流風(fēng)速極大，但該風(fēng)況下的風(fēng)力機(jī)變槳故障研究鮮見相關(guān)報道。因此，本文針對NREL 5 MW風(fēng)力機(jī)在切出風(fēng)速下變槳故障葉片展開研究，對比分析變槳成功及各典型方位角下變槳故障葉片氣動特性，然后開展變槳故障葉片準(zhǔn)靜態(tài)結(jié)構(gòu)響應(yīng)分析，以期為葉片結(jié)構(gòu)設(shè)計提供更全面的參考數(shù)據(jù)。

1 研究對象及數(shù)值計算模型

1.1 葉片參數(shù)化建模及鋪層

葉片是風(fēng)力機(jī)最為關(guān)鍵的部件。本文利用Unigraphics NX軟件建立NREL 5 MW風(fēng)力機(jī)葉片三維實(shí)體模型，如圖1所示，其中對關(guān)鍵部位進(jìn)行了放大顯示。

圖1 NREL 5 MW風(fēng)力機(jī)葉片模型Fig. 1 Blade model of NREL 5 MW wind turbine

由于本文結(jié)構(gòu)側(cè)計算涉及有限元分析，因此在Unigraphics NX軟件中完成NREL 5 MW風(fēng)力機(jī)三維實(shí)體建模后，還需對其賦予相關(guān)材料屬性。共選取9種單一材料，通過對材料種類、纖維方向以及鋪層順序進(jìn)行組合便可得到適用于風(fēng)力機(jī)葉片各部位的復(fù)合材料。例如，葉片主梁帽需要承擔(dān)彎曲載荷，可以選擇抗拉性強(qiáng)的碳纖維和玻璃纖維并使纖維方向順著葉片展向進(jìn)行鋪層，便可滿足梁帽的結(jié)構(gòu)力學(xué)性能要求，其他部位鋪層方法類似。最終葉片鋪層方案如圖2所示。通過Abaqus軟件Property模塊的Composite Layup功能進(jìn)行NREL 5 MW風(fēng)力機(jī)葉片復(fù)合材料鋪層。這與真實(shí)風(fēng)力機(jī)葉片結(jié)構(gòu)設(shè)計和制造工藝相同，從而可保證計算結(jié)果能夠直接為相關(guān)設(shè)計制造及實(shí)際工程問題提供參考。

圖2 葉片鋪層方案Fig. 2 Blade layout

1.2 氣動側(cè)模型

氣動載荷分析需通過與葉片同步旋轉(zhuǎn)的局部坐標(biāo)系定義各方向力矩。葉片坐標(biāo)系及各方向力矩如圖3所示。

圖3 葉片坐標(biāo)系及各方向力矩Fig. 3 Blade coordinate system and torque in all directions

為降低計算域邊界對數(shù)值模擬的影響并減小風(fēng)洞阻塞效應(yīng)，需盡可能擴(kuò)大計算域尺寸以增加風(fēng)力機(jī)葉片與各邊界間的距離。本文中計算域以及邊界條件如圖4所示，圖中右下角坐標(biāo)系為數(shù)值模擬全局坐標(biāo)系，為風(fēng)輪直徑126 m。入口邊界條件為指數(shù)率風(fēng)剪切速度入口，參高度處是NREL 5 MW風(fēng)力機(jī)25 m·s的切出風(fēng)速，地面粗糙度指數(shù)選為0.16，出口為壓力出口，其值為1個準(zhǔn)大氣壓.地面為無滑移壁面，頂面是對稱平面。兩側(cè)面采用周期性邊界條件，即計算域一側(cè)面網(wǎng)格節(jié)點(diǎn)與另一側(cè)面網(wǎng)格節(jié)點(diǎn)一一對應(yīng)，某側(cè)周期邊界計算域外“鏡像單元”信息由緊鄰另一側(cè)周期邊界計算域內(nèi)的單元提供。CCM+前處理創(chuàng)建高質(zhì)量多面體網(wǎng)格，且壁面邊界層采用低處理（＜1）。近壁面網(wǎng)格高度為5×10m，邊界層網(wǎng)格共24層，總厚度為0.1 m。同時為使風(fēng)力機(jī)數(shù)值模擬準(zhǔn)確度更高，對其尾跡

圖4 計算域以及邊界條件Fig. 4 Calculation domain and boundary conditions

計算域最終網(wǎng)格分布如圖5所示。由圖中可看出，相較于采用四面體網(wǎng)格，全部采用多面體網(wǎng)格類型能以更少的網(wǎng)格數(shù)量達(dá)到相同的計算精度。本文中氣動側(cè)均利用計算流體力學(xué)軟件STAR-區(qū)域作了加密處理。

圖5 網(wǎng)格分布Fig. 5 Grid distribution

數(shù)值計算采用隱式非定常算法（implicit unsteady），時間步長為0.001 s，對流項與時間項的離散均使用二階格式，壓力-速度解耦基于SIMPLE算法，湍流處理采用SST模型，湍流強(qiáng)度取0.05，湍流長度尺度為1.0 m。當(dāng)數(shù)值模擬運(yùn)行至載荷處于相對穩(wěn)定或呈周期性振蕩時，則判定為已收斂。

1.3 結(jié)構(gòu)側(cè)模型

圖6為流固耦合結(jié)構(gòu)側(cè)計算模型。風(fēng)力機(jī)葉片幾何模型由片體（無厚度曲面）構(gòu)建，這是由于有限元分析中若三維結(jié)構(gòu)某一方向尺寸遠(yuǎn)小于另外兩個方向尺寸即可采用片體模型，如此既能保證精度又可節(jié)省相當(dāng)多的計算資源。同樣為節(jié)省計算資源，將風(fēng)力機(jī)葉片劃分為結(jié)構(gòu)化網(wǎng)格。與非結(jié)構(gòu)化網(wǎng)格相比，結(jié)構(gòu)化網(wǎng)格生成速度快且質(zhì)量高，基于結(jié)構(gòu)化網(wǎng)格數(shù)值計算的數(shù)據(jù)存儲結(jié)構(gòu)更簡單且占用計算內(nèi)存較少。有限元分析中片體結(jié)構(gòu)化網(wǎng)格適合使用S4R單元，其為一種通用殼單元類型，既可應(yīng)用于厚殼問題求解，亦能應(yīng)用于薄殼模型計算。此外，將葉根一圈通過運(yùn)動耦合約束在葉根圓形的中心點(diǎn)（即約束點(diǎn)）上，然后在該約束點(diǎn)施加ENCASTRE邊界條件對六個自由度同時進(jìn)行約束。

圖6 結(jié)構(gòu)側(cè)計算模型Fig. 6 Calculation model for blade structure

2 準(zhǔn)靜態(tài)結(jié)構(gòu)分析方法

2.1 雙向弱流固耦合

雙向弱流固耦合流程如圖7所示。該流程共有8個步驟，對各步驟的解釋分別為：①對固體域進(jìn)行網(wǎng)格劃分及物理模型設(shè)置，并參考?xì)鈩佑嬎憬Y(jié)果施加一個近似載荷條件，然后單獨(dú)運(yùn)行一段時間以驗(yàn)證結(jié)構(gòu)側(cè)相關(guān)設(shè)置的準(zhǔn)確性及可行性；②流體域保留氣動計算數(shù)據(jù)結(jié)果，并對多面體網(wǎng)格應(yīng)用Morphing功能以控制其變形；③判斷流固兩側(cè)柔性結(jié)構(gòu)是否完全匹配，若一切正常便可通過Import mode程序?qū)С隽黧w域柔性結(jié)構(gòu)所受壓力和黏性剪力，并通過abaqus關(guān)鍵字加載至固體域柔性結(jié)構(gòu)上；④對結(jié)構(gòu)側(cè)進(jìn)行計算，柔性結(jié)構(gòu)變形等數(shù)據(jù)保存至ODB文件；⑤首先將結(jié)構(gòu)側(cè)的ODB文件加載至Import mode程序，其次由Import mode程序提取其中變形數(shù)據(jù)傳遞給STAR-CCM + ，最后STAR-CCM +對柔性結(jié)構(gòu)形狀進(jìn)行更新，同時利用Morphing功能控制計算域網(wǎng)格變形以適應(yīng)當(dāng)前柔性結(jié)構(gòu)形狀；⑥對流體域進(jìn)行計算并得到柔性結(jié)構(gòu)新的壓力以及黏性剪力；⑦重復(fù)前面第③～⑥步直至流固兩側(cè)相關(guān)特征不再有明顯變化為止（本文中流固兩側(cè)數(shù)據(jù)共交換4次）；⑧雙向弱流固耦合數(shù)值模擬完成。

圖7 雙向弱流固耦合流程Fig. 7 Flow chart of two-way weak fluid-structure coupling

2.2 屈曲分析

細(xì)長結(jié)構(gòu)體由于長度方向尺寸遠(yuǎn)大于另外兩個方向尺寸而表現(xiàn)出柔性特征，當(dāng)受到微小載荷或擾動后便易發(fā)生失穩(wěn)。當(dāng)這些外部微小影響消失后，細(xì)長彈性體能夠恢復(fù)至穩(wěn)態(tài)，但若載荷或擾動較大時，細(xì)長彈性體則很可能發(fā)生屈曲。目前風(fēng)力機(jī)大型化趨勢非常明顯，風(fēng)力機(jī)越來越長導(dǎo)致其柔性增加，NREL 5 MW風(fēng)力機(jī)葉片長達(dá)61.5 m，屬于標(biāo)準(zhǔn)的細(xì)長彈性體，因此有必要對其進(jìn)行屈曲分析。屈曲分析一般分為線性與非線性，基于線性屈曲所研究的結(jié)構(gòu)并不會發(fā)生塑性變形，而是在結(jié)構(gòu)彈性階段產(chǎn)生彈性失穩(wěn)。線性屈曲分析最終所得臨界載荷可能會稍大于結(jié)構(gòu)發(fā)生屈曲時的實(shí)際值，但線性屈曲分析計算效率非常高，且可作為結(jié)構(gòu)分析時的前期預(yù)估值，故本文采用線性屈曲分析。

3 結(jié)果分析

3.1 單葉片模型可行性驗(yàn)證

首先對比全葉輪模型（full rotor model，F(xiàn)RM）與單葉片模型（single blade model，SBM）的計算結(jié)果，以便為減少計算量的策略提供理論依據(jù)。圖8為SBM與FRM揮舞力矩和擺振力矩時域圖。對比圖8（a）、（b）中的力矩可知：SBM與FRM揮舞力矩和擺振力矩時域曲線基本重合，初步證明了基于SBM計算的可行性。

圖8 SBM與FRM揮舞力矩和擺振力矩時域圖Fig. 8 Time domain spectra of flapwise moment and edgewise moment for SBM and FRM

將SBM和FRM揮舞力矩和擺振力矩脈動時域序列進(jìn)行快速傅里葉變換，可得到表1中的SBM與FRM載荷脈動主要刺激頻率。由表中可以看出，SBM與FRM揮舞力矩和擺振力矩主要刺激頻率相同或極其接近，這更進(jìn)一步證明了水平軸風(fēng)力機(jī)SBM與FRM載荷脈動差異非常小。

表1 SBM與FRM載荷脈動主要刺激頻率Tab. 1 Main excitation frequency of load pulsation for SBM and FRM (單位：Hz)

經(jīng)過葉片頂點(diǎn)且與葉片坐標(biāo)系軸垂直平面上SBM與FRM在軸方向的葉尖渦量場如圖9所示。從圖中可觀察到，SBM與FRM葉尖渦大小及位置幾乎完全相同，說明SBM與FRM流場流動狀況極為接近，此亦為SBM與FRM載荷頻域主要刺激頻率以及力矩時域曲線差別甚微的根本原因。綜上可知，使用SBM進(jìn)行數(shù)值模擬具有較高的準(zhǔn)確性及可行性，故后文均基于SBM開展計算。

圖9 SBM與FRM在Z軸方向的葉尖渦量場Fig. 9 Blade tip vortex of SBM and FRM in Z-axis direction

3.2 變槳故障葉片氣動特性

圖10為NREL 5 MW風(fēng)力機(jī)在切出風(fēng)速下變槳故障與變槳成功葉片揮舞力矩和擺振力矩。對比圖10（a）、（b）可看出，揮舞力矩比擺振力矩大1個數(shù)量級，但都表現(xiàn)出變槳故障葉片力矩大于變槳成功葉片的特征，其中變槳故障葉片揮舞力矩平均值為變槳成功葉片的13.8倍。

為進(jìn)一步分析變槳故障葉片與變槳成功葉片的差異，對周圍流場速度云圖進(jìn)行分析，結(jié)果如圖11所示。由圖中可知，變槳故障葉片的尾跡較變槳成功葉片的更為明顯。這是因?yàn)樽儤收先~片從來流風(fēng)中吸收了更多能量，因此變槳故障葉片的受載比變槳成功葉片的大，這一結(jié)果與圖10中的數(shù)據(jù)非常吻合。綜上，后文僅分析變槳故障葉片特征。

圖10 變槳故障與變槳成功葉片揮舞力矩和擺振力矩Fig. 10 Flapwise moment and edgewise moment of a blade under fault and normal status of pitch system

圖11 變槳故障葉片與變槳成功葉片周圍流場速度云圖Fig. 11 Velocity contour around a blade under fault and normal status of pitch system

NREL 5 MW風(fēng)力機(jī)達(dá)到切出風(fēng)速停機(jī)后葉輪方位角是隨機(jī)的。圖12為0°、60°、120°、180°、240°和300°共6個方位角下變槳故障葉片揮舞力矩和擺振力矩時域圖。從圖12（a）中可看出，0°方位角下?lián)]舞力矩最大，其次是60°和300°方位角下的，然后是120°、240°、180°方位角下?lián)]舞力矩最小。揮舞力矩在各方位角下的差異主要是由風(fēng)剪切造成風(fēng)速從地面沿豎直方向越來越大所致，故0°方位角葉片平均來流風(fēng)速最大，180°方位角葉片平均來流風(fēng)速最小。此外，60°與300°方位角以及120°與240°方位角葉片平均來流風(fēng)速兩兩相同，因此它們的力矩大小也一樣，故后續(xù)研究僅分析60°和120°方位角，對300°和240°方位角的結(jié)果不再贅述。

圖12 各方位角下變槳故障葉片揮舞力矩及擺振力矩時域圖Fig. 12 Time domain spectra of flapwise moment and edgewise moment of a blade with pitch fault at different azimuths

3.3 準(zhǔn)靜態(tài)結(jié)構(gòu)響應(yīng)分析

圖13為雙向弱流固耦合葉片應(yīng)力及葉尖位移。由圖中可知，切出風(fēng)速下4個典型方位角的變槳故障葉片均出現(xiàn)了應(yīng)力集中現(xiàn)象，其中0°方位角下變槳故障葉片應(yīng)力集中最強(qiáng)，為9.886 MPa，而180°方位角下應(yīng)力集中最弱，為6.932 MPa，后者較前者降低了29.8%。葉尖位移變化趨勢與應(yīng)力集中變化趨勢相同，從0°、60°、120°到180°方位角下葉尖位移逐漸減小，其中180°方位角下變槳故障葉片葉尖位移較0°方位角下的減少了32.7%。此外，還發(fā)現(xiàn)雙向弱流固耦合結(jié)構(gòu)側(cè)結(jié)果與圖12中氣動力矩情況吻合較好，證明了本文計算結(jié)果的準(zhǔn)確性。

圖13 雙向弱流固耦合葉片應(yīng)力及葉尖位移Fig. 13 Stress and tip displacement of a blade under two-way weak fluid-structure coupling

圖14為0°方位角下變槳故障葉片前六階屈曲因子和屈曲模態(tài)。屈曲因子的物理意義是載荷與其之乘積即為葉片發(fā)生屈曲的臨界載荷。從圖中可看出，最小的屈曲因子是一階屈曲因子1.234 5，說明葉片并未發(fā)生屈曲。此外，仔細(xì)觀察圖14能夠發(fā)現(xiàn)，隨著屈曲階數(shù)的增加，葉片發(fā)生屈曲部位的面積越來越大，且屈曲因子逐漸提升，但實(shí)際工程中最關(guān)注的是一階屈曲因子，因?yàn)橐浑A屈曲發(fā)生之后其他階結(jié)構(gòu)屈曲將無意義。此外，其他方位角結(jié)構(gòu)屈曲情況與0°方位角變化趨勢一樣，且隨著方位角增加葉片屈曲越弱，其中180°方位角下變槳故障葉片的一階屈曲因子較0°方位角下的增大20.2%，結(jié)構(gòu)安全性更高。需要說明的是，本文結(jié)構(gòu)側(cè)屈曲分析載荷來自雙向弱流固耦合最后一次流固兩次數(shù)據(jù)交換并計算結(jié)束后氣動側(cè)的載荷數(shù)據(jù)。

圖14 0°方位角下變槳故障葉片前六階屈曲Fig. 14 The first six order buckling of a blade with pitch fault at 0°azimuth

4 結(jié) 論

本文基于計算流體力學(xué)方法對NREL 5 MW風(fēng)力機(jī)變槳故障/成功葉片氣動側(cè)狀態(tài)進(jìn)行對比分析，并利用雙向弱流固耦合及曲屈分析對典型方位角下變槳故障葉片特征展開研究。主要結(jié)論有：

（1）SBM與FRM揮舞力矩和擺振力矩時域序列、頻域特征及渦量場強(qiáng)度和分布差別均非常小或相同，故采用節(jié)省大量計算資源的SBM模型可行。

（2）切出風(fēng)速下變槳故障葉片較變槳成功葉片受載更大，且前者流場尾跡比后者亦更明顯，說明變槳故障葉片從來流風(fēng)中吸收了更多能量。

（3）180°方位角下變槳故障葉片較0°方位角下變槳故障葉片應(yīng)力及葉尖位移分別減小29.8%和32.7%，一階屈曲因子增加20.2%。