廣義α算法在拖曳陣列動態仿真中的應用

張廣磊，張維競，劉濤

(上海交通大學 海洋工程國家重點實驗室，上海 200240)

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

本人在前人研究的基礎上，對三維拖纜非線性動態方程在時間域上采用廣義α算法離散，分析證明了該方法求解拖纜方程在時間上具有二階精度，將該方法應用到三維拖纜方程的求解，并通過將該方法與幾種常用的方法進行仿真實驗對比，分析了其求解的穩定性及精確性.

1 拖纜的控制方程

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

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

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

2 拖纜方程的離散

2.1 空間的離散

這里拖纜方程在空間上的離散采用Ablow和Schechter所用的基于二階精度的有限差分 Box method方法，假設把纜劃分為n個節點，那么拖纜的動態方程就可以離散為n－1個矩陣方程:

2.2 時間的離散

2.2.1 廣義α算法

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

其中，

式中:Δt表示時間步長，αm、αk和γ則構成了廣義α算法.

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

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

假設y(ti)是常微分方程(17)的一個精確解，那么截斷誤差可以定義為

其中，Bn、Cn為式(16)中yi和fi的系數，將yi、fi用泰勒級數展開，并代入式(17)經過運算化簡可得

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

式中增廣矩陣定義為

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

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

由以上廣義算法的形式可知，有限差分 Box method、廣義梯形法以及向后差分法都是廣義α算法的特殊形式，只是αm、αk和γ的取值不同而已，表1給出了幾個不同算法α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內的幾種廣義α算法的譜半徑圖，從圖中可知，有限差分Box method(λ∞=－1.0)就像一個全通濾波器，它沒有任何的數值耗散.向后差分算法子在較高頻段衰減曲線雖然比較光滑，但是這種算法只具有一階精度.當λ∞=0.0時廣義α算法具有二階精度，但是其譜半徑圖上有個拐點，這說明增廣矩陣的一個特征值從拐點開始轉變成2個共軛的特征值，與假設不符.廣義梯形法對中低頻濾波效果明顯，然而在較高的頻率段，梯形法則幾乎沒有數值耗散.Gobat等［2-3，12］經過仿真研究對比發現，對于大多數拖纜問題，當λ∞在［－0.3，－0.7］取值時，雖然增廣矩陣具有2個不相等的實特征值，但其譜半徑隨著頻率的增加逐漸減小，或者僅有微小的增加，具有較好的穩定性.這里在仿真計算時取λ∞=－0.5，其譜半徑見圖1所示.

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

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

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

3 數值仿真及對比分析

3.1 邊界條件

1)拖纜首端邊界條件:

首端纜上拖點的速度與拖船速度相同，即

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

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

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

3.2 數值仿真及對比分析

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

表2 拖纜陣列參數Table 2 Towed array parameters

圖2 幾種常用方法計算結果Fig.2 Calculation results of several commonly used methods

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

表3 實驗與幾種計算方法結果的比較Table 3 Comparison of experimental and numerical resultsm

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

表4 幾種方法的運算時間比較Table 4 Operation time of several methods

圖3 幾種計算收斂性圖Fig.3 Convergence of several methods

通過以上對廣義α算法、向后差分、廣義梯形和有限差分Box method仿真計算對比分析研究，可以發現廣義α算法具有收斂效果好，計算時間短，計算精度高等優點，為求解非線性拖纜方程提供了一個較好的解決方法.

4 結束語

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

［1］張維競，張小卿，陳峻.基于嵌入式水鳥的海洋地震拖纜運動狀態仿真研究［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，FU 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］馬秀騰，陳立平，張云清.約束力學系統運動方程積分的數值耗散研究［J］.系統仿真學報，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］王飛，陳錦標，涂興華.低速低應力時拖攬運動仿真［J］.上海交通大學學報，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.