0 引　言

在海洋中，由于盐度和温度的垂向差异，导致海水密度垂向分布不均匀，最终在海洋内部形成稳定连续分布的密度跃层。且相对于海水与空气接触的自由面上下密度差来说，海水内部密度跃层的密度差极小，相当于将分层介质置于微重力场中，其恢复力极小，仅为水面波的千分之一^[1]。因此极小的扰动也会在海水内部产生较大振幅的波动，从而在密度跃层附近产生内波。

内波多由强潮流越过陡峭地形激发，具有很强的随机性，周期在几分钟到几十小时，波长可从几米到几百千米，振幅可达几百米。且由于内波通常携带巨大能量，在传播过程中常会引起海水层结剧烈波动，使密度跃层上下流速相反，导致海水的辐聚和突发性强流^[2]，对水下航行体的航行安全有巨大威胁。

海洋内波对水下航行体的水动力影响，国内开展了许多研究。其中针对内波的数值水槽模拟主要基于两层流体假设下的内孤立波理论解，包括KdV（Korteweg-de Vries）、mKdV（Modified KdV）、eKdV（Extended KdV）、MCC（Miyata-Choi-Camassa）理论等^[3]。其中，KdV理论主要适用于浅水弱非线性、小振幅内波。付东明等^[4]基于双推板的内孤立波数值造波方法，数值模拟了与KdV理论解吻合良好的数值造波结果，并计算了实尺度光体SUBOFF模型约束状态下遭遇内波时的受力载荷。黄苗苗等^[5]基于KdV理论，采用速度入口造波法，分析研究了约束状态SUBOFF模型遭遇内波时的受力载荷和流场变化。但随着内波振幅的增大，KdV理论不再适用，mKdV理论和eKdV理论通过引入高阶非线性项对KdV方程进行修正，对稍大振幅内波进行模拟，但存在极限波幅。姚志崇等^[6]基于二层流体的eKdV理论进行内波模拟，对水下航行体遭遇内波受力变化机制进行分析。发现速度与密度场分布对垂向力起决定作用，水平力主要受水平速度影响。且会出现俯仰力矩的极值现象。汪超等^[7]基于强分层假设下的内波mKdV理论，结合速度入口造波和重叠网格方法，模拟不同内波环境下内孤立波与固定、悬浮航行体耦合作用过程，通过数值求解得到水下航行体的水动力载荷变化规律。并用实验结果验证了数值模拟的准确性。MCC理论可以很好地模拟强非线性、大振幅内波。高雪原等^[8]基于2层流体模型和 MCC 方程，采用 Fluent 研究了内孤立波作用下深水中立管的荷载变化。

而上述研究都基于2层流体假设，用2层不同密度的均质流体来表征密度跃层，忽略了密度跃层厚度的影响。实际海洋中浅水区密度跃层厚度在10～20 m，深水区密度跃层厚度最大可达150 m^[9]，无法忽略。且2层流体假设下，内波交界面处会出现水平速度、密度突变，在水下航行体遭遇内波时的受力载荷计算中会存在一定误差，与实际情况有所差距。李景远等^[10]基于DJL（Dubriel-Jacobin-Long）方程，采用初始场造波方法构建了连续密度分层的内孤立波，与eKdV理论对比，发现忽略密度跃层的厚度会高估内孤立波的传播速度。上述研究均基于理论解进行造波，没有引入地形变化、背景潮流对内波环境场的影响。而麻省理工非静压海洋环流模型MITgcm（Massachusetts Institute of Technology General Circulation Model）可以模拟内潮与地形相互作用生成内波的全过程以及陆架等地形变化对内波的影响，能构建大范围较实际环境场。且该模型可得到连续内波波列分布，相较于理论解中单一的内孤立波作用，更能探讨实际内波环境场中水下航行体运动的累积效应。

本文采用MITgcm模型，模拟了连续密度分层背景下潮地作用生成内波的全过程。并通过自定义边界条件，将MITgcm模型中稳定传播的内波波列导入STAR-CCM+软件中，实现小范围水槽内连续密度分层环境下内波的精准模拟。在此基础上，针对约束状态的大尺度SUBOFF模型开展了连续内波波列作用下的数值模拟及水动力载荷分析，并分析归纳了不同潜深工况下的水动力载荷变化规律。

1 数学模型 1.1 MITgcm模型

本文采用的MITgcm模型是一种基于不可压缩的N-S方程建立的非静压海洋环流数值模型，控制方程如下：

连续方程：

$ \nabla \vec u = 0 。$

(1)

动量方程：

$ \frac{{\partial \vec u}}{{\partial t}} + \vec u \cdot \nabla \vec u + f{\vec e_3} \times \vec u = - \frac{1}{\rho }\nabla p + \nabla \cdot (\boldsymbol A\nabla \vec u)。$

(2)

热扩散方程：

$ \frac{\partial T}{\partial t}+\vec{u}\cdot\nabla T=\nabla\cdot\left(\boldsymbol{K}^{{T}}\nabla T\right)。$

(3)

盐量扩散方程：

$ \frac{{\partial S}}{{\partial t}} + \vec u \cdot \nabla S = \nabla \cdot \left( {{{\boldsymbol{K}}^S}\nabla S} \right)。$

(4)

海水状态方程：

$ \rho = \rho (T,S,p) 。$

(5)

式中：速度$ \vec u = ({{u,v,w)}} $；$ \rho $为海水密度；$T$为海水温度；$S$为海水盐度，海水压强$P$为未知数；$f$为科里奥利参数，$ {\vec e_3} $为垂直方向单位矢量，$ \boldsymbol A、{\boldsymbol K}^{ T}、{\boldsymbol K}^{S} $分别为粘性、温度、盐度扩散矩阵，其中下标$ h、v $分别表示水平和垂直方向：

$ \begin{split} & \boldsymbol{A}=\left(\begin{array}{*{20}{c}}A_h & 0 & 0 \\ 0 & A_h & 0 \\ 0 & 0 & A_v\end{array}\right),\quad\boldsymbol{K}^{{T}}=\left(\begin{array}{*{20}{c}}K_h^{{T}} & 0 & 0 \\ 0 & K_h^{{T}} & 0 \\ 0 & 0 & K_v^{{T}}\end{array}\right), \\ & \boldsymbol{K}^S=\left(\begin{array}{*{20}{c}}K_h^S & 0 & 0 \\ 0 & K_h^S & 0 \\ 0 & 0 & K_v^S\end{array}\right)。\end{split} $

(6)

1.2 连续密度分层数值模拟方法

本文采用流体体积法（VOF）来模拟密度跃层的连续分布，该方法通过计算各相不混溶流体在单元中的体积分数来确定该位置流体的性质。假设第$q$种流体在单元格中的体积分数为${\alpha _q}$，对第$q$种流体有：

$ \frac{{\partial {\alpha _{\text{q}}}}}{{\partial t}} + {u_i}\frac{{\partial {\alpha _{\text{q}}}}}{{\partial {x_i}}} = 0。$

(7)

且各相流体体积分数之和需满足：

$ \sum\limits_{q = 1}^n {{\alpha _q}} = 1。$

(8)

单元格内流体的密度为：

$ \rho = \sum\limits_{q = 1}^n {{\alpha _q}} {\rho _q}。$

(9)

2 计算模型

本文研究对象为放大20倍的SUBOFF全附体标准模型，主尺度数据见表1。

如图1所示，将模型置于长4 km，宽300 m，深300 m的数值水槽中。模型头部距速度入口边界1 km。SUBOFF近壁面区域采用标准壁面函数模型，根据y+值（范围为30～60），划分水下航行体边界层网格，并在附体附近进行网格加密，如图2所示。

对网格无关性进行验证，设计了网格总数为395万、448万、694万、1395万等4套网格，数值模拟计算水下航行体总阻力系数，结果如图3所示。发现448万网格和694万网格之间的阻力结果差在1%以内，因此认为进一步细化网格对结果影响较小，选择448万网格方案进行计算。

将不同雷诺数下的水下航行体表面摩擦力系数数值计算结果与ITTC公式计算的摩擦系数理论值对比，如表2所示，误差在4%以内。

3 数值模拟结果与分析 3.1 内波生成传播模拟

基于内波生成原理，本文首先在MITgcm中模拟二维内波的生成与传播过程。所建立的二维模型中x轴放置在水平面上，z轴向上为正，总水深300 m，水平计算域300 km。水平方向网格尺寸为50 m，垂直方向网格尺寸为2 m，地形设置为单一高斯海山地形，高180 m。初始时刻流场垂向密度剖面采用双曲正切形式，如图4(a)所示。

$ \bar \rho (z) = \frac{{{\rho _1} + {\rho _2}}}{2} - \frac{{{\rho _2} - {\rho _1}}}{2}\tanh \left( {\frac{{z - {z_{pyc}}}}{{{d_{pyc}}}}} \right)。$

(10)

式中：$ {\rho _1} $和$ {\rho _2} $为上下层流体密度，分别为1022、1028 kg/m³；$ {z_{pyc}} $为密度跃层中心位，位于水下50 m处；$ {d_{pyc}} $为密度跃层厚度，根据海洋实际环境，设为20 m。

模式左右入口为开边界，用M2潮流进行驱动，并设置海绵边界层进行消波。海表面采用缸盖假设，数值模拟结果见图4。在7 h时刻，内潮与山脊相互作用，在海山处出现大幅密度凹陷；在14 h时刻，海山右侧形成内波波列；在29 h时刻，内波波列中首波传播至180 km处，振幅可达83 m，相速度c约为1.263 m/s，次波与首波相距2 km，振幅为64 m，如图4(e)所示；与此同时下一潮周期生成的内波波列传播至100 km处。图4(e)中的内孤立波波列与邝芸艳等人给出的实际海洋中内孤立波包的地震剖面资料图类似，如图5所示，首波振幅最大，首波次波之间相距约2 km，振幅逐渐衰减。对比验证了该模型生成的内孤立波波列有一定普适性。对波形进行分析，传播至180 km处的首波波形进行分析。将首波等密度线与全非线性内波DJL理论解对比，如图6所示，发现MITgcm中密度跃层上边缘等密度线在波谷处吻合良好，下边缘等密度线振幅略大于理论解，整体上MITgcm中内波波宽小于DJL理论解。该部分差异主要来源于背景潮流的作用，使内波波形变窄。

3.2 内波数值水槽

在STAR-CCM+中建立长4 km、宽300 m、深300 m的三维水槽。水平方向网格尺寸为12.5 m，垂直方向网格尺寸均为1 m。上下边界采用固壁边界，左右边界均采用速度入口边界条件。

为在小范围内精确模拟MITgcm中所生成的大振幅内波首波，相当于将该4 km长的数值水槽置于MITgcm中的相应位置，如图4(e)中红框所示。与2层流体假设下静止稳定的初始条件不同，由于受到背景潮流的影响，以27 h时刻该位置处的密度、速度分布情况作为4 km长数值水槽的初始条件。记录MITgcm模型中水槽左右边界位置处的密度、速度时历数据，作为数值水槽左右边界的速度入口条件输入。在小范围数值水槽中生成内波波列，如图7所示，首部与次波之间相差约1850 m。针对该内波波列的首波，对比STAR-CCM+与MITgcm中同一时刻、同一位置处的密度等值线，如图8所示，波形误差在2%以内，可认为在STAR-CCM+中生成了与MITgcm中相同的内波。该方法实现了大小范围海洋环境的衔接，对内波波列的首波和次波进行了精确模拟，为后续水下航行体遭遇内波奠定了更实际的环境基础。

3.3 约束状态水下航行体穿越连续内波波列水动力载荷研究

将水下航行体固定在水下110 m处，即密度跃层下边缘以下位置。水下航行体遭遇内孤立波波列时所受水动力载荷时历变化曲线如图9所示。

力与力矩通过无因次形式给出，无因次形式为：

$ F_x^\prime = \frac{{{F_x}}}{{0.5\bar \rho {c^2}{L^2}}}，$

(11)

$ F_z^\prime = \frac{{{F_z}}}{{0.5\bar \rho {c^2}{L^2}}} ，$

(12)

$ {\boldsymbol M}_y^\prime = \frac{{{{\boldsymbol M}_y}}}{{0.5\bar \rho {c^2}{L^3}}}。$

(13)

式中：${F_x}$、${F_z}$和${{\boldsymbol M}_y}$分别为水下航行体所受水平力、垂向力、俯仰力矩；水平力向右为正；垂向力向上为正；以重心为转动中心，力矩顺时针为正；$\bar \rho $为平均密度，取1025 kg/m³；c为内波相速度。

由图9所示，随着内波的接近，固定水下航行体在t=834 s时遭遇首波，t=3634 s时遭遇次波，遭遇内孤立波过程中水动力载荷变化趋势类似，水平力都先达到负向最大值，再达到正向最大值；垂向力会在短时间内达到负向最大值；俯仰力矩先达到逆时针最大值，再达到顺时针最大值。而在首波、次波中间时，水动力载荷趋于0。但水动力载荷极值大小有一定差异，如表3所示。水下航行体遭遇内波右侧波面时受到的水平力、俯仰力矩极值都大于水下航行体穿过内波左侧波面时。且由于振幅减小，水下航行体在穿越次波时所受水平力和俯仰力矩极值略小于水下航行体穿越首波时，水平力最大值可达44.6 KN，俯仰力矩最大可达6757 kNm。水下航行体穿越主波时的垂向力极值约为穿越次波时的2倍，垂向力水动力载荷可达238 kN，对因为水下航行体穿越密度跃层时，垂向力变化主要由浮力决定，$ F_浮=\displaystyle\int_z\rho_zgV\mathrm{d}z-G $。由于振幅的衰减，相同潜深水下航行体穿越波面时相对于密度跃层的位置发生了改变。水下航行体在未遭遇内波时均处于下层流体中，所受垂向力平衡。在穿越内波右侧波面后位于波谷中时，如图10所示，在次波波谷中时水下航行体位于密度跃层中，而在首波波谷中时水下航行体位于密度跃层上边缘，$ {\rho }_{z首波} $<$ {\rho }_{z次波} $，因此$ {F}_{浮首波} $变化幅值要大于$ {F}_{浮次波} $。

3.4 潜深影响

对不同潜深工况下固定水下航行体遭遇内波波列中的大振幅首波进行数值模拟计算，模拟工况潜深相对位置如表4和图11所示，其中d为矩密度跃层中心的距离，密度跃层中心以上为正。

如图12(a)所示，水平力变化规律主要分为2类。从水平速度诱导的动压场角度分析（见图13），位于正向动压场的工况（d=+20，−20 m）由于水平速度先增大后减小，在波谷处达到最大，因次水平力先达到正向最大再反向增加达到负向最大，最后衰减趋于0。位于负向动压场的工况（d=−60，105 m）变化规律则相反。重点对比分析波谷之前时刻水平力变化幅值（见表4）。发现工况4水下航行体水平力幅值最大，且峰值出现时间最晚，接近波谷时在最短时间内水平力趋于0。越接近密度跃层下边缘，水平力幅值越小。

由图12(b)所示，位于密度跃层附近的工况（d=+20，−20 m）和穿越密度跃层的工况（d=−60 m），由于密度的变化，垂向力变化趋势类似。垂向力主要由浮力决定，随着内波的接近，水下航行体周围流体密度减小，垂向力迅速达到负向最大值。位于波谷以下的工况（d=−105 m），水下航行体一直位于下层流体中未穿越密度跃层，因此垂向力仅受内波诱导的垂直速度场影响。如表4所示，工况1最早达到垂向力最大值，因为随着内波的接近，工况1中水下航行体下方密度跃层下压，水下航行体快速全部进入上层液体中，达到垂向力最大值。工况2和工况3垂向力最大值出现时间相近，都为水下航行体位于波谷时刻，且垂向力幅值远大于工况1、工况4，比同工况下水平力约大1个量级。

由图12(c)所示，位于内波波高范围内的工况（d=+20，−20，−60 m），俯仰力矩变化规律类似，都先达到逆时针力矩最大值，再达到顺时针力矩最大值。而波谷以下工况（d=−105 m），变化趋势相反。但在该过程中位于波谷时刻，力矩都为0。对比各工况力矩极值点，如表4所示，越接近密度跃层，俯仰力矩极值越小；密度跃层以下波高范围内工况（d=−60 m）水下航行体俯仰力矩极值最大，比其他工况大1～2个量级；波谷以下工况未穿越密度跃层，俯仰力矩稍小，但也比工况1、工况2大1个量级多。

4 结　语

本文将大范围非静压海洋环流模型与小范围数值水槽相结合，完成了连续密度分层背景场下，内孤立波波列的生成传播与小范围高精度重现。采用数值模拟的方法计算了连续内波波列作用下的实尺度水下航行体约束状态所受到的水动力载荷，并进一步对比研究了不同潜深工况下水下航行体遭遇大振幅内波时的水动力载荷变化规律与极值。获得的主要结论如下：

1） 本文基于连续密度分层假设的非静压大范围环境场模拟内波波形与全非线性DJL内波理论解结果相当，波列形式与实际内波波列观测结果一致。并在小范围上对实际环境场内波波列进行了高精度重现，考虑了背景潮流对内波的影响，避免了2层流体假设下的垂向水平速度突变现象，为后续水下航行体载荷计算建立了更实际可靠的环境场基础。

2） 水下航行体在穿越连续内波波列首波、次波时，水平力、垂向力、俯仰力矩变化规律类似。但由于首波振幅大于次波振幅，穿越次波时的水平力和俯仰力矩极值略小于穿越首波时的。由于振幅不同，水下航行体穿越内孤立波时相对于密度跃层的位置不同，穿越次波时的垂向力幅值是穿越首波时的一半。

3） 不同潜深工况下水下航行体水动力载荷变化规律不同。水平力变化规律与内波诱导的水平速度相关，以水平正负动压场为界，变化趋势相反。垂向力变化主要分为两大类：水下航行体周围存在密度变化和一直处于同一密度流体中，前者垂向力比水平力大约一个量级。潜深位于波谷以下水下航行体俯仰力矩变化规律与其他工况相反，且越靠近密度跃层，俯仰力矩极值越小。潜深在波高范围内的水下航行体所受俯仰力矩极值最大，约为密度跃层处的20倍。在连续密度分层的内波场中，考虑了密度跃层厚度的影响，更准确地归纳了不同潜深水下航行体的水动力载荷变化规律，更有利于后续水下行体的运动研究。