基于CFD的剛性養(yǎng)殖網(wǎng)衣流場數(shù)值模擬及不確定度分析

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

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

漁業(yè)過度捕撈使?jié)O業(yè)資源量下降，單純依靠海洋捕撈早已無法滿足人類對水產(chǎn)品日益增長的需求，現(xiàn)在迫切需要海水養(yǎng)殖來彌補(bǔ)水產(chǎn)品消費市場的不足。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)衣對海水有阻流作用，這不僅會影響?zhàn)B殖對象生長，還會阻礙水體交換，影響網(wǎng)箱內(nèi)部的養(yǎng)殖環(huán)境，因此網(wǎng)衣內(nèi)部的流場特性對網(wǎng)箱養(yǎng)殖非常重要。

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

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

1 理論模型

1.1 控制方程

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

(1)

(2)

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

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

(3)

(4)

式(3)～(4)中：k和ε分別為湍動能和耗散率；μt為湍流黏性系數(shù)；Gk和Gb分別為平均速度梯度和浮力引起的湍動能產(chǎn)生項；σk和σε分別為k和ε的普朗特數(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ù)值模擬中，以動量方程的右邊增加源項來模擬多孔介質(zhì)，動量方程[13]變?yōu)椋?/p>

(5)

式(5)中：Si為源項，其表達(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)衣平面垂直時(攻角α=90°，如圖1所示)，D和L的表達(dá)式[13]分別為：

(8)

(9)

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

(10)

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

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

S+χ=1

(11)

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

2 數(shù)值計算方法

2.1 計算模型

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

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

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

2.2 網(wǎng)格生成

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

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

2.3 計算方法

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

對網(wǎng)衣試驗結(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ù)值計算結(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)格的計算結(jié)果

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

1)驗證

數(shù)值不確定度由迭代不確定度(UI)、網(wǎng)格不確定度(UG)以及其他因素的不確定度(UP)組成。本文數(shù)值模擬采用定常計算，網(wǎng)衣阻力曲線幾乎無波動，迭代誤差可以忽略不計，驗證過程主要為計算數(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)是利用試驗數(shù)據(jù)來評估數(shù)值模型不確定度USN的過程，將數(shù)值模擬值(S)與試驗值(D)進(jìn)行比較。

(19)

(20)

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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