基于类物理数值模拟进行海洋系统的操纵性能预报，最早由美国海军在20世纪90年代提出。这种数值模拟由于直接模拟高速旋转的推进器运动，导致模型非常复杂，网格尺度不一，高速与低速相对运动耦合，使数值模拟的计算时间步长较小，大大增加了数值计算量。同时各个部件交界面复杂，降低了数值精度。这限制了类物理数值模拟的发展。直到2010年以后，CFD(Computational Fluid Dynamics)软件和超算中心的快速发展，使类物理数值模拟有了实质的进展，主要采用重叠网格法在超算中心上运行。如Carrica等[2-4]采用离散螺旋桨对水面船舶的自航运动、自由深沉运动、横漂运动进行了操纵性模拟；Chase和Carrica[5]对全附体潜艇自航运动采用动态重叠网格法进行了数值模拟；Mofidi和Carrica[6]对集装箱船(KRISO Container Ship，KCS)带全附体桨舵的ZigZag操纵运动，采用移动舵和直接螺旋桨离散的方法进行数值模拟；沈志荣[7]对水面船舶操纵运动进行数值模拟。
重叠网格法和移动网格法都是运动边界模拟的方法。移动网格法发展较早，广泛应用于边界移动的数值模拟[8-12]，如机载外挂物分离运动，二维机翼和三维机翼的振荡运动，火箭发射的二维模拟，鱼群中等幅度的摆动，AUV从鱼雷管发射和AUV水下对接等运动模拟。但是，相对于重叠网格的快速发展，已经应用于类物理数值模拟，移动网格则发展较慢，甚至处于停滞不前的状态。主要是因为移动网格技术随着边界运动，网格跟着变形，导致网格畸变和网格大量更新，网格数增加，由此导致计算精度变差和计算时间增长，这限制了移动网格法在三维大位移高频响应问题中的应用[10]，使移动网格技术主要应用于小幅运动、单自由度运动和二维运动等的模拟，还未曾出现关于高速螺旋桨推进的载体操纵运动的模拟。解决移动网格技术的难点是构建合理的网格，采用有效的动网格更新技术。笔者将探索一种多块混合网格结合动态层方法[13-14]，采用与重叠网格法相似的子区域运动方法，在台式机上实现AUV自航对接的类物理数值模拟，从网格技术方面提高移动网格法的计算精度和计算效率。本文将为航空和水下载体复杂操纵运动的类物理数值模拟提供方法，提高复杂载体系统的安全操纵范围。
1 模型和方法 1.1 数学模型 AUV水下对接系统模型如图 1所示。AUV为改进型类REMUS AUV。AUV直径d=0.191m, 长度L=1.74m。螺旋桨为MAU4-40桨，直径D=0.156m。dock包含开孔的圆锥结构nozzle和框架式的圆管结构tube两部分[12]。dock总长为0.84L, nozzle开口直径为5.6d, tube直径为1.8d。起始时刻AUV与dock的相对距离约为5L。其中：A、B和C点分别是对接入口点A、井口B和终点C。大地坐标系OXYZ的原点建立在初始时刻AUV的首部端点，OX指向AUV对接方向，OY垂直向上，OZ指向右舷。
 |
| 图 1 三维对接系统模型 Fig. 1 Three-dimensional model of docking system |
| 图选项 |
1.2 多块动态混合网格法 网格系统为多块混合网格, 由多块结构网格和非结构网格构成[15]，如图 2所示。整个流场共划分为5个区域，其中区域①、③和④为结构网格区域，分别对应螺旋桨尾流拉伸区域，AUV及其边界层区域，AUV和dock之间的压缩区域。
 |
| 图 2 AUV对接的对称面网格 Fig. 2 Symmetrical mesh for AUV docking |
| 图选项 |
区域②为螺旋桨及其尾流场的非结构网格区域。区域⑤为包括dock的外围非结构网格区域。由于AUV尾部流场包括高速旋转的螺旋桨，以及AUV和螺旋桨耦合的复杂尾迹流，为捕捉其运动，对这部分区域进行网格加密，如图 3所示。整个流场网格总数为1.073×108：结构网格总数为5.03×107，非结构网格总数为5.7×107，其中螺旋桨区域②为4×107。图 4为AUV带舵翼和螺旋桨的网格图。AUV对接dock过程中，区域①、②、③和④随AUV做直航运动，区域②同时随螺旋桨作旋转运动。交界面Interface6和交界面Interface7静止，紧邻动区域一侧分别进行拉伸和压缩。这种多块混合网格法，采用结构网格划分近场区域，应用子区域移动代替边界移动，有利于网格移动过程中近场网格质量保持，提高数值求解精度[16]；同时采用动态层方法进行网格压缩和拉伸，有利于网格快速更新和网格数量的保持，提高数值求解效率。
 |
| 图 3 AUV和螺旋桨尾部网格放大图 Fig. 3 Zoomed-in mesh astern of AUV and propeller |
| 图选项 |
 |
| 图 4 AUV网格 Fig. 4 AUV mesh |
| 图选项 |
1.3 用户自定义函数 AUV自航对接过程采用用户自定义函数(User Defined Function，UDF)来实现AUV、螺旋桨和流域之间的力和速度传递。这包含3个UDF模块：螺旋桨模块、AUV模块和运动流域模块，具体流程如图 5所示。
 |
| 图 5 UDF流程图 Fig. 5 Flowchart of UDF |
| 图选项 |
1) 在螺旋桨模块中，螺旋桨以一定的初始速度V和转速n运动，产生推力T，将推力T存储到文本文件Fp中。
2) 在AUV模块中，读取Fp文件中的推力T，并传递给AUV，根据推力T和流场求解的阻力R，求解AUV 6DOF方程获得AUV新的速度V，将V、R存入文本文件Fa中。同时AUV以速度V航行。
3) 在运动流域模块，区域①、②、③、④以速度V移动。
4) 交界面Interface6和Interface7紧邻动区域的网格采用动态层法分别进行拉伸和压缩。
5) 区域⑤的网格进行变形，各交界面传递通量。
重新更新流场、边界位置和速度，螺旋桨从Fa文件中读入新的速度V，从步骤1)开始循环，直到AUV航行到dock的终端C，此时航行距离LS=LEND(LEND为AUV起点到dock终点距离), 对接过程终止，保存所有中间变量T、R、V和LS等。
1.4 数值方法和验证 本文采用有限体积法求解三维非定常不可压缩雷诺平均Navier-Stokes方程。空间离散采用基于压力方法，2阶精度。时间离散采用隐式离散法，外循环时间步长为定时间步长，每步时间为螺旋桨旋转1°对应的时间，内循环迭代20步。湍流模型为两方程k-ω SST湍流模型。物面边界条件为无滑移条件。数值设计的难点在于数值模型的建立、多块混合网格的网格划分和组装、UDF程序的编写。在进行非定常AUV自航数值模拟之前，将定常收敛结果作为初始值进行迭代计算。本文的载体和螺旋桨的网格无关性验证见文献[12]。
为了验证本文数值方法的准确性，将螺旋桨敞水试验的水动力计算结果和试验结果对比，如图 6所示, KT、KQ、J和η0分别为螺旋桨的推力系数、扭矩系数、进速系数和敞水效率。可见数值模拟和试验结果吻合较好。将AUV自航达到匀速下的推进系数与多参考系坐标(Multiple Frames of References, MFR)的进行对比，两者误差在5%以内，见文献[16]。此外，对AUV由静止开始，在螺旋桨定转速300r/min推进作用下进行自航试验，将航速变化的数值模拟结果和试验结果对比，如图 7所示，两者吻合较好。这些验证了本文方法进行数值模拟的精度。
 |
| 图 6 敞水试验的数值和试验结果对比 Fig. 6 Comparison between numerical and experimental results in open water tests |
| 图选项 |
 |
| 图 7 AUV自航速度曲线 Fig. 7 AUV velocity curves during self-propulsion |
| 图选项 |
2 数值结果分析 数值模拟AUV自航对接过程，AUV在恒定转速300r/min作用下，从静止开始加速，总迭代约28000步，历时约16s，最终速度达0.75m/s。对接速度满足对接碰撞需求(文献[17]要求对接速度小于1m/s)。根据AUV相对dock的位置变化，将AUV速度分为3段。远离dock阶段1：0~12s；近距离dock阶段2：12~14.4s；进入dock阶段3：14.4~16s。阶段1不受dock的影响。因此，可以认为阶段1为AUV的自航加速段，将此段与AUV单独在水池中的自航试验结果对比，如图 7所示，AUV起始加速度较大，最终速度渐趋于稳定。
图 8为AUV自航对接过程推力和阻力随速度变化。初始时刻，AUV速度为0，AUV阻力为0，螺旋桨推力最大。这与螺旋桨敞水性能一致。由于AUV推力大于阻力，AUV做加速运动。随着AUV速度增加，阻力增加，推力下降，速度仍然增加，但是加速度逐渐变小。最后，当推力和阻力接近时，速度达到近似匀速状态(见图 7速度曲线)。全对接过程中，当AUV进入dock后，阻力和推力有微小振荡。
 |
| 图 8 推力和阻力曲线 Fig. 8 Curves of thrust and resistance |
| 图选项 |
AUV进入dock阶段3的受力和速度变化如图 9和图 10所示。AUV进入A点后，其阻力稍微增加，AUV速度仍然增加，曲率稍微减小；经过连接点B，AUV阻力开始下降，直到阻力降到最低点D点，对应此区域的AUV的速度曲率增大；当AUV经过D点，直到终端C，阻力趋于稳定值，此时速度曲率又变小。dock对AUV影响的关键点是B点，在B点之前，dock阻碍AUV前进；在B点之后，dock吸附AUV对接。终止点C相对初始点A，AUV阻力增加2.4%，速度增加1.75%。整体而言，dock对AUV阻力和速度影响并不大。整个过程推力呈现微幅振荡，大小随速度增加略有下降。推力振荡的原因是由于尾部舵翼的非轴对称分布，其尾流场影响了螺旋桨的进流场，使推力出现周期性振荡[2]。
 |
| 图 9 AUV推力和阻力(在A和C点之间) Fig. 9 AUV thrust and resistance (between point A and C) |
| 图选项 |
 |
| 图 10 AUV速度(在A和C点之间) Fig. 10 AUV velocity (between point A and C) |
| 图选项 |
图 11为AUV舵翼和螺旋桨尾迹的速度矢量图，螺旋桨进流场可见舵翼后的回流，其速度与AUV航行方向相同。螺旋桨叶片后方可见梢涡和毂涡的旋转方向。图 12为AUV对接非定常过程6个典型时刻(t=0.6, 3.2, 5.8, 10, 14.8, 16s)轴向速度云图。左侧为螺旋桨尾迹局部放大图，右侧为全图。由于螺旋桨尾部速度场变化明显，速度范围跨度较大，因此采用2个图例显示尾迹场。在t=14.8s和t=16s的2个云图中，采用二维图显示dock的速度流场。AUV自航对接过程中，螺旋桨叶梢曳出梢涡，梢涡强度随AUV速度增加而变弱。螺旋桨尾迹中靠近桨毂后方存在毂涡，与AUV速度方向一致，涡强随速度增加而增大。随着AUV航行距离增加，AUV尾迹线不断延伸，强度逐渐变小。当AUV进入dock(t=14.8s)，AUV的首部速度场范围扩大，流体获得加速，AUV阻力增加。
 |
| 图 11 AUV尾部速度矢量图 Fig. 11 Velocity vectors astern of AUV |
| 图选项 |
 |
| 图 12 AUV自航对接轴向速度云图 Fig. 12 Axial velocity contours in AUV docking by self-propulsion |
| 图选项 |
图 13为AUV对接非定常过程的6个典型时刻(t=0.6, 3.2, 5.8, 10, 14.8, 16s)的三维压力云图。由图 13可见，螺旋桨刚开始定转速旋转时，AUV的速度还未显著增加，AUV表面的压力没有变化(t=0.6s)。随着AUV速度增加，螺旋桨叶面压力变小，如左侧图所示。而AUV表面出现了3个明显的压力区：AUV首部肩颈、平行中体和AUV尾部肩颈。AUV的平行中体压力稍大。首部肩颈和尾部肩颈处为低压区，随着AUV速度增加而范围扩大，压力下降(t=3.2, 5.8，10s)。AUV首部肩颈处的低压区表明：AUV首部驻点压力在增加。在大地坐标系下，AUV首部高压作用于流体，流体会向AUV尾部流动，高压越大，靠近高压区的AUV首部肩颈处的低压值越小，此处对应的速度越大，如图 12所示，摩擦阻力增加，则AUV阻力增加。当AUV首部进入dock(t=14.8s)，nozzle上呈现高压分布。当AUV平行中体进入dock (t=16s)，nozzle上呈现低压分布。这表明dock在AUV进入B点前后呈现相反的作用。B点之前，阻碍AUV航行；B点之后，吸附AUV对接。
 |
| 图 13 AUV自航对接的三维压力云图 Fig. 13 Three-dimensional pressure contours of AUV docking by self-propulsion |
| 图选项 |
3 结论 本文应用多块混合网格结合动态层方法实现AUV自航对接dock的类物理数值模拟。该数值模拟在台式机(i5-6400 CPU @2.70GHz, 2.70GHz, 内存16.0GB)上以4个处理器并行计算，历时约14天(文献[2-5]采用超算中心计算，需要1个月)。数值仿真得出如下结论：
1) AUV从静止开始，在螺旋桨恒定转速300r/min的推进下，自航对接历时约16s，最终航速达到0.75m/s。
2) 螺旋桨旋转推进AUV自航中，桨叶叶梢向下游曳出梢涡，梢涡强度随AUV速度增加而变弱。靠近桨毂存在毂涡，与AUV的速度方向一致，强度随AUV速度增加而增大。AUV尾迹线随AUV航行距离增加不断拉长，强度逐渐变小。
3) AUV进入dock后，阻力和推力出现振荡。dock对AUV的影响以B点为分界点，在B点之前起阻碍作用；在B点之后起吸附作用。对接终止点C相对对接入口点A，AUV阻力增加2.4%，速度增加1.75%。
研究结果指出了dock对AUV对接阻碍作用是较小的，AUV无需增加额外推力就可以克服阻力航行到对接终端。此外，外部海洋环境会对AUV航行姿态、对接航向和dock开口朝向产生影响，下一步将重点研究，以提高真实海洋环境下AUV的对接成功率。
参考文献
| [1] | 燕奎臣, 吴利红. AUV水下对接关键技术研究[J]. 机器人, 2007, 29(3): 267-273. YAN K C, WU L H. A survey on the key technologies of AUV underwater docking[J]. Robot, 2007, 29(3): 267-273. (in Chinese) |
| [2] | CARRICA P M, CASTRO A M, STERN F. Self-propulsion computations using a speed controller and a discretized propeller with dynamic overset grids[J]. Journal of Marine Science and Technology, 2010, 15: 316-330. |
| [3] | CARRICA P M, FU H P, STERN F. Computations of self-propulsion free to sink and trim and of motions in head waves of the KRISO container ship(KCS)model[J]. Applied Ocean Research, 2011, 33(4): 309-320. |
| [4] | CARRICA P M, HOSSEINI H S, STERN F. CFD analysis of broaching for a model surface combatant with explicit simulation of moving rudders and rotating propellers[J]. Computers & Fluids, 2012, 53: 117-132. |
| [5] | CHASE N, CARRICA P M. Submarine propeller computations and application to self-propulsion of DARPA SUBOFF[J]. Ocean Engineering, 2013, 60: 68-80. |
| [6] | MOFIDI A, CARRICA P M. Simulation of ZigZag maneuvers for a container ship with direct moving rudder and propeller[J]. Computers & Fluids, 2014, 96: 191-203. |
| [7] | 沈志荣.船桨舵相互作用的重叠网格技术数值方法研究[D].上海: 上海交通大学, 2014: 150-162. SHEN Z R.Development of overset grid technique for hull-propeller-rudder interactions[D].Shanghai: Shanghai Jiao Tong University, 2014: 150-162(in Chinese). |
| [8] | 段旭鹏, 常兴华, 张来平. 基于动态混合网格的多体分离数值模拟方法[J]. 空气动力学学报, 2011, 29(4): 447-452. DUAN X P, CHANG X H, ZHANG L P. A CFD and 6DOF coupled solver for multiple moving object problems based on dynamic hybrid grids[J]. Acta Aerodynamica Sinica, 2011, 29(4): 447-452. (in Chinese) |
| [9] | 张玉东, 纪楚群. 包含动边界的非定常流场动网格数值模拟[J]. 计算物理, 2006, 23(2): 166-170. ZHANG Y D, JI C Q. Numerical methods with dynamic structured meshes for unsteady flows with moving boundaries[J]. Chinese Journal of Computational Physics, 2006, 23(2): 166-170. (in Chinese) |
| [10] | 张来平, 邓小刚, 张涵信. 动网格生成技术及非定常计算方法进展综述[J]. 力学进展, 2010, 40(4): 424-446. ZHANG L P, DENG X G, ZHANG H X. Reviews of moving grid generation techniques and numerical methods for unsteady flow[J]. Advances in Mechanics, 2010, 40(4): 424-446. (in Chinese) |
| [11] | 吴利红, 封锡盛, 胡志强. 三维动态混合网格在AUV发射过程中的应用[J]. 船舶力学, 2010, 14(7): 717-722. WU L H, FENG X S, HU Z Q. Application of 3-D hybrid dynamic grids to simulate the flow in AUV swim-out[J]. Journal of Ship Mechanics, 2010, 14(7): 717-722. (in Chinese) |
| [12] | WU L H, LI Y P, SU S J, et al. Hydrodynamic analysis of AUV underwater docking with a cone-shaped DOCK under ocean currents[J]. Ocean Engineering, 2014, 85: 110-126. |
| [13] | ZHAN J M, CAI W H, HU W Q, et al. Numerical study on the six-DOF anchoring process of gravity anchor using a new mesh update strategy[J]. Marine Structures, 2017, 52: 173-187. |
| [14] | 郝继光, 姜毅, 韩书永, 等. 一种新的动网格更新技术及其应用[J]. 弹道学报, 2007, 19(2): 88-92. HAO J G, JIANG Y, HAN S Y, et al. A new dynamic mesh update method and its applications[J]. Journal of Ballistics, 2007, 19(2): 88-92. (in Chinese) |
| [15] | WU L H, LI Y P, ZHANG H C, et al.Meshing impact on numerical simulation of marine systems performance[C]//Proceedings of Oceans 2016.Piscataway, NJ: IEEE Press, 2016: 1-6. |
| [16] | 吴利红, 李一平, 刘开周, 等. 基于多块动态混合网格的AUV自航类物理数值模拟[J]. 机器人, 2019, 41(6): 706-712. WU L H, LI Y P, LIU K Z, et al. Physics-based numerical simulation of AUV self-propulsion using multi-block hybrid dynamic mesh method[J]. Robot, 2019, 41(6): 706-712. (in Chinese) |
| [17] | ROBERT S M, BRETT W H, LANCE M, et al. Docking control system for a 54-cm-diameter(21-in)AUV[J]. IEEE Journal of Oceanic Engineering, 2008, 33(4): 550-562. |
