水流作用下的漁網水動力數值方法研究

胡 克，耿寶磊，趙 旭，沈文君

(交通運輸部天津水運工程科學研究所 港口水工建筑技術國家工程實驗室 工程泥沙交通行業(yè)重點實驗室, 天津 300456)

*通訊作者：耿寶磊(1980-), 男，副研究員，主要從事港口及海洋工程水動力研究。E-mail：tksgengbaolei@163.com

隨著人口的急劇增長，對于食品的需求也在不斷增加。海洋面積約占地球面積的70%，目前，人類食品的主要來源仍來自陸地，海洋所提供的食品來源比重很小，未來海洋養(yǎng)殖業(yè)將會成為人類食品的主要來源之一。因此，網箱在波流載荷下的水動力問題將會成為未來研究的重點。

網箱通常放置在有水流的地方，以便改善網箱內的水質，同時進行網箱內氧氣的更新，所以流速對于網衣的水動力和容積影響是研究人員關注的重點。Aarnses[1]通過對重力式網箱進行的拖曳試驗研究了作用于網箱上的阻力以及網衣的變形和箱體內的流速衰減規(guī)律。目前，網衣數值結果比較多的試驗是Lader[2]進行的三維網衣在不同流速和不同配重下的水動力試驗，通過研究發(fā)現(xiàn)，網衣的受力與網衣的變形相互影響；網衣的受力受到雷諾數的影響很大；在計算柔性網衣的受力時，如果使用基于彈性很小的單片網衣得到的阻力公式計算，將會產生很大的誤差。此外，還通過試驗測量了不同流速和配重下的網衣升力和阻力，并擬合出升力和阻力相對雷諾數的經驗公式。除了試驗研究外，Schlichting[3]最早將尾流場的概念描述成“尾流遠場的平均速度損失”。 L?land[4]通過研究發(fā)現(xiàn)，水流經過網平面后，在網衣的寬度范圍內是呈均勻分布的，而且衰減后的流速并不隨著遠離網衣的距離增大而明顯減小。此外，他還根據研究成果推導出了速度衰減的經驗公式[5]。

與大尺度結構波浪力計算方法不同[6]，網箱與牧場周圍圍欄設施[7]等屬于小尺度結構物。計算網衣水動力的模型主要有兩種：Morison方程模型和網平面模型。對于Morison方程模型來說，是將每根網線看作圓柱體，根據圓柱的雷諾數和KC數來計算網衣的水動力系數，進而根據Morison方程求解水動力。對于網平面模型，是先將網衣按照不同密實度(Solidity Ratio)和與來流攻角測量出網衣的法向和切向力系數，再將網箱上每個網目按照響應的法向和切向力系數計算出來。

目前，使用Morison方程計算網衣的方法較多，例如Huang[8]、趙云鵬[9-10]和許條建[11]等都使用了這一方法計算波流作用下的重力式網箱。在他們的水動力計算模型中都是在每個迭代步后根據雷諾數更新水動力系數的，在對網衣結構的模擬是基于集中質量點法的，這種方法假定網衣的質量集中在節(jié)點處，相鄰節(jié)點之間使用彈簧連接并模擬剛度。另一方面，Moe[12]和Li[13]則使用有限元方法來模擬網衣合股線結構，但是水動力模型進行了簡化，即根據流速和網衣合股線直徑計算出雷諾數再選取固定的阻力系數計算網衣的水動力。除了Morison方程以外， Kristiansen[14]則基于網平面模型用于計算網衣的水動力。在他對網衣研究成果中，首先根據Rudi[15]進行的單片網衣試驗結果，使用傅里葉級數的方式擬合出升力系數和阻力系數關于網衣密實度和攻角的經驗公式，再計算出網衣的水動力。在他的網衣模型中也是集中質量點法來模擬網衣的結構。本文嘗試使用一種新的方法：根據網衣的物理特性，基于有限元方法對網衣結構進行離散，建立相應的系統(tǒng)控制方程，計算求解網衣的動力響應，并基于文獻[14]中的水動力模型，準確求解出網衣的水動力。

為了實現(xiàn)這一方法，本文基于ABAQUS用戶子程序的二次開發(fā)技術，采用FORTRAN語言，編寫了適合于網衣的水動力計算子程序。根據這一方法，先通過每個迭代步輸出的節(jié)點信息(包含速度和加速度)計算出網目雷諾數以及與來流的夾角，然后根據試驗所得的單片網衣阻力和升力系數，計算出每個網目的阻力和升力，進而得到各個時刻下的整個漁網的水動力。

1 基本理論介紹

1.1 運動方程

當網箱處于波流作用下，其結構的動力平衡方程可以表示成為

(1)

Fext=G+fb+fw+fc

(2)

式中：[m]為質量矩陣，[c]為阻尼矩陣，[k]為剛度矩陣。結構的外載荷用Fext表示，如式(2)所示，共包含四個部分：重力G，浮力為fb，波浪力為fw，水流阻力為fc。

式(1)和式(2)中，重力可以由靜態(tài)浮力計算得到，浮力fb由結構的瞬時浸沒的浮圈結構計算得出。

1.2 水動力模型

對于網衣結構，本文使用的是網平面模型計算。每個網目上的升力和阻力都是按照如下公式進行計算的

(3)

(4)

Rudi等人[15]首先通過試驗方法對平面網衣平面下不同來流攻角下，不同密實度的網衣進行阻力系數和升力系數測量，在試驗中的網衣密實度Sn包括0.130,0.144,0.184,0.243和0.317，試驗流速包括0.159，0.316和0.966 m/s。Kristiansen[14]對這次試驗中的阻力和升力系數使用正弦行數曲線在不同攻角(θ)和網衣密實度(Sn)的網衣平面的阻力系數和升力系數進行擬合。得到的結果如下

CD=cd(α1cosθ+α3cos3θ)

(5)

CL=cL(b2sinθ+b4sin3θ)

(6)

在式(5)和式(6)中，α1和α3分別等于0.9和0.1，b2和b4分別等于1和0.1。cd、cL則參照文獻[14]中的計算方法進行計算。

1.3 幾何非線性理論

基于小變形假設，建立平衡條件時不必考慮物體位形的變化，因而應變可以采用位移的一階偏導進行度量。在本文所處理的網箱變形中由于大撓度和大轉動的存在，使其不滿足小變形假設，在有限元求解中必須考慮幾何非線性的影響。

1.3.1 應變-位移之間的關系

當結構需要考慮幾何非線性時，應力與應變的關系可以寫成

(7)

(8)

式中：[B0]是線性應變項，與位移{δ}e呈線性關系；[BL]是大位移應變項，是位移{δ}e的函數。

1.3.2 應力-位移之間的關系

單元內部應力增量和應變增量關系可以表示成為：

d{σ}e=[D]d{ε}e

(9)

式中：[D]為材料的本構關系矩陣。聯(lián)合式(4)、(5)和(6)可以得到如下形式

(10)

式中：{ε*}e和{σ}e分別代表單元的虛應變和虛應力，{δ*}e為單元節(jié)點虛位移，{F}e為單元節(jié)點力。

聯(lián)合式(7)和(10)可以得到

(11)

再對式(11)進行微分后可得

(12)

式(12)又可以寫成

([k0]+[kσ]+[kL])d{δ}e=d{F}e

(13)

式(13)為幾何非線性問題的求解方程。其中，[k0]是小位移線性剛度陣；[kL]為大位移非線性剛度陣；[kσ]是初始應力非線性剛度陣；以上三個剛度陣分別表示為

(14)

(15)

(16)

1.4 等效網目計算方法介紹

網箱的網衣在計算水動力時，由于實尺度網衣上的單元較多，因此很難將計算模型的單元數與實尺度的網衣單元保持一致，事實上，網衣計算模型需要簡化。對于簡化后的模型，應該保持質量和剛度與簡化前一致，如式(17)～(18)所示

Mtruss=M

(17)

(EAsection)truss=ΣEkAsection＿k

(18)

式中：A和Asection代表網衣合股線的投影面積和橫截面積，M代表網衣的質量，E代表網衣的彈性模量。蘇煒[16]和Moe[12]曾驗證過使用等效網格將網目合并后計算網衣模型與試驗模型中網衣的變形，并得到很好的計算結果，本文將基于這一方法簡化網衣模型。

圖1 浮筒分段浮力計算示意圖Fig.1 Buoyancy calculation of floating collar′s segment

1.5 網箱體積計算方法介紹

網箱的體積變化直接影響到養(yǎng)殖物的生存空間，因此，對于網箱容積的變化是網箱水動力研究的主要部分，為了對網箱容積變化進行研究，首先將網箱在波流作用后的容積與靜止時的網箱比值定義為網箱的容積變化率，如式(19)所示

(19)

式中：Cv為水流條件下網箱的容積變化率，Vc為波流作用下的網箱容積，V0為網箱在靜止狀態(tài)下的體積。

對于網箱的容積計算，本文使用的是如下計算方法：首先，將網箱沿中軸線等分為若干楔形體，而每個楔形體又近似看做為棱柱結構，將每個棱柱結構可以分為3個四面體結構，如圖1所示。

圖2 四面體示意圖Fig.2 Sketch of tetrahedron

對于任意四面體OABC的體積計算方法(如圖2所示)，其計算公式可以寫為

(20)

按照式(20)將圖1中的三個四面體的體積計算出來以后再相加，就得到了一個棱柱的體積，最后再將網箱上每個棱柱的體積相加，就得到了此時網箱的容積。

1.6 流速衰減計算方法介紹

水流經過網衣平面后會出現(xiàn)流速減小的情況。L?land[5]曾經對流速衰減進行了研究，他將網箱分解成若干網片，并對流速通過網衣后的阻力進行了總結，得到了不同網衣密實度Sn下的網衣的流速衰減系數r和阻力系數之間的關系式，得出公式(21)的表達式

r=1.0-0.46Cd

(21)

圖3 網衣水動力施加方法示意圖Fig.3 Sketch of the method of hydrodynamic force acting on the twine

2 水動力施加方法介紹

本文使用的方法水動力模型是基于試驗結果得到的，其水動力系數計算方法已經在1.2節(jié)中進行了闡述，在這里只介紹在ABAQUS中的網衣水動力施加方法。

3 計算結果分析與討論

3.1 網衣疏密變化的比較

根據Moe[12]的介紹，計算的網衣模型的整體直徑為1.435 m，深度為1.44 m，網目的邊長為1.76 cm，合股線直徑為2 mm，網衣的密實度是按照如下公式進行計算的

(22)

在式(22)中，d代表網衣合股線直徑，λ表示網目邊長。計算得出Sn=0.214，此外，根據Moe[12]提供的參數，網衣的剛度EI=8.2×107Pa。在網衣的最上端為剛性圓圈，且由均勻分布四個簡支的支點固定。沉子為16個800 g、600 g和400 g的重塊，根據Moe[12]的介紹，重塊在水中的重量分別為6 N、4.5 N和3 N，在本文的計算模型中將以集中力的方式施加到網衣底部。為了研究不同網格數對于網箱水動力計算結果的影響，在本文中使用了兩種網格數的網衣，分別為64×21和32×12。與Lader的試驗模型相同，本文建立的網衣模型也是沒有底部網衣的，兩種網格的網衣模型如圖4所示。

圖5給出了兩種數量的網格模型計算結果。在本文計算的流速中，共進行了三個流速：0.2 m/s，0.34 m/s和0.5 m/s。從總體上看，兩種網格的計算結果在不同流速下較為相近，在流速為0.5 m/s時，計算值小于試驗值，此時的相對誤差為9.5%，這可能是由于此時的網衣變形最大，所以相對誤差要大于流速較小的情況。對于升力，網箱的升力是通過網箱頂部固定的四個節(jié)點在垂直方向的合力減去沉子的重力后得出的。從趨勢來看，數值模擬的結果與網箱的實驗結果都是隨著流速的增大而增大的。在流速為0.5 m/s時計算值與實驗結果的相對誤差較大，為10.8%。

4-a 64×214-b 32×12圖4 網格對比Fig.4 Comparison of different grids of the net cage圖5 網箱試驗和計算的阻力和升力結果對比圖(16×800 g)Fig.5 Comparisons of drag force and lift force between model tests and calculated result in numerical analysis(16×800 g)

圖6中顯示的是在流速為0.34 m/s情況下，通過不同網格數下網箱的變形情況與試驗結果的比較可知，疏網格(32×12)與密網格(62×21)的變形大體上是一致的，同時結合阻力與升力的計算結果，可以使用疏網格計算網衣的阻力和變形情況。此外，為了節(jié)省計算時間，下面將使用疏網格(32×12)對網箱的水動力進行計算。

6-a 64×216-b 32×126-c 模型試驗[5]圖6 在流速為0.34 m/s時不同網格的數值模擬變形與模型試驗結果對比Fig.6 Comparison of deformed shape from numerical models with detailed and coarse mesh and physical model(current velocity of 0.34 m/s)

3.2 不同沉子重量下網箱的水動力

為了研究沉子的重量對于網箱的水動力的影響，本文對比Lader[2]的試驗結果，計算不同沉子質量下的網箱的水動力。在3.1節(jié)計算了16×800 g沉子的情況，在這里顯示的是其余兩個配重的計算結果。

對比圖6～8，可以看出，計算得到的阻力和升力隨流速的變化趨勢與試驗結果是相同的，都是隨著流速的增大而增大的。同時通過對比不同沉子下的阻力和升力的計算值與試驗值結果，發(fā)現(xiàn)兩者較為吻合。此外，隨著流速的增大，計算結果相對實驗的誤差也是在逐漸增大的。在本文計算的幾個工況中，最大的阻力和升力誤差為沉子16×400 g且流速為0.56 m/s的情況(如圖8所示)，其與試驗結果的相對誤差分別為13.1%和10.7%。Kristiansen[2]曾對該試驗進行了誤差分析：由于漁網具有吸水的特性，使得試驗時的漁網直徑相對實際直徑較小。同時，該試驗中的沉子在測量水中的質量時也存在一定的誤差。所以，由于試驗存在一定不確定性因素可能使得計算結果較試驗結果存在一定的誤差。

3.3 不同流速下網箱的容積變化規(guī)律研究

為了研究流速對于網箱容積的影響，本文首先根據靜止狀態(tài)下的網衣各節(jié)點坐標，利用公式(20)計算出網箱初始容積。為了更加方便的導出計算結果，本文基于Python語言，編寫了可以直接將ABAQUS結果導出的命令流，使用這一方法將穩(wěn)定時刻網衣上各個節(jié)點的位移導出，再加上網衣初始時刻的節(jié)點坐標，即可得到網箱各個節(jié)點的位置。在水流作用下由于沒有相關的網箱體積的相關試驗結果，本文對比的是Moe[12]的計算結果，如圖9所示。

圖7 網箱試驗和計算的阻力和升力結果對比圖(16×600 g)Fig.7 Comparisons of drag force and lift force between model tests and calculated result in numerical analysis(16×600 g)圖8 網箱試驗和計算的阻力和升力結果對比圖(16×400 g)Fig.8 Comparisons of drag force and lift force between model tests and calculated result in numerical analysis(16×400 g)圖9 不同沉子重量下的網箱相對容積隨流速變化情況Fig.9 Relative volume in net cage under different weight sinker corresponding to current velocity

在圖9中，本文計算的網箱是沒有網箱底部的，所以在計算時只計算了網箱網衣圍住的部分，這一點與試驗[5]和Moe[12]的計算模型是一致的。網箱的初始體積是根據網箱上各個點的初始坐標計算出來，使用的是32×12的網格。可以看出，隨著流速的增大，網箱的體積是呈遞減的趨勢。當流速增大到0.5 m/s時，網箱的體積減少了約50%。此外，還可以明顯地看到，網箱的容積與流速之間是非線性關系。

4 結論

本文基于ABAQUS的二次開發(fā)技術對網箱在水流作用下的受力進行了數值模擬。在計算中，本文使用子程序編譯了基于網平面方法的水動力模型，并與試驗結果進行了比較。研究結果顯示，網衣的升力和阻力均隨著流速的增大而增大，且計算與試驗結果的最大相對誤差為10.8%。而網箱的容積是隨著流速的增大而減小的，在本文計算的模型中，當流速達到0.5 m/s時網箱的容積減小了約50%。網箱容積的計算結果與Moe[12]的基本一致。所以，通過本文的計算結果與試驗結果對比可以看出，該方法可以應用于的網衣在水流作用下的數值模擬。

目前，該方法主要的研究對象是水流作用下的網箱的受力以及網箱體積的變化情況，未來該方法還可用于計算海洋柔性結構物中需自定義水動力載荷的情況。