陆争光，夏 志，阮 哲，黄冬云，王祺来

(中海油研究总院有限责任公司，北京 100028)

0 引言

近30年内，中国海域内已投运315条海底管道累计发生51起事故[1]，主要包括管道变形、穿孔小漏、裂缝大漏和断裂4大类，其中，在海洋平台吊机作业和补给船装卸过程中的落物风险是导致水下设施(海底管道、井口、采油树以及平台结构等)变形损伤的主要原因之一。为有效降低海洋平台落物风险、指导海洋平台方位设计，准确并快速预测海洋平台落物轨迹、落点和触底冲击能量的方法和技术显得至关重要。但目前，行业内依旧参考2017年DNVGL-RP-F107推荐性规范中的角度偏差经验值方法，经验化、平面化、过于简化、单一化等特点导致其计算能力有限(触底冲击能量)或计算结果(轨迹和落点)与现场实验结果存在较大误差[2-4]，在保证方法可行、结果可信的前提下，科学化、立体化、完整化、多因素概化的落物轨迹、落点分析方法和软件或工具有待攻克，以科学完善DNVGL-RP-F107推荐性规范。

本质上，落物水下运动属于典型的综合性工程学科问题，涉及到非线性动力学、流体动力学、运动理论和概率与统计方法等基础学科理论，国外早在20世纪80年代开始了大量研究，主要集中在模型化现场实验分析[5～11]、2D运动理论模型搭建与应用[5,6,12]、3D运动理论模型搭建与校验改正[8,12,14-15]等方面，相关应用主要体现在潜布无动力水雷运动受力、攻击型水雷下潜弹道预测[9-10,13]、海洋平台落物运动危险区和海底管线损伤分析[2,11,15]、海上失事目标搜救应急预报[16]等多个重点领域。其中，海洋平台圆柱体坠落运动理论模型研究的提升空间依旧很大，其研究应用停留在2D模型层面，方法可行性和模拟结果可信度均有待提高，研究趋势将主要集中在3D运动理论模型方面，从平面化向立体化、受力简化向完整化、多因素确定向多因素不确定、单一均质简单形状向多类非均质复杂形状不断改进或细化。

为此，笔者将以海洋平台圆柱体落物(管子、钻杆等)为研究对象，结合势流和粘性流等流体力学理论、弗劳德模型相似理论，考虑3维、6度受力与运动参数，建立具有普适性、可行性、可信性的的落物水下2D/3D轨迹预测模型，以科学化DNVGL-RP-F107落物轨迹预测及落点取值。

1 2D运动理论模型建立及验证

1.1 坐标系建立

2D右手直角坐标系与关键运动参数示意如图1所示，全局坐标系OXZ即地坐标系，局部坐标系oxz以落物质心为坐标圆心、落物轴向为x轴、落物径向为z轴。

图1 2D右手直角坐标系与关键运动参数示意Fig.1 Schematic diagram of 2D right hand coordinate system and key motion parameters

1.2 模型建立

为合理简化落物运动理论，假设圆柱体落物为均质、细长型刚性封闭体，其质心与几何中心重合，建立只考虑oxz平面坐标系内运动的2D理论模型，主运动方程如式(1)～(3)所示。

(1)

(2)

(3)

式中：m，ΔV分别为圆柱体质量和封闭体积；ρ为海水密度；β为2D初始倾斜角(x轴与X轴间的起始夹角)；U1，U3分别为圆柱体沿x轴、z轴的平移速度；W2为圆柱体绕y轴的旋转角速度；xt，mt分别为有效轴向边缘位置及其2D附加质量系数；m33，m55为不同维度的附加质量；I55为圆柱体绕y轴的惯性力矩；式(2)和(3)右侧中括号内公式项分别为势流理论中的附加水动力、力矩项，可由Lighthill理论公式推导计算。主运动方程中涉及到的其余参数项如下所示：

1)x轴粘性拖曳力Fdx

粘性拖曳力一般由摩擦拖曳力和压差拖曳力两部分组成，根据Gong Xiang等[14]相关研究，前者可通过平板边界层理论近似推算得到(与绕流雷诺数有关)，后者则由光滑圆柱体的绕流拖曳力公式计算，如式(4)所示。

(4)

2)z轴粘性拖曳力Fdz、y轴力矩Mdy

利用莫里森(Morison)方程计算Fdz，Mdy，如式(5)和式(6)所示。

(5)

(6)

式中：Cdz为z轴拖曳力系数；Uz(x)为z轴方向水流-圆柱体落物相对速度。

1.3 模型误差控制

按照细长体理论的基础假设，纵向平移速度应远大于其他方向的平移速度且基本保持不变。但是，实际落物运动过程并不满足该假设，为此模型运动方程引入了粘性力项，并以某一时间步长内起点纵向平移速度、平均横向平移速度计算位移递进过程。另外，为从算法上有效抑制模拟误差，模型节点平移速度递进可采用高精度单步算法R-K4(即四阶龙格库塔法)，一阶递进时间步长为0.01 s、四阶递进时间步长0.01 s/0.02 s。

1.4 模型验证

模型验证选取现场实验管子#1为对象，主要参数和现场边界如表1所示。

表1 现场实验管子#1参数与边界Table 1 Property of the Pipe #1 and experiment boundary

以Matlab 2016为程序设计平台、采用R-K4算法求解(未采用并行的单次计算时间未超过0.40 s)，模拟不同初始倾斜角、不同有效边缘位置系数、相同拖曳力系数下的运动轨迹及落点，验证分析结论如下：1)2D运动理论模型模拟轨迹、落点及下落时间等关键参数与现场实验结果基本吻合，关键参数可信；2)相对于Aanesland模型模拟结果[5]，修改后的2D运动理论模型在偏小初始倾斜角(β<60°)下的模拟落点和下落时间更趋近于实验结果，更合理；3)以θ0=60°为例，如图2所示，无论是Aanesland模型还是修改后的2D运动理论模型，均未能同时实现模拟运动轨迹和落点均在实验包络线内，存在一定程度不可避免的误差，其根本原因在于2D运动理论模型无法完全逼近实际运动过程，因此3D运动理论模型研究势在必行。

表2 现场实验与模拟运动轨迹关键参数对比Table 2 Comparison of key parameters between experiments and simulated trajectories using 2D model

注：XA，TA分别代表偏离起始落点垂直线最远处的水平位移及其坠落时间；XB，TB分别代表最终落点的水平位移及其坠落时间；Experiment-V即Aanesland现场实验数据；NS-V即Aanesland数值模型模拟结果；NS-L即本文基于改进的2D运动理论模型模拟结果。

图2 管子#1-现场实验与2D理论修正后的运动轨迹模拟对比(θ0=60°，Cdx=1.2，Cdz=1.0)Fig.2 Pipe #1- comparison between experiment and simulation trajectory based on modified 2D model at drop angle 60° with Cdx=1.2，Cdz=1.0

2 3D运动理论模型建立及验证

2.1 坐标系建立

DNVGL-RP-F107和早期的现场实验只注重其OXZ平面垂直、横向方向的运动，以横向X轴方向的平移距离来判定落点，但是在实际运动过程中，落物的轨迹是立体化、3D化的，应以在OXY平面内偏移初始落点垂直中心线的距离来判定落点。3D右手直角坐标与关键运动参数如图3所示。

图3 3D右手直角坐标系与关键运动参数示意Fig.3 Schematic diagram of 3D right hand coordinate system and key motion parameters

2.2 模型建立

基于改进的2D运动理论模型，在3D右手直角坐标系中，重点考虑落物轴向x轴、横向y轴、径向z轴3个方向上的平移速度和旋转角速度，引入粘性拖曳力系数经验公式、升力和相应力矩项、海流运动影响项，对式(1)～(3)扩展得到3D主运动方程如式(10)～(15)。

(7)

(8)

(9)

(10)

(11)

(12)

式中：φ，θ，ψ均为瞬时欧拉角，分别表示X轴与N轴(OXY与oxy平面的交线)间的夹角、Z轴与z轴间的夹角、x轴与N轴间的夹角；U2分别为圆柱体沿y轴的平移速度；W1，W3分别为圆柱体绕x轴、z轴的旋转角速度；cd为W1随时间变化的平均折损率；mt2，mt3分别为轴向边缘位置y轴、z轴附加质量系数；m22，m66分别为因沿y轴摇摆、绕z轴偏航引起的附加质量，m22=m33，m55=m66；I44，I66分别为圆柱体绕x轴、z轴的惯性力矩；Fdy，Mdz分别为y轴粘性拖曳力和力矩，利用莫里森(Morison)方程计算即可。主运动方程中涉及到的其余未定义参数项如下所示：

1)因绕x轴旋转W1产生的升力Fly/Flz，力矩Mly/Mlz

(13)

(14)

(15)

(16)

2)粘性拖曳力系数经验公式

Cdy，Cdz分别为y轴、z轴拖曳力系数，采用H.Rouse经验公式进行计算，如式(17)所示。

(17)

2.3 模型验证

模型验证以管子#1为模拟对象，同样以Matlab 2016为程序设计平台、对2D源程序语言进行扩展修改，采用R-K4算法或分段化升力参数项后采用一阶算法求解(单次计算时间均未超过0.40 s)，模拟不同倾斜角下的运动轨迹，以θ0=60°为例，基于3D模型模拟的平面轨迹如图4。

图4 3D模拟平面轨迹(θ0=60°，W10=0.05 rad/s，W20=W30=0 rad/s)Fig.4 Simulated OXZ plane and OXY plane trajectories at drop angle 60° with W10=0.05 rad/s and W20=W30=0 rad/s

通过对比分析不同初始欧拉角θ0下的模拟结果与现场实验结果，3D运动理论模型验证分析结论如下：

1)由OXZ平面轨迹可以看出，对标现场实验数据，相对于2D运动理论模型，在合理的Xt条件下，模拟轨迹和落点同时满足在实验包络线内，总体结果可信、与现场实验结果更加吻合；

2)由OXY平面轨迹可以看出，在考虑落物轴向旋转时，落物沿Y轴方向的位移是不可忽略的(介于0～0.88 m)，以在OXY平面内偏移初始落点垂直中心线的距离来判定落点更合理；

3)对于不同的初始欧拉角θ0，合理的Xt是不同的，整体呈现正相关关系，在相同初始欧拉角θ0下，Xt增大会使平面轨迹向轨迹曲线外开口侧偏移，但落点变化趋势不一。

3 落点规律分析

在模型验证可行、可信、更加合理的前提下，以在OXY平面内偏移初始落点垂直中心线的距离来判定落点，落点具体参数如表3所示，综合分析得出：

1)不同θ0下的落点参数差别较大，偏移距离R介于0～5.83 m、X轴向偏移距离XB介于0～5.81 m；随着的θ0增大，R和XB均呈现先增后降的趋势，于θ0=62°时达到最大值，运动轨迹整体呈现螺线曲率减小、向轨迹曲线外开口侧偏移趋势；

2)落点分布变化可大致分为以下3个阶段：①当θ0≤20°时，落点分布较为集中，为Y轴偏移稳定主导阶段，Y轴偏移距离YB对R贡献率较大(60%～100%)；②当20<θ0≤62°时，落点分布相对较为分散，为X轴偏移实力型快速主导阶段，YB对R贡献率迅速减少；③当62°<θ0≤90°时，落点分布相对更为分散，为X轴偏移机会型主导阶段，XB和YB均下降较快，但YB对R贡献率始终很小(小于1%)；

3)就目前掌握现场实验数据而言，当45≤θ0≤60°时，模拟落点分布参数与实验XB均值之间差异性较大。需要注意的是，每次现场实验的起始条件存在不确定性，因此XB并不是唯一值而是范围值；

4)与2017年DNV推荐规范经验性结果相比，落点R相对较大，但是与现场实验结果和模拟结果均更接近，这说明DNV推荐规范经验性结果存在一定的偏差或不合理性。究其根本原因，一是DNV推荐规范仅以横向X轴方向的平移距离来判定落点距离，忽略了Y轴方向的平移距离；二是DNV推荐规范未考虑海流运动的影响，而海流运动会对落点分布产生较为显著的影响。

表3 3D运动理论模拟、现场实验和DNV规范推荐的落点分布对比Table 3 Comparison of excursion distribution between experiments, simulated trajectories using 3D model and recommended results of DNV specification

图5 3D运动理论模型模拟不同θ0下的落点分布(W10=0.05 rad/s，W20=W30=0 rad/s)Fig.5 Excursion distribution using 3D model with drop angle from 0° to 90°, W10=0.05 rad/s，W20=W30=0 rad/s

4 结论及建议

4.1 结论

1)相对于现场实验数据，2D运动理论模型模拟关键参数可信，但并未实现全过程满足实验包络线条件；3D运动理论模型可行、整体优化结果可信、与现场实验结果更加吻合，以在OXY平面内偏移初始落点垂直中心线的距离来判定落点更合理。

2)落点分布变化可大致分为Y轴偏移稳定主导阶段、X轴偏移实力型快速主导阶段和X轴偏移机会型主导阶段3个阶段。随着的θ0增大，R和XB均呈现先增后降的趋势，运动轨迹整体呈现螺线曲率减小、向轨迹曲线外开口侧偏移趋势。

3)自主开发工具MREDP可实现不同长度直径比、不同海流运动下的圆柱体坠落全过程轨迹和触底冲击能量预测，具有工程应用性和推广性。

4.2 建议

1)海洋表面波浪的影响、非刚性体的弯曲特性可在3D运动理论模型中加以考虑、细化。

2)针对落物起始运动参数、入水撞击产生的运动参数变化等不确定性，可以借助Monte Carlo方法进行表征，模拟分析相应的落点分布、确定危险区。