收稿日期：2023-10-24

基金項目：國家重點研發(fā)計劃課題（2018YFB1501304）

通信作者：田 德（1958—），男，博士、教授、博士生導(dǎo)師，主要從事風(fēng)力發(fā)電系統(tǒng)理論和技術(shù)方面的研究。tdncepu@163.com

DOI：10.19912/j.0254-0096.tynxb.2023-1751 文章編號：0254-0096（2024）02-0198-08

摘 要：以IEA Wind 15 MW 參考風(fēng)電機(jī)組為研究對象，基于非線性時域氣彈耦合分析工具，研究風(fēng)剪切、偏航入流、質(zhì)量不平衡3種不平衡載荷單一或共同作用時，超大型風(fēng)電機(jī)組葉片的顫振極限。結(jié)果表明，單一不平衡載荷作用時，風(fēng)剪切和±15°以內(nèi)的偏航角會增大臨界顫振速度，而質(zhì)量不平衡會降低臨界顫振速度，且質(zhì)量偏重時顫振速度下降更加明顯。3種不平衡載荷共同作用時，臨界顫振速度的最大值和最小值都出現(xiàn)在剪切系數(shù)為0.3時；最大值總是在偏航角為非負(fù)值時取得，最小值在偏航角為[-20°]時取得；質(zhì)量不平衡對臨界顫振速度的影響與其單一作用時的規(guī)律一致。因此，在進(jìn)行氣彈穩(wěn)定性分析時，應(yīng)考慮多種風(fēng)輪不平衡載荷共同作用，尤其是葉片質(zhì)量偏重和偏航角度為[-20°]時這種最不利的情況。

關(guān)鍵詞：風(fēng)電機(jī)組葉片；顫振；穩(wěn)定性；氣動彈性；風(fēng)輪不平衡載荷

中圖分類號：TK83"""""""""nbsp;""""""""""""""""""""""" 文獻(xiàn)標(biāo)志碼：A

0 引 言

為降低風(fēng)電機(jī)組度電成本，風(fēng)電機(jī)組單機(jī)容量逐漸增大，裝備超長葉片的10 MW級超大型風(fēng)電機(jī)組將成為未來主流機(jī)型。據(jù)全球風(fēng)能理事會估計，2025年全球新增海上風(fēng)電機(jī)組的平均功率將達(dá)到10～12 MW［1］。另一方面，為突破平方-立方定律的限制，葉片的質(zhì)量和剛度會進(jìn)一步縮減，柔性變大。這使得葉片可能在運(yùn)行過程中有發(fā)生顫振的風(fēng)險。文獻(xiàn)［2］強(qiáng)調(diào)超長柔性葉片穩(wěn)定性分析的必要性，以確保顫振的安全裕度。葉片顫振是其在氣動力、彈性力和慣性力的耦合作用下發(fā)生的動態(tài)不穩(wěn)定［3］，而不平衡載荷可能會加劇這種不穩(wěn)定，甚至降低臨界顫振速度。因此，分析風(fēng)輪不平衡載荷下，風(fēng)電機(jī)組葉片的顫振特性，對風(fēng)電機(jī)組的安全運(yùn)行具有重要意義。

針對風(fēng)電機(jī)組葉片的顫振問題，文獻(xiàn)［4］考慮幾何非線性，研究3種風(fēng)電機(jī)組葉片：NREL 5 MW、SNL 61.5和SNL 100-00的顫振不穩(wěn)定性；文獻(xiàn)［5］使用多種氣彈工具分析了IWT 7.5 MW風(fēng)電機(jī)組失控狀態(tài)下的氣彈穩(wěn)定性；文獻(xiàn)［6］基于907自由度氣動-水動-控制-彈性風(fēng)電機(jī)組模型，研究DTU 10 MW海上風(fēng)電機(jī)組的顫振性能。

針對風(fēng)電機(jī)組不平衡載荷的問題，學(xué)者們主要關(guān)注于不平衡載荷的檢測和評估。文獻(xiàn)［7］提出一種檢測葉片質(zhì)量不平衡的方法，結(jié)果表明該方法在多個風(fēng)速和湍流強(qiáng)度下有效；文獻(xiàn)［8］研究風(fēng)剪切、湍流、偏航、塔影效應(yīng)、質(zhì)量不平衡、空氣動力不平衡對于葉片攻角和升力系數(shù)波動的影響；文獻(xiàn)［9-11］基于不同控制策略設(shè)計了獨立變槳控制器來抑制風(fēng)輪不平衡載荷；對于不平衡引起的顫振速度變化，文獻(xiàn)［12］發(fā)現(xiàn)湍流強(qiáng)度會降低臨界顫振速度。

上述研究中關(guān)于風(fēng)輪不平衡載荷下葉片顫振特性的研究較少，而且研究對象多為容量在10 MW及以下的風(fēng)電機(jī)組，葉片長度不超過100 m。隨著葉片長度的增加，載荷作用下的葉尖撓度將會更大，風(fēng)輪不平衡載荷作用更加明顯。因此，建立風(fēng)電機(jī)組葉片氣彈耦合模型，基于軟件HAWC2［13］，采用時域方法，分析風(fēng)剪切、偏航入流、質(zhì)量不平衡3種不平衡載荷單一作用或共同作用時，超大型風(fēng)電機(jī)組葉片臨界顫振速度的變化規(guī)律。

1 研究對象和研究工具

1.1 風(fēng)電機(jī)組介紹

風(fēng)電機(jī)組仿真模型為IEA Wind 15 MW Reference Wind Turbine（RWT）［14］，由美國國家可再生能源實驗室和丹麥技術(shù)大學(xué)共同開發(fā)，是目前世界上單機(jī)容量最大的參考型風(fēng)電機(jī)組［15］，表1中列出了風(fēng)電機(jī)組的關(guān)鍵參數(shù)。

1.2 研究工具

HAWC2［13］是由DTU開發(fā)并持續(xù)發(fā)展的非線性時域氣彈

工具，能處理風(fēng)電機(jī)組所有結(jié)構(gòu)構(gòu)件的大旋轉(zhuǎn)和大撓度，可用于模擬風(fēng)電機(jī)組葉片發(fā)生顫振時的動態(tài)響應(yīng)。

2 氣彈模型

2.1 氣動模型

2.1.1 極坐標(biāo)網(wǎng)格葉素動量理論

由于風(fēng)剪切、偏航入流和質(zhì)量不平衡引起的風(fēng)輪不平衡載荷，使得整個風(fēng)輪平面載荷均勻分布的假設(shè)不合理，原始的葉素動量理論不再適用。因此，引入極坐標(biāo)網(wǎng)格葉素動量理論，其基本思想是分別計算各個葉片的推力系數(shù)和轉(zhuǎn)矩系數(shù)，不在葉片上的網(wǎng)格點的推力系數(shù)和轉(zhuǎn)矩系數(shù)由網(wǎng)格點和相鄰兩個葉片之間的方位角加權(quán)得到，示意圖如圖1所示。在文獻(xiàn)［16］中詳細(xì)介紹了極網(wǎng)格葉素動量理論具體計算過程。此外，為準(zhǔn)確計算氣動載荷，進(jìn)行顫振分析，采用Prandtl葉尖損失模型［17］和Beddoes-Leishman模型［18］進(jìn)行修正。

2.1.2 偏航入流模型［19］

風(fēng)電機(jī)組處于偏航入流狀態(tài)下時，尾流將向一側(cè)傾斜，導(dǎo)致誘導(dǎo)速度的大小和方向發(fā)生變化，示意圖如圖2所示。圖中[ψ]表示偏航角，角度偏向負(fù)[x]軸為正；[χ]表示尾跡偏斜角；[U0，l]表示局部平均風(fēng)速；[Ui]表示局部平均誘導(dǎo)速度。

偏航入流下的推力系數(shù)與誘導(dǎo)因子[a]的關(guān)系式為：

[CT=4a（1+a2-2acosψ）12]"""" （1）

式中：[ψ]——偏航角，rad。

同樣采用多項式的表達(dá)方式，用折減系數(shù)[Ra]和誘導(dǎo)因子[a]作為[CT]的函數(shù)，可導(dǎo)出：

[a=Ra（0.0883CT3+0.0586CT2+0.2460CT）] （2）

折減系數(shù)[Ra]與偏航角和推力系數(shù)有關(guān)，采用多項式的方式擬合，表示為：

[Ra=Ra，3min（CT，0.9）3+Ra，2min（CT，0.9）2+Ra，1min（CT，0.9）+1]"""""""""""""""""""" （3）

[Ra，i]是偏航角[ψ]的函數(shù)，表示為：

[Ra，i=Ai，3ψ3+Ai，2ψ2+Ai，1ψ]" （4）

其中系數(shù)矩陣[A]表示為：

[A=-0.51360.4438-0.16402.1735-2.61450.8646-2.07052.1667-0.6481]""" （5）

偏航導(dǎo)致的尾跡偏斜角[χ]可表示為：

[tanχ=U0，lsinψU0，lcosψ-Ui]""" （6）

由于尾跡在風(fēng)輪后偏斜，尾跡偏斜的一側(cè)誘導(dǎo)更高，且誘導(dǎo)隨葉片方位角變化，誘導(dǎo)速度表示為：

[Ui?=Ui（1+kxrsin?+kzrcos?）]""""" （7）

式中：[?]——方位角，rad；[r]——無量綱半徑；[kx]、[kz]——常數(shù)。根據(jù)文獻(xiàn)［19］，實現(xiàn)最佳相關(guān)性的取值為[kx=tan（0.4χ）]；[kz=0]。

2.2 結(jié)構(gòu)模型

單一的經(jīng)典線性Timoshenko或BernoullieEuler梁模型無法捕捉大型風(fēng)電機(jī)組葉片的非線性行為，HAWC2中風(fēng)電機(jī)組結(jié)構(gòu)模型基于多體動力學(xué)，將葉片離散成多個Timoshenko梁單元［20］。HAWC2中坐標(biāo)系示意圖如圖3所示。

柔性體上某一點位置可表示為：

[rp=R+Aup=R+Au0+uf]"""""" （8）

式中：[R]——慣性坐標(biāo)系原點到浮動坐標(biāo)系原點的矢量；[A]——浮動坐標(biāo)系變換到慣性坐標(biāo)系的旋轉(zhuǎn)變換矩陣；[up]——點在浮動坐標(biāo)系中的位置矢量，它是柔性體未變形之前的位置[u0]和[uf]的和，而[uf]可通過形函數(shù)[S]和自由度向量[qf]表示。

這點的速度可表示為：

[rp=R+Aup+ASqf=R+Aω×up+ASqf]"""""" （9）

將式（9）寫成如下矩陣的形式：

[rp=I-Aup×IASRωqf=Hq]" （10）

式中：[q]——廣義速度矢量。

這點的加速度可表示為：

[rp=Hq+Hq=I-Aup×IASRωqf+"""""" Aω×ω×up+2ω×up]"" （11）

對于各向異性梁單元，形函數(shù)[S]可表示為：

[S=TESNTNZTTES4]"nbsp; （12）

其中，

[TTES4=TTESTTESTTESTTES]" （13）

式中：[TES]——單元坐標(biāo)系到浮動坐標(biāo)系的標(biāo)準(zhǔn)正交轉(zhuǎn)換矩陣；[NT]——不在中心線上的點的旋轉(zhuǎn)矩陣；[NZ]——中心線的偏移和旋轉(zhuǎn)矩陣。

虛功和虛位移分別表示為：

[δW=Vδr?ρrdV=VHδq?ρrdV]"""""" （14）

[δr=δR+δAu0+ASδqf=I-Aup×IASδRδθδqf=Hδq]"""""""""""""""" （15）

將式（11）和（15）代入式（14）得到：

[δW=δqTMRωqf+Q]"""""" （16）

其中[M]為質(zhì)量矩陣，[Q]為載荷向量，表示為：

[M=VρI-Aup×IAS-Aup×Iup×ITup×I-up×ITSAS-up×ITSSTS dV] （17）

[Q=VρAω×ω×up+2ω×up-up×ITω×ω×up+2ω×upSTω×ω×up+2ω×updV]"" （18）

3 失控狀態(tài)下的穩(wěn)定性分析

仿真基礎(chǔ)條件設(shè)置如下：初始風(fēng)速為4 m/s，500～2100 s內(nèi)，風(fēng)速從4 m/s均勻增加到16 m/s，葉片槳距角始終設(shè)置為0°，發(fā)電機(jī)扭矩也為0，模擬風(fēng)電機(jī)組“失控”狀態(tài)下葉片的氣彈響應(yīng)?！笆Э亍睜顟B(tài)下，葉片上的氣流為附著流，風(fēng)輪轉(zhuǎn)速由風(fēng)速控制，是現(xiàn)實中可能發(fā)生顫振的情況［21］。

3.1 基 準(zhǔn)

圖4顯示了不考慮風(fēng)輪不平衡載荷時，風(fēng)輪轉(zhuǎn)速隨風(fēng)速的時間變化?？煽吹剑L(fēng)輪轉(zhuǎn)速在初始的500 s內(nèi)達(dá)到穩(wěn)定狀態(tài)，這個轉(zhuǎn)速稱為最終風(fēng)輪轉(zhuǎn)速［21］。500 s后，隨著風(fēng)速增大，更多風(fēng)能轉(zhuǎn)化為風(fēng)輪旋轉(zhuǎn)的動能，風(fēng)輪轉(zhuǎn)速逐漸增加。直到1428 s，轉(zhuǎn)速為13.058 r/min時，風(fēng)輪轉(zhuǎn)速出現(xiàn)跌落。

截取1350～1850 s內(nèi)葉尖的揮舞位移和攻角（attack of angle，AOA），如圖5所示?？煽吹?，風(fēng)輪轉(zhuǎn)速跌落前后，葉尖的揮舞位移和攻角都存在顯著差異，且轉(zhuǎn)速跌落后揮舞位移幅值和攻角3變化范圍持續(xù)增大，即出現(xiàn)發(fā)散現(xiàn)象。通過觀察1450～1451 s內(nèi)的揮舞位移和攻角，可發(fā)現(xiàn)兩者之間存在相位差，說明揮舞和扭轉(zhuǎn)運(yùn)動之間存在耦合。因此，在不考慮風(fēng)輪不平衡載荷的情況下，IEA 15 MW RWT的臨界顫振速度為13.058 r/min。

3.2 風(fēng)剪切

風(fēng)剪切對風(fēng)電機(jī)組的氣動載荷有重要影響，且隨著風(fēng)輪直徑和塔架高度的增加，風(fēng)剪切引起的整個風(fēng)輪旋轉(zhuǎn)平面的風(fēng)速空間變化更加明顯［22］。采用指數(shù)率描述風(fēng)剪切效應(yīng)，風(fēng)剪切系數(shù)分別取0.1、0.2、0.3、0.4。

圖6顯示了不同剪切系數(shù)下風(fēng)輪轉(zhuǎn)速隨時間的變化規(guī)律。可看出，隨著剪切系數(shù)的增大，臨界顫振速度升高。風(fēng)剪切系數(shù)增大至0.1、0.2、0.3、0.4時，臨界顫振速度分別為13.061、13.070、13.102、13.103 r/min，相較于基準(zhǔn)情況，分別提升了0.023%、0.092%、0.337%、0.345%。這說明風(fēng)剪切系數(shù)的增大會提高臨界顫振速度，可能與風(fēng)剪切系數(shù)的增大而減小了風(fēng)輪平面的平均推力有關(guān)。

3.3 偏航入流

由于風(fēng)向的隨機(jī)性以及風(fēng)電機(jī)組的偏航控制存在響應(yīng)時間，風(fēng)電機(jī)組經(jīng)常處于偏航狀態(tài)下運(yùn)行［23］?？紤]±20°以內(nèi)的偏航誤差，偏航角的正負(fù)根據(jù)圖2中的坐標(biāo)和方向來定義，不同偏航角下的風(fēng)輪轉(zhuǎn)速變化趨勢如圖7a、圖7b所示。

可看到，多個偏航角下的臨界顫振速度比較接近，變化較小。由圖7a可知，偏航角[ψ]為5°、10°、15°、20°時，臨界顫振速度為13.063、13.062、13.059、13.046 r/min。相較于基準(zhǔn)，[ψ]為5°、10°、15°時，臨界顫振速度分別提升了0.038%、0.038%、0.008%，而[ψ]為20°時，顫振速度降低了0.092%。從

圖6b可看出，偏航角為[-5°、-10°、-15°、-20°]時顫振速度為13.061、13.061、13.074、13.043 r/min，相較于基準(zhǔn)，[ψ]為[-5°、][-10°、][-15°]時，臨界顫振速度分別提升了0.023%、0.023%、0.123%，而[ψ]為[-20°]時，顫振速度降低了0.115%。這說明±15°以內(nèi)的偏航角能增大臨界顫振速度，而且在[-15°]時取得最大值，在[-20°]時取得最小值。

3.4 質(zhì)量不平衡

由于制造誤差或者葉片覆冰，風(fēng)電機(jī)組可能會在質(zhì)量不平衡的情況下運(yùn)行［7］。定義葉片質(zhì)量變化率[k=（m-m0）/m0]，[m0]表示葉片1初始質(zhì)量，那么[k]為正表示葉片1質(zhì)量偏重，[k]為負(fù)表示質(zhì)量偏輕，[k]的絕對值則表示不平衡程度的大小。

圖8a、圖8b分別顯示了葉片1質(zhì)量變化率為0（基準(zhǔn)）、0.01、0.02、0.03、0.04和[-0.01、-0.02、-0.03、-0.04]時的風(fēng)輪轉(zhuǎn)速。可看到由于葉片質(zhì)量不平衡，轉(zhuǎn)速跌落時間提前，臨界顫振轉(zhuǎn)速下降。同時，風(fēng)輪轉(zhuǎn)速波動明顯增加，且波動隨質(zhì)量不平衡程度增大而增大。[k]為0.01、0.02、0.03、0.04時，轉(zhuǎn)速極值分別為12.944、12.924、12.924、12.930 r/min。[k]為[-0.01、-0.02、-0.03、-0.04]時，轉(zhuǎn)速極值分別為12.988、12.994、13.020、13.045 r/min。

注意到，當(dāng)[kgt;0]時，隨著[k]的增大，轉(zhuǎn)速曲線整體趨勢更加偏下。但由于質(zhì)量不平衡程度更大，轉(zhuǎn)速波動更劇烈，導(dǎo)致[k]為0.02和0.03時轉(zhuǎn)速極值相等，甚至在[k]為0.04時轉(zhuǎn)速極值更大。因此為了合理地考慮轉(zhuǎn)速的整體變化趨勢以及更為保守地估計臨界顫振速度，對質(zhì)量不平衡狀態(tài)下的轉(zhuǎn)速進(jìn)行算數(shù)平均處理，轉(zhuǎn)速振動周期約為5 s，因此窗口大小設(shè)置為5 s，結(jié)果如圖8c、圖8d所示。可看到，平均處理后，曲線與原始數(shù)據(jù)整體趨勢一致，且曲線光滑，轉(zhuǎn)速波動較小，因此，可將處理后的局部轉(zhuǎn)速極值作為臨界顫振轉(zhuǎn)速。[k]為0.01、0.02、0.03、0.04時，臨界轉(zhuǎn)速為12.935、12.913、12.905、12.904 r/min。相較于基準(zhǔn)，臨界顫振速度降低了0.942%、1.110%、1.172%、1.179%。而[k]為[-0.01、-0.02、-0.03、][-0.04]時，臨界轉(zhuǎn)速為12.98、12.985、13.006、13.026 r/min。相較于基準(zhǔn)，臨界顫振速度降低了0.597%、0.559%、0.398%、0.245%。這說明質(zhì)量不平衡會降低臨界顫振速度，當(dāng)葉片質(zhì)量較輕，即[klt;0]時，顫振速度雖然相較于質(zhì)量平衡時有所下降，但隨著不平衡程度的增大，葉片質(zhì)量也減少，在假設(shè)剛度不變的前提下，葉片質(zhì)量減少能提高臨界顫振速度［24］。因此盡管轉(zhuǎn)速波動更大，但顫振速度卻不降反升。然而，當(dāng)葉片質(zhì)量較重，即[kgt;0]時，隨著不平衡程度的增大，葉片質(zhì)量增加，顫振速度下降。而且在同等不平衡程度下，相對于葉片質(zhì)量較輕（[klt;0]），葉片質(zhì)量較重（[kgt;0]）時的臨界顫振速度下降更為明顯。

3.5 組合工況

風(fēng)電機(jī)組運(yùn)行時可能會承受多種不平衡載荷，因此對前面所述的3種不平衡載荷進(jìn)行組合，研究風(fēng)電機(jī)組“失控”狀態(tài)下葉片的氣彈響應(yīng)，得到臨界顫振速度[Ωf]。由于風(fēng)剪切是客觀存在的，因此根據(jù)風(fēng)剪切系數(shù)的取值分為4組，共計324個組合工況，結(jié)果如圖9a～圖9d所示。

從圖9a中可看出，當(dāng)剪切系數(shù)為0.1時，在偏航角相同的情況下，臨界顫振速度[Ωf]隨質(zhì)量變化率的變化規(guī)律與只存在質(zhì)量不平衡時的規(guī)律一致，臨界顫振速度最大值[Ωf-max]在k=0時取得，最小值[Ωf-min]在k=0.04時取得。當(dāng)質(zhì)量變化率相同，臨界顫振速度[Ωf]隨偏航角的變化規(guī)律與質(zhì)量變化率的取值有關(guān)，不同程度的質(zhì)量不平衡狀態(tài)下的[Ωf-max]在0°～10°的偏航角下取得，[Ωf-min]在偏航角為-20°時取得，且偏航角為同等數(shù)值時，正偏航角下的臨界顫振速度更高。就整體而言，最大值[Ωf-max]為13.081 r/min，此時的偏航角為10°，葉片質(zhì)量變化率為0；最小值[Ωf-min]為12.882 r/min，此時的偏航角為[-20°]，葉片質(zhì)量變化率為0.04。相較于基準(zhǔn)，[Ωf-max]增大了0.191%，[Ωf-min]減小了1.333%。

圖9b顯示了剪切系數(shù)為0.2時，不同偏航角和質(zhì)量變化率下的臨界顫振速度。當(dāng)偏航角相同時，臨界顫振速度[Ωf]隨質(zhì)量變化率的變化規(guī)律與只存在質(zhì)量不平衡時的規(guī)律一致。當(dāng)質(zhì)量變化率相同時，[Ωf-max]在10°～20°的偏航角下取得，[Ωf-min]在偏航角為-20°時取得，且偏航角為同等數(shù)值時，正偏航角下的臨界顫振速度更高?？偟膩砜矗?dāng)偏航角為[10°]，葉片質(zhì)量變化率為0時，[Ωf]達(dá)到最大，為13.100 r/min；偏航角為[-20°]，葉片質(zhì)量變化率為0.04時，[Ωf]有最小值，為12.868 r/min。相較于基準(zhǔn)，[Ωf-max]增大了0.337%，[Ωf-min]減小了1.440%。

由圖9c可知，當(dāng)剪切系數(shù)為0.3時，相同偏航角下的臨界顫振速度[Ωf]隨質(zhì)量變化率的變化規(guī)律與只存在質(zhì)量不平衡時的規(guī)律一致。當(dāng)質(zhì)量變化率相同時，[Ωf-max]在10°～20°的偏航角下取得，[Ωf-min]在偏航角為[-20°]時取得，且偏航角為同等數(shù)值時，正偏航角下的臨界顫振速度更高。總的來看，當(dāng)偏航角為15°，葉片質(zhì)量變化率為0時，[Ωf]達(dá)到最大，為13.176 r/min；偏航角為[-20°]，葉片質(zhì)量變化率為0.04時，[Ωf]有最小值，為12.863 r/min。相較于基準(zhǔn)，[Ωf-max]增大了0.919%，[Ωf-min]減小了1.478%。

由圖9d所示，在剪切系數(shù)為0.4的情況下，當(dāng)偏航角相同時，臨界顫振速度[Ωf]隨質(zhì)量變化率的變化規(guī)律與只存在質(zhì)量不平衡時的規(guī)律一致。當(dāng)質(zhì)量變化率相同時，[Ωf-max]在10°～20°的偏航角下取得，[Ωf-min]在偏航角為[-20°]時取得，且偏航角為同等數(shù)值時，正偏航角下的臨界顫振速度更高?？偟膩砜矗琜Ωf]在偏航角為10°，葉片質(zhì)量變化率為0時，取得最大值13.169 r/min；在偏航角為[-20°]，葉片質(zhì)量變化率為0.04時，取得最小值12.865 r/min。相較于基準(zhǔn)，[Ωf-max]增大了0.866%，[Ωf-min]減小了1.463%。

4 結(jié) 論

基于葉片氣彈耦合分析模型，采用時域分析方法，以IEA 15 MW RWT為對象，研究3種風(fēng)輪不平衡載荷：風(fēng)剪切、偏航入流、質(zhì)量不平衡，單一作用或共同作用時，超大型風(fēng)電機(jī)組葉片的顫振特性。主要結(jié)論如下：

1）3種風(fēng)輪不平衡載荷單一作用時，風(fēng)剪切、偏航入流對于葉片臨界顫振速度影響較小，質(zhì)量不平衡對于葉片臨界顫振速度影響較大。臨界顫振速度隨著風(fēng)剪切系數(shù)的增大而增大，風(fēng)剪切系數(shù)為0.4時，臨界顫振速度[Ωf]增大了0.345%?！?5°以內(nèi)的偏航角能增大臨界顫振速度，且在[-15°]時取得最大值[Ωf-max]，[-20°]時取得最小值[Ωf-min]，相比于基準(zhǔn)情況，[Ωf-max]增大了0.123%，[Ωf-min]減小了0.115%。質(zhì)量不平衡會降低臨界顫振速度，且質(zhì)量偏重時，臨界顫振速度下降更明顯，當(dāng)質(zhì)量變化率[k=0.04]時，臨界顫振速度降低了1.179%。

2）3種風(fēng)輪不平衡載荷共同作用時，顫振速度最大值和最小值都出現(xiàn)在剪切系數(shù)為0.3時。任意剪切系數(shù)和偏航角下，臨界顫振速度[Ωf]隨質(zhì)量變化率的變化規(guī)律與只存在質(zhì)量不平衡時的規(guī)律一致，最大值[Ωf-max]在[k=0]時取得，最小值[Ωf-min]在[k=0.04]時取得。當(dāng)剪切系數(shù)和質(zhì)量變化率相同時，[Ωf]隨偏航角變化的規(guī)律與兩者的取值有關(guān)。但共同點是，取得最大值[Ωf-max]的偏航角都為非負(fù)值，而最小值[Ωf-min]都在偏航角為[-20°]時取得，且偏航角為同等數(shù)值時，正偏航角下的臨界顫振速度更高。共同作用時的[Ωf-max]為13.176 r/min，[Ωf-min]為12.863，相比于基準(zhǔn)情況，[Ωf-max]增大了0.919%，[Ωf-min]減小了1.478%。

3）在進(jìn)行超大型風(fēng)電機(jī)組氣彈穩(wěn)定性分析時，應(yīng)考慮多種風(fēng)輪不平衡載荷共同作用，特別是葉片質(zhì)量偏重以及偏航角度為-20°時這種最不利的情況，確保風(fēng)電機(jī)組在正常運(yùn)行期間不發(fā)生顫振失穩(wěn)。

［參考文獻(xiàn)］

［1］"""" GWEC. Global Offshore Wind Report 2021［R］. 2021.

［2］"""" VEERS P， DYKES K， LANTZ E， et al. Grand challenges in the science of wind energy［J］. Science，" 2019， 366（6464）： eaau2027.

［3］"""" CHANG L， YU Y J， LIU T R. Aeroelastic flutter and sliding mode control of wind turbine blade［J］. Shock and vibration， 2020， 2020： 1-8.

［4］"""" SHAKYA P， SUNNY M R， MAITI D K. Nonlinear flutter analysis of a bend-twist coupled composite wind turbine blade in time domain［J］. Composite structures， 2022， 284：115216.

［5］"""" HACH O， VERDONCK H， POLMAN J D， et al. Wind turbine stability： comparison of state-of-the-art aeroelastic simulation tools［J］. Journal of physics： conference series， 2020， 1618（5）： 052048.

［6］"""" HUA X G， MENG Q S， CHEN B， et al. Structural damping sensitivity affecting the flutter performance of a 10-MW offshore wind turbine［J］. Advances in structural engineering， 2020， 23（14）： 3037-3047.

［7］"""" HüBNER G R， PINHEIRO H， DE SOUZA C E， et al. Detection of mass imbalance in the rotor of wind turbines using" support" vector" machine［J］." Renewable" energy， 2021， 170：49-59.

［8］"""" REZAEIHA A， PEREIRA R， KOTSONIS M. Fluctuations of angle of attack and lift coefficient and the resultant fatigue loads for a large horizontal axis wind turbine［J］. Renewable energy， 2017， 114： 904-916.

［9］"""" TANG S Z， TIAN D， FANG J J， et al. Individual pitch controller characteristics analysis and optimization under aerodynamic imbalanced loads of wind turbines［J］. Energy reports， 2021， 7： 6489-6500.

［10］""" 田德， 方建駒， 劉楓， 等. 大型風(fēng)電機(jī)組模型預(yù)測獨立變槳控制器設(shè)計［J］. 太陽能學(xué)報， 2022， 43（4）： 461-467.

TIAN D， FANG J J， LIU F， et al. Design of model predictive individual pitch controller for large-scale wind turbine［J］. Acta energiae solaris sinica， 2022， 43（4）： 461-467.

［11］""" 馬磊明， 肖玲斐， 姜斌. 基于自適應(yīng)非奇異智能終端滑模觀測器的風(fēng)力機(jī)載荷增廣預(yù)測控制［J］. 太陽能學(xué)報， 2022， 43（11）： 259-268.

MA L M， XIAO L F， JIANG B. Augmented predictive intelligent control of wind turbine load based on adaptive nonsingular terminal sliding mode observer［J］. Acta energiae solaris sinica， 2022， 43（11）： 259-268.

［12］""" PIRRUNG G R， MADSEN H A， KIM T. The influence of trailed vorticity on flutter speed estimations［J］. Journal of physics： conference series， 2014， 524（1）： 012048.

［13］""" LARSEN T J， HANSEN A M. How 2 HAWC2， the user's manual［R］. Roskilde， Denmark：DTU Wind Energy， 2019.

［14］""" GAERTNER E， RINKER J， SETHURAMAN L， et al. Definition of the IEA wind 15-megawatt offshore reference wind turbine［R］. National Renewable Energy Laboratory， 2020.

［15］""" RINKER J， GAERTNER E， ZAHLE F， et al. Comparison of loads from HAWC2 and OpenFAST for the IEA wind 15 MW reference wind turbine［J］. Journal of physics： conference series， 2020， 1618（5）： 052052.

［16］""" 李貝， 田德， 唐世澤， 等. 超大型風(fēng)電機(jī)組葉片顫振分析及參數(shù)靈敏度研究［J］. 太陽能學(xué)報， 2023， 44（9）： 295-301.

LI B， TIAN D， TANG S Z， et al. Flutter analysis and parameter" sensitivity" study" of" ultra-large" wind" turbine blades［J］. Acta energiae solaris sinica， 2023， 44（9）： 295-301.

［17］""" EL KHCHINE Y， SRITI M. Tip loss factor effects on aerodynamic performances of horizontal axis wind turbine［J］. Energy procedia， 2017， 118：136-140.

［18］""" HANSEN M H. A Beddoes-Leishman type dynamic stall model in state-space and indicial formulations［R］. Riso National Laboratory， Denmark， 2004.

［19］""" MADSEN H A， LARSEN T J， PIRRUNG G R， et al. Implementation of the blade element momentum model on a polar grid and its aeroelastic load impact［J］. Wind energy science， 2020， 5（1）： 1-27.

［20］""" KIM T， HANSEN A M， BRANNER K. Development of an anisotropic beam finite element for composite wind turbine blades in multibody system［J］. Renewable energy， 2013， 59： 172-183.

［21］""" BAK C， ZAHLE F， BITSCHE R， et al. The DTU 10-MW reference wind turbine［R］. Roskilde， Denmark： DTU Wind Energy， 2013.

［22］""" LIU Y Q， QIAO Y H， HAN S， et al. Rotor equivalent wind speed calculation method based on equivalent power considering wind shear and tower shadow［J］. Renewable energy， 2021， 172： 882-896.

［23］""" DAI J C， YANG X， HU W， et al. Effect investigation of yaw on wind turbine performance based on SCADA data［J］. Energy， 2018， 149（APR.15）： 684-696.

［24］""" GRIFFITH D T， CHETAN M. Assessment of flutter prediction and trends in the design of large-scale wind turbine rotor blades［J］. Journal of physics： conference series， 2018， 1037： 042008.

STUDY ON FLUTTER LIMIT OF ULTRA-LARGE WIND TURBINE

BLADES UNDER ROTOR IMBALANCE LOADS

Tian De，Li Bei，Wu Xiaoxuan，Meng Huiwen，Wang Haodong，Su Yi

（State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources，

North China Electric Power University， Beijing 102206， China）

Abstract：This study focuses on the IEA Wind 15 MW Reference Wind Turbine and investigates the flutter limit of ultra-large wind turbine blades under imbalanced loads such as wind shear， yaw inflow， and mass imbalance. The results show that under a single unbalonced load， both wind shear and yaw misalignment within ±15° increases the critical flutter speed. Conversely， mass imbalance reduces the critical flutter speed， and the flutter speed is decreased more obviously under heavy mass. When all three imbalanced loads act simultaneously， the critical flutter speed exhibit maximum and minimum values at a shear coefficient of 0.3. The maximum value consistently occurs at a non-negative yaw angle， while the minimum value is observed at a yaw angle of -20°. The impact of mass imbalance on the critical flutter speed is consistent with that of single unbalanced load action. Therefore， it is crucial to consider the combined effects of multiple rotor imbalanced loads， particularly in the unfavorable scenario of heavier blade mass and a yaw angle is-20°， when conducting aeroelastic stability analysis.

Keywords：wind turbine blades； flutter； stability； aeroelasticity； rotor imbalance loads