内波是发生在密度稳定分层的流体内部的一种波动, 而内孤立波是一种特殊的内波。内孤立波频繁地发生在海洋中, 在世界海洋的许多大陆边缘斜坡上都曾观测到内孤立波痕迹(Filonov et al, 2000; Azevedo et al, 2006; Vázquez et al, 2008)。内孤立波具有强垂向流和水平流, 受大陆或近岸斜坡影响的内孤立波会浅化、破碎, 并造成强烈的湍流混合(Forgia et al, 2018), 进而影响温度(Leichter et al, 1996; Davis et al, 2011; Walter et al, 2012)、氧气和营养物质的空间分布(Leichter et al, 1996; Omand et al, 2011; Walter et al, 2014)。

内孤立波对于结构物具有潜在的破坏性。在非线性效应和频散效应的平衡作用下, 在传播过程中, 内孤立波的波形和速度可以保持数百公里不变(Grue et al, 2000; Grimshaw et al, 2002)。在中国南海(South China Sea, SCS)曾记录到振幅最大为170 m的内孤立波, 且上、下层之间的水平速度差超过3.4 m/s (Warn-Varnas et al, 2010)。正是由于内孤立波具有振幅尺度大、能量集中的特点, 其引起流场和密度场的迅速变化可能对海洋工程的结构物以及水下潜体造成严重威胁(邹丽等, 2020)。由此可见, 内孤立波研究在海洋工程、海洋观测等方面具有十分重要的应用意义。因此研究不同造波条件下生成的内孤立波运动的流场特征具有重要的学术意义和实际应用价值。

在Chen等(2007)的实验室实验中, 内孤立波是在水闸两侧流体界面高度不同的坍塌机制下生成的(由水闸划分为左侧环境区和右侧扰动区)。目前在内孤立波的实验室研究过程中, 在分层水体制取方面, 要配制出质量高的稳定的分层水体, 需保证流入水槽流体运动的平稳性及均匀性, 使得配制过程要耗用很长的时间, 而且试验成功率低, 需经多次重新配制才能达到要求(徐鑫哲, 2012)。此外, 由实验仪器收集到的信息也往往比较有限, 即使使用了高频相机和超声波测试仪, 依然难以得出流场内部精确的密度、涡度的时空变化。

数值模拟是一种可以获得精确的流场内部变化信息的有效方法。对内孤立波的数值模拟研究基本都是从数值求解Navier-Stokes (N-S)方程入手, 根据求解N-S方程的方法可进一步分为: (1)直接求解N-S方程(直接数值模拟方法, direct numerical simulation, DNS); (2)求解雷诺平均的N-S方程(reynolds-averaged Navier-Stokes simulation, RANS) (Vlasenko et al, 2002; Lin et al, 2012); (3)求解空间滤波的N-S方程(大涡数值模拟法, large-eddy simulation, LES) (沈治等, 2009)。然而, RANS的模拟方法只能提供湍流的平均信息, 在模拟结果中, 上凸型内孤立波基本不呈现涡旋结构, 这与实验室实验结果不符。LES的模拟方法适用于研究运动尺度较大的湍流。目前有关直接数值模拟方法的内孤立波研究相对较少, 理论上来说, 直接数值模拟方法能够更好地模拟出内孤立波发展和传播过程中的细微结构, 如因断裂形成的湍流涡旋结构。

本文参考Chen等(2007)重力塌陷激发内孤立波的水槽实验, 构建了重力塌陷激发内孤立波的数学模型。根据Druzhinin等(2015)的密度跃迁公式计算初始密度场。以Boussinesq假定下的N-S方程组及盐分对流扩散方程作为控制方程, 并采用有限差分法进行离散, 实现了重力塌陷造波方式下形成的内孤立波运动的数值模拟。通过与Chen等(2007)的实验结果对比, 验证了本文方法的有效性, 而后本文进一步与Lin等(2012)对比不同数值模拟方法下呈现的涡旋结构。之后, 通过改变造波参数, 即流体密度差、上下层水深比和台阶深度(水闸两侧初始密度界面的高度差), 对一系列算例进行数值模拟, 得到流场的密度界面、涡旋结构、涡度等的时空演化结果, 比较不同造波条件对内波后续发展过程中的波型、涡度、水平速度、流场动能和特征波长的影响。

1 数学模型

本文参考Chen等(2007)重力崩塌激发内孤立波的水槽实验, 构建了含有斜坡、重流体、轻流体等的数学模型, 如图 1所示。模型设置为: 在H × L的水箱中(H为水箱高度、L为水箱长度), 设置左右两侧不同初始密度界面高度的两层流体, 其中上层流体密度ρ₁, 下层流体密度ρ₂, 左侧环境区域的上层流体水深H₁, 右侧扰动区域的上层流体水深h₁, H₂、h₂为对应的下层流体水深, 其中H₁+H₂=h₁+h₂=H, H₂=h₂-_η0, _η0为台阶深度。扰动区长度L_p=0.3 m, 斜坡长度L_s=1.882 m。L_w为波长, a为波高。两种流体的运动黏度均为v。计算时的各算例参数设置参见表 1, 时间步长Δt=0.005Δx, 盐分施密特数(salt Schmidt number) Sc_s为100。

1.1 控制方程

本文以不可压缩N-S方程组(连续性方程及动量守恒方程)和对流扩散方程为控制方程, 三个方程分别为:

其中, 为梯度算子; 为流体速度; 为时间; 为流体密度; 为压力; 为流体的运动黏度; e_z为重力方向无量纲单位矢量, e_z=(0, 0, -1); 为盐度; 为盐分的有效扩散系数; 为水的密度; 为盐分的存在造成的流体密度增加系数; 控制方程中变量的上标"^"符号表示有量纲的物理量。根据Boussinesq假定, 将控制方程进行无量纲化(He et al, 2018)可得(具体可参见表 2):

其中, 为无量纲梯度算子; u为无量纲流体速度; t为无量纲时间; p为无量纲压力; s为无量纲盐度; 为特征运动黏度; 为特征长度(验证算例取), 特征长度为0.2 m; 为特征速度, 除工况6和工况7外, 特征速度均为0.2465 m/s, 工况6的特征速度为0.1468 m/s, 工况7的特征速度为0.2029 m/s; 为特征时间。

控制方程中的盐分对流扩散方程用迎风紧致差分格式(upwinding combined compact difference, UCCD)离散(Yu et al, 2019)。N-S方程, 利用二阶显式Adams-Bashforth法离散时间项, 二阶迎风格式离散对流项, 以二阶中心差分格式离散黏性扩散项, 压力项通过超松弛迭代法(successive over relaxation method, SOR)求解(Yu et al, 2019)。

1.2 初始及边界条件

本文采用了密跃层假设(Druzhinin et al, 2015), 认为流体间的密度是按双曲正切函数变化。相对于简单的线性分层, 双曲正切函数具有阶跃突变的特点, 这与实际海洋中的密跃情况相接近。本文设置计算流体的初始盐度分布如下:

式中, s(t=0)为初始盐度; ϕ为网格点到密跃层中间界面ψ的符号距离函数(王德军等, 2003); Ω₁、Ω₂分别为上层流体和下层流体; x=(x, y)为网格点; d为任意网格点到密跃层中间界面的最短距离, 最短距离沿法线方向; ε为密跃层厚度。初始密度可通过式(12)计算:

式中, ρ(t=0)为初始密度, 为特征密度, s(t=0)为初始盐度。系统由静止释放, 初始速度为0。

对于速度边界条件, 在顶底及两侧界面处均使用了自由滑移边界条件(free-slip boundary condition)。以顶部为例:

对于盐度边界条件, 在界面处均使用了无通量条件。

采用浸入边界法, 求解斜坡边界处的速度及盐分浓度变化。通过在动量方程中添加虚拟力F, 用以反映不规则地形对流体的作用力:

其中, u_s为固体速度, u_s=0, u_l为流体速度, η为固体在网格单元中的体积分数, Δt为模拟时的迭代时间间隔。

2 模型检验

为检验本文数值方法的有效性, 将本文模拟结果和前人的实验(Chen et al, 2007)及模拟结果(Lin et al, 2012)进行比较。

本文与Chen等(2007)实验结果进行了两方面的比对, 分别为界面波形的演化过程和定点波动时间序列。图 2为界面波形的演化过程, 图 2a为Chen等(2007)实验的结果, 图 2b为本文的模拟结果。由图 2可以看出, 实验与模拟方式之间最明显的差异是: 在实验水槽内填充阶段引入了厚度较小的水闸, 以便在左右两侧流体区域之间产生不同的分层。实验中水闸的去除会引起垂直剪切应力, 并将力传递到与其直接接触的流体。但考虑到水闸的厚度较小, 其所产生的干扰应是相当有限的。图 2证实, 在去除水闸不久后, 扰动区附近模拟得到的密度界面便与通过实验观察到的染料的空间分布非常吻合: 在无量纲时间t=3.5, 最初放置在扰动区内的较重的流体开始呈现出内孤立波的典型形状; 在t=6.22至t=10.5, 内孤立波向左传播, 同时保留其形状; 在t=6.22至t=15, 在实验和模拟中, 均存在界面尾随内波。

总体而言, 本文的直接数值模拟结果与Chen等(2007)实验结果所得的流场界面波形具有较好的一致性。

图 3为定点波动时间序列对比图, 参照实验选取定点为(8.7 m, 0.3 m)。结果表明, 本文模拟的定点波动较实验结果波动幅度更大。但总体而言, 本文模拟结果的定点波动与实验结果较为吻合, 能够基本体现定点的波动特征。图 3从微元层面反映了稳定传播的下凹型内孤立波(h₁>H₁, H₁＜H₂)。

本文与Lin等(2012)的RANS模拟结果进行了界面涡旋结构对比。如图 4所示, 在扰动区水体刚开始下落时, 本文与Lin等(2012)模拟结果中的密度界面位移变化基本一致。随着时间的推移, 两种模拟方法开始有所区别: 在Lin等(2012)模拟结果中, 上凸型内波基本不呈现涡旋结构, 这与实验室实验结果不符; 本文模拟结果则呈现了因断裂形成湍流涡旋结构的详细过程, 这为后续对比内孤立波的起涡能力提供了充分条件。

综上, 对验证算例的模拟结果表明, 本文所采用的直接模拟数值模型有效地模拟了重力塌陷式内孤立波的波动变化过程, 相比于Lin等(2012)的RANS模拟结果呈现出更为明显的涡旋结构。

3 模拟结果分析

为研究造波条件对重力塌陷式内孤立波发展的影响, 本文设计了8组算例(工况1—工况8), 各算例的主要参数参见表 1。由于本文主要考虑的是造波条件对内孤立波后续发展过程中的波型、涡度、振幅和水平速度的影响, 而根据验证算例和Chen等(2007)实验可知, 如果斜坡与扰动区距离较远, 在内孤立波没有传递到斜坡之前, 斜坡的影响基本可以忽略, 因此本文没有在这些算例中加入斜坡。

工况1、工况2、工况3的H₁相同, η₀不同, 可研究台阶深度η₀对内孤立波发展的影响。工况1、工况6、工况7可研究流体密度差对内孤立波发展的影响。工况1、工况4、工况5、工况8可研究环境区上下层水深比(H₁/H₂)对内孤立波发展的影响, 四个算例的η₀相同, 且均设置为H₁>h₁, 即初始波型为上凸型。其中工况1、工况4、工况5设置为H₁>H₂; 工况8设置为H₁＜H₂, 可研究初始波型为上凸型内孤立波的反转现象。

3.1 流场涡度

涡度, 由于本文为二维的数值模拟, 因此涡度为

其中, k为z方向单位长度分量。

在研究过程中发现流体密度差对无量纲速度几乎没有影响, 因此本节将不给出工况6、工况7流场涡度对比图。图 5a、5b、5c为不同台阶深度的流场涡度时间演化的对比图。结果表明, 台阶深度η₀对流场涡度的影响较为显著。在无量纲时间t=2.5时, 工况3几乎没有出现正涡度, 在之后的时间t=3.5至t=5.5, 工况1、工况2均表现出正涡度先增大后衰减的过程, 工况3的衰减过程则相对不明显。这表明台阶深度越大, 对尾部的湍流振荡影响越显著, 湍流振荡越剧烈, 且振荡持续时间相对越长。除此之外, 台阶深度越大, 产生的内孤立波波高越大, 内孤立波传播的平均速度越大。

图 5a、5d、5e为不同环境区上下层水深比(H₁/H₂)的流场涡度时间演化的对比图。结果表明, 上下层水深比(H₁/H₂)对流场涡度的影响较小。在t=2.5时, 3个算例均开始出现正涡度, 在之后的时间t=3.5至t=5.5, 工况5的负涡度衰减的更快, 且正涡度相对工况1、工况4更明显。这表明在H₁>h₁的条件下, 上下层水深比(H₁/H₂)越小, 前部波衰减的相对更快, 尾部的湍流振荡相对更剧烈。

图 5a、5f为不同水层条件下(分别为稳定和不稳定的内孤立波)的流场涡度时间演化的对比图。结果表明, 工况8的负涡度衰减的更快, 且正涡度相对工况1数值更大。这定性地表明了不稳定内孤立波在向前传播的过程中, 存在能量的向后部传递。

3.2 最大水平速度

图 6为不同算例的流场最大水平速度随水平位置的变化, 两者均为有量纲形式(由于流体密度差对无量纲速度U_max影响不大, 但密度差会影响到特征速度, 考虑到图形重叠问题, 这里采用有量纲形式), 其中为不同时间不同垂向深度在同一水平位置处所能达到的最大水平速度分量负值(在波传播方向上向前)。结果表明:

(1) 增大流体密度差, 流场最大水平速度将随之增大;

(2) 随着台阶深度的增大, 流场最大水平速度将随之增大; 台阶深度较大的算例(工况3)初始会有一段最大水平速度保持不变的阶段, 这可能是因为增大台阶深度会加长初始水体顺时针涡流运动的持续时间;

(3) 在水层条件h₁＜H₁, H₁>H₂下, 产生稳定的上凸型内孤立波, 随着上下层水深比的增大, 流场最大水平速度将随之增大;

(4) 在水层条件h₁＜H₁, H₁＜H₂下, 产生不稳定的内孤立波, 在m位置处, 不稳定内孤立波的最大水平速度发生突然下降, 后继续恢复平稳。

由于不稳定内孤立波的前部上凸型内孤立波处为顺时针涡流, 后部下凹型内孤立波处为逆时针涡流。这表明在位置处, 下凹型内孤立波的上层流体负速度的绝对值超过了前部上凸型内孤立波的下层流体负速度的绝对值, 证实了不稳定内孤立波在向前传播的过程中, 存在能量的向后部传递。

3.3 最大波高及振幅

图 7为不同算例的最大波高H_max随时间t的变化, 两者均为无量纲形式, 其中H_max为内孤立波运动状态趋于稳定时, 流场内的最大波高。由图可知, 最大波高的变化过程可大致分为三个阶段: 第一阶段为快速衰减阶段, 在这一阶段, 由于两侧水体的密度交互作用, 水体倾覆, 最大波高快速衰减; 第二阶段为振荡阶段, 在这一阶段, 水体倾覆运动基本完成, 密度界面形成稳定的前部高程波, 但由于尾随界面振荡的存在, 表现出最大波高下降再升高的反复现象, 但总体趋势是呈相对快速衰减; 第三阶段为缓慢衰减阶段, 达到一定时间后, 前部波以较为稳定的状态向左传播, 由于能量耗散, 内孤立波波高逐渐缓慢减小。对比各算例结果可知:

(1) 流体密度差对最大波高基本没有影响;

(2) 随着台阶深度的增大, 最大波高第二阶段的反复现象越明显, 即振荡越剧烈;

(3) 在台阶深度一定, 且产生稳定内孤立波的水层条件下(h₁＜H₁, H₁>H₂), 内孤立波表现出极为相似的波高变化。

根据上述的内孤立波传播特点, 本文选取无量纲时间t=9—10内的界面波动振幅的平均值作为内孤立波的振幅A_w。计算公式如下:

图 8为内孤立波的振幅A_w随各项造波参数的变化曲线。结果表明: 台阶深度对内孤立波的振幅影响显著, 上下层水深比对振幅影响相对较小, 流体密度差对振幅基本不产生影响。

3.4 波形反转

图 9为工况1流场界面内孤立波的演化过程。稳定的内孤立波在传播过程中, 其前部内孤立波处的顺时针涡流始终维持较大速度, 这表明能量集中在前部波内。图 10为工况8流场界面内孤立波波型反转的演化过程。与图 9对比, 可清楚观察到工况8的内孤立波在传播过程中发生波型不稳定转换的现象:

(1) 在H₁>h₁的条件下, 工况8内孤立波的初始波型表现为上凸型。随着时间的递进, 前部的上凸型内孤立波振幅逐渐减小, 反之, 其后部的下凹型内孤立波振幅逐渐增大;

(2) 不稳定的内孤立波在传播过程中, 前部上凸型内孤立波处的顺时针涡流速度逐渐减小; 反之, 后部的下凹型内孤立波处的逆时针涡流速度逐渐增大。这证实了不稳定内孤立波在传播过程中会发生内部能量传递, 能量由前部上凸型内孤立波, 逐渐传递到后部下凹型内孤立波。

以上均是定性判断了内孤立波运动过程中的能量变化, 下面将通过给出内孤立波动能的定量计算, 分析内孤立波运动过程中的能量变化。

图 11为稳定内孤立波(工况1)与不稳定内孤立波(工况8)的流场动能水平分布对比图。其中, E_kr为不同水平位置处的流体动能占整个流场动能的比例, 计算公式如下:

其中, ρ为微元的密度, u_Ω为微元的速度, x₀为微元的水平位置, dV为微元的体积, ds为微元的剖面积(二维时可认为是高度dy), Ω为计算域。

结果表明, 在无量纲时间t=12.5至t=16的这段时间内, 稳定内孤立波(工况1)的流场动能水平分布形态基本不发生改变, 第一个转折点处判断为前部波波峰位置; 不稳定内孤立波(工况8)的流场动能水平分布形态发生明显改变, 对比t=12.5与t=16的第一个转折点处的动能占比, 可看到随时间的推移, 前部流体动能占比减小。这表明在不稳定内孤立波在向前传播的过程中, 存在能量的向后部传递。除此之外, 工况8的后部动能占比明显大于工况1, 这可能是由于工况8的上下层水深比接近1, 不稳定内孤立波在主波后方产生了连续的内波波列, 造成能量分散。

4 结论

本文通过建立重力塌陷激发内孤立波运动的数值模型, 采用直接模拟数值方法, 研究了不同造波条件对内孤立波发展的影响, 结果表明:

(1) 通过直接数值模拟方法得到的流场界面呈现更为明显的涡旋结构。此外, 从定性和定量的角度, 证实了不稳定内孤立波传播过程中存在能量的向后传递。

(2) 增大台阶深度、流体密度差、上下层水深比, 流场最大水平速度都将随之增大, 其中台阶深度的影响最为显著。

(3) 对于相同的台阶深度, 初始涡流保持相同, 但是随着上下层水深比的减小, 其强度明显下降; 台阶深度对初始涡流的垂直结构的影响要大于上下层水深比。

(4) 台阶深度对内孤立波的振幅影响显著, 上下层水深比对振幅影响相对较小, 流体密度差对振幅基本不产生影响。