張 為，郭 軍*，扈 喆，張曉瑩，丁 蘭

(1. 集美大學(xué)輪機(jī)工程學(xué)院，福建 廈門 361021；2. 福建省水產(chǎn)研究所，福建省海洋生物增養(yǎng)殖與高值化利用重點(diǎn)實(shí)驗(yàn)室，福建 廈門361013)

漁業(yè)過度捕撈使?jié)O業(yè)資源量下降，單純依靠海洋捕撈早已無法滿足人類對(duì)水產(chǎn)品日益增長的需求，現(xiàn)在迫切需要海水養(yǎng)殖來彌補(bǔ)水產(chǎn)品消費(fèi)市場(chǎng)的不足。2021年，我國水產(chǎn)品總產(chǎn)量達(dá)6 690×104t，養(yǎng)殖產(chǎn)量與捕撈產(chǎn)量之比約為4∶1，網(wǎng)箱養(yǎng)殖為人們提供了豐富的海鮮食品。隨著國內(nèi)外各種大型養(yǎng)殖設(shè)施的建造，養(yǎng)殖業(yè)逐漸向深遠(yuǎn)海方向發(fā)展[1]，養(yǎng)殖的種類有鮑魚[2]、大黃魚[3]和云龍石斑魚[4]等。在養(yǎng)殖過程中，網(wǎng)衣對(duì)海水有阻流作用，這不僅會(huì)影響?zhàn)B殖對(duì)象生長，還會(huì)阻礙水體交換，影響網(wǎng)箱內(nèi)部的養(yǎng)殖環(huán)境，因此網(wǎng)衣內(nèi)部的流場(chǎng)特性對(duì)網(wǎng)箱養(yǎng)殖非常重要。

研究網(wǎng)衣流場(chǎng)特性通?？梢圆捎媚Ｐ驮囼?yàn)和數(shù)值模擬兩種方法。雖然模型試驗(yàn)可以進(jìn)行養(yǎng)殖網(wǎng)箱網(wǎng)衣流場(chǎng)的分析，但成本高，且受試驗(yàn)場(chǎng)地的限制。隨著計(jì)算流體力學(xué)(Computational fluid dynamics，CFD)的不斷發(fā)展，CFD技術(shù)已能較為準(zhǔn)確地分析一些流體力學(xué)問題，不僅能得到養(yǎng)殖網(wǎng)箱整體受力情況，還可以捕捉養(yǎng)殖網(wǎng)箱內(nèi)外流場(chǎng)的細(xì)節(jié)。桂福坤[5]對(duì)深水重力式網(wǎng)箱的水動(dòng)力特性開展了模型試驗(yàn)，給出了重力式、碟形及擬碟形網(wǎng)箱的受力及運(yùn)動(dòng)特性；Patursson O等[6-7]對(duì)平面網(wǎng)衣的流場(chǎng)進(jìn)行了模型試驗(yàn)和數(shù)值模擬，測(cè)量了網(wǎng)衣的阻力系數(shù)及升力系數(shù)，模擬了網(wǎng)衣后面的流場(chǎng)速度分布；趙云鵬等[8-9]借助商業(yè)軟件FLUENT對(duì)單片網(wǎng)衣周圍的流場(chǎng)特性進(jìn)行了二維數(shù)值模擬，結(jié)果表明：FLUENT中的多孔介質(zhì)模型可以用于網(wǎng)衣周圍流場(chǎng)的數(shù)值模擬；劉興[10]采用多孔介質(zhì)模型，分別對(duì)平面網(wǎng)衣的流場(chǎng)進(jìn)行了二維和三維數(shù)值模擬，結(jié)果表明：三維數(shù)值模擬的誤差更小；劉春宏等[11]對(duì)網(wǎng)衣和魚類共同作用下的網(wǎng)箱周圍流場(chǎng)進(jìn)行了數(shù)值模擬，發(fā)現(xiàn)需同時(shí)考慮魚類和網(wǎng)衣才能準(zhǔn)確地模擬養(yǎng)殖網(wǎng)箱周圍的流場(chǎng)。

為研究及分析養(yǎng)殖網(wǎng)衣周圍的流場(chǎng)特性，本文基于CFD理論，采用多孔介質(zhì)模型對(duì)單片網(wǎng)衣的流場(chǎng)進(jìn)行數(shù)值模擬，首先進(jìn)行不確定度分析，再對(duì)網(wǎng)衣在不同流速下的流場(chǎng)進(jìn)行數(shù)值計(jì)算，并與試驗(yàn)結(jié)果進(jìn)行了比較，最后對(duì)網(wǎng)衣在不同攻角下的流場(chǎng)進(jìn)行數(shù)值計(jì)算，分析了網(wǎng)衣的流場(chǎng)特性，給出網(wǎng)衣對(duì)流場(chǎng)內(nèi)速度影響的量化區(qū)域及范圍，旨在能夠較好地為實(shí)際養(yǎng)殖生產(chǎn)活動(dòng)提供參考依據(jù)。

1 理論模型

1.1 控制方程

流體力學(xué)問題要遵守質(zhì)量守恒、動(dòng)量守恒和能量守恒定律，對(duì)于不可壓縮流體，其連續(xù)性方程和動(dòng)量守恒方程分別如下[12]：

(1)

(2)

式(1)～(2)中：ui和uj為速度；xi和xj為坐標(biāo)分量；t為時(shí)間；p為流體的壓力；ρ為流體的密度；μ為流體的黏性系數(shù)。

控制方程不封閉，因而需引入湍流模型，使方程能夠求解。本文計(jì)算采用標(biāo)準(zhǔn)k-ε湍流模型[12]，該模型具有較髙的穩(wěn)定性、經(jīng)濟(jì)性和計(jì)算精度，應(yīng)用廣泛，其表達(dá)形式如下：

(3)

(4)

式(3)～(4)中：k和ε分別為湍動(dòng)能和耗散率；μt為湍流黏性系數(shù)；Gk和Gb分別為平均速度梯度和浮力引起的湍動(dòng)能產(chǎn)生項(xiàng)；σk和σε分別為k和ε的普朗特?cái)?shù)。常數(shù)σk=1.0、σε=1.3、Cε1=1.44、Cε2=1.92。

1.2 多孔介質(zhì)模型

將養(yǎng)殖網(wǎng)衣視為多孔透水板，采用多孔介質(zhì)模型進(jìn)行數(shù)值模擬。在數(shù)值模擬中，以動(dòng)量方程的右邊增加源項(xiàng)來模擬多孔介質(zhì)，動(dòng)量方程[13]變?yōu)椋?/p>

(5)

式(5)中：Si為源項(xiàng)，其表達(dá)式為：

Si=-(Pνu+Piuu)

(6)

式(6)中：Pν和Pi分別為多孔介質(zhì)黏性力系數(shù)矩陣和慣性力系數(shù)矩陣。

水流作用在多孔介質(zhì)區(qū)域的合力(F)[13]為：

F=SitA

(7)

式(7)中：t為多孔介質(zhì)的厚度；A為多孔介質(zhì)的面積。

可將F分解為與水流方向相平行的阻力(D)和與水流方向相垂直的升力(L)，當(dāng)水流方向與網(wǎng)衣平面垂直時(shí)(攻角α=90°，如圖1所示)，D和L的表達(dá)式[13]分別為：

(8)

(9)

網(wǎng)衣密實(shí)度(S)的表達(dá)式[13]為：

(10)

式(10)中：d為網(wǎng)線直徑；λ為網(wǎng)目長度。

網(wǎng)衣密實(shí)度與多孔介質(zhì)孔隙率的關(guān)系為：

S+χ=1

(11)

式(11)中：χ為多孔介質(zhì)的孔隙率。

2 數(shù)值計(jì)算方法

2.1 計(jì)算模型

本文研究對(duì)象為Patursson單片網(wǎng)衣試驗(yàn)過程中的平面網(wǎng)衣[6]，網(wǎng)衣尺寸的具體參數(shù)見表1。

表1 網(wǎng)衣的尺寸

網(wǎng)衣試驗(yàn)的水槽長37.00 m，寬3.66 m。網(wǎng)衣距離水槽兩側(cè)1.33 m，位于水槽液面0.73 m，以減少興波的影響。網(wǎng)衣通過鋼管固定在水槽中心，鋼管直徑為0.033 m，鋼管可繞中心軸轉(zhuǎn)動(dòng)，從而改變網(wǎng)衣的攻角[6]。網(wǎng)衣布置圖如圖2所示。

2.2 網(wǎng)格生成

數(shù)值計(jì)算的計(jì)算域和邊界條件如圖3所示。計(jì)算域的長度為5 m，寬度和高度與試驗(yàn)水槽相同，多孔介質(zhì)模型距離進(jìn)流面入口1.5 m，多孔區(qū)域厚度為0.05 m。計(jì)算域的邊界條件設(shè)置如下：進(jìn)流面設(shè)置為速度入口條件；出流面為壓力出口條件；水槽四周為無滑移壁面條件；水槽與多孔區(qū)域的交界面為內(nèi)部界面條件。

數(shù)值計(jì)算采用的網(wǎng)格類型為非結(jié)構(gòu)切割體網(wǎng)格。在網(wǎng)衣周圍區(qū)域和網(wǎng)衣后側(cè)區(qū)域進(jìn)行了網(wǎng)格細(xì)化，用于捕捉網(wǎng)衣后側(cè)的流場(chǎng)。計(jì)算域的總網(wǎng)格數(shù)為31×104，其中加密區(qū)域1和加密區(qū)域2的網(wǎng)格為18×104和8×104。整體計(jì)算域網(wǎng)格和中間水平截面網(wǎng)格如圖4所示。

2.3 計(jì)算方法

在數(shù)值計(jì)算中，求解器采用基于分離流的黏性求解器，壓力速度耦合算法為SIMPLE法，對(duì)流項(xiàng)采用二階離散格式。湍流模型采用標(biāo)準(zhǔn)k-ε湍流模型，壁面函數(shù)采用高y+壁面處理。采用三維定常計(jì)算，最大迭代步數(shù)為1 000步。

對(duì)網(wǎng)衣試驗(yàn)結(jié)果進(jìn)行擬合，多孔介質(zhì)的法向慣性系數(shù)和切向慣性系數(shù)分別為2 492.5、830.0kg/m4，法向黏性系數(shù)和切向黏性系數(shù)分別為75.08、38.31kg/(m3·s)[9]，多孔介質(zhì)的孔隙率為0.802。

3 數(shù)值計(jì)算結(jié)果

3.1 不確定度分析

網(wǎng)衣阻力系數(shù)和升力系數(shù)的定義[13]為：

(12)

(13)

式(12)～(13)中：Cd為網(wǎng)衣阻力系數(shù)；D為網(wǎng)衣所受的阻力；Cl為網(wǎng)衣升力系數(shù)；L為網(wǎng)衣所受的升力。

表2 網(wǎng)格的計(jì)算結(jié)果

可以看出，3套網(wǎng)格在數(shù)值迭代計(jì)算75步后，網(wǎng)衣阻力都已經(jīng)完全收斂，將數(shù)值模擬的總迭代步數(shù)設(shè)置為75～100，可節(jié)省計(jì)算工作量。1號(hào)網(wǎng)格、2號(hào)網(wǎng)格和3號(hào)網(wǎng)格的網(wǎng)衣阻力系數(shù)分別為0.263、0.262和0.260。采用國際拖曳水池(ITTC)的不確定計(jì)算方法對(duì)網(wǎng)衣阻力進(jìn)行不確定度分析，不確定度分析分為驗(yàn)證與確認(rèn)兩個(gè)過程[14]。

1)驗(yàn)證

數(shù)值不確定度由迭代不確定度(UI)、網(wǎng)格不確定度(UG)以及其他因素的不確定度(UP)組成。本文數(shù)值模擬采用定常計(jì)算，網(wǎng)衣阻力曲線幾乎無波動(dòng)，迭代誤差可以忽略不計(jì)，驗(yàn)證過程主要為計(jì)算數(shù)值模擬中的網(wǎng)格不確定度UG。

網(wǎng)格的收斂半徑RG為：

(14)

式(14)中：ε21為細(xì)網(wǎng)格與中網(wǎng)格的差值；ε32為中網(wǎng)格與粗網(wǎng)格的差值。

網(wǎng)格收斂半徑為0.281，0

(15)

(16)

(17)

(18)

2)確認(rèn)

確認(rèn)是利用試驗(yàn)數(shù)據(jù)來評(píng)估數(shù)值模型不確定度USN的過程，將數(shù)值模擬值(S)與試驗(yàn)值(D)進(jìn)行比較。

(19)

(20)

式(19)～(20)中：E為比較誤差；UV為數(shù)值結(jié)果確認(rèn)過程的不確定度；UD為試驗(yàn)結(jié)果的不確定度。

當(dāng)E>UV，說明數(shù)值結(jié)果與試驗(yàn)值偏離較遠(yuǎn)，需要改進(jìn)數(shù)值計(jì)算方法。

3.2 不同流速下的計(jì)算結(jié)果

在來流為0.125、0.250、0.500和0.750 m/s等4個(gè)速度下，對(duì)攻角為90°的網(wǎng)衣流場(chǎng)進(jìn)行數(shù)值模擬，并與Patursson O的網(wǎng)衣阻力試驗(yàn)結(jié)果[6]進(jìn)行比較，以驗(yàn)證數(shù)值模擬方法的準(zhǔn)確性及有效性。

網(wǎng)衣阻力系數(shù)計(jì)算結(jié)果和試驗(yàn)結(jié)果對(duì)比如圖6所示。由圖可見，當(dāng)來流為0.125 m/s時(shí)，網(wǎng)衣阻力的數(shù)值模擬結(jié)果偏大，誤差為13.4%，主要原因?yàn)榈退贂r(shí)網(wǎng)衣阻力的數(shù)值較小，相對(duì)誤差較大；當(dāng)來流分別為0.250、0.500和0.750 m/s時(shí)，網(wǎng)衣阻力的誤差分別為-1.85%、1.75%和-1.11%。整體上，網(wǎng)衣阻力的數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果吻合良好，說明本文數(shù)值模擬方法有效，可用于養(yǎng)殖網(wǎng)箱流場(chǎng)的數(shù)值預(yù)報(bào)。

在4個(gè)速度下，網(wǎng)衣中心截面的速度分布如圖7所示。由圖7可見，網(wǎng)衣邊緣端部的流體速度大約會(huì)增加2%；由于網(wǎng)衣的阻擋作用，網(wǎng)衣前后的流體速度都會(huì)降低，后方的流體速度大約會(huì)降低10%；隨著來流速度的增加，網(wǎng)衣兩端流體速度增加的區(qū)域逐漸變小，網(wǎng)衣后方對(duì)來流速度的阻擋作用逐漸減弱。

在網(wǎng)衣中心法方向設(shè)置50個(gè)監(jiān)測(cè)點(diǎn)以監(jiān)測(cè)網(wǎng)衣前后流場(chǎng)的速度變化，在4個(gè)速度下，網(wǎng)衣中心前后的速度分布如圖8所示。由圖8可見，在網(wǎng)衣前后(網(wǎng)前0.5倍網(wǎng)衣長度到網(wǎng)后1.0倍網(wǎng)衣長度的范圍內(nèi))，流體速度下降較快，網(wǎng)衣對(duì)流體有明顯的阻擋作用；超過1.0倍網(wǎng)衣長度，流體速度基本不再變化。

3.3 不同攻角下的計(jì)算結(jié)果

在來流為0.5 m/s的速度下，對(duì)攻角分別為90°、60°、45°、30°和0°的網(wǎng)衣流場(chǎng)進(jìn)行數(shù)值模擬，在不同攻角情況下，網(wǎng)衣阻力系數(shù)和升力系數(shù)的計(jì)算結(jié)果和試驗(yàn)結(jié)果對(duì)比分別如圖9和圖10所示。由圖9可見，不同攻角下網(wǎng)衣阻力的數(shù)值模擬結(jié)果偏大；隨著攻角角度的減小，網(wǎng)衣阻力系數(shù)不斷減小，而誤差逐漸增大。由圖10可見，不同攻角下網(wǎng)衣升力的數(shù)值模擬與試驗(yàn)結(jié)果吻合較好；隨攻角角度的增加，網(wǎng)衣升力系數(shù)先增大后減小；當(dāng)攻角為45°時(shí)，網(wǎng)衣受到的升力最大。

在不同攻角下，網(wǎng)衣中心截面的速度分布如圖11所示。由圖11可見，當(dāng)攻角變化后，網(wǎng)衣后面減流區(qū)域的寬度也發(fā)生變化，寬度基本上為網(wǎng)衣的投影寬度；網(wǎng)衣上端部的流體速度大于下端部的流體速度，網(wǎng)衣下部對(duì)流體的阻擋作用大于上部。

在不同攻角下，網(wǎng)衣中心前后的速度分布如圖12所示。由圖12可見，當(dāng)攻角在30°至90°變化時(shí)，網(wǎng)衣中心前后速度變化的規(guī)律基本相同。當(dāng)攻角為0°(來流與網(wǎng)衣平行)時(shí)，受網(wǎng)衣的阻擋，在網(wǎng)后1.0倍網(wǎng)衣長度范圍內(nèi)，流體速度小于其他攻角的速度；在網(wǎng)后大于1.0倍網(wǎng)衣長度時(shí)，發(fā)生繞射現(xiàn)象，最終流體速度反而大于其他攻角的速度。

4 結(jié)論

本文基于CFD理論，采用多孔介質(zhì)模型對(duì)單片網(wǎng)衣的流場(chǎng)進(jìn)行了數(shù)值模擬，首先進(jìn)行了不確定度分析，然后對(duì)網(wǎng)衣在不同流速和不同攻角下的流場(chǎng)進(jìn)行了數(shù)值計(jì)算，再將數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行了對(duì)比，并對(duì)網(wǎng)衣流場(chǎng)特性進(jìn)行了分析。結(jié)果如下：

1)網(wǎng)衣阻力和升力的數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果吻合良好，說明本文數(shù)值模擬方法有效，可用于養(yǎng)殖網(wǎng)箱流場(chǎng)的數(shù)值預(yù)報(bào)。

2)在網(wǎng)前0.5倍網(wǎng)衣長度到網(wǎng)后1.0倍網(wǎng)衣長度的范圍內(nèi)，流體速度下降較快，網(wǎng)衣對(duì)流體有明顯的阻擋作用；當(dāng)超過1.0倍網(wǎng)衣長度時(shí)，流體速度基本不再變化。

3)隨著來流速度的增加，網(wǎng)衣后方對(duì)來流速度的阻擋作用逐漸減弱。

采用本文數(shù)值模擬方法可以準(zhǔn)確獲取網(wǎng)衣周圍的流場(chǎng)數(shù)據(jù)，在后續(xù)工作中，將進(jìn)行本文數(shù)值模擬方法與網(wǎng)衣變形程序之間相互耦合的研究，從而可對(duì)具有百萬級(jí)網(wǎng)目的柔性網(wǎng)衣進(jìn)行求解。