深水SPAR風機系統(tǒng)全耦合動力響應分析研究

閆發(fā)鎖，門驥遠，彭 成

（1.哈爾濱工程大學 船舶工程學院，哈爾濱 150001；2.美國德州農(nóng)工大學 土木工程系，大學城 77843）

深水SPAR風機系統(tǒng)全耦合動力響應分析研究

閆發(fā)鎖1，門驥遠1，彭 成2

（1.哈爾濱工程大學 船舶工程學院，哈爾濱 150001；2.美國德州農(nóng)工大學 土木工程系，大學城 77843）

文章采用聯(lián)合開發(fā)的計算程序對深水SPAR風機的浮體、錨泊和風機各子系統(tǒng)進行了水—氣動力的全耦合數(shù)值分析，研究了深水浮式風機系統(tǒng)的動力響應特點。浮體水動力計算采用基于二階精度的混合波浪模型（Hybrid Wave Model）的MORISON公式，錨泊系統(tǒng)采用細長桿理論通過非線性有限元方法實現(xiàn)，風機系統(tǒng)的空氣動力分析采用基于多體氣動彈性理論的FAST模塊。以浮體控制方程為主體，通過模塊間的載荷與位移傳遞在每個時間步上迭代求解，形成完全耦合的時域分析方法。通過對NREL的5MW SPAR風機系統(tǒng)在隨機海況下的水動力響應分析，驗證了該方法的有效性，并分析了浮式風機子系統(tǒng)間的混合動力作用。

浮式風機；動力響應；混合波浪模型；系統(tǒng)耦合

0 引 言

海洋風能儲量豐富，在倡導綠色環(huán)保、節(jié)能減排的可持續(xù)發(fā)展的全球背景下，海上風電的開發(fā)利用是資源環(huán)境與社會可持續(xù)發(fā)展模式的典型示范，已成為國際可再生能源開發(fā)的主力方向。由于技術和成本的限制，直至上世紀90年代海上風電場的建設才取得了長足的發(fā)展。海上漂浮風機系統(tǒng)由葉輪、塔架、浮體和錨泊系統(tǒng)構成，是深水風電開發(fā)的基礎裝備。目前在深水浮式風機技術方面，歐美和日本等發(fā)達國家已積累了豐富的理論成果和實踐經(jīng)驗，完成了多個概念設計和試驗樣機，其中影響較大的有SWAY、HYWIND、WindSea AS、Blue H和WindFloat等（圖1）。

圖1 幾種海上浮式風機系統(tǒng)（由左至右分別為：SWAY-2011、HYWIND-2009、WindSea AS-2008、Blue H-2008和WindFloat-2008）Fig.1 Some types of offshore floatingwind turbine systems

與固定式風力機相比，漂浮式風力機浮式基礎在海洋環(huán)境載荷作用下有一定幅度的運動。風力機位于距海面近百米的高空，浮式基礎的微幅運動即可造成風力機姿態(tài)較顯著的變化，對風力機氣動載荷和功率輸出會產(chǎn)生顯著影響。因此，研究浮式風機系統(tǒng)的運動響應對發(fā)揮海上風力機的性能特點，進行合理的設計優(yōu)化具有重要的指導意義。

風機葉輪的氣動載荷是浮式風機系統(tǒng)重要的環(huán)境激勵。風機氣動力性能的分析目前主要有兩種方法：傳統(tǒng)的動量葉素理論（BEM）和基于CFD的渦面元方法[1]。BEM方法基于二維假設，理論模型簡潔直觀、計算快速，但葉素單元的受力計算需要合理的試驗或經(jīng)驗系數(shù)和端部修正。但從實用意義考慮，針對其局限性已發(fā)展了豐富的修正處理方法，BEM仍然是風機葉片設計載荷的主要計算手段。Gutierrez等[2]采用基于BEM的FAST軟件開發(fā)了浮式風力機的優(yōu)化設計工具FASTLognoter，具備風力機的參數(shù)設計能力。渦面元法考慮了葉片流場旋渦尾流特性，理論模型精確，但數(shù)值實現(xiàn)技術難度高，計算量大，目前還主要處于研究階段。此外，氣動力分析方法還要結合風機控制系統(tǒng)的控制模式和葉輪、塔架等的結構特征，才能全面考慮系統(tǒng)的動力性質。基于渦面元法，Sebastian等[3]發(fā)展了自由渦放尾流（FVM）方法進行了風力機的氣動響應分析。該方法相對BEM在葉片尾流演變和氣動載荷方面體現(xiàn)了較高的精度，但計算資源耗費較大。

借助于大型海洋平臺的設計分析技術，海上風機系統(tǒng)中浮體的選型與優(yōu)化、錨泊系統(tǒng)設計以及二者的耦合水動力響應分析已具備一定的基礎。Jonkman等[4-5]基于FAST對已有的半潛式、張力腿式和SPAR式三種類型風機浮體的水動力進行了分析。唐友剛等[6]針對一座5MW海上風機設計了SPAR浮體和三立柱—撐桿形式的半潛浮體，并進行了性能和結構分析。Sethuraman等[7]對一個階梯型SPAR平臺風機的水動力響應進行了模型試驗和數(shù)值預報。Karimirad等[8]采用兩種分析代碼針對一個單柱張力筋系泊的浮式風機系統(tǒng)進行了平臺運動和系泊張力的水彈性時域模擬，并對系統(tǒng)耦合進行了模型試驗研究。

綜上所述，海上風機各組成系統(tǒng)綜合了水、空氣結構之間的多域動力耦合分析技術。以浮式平臺的相關理論為基礎，利用現(xiàn)有浮體水動力分析方法結合風力機械的氣動力載荷模塊，進行浮式風機系統(tǒng)的動力響應分析成為合理可行的研究方式。目前已有的分析工具包括業(yè)界具有較大影響的DeepC、OrcaFlex、Charm3D等水動力分析程序和FAST、GH Bladed、TurbSim等風機氣動分析模塊[9-11]。針對風機—浮體—系泊三者之間的氣動力與水動力的多域耦合，目前的做法需要采用各自獨立的分析模塊在時域上迭代耦合，分析結果和計算效率也有待檢驗和提高。本文將水動力耦合分析程序與現(xiàn)有的風機氣動力程序結合，完成了適用于深水浮式風機系統(tǒng)的全耦合時域動力響應分析。浮體水動力的分析采用具有二階精度的混合波浪模型（Hybrid Wave Model）與基于細長體假設的MORISON方法；錨泊系統(tǒng)采用細長桿理論通過非線性有限元方法實現(xiàn)；風機系統(tǒng)的空氣動力分析采用基于BEM和多體氣動彈性理論的FAST模塊。通過對NREL 5MW Hywind OC3風機動力響應的數(shù)值模擬與現(xiàn)有結果的比較，分析了浮式風機子系統(tǒng)間的混合動力作用。

1 數(shù)值模型

1.1 風機氣動載荷分析

系統(tǒng)中風力機的氣動載荷采用了FAST[12]（Fatigue,Aerodynamics,Structures,Turbulence）。FAST是由美國國家可再生能源實驗室（NREL）開發(fā)的一個開放資源程序，用以模擬傳統(tǒng)兩葉或三葉水平軸風機的氣動彈性載荷。FAST中計算風載荷的主要模塊是AeroDyn，其主要采用了經(jīng)典的動量葉素理論[13]。

基于動量葉素理論，風輪氣動性能的計算

式中：Q和T分別為作用在風輪上的扭矩和盤面推力；ρ為空氣密度；CL和CD分別為升力和阻力系數(shù)；φ為入流角；RH和R分別為輪轂及葉片半徑；B為葉片數(shù)目；c為葉片單元弦長；r為葉片單元的展向半徑；Vr為葉片單元處的氣流相對速度，其中U為遠方來流風速，a和a′分別

r為軸向和周向誘導因子；F為葉片端部尾流的修正系數(shù)，源于Prandtl的氣動理論方法[13]，其中

1.2 浮體動力分析

采用圖2（a）所示坐標系，對某一剛體其上任意一點的6自由度非線性運動方程可以用（3）式和（4）式表示。

圖2 分析所采用的坐標系統(tǒng)Fig.2 Coordinate systems for theorymodels

式中：上標‘t’表示轉換矩陣，m為結構質量，ξ=（ξ1,ξ2,ξ3）t為坐標系（即大地坐標系）o點坐標值，=d2ξ/ d t2為線性加速度，ω= （ω1,ω2,ω3）t為角速度，rg=（xg,yg,zg）t為中心坐標值，I0為剛體的慣性矩。此外，u、rg和I在隨體坐標系即坐標系中定義。要注意到和M分別是作用于剛體上的總力和總力矩，后者是oxyz

00在坐標系中定義的，并且T是（oxyz）坐標系到（OXYZ）坐標系的轉換矩陣，其為正交矩陣，正交變換與自身相同 （Tt=T-1）。

當柱體直徑相對波長較小時，選取合理的阻力系數(shù)和附加質量系數(shù)，莫里森方程可較準確地用于計算柱體上的波浪載荷。沿著與柱體軸向垂直的方向n作用于柱體單元上的受力（慣性力和阻力）可以用下式表示：

式中：下標‘n’表示法向分量，Dj為等效柱體直徑，zj為第j個單元的濕長，ρf為環(huán)境流體的密度，Cm為附加質量的系數(shù)，CD為阻力系數(shù)。vr和ar分別為環(huán)境流體相對于平臺結構的速度和加速度，且an為柱體加速度。作用于某一穿過自由液面的柱體的軸向水動力，包括作用于單元底部面積SBj上的動壓力，阻力和附加質量力。軸向力可以表示為：

式中：下表‘t’表示軸向分量，φ（1）和φ（2）為一階和二階入射波勢，CDt和Cmt為垂向阻力系數(shù)及附加質量系數(shù)，且vrt和art是柱體相對于靠近柱體底部附近的環(huán)境流體相對速度和加速度的軸向分量。

1.3 系泊系統(tǒng)分析

當沒有額外轉矩和力矩作用于系泊線/立管，在圖2（b）坐標系下運動的方程可以表示為

式中：r（s,t）為系泊線/立管的位移矢量，是弧長s和時間t的函數(shù)；B表示彎曲剛度；q表示單位長度外力，其可通過莫里森方程算得；r為單位長度質量；λ=T-Bκ2為一個標量變量，其中κ為系泊線/立管的曲率，而T為局部張力。“’”表示對s的導數(shù)，“·”表示對時間t的導數(shù)=?2t/?t2為局部加速度。此外，對于不計軸向拉伸的纜索，r必須滿足幾何約束方程

作用于桿上的外力包括重力、流體靜水力和流體動水力。作用于桿上的流體動水力包括附加質量力、阻力和Froude-Krylov力。方程（7）可以通過非線性FEM法數(shù)值求解。使用Spar及其系泊系統(tǒng)間的配置條件，柱體及系泊系統(tǒng)的動力方程可以通過Newmark-β方法在時域上同時求解。

2 風機系統(tǒng)參數(shù)

為了驗證浮體與系泊系統(tǒng)間的耦合作用，本文選取NREL 5MW OC3 SPAR[14]風機進行了系統(tǒng)動力響應分析，并將計算結果與現(xiàn)有的三維勢流理論軟件分析結果進行了比較。

風機系統(tǒng)包括上部傳統(tǒng)的三葉片變螺距風力機及塔架與下部Hywind-OC3 Spar平臺及錨泊兩部分。Hywind-OC3 Spar風機設計水深320m，系泊系統(tǒng)采用3根系泊線，它們在靜水面投影之間的夾角為120°。導纜器位置在靜水面以下70m處，導纜口距平臺柱體中心線之間的半徑為5.2m。海底錨點分布在系統(tǒng)靜平衡時距SPAR中心線半徑為853.9 m的圓周上。每根系泊線長度為902.2 m，直徑0.09m。對系泊系統(tǒng)的模擬做了以下簡化：系泊線與平臺的三個連接點用鉸鏈連接代替，艏搖方向加入彈簧以模擬原結構剛度；實際上多成分的系泊線近似為均勻濕重（單位長度）和剛度的錨纜。風機結構的主體尺度與系泊系統(tǒng)的布置如圖3所示，系統(tǒng)總體的結構參數(shù)見表1。

圖3 Hywind-OC3 5MW浮式風機系統(tǒng)組成Fig.3 Sketch of Hywind-OC3 5MW floating wind turbine

表1 Hyw ind-OC3 5MW浮式風機系統(tǒng)結構參數(shù)Tab.1 Characteristics of Hyw ind-OC3 5MW w ind turbine system

3 系統(tǒng)頻域響應分析

為例驗證各子系統(tǒng)理論方法及數(shù)值實現(xiàn)，首先對風機浮體進行頻域水動力分析，得到自由浮體的線性頻響函數(shù)；之后，考慮錨泊系統(tǒng)的影響，得到系統(tǒng)的水動力耦合頻響函數(shù)。二者在分析時假設風機塔架和葉片是靜止的，未計及上部風機的氣動力載荷，但考慮了風機結構的質量及分布。圖4為浮體系統(tǒng)縱蕩、垂蕩、縱搖和首搖4個自由度的波激運動頻響函數(shù)。浮體在自由狀態(tài)下，垂蕩運動的固有周期為31.4 s.，縱搖/橫搖為31 s。

圖4 浮體系統(tǒng)的運動頻響函數(shù)Fig.4 Motions frequency response function of floating body including turbine structure

SPAR平臺系統(tǒng)通常使用半張緊式的錨泊系統(tǒng)，錨泊系統(tǒng)的影響主要體現(xiàn)在平動自由度上。基于自由衰減法進行風機系統(tǒng)的時域耦合分析，表2列出了自由衰減模擬的前五個周期的平均值。由表2可見，有了錨泊的約束，系統(tǒng)平面運動固有周期明顯減少，垂向及相關運動變化不顯著。為了便于驗證比較，表中同時給出了文獻[9]水動力耦合計算的結果。二者在橫蕩、縱蕩和垂蕩運動的周期幾乎一致，而在縱搖、橫搖和艏搖有較小差別。本文方法計算時較準確地考慮了錨泊系統(tǒng)的影響，與自由浮體相比，垂蕩、縱搖固有周期有較小降低，縱搖、橫搖因錨泊布置的差別稍有不同，艏搖周期相對有稍大差別。

表2 平臺系統(tǒng)的運動固有周期Tab.2 Natural frequencies of platform motions

4 全耦合數(shù)值分析與比較

應用數(shù)值技術將浮式風機系統(tǒng)的水動力模塊與氣動彈性模塊耦合，得到風機系統(tǒng)結構中心處的運動響應及譜密度，如圖5所示。風機作業(yè)海域水深320m，分析環(huán)境條件選用JONSWAP表示的隨機不規(guī)則海況，有義波高5.0m，峰值周期8.69 s，參數(shù)γ取值3.6，環(huán)境風速13m/s。風、浪的傳播方向沿x軸正向。因為計算步長小（0.02 s）,且錨泊系統(tǒng)應用細長桿理論需要大量的迭代計算，計算耗時較長。為了合理控制時間和充分利用計算資源，分析采用了分段進行的方式。本例對于3 600 s的計算分三段進行，每段計算1 300 s。為了避免每個分段起始部分重復，各分段之間計算預留了150 s的共用重合區(qū)，該區(qū)間內應用適當方法進行了處理。為了清晰顯示，圖中的運動時歷的時間為1 300 s，對應的功率譜密度（PSD）的統(tǒng)計是基于3 600 s的計算結果。

圖5 風機系統(tǒng)全耦合運動響應與能量譜密度Fig.5 Motion responses and PSD by the full coupledmethod

由于全耦合分析計入了風力機葉片的氣動彈性載荷，系統(tǒng)在風浪方向的運動（縱蕩）具有較大的運動幅度和平均偏移。根據(jù)六自由度運動的時域統(tǒng)計，縱蕩運動平均值為9.369 m，最大值為17.708 m，最小值為4.125m，標準差為1.747m。因為在耦合分析中，由于系統(tǒng)上部風力機的葉片的旋轉，再加上SPAR主體縱搖速度影響葉片迎面風速，系統(tǒng)縱向受力幅值相比忽略風力作用（圖5）時變化大，響應的運動幅值也較大。由此表明，風力機縱向載荷對水動力載荷在頻域方面的影響是較大的。本文方法要考慮圓柱體繞流旋渦脫落導致的橫向力，而且當風機工作時也會產(chǎn)生一個相當小的橫向力。因此，盡管風和浪在x軸方向，主體仍具有一定的橫蕩和橫搖。耦合分析得到的艏搖運動比較小。這實際上當主體迫使風機做艏搖偏轉時，風機便可以通過監(jiān)控機械來實現(xiàn)控制，這種作用于風機上主體首搖運動的反饋在分析時應該被考慮在內。

采用細長桿理論與浮體運動的耦合求解，得到的錨泊系統(tǒng)的受力如圖6所示。圖6（a）為三根系泊線的頂部張力，圖6（b）為底部水平張力。每根系泊線的頂部張力平均比底端張力大約200 kN。由于對稱性，理論上#1系泊線的張力應該與#2的張力相同。實際數(shù)值計算時，兩根系泊線是獨立迭代計算的，并且由于系統(tǒng)橫向載荷和運動的影響，二者的張力在時域上有較小的局部差別，但總體的吻合度很好。#3的布置與#1和#2相反，其張力大小變化與后兩者呈相反趨勢，結果也得到了很好的體現(xiàn)。

圖6 錨泊線張力Fig.6 Tensions inmooring lines

全耦合在上部風力機氣動載荷計算時，考慮了浮體SPAR的運動狀態(tài)。已有的部分研究在計算風力時假設浮體SPAR是固定的，即認為浮體的姿態(tài)和運動對風力機的影響較小。通過將本文全耦合方法與后者假設情況下，在上述隨機海況下進行計算分析對比表明，系統(tǒng)的垂蕩響應具有較明顯的差別（圖7）。計算氣動力載荷時不考慮浮體的運動狀態(tài)，系統(tǒng)的最大垂蕩響應接近全耦合計算時的2倍。因為系統(tǒng)的縱搖的自然周期（30.7 s）非常接近于升沉運動固有周期（30.9 s），風機所產(chǎn)生的外部縱搖力矩與垂向波浪力的耦合加劇了垂向運動。而在全耦合分析時考慮了SPAR浮體的每時刻的運動和空間姿態(tài)，垂蕩和縱搖運動的固有頻率產(chǎn)生了偏移，二者的諧振效應明顯減弱。

圖7 全耦合與近似耦合垂蕩響應的比較Fig.7 Comparison of heavemotion between full coupling and simple coupling

5 結 論

本文研究了海上浮式風機系統(tǒng)全時域耦合分析方法，將水動力載荷模塊與空氣動力載荷模塊進行聯(lián)合求解，并對NRELHywind-OC3 5MW風機系統(tǒng)的運動響應進行了全耦合分析。計算結果表明，全耦合方法有效地對風機系統(tǒng)的自由浮體和具有張緊式系泊狀態(tài)下的運動RAO進行分析。在隨機海況下，由于全耦合分析計入了風力機葉片的氣動彈性載荷，系統(tǒng)在風浪方向的運動（縱蕩）具有較大的運動幅度和平均偏移。由于全耦合考慮了風力機迎風速度的變化，全耦合時系統(tǒng)的動力響應與不考慮浮體運動狀態(tài)的近似耦合相比，某些自由度的運動可能會有較大差別。針對本文算例所采用SPAR類型平臺，因為系統(tǒng)的垂蕩和縱搖固有頻率相近，兩種方法得到的垂蕩運動在時域和頻域結果有顯著的差別。本文的研究表明，針對浮式風機系統(tǒng)發(fā)展全耦合分析方法是必要的。

[1]Hansen M O L,Sorensen JN,Voutsinasb S,Sorensen N,Madsen H Aa.State of the art in wind turbine aerodynamics and aeroelasticity[J].Progress in Aerospace Sciences,2006,42:285-330.

[2]Gutierrez JE,Zamora B,García J,Peyrau M R.Tool development based on FAST for performing design optimization of offshorewind turbines:FASTLognoter[J].Renewable Energy,2013,55:69-78.

[3]Sebastian T,Lackner M A.Developmentof a free vortex wakemethod code for offshore floating wind Turbines[J].Renewable Energy,2012,46:269-275.

[4]Jonkman JM,Butterfield S,MusialW,Scott G.Definition of a 5-MW reference wind turbine for offshore system development[R].NREL/TP-500-38060,National Renewable Energy Laboratory.Colorado,USA,2009.

[5]Jonkman JM,Matha D.Dynamics of offshore floating wind turbines-analysis of three concepts[J].Wind Energy,2011, 14:557-569.

[6]唐友剛,王 涵,陶海成,劉中柏.海上風機半潛型浮式基礎結構設計及整體強度分析[J].中國造船,2013,54(3):85-93.Tang Yougang,Wang Han,Tao Haicheng,Liu Zhongbai.Structure design and global strength analysis for semi-submersible floating foundation of offshorewind turbine[J].Shipbuilding of China,2013,54(3):85-93.

[7]Sethuraman L,Venugopal V.Hydrodynamic response of a stepped-spar floating wind turbine:Numericalmodelling and tank testing[J].Renewable Energy,2013,52:160-174.

[8]Karimirad M,Meissonnier Q,Gao Zhen,Moan T.Hydroelastic code-to-code comparison for a tension leg spar-type floating wind turbine[J].Marine Structures,2011,24:412-435.

[9]Bae Y H,Kim M H,Im SW,et al.Aero-Elastic-Control-Floater-Mooring coupled dynamic analysis of floating offshore wind turbines[C]//International Offshore and Polar Engineering Conference,June 19-24,2011.Maui,Hawaii,USA,2011.

[10]Masciola M,Robertson A,Jonkman JM,Driscoll F.Investigation of a FAST-OrcaFlex couplingmodule for integrating turbine andmooring dynamics of offshore floatingwind turbines[C]//International Conference on OffshoreWind Energy and Ocean Energy.Beijing,China,2011.

[11]Karimirad M.Dynamic response of floatingwind turbine[J].Mechanical Engineering,2010,17(2):146-156.

[12]Jonkman JM.,Jr,BuhlM L.FASTUser’s Guid[R].NREL/EL-500-38230,2005.

[13]Moriarty P J,Colorado G,Craig Hansen A.AeroDyn Theory Manual[R].NREL/EL-500-36881.National Renewable Energy Laboratory.Colorado,USA,2005.

[14]Jonkman J.Definition of the floating system for phase IV of OC3[R].NREL/TP-500-47535.National Renewable Energy Laboratory.Colorado,USA,2010.

Study on full coup ling dynam ic responses of a deepwater SPAR w ind turbine system

YAN Fa-suo1,MEN Ji-yuan1,PENG Cheng2
(1.Harbin Engineering University,Harbin 150001,China;2.Texas A&M University,College Station 77843,USA)

Themethod to perform coupling analysis for offshore floating wind turbines(OFWT)is developted through integrating hydrodynamic and aerodynamicmodules in time domain.Morisonmethod is used for hydrodynamic computation of floating body and itsmooring system,in which relative velocity between structure elements and waves is implemented by the Hybrid Wave Modelwith second order accuracy.Slender rods theories are applied on themooring systems,and the aerodynamics load induced by wind turbine is simulated by NREL’s code-FAST.Loads and displacements are transferred between the submodules based mainly on floating body control equations in every time step by Newmark-βmethod.Motion response of a 5MW 3 blades spar type OFWT is predicted with and without FAST to validate the combined program.A comparison with results from available 3D linear potential flow method in a random sea condition shows that the code is capable of hydro-aero dynamic analysis for OFWT.

OFWT;dynamic responses;hybrid wavemodel;subsystems coupling analysis

P752 TK89

：Adoi:10.3969/j.issn.1007-7294.2017.02.005

2016-09-30

國家自然科學基金項目（51379051，50909022）；上海交通大學海洋工程國家重點實驗室開放基金（1208）；國家國際科技合作專向項目（2013DFE73060）

閆發(fā)鎖（1977-），男，博士，副教授，E-mail:yanfasuo@hrbeu.edu.cn；

門驥遠（1990-），男，博士研究生。

1007-7294（2017）02-159-09