系泊缆索动力分析中Newmark-β迭代解法的计算效率

， ， ，

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

0 引 言

在深海系泊浮体分析中，系泊缆索和立管的动态响应对浮体有很大影响，需采用动力耦合方法求解浮体与系泊缆索及立管的相互作用。系泊浮体耦合运动包括3个方面：浮体水动力建模、系泊缆索力学分析建模和系泊缆索的耦合算法。按照JACOB等[1]的讨论，耦合算法有“强耦合”方法和“弱耦合”方法。“强耦合”方法将浮体和系泊线动力学方程通过弹簧反力(矩)和阻尼力(矩)耦合在一起，形成一个浮体和细长杆件动力学系统同时求解的多自由度方程组，在这种方法中浮体和系泊缆索是整体求解的。“弱耦合”方法将浮体动力学模块和柔性杆件动力学模块独立求解，两者通过位移和受力的传递进行耦合。RAN[2]和GARRET[3]分别开发了时域全耦合动态分析程序。KIM等[4]将聚酯缆材料的非线性数学模型融合进系泊缆索的运动和控制方程，并与深海平台进行耦合分析。ZHANG等[5]分别采用准静态耦合、半耦合和全耦合3种方法对某Spar平台进行数值模拟，并将结果与试验进行对比。袁梦[6]研究了聚酯缆的动力特性，并对某张紧式Spar平台进行耦合动力分析。杨敏冬[7]对风、浪、流共同作用下的Truss Spar平台分别进行了静态耦合模拟和动态模拟。基于JING等[8]提出的异步耦合方法，马山等[9-11]开发了深海浮体与系泊、立管的动力耦合分析程序，在对单点系泊FPSO、Spar平台的动力响应分析中得到了很好的验证。该方法属于“弱耦合”方法，由于系泊系统作为独立模块求解，适合采用并行计算。

在浮体动力耦合分析数值模拟中系泊缆索、立管的动力分析时间步长较小，耗时较长，因此，提高系泊缆索动力分析的计算效率对于减小系泊浮体耦合分析的计算耗时至关重要。JACOB等分别采用“强耦合”和“弱耦合”方法对某半潜式平台进行动力分析对比，发现“弱耦合”方法结合并行计算能大幅提高计算效率。丁佐鹏等[12]基于细长杆理论将Adams-Moulton法和Newmark-β迭代法两种算法用于系泊缆索动力分析，并对这两种算法进行对比研究工作，发现在相同时间步长的情况下Adams-Moulton法计算效率更高，但是其只适用于小步长，而Newmark-β迭代法稳定性更好，适用于大步长。

在采用系泊浮体与动力分析异步耦合算法时，针对系泊缆索、立管动力分析常采用细长杆理论。该方法最早由GARRET[13]提出，适用于不可拉伸细长杆件，后经MA等[14]改进，考虑了杆件的拉伸变形。唐友刚等[15]基于细长杆理论建立了三维非线性细长杆单元刚度矩阵，实现有限元软件ABAQUS对缆索单元的调用，该方法非常适用于工程实际。马刚[16]结合细长杆理论分析了某单点系泊FPSO在波浪中的运动和受力。袁梦等[17]对控制方程进行修改，将细长杆理论应用于深水轻质缆索的动力分析。细长杆理论动力求解中经常采用的一种方法是Newmark-β迭代法[18]，本文采用该方法，以提高系泊缆索动力分析计算效率为目标进行以下工作：(1)研究单元数、时间步长、控制精度和每个时间步的迭代次数等要素对计算效率的影响。(2)将OpenMP并行计算应用中系泊缆索的动力分析中，研究其对计算效率的影响。发现适当减小单元数、增加时间步长、降低动力分析迭代求解控制精度以及减小每个时间步的迭代次数有利于提高计算效率，采用并行计算对提高计算效率也是有利的。

图1 缆索分析空间坐标示例

1 基于细长杆理论的Newmark-β迭代法

1.1 基于细长杆理论的动力学响应控制方程简介

图1给出了缆索动力分析的弧长坐标系，系泊缆索的运动及控制方程如下：

r′·r′=(1+ε)2(2)

引入插值函数，利用Galerkin法离散后的有限元运动及拉伸控制方程如下：

μimqmn+fin(3)

1.2 Newmark-β迭代法求解格式

系泊缆索有限元离散动力学方程组采用Newmark-β迭代法求解。具体求解格式如下：

在t=0时，或者时间步k=0时，有

(5)

假设初始时刻结构静止，初始条件为

(6)

(7)

时间步k时各变量的值用前一时间步k-1步的值通过式(9)进行预估：

(9)

式中：Δt为时间步长。

M(k)，q(k)可以通过以上各值求得。

式中：ρ为水的密度;A为材料截面积;FEA为缆索抗拉刚度。

求解出以上方程后，通过式(12)修正各变量的值：

(12)

式中：γ、β为常数，γ=1/2，β=1/4。

1.3 海底边界条件

海底对系泊缆索同时有支持力和摩擦力的作用，本文忽略摩擦力的作用，只考虑支持力的作用。将海底对缆索支撑力作用看成线弹性弹簧，参考文献[18]，支持力的公式为

(13)

参数运动及控制方程，对支持力进行离散积分后得到其表达式：

(14)

对于完全触底的单元，式(14)中的l1即为单元长度，而μim、γikm仍为积分常数；但是对于一部分触底、一部分悬空的单元，式(14)中的l1不是单元长度，μim、γikm也不再是常数，积分只在触底的那一部分区间上进行。

对于部分触底的单元需先求解出触地点，即求解式(15)中的零点：

(15)

2 各要素对计算效率的影响

基于以上理论，本文编制了系泊缆索动力求解程序。现以单一成分系泊线的强迫运动动力分析为例，研究有限单元数、时间步长、控制精度和每个时间步的迭代次数对计算效率的影响。系泊线的具体参数见表1。

表1 单一成分系泊线参数

图2 单一成分系泊线强迫运动示例

系泊线作业水深1 000 m，在系泊线顶端施加一水平强迫正弦运动，如图2所示，运动位移关系：

x(t)=x0+a·sin(ωt)(16)

式中：x0为系泊线顶端初始时刻坐标；a为施加的强迫运动的幅值；ω=2π/T为振荡圆频率，T为振荡周期。

坐标系原点在静水面上，取x0=0，a=10 m，T=15 s，系泊线顶端位于水平面下8 m处。

2.1 单元数和时间步长的影响

现将整个系泊线分别划分为10、20、40个单元，时间步长分别取0.1 s、0.2 s、0.3 s，共9种工况进行数值模拟，控制精度取10-2，并将系泊线顶端张力的时历曲线与商业软件OrcaFlex的计算结果进行比较，如图3～图5所示。

图3 单元数为10的计算结果 图4 单元数为20的计算结果

由图3～图5可以看出:划分不同的单元数，取不同的时间步长，采用Newmark-β迭代法进行数值计算得出的结果总体与OrcaFlex模拟结果是一致的，说明本文采用的数值方法是可靠的。上述比较也说明：在相同单元数目情况下，不同时间步长下计算结果较为一致，本文采用的Newmark-β迭代法数值稳定性较好，在较大的时间步长下缆索动力分析也能够较为稳定地进行下去。

图5 单元数为40计算结果

现按以上各工况进行数值模拟，数值模拟时长均3 000 s，得到计算耗时见表2。

表2 各工况模拟3 000 s的计算耗时对比 s

由表2可以发现：在固定控制精度为10-2，单元数相同时，随着时间步长的增加，计算耗时逐渐减小，随着单元数的增加，这一趋势越明显。在时间步长相同时，单元数越多，计算耗时越长，而且随着单元数的增加，耗费时间的增加速度越快。从这2点可知，单元数的划分对模拟的耗时非常敏感，虽然增加单元数可以提高计算精度，但是单元数增加到一定程度时，数值计算的精度已不能提高很多，如在单元数取20时计算精度已经很高，增加到40并没有带来太大的精度提升，反而由于分析的单元数增加了，扩大了整个系泊缆索求解的矩阵，导致计算耗时的增加。因此在采用Newmark-β迭代法进行数值计算时，在满足计算精度的要求前提下尽量减少单元数，对于提高计算效率有很大的帮助。同时可以看到，提高时间步长也减小了计算耗时，这一现象与丁佐鹏等所述情况有所不同，主要是因为此处控制精度取的是10-2，比较容易达到要求，增加时间步长对迭代次数的影响并不是很大。而丁佐鹏等在进行缆索动力分析计算效率研究时，由于针对系泊缆索动力分析迭代控制精度比较高，不容易达到，增加时间步长导致迭代次数明显增加，导致大时间步长反而比小步长更耗时。

2.2 控制精度的影响

将缆索划分为20个单元，时间步长取0.2 s，控制精度分别取10-2、10-3、10-4、10-5、10-6进行数值模拟，计算得到缆索顶端张力时历曲线，如图6所示。数值模拟时间取3 000 s，将各控制精度下每个时间步的平均迭代次数及计算耗时进行对比，见表3。

图6 各控制精度下顶端张力时历曲线对比

表3 各控制精度下每个时间步平均迭代次数及计算耗时

从图6中可以看出，各控制精度下顶端张力的时历曲线基本一致，表3显示每个时间步的平均迭代次数随着控制精度的提高而迅速增加，迭代次数的增加必然导致计算耗时的增加。因此，一味地提高控制精度并没有太大的意义，反而会增加耗时，得不偿失。在计算精度得到保证的前提下，降低控制精度有利于减少每个时间步的迭代次数，从而减小计算耗时，提高计算效率。

2.3 每个时间步迭代次数的影响

图7 各迭代次数下位移的精度 图8 各迭代次数下张力的精度

图9 每个时间步迭代次数为3次和20次时 顶端张力时历曲线(1 000 s)

图10 每个时间步不同迭代次数计算耗时

从图7和图8可以看出：各迭代次数下位移的精度远高于张力的精度，可以推断控制迭代精度时其实主要控制的是张力的迭代精度。但无论是位移还是张力，随着每个时间步迭代次数的增加，计算的精度均在不断提高。在迭代次数较少的时候，精度提高的速度很快，但是当迭代次数增多之后，精度的提高越来越慢，迭代次数达到10次之后精度的提高已不很明显。即使每个时间步的迭代次数仅为3次，张力的精度相较于张力的幅值已经很高了，且模拟时长达到1 000 s计算的结果如图9所示，也能保持良好的稳定性，并与迭代20次的结果吻合得很好。

从图10可以看出：计算耗时随着每个时间步的迭代次数的增加，基本呈线性增加趋势，因此控制每个时间步的迭代次数是提高计算效率的一种有效方法。

图11 并行计算

图12 并行计算耗时

3 并行计算

系泊浮体的系泊系统通常由多根系泊线组成。在对系泊浮体动力耦合分析开展研究的“弱耦合”方法中，在系泊缆索求解时各根系泊缆索动力分析是独立进行的，允许对系泊缆索并行求解提高计算效率，如图11所示。本文应用Intel Fortran自带的并行运算方法OpenMP实现多根系泊缆索并行计算，该方法相对简单，容易实现。

为研究并行计算的效率，取12根系泊线，简便起见，其参数均与第2节中的一致，且均做与上节中同样的强迫运动。将泊线划分为20个单元，时间步长取0.2 s，模拟时长取500 s。运行用计算机参数为Intel Core i5-4460，四核。并行计算时分别采用2线程、3线程和4线程进行计算，将并行计算与串行计算(即线程数为1)的耗时进行对比，如图12所示。

从图12可以看出：随着线程数的增加并行计算的耗时显著减小，说明将并行计算运用于系泊缆索的动力分析是十分有效的。

4 结 语

本文在系泊缆索动力分析细长杆理论的基础上，建立系泊缆索有限元模型，采用Newmark-β迭代法求解，编制动力分析计算程序，并利用商业软件OrcaFlex验证程序的正确性。以此为前提，本文着重讨论了各要素对计算效率的影响。其中系泊线单元数的划分对计算耗时的影响比较敏感，随着单元数的增加，计算耗时呈非线性增加趋势，适当增加时间步长有利于减小计算耗时。随着控制精度的提高，每个时间步的迭代次数增加速度越来越快，一味提高控制精度对计算结果精度的提高并没有显著的作用，反而造成计算耗时的大幅增加。控制每个时间步的迭代次数，能有效减小计算耗时，而且计算结果的精度也并未有太大损失。对于实际应用中深水悬链式系泊，系泊缆索的长度常在千米以上，根据前文算例分析及笔者经验，动力分析的时间步长可取0.2 s，单元数可取20，每个时间步的迭代次数可取3，动力分析结果精度能够满足实际工程的要求。另外，针对系泊浮体动力耦合分析的“弱耦合”方法，在对多根系泊线进行动力分析时可采用并行计算。本文针对系泊线动力分析开展的并行计算表明：随着线程数的增加，计算耗时呈线性减少趋势，采用并行计算能够有效降低系泊缆索动力求解的耗时。下一步将集成系泊缆索与浮体动力耦合分析程序，在此基础上开展系泊缆索并行计算下系泊浮体动力耦合算法计算效率研究，以便进一步提高计算效率。