自主水下航行器(autonomous underwater vehicle，AUV)作为可以在水下自主航行的移动节点，主要完成监视侦察、海洋环境监测、通讯导航节点、反水雷、反潜作战、探测识别、载荷投放、时敏打击等任务，无论在民用领域还是在军用领域都具有重要的作用，尤其是在未来的深海对抗中将发挥至关重要的作用[1, 2]。但随着科技的进步以及军民等方面的实际需要，考虑到远距离、长时效等因素，研制带有多个载荷的航行器技术已经成为未来的发展方向。这种航行器是由载荷和运载段组成的一个多体系统。当运载段将载荷运至指定区域后，载荷从顶端逐级分离。这一思想对于新型运载技术的研究具有重要的现实意义和参考价值，多个声呐的布放，大间距传感器的安置，1个运载段携带多个载荷1次完成任务可以大大缩短工作时间，减小工程成本等优点[3]。多载荷AUV在分离运动时，通过进水以增加其负浮力，从而实现分离。如此一来，分离过程中的载荷段及剩余推进段都会受到由于进水引起的附加外力干扰，稍有不慎就存在分离后的载荷和运载体发生撞击而导致分离失败的可能。因此，研究不同工况下多载荷AUV的分离进水过程对其运动姿态以及安全、高效分离具有重要意义。
多载荷AUV的分离进水过程中，水在压力作用下将空腔内空气排出，属于典型的气液多相流问题。李书斌[4]数值分析了水利工程立轴漩涡的进水孔对进水量、进水速度等影响。马斗[5]采用多相流技术对气泡在液体中上升的过程进行了数值模拟，计算结果与试验结果吻合较好。徐玲君等[6]对不同直径气泡在静水中运动特性进行了数值分析与试验研究。上述研究均针对气泡在静水中的运动特性进行研究，与本文中所研究的动态进水过程不同。基于此，本文建立了适合描述分离过程的气液两相流数学模型，利用计算流体力学技术对后平衡舱腔在静止情况下的进水过程进行数值分析，采用大涡模拟(LES,large eddy simulation)与流体体积函数法(VOF,volume of fluid)相结合的方法研究进水孔直径、位置以及航行深度对进水过程的影响规律，得到进水体积随时间的变化过程，避免加工模型水池试验所带来的周期长、费用高等问题，相关研究结论可用于多载荷AUV的设计。
1 控制方程 1.1 VOF方法多载荷AUV的进水过程属于涉及空气、水的气液两相流问题，在进水过程中会出现漩涡、气泡的形成及脱离等现象，过程较为复杂。需要采用合适的界面追踪方法对气液界面进行捕捉。目前已提出多种追踪自由边界的方法，按其基本思想归结起来主要有以下几种：静压假定、刚盖假定、标高函数法、MAC法和VOF法[4]。
由于VOF方法追踪的是网格中的流体体积，不是追踪流体质点的运动，因而具有容易实现、计算量小和精度高等优点。同时该方法能清晰捕捉界面、有效地预测气泡的形状、速度和轨迹[5, 6]，本文采用该方法对多载荷AUV的进水过程进行模拟。
VOF方法通过求解单独的动量方程和处理穿过区域的每一流体的体积分数来模拟2种或2种以上不相混合的流体流动。两相间的界面通过计算每个网格中第k相所占的体积分数φk来追踪：
且在任意网格中存在：
φk通过连续性方程来求解，即：
在多载荷AUV进水过程的两相混合流动中，动量方程对各个流体相均适用，其速度场通过求解Navier-Stokes方程获得：
式中，p为静态压力;u=(u,v,w)为流体的速度场；F为动量源项，由表面张力产生；ρ和μ分别代表容积平均密度和动力黏度：
式中，下标g和l分别表示气相和液相。
动量源项的处理采用CSF(continuum surface force)模型[7]：
计算区域内界面重构和界面推进采用分段性重构方法，界面的位置由界面法向量n和定义在网格中心的体积分数φk确定，并采用拉格朗日方法跟踪随流动传播的界面。界面曲率κk为：
式中，n=∇φk。
1.2 大涡模拟方法大涡模拟作为湍流基础性研究的高级模拟手段，已随计算机硬件的发展逐渐成为工业湍流流动问题计算的可选方法之一。经各向同性滤波器滤波后的大涡模拟基本控制方程可表示为[8]：
式中，上划线代表经过空间滤波的变量；ui为速度；p为压力；τij是一个由非线性项产生的未知量，被定义为亚格子尺度应力，表征过滤掉的小尺度脉动和可解尺度湍流间的动量输运：
采用Smagorinsky涡黏模型来模化亚格子应力，其张量中的偏量部分为
式中，vt为涡黏系数；Cs为Smagorinsky常数，Δ为滤波尺度，通常取0.1～0.23，sij为经过滤波后的速度变形张量，定义如下
2 物理模型和计算方法 2.1 物理模型
以无推力推动的多载荷AUV为研究对象，该模型具有前、后2个平衡舱段和中间载荷段3个部分。其中前平衡舱在航行初期通过进水来增加多载荷AUV的重力，依靠重心位置的变化使其在航行中逐步旋转下潜至分离方位；然后，后平衡舱开始进水，通过重力的增加和重心的偏移，后平衡舱在重力作用下分离将载荷释放，载荷脱离后自由航行到达预定位置，完成既定任务。
由于前、后平衡舱的进水过程具有相似性且主要依靠后平衡舱的脱离释放载荷，以后平衡舱为研究对象，如图 1a)所示，该模型为圆柱形回转体，长L=1 440 mm，最大直径D=534.4 mm，后平衡舱空腔的体积V=0.554 2mm。图 1b)为模型的轴对称面，如图所示模型上有2个进水孔，分别为进水孔1和进水孔2，2个进水孔直径相等且进水孔的中心线与模型的轴线在同一平面上。通过数值仿真展示其进水全过程，重点研究进水孔直径、进水孔位置以及航行深度对进水过程的影响，为多载荷AUV的安全快速分离提供技术支持。
 |
| 图 1 计算模型 |
采用有限体积法，压力速度的耦合采用PISO算法。为了得到足够高的网格分辨率来精确模拟进水过程，计算域全部采用结构化多块网格划分的思路，所有的变量都在网格中心给出，最终生成了高质量的内外流场网格，网格总数为145万左右，如图 2所示。
 |
| 图 2 计算域网格 |
图 3所示的模拟区域各个面都采用对称平面来处理，底面给定无滑移边界条件，上表面设为压力出口，整个计算区域做成一个开放系统。因此，在水体上的气相空间为大气压强，计算中方程的对流项采用混合差分格式，对时间积分采用迎风差分格式，利用geo-reconstruct方法进行网格重构。通过试算发现，当时间步长取为1×10-4 s时既能获得较好的收敛结果，计算量也不会过大。
 |
| 图 3 边界条件设置 |
为了详细研究进水孔直径、进水孔位置以及航行深度对进水过程的影响，计算生成了一系列后平衡舱空腔进水过程的进水云图，通过计算模拟的方式展示了后平衡舱空腔进水的全过程。图 4展示了不同时刻下的计算结果，从图中可以看出，水从进水孔1沿着平衡舱底部进入后平衡舱空腔，在压力的作用下将空气从进水孔2排出。排出的气泡呈球冠形生长，气泡在初始时刻，内外压强保持平衡，保持气泡形状的稳定。随着气体的不断排出，在表面张力作用下孔口附近三相接触线出现外扩现象，气泡越来越大。随着气体体积继续增加，气泡受到的剪切力以及惯性力也越来越大，气泡不再呈规则的圆球形长大，而是出现凹凸不平的非球形形状，在浮力和绕流曳力共同作用下气泡缩颈断裂，最后脱离孔口。
 |
| 图 4 空腔内空气排出的变化过程 |
进水孔的直径对平衡舱空腔的进水过程有着重要的影响。图 5为水深10 m、L1=150 mm和L2=850 mm时，进水孔直径分别为60 mm、70 mm和80 mm时进水体积随时间的变化曲线。从图中可以看出，进水孔直径对进水速率产生显著影响。直径越大，平衡舱的进水速率越大，但是到达气液两相平衡状态的时间基本相等，均在0.768 s左右。
 |
| 图 5 不同直径时进水体积随时间变化曲线 |
为了便于分析，图 6给出了在1 s时不同直径下空腔内空气排出状态的气液两相云图(轴对称面)。从图中可以看出，当进水孔直径较小时，形成的气泡较小，气泡在轴向上所受的剪切力及惯性力也较小，因此在该方向上气泡的倾斜角度很小。在气泡长大过程中，一方面由于轴向上的剪切与挤压作用，另一方面由于铅垂方向上重力的作用，使得气泡在铅垂方向慢慢上升，连接气泡与孔口底部的气泡颈部越来越细，最后气泡从此处断开，脱离孔口。对于较大直径的进水孔，连接孔口与气泡主体的气泡颈部较粗，需排出大量气体以增大气泡表面积来提高轴向上气泡所受到的剪切力及惯性力，增加了气泡脱离孔口的时间。直径太小，影响空腔的进水速率；直径太大，气泡脱离孔口的时间较长，易引起孔口堵塞现象。这表明，在实际的进水孔设计中，需要根据进水速率以及气泡脱离孔口时间的需求，选择合适的孔口直径。
 |
| 图 6 1 s时不同直径的气液两相云图 |
由图 1可见，在静水压力作用下，水从进水孔1进入平衡舱空腔，将平衡舱空腔内的气体从进水孔2排出，进水孔的相对位置对平衡舱空腔的进水过程有影响。图 7为水深10 m、L1=150 mm和Φ=80 mm时，进水孔2与进水孔1水平间距L2分别为190 mm、520 mm和850 mm时的进水体积随时间的变化曲线。从图 7中可以看出，水平间距对平衡舱进水速率及到达平衡状态时的平衡舱进水体积均有影响。水平间距越小，平衡舱进水体积随时间变化曲线斜率越大，即平衡舱的进水速率越大。同时，到达气液两相平衡状态时空腔内进水体积越大，其进水体积分别为99.94%、98.09%和94.41%，但到达气液两相平衡状态的时间基本相等。
 |
| 图 7 不同水平间距时进水体积随时间变化曲线 |
为了便于分析，图 8为在0.5 s时不同水平间距下空腔内空气排出状态的气液两相云图(轴对称面)。从图中可以看出，排出的气体在表面张力下维持气泡的状态，在压力梯度作用下气泡向上运动，气泡在运动过程中，内部压力保持恒定。进水孔水平间距越小，连接孔口与气泡主体的气泡颈部较粗，气泡脱离孔口所需时间较长。随着时间的推进，在空腔内外压强作用下，气液两相平衡状态下空腔内残留的空气体积与进水孔2的位置相关。由于平衡舱尾部为一收缩段，进水孔的水平间距越大，进水孔2的垂向距离越小，使得残留在空腔内的空气越多。在实际的进水孔设计中，需要根据使用条件合理设计进水孔的相对位置。
 |
| 图 8 0.5 s时不同水平间距的气液两相云图 |
水深越深平衡舱承受的静水压力越大，而静水压力影响着空腔的进水过程。为了分析水深对进水过程的影响规律，数值计算了在Φ=80 mm、L1=150 mm和L2=850 mm时，水深分别为10 m、20 m和30 m时的进水过程。图 9为不同水深时进水体积随时间的变化曲线。
 |
| 图 9 不同水深时进水体积随时间变化曲线 |
从图中可以看出，水深越大，平衡舱的进水速率越大，到达气液两相平衡状态的时间越短，其达到平衡的时间分别为0.808 s、0.69 s和0.624 s。这表明，增加水深可以加快水从进水孔1流入平衡舱空腔的速度，进而加剧了空腔内空气的排出速度。
4 结 论本文运用LES与VOF相结合的方法对多载荷AUV后平衡舱空腔在静止情况下的进水过程进行了数值模拟，研究了不同进水孔直径、进水孔位置和航行深度对进水过程的影响规律，得到了不同工况下进水体积及质量随时间的变化过程，深入分析结果可以得到以下结论：
1)进水孔直径越大，进水速率越大，但气泡脱离孔口所需的时间越长，且生成的气泡越大；
2)进水孔水平间距越小，由于尾部收缩段的存在，在内外压强作用下使得进水速率增大，且到达气液两相平衡状态时空腔内进水体积越大；
3)航行深度越大，承受的静水压力越大使得后平衡舱空腔到达气液两相平衡状态的时间越短
| [1] | The Navy Unmanned Undersea Vehicle (UUV) Master Plan[R]. Assistant Secretary of the Navy for Research Development & Acquisition, 2004 |
| [2] | FY2011-2036 Unmanned Systems Integrated Roadmap[R]. United States of America, Department of Defense, 2011, 11-S-3613 |
| [3] | 潘光, 施瑶, 杜晓旭, 杨智栋. 凸台对多载荷AUV运载段阻力和流噪声的影响[J]. 西北工业大学学报, 2013, 31(6): 841-847 Pan Guang, Shi Yao, Du Xiaoxu, Yang Zhidong. Exploring Effect of Carrier Appendage on Drag and Flow Noise of Carrier of Multi-Load AUV[J]. Journal of Northwestern Polytechnical University, 2013, 31(6): 841-847 (in Chinese) |
| Cited By in Cnki (1) | Click to display the text | |
| [4] | 李书斌. 进水口漩涡影响因素研究[D]. 天津: 天津大学,2009 Li Shubin. Research on the Influencing Factors of Vortex at Hydraulic Intakes[D]. Tianjin: Tianjin University, 2009 (in Chinese) |
| Cited By in Cnki (2) | |
| [5] | 马斗. 单孔气泡动力学行为的VOF数值模拟[D]. 天津:天津大学,2009 Ma Dou. Numerical Simulation of Single Orifice Gas Bubble Motion Behavior with Volume of Fluid Method[D]. Tianjin: Tianjin University, 2009 (in Chinese) |
| Cited By in Cnki (7) | |
| [6] | 徐玲君, 陈刚, 邵建斌, 薛阳. 不同直径气泡在静水中运动特性的研究[J]. 水动力学研究与进展A辑, 2012, 27(5): 582-588 Xu Lingjun, Chen Gang, Shao Jianbin, Xue Yang. Study on Motion Behavior of Bubbles with Different Diameters in Still Water[J]. Chinese Journal of Hydrodynamics, 2012, 27(5): 582-588 (in Chinese) |
| Cited By in Cnki (4) | Click to display the text | |
| [7] | Brackbill J U, Kothe D B, Zemach C. A Continuum Method for Modeling Surface tension[J]. Journal of Computational Physics, 1992, 100(2): 335-354 |
| Click to display the text | |
| [8] | Smagorinsky J. General Circulation Experiments with the Primitive Equations[J]. Monthly Weather Review, 1963, 91(3): 99-164 |
| Click to display the text |