摘要： 風(fēng)力機(jī)葉片預(yù)彎設(shè)計(jì)大多采用靜氣彈分析方法，忽略了葉片氣動(dòng)力、慣性力和彈性力三者相互作用而引發(fā)的氣彈耦合失穩(wěn)問題，而預(yù)彎對(duì)百米量級(jí)超長(zhǎng)柔性葉片的顫振性能影響尤為顯著。為對(duì)比分析不同預(yù)彎尺寸對(duì)葉片顫振臨界狀態(tài)的影響，基于主梁剛度等效原則進(jìn)行葉片氣彈模型設(shè)計(jì)，通過風(fēng)洞試驗(yàn)發(fā)現(xiàn)了15 MW兩種預(yù)彎葉片顫振區(qū)間及其臨界風(fēng)速的差異；進(jìn)一步基于修正的葉素動(dòng)量理論和幾何精確梁理論（Blade Element Momentum Theory?Geometrically Exact Beam Theory，BEM?GEBT）耦合計(jì)算方法對(duì)4種預(yù)彎葉片進(jìn)行分析，對(duì)比了不同預(yù)彎尺寸葉片的顫振臨界風(fēng)速、氣動(dòng)力分布和位移頻譜特性，并揭示了顫振耦合模態(tài)機(jī)理。研究表明：BEM?GEBT耦合計(jì)算方法結(jié)果與風(fēng)洞試驗(yàn)結(jié)果吻合較好；隨著預(yù)彎尺寸的增大，揮舞?擺振耦合顫振臨界風(fēng)速提高，顫振區(qū)間范圍基本一致；不同預(yù)彎尺寸葉片升力系數(shù)和俯仰力矩系數(shù)發(fā)散速率與位移發(fā)散速率呈正相關(guān)，平均風(fēng)壓曲線在預(yù)彎3～4 m范圍內(nèi)出現(xiàn)顯著變化，揮舞?擺振耦合效應(yīng)大于揮舞?扭轉(zhuǎn)耦合效應(yīng)，其顫振耦合頻率由一階揮舞頻率主導(dǎo)。

關(guān)鍵詞： 超長(zhǎng)柔性葉片； 預(yù)彎效應(yīng)； 風(fēng)洞試驗(yàn)； 氣動(dòng)彈性耦合； 顫振失穩(wěn)

中圖分類號(hào)： TK83 文獻(xiàn)標(biāo)志碼： A 文章編號(hào)： 1004-4523（2024）06-0986-11

DOI：10.16385/j.cnki.issn.1004-4523.2024.06.009

引 言

預(yù)彎［1］作為大兆瓦級(jí)風(fēng)力機(jī)葉片設(shè)計(jì)的關(guān)鍵環(huán)節(jié)，是防止超長(zhǎng)柔性葉片出現(xiàn)掃塔事故的主要措施。預(yù)彎改變了葉片氣動(dòng)外形和結(jié)構(gòu)特性，影響其氣動(dòng)力、慣性力和彈性力的相互作用和不同模態(tài)間的耦合振動(dòng)，易造成葉片氣動(dòng)彈性耦合顫振失穩(wěn)?，F(xiàn)有預(yù)彎設(shè)計(jì)大多采用靜氣彈設(shè)計(jì)方法［2?3］，未考慮動(dòng)氣彈耦合影響。然而在強(qiáng)/臺(tái)風(fēng)作用下，風(fēng)力機(jī)葉片因預(yù)彎效應(yīng)發(fā)生氣彈失穩(wěn)破壞的事故［4?7］屢見不鮮，如2008年臺(tái)風(fēng)“薔薇”導(dǎo)致風(fēng)力機(jī)葉片發(fā)生嚴(yán)重破壞，2013年臺(tái)風(fēng)“鲇魚”造成風(fēng)力機(jī)葉片折斷等。系統(tǒng)研究預(yù)彎尺寸效應(yīng)對(duì)葉片氣彈失穩(wěn)的發(fā)生機(jī)制具有重要意義。

目前針對(duì)風(fēng)力機(jī)葉片預(yù)彎效應(yīng)的研究大多采用數(shù)值模擬方法，主要內(nèi)容涉及葉片穩(wěn)定性能［8］、預(yù)彎靜力響應(yīng)［9?10］、機(jī)組輸出功率［3，11］等。研究發(fā)現(xiàn)，預(yù)彎有助于葉片扭轉(zhuǎn)載荷的降載和穩(wěn)定，預(yù)彎葉片截面的扭轉(zhuǎn)變形則使風(fēng)力機(jī)的實(shí)際氣動(dòng)功率減小，其軸向彈性變形有利于實(shí)際功率的提高。但靜氣彈分析方法忽略了結(jié)構(gòu)附加變形效應(yīng)、結(jié)構(gòu)運(yùn)動(dòng)位移和速度與氣動(dòng)力耦合效應(yīng)，很難精確反映預(yù)彎效應(yīng)的顫振失穩(wěn)影響。氣彈風(fēng)洞試驗(yàn)是研究顫振預(yù)彎效應(yīng)最有效的手段之一，但由于風(fēng)力機(jī)葉片翼型不規(guī)則，其截面、剛度、剪心等沿展長(zhǎng)不規(guī)則分布使得氣彈模型設(shè)計(jì)難，又因其大縮尺比帶來的測(cè)點(diǎn)布置難、采集干擾性強(qiáng)和測(cè)量精度低等試驗(yàn)因素，導(dǎo)致國內(nèi)外缺少風(fēng)力機(jī)葉片的三維氣彈風(fēng)洞試驗(yàn)?，F(xiàn)主要采用節(jié)段翼型［12?15］進(jìn)行測(cè)壓、測(cè)力試驗(yàn)，探討二維翼型測(cè)壓、測(cè)力的一致性，但二維翼型試驗(yàn)難以真實(shí)反映三維葉片預(yù)彎效應(yīng)對(duì)于顫振的影響。

鑒于此，本文以NREL?15 MW風(fēng)力機(jī)超長(zhǎng)柔性葉片為研究對(duì)象，首先基于主梁剛度等效原則設(shè)計(jì)超長(zhǎng)柔性葉片三維氣彈模型，并采用高速攝像技術(shù)和高頻六分量天平進(jìn)行同步測(cè)振、測(cè)力氣彈風(fēng)洞試驗(yàn)，提煉預(yù)彎尺寸對(duì)顫振區(qū)間及臨界風(fēng)速的組合規(guī)律；然后基于BEM?GEBT耦合計(jì)算方法研究預(yù)彎尺寸對(duì)葉片顫振性能的影響，最后揭示預(yù)彎尺寸對(duì)氣彈耦合模態(tài)的影響機(jī)理，研究結(jié)論可為超長(zhǎng)柔性風(fēng)力機(jī)葉片預(yù)彎設(shè)計(jì)提供參考依據(jù)。

1 氣彈模型風(fēng)洞試驗(yàn)

1.1 氣彈模型設(shè)計(jì)

以美國國家可再生能源實(shí)驗(yàn)室NREL?15 MW風(fēng)力機(jī)葉片［16］為研究對(duì)象，該葉片全長(zhǎng)117 m，葉片設(shè)計(jì)采用DTU FFA?W3翼型族，葉片根部直徑為5.2 m，在27.2 m（23.3%）處弦長(zhǎng)最大為5.77 m，單葉片質(zhì)量為65.252 t，葉片質(zhì)心為26.8 m，最大葉尖預(yù)彎位移量為4 m。葉片具體幾何參數(shù)如表1所示。表1中，r為以葉根處為起點(diǎn)，沿葉根到葉尖方向的實(shí)際位置；R為葉片總長(zhǎng)度；x為截面槳距中心位置；c為截面弦長(zhǎng)。

氣動(dòng)外形相似和結(jié)構(gòu)動(dòng)力學(xué)相似是氣彈模型設(shè)計(jì)的基本原則，風(fēng)洞試驗(yàn)葉片氣彈模型需要模擬幾何尺寸、風(fēng)場(chǎng)特性和氣動(dòng)彈性特性。氣彈相似條件物理量可采用無量綱參數(shù)來表示，如Reynolds數(shù)、Froude數(shù)、Cauchy數(shù)、密度比、阻尼比等。而風(fēng)力機(jī)葉片所在流場(chǎng)空氣為低速、不可壓縮、牛頓黏性流，其流體運(yùn)動(dòng)方程與結(jié)構(gòu)運(yùn)動(dòng)方程分別為：

（1）

（2）

式中 u為流體運(yùn)動(dòng)廣義速度；f為流體廣義外力；ρ為空氣密度；P為壓強(qiáng)；為空氣動(dòng)力黏度；x為坐標(biāo)系主軸；ω為鎖頻風(fēng)振頻率；M為廣義質(zhì)量；K為廣義剛度；g為阻尼系數(shù)；V為流場(chǎng)速度；b1為參考長(zhǎng)度；A為廣義空氣動(dòng)力系數(shù)；q為廣義坐標(biāo)。

引用符號(hào)λX表示模型的Xm同實(shí)物的Xa之比，即λX=Xm/Xa，可以用實(shí)物參數(shù)寫為：

（3）

（4）

對(duì)流體運(yùn)動(dòng)方程所有項(xiàng)乘以λl/λ，對(duì)結(jié)構(gòu)運(yùn)動(dòng)方程考慮減縮頻率相同，即λbλω/λV=1。為保證原型和模型流體運(yùn)動(dòng)與結(jié)構(gòu)運(yùn)動(dòng)的相似性，物理量的比值需滿足下式：

（5）

（6）

考慮風(fēng)洞試驗(yàn)條件限制，不易保證Re數(shù)相等，且模型與原型Re數(shù)均處于自模區(qū)［17］，不是決定性相似準(zhǔn)則，因此忽略Re數(shù)相等準(zhǔn)則。根據(jù)風(fēng)洞阻塞率要求，預(yù)彎與未預(yù)彎氣彈模型縮尺比均選為1∶70，其中，未預(yù)彎氣彈模型沿展長(zhǎng)預(yù)彎尺寸設(shè)為零，縮尺后的模型葉片底部最大直徑為0.074 m，高度為1.67 m，滿足阻塞率要求。因添加質(zhì)量塊和安裝葉片底座連接件產(chǎn)生了附加質(zhì)量，使得模型質(zhì)量與目標(biāo)質(zhì)量偏差稍大，模型動(dòng)力相似通過控制試驗(yàn)?zāi)Ｐ偷墓逃姓駝?dòng)特性來保證。其余無量綱參數(shù)由相似準(zhǔn)則推算得出，葉片氣彈模型主要相似參數(shù)如表2所示。

結(jié)構(gòu)動(dòng)力特性是進(jìn)行結(jié)構(gòu)動(dòng)力試驗(yàn)和動(dòng)力響應(yīng)分析的基礎(chǔ)，氣彈模型需考慮結(jié)構(gòu)頻率分布特點(diǎn)。為了便于調(diào)節(jié)葉片揮舞剛度、擺振剛度和扭轉(zhuǎn)剛度，基于變分漸進(jìn)梁截面法（VABS）提出等效梁截面法進(jìn)行風(fēng)力機(jī)葉片的氣彈模型制作，葉片通過變截面異形梁提供全部剛度，每段梁截面具體尺寸根據(jù)模型控制截面的目標(biāo)剛度設(shè)計(jì)，梁截面尺寸由葉根到葉尖逐漸減小，氣彈模型主梁理論剛度與實(shí)際剛度對(duì)比曲線及主梁截面形狀如圖1所示。

氣彈模型的氣動(dòng)外形采用“維形框段”結(jié)構(gòu)形式，維形框段保證模型的氣動(dòng)外形，框段和主梁連接方式采用單點(diǎn)連接。為保證框段不對(duì)整體結(jié)構(gòu)提供額外剛度，框段采用分段式打印，相鄰框段間有3 mm間隙，分段縫隙采用高密泡沫進(jìn)行填充，前后緣和檁條之間的鏤空區(qū)域采用輕質(zhì)木片填充打磨，以保證氣動(dòng)外形封閉。風(fēng)力機(jī)葉片整體采用聚酰胺纖維（未來7500尼龍）分段式3D打印，重心和轉(zhuǎn)動(dòng)慣量通過質(zhì)量塊調(diào)整配重進(jìn)行模擬。預(yù)彎與未預(yù)彎氣彈模型整體設(shè)計(jì)制作示意圖如圖2所示。

1.2 動(dòng)力特性分析

為驗(yàn)證氣彈模型與原型結(jié)構(gòu)的相似性，基于錘擊法測(cè)量氣彈模型固有頻率，采用隨機(jī)減量法［18］處理得到結(jié)構(gòu)模態(tài)參數(shù)。表3為氣彈模型與原型結(jié)構(gòu)動(dòng)力特性對(duì)比分析列表，其中NREL?15 MW研究報(bào)告［16］僅給出結(jié)構(gòu)前2階模態(tài)頻率。分析發(fā)現(xiàn)：預(yù)彎氣彈數(shù)值模型和預(yù)彎氣彈模型與原型結(jié)構(gòu)的前2階模態(tài)固有頻率基本吻合，誤差均在10%以內(nèi)，表明氣彈模型的動(dòng)力學(xué)特性與原型結(jié)構(gòu)匹配較好，保證了氣彈模型的顫振特性與原型結(jié)構(gòu)的相似精度；未預(yù)彎氣彈模型前2階頻率明顯大于預(yù)彎氣彈模型。

1.3 風(fēng)洞測(cè)量系統(tǒng)

風(fēng)洞試驗(yàn)段長(zhǎng)20 m，高2 m，寬2.5 m，最大試驗(yàn)風(fēng)速為30 m/s，風(fēng)洞上游來流風(fēng)速通過眼鏡蛇三維脈動(dòng)風(fēng)速測(cè)量?jī)x測(cè)量。葉片模型豎直放置在風(fēng)洞的轉(zhuǎn)盤中心處，通過連接件與高頻六分量天平固定，其底端與風(fēng)洞轉(zhuǎn)盤固結(jié)，模型可以隨轉(zhuǎn)盤360°轉(zhuǎn)動(dòng)，通過采集系統(tǒng)實(shí)時(shí)測(cè)量根部六分力。葉尖位移采集通過補(bǔ)光燈照射葉尖反光點(diǎn)，并使用高速攝像系統(tǒng)記錄葉尖反光點(diǎn)移動(dòng)軌跡，高速攝像采樣頻率為330 Hz，每個(gè)試驗(yàn)工況采集60 s葉尖運(yùn)動(dòng)高幀視頻，再使用編程處理高速攝像視頻，得到葉尖揮舞和擺振方向位移。風(fēng)洞氣彈模型測(cè)量系統(tǒng)示意圖如圖3所示。

鑒于風(fēng)洞試驗(yàn)中葉片槳距角不能主動(dòng)調(diào)節(jié)，又考慮到實(shí)際風(fēng)向的突變性［19?20］，葉片的風(fēng)攻角區(qū)間為0°～360°，以每10°為間隔共設(shè)置36個(gè)槳距角，每個(gè)槳距角測(cè)量7個(gè)風(fēng)速工況，逐級(jí)加載，在風(fēng)振響應(yīng)顯著的工況下每隔0.5°進(jìn)行內(nèi)插，提取風(fēng)振響應(yīng)敏感槳距角。定義葉片正常運(yùn)行狀態(tài)下的轉(zhuǎn)盤為0°，轉(zhuǎn)盤以面向來流方向順時(shí)針旋轉(zhuǎn)為正，順風(fēng)向?yàn)槿~片揮舞方向，垂直順風(fēng)向?yàn)槿~片擺振方向，風(fēng)洞試驗(yàn)來流風(fēng)向角定義示意圖如圖4所示。

1.4 試驗(yàn)結(jié)果分析

圖5給出了全槳距角下風(fēng)速為8.7 m/s（實(shí)際風(fēng)速為72.8 m/s）時(shí)，預(yù)彎與未預(yù)彎氣彈模型的葉尖位移均方根對(duì)比圖。由圖5可知，在槳距角0°～180°范圍內(nèi)，預(yù)彎與未預(yù)彎氣彈模型顫振區(qū)間分別為93°～96°和92°～94°；在槳距角180°～360°范圍內(nèi)，預(yù)彎與未預(yù)彎氣彈模型顫振區(qū)間分別為284°～287°和282°～285.5°；葉尖揮舞位移均方根顯著大于葉尖擺振位移均方根。

將風(fēng)速從4 m/s開始每隔0.1 m/s逐級(jí)加載，預(yù)彎與未預(yù)彎氣彈模型的顫振臨界風(fēng)速結(jié)果如圖6所示。對(duì)比發(fā)現(xiàn)：兩葉片在顫振區(qū)間內(nèi)顫振臨界風(fēng)速均呈現(xiàn)先減小后增大的趨勢(shì)；在槳距角0°～180°范圍內(nèi)，預(yù)彎與未預(yù)彎氣彈模型最小顫振臨界風(fēng)速分別為5.4和4.5 m/s，最大顫振臨界風(fēng)速分別為8.7和5.8 m/s；在槳距角180°～360°范圍內(nèi)，預(yù)彎與未預(yù)彎氣彈模型最小顫振臨界風(fēng)速分別為5.4和4.1 m/s，最大顫振臨界風(fēng)速分別為8.7和6.4 m/s；在二者相交的顫振區(qū)間內(nèi)，預(yù)彎模型的顫振臨界風(fēng)速比未預(yù)彎模型的大。

圖7給出了顫振區(qū)間臨界風(fēng)速下葉尖位移振動(dòng)幅值隨槳距角變化的曲線。對(duì)比顫振區(qū)間臨界風(fēng)速下葉尖位移振動(dòng)幅值的變化可知：預(yù)彎與未預(yù)彎氣彈模型的揮舞位移明顯大于擺振位移，揮舞位移最大值和最小值的絕對(duì)值基本相等，擺振位移亦呈現(xiàn)相同特點(diǎn)。結(jié)合圖7分析可知：預(yù)彎與未預(yù)彎氣彈模型在顫振槳距角0°～180°范圍內(nèi)，揮舞/擺振振幅與顫振臨界風(fēng)速變化趨勢(shì)一致；在顫振槳距角180°～360°范圍內(nèi)，揮舞/擺振振幅與顫振臨界風(fēng)速呈非線性關(guān)系。

圖8給出了風(fēng)速為8.7 m/s，槳距角分別為94°，280°和285°時(shí)，預(yù)彎與未預(yù)彎氣彈模型的揮舞/擺振葉尖位移頻譜曲線。由圖8可知：預(yù)彎與未預(yù)彎氣彈模型發(fā)生顫振時(shí)，揮舞位移主導(dǎo)頻率與擺振位移主導(dǎo)頻率基本相等，耦合主導(dǎo)頻率與各自的一階揮舞頻率（4.682和5.463 Hz）基本一致；預(yù)彎與未預(yù)彎氣彈模型未發(fā)生顫振時(shí)，葉尖位移頻譜曲線呈平緩波動(dòng)趨勢(shì)，均沒有主導(dǎo)頻率。

2 BEM-GEBT耦合計(jì)算方法

葉片發(fā)生大幅振動(dòng)時(shí)結(jié)構(gòu)變形較大，超長(zhǎng)預(yù)彎葉片的柔性變形具有明顯的幾何非線性特點(diǎn)，且沿展長(zhǎng)方向各翼型攻角隨時(shí)間發(fā)生周期性變化，從而影響翼型氣動(dòng)力。為了研究預(yù)彎尺寸的顫振機(jī)理，結(jié)合修正的葉素動(dòng)量理論與幾何精確梁理論分析葉片非定常效應(yīng)。

2.1 BEM-GEBT耦合方法

2.1.1 氣動(dòng)控制方程

基于葉素動(dòng)量理論計(jì)算風(fēng)力機(jī)葉片氣動(dòng)力，假設(shè)風(fēng)輪所受氣動(dòng)載荷沿葉片展向均勻分布，且作用在葉片截面弦上距離前緣1/4即氣動(dòng)中心處。為了獲得較為準(zhǔn)確的計(jì)算結(jié)果，采用Prandtl葉尖損失修正［21］、Pitt/Peters模型［22］、動(dòng)態(tài)失速模型［23］對(duì)該算法進(jìn)行修正，進(jìn)行非線性迭代求解。

由葉素定理和動(dòng)量定理可得，作用在風(fēng)輪半徑r1處的微元環(huán)所受推力和轉(zhuǎn)矩分別為：

（7）

（8）

式中 B為葉片數(shù)；V0為自由來流速度；Cn和Ct分別為法向力系數(shù)和切向力系數(shù)；?為入流角；askew為修正后的軸向誘導(dǎo)因子；a′為切向誘導(dǎo)因子；F為葉尖損失修正因子。

2.1.2 結(jié)構(gòu)控制方程

幾何精確梁理論［24］適用于幾何非線性及大撓度變形結(jié)構(gòu)，具有伸展、彎曲、剪切和扭轉(zhuǎn)自由度，使用完整的6×6質(zhì)量和剛度矩陣，其運(yùn)動(dòng)控制方程為：

（9）

（10）

式中 和分別為慣性系中的線動(dòng)量和角動(dòng)量；和分別為梁的截面力和力矩；為梁參考線上一點(diǎn)的一維偏移量；為沿著梁參考線上一點(diǎn)的方向向量；和分別為作用在梁結(jié)構(gòu)上的分布力和力矩；其中，變量下方加下劃線標(biāo)記表示向量或相應(yīng)矩陣的維度，加雙下劃線標(biāo)記表示張量；另外，上角標(biāo)“'”表示對(duì)梁的軸線求導(dǎo)數(shù)，變量上方加“?”表示對(duì)時(shí)間求導(dǎo)數(shù)，變量上方加“～”表示給定向量的反對(duì)稱張量。

葉片的阻尼力與應(yīng)變率成正比，針對(duì)6個(gè)自由度（3個(gè)平移，3個(gè)旋轉(zhuǎn)）給定葉片6個(gè)阻尼系數(shù)。通過每個(gè)截面的剛度矩陣和阻尼系數(shù)矩陣以及應(yīng)變率可以得到阻尼力為：

（11）

式中 為阻尼力；為應(yīng)變率；為6×6剛度矩陣；為阻尼系數(shù)矩陣。

2.2 OpenFAST二次開發(fā)建模

為進(jìn)一步研究不同預(yù)彎尺寸對(duì)顫振性能的影響，設(shè)計(jì)了4種預(yù)彎葉片，其尺寸分別為0，3，4和5 m，葉片翼型選用與NREL?15 MW風(fēng)力機(jī)葉片相同數(shù)據(jù)，基于積疊線形狀設(shè)計(jì)方法［2］，采用與葉片展向線性和二次曲線結(jié)合的形式對(duì)葉片進(jìn)行預(yù)彎，預(yù)彎起于葉片根部。預(yù)彎表達(dá)式如下式所示：

（12）

式中 為沿展向的預(yù)彎位移；R0為葉片長(zhǎng)度；a和b分別為線性和二次項(xiàng)系數(shù)。

綜合考慮風(fēng)力機(jī)輸出功率隨葉片預(yù)彎尺寸增大呈非線性減小［11］及葉片掃塔因素，預(yù)彎尺寸參數(shù)分別取a=2.5，b=5.5（預(yù)彎3 m）；a=3.3，b=7.3（預(yù)彎4 m）；a=4.1，b=9.1（預(yù)彎5 m）進(jìn)行設(shè)計(jì)。

根據(jù)公式（12）得到4種葉片預(yù)彎曲線，并將葉片與NREL?15 MW原型結(jié)構(gòu)進(jìn)行對(duì)比，結(jié)果如圖9所示。由圖9可知：基于積疊線形狀設(shè)計(jì)方法得到的預(yù)彎4 m曲線與NREL?15 MW原型結(jié)構(gòu)擬合較好。

對(duì)4種預(yù)彎葉片進(jìn)行模態(tài)計(jì)算，表4給出了模態(tài)計(jì)算結(jié)果與NREL?15 MW原型結(jié)構(gòu)（預(yù)彎4 m）前2階固有頻率［16］的對(duì)比。對(duì)比4種預(yù)彎葉片固有頻率可知：所設(shè)計(jì)的預(yù)彎葉片具有較高的精度；葉片前4階振型以揮舞或擺振方向?yàn)橹?；?階振型以揮舞方向?yàn)橹鳎晕殡S扭轉(zhuǎn)方向振動(dòng)；第6階振型主要表現(xiàn)為扭轉(zhuǎn)方向振動(dòng)；葉片前3階固有頻率隨預(yù)彎尺寸的增大均有所降低，且一階揮舞頻率與一階擺振頻率偏離程度逐漸增大。

2.3 有效性驗(yàn)證

為驗(yàn)證數(shù)值模擬與風(fēng)洞試驗(yàn)數(shù)據(jù)結(jié)果的一致性，選取風(fēng)洞試驗(yàn)風(fēng)速為8.7 m/s及槳距角為94°的工況結(jié)果進(jìn)行驗(yàn)證對(duì)比。考慮數(shù)值模擬采用足尺模型，風(fēng)洞試驗(yàn)相應(yīng)數(shù)據(jù)按縮尺比進(jìn)行調(diào)整，其足尺模型對(duì)應(yīng)風(fēng)速為72.8 m/s，葉尖位移和葉根反力動(dòng)態(tài)響應(yīng)頻譜曲線如圖10所示。對(duì)比葉尖位移和葉根反力動(dòng)態(tài)響應(yīng)頻譜曲線可知：數(shù)值頻率分別為0.578和0.560 Hz，與一階揮舞頻率（0.555 Hz）基本對(duì)應(yīng)，風(fēng)洞試驗(yàn)頻率分別為4.900和4.859 Hz，按照縮尺比換算全尺寸模型為0.586和0.581 Hz；數(shù)值揮舞/擺振位移及揮舞方向葉根反力的主導(dǎo)頻率與風(fēng)洞試驗(yàn)結(jié)果基本一致，誤差均小于5%，表明BEM?GEBT耦合計(jì)算方法結(jié)果與風(fēng)洞試驗(yàn)結(jié)果吻合較好。

3 預(yù)彎效應(yīng)分析

3.1 顫振臨界風(fēng)速

圖11給出了由BEM?GEBT耦合計(jì)算方法計(jì)算來流風(fēng)速為72.8 m/s時(shí)，4種預(yù)彎葉片的顫振區(qū)間及顫振臨界風(fēng)速。由圖11可知：4種預(yù)彎葉片在槳距角0°～180°和180°～360°范圍內(nèi)均會(huì)發(fā)生顫振，定義槳距角0°～180°范圍內(nèi)的顫振區(qū)間為顫振Ⅰ區(qū)間，180°～360°范圍內(nèi)的顫振區(qū)間為顫振Ⅱ區(qū)間。4種預(yù)彎葉片的顫振區(qū)間有較大重合部分，顫振Ⅰ，Ⅱ區(qū)間的臨界風(fēng)速均呈拋物線趨勢(shì)，初始顫振槳距角隨預(yù)彎尺寸的增大呈后移現(xiàn)象，與風(fēng)洞試驗(yàn)得到的規(guī)律一致；同一顫振槳距角下，顫振臨界風(fēng)速隨預(yù)彎尺寸的增大而增大；隨著預(yù)彎尺寸改變，顫振區(qū)間未發(fā)生明顯變化。

3.2 氣動(dòng)力特性

為分析預(yù)彎尺寸對(duì)葉片顫振氣動(dòng)力的影響，選取未預(yù)彎葉片和預(yù)彎4 m葉片兩種典型工況進(jìn)行分析。圖12給出了來流風(fēng)速為72.8 m/s及槳距角為94°時(shí)，兩種工況葉片的葉尖截面升力系數(shù)Cl和俯仰力矩系數(shù)CT的歸一化曲線。由圖12可知：升力系數(shù)隨時(shí)間變化呈發(fā)散狀態(tài)，葉尖截面翼型攻角隨葉片大幅發(fā)散呈現(xiàn)周期性變化；升力系數(shù)與俯仰力矩系數(shù)幅值增長(zhǎng)速度隨預(yù)彎尺寸的增大而減緩；同一預(yù)彎尺寸的升力系數(shù)曲線與俯仰力矩系數(shù)曲線波動(dòng)相位一致，可以看出不同預(yù)彎尺寸的升力特性主要受葉片扭轉(zhuǎn)的影響。

圖13給出了來流風(fēng)速為72.8 m/s及槳距角為94°時(shí)，4種預(yù)彎葉片沿展長(zhǎng)平均風(fēng)壓變化的曲線。由圖13可知：4種預(yù)彎葉片平均風(fēng)壓曲線變化趨勢(shì)基本一致，沿展長(zhǎng)呈現(xiàn)先減小后增大的趨勢(shì)，且平均風(fēng)壓變化速率隨預(yù)彎尺寸的增大而減?。辉谌~片展長(zhǎng)0～47 m范圍內(nèi)，平均風(fēng)壓相差不明顯，展長(zhǎng)超過47 m后，平均風(fēng)壓隨預(yù)彎尺寸的增大呈減小趨勢(shì)且在3～4 m范圍內(nèi)減小趨勢(shì)顯著；沿展長(zhǎng)方向平均風(fēng)壓呈非線性增長(zhǎng)且隨預(yù)彎尺寸的增大趨于平緩分布。

3.3 氣彈耦合模態(tài)

圖14給出了4種預(yù)彎葉片的典型響應(yīng)對(duì)比曲線。由位移時(shí)程曲線可知：顫振位移發(fā)展階段隨預(yù)彎尺寸的增加而逐漸延長(zhǎng)，但位移不會(huì)無限發(fā)散，最終發(fā)展為鎖幅振動(dòng)且4種預(yù)彎葉片振動(dòng)幅值基本一致。由頻譜對(duì)比曲線可知：不同預(yù)彎尺寸葉片揮舞、擺振和扭轉(zhuǎn)的主導(dǎo)頻率基本相等并對(duì)應(yīng)各自一階揮舞頻率，即葉尖顫振變形主要受一階揮舞振動(dòng)的影響。由振動(dòng)幅值可知：耦合頻率處的振動(dòng)幅值隨預(yù)彎尺寸的增大而減小，葉尖運(yùn)動(dòng)能量隨預(yù)彎尺寸的增大呈減小趨勢(shì)，4種預(yù)彎葉片的揮舞?擺振耦合效應(yīng)大于揮舞?扭轉(zhuǎn)耦合效應(yīng)。

圖15給出了4種預(yù)彎葉片的揮舞?擺振雙軸同步振動(dòng)葉尖位移軌跡曲線。由圖15可知：擺振方向位移隨預(yù)彎尺寸的增大呈現(xiàn)先減小后增大再減小的趨勢(shì)，揮舞方向位移隨預(yù)彎尺寸的增大出現(xiàn)先增大再減小的趨勢(shì)；4種預(yù)彎葉片雙軸軌跡圖由初始振動(dòng)逐步增大，發(fā)展到一定階段達(dá)到穩(wěn)態(tài)，隨后進(jìn)入鎖幅振動(dòng)狀態(tài)，雙軸位移軌跡表現(xiàn)為類橢圓閉合曲線。

4 結(jié) 論

（1）基于主梁剛度等效原則設(shè)計(jì)了三維葉片氣彈模型，發(fā)展了基于修正的葉素動(dòng)量和幾何精確梁理論耦合的顫振分析方法，與風(fēng)洞試驗(yàn)結(jié)果對(duì)比，表明BEM?GEBT耦合計(jì)算方法能有效模擬風(fēng)力機(jī)葉片的顫振預(yù)彎效應(yīng)，最大誤差為5%。

（2）風(fēng)力機(jī)葉片固有頻率隨預(yù)彎尺寸的增大逐漸減小，且一階揮舞頻率與一階擺振頻率偏離程度呈加大趨勢(shì)；顫振區(qū)間隨預(yù)彎尺寸的增大呈現(xiàn)小幅度后移趨勢(shì)，相同槳距角下顫振臨界風(fēng)速亦隨預(yù)彎尺寸的增大而增大，顫振區(qū)間范圍基本一致。

（3）升力系數(shù)與俯仰力矩系數(shù)發(fā)散速率隨預(yù)彎尺寸的增大而減緩，并與位移發(fā)散速率呈正相關(guān)，且不同預(yù)彎尺寸的升力特性主要受葉片扭轉(zhuǎn)的影響；4種預(yù)彎葉片平均風(fēng)壓曲線變化趨勢(shì)基本一致，沿展長(zhǎng)呈現(xiàn)先減小后增大的趨勢(shì)。

（4）葉尖運(yùn)動(dòng)能量積累隨預(yù)彎尺寸的增大呈減小趨勢(shì)；不同預(yù)彎超長(zhǎng)柔性葉片的揮舞?擺振耦合效應(yīng)大于揮舞?扭轉(zhuǎn)耦合效應(yīng)，其顫振耦合頻率由一階揮舞頻率主導(dǎo)，雙軸位移軌跡表現(xiàn)為規(guī)則的類橢圓閉合曲線。

參考文獻(xiàn)：

［1］ OLESEN E E， KOEFOED M S， PEDERSEN S H， et al. Method of manufacturing pre-bent wind turbine blades： US13176419［P］. 2012-07-10.

［2］ HANSEN M H. Aeroelastic properties of backward swept blades［C］. 49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition. Orlando， Florida， 2011： 260.

［3］ 郭小鋒， 陳進(jìn)， 李松林. 風(fēng)力機(jī)葉片預(yù)彎設(shè)計(jì)機(jī)理與方法［J］. 東北大學(xué)學(xué)報(bào)（自然科學(xué)版）， 2016， 37（6）： 851-855.

GUO X F， CHEN J， LI S L. Mechanism and method of pre-bending design for the wind turbine blades［J］. Journal of Northeastern University （Natural Science）， 2016， 37（6）： 851-855.

［4］ CHOU J S， CHIU C K， HUANG I K， et al. Failure analysis of wind turbine blade under critical wind loads［J］. Engineering Failure Analysis， 2013， 27： 99-118.

［5］ 王浩， 柯世堂， 王同光. 臺(tái)風(fēng)過境全過程大型風(fēng)力機(jī)風(fēng)荷載特性［J］. 空氣動(dòng)力學(xué)學(xué)報(bào)， 2020， 38（5）： 915-923.

WANG H， KE S T， WANG T G. Wind loads characteristic of large wind turbine considering typhoon transit process［J］ Acta Aerodynamica Sinica， 2020， 38（5）： 915-923.

［6］ 柯世堂， 徐璐， 王同光， 等. 考慮停機(jī)位置風(fēng)力機(jī)塔架-葉片耦合體系風(fēng)致穩(wěn)定性能分析［J］. 太陽能學(xué)報(bào)， 2020， 41（1）： 264-272.

KE S T， XU L， WANG T G， et al. Analysis of wind-induced stability tower-blade coupling system of wind turbine considering stopping positions［J］. Acta Energiae Solaris Sinica， 2020， 41（1）： 264-272.

［7］ CHEN X， XU J Z. Structural failure analysis of wind turbines impacted by super typhoon Usagi［J］. Engineering Failure Analysis， 2016， 60： 391-404.

［8］ BAZILEVS Y， HSU M C， KIENDL J， et al. A computational procedure for prebending of wind turbine blades［J］. International Journal for Numerical Methods in Engineering， 2012， 89（3）： 323-336.

［9］ RIZIOTIS V A， VOUTSINAS S G， MANOLAS D I， et al. Aeroelastic analysis of pre-curved rotor blades［C］. Proceedings of the European Wind Energy Conference and Exhibition 2010. Warsaw， Poland， 2010： 1173-1202.

［10］ 田德， 羅濤， 梁俊宇， 等. 風(fēng)電機(jī)組預(yù)彎葉片模型的有限元分析［J］. 太陽能學(xué)報(bào)， 2017， 38（6）： 1482-1488.

TIAN D， LUO T， LIANG J Y， et al. Finite element analysis of pre-bending blade model for wind turbine generation system［J］. Acta Energiae Solaris Sinica， 2017， 38（6）： 1482-1488.

［11］ 郭婷婷， 吳殿文， 王成蔭， 等. 風(fēng)力發(fā)電機(jī)葉片預(yù)彎設(shè)計(jì)及其數(shù)值研究［J］. 動(dòng)力工程學(xué)報(bào)， 2010， 30（6）： 450-455.

GUO T T， WU D W， WANG C Y， et al. Pre-bend design and numerical simulation of wind turbo-generator［J］. Journal of Chinese Society of Power Engineering， 2010， 30（6）： 450-455.

［12］ 程鵬. 水平軸風(fēng)力機(jī)葉片氣動(dòng)彈性的風(fēng)洞試驗(yàn)研究［D］. 哈爾濱： 哈爾濱工業(yè)大學(xué)， 2009.

CHENG P. Experimental study of aeroelasticity for the blades of horizontal axis wind turbines［D］. Harbin： Harbin Institute of Technology， 2009.

［13］ RAZAK N A， ANDRIANNE T， DIMITRIADIS G. Flutter and stall flutter of a rectangular wing in a wind tunnel［J］. AIAA Journal， 2011， 49（10）： 2258-2271.

［14］ DIMITRIADIS G， LI J. Bifurcation behavior of airfoil undergoing stall flutter oscillations in low-speed wind tunnel［J］. AIAA Journal， 2009， 47（11）： 2577-2596.

［15］ 張超宇. 后掠葉段氣動(dòng)及氣彈特性的風(fēng)洞實(shí)驗(yàn)研究［D］. 北京： 中國科學(xué)院大學(xué)（中國科學(xué)院工程熱物理研究所）， 2020.

ZHANG C Y. Study on aerodynamic and aeroelastic characteristic of swept-back airfoil model based on wind tunnel experiments［D］. Beijing： University of Chinese Academy of Sciences （Institute of Engineering Thermophysics， Chinese Academy of Sciences）， 2020.

［16］ GAERTNER E， RINKER J， SETHURAMAN L， et al. IEA wind TCP task 37： definition of the IEA wind 15-megawatt offshore reference wind turbine： NREL/TP-5000-75698［R］. Golden， C.O.， United States： National Renewable Energy Laboratory （NREL）， 2020.

［17］ 黃本才， 汪叢軍. 結(jié)構(gòu)抗風(fēng)分析原理及應(yīng)用［M］. 2版. 上海： 同濟(jì)大學(xué)出版社， 2008： 112-113.

HUANG B C， WANG C J. Principle and Application of Structural Wind Resistance Analysis［M］. 2nd ed. Shanghai： Tongji University Press， 2008： 112-113.

［18］ IBRAHIM S R. Random decrement technique for modal identification of structures［J］. Journal of Spacecraft and Rockets， 1977， 14（11）： 696-700.

［19］ 吳金城， 張容焱， 張秀芝. 海上風(fēng)電機(jī)的抗臺(tái)風(fēng)設(shè)計(jì)［J］. 中國工程科學(xué)， 2010， 12（11）： 25-31.

WU J C， ZHANG R Y， ZHANG X Z. Anti-typhoon design for offshore wind turbines［J］. Strategic Study of CAE， 2010， 12（11）： 25-31.

［20］ 宋麗莉， 毛慧琴， 錢光明， 等. 熱帶氣旋對(duì)風(fēng)力發(fā)電的影響分析［J］. 太陽能學(xué)報(bào)， 2006， 27（9）： 961-965.

SONG L L， MAO H Q， QIAN G M， et al. Analysis on the wind power by tropical cyclone［J］. Acta Energiae Solaris Sinica， 2006， 27（9）： 961-965.

［21］ SHEN W Z， MIKKELSEN R， S?RENSEN J N， et al. Tip loss corrections for wind turbine computations［J］. Wind Energy， 2005， 8（4）： 457-475.

［22］ NING S A， HAYMAN G， DAMIANI R， et al. Development and validation of a new blade element momentum skewed-wake model within AeroDyn［C］//33rd Wind Energy Symposium. Kissimmee， Florida， 2015： 0215.

［23］ LEISHMAN J G， BEDDOES T S. A semi-empirical model for dynamic stall［J］. Journal of the American Helicopter Society， 1989， 34（3）： 3-17.

［24］ JúNIOR C J F， CARDOZO A C P， JúNIOR V M， et al. Modeling wind turbine blades by geometrically-exact beam and shell elements： a comparative approach［J］. Engineering Structures， 2019， 180： 357-378.

Analysis of flutter effect in the pre?bending state of wind turbine blade

GAO Mu-en1， KE Shi-tang1，2， WU Hong-xin1，2， ZHANG Chun-wei1，2，

TIAN Wen-xin1， LU Man-man1

（1.Department of Civil and Airport Engineering， College of Civil Aviation， Nanjing University of Aeronautics and Astronautics， Nanjing 211106， China； 2.Jiangsu Key Laboratory of Hi-Tech Research for Wind Turbine Design，Nanjing University of Aeronautics and Astronautics， Nanjing 210016， China）

Abstract： Most wind turbine blade pre-bending designs use the static aeroelastic analysis method. This approach often overlooks the aeroelastic coupling instability caused by the interaction of blade aerodynamic force， inertial force and elastic force. This oversight is particularly significant when considering flutter performance of ultra-long flexible blades of around 100 meters. To analyze the influence of different pre-bending sizes on flutter critical state of blade， aeroelastic model of the blade was designed based on the stiffness equivalence principle of the main beam. Wind tunnel tests revealed differences between the flutter interval and the critical wind speed of two pre-bending blades of a 15 MW wind turbine. Further analysis was conducted on four pre-bending blades using the corrected Blade Element Momentum Theory-Geometrically Exact Beam Theory （BEM-GEBT） coupling calculation method. This analysis compared and analyzed the flutter critical wind speed， aerodynamic force distribution and displacement spectrum characteristics of blades with different pre-bending sizes， revealing the flutter coupling modal mechanism. The research shows that the results of BEM-GEBT coupling calculation method align well with those of wind tunnel test. As the pre-bending size increases， the flutter critical wind speed of flap-edge coupling increases， and the flutter interval range remains essentially the same. The divergence rates of lift coefficient and pitching moment coefficient of different pre-bending blades are positively correlated with the displacement divergence rate. The average wind pressure curve shows significant changes in the pre-bending range of 3～4 m. The flap-edge coupling effect is larger than the flap-torsion coupling effect， and the flutter coupling frequency is dominated by the first-order flapwise frequency.

Key words： ultra-long flexible blades；pre-bending effect；wind tunnel test；aeroelastic coupling；flutter instability

作者簡(jiǎn)介： 高沐恩（1997―），男，碩士。E-mail：gaomuen@nuaa.edu.cn。

通訊作者： 柯世堂（1982―），男，博士，教授。E-mail：keshitang@163.com。