內孤立波數值模擬方法研究?

苗得勝， 郭海燕， 趙　婧， 王　飛， 王　偉

(1.中國海洋大學工程學院，山東 青島 266100； 2.青島酒店管理職業技術學院，山東 青島 266100)

?

內孤立波數值模擬方法研究?

苗得勝1， 郭海燕1， 趙婧2， 王飛1， 王偉1

(1.中國海洋大學工程學院，山東 青島 266100； 2.青島酒店管理職業技術學院，山東 青島 266100)

本文選用RNGk-ε湍流模型封閉二維不可壓縮N-S方程組，選用KdV、mKdV理論作為輸入條件進行內孤立波造波模擬。運用流體模擬軟件Fluent建立二維內波數值水槽，用動網格手段實現網格變形模擬物理邊界造波法，用VOF方法進行密度分層流體的內界面追蹤，對不同非線性強度的二相流內孤立波進行了數值模擬。并將數值模擬結果與KdV、mKdV理論波形以及內波試驗水槽所造內孤立波進行對比，驗證了此套造波方法在Fluent中進行內孤立波模擬的可行性。

內孤立波；Fluent；動網格；造波方法；數值模擬

引用格式：苗得勝， 郭海燕， 趙婧， 等. 內孤立波數值模擬方法研究 [J].中國海洋大學學報(自然科學版), 2016, 46(10):123-128.

MIAO De-Sheng， GUO Hai-Yan， ZHAO Jing， et al. Study of numerical simulation method of internal solitary waves [J].Periodical of Ocean University of China, 2016, 46(10):123-128.

內波是發生在密度穩定的分層水體中、最大振幅出現在水體內部的波動。1902年F. Nansen通過研究船舶在遇冰溶解后的水面上航行時速度突然變慢的現象首次發現了海洋內波現象[1]。1970年代末美國首次利用衛星觀測內波的合成孔徑雷達(SAR)影像,開辟了海洋內波研究的新時代[2]。海洋內波振幅出現在海洋內部,活動不易察覺。密度變化面上產生的內波對海洋結構物及人員的安全構成了巨大威脅。我國南海水深、風大、浪高、海水密度層化嚴重，海洋內波活動頻繁，南海海域頻繁發生的內波活動已經成為影響油氣田開發工程結構的災害性因素。這就決定了內波研究在工程中具有重要的現實意義。

在工程結構內波作用研究中關鍵問題是獲取內波數據, 獲取內波數據的途徑主要有實地觀測、實驗室造波、理論研究以及數值模擬造波技術。徐肇延等人在分層流體內波水槽中進行了內波產生、內波尾跡以及內波增阻等方面的試驗[3]。KdV和mKdV內波理論是目前進行理論研究比較常用的兩個內孤立波理論。尤云祥、李巍[4]等運用CFD方法建立內波數值水槽，研究內波作用。實地觀測內波費用昂貴,實驗室造波有尺寸限制,而數值模擬耗費小、用時短,沒有尺寸限制,同時可以根據需求調整參數來制造符合要求的內波,因此越來越多的內波研究采用數值模擬的方法進行。

目前進行造波數值模擬常用的工具軟件主要有Fluent和ANSYS中的CFX模塊。兩者均采用有限體積法，并可通過C語言或C++等進行二次開發。一般而言，CFX在求解算法上的優勢使得其在求解收斂速度上有一定優勢，Fluent在誤差控制和算法穩定性上更有優勢。另外，Fluent具有網格重構功能，處理動邊界問題有較大的優勢。本文的水槽模型涉及到邊界運動問題，因此本文采用Fluent作為工具來建立數值水槽對內孤立波進行模擬。

本文的主要工作：選用恰當的湍流模型來封閉N-S方程組，推導造波推板運動方程，選取合適的內孤立波理論寫入造波推板運動方程函數中作為輸入條件，然后用Fluent建立內波數值水槽，通過用戶自定義函數(UDF)將造波推板運動程序導入，用動網格手段模擬物理造波法，用VOF方法進行界面追蹤，模擬兩層密度分層水體中內孤立波的產生與傳播。最后將數值造波結果與理論波形以及項目組在內波水槽所做內孤立波試驗結果進行對比，驗證Fluent數值造波技術的可行性。

1　Fluent造波模擬

本文將三維內波問題簡化為二維平面內波問題，將兩層流體均視為有黏不可壓縮流體，選用二維不可壓縮Navier-Stokes方程作為流體控制方程。

Fluent進行內孤立波波模擬主要有水平推板造波法、雙推板造波法、速度入口法和源項造波法。雙推板造波法容易造成流體上表面波動。速度入口法在計算過程中由于誤差的積累會導致進出流量不平衡而計算出錯或造波結果誤差較大。水平推板造波法不會引起流體上表面的明顯波動，同時避免了速度入口法進出水量不平衡的問題。因此本文采用水平推板法造波。

水槽參數：為方便跟試驗進行對比，本文的數值水槽參數參照中國海洋大學內波試驗水槽進行選取。內波試驗水槽總高度60cm，長1500cm，額定水深58cm。數值模擬中水槽高度參照試驗水槽的額定水深取58cm。其他尺寸同內波試驗水槽。

本文在參考運用水平推板法進行內周期波數值模擬的基礎上，對水平推板的長度、放置的位置以及隔板的設置進行改進，實現對內孤立波的數值模擬。考慮所造內孤立波波形，取造波推板長度為100cm。為防止造波推板上下運動引起上下層流體的摻混，在造波板右側加豎向隔板。調整隔板的高度保證其不會影響內孤立波的生成。綜合考慮分層流體的厚度、內孤立波波高和造波推板運動的范圍，取隔板長度為19cm，隔板最上端高出流體內分界面1cm。水槽布置見圖1。

圖1　水槽示意圖Fig.1　Sketch of flume

邊界條件設定設置水槽上邊界為壓力入口邊界。下邊界和左右邊界均設為無滑移固壁邊界。隔板設為有滑移固壁邊界。設置水平造波推板為移動邊界。

造波推板控制推板造波法的原理是通過控制推板上下移動，實現分層流體的反向流動，從而激發內波。為保證所造波形與理論波形具有對比性，需選取合適的理論寫入用戶自定義函數(UDF)來控制造波推板的實時運動。目前，內孤立波理論模型主要有KdV、eKdV、mKdV 和 MCC等。Helfrich和Melville[5]通過實驗研究發現KdV理論對于波高水深比小于0.1的兩層流體模型適用性較好。尤云翔等的研究表明mKdV模型適用于波高水深比大于0.1的內孤立波模型[6]。因此本文選取這兩種內孤立波理論作為輸入條件。

KdV理論將實際海洋中的密度分層簡化為兩層模型，取上下層流體深度和密度分別為h1、ρ1和h2、ρ2，總水深為h。內孤立波沿正向傳播，振幅為η0，波面方程的KdV理論解[4]為：

(1)

(2)

式中：

(3)

波形為下凹時η0取負值，上凸時η0取正值。

內孤立波的mKdV理論解[5]為：

(4)

其中：

(5)

(6)

(7)

(8)

波形為下凹時η0取負值,上凸時η0取正值。

(9)

因為：

Δx=cΔt

(10)

可得推板運動速度：

(11)

將KdV理論解的波面方程(1)和波速方程(2)帶入(11)式,得到KdV作為輸入條件的UDF。將mKdV理論解的波面方程(4)和波速方程(9)帶入(11)式，得到mKdV作為輸入條件的UDF。

網格劃分在造波區，由于造波板的運動要用到動網格，為保證網格更新質量，此部分長寬方向劃分較細的均勻網格。在工作區需要準確捕捉波面位置，因此將波面經過的區域在垂向方向進行局部加密。下層流體厚度較大，底部速度很小，為了兼顧計算速度，從加密區到水槽底部采用逐漸增大的漸變網格，工作區網格劃分如圖2所示。

圖2　工作區網格劃分Fig.2　Mesh of workspace

湍流模型選取湍流模型的選擇對內孤立波傳播過程中能量的衰減有重要影響。針對本文造波工況，選用k-ε湍流模型來封閉N-S方程組。Fluent提供了3種k-ε模型，其中標準k-ε模型魯棒性最好，適用于初始迭代、參數選型和參數研究；Realizablek-ε模型計算精度較高，耗時相對長，適用于漩渦較大的工況；RNGk-ε模型在處理高應變率及流線彎曲程度較大的流動時表現較好，在結構較為復雜時計算結果較好。本文所造內孤立波流線彎曲程度較大，故選取RNGk-ε湍流模型來封閉N-S方程組。

波面追蹤方法本文采用VOF方法[7]對兩層流體的分界面進行追蹤。VOF模型是一種固定的歐拉網格下的表面追蹤方法，通過求解一組動量方程和追蹤計算域中每一相流體的體積分數來模擬兩種互不相溶的流體之間的自由界面。構造網格體積函數F，若在某時刻網格單元中F=l，則說明該單元全部為指定相流體所占據，為流體單元。若F=0，則該單元全部為另一相流體所占據，為空單元。當0

求解設置建好模型后導入Fluent進行求解設置。選擇基于壓力的二階隱式時間離散的非穩態求解器。定義基本相流體密度為1035kg/m3，第二相流體密度為1055kg/m，兩者動力黏度均取1.003e-03kg/m-s。選用動網格模型來實現造波板的運動，網格重構方法選用動態層模型，網格分割因子取0.4。壓力速度耦合方式選PISO。壓力插值選擇體積力加權法。對流相及輸運方程采用一階迎風格式進行離散。

消波方法本文采用阻尼消波法進行消波。通過在消波區動量方程中添加附加源項來改變消波區流體的阻尼進行消波?？紤]內孤立波波長，消波區長度取400cm。為了避免流場內粘度突變，消波區阻尼采用線性漸變的增大方式。修改后的消波區動量方程[8]如下：

(12)

其中：

(13)

其中a為阻尼消波系數，本文取a=8.0。

求解計算設置好求解控制條件后進行迭代求解社設置。迭代時間步長取0.005s，為保證計算精度，每個時間步內最大迭代次數40次。然后進行迭代求解。

圖3是運用上述方法進行數值模擬得到的內孤立波波形圖，橫軸為內孤立波傳播正方向，水平向左。縱軸為水槽深度方向，豎直向上。內孤立波波形光滑，對稱性良好。

圖3　內孤立波波形圖Fig.3　Shape of internal solitary wave

圖4是數值模擬得到的內孤立波流場速度矢量圖，箭頭方向為流體質點瞬時速度指向，長度為速度大小。此圖描述了內孤立波傳播過程中波谷周圍的流速分布情況。橫軸水平向左，縱軸豎直向上。從圖中可以看出，在振幅最大值附近，上層流體水平流速與內波傳播方向相同，流速較大，下層流體水平流速與內波傳播方向相反，流速相對較小。

圖4　內孤立波流場速度矢量圖Fig.4　Velocity vector of internal solitary wave

2　結果分析

運用Fluent進行內孤立波造波模擬得到內孤立波結果的準確性需要用等比尺的內孤立波試驗進行驗證。

在中國海洋大學物理海洋實驗室內波試驗水槽中進行等比尺內孤立波造波試驗。水槽高60cm，寬35cm，長1500cm，額定水深58cm。淡水密度1035kg/m3，厚度9.5cm，鹽水密度1055kg/m3，厚度48.5cm。

此次試驗采用重力塌陷法制造內孤立波。塌陷區長40cm。為得到滿足KdV理論和mKdV理論適用范圍的不同非線性強度內孤立波，設置4個試驗工況進行試驗。四個工況分別取塌陷區流體內分界面與工作區流體內分界面之間的高度差10，15，20和25cm。由此得到四個不同波高的內孤立波，波高分別為4.060、5.877、7.378和8.108cm。

試驗過程中抽離隔板后內液面差造成的擾動并未激起可見表面波，因此不考慮其對內波試驗的影響。水槽依靠兩塊在末端傾斜放置的多孔鐵板進行消波，大部分內波能量在經過多孔鐵板時被耗去，余下的小振幅反射波傳回工作區時各項數據的觀測工作已經完成，因此不考慮反射波對內波試驗的影響。

將試驗工況1和工況2得到的小波高水深比的內孤立波波高作為輸入條件寫入基于KdV理論編寫的UDF函數，導入Fluent控制水平推板的運動；將試驗工況3和工況4得到的波高水深比大于0.1的內孤立波波高作為輸入條件導入Fluent控制水平推板的運動；進行造波模擬。得到4個工況下的波高結果見表1。

表1　數值模擬與試驗波高對比Table 1　wave height contrast of numerical solution and experimental solution

Note：①Experimental wave height；②Numerical wave height；③Percentage error

分析表1，前3組數值解波高略大于試驗值，工況4數值解波高略小于試驗值。四組數據的誤差均在5%以內，數值模擬結果較好。

為驗證數值造波波形的準確性，將數值解波形與實驗結果以及理論結果進行對比。對4個工況得到的數值解和試驗結果波形監測數據進行數據處理并導入Matlab作圖。同時將試驗得到的4個波高帶入KdV和mKdV方程得到4個工況對應的理論波面解，將三者放到一起進行對比，如圖5～圖8所示。橫軸為時間軸，縱軸為波高水深比。

圖5　工況1試驗結果、數值解及理論解對比Fig.5　Contrast of experimental solution，numerical solution and theoretical solution in condition 1

從圖5可以看出，數值模擬所造內孤立波(KdV)波形光滑、尾波幅度較小且穩定，試驗所造內波不太光滑，對稱性稍差。理論波形光滑度和對稱性都很好。試驗結果與數值解在波形上吻合度很高，兩者均比KdV理論波形稍寬。由于工況1的波高水深比較為接近KdV理論最適用的波高水深比，此時得出的理論波形基本與試驗結果和理論解吻合，因此圖中3個波只有細微的差異?？紤]到表1中工況1試驗波高和理論波高要比數值模擬波高小5%，因此數值模擬結果得到的波形與試驗結果和理論解的吻合程度比圖中展示出的吻合度更高。

分析圖6可知，數值模擬所造內孤立波(KdV)與試驗所造內波在波形上吻合度較高，數值解波形略寬。KdV理論波形比數值解和試驗結果波形要窄，偏離的程度較為明顯。這是由于工況2的波高水深比接近0.1，此時KdV理論對應的理論波形波長偏小，所以做出的波形圖會顯示出比數值解和試驗結果稍窄的情況，因此展示出如圖所示的差異。

分析圖7可知，數值模擬所造內孤立波(mKdV)與試驗所造內波在波形上吻合度很高，數值解波形更為平滑，在尾端更符合理論解。試驗結果尾部波形失真較大，對比前幾組試驗，原因是測量誤差。mKdV理論波形比數值解和試驗波形寬一些，這與上兩組工況采用KdV理論的造波結果剛好相反。由于工況3中內孤立波的波高水深比接近0.1，此時的mKdV理論波形波長偏大，因此圖中會顯示出數值解與試驗波形高度吻合而兩者均比理論波形稍窄的情況。

圖6　工況2試驗結果、數值解及理論解對比Fig.6　Contrast of experimental solution，numerical solution and theoretical solution in condition 2

圖7　工況3試驗結果、數值解及理論解對比Fig.7　Contrast of experimental solution，numerical solution and theoretical solution in condition 3

分析圖8可知，數值擬所造內孤立波(mKdV)波形對稱性良好，波形光滑。試驗結果波形則表現出較強的不對稱性，波形也較為粗糙。這是因為工況4在最后一組進行試驗，此時兩層液面的分界面因相互滲透而變得模糊，導致數據采集出現較大誤差。數值解大體上與試驗結果吻合，而mKdV理論波形比數值解和試驗結果都要寬。此差異產生原因與工況3中理論波形與數值解和試驗波形存在差異的原因相同。

圖8　工況4試驗結果、數值解及理論解對比Fig.8　Contrast of experimental solution，numerical solution and theoretical solution in condition 4

3　結論

(1)本文用Fluent軟件建立數值水槽，選取KdV、mKdV理論作為輸入條件，運用改進的水平推板法實現了內孤立波的數值模擬。將在多個工況下得到的內孤立波與理論波形以及試驗結果進行了對比，發現數值模擬結果與試驗結果吻合度較高，驗證了選用RNGk-ε湍流模型封閉N-S方程組進行內孤立波數值造波的適用性和準確性。

(2)對比發現，工況1、2數值模擬波形比KdV理論波形稍寬，工況3、4數值模擬波形比mKdV理論波形稍窄。4個工況均采用理論方程解作為輸入條件，得到的波形與試驗得到的波形更為吻合，說明內波在傳播過程中發生了演化并最終達到一個穩定狀態，說明內孤立波最終的波形取決于流場自身而非外在激勵。同時，試驗結果與數值結果的吻合度明顯高于兩者與理論值的吻合度，說明KdV理論和mKdV理論的波面方程在本文選定波高水深比的工況下存在一定的誤差。

[1]方欣華, 杜濤. 海洋內波基礎和中國海內波[M].青島：中國海洋大學出版社，2004: 69-281.

Fang X H，Du T. Fundamentals of Oceanic Internal Waves and Internal Waves in the China Seas[M]. Qingdao：China Ocean University Press， 2004: 69-281.

[2]魏崗. 分層流體中運動潛體生成的內波以及內波的垂向結構研究[D]. 上海：上海大學, 2003.

Wei G. Waves Generated by a Submerged Body Moving in Stratified Fluids and Vertical Structures of Internal Waves[D]. Shanghai: Shanghai University， 2003.

[3]徐肇廷, 姚鳳朝, 隋紅波. 分層海洋中運動物體生成的內波[J]. 中國海洋大學學報(自然科學版)，2001, 31(4): 461-466.

Xu Z T，Yao F C，Sui H B. Internal waves generated by moving body in the stratified ocean[J]. Journal of Ocean University of Qingdao，2001， 31(4): 461-466.

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

Gao Y X, You Y X, Wang X，et al. Numerical simulation for the internal solitary wave based on MCC theory[J].The Ocean Engineering，2012(4): 29-36.

[5]Helfrich K R，Melville W K. On long nonlinear internal waves over slope-shelf topography [J]. Journal of Fluid Mechanics， 1986， 167: 285-308.

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

Hen W H，You Y X，Wang X，et al. Wave-making experiments and theoretical models for internal solitary waves in a two-layer fluid of finite depth[J]. Acta Phys， 2013， 62(8)： 084705： 12-13.

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

Dong Z， Zhan J M. Comparison of existing methods for wave generating and absorbing in VOF-based numerical tank[J]. Journal of Hydrodynamics， 2009(1): 15-21.

[8]徐鑫哲. 內波生成機理及二維內波數值水槽模型研究[D]. 哈爾濱：哈爾濱工程大學, 2012.

Xu X Z. Study on Generational Mechanism and 2-D Numerical Flume Model of Internal Waves[D]. Harbin：Harbin Engineering University， 2012.

[9]Wessels F，Hutter K. Interaction of internal waves with a topographic still in a two-layered fluid[J]. J Phys Ocean Ogr， 1996， 26 (1): 5-20.

[10]Walker S A， Martin A J， Easson W J， et al. Comparison of laboratory and theoretical internal solitary wave kinematics[J]. Journal of waterway， port， Coastal and Ocean Engineering， 2003， 129(5): 210-218.

[11]Cai S Q， Gan Z J， Long X M. Some characteristics and evolution of the internal soliton in the northern South China Sea[J]. Chinese Science Bulletin， 2002， 47(1): 21-26.

責任編輯陳呈超

Study of Numerical Simulation Method of Internal Solitary Waves

MIAO De-Sheng1， GUO Hai-Yan1， ZHAO Jing2， WANG Fei1， WANG Wei1

(1.College of Engineering， Ocean University of China， Qingdao 266100， China； 2.Qingdao Vocational and Technical College of Hotel Management, Qingdao 266100， China)

A two-dimensional flume numerical simulation model is established with Fluent to study generation and propagation of internal solitary wave. RNG k-epsilon turbulent model is adopted to close Navier-Stokes equations. Internal solitary wave is generated by wave-making board which is controlled by user code written in C language based on KdV and mKdV theories. Dynamic grid method is adopted to adapt the grid deformation around the wave-making board during the simulation. VOF method is adopted to catch the interface between density-stratified fluids, which represents the waveform of internal solitary wave. Internal solitary waves with diverse nonlinear strength are simulated in this paper to study the application range of KdV and mKdV theories. Data comparison of numerical solution, theoretical solution and experiment solution is carried out after the simulation. The feasibility and accuracy of internal solitary wave simulation model raised in this paper is verified through the comparison. The application range of KdV and mKdV theories is also discussed based on the comparison.

internal solitary wave; Fluent; dynamic grid; VOF method; numerical simulation

國家自然科學基金項目(51279187)；中央高?；究蒲袠I務費項目(201262005)

2014-10-11；

2015-07-10

苗得勝(1989-)，男，碩士。E-mail:mdsouc@163.com

P751

A

1672-5174(2016)10-123-06

10.16441/j.cnki.hdxb.20140311

Supported by the Natural Science Foundation of China(51279187);the Fundamental Research Funds for the Central Universities(201262005)