小尺度平面網(wǎng)衣流場特性數(shù)值研究

李鵬，韓鑫悅，秦洪德，鄧忠超

1.哈爾濱工程大學(xué) 船舶工程學(xué)院, 黑龍江 哈爾濱 150001

2.哈爾濱工程大學(xué) 煙臺研究 (生)院, 山東 煙臺 264006

水產(chǎn)養(yǎng)殖為人類的日常生活提供食物和生物能源，據(jù)2018年聯(lián)合國糧食及農(nóng)業(yè)組織報告，2030年，海產(chǎn)品需求可能達(dá)到約1.8億噸[1]。網(wǎng)箱養(yǎng)殖具有清污能力強(qiáng)、養(yǎng)殖產(chǎn)品質(zhì)量高、對環(huán)境污染小等特點(diǎn)，網(wǎng)衣能限制海洋生物活動范圍，保證網(wǎng)箱內(nèi)外的水質(zhì)交換，抗風(fēng)浪破壞。而金屬網(wǎng)衣具有耐腐蝕、抗風(fēng)浪能力強(qiáng)、不易附著物等優(yōu)點(diǎn)，應(yīng)用逐漸廣泛。研究金屬網(wǎng)衣的阻力特性及其周圍的流場特性關(guān)系著魚群養(yǎng)殖的質(zhì)量和網(wǎng)箱安全性。

在過去的幾十年里，關(guān)于網(wǎng)衣的阻力特性及其周圍的流場特性一直是熱門課題，不同的學(xué)者進(jìn)行了大量的研究。L?land[2]采用線性自由尾跡方程結(jié)合漩渦粘公式模擬了平面網(wǎng)周圍的流場，研究了近場和遠(yuǎn)場尾流的運(yùn)動特性，并與模型試驗結(jié)果進(jìn)行了比較。Tsukrov等[3]利用經(jīng)驗公式計算求得的網(wǎng)衣所受力的大小，利用等效網(wǎng)單元模型對網(wǎng)衣在波浪水流聯(lián)合作用下的運(yùn)動相應(yīng)進(jìn)行模擬。詹杰民等[4]采用物理實(shí)驗和理論分析的方法，研究了網(wǎng)衣的水動力性能與雷諾數(shù)、攻角及網(wǎng)衣特征參數(shù)之間的關(guān)系。桂福坤[5]基于有效繞流面積系數(shù)假設(shè)理論，得出了平面網(wǎng)衣系統(tǒng)流速衰減的計算方法。Patursson[6]利用多孔介質(zhì)模型來模擬平面網(wǎng)衣，研究分析平面網(wǎng)衣后的水流流速衰減情況。趙云鵬等[7?8]建立三維數(shù)值模型來模擬水流中平面網(wǎng)衣周圍的流場，分析了不同攻角、不同高度和不同距離的兩片網(wǎng)衣以及不同數(shù)量的平面網(wǎng)衣周圍的流場特性。畢春偉等[9]運(yùn)用集中質(zhì)量法和多孔介質(zhì)法數(shù)值模擬了流場和網(wǎng)衣之間的相互作用，然后計算分析了單片柔性網(wǎng)衣與不同間距的雙片柔性網(wǎng)衣周圍的流場特性。Liu等[10]基于有限元分析方法建立了動力有限元數(shù)值模型，對有周期性外力作用在網(wǎng)衣上的動力響應(yīng)問題分析研究，研究發(fā)現(xiàn)雷諾數(shù)的大小和振動頻率的高低都會對網(wǎng)衣的水動力系數(shù)產(chǎn)生影響。Hao等[11?13]提出了一種基于Morison型荷載模型轉(zhuǎn)換的多孔阻力系數(shù)計算方法。

本文采用數(shù)值模擬方法，基于計算流體力學(xué)方法，利用STAR-CCM+軟件，對不同雷諾數(shù)下、不同攻角、不同密實(shí)度以及不同網(wǎng)眼形狀的平面網(wǎng)衣的阻力性能和網(wǎng)衣周圍流場特性進(jìn)行了研究分析。本文對明確平面網(wǎng)衣的流場特性具有一定的指導(dǎo)價值。

1 數(shù)值模型

1.1 控制方程和湍流方程

本文使用CFD方法進(jìn)行數(shù)值模擬時，基于不可壓縮粘性流體的運(yùn)動，通過有限體積法求解流場的控制方程，不可壓縮流的運(yùn)動可以由雷諾平均Navier-Stokes方程系統(tǒng)描述，連續(xù)方程和動量方程可以表示為

式中：u、v、w為速度U在3個方向上的分量，t為時間，p為流體受到的壓力，μ為流體的動力粘性系數(shù)，ui為x、y、z方向的速度分量（i=1,2,3），u'為脈動速度。

在湍流運(yùn)動包含著脈動應(yīng)力項，該項使方程不封閉，因而必須對該項進(jìn)行假設(shè)，即湍流模型。因為網(wǎng)衣對水的流動有著一定的阻攔作用，故網(wǎng)衣前后會存在不可忽略的壓力梯度，Realizablek-ε模型能夠?qū)ι鲜霈F(xiàn)象進(jìn)行很好地模擬分析[14]，所以本文選擇 Realizablek-ε模型進(jìn)行對湍流流動的模擬分析。

k方程：

ε方程：

式中：

k方程、ε方程中：Gk為由平均速度梯度產(chǎn)生的湍流動能，Gb為由浮力產(chǎn)生的湍流動能，YM為可壓縮湍流中脈動擴(kuò)散分布的整體耗散率，C2、C1ε和C3ε為常量，S為應(yīng)變模量，σk和σε為k方程和ε方程的湍流Prandtl數(shù)，Sk和Sε為用戶自定義的源項。Gb=0，YM=0，Sk=0，Sε=0。對于 Realizablek-ε模型，C1ε=1.44，C2=1.9，σk=1.0，σε=1.2。

1.2 相關(guān)參數(shù)

與本文中流體運(yùn)動相關(guān)的無量綱系數(shù)為

雷諾數(shù)：

阻力系數(shù)：

式中：u為流體的速度，d為網(wǎng)線的直徑，S為投影面積，v為流體運(yùn)動粘性系數(shù)，F(xiàn)D為物體阻力。

攻角 (angle of attack, AOA)α定義為流體正方向與網(wǎng)衣方向的夾角，如圖1所示。

圖1 攻角的定義

網(wǎng)眼密實(shí)度Sn定義為網(wǎng)線所覆蓋的面積與平面網(wǎng)的總面積之比（L?land[2]）。

對于方形網(wǎng)衣，參考文獻(xiàn)[15]文章，其網(wǎng)眼密實(shí)度為

對于菱形網(wǎng)衣，參考文獻(xiàn)[3]。如圖2所示，在a×b矩形網(wǎng)衣有X×Y個菱形網(wǎng)格的情況下，對菱形網(wǎng)衣密實(shí)度定義為

圖2 菱形網(wǎng)眼網(wǎng)衣

式中：a為菱形網(wǎng)衣長度，b為菱形網(wǎng)衣寬度，X為長度方向為菱形網(wǎng)格個數(shù)，Y為寬度方向上菱形網(wǎng)格個數(shù)，d為網(wǎng)衣網(wǎng)線直徑。

目前沒有對六邊形網(wǎng)衣的密實(shí)度有確定的計算公式，故六邊形網(wǎng)衣密實(shí)度Sn參考文獻(xiàn)[2]，為網(wǎng)線所覆蓋的面積與平面網(wǎng)的總面積之比。

2 網(wǎng)格無關(guān)性驗證

數(shù)值模擬使用的網(wǎng)衣為無結(jié)金屬網(wǎng)，網(wǎng)線可簡化為剛體，Tsukrov等[3]運(yùn)用物理試驗法研究分析純流作用下不同材料網(wǎng)衣的阻力系數(shù)，本文選取其中的金屬網(wǎng)進(jìn)行數(shù)值模擬，考慮到后續(xù)計算的簡便性，本文選取3×3小尺度網(wǎng)衣進(jìn)行研究分析，直徑dw為 0.001 35 m，網(wǎng)眼長度lw為 0.018 m，如圖3所示。

圖3 無結(jié)金屬網(wǎng)衣

定義寬度方向為y軸方向，長度方向為x軸方向，為水流前進(jìn)方向，垂直方向為z軸方向。設(shè)3×3網(wǎng)衣的總長度為l，流體進(jìn)口距網(wǎng)衣中心4l，流體出口距網(wǎng)衣中心8l，左右上下各距網(wǎng)衣中心2l。流體域左側(cè)為速度進(jìn)口，右側(cè)為壓力出口，其余邊界為對稱平面，用以模擬無限流體域。網(wǎng)格采用切割體網(wǎng)格，對網(wǎng)衣和其壁面區(qū)域、網(wǎng)衣后尾流區(qū)域進(jìn)行加密。圖4為建立好的計算域和網(wǎng)格模型。采用4種不同尺寸的網(wǎng)格進(jìn)行網(wǎng)格無關(guān)性分析，得到了不同數(shù)量網(wǎng)格網(wǎng)衣的阻力系數(shù)Cd。表1為網(wǎng)格無關(guān)性驗證結(jié)果。

圖4 計算域及網(wǎng)格

表1 網(wǎng)格無關(guān)性

當(dāng)網(wǎng)格數(shù)量為127萬和172萬時，兩者阻力系數(shù)相對誤差為0.086%，說明已收斂，驗證了網(wǎng)格無關(guān)性?？紤]數(shù)值的準(zhǔn)確性，網(wǎng)格數(shù)量采用127萬。

3.1 雷諾數(shù)的影響

圖5為α=90°時平面網(wǎng)衣在不同流速下的阻力系數(shù)，可以看出，隨著雷諾數(shù)的增大，阻力系數(shù)逐漸減小。當(dāng)速度從 0.1 m/s增至 0.5 m/s時，阻力系數(shù)平均降低了9.6%；當(dāng)速度從1 m/s增至3 m/s時，阻力系數(shù)平均降低了3.7%，隨著雷諾數(shù)的增加，平面網(wǎng)衣的阻力系數(shù)趨于平緩。

圖5 不同流速 (雷諾數(shù))下阻力系數(shù)

圖6為α=90°時平面網(wǎng)衣在不同流速下的速度場分布?？梢钥闯?，流體流經(jīng)網(wǎng)孔時流速增大，流過網(wǎng)衣后流速慢慢降低，但其尾流區(qū)域持續(xù)較長的時間。然而，對比不同雷諾數(shù)下速度分布的變化，可以發(fā)現(xiàn)雷諾數(shù)的變化并不影響尾流的寬度和區(qū)域。

圖6 不同雷諾數(shù)下網(wǎng)衣速度分布

為了更準(zhǔn)確地描述速度的分布，繪制沿著網(wǎng)中心線的x方向流速衰減情況圖，如圖7所示。

圖7 不同雷諾數(shù)下速度衰減

其流速衰減幅度采用式(1)來計算：

式中：u0為速度入口處速度，u為網(wǎng)衣中心線上網(wǎng)孔處速度。

從圖7可看出，不同雷諾數(shù)下速度衰減范圍為5%~6.5%，可以近似認(rèn)為雷諾數(shù)并不影響速度衰減。

3.2 AOA 的影響

圖8為平面網(wǎng)衣在不同攻角下的阻力系數(shù)。可以看出，隨攻角增大，阻力系數(shù)逐漸增大，當(dāng)攻角在45°附近時，阻力系數(shù)增大相對明顯。

圖8 不同攻角下阻力系數(shù)

圖9為平面網(wǎng)衣在不同攻角下的速度分布?？梢钥闯?，隨著攻角的增大，近場區(qū)域出現(xiàn)流速衰減范圍越大，遠(yuǎn)場區(qū)域流速增大的范圍也越大；網(wǎng)衣后尾流的寬度也逐漸增大，且流場越來越紊亂。

圖9 不同攻角下網(wǎng)衣速度分布

為了更準(zhǔn)確地描述速度的分布，繪制了沿著網(wǎng)中心線的x方向速度分布情況。圖10為不同攻角下網(wǎng)衣中心軸線速度。

圖10 不同攻角下網(wǎng)衣中心軸線速度

通過圖10可以看出，當(dāng)α≥30°時，水流流經(jīng)網(wǎng)衣的過程中，速度先減小后增大，當(dāng)攻角越大時，速度變化的幅值越小；水流流過網(wǎng)衣后，水流速度會恢復(fù)至0.5 m/s，攻角越小越快出現(xiàn)恢復(fù)的趨勢。當(dāng)α=15°時，水流流過網(wǎng)衣時，速度出現(xiàn)波動明顯，呈現(xiàn)出先大幅減小，后小幅增大，而后大幅減小的趨勢，最后趨于平緩。

3.3 Sn 的影響

圖11為網(wǎng)衣在不同密實(shí)度下的阻力系數(shù)，隨密實(shí)度的增大，阻力系數(shù)逐漸減??；當(dāng)密實(shí)度由0.15增至0.2，阻力系數(shù)平均降低了10.1%，密實(shí)度由0.25增至0.3，阻力系數(shù)平均降低了4.2%，隨著密實(shí)度的增加，阻力系數(shù)衰減的速率減緩。

圖11 不同密實(shí)度下阻力系數(shù)

圖12為網(wǎng)衣在不同密實(shí)度下的速度分布，可以看出，隨著密實(shí)度的增大，近場和遠(yuǎn)場區(qū)域出現(xiàn)流速衰減的范圍都越大，網(wǎng)衣后尾流的寬度逐漸增大，并且尾流越來越紊亂。流體流經(jīng)網(wǎng)孔時，流速明顯增大，而隨著密實(shí)度的增大，其流速增大區(qū)域也越大。為了更準(zhǔn)確地描述隨密實(shí)度變化速度的分布，繪制沿著網(wǎng)中心線的x方向速度的分布情況。圖13為不同攻角下網(wǎng)衣中心軸線速度。

圖12 不同密實(shí)度下網(wǎng)衣速度分布

圖13 不同攻角下網(wǎng)衣中心軸線速度

通過圖13可以看出，在水流流經(jīng)網(wǎng)衣的過程中，速度呈現(xiàn)先減小后增大的趨勢。隨著密實(shí)度的增大，速度變化的趨勢也越大。即網(wǎng)眼密實(shí)度越大，網(wǎng)衣網(wǎng)線對軸線上的速度影響越明顯。

為了更準(zhǔn)確地描述速度的分布，繪制沿著網(wǎng)中心線的x方向流速衰減情況圖，速度衰減幅度定義采用式(1)。

圖14為不同密實(shí)度下速度衰減變化，可以看出，當(dāng)密實(shí)度Sn＜0.2時，速度衰減Δu小于6%，當(dāng)密實(shí)度Sn=0.3時，其速度衰減Δu超過了10%，即密實(shí)度對速度衰減的影響較大。隨著密實(shí)度的增大，速度衰減逐漸增大，速度衰減增大的速率基本不變。

3.4 網(wǎng)眼形狀的影響

圖15為平面網(wǎng)衣在不同密實(shí)度下的阻力變化，可以看出，當(dāng)水流速度相同時，菱形網(wǎng)眼網(wǎng)衣所受的阻力最大，方形網(wǎng)網(wǎng)眼衣所受阻力最小。當(dāng)水流速度較小時，網(wǎng)眼形狀對網(wǎng)衣所受阻力的大小影響較小；隨著水流速度的增大，不同網(wǎng)眼形狀網(wǎng)衣所受的阻力相差越大。

圖16為平面網(wǎng)衣在不同密實(shí)度下的速度分布情況，可以看出，水流流經(jīng)網(wǎng)衣的過程中，近場區(qū)域流速都發(fā)生了衰減，網(wǎng)眼類型近場區(qū)域流速衰減程度影響不大。網(wǎng)眼類型對遠(yuǎn)場區(qū)域流速衰減現(xiàn)象影響較大，六邊形網(wǎng)眼網(wǎng)衣后方流速衰減的幅值比較小，且流速越快恢復(fù)至0.5 m/s；菱形網(wǎng)眼網(wǎng)衣后方流速衰減幅值較大，且出現(xiàn)流速衰減的區(qū)域較大。

圖16 不同網(wǎng)眼形狀網(wǎng)衣速度分布

為了更準(zhǔn)確地描述流場的狀態(tài)，繪制了沿著網(wǎng)中心線的x方向速度分布情況，如圖17所示。通過圖17可以看出，在水流流經(jīng)網(wǎng)衣的過程中，速度都呈現(xiàn)先減小后增大的趨勢，且速度最后都會恢復(fù)至0.5 m/s。在水流流經(jīng)網(wǎng)衣的最初階段，水流速度會出現(xiàn)衰減，其中六邊形網(wǎng)眼網(wǎng)衣較慢出現(xiàn)減小的趨勢，且減小幅度為0.83%，方形網(wǎng)眼網(wǎng)衣流速衰減幅度為5.2%，即網(wǎng)眼形狀對流速衰減影響較大；水流流經(jīng)網(wǎng)衣時，水流速度會出現(xiàn)驟增，其中方形網(wǎng)眼網(wǎng)衣流速增大幅度為9.7%，六邊形網(wǎng)眼網(wǎng)衣流速增幅為1.6%，網(wǎng)眼形狀對網(wǎng)孔后流速增加的影響不可忽略；水流流經(jīng)網(wǎng)衣的最后階段，六邊形網(wǎng)眼網(wǎng)衣的水流速度較快恢復(fù)到0.5 m/s。對水流流經(jīng)網(wǎng)衣中心線速度進(jìn)行對比分析可以發(fā)現(xiàn)，六邊形網(wǎng)眼形狀的網(wǎng)衣中心線速度變化最小。

圖17 不同網(wǎng)眼形狀網(wǎng)衣中心軸線速度

4 結(jié)論

本文基于CFD數(shù)值方法，研究均勻流作用下小尺度無結(jié)網(wǎng)衣在不同雷諾數(shù)、攻角、密實(shí)度和網(wǎng)眼形狀下的阻力性能和速度場分布。結(jié)論如下：

1）隨著雷諾數(shù)的增大，網(wǎng)衣阻力系數(shù)逐漸減小，且減小速率逐漸降低，雷諾數(shù)對速度分布中網(wǎng)衣后尾流區(qū)域和流經(jīng)網(wǎng)衣的速度衰減影響較小。

2）隨著AOA的增大，阻力系數(shù)逐漸增大，網(wǎng)衣后尾流的區(qū)域越大且越紊亂；攻角較小時，流經(jīng)網(wǎng)衣的水流速度的波動相對較明顯。

3）隨著密實(shí)度的增大，阻力系數(shù)逐漸減小阻力系數(shù)衰減的速率減緩，網(wǎng)衣后尾流的寬度越大，尾流越紊亂，網(wǎng)線對網(wǎng)衣軸線上的速度影響越明顯。

4）當(dāng)水流速度相同時，菱形網(wǎng)衣所受的阻力最大，方形網(wǎng)衣所受阻力最小；網(wǎng)眼類型對速度分布中網(wǎng)衣后尾流范圍影響不大，六邊形網(wǎng)眼形狀的網(wǎng)衣中心線速度變化最小。