于定勇， 楊遠(yuǎn)航， 李宇佳

(中國(guó)海洋大學(xué)工程學(xué)院，山東 青島 266100)

近年來隨著人類進(jìn)軍海洋步伐的加快，許多魚類棲息地遭到嚴(yán)重破壞，海洋漁業(yè)資源急劇衰退，因此許多國(guó)家提出了海洋牧場(chǎng)的建設(shè)構(gòu)想。根據(jù)《山東省人工魚礁建設(shè)規(guī)劃(2014—2020年)》，到2020年，山東省沿海區(qū)域?qū)⒔ㄔO(shè)9大人工魚礁帶，40個(gè)人工魚礁群。

人工魚礁體作為海洋牧場(chǎng)的重要組成部分，是設(shè)置于預(yù)定海域的人工構(gòu)筑物，可為魚類創(chuàng)建適宜的繁育、生長(zhǎng)、避敵場(chǎng)所，目前正受到越來越廣泛的關(guān)注。人工魚礁體投放入海后產(chǎn)生的上升流可引起海底營(yíng)養(yǎng)物質(zhì)卷起和擴(kuò)散，礁體附近形成較理想的索餌區(qū)，能改善魚群覓食條件。人工魚礁流場(chǎng)效應(yīng)是魚礁體影響漁業(yè)資源增殖的重要因素。

M. Falc?o等[1]通過對(duì)比人工魚礁區(qū)和非魚礁區(qū)營(yíng)養(yǎng)鹽通量和有機(jī)顆粒物的不同，研究了礁體影響營(yíng)養(yǎng)鹽通量和懸浮泥沙變化的水動(dòng)力機(jī)制。

Jinho Woo等[2]通過數(shù)值模擬給出了包括透空立方體在內(nèi)的24種不同形狀人工魚礁阻力系數(shù)值，結(jié)果表明阻力系數(shù)不隨初始流速的變化而變化，但其受礁體迎流角度的影響較為顯著。

Dongha Kim等[3]通過數(shù)值模擬研究了拱型、圓屋頂型等6類不同形狀礁體尾渦長(zhǎng)度與阻力系數(shù)的關(guān)系，表明礁體尾渦長(zhǎng)度與阻力系數(shù)沒有明顯的線性關(guān)系。

唐衍力[4]通過水槽實(shí)驗(yàn)測(cè)得邊長(zhǎng)為0.15 m、開口比為0.22的有蓋和無蓋兩種圓形開口魚礁模型在不同來流速度、不同迎流角度時(shí)所受阻力，得到礁體阻力系數(shù)。

劉彥[5]通過數(shù)值模擬研究了中空正方體魚礁單體及組合礁體在不同來流速度下的流場(chǎng)效應(yīng)，計(jì)算了開口比為0.36的正方體人工魚礁的抗滑移及抗傾覆安全系數(shù)。

趙云鵬等[6]建立三維數(shù)值波浪水槽研究波浪作用下三角型鏤空人工魚礁的受力情況，得到了三角型鏤空魚礁體阻力系數(shù)、慣性力系數(shù)隨Kc數(shù)、Re數(shù)的變化情況。

由上述研究工作可知目前國(guó)內(nèi)外學(xué)者對(duì)于人工魚礁體的研究主要集中于單一開口比情況水流作用下礁體周圍流場(chǎng)形態(tài)隨海流流速、礁體形狀、礁體迎流角度的變化情況及其生態(tài)效應(yīng)方面，但是對(duì)于礁體流場(chǎng)效應(yīng)、阻力系數(shù)、礁體受力及其穩(wěn)定性隨開口比的變化情況尚缺乏研究，已有研究工作[7]表明礁體開口比變化對(duì)上述流場(chǎng)特性有一定的影響。本文利用Fluent軟件模擬方型人工魚礁體水動(dòng)力性能，旨在解決礁體流場(chǎng)效應(yīng)、阻力系數(shù)、抗滑移及抗傾覆安全系數(shù)隨開口比的變化問題，為實(shí)際礁體結(jié)構(gòu)的設(shè)計(jì)提供參考。

1 模型

1.1 控制方程

假設(shè)礁體附近的流動(dòng)為粘性不可壓縮流體的湍流運(yùn)動(dòng)，溫度變化不大，能量方程可以忽略。

連續(xù)方程：

(1)

動(dòng)量方程：

(2)

式中：ui(i=1，2，3)分別為x、y、z方向的雷諾平均速度；ρ為流體密度；p為壓強(qiáng)；ν為運(yùn)動(dòng)粘性系數(shù)；fi為體積力。

1.2 湍流模型

本文在計(jì)算粘性流體運(yùn)動(dòng)時(shí)采用RNGκ-ε兩方程模型。此模型可以有效模擬分布較均勻、湍流結(jié)構(gòu)較小的湍流流動(dòng)[7]，適合人工魚礁體流場(chǎng)效應(yīng)的研究。

湍動(dòng)能κ方程：

(3)

湍流耗散率ε方程：

(4)

本文采用有限體積法離散控制方程，壓力-速度耦合采用SIMPLEC算法，壓力項(xiàng)處理采用標(biāo)準(zhǔn)差分格式，各方程空間離散均采用二階迎風(fēng)格式，計(jì)算殘差值取10-5，采用非結(jié)構(gòu)化網(wǎng)格，通過數(shù)值模擬方型礁體水動(dòng)力特性，分析其流場(chǎng)效應(yīng)、阻力系數(shù)隨開口比的變化情況并計(jì)算礁體抗滑移及抗傾覆安全系數(shù)。

1.3 魚礁結(jié)構(gòu)及模擬計(jì)算區(qū)域

4個(gè)側(cè)面具有相同開口形式的礁體，其迎流面在垂直于水流方向上開口的投影面積與迎流面投影面積之比稱為開口比(φ)[7]。圖1表示的是開口比為0.4的方型人工魚礁體結(jié)構(gòu)示意圖。

本文選取邊長(zhǎng)L為3 m、開口比φ分別為0、0.1、0.2、0.3、0.4、0.5、0.6的方形開口礁體，流速取為0.8 m/s，分析其周圍流場(chǎng)形態(tài)、阻力系數(shù)隨開口比的變化情況并計(jì)算其抗滑移、抗傾覆安全系數(shù)。仿真計(jì)算區(qū)域如圖2所示。

圖1 開口比為0.4的方型人工魚礁體結(jié)構(gòu)Fig.1 The sketch of cubic artificial reef with 0.4 opening ratio

圖2 計(jì)算區(qū)域Fig.2 Computational domain

文中所設(shè)定的初邊界條件如下：

(1)入口邊界設(shè)置為速度入口邊界條件(Velocity inlet)，來流速度為0.8 m/s，設(shè)置邊界上各方向的速度矢量分量，并給出邊界上湍動(dòng)能κ和湍動(dòng)耗散率ε。

(2)出口邊界設(shè)置為自由出流邊界條件(Outflow)。

(3)計(jì)算域的兩側(cè)面設(shè)置為對(duì)稱邊界(Symmetry)。

(4)計(jì)算域的頂面設(shè)置為具有與入口水流相同速度的可移動(dòng)壁面，剪切力為零，底面和礁體表面設(shè)置為無滑移壁面(Wall)。

2 模型可靠性驗(yàn)證

本文選擇的是渦粘模型中的RNGκ-ε兩方程模型，屬于湍流非直接數(shù)值模擬方法中的Reynolds平均法(RANS)，使用ANSYS Workbench 【Mesh】模塊對(duì)計(jì)算域進(jìn)行四面體單元非結(jié)構(gòu)化網(wǎng)格劃分。為了驗(yàn)證本文選擇的湍流模型、參數(shù)設(shè)置、網(wǎng)格劃分的準(zhǔn)確性，選擇唐衍力建立的邊長(zhǎng)為0.15 m，頂面有一個(gè)、每個(gè)側(cè)面有四個(gè)直徑均為0.04 m圓形開口，開口比為0.22的無底立方體人工魚礁[4]進(jìn)行數(shù)值模擬，如圖3所示。

數(shù)值模擬時(shí)最大網(wǎng)格尺寸設(shè)置為0.025 m，礁體表面第一層邊界層網(wǎng)格高度為0.002 m，增長(zhǎng)率為1.1，共設(shè)置10層。設(shè)定來流速度為0.5 m/s，將模擬結(jié)果與已有學(xué)者進(jìn)行的水槽實(shí)驗(yàn)和數(shù)值模擬結(jié)果進(jìn)行對(duì)比，測(cè)點(diǎn)相對(duì)位置分布如圖4所示[4]，結(jié)果如圖5及表1所示。

由對(duì)比結(jié)果可以看出，利用本文選擇的數(shù)值模型模擬得到的礁體前后測(cè)點(diǎn)的流速及礁體阻力系數(shù)與其他學(xué)者通過物理模型實(shí)驗(yàn)和數(shù)值模擬得到的結(jié)果相近。本文利用Fluent軟件模擬人工魚礁水動(dòng)力特性的研究方法是可行的。

圖3 圓形開口的人工魚礁體模型Fig.3 Artificial reef model with circled cut-openings

圖4 測(cè)點(diǎn)相對(duì)位置分布示意圖“○”—測(cè)點(diǎn)Fig.4 Sketch of relative distribution of measuring points “○”—measuring points

圖5 測(cè)量點(diǎn)流速的實(shí)驗(yàn)值與模擬值的比較Fig.5 Comparison of the simulated and measured velocities of the measuring stations

實(shí)驗(yàn)值Experimental result[4]本文模擬值Simulated value相對(duì)誤差Relative error阻力系數(shù)Drag coefficient1.6551.7566.09%

3 不同開口比礁體水動(dòng)力特性研究

為了減小數(shù)值模擬過程中由于網(wǎng)格尺度所產(chǎn)生的誤差，以邊長(zhǎng)為3 m、開口比為0.2的方型礁體為例，對(duì)計(jì)算域進(jìn)行不同尺寸網(wǎng)格劃分，網(wǎng)格收斂性驗(yàn)證以礁體阻力系數(shù)作為變量，結(jié)果如表2所示。

表2 不同網(wǎng)格尺寸模擬結(jié)果Table 2 Simulation results with different gird sizes

由表2可以看出邊長(zhǎng)為3 m方型礁體的計(jì)算域最大網(wǎng)格尺寸為0.500 m時(shí)，礁體阻力系數(shù)受網(wǎng)格劃分尺寸影響較小，網(wǎng)格收斂性較好，因此本文進(jìn)行模擬時(shí)最大網(wǎng)格尺寸均設(shè)置為0.500 m，礁體表面第一層邊界層網(wǎng)格高度為0.002 m，增長(zhǎng)率為1.1，共設(shè)置10層。

基于上述數(shù)值模型，通過數(shù)值模擬研究邊長(zhǎng)為3 m、開口比為0～0.6之間7種礁體流場(chǎng)效應(yīng)、阻力系數(shù)隨開口比變化情況，計(jì)算不同開口比礁體抗滑移及抗傾覆安全系數(shù)。

3.1 礁體流場(chǎng)效應(yīng)分析

本研究中上升流區(qū)含義采用了黃遠(yuǎn)東等[9]提出的定義，即礁體附近豎直方向速度分量大于或等于5%來流速度的區(qū)域。上升流區(qū)范圍越大，礁體的流場(chǎng)效應(yīng)越顯著，集魚效果越好。

圖6表示的是來流速度為0.8 m/s時(shí)，邊長(zhǎng)為3 m的不同開口比礁體在y=1.5 m截面上的速度矢量分布。由圖可以看出當(dāng)開口比較小時(shí)，由于礁體的阻水作用，水流在礁體背流面后端形成速度很小、范圍較大的漩渦，稱為流向渦[10]，此區(qū)域稱為背渦區(qū)。當(dāng)礁體開口比為0時(shí)，流向渦長(zhǎng)度約為礁高的2.3倍，這與李曉磊等[10]得出的流向渦長(zhǎng)度約為礁高的2.42倍相近，同時(shí)也與劉同渝[11]得出的流向渦長(zhǎng)度為礁長(zhǎng)的2～3倍相吻合；當(dāng)礁體開口比為0.1時(shí)，流向渦長(zhǎng)度略大，約為礁寬的3.7倍；當(dāng)開口比大于0.3時(shí)，由于礁體開口透水作用增大，流經(jīng)礁體中心至礁體后側(cè)水流流量增大，因而觀察不到明顯的背渦區(qū)。

圖6 不同開口比礁體y=1.5 m平面速度矢量分布Fig.6 The velocity vector diagrams of reefs with different opening ratios in the plane of y=1.5 m

φHmax/HWmax/WZmax/V02.573.240.410.12.753.410.370.22.332.820.370.32.122.420.370.41.912.080.290.51.781.930.270.61.611.660.27

表3中Hmax/H表示礁體產(chǎn)生的上升流最大高度/礁高，Wmax/W表示上升流水平跨度/礁寬，Zmax/V表示上升流區(qū)豎直方向最大速度分量/來流速度。與開口比為0.1時(shí)相比，礁體開口比為0時(shí)產(chǎn)生的上升流相對(duì)高度及相對(duì)寬度均略小，這與黃遠(yuǎn)東等[12]得到的規(guī)律相同。當(dāng)開口比大于0.1時(shí)，隨著開口比的增大，來流受到礁體開口的分流增多，礁體上側(cè)流管效應(yīng)減弱，魚礁引起的上升流最大高度、水平跨度及豎直方向最大速度分量均逐漸減小。

礁體阻力系數(shù)是表征人工魚礁體穩(wěn)定性的重要參數(shù)。已知礁體的阻力系數(shù)，其在水下抗滑移、抗傾覆安全系數(shù)便可通過計(jì)算求得。本文數(shù)值模擬中僅考慮礁體在來流速度不變時(shí)的受力情況，阻力系數(shù)可通過下式求得：

(5)

式中：F為礁體沿水流方向受力(N)；ρ為海水密度(kg/m3)；A為礁體迎流面積(m2)；u為水流速度(m/s)。

圖7 阻力系數(shù)與開口比的關(guān)系Fig.7 The relation between the drag coefficient and opening ratio

圖7給出了礁體阻力系數(shù)隨其開口比的變化關(guān)系。在礁體開口比為0時(shí)，其阻力系數(shù)為1.052，與Jinho Woo等[2]通過數(shù)值模擬得到的1.069及Fox等[2]得到的1.05差異較小，與White和Young等[13]得到的立方型三維實(shí)體結(jié)構(gòu)在Re≥104、垂直于水流放置時(shí)阻力系數(shù)為1.05～1.07相吻合。通過最小二乘法擬合得到邊長(zhǎng)為3 m、迎流面中心開口為方形的礁體阻力系數(shù)Cd與開口比φ的關(guān)系式(R2=0.963，P<0.01)為：

Cd=0.875φ+1.088。

(6)

可以看出在顯著性水平α<0.01時(shí)，該回歸方程回歸效果極顯著。方型人工魚礁開口比小于0.6時(shí)，隨著礁體開口比的增大，礁體阻力系數(shù)逐漸增大。

為了研究關(guān)系式(6)對(duì)于方型不同開口形狀礁體阻力系數(shù)計(jì)算的適用性，選擇圓形開口礁體進(jìn)行數(shù)值模擬，礁體尺寸為10 cm×10 cm×10 cm，中空無底，壁厚為1.2 cm，每個(gè)側(cè)面均開有直徑為3.5 cm的圓孔，如圖8所示。將本文模擬結(jié)果與姜昭陽的實(shí)驗(yàn)結(jié)果[14]進(jìn)行對(duì)比，如表4所示。

通過比較礁體阻力系數(shù)實(shí)測(cè)值與關(guān)系式計(jì)算值可以看出二者之間存在差異。利用實(shí)測(cè)阻力計(jì)算礁體阻力系數(shù)時(shí)，流速采用的是流速儀測(cè)量的平均值，與實(shí)驗(yàn)采集測(cè)力值時(shí)對(duì)應(yīng)的瞬時(shí)流速不同；此外，由于水流的作用，在測(cè)力鋼桿的后方會(huì)產(chǎn)生小的漩渦，引起鋼桿的輕微振動(dòng)，導(dǎo)致采集的測(cè)力值有波動(dòng)，使得礁體受力測(cè)量值與計(jì)算值有差異。

圖8 礁體模型尺寸Fig.8 Dimensions of the reef model

圖9 圓形開口礁體阻力系數(shù)與開口比的關(guān)系Fig.9 The relation between the drag coefficient and opening ratio of circular opening reef

阻力Resistance/N阻力系數(shù)Drag coefficient實(shí)測(cè)值[14]Experimental result[14]0.2151.056模擬值[14]Simulated value[14]0.2741.345本文模擬值Simulated value0.2411.194計(jì)算值Calculated value0.2391.172計(jì)算值與實(shí)測(cè)值相對(duì)誤差Relative error/%11.1610.98

如果考慮阻力值的數(shù)量級(jí)，關(guān)系式計(jì)算值與實(shí)測(cè)值結(jié)果較為一致，說明阻力系數(shù)與開口比的關(guān)系式可以應(yīng)用于無底圓形開口礁體阻力系數(shù)的計(jì)算。

3.2 礁體穩(wěn)定性計(jì)算

日本學(xué)者中村充[5]研究了人工魚礁體在波流共同作用時(shí)的受力情況，將波流速度看成是海流速度u0和波浪速度u1的疊加。

賈曉平等[15]認(rèn)為人工魚礁體在流速u下受到的力F可根據(jù)Morison方程計(jì)算：

(7)

式中Cd為阻力系數(shù)，CMA為附加質(zhì)量系數(shù)，A為礁體迎流面積(m2)，V為礁體實(shí)體體積(m3)。

x4+2αx3+(α2+β2-1)x2-2αx-α2=0。

(8)

(9)

通過編寫MATLAB程序?qū)崿F(xiàn)牛頓迭代法求解式(3.3)，可求得F最大時(shí)的sinθ和cosθ值，進(jìn)而求得魚礁體受力最大值Fmax。

本文選擇日照近海某人工魚礁區(qū)波浪參數(shù)[16]進(jìn)行計(jì)算，見表5。

表5 日照近海某人工魚礁區(qū)水文資料Table 5 Hydrological data of an offshore artificial reef region in Rizhao

3.2.1 礁體不滑移的安全性校核 人工魚礁體投放于海底，不發(fā)生滑動(dòng)的條件為礁體所受最大靜摩擦力大于波流作用力，即抗滑移安全系數(shù)SF需滿足下式：

(10)

式中：σ為礁體材料的單位體積重量，對(duì)于混凝土人工魚礁來說一般取2 000 kg/m3[17]；μ為礁體與海底的靜摩擦系數(shù)，取為0.5[17]。當(dāng)SF大于1時(shí)，礁體不會(huì)發(fā)生滑移。為安全起見，SF應(yīng)取1.2以上[18]。計(jì)算可得此海區(qū)波幅um為0.507 m/s，進(jìn)而算得α為1.58。

表6 不同開口比礁體受力及抗滑移安全系數(shù)Table 6 The stress and safety factor against slippage of reefs with different opening ratios

從表6中可以看出，當(dāng)開口比較小時(shí)，方型礁體所受的最大波流作用力中，速度力所占比重較小，對(duì)波流作用力產(chǎn)生較大影響的是礁體所受的慣性力。隨著開口比的增大，速度力占比逐漸增大，礁體所受最大波流作用力及抗滑移安全系數(shù)均逐漸減小，但礁體不會(huì)發(fā)生滑移，滿足安全性要求。

(11)

礁體抗傾覆安全系數(shù)為：

(12)

式中：lw為傾覆中心到礁體重心的水平距離。當(dāng)SF大于1時(shí)，礁體不會(huì)發(fā)生傾覆。為安全起見，SF應(yīng)大于1.2。

圖10 方型人工魚礁體迎流面分塊示意圖Fig.10 Schematic diagram of incident flow area of cube artificial reef

φq/kN·m-2M1/kNM2/kNSF04.001386.97854.0197.160.13.520295.13842.7596.900.23.103224.33033.5306.690.32.745166.47925.9816.410.42.415117.92819.5266.040.52.19980.35714.8645.410.61.93850.63310.5144.82

由表7可以看出隨著礁體開口比的增大，礁體迎流面均布荷載、波流最大作用力矩及抗傾覆安全系數(shù)均逐漸減小，但在日照近海人工魚礁區(qū)海浪條件下，礁體不會(huì)發(fā)生傾覆，滿足安全性要求，可以投入實(shí)際應(yīng)用。

4 結(jié)論

本文利用Fluent軟件中的κ-ε湍流模型，模擬了不同開口立方型魚礁體周圍水流場(chǎng)，通過分析其流場(chǎng)效應(yīng)、阻力系數(shù)、礁體穩(wěn)定性隨開口比變化的差異，得到了如下結(jié)論：

(1)在開口比變化的情況下，人工魚礁體流場(chǎng)效應(yīng)等方面的水動(dòng)力特性差異顯著，而礁體穩(wěn)定性的差異較小。

(2)當(dāng)開口比為0～0.2時(shí)，礁體后方流向渦和背渦區(qū)范圍較大，礁體流場(chǎng)效應(yīng)較明顯，開口比為0.3～0.6時(shí)，礁體上側(cè)流管效應(yīng)較弱，后方?jīng)]有明顯的背渦區(qū)；礁體產(chǎn)生的上升流高度、水平跨度及豎直方向最大速度分量均隨著礁體開口比的增大而減小。

(3)當(dāng)礁體開口比在0～0.6之間時(shí)，六個(gè)面中心均布有方形開口的立方型礁體的阻力系數(shù)與開口比間存在如下關(guān)系(R2=0.963，P<0.01)：Cd=0.875φ+1.088。

(4)當(dāng)開口比變化時(shí)，礁體穩(wěn)定性安全系數(shù)有一定的差異。對(duì)于方形開口礁體，隨著開口比的增大，速度力占比逐漸增大，礁體所受最大波流作用力、波流最大作用力矩、抗滑移和抗傾覆安全系數(shù)均逐漸減小，但礁體仍處于穩(wěn)定狀態(tài)。