吳 浩 高志亮 蘇 焱

(武漢理工大學(xué)高性能船舶技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室1) 武漢 430063) (武漢理工大學(xué)交通學(xué)院2) 武漢 430063)

0 引 言

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

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

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

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

1 數(shù)學(xué)模型

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

(1)

(2)

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

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

(3)

(4)

式中：ω為耗散率；Γk和Γω為k和ω的有效擴(kuò)散項(xiàng)；Pk和Pω為k和ω的產(chǎn)生項(xiàng)；fβ*為自由剪切修正系數(shù)；fβ為漩渦拉伸修正系數(shù)；ω1和k1為環(huán)境湍流值.上述項(xiàng)的計(jì)算和系數(shù)取值見STAR-CCM+用戶手冊.

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

(5)

船舶在流體中的運(yùn)動滿足

(6)

(7)

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

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

2 數(shù)值預(yù)報(bào)及分析

2.1 研究對象

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

圖1 船型示意圖

表1 DTMB 5415船型參數(shù)

表2 DTMB 5415破損艙室信息

2.2 數(shù)值波浪水池建立及驗(yàn)證

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

(8)

(9)

η=acos(Kx1-ω0t)

(10)

式中：a為波幅；ω0為波浪頻率；K為波數(shù)；x1方向?yàn)椴ɡ藗鞑シ较?；x3為距靜水面的垂直距離；d為水深.

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

圖2 數(shù)值水池示意圖

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

(11)

(12)

式中：xsd和xed分別為阻尼區(qū)的起始點(diǎn)和結(jié)束點(diǎn).

數(shù)值水池的波面附近網(wǎng)格劃分見圖3，在波面區(qū)域內(nèi)沿波高和波長方向進(jìn)行網(wǎng)格加密；波高范圍內(nèi)網(wǎng)格均勻分布，并向水池頂部和底部逐漸變疏；波長方向水池前端6λ區(qū)域采用均勻網(wǎng)格布置形式，水池后端2λ區(qū)域網(wǎng)格逐漸變疏.此處模擬二維波浪傳播，故在計(jì)算域?qū)挾确较蛑辉O(shè)置了兩層網(wǎng)格.

圖3 數(shù)值水池網(wǎng)格劃分示意圖

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

表3 網(wǎng)格依賴性分析結(jié)果(距入口邊界4λ處測點(diǎn))

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

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

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

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

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

2.3 破損船舶在規(guī)則橫浪中運(yùn)動模擬

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

圖6 破損船舶在規(guī)則橫浪中運(yùn)動計(jì)算區(qū)域

圖7 船體附近網(wǎng)格布置

在計(jì)算域末端2λ區(qū)域采用3.2中描述的阻尼消波處理.此外，為了消除波浪在計(jì)算域入口邊界的二次反射，在方程(2)中源項(xiàng)加入力源項(xiàng)，力源作用區(qū)域?yàn)橛?jì)算域前端2λ區(qū)域，力源項(xiàng)表達(dá)式為

(13)

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

表5 計(jì)算橫搖角與試驗(yàn)值比較

圖8 橫搖頻率響應(yīng)函數(shù)比較

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

圖9 船舶運(yùn)動歷時曲線

圖10 艙內(nèi)水體運(yùn)動情況

2.4 破損船舶在隨機(jī)橫浪中運(yùn)動預(yù)報(bào)

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

(14)

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

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

表6 破損船舶在隨機(jī)橫浪中的運(yùn)動響應(yīng)統(tǒng)計(jì)值

3 結(jié) 束 語

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