基于Jourdain原理和有限元離散的浮式風機動力響應分析

劉利琴，肖昌水，郭穎，李焱

(1.天津大學 水利工程仿真與安全國家重點實驗室，天津 300072； 2.天津大學 天津市港口與海洋工程重點實驗室，天津 300072)

近年來，化石能源帶來的環境問題日趨嚴重，風能作為可再生清潔能源受到人們的關注。陸上風機技術發展較為成熟，海上風機也因其優勢得到迅猛發展，風機逐漸由陸地向海洋發展。出于經濟效益考慮，水深超過50 m后宜采用浮式基礎，包括半潛型、張力腿型(tension leg platform，TLP)、單柱式(SPAR)等。

針對浮式風機已有大量研究。Koo等[1]針對3種基礎型式的風機進行模型試驗。Shin等[2]針對OC3-Hywind浮式風機進行了模型試驗并分析了風機系統運動特性。Jonkman等[3]開發了氣動-水動-伺服-彈性全耦合浮式風機仿真程序，并分析了5 MW風機動力特性，Robertson等[4]運用該仿真程序對幾種不同基礎型式的浮式風機進行環境載荷分析。Ren等[5]基于計算流體動力學(computational fluid dynamics, CFD)方法計算并分析了風浪聯合作用下時域運動響應。Li等[6]考慮下沉運動研究了二階波浪力和氣動載荷對TLP型浮式風機動力響應的影響。

目前大部分研究將浮式風機系統考慮為單剛體。本文將浮式風機系統處理為剛-柔耦合的多體系統，其中基礎、輪轂為剛體，塔柱和葉片為柔性體，并考慮了葉片軸向變形對彎曲變形(揮舞和擺振)的影響。采用Jourdain速度變分原理，采用有限單元法離散柔性結構，推導浮式風機剛-柔耦合動力學方程，采用自動變步長的向后差分法求解微分-代數方程組，并依此編制Matlab仿真程序。建立了剛-柔耦合多體系統動力學方程，分析了浮式風力機系統的動力響應。

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

1.1 坐標系的建立

為了分析海上浮式風機系統耦合動力響應，本文簡化整機系統，將機艙考慮為塔柱末端質量模型，整機系統分為浮式基礎剛體(B1)、塔柱柔體(B2)、輪轂剛體(B3)、葉片柔體(B4、B5、B6)。其中剛體部分考慮6個自由度，引入如圖1所示坐標系：1)大地坐標系(O0,e(0))，慣性坐標系；2)浮式基礎隨體坐標系(O1,e(1))，固結于浮式基礎重心；3)塔柱浮動坐標系(O2,e(2))，固結塔柱底端，描述塔柱變形；4)輪轂隨體坐標系(O3,e(3))，固結于輪轂中心，隨輪轂轉動；5)葉片浮動坐標系(O4,e(4))、(O5,e(5))、(O6,e(6))，固結于葉片根部，描述葉片變形。

圖1 浮式風機示意Fig.1 Schematic diagram of floating wind turbine

風機系統各體之間的轉動位移采用卡爾丹角描述，即繞x軸轉動角α、繞y軸轉動角β、繞z軸轉動角γ，表征轉動關系的方向余弦矩陣為[7]：

(1)

式中：c表示cos；s表示sin。

角速度ω表示在隨體坐標下坐標列陣，可用相對慣性基轉動角的導數表示，即：

(2)

(3)

1.2 柔性梁運動學關系描述

設r0為柔性梁浮動基到慣性基的矢徑，ρ為柔性梁上任意點k相對浮動基的矢徑，ρ0為未變形時k點關于浮動基的矢徑，u為柔性梁關于浮動基的彈性位移，有：

r=r0+ρ

(4)

ρ=ρ0+u

(5)

因此，柔性梁上k點關于慣性基的速度為:

(6)

式中：0()表示慣性基；b()表示浮動基；°表示對浮動基的局部導數。

柔性梁上k點關于慣性基的加速度為:

(7)

1.3 柔性梁變形位移描述

設ui(i=1,2,3)為k點的變形位移在浮動基的坐標分量，根據平斷面假設ui(i=1,2,3)可以用中線上對應點k0的變形位移wi(i=1,2,3)表示:

(8)

式中，在小變形假設下，考慮k點的縱向變形位移的二次耦合項[8]：

對于作大范圍運動的剛-柔耦合動力學系統，二次耦合項對剛度項的貢獻不可不計。就風力機旋轉風輪而言，不考慮二次耦合項可能導致數值仿真結果發散。

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

(9)

不計剪切及扭轉效應的柔性梁彈性虛功率為:

(10)

1.4 空間柔性梁有限元離散

目前，變形體的離散技術大致可分為3類：假設模態法、有限元法和部件模態綜合法[9]。為了反映旋轉葉片的真實振型和頻率，本文采用有限元法離散塔柱和葉片。

(11)

其中：

因此，中線上對應點k0的變形位移wi表示為:

(12)

每個單元的節點變形位移列陣pj可通過布爾矩陣Bj及邊界條件轉化矩陣R從總體變形位移列陣p中“拾取”，即：

pj=BjRp

(13)

式中：布爾矩陣Bj和邊界條件轉化矩陣R分別為：

本文柔性梁均為一端固支的懸臂梁，因此節點1的位移全為0，總體變形位移列陣p表示為：

p=[w12w22θ22w32θ32…w1n+1w2n+1

θ2n+1w3n+1θ3n+1]T

(14)

式中：元素下標第1位表示方向，第2位表示節點數；n為單元數。

將式(11)～(13)代入式(8)，并令Ni=Nj,iBjR，則k點相對浮動基的變形位移可表示為

(15)

式中耦合形函數為[10]：

1.5 空間梁剛-柔耦合動力學方程

本文通過Jourdain速度變分原理推導空間動力學方程，其形式為：

(16)

式中δWe為外力虛功率。

設廣義坐標列陣為：

(17)

現代大型風力機都安裝伺服控制系統，對于變槳調節的風力機在未達到額定風速之前，以輸出功率最大的轉速運行，當達到額定風速之后便以額定轉速運行。因此需要在輪轂與塔柱之間施加約束條件形如下，其他連接處設為剛性固定：

(18)

式中：下標t代表塔柱；h代表輪轂；At為塔柱浮動基相對慣性基的余弦矩陣；ρt為考慮塔柱變形的塔柱浮動基到輪轂浮動基的矢徑；ω為風輪轉速。

因此，約束方程可表示為:

(19)

將式(6)、(7)、(11)、(15)代入動力學方程式(16)，對于剛體方程中的變形位移和彈性虛功率為零，并考慮約束方程得到形如下式的第一類拉格朗日方程：

(20)

式中：()q表示Jacobian矩陣；()t表示對t求偏導數；λ為拉格朗日乘子；M為廣義質量陣；Q為廣義力陣。式(20)是指標為1的微分-代數方程組，采用傳統增廣法求解，消去拉格朗日乘子，考慮位置約束和速度約束的違約問題，并采用向后差分法(backward differentiation formulas, BDFs)求解增廣后的常微分方程組。本文數值求解為變步長積分，積分精度為10-6。

2 系統載荷計算

2.1 風機葉片非定常氣動載荷

本文將經典葉素動量理論(blade element momentum, BEM)擴展到局部，計算每個葉素的非定常氣動載荷，考慮風剪切、塔影效應、葉尖輪轂損失因子及Glauert修正，在此基礎上引入尾流傾斜模型、B-L動態失速模型，經過修正的BEM能計算得到較為準確的氣動載荷，滿足本文計算要求。

由圖2葉素速度三角關系可得[11]：

圖2 葉素速度及作用力三角關系Fig.2 The trigonometric relationship of velocity and force of blade element

(21)

式中：Vrel,x、Vrel,y分別為來流風速疊加葉片運動速度的軸向和周向速度分量。

針對每個葉素，軸向誘導因子和切向誘導因子可表示為：

(22)

(23)

式中：B為葉片數；Cl、Cd為升阻系數；φ為入流角；F為葉尖/輪轂損失因子：

(24)

式中：R為風輪半徑；Rhub為輪轂半徑。

Buhl基于Glauert修正模型，提出當a≥0.4時，推力系數及軸向誘導因子的表達式[12]為：

(25)

(26)

以上推導過程基于來流風速正對風輪平面，當風輪發生偏航時，尾流將發生傾斜。對于海上浮式風機，當浮式基礎發生艏搖，風輪平面發生偏航。因此，必須引入尾流傾斜模型對軸向誘導因子修正。

軸向誘導因子偏航修正可表示為：

(27)

式中：γ為偏航角；θ為風輪方位角；θ0為葉片指向尾流最深處時的方位角。

尾流傾角可用下式近似計算：

χ=(0.6a+1)γ

(28)

考慮風剪切及塔影效應，計算入流速度。風剪切模型可表示為：

(29)

式中：H為輪轂高度；vH為輪轂處風速；η為風切變指數，通常取0.1～0.25。

塔柱的存在影響了尾流分布，如圖3所示，對于上風型風力機可采用潛流模型[13]：

圖3 風剪切和塔影效應示意Fig.3 Schematic diagram of wind shear and tower shadow

(30)

式中：DT為葉素所在高度塔柱直徑；x、y為葉素與塔柱中心的距離。式(30)在葉片處于輪轂中心以下，相對塔柱中心±60°區域內成立。當葉片處于輪轂中心之上，相對塔柱中心±60°區域內無塔影效應，其他區域修正為A(0.5-cosθ)+(0.5+cosθ)，使得各區域平滑過渡。

因此，考慮風剪切和塔影效應后的來流速度可表示為：

U∞=Av(z)

(31)

本文使用B-L動態失速模型進行動態失速的模擬，分為非定常附著流、尾緣流動分離及動態失速3個過程[14]。編制Matlab程序，模擬S809翼型動態失速過程，翼型攻角以正弦規律振動α=15+10sin(ωt)，衰減頻率k=0.051 ，馬赫數M=0.11，并與靜態風洞試驗數據對比[15],如圖4所示。

圖4 動態失速對升阻系數的影響Fig.4 The effect of dynamic stall on lift and drag coefficient

從圖4可以發現，動態失速體現在升力失速延遲，在失速區動態氣動系數與靜態數據有明顯差異。攻角增加變化時，升力持續增加并且要大于靜態升力數據；攻角回程階段，升力小于靜態升力數據，阻力系數也有類似的現象，攻角持續變化氣動系數形成遲滯環。氣動仿真時，基礎的運動以及葉片自身變形使得翼型攻角持續變化，因此有必要考慮動態失速的影響。

2.2 浮式基礎系泊力/矩

本文以外載荷的方式考慮系泊系統對浮式風機系統動力學方面的貢獻。以浮式風力機OC4-DeepCwind為例，該風機為懸鏈線式系泊，如圖5。懸鏈線式系泊根據其形狀，忽略水動力的影響，可通過懸鏈線方程快速估計系泊張力及系泊恢復力(矩)[16]。浮式基礎六自由度系泊恢復力(矩)隨縱蕩位移變化趨勢如圖6所示。

圖6 系泊六自由度恢復力和力矩Fig.6 Restoring force and moment of mooring system

圖5 系泊系統示意Fig.5 Schematic diagram of mooring system

2.3 浮式基礎水動力

在海洋工程中，小尺度構件用Morison方程計算，大尺度結構物通過三維勢流理論計算。本文采用SESAM/Wadam模塊進行水動力方面參數的計算，包括一階波浪力傳遞函數、附連水質量、勢流阻尼、靜水恢復剛度。勢流理論不計及粘性阻尼，工程中認為粘性阻尼可取臨界阻尼β0的5%～8%，即：

(32)

式中：M0為浮式基礎質量；Ma為浮式基礎附連水質量；Ci為浮式基礎第i自由度的靜水恢復剛度。

浮式基礎時域控制方程為：

(33)

3 半潛式浮式風機算例分析

本文以OC4-DeepCwind浮式風機系統為計算模型，相關參數可參考文獻[17-18]。綜合上述建模過程，可得到如圖7所示程序框圖，并編制Matlab仿真程序。

圖7 仿真程序框圖Fig.7 Simulation block diagram

3.1 柔性葉片固有振動特性

本文采用一次近似耦合模型，考慮二次耦合變形量，體現了葉片離心力和軸向慣性力對揮舞和擺振方向變形的影響，在動力剛度項中產生的附加耦合剛度項。考慮式(16)中δWe=0，且大范圍運動已知，則方程退化為非慣性系下自由振動方程。圖8為單根葉片在轉速為12.1 r/min時，劃分不同單元數目計算一階揮舞固有頻率的結果。

圖8 驗證本文方法單元收斂性Fig.8 Verification of the convergence of element

圖8結果顯示，葉片劃分至少9個單元就能得到收斂的計算結果。采用本文程序計算零轉速下葉片自由振動，葉片劃分10個單元，分析其固有頻率，并與FAST和ADAMS計算結果對比[17]，如表1所示，本文程序結果與二者吻合較好。

表1 本文程序結果驗證Table 1 Verification of results by computer program Hz

為了揭示葉片大范圍運動產生的“動力剛化”現象，將浮式基礎固定并且不考慮塔柱柔性變形，采用本文程序給定轉速計算葉片自由振動，固有頻率隨轉速變化如圖9所示。

圖9 葉片固有頻率隨轉速變化關系Fig.9 The frequencies of blade vary with rotor speed

由圖9可以發現：一次近似耦合模型計算的揮舞和擺振方向固有頻率均隨著轉速的增大而增加，并且轉速越高，固有頻率增大越快。固有頻率的增大，反映出結構剛度的增大。而傳統零次近似模型的擺振方向固有頻率隨轉速的增大呈減小的趨勢，揮舞方向的剛度不隨轉速變化。這是因為揮舞方向在旋轉平面外，擺振方向在旋轉平面內，傳統零次近似模型忽略始終為正的附加耦合剛度項。一次近似耦合模型考慮了附加耦合剛度項，保證了旋轉葉片各方向的剛度不隨轉速增大而減小，旋轉葉片出現了“動力剛化”現象。

3.2 浮式風機系統動力響應

在額定工況下(轉速12.1 r/min，風速11.4 m/s)，選定海況規則波波高6 m、周期10 s，波浪和風向均為0°，采用本文程序計算浮式風機系統動力響應。

圖10為葉尖揮舞和擺振方向變形，從圖中可以看出，一次近似耦合模型的變形量明顯小于傳統零次近似模型，結果與上述分析的旋轉葉片特性相對應，考慮了葉片的“動力剛化”效應。如果采用傳統零次近似模型，隨著轉速不斷增大，葉片擺振方向的剛度不斷減小，有可能導致葉片變形的仿真失敗，與實際情況不符。

圖10 葉尖變形時域曲線Fig.10 Blade tip deformation in time domain

通過頻譜分析可以發現，葉片軸向變形二次耦合項只對葉片變形幅值有所影響，并不改變其運動規律。葉片揮舞方向振動頻率以波浪和1P(風輪轉速)為主，這是由于氣動載荷耦合了基礎的運動，受波浪影響較大。此外，揮舞變形還存在較為明顯的波浪頻率與1P的耦合頻率、2P以及3P等頻率成分。與葉片揮舞方向不同，擺振方向受到的切向氣動載荷較小，因此波浪頻率響應較小，而擺振方向在風輪旋轉過程中受交變重力影響較大，因此擺振方向變形主要以重力引起的1P響應為主導。

圖11 葉尖變形頻譜分析Fig.11 Spectral analysis of blade tip deformation

由于篇幅關系，本文只列出風浪作用平面的3個基礎自由度響應，即縱蕩、垂蕩、縱搖?；A縱蕩、縱搖和垂蕩的運動響應，如圖12所示。

從圖12中可以發現，柔性葉片的2種近似模型對浮式基礎的三自由度運動影響不大，這是因為考慮葉片軸向變形二次耦合項對葉片的振動速度影響不大，如圖13所示，因此2種模型計算的氣動載荷相差無幾。相比較而言，零次近似模型縱搖平衡位置為2.70°，而一次近似模型縱搖的平衡位置為2.71°，一次近似耦合模型的縱搖平衡位置略有增大。

圖13 葉尖揮舞速度時歷曲線Fig.13 Blade tip flap velocity in time domain

圖12 浮式基礎動力響應Fig.12 The dynamic response of floating foundation

此外，對三自由度作頻譜分析，如圖14所示，可以發現，浮式基礎運動響應頻率基本上以波頻為主，氣動載荷對基礎運動響應頻率成分貢獻非常小。這是由于水平軸3根葉片在同一旋轉平面內，浮式基礎所受氣動推力相當于葉片沿展向積分求和，葉片之間相位120°，求和之后氣動載荷幅值相互抵消大大減小，與波浪載荷波動幅值相比幾乎可忽略不計。

圖14 浮式基礎運動響應頻譜分析Fig.14 Spectral analysis of response of foundation

4 結論

1) 通過計算額定轉速下葉片固有頻率，驗證了本文建模方法的單元收斂性。針對5 MW風機葉片至少劃分9個單元可得到收斂的計算結果。僅考慮葉片運動，計算零轉速下葉片固有頻率，與FAST和ADAMS結果對比吻合較好。

2) 僅考慮葉片運動，計算不同轉速下葉片振動，分析其固有頻率。傳統零次近似模型計算的葉片擺振方向固有頻率隨轉速增大而減小，揮舞方向不發生變化。一次近似耦合模型考慮旋轉葉片的“動力剛化”效應，剛度隨轉速增大而增大，與實際情況相符。

3) 風浪聯合作用下，由于一次近似耦合模型考慮了旋轉葉片的“動力剛化”效應，葉片振動幅值較小，對葉片振動速度影響較小。葉片揮舞主要以波浪頻率和氣動1P響應為主，此外還存在二者的耦合頻率以及2P、3P響應，擺振方向變形主要受交變重力影響。

4) 在定常風與規則波工況下，氣動載荷對基礎運動響應頻率成分貢獻非常小，基礎運動響應以波浪頻率為主。事實上在某些復雜海況下，氣動載荷對基礎運動的影響明顯增大，因而考慮風機與浮體的耦合效應，對評估風機系統的運動性能至關重要。