0 引　言

内波是指发生在密度稳定层化的海水内部的一种波动现象，其波长和周期分布在很宽的范围内，常见波长为几十米至几百千米，周期为几分钟至几十小时，最大振幅可达百米，且出现在海洋内部，不易察觉。大振幅内孤立波是海洋中活动最为频繁的一种强非线性潮成内波，其不仅能使等密度面大振幅垂向起伏，更会产生具有强垂向剪切的往复水平流，对潜水器的水下适航性和安全性构成严重威胁。我国南海是内孤立波多发海域，因此对内孤立波作用下潜水器的运动响应研究具有现实意义。

国内外关于内孤立波的理论、实验及数值造波的研究已经相当深入，在小振幅浅水假定条件下，D. J. Korteweg及G. de Vrie首次提出了描述内孤立波传播的KdV数学模型，在理论上对内孤立波现象进行了解释。BENJAMIN^[1]推导了浅水中传播的内孤立波KdV方程。由于内孤立波的非线性较强，缺乏非线性项的KdV模型并不能广泛适用，MELVILLE等^[2]提出了eKdV理论，该理论可以较好地解释浅水中多数工况下的内孤立波。李效民^[3]基于KdV，mKdV和eKdV三种理论模型，分别采用双推板、平板拍击和速度入口边界3种数值造波方法，对内孤立波进行了模拟，并将数值模拟结果与理论及实验结果进行对比，结果表明速度入口边界法造波效果及造波效率最优；平板拍击法造波效果好，但造波效率略差；双推板法造波效果及造波效率均较差。

有关内孤立波与水下结构物的相互作用，相关研究人员在数值和实验上已有大量研究成果。武军林^[4]运用重力式分层流水槽开展横置细长潜体遭遇内孤立波的流场和作用力特性实验，发现潜体所受水平作用力主要由内孤立波诱导流场决定，垂向作用力由内孤立波引起的密度场和流场变化的耦合作用决定。关晖^[5]采用数值方法模拟了Suboff约束模型穿越内孤立波时的流场形态和受力载荷特性，研究表明Suboff在穿越内孤立波时在垂向和水平方向的受力都会在短时间内发生巨变，力矩的突变也有可能使潜艇发生倾覆。现有的研究主要集中在约束模型遭遇内孤立波时的受力特性，而内孤立波作用下潜水器的运动响应问题却鲜有研究，但潜水器在遭遇内孤立波时的“掉深”及姿态角变化才是工程实际中最需要解决的问题。

本文采用CFD方法，依据eKdV理论采用速度入口边界生成内孤立波，将数值结果与理论值进行对比，并分析了内孤立波的流场特点；然后模拟了4种不同潜深的潜水器在内孤立波作用下的运动响应，并对不同潜深的潜水器在内孤立波作用下的运动机理进行了深入分析，为后续内孤立波作用下潜水器的安全操控策略研究提供支撑。

1 内孤立波造波数值模拟 1.1 数值水池建立

内孤立波数值水池示意图如图1所示，采用2层流体模型，水池长46 m，高1.6 m，宽2 m。X轴为2层流体分界面，向右为正方向，上层定义为淡水，密度ρ₁=998 kg/m³，层深h₁=0.48 m；下层定义为海水，密度ρ₂=1025 kg/m³，层深h₂=1.12 m。水池前26 m为造波区，后20 m为消波区。Z轴沿水深方向，垂直向上为正方向，Y轴正方向依据右手坐标系确定。

1.2 内孤立波eKdV理论模型

根据内孤立波eKdV理论^[6]，内孤立波波形表达式如下式：

$ \eta =\frac{a}{B+(1-B){\mathit{cos}h}^{2}\left(\dfrac{x-ct}{\lambda }\right)} ，$

(1)

$ c={c}_{0}+\frac{a}{3}\left({c}_{1}+\frac{1}{2}{c}_{3}a\right) ，$

(2)

$ {c}_{0}=\sqrt{\frac{g{h}_{1}{h}_{2}({\rho }_{2}-{\rho }_{1})}{{\rho }_{1}{h}_{2}+{\rho }_{2}{h}_{1}}} ，$

(3)

$ {c}_{1}=-\frac{3{c}_{0}}{2}\frac{{\rho }_{1}{h}_{2}^{2}-{\rho }_{2}{h}_{1}^{2}}{{\rho }_{1}{{h}_{1}h}_{2}^{2}+{\rho }_{2}{{h}_{2}h}_{1}^{2}} ，$

(4)

$ {c}_{2}=\frac{{c}_{0}}{6}\frac{{\rho }_{1}{{h}_{2}h}_{1}^{2}+{\rho }_{2}{{h}_{1}h}_{2}^{2}}{{\rho }_{1}{h}_{2}+{\rho }_{2}{h}_{1}} ，$

(5)

$ {c}_{3}=\frac{7{c}_{1}^{2}}{18{c}_{0}}-\frac{{c}_{0}{(\rho }_{1}{h}_{2}^{3}+{\rho }_{2}{h}_{1}^{3})}{{h}_{1}^{2}{h}_{2}^{2}({\rho }_{1}{h}_{2}+{\rho }_{2}{h}_{1})} ，$

(6)

$ \lambda =\sqrt{\frac{24{c}_{2}}{a(2{c}_{1}+{c}_{3}a)}} ，$

(7)

$ B=\frac{-a{c}_{3}}{2{c}_{1}+a{c}_{3}} ，$

(8)

$ {U}_{1}=-c\frac{\eta }{{h}_{1}-\eta } ，$

(9)

$ {U}_{2}=c\frac{\eta }{{h}_{2}+\eta } 。$

(10)

式中：a为内孤立波波幅；λ为一个表征内孤立波长度的物理量，整个波形基本包含于约4λ的长度范围内；c为内孤立波传播的相速度；g为重力加速度，9.81 m/s²。

1.3 CFD计算方程

以RANS方程和时均连续性方程为基本方程，时均连续性方程具体形式为：

$ \frac{\partial {u}_{i}}{\partial {x}_{i}}=0 ，$

(11)

RANS方程具体形式为：

$ \frac{{\partial u}_{i}}{\partial t}+\frac{\partial }{\partial {x}_{j}}\left({u}_{i}{u}_{j}\right)=-\frac{1}{\rho }\frac{\partial p}{\partial {x}_{i}}+{S}_{i}+\frac{1}{\rho }\frac{\partial }{\partial {x}_{j}}\left(\mu \frac{\partial {u}_{i}}{\partial {x}_{j}}-\rho \overline{{u}_{i}'{u}_{j}'}\right) 。$

(12)

式中：ρ为流体密度；μ为流体粘度；p为静水压力；S_i为质量力；u_i和u_j为速度分量。

1.4 计算域设置及网格划分

数值水池顶部依据“刚盖”假定，设置为对称边界；底部设置为滑移壁面边界；左侧为入口边界，边界条件为速度入口，入口速度场函数采用自编函数表达式，表达式见式（9）和式（10）；右侧为出口边界，边界条件为压力出口；沿y轴方向的2个面均设置为对称边界。2层流体的交界面采用VOF法捕捉，消波方式采用阻尼消波，湍流模型采用可实现k-ε模型，采用隐式不定常求解器求解，求解时间步长为0.02 s。采用六面体网格对整个数值水池进行网格划分，总网格数为385万。

1.5 数值造波结果及分析

为了验证内孤立波数值水池的准确性，采用上述内孤立波数值水池生成设计波幅a=0.17 m的内孤立波，并将波谷传播至X=14 m处的数值模拟结果与设计波幅下的eKdV理论解进行比较，同时对比波谷传播至X=14 m，X=16 m，X=18 m的结果，探究内孤立波传播过程中的波浪衰减问题，结果如图2所示。

可以看出，波谷传播至X=14 m处的内孤立波形和波幅与理论解基本吻合，但数值模拟的内孤立波具有明显的连续尾波，而且尾波是上凸型和下凹型交替出现。内孤立波由X=14 m传播到X=18 m的过程中波幅值有所衰减，X=18 m的波幅值为0.164 m，与设计波幅相比，衰减了3.5%，衰减幅度不大。同时从图中可以看出内孤立波传播过程中，靠近出口的波面始终在Z=0处，说明消波效果良好。

由图3可以看出，内孤立波速度矢量轨迹类似于一个顺时针的椭圆，以波面为分界面，波面上方与下方的水平流速相反，上方的水平流速与内孤立波传播方向一致，水平向右，而下方的水平流速与内孤立波传播方向相反，水平向左，上层流体的流速大小大于下层流体。内孤立波的垂向速度，以波谷为分界，波谷左侧向上，右侧向下，左右两侧流速大小相当。

2 内孤立波作用下潜水器的运动响应预报 2.1 计算模型

采用的计算模型为Suboff缩比模型，如图4所示。该模型不仅具有围壳，还带有十字尾舵，其主要参数见表1，模型的重量和重心纵向位置与模型在数值水池中所受到的浮力和浮心纵向位置一致。

2.2 计算设置及网格划分

潜水器在内孤立波中的运动是在上述的内孤立波数值水池中模拟，潜水器初始时刻所处位置在X=17～18 m之间。模型在内孤立波作用下的大幅度运动采用重叠网格技术模拟，数值水池各个边界的边界条件和数值造波时保持一致，潜水器的边界条件设置成无滑移壁面。湍流模型采用可实现k-ε模型，采用隐式不定常求解器求解，求解时间步长为0.02 s。模型壁面网格经过细化，保证Y+≈1，总网格数约426万，模型运动区域的网格如图5所示。计算初始阶段对模型施加约束，当内孤立波即将抵达模型首端时释放模型，即t=58 s时。

2.3 计算工况设置

潜水器在内孤立波作用下的运动是由内孤立波诱导的流场以及密度分层引起的静力变化共同引起的。潜水器在水池中潜深不同，则其所处流场和密度场也不同，因此有必要探究潜深对潜水器运动的影响。因此计算4种典型的潜深工况，即潜水器位于淡水层、位于2层流体分层交界面处、位于内孤立波波谷以下。计算过程中，模型处于随波逐流状态，由于内孤立波近似二维波，数值计算中只模拟了纵摇、垂荡和纵荡三自由度运动。潜深d为潜水器圆柱段的水平对称面距两层流体分层交界面的垂直距离，工况设置如表2所示。

2.4 计算结果及分析

1）工况1下潜水器运动响应

由图6（a）垂荡时历曲线可知，内孤立波传播至潜水器时，流场垂向速度向下，潜水器开始下沉，80 s左右，达到下沉最大值−0.185 m，80～86 s之间，潜水器平卧在内孤立波波谷处。通过监测潜水器的运动场景，发现潜水器通过内孤立波的过程中始终平稳坐底在流体分层交界面上方，而内孤立波的波幅值为−0.164 m，与潜水器垂荡值接近。由此可以看出，处于淡水层且靠近流体分层交界面的潜水器在内孤立波作用下的垂荡值由内孤立波的波幅决定。86 s后，潜水器处于波谷左侧，流场垂向速度向上，潜水器开始上浮，90 s之后在尾波作用，潜水器持续上浮；

由图6（b）纵荡时历曲线可知，内孤立波刚传播至潜水器时，流场水平速度向右，潜水器向右漂去，90 s左右纵荡达到最大值0.663 m，此时内孤立波主波基本上经过潜水器。之后在上凸与下凹交替的尾波作用下，潜水器开始向左漂，虽然尾波水平流速在不断变化，但在惯性的影响下，潜水器始终向左漂。

由图6（c）纵摇时历曲线可知，内孤立波经过潜水器时，潜水器开始首倾，然后开始扶正，后又开始首倾至最大艏倾角1.5°。75 s之后，潜水器接近内孤立波波谷，波谷较平缓，波倾角小，潜水器快速扶正，75～86 s，潜水器平卧在波谷附近，纵摇角在0°附近。86 s之后，潜水器运动至波谷左侧，开始抬首，90 s左右潜水器抬首至最大尾倾角3.7°，此时潜水器通过了内孤立波主波，经测量，内孤立波波谷左侧的波倾角约3.5°，波谷右侧的波倾角约2.8°，由此可看出处于淡水层且靠近流体分层交界面的潜水器最大纵摇角与内孤立波主波最大波倾角接近。90 s后，潜水器进入尾波区域，纵摇幅值更大，在93 s出现首倾最大角度4.7°，甚至超过了潜水器经过主波的幅值，这可能是因为尾波波长较短，接近于潜水器的船长，引起了潜水器谐振。

由图7可以看出，处于淡水层且靠近流体分层交界面的潜水器随波逐流状态下作逆时针运动，先向右下运动，到达谷底后，在内孤立波垂向流作用下迅速向上运动，然后受尾波流场作用，向左上运动，最后返回至初始位置的上方。

2）工况2下潜水器运动响应

由图8（a）的垂荡时历曲线可知，内孤立波传播至潜水器时，流场垂向速度向下，同时流体分层交界面向下塌陷，在流体垂向速度和分层密度差的作用下，潜水器开始下沉，78 s时，达到下沉最大值−0.166 m。此时潜水器刚好处于内孤立波波谷处，通过监测潜水器的运动场景，发现潜水器通过内孤立波的过程中始终处于流体分层交界面处，内孤立波78 s时波谷值为−0.164 m，与潜水器垂荡值大致相同。由此可以看出，处于分层交界面处的潜水器在内孤立波作用下的垂荡值由内孤立波波幅决定。78 s后，潜水器处于波谷左侧，流场垂向速度向上，潜水器开始上浮，89 s之后在尾波作用，潜水器在平衡位置处上下振荡。

由图8（b）纵荡时历曲线可知，潜水器大部分处于海水中，内孤立波刚传至潜水器时，海水层水平速度向左，潜水器开始向左漂，89 s左右纵荡值为−0.873 m，此时内孤立波主波刚刚经过潜水器。之后潜水器在交替出现的尾波作用下，先停留一段时间后，继续向左漂去，直到112 s，潜水器达到最大纵荡值−0.969 m，之后潜水器开始向右漂，这可能是因为尾波流场的反复变化，导致潜水器的纵荡运动更加复杂。

由图8（c）纵摇时历曲线可知，由于潜水器在经过内孤立波过程中始终处于分层交界面处，因此潜水器在经过波谷之前是以埋首为主，73 s时达到最大首倾角2.2°。经过波谷之后，潜水器开始抬首，82 s时达到最大尾倾角2.2°，经测量，内孤立波波谷左侧的波倾角约3.5°，波谷右侧的波倾角约2.8°，由此可看出处于两层流体分层交界面处的潜水器的最大纵摇角与内孤立波最大波倾角接近。89 s之后潜水器进入尾波区域，抬艏和埋首变化时间更短，且纵摇幅值也较大，这可能是因为尾波波长较短，接近于潜水器的船长，引起了潜水器谐振。

由图9可以看出，处于2层流体分层交界面处的潜水器随波逐流状态下作顺时针运动，这与Case1中潜水器的运动方向刚好相反，潜水器先向左下运动，到达谷底后，在内孤立波垂向流作用下迅速向左上运动，然后受尾波流场作用，向左上运动，最后返回至分层交界面处。

3）工况3下潜水器运动响应

由图10（a）的垂荡时历曲线可知，内孤立波传播至潜水器时，流场垂向速度向下，潜水器开始下沉，79 s时，达到下沉最大值−0.134 m，此时潜水器刚好经过内孤立波波谷处。通过监测潜水器的运动场景，发现潜水器通过内孤立波的过程中始终处于流体分层交界面下方。79 s后，潜水器处于波谷左侧，流场垂向速度向上，潜水器开始上浮，直至119 s左右，上浮至最大值0.209 m，此后潜水器在紧邻分层交界面的位置振荡。

由图10（b）纵荡时历曲线可知，潜水器刚释放开始向右漂，然后在水平向左流场的作用下，开始向左漂，89 s左右，达到向左的最大纵荡值−0.274 m，此时内孤立波主波刚刚经过潜水器。89 s后潜水器开始在尾波的作用下再次向右漂，122 s左右达到向右的最大纵荡值0.155 m。

由图10（c）纵摇时历曲线可知，潜水器始终在抬首和埋首之间不停振荡，71 s左右达到首倾最大角度2°，76 s左右达到尾倾最大角度2.1°。90 s左右，潜水器通过了内孤立波主波，开始进入尾波区域。之后潜水器在106 s出现尾倾最大角度2°，129 s出现首倾最大角度2.5°，这可能是因为尾波波长较短，接近于潜水器的船长，引起了潜水器谐振。

由图11可以看出，处于内孤立波波谷以下的潜水器随波逐流状态下作顺时针运动，且通过内孤立波后，潜水器运动到2层流体分层交界面附近。

4）不同潜深工况下潜水器运动响应对比

由图12（a）垂荡运动曲线对比可知，初始潜深不同的潜水器在内孤立波作用下的运动差异较大。工况1和工况3下的潜水器通过内孤立波主波的运动过程基本一致，均是随波面起伏。但在完全通过内孤立波之后，工况1下的潜水器上浮至自由面处，而工况2和工况3下的潜水器最终都停留在流体分层交界面附近。

由图12（b）纵荡运动曲线对比可知，工况1和工况3的纵荡运动差异较大，这主要是因为潜水器的潜深不同导致潜水器所处的速度场不同，因此在经历内孤立波的过程中潜水器纵荡运动方向和幅值各有差异。

由图12（c）纵摇运动曲线对比可知，3种工况下，在经过内孤立波主波的过程中，潜水器都是在埋首和抬首之间不断振荡，但埋首和抬首的纵摇幅值均未超过内孤立波的最大波倾角，且最大纵摇幅值出现的位置并不一定是在内孤立波主波经过潜水器的时候，也可能发生在内孤立波的尾波区域内，这说明内孤立波的尾波对潜水器的运动有相当大的影响，不可忽视，这可能是因为内孤立波的尾波波长较短，与潜水器的艇长接近，造成了潜水器谐振。

图13为3种工况下潜水器在内孤立波作用下的运动轨迹。可以看出，工况2和工况3的潜水器在内孤立波作用下运动特点相似，均是顺时针旋转，最终均运动到分层交界面附近，而工况1的潜水器运动就刚好相反，是逆时针运动，最终快速上浮至自由面附近。

3 结　语

采用CFD方法建立了内孤立波数值池，对Suboff模型受内孤立波作用的运动响应进行了数值模拟，分析了不同潜深下潜水器在内孤立波作用下的运动响应规律，阐明了内孤立波水动力和分层密度差在潜水器运动响应过程中所起的作用，比较了不同潜深下运动响应特性的差异，得出以下结论：

1）CFD方法可以有效模拟内孤立波，所造波形与eKdV理论解吻合，误差较小，且数值模拟的内孤立波具有明显的连续尾波。

2）模型处于分层交界面时，受分层密度差静力作用，模型“漂浮”在内界面上，并随内孤立波波面起伏，呈“随波逐流”状态，模型在内孤立波作用下的垂荡值与内孤立波的波幅基本一致，通过内孤立波时的纵摇角由内孤立波的波倾角决定；模型处于分层交界面上或分层交界面下时，受内孤立波流场水动力作用，模型也产生了明显的垂荡、纵荡和纵摇。

3）模型通过内孤立波时，最大纵摇角不一定出现在经过主波期间，也可能出现在经过内孤立波尾波期间，因此在探究内孤立波对潜水器运动响应尤其是纵摇运动的影响时，不能只考虑内孤立波主波，尾波也要考虑。

4）处于分层交界面上方的模型在内孤立波作用下作逆时针运动，最终会持续上浮至自由面附近；而处于分层交界面处或分层交界面下方的模型运动轨迹刚好相反，为顺时针运动，且最终均会运动至分层交界面附近。

本文仅对不同潜深下模型受内孤立波作用的运动响应特性进行了分析，后续将继续开展航速、内孤立波波幅、潜水器的水下稳性高、应急上浮等操控策略对潜水器在内孤立波作用下运动的影响分析，并开展试验验证研究。