考慮葉片旋轉(zhuǎn)與土構(gòu)耦合作用的海上風電機組 塔筒動態(tài)特性研究

樊 昂， 李錄平， 歐陽敏南， 陳尚年

(長沙理工大學 能源與動力工程學院，長沙 410114)

目前，單樁式海上風電機組應(yīng)用越來越廣泛[1]。柔性風電機組由于較大的尺寸更容易受到外部振動源的影響，這可能會影響風電機組功率的輸出，造成塔筒的疲勞損傷，甚至在極端條件下直接引發(fā)倒塔事故。為了保證海上風電機組的安全、可靠運行，準確了解風力機塔筒在外部振動載荷作用下的動態(tài)特性是非常重要的。對此，研究人員進行了大量的研究工作，以獲得塔筒的動力學特性。

為了簡化分析，通常假定風電機組處于停車狀態(tài)，忽略葉片的幾何構(gòu)型以及塔筒與葉片之間的相互作用[2]，將葉片簡化為位于塔頂?shù)募匈|(zhì)量[3]或簡化為兩自由度系統(tǒng)。而在工程實際中，葉片的幾何特性和轉(zhuǎn)速直接影響作用在葉片上的風荷載[4]。此外,在運行條件下，葉片產(chǎn)生的離心力會改變?nèi)~片的剛度和模態(tài)頻率[5],這反過來又會間接地影響塔筒的動態(tài)響應(yīng)。因此，簡化的集中質(zhì)量模型與停機狀況假設(shè)可能會導致結(jié)構(gòu)響應(yīng)估計不準確。而對于動力剛化效應(yīng)的研究大多局限于葉片及風輪本身，并沒有將其施加在塔筒動態(tài)特性分析中[6]。

另外，對于單樁這樣一個細長的柔性基礎(chǔ)，其與周圍土體的相互作用是不可避免的，土構(gòu)耦合(Soil-Structure Interaction，SSI)可以改變海上風電機組塔筒的振動特性和動力響應(yīng)[7]，會降低塔筒的模態(tài)頻率[8]。Bordón等[9]采用數(shù)值模擬和實驗方法研究了土構(gòu)耦合對風電機組振動特性的影響。然而在這些研究中，要么假設(shè)風電機組處于停車狀態(tài)，風輪質(zhì)量集中在塔筒頂部[10]，要么只注重于單樁基礎(chǔ)性能的分析[11-12]，而對單樁上方的塔筒振動特性研究則相對較少。在實際工程中，海上風電機組的塔筒結(jié)構(gòu)也很容易遭受破壞作用[13]，因此有必要研究塔筒在土構(gòu)耦合作用下的振動特性。

筆者采用有限元軟件Ansys對NREL 5 MW單樁式海上風電機組進行整體建模，模擬分析塔筒在葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下的動力響應(yīng)。在建立葉片與土體三維實體模型的基礎(chǔ)上，將葉片旋轉(zhuǎn)時產(chǎn)生的離心力與離心剛度作為初始狀態(tài)；在單樁周圍建立分層土體模型，模擬土壤與塔筒單樁的相互作用，定量研究葉片旋轉(zhuǎn)工況與土構(gòu)耦合作用對塔筒動力學行為的影響。

1 風電機組有限元模型

1.1 NREL 5 MW單樁式海上風電機組模型

NREL 5 MW單樁式海上風電機組整體結(jié)構(gòu)模型如圖1所示，風電機組主要參數(shù)[14]如表1所示。

圖1 NREL 5 MW單樁式海上風電機組三維實體模型

表1 NREL 5 MW單樁式海上風電機組主要參數(shù)

該風電機組的葉片采用聚酯材料，塔和單樁材料為鋼。筆者不考慮海水振動引起的附加質(zhì)量。風電機組結(jié)構(gòu)材料屬性見表2。

表2 風電機組結(jié)構(gòu)材料屬性

1.2 葉片旋轉(zhuǎn)工況的施加

考慮到葉片旋轉(zhuǎn)對塔筒振動特性的影響，將風電機組額定轉(zhuǎn)速(12.1 r/min)工況下葉片旋轉(zhuǎn)產(chǎn)生的離心力作為預應(yīng)力形成離心剛度矩陣，然后將其疊加到模態(tài)和瞬態(tài)剛度矩陣上，以此來考慮葉片旋轉(zhuǎn)產(chǎn)生的動力剛化效應(yīng)。動力剛化效應(yīng)下葉片的應(yīng)力-應(yīng)變關(guān)系[15]為：

σ=σr+Hε

(1)

式中：σ為葉片應(yīng)力，MPa；σr為初始應(yīng)力，MPa；H為應(yīng)力-應(yīng)變關(guān)系方程；ε為葉片應(yīng)變張量。

葉片在旋轉(zhuǎn)過程中產(chǎn)生的離心力總和Nr(i)為：

(2)

式中：r為葉片i點到轉(zhuǎn)軸中心的距離，m；m(r)為單位長度質(zhì)量，kg/m；Ωr為葉片旋轉(zhuǎn)速度，r/min；Δr為單位長度;R為葉輪半徑,m。

由葉片旋轉(zhuǎn)引起的離心剛度矩陣K*為：

(3)

式中：Li為i單元長度，m；Ni為i單元離心力，N。

1.3 土-結(jié)構(gòu)耦合作用模型

選擇風電機組單樁土體模型的直徑為80 m，高60 m，模型上下共分為3層：上層厚度為15 m，土質(zhì)為淤泥質(zhì)粉質(zhì)黏土；中層厚度為15 m，土質(zhì)為粉砂；底層厚度為30 m，土質(zhì)為粉質(zhì)黏土。土質(zhì)參數(shù)采用江蘇某海上風電場實測參數(shù)[16]，見表3。

表3 海上風電場土質(zhì)實測參數(shù)

樁土接觸底部采用硬接觸形式，側(cè)向采用摩爾庫倫摩擦罰函數(shù)形式，通過摩擦因數(shù)定義接觸面關(guān)系，表達式為：

τc=μp

(4)

式中：τc為接觸面滑移臨界切應(yīng)力，MPa；μ為摩擦因數(shù)；p為接觸壓力，MPa。

土體邊界條件為：底邊三向固定約束，外側(cè)徑向位移約束。模型中將單樁內(nèi)部填充物視為混凝土。單樁土體模型內(nèi)部截面示意圖如圖2所示。

圖2 單樁土體模型內(nèi)部截面示意圖

為減小計算誤差，在進行動態(tài)分析之前對單樁土體模型進行初始應(yīng)力平衡。先向單樁土體模型施加重力載荷，再將得到的土體內(nèi)部應(yīng)力作為初始條件寫入模型。

1.4 單樁式海上風電機組結(jié)構(gòu)運動數(shù)學模型

單樁式海上風電機組結(jié)構(gòu)運動微分方程[17]為：

(5)

附加水質(zhì)量和離心剛度矩陣后，風電機組結(jié)構(gòu)運動微分方程變?yōu)椋?/p>

(6)

式中：M0為考慮附加水質(zhì)量的結(jié)構(gòu)質(zhì)量矩陣；K0為結(jié)構(gòu)整體剛度矩陣。

M0可以表示為：

M0=M+M*

(7)

式中：M*為附加質(zhì)量矩陣。

C采用Rayleigh阻尼：

C=αM+βK

(8)

(9)

式中：ωi、ωj分別為第i階和第j階振型的模態(tài)頻率，Hz；ξi、ξj分別為第i階和第j階振型阻尼比。

結(jié)構(gòu)整體剛度矩陣K0為：

K0=K+K*

(10)

(11)

式中：L為線性應(yīng)變-位移關(guān)系矩陣；ψ為振動幅值矩陣。

1.5 網(wǎng)格劃分

采用Ansys Workbench自帶網(wǎng)格劃分模塊進行網(wǎng)格劃分，模型單元類型為實體單元，選擇占用計算資源較少的四面體網(wǎng)格與質(zhì)量較高的六面體網(wǎng)格相結(jié)合，在單樁與土體接觸處細化網(wǎng)格[18]。由于單樁土體模型結(jié)構(gòu)簡單、尺寸較大，且筆者著重分析塔筒性能，因此設(shè)置單樁土體模型網(wǎng)格尺寸大于風電機組結(jié)構(gòu)網(wǎng)格尺寸。同時進行網(wǎng)格劃分無關(guān)性驗證，以葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下塔筒瞬態(tài)響應(yīng)前后方向位移最大值為算例，無關(guān)性驗證結(jié)果如表4所示。

表4 網(wǎng)格劃分無關(guān)性驗證結(jié)果

由表4可知，方案三計算量合理，網(wǎng)格質(zhì)量滿足要求，算例結(jié)果與方案四差別不大，為合理利用計算資源，保證計算效率，故選擇方案三作為本文網(wǎng)格劃分方案。

2 風力機結(jié)構(gòu)與載荷計算模型

2.1 脈動風載荷

湍流風的風速時程由平均風速和脈動風速組成，選取風電機組的額定風速(11.4 m/s)作為平均風速，脈動風功率譜選取IEC 61400-3標準[19]中的Kaimal模型，風速譜Sk(f)的表達式為：

(12)

(13)

式中：f為風速頻率，Hz；σk為風速均方根值，m/s；Lk為整體尺度參數(shù)；uhub為輪轂高度處平均風速，m/s；Lu為湍流尺度參數(shù)，Lu=0.7min(60，hub)，其中hub為輪轂高度，m;k取a、b、c3個方向。

平均風速為11.4 m/s時100 s內(nèi)的脈動風速分布如圖3所示。

圖3 平均風速為11.4 m/s時100 s內(nèi)的脈動風速分布

作用在風電機組上的風載荷值按式(14)～式(19)確定。

作用在塔筒的風載荷Fz[20]為：

(14)

式中：ρ為空氣密度，kg/m3；At為塔筒受風面積，m2；V為風速，m/s；Ce為阻力系數(shù)。

風剪切采用指數(shù)模型：

(15)

式中：Vh為距地面高度h處風速，m/s；α為切變系數(shù)。

風輪水平軸向推力Ft為：

(16)

式中：Ab為掠葉面積，m2；Ct為推力系數(shù)。

葉片氣動載荷切向力Fy為：

(17)

式中：CP為風電機組風能利用系數(shù)；A為葉片切向面積,m2。

葉輪旋轉(zhuǎn)力矩Mr為：

Mr=9 950Pγ/ω

(18)

式中：P為風電機組功率，W；γ為發(fā)電效率；ω為葉片轉(zhuǎn)速，r/min。

葉輪俯仰力矩Mp為：

(19)

式中：B為葉片數(shù)；V1、V2分別為葉輪中心上、下2/3葉片長度處風速，m/s。

圖4為塔筒頂部脈動風載荷100 s時程曲線。

圖4 塔筒頂部脈動風載荷100 s時程曲線

2.2 海浪載荷

海浪會對單樁產(chǎn)生顯著的水平方向作用力，且海浪的高度是不規(guī)則的，可以從各個方向影響風電機組結(jié)構(gòu)。采用線性波浪理論描述海浪運動[21]：

(20)

式中：φ為速度勢；Hw為海浪高度，m；ζ為波數(shù)；T為波浪周期，s；x為水平坐標；z為豎向坐標；d為水深，m；ωw為波浪圓頻率，rad/s。

海浪高程η為：

(21)

(22)

(23)

采用Morison方程[22]計算海浪載荷Fwave：

Fwave=Fm+Fd

(24)

(25)

(26)

式中：Fm為慣性力，kN；Fd為拖曳力，kN；Cm為慣性力系數(shù)；ρw為水密度，kg/m3；D為結(jié)構(gòu)直徑，m；Cd為拖曳力系數(shù)。

單樁式海上風電機組是典型的高聳柔性結(jié)構(gòu)，其自身的變形與位移不可忽略，海浪中水質(zhì)點會受到影響，其速度和加速度會發(fā)生改變，因此要對式(25)和式(26)進行修正：

(27)

(28)

海浪載荷100 s時程曲線如圖5。

圖5 海浪載荷100 s時程曲線

下面的模擬計算中，在2個計算工況下均施加了海浪載荷，施加位置在單樁高于海床之上的部分。

3 塔筒振動特性分析

3.1 不考慮土構(gòu)耦合時塔筒振動特性分析

不考慮土構(gòu)耦合時，風電機組模型不包含單樁土體模型，塔筒底部與單樁頂部綁定接觸，單樁底部固定約束。將風電機組受載運行狀況下的葉片應(yīng)力矩陣與剛度矩陣施加在模型上，以此作為初始狀態(tài)。

3.1.1 風電機組結(jié)構(gòu)的模態(tài)特性

通過有限元計算，得到風電機組塔筒結(jié)構(gòu)前六階模態(tài)振型和模態(tài)頻率，分別見圖6和表5。由表5可知，葉片旋轉(zhuǎn)時所產(chǎn)生的離心力和離心剛度會增大塔筒的模態(tài)頻率，其中對前兩階模態(tài)影響最大，模態(tài)頻率變化率分別為7.52%和7.10%。

表5 葉片旋轉(zhuǎn)工況下塔筒結(jié)構(gòu)模態(tài)頻率

3.1.2 風電機組結(jié)構(gòu)的位移響應(yīng)特性

通過模擬得到的風電機組結(jié)構(gòu)瞬態(tài)響應(yīng)位移云圖如圖7所示。風電機組塔筒頂部前后方向與側(cè)向位移時程曲線如圖8所示。

由圖7可知，沿塔筒高度的瞬態(tài)響應(yīng)位移不同，最大位移響應(yīng)出現(xiàn)在塔頂，因此僅討論塔筒頂部位移。由圖8可知，在葉片旋轉(zhuǎn)工況下，塔筒頂部前后方向的位移主要在250～700 mm范圍內(nèi)波動，最大位移為846 mm。

(a) 機組總位移

(a) 前后方向

塔筒頂部側(cè)向位移主要在75～250 mm范圍內(nèi)波動，最大位移為288.78 mm；側(cè)向位移波動幅度較小，遠不如前后方向位移波動幅度大。

計算結(jié)果表明，塔筒的側(cè)向位移響應(yīng)明顯小于前后方向位移響應(yīng)，出現(xiàn)這種現(xiàn)象的主要原因為：第一是數(shù)值模擬中，脈動風載荷只沿前后方向作用在塔筒上；第二是因為當葉片處于旋轉(zhuǎn)狀態(tài)時，葉片揮舞方向(對應(yīng)于塔筒前后方向)上的風載荷遠大于擺振方向(對應(yīng)于塔筒側(cè)向)上的風載荷，葉片上較大的風載荷導致塔架與葉片之間更劇烈的相互作用。

圖9給出了葉片旋轉(zhuǎn)工況下塔筒頂部前后方向與側(cè)向加速度響應(yīng)的功率譜密度。

(a) 前后方向

由圖9可知，在0.247 6 Hz處出現(xiàn)明顯峰值，這對應(yīng)了結(jié)構(gòu)一階模態(tài)頻率，說明主要能量集中在前兩階模態(tài)中。從圖9可以很明顯地看出能量在側(cè)向上小得多，這導致在該方向上較小的塔筒位移響應(yīng)。

3.1.3 風電機組結(jié)構(gòu)應(yīng)力分布特性

通過模擬得到的風電機組結(jié)構(gòu)瞬態(tài)響應(yīng)等效應(yīng)力云圖和最大等效應(yīng)力時程曲線如圖10所示。

(a) 等效應(yīng)力云圖

在此工況下，塔筒瞬態(tài)響應(yīng)最大等效應(yīng)力分布在塔筒與機艙連接處以及塔筒底部，最大等效應(yīng)力出現(xiàn)在1.6 s，其值為703.65 MPa。

3.2 土構(gòu)耦合作用下塔筒振動特性分析

3.2.1 風電機組結(jié)構(gòu)的模態(tài)特性

通過有限元計算，得到葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下風電機組前六階模態(tài)振型和模態(tài)頻率，分別見圖11和表6。

由表6可知，土構(gòu)耦合作用會顯著降低風電機組塔筒結(jié)構(gòu)前兩階模態(tài)頻率，前兩階模態(tài)頻率的相對變化量分別為-31.68%和-31.74%。土構(gòu)耦合對塔筒結(jié)構(gòu)三階、四階模態(tài)頻率的影響略小，其相對變化量分別為-4.75%和-4.41%。這是因為當考慮土構(gòu)耦合時，單樁插入海床并被土壤包圍，這使得塔筒比底部完全固定的情況下更靈活，土壤的彈性使整個結(jié)構(gòu)系統(tǒng)的剛度降低，同時又增加了塔筒的參振質(zhì)量，因此塔筒的模態(tài)頻率降低。

表6 考慮土構(gòu)耦合時風電機組塔筒結(jié)構(gòu)模態(tài)頻率

3.2.2 風電機組結(jié)構(gòu)的位移響應(yīng)特性

通過模擬得到的葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下風電機組結(jié)構(gòu)瞬態(tài)響應(yīng)位移云圖如圖12所示，塔筒頂部前后方向與側(cè)向位移時程曲線如圖13所示。

由圖12和圖13可知，在葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下，塔筒前后方向位移瞬態(tài)響應(yīng)量明顯增大，響應(yīng)幅值在800～1 500 mm內(nèi)波動；在41.5 s時，塔筒頂部前后方向最大位移為1 581.52 mm，與不考慮土構(gòu)耦合時的最大位移846 mm相比，增幅為86.94%。這是因為考慮土構(gòu)耦合后，土體結(jié)構(gòu)的剛度、阻尼以及單樁周圍土體的變形將會對風電機組結(jié)構(gòu)的響應(yīng)產(chǎn)生較大影響，使其位移和應(yīng)力大大超過剛性基礎(chǔ)下的結(jié)果，并且塔筒是高柔結(jié)構(gòu)，該現(xiàn)象更為明顯。

(a) 前后方向

與不考慮土構(gòu)耦合時相比，考慮土構(gòu)耦合后塔筒頂部側(cè)向位移也有所增大,其響應(yīng)幅值在175～375 mm內(nèi)波動；塔筒頂部側(cè)向最大位移出現(xiàn)在37.5 s，達到423.33 mm，與不考慮土構(gòu)耦合時的最大位移288.78 mm相比，增幅為46.59%。塔筒頂部側(cè)向位移波動幅度遠小于其前后方向位移波動幅度。

圖14給出了葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下塔筒頂部前后方向與側(cè)向加速度功率譜密度。在0.157 Hz處塔筒頂部前后方向與側(cè)向加速度功率譜密度出現(xiàn)明顯的峰值，這對應(yīng)了土構(gòu)耦合作用下結(jié)構(gòu)一階模態(tài)頻率。將圖14與圖9進行比較發(fā)現(xiàn)，考慮土構(gòu)耦合時塔筒頂部加速度功率譜密度明顯增大。

(a) 前后方向

3.2.3 風電機組結(jié)構(gòu)應(yīng)力分布特性

通過模擬計算，得到風電機組結(jié)構(gòu)在葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下瞬態(tài)響應(yīng)等效應(yīng)力云圖和最大等效應(yīng)力時程曲線，如圖15所示。由圖15(b)可知，考慮土構(gòu)耦合時塔筒瞬態(tài)響應(yīng)最大等效應(yīng)力出現(xiàn)在1.6 s，其值為1 107.4 MPa，相比于不考慮土構(gòu)耦合時的最大等效應(yīng)力703.65 MPa，增幅為57.38%。這是因為土體結(jié)構(gòu)阻尼較小，對風電機組結(jié)構(gòu)應(yīng)力響應(yīng)影響較大，且單樁周圍的土體結(jié)構(gòu)在擠壓下產(chǎn)生位移，導致上部風電機組結(jié)構(gòu)產(chǎn)生應(yīng)力重分布。同時最大應(yīng)力分布位置主要出現(xiàn)在塔筒底部和單樁處，與實際情況相符合。

(a) 等效應(yīng)力云圖

3.3 風電機組結(jié)構(gòu)共振分析

為了避免風電機組結(jié)構(gòu)在復雜的外部激勵作用下產(chǎn)生共振，利用坎貝爾圖來甄別風電機組結(jié)構(gòu)潛在的共振點，并通過風電機組結(jié)構(gòu)的模態(tài)頻率與葉片激勵頻率的重合點來識別。其中，風電機組單葉片旋轉(zhuǎn)一周的激勵頻率稱為1P，三葉片旋轉(zhuǎn)一周的激勵頻率稱為3P。當風電機組運行時，坎貝爾圖中線條產(chǎn)生的任何交點即為共振點。

圖16(a)為只考慮葉片旋轉(zhuǎn)作用時的風電機組坎貝爾圖，圖16(b)為考慮葉片旋轉(zhuǎn)和土構(gòu)耦合共同作用時的風電機組坎貝爾圖。

(a) 葉片旋轉(zhuǎn)作用

由圖16(a)可知，在考慮葉片動力剛化效應(yīng)下，葉片以不同轉(zhuǎn)速旋轉(zhuǎn)時，在1P轉(zhuǎn)頻范圍內(nèi)，塔筒的一階及二階模態(tài)頻率會在轉(zhuǎn)速14.8 r/min附近與1P轉(zhuǎn)頻產(chǎn)生交點，發(fā)生共振。而三階及四階模態(tài)頻率分別在轉(zhuǎn)速14.5 r/min附近和15.8 r/min附近與3P轉(zhuǎn)頻產(chǎn)生交點，說明風電機組結(jié)構(gòu)在該3P轉(zhuǎn)頻內(nèi)存在2個共振點，15 r/min附近是風電機組易共振區(qū)域。風電機組以額定轉(zhuǎn)速(12.1 r/min)運行時，在1P、3P轉(zhuǎn)頻激勵下不會產(chǎn)生共振。

由圖16(b)可知，在葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下，塔筒的一階及二階模態(tài)頻率會在轉(zhuǎn)速9.5 r/min附近與1P轉(zhuǎn)頻產(chǎn)生交點，發(fā)生共振，該點位于在風電機組達到額定轉(zhuǎn)速之前，因此應(yīng)在運行到此轉(zhuǎn)速處快速通過，避免共振。而三階及四階模態(tài)頻率分別在轉(zhuǎn)速14 r/min附近和14.8 r/min附近與3P轉(zhuǎn)頻產(chǎn)生交點，該點與圖16(a)差別不大，都位于額定轉(zhuǎn)速外。

4 結(jié) 論

(1) 考慮土構(gòu)耦合時，風電機組塔筒結(jié)構(gòu)的前兩階模態(tài)頻率顯著降低，相對變化率分別為-31.68%和-31.74%,而對三階～六階模態(tài)頻率的影響較小。

(2) 當風電機組以額定轉(zhuǎn)速(12.1 r/min)旋轉(zhuǎn)時，在葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下，塔筒頂部在前后方向和側(cè)向上的最大位移比底部固定約束時大86.94%和46.59%。若忽略塔筒底部的土構(gòu)耦合作用，可能導致非保守的結(jié)構(gòu)響應(yīng)估計和結(jié)構(gòu)部件的不安全設(shè)計。

(3) 在葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下，塔筒的一階及二階模態(tài)頻率會在轉(zhuǎn)速9.5 r/min附近與1P轉(zhuǎn)頻產(chǎn)生共振點，在額定轉(zhuǎn)速以內(nèi)不會與3P轉(zhuǎn)頻產(chǎn)生共振。

(4) 在葉片旋轉(zhuǎn)與土構(gòu)耦合共同作用下，塔筒等效應(yīng)力顯著增大，最大等效應(yīng)力相對于不考慮土構(gòu)耦合時增大了57.38%，主要分布位置為塔筒底部和單樁處。