有限深两层流体内孤立波非线性数值模拟入口边界条件研究

王 旭，尤云祥，黄文昊

(上海交通大学海洋工程国家重点实验室，上海 200240)

内孤立波是最大振幅发生在密度分层海水内部的一种常见海洋动力学现象，不仅是海洋能量级串中的一个重要环节，也是对海洋工程结构物安全性产生重要影响的一类海洋环境因素。在南海流花油田的早期延长测试期间，就曾发生过因内孤立波产生的突发性强流而导致缆绳拉断、船体碰撞，甚至拉断和挤破漂浮软管等事故［1］。在南海陆丰油田的早期延长测试期间，也曾发生过内孤立波产生的突发性强流使半潜式钻井船与锚定油轮在连接输油管道时发生困难等问题［2］。

非线性和频散效应是密度分层海洋中影响内孤立波传播演化特性的两个基本因素，非线性效应趋向于使波形变陡，而频散效应则趋向于使波形变得平坦，内孤立波的稳定传播是非线性与频散效应在一定尺度上达到平衡的结果，可以用KdV、eKdV和MCC等理论模型来描述［3］。设ε=a/h为非线性参数，μ=( h /λ)2为色散参数(其中a 和λ分别为内孤立波振幅和特征宽度，h为垂向特征尺度)，KdV理论是建立在弱非线性、弱色散且两者平衡(即ε=O(μ)≪1)条件下的一类内孤立波理论［4］。研究表明，在小振幅情况下KdV理论与实验结果相符，但在大振幅情况KdV理论明显与实验结果不符［5-8］。

随着内孤立波振幅的增大，高阶非线性项的影响会逐渐增大，而且非线性和色散性不再保持平衡［9］，这正是KdV理论不能适用于大振幅内孤立波的主要原因之一。这个缺陷可以通过在KdV方程中加入立方非线性项等途径解决，所得修正方程称为eKdV方程［10］。但eKdV理论依然是建立在弱非线性和弱色散条件下的，为克服这两类理论中需要弱非线性这个限制性条件的缺陷，从完全非线性欧拉方程出发，Choi和Camassa［11］建立了一个强非线性和弱色散的两层Green-Naghdi模型，在定态条件下该理论可以简化为一个常微分方程形式的非线性方程，称为MCC理论。这是建立在强非线性和弱色散条件下的一类定态内孤立波理论，研究表明在大振幅内孤立波的情况，MCC理论与实验结果相符［12］。

上述KdV、eKdV和MCC理论均是在弱色散条件下建立的，而对KdV和eKdV理论还需要弱非线性这个限制性条件，但关于弱非线性和弱色散这两个限制性条件迄今仍是定性的。最近，黄文昊等［13］利用大型密度分层水槽，采用双推板内孤立波造波方法，研究了弱非线性和弱色散这两个条件的定量表征问题，结果表明在给定密度分层及内孤立波振幅条件下，KdV、eKdV和MCC理论的适用性与非线性和色散参数的某种组合有关。

各种海洋环境条件下深海浮式平台载荷特性及其运动响应的合理确定，是其设计与应用中的一项关键性指标。程友良等［14］和蔡树群等［15-16］将Morison公式与KdV理论结合，而XIE等［17］则将Morison公式与MCC理论结合，研究了内孤立波作用下直立贯底圆柱体的载荷特性问题。最近，尤云祥等［18-19］将Morison公式与eKdV理论结合，研究了内孤立波作用下张力腿和半潜式平台的载荷与动力响应问题，而宋志军等［20］将Morison公式与KdV理论结合，研究了内孤立波作用下Spar平台的载荷与动力响应问题。在这些研究中，Morison公式中的惯性力和拖曳力系数是参照表面波的方法选取的，这种选取方法尚缺乏理论和实验依据。

随着计算流体力学(CFD)技术的发展，采用CFD方法来研究内孤立波的生成传播及其与海上结构物的相互作用问题，已成为目前国际上的研究热点之一。由于采用CFD方法可以直接获得内孤立波与海洋结构物相互作用过程中的速度场及压力场的变化特性，因此可以直接获得内孤立波作用下海上结构物的载荷等水动力特性。在这种方法中，通常在入口边界条件中以定态内孤立波的解析解(包括KdV、eKdV和MCC等)来驱动完全非线性模型，进而模拟内孤立波的生成传播及其与海上结构物的作用特性。

陈钰等［21］、关辉等［22］和李水娟等［23］利用 KdV 理论解作为入口边界条件，而高原雪等［24］利用 MCC 理论解作为入口边界条件，采用CFD方法研究了内孤立波的生成传播问题。付东明等［25］利用KdV理论解作为入口边界条件，而陈杰等［26］和关辉等［27］利用eKdV理论解作为入口边界条件，采用CFD方法研究了内孤立波作用下水下潜体的载荷特性问题。刘碧涛等［28］利用eKdV理论解作为入口边界条件，采用CFD方法研究了内孤立波作用下深海立管的载荷特性问题。

根据文献［13］的实验结果，在给定密度分层及内孤立波振幅条件下，由于KdV、eKdV和MCC理论的适用性与非线性和色散参数的某种组合有关，因此如果选择不合适的内孤立波理论解作为CFD数值模拟的入口边界条件，则有可能会导致振幅及波形都不可控的数值模拟结果。因此，在给定密度分层及内孤立波振幅条件下，如何选择合适的内孤立波理论解作为CFD数值模拟的入口边界条件是利用CFD方法来研究内孤立波的生成传播及其与海上结构物相互作用问题时需要研究和解决的关键之一。

有鉴于此，以有限深两层流体中内孤立波诱导上下层深度平均水平速度为入口边界条件，以理想流体完全非线性欧拉方程为控制方程，结合VOF动态界面跟踪方法，对内孤立波的生成特性进行完全非线性数值模拟。在此基础上，研究在给定密度分层及内孤立波振幅条件下如何选择合适的内孤立波理论解作为完全非线性数值模拟入口边界条件等问题。

1 数值模拟方法

设两层流体均为理想不可压缩流体，当流体处于静平衡状态时，上层流体深度与密度分别为h1和ρ1，而下层流体深度与密度分别为h2和ρ2，水深为h=h1+h2。建立直角坐标系oxz，其中ox轴位于流体静止时两层流体的界面上，oz轴垂直向上为正。对水表面作刚盖假定，假定水底为不可渗透平坦的刚性固壁。设ζ为两层流体界面位移，ui、wi为流场中的速度矢量，其中i=1、2分别表示上下层流体。

在上下层流体中，流场控制方程如下:

其中，g为重力加速度，各物理量的下标t、x和z分别表示关于时间和空间坐标的导数，pi为压力。

在两层流体界面上，要求满足如下连续性条件:

在水面和水底，要求满足如下固壁条件:

流场计算的控制区域如图1所示，包括内孤立波生成传播区和消波区两个区域。采用速度入口方法生成内孤立波，当形成稳定的内孤立波后，对所生成内孤立波的传播特性进行监测分析。

图1 流场计算控制区域Fig.1 Governing domain for flow computation

速度入口方法为:在上层流体入口处，入口速度取为u1;在下层流体入口处，入口速度取为。内孤立波在生成与传播过程中，两层流体的界面会发生变化，采用VOF(volume of fluid)方法追踪两层流体界面的变化［29］。在VOF方法中，用体积分数aq表示单元内第q相流体所占体积与该单元的体积之比，其中q=1表示上层流体，q=2表示下层流体。当aq=1时，表示单元内全部为第q相流体，而当aq=0时，表示单元不存在第q相流体，对含界面的单元则有0＜aq＜1。在每一个流体单元中各相的体积分数之和为1，即=1，而且对于任一相的体积分数aq，要求满足如下连续方程:

采用海绵层方法对水槽末端的内孤立波进行消波处理，该方法通过添加人工粘性项实现阻尼消波之目的。为此，对消波区，在动量方程的右端添加源项，将其修改为

其中，海绵层衰减系数μ依据文献［30］中的方法确定。

采用KdV、eKdV和MCC理论计算内孤立波诱导上下层流体中的层深度平均水平速度。为此，设a和λ分别为内孤立波振幅与特征宽度，定义非线性参数ε=a/h和色散参数μ=(h/λ)2，那么在弱色散和弱非线性(即ε≪1而且μ≪1)条件下，界面位移ζ满足如下的eKdV方程［3］

其中

当内孤立波为定态平面前进波时，eKdV方程(10)有如下形式的解［3］:

其中

其中，cekdv为内孤立波eKdV理论解的相速度。

当ε=O(μ)≪1，即当内孤立波是弱非线性、弱色散且两者平衡时，在eKdV方程中的立方非线性项可以忽略，即可取c3=0，这时eKdV理论解(13)即退化为如下的KdV理论解［3］

其中，ckdv为内孤立波KdV理论解的相速度。

在上述理论中，均要求内孤立波是弱非线性和弱色散的，对KdV理论还要求两者是平衡的。因此，这两类内孤立波理论一般只适用于振幅较小的情况。为克服这个缺陷，Choi和Camassa［11］建立了一个强非线性和弱色散的内孤立波理论模型，称为MCC理论模型，其定态解可表示为

其中，cMCC为内孤立波MCC理论解的相速度，

a-和a+(a-＜a+)分别为下面方程的两个根

式中:

对于KdV理论解来说，只要h1/h2≠ ( h1/h2)c，则不存在极限振幅，即无论振幅有多大，该理论均能给出对应的内孤立波解。Camassa等［12］指出，KdV理论无极限振幅的这个缺陷，正是其不能适用于大振幅内孤立波的原因之一。与KdV理论的一个重要不同之处是，其它两类内孤立波理论均存在极限振幅，其中eKdV和MCC理论解的极限振幅分别为［3］

2 结果与分析

文中将结合该文献中的相关实验结果进行数值模拟与分析。为此，数值水槽主尺度、上下层流体密度及其深度比均与该文一致。在数值模拟中，内孤立波生成传播区的长度为24 m，而消波区的长度为6 m。生成传播区内网格的纵向间距为0.03 m，垂向网格划分以底部向上0.4 m位置为界，上部区域取网格垂向间距为0.005 m，下部区域首层垂向间距取为0.005 m，为保证网格平滑过渡，后续网格间距按照1.02的比例逐渐加宽;消波区内垂向网格间距同生成传播区的情况一致，而纵向首层网格间距取0.03 m，后续网格间距按1.04的比例逐渐加宽;整个计算域的网格单元共计60 800个。

利用商业软件Fluent进行数值模拟，初始时水槽中流场物理量均取为零，时间步长为0.01 s。计算区域采用结构化网格进行离散，在两层流体界面附近，网格需进行加密处理。控制方程采用有限体积法(FVM)进行离散，对流项采用二阶迎风格式离散，时间步推进选用二阶隐式时间匹配格式，两相界面的构造方法选用几何重构法，应用非定常压力的隐式算子分割法(PISO)进行压力-速度的耦合计算。

文献［13］的实验结果表明，存在一个临界色散参数μ0=0.1，以此为基准对弱色散条件μ≪1进行定量化，可将其分为弱色散(μ＜μ0)和强色散(μ≥μ0)两种情况;在弱色散(μ＜μ0)条件下，对非线性条件也进行定量化，可将其分为弱非线性(ε≤μ)、中等非线性(μ＜ε≤)和强非线性(ε＞)三种情况;在系列实验基础上，获得了KdV、eKdV和MCC理论的如下适用性条件:KdV理论适用于弱非线性和弱色散情况，eKdV理论适用于中等非线性和弱色散情况，而MCC理论适用于强非线性和弱色散或强色散情况。

在三种上下层流体深度比情况下，文献［13］中系列实验所得内孤立波振幅如表1所示。在本文数值模拟中，初始内孤立波的设计振幅ad也取为表中所列相应的值。由式(21)可知，当h1∶h2=3∶7和1∶9时，h=-0.182 6和 -0.211 8，而/h=-0.196 7和 -0.396 7;当 h1∶h2=2∶8时，/h =-0.234 1，而/h=-0.296 7。由此可知，当h1∶h2=3∶7和1∶9时，各设计振幅均在eKdV和MCC理论解的极限振幅内;而当h1∶h2=2∶8时，各设计振幅虽均在MCC理论解的极限振幅内，但最后三个振幅已超出eKdV理论解的极限振幅。

表1 数值模拟中的无因次设计振幅ad/hTab.1 Non-dimensional design amplitudes for numerical simulations ad/h

在数值模拟中，对给定的设计振幅ad，利用KdV、eKdV和MCC理论解分别计算其诱导的上下层流体中的层深度平均水平速度(见式(6))，将其作为入口边界条件来驱动完全非线性数值模型(见式(1)-(5))。在上述初始内孤立波设计振幅ad下，这里重点研究如何选择合适的定态内孤立波理论解来驱动完全非线性数值模型的问题。

为此，设εm和μm分别为数值模拟所得波的非线性和色散参数，而εd和μd分别为在设计振幅ad下利用内孤立波理论解计算得到的非线性和色散参数。在图2中，给出了当h1/h2=3∶7时，(εm，μm)和(εd，μd)的结果比较。其中:“▲”表示(εm，μm)，“2”表示(εd，μd)。在给定设计振幅ad下，点▲和2之间的横向距离表示μm和μd之间的绝对误差，而点▲和2之间的垂向距离表示εm和εd之间的绝对误差。由图可知，在给定设计振幅ad下，利用eKdV理论解作为入口边界条件，所得μm和μd以及εm和εd之间的绝对误差最小。因此，在表1所示各设计振幅下，利用eKdV理论解作为入口边界条件，数值模拟所得内孤立波振幅及设计振幅误差最小，而且两者波形最为吻合。

图2 当 h1∶h2=3∶7 时 (εm，μm)和 (εd，μd)结果比较Fig.2 Comparisons between(εm，μm)and(εd，μd)when h1∶h2=3∶7

由图2(a)可知，在各设计振幅ad下，KdV理论解和数值模拟所得波的非线性和色散参数均属于弱色散和中等非线性区域，而文献［13］的实验表明，KdV理论适用于弱色散和弱非线性的情况，因此采用KdV理论解作为入口边界条件本身就是不合适的。由图2(b)可知，在各设计振幅ad下，eKdV理论解和数值模拟所得波的非线性和色散参数也均属于弱色散和中等非线性区域，而文献［13］的实验表明，eKdV理论适用于弱色散和中等非线性的情况，两者正好是一致的，即采用eKdV理论解作为入口边界条件是合适的。由图2(c)可知，MCC理论解和数值模拟所得波的非线性和色散参数均属于弱色散和中等非线性区域，而文献［13］的实验表明，MCC理论适用于弱色散和强非线性或强色散的情况，因此采用MCC理论解作为入口边界条件本身就是不合适的。

在图3中，给出了当h1/h2=3∶7时，数值模拟结果与初始内孤立波及相关实验结果的比较。由图中第一列可知，当利用KdV理论解作为入口边界条件时，数值模拟所得波的振幅及其波形与初始输入波的振幅及其波形差异都很大，而且随着设计振幅的增大，数值模拟结果还会出现显著的尾波现象。由图中第二列可知，当利用eKdV理论解作为入口边界条件时，数值模拟所得波的振幅及其波形与初始输入波和文献［13］中实验所得波振幅及其波形之间的差异都很小，而且数值模拟结果没有明显的尾波现象。由图中第三列可知，与利用KdV理论解作为入口边界条件的结果相比，当利用MCC理论解作为入口边界条件时，数值模拟所得波的振幅及其波形与初始输入波的振幅及其波形之间的差异明显减小，但这种差异要比利用eKdV理论解作为入口边界条件的结果大，因此采用eKdV理论解作为入口边界条件更为合适。

由图2可知，当h1/h2=3∶7时，在各设计振幅ad下，利用KdV理论解作为入口边界条件，数值模拟所得波均属于弱色散和中等非线性的。因此，根据文献［13］中的实验结果，如果该波为内孤立波，则应该属于eKdV型的。为此，在图4中，给出了当ad/h=0.034 0、0.105 5和0.158 0时，利用KdV理论解作为入口边界条件所得数值模拟结果与等振幅eKdV理论解的比较。由图4(a)可知，在较小无因次振幅的情况，数值模拟结果的波形与等振幅eKdV理论解波形一致，这表明数值模拟所得波为eKdV型内孤立波。由图4(b)和4(c)可知，在较大无因次振幅的情况，数值模拟所得波的波形与等振幅eKdV理论解的波形差异很大，这表明在这种情况下利用KdV理论解作为入口边界条件是不合适的。

图3 当h1∶h2=3∶7时数值模拟结果与初始内孤立波及相关实验结果比较Fig.3 Comparisons of numerical results with theoretical and experimental ones when h1∶h2=3∶7

图4 当h1∶h2=3∶7时利用KdV理论解作为入口边界条件所得数值模拟结果与eKdV理论解比较Fig.4 Comparisons of numerical results with eKdV theoretical solutions based on the KdV input conditions when h1∶h2=3∶7

图5给出了当h1/h2=2∶8时，(εm，μm)和(εd，μd)的结果比较。由图可知，对表1中前3个设计振幅ad，利用KdV理论解作为入口边界条件，所得μm和μd以及εm和εd之间的误差最小，而且利用KdV理论解和数值模拟所得波的非线性和色散参数均属于弱色散和弱非线性区，与文献［13］实验结果一致，即这时采用KdV理论解作为入口边界条件最为合适。对表1中第4到9个设计振幅ad，利用eKdV理论解作为入口边界条件，所得μm和μd以及εm和εd之间的误差最小，而且利用eKdV理论解和数值模拟所得波的非线性和色散参数均属于弱色散和中等非线性区，与文献［13］实验结果一致，即这时采用eKdV理论解作为入口边界条件最为合适。对表1中最后4个设计振幅ad，利用MCC理论解作为入口边界条件，所得μm和μd以及εm和εd之间的误差最小，而且利用MCC理论解和数值模拟所得波的非线性和色散参数均属于弱色散和强非线性区域，与文献［13］实验结果一致，即这时采用MCC理论解作为入口边界条件最为合适。

图5 当 h1∶h2=2∶8 时 (εm，μm)和 (εd，μd)结果比较Fig.5 Comparisons between(εm，μm)and(εd，μd)when h1∶h2=2∶8

在图6中，给出了当h1/h2=2∶8时，数值模拟结果与初始内孤立波及相关实验结果的比较。由图可知，当ad/h=0.062 1时，与利用eKdV和MCC理论解作为入口边界条件的结果相比，利用KdV理论解作为入口边界条件所得数值模拟波的振幅及其波形与初始输入波最为吻合，而且与文献［13］中的实验一致。当ad/h=0.130 4时，与利用KdV和MCC理论解作为入口边界条件的结果相比，利用eKdV理论解作为入口边界条件所得数值模拟波的振幅及其波形与初始输入波最为吻合，而且与文献［13］中的实验一致。当ad/h=0.243 3时，与利用KdV和eKdV理论解作为入口边界条件的结果相比，利用MCC理论解作为入口边界条件所得数值模拟波的振幅及其波形与初始输入波最为吻合，而且与文献［13］中的实验一致。

图6 当h1∶h2=2∶8时数值模拟结果与初始内孤立波及相关实验结果比较Fig.6 Comparisons of numerical results with theoretical and experimental ones when h1∶h2=2∶8

由图5可知，当h1/h2=2∶8时，对表1中最后10个设计振幅ad，利用KdV理论解作为入口边界条件，数值模拟所得波均属于弱色散和中等非线性的。因此，根据文献［13］中实验结果，如果该波为内孤立波，则应该属于eKdV型的。为此，在图7中，给出了当ad/h=0.130 4和0.243 3时，利用KdV理论解作为入口边界条件所得数值模拟结果与等振幅eKdV理论解的比较。由图可知，数值模拟结果的波形与等振幅eKdV理论解波形一致，这表明数值模拟所得波为eKdV型内孤立波。

图7 当h1∶h2=2∶8时，利用KdV理论解作为入口边界条件所得数值模拟结果与eKdV理论解比较Fig.7 Comparisons of numerical results with eKdV theoretical solutions based on the KdV input conditions when h1∶h2=2∶8

图8 当 h1∶h2=1∶9 时 (εm，μm)和 (εd，μd)结果比较Fig.8 Comparisons between(εm，μm)and(εd，μd)when h1∶h2=1∶9

在图8中，给出了当h1∶h2=1∶9时，(εm，μm)和 (εd，μd)的结果比较。由图可知，在给定设计振幅 ad下，利用MCC理论解作为入口边界条件，所得μm和μd以及εm和εd之间的误差最小。因此，在表1所示各设计振幅下，利用MCC理论解作为入口边界条件，数值模拟所得内孤立波振幅与设计振幅误差最小，而且两者波形最为吻合。由图8(a)和8(b)可知，在各设计振幅ad下，KdV和eKdV理论解和数值模拟所得波的非线性和色散参数均属于强弱色散区域。文献［13］的实验表明，KdV理论适用于弱色散和弱非线性的情况，而eKdV理论适用于弱色散和中等非线性的情况，因此采用KdV和eKdV理论解作为入口边界条件本身就是不合适的。由图8(c)可知，在各设计振幅ad下，利用MCC理论解和数值模拟所得波的非线性和色散参数均属于强弱色散区域，而文献［13］的实验表明，MCC理论适用于强弱色散的情况，两者正好是一致的，即采用MCC理论解作为入口边界条件是合适的。

在图9中，给出了当h1∶h2=1∶9时，数值模拟结果与初始内孤立波及相关实验结果的比较。由图中第一列和第二列可知，当利用KdV和eKdV理论解作为入口边界条件时，随着设计振幅的增大，数值模拟所得波的振幅及其波形与初始输入波的振幅及其波形之间的差异也随着增大。由图中第三列可知，当利用MCC理论解作为入口边界条件时，数值模拟所得波的振幅及其波形与初始输入波和文献［13］中实验所得波振幅及其波形之间的差异都很小。

图9 当h1∶h2=1∶9时数值模拟结果与初始内孤立波及相关实验结果比较Fig.9 Comparisons of numerical results with theoretical and experimental ones when h1∶h2=1∶9

图10 当h1∶h2=1∶9时利用KdV和eKdV理论解作为入口边界条件所得数值模拟结果与MCC理论解比较Fig.10 Comparisons of numerical results with MCC theoretical solutions based on the KdV and eKdV input conditions when h1∶h2=1∶9

由图8可知，当h1/h2=1∶9时，在各设计振幅ad下，利用KdV和eKdV理论解作为入口边界条件，数值模拟所得波均属于强色散的。因此，根据文献［13］中的实验结果，如果该波为内孤立波，则应该属于MCC型的。为此，在图10中，给出了当ad/h=0.041 8，0.137 3和0.151 3时，利用KdV和eKdV理论解作为入口边界条件所得数值模拟结果与等振幅MCC理论解的比较。由图可知，数值模拟结果的波形与等振幅MCC理论解波形一致，这表明数值模拟所得波为MCC型内孤立波。

3 结语

这里为有限深水内孤立波非线性数值模拟中如何选择合适的理论解作为入口边界条件提供了理论依据。所得入口边界条件选择方法仅限于有限深水的情况，针对无限深水条件下应如何选择合适的内孤立波理论解作为入口边界条件等问题，尚需做进一步研究。

［1］ 陈景辉.南海流花11-1深海油田开发工程［J］.中国海洋平台，1996，11(1):43-45.(CHEN Jing-hui.Liuhua 11-1 deep sea oil-field development project in South China Sea［J］.China Offshore Platform，1996，11(1):43-45.(in Chinese))

［2］ Bole J B，Ebbesmeyer C C，Romea R D.Soliton currents in the South China Sea:measurements and theoretical modeling［C］//The 16thOffshore Technology Conference.1994:367-376.

［3］ Helfrich K R，Melville W K.Long nonlinear internal waves［J］.Ann.Rev.Fluid Mech.，2006，38:395-425.

［4］ Ostrovsky L A，Stepanyants Yu A.Do internal solitons exist in the ocean［J］.Rev.Geophys.，1989，27(3):293-310.

［5］ Michallet H，Barthelemy E.Experimental study of interfacial solitary waves［J］.J.Fluid Mech.，1998，366:159-177.

［6］ Grue J，Jensen A，Rusas P E，et al.Properties of large-amplitude internal waves［J］.J.Fluid Mech.，1999，380:257-278.

［7］ Sveen J K，Guo Y，Davies P A，et al.On the breaking of internal solitary waves at a ridge［J］.J.Fluid Mech.，2002，469:161-188.

［8］ Walker S A，Martin A J，Easson W J，et al.Comparison of laboratory and theoretical internal solitary wave kinematics［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，2003，129(5):210-218.

［9］ Brandt P，Rubino A，Alpers W，et al.Internal waves in a strait of messina studied by a numerical model and synthetic aperture radar form the ERS1/2 satellites［J］.J.Phys.Oceanography.，1997，27(5):648-663.

［10］Djordjevic V D.The fission and disintegration of internal solitary waves moving over two-dimensional topography［J］.J.Phys.Oceanography，1978，8(6):1016-1024.

［11］Choi W，Camassa R.Fully nonlinear internal waves in a two-fluid system［J］.J.Fluid Mech.，1999，396:1-36.

［12］Camassa R，Choi W，Michallet H，et al.On the realm of validity of strongly nonlinear asymptotic approximations for internal waves［J］.J.Fluid Mech.，2006，549:1-23.

［13］黄文昊，尤云祥，王旭，等.有限深两层流体中内孤立波造波实验及其理论模型［J］.物理学报，2013，62(8):1-14.(HUANG Wen-hao，YOU Yun-xiang，WANG Xu，et al.Wave-making experiments and theoretical models for internal solitary waves in a two-layer fluid of finite depth［J］.Acta Phys.Sin.，2013，62(8):1-14.(in Chinese))

［14］CHENG Y L，LI J C，LIU Y F.The induced flow field by internal solitary wave and its action on cylindrical piples in the stratified ocean［C］//The 4th International Conference on Fluid Mechanics.2004:296-299.

［15］CAI S Q，LONG X M，GAN Z J.A method to estimate the forces exerted by internal solitons on cylinder piles［J］.Ocean Engineering，2003，30(5):673-689.

［16］CAI S Q，LONG X M，WANG S G.Force and torques exerted by internal solitions in shear flows on cylinder piles［J］.Applied Ocean Research，2008，30(1):72-77.

［17］XIE J S，JIANG Y J，YANG L G.Strongly nonlinear internal solution load on a small vertical circular cylinder in two-layer fluids［J］.Applied Mathematical Modelling，2010，34(8):2089-2101.

［18］尤云祥，李巍，何景异，等.海洋内孤立波中张力腿平台的水动力特性［J］.上海交通大学学报，2010，44:12-17.(YOU Yun-xiang，LI Wei，HE Jing-yi，et al.Hydrodynamic characteristics of tension leg platforms in ocean internal solitary waves［J］.Journal of Shanghai Jiao Tong University，2010，44:12-17.(in Chinese))

［19］尤云祥，李巍，胡天群，等.内孤立波中半潜平台动力响应特性［J］.海洋工程，2012，30(2):1-7.(YOU Yun-xiang，LI Wei，HU Tian-qun，et al.Dynamic responses of a semi-submersible platform in internal solitary waves［J］.The Ocean Engineering，2012，30(2):1-7.(in Chinese))

［20］宋志军，勾莹，滕斌，等.内孤立波作用下Spar平台的运动响应［J］.海洋学报，2010，32:12-19.(SONG Zhi-jun，GOU Ying，TENG Bin，et al.The motion responses of a Spar platform under internal solitary wave［J］.Acta Oceanologica Sinica，2010，32(2):12-19.(in Chinese))

［21］陈钰，朱良生.基于 FLUENT的海洋内孤立波数值水槽模拟［J］.海洋技术，2009，28(4):72-75.(CHEN Yu，ZHU Liang-sheng.Numerical wave tank simulation of oceanic internal solitary waves based on FLUENT［J］.Ocean Technology，2009，28(4):72-75.(in Chinese))

［22］关晖，苏晓冰，田俊杰.分层流体中内孤立波数值造波方法及其应用［J］.力学季刊，2011，32(2):218-223.(GUAN Hui，SU Xiao-bing，TIAN Jun-jie.Numerical method of solitary wave generation in a stratified flow and its applications［J］.Chinese Quarterly Of Mechanics，2011，32(2):218-223.(in Chinese))

［23］李水娟，余真真，罗龙洪，等.三维数值波浪水槽内孤立波特性分析［J］.水电能源科学，2012，30(7):96-99.(LI Shui-juan，YU Zhen-zhen，LUO Long-hong，et al.Analysis of internal solitary waves characteristics of 3d numerical wave flume［J］.Water Resources and Power，2012，30(7):96-99.(in Chinese))

［24］高原雪，尤云祥，王旭，等.基于 MCC理论的内孤立波数值模拟［J］.海洋工程，2012，30(4):29-36.(GAO Yuan-xue，YOU Yun-xiang，WANG Xu，et al.Numerical simulation for the internal solitary wave based on MCC theory［J］.The Ocean Engineering，2012，30(4):29-36.(in Chinese))

［25］付东明，尤云祥，李巍.两层流体中内孤立波与潜体相互作用数值模拟［J］.海洋工程，2009，27(3):38-44.(FU Dong-ming，YOU Yun-xiang，LI Wei.Numerical simulation of internal solitary waves with submerged body in a two-layer fluid［J］.The Ocean Engineering，2009，27(3):38-44.(in Chinese))

［26］陈杰，尤云祥，刘晓东，等.内孤立波与有航速潜体相互作用数值模拟［J］.水动力学研究与进展，2010，25(3):344-350.(CHEN Jie，YOU Yun-xiang，LIU Xiao-dong，et al.Numerical simulation of interaction of internal solitary waves with a moving submarine［J］.Chinese Journal of Hydrodynamics，2010，25(3):344-350.(in Chinese))

［27］关晖，魏岗，杜辉.内孤立波与潜艇相互作用的水动力学特性［J］.解放军理工大学学报:自然科学版，2012，13(5):577-582.(GUAN Hui，WEI Gang，DU Hui，et al.Hydrodynamic properties of interactions of three-dimensional internal solitary waves with submarine［J］.Journal of PLA University of Science and Technology，Natural Science Edition，2012，13(5):577-582.(in Chinese))

［28］刘碧涛，李巍，尤云祥，等.内孤立波与深海立管相互作用数值模拟［J］.海洋工程，2011，29(4):1-7.(LIU Bi-tao，LI Wei，YOU Yun-xiang，et al.Numerical simulation of interaction of internal solitary waves with deep-sea risers［J］.The Ocean Engineering，2011，29(4):1-7.(in Chinese))

［29］Hirt C W，Nichols B D.Volume of fluid(VOF)method for the dynamics of free boundaries［J］.J.Computational Physics，1981，39(1):201-225.

［30］韩鹏.基于VOF方法的不规则波阻尼消波研究［D］.大连:大连理工大学，2008.(HAN Peng.The study of damping absorber for irregular waves based on VOF method［D］.Dalian:Dalian University of Technology，2008.(in Chinese))