基于AUTODYN的气泡与固定壁面相互作用数值模拟

张 伟 岳永威 张阿漫 孙龙泉

1 中国舰船研究设计中心，湖北武汉 430064

2 哈尔滨工程大学船舶工程学院，黑龙江哈尔滨 150001

0 引 言

水下爆炸气泡引起的结构破坏可分为3种：爆炸气泡脉动激发船体梁总体振动，造成整体失稳甚至断裂失效；远场爆炸时，气泡脉动引起舰船上较敏感设备的共振，造成设备破坏；当炸药近场爆炸时，气泡受舰船结构边界的影响，形成冲击射流，造成舰艇结构局部损伤。第3种情况属气泡近壁面运动规律问题，进行理论研究的依据主要是以势流理论建立的水平及垂直刚性面附近在浮力作用下运动的气泡理论模型。该模型基本能反映水下爆炸气泡和周围流体介质的运动规律，但其忽略了边界对气泡形状的影响，较适于远场气泡脉动分析。在试验研究方面，关于水下爆炸气泡对结构的毁伤作用试验研究多采用规则结构或缩比模型，鲜有实船试验。

近年来，由水下爆炸引起的气泡动力学问题成为海军舰船生命力技术领域关注的重点。然而，水下爆炸气泡从形成、膨胀到最终溃灭是一个复杂的物理演化过程，尤其是气泡在运动过程中与周围结构的作用受许多因素的影响，研究难度较大。目前，我国学者主要是以高速摄像的方法对电火花生成的气泡进行观测，进而对气泡的运动规律进行研究［1-2］，而基于数值平台的仿真模拟则开展得相对较少。

AUTODYN［3］是一种显式有限元分析程序，主要用于解决固体、流体、气体及其相互作用的高度非线性动力学问题。目前，我国关于AUTODYN软件在水下爆炸领域的应用与其他商业软件比相对较少。就水下爆炸［4-6］而言，由于爆轰产物和水都属于流体，所以相对于Lagrange描述方法，在AUTODYN中采用Euler方法描述更为方便、有效。数值模拟的过程就是对守恒方程，包括连续性方程、动量方程及能量方程的离散，单纯的守恒方程无法求解动力学问题，必须与状态方程（EOS）联立，才能构成封闭方程组，继而对流体动力学问题进行求解。

1 有效性验证

为了验证AUTODYN软件的有效性，使用球对称计算模块分析0.229 kg的TNT在178.6 m水深处爆炸时的相关数据，并将计算值与文献［7］中的实验值进行了对比。在深水爆炸过程中，静水压力梯度可以忽略不计，取气泡周围的静水压力一定，水的计算域取为50 m。表1和表2所示为流场中距药包中心0.5 m和1 m处冲击波和气泡压力峰值、气泡最大半径及脉动周期的计算值与实验值的对比。

表1 冲击波及气泡压力峰值的AUTODYN计算值与实验值对比Tab.1 Comparison between experimental data and the calculated results about shock and bubble pressure peak

表2 气泡最大半径及脉动周期的计算值与实验值对比Tab.2 Comparison between experimental data and the calculated results about the maximum bubble radius and pulse cycle

从表中可以看出，数值模拟的气泡脉动最大半径和脉动周期与实验值间的误差约为10%，误差在可接受的范围内。一维球对称求解器可有效模拟水下爆炸气泡的脉动。

下面将以重力场为例进行对比分析。在模拟重力场中水下爆炸气泡的运动时，采用AUTODYN软件中独有的映射技术，将一维球对称计算结果映射至二维轴对称求解器，从而解决了网格尺寸过小、计算时间过长的问题。设置长、宽、高分别为18 m，18 m和7 m的流场，以35 g药量在流场中心下3.5 m处引爆，观察该工况下气泡在重力场中的运动规律。在流场中预设A，B，C等3个测点以便测量气泡在运动过程中的流场压力，它们分别位于爆心水平方向0.7 m处；爆心下方水平方向0.7 m、垂向0.71 m处；爆心上方水平方向0.7 m、垂向1.095 m处。图1所示为测定A，B，C的数值模拟压力时历曲线，通过与文献［7］中的相似工况及测点压力曲线进行对比，发现曲线的时间发展趋势以及压力峰值基本吻合，进一步验证了AUTODYN在模拟气泡在重力场中运动的精确性。

图1 流场中测点的压力数值模拟曲线Fig.1 Bubble pressure of flow field numerical simulation curve of three measuring points

图1给出了80～83 ms时间段的压力值，此时，气泡收缩至最小体积，流场中辐射气泡二次压力波。由图可见，气泡二次压力波的峰值和持续时间模拟值与实验值吻合良好，平均误差在10%以内，但压力峰值发生的时间略有提前。究其原因，可能是数值上虽然采用了边界处理，但有限的边界仍会对气泡脉动产生影响，使气泡的周期较实验值略小。

2 近固壁面的水下爆炸气泡射流

2.1 概 述

气泡与壁面的相互作用一直是研究人员关注的问题。处于固壁面附近的气泡在受到壁面Bjerknes力［8］的同时还受重力的作用，为此，设置不同的无量纲距离参数［9］（爆心距壁面的距离与气泡最大半径的比值），将炸药置于刚固平板的下方，如图2所示。

图2 气泡在固壁面下的计算模型Fig.2 Model of bubble under a solid boundary

水下爆炸气泡的模拟涉及到的流体材料包括空气、炸药和水，下面将分别介绍这3种材料的状态方程。

2.1.1 空气的状态方程

计算中可能会涉及到自由面，因此，需要对空气进行模拟。对空气进行模拟采用理想气体状态方程（Ideal Gas）：

式中，ρ=1.225×10-3g/cm3；γ为绝热指数，γ=1.4；e为比内能，e=2.06785e5kJ/kg。

2.1.2 炸药的状态方程

此处，选取TNT炸药材料模型作为水下爆炸气泡形成的条件。采用JWL状态方程描述炸药的爆轰过程：

需要说明的是，当炸药膨胀到相当的体积时，JWL方程右端的前两项可以忽略，此时，可以用理想气体的状态方程模拟炸药的行为：p=ρ(γ-1)e，其中γ=ω+1。

2.1.3 水的状态方程

AUTODYN的材料库中，水的状态方程有2种，多项式（Polynomial）状态方程和冲击（shock）状态方程。由于需要考虑静水压力，因此，本文选用多项式状态方程进行计算：

2.2 不同距离固壁面下的气泡形状

为考察固壁面下不同距离的气泡形状，计算50 kg TNT在75 m水深处爆炸时的相关数据，使用STEEL材料模拟固壁，计算中采用映射技术。设置无量纲距离γf=0.6，0.8，1.0，1.5，2.0，2.5，3.0，5.0。图3～图6所示为距离参数γf=0.6，1.0，2.0，3.0时气泡形状的演变。

图3所示为γf=0.6时的气泡形状图。在膨胀阶段，气泡受到壁面的阻碍，在103.5 ms时气泡体积达到最大，其上端成扁平状；在收缩阶段，气泡上部仍被吸附在固壁表面，在183.5 ms时下端产生指向壁面的射流，196.5 ms时射流冲破气泡，213 ms时体积达到最小。

图4所示为γf=1.0时的气泡形状图。值得注意的是，在收缩阶段，气泡上部距固壁表面的距离变近，这主要是受浮力的影响。

图5～图6所示为γf=2.0，3.0时的气泡形状图。可以看出，距离壁面较远时，壁面的影响减弱，γf=2.0与γf=3.0时的气泡形状并没有明显的差别。因γf=5.0时的气泡形状变化与自由场中的相似，因而在此不再给出。

图3 γf=0.6时的气泡形状演变Fig.3 Evolution of bubble shape with the distance parameterγf=0.6

图5 γf=2.0时的气泡形状演变Fig.5 Evolution of bubble shape with the distance parameterγf=2.0

图6 γf=3.0时的气泡形状演变Fig.6 Evolution of bubble shape with the distance parameterγf=3.0

将自由场气泡与壁面下方的气泡进行对比可见，在距离参数较大（γf=2.0，3.0，5.0）的工况下，气泡形状的演变过程与自由场相似；距离壁面越近，射流的宽度越大。此外，距离参数不同，气泡的脉动周期及射流发生的时间也不同，下面，将对脉动周期及射流时间进行分析。

2.3 固壁面对气泡最大半径、脉动周期及射流时间的影响

药包距壁面的距离不同，受到壁面的影响也不同，最直观的体现可参见气泡第1次脉动周期及最大半径，此外，射流冲击时间也有所不同。第1次脉动周期指的是气泡第1次收缩至最小体积，将要进行第2次膨胀的时刻，射流冲击时刻指的则是射流将气泡的上下表面击穿的时刻。

50 kg TNT在75 m水深处爆炸，按照经验公式计算，可得气泡第1次脉动周期为182.6 ms，第1次膨胀的最大半径为2.85 m。下面，将给出不同距离参数壁面下气泡的第1次膨胀最大半径、第1次脉动周期及射流冲击时刻的数值模拟值，如表3所示。

表3 不同距离参数壁面下的气泡数值模拟值Tab.3 Numerical simulation results with different distance parameters

近距离（γf=0.6，0.8）的气泡由于其形状不再保持为球形，因而其最大半径为最大体积的等效半径。由表3可见，距离壁面越近，气泡的最大半径越小，距离壁面越远，其最大半径就和自由场的相差越小。在γf=3.0的情况下，气泡最大半径已和自由场的相差无几。总之，壁面的存在对气泡的最大半径影响不大，在γf=0.6时仅与自由场相差（2.842-2.821）/2.842=0.7%。

气泡第1次脉动周期受壁面的影响较大，距离参数越小，脉动周期越大。在γf=3.0时，气泡脉动周期已与自由场的值相差不大，这是因为壁面在气泡膨胀过程中对其造成了阻碍作用，而在收缩时，因气泡依附在壁面上，因此距壁面越近，气泡脉动周期便越大。

为研究壁面的存在对射流时间的影响，使用气泡脉动周期作为参考量将射流冲击时间进行无量纲化：Tj′=Tj/T，绘制其随无量纲距离的变化而发生的变化，如图7所示。图中，圆圈为数值解；实线为自由场中的无量纲射流，时间Tj′=1.083；虚线代表射流冲击时刻与气泡脉动周期相同，Tj′=1 。

图7 无量纲射流冲击时间随距离参数的变化曲线Fig.7 Variation of dimensionless jet shock time with respect to the distance parameter

由图7可见，射流冲击时刻与气泡第1次脉动周期并非总是一致，但两者相差不大，即射流冲击是在气泡第1次体积达到最小时刻附近发生的，此时也是气泡开始辐射2次压力波的时刻；对于壁面下的气泡，其无量纲射流冲击时间是随着距离参数的增加而增加，这是因为壁面效应是随距壁面的距离增加而减弱，当距离参数约为2.0（γf≈2.0）时，射流冲击时刻与气泡第1次脉动周期基本相等；当距离参数大到一定（γf≥3.0）时，无量纲射流冲击时间与自由场的就已非常接近。

2.4 近固壁面气泡射流速度及压力

首先，研究气泡射流速度随距离参数的变化规律。射流速度在射流形成过程中不断提高。当射流冲击后，射流顶端的速度会有一定程度的损失，而气泡射流冲击后的射流速度是人们最为关心的。因此，定义射流冲击完成后的射流顶端速度为Vj。图8给出了γf=0.6，1.0，2.0，3.0，5.0时的射流顶端速度时历曲线及重力场中的射流顶端速度的时历曲线。

图8 不同距离参数下的射流顶端速度时历曲线Fig.8 Time history of the jet top velocity with different distance parameters

由图8可见，在射流冲击时刻，射流顶端速度急速增加，在达到最大值后开始下降。除γf=0.6的情况外，其余工况的射流顶端速度均趋于20 m/s。同时还可看出，距离参数越小，射流冲击持续的时间越短，这是因为壁面的存在使得气泡的射流在更短的时间内完成。为研究距离参数对射流顶端速度的影响，绘制了射流顶端速度最大值随距离参数的变化曲线，如图9所示。

图9 射流顶端速度峰值随距离参数的变化曲线Fig.9 Variation of max top velocity with respect to the distance parameter

由图9可见，射流顶端速度峰值随距离参数的变化规律较复杂。在距离参数较小时，射流顶端速度峰值随距离参数的增加而增加，在γf=2.0附近达到最大值，随后又随距离参数的增加而减小。当距离参数较大时，射流顶端速度峰值与自由场的相差不大。据分析，距离参数不同，壁面对气泡的作用也不同，气泡射流冲击发生的时间以及射流的形状也不同。因此，当距离参数很小时，虽然壁面的作用非常强烈，但射流冲击发生时间相对较早，而且气泡紧贴壁面，射流无法完全发展，形成的射流顶端速度峰值并不是最大的。随着距离参数的略微增大，射流时间稍微推迟，因气泡离壁面有一定的距离，此时射流能更进一步发展，因此，射流顶端速度峰值不断增加。当γf=2.0时，射流时间与气泡第1次脉动时间基本一致，此时射流发展最为完全，射流顶端速度峰值达到最大值。随着距离参数的进一步增大，壁面的影响削弱，在气泡第1个脉动周期内不再发生射流，因而射流顶端速度峰值也就不断减小。

气泡射流冲击后，射流顶端速度经过气泡顶端与壁面间流体的衰减作用，即为作用在壁面上的射流速度Vjp。图10所示为壁面附近0.1 m处射流速度峰值Vjp(max)随距离参数的变化曲线。

图10 壁面上射流速度峰值随距离参数的变化曲线Fig.10 Variation of the peak velocity of jet near a solid boundary with respect to the distance parameter

由图10可见，经过流体的衰减作用，壁面射流速度峰值与气泡顶端射流速度峰值随距离参数的变化规律又有所不同。在小距离参数（γf=0.6，0.8）情况下，无流体的衰减作用，因此射流顶端速度即为壁面射流速度。但距离参数越大，流体的衰减作用便也越大，最后作用在壁面的射流速度就越小。当γf＞3时，便可以认为作用在壁面的已不再是射流，而是气泡第2次膨胀引起的流体运动。可见，壁面上的气泡射流载荷是非常复杂的，气泡攻击壁面的威力也与多种因素有关。

图11所示为γf=0.6时壁面附近0.1 m处的压力时历曲线。

图11 γf=0.6时壁面附近0.1 m处的压力时历曲线Fig.11 Time history of pressure near a solid boundary about 0.1 m with the distance parameterγf=0.6

由图11可见，当水射流跃击到壁板上时，速度降至零，同时在水中产生压缩波，作用在平板上的压力即为水锤压力。水锤压力为：

式中，Ph为水锤压力；ρ为水的密度；C为水的声速；Vj为作用在平板上的射流速度。按式（4）计算，γf=0.6时壁面的水锤压力应为145 MPa，数值模拟值为170 MPa，考虑到此测点并非壁面上的点，数值模拟值是在误差允许的范围内。同时注意到，在射流压力之后，程序可以捕捉到气泡的二次压力波，其峰值为20.56 MPa。

气泡在发生射流之后，仍然会第2次膨胀，图12所示为γf=0.6，1.5，3.0时爆心水平方向2.6 m处的气泡脉动压力时历曲线，图13所示为爆心水平方向2.6 m处气泡脉动压力峰值随距离参数的变化曲线。

图12 爆心水平方向2.6 m处气泡脉动压力时历曲线Fig.12 Time history of bubble pulse pressure at 2.6 m from explosion center horizontally

图13 爆心水平方向2.6 m处气泡脉动压力峰值随距离参数变化曲线Fig.13 Variation of bubble pulse pressure peak with respect to distance parameter at 2.6 m from explosion center horizontally

由图12和图13可见，水平方向的气泡脉动压力峰值随距离参数的变化规律与壁面附近射流速度的规律是相反的，在γf=0.8附近最小。这是因为气泡射流消耗的能量越多，气泡再次膨胀的能量就会越发缩小，相应的脉动压力也就越小。这也再次验证了在γf=0.8时气泡对壁面的射流作用最大。

3 结 论

本文运用AUTODYN软件分别模拟了球对称气泡、重力场中气泡以及刚性和弹性壁面与气泡之间的相互作用，计算结果与实验数据的对比和分析表明：

1）运用AUTODYN软件计算得出的一维及重力场中水下爆炸气泡的平均误差在10%以内，完全适用于气泡动力学的研究。

2）刚性壁面的存在对气泡的运动特性具有较大影响。对于固壁面下方的水下爆炸气泡，距离参数γf越小，壁面的影响就越大，气泡射流的宽度越大，第1次脉动周期越大。

3）气泡的无量纲射流冲击时间随距离参数γf的增加而增加，在γf≈2.0时，无量纲射流冲击时间为1，气泡射流发展最为完全，射流顶端速度峰值达到最大值。距离参数γf越小，气泡射流冲击持续的时间就越短。

4）壁面上的气泡射流载荷非常复杂，与气泡顶端射流速度、气泡射流形状以及射流顶端距离壁面的距离均有关。

［1］张阿漫，王诗平，白兆宏，等.不同环境下气泡脉动特性实验研究［J］.力学学报，2011，43（1）：71-83.ZHANG A M，WANG S P，BAI Z H，et al.Experimental study on bubble pulse features under different circumstances［J］.Chinese Journal of Theoretical and Applied Mechanics，2011，43（1）：71-83.

［2］张阿漫，姚熊亮.近自由面水下爆炸气泡的运动规律研究［J］.物理学报，2008，57（1）：339-353.ZHANG A M，YAO X L.The law of the underwater explosion bubble motion near free surface［J］.Acta Physica Sinica，2008，57（1）：339-353.

［3］VUYST T D，VIGNJEVIC R，CAMPBELL J C.Coupling between meshless and finite element methods［J］.International Journal of Impact Engineering，2005，31（8）：1054-1064.

［4］COLE R H.Underwater explosion［M］.Princeton：Princeton University Press，1948.

［5］GEERS T L，HUNTER K S.An integrated waveeffects model for an underwater explosion bubble［J］.The Journal of the Acoustical Society of America，2002，111（4）：1584-1601.

［6］BRUJAN E A，NAHEN K，SCHMIDT P，et al.Dynamics of laser-induced cavitation bubbles near an elastic boundary：influence of the elastic modulus［J］.Journal of Fluid Mechanics，2001，433：251-281.

［7］KLASEBOER E，HUNG K C，WANG C，et al.Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/rigid structure［J］.Journal of Fluid Mechanics，2005，537：387-413.

［8］BJERKNES.Fields of force［M］.New York：Columbia University Press，1966.

［9］张阿漫.水下爆炸气泡三维动态特性研究［D］.哈尔滨：哈尔滨工程大学，2006.