一種半潛式商用6 MW級風力機載荷數值模擬研究

方龍, 翟恩地, 李榮富, 寧巧珍, 趙斌, 李曄, 周雅容, 章麗駿

(1.天津大學 建筑工程學院, 天津 300350; 2.北京金風科創風電設備有限公司, 北京 100176; 3.南方海洋科學與工程廣東省實驗室(湛江), 廣東 湛江 524002; 4.上海交通大學 船舶海洋與建筑工程學院, 上海 200240)

海上風電具有風資源儲量高、遠離居民區、選址方案限制低等諸多優勢,近年來發展迅猛。我國近海海域的海上可開發風資源豐富,海上年均風速高,尤其是南海的部分海域,其年均風速值在5.5～10.7 m/s區間[1],適合安裝大型海上風力機裝置。根據支撐基礎的形式,海上風力機可分為固定式和漂浮式。隨著海上風電產業逐漸走向深遠海,固定式風力機在其結構可靠性及經濟性方面均難以滿足需求。漂浮式風力機由于其在深遠海海域的極強的適用性,已經成為全球海上風電領域新的研究熱點[2]。我國也在浮式風電領域積極探索,國內首臺漂浮式風電機組“三峽引領號”[3]于2021年12月在廣東陽江海上風電場成功并網發電。2022年5月,國內首臺深遠海浮式風力機“扶搖號”[4]在廣東湛江羅斗沙海域進行示范應用。這標志著我國在全球率先具備大容量抗臺風型深遠海漂浮式海上風電機組自主研發、制造、安裝及運營能力。

浮式風力機工作在復雜的風浪流條件下,風力機、浮式平臺以及系泊系統之間的作用相互耦合[5]。作用在葉輪上的力對系統的載荷有貢獻,浮式平臺的六自由度運動又改變了葉輪的氣動性能,系泊系統進一步復雜化了整個漂浮式風力機系統的載荷與運動響應,漂浮式風力機的載荷分析是典型的多場-多體耦合問題[6]。同時,深遠海環境的極端天氣給漂浮式風力機的安全運行帶來了巨大的挑戰,因此,系統地評估漂浮式風力機在各種風浪流條件下的運動響應和載荷變化至關重要。由于漂浮式風力機成本昂貴,目前針對漂浮式風力機的研究以模型實驗和數值分析居多。

在實驗研究方面,國內外一些學者對不同形式的漂浮式風力機在不同環境下的動力響應進行了大量的實驗研究。Koo等[7-8]對比了Spar和半潛式等不同形式的漂浮式風力機在不同波浪條件下的水動力響應,研究了不同浮式基礎對塔筒固有頻率的影響;Shin[9]研究了一種新型張力腿浮式平臺在不同波浪條件下的水動力響應;Duan[10]評估了風載對浮式風力機縱蕩和縱搖的影響,以及風載對系泊系統的動態響應的影響,揭示耦合響應的本質;Chen[11]針對模型浮式風力機的推力遠小于實驗弗勞德縮放值的問題,發展了一種葉片的高階再設計方法,可以顯著提高轉子推力,提高氣動力相似水平。

在數值分析方面,美國國家可再生能源實驗室還與DNV GL、丹麥可持續能源國家實驗室和挪威理工大學合作開發了數個數值模擬模塊,如FAST[12]、HAOC2[13]等,并進行了代碼比較協作[14-15]。張立等[16-17]結合葉素-動量理論、輻射/繞射理論及有限元方法,研究了不同風浪條件下對各種浮體形式的漂浮式風力機平臺運動響應的影響;丁勤衛[18]基于輻射/繞射理論對比研究了垂蕩板及其安裝位置對漂浮式風力機Spar平臺動態響應特性的影響。這些方法通常使用莫里森方程結合勢流理論方法來計算水動力載荷,該方法一般忽略粘性項,引入阻力系數彌補粘性損失。此外,勢流理論中的阻尼模型不能考慮與渦旋脫落有關的橫向力或升力,這將影響了浮體在橫向方向上的預測運動響應的精度[19]。蔡恒等[20]結合勢流理論和葉素動量理論對比研究了單純波浪以及風浪聯合作用對浮式風力機系泊力的影響。姜勁等[21]對一種5 MW半潛式垂直軸風力機進行了氣動、水動及系泊載荷的全耦合數值模擬,探究了速比對平臺縱蕩運動響應以及葉輪氣動載荷的影響規律。

本文以一種商用6 MW級漂浮式風力機為研究對象,在荷蘭MARIN實驗室對其縮尺模型進行了風浪耦合試驗研究。同時,基于Spalart-Allmaras湍流模型、流體-剛體動態相互作用模型和重疊網格技術,對該模型進行了氣動力-水動力-系泊力的全耦合數值模擬。研究了波浪狀況對浮式風力機平臺運動狀態以及葉輪氣動性能的影響,并與實驗值進行了對比,分析了引起浮式風力機氣動載荷波動的原因。

1 研究對象

本文所述仿真模型為一種商用6 MW級浮式風力機的水池試驗縮比模型,該機組根據我國南海的實際海況設計,縮比后的模型如圖1所示。在該漂浮式風力機的水池實驗[22]中,試驗模型整體滿足弗勞德數相似準則,該相似準則可以盡可能地滿足模型機組與原機組之間的幾何相似、動力相似和剛度相似。在滿足弗勞德數相似的條件下,葉輪空氣動力的雷諾數相似難以實現,但可以通過調整葉片的外形參數,達到模型機組與原機組在推力系數隨葉尖速比變化趨勢上的一致。限于實驗條件,確定縮尺比例為1∶55,縮尺后的模型機組的基本尺寸如表1所示。

表1 實驗模型浮式風力機基本物理量參數Table 1 Basic physical parameters of the scaled test model

圖1 水池實驗縮比模型Fig.1 Scaled model of the test in offshore basin

2 數值模型

本文基于非穩態的Spalart-Allmaras 雷諾平均湍流模型,利用重疊網格技術和滑移網格技術,以及DFBI六自由度運動模塊,借助Star-CCM+軟件模擬漂浮式風力機在風浪作用下的運動響應與載荷變化。SA(Spalart-Allmaras)湍流模型是工程上廣泛應用的湍流模型之一,具有計算收斂性好和可靠性高等優點。對于漂浮式風力機,其運動響應受到水動力載荷和空氣動力載荷的影響,同時,其運動狀態也同樣會影響作用在機組上的載荷。因此,漂浮式風力機載荷和運動響應是一個動態的相互作用過程。流體-剛體動態相互作用模型(dynamic fluid body interaction, DFBI)可以通過模擬剛體在流體域中的壓力和剪切力作用下的運動,以及系泊線的恢復力,該模型可以很好地模擬浮式風力機與流體的相互作用。

2.1 計算域設置

仿真計算域的設置參考實驗條件,圖2給出了計算域示意圖,入口邊界和出口邊界到浮式風力機的距離分別為5.2 m和12 m,兩側邊界到浮式風力機的距離均為5 m,水深為2.46 m,頂部邊界到水面的距離為7 m。為了盡量減少網格對于風速及波浪的耗散,在入口邊界(圖2中左側邊界)、頂部邊界和兩側邊界均設為速度入口。出口邊界(圖2中右側邊界)設為壓力出口,以保證計算域內外壓力平衡。在本文中,不考慮海床的粗糙度,將海底視作一個光滑的壁面,因此計算域的底部邊界設為壁面。

圖2 計算域尺寸定義及邊界設置Fig.2 Size definition and boundary setting of the fluid domain

2.2 網格劃分

本文選用切割體網格對計算域進行離散,如圖3,利用重疊網格和滑移網格技術,將計算域劃分為背景網格域、重疊網格域以及滑移網格域。重疊網格域將浮式風力機包絡在其中,用來模擬機組的六自由度運動,滑移網格域為葉輪所在的區域,用該域網格的旋轉來模擬葉輪的旋轉。為了提高模擬精度,對葉輪的尾流區和水面線附近區域進行加密。葉片近壁面第1層網格高度為5×10-5m,經過驗證,葉片大部分位置的y+值接近1,滿足SA湍流模型的近壁面網格計算要求。

圖3 網格劃分示意Fig.3 Diagram of the mesh distribution

3 數值驗證

仿真采用隱式非穩態計算,需要設置時間步長來實現時間上的離散。一般情況下,時間步長越小,模擬精度越高,計算量越大。反之,時間步長越大,計算量越低,模擬精度也越差。對于漂浮式風力機,空氣動力和水動力的模擬對時間步長的要求不同,為了在控制計算成本時不影響模擬精度,本節分別從葉輪推力和浮體自由衰減2個方面來討論時間步長對仿真結果的影響。

3.1 時間步長無關性驗證

為了驗證空氣動力模擬所需的時間步長要求,本文在保持浮體固定的情況下模擬了葉輪的氣動載荷,分別采用了0.05、0.03、0.02、0.01、0.005和0.003 s 6個時間步長。圖4為葉輪推力系數隨時間步長變化的結果,可以看出,時間步長越小葉輪的推力越大,并逐漸趨于平穩,時間步長為0.005 s和0.003 s的推力系數值的相對誤差為0.17%,說明當時間步長滿足Δt≤0.005 s時,可以達到空氣動力的模擬精度要求。圖5給出了4個時間步長下浮式平臺在垂蕩方向上的自由衰減模擬值,時間步長越大,垂蕩自由衰減得越快,隨著時間步長的減小,平臺垂蕩的自由衰減值趨向于收斂,時間步長為0.01 s和0.005 s的衰減周期的相對誤差為2.05%,即當時間步長滿足Δt≤0.01 s時,可以達到水動力計算的精度要求。

圖4 葉輪推力系數時間步長驗證Fig.4 Time step validation of the thrust coefficient

圖5 自由衰減時間步長驗證Fig.5 Time step validation for free decay

綜合考慮水動力和空氣動力的模擬精度要求,即當時間步長滿足Δt≤0.005 s時,可滿足模擬漂浮式風力機的精度要求,故后文模擬均采用0.005 s的時間步長。

3.2 土地利用動態度。在人和自然的共同作用下，區域內各種土地利用類型的數量在不同時間段變化的速度是不同的。土地利用類型的變化率與土地利用動態度可定量描述區域土地利用變化的速度。單一土地利用類型動態度可運用土地利用動態度模型，分析土地利用類型的動態變化，進而真實反映區域土地利用類型的變化速度和劇烈程度。公式為[13]：

3.2 浮體自由衰減對比

自由衰減通常用于確定浮體的自然周期,設置波浪條件為靜水,給定浮體一個初始位移,然后釋放,使其從初始位置自由移動,為了驗證浮式風力機水動力模擬的準確性,對浮體在縱搖、橫揺、縱蕩和橫蕩上的自由衰減進行了模擬,并與實驗測量值進行對比。圖6和圖7分別給出了浮式平臺在縱搖和橫揺上的自由衰減模擬值與實驗值的對比。結合表2可以看出,模擬值在縱搖、橫揺和垂蕩上與實驗值吻合程度較高,縱蕩和橫蕩的模擬值在衰減周期上大于實驗值,艏搖衰減周期的模擬值低于實驗值,其原因為數值模擬中采用懸鏈線模型對系泊系統進行了簡化,該模型計算得出的系泊系統作用在浮體上的載荷與實際值存在一定誤差。因此,全耦合數值模擬的方法對浮式風力機水動力的模擬較為準確。

表2 浮式風力機自然周期Table 2 Natural period of floating wind turbine s

圖6 縱搖自由衰減對比Fig.6 Comparison of free decay results of pitching motion

圖7 橫揺自由衰減對比Fig.7 Comparison of free decay results of rolling motion

4 浮式風力機的運動響應與載荷分析

本文分別對靜水、規則波和非規則波作用下的漂浮式風力機進行了模擬,研究漂浮式風力機在不同外部條件作用下的載荷與運動響應,并在時域和頻域上分析了漂浮式風力機的載荷波動與運動響應的關系。

4.1 靜水工況

本文模擬了靜水中的浮式風力機在來流風作用下的載荷與運動響應,來流風速設為U0=1.53 m/s,葉輪轉速設為ω=9.32 rad/s,對應風力機葉尖速比為RTS=9.47。模擬過程中監測了浮式風力機的六自由度運動,通過轉換可得到浮式風力機機艙的位移、速度和加速度。本文僅對來流風作用下機艙沿X軸方向(假設X軸正方向與風的來流方向相同)的運動進行分析。風力機機艙的沿X軸方向的位移xhub、速度vhub和加速度ahub的表達式為:

(1)

(2)

(3)

式中:xsurge代表漂浮式風力機平臺沿X方向的線位移;h為機艙中心距離風力機重心的高度;θpitch代表漂浮式風力機繞Y軸的角位移。機艙相對速度指機艙運動速度與來流風速的比值,其表達式為:

(4)

圖8給出了浮式風力機葉輪推力系數與機艙加速度波動的時歷曲線,可以看出機艙加速度和葉輪推力呈正相關,當葉輪推力增大時,浮式風力機的傾覆力矩越大,機艙沿來流方向的加速度也越大;反之,當葉輪推力減小時,浮體回復力矩和系泊力的作用大于葉輪的推力作用,機艙加速度減小為負值。

圖8 定常風靜水工況下機艙加速度與葉輪推力系數Fig.8 Nacelle acceleration and turbine thrust coefficient curves under the condition of steady wind and still water

為了分析浮式風力機葉輪載荷波動的原因,圖9給出了葉輪推力系數的頻域分布,可以看出,推力系數的峰值①的頻率為4.250 Hz,接近浮式風力機葉片的通過頻率4.449 Hz,說明在純風作用下,浮式風力機氣動載荷的波動主要受葉輪轉動的影響,一般情況下,葉片通過頻率附近氣動載荷的影響由風剪切和塔影效應造成,由于模擬中的采用的是均勻來流而非剪切來流,因此排除風剪切的影響,說明葉輪推力系數的波動是由塔影效應造成的。

圖9 定常風靜水工況下葉輪推力系數頻域Fig.9 Frequency domain analysis of rotor thrust coefficient under the condition of steady wind and still water

4.2 規則波工況

對規則波作用下漂浮式風力機進行模擬,模擬與實驗的來流、波浪以及風力機的運行參數一致,具體參數如下:來流風速為1.53 m/s,波高0.091 m,波周期1.079 s,浪向和來流風向一致,葉輪轉速為9.32 rad/s,對應風力機葉尖速比為9.47。

由于縱搖和縱蕩直接影響風力機的實際來流速度,推測這2種運動對風力機氣動載荷的影響最大,通過合成縱搖和縱蕩在機艙高度處沿來流方向的運動研究機艙運動與葉輪推力之間的關系。圖10給出了機艙沿X軸方向相對速度與推力系數波動的時歷曲線,模擬推力系數的平均值為0.703,實驗測量的平均推力系數為0.671,相對誤差為4.77%,說明全耦合數值模擬的方法較好地模擬了漂浮式風力機在風浪作用下的氣動載荷。規則波作用下機艙運動的波動十分規律,呈現出正弦曲線的變化趨勢。漂浮式風力機縮尺模型的縱搖運動時歷曲線的數值模擬和水池模型試驗結果如圖11所示,可見兩者在波動頻率上幾乎一致,而水池模型試驗的縱搖幅值略大于模型試驗,這主要是由于模型實驗中的物面粗糙度和系泊結構的簡化處理以及數值模擬時沒有環境噪聲的干擾。因此可以認為,在風浪耦合工況下,本文所述研究方法對于半潛式風力機的數值模擬結果可靠。此外,推力系數和機艙相對速度在1.078 s的周期上均有較大的波動,且波動峰值出現的時間節點接近,當機艙位移速度最小時,速度值為負,即機艙向來流的反方向運動,葉輪與空氣的相對速度最大,此時葉輪的推力系數也達到了最大。相反,當機艙位移速度最大時,速度值為正,即機艙和來流同方向運動,葉輪與空氣的相對速度最小,此時的推力系數也達到最小值,即葉輪推力系數變化和機艙沿來流方向的速度呈負相關。說明浮式風力機縱搖和縱蕩運動對風力機的氣動載荷影響較大,其原因為縱蕩和縱搖運動使得機艙在來流方向上前后運動,進而影響了葉輪平面的實際來流速度。

圖10 規則波下浮式機組機艙速度與葉輪推力系數Fig.10 Nacelle velocity and rotor thrust coefficient under regular wave

圖11 漂浮式風力機縱搖運動曲線Fig.11 Pitch motion curves of the floating wind turbine

圖12給出了浮式風力機推力系數的頻域分布,可以看出,推力系數峰值①的頻率為0.924 Hz,該頻率接近波浪頻率0.927 Hz,說明規則波工況下,浮式風力機氣動載荷的波動主要受波浪的影響。進一步證明在波浪載荷的影響下,浮式風力機發生縱搖和縱蕩引起葉輪實際來流速度的波動,從而導致了葉輪氣動載荷的波動。推力系數的峰值①頻率為4.276 Hz,接近葉片的通過頻率4.449 Hz。根據4.1節的分析,該頻率氣動載荷的波動是由塔影效應造成的。圖13和圖14分別給出了規則波下浮式風力機縱蕩和縱搖的頻域分布,可以看出,浮式風力機縱蕩和縱搖的峰值頻率均出現在波浪頻率0.927 Hz附近,說明規則波工況下,浮式風力機的運動的主要受波浪載荷的影響,且縱蕩和縱搖運動的頻率與波浪頻率一致。圖13與圖14中在②號峰值頻率左側分別出現了一個低頻頻峰,該低頻運動為平臺的固有運動頻率導致的。在圖13中,②號峰值頻率右側分別出現了多個高頻頻峰,其分別為②號峰值頻率的2、3和4倍頻。由于葉輪的氣動載荷相較于浮體的水動力載荷小得多,浮式風力機的縱蕩和縱搖運動在葉片通過頻率處的峰值不明顯,說明波浪載荷對浮式風力機運動的影響比氣動載荷大得多。

圖12 規則波下葉輪推力系數頻域圖Fig.12 Frequency domain of rotor thrust coefficient under the regular wave

圖13 規則波下浮式風力機縱蕩運動頻域圖Fig.13 Frequency domain of the floating wind turbine′s surging motion under the regular wave

圖14 規則波下浮式風力機縱搖運動頻域圖Fig.14 Frequency domain of the floating wind turbine′s pitching motion under the regular wave

4.3 非規則波工況

對非規則波作用下的漂浮式風力機進行模擬,來流風速為3.371 m/s,有義波高0.12 m,譜峰周期1.48 s,波譜采用JONSWAP譜,浪向與來流同向,葉輪轉速9.32 rad/s,對應風力機的葉尖速比為4.29。

圖15給出了非規則波下機艙X軸方向相對速度與推力系數波動的時程曲線,相較于規則波工況,非規則波工況下浮式風力機機艙的運動速度以及推力系數的波動更加無序。風力機推力系數和機艙運動速度直接相關,當機艙速度為正,即機艙向來流同方向運動時,葉輪與空氣的相對來流速度減小,推力系數整體上低于平均值;反之,當機艙速度為負,即機艙向來流的反方向運動時,葉輪與空氣的相對來流速度增加,推力系數整體上高于平均值,即葉輪推力系數變化和機艙沿來流方向的速度呈負相關,浮式風力機縱搖和縱蕩運動對風力機的氣動載荷影響較大。

圖15 非規則波下浮式機組機艙速度與推力系數Fig.15 Nacelle velocity and thrust coefficient under the irregular wave

對非規則波工況下浮式風力機的推力系數進行頻域的分析,如圖16,推力系數頻域峰值①的頻率為0.623 Hz,接近非規則波譜峰周期對應的頻率0.674 Hz,說明非規則波工況下,浮式風力機氣動載荷的波動主要受波浪載荷的影響。推力系數頻域峰值②的頻率為4.468 Hz,接近風力機的葉片通過頻率4.449 Hz,同理,該頻率下風力機氣動載荷的波動是由塔影效應造成的。圖17和圖18給出了浮式風力機在非規則波下縱蕩和縱搖的頻域分布,縱蕩和縱搖峰值的頻率均集中在0.674 Hz附近,而在葉片通過頻率和葉輪轉頻上未出現峰值,同樣說明在非規則波工況下,浮式風力機縱搖和縱蕩主要受波浪的影響,葉輪的氣動推力的波動對其運動的影響可以忽略不計。

圖16 非規則波下葉輪推力系數頻域圖Fig.16 Frequency domain of rotor thrust coefficient under the irregular wave

圖17 非規則波下浮式風機縱蕩運動頻域圖Fig.17 Frequency domain of the floating wind turbine′s surging motion under the irregular wave

圖18 非規則波下浮式風機縱搖運動頻域圖Fig.18 Frequency domain of the floating wind turbine′s pitching motion under the irregular wave

5 結論

1)模擬值在縱搖、橫揺和垂蕩上與實驗值吻合程度較高,縱蕩和橫蕩的模擬值在衰減周期上略大于實驗值,艏搖衰減周期的模擬值低于實驗值,其原因為本文所述研究方法中采用懸鏈線模型對系泊系統進行了簡化,難以完全匹配浮式平臺復雜的六自由度運動狀態,計算得出的載荷與試驗測得的載荷數值存在一定誤差。

2)在靜水工況下,浮式風力機僅受定常風作用,葉輪氣動載荷的波動主要由塔影效應產生,機艙加速度的波動與葉輪推力系數變化呈正相關。

3)在規則波和非規則波工況下,葉輪推力的波動受浮式基礎縱蕩和縱搖的影響較大,其主要原因為這2種運動使得機艙沿來流方向前后運動,進而導致了葉輪的實際來流的速度不斷發生變化,葉輪推力系數變化和機艙速度呈負相關。另外氣動載荷的波動還與葉片的轉動相關,這主要是塔影效應造成的。

4)在規則波和非規則波工況下,浮式風力機的運動主要受波浪載荷的影響,氣動載荷對其影響相對較小。