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

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

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

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

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

1 建立數值水槽

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

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

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

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

式中:

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

式中:

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

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

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

圖2 數值水槽網格劃分示意圖Fig.2 The gridding of the numerical flume

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

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

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

表1 數值模擬工況設置Tab.1 Setting of amplitude for the numerical simulation

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

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

2 數值模擬結果分析

2.1 內波物理實驗方法

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

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

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

2.2 數值造波結果驗證

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

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

綜上, 本文建立的數值水槽能夠實現設計振幅的內孤立波數值造波, 且數值結果與實驗結果較吻合, 能有效實現對內孤立波的數值模擬。

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

3 討論

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

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

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

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

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

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

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

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

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

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

4 結論

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

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

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

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

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

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

方欣華, 杜 濤, 2005. 海洋內波基礎和中國海內波. 青島:中國海洋大學出版社, 258—281

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

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

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

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

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

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

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

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

韓 朋, 任 冰, 李雪臨等, 2009. 基于VOF方法的不規(guī)則波數值波浪水槽的阻尼消波研究. 水道港口, 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