考慮葉片旋轉與土構耦合作用的海上風電機組 塔筒動態特性研究

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

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

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

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

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

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

1 風電機組有限元模型

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

NREL 5 MW單樁式海上風電機組整體結構模型如圖1所示，風電機組主要參數[14]如表1所示。

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

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

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

表2 風電機組結構材料屬性

1.2 葉片旋轉工況的施加

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

σ=σr+Hε

(1)

式中：σ為葉片應力，MPa；σr為初始應力，MPa；H為應力-應變關系方程；ε為葉片應變張量。

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

(2)

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

由葉片旋轉引起的離心剛度矩陣K*為：

(3)

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

1.3 土-結構耦合作用模型

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

表3 海上風電場土質實測參數

樁土接觸底部采用硬接觸形式，側向采用摩爾庫倫摩擦罰函數形式，通過摩擦因數定義接觸面關系，表達式為：

τc=μp

(4)

式中：τc為接觸面滑移臨界切應力，MPa；μ為摩擦因數；p為接觸壓力，MPa。

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

圖2 單樁土體模型內部截面示意圖

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

1.4 單樁式海上風電機組結構運動數學模型

單樁式海上風電機組結構運動微分方程[17]為：

(5)

附加水質量和離心剛度矩陣后，風電機組結構運動微分方程變為：

(6)

式中：M0為考慮附加水質量的結構質量矩陣；K0為結構整體剛度矩陣。

M0可以表示為：

M0=M+M*

(7)

式中：M*為附加質量矩陣。

C采用Rayleigh阻尼：

C=αM+βK

(8)

(9)

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

結構整體剛度矩陣K0為：

K0=K+K*

(10)

(11)

式中：L為線性應變-位移關系矩陣；ψ為振動幅值矩陣。

1.5 網格劃分

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

表4 網格劃分無關性驗證結果

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

2 風力機結構與載荷計算模型

2.1 脈動風載荷

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

(12)

(13)

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

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

圖3 平均風速為11.4 m/s時100 s內的脈動風速分布

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

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

(14)

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

風剪切采用指數模型：

(15)

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

風輪水平軸向推力Ft為：

(16)

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

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

(17)

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

葉輪旋轉力矩Mr為：

Mr=9 950Pγ/ω

(18)

式中：P為風電機組功率，W；γ為發電效率；ω為葉片轉速，r/min。

葉輪俯仰力矩Mp為：

(19)

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

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

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

2.2 海浪載荷

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

(20)

式中：φ為速度勢；Hw為海浪高度，m；ζ為波數；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為慣性力系數；ρw為水密度，kg/m3；D為結構直徑，m；Cd為拖曳力系數。

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

(27)

(28)

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

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

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

3 塔筒振動特性分析

3.1 不考慮土構耦合時塔筒振動特性分析

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

3.1.1 風電機組結構的模態特性

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

表5 葉片旋轉工況下塔筒結構模態頻率

3.1.2 風電機組結構的位移響應特性

通過模擬得到的風電機組結構瞬態響應位移云圖如圖7所示。風電機組塔筒頂部前后方向與側向位移時程曲線如圖8所示。

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

(a) 機組總位移

(a) 前后方向

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

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

圖9給出了葉片旋轉工況下塔筒頂部前后方向與側向加速度響應的功率譜密度。

(a) 前后方向

由圖9可知，在0.247 6 Hz處出現明顯峰值，這對應了結構一階模態頻率，說明主要能量集中在前兩階模態中。從圖9可以很明顯地看出能量在側向上小得多，這導致在該方向上較小的塔筒位移響應。

3.1.3 風電機組結構應力分布特性

通過模擬得到的風電機組結構瞬態響應等效應力云圖和最大等效應力時程曲線如圖10所示。

(a) 等效應力云圖

在此工況下，塔筒瞬態響應最大等效應力分布在塔筒與機艙連接處以及塔筒底部，最大等效應力出現在1.6 s，其值為703.65 MPa。

3.2 土構耦合作用下塔筒振動特性分析

3.2.1 風電機組結構的模態特性

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

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

表6 考慮土構耦合時風電機組塔筒結構模態頻率

3.2.2 風電機組結構的位移響應特性

通過模擬得到的葉片旋轉與土構耦合共同作用下風電機組結構瞬態響應位移云圖如圖12所示，塔筒頂部前后方向與側向位移時程曲線如圖13所示。

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

(a) 前后方向

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

圖14給出了葉片旋轉與土構耦合共同作用下塔筒頂部前后方向與側向加速度功率譜密度。在0.157 Hz處塔筒頂部前后方向與側向加速度功率譜密度出現明顯的峰值，這對應了土構耦合作用下結構一階模態頻率。將圖14與圖9進行比較發現，考慮土構耦合時塔筒頂部加速度功率譜密度明顯增大。

(a) 前后方向

3.2.3 風電機組結構應力分布特性

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

(a) 等效應力云圖

3.3 風電機組結構共振分析

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

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

(a) 葉片旋轉作用

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

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

4 結 論

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

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

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

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