内孤立波是一种最大振幅发生在密度稳定层化海洋内部的波动，通常由洋流通过不平坦的海底而产生，是海洋中频繁发生的现象^[1-2]。我国南海是内孤立波的频发海域, 被评价为研究大幅内波的“天然实验室”，Huang等^[3]在南海观测到波幅240 m，流速2.55 m/s的大幅内孤立波。大幅内孤立波在传播中携带着巨大能量，对海洋结构物存在严重威胁。因此，内孤立波的研究对于南海海洋环境的认识和开发具有极大的现实意义。在内波问题的研究上，当密度跃层的厚度较小时，可以简化为两层流体间的内波问题，即每一层流体的密度均为常数，两层流体界面处存在密度跳跃。在两层内孤立波的水池实验研究上，文献[4-7]采用重力塌陷法对两层流体间内孤立波的波形、波速以及速度场等特性进行研究。黄文昊等^[8]采用双推板法生成内孤立波，对实验中内孤立波的振幅实现了自主控制。国际上流行的内孤立波理论模型主要有Korteweg-deVries (KdV)模型、MCC模型以及High-Level Green-Naghdi (HLGN)模型等。KdV模型是当前应用最广泛的内波模型，但是对大幅内孤立波进行研究时，KdV模型给出的结果往往不尽如人意^[4]。

考虑到大幅内孤立波的强非线性，文献[9-10]假设特征波长较总水深更长，振幅与总水深相当，从完全非线性欧拉方程出发，采用层平均水平速度方法，建立了适用于刚盖近似下大幅内孤立波的MCC-RL模型。在MCC-RL模型中，由于没有引入小幅波假定，同时其形式简单，因此被广泛地应用到浅水模式下的大幅内孤立波问题的研究中，并展示了良好的性能^{[7-8, 10-11]}。

当上、下层流体密度相差较大时，自由面对内孤立波特性的影响不可忽略，引入刚盖近似假设是不合理的。文献[5, 12]发现考虑自由面效应，MCC-FS模型得到的波形结果相对MCC-RL模型要更窄，同样更接近实验值。然而，Kodaira等^[5]求解MCC-FS模型时，在MCC-FS模型的4个方程中消去上、下层平均速度，推导得到只保留自由面以及内波界面2个未知数的2个非线性微分方程，求解出自由面和内波波面之后，可计算上、下层的层平均速度。本文中则是对MCC-FS模型4个方程进行直接求解。另外，MCC-FS模型的速度场结果从未被展示以及验证。

本文通过与实验值以及他人数值结果的对比，对本文算法以及速度场结果进行了验证。推导出MCC-FS模型的算法以及基于MCC-FS模型速度场结果，并研究了不同密度比对浅水模式下内孤立波特性的影响。

1 保留自由面的MCC模型

考虑自由面的两层流体间内孤立波传播示意图如图 1所示。假定两层流体均无粘、不可压缩。在未受扰动时，上层流体厚度为h₁，密度为ρ₁；下层流体厚度为h₂，密度为ρ₂。建立大地直角坐标系oxz，其中，ox轴与未扰动的内界面重合，指向内孤立波传播方向。η₁、η₂分别为上层流体与下层流体的厚度。

保留自由面的MCC-FS模型方程为^[5]：

$ {\eta _{1,t}} + {({\eta _1}{\bar u_1})_x} = 0 $

(1)

$ {\eta _{2,t}} + {({\eta _2}{\bar u_2})_x} = 0 $

(2)

$ \begin{array}{c} {{\bar u}_{1,t}} + {{\bar u}_1}{{\bar u}_{1,x}} + {\rm{ }}g{({\eta _1} + {\eta _2})_x} = \frac{1}{{{\eta _1}}}{\left( {\frac{1}{3}\eta _1^3{G_1}} \right)_x} - \\ \frac{1}{{{\eta _1}}}{\left( {\frac{1}{2}\eta _1^2D_1^2{\eta _2}} \right)_x} + \left( {\frac{1}{2}{\eta _1}{G_1} - D_1^2{\eta _2}} \right){\eta _{2,x}} \end{array} $

(3)

$ \begin{gathered} \bar{u}_{2, t}+\bar{u}_{2} \bar{u}_{2, x}+g\left(\frac{\rho_{1}}{\rho_{2}} \eta_{1}+\eta_{2}\right)_{x}= \\ \frac{1}{\eta_{2}}\left(\frac{1}{3} \eta_{2}^{3} G_{2}\right)_{x}+\frac{\rho_{1}}{\rho_{2}}\left(\frac{1}{2} \eta_{1}^{2} G_{1}-\eta_{1} D_{1}^{2} \eta_{2}\right) \end{gathered} $

(4)

式中：ζ₁表示自由面位移；ζ₂表示内波界面位移。

$ {{\eta _1} = {\rm{ }}{h_1} + {\zeta _1} - {\zeta _2}} $

(5)

$ {{\eta _2} = {\rm{ }}{h_2} + {\zeta _2}} $

(6)

其中：u₁和u₂分别是上层流体和下层流体的层平均水平速度为：

$ {{{\bar u}_1}(x,t) = \frac{1}{{{\eta _1}}}\int_{{\zeta _2}}^{{h_1} + {\zeta _1}} {{u_1}(x,z,t){\rm{d}}z} } $

(7)

$ {{{\bar u}_2}(x,t) = \frac{1}{{{\eta _2}}}\int_{ - {h_2}}^{{\zeta _2}} {{u_2}(x,z,t){\rm{d}}z} } $

(8)

G_i作为一个运算规则，其定义为：

$ \begin{array}{*{20}{l}} {{G_i}(x,t) = - \frac{{({D_i}^2{\eta _i})}}{{{\eta _i}}}} \end{array} $

(9)

其中D_i≡∂_t+u_i∂_x。

通过消去速度项，得到了简化后的含二阶导数项的非线性微分方程^[5]：

$ {\alpha _{j1}}{\eta^{''}{}_1} + {\alpha _{j2}}{\eta^{''}{}_2} + {\alpha _{j3}}\eta _{12}^\prime + {\alpha _{j4}}\eta _{22}^\prime + {\alpha _{j5}}\eta _1^\prime \eta _2^\prime = {\rm{ }}{\alpha _{j6}} $

(10)

$ \bar{u}_{i}=c\left(1-\frac{h_{i}}{\eta_{i}}\right) $

(11)

式中：i, j=1, 2；式子中的上标撇号代表对水平方向上的X的求导，需给出波速c求解上式。该方程组中包含η₁、η₂2个未知数，封闭可解。

更多关于MCC-FS模型原始算法的内容参考文献[5]。

通过式(11)得到层平均水平速度，然而，MCC-FS模型方程是可以直接求解的，在方程(1)~(4)中，包含η₁、η₂、u₁、u₂等4个未知数，因此，方程组是封闭的，可以直接数值求解，并且可以直接求解出自由面、内波界面以及上、下2层平均水平速度，为之后计算MCC-FS模型的速度场提供了便利。图 2中，具体的给出了2种算法的不同之处。

2 MCC-FS模型的数值算法

在MCC-FS模型的本文算法中，将大地坐标系转换到随着内波面相速度平动的稳态坐标系oxz：

$ \eta_{i, t}=-c \eta_{i, X}, \bar{u}_{i, t}=-c \bar{u}_{i, X}, i=1,2 $

(12)

式(12)代入MCC-FS模型的方程(1)~(4)，可得到稳态坐标系下的MCC-FS模型方程。采用差商代替导数的方法进行空间离散，在这里以内波界面位移ζ₂为例，其他各量的表示方法类似，各阶导数可表示为：

$ \zeta_{2, X}(X)=\frac{\zeta_{2}(X+\Delta X)-\zeta_{2}(X-\Delta X)}{2 \Delta X} $

(13)

$ \begin{aligned} \zeta_{2, X X}(X)=&\frac{1}{\Delta X^{2}}\left[\zeta_{2}(X+\varDelta X)-2 \zeta_{2}(X)+\right. \\ &\left.\zeta_{2}(X-\Delta X)\right] \end{aligned} $

(14)

$ \begin{aligned} \zeta_{2, X X X}(X)=&\frac{1}{2 \Delta X^{3}}\left[2 \zeta_{2}(X-\Delta X)-\zeta_{2}(X-2 \Delta X)-\right. \\ &\left.2 \zeta_{2}(X+\Delta X)+\zeta_{2}(X+2 \Delta X)\right] \end{aligned} $

(15)

将内孤立波的波峰设在X=0处，无论是内波界面位移、自由面位移还是上、下层平均速度均满足对称性。同时，在无穷远处，内波界面位移、自由面位移以及上、下层平均速度均为0，同样以内波界面位移ζ₂为例，边界条件可以表示为：

$ \zeta_{2}(-\Delta X)=\zeta_{2}(\Delta X) $

(16)

$ \zeta_{2}(-2 \Delta X)=\zeta_{2}(2 \Delta X) $

(17)

$ \zeta_{2}(X)=0(X \rightarrow \infty) $

(18)

数值求解采用牛顿迭代法，当相邻迭代步得到变量的绝对值之差小于10^-7时，迭代停止。

Camassa等^[11]给出了MCC-RL模型中计算流体质点水平速度分布的近似方法，该方法可以应用于MCC-FS模型上。

1) 对于上层流体：

$ \begin{aligned} u_{1}(X, z)=& c\left[1-\frac{h_{1}}{\eta_{1}}+\left(\frac{\eta_{1}^{2}}{6}-\frac{\left(h_{1}-z\right)^{2}}{2}\right)\cdot\right.\\ &\left.\left(\frac{h_{1} \eta_{1, X X}}{\eta_{1}^{2}}-\frac{2 h_{1}\left(\eta_{1, X}\right)^{2}}{\eta_{1}^{3}}\right)\right] \end{aligned} $

(19)

2) 对于下层流体：

$ \begin{aligned} u_{2}(X, z)=& c\left[1-\frac{h_{2}}{\eta_{2}}+\left(\frac{\eta_{2}^{2}}{6}-\frac{\left(h_{2}-z\right)^{2}}{2}\right)\cdot\right.\\ &\left.\left(\frac{h_{2} \eta_{2, X X}}{\eta_{2}^{2}}-\frac{2 h_{2}\left(\eta_{2, X}\right)^{2}}{\eta_{2}^{3}}\right)\right] \end{aligned} $

(20)

根据式(19)、(20)以及MCC-FS模型得到的自由面以及内波波面结果，可以得到MCC-FS模型的水平速度分布。在垂向速度的计算上，Choi等^[10]给出了通过层平均水平速度得到垂向速度为：

$ w_{i}(X, z)=-\bar{u}_{i, X}\left(z-\zeta_{2}\right)+D_{i} \zeta_{2} $

(21)

式中：i=1代表上层流体；i=2代表下层流体。本文算法可以直接得到上、下层流体的层平均速度，可以更方便地求解MCC-FS模型的垂向速度。

在数值计算中，MCC-FS模型本文算法更为直接方便，在MATHMATICA中，得到数值结果的时间要小于1 min，计算耗时很短。

3 MCC-FS模型速度场算法结果与分析

本文基于Kodaira等^[5]的算例，比较了本文算法、原始算法、实验值得到的内孤立波的波形以及波速结果，以验证本文算法的准确性。并基于邹丽等^[7]的PIV实验结果，对本文算法以及给出的MCC-FS模型速度场求解方法进行进一步验证。最后，研究了不同密度比对内孤立波特性的影响。

3.1 内孤立波速度场研究

Kodaira等^[5]的算例中，上层流体的密度为ρ₁=856 kg/m³，厚度为h₁=0.05 m；下层流体的密度为ρ₂=996 kg/m³，厚度为h₂=0.25 m，对幅值为a/h₁为-0.50、-0.77、-0.99以及-1.21的内孤立波进行计算。内孤立波波形结果如图 3所示。图 3中，本文算法与Kodaira等^[5]的原始算法得到的波形结果完全重合，与实验值相比，波形结果也很准确。

图 4为不同波幅内孤立波的波速曲线，其中，速度无量纲化参数c₀为：

$ c_{0}=\sqrt{g h_{1} h_{2}\left(\rho_{2}-\rho_{1}\right) /\left(\rho_{1} h_{2}+\rho_{2} h_{1}\right)} $

(22)

可以看出，本文算法与Kodaira等^[5]得到的内孤立波波速结果也完全重合。需要指出的是，对于该算例，Kodaira等^[5]并未给出速度场的数值结果与实验结果。

邹丽等^[7]的算例中，上层流体的密度为ρ₁=941 kg/m³，厚度为h₁=0.05 m；下层流体的密度为ρ₂=1 003 kg/m³，厚度为h₂=0.25 m；幅值为a/h₁=-0.875, -1.323的内孤立波波形的实验结果，MCC-FS模型本文算法与实验结果对比如图 5所示。在该实验布置下，MCC-FS模型本文算法得到的波形结果同实验值基本一致。

邹丽等^[7]基于PIV技术得到了内孤立波的速度场结果。其中，实验中内孤立波幅值分别为a/h₁=-0.24, -0.59, -0.70, -0.96。需要说明的是，邹丽等^[7]并未直接给出内孤立波的波幅，本文中内孤立波的波幅是对文章中给出的内孤立波波面描点后得到的。图 6中，基于MCC-FS模型本文算法得到的速度场与邹丽等^[7]得到的实验速度场基本一致。

针对幅值为a/h₁=-0.59的内孤立波，邹丽等^[7]给出了不同深度下水平流速时历曲线以及波峰竖直剖面水平速度分布的实验结果。

图 7中，ρ₁/ρ₂=0.938, h₁/h₂=1/5, a/h₁=-0.59时，MCC-FS模型本文算法预报的不同深度的流体质点水平速度时历曲线与邹丽等^[7]的实验结果很接近。

图 8中，MCC-FS模型本文算法得到的波峰所在竖直剖面水平速度分布的数值结果同实验值吻合的很好，尤其是下层流体速度的预报上，基本与实验值重合。

3.2 不同密度比的浅水模式内孤立波特性研究

在不同密度比内孤立波特性的研究上，算例的布置为上层流体的厚度h₁=0.05 m，下层流体的厚度h₂=0.25 m，下层流体的密度ρ₂=1 000 kg/m³，改变上层流体的密度ρ₁研究不同密度比的影响，上层密度ρ₁=800, 900, 999 kg/m³，内孤立波幅值a=-0.025 m。这里需要说明的是，由于速度的无量纲化参数c₀在不同的上、下层密度比是不同的，因此，接下来展示的结果均有量纲。

图 9中，h₁/h₂=1/5, a/h₁=-0.5时，上、下层密度越接近，得到的内孤立波波形越宽，上、下层密度相差越大，得到的内孤立波波形越窄。图 10中，h₁/h₂=1/5, a/h₁=-0.5时，在上、下层密度相差不大时，自由面波形幅值很小，随着上、下层密度相差的越大，得到的自由面幅值越大。上、下层密度比为0.999、0.9以及0.8时，自由面最大位移ζ₁/h₁分别为0.003、0.033 2、0.071 9。这一结论同样也解释了对于上、下层流体密度接近的情况，引入刚盖假设是合理的，而在上、下层密度相差较大时需要考虑自由面的影响。

图 11中，h₁/h₂=1/5, a/h₁=-0.5时随着密度相差越大，上、下层流体质点的水平速度变得更大，内波界面上的速度跳跃值也越大。图 12中，h₁/h₂=1/5, a/h₁=-0.5时，随着上、下层密度比的增大，内孤立波的波速变小，上下层密度比为0.8到0.95时，波速随密度比近似成线性变化，密度比为0.95~1时，随着密度比增大，波速减小的较快。

4 结论

1) 通过与原算法的结果、文献[5, 7]实验的波形结果进行比较，表明了本文提出的本文算法准确可靠，本文算法可同时求得内波波形、自由表面波波形以及上、下层平均速度，更为直接、方便。

2) 通过与文献[7]实验的速度场结果对比，证明了MCC-FS模型可以准确的预报浅水模式内孤立波的速度场。

3) 基于MCC-FS模型的速度场结果，发现内孤立波波峰位置的竖直剖面上，随着上、下层密度差增大，上、下层流体质点的水平速度变得更大，内波界面上的速度跳跃值也增大。