半潛式海上風力機流固耦合特性分析

胡丹梅，曾理，陳云浩

（上海電力大學能源與機械工程學院，上海市 楊浦區 200090）

0 引言

隨著化石能源的日益消耗，人類對新能源的研究日益加深。風能作為一種清潔無污染的可再生能源，受到人們的廣泛關注[1]。隨著技術的進步，海上風電以其土地資源零占用、風速更高且更穩定、流場中湍流強度小、視覺及噪聲污染對環境影響小等優勢，近年來得到了許多國家的重視[2]。與陸上風電相比，海上風力發電面臨的發電環境更加復雜[3-5]，不僅會受到鹽霧對機組的腐蝕，還會受到波浪、浮冰、海流等對風力發電機組的沖擊[6]。

海上風力發電機按基礎劃分一般分為固定式基礎和漂浮式基礎。浮式基礎一般分為張力腿式(tension leg platform，TLP)、Spar 式、半潛式[7]，其適用范圍大多為水深超過50 m的海域。張力腿式對系泊纜材料要求較高。Spar 式基礎水平面外運動響應較大，較TLP 式更適應深水。半潛式基礎吃水小，在運輸和安裝時具有良好的穩定性，因此更具有經濟性[8]。

國內外很多專家學者對海上半潛式風力機進行了研究。Masciola[9]等人研究了系泊系統對半潛式海上風力機模型運動響應的影響，主要分析了3 種不同波浪環境條件對整個風力機模型運動的影響，研究發現系泊線動力學在影響浪涌和升沉半潛式運動中的作用有限。Cheng[10]等人建立了海上浮式風力機的水動力模型，分析了風力機的氣動性能，研究了氣動力對浮式平臺水動力響應的影響。劉海鋒等[11]基于流固耦合理論，采用雙向流固耦合方法和滑移網格技術對5 MW 風力機的尾跡特性進行了系統分析。鄧新麗等[12]利用Blade 軟件分析并計算了1.5 MW 風力機葉片在三維湍流下的載荷，發現葉片沿翼展方向應力最大且極限工況下的最大應力出現在離葉根1/3 的位置。Chen[13]等人基于2 種不同概念的半潛式海上風力機動力特性對比，從六自由度運動、氣動阻尼效應、陀螺效應、塔頂動力響應以及系泊系統的動力學入手分析其性能。

以上研究大多對半潛式風力機模型運動進行動力響應分析，不能全面地觀測風力機整機運行時的真實流場情況。本文除了對風力機進行整體的雙向流固耦合分析，還利用仿真軟件模擬實際運行時的環境載荷，在仿真軟件中展現風及海水流過風力機的實際情況，以此對比整機的形變，流線圖、形變云圖等的變化，研究結果可為海上半潛式風力機設計提供參考。

1 半潛式風力機模型及計算域模型

1.1 半潛式風力機模型

根據葉素理論將美國國家可再生能源實驗室(NREL)提供的參數進行坐標變換，進而將變換后的葉片坐標結果導入ProE構建葉輪模型。浮式平臺的選取來源于Marco Masciola[5]所提供的平臺數據，進而根據給定的東海水文數據對該半潛式浮式風力機進行建模。其中有關風力機整體構造數據及運行工況數據的詳細參數如表1所示。

根據提供的模型數據對半潛式漂浮式風力機進行整體建模，整機模型如圖1所示。

圖1 半潛式浮式風力機結構模型Fig.1 Semi-submersible floating wind turbine structure model

1.2 計算域模型

設置輪轂中心為坐標原點，風力機逆風方向為Y軸正向，軸心至塔筒方向為Z軸正向，風力機風輪繞Y軸旋轉，X軸方向由Y軸和Z軸右手法則決定。空氣場上部半圓半徑為150 m，氣場高度為242 m；具體計算域尺寸設置為：空氣場中的旋轉域以風輪輪轂中心為原點，半徑為70 m，厚度為5 m；下部分海水域高度設計為70 m；流體域入口距風力機150 m，風力機距流體域出口300 m。計算域整體模型如圖2所示。

圖2 計算域整體模型Fig.2 Overall model of computational domain

2 流固耦合理論

流固耦合作用是自然界客觀存在的一種特殊現象，是指流體與固體之間的相互作用。流固耦合問題可由其耦合方程定義，這組方程的定義域同時有流體域與固體域。而未知變量含有描述流體現象的變量和含有描述固體現象的變量，一般具有以下特征[14]：

1）流體域與固體域均不可單獨求解；

2）無法顯式地削去描述流體運動的獨立變量及描述固體現象的獨立變量。

在流固耦合過程中，為保證流體與固體間能量的守恒，則在流體與固體的交界面處二者應滿足以下方程：

式中：τ為應力；n為方向余弦；d為位移量；q為熱流量；T為溫度；下標f表示流體，s表示固體。

同時，可以建立控制方程的通式，對耦合簡化分析需要結合模型所處的初始條件參數及相關邊界條件。此時將流、固方程耦合到一個方程矩陣中進行求解，即在一個求解器中同時對流體和固體方程進行求解。方程如下：

式中：k為迭代時間步；Aff為流場的系統矩陣；Ass為固體的系數矩陣；ΔXfk為流體待求解力；ΔX ks為固體待求解力；Bf為流體外部作用力；Bs為固體外部作用力；Afs和Asf是流固耦合矩陣。

ANSYS 中采用的CFX 模塊和瞬態模塊中流體域計算與固體域計算之間流固耦合流程如圖3所示。

圖3 流固耦合流程圖Fig.3 Fluid-structure coupling flow chart

3 網格劃分及網格無關性驗證

3.1 計算域網格劃分及相關設置

由于需要確保雙向耦合作用下流場域與固體域的數據傳遞準確性，劃分時需要將葉片網格與旋轉網格接觸面的網格節點保持一致[15]。同時，由于風力機葉輪部分進行旋轉，需要對此區域網格進行加密，使旋轉時網格可以平穩過渡，確保數據傳遞的準確性。固體域網格如圖4所示，并在Transient structural 中設置整個風力機表面為流固耦合面，3 條系泊纜繩底面設為固定面。

圖4 風力機固體域網格Fig.4 Wind turbine structural domain grid

在風力機外部設置流體域，分別為風輪外部的旋轉域，包含塔筒、風輪及一半浮筒的空氣域，最后是下方的海水域。旋轉域以風輪輪轂中心為原點，半徑為70 m，厚度為5 m。空氣域上部半圓半徑為150 m，高度為300 m，下部海水域采用設計水深70 m。從流體域入口至風力機處為200 m，風力機距流體域出口為400 m，總長度為10R(R=63 m，為風輪半徑)。整個流場域建好后劃分網格，外部流場域采用六面體網格，內部旋轉域采用四面體，如圖5所示。

圖5 流體域網格Fig.5 Fluid domain grid

3.2 數值計算方法與邊界條件

本文基于Navier-Stocks 方程和RNGk-ε湍流模型[16]，利用滑移網格及網格重構對NREL 5 MW半潛式風力機進行數值模擬，采用二階迎風格式[17]。

控制方程為

式中：ρ為流體密度；φ為通用變量；t為時間；u為來流風速；Γ為廣義擴散系數；S為廣義源項。

RNGk-ε湍流模型為：

式中：k和ε分別為湍動能及其耗散率；ui為流體速度分量；xi，xj為流體坐標分量；μeff為風流黏度與湍流黏度之和；Gk為平均速度梯度產生的湍動能；Gb浮力所產生的湍動能；YM為可壓縮湍流中波動膨脹對總耗散率的貢獻；αk和αε分別為湍動能及耗散率普朗特數的倒數；C1ε，C2ε和C3ε為模 型默 認常 數；Sk，Sε和Rε分 別為 用 戶 定 義 的源項[18]。

CFX 定義的初始邊界條件：入口為速度進口，風速為11.4 m/s，水流速度為0.5 m/s；出口為自由出流；外圍流場區域設置為靜止域wall，內部旋轉域設置為interface；空氣域與海水域交界處設為interface面，整個風力機設為wall。設置旋轉域旋轉角速度為1.267 rad/s，風輪跟隨旋轉域一起轉動，設置網格重構。

3.3 功率驗證

為了驗證結果的準確性，采用不同網格數在機組的額定工況下進行數值計算。設置外流場入口為速度入口，出口為自由出流。額定風速為11.4 m/s，對應風輪轉速為12.1 r/min，風輪表面以及外流場邊界為無滑移壁面。

模擬采用SSTk-ω湍流模型，內外流場接觸面設為interface，采用二階迎風格式，收斂殘差設置為1×10-4。分別設置風速為5、8、11.4 m/s，求對應的計算輸出功率。

通過數值分析可以得到風輪轉矩M，由式(6)可計算得到風力機的功率。

式中：P為風力機輸出功率，MW；n為風力機額定工況下轉速12.1 r/min。

圖6 表明隨著輸入風速的增大，計算所得輸出功率越來越接近于設計功率5 MW。當風速設置為額定風速11.4 m/s 時，計算得到輸出功率為4.75 MW，相對于設計功率5 MW，誤差僅為5%，在可接受的誤差范圍內。

圖6 輸出功率驗證Fig.6 Computational domain model

4 風浪聯合作用下風力機耦合分析

4.1 雙向流固耦合計算域相關設置

利用Workbench平臺構建雙向流固耦合分析，在Mechnical 中構建結構域，在CFX 中構建流場域，在DM 中建立整個模型，并設置好相關邊界條件。空氣域與海水域交界面為interface 面，分別設置風速和水速，將整個半潛式風力機設置為耦合面。上部為空氣域，下部為海水域。設置固體域與流體域運行時間均為10 s，即總時間為20 s，時間步長設置為0.1 s。同時，設置動網格開啟，將各區域相連，實現數據傳遞。在計算中選擇二階歐拉方程，湍流選擇二階迎風格式，收斂殘差設置為1.0×10-4。圖7為計算域模型圖。

圖7 計算域模型圖Fig.7 Computational domain model

在流體域中劃分好流場網格，在結構域中劃分好整機網格，導入流體域并進行雙向流固耦合計算。在結構域中，風輪沿用NREL 5 MW 風力機中葉片設計的原材料，將玻璃纖維環氧樹脂設定為葉片的復合材料，該材料具有成型后尺寸穩定、抗腐蝕以及承受載荷較大等優點，浮式平臺的系泊纜選用聚酯纖維纜新型材料[19-20]。在流體域中分別設置空氣域、海水域的入口、出口及壁面，空氣域中的旋轉域需設為interface 面，再將計算結果以壓力載荷的方式導入結構，對風力機進行穩定性分析。

4.2 結構域應力應變分析

由文獻[21]所知，風力機的轉速與當時風速相關，本文所用NREL 5 MW 風力機轉速與風速的關系如表2所示。

表2 不同風速下風力機轉速Tab.2 Wind turbine speed at different wind speeds

根據表2，分別取5、8、11.4 m/s 風速下對應的風輪角速度以及風輪停止轉動時的角速度0 rad/s，設置系泊纜底部為固定約束面，同時設置半潛式浮式風力機表面為流固耦合面，輸出結構域中葉片的最大變形量。

圖8為不同風力機轉速下，葉片的最大變形量隨時間的變化規律圖。從圖8中看出，當風輪轉速為0 rad/s 時，風力機風輪停止轉動，處于自存狀態，葉片變形量隨時間變化幾乎不變，圖上顯示為一條直線，因此在自存狀態下，風輪變形量可忽略不計。在1.267 rad/s 轉速下，最大變形量峰值為21.038 mm；在0.959 rad/s 轉速下，最大變形量峰值為11.866 mm；在0.773 rad/s轉速下，最大變形量峰值則為7.757 mm。可見，隨著風輪轉速的增大，風輪的最大變形量也隨之增大，且其峰值隨轉速的增大出現左移，即風輪轉速越大，越早出現最大變形量。最大變形量在運行開始階段變化較大，隨著時間逐漸變小，在一定數值內趨于平緩。葉片變形主要受風力機轉速(即輸入風速)、材料影響，風輪轉速增大(即輸入風速增大)，變形量明顯提高。這一結論與張建平[20]所得結果一致，反映了葉片響應規律。在不同轉速下，最大變形量隨時間變化曲線的整體走勢基本一致。

圖8 不同轉速下葉片最大變形量對比Fig.8 Comparison of maximum blade deformation at different speeds

從圖8 可看出，在風輪停止運轉時最大變形量隨時間變化曲線幾乎為直線，因此在等效應力曲線圖中只討論3種風輪轉速，如圖9所示。由圖9可知，在風輪旋轉初期，最大應力的變化很大，隨著時間的推移，風輪最大等效應力逐步趨向平穩，并穩定在一定數值內變化，此時可判斷風輪旋轉,趨于穩定。且隨著輸入風速的增大，風輪轉速隨之增大，最大等效應力也隨之增大。當風力機風輪轉速設置為0.727 rad/s時，風輪的最大等效應力為8.01 MPa；風輪轉速為0.959 rad/s時，其最大等效應力為12.44 MPa；風輪轉速為1.267 rad/s時，其最大等效應力為21.46 MPa。由此可見，風輪等效應力值明顯小于設計的材料許用應力值250 MPa。

圖9 風輪最大等效應力隨時間變化圖Fig.9 Maximum equivalent stress of the wind wheel changes with time

在流固耦合設置中將風速都設置為額定風速11.4 m/s。在額定風速下，將海水速度分別設置為0.5 m/s與2 m/s，對半潛式風力機浮式平臺進行雙向流固耦合特性分析。將雙向流固耦合結果導出，在CFD-POST 中以云圖的方式展示結果，如圖10所示。

圖10 風力機浮筒位移圖Fig.10 Displacement cloud diagram of wind turbine foundation

浪速為0.5 m/s 時，浮筒最大位移為0.09 m。浪速為2 m/s時，浮筒最大位移為0.37 m。隨著浪速的提升，浮筒位移顯著增大，且位移最大處基本出現在下浮筒與系泊纜連接處附近。在浮筒部分，位移主要體現在系泊纜連接處，且系泊纜由于與浮筒相連，隨浮筒一起運動。

為了更好地觀察浮筒的運動情況，在CFX中導出浮筒及系泊纜位移參數，畫出浮筒及系泊纜位移隨時間變化曲線，如圖11所示。

圖11 浮筒及系泊纜位移隨時間變化曲線圖Fig.11 Float foundation and mooring line displacement curve diagram with time

由圖11可以看出，在5.5 s左右位移達到最大值。可以合理地推測當時間進行到11 s 時，浮筒回到初始位置，完成一個周期運動。從原因分析，由于平臺與系泊纜相連，平臺受風浪聯合作用，同時受到系泊纜拉力的作用，在變形過大時及時將平臺拉回初始位置。

4.3 流體域速度流線分析

在額定風速11.4 m/s，海浪速度2 m/s工況下，在CFD-POST 中導出半潛式風力機流場域速度矢量圖，如圖12、13所示。

圖12 風力機風輪附近速度矢量圖Fig.12 Velocity vector diagram around wind turbine wheel

圖13 流體域速度矢量圖Fig.13 Velocity vector diagram of fluid domain

由圖12 可知，風輪處速度由葉尖向葉根遞減，在半潛式風力機中速度最大存在于葉尖部分，且葉尖處速度矢量密度最大，速度約為87.5 m/s。在風輪部分可以清晰地觀察到風輪旋轉時的速度矢量。在半潛式平臺部分，發現圍繞系泊纜區域，速度矢量明顯高于海水域其他部分，易形成繞流。圖14為在額定風速11.4 m/s、海浪速度2 m/s 工況下風力機的流線圖，其中用粗實線標出的為風力機葉輪所在位置。從圖14中可以看出，風流經風力機后發生偏轉，流線有一定紊亂，但經過旋轉域后逐漸恢復穩定。

圖14 風力機流線圖Fig.14 Wind turbine streamline diagram

5 結論

對半潛式風力機進行雙向流固耦合特性分析，結論如下：

1）在不同風速對應的風輪轉速下，隨轉速的增大，風輪形變增大，且峰值出現左移，即轉速越大，越早出現最大變形量。最大變形量的大小隨著時間的變化，逐漸趨于平緩，并最終穩定。葉片變形主要受風力機轉速(即輸入風速)、材料影響，其中受風速影響較大。在運行開始后0.3 s，風輪變形量達到最大值。

2）隨著時間推移，風力機半潛式平臺位移與等效應力值都趨于穩定。風力機等效應力值在所設工況中符合設計要求，明顯小于設計的材料應力值250 MPa。且在結構分析中未出現較大程度變形，符合設計結構穩定性要求。

3）在風力機運行過程中，葉片的變形集中于葉尖附近，且越靠近葉尖越明顯。在整機中，形變主要集中于系泊纜連接處。因此，葉片及浮筒系泊纜的材料選取尤為重要。經測算，浮筒與系泊纜的振動周期約為11 s。因此在系泊纜處應進行加固或選用更結實的材料。

4）根據半潛式風輪速度矢量圖，可看出葉尖至葉片中段速度較大，顏色較深。風力機流線圖顯示風經過風輪產生一定擾動，后逐漸穩定，流場流線基本穩定。