波浪中接驳小艇船体和浮冰块的耦合动力学特性数值研究

2022-08-20刘永涛封培元

船舶 2022年4期

刘永涛封培元

（1. 江苏科技大学船舶与海洋工程学院镇江 212003；2.中国船舶及海洋工程设计研究院上海 200011）

0 引言

伴随着全球气候逐渐变暖，极地海区夏季融冰速度加快、覆冰范围减少，使当地具备了通航、观光旅行的条件。航行于北冰洋的船舶，将面临大面积漂浮海冰引起的航行阻力问题、船体局部结构冰载荷问题等。对航行于极地海区的邮轮，常配备接驳小艇，用于满足游客登临冰川的需求。在接驳小艇航行于浮冰区时，受波浪作用船体极易碰触浮冰块、爬升并形成骑冰状态，但不致于使浮冰发生破坏。此时船体的运动响应、船体和浮冰块的作用力等问题，少见公开发表的研究成果。

上述问题涉及到船体、波浪和浮冰块三者的相互耦合作用，既有船体和波浪、波浪和浮冰块的流固耦合作用，还包括船体和浮冰块2种固体间的耦合作用。针对该问题，本文应用光滑粒子水动力-离散单元法（基于DualSPHysics代码开发的SPH-DEM法），应用SPH (smoothed particle hydrodynamics) 方法模拟波浪运动，基于DEM(discrete element method)建立并求解船体、浮冰块相互作用的动力学方程，从而建立波浪、船体和浮冰块三者耦合的动力学方程，最后针对接驳小艇计算分析了三者耦合的动力学特性。

1 计算原理与方法

求解波浪、船体和浮冰块三者耦合的动力学问题，需分别对流体、船体和浮冰块2种固体建立数值计算模型。

1.1 流体运动数值计算模型

对于不可压缩流动，控制方程包括连续方程和动量方程，分别为：

式中：为流动相关的速度，m/s；为流动相关的密度；为流动相关的压强，Pa；为流动相关的质量力，m/s；为流动相关的黏度，m/s。

计算流场进行粒子离散化后，流体粒子就拥有位置、速度、密度和压强等物理量，在其影响半径内和临近粒子相互作用并交换物理量。这种粒子间相互作用、交换物理量的规律可用核函数来描述，具体为：

式中：=1，2，3，...，，并且为粒子影响半径范围Ω内的粒子数量；且m为粒子的质量，kg；ρ为粒子的密度，kg/m；r为粒子的位置；为影响半径，m。

本文采用Wendland核函数，其表达式为：

式中：r=|r-r|为粒子和粒子间距离，m；常数α在二维和三维取值分别为7/4π和7/4π，变量为q=r/h。

基于SPH方法的核函数近似，控制方程为：

基于若可压缩假定，流场内流体粒子处的压强可用状态方程表示为：

最终，每一瞬时流体粒子的位置为：

1.2 固体运动数值计算模型

固体运动求解，需要建立其运动方程、确定其所受的流体作用力和固体接触力，分别介绍如下。

对于某一固体，根据牛顿第二运动定律建立运动方程如下：

式中：M为固体的质量；为固体的重心位置；θ为固体的转动角度；F

（）、M （）为固体的重力和力矩；F （）、M （）为流体对固体的作用力和力矩；（t）、（t）为临近固体对固体的接触力和力矩。

固体的受力，在离散为大量弹性固体颗粒的基础上，则其所受的流体作用力和力矩、固体接触力和力矩可分别由每一固体颗粒所受的流体粒子作用力和力矩、临近固体颗粒接触力和力矩累加得到。其中，前者由SPH方法通过边界流体粒子作用力求和得到；后者用基于软球模型的DEM方法，考虑固体材料属性、颗粒间相对速度和重合量确定。考虑固体边界上的某一离散固体颗粒，由离散单元法，其受力和运动方程为：

1.3 流体-固体耦合的SPH-DEM计算模型

考虑流体和固体的耦合作用，根据流体在固体表面满足不可穿透物面条件，固体所受的流体作用力F 和力矩M ，为边界上固体颗粒对应受力F 和力矩M 分别求和得到：

式中：为固体边界的某一离散固体颗粒；为与其临近的流体粒子；R为固体颗粒与所在固体重心的距离矢量。

固体之间的接触力和力矩，考虑发生接触的固体和固体，对应边界接触颗粒分别为固体颗粒和，假定单个固体内颗粒无相对运动，通过累加计算得到发生接触的固体颗粒受力和力矩：

固体颗粒之间的法向接触力：

式中：k为颗粒法向弹性系数；η为颗粒法向阻尼系数；δ为颗粒间重合量；u为颗粒间法向相对速度，m/s。

固体颗粒之间的切向接触力f基于库伦摩擦定律确定：

式中：k为切向弹性系数；δ为切向重合量；η为切向阻尼系数；u为切向相对速度，m/s。

2 数值算例

计算船体选为接驳小艇，其主要参数为：长度3.1 m、吃水0.3 m、水线长2.9 m、重量为7 590 N。考虑2种重心位置，纵向分别为距舯-0.45 m和 0.05 m，垂向距船底0.3 m。浮冰块主体为梯形，中部左右侧有一水平台阶，整体左右对称，主要参数为：底部长度为5.0 m、顶部长度为1.5 m、高度为1.25 m、密度为930 kg/m、重量为41 700 N、重心距梯形底部0.51 m，如图2所示。

图2 船体和浮冰块尺度图

建立数值波浪水池，其长度为35 m，水深 4.0 m。水池右侧5 m长度范围为消波区，左侧布置造波机，产生周期为3.0 s、波高为0.15 m的规则波。计算粒子直径为0.02 m，共有22.56万个粒子。针对二维问题，考虑0°、180° 2个浪向，以及舯前、舯后2个船体重心位置，共有4种计算工况。其中，0°、180°浪向对应船体位于浮冰块上方的状态，此时船体和浮冰块在波浪中的纵摇、垂荡运动位移及二者相互作用力较大，也是实际船体与浮冰块常见的接触状态，具体如表1所示。

表1 计算工况表

为研究船体与浮冰块从初始接触到达到稳定相互作用的过程，假定初始时刻浮冰块顶部位于静水面，船体正浮且底部恰好接触中部台阶。此后，浮冰块将上浮并与船体接触。采用大地坐标系描述船体、浮冰块2个物体的运动和受力，按照右手螺旋法则建立坐标轴，沿波浪水池长度方向水平向右为轴正向，垂直向上为轴正向，倾角以逆时针为正。船体和浮冰块之间的接触力，与材料属性和局部变形量有关，船体材料为钢材，浮冰块材料为冰，假定船体、浮冰块受力均不破损，其局部变形量由2种颗粒间的瞬时重合量确定。其材料性能参数如表2所示。

表2 材料性能参数表

上述4种工况的数值模拟结果，分别给出典型时刻船体、浮冰块和波浪的运动状态图。船体运动，及所受流体力、重力、浮冰块对船体作用力和上述外力合力；浮冰块运动，及所受流体力、重力、船体对浮冰块作用力和上述外力合力，对应时历图及说明如图3～19所示。

图3 典型时刻船体、浮冰块、波浪运动状态图

Case1结果（0°浪向、纵向重心距舯0.05 m）

图4 船体和浮冰块运动时历图

图5 船体受力时历图

图6 浮冰块受力时历图

该工况船体重心距舯0.05 m，使船首依靠自身重量与上浮浮冰块在右侧梯形上斜面挤压，造成浮冰块右倾2.2°。并且，船体和浮冰块在波浪中的周期性纵摇、纵荡、垂荡运动幅值和相位较为一致。对比两者纵荡相对运动，可知船体相对浮冰块向右的纵荡运动量大，故受到水平向左的浮冰块作用力。在右侧梯形上斜面，由于浮冰块姿态右倾，故对船体施加了垂直向下的作用力，如图3～6所示。

Case2结果（0°浪向、纵向重心距舯-0.45 m）

图7 典型时刻船、冰、浪运动状态图

图9 船体受力时历图

图10 浮冰块受力时历图

该工况船体重心在舯后0.45 m，形成12.5°的船体右侧尾倾，并绕此角做波浪中的周期性纵摇运动，且船体纵摇峰谷值和浮冰块相应值反向出现，即发生船体纵摇角到达峰值、浮冰块纵摇角到达谷值现象，使浮冰块受到垂直向上的作用力。由纵荡时历曲线图可知，船体相对浮冰块向右侧运动，故对其施加了水平向右的拖曳力，在23 s后两者基本脱离接触，水平和垂直相互作用力接近于零，如图8～10所示。

图8 船体和浮冰块运动时历图

Case3结果（180°浪向、纵向重心距舯0.05 m）

图12 船体和浮冰块运动时历图

图13 船体受力时历图

图14 浮冰块受力时历图

该工况船体重心在舯前0.05 m，运动开始以后浮冰块上浮，造成其在左侧梯形上斜面与船首紧密接触、部分船体重量施加于浮冰块左侧，形成2.2°的左倾角。两者在波浪中的运动幅值和相位较为一致周期性纵摇运动。对比而言，浮冰块的纵荡运动量大于船体相应运动量，故受到水平向左的船体作用力，而在左侧的梯形上斜面，船体对浮冰块施加了垂直向上的作用力。如图11～14所示。

图1 SPH流体粒子和DEM固体颗粒临近关系图

图11 典型时刻船、冰、浪运动状态图

Case4结果（180°浪向、纵向重心距舯-0.45 m）

图15 典型时刻船、冰、浪运动状态图

图17 船体受力时历图

图18 浮冰块受力时历图

该工况船体重心在舯后0.45 m，模拟开始以后船体即产生较大尾倾，在波浪中绕12.5°尾倾角做周期性纵摇运动。相对于浮冰块，船体纵荡运动量在12 s后较快增大，使船首逐渐爬升至浮冰块左上侧梯形斜边顶部，随其纵摇运动对浮冰块施加方向为右下方的作用力，如图16～19所示。

图16 船体和浮冰块运动时历图

对比4种工况船体所受浮冰块作用力，可见垂直方向作用力极值大于水平方向作用力极值， 180°浪向工况Case3和Case4浮冰块作用力极值大于0°浪向工况Case1和Case2作用力极值。180°浪向船体位于浮冰块上游，船体和浮冰块相对运动量大；0°浪向船体位于浮冰块下游受其遮蔽效应影响，2者相对运动量较小，因此所受浮冰块作用力极值较小。两种重心工况Case1、Case3和Case2、Case4结果对比可见，Case1、Case3船体重心在舯前，船首处其重量施加于浮冰块形成两者的稳定接触，两者运动周期、幅值较为一致，整个运动过程浮冰块作用力存在稳定的时均值，其水平分力分别占自重的2.1%和4.7%，垂直分力分别占自重的44.6%和88.9%；而Case2、Case4船体重心在舯后形成稳定的尾倾，浮冰块、船体两者在23 s前期存在接触、23 s后期相对独立，浮冰块作用力在前期达到极值而在后期快速下降至接近于零值，见图19所示。