摘 要：風(fēng)力機(jī)的大型化、深水化發(fā)展使得葉片尺寸急劇增加，長(zhǎng)柔葉片的氣彈效應(yīng)日趨顯著。該文以某95 m葉片為研究對(duì)象，基于BEM方法和多體動(dòng)力學(xué)理論構(gòu)建氣動(dòng)結(jié)構(gòu)耦合方程并求解極端風(fēng)載條件下的葉片氣彈響應(yīng)。系統(tǒng)地研究了發(fā)生顫振現(xiàn)象后葉片的振動(dòng)特征及不同轉(zhuǎn)速條件下的葉片阻尼特性，發(fā)現(xiàn)擺振負(fù)阻尼是誘導(dǎo)葉片發(fā)生顫振失穩(wěn)的關(guān)鍵因素；同時(shí)通過(guò)平行對(duì)照分析了空氣密度、變槳角度等因素對(duì)大型葉片氣彈穩(wěn)定性的影響。

關(guān)鍵詞：風(fēng)電葉片；氣動(dòng)彈性；阻尼；穩(wěn)定性分析；葉素動(dòng)量方法

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

DOI：10.19912/j.0254-0096.tynxb.2023-0658

文章編號(hào)：0254-0096（2024）08-0572-09

1. 廈門雙瑞風(fēng)電科技有限公司，廈門 361000；

2. 中國(guó)船舶集團(tuán)海裝風(fēng)電股份有限公司，重慶 401121；

3. 廣東海裝海上風(fēng)電研究中心有限公司，湛江 524000

0 引 言

近些年來(lái)，風(fēng)力機(jī)相關(guān)技術(shù)發(fā)展迅速，海上風(fēng)電迅速崛起，風(fēng)電的平價(jià)化、深水化發(fā)展已成為必然趨勢(shì)［1-2］。隨著海上浮式風(fēng)力機(jī)的大力發(fā)展，機(jī)組發(fā)電功率在不斷提高，為了滿足大型風(fēng)電設(shè)備對(duì)于風(fēng)能捕捉的需求，風(fēng)力機(jī)尺寸呈現(xiàn)出大型化發(fā)展趨勢(shì)，風(fēng)電葉片設(shè)計(jì)尺寸激增［3］。當(dāng)前葉片長(zhǎng)度已突破120 m量級(jí)，正朝向150 m以上的超大尺度快速前進(jìn)，為保證葉片總體重量不會(huì)急劇增加以及受限于成本制約，葉片普遍采用大柔性設(shè)計(jì)，這種超大柔性葉片在氣動(dòng)載荷作用下的彈性變形響應(yīng)愈發(fā)嚴(yán)重，風(fēng)力機(jī)的氣彈不穩(wěn)定問(wèn)題便尤為顯著，通常需要精細(xì)化的分析和安全性評(píng)估。

當(dāng)大柔性葉片在近地面大氣環(huán)境下運(yùn)行時(shí)，復(fù)雜而隨機(jī)的氣動(dòng)力、慣性力和彈性力等相互耦合使大柔性葉片的非線性氣彈效應(yīng)極為突出，這其中主要有失速誘導(dǎo)的振動(dòng)和顫振問(wèn)題［4-5］，第一種又被稱為失速顫振?，F(xiàn)有風(fēng)力機(jī)設(shè)計(jì)標(biāo)準(zhǔn)，要求風(fēng)力機(jī)在停機(jī)狀態(tài)下，能抵抗50年一遇的大風(fēng)且風(fēng)向并不固定，其中最極端的情況就是某風(fēng)向產(chǎn)生了失速誘發(fā)的振動(dòng)問(wèn)題，所以需要考慮葉片部分截面來(lái)流的流動(dòng)分離問(wèn)題，即失速狀態(tài)。同時(shí)風(fēng)力機(jī)在工作狀態(tài)葉片持續(xù)旋轉(zhuǎn)，此時(shí)在氣動(dòng)力、重力和離心力作用下，葉片不斷振動(dòng)誘發(fā)動(dòng)態(tài)失速，從而使氣動(dòng)力脈動(dòng)加劇，進(jìn)而結(jié)構(gòu)所受到的氣動(dòng)載荷變得復(fù)雜，造成結(jié)構(gòu)振動(dòng)加劇。動(dòng)態(tài)載荷與大柔性葉片耦合會(huì)帶來(lái)劇烈的隨機(jī)和自激振動(dòng)，被認(rèn)為是當(dāng)前葉片折斷和整機(jī)倒塌的主因，給機(jī)組安全運(yùn)行帶來(lái)嚴(yán)重挑戰(zhàn)。國(guó)內(nèi)外學(xué)者做了大量有關(guān)試驗(yàn)和理論方面的研究，以研究風(fēng)電葉片的氣彈穩(wěn)定性。在數(shù)值模擬方面，Yu等［6］基于CFD-CSD耦合方法，發(fā)現(xiàn)葉片氣動(dòng)扭轉(zhuǎn)變形對(duì)非定常氣動(dòng)載荷影響顯著，但CFD-CSD仿真采用松耦合精度較低，存在明顯的滯后效應(yīng)；陳川等［7］通過(guò)有限單元法求解不同方位角和不同槳距角下的葉片氣動(dòng)阻尼，總結(jié)了方位角對(duì)葉片氣動(dòng)阻尼的影響，進(jìn)而歸納出獨(dú)立變槳對(duì)大型葉片氣彈穩(wěn)定性的提升作用；黃俊東［8］應(yīng)用多體動(dòng)力學(xué)方法建模采用“超級(jí)單元”構(gòu)建了葉片等效模型考慮了剛?cè)狁詈霞胺嵌ǔ鈴楍詈闲?yīng)，發(fā)現(xiàn)顫振失穩(wěn)時(shí)振型由單向振動(dòng)演變?yōu)槎喾较蛘駝?dòng)，風(fēng)力機(jī)氣彈失穩(wěn)主要表現(xiàn)在擺振方向上。風(fēng)洞試驗(yàn)是揭示葉片強(qiáng)非線性運(yùn)動(dòng)的最直觀有效手段之一，但由于風(fēng)電葉片翼型復(fù)雜不規(guī)則，葉片預(yù)彎、后掠、扭轉(zhuǎn)角的存在使得葉片剛度質(zhì)量分布縮比設(shè)計(jì)難度大，測(cè)點(diǎn)布置難、采集干擾性強(qiáng)、洞壁效應(yīng)修正精度不足等試驗(yàn)困難導(dǎo)致國(guó)內(nèi)外較少開(kāi)展超長(zhǎng)柔性葉片三維顫振彈性模型試驗(yàn)研究?？率捞玫龋?］基于變截面梁設(shè)計(jì)某15 MW風(fēng)電葉片三維彈性縮比模型，進(jìn)行了同步測(cè)振、測(cè)力、氣彈風(fēng)洞試驗(yàn)，研究表明：風(fēng)電葉片的槳距角存在風(fēng)振敏感區(qū)間，在該區(qū)間內(nèi)超過(guò)臨界風(fēng)速即可發(fā)生大幅鎖頻振動(dòng)，同時(shí)存在能量積累突變界線，超過(guò)該界線對(duì)應(yīng)風(fēng)速后的能量積累尤為顯著，表現(xiàn)為風(fēng)致振動(dòng)能量隨時(shí)間呈現(xiàn)顯著的非平穩(wěn)特性，氣動(dòng)負(fù)阻尼是結(jié)構(gòu)系統(tǒng)發(fā)散從而誘發(fā)顫振的主要原因。白學(xué)宗等［10］基于試驗(yàn)方法驗(yàn)證了突變風(fēng)速載荷沖擊會(huì)帶來(lái)葉片疲勞損傷的大幅增，總結(jié)了突變載荷下葉片剛度的退化規(guī)律。楊景云等［11］分別構(gòu)建低速和高速旋轉(zhuǎn)狀態(tài)下的葉片結(jié)構(gòu)方程，發(fā)現(xiàn)轉(zhuǎn)速對(duì)葉片剛度的具有明顯的影響，其中低速轉(zhuǎn)動(dòng)會(huì)使葉片剛度降低而高速轉(zhuǎn)動(dòng)會(huì)使葉片剛度明顯提高，且對(duì)于大型葉片影響更為明顯。針對(duì)長(zhǎng)柔葉片的氣彈耦合問(wèn)題，本文采用葉素動(dòng)量方法結(jié)合模態(tài)分析法分別對(duì)葉片的氣動(dòng)、結(jié)構(gòu)完成建模，研究葉片極限工況下的氣彈響應(yīng)，對(duì)葉片氣彈失穩(wěn)現(xiàn)象及其影響因素進(jìn)行歸納，為長(zhǎng)柔葉片的設(shè)計(jì)和風(fēng)電場(chǎng)運(yùn)維提供參考，具有工程價(jià)值。

1 理論方法

1.1 葉素動(dòng)量理論

葉素動(dòng)量理論（bladed element momentum， BEM）方法是一維動(dòng)量理論和二維葉素理論的結(jié)合。1920年，由Betz［12］首次將動(dòng)量理論應(yīng)用于風(fēng)力機(jī)領(lǐng)域，認(rèn)為風(fēng)輪吸收的能量是入流通過(guò)風(fēng)輪后減少的能量。而后Glauert［13］于1935年在此基礎(chǔ)上加入了葉素理論，從而形成了一個(gè)相對(duì)完整的葉素動(dòng)量理論體系。由于BEM方法采用了一系列簡(jiǎn)化的假設(shè)，因此可能會(huì)導(dǎo)致計(jì)算結(jié)果出現(xiàn)一些誤差。為了提高BEM方法的計(jì)算精度，使葉素動(dòng)量方法計(jì)算結(jié)果更加可靠，研究人員提出多種基于半經(jīng)驗(yàn)公式的修正模型。Prandtl［14］于1963年提出一種葉尖損失修正模型，引入葉尖損失系數(shù)，該系數(shù)由半經(jīng)驗(yàn)公式導(dǎo)出，以修正BEM理論。1981年，Pitt等［15］基于旋翼飛機(jī)的斜尾流模型，改進(jìn)了適用于水平軸風(fēng)力機(jī)的斜尾修正模型，使BEM模型更準(zhǔn)確地計(jì)算風(fēng)機(jī)在偏航狀態(tài)下的氣動(dòng)性能。1986年，Leishman等［16］提出動(dòng)態(tài)失速效應(yīng)，并將其引入風(fēng)電葉片的氣動(dòng)計(jì)算中。采用半經(jīng)驗(yàn)公式對(duì)升力和阻力系數(shù)進(jìn)行修正，有效地提高了非穩(wěn)態(tài)入流下BEM計(jì)算結(jié)果的準(zhǔn)確性。

葉素動(dòng)量理論的基本思想是將風(fēng)電葉片沿展向離散為若干段長(zhǎng)度為[r]的葉素。葉輪旋轉(zhuǎn)使得每一段葉素都會(huì)掃掠出一個(gè)圓環(huán)，并假設(shè)每一段葉素之間不存在相互影響，通過(guò)臨近葉素掃掠圓環(huán)的入流不會(huì)對(duì)其產(chǎn)生誘導(dǎo)作用。葉素理論將每一段葉素視為二維翼型，葉片整體的受力情況可由每一段葉素上二維翼型的受力沿展向積分來(lái)近似表示，葉素受力分析如圖1所示。

圖1中，[D]為阻力；[L]為升力； [V01-a]為軸向速度分量；[ωr1+a′]為速度切向分量；[V0]為遠(yuǎn)方入流風(fēng)速；[ω]為風(fēng)輪旋轉(zhuǎn)速度；[a]為軸向誘導(dǎo)因子，[a′]為切向誘導(dǎo)因子；[α]為攻角；[β]為扭轉(zhuǎn)角；[φ]為入流角。

軸向誘導(dǎo)因子[a]和切向誘導(dǎo)因子[a′]分別定義為：

[a=ΔVV0=V0-VV0] （1）

[a′=Vt2ωr] （2）

式中：[V]——風(fēng)輪前風(fēng)速；[Vt]——切向速度分量。

軸向誘導(dǎo)速度分量為：

[vd=v∞1-a] （3）

切向誘導(dǎo)速度分量為：

[vt=1+a′ωr] （4）

合成相對(duì)入流速度為：

[w=V2d+V2t=v2∞1-a2+ωr21+a′2] （5）

入流角為：

[φ=tan-1V∞1-a1+a′ωr] （6）

根據(jù)動(dòng)量守恒定理，軸向推力與風(fēng)輪轉(zhuǎn)矩可表達(dá)為：

[dF=4πρv2∞a1-ardr] （7）

[dM=4πρωv∞1-aa′r3dr] （8）

再根據(jù)葉素理論方程，推力與風(fēng)輪轉(zhuǎn)矩可分別表達(dá)為：

[dF=NLcosφ+Dsinφ=12ρw2NcCLcosφ+CDsinφdr] （9）

[dM=NLsinφ-Dcosφr=12ρw2NcCLsinφ-CDcosφdr] （10）

式（7）和式（9）分別表示從動(dòng)量守恒和葉素理論導(dǎo)出的風(fēng)輪推力平衡方程。式（8）和式（10）分別表示從動(dòng)量守恒和葉素理論導(dǎo)出的風(fēng)輪扭矩方程。聯(lián)立式（7）～式（10），可獲得軸向誘導(dǎo)因子和切向誘導(dǎo)因子與葉片弦長(zhǎng)和入流角之間的關(guān)系為：

[a1-a=Nc8πr·CLcosφ+CDsinφ（sinφ）2] （11）

[a′1+a′=Nc8πr·CLsinφ-CDcosφsinφcosφ] （12）

式中：[CL]——翼型的升力系數(shù)；[CD]——翼型的阻力系數(shù)，其數(shù)值可由翼型參數(shù)數(shù)據(jù)表獲得。

在計(jì)算時(shí)，首先初始化[a]和[a′]，假設(shè)[a=a′=0]；然后根據(jù)式（6）計(jì)算入流角[φ]。在得到入流角后，翼型對(duì)應(yīng)的攻角即可獲得，參照翼型數(shù)據(jù)表可得到對(duì)應(yīng)攻角下的升阻力系數(shù)。

[α=φ-β] （13）

將升阻力系數(shù)和入流角帶入式（11）和式（12）即可求得軸向誘導(dǎo)因子[a]和切向誘導(dǎo)因子[a′]，至此完成一次迭代過(guò)程。重復(fù)迭代此過(guò)程，直至[a]和[a′]均收斂。再將收斂的[a]和[a′]值回代至式（9）和式（10），即可求得葉片上的氣動(dòng)載荷。葉素動(dòng)量方法被認(rèn)為是計(jì)算氣動(dòng)載荷問(wèn)題的有效方法，在工程上得到廣泛應(yīng)用。

1.2 葉片結(jié)構(gòu)動(dòng)力學(xué)方程

在風(fēng)力機(jī)發(fā)展初期葉片塔架等結(jié)構(gòu)往往采用剛體建模，忽略了彈性變形所帶來(lái)的影響，20世紀(jì)70年代末開(kāi)始有研究人員使用彈性梁模型來(lái)建立葉片模型［17］，結(jié)合了更為可靠的動(dòng)態(tài)分析方法。目前針對(duì)大型風(fēng)電柔性葉片的結(jié)構(gòu)動(dòng)力學(xué)分析方法主要有兩種：模態(tài)分析法和有限元法。其中模態(tài)分析方法通過(guò)前處理手段將二維截面特性施加在一維梁?jiǎn)卧Ｐ蜕弦源藖?lái)達(dá)到等效三維葉片建模目的，再配合多體動(dòng)力學(xué)方程建?？蓸O大程度上減少計(jì)算模型的自由度，對(duì)計(jì)算效率提升明顯，此外模態(tài)分析方法中可使用模態(tài)阻尼直觀地描述柔性葉片的結(jié)構(gòu)阻尼，方便地計(jì)算葉片單元的動(dòng)能和應(yīng)變，這對(duì)于大型柔性葉片的氣彈分析有重要意義。

如圖2所示，在葉片梁模型中一個(gè)單元共2個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)包含6個(gè)自由度，即節(jié)點(diǎn)位移[x、y、z]，節(jié)點(diǎn)扭轉(zhuǎn)角[θx、θy、θz]，梁軸線和截面的位移模式遵循小變形假設(shè)。單元坐標(biāo)系以內(nèi)側(cè)第一個(gè)節(jié)點(diǎn)為原點(diǎn)。單元坐標(biāo)系以第三軸z軸從內(nèi)側(cè)節(jié)點(diǎn)指向外側(cè)節(jié)點(diǎn)，其余兩軸線[x]軸和[y]軸與截面彎曲第一主軸和第二主軸平行，且正方向分別指向前緣和吸力面。

葉片結(jié)構(gòu)動(dòng)力學(xué)模型可通過(guò)式（14）動(dòng)力學(xué)方程描述：

[McrrMcrεDcTrMcTrεMcεεDcTεDcrDcε0vcεcλc+0CCεε0εc+0KCεε0εc=fca+fc0σca+σc0ac2-fciσci0] （14）

式中：[vc]——節(jié)點(diǎn)的速度矢量；[εc]——廣義應(yīng)變矢量；[λc]——拉格朗日乘子矢量；[fci]和[σci]——非線性慣性力矢量和節(jié)點(diǎn)速度、應(yīng)變率對(duì)應(yīng)的應(yīng)力矢量；[fca]——施加的節(jié)點(diǎn)力矢量；[σca]——由于廣義應(yīng)變率施加的廣義應(yīng)力矢量；[fc0]——由于節(jié)點(diǎn)速度帶來(lái)的鉸鏈反力矢量；[σc0]——對(duì)應(yīng)于制定應(yīng)變的廣義限制應(yīng)力矢量；[ac2]——加速度轉(zhuǎn)換的傳導(dǎo)項(xiàng)；[Mcrr]，[Mcrε]和[Mcεε]——節(jié)點(diǎn)位移和應(yīng)變率對(duì)應(yīng)的質(zhì)量矩陣子塊；[CCεε]——對(duì)應(yīng)于應(yīng)變率的結(jié)構(gòu)阻尼矩陣；[KCεε]——對(duì)應(yīng)于應(yīng)變率的結(jié)構(gòu)剛度矩陣；[Dcr]和[Dcε]——與節(jié)點(diǎn)位移和應(yīng)變率有關(guān)的變量。將所有梁?jiǎn)卧馁|(zhì)量、阻尼和剛度矩陣以及右側(cè)項(xiàng)進(jìn)行組裝，即可得到整個(gè)柔性葉片的結(jié)構(gòu)運(yùn)動(dòng)方程。

[Mq+Cq+Kq=p] （15）

式中：[M]——葉片的質(zhì)量矩陣；[C]——葉片的阻尼矩陣；[K]——葉片的剛度矩陣；[q]、[q]、[q]——葉片單元的加速度項(xiàng)、速度項(xiàng)以及位移項(xiàng)；[p]——外力載荷即氣動(dòng)載荷及慣性載荷。直接求解運(yùn)動(dòng)方程，可得到梁的瞬態(tài)時(shí)歷響應(yīng)；求解相應(yīng)的特征值問(wèn)題，可得到振動(dòng)頻率（特征值）和振型（特征向量）；用特征向量左乘運(yùn)動(dòng)方程，利用結(jié)構(gòu)模態(tài)的正交性質(zhì)，可將系統(tǒng)降階為一個(gè)簡(jiǎn)化的模態(tài)系統(tǒng)。進(jìn)而可完成瞬態(tài)氣彈響應(yīng)仿真等內(nèi)容。如果無(wú)外加氣動(dòng)力，則可進(jìn)行結(jié)構(gòu)模態(tài)分析，從多個(gè)角度分析風(fēng)電葉片的動(dòng)力學(xué)穩(wěn)定特性。

2 葉片計(jì)算模型

本文選取計(jì)算模型為某95 m葉片模型，葉片長(zhǎng)度為95 m，匹配某6.x MW機(jī)組，葉輪模型見(jiàn)圖3。設(shè)計(jì)額定轉(zhuǎn)速9 r/min，額定風(fēng)速為9.5 m/s，匹配風(fēng)區(qū)等級(jí)為IEC ⅢC類風(fēng)區(qū)。具體參數(shù)如表1所示。

2.1 葉片模態(tài)分析

通過(guò)試驗(yàn)和仿真分析得到葉片的固有頻率，其中現(xiàn)場(chǎng)實(shí)測(cè)時(shí)在葉尖激勵(lì)葉片擺動(dòng)，保持葉尖擺動(dòng)幅值不超過(guò)傳感器采集信號(hào)的最大量程，隨后釋放使其發(fā)生自由振動(dòng)，在葉片不同截面處布置位移傳感器和加速度傳感器測(cè)得其位移信號(hào)，進(jìn)而求得葉片固有頻率，傳感器安裝位置及現(xiàn)場(chǎng)實(shí)測(cè)見(jiàn)圖4。

表2中給出了該葉片的仿真模態(tài)頻率及現(xiàn)場(chǎng)實(shí)測(cè)頻率，仿真模型的前4階頻率以及扭轉(zhuǎn)頻率與葉片實(shí)測(cè)頻率吻合良好。其中Bladed模型計(jì)算頻率與樣片測(cè)試頻率最大誤差為5.48%對(duì)應(yīng)一階擺振方向，出現(xiàn)誤差的主要原因在于測(cè)試樣片存在一定的測(cè)試補(bǔ)強(qiáng)，以避免在進(jìn)行靜力測(cè)試和疲勞測(cè)試過(guò)程中出現(xiàn)加載點(diǎn)損傷，因此樣片測(cè)試頻率略低。圖5中分別為對(duì)應(yīng)的各階頻率模態(tài)振型。

2.2 穩(wěn)態(tài)計(jì)算結(jié)果

通過(guò)BEM穩(wěn)態(tài)計(jì)算得到葉輪的穩(wěn)態(tài)運(yùn)轉(zhuǎn)功率、轉(zhuǎn)速隨風(fēng)速變化的曲線。圖6給出了機(jī)組在切入、切出風(fēng)速區(qū)間內(nèi)的穩(wěn)態(tài)運(yùn)行功率、轉(zhuǎn)速隨風(fēng)速的變化關(guān)系。葉片運(yùn)轉(zhuǎn)額定轉(zhuǎn)速為9 r/min，額定風(fēng)速為9.5 m/s。

3 極限工況下的風(fēng)力機(jī)時(shí)域響應(yīng)分析

為研究出現(xiàn)失穩(wěn)時(shí)風(fēng)力機(jī)的響應(yīng)特征，本文選取極端工況計(jì)算極限工況下的風(fēng)力機(jī)動(dòng)態(tài)響應(yīng)。時(shí)域數(shù)值仿真部分分為飛車工況及停機(jī)偏航工況下的響應(yīng)計(jì)算。

3.1 飛車工況響應(yīng)分析

此部分模擬風(fēng)輪在不接發(fā)電機(jī)的情況下發(fā)生飛車時(shí)的狀態(tài)。在特定的風(fēng)速和槳距角下，風(fēng)輪的轉(zhuǎn)速最終會(huì)達(dá)到一個(gè)臨界值，風(fēng)輪在此轉(zhuǎn)速下所受到的氣動(dòng)力矩為零，維持勻轉(zhuǎn)速運(yùn)動(dòng)，隨著風(fēng)速的不斷增大風(fēng)輪轉(zhuǎn)速不斷提升直至發(fā)生顫振的狀態(tài)。

在穩(wěn)態(tài)遞增的風(fēng)速條件下風(fēng)輪運(yùn)轉(zhuǎn)情況如圖7所示，模擬過(guò)程中保持葉片槳距角為0°，結(jié)果顯示在200 s時(shí)葉尖揮舞方向扭轉(zhuǎn)方向振動(dòng)幅值出現(xiàn)突增，隨后擺振方向的振動(dòng)也開(kāi)始逐步增大開(kāi)始發(fā)散。風(fēng)輪在14.8 r/min轉(zhuǎn)速附近發(fā)生顫振，此時(shí)葉尖的臨界顫振速度在約147 m/s，遠(yuǎn)超出機(jī)組的正常運(yùn)行范圍。

分別計(jì)算了空氣密度為1.000、1.100、1.225和1.3 kg/m3時(shí)的風(fēng)輪顫振特性，表3統(tǒng)計(jì)了這幾種情況下的臨界顫振轉(zhuǎn)速和葉尖速度。在相同的穩(wěn)態(tài)遞增風(fēng)速工況條件下，隨著空氣密度的增大風(fēng)輪出現(xiàn)失穩(wěn)的時(shí)間節(jié)點(diǎn)提前，葉片的臨界顫振轉(zhuǎn)速隨著空氣密度的增加會(huì)略有降低，說(shuō)明在空氣密度較大的風(fēng)場(chǎng)運(yùn)行時(shí)葉片的顫振風(fēng)險(xiǎn)有增大的趨勢(shì)。各組空氣密度條件下，葉片動(dòng)態(tài)表現(xiàn)如圖8所示。

3.2 停機(jī)狀態(tài)響應(yīng)分析

不同于飛車工況，在出現(xiàn)極端風(fēng)速陣風(fēng)時(shí)風(fēng)力機(jī)往往會(huì)采取停機(jī)順槳策略，以規(guī)避極端載荷，避免發(fā)生危險(xiǎn)。

本節(jié)評(píng)估了風(fēng)輪在停機(jī)狀態(tài)下遭遇極端風(fēng)速的響應(yīng)特性。計(jì)算了停機(jī)狀態(tài)風(fēng)輪偏航角分別為0°、30°和60°這3種情況的動(dòng)態(tài)響應(yīng)。模擬過(guò)程中葉片始終保持90°順槳狀態(tài)，輪轂中心風(fēng)速?gòu)? m/s 增加到120 m/s。圖9為風(fēng)輪動(dòng)態(tài)響應(yīng)結(jié)果。

結(jié)果顯示，在停機(jī)順槳狀態(tài)下，風(fēng)輪偏航角為0°時(shí)風(fēng)輪在102 m/s以上的風(fēng)速時(shí)會(huì)出現(xiàn)不穩(wěn)定狀態(tài)；風(fēng)輪偏航角為30°時(shí)風(fēng)輪在50 m/s風(fēng)速附近會(huì)出現(xiàn)不穩(wěn)定狀態(tài)；風(fēng)輪偏航角為60°時(shí)風(fēng)輪在84 m/s以上的風(fēng)速時(shí)會(huì)出現(xiàn)不穩(wěn)定狀態(tài)。相比于正面對(duì)風(fēng)工況，存在入流偏航時(shí)會(huì)大大降低風(fēng)力機(jī)的安全裕量，因此在停機(jī)順槳狀態(tài)下需要避免長(zhǎng)時(shí)間的偏航入流，機(jī)組需要主動(dòng)對(duì)風(fēng)以規(guī)避偏航狀態(tài)下的不穩(wěn)定狀態(tài)。0°偏航時(shí)風(fēng)輪在102 m/s以上的風(fēng)速時(shí)才會(huì)出現(xiàn)不穩(wěn)定狀態(tài)，遠(yuǎn)高于Ⅲ類風(fēng)場(chǎng)52.5 m/s的極端風(fēng)速條件。在保證偏航對(duì)風(fēng)時(shí)停機(jī)順槳狀態(tài)下機(jī)組不會(huì)出現(xiàn)顫振失穩(wěn)風(fēng)險(xiǎn)。

4 葉片頻域響應(yīng)分析

本節(jié)首先重點(diǎn)分析在葉片0°槳距角下風(fēng)輪各階模態(tài)頻率和阻尼隨轉(zhuǎn)速的變化情況，隨后評(píng)估槳距角變化對(duì)風(fēng)輪系統(tǒng)穩(wěn)定性的影響。

4.1 風(fēng)輪穩(wěn)定性分析

計(jì)算分析了0°槳距角下風(fēng)輪系統(tǒng)頻率小于2 Hz的模態(tài)特征，圖10中給出了風(fēng)輪前11階模態(tài)頻率隨轉(zhuǎn)速的變化曲線。其中一階揮舞頻率隨著轉(zhuǎn)速的增加有所增大，這是由于氣彈剛化的效應(yīng)，此時(shí)一階揮舞模態(tài)和扭轉(zhuǎn)相耦合，葉片加速時(shí)截面會(huì)產(chǎn)生抬頭動(dòng)作，導(dǎo)致截面氣動(dòng)攻角增大，將在加速度方向上產(chǎn)生更大的升力，而此時(shí)氣動(dòng)力和葉片變形的相位相反，使得葉片剛度增加。

圖11為風(fēng)輪前11階模態(tài)阻尼隨轉(zhuǎn)速的變化曲線，在機(jī)組額定轉(zhuǎn)速內(nèi)未發(fā)現(xiàn)風(fēng)輪各階模態(tài)有負(fù)阻尼產(chǎn)生，說(shuō)明在正常運(yùn)轉(zhuǎn)條件下機(jī)組穩(wěn)定不存在失穩(wěn)狀態(tài)。

觀察所有計(jì)算模態(tài)阻尼，其中有2階模態(tài)在9.8 r/min以后出現(xiàn)負(fù)阻尼現(xiàn)象，有二階模態(tài)在12.1 r/min時(shí)出現(xiàn)負(fù)阻尼情況，有一階模態(tài)在13 r/min附近出現(xiàn)負(fù)阻尼情況，說(shuō)明該風(fēng)輪系統(tǒng)在轉(zhuǎn)速高于9.8 r/min后將開(kāi)始可能產(chǎn)生不穩(wěn)定性風(fēng)險(xiǎn)。各階不穩(wěn)定模態(tài)出現(xiàn)負(fù)阻尼時(shí)的風(fēng)輪轉(zhuǎn)速和對(duì)應(yīng)頻率見(jiàn)表4，風(fēng)輪的不穩(wěn)定頻率在1.64 Hz以內(nèi)。

產(chǎn)生負(fù)阻尼的模態(tài)集中在擺振方向，表4中風(fēng)輪系統(tǒng)的各階模態(tài)又是風(fēng)輪的結(jié)構(gòu)模態(tài)和氣動(dòng)狀態(tài)的組合。利用Bladed坎貝爾圖分析工具將風(fēng)輪系統(tǒng)的各階模態(tài)主要分量分離。例如在轉(zhuǎn)速9.82 r/min下，“一階擺振滯后模態(tài)”的模態(tài)出現(xiàn)負(fù)阻尼，頻率為0.362 Hz，風(fēng)輪系統(tǒng)的該階模態(tài)主要分量構(gòu)成在表5中給出，是風(fēng)輪的結(jié)構(gòu)模態(tài)和氣動(dòng)模態(tài)的耦合結(jié)果，每一分量的貢獻(xiàn)比例不同。

其中貢獻(xiàn)最大的是轉(zhuǎn)子一階擺振正弦循環(huán)分量和一階擺振余弦循環(huán)分量模態(tài)，這兩種模態(tài)振型為葉片擺振方向主導(dǎo)陣型如圖12所示，因而葉片擺振方向的阻尼對(duì)風(fēng)輪的穩(wěn)定性具有主要影響。

頻域分析結(jié)果顯示大柔性葉片的不穩(wěn)定性主要來(lái)源于擺振方向的振動(dòng)及其對(duì)應(yīng)的負(fù)阻尼，因此可嘗試通過(guò)提高擺振方向的結(jié)構(gòu)阻尼用以提高葉片整體穩(wěn)定性。

4.2 槳距角影響

通常在風(fēng)速超過(guò)額定風(fēng)速后葉片會(huì)進(jìn)行變槳，以保持相對(duì)穩(wěn)定的輸出功率同時(shí)減小葉片沖擊載荷。補(bǔ)充分析了葉片在4°槳距角和8°槳距角狀態(tài)下風(fēng)輪的阻尼表現(xiàn)，計(jì)算結(jié)果如圖13所示。

4°槳距角時(shí)，一階擺振模態(tài)出現(xiàn)負(fù)阻尼對(duì)應(yīng)轉(zhuǎn)速提高約10%，二階擺振模態(tài)在仿真過(guò)程中未表現(xiàn)出負(fù)阻尼不穩(wěn)定狀態(tài)，風(fēng)輪整體穩(wěn)定性有明顯提高；在8°槳距角下各階模態(tài)已經(jīng)完全不存在負(fù)阻尼情況，且隨著轉(zhuǎn)速的增大各階擺振阻尼有提高現(xiàn)象?？傮w來(lái)講隨著葉片槳距角的增加，風(fēng)輪穩(wěn)定性得到提升，特別是轉(zhuǎn)速較大時(shí)風(fēng)輪的阻尼迅速提高。因此可考慮適當(dāng)犧牲一定發(fā)電量在額定風(fēng)速之前提前變槳，以提高風(fēng)力機(jī)系統(tǒng)穩(wěn)定性同時(shí)可一定程度上降低載荷。

5 結(jié) 論

本文中使用Bladed軟件對(duì)大型柔性風(fēng)電葉片建模仿真計(jì)算，結(jié)合時(shí)域仿真和頻域分析方法系統(tǒng)分析大型機(jī)組葉片的氣彈動(dòng)態(tài)響應(yīng)過(guò)程及其顫振失穩(wěn)影響因素，主要結(jié)論如下：

1）在正常運(yùn)轉(zhuǎn)工況轉(zhuǎn)速范圍內(nèi)，該95 m葉片各階模態(tài)未出現(xiàn)負(fù)阻尼情況。同時(shí)時(shí)域仿真結(jié)果表明0°槳距角下機(jī)組在約14.8 r/min出現(xiàn)失穩(wěn)，此時(shí)葉尖線速度接近147 m/s，與正常運(yùn)行范圍之間具有約70%裕量。因此葉片在正常運(yùn)行的工況范圍內(nèi)，其動(dòng)力學(xué)特性是穩(wěn)定的，無(wú)氣彈失穩(wěn)風(fēng)險(xiǎn)。

2）葉片在轉(zhuǎn)速超過(guò)9.82 r/min后，擺振方向部分模態(tài)出現(xiàn)負(fù)阻尼，風(fēng)輪開(kāi)始產(chǎn)生不穩(wěn)定風(fēng)險(xiǎn)。葉片出現(xiàn)氣彈失穩(wěn)主要源于擺振方向的振動(dòng)及其較低的阻尼特性。

3）葉片槳距角增加使得風(fēng)輪整體穩(wěn)定性明顯提升，特別是在大轉(zhuǎn)速下風(fēng)輪的各階模態(tài)阻尼迅速提高，可選擇適當(dāng)提前變槳來(lái)提升風(fēng)力機(jī)穩(wěn)定性。

4）研究了空氣密度影響，目標(biāo)風(fēng)場(chǎng)的空氣密度降低，對(duì)風(fēng)力機(jī)系統(tǒng)的穩(wěn)定性提升有正面作用。

5）停機(jī)順槳工況下，正對(duì)風(fēng)向時(shí)在風(fēng)速達(dá)到102 m/s以上才會(huì)產(chǎn)生不穩(wěn)定性顫振現(xiàn)象，在存在偏航角度時(shí)，安全風(fēng)速明顯降低。因此應(yīng)避免長(zhǎng)時(shí)間的停機(jī)偏航狀態(tài)，以避免風(fēng)輪出現(xiàn)不穩(wěn)定狀態(tài)。

［參考文獻(xiàn)］

［1］ 許立國(guó)， 任建宇. 海上風(fēng)力發(fā)電的現(xiàn)狀及展望［J］. 港工技術(shù)， 2022， 59（6）： 45-48.

XU L G， REN J Y. Present situation and prospects of offshore wind power［J］. Port engineering technology， 2022， 59（6）： 45-48.

［2］ 李志川， 胡鵬， 馬佳星， 等. 中國(guó)海上風(fēng)電發(fā)展現(xiàn)狀分析及展望［J］. 中國(guó)海上油氣， 2022， 34（5）： 229-236.

LI Z C， HU P， MA J X， et al. Analysis and prospect of offshore" "wind" power" development" "in" "China［J］." China offshore oil and gas， 2022， 34（5）： 229-236.

［3］ 李會(huì)康， 束志鵬， 劉升貴. 大型風(fēng)力發(fā)電機(jī)葉片強(qiáng)度分析及發(fā)展趨勢(shì)研究［J］. 電子技術(shù)與軟件工程， 2023（1）： 100-103.

LI H K， SHU Z P， LIU S G. Strength analysis and development trend research of large wind turbine blades［J］. Electronic technology amp; software engineering， 2023（1）： 100-103.

［4］ 王明軍. 風(fēng)電機(jī)組葉片設(shè)計(jì)與氣動(dòng)彈性問(wèn)題［C］//第六屆中國(guó)風(fēng)電后市場(chǎng)交流合作大會(huì)論文集. 天津， 2019： 32-38.

WANG M J. The design and aeroelasticity problem of wind turbine blade［C］//Proceedings of the 6th China Wind Power Aftermarket Exchange and Cooperation Conference.Tianjin， 2019： 32-38.

［5］ 孫麗萍， 王昊， 丁嬌嬌. 風(fēng)力發(fā)電機(jī)葉片的氣動(dòng)彈性及顫振研究綜述［J］. 液壓與氣動(dòng)， 2012（10）： 1-5.

SUN L P， WANG H， DING J J. A summary on the aeroelasticity and flutter of wind turbine blade［J］. Chinese hydraulics amp; pneumatics， 2012（10）： 1-5.

［6］ YU D O， KWON O J. Predicting wind turbine blade loads and aeroelastic response using a coupled CFD-CSD method［J］. Renewable energy， 2014， 70： 184-196.

［7］ CHEN C， ZHOU J W， LI F M， et al. Stall-induced vibrations analysis and mitigation of a wind turbine rotor at idling state： theory and experiment［J］. Renewable energy， 2022， 187： 710-727.

［8］ 黃俊東. 基于特征值分析的后掠葉片的氣彈阻尼分析方法研究［D］. 廣州： 廣東工業(yè)大學(xué)， 2020.

HUANG J D. Research of aeroelastic damping of swept-blade" "based" "on" "eigenvalue" "method［D］." "Guangzhou： Guangdong University of Technology， 2020.

［9］ 柯世堂， 陸曼曼， 吳鴻鑫， 等. 基于風(fēng)洞試驗(yàn)15 MW風(fēng)電葉片顫振后形態(tài)與能量圖譜研究［J］. 空氣動(dòng)力學(xué)學(xué)報(bào)， 2022， 40（4）： 169-180.

KE S T， LU M M， WU H X， et al. Experimental study on the post-flutter morphological chracteristics and energy dissipation of a 15 MW wind turbine blade［J］. Acta aerodynamica sinica， 2022， 40（4）： 169-180.

［10］ 白學(xué)宗， 安宗文， 侯運(yùn)豐， 等. 低速?zèng)_擊載荷下的風(fēng)電葉片剛度退化規(guī)律［J］. 太陽(yáng)能學(xué)報(bào)， 2022， 43（1）： 132-139.

BAI X Z， AN Z W， HOU Y F， et al. Stiffness degradation rule of wind turbine blade under low-velocity impact load［J］. Acta energiae solaris sinica， 2022， 43（1）： 132-139.

［11］ 楊景云， 王文韞， 戴巨川. 基于攝動(dòng)模態(tài)分析的風(fēng)電葉片動(dòng)力學(xué)特性研究［J］. 太陽(yáng)能學(xué)報(bào)， 2023， 44（11）： 231-238.

YANG J Y， WANG W Y， DAI J C. Research on dynamic characteristics of wind power blades based on perturbation mode analysis［J］. Acta energiae solaris sinica， 2023， 44（11）： 231-238.

［12］ BETZ A. Das maximum der theoretisch mglichen Ausnützung des windes durch windmotoren［J］. Zeitschrift für das gesamte turbinewesen， 1920， 26： 307-309.

［13］ GLAUERT" "H." "Airplane" "propellers［M］." "Aerodynamic Theory. Heidelberg： Springer， 1935： 169-360.

［14］ PRANDTL L. The essentials of fluid dynamics［M］. New York： Hafner Publishing Company， 1963.

［15］ PITT D M， PETERS D A. Theoretical prediction of dynamic inflow derivatives［J］. Vertica， 1981， 5（1）： 21-34.

［16］ 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.

［17］ SHABANA A A. Dynamics of multibody systems， second edition［M］. Cambridge： Cambridge University Press， 1998.

AEROELASTIC RESPONSE ANALYSIS OF LARGE FLEXIBLE

WIND TURBINE BLADES

Gao Long1，Lin Lihui1，Yang Wansheng2，Gao Erjie3，Zhu Lei1

（1. Xiamen SunRui Wind Power Technology Co.， Ltd.， Xiamen 361000， China；

2. CSSC Haizhuang Windpower Equipment Co.， Ltd.， Chongqing 401121， China；

3. Guangdong Haizhuang Offshore Windpower Research Center Co.， Ltd.， Zhanjiang 524000， China）

Abstract：The aeroelastic effect of flexible blades is becoming increasingly significant with the rapid development of wind turbines and the sharp increase in blade length. Based on blade element momentum（BEM） theory and multi-body dynamics theory， the aerodynamic structure coupling equation of 95 m blade is established， and the aerodynamic response of blades under extreme wind load is analyzed which include a detailed study of the vibration characteristics of the blades after flutter and the damping characteristics of the blades under different speed conditions. According to the research， negative damping oscillation is a key factor causing blade flutter instability. At the same time， the influence of density of air， pitch angle and other factors on the aeroelastic stability of large blades are analyzed by parallel comparison method.

Keywords：wind turbine blades； aeroelasticity； damping； stability analysis； BEM