席仁強，許成順，杜修力，許 坤

(1. 常州大學(xué)機械與軌道交通學(xué)院，江蘇，常州 213164；2. 北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點實驗室，北京 100124)

為應(yīng)對環(huán)境挑戰(zhàn)，風(fēng)能在現(xiàn)代能源體系中所占比重不斷擴大[1]。20 世紀90 年代，風(fēng)力發(fā)電率先在歐洲興起；進入21 世紀，隨著兆瓦級風(fēng)力機的應(yīng)用，風(fēng)力發(fā)電從北歐擴展至世界各地[2-3]。近年來，越來越多的風(fēng)電場在高地震危險性區(qū)域建設(shè)，風(fēng)力機抗震受到研究者的關(guān)注[4]。對于多兆瓦風(fēng)力機，葉片柔度大，當(dāng)葉尖位移峰值為1 m 時，其法向振動速度峰值可達6 m/s[5]。因此，葉輪氣動荷載計算需考慮氣彈性效應(yīng)[6]。根據(jù)氣彈性模型，風(fēng)力機結(jié)構(gòu)地震響應(yīng)有耦合、解耦兩類分析方法。

耦合方法直接考慮風(fēng)力機動力響應(yīng)的氣動-伺服-彈性耦合行為。Witcher[7]利用BLADED 軟件，開展了風(fēng)-地震聯(lián)合激勵下風(fēng)力機動力反應(yīng)分析，結(jié)果表明：該軟件能夠模擬氣彈性和控制器的影響。隨后，Prowell 等[8]對開源軟件FAST 進行二次開發(fā)，使其能夠?qū)嵤╋L(fēng)-地震作用下風(fēng)力機動力反應(yīng)分析，他們還通過振動臺試驗和數(shù)值模擬驗證該方法的可靠性[9]。改進的FAST 程序被廣泛用于分析風(fēng)力機地震反應(yīng)[10-12]，相關(guān)工作[13-14]表明風(fēng)荷載與地震作用效應(yīng)具有顯著的耦合行為。然而，現(xiàn)有風(fēng)力機氣彈性軟件的建模能力有限，難以模擬復(fù)雜基礎(chǔ)結(jié)構(gòu)和巖土體非線性行為；同時，由于風(fēng)速場、地震動具有隨機性，耦合方法計算量較大。這些困難推動了風(fēng)力機解耦動力響應(yīng)分析方法的發(fā)展。

氣動-彈性解耦地震響應(yīng)分析方法的應(yīng)用要早于耦合方法。針對千瓦級風(fēng)力機，葉輪-機艙(RNA)被簡化為塔頂處的集中質(zhì)量，氣動力則被當(dāng)做靜荷載[15]。根據(jù)葉素動量(BEM)理論，Zhao和Mai?er[16]將葉輪氣動荷載分解為平均推力和脈動推力，后者可等效為氣動阻尼。鑒于多兆瓦風(fēng)力機葉片結(jié)構(gòu)模型對系統(tǒng)的地震響應(yīng)至關(guān)重要[17]，Asareh 等[18]通過梁單元模擬葉片，將葉輪氣動阻尼力等效為塔架模態(tài)氣動阻尼比，提出了一種風(fēng)力機地震響應(yīng)解耦分析方法，并開展NREL 5 MW風(fēng)力機地震易損性分析。Santangelo 等[19]通過耦合、解耦方法計算結(jié)果的比較，評價了這種解耦方法的可靠性。解耦方法計算效率高，被廣泛用于風(fēng)力機地震響應(yīng)分析。然而，該解耦方法將不同模態(tài)和風(fēng)速的氣動阻尼比取為相同值。這一假定與水平軸風(fēng)力機氣動阻尼研究的最新研究結(jié)論并不一致[20]。同時，Xi 等[21]選取一組包括近場、遠場和近場脈沖地震動的強震記錄，發(fā)現(xiàn)該方法不具有普適性，可能產(chǎn)生較大誤差。

考慮到現(xiàn)有解耦方法的氣動阻尼模型存在錯誤、其仿真結(jié)果也不滿足工程精度要求，本文基于風(fēng)力機模態(tài)氣動阻尼理論模型，通過最小二乘法，簡化氣動阻尼比與平均風(fēng)速關(guān)系，形成一種改進的風(fēng)力機地震響應(yīng)解耦分析方法，通過耦合、解耦方法計算結(jié)果的比較，驗證該方法的可靠性，并評估解耦方法的計算效率。

1 氣動阻尼簡化模型

圖1(a)為典型的三葉片水平軸風(fēng)力機，是工業(yè)界應(yīng)用最廣泛的機型。圖示的慣性參考系xyz是固定于地面的笛卡爾坐標(biāo)系，x軸沿風(fēng)力機前后向，指向風(fēng)速場下風(fēng)向，z軸豎直向上，y軸沿風(fēng)力機側(cè)向。

1.1 模態(tài)氣動阻尼理論

根據(jù)文獻[22]建議，不考慮葉片變形，依據(jù)葉素動量理論，半徑r處葉素的氣動荷載為：

葉輪旋轉(zhuǎn)使得風(fēng)力機成為時變體系，為建立系統(tǒng)氣動阻尼模型，參照文獻[22]，將葉片近似為剛體，塔架變形引起的葉片振動速度可表示為：

這些氣動荷載增量與塔頂平動和轉(zhuǎn)動速度大小成正比，方向相反，對風(fēng)力機具有阻尼效應(yīng)，因此，被稱為氣動阻尼力。Xi 等[20]將式(9)～式(10)所示的葉片分布氣動阻力向塔頂截面形心處簡化，并根據(jù)風(fēng)力機前后向、側(cè)向模態(tài)，得到水平軸風(fēng)力機模態(tài)氣動阻尼比，具體的推導(dǎo)過程可參照該文獻。

NREL 5 MW 風(fēng)力機是應(yīng)用最廣泛的參考機型，主要技術(shù)參數(shù)如表1 所示，更多的結(jié)構(gòu)、控制系統(tǒng)參數(shù)可參考其研究報告[23]。本文采用歐拉梁單元模擬葉片和塔架動力行為，如圖1(b)所示，將塔架等分為50 個梁單元。葉片劃分為17 個梁單元，如圖1(c)所示，長度依次為：自葉根起，前3 個2.73 m，中間11 個4.1 m，其余3 個2.73 m。該風(fēng)力機前后向、側(cè)向前二階振型的模態(tài)氣動阻尼比如圖2 所示。對于前后向模態(tài)，風(fēng)速小于額定風(fēng)速時，氣動阻尼比單調(diào)遞增；風(fēng)速超過額定風(fēng)速后，氣動阻尼比變化較小，呈單調(diào)遞減趨勢。對于側(cè)向模態(tài)，氣動阻尼比單調(diào)遞增。同時，前后向和側(cè)向二階模態(tài)的氣動阻尼比顯著小于一階模態(tài)。因此，現(xiàn)有解耦方法將各模態(tài)的氣動阻尼比取為常數(shù)是不恰當(dāng)?shù)摹?/p>

表1 NREL 5 MW 風(fēng)機的主要性質(zhì)Table 1 Properties of NREL 5 MW wind turbines

圖2 NREL 5 MW 風(fēng)力機模態(tài)氣動阻尼比Fig. 2 Modal aerodynamic damping ratios of NREL 5 MW wind turbine

1.2 雙線性氣動阻尼模型

水平軸風(fēng)力機的模態(tài)氣動阻尼比與平均風(fēng)速有關(guān)，呈非線性變化，不利于其在工程中的應(yīng)用。因此，本節(jié)通過最小二乘法，對不同風(fēng)速條件下的風(fēng)力機模態(tài)氣動阻尼比進行簡化。

根據(jù)圖2，對于前后向模態(tài)，平均風(fēng)速介于切入風(fēng)速、額定風(fēng)速之間時，假定氣動阻尼比呈線性變化；平均風(fēng)速介于額定風(fēng)速、切出風(fēng)速之間時，假定氣動阻尼比為常數(shù)。因此，風(fēng)力機前后向一階、二階模態(tài)氣動阻尼比均可表示為：

表2 擬合曲線參數(shù)Table 2 Parameters of fitting curves

圖3 和圖4 分別為NREL 5 MW 風(fēng)力機前后向、側(cè)向模態(tài)氣動阻尼比。擬合值與理論解具有良好的一致性，因此，可采用這一簡化模型描述氣動阻尼比隨輪轂高度處10 min 平均風(fēng)速的變化規(guī)律。

圖3 NREL 5 MW 風(fēng)力機前后向模態(tài)氣動阻尼比Fig. 3 Modal aerodynamic damping ratios of fore-aft modes for NREL 5 MW wind turbine

圖4 NREL 5 MW 風(fēng)力機側(cè)向模態(tài)氣動阻尼比Fig. 4 Modal aerodynamic damping ratios of side-side modes for NREL 5 MW wind turbine

2 風(fēng)力機地震響應(yīng)的解耦分析方法

考慮氣動-彈性效應(yīng)，半徑r處葉素的法向氣動荷載為：

根據(jù)式(15)～式(19)，葉片振動速度影響葉輪氣動荷載，風(fēng)-地震激勵下的風(fēng)力機動力響應(yīng)屬于典型的氣彈性問題。現(xiàn)代風(fēng)力機廣泛采用了變速-變槳距控制技術(shù)[2]。因此，需要采用氣動-伺服-彈性耦合模型分析風(fēng)-地震共同激勵下的風(fēng)力機動力響應(yīng)。以FAST 為例，該軟件基于Kane 方法，建立的風(fēng)力機系統(tǒng)運動控制方程為：

3 算例與分析

研究表明：現(xiàn)有解耦方法的計算精度與地震動有關(guān)[21]，因此，參照FEMA 建議[25]，選取一組包括近場、遠場和脈沖地震動的強震記錄作為輸入地震動?？紤]到本文旨在探討氣彈性解耦地震響應(yīng)分析方法的可靠性，將這些原始記錄作為輸入地震動。

一個完整的強震記錄包含兩個水平分量和一個豎向分量。風(fēng)力機的氣動-彈性行為主要在其前后向，包括其他兩個分量會低估解耦模型的誤差。因此，對于所選強震記錄，采用附錄A 所列的水平分量為輸入地震動，并假定地震動和風(fēng)速方向一致，不考慮偏航，這是氣彈性效應(yīng)最為突出的工況?？偰M時長取600 s，時間步長取0.002 s，為避免初始條件的影響，地震動400 s 時開始輸入。驗證性算例表明計算參數(shù)的取值能夠保證數(shù)值模擬的精度與穩(wěn)定性。

3.1 穩(wěn)態(tài)風(fēng)-地震激勵

鑒于1.1 節(jié)所述氣動阻尼模型針對均勻穩(wěn)態(tài)風(fēng)，本節(jié)探討了穩(wěn)態(tài)風(fēng)-地震動激勵下新解耦方法分析風(fēng)力機地震響應(yīng)的可靠性。附表A1 所選的2 號和4 號地震動分量分別來自1999 年土耳其Duzce 地震的Bolu 記錄和臺灣Chi-Chi 地震的Chy101 記錄，各自的加速度時程和反應(yīng)譜如圖5所示。其中，無量綱加速度反應(yīng)譜值q通過式(23)獲得：

圖5 Bolu 強震記錄的一個水平分量Fig. 5 One horizontal component of Bolu seismic record

將風(fēng)速取為11.4 m/s，首先，以Bolu 波為輸入地震動。穩(wěn)態(tài)風(fēng)、地震動單獨作用引起的塔底彎矩時程如圖6 所示。穩(wěn)態(tài)風(fēng)單獨激勵下，由于氣彈性效應(yīng)，塔底彎矩均值為60 MN·m，塔頂加速度峰值為0.01 m/s2，因此，其對塔頂加速度的貢獻可以忽略。圖7 為耦合、新解耦方法得到的塔頂加速度和塔底彎矩時程。整體而言，兩種方法所得結(jié)果具有較好的一致性。耦合、新解耦方法得到的塔頂加速度峰值分別為2.21 m/s2和2.23 m/s2。對于塔底彎矩，耦合、新解耦方法得到的峰值分別為156 MN·m 和158 MN·m。相對誤差不超過5 %，滿足工程精度要求。

圖6 Bolu 強震記錄單獨激勵下的風(fēng)力發(fā)電機塔底彎矩Fig. 6 Tower-base bending moment of wind turbine excited by Bolu seismic record only

圖7 穩(wěn)態(tài)風(fēng)-Bolu 波共同激勵下的風(fēng)力發(fā)電機動力響應(yīng)Fig. 7 Dynamic response of wind turbine excited by combined steady wind and Bolu wave

隨后，以Chy101 波為輸入地震動，風(fēng)速仍取為11.4 m/s。地震動單獨作用引起的塔底彎矩時程如圖8 所示。圖9 為耦合、新解耦方法得到的塔頂加速度和塔底彎矩時程。耦合方法和新解耦方法得到的塔頂加速度峰值分別為3.14 m/s2和3.29 m/s2。對于塔底彎矩，耦合方法和解耦方法得到的峰值分別為173 MN·m 和185 MN·m。相對誤差小于10 %，滿足工程精度要求。

圖8 Chy101 強震記錄單獨激勵下的風(fēng)力發(fā)電機塔底彎矩Fig. 8 Tower-base bending moment of wind turbine excited by Chy101 seismic record only

圖9 穩(wěn)態(tài)風(fēng)-Chy101 波激勵下的風(fēng)力發(fā)電機動力響應(yīng)Fig. 9 Dynamic response of wind turbine excited by steady wind and Chy101 wave

考慮篇幅，本節(jié)僅列出上述兩個算例。實際上，對于附表A1 的所有輸入地震動，穩(wěn)態(tài)風(fēng)-地震激勵下，新解耦方法分析風(fēng)力機動力響應(yīng)具有足夠的精度，滿足工程要求。

3.2 湍流風(fēng)-地震激勵

風(fēng)具有顯著的湍流特性，脈動風(fēng)功率譜密度函數(shù)采用Kaimal 譜：

式中：Δr為兩節(jié)點的距離；b為衰減系數(shù)，取為12；LC為湍流積分尺度，取為340.2 m。將湍流風(fēng)速場總時長取為650 s，利用Turbsim 軟件，生成風(fēng)速場樣本。

根據(jù)已有研究，輪轂高度處平均風(fēng)速等于額定風(fēng)速時，風(fēng)力機受到的氣動載荷最大。首先，輸入地震動為Bolu 波，將輪轂高度處平均風(fēng)速取為11.4 m/s，采用上述模型，生成5 個隨機風(fēng)速場，分別采用耦合、新解耦和現(xiàn)有解耦方法分析風(fēng)力機地震響應(yīng)，塔架響應(yīng)峰值如圖10 所示。由于生成了5 個風(fēng)速場，此處的峰值是關(guān)于風(fēng)速場樣本的平均。對于塔架位移、加速度、剪力和彎矩峰值，新解耦方法與耦合方法的誤差不超過10%，滿足工程精度要求；而現(xiàn)有解耦方法具有顯著的誤差。由于脈動風(fēng)的影響，湍流風(fēng)速隨時間變化。現(xiàn)有解耦方法采用的模態(tài)阻尼和瑞利阻尼模型均為時不變的，也為了保持解耦方法的易實施性，因此，根據(jù)輪轂高度處10 min 平均風(fēng)速值確定模態(tài)氣動阻尼比。

圖10 Bolu 波-湍流風(fēng)激勵下的塔身響應(yīng)峰值Fig. 10 Response amplitudes of the tower when the wind turbine excited by Bolu wave and turbulent wind

為涵整個運行狀態(tài)，將輪轂高度處10 min 平均風(fēng)速分別取為5 m/s、7 m/s、9 m/s、11.4 m/s、13 m/s、15 m/s、18 m/s、21 m/s 和24 m/s，圖11為塔頂加速度和塔底彎矩幅值。新解耦模型與耦合模型的相對誤差不超過10%，滿足工程精度要求；現(xiàn)有解耦方法則存在較大誤差。

圖11 Bolu 波-風(fēng)激勵下風(fēng)力發(fā)電機響應(yīng)峰值Fig. 11 Response amplitude of wind turbines excited by Bolu seismic wave and wind

隨后，輸入地震動為Chy101 波，輪轂高度處平均風(fēng)速等于11.4 m/s 的塔架響應(yīng)峰值如圖12 所示。新解耦方法得到的塔架位移、加速度、剪力和彎矩峰值與基準方法的誤差不超過10%。圖13為不同平均風(fēng)速條件下的塔頂加速度和塔底彎矩峰值。新解耦模型與耦合模型的相對誤差不超過10%，滿足工程精度要求。對于Bolu 波和Chy101波，湍流風(fēng)-地震動激勵下，新解耦模型分析風(fēng)力機動力響應(yīng)的精度滿足工程需要，且顯著高于現(xiàn)有解耦方法。

圖12 Chy101 波-湍流風(fēng)激勵下的塔身響應(yīng)峰值Fig. 12 Response amplitude of the tower induced by Chy101 wave and turbulent wind

圖13 Chy101 波-風(fēng)激勵下的風(fēng)力發(fā)電機響應(yīng)峰值Fig. 13 Response amplitude of wind turbines excited by Chy101 seismic wave and wind

圖14 為葉片各單元法向氣動荷載，分布特征與輪轂高度處平均風(fēng)速有關(guān)。對于葉根起的前3 個單元，氣動荷載較小，因此，葉輪氣動合力受它們的影響較小。輪轂高度處10 min 平均風(fēng)速等于切入風(fēng)速(3 m/s)時，相對風(fēng)速最小，此時，地震引起的葉片振動速度平方與相對風(fēng)速平方之比η 最大：

圖14 葉片各單元法向氣動荷載Fig. 14 Aerodynamic forces of blade element in the normal direction

輸入地震動為Bolu 波、Chy101 波的葉片振動速度平方比如圖15 所示。對于前4 個節(jié)點，相對風(fēng)速較小，因此，速度平方比η 較大。此時，盡管式(22)會產(chǎn)生一定的截斷誤差，然而，這些節(jié)點所在單元對葉輪氣動荷載貢獻小，因此，對于支撐結(jié)構(gòu)動力響應(yīng)的影響較小。對于其他地震動，地震引起的葉片振動速度平方與相對風(fēng)速平方之比η 也具有這一特征。實際上，新解耦方法和耦合方法仿真結(jié)果的比較也間接表明了這一假定的合理性。

圖15 葉片振動速度平方比Fig. 15 Ratio for squared blade vibration velocity

3.3 解耦模型誤差分析

為評估解耦模型分析風(fēng)力機地震響應(yīng)的可靠性，支撐結(jié)構(gòu)響應(yīng)誤差δ 定義為：

輪轂高度處平均風(fēng)速等于11.4 m/s 時，塔頂加速度和塔底彎矩誤差如圖16 所示。對于所有地震動，風(fēng)力機支撐結(jié)構(gòu)響應(yīng)峰值誤差均小于15%，與Xi 等[21]評估的解耦方法相比，本文改進的解耦模型具有更高的精度。

圖16 平均風(fēng)速等于11.4 m/s 的解耦模型誤差Fig. 16 Errors of uncoupled model when mean wind speed is 11.4 m/s

輪轂高度處平均風(fēng)速等于5 m/s 和18 m/s 時，塔頂加速度和塔底彎矩誤差分別如圖17 和圖18所示。對于所有地震動，支撐結(jié)構(gòu)響應(yīng)誤差均小于15 %。實際上，平均風(fēng)速取其他值時，新解耦模型預(yù)測風(fēng)力機支撐結(jié)構(gòu)響應(yīng)峰值的誤差也不超過15 %，顯著小于Xi 等[20]評估的解耦方法。

圖17 平均風(fēng)速等于5 m/s 的解耦模型誤差Fig. 17 Errors of uncoupled model when mean wind speed is 5 m/s

圖18 平均風(fēng)速等于18 m/s 的解耦模型誤差Fig. 18 Errors of uncoupled model when mean wind speed is 18 m/s

3.4 計算效率的比較

GL 標(biāo)準[27]要求將地震作用與所有可能出現(xiàn)的環(huán)境荷載組合?？紤]到這些作用的隨機性，風(fēng)力機地震響應(yīng)分析涉及多種荷載工況。以3.3 節(jié)的湍流風(fēng)-地震動激勵下風(fēng)力機動力響應(yīng)分析為例，涉及風(fēng)速場樣本45 個、強震記錄50 個。耦合方法需要開展2 250 次時域地震響應(yīng)分析，解耦方法僅需進行95 次時域分析。當(dāng)然，解耦方法還需將風(fēng)作用效應(yīng)與地震作用效應(yīng)時程線性疊加。整體而言，解耦方法的計算量小于耦合方法。對于本算例，使用的計算機置為：CPU Intel I7 7700k，內(nèi)存16 GB。耦合方法消耗的CPU 時長為11 250 min，解耦方法消耗的時長為670 min。因此，解耦方法計算效率高于耦合方法。

4 結(jié)論與展望

本文基于水平軸風(fēng)力機氣動阻尼理論，通過最小二乘法將氣動阻尼模型簡化，改進風(fēng)力機地震響應(yīng)的解耦分析方法；隨后，以NREL 5 MW風(fēng)力機為例，探討穩(wěn)態(tài)風(fēng)、湍流風(fēng)分別與地震聯(lián)合激勵下風(fēng)力機地震響應(yīng)?？傻靡韵陆Y(jié)論：

(1) 線性化模態(tài)氣動阻尼模型與風(fēng)力機氣動阻尼理論解具有良好的一致性，因此，可采用這一簡化模型描述水平軸風(fēng)力機模態(tài)氣動阻尼比隨風(fēng)力機輪轂高度處10 min 平均風(fēng)速的變化規(guī)律。

(2) 對于本文選取的所有輸入地震動，解耦方法分析風(fēng)-地震聯(lián)合激勵下風(fēng)力機動力響應(yīng)時，塔頂加速度和塔底彎矩峰值誤差均不超過15 %，基本滿足工程精度需要。

(3) 對于本文3.3 節(jié)的驗證性算例，風(fēng)力機耦合地震響應(yīng)分析方法需開展的時域模擬數(shù)量超過解耦方法的20 倍，由于解耦方法需要在時域?qū)L(fēng)、地震作用效應(yīng)疊加，耦合方法消耗的CPU 時長約解耦方法的16 倍。

為保證解耦、耦合兩類模型的差異僅在于氣動-彈性處理方法，本文的數(shù)值模擬工作均在FAST 軟件中實施。實際上，這種解耦方法可用于商用有限元軟件。需要指出，本文將風(fēng)力機氣動阻尼效應(yīng)以模態(tài)氣動阻尼施加于塔架，還可以通過塔頂集中阻尼器或者葉片分布阻尼器的方式實現(xiàn)，這些方法待進一步探討。同時，NREL 發(fā)布的FAST 軟件尚不能考慮土-結(jié)相互作用，后續(xù)工作應(yīng)基于FAST 程序開發(fā)樁-土非線性相互作用模型，進一步驗證上述解耦方法在考慮樁-土非線性相互作用后的適用性。

附錄：

附表A1 地震地面運動Attached table A1 Earthquake Ground motions