劉永濤 封培元

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

0 引 言

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

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

1 計算原理與方法

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

1.1 流體運動數(shù)值計算模型

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

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

計算流場進(jìn)行粒子離散化后，流體粒子就擁有位置、速度、密度和壓強等物理量，在其影響半徑內(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ù)近似，控制方程為：

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

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

1.2 固體運動數(shù)值計算模型

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

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

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

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

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

1.3 流體-固體耦合的SPH-DEM計算模型

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

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

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

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

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

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

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

2 數(shù)值算例

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

圖2 船體和浮冰塊尺度圖

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

表1 計算工況表

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

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

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

圖3 典型時刻船體、浮冰塊、波浪運動狀態(tài)圖

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

圖4 船體和浮冰塊運動時歷圖

圖5 船體受力時歷圖

圖6 浮冰塊受力時歷圖

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

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

圖7 典型時刻船、冰、浪運動狀態(tài)圖

圖9 船體受力時歷圖

圖10 浮冰塊受力時歷圖

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

圖8 船體和浮冰塊運動時歷圖

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

圖12 船體和浮冰塊運動時歷圖

圖13 船體受力時歷圖

圖14 浮冰塊受力時歷圖

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

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

圖11 典型時刻船、冰、浪運動狀態(tài)圖

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

圖15 典型時刻船、冰、浪運動狀態(tài)圖

圖17 船體受力時歷圖

圖18 浮冰塊受力時歷圖

該工況船體重心在舯后0.45 m，模擬開始以后船體即產(chǎn)生較大尾傾，在波浪中繞12.5°尾傾角做周期性縱搖運動。相對于浮冰塊，船體縱蕩運動量在12 s后較快增大，使船首逐漸爬升至浮冰塊左上側(cè)梯形斜邊頂部，隨其縱搖運動對浮冰塊施加方向為右下方的作用力，如圖16～19所示。

圖16 船體和浮冰塊運動時歷圖

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

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

3 結(jié) 論

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

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

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

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