0 引　言

随着海洋科技的不断发展，无人水下航行器（Unmanned Underwater Vehicle，UUV）作为一种高效、自主的水下探测，受到了广泛关注与研究。UUV不仅能够在复杂多变的水下环境中执行各种任务，还因其使用费用低、作战半径大的特点广泛应用于探测^[1]、通信^[2]、环境监测和目标识别^[3]等商业或军事领域。在双UUV协同作业场景中^[4]，如何优化其流场特性，提高作业效率，成为了亟待解决的问题。机械臂作为UUV的重要执行机构，即在UUV的基础上加装机械手构成水下机器人-机械手系统（Underwater Vehicle-Manipulator System，UVMS），在双UUV协同作业中扮演着至关重要的角色。当2个UUV搭载机械臂进行水下作业时，它们之间的流场特性将直接影响到机械臂的操作精度和作业效率。因此，对带机械臂的双UUV流场特性进行数值模拟，分析水动力特性，对潜器设计研发有重要的理论和实际意义。

近年来，国内外学者针对无人潜器的流场仿真分析开展了广泛研究^{[5 − 8]}。吴利红等^[9]采用“移动网格法”对不同航速潜器的直航运动进行水动力预报。王佳茂等^[10]对潜航器进行了水动力分析计算，揭示了UUV在不同浪向下的运动响应特征。Korol等^[11]利用CFD方法，分析了潜器的水动力性能随距海底距离和攻角的变化特性，并优化了潜器受海底干扰条件下的水动力参数计算方法。此外，Mitra 等^[12]采用试验与数值模拟结合的方法，探究了倾斜河床对潜器工作性能的影响。然而，上述研究主要聚焦于单潜器的水动力学特性分析，在潜器集群化作业日益成为发展趋势的背景下，多潜器协同作业工况下的水动力特性分析具有广泛的现实意义^{[13 − 15]}。Rehman等^[16]考察了一种半经验方法获得由2艘航行器组成的运输系统水动力系数的适用性，对阻力进行验证，半经验计算结果与CFD计算结果一致。多潜器协同工作往往存在于复杂海洋环境，相关研究具有特殊性，机械臂作为无人潜器的重要执行机构，其对流场的干扰研究相对较少。Du等^[17]建立了水下潜器近海底运动的数学模型，得出了阻力、升力、俯仰力矩系数的绝对值随攻角的增大而增大的结论。

由于机械臂的加装以及潜器间流场的相互扰动，使得作业过程中，流场环境更为复杂多变。本文采用计算流体动力学（CFD）方法开展了双UUV水下协同作业的流域分析，旨在分析不同间距下UUV的水动力系数变化规律。在有无机械手系统加装的不同情况下，开展潜器流场模拟仿真分析，以获取不同工况下潜器的阻力、升力、压力及速度分布情况。

1 数学模型和数值方法 1.1 控制方程

流体的动力学特性遵守三大守恒定律，即质量守恒定律、动量守恒定律和能量守恒定律^[18]。对于不可压缩流体，流体的密度恒定不变，密度的物质导数为0，连续性方程可表示为：

$ \nabla \cdot U = 0。$

(1)

式中：$ U $为流体速度。

动量守恒原理表明，若没有外力施加在物体上，该物体保持质量与速度的乘积不变，即N-S方程为：

$ \rho\frac{\mathrm{d}U}{\mathrm{d}t}=\rho F_b-\nabla p+\nabla\cdot(2\mu S)-\frac{2}{3}\nabla(\mu\nabla\cdot U)。$

(2)

式中：$ \rho $为密度；t为时间；$ {F_b} $为体力；p为压强；$ \mu $为动力粘度系数，S为应变率张量。该式等号左端为单位体积流体的惯性力，等式右端的第一项$ \rho {F_b} $为作用在单位体积流体上的质量力，第二项为作用在单位体积流体上的压强梯度力，第三项为粘性变形应力，最后一项为粘性体膨胀应力。

1.2 湍流模型

湍流是种常见的复杂流动，目前数值计算可以对湍流进行模拟。湍流模型分为涡粘和雷诺应力模型。涡粘模型中的$ k - \varepsilon $模型和$ k - \omega $模型被广泛应用在实际工程。参考相关文献及相应仿真模拟^{[19 − 20]}，本文采用SST $ k - \omega $模型。SST $ k - \omega $模型在自由剪切层和边界层边缘处使用$ k - \varepsilon $模型，在近壁面自动转换成Wilcox $ k - \omega $模型。因此，SST $ k - \omega $模型既擅长远场计算，也能在近壁面进行准确的计算，使得计算结果更为精确可信。SST $ k - \omega $模型输运方程式为：

$ \frac{\partial }{{\partial t}}(\rho k) + \frac{\partial }{{\partial {x_i}}}(\rho k{u_i}) = \frac{\partial }{{\partial {x_j}}}\left({\Gamma _k}\frac{{\partial k}}{{\partial {x_j}}}\right) + {G_k} - {Y_k} + {S_k}, $

(3)

$ \frac{\partial }{{\partial t}}(\rho \omega ) + \frac{\partial }{{\partial {x_i}}}(\rho \omega {u_i}) = \frac{\partial }{{\partial {x_j}}}\left({\Gamma _\omega }\frac{{\partial k}}{{\partial {x_j}}}\right) + {G_\omega } - {D_\omega } + {S_\omega }。$

(4)

式中：i和j取值范围为（1,2,3）；$ {\Gamma _k} $和$ {\Gamma _\omega } $为k和$ \omega $的有效扩散项；$ {G_k} $和$ {G_\omega } $为平均速度导致的湍动能k和$ \omega $的产生项；$ {Y_k} $为湍流扩散项；$ {D_\omega } $为交叉扩散项；$ {S_k} $和$ {S_\omega } $为自定义源项。$ {\Gamma _k} $、$ {G_k} $和$ {\Gamma _\omega } $的表达式如下：

$ {\Gamma _k} = \mu + \frac{{{\mu _t}}}{{{\sigma _k}}}，$

(5)

$ {\Gamma _\omega } = \mu + \frac{{{\mu _t}}}{{{\sigma _\omega }}}，$

(6)

$ {G_k} = \min ({G_k},10\rho {\beta ^ * }k\omega )。$

(7)

式中：$ {\sigma _k} $和$ {\sigma _\omega } $为普朗特数；$ {\beta ^ * } $为修正系数。

本文采用的湍流模型为SST $ k - \omega $模型。

2 物理模型

海底洋流环境较复杂，实际海底作业中有诸多不确定因素，因此，单一水下潜器并不能顺利地完成水下任务，这就要再引入一个潜器协同作业。为便于区分，本文对潜器编号，分别称为潜器1和潜器2。

潜器1长为2.68 m，宽为1.60 m，高为0.69 m，由首部、尾部和平行中体3个部分构成。平行中体长为1.48 m，首部和尾部的长均为0.6 m，平行中体通过顺滑的曲线与首尾部过渡。为减少网格数便于计算，需对机械臂进行合理简化。如图1所示，在艇体两侧分别有一个长为1 m，半径为0.05 m类圆柱体的结构。

潜器2艇身长为2.62 m，宽为1.10 m，高为0.73 m，首部长为0.6 m，通过顺滑曲线过渡到中段，尾部宽为0.5 m，推进器在净深为0.2 m的开口处。为顺利进行仿真，不考虑如螺旋桨一类零散结构的影响来进行模型简化。潜器两侧机械臂简化为长0.75 m、半径0.05 m的类圆柱体，潜器2模型如图2所示。

3 数值模拟和结果讨论 3.1 计算域设置和网格划分

计算域是流体流经的区域，本文将计算域长度设为潜器长度的6倍，宽和高均为潜器长度的3倍。如图3，潜器长度记为L，计算域尺寸为$ 6L \times 3L \times 3L $。对计算域离散化的过程就是划分网格。在STAR-CCM+中有3种不同的体网格，本文选择的切割体网格适用于表面曲率变化大、外形复杂的模型。同时，为精确地观测潜器周围流场变化，本文采用了棱柱层网格。网格划分如图4所示。

3.2 网格收敛性验证

为计算结果的准确，需对网格数量进行收敛性分析。本文模拟潜器1在来流速度为5.144 m/s的情况下，潜器受到的总阻力。生成了154万、218万、304万3种不同网格数量模型计算潜器的阻力值，这里选取的网格收敛率为$ \sqrt 2 $，调整网格基础尺寸，生成的网格数量以$ \sqrt 2 $倍增加，如表1所示。其中，总阻力为压差与剪切阻力的线性叠加和，即:

$ R = {R_p} + {R_f}。$

(8)

式中：R为总阻力；$ {R_p} $为压差阻力；$ {R_f} $为剪切阻力。

网格的精细度随网格数量的增加而增加，由表可知，218万和304万网格的计算结果基本一致，相差1.24 N，这表明网格的增长对计算结果影响极小，计算结果已经收敛。但154万和304万网格的计算结果相差36.85 N，差值较大。考虑计算成本及结果精确度的问题，本文采用218万网格进行仿真计算。

3.3 单一潜器特征数值分析 3.3.1 单一潜器流场分析

针对潜器1，图5(a)是在来流速度为5.144 m/s流场中的水平截面速度分布云图，潜器1的头部、机械臂与艇身连接的凸起处、机械臂与艇身的夹角处、潜器的尾部以及机械臂的尾部均出现了低速区。其中潜器头部和机械臂与艇身连接凸起处的低速区是由于流体与潜器模型壁面发生直接碰撞形成的。流体的粘性使臂和艇身连接拐角处的摩擦力增大，这种摩擦效应使流体的流速变小，形成了一片低速区。潜器和机械臂的尾部出现低速区是由于该区域曲率变化较大，流经尾部的流体具有不同速度方向，不同速度方向相互抵消，导致尾部区域的速度减小形成低速区。而机械臂与艇身的夹角处是由于流体具有惯性、粘性使流速分布不均匀，在夹角内形成了一片低速区。

如图5 (b)所示，中纵剖面是潜器模型的对称平面。由图可知，潜器1头部和尾部流速减慢，尾流场在出现一片低速区后，其他区域流体汇入到尾流区，使尾流速度有所增加，并持续增大到与初始来流速度相近。由于潜器模型的头部和尾部区域表面曲率变化较大，流速的不均匀性较明显，形成了驼峰状的高速区，该区域流速的增加值是最大的，称为“M”型速度剖面，即在潜器的中部出现了速度凹陷的现象。

潜器2模型上下不对称，导致纵剖面的流场环境更加复杂。如图6(a)所示，速度云图的上下两侧不再对称分布。其头部和尾部与潜器1一样产生了低速区。与潜器1不同的是，潜器2模型上下两侧及尾部均有较大的折角出现，尤其在潜器下侧，曲率变化较大，因此在这些区域出现了不同程度的低速区。云图中同样出现了“M”型速度剖面，且由于潜器头部曲率变化相较于潜器尾部曲率变化更显著，导致潜器头部的“驼峰”明显大于潜器尾部的“驼峰”。

潜器2水平截面的速度云图如图6(b)所示，潜器2是左右对称模型，其流场中的速度呈对称分布。高速区是很明显的“M”型分布在潜器的左右两侧。但潜器2尾部的低速区的面积明显要比潜器1大，这是潜器2尾部两侧凸起的形状导致的，尾部U形凸起阻碍了更多流体去向，使得低速区面积增大。

无臂的潜器1计算结果如图7所示。没有机械臂的干扰，从水平截面的速度云图中可知，尾流场的低速区面积明显减小了，且尾流区域速度有所增加。无臂潜器2的结果如图8所示，高速区的“M”型分布区域面积有所减小，且流速有所增加。这是由于拆除机械臂后，潜器外观形状变化所导致的。

3.3.2 单一潜器的阻力分析

通过仿真计算，在来流速度为5.144 m/s的条件下，潜器1、潜器2、无机械臂的潜器1及潜器2的阻力值如表2所示。可知，压差阻力在合阻力中所占的比例约80%，起主导地位。潜器1在加装机械臂后，受到的压差阻力大幅减小了322.57 N，而剪切阻力增加了69.03 N。因此在加装机械臂后，合阻力约减少了18.46%，但升力增加了1.47%，合阻力的变化幅度比升力大，即机械臂加装对阻力有较大影响，对升力影响较微弱。阻力值与潜器表面的形状密切相关，加装机械臂后，潜器外部形状更趋近于流线型，减小了阻力值。因此合理优化机械臂的外形，可提升潜器的快速性。潜器2在加装机械臂后，压差阻力减小了139.93 N，剪切阻力增加了48.72 N，因此加装机械臂后合阻力减少了4.66%。且升力减少了0.18%。阻力的变化比升力小。因此，相较于潜器1，潜器2机械臂的加装对潜器性能的影响不大。

对比2个潜器的仿真结果可知，潜器1的合阻力小于潜器2，但升力却大于潜器2，这是2个模型不同外型差异导致的。优化潜器主体外型能有效减小潜器的阻力值，且潜器1比潜器2的效果显著提高了13.80%。

3.4 双潜器计算分析 3.4.1 计算模型

为顺利进行海底作业，本文对双潜器的流场干扰和水动力特性进行研究。以潜器1和潜器2为研究对象，改变潜器间距离，探究不同间距潜器周围流场变化。计算模型如图9所示，以x轴为基准，潜器1和潜器2分别沿y轴负向和正向移动相同距离。潜器1和潜器2间实际距离记作d，平移所需的距离记作s：

$ s = - \frac{{d + D}}{2}。$

(9)

式中：$ D $为2个潜器最大直径的平均值，值为1.35 m。

此外，本文采用的物理模型是不定常的，为提升工作效率，对比了1.543 m/s、3.086 m/s、5.144 m/s以及6.173 m/s不同流速的收敛时间快慢，本文选择了收敛较快的5.144 m/s，该流速较3.086 m/s和6.173 m/s的收敛时间缩短了约10%，较1.543 m/s的收敛时间缩短了约16%。综上，设置了如表3所示的5组工况。为便于分析，引入了无量纲间距$ d/D $。

3.4.2 水平剖面速度

由图10可知，在无量纲间距$ d/D $小于2时，2个潜器首部前椭圆状的低速区连成了一片，且随着潜器间距离的增大而分离。当$ d/D $=3时，2个潜器首部前的低速区已完全分离，且潜器1低速区面积略大于潜器2，这是由于潜器1的迎流面积大于潜器2导致的。当$ d/D $小于1时，2个潜器间高速区融合在一起，两者的相互干扰情况较显著。当$ d/D $=1.5时，两者间高速区虽已分离，但图中仍能明显看出两者会相互影响。直到$ d/D $=2时，高速区才分离开。但高速区外围的部分，在$ d/D $=3时仍连成了一片，这表明2个潜器并不能够摆脱彼此造成的流场扰动。针对潜器尾部的低速区，潜器2低速区域面积大于潜器1，这是由于潜器2尾部凸起位置阻挡了更多流体的去路，导致流速变缓的情况更加显著。当$ d/D $<1.5时，2个潜器尾流场相互干扰情况比较强烈。尤其在$ d/D $=0.5时，2个潜器尾流场的低速区有相互融合现象。当$ d/D $=1时，2个潜器尾流场低速区呈现出分离的现象，即随着间距增大，潜器间的相互干扰逐渐减弱。综上，双潜器在水下协同作业时，要考虑到由于他们之间相互干扰进而导致流场环境更加复杂的情况。

3.4.3 中纵剖面速度

由图11可知，在$ d/D $=0.5时，潜器两侧外围的速度等值线在艇身中间处收缩成驼峰状。当$ d/D $>0.5后，首部和尾部上下两侧的高速区分离开，不再连成一片。这表明，在$ d/D $大于0.5后，2个潜器在中纵剖面方向上，流场相互干扰现象是微弱的。但不同间距下潜器2中纵剖面的速度云图差别并不大，即在中纵剖面方向上，潜器2受到流场干扰的影响并不强烈。

3.4.4 潜器表面压力

流场中压力与速度的变化是相关的。为便于观察分析，本文选取在$ d/D $=0.5协同作业的工况和单潜器作业的工况进行对比。当壁面周围流体流速减小时，压强会增大，反之亦然。2个模型间流体的速度随间距的减小而增大，导致模型表面压力分布具有不对称性。2个模型间相邻侧，即潜器1右侧和潜器2左侧表面压力较小。如图12所示，潜器1单独作业时，潜器表面左右两侧压力分布是完全对称的，而在$ d/D $=0.5的协同作业情况下，表面压力分布不再对称，且右侧的压力明显比单一作业时小，而潜器左侧表面的压力分布差别不是很大。针对潜器2，如图13所示，这种现象更加显著。相较于单一作业，协同作业情况下，潜器2左侧表面靠近艇臂接触区域的高压区几乎消失了，即潜器表面压力减小。而潜器右侧表面的压力分布差别不是很大。此外，这种两侧差异现象会随着间距的增加而逐渐消除。

3.4.5 阻力分析

数值仿真采用了湍流模型以及非定常流，所以阻力的数值随着时间的变化而变化，最终表现在一定范围内小浮动波动的状态，即为阻力值的收敛状态。具体的阻力值如表4所示。

显而易见，剪切阻力占据总阻力的份额较少，压差阻力的大小是影响总阻力数值大小的关键因素，这反映出潜器在水下的压力差极大地阻碍潜器向前运动。随着潜器1和潜器2之间距离增加，剪切、压差以及总阻力的值随之减小，如图14所示。随着两者间的距离增大，各项阻力值趋近于单一潜器所受的阻力值。因此，在涉及多个水下航行器的快速性设计时，要充分考虑到潜器间的间距这一重要因素。

4 结　语

1）在来流速度为5.144 m/s的流场环境中，潜器1和潜器2的速度剖面云图均出现了“M”型的高速区，潜器1受到的合阻力值为1120.29 N，潜器2受到的合阻力值为1864.44 N。由于模型外形的差异，潜器1受到的合阻力值比潜器2受到的合阻力小744.15 N。即同样的工作环境下，潜器1会更加灵活易操纵。

2）潜器1在加装机械臂后合阻力值减小了253.55 N，加装机械臂的潜器2合阻力减小了91.21 N。机械臂的加装使得潜器的外形更加趋近于流线型，使得合阻力值有所减小，且潜器1机械臂的加装在提升潜器性能方面比潜器2显著提高了13.80%，即潜器1机械臂的加装效果更好。

3）双潜器在不同间距工况下的计算结果表明，随着2个潜器之间距离的增加，各项阻力数值大小随之减小，并趋近于单一潜器所受到的阻力值，这关乎潜器快速性设计等相关问题。