0 引　言

水下拖曳体是进行水下资源探测和海洋资源开发的重要载体，在现代海洋开发、水声探测、军事和民用等诸多领域有着十分广泛的应用，为了满足各种拖曳航速和应用环境的需求，对水下拖曳体的性能提出越来越高的要求^{[1 − 2]}。拖曳系统通常由水面拖曳母船、牵引缆、水下拖体、尾绳等设备组成，拖曳形式也因用途的不同而各不相同^[3]。水下拖曳系统动力学的理论研究不断完善，国内外诸多学者针对各自的研究，建立了运动的数学模型，采用数值模拟等方法研究其在不同情况下的稳态、动态运动^[4]。Chen等^[5]基于CFD模拟和NSGA-II多目标优化算法的结合，对牵引式和自航式的自主水下航行器（Autonomous Underwater Vehicle, AUV）外形参数进行多目标优化，实现了航行时阻力低、运动稳定性好的设计，所提出的多目标优化方法有效地优化了AUV的外形，为AUV的设计提供了重要参考；Yang等^[6]采用CFD的方法对拖船和水下拖曳体的运动进行模拟，分析了水下拖曳体在波动运动模式下的水动力特性，结果表明，在周期控制作用下，系统的外力、内力以及拖曳系统各部件的位置和姿态呈现周期性变化，作用在拖曳体的压力幅值随着拖曳力的增大而增大，此外，当移动配重从后到前经过其平衡位置时，作用在牵引缆线上的流体动力分别达到最小值和最大值，反之亦然；孙烨等^[7]为了提升拖曳系统深度的调节性能，设计了一种通过内部调节机构改变运动姿态和深度的水下拖曳体，通过对水动力学数值计算，分析水下拖曳体重心位置与运动姿态的关系，验证内部调节机构的可行性，数值计算结果表明，不同俯仰角下水下拖曳体深度变化的时间有所不同，在俯仰角40°左右范围内，拖曳体的升力系数最大，并呈现较好的深度调节运动姿态；杨显原等^[8]采用重叠网格嵌套滑移网格技术对不同周期的水下拖曳系统的水动力特性进行数值分析，研究了水下拖曳系统控制动作与控制力、运动姿态、阻力以及拖曳缆绳拖曳力等特性及其相互之间的关系，结果表明，水下拖曳系统在控制动作下呈现出复杂的运动学及动力学特性，文中的水动力模型可对水下拖曳系统水动力特性进行有效的预测。

本文采用Fluent软件Realizable $k - \varepsilon$湍流模型对搭载声学设备的水下拖曳体进行了数值模型的建立，采用数值模拟方法，研究了在3种砰击载荷工况下水下拖曳体的水动力特性，重点研究水下拖曳体各部分在运动过程中的受力情况，以此来满足实际使用的需求。

1 数值模拟方法 1.1 控制方程

本文选取Realizable $k - \varepsilon$湍流模型^[9]，该模型中涡粘系数引入与旋转和曲率相关内容，已被有效地应用于包括带有分离的流动中。它是两方程模型，k和$\varepsilon$是2个未知量，其输送方程如下：

$\frac{{\partial (\rho k)}}{{\partial t}} + \frac{{\partial (\rho k{u_i})}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}\left[\left(\mu + \frac{{{\mu _i}}}{{{\partial _k}}}\right)u\left(\frac{{\partial k}}{{\partial {x_j}}}\right)\right] + {G_K} - \rho \varepsilon,$

(1)

$\begin{split} \frac{{\partial (\rho \varepsilon )}}{{\partial t}} + \frac{{\partial (\rho \varepsilon {u_i})}}{{\partial {x_i}}} = &\frac{\partial }{{\partial {x_j}}}\Bigg[\Bigg(\mu + \frac{{{\mu _i}}}{{{\partial _\varepsilon }}}\Bigg)u\Bigg(\frac{{\partial \varepsilon }}{{\partial {x_j}}}\Bigg)\Bigg] +\\& \rho {C_1}E\varepsilon - \rho {C_2}\frac{{{\varepsilon ^2}}}{{k + \sqrt {v\varepsilon } }}, \end{split}$

(2)

${C_1} = \max\Bigg (0.43,\frac{\eta }{{\eta + 5}}\Bigg),$

(3)

$\eta = {(2{E_{ij}} \cdot{E_{ij}})^{1/2}}\frac{k}{\varepsilon },$

(4)

${E_{ij}} = \frac{1}{2}\Bigg(\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}\Bigg),$

(5)

${\mu _{\text{t}}} = \rho {C_\mu }\Bigg(\frac{{{k^2}}}{\varepsilon }\Bigg),$

(6)

${C_\mu } = \frac{1}{{{A_0} + {A_S}{U^*}k/\varepsilon }},$

(7)

${A_S} = \sqrt 6 \cos \phi ,$

(8)

${U^{*}} = \sqrt {{{E}_ij}{E_{ij}} + {{\tilde \Omega }_{ij}}{\tilde \Omega }_ij} ,$

(9)

$\phi = \frac{1}{3}{\cos ^{ - 1}}(\sqrt 6 W),$

(10)

$W = \frac{{{E_{ij}}{E_{jk}}{E_{kj}}}}{{{{({E_{ij}}{E_{ij}})}^{1/2}}}} ,$

(11)

$\tilde \Omega _ij = {\Omega _{ij}} - 2\varepsilon {i_{jk}}{\omega _k},$

(12)

${\Omega _{ij}} = {\overline \Omega _{ij}} - {\varepsilon _{ijk}}{\omega _k}。$

(13)

式中：$\rho$为密度；t为时间，$\mu$为动力粘度；k为流体的传热系数；${\sigma _k}$=1.0；${\sigma _\varepsilon }$=1.0；${C_2}$=1.9；${A_0}$=4。

使用VOF方法描述波浪运动，该方法可通过体积分数F来描述波浪的运动，如下式：

$\frac{{\partial F}}{{\partial {t}}} + (\nu \cdot \nabla )F = 0。$

(14)

式中：当F =0时，表示网格为与空气相中；当0< F <1时，表示网格处于液相与空气相之间的自由表面；当F =1时，表示网格位于液相。

1.2 物理模型

拖体数值计算模型如图1所示，分别对头壳、上部浮力块导流罩、换能器导流罩、基阵安装面、水平翼板、水平襟翼、尾部组件及拖体整体进行受力监测区域划分。

1.3 边界条件及网格 1.3.1 边界条件

1）拖体边界条件，如图2所示，由于拖体模型对称，为减少计算量，对计算区域设置对称边界。在高海况下拖体上下往复运动砰击波浪涉及大幅运动模拟，在拖体上下往复运动方向设置足够运动空间，于拖体前方45 m（1.5倍波长）位置设置为流体域速度进口，拖体后方90 m（3倍波长）位置设置为流体压力出口。

2）砰击工况，根据拖体使用环境要求，取5级海况对拖体砰击载荷进行计算，其中波高2.1 m，波长30 m。设5级海况下拖曳平台（船）捕捉拖体后上下往复升沉运动满足下式：

$Z = A\sin \Bigg(\frac{{2{\text π} t}}{T}\Bigg) {Z}={A}\sin (2{\text π}{t}/{T}) 。$

(15)

式中：A为升沉运动幅值；T为运动周期。

拟分别计算拖体初始位置为拖曳点位于水线、拖曳点距离水线以上半波高、拖曳点距离水线以下半波高，在5级海况下以运动幅值1.5 m、周期5 s做升沉运动时受波浪的砰击载荷。

1.3.2 网格生成

1）网格无关性验证

基于流体动力仿真软件对拖体不同姿态不同襟翼攻角状态进行数值计算时，其结果准确性很大程度上取决于网格的数量和质量，因此有必要对拖体模型开展网格无关性验证，选取航速6 kn、仰首3°、水平襟翼攻角12°状态拖体开展网格无关性验证。

采用减小第一层网格高度和减小表面网格大小2种方式增加网格数量，当拖体流体升力随网格数量增加变化值小于1%且变化减缓时，确定最优的网格数量。仰首3°、水平襟翼攻角12°状态拖体边界层和表面网格划分如图3所示，网格无关性计算结果见表1。

分析表1不同网格数量升力变化值数据可知，当网格数量增大至236万后（表面网格基准值1.2 m，y+设置值100时），继续增大网格数量升力变化值仅为−0.7%，对于解决工程问题，数值计算结果已满足要求。网格数量为236万时，仰首3°、水平襟翼攻角12°状态拖体y+分布如图3所示，其y+值主要集中在60左右，满足壁面函数使用要求。

2）网格生成

综合考虑计算精度及计算周期，将236万网格对应网格参数应用于后续拖体流体动力计算中。采用切割体重叠网格，具体网格划分如图4所示。

2 计算结果及分析

共对拖体进行3种砰击载荷工况计算，其中拖体初始位置分别为拖曳点距离水线以上半波高、拖曳点位于水线、拖曳点距离水线以下半波高，如图5所示。对拖体升沉运动进行实时监测，如图6所示。

由图6可知，拖体按照仿真输入要求进行了上下往复升沉运动。在高海况下拖体砰击水面过程，如图7所示。

对拖体砰击载荷各工况进行32 s时间历程仿真计算，得到各监测区域受力历时曲线，如图8所示。其中，各监测区域受力仅为流体力，不包含惯性力及重力。

在流体力中添加重力后进行统计，各监测区域受力极值结果，如表2所示。

拖体实际工作中，在各受力点处测出试验结果，试验结果如表3所示。

由于水下拖曳体为对称分布，故数值模拟计算中采用一半对称计算，因此表2监测到的受力为实际受力数据的一半，与表3实际试验结果相比出现了略微大的情况，这是由于在数值模拟中未考虑惯性力的影响。

3 结　语

本文对高海况拖体收放抨击载荷进行数值计算和试验测试，得到以下结论：

1）对拖体进行3种抨击载荷工况计算，其中初始位置分别为拖曳点距离水线以上半波高，拖曳点位于水线，拖曳点距离水线以下半波高。

2）对拖体受力监测区域进行仿真计算和试验测试，分别得到各监测区域受力情况，相同条件下的实验与数值模拟结果基本吻合，验证了数值模拟结果的有效性。

3）本文研究内容为水下拖曳体的设计提供了理论支撑和依据。