漁船模擬器中拖網(wǎng)放網(wǎng)過(guò)程的數(shù)值模擬*

高 帥，尹 勇, 孫霄峰

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

在漁船操作模擬器的研究方面，英國(guó)Transas的NTPRO 5000[1]、俄羅斯Vector Marine Electronics公司的NFS-3000[2]等產(chǎn)品都可以用于漁船船員培訓(xùn)，能夠進(jìn)行漁船拖網(wǎng)和圍網(wǎng)作業(yè)仿真。利用漁船模擬器對(duì)相關(guān)人員進(jìn)行培訓(xùn)，可以讓船員在培訓(xùn)中熟練地操作各種船舶設(shè)備，提高應(yīng)對(duì)危機(jī)和海上安全管理的經(jīng)驗(yàn)?zāi)芰?。目前我?guó)漁船模擬器領(lǐng)域尚處于起步階段，還沒(méi)有可用于培訓(xùn)的產(chǎn)品。

漁網(wǎng)作為漁業(yè)模塊的重要組成部分，近年來(lái)其水動(dòng)力特性的研究已取得很大進(jìn)展。大多數(shù)學(xué)者采用集中質(zhì)量法建立漁網(wǎng)的數(shù)學(xué)模型，來(lái)進(jìn)行漁網(wǎng)的動(dòng)態(tài)模擬。Takagi、Suzuki[3-6]等針對(duì)特定的網(wǎng)片進(jìn)行仿真分析，將網(wǎng)片的仿真與水槽實(shí)驗(yàn)進(jìn)行比較，實(shí)現(xiàn)網(wǎng)片的三維顯示；Chun-Woo Lee[7-8]開(kāi)發(fā)了漁網(wǎng)自動(dòng)設(shè)計(jì)建模軟件，該軟件對(duì)拖網(wǎng)和圍網(wǎng)都適用；在國(guó)內(nèi)，李玉成[9]研究了網(wǎng)箱在水流作用下網(wǎng)衣的變形和應(yīng)力分布情況；孫霄峰[10]并實(shí)現(xiàn)了單船中層拖網(wǎng)的實(shí)時(shí)仿真；黃小華[11]研究了不同配重和流速下網(wǎng)的受力和變形；陳英龍[12]考慮了拖網(wǎng)網(wǎng)板及升力帆布水動(dòng)力的作用建立整個(gè)網(wǎng)具系統(tǒng)的仿真模型，通過(guò)海上試驗(yàn)驗(yàn)證了模型的準(zhǔn)確性。漁網(wǎng)各質(zhì)量點(diǎn)的運(yùn)動(dòng)方程相互聯(lián)系，構(gòu)成數(shù)目巨大的非線性微分方程組，這就給模型的數(shù)值解法提出了很高的要求，常用的求解方法主要有：四階Runge-Kutta法、六級(jí)五階Runge-Kutta法及Newmark-β法。Runge-Kutta是一種顯式的數(shù)值解法, 其穩(wěn)定求解的時(shí)間步長(zhǎng)受到柔性體剛度的極大限制；而隱式的Newmark-β法其計(jì)算精度和模型解算時(shí)間有待于進(jìn)一步提高。

鐘萬(wàn)勰[13]提出線性定常結(jié)構(gòu)的精細(xì)時(shí)程積分方法，以其高精度、無(wú)條件穩(wěn)定等優(yōu)點(diǎn)得到廣泛應(yīng)用。而非線性問(wèn)題研究的關(guān)鍵是如何求解非齊次項(xiàng)產(chǎn)生的Duhamel積分。譚述君[14]在Duhamel項(xiàng)精細(xì)積分方法的基礎(chǔ)上構(gòu)造了單步法和多步法進(jìn)行非線性微分方程的求解；高強(qiáng)[15]則提出了一種針對(duì)大規(guī)模動(dòng)力系統(tǒng)的改進(jìn)的快速精細(xì)積分方法；鄒洋[16]通過(guò)精細(xì)積分方法對(duì)起重船吊物系統(tǒng)的三維非線性動(dòng)力方程進(jìn)行求解。

本文將采用集中質(zhì)量法建立漁網(wǎng)網(wǎng)衣的數(shù)學(xué)模型，采用精細(xì)積分法進(jìn)行模型的數(shù)值求解，并與模型試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，分別對(duì)水流作用下網(wǎng)衣所受的流體力、網(wǎng)衣形狀和位移隨時(shí)間的變化進(jìn)行了動(dòng)態(tài)模擬，驗(yàn)證了本文數(shù)值解法的可靠性。

1 漁網(wǎng)網(wǎng)衣的數(shù)學(xué)模型

根據(jù)集中質(zhì)量法將網(wǎng)衣離散為通過(guò)無(wú)質(zhì)量彈簧相連的質(zhì)量點(diǎn)集合。假設(shè)結(jié)節(jié)為球體，目腳為圓柱體，結(jié)節(jié)和目腳的質(zhì)量分別集中于各自的形心，如圖1所示。球體的流體阻力在各個(gè)方向上均相同，而作用在圓柱體的流體阻力與來(lái)流方向有關(guān)，以下分別建立結(jié)節(jié)和目腳的數(shù)學(xué)模型。

圖1 網(wǎng)衣模型示意圖Fig.1 Schematic diagram of the mesh

1.1 結(jié)節(jié)的運(yùn)動(dòng)方程

漁網(wǎng)網(wǎng)衣假設(shè)在均勻流作用下，由于球體的流體阻力系數(shù)在各方向上都相同，在空間坐標(biāo)系下對(duì)第i個(gè)結(jié)節(jié)進(jìn)行受力分析，如圖2所示。結(jié)節(jié)運(yùn)動(dòng)方程可表示如下：

(Mi+ΔMi)a=T+F+W+B。

(1)

其中：a是結(jié)節(jié)i的加速度；T、F分別表示結(jié)節(jié)受到的彈力和流體阻力，其中彈力是結(jié)節(jié)i受到與其相連的所有質(zhì)量點(diǎn)對(duì)其作用力的合力；W為結(jié)節(jié)的重力；B為結(jié)節(jié)所受的浮力；Mi和ΔMi表示第i個(gè)結(jié)節(jié)的質(zhì)量和附加質(zhì)量。

圖2 結(jié)節(jié)模型的受力示意圖Fig.2 Schematic diagram of the model for mesh knots

結(jié)節(jié)i與目腳j相連，其受的彈性力的大小可表示如下：

(2)

其中：E,Aij分別為質(zhì)量點(diǎn)i,j之間目腳的彈性模量和橫截面積；lij,l0ij分別為質(zhì)量點(diǎn)i,j之間的實(shí)際長(zhǎng)度和未伸長(zhǎng)時(shí)的長(zhǎng)度。作用在結(jié)節(jié)i上的彈性力在空間坐標(biāo)系三個(gè)坐標(biāo)軸上的分量可表示如下：

(3)

其中N表示同結(jié)節(jié)i相連的質(zhì)量點(diǎn)的集合。

考慮結(jié)節(jié)在均勻流作用下，其流體阻力在空間坐標(biāo)系三個(gè)坐標(biāo)軸上的分量可表示如下：

(4)

其中：ρ為流體密度；Cd為結(jié)節(jié)的阻力系數(shù)；S為結(jié)節(jié)的投影面積；vcx,vcy,vcz分別為流速在空間坐標(biāo)系三個(gè)坐標(biāo)軸上的分量。

結(jié)節(jié)在水中的重力表示為：

W+B=(mi-ρVi)g。

(5)

其中：Vi是結(jié)節(jié)的體積；g為重力加速度。

把以上(2)～(5)式代入(1)式可得結(jié)節(jié)i在空間坐標(biāo)系下的運(yùn)動(dòng)方程：

(6)

1.2 目腳的運(yùn)動(dòng)方程

目腳的數(shù)學(xué)模型與結(jié)節(jié)相似，只是結(jié)節(jié)是在空間坐標(biāo)系下進(jìn)行受力分析，而目腳在局部坐標(biāo)系下分析其受力(見(jiàn)圖3)。將目腳視作圓柱體，流體阻力系數(shù)與來(lái)流方向有關(guān)，目腳i的運(yùn)動(dòng)方程可表示如下：

圖3 目腳模型的受力示意圖Fig.3 Schematic diagram of the model for mesh bars

(7)

目腳i所受的流體阻力表示為：

(8)

其中：Sτ,Sη,Sξ分別為目腳在τ,η,ξ方向上的投影面積；vcτ,vcη,vcξ分別為流速在τ,η,ξ方向上的速度分量；Cdτ,Cdη,Cdξ分別為目腳在τ,η,ξ方向上的阻力系數(shù)。

方程(4)和(8)中的阻力系數(shù)Cd在雷諾數(shù)小于200時(shí)會(huì)發(fā)生急劇的改變，因此應(yīng)將Cd看作是雷諾數(shù)Re的函數(shù)[5]。

對(duì)于結(jié)節(jié)：Cd=101.2Re-0.6，

對(duì)于目腳：Cdτ=0.1，Cdη=Cdξ=100.7Re-0.3。

2 數(shù)值解法

精細(xì)積分法自鐘萬(wàn)勰提出來(lái)后，在結(jié)構(gòu)動(dòng)力分析中獲得了廣泛應(yīng)用。在此基礎(chǔ)上進(jìn)一步研究，使得非線性動(dòng)力方程的高精度計(jì)算成為可能。結(jié)構(gòu)系統(tǒng)非線性動(dòng)力方程的一般形式為：

(9)

(10)

其中：

B=-M-1K，G=-M-1C，I是單位矩陣。

方程(10)的解可寫(xiě)成如下形式：

(11)

對(duì)式(11)進(jìn)行數(shù)值離散，并假設(shè)時(shí)間步長(zhǎng)Δt=tk+1-tk，則第k+1步時(shí)的響應(yīng)公式可寫(xiě)成：

(12)

為使(12)式能進(jìn)行數(shù)值計(jì)算，選用計(jì)算格式較為簡(jiǎn)單同時(shí)精度又較高的辛普生積分法，也即：

(13)

(14)

式(14)即是計(jì)算系統(tǒng)響應(yīng)的迭代公式，因?yàn)門，T1/2都可以由數(shù)值方法精確算得，而F(t)又是已知載荷，因此通過(guò)這個(gè)迭代公式就能直接求出系統(tǒng)的響應(yīng)。

(15)

(16)

將X+f(t)看作外部激勵(lì)，對(duì)結(jié)節(jié)運(yùn)動(dòng)方程來(lái)說(shuō)：

(17)

而對(duì)目腳運(yùn)動(dòng)方程，其解算在局部坐標(biāo)系下，F(xiàn)(t)如下式：

(18)

(14)式右邊各項(xiàng)是矩陣相乘，對(duì)于本文來(lái)講，矩陣最高是6×6階，而矩陣相乘的計(jì)算效率很大程度上的影響了整個(gè)程序的運(yùn)行速度，所以對(duì)矩陣相乘算法進(jìn)行一些改進(jìn)是必要的。矩陣T，T1/2的形式相同，其構(gòu)造如下：

(19)

由此，將(14)式中的矩陣相乘寫(xiě)成矩陣各元素乘積的形式可得：

(20)

Vk=[x,y,z,vx,vy,vz]T，

F(t)=[0,0,0,Ftx,Fty,Ftz]T。

其中帶“′”的元素表示T1/2的相應(yīng)元素。

誠(chéng)然，政治究竟過(guò)不過(guò)硬，不由我們自己說(shuō)了算。黨的十九大報(bào)告指出，一個(gè)政黨，一個(gè)政權(quán)，其前途命運(yùn)取決于人心向背。民心是最大的政治。人民群眾反對(duì)什么、痛恨什么，我們就要防范和糾正什么。對(duì)十八大以來(lái)全面從嚴(yán)治黨取得的卓著成效，人民群眾給予了很高評(píng)價(jià)，之所以如此，追根溯源就在于我們緊緊抓住了人民群眾反對(duì)和痛恨的作風(fēng)問(wèn)題、腐敗問(wèn)題不放松，下大氣力、用真功夫加以解決和遏制，辦成了許多大事、解決了不少難題，從而贏得了黨心民心。

圖4是用精細(xì)積分法求解網(wǎng)衣數(shù)學(xué)模型的計(jì)算流程圖，Bar和Knot分別表示目腳、結(jié)節(jié)。在實(shí)際程序中直接用(20)式求解結(jié)節(jié)和目腳的位置及速度，模型的解算時(shí)間比直接用(14)式矩陣相乘提高很大。

3 模型驗(yàn)證

日本學(xué)者Takagi基于集中質(zhì)量法建立了漁網(wǎng)網(wǎng)衣三維動(dòng)態(tài)模擬[4-6]，采用六級(jí)五階Runge-Kutta法進(jìn)行數(shù)值求解，并且通過(guò)水槽試驗(yàn)驗(yàn)證了數(shù)值計(jì)算的可靠性。其水槽試驗(yàn)方法如下：將網(wǎng)片垂直放置于水槽中，水深為1.0 m。網(wǎng)片上緣的6個(gè)結(jié)節(jié)等距固定在1.5 m長(zhǎng)的鋼桿上，其左下角和右下角用沉子固定于水池底部，以不同速度的水流(8，20，41 cm/s)沿垂直于網(wǎng)片的方向流向網(wǎng)片，如圖5所示。

圖4 精細(xì)積分計(jì)算流程圖Fig.4 Precise integration computational flow chart

圖5 網(wǎng)片水槽試驗(yàn)示意圖Fig.5 Schematic diagram of the devices used in the flume tank experiments

本文選取其中一類網(wǎng)片B進(jìn)行數(shù)值模擬，該網(wǎng)片目大50 mm，橫向45目，縱向48目，其具體參數(shù)如表1所示。圖6是Takagi在試驗(yàn)中，用數(shù)碼相機(jī)拍攝到的網(wǎng)片在流速41.0 cm/s下達(dá)到穩(wěn)定后圖像。

圖6 水槽試驗(yàn)拍攝的網(wǎng)片圖像[5]Fig.6 Video image of the net configuration in the flume tank[5]

表1 測(cè)試網(wǎng)片B的計(jì)算參數(shù)Table 1 Calculation parameters of testing Net B

Note：①Number of mass points；②Mass of knot；③Mass of bar；④Bar；⑤Diameter of knot；⑥Length；⑦Diameter.

采用本文的精細(xì)積分方法進(jìn)行網(wǎng)片數(shù)學(xué)模型的求解，在均勻流41.0 cm/s條件下網(wǎng)片穩(wěn)定后得到圖7所示的仿真結(jié)果，圖8和9則是本文仿真結(jié)果與文獻(xiàn)[4]計(jì)算值及其水槽試驗(yàn)值的比較，可以看出水平張力載荷和網(wǎng)片穩(wěn)定后的L/H值，在不同水流作用下趨勢(shì)是相同的，誤差在可以接受范圍內(nèi)，仿真結(jié)果比較理想。

圖7 網(wǎng)片B的仿真結(jié)果Fig.7 The simulated result of Net B

圖8 不同流速下水平張力載荷Fig.8 Horizontal tension loading at different flow

圖9 不同流速下L/H的值Fig.9 The value of L/H at different flow

4 拖網(wǎng)仿真模擬

我國(guó)東海水產(chǎn)研究所在漁具模型靜水槽內(nèi)進(jìn)行了中層拖網(wǎng)及底拖網(wǎng)共九項(xiàng)網(wǎng)具的模型試驗(yàn)，并給出了試驗(yàn)結(jié)果。本文選取其中一類型中層拖網(wǎng)ZT8909(52×26 m) 進(jìn)行了仿真研究，參數(shù)如表2所示，將仿真結(jié)果同其水槽試驗(yàn)結(jié)果進(jìn)行了比較。

表2中層拖網(wǎng)參數(shù)
Table 2 Parameters of a midwater trawl

網(wǎng)段編號(hào)Number of net長(zhǎng)度/mLength目數(shù)/個(gè)Quantity直徑/cmDiameter目長(zhǎng)/mLength of net mesh1130.514262130.514263130.514264130.512265130.510266100.510207100.51020

續(xù)表2

網(wǎng)段編號(hào)Number of net長(zhǎng)度/mLength目數(shù)/個(gè)Quantity直徑/cmDiameter目長(zhǎng)/mLength of net mesh87.50.5815950.5810103.20.586.41118563.6129551.8138.410.540.8144.81230.41513.46730.2161212030.11711.2515030.075

圖10是網(wǎng)具阻力隨拖曳速度的變化曲線，阻力大小與拖曳速度基本呈線性增長(zhǎng)，本文仿真值比實(shí)測(cè)值偏小，原因分析如下：本文計(jì)算時(shí)假設(shè)漁具不受洋流的影響，是在理想環(huán)境下計(jì)算的結(jié)果。

表3中網(wǎng)口高度隨拖速的變化可以看出本文精細(xì)積分法求解拖網(wǎng)運(yùn)動(dòng)方程的結(jié)果與試驗(yàn)值最大誤差是6.6%，這說(shuō)明本文數(shù)值方法的精確度較高。圖11是采用本文方法進(jìn)行的拖網(wǎng)仿真，拖速為41.0 cm/s在不同時(shí)刻漁網(wǎng)網(wǎng)口的形狀。圖12是采用本文方法實(shí)現(xiàn)的拖網(wǎng)三維仿真。

圖10 網(wǎng)具阻力隨拖速的變化Fig.10 Drag force of nets with different towing speed

表3 網(wǎng)口高度隨拖速變化值Table 3 Height of trawl mouth with different towing speed

圖11 不同時(shí)刻漁網(wǎng)網(wǎng)口形狀Fig.11 The trawl mouth at different time

圖12 漁網(wǎng)穩(wěn)定后的形狀Fig.12 The net shape after stability

漁船模擬器中對(duì)漁網(wǎng)的仿真不僅需要其物理模型的準(zhǔn)確性，更需要保證整個(gè)視景系統(tǒng)的實(shí)時(shí)性。前者也即仿真的效果和穩(wěn)定性，后者則對(duì)應(yīng)于仿真的速度，這兩者之間的矛盾在計(jì)算機(jī)仿真中必須仔細(xì)衡量。圖13是采用不同數(shù)值解法進(jìn)行網(wǎng)衣數(shù)學(xué)模型解算一次的耗時(shí)時(shí)間比，網(wǎng)衣質(zhì)量點(diǎn)數(shù)是21 720時(shí)，六級(jí)五階龍格庫(kù)塔法解算時(shí)間為95 ms，Newmark法耗時(shí)68 ms，四階龍格庫(kù)塔法耗時(shí)56 ms，而精細(xì)積分法耗時(shí)30 ms；可以看出在這幾種方法中，精細(xì)積分法的模型解算耗時(shí)是最少的。

圖13 不同數(shù)值解法的解算時(shí)間Fig.13 Solution time of different numerical methods

5 結(jié)語(yǔ)

本文在漁網(wǎng)網(wǎng)衣建模理論和數(shù)值解法的研究基礎(chǔ)上，基于集中質(zhì)量法建立了漁網(wǎng)網(wǎng)衣的運(yùn)動(dòng)數(shù)學(xué)模型，采用精細(xì)積分的方法對(duì)數(shù)學(xué)模型進(jìn)行求解，通過(guò)和相關(guān)文獻(xiàn)中的水槽試驗(yàn)數(shù)據(jù)的對(duì)比證明該數(shù)值解法的準(zhǔn)確性。整個(gè)系統(tǒng)的仿真速率也得到保證，從而為漁船模擬器中拖網(wǎng)模型建立和視景可視化研究打下了基礎(chǔ)。