基于离散模块梁单元理论的超大型浮体在非稳态载荷下的水弹性响应

陈永强, 宋 炜, 张显涛

(1. 上海交通大学 海洋工程国家重点实验室, 上海 200240;2. 中国水产科学研究院 东海水产研究所, 上海 200090;3. 上海交通大学 三亚崖州湾深海科技研究院, 海南 三亚 572024)

0 引 言

大型浮式结构物是一种水平尺度远大于垂向尺度的特殊海洋工程结构物，一般依托于沿海岛屿设置，将对其附近区域产生巨大深刻的经济和军事影响[1]。大型浮式结构物的动力响应问题是一种水弹性问题，即一种惯性力、水动力和弹性变形力相互作用的力学行为。模态叠加法和直接法是处理水弹性问题的两大主要方法。模态叠加法需要先对结构进行模态分析，确定最佳模态组合后对各模态进行叠加[2-3]。直接法则针对结构和流域均进行单元划分，不借助模态函数直接求解水弹性方程，计算成本极大[4]。

近年来，一种全新的离散模块-梁单元(Discrete-Module-Beam，DMB)水弹性方法[5]逐渐应用于海洋工程中。DMB方法将浮体划分为若干子模块，并引入梁单元连接各子模块重心，结合三维线性势流理论和梁的弯曲理论进行水弹性分析。诸多学者使用DMB方法解决若干海洋工程问题。ZHANG等[6]和SUN等[7]应用DMB方法准确预测某铰接超大型浮体在规则波下的水弹性响应，展示采用DMB方法求解具有复杂连接形式的大型浮式结构物动力响应的便捷性；陈永强等[8]在DMB框架中提出针对铰接的建模方法并求解铰接超大型浮体在规则波下的水弹性响应；BAKTI等[9]将前进速度因素融入DMB框架中，研究其对水弹性响应的影响；JIN等[10]考虑某系泊大型浮式结构物，使用DMB方法解决连接体-浮体-锚链耦合系统的水弹性响应问题。

本文基于DMB时域分析模型，针对学界较为知名的MF-300超大型浮体，计算其在非稳态载荷(货物移动、飞机降落)下的水弹性响应，给出MF-300的位移响应和弯矩响应，方便其他学者进行对比。

1 离散模块梁单元方法

1.1 水弹性方程建立

以箱式超大型浮体为例说明DMB水弹性方法的处理思路。如图1所示，某箱式超大型浮体受到入射角为φ的波浪作用，浮体几何尺寸为L×b×d，吃水为T。

图1 离散模块梁单元方法分析示例

DMB水弹性分析理论主要包含水动力分析和结构分析两部分。在水动力分析方面，将超大型浮体划分为若干子模块，将每个子模块视作刚体，基于三维势流理论进行多刚体水动力分析，求解附加质量系数矩阵A(ω)、辐射阻尼系数矩阵B(ω)、静水恢复力系数矩阵C和波浪激励力矩阵FE。在结构分析方面，将每个子模块抽象为一个位于相应重心处的集中质量，引入线性欧拉-伯努利梁连接相邻集中质量，基于梁弯曲理论建立节点受力平衡方程，得到DMB方法框架下的频域水弹性方程：

{-ω2[M+A(ω)]-iωB(ω)+(C+Kst)}ξ=FE

(1)

式中：ω为入射波频率；M为结构质量矩阵；Kst为结构总刚度阵，ξ为集中质量的六自由度位移。需要注意的是，N个子模块对应N个集中质量，但为了后续刚度阵的推导，定义左侧自由端为集中质量0，右侧自由端为集中质量N+1。

基于矩形截面规则均匀梁刚度阵的理论解，按照有限元标准程序组装总结构刚度阵Kst(6N×6N)，具体细节参见ZHANG等[6]。设梁单元e连接集中质量i与j，则梁单元两端的结构变形力Fst，i和Fst，j与2个集中质量的位移ξi和ξj可用梁单元e的结构刚度阵ke连接：

(2)

式中：ke，(i，i)、ke，(i，j)、ke，(j，i)和ke，(j，j)为其子矩阵，维度均为6×6，其中，0≤i

图2 梁单元e0示例

因此，式(2)变为

(3)

将式(3)展开可得

0=k0，(0，0)ξ0+k0，(0，1)ξ1⟹ξ0=
-(k0，(0，0))-1k0，(0，1)ξ1

(4)

Fst，1=-k0，(1，0)ξ0-k0，(1，1)ξ1

(5)

将式(4)代入式(5)计算可得：

(6)

同样，对于梁单元eN可得：

(7)

ZHANG等[6]基于Cummins方程将DMB方法扩展至时域，提出处理非稳态集中力的准静态位移增量法。假设超大型浮体受集中力fc(t)作用，则超大型浮体的时域运动方程为

(8)

式中：A(∞)为无穷大频率下结构的附加质量；Γ(t)为时延函数，表达流体的记忆效应。需要注意的是，式(8)左侧是针对所有子模块重心建立的矩阵方程，而右侧的集中力仅为实数。因此式(8)使用“←”符号连接，表示右侧激励引起左侧的水动力响应。

ZHANG等[6]提出准静态位移等效法，假设作用在每个集中质量上的等效分布力FUN(t)在任意时刻引起的水弹性响应与集中力fc(t)相同，这一转变使式(8)变为

(9)

图3 集中力等效为均布载荷的过程示例

(10)

(11)

(C+Kst)Δξs(tn)←fc(tn)

(12)

(C+Kst)Δξs(tn)=FUN(tn)

(13)

(14)

1.2 位移响应还原

图4为离散之后的超大型浮体结构，结构被划分为N个子模块，由此产生N+1个交界面；每个子模块被抽象为一个集中质量(图4中节点)；每个梁单元均具有左右2个端面，记为L(left)和R(right)。

图4 大型浮式结构物离散模型示例

求解式(1)或式(9)可获得所有集中质量位置的位移响应，采用子梁法还原梁单元上任意位置处的位移响应。如图5所示，任一给定位置xany位于梁单元e上，该梁单元连接集中质量i和j，具有端面p和q，集中质量i(j)与端面p(q)位移相同。位置xany与端面p构成子梁f，根据子梁f的刚度阵可构建以下方程：

(15)

图5 子梁法还原位移响应示例

进一步运算可得：

ξany=(kf，(p，any))-1(Fp-kf，(p，p)ξp)

(16)

Fany=kf，(any，p)ξp+kf，(any，any)(kf，(p，any))-1·

(Fp-kf，(p，p)ξp)

(17)

式(16)即给出了任意位置处的位移响应。

1.3 弯矩响应还原

Kstξ(t)=fc(t)+Fa(t)+Fadm(t)+
Fdamp(t)+G+B(t)

(18)

将超大型浮体视作无质量、无支座支撑和弹性基础的两端自由梁，将外载荷加载过程中所有作用在结构上的力均视作外力。求解结构力响应的本质就是对外载荷进行积分，当超大型浮体受到非稳态集中力时，外载荷可分为以下2个部分：

(19)

式中：q1(t)的分布可准确求解，q2(t)则相反。大型浮式结构物的剪力Nstr和弯矩Mstr分布即为两者引起的叠加：

(20)

Nq1(t)和Mq1(t)的求解步骤如下。如图6所示，以ΔL为间距沿浮体纵向取N″个位置点，每个位置点的位移均可通过式(16)求解。设ξj(t)=[ξj，1(t)ξj，2(t)ξj，3(t)ξj，4(t)ξj，5(t)ξj，6(t)]T为位置点j的位移，则作用于位置点j的重力Gj(t)和浮力Bj(t)为

(21)

图6 Nq1(t)和Mq1(t)的求解思路示例

式中：ρstru为浮体密度；ρwater为海水密度。式(21)给出作用于任意位置点的重力和浮力，此外，假设重力和浮力在相邻位置点间呈线性变化。当ΔL变小时，即可得到重力分布曲线和浮力分布曲线，非稳态集中力fc(t)作用于xc处。

基于结构任意位置处位移求解所受浮力B(t)，结构重力G(t)分布显然已知，因此两者在任意位置处引起的剪力Nq1，any(t)和弯矩Mq1，any(t)为

(22)

Nq2(t)和Mq2(t)的求解步骤如下。对于首端梁单元e0，q2(t)在位置点1必然引起剪力弯矩响应，但由于位置点1为自由端，q1(t)和q2(t)在该位置引起的结构力之和必为零。位置点N″亦是如此。

结合式(20)可得：

(23)

相应地，式(3)变为

(24)

可得：

(25)

相似地，

(26)

(27)

(28)

M(x)=a+bx+cx2+dx3

(29)

则必有

(30)

求解式(30)中4个待定系数，进而求解任意位置处的弯矩。最后，按式(20)进行叠加即可获得超大型浮体在非稳态集中力作用下任意时刻的结构力分布。

2 模型验证与对比

研究对象为日本学者ENDO等[11]提出的MF-300超大型浮体，其真实尺寸如表1所示。众多学者对MF-300进行详细的试验研究和数值研究。在使用DMB方法时，将浮体分为8个子模块建立水弹性方程；将浮体划分为300个子模块实施准静态位移等效法求解等效分布力FUN；采用4阶龙格-库塔法求解时域水弹性方程(9)(具体见ZHANG等[6])。

表1 MF-300真实尺寸及环境参数

2.1 货物移动

图7为ENDO等[11]进行的货物移动试验，试验缩尺比为30.77。模型先静止于静水面，以0.125L为间距共取9个测点Z1～Z9(图中未完全展示)，记录货物移动过程中Z1、Z5和Z7等3个测点的垂向位移变化。一个0.9 kg重物从0.19L处，以0.61 m/s的速度均匀移动至0.9L位置。

图7 货物移动试验示例

货物移动过程中的位移模拟值与试验值对比如图8所示。总体来看，DMB的模拟结果与试验值吻合较好，但是与试验值相比呈现出较大的波动，一个可能的原因是试验值进行滤波处理，过滤了高频成分。

图8 货物移动工况下试验值与模拟值对比

在弯矩响应方面，目前未见学者给出该工况下弯矩响应情况，因此仅给出DMB方法的计算结果。图9给出不同时刻MF-300上的弯矩分布曲线，水平辅助线表示0刻度线，竖直辅助线表示重物位置。由图9可知，弯矩分布曲线在两自由端为零，且与水平辅助线相切(剪力为零)。弯矩在重物位置存在突变，这也是合理的。

图9 货物移动过程中不同时刻弯矩分布曲线

2.2 飞机降落

在计算飞机降落工况时，采用以下假设[12]：

(1)将飞机视作沿浮体纵向中线移动的质点，对浮体仅有垂向力作用；

(2)在降落过程中飞机所受升力与速度的平方成正比；

(3)不考虑飞机与浮体的耦合作用，忽略由飞机自身垂向运动引起的垂向力变化。

基于上述假设，飞机降落载荷被理想化为一个作用位置与载荷大小均随时间变化的竖直方向的集中力。

图10给出飞机降落过程：超大型浮体MF-300的飞机跑道设置于纵向中部位置处，且在飞机降落前保持静止；飞机重达3 t，以150 km/h速度降落于浮体上0.19L位置处，滑行过程中的加速度为-5.79 m/s2；计算可得飞机滑行时间为7.2 s，滑行距离为150 m，最终静止于浮体0.69L处。基于假设，飞机降落载荷如图11所示。需要强调的是，在ENDO[12]的计算中，跑道平行于纵向中线且距离其3.3 m。

图10 飞机降落过程示例

图11 飞机降落载荷时历(基于本节假设，将飞机降落载荷理想化为垂向集中力)

图12给出不同时刻下浮体位移响应曲线的计算值与ENDO[12]结果的对比，两者吻合度非常好。其中竖直细实线表示飞机位置，注意飞机在7.2 s时刻停留在0.69L(207 m)位置处并保持静止。

图12 飞机降落过程中不同时刻位移分布曲线对比

图13给出不同时刻(对应图12)下的弯矩分布曲线，与货物移动工况类似，DMB方法计算的结果反映了集中力作用引起的弯矩突变。其中，水平细实线为零刻度线，竖直细实线为飞机位置。

图13 不同时刻(对应图12)弯矩分布曲线对比

3 结 论

使用DMB水弹性时域分析方法计算超大型浮体MF-300在货物移动和飞机降落两种非稳态激励下的水弹性响应，尤其是提出改进后的三次曲线插值法计算弯矩响应分布。改进方法可体现集中力作用引起的弯矩突变，且满足弯矩分布曲线在2个自由端函数值和一阶导数值(即剪力)均为零，拓展三次曲线插值法的应用范围。所得位移响应和试验值与其他学者模拟值吻合较好。本文也给出超大型浮体在非稳态加载过程中的弯矩分布曲线，供其他方法对比。