高帥，尹勇，孫霄峰，神和龍

(大連海事大學航海學院，遼寧大連116026)

曳綱作為水下拖纜的一種，其水動力研究可借鑒現(xiàn)有的水下拖纜相關理論成果，根據(jù)目前國內外發(fā)表的相關文獻可知，在拖曳系統(tǒng)的數(shù)學模型中應用最為廣泛的有兩類:有限差分法和集中質量法。Ablow等［1］首次用有限差分法求解拖纜的三維運動，Howell［2］、張潞怡［3］、Sun 等［4］、李英輝［5］、孫霄峰［6-7］對有限差分法進行了進一步修正;Walton等［8］首先采用集中質量法研究了海軍武器試驗中水下錨鏈的二維運動響應，Delmer等［9］、Huang［10］、王飛［11］、陳英龍［12］先后對集中質量法進行了改進。國內外學者在進行水下拖纜建模時，大多只考慮模型的精準性及其與模型試驗的驗證，沒有進行可視化和實時性方面的研究。而模擬器中的視景系統(tǒng)為操作者提供一個訓練仿真環(huán)境，其可視化效果對環(huán)境真實感和訓練效果具有重要意義。本研究中，在現(xiàn)有拖曳系統(tǒng)水動力理論研究的基礎上，忽略曳綱的扭轉、抖動、彎矩等，基于集中質量法建模思想，針對拖曳整體運動進行分析，建立拖網(wǎng)漁船曳綱的動態(tài)運動模型，并考慮視景顯示的實時性要求，利用三維圖形引擎OSG(OpenScene-Graph)實現(xiàn)拖曳系統(tǒng)的三維可視化。

1 曳綱穩(wěn)態(tài)運動

當拖網(wǎng)漁船在水中穩(wěn)定直航時，系統(tǒng)處于相對靜止的狀態(tài)，此時曳綱的運動狀態(tài)稱為曳綱的穩(wěn)態(tài)解。穩(wěn)態(tài)運動方程是一個復雜的非線性偏微分方程組，通過解析解很難得到，只有通過數(shù)值方法來近似求解。穩(wěn)態(tài)解是動態(tài)運動方程的初始條件。

1.1 坐標系統(tǒng)

為進行曳綱的運動分析計算，采取圖1兩種坐標系:曳綱慣性坐標系o-ijk和局部坐標系ptnb。在o-ijk中，o為拖點位置，oi、oj軸分別平行于船艏向和船舶右舷正橫方向;ok軸垂直水平面向下;在p-tnb中，p為曳綱上一點，pt軸向為曳綱在p點處的切線方向，pn、pb方向分別為在p點處的兩個法線方向，pb方向平行于oi、oj所組成的平面。兩個坐標系統(tǒng)可以通過拖纜姿態(tài)角(θ，φ)相互轉換，其中θ為偏航角，φ為俯仰角。它們的轉換關系以矩陣形式表示為

其中，

1.2 曳綱穩(wěn)態(tài)運動數(shù)學模型

曳綱為細長、柔性的圓柱體，忽略彎矩和扭矩

圖1 拖曳系統(tǒng)坐標示意圖Fig.1 Schematic diagram of a towed system

其中:T為微元段曳綱的張力，指向曳綱切向;B、G分別為單位長度曳綱的浮力和重力;F為曳綱受到的流體阻力;S為微元段曳綱拉伸后的弧長。

將式 (2)在曳綱局部坐標系下展開，則平衡方程可寫成標量形式［11］:

其中:w為曳綱水中質量;ρ為海水密度;d為曳綱截面的直徑;s表示未拉伸的曳綱弧長;ε為曳綱的應變;Ct和 Cn為切向和法向阻力系數(shù);ut、un、ub為局部坐標系下的速度分量，由系統(tǒng)相對于水流的拖曳速度轉換到局部坐標系下得到，如下式:

式中，v1為拖曳速度，v2為流速。而曳綱在慣性坐標系下的坐標由下列關系式給出:

給定曳綱在端點處或者任意位置上的狀態(tài)值(T0，θ0，φ0)，便可由式 (3)通過四階龍格庫塔法進行積分求解，再由式 (5)求出曳綱分段的位置信息。對于兩點邊值的穩(wěn)態(tài)問題，邊界條件不能直接應用到數(shù)值求解中，聯(lián)立方程將轉換為關于初始值的隱式方程組，其中初始值(T0，θ0，φ0)為待求量。

1.3 穩(wěn)態(tài)初始值及方程求解

穩(wěn)態(tài)運動控制方程的求解，需要在曳綱首尾兩端加上一定的邊界條件，將其轉變?yōu)槌跏贾祮栴}再進行求解。考慮拖船以某一航速勻速直航，在拖曳首端，拖船拖點和曳綱的速度與位置始終一致。此時，ut、un、ub由式 (4)可以直接求出。為驗證穩(wěn)態(tài)模型的合理性，先假設均勻海流下拖曳尾端不受漁網(wǎng)的拉力而是自由端，此時，T0=0(T=0對方程 (3)來說為一奇點，可取一接近零的數(shù)來代替)和θ0=φ(φ為船艏向)。將T0=0代入式(3)可求得φ0值，如下式所示:

自由端的邊界條件表示成一個顯式的初始值方程，將(T0，θ0，φ0)代入方程 (3)通過四階龍格庫塔法進行積分求解，而曳綱的首端與拖船相連，拖點的位置 (x0，y0，z0)可認為已知，這樣由式(5)求出曳綱分段的位置信息。

為了比較的一致性，考慮均勻流下穩(wěn)定直航狀態(tài)，在進行曳綱的分析時采用了文獻 ［4］的物理參數(shù)，如表1所示。表2為本研究結果與文獻［4］的參數(shù)對照，可以看出，拖點的俯仰角φ、張力T和曳綱尾端的深度基本一致。曳綱長度不變，尾端為自由端時，仿真計算結果 (圖2)表明，隨著拖速的增加，拖曳的深度降低;圖3表明，隨著拖速的增加，拖點張力增大。

表1 曳綱物理參數(shù)Tab.1 Physical parameters of the warp

表2 拖曳穩(wěn)定后相關參數(shù)對照Tab.2 Parameter contrast of the trawling stabilization

圖2 曳綱陣形隨拖速的變化Fig.2 Displayed warp at different towing speeds of trawling

圖3 均勻海流不同拖速的張力分布Fig.3 Distribution of strain at different towing speeds in uniform current

2 曳綱動態(tài)運動

2.1 曳綱動態(tài)運動方程

在曳綱慣性坐標系下對曳綱動態(tài)運動進行分析，不考慮彎矩和扭矩的影響，每段的作用力 (重力G、浮力B、流體阻力F和張力T)均集中作用在節(jié)點上。對曳綱第i個節(jié)點應用牛頓第二定律，建立在流體中運動的曳綱節(jié)點i的動力方程:

其中:Mi和Δ Mi分別為節(jié)點i的質量和附加質量;為節(jié)點i的加速度。

(1)附加質量。流體中的物體加速運動會引起周圍流體也做加速運動，流體的質量會引起一個反作用力，相當于物體具有了一個附加質量，即

其中:Cmt、Cmn、Cmb分別為局部坐標系各方向的附加質量系數(shù)，參考文獻 ［12］中的值，Cmt=0，Cmn=Cmb=1．0;Vi表示節(jié)點i段的體積。

(2)張力。繩索在拉伸后具有一定的張力。當曳綱尾部為自由端時，尾部節(jié)點只受到上一節(jié)點的張力作用;當曳綱尾部為拖網(wǎng)時，尾部節(jié)點與漁網(wǎng)手綱相連，受到網(wǎng)具的阻力作用。

其中:Tij為曳綱在慣性坐標系受到的張力矢量;E為彈性模量;Aij為i段的橫截面積;lij和loij分別為微元段的實際長度和初始長度;e=(lij-loij)/loij。

(3)流體阻力。參照Albow等［1］使用的方法，將拖纜阻力進行簡化，分為切向和法向兩部分處理:

在曳綱慣性坐標系下，只在垂直方向受到重力和浮力作用，由此，曳綱的動態(tài)運動方程 (7)可寫成下式:

2.2 邊界條件

(1)拖點邊界條件。在拖曳的上端點，漁船對曳綱的影響主要是改變了曳綱首節(jié)點的位置和速度。設導纜孔在隨船坐標系下的坐標為(xw，yw)，轉換到空間坐標系下可表示為

其中:Xw、Yw為導纜孔在空間坐標系下的坐標;x、y為漁船重心在空間坐標系下的坐標;φ為船艏向。

導纜孔在隨船運動坐標系下的速度可表示為

其中:uw、vw分別為導纜孔處的縱向和橫向速度;u、v、r分別為漁船重心處的縱向速度、橫向速度和轉艏角速度;α=tan-1(yw/xw)。

(2)尾端邊界條件。曳綱末端與漁網(wǎng)手綱相連，為簡化處理，漁網(wǎng)的阻力采用日本小山武夫網(wǎng)具阻力公式［13］:

其中:Fz為網(wǎng)具阻力 (N);d/a為網(wǎng)具線面積系數(shù)，即網(wǎng)線直徑與目腳長度之比;LW0為網(wǎng)具拉直的總長 (m);CW0為網(wǎng)具網(wǎng)口拉直的周長 (m);v3為相對拖速 (m/s)。

將漁網(wǎng)視為一個點融入到曳綱尾部結點i=0的控制方程中，此時邊界條件可近似為

其中，mZ為網(wǎng)具質量。

2.3 曳綱動態(tài)方程數(shù)值求解

聯(lián)立控制方程 (7)和相應的邊界條件及速度v=dx/dt，便組成一個完整的微分方程組:

采用四階龍格庫塔方法求解此微分方程，便可得到t+Δt時刻的位置和速度，其計算公式如下:

在求解此方程時，張力T是位移的隱式函數(shù)，由式 (17)可以看出，單次求解過程需要4次循環(huán)后得到節(jié)點的位移和速度，在每次循環(huán)后更新一次數(shù)據(jù)，用得到的臨時位置來計算各點的張力和流體阻力。用表1的參數(shù)進行動態(tài)仿真，尾部設為自由端，拖點在30 s內速度從0 m/s加速到 0．8 m/s，此后以0．8 m/s的速度前進，所得到的纜形如圖4所示，該纜形與文獻 ［4］采用有限差分法所得纜形基本一致。

在實際中，曳綱都是與漁網(wǎng)手綱相連的，而漁網(wǎng)的阻力在不考慮漁獲物的情況下受力較大，此時曳綱基本呈直線狀態(tài)［14］。

對曳綱進行空間離散，離散的段數(shù)直接決定了運動方程的復雜程度，即離散段數(shù)越多，預測結果精度就越高，但同時計算量就越大。在可視化過程中，既要保證仿真的準確性，又要考慮系統(tǒng)的實時性，從而選擇合適的空間步長。圖5為本研究中采用集中質量法和文獻 ［6］利用有限差分法建立模型的分段數(shù)與解算時間的關系，可以看出，集中質量法在實時性方面具有更高的效率。

圖4 自由端拖纜陣形的變化Fig.4 Displayed configurations of towed cable with free end

圖5 模型分段數(shù)與解算時間關系Fig.5 Relationship between solution time and model segments

至此，建立了曳綱的穩(wěn)態(tài)運動模型和動態(tài)控制方程，動態(tài)運動是以穩(wěn)態(tài)運動作為初始值，進行積分求解得到每個節(jié)點的速度和位置，由于采用龍格庫塔法求解集中質量模型需滿足一定條件才能穩(wěn)定，即時間步長Δt必須滿足一定的穩(wěn)定性條件，由文獻［11］可知，數(shù)值穩(wěn)定性的充要條件為

Eσ/m·Δ t2/Δ s2＜1(σ為橫截面積)。

3 拖曳系統(tǒng)的可視化

利用三維圖形引擎OSG開發(fā)了拖網(wǎng)曳綱的視景驅動程序，根據(jù)曳綱的數(shù)學模型解算得到各節(jié)點的位置信息進行實時繪制，實現(xiàn)曳綱的三維動態(tài)可視化，整個系統(tǒng)的仿真流程如圖6所示。

圖6 仿真流程圖Fig.6 Flowing chart of simulation

為增強可視化效果，加入海浪、天空、島嶼等進行渲染，但并不考慮曳綱受到的波浪力。漁船在航行中，船上物體會隨漁船發(fā)生振蕩，這要求曳綱的物理模型要進行相應的矩陣變換。假定漁船的位置矩陣為M，M是一個3×3的矩陣，海面波浪的坐標矩陣為A(posX posY h)，這里(posX，posY)為船所在的海平面坐標，h為該點的波高，則船在波浪里的位置矩陣為M':M'=M·A，船上物體相對船中心的位置為M0(X0Y0Z0)，物體隨船振蕩的坐標矩陣表示為P=M0·M'。這樣便可使得物體隨船振蕩，如圖7所示。

在仿真過程中，將漁網(wǎng)視作水下拖體，不考慮漁網(wǎng)的數(shù)學模型。OSG采用一種自上而下的、分層的樹狀數(shù)據(jù)結構來組織空間數(shù)據(jù)集，以提升渲染效率。將漁船、曳綱和漁網(wǎng)分別作為獨立的節(jié)點，這些節(jié)點中包含了幾何信息和用于控制其外觀的渲染狀態(tài)信息，通過回調函數(shù)進行實時渲染。

將曳綱離散為許多相連的圓柱體，并將紋理對象映射到每個圓柱曲面上，選取紋理對象時首先應注意單個紋理映射到圓柱體上的融合，其次為使相鄰兩個圓柱體連接時不會出現(xiàn)紋理圖像的斷裂，紋理圖片的上側與下側應保證對稱［7］。曳綱紋理可視化結果如圖8～圖10所示。

圖7 物體隨船振蕩Fig.7 Object shaking with ship

圖8 曳綱局部可視化效果Fig.8 Local visualization of warp

圖9 雙船拖曳運動 (俯視)Fig.9 Drag motion of a pair trawler(vertical view)

圖10 雙船拖曳運動 (側視)Fig.10 Drag motion of a pair trawler(side view)

整個場景的渲染過程中，除去曳綱后系統(tǒng)的渲染幀率基本維持在60幀，這樣曳綱模型的解算時間成為影響整個系統(tǒng)實時性的主要因素。分析圖5選擇合適的空間步長，以滿足視景系統(tǒng)實時性的要求。

4 結語

本研究中在現(xiàn)有拖曳系統(tǒng)水動力理論研究的基礎上，采用集中質量法建立了曳綱的運動數(shù)學模型，最后利用三維圖形引擎OSG實現(xiàn)拖曳系統(tǒng)的三維可視化。通過和已有的試驗數(shù)據(jù)對比，證明該數(shù)學模型是合理的，在三維可視化中，其顯示的速率也得到保證，從而為漁具系統(tǒng)的水動力模型的建立和視景可視化奠定了基礎。

［1］Ablow C M，Schechter S．Numerical simulation of undersea cable dynamics［J］．Oceanic Engineering，1983，10(6):443-457．

［2］Howell C T．Investigation of the dynamics of low-tension cables［D］．Massachusetts:Massachusetts Institute of Technology and Woods Hole Oceanographic Institution，1992．

［3］張潞怡．6000米深拖系統(tǒng)非線性運動理論研究及仿真［D］．上海:上海交通大學，1996．

［4］Sun Yang，Lenard J W．Dynamics of ocean cables with local low tension regions［J］．Ocean Engng，1998，25(6):443-463．

［5］李英輝．合成孔徑聲納拖曳系統(tǒng)運動的理論與實驗研究［D］．哈爾濱:哈爾濱工程大學，2002．

［6］孫霄峰．單船中層拖網(wǎng)系統(tǒng)的建模與仿真［D］．大連:大連海事大學，2008．

［7］孫霄峰，高帥，尹勇，等．漁船模擬器中中層拖網(wǎng)的建模與仿真［J］．大連海洋大學學報，2012，27(3):284-288．

［8］Walton T S，Polachek H．Calculation of transient motion of submerged cables［J］．Mathematics of Computation，1960，14:27-46．

［9］Delmer T N，Stephens T C，Coe J M．Numerical simulation of towed cables［J］．Ocean Engineering，1983，10:119-132．

［10］Huang S．Dynamic analysis of three-dimensional marine cables［J］．Ocean Engineering，1994，21:587-605．

［11］王飛．海洋勘探拖曳系統(tǒng)運動仿真與控制技術［D］．上海:上海交通大學，2006．

［12］陳英龍．拖網(wǎng)捕撈拖曳系統(tǒng)的建模及控制研究［D］．杭州:浙江大學，2013．

［13］崔建章．漁具與漁法學［M］．北京:中國農業(yè)出版社，1997．

［14］陳英龍，趙勇剛，周華，等．大型中層拖網(wǎng)網(wǎng)具系統(tǒng)的仿真研究［J］．浙江大學學報:工學版，2014，48(4):625-632．