内孤立波数值模拟方法研究❋

苗得胜， 郭海燕， 赵　婧， 王　飞， 王　伟

(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)；中央高校基本科研业务费项目(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)