三维缆索形态计算平台的开发与应用

,,

(华中科技大学 船舶与海洋工程学院,武汉 430074)

1 缆索三维流场理论模型

基于以下假设：缆索不可深伸长；缆索只能承受拉力；不能承受弯矩及压力；作用在缆索上的流体动力可以分解成切线方向的分力和法线方向的分力。设u，v，w为缆索上任意点P上的直角坐标单位矢量;u为P点处缆索方向。同样设i，j和k为大地坐标系；dsu为P点处的长为无穷小的缆索元。φ为水平面与矢量u之间的夹角；φ为u在水平面的投影与矢量i的夹角，则

dx=dsu·i=dscosφcosφ

(1)

dy=dsu·j=dscosφsinφ

(2)

dz=dsu·k=dssinφ

(3)

设三维速度场(缆索相对于海水的速度)为V=Vxi+Vyj+Vzk，阻力的切向分力F和法向分力G、H的表达式就为

F=Fdsu=ρCdtπd(V·u)2dsu/2

(4)

G=Gdsv=ρCdnd(V·u)2dsv/2

(5)

H=Hdsw=ρCdnd(V·w)2dsw/2

(6)

式中：

V·u=Vxcosφcosφ+Vycosφsinφ+Vzsinφ

(7)

V·v=-Vxsinφ+Vycosφ

(8)

V·w=-Vxsinφcosφ-Vysinφsinφ+Vzcosφ

(9)

Cdt——缆索的切向阻力系数,Cdt=γCdn;

其中：Cdn——缆索的法向阻力系数。

缆索元上的作用力之和为零，平衡条件为

-Tu+(T+dT)(u+du)+Fdsu+Gdsv+
Hdsw-kpdsk=0

(10)

忽略高阶无穷小量，得

(dT+Fds)u+(Tcosφdφ+Gds)v+
(Tdφ+Hds)w-pdsk=0

(11)

在u方向上：

dT+Fds-psinφds=0

(12)

在v方向上：

Tcosφdφ+Gds=0

(13)

在w方向上：

Tdφ+Hds-Pdscosφ=0

(14)

综上可得缆索的平衡微分方程：

(15)

2 求解方法

上述微分方程组显然是一阶常微分方程组，根据特定的边界条件对这些方程积分，得到缆索的平衡方程和拉力沿缆索径迹的稳态分布。

本文采用变步长四阶五级龙格-库塔-芬尔格(Runge-Kutta-Fehlberg)方法[1-3]来求解。针对拖曳缆索的初值问题，即已知初始条件为

X0=[T0，φ0，φ0，x0,y0,z0]

将X0作为输入条件沿着S=[0,S0]以一定步长进行迭代计算，其中S0为已知缆索总长度或者已知工作深度，就可以得出终端参数为

X1=[T1，φ1，φ1，x1,y1,z1]

同时步长点上的Xi也可以得出，这样就可以得出整个拖曳缆索的形态。

3 GUI计算平台

基于上述理论模型，利用Matlab图形用户界面功能，简称为GUI(graphic user interface)，编制计算平台。

程序核心部分采用变步长四阶五级龙格-库塔-芬尔格方法，利用ODE模块编写计算模块，同时利用GUI的I/O交互功能，对于缆索的一般初值问题，输入初始参数即可求解。

平台有文本框输入初始参数，并显示输出结果，根据计算结果在图像窗口自动绘制缆索三维形态。同时，实际工程中主要是已知绳长求深度和已知深度求绳长两种问题，故平台设置了工况选择菜单，可以根据实际需要选择不同工况求解。

4 实例

某拖曳器参数见表1。

表1 实例参数

4.1 无海流情况

如图1，调用GUI界面进行计算，选择“已知绳长求深度”工况，输入初始参数，计算结果：缆索放出800 m时，潜水器的深度为304.7 m，此结果与原文献结果十分接近。

4.2 有海流情况

有海流的情况，实质上是缆索相对于海水的速度矢量V=Vxi+Vyj+Vzk发生了变化，故求解方法与无海流情况无异。

对于有海流的情况，定义θ为海流速度矢量与拖曳速度矢量的夹角。当海流速度大小一定时，定性分析，θ从0～180°的变化过程中，海流阻碍拖曳的速度分量逐渐变大，那么必然使得拖曳缆索的深度逐渐变小。

图1 GUI计算界面

假设存在与拖曳方向不一致的海流，其流速为0.5 kn；并且假定拖曳器受的阻力和净浮力不变。当夹角为0、30、60、90、120、150及180°时，调用GUI界面进行计算，可以得出各角度下的缆索空间形态空间分布，见图2，计算结果见表2。

图2 不同海流夹角下缆索空间分布

表2 缆索深度随海流夹角变化值

可以看出，在同一缆长和假定拖曳器受的阻力和净浮力不变时，海流的存在对模型计算的影响是存在的，当海流速度为定值时，随着夹角θ增大，缆索深度变小，也即意味着，在同一深度时，随着夹角增大，所需缆长将增加。这与上述定性分析相吻合。

同时，可以对比θ与(180°-θ)的情况，事实上，这两种情况，海流速度在拖曳方向的法向分量大小相同，故缆索的空间形态应在同一空间平面(曲面)内，只是由于在拖曳方向的分量方向相反，造成缆索的深度不同。

在不同角度下，缆索终端的拉力T也不同，故工程实际中可以根据计算结果选择合适的拖曳方向，以在不影响拖曳目标的情况下实现拉力最小，从而达到节能目的。

4.3 考虑轴向变形的情况

上述理论模型假设缆索在拉力的作用下不存在轴向变形，下面进一步探讨缆索在拉力的作用下伸长，存在轴向变形的情况

ds′=ds(1+ε)

(16)

(17)

(18)

式中：ε——应变；

ds′——拉伸后的微元长；

ds——拉伸前的微元长；

A——拉伸后缆索的横截面；

d0——拉伸前缆索微元段的直径；

d1——拉伸后缆索微元段的直径；

E——材料的弹性模量。

根据式(16)～(18)，可推导出应变ε的表达式为

(19)

于是考虑弹性后的缆索三维模型微分方程组为

(20)

根据上述方程组，对上例平台中的方程组做出相应修改，调用GUI界面计算，对于铜制(或黄铜，青铜)缆芯，E=70～130 GPa，此例取为70 GPa，计算结果见图3。

图3 考虑弹性影响下的缆索曲线

可以看出，考虑弹性和不考虑弹性的情况下，两条曲线的吻合程度非常好，这是由于对于缆芯为铜制的缆索，弹性模量E一般较大，故应变就很小，所以对于常用缆索材料，弹性对于缆索形态的影响可以忽略不计。

[1] RECKTENWALD G.数值分析和MATLAB实现与应用[M].伍卫国,万 群,张 辉,等.译.北京：机械工业出版社，2004.

[2] 薛定宇.高等应用数学问题的Matlab求解[M].北京：清华大学出版社,2008.

[3] 张圣君.浮标拖曳缆索形状及其力学计算[D].武汉：华中科技大学，2003.