半潜式平台系泊运动数学模型及其仿真视景构建

张磊， 贾鹏， 苏峰， 王恒， 姚广旭

（1.海洋石油工程股份有限公司，天津 300451；2.哈尔滨工程大学机电工程学院，哈尔滨 150001）

0 引言

随着近海油气的开采日渐减少，对深海油气资源的开发已成必然的趋势，海洋平台相关技术成为海洋工程界的研究热点[1]，半潜式采油平台的应用越来越广泛。因此使用数学方法对半潜式平台锚泊运动进行研究，对于平台的安全性评估、保证平台正常运行与生产，以及对远洋深海工程锚泊定位技术的发展具有重大意义。

现阶段随着越来越多的深海浮式生产系统的服役，对各种浮式生产系统的运动方面的仿真也越来越多。对于油气田浮式生产系统的仿真，基本上分为对各种类型的半潜式平台的锚泊运动的仿真以及对与其相连的锚泊锁链的仿真。其中，由多点系泊锚泊锁链组成的锚泊系统，在环境的作用下为平台提供了巨大的回复力，使平台固定在海面的指定范围内[2]，半潜式平台在海洋环境中，在环境风、浪、流载荷的作用以及锚泊系统的回复力作用下，能够产生6自由度方向上的运动，基于以上原理，平台和锚泊系统组成了一个具有复杂运动的耦合运动系统[3-6]。

目前，对于半潜式锚泊系统的仿真技术已经相当成熟，C.M.Leech使用绳子模拟，在理论和实际方面分别对悬链线方法做了分析，能够分析出悬链线扭曲的状态，也能够对输水软管做出模拟；Joel Brown、Jean Claude、Latombel等基于绳子提出了一种计算方法，能够模拟出绳子在运动状态下打结的效果；Gunnar Teichelmann、Mike Schaub、Bernd Simeon等基于铁路电缆线的研究，在稳定受力平衡和动态振荡的效果下，研究了对悬链线仿真方法。由于半潜式平台结构相对复杂且种类较多，现阶段并没有普遍的数学模型和经验公式对其6自由度运动进行模拟，对其运动的模拟，更多的是采用有限元分析软件，基于有限元模型对其频率特性、运动特性等做出分析[7-9]。

为了模拟复杂海况中各种环境载荷对平台锚泊系统所产生的影响，本文根据所给定的平台尺寸构建合适的半潜式平台动力学模型，然后针对所得到的数学模型，研究系统锚泊运动特性，在Unity3d中开发虚拟海洋环境，使用3dsmax建立半潜式平台的三维模型，并将三维模型导入到Unity3d的虚拟场景中，实现对半潜式平台仿真实景的构建。

1 半潜式平台锚泊运动方程

在实际海洋油田环境中，半潜式平台与众多锚泊线连接组成耦合系统[10]。对半潜式平台的运动分析需建立2个方程：一为半潜式平台本身的6自由度运动方程，另一为锚泊线拉力方程。

1.1 半潜式平台运动方程

建立半潜式平台运动方程，将平台考虑为刚体，以平台自身重心为原点、以平台侧弦方向为X轴、以头部方向为Y轴、以垂直于甲板方向为Z轴建立惯性的空间笛卡尔坐标系。在惯性空间坐标系中应用牛顿第二运动定律建立平台运动方程式如下：

式中：Mb为平台质量矩阵；Cb为平台辐射阻尼矩阵；Kb为平台刚度矩阵；U为平台的位移矢量；为平台的速度矢量为平台的加速度矢量；Fb为平台受到的环境载荷。

1.2 锚泊线方程

对于锚泊系统的分析，需分析自然状态下的单根锚泊线对平台的拉力，然后基于悬链线方程建立单根锚泊线数学模型，悬链线方程是一根固定于两点间的自然悬垂状态下的曲线。在笛卡尔坐标系中，悬链线方程为：

方程三式分别代表在笛卡尔坐标系中，悬链线方程纵坐标y、长度L、拉力T随横坐标x的变化，其中a=H/（ρg）；H为悬链线横向拉力；ρ为悬链线在水中单位重量；θ为悬链线初始张角。

将悬链线初始张角θ近似取零，则可得到简化后的悬链线锚泊方程[11]：

2 环境载荷

对于海洋上的受约束低速运动结构，其受到的环境载荷包括环境风、浪和海流的作用。对于海洋波浪力载荷的确定，一般基于线性波理论，确定波浪力数学模型。风载荷和流载荷作用对平台微幅运动的效果影响并不明显，在实际计算中，一般将风载荷和流载荷考虑为线性作用[12]。

2.1 波浪干扰力数学模型建立

海洋中的波浪具有高度的随机性，海浪并没有固定的参数、形状等供人们研究。实际工程中，波浪一般被分为规则波和不规则波，规则波与不规则波又根据其对漂浮物的作用特性，被分为主体力与横漂力。本文分别建立规则波和不规则波的主体力及横漂力数学模型，来描述整体的波浪干扰力。

2.1.1 规则波浪干扰力数学模型

对于被视作刚体的海洋漂浮物，其受到的规则波主干扰力可被分解为作用在其6个自由度上对其结构体积的三重积分式：

式中：k为2π内波的个数，k=2π/λ；λ为波浪平均波长；a为波浪运动幅度；c为波浪的传播速度，其与频率有关系，c=ω/k；ψ为平台艏摇角。

以上对规则波浪主干扰力的计算公式，是基于对半潜式平台浸水体积的积分。在计算时，将需要计算的部分考虑为若干立方体的组合以简化计算。

2.1.2 规则波浪横漂力数学模型

在波浪力中，横漂力是作用在水平面方向上的作用力。实际计算时，通过已知条件波浪横漂力与波幅的平方成正比来简化计算[13]。规则波浪横漂力标准形式如下式：

式中：ρ为海水密度；L为半潜式平台长度；a为平均波幅；CXwD、CYwD、CNwD分别为纵向、横向以及波浪力矩飘移系数，χ为波向角。

其中，CXwD、CYwD、CNwD可通过下式计算：

2.1.3 不规则波主干扰力数学模型

对于不规则波的数值计算，一般采用统计学对其进行分析。从统计理论方面考虑，不规则波可以视为是对若干个参数不同的单位规则波的叠加，基于以上理论，海浪特征可用波能谱的概念进行描述。Jonswap波能谱函数为

式中：ωP为峰值频率，rad/s，参数β和σ一般取：

γ为谱峰提升因子，一般取平均值3.3。α取下值：

其中Hs为有义波高，一般取平均波高的1.6倍。基于波能谱函数，不规则波主干扰力数学模型如下式所示：

对于不规则波主干扰力的数学模型，也是6自由度方向上对不同频率的规则波在平台浸水体积的三重积分的叠加[4]，为简化计算，在波能谱函数中，仅取影响较为明显的规则波叠加到数学模型中。

2.1.4 不规则波横漂力数学模型

不规则波浪中波浪飘移力的计算，仍然考虑为对各种参数的规则波波浪飘移力的叠加[14]，不规则波二阶波浪干扰力如下式：

2.2 风干扰力数学模型

对于半潜式平台而言，受风部分主要是平台上层建筑以及水线面以上的平台浮桩，风力大小受到受风面积大小、受风点高度、受风面形状以及风速的影响。半潜式平台受到的简化水平风力如下式：

式中：Cw为半潜式平台风力形状系数，对于平台的表面结构取值为1；AT为平台横向受风面面积；AL为平台纵向受风面面积；Lpp为平台两立柱间距离；ρw为空气密度；VwR为平面上10 m处相对风速。

2.3 海流干扰力数学模型

海洋中海流流速一般随海水深度的增大而减小，对于半潜式平台，其整个水下部分基本处于浅水区，垂直方向流速变化不明显，故将流速做常数值处理。对于半潜式平台，海流力可通过如下公式计算：

式中：Cc为半潜式平台流力系数，对于平台的表面结构取值为1.5；Ap为平台横向浸水面面积；Aq为平台纵向浸水面面积；Lpp为平台两立柱间距离；ρc为海水密度，VcR为平均相对流速。

3 半潜式平台6自由度方程水动力系数分析

半潜式平台6自由度方程水动力参数包括质量矩阵系数Mb、辐射阻尼矩阵系数Cb与刚度矩阵系数Kb，其中质量矩阵包含半潜式平台主体结构质量与附加质量，附加质量在方程中代表与平台接触的水体速度变化时对平台的作用力；辐射阻尼矩阵为与平台接触的水体对平台的黏性阻力作用；刚度矩阵是由平台结构及锚泊系统对平台产生的回复力作用。

3.1 附加质量与辐射阻尼的计算

3.1.1 基于有限元分析的平台水动力参数计算

使用AQWA对半潜式平台的附加质量及辐射阻尼做出分析。对于海洋浮体的水动力分析，一般只考虑浮体的外壳受到的影响而忽略浮体内部结构。AQWA默认的挠射单元为处于水线面以下的部分，而水线面以上的部分是不需要做水动力分析的。如图1所示为平台主体的有限元模型。

经计算后，在AQWA中打开平台有限元模型如图2所示。

AQWA会生成平台的附加质量矩阵以及辐射阻尼矩阵中的全部参数，其中数量较大的为矩阵右对角线的6个参数，矩阵中其余参数全部视作0。

3.1.2 半潜式平台附加质量矩阵

图1 平台有限元模型

图2 平台AQWA模型

在平台的二阶运动微分方程中，Mb是由平台主体质量矩阵和附加质量矩阵组成的，其中附加质量矩阵为6×6的矩阵，在AQWA计算得到矩阵右斜对角线上的参数如图3所示。

所以平台的质量矩阵为：

其中，M为平台主体质量矩阵。

3.1.3 半潜式平台辐射阻尼矩阵

同理，平台的二阶运动微分方程中的阻尼矩阵Cb与Mb有同样的形式，在AQWA计算得矩阵右斜对角线上的参数如图4所示。

所以平台的阻尼矩阵为

3.2 半潜式平台刚度计算

半潜式平台的刚度受到锚泊系统以及平台自身结构的影响，其中锚泊系统提供平台在横荡、纵荡及艏摇方向上的刚度，平台自身结构提供平台在横摇、纵摇及垂荡方向的刚度。

3.2.1 锚泊系统回复力计算

平台的锚泊系统由11根锚链组成，如图5所示。求解11根锚泊线在水平方向的和作用，通过分析平台在水平面上受到的锚泊作用随平台水平面上位移的变化，求解锚泊系统在横荡、纵荡和艏摇方向上的刚度系数。

对于锚链在水平方向上拉力的计算：根据其在水平面的俯视长度和锚泊线长度方程，解得参量a；根据横向拉力公式H=ρga解得单根锚链横向拉力；再通过向量加法计算合力及合力矩，如图5所示。经计算，平台在水平面上指定位置受到的锚泊作用如表1所示。

图3 平台附加质量系数

图4 平台阻尼系数

图5 锚泊系统刚性计算

以上数据点对半潜式平台在横荡、纵荡和艏摇方向上受到锚泊作用进行线性拟合，得到图6～图8。

3.2.2 半潜式平台自身回复力计算

在半潜式平台的锚泊运动中，平台自身回复力出现在横摇、纵摇与升沉3个方向上，且主要受平台结构的影响。以下分别对3个方向的回复力作出分析计算：

1）升沉方向回复力计算。半潜式平台漂浮在海面上时，其在垂直方向上具有上浮和下沉的效果，平台的上浮和下沉均会引起平台排水体积的变化，也就会引起其在垂直方向上受力的变化。

图6 平台横荡方向上锚泊作用拟合图

由于浮筒在水线面上的运动面积是不发生变化的，如图9所示，所以对于平台在升沉方向的刚度可做简单的线性处理，如下所示：

Ka=ρgAσ。

其中Aσ为水线面的面积，式中的算术值代表平台在升沉方向上具有单位位移时，其受到浮力的变化。

表1 平台在横荡方向上受锚泊作用点数据

表2 平台在纵荡方向上受锚泊作用点数据

表3 平台在艏摇方向上受锚泊作用点数据

2）横摇与纵摇方向回复力计算。当平台具有横摇或纵摇运动时，其受到的重力和浮力的合作用会产生反方向的倾覆力矩，则平台在横摇和纵摇上的刚度，可用下式表示：

其中zb为平台重心与浮心间的垂直距离，上式代表平台具有单位角位移时，其在相应方向上倾覆力矩的变化，由于平台的摇动角度相对较小，所以忽略θ角的影响。

4 半潜式平台锚泊运动仿真视景构建

半潜式平台锚泊运动的仿真视景，一方面是支撑三维视景系统运行的硬件部分；另一方面为视景软件部分，包括虚拟海洋场景、三维模型以及视景的后台程序。

图7 平台纵荡方向上锚泊作用拟合图

图8 平台艏摇方向上锚泊作用拟合图

4.1 视景硬件系统

视景系统的硬件部分包括信号采集设备及视景显示设备两部分。信号采集设备用以探测海空的各类环境参数，并将其处理成计算机能够识别的信息，并在视景机端通过数学模型解算，将解算的结果以视景的形式在视景显示设备上显示出来。

4.2 半潜式平台视景构建

图9 平台自身回复力计算原理图

在Unity3d中开发虚拟海洋环境，使用3dsmax建立半潜式平台的三维模型，并将半潜式平台的三维模型的.FBX格式导入到Unity3d的虚拟场景中，组成半潜式平台的虚拟视景。

4.3 差分法解平台运动微分方程组

在视景系统中，使用VS2012作为开发环境，在VS2012中，使用C#语言开发脚本，编辑差分法解算平台的运动微分方程组，并将解算得到的平台位移通过脚本实时赋予视景中的平台，从而赋予平台相应的运动效果。

4.3.1 差分法解微分方程组

半潜式平台的运动微分方程是六元二阶微分方程。求解二阶微分方程最终目的是求解半潜式平台6自由度方向上的实时位移，在二阶微分方程中可看出，方程结果受到海况参数和环境参数影响。系统中，海况参数和环境参数通过通信接口实时获取，在具备这些参数的情况下，半潜式平台6个自由度方向上的位移也就是关于时间的函数。根据以上分析，将半潜式平台运动微分方程简化成下式：

其中，p（t）、q（t）、f（t）均为[a,b]上的连续函数，则根据二阶微分方程解的条件，可得微分方程存在唯一解。

二阶微分方程求解往往伴随着初值问题，方程的初值［u（a）,u′（a）］代表半潜式平台在某一方向上的位移与速度，如下式所示：

图10 视景硬件系统

差分法求解微分方程是基于连续方程离散化的原理，首先将区间a,［］b离散为N等份，则有下式：

图11 半潜式平台锚泊运动仿真视景

用ui表示u在点xi处的取值，对连续可导函数u（t），由泰勒级数展开式可得：

则将式（1）、式（4）与式（5）联合，省略局部截断误差，可得二阶微分方程三点差分格式式：

上式可改写为：

其中：ai=2/h-pi,bi=-4/h+2hqi,ci=2/h+pi,di=2hfi。

以上即为二阶微分方程差分离散形式。

由数值微分公式u′（a）≈u1-u0/h以及初值条件u′（a）=β可得：

则根据式（7）可依次求解各差分点值。

本文需求解的二阶微分方程如下：

则对于平台的单一自由度，其标准差分形式如下式：

式中的各参量中，仅有Fb为随时间变化，其余参量均与平台水动力系数有关，半潜式平台在一定运动频率范围内时，取平台水动力系数为常数进行迭代计算。

4.3.2 差分法解微分方程组脚本编辑

在视景系统中，使用C#语言编辑差分法解算平台的运动微分方程组。半潜式平台的6自由度微分方程经处理后，其最后能够简化成6个平台位移、平台受风向角、计算风速、波浪周期、波浪平均波长、波浪运动频率、波浪传播速度、洋流速度、平均浪高以及时间有关的6个二阶微分方程，通过通信接口或者系统数据中的环境参数和海况参数，并通过计算获得波浪运动频率，则平台的二阶运动微分方程化简为6个位移随时间变化的二阶微分方程。

在系统中，定义6个自定义函数，使用差分法分别对平台6个自由度方向的运动微分方程进行求解，并将解得的运动数据赋予平台。在单个自定义函数中，赋予平台运动的流程如图12所示。

图12 赋予平台运动的流程

图中的程序中，每隔0.5 s取得平台6自由度中单一方向上的位移量，并赋予平台驱动平台产生运动。其中子函数实现了对平台运动微分方程的解算，其流程如图13所示。

分别在上文的6个自定义函数中编辑图的程序，同时计算出平台6自由度的实时位移，图中设置迭代步长 为 0.5，也就是系统每间隔0.5 s输出一个位置点，根据精度的要求，可缩小迭代步长。

图13 平台运动微分方程的解算流程

5 半潜式平台锚泊运动仿真视景系统测试

使用Matlab对半潜式平台运动微分方程组加以解算，得到平台的位移随时间变化的关系；并在视景系统中记录平台位移随时间变化的关系，验证视景系统性能，并验证视景系统对平台运动的准确性。

5.1 运动微分方程组解算

取平台6自由度上的初始位移与初始速度均为0，使用表4的海况数据及表5的平台结构尺寸，分析平台在6自由度方向上的运动，如图14所示。

使用表4海况数据对平台6自由度方向上运动的解算结果如图14（a）～（f）。

本文对于微分方程的求解，赋予方程的初始值为[0,0]，其物理意义代表平台在单一自由度方向上的位移和速度均为0，在宏观上，其相当于在三维空间的（0,0,0）位置，赋予平台初始速度为0，让其在海水作用、环境载荷和锚泊作用下做进一步的运动。对以上各图中平台的运动结果进行分析，在海水作用、环境载荷以及锚泊作用下，由于三维空间中的[0,0,0]位置并不是平台受力平衡的位置，在平台的横荡和纵荡以及艏摇方向上，平台会在风力、海流力以及二阶波浪力的作用下产生偏移，偏移达到一定程度时，会因为锚泊系统的回复力作用，将平台拉回反方向，如此反复直到平台到达平衡位置，在运动图像中，其反应为平台的振荡运动，在垂荡、横摇和纵摇方向上，由于平台自身的回复力作用，其具备相同的振荡效果。当各自由度方向振荡结束后，在规则及不规则波浪干扰力的作用下，平台会在6自由度方向上产生具有一定规则性的小周期运动，这部分运动即是平台的主体运动。通过对以上各图的分析，在海况较为不稳定的情况下，平台在6自由度方向上的小周期振动也会比较强烈。

表4 计算用海况数据

表5 平台主要尺寸表 m

图14 平台历时运动效果图

图15 平台运动验证

为验证半潜式平台锚泊运动的准确性，需要将计算结果与实测数据进行比较。本文在计算数据与实测数据的比较上，忽略掉其运动的振荡环节，取平台运动规律性较强的平稳运动时段进行比较与分析。图15（a）～（f）为对平台运动计算数据与实测数据的比较。

通过以上的分析，应用数学模型计算的结果，与实际海洋中的大型漂浮物的锚泊运动规律是大体类似的。将平台在6自由度方向上运动的计算结果同实测数据进行比较，可发现数学模型计算结果同实测数据是基本吻合的，其中若干数据点的运动数据脱离数学模型的控制，经分析是由于海洋海况不稳定性造成的，如果不考虑这些不稳定性因素的影响，使用本文中的平台运动微分方程，是可以在很大程度上描述出平台的运动规律的。并能够应用于三维视景系统，描述平台的6自由度运动。

5.2 视景系统性能测试

在视景系统软件端，是采用差分法[15]，每隔0.5 s计算一个平台位置点，从而实现对平台位移的解算。其结果受差分间隔和视景机性能的影响。

图16 平台运动视景效果图

表6 初始时刻平台垂荡方向运动数据 m

表7 系统运行20 min后平台垂荡方向运动数据 m

表8 系统运行40 min后平台垂荡方向运动数据 m

表9 系统运行60 min后平台垂荡方向运动数据 m

表10 系统运行80 min后平台垂荡方向运动数据 m

表11 系统运行100 min后平台垂荡方向运动数据 m

表12 系统运行120 min后平台垂荡方向运动数据 m

在视景系统接口接入表4中海况信息，从而在系统中实现以差分法对平台运动微分方程组的解算。在系统中，平台运动效果可通过视景输出及数据输出表现出来，如图16所示。

表13 系统运行140 min后平台垂荡方向运动数据 m

表14 系统运行160 min后平台垂荡方向运动数据 m

表15 系统运行180 min后平台垂荡方向运动数据 m

如图16所示，当三维视景系统中具有海况信息时，系统会根据海况数据信息，通过数学模型解算出平台的6自由度方向上的运动规律。图中为根据数学模型解算出的平台运动规律在视景、图像和数据方面的表现。本文中监控的海况波浪平均周期为8.7 s，取图像的监控时长为60 s，监控时长可根据实际需要手动调节其大小。

由于数据量较大，本文取平台垂荡方向运动数据进行分析，并同第3节中对方程的计算结果进行比较，分析三维视景系统对平台运动描述的准确性，如表6～表15所示，每隔20 min取一组一周期上的运动数据，并对系统进行连续3 h的观测。

对表中平台的运动信息进行分析，由于系统运行的迟滞性，当系统长时间运行后，系统对平台运动数据的表述也会出现迟滞。在60 min以内，系统能够对平添运动信息进行准确表述，并进一步表明了平台运动微分组对平台运动描述的准确性。系统运行超过60 min后，系统内计算数据与理论数据的偏差会越来越大。

6 结论

建立了六元的二阶微分方程组，并分析了方程组中的质量矩阵、阻尼矩阵、刚度矩阵及环境载荷，并建立了半潜式平台的运动仿真视景，验证了方程组描述半潜式平台运动的合理性。

［参考文献］

[1] 孙涛,桂文彬,俞志刚.半潜式平台定位系泊控制试验系统设计与应用[J].船舶工程,2012,34(2)：84-86.

[2] 张维.深水平台锚泊系统锚索链受力分析及视景仿真[D].武汉：武汉理工大学,2014.

[3] XIE S.Design of Corrosion of Marine Oil and Gas Fields Underwater Manifold Pry[J].Guangdong Chemical Industry,2016,43(7)：158-159.

[4] 韦红术,张伟国,陈国明,et al.Position determining method and position determining device of Christmas tree of marine oiland-gas field producing platform：CN 103132991 A[P].2013-06-05.

[5] 王俊荣,谢彬.半潜式平台水动力性能及运动响应研究综述[J].中国造船,2009(a11)：255-261.

[6] 史琪琪,杨建民.半潜式平台运动及系泊系统特性研究[J].海洋工程,2010,28(4)：1-8.

[7] KOO B J,KIM M H.Hydrodynamic interactions and relative motions of two floating platforms with mooring lines in sideby-side offloading operation[J].Applied Ocean Research,2005,27(6)：292-310.

[8] YANG M,TENG B,NING D,et al.Coupled dynamic analysis for wave interaction with a truss spar and its mooring line/riser system in time domain[J].Ocean Engineering,2012,39(1)：72-87.

[9] 朱忠显,尹勇,申和龙.基于集中质量法的锚泊系统动力学仿真[J].中国航海,2015,38(4)：38-42.

[10]童波,杨建民,李欣.深水半潜式平台系泊系统动力特性研究[J].海洋工程,2009,27(1)：1-7.

[11] 王冬姣.索-链-浮子/沉子组合锚泊线的静力分析[J].中国海洋平台,2007,22(6)：16-20.

[12] 罗强.981半潜式钻井平台锚泊系统研究[D].荆州：长江大学,2012.

[13]韩鑫,蓝尉冰.二阶波浪载荷计算方法探究——以深海船式平台动力定位系统为例[J].钦州学院学报,2014,29(5)：1-3.

[14] 姜宗玉,崔锦,董刚,等.不规则波中动力定位半潜平台慢漂运动研究[J].中国海洋平台,2015,30(3)：89-94.

[15] 张文生.科学计算中的偏微分方程有限差分法[M].北京：高等教育出版社,2006.