海洋内波是指在海水稳定层化的海洋中产生的、最大振幅出现在海洋内部的波动。海洋内波的产生应具备两个条件:一是海水密度稳定分层, 二是要有扰动能源, 缺一不可。由于海水的密度分布经常处于不均匀状态, 因此海洋内波是一种比较普遍的现象(杜涛等, 2001; 柯自明等, 2009; Xu et al, 2012)。其中, 内孤立波是一种强非线性的潮成内波, 非线性与色散效应的平衡使得其在传播过程中可以保持稳定的波形和不变的速度(方欣华等, 2005)。内孤立波通常具有很大的振幅和波长, 在传播过程中携带有巨大的能量, 对海洋结构物构成很大的威胁(Roberts, 1975; Osborne et al, 1980; 张波等, 1998)。

现在广泛使用的内孤立波理论有KdV理论、mKdV理论、eKdV理论及MCC理论等(Michallet et al, 1998; Choi et al, 1996, 1999; Helfrich et al, 2006)。近些年来, 随着计算机技术及计算流体力学的发展, 采用计算流体力学(CFD)的方法研究内孤立波的生成传播及其与结构物的相互作用, 已成为重要的研究方法之一。内孤立波数值模拟常用的方法有两类--仿物理造波法和纯数值造波法。仿物理造波法包括平板拍击法、双推板法、摇板法等, 纯数值造波法包括速度入口法和源项造波法。付东明等(2009)结合KdV理论发展了双推板数值造波方法; 陈钰等(2009)基于KdV理论, 采用速度入口法, 建立了可有效模拟弱非线性海洋内孤立波的二维分层数值水槽; 高原雪等(2012)基于MCC理论, 采用速度入口法, 建立二维内孤立波数值水槽, 实现了振幅可控的内孤立波数值模拟; 李水娟等(2012)基于KdV理论, 采用速度入口法建立了三维内孤立波分层数值水槽; 王旭等(2014)结合eKdV理论, 从理想流体完全非线性欧拉方程出发, 发展了内孤立波质量源数值造波方法, 所得结果跟理论与实验吻合较好; 王伟等(2016)基于KdV、mKdV理论, 采用平板拍击法进行内孤立波数值模拟, 并详细分析了内孤立波波致流场变化。

自从Brorsen等(1987)基于“水下爆炸”原理, 首次提出基于势流理论的适合于边界积分方程法的质量源函数造波方法后, 质量源方法已被广泛应用于表面波数值模拟(Lin et al, 1999; 刘秀丽等, 2015; 周玲玲等, 2015), 但目前采用质量源法模拟内孤立波的论文很少。王旭等(2014)采用线源形式的质量源, 结合eKdV理论成功模拟了内孤立波。本文采用两个点源形式的质量源, 基于FLUENT软件, 从黏性不可压缩流体的Navier-Stokes方程出发, 结合KdV、eKdV理论建立有限水深两层流体中的内孤立波质量源数值模拟方法, 并将所得结果与理论和实验相对比, 验证方法的可行性。

1 数学模型 1.1 造波方法

流体力学中的源方法是在流体力学控制方程中加入特定的源项, 以此激发产生所需要的流场。源项是一个广义量, 它代表了那些不能包括到控制方程的非稳态项、对流项与扩散项中的所有其它各项之和, 包括质量源项、动量源项、能量源项和湍流源项。其中, 质量源方法是在连续性方程中加入附加质量源项, 通过对附加质量源项进行特定的设置, 以此来产生目标流场。Brorsen等(1987)对质量源方法进行了详细的数学物理推导, Kim等(2010)从雷诺输运定理出发, 详细推导了质量源表达式。

本文用质量源法建立二维内孤立波数值水槽, 如图 1所示。上层流体的深度为h₁, 密度为ρ₁, 下层流体的深度为h₂, 密度为ρ₂, 建立如图 1中所示的坐标系oxz, ox轴位于两层流体尚未扰动时的分界面上, oz轴垂直向上。在左边界处, 上下层流体中各放置一个点源形式的质量源, 上层中的源Ω₁的强度为S₁(x, z, t), 下层中的源Ω₂的强度为S₂(x, z, t)。源Ω₁不断释放质量, 源Ω₂不断吸收质量, 从而产生下凹型的内孤立波。

在质量源区域, 加入源项后, 流体的控制方程变为如下形式:

式中, S(x, z, t)为源强, μ为动力黏度系数。

假定质量源引起的流体质量的增减对应入射的孤立波波面的变化, 则对于源Ω₁、Ω₂有:

式中, C为内孤立波相速度, η(t)为内孤立波波面。

质量源为点源, 在质量源区域内均匀分布, 可认为质量源强度仅为时间t的函数, 即S_i(x, z, t)=S_i(t), 则质量源Ω₁、Ω₂的源强的表达式如下:

式中, A₁、A₂分别为源Ω₁、Ω₂的面积。

1.2 消波方法

采用海绵层阻尼消波法(韩朋等, 2009)进行消波。该方法的原理是, 通过在特定的计算区域的动量方程中添加阻尼源项来削弱传入该区域的波动。阻尼消波的效果与消波段的长度有关, 该长度应至少为入射波波长的一倍。消波区加入阻尼源项后的动量方程为:

其中, ν为运动黏度系数, φ(x)为阻尼消波系数, 是在阻尼段起点为零的单调递增函数。阻尼消波系数取线性分布的形式, 如下式所示:

式中: α为阻尼系数内的经验参数, x₁为海绵层的起始坐标, x₂为海绵层末端坐标; L_S为海绵层的长度, 在数值计算中通常取1-2倍波长。

2 内孤立波数值模拟 2.1 模型建立与网格划分

王飞(2015)在中国海洋大学物理海洋实验室, 利用密度分层水槽, 采用重力塌陷法制造不同振幅的内孤立波, 通过实验研究内孤立波的生成及流场特性。水槽长15m, 宽0.35m, 高0.7m, 额定水深0.6m。采用Oster双缸法制取密度分层流体, 上层流体高度9.5cm, 密度1035kg/m³, 下层流体高度48.5cm, 密度1054kg/m³; 采用PIV技术监测流场信息。为与实验结果进行对比, 本文选取上述流体参数, 以实验中得到的三个波高作为设定波高进行数值模拟。

黄文昊等(2013)通过实验, 研究了四种内孤立波理论的适用范围, 利用内孤立波的非线性参数ε和色散参数ξ来定量表征各内孤立波理论的适用区间, 结果表明: ξ < 0.1, ε≤ξ时, KdV理论适用; ξ < 0.1, ξ < ε < 时, eKdV理论适用。以此为依据, 对三组工况分别选用适用的内孤立波理论作为数值模拟输入条件。数值模拟工况如表 1所示:

网格划分:数值水槽采用结构化网格进行划分, x方向网格尺寸取L_X=2cm, z方向网格尺寸取L_Z=0.5cm, 共87000个网格。消波区的长度取为1-2倍波长。

源的设置:两质量源的大小应满足点源假设, 即其高度H_S应远小于数值水槽的水深H, 宽度L_Ω应远小于内孤立波的波长λ, 通常质量源宽度应小于波长的5% (Lin et al, 1999; Hafsia et al, 2009), 故将源Ω₁、Ω₂分别放置在左边界处上、下层流体中的一个网格内。Lin等(1999)指出质量源位置距离分界面1/3层厚处模拟结果与理论最相符, 并将此位置称为“standard elevation”。故参考Lin等(1999)的结论, 将两源距离两层流体分界面的距离各设为相应层厚度的1/3, 即Ω₁距离两层流体分界面3.5cm, Ω₂距离分界面16cm。从后面的模拟结果来看, 这样的设置取得了良好的效果。

2.2 FLUENT设置

边界条件:水槽左边界设为对称边界; 水表面采用刚盖假设(方欣华, 2005), 设为固壁边界; 水底设为固壁边界; 右边界设置为压力出口。

液面追踪方法:采用VOF方法(韩朋等, 2009)来追踪两层流体分界面处的内孤立波波面。流体体积函数F定义为单元内流体所占有的体积与该单元可容纳流体体积之比。若单元体被指定相流体占满, 则该单元F值为1;若单元体被另一相流体占满, 则单元的F值为0;若单元体的F值在0与1之间, 则表示该单元为包含两相流体的交界面单元。

UDF设置:利用FLUENT中的DEFINE_ SOURCE宏定义上下两源的源强, 添加到FLUENT中造波源区域的源项加载区, 嵌入到造波区控制方程中; 利用DEFINE_PROFILE宏在压力出口边界给定计算域静水压强, 使得出口边界处水体能够自由流出, 保证计算域内压力出口附近上下层水深为定值(张齐焰, 2008; 陈钰等, 2009)。

求解设置:控制方程的对流项采用一阶迎风格式, 扩散项采用中心差分格式, 压力速度耦合方式采用PISO算法, 压力插值选用体积力加权法, 两相界面的构造方法选用几何重构; 初始时水槽中流场速度取为0, 时间步长设为0.01s, 每个时间步最大迭代次数20次。

3 数值造波结果及分析 3.1 质量源法中内孤立波的生成

两质量源的强度随模拟时间不断发生变化, 使得造波源附近的流场速度发生变化。图 2为工况A中两造波源附近两个时刻的速度矢量图, 图中白线分别为数值水槽左边界和上边界, 黑色矩形分别为质量源Ω₁和Ω₂。从速度矢量图中可以看出, 随着源Ω₁不断释放质量、源Ω₂不断吸收质量, 源附近上下层流体的速度场均发生变化, 上层流体产生沿x正向的水平速度, 下层流体产生沿x负向的水平速度, 这与内孤立波波致流场中的速度趋势是相同的。

图 3为A工况造波初期四个不同时刻质量源区域附近两相界面图, 展示了造波开始后左侧边界附近两相界面的变化, 其中灰色部分为上层流体, 黑色部分为下层流体。从图中可以看出, 两层流体的界面在源Ω₁不断释放质量、源Ω₂不断吸收质量的驱动下, 缓慢而平稳地发生下凹, 并最终形成一个完整的下凹形内孤立波。

3.2 数值造波结果及对比

图 4为三个工况下内孤立波数值模拟所得波形(CFD)与理论波形(KDV)和实验波形的对比。从图中可以看出, 三个工况下数值模拟波形均光滑对称, 且与理论及实验波形吻合良好; 各工况数值模拟所得波高分别为-3.988cm、-7.320cm和-8.092cm, 与设定波高吻合良好。从对比中可以看出, 本文所述质量源方法模拟内孤立波所得波形和波高同理论解和实验结果是一致的。

图 5为波谷断面处波致水平流速沿垂向分布数值模拟、理论解及实验结果的对比图。从图中可以看出, 工况A中数值模拟与理论解及实验均吻合良好。工况B、C中, 上层流体中波致水平流速较大, 致使实验监测数据失效, 但在该层数值模拟与理论解吻合较好; 在下层流体中, 数值模拟与理论解及实验均吻合较好。

图 6为内孤立波波致水平流速数值模拟、理论解及实验的对比图。从对比中可以看出, 工况A中模拟峰值比理论和实验结果略小, 图像介于理论解和实验之间, 整体上与二者吻合良好; 工况B、C中同样由于上层流体的流速较大, 超出系统的测量范围, 造成实验监测数据失真, 仍可以看出, 上层流体中数值模拟与理论解吻合良好, 下层流体中三者的吻合度均良好。

4 结论

采用两个点源形式的质量源, 分别放置于两层流体中的上下层中, 上层流体中的质量源释放流体, 下层流体中的质量源吸收流体, 激发产生下凹型内孤立波。推导质量源强度表达式, 并利用商业软件FLUENT, 基于其UDF编译功能, 以内孤立波KdV和eKdV理论为依据, 从不可压缩流体的Navier-Stokes方程出发, 发展了一种内孤立波数值造波方法, 并将数值模拟结果与理论及实验进行对比。结果表明:

(1)造波过程中, 在上层流体中的源释放质量、下层流体中的源吸收质量的驱使下, 上下层流体速度呈现出相反的情形, 使得两层流体分界面逐渐下凹, 最终形成一个不断向前传播的下凹形内孤立波。

(2)三个工况下波形、波高及波致水平速度的数值模拟结果均与理论及实验的吻合较好, 验证了此方法的有效性, 表明本文所述质量源法能够有效模拟内孤立波的生成与传播。

(3)该方法作为一种纯数值造波方法, 避免了动网格的使用, 计算效率高, 耗时短, 另外不需要已知速度等边界条件, 可用于研究速度边界未知的情况, 如波流相互作用等。

(4)该方法同样可以模拟上凸型内孤立波。对于上凸型内孤立波, 其波面η(t)始终为正, 则(5)、(6)两式中的源强表达式S₁(t)始终为负, S₂(t)始终为正, 表示上层流体中的源吸收质量, 下层流体中的源释放质量, 从而驱动两层流体界面发生上凸, 形成上凸型内孤立波。

实际海洋中环境复杂, 由于表面波、流以及海底地形等因素的影响, 内孤立波波形通常不如理论波形和数值模拟波形那样光滑对称(Xu et al, 2012), 在以后的内孤立波数值模拟中可以考虑波流相互耦合以及复杂海底地形下的内孤立波的数值模拟, 这样得到的波形将更趋近于实际海洋中的内孤立波的波形。