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

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

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

1.1 KdV理论

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

其中:

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

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

1.2 eKdV理论

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

其中:

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

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

1.3 mKdV理论

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

其中:

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

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

2 数值造波方法

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

2.1 双推板法

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

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

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

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

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

实际设置时, 应使d₁略大于h₁+|η₀|, 以防上下层流体掺混(付东明等, 2009)。

2.2 平板拍击法

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

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

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

2.3 速度入射边界法

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

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

3 数值模拟与结果分析 3.1 内孤立波物理实验

王飞(2015)在中国海洋大学物理海洋实验室分层流体水槽中进行了内孤立波物理实验。水槽长15m, 宽0.35m, 高0.7m, 额定水深0.6m。采用Oster双缸法制取密度分层流体, 上层流体高度9.5cm, 密度1035kg/m³, 下层流体高度48.5cm, 密度1054kg/m³。采用重力塌陷法制造不同振幅的内孤立波, 使用两个斜面在水槽末端进行消波, 利用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的比例逐渐变大。

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

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

3.3 计算结果与分析

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

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

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

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

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

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

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

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

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

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

4 结论

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

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

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

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