MPS-DEM耦合方法的多段系泊浮体波浪响应分析

孙一颉,席光,孙中国

(西安交通大学能源与动力工程学院,710049,西安)

海洋蕴含丰富的能量和资源,海洋资源和能源的开发离不开各种各样的海上设施,其中海上浮体平台作为一种常见的海上设施,在海洋工程中发挥着重要角色。为了保证大型浮体平台运动的稳定性和平台上的人员安全,平台结构常采用分段多模块式连接。海浪中多模块分段式耦合运动与单一整体自由运动存在明显差异,了解分段式平台的运动以及多模块之间的动态响应关系对合理设计大型浮体平台的结构十分重要。

海上浮体平台与海浪的相互作用是一个复杂的强非线性问题,包含浮体结构受到的非线性水动力载荷和非线性变形。浮体平台在海浪的冲击下会发生变形,平台的反复形变不仅可能会导致结构疲劳失效,还会引起海水与平台相互作用力的变化,这类复杂问题通常采用无网格数值方法对浮体平台的运动响应与结构受力进行预测。无网格法将计算区域离散为一系列的粒子,通过求解粒子的速度、位置、压力及时间层推进来获得流场的动态信息,与传统网格方法相比,无网格法避免了网格重构和畸变等问题,在捕捉海浪破碎、甲板上浪等流体大变形现象和模拟固体变形断裂等问题时有较大优势。

近年来,许多学者将无网格法或无网格法与其他方法的耦合算法应用于海上浮体结构物在海浪作用下的运动响应。Shibata等建立了三维刚性船体模型,并利用耦合移动粒子半隐式法(MPS)方法数值模拟了船舶在高海浪下的运动响应[1];Sun等采用可压缩SPH方法模拟了刚体和波浪自由面的相互作用[2];Zhang等利用MPS-FEM方法,模拟了柔性浮体和波浪的相互作用[3];Sun等将MPS方法与模态叠加法耦合,计算了溃坝冲击下浮体的变形与运动[4]。多数研究采用无网格法与传统网格法相耦合的方法来计算海浪与浮体间的相互作用,使用传统网格法离散的浮体结构区域仍存在网格,当浮体出现大变形或断裂时网格易出现畸变,影响计算结果精度。

因此,本文耦合了MPS、离散单元法(DEM)两种无网格方法,完全摆脱了数值计算对网格的依赖,充分发挥无网格法解决海浪与浮体非线性相互作用的优势,并利用MPS-DEM耦合算法,采用MPS方法模拟海浪运动,利用DEM方法离散分段浮体平台。首先对二维活塞式造波模型进行了修正,获得更加规则的波浪便于准确分析浮体响应;其次改进了消波模型并与解析解对比验证了二维数值水槽造波的准确性;最后数值模拟了分段浮体平台在规则波冲击下的运动响应,分析了子平台间的动态关系,讨论了子平台纵荡、垂荡和纵摇这3类运动响应。

1 数值方法

1.1 MPS方法简介

MPS方法用于模拟液体的非定常大变形运动,液体粒子在拉格朗日框架下的控制方程为

(1)

(2)

式中:ρ为流体密度;μ为动力黏性系数;u为流体速度;p为压力;f为体积力。

粒子之间的相互作用强弱由以下核函数衡量

(3)

式中:w(r)为核函数值;r=|ri-rj|为任意两个粒子i、j之间的距离;re为粒子相互作用半径。梯度算子采用文献[5]改进的模型,即

〈φ〉

(4)

压力泊松方程计算公式改进[6]为

(5)

参数α取值范围为(0,1),为了有效减小压力计算误差[7],本文取值为0.008。

本文液体自由表面的处理方式是将粒子数密度满足n<βn0(β=0.97)条件的粒子判定为表面粒子,并将液体表面粒子压力取值为0。

1.2 DEM方法简介

DEM方法用于计算固体的运动、变形及相互作用,其中关于固体运动的流固耦合模拟,固体粒子运动控制方程[8]为

(6)

(7)

式中:m、I、u、ω分别为每个固体粒子的质量、转动惯量、速度和角速度;Fc、Tc分别为固体粒子间相互接触产生的力和力矩;Fo、To分别为固体粒子受到的外力和外力矩;Fls、Tls为液体粒子对固体粒子的力和力矩。流体对固体的作用是通过压力实现的,压力作用于固体表面,方向沿固体表面法线方向指向固体内部,相关计算公式为

Fls=-∬SpdS=-∭V

(8)

Tls=-∬SrralpdS=-∭VrralpdV=

(9)

式中:rral为粒子对于固体质心的相对位置;l0为粒子直径。

关于固体的变形模拟,即固体粒子间的相互作用与变形关系,离散单元法接触模型如图1所示。粒子间接触力由刚性项和阻尼项构成,其法向分量和切向分量计算公式[9]为

Fcn=-knδn-dnun

(10)

(11)

式中:k、d、δ、u分别为刚度系数、阻尼系数、接触位移和粒子速度;f为粒子间摩擦系数;下标n表示法向;t表示切向。

图1 DEM接触模型

(a)粒子-粒子连接

(b)平行键连接图2 DEM粒子连接模型

2 造波数值模型及改进

2.1 造波修正与精确造波

造波池造波有两种方式:①直接设定造波板正弦运动[12];②由目标波形反演造波板运动规律[13]。第1种方式简单易行,直接给定造波板的运动幅值与周期即可,但此方法造出波浪的实际参数与造波板运动规律并不完全一致,不适用于对波形参数要求较高的情况。斜坡式消波的活塞造波池几何模型如图3所示,粒子初始间距为0.025 m,粒子总数为16 095,水池左侧为活塞式造波板,右侧布置斜坡式消波段,防止反射波对波浪运动的影响。采用第1种方式进行造波,给定造波板运动

(12)

式中:l为造波板振幅;T为造波板运动周期。

图3 斜坡式消波的活塞造波池几何模型

直接造波法的自由面与解析解比较和反演方式造波的自由面与解析解比较如图4、图5所示。与振幅为0.1 m、周期为1.4 s的正弦波相比,虽然波浪周期与造波板运动周期相同,但是造波板振幅为0.1 m时,造出的实际波高更大且存在波动。图5中采用反演方式造出的波形振幅和周期与解析解基本吻合。

图4 直接造波法的自由面与解析解比较

图5 反演方式造波的自由面与解析解比较

采用第2种方式造波时,活塞式造波板的水平方向速度运动规律为

(13)

(14)

ηp=Hcos(ωt)

(15)

式中:ω为波浪频率;W为水力传递函数;ηp为目标波形;k为波数;d为水深;H为波高。

两种造波方式下造波板的运动速度和位移比较如图6所示。由图6可知,在相同波形输出条件要求下,反演式造波板的运动幅值和速度均小于直接造波。直接式造波方式虽然简单,但是造出的波振幅被放大,而反演式造波通过目标波形条件来控制造波板运动,实现波浪参数实时反馈修正,适合应用于波浪条件已知、对波形参数要求较高的情况。

(a)造波板运动速度比较

(b)造波板横向位移比较图6 两种造波方式造波板运动速度和位移比较

2.2 消波段改进

斜坡式消波段通过将液体横向运动转为竖直方向运动来逐渐降低波浪传播速度,从而消除二次反射波,阻尼式消波的活塞造波池几何模型如图7所示。斜坡式消波效果较好,但在波长较长时,消波段所需长度较长,会大大增加计算量。本文采用阻尼式消波方式[13],在相同粒子尺寸时,粒子总数较斜坡式消波明显降低,相对减少了12.1%,有效提高了计算效率。

图7 阻尼式消波的活塞造波池几何模型

消波段内采用对粒子加速度的修正实现对波浪能的衰减,修正公式[13]为

(16)

阻尼式消波自由面高度与解析解比较如图8所示,可知消波段改进后液面高度的监测值和解析解吻合较好,说明消波形式的改变并未影响到池内波形传播,消波效果与斜坡式消波基本相同。

图8 阻尼式消波自由面高度与解析解比较

3 算例与结果分析

本文基于MPS-DEM耦合方法,建立了双向流固耦合模型,并通过数值计算与其他方法和实验做了对比[8],验证了方法的准确性。在此基础上,数值模拟了分段浮体平台在规则波冲击下的运动,分段浮体平台的首段和尾段分别系泊于造波池底部,不计系泊缆绳的质量,如图9所示。相邻两个子平台间设立连接部件,平台相对连接部件的运动有3个自由度,连接部件与子平台内设置连接粒子,连接粒子间作用包含弹性作用和阻尼作用,通过改变连接粒子间的弹性系数与阻尼系数,达到控制平台最大旋转角度与最大平移的目的,平台间约束作用示意图如图10所示。

图9 分段浮体平台几何模型

(a)平台间约束作用

(b)连接粒子间作用关系图10 平台间约束作用示意图

二维情况下每段子平台相对连接部件有纵荡、垂荡和纵摇3个自由度。本文针对含有2～5段子平台的4种分段浮体平台在波浪中的运动进行了数值模拟,并分析子平台间的运动响应。通过分析不同系泊方式对浮体运动的影响[14],本文中的4组算例选取了定位效果更好的首尾双缆绳系泊形式。分段浮体平台算例参数如表1所示。

表1 分段浮体平台算例参数

含有5节子平台的浮体在规则波冲击下运动的局部时序图如图11所示,深色部分为分段浮体平台,浅色部分为海水。由图11可知:t=2.0 s时规则波刚到达平台前端,首段子平台前仰带动后部平台;t=2.2 s时波浪冲击首段子平台,由于子平台前端系泊,缆绳牵拉首段子平台,平台仰角减小浸入海水,少量海水跃上平台;t为2.4～2.6 s时分段平台主体受到波浪冲击,分段结构呈弯曲状态,每个子平台下部与海水紧密接触,削弱了平台主体的整体转动;t为2.8～3.0 s时,波浪离开平台,平台恢复初始平展状态。

(a)t=2.0 s (b)t=2.2 s

(c)t=2.4 s (d)t=2.6 s

(e)t=2.8 s (f)t=3.0 s图11 5段浮体在规则波冲击下运动的局部时序图

(a)2段浮体

(b)3段浮体

(c)4段浮体

图12为含有2～5节子平台的4种浮体平台的纵荡响应,纵荡指平台沿最长延展方向的水平线性振荡运动,即平台沿x轴的振荡。由于每个子平台长度相同,因此相邻子平台间纵荡曲线起始点(t=0 s)间隔相同;子平台位置越靠近浮体尾部,纵荡振幅越大,具体数据如表2所示,尾部的子平台在波浪冲击和连接部件约束下呈正弦运动,正弦运动幅值越大,说明尾部子平台较前部子平台需要连接部件更大的约束力以保证连接稳定。在极端海浪冲击下,大型分段浮体平台需注意尾部子平台水平方向的振荡,防止子平台间的相互碰撞和子平台脱离分段浮体平台主体。

(d)5段浮体图12 4种浮体平台的纵荡响应

(a)2段浮体

(b)3段浮体

(c)4段浮体

(d)5段浮体图13 4种浮体平台的垂荡响应

4种浮体平台的垂荡响应如图13所示。由图13可知:子平台沿波浪传播方向的位置先后,使得垂荡曲线存在一定相位差;随着波浪沿平台尾部传递,能量逐段递减,尾段子平台垂荡幅度最小;t=2.5 s时第2段子平台垂荡响应远大于首段子平台,浮体平台在一个周期内的平均运动响应如表2所示,这是由于首段子平台前端系泊,当海浪浪高较大时系泊缆绳会对子平台产生位移约束,避免子平台运动幅值过大而导致分段浮体平台主体振荡不稳。

(a)2段浮体

(b)3段浮体

(c)4段浮体

4种浮体平台与水平线的夹角变化,反映了各个子平台的纵摇响应,4种浮体平台的纵摇响应如图14所示,纵摇是平台绕横轴的回转振荡运动,文中定义平台下俯时,与水平线的夹角为正。由图14可知,波浪刚冲击到首段子平台时,平台前端上仰,纵摇角度急剧下降由正转负,随后通过波峰,平台前端下俯,纵摇由负转正,整个过程纵摇曲线呈三角函数趋势。由于尾部子平台最后受到波浪冲击,波浪能量衰减,因而纵摇振幅最小,纵摇角度变化幅度随着子平台数量增加而逐段递减。

由表2可知:子平台越靠近尾部,平台的平均纵荡幅值越大,平均纵摇幅值越小;由于第1段子平台前端系泊于造波池底,运动范围受到系泊缆绳的限制,因而平均垂荡幅值的最大值出现于第2段子平台。随着子平台数量的增加,纵荡幅值逐渐降低,整体平台更加稳定。

表2 浮体平台在一个周期内的平均运动响应

4 结 论

本文利用无网格法计算波浪运动大变形的优势和离散元法模拟多体运动的特长,将无网格MPS法和离散元DEM法进行耦合,建立了流固(多体)耦合数值模拟方法;采用反演式造波形式,提升了造波水池的造波精度和稳定性;利用阻尼式消波的技术,缩减了计算区域,提升了计算效率;基于系泊浮体数值模型,设置并实现了首尾系泊的多段分段浮体平台在规则波冲击下的运动的数值模拟。

研究结果表明,分段浮体平台的子平台位置越靠近整体平台下游,垂荡和纵摇响应逐渐削弱,纵荡响应成为主导特征。大型分段浮体平台在设计时可通过增加系泊缆绳等措施来降低前部子平台的位移,进而减小前部子平台的垂荡和纵摇;其次,还应增加连接部件对尾部子平台的约束力,防止极端海况下尾部子平台间的相互碰撞以及尾部子平台与主体平台的脱离;同时,适当增加子平台数量有利于减小整体平台的纵荡幅值,保持稳定性。研究结果可为大型浮体平台的设计及优化提供参考和依据。