1 引言
水下爆炸模型试验是研究水下爆炸的重要手段之一。然而目前模型试验的推广还存在一些困难，主要是因为水下爆炸冲击波、气泡脉动、迁移难以满足统一的相似律，以及边界效应的消除问题[1]。
从流体动力学角度来说，刚性边界或弹性边界的存在会阻碍气泡、边界之间的水在气泡径向方向上的运动，库尔[2]把这种阻碍作用定性地表述为“壁面使得流体‘惯性’增大”；而自由面的存在则导致气泡与自由面间水的“惯性”相比自由场要小。这种定性的理论，可以粗略地阐释自由面附近气泡脉动周期缩短而刚性壁面附近气泡周期延长这一现象。除此以外，边界面的存在也是气泡中心发生迁移的原因之一。气泡内压力大于流体静压时，气泡、壁面之间水较难发生位移，气泡中心远离壁面；气泡内压力小于流体静压时，远离壁面的流体加速度远大于气泡、壁面间的流体加速度，气泡中心产生向壁面方向运动的趋势。事实上，在气泡脉动第一个周期内约80%时间里，气泡内的压力都小于流体静压，因而前者效应相对而言十分微小，气泡总体上表现为受到壁面的吸引。与之相对的，因为自由面对气泡-边界之间流体“惯性”的影响与刚性壁面截然相反，所以在自由面附近，不考虑重力的影响，总体来说气泡中心向远离自由面方向迁移。
关于边界对气泡运动的影响，研究人员开展了很多工作。张阿漫等[3]利用高速摄影系统对气泡与自由液面的相互作用进行了试验研究，发现距离自由液面越近，气泡脉周期越小，最大半径越大，且射流形成时间越早、宽度越大。牟金磊等[4]使用MSC.Dytran对常压下自由场、近自由液面、近刚性壁和弹性边界等不同条件下水下爆炸气泡的动态特性进行了仿真计算，结果表明：当爆心距离边界约等于气泡最大半径时，自由表面会使气泡的周期缩短，刚性壁的影响与之相反，弹性边界的影响介于二者之间；与自由面距离大于2倍气泡半径，与弹性结构和刚性壁边界距离大于3倍气泡半径时，水下爆炸气泡可近似认为做球状脉动。Xu Liuyi等[5]使用有限元法对可变形海床与自由水面之间的水下爆炸进行了数值模拟，计算结果表明随着初始装药深度增大，气泡脉动周期延长，气泡膨胀阶段非球状脉动现象愈加明显。接近海床时气泡收缩阶段由于Bjerknes效应产生向下的诱导射流，与气泡在刚性边界附近的现象相似。
前人的研究[2]中，多只考虑不同边界对气泡动态特性的影响，对边界附近气泡运动过程中一系列复杂动力学行为背后的力学机理分析较少。本研究采用LS-DYNA的S-ALE方法进行数值模拟，分别研究单一自由面及单一刚性壁面附近气泡的动态特性，并结合压力场云图分析单一自由面及单一刚性壁面附近气泡运动特性背后的力学机理。
2 数值模型
采用LS-DYNA的S-ALE方法建立三维1/4流场模型，求解水下爆炸气泡运动，水域横截面尺寸为0.25 m×0.25 m。使用*INITIAL_HYDROSTATIC_ALE关键字设置流体静压及重力作用下形成的静水压力梯度[12]，环境压力值通过*CONTROL_ALE关键字中的PREF项设置。在所有工况中，TNT药包质量均为1 g。由于要考虑边界对气泡运动的影响，因此所有水域网格均匀划分，网格尺寸0.25 cm，约为TNT药包半径的1/2。
考虑自由面效应时，自由面以上空气划分为重力方向上变密度的网格，空气高度0.2 m。考虑刚性壁面影响时，设置刚性面处的位移约束。在除了对称面之外的水域其他外表面上施加无反射边界条件。
计算结果以无量纲形式展现，用以无量纲化的基本物理量[13]有：气泡最大半径am为特征长度，气泡初始位置的初始流体静压p0为特征压力，气泡外部流体密度ρ0为特征密度。涉及到的各种物理量的无量纲化因子如表1所示[13]。
表1 相关物理量的无量纲化因子
Table 1 Dimensionless factors of physical quantities
物理量长度压力密度时间速度加速度 无量纲化因子amp0ρ0amρ0p0p0ρ0amp0ρ0
采用表1中的无量纲因子将气泡初始条件无量纲化，本文涉及到的无量纲初始参数有：
其中：γf为初始时刻气泡中心与自由面的距离；γs为初始时刻气泡中心与壁面的距离；δ为浮力参数，用以衡量气泡受到的浮力效应。
TNT炸药采用JWL状态方程表示，其方程的形式如下：
式中：V 为爆轰产物相对体积；E 为爆轰产物单位体积的内能；A 、B 、R1、R2、ω为材料经验参数。TNT材料参数如表2所示[12]。
表2 TNT的材料及状态方程参数
Table 2 Material parameters for the JWL equation
ρ/(kg·m-3)pCJ/GPaDCJ/(m·s-1)ΔQ/(kJ·kg-1)1 57019.16 8504 400A/GPaB/GPaR1R2ω374.23.234.150.950.30
通常采用MAT_NULL材料模型及EOS_LINEAR_POLYNOMIAL多项式状态方程描述空气，形式如下：
p=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E
式中：μ=ρ/ρ0-1；ρ/ρ0为当前密度与初始密度的比值，ρ0在*MAT_NULL中定义；E为空气单位体积的内能；Ci为有量纲的常数。空气材料参数如表3[12]，当设定C0=C1=C2=C3=C6=0以及C4=C5=γ-1时，可描述理想气体：
表3 空气材料及状态方程参数
Table 3 Material parameters for the linear polynomial equation
ρ/(kg·m-3)C0～C3C4C5C6V01.184 500.40.401.0
采用MAT_NULL材料模型及EOS_GRUNEISEN状态方程描述水，其定义材料压力的形式为：
式中：
为介质声速；γ0为GRUNEISEN系数；S1、S2、S3为介质冲击波-粒子速度曲线的斜率系数；a为体积修正系数。水的材料参数如表4[12]。
表4 水的Gruneisen状态方程参数
Table 4 Material parameters for the Gruneisen equation
ρ/(kg·m3)C/(m·s-1)S1S2S3γ0a998.211 6471.921-0.0960.00.350.0
3 计算结果及分析
3.1 近加压自由面气泡运动特性及力学机理
本节通过数值方法来计算加压自由面附近气泡的运动特性，其中自由面压力pa=29.98 atm(水深300 m处流体静压)，气泡由1 g TNT炸药爆炸生成，初始位置距离自由面d=50 mm。于是可得到无量纲距离γf≈0.95，浮力参数δ≈0.013，在此条件下浮力作用可以忽略不计[11]。
给出若干典型时刻气泡及附近流场的压力云图的演化，如图1，从而在力学层面阐述气泡各阶段各种物理现象的成因，明晰气泡动力学行为背后的力学机理。图1中每幅图以初始药包中心为原点，横向坐标范围为-1.90≤r≤1.90；每幅图下的标注为对应的无量纲时刻；无量纲时间因子为955 μs，无量纲长度因子为52.7 mm。
图1 d=50 mm时的无量纲压力场云图
Fig.1 Dimensionless pressure field cloud maps at d=50 mm
在释放出冲击波之后，气泡内部处于高压状态，开始迅速向外膨胀，自由液面随之开始微微隆起。冲击波在自由面处反射，反射稀疏波与衰减后的入射波叠加可能导致自由面附近某些区域发生压力截断现象，部分区域压力甚至降到水的饱和蒸气压以下，促使空化效应的产生。但无论是否考虑自由面处冲击波的反射，气泡外各处水介质的压力相较于此时气泡内的高压来说都很微小，可以认为气泡上方自由面处的流体静压与相同深度无限远处水域中流体静压相差无几，而高压气泡与上方自由面相距最近，因此初生时刻气泡与上方自由液面之间的正向压力梯度最大。这导致相对其他方向来说，气泡顶部膨胀速度最快，气泡中心向自由面不断靠近。如第3幅图所示，此阶段气泡整体呈现卵形，顶部曲率大，其他方向曲率小。
当膨胀到一定程度时，气泡内压力开始低于外部流体压力，惯性效应推动着气泡继续膨胀(过度膨胀阶段)，膨胀速度逐渐减小，直到约0.73时刻，整个流场动能几乎完全转化为气泡外部水介质的势能，气泡整体达到最大体积。但如果从局部来看，约0.63时刻气泡顶部就率先开始收缩，此时气泡其他方向仍在膨胀。这意味着在外界不均匀的静水压力作用下，气泡表面各个方向的膨胀收缩过程并不统一。过度膨胀阶段，气泡中心向自由面方向迁移的速度逐渐减小到0，之后开始反向迁移，气泡中心与自由面距离不断增大。气泡达到最大体积时，形状依然为卵形。
文献0指出：对于坍塌初始阶段而言，后续气泡表面收缩速度与该处的曲率成正比。
事实上，可以将坍塌初始时刻具有不同曲率的气泡表面各处分离开来，想象为不同压力自由场中具有不同最大半径(最大曲率)的球形气泡，那么由气泡脉动周期及最大半径的经验公式可得自由场中气泡表面径向平均速度约为
其中T为气泡脉动周期，J、K、K1为常数，d0为大气压力贡献的压头)，即自由场压力越大，气泡达到最大体积时曲率越高(最大半径越小)，并且气泡表面径向平均速度越大，从而解释非自由场中气泡表面膨胀收缩速度与该处的曲率之间的正比关系。
本文中的仿真结果与以上论述吻合。由第4幅图可见，刚开始收缩时，气泡顶部曲率最高，收缩速度最快，同时气泡与自由面间逐渐形成压力不断提高的局部高压区。由第4～8幅图可见，局部高压区反过来又导致气泡顶部收缩速度进一步加快(水射流的加速)。根据流体力学质量守恒定律，水射流加速导致向自由面-气泡顶部之间的驻点处汇聚的流体速度加快，而此时驻点与局部高压区基本重合，因此汇聚的流体使得局部高压区压力逐渐增大，更大压力梯度的高压区又进一步推动了水射流的加速。
气泡坍塌阶段，水射流不断发展，表现为速度的不断提升、直径的逐渐增大以及射流头部曲率的逐渐减小。水射流于1.55时刻击穿气泡，产生砰击载荷。此后气泡进入环形脉动阶段，气泡剖面形状呈现橘瓣形。约1.60时刻，气泡收缩到最小，在这一时刻前后气泡处于过度压缩状态，释放出脉动压力波，无量纲压力峰值约为17.43。
收缩阶段气泡中心始终向远离自由面方向迁移。气泡第二次脉动初期，气泡中心才略微有朝向自由面的迁移，但是随着气泡内压力迅速减小到周围水介质流体静压以下，气泡中心再次向远离自由面方向迁移。气泡整个运动过程的迁移曲线(以无量纲化)见图3所示。由图3可知，自由面边界对高压气泡具有吸引作用，对低压气泡具有排斥作用。
观察自由面的演化，可以发现整个过程水冢不断升高。气泡初生阶段自由面仅仅略微隆起。0.33时刻，区域空化现象基本消失，后续整个流场的流动完全随着气泡与外部水介质的能量转化过程进行，气泡顶部与初始自由液面基本持平，自由液面隆起曲率增大。接近0.73时刻，气泡达到最大体积。这一时刻前后气泡上方自由面形成突起，可以观察到自由液面与气泡间明显的压力梯度。膨胀阶段自由面的隆起主要是气泡的推动造成的。坍塌阶段，流场其他区域的水不断向着气泡与自由面之间汇聚，靠近气泡的水介质汇聚形成水射流，而靠近自由面的水则汇聚成为水冢。之后局部高压区不断向下移动，但驻点位置移动较小，水冢依然不断升高。环形脉动阶段，水冢直径逐渐减小，水冢根部形成环形的水幕，但是不明显，一般称这种水冢形态为皇冠形[11]。
3.2 水平刚性壁面附近气泡运动特性及力学机理
设置静水压力等效300 m水深，TNT药量1 g，药包中心与刚性壁面距离b=50 mm。因此无量纲浮力参数δ约为0.013，气泡初始位置与刚性壁面的无量纲距离γs约为1.1。
本文在γs≈1.1、δ≈0.013条件下分别进行过气泡处于水平刚性壁面上方、下方2种工况的数值模拟，发现气泡运动特性并无明显差别，即这2种工况中浮力的影响可以忽略不计。以下结果出自气泡处于水平刚性壁面上方的数值模拟工况。气泡的形态演变过程及水中压力场云图如图2，每幅图以初始药包中心为原点，横向坐标范围为-1.67≤r≤1.67；无量纲时间因子为813 μs，无量纲长度因子为44.9 mm。
图2 b=50 mm时的无量纲压力场云图
Fig.2 Dimensionless pressure field cloud maps at b=50 mm
气泡初生时刻呈球形，内部处于高压状态，迅速向外膨胀，气泡内部压力随之迅速减小。在很短的时间内气泡内压力就减小到周围流体静压以下(如第2幅图所示，靠近气泡表面的流体静水压力小于更外部的流体)，此后完全依靠惯性继续膨胀，膨胀速度不断减小。膨胀到最大体积时，气泡形状依然接近球形，只是靠近壁面的一侧由于受到壁面的阻碍作用，较为扁平。由第2～4幅图可知，随着气泡膨胀，气泡与壁面之间的流体逐渐变得狭窄，这一区域逐渐形成气泡-壁面间连通的低压区，而气泡上方流场的压力较高，导致气泡受到朝向壁面的合力，使其向壁面缓缓靠近。这便是刚性壁面吸引气泡的力学机理，一般称这种使得气泡趋近壁面的力为Bjerknes力。在膨胀时，壁面的Bjerknes效应尚不显著。
与自由面附近观察到的现象一致，在外部不均匀的流体静压作用下，气泡表面各处并非同时停止膨胀的：0.86时刻，气泡上端就率先开始收缩，而气泡下端在整个第一次脉动期间完全没有收缩，气泡整体则是于约0.98时刻才达到最大体积。
坍塌阶段初期，气泡表面各处曲率相差不大。由第4～8幅图可见，坍塌前半段，气泡下端在低压区作用下完全不收缩，气泡逐渐变为椭球形，顶部曲率最大。最大的曲率使得气泡顶部收缩速度开始加快，对流场产生扰动最大，导致上方流场中出现范围逐渐缩小、与气泡距离不断减小、压力梯度不断升高的局部高压区。第7幅图，局部高压区与气泡顶部产生接触。第8幅图，气泡顶部出现内凹，形成水射流雏形。虽然水射流出现的很晚，但水射流迅速加速，并于2.35时刻击穿气泡。坍塌时在气泡下方低压区以及气泡上方局部高压区的共同作用下，气泡中心向壁面方向迅速迁移，这就是坍塌阶段更为强烈的壁面吸引效应的成因。
2.40时刻气泡压缩到最小，之后气泡下端逐渐贴壁。气泡中轴线上的水射流不断冲击刚性壁面，并在贴近壁面处向四周扩散，导致水射流形状类似一个中间细两端粗的酒杯。砰击射流不断外溢并挤压靠近刚性壁面的气泡下半部，使得环形气泡下半部不断向外扩张，与此同时气泡中部则出现向内的凹陷，并使环形气泡有撕裂的趋势。水射流的砰击过程中，载荷先迅速升高到水锤压力峰值，之后又逐渐衰减为驻压。
将刚性壁面工况下的气泡中心迁移曲线(已无量纲化)同样绘制在图3中。
图3 气泡中心迁移曲线
Fig.3 Migration curve of bubble center
为了更好地表述边界对气泡的吸引或排斥效应，做了以下处理：对于自由面附近工况，纵轴原点设在气泡初始中心处，正方向竖直向上；对于刚性壁面附近工况，纵轴原点设在气泡初始中心处，正方向竖直向下。于是，曲线各点纵坐标为正表示边界对气泡的吸引作用，纵坐标为负表示边界对气泡的排斥作用。由图可知，刚性壁面边界对高压气泡具有排斥作用，对低压气泡具有吸引作用。
4 结论
本研究采用LS-DYNA的S-ALE方法对自由面与刚性壁面附近气泡的运动进行了数值模拟，通过典型时刻气泡外流场压力云图的演化，分析了不计浮力条件下自由面、刚性壁面附近气泡动力学行为及其背后的力学机理，验证了S-ALE方法在水下爆炸三维数值模拟中的准确性，并取得以下结论：
1)不计浮力条件下，在自由面附近：初始膨胀阶段气泡与上方自由面间压力梯度最大，因此气泡顶部表面膨胀速度相比其他方向较快，气泡整体呈现卵形，顶部曲率较大。坍塌阶段气泡顶部与自由面间逐渐形成局部高压区，该局部高压区促使气泡顶部背离自由面方向水射流的形成。
2)不计浮力条件下，在刚性壁面附近：膨胀阶段气泡与壁面之间形成连通的低压区，气泡上方流场压力较大，这便是刚性壁面对气泡吸引作用的力学机理。坍塌阶段气泡顶部附近流场中形成的局部高压区与气泡下方低压区共同作用，使得气泡向壁面方向的迁移相比膨胀阶段更加剧烈。与自由面附近现象相似的是，该局部高压区直接导致朝向壁面方向水射流的形成。但相比自由面，刚性壁面附近水射流形成时间较晚。
3)自由面边界对高压气泡具有吸引作用，对低压气泡具有排斥作用。刚性壁面边界的作用与之相反，对高压气泡具有排斥作用，对低压气泡具有吸引作用。
[1] 张阿漫,王诗平,彭玉祥,等.水下爆炸与舰船毁伤研究进展[J].中国舰船研究,2019,14(03):1-13.
Zhang A M,Wang S P,Peng Y X,et al.Research progress in underwater explosion and its damage to ship structures[J].Chinese Journal of Ship Research,2019,14(3):1-13.
[2] Cole H R.Underwater explosion[M].Princeton,NJ,US:Princeton University Press,1948.
[3] 张阿漫,王超,王诗平,等.气泡与自由液面相互作用的实验研究[J].物理学报,2012,61(08):300-312.
Zhang A M,Wang C,Wang S P,et al.Experimental study of interaction between bubble and free surface[J].Acta Physica Sinica,2012,61(8):300-312.
[4] 牟金磊,朱石坚,刁爱民,等.边界条件对水下爆炸气泡运动特性的影响分析[J].振动与冲击,2014,33(13):92-97.
Mou J L,Zhu S J,Diao A M,et al.Analysis on the characteristics of UNDEX bubbles under different boundary conditions[J].Journal of Vibration and Shock,2014,33(13):92-97.
[5] Xu Liuyi,Wang Shiping,Liu Yunlong,et al.Numerical simulation on the whole process of an underwater explosion between a deformable seabed and a free surface[J].Ocean Engineering,2021,219:108311.
[6] 王树山,李梅,马峰.爆炸气泡与自由水面相互作用动力学研究[J].物理学报,2014,63(19):239-248.
Wang S S,Li M,Ma F.Dynamics of the interaction between explosion bubble and free surface[J].Acta Physica Sinica,2014,63(19):239-248.
[7] TIAN Zhaoli,LIU Yunlong,ZHANG A’man,et al.Jet development and impact load of underwater explosion bubble on solid wall[J].Applied Ocean Research,2020,95.
[8] Phan T J,Nguyen V T,Park W G.Numerical study on strong nonlinear interactions between spark-generated underwater explosion bubbles and a free surface[J].International Journal of Heat and Mass Transfer,2020,163.
[9] Nguyen A T.The Application of CEL Technique to Simulate the Behavior of an Underwater Explosion Bubble in the Vicinity of a Rigid Wall[J].Applied Mechanics and Materials,2020,4796.
[10] CHEN Yingyu,YAO Xiongliang,CUI Xiongwei,et al.Study on the quantitative law of wall pressure caused by mini-charge underwater explosion bubble[J].Ocean Engineering,2020,213.
[11] 崔杰.近场水下爆炸气泡载荷及对结构毁伤试验研究[D].哈尔滨：哈尔滨工程大学,2013:148-179.
Cui J.Experimental study on underwater explosion bubble loads and damage on the structure nearby[D].Harbin:Harbin Engineering University,2013:148-179.
[12] 安丰江,施惠基,刘琼,等.水下爆炸近场压力特性及其与结构的流固耦合作用[J].兵工学报,2015,36(S1):13-24.
An F J,Shi H J,Liu Q,et al.Load characteristics of near-field of underwater explosion and its fluid-structure interaction[J].Acta Armamentarii,2015,36(S1):13-24.
[13] 李帅,张阿漫,韩蕊.水中高压脉动气泡水射流形成机理及载荷特性研究[J].力学学报,2019,51(06):1666-1681.
Li S,Zhang A M,Han R.The mechanism of jetting behaviors of an oscillating bubble[J].Chinese Journal of Theoretical and Applied Mechanics,2019,51(6):1666-1681.
[14] Lauterborn W.Cavitation bubble dynamics — new tools for an intricate problem[Z].Springer Netherlands,1982.