規則和隨機橫浪中破損船舶運動響應研究*

吳 浩 高志亮 蘇 焱

(武漢理工大學高性能船舶技術教育部重點實驗室1) 武漢 430063) (武漢理工大學交通學院2) 武漢 430063)

0 引 言

隨著船舶尺度的增大，船舶發生碰撞或觸礁的概率增加.2018年初，油船“桑吉”號與散貨船“長峰水晶”號發生碰撞，“桑吉”號油船燃爆失火并最終沉沒；同年，挪威護衛艦“海爾格·英斯塔”號與油輪“索拉”號相撞，事故發生后護衛艦雖采用了海岸擱淺、損管作業，以及鋼索固定等補救措施，最終仍然沉沒.對大量的船舶碰撞事故事后分析表明，船舶破艙進水后正確的應急決策能將人員傷亡、財產損失，以及環境污染降至最低.破損船舶在海浪中的運動行為十分復雜，對其準確判斷決定了應急決策系統的正確性.

船模試驗方法難以對破損船舶運動的影響因素進行系統分析，不利于深入探究破損船舶與流體相互作用的動力學機理.數值模擬方法為上述問題研究提供了可行的途徑.早期，科研工作者將傳統耐波性理論與基于準靜態假設的破艙進水模型相結合，開展船舶在波浪中的運動研究[1-3].該方法采用勢流理論計算船舶波浪力，基于伯努利方程計算破口處的流量，并假設艙室內水面始終保持水平.上述方法計算簡單，但忽略了水體流經破口時的動態行為及艙室內部自由液面運動等瞬時因素，因此無法準確預報水流經破口和在艙內運動時的動態作用力.采用計算流體動力學(computational fluid dynamics，CFD)方法則可合理地計入破艙進水問題中流體的動態效應.

近年來，學者們采用基于求解納維-斯托克斯方程的CFD方法開展了破損船舶在波浪中的運動研究.Gao等[4-5]采用勢流理論求解船波干擾問題，并與CFD方法求解的破艙進水問題進行耦合，分析了破損滾裝客船在橫浪中的運動和受力.Sadat等[6]采用CFD方法模擬了破損客船在規則橫浪中的橫搖運動.Haro等[7]采用類似方法對同一艘船在破損狀態下規則迎浪航行時的快速性問題進行了研究.上述研究分析了破損船舶在規則波中的運動，而破損事故通常在真實海況中發生，因此，開展破損船舶在隨機波浪中的運動預報更加符合實際情況，同時也可為事故應急決策系統提供有用參考[8-9].

本文將基于求解雷諾平均納維-斯托克斯(reynolds averaged navier stokes，RANS)方程的CFD方法與譜分析方法相結合，研究了破損艦船DTMB 5415在規則橫浪和隨機橫浪中的運動響應.通過CFD方法建立數值波浪水池，進行了數值水池造波精度的驗證.采用所建立的數值水池模擬了破損船舶在規則橫浪中的運動，并將船舶橫搖頻率響應函數與試驗結果進行對比驗證.在此基礎上，采用頻譜分析方法預報了該艘破損船在四到六級海況的隨機橫浪中橫搖運動響應.

1 數學模型

本文采用CFD商用軟件STAR-CCM+對船舶運動流場進行模擬，將空氣和水看作不可壓縮流體，流體運動由連續性方程和RANS方程控制.在笛卡爾坐標系下，方程描述為

(1)

(2)

式中：t為時間；xi為坐標分量(i=1,2,3)；ui為xi方向的速度分量；ρ為流體密度；μ為黏性；I為單位矩陣；k為湍動能；p為壓力；F為體積力；Si為自定義源項.

采用SSTk-ω湍流模型，方程描述如下.

(3)

(4)

式中：ω為耗散率；Γk和Γω為k和ω的有效擴散項；Pk和Pω為k和ω的產生項；fβ*為自由剪切修正系數；fβ為漩渦拉伸修正系數；ω1和k1為環境湍流值.上述項的計算和系數取值見STAR-CCM+用戶手冊.

本文采用流體體積(volume of fluid，VOF)方法對水和空氣交界面進行捕捉，該方法通過計算網格單元中流體所占的體積分數(α)來捕捉流體的交界面，α滿足以下方程.

(5)

船舶在流體中的運動滿足

(6)

(7)

式中：vi和ωi分別為船體運動線速度和角速度；m為船體質量；Ji為船體轉動慣量；fi和Mi分別為作用于船體上的合力和合力矩.

STAR-CCM+軟件基于有限體積法對方程(1)～(5)進行離散，其中時間項采用二階隱式格式，方程(5)對流項采用二階高分辨率交界面捕捉(high-resolution interface capturing，HRIC)格式計算，其余方程中的對流項采用二階迎風格式進行計算，擴散項則采用中心差分格式進行計算.采用半隱式速度壓力耦合方程(semi-implicit method for pressure linked equation，SIMPLE)方法對流場的速度和壓力進行耦合求解.

2 數值預報及分析

2.1 研究對象

Begovic等[10]對DTMB 5415船模破損情況下在波浪中的運動開展了試驗研究.本文的研究對象為該船模(見圖1)，其主要船型參數和破損艙室參數見表1～2.破損艙室為位于船中附近的兩個艙室，兩艙之間艙壁位于破口中間.數值模擬中，在破損艙室頂部布置了與試驗相同的兩個通氣管，以便船舶在運動過程中艙室內的空氣可以順利排出.

圖1 船型示意圖

表1 DTMB 5415船型參數

表2 DTMB 5415破損艙室信息

2.2 數值波浪水池建立及驗證

STAR-CCM+軟件提供了基于Stokes波理論的邊界造波方法，一階Stokes波的流場水平速度(u1)、垂向速度(u3)和波面抬高(η)方程為

(8)

(9)

η=acos(Kx1-ω0t)

(10)

式中：a為波幅；ω0為波浪頻率；K為波數；x1方向為波浪傳播方向；x3為距靜水面的垂直距離；d為水深.

采用上述邊界造波方法生成一階Stocks波，波長λ取3.2 m，波高取波長的1/50.圖2為數值波浪水池的示意圖，其中計算域長8λ、寬0.5 m，水深2.15 m，該水深與文獻[10-11]試驗水池深度一致，水面以上部分高2 m.水池左側設為速度入口邊界，右側設為壓力出口邊界，前后兩側設為對稱面邊界，上下兩側設為不可滑移固壁邊界.

圖2 數值水池示意圖

為了減少波浪在水池末端反射影響，在方程(2)中施加阻尼源項，源項作用范圍為距離計算域末端2λ區域，源項表達式為

(11)

(12)

式中：xsd和xed分別為阻尼區的起始點和結束點.

數值水池的波面附近網格劃分見圖3，在波面區域內沿波高和波長方向進行網格加密；波高范圍內網格均勻分布，并向水池頂部和底部逐漸變疏；波長方向水池前端6λ區域采用均勻網格布置形式，水池后端2λ區域網格逐漸變疏.此處模擬二維波浪傳播，故在計算域寬度方向只設置了兩層網格.

圖3 數值水池網格劃分示意圖

首先對波浪傳播模擬進行了網格依賴性分析，通過對波高和波長方向的網格調整生成了三種疏密不同的網格，分別記為Mesh1、Mesh2和Mesh3，計算時間步長取波浪周期的1/1 000，在距入流口4λ處記錄波面高度，表3為計算結果與理論結果的比較.可以看出，隨著網格加密，計算所得的波高和波周期誤差降低.當網格達到Mesh2的分辨率時，誤差已小于1%，故后面計算根據Mesh2的網格布置形式來生成網格.

表3 網格依賴性分析結果(距入口邊界4λ處測點)

采用Mesh2的網格布置形式，通過選取三種時間步長，分別記為Timestep1、Timestep2和Timestep3，對計算模型進行了時間步長依賴性測試，時間步長依賴性分析結果見表4.隨著時間步長變小，計算所得的波高和波周期誤差降低；當時間步長小于波周期的1/800時，計算誤差小于1%，所以后面計算時間步長取波周期的1/800.

表4 時間步長依賴性分析結果(距入口邊界4λ處測點)

圖4為時間步長為波周期的1/800時，距入口4λ處測得的波高歷時曲線與理論值比較.圖5為波傳播20個周期后數值水池波面與理論值比較.可以看出，在工作區內，數值水池造波結果與理論結果基本一致；在消波區內，波浪朝著出口邊界方向逐漸衰減為零.

圖4 距入口4λ處計算波高與理論波高比較(Mesh2，Timestep2)

圖5 波傳播20周期時計算波高與理論波高比較(Mesh2，Timestep2)

2.3 破損船舶在規則橫浪中運動模擬

破損船舶在規則橫浪中運動模擬的計算區域及邊界設置見圖6，船體左右各取4λ，前后各取1倍船長，水深為2.15 m，水面以上高2 m.將整個計算域分為背景區域和重疊區域，在計算中背景區域保持靜止，而重疊區域隨船體一起運動.水面附近網格生成參考上述數值波浪水池Mesh2網格布置形式，艙室內網格尺寸為0.01 m.船體附近網格采用了局部加密的布置形式，遠離該區域網格逐漸變稀疏，見圖7.

圖6 破損船舶在規則橫浪中運動計算區域

圖7 船體附近網格布置

在計算域末端2λ區域采用3.2中描述的阻尼消波處理.此外，為了消除波浪在計算域入口邊界的二次反射，在方程(2)中源項加入力源項，力源作用區域為計算域前端2λ區域，力源項表達式為

(13)

采用上述模型模擬了破損船舶在不同波長的規則橫浪中的運動，波高取波長的1/50，破口面向來波方向.計算中僅考慮船體的橫搖、縱搖和垂蕩三個自由度的運動，初始時刻艙內水面與靜水面持平.表5為計算所得的船舶運動穩定后的橫搖幅值與試驗結果[12]比較，圖8為橫搖頻率響應函數計算結果與試驗結果比較.CFD計算結果與試驗結果吻合較好，當波長較小時由于船舶運動幅度較小，數值計算結果與試驗結果的相對誤差較大，但橫搖幅值絕對誤差小于1°，說明本文計算所得結果比較可靠.

表5 計算橫搖角與試驗值比較

圖8 橫搖頻率響應函數比較

圖9為波浪周期為1.47 s時船舶的運動歷時曲線圖，該波浪周期與船舶破損情況下橫搖固有周期相近.由圖9可知，計算14 s后船舶的橫搖運動趨于穩定，其橫搖角幅值為16.67°；船舶的縱搖運動幅度較小，其幅值約為0.2°，平均縱傾角約為1.1°；船舶的垂蕩運動幅值約為0.03 m.圖10為船舶運動不同時刻艙內水體的分布情況，由圖10可知，艙室內水體運動較為劇烈，出現自由液面爬升、翻卷和觸碰艙頂等現象.

圖9 船舶運動歷時曲線

圖10 艙內水體運動情況

2.4 破損船舶在隨機橫浪中運動預報

將上述CFD計算所得的船舶在規則橫浪中的頻率響應函數換算成實船尺度，在已知波浪譜密度的條件下，采用譜分析方法預報該實尺度破損船舶在隨機橫浪中的橫搖運動.波浪譜采用國際拖曳水池會議(international towing tank conference，ITTC)建議的雙參數譜，譜密度函數公式為

(14)

式中：H為三一波高；T為波浪特征周期.

表6為譜分析法得到的船舶橫搖運動的統計值.由表6可知，對于四級和五級海況，該破損船舶的三一橫搖幅值分別為1.47°和4.69°，對應的單幅有義值分別為1.4°和3.01°，均滿足北歐合作研究計劃建議的軍船耐波性衡準指標橫搖角單幅有義值小于4°的要求[13]，說明該船出現兩艙破損情況時在四級和五級海況下橫搖運動量較小；在這兩種海況下，船舶的百一橫搖幅值分別為2.46°和7.83°，表明當前破損條件下該船基本上不會出現較大的橫搖角度和傾覆行為.對于六級海況，該破損船舶的三一橫搖幅值為13.55°，對應的單幅有義值為5.06°，該值略大于北歐合作研究計劃建議的衡準值；船舶的百一橫搖幅值為22.63°，表明在六級海況下該船有可能出現較大的橫搖角度，需要防范傾覆行為的發生.此外，船舶破損后在上述三種海況下的橫搖運動的平均周期與船舶在完整狀態下的橫搖運動固有周期9.77 s較為接近，且均大于各自海況對應的波浪特征周期.

表6 破損船舶在隨機橫浪中的運動響應統計值

3 結 束 語

本文采用基于求解RANS方程的CFD方法對破損艦船DTMB 5415在規則橫浪中的運動進行了時域模擬.研究中建立了數值波浪水池并對其進行了網格和時間步長依賴性分析；在此基礎上模擬了破損船舶在規則橫浪中的運動，計算所得的橫搖運動響應與試驗結果吻合較好，表明本文采用的計算模型是有效的.對該船破損情況下在規則橫浪中發生橫搖諧搖時的運動進行了分析，得出其橫搖諧搖角幅值約為17°，平均縱傾角約為1°，垂蕩運動的幅值約為0.03 m.隨后，基于ITTC雙參數譜，采用譜分析方法預報了該破損船在隨機橫浪中的橫搖運動統計值.統計結果表明，船舶在四、五和六級海況下，單幅有義值分別約為1°，3°和5°，橫搖運動平均幅度較小；在六級海況下，船舶百一橫搖幅值約為23°，表明該船出現兩艙破損情況時在該海況下會出現較大的橫搖運動，需要注意其穩性相關的安全性問題.今后工作將在本研究基礎上開展船舶在波浪中發生瞬間破損后的運動行為研究.