劉永濤 封培元

（1. 江蘇科技大學(xué) 船舶與海洋工程學(xué)院 鎮(zhèn)江 212003；2.中國(guó)船舶及海洋工程設(shè)計(jì)研究院 上海 200011）

0 引 言

伴隨著全球氣候逐漸變暖，極地海區(qū)夏季融冰速度加快、覆冰范圍減少，使當(dāng)?shù)鼐邆淞送ê?、觀光旅行的條件。航行于北冰洋的船舶，將面臨大面積漂浮海冰引起的航行阻力問題、船體局部結(jié)構(gòu)冰載荷問題等。對(duì)航行于極地海區(qū)的郵輪，常配備接駁小艇，用于滿足游客登臨冰川的需求。在接駁小艇航行于浮冰區(qū)時(shí)，受波浪作用船體極易碰觸浮冰塊、爬升并形成騎冰狀態(tài)，但不致于使浮冰發(fā)生破壞。此時(shí)船體的運(yùn)動(dòng)響應(yīng)、船體和浮冰塊的作用力等問題，少見公開發(fā)表的研究成果。

上述問題涉及到船體、波浪和浮冰塊三者的相互耦合作用，既有船體和波浪、波浪和浮冰塊的流固耦合作用，還包括船體和浮冰塊2種固體間的耦合作用。針對(duì)該問題，本文應(yīng)用光滑粒子水動(dòng)力-離散單元法（基于DualSPHysics代碼開發(fā)的SPH-DEM法），應(yīng)用SPH (smoothed particle hydrodynamics) 方法模擬波浪運(yùn)動(dòng)，基于DEM(discrete element method)建立并求解船體、浮冰塊相互作用的動(dòng)力學(xué)方程，從而建立波浪、船體和浮冰塊三者耦合的動(dòng)力學(xué)方程，最后針對(duì)接駁小艇計(jì)算分析了三者耦合的動(dòng)力學(xué)特性。

1 計(jì)算原理與方法

求解波浪、船體和浮冰塊三者耦合的動(dòng)力學(xué)問題，需分別對(duì)流體、船體和浮冰塊2種固體建立數(shù)值計(jì)算模型。

1.1 流體運(yùn)動(dòng)數(shù)值計(jì)算模型

對(duì)于不可壓縮流動(dòng)，控制方程包括連續(xù)方程和動(dòng)量方程，分別為：

式中：為流動(dòng)相關(guān)的速度，m/s；為流動(dòng)相關(guān)的密度；為流動(dòng)相關(guān)的壓強(qiáng)，Pa；為流動(dòng)相關(guān)的質(zhì)量力，m/s；為流動(dòng)相關(guān)的黏度，m/s。

計(jì)算流場(chǎng)進(jìn)行粒子離散化后，流體粒子就擁有位置、速度、密度和壓強(qiáng)等物理量，在其影響半徑內(nèi)和臨近粒子相互作用并交換物理量。這種粒子間相互作用、交換物理量的規(guī)律可用核函數(shù)來描述，具體為：

式中：=1，2，3，...，，并且為粒子影響半徑范圍Ω內(nèi)的粒子數(shù)量；且m為粒子的質(zhì)量，kg；ρ為粒子的密度，kg/m；r為粒子的位置；為影響半徑，m。

本文采用Wendland核函數(shù)，其表達(dá)式為：

式中：r=|r-r|為粒子和粒子間距離，m；常數(shù)α在二維和三維取值分別為7/4π和7/4π，變量為q=r/h。

基于SPH方法的核函數(shù)近似，控制方程為：

基于若可壓縮假定，流場(chǎng)內(nèi)流體粒子處的壓強(qiáng)可用狀態(tài)方程表示為：

最終，每一瞬時(shí)流體粒子的位置為：

1.2 固體運(yùn)動(dòng)數(shù)值計(jì)算模型

固體運(yùn)動(dòng)求解，需要建立其運(yùn)動(dòng)方程、確定其所受的流體作用力和固體接觸力，分別介紹如下。

對(duì)于某一固體，根據(jù)牛頓第二運(yùn)動(dòng)定律建立運(yùn)動(dòng)方程如下：

式中：M為固體的質(zhì)量；為固體的重心位置；θ為固體的轉(zhuǎn)動(dòng)角度；F

（）、M （）為固體的重力和力矩；F （）、M （）為流體對(duì)固體的作用力和力矩；（t）、（t）為臨近固體對(duì)固體的接觸力和力矩。

固體的受力，在離散為大量彈性固體顆粒的基礎(chǔ)上，則其所受的流體作用力和力矩、固體接觸力和力矩可分別由每一固體顆粒所受的流體粒子作用力和力矩、臨近固體顆粒接觸力和力矩累加得到。其中，前者由SPH方法通過邊界流體粒子作用力求和得到；后者用基于軟球模型的DEM方法，考慮固體材料屬性、顆粒間相對(duì)速度和重合量確定。考慮固體邊界上的某一離散固體顆粒，由離散單元法，其受力和運(yùn)動(dòng)方程為：

1.3 流體-固體耦合的SPH-DEM計(jì)算模型

考慮流體和固體的耦合作用，根據(jù)流體在固體表面滿足不可穿透物面條件，固體所受的流體作用力F 和力矩M ，為邊界上固體顆粒對(duì)應(yīng)受力F 和力矩M 分別求和得到：

式中：為固體邊界的某一離散固體顆粒；為與其臨近的流體粒子；R為固體顆粒與所在固體重心的距離矢量。

固體之間的接觸力和力矩，考慮發(fā)生接觸的固體和固體，對(duì)應(yīng)邊界接觸顆粒分別為固體顆粒和，假定單個(gè)固體內(nèi)顆粒無相對(duì)運(yùn)動(dòng)，通過累加計(jì)算得到發(fā)生接觸的固體顆粒受力和力矩：

固體顆粒之間的法向接觸力：

式中：k為顆粒法向彈性系數(shù)；η為顆粒法向阻尼系數(shù)；δ為顆粒間重合量；u為顆粒間法向相對(duì)速度，m/s。

固體顆粒之間的切向接觸力f基于庫(kù)倫摩擦定律確定：

式中：k為切向彈性系數(shù)；δ為切向重合量；η為切向阻尼系數(shù)；u為切向相對(duì)速度，m/s。

2 數(shù)值算例

計(jì)算船體選為接駁小艇，其主要參數(shù)為：長(zhǎng)度3.1 m、吃水0.3 m、水線長(zhǎng)2.9 m、重量為7 590 N?？紤]2種重心位置，縱向分別為距舯-0.45 m和 0.05 m，垂向距船底0.3 m。浮冰塊主體為梯形，中部左右側(cè)有一水平臺(tái)階，整體左右對(duì)稱，主要參數(shù)為：底部長(zhǎng)度為5.0 m、頂部長(zhǎng)度為1.5 m、高度為1.25 m、密度為930 kg/m、重量為41 700 N、重心距梯形底部0.51 m，如圖2所示。

圖2 船體和浮冰塊尺度圖

建立數(shù)值波浪水池，其長(zhǎng)度為35 m，水深 4.0 m。水池右側(cè)5 m長(zhǎng)度范圍為消波區(qū)，左側(cè)布置造波機(jī)，產(chǎn)生周期為3.0 s、波高為0.15 m的規(guī)則波。計(jì)算粒子直徑為0.02 m，共有22.56萬(wàn)個(gè)粒子。針對(duì)二維問題，考慮0°、180° 2個(gè)浪向，以及舯前、舯后2個(gè)船體重心位置，共有4種計(jì)算工況。其中，0°、180°浪向?qū)?yīng)船體位于浮冰塊上方的狀態(tài)，此時(shí)船體和浮冰塊在波浪中的縱搖、垂蕩運(yùn)動(dòng)位移及二者相互作用力較大，也是實(shí)際船體與浮冰塊常見的接觸狀態(tài)，具體如表1所示。

表1 計(jì)算工況表

為研究船體與浮冰塊從初始接觸到達(dá)到穩(wěn)定相互作用的過程，假定初始時(shí)刻浮冰塊頂部位于靜水面，船體正浮且底部恰好接觸中部臺(tái)階。此后，浮冰塊將上浮并與船體接觸。采用大地坐標(biāo)系描述船體、浮冰塊2個(gè)物體的運(yùn)動(dòng)和受力，按照右手螺旋法則建立坐標(biāo)軸，沿波浪水池長(zhǎng)度方向水平向右為軸正向，垂直向上為軸正向，傾角以逆時(shí)針為正。船體和浮冰塊之間的接觸力，與材料屬性和局部變形量有關(guān)，船體材料為鋼材，浮冰塊材料為冰，假定船體、浮冰塊受力均不破損，其局部變形量由2種顆粒間的瞬時(shí)重合量確定。其材料性能參數(shù)如表2所示。

表2 材料性能參數(shù)表

上述4種工況的數(shù)值模擬結(jié)果，分別給出典型時(shí)刻船體、浮冰塊和波浪的運(yùn)動(dòng)狀態(tài)圖。船體運(yùn)動(dòng)，及所受流體力、重力、浮冰塊對(duì)船體作用力和上述外力合力；浮冰塊運(yùn)動(dòng)，及所受流體力、重力、船體對(duì)浮冰塊作用力和上述外力合力，對(duì)應(yīng)時(shí)歷圖及說明如圖3～19所示。

圖3 典型時(shí)刻船體、浮冰塊、波浪運(yùn)動(dòng)狀態(tài)圖

Case1結(jié)果（0°浪向、縱向重心距舯0.05 m）

圖4 船體和浮冰塊運(yùn)動(dòng)時(shí)歷圖

圖5 船體受力時(shí)歷圖

圖6 浮冰塊受力時(shí)歷圖

該工況船體重心距舯0.05 m，使船首依靠自身重量與上浮浮冰塊在右側(cè)梯形上斜面擠壓，造成浮冰塊右傾2.2°。并且，船體和浮冰塊在波浪中的周期性縱搖、縱蕩、垂蕩運(yùn)動(dòng)幅值和相位較為一致。對(duì)比兩者縱蕩相對(duì)運(yùn)動(dòng)，可知船體相對(duì)浮冰塊向右的縱蕩運(yùn)動(dòng)量大，故受到水平向左的浮冰塊作用力。在右側(cè)梯形上斜面，由于浮冰塊姿態(tài)右傾，故對(duì)船體施加了垂直向下的作用力，如圖3～6所示。

Case2結(jié)果（0°浪向、縱向重心距舯-0.45 m）

圖7 典型時(shí)刻船、冰、浪運(yùn)動(dòng)狀態(tài)圖

圖9 船體受力時(shí)歷圖

圖10 浮冰塊受力時(shí)歷圖

該工況船體重心在舯后0.45 m，形成12.5°的船體右側(cè)尾傾，并繞此角做波浪中的周期性縱搖運(yùn)動(dòng)，且船體縱搖峰谷值和浮冰塊相應(yīng)值反向出現(xiàn)，即發(fā)生船體縱搖角到達(dá)峰值、浮冰塊縱搖角到達(dá)谷值現(xiàn)象，使浮冰塊受到垂直向上的作用力。由縱蕩時(shí)歷曲線圖可知，船體相對(duì)浮冰塊向右側(cè)運(yùn)動(dòng)，故對(duì)其施加了水平向右的拖曳力，在23 s后兩者基本脫離接觸，水平和垂直相互作用力接近于零，如圖8～10所示。

圖8 船體和浮冰塊運(yùn)動(dòng)時(shí)歷圖

Case3結(jié)果（180°浪向、縱向重心距舯0.05 m）

圖12 船體和浮冰塊運(yùn)動(dòng)時(shí)歷圖

圖13 船體受力時(shí)歷圖

圖14 浮冰塊受力時(shí)歷圖

該工況船體重心在舯前0.05 m，運(yùn)動(dòng)開始以后浮冰塊上浮，造成其在左側(cè)梯形上斜面與船首緊密接觸、部分船體重量施加于浮冰塊左側(cè)，形成2.2°的左傾角。兩者在波浪中的運(yùn)動(dòng)幅值和相位較為一致周期性縱搖運(yùn)動(dòng)。對(duì)比而言，浮冰塊的縱蕩運(yùn)動(dòng)量大于船體相應(yīng)運(yùn)動(dòng)量，故受到水平向左的船體作用力，而在左側(cè)的梯形上斜面，船體對(duì)浮冰塊施加了垂直向上的作用力。如圖11～14所示。

圖1 SPH流體粒子和DEM固體顆粒臨近關(guān)系圖

圖11 典型時(shí)刻船、冰、浪運(yùn)動(dòng)狀態(tài)圖

Case4結(jié)果（180°浪向、縱向重心距舯-0.45 m）

圖15 典型時(shí)刻船、冰、浪運(yùn)動(dòng)狀態(tài)圖

圖17 船體受力時(shí)歷圖

圖18 浮冰塊受力時(shí)歷圖

該工況船體重心在舯后0.45 m，模擬開始以后船體即產(chǎn)生較大尾傾，在波浪中繞12.5°尾傾角做周期性縱搖運(yùn)動(dòng)。相對(duì)于浮冰塊，船體縱蕩運(yùn)動(dòng)量在12 s后較快增大，使船首逐漸爬升至浮冰塊左上側(cè)梯形斜邊頂部，隨其縱搖運(yùn)動(dòng)對(duì)浮冰塊施加方向?yàn)橛蚁路降淖饔昧?，如圖16～19所示。

圖16 船體和浮冰塊運(yùn)動(dòng)時(shí)歷圖

對(duì)比4種工況船體所受浮冰塊作用力，可見垂直方向作用力極值大于水平方向作用力極值， 180°浪向工況Case3和Case4浮冰塊作用力極值大于0°浪向工況Case1和Case2作用力極值。180°浪向船體位于浮冰塊上游，船體和浮冰塊相對(duì)運(yùn)動(dòng)量大；0°浪向船體位于浮冰塊下游受其遮蔽效應(yīng)影響，2者相對(duì)運(yùn)動(dòng)量較小，因此所受浮冰塊作用力極值較小。兩種重心工況Case1、Case3和Case2、Case4結(jié)果對(duì)比可見，Case1、Case3船體重心在舯前，船首處其重量施加于浮冰塊形成兩者的穩(wěn)定接觸，兩者運(yùn)動(dòng)周期、幅值較為一致，整個(gè)運(yùn)動(dòng)過程浮冰塊作用力存在穩(wěn)定的時(shí)均值，其水平分力分別占自重的2.1%和4.7%，垂直分力分別占自重的44.6%和88.9%；而Case2、Case4船體重心在舯后形成穩(wěn)定的尾傾，浮冰塊、船體兩者在23 s前期存在接觸、23 s后期相對(duì)獨(dú)立，浮冰塊作用力在前期達(dá)到極值而在后期快速下降至接近于零值，見圖19所示。

圖19 4種工況船體所受浮冰塊作用力對(duì)比圖

3 結(jié) 論

本文應(yīng)用SPH-DEM方法，以接駁小艇為對(duì)象，針對(duì)船體、浮冰塊在波浪中的運(yùn)動(dòng)響應(yīng)、相互作用力問題，考慮浪向、重心位置因素的影響，開展了多工況數(shù)值計(jì)算，得到有關(guān)結(jié)論如下：

（1）考慮相同船體重心位置、不同浪向，船體位于浮冰塊上游的180°浪向條件下，船體首部受波浪力作用向上爬升至浮冰塊頂部；而在船體位于浮冰塊下游的0°浪向條件下，船體位于浮冰塊下游遮蔽區(qū)運(yùn)動(dòng)量減??；因此，船體和浮冰塊作用力極值在180°浪向結(jié)果大于0°浪向相應(yīng)結(jié)果。

（2）考慮相同浪向、不同船體重心位置，在船體、浮冰塊兩者運(yùn)動(dòng)穩(wěn)定后，船體重心位于舯前，船首處兩者緊密接觸，形成穩(wěn)定而明顯的時(shí)均作用力；在船體重心位于舯后，船體尾傾明顯，船體、浮冰塊處于相對(duì)自由狀態(tài)，作用力下降至接近于零值。

（3）初始自由漂浮條件下，船體和浮冰塊作用力的4種工況結(jié)果對(duì)比可知，作用力的垂直方向分量占主要部分，最大穩(wěn)定時(shí)均力約為船體重量的88.9%，是需要關(guān)注的主要載荷成分。