波流數(shù)值水池模擬研究

吳梓鑫，朱仁慶，繆志剛，劉 珍，顧思琪

(江蘇科技大學(xué)船舶與海洋工程學(xué)院，江蘇鎮(zhèn)江212003)

在復(fù)雜的海洋環(huán)境中，波浪與水流通常是共同存在的.潮流、洋流、河口處徑流以及與波浪破碎有關(guān)的沿岸流、裂流及其所形成的近岸環(huán)流等水流會(huì)與波浪產(chǎn)生相互作用，使波幅、波高、周期等波浪要素發(fā)生變化，并使波浪發(fā)生折射甚至破碎.因此，研究波流的相互作用是開發(fā)海洋需要解決的關(guān)鍵問題，受到了國內(nèi)外科研人員的重視.1911年，Rayleigh發(fā)現(xiàn)在潮汐中存在著波流的相互作用.文獻(xiàn)［1-3］中通過實(shí)驗(yàn)研究了水流對(duì)波浪要素的影響;文獻(xiàn)［4-5］中在數(shù)值理論方面做出了一定的貢獻(xiàn).但是，波流相互作用是極其復(fù)雜的非線性問題，其流動(dòng)機(jī)理和作用特性仍需深入研究.隨著計(jì)算機(jī)技術(shù)的成熟，數(shù)值水池已具備了信息量豐富、計(jì)算結(jié)果精確的特點(diǎn)，加上其成本低、改變條件方便、重復(fù)性好等優(yōu)勢(shì)，數(shù)值水池已經(jīng)成為國際上研究海洋波浪問題的熱門方法.目前，眾多學(xué)者已經(jīng)通過推板造波［6］、源項(xiàng)造波［7］、速度入口［8］等數(shù)值方法模擬出理想波浪.然而針對(duì)波流相互影響的研究還很少，此外，波浪與剪切流作用的數(shù)值模擬更是被很多人忽視.文中基于軟件平臺(tái)FLUENT二次開發(fā)構(gòu)造二維數(shù)值波流水池，采用速度入口造波技術(shù)，不可壓縮粘性流體Navier-Stokes方程為控制方程，通過流體體積法(volume of fluid，VOF)追蹤自由液面，采用阻尼消波法在水池后方進(jìn)行消波處理，成功模擬了波浪與均勻流及與剪切流共同作用的數(shù)值水池，并分析了水流對(duì)波浪的影響.

1 數(shù)學(xué)模型

1.1 控制方程

1)連續(xù)性方程:質(zhì)量守恒方程

對(duì)于不可壓縮流體，式(1)簡化為:

式中:ρ為流體密度;v為流場(chǎng)速度矢量.

2)N-S方程:動(dòng)量守恒方程

假設(shè)流體粘性系數(shù)為常數(shù)，則不可壓縮流體的

運(yùn)動(dòng)方程，即N-S方程為:

式中:μ為動(dòng)力粘性系數(shù);Fb為體力;p為壓強(qiáng);S為應(yīng)變率張量.式(3)左端為單位體積流體的慣性力;右端第1項(xiàng)為單位體積的質(zhì)量力;第2項(xiàng)為作用于單位體積流體的壓強(qiáng)梯度力;第3項(xiàng)為粘性變形應(yīng)力;第4項(xiàng)為粘性體膨脹應(yīng)力.

對(duì)于不可壓縮流體，N-S方程簡化為:

1.2 邊界條件及自由液面追蹤方法

數(shù)值水池底部法向速度為零，即:

文中采用VOF法追蹤氣液交界面.定義體積分?jǐn)?shù)函數(shù)αq為單元內(nèi)第q相流體所占有全部的體積與單元總體積之間的比.若αq=0，表明該單元內(nèi)全部為空氣;若αq=1，則表明該單元內(nèi)都為液體;若0＜αq＜1，即為交界面.

每個(gè)控制體中，兩相的體積分?jǐn)?shù)之和為1，即:

兩相之間的流體界面通過求解體積分?jǐn)?shù)連續(xù)性方程來追蹤，第q相流體的體積分?jǐn)?shù)滿足連續(xù)性方程:

1.3 線性小振幅波理論

線性波又稱Airy波，對(duì)于無限水深平面線性波，其解為:

波面方程:

速度勢(shì):

x方向速度:

z方向速度:

式中:H為波高;k=2π/λ;ω =2π/T;λ為波長;T為波浪周期;z為垂向坐標(biāo)(向上為正);z=0為靜水面.

1.4 消波理論

物理水池中，波浪傳播到另一端池壁，將會(huì)產(chǎn)生反射，因而需在水池末端設(shè)置消波裝置以消除反射波.對(duì)于數(shù)值水池，同樣需要在水池的后部進(jìn)行消波處理.數(shù)值消波主要有阻尼消波法和多孔介質(zhì)消波法，這里采用阻尼消波法.

消波區(qū)內(nèi)，動(dòng)量方程(二維)為:

式中:C(x)為消波系數(shù);υ為運(yùn)動(dòng)粘性系數(shù).

式中:c為系數(shù)，不同的數(shù)值水池其值不同，根據(jù)經(jīng)驗(yàn)及試算，算例中c取10;λ為波長;x0，xL分別為消波區(qū)左、右邊界的x坐標(biāo)值;ρ為流體密度.

2 波流數(shù)值水池模擬

2.1 建模與網(wǎng)格劃分

采用ANSYS 14.5中GEOMETRY模塊建模.數(shù)值水池尺寸為30 m×5 m，水深3 m，水面以上為空氣，水池后部10 m為消波區(qū).用MESH模塊劃分網(wǎng)格，在波長方向上，單位網(wǎng)格長度為波長的1/100;在水深方向，為了更精確地捕捉到自由液面，在靜水面上下一個(gè)波高內(nèi)加密網(wǎng)格，網(wǎng)格高度為波高的1/60，在遠(yuǎn)離自由液面處網(wǎng)格漸疏.模型及網(wǎng)格劃分情況如圖1，2所示.

圖1 數(shù)值水池Fig.1 Schematic diagram of the numerical tank

圖2 網(wǎng)格劃分Fig.2 Schematic diagram of the grid

2.2 邊界條件及計(jì)算參數(shù)設(shè)置

水池左端AB設(shè)置為速度入口(velocity-inlet)，上方邊界AE為壓力入口(pressure-inlet)，底部BF為固壁(wall)，右端EF為壓力出口(pressure-outlet)，區(qū)域CEFD為消波區(qū).

采用FLUENT中二維非定常分離隱式求解器求解，壓力方程選用加權(quán)體積力格式(body force weighted)，壓力速度耦合方式采用PISO算法，體積分?jǐn)?shù)方法為幾何重構(gòu)(geo-reconstruct)，壓力參考值為一個(gè)標(biāo)準(zhǔn)大氣壓，重力加速度為-9.81m/s2，時(shí)間步長為0.005s.

2.3 波浪數(shù)值水池驗(yàn)算

為了驗(yàn)證波流數(shù)值水池的正確性，首先對(duì)具有如下波浪要素的余弦波進(jìn)行了數(shù)值模擬.

表1 余弦波Table 1 Cosine wave

圖3為波浪生成和傳播過程示意圖，上部分為空氣，下部分為水.可以看出，在t=30s時(shí)，水池已形成穩(wěn)定的波形.

圖3 波浪的生成和傳播過程Fig.3 Generation and propagation process of the wave

在數(shù)值波浪水池中x=1/5λ，λ(即x=5m)處設(shè)置波高儀以實(shí)時(shí)監(jiān)測(cè)波高.圖4為所監(jiān)測(cè)波高儀隨時(shí)間t變化時(shí)程曲線，藍(lán)色為數(shù)值計(jì)算結(jié)果曲線，紅色為理論曲線，可以得到其運(yùn)動(dòng)的周期為1.79s，與理論波形的周期基本吻合，波幅為0.148 m，與理論值誤差為1.33%.

圖4 自由液面坐標(biāo)時(shí)程曲線Fig.4 Time curve of the surface

圖5為消波區(qū)波形隨時(shí)間走勢(shì)圖，從圖中可以清晰看出波形衰減的過程，在水池最右端，波浪已基本完全消失，驗(yàn)證了文中所用阻尼消波的可行性.

圖5 不同時(shí)刻消波區(qū)波形Fig.5 Waves′shape in the damping domain at different time

圖6為數(shù)值水池在靜水狀態(tài)和波浪狀態(tài)下的靜壓力云圖.a)圖顯示的是靜水狀態(tài)的壓力云圖，圖中壓力狀態(tài)呈分層分布，壓力從上至下逐漸變大，且最大靜壓力的值與理論值相當(dāng).b)圖展示的是波浪狀態(tài)下的壓力云圖，可以看出波谷處的靜水壓力值小于波峰處的靜水壓力值.從而進(jìn)一步通過靜水和波浪兩個(gè)狀態(tài)分別驗(yàn)證了數(shù)值模擬波浪的正確性.

圖6 靜水壓力云圖Fig.6 Hydrostatic pressure contours

2.4 波浪與均勻流數(shù)值模擬及分析

目前，為兼顧準(zhǔn)確性和工程需要，波浪和水流共同作用一般采取文獻(xiàn)［9］中提出的波流共同作用的流速場(chǎng)理論，即波流水平速度疊加法，其速度關(guān)系如下:

式中:uwc(x，z，t)為波流聯(lián)合場(chǎng)水平速度;uw(x，z，t)為波浪場(chǎng)水平速度;uc為均勻流場(chǎng)水平速度.

文中對(duì)流速uc=uw/4，uc=uw/8和uc=-uw/8的3種工況波流水池進(jìn)行數(shù)值模擬，分別命名為算例1，2，3，uw=2.79m/s.為保證波流數(shù)值水池與上述單獨(dú)波浪水池結(jié)果的可比性，其幾何模型、邊界條件、物理屬性、計(jì)算參數(shù)設(shè)置以及消波區(qū)域設(shè)置與單獨(dú)波浪水池保持相同.

圖7為算例0～3波長方向x=15 m處波高時(shí)程曲線，通過比較可看出，由于水流的作用，形成穩(wěn)定波形所需時(shí)間發(fā)生了變化，當(dāng)波流同向時(shí)，流體對(duì)波浪具有推動(dòng)作用，其傳播速度加快，形成穩(wěn)定波形較單獨(dú)波浪所需時(shí)間要短;反之，波流相向時(shí)，流體對(duì)波浪傳播具有阻礙作用，波流場(chǎng)形成穩(wěn)定波形則所需時(shí)間變長.

圖7 x=15m處波幅時(shí)程曲線Fig.7 Time curve of the wave amplitude at x=15m

對(duì)于無限水深，設(shè)L為波流場(chǎng)波長;A為波流場(chǎng)波幅;L0為單獨(dú)波浪場(chǎng)波長;A0為單獨(dú)波浪場(chǎng)波幅，根據(jù)文獻(xiàn)［10］，有

表2，3為波長和波幅計(jì)算結(jié)果分析表，數(shù)據(jù)顯示數(shù)值模擬很好地驗(yàn)證了波流場(chǎng)中水流對(duì)波浪波長及波幅的影響.當(dāng)波流同向時(shí)，由于波速增大，波形被拉長，波幅減小;反之，波流相向時(shí)，波速減小，波形變窄，波幅增加.圖8為波長和波幅計(jì)算結(jié)果與理論值對(duì)比曲線，數(shù)值結(jié)果與理論結(jié)果基本吻合.

表2 波長計(jì)算結(jié)果Table 2 Calculation results of wavelength

表3 波幅計(jì)算結(jié)果Table 3 Calculation results of amplitude

圖8 計(jì)算結(jié)果與理論值對(duì)比曲線Fig.8 Contrast curve of the calculations and the theoretical values

2.5 波浪與剪切流聯(lián)合作用數(shù)值模擬

在實(shí)際工程應(yīng)用中，基本不存在完全的均勻流，海底河床等因素必然形成剪切流等非均勻流.文中對(duì)r=0.03，0.06，0.09和0.12的4種剪切率進(jìn)行波浪與剪切流的數(shù)值模擬，分別記為算例4～7，水池幾何模型、邊界條件、物理屬性、計(jì)算參數(shù)設(shè)置以及消波區(qū)域設(shè)置與單獨(dú)波浪水池保持相同，在水深1.5m處，保持所加流速為波速的1/8，即uc=1/8uw.圖9為t=50s時(shí)4種工況波形圖，波長及波幅變化情況如表4所示.

圖9 t=50s時(shí)波形Fig.9 Waves′shape at t=50s

表4 波長、波幅計(jì)算結(jié)果Table 4 Calculation results of wavelength and amplitude

與單獨(dú)波浪場(chǎng)相比，疊加剪切流后，由于流的作用，波浪波長增加了，波幅則減小，變化趨勢(shì)與疊加均勻流一致.波長及波幅隨剪切率變化曲線如圖10所示，根據(jù)數(shù)值模擬結(jié)果，文中對(duì)波長和波幅進(jìn)行函數(shù)擬合，結(jié)果如下:

圖10 λ/λ0及A/A0隨剪切率變化曲線Fig.10 Curve of λ/λ0and A/A0changing with the rate

3 結(jié)論

以商業(yè)軟件ANSYS中的FLUENT為平臺(tái)，通過加載UDF程序進(jìn)行二次開發(fā)實(shí)現(xiàn)速度入口造波，采用VOF方法追蹤氣液兩項(xiàng)自由面，在水池后方動(dòng)量方程中添加源項(xiàng)形成阻尼消波，成功建立了二維波浪、波流數(shù)值水池.

1)將設(shè)置造波邊界法應(yīng)用于FLUENT軟件，成功模擬了效果良好的單獨(dú)波浪以及波浪—均勻流和波浪—剪切流數(shù)值水池，驗(yàn)證了該造波方式和所用阻尼消波技術(shù)的可行性;

2)通過在入口處定義波流疊加的邊界條件，模擬了波浪與不同流速聯(lián)合作用的波浪—均勻流數(shù)值水池.表明波流同向時(shí)，水流對(duì)波浪傳播具有推動(dòng)作用，使波形變長，波幅減小;波流相向時(shí)，水流阻礙波浪傳播，使波形變短，波幅增加，與理論結(jié)果基本吻合;

3)模擬了波浪與不同剪切率剪切流聯(lián)合作用的波浪——剪切流數(shù)值水池，探究了剪切率對(duì)波流數(shù)值水池的影響，通過擬合，得到了波長及波幅與剪切率的函數(shù)關(guān)系，以便工程應(yīng)用.

References)

［1］ Thomas G P.Wave-current interactions:an experimental and numerical study.part 1.linear waves［J］.Journal of Fluid Mechanics，1981，110:457 -474.

［2］ Kemp P H，Simons R R.The interaction between waves and a turbulent current:waves propagating against the current［J］.Journal of Fluid Mechanics，1983，130:73 -89.

［3］ Hughes B A，Stewart R W.Interaction between gravity waves and as hear flow［J］.Journal of Fluid Mechanics，1961，10:385 -402.

［4］ Longuet M S，Stewart R W.Changes in the form of short gravity waves on long waves and tidal currents［J］.Journal of Fluid Mechanics，1960，8:565 -583.

［5］ Grant W D，Madsen O S.Combined wave and current interaction with a rough bottom［J］.Journal of Geophysical Research，1979，84:1797 -1808.

［6］ 趙艷，朱仁慶，劉珍，等.三維數(shù)值波浪水池的構(gòu)建和粘性影響研究［J］.艦船科學(xué)技術(shù)，2014，36(5):42-48.Zhao Yan，Zhu Renqing，Liu Zhen，et al.Simulation of 3D numerical wave tank and viscosity research［J］.Ship Science and Technology，2014，36(5):42 - 48.(in Chinese)

［7］ 楊錦凌，孫大鵬，吳浩，等.斜坡堤波浪爬高和越浪數(shù)值模擬［J］.海洋工程，2013，31(2):45-52.Yang Jinlin，Sun Dapen，Wu Hao，et al.Numerical simulation of wave run-up and overtopping on sloping seawall［J］.The Ocean Engineering，2013，31(2):45 -52.(in Chinese)

［8］ 劉莎莎，顧煜炯，惠萬馨，等.基于邊界造波法的波浪數(shù)值模擬［J］.可再生能源，2013，31(2):100-103.Liu Shasha，Gu Yujiong，Hui Wanxin，et al.Wave numerical simulation based on wave-generation method of defining inlet boundary conditions.［J］.Renewable Energy Resources，2013，31(2):100 -103.(in Chinese)

［9］ 李玉成.波浪與水流共同作用下的流速場(chǎng)［J］.海洋工程，1983，(4):12 -23.

［10］ 鄒志利.水波理論及其應(yīng)用［M］.北京:科學(xué)出版社，2005:472-482.