張廣磊，張維競(jìng)，劉濤

(上海交通大學(xué) 海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

水下拖纜系統(tǒng)在海洋上有著廣泛的應(yīng)用:軍事上，海軍利用水面拖曳陣列聲吶系統(tǒng)來執(zhí)行勘探及反潛任務(wù);民用上，海上資源勘探拖纜系統(tǒng)是目前世界上進(jìn)行海洋油氣資源勘探的主要手段.拖曳陣列系統(tǒng)通常由拖船(或軍艦)、引導(dǎo)拖纜(重力纜)、零浮力拖纜陣列、水下拖體及各種纜載探測(cè)控制設(shè)備等組成.拖纜作為整個(gè)拖曳系統(tǒng)的重要組成部分，而我國在海洋拖纜系統(tǒng)方面的研究起步較晚，技術(shù)上明顯落后于美國等發(fā)達(dá)國家［1］，對(duì)于拖纜動(dòng)力學(xué)方程求解方法的研究相對(duì)較少.基于二階精度的有限差分Box method［2-3］是目前拖纜方程求解最普遍適用的方法，它把拖纜的控制方程在時(shí)間和空間上采用半網(wǎng)格點(diǎn)方法離散，自從Ablow和Schechter［4］首次采用它來求解拖纜的方程以來，由于它具有精度高和易實(shí)現(xiàn)的特點(diǎn)，相繼被Park、Ioannis、HUO Cunfeng和Mark等［5-8］很多研究人員采用，但該方法在時(shí)間上離散求解存在著穩(wěn)定性較差的缺點(diǎn)，Koh等［9］在其論文中證明了有限差分Box method時(shí)域上穩(wěn)定性差的缺點(diǎn)，在空間上對(duì)變量離散保留了Box method的方法，而時(shí)間上提出采用向后差分的方法避免了時(shí)間上離散求解的不穩(wěn)定性，然而這種方法時(shí)間上只有一階精度.Sun等［10］在時(shí)間上采用廣義梯形方法來優(yōu)化求解的穩(wěn)定性，具有一定的數(shù)值耗散［11］性，但這種方法只有在其數(shù)值耗散性能最差的時(shí)候才具有二階精度.Gobat等［2-3，12］提出了具有二階精度的廣義α算法，并將之應(yīng)用到二維拖纜的求解過程中.

本人在前人研究的基礎(chǔ)上，對(duì)三維拖纜非線性動(dòng)態(tài)方程在時(shí)間域上采用廣義α算法離散，分析證明了該方法求解拖纜方程在時(shí)間上具有二階精度，將該方法應(yīng)用到三維拖纜方程的求解，并通過將該方法與幾種常用的方法進(jìn)行仿真實(shí)驗(yàn)對(duì)比，分析了其求解的穩(wěn)定性及精確性.

1 拖纜的控制方程

采用Ablow所建立的模型，假設(shè)拖纜為圓柱型剛性纜，忽略拖纜彎曲應(yīng)力［13］和海流的影響得到的方程如下:

式中:s為未伸長(zhǎng)纜單位長(zhǎng)度，t為時(shí)間變量，m為單位長(zhǎng)度未伸長(zhǎng)纜的質(zhì)量，m1為纜的附加質(zhì)量與纜的質(zhì)量的和，A為未伸長(zhǎng)纜的截面積，E為彈性模量，ρ為水的密度，ω0為單位長(zhǎng)度纜在水中的重量，θ、φ為歐拉角度，V(Vt，Vn，Vb)局部坐標(biāo)中的速度矢量，Rd(Rd1，Rd2，Rd3)為水動(dòng)力.

定義Y(s，t)=(T，Vt，Vn，Vb，θ，φ)，因此方程(1)～(6)寫成矩陣的形式為

2 拖纜方程的離散

2.1 空間的離散

這里拖纜方程在空間上的離散采用Ablow和Schechter所用的基于二階精度的有限差分 Box method方法，假設(shè)把纜劃分為n個(gè)節(jié)點(diǎn)，那么拖纜的動(dòng)態(tài)方程就可以離散為n－1個(gè)矩陣方程:

2.2 時(shí)間的離散

2.2.1 廣義α算法

在時(shí)間上的離散采用廣義 α算法，Gobat和 Grosenbaugh［11］證明了有限差分Box method在時(shí)間離散上存在著求解的不穩(wěn)定性，并把廣義α算法應(yīng)用到二維拖纜的求解問題.廣義α算法在時(shí)間上離散形式為

其中，

式中:Δt表示時(shí)間步長(zhǎng)，αm、αk和γ則構(gòu)成了廣義α算法.

為了比較方便的證明廣義α算法的精度及穩(wěn)定性，考慮單一自由度問題，則有

其中，ω為頻率，式(16)即為下列方程的差分形式:

假設(shè)y(ti)是常微分方程(17)的一個(gè)精確解，那么截?cái)嗾`差可以定義為

其中，Bn、Cn為式(16)中yi和fi的系數(shù)，將yi、fi用泰勒級(jí)數(shù)展開，并代入式(17)經(jīng)過運(yùn)算化簡(jiǎn)可得

可以看出當(dāng)αm－αk+γ=0.5時(shí)，廣義α算法具有二階精度.為了便于討論廣義α算法的穩(wěn)定性，將廣義α算法寫成增廣矩陣的形式:

式中增廣矩陣定義為

下面，通過分析矩陣的特征值λ1，2以及譜半徑ρ來檢測(cè)算法的數(shù)值耗散性能，增廣矩陣的譜半徑定義為

假設(shè)當(dāng)ωΔt→∞增廣矩陣的特征值λ1，2為2個(gè)相等的實(shí)數(shù)值，那么式(21)就可以化為只含有1個(gè)參數(shù)的方程，定義λ∞為ωΔt∞時(shí)矩陣的雙重特征值，則αm、αk和γ可由以下方程來計(jì)算:

由以上廣義算法的形式可知，有限差分 Box method、廣義梯形法以及向后差分法都是廣義α算法的特殊形式，只是αm、αk和γ的取值不同而已，表1給出了幾個(gè)不同算法αm、αk和γ的值.

表1 廣義α算法的幾種特例Table 1 Several special cases for generalized-α algorithm

圖1 廣義α算法及其幾種特例的譜半徑Fig.1 The spectral radius of generalized-α algorithm and several special cases

圖1所示為表1內(nèi)的幾種廣義α算法的譜半徑圖，從圖中可知，有限差分Box method(λ∞=－1.0)就像一個(gè)全通濾波器，它沒有任何的數(shù)值耗散.向后差分算法子在較高頻段衰減曲線雖然比較光滑，但是這種算法只具有一階精度.當(dāng)λ∞=0.0時(shí)廣義α算法具有二階精度，但是其譜半徑圖上有個(gè)拐點(diǎn)，這說明增廣矩陣的一個(gè)特征值從拐點(diǎn)開始轉(zhuǎn)變成2個(gè)共軛的特征值，與假設(shè)不符.廣義梯形法對(duì)中低頻濾波效果明顯，然而在較高的頻率段，梯形法則幾乎沒有數(shù)值耗散.Gobat等［2-3，12］經(jīng)過仿真研究對(duì)比發(fā)現(xiàn)，對(duì)于大多數(shù)拖纜問題，當(dāng)λ∞在［－0.3，－0.7］取值時(shí)，雖然增廣矩陣具有2個(gè)不相等的實(shí)特征值，但其譜半徑隨著頻率的增加逐漸減小，或者僅有微小的增加，具有較好的穩(wěn)定性.這里在仿真計(jì)算時(shí)取λ∞=－0.5，其譜半徑見圖1所示.

2.2.2 廣義α算法求解非線性問題

廣義α算法在求解線性問題時(shí)，由于系數(shù)矩陣M、K是不變的，所以可以用式(12)進(jìn)行計(jì)算，但是對(duì)于非線性問題系數(shù)矩陣則隨時(shí)間的變化而改變，那么算法的穩(wěn)定性將隨時(shí)間有條件收斂.

為了避免這個(gè)問題，提高算法的穩(wěn)定性，系數(shù)矩陣必須取固定值，普遍采用相同權(quán)重的方法來計(jì)算平均系數(shù)矩陣，其是用速度和位移矢量的平均值來計(jì)算.廣義α算法矩陣M和矩陣K的系數(shù)分別用αm、αk表示，這樣算法的表達(dá)形式為

3 數(shù)值仿真及對(duì)比分析

3.1 邊界條件

1)拖纜首端邊界條件:

首端纜上拖點(diǎn)的速度與拖船速度相同，即

式中:Vx、Vy、Vz為拖船航速在固定坐標(biāo)系X、Y、Z方向上的速度分量.

2)尾端自由端邊界條件:

自由端張力為零，即T=0.

3.2 數(shù)值仿真及對(duì)比分析

為了分析廣義α算法的精度和穩(wěn)定性，這里僅考慮回轉(zhuǎn)運(yùn)動(dòng)情況.計(jì)算的狀態(tài)為:首先拖船沿某一固定方向以18.5 kn的速度直線運(yùn)動(dòng)1 s，然后進(jìn)入半徑為0.64 km的回轉(zhuǎn)運(yùn)動(dòng)440 s，最后，在完成375°的回轉(zhuǎn)后，沿圓周的切線直線運(yùn)動(dòng)300 s.在計(jì)算時(shí)，所用的時(shí)間步長(zhǎng)和整個(gè)系統(tǒng)的分段數(shù)均與文獻(xiàn)［4］相同.這里參數(shù)設(shè)定部分參考文獻(xiàn)［4］，表2如下.

表2 拖纜陣列參數(shù)Table 2 Towed array parameters

圖2 幾種常用方法計(jì)算結(jié)果Fig.2 Calculation results of several commonly used methods

圖2給出了可用廣義α算法表示的幾種算法計(jì)算拖纜回轉(zhuǎn)運(yùn)動(dòng)的結(jié)果圖，利用表3對(duì)仿真結(jié)果進(jìn)行對(duì)比分析比較，其中Rispin的數(shù)據(jù)為實(shí)驗(yàn)數(shù)據(jù)，來自于文獻(xiàn)［4］.考慮到實(shí)驗(yàn)的誤差、海流、船舶等的影響，廣義α算法和差分Box method結(jié)果基本和實(shí)驗(yàn)相符，而向后差分、廣義梯形法由于算法精度的影響與實(shí)驗(yàn)有少許的偏差.

表3 實(shí)驗(yàn)與幾種計(jì)算方法結(jié)果的比較Table 3 Comparison of experimental and numerical resultsm

為了說明算法收斂性的優(yōu)越，圖3給出了幾種算法計(jì)算拖纜回轉(zhuǎn)運(yùn)動(dòng)的每一點(diǎn)的平均誤差圖，表4記錄了各種算法的運(yùn)算時(shí)間，這里空間步長(zhǎng)取2 m，時(shí)間步長(zhǎng)與文獻(xiàn)［4］相同均為20 s，仿真總時(shí)間為37 s，所用的計(jì)算機(jī)參數(shù)為:CPU為 Inter Pentium Dual E2160主頻為1.79 Hz，內(nèi)存為1.99 GB.仿真實(shí)驗(yàn)分析表明，廣義α算法、向后差分以及廣義梯形法比差分Box method具有較好的收斂性，因此計(jì)算它們的計(jì)算時(shí)間相對(duì)于差分Box method有很大的優(yōu)勢(shì).

表4 幾種方法的運(yùn)算時(shí)間比較Table 4 Operation time of several methods

圖3 幾種計(jì)算收斂性圖Fig.3 Convergence of several methods

通過以上對(duì)廣義α算法、向后差分、廣義梯形和有限差分Box method仿真計(jì)算對(duì)比分析研究，可以發(fā)現(xiàn)廣義α算法具有收斂效果好，計(jì)算時(shí)間短，計(jì)算精度高等優(yōu)點(diǎn)，為求解非線性拖纜方程提供了一個(gè)較好的解決方法.

4 結(jié)束語

本文通過數(shù)學(xué)運(yùn)算證明了廣義α算法具有二階精度，對(duì)比分析了有限差分Box method、廣義梯形算法，向后差分算法以及不同的值廣義α算法的穩(wěn)定性和精確性，研究表明，在－0.3到－0.7之間取值時(shí)，廣義α算法對(duì)拖曳陣列有較好的求解穩(wěn)定性.將廣義α算法用于拖纜的非線性動(dòng)態(tài)方程求解，在空間上保留有限差分Box method的變量離散的方法，在時(shí)間上采用廣義α算法變量離散方法，并與向后差分、廣義梯形和有限差分Box method進(jìn)行對(duì)比仿真實(shí)驗(yàn)，仿真結(jié)果表明，廣義α算法不僅有較好的求解準(zhǔn)確性，而且有較短的仿真運(yùn)算求解時(shí)間，為提高我國拖曳陣列的控制技術(shù)提供了一個(gè)較好的解決途徑.

［1］張維競(jìng)，張小卿，陳峻.基于嵌入式水鳥的海洋地震拖纜運(yùn)動(dòng)狀態(tài)仿真研究［J］.海洋工程，2009，27(4):81-86.

ZHANG Weijing，ZHANG Xiaoqing，CHEN Jun.Simulation study on the movement states of marine seismic streamer based on embedded birds［J］.The Ocean Engineering，2009，27(4):81-86.

［2］GOBAT J I，GROSENBAUGH M A.Applications of the generalized-α method to the time integration of the cable dynamics equations［J］.Computer Methods in Applied Mechanics and Engineering，2001，190(6):4817-4829.

［3］GOBAT J I，GROSENBAUGH M A.Dynamics in the touchdown region of catenary moorings［J］.International Journal of Offshore and Polar Engineering，2001，11(1): 273-281.

［4］ABLOW C M，SCHECHTER S.Numerical simulation of undersea cable dynamics［J］.Ocean Engineering，1983，10 (6):443-457.

［5］PARK H I，JUNG D H，KOTERAYAMA W.A numerical and experimental study on dynamics of a towed low tension cable［J］.Applied Ocean Research，2003，25(10):289-299.

［6］HUO Cunfeng，YAO Baoheng，F(xiàn)U Bin，et al.Investigation on transient dynamic behaviors of low-tension undersea cables［J］.Journal of Shanghai Jiaotong University，2011，16 (1):34-39.

［7］CHATJIGEORGIOU I K.3 dimensional nonlinear dynamics of submerged，extensible catenary pipes conveying fluid and subjected to end-imposed excitations［J］.International Journal of Non-Linear Mechanics，2010，45(4):667-680.

［8］GROSENBAUGH M A.Transient behavior of towed cable systems during ship turning maneuvers［J］.Ocean Engineering，2007，34(8):1532-1542.

［9］KOH C G，ZHANG Y，QUEK S T.Low-tension cable dynamics:numerical and experimental studies［J］.Journal of Engineering Mechanics，1999，125(3):347-354.

［10］SUN Y，LEONARD J W，CHIOU R B.Simulation of unsteady oceanic cable deployment by direct integration with suppression［J］.Ocean Engineering，1994，21(6):243-256.

［11］馬秀騰，陳立平，張?jiān)魄?約束力學(xué)系統(tǒng)運(yùn)動(dòng)方程積分的數(shù)值耗散研究［J］.系統(tǒng)仿真學(xué)報(bào)，2009，21(20): 6373-6377.

MA Xiuteng，CHEN Liping，ZHANG Yunqing.Research on numerical dissipation of integration methods for motion equations in constrained mechanical systems［J］.Journal of System Simulation，2009，20:6373-6377.

［12］GOBAT J I，GROSENBAUGH M A.Time-domain numerical simulation of ocean cable structures［J］.Ocean Engineering，2006，33(10):1373-1400.

［13］王飛，陳錦標(biāo)，涂興華.低速低應(yīng)力時(shí)拖攬運(yùn)動(dòng)仿真［J］.上海交通大學(xué)學(xué)報(bào)，2010，44(6):828-832.

WANG Fei，CHEN Jingbiao，TU Xinghua.Numerical simulation of towed low-tension cable［J］.Journal of Shanghai Jiaotong University，2010，44(6):828-832.