摘 要：將改進(jìn)的流固耦合方法用于不同風(fēng)速下風(fēng)力機(jī)的載荷和響應(yīng)特性研究，其準(zhǔn)確性得到NREL數(shù)據(jù)和風(fēng)力機(jī)葉片模態(tài)振動(dòng)實(shí)驗(yàn)的驗(yàn)證。與傳統(tǒng)的流固耦合方法相比，該方法減少了動(dòng)網(wǎng)格使用數(shù)量、避免負(fù)體積網(wǎng)格，提高了計(jì)算穩(wěn)定性，將計(jì)算時(shí)間縮短約50%。仿真結(jié)果表明，將風(fēng)速由25 m/s提升至35 m/s，葉片的推力和扭矩載荷的波動(dòng)振幅分別增加了6.8倍和9.8倍。同時(shí)，在葉片吸力面發(fā)現(xiàn)存在結(jié)構(gòu)屈曲，這與臺(tái)風(fēng)現(xiàn)場(chǎng)葉片斷裂位置接近，可能是風(fēng)力機(jī)葉片結(jié)構(gòu)在臺(tái)風(fēng)環(huán)境下的強(qiáng)度薄弱點(diǎn)，同時(shí)，在葉片吸力面發(fā)現(xiàn)存在結(jié)構(gòu)屈曲，這與臺(tái)風(fēng)現(xiàn)場(chǎng)葉片斷裂位置接近，是可能的風(fēng)力機(jī)葉片結(jié)構(gòu)在臺(tái)風(fēng)環(huán)境下的強(qiáng)度薄弱點(diǎn)。

關(guān)鍵詞：風(fēng)力機(jī)葉片；臺(tái)風(fēng)環(huán)境；振動(dòng)環(huán)境；流固耦合；數(shù)值模擬；葉片強(qiáng)度

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

0 引 言

近年來，由于臺(tái)風(fēng)等極端天氣的影響，許多海上和沿海風(fēng)電場(chǎng)都遭受了重大的經(jīng)濟(jì)損失。隨著風(fēng)力機(jī)容量和尺寸的不斷增加，在臺(tái)風(fēng)環(huán)境下的葉片斷裂和損壞等事故時(shí)有發(fā)生。例如，在2013年，超級(jí)臺(tái)風(fēng)天兔對(duì)沿海風(fēng)電場(chǎng)造成災(zāi)難性破壞，導(dǎo)致35個(gè)葉片斷裂和8個(gè)塔架坍塌［1］。2015年8月，強(qiáng)臺(tái)風(fēng)蘇迪羅導(dǎo)致6臺(tái)2 MW的風(fēng)力機(jī)失效，1臺(tái)風(fēng)電機(jī)葉片斷裂［2］。

為減少事故率［3］，減少運(yùn)維成本［4］，學(xué)者們針對(duì)臺(tái)風(fēng)造成的風(fēng)力機(jī)破壞事故進(jìn)行了調(diào)查研究。陳曉等［1］對(duì)2013年天兔臺(tái)風(fēng)造成的風(fēng)電場(chǎng)破壞事故進(jìn)行全面調(diào)查，收集并分析當(dāng)?shù)仫L(fēng)場(chǎng)、塔架和葉片的數(shù)據(jù)，通過將風(fēng)載荷計(jì)算與結(jié)構(gòu)建模相結(jié)合，建立系統(tǒng)的計(jì)算方法，并給出風(fēng)力機(jī)相關(guān)設(shè)計(jì)標(biāo)準(zhǔn)的修改建議；文獻(xiàn)［2］使用有限元方法分析了塔架倒塌和葉片斷裂的機(jī)械機(jī)理，確定了風(fēng)力機(jī)的失效位置和弱點(diǎn)，并為未來風(fēng)力發(fā)電機(jī)的風(fēng)險(xiǎn)防范和抗風(fēng)設(shè)計(jì)提供了依據(jù)；柯世堂等［5］研究臺(tái)風(fēng)“威馬遜”環(huán)境下海上風(fēng)電場(chǎng)風(fēng)荷載特性，并找到最不利停機(jī)位置；王立忠等［6］對(duì)近海風(fēng)力機(jī)發(fā)電機(jī)關(guān)鍵結(jié)構(gòu)在臺(tái)風(fēng)中的災(zāi)變研究進(jìn)行綜述，認(rèn)為準(zhǔn)確模擬臺(tái)風(fēng)及其誘發(fā)的極端海浪特征，清楚風(fēng)力機(jī)組在極端臺(tái)風(fēng)浪載荷作用下的關(guān)鍵結(jié)構(gòu)累積變形物理機(jī)制，是保證風(fēng)機(jī)結(jié)構(gòu)安全的重點(diǎn)。

所以，為深入研究風(fēng)力機(jī)振動(dòng)、形變和失效的根本原因，采用非線性動(dòng)力學(xué)或有限元的方法對(duì)風(fēng)力機(jī)自身結(jié)構(gòu)強(qiáng)度進(jìn)行分析是必要的。例如通過比較在擺動(dòng)載荷下葉片的幾何非線性結(jié)構(gòu)響應(yīng)和線性結(jié)構(gòu)響應(yīng)［7］，可確定葉片在極限強(qiáng)度和穩(wěn)定極限狀態(tài)下的抗力，揭示實(shí)際的故障機(jī)制和關(guān)鍵模式。在此基礎(chǔ)上，Nelson等［8］使用有限元方法建立一個(gè)風(fēng)力機(jī)的鋼框架結(jié)構(gòu)模型，并在極端環(huán)境載荷下進(jìn)行測(cè)試，分析了臺(tái)風(fēng)條件下近海風(fēng)力機(jī)的基礎(chǔ)強(qiáng)度。該研究探索了風(fēng)力機(jī)結(jié)構(gòu)動(dòng)態(tài)響應(yīng)與風(fēng)載荷持續(xù)時(shí)間之間的強(qiáng)相關(guān)性。發(fā)現(xiàn)載荷持續(xù)時(shí)間減少25%會(huì)導(dǎo)致動(dòng)態(tài)應(yīng)力增加10%。安利強(qiáng)等［4］基于葉素動(dòng)量理論研究臺(tái)風(fēng)作用下風(fēng)力機(jī)葉根和塔底的極限載荷，結(jié)果表明，風(fēng)力機(jī)和外部環(huán)境關(guān)鍵參數(shù)的隨機(jī)性是造成臺(tái)風(fēng)過境后風(fēng)電場(chǎng)機(jī)組出現(xiàn)不一致破壞狀態(tài)的主要原因。然而，上述研究中的風(fēng)力機(jī)氣動(dòng)載荷計(jì)算部分均存在大量假設(shè)，葉片表面的流動(dòng)分離和脫落渦被簡化甚至忽略。然而，參考熱帶氣旋區(qū)風(fēng)力機(jī)設(shè)計(jì)的湍流旋風(fēng)模型，高風(fēng)速、湍流強(qiáng)度和風(fēng)剪切是影響風(fēng)力機(jī)生存的重要因素［9］。所以，進(jìn)行更為精細(xì)化的風(fēng)力機(jī)流場(chǎng)氣動(dòng)載荷仿真，用于風(fēng)力機(jī)結(jié)構(gòu)應(yīng)力應(yīng)變分析，對(duì)于提高風(fēng)力機(jī)在臺(tái)風(fēng)環(huán)境下的可靠性非常必要。

將高精度的計(jì)算流體力學(xué)方法（computational fluid dynamics，CFD）與有限元方法相結(jié)合的流固耦合計(jì)算方法（fluid-structure interaction，F(xiàn)SI）是一種更準(zhǔn)確的風(fēng)力機(jī)結(jié)構(gòu)響應(yīng)求解方法［10-11］。Carrión等［12］利用FSI的概念，使用靜態(tài)氣彈性方法和動(dòng)態(tài)氣彈性方法分析了兩臺(tái)風(fēng)力機(jī)的氣動(dòng)彈性，并發(fā)現(xiàn)葉片擺動(dòng)在風(fēng)力機(jī)的結(jié)構(gòu)整體變形中占主導(dǎo)位置。另一項(xiàng)研究中，繆維跑等［13］使用高精度的FSI方法發(fā)現(xiàn)扭轉(zhuǎn)葉片結(jié)合耦合材料可有效地減少由于葉片形變引起的氣動(dòng)載荷和最大應(yīng)力，進(jìn)而顯著降低臺(tái)風(fēng)環(huán)境下的負(fù)載。然而，使用FSI會(huì)消耗大量計(jì)算資源［14］，雙向FSI還會(huì)引起網(wǎng)格負(fù)體積問題，造成求解困難。因此，需對(duì)傳統(tǒng)FSI的使用策略進(jìn)行優(yōu)化改進(jìn)，使其具有更高的仿真效率和計(jì)算穩(wěn)定性。

綜上，本文以在臺(tái)風(fēng)期間發(fā)生葉片結(jié)構(gòu)破壞的5 MW風(fēng)力機(jī)（如圖1所示）為研究對(duì)象，以一種改進(jìn)的雙向流固耦合方法（improved fluid-structure interaction，IFSI）對(duì)其風(fēng)輪結(jié)構(gòu)進(jìn)行瞬態(tài)結(jié)構(gòu)響應(yīng)分析。該方法可顯著減少計(jì)算時(shí)間，提高計(jì)算穩(wěn)定性。通過與葉片振動(dòng)實(shí)驗(yàn)和美國國家可再生能源實(shí)驗(yàn)室（National Renewable Energy Laboratory，NREL）數(shù)據(jù)比對(duì)，證明了IFSI方法的可行性?；诖朔椒?，研究臺(tái)風(fēng)環(huán)境下風(fēng)力機(jī)的載荷和響應(yīng)特性，并發(fā)現(xiàn)了風(fēng)力機(jī)結(jié)構(gòu)的薄弱點(diǎn)，為后續(xù)風(fēng)力機(jī)設(shè)計(jì)優(yōu)化提供了參考。

1 風(fēng)力機(jī)模型

以一臺(tái)5 MW風(fēng)力機(jī)作為研究對(duì)象，葉片長度[R]為63 m。用兩個(gè)重疊的流體域（見圖2）來模擬風(fēng)流過風(fēng)輪轉(zhuǎn)子的整個(gè)過程。外流域?yàn)閳A柱狀，風(fēng)的進(jìn)口和出口分別位于轉(zhuǎn)子上游[5R]和下游[10R]處，半徑為[2R]。內(nèi)流域長度為[2R]，半徑為[1.2R]。風(fēng)力機(jī)模型位于內(nèi)流域中心。為盡量減小流場(chǎng)尺寸對(duì)風(fēng)力機(jī)的影響，流場(chǎng)外部邊緣面被設(shè)定為對(duì)稱面，并將內(nèi)部交界面為重疊截面，進(jìn)行數(shù)據(jù)傳輸。

在本研究中，默認(rèn)風(fēng)速是影響風(fēng)力機(jī)在極端環(huán)境運(yùn)行的主要因素［15］。選擇3個(gè)入口風(fēng)速，分別代表熱帶低氣壓風(fēng)速（15 m/s）、強(qiáng)熱帶風(fēng)暴風(fēng)速（25 m/s）和臺(tái)風(fēng)風(fēng)速（35 m/s）。在本研究中，風(fēng)速被假定為均勻的，中等湍流水平，5%的湍流強(qiáng)度和10%的湍流粘性比。設(shè)計(jì)6個(gè)計(jì)算仿真案例，分別為使用剛性的CFD方法（R1、R2和R3）和柔性的IFSI方法（F1、F2和F3），具體細(xì)節(jié)見表1。

理論上，葉片的槳距角應(yīng)隨風(fēng)速變化而改變，以保證發(fā)電功率穩(wěn)定。然而，為減小葉片轉(zhuǎn)動(dòng)引起的網(wǎng)格畸變問題，本研究中以限制葉片自由度的方法將槳距角設(shè)為13.308°。此角度是根據(jù)文獻(xiàn)［16］的A2設(shè)定，并與圖1中風(fēng)力機(jī)的槳距角相似。

2 IFSI方法

2.1 變形和運(yùn)動(dòng)分解原理

IFSI方法與傳統(tǒng)的雙向流固耦合計(jì)算具有近似的流程步驟，即CFD計(jì)算流場(chǎng)氣動(dòng)載荷、有限元計(jì)算結(jié)構(gòu)應(yīng)力應(yīng)變、流場(chǎng)網(wǎng)格更新和時(shí)間步更新。但傳統(tǒng)FSI通過動(dòng)網(wǎng)格方法實(shí)現(xiàn)流場(chǎng)網(wǎng)格運(yùn)動(dòng)和變形（見圖3a），隨著時(shí)間步增加或運(yùn)動(dòng)幅度過大時(shí)，網(wǎng)格質(zhì)量會(huì)嚴(yán)重下降，極易出現(xiàn)負(fù)體積現(xiàn)象，導(dǎo)致計(jì)算中止［17］。同時(shí)，每個(gè)時(shí)間步中大量的動(dòng)網(wǎng)格更新，也會(huì)增加仿真，提高研究成本。

針對(duì)以上問題，IFSI采用變形和運(yùn)動(dòng)分解的思路，根據(jù)不同的變形和運(yùn)動(dòng)形式劃分流場(chǎng)網(wǎng)格，可減少動(dòng)網(wǎng)格使用，進(jìn)而減少仿真計(jì)算時(shí)長和避免負(fù)體積網(wǎng)格出現(xiàn)。如圖3b所示，在IFSI中，風(fēng)力機(jī)風(fēng)輪整體的變形和運(yùn)動(dòng)可看作是整體的勻速旋轉(zhuǎn)運(yùn)動(dòng)和葉片局部變形的疊加。其中，風(fēng)力機(jī)風(fēng)輪的變形主要集中在葉片結(jié)構(gòu)上，有揮舞、擺振和扭轉(zhuǎn)3種形式，通常為局部變形，所以，動(dòng)網(wǎng)格的使用區(qū)域較小。且由于葉片變形幅度有限，有效避免了由于網(wǎng)格過量拉伸與壓縮造成的質(zhì)量惡化，減少了計(jì)算迭代過程中出現(xiàn)負(fù)體積網(wǎng)格的可能性。

風(fēng)輪的旋轉(zhuǎn)運(yùn)動(dòng)可看作風(fēng)力機(jī)風(fēng)輪整體運(yùn)動(dòng)和變形的背景，通過滑移網(wǎng)格法在圖2的內(nèi)流域設(shè)置風(fēng)輪轉(zhuǎn)速實(shí)現(xiàn)。由于每個(gè)位置的轉(zhuǎn)速均相同，故風(fēng)輪的旋轉(zhuǎn)運(yùn)動(dòng)可看作是一種剛性運(yùn)動(dòng)，在計(jì)算中無需使用動(dòng)網(wǎng)格方法。

綜上，IFSI方法可通過風(fēng)力機(jī)風(fēng)輪變形和運(yùn)動(dòng)的分解原理，減少動(dòng)網(wǎng)格的使用范圍，保證了網(wǎng)格更新后質(zhì)量，有利于提高仿真效率和計(jì)算精度。

2.2 CFD方法

整個(gè)流場(chǎng)的物理特性符合N-S方程的基本數(shù)學(xué)描述：

[?ui?xi=0]""""""""" （1）

[ρ?ui?t+?uiuj?xi=??xiμ+μt?ui?xi-?p?xi]" （2）

式中：[t]——空氣流動(dòng)時(shí)間，s；[u]——流速，m/s；[p]——空氣壓力，Pa；[ρ]——流體密度，kg/m3；[μ]——?jiǎng)恿︷ざ认禂?shù)，（N·s）/m2；[μt]——湍流黏度，m2/s。在此基礎(chǔ)上使用雷諾平均（Reynolds averaging Navier-Stokes，RANS）方程對(duì)湍流波動(dòng)進(jìn)行時(shí)間平均而得到多種湍流模型：

[?ui?t+uj?ui?xj=fi-1ρ?p?xi+v?2ui?xj?xi-?u′iu′j??ui?xj=0]" （3）

SST k-ω模型綜合了標(biāo)準(zhǔn)k-ω和標(biāo)準(zhǔn)k-ω模型的優(yōu)點(diǎn)。它通過在黏性子層和遠(yuǎn)場(chǎng)流動(dòng)之間進(jìn)行動(dòng)態(tài)的模型轉(zhuǎn)換，使得遠(yuǎn)場(chǎng)流動(dòng)動(dòng)力學(xué)的求解成為可能。因此，SST k-ω模型在風(fēng)力機(jī)數(shù)值模擬研究中被廣泛使用［18，13］。在本文同樣選擇SST k-ω模型來模擬風(fēng)力機(jī)的流場(chǎng)分布。

2.3 流場(chǎng)網(wǎng)格

使用重疊網(wǎng)格技術(shù)實(shí)現(xiàn)內(nèi)外流域的數(shù)據(jù)交互（圖4a）。外流域負(fù)責(zé)風(fēng)力機(jī)的遠(yuǎn)場(chǎng)模擬，其具有較大的空間尺寸，因此使用較為稀疏的六面體網(wǎng)格以減少計(jì)算時(shí)間。另一方面，內(nèi)流域需詳細(xì)捕捉葉片周圍流場(chǎng)的細(xì)節(jié)，并滿足動(dòng)網(wǎng)格的兼容性要求。因此，內(nèi)流域使用非結(jié)構(gòu)體網(wǎng)格，其網(wǎng)格大小受到葉片表面網(wǎng)格大小的限制。離葉片表面越近網(wǎng)格約密集（如圖4b所示）。但由于使用風(fēng)輪變形和運(yùn)動(dòng)分解后，內(nèi)流域的動(dòng)網(wǎng)格區(qū)域較少，所以對(duì)內(nèi)部精細(xì)網(wǎng)格的需求減少，網(wǎng)格總數(shù)量降低，同時(shí)更易保證網(wǎng)格質(zhì)量。

建立葉片表明的邊界層網(wǎng)格需估算在葉片表面的流體速度，葉片表面邊界層的第1層網(wǎng)格高度設(shè)置為0.008 mm，[y+]約為1。該高度經(jīng)過模擬驗(yàn)證，滿足計(jì)算要求。

外流域網(wǎng)格保持不變，對(duì)內(nèi)流域網(wǎng)格進(jìn)行無關(guān)性驗(yàn)證。通過修改葉片的線網(wǎng)格尺寸[L]來控制網(wǎng)格的生成（如圖5a），取[L]的值分別為0.40、0.30、0.20、0.10、0.08、0.07、0.06和0.04 m。在標(biāo)準(zhǔn)工況下，風(fēng)速為11.4 m/s，葉片旋轉(zhuǎn)速度為12.1 r/min，比較不同網(wǎng)格情況下計(jì)算出的推力和扭矩值。綜合考慮模擬的準(zhǔn)確性和網(wǎng)格計(jì)算的成本，在本研究中選擇[L=0.07 m]，網(wǎng)格數(shù)為868萬，推力和扭矩值分別為785 kN和4390 kN·m（如圖5b所示）。

NREL使用FAST軟件得到的標(biāo)準(zhǔn)條件下5 MW風(fēng)力機(jī)模擬結(jié)果［16］可用作準(zhǔn)確性對(duì)比，推力和扭矩值分別為813.6 kN和4190.4 kN·m。與FAST計(jì)算結(jié)果相比較，推力和扭矩的差異分別為3.52%和4.76%，模擬計(jì)算的方法和風(fēng)力機(jī)結(jié)構(gòu)差異是導(dǎo)致計(jì)算結(jié)果不同的主要原因。FAST采用葉素動(dòng)量（blade element momentum，BEM）理論簡化了葉片在展向方向上的結(jié)構(gòu)［19］，將風(fēng)力視為載荷，并通過疊加獲得最終結(jié)果。相反，本研究中采用的CFD方法未在流體部分進(jìn)行許多假設(shè)和簡化，充分考慮了空氣的展向流動(dòng)、葉片表面的黏性和邊界層分離。此外，NREL計(jì)算整個(gè)風(fēng)力機(jī)結(jié)構(gòu)，而本研究只關(guān)注風(fēng)力機(jī)葉片轉(zhuǎn)子系統(tǒng)，忽略了塔架和機(jī)艙。

2.4 葉片結(jié)構(gòu)仿真和實(shí)驗(yàn)

風(fēng)力機(jī)的結(jié)構(gòu)響應(yīng)是使用有限元法（finite element method，F(xiàn)EM）計(jì)算的。風(fēng)力機(jī)的葉片和腹板采用單向E-玻璃纖維/環(huán)氧復(fù)合材料制成，材料性質(zhì)與文獻(xiàn)［13］中一致。這種結(jié)構(gòu)與材料確保了葉片的機(jī)械強(qiáng)度，同時(shí)降低了重量（圖6）。

通過葉片振動(dòng)特性實(shí)驗(yàn)驗(yàn)證FEM方法的準(zhǔn)確性。使用玻璃纖維和碳纖維制造3個(gè)成比例縮小的實(shí)驗(yàn)葉片（如圖7a，[R]為2.5 m），使用DH5923N收集加速度信號(hào)和DYTRAN以瞬態(tài)激勵(lì)方法測(cè)量擺動(dòng)方向的模態(tài)頻率，振動(dòng)頻率采集范圍設(shè)置在0～500 Hz。激勵(lì)點(diǎn)在葉片頂端，激勵(lì)方向垂直于激勵(lì)點(diǎn)處的葉片表面。使用蜂蠟對(duì)稱粘貼加速度傳感器，測(cè)量點(diǎn)如圖7b所示。使用力錘手動(dòng)激勵(lì)葉片，并收集、保存、分析和模態(tài)分析所得的加速度數(shù)據(jù)。使用有限元方法對(duì)葉片模型進(jìn)行了仿真，并將前五階振動(dòng)特性與振動(dòng)實(shí)驗(yàn)獲得的振型進(jìn)行了比較見表2，有限元法所得到的葉片模型振動(dòng)特性與實(shí)驗(yàn)結(jié)果相似，證明了有限元方法的準(zhǔn)確性。

2.5 IFSI計(jì)算

在剛性案例下（R1、R2、R3），風(fēng)力機(jī)整體形狀保持不變，只發(fā)生勻速轉(zhuǎn)動(dòng)，因此采用滑動(dòng)網(wǎng)格模型（SMM）處理內(nèi)部流體場(chǎng)中的網(wǎng)格運(yùn)動(dòng)。而在柔性計(jì)算案例（F1、F2、F3）中則需采用流固耦合方法，模擬由于空氣載荷導(dǎo)致的葉片變形和振動(dòng)。所以方法的關(guān)鍵在于對(duì)流場(chǎng)網(wǎng)格進(jìn)行及時(shí)更新，同時(shí)要保證動(dòng)網(wǎng)格的質(zhì)量，避免由于網(wǎng)格過于拉伸或壓縮導(dǎo)致負(fù)體積。通常有3種不同的網(wǎng)格更新方法：Smoothing、Layering和Remeshing。

Smoothing特別適用于移動(dòng)邊界的平移，網(wǎng)格的運(yùn)動(dòng)是通過求解擴(kuò)散方程來確定的。

[▽?γ▽u=0γ=1dα]"" （4）

式中：[u]——網(wǎng)格運(yùn)動(dòng)的速度，m/s；[γ]——擴(kuò)散系數(shù)，用于控制邊界運(yùn)動(dòng)對(duì)內(nèi)部網(wǎng)格變形的影響；[d]——法向邊界距離，m；[α]——擴(kuò)散參數(shù)，通過許多研究嘗試，[α]被設(shè)定為1.7較為合理。

Layering是更新六面體或楔形網(wǎng)格的常用方法。它基于邊界運(yùn)動(dòng)來進(jìn)行網(wǎng)格合并和分割。當(dāng)與邊界相鄰的網(wǎng)格層高度達(dá)到指定閾值時(shí)，網(wǎng)格被分成兩層。同樣地，當(dāng)網(wǎng)格層高度降低到一定水平時(shí)，與邊界相鄰的兩層網(wǎng)格合并成為一層。該技術(shù)特別適用于邊界層網(wǎng)格更新。網(wǎng)格大小為：

[hmaxgt;1+αshidealhminlt;αchideal]"""""" （5）

式中：[hmax]——網(wǎng)格層的最大高度，m；[hmin]——網(wǎng)格層的最小高度，m；[hideal]——理想網(wǎng)格高度，m；[αs]——分裂因子；[αc]——折疊因子。

在大變形情況下，僅使用Smoothing和Layering可能會(huì)導(dǎo)致網(wǎng)格質(zhì)量不佳，甚至在網(wǎng)格更新迭代過程中出現(xiàn)負(fù)體積網(wǎng)格，導(dǎo)致計(jì)算中斷。為解決這個(gè)問題，對(duì)于內(nèi)部流場(chǎng)的大尺寸網(wǎng)格位移，采用Remeshing方法。然而，它僅適用于四面體非結(jié)構(gòu)網(wǎng)格。在本研究中，為提高網(wǎng)格更新質(zhì)量，同時(shí)使用了這3種方法。

在IFSI中，充分考量了常規(guī)FSI方法和CFD方法的優(yōu)缺點(diǎn)，對(duì)以下3個(gè)內(nèi)容進(jìn)行改進(jìn)提升：1）使用CFD穩(wěn)態(tài)和瞬態(tài)數(shù)值計(jì)算結(jié)果作為流固耦合初始化，以減少流場(chǎng)局部的速度和壓力畸變，是保證流固耦合計(jì)算初始階段穩(wěn)定的關(guān)鍵；2）3種動(dòng)網(wǎng)格方法配合使用保證了計(jì)算過程中的網(wǎng)格質(zhì)量；3）將風(fēng)力機(jī)轉(zhuǎn)子的運(yùn)動(dòng)分解為剛性的轉(zhuǎn)動(dòng)和柔性的振動(dòng)，參考剛性案例以滑移網(wǎng)格處理轉(zhuǎn)子圓周運(yùn)動(dòng)，減少動(dòng)網(wǎng)格使用頻率。IFSI的工作流程如圖8所示。

表3展示了CFD、傳統(tǒng)FSI和IFSI這3種方法的統(tǒng)計(jì)信息。CFD的計(jì)算耗時(shí)最少，傳統(tǒng)FSI方法的計(jì)算時(shí)間成本最高。經(jīng)過多次比較，IFSI方法的計(jì)算時(shí)間比FSI方法降低了52%，且不存在負(fù)體積網(wǎng)格，極大減少了研究人員和工程師的工作量和調(diào)試程序的難度。

3 結(jié)果與討論

3.1 風(fēng)力機(jī)載荷分析

6個(gè)案例中風(fēng)力機(jī)的推力和扭矩如圖9所示。在圖9a中，當(dāng)風(fēng)速為15 m/s時(shí)，推力和扭矩的趨勢(shì)幾乎相同。在R1的初始計(jì)算階段，推力和扭矩的振幅最大，但隨著時(shí)間的推移逐步穩(wěn)定在922.5和6825 kN·m。在風(fēng)速為25 m/s的情況下（圖9b），推力和扭矩相對(duì)于R1有顯著增加，推力和扭矩的最大值分別為1046.6和7534.8 kN·m。推力和扭矩的趨勢(shì)分布相似，在短時(shí)間內(nèi)均有明顯幅值脈動(dòng)。在風(fēng)速為35 m/s時(shí)，R3（圖9c）具有更高的推力值（最大值1463.7 kN；最小值1429.6 kN），但扭矩相對(duì)于風(fēng)速為25 m/s時(shí)有所降低（最大值7853.6 kN·m；最小值7211.8 kN·m）。

與剛性情況相比，使用IFSI方法計(jì)算的推力和扭矩都有所降低。在圖9d中，推力和扭矩的趨勢(shì)相似，振動(dòng)波形呈遞減的山峰形狀（箭頭標(biāo)識(shí)）。在25 m/s的風(fēng)速下（圖9e），起始段的振動(dòng)形狀與F1的趨勢(shì)相同，但極值更高。中間段和后段的推力和扭矩都表現(xiàn)出明顯的周期性（由矩形方框和橢圓形框標(biāo)識(shí)）。在35 m/s的風(fēng)速下，F(xiàn)3情況表現(xiàn)出最強(qiáng)烈的波動(dòng)，如圖9f所示。在顯著的上升后，推力和扭矩的值逐漸穩(wěn)定在一定范圍內(nèi)，并呈周期性波動(dòng)。推力和扭矩的振幅分別達(dá)到11.93 kN和127.68 kN·m。

定義[K]來量化比較風(fēng)力機(jī)推力和扭矩在短時(shí)間內(nèi)波動(dòng)的劇烈程度。

[K=maxtmax/tminτend-τstart]"" （6）

式中：[tmax]、[tmin]——在最劇烈振動(dòng)部分（在圖中用箭頭標(biāo)注）的推力或扭矩的最大值和最小值，kN和kN·m；[τstart]、[τend]——最劇烈振動(dòng)部分的開始時(shí)間和結(jié)束時(shí)間，s。推力和扭矩的[K]值（簡稱[Ktr]和[Kto]，單位均為[s-1]，只需考慮相對(duì)大小）如表4所示。可發(fā)現(xiàn)R1的[K]值最?。╗Ktr=1.62;Kto=1.22]），其次是R2（[Ktr=1.49;Kto=2.58]）和R3（[Ktr=3.52;Kto=3.18]）。結(jié)果表明，R3經(jīng)歷了最顯著的流場(chǎng)波動(dòng)，這可能是由于風(fēng)力機(jī)后方的氣流出現(xiàn)了大規(guī)模的流動(dòng)分離［20］，引起上游流場(chǎng)波動(dòng)，導(dǎo)致強(qiáng)烈的負(fù)載變化。

流固耦合的情況下，[Ktr]和[Kto]在15 m/s風(fēng)速時(shí)分別降至0.92和0.64。F2的[Kto]為1.68，在25 m/s風(fēng)速下與R2相比略有下降，[Ktr]基本相同。F3在35 m/s風(fēng)速下可能是6種情況中最危險(xiǎn)的條件，具有最大的[Ktr]和[Kto]（[Ktr=7.20]；[Kto=6.78]）。與低風(fēng)速的15 m/s相比，35 m/s的推力和扭矩振動(dòng)強(qiáng)度[K]分別增加了6.8倍和9.8倍。

3.2 空氣流動(dòng)特性分析

使用λ2標(biāo)準(zhǔn)等值面顯示不同風(fēng)速下風(fēng)力機(jī)周圍的流場(chǎng)特征（如圖10），可看到葉片與尾流渦街之間的相互作用。在低風(fēng)速下，脫落渦體積較少，主體為葉片的翼尖渦。在25 m/s時(shí)，可在圖10b中觀察到葉片表面的流動(dòng)分離。同時(shí)由于主流方向的高流速導(dǎo)致大的螺旋間距，使得尾流渦相對(duì)較早地到達(dá)流場(chǎng)網(wǎng)格的粗糙部分并過早地消散［12］。在35 m/s（圖10c）下，螺旋流的轉(zhuǎn)數(shù)進(jìn)一步減少，并在葉片后形成了一個(gè)類似于扇形的流動(dòng)面（圖中箭頭），表明流動(dòng)葉片表明的流動(dòng)分離進(jìn)一步惡化。

流固耦合條件下，葉片表面的渦結(jié)構(gòu)，如圖11所示。在15 m/s的情況下，渦旋體積和葉片彎曲曲率較小，幾乎無流動(dòng)分離的跡象。風(fēng)力機(jī)葉片保持在最佳性能狀態(tài)下。在25和35 m/s的情況下，葉片的曲率小幅增加，流動(dòng)分離在高風(fēng)速下變得越來越明顯，特別是在35 m/s的情況下流動(dòng)分離現(xiàn)象嚴(yán)重。

3.3 葉片結(jié)構(gòu)響應(yīng)

風(fēng)力機(jī)葉片的結(jié)構(gòu)響應(yīng)如圖12所示。3個(gè)案例擁有一致的葉片擺振過程，但隨風(fēng)速增加最大變形增大。風(fēng)力機(jī)葉片的擺振變形可分為以下3個(gè)階段：1）在風(fēng)壓作用下，葉片順著下風(fēng)方向彎曲和變形。隨著葉片的彎曲，受風(fēng)壓力作用的面積減少，同時(shí)根據(jù)胡克定律，葉片的彈性增加；2）葉尖快速向后移動(dòng)，直到最大位移處。此時(shí)，葉片的變形速度降至零，葉片回彈力大于風(fēng)壓力，導(dǎo)致葉片回彈；3）當(dāng)回彈速度降至零時(shí)，回彈結(jié)束。在風(fēng)壓力的驅(qū)動(dòng)下，下一個(gè)葉片擺振周期開始。

葉片表面上的等效應(yīng)力如圖13所示。葉片上的等效應(yīng)力分布趨勢(shì)在3種風(fēng)速下是一致的，在吸力面最大弦長處出現(xiàn)應(yīng)力集中點(diǎn)。隨著風(fēng)速的增加，葉片表面上的應(yīng)力集中現(xiàn)象也越明顯。

應(yīng)力集中點(diǎn)可能是風(fēng)力機(jī)葉片在極端條件下（如臺(tái)風(fēng)）失效的關(guān)鍵位置點(diǎn)。如圖14所示，吸力面的應(yīng)力集中位置存在明顯的結(jié)構(gòu)屈曲，這一位置與圖1中的葉片斷裂位置基本一致。在風(fēng)力機(jī)葉片的設(shè)計(jì)和制造過程中加強(qiáng)這一區(qū)域可能會(huì)提高風(fēng)力機(jī)的抗風(fēng)能力和耐久性。

4 結(jié) 論

本研究使用改進(jìn)的流固耦合方法（IFSI）研究在臺(tái)風(fēng)環(huán)境下風(fēng)力機(jī)的載荷和響應(yīng)特性。以一臺(tái)在臺(tái)風(fēng)中葉片失效的5 MW風(fēng)力機(jī)為研究對(duì)象，并通過NREL數(shù)據(jù)和葉片振動(dòng)實(shí)驗(yàn)進(jìn)行對(duì)比驗(yàn)證，主要結(jié)論如下：

1）與傳統(tǒng)的流固耦合方法相比較，IFSI方法基于風(fēng)力機(jī)變形和運(yùn)動(dòng)分解原理，通過減少動(dòng)網(wǎng)格使用，保證了網(wǎng)格更新質(zhì)量，有效地避免了負(fù)體積網(wǎng)格發(fā)生，顯著提高了仿真效率。每次迭代所需的平均時(shí)間從590 s減少到308.8 s，約為傳統(tǒng)流固耦合方法的52%。

2）模擬結(jié)果表明，高風(fēng)速（35 m/s）會(huì)增加葉片上的推力和轉(zhuǎn)矩，導(dǎo)致葉片變形增大，提高了結(jié)構(gòu)失效的風(fēng)險(xiǎn)。并通過引入[K]值來定量化比較風(fēng)力機(jī)推力和扭矩在短時(shí)間內(nèi)波動(dòng)的劇烈程度，發(fā)現(xiàn)在高風(fēng)速時(shí)葉片的推力和轉(zhuǎn)矩載荷的波動(dòng)幅度分別增加了6.8和9.8倍。

3）每個(gè)葉片的吸力面在靠近根部附近均有一個(gè)應(yīng)力集中點(diǎn)，與在風(fēng)電場(chǎng)的實(shí)際調(diào)查結(jié)果中葉片斷裂位置對(duì)應(yīng)，故可能是葉片在臺(tái)風(fēng)環(huán)境下的強(qiáng)度薄弱位置。這一發(fā)現(xiàn)可為后續(xù)的風(fēng)力機(jī)結(jié)構(gòu)設(shè)計(jì)和優(yōu)化提供參考。

［參考文獻(xiàn)］

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

［2］"""""" CHOU J S， OU Y C， LIN K Y， et al. Structural failure simulation of onshore wind turbines impacted by strong winds［J］. Engineering structures， 2018， 162： 257-269.

［3］"""""" DAO C， KAZEMTABRIZI B， CRABTREE C. Wind turbine reliability data review and impacts on levelised cost of energy［J］. Wind energy， 2019， 22（12）： 1848-1871.

［4］"""""" 安利強(qiáng)， 孫陽， 王鵬， 等. 考慮隨機(jī)參數(shù)的風(fēng)電機(jī)組臺(tái)風(fēng)載荷特征研究［J］. 太陽能學(xué)報(bào)， 2021， 42（4）： 417-423.

AN L Q， SUN Y， WANG P， et al. Study on typhoon load characteristics of wind turbines considering random parameters［J］. Acta energiae solaris sinica， 2021， 42（4）： 417-423.

［5］"""""" 柯世堂， 王碩， 趙永發(fā)， 等. 臺(tái)風(fēng)-浪-流耦合作用下海上10 MW級(jí)特大型風(fēng)力機(jī)風(fēng)荷載特性分析［J］. 振動(dòng)工程學(xué)報(bào)， 2023， 36（2）： 299-310.

KE S T， WANG S， ZHAO Y F， et al. Wind load characteristics of 10 MW-level super-large offshore wind turbine under the coupling effect of typhoon-wave-current［J］. Journal of vibration engineering， 2023， 36（2）： 299-310.

［6］"""""" 王立忠， 洪義， 高洋洋， 等. 近海風(fēng)電結(jié)構(gòu)臺(tái)風(fēng)環(huán)境動(dòng)力災(zāi)變與控制［J］. 力學(xué)學(xué)報(bào)， 2023， 55（3）： 567-587.

WANG L Z， HONG Y， GAO Y Y， et al. Dynamic catastrophe and control of offshore wind power structures in typhoon"" environment1）［J］.""" Chinese"" journal""" of theoretical and applied mechanics， 2023， 55（3）： 567-587.

［7］"""""" ROSEMEIER M， BERRING P， BRANNER K. Non-linear ultimate strength and stability limit state analysis of a wind turbine blade［J］. Wind energy， 2016， 19（5）： 825-846.

［8］"""""" NELSON B， LIN T Y， QUéMéNER Y， et al. Extreme typhoon loads effect on the structural response of an offshore wind turbine［C］//Proceedings of 7th PAAMES and AMEC2016， Hong Kong， China， 2016， 13： 14.

［9］"""""" HAN T， MCCANN G， MüCKE T A， et al. How can a wind turbine survive in tropical cyclone？［J］. Renewable energy， 2014， 70： 3-10.

［10］"""" ZHANG R X， HE L， AN L Q. A dynamic probabilistic analysis method for wind turbine rotor based on the surrogate model［J］. Journal of renewable and sustainable energy， 2023， 15（1）： 013304.

［11］""nbsp;" XU M， WEI M J， YANG T， et al. An embedded boundary approach for the simulation of a flexible flapping wing at different density ratio［J］. European journal of mechanics - B， 2016， 55： 146-156.

［12］"""" CARRIóN M， STEIJL R， WOODGATE M， et al. Aeroelastic analysis of wind turbines using a tightly coupled"" CFD-CSD"" method［J］." Journal" of" fluids"" and structures， 2014， 50： 392-415.

［13］"""" MIAO W P， LI C， WANG Y B， et al. Study of adaptive blades in extreme environment using fluid-structure interaction method［J］. Journal of fluids and structures， 2019， 91： 102734.

［14］"""" 王蕤， 仲繼澤， 徐自力， 等. 動(dòng)網(wǎng)格區(qū)域?qū)θ~片顫振流固耦合計(jì)算效率及精度的影響［J］. 推進(jìn)技術(shù)， 2017， 38（9）： 2086-2092.

WANG R， ZHONG J Z， XU Z L， et al. Effects of coverage of dynamic mesh region on efficiency and accuracy of coupled fluid structure simulation for blade flutter［J］. Journal of propulsion technology， 2017， 38（9）： 2086-2092.

［15］"""" BAO M L， DING Y， SANG M S， et al. Modeling and evaluating nodal resilience of multi-energy systems under windstorms［J］. Applied energy， 2020， 270： 115136.

［16］"""" JONKMAN J， BUTTERFIELD S， MUSIAL W， et al. Definition of a 5-MW reference wind turbine for offshore system development［R］. National Renewable Energy Lab（NREL）， Golden， CO （United States）， 2009.

［17］"""" 劉樂璠， 鄭直， 閔為， 等. 基于動(dòng)態(tài)重疊網(wǎng)格的閥芯振蕩空化的研究［J］. 液壓與氣動(dòng)， 2023， 47（5）： 34-41.

LIU L F， ZHENG Z， MIN W， et al. Research on cavitation of valve core oscillation based on dynamic overset"" mesh［J］."" Chinese"" hydraulics"" amp;" pneumatics， 2023， 47（5）： 34-41.

［18］"""" RIYADH A A， ZISHAN A M， KHALIFA O F. Performance enhancement of a small-scale wind turbine featuring"" morphed" trailing" edge［J］." Sustainable" energy technologies and assessments， 2021， 46： 101229.

［19］"""" 胡超， 周丙浩， 馬勇， 等. 基于FAST和Simulink的海上風(fēng)力機(jī)葉片扭角和弦長優(yōu)化設(shè)計(jì)研究［J］. 太陽能學(xué)報(bào)， 2021， 42（3）： 135-141.

HU C， ZHOU B H， MA Y， et al. Optimization research of twist angle and chord length of offshore wind turbine blade based" on" fast"" and"" simulink［J］." Acta"" energiae"" solaris sinica， 2021， 42（3）： 135-141.

［20］"""" TANG D， BAO S Y， LUO L J， et al. A CFD/CSD coupled method with high order and its applications in flow induced vibrationsnbsp;" of"" tube" arrays" in" cross" flow［J］." Annals"" of nuclear energy， 2019， 130： 347-356.

STUDY ON DYNAMIC RESPONSE CHARACTERISTICS OF

WIND TURBINE ROTOR FLUID-STRUCTURE INTERACTION UNDER

TYPHOON WIND SPEED

Zhang Ruixing1，2，An Liqiang1，2，He Lun1，2，Zhang Ying3

（1. Department of Mechanical Engineering， North China Electric Power University， Baoding 071003， China；

2. Hebei Key Laboratory of Electric Machinery Health Maintenance amp; Failure Prevention， Baoding 071003， China；

3. College of Civil Engineering and Architecture， Hebei University， Baoding 071002， China）

Abstract：High-precision fluid-structure interaction simulation is an effective method for studying the structural reliability of wind turbines under typhoon conditions. However， the substantial computational costs and issues related to negative volume grids significantly affect the efficiency and stability of simulation. In this study， an improved fluid-structure interaction method was applied to investigate the load and response characteristics of wind turbines at different wind spe eds. Its accuracy was validated using NREL data and wind turbine blade modal vibration experiments. Compared to traditional fluid-structure coupling methods， this approach reduces the number of dynamic grids used， avoids negative volume grids， enhances computational stability， and shortens computation time by approximately 50%. Simulation results indicate that increasing the wind speed from 25 m/s to 35 m/s leads to a 6.8-fold increase in fluctuation amplitude for blade thrust and a 9.8-fold increase for torque loads. Additionally， structural buckling was observed on the suction side of the blade， which is close to the location of blade fractures observed in actual typhoon conditions. This may indicate a structural weak point of wind turbine blades under typhoon environments， offering valuable insights for subsequent wind turbine design optimization.

Keywords：wind turbine blades； typhoon environment； vibration analysis fluid structure interaction； numerical models； blade strength