鄒友家，沈 淳，奚祥英（上海海事大學 商船學院，上海 0035；武漢理工大學 管理學院，武漢430063）

鎳礦運輸船沉沒機理的數值模擬分析

鄒友家1，沈 淳1，奚祥英2
（1上海海事大學 商船學院，上海 200135；2武漢理工大學 管理學院，武漢430063）

近幾年來，多條裝運紅土鎳礦的船舶在南中國海傾覆沉沒，引起業(yè)界的震驚。紅土鎳礦是一種非常特殊的礦物質，在運輸過程中，如果船舶遭遇風浪，或船舶本身長時間的振動，都會導致該礦石的表面液化。液化后的混合液體具有較強粘性，可以形成自由液面。但由于現存的用于計算船舶自由液面的方法不適合于具有較大粘性的液體，因而無法解釋運輸該種貨物船舶傾覆沉沒的原因。所以，該文采用數學建模的方法，對運輸紅土鎳礦的船舶進行了模擬計算。得出的結論是當流體的慣性力大于流體的粘性剪力時，粘性流體對船舶的作用力矩隨時間的變化總趨勢是逐漸增大。當流體產生的橫向慣性力遠大于流體的剪切阻力時，流體將會掙脫剪切阻力的束縛，迅速滑向船舶的一側，使船舶的橫穩(wěn)性迅速消失，最終導致船舶很快傾覆沉沒。

紅土鎳礦；粘性液體；液化；粘性剪力崩潰；傾覆沉沒

0 引 言

隨著我國經濟的高速發(fā)展，國內對鎳礦的需求急劇增加，不可避免地導致了鎳礦運輸船的迅速增加。但是，自從2010年以來，相繼有多條載重噸在5萬左右的鎳礦運輸船在南中國海附近翻沉，造成了重大的人員傷亡與財產損失，引起了國際航運界的高度關注。

研究表明，這些事故與船舶所承載的貨物的特殊屬性有關。含有高氧化鐵的粘土類土質的鎳礦（俗稱紅土鎳礦）所含水分達到30%～40%，這種鎳礦物質液化后形成的流體具有較大的粘性[1]。傳統(tǒng)的船舶穩(wěn)性理論認為，粘性流體對船舶的穩(wěn)性具有正面影響，亦即當船舶向一舷發(fā)生傾斜時，粘性流體可以產生阻力，阻止船舶繼續(xù)向該側傾斜。因而，所有有關船舶穩(wěn)性理論的教科書都只研究非粘性流體對船舶穩(wěn)性產生的影響，即船舶的自由液面理論。但是，隨著多條鎳礦運輸船的翻沉，出現了船舶的自由液面理論無法解釋的現象。因此，有必要通過對船舶翻沉的過程建立數學模型，動態(tài)數值模擬鎳礦運輸船的沉沒，從而可能找到粘性液貨的粘性系數與船舶的橫搖周期、振幅等變量之間的關系。通過對鎳礦運輸船舶的翻沉過程進行數值模擬，也許可以找到鎳礦運輸船沉沒的真正原因。越來越多的研究表明，粘性流體可能并不總是對船舶的穩(wěn)性具有正面影響。

1 數值模擬研究

通過對海難事故船舶的調查了解得知，這些船舶有些共同點：出事船舶都是滿載紅土鎳礦，船上所有導航設備、主機副機及舵設備在當時都處于良好狀態(tài)，出事海域當時都有6～8級的橫風及2～3 m左右的橫浪，但傾覆沉沒的船舶狀況有好有壞，總的說來，大部分船況較好，有的船甚至是剛下水不到一年的新船，這就基本上可以排除是由于船舶縱向強度不夠造成船舶瞬間斷裂的可能性。

根據船舶沉沒之前的最后通訊記錄及幸存者事后回憶，船舶在沉沒之前都向一側發(fā)生了傾斜，隨即傾覆沉入海底。由此可以推斷，船舶的穩(wěn)性消失可能是造成船舶遇難的主要原因。那么，散裝礦石船的重心通常都是非常低的，其穩(wěn)性為什么會一下子消失呢？

1.1 自由液面理論

紅土鎳礦粉是一種被水分凝聚著的礦粉。當含水量超過35.7%時，遇到船舶搖擺和振動都會發(fā)生水份析出礦體表面的現象，稱作礦粉的“液化流動”（Liquefaction）。初始產生礦粉液化現象的含水量，被稱作流動水分點FMP（Flow Moisture Point）。液化析出的水分會聚集在紅土鎳礦的表面，形成自由液面，從而對船舶的初穩(wěn)性高度GM值產生負面影響。國內有研究人員利用自由液面對穩(wěn)性影響的計算公式對穩(wěn)性損失進行了計算。貨物表面自由液面對穩(wěn)性影響的計算公式為：LB（3×）γ×2/12×D。式中：貨艙長為L，寬為B，排水量為D，液體密度為γ。以一條排水量為21 000 t，5個貨艙的散貨船為例來進行穩(wěn)性損失的計算，穩(wěn)性變化值為：ΔGMi=-0.38 m。

研究人員假定船舶開航時GM為0．4 m，那么船舶在風浪中航行，紅土鎳礦粉游離出水份，導致穩(wěn)性損失達0.38 m，可能會出現初穩(wěn)性高度為負值，在外力的作用下可能會導致船舶向一側傾斜，以致傾覆[2]。

但是，以上的分析方法值得商榷。一般說來，運輸散裝鎳礦石的船舶，其貨物堆裝高度都在大艙深度的一半以下，保守地估算，初穩(wěn)性高度值一般都會大于1 m，正常情況下是1.5-2.5 m，有的甚至更大[3]。因此，由于自由液面的存在而使船舶的初穩(wěn)性高度損失0.38 m時，不太可能使船舶的穩(wěn)性高度變?yōu)榱慊蜇撝?。那么，船舶傾覆沉沒又該如何解釋呢？

圖1 紅土鎳礦在液化前與液化后的對比Fig.1 Lateritic nickel ore before and after liquefaction

1.2 數值模擬的建模

通過大量的調查研究，我們發(fā)現紅土鎳礦是一種非常特殊的貨物。當其水份含量超過35.7%時（還需進一步測試數據的可靠性），如果船舶長時間振動或在風浪的作用下就會發(fā)生液化現象，貨物的表面就會析出“一層水”。由于紅土鎳礦是一種有較強粘性的物質，析出的水份不會通過貨物流到艙底，而是附著在貨物的表面。而且，析出的水份不完全是水，而是一種水與腐泥土的紅色混合流體，如圖1。根據各礦區(qū)產出的礦石成分及化學組成不同，混合液體的粘度有很大的不同。

目前，計算自由液面對穩(wěn)性的影響是有條件的，即液體的粘度可以忽略不計。當液體的粘度不能忽略的時候，現存的計算自由液面對穩(wěn)性高度的影響的公式就不再適用。這也許可以解釋為什么上述算法不能解釋船舶穩(wěn)性消失的原因。

當船舶裝載粘性液貨在風浪中航行時，粘性液貨必然產生晃蕩現象。大幅晃蕩時，自由液面的運動具有強烈的非線性。圖2所示為外力作用下三維貨艙坐標系統(tǒng)。

二維不可壓縮帶自由表面的粘性流體運動，其控制方程分別由連續(xù)方程[4]

圖2 外力作用下3維貨艙坐標系統(tǒng)Fig.2 Coordinate systems for a three-dimensional tank under forced excitation

式中：a為質點o相對于固定坐標系的平均加速度，ω為船舶橫搖的角速度，rG為從o點出發(fā)的矢徑。

盡管當今有多種方法處理邊界問題，但本文主要采用數值方法解決流體在艙內運動產生的晃蕩。由于流體在艙內的晃蕩表現為強烈的非線性，特別是接近共振區(qū)時。這些非線性包括破碎波、飛濺、噴射流以及沖擊波等。若要將這些因素全部考慮進去是很困難的，而且，有些因素對本研究沒有重大的意義。因而，流體的自由液面邊界可以認為是一個單值函數。運動自由液面邊界條件可以寫成[5]：

式中：η為自由液面高度，u為速度矢量。

1.3 數值方法

由于有些局部的流體不具有普遍意義，因而我們可以使用一種快速有效的數值方法來模擬流體的運動[6]。國外學者Kim利用由Hirt提供的SOLA方法對不帶斜面的三維容器（四角均為直角）成功進行了模擬。而我們則利用同樣的方法來處理帶有斜面（靠近艙底的一段側面艙壁）的三維容器。具體的方法如下：

圖3 棱型容器笛卡爾網格Fig.3 Cartesian grids for a prismatic tank with internal members

首先，將容器的體積離散為三維網格（I，J，K），采用笛卡爾交錯網格，速度分量定義在網格的邊界，而壓力則在網格的中心，如圖3。在SOLA方法中，瞬時壓力和速度的計算采用迭代法。迭代法的第一步是計算偽速度[5]：

式中：Δt是時間分量，＾表示離散梯度，下標i，j，k表示在x，y，z方向的網格指數，上標n表示第n步。將（5）式代入連續(xù)性方程，我們就可以計算出壓力改正值：

式中：Δ（x，y，）

z為在x，y，z方向的空間分量，因此，實際速度可以寫成：

利用（6）式以及（7）式進行的壓力與速度改正應該在每一個網格中都要做。改正的速度可以代入（6）式進行下一步的壓力改正，這種迭代方法要一直運用下去，直到壓力改正值在容許的誤差范圍內。

上式對有限差分中的對流項靈敏度較高。Hirt提出了將二階中心差分法與一階逆風差分法結合的方法[7]。采用分量參數α，則可以合成為：

很明顯，在一階逆風差分中，如果參數α增大，例如增加到接近1時，將會造成數值衰減。方程（4）的離散形式可以寫成：

式中：wη，uη，vη為自由液面的速度分量，可以通過內插求得。液面坡度項可以用方程（8）的有限差分形式進行估算。特別要注意的是，為了克服鋸齒波對方程穩(wěn)定性造成的影響，需要使用數值濾波[9]。

由Longuet Higgins以及Cokelet提出的5點式公式對解決水波問題非常有用，對本研究也仍然適用[10]：

通過運用Miyata et al提出的3維不規(guī)則星方法[11]，可以滿足方程（2）中自由液面交叉點各點的動態(tài)壓力。該方法是基于6個相鄰點的泰勒壓力級數展開式，如圖4。

在每一個自由液面網格，其中心壓力可以寫成：

圖4 3維不規(guī)則星方法的6點示意圖Fig.4 Six points for the three-dimensional irregular-star method

式中：li是到第i個相鄰點的距離。如果第i個相鄰點位于自由液面上，則pi就應該等于大氣壓力。對于單值自由液面的剖面，可以認為5點中的最大值等于大氣壓力[10]。

1.4 模擬船舶橫搖的影響

假定波浪對船舶產生一個周期性的擾動，其擾動力矩可以表示為[11-14]：

船舶運動的方程為：

2

式中：J為船舶相對于縱軸的船舶慣性及附加質量慣性，φ是船舶的橫搖角，B是阻尼系數，C是復原力矩，Mω是每一種波浪頻率引起的擾動力矩，φ是擾動相位，Mt為流體的力矩振幅，φ′是相位滯后角。

2 數值模擬結果與分析

根據以上建立的數學模型，我們可以動態(tài)地模擬出粘性液體對艙壁的壓強隨時間的變化關系曲線圖?，F在以一條載重噸為5萬噸的散貨船為例（事故船的載重噸均在5萬噸左右），滿載紅土鎳礦，航速12 Kn，吃水11 m，風力6～8級，左舷受風，橫浪，浪高2～3 m，橫搖角6°，周期10 s。從圖5可以看出，在流體粘度很小的情況下，艙壁所受壓強呈現準周期變化。

船舶在外力的作用下，一定會產生一個外力矩，這個外力矩由風浪作用在船舶上的力矩與粘性流體作用在船舶上的力矩組成。由圖6可以看出，風浪作用在船舶上的力矩呈現準周期性變化。

圖5 艙壁所受壓強與時間的關系Fig.5 Pressure on ship’s wall vs time

圖6 風浪作用在船舶上的力矩與時間的關系Fig.6 Moments of winds and wave acting on ships vs period

圖7 自由液面位移與外激勵頻率的關系Fig.7 Free surface displacement versus excitation frequency

圖8 速度矢量與壓力場在3個時刻的模擬圖：h/B=0.026，s/B=0.1，Te=10 sFig.8 Velocity vector and pressure contour at three time steps:h/B=0.026,s/B=0.1,Te=10 s

模擬實驗表明，當激勵頻率與船舶的固有頻率接近或相等時（f=0.093 1 rad/s），自由液面產生的位移（液面離開原來的水平位置）最大。外力（風浪）所做的功主要消耗在使船舶產生橫搖運動以及使自由液面勢能增加上（其它的消耗暫時忽略不計）[11]，見圖7。

由于鎳礦的比重大，完全液化后其堆裝高度通常只占大艙高度的五分之一左右，因而做低液位模擬更符合實際情況。圖8是對液位高度、速度場和壓力場的數值模擬。h為液面高度，B為船寬，s為橫蕩位移，Te為激勵周期。

但是，粘性流體（本例為鎳礦混合液體）在外力的激勵下，其力矩隨著時間的變化與風浪對船舶的作用力矩隨時間的變化有很大的不同。風浪對船舶的作用力矩隨時間的變化具有準周期性。但粘性流體對船舶的作用力矩可以分為2個階段：當慣性力小于粘性剪力時，粘性流體對船舶的作用力矩也呈準周期變化；當慣性力大于粘性剪力時，粘性流體對船舶的作用力矩隨時間的變化總趨勢是逐漸增大。增大的幅度與流體的粘度系數、船舶的橫搖振幅及風浪周期有關。如圖9.當船舶的橫搖振幅增加時，流體產生的慣性力就會增加。而當流體產生的橫向慣性力大于流體的剪切阻力時，流體將會掙脫剪切阻力的束縛，迅速滑向船舶的一側。使得粘性流體作用在船舶上的力矩迅速增大，船舶的橫穩(wěn)性可能迅速消失，從而引起船舶的傾覆。

3 結 論

本文利用數學模型方法對載運粘性流體的船舶進行了穩(wěn)性分析與計算。得出的結論是：當流體的慣性力大于流體的粘性剪力時，粘性流體對船舶的作用力矩隨時間的變化總趨勢是逐漸增大。當流體產生的橫向慣性力遠大于流體的剪切阻力時，流體將會掙脫剪切阻力的束縛，迅速滑向船舶的一側，使船舶的橫穩(wěn)性迅速消失，最終導致船舶很快傾覆沉沒。

圖9 粘性流體作用在船舶上的力矩與時間的關系（當慣性力大于粘性剪力時）Fig.9 Moments of viscous liquid acting on ships vs period(when inertial force is greater than shear stress)

上述結論還需要做進一步的研究，還需要進行更多的實驗，以便確定較為準確的船舶橫搖振幅，周期與流體粘性系數的關系。另外，風浪對船舶的作用具有非對稱性，往往是從船舶的一側作周期性的打擊，而不是兩側交替打擊，這樣，我們的數學模型還需要進一步修改和完善，以便與實際觀測的結果相一致。

[1]崔建輝,王建平.海運紅土鎳礦的安全運輸[J].航海技術,2011,6:30-33. Cui Jianhui,Wang Jianping.Safe transport for lateritic nickel ores[J].Marine Technology,2011,6:30-33.

[2]雷 海.紅土鎳礦粉的安全運輸[J].航海技術,2011,1:27-28. Lei Hai.Safe transport for lateritic nickel ore[J].Marine Technology,2011,1:27-28.

[3]李榮輝,徐邦禎.基于臨界橫傾角的散裝谷物船舶穩(wěn)性衡準[J].航海技術,2009,6:26-27. Li Ronghui,Xu Bangzhen.Stability norm for grain carriers based on critical heeling angle[J].Marine Technology,2009,6: 26-27.

[4]朱仁慶,方智勇,張照鋼,等.Level-set法預報液艙內劇烈晃蕩引起的沖擊壓強[J].船舶力學,2008,12(6):345-451. Zhu Renqing,Fang Zhiyong,Zhang Zhaogang,et al.Level-set method for predicting impact pressure induced by violent sloshing in a tank[J].Journal of Ship Mechanics,2008,12(3):345-451.

[5]Kima Yonghwan,Shin Yungsup,Lee Kwanghyun.Numerical study on slosh-induced impact pressures on three-dimensional prismatic tanks[J].Applied Ocean Research,2004,26:213-226.

[6]賀俊松,陳 震,肖 熙.五體船大傾角穩(wěn)性研究[J].船舶工程,2007,29(5):12-15. He Junsong,Chen Zhen,Xiao Xi.Study on the stability at large angle of the pentamra[J].Ship Engineering,2007,29(5): 12-15.

[7]Hirt C W,Nicholas B D,Romero N C.SOLA-a numerical solution algorithm for transient fluid flows[R].Los Alamos Scientific Laboratory Report,1975,LA-5852.

[8]魯 江,馬 坤,黃武剛.斜浪長峰不規(guī)則波中船舶復原力變化計算[J].中國造船,2010,51(3):11-17. Lu Jiang,Huang Wugang.Calculations on the righting forces of ships in irregular waves with long inclined sea[J].China Shipbuilding,2010,51(3):11-17.

[9]丁仕風,唐文勇,張圣坤.速度突變情況下液化天然氣船液艙內晃蕩問題的仿真[J].上海交通大學學報,2008,42(6): 920-923. Ding Shifeng,Tang Wenyong,Zhang Shengkun.Simulation of liquid sloshing in the cabin caused by liquefied natural gas ship’S variable speed[J].Journal of Shanghai Jiaotong University,2008,42(6):920-923.

[10]Longuet-Higgins M S,Cokelet E D.The deformation of steep surface waves on water.I.A numerical method of computation[J].Proc R Soc,1976,A350:1-26.

[11]Miyata H,Inui T,Kajitani H.Free surface shockwaves around ships and their effects on ship resistance[J].J Soc Naval Architect Jpn,1980,147:1-9.

[12]雒高龍,任慧龍.超大型LNG船晃蕩壓力直接計算法[J].造船技術,2010,297(5):7-10. Luo Gaolong,Ren Huilong.Direct calculations of sloshing forces for the super LNG[J].Shipbuilding Technology,2010, 297(5):7-10.

[13]劉新立,嚴仁軍,史 政,等.基于MSC.Dytran的液艙晃蕩分析[J].船海工程,2010,39(2):44-47. Liu Xinli,Yan Renjun,Shi Zheng,et al.Sloshing analysis of liquid tanks based on MSC.Dytran[J].Ocean Engineering, 2010,39(2):44-47.

[14]陸志妹,范佘明.船舶液艙晃蕩研究進展[J].上海造船,2010,2:14-30. Lu Zhimei,Fan Sheming.Outlooks of sloshing research[J].Shanghai Shipbuilding,2010,2:14-30.

Numerical simulations on the mechanisms of the capsizing of bulk nickel ore carrier

ZOU You-jia1,SHEN Chun1,XI Xiang-ying2
(1 Shanghai Maritime University,Shanghai 200135,China;2 Wuhan University of Technology,Wuhan 430063,China)

A number of bulk nickel ore carriers capsizing in South China Sea have made shipping industry shocked.Clay lateritic nickel ore is characterized by liquefaction on the surface under the influence of waves or long time vibrations during transportation.The blending liquid after liquefied is highly viscous and can generate free surface.However,the existing methods based on the calculation of GM values fail in explaining the capsizing of the carriers during transportation.In this paper,the mathematical modeling is adopted to make a thorough analysis to the mechanisms of capsizing.The results show that the imbalance of the shear force leads to the free movements of the liquid cargoes as the inertial force of the liquid cargoes is greater than that of the shear stress,and thus caused a loss of stability.Finally,the increasing external moments result in the capsizing.

clay lateritic nickel ore;viscous liquid;liquefaction;shear stress;collapse;capsizing

U661.32

A

10.3969/j.issn.1007-7294.2015.08.004

1007-7294（2015）08-0912-07

2015-01-28

鄒友家（1964-），男，博士，高級船長，E-mail:marscar@126.com；

沈 淳（1966-），男，副教授；奚祥英（1966-），女，副教授。