徐子鳴,林志良, 2,謝 彬, 2
(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國家重點實驗室,上海 200240; 2. 上海交通大學(xué) 船海工程數(shù)值試驗中心,上海 200240)
水產(chǎn)養(yǎng)殖一直是世界糧食生產(chǎn)的重要方式,在全球范圍內(nèi)處于穩(wěn)步擴張階段。由于人口密度增加、過度捕撈和不合理開發(fā)等因素影響,近海水產(chǎn)資源逐漸枯竭。為了尋求更優(yōu)越的養(yǎng)殖環(huán)境,把網(wǎng)箱移至外海的深海養(yǎng)殖應(yīng)運而生[1-4]。在深水養(yǎng)殖過程中,網(wǎng)衣的減流效應(yīng)備受關(guān)注,不僅決定了整個網(wǎng)箱的水動力學(xué)特性,而且對網(wǎng)箱內(nèi)水質(zhì)環(huán)境的影響也非常顯著,因此需要更準(zhǔn)確的水動力分析和計算方法。
網(wǎng)衣屬于小尺度、大變形的柔性多孔結(jié)構(gòu)體,周圍的流場分布極其復(fù)雜,許多學(xué)者對此提出了不同研究方法。Lφland[5]通過Rudi等[6]的試驗數(shù)據(jù)擬合出網(wǎng)衣阻力和升力系數(shù)與來流速度的關(guān)系;宋偉華等[7]通過模型試驗研究了方形網(wǎng)箱外部的流速分布;Tsukrov等[8]和Tang等[9]通過大量試驗發(fā)現(xiàn),網(wǎng)衣水動力特性主要依賴于雷諾數(shù)Re和堅固性兩個無量綱變量。網(wǎng)衣的數(shù)值模擬難度較大,通常將流場導(dǎo)致的網(wǎng)衣變形和網(wǎng)衣對流場的作用分別進行求解。其中網(wǎng)衣變形模型一般使用經(jīng)驗公式,經(jīng)驗公式有兩種,一種是Morison模型[10-11],把每一根網(wǎng)線當(dāng)作獨立的圓柱并將其作用力通過求解Morison方程計算;第二種為Screen模型[12],將網(wǎng)衣劃分為多個平面網(wǎng)板單元并根據(jù)密實度和攻角計算網(wǎng)衣的阻力和升力系數(shù)。雷江濤等[13]和Cheng等[14]分別通過Morison模型和Screen模型對網(wǎng)衣在定流速情況下的變形和受力情況進行分析。而網(wǎng)衣周圍的流場分布可采用計算流體力學(xué)(CFD)方法進行求解。這種方法將網(wǎng)衣表示成多孔介質(zhì)[15]并通過在Navier-Stokes方程中添加源項來表示網(wǎng)衣對流體的影響,其特點是不需要對網(wǎng)衣的詳細幾何形狀進行建模,因此具有計算簡便和效率較高的優(yōu)點。趙云鵬等[16]結(jié)合多孔介質(zhì)模型模擬網(wǎng)衣,探究了在不同流速、不同攻角、不同網(wǎng)衣片數(shù)等工況下水流作用下平面網(wǎng)衣周圍流場的情況。在網(wǎng)衣和流場的耦合模擬方面,Bi等[17-18]用軟件Fluent將多孔介質(zhì)模型和集中質(zhì)量模型耦合,分別模擬了平面網(wǎng)衣和圓形網(wǎng)衣在定常流速下的流場情況,和試驗對比了網(wǎng)衣受力。Chen和Christensen[19-20]通過在軟件OpenFOAM框架下結(jié)合多孔介質(zhì)模型和集中質(zhì)量模型,開發(fā)出一種新的耦合接口用以流體和結(jié)構(gòu)求解器間的數(shù)據(jù)交換,模擬結(jié)果與試驗數(shù)據(jù)基本一致。但這種耦合方式需要在每個時間步內(nèi)迭代至穩(wěn)態(tài),計算量很大,在實際應(yīng)用中存在一定的局限性。Cheng等[21]利用code_aster和OpenFOAM進行網(wǎng)衣的流固耦合分析。新的耦合算法通過消除多孔介質(zhì)系數(shù)的數(shù)據(jù)擬合過程并通過Screen模型計算網(wǎng)衣的水動力,從而達到簡化計算過程提高模擬精度的效果。
下文提出一種更簡便的計算方法,通過Screen模型求解網(wǎng)衣水動力,然后運用桁架模型建立其網(wǎng)衣節(jié)點位移的偏微分方程,在code_aster中計算網(wǎng)衣節(jié)點的位移值。當(dāng)解得網(wǎng)衣變形后各節(jié)點的坐標(biāo)后,再將其導(dǎo)入OpenFOAM中生成多孔介質(zhì)區(qū)域并通過求解流體的控制方程來模擬流場情況。計算結(jié)果表明,網(wǎng)衣變形、流速分布和網(wǎng)阻力等特征參數(shù)與試驗數(shù)據(jù)對比良好。研究結(jié)果為養(yǎng)殖網(wǎng)箱周圍的流場分析提供了參考。
1.1.1 概述
柔性網(wǎng)衣在桁架模型中被視為一組桿件和網(wǎng)格點。對于三維的柔性網(wǎng),如圖1(a)所示,單根桿件包含兩個節(jié)點,并且每個節(jié)點都能在三維空間里任意運動。對于變形的柔性網(wǎng),圖1(a)展示柔性網(wǎng)的配置情況,作用在桿件的力被均勻地分布到其所關(guān)聯(lián)的節(jié)點中,并且在每個節(jié)點上列出運動偏微分方程。這里結(jié)構(gòu)求解器將根據(jù)Cheng等[14]提出的原理進行編寫。
1.1.2 節(jié)點受力
作用在每個節(jié)點上的力包括了水動力Fh、相鄰的兩個節(jié)點的張力Fs、節(jié)點重力Fg和浮力Fb。Fb和Fg是恒定的,由網(wǎng)衣的材料所決定。在求解中,水動力采用Screen模型公式計算,如圖1(b)所示,點A、B、C、D組成三維空間任意的四邊形,一般來說A、B、C、D不在同一個平面上,為了保證節(jié)點受力容易計算,將四邊形ABCD拆分為四個小三角形(ΔABC,ΔABD,ΔACD,ΔBCD),以任意的一個小三角形ABC為例,水動力Fh=FD+FL,F(xiàn)D和FL的計算公式為:
圖1 桁架模型Fig. 1 Sketch of the truss model
(1)
(2)
其中,At是ΔABC網(wǎng)片的面積,ur為網(wǎng)片與來流的相對速度,iD,iL是阻力、升力的單位力方向矢量,分別與相對速度方向相同和垂直,對于節(jié)點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為網(wǎng)衣密實度,dw為網(wǎng)衣直徑,L為目腳長度(如圖2所示),θ為網(wǎng)片相對來流的迎角。
圖2 網(wǎng)衣模型Fig. 2 Net model
節(jié)點上的重力和浮力如圖1(c)所示,將桿的重力和浮力平均分給每個節(jié)點,以桿件AB為例,A節(jié)點的重力和浮力為:
(6)
式中:dw為網(wǎng)衣直徑,L為目腳長度,ρ和ρnet分別為流體密度和網(wǎng)的密度。
1.1.3 尾流效應(yīng)
尾流效應(yīng)是分析繞流的重要機制。在水動力模型Screen模型中,后網(wǎng)的來流速度衰減,因此尾流效應(yīng)必須考慮。一般來說,尾流效應(yīng)包含三種,如圖3所示,網(wǎng)線之間、單個網(wǎng)衣之間、網(wǎng)籠之間,網(wǎng)線之間的尾流效應(yīng)已經(jīng)包含在Screen模型之中[14], 對于單片網(wǎng)衣之間,若前網(wǎng)的速度為u,后網(wǎng)速度為ru,r為衰減因子,Lφland[5]給出的經(jīng)驗公式為r=1-0.46CD,θ=0,CD,θ=0為迎角為0的阻力系數(shù),網(wǎng)籠之間的尾流效應(yīng)一般與單片網(wǎng)衣的經(jīng)驗公式相同。
圖3 尾流效應(yīng)Fig. 3 Wave effect
1.1.4 運動方程
對網(wǎng)衣節(jié)點位移建立偏微分方程,與Cheng等[14]研究有一樣的方程為:
(7)
式中:M為質(zhì)量矩陣,K為剛度矩陣,F(xiàn)h、Fg和Fb分別為節(jié)點的水動力、重力和浮力,已經(jīng)在1.1.2節(jié)中敘述,x為瞬態(tài)的節(jié)點位移。圖4為由網(wǎng)衣節(jié)點1和2組成的桿件,在局部坐標(biāo)系OS中可得:
圖4 單桿單元Fig. 4 A rod element
(8)
由于節(jié)點1和節(jié)點2在全局坐標(biāo)系oxyz都是任意的運動,都具有三個自由度,轉(zhuǎn)換矩陣為:
(9)
式中:cos(x,s),cos(y,s),cos(z,s)分別為os與ox、oy、oz夾角的余弦值。因此,在全局坐標(biāo)系oxyz下,M=M1·T,K=K1·T。
同時,對于單個節(jié)點的位移和速度推進運用HHT-α方法[14]:
(10)
(11)
(12)
其中,取α=0.3,Δt=0.01 s,在code_aster中模擬時間為10 s。
1.2.1 控制方程
Realizable k-ε模型是一種具有廣泛適用性的湍流模型,適合解決涉及強壓力梯度的流動問題。對于純水流繞流于平面網(wǎng),由于網(wǎng)的屏蔽作用,存在很強的壓力梯度。因此,選擇Realizable k-ε湍流模型來解決這個問題??刂品匠贪ú豢蓧嚎s流體的連續(xù)性方程和動量守恒方程:
?·u=0
(13)
(14)
其中,t是時間,p是壓力,u是流體速度,μ和ρ分別是流體的動力黏度和密度,g是重力場。S是動量方程的源項。在多孔介質(zhì)理論中,
(15)
其中,C是多孔介質(zhì)系數(shù)矩陣。在局部坐標(biāo)系中可表示為:
(16)
其中,Cn是法向慣性阻力系數(shù),Ct是切向慣性阻力系數(shù),對于矩陣C來說,可以由試驗數(shù)據(jù)通過最小二乘法[16]得到,也可以通過Chen等[19]推導(dǎo)的公式獲得,如果局部坐標(biāo)系不與全局坐標(biāo)系對齊,則需要對系數(shù)矩陣進行轉(zhuǎn)換[19]。
Realizable k-ε模型方程為:
(17)
(18)
1.2.2 多孔介質(zhì)區(qū)域的生成
根據(jù)Chen和Christensen[20]的方法,對于每個小的網(wǎng)片,如圖1(b)所示,是由四個節(jié)點組成。但是四個節(jié)點不一定在一個平面上,因此將方形網(wǎng)格拆分成兩個三角形面板,相應(yīng)的多孔介質(zhì)區(qū)域是由質(zhì)量點組成三角形延伸到厚度t0的棱柱來生成的。
圖5描述了如何通過全局網(wǎng)格生成多孔介質(zhì)區(qū)域,圖中小正方體代表一個網(wǎng)格,x1,x2,x3分別代表一個三角形面板的坐標(biāo),x代表全局網(wǎng)格體心的坐標(biāo),d代表全局網(wǎng)格體心到三角形面板的距離,n是三角形面板的單位法向量,x0為全局網(wǎng)格體心在三角形面板投影的坐標(biāo),其中d和x0的計算為:
圖5 三維情況下流體模型與結(jié)構(gòu)模型的映射圖解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)
如果網(wǎng)格坐標(biāo)點滿足下列兩個條件:1)d<1/2t0;2)由點x0x1x2,x0x2x3,x0x3x1組成的三個小三角形面積S1,S2,S3之和等于x1x2x3組成的三角形面積S,即投影點x0應(yīng)該在x1x2x3組成的三角形內(nèi)部,則將其標(biāo)記為多孔區(qū)域的網(wǎng)格。多孔區(qū)域由所有滿足判定條件的網(wǎng)格組成。
采用Bi等[18]所得的試驗結(jié)果進行模型驗證,試驗在波浪水槽中進行,測試網(wǎng)為圓形網(wǎng),在周向有40目,在垂直方向有8目,網(wǎng)衣直徑為dw=1.2 mm,目腳長度為λ=20 mm,圓形框架底部直徑為D=0.254 m,網(wǎng)衣高度為H=0.15 mm,框架為m=8 g的配重,在本算例中,8 g的配重除去浮力的部分將會平均分給40個底部的網(wǎng)衣節(jié)點,每個節(jié)點的受力為0.001 73 N。物理模型和測量點如圖6所示,圓形網(wǎng)衣需要考慮容積損失,通過Zhao等[22]的方法來計算出容積損失率。由于多孔介質(zhì)的厚度影響不大[15],選擇多孔介質(zhì)為20 mm的厚度,多孔介質(zhì)系數(shù)Cn=9.16 m-1,Ct=5.67 m-1。
圖6 圓形網(wǎng)衣的數(shù)值模型和測量點的設(shè)置Fig. 6 Physical model of a circular net and general setting of the measurement points
數(shù)值模型的網(wǎng)格劃分如圖7所示。x為來流方向,y為水平面上與x垂直的方向,z為水平面上法向方向。水槽左端定義為速度入口邊界;右端定義為出流邊界;水槽側(cè)壁和水面定義為滑移邊界;水槽底面定義為固壁邊界。網(wǎng)格類型為非結(jié)構(gòu)化六面體網(wǎng)格,在多孔介質(zhì)區(qū)域向外加密兩層網(wǎng)格。采用有限體積法(FVM)離散控制方程,時間項目采用一階歐拉隱式格式,對流項采用二階迎風(fēng)差分格式,壓力—速度耦合選擇瞬態(tài)的PISO算法,計算時間步長取0.005 s,求解時長為40 s。
圖7 計算網(wǎng)格示例Fig. 7 Example of computational grids
為了保證在OpenFOAM計算準(zhǔn)確,需要進行網(wǎng)格收斂性研究,在流速U=0.242 m/s的情況下,使用三種網(wǎng)格來計算網(wǎng)衣的阻力。圖8是在運用1.2.2節(jié)所述的多孔介質(zhì)生成方法,在OpenFOAM得到的多孔區(qū)域,從多孔區(qū)域看出,其中有部分網(wǎng)格缺失,原因是判斷條件算法的局限性和靜態(tài)網(wǎng)格的使用,多孔區(qū)域只能近似變形網(wǎng)的幾何形狀。表1是在三種不同網(wǎng)格下的計算結(jié)果,計算結(jié)果表明,在不同的全局網(wǎng)格數(shù)和多孔介質(zhì)區(qū)域網(wǎng)格數(shù)下,計算流體力學(xué)(CFD)模擬的網(wǎng)衣阻力隨網(wǎng)格數(shù)變化很小,與試驗值基本吻合,網(wǎng)格收斂性良好,圖9為在三種網(wǎng)格下的切面y=-0.05 m的速度切面云圖,來流通過前網(wǎng)和后網(wǎng)都有明顯的速度降低,兩側(cè)的繞流增加,有明顯的速度過渡區(qū)域。
圖8 OpenFOAM與試驗結(jié)果之間凈變形比較Fig. 8 Comparison of the net deformation between OpenFOAM and the experimental result
表1 OpenFOAM中網(wǎng)格收斂性研究Tab. 1 Mesh convergence study in OpenFOAM
圖9 三種網(wǎng)格下切面y=-0.05 m的速度云圖Fig. 9 Velocity distribution of three kinds of meshes under the slice y=-0.05 m
3.2.1 在code_aster中生成圓形網(wǎng)衣的變形
在圓形網(wǎng)衣的數(shù)值模擬中,三種不同的來流速度(U=0.122、0.178、0.242 m/s)通過底部配重8 g的單片圓形網(wǎng)衣,圖10給出了code_aster計算中網(wǎng)衣在初始時刻和穩(wěn)定狀態(tài)下的位置和形狀,與試驗結(jié)果吻合較好,隨著流速從0.122 m/s增加至0.242 m/s,圓形網(wǎng)衣的變形加劇,網(wǎng)衣的前端變形未能很好地與試驗匹配,主要是因為在試驗中,配重片是以圓環(huán)的形式存在的,而在數(shù)值計算中,是將圓環(huán)的配重平均分給了網(wǎng)衣節(jié)點,這導(dǎo)致了數(shù)值模型和試驗?zāi)P偷牟町?。?為網(wǎng)衣的容積損失率,隨著流速從0.122 m/s增加至0.242 m/s,容積損失率從20.9%增加至48.5%,這與圖10的結(jié)果是符合的。
表2 容積損失率Tab. 2 Volume loss rate
圖10 數(shù)值模擬和模型試驗在不同流速下圓形網(wǎng)衣的變形對比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的水流通過圓形網(wǎng)衣的數(shù)值模擬結(jié)果。從圖11(a)來看,網(wǎng)衣內(nèi)部有一段區(qū)域流速衰減,在網(wǎng)的后方有相當(dāng)長的流速衰減區(qū)域,流速衰減最小達到了0.12 m/s,出現(xiàn)在網(wǎng)兩側(cè)的正后方,繞流速度最大值達到了0.26 m/s,出現(xiàn)區(qū)域為網(wǎng)的兩側(cè)開外區(qū)域。 從圖11(b)來看,網(wǎng)的內(nèi)部和后方流速衰減明顯,后網(wǎng)的后方是一片很大的流速衰減區(qū)域,網(wǎng)的底部繞流明顯,速度分層明顯。從圖11整體來看,網(wǎng)后方的流速降低主要是由于網(wǎng)的阻塞作用,隨著x方向距離柔性網(wǎng)的距離增加,兩側(cè)和底部均有繞流,網(wǎng)正后方速度有明顯衰減。由于使用了Realizeable k-ε湍流模型,在速度的過渡層區(qū)域,急劇的速度梯度被減小,不同速度層之間的動量交換,速度分布比較平滑。
圖11 速度0.242 m/s數(shù)值模擬計算網(wǎng)周圍的流速分布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為在兩個求解器中計算的阻力與試驗數(shù)據(jù)互相對比,從數(shù)值來看吻合較好。在圓形網(wǎng)衣的計算中,結(jié)構(gòu)求解器code_aster的后網(wǎng)來流速度是根據(jù)前網(wǎng)的速度衰減[14]估算得到的,將流速衰減因子r看成常數(shù)是不合理的,這與Bi等[23]的測量結(jié)果不符,因此在code_aster中計算的受力與試驗數(shù)據(jù)相比存在偏差且最大誤差為15.3%,出現(xiàn)在來流速度U=0.178 m/s。OpenFOAM將網(wǎng)衣節(jié)點的坐標(biāo)位置導(dǎo)入計算域并采用多孔介質(zhì)模型計算出流場分布和受力情況,大小與試驗值較為符合。
圖12 不同流速下數(shù)值模擬和試驗數(shù)據(jù)之間的阻力比較Fig. 12 Comparison of the drag force between numerical simulations and experimental data at different current velocities
3.3.1 在code_aster中生成圓形網(wǎng)衣的變形
對具有不同配重m(m=8、45、367 g)的網(wǎng)衣進行數(shù)值模擬,來流速度均為U=0.242 m/s,圖13表示了code_aster的計算結(jié)果和試驗結(jié)果對比,給出網(wǎng)衣的初始狀態(tài)和變形狀態(tài),數(shù)值模擬結(jié)果和試驗結(jié)果吻合較好,隨著配重從m=8 g增加至m=367 g,圓形網(wǎng)衣的變形減小,在w=367 g的變形十分小。表3為網(wǎng)衣的容積損失率,隨著配重從m=8 g增加至m=367 g,容積損失率從48.5%減小至2.0%,這與圖13所得的結(jié)果相符合。
表3 容積損失率Tab. 3 Volume loss rate
圖13 數(shù)值模擬和模型試驗在同流速不同配重下圓形網(wǎng)衣的變形對比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中的模擬結(jié)果,結(jié)果表明,網(wǎng)的內(nèi)部和后方都有明顯的流速衰減,底部有明顯的繞流,過渡區(qū)域明顯。相較于配重大的網(wǎng)衣,配重m=8 g的網(wǎng)衣后方流速衰減更大,這表明,網(wǎng)衣的屏蔽作用隨著配重的增加而降低,尤其在網(wǎng)衣的正后方。
圖14 不同配重下數(shù)值計算切面y=0的流速分布Fig. 14 Velocity distribution of slice y=0 simulated by numerical simulation in different weights
圖15為沿著圖6監(jiān)測點的線中平均流速的分布,可以看出,流速與試驗數(shù)據(jù)大小相近,變化趨勢一致,最大的相對誤差為 5.6%,模擬結(jié)果較好。在配重m=8 g的圓形網(wǎng)衣模擬中,在x=1.25 m左右流速趨向穩(wěn)定,而在配重m=45 g和m=367 g的情況下,在x=0.75 m左右流速趨向穩(wěn)定,配重m=8 g的穩(wěn)定流速要比配重m=45 g和m=367 g更小。
圖15 數(shù)值模擬和試驗中測量點的平均流速的比較Fig. 15 Comparison of the average flow-velocity between numerical simulations and experimental measurement points
3.3.3 OpenFOAM和code_aster受力對比
由于在模型試驗中并未提供試驗數(shù)據(jù),本文只進行兩個求解器中網(wǎng)衣阻力的對比,圖16為兩個求解器的阻力對比,在OpenFOAM中的計算結(jié)果要比code_aster中大,這與圖12的結(jié)果較為一致。隨著配重的增加,阻力也持續(xù)加大。
圖16 不同流速下數(shù)值模擬和試驗數(shù)據(jù)之間的阻力比較Fig. 16 Comparison of the drag force between numerical simulations and experimental data at different current velocities
利用結(jié)構(gòu)求解器code_aster和流體求解器OpenFOAM提出了一種更簡便的柔性網(wǎng)衣和流場耦合計算模型,通過圓形網(wǎng)衣進行了驗證。具體做法是,首先采用Screen模型計算網(wǎng)衣的水動力,通過桁架模型code_aster計算網(wǎng)衣的結(jié)構(gòu)變形,模擬結(jié)果與試驗相符,然后導(dǎo)出網(wǎng)衣節(jié)點坐標(biāo)至OpenFOAM中,根據(jù)判斷條件在流場的網(wǎng)格區(qū)域內(nèi)生成多孔區(qū)域,通過加多孔介質(zhì)模型源項的方法,利用OpenFOAM在不同來流和不同配重情況下進行了流場分析,網(wǎng)衣后方有明顯的流速衰減,兩側(cè)和底部繞流明顯,流速衰減與試驗數(shù)據(jù)趨勢一致,大小相近,說明多孔介質(zhì)模型的有效性,最后比較了試驗所得網(wǎng)衣阻力與code_aster、OpenFOAM計算的網(wǎng)衣阻力,誤差均在合理范圍之間。上述結(jié)果均顯示,與Bi等[18]的模擬結(jié)果和試驗數(shù)據(jù)吻合,說明了桁架模型和多孔介質(zhì)模型耦合方法的準(zhǔn)確性。
在柔性網(wǎng)衣與流場的耦合計算中,時間開銷主要由流體求解器承擔(dān),結(jié)構(gòu)求解器占比較小。由于本模型采用的是單向耦合方式,因此不必考慮兩個求解器的數(shù)據(jù)交換和時間同步。其次,Screen水動力模型是由經(jīng)驗公式得到,在流速不高的條件下求解網(wǎng)衣變形是比較可靠的。因此本文提出的基于桁架模型和多孔介質(zhì)模型建立的柔性網(wǎng)衣和流場的單向耦合方法對于實際問題的研究是經(jīng)濟、有效的。