基于CFD的風電安裝船首部砰擊載荷分析

葉 晗 黃亞南 李飛虎 馮 珺 康 琪

（大連海洋大學 航海與船舶工程學院 大連 116023）

0 引 言

在全球可再生能源產(chǎn)業(yè)日益增長的背景下[1]，風能因其巨大的環(huán)境效益和商業(yè)潛力而獲得快速發(fā)展。為實現(xiàn)2050年零排放的目標[2]，有必要進一步發(fā)展風能，特別是海上風能。為此，各國加快了海上風電項目的建設。由于深遠海域風能資源更加豐富，因此走向深遠海是海上風電項目的未來趨勢。然而，海上風電場的基礎、建設和安裝成本占開發(fā)成本的比例較高，比陸地上的相應成本高出近3.6倍。鑒于海上風電場所處的風況、海況和海床的具體特點，水波和風壓將極大地影響風電場開發(fā)的安全性和經(jīng)濟性。艏部位置是風電安裝船重要受力處，因此研究艏部砰擊載荷尤為重要。

波浪砰擊是船舶海洋工程領域的一個重要問題，也是造成船體結構破壞的主要原因之一。自20世紀60年代初以來，國內外學者相繼開展大量砰擊研究，包括理論研究、試驗研究和計算仿真研究等。近些年來，有關砰擊問題的理論研究主要集中于船體局部垂直入水砰擊（即二維砰擊理論），該研究相對簡單且更為成熟。COINTE[3]研究了二維情況下水波沖擊剛性結構物的問題，采用一種基于波浪沖擊力的分析方法，通過計算分析得出了波浪砰擊的時間歷程、瞬時沖擊力的分布和砰擊面積等參數(shù)；同時，該研究還討論了波浪砰擊對結構物穩(wěn)定性和疲勞壽命的影響，并提出了相應的建議和措施。HOWISON等[4]研究了小入水角度情況下，不可壓結構物的入水問題，并進行了數(shù)學建模和理論分析。研究發(fā)現(xiàn)，入水問題涉及到很多物理現(xiàn)象（如空氣被夾帶進入水中，水面的形變和擾動等），論文提出了一種基于邊界元方法的數(shù)值計算模型，可以較為準確地描述這些物理現(xiàn)象的發(fā)生和演化過程，通過計算分析便可得出不同入水角下的入水形態(tài)、水下壓力分布等參數(shù)。

研究船舶波浪砰擊問題的最終目的是為了在實際工程中應用。陳占陽等[5]研究了外張砰擊載荷的時間分布規(guī)律，其通過模型試驗，采用高速壓力傳感器測量砰擊載荷的時間歷程，并進行分析和擬合，得到了外張砰擊載荷時間分布的數(shù)學模型。洪堯等[6]對圓盤入水砰擊載荷進行試驗研究，在不同入射速度和角度下，對入水的砰擊載荷進行測試和分析，研究結果表明入射速度和角度都對砰擊載荷有顯著影響，會產(chǎn)生劇烈的波浪和噴濺現(xiàn)象。

隨著科技的飛速發(fā)展，海上工程越發(fā)普遍，也引發(fā)更多實際問題。傳統(tǒng)的理論方法已經(jīng)難以解決各種實際問題，其局限性變得越來越嚴重。而實體試驗研究不僅需要耗費大量時間和成本，還存在很多不可抗因素，使得人們開始將目光轉向數(shù)值仿真計算，以解決砰擊問題。

司海龍等[7]結合計算流體力學（computational fluid dynamics, CFD）數(shù)值模擬計算技術和船舶的時域運動預測方法，得出可適用于一般船型的砰擊載荷預測方法。研究表明，砰擊載荷系數(shù)的大小，與船舶入水速度、浪向角以及船舶首部在波峰、波谷位置都有非常大的關系。劉正國等[8]為了計算某江海直達船整船在波浪中的運動，使用AQWA軟件分析得到船體表面三維形狀和船舶航行速度對船舶所受砰擊載荷的影響，并提出可行性優(yōu)化方案。吳小平等[9]利用Nastran軟件，基于頻域法對船舶尾部與首部偏大的客滾船進行計算模擬，將砰擊載荷反應到有限元計算模型的各部分上，得到客滾船砰擊強度的安全評估報告。張文挺等[10]對計算船舶的運動響應進行研究，基于時域非線性方法并結合三維Rankine源法，得到了船舶的垂向相對運動速度，進而結合砰擊力系數(shù)計算出船舶的砰擊載荷。

本文基于CFD并使用STAR-CCM+商業(yè)軟件與重疊網(wǎng)格技術，結合六自由度剛體運動模型，自由液面捕捉采用流體體積（volume of fluid,VOF）方法，在規(guī)則波（斯托克斯五階波）作用下，對不同波高和波長下的風電安裝船首部砰擊載荷進行計算。

1 CFD理論基礎

1.1 控制方程

基于雷諾平均N-S方程，考慮流體不可壓縮，連續(xù)性方程和雷諾平均N-S方程如下：

式中：iμ為時均速度，m/s；ρ為流體密度，kg/m3；t為時間，s；ν為流體的運動黏度，mPa·s；fi為質量力，N；P為流體壓強，Pa；μ'i為脈動速度，為雷諾應力項。自由液面捕捉采用VOF方法。

1.2 k-ε湍流模型

本文選用標準k-ε模型，其穩(wěn)定性較好并且計算精度高，湍流黏度方程表示為：

式中：Cμ代表經(jīng)驗常數(shù)，取值為0.09。

k和ε對應運算方程如下所示：

式中：C1ε、C2ε、σk、σε分別取值1.44、1.92、1.0、1.3；Gk為湍動能產(chǎn)生項，表達式如下：

1.3 造波方法

本文采用斯托克斯五階波，其中包含了多個非線性項和耗散項。斯托克斯五階波的方程如式（7）至式（10）所示。

波浪剖面方程η：

水平速度u：

垂直速度v：

速度勢φ：

式中：波長為L，m；波浪周期為T，s；波數(shù)為k；波速為c，m/s；波形系數(shù)為nφ；速度勢系數(shù)為nη。

2 船型參數(shù)與工況選取

2.1 幾何模型與邊界條件

目標船型吃水較淺、甲板面積大，可存放數(shù)量較多的風機設備。船上配備1臺主起重機、2臺輔吊機以及船首定位裝置，可在最大作業(yè)水深為40 m的近海進行作業(yè)。風電安裝船主尺度如表1所示，全船表面模型如下頁圖1所示，同時在船首位置布置了16個觀測點，參見下頁圖2。

圖1 風電安裝船模型

圖2 船首觀測點

表1 風電安裝船主尺度

由于設置船舶為迎浪，船橫搖角度較小，并且流場具有對稱性，所以采用半船計算域。計算域中船首方向為速度入口，船尾方向為壓力出口，船首距離速度入口為1倍船長，船尾為2倍船長，如圖3所示，為了防止出口處波浪反射對流場造成影響，在出口處激活阻尼消波。

圖3 計算域及邊界條件

2.2 工況選取

為了探究船首處波浪載荷與波浪參數(shù)的規(guī)律，根據(jù)風電安裝船工作特性設置船舶為迎浪運動，船體滿載吃水5.5 m。本文設計了3組不同波高，分別為2 m、3 m、4 m；5組不同波長，分別為66.5 m（0.5倍船長）、106.4 m（0.8倍船長）、133.0 m（1倍船長）、159.6 m（1.2倍船長）、266.0 m（2倍船長）。共計15組計算工況，如表2所示。

表2 各組工況參數(shù)

為了計算分析便利，本文采用模型縮尺比1∶20，三維數(shù)值水池中對應各波高分別為0.10 m、0.15 m和0.20 m，波長分別為3.325 m、5.320 m、6.650 m、7.980 m和13.300 m。

3 網(wǎng)格劃分與數(shù)值方法驗證

本文采用VOF方法進行自由液面捕捉，將計算區(qū)域分為若干個小單元，每個小單元內的流體狀態(tài)用1個分數(shù)來描述，表示該單元內流體的體積分數(shù)。通過對這些小單元體積分數(shù)的計算，可以得到流體的總體積和質量。本文為更好模擬波浪運動的變化，劃分數(shù)值水池網(wǎng)格時，在自由液面上下一個波峰處進行了網(wǎng)格細化，而其余網(wǎng)格則采用線性遞增的方式進行劃分，以節(jié)省計算資源。這種網(wǎng)格劃分方法可以更好地捕捉自由液面變化，提高模擬準確度。自由液面處的網(wǎng)格劃分方案如圖4所示。

圖4 水和空氣體積分數(shù)示意圖

為驗證數(shù)值方法的穩(wěn)定性，采用速度入口邊界造波，且對自由液面處劃分不同尺度的網(wǎng)格，討論網(wǎng)格密度對造波產(chǎn)生的影響?；诎惭b船模型，選定1個計算工況（波高0.2 m、波長6.65 m），在波浪場中，X軸方向上，由原點向正負方向等距離（4 m）選取觀測點，如圖5所示。

圖5 波浪場示意圖

隨著波浪衰弱進行，離速度入口越遠、衰弱越明顯，所以著重觀測速度入口遠端點（-8觀測點）。計算時間步長取0.01 s，1個波高方向網(wǎng)格取10 ～30個進行比較。

當時間步為0.01 s時，1個波高方向網(wǎng)格設置為10、20和30個，并進行加密；網(wǎng)格高分別為0.020 m、0.010 m、0.007 m，且自由液面處單個網(wǎng)格的橫向間距x與垂向間距z的比例一般選取為1∶2、1∶4或者1∶8。具體比例可根據(jù)波陡參數(shù)進行選取，本文選取x/z=1∶4，橫向網(wǎng)格密度與縱向保持一致。數(shù)值水池網(wǎng)格具體劃分形式如表3所示。

表3 數(shù)值水池網(wǎng)格劃分方案

將不同網(wǎng)格劃分下的觀測點距原點時歷曲線與理論值相比較（如圖6所示），可發(fā)現(xiàn)在計算時間內：當網(wǎng)格為10個時，波浪衰弱較明顯；當網(wǎng)格為20個和30個時，波浪衰弱都不明顯，符合計算要求。綜合計算時間成本考慮，當網(wǎng)格大于20個再增加時，計算結果基本無明顯變化，因此選用網(wǎng)格為20個、網(wǎng)格高為0.01 m，作為本文計算網(wǎng)格劃分方案。

圖6 不同網(wǎng)格高度下，觀測點距自由液面高度

4 數(shù)值計算結果及分析

本文通過邊界造波方法，實現(xiàn)了斯托克斯五階波的數(shù)值模擬，該方法以不可壓縮黏性流體的N-S方程和VOF自由液面捕捉法為基礎，借助Star-CCM+軟件定義了入射邊界處波浪解析表達式，從而模擬了五階斯托克斯波的產(chǎn)生與傳播。

風電安裝船在波浪中航行，定義其設計航速為12.0 kn?？s尺比1∶20下，設計航速為0.6 kn，同時波高分別為0.10 m、0.15 m和0.20 m，波長分別為3.325 m、5.320 m、6.650 m、7.980 m和13.300 m。如圖2所示，共取16個觀測點。當波高為0. 15 m、波長為6.65 m時，比較觀測點1、5、9、13（x不同，z相同）；點2、6、10、14（x不同，z相同）；點3、7、11、15（x不同，z相同）；點4、8、12、16（x不同，z相同）。這4組受力情況見圖7至圖10，圖11為船體壓強云圖。

圖7 觀測點1、5、9、13砰擊載荷對比曲線

圖8 觀測點2、6、10、14砰擊載荷對比曲線

圖9 觀測點3、7、11、15砰擊載荷對比曲線

圖10 觀測點4、8、12、16砰擊載荷對比曲線

圖11 船體壓強云圖

可見，此4組總體受壓情況為離水線越近，所受砰擊載荷變化越大。本文重點研究砰擊載荷的變化，所以重點分析距水線較近的1組，參見圖7。

由此可知，離船首越近時砰擊載荷越大，砰擊載荷曲線變化劇烈，短時間內即可到達極值。此種砰擊會導致船舶加速度出現(xiàn)突變，對船舶結構破壞較大；同時也發(fā)現(xiàn)，高度相同時，砰擊載荷的周期幾乎相同。

當波高為0.15 m、波長為6.65 m時，比較觀測點1、2、3、4（x相同，z不同）；點5、6、7、8（x相同，z不同）；點9、10、11、12（x相同，z不同）；點13、14、15、16（x相同，z不同）。這4組的受力情況，見圖12至圖15。

圖12 觀測點1、2、3、4砰擊載荷對比曲線

圖13 觀測點5、6、7、8砰擊載荷對比曲線

圖14 觀測點9、10、11、12砰擊載荷對比曲線

圖15 觀測點13、14、15、16砰擊載荷對比曲線

由此可知：離水線面越遠，砰擊載荷越小，距水線面較遠的2組觀測點，由于幾乎沒接觸水面，所以砰擊載荷較小，因此重點比較離水線較近的觀測點。

通過上述比較發(fā)現(xiàn)，船首最前端靠近水線面處的砰擊載荷為最大值。為了探究成因，通過比較船?？v搖和砰擊載荷曲線發(fā)現(xiàn)，船模垂蕩周期與砰擊載荷周期相同。如圖16所示，由此可以推斷，船首砰擊載荷產(chǎn)生的主要原因是船首的出水與入水所造成；波浪沖擊船首外飄區(qū)造成的外飄砰擊較小，可能是由于外飄面積大且作用時間長。由于點13處的砰擊載荷為最大值，所以下文將主要對點13處的砰擊載荷規(guī)律展開研究。

圖16 船?？v搖和砰擊載荷曲線

當波高為0.15 m時，比較不同波長下（波浪周期），點13處所受的砰擊載荷，如圖17所示。當波長等于船長時，所受砰擊載荷最大。由船舶耐波性可知，在規(guī)則波中，當L/λ在1附近時，波浪的擾動力最大，縱搖和垂蕩都十分劇烈。由此又證明了船首的波浪載荷主要是由于船舶在波浪中出水、入水所導致。

圖17 觀測點13在不同波長下的砰擊載荷時歷曲線

由以上結論可以判斷出，當波長等于船長時，船首靠近水線面處所受砰擊壓力最大。由此對波長等于6.65 m時，點13處在不同波高的砰擊載荷進行研究，參見下頁圖18。

圖18 觀測點13在不同波高下的砰擊載荷時歷曲線

由此可見，當波高為0.2 m時，點13處所受砰擊載荷最大。

5 結 語

綜上所述，風電安裝船作業(yè)中：隨著波高增大，船舶首部砰擊載荷越大；隨著波長增長，砰擊載荷呈先上升后下降的趨勢；當波長等于船長時，砰擊載荷達到最大值。船首砰擊主要由船舶在波浪中出水、入水所導致，其劇烈程度與船舶縱搖幅值相關，減少砰擊載荷就是減少縱搖幅值，因此船舶應當加強首部結構強度。并且，因距水線面越近的點所受砰擊壓力越大，所以風電安裝船首部距水線面較近位置的結構強度尤其需要重點加強。

本文分析了船首位置砰擊載荷分布以及載荷峰值最大值，為風電安裝船首部強度設計提供了參考。