畢東魁，王翔宇，孫樹政，陳睿童，郁 峰，尤婷婷

(1. 中國船級社，山東 威海 264200；2. 哈爾濱工程大學船舶工程學院，黑龍江 哈爾濱 150000；3. 煙臺哈爾濱工程大學研究院，山東 煙臺 264000；4. 中國船舶集團有限公司，上海 200011)

0 引 言

我國處于北半球，船運走北極航線必須經過邊緣冰帶（MIZ），這個區(qū)域處于極地和溫帶氣候系統的交匯處。該區(qū)域浮冰的大小從數米到數百米不等[1]。這些充滿隨機性的浮冰會對船舶的總體結構與人員生命安全造成嚴重威脅。所以，船舶在浮冰區(qū)域的冰阻力預測成為船舶安全保障的關鍵。

對于船舶浮冰阻力的研究，有實船試驗法、水池模型船試驗法、計算機仿真模擬等，綜合對比水池試驗法在耗費人力物力、試驗周期、結果精確性等方面有較高優(yōu)勢，也是現代學者研究冰載荷的一種較為認可的方式。陳昭煬[2]在拖曳水池中進行了模型冰船阻力試驗，國內首次應用試驗方法研究浮冰的形狀對冰阻力的影響，觀察浮冰的堆積和貼體情況，分析工況、浮冰運動、阻力之間的影響規(guī)律。Jeong[3]在冰池中建立了一條碎冰航道，并根據船舶在該碎冰航道中航行時拖曳力與推進阻力的關系，推導出帶有螺旋槳模型船的性能。隨著計算機的發(fā)展，可以做到對水池試驗進行準確仿真，不僅可以保證數據精準，而且有容易控制變量、靈活度高等優(yōu)勢。徐棟梁[4]在研究船舶舷側與浮冰的碰撞過程中考慮了溫度場對船體鋼材料的影響，總結溫度、速度、位置等因素對加強筋結構碰撞的影響，根據仿真結果擬合出相對應的經驗公式并驗證。

本文應用光滑粒子流體動力學方法(smoothed particle hydrodynamics, SPH)與有限單元法(finite element method，FEM)相耦合模擬冰水池中船與浮冰碰撞。SPH 模擬流體具有自適應的優(yōu)勢，本質上仍然是拉格朗日方法，但在處理高速碰撞中不會出現傳統FEM 中網格嚴重變形導致結果不精確的問題。SPH 粒子在耦合中被視為特殊的節(jié)點，控制參數為節(jié)點編號、質量以及空間位置，粒子之間自行耦合，與有限單元點面接觸耦合，可以較為準確做到冰-水-船耦合。

1 數值模型

SPH 的基本原理為拉格朗日法。將連續(xù)的海水與浮冰離散成諸多粒子質點，這些粒子擁有獨立的質量、位置坐標、速度、加速度等物理量。隨后求解整體的動力學方程，跟隨每個質點的各項物理參數隨時間的變化情況，最終綜合所有質點得到整個流場與冰體的物理情況。

SPH 方程構造首先采用核函數逼近的方法，將描述場的函數轉化成積分表達式，然后粒子近似，實現了對核函數近似積分表達式的粒子離散。

粒子近似式方程表達為：

式中，W為核函數。

核函數W由函 數 θ定義：

式中：d為空間維數；h為光滑長度。光滑長度隨時間空間而變化。

SPH 方程最常用的光滑內核是三次B-樣條函數，由下式定義：

式中：C為標準化常數，由空間維數決定。

邊界采用剛性墻與虛粒子2 種方法結合處理。方法1：用有限或無限大小定義1 個平面剛性壁，保證內部粒子無法穿過壁面。方法2：創(chuàng)建虛粒子來定義對稱面，虛粒子為靠近邊界兩倍初始光滑長度范圍內實粒子的鏡像，與對稱實粒子具有相同的質量、速度等物理量。不僅提高了計算效率，并且仿真結果的精度得到了保證。

接觸算法采用罰函數進行求解[5]，粒子的接觸力分解為法向接觸力fn和切向接觸力，計算公式為：

2 計算模型

2.1 LNG 船體模型

本文研究一艘航行于北極航道的LNG 船，破冰能力為ARC-7，自主破冰厚度可達1.7 m，主要參數如表1所示，仿真模型縮尺比例為1:100, 材料設為剛體。忽略船體與浮冰碰撞時的變形，簡化船體結構，保留船體外殼與內部主要骨架，用添加質量點的方式平衡船體的質量、重心位置和轉動慣量。

表1 船體主要參數Tab. 1 Main parameters of the hull

2.2 浮冰域

由于北極航道環(huán)境復雜多變，浮冰的大小與形狀也不盡相同，忽略不同大小與形狀帶來的結果偏差，只保證浮冰的密集度。將海面漂浮的浮冰統一設置成大小相同的長方體，為模擬真實環(huán)境，長方體的排布為隨機布置。浮冰為SPH 建立，粒子直徑為0.01 m，材料的本構模型選擇各向同性彈塑性斷裂模型，失效準則滿足VON-MISES 屈服準則，以最大塑性應變模式作為材料的破壞模式，材料的分離模式為恒定最小壓力模式[6]，具體參數見表2。

表2 冰材料模型參數Tab. 2 Ice material model parameters

由于船模寬度相較于拖曳水池寬度來說為小值，適當縮小浮冰域寬度對實驗結果影響不大，用剛性墻模擬冰水池試驗中浮式圍欄的影響，浮冰域設置為長9 m，寬2 m。

海水模型同樣使用SPH 建立，粒子半徑為0.02 m。水的狀態(tài)方程采用Gruneisen 狀態(tài)方程，壓力方程為[7]：

式中：C1為 水中聲音的傳播速度； α為 對系數 γ0的一階修正；S1-S3為 μs-μp曲線斜率無量綱系數；E為初始材料內能； μ為體積變化率。

圖2 為船在水中靜置10 s 的浮力隨時間變化的曲線圖。前2 s 船與水開始接觸，浮力數值發(fā)生較大變化，后8 s 趨于穩(wěn)定，浮力穩(wěn)定在1 275 N 附近，而船縮尺后滿載排水量為127 kg，可以認為水的浮力模擬結果較為準確。

圖2 浮力-時間曲線Fig. 2 Buoyancy-time curve

2.3 浮冰工況與結果對比

北極地區(qū)冬季最低溫可至-50℃以下，夏季平均氣溫在-10℃至10℃之間，某些地區(qū)最高溫度可達30℃以上，這就導致同一條航道中浮冰情況會有很大差別。綜合實際情況并縮尺選取浮冰厚度均為0.01 m，密集度為80%，60%，40%，船舶釋放六自由度，仿真時長9 s，前1 s 船舶靜置在水中達到正浮，后8 s 以0.448 m/s2沿船長方向直線勻速運動。水域四周與底部設置剛性墻防止粒子穿透，船體縱剖面設置虛粒子邊界。得出的阻力值與Du Brovin 在1950-1955 年通過4 次船模水池試驗得出的結果推導出針對船模試驗冰阻力預報的經驗公式[8]對比。計算船模試驗冰阻力Rice的經驗公式為：

式中：Rice為船舶所受冰阻力；p1，p2為經驗修正系數，與浮冰密集度和浮冰域寬度與船寬的比值有關；Fn為 傅汝德數；n為 功率系數。參數A與 φ的計算公式為：

式中：B為船型寬；L為船總長；r為浮冰密集度；hice為浮冰厚度； ρice為浮冰密度；f0為船體與浮冰的摩擦因數； αH為船舶前體棱形系數； α0為水線面入射角的1/2。

3 計算結果分析

3.1 浮冰應變云圖分析

各時間浮冰應變云圖如圖3～圖5 所示。

圖3 40%浮冰密集度的應變云圖Fig. 3 Strain cloud map of 40% ice floe density

圖4 60%浮冰密集度的應變云圖Fig. 4 Strain cloud map of 60% ice floe density

圖5 80%浮冰密集度的應變云圖Fig. 5 Strain cloud map of 80% ice floe density

從時間云圖可以看到，船首先與浮冰碰撞，接觸部位的浮冰應變達到0.001 失效應變后發(fā)生破碎，隨后被船首擠壓到兩邊位置，對周圍浮冰產生擠壓碰撞。碰撞主要發(fā)生在船體的首部與肩部，平行中體部位與尾部附近的浮冰應變幾乎沒有變化。但由于波浪和浮冰之間碰撞浮冰會向船首方向運動，密集度越低的浮冰之間運動的相對距離會遠，密集度高的浮冰會聚在一起，整體向前。當船全部進入浮冰區(qū)域后，船尾被碰撞的碎冰會重新聚合。

3.2 浮冰運動狀態(tài)分析

如圖6 所示，在碰撞時浮冰的運動狀態(tài)分為3 種情況。

圖6 浮冰不同運動狀態(tài)Fig. 6 Floating ice in different motion states

在仿真過程中可以明顯看到模型浮冰有推積、貼體、平移等現象。浮冰在與船首接觸時發(fā)生破碎且周圍浮冰的擠壓，導致部分碎冰無法從兩邊沿著船身向船尾方向移動，于是這些碎冰會堆積到船首的位置。船首的型線會導致在碰撞浮冰時，浮冰不僅會受到船長方向的碰撞力，也會受到來自船首向下的擠壓力。這些擠壓力將浮冰擠壓成大小不一的碎冰，碎冰受到水中浮力作用緊緊貼附在水下船殼表面。平移現象在密集度小的工況中更為明顯。因為浮冰在船頭的碰撞與船身的摩擦之下向前運動，由于前方水域較為開闊，無浮冰阻攔，所以僅收到水阻力的情況下向前漂浮一段距離后靜止。因為船速相對較低，所以整體工況中反轉現象并不明顯。這些浮冰的不同運動情況是使冰阻力增加的一個重要因素。

3.3 不同密集度下冰力隨時間變化曲線

船體分別在40%，60%，80%密集度浮冰區(qū)域航行時受到的冰阻力隨時間變化曲線如圖7～圖9 所示。

圖7 40%浮冰密集度冰力-時間曲線Fig. 7 40% Ice density ice force time curve

圖8 60%浮冰密集度冰力-時間曲線Fig. 8 60% Ice density ice force time curve

圖9 80%浮冰密集度冰力-時間曲線Fig. 9 80% Ice density ice force time curve

實線為冰阻力-時間曲線，The mean 劃線為4-6 s 船首進入浮冰區(qū)后數據平穩(wěn)段阻力的平均值。Experience 劃線為DuBrovin 經驗公式值。各項詳細數值如表3 所示。

表3 各項冰力數值表Tab. 3 Table of various ice force values

各浮冰密集度冰力-時間曲線具有明顯的動態(tài)非線性特征，趨勢分為3 部分，第1 部分0-2 s 時船開始向前運動，與浮冰未接觸，冰力為0。在2-4 s 中船首部與浮冰接觸，冰力開始增大。4-10 s 中船身與浮冰接觸，直到整個船身完全進入浮冰區(qū)域。在這個階段冰力大浮動波動，是因為船與浮冰接觸時冰力增大。浮冰收到碰撞力向力的相反方向移動冰力卸載，浮冰又受到水阻力與其他浮冰的影響速度減小與船體發(fā)生二次碰撞冰力增大，直至破碎或沿著船身移動。但浮動范圍變化很小，說明船身對冰力大小的影響較小。

冰力極值方面密集度為40%時峰值點較多且出現時間較晚，說明船舶航行中波浪對于浮冰的運動起重要作用。

由圖10 可知，隨著浮冰密集度的增加，浮冰阻力均值與經驗公式值均程增加趨勢，近似呈線性關系。仿真數據比經驗公式在60%，80% 的工況中數值偏大，是因為浮冰隨著密集度增加會出現破碎、堆積、貼體等不確定因素導致的阻力增加。

圖10 浮冰阻力均值和經驗公式值與浮冰密集度的散點圖Fig. 10 Scatter plot of mean and empirical formula values of floating ice resistance and floating ice density

4 結 語

本文基于SPH-FEM 耦合算法建立模型船在水中與不同密集度浮冰域碰撞的仿真場景并進行模擬，分析船體航行冰阻力變化情況，主要結論如下：

1）基于FEM-SPH 耦合算法對LNG 船在水中與浮冰碰撞過程有很好的模擬效果。

2）在LNG 船航行過程中，浮冰與船體接觸被破壞區(qū)域主要在船首與肩部，船身處浮冰應變變化不大。浮冰運動狀態(tài)主要有推積、貼體和平移，密集度小的工況中推積與貼體較少，浮冰平移較多。

3）冰力-時間曲線可反映船體與浮冰碰撞全過程。浮冰密集度小的工況中船體所受平均冰力較小，與DuBrovin 經驗公式對比，整體變化趨勢一致，在高密集度工況下仿真結果偏大。