大網目底拖網在穩定流場中的數值模擬

王永進，萬 榮，張 勛，張 禹，劉龍騰，王魯民

（1.中國水產科學研究院東海水產研究所，上海 20090；2.上海海洋大學，上海 201306；3.中國水產科學研究院，北京 100141）

底拖網是一種復雜的錐形網具，為了設計或優化拖網漁具，研究方法主要包括模型試驗、海上測試和數值模擬，來分析和評估在實際作業中現有漁具（或新設計漁具）的水動力性能［1］。但是，漁具模型試驗或海試需要花費較多的人力和成本，而計算機技術的進步推動了模擬柔性漁具數值方法的發展。

柔性漁具數值計算方法之一是將網目看作一組集中質量的點和彈簧桿件，質量點代表節點，桿件代表目腳。LEE等［2-3］首先采用集中質量法建立網片、拖網和圍網網具系統的數學模型，計算水下柔性網具的形狀和運動狀態；BESSONEAU等［4］建立了網衣的力學平衡方程進行求解，開展了水流中拖網變形和受力情況的分析；還有學者對網狀和實體錐形結構［5］以及特定的網片［6］等進行了仿真分析，以便能通過可視化技術呈現網具的形狀和受力情況；萬榮等［7-10］將其應用于穩流中網片、中層拖網以及網箱的形狀變化和受力情況研究，得到了良好的模擬效果，并發展了適用于網箱網目群化計算以及漁具模型空間形狀的計測方法；陳英龍等［11］構建大型中層拖網模型，并利用三維仿真呈現出瞬時作業狀態，分析了拖網網具的張力分布。另一種方法是通過將網片離散成三角形單元進行網具形狀計算［12］，克服了目腳單元的目腳必須與網線平行的缺陷，通過此方法，DANIEL PRIOUR［13］以三角形單元離散法計算了網口阻力；在這些研究基礎上，AMELIA等［14］以1D桿單元模擬繩索，2D三角形單元模擬網片，比較分析了底拖網漁具在靜態平衡時形狀的兩種數值方法 牛頓-拉斐爾迭代法和動態仿真法，研究表明后者計算結果較穩定但計算過程較慢。

根據以上研究，雖然數值模擬在計算網具形狀和受力等方面取得了一定的效果，但在建模過程中，由于拖網網具系統規模龐大，結構復雜，網目繁多，或因計算機能力的限制，如果在有限元分析中單純以實際單元進行數值計算，將導致計算量巨大而難以實現。為此，BESSONEAU等［4］在拖網數值計算中提出了網目群化的概念，即根據實際需要可以用1個虛擬的計算網目代替一定數量的實際網目，從而減少拖網的計算單元數。IGOR等［15］也提出了等效桿單元的概念，試圖以有限個等效桿單元替代部分網衣單元，以提高計算效率。萬榮等［10］提出了一種適用于網箱耐流特性有限元分析的網目群化方法，驗證了其可行性，并討論了不同流速和配重條件下網目群化數目對計算結果和計算效率的影響。DANIEL PRIOUR等［12］提出的三角形有限單元離散網片實質也是一種網目群化的概念。網目群化是在網箱、拖網等漁具數值模擬上提出并應用，目前，除萬榮等［10］較詳細介紹了適用網箱的群化方法，其他學者均未對網目群化提出具體的描述。另外，萬榮等的方法并不完全適用于拖網、張網等更為復雜的網漁具。

大網目拖網網目尺寸變化快，網身部位分段多，手工編織的網片縱向目數少且存在較多的掛目或合并目。本文采用集中質量法，將網衣網目離散為一系列在兩端節點處通過無摩擦鉸接而成的柔性集合體，建立大網目底拖網的數值模型，模擬其在穩流狀態下受力及形狀變化。在保證準確性的前提下，為了減少計算量以提高計算效率，提出一種適用于拖網漁具數值計算的網目群化方法。同時，利用此群化方法，對同一頂拖網進行不同網目數量的群化，并對比水槽試驗結果，驗證此方法的實用性以及拖網漁具數值模擬的準確性。

1 材料與方法

1.1 物理模型

所選母型網為中國東海廣泛使用的大網目底拖網，屬于單船有翼單囊拖網。物理模型根據田內準則設計制作，大尺度比設為1∶35，模型圖見圖1。物理模型網全網網線直徑均為0.5 mm。網口部分及袖網網目尺寸為514mm，其他網身各段網目尺寸為450 mm到10 mm不等。上綱長度為6 214 mm，下綱長度為6 271mm。沉子綱空氣中重175 g，總浮力為126 g，有8個14 g浮球，2個7 g的浮球。網身部位網衣具體規格見表1。

圖1 拖網物理模型網衣展開圖Fig.1 Unfolded draw ing of traw lmodel

表1 拖網物理模型網身部位網衣規格Tab.1 Specifications of body netting of traw lmodel

1.2 儀器與設備

物理模型試驗在中國水產科學研究院東海水產研究所的靜水槽中實施。靜水槽規格為長90 m×寬6 m×深3 m。試驗測試時的水深為2.81 m，水槽底部為光滑的陶瓷磚。

水槽拖車車速范圍為0～4.0 m·s-1，配有微處理機調速系統，相對精度為P≤1%；光電測速儀精度為±0.01%；空綱前端連接Lu-A型測力傳感器（日本），進行測力，量程為100.00 N，數據單位為g，測力儀器的線性誤差小于滿量程的0.05%；小型魚探儀StructureScan（墨西哥）用于測量網口的垂直擴張，距離分辨率為30.00 mm［16］。

1.3 數值方法

1.3.1 基本假設

本文將拖網網目單元作為空間索桿單元進行處理，在對網目單元解決幾何非線性平衡方程時，對網線材料、變形特性等方面作如下基本假定：

1）網線絕對柔軟，且只能承受拉伸應力；

2）繩索和網線材料各向同性，并且符合應力與應變的彈性虎克定律，且屬于彈性范圍內的小應變；

3）單元無摩擦鉸合，單元張力僅作用于單元的軸線上，并且在單元的整個橫截面上均勻分布且無衰減；

4）外載荷作用力視作作用在桿件的節點上，即桿件單元不發生彎曲，桿件內僅有軸向內力。

底拖網在作業時，其下綱與海底接觸，下綱與海底之間會產生摩擦力，同時由于邊界層的存在，使下綱及下部網衣處的流場會產生變化。本研究中，對“海底”—— 即靜水槽槽底作如下假設：

1）水槽底部是等深平整的，底部摩擦系數相同；

2）忽略邊界層產生的流場變化對拖網形狀和張力的影響；

3）流場中沒有波浪的影響，拖網是在定常流中運動。

1.3.2 基本方程

當網具在穩流中運動時，構成網衣的單元會產生相對位移，從初始位置到形狀變化后的狀態。對于整個網具系統，其總的勢能可描述為［7］：

式（1）中，П是網具系統總勢能，Fi是作用在第i個節點上的等效節點力，Di是單元i端的節點位移，Tg是第g個單元的軸向力，Lg0是第g個單元的初始長度，Lg是第g個單元變形后的長度，Ag是第g個單元或目腳的橫截面積，E是材料的楊氏模量，f是節點自由度，m是單元數量。

由于張力或在大變形中單元有限位移，第g個單元的伸長量［Lg（Di）-Lg0］可描述為：

式（2）中，x、y、z是單元分別在x、y、z方向的節點坐標。x軸和z軸分別表示水流方向和水深方向的反向。y軸垂直于xz平面。變量u、v和w表示單元在x、y和z方向的節點位移。下標i和j表示單元的兩端。Di是單元i端的節點位移，Lg是第g個單元變形后的長度，Lg0是第g個單元的初始長度。

根據最小勢能原理，當這個離散系統處于平衡狀態時，也就是整個拖網系統在水流作用下處于平衡狀態時，系統的勢能最小。對于具有兩個獨立變量的能量泛函數表達式，其最小值就是泛函數的駐點。由式（1）對Tg求變分，可以得到張力的駐值表達式：

根據上式，各個單元的應力-應變關系表達式為：

式中，П是網具系統總勢能，Tg是第g個單元的軸向力，m是單元數量，E是材料的楊氏模量，f是節點自由度，Lm是第m個單元變形后的長度，Df是單元f端的節點位移，Ff是作用在單元f端等效節點力，Ag是第g個單元或目腳的橫截面積，Lg是第g個單元變形后的長度，Lg0是第g個單元的初始長度，Di是單元i端的節點位移。

根據作用在目腳上的水動力公式：

式（5）中，R是作用在目腳上的水動力，CR是水動力系數，d是單元直徑，L是單元長度，V是水流的相對速率，ρ是流體密度。

當l、m和n分別表示變形后單元對于x、y和z軸的方向余弦，則水動力系數CR（Cx、Cy和Cz為3個方向的水動力系數）表示為：

式（6）中，C0為單元與水流垂直時的水動力系數，通常來講，水動力系數取決于根據單元直徑計算出來的雷諾數（Re），在一般的網箱和拖網研究中，水動力系數一般為2×103～6×103，故取C0為1.3［17］。

本文中摩擦力計算公式為：

式（7）中，fm為底綱與水槽底部的摩擦力，μ為底綱（鐵鏈）與水槽底部（瓷磚）的摩擦系數，其值約為0.1，FN為漁具模型對水槽底部的正壓力。

1.3.3 網目群化方法

為了精確地模擬拖網的形狀和受力，考慮作用于拖網網衣及綱索所有單元的水動力是非常必要的。但由于拖網系統規模大，導致其計算量巨大或耗費較長時間或難以實現。為此，有的學者引入“網目群化”概念，并提出群化的原則為［3，10］：

1）保持群化前后網片投影面積相等。在保證幾何尺度相似時，假若網形相似，則網片與水流的沖角相似，則需要有相同的線面積（將結節分散到各個單腳）；

2）保持群化前后網片的水動力相等。網片的水動力與網片的線面積、網片與水流的沖角相關，若沖角和線面積相等，則網片受到的水動力即相等；

3）保持群化前后網衣的重量相等。

水動力基本公式為［18］：

式（8）中，R為網片受到的水動力，ρ為流體密度，CDN為水動力系數，ST為網片的投影面積，V為流體流動速度。

根據群化前后單元所承受的水動力相等的原則，在水動力計算時，ρ、CDN、ST、V在網目群化前后保持一致。因為需要群化前后網具的形狀基本保持一致，則群化前后網衣與水流的沖角基本相似，則線面積保持一致（不考慮結節，或將結節平均分到相連的各個單腳）。

根據線面積計算公式：

式（9）中，S為網片的線面積，N為網目數量，a為目腳長度，d為目腳直徑。則有：

式（10）中，d1、N1為群化之前的目腳直徑和網目數量，d2、N2為群化之后的目腳直徑和網目數量，n為群化目數，即1個群化的虛擬計算網目的目腳所包含的真實網目的網目腳數量，K為群化系數。

根據群化前后所有計算的單元的總重量相等原則：

從而得出：

式（11）～（13）中，V1、V2為網目群化前、后網衣的體積，ρ1、ρ2為群化前、后的網線密度，a1、d1、N1為群化之前的目腳長度、目腳直徑和網目數量，a2、d2、N2為群化之后的目腳長度、目腳直徑和網目數量，K為群化系數。

由于在網目群化中，網衣的重量不變，其重力即不變，而其體積發生了變化，必然導致網衣受到的浮力產生變化，這在研究文獻中并沒有提及，若未考慮則會導致群化后的網衣所受到的浮力和重力產生較大的影響，因此本文提出在計算網具系統的水動力時，流體的密度不變，而在計算群化網目浮力時，將流體的密度進行修正，使得網目群化前后網衣所受到的浮力和重力的合力相等，由此有公式：

式（14）、式（15）中，Bf1、Bf2分別為網目群化前、后的浮力與重力的合力，ρ水1、ρ水2分別為群化前、后的流體密度，V1、V2為網目群化前、后網衣的體積，g為重力加速度。

由于體積不相等，所以得出：

或者為：

式（16）、式（17）中，ρ水1、ρ水2分別為群化前、后計算網衣浮力的水體密度，V1、V2分別為群化前、后網衣體積，K為群化系數，ρ1、ρ2為群化前、后的網線密度。

由于此網的特點是網目較大，網口及網袖處的網目數較少，為了保持網口部分結構相似，本文只將網身部分進行網目群化，并且群化為兩種規格，model-1為單元數為6 513的拖網數學模型，model-2為單元數為1 787的拖網數學模型，兩種群化后的拖網數學模型結構規格如表2、表3所示。

表2 model-1拖網數學模型結構參數Tab.2 Structure parameters ofmathematicalmodel ofmodel-1 traw l

表3 model-2拖網數學模型結構參數Tab.3 Structure parameters ofmathematicalmodel ofmodel-2 traw l

1.3.4 數學模型初始位置

拖網是一種由網片、繩索、浮球以及沉子等構成的大規格漁具。為了建立數學模型，將整個拖網網目及繩索離散成有限單元，每個單元包括一個彈性桿件和兩個集中質量的點。AMELIA等［14］認為初始位置不必是準確的，但需要給予一定的預加張力。本文將兩種規格群化后的拖網初始位置設置為如圖2和圖3。預加張力為0.1 g。程序編碼在MATLAB平臺上完成。

圖2 m odel-1拖網數學模型初始形狀Fig.2 Initially assumed geometry ofmodel-1 traw l net

圖3 model-2拖網數學模型初始形狀Fig.3 Initially assum ed geometry ofmodel-2 traw l net

2 結果與分析

2.1 流場中拖網形態變化

圖4為流速分別在42.9、51.4、60.0、68.6、77.2 cm·s-1時model-1拖網模型展開形態的計算機模擬結果。如圖4所示，當拖網受到水流作用時，網口部位從初始位置整體向后方移動，上、下綱和空綱向中間偏移，上中綱向上鼓起，當拖網受力均衡之后，上、下綱均成懸鏈線狀，網衣由初始設置的矩形筒變成圓筒狀。袖端及網身網目展開良好，網身后部及網囊部位網衣沒有充分展開，網型收縮較快。隨著流速的顯著增大，網口高度大幅降低，網目張開也越小，網型整體趨于流線型。

圖5分別為流速在42.9、51.4、60.0、68.6、77.2 cm·s-1時，model-2拖網模型展開形態的計算機模擬結果。如圖5所示，拖網在水流作用下，網口網位從設定的初始位置向后方移動，上、下綱和空綱向內偏移，上中綱向上鼓起，下空綱和下綱Z軸方向沒有位移。拖網受力均衡時，上、下綱均成懸鏈線狀，網衣由初始設置的矩形筒變成圓筒狀。袖端及網身背部網目展開良好，網身腹部網衣略向上浮起，網囊末端有上翹現象，這可能由于群化后網囊末端網目數極少（僅1目），網目不能很好的展開，受浮力略大的影響而產生上翹。隨著流速的顯著增大，網口高度大幅降低，網目張開也越小。model-2模型網相較于model-1模型網，網目數少，網目尺寸大，從兩者流場中形態的模擬結果看，在相同工況下，兩者的形態變化趨勢基本相似，特別是網袖及網身部位網衣展開及其形態變化一致。

圖4 不同流速下model-1拖網模型展開形態Fig.4 Simulated model-1 traw l shape w ith different current velocities

2.2 試驗與計算網具阻力比較

圖6為水槽模型試驗（model）和計算機模擬（model-1和model-2）阻力結果的比較。由圖6中可以看出，隨著流速由42.9 cm·s-1增加到77.2 cm·s-1，拖網受到的阻力顯著增加，而不同速度下，模擬結果均比物理模型試驗結果大，且速度越大，計算機模擬的阻力值越大。從圖6中也可以看出，群化單元較多的mode-1的數學拖網模型的計算結果與物理模型model的結果非常相近，只是在較大流速下，兩者略有差異，群化單元較少的model-2的數學拖網模型的計算結果略大，不同速度下兩者相差16%～26%。結果說明，群化單元越多的數學模型的網具阻力計算值與物理模型更相近，群化單元由6 513個減少到1 787個，將會產生20%左右的誤差。

圖5 不同流速下model-2拖網模型展開形態Fig.5 Simulated model-2 traw l shape w ith different current velocities

2.3 試驗與計算網口垂直擴張比較

圖7為水槽模型試驗（model）和計算機模擬（model-1和model-2）拖網網口垂直擴張結果的比較。由圖7中可以看出，隨著流速由42.9 cm·s-1增加到77.2 cm·s-1，網口垂直擴張具有明顯的下降趨勢，而3者之間的差距不大。當流速較低時，物理模型的網口垂直擴張與群化單元較多的model-1的數值模擬結果相近，而大于群化單元較少的model-2數學模型計算結果；當流速較高時，兩種數學模型的模擬結果基本一致，略大于物理模型試驗結果。整體上數學模型計算結果與物理模型結果相差5%左右，速度最大時相對誤差也僅在10%左右，結果符合良好。

圖6 試驗與計算的網具阻力比較Fig.6 Com parison of resistance between experiment and calculation

圖7 試驗與計算的網口垂直擴張比較Fig.7 Com parison of vertical expansion of net mouth between experiment and calculation

2.4 試驗與計算袖端水平擴張比較

將模型試驗拖網下袖端距離與計算機模擬結果進行比較（圖8），圖8中實線表示模型試驗結果。在拖曳過程中，為了保證水平擴張比為0.35，在試驗之前，對網板（拖車牽引桿）間距進行調試，并固定值為3.4 m，通過觀測，下袖端間距幾乎不變，其距離為2.19 m左右。圖中虛線分別為兩種數學拖網模型的計算結果，為了與試驗值相比較，將數學模型網板間距設置為3.4 m，以保證相同的試驗工況。由圖8可以看出，兩種數學模型的計算結果與模型試驗結果相差不大，隨著拖速的變化，袖端水平擴張均略有浮動，浮動程度較小。由圖8可以看出，不同拖速下，計算值均較試驗值小，但兩者相差僅5%左右，說明數值模擬結果與水槽模型試驗結果符合較好。

圖8 試驗與計算的袖端水平擴張比較Fig.8 Comparison of horizontal expansion of the end of net w ings between experiment and calculation

3 討論

本文針對大網目底拖網建立了數學模型，提出了適用于拖網類漁具的網目群化方法。萬榮等［10］和LIU等［19］認為在網目群化前后，需要網衣水中重量和投影面積相等以保證網箱網衣總體水動力相同。然而，材料密度的變化（文中分析密度是減小的）必然導致網衣體積增大，浮力增大，這在萬榮等［10］的研究中并未考慮，原因可能是網箱下緣的重力足夠大從而可以對浮力的增量忽略不計。但是，在拖網、張網等此類漁網具中，網身部分沒有配重，需要對體積變化加以考慮，否則可能造成網衣整體浮力大于重力，出現網體上翹現象。

兩種數值模型的網具阻力計算結果有一定差距，model-1計算值接近試驗值，而model-2數值模擬的網具阻力較模型試驗大16%～26%，但這與QUEIROLO等［20］發現數值模擬和物理模型試驗結果之間存在差異（13%～23%）的結果相近，可能與群化后的數值模型網衣單元較少，虛擬單元長度較長有一定關系。兩種群化模型在穩定流場中的形態變化趨勢相近，特別是網口及網袖部分，網口垂直擴張和袖端間距結果相差較小，這說明大網目底拖網網身網目群化的程度對計算結果的影響可控，可依據試驗條件適當減少網身部位的網衣單元，在保證一定精確率下提高計算效率。

在數值模型構建過程中，虛擬目腳單元最大長度達到514mm，同時選用的拖網網身網衣具有較多的工藝性掛目或并目，都在一定程度上對計算結果造成了影響，網衣在局部未達到理想的流線型狀態，在今后的優化計算中需要加以分析和改進，以提高計算模擬的精度，更好的呈現拖網在流場中的狀態。