基于無網格MPS方法數值分析方艙破艙進水問題

盧澤宇，萬德成

(上海交通大學 船海計算水動力學研究中心(CMHL)船舶海洋與建筑工程學院，上海 200240)

伴隨著經濟全球化和地區貿易自由化，航運業也得到快速發展，然而多種多樣的海難事故，例如著火或爆炸、進水、碰撞、擱淺、傾斜或傾覆、沉沒、失控及漂浮等，造成了嚴重的經濟損失，對人員的生命安全也造成了嚴重威脅[1]。海難事故導致的受損船舶水密艙室進水也一直是行業關注的焦點，眾多學者對其進行了研究。受損船舶的水密艙室進水可能會嚴重影響船舶的穩性并威脅其安全。水的快速流入可能導致大振幅的船舶運動從而導致傾覆，研究的目的是對進水過程中二維受損艙段的運動特性進行數值模擬分析。

按照時間順序，船舶受損后，如圖1所示，進水過程主要分為3個階段，瞬時進水階段、持續性進水階段和穩定階段[2]。在大多數損壞情況下，瞬時進水階段會引起受損船舶大幅度橫搖，但持續時間較短，這使得該階段往往是最危險的階段；穩定階段持續時間較長，但船舶橫搖過程較平緩；持續性進水階段需要更長的時間才能使受損船舶達到新的平衡狀態。

圖1 進水過程的主要階段

因此，目前的研究大多數集中在持續性進水階段[3]。由于靜態方法過于簡單，無法應用于持續性進水階段，許多學者采用準靜態方法進行持續性進水階段的模擬[4]。Papanikolaou等[5]和 Jasionowski[6]考慮了進水的影響。為了研究破艙的過程，Palazzi和De Kat[7]將空氣和水看作是在一定范圍內機械運動的集中質量塊，一定程度上，這種假設更接近于受損船舶運動的模型試驗結果。Ruponen[8]通過時域浸水分析，研究了非水密結構對典型現代大型客船持續進水和傾覆時間的影響。Lee[9]還提出了通風隔間的新模型和蓄能器模型，以研究復雜的開口系統對持續進水過程中水流和空氣的影響。以上研究都基于準靜態假設，計算成本相對較低。然而，進一步的研究表明，艙內的流動確實會影響船舶的瞬態運動特性，使用簡化的準靜態模型難以對問題進行綜合考慮。因此，研究人員開始探索其他方法。

Manderbacka等[10]提出了一種全耦合方法來研究船舶對瞬時進水的響應，其中考慮了傳遞動量的進水入口，并通過進水船舶的橫搖阻尼和瞬時進水試驗數據驗證了該方法。Acanfora和Cirillo[11]通過將準靜態方法與完全耦合方法耦合開發了一種中間方法，不再假設自由表面是水平的，并且平面的傾斜角取決于船的水平加速度，將模擬數據與駁船模型的試驗測量值進行比較并得到了驗證。Manderbacka和Ruponen[12]通過比較全動態和準靜態洪水模擬，對瞬時進水的水動力學進行了研究，通過改變艙壁開口的大小來評估進水射流對船舶響應的影響。一些研究人員依靠模型試驗來評估受損船舶運動的高度非線性進水過程。Zhang等[13]的試驗考慮了兩種帶有舷側開口和底部開口的船舶模型，以證明損壞對船舶失穩和沉沒時間的影響。Gu等[14]試驗研究了完整和損壞的戰艦 DTMB-5415 在波束中的性能，表明海水的進出及其與船舶運動的相互作用顯著影響船舶上受到的沖擊載荷。

當前試驗方法存在局限性，例如難以對整個流入區域的數據進行全面監測，試驗準備成本高且復雜等。近年來，計算機計算能力的提高使得船舶水動力響應的數值模擬變得容易。Gao和Vassalos[15]使用基于雷諾平均N-S方程(RANS)的計算流體動力學(CFD)求解器和基于自由表面的流體體積(VOF)方法來研究液艙晃蕩和進水對損壞船舶水動力學的影響。Cao等[16]開發了一個三維并行光滑粒子流體動力學(SPH)程序來模擬具有舷側和底部開口的模型下沉過程，以研究船舶的穩定性和下沉時間。為了揭示破艙水動力學與波浪中受損船舶之間的相互作用，Gao和Vassalos[17]使用結合了耐波求解器和 Navier-Stokes 求解器的集成方法來模擬滾裝渡輪在規則波中的運動，發現受損船舶的自然橫搖頻率變化和橫搖響應降低主要是由浸水艙室內的水晃動所致。

近些年的研究主要集中在持續性進水階段，很少有研究考慮過瞬時進水。事實上，瞬時進水期間的大角度橫搖可能會導致受損船舶傾覆。因此，為了評估這一瞬時進水階段，對受損船舶進行可靠的數值模擬是必要的。Manderbacka等[10]提出了一種用于受損船舶運動的非線性時域仿真方法，并研究了兩種不同的初始穩定性條件、兩種不同的艙室布局以及黏性阻尼的影響。此外，現有研究還表明，不應忽視破艙進水流量的影響。而使用粒子法的優點是可以準確地捕捉自由面的劇烈變形，因此，已被廣泛用于模擬液艙晃蕩、潰壩流、物體入水、甲板上浪、波物相互作用等劇烈流動問題。張雨新等[18-19]采用MLParticle-SJTU程序對二維和三維液艙晃蕩問題進行了模擬，計算的壁面砰擊壓力和自由面波形與試驗和其他方法結果吻合。Zhang等[20]使用改進的MPS方法對二維損壞船段進行了不同開口和擋板的數值模擬，并對進水量等物理量進行了分析。田鑫和萬德成[21]采用MLParticle-SJTU程序對逐漸潰壩的問題進行數值模擬，研究了壩體不同失效時間對潰壩流動以及砰擊壓力的影響。Xie等[22]基于自主開發的軟件MPSGPU-SJTU研究了方形液艙內晃蕩的三維效應。Wen等[23]在改進的MPS方法基礎上，進一步發展形成了具有高精度和高穩定性的多相流MPS方法，并結合多CPU 并行加速技術和GPU并行加速技術，開發了擁有獨立自主知識產權的高性能多相流求解器軟件MMPS，隨后，成功將文中方法和求解器應用于多個典型多相流問題，包括含復雜界面的氣泡流問題，船海工程中內孤立波問題，以及劇烈兩相晃蕩和潰壩問題。

本文簡要介紹了改進的MPS方法數值模型，MPS方法通過與SPH方法和VOF方法進行對比得到驗證，最后給出數值模擬模型的參數和不同的橫截面模型，比較了受損艙段截面各種模型的破艙進水情況。

1 數值方法

從核函數、壓力梯度模型、壓力泊松方程源項以及自由面判斷等方面入手，對MPS方法進行了一系列改進，并在單相流問題中取得了比較好的應用成果。這里將對MPS方法基本格式進行介紹。

1.1 控制方程

對于不可壓縮黏性流體，MPS方法控制方程包括連續性方程和N-S方程，其形式為：

(1)

(2)

式中：D/Dt表示物質導數，t為時間，ρ為流體密度，P為壓力，μ為動力黏性系數，取1.01×10-5，V為速度矢量，FV為體積力。可以注意到，在MPS方法中，動量方程左側的時間導數項以物質導數形式給出，因此有效避免了對流項帶來的數值耗散。

在無網格粒子類方法中，整個計算域通過對一系列空間粒子進行離散，并根據粒子間相互作用對上述控制方程進行求解。但由于此類方法中的粒子可自由運動，使粒子間拓撲關系變得不固定，因此需在每步計算中進行鄰居粒子搜尋和配對，并對不同粒子間的相互作用大小進行評估。為此，MPS方法中引入了核函數及粒子相互作用模型的概念。

1.2 MPS核函數

核函數(kernel function)在無網格粒子類方法中起到權重函數的作用。一般而言，當粒子間距離大于一定范圍時，核函數的值為0，即該對粒子之間不發生相互作用。隨著兩個粒子不斷靠近，核函數的值逐漸變大，粒子間的相互作用也變得更強，反之亦然。與SPH 方法不同，MPS方法中并不要求核函數光滑可導，因此核函數形式更簡單，也更多樣。張雨新和萬德成[18]提出一種無奇點的核函數，形式為：

(3)

式中：re為粒子作用半徑，r為相鄰兩個粒子的距離。通過Koshizuka和Oka[24]的數值研究，通常對拉普拉斯模型的作用半徑re1取4.01l0，對粒子數密度模型、梯度模型、散度模型和鄰居粒子分布函數模型的作用半徑re2取2.1l0，其中l0為初始粒子間距。

1.3 粒子相互作用模型

基于上述核函數，可對某一目標粒子受到的所有鄰居粒子作用進行加權平均，并推導得到一系列的粒子相互作用模型。利用這些模型，最終實現對控制方程中多種微分算子的離散和求解。MPS方法中所采用的粒子相互作用模型，主要包括梯度模型、散度模型和拉普拉斯模型，分別為：

(4)

(5)

(6)

式中：P為壓力，V為速度矢量，φ代表任意標量，D為空間維度，r為粒子位置矢量。n0為初始粒子數密度，用來表征流體粒子分布的疏密程度，如式(7)。λ是對用有限范圍的核函數來代替無限范圍的高斯函數帶來誤差的一種補償，可用式(8)計算得到。

(7)

(8)

1.4 壓力泊松方程

MPS方法的流場壓力是通過壓力泊松方程計算得到，Tanaka 和Masunaga[25]結合粒子數密度和速度散度，提出了一種混合源項的壓力泊松方程，如式(9)。

(9)

式中：Δt為時間步長，n*為臨時粒子數密度，V*為粒子臨時速度矢量。根據Lee等[26]的研究，γ通常取值為0.01～0.05。

1.5 固體壁面和自由面邊界條件

在泊松方程求解時，自由面粒子和第二類邊界粒子的壓力被賦值為0作為方程的邊界條件，因此需要正確判斷出自由面粒子。如圖2所示，對于單相流MPS方法，并不考慮氣相和液相之間的相互作用，沒有粒子表示氣相介質，因此原本的氣液兩相交界面就是流體相的自由液面，并且自由面附近流體粒子的鄰居粒子基本分布在一側。基于這種粒子分布的不對稱性，張雨新和萬德成[18]定義了一個衡量鄰居粒子不對稱性的矢量：

圖2 鄰居粒子分布

(10)

當粒子滿足式(11)，該粒子則被判斷為自由面粒子。

|F|>α|F|0

(11)

式中：|F|0為初始時刻自由面粒子的|F|值，α為常數。經過張雨新和萬德成[18]的數值研究，α建議取值0.9。

MPS方法采用布置多層粒子的策略作為物面邊界條件，在物面上布置一層粒子稱為第一類邊界粒子，然后在第一類邊界粒子之外再布置兩層第二類邊界粒子，如圖3所示。其中，第一類邊界粒子的壓力和流體粒子一起參與壓力泊松方程求解，第二類邊界粒子的壓力通過第一類邊界粒子和流體粒子的壓力插值得到。這種多層粒子布置方式克服了在物面附近發生積分域的截斷，引起數值結果失真問題，還避免了流體粒子在邊界粒子運動時穿透物面。

圖3 邊界粒子分布

2 結果與討論

2.1 數值模型驗證

為了證明模型的可行性，采用基于GPU的MPS方法來模擬二維方艙舷側上部開口的破艙進水過程，與Shen等[27]使用VOF和SPH方法的模擬結果進行比較。

如圖4所示，模型為邊長0.1 m的正方形艙，舷側上方有一個0.02 m的開口，在計算域內模擬靜水下的破艙進水，其參數如表1所示。

圖4 模型尺寸示意

表1 方艙參數

模擬過程如圖5所示，從上到下分別是VOF、SPH、MPS方法，從左到右分別是在0.1 s、0.2 s、0.3 s、0.4 s時刻的模擬。可以看出在各個時刻3種方法的模擬現象基本一致。

圖5 3種方法在不同時刻的模擬

將模擬結果與SPH方法對比如圖6所示，可以看出文中所用MPS方法在橫蕩和垂蕩上的加速度與SPH方法的結果擬合較好。

圖6 兩種方法的縱蕩和垂蕩模擬結果對比

2.2 橫截面模型

為了系統地研究不同破艙場景下的影響，文中模擬了受損艙段在破艙進水后的運動響應，在其舷側上的關鍵位置開孔來表示船體損壞。在受損船段內使用非穿透性擋板進行進一步的模擬。艙的尺寸在上一節中有介紹，開孔位置以及艙內擋板分布如圖7所示。初始粒子間距為0.05 mm，時間間隔為0.000 01 s。對于A1模型，粒子總數為328 738個，其中流體粒子289 750個，壁面粒子為13 500個，每個時間步的計算時間約為30 s。其他模型的粒子數量和計算時間大致相同。

圖7 橫截面模型的9種配置

2.3 開口位置對破艙進水的影響

圖8(a)繪制了模型A1、A2、A3的進水量N的時間歷程。可以看出，損壞孔的位置對進水速率的影響較為顯著。一般來說，水面以下距離水面較遠的受損孔洞會導致更大的進水量，這是由于距離水面越遠水壓越大，會迫使更多的海水涌入艙室。

圖8(b)可以看出進水導致舷側開口模型A1、A2和A3向損壞側水平移動。位移在達到峰值后，由于右側反射回來的流體沖擊左側艙壁導致位移回調。從圖8(c)可以看出，開口位置越低，下沉速率越快，與進水速率成正相關。

圖8(d)表明，開孔的位置顯著影響了橫搖運動。A1位置的開口由于壓力缺口導致艙左傾，而A2、A3開口由于入射水柱的沖擊力更強，進水量更大，導致初始階段產生右傾。同時可以看出開口位置更靠下部的A3比A2對右傾的影響更大，和進水量成正相關。

圖8 舷側不同開口位置模型的進水量、橫蕩、垂蕩和橫搖時域歷程

圖9展示了A2模型破艙進水的過程，從圖中可以看出，損壞開孔位置對流體的流動特性也有顯著的影響，并在艙室內產生渦流。

圖9 模型A2的舷側開口引起的進水過程屏幕截圖

2.4 水平擋板對破艙進水的影響

為了研究障礙物對受損艙段截面模型相關運動響應的影響，在船截面模型B1、B2、B3和C1、C2、C3內分別設置了非穿透性水平和垂直擋板，如圖7所示。

圖10(a)繪制了帶有水平擋板的模型進水量N的時間歷程。由圖可知，這些模型的進水量在同一時刻幾乎沒有差別。也就是說，流入速率之間的差異僅隨時間略微增加。說明水平擋板位置對進水量的影響有限。

圖10(b)、10(c)、10(d)分別繪制了帶有擋板的模型B1、B2、B3引起的橫蕩、垂蕩和橫搖運動時間歷程。如圖所示，在不同位置布置水平擋板對橫蕩和垂蕩的影響差異很小而對橫搖的影響差異很大，且隨著時間增加這種差異會越來越大。在遠離開口位置附近布置水平擋板對橫搖的改善明顯要優于開口位置附近布置水平擋板。在t=0.8 s時，上部水平擋板模型(B1)的橫搖角度與下部水平擋板模型(B2)相差近40%。

圖10 艙內不同位置水平擋板模型的進水量、橫蕩、垂蕩和橫搖時域歷程

2.5 垂直擋板對破艙進水的影響

圖11(a)繪制了帶有垂直擋板的模型進水量N的時間歷程，如圖所示，垂直擋板遠離開口位置的模型進水速率更快。產生這種現象的原因是垂直擋板阻止了初始入射流擴散到遠離損壞側的區域，導致更大的橫搖力矩和更劇烈的橫搖運動(如圖12所示)。

圖11(b)、11(c)、11(d)分別繪制了帶有擋板的模型C1、C2、C3引起的橫蕩、橫搖和垂蕩運動的時間歷程，如圖所示，在不同位置布置豎直擋板對橫蕩和垂蕩的影響差異很小而對橫搖的影響差異很大，且隨著時間推移這種差異會越來越大。垂直擋板影響了初始入射流的分散，導致更劇烈的橫搖運動。在t=0.8 s時，最右邊垂直擋板模型(C3)的橫搖角度與最左邊水平擋板模型(C1)相差近30°。

圖12展示了模型B1、B2、B3(帶水平擋板)和模型C1、C2和C3(帶垂直擋板)的破艙進水過程中在時間t=1.1 s時的壓力云圖。可以看出，垂直擋板對破艙內流動的影響大于水平擋板，表明擋板的位置和形式會對破艙后的運動響應有影響，特別是對橫搖運動產生很大影響。

圖12 帶擋板方艙破艙進水過程在t=1.1 s時刻的壓力云圖

3 結 論

本文數值模擬了靜水條件下沉沒過程中受損二維方艙艙段模型的進水及其相關運動。首先，為了滿足數值計算的需要，選用基于GPU加速技術的MPS方法。通過與以往其他數值結果的對比，證明該模型的準確性，適用于研究自由漂浮二維破損艙段的瞬時進水過程。其次，對該方艙的各種模型進行了數值模擬。這些模型由方艙艙段橫截面組成，其舷側不同位置有開孔，同時在內部布置水平和垂直的非穿透性擋板。通過數值模擬分析，得出以下主要結論：

1)驗證了MPS方法對船舶破艙進水數值模擬的適用性，研究表明，對于不同類型的船體截面模型，可以很好地預測船體截面的瞬時進水和相關運動響應。

2)受損孔的位置影響了破艙進水的進水量。開孔距靜水面以下越遠，進水量越大。

3)受損孔的位置也會影響模型的運動。開孔距水面以下越遠，橫蕩運動和垂蕩運動就越大，同時進水量也越大，進水速率更快。開孔是否全部浸沒在水面以下對于橫搖的方向有很大影響，同時浸沒的深度會影響橫搖的角度大小。

4)內部擋板顯著影響了模型的運動。垂直擋板比水平擋板對破艙進水的運動響應影響更大，尤其是橫搖運動。這主要是由于垂直擋板對入射流擴散的影響較大，導致艙內流體分布不均，間接增加了橫搖力矩。

本文的研究有助于了解不同位置的受損孔和擋板位置類型對船舶破艙進水過程中瞬時進水階段船舶運動特性的影響，為三維模擬的進行奠定基礎，對艙內的艙室布置有一定的指導意義。