海上浮式風機運動響應的時域耦合計算方法

金飛，滕斌

（大連理工大學海岸和近海工程國家重點實驗室，遼寧大連 116024）

海上浮式風機運動響應的時域耦合計算方法

金飛，滕斌

（大連理工大學海岸和近海工程國家重點實驗室，遼寧大連 116024）

建立了一種用于計算海上浮式風機運動響應的時域耦合方法，該算法主要由氣動力模塊、水動力模塊、系泊模塊和系統(tǒng)運動模塊構(gòu)成。氣動力模塊采用葉素動量法；水動力模塊采用一階勢流理論，通過邊界元法求解；系泊模塊采用懸鏈線模型，用Chebyshev多項式進行擬合計算。對于系統(tǒng)運動模塊，采用Runge-Kutta法求解。對OC3-Hywind spar風機進行了建模，對各模塊及耦合模型的進行對比研究，驗證了此耦合計算方法的準確性。并利用該方法計算和分析了此風機的運動響應及其對風機功率的影響。

浮式風機；時域耦合方法；葉素動量法；邊界元法；懸鏈線模型

風力發(fā)電作為一種利用可再生、清潔能源的技術(shù)，已經(jīng)得到國內(nèi)外廣泛研究和使用[1]。目前風機的主要形式可分為陸上風機和海上風機2種。而海上風機又可以分為海上固定式和海上浮式風機。相比于陸上風機，海上風機擺脫了土地利用以及風場大小的限制[2]。而由于近岸區(qū)域海岸線利用的限制，海上風機將逐漸走向深水區(qū)域，此時海上固定式風機的建造費用將劇增，不利于經(jīng)濟效益。海上浮式風機的概念應運而生，它能適應深水環(huán)境，具有發(fā)電穩(wěn)定、風能利用率高等優(yōu)勢，具有廣闊的應用前景。

海上浮式風機由風機、浮體平臺、錨鏈等多部分組成，受到風場、波浪、潮流等共同作用，受力非常復雜，在運動過程中可能出現(xiàn)失穩(wěn)、傾覆等問題，所以需要對系統(tǒng)作整體的耦合分析。耦合分析可以充分考慮結(jié)構(gòu)各部分以及各力之間的相互影響，在計算結(jié)構(gòu)運動、動力響應等方面能得到更準確的結(jié)果。Withee（2004年）[3]使用FAST和ADAMS開發(fā)了全耦合動力響應程序，用來計算浮式風機在風場和波浪作用下的響應。該程序同時考慮了浮體受到的非線性波浪荷載和風機受到的空氣動力荷載。Nielsen（2006年）[4]使用空氣動力學程序HAWC2和計算海洋結(jié)構(gòu)物動力響應的程序SIMO/RIFLEX，開發(fā)了計算浮式風機在風場、波浪作用下動力響應的程序。唐耀（2013年）[5]對Spar型浮式風機平臺做了耦合動力響應分析。對于風力荷載，采用Knauer等提出的經(jīng)驗公式進行簡化計算。對于水動力計算，使用SESAM軟件中的DeepC模塊計算系泊狀態(tài)下浮式風機系統(tǒng)的運動響應及錨鏈的張力。

以往在進行時域內(nèi)浮式風機水動力計算時常采用頻域轉(zhuǎn)時域的方法，如商業(yè)軟件SESAM中的一些計算模塊就采用這樣的方法。雖然此種方法計算速度快，但是需要先求出高頻下結(jié)構(gòu)物的水動力系數(shù)，這對于復雜結(jié)構(gòu)物是比較困難的。另外，頻域理論中通常假設浮體的運動幅度較小，而浮式風機由于受到風的作用，水平方向的漂移較大，這時此種方法將帶來誤差。針對這個問題，本文采用直接時域方法進行浮體水動力計算。這種方法對有限振幅運動的時域計算將更加準確和有效，同時也可以與氣動力計算和系泊計算方便地進行耦合。這樣建立的耦合計算方法將更真實地刻畫時域內(nèi)海上浮式風機的運動響應。

1 耦合計算方法

本文將建立一種用于計算海上浮式風機運動響應的時域耦合方法，并形成一套完整而獨立的時域耦合計算程序。對于風機所受的氣動力荷載，采用經(jīng)典的葉素動量法計算。對于浮體平臺的水動力荷載，在一階勢流理論下采用邊界元法進行計算。對于錨鏈模型，采用Chebyshev多項式擬合的懸鏈線模型計算。對于系統(tǒng)運動模塊，采用Runge-Kutta法求解。

1.1 風機氣動力理論

本文采用葉素動量法（blade element momentum method）進行風機氣動力計算[6]。葉素動量法由動量理論和葉素理論構(gòu)成。

根據(jù)動量理論，風機平面dr圓環(huán)上的軸向力可以表示為：

式中：v1為來流風速；r為圓環(huán)至圓心距離；ρ為空氣密度；a為軸向誘導因子。環(huán)形單元的扭矩為：

式中：ω為風機旋轉(zhuǎn)的角速度；a′為切向誘導因子。

根據(jù)葉素理論，將葉片分為若干個葉素段，每一段內(nèi)采用一種翼型數(shù)據(jù)。對于給定的葉素段，通過推導，風機平面dr圓環(huán)上的軸向力和扭矩可以表示為：

式中：B為葉片數(shù)；φ為相對來流風速與風機平面的夾角；c為此翼型的弦長；Cn為垂向分力系數(shù)。

讓式（1）與式（3）相等、式（2）與式（4）相等，即有：

具體計算時，不斷迭代軸向和切向誘導因子直至收斂，即可計算出dr圓環(huán)上的軸向力和扭矩，積分后即可得到整個風機的受力。

以上即為葉素動量法，為了獲得更精確的結(jié)果，這里進行了兩項修正。第一項為Prandtl葉尖損失修正。第二項為Glauert修正。

浮式風機在波浪的作用下將產(chǎn)生縱蕩、縱搖等運動形式，此時風機平面將發(fā)生變化，導致相對來流風速vrel和入流角度φ的改變，從而引起風機軸向力和功率的變化。而浮式風機在運動過程中還具有速度，也會導致相對來流風速的變化。比如當浮式風機發(fā)生縱搖時，風機平面發(fā)生偏轉(zhuǎn)，首先導致入流角度發(fā)生變化，這時垂直于風機平面的來流風速分量也將發(fā)生變化，因此迭代計算結(jié)果也會改變，即軸向和切向誘導因子會改變，最終風機軸向力和功率發(fā)生變化。因此計算浮式風機氣動力時需要考慮系統(tǒng)的位移和速度。

1.2 波浪與浮體作用的時域分析方法

由于浮式風機一般布置在深水區(qū)，其截面尺度相對于波浪要素基本上是小尺度，因此采用勢流理論計算一階波浪力就可以較準確地計算出其受力[7]。

一階勢流理論下速度勢函數(shù)Φ和波面函數(shù)η可以分為已知的入射分量和未知的散射分量：

式中：散射勢ΦS包含所有由于波浪和物體擾動而產(chǎn)生的散射分量。

散射勢在域內(nèi)滿足Laplace方程：

并滿足如下邊界條件:

式中：n為物面的法線方向，規(guī)定指出物體為正。

為了構(gòu)成定解問題，還需增加一個遠場輻射邊界條件。在數(shù)值計算中，采用增加人工阻尼層來保證散射波向外傳播而不反射回來。

利用簡單格林函數(shù)法，選取Rankine源和它關于海底的像作為格林函數(shù)：

式中：x為場點；ξ為源點。利用格林第二定律，可以得到關于計算域邊界上散射勢的邊界積分方程：

式中：邊界S包括淹沒在水中的物體平均濕表面SB和從物體到阻尼層外邊界的有限靜水面SF；α為固角系數(shù)。

根據(jù)散射勢的自由水面條件式（11）、式（12），采用高階邊界元法求解邊界積分方程即可計算出速度勢。得到速度勢后，就可以計算出物面的壓強及受力。對于線性波浪，作用在物體上的一階波浪力可以表示為：

1.3 錨鏈模型

對于海洋中的浮體結(jié)構(gòu)，需要錨鏈對其進行定位和約束。進行錨鏈力的計算時，常采用懸鏈線模型。具體計算時，常采用分段外推-校正法計算，即首先對錨鏈分段，對每一段求解，然后外推，經(jīng)過多次迭代校正后得到最終的結(jié)果。

但在時域內(nèi)的浮體—錨鏈系統(tǒng)中，每個時間步上浮體的位置都在變化，因此迭代的計算量會顯著增加。針對這個問題，本文采用Chebyshev多項式對錨鏈頂端位置和頂端拉力函數(shù)進行擬合，求出用于計算錨鏈頂端位移與受力關系的擬合公式[8]。

應用Chebyshev多項式進行擬合能夠根據(jù)錨鏈的頂端位置，方便、快速地求出錨鏈頂端的拉力，能夠應用在浮體—錨鏈系統(tǒng)的耦合計算中。

1.4 時域耦合模[10]

海上浮式風機的運動受到波浪力、錨鏈力、風機軸向力的相互耦合作用。對于一個給定的系統(tǒng)狀態(tài)（已知浮式風機的位移和速度），可以計算出此時各部分的外力。系統(tǒng)受到這些外力的作用，下一時刻的狀態(tài)就會改變。而波浪與結(jié)構(gòu)物、錨鏈與浮體、風與浮式風機各自都存在相互作用，而各部分間又有相互作用，因此使得整個問題成為一個非常復雜的耦合問題。

本文在耦合計算時采用的思路可以由圖1體現(xiàn)。系統(tǒng)運動模塊將系統(tǒng)的位移、速度傳遞給對應的計算模塊，各模塊計算后將力、力矩返回給運動方程，周而復始完成時域計算。圖1中的風機參數(shù)包括風機模型尺寸、葉片氣動力特性等參數(shù)，波浪參數(shù)包括波高、周期或者波浪譜的參數(shù)，錨鏈參數(shù)包括錨鏈強度、長度、類型等參數(shù)。

圖1 耦合計算流程圖Fig.1 Flow chart of the coupled calculation

對于氣動力模塊，系統(tǒng)運動模塊將為其提供浮體的位移和運動速度，根據(jù)給定的風機參數(shù)計算出當前時刻的風機軸向力和力矩；對于水動力模塊，系統(tǒng)運動模塊也將為其提供浮體的位移和運動速度，根據(jù)給定的波浪參數(shù)計算出當前時刻的水動力和力矩；對于系泊模塊，系統(tǒng)運動模塊只需為其提供浮體的位移，因為在錨鏈力的計算中，不涉及浮體的運動速度，根據(jù)給定的錨鏈參數(shù)計算出當前時刻的錨鏈力和力矩。這些力和力矩將返回到運動模塊中，根據(jù)運動方程計算出下一時刻的位移和速度。系統(tǒng)運動模塊分別與其余3個模塊有變量交換，完成了相互之間的耦合作用，因此這樣的模型是在時域上全部耦合的，能夠充分考慮各部分外力與整體系統(tǒng)的相互作用。

系統(tǒng)運動模塊中的運動方程可以表述為：

式中：Fk為總的廣義水動力荷載分量（包括力和力矩）；Gk為錨鏈作用力和力矩；Tk為風機軸向力和力矩。對于以上3個參數(shù)，k=1～3時分別表示x，y，z方向的力；k=4～6時分別表示對x，y，z軸的力矩。Mkj為物體的質(zhì)量矩陣；Bk為系統(tǒng)阻尼矩陣；Ckj為恢復力矩陣。

具體數(shù)值求解時，把式（18）的二階微分方程寫成如下的通用格式：

采用四階Runge-Kutta法求解此二階微分方程，計算時首先根據(jù)t時刻物體的位移ξ（t）和速度ξ˙（t），由水動力模塊確定波浪激振力、恢復力和阻尼力，由系泊模塊確定錨鏈力，由氣動力模塊確定風機軸向力，從而求得t時刻函數(shù)，再由Runge-Kutta法計算出t+Δt時刻新的位移ξ（t）和速度。重復該過程直至計算結(jié)束。

2 模型建立與驗證

計算模型采用的是2010年Jonkman提出的OC3-Hywind spar型海上浮式風機[11]。上部風機為NREL 5-MW風機。浮體部分為spar型圓柱結(jié)構(gòu)。采用3根對稱分布的錨鏈進行系泊。部分重要參數(shù)見表1—表2，錨鏈布置形式見圖2。

表1 風機與塔架參數(shù)Tab.1 Parameters of the wind turbine and tower

2.1 氣動力驗證

利用1.1中的氣動力理論編寫程序，對于不同風速下的固定式風機進行計算（不考慮風機實際的控制策略），并與其他程序和結(jié)果進行了對比驗證。

2.1.1 功率驗證

利用Blade Element Method計算不同風速下、風機固定時的功率，并與FAST程序結(jié)果進行比較，對比結(jié)果如圖3所示。

表2 浮體平臺與錨鏈參數(shù)Tab.2 Parameters of the platform and anchor chain

圖2 錨鏈布置形式Fig.2 Layout of the anchor chains

圖3 不同風速下固定式風機功率Fig.3 Powers of the fixed wind turbine in different wind velocities

對比顯示，在風速7～14 m/s時，程序計算結(jié)果與FAST結(jié)果基本吻合，當風速較大時，存在一定誤差。

2.1.2 推力驗證

11.4 m/s是該風機的額定風速，需要對此風速下風機受到的推力重點驗證。這里計算了該風速及其附近風速下風機受到的推力，并與文獻結(jié)果進行了比較，結(jié)果如表3所示。

對比顯示，額定風速及其附近風速下風機受到的推力（即軸向力）與文獻結(jié)果吻合。

綜合以上2個驗證，表明風機參數(shù)設置正確，氣動力模塊計算結(jié)果準確。

表3 風機推力結(jié)果Tab.3 Results of the axial force of the wind turbine

2.2 水動力驗證

利用1.2和1.3中的理論和模型，對浮式風機系統(tǒng)進行水動力計算?？紤]入射波為規(guī)則波，波高2.56 m，周期10 s，入射方向朝向x正方向。錨鏈采用懸鏈線模型。此時不考慮風的作用。由于浮體平臺和錨鏈布置的對稱特性，當來浪方向平行于x軸時，其運動響應只有在縱蕩、垂蕩、縱搖方向上有值，其余3個方向上為0。圖4—圖6是浮體縱蕩、垂蕩、縱搖的計算結(jié)果。

圖4 浮體縱蕩歷程及其幅值譜Fig.4 The history and the amplitude spectrum of surge of the floating platform

從運動響應可以看出，經(jīng)過一段時間后，浮體在平衡位置附近做規(guī)則的往復運動。從譜分析的結(jié)果也可以看出，譜峰頻率是入射規(guī)則波的頻率，說明浮體主要受波浪的影響。為了驗證水動力計算的準確性，將穩(wěn)定后的運動響應幅值與時域下錨鏈采用等效剛度以及頻域下錨鏈采用等效剛度的結(jié)果進行了比較，如表4所示。

圖5 浮體垂蕩歷程及其幅值譜Fig.5 The history and the amplitude spectrum of heave of the floating platform

圖6 浮體縱搖歷程及其幅值譜Fig.6 The history and the amplitude spectrum of pitch of the floating platform

從上述比較結(jié)果中可以看出，采用懸鏈線的時域計算結(jié)果與其他2種方法非常吻合，說明浮體模型參數(shù)設置準確，水動力模塊計算準確。

表4 浮式風機運動響應幅值比較Tab.4 Comparison among amplitudes of the motion response of the floating wind turbine

2.3 耦合計算驗證

浮式風機受到風、波浪、錨鏈的耦合作用，除了需要對每個模塊驗證準確性，還需對整體耦合計算模型進行準確性驗證。為了與已有的計算結(jié)果進行比較，本文進行了浮式風機平衡位置的計算，以此作為耦合模型的驗證。

風速為11.4 m/s，朝向x正方向。入射波波幅為0（存在輻射波，仍考慮水動力模塊），錨鏈采用懸鏈線模型。其縱蕩計算結(jié)果如圖7所示。

圖7 浮式風機縱蕩歷程Fig.7 The history of surge of the floating wind turbine

從圖7中可以看出，經(jīng)過一段時間后，縱蕩達到了穩(wěn)定，平衡在27.8 m。文獻[5]中給出的結(jié)果為25.7 m。考慮到所用的計算模型存在誤差，計算結(jié)果還是比較吻合的。平衡位置的耦合計算雖然沒有計入入射波的作用，但是依然考慮了水動力模塊，浮體的運動會產(chǎn)生輻射波浪力，也與氣動力形成了相互作用。因此，仍然是氣動力模塊、水動力模塊、系泊模塊的相互耦合計算。上述驗證說明整體系統(tǒng)模型參數(shù)設置準確，各部分模塊及耦合程序計算準確。

3 結(jié)果分析

3.1 規(guī)則波計算

入射波為規(guī)則波，波高2.56 m，周期10 s，入射方向朝向x正方向。錨鏈采用懸鏈線模型。風速11.4 m/s，朝向x正方向。圖8給出了2種計算中縱蕩結(jié)果的比較。其中wave-induced指的是只考慮水動力作用，不計入氣動力影響。wind-wave指的是氣動力和水動力的聯(lián)合作用。

圖8 規(guī)則波作用下浮式風機縱蕩歷程Fig.8 The history of surge with the regular wave of the floating wind turbine

從圖8（a）、圖8（b）中可以看出，在此計算條件下，氣動力的作用主要是使系統(tǒng)在縱蕩方向漂移了很長一段距離，穩(wěn)定后，兩者運動幅值與周期幾乎一致，說明此時系統(tǒng)主要受到波浪的作用。表5中給出了3個方向上運動的具體數(shù)值。

表5 浮式風機運動的平衡位置與幅值Tab.5 Balance position and amplitude of motion of the floating wind turbine

3.2 運動響應對風機功率的影響

計算條件同3.1中的耦合計算一致，這里計算了浮式風機的功率，如圖9所示。

由于初始效應，風機功率在0～300 s時有所波動，300 s后在一個平衡位置附近作周期性振蕩。平衡后的平均功率為5 216 kW。而同樣風速下，固定式風機的功率為5 363 kW。由此表明，由于浮體的運動，浮式風機的平均功率有所減小，減小約3%。

圖9 風機功率歷程Fig.9 The history of power of the wind turbine

4 結(jié)論

本文建立了一種用于計算時域內(nèi)海上浮式風機運動響應的耦合方法，主要由氣動力模塊、水動力模塊、系泊模塊、系統(tǒng)運動模塊構(gòu)成。前3個模塊通過與系統(tǒng)運動模塊相互傳遞變量進行耦合計算。通過氣動力驗證、水動力驗證以及耦合驗證，表明各模塊計算準確、各模塊之間的耦合也是準確的。這樣建立的時域內(nèi)海上浮式風機運動響應的耦合計算方法是準確而有效的。此耦合計算方法計算速度快，準確性高，可以用于海上浮式風機的設計及研究。

本文計算了OC3-Hywind spar型海上浮式風機在11.4 m/s的定常風，波高2.56 m、周期10 s的規(guī)則波作用下的運動及其功率變化。初始階段，浮式風機主要收到風機推力的影響，在縱蕩方向會漂移很大一段距離，穩(wěn)定值為27.8 m。穩(wěn)定后，主要受波浪作用，在平衡位置處做往復運動。對于功率而言，由于浮體的運動，浮式風機的平均功率減小了約3%。

[1]張志春，劉洪偉.我國發(fā)展風電產(chǎn)業(yè)的戰(zhàn)略性分析[J].電網(wǎng)與清潔能源，2013，29（5）:67-72.ZHANG Zhichun，LIU Hongwei.Strategic analysis of wind power industry development in China[J].Power System and Clean Energy，2013，29（5）:67-72（in Chinese）.

[2]鄭小霞，葉聰杰，符楊.海上風電場運行維護的研究與發(fā)展[J].電網(wǎng)與清潔能源，2012，28（11）:90-94.ZHENG Xiaoxia，YE Congjie，F(xiàn)U Yang.Research and development of operation and maintenance for offshore wind farms[J].Power System and Clean Energy，2012，28（11）:90-94（in Chinese）.

[3]WITHEE J E.Fully coupled dynamic analysis of a floating wind turbine system[D].Monterey，California: Naval Postgraduate School，2004.

[4]NIELSEN F G，HANSON T D，SKAARE B.Integrated dynamic analysis of floating offshore wind turbines[C]//American Society of Mechanical Engineers.25th International Conference on Offshore Mechanics and Arctic Engineering，June 4-9，2006，Hamburg，Germany:ASME，2006:671-679.

[5]唐耀.Spar型浮式風機平臺動力響應分析[D].上海:上海交通大學，2013.

[6]HANSEN M O L.Aerodynamics of wind turbines[J].Rotors Loads&Structure James，2015，5（2/3）:141-167.

[7]李玉成，滕斌.波浪對海上建筑物的作用[M].北京:海洋出版社，2002:95-98.

[8]滕斌，郝春玲，韓凌.Chebyshev多項式在錨鏈分析中的應用[J].中國工程科學，2005（1）:21-26.TENG Bin，HAO Chunling，HAN Ling.Numerical simulation of static behavior of the single anchor cable[J].Engineering Science，2005（1）:21-26（in Chinese）.

[9]SHIM S.Coupled dynamic analysis of floating offshore wind farms[D].Texas:Texas A&M University，2007.

[10]WAYMAN E N.Coupled dynamics and economic analysis of floating wind turbine systems[D].Massachusetts:Massachusetts Institute of Technology，2006.

[11]BUTTERFIELD S，MUSIAL W，SCOTT G.Definition of a 5-MW reference wind turbine for offshore system development[M].Golden，CO:National Renewable Energy Laboratory，2009:5-16.

[12]JONKMAN J M.Definition of the floating system for phase IV of OC3[M].Golden，CO:National Renewable Energy Laboratory，2010:2-21.

[13]馬鈺.單柱式浮式風機動力性能機理研究[D].上海:上海交通大學，2014.

（編輯 馮露）

A Coupled Time-Domain Calculation Method for Motion Response of Floating Offshore Wind Turbine

JIN Fei，TENG Bin
（State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian 116024，Liaoning，China）

In this paper，a coupled time-domain calculation method is established for the motion response of floating offshore wind turbine，mainly including aerodynamic，hydrodynamic，mooring and system motion modules.The blade element momentum method is used in the aerodynamic module and in the hydrodynamic module；the boundary element method is used for calculation under the first-order potential theory.In the mooring module，the catenary model is used with fitting calculation by the Chebyshevpolynomial.As for the motion module of the whole system，Runge-Kutta method is used to solve it.The OC3-Hywind spar wind turbine is modeled and the accuracy of this coupled calculation method is verified through the comparisons of each module and the coupled model.The method is used to calculate and analyze the motion response and its influence on the power of the wind turbine.

floating wind turbine；coupled time-domain method；blade element momentum method；boundary element method；catenary model

國家自然科學基金項目（51379032，51490672）。

Project Supported by National Natural Science Foundation of China（51379032，51490672）.

1674-3814（2016）08-0093-07

P751

A

2016-01-18。

金 飛（1991—），男，碩士研究生，研究方向主要為海上浮式風機的耦合計算。