高 帥，尹 勇, 孫霄峰

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

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

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

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

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

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

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

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

1.1 結節(jié)的運動方程

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

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

(1)

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

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

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

(2)

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

(3)

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

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

(4)

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

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

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

(5)

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

把以上(2)～(5)式代入(1)式可得結節(jié)i在空間坐標系下的運動方程：

(6)

1.2 目腳的運動方程

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

圖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時會發(fā)生急劇的改變，因此應將Cd看作是雷諾數(shù)Re的函數(shù)[5]。

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

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

2 數(shù)值解法

精細積分法自鐘萬勰提出來后，在結構動力分析中獲得了廣泛應用。在此基礎上進一步研究，使得非線性動力方程的高精度計算成為可能。結構系統(tǒng)非線性動力方程的一般形式為：

(9)

(10)

其中：

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

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

(11)

對式(11)進行數(shù)值離散，并假設時間步長Δt=tk+1-tk，則第k+1步時的響應公式可寫成：

(12)

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

(13)

(14)

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

(15)

(16)

將X+f(t)看作外部激勵，對結節(jié)運動方程來說：

(17)

而對目腳運動方程，其解算在局部坐標系下，F(xiàn)(t)如下式：

(18)

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

(19)

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

(20)

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

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

其中帶“′”的元素表示T1/2的相應元素。

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

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

3 模型驗證

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

圖4 精細積分計算流程圖Fig.4 Precise integration computational flow chart

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

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

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

表1 測試網(wǎng)片B的計算參數(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.

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

圖7 網(wǎng)片B的仿真結果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)仿真模擬

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

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

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

續(xù)表2

網(wǎng)段編號Number of net長度/mLength目數(shù)/個Quantity直徑/cmDiameter目長/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)具阻力隨拖曳速度的變化曲線，阻力大小與拖曳速度基本呈線性增長，本文仿真值比實測值偏小，原因分析如下：本文計算時假設漁具不受洋流的影響，是在理想環(huán)境下計算的結果。

表3中網(wǎng)口高度隨拖速的變化可以看出本文精細積分法求解拖網(wǎng)運動方程的結果與試驗值最大誤差是6.6%，這說明本文數(shù)值方法的精確度較高。圖11是采用本文方法進行的拖網(wǎng)仿真，拖速為41.0 cm/s在不同時刻漁網(wǎng)網(wǎng)口的形狀。圖12是采用本文方法實現(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 不同時刻漁網(wǎng)網(wǎng)口形狀Fig.11 The trawl mouth at different time

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

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

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

5 結語

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