基于桁架模型和多孔介質模型的柔性網衣和周圍流場單向耦合方法研究

徐子鳴，林志良, 2，謝 彬, 2

(1. 上海交通大學 船舶海洋與建筑工程學院 海洋工程國家重點實驗室，上海 200240; 2. 上海交通大學 船海工程數值試驗中心，上海 200240)

水產養殖一直是世界糧食生產的重要方式，在全球范圍內處于穩步擴張階段。由于人口密度增加、過度捕撈和不合理開發等因素影響，近海水產資源逐漸枯竭。為了尋求更優越的養殖環境，把網箱移至外海的深海養殖應運而生[1-4]。在深水養殖過程中，網衣的減流效應備受關注，不僅決定了整個網箱的水動力學特性，而且對網箱內水質環境的影響也非常顯著，因此需要更準確的水動力分析和計算方法。

網衣屬于小尺度、大變形的柔性多孔結構體，周圍的流場分布極其復雜，許多學者對此提出了不同研究方法。Lφland[5]通過Rudi等[6]的試驗數據擬合出網衣阻力和升力系數與來流速度的關系；宋偉華等[7]通過模型試驗研究了方形網箱外部的流速分布；Tsukrov等[8]和Tang等[9]通過大量試驗發現，網衣水動力特性主要依賴于雷諾數Re和堅固性兩個無量綱變量。網衣的數值模擬難度較大，通常將流場導致的網衣變形和網衣對流場的作用分別進行求解。其中網衣變形模型一般使用經驗公式，經驗公式有兩種，一種是Morison模型[10-11]，把每一根網線當作獨立的圓柱并將其作用力通過求解Morison方程計算；第二種為Screen模型[12]，將網衣劃分為多個平面網板單元并根據密實度和攻角計算網衣的阻力和升力系數。雷江濤等[13]和Cheng等[14]分別通過Morison模型和Screen模型對網衣在定流速情況下的變形和受力情況進行分析。而網衣周圍的流場分布可采用計算流體力學(CFD)方法進行求解。這種方法將網衣表示成多孔介質[15]并通過在Navier-Stokes方程中添加源項來表示網衣對流體的影響，其特點是不需要對網衣的詳細幾何形狀進行建模，因此具有計算簡便和效率較高的優點。趙云鵬等[16]結合多孔介質模型模擬網衣，探究了在不同流速、不同攻角、不同網衣片數等工況下水流作用下平面網衣周圍流場的情況。在網衣和流場的耦合模擬方面，Bi等[17-18]用軟件Fluent將多孔介質模型和集中質量模型耦合，分別模擬了平面網衣和圓形網衣在定常流速下的流場情況，和試驗對比了網衣受力。Chen和Christensen[19-20]通過在軟件OpenFOAM框架下結合多孔介質模型和集中質量模型，開發出一種新的耦合接口用以流體和結構求解器間的數據交換，模擬結果與試驗數據基本一致。但這種耦合方式需要在每個時間步內迭代至穩態，計算量很大，在實際應用中存在一定的局限性。Cheng等[21]利用code_aster和OpenFOAM進行網衣的流固耦合分析。新的耦合算法通過消除多孔介質系數的數據擬合過程并通過Screen模型計算網衣的水動力，從而達到簡化計算過程提高模擬精度的效果。

下文提出一種更簡便的計算方法，通過Screen模型求解網衣水動力，然后運用桁架模型建立其網衣節點位移的偏微分方程，在code_aster中計算網衣節點的位移值。當解得網衣變形后各節點的坐標后，再將其導入OpenFOAM中生成多孔介質區域并通過求解流體的控制方程來模擬流場情況。計算結果表明，網衣變形、流速分布和網阻力等特征參數與試驗數據對比良好。研究結果為養殖網箱周圍的流場分析提供了參考。

1 數學模型

1.1 桁架模型

1.1.1 概述

柔性網衣在桁架模型中被視為一組桿件和網格點。對于三維的柔性網，如圖1(a)所示，單根桿件包含兩個節點，并且每個節點都能在三維空間里任意運動。對于變形的柔性網，圖1(a)展示柔性網的配置情況，作用在桿件的力被均勻地分布到其所關聯的節點中，并且在每個節點上列出運動偏微分方程。這里結構求解器將根據Cheng等[14]提出的原理進行編寫。

1.1.2 節點受力

作用在每個節點上的力包括了水動力Fh、相鄰的兩個節點的張力Fs、節點重力Fg和浮力Fb。Fb和Fg是恒定的，由網衣的材料所決定。在求解中，水動力采用Screen模型公式計算，如圖1(b)所示，點A、B、C、D組成三維空間任意的四邊形，一般來說A、B、C、D不在同一個平面上，為了保證節點受力容易計算，將四邊形ABCD拆分為四個小三角形(ΔABC，ΔABD，ΔACD，ΔBCD)，以任意的一個小三角形ABC為例，水動力Fh=FD+FL，FD和FL的計算公式為：

圖1 桁架模型Fig. 1 Sketch of the truss model

(1)

(2)

其中，At是ΔABC網片的面積，ur為網片與來流的相對速度，iD,iL是阻力、升力的單位力方向矢量，分別與相對速度方向相同和垂直，對于節點A的水動力，為Fh/3。在式(1)、式(2)中的CD,CL由式(3)～(5)計算：

CD=0.04+(-0.04+Sn-1.24Sn2+13.7Sn3)cosθ

(3)

CL=(0.57Sn-3.54Sn2+10.1Sn3)sin2θ

(4)

(5)

式中：Sn為網衣密實度，dw為網衣直徑，L為目腳長度(如圖2所示)，θ為網片相對來流的迎角。

圖2 網衣模型Fig. 2 Net model

節點上的重力和浮力如圖1(c)所示，將桿的重力和浮力平均分給每個節點，以桿件AB為例，A節點的重力和浮力為:

(6)

式中：dw為網衣直徑，L為目腳長度，ρ和ρnet分別為流體密度和網的密度。

1.1.3 尾流效應

尾流效應是分析繞流的重要機制。在水動力模型Screen模型中，后網的來流速度衰減，因此尾流效應必須考慮。一般來說，尾流效應包含三種，如圖3所示，網線之間、單個網衣之間、網籠之間，網線之間的尾流效應已經包含在Screen模型之中[14], 對于單片網衣之間，若前網的速度為u，后網速度為ru，r為衰減因子，Lφland[5]給出的經驗公式為r=1-0.46CD,θ=0，CD,θ=0為迎角為0的阻力系數，網籠之間的尾流效應一般與單片網衣的經驗公式相同。

圖3 尾流效應Fig. 3 Wave effect

1.1.4 運動方程

對網衣節點位移建立偏微分方程，與Cheng等[14]研究有一樣的方程為：

(7)

式中：M為質量矩陣，K為剛度矩陣，Fh、Fg和Fb分別為節點的水動力、重力和浮力，已經在1.1.2節中敘述，x為瞬態的節點位移。圖4為由網衣節點1和2組成的桿件，在局部坐標系OS中可得：

圖4 單桿單元Fig. 4 A rod element

(8)

由于節點1和節點2在全局坐標系oxyz都是任意的運動，都具有三個自由度，轉換矩陣為：

(9)

式中：cos(x,s)，cos(y,s)，cos(z,s)分別為os與ox、oy、oz夾角的余弦值。因此，在全局坐標系oxyz下，M=M1·T,K=K1·T。

同時，對于單個節點的位移和速度推進運用HHT-α方法[14]:

(10)

(11)

(12)

其中，取α=0.3，Δt=0.01 s，在code_aster中模擬時間為10 s。

1.2 多孔介質模型

1.2.1 控制方程

Realizable k-ε模型是一種具有廣泛適用性的湍流模型，適合解決涉及強壓力梯度的流動問題。對于純水流繞流于平面網，由于網的屏蔽作用，存在很強的壓力梯度。因此，選擇Realizable k-ε湍流模型來解決這個問題。控制方程包括不可壓縮流體的連續性方程和動量守恒方程：

?·u=0

(13)

(14)

其中，t是時間，p是壓力，u是流體速度,μ和ρ分別是流體的動力黏度和密度，g是重力場。S是動量方程的源項。在多孔介質理論中，

(15)

其中，C是多孔介質系數矩陣。在局部坐標系中可表示為：

(16)

其中，Cn是法向慣性阻力系數，Ct是切向慣性阻力系數，對于矩陣C來說，可以由試驗數據通過最小二乘法[16]得到，也可以通過Chen等[19]推導的公式獲得，如果局部坐標系不與全局坐標系對齊，則需要對系數矩陣進行轉換[19]。

Realizable k-ε模型方程為：

(17)

(18)

1.2.2 多孔介質區域的生成

根據Chen和Christensen[20]的方法，對于每個小的網片，如圖1(b)所示，是由四個節點組成。但是四個節點不一定在一個平面上，因此將方形網格拆分成兩個三角形面板，相應的多孔介質區域是由質量點組成三角形延伸到厚度t0的棱柱來生成的。

圖5描述了如何通過全局網格生成多孔介質區域，圖中小正方體代表一個網格，x1,x2,x3分別代表一個三角形面板的坐標，x代表全局網格體心的坐標，d代表全局網格體心到三角形面板的距離，n是三角形面板的單位法向量，x0為全局網格體心在三角形面板投影的坐標，其中d和x0的計算為：

圖5 三維情況下流體模型與結構模型的映射圖解Fig. 5 Illustration on mapping of the element between the fluid and structure model in 3D case

(19)

d=(x-x2)·n

(20)

x0=x-d·n

(21)

如果網格坐標點滿足下列兩個條件：1)d<1/2t0；2)由點x0x1x2,x0x2x3,x0x3x1組成的三個小三角形面積S1,S2,S3之和等于x1x2x3組成的三角形面積S，即投影點x0應該在x1x2x3組成的三角形內部，則將其標記為多孔區域的網格。多孔區域由所有滿足判定條件的網格組成。

2 數值模型驗證

2.1 圓形網衣模型驗證

采用Bi等[18]所得的試驗結果進行模型驗證，試驗在波浪水槽中進行，測試網為圓形網，在周向有40目，在垂直方向有8目，網衣直徑為dw=1.2 mm，目腳長度為λ=20 mm，圓形框架底部直徑為D=0.254 m，網衣高度為H=0.15 mm，框架為m=8 g的配重，在本算例中，8 g的配重除去浮力的部分將會平均分給40個底部的網衣節點，每個節點的受力為0.001 73 N。物理模型和測量點如圖6所示，圓形網衣需要考慮容積損失，通過Zhao等[22]的方法來計算出容積損失率。由于多孔介質的厚度影響不大[15]，選擇多孔介質為20 mm的厚度，多孔介質系數Cn=9.16 m-1,Ct=5.67 m-1。

圖6 圓形網衣的數值模型和測量點的設置Fig. 6 Physical model of a circular net and general setting of the measurement points

2.2 網格劃分和計算條件

數值模型的網格劃分如圖7所示。x為來流方向，y為水平面上與x垂直的方向，z為水平面上法向方向。水槽左端定義為速度入口邊界；右端定義為出流邊界；水槽側壁和水面定義為滑移邊界；水槽底面定義為固壁邊界。網格類型為非結構化六面體網格，在多孔介質區域向外加密兩層網格。采用有限體積法(FVM)離散控制方程，時間項目采用一階歐拉隱式格式，對流項采用二階迎風差分格式，壓力—速度耦合選擇瞬態的PISO算法，計算時間步長取0.005 s，求解時長為40 s。

圖7 計算網格示例Fig. 7 Example of computational grids

3 結果與討論

3.1 收斂性研究

為了保證在OpenFOAM計算準確，需要進行網格收斂性研究，在流速U=0.242 m/s的情況下，使用三種網格來計算網衣的阻力。圖8是在運用1.2.2節所述的多孔介質生成方法，在OpenFOAM得到的多孔區域，從多孔區域看出，其中有部分網格缺失，原因是判斷條件算法的局限性和靜態網格的使用，多孔區域只能近似變形網的幾何形狀。表1是在三種不同網格下的計算結果，計算結果表明，在不同的全局網格數和多孔介質區域網格數下，計算流體力學(CFD)模擬的網衣阻力隨網格數變化很小，與試驗值基本吻合，網格收斂性良好，圖9為在三種網格下的切面y=-0.05 m的速度切面云圖，來流通過前網和后網都有明顯的速度降低，兩側的繞流增加，有明顯的速度過渡區域。

圖8 OpenFOAM與試驗結果之間凈變形比較Fig. 8 Comparison of the net deformation between OpenFOAM and the experimental result

表1 OpenFOAM中網格收斂性研究Tab. 1 Mesh convergence study in OpenFOAM

圖9 三種網格下切面y=-0.05 m的速度云圖Fig. 9 Velocity distribution of three kinds of meshes under the slice y=-0.05 m

3.2 不同來流速度通過單片網衣計算結果

3.2.1 在code_aster中生成圓形網衣的變形

在圓形網衣的數值模擬中，三種不同的來流速度(U=0.122、0.178、0.242 m/s)通過底部配重8 g的單片圓形網衣，圖10給出了code_aster計算中網衣在初始時刻和穩定狀態下的位置和形狀，與試驗結果吻合較好，隨著流速從0.122 m/s增加至0.242 m/s，圓形網衣的變形加劇，網衣的前端變形未能很好地與試驗匹配，主要是因為在試驗中，配重片是以圓環的形式存在的，而在數值計算中，是將圓環的配重平均分給了網衣節點，這導致了數值模型和試驗模型的差異。表2為網衣的容積損失率，隨著流速從0.122 m/s增加至0.242 m/s，容積損失率從20.9%增加至48.5%，這與圖10的結果是符合的。

表2 容積損失率Tab. 2 Volume loss rate

圖10 數值模擬和模型試驗在不同流速下圓形網衣的變形對比Fig. 10 Deformation comparison of the circular net under different current velocities by the numerical simulations and the physical model tests

3.2.2 在OpenFOAM中流場情況

圖11為速度U=0.242 m/s的水流通過圓形網衣的數值模擬結果。從圖11(a)來看，網衣內部有一段區域流速衰減，在網的后方有相當長的流速衰減區域，流速衰減最小達到了0.12 m/s，出現在網兩側的正后方，繞流速度最大值達到了0.26 m/s，出現區域為網的兩側開外區域。 從圖11(b)來看，網的內部和后方流速衰減明顯，后網的后方是一片很大的流速衰減區域，網的底部繞流明顯，速度分層明顯。從圖11整體來看，網后方的流速降低主要是由于網的阻塞作用，隨著x方向距離柔性網的距離增加，兩側和底部均有繞流，網正后方速度有明顯衰減。由于使用了Realizeable k-ε湍流模型，在速度的過渡層區域，急劇的速度梯度被減小，不同速度層之間的動量交換，速度分布比較平滑。

圖11 速度0.242 m/s數值模擬計算網周圍的流速分布Fig. 11 Flow-velocity distribution around a net calculated by numerical simulation for an incoming velocity of 0.242 m/s

3.2.3 OpenFOAM和code_aster受力對比

圖12為在兩個求解器中計算的阻力與試驗數據互相對比，從數值來看吻合較好。在圓形網衣的計算中，結構求解器code_aster的后網來流速度是根據前網的速度衰減[14]估算得到的，將流速衰減因子r看成常數是不合理的，這與Bi等[23]的測量結果不符，因此在code_aster中計算的受力與試驗數據相比存在偏差且最大誤差為15.3%，出現在來流速度U=0.178 m/s。OpenFOAM將網衣節點的坐標位置導入計算域并采用多孔介質模型計算出流場分布和受力情況，大小與試驗值較為符合。

圖12 不同流速下數值模擬和試驗數據之間的阻力比較Fig. 12 Comparison of the drag force between numerical simulations and experimental data at different current velocities

3.3 固定來流速度通過不同配重的單片網衣計算結果

3.3.1 在code_aster中生成圓形網衣的變形

對具有不同配重m(m=8、45、367 g)的網衣進行數值模擬，來流速度均為U=0.242 m/s，圖13表示了code_aster的計算結果和試驗結果對比，給出網衣的初始狀態和變形狀態，數值模擬結果和試驗結果吻合較好，隨著配重從m=8 g增加至m=367 g，圓形網衣的變形減小，在w=367 g的變形十分小。表3為網衣的容積損失率，隨著配重從m=8 g增加至m=367 g，容積損失率從48.5%減小至2.0%，這與圖13所得的結果相符合。

表3 容積損失率Tab. 3 Volume loss rate

圖13 數值模擬和模型試驗在同流速不同配重下圓形網衣的變形對比Fig. 13 Deformation comparison of the circular net under the same current velocity and different weights by the numerical simulations and the physical model tests

3.3.2 在OpenFOAM中流場情況

圖14為OpenFOAM中的模擬結果，結果表明，網的內部和后方都有明顯的流速衰減，底部有明顯的繞流，過渡區域明顯。相較于配重大的網衣，配重m=8 g的網衣后方流速衰減更大，這表明，網衣的屏蔽作用隨著配重的增加而降低，尤其在網衣的正后方。

圖14 不同配重下數值計算切面y=0的流速分布Fig. 14 Velocity distribution of slice y=0 simulated by numerical simulation in different weights

圖15為沿著圖6監測點的線中平均流速的分布，可以看出，流速與試驗數據大小相近，變化趨勢一致，最大的相對誤差為 5.6%，模擬結果較好。在配重m=8 g的圓形網衣模擬中，在x=1.25 m左右流速趨向穩定，而在配重m=45 g和m=367 g的情況下，在x=0.75 m左右流速趨向穩定，配重m=8 g的穩定流速要比配重m=45 g和m=367 g更小。

圖15 數值模擬和試驗中測量點的平均流速的比較Fig. 15 Comparison of the average flow-velocity between numerical simulations and experimental measurement points

3.3.3 OpenFOAM和code_aster受力對比

由于在模型試驗中并未提供試驗數據，本文只進行兩個求解器中網衣阻力的對比，圖16為兩個求解器的阻力對比，在OpenFOAM中的計算結果要比code_aster中大，這與圖12的結果較為一致。隨著配重的增加，阻力也持續加大。

圖16 不同流速下數值模擬和試驗數據之間的阻力比較Fig. 16 Comparison of the drag force between numerical simulations and experimental data at different current velocities

4 結 語

利用結構求解器code_aster和流體求解器OpenFOAM提出了一種更簡便的柔性網衣和流場耦合計算模型，通過圓形網衣進行了驗證。具體做法是，首先采用Screen模型計算網衣的水動力，通過桁架模型code_aster計算網衣的結構變形，模擬結果與試驗相符，然后導出網衣節點坐標至OpenFOAM中，根據判斷條件在流場的網格區域內生成多孔區域，通過加多孔介質模型源項的方法，利用OpenFOAM在不同來流和不同配重情況下進行了流場分析，網衣后方有明顯的流速衰減，兩側和底部繞流明顯，流速衰減與試驗數據趨勢一致，大小相近，說明多孔介質模型的有效性，最后比較了試驗所得網衣阻力與code_aster、OpenFOAM計算的網衣阻力，誤差均在合理范圍之間。上述結果均顯示，與Bi等[18]的模擬結果和試驗數據吻合，說明了桁架模型和多孔介質模型耦合方法的準確性。

在柔性網衣與流場的耦合計算中，時間開銷主要由流體求解器承擔，結構求解器占比較小。由于本模型采用的是單向耦合方式，因此不必考慮兩個求解器的數據交換和時間同步。其次，Screen水動力模型是由經驗公式得到，在流速不高的條件下求解網衣變形是比較可靠的。因此本文提出的基于桁架模型和多孔介質模型建立的柔性網衣和流場的單向耦合方法對于實際問題的研究是經濟、有效的。