陳 誠，宋煒，謝正麗,王磊，劉永利

(1 中國水產(chǎn)科學研究院東海水產(chǎn)研究所，上海 200090；2 上海交通大學海洋工程國家重點實驗室，上海 200240；3 中國水產(chǎn)科學研究院漁業(yè)機械儀器研究所，上海 200090)

近年來，海上養(yǎng)殖業(yè)逐漸由近海向深遠海發(fā)展，但深遠海海域的海況復雜惡劣且難以預測，在極端條件下會對結構物產(chǎn)生巨大的載荷作用。網(wǎng)衣作為養(yǎng)殖裝備的重要組成部分，對其水動力特性的研究對于了解養(yǎng)殖裝備的水動力響應具有重要意義。

研究人員對網(wǎng)衣的水動力特性進行了大量的實驗探究，如Zhan等[1]，李玉成等[2]，Lader等[3-4]，宋偉華等[5]，Patursson等[6]，這些研究主要聚焦于雷諾數(shù)、密實度、來流攻角等參數(shù)對于網(wǎng)衣阻力和升力的影響，并根據(jù)實驗結果擬合出一系列網(wǎng)衣阻力和升力的經(jīng)驗公式。此外，Kristiansen等[7]〗開展了風浪下網(wǎng)箱系泊特性的研究，Bi等[8]，曾啟東等[9]從波浪和網(wǎng)衣參數(shù)的角度對網(wǎng)箱的阻尼特性開展了研究。

隨著CFD的發(fā)展，多種數(shù)值模型逐漸發(fā)展并被應用于網(wǎng)衣水動力的研究上。Van等[10]最早發(fā)展了基于Navier-Stokes方程的多孔介質(zhì)模型，并將其應用于近岸多孔結構的研究上，此后，多孔介質(zhì)模型在其他學者如Nakayama等[11]、趙云鵬等[12-13]、Jensen等[14]、Bi等[15]的研究中不斷完善與擴展。與此同時，Morison模型和Screen模型被應用于網(wǎng)衣水動力模型的構建，Tsukrov等[16]構建了基于修正Morison方程的網(wǎng)衣水動力的有限元模型，Chen等[17]提出了一種基于Morison載荷模型的多孔介質(zhì)阻力系數(shù)計算方法，Kristiansen等[18]提出并討論了粘性水動力的Screen模型。此外，集中質(zhì)量法(Huang等[19]，Martin等[20])、有限體積法(Zhao等[21-22])、流體體積法(Jafari等[23])、有限元法(崔勇等[24])等也被廣泛應用于網(wǎng)衣水動力研究中。在數(shù)值計算的載體上，由Jensen等[14]，F(xiàn)eichtner等[25]在OpenFOAM軟件中開發(fā)了基于多孔介質(zhì)模型的計算框架，以此為基礎開展數(shù)值研究。

然而，在現(xiàn)有的網(wǎng)衣水動力研究中，絕大多數(shù)以規(guī)則波作為輸入，無法體現(xiàn)出深遠海中所可能面對的極端海況。為此，需要采用一些有代表性的不規(guī)則波來模擬極端海況。Xu等[26]在波浪水槽中進行了一系列極端波的實驗，得出密實度和KC數(shù)是影響網(wǎng)板作用力的主要參數(shù)。俞嘉臻等[27]將聚焦波與網(wǎng)衣相互作用，分析了不同網(wǎng)衣參數(shù)對升阻力特性的影響及對波浪場的擾動規(guī)律。新波理論(NewWave)從概率統(tǒng)計的意義上代表了波峰附近最可能發(fā)生的波浪時歷，既包含了波浪譜的特征信息，又能有效減少數(shù)值模擬的時間[28]。

針對波浪與網(wǎng)衣結構的相互作用展開研究，將多孔介質(zhì)模型引入網(wǎng)衣系統(tǒng)，使用Morison模型對網(wǎng)衣進行力學分析，并通過等效分析得到了多孔介質(zhì)阻力系數(shù)的直接估計方法。數(shù)值模擬過程在開源CFD軟件OpenFOAM中進行，并使用造波工具箱waves2Foam實現(xiàn)規(guī)則波和聚焦波的模擬?；跀?shù)值模擬的結果，通過參數(shù)化分析研究了影響網(wǎng)衣所受阻力的因素以及網(wǎng)衣對波浪場的影響。最后，為尋求一種對極端波浪下網(wǎng)衣多孔介質(zhì)阻力的有效預報，根據(jù)已有的數(shù)據(jù)總結給出了網(wǎng)衣水動力載荷的預測模型?；诒疚牡难芯靠梢詫W(wǎng)衣的設計與裝配提供理論參考。

1 網(wǎng)衣水動力模型

1.1 多孔介質(zhì)模型

本數(shù)值模擬研究在開源CFD軟件OpenFOAM中展開，采用了有限體積離散化的方法，求解了體積平均、雷諾平均的N-S方程(即VARANS方程)。Jensen等[14]給出了詳細的推導過程，其最終的連續(xù)性方程和動量方程如下：

(1)

(2)

(3)

式中：cij為非線性阻力系數(shù)矩陣。

流體對網(wǎng)衣的作用力即為(3)式所表征流體所受阻力的反作用力，其合力可由每個微元的受力的積分得到：

(4)

式中：PV為多孔介質(zhì)區(qū)域的體積，cij是未知的，將在1.2節(jié)中使用Morison載荷模型進行等效得到。

1.2 Morison模型

在Morison模型中，將每一根網(wǎng)線視為一個獨立的圓柱體進行分析，不考慮結節(jié)的影響，并利用Morison方程來計算其上所受的作用力。由于網(wǎng)衣的孔隙率較高，在研究中忽略了網(wǎng)線間的水動力影響。因而，作用在整個網(wǎng)面上的合力可由每根網(wǎng)線上的力疊加得到。

如圖1所示，假設網(wǎng)衣在y-z平面內(nèi)，網(wǎng)面的法線沿x方向，網(wǎng)目的基本單元由橫向的網(wǎng)線1及垂向的網(wǎng)線2組成。

注：U∞為來流速度，Uxy為來流速度在xoy平面內(nèi)的分量，Uz為來流速度在z方向上的分量

來流速度為U∞，其分量Uxy=U∞cosβ，Uz=U∞sinβ。x-y平面內(nèi)的來流與網(wǎng)衣相對方位圖如圖2所示。

注：Uxy為來流速度在xoy平面內(nèi)的分量，F(xiàn)xy為波浪對網(wǎng)線的作用力

對于網(wǎng)線1，流速Uxy可分解為切向分量和法向分量，其中切線分量引起的摩擦力是一個小量，可忽略不計，流速Uz引起一個z方向的力，因此網(wǎng)線1主要受法向力作用，則其所受力大小為：

(5)

Fy,1=0

(6)

(7)

對于網(wǎng)線2，流速Uxy沿網(wǎng)線的法向，流速Uz與網(wǎng)線平行，則其所受力大小為：

(8)

(9)

Fz,2=0

(10)

式中：Cd為網(wǎng)線的阻力系數(shù)，A1,A2等于網(wǎng)線1在y-z平面內(nèi)的總投影面積。

則作用在網(wǎng)衣上的合力就是每根網(wǎng)線受力的疊加

(11)

(12)

(13)

式中：M,N分別為平行于網(wǎng)線1、2的網(wǎng)線根數(shù)。

又由1.1節(jié)得，作用在多孔介質(zhì)結構上的力為：

(14)

(15)

(16)

式中：c1,c2,c3分別為x,y,z方向上的阻力系數(shù)，V為網(wǎng)衣的體積，網(wǎng)衣是一種高孔隙率的多孔結構，波浪流經(jīng)網(wǎng)衣前后，流速的降低并不明顯，因此可以認為網(wǎng)衣附近的流速Uxy等于遠處的來流速度U∞。

由于作用在網(wǎng)衣上的力F等于每根網(wǎng)線受力Q的疊加，將二者列為等式，即可得到阻力系數(shù)C1,C2,C3的表達式：

(17)

(18)

(19)

則基于多孔介質(zhì)模型的網(wǎng)衣受力如下：

(20)

(21)

(22)

Chen等[17]通過對前人實驗結果的總結，給出了在一定密實度范圍內(nèi)，系數(shù)a和b的參考值如下：

(23)

(24)

在本研究中，波浪沿x軸正方向垂直于網(wǎng)衣入射，即攻角θ=90°，則網(wǎng)衣受到沿x方向的阻力Qx以及z方向的升力Qz，y方向受力為0。

1.3 阻力系數(shù)等效選取

Lader等[3]通過多組網(wǎng)衣實驗得到了不同波浪與網(wǎng)衣相互作用后網(wǎng)衣阻力數(shù)據(jù)。Chen等[17]對其進行了整理，得到了對應的網(wǎng)線阻力系數(shù)cd的數(shù)據(jù)。圓柱形網(wǎng)線的阻力系數(shù)受波浪的流態(tài)影響，可以認為主要與雷諾數(shù)Re和KC數(shù)有關，二者的表達式如下：

(25)

(26)

式中：U為波浪的特征速度U=Aω，即波幅A與圓頻率ω之積，v為運動粘度，T為波浪周期。

為此，可將網(wǎng)線的阻力系數(shù)視為Re和KC的函數(shù)，用多參數(shù)擬合的方法擬合數(shù)據(jù)，并用于內(nèi)插法取值，擬合結果如圖3所示。

注：cd為網(wǎng)線阻力系數(shù)，Re為雷諾數(shù)，KC為KC數(shù)

2 數(shù)值模擬方案

2.1 數(shù)值水池

二維數(shù)值水池的參數(shù)根據(jù)波浪和網(wǎng)衣的參數(shù)來決定。如圖4所示，其中水池總長為5.5λ，水深1λ，滿足深水條件，其中λ為波長。左側邊界為造波板，所產(chǎn)生的波浪向x軸正方向傳播，在4～5.5λ的范圍內(nèi)設置了消波區(qū)，用于吸收波浪，減少波浪反射的影響。

圖4 數(shù)值水池示意圖

網(wǎng)衣高1 m，厚0.05 m，被豎直放置在距造波板2λ的位置，一半被浸沒在水中，在數(shù)值模擬中測量網(wǎng)衣前后1λ范圍內(nèi)的波浪時歷以及網(wǎng)衣處的多孔介質(zhì)阻力及升力。

1.5 模型有效性驗證

為驗證所建立模型的有效性，在OpenFOAM中建立與Lader等[3]所做試驗相同的水池模型，使用相同的波浪、網(wǎng)衣參數(shù)進行數(shù)值模擬，對比分析波浪時歷η、網(wǎng)衣處的阻力Qx、升力Qz等數(shù)據(jù)，結果如下。

由圖5可知，數(shù)值模擬的結果與試驗結果比較接近，阻力誤差的主要來源可能在于試驗中網(wǎng)衣的變形等不確定性在數(shù)值模擬中并未被考慮，但由結果來看模型的準確性已基本滿足要求。

注：周期T=0.70 s，波幅A=0.022 m

3 數(shù)值模擬研究

3.1 規(guī)則波數(shù)值模擬

3.1.1 密實度的影響

密實度Sn描述著網(wǎng)線投影面積占網(wǎng)衣總面積的比值，在海洋環(huán)境中，隨著時間的推移，各種藻類、砂礫等會附著在網(wǎng)衣上，這將改變網(wǎng)衣原有的密實度，對網(wǎng)衣的水動力特性產(chǎn)生影響，因而研究密實度對網(wǎng)衣載荷的影響有重要意義。密實度主要從兩個方面來影響波浪的傳播。一方面，隨著密實度的增加，漁網(wǎng)在波浪方向上的有效面積增大，從而增加了對波浪傳播的阻礙作用；另一方面，隨著密實度的增加，波浪與網(wǎng)衣的相互作用面積增大，因此，摩擦所引起的波浪能量耗散也越大。

采用同一波浪與四種不同密實度的網(wǎng)衣相互作用，數(shù)值模擬參數(shù)和結果見下。

表1 網(wǎng)衣參數(shù)

注：周期T=0.80 s，波幅A=0.02 m

由圖7可知，在一定的密實度范圍內(nèi)，網(wǎng)衣所受的阻力隨著密實度的增大而增大，二者近似呈線性關系。當密實度大到一定程度，網(wǎng)衣阻力與密實度之間顯示出非線性關系。

圖7 阻力Qx與密實度Sn關系圖

3.1.2 特征速度的影響

流速是影響網(wǎng)衣阻力的重要參數(shù)，然而，在實際情況下網(wǎng)衣各處的流速存在差異且難以捕捉，為此，需要找出一個具有代表性的速度來進行網(wǎng)衣的阻尼特性研究。

流體力學中的很多無量綱常數(shù)計算中都包含著特征速度U=Aω，即波幅A與波浪圓頻率ω的積，它能夠從幅值和頻域兩個角度來表示流場的信息，并且不受結構物的影響，容易得到。為研究特征速度對阻力的影響，使用相同圓頻率、不同波高的規(guī)則波與網(wǎng)衣相互作用，數(shù)值模擬的參數(shù)如表2所示，測得不同特征速度下的阻力時歷如圖8所示。

注：網(wǎng)目尺寸D=0.016 m，網(wǎng)線直徑dω=0.001 8 m，密實度Sn=0.225

表2 波浪參數(shù)

由圖9可知，隨著特征速度的增大，網(wǎng)衣的阻力逐漸增大。當波長一定時，隨著波高的增加，水粒子的速度增大，從而導致了粘性阻力的增大，波浪力的增加表明了波浪與網(wǎng)衣的相互作用更為劇烈，能量耗散得越快，在3.1.3節(jié)中將會對波浪透射系數(shù)的影響因子做進一步探究；此外，隨著波高的增加，網(wǎng)面對波浪的有效面積增大，也會使網(wǎng)衣對波浪的阻尼作用增強。

圖9 阻力Qx與特征速度U關系圖

3.1.3 網(wǎng)后波浪透射系數(shù)

網(wǎng)衣對波浪的阻尼作用會導致波能的部分消耗，由水波理論可知，諧波的能量正比于波幅的平方，波幅的大小體現(xiàn)了能量的大小。定義了一個波浪透射系數(shù)β[29]，其數(shù)值等于規(guī)則波與網(wǎng)衣作用穩(wěn)定后網(wǎng)衣后1λ處與無網(wǎng)衣時同一位置處的平均幅值之比。在流體與結構物的相互作用中，波長與結構物特征長度的比值具有重要影響，當波長與特征長度差不多時，波浪流經(jīng)物體主要體現(xiàn)為繞射；當特征長度遠小于波長時，結構物的存在對波浪的影響很小。對于網(wǎng)衣來說，網(wǎng)線的直徑dW為特征長度，現(xiàn)對波浪透射系數(shù)同網(wǎng)繩直徑與波長的比值dW/λ的關系進行研究，數(shù)值模擬的參數(shù)如表3所示，測得不同dW/λ下的波浪時歷如圖10所示。

注：網(wǎng)目尺寸D=0.016 m，網(wǎng)線直徑dω=0.002 m，密實度Sn=0.25

表3 波浪參數(shù)

由圖11可知，隨著網(wǎng)線直徑相對值的增大，透射系數(shù)有下降的趨勢，表明結構物尺度與波長越接近，對波浪的阻尼作用越強，這也與特征長度與波浪作用之間的關系相吻合。在本研究中，網(wǎng)線直徑與波長的比值均為小量，因而結果呈現(xiàn)出一定的線性。

圖11 透射系數(shù)β與網(wǎng)線直徑與波長比值dW/λ的關系

3.2 聚焦波數(shù)值模擬

3.2.1 能量譜特性研究

在自然情況下，波浪總是由多種頻率成分組成的，聚焦波就可看作是一系列不同頻率和波幅的正弦波相疊加而成，即：

(27)

式中：N為所取正弦波的個數(shù)，kj為波數(shù)，ωj為圓頻率，φj為相位角，x0為聚焦位置，t0為聚焦時刻，Aj為波幅，其表達式為：

(28)

式中：Ac為譜峰波幅，δω為所取的頻率間隔，S(ωj)為海浪譜，本研究中采用的是Jonswap譜，其表達式為:

(29)

式中：H1/3為有義波高，γ取3.3。

為研究網(wǎng)衣與波浪作用后波浪的頻域變化情況，使用聚焦波為輸入，令其在網(wǎng)衣處發(fā)生聚焦，分別測量在無網(wǎng)衣和有網(wǎng)衣情況下聚焦位置后1λ處的波浪時歷，數(shù)值模擬的參數(shù)及在有無網(wǎng)衣條件下測得的波浪時歷如圖12所示。

注：譜峰波幅Ac=0.02，譜峰波長λp=1 m，網(wǎng)目尺寸D=0.01 m，網(wǎng)線直徑dω=0.001 5 m，密實度Sn=0.30

由圖12的結果可以發(fā)現(xiàn)，二者在有網(wǎng)衣和無網(wǎng)衣時在聚焦位置后1λ處的波浪時歷幾乎相同，但在有網(wǎng)衣情況下的幅值小幅降低了3.4%左右，原因是波浪與網(wǎng)衣作用后的能量損耗，此外，在作用后的尾流也有所變化。使用FFT對兩個時歷做頻域分析，得到二者的頻域譜如圖13所示。由圖13可知，在小密實度條件下，波浪在與網(wǎng)衣作用后，波浪整體的頻域分布變化不大，但有部分能量向高頻發(fā)生了轉移。

圖13 有網(wǎng)衣和無網(wǎng)衣情況下網(wǎng)后1m處波浪時歷的FFT分析結果

3.2.2D/λp的影響

在3.1.3節(jié)中，研究了特征長度與網(wǎng)線直徑的比值對網(wǎng)后波浪透射系數(shù)的影響，但對于海洋養(yǎng)殖網(wǎng)來說，網(wǎng)目尺寸是一個更重要的參數(shù)，既要滿足一定的間距，保證充分的水交換，還要避免所養(yǎng)殖魚類的逃逸，并在此基礎上，滿足基本的安全性能。

為研究網(wǎng)目尺寸D對網(wǎng)衣阻力的影響，以聚焦波為輸入，作用在四種網(wǎng)目尺寸不同的網(wǎng)衣上，測量網(wǎng)衣處的阻力，數(shù)值模擬的參數(shù)如表4所示，測得不同D/λp下的阻力時歷如圖14所示。

注：譜峰波幅Ac=0.02 m，譜峰波長λp=1 m

表4 網(wǎng)衣參數(shù)

由圖15可知，在波長相等的情況下，隨著網(wǎng)目尺寸的增大，網(wǎng)衣的阻力逐漸減小。這意味著從網(wǎng)衣出發(fā)有兩種角度來降低載荷，一是增大網(wǎng)目尺寸，即在網(wǎng)衣大小不變的條件下減少網(wǎng)線的根數(shù)，二是在網(wǎng)線直徑不變的情況下降低網(wǎng)線直徑，但同時需要保證網(wǎng)孔的尺度不至于讓養(yǎng)殖的魚類通過，并且保證一定的直徑來保證安全強度。

圖15 阻力Qx與網(wǎng)目尺寸與譜峰波長比值D/λp關系圖

3.3 阻力預報

在現(xiàn)實海況下，想要實時獲取或預測網(wǎng)衣上所受的載荷是十分困難的，如果能夠通過一些容易得到的數(shù)據(jù)，來對網(wǎng)衣載荷進行及時的預報，或者在預設的海況下能夠事先了解網(wǎng)衣上的載荷分布情況，并為網(wǎng)衣的設計提供一定的參考，這是研究網(wǎng)衣水動力特性的主要目的。

通常計算網(wǎng)衣所受的阻力的方法是采用前人總結得到的經(jīng)驗阻力系數(shù)，改方法對于均勻穩(wěn)定的流場比較有效，但目前針對極端波浪的網(wǎng)衣阻力計算和預測方法還很少。從作用的效果上看，阻力是網(wǎng)衣上的主要載荷，會引起網(wǎng)衣的形變，也是過往研究中的重點，根據(jù)之前的研究可以知道，阻力的大小與波浪特性與網(wǎng)衣特性均有緊密的關系?？紤]來流攻角為90°的危險情況時，根據(jù)Morison模型推出的多孔介質(zhì)阻力的表達式為:

(30)

從表達式看，當來流攻角不變時，對網(wǎng)衣阻力起關鍵作用的變量有多孔阻力系數(shù)、流速、網(wǎng)線總面積，其中密度ρ、系數(shù)a和網(wǎng)線總面積S是一個定值，則只需要知道多孔阻力系數(shù)Cd和流速Uxy平方的積分。

(31)

其中U為3.1.2節(jié)中定義的特征速度U=Aω，由對應海域的波浪譜計算得到，如對于本研究中的聚焦波可選取譜峰波長所對應的流速，則簡化后的新多孔介質(zhì)阻力公式為

(32)

在2.3節(jié)的分析中，將網(wǎng)線的阻力系數(shù)cd看作Re和KC的函數(shù)。

由于U即為Re和KC表達式中的U，因此同樣可以將參數(shù)C看作是Re和KC的函數(shù)，通過之前的規(guī)則波和聚焦波的數(shù)值模擬，已經(jīng)得到了一系列參數(shù)及其對應的阻力模擬結果，由此可以計算出相應的阻力系數(shù)C，對Re和KC，同樣采用多參數(shù)擬合的方法，得到的擬合模型如圖16所示。

注：C為新阻力系數(shù)，Qx為阻力，ρ為密度，a為無因次參數(shù)，S為網(wǎng)線總面積，U為特征速度，V為多孔介質(zhì)區(qū)域體積，Re為雷諾數(shù)，KC為KC數(shù)

4 結論

本研究基于多孔介質(zhì)模型和Morison模型建立了網(wǎng)衣模型，使用OpenFOAM對網(wǎng)衣的水動力性能開展數(shù)值模擬研究。研究中發(fā)現(xiàn)，對同一規(guī)則波，隨著密實度的增加，網(wǎng)衣對波浪的阻尼作用逐漸增強，在一定范圍內(nèi)兩者近似呈線性關系；對不同的規(guī)則波，隨著特征速度的增大，網(wǎng)衣的多孔介質(zhì)阻力逐漸增大；隨著特征長度的增大，網(wǎng)后的波浪透射系數(shù)逐漸增大。在以聚焦波為輸入的研究中發(fā)現(xiàn)，在與網(wǎng)衣作用后，波浪的頻域分布沒有明顯的變化，但有部分能量向高頻區(qū)域發(fā)生了轉移；譜峰波長一定時，隨著網(wǎng)目尺寸的增大，網(wǎng)衣所受的多孔介質(zhì)阻力逐漸降低。

最后，提出了一種針對極端波浪的網(wǎng)衣多孔介質(zhì)阻力預報方法，對多孔介質(zhì)阻力公式進行了簡化，引入了波浪譜的特征參數(shù)作為輸入，并定義了新的阻力系數(shù)，將其表示為Re數(shù)和KC數(shù)的函數(shù)，結果顯示，擬合結果在取值范圍內(nèi)收斂，進而可以計算得到阻力，為網(wǎng)衣阻力的預測分析、網(wǎng)衣參數(shù)的設計提供一定的參考。

然而，在當前的研究中，網(wǎng)衣密實度的研究范圍普遍在0.1～0.3的范圍內(nèi)，對大密實度網(wǎng)衣的研究還存在著很多空白；此外，對于一個典型的水產(chǎn)養(yǎng)殖網(wǎng)來說，Re數(shù)的范圍通常在100～10 000[30]，而本研究僅限于188～471，因此還需要有更多的實驗來得到更廣泛的多孔介質(zhì)阻力數(shù)據(jù)以完善擬合模型，得到更多定量的結論；本研究中始終保持波浪傳播方向沿網(wǎng)衣的法向，在后續(xù)的研究中還可以將攻角的影響引入模型中。