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

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

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

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

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

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

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

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

1 基本理論介紹

1.1 運(yùn)動方程

當(dāng)網(wǎng)箱處于波流作用下，其結(jié)構(gòu)的動力平衡方程可以表示成為

(1)

Fext=G+fb+fw+fc

(2)

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

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

1.2 水動力模型

對于網(wǎng)衣結(jié)構(gòu)，本文使用的是網(wǎng)平面模型計算。每個網(wǎng)目上的升力和阻力都是按照如下公式進(jìn)行計算的

(3)

(4)

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

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則參照文獻(xiàn)[14]中的計算方法進(jìn)行計算。

1.3 幾何非線性理論

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

1.3.1 應(yīng)變-位移之間的關(guān)系

當(dāng)結(jié)構(gòu)需要考慮幾何非線性時，應(yīng)力與應(yīng)變的關(guān)系可以寫成

(7)

(8)

式中：[B0]是線性應(yīng)變項(xiàng)，與位移{δ}e呈線性關(guān)系；[BL]是大位移應(yīng)變項(xiàng)，是位移{δ}e的函數(shù)。

1.3.2 應(yīng)力-位移之間的關(guān)系

單元內(nèi)部應(yīng)力增量和應(yīng)變增量關(guān)系可以表示成為：

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

(9)

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

(10)

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

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

(11)

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

(12)

式(12)又可以寫成

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

(13)

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

(14)

(15)

(16)

1.4 等效網(wǎng)目計算方法介紹

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

Mtruss=M

(17)

(EAsection)truss=ΣEkAsection＿k

(18)

式中：A和Asection代表網(wǎng)衣合股線的投影面積和橫截面積，M代表網(wǎng)衣的質(zhì)量，E代表網(wǎng)衣的彈性模量。蘇煒[16]和Moe[12]曾驗(yàn)證過使用等效網(wǎng)格將網(wǎng)目合并后計算網(wǎng)衣模型與試驗(yàn)?zāi)Ｐ椭芯W(wǎng)衣的變形，并得到很好的計算結(jié)果，本文將基于這一方法簡化網(wǎng)衣模型。

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

1.5 網(wǎng)箱體積計算方法介紹

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

(19)

式中：Cv為水流條件下網(wǎng)箱的容積變化率，Vc為波流作用下的網(wǎng)箱容積，V0為網(wǎng)箱在靜止?fàn)顟B(tài)下的體積。

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

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

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

(20)

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

1.6 流速衰減計算方法介紹

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

r=1.0-0.46Cd

(21)

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

2 水動力施加方法介紹

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

3 計算結(jié)果分析與討論

3.1 網(wǎng)衣疏密變化的比較

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

(22)

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

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

4-a 64×214-b 32×12圖4 網(wǎng)格對比Fig.4 Comparison of different grids of the net cage圖5 網(wǎng)箱試驗(yàn)和計算的阻力和升力結(jié)果對比圖(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情況下，通過不同網(wǎng)格數(shù)下網(wǎng)箱的變形情況與試驗(yàn)結(jié)果的比較可知，疏網(wǎng)格(32×12)與密網(wǎng)格(62×21)的變形大體上是一致的，同時結(jié)合阻力與升力的計算結(jié)果，可以使用疏網(wǎng)格計算網(wǎng)衣的阻力和變形情況。此外，為了節(jié)省計算時間，下面將使用疏網(wǎng)格(32×12)對網(wǎng)箱的水動力進(jìn)行計算。

6-a 64×216-b 32×126-c 模型試驗(yàn)[5]圖6 在流速為0.34 m/s時不同網(wǎng)格的數(shù)值模擬變形與模型試驗(yàn)結(jié)果對比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 不同沉子重量下網(wǎng)箱的水動力

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

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

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

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

圖7 網(wǎng)箱試驗(yàn)和計算的阻力和升力結(jié)果對比圖(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 網(wǎng)箱試驗(yàn)和計算的阻力和升力結(jié)果對比圖(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 不同沉子重量下的網(wǎng)箱相對容積隨流速變化情況Fig.9 Relative volume in net cage under different weight sinker corresponding to current velocity

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

4 結(jié)論

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

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