内孤立波是一种发源于海洋密度跃层中的非线性大振幅波动，它通常由强潮流经过地形陡峭的海底陆坡、海脊时激发产生。在传播过程中，内孤立波携带着巨大能量，能将能量和动量从上层流体传至值较低的下层流体，是能量和动量垂向传输的重要载体^[1]。尤其是强非线性内孤立波，甚至会对穿越波面的水下潜器产生巨大的剪切作用，引起海水强烈的幅聚、幅散和突发性的强流，极易导致海洋工程或者海底潜器出现疲劳甚至损坏^[2-4]。

目前，针对结构物在内孤立波在流场中载荷特性的研究主要有模型试验、理论分析以及CFD数值模拟这三种方法。众多学者设计了内孤立波实验水槽，通过实验手段研究了内孤立波作用下结构的载荷特性。Wang等^[5-7]在大型实验水槽中开展了内孤立波流场中固定式水下细长体的受力特性研究；谷梦梦等^[8]和黄文昊等^[9]用大型重力式分层流试验水槽研究了内孤立波对细长潜体和圆柱结构的载荷特性；Wei等^[10]和Du等^[11]在分层流水池中对比研究了细长圆柱体和潜器与内孤立波的相互作用；Wang等^[12]在考虑密度连续分层的内孤立波水槽中开展了水下结构物内孤立波载荷实验研究，获取了内孤立波作用下不同浸没深度圆桩的水动力系数。

此外，大量学者基于内孤立波理论与Morison公式发展了内孤立波场作用下结构载荷计算的解析方法，能够实现在理想条件下计算规则结构所受到的内孤立波流场载荷。Cai等^[13-14]将内孤立波理论与Morison公式相结合，提出了一种简单细长结构物内孤立波载荷的计算方法，发现内孤立波诱导产生的强剪切流将对水面、水下结构物的受力产生显著影响；文献[15-19]的研究中均采用Morison公式计算垂直圆柱结构和潜器表面的内孤立波载荷；刘忻丹等^[20]利用Morison方程和Froude-Krylov计算潜器受到的内孤立波载荷，并结合三自由度运动方程，求解了潜器的运动姿态和轨迹。

总体而言，实验方法成本较高并且存在尺度效应问题(模型与实物之间相似条件不能全部满足而造成的换算误差)。解析方法有一定的局限性，忽略了结构物对内孤立波流场的干扰效应，并且很难捕获内孤立波与复杂结构形式之间的耦合效应。因此，近年来计算流体力学(Computational fluid mechanics, CFD)方法多凭借其直观、可重复性好的特点，且已逐渐成为研究内孤立波与结构物耦合作用问题的重要手段。基于两层流体间断分层模型，汪逢港^[21]和刘倩等^[22-23]研究了内孤立波流场分布和立管所受到的内孤立波载荷；付东明等^[24]和汪超等^[25-29]对内孤立波作用下潜器的受力载荷和运动响应进行了一系列的研究。对于固定潜器，他们研究发现，当潜器位于内孤立波界面上、下时，其所受的水平力不可忽视。潜器在穿越波面时，垂向力和水平力都会在短时间内发生巨变；李国迎^[30]和关晖等^[31]研究了不同内孤立波类型下潜器的载荷特性。基于流体连续分层模型，杨永春等^[32]和Li等^[33]研究了内孤立波的生成及其与结构物的相互作用。王展等^[34]对比研究了连续分层和间断分层内孤立波模型各自的适用范围。Liu等^[35]和杨帆^[36]分析了在不同潜深、流体分层比以及不同内孤立波波幅作用下固定潜器的载荷特性，认为潜器潜深是影响内孤立波作用下潜器受力载荷的关键性因素，尤其是对于能否穿越内孤立波波面的潜器，其垂向力载荷存在较大的区别。本研究团队以往研究了自航潜器在内孤立波流场中的流固耦合效应，重点分析了高航速潜器在内孤立波流场中大幅抬艏的运动现象^[37-39]。

尽管前人已经开展了大量的数值研究来分析结构物在内孤立波中的受力特性，但对于潜器在内孤立波流场中的力学机理及其受力突变的原因，目前仍未完全揭示。鉴于此，本文以固定潜器在内孤立波中受力特性这一典型工况为例，系统研究潜器潜深、内孤立波波幅、流体分层比和内孤立波类型等多种参数对固定潜器载荷特性的影响规律，寻找影响其受力突变的关键因素，以解释潜器在内孤立波流场中失控和掉深的原因，旨在为潜器在内孤立波流场中安全航行提供理论与技术支撑。

1 数值方法

内孤立波流场可被视为不可压缩黏性流体，流场的控制方程为质量连续性方程和雷诺平均化的Navier-Stokes (RANS)方程^[40]

$ \frac{\partial \bar{u}_{i}}{\partial x_{i}}=0, $

(1)

$ \begin{gather*} \frac{\partial \bar{u}_{i}}{\partial t}+\sum\limits_{j} \bar{u}_{j} \frac{\partial \bar{u}_{i}}{\partial x_{j}}=-\frac{1}{\rho} \frac{\partial P}{\partial x_{i}}+F_{i}+ \\ v \sum\limits_{j} \frac{\partial^{2} \bar{u}_{i}}{\partial x_{i}^{2}}-\sum\limits_j \frac{\partial \overline{u_i^{\prime} \bar{u}_j^{\prime}}}{\partial x_j} 。\end{gather*} $

(2)

式中：$\bar{u}_{i}$ 为平均速度的第$i$ 个分量（$i=1, 2, 3$，分别对应$x、y、z$ 方向）；$x_{i}$ 为空间直角坐标系的第$i$ 个坐标分量；$x_{j}$ 为空间直角坐标系的第$j$ 个坐标分量；$P$ 为压力；$F_{i}$ 为质量力；$v$ 为运动黏滞系数；$t$ 为时间，表征流场的时间演化；$\overline{u_{i}^{\prime} u_{j}^{\prime}}$ 为雷诺应力项；$u_{i}^{\prime}$ 和$u_{j}^{\prime}$ 分别为脉动速度的第$i$ 个和第$j$ 个分量。

为封闭求解控制方程(式(2))，还需引入湍流模型。RNG k-ε湍流模型^[41]在计算复杂几何问题中能展现出较好的适应性与计算收敛性。考虑到内孤立波波面流场附近的流线具有较大的曲率，并且潜器结构表面形状复杂，本文选取RNG k-ε模型

$ \begin{gather*} \frac{\partial}{\partial t}(\rho k)+\frac{\partial}{\partial x_{i}}\left(\rho k u_{i}\right)=\frac{\partial}{\partial x_{j}}\left(\alpha_{k} \mu_{\text {eff }} \frac{\partial k}{\partial x_{j}}\right)+ \\ G_{k}+G_{b}-\rho \varepsilon-Y_{M}+S_{k} , \end{gather*} $

(3)

$ \begin{gather*} \frac{\partial}{\partial t}(\rho \varepsilon)+\frac{\partial}{\partial x_{i}}\left(\rho \varepsilon u_{i}\right)=\frac{\partial}{\partial x_{j}}\left(\alpha_{\varepsilon} \mu_{\text {eff }} \frac{\partial \varepsilon}{\partial x_{j}}\right)+ \\ C_{\varepsilon 1} \frac{\varepsilon}{k}\left(G_{k}+C_{\varepsilon 3} G_{b}\right)-C_{\varepsilon 2} \rho \frac{\varepsilon^{2}}{k}-R_{\varepsilon}+S_{\varepsilon} 。\end{gather*} $

(4)

式中：ρ为流体密度；k和ε分别为湍流动能和湍流耗散率；u_i为流体平均速度；G_k为平均速度梯度产生的湍流动能；G_b为浮力b产生的湍流动能；Y_M为可压缩湍流波动膨胀对整体耗散率的贡献；μ_eff为有效黏滞系数；α_k和α_ε分别为k和ε的有效普朗特数的倒数；C_ε1、C_ε2和C_ε3均为耗散方程的经验参数；S_k与S_ε分别为用户自定义的k和ε的源项。R_ε的表达式为

$ R_{\varepsilon}=\frac{C_{\mu} \rho \eta^{3}\left(1-\eta / \eta_{0}\right)}{1+\beta \eta^{3}} \times \frac{\varepsilon^{2}}{k} 。$

(5)

式中：C_μ为模型修正系数；η为无量纲应变或平均流时间尺度与湍流时间尺度的比值；η₀为应变强度的临界常数(η₀=4.38)，核心作用是界定η的物理合理范围，进而保证模型修正项R_ε的有效性；β为修正项的权重调节常数(β=0.012)，其核心作用是量化η³对湍流耗散修正的影响权重^[37]。

本文选择eKdV理论描述内孤立波传播特征。eKdV理论在KdV理论基础上引入了一个高阶非线性项，可以模拟中等非线性大振幅的内孤立波。eKdV理论中，内孤立波界面位移ζ满足

$ \frac{\partial \zeta}{\partial t}+c_{0} \frac{\partial \zeta}{\partial x}+c_{1} \zeta \frac{\partial \zeta}{\partial x}+c_{2} \frac{\partial^{3} \zeta}{\partial x^{3}}+c_{3} \zeta^{2} \frac{\partial \zeta}{\partial x}=0 。$

(6)

其中

$ \begin{gathered} c_{0}^{2}=g h_{1} h_{2}\left(\rho_{2}-\rho_{1}\right) /\left(\rho_{1} h_{2}+\rho_{2} h_{1}\right), \\ c_{1}=-3 c_{0}\left(\rho_{1} h_{2}^{2}-\rho_{2} h_{1}^{2}\right) /\left[2\left(\rho_{1} h_{1} h_{2}^{2}+\rho_{2} h_{1}^{2} h_{2}\right)\right], \\ c_{2}=c_{0}\left(\rho_{2} h_{1} h_{2}^{2}+\rho_{1} h_{1}^{2} h_{2}\right) /\left[6\left(\rho_{1} h_{2}+\rho_{2} h_{1}\right)\right], \\ c_{3}=\frac{3 c_{0}}{h_{1}^{2} h_{2}^{2}}\left(\frac{7}{8}\left(\frac{\rho_{1} h_{2}^{2}-\rho_{2} h_{1}^{2}}{\rho_{1} h_{2}+\rho_{2} h_{1}}\right)^{2}-\frac{\rho_{1} h_{2}^{3}+\rho_{2} h_{1}^{3}}{\rho_{1} h_{2}+\rho_{2} h_{1}}\right) 。\end{gathered} $

(7)

式中：c₀、c₁、c₂、c₃均为中间参数；ρ₁和ρ₂分别为上、下层流体的密度；h₁和h₂分别为上、下层流体的厚度。由式(7)可以解得

$ \zeta(x, t)=\zeta_{0} /\left[B+(1-B) \cosh ^{2}\left(\left(x-c_{\mathrm{eKdV}} t\right) / \lambda_{\mathrm{eKdV}}\right)\right] \text { 。} $

(8)

式中：$B=-c_{3} \zeta_{0} /\left(2 c_{1}+c_{3} \zeta_{0}\right); \zeta_{0}$ 为内孤立波波幅；$c_{\mathrm{eKdV}}$ 与$\lambda_{\mathrm{eKdV}}$分别为内孤立波传播相速度及特征波长，二者表达式分别为

$ \begin{gather*} c_{\mathrm{eKdV}}=c_{0}+\zeta_{0}\left(c_{1}+c_{3} \zeta_{0} / 2\right) / 3 , \\ \lambda_{\mathrm{eKdV}}=\sqrt{12 c_{2} /\left[\left(c_{1}+c_{3} \zeta_{0} / 2\right) \zeta_{0}\right]} 。\end{gather*} $

(9)

本文采用速度入口法造波，在速度入口边界给定波面速度方程，速度入口处上、下层流体速度设为u₁和u₂，二者表达式分别为

$ \begin{align*} & u_{1}=-c_{\mathrm{eKdV}} \frac{\zeta_{0}}{h_{1}}, \\ & u_{2}=c_{\mathrm{eKdV}} \frac{\zeta_{0}}{h_{2}} 。\end{align*} $

(10)

为捕捉两层流体交界面的时历变化，采用流体体积(Volume of fluid，VoF)模型来捕捉界面，从而确定内孤立波界面的位置。VoF模型通过计算网格单元中流体所占网格体积的比例，即体积分数函数α_j实现对波面的跟踪捕捉。当α_j=0时，表示该单元内没有指定的第j相流体；当α_j=1时，则表示该单元内全部是第j相流体；当0 < α_j < 1时，表示该单元内同时占有两相流体。在每一个控制体中，各相的体积分数之和为1，则

$ \sum\limits_{j=1}^{2} \alpha_{j}=1。$

(11)

通过求解体积分数连续性方程捕捉两相界面，第j相的体积分数应满足连续性方程

$ \frac{\partial \alpha_{j}}{\partial t}+\frac{\partial\left(\alpha_{j} u_{i}\right)}{\partial x_{i}}=0 。$

(12)

式中i=1, 2。

2 模型的建立

本文以中国南海海洋环境为研究背景，选取1∶20的缩尺比，利用ANSYS Fluent软件建立数值模型，将数值水槽的长×宽×高设置为150 m×15 m×20 m，上、下流体层厚度分别设置为h₁ = 5 m和h₂ = 15 m，上、下流体层密度分别设置为ρ₁ = 998 kg/m³和ρ₂=1 025 kg/m³。内孤立波数值水槽的顶部依据“刚盖近似”定义为对称边界，数值水槽左、右侧边界分别定义为速度入口边界以及压力出口边界，数值水槽的底面被定义无滑移壁面，前、后边界均被设置为对称边界，如图 1所示。

本文采用美国国防高级研究计划局Suboff标准潜器模型作为研究对象，如图 2所示，潜器的结构由前体、平行中体、后体以及附体四部分组成，其中附体包括指挥塔与尾翼。潜器模型的主尺度如表 1所示。

在内孤立波数值水槽中引入潜器模型，建立内孤立波流场中的固定潜器模型，如图 3所示，潜器位置距速度入口为50 m。d为初始时刻潜器与内孤立波界面的相对距离，即初始潜深。

采用重叠网格技术计算潜器在内孤立波流场中的受力特性，整个计算域网格分为前景网格与背景网格两部分，二者通过构建重叠交界耦合面，采用插值的方式实现数据传递。背景网格仅用于数值造波，在波面附近进行网格的垂向加密，使局部网格尺寸为0.2 m，采用六面体结构化网格离散背景网格，背景网格单元总数量约为131万；前景网格主要用于描述潜器的受力情况，采用影响体(Body of influence，BoI)法对潜器附近的流体域进行局部加密，能够更准确地捕捉潜艇周围流场特征，采用多面体非结构化网格离散前景网格，前景网格单元总数量约为80万，两者组装后，网格单元总数量为211万，如图 4所示。

在本人团队以往工作^[37-39]中，已将数值结果与文献[7]中固定圆柱在内孤立波流场中的实验结果进行了对比验证。对比内容包括流场特征监测点的水平速度和垂直速度，以及水下圆柱体的水动力系数C_x和C_z时历曲线，数值结果和实验结果吻合良好。从而证明本文数值模型能够有效模拟内孤立波的传播过程，并且能够准确地计算内孤立波环境下水下潜器的水动力载荷，验证了模型的可靠性，故本文不再重复验证。本文将选取网格尺寸为0.02 m，时间步长Δt = 0.02 s。

3 内孤立波作用下固定潜器载荷特性研究 3.1 潜深对内孤立波作用下潜器受力特性的影响

本文以水槽深度H、位于下层流体的潜器所受重力(ρ₂gV)与其恢复力矩(ρ₂gVL_pp)分别为长度、力与力矩的特征量，对相关物理量进行无量纲变换，具体表达式如下：

$ \begin{aligned} & X^*=\frac{X}{H}, Y^*=\frac{Y}{H}, Z^*=\frac{Z}{H}, \zeta_0^*=\frac{\zeta_0}{H}, d^*=\frac{d}{H}, \\ & F_x^*=\frac{\boldsymbol{F}_x}{\rho_2 \boldsymbol{g} V}, F_y^*=\frac{\boldsymbol{F}_y}{\rho_2 \boldsymbol{g} V}, F_z^*=\frac{\boldsymbol{F}_z}{\rho_2 \boldsymbol{g} V}, \\ & M_x^*=\frac{\boldsymbol{M}_x}{\rho_2 \boldsymbol{g} V L_{\mathrm{pp}}}, M_y^*=\frac{\boldsymbol{M}_y}{\rho_2 \boldsymbol{g} V L_{\mathrm{PP}}}, M_z^*=\frac{\boldsymbol{M}_z}{\rho_2 \boldsymbol{g} V L_{\mathrm{pp}}} 。\end{aligned} $

(13)

式中：上标星号的物理量为无量纲物理量；X、Y和Z分别表示潜器的纵荡、横荡和垂荡位移；d代表初始时刻潜器距离内孤立波界面的相对距离；ζ₀表示内孤立波的振幅；F_x、F_y和F_z分别为潜器所受到的x、y和z方向的作用力；M_x, M_y和M_z分别为潜器在横摇、纵摇和艏摇方向所受到的力矩。

为探究初始潜深对固定潜器受力特性的影响，本文设置8个算例，分别模拟位于上层流体(d^* =0.1, 0.05)、流体分界面(d^* =0)、下层流体(d^* =-0.05, -0.1, -0.15, -0.25, -0.35)的固定潜器在相同振幅(ζ₀^* = -0.15)时内孤立波作用下的受力，如表 2所示。图 5为固定潜器不同潜深位置示意图。为保证潜器浮力与重力平衡，针对潜器的密度进行修改，使其与所在位置处的初始流场流体密度保持一致。

本文重点关注固定潜器位于垂直平面的载荷变化情况，即无量纲水平载荷F_x^*、无量纲垂直载荷F_z^*、无量纲纵摇力矩M_y^*。

图 6(a)展示了不同潜深潜器的水平载荷时历存在差异，依据其变化特征可分为以下三类。第一类是潜器位于内孤立波界面及其上方的工况(d^* ≥ 0)。潜器在内孤立波传播过程中始终位于分界面或上层流体，在上层流体场的诱导下潜器的无量纲水平载荷F_x^*先增加到正幅值后再减小到负幅值，随着内孤立波远离，其无量纲水平载荷F_x^*缓慢减小到零附近，并伴随着尾波作用小幅波动。第二类是潜器位于分界面以下且穿过内孤立波波面的工况(ζ₀^*≤d^*＜0)。在内孤立波传播过程中，潜器从下层流体穿过右侧波面短暂地来到上层流体，最终再次穿过波面回到下层流体，如图 7(b)和(d)所示，其水平载荷F_x^*在右侧波面的诱导下先小幅波动后增加到正幅值，由于内孤立波流场为一个椭圆形双剪切流场，穿过波面来到上层流体的潜器水平载荷持续减小到负幅值，随后缓慢趋于0，而后再次穿过波面，并在左侧波面流场影响下载荷再次增加到正幅值。距分界面越远的潜器，其水平载荷正幅值越小，负幅值越大，且到达第一个负幅值的时间越晚。第三类是潜器位于分界面以下且不穿过内孤立波波面的工况(d^* < ζ₀^*)。由于内孤立波上、下层流体场水平流速相反，位于下层流体且不穿过波面潜器的水平载荷时历(d^* < ζ₀^* < 0)正好与第一类工况(潜器位于分界面及其以上的工况，d^* ≥ 0)的载荷时历变化相反。在下层流体场的诱导下，潜器水平载荷先减小到负幅值又增加到正幅值，随着波面远离，其载荷慢慢恢复到零。

水下潜器在初始时刻处于静力平衡状态，随着内孤立波的传播，潜器所处环境流场的密度将发生改变，潜器与流场之间密度差是导致潜器无量纲垂向载荷F_z^*异变的关键因素。与无量纲水平载荷F_x^*类似，将图 6(b)中不同潜深潜器的垂向载荷时变特性也分为三类。第一类是潜器位于内孤立波界面以上的工况(d^*>0)。潜器位于上层流体且不穿过波面，潜器与周围流场之间的密度差不会发生明显变化，其垂向载荷随流场的传播只会产生小幅波动。第二类是潜器位于分界面及其以下且穿过波面的工况(ζ₀^*≤d^*≤ 0)。在此工况下的内孤立波传播过程中，潜器所在的密度场会发生变化。上、下层流体的密度差会导致潜器无量纲垂向载荷F_z^*产生大幅度波动，其变化与潜器进入上层流体的体积、在上层流体滞留时间以及密度场变化等因素都存在紧密联系。对于穿过波面时正好位于波谷的潜器(d^* = -0.15)，由于其进入上层流体的体积有限，其垂向载荷幅值远小于能够完全浸没于上层流体的潜器(d^*= -0.05和d^*= -0.1)。当潜器完全穿过波面(d^*= -0.05和d^*= -0.1)时，垂向载荷幅值基本相等，且无量纲垂向载荷F_z^*的两次幅值跃迁时间间隔随着潜器在上层流体场滞留时间的增加而增加。初始时刻位于流体分界面潜器(d^*= 0)的垂向载荷变化随波面波动，随着内孤立波波面的接近，潜器深入到上层流体，潜器垂向载荷减小，其幅值小于位于下层流体且能够穿过波面完全浸没在上层流体潜器(d^*= -0.05和d^*= -0.1)的垂向载荷幅值。特别地，当内孤立波远离时，潜器部分体积深入到下层流体，由于潜器密度小于下层流体场密度，其垂向载荷持续增加到零以上。第三类是潜器位于分界面以下且不穿过内孤立波波面的工况(d^* < ζ₀^*)。由于潜器没有受到上层流体密度场的干扰，潜器受到的垂向载荷只产生小幅波动。对于d^*= -0.35的工况，潜器垂向载荷的幅值约为d^*=-0.05和d^*=-0.1时工况的3%。

图 6(c)展示了固定潜器无量纲纵摇力矩M_y^*的时变特性，内孤立波载荷在潜器表面的分布情况决定了固定潜器无量纲纵摇力矩M_y^*的时变规律，具体可以分为以下三种类型。第一类是潜器位于分界面及其以上的工况(d^* ≥ 0)。对于d^*= 0.1的工况，潜器纵摇力矩M_y^*先在右侧波面上层流场向下速度的诱导作用下减小到最小值，在左侧波面向上流场速度诱导作用下逐渐增加到最大值，随着内孤立波的远离，力矩逐渐减小到零附近。对于d^*= 0的工况，潜器位于分界面处，受界面影响无量纲纵摇力矩M_y^*波动比较剧烈。第二类是潜器位于分界面以下且穿过内孤立波波面的工况(ζ₀^*≤d^*<0)。对于d^*=-0.05、d^*=-0.1和d^*=-0.15的工况，潜器穿过波面时垂向力会出现大幅值跃迁，导致潜器纵摇力矩时历出现较为显著的波谷与波峰。与水平载荷类似，距分界面越远的潜器，第一个波谷到达的时间越晚。第三类是潜器位于下层流体且不穿过波面的工况($d^* < \zeta_0^* < 0$)。对于d^*=-0.3的工况，由于潜器未受到上层流场的干扰，其纵摇力矩变化规律与第一类(位于上层流体及其分界面潜器)正好相反，即在右侧波面下层流场的作用下增加到最大值，随后在左侧波面附近下层流场的诱导作用下减小到最小值，随着内孤立波的远离，力矩逐渐减小到零附近。

本研究发现，潜器受到内孤立波载荷的作用，与此同时，内孤立波经过潜器附近会产生绕射现象，其波面会产生波动。当内孤立波远离潜器时，波动效果会更为显著，如图 7(d)所示。并且随着潜器与波面距离减小，潜器对波面的干扰效应愈加显著。由于本文潜器的特征长度相比于内孤立波波长要小一个量级，潜器运动对流场干扰比较有限，在实际物理背景下，潜器对流场的干扰效应可以被忽视。但在实验中，其结构物的尺度被过分高估，导致模拟结果存在一定的误差。

3.2 波幅对内孤立波作用下潜器受力特性的影响

内孤立波的振幅是影响潜器受力特性的重要参数，内孤立波振幅越大，诱导产生的流场速度越高，设置4种内孤立波振幅($\xi_0^*$=-0.05, -0.1, -0.15, -0.2)的工况(对应实际海洋中内孤立波波幅为20、40、60和80 m)，探讨其对固定潜器(d^* = -0.05)受力特性的影响，如图 8所示。

图 9展示了位于下层流体(d^* = -0.05)固定潜器在不同波幅内孤立波作用下的载荷时历：无量纲水平载荷F_x^*随内孤立波波幅增大而增大，且到达水平载荷幅值的时间随着内孤立波波幅的增加而提前。同样，无量纲纵摇力矩M_y^*的负幅值也随着内孤立波波幅的增加而增大。高幅值的内孤立波($\xi_0^*$ = -0.2)诱导产生的无量纲纵摇力矩M_y^*正幅值大于低幅值内孤立波(ξ₀^*=-0.05)作用下的幅值。且随着内孤立波波幅增加，内孤立波波面非线性增加，其产生的尾波序列愈加显著，尾波诱导的固定潜器载荷幅值波动也愈加明显。

影响潜器无量纲垂向载荷F_z^*变化的主控因素是潜器进入上层流体的体积，内孤立波发挥的作用并不显著。当内孤立波波幅与潜器初始潜深相等($\xi_0^*$= d^* =-0.05)时，内孤立波通过时固定潜器恰好位于波谷，潜器仅有部分体积进入上层流体，其垂向载荷F_z^*幅值明显小于其他工况($\xi_0^*$=-0.1, -0.15, -0.2)。当潜器完全浸没于上层流体时，内孤立波幅值的变化几乎不改变潜器垂向载荷幅值。值得注意的是，随着波幅的增加，潜器穿过波面滞留在上层流体时间变长，其垂向载荷出现的两次峰值的时间间距变长。且波面传播相速度随着波幅的增加而增加，固定潜器垂向载荷的第一次峰值发生时间随波幅的增加而提前。

3.3 流体分层比对内孤立波作用下潜器受力特性的影响

不同上、下流体厚度分层比的内孤立波流场流速存在差异，本研究设置3种上、下层厚度比(h₁∶h₂=1∶4, 1∶3, 3∶7)来模拟实际海洋中的分层情况，对应实际海洋中的上、下层厚度组合分别为80 m & 320 m、100 m & 300 m、120 m & 280 m，以探讨其对固定潜器(d^* = -0.05)受力特性的影响。

图 10展示了位于下层流体(d^*= -0.05)的固定潜器在不同分层比内孤立波流场作用下的受力时历曲线。分析可知，当分层比不同时，受力曲线变化趋势基本一致：因潜器两次穿过波面，其水平载荷时历曲线存在显著的波峰与波谷各两个；随着内孤立波接近潜器来到上层流体，其垂直载荷迅速减小并维持一段时间的负幅值，待其回到下层流体后，垂直载荷增大并趋于零。无量纲纵摇力矩$M_y^*$随垂直载荷的变化而变化，其时历曲线存在显著的波峰与波谷各一个。

随着上、下层流体分层比增大，内孤立波波长增大，流场速度反而减小。因此，固定潜器无量纲水平载荷F_x^*以及无量纲纵摇力矩M_y^*的幅值随着流体分层比的增大而减小，各个方向上到达幅值的时间随着流体分层比的增大而相对延后。值得注意的是，无量纲垂向载荷F_z^*幅值几乎不受流体分层比的影响，但随着分层比的增加，潜器滞留在上层流体的时间也增加，潜器垂向载荷维持负幅值的时间也会变长。且由于内孤立波引起的流速随着流体分层比的增大而降低，固定潜器垂向载荷峰值会略微减小。由此可以推测，当潜器位于其他位置且不受上、下层流体密度干扰时，分层比对垂向载荷的影响会愈加显著。

3.4 内孤立波类型对潜器受力特性的影响

在实际海洋中存在上凸形与下凹形两种形态的内孤立波，为探讨内孤立波类型对固定潜器载荷的影响，将针对上凹型或下凹形内孤立波作用下不同位置处(d^* = 0.05, -0.05)潜器的载荷特性进行分析。下凹形内孤立波流体分层比h₁∶h₂ =1∶3，内孤立波波幅$\xi_0^*$= -0.15，上凸形内孤立波流体分层比h₁∶h₂ =3∶1，内孤立波波幅ξ^*₀= 0.15。

如图 11所示，由于上凸形内孤立波上层流体和下凹形内孤立波下层流体在流速大小和方向上几乎完全相同，同时上凸形内孤立波下层流体和下凹形内孤立波上层流体也具有高度一致的流速特征，因此，位于上凸形内孤立波流场的上层流体固定潜器(d^* = 0.05)与位于下凹形内孤立波流场的下层流体固定潜器(d^* =-0.05)的水平载荷时历几乎完全重合。同样的，位于上凸形内孤立波流场的下层流体固定潜器(d ^*=-0.05)与位于下凹形内孤立波流场的上层流体固定潜器(d^* = 0.05)的水平载荷时历也几乎同步。

无量纲垂向载荷F_z^*时历曲线如图 12 (b)所示，位于上凸形内孤立波流场的上层流体固定潜器(d^*=0.05)与位于下凹形内孤立波流场的下层流体固定潜器(d^* = -0.05)均2次穿过波面。不同的是，上凸形内孤立波流场是从上层流体穿过波面来到下层流体，然后再穿过波面，并最终回到上层流体，如图 12(a)—(f)，因此其垂向载荷存在明显的正峰值。下凹形内孤立波流场是从下层流体穿过波面来到上层流体，然后再穿过波面回到下层流体，如图 7(a)—(e)，其垂向载荷存在明显的负峰值，两者的时历曲线几乎沿着F_z^*= 0镜像对称。当位于上层或下层流体的固定潜器不穿过波面时，两者垂向载荷仅在流场影响下小幅波动。由于在上凸形与下凹形内孤立波流场作用下，不同流体层固定潜器垂直载荷时历特征恰好相反，因此位于上凸形与下凹形内孤立波流场作用下不同流体层固定潜器纵摇力矩时历也几乎沿着M_y^*= 0镜像对称。

4 结论

本文基于CFD方法建立了内孤立波数值水槽，研究了内孤立波作用下固定潜器的受力特性，并分析讨论了潜器潜深位置、内孤立波波幅、内孤立波类型和上、下流体厚度分层比等参数对固定潜器载荷特性的影响规律，得出以下结论：

(1) 潜器相对于内孤立波的潜深位置是影响其载荷的关键参数，当潜器能穿过内孤立波时，距分界面越远的潜器，其受到的无量纲水平载荷F_x^*的正幅值越小，负幅值越大，且到达第一个负幅值的时间越晚，其无量纲垂向载荷F_z^*和无量纲纵摇力矩M_y^*会出现大幅波动，与不穿过波面的潜器形成鲜明对比。尤其是位于分界面以下且不穿过波面的潜器，其垂向载荷幅值仅为穿过波面且能够完全浸没在上层流体的潜器垂向载荷幅值的3%。

(2) 内孤立波波幅是影响固定潜器受力特性的重要参数，波幅越大，诱导产生的流场速度越高，传播相速度就越大。但当潜器能够完全进入内孤立波时，其无量纲垂向载荷F_z^*幅值受波幅变化的影响并不明显，但载荷作用时长存在明显差异。

(3) 上、下层流体分层比也是影响潜器载荷变化的重要参数，随着上、下层流体分层比增加，内孤立波波长增加，流场速度反而减小。固定潜器无量纲水平载荷F_x^*以及无量纲纵摇力矩M_y^*的幅值随着流体分层比的增大而减小，到达幅值的时间也随着流体分层比的增大而相对延后。无量纲垂向载荷F_z^*幅值几乎不受流体分层比的影响，随着分层比增加，潜器滞留在上层流体的时间增加，潜器垂向载荷维持负幅值的时间也会变长，但其载荷幅值并没有明显的变化。

(4) 在流体分界面镜面对称位置处，上凸形内孤立波流场的流速大小、方向几乎与下凹形的完全相同，因此二者的固定潜器水平载荷几乎完全重合。与此同时，由于上、下流场密度场在两种波型中变化趋势恰好相反，导致固定潜器的无量纲垂向载荷F_z^*和无量纲纵摇力矩M_y^*时历曲线几乎沿着水平零线镜像对称。