拖網(wǎng)漁船曳綱的運動建模及仿真

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

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

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

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

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

1.1 坐標(biāo)系統(tǒng)

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

其中，

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

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

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

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

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

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

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

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

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

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

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

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

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

表2 拖曳穩(wěn)定后相關(guā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)運動方程

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

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

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

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

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

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

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

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

2.2 邊界條件

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

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

導(dǎo)纜孔在隨船運動坐標(biāo)系下的速度可表示為

其中:uw、vw分別為導(dǎo)纜孔處的縱向和橫向速度;u、v、r分別為漁船重心處的縱向速度、橫向速度和轉(zhuǎn)艏角速度;α=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)視為一個點融入到曳綱尾部結(jié)點i=0的控制方程中，此時邊界條件可近似為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

將曳綱離散為許多相連的圓柱體，并將紋理對象映射到每個圓柱曲面上，選取紋理對象時首先應(yīng)注意單個紋理映射到圓柱體上的融合，其次為使相鄰兩個圓柱體連接時不會出現(xiàn)紋理圖像的斷裂，紋理圖片的上側(cè)與下側(cè)應(yīng)保證對稱［7］。曳綱紋理可視化結(jié)果如圖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 雙船拖曳運動 (側(cè)視)Fig.10 Drag motion of a pair trawler(side view)

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

4 結(jié)語

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

［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］．上海:上海交通大學(xué)，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］．哈爾濱:哈爾濱工程大學(xué)，2002．

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

［7］孫霄峰，高帥，尹勇，等．漁船模擬器中中層拖網(wǎng)的建模與仿真［J］．大連海洋大學(xué)學(xué)報，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)運動仿真與控制技術(shù)［D］．上海:上海交通大學(xué)，2006．

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

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

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