內(nèi)孤立波波致流場數(shù)值模擬研究*

王 偉 郭海燕① 王 飛 苗得勝 馬 東

(1. 中國海洋大學(xué)工程學(xué)院 青島 266100; 2. 山東科技大學(xué)土木工程與建筑學(xué)院 青島 266590)

隨著南海資源開發(fā)的不斷深入, 其復(fù)雜的海洋環(huán)境因素越來越受到專家學(xué)者的重視。南海內(nèi)孤立波發(fā)生頻繁, 對水下結(jié)構(gòu)的正常工作產(chǎn)生嚴重影響, 給南海油氣資源開發(fā)帶來了極大挑戰(zhàn), 是南海海洋工程結(jié)構(gòu)建設(shè)中不可忽略的環(huán)境荷載之一。

內(nèi)孤立波發(fā)生在密度穩(wěn)定分層的海洋內(nèi)部, 其振幅遠大于表面波(杜濤等, 2001), 主要通過非線性效應(yīng)與色散效應(yīng)達到一定程度的平衡實現(xiàn)穩(wěn)定傳播(Osborne et al, 1980), 目前應(yīng)用比較廣泛的內(nèi)孤立波理論模型主要有 KdV、mKdV理論等(Choi et al,1999;Helfrich et al, 2006)。近幾年, 內(nèi)孤立波的研究方法有現(xiàn)場觀測、理論研究、物理實驗、數(shù)值模擬等。Liu等(2014)利用在中國南海獲得的內(nèi)孤立波SAR影像, 探索了海底地形以及水深變化對內(nèi)孤立波傳播速度的影響。2005年5月Xu等(2010)在南海西北陸架監(jiān)測到高非線性的內(nèi)孤立波, 通過將內(nèi)孤立波數(shù)據(jù)與KdV、mKdV理論進行對比, 發(fā)現(xiàn)對于高非線性內(nèi)孤立波, eKdV理論吻合效果更好?？伦悦鞯?2009)在南海文昌內(nèi)波實驗中, 通過單點錨系測量觀測到典型內(nèi)孤立波, 并分析內(nèi)孤立波波形、溫度結(jié)構(gòu)及波致斜壓流。Michallet等(1998)通過將波形、相速度、振幅的理論模型計算結(jié)果與實驗結(jié)果對比, 提出在振幅較小時 KdV理論與實驗結(jié)果吻合較好, 而振幅較大時KdV-mKdV理論與實驗結(jié)果較接近; Buick等(2003)通過實驗室造波, 采用PIV技術(shù)監(jiān)測內(nèi)波流場,驗證了理論計算結(jié)果與實驗結(jié)果。近幾年, 利用數(shù)值水槽進行模擬的計算流體力學(xué)(CFD)方法迅速發(fā)展起來。Zhang等(2012)使用CFD方法建立三維數(shù)值波浪水槽, 所得結(jié)果與 KdV理論計算結(jié)果較吻合, 證明該方法可有效模擬非線性內(nèi)孤立波。陳鈺等(2009)利用Fluent軟件, 通過設(shè)置造波邊界法, 成功建立了可有效模擬弱非線性內(nèi)孤立波的分層流數(shù)值水槽; 高原雪等(2012)基于 MCC理論, 利用速度入口數(shù)值造波方法, 建立內(nèi)孤立波數(shù)值水槽, 并實現(xiàn)了振幅可控的內(nèi)孤立波數(shù)值模擬。目前國內(nèi)外主要將 CFD數(shù)值模擬結(jié)果與理論結(jié)果進行對比分析, 而結(jié)合物理實驗結(jié)果對比分析的研究較少。

本文采用CFD數(shù)值模擬方法, 基于Fluent軟件,通過“平板拍擊”方法, 建立具有多種設(shè)計振幅的內(nèi)孤立波二維數(shù)值水槽。通過VOF方法(董志等, 2009)追蹤兩層流體界面變化, 利用軟件自帶的監(jiān)視器功能觀測不同高度處的波致水平流速, 并與物理實驗結(jié)果進行對比, 驗證數(shù)值模擬的有效性。利用數(shù)值造波結(jié)果, 研究內(nèi)孤立波波致流場特性, 分析不同水深處的水平流速, 以及波谷經(jīng)過斷面時波致水平流速垂向分布。

1 建立數(shù)值水槽

采用兩層流體模型, 使用Navier-Stokes方程作為流場的控制方程。建立內(nèi)孤立波二維數(shù)值水槽。數(shù)值水槽尺寸如下: 長 1500cm, 高 58cm, 其中上層流體高度9.5cm, 下層流體高度48.5cm。利用直角坐標(biāo)系對水槽進行定位, x軸與兩層流體的交界面重合, 正方向沿水槽長度方向向右; y軸與水槽左邊界重合, 且正方向沿水槽高度垂直向上。數(shù)值水槽模型如圖1所示, 其中造波板長度為80cm, 擋板長度為15cm。

圖1 數(shù)值水槽示意圖Fig.1 The sketch of numerical wave tank

其邊界條件設(shè)置為: 造波板設(shè)置為移動邊界, 擋板、水槽的下邊界及左、右邊界設(shè)為壁面邊界(wall)。根據(jù)剛蓋近似原理(方欣華等, 2005), 內(nèi)孤立波在兩層流體的自由表面上引起的波動較小, 故將水槽的上表面設(shè)為壁面邊界(wall)。

“平板拍擊”造波方法原理為: 選取合適的內(nèi)波理論, 通過 UDF用戶自定義函數(shù)來控制造波板的運動速度, 使造波板作上下運動, 并利用動網(wǎng)格技術(shù),在兩層流體界面產(chǎn)生內(nèi)孤立波。使用 Navier-Stokes方程作為流場的控制方程, 在兩層流體模式下, 描述內(nèi)孤立波運動的方程主要有KdV、mKdV等。設(shè)上下層流體水深、密度分別為:1h、1ρ和2h、2ρ, (,)x tη為波面方程。KdV理論下方程解的形式(Michallet et al, 1998)如下:

式中:

mKdV理論下方程解的形式(Michallet et al,1998)為:

式中:

a為內(nèi)孤立波振幅, 負號表示內(nèi)孤立波波面是向下凹的。L表示內(nèi)孤立波的特征波長。根據(jù)KdV、mKdV理論對不同振幅水深比(a/H)的使用條件, 可分別選用 KdV或 mKdV理論編寫控制造波板運動速度的UDF程序。其中造波板的運動速度 v(t)可定義為(徐鑫哲, 2012):D為造波板的長度, c為內(nèi)孤立波的相速度。其數(shù)學(xué)原理為: 造波板上下運動過程中, 根據(jù)體積守恒, 單位時間內(nèi)造波板排開的流體體積與造波板右端兩層流體交界面沿 x軸前進的體積相等, 即:

其中, Δx = cΔ t 。因此:

數(shù)值水槽的網(wǎng)格劃分: 采用結(jié)構(gòu)化網(wǎng)格對水槽進行網(wǎng)格離散。為了方便、準(zhǔn)確的捕捉內(nèi)孤立波波面, 在內(nèi)孤立波形成的高度范圍內(nèi), 即兩層流體的交界面附近, 對網(wǎng)格劃分進行加密處理。其中在垂直方向y= –8.5cm至y=0.5cm范圍內(nèi), 劃分的網(wǎng)格大小為 0.5cm, 而為了提高計算效率, 在加密區(qū)兩側(cè)的高度范圍內(nèi), 分別從加密區(qū)開始, 向兩側(cè)方向采用漸變網(wǎng)格劃分方法, 網(wǎng)格尺寸的漸變比例為1.02。由于水平方向即x向的網(wǎng)格劃分對內(nèi)孤立波波面的形成影響較小, 故沿x向的網(wǎng)格劃分可以取為2cm一個網(wǎng)格。數(shù)值水槽的網(wǎng)格劃分情況如下圖所示:

圖2 數(shù)值水槽網(wǎng)格劃分示意圖Fig.2 The gridding of the numerical flume

求解設(shè)置: 將建立好的數(shù)值水槽模型導(dǎo)入Fluent中進行求解設(shè)置, 設(shè)置上層流體密度1ρ=1035kg/m3,下層流體密度2ρ=1054kg/m3。選取Fluent自帶的kε-湍流模型, 采用VOF方法追蹤兩層流體的交界面, 對控制方程的離散使用有限體積法, 壓力速度耦合方法選用 PISO算法, 對流相及輸運方程的離散采用一階迎風(fēng)格式。

由于數(shù)值水槽的右邊界為固壁邊界, 內(nèi)孤立波傳播到該位置處時易產(chǎn)生反射波, 為避免反射波對流場產(chǎn)生干擾, 采用阻尼消波(孫大鵬等, 2000; 韓朋等, 2009)的方法在水槽末端進行消波處理: 取水槽末端 1—2倍波長處作為消波段, 通過在消波段動量方程中添加阻尼源項使波浪傳播到該消波段時能量逐漸被吸收。

按照振幅大小分為四個工況, 將四個工況下的內(nèi)孤立波振幅分別作為數(shù)值水槽的設(shè)置波高, 進行數(shù)值造波模擬計算:

表1 數(shù)值模擬工況設(shè)置Tab.1 Setting of amplitude for the numerical simulation

根據(jù) KdV、mKdV理論對不同振幅水深比(a/H)的使用條件, 由于工況1、2中a/H＜0.1, 故選用KdV理論定義造波板的速度; 而工況三、四中a/H＞0.1, 故選用mKdV理論定義造波板的速度(黃文昊等, 2013)。

迭代計算: 各項參數(shù)設(shè)置完成后進行流場的迭代計算, 設(shè)置迭代時間步長為 0.01s, 每個時間步的最大迭代計算次數(shù)為20次。

2 數(shù)值模擬結(jié)果分析

2.1 內(nèi)波物理實驗方法

開展內(nèi)波物理實驗檢驗數(shù)值模擬結(jié)果, 該實驗在中國海洋大學(xué)物理海洋實驗室分層流水槽中進行,造波方法采用重力塌陷法。水槽尺寸為: 15m×0.35m×0.7m(長、寬和高), 額定水深為0.6m。密度分層采用Oster雙缸法制取, 上層流體高h1=9.5cm, 密度ρ1=1035kg/m3; 下層流體高度h2=48.5cm, 密度ρ2=1054kg/m3。圖3為物理造波示意圖。

使用PIV(Wangetal, 2015)技術(shù), 將體積很小、密度與流體相當(dāng)、反光性能良好的鍍銀空心微珠粒子按一定濃度布置于液體中, 粒子能隨流體運動, 代表流體質(zhì)點的運動。待內(nèi)波穩(wěn)定傳播后, 在水槽中部采用CCD系統(tǒng)對粒子的位置進行拍攝, 采樣頻率為50Hz,每張圖像分辨率為1920×1080像素。運用圖像處理軟件進行分析, 可以得到斷面處粒子的運動狀態(tài), 即該處的流場狀態(tài)。

圖3 造波示意圖Fig.3 Experimental wave generation

2.2 數(shù)值造波結(jié)果驗證

首先分析內(nèi)孤立波數(shù)值模擬得到的波形, 并與實驗結(jié)果進行對比, 以驗證本次數(shù)值造波的合理性:

圖5為數(shù)值造波得到的波形數(shù)據(jù), 并分別與內(nèi)孤立波物理實驗結(jié)果進行對比。由圖中可以看出, 工況一、二的數(shù)值造波結(jié)果與實驗結(jié)果吻合較好, 只是在波形尾端數(shù)值造波結(jié)果存在較小的波動, 這是內(nèi)孤立波在生成過程中會產(chǎn)生尾波所致, 但在內(nèi)波傳播過程中會逐漸衰減, 并不會對內(nèi)波產(chǎn)生影響。工況三的數(shù)值造波結(jié)果在波形方面與實驗一致, 且尾波較小。工況四的數(shù)值造波結(jié)果在波形的后半段與實驗一致, 在波形前半段數(shù)值造波結(jié)果優(yōu)于實驗結(jié)果, 因此工況四的造波結(jié)果也是合理的。

綜上, 本文建立的數(shù)值水槽能夠?qū)崿F(xiàn)設(shè)計振幅的內(nèi)孤立波數(shù)值造波, 且數(shù)值結(jié)果與實驗結(jié)果較吻合, 能有效實現(xiàn)對內(nèi)孤立波的數(shù)值模擬。

圖4 數(shù)值造波波形與實驗對比Fig.4 Comparison in shape of the waves generated between numerical simulation (red lines) and experiment (black lines) in different amplitudes

3 討論

3.1 上下層流體中波致水平流速分布

在數(shù)值造波過程中, 針對四個工況, 分別設(shè)置監(jiān)視器觀測上下層流體中間位置處波致水平流速的變化。觀測位置高度: (1)y=5cm, 即上層流體中間位置處; (2)y= –24cm, 即下層流體中間位置處。觀測結(jié)果如下圖所示, 并與實驗結(jié)果進行對比:

圖5 內(nèi)孤立波波致水平流速變化Fig.5 Comparison in horizontal velocity magnitude induced by generated internal solitary waves between numerical simulation(red lines) and experiment (black lines) at different depth

圖6為內(nèi)孤立波波致水平流速變化圖, 并與實驗結(jié)果進行對比, 經(jīng)過對比發(fā)現(xiàn), 工況一、二中, 波致水平流速的數(shù)值模擬結(jié)果的變化與實驗一致; 工況三、四中, 在下層流體中, 數(shù)值模擬結(jié)果的變化與實驗一致, 而在上層流體中差異較大。原因為: 在工況三和工況四的實驗中, 上層流體的流速值較大, 超出流場測量系統(tǒng)的測量范圍, 造成監(jiān)測數(shù)據(jù)失真。但是從工況一、二的對比以及工況三、四下層流體中流速的對比可得出結(jié)論: 數(shù)值模擬結(jié)果與實驗結(jié)果吻合較好。

由圖6可看出, 隨時間的推移, 內(nèi)孤立波波致水平流速在上下層流體中的方向相反, 上層流體中的速度與波形傳播的方向一致, 下層流體中的速度與波形傳播方向相反; 上下層流體中的水平速度均呈現(xiàn)先增大后減小的趨勢; 在波谷經(jīng)過的時刻, 上下層流體中的水平速度均達到最大, 且上層流體的速度最大值大于下層流體。

表2為四個工況下, 由數(shù)值造波得到的上下層流體中間位置處內(nèi)孤立波波致水平流速的最大值。由表2而得出結(jié)論: 隨振幅的增加, 波致水平流速逐漸增大。

表2 內(nèi)孤立波波致水平流速的最大值Tab.2 The maximum horizontal velocity induced by internal solitary waves

3.2 波致水平流速沿垂向分布

圖 7為各工況下波谷經(jīng)過斷面處波致水平流速沿垂向的分布, 并與實驗結(jié)果進行對比。其中工況一、二、三中, 數(shù)值模擬結(jié)果與實驗結(jié)果一致; 在工況四中, 由于上層流體中波致水平流速較大, 致使實驗監(jiān)測失效, 因此實驗結(jié)果失真。

從圖中可以看出: (1)內(nèi)孤立波波致水平流速在上下層流體中的方向相反; (2)內(nèi)孤立波波致水平流速在上層流體中沿垂向分布較均勻, 在波面以下的下層流體中沿垂向分布有衰減的趨勢, 但衰減很小,而在兩層流體交界面與波谷之間的高度范圍內(nèi), 水平流速由正變負, 衰減明顯; (3)上層流體中的波致水平流速大于下層流體; (4)定義兩層流體界面與波谷的水深范圍為過渡水深范圍, 水平流速在該水深范圍內(nèi)沿垂向衰減明顯, 首先由正向衰減至零, 繼而反向增大。(5)對比不同振幅下內(nèi)孤立波波致水平流速, 發(fā)現(xiàn)隨內(nèi)孤立波振幅的增大, 過渡水深范圍有所增大。

圖6 內(nèi)孤立波波致水平流速沿垂向分布Fig.6 Profiles of magnitude of the horizontal velocity induced by internal solitary waves

4 結(jié)論

利用 Fluent軟件, 采用“平板拍擊”方法, 利用動網(wǎng)格技術(shù), 基于UDF編譯功能, 使用KdV、mKdV理論控制造波板的運動速度, 對內(nèi)孤立波進行數(shù)值模擬, 并研究了內(nèi)孤立波波致流場特性。研究結(jié)果表明:

(1) 該方法可以成功模擬內(nèi)孤立波的生成和傳播, 其生成的內(nèi)孤立波波形與實驗結(jié)果一致。

(2) 在內(nèi)孤立波傳播過程中, 上層流體的流速始終與內(nèi)孤立波傳播方向一致, 下層流體的流速始終與內(nèi)孤立波傳播方向相反, 說明內(nèi)孤立波在上、下層產(chǎn)生的流速是單方向的沖擊流, 而非振蕩流, 且流速的方向在上、下層始終相反形成剪切流。

(3) 對于一個空間固定位置, 內(nèi)孤立波經(jīng)過時,上下層流體的波致水平流速先增大后減小, 在波谷經(jīng)過時, 波致水平流速最大且上層流體中的水平流速大于下層流體。

(4) 波谷經(jīng)過斷面處的內(nèi)孤立波波致水平流速,在上層流體中沿垂向接近均勻分布, 在波面以下的下層流體中沿水深有較小的衰減趨勢。而兩層流體界面與波谷的水深范圍為過渡水深范圍, 水平流速在該水深范圍內(nèi)沿垂向衰減明顯, 首先由正向衰減至零, 繼而反向增大。

(5) 比較不同振幅下的內(nèi)孤立波波致水平流速, 發(fā)現(xiàn)隨內(nèi)孤立波振幅的增大, 過渡水深范圍有所增大。

方欣華, 杜 濤, 2005. 海洋內(nèi)波基礎(chǔ)和中國海內(nèi)波. 青島:中國海洋大學(xué)出版社, 258—281

孫大鵬, 李玉成, 2000. 數(shù)值水槽內(nèi)的阻尼消波和波浪變形計算. 海洋工程, 18(2): 46—50

杜 濤, 吳 巍, 方欣華, 2001. 海洋內(nèi)波的產(chǎn)生與分布. 海洋科學(xué), 25(4): 25—28

陳 鈺, 朱良生, 2009. 基于Fluent的海洋內(nèi)孤立波數(shù)值水槽模擬. 海洋技術(shù), 28(4): 72—75, 100

柯自明, 尹寶樹, 徐振華等, 2009. 南海文昌海域內(nèi)孤立波特征觀測研究. 海洋與湖沼, 40(3): 269—274

徐鑫哲, 2012. 內(nèi)波生成機理及二維內(nèi)波數(shù)值水槽模型研究.哈爾濱: 哈爾濱工程大學(xué)碩士學(xué)位論文, 32—34

高原雪, 尤云祥, 王 旭等, 2012. 基于MCC理論的內(nèi)孤立波數(shù)值模擬. 海洋工程, 30(4): 29—36

黃文昊, 尤云祥, 王 旭等, 2013. 有限深兩層流體中內(nèi)孤立波造波實驗及其理論模型. 物理學(xué)報, 62(8): 084705

董 志, 詹杰民, 2009. 基于VOF方法的數(shù)值波浪水槽以及造波、消波方法研究. 水動力學(xué)研究進展, 24(1): 15—21

韓 朋, 任 冰, 李雪臨等, 2009. 基于VOF方法的不規(guī)則波數(shù)值波浪水槽的阻尼消波研究. 水道港口, 30(1): 9—13

Buick J M, Martin A J, Cosgrove J A et al, 2003. Comparison of a lattice Boltzmann simulation of steep internal waves and laboratory measurements using particle image velocimetry.European Journal of Mechanics-B/Fluids, 22(1): 27—38

Choi W, Camassa R, 1999. Fully nonlinear internal waves in a two-fluid system. Journal of Fluid Mechanics, 396: 1—36

Helfrich K R, Melville W K, 2006. Long nonlinear internal waves. Annual Review of Fluid Mechanics, 38(1): 395—425

Liu B Q, Yang H, Zhao Z X et al, 2014. Internal solitary wave propagation observed by tandem satellites. Geophysical Research Letters, 41(6): 2077—2085

Michallet H, Barthélemy E, 1998. Experimental study of interfacial solitary waves. Journal of Fluid Mechanics, 366:159—177

Osborne A R, Burch T L, 1980. Internal solitons in the Andaman Sea. Science, 208(4443): 451—460

Wang F, Wu K F, Guo H Y et al, 2015. Experimental study on internal solitary wave induced flow field. In: Xie L Q ed.Advanced Engineering and Technology II. Hong Kong,China: CRC Press, 111

Xu Z H, Yin B S, Hou Y J, 2010. Highly nonlinear internal solitary waves over the continental shelf of the northwestern South China Sea. Chinese Journal of Oceanology and Limnology, 28(5): 1049—1054

Zhang L, Wang L L, Yu Z Z et al, 2012. Characteristics of nonlinear internal waves in a three-dimensional numerical wave tank. Applied Mechanics and Materials, 212—213: 1123—1130