于定勇， 楊遠航， 李宇佳

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

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

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

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

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

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

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

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

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

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

1 模型

1.1 控制方程

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

連續(xù)方程：

(1)

動量方程：

(2)

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

1.2 湍流模型

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

湍動能κ方程：

(3)

湍流耗散率ε方程：

(4)

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

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

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

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

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

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

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

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

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

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

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

2 模型可靠性驗證

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

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

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

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

圖4 測點相對位置分布示意圖“○”—測點Fig.4 Sketch of relative distribution of measuring points “○”—measuring points

圖5 測量點流速的實驗值與模擬值的比較Fig.5 Comparison of the simulated and measured velocities of the measuring stations

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

3 不同開口比礁體水動力特性研究

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

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

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

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

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

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

圖6表示的是來流速度為0.8 m/s時，邊長為3 m的不同開口比礁體在y=1.5 m截面上的速度矢量分布。由圖可以看出當開口比較小時，由于礁體的阻水作用，水流在礁體背流面后端形成速度很小、范圍較大的漩渦，稱為流向渦[10]，此區(qū)域稱為背渦區(qū)。當礁體開口比為0時，流向渦長度約為礁高的2.3倍，這與李曉磊等[10]得出的流向渦長度約為礁高的2.42倍相近，同時也與劉同渝[11]得出的流向渦長度為礁長的2～3倍相吻合；當礁體開口比為0.1時，流向渦長度略大，約為礁寬的3.7倍；當開口比大于0.3時，由于礁體開口透水作用增大，流經(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時相比，礁體開口比為0時產(chǎn)生的上升流相對高度及相對寬度均略小，這與黃遠東等[12]得到的規(guī)律相同。當開口比大于0.1時，隨著開口比的增大，來流受到礁體開口的分流增多，礁體上側(cè)流管效應(yīng)減弱，魚礁引起的上升流最大高度、水平跨度及豎直方向最大速度分量均逐漸減小。

礁體阻力系數(shù)是表征人工魚礁體穩(wěn)定性的重要參數(shù)。已知礁體的阻力系數(shù)，其在水下抗滑移、抗傾覆安全系數(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ù)為1.052，與Jinho Woo等[2]通過數(shù)值模擬得到的1.069及Fox等[2]得到的1.05差異較小，與White和Young等[13]得到的立方型三維實體結(jié)構(gòu)在Re≥104、垂直于水流放置時阻力系數(shù)為1.05～1.07相吻合。通過最小二乘法擬合得到邊長為3 m、迎流面中心開口為方形的礁體阻力系數(shù)Cd與開口比φ的關(guān)系式(R2=0.963，P<0.01)為：

Cd=0.875φ+1.088。

(6)

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

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

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

圖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實測值[14]Experimental result[14]0.2151.056模擬值[14]Simulated value[14]0.2741.345本文模擬值Simulated value0.2411.194計算值Calculated value0.2391.172計算值與實測值相對誤差Relative error/%11.1610.98

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

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

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

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

(7)

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

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

(8)

(9)

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

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

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

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

(10)

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

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

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

(11)

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

(12)

式中：lw為傾覆中心到礁體重心的水平距離。當SF大于1時，礁體不會發(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ū)海浪條件下，礁體不會發(fā)生傾覆，滿足安全性要求，可以投入實際應(yīng)用。

4 結(jié)論

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

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

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

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

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