基于CFD的剛性養殖網衣流場數值模擬及不確定度分析

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

(1. 集美大學輪機工程學院，福建 廈門 361021；2. 福建省水產研究所，福建省海洋生物增養殖與高值化利用重點實驗室，福建 廈門361013)

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

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

為研究及分析養殖網衣周圍的流場特性，本文基于CFD理論，采用多孔介質模型對單片網衣的流場進行數值模擬，首先進行不確定度分析，再對網衣在不同流速下的流場進行數值計算，并與試驗結果進行了比較，最后對網衣在不同攻角下的流場進行數值計算，分析了網衣的流場特性，給出網衣對流場內速度影響的量化區域及范圍，旨在能夠較好地為實際養殖生產活動提供參考依據。

1 理論模型

1.1 控制方程

流體力學問題要遵守質量守恒、動量守恒和能量守恒定律，對于不可壓縮流體，其連續性方程和動量守恒方程分別如下[12]：

(1)

(2)

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

控制方程不封閉，因而需引入湍流模型，使方程能夠求解。本文計算采用標準k-ε湍流模型[12]，該模型具有較髙的穩定性、經濟性和計算精度，應用廣泛，其表達形式如下：

(3)

(4)

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

1.2 多孔介質模型

將養殖網衣視為多孔透水板，采用多孔介質模型進行數值模擬。在數值模擬中，以動量方程的右邊增加源項來模擬多孔介質，動量方程[13]變為：

(5)

式(5)中：Si為源項，其表達式為：

Si=-(Pνu+Piuu)

(6)

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

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

F=SitA

(7)

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

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

(8)

(9)

網衣密實度(S)的表達式[13]為：

(10)

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

網衣密實度與多孔介質孔隙率的關系為：

S+χ=1

(11)

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

2 數值計算方法

2.1 計算模型

本文研究對象為Patursson單片網衣試驗過程中的平面網衣[6]，網衣尺寸的具體參數見表1。

表1 網衣的尺寸

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

2.2 網格生成

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

數值計算采用的網格類型為非結構切割體網格。在網衣周圍區域和網衣后側區域進行了網格細化，用于捕捉網衣后側的流場。計算域的總網格數為31×104，其中加密區域1和加密區域2的網格為18×104和8×104。整體計算域網格和中間水平截面網格如圖4所示。

2.3 計算方法

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

對網衣試驗結果進行擬合，多孔介質的法向慣性系數和切向慣性系數分別為2 492.5、830.0kg/m4，法向黏性系數和切向黏性系數分別為75.08、38.31kg/(m3·s)[9]，多孔介質的孔隙率為0.802。

3 數值計算結果

3.1 不確定度分析

網衣阻力系數和升力系數的定義[13]為：

(12)

(13)

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

表2 網格的計算結果

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

1)驗證

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

網格的收斂半徑RG為：

(14)

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

網格收斂半徑為0.281，0

(15)

(16)

(17)

(18)

2)確認

確認是利用試驗數據來評估數值模型不確定度USN的過程，將數值模擬值(S)與試驗值(D)進行比較。

(19)

(20)

式(19)～(20)中：E為比較誤差；UV為數值結果確認過程的不確定度；UD為試驗結果的不確定度。

當E>UV，說明數值結果與試驗值偏離較遠，需要改進數值計算方法。

3.2 不同流速下的計算結果

在來流為0.125、0.250、0.500和0.750 m/s等4個速度下，對攻角為90°的網衣流場進行數值模擬，并與Patursson O的網衣阻力試驗結果[6]進行比較，以驗證數值模擬方法的準確性及有效性。

網衣阻力系數計算結果和試驗結果對比如圖6所示。由圖可見，當來流為0.125 m/s時，網衣阻力的數值模擬結果偏大，誤差為13.4%，主要原因為低速時網衣阻力的數值較小，相對誤差較大；當來流分別為0.250、0.500和0.750 m/s時，網衣阻力的誤差分別為-1.85%、1.75%和-1.11%。整體上，網衣阻力的數值模擬結果與試驗結果吻合良好，說明本文數值模擬方法有效，可用于養殖網箱流場的數值預報。

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

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

3.3 不同攻角下的計算結果

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

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

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

4 結論

本文基于CFD理論，采用多孔介質模型對單片網衣的流場進行了數值模擬，首先進行了不確定度分析，然后對網衣在不同流速和不同攻角下的流場進行了數值計算，再將數值計算結果與試驗結果進行了對比，并對網衣流場特性進行了分析。結果如下：

1)網衣阻力和升力的數值模擬結果與試驗結果吻合良好，說明本文數值模擬方法有效，可用于養殖網箱流場的數值預報。

2)在網前0.5倍網衣長度到網后1.0倍網衣長度的范圍內，流體速度下降較快，網衣對流體有明顯的阻擋作用；當超過1.0倍網衣長度時，流體速度基本不再變化。

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

采用本文數值模擬方法可以準確獲取網衣周圍的流場數據，在后續工作中，將進行本文數值模擬方法與網衣變形程序之間相互耦合的研究，從而可對具有百萬級網目的柔性網衣進行求解。