不同動力學建模方法對浮式風機系統動力響應的影響研究

楊禮東，姜玉挺，李 昊，李 焱，劉利琴

(1. 中國電力工程顧問集團東北電力設計院有限公司, 吉林 長春 130021；2. 天津大學 水利仿真與安全國家重點實驗室, 天津 300072)

0 引 言

風力發電機是將風能轉化為電能的大型工程設施，海上浮式風力機是風力發電的一種重要形式。由于海洋環境相較陸地更為復雜，與固定式風機相比，浮式風機會發生大范圍運動，從而影響風機葉片氣動載荷。因此采用耦合的分析方法對整個浮式風力機系統進行動力分析很有必要。

目前有許多有關浮式風力機系統的建模及分析方法，主要包括剛體動力學方法和剛-柔耦合的多體動力學方法。李焱等[1]基于剛體動力學理論，考慮氣動力與浮式基礎運動的耦合，建立了Spar 型風機系統動力學模型，并分析了浮式風機系統動力特性。Minu 等[2]與Shen 等[3]分別使用渦格法與勢流理論建立剛體模型計算了風機的非定常空氣動力學響應。Jonkman 等[4]采用氣動-水動-伺服-彈性全耦合方法初步分析了5 MW風機動力響應特性。Roberson 等[5]在此基礎上應用該仿真程序對幾種不同形式的浮式風機進行分析，對比研究了不同環境載荷下系統的動力響應。肖昌水等[6]基于Jourdain 原理和有限元離散方法研究了風機系統的動力響應。葉江舟等[7]研究了剛-柔耦合理論模型及其數值原理。陳嘉豪等基于OC3 平臺建立了時域耦合分析程序并加以驗證[8]，并對比研究了半潛式風機的氣動阻尼響應特性[9]。不同的建模方法具有不同的精度、速度與穩定性，單剛體模型的計算速度快，但計算精度相對較低。相對而言，剛-柔耦合模型的計算精度高，同時其計算成本也很高。如何權衡并合理使用2 種不同的方法成為亟待解決的問題。

基于此，本文采用單剛體和剛-柔耦合[6]2 種不同建模方法，在保證氣動力、水動力、系纜等外載荷相同算法的前提下，研究浮式分機系統動力響應，對比分析浮式風機系統剛體模型和剛-柔耦合模型對計算產生的影響。

1 浮式風機系統動力學方程

本文以5 MW 浮式風機系統為例進行分析，其中，浮式基礎為OC4 基礎，上部為NERL-5 MW 風力機，具體參數見文獻[10]。

1.1 單剛體動力學方程

將浮式風機系統處理為一個剛體，建立系統動力學方程。采用笛卡爾-右手坐標系，取其整體重心為慣性系原點，z軸由塔柱中垂線豎直向上，x軸由盤面法向在z軸法平面投影方向，與無窮遠處來流方向保持一致，如圖1 所示。

圖1 剛體模型坐標系Fig. 1 Coordinate system of rigid body

考慮浮式風機系統葉片受到的氣動載荷(Fblade，考慮變槳)、塔柱上的風壓載荷(Ftower)以及作用在浮式基礎上的波浪力(Fwave，包括一階力、二階力)及系泊力(Fmoor)，寫出系統單剛體動力學方程為：

式中：M為浮式風機質量矩陣；m∞為無窮頻率下浮式基礎附加質量矩陣；C為阻尼矩陣，浮式基礎粘性阻尼的大小取臨界阻尼的5%～8%；K為靜水恢復剛度矩陣；R為 速度脈沖函數矩陣，其與輻射阻尼矩陣B和附加質量矩陣m的關系如下：

1.2 剛-柔耦合多體動力學方程

基于Jourdain 速度變分原理，考慮剛-柔耦合建立浮式風機系統動力學方程，其中浮式基礎處理為剛體，考慮六自由度運動；將塔柱和葉片處理成彈性體，采用有限元方法進行離散；機艙處理為塔柱末端的集中質量；輪轂處理為剛體。浮式風力機系統模型如圖2 所示。

圖2 剛-柔耦合模型坐標系Fig. 2 Coordinate system of R-F coupling body

建立多體系統坐標系如下：大地坐標系(O0,e(0))，即慣性坐標系；浮式基礎隨體坐標系(O1,e(1))，坐標原點位于浮式基礎重心；塔柱浮動坐標系(O2,e(2))，坐標原點位于塔柱底端，描述塔柱變形；輪轂隨體坐標系(O3,e(3))，坐標原點位于輪轂中心，隨輪轂轉動；葉片浮動坐標系(O4,e(4))、(O5,e(5))、(O6,e(6))，固結于葉片根部，描述葉片變形。表示在相應坐標系上的基矢量。

根據卡爾丹角坐標轉換定理，物體的旋轉可分解為先繞x軸，再繞y軸，最后繞z軸的三維旋轉運動，故從一坐標系轉換到另一坐標系時，需左乘卡爾丹角方向余弦矩陣：

不同坐標系上的角速度相互轉換時也需要左乘相應的角速度轉換矩陣：

式中：與分別代表對相應角度的正、余弦計算。

利用有限元原理對塔柱與葉片進行離散。根據材料平斷面假定，柔性梁上任一點k處的變形位移 uk可由中線上對應點的變形位移得到，即：

式中，ug為考慮該點縱向變形位移的二次耦合項。

據連續介質力學理論，k點處的縱向正應變為：

由上述二式即可得到柔性梁各方向應變及變形。通過Jourdain 速度變分原理推導出系統空間變形動力學方程。對于柔性體變形為：

對于剛體，根據動量矩定理有：

得到浮式風機系統的動力學方程為：

其中：δWe為外力虛功；δP為由柔性梁彈性虛功率，表達式為

考慮不同結構之間的約束條件，最終得到離散后的剛-柔耦合動力學方程：

其中： ()q表示Jacobian 矩陣； ()t表示對t求偏導數；M為廣義質量陣；Q為廣義力陣； Φq,λ和γ分別為約束方程的雅克比矩陣，拉格朗日乘子，以及加速度約束方程右項；q為廣義坐標，包括浮式基礎運動，塔柱變形，以及葉片質量點的運動。

廣義力列陣可表示為：

其中，平臺外力列陣Fplat根據公式(1)得到，包含波浪力Fwave與系纜力Fmoor兩部分。輪轂外力列陣是通過位移邊界條件，依照約束方程(11)計算得到的。塔柱與葉片由于采用有限元方法離散，其各個分段ni上所受均布力需根據有限元離散規則積分得到，如下式：

式中，fi為相應位置的均布載荷強度。

2 外載荷計算

在海洋環境中，影響風機運動響應的外載荷主要為風載荷、波浪載荷、系纜力等。

2.1 風載荷

所有位于水面以上的浮體結構均受到風載荷的影響。對于塔柱及停機狀態的葉片，采用繞流理論計算風壓，即：

式中：fL,D為升力、阻力風壓，CL,D為截面升力與阻力系數，ρ為氣流密度，vw為來流風速。對于剛柔耦合問題，fL,D經過坐標轉換后就可得到式(13)中的均布載荷強度fi；對于剛體問題，還需積分求得整根梁上的繞流風力，包括升力FL與阻力FD：

對于一般作業工況下的葉片，風載荷主要表現為氣動載荷，本文使用經典葉素動量理論來解決此問題，具體求解過程與理論可參考文獻[4]。

2.2 波浪載荷

對于波浪載荷，基于三維勢流理論計算，采用Sesam/Wadam 模塊進行模擬。在Sesam 軟件中建模計算，得到一階波浪力、平均漂移力、附加質量、靜水恢復剛度、勢流阻尼的水動力傳遞函數后，將頻域結果轉換為時域結果，并將其分別代入剛體程序與剛-柔耦合程序中，之后進行浮式風機系統動力響應的時域計算。

2.3 系纜載荷

風機的系泊系統由3 根系泊纜組成，間隔120°設置，如圖3 所示。錨泊點處于水下200 m，導纜孔的位置在浮體水下14 m 的位置處，整體的錨鏈未被拉伸時長837.6 m，具體的系纜參數見文獻[10]。

圖3 系纜布置方式Fig. 3 Mooring arrangement

使用準靜態懸鏈線理論[3]求解系纜載荷，通過分類討論判別錨鏈躺底狀態并使用迭代方法求解。

2.4 變槳控制策略

依照參考規范[10]中的槳距角-風速控制策略得到低風速下的槳距，并根據計算結果優化擬合得到高風速下的槳距，從而得到各個風速下的最佳槳距，結果如表1 所示。

表1 槳距角控制Tab. 1 Pitch control

通過使用槳距角控制，浮式風機可以在較低風速下達到最高輸出功率的同時，在較高風速下又不至于葉片損壞、系纜斷裂或基礎位移過大，從而實現在各種環境條件下浮式風機系統的數值模擬。

2.5 浮式風機系統運動響應計算過程

對于單剛體計算，每一迭代步計算風機氣動載荷，之后考慮氣動載荷影響計算浮式風機整體運動，而浮式基礎運動又影響下一時刻步的氣動載荷計算，從而實現浮式風機系統氣動力和浮式基礎運動的相互耦合。

對于剛-柔耦合計算，每一步迭代計算風機氣動載荷，之后考慮氣動載荷影響計算浮式基礎運動、葉片/塔柱振動，而浮式基礎運動、葉片/塔柱振動又影響下一時刻的氣動載荷計算，實現浮式風機系統氣動力-浮式基礎運動-葉片/塔柱振動的相互耦合。

單剛體計算程序和多剛體計算程序都采用相同的風載荷、波浪載荷、系纜力計算方法及變槳控制策略。對2 種程序的正確性做了大量的驗證，見參考文獻[6,11]，本文直接采用2 種程序開展計算分析。

3 響應結果與分析

參考規范[10]中給出的切入、額定、切出、極限4 個工況對單剛體建模方法與剛-柔耦合建模方法分別得到的系統運動響應進行對比分析，各個工況的環境參數如表2 所示。

表2 工況環境參數Tab. 2 Working conditions

風、浪作用方向同系纜1 的布置方向，即作用在浮式風機系統的縱蕩方向，因此主要針對浮式風機整體運動較為明顯的縱蕩、縱搖自由度以及直接影響輸出功率的氣動轉矩進行分析。

3.1 不同工況的響應對比

分別采用單剛體算法和剛-柔耦合算法對不同工況下浮式風機系統響應進行計算。其中，對于切出工況，風機從穩定運行狀態到降速乃至停機狀態的轉捩點，但為防止葉片受損或因機械功率過大導致電機受損，風機葉片需要進行較大的變槳，依照規范，在23°槳距下對單剛體及耦合2 種程序進行模擬計算；對于極限工況，風輪停轉，槳距角調整為90°附近（順槳）。得到浮式風機系統縱蕩、縱搖及氣動轉矩的時間歷程和頻譜圖，如圖4～圖15 所示。對不同工況的響應結果進行統計，得到其平均值和標準差，并進行對比，結果如圖16～18 所示。

圖4 切入工況縱蕩響應Fig. 4 Surge, LC1

圖5 切入工況縱搖響應Fig. 5 Pitch, LC1

圖6 切入工況氣動轉矩Fig. 6 Torque, LC1

圖7 額定工況縱蕩響應Fig. 7 Surge, LC2

圖8 額定工況縱搖響應Fig. 8 Pitch, LC2

圖9 額定工況氣動轉矩Fig. 9 Torque, LC2

圖10 切出工況縱蕩響應Fig. 10 Surge, LC3

圖12 切出工況氣動轉矩Fig. 12 Torque, LC3

圖13 極限工況縱蕩響應Fig. 13 Surge, LC4

圖14 極限工況縱搖響應Fig. 14 Pitch, LC4

圖15 極限工況轉矩Fig. 15 Torque, LC4

圖16 各工況縱蕩統計值對比Fig. 16 Surge comparison

圖17 各工況縱搖統計值對比Fig. 17 Pitch comparison

圖18 各工況轉矩統計值對比Fig. 18 Torque comparison

針對不同工況進行分別分析，分析中用到的相對差異均指剛-柔耦合方法相對于單剛體方法的差異，即以上計算表明：

1）對于切入工況，風浪環境條件較為溫和，風機以低功率轉動，槳距角為0°，浮式風機的縱蕩、縱搖運動較小，2 種方法所得氣動轉矩略有差異，而縱蕩和縱搖運動一致性較好，以下重點分析另外3 個工況。

2）浮式風機縱蕩/縱搖運動及風機氣動轉矩的平均值采用兩種不同模型得到的結果非常接近；對于極限工況，2 種模型得到的浮式風機縱搖平均值相對差異較大，但絕對差異僅為0.5°。因此，浮式風機系統的縱蕩/縱搖運動及風機氣動轉矩的平均值非常接近，兩種不同模型的計算結果差異很小。

3）2 種不同模型計算得到的浮式風機縱蕩/縱搖及氣動功率標準差的相對差異較為顯著。對于浮式風機的縱蕩運動，在額定和切出工況，2 種模型的縱蕩標準差相對差異較小，剛體模型結果略小于剛-柔耦合模型結果；而極限工況結果相反，和剛體模型結果相比，剛-柔耦合模型得到的縱蕩標準差減小約20%。對于浮式風機的縱搖運動，剛體模型得到的縱搖運動標準差均大于剛-柔耦合模型的結果，相對差異最大為極限工況約為30%。這是由于柔性葉片和塔柱的振動吸收了部分外載荷能量，使得剛-柔耦合模型的結果小于剛體模型的結果。對于氣動轉矩，2 種模型得到的標準差差異較小，額定工況的相對差異最大為10%左右。

3.2 計算效率統計

為研究不同計算方法的計算效率，本文各工況均選用Matlab 算法進行數值模擬，模擬時長均為1 000 s，最終得到的計算效率統計對比如表3 所示。

表3 計算效率對比Tab. 3 Efficiency comparison

可知：切入工況耗時最長，這是由于在低風速時BEM 方法迭代收斂需要的時間較長；剛體模型相對剛-柔耦合模型的計算效率優勢極為明顯，如切入工況，剛體模型的計算時間僅為剛-柔耦合模型計算時間的10%左右。這是由于剛-柔耦合程序不僅計算了整體載荷、位移，還在每個時間步中計算塔柱與葉片變形，并將此變形代入氣動模塊進行耦合。這一過程使得結構動力學矩陣維度陡增，從而求解難度與消耗時長激增。

實際設計中，應充分考慮2 種計算模型的特點，合理使用2 種計算模型。在浮式風機系統的動力學參數設計初始階段，由于剛體模型計算的快速性，可以用來獲得浮式風機系統運動及氣動轉矩，通過大量參數組合計算及方案篩選，確定初步的浮式風機系統動力學設計參數。之后，針對已有的浮式風機系統設計參數，采用剛-柔耦合模型開展進一步的校核，以獲得更為精確的計算結果。

4 結 語

本文以OC4-NERL5 MW 浮式風機為例，建立浮式風機系統剛體動力學模型及剛-柔耦合動力學模型；考慮不同工況，對比研究了剛體模型及剛-柔耦合模型2 種方法在計算浮式風機系統縱蕩/縱搖及氣動轉矩的差異，分析了2 種模型的實用性。研究結論如下：

1）對于浮式風機系統縱蕩/縱搖運動，切入工況下浮式風機的運動較小，2 種方法所得的風機系統縱蕩/縱搖的平均值和標準差都非常接近；對于其他工況，采用2 種不同模型得到的浮式風機縱蕩/縱搖的平均值結果也非常接近。

2）除去切入工況，2 種模型計算得到的浮式風機運動標準差的差異較為顯著，相對差異最大為極限工況。和剛體模型結果相比，剛-柔耦合模型的縱蕩標準差最大減小約20%，縱搖標準差最大減小約30%。彈性葉片和塔柱振動可吸收部分外載荷能量，使得剛-柔耦合模型得到的運動標準差小于剛體模型結果。

3）對于所有的工況，2 種模型得到的氣動轉矩均值幾乎相同；標準差略有差異，額定工況氣動轉矩標準差的相對差異最大，約10%左右。

(4)實際設計中應合理使用2 種計算模型，剛體模型可用于浮式風機系統的動力學參數設計初始階段，進行多組參數組合及方案篩選的大規模計算；剛-柔耦合模型可用于后續的動力學參數校核，以獲得更為精確的計算結果。

本文計算中，為防止槳距角不同對2 種模型計算產生的影響，采用了定槳計算，與實際情況存在一定差異。另外，本文計算中，只考慮了風浪同向的情況，后續將進一步研究風浪異向的情況。