MPS-FEM 方法模拟弹性船体在规则波中的运动1)

黄聪祎 赵伟文 万德成

(上海交通大学 船舶海洋与建筑工程学院,船海计算水动力学研究中心,上海 200240)

引言

船舶在海上航行时,经常会遭遇波浪的作用.在波浪的作用下,船体可能会发生较大幅度的运动响应,给船上人员和生活带来诸多不便.另外,船体的大幅度运动也可能会造成船体表面与水面的剧烈砰击,作用在船体表面的砰击压力可能会破坏船体结构,危害船体和所载人员的安全.因此,准确地预报船体的运动响应,并估计船体受到的砰击载荷,对于提高船舶的安全系数,提升操作性能都具有重要的意义.

船体在波浪中的运动属于一典型的流固耦合问题.当船体运动幅度较小时,船体的弹性可以忽略,船体运动可以简单的视为刚体浮体的运动.但当波高较高,波浪砰击载荷较大时,船体可能会在波浪的作用下发生变形,此时船体的弹性则不可忽略.

近年来,国内外诸多学者都在研究浮体在波浪中的运动方面完成了许多工作.实验方法在处理此类问题时成本较高,实验难度也较大.在分析结构受力时,只能通过在结构表面附着压力测点或应变片的方式来实现.这样不能获得结构表面连续的压力分布,很难准确地捕捉到船体表面的最大压力,更难以分析流场中的压力分布情况.随着近些年高性能计算技术的发展以及软硬件性能的提高,数值模拟凭其高效、环保、适应性强的特点,被越来越多的研究人员选择.

在求解流固耦合问题时,由于流体和固体的物理差异巨大,因此,往往对流体域和结构域采用不同的方法分别计算,然后通过流固界面上信息的传递实现二者的耦合.目前已经有很多研究人员采用耦合的方法对流固耦合问题进行模拟.王加夏等[1]通过计算流体力学分析软件STAR-CCM+模拟流场,使用有限单元分析软件Abaqus 模拟结构,通过二者之间的耦合,模拟了弹性船体在波浪中的运动,结果与试验结果吻合较好.Oberhagemann等[2]使用黏流求解器模拟流场,使用有限元求解器ANSYS 模拟结构,模拟了LNG 船在波浪中的运动,并计算了船体受到的砰击载荷.Lakshmynarayanana等[3]使用RANS 方法,在StarCCM+平台上模拟流体,使用FEM方法使用ABAQUS 软件求解结构,通过二者的耦合模拟了弹性船体在规则波浪中运动响应.Kim等[4]将BEM 方法与FEM 方法耦合,模拟了大型船舶在迎浪时的变形情况.Malenica等[5]将FEM 方法与格林函数法结合,模拟了集装箱船在波浪中运动的流固耦合问题.

以上例子证明了在求解流固耦合问题时,耦合类的方法具有较好的表现.本文使用MPS-FEM 耦合方法,流体域的计算采用MPS 方法.MPS 方法作为一种无网格方法,追踪粒子运动,具有拉格朗日特性,其在处理自由表面变形较大的问题时,具有天然的优势.目前,无网格粒子法也被广泛地应用到对流固耦合问题的模拟中.Gao等[6]使用SPH 方法在二维数值水池中,模拟了规则波与结构的相互作用,分析了规则波与结构底面相互作用附近的压力分布和速度分布;Omidvar等[7]研究了二维浮体与波浪相互作用问题.也有一些研究人员利用粒子法模拟了结构与规则波之间的相互作用.Sueyoshi等[8-9]首先模拟了二维的波浪中的浮体运动问题,计算得到的浮体的运动曲线与试验数据相比吻合较好.随后又增大规则波的波幅,模拟船体剖面在规则波中的剧烈运动.饶成平等[10]使用MPS 耦合方法计算了孤立波对弹性结构的砰击载荷.Zhang等[11]使用MPSFEM 耦合方法,模拟了孤立波对弹性平板的砰击过程.Lind等[12]考虑砰击过程中空气的影响,采用两相流模型模拟了波浪对固定平板的砰击.Sun等[13]使用两相SPH 模型,模拟了畸形波对固定方柱的砰击作用.Zhang等[14]使用MPS-FEM 方法,模拟了弹性浮体和弹性系泊平台在孤立波中的受力及运动.Zhang等[15]用MPS 方法模拟了二维浮体在规则波中的运动,研究了规则波频率对浮体的横摇响应的影响.Xie等[16]采用MPS-DEM耦合方法,模拟了流体与弹性薄壁结构的相互作用.Zhang等[17]采用MPS-FEM 方法,模拟了三维溃坝流域弹性闸门的相互作用.Sun等[18]采用MPS-DEM耦合方法,模拟了弹性系泊平台在遭遇孤立波时受到的砰击压力及变形情况.Sun等[19-20]使用MPS 方法,模拟了二维和三维浮式平台在剧烈波浪作用下的变形和运动.

本文使用MPS-FEM 方法模拟弹性船体在规则波中的运动.在第1 节中,分别介绍了MPS 方法,FEM 方法的控制方程和求解流程,以及二者之间的耦合方法.在第2 节中,分别验证了MPSFEMSJTU 求解器的造波消波模块,以及模拟弹性体在规则波中的变形情况的准确性和稳定性.在第3 节中,首先模拟了刚性船体在规则波中的运动,研究波长对刚性船体运动响应的影响;接着模拟了弹性船体在规则波中的运动,并与刚性船体的结果对比,分析船舶的弹性对运动的影响.最后,第4 节总结了本文得到的结论.

1 数值方法

本文采用MPS-FEM 耦合方法,模拟弹性船体在规则波中的运动.在流体域使用MPS 方法,在结构域使用FEM 方法.在本节中,将分别介绍MPS 方法和FEM 方法的求解流程,以及二者的耦合方法.

1.1 控制流场的MPS 方法

本文采用一种无网格方法,移动粒子半隐式方法(moving particle semi-implicit method,MPS)来模拟流体的运动.传统的网格类数值模拟方法将计算域离散为一系列连续的网格,这种方法处理具有自由面大变形特性的问题时,很容易出现网格畸变的现象,导致计算的失效.而具有拉格朗日特性的无网格方法则不会出现这样的问题.MPS 方法通过一系列均匀分布的粒子来离散计算域.每个粒子都携带着自身的物理信息,包括粒子类型,位置,压力,速度等.每个粒子之间不存在固定的拓扑关系,而是通过一个权函数—核函数,来传递压力或位移等物理信息.通过追踪流场内每个粒子的运动来实现对流动的模拟.基于这一特点,无网格方法可以非常容易地追踪自由表面,这是网格类方法无可取代的优势.因此,无网格方法被广泛地应用于处理自由表面变形较剧烈的问题.

在MPS 方法中,控制方程为模拟不可压缩黏性流体的连续性方程和N-S(Navier-Stokes)方程.拉格朗日形式的控制方程可写为

式中,ρ为流体密度,P为压力,V为速度向量,g是重力,v是流体的运动黏性系数.

接下来介绍MPS 方法对计算域的离散方法.MPS 方法在离散计算域时不需要网格,而是将计算域离散为一系列无固定拓扑关系的粒子,粒子的运动由其邻居粒子的相应信息进行加权平均得到,这一过程通过核函数来实现.本文中使用的核函数如下所示

式中,r表示两个粒子之间的距离,re表示粒子的影响半径.式(3)所示的核函数是一种改进了的核函数,它既保持了原有核函数的简洁形式,又不存在奇点,避免了当两个粒子距离较近时出现的压力较大的现象,从而增强了压力场的稳定性.

在MPS 方法中,使用粒子数密度来表示粒子分布的疏密程度,其定义为粒子作用域内所有粒子的核函数之和,定义式如下

对于不可压缩流体,流体内部的粒子数密度是常数.

对于控制方程中压力梯度的离散,采用梯度模型来完成.在MPS 方法中,保持压力场的稳定性至关重要.流场压力的振荡不仅会带来计算误差,而且可能会导致计算的崩溃.因此,选择准确稳定的压力梯度模型至关重要.近年来,众多研究者按照不同的思路改进压力梯度模型[21-23],改进的总的原则就在于要保证粒子间的作用力始终为排斥力,从而保持动量守恒,保证计算的稳定性.本文中使用的压力梯度模型为Tanaka等[23]提出的梯度模型,表达式为

式中,D表示问题的维度,n0表示初始粒子数密度.张雨新[24]设计了一系列数值模拟试验,验证了式(5)所示的压力梯度模型可以较好地保证系统的动量守恒,改善压力控制方程中的速度散度项通过散度模型压力离散,散度模型表达式为

在MPS 方法中,粒子的压力通过求解压力泊松方程得到.因此,压力泊松方程模型的选择也关系到计算的稳定性和准确性.许多研究人员,尝试从不同的角度来对压力泊松方程进行改善[25-26],其中包括对拉普拉斯算子的改善以及对源项的改善.传统的MPS 方法中的泊松方程的源项完全由粒子数密度构成,但当流动较为剧烈时,粒子分布不均匀,必然会造成流场压力计算的不准确,导致流场压力振荡.于是又有部分学者提出,利用不可压流动中速度散度为0的特点,使用速度散度作为泊松方程源项.以速度散度为源项的泊松方程可以得到较为稳定光滑的压力场,但是由于计算速度散度时误差的积累,可能会导致流体体积的不守恒.为了综合以上两种泊松方程源项的优点,Tanaka等[23]提出了一个混合源项法,即速度散度和粒子数密度在计算中各占一定的比例,如式(7)所示

上角标k,k+1和*分别表示当前时间步、下一时间步和中间时刻,γ表示粒子数密度源项中所占的比例.Lee等[27]研究得出γ取值为[0.01,0.05]区间时,流场压力的稳定性和流体体积的守恒性较好.

在求解压力泊松方程时,将自由表面粒子压力为0作为边界条件.因此,准确地捕捉自由表面粒子对于稳定地求解压力泊松方程非常重要.在MPS 方法中,使用粒子数密度来判断粒子分布的疏密程度.因此可以利用自由表面粒子周围的邻居较少的特征,使用粒子数密度来判断自由表面粒子.在本文中,当 〈ni〉 ＜0.8n0时,n0为初始流体内部粒子数密度,认为该粒子位于自由表面,当 〈ni〉 ＞0.97n0,粒子则被判断为流体内部粒子.若仅以这种方式来判断,则在流动较为剧烈,粒子分布不均匀时,非常容易造成自由表面粒子的误判,导致在流体内部出现压力为0的位置,导致流场压力的振荡.因此,Zhang等[28]引入了一相对位置矢量函数F,利用自由表面处粒子两侧的粒子分布不均匀的特征来判断粒子是否位于自由表面上.相对位置矢量函数的定义如下

当〈 |F|〉i>0.9|F|0时,认为粒子位于自由表面.其中,|F|0表示初始时刻自由面粒子的 |F| 值.张雨新[24]通过数值模拟,验证了这种自由表面粒子判断方法可以有效地避免当流动剧烈时,将流体内部的粒子误判为自由表面粒子,从而提高了流场压力的稳定性.在处理固定边界时,为了保证边界附近粒子支持域的完整性,将边界设成一层第一类边界粒子和两层第二类边界粒子,如图1 所示.其中,第一类边界粒子的压力通过和流体粒子一起参与压力泊松方程的求解来得到.第二类边界粒子的压力通过第一类边界粒子和流体粒子的压力插值得到.这种多层边界粒子的方式既可以避免粒子的穿透,又可以保证边界附近流体粒子支持域的完整.

图1 固壁边界粒子示意图Fig.1 Schematic diagram of wall boundary particles

1.2 求解结构的FEM 方法

在这一小节中,介绍用于计算结构受力和变形的FEM 方法.FEM 方法的控制方程为动力平衡方程,即建立节点单元的平衡条件.控制方程可以写为

式中,M表示质量矩阵,K表示刚度矩阵,C表示阻尼矩阵.其中,为了简化计算,阻尼矩阵可以表示为质量矩阵M和刚度矩阵K的线性叠加,可以写为

式中,α1和α2分别表示瑞利阻尼系数.F(t)表示节点处随时间变化的外力矢量.本文采用Newmark-β[29]方法求解结构动力平衡方程,通过对速度和位移进行泰勒展开,可以得到在下一时间步t=t+Δt的节点位移,可以写为

式中,β和γ是保证计算稳定性的控制参数,需要满足 γ ≥0.5,β ≥0.25(0.5+γ)2条件,在本文模拟中,β=0.25,γ=0.5.将式(11)和式(12)代入控制方程中,就可以求解得到下一时间步的结构节点位移.求解流程表达式如式(13)～式(16)所示

1.3 MPS 与FEM 的耦合策略

在本文中,MPS 与FEM 求解模块之间的信息传递通过串行交错算法(conventional serial staggered,CSS)的分区耦合策略来实现.结构域和流体域分别在各自的区域内独立求解,之后在流固交界面处通过插值完成信息的传递.当结构单元的离散形式不同时,流固界面不同,插值的方式也不同.在本文中,采用核函数插值技术,利用粒子法的核函数具有的插值功能,实现了流固异构界面上的位移和载荷信息的传递.图2 中分别展示了流固界面之间的位移和载荷的插值技术.

如图2(a),流体粒子m的位移为

图2 流固界面插值技术Fig.2 The schematic diagram of coupling fluid-structure interface

式中,i表示以rei为插值半径的流体粒子m的插值域内的计算结构节点,δi表示节点i的位移.节点i的载荷为

式中,i表示计算结构节点插值域内的流体粒子,Pi为流体粒子i的流体压力,l0为流体域初始粒子间距.

2 数值方法验证

在这一小节中,将对本文中使用的数值方法以及采用的求解器MPSFEM-SJTU 进行验证.验证分为造波模块验证和弹性体变形验证两个部分.

2.1 数值造波模块验证

首先在三维数值水池中生成三维规则波.数值波浪水池模型如图3 所示,水池长L=12 m,宽D=1.5 m,水深h=1.8 m.水池左侧设置了造波板,通过推板造波的方式生成波浪.造波板的位移可基于线性造波理论推导得出,表达式为

图3 数值波浪水池示意图Fig.3 Schematic diagram of numerical wave pool

式中,H表示波高,k表示波数.在水池的右侧设置了消波区,消波区长l=3 m.消波是通过在方程中增加源项来实现的,源项可以写为

式中,xs是消波区的起始位置,Ls是消波区的长度,αs是一个无因次化的人工黏性系数,用于控制消波强度;u表示流体速度矢量.生成目标规则波的波长为λ=3m,通过Airy波的色散关系ω=以及T=2π/ω,求得规则波周期为T≈1.407 8 s.波高H=0.3 m,波陡 δ=H/λ=0.1.

取三种初始粒子间距对造波模块进行收敛性验证,分别为Δx=0.04 m,Δx=0.03 m和Δx=0.025 m,分别模拟规则波的生成和传播过程.在水池底部分别布置了波高测点,波高测点的横坐标分别为x=4.5 m,x=6.0m,x=7.5 m和x=11.0m.其中,x=11.0m 的测点位于消波区内,用以验证消波区的消波效果.四个测点处的波高时历曲线如图4 所示.

从图4 中可以看出,在前面三个测波点处,波高的时历曲线均与解析解吻合较好.当Δx=0.04 m 时,模拟得到的规则波的波高不是很稳定,与解析解相差较大.并且随着波浪的传播,在距离造波板较远的位置,波高略低于理论值,这可能是波浪在船舶运动过程中会出现能量的衰减.但随着初始粒子间距减小,时历曲线逐渐贴近解析解,证明了求解器的造波模块具有较好的收敛性.并且从图4 中也可以看出,对于粒子间距为Δx=0.03 m和Δx=0.025 m的两个算例,已经不存在明显的波高衰减的现象.由此可以说明,减小初始粒子间距可以较好地避免出现波高衰减的现象.对于位于消波区内的测点,从图3(d)中可以看出,波高测点的时历曲线非常平稳,证明了求解器消波模块的有效性.

图4 四个测点处的波高时历曲线Fig.4 Time history curves of wave height at four measuring points

图5为粒子间距为Δx=0.03 m 时,在一个周期内,波浪场的波面分布.从图5 中可以看出,生成的规则波波面具有较好的正弦特性,波面连续,无明显的非线性特征.沿着波浪的传播方向不存在波高的衰减现象.因此可以证明本文使用求解器可以较好地实现造波和消波功能.

图5 一个周期内规则波的波面分布情况Fig.5 Wave surface of regular wave in a period

2.2 弹性结构变形验证

在这一节中,模拟二维圆柱绕流后弹性悬臂梁的运动,验证求解器弹性结构在流体中的变形的准确性.

本文的计算模型按照Turek等[30]的模型实验尺度建立,如图6 所示.通道内有一个固定圆柱,圆柱后连接着一个弹性板.其中,通道高H2=0.41 m,通道长L2=2.5 m,圆柱直径为D=0.1 m.以通道左下角为坐标原点建立坐标系,圆柱圆心坐标为(0.2,0.2),圆柱后的悬臂梁厚H1=0.02 m,长L1=0.35 m,悬臂梁使用梁单元来离散,悬臂梁端点的坐标为(0.6,0.2),弹性悬臂梁的弹性模量为Es=5.6×106Pa.通道内来流平均速度为U=1 m/s,方向沿着通道方向,流体的黏度为μ=1.0×10-3m2/s.粒子间距取为Δx=0.002 5 m,D/Δx=40.

图6 计算模型示意图Fig.6 Schematic diagram of the calculation model

来流流经静止的圆柱后,会在圆柱后产生均匀规律的涡脱落现象.圆柱后的弹性悬臂梁受到涡脱落的影响产生受迫振动,其振动频率与漩涡脱落频率相同,周期约为T=0.8 s.图7 展示了在半个周期内,弹性悬臂梁的变形情况与其他数值方法得到的模拟结果的对比.

图7 t=t0,t=t0+T/8,t=t0+T/4,t=t0+3T/8 时刻弹性悬臂梁变形情况Fig.7 Deformation of elastic beam at t=t0,t=t0+T/8,t=t0+T/4,t=t0+3T/8

从图7 中可以看出,圆柱后的弹性悬臂梁,表现出了较好的弹性特征,变形符合物理规律,并与基于SPH[31]方法模拟得到的结果吻合较好.为进一步验证对于结构变形求解的准确性,记录了两个周期内悬臂梁尾端的位移随时间的变化,并与文献[30-31]得到的计算结果进行对比,如图8 所示.从图8 可以看出,MPSFEM-SJTU 求解器模拟得到的悬臂梁端点位移时历曲线,与文献[30-31]得到的数值计算结果在t=1.0T之前稍有差异,可能的原因是本文中选取的初始粒子间距较大,为D/Δx=40,而文献[30]选用的初始粒子间距为D/Δx=100.但时历曲线的振幅较为接近,可以证明本文使用的求解器在模拟弹性结构在流体中的变形时,具有较好的准确性和稳定性.

图8 悬臂梁端点位移的时历曲线Fig.8 The time history curves of the displacement of the cantilever beam endpoint

3 数值模拟

在本小节中,模拟三维船体在规则波中的运动.首先模拟了刚性船体在规则波中的运动,研究了规则波的波长对船体运动的影响.接着模拟了弹性船体在规则波中的运动,并与刚性船体的运动进行对比,研究了船体的弹性对船体运动的影响.

本文中选用的船模为船长3.2 米的kvlcc2 船,缩尺比为1:100.kvlcc2 船的实船尺度和模型尺度如表1 所示.数值水池的尺度与2.1 节中的造波水池一致,船艏与造波板的距离为3 m.

表1 kvlcc2 船实船和模型船参数Table 1 Characteristic parameters of the kvlcc2 ship

在对船体模型进行离散时,可以离散为全船有限元模型,也可以离散为非均匀梁模型.前者需要大量的建模工作和计算资源,而采用变截面梁模型可以大大简化计算,提高工作效率.在本文中选择非均匀梁模型,将船体沿船长方向离散为一系列变截面梁单元,船体模型与梁单元的等效关系如图8 所示.图9(a)展示了船体梁模型分组插值的方法,即将位于同一纵剖面处的船体粒子分为一组,这样一组粒子被等效为梁单元上的一个节点.从图9(b)可以看出,通过映射关系,将船体梁单元上的节点的位移和变形信息传递到船体边界粒子上,船体边界粒子的载荷信息也可以通过这样的映射关系传递到梁单元的节点上.

图9 船体粒子与梁模型的等效关系Fig.9 Equivalent relationship between hull particles and beam model nodes

船体模型的初始粒子间距为dp=0.025 m,根据船体模型的特征参数,将船体密度设为ρs=945 kg/m3.非均匀梁的横截面积、质量分布、横截面惯性矩在各个截面各不相同,在本文中的计算式为

式中,A,M,I分别表示横截面积、截面质量和截面惯性矩.下标i,m和S分别表示船体边界粒子、梁单元的节点和船体结构.图10为统计得到的船体梁截面质量和截面惯性矩沿船长方向的分布.将图10(a)曲线进行积分,可以得到统计的船体梁模型总质量为311.6 kg,与表1 中数据一致.由此可以说明本文中采用的离散方式的合理性.

图10 非均匀船体梁模型特征参数沿船长的变化Fig.10 Variation of characteristics along the length of a nonuniform hull beam model

3.1 模拟刚性船体在三维规则波中的运动

首先对船体在规则波中的运动响应进行收敛性验证,选取三组粒子模型,粒子间距分别为:Δx=0.03 m,Δx=0.035 m和Δx=0.04 m,三组计算模型对应的粒子总数分别为101 万和64 万和43 万.图11展示了初始粒子间距Δx=0.03 m 的模型,在一个波浪周期内船体的运动形态和自由面速度分布.

从图11 中可以看出,MPSFEM-SJTU 求解器可以较好地模拟三维船体在规则波中的运动情况.模拟得到的船舶运动形态以及自由表面粒子速度分布都符合物理规律.在规则波的作用下,船舶会发生周期性的运动,在本文中只考虑垂荡和纵摇的耦合运动.

图11 三维船体在波浪中的运动Fig.11 The ship motion in regular wave

图12 是10～ 20s 内,不同的粒子间距下船体垂荡和纵摇的时历曲线.从图中可以看出,在三组粒子间距下,船体的运动响应都表现出了较好的周期性,符合规则波激励下的响应特征.但粒子间距对船体运动响应影响较大.当Δx=0.04 m 时,船体的垂荡和纵摇运动都较小,这是因为当粒子间距较大时,粒子在运动过程中会出现较为明显的能量衰减.但当初始粒子间距减小时,衰减程度减小,验证了本文的算法具有较好的收敛性.关于船体运动响应时历曲线数据的分析如表2 所示.可以看出,粒子间距增大,船体运动响应幅值也会大幅度减小.从图4 中的分析可知,当粒子间距较大时,规则波的波高也会低于理论值.因此,此处的船体运动响应曲线振幅的差别不仅由粒子间距造成,规则波波高的不同也会影响船体的响应幅值.从图12(b)和表2 中可以看出,粒子间距对于模拟船体的纵摇运动影响较小.与垂荡响应曲线相比,三种粒子间距下船体纵摇的时历曲线相差较小.

表2 船舶运动响应参数Table 2 Ship motion response parameters

图12 船体在规则波中运动的时历曲线Fig.12 The time histories curves of the ship motion

接下来分析规则波的波长对船体运动响应的影响.在这一小节均采用初始粒子间距Δx=0.03 m 的模型来进行分析和计算,船型和水池尺度与上一小节相同.取波长分别为λ=2 m,λ=3 m,λ=3.2 m,λ=3.5 m,λ=4 m进行实验.已知船长L=3.2 m,即3 组实验的波长分别为λ=0.625L,λ=0.9375L,λ=L,λ=1.09375L和λ=1.25L,规则波的波长均在船长附近.在本文中计算得到的船体的垂荡和纵摇时历曲线如图13 所示.

由于本文波浪基于线性波浪理论,规则波的波长与频率成反比,所以随着波长的增大,规则波的频率降低.由于受迫振动的响应频率与激励频率相同,所以在相同波长规则波的激励下,船体的垂荡和纵摇响应频率相同.随着规则波波长增大,船体运动响应的频率降低.在不同的规则波激励频率下,船体运动响应的幅值也不同,统计了图13 中几个周期内的船体运动响应平均幅值,绘制了船体运动响应随波长的变化曲线,如图14 所示.

图13 不同波长下船体运动响应的时历曲线Fig.13 Ship motion response at different wavelengths

图14 不同波长下的船体运动响应平均幅值Fig.14 Motion response average amplitude at different wavelengths

从图14 中可以看出,当规则波的波长由2 m 增大到3 m 时,纵摇平均幅值增大,而垂荡平均幅值减小,当波长从 λ=3 m 继续增大时,垂荡和纵摇曲线的变化趋势大致相同.当波长 λ=3.2 m=L时,船体的垂荡运动和纵摇运动的幅值都最小.随着波长增大,船体的运动幅值逐渐增大.当波长 λ=4 m=1.25L时,纵摇和垂荡的幅度同时达到最大值,并且此时的运动响应幅值远大于其他波长下的平均响应幅值.由以上分析可以认为,对于本文中所研究的波长范围内(0.625L～ 1.25L),当规则波的波长大于船长并逐渐增大时,船体的运动响应也逐渐增大.当规则波的波长小于船长时,垂荡和纵摇曲线出现极值时对应的规则波波长可能不相同.

3.2 模拟弹性船体在三维规则波中的运动

在这一小节中,模拟弹性船体在规则波中的运动,水池尺度、船体模型以及造波模式等均与3.1节相同.图15 展示了弹性船体一个周期内的运动情况.

从图15 中可以看出,弹性船体在规则波的作用下发生了周期性的纵摇运动,但由于弹性船体的弹性模量较大,刚性较强,因此难以观察到明显的变形.图16 展示了在图15(c)情况下,即当波峰位于船舯附近时,船体表面的压力分布情况.

图15 弹性船体在规则波中的运动Fig.15 Motion of elastic hull in regular waves

从图16 中可以看出,与刚性船体相比,弹性船体在船舯部的压力分布颜色更深,压力更大,弹性船体的压力大于刚性船体.弹性船体表面最大压力为2393.52 Pa,刚性船体底部最大压力为2202.61 Pa,约为弹性船体最大压力的92%.并且较大的压力主要分布在船舯附近.在这一时刻下,造成这一现象的主要原因有以下两个,第一是波峰位于船舯附近,因此处水深较大,受到的水压力也较大.第二是弹性船体在这种波浪条件下会发生中拱,即使在船体弹性模量较大的情况下.船体梁的弯曲也会增大船底粒子受到的压力.

图16 弹性和刚性船体的表面压力分布Fig.16 Pressure distribution of elastic and rigid hull

为研究弹性船体和刚性船体在规则波中的运动响应的区别,分别绘制刚性船体和弹性船体的垂荡和纵摇的时历曲线,如图17 所示.

从图17 中可以看出,刚性船体的垂荡和纵摇运动幅值都大于弹性船体.8 s 之后,船体的垂荡曲线的幅值还存在着较小的波动,但船体纵摇的响应曲线幅值基本稳定.统计8 s 之后图17 所示的9 个周期内的船体运动响应的平均幅值,数据如表3 所示.

图17 弹性和刚性船体的运动响应时历曲线Fig.17 Time history curves of motion response of elastic and rigid hull

从表3 中的数据可以看出,刚性船体的纵摇平均幅值约为5.676°,弹性船体的纵摇幅值约为4.253°,相比刚性船体减小了25.1%.而对于垂荡响应,9 个周期内刚性船体的垂荡平均幅值为0.014 8 m,弹性船体的垂荡平均幅值约为0.013 1 m,相比刚性船体减小了11.5%.由此可见刚性船体在规则波中的运动比弹性船体的运动更剧烈,不考虑船体的弹性时,计算得到的船体运动响应会略大于实际的船体运动响应.

表3 弹性刚性船体运动响应参数对比Table 3 Comparison of motion response parameters of elastic and rigid hulls

接下来沿船长方向,分别取了5 个截面,分别位于船艏,距船艏Lpp/3 处,船舯,距离船艏2Lpp/3处,以及船尾.五个截面受到的压力随时间的变化曲线如图18 所示.

图18 刚性和弹性船体不同截面处的压力时历曲线Fig.18 Time history curves of pressure at different sections of rigid and elastic hull

从图18(a)中可以看出,在船艏处,弹性船体和刚性船体的受力情况有明显的不同.虽然二者的压力的时历曲线都具有明显的周期性和脉冲特性.但对于刚性船体,压力大多数情况下都表现为正值,其主要成分应该是浮力.而对于弹性船体,压力曲线在零附近有规律地波动,正值和负值同时存在.这种情况下,压力的组成成分主要来自于浮力和甲板上浪引起的砰击压力.从图15 中可以看出,刚性船体的纵摇幅度较大,在发生甲板上浪时,船艏已被完全淹没.因此,刚性船体的船艏处难以捕捉到由甲板上浪带来的砰击压力,这部分砰击压力可能会作用到船艏之后的某一截面上.从图18(b)～ 图18(d)可以看出,在1/3 船长、船中和2/3 船长这三个截面处,砰击载荷同样呈周期性变化,并且刚性船体与弹性船体的载荷曲线较为接近.对位于船尾截面受到的压力,刚性船体受到的压力幅值大于弹性船体.并且弹性船体和刚性船体都存在着一部分压力为零的时刻,这是因为此时船体艏倾,船尾被抬起,离开了水面,因此此时压力为零.仔细观察还可以发现,弹性船体压力为零的时间要短于刚性船体压力为零的时间,这同样是因为弹性船体的运动幅度更小.

4 结论

在本文中,采用MPS-FEM 耦合方法,使用自主开发的求解器MPSFEM-SJTU,模拟了刚性和弹性船体在规则波中的运动.首先在三维数值水池中进行空场造波,验证了求解器的造波模块和消波模块.接着通过模拟二维圆柱绕流后弹性悬梁的振动,验证求解器模拟弹性流固耦合问题的准确性.验证了求解器后,本文模拟了弹性船体和刚性船体在规则波中的运动,得到了以下结论.

(1) 船体在规则波中的纵摇运动比垂荡运动更稳定.当刚性或弹性船体在规则波中的运动响应发展到较为稳定的阶段时,船体纵摇运动的幅度基本保持不变,但船体垂荡运动的幅值还有轻微的波动.

(2) 弹性船体的运动幅度比刚性船体小,从模拟得到的弹性和刚性船体的运动响应时历曲线可以看出,无论是垂荡运动还是纵摇运动,刚性船体的运动幅度都较大.不考虑船体弹性时计算得到的运动响应可能偏大.

(3) 由于弹性船体会发生弯曲,因此在船体底部,弹性船体船舯附近的压力大于刚性船体.而由于刚性船体的运动幅度大于弹性船体,所以刚性船体的船艏和船尾处受到的砰击压力大于弹性船体.