波流數值水池模擬研究

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

(江蘇科技大學船舶與海洋工程學院，江蘇鎮江212003)

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

1 數學模型

1.1 控制方程

1)連續性方程:質量守恒方程

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

式中:ρ為流體密度;v為流場速度矢量.

2)N-S方程:動量守恒方程

假設流體粘性系數為常數，則不可壓縮流體的

運動方程，即N-S方程為:

式中:μ為動力粘性系數;Fb為體力;p為壓強;S為應變率張量.式(3)左端為單位體積流體的慣性力;右端第1項為單位體積的質量力;第2項為作用于單位體積流體的壓強梯度力;第3項為粘性變形應力;第4項為粘性體膨脹應力.

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

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

數值水池底部法向速度為零，即:

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

每個控制體中，兩相的體積分數之和為1，即:

兩相之間的流體界面通過求解體積分數連續性方程來追蹤，第q相流體的體積分數滿足連續性方程:

1.3 線性小振幅波理論

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

波面方程:

速度勢:

x方向速度:

z方向速度:

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

1.4 消波理論

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

消波區內，動量方程(二維)為:

式中:C(x)為消波系數;υ為運動粘性系數.

式中:c為系數，不同的數值水池其值不同，根據經驗及試算，算例中c取10;λ為波長;x0，xL分別為消波區左、右邊界的x坐標值;ρ為流體密度.

2 波流數值水池模擬

2.1 建模與網格劃分

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

圖1 數值水池Fig.1 Schematic diagram of the numerical tank

圖2 網格劃分Fig.2 Schematic diagram of the grid

2.2 邊界條件及計算參數設置

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

采用FLUENT中二維非定常分離隱式求解器求解，壓力方程選用加權體積力格式(body force weighted)，壓力速度耦合方式采用PISO算法，體積分數方法為幾何重構(geo-reconstruct)，壓力參考值為一個標準大氣壓，重力加速度為-9.81m/s2，時間步長為0.005s.

2.3 波浪數值水池驗算

為了驗證波流數值水池的正確性，首先對具有如下波浪要素的余弦波進行了數值模擬.

表1 余弦波Table 1 Cosine wave

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

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

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

圖4 自由液面坐標時程曲線Fig.4 Time curve of the surface

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

圖5 不同時刻消波區波形Fig.5 Waves′shape in the damping domain at different time

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

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

2.4 波浪與均勻流數值模擬及分析

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

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

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

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

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

對于無限水深，設L為波流場波長;A為波流場波幅;L0為單獨波浪場波長;A0為單獨波浪場波幅，根據文獻［10］，有

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

表2 波長計算結果Table 2 Calculation results of wavelength

表3 波幅計算結果Table 3 Calculation results of amplitude

圖8 計算結果與理論值對比曲線Fig.8 Contrast curve of the calculations and the theoretical values

2.5 波浪與剪切流聯合作用數值模擬

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

圖9 t=50s時波形Fig.9 Waves′shape at t=50s

表4 波長、波幅計算結果Table 4 Calculation results of wavelength and amplitude

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

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

3 結論

以商業軟件ANSYS中的FLUENT為平臺，通過加載UDF程序進行二次開發實現速度入口造波，采用VOF方法追蹤氣液兩項自由面，在水池后方動量方程中添加源項形成阻尼消波，成功建立了二維波浪、波流數值水池.

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

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

3)模擬了波浪與不同剪切率剪切流聯合作用的波浪——剪切流數值水池，探究了剪切率對波流數值水池的影響，通過擬合，得到了波長及波幅與剪切率的函數關系，以便工程應用.

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］ 趙艷，朱仁慶，劉珍，等.三維數值波浪水池的構建和粘性影響研究［J］.艦船科學技術，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］ 楊錦凌，孫大鵬，吳浩，等.斜坡堤波浪爬高和越浪數值模擬［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］ 劉莎莎，顧煜炯，惠萬馨，等.基于邊界造波法的波浪數值模擬［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］ 李玉成.波浪與水流共同作用下的流速場［J］.海洋工程，1983，(4):12 -23.

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