内孤立波数值造波方法及其与理论和实验结果的比较*

李效民　张　林　郭海燕　姜　海　王　伟

(中国海洋大学工程学院　青岛　266100)

经过长期的开采,陆上油气资源日益减少,人们把目光转向了资源丰富的海洋; 而复杂的海洋环境存在许多制约海洋开发的因素,内波即是其中之一。内波会对海洋立管、海洋平台等海洋结构物产生较强的破坏作用,给海洋油气的开发带来极大的挑战。

海洋内波是发生在密度稳定层化的海水内部的一种波动。海洋内波的频率和波速远小于表面波,但波长和振幅远大于表面波。据记载,最大垂向振幅高达180m(方欣华等,2005)。内孤立波是一种强非线性潮成内波,主要通过非线性效应与色散效应达到一定程度的平衡,实现稳定传播(Osborneet al,1980)。

目前广泛使用的内孤立波理论模型有 KdV理论、eKdV理论、mKdV理论、MCC理论等。研究内孤立波的手段有理论分析、实地观测及实测图像分析、物理实验、数值模拟,其中,物理实验和数值模拟实施简便、结果直观可靠、适用性较强,深受研究者的青睐,而数值模拟更是凭借其成本低、可重复性好、快速的优势成为最常用的内孤立波研究方法。Wei等(2009)采用质量源数值造波方法模拟了两层流体中内孤立波的生成与传播。付东明等(2009)结合KdV和mKdV理论发展了双推板数值造波方法。陈钰等(2009)利用 Fluent软件,采用速度入射边界法,成功建立了模拟弱非线性内孤立波的分层数值水槽。王旭等(2014)结合内孤立波理论进一步研究了质量源数值造波方法,与理论结果较符合。Zhang等(2012)使用 CFD方法建立了三维数值水槽,所得结果与KdV理论结果吻合较好。高原雪等(2012)基于 MCC理论,采用速度入射边界法建立数值水槽,实现了振幅可控的内孤立波数值模拟。目前国内外主要将数值模拟结果与理论结果进行对比分析,而将理论、实验、数值模拟三者进行对比分析的研究较少。此外,对各种数值造波方法优劣的比较也较为缺乏。

本文针对三种内孤立波理论模型与数值造波方法,基于KdV、mKdV和eKdV三种理论模型与双推板、平板拍击和速度入射边界三种数值造波方法所组成的一系列工况,利用Fluent软件进行计算,并将数值模拟结果与理论公式及实验结果进行对比分析,验证了数值造波的合理性,评估了这三种数值造波方法。

1　内孤立波理论模型

采用有限水深两层流体简化模型,假设上下层流体为无旋且不可压缩的理想流体,上层流体密度为ρ1,深度为h1,下层流体密度为ρ2,深度为h2,总水深为h=h1+h2。如图1,建立直角坐标系oxyz。其中oxy平面为未扰动的两层流体界面,ox轴正向为内孤立波传播方向,oz轴垂直于oxy面向上。令η(x,t)表示内孤立波的波面函数;为振幅大小,符号满足上凸波为正,下凹波为负;λ为内孤立波特征波长;c为内孤立波相速度。

图1　内孤立波及其参数Fig.1　Schema of internal solitary wave and the parameterization

1.1　KdV理论

波面方程的KdV理论解为(Choiet al,1999):

其中:

式中cKdV为 KdV理论下内孤立波的相速度;lKdV为KdV理论下内孤立波的特征长度。

KdV理论不存在极限振幅。但实验研究发现,振幅与水深比小于0.1时,KdV理论完全适用,而对于大振幅内孤立波不再适用。

1.2　eKdV理论

波面方程的eKdV理论解为(Helfrichet al,2006):

其中:

式中ceKdV为eKdV理论下内孤立波的相速度。

eKdV理论适用于振幅与水深比超过0.1的广泛振幅范围,且存在极限振幅,其大小为:

1.3　mKdV理论

波面方程的 mKdV理论解为(Michalletet al,1998):

其中:

式中cmKdV为mKdV理论下内孤立波的相速度;为极限振幅大小。

mKdV理论适用于振幅与水深比超过0.1且上下层水深比接近临界水深比的情况。

2　数值造波方法

目前,数值造波方法可大致归为两类: 一类是仿物理造波法,该方法源于实验室造波技术,即仿造物理造波机的工作原理。如: 推板法、平板拍击法等。另一类是纯数值造波法,如: 速度入射边界法等。

2.1　双推板法

图2为双推板法的示意图,水槽上、下边界及右边界设置为固壁边界,左侧上、下推板设置为移动边界。其造波原理为: 选取合适的内孤立波理论,采用动网格技术,通过用户自定义函数(UDF)的DEFINE_CG_MOTION宏控制两块推板的水平反向运动实现两层流体界面处内孤立波的模拟。

图2　双推板法造波示意图Fig.2　Wave generation by double push-pedals method

令上板高度为d1,下板高度为d2,上下板的水平运动速度分别为:U1、U2。根据上推板运动排开的流体体积与两层流体界面处水位在水平方向的积分相等,即:

可得上推板运动速度为(徐鑫哲,2012):

再根据流体质量守恒(徐鑫哲,2012),得下推板运动速度为:

式中c为内孤立波的相速度,η(x,t)为内孤立波波面函数(上凸波为正,下凹波为负)。

实际设置时,应使d1略大于以防上下层流体掺混(付东明等,2009)。

2.2　平板拍击法

图3为平板拍击法的示意图,水槽上、下、左、右边界及隔板均设置为固壁边界,造波板设置为移动边界。其造波原理为: 选取合适的内孤立波理论,采用动网格技术,通过用户自定义函数(UDF)的DEFINE_CG_MOTION宏控制造波板做上下运动实现两层流体界面处内孤立波的模拟。

图3　平板拍击法造波示意图Fig.3　Wave generation by moving pedal method

令造波板的长度为D,运动速度为U。根据造波板移动所排开的流体体积与造波板右端两层流体界面处水位在水平方向的积分相等,即:

可得造波板运动速度为(徐鑫哲,2012):

2.3　速度入射边界法

图4为速度入射边界法的造波示意图,水槽上、下边界及右边界设置为固壁边界,左侧边界设置为速度入射边界。该方法是以边界为扰动源,但边界是固定不运动的。其造波原理为: 选取合适的内孤立波理论,通过用户自定义函数(UDF)的 DEFINE_PROFILE宏在入射边界上给定该理论下的水流速度实现两层流体界面处内孤立波的模拟。

图4　速度入射边界法造波示意图Fig.4　Wave generation by velocity-inlet boundary method

入射边界处水流速度设为:

3　数值模拟与结果分析

3.1　内孤立波物理实验

王飞(2015)在中国海洋大学物理海洋实验室分层流体水槽中进行了内孤立波物理实验。水槽长15m,宽0.35m,高0.7m,额定水深0.6m。采用Oster双缸法制取密度分层流体,上层流体高度 9.5cm,密度1035kg/m3,下层流体高度48.5cm,密度1054kg/m3。采用重力塌陷法制造不同振幅的内孤立波,使用两个斜面在水槽末端进行消波,利用PIV技术对内孤立波致流场进行测量,使用高速相机进行图像采集,利用图像处理软件进行数据分析。物理实验结果如图6、图8所示。

3.2　建立数值水槽

为便于与实验结果对比,本文参照内波实验建立数值水槽模型。采用两层流体简化模型,利用Fluent进行内孤立波数值模拟。采用结构化网格离散,在两层流体交界面附近对网格划分进行加密。工作区网格划分如图5(a)所示: 底部向上38cm至50cm高度范围内,网格尺寸为0.5cm,其余为1cm,x轴方向网格尺寸为2cm。选用二维不可压缩Navier-Stokes方程作为流体控制方程,使用有限体积法对控制方程进行离散。选用κ-ε湍流模型,采用VOF方法追踪两层流体的交界面,两层流体界面的构造采用几何重构法,压力速度耦合选用 PISO算法,压力插值采用体积力加权法,对流相及输运方程的离散采用一阶迎风格式。

由于水槽的右边界为固壁边界,易产生反射波,本文采用阻尼消波法在水槽末端进行消波。消波段长度取为 1—2倍波长。同时,在消波区沿x轴方向采用逐渐增大的渐变网格来耗散波浪的能量。消波区网格划分如图 5(b)所示:z轴方向网格划分与工作区相同,x轴方向首层网格尺寸为 2cm,后续网格尺寸按1.04的比例逐渐变大。

图5　数值水槽网格划分示意图Fig.5　Grids of the numerical flume

根据内孤立波理论对振幅水深比的使用条件(黄文昊等,2013),按照不同的振幅,选取合适的理论模型,采用不同的数值造波方法进行数值模拟。具体工况如表1所示。

设置各项参数后进行流场的迭代计算,迭代时间步长取 0.01s,每个时间步长的最大迭代次数取 20次。数值计算结果如图6、图8所示。

3.3　计算结果与分析

图 6为各工况下数值模拟所得波形与相应理论及实验结果的对比。结果表明: 三种数值造波方法均能生成内孤立波,但模拟效果有差异。由工况 A1、A2、A3(B1、B2、B3 和 C1、C2、C3)可知,对于同一种振幅、同一种内孤立波理论,速度入射边界法所模拟波形与理论及实验波形吻合最好。平板拍击法所模拟波形与理论及实验波形吻合较好,只是在波形尾端存在较小的尾波,但在传播过程中尾波会逐渐衰减,不会对内波产生影响。双推板法所得模拟波形与理论及实验波形趋势相同,但最大振幅达不到设计振幅。

表1　数值模拟工况设置Tab.1　Condition setting of numerical simulation

由工况 A1、B1、C1(A2、B2、C2和 A3、B3、C3)可知,随着振幅的增加,速度入射边界法和平板拍击法所得波形与理论及实验波形吻合仍然较好。双推板所得波形与理论及实验波形趋势仍然一致,但最大振幅与设计振幅的差值越来越大,其原因为流体阻碍了推板的运动造成耗散。由图7(a)可知A1工况下加大推板的运动幅值后所得波形与理论及实验波形吻合较好。因此,可通过加大推板的运动幅值使所生成内孤立波的幅值满足要求。此外,双推板法数值造波所得幅值还受上下层流体深度比的影响。实验研究发现: 实测振幅均比设计振幅小,且两者之间存在线性关系(黄文昊等,2013),因此,实际使用时需拟合相同条件下不同设计振幅与所得振幅之间的关系,然后根据实际所需振幅反算出设计振幅。

图 6显示三种数值造波方法结果与实验结果存在差异,尤其随着波高的增大,两者的波形在末端差异越大。其原因有: (1)内波在传播过程中粘滞效应影响会导致能量损耗。(2)实验造波方法为重力塌陷法,波高与塌陷高度差呈正相关,塌陷高度差越大,对流体产生的激荡也越强,引发的尾波越明显。(3)实验的测量设备、边界条件等的限制导致实验结果与数值模拟结果不能完全一致。

图6　数值造波波形与理论及实验结果对比Fig.6　Comparison in numerical waveform between experimental and theoretical results

图7　A1工况下修正后的数值造波结果与理论及实验结果的对比Fig.7　Comparison of calibrated numerical result with experimental and theoretical ones for Case A1

由工况B1、B2、B3可知,mKdV理论波形明显要宽于实验波形与数值模拟波形。这一偏差说明mKdV理论对这种波高水深比情况下的内孤立波的描述存在一定偏差,但基于该理论的数值造波仍然能给出比较接近实际情况的波形。

图 8为各工况下数值模拟所得波谷经过断面处波致水平流速沿垂向的分布,并与理论及实验结果进行了对比。结果表明: 三种数值造波方法所得波致水平流速沿垂向分布趋势基本一致,但有所差异。由工况 A1、A2、A3(B1、B2、B3)可知,平板拍击法和速度入射边界法所得波致水平流速均与理论及实验结果吻合较好。双推板法所得波致水平流速相对于理论及实验结果明显偏小。而对于工况C1、C2、C3,由于实验时上层流体中波致水平流速较大,致使实验监测数据失效,实验结果失真。但仍然可以看出速度入射边界法和平板拍击法所得波致水平流速与理论结果吻合较好,双推板法与理论结果存在一定差异。

图8　内孤立波波致水平流速沿垂向分布Fig.8　Horizontal velocity magnitude induced by internal solitary waves along vertical section

由工况 A1、B1、C1(A2、B2、C2和 A3、B3、C3)可知,随着振幅的增加,速度入射边界法和平板拍击法所得结果与理论及实验结果吻合仍然较好。双推板法所得结果与理论及实验结果趋势相同,但上下层流体水平速度与理论及实验结果的差值越来越大。加大推板的运动幅值,所得内孤立波的上下层流体水平速度与理论及实验结果吻合较好,如图7(b)所示。

图8显示实验过程中通过PIV测量得到的上层流体水平流速变化情况比较混乱,没有规律,与数值模拟结果存在差异,其原因为上层流体中波致水平流速较大,超出了PIV系统的测量范围,故测得的数据存在一定失真。

内界面与波面之间区域的误差主要原因是实验与数值模拟存在密度跃层,密度跃层内会有流速变化,而理论模型是简单地分为两层处理,没有中间过渡阶段。

综合上述分析,可以发现三种造波方法均能生成内孤立波,但双推板法需拟合设计振幅与所需振幅之间的关系,通过加大推板的运动幅值来生成所需振幅的内孤立波,无法直接有效实现对内孤立波的数值模拟。而其余两种数值造波方法均能直接有效模拟内孤立波的生成与传播,且所得结果与理论及实验结果吻合较好。

表2为三种数值造波方法模拟效率的对比。从中可以看出: 在相同电脑配置(windows 7 64位,16G内存,i7处理器)、相同振幅、相同理论模型、相同模拟时间的情况下,平板拍击法与双推板法所用时间相当,速度入射边界法所用时间明显少于前面两种方法。其原因是仿物理造波方法使用动网格,需占用更多的计算资源。

表2　三种数值造波方法模拟效率的比较Tab.2　Efficiency of three numerical wave-generating methods

4　结论

利用 Fluent软件计算了基于 KdV、mKdV和eKdV三种理论模型与双推板、平板拍击和速度入口三种数值造波方法所组成的九种工况,并将模拟结果与理论及实验结果进行对比和分析,得到如下结论:

(1) 三种数值造波方法均能实现对内孤立波的数值模拟,除双推板法所得结果与理论及实验结果存在一定差异,需通过加大推板的运动振幅来实现所需的造波效果外,其余两种数值造波方法均能直接有效模拟内孤立波的生成与传播,且所得结果与理论及实验结果吻合较好。

(2) 双推板法与平板拍击法均属于仿物理造波法,但平板拍击法所得结果明显优于双推板法,因此实际使用时应优先选用平板拍击法。

(3) 同工况同配置条件下,速度入射边界法模拟效率最高,所得结果更准确,因此实际使用时应优先选用速度入射边界法。

王飞,2015. 内孤立波作用下小尺度竖直圆柱体的水动力特性研究. 青岛: 中国海洋大学博士学位论文,21—32

王旭,林忠义,尤云祥,2014. 两层流体中内孤立波质量源数值造波方法. 上海交通大学学报,48(6): 850—855

方欣华,杜涛,2005. 海洋内波基础和中国海内波. 青岛:中国海洋大学出版社,1—2

付东明,尤云祥,李巍,2009. 两层流体中内孤立波与潜体相互作用数值模拟. 海洋工程,27(3): 38—44

陈钰,朱良生,2009. 基于Fluent的海洋内孤立波数值水槽模拟. 海洋技术,28(4): 72—75,100

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

徐鑫哲,2012. 内波生成机理及二维内波数值水槽模型研究.哈尔滨: 哈尔滨工程大学硕士学位论文,35—36

黄文昊,尤云祥,王旭等,2013. 有限深两层流体中内孤立波造波实验及其理论模型. 物理学报,62(8): 084705

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

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

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

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

Wei G,You Y X,Su X B,2009. Two-dimension numerical internal wave tank for navier-stokes equation model in the stratified fluid. In: Proceedings of the Fifth International Conference on Fluid Mechanics. Berlin Heidelberg: Springer,364—367

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