水面拖曳系统龙须缆静态构型算法

孙洪波， 施朝健， 翁跃宗

(1. 上海海事大学 商船学院， 上海 201306； 2. 集美大学 航海学院， 福建 厦门 361021)

水面拖曳系统龙须缆静态构型算法

孙洪波1, 2， 施朝健1， 翁跃宗2

(1. 上海海事大学 商船学院， 上海 201306； 2. 集美大学 航海学院， 福建 厦门 361021)

为快速确定水面拖曳系统龙须缆的稳态运动参数，建立拖缆三维稳态平衡微分方程。根据拖缆两端的边界条件，利用二分法和龙哥库塔方法，对拖缆微分方程进行积分求解，确定水面拖曳系统中龙须缆的构型及稳态运动参数。以具体的拖曳实例为例进行数值仿真计算，分析龙须缆物理参数变化对系统稳态运动的影响。结果表明，运用这种算法，可以在拖曳系统初步设计阶段快速确定系统索具的参数，并能满足系统设计要求。

水路运输； 水面拖曳； 龙须缆； 静态构型； 二分法

水面拖曳系统作为一种水上运输及失事船舶救助系统，应用于诸多领域。其中的一种龙须缆拖曳方式，由于能够有效抑制整个拖曳系统在拖航过程中的偏荡，应用更加广泛。特别是在大型远洋船舶拖航时，采用的多是龙须缆。

龙须缆拖曳系统坐标系见图1，其中龙须缆PLOC和PROC的长度由被拖物拖力点距离确定，其夹角α应在35°～60°。沿海及遮蔽航区可用龙须缆；导缆孔磨擦区域必须使用龙须链，不得使用钢缆[1]；主拖缆OCPT为钢缆。

拖曳系统的静态构型指的是拖曳系统匀速直线运动时的运动姿态，在该系统的计算中具有十分重要的作用，不仅可以在系统设计初期直接帮助确定系统的物理参数，而且可以作为拖曳系统动态运算初始值，以加速运算收敛。因此，提出的龙须缆拖曳系统的静态构型计算方法不仅对龙须缆的设置有现实指导意义，而且可为系统的动态分析提供计算基础。

图1 龙须缆拖曳系统坐标系

拖缆水动力模型建模方法主要包括: 集中质量法、有限元法、有限差分法和直接积分法等。WALTON等[2]采用集中质量法研究了海军武器试验中水下锚链的二维运动响应，并给出了详细的求解算法和公式。中岛俊夫等[3]建立了船舶系留锚链的三维集中质量法模型，并给出了详细的算例。ABLOW等[4]在建立拖缆运动控制微分方程的基础上, 采用有限差分法求解其三维动态运动。虽然有限元法和直接积分法更适于分析拖缆静力，但计算繁琐。[5]国内的王飞、朱军等[6-8]在拖缆的建模与仿真方面，分别对集中质量法和有限差分法进行了详细的总结和研究。上述研究多集中于单缆及水下拖曳系统，关于水面拖曳特别是龙须缆拖曳的研究比较少。

基于以上分析，首先对拖缆进行空间离散；然后建立拖缆微元在流体中稳定运动时，在重力、浮力和流体阻力作用下的平衡方程；最后将拖缆的边界条件转换为一组初始值，采用二分法进行积分求解。

1 拖缆三维稳态运动数学模型

1.1拖缆系统坐标系转换关系

图1中，惯性坐标系Oxyz的3个坐标轴单位矢量为(i,j,k)，局部坐标系Otnb的3个坐标轴单位矢量为(t,n,b)，2个坐标系均为右手坐标系,其中：t为缆的切线方向；n，b为缆的2个法线方向。根据欧拉角的定义，可得2个坐标系的变换关系为

(i,j,k)=BV(t,n,b)

(1)

(2)

式(2)中:θ为水平面方位角；φ为垂向仰角。

1.2拖缆微元受力平衡方程式

设拖缆细长、有柔性、圆柱形状，受重力、水动力和惯性力作用，由拖船拖动前进，不考虑所受扭矩。拖缆微元Δs在t时刻所受张力、惯性力、水中重力和流体力分别为T，B，W和F。s为未拉伸的拖缆长度，S为拉伸后的拖缆长度。考虑一般拖缆材料的拉伸特性，有S′=∂S/∂s=1+eT，其中：e=1/EA，E为拖缆弹性模量，A为拖缆横截面积。根据牛顿第二定律，拖缆微元的动力平衡方程为

(3)

在计算拖缆稳态构型时，惯性力B=0，因此将式(3)与重力、浮力及阻力代入平衡方程中，并在拖缆局部坐标系下沿各个坐标轴方向展开，所得平衡

方程可写为标量形式：

(4)

式(4)中：ρ为海水密度；g为重力加速度；A为拖缆横截面积；m为每米拖缆质量；w=(m-ρA)g为水中每米缆重量；d=(1+eT)1/2d0为拉伸后缆的直径，d0为未拉伸缆的直径；Ct，Cn和Cb为拖缆的水动力系数；ut，un，ub分别为拖缆微元在3个局部坐标轴方向上的速度。

拖缆在惯性系下的坐标见式(5)。

(5)

2 边界条件及稳态初始值

为简化计算，设定拖船和被拖船的初始航向为零，这样主拖缆部分可简化到二维空间来考虑，而龙须缆的两根支缆存在以下关系：

(6)

式(6)中：下标L,R分别为左支缆和右支缆。

2.1不连续点的处理

经过上述简化，不连续点(即OC点)处的约束方程可根据经典力学平衡关系得到。

(7)

联立式(6)可得

(8)

式(8)中：下标M为主拖缆；下标MC,LC,RC分别为主拖缆连接点、左支缆连接点和右支缆连接点。

2.2端点的处理

首先考虑被拖船端，当给定了稳定拖航速度V后，被拖船端点PL和PR所受的水平拉力FL=FR=0.5XuuV2，Xuu为被拖船纵向运动的总阻力系数。根据拖缆两端的空间位置以及两根支缆的对称关系，拖缆的边界条件可确定为

(9)

式(9)中：下标PL,PR,PT分别为左支缆拖带点、右支缆拖带点和主拖缆拖带点；sC和B分别为支缆长度、两支缆拖带点间的距离；Z为水中浸深，如龙须缆两端都在水面处可取为0。

根据实际情况，可确定φPL,φPR和θPL,θPR的取值范围为：φPL,φPR∈(-π/2,π/2)；θPL∈[0,π/2)；θPR∈(-π/2,0]。

式(9)可采用二分法迭代求解，直至满足式(10)要求。其中，εz和εy分别为龙须缆两端位置点垂向误差和两支缆连接点位置误差，计算时均取0.001。

(ZPL-ZPT)|lt;εz

(10)

3 稳态运动计算仿真

3.1计算方法及步骤

通过联立式(6)和式(9)，可将式(9)简化为只含有2个变量的方程组形式：

(11)

式(11)可以采用二分法迭代求解。虽然二分法在扩展后可以解决二维、三维甚至更高维的问题[7]，但通过对龙须缆拖带的实际情况进行分析，可以将式(11)变为一维问题进行求解。由于左右两根龙须缆对称，仅需以其中一根龙须缆为例进行求解。以左侧龙须缆为例，根据φP和θP的取值范围，具体的计算求解步骤如下。

运用该方法迭代求解龙须缆静态构型的计算过程见图2，从中可以看出该方法每次循环的计算结果和龙须缆的收敛过程。需指出，φP的计算值在计算过程中是临时的或无用的。如果在计算过程中动态地设置φP的取值范围(即逐渐减小φP的取值范围)，还可以提高计算效率。

3.2计算仿真

为验证该算法的可行性与精确性，以直径为70 mm的拖缆为例进行计算仿真。拖缆主尺度见表1。

表1 拖缆参数

以300 m长的拖缆为例，进行不同拖力下的计算仿真，并与实验结果[9]进行比对，以验证算法的精确性，结果见图3。

通过与实验数据比较可知，本算法不仅可行，而且在精度上完全满足实际工程计算要求。

图3 300 m拖缆构型及悬垂深度

为进一步验证本算法在龙须缆构型计算中的可行性与精确性，以海洋石油支持船“海洋石油681”号拖船以225 t系柱拉力海上拖带“海洋石油981”号半潜平台为例进行计算仿真。表2为海上拖航的索具配置，图4为拖船和平台在拖航作业工况下静水中的拖航阻力曲线。[9]平台龙须缆悬挂点宽度B=56.56 m。

表2 海上拖航索具配置

图4 拖航阻力曲线

图5和图6为主拖缆长1 200 m，龙须缆长60 m，在最大拖航速度下的龙须缆构型及沿缆长的张力分布情况。

图5 拖缆构型

图6 沿缆长的张力分布

根据拖船与平台的阻力曲线，可以反推出此时的拖航速度为7 kn。但分析图5和图6可知，虽然此时满足拖缆安全工作负荷要求，但拖缆的悬垂深度并不能满足《海上拖航技术要求》中“远洋拖带移动式平台的悬垂值≥15 m”的规定[1]。其原因是在计算拖曳系统拖航阻力时，既考虑了拖船与被拖平台的拖航阻力，又计算了拖缆的拖带阻力。因此，计算方法不仅可行，而且更加精确。

图7和图8为不同拖航速度下的龙须缆构型及沿缆长的张力分布情况。

图7 拖缆构型

图8 沿缆长的张力分布

图9和图10为3 kn拖航速度、不同主拖缆长度下的龙须缆构型及沿缆长的张力分布情况。

图9 拖缆构型

图10 沿缆长的张力分布

4 结 语

在建立拖缆三维稳态平衡微分方程的基础上，利用二分法和龙哥库塔方法对拖缆微分方程进行了积分求解，最终确定了龙须缆拖曳系统稳态拖曳下的拖缆构型及拖缆上的张力分布，通过实际算例验证了算法的可行性与精确性。此外，还进行了不同拖航速度及拖缆长度下的计算仿真，结果充分表明

了拖曳系统参数对龙须缆拖曳系统的影响。因此，这种算法不仅对龙须缆拖曳系统设计初期确定缆索参数具有较强指导意义，而且由其得到的拖曳系统的稳态构型可为该系统的动态分析提供计算基础(即可加速龙须缆动态计算仿真的收敛性)。

[1] 中华人民共和国交通部.JT/T 214—1995,海上拖航技术要求[S].

[2] WALTON T S, POLACHECH H. Calculation of Transient Motion of Submerged Cables[J]. Mathematics of Computation, 1960,14:27-46.

[3] 中嶋俊夫，元良誠三，藤野正隆. 質点系モデルによる係留ラインの3次元動的解析法[C]. 东京：日本造船学会論文集, 1983.

[4] ABLOW C M, SCHECHTER S. Numerical Simulation of Undersea Cable Dynamics[J]. Ocean Engineering, 1983, 10(6):443-457.

[5] VANDANA R K. Finite Element Analysis of Under Water Towed Cables[C]. Proceedings of International Conference on Energy and Environment, 2013.

[6] 王飞, 黄国梁，邓德衡. 水下拖曳系统的稳态运动分析与设计[J]. 上海交通大学学报, 2008, 42(4):679-684.

[7] 王飞. 海洋勘探拖曳系统运动仿真与控制技术研究 [D].上海：上海交通大学，2007.

[8] 朱军，熊鹰，王志国. 拖缆系统直线定常运动仿真计算[J]. 海军工程大学学报, 2001,13(2):17-20.

[9] 林文锦, 蔡钿, 邵哲平. 海洋石油支持船拖航作业技术及应用指导书[M].北京:人民交通出版社,2011：10.

SteadyStateSolutionofTowedSurfaceSystemwithBridleConfiguration

SUNHongbo1,2,SHIChaojian1,WENGYuezong2

(1. Merchant Marine College, Shanghai Maritime University, Shanghai 201306, China; 2. Navigation College, Jimei University, Xiamen 361021, China)

In order to quickly determine the steady state motion parameters of the surface towed system with bridle configuration, the towing steady-state equilibrium differential equation is established. The differential equations are solved under the boundary conditions at the ends of the towing lines by both bisection method and 4th order Runge-Kutta method, to determine the towing bridle configuration and the steady state motion parameters. The numerical simulation was carried out and the effect of physical parameters of the towing bridle on the steady state motion of the system is analyzed. The simulation results prove that this algorithm can be used to determine the parameters of the systems quickly for preliminary design.

waterway transportation; surface towing; towing bridle; static configuration; bisection method

2014-08-12

国家自然科学基金(51109090)；李尚大基金资助项目(ZC2010011)

孙洪波(1977—)，男，黑龙江尚志人，讲师，博士生，研究方向为船舶运动建模及仿真。E-mail: sunhongbo1977@126.com

1000-4653(2014)04-0059-05

U661.3

A