基于数值仿真的球体垂直入水空泡演化过程研究

2024-03-05张晨星谷立祥权晓波魏海鹏程少华

导弹与航天运载技术 2024年1期

张晨星，谷立祥，权晓波，魏海鹏，程少华

（北京宇航系统工程研究所，北京，100076）

0 引言

物体由空气介质穿越气水界面运动至水中的过程称为入水过程。最早对物体入水开展研究的是英国学者Worthington A M和Cole R S［1］，他们通过高速闪光相机对球体落入液体中所产生的喷溅进行了拍摄，得到了大量的球体入水图像。此后，随着鱼雷空投入水、水上飞机水上着陆、火箭助推器以及航天器海面回收等工程问题的出现，大量的学者对入水问题开展了相应的研究。

作为最早得到研究的入水问题，球体入水一直受到学者们的关注。球体入水的过程往往伴随着入水空泡的产生和演化。入水空泡的演化过程涉及到固液冲击和多相流动，是球体入水问题的一项重要研究内容，国内外学者均对球体入水空泡开展了相应的研究工作。

早期对球体入水空泡的研究主要以试验和理论为主，Gilbarg D 和Anderson R A［2］通过试验研究了液面压力和入水速度对球体入水空泡闭合时间的影响，发现液面压力越高，入水速度越快，球体入水空泡表面闭合的速度就越快；May A［3］研究了球体表面状态对入水空泡的影响，发现球体表面污染越严重、表面沾有的液体黏性越强，球体入水产生空泡所需的速度越小；此外，May A［4］在理论和试验的基础上，提出了适用于不同外形物体入水空泡早期形态的计算公式；Duez 等［5］提出了用于判断物体入水是否产生空泡的依据，指出物体入水空泡的产生受到入水速度和物体表面亲/疏水性两个因素的影响；在Duez 研究的基础上，Truscott T T［6］开展了旋转球体的垂直入水试验，发现除了入水速度和球体表面亲/疏水性外，旋转也会对空泡的形成产生影响；马庆鹏等［7］对球体垂直入水开展了试验研究，得到了球体入水空泡的演化图像，分析了入水速度和表面沾湿状态对入水空泡流动的影响。

在20 世纪80 年代以后，随着计算机技术的发展和流体力学数值计算方法的成熟，数值仿真成为研究球体入水问题新的有效方法。Lee M等［8］为验证其提出的高速入水空泡动力学模型，对球体的高速入水过程进行了数值模拟，得到了球体高速入水时的空泡形态；张伟伟等［9］采用多物质ALE 算法对球体的入水空泡进行了计算，发现该算法能够较好地对球体入水过程进行仿真模拟；余彧等［10］研究了不同的计算参数对亲/疏水性球体入水过程仿真结果的影响；孙钊等［11］对不同亲/疏水性球体的入水空泡进行了仿真计算，得到了4种不同的空泡形态。

综上所述，在对球体入水空泡的研究方面，国内外学者在试验、理论和数值仿真方面均取得了一定的成果，但目前对球体入水空泡的数值仿真研究多集中于对算法的研究以及对不同因素对入水空泡的影响规律的研究，对入水空泡演化过程的详细分析研究较少。而试验研究会受到试验条件和测量设备的限制，获取的球体入水流场信息以图像和压力数据为主，难以对入水空泡的演化进行细致的分析。因此，本文采用数值仿真的方法对球体垂直入水过程进行了仿真计算，基于仿真结果对球体入水空泡的演化过程进行全面系统地分析。

1 数学模型

1.1 基本控制方程

球体入水涉及到多相流动，因此本文采用均质多相流的N-S 方程作为流动的基本控制方程。多相流模型选取VOF 模型，即将多相流体视为单一的流体介质混合物，各相具有相同的压力场和速度场。αl和αg分别表示液相和气相的体积分数，两者满足关系式：

式中ρm为混合介质的密度，且ρm=αlρl+αgρg，ρl和ρg分别为液相和气相的密度；ui为速度分量。

动量守恒方程为

式中μm为混合介质的黏性系数，且μm=αlμl+αgμg，μl和μg分别为液相和气相的黏性系数；μt为湍流黏性系数。

在连续性方程和动量守恒方程的基础上，VOF模型增加了第n相体积分数的输运方程：

通过求解该方程得到控制单元内第n相的体积分数αn来表征该流体单元内的流体相分布。

1.2 湍流模型

本文采用标准k-ε湍流模型，其湍流动能方程为

式中k为湍流动能；ε为湍流能量耗散率；Gk和Gb为湍流生成项；YM为湍流耗散项；Sk和Sε为源项，C1ε=1.44，C2ε=1.92，σk= 1.0，σε=1.3。

1.3 网格划分及边界条件设置

Aristoff J M 等［12］对直径为25.4 mm，不同密度比的小球以2.17 m/s 左右速度的入水过程开展试验研究，本文采用二维轴对称模型对其试验中密度比为2.3的球体的入水过程展开仿真计算。如图1所示，空气域的高度L2=160 mm，水域高度L3=400 mm，径向直径L1=400 mm。球体壁面第一层网格高度取为0.002 mm，网格均采用结构化网格且在球体附近加密，网格数量为158 400。入口条件为压力入口，入口压力设置为标准大气压强p0=101 325 Pa，出口条件为压力出口的压力为p=p0+ρlgh，其中g为重力加速度，h为水深。本文中压力出口处水深h=L3=400 mm，因此压力出口处的压力设置为105 242 Pa。球体表面和计算域侧面定义为壁面条件。计算初始时刻，球体距水面1 mm，设置球体的初始速度为2.17 m/s。

图1 计算域、边界条件及网格划分示意Fig.1 Schematic of computational area, boundary conditions and mesh

为方便对球体的入水过程进行描述，如图1 所示，记小球下落的方向为轴向x，记垂直于下落的方向为径向y，对称轴和自由水平液面的交点为坐标原点。定义角度θ为球心和球面一点连线与轴向的夹角，逆时针为正。此外，入水深度为球心距水面的距离，取球心到达水面的时刻为零时刻。

1.4 数值计算方法

本文应用流体计算软件对球体入水过程进行计算。采用有限体积法对控制方程进行离散，采用SⅠMPLE算法对压力-速度场的耦合进行求解。压力场的离散选用PRESTO格式，体积分数的离散采用Geo-Reconstruct格式，动量方程和能量方程的离散采用二阶迎风格式。由于入水过程中球体在计算域中发生运动，因此采用动网格技术来实现网格的移动和更新。同时，编写用户自定义函数来计算球体的运动速度和位移。

1.5 仿真结果验证

如图2所示，本文通过对仿真得到的入水空泡形态和文献［2］的试验结果进行对比来验证数值仿真的正确性。可以看到，仿真较好地模拟出了球体入水空泡的扩张、收缩以及闭合等演化过程。仿真得到的空泡形态在球体入水初期和试验结果的一致性较好。仿真中空泡的闭合时间相比试验提前了约6.7 ms，存在10.14%的时间误差。出现误差的原因可能包括网格离散引起计算误差、试验中对高速摄像进行时间标定时存在测量误差以及仿真模型中球体的物理特性与试验中不完全相同等。

图2 球体入水空泡形态试验和仿真结果对比Fig.2 Comparision of experimental and numerical results of cavity's shape which was formed during sphere's waterentry

球体的表面特性是球体的一项重要物理特性，球体表面出现沾湿、沾灰等污染现象都会影响球体的表面特性，即会改变球体表面的润湿性［11］。表面接触角是表征球体表面润湿性的一个常用物理量，在Aristoff J M 试验中，小球的表面接触角为（120±5）°，考虑球体表面受到污染，取球体表面接触角为150°，对球体的入水过程进行仿真计算，得到球体入水空泡的闭合时间为62.7 ms，空泡闭合时间的误差降至5.14%。可见球体表面的污染特性会对球体入水空泡的闭合时间产生一定的影响。

针对Aristoff J M 的球体入水试验，张伟伟等［9］采用多物质ALE仿真算法进行了仿真模拟，重点分析了球体密度和入水速度对入水空泡闭合时间的影响。得到的空泡闭合时间与试验值也存在约10%的偏差，可见实现与试验结果高度吻合的仿真较为困难。

2 球体入水空泡演化过程分析

2.1 入水空泡演化过程阶段划分

图3 给出了球体入水过程中空气体积分数的云图。可以看到，在球体撞击水面之后，球体底部的水体被排开，球体底部被沾湿。部分水体沿着球体表面的切向斜向上涌起，形成高于液面的射流。随着球体的下落，球体的沾湿面积迅速增大，直至水体从球体表面分离。水体从球体表面分离之后，球体后方的空气源源不断地进入水面以下，形成一个开口的空泡，通常称之为入水空泡。空泡形成之后，首先随着球体的下落被不断地拉伸，同时向径向扩张，空泡呈现为倒圆台的形态。随着球体入水深度的不断增加，水面和球体中间的部分空泡截面开始由扩张变为收缩，收缩最快的空泡截面率先收缩为一点，将入水空泡分割为两个部分：上方的敞口空泡和下方附着在球体上的尾空泡。这种现象通常称为空泡的深度闭合，空泡闭合时，上方的敞口空泡呈现为倒锥形，下方的尾空泡呈现为锥形，且尾空泡的体积明显小于敞口空泡。空泡的深度闭合阻断了空气进入下方的尾空泡，同时产生上下两股高速射流。向上射流的上涌使得敞口空泡逐渐消失，向下的射流则击打在球体表面，使得水体从球体的上表面沾湿球体。随着球体进一步下落，向上的射流涌出水面，而尾空泡则不断被压缩，泡内气体不断泄漏，空泡溃灭。

图3 球体入水空泡演化过程示意Fig.3 Schematic of the evolution of cavity formed during sphere's water-entry

根据球体沾湿面积的变化规律对球体入水初期的空泡演化过程进行阶段划分。球体的沾湿角度随时间的变化曲线如图4所示，这里的沾湿角度θw是指固液接触线处的角度值。记球体入水至沾湿角度快速增大的阶段为冲击及流动形成阶段，可以看到在此阶段球体沾湿角度在6.6 ms的时间内由0°增加至约78.5°。流动形成之后，球体沾湿面积增大的速度明显放缓，空泡壁面各处沿径向均在不断扩张，这个阶段即为空泡的发展阶段。而在空泡深度闭合阶段，空泡壁面则部分沿径向扩张，部分向内收缩，直至壁面某处在空间收缩为一点。这两个阶段是空泡发展演化最重要的阶段，即空泡发展及深度闭合阶段。空泡闭合后，敞口空泡的长度不断减小并最终消失。尾空泡内则不断被压缩，气体不断泄漏，直至球体表面完全浸湿，这个阶段即为空泡的溃灭阶段。

图4 沾湿角度θw随时间变化曲线Fig.4 Curve of wet angle θw versus time

2.2 冲击及流动形成阶段空泡演化分析

入水冲击是指球体与水接触，球体压缩水体产生一个高于水中声速的冲击波的过程，这个过程时间很短，通常为微秒级［14］，球体会受到极大的冲击载荷。冲击及流动形成阶段流场的压力云图和气水界面图如图5所示。可以发现，球体经历入水冲击后，球体沾湿部分下方的流场形成高压区。随着时间的推移，球体的沾湿面积不断增大，流场的高压区在不断地扩展，同时高压区的压力在逐渐降低。

图5 冲击及流动形成阶段流场压力云图及气水界面图Fig.5 Contours of pressure and air-water interface of the flow field during impact and flow formation stage

进一步对作用在球体表面上的压力进行分析，不同时刻球体表面的压力随θw的变化曲线如图6 所示，图中虚线依次表示不同时刻球体的沾湿角度。可以看到，在入水初期的某一固定时刻，球体表面的压力分布规律如下：球体正下方出现高压，且高压向球体外侧不断拓展、增大，在固液接触线附近达到了最大值。压力出现上述分布的原因在于，入水初期，沾湿部分的水体在冲击作用下开始流动，而在固液接触线处的水体流动则尚未形成，因此与球体的冲击作用最为剧烈，从而压力在此处最高。随着时间推移，球体的沾湿面积不断增大，球体附近水体的流动逐渐稳定，球体受到的冲击力逐渐转变为稳定的流体阻力，因此球体下方的压力不断减小，同时最大压力的位置逐渐回归到球体的正下方。

图6 不同时刻球体表面无量纲压力随角度变化曲线Fig.6 Curve of dimensionless pressure versus θw at different times

图7为冲击及流动形成阶段流场的速度云图和气水界面图。可以发现，在球体临近入水时刻，球体下方的空气被高速排出。随着球体的入水，固液接触线处的水体在冲击高压产生的压差力的作用下沿着球体表面迅速向外延伸，形成高速射流。同时，球体的动能不断转化为水的动能，当球体表面的水体获得足够大的径向速度时，水体脱离球体表面开始形成入水空泡，流动也逐渐稳定。

2.3 空泡发展及深度闭合阶段空泡演化分析

空泡发展阶段t=17.15 ms 时流场的压力云图、速度云图、流线图和气水界面图如图8所示。由压力云图可知，由于开口空泡和水面大气相连通，因此泡内的压力接近大气压力。而空泡外部的压力则沿径向逐渐增加至水体静压，因此在空泡发展阶段，空泡壁面会受到周围水的压力，且水深越深，压力就越大。同时可以看到在固液接触线的下方出现了低压区。而由速度云图及流线图可以发现，随着球体的下落，球体后方的空气会高速进入空泡内部，泡内空气的流动速度要明显大于空泡周围水体的流动速度。

图8 空泡发展阶段流场示意Fig.8 Schematic of the flow field during the cavity development stage

由于空泡某一截面的径向位置是径向速度对时间的积分，因此空泡壁面的径向速度是影响入水空泡形态的重要参数。空泡壁面上流体的径向速度随轴向位置和时间的变化曲线如图9所示。

图9 空泡发展阶段空泡壁面径向速度变化Fig.9 Curve of radial velocity of the cavity's wall versus time and radial position during the cavity development stage

由图9a可知，对于某一固定空泡截面，空泡壁面的径向速度在球体经过时迅速增大，而后空泡壁面的径向速度在压差力的作用下经历了一个缓慢的衰减过程。轴向位置越大，即水深越深的截面处，由于球体经过该位置时的速度更低，且压差力也更大，因此该处空泡壁面能达到的最大径向速度越小，且速度衰减得越快。

由图9b可知，对于某一固定时刻，在空泡发展阶段的早期，空泡壁面的径向速度沿着轴向降低，球体表面处空泡的径向速度最小。但在空泡发展阶段的后期，空泡壁面最小径向速度的位置出现在了水面和球体之间的一点。这是由于该点上方的空泡壁面受到的压差力相对较小，因此径向速度衰减较慢。而该点下方的空泡壁面虽然受到较大的压差力，但由于其仍受到球体运动的影响，因此径向速度衰减也相对较慢。

在空泡发展阶段，空泡壁面的径向速度均为正值，空泡壁面各点均保持扩张。而一旦空泡壁面的径向速度出现负值，则表明空泡出现了部分收缩，空泡演化进入深度闭合阶段。

空泡深度闭合阶段流场的压力云图及气水界面图如图10 所示。可以看出，在深度闭合阶段的初期，流场的压力场特征和空泡发展阶段相类似，空泡内的压力和大气压接近。随着空泡的进一步收缩，收缩位置下方的空泡受到压缩，泡内的压力开始上升。空泡闭合的瞬间，径向运动的水体在闭合点处撞击，形成高压区。

图10 空泡深度闭合阶段流场的压力云图及气水界面图Fig.10 Contour of pressure and air-water interface of the flow field during the cavity deep closure stage

图11为空泡在深度闭合阶段的速度云图、流线图及气水界面图。由图11b和11c可知，空泡开始收缩后，收缩位置的空气在水体的挤压下，向上、向下排出。由于空气向下排出时空间有限，会造成下方空泡的压力上升，流动受阻，因此下方空泡内的空气流速较低。而空气向上排出时，由于无空间限制，因此流动较为顺畅，流速也更高。同时，在空泡收缩的位置，存在速度滞止的区域。空泡闭合时如图11d所示，冲击产生的高压驱使水体向上、向下运动，形成高速射流。

图11 空泡深度闭合阶段流场速度云图、流线图及气水界面图Fig.11 Contour of velocity, streamlines and air-water interface of the flow field during the cavity deep closure stage

对空泡形态的演化进行进一步分析，提取深度闭合阶段空泡壁面的径向速度随轴向位置的变化曲线如图12 所示。可以看出，在深度闭合阶段初期，空泡中间偏向球体位置处的壁面径向速度在周围水体压力的作用下开始减小为负值，空泡在此处收缩。但收缩位置上下方的空泡壁面仍保持扩张，这种现象通常称为空泡的“颈缩”。随着空泡闭合的进行，空泡壁面越来越多的部分开始处于收缩状态，甚至刚刚脱离球体表面的空泡壁面也在收缩。这表明相比于球体向水体转移的动能，压力已在空泡的发展中起主要的影响作用。此外还可以发现，水面处的空泡壁面一直在保持扩张，这是由于水面的压力过小。

图12 空泡深度闭合阶段空泡壁面径向速度随轴向位置变化曲线Fig.12 Curve of radial velocity of the cavity's wall versus radial position during the cavity deep closure stage

2.4 空泡溃灭阶段空泡演化分析

空泡溃灭阶段流场的压力云图和气水界面图如图13所示。由图13a可知，空泡闭合之后，闭合点上下方的空泡壁面继续收缩，闭合点扩展为闭合线。径向运动的水体在闭合位置继续相撞形成高压，高压区沿着闭合线向上、向下扩展，压力等值线呈现“8”字形。提取t=59.75 ms 时刻空泡壁面的径向速度，发现闭合线上方的径向速度的量值要比下方的大，因此敞口空泡底部的冲击更为剧烈，因而敞口空泡底部的冲击高压区要比尾空泡上方的高压区大。