不同動力學(xué)建模方法對浮式風(fēng)機系統(tǒng)動力響應(yīng)的影響研究

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

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

0 引 言

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

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

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

1 浮式風(fēng)機系統(tǒng)動力學(xué)方程

本文以5 MW 浮式風(fēng)機系統(tǒng)為例進行分析，其中，浮式基礎(chǔ)為OC4 基礎(chǔ)，上部為NERL-5 MW 風(fēng)力機，具體參數(shù)見文獻[10]。

1.1 單剛體動力學(xué)方程

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

圖1 剛體模型坐標(biāo)系Fig. 1 Coordinate system of rigid body

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

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

1.2 剛-柔耦合多體動力學(xué)方程

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

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

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

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

不同坐標(biāo)系上的角速度相互轉(zhuǎn)換時也需要左乘相應(yīng)的角速度轉(zhuǎn)換矩陣：

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

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

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

據(jù)連續(xù)介質(zhì)力學(xué)理論，k點處的縱向正應(yīng)變?yōu)椋?/p>

由上述二式即可得到柔性梁各方向應(yīng)變及變形。通過Jourdain 速度變分原理推導(dǎo)出系統(tǒng)空間變形動力學(xué)方程。對于柔性體變形為：

對于剛體，根據(jù)動量矩定理有：

得到浮式風(fēng)機系統(tǒng)的動力學(xué)方程為：

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

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

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

廣義力列陣可表示為：

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

式中，fi為相應(yīng)位置的均布載荷強度。

2 外載荷計算

在海洋環(huán)境中，影響風(fēng)機運動響應(yīng)的外載荷主要為風(fēng)載荷、波浪載荷、系纜力等。

2.1 風(fēng)載荷

所有位于水面以上的浮體結(jié)構(gòu)均受到風(fēng)載荷的影響。對于塔柱及停機狀態(tài)的葉片，采用繞流理論計算風(fēng)壓，即：

式中：fL,D為升力、阻力風(fēng)壓，CL,D為截面升力與阻力系數(shù)，ρ為氣流密度，vw為來流風(fēng)速。對于剛?cè)狁詈蠁栴}，fL,D經(jīng)過坐標(biāo)轉(zhuǎn)換后就可得到式(13)中的均布載荷強度fi；對于剛體問題，還需積分求得整根梁上的繞流風(fēng)力，包括升力FL與阻力FD：

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

2.2 波浪載荷

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

2.3 系纜載荷

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

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

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

2.4 變槳控制策略

依照參考規(guī)范[10]中的槳距角-風(fēng)速控制策略得到低風(fēng)速下的槳距，并根據(jù)計算結(jié)果優(yōu)化擬合得到高風(fēng)速下的槳距，從而得到各個風(fēng)速下的最佳槳距，結(jié)果如表1 所示。

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

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

2.5 浮式風(fēng)機系統(tǒng)運動響應(yīng)計算過程

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

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

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

3 響應(yīng)結(jié)果與分析

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

表2 工況環(huán)境參數(shù)Tab. 2 Working conditions

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

3.1 不同工況的響應(yīng)對比

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

圖4 切入工況縱蕩響應(yīng)Fig. 4 Surge, LC1

圖5 切入工況縱搖響應(yīng)Fig. 5 Pitch, LC1

圖6 切入工況氣動轉(zhuǎn)矩Fig. 6 Torque, LC1

圖7 額定工況縱蕩響應(yīng)Fig. 7 Surge, LC2

圖8 額定工況縱搖響應(yīng)Fig. 8 Pitch, LC2

圖9 額定工況氣動轉(zhuǎn)矩Fig. 9 Torque, LC2

圖10 切出工況縱蕩響應(yīng)Fig. 10 Surge, LC3

圖12 切出工況氣動轉(zhuǎn)矩Fig. 12 Torque, LC3

圖13 極限工況縱蕩響應(yīng)Fig. 13 Surge, LC4

圖14 極限工況縱搖響應(yīng)Fig. 14 Pitch, LC4

圖15 極限工況轉(zhuǎn)矩Fig. 15 Torque, LC4

圖16 各工況縱蕩統(tǒng)計值對比Fig. 16 Surge comparison

圖17 各工況縱搖統(tǒng)計值對比Fig. 17 Pitch comparison

圖18 各工況轉(zhuǎn)矩統(tǒng)計值對比Fig. 18 Torque comparison

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

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

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

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

3.2 計算效率統(tǒng)計

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

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

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

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

4 結(jié) 語

本文以O(shè)C4-NERL5 MW 浮式風(fēng)機為例，建立浮式風(fēng)機系統(tǒng)剛體動力學(xué)模型及剛-柔耦合動力學(xué)模型；考慮不同工況，對比研究了剛體模型及剛-柔耦合模型2 種方法在計算浮式風(fēng)機系統(tǒng)縱蕩/縱搖及氣動轉(zhuǎn)矩的差異，分析了2 種模型的實用性。研究結(jié)論如下：

1）對于浮式風(fēng)機系統(tǒng)縱蕩/縱搖運動，切入工況下浮式風(fēng)機的運動較小，2 種方法所得的風(fēng)機系統(tǒng)縱蕩/縱搖的平均值和標(biāo)準(zhǔn)差都非常接近；對于其他工況，采用2 種不同模型得到的浮式風(fēng)機縱蕩/縱搖的平均值結(jié)果也非常接近。

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

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

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

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