基于SPH-FEM 算法的浮冰-水-船耦合作用數(shù)值模擬方法

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

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

0 引 言

我國(guó)處于北半球，船運(yùn)走北極航線必須經(jīng)過(guò)邊緣冰帶（MIZ），這個(gè)區(qū)域處于極地和溫帶氣候系統(tǒng)的交匯處。該區(qū)域浮冰的大小從數(shù)米到數(shù)百米不等[1]。這些充滿(mǎn)隨機(jī)性的浮冰會(huì)對(duì)船舶的總體結(jié)構(gòu)與人員生命安全造成嚴(yán)重威脅。所以，船舶在浮冰區(qū)域的冰阻力預(yù)測(cè)成為船舶安全保障的關(guān)鍵。

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

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

1 數(shù)值模型

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

SPH 方程構(gòu)造首先采用核函數(shù)逼近的方法，將描述場(chǎng)的函數(shù)轉(zhuǎn)化成積分表達(dá)式，然后粒子近似，實(shí)現(xiàn)了對(duì)核函數(shù)近似積分表達(dá)式的粒子離散。

粒子近似式方程表達(dá)為：

式中，W為核函數(shù)。

核函數(shù)W由函 數(shù) θ定義：

式中：d為空間維數(shù)；h為光滑長(zhǎng)度。光滑長(zhǎng)度隨時(shí)間空間而變化。

SPH 方程最常用的光滑內(nèi)核是三次B-樣條函數(shù)，由下式定義：

式中：C為標(biāo)準(zhǔn)化常數(shù)，由空間維數(shù)決定。

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

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

2 計(jì)算模型

2.1 LNG 船體模型

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

表1 船體主要參數(shù)Tab. 1 Main parameters of the hull

2.2 浮冰域

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

表2 冰材料模型參數(shù)Tab. 2 Ice material model parameters

由于船模寬度相較于拖曳水池寬度來(lái)說(shuō)為小值，適當(dāng)縮小浮冰域?qū)挾葘?duì)實(shí)驗(yàn)結(jié)果影響不大，用剛性墻模擬冰水池試驗(yàn)中浮式圍欄的影響，浮冰域設(shè)置為長(zhǎng)9 m，寬2 m。

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

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

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

圖2 浮力-時(shí)間曲線Fig. 2 Buoyancy-time curve

2.3 浮冰工況與結(jié)果對(duì)比

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

式中：Rice為船舶所受冰阻力；p1，p2為經(jīng)驗(yàn)修正系數(shù)，與浮冰密集度和浮冰域?qū)挾扰c船寬的比值有關(guān)；Fn為 傅汝德數(shù)；n為 功率系數(shù)。參數(shù)A與 φ的計(jì)算公式為：

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

3 計(jì)算結(jié)果分析

3.1 浮冰應(yīng)變?cè)茍D分析

各時(shí)間浮冰應(yīng)變?cè)茍D如圖3～圖5 所示。

圖3 40%浮冰密集度的應(yīng)變?cè)茍DFig. 3 Strain cloud map of 40% ice floe density

圖4 60%浮冰密集度的應(yīng)變?cè)茍DFig. 4 Strain cloud map of 60% ice floe density

圖5 80%浮冰密集度的應(yīng)變?cè)茍DFig. 5 Strain cloud map of 80% ice floe density

從時(shí)間云圖可以看到，船首先與浮冰碰撞，接觸部位的浮冰應(yīng)變達(dá)到0.001 失效應(yīng)變后發(fā)生破碎，隨后被船首擠壓到兩邊位置，對(duì)周?chē)”a(chǎn)生擠壓碰撞。碰撞主要發(fā)生在船體的首部與肩部，平行中體部位與尾部附近的浮冰應(yīng)變幾乎沒(méi)有變化。但由于波浪和浮冰之間碰撞浮冰會(huì)向船首方向運(yùn)動(dòng)，密集度越低的浮冰之間運(yùn)動(dòng)的相對(duì)距離會(huì)遠(yuǎn)，密集度高的浮冰會(huì)聚在一起，整體向前。當(dāng)船全部進(jìn)入浮冰區(qū)域后，船尾被碰撞的碎冰會(huì)重新聚合。

3.2 浮冰運(yùn)動(dòng)狀態(tài)分析

如圖6 所示，在碰撞時(shí)浮冰的運(yùn)動(dòng)狀態(tài)分為3 種情況。

圖6 浮冰不同運(yùn)動(dòng)狀態(tài)Fig. 6 Floating ice in different motion states

在仿真過(guò)程中可以明顯看到模型浮冰有推積、貼體、平移等現(xiàn)象。浮冰在與船首接觸時(shí)發(fā)生破碎且周?chē)”臄D壓，導(dǎo)致部分碎冰無(wú)法從兩邊沿著船身向船尾方向移動(dòng)，于是這些碎冰會(huì)堆積到船首的位置。船首的型線會(huì)導(dǎo)致在碰撞浮冰時(shí)，浮冰不僅會(huì)受到船長(zhǎng)方向的碰撞力，也會(huì)受到來(lái)自船首向下的擠壓力。這些擠壓力將浮冰擠壓成大小不一的碎冰，碎冰受到水中浮力作用緊緊貼附在水下船殼表面。平移現(xiàn)象在密集度小的工況中更為明顯。因?yàn)楦”诖^的碰撞與船身的摩擦之下向前運(yùn)動(dòng)，由于前方水域較為開(kāi)闊，無(wú)浮冰阻攔，所以?xún)H收到水阻力的情況下向前漂浮一段距離后靜止。因?yàn)榇傧鄬?duì)較低，所以整體工況中反轉(zhuǎn)現(xiàn)象并不明顯。這些浮冰的不同運(yùn)動(dòng)情況是使冰阻力增加的一個(gè)重要因素。

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

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

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

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

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

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

表3 各項(xiàng)冰力數(shù)值表Tab. 3 Table of various ice force values

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

冰力極值方面密集度為40%時(shí)峰值點(diǎn)較多且出現(xiàn)時(shí)間較晚，說(shuō)明船舶航行中波浪對(duì)于浮冰的運(yùn)動(dòng)起重要作用。

由圖10 可知，隨著浮冰密集度的增加，浮冰阻力均值與經(jīng)驗(yàn)公式值均程增加趨勢(shì)，近似呈線性關(guān)系。仿真數(shù)據(jù)比經(jīng)驗(yàn)公式在60%，80% 的工況中數(shù)值偏大，是因?yàn)楦”S著密集度增加會(huì)出現(xiàn)破碎、堆積、貼體等不確定因素導(dǎo)致的阻力增加。

圖10 浮冰阻力均值和經(jīng)驗(yàn)公式值與浮冰密集度的散點(diǎn)圖Fig. 10 Scatter plot of mean and empirical formula values of floating ice resistance and floating ice density

4 結(jié) 語(yǔ)

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

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

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

3）冰力-時(shí)間曲線可反映船體與浮冰碰撞全過(guò)程。浮冰密集度小的工況中船體所受平均冰力較小，與DuBrovin 經(jīng)驗(yàn)公式對(duì)比，整體變化趨勢(shì)一致，在高密集度工況下仿真結(jié)果偏大。