考慮流固耦合的近海風機動力響應數值計算*

鄧 露，肖志穎，黃民希,宋曉萍，吳海濤

(1.湖南大學 土木工程學院，湖南 長沙 410082; 2. 湘電風能有限公司，湖南 湘潭 411102)

考慮流固耦合的近海風機動力響應數值計算*

鄧 露1?，肖志穎1，黃民希1,宋曉萍2，吳海濤2

(1.湖南大學 土木工程學院，湖南 長沙 410082; 2. 湘電風能有限公司，湖南 湘潭 411102)

研究了固定式海上風機的流固耦合問題.選取美國可再生能源實驗室(NREL)建立的5 MW海上單樁式風機模型，采用Kaimal譜和風剪切指數理論模型，通過NREL編寫的程序—TurbSim，生成脈動風場，進而使用氣彈程序FAST計算風力機的氣動荷載，輸出包括剪力和彎矩在內的塔頂荷載時程.選取Pierson-Moskowits波浪譜，并利用諧波疊加法模擬隨機波浪.采用Morison公式計算波浪力，分別建立了海上風機塔架考慮與不考慮流固耦合效應時的動力方程及數值模擬方法，并對比分析了這兩種情況下的波浪力和結構響應的差異.結果表明，兩種情況下計算得到的波浪力和風機響應有明顯差異，并且這種差別跟風浪大小等因素有關.

海上風機；隨機風浪；流固耦合；Morison方程

在海洋工程中，小尺度(直徑相對于波長較小)圓柱結構應用相當廣泛，且波浪荷載是其設計的控制荷載.因此，研究波浪荷載對此類結構物的作用具有重要的工程意義.1950年Morison[1]等提出了一個計算小尺度結構所受波浪力的經驗公式，該公式形式簡單，能夠滿足一般工程需要，至今仍廣泛使用.但當結構本身的振動較大時，結構的運動會改變波浪力的大小，結構與波浪之間的流固耦合效應不能忽略.Sarpkaya[2]等人針對這一問題將Morison方程進行了修改，提出了考慮相對運動的修正Morison方程.Fish[3]等人針對修正Morison方程中的流固耦合問題進行了試驗和理論研究，提出了相應的簡化計算方法.國內方面，王元戰[4]等研究了離岸深水碼頭樁基所受波浪力的流固耦合問題；徐亞洲[5]等提出了一階預估校正法來研究結構流固耦合波浪動力響應.

對于深海結構物，結構的柔度較大，其自身運動對波浪的影響較大.因此，計算波浪力時應該考慮結構自身的運動，即流固耦合問題.比如Chandrasekaran[6]等應用修正Morison方程進行深海結構的波浪力計算，采用Monte-Carlo方法模擬隨機波浪荷載，研究了張拉腿平臺結構在隨機波浪作用下的結構響應.而對于近海結構物，結構剛度較大，計算波浪力時一般不考慮流固耦合.

然而，修正前和修正后的Morison公式兩者的應用范圍并沒有明確的界限.對于固定式海上風機結構，計算其所受的波浪力時可按DNV[7]規范的規定不考慮流固耦合.Shirzadeh[8]等對比研究了評估海上風機阻尼的實驗和數值方法，計算波浪力時沒有考慮流固耦合.但是也有學者考慮了流固耦合：Agarwal[9]等研究海上風機aero-servo-hydro-elastic(空氣動力-水動力-伺服)耦合模型時，考慮了流固耦合問題計算不規則波浪力；李德源[10]等研究了固定式海上風力機塔架在隨機風、浪載荷作用下的動力響應，考慮了流固耦合問題計算不規則波浪力.張學志[11]等人對海洋平臺的流固耦合問題進行了相關研究，發現是否考慮流固耦合對結構的位移和加速度等響應有較大影響，但其研究對象沒有涉及海上風機.

海上風機結構在振動控制方面有較高要求，振動過大會對風機造成損害，導致運行故障.GL[12]風機認證規范專門規定了振動限值.因此準確計算風機結構的動力響應不僅能夠指導工程設計，亦能對風機工作性能進行準確評估.由于風機自身結構高聳，且受到風浪流等多種復雜荷載作用，其結構的力學特性與海洋平臺、浮式風機有較大差異，是否需要考慮流固耦合問題，以及流固耦合問題對計算波浪荷載和結構的動力響應的影響程度，都值得深入研究.

針對上述問題，本文對比分析了固定式海上風機結構考慮與不考慮流固耦合問題的波浪力計算方法，并分別建立了兩種情況下的海上風機塔架動力方程及有限元數值模擬方法，以此探究海上風機波浪力計算中的流固耦合問題.并以NREL-5 MW(5兆瓦)風機模型為算例，對比分析了計算結果的差異，獲得了一些對工程實踐有意義的結論.

1 理論介紹

1.1 風機模型

固定式海上風機主要由基礎、塔架、發電機組構成.本文采用Ansys軟件建立有限元模型，將機艙和葉片等效成為塔頂的質點；用梁單元模擬塔架以及基礎結構；為簡化計算，將基底邊界條件假設為固接.風機結構受到的隨機風荷載和隨機波浪荷載作用在相應的質點.風機模型圖如圖1所示.

圖1 風機模型圖

1.2 風機結構動力學方程

本文計算風機x方向的動力響應，其動力學方程如下：

[Fwave]+[Fwind]

(1)

式中[M],[C]和[K]分別為風機結構的質量矩陣、阻尼矩陣和剛度矩陣；[Fwave]，[Fwind]分別為波浪荷載和風荷載.[C]采用Rayleigh阻尼，計算公式如(2)～(4)所示.動力學方程求解采用Newmark-β法求解，利用自編Matlab程序計算.

(2)

(3)

(4)

式中α1,β1為比例系數；ω1,ω2分別為第一、第二階自振頻率(rad)；ζ1，ζ2為對應的阻尼比.

1.3 波浪荷載模擬

1.3.1 隨機波浪模型

海上風機結構承受多種復雜荷載，其中波浪荷載是風機塔架的主要設計荷載[13].工程上有多種方法描述波浪，比如Stokes波理論、流函數理論等.由于海浪具有隨機性, 設計過程中必須考慮結構響應的隨機性.因此，本文采用隨機波浪理論描述波浪.工程上一般采用Jonswap,Neumann或Pierson-Moskowits[14]等譜來模擬隨機波浪.首先，根據近海風電場的風況和波候等條件得到海浪的統計參數，如有效波高、跨零周期等，以產生波浪的功率譜密度信號，然后利用諧波疊加法得出波高的時域信號，最后選取合適的波浪理論得到水質點的速度、加速度等運動參數.本文采用Pierson-Moskowits譜(簡稱P-M譜)進行計算，具體表達式如下：

(5)

式中Sη(ω)為波浪譜密度函數；ω為波浪頻率(rad)；α2，β2為常數，α2=0.008 1，β2=0.74；U為海面以上19.5 m的風速(m/s)；g為重力加速度(m/s2).

此波浪譜公式含有風速項，說明其蘊含了風浪之間的關系，即確定風速后就可以確定波浪譜信息.根據設計資料，可以確定式(5)中的具體參數，代入公式后，即可以得到能量譜分布，典型P-M譜函數曲線如圖2所示.

隨機波浪理論[15]將波浪視為由無限多個振幅不等、頻率不等、初相位隨機的余弦波疊加而成，利用諧波疊加法得到波高，理論計算公式如下：

(6)

式中η(t)為波高函數；t為時間；εn為0～2π之間均勻分布的初相位.

ω/rad

為進行數值計算，將自變量波浪頻率平均分成N個區間，每個區間長度為Δωn，用區間邊界的平均值作為該區間的代表值ωn，即可用數值計算公式(7)得到波高.

(7)

計算得到的典型波高時程曲線如圖3所示.

t/s

假設流體為無粘性不可壓縮的均勻流體，采用airy線性波浪理論來求解波浪.目前工程上隨機波的模擬都是基于線性波來模擬的，即Longuet-Higgins隨機波浪模型，這樣可得到波浪的速度、加速度等運動特性.airy波浪理論的計算公式表述如下.

波浪彌散關系：

ωn2=gkntanhknd

(8)

波浪速度場：

(9)

加速度場：

(10)

式中Hn為波高；Tn為波浪周期；d為水深；kn為波數；z為質點z方向坐標.通過式(9),(10)，可以得到相關質點波浪運動參數.

1.3.2 隨機波浪荷載

在得到波浪速度、加速度之后，即可計算海上風機的波浪荷載.采用Morison公式，將水平波浪荷載分為水平慣性力和水平拖拽力，當結構振動比較小時，風機任意質點波浪力計算公式如下：

(11)

式中fD為拖拽力；fI為慣性力；CD為拖拽力系數；CM為慣性力系數；ρ為海水密度；V0為單位柱體高度的排水體積.

當考慮流固耦合時，結構質點波浪力計算為：

(12)

將式(12)代入式(1)，得到質點動力學方程：

(13)

(14)

不難發現方程(14)等號右邊的波浪力項與方程(11)是一致的.考慮耦合的時候，只需要將質量矩陣相應的自由度加上一個附加質量系數即可，相對于張學志[11]等人的方法，計算相對簡單.最后將計算出風機塔架的速度，加速度代入式(12)就可以得到考慮流固耦合情況時的波浪力.

1.4 風荷載模擬

1.4.1 隨機風模型

風場中任一點的風速可分為兩部分：平均風速和脈動風速.前者需考慮在大氣邊界層范圍內的風剪切效應來模擬平均風速隨高度的變化；后者可視為零均值、平穩各態歷經的隨機過程，根據功率譜函數進行數值模擬.

風剪切模型有對數和指數模型，分別假定風速隨高度呈對數和指數變化.常用的脈動風功率譜有Kaimal譜、Von Karman譜等.在模擬空間網格上每點的脈動風速時間歷程時，必須考慮到這些時間歷程是相干的，且這種相關性隨著距離的增大而減小，在高頻時比低頻小，用相干函數來描述這種關系.

通過NREL開發的前處理程序——TurbSim生成風速場.該程序能夠根據不同的參數生成三維、隨機的湍流風速場.選用指數模型考慮風剪切效應：

(15)

在TurbSim中選用IEC61400-1中的Kaimal譜[16]，計算公式如下：

(16)

根據IEC61400-1標準，模擬正常湍流風況；選用該標準的相干方程考慮縱向風速u的相關性：

(17)

但該標準沒有定義橫向風速v和垂向風速w的相干模型，故在選用了IEC標準的前提下，TurbSim不考慮橫向和垂向風速的相干性.

1.4.2 隨機風荷載計算

基于隨機三維全域風場計算軸向和切向誘導因子在葉片上的分布，從而得到誘導速度，進而求解葉片所經歷的風速及攻角.圖4為風輪平面內的速度.

常用計算誘導速度的理論有葉素動量理論(BEM)和動態尾流理論(GDW).動態尾流理論是基于拉普拉斯方程的勢流解，相比于葉素動量理論，能夠考慮轉子平面上的更廣義的壓力分布，且能模擬流經葉片的風速變化時葉片受到的氣動載荷變化的時間滯后效應.在此基礎上，利用由風洞實驗或者CFD計算獲得的翼型空氣動力學參數，包括升、阻力系數與攻角的關系和動態失速參數，最終計算作用在葉輪上的氣動荷載.局部的升力和阻力按如下公式計算：

(18)

(19)

采用NREL開發的程序——FAST計算葉輪上的氣動荷載，在該程序的氣動荷載計算子模塊AeroDyn對應的輸入參數中，選用Peters與He[17]提出的動態尾流模型來計算誘導速度.在FAST的輸入文件中，需要配置好所用到的翼型數據文件，每個翼型數據文件中包含了在±180°范圍內的攻角所對應的升力系數 、阻力系數和變槳系數.表1為DU系列和NACA系列翼型族的葉素參數.由于風剪切、偏航和湍流等因素，葉片將不可避免地出現動態失速效應，故在AeroDyn中采用Leishman和Beddoes[18]提出的半經驗模型動態失速模型，并在翼型數據文件中配置好動態失速參數.最后輸出3個方向的剪力和彎矩時程，將風荷載作用在風機塔頂節點.圖5和圖6分別為較大風浪時塔頂x方向剪力和繞y軸彎矩的時程圖.

圖4 風輪平面內的速度

表1 葉素參數

時間/s

時間/s

2 算 例

2.1 單樁式基礎海上風機

本文使用美國可再生能源實驗室(NREL)建立的5 MW海上單樁風機模型進行分析，該模型整體參數見表2[19].樁頂高出海底30 m，塔架安裝在樁頂.塔架呈錐形，頂部直徑3.87 m，壁厚0.019 m，底部直徑6 m，壁厚0.027 m.單樁呈圓柱形，直徑6 m，壁厚0.06 m.

表2 風機總體參數

2.2 計算工況

本文選取5 WM風機模型，以水深、風浪大小、波高影響為變量，計算了24種組合工況.如表3所示.

表3 計算工況

其中“風浪大小”可以用輪轂風速控制，選取不同的輪轂風速，然后可以根據式(16)得到海上19.5 m處的風速，繼而代入P-M譜公式(5)就可以得到波浪參數.3個取值分別代表了“較小風浪”、“中等風浪”和“較大風浪”.

為了計算方便，工程上通常將波浪力從海底計算到靜水面.但是對于淺水區非線性波浪, 這種取法可能產生較大的誤差，因此也有學者將其從海底計算到自由波面.為了考慮波面變化對流固耦合問題的影響，將其作為一個自變量，不考慮其影響時，將任意時刻波浪力從海底計算到靜水面，考慮其影響時，將任意時刻波浪力從海底計算到自由波面.本文是按照線性外延法考慮自由液面的影響，即把線性波理論延伸到自由液面.波浪質點速度的修正采用Chakrabarti[20]提出的等效水深方法，具體計算公式如下：

(20)

式中H為波高；T為波浪周期；d為水深；k為波數；z為質點z方向的坐標；η為波高.

2.3 計算結果

為進行分析，定義相對誤差δ的計算公式：

(21)

式中Max可為波浪力、結構位移或加速度最大值.

對于波浪力，在所有工況的計算結果中，考慮耦合與不考慮耦合計算結果的差別較大.圖7為工況“d=35 m，較大風浪，考慮波面影響”時，結構上z=30 m處的質點所受波浪力時程圖.當發生最大值時，計算兩種情況的相對誤差值，所有工況中的最大相對誤差達到41.7%.將圖7所示波浪力進行頻譜分析發現考慮流固耦合時，波浪力出現了高頻的成分，而不考慮耦合時波浪力沒有高頻成分，如圖8所示.

為研究高頻成分出現的原因，對結構進行模態分析.考慮到結構的截面是xy平面的圓形截面，其在x和y方向振動的模態相同.本文研究的是結構x方向運動，只需要研究第1,4階模態，第2,3階是與1,4階對應的y方向振動，如圖9所示.可以發現波浪力中的高頻成分與結構自振頻率基本吻合.

t/s

頻率/Hz

圖9 第1，4階振型圖

為了研究波浪力的差別與波面因素的關系，計算所有質點的總波浪力，在結果中發現波面因素對波浪力的影響較大.比如在工況“水深20 m，中等風浪”中，考慮波面影響時的波浪力峰值較不考慮時大47.1%，所以計算波浪力時，建議考慮波面因素的影響，并且波浪力峰值的差別隨水深加大而變大.圖10為考慮波面影響時，水面以下5 m處質點波浪力峰值相對誤差與水深的關系圖，圖11為波浪合力的峰值相對誤差與水深的關系圖.從圖10和圖11能看出，在“風浪較小”時的波浪力相對誤差較大，這是由于考慮了風機變漿的緣故.在一定范圍內，隨著風速的增大，雖然波浪荷載變大但風荷載反而變小，導致“風浪較小”的情況下波浪力的相對誤差較大，由此推斷風浪大小對波浪力有一定的影響.

水深/m

水深/m

對于動力響應的差別，在所有工況的計算結果中，考慮耦合與不考慮耦合計算的位移、速度的差別較小，而加速度差別較大.所有工況中，塔頂位移最大相對誤差為7.6%，塔頂速度最大相對誤差為4.6%，而塔頂加速度最大相對誤差為35.5%.圖12～14為工況“d=35 m，較大風浪，考慮波面影響”時，塔頂處的動力響應時程圖.并且最大加速度相對誤差值跟風浪和水深也有一定的關系.圖15為考慮波面影響時，最大加速度相對誤差與風浪、水深的關系圖.

時間/s

時間/s

時間/m

水深/m

3 結 論

本文研究了固定式海上風機的流固耦合問題，分別建立了海上風機塔架考慮與不考慮流固耦合效應時的動力方程及數值模擬方法.計算結果表明：

1)考慮流固耦合問題時比不考慮流固耦合時計算的風機波浪力峰值較大，并且從頻譜分析可以發現：考慮流固耦合時波浪力會出現高頻的成分，波浪力高頻的成分可能引起結構高頻振動.并且，波浪力峰值差別的大小隨水深加大而變大，且與風浪大小有一定關系.此外，計算波浪力時應該考慮波面的影響，并需要綜合考慮風浪大小.

2)對比考慮和不考慮流固耦合對風機塔架動力響應的影響時發現，兩者情況下獲得的位移和速度的差別較小，加速度差別較大，并且加速度的差別與水深和風浪大小有一定關系.

加速度是描述風機結構振動的重要指標，振動過大會對風機正常運行造成嚴重影響.基于本研究結果，對于固定式海上風機，建議在水深和風浪較大的情況下分析風機塔架的動力響應時應考慮流固耦合問題，且同時考慮自由波面的影響，從而減小計算結果的誤差.這對風機的設計及工作性能評估均具有現實意義.

[1] MORISON J, JOHNSON J, SCHAAF S. The force exerted by surface waves on piles[J]. Journal of Petroleum Technology, 1950,2(5): 149-154.

[2] SARPKAYA T, ISAACSON M. Mechanics of wave forces on offshore structures[M]. New York: Van Nostrand Reinhold Company,1981:25-32.

[3] FISH P, DEAN R, HEAF N. Fluid-structure interaction in Morison's equation for the design of offshore structures[J]. Engineering Structures, 1980,2(1):15-26.

[4] 王元戰, 龍俞辰, 王朝陽. 考慮流固耦合影響的樁基波浪力簡化計算方法[J]. 水道港口, 2014，35(2):93-98.

WANG Yuan-zhan, LONG Yu-chen, WANG Zhao-yang. Simplified calculation method for wave force of piles under effect of fluid-structure interaction[J]. Journal of Waterway and Harbor, 2014，35(2):93-98.(In Chinese)

[5] 徐亞洲, 李杰. 結構流固耦合波浪動力響應的一階預估校正法[J]. 船舶力學, 2012，16(7):781-786.

XU Ya-zhou, LI Jie. First order prediction correction method for dynamic wave responses considering fluid-structure-interaction[J]. Journal of Ship Mechanics, 2012，16(7):781-786. (In Chinese)

[6] CHANDRASEKARAN S, JAIN A. Triangular configuration tension leg platform behaviour under random sea wave loads[J]. Ocean Engineering, 2002,29(15):1895-1928.

[7] DNV-OS-J101 Design of offshore wind turbine structures[S]. Oslo,Norway:Det Norske Veritas,2010:48-50.

[8] SHIRZADEH R, DEVRIENDT C, BIDAKHVIDI MA,etal. Experimental and computational damping estimation of an offshore wind turbine on a monopile foundation[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2013,120:96-106.

[9] AGARWAL P, MANUEL L. Incorporating irregular nonlinear waves in coupled simulation and reliability studies of offshore wind turbines[J]. Applied Ocean Research, 2011,33(3):215-227.

[10]李德源, 劉勝祥, 張湘偉. 海上風力機塔架在風波聯合作用下的動力響應數值分析[J]. 機械工程學報, 2009,45(12):46-52.

LI De-yuan, LIU Sheng-xiang, ZHANG Xiang-wei. Dynamical response numerical analysis of the offshore wind turbine tower under combined action of wind and wave[J]. Journal of Mechanical Engineering, 2009,45(12):46-52. (In Chinese)

[11]張學志, 黃維平, 李華軍. 考慮流固耦合時的海洋平臺結構非線性動力分析[J]. 中國海洋大學學報 :自然科學版, 2005，35(5):823-826.

ZHANG Xue-zhi, HUANG Wei-ping, LI Hua-jun. Nonlinear dynamic analysis of offshore platform considering fluid-structure interaction[J]. Periodical of Ocean University of China, 2005,35(5):823-826. (In Chinese)

[12]GL2010 Guideline for the certification of offshore wind turbines[S]. Uetersen,Germany: Germanischer Lloyd Wind Energie GmbH,2010: 51-58.

[13]徐亞洲, 李杰. 近海風力發電高塔波浪隨機動力響應分析[J]. 振動工程學報, 2011，24(3):315-322.

XU Ya-zhou, LI Jie. Dynamic responses of offshore wind turbines subjected to random waves[J]. Journal of Vibration Engineering, 2011，24(3):315-322. (In Chinese)

[14]PIERSON W, MOSKOWITZ L. A proposed spectral form for fully developed wind seas based on the similarity theory of SA Kitaigorodskii[J]. Journal of Geophysical Research, 1964,69(24):5181-5190.

[15]余聿修. 隨機波浪及其工程應用[M]. 大連: 大連理工大學出版, 2003:132-136.

YU Yu-xiu. Random wave and its applications to engineering[M]. Dalian:Dalian University of Technology Press, 2003:132-136. (In Chinese)

[16]IEC 61400-1 Wind turbines-Part 1: Design requirements[S]. Geneva,Switzerland: International Electrotechnical Commission, 2005:70-71.

[17]PETERS D, HE C. Correlation of measured induced velocities with a finite-state wake model[J]. Journal of the American Helicopter Society, 1991,36(3):59-70.

[18]LEISHMAN J, BEDDOES T. A semi-empirical model for dynamic stall[J]. Journal of the American Helicopter Society, 1989,34(3):3-17.

[19]JONKMAN J, BUTTERFIELD S, MUSIAL W,etal. Definition of a 5-MW reference wind turbine for offshore system development[R]. Colorado,US:National Renewable Energy Laboratory,2009:2-16.

[20]CHAKRABARTI S K. Hydrodynamic of offshore structures[M]. Berlin: Springer, 1987:45-50.

Numerical Simulation of Dynamic Response for Offshore Wind Turbines including Fluid-Structure Interaction

DENG Lu1?, XIAO Zhi-ying1, HUANG Ming-xi1, SONG Xiao-ping2, WU Hai-tao2

(1.College of Civil Engineering, Hunan Univ, Changsha, Hunan 410082, China;2.XEMC Windpower Co Ltd, Xiangtan,Hunan 411102, China)

This paper investigated the issue of fluid-structure interaction (FSI) for fixed offshore wind turbines. The model for a 5 MW monopile wind turbine built by the National Renewable Energy Laboratory (NREL) was selected as an example. The stochastic inflow turbulence code written by NREL-TurbSim, the Kaimal spectrum and exponential wind profile were used to generate stochastic and full field turbulent wind. The aeroelastic code FAST developed by NREL was then used to calculate the aerodynamic load on the wind turbine. The time histories of shear force and bending moment at the top of the tower were obtained. The random wave was simulated using the Pierson-Moskowits spectrum through the harmony superposition method. Two cases were considered in this study, namely, the case that considered the effect of FSI and the case in which the effect of FSI was neglected. Kinematic equations and numerical models under the two cases were developed respectively, and the comparison between the two cases was made to investigate the effect of FSI. The results indicate that distinct deviations affected by the intensity of the wind and wave exists in both wave force and dynamic responses in the two situations.

offshore wind turbine; random wind and wave; Fluid-Structure interaction; Morison equation

1674-2974(2015)07-0001-08

2014-10-12

國家863計劃資助項目(2013AA050603)；湖湘青年創新創業平臺

鄧 露(1984-)，男，湖南婁底人，湖南大學教授，博士

?通訊聯系人，E-mail:denglu@hnu.edu.cn

TV139.2

A