基于CFD的漁船破艙進水時域模擬及破艙穩(wěn)性和壓強分布研究

于 博 文，張 維 英，陳 靜，于 欣，于 洋，周 俊 秋

(大連海洋大學 航海與船舶工程學院，遼寧 大連 116023 )

0 引 言

據中國漁業(yè)互保協會統計，1994～2015年的漁船事故中，全損(無法修復)漁船2 996艘，部分(可修復)損失漁船76 454艘．事故原因方面，由人員操作不當造成的碰撞事故占比44.81%，觸損事故占比39.12%，觸礁事故占比49.23%，擱淺事故占比36.04%[1]．所以探討漁船在破艙條件下的穩(wěn)性，尤其是非對稱進水過程的六自由度變化，可以判斷漁船的分艙合理性和給出破艙應急操作的建議，以達到避免或減少漁船事故損失的目的，對漁船的安全性具有實用價值[2]．

不斷發(fā)生的海難事件促使人們不斷發(fā)展船舶破艙穩(wěn)性的研究，在模型試驗之外對于破艙進水的研究方法可以分為兩種．準靜態(tài)法：將整體分成破艙管涌、破損艙室內的水運動以及波浪運動下的船舶運動[3-9]．在這些研究中，基于伯努利方程水動力模型被用來確定通過開口的流速．假定涌入艙內的水是瞬間沉降到一個水位[3，5-7]，并作為塊狀物在特定路徑表面上自由移動[8-9]，或者晃蕩并且遵循淺水方程[4]．勢流理論法：應用切片理論[3，4，6]或面元法[5，7-9]計算船-波相互作用動力學．計算機硬件的飛速發(fā)展，增強了數值模擬用于解決非線性問題的能力．破艙進水過程中水的涌入、沖擊、飛濺等都是典型的非線性狀況，通過數值模擬方法解決非線性問題的應用日趨成熟，對基于計算流體動力學方法的高保真度模擬的需求也在增加．

劉強等[10]基于Fluent軟件用Navier-Stokes(N-S)方程的CFD求解器對進水以及二維艙室耦合運動進行模擬，從艙內空氣對水進入艙內的影響分析，提出了減少外域計算量的方法，并探討了空氣壓縮性對于進水過程的影響．盧俊尹[11]用模型試驗法和一種雷諾平均的N-S(RANS)方程對船舶固定狀態(tài)下不考慮外流體域的破艙進水問題進行了研究，兩者結果吻合較好．Manderbacka等[12]為了增強水流和晃動的耦合效應，對集總質量模型進行了改進，使其能夠考慮入流沖量和流量對分隔室內湍流的影響，并且具有計算效率高和實用性強的優(yōu)點．Sadat-Hosseini等[13]使用CFD中的Ship-Iowa代碼(RANS方程求解器)，通過水平設置方案和動態(tài)重疊網格技術，模擬了破浪船在橫波海況中的淹沒過程．Ming等[14]利用光滑粒子流體動力學(SPH)研究相對于入射波的損壞位置對船艙行為的影響．Cao等[15]在2019年建立了一個基于條帶理論的振動模型，用于計算不同損傷和荷載條件下受損艙室的運動響應．Gao等[16]在2020年基于RANS方程結合流體體積法(VOF)、滑動界面和動力數值模擬，并采用分層技術處理網格更新，對護衛(wèi)艦DTMB-5415在靜水、規(guī)則波中的運動進行了模擬，理論計算和實驗結果高度吻合．

CFD逐漸與實驗流體力學一起成為船舶穩(wěn)性研究中的重要手段，其動網格和重疊網格技術在對于非穩(wěn)態(tài)運動物體的運動和強烈的非線性動力學的模擬中，有著很好的應用．

本文對一條養(yǎng)殖看護船DLHT 892 L288先用COMPASS進行穩(wěn)性計算，再選取一個艙室并用UG NX建模，在考慮內外流體影響下，用建立的STAR-CCM+數值模擬方法計算船體和進水之間橫搖的耦合運動，并對船模破艙進水進行時域模擬．通過對破損船和波浪的橫向耦合，盡可能還原作業(yè)環(huán)境下漁船的穩(wěn)性變化和壓力分布，對船舶破艙后扶正及減損提供參考．

1 數學模型

馬崢等[17]基于Fluent軟件分別用Standardk-ε、RNGk-ε、SSTk-ω、Realizablek-ω和Standardk-ω湍流模型在同一條件下對不同船舶進行模擬并分析．對于方形系數較小、航速較高的船，Standardk-ε和SSTk-ω模型有相對較好的預測精度．SSTk-ω模型對方形系數較大、航速較低的船的總阻力也具有較高的預測精度．而RNGk-ε模型對剩余阻力具有較強的預測能力．綜上所述，本文選取Standardk-ε模型．

1.1 湍流模型

1.1.1 不可壓縮黏性流體控制方程 本文以三維不可壓縮的黏性流體瞬態(tài)運動方程為理論基礎，流體密度和黏性系數為常數．

質量守恒方程(連續(xù)性方程)為

(1)

式中：u、v、w分別為速度矢量v在x、y、z軸上的分量．

動量守恒方程(運動方程)為

(2)

式中：F為質量力，p為壓強，μ為流體動力黏度．

有限體積法是在控制體積內對一般形式的控制微分方程的積分，即求解積分形式的守恒方程：

(3)

式中：φ為通用變量，V為控制體積，Γ為廣義擴散系數，S為廣義源項．

1.1.2 湍流方程 湍流的物理模型基礎是三維非定常的，這使得模擬湍流的計算量非常大．前人為了在有限計算資源下模擬湍流，提出了大渦模擬(large eddy simulation，LES)和RANS兩種方法．RANS法適用于本文．

連續(xù)方程

(4a)

動量方程

式中：u′i表示略去平均符號的雷諾平均速度分量，ρ為密度，u′j為脈動速度，σij為應力張量分量；i=1，2，3．

k-ε模型湍流控制方程如下：

湍流動能方程

Gb-ρε-YM+Sk

(5a)

擴散方程

式中：Gk為速度梯度產生的湍動能項；Gb為浮力產生的湍動能項；i，j=1，2，3分別表示x、y、z方向；YM為脈動擴張項，C1ε、C2ε、C3ε為經驗常數；σk、σε分別為與湍動能k和耗散率相對應的Prandtl數；Sk和Sε為用戶自定義的源項．

1.2 運動控制方程

船舶在海上運動時，把船體運動看成剛體運動，其平移運動通過高斯-賽德爾(GS)描述的線性方程求解：

(6)

船舶的旋轉運動受以下固定于質心的相對坐標系(BS)中的動量方程控制：

(7)

式中：Jc是船舶相對于中心的慣性矩張量，并且相對于BS保持恒定；Ω是船的角速度矢量；M′c是相對于中心作用在船上力矩的合成矢量．

1.3 入射波的產生與消散

通過背景區(qū)域的邊界條件規(guī)定入射波方向，入射波由斯托克斯(Stokes)波理論確定．

水平速度方程

u=Aωcos(K·x-ωt)ekz

(8)

垂直速度方程

w=Aωsin(K·x-ωt)ekz

(9)

表面高度方程

η=Acos(K·x-ωt)

(10)

式中：A為波幅值，ω為波頻率，K為波矢量，k為波矢量的幅值，z為與平均水位的垂直距離．

為了減少波在反射邊界附近的振蕩，可以在選定邊界附近添加阻尼衰減．阻尼將垂直阻力引入垂直運動，以避免其在邊界反射．Choi等的方法[18]為垂直速度w方程添加了一個阻力項：

(11)

其中

(12)

式中：xsd為波阻尼(在x方向上傳播)的起點，xed為波阻尼的終點(邊界)，l為阻尼波長度，f1、f2和nd為阻尼模型的參數．

2 模型建立與計算域選取

2.1 實驗船和破損艙室建模

在本研究中，養(yǎng)殖看護船DLHT 892 L288的主要參數如表1所示．

表1 養(yǎng)殖看護船DLHT 892 L288的主要參數Tab.1 Main parameters of breeding care boat DLHT 892 L288

依據型值表給出的船型參數，在建模軟件UG-NX中，以尾垂線與BL的交點為坐標原點，船頭為x軸正方向，左舷為y軸正方向，鉛垂向上為z軸正方向建立基準笛卡兒坐標系．依據型線，采用1∶1進行網格曲面的建立，各曲面之間采用G1連續(xù)相切，如圖1(a)所示；在同一個坐標系中，選取該船的工具艙(艙容為28.00 m3)并對該艙室進行建模，舷側板厚為0.01 m，如圖1(b)所示．

將模型導入STAR-CCM+中，并構建該模型破損艙室頂部的通風孔，其直徑為1.50 m，保證船艙進水時，排除內部氣壓對進水的影響；其破損進水口的直徑為0.30 m，如圖2所示．

圖2 受損船DLHT 892 L288的數值模擬模型Fig.2 Model of damaged boat DLHT 892 L288 in the numerical simulation

(a)船體

船艙的內表面與船外表面均為混合壁面，其切向速度固定，剪應力無滑移，壁面規(guī)格平滑．

為較好地捕捉進水流量特征，進水口的最小尺寸為0.005 m，最大尺寸為0.010 m，對艙內區(qū)域的空氣和水體積采用各向異性加密的方法，其網格尺寸的z方向為0.050 m，x和y方向均為0.200 m．為合理表現出船外表面高曲率部分的曲面特征，表面曲率最大點數為200，其面網格最小尺寸為0.010 m，最大尺寸為0.400 m，劃分結果如圖3所示．

(a)局部加密

本文選取DLHT 892 L288的裝載情況為空載到港，部分裝載情況為No.1壓載艙100%．No.2壓載艙60%，船員100%，淡水10%，燃油10%，食品及備品10%．該裝載下排水量為109.4 t，平均型吃水1.575 m．

通過COMPASS計算該船空載返港狀況下船舶復原力臂曲線，如圖4所示，橫坐標為橫傾角θ，縱坐標為靜穩(wěn)性力臂Ls．

圖4 DLHT 892 L288空載返港狀態(tài)下的復原力臂曲線Fig.4 DLHT 892 L288 restoring force arm curve under no-load return to port state

2.2 外計算域

2.2.1 外計算域劃分 在本文中，監(jiān)測破損船分別在靜水、常規(guī)海況和危險海況作業(yè)下的六自由度變化．為更好地響應復雜運動，采用嵌套網格技術和VOF法用于破損船六自由度運動的網格更新．

將整個區(qū)域分為6個計算域，如圖5所示．其中區(qū)域1為船舶運動的背景域；區(qū)域2為船舶運動的滑道域；區(qū)域3為嵌套網格域，其與船舶和艙室區(qū)域同步運動，產生相同的平移和旋轉；區(qū)域4為波浪域；區(qū)域5是船舶本體域；區(qū)域6為破損艙室域；破損船的外計算域尺寸見表2．

圖5 計算域劃分Fig.5 Partition of the computational domain

表2 計算域的大小Tab.2 Size of computational domain

2.2.2 外計算域網格 區(qū)域1的網格基礎尺寸為0.5 m，應用尺寸為2 m，并設置為各向同性．為更好地響應破損船運動，添加嵌套網格為區(qū)域3，基礎尺寸為0.5 m，其面網格最小尺寸為0.05 m，最大尺寸為2 m，并設置為各向同性．在嵌套網格和背景網格之間添加一個區(qū)域2為滑道域，其尺寸為0.25 m，作為區(qū)域1和區(qū)域3的過渡層，使網格變化遞進協調．

在網格控制中，體積增長率設為非常慢，面網格增長率為1.1，表面曲率為72，最小接近修復值為0.01．經以上設置，取該網格設置方案y向截面，如圖6所示，各區(qū)域網格過渡較好．

經網格檢測，在面網格中，高質量面網格數量為8 329 281，占比為100%；體網格中，較好質量體網格數量為18 464，占比為0.222%，高質量體網格數量為8 310 818，占比為99.778%．

2.2.3 波浪區(qū)網格 該船現作業(yè)區(qū)域為渤海范圍，結合渤海海域波浪特征[19]，采用以下3種不同的波浪形態(tài)：Wave 1(靜水)、Wave 2(常規(guī))、Wave 3(惡劣)，分別對破損失速養(yǎng)殖看護船DLHT 892 L288進行穩(wěn)性分析，其波浪參數見表3．

表3 數值模擬中入射波的參數Tab.3 Parameters of incident waves in the numerical simulation

為保證所加載波浪的準確性，采用二維潰壩方式對Wave 2和Wave 3進行了無風條件下波浪形態(tài)檢驗，網格尺寸的長度分別為對應波長的1/80，高度為對應波高的1/20；時間步Ts為0.01 s(時間步≤P/2.4n，其中P為波周期，n為用于劃分波長的段數，為方便與破損船的同步流固耦合運動，取0.01 s)．

所加載的波浪運動的殘差圖如圖7所示，其殘差值ξ不大，趨于穩(wěn)定，有較好的收斂性，該波浪網格方案合理可用．

(a)Wave 2

模擬波所呈現出來的形態(tài)分別如圖8和9所示，其模擬波浪的形態(tài)與阻尼消波情況較好．

圖8 阻尼條件下模擬Wave 2的波浪形態(tài)Fig.8 Wave shape of simulated Wave 2 under damping condition

圖9 阻尼條件下模擬Wave 3的波浪形態(tài)Fig.9 Wave shape of simulated Wave 3 under damping condition

3 模擬結果與討論

3.1 靜水狀態(tài)下破艙失速狀態(tài)

Wave 1海況下的下模擬網格運動殘差分析如圖10所示．其殘差值較小，趨勢具有較好的穩(wěn)定性，該網格劃分具有較好的合理性并滿足計算需要．

圖10 Wave 1條件下的殘差圖Fig.10 Residual plot under Wave 1 condition

在Wave 1條件下，外域湍流由破口涌入破損艙室，在考慮進艙水晃蕩影響下，其進水過程如圖11所示．

(a)t=0.01 s

Wave 1條件下，破口側的舷側壓力分布如圖12所示．

(a)t=0.01 s

Wave 1條件下，相對于初始狀態(tài)，受損船DLHT 892 L288的六自由度隨時間變化如圖13所示．

(a)橫搖

Wave 1條件下的破口流量如圖14所示．

圖14 Wave 1條件下破口質量流量隨時間變化Fig.14 The mass flow rate changes of the break under Wave 1 condition with time

在Wave 1條件下，在進水階段的0.56 s時，從破口涌入的水噴射到艙壁，對艙壁產生沖擊并開始形成新的自由液面，水流沖擊和新的湍流運動使靜止破損船產生橫搖并在10.50 s達到最大，隨后開始衰減．破口進水在0.81 s趨于穩(wěn)定，艙內進水形成渦流，初期湍流的涌入對船舶的六自由度影響也趨于穩(wěn)定，隨后產生規(guī)律穩(wěn)定的響應．破口舷側的壓力分布穩(wěn)定．船尾為壓強變化最大區(qū)域，在圖12中，從初始的18 459 Pa在20 s內變化到19 357 Pa，變化幅度很小，對船體影響不大．破損船DLHT 892 L288在靜水條件下的六自由度變化很小，橫搖響應周期較長，舷側壓力沒有較大變化，此狀態(tài)下的船舶穩(wěn)性較好．

3.2 Wave 2波浪狀態(tài)下的破艙失速狀態(tài)

Wave 2海況下的下模擬網格運動殘差分析如圖15所示．其殘差值較小，趨勢具有較好的穩(wěn)定性，該網格劃分具有較好的合理性并滿足計算需要．

圖15 Wave 2條件下的殘差圖Fig.15 Residual plot under Wave 2 condition

假定船初始位于Wave 2的波峰處，此時橫波在船處浪位最高，其船隨波浪產生的剛體運動變化如圖16所示．Wave 2波浪載況下破口側舷側壓力分布如圖17所示．Wave 2工況下，破口涌入的海水響應比Wave 1條件下要強很多，故在此捕捉進水的瞬時狀態(tài)，在考慮進艙水的晃蕩影響下，進水過程的瞬時變化如圖18所示．

(a)t=0.50 s

(a)t=0.50 s

(a)t=0.01 s

Wave 2條件下，相對于初始狀態(tài)，受損船DLHT 892 L288的六自由度隨時間變化如圖19所示．

Wave 2條件下的破口質量流量如圖20所示．

圖20 Wave 2條件下破口質量流量隨時間變化Fig.20 The mass flow rate changes of the break under Wave 2 condition with time

在Wave 2條件下，破艙管涌形成蘑菇狀射流沖擊艙壁，結合波浪作用，其艏搖在1.30 s達到6.8°，橫搖在第一周期達到7.0°，并在持續(xù)沖擊狀態(tài)中于8.00 s前后達到15.0°的傾斜位置，隨后逐漸收斂．縱蕩以及橫蕩在波浪與破損船的耦合作用下產生了較大的變化．其壓力最大處位于龍骨，隨縱搖前后移動，最大壓強變化區(qū)間為22 158～15 979 Pa，如圖17所示，變化幅度相對平緩，無中拱中垂現象．此條件下破艙后的六自由度變化雖然較大，但船舶無傾覆趨勢，最大橫搖為最大橫傾角一半，穩(wěn)性較好，比較安全．

(a)橫搖

3.3 Wave 3波浪狀態(tài)下的破艙失速狀態(tài)

Wave 3海況下的下模擬網格運動殘差分析如圖21所示．其殘差值較小，趨勢具有較好的穩(wěn)定性，該網格劃分具有較好的合理性并滿足計算需要．

圖21 Wave 3條件下的殘差圖Fig.21 Residual plot under Wave 3 condition

假定船初始位于Wave 3的波峰處，此時橫波在船處浪位最高，其船隨波浪剛體運動變化如圖22所示．

Wave 3波浪載況下破口側舷側壓力分布如圖23所示．

在Wave 3海況下，破口涌入的海水響應比Wave 2海況下的要劇烈很多，故捕捉其短時間內進水的瞬時狀況，在考慮進艙水的晃蕩影響下，進水過程的瞬時變化如圖24所示．

(a)t=0.20 s

(a)t=0.20 s

(a)t=0.01 s

Wave 3條件下，相對于初始狀態(tài)，受損船DLHT 892 L288的六自由度隨時間變化如圖25所示．

Wave 3條件下的破口質量流量如圖26所示．

(a)橫搖

圖26 Wave 3條件下破口質量流量隨時間變化Fig.26 The mass flow rate changes of the break under Wave 3 condition with time

在Wave 3條件下，破口處初始質量流量高達42 000 kg/s，破艙管涌對艙壁產生劇烈沖擊，新自由液面的晃蕩再結合外域波浪的推動，使其橫搖在1.10 s達到75°，縱搖在3.10 s達到12°，并且產生了較為嚴重的甲板上浪情況．雖然產生了巨大的橫傾后趨于收斂，這得益于艏搖的快速轉變，使得破損船舷側受波浪沖擊的相對面積大幅度減小，但如此夸張的橫傾角已經遠超該船的最大橫傾角，依賴波浪作用才不至傾覆．由于出現的艏搖運動使該船從遭受橫浪影響轉為遭受縱浪影響，且在圖23(b)中，龍骨中段最大壓強為40 180.0 Pa，兩端最小處為8 042.7 Pa，容易造成中拱現象．運動過程中，龍骨處最大壓強變化區(qū)間為40 180.0～20 810.0 Pa，該變化劇烈且強度大，容易給船體結構帶來一定損傷．此條件下破艙后的六自由度變化劇烈，并且伴隨較大的甲板上浪和中拱現象，雖然不至于傾覆，但是極度危險．

4 結 語

沿海漁船作為高發(fā)事故群體，其破艙后失穩(wěn)容易造成船體以外的各項損失．通過對比該船在Wave 1、Wave 2、Wave 3條件下的六自由度變化，可以發(fā)現在進水階段的過渡相中，破損船會向破口側傾斜，破艙管涌沖擊艙壁使船反向側傾，憑借復原力矩產生橫搖運動，其運動劇烈程度與破艙管涌的強弱相關．在初始進水流量中Wave 1、Wave 2、Wave 3流量比約為1∶20∶40，結合風浪荷載產生的最大橫搖比約為1∶10∶30，后者甚至達到了夸張的75°，按國際拖船大會的建議[20]，該船已經屬于傾覆狀態(tài)．幸運的是該破損艙室位于船頭，通過觀察舷側壓強分布變化，發(fā)現產生的艏搖運動在波浪的推動下使船頭朝向改變?yōu)轱L浪前進方向，進而較大地減少了橫搖運動的幅度．

綜上所述，沿海漁船的靜穩(wěn)性對于作業(yè)時破艙進水后物理行為的參考意義不大．考慮船上人員對于進水的第一階段難以及時反應，發(fā)現船艙破損時，進水階段往往處于漸進相，此時可通過調轉船頭方向，使船盡可能處于縱波狀態(tài)，以減小風浪和破艙管涌帶來的影響．

對于漁船艙室，盡量設置在漁船中橫剖線兩側，當發(fā)生破艙時，借助破艙進水的湍流沖擊和波浪耦合作用，推動破損船盡快從橫波狀態(tài)變?yōu)榭v波狀態(tài)，以減小橫搖幅度，減少損失．