基于非線(xiàn)性有限元法的二維水下拖纜研究

邵 校，劉祚秋

（中山大學(xué)應(yīng)用力學(xué)與工程系，廣東廣州 510275）

0 引言

拖曳系統(tǒng)通常由拖曳電纜、拖曳體及專(zhuān)用絞車(chē)組成,拖曳體內(nèi)可根據(jù)不同的用途搭載溫、鹽、壓、營(yíng)養(yǎng)鹽及聲納等各種海洋化學(xué)元素探測(cè)傳感器或聲、光等物理探測(cè)傳感器。拖曳船上的操作者可以通過(guò)一定的控制方式實(shí)施對(duì)拖曳體的軌跡與姿態(tài)控纜作為整個(gè)系統(tǒng)的重要組成部分，起著連接拖曳母船（潛艇）與拖體設(shè)備的作用，它具有遠(yuǎn)距離傳輸能源和信號(hào)的能力，并為水下設(shè)備提供前進(jìn)的動(dòng)力。

國(guó)內(nèi)外學(xué)者對(duì)此進(jìn)行了廣泛而深入的研究，目前描述拖纜的理論方法大致可分為[3-4]：解析法[5]、有限元法[6-7]、集中質(zhì)量法[8]、有限差分法[9]。解析法是指那種可以給出封閉解的理論模型,由于該模型只能局限于求解簡(jiǎn)單的靜態(tài)問(wèn)題,其使用范圍有限；集中質(zhì)量法是將纜索視為由一系列無(wú)質(zhì)量彈簧連接著離散的集中質(zhì)量點(diǎn)。顧懋祥等人[10]的計(jì)算經(jīng)驗(yàn)表明：在使用集中質(zhì)量法時(shí)，纜繩模型分段太多將引起數(shù)值發(fā)散；有限差分法和集中質(zhì)量法在應(yīng)用過(guò)程中受到一些限制，且在強(qiáng)幾何非線(xiàn)性和材料非線(xiàn)性問(wèn)題上計(jì)算不易控制；基于小變形線(xiàn)彈性理論的線(xiàn)性有限元法忽略了纜索的及和非線(xiàn)性和材料非線(xiàn)性，在實(shí)際纜索計(jì)算中存在誤差。

為了克服其他計(jì)算方法的不足，本文采用非線(xiàn)性有限元方法計(jì)算纜索幾何形態(tài)及受力特性，考慮了纜索的各種受力及纜索的彈性變形，可有效對(duì)拖曳模型進(jìn)行模擬。通過(guò)Matlab語(yǔ)言，設(shè)計(jì)非線(xiàn)性有限元算法程序，結(jié)合懸鏈線(xiàn)法能使計(jì)算結(jié)果較快收斂。

圖1 水下拖纜模型

1 拖纜模型計(jì)算

如圖1中所示，由于纜索運(yùn)動(dòng)機(jī)理非常復(fù)雜，故作如下假設(shè)：1）忽略波浪對(duì)拖體設(shè)備的作用；2）只考慮在材料的線(xiàn)彈性階段工作，即索始終處于彈性工作階段，符合虎克定律；3）索是一維體；4）索是理想柔性的，只受拉力而不能受壓和抗彎；5）大位移小應(yīng)變。

1.1 懸鏈線(xiàn)方程

懸鏈線(xiàn)方法對(duì)松弛的（有剩余索鏈躺在水底的）和張緊的（索鏈全部提起其下端拉力傾斜的）兩種狀態(tài)都適用。由于張緊狀態(tài)懸鏈線(xiàn)適用范圍更廣，因此本文著重討論張緊式懸鏈線(xiàn)。為便于公式推導(dǎo)，分別建立笛卡爾坐標(biāo)系及拉格朗日坐標(biāo)系。以纜索右下端為原點(diǎn)定義笛卡爾坐標(biāo)系。以纜索右下端為原點(diǎn)定義為纜索弧坐標(biāo)。

圖2 靜水中纜索受力示意圖

如圖2所示使用截?cái)喾ㄍ茖?dǎo)懸鏈線(xiàn)方程，S0、Sd分別表示自然無(wú)伸張狀態(tài)，動(dòng)態(tài)纜索的弧坐標(biāo)。區(qū)間[0,l]上靜水中纜索只受端部力與重力的作用如圖1所示，S0=l處拉力為其中，水平方向的分量Th=Th0-Rwx，豎直方向分量Tv=Tv0+(W-B)-Rwy；Th0和Tv0分別為S0=0處作用在纜索上的水平力和豎直力，對(duì)于張緊狀態(tài)懸鏈線(xiàn)有：Tv0＞0。

纜索的拉格朗日應(yīng)變?yōu)椋害?(dsd-ds0)/ds0，即：dsd=(1+ε)ds0，s0=l處纜索微段與x軸正方向的夾角為?，那么：

以纜索微段切線(xiàn)方向與法線(xiàn)方向?yàn)樽鴺?biāo)軸建立局部坐標(biāo)系，如圖3。為簡(jiǎn)化計(jì)算將海水與纜索之間的相對(duì)速度vr視做勻速且只有水平方向的相對(duì)速度。那么纜索微段與水流切向相對(duì)速度vrτ、纜索微段與水流法向相對(duì)速度vrn分別為：

圖3 纜索微段受力示意圖

根據(jù)水動(dòng)力學(xué)Morison公式，流體施于微段纜索上的作用力為（圖3中所示）：

式中：ρw為海水密度；D為纜索直徑；CDτ、CDn分別為纜索微段切向和法向阻力系數(shù)。

將求出的fτ、fn再轉(zhuǎn)化到整體坐標(biāo)系下，結(jié)合式（1）即可得到：

在區(qū)間[0,s]上，纜索在水中的凈重：

式中：ρc為纜索在空氣中的密度；A為纜索橫截面積；g為重力加速度。若不考慮材料非線(xiàn)性，則：式中：E為彈性模量，A為纜索橫截面積。將以上各式代入方程的積分形式便可得到纜索各點(diǎn)的橫縱坐標(biāo)：

在區(qū)間[0,l]上取定積分可得：

顯然，這種情況下纜索的控制方程為非線(xiàn)性，無(wú)法求得解析解。為簡(jiǎn)化計(jì)算，令vr=0（即，使纜索處于靜水中），積分后可得下式：

1.2 非線(xiàn)性有限元基本方程

由上述推導(dǎo)過(guò)程不難看出拖曳纜索表現(xiàn)出很強(qiáng)的幾何非線(xiàn)性，在外部荷載作用下纜索產(chǎn)生大位移但應(yīng)變?nèi)匀缓苄。鋺?yīng)力應(yīng)變關(guān)系仍處于彈性范圍內(nèi)[11]，其平衡方程無(wú)解析解，須采用數(shù)值方法進(jìn)行計(jì)算。

圖4 纜索單元構(gòu)形及變形構(gòu)形

根據(jù)虛位移原理采用完全拉格朗日表述，推導(dǎo)得到連續(xù)體增量運(yùn)動(dòng)平衡方程[12]：

式中：

在大位移小應(yīng)變的假定下，柔性索單元的應(yīng)變?chǔ)胖恍杩紤]格林應(yīng)變的第一項(xiàng)，即：

以上各式中：E為彈性模量；ρc為纜索密度；為時(shí)間t+Δt時(shí)刻的加速度；為第二類(lèi)Piola-Kirchhoff應(yīng)力；和分別是Green應(yīng)變張量增量的線(xiàn)性項(xiàng)和二次項(xiàng)；時(shí)間t+Δt位形外力所做的虛功；為t時(shí)刻增量位移；為t時(shí)刻坐標(biāo)。

本文采用二節(jié)點(diǎn)索單元離散纜索，依照等參變換的概念，可直接得到單元節(jié)點(diǎn)坐標(biāo)向量及單元節(jié)點(diǎn)位移向量，各式有如下關(guān)系：

式中：[N]為插值形函數(shù)矩陣。忽略慣性項(xiàng)推導(dǎo)出拖曳纜索平衡方程：

在應(yīng)用有限元方法時(shí)需將作用在纜索上的水流力等效地作用在單元節(jié)點(diǎn)上。假定水流力所做的功為：

如圖3及所示，將式(3)代入式(25)可得：

1.3 非線(xiàn)性有限元方程求解方法

解決非線(xiàn)性有限元問(wèn)題常常使到

Newton-Raphson迭代法（N-R法）或修正的

Newton-Raphson迭代法（mN-R法）[13]，為了更好地得到收斂結(jié)果，這里使用N-R法。使用迭代法求解非線(xiàn)性問(wèn)題需要給定初始節(jié)點(diǎn)坐標(biāo)，若采用偏離真實(shí)值的較遠(yuǎn)的時(shí)極容易導(dǎo)致迭代，發(fā)散或漂移。使用上述的懸鏈線(xiàn)方程即可求得較好。計(jì)算的基本步驟為：1）假定纜索初始張力Th0、Tv0，使用本文中的懸鏈線(xiàn)法求得纜索上節(jié)點(diǎn)坐標(biāo)向量并假定纜索初始位移為零；2）將代入非線(xiàn)性有限元解法中進(jìn)行N-R迭代，直至收斂。

2 數(shù)值算例

為了驗(yàn)證本方法對(duì)拖曳纜索研究的適用性。以拖船端為拖纜起始端，分別取航速為3kn、5kn、7kn、9kn、11kn、13kn、15kn、17kn和19kn來(lái)模擬拖纜形態(tài)，計(jì)算纜索張力并作出纜索收斂圖。拖纜各參數(shù)如下表：

表1 拖纜及拖曳體參數(shù)

以(xi,yi)為坐標(biāo)做收斂圖（見(jiàn)圖5）觀(guān)察其收斂性。ne為單元個(gè)數(shù)，其中：

式中[U]為節(jié)點(diǎn)位移增量向量；[Ψ]為節(jié)點(diǎn)不平衡力向量。

圖5 各航速收斂圖

由圖5可以看出以(xi,yi)為坐標(biāo)的點(diǎn)隨著迭代步數(shù)的增長(zhǎng)都落在在了收斂區(qū)（虛線(xiàn)矩形框內(nèi)），說(shuō)明本文設(shè)計(jì)的方法收斂性較好。隨后將各航速下拖纜形態(tài)及受力情況做了對(duì)比（如圖6和圖7）。

圖6 不同航速纜索形態(tài)圖

圖7 不同航速纜索張力圖

由圖6和圖7可知雖然拖纜的形態(tài)及受力情況在不同航速條件下很大差別，但也呈現(xiàn)出一定的連續(xù)性。

3 結(jié)語(yǔ)

1）近年來(lái)非線(xiàn)性有限元方法常用于海洋系泊纜索的理論研究，本文將非線(xiàn)性有限元法配合懸鏈線(xiàn)法應(yīng)用于拖曳纜索的研究達(dá)到了較為理想的效果。以懸鏈線(xiàn)法的結(jié)果作為Newton-Raphson迭代法的初始節(jié)點(diǎn)坐標(biāo)能達(dá)到較好的收斂效果。

2）航速會(huì)影響拖曳體水平力與豎向力的比值進(jìn)而影響纜形，忽略拖纜振動(dòng)情況下，在拖船勻速行駛時(shí)拖纜呈現(xiàn)相對(duì)固定形態(tài)。

3）隨航速的遞增纜索張力越呈現(xiàn)非線(xiàn)性增長(zhǎng)，增長(zhǎng)越緩慢，在航速較大時(shí)，航速對(duì)拖纜張力影響不大。

[1]吳家鳴, 葉家瑋, 李寧.拖曳式多參數(shù)剖面測(cè)量系統(tǒng)水動(dòng)力與控制性能研究述評(píng)[J].海洋工程,2004.22(1):111-120.

[2]王飛.海洋勘探拖曳系統(tǒng)運(yùn)動(dòng)仿真與控制技術(shù)研究[D].上海: 上海交通大學(xué),2007.

[3]Vineetk.Srivastava.Analyzing Parabolic Profile Path for Underwater Towed-Cable[J].Journal of Marine Science and application, 2014, 02: 185-192.

[4]張攀.拖曳系統(tǒng)仿真計(jì)算[D].武漢:武漢理工大學(xué),2005.

[5]Irvine H.M, Caughey T.K.The linear theory of free vibrations of a suspended cable[A].Proceedings of royal society of London, series A, 1974, 341: 299 -315.

[6]Leonard J W, Recher W W.Nonlinear dynamics of cables with low initial tension[J].Journal of the Engineering Mechanics Division, ASCE.1972, 98 (EM 2): 293 -309.

[7]Ma David, Leonard John, Chu Kuang-Han.Slack-elasto-plastic dynamics of cable systems[J].Journal of the EngineeringMechanics Division ASCE.1979, 105 (EM2):207 -222.

[8]Walton T.S, Polacheck H.Calculation of transient motion of submerged cables.Mathematics of Computation,1960, 14, 27-46.

[9]Ablow C M, Schechter S.Numerical simulation of undersea cable dynamics.Ocean Engineering, 1983,10(6),443-457.

[10]顧懋祥,邵建南.拖曳系統(tǒng)動(dòng)力學(xué)計(jì)算[J].艦船性能研究, 1988, (3): 31-43.

[11]楊敏冬, 滕斌, 勾瑩，姜?jiǎng)俪?海洋系泊纜索非線(xiàn)性有限元靜力分析[J].海洋工程, 2009,27(2):14-20.

[12]Bathe K J.Finite element procedure[M].Englewood Cliffs: Prentice-hall, 1996.

[13]王興剛, 深海浮式結(jié)構(gòu)物與其系泊纜索的耦合動(dòng)力分析[D].大連: 大連理工大學(xué), 2011.