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

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

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

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

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

1 建立數(shù)值水槽

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

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

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

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

式中:

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

式中:

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

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

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

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

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

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

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

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

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

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

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

2.1 內(nèi)波物理實(shí)驗(yàn)方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

從圖中可以看出: (1)內(nèi)孤立波波致水平流速在上下層流體中的方向相反; (2)內(nèi)孤立波波致水平流速在上層流體中沿垂向分布較均勻, 在波面以下的下層流體中沿垂向分布有衰減的趨勢, 但衰減很小,而在兩層流體交界面與波谷之間的高度范圍內(nèi), 水平流速由正變負(fù), 衰減明顯; (3)上層流體中的波致水平流速大于下層流體; (4)定義兩層流體界面與波谷的水深范圍為過渡水深范圍, 水平流速在該水深范圍內(nèi)沿垂向衰減明顯, 首先由正向衰減至零, 繼而反向增大。(5)對(duì)比不同振幅下內(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軟件, 采用“平板拍擊”方法, 利用動(dòng)網(wǎng)格技術(shù), 基于UDF編譯功能, 使用KdV、mKdV理論控制造波板的運(yùn)動(dòng)速度, 對(duì)內(nèi)孤立波進(jìn)行數(shù)值模擬, 并研究了內(nèi)孤立波波致流場特性。研究結(jié)果表明:

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

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

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

(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)的阻尼消波和波浪變形計(jì)算. 海洋工程, 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)波生成機(jī)理及二維內(nèi)波數(shù)值水槽模型研究.哈爾濱: 哈爾濱工程大學(xué)碩士學(xué)位論文, 32—34

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

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

董 志, 詹杰民, 2009. 基于VOF方法的數(shù)值波浪水槽以及造波、消波方法研究. 水動(dòng)力學(xué)研究進(jìn)展, 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