系泊系统耦合动力分析不同解法的比较研究

丁佐鹏，马 山，段文洋，刘 昊

(哈尔滨工程大学船舶工程学院，黑龙江哈尔滨150001)

深海浮式平台总体分析中，系泊缆索和立管的水动力阻尼效应对平台低频运动有重要影响，同时浮体摇荡运动对系泊线诱导动态张力作用，需采用动力耦合分析方法同时求解浮体和系泊缆索的动力响应［1－2］.由于系泊浮体动力耦合分析涉及的学科领域较多，力学机理和分析技术复杂，尚有诸多问题需进一步探索和完善，其中耦合分析效率和算法的稳定性研究便是一个亟待解决的关键问题.为探索在浮体和系泊/立管耦合动力分析中提高计算效率，文中采用“弱耦合”模型，基于异步耦合动力分析［3］的求解思想，开发了深海浮体与系泊、立管的动力耦合程序，在单点系泊FPSO，SPAR平台动力响应模拟上，其结果同商业软件Aqwa获得了很好验证［4－5］.同时在数值模拟中发现针对浮体与系泊/立管动力耦合分析过程，通常柔性构件的稳定求解需要较小的时间步长，导致泊缆索/立管动力响应分析相对于浮体动力求解耗时最多.但是不能为了节约计算时间而随意增大柔性杆件求解步长，这样容易影响算法稳定性，因此有必要对系泊线/立管动力响应分析方法进行深入研究以提高耦合计算效率.

细长杆理论用于系泊线/立管动力响应分析，最早由 Garret［6］提出，然后 Webster等［7］对其改进，考虑杆件小拉伸变形效应.基于细长杆理论进行缆索动力分析时常用的有Adams-Moulton单步积分法［8］和 Newmark-β 迭代法［9］.文中针对这两种方法，基于细长杆理论，编制了有限元数值求解的动力响应分析程序，分析了两种算法的数值精度，在此基础上开展了两种算法计算稳定性和计算效率的对比分析.

1 细长杆理论模型简介

1.1 缆索控制方程

图1所示系泊线典型杆单元的运动控制方程可以写成如下形式:

图1 缆索分析空间坐标系示意Fig.1 Coordinates for cable analysis

式中:B为抗弯刚度;r(s，t)为缆索中心线位置向量，s为弧长变量，t为时间变量;λ为拉格朗日乘子，λ=T-Bκ2，κ 为曲率，T(s，t)=F·r′为轴向张力;ε为应变.借助Galerkin法对缆索控制方程中的高阶偏导数降阶，同时利用插值函数将沿杆单元的变量用节点变量近似表示，可得如下所示的控制方程:

式中:M为质量矩阵;U为节点位移坐标;q为外载荷;γikm，αikm，βikm，μim，τm，ηlm为积分常数;i=1，2，3，4;k=1，2，3，4;m=1，2，3;n=1，2，3.

1.2 水动力载荷和海底边界条件的模拟

作用在缆索单元上的力主要有附加质量力qI，拖曳力qD和Froude-Krylov力qF-K.

上式中的前两项可通过莫里森公式求得，在知道了水质点的加速度之后，Froude-Krylov力也可求得.

海底刚度对杆单元的刚度矩阵有影响，同时海底支持力也应该作为单元外载荷的一部分来考虑.文中通过在海底布置弹簧来模拟海底弹性基础的作用，海底支持力为:

式中:c为海底弹簧刚度系数，文中取c=w/d2，w为缆索单位长度浸水重量，d为缆索直径;r2为缆索单元垂向坐标，h为水深，e2为垂向单位向量.

2 两种不同的动力学算法分析

2.1 Adams-Moulton 积分法

将式(2)中的二阶微分方程用关于变量的两个一阶微分方程表示:

对上式在时间步内进行积分，并利用梯形公式进行积分运算:

通过联立求解方程(7)和拉伸控制方程，即可求得节点位移和张力.由数值分析的知识可知，梯形公式的截断误差满足以下关系:

由式(8)可知，Adams-Moulton法的截断误差为O((Δt)3)，具有二阶精度，可以满足工程计算需要.

2.2 Newmark-β 法

文中采用平均加速度方法，假设加速度在每个时间步内为常数，其值为步长始末两点的平均值，在步长内通过积分获得速度和位移的表达式［10］.在初始时刻，运动控制方程如下:

第k步的位移和速度表示为:

M(k)和q(k)的值可以通过u(k)，·u(k)和u¨(k)的值计算得到.

将式(10)代入式(9)，经整理可将第k步的增量方程变成如下形式:

将其与拉伸控制方程联立求解可获得第k步内的位移增量δU(knk)和拉格朗日乘子增量δλ(mk).取γ=1/2，β=1/4，并利用加速度为常数的假设可将式(10)中位移写成如下形式:

可见，平均加速度假设下的位移具有二阶泰勒展开形式，即具有二阶精度.平均加速度法得到的加速度计算结果在节点处不连续，但是它是无条件稳定的.文中使用其迭代形式，即将前一步的结果作为初值代入控制方程(11)进行迭代求解直到δu(k)和δλ(k)足够小，停止迭代，进行下一时间步的计算.初步比较 Adams-Moulton积分法和Newmark-β法，可知在进行缆索动力响应步进求解时，两种方法具有相同的计算精度.下面将两种方法分别应用到实际海洋工程系泊系统的动力响应求解中，对两者结果进行比较.

2.3 两种算法计算精度对比

本算例中分别利用Adams-Moulton法和Newmark-β迭代法研究了顶端受迫振动情况下系泊缆索顶端张力的时历变化，并将结果与商用软件Orcaflex［10］进行对比.系泊缆索为钢链，工作水深223.5 m，链长762 m，名义直径0.14 m，抗拉刚度2.21×109N，单位长度质量235.2 kg/m，法向附加质量系取2.8，法向拖曳力系数取3.2，顶端预张力1558.8 kN.系泊线顶端受迫运动形式如下:

式中:x0为缆索顶端初始横坐标值，t为时间变量，G(t)为缓载函数.文中为如下形式:G(t)=3α2-2α3，α=t/Tload，Tload为缓载区间，a为受迫运动幅值(m)，ω为受迫运动频率(rad/s)，T为受迫运动周期(s).图2，3给出了文中采用的两种算法计算得到的缆索顶端动态张力结果与商业软件Orcaflex结果的比较.

图2 振幅0.9144 m、周期4 s下顶端张力对比Fig.2 Comparison of top tension(a=0.9144 m，T=4 s)

图3 振幅2.75 m、周期8 s下顶端张力验证Fig.3 Comparison of top tensions(a=2.75 m，T=8 s)

从图2，3中可以看出对不同振幅和频率下的受迫振动Adams-Moulton法和Newmark-β迭代法所进行的系泊缆索张力响应分析结果与Orcaflex的结果具有很好的一致性，从而证明在实际的动力响应求解过程中，这两种方法也具有相似的计算精度.

2.4 动力学算法稳定性和计算效率比较研究

在深海浮体动力耦合响应分析时，动力分析时间步长是一个关键参数，其对动力分析计算效率和稳定性有显著影响.下面通过改变时间步长研究两者收敛性和计算效率的比较规律.

选取某系泊线，其材料和工作水深见表1，EA为抗拉(压)刚度，锚链和钢缆的附加质量系数取1.0，拖曳力系数取2.0，按照式(13)分别选取了小幅高频和大幅低频两种典型的受迫运动形式.

图4，5分别给出了在幅值为1 m，振动周期为5 s以及幅值为4 m，振动周期为20 s两种简谐受迫运动情况下，两种动力分析算法获得的缆索顶端张力随时间步长变化的结果.

表1 系泊线参数Table 1 Paramaters of the mooring line

图4 小幅高频运动时收敛性比较(振幅1 m，周期5 s)Fig.4 Convergence analysis of the two algorithms under different time step(a=1 m，T=5 s)

图5 大幅低频运动时收敛性比较(振幅4 m，周期20 s)Fig.5 Convergence analysis of the two algorithms under different time step(a=4 m，T=20 s)

可以看出，在保证计算结果收敛的前提下，改变动力分析的时间步长对结果的影响很小，说明这两种算法都具有很好的收敛性.另外，比较两个算例最大时间步长可以发现，Newmark-β法适用的步长范围明显比Adams-Moulton法大.通过分析可知，由于每次迭代过程相当于将结果重新代入平衡方程对加速度的结果进行误差修正，减小了误差累积，稳定性因而更高.以此为基础，进一步研究了在不同的时间步长下系泊系统动力响应的数值求解效率.

文中算例通过模拟相同时间的缆索动力响应过程，对两种算法计算效率进行比较(图6，7).

图6 不同时间步长下两种算法计算效率比较Fig.6 Comparison of computational efficiency between the algorithms under different time steps

图7 不同时间步长下两种算法计算效率比较Fig.7 Comparison of computational efficiency between the algorithms under different time steps

从图6，7中可看出随着时间步长的增大，Adams-Moulton法的数值计算时间显著降低，而Newmark-β法则未能表现出与Adams-Moulton法相似的规律，甚至在一些大步长下计算效率反而下降了.通过仔细研究发现，迭代过程其实对Newmarkβ法的计算时间影响很大.尤其是在振动初始阶段，由于结构的响应比较大，为了满足精度要求往往在一个时间步需要迭代几十次，即使到了稳定计算阶段，每步也需要数次的迭代，从而影响了整体的计算效率.

3 结论

文中利用细长杆理论建立了系泊缆索的有限元模型，分别根据Adams-Moulton法和Newmark-β法这两种动力学算法编写了相应的缆索动力学求解程序.通过与商业软件的结果进行对比，表明在考虑了实际环境下的各种非线性因素(系泊线几何参数，水动力载荷和海底接触条件)之后，这两种算法依然能够得到满足工程应用的结果，具有较好的数值精度.

通过算例表明，Newmark-β法能适应较大时间步长，算法稳定性更好;Adams-Moulton单步法的计算效率随时间步长的增大提高很快，而增大步长对Newmark-β迭代方法的效率影响不大，步长太大了效率甚至会下降.分析发现，增大Newmark-β迭代法时间步长之后，为获得稳定的结果每一步所需的迭代次数增多，每一步的计算时间就增加了，从而影响了整体的计算效率.总体来讲，Newmark-β迭代方法算法稳定性好，更适合环境载荷瞬息万变的海洋工程柔性杆件分析，但需要进一步研究如何在保证计算精度的前提下减小每一时间步的迭代次数，以提高求解效率.

References)

［1］ Ormberg H，Larsen K.Coupled analysis of floater motion and mooring dynamics for a turret moored ship［J］.Applied Ocean Research，1998，20:55 －67.

［2］ 窦培林，朱建，陈刚，等.风浪流成长全过程单点系泊FPSO运动低频响应极值研究［J］.江苏科技大学学报:自然科学版，2014，28(2):114 －119.Dou Peilin，Zhu Jian，Chen Gang，et al.Low frequency extreme response of single-point moored FPSO during the whole process of growing of wind/wave/current［J］.Journal of Jiangsu University of Science and Technology:Natural Science Edition，2014，28(2):114 －110.(in Chinese)

［3］ Jing X，Webster W C，Xu Q，et al.Coupled dynamic modeling of a moored floating platform with Risers［C］∥Proceedings of the ASME30th Inter:Conference on Ocean，Offshore and Arctic Engineering.OMAE2011 －49553.Rotterdam，Netherlands:ASME，2011.

［4］ Ma Shan，Duan Wenyang.Asynchronous coupling analysis of dynamic response for deepwater moored floating platform［C］∥Proceedings of the Eighth International Workshop on Ship Hydrodynamics.IWSH2013 －P3B －3.Seoul，Korea:September，2013:23 －25.

［5］ 马山，段文洋.深海移动式平台与系泊/立管系统动力响应异步耦合分析研究［C］∥第二十五届全国水动力学研讨会暨第十二届全国水动力学学术会议论文集北京，2013:786－792.

［6］ Garret D L.Dynamic analysis of slender rods［J］∥Journal of Energy Resource Technology，1982，104(4)，302－306.

［7］ Pauling J R，Webster W C.A consistent larger amplitude analysis of the coupled response of a TLP and tendon system ［C］∥5th OMAE Conference，Tokyo，Japan，1986，3:126 －133.

［8］ Ran Z.Coupled dynamic analysis of floating structures in waves and current［D］.Texas:Texas A＆M University，2000.

［9］ Chen Xiaohong.Studies on dynamic integration between deep-water floating structures and their mooring/tendon systems［D］.Texas:Texas A＆M University，2002.

［10］ 李鸿晶，王通，廖旭.关于 Newmark-β法机理的一种解释［J］.地震工程与工程振动，2011，31(2):55－62.Li Hongjing，Wang Tong，Liao Xu.An interpretation on Newmark beta methods in mechanism of numerical analysis［J］.Journal of Earthquake Engineering and Engineering Vibration，2011，31(2):55 －62.(in Chinese)