水下爆炸异相气泡动力学特性的Euler有限元数值模拟研究

唐 皓, 刘云龙, 冯集团, 鞠欣洋, 张阿漫,2

(1. 哈尔滨工程大学 船舶工程学院, 哈尔滨 150001;2. 哈尔滨工程大学南海研究院, 海南 三亚 572024)

0 引 言

气泡是自然界中常见的物理现象,气泡动力学在人类的生产生活中有着广泛的应用[1-4]．超声波气泡可以用来清洗[5],气枪产生的高压气泡可以用来勘探海洋资源[6]．除此之外,螺旋桨产生的空泡会对结构产生剥蚀[7],水下爆炸气泡会对舰艇造成严重毁伤[8-9]．水下爆炸不仅可以通过冲击波毁伤舰艇,还可以通过气泡脉动和气泡射流对舰艇造成进一步的毁伤[10-11]．关于气泡动力学的研究不仅可以提高生产生活的效率[12],而且对国防安全和海军装备具有重要的意义．

气泡往往不是以单一气泡的形式出现,而是以多气泡的形式出现．在海战中,舰艇受到水下和水面武器攻击时,通常会面临多发武器的威胁,由于鱼雷和水雷等武器触发时间不同,所产生的气泡往往是不同相位的．单气泡的动力学过程本身就包含射流和气泡撕裂等复杂拓扑结构变化,异相多气泡之间耦合过程更为复杂,气泡之间的穿透、撕裂和融合存在很强的几何非线性,异相双气泡问题是多气泡问题中较为典型的一种．关于多气泡的耦合过程很难用单一理论来进行探究,研究人员主要通过实验和数值方法对多气泡的耦合进行探究．利用激光气泡实验,Tomita等[13]发现了同相双气泡的运动和气泡间距的关系,并且研究了刚性壁面附近的双气泡运动特性．Cho等[14]通过实验发现在一定参数下,气泡附近的小气泡可以使气泡射流的冲量达到原来的两倍．Tomita等[15]研究了在不同时间产生的两个激光诱导气泡的射流行为,并发现通过控制控制参数可以实现强化射流冲击．Chew等[16]利用电火花气泡探究了异相气泡初始距离和相位差的影响,并将异相双气泡的耦合类型分为融合、弹射、对向射流和反向射流四种．Liang等[17]进一步利用电火花实验对异相双气泡进行探究,增加了未弹射的类别．Cui等[18]利用电火花气泡探究了双气泡破冰问题,开拓了双气泡应用的新方向．

由于实验研究难以探究双气泡耦合的详细过程,气泡撕裂融合的细节以及气泡内部的情况很难通过摄像机记录到,气泡的能量等物理量也很难检测．因此许多研究学者利用数值计算方法去完善双气泡的研究,边界积分方法(BIM)是最为典型的计算方法[19-20]．边界积分方法具有精度高、计算速度快的优势,但是其在计算气泡撕裂和融合等复杂拓扑结构变化时需要人工干预．采用VOF方法[21-22]和LS方法[23-25]的多相流求解器是求解气泡动力学的另外一种途径, 在处理界面大变形问题时具有很大的优势．LS方法和VOF方法都是Euler界面捕捉方法,通过对示踪函数进行处理从而确定界面的位置和形状．除了Euler界面捕捉方法之外还有Lagrange追踪方法(界面追踪方法),例如front tracking方法[26-28]等．Barras等[29]采用任意Euler-Lagrange方法(ALE)针对水下爆炸在无限介质中产生的气泡动力学问题进行了建模和研究．Li等[30]利用OpenFoam结合VOF方法对同相双气泡耦合问题进行了模拟计算,对双气泡的耦合方式进行了模式分类．He和Liu等[31-35]提出了Euler有限元方法(EFEM),将Euler方程拆分为Lagrange计算步和Euler计算步进行求解,在模拟强间断冲击波的同时,也能很好地处理水下爆炸气泡动力学问题中多项流界面变化问题．在此基础上,Liu等[36]基于EFEM结合VOF方法实现了对多相界面的追踪,建立了同相水下爆炸双气泡耦合的轴对称力学模型,模拟了强间断冲击波传播与双气泡相互作用的连续变化过程,探究了距离参数和浮力参数对同相双气泡的影响．

本文基于EFEM建立了水下爆炸异相气泡动力学模型,实现了异相水下爆炸的全过程数值模拟．首先本文介绍了EFEM的基本理论,然后将数值计算结果和水下爆炸实验以及气泡统一理论进行对比,验证了数值计算模型的正确性．在此基础上本文对比了异相气泡和单气泡动力学的区别,探究了相位和距离参数对异相双气泡耦合过程以及流程压力载荷的影响．

1 Euler有限元数值计算模型

1.1 问题描述和变量声明

本文用轴对称模型计算异相爆炸问题,如图1所示．以上气泡的初始位置作为坐标原点O,并建立坐标系．选取上气泡作为参考气泡,上气泡在自由场的最大半径为Rm,第一脉动周期为T1,起爆时刻为t1．下气泡的起爆时刻为t2,第一脉动周期为T2,上气泡和下气泡的间距为dbb,无量纲距离参数定义为γbb=dbb/Rm,异相双气泡的相位差定义为Δθ=2π·(t2-t1)/Ti,如果上气泡先起爆,则i=1,反之为2．计算域的大小为4Rm×10Rm,边界为无反射边界条件[32],无反射边界条件可以有效地减弱压力波在边界的反射．

图1 异相双气泡问题示意图Fig. 1 Schematic diagram of the out-of-phase double-bubble problem

1.2 EFEM

水下爆炸气泡的尺度较大,运动过程属于高Reynolds数问题,计算域流体的黏性和表面张力可以忽略．假设水下爆炸问题涉及流体是无黏性且可压缩的,因此可以用Euler方程来描述流场的运动过程[37],体积分数、质量、动量和能量方程如下:

(1)

式中i为流体类别;f为体积分数,当值为1时表示完全被该流体占有,值为0时表示不含该流体;t为时间;u为流场速度矢量;ρ为流场密度;p为压力;g为重力加速度;ein代表单位质量的内能;K=ρc2为流体体积模量,c表示单元声速．当两种流体在一个单元混合时,平均密度为

(2)

平均模量可以表示为

(3)

混合单元内的平均声速为

(4)

对于可压缩流体而言,需要引入流体状态方程来封闭Euler方程组．本文数值计算水的状态方程用的是Tammann方程[38]:

p=ρein(γ-1)-γPw,

(5)

式中,γ取值为7.15,Pw取值为330.9 MPa,声速可以表示为

(6)

本文涉及高压气泡的计算,气泡内部气体采用理想气体状态方程[38],具体如下:

p=ρein(γg-1),

(7)

式中γg为气体的热容比,是一个无量纲数,取值为1.25．理想气体状态方程也可以写成Tammann方程[38]．针对水下爆炸气泡,爆轰气体采用JWL状态方程[39],具体如下:

(8)

本文数值计算时炸药种类为TNT炸药,A为371.2 GPa,B为3.2 GPa,R1为4.15,R2为0.95,ω为0.25,ein初始值为4.3 MJ/kg,ρ0为1 630 kg/m．

在Lagrange步,计算域网格节点跟随流体运动．本文采用显式有限元方法求解,通过分部积分和GaussGreen公式得到不包含对流项的半离散形式的动量方程:

(9)

(10)

(11)

式中M为质量矩阵,M′是M对时间的导数,U表示节点速度对应的列向量,F表示节点力．依据数值分析的CFL条件,EFEM的稳定时间步长为

(12)

式中NCFL本文取值0.5,上标“min”表示计算域内所有单元的最小值,Le为单元的最小尺寸．

在Euler步,需要将网格移动回原始位置,流体与网格之间发生相对运动．通过变形之后的网格和原始网格之间的守恒关系,计算出相邻单元之间的物理量输运．具体为利用PLIC算法重建多介质单元中的介质界面,之后应用VOF方法计算输运体积,使用MUSCL解单元中心变量输运的质量和能量．使用半指数位移法输运节点中心动量,密度以及内能更新之后利用流体的状态方程求解单元中心压力,最后根据单元声速计算时间步增量．将水下爆炸过程简化为轴对称问题,可以极大地提高计算效率．在柱坐标系的Euler方程中,散度算子和梯度算子如下:

(13)

(14)

1.3 数值计算模型验证及收敛性分析

本文首先将数值计算模型和气泡统一理论进行对比并进行收敛性分析．气泡统一理论是由Zhang(张阿漫)等[40]推导出的,其考虑气泡迁移特性的脉动方程为

(15)

式中,C为声速,t为时间,R为气泡半径．方程左侧相当于对应问题的驱动力,dF/dt可以根据具体的物理问题确定．气泡统一理论可以根据不同的物理问题进行扩展,同时保持统一、简洁和优雅的数学形式．其展开形式为

(16)

对比工况为两个高压异相气泡,相位差为π,上气泡和下气泡的能量比为E1/E2=1.6．不考虑重力,环境压力P∞为394 kPa．上气泡先开始膨胀,初始内压为100倍P∞,初始半径为0.02 m．下气泡位于上气泡下方0.567 7 m的位置, 初始内压为100倍P∞, 初始半径为0.018 m, 相比于上气泡延迟0.005 677 s开始膨胀．本文除了对比两种方法计算的气泡的等效半径变化, 也对比了流场的压力, 流场中压力测点为(0.683 0 m, -0.339 7 m)．本文利用Euler有限元数值计算模型进行了4种网格精度的计算．当下气泡在第一次收缩达到最小体积时,气泡统一理论计算的气泡半径为2.079 cm,和气泡统一理论的计算结果相比,网格精度为0.08Rm,0.04Rm,0.02Rm和0.01Rm时的EFEM计算结果相对误差分别为12.2%,1.4%,1.4%和2.7%．从图2(a)中可以看出当网格精度从0.02Rm提高到0.01Rm时,下气泡等效半径变化很小,说明用0.01Rm的网格计算满足计算精度要求．

(a) 气泡等效半径变化对比(b) 流场压力对比(a) Comparison of bubble equivalent radius changes(b) Comparison of flow field pressures图2 EFEM计算结果与气泡统一理论[40]计算结果对比Fig. 2 Comparisons between the EFEM calculation results and the unified theory calculation results

EFEM的计算结果和气泡统一理论在上气泡和下气泡的第一个周期吻合较好,从第二个计算周期相差开始明显,并且EFEM计算的气泡等效半径略小于气泡统一理论,这主要是由于EFEM计算时的能量耗散较大．在图2(b)中上气泡冲击波阶段和下气泡冲击波阶段两种计算方法的结果几乎完全一致,包括下气泡的冲击波在上气泡表面的反射波扰动细节都很一致．气泡载荷除在时间上两种方法有差别,在载荷峰值和脉宽上相差很小．

图3 EFEM计算结果与试验结果[41]的气泡形态对比Fig. 3 Comparison of bubble shapes between the EFEM calculation results and the experimental results[41]

除了与理论方法进行对比,本文还将数值计算结果和试验进行对比,从气泡形态上验证数值计算模型．试验工况选自Liu等[41]关于异相爆炸的论文,试验是在一个4 m×4 m×4 m的水箱中进行的．试验中用的炸药当量为10 g,转换成TNT需要乘以转换系数1.319,等效为13.19 g TNT．试验中上气泡爆炸水深为1.55 m,相位差Δθ=-π．当Δθ=-π,下气泡膨胀到最大体积时上气泡产生．如图3中t=30.34 ms所示,在上气泡开始膨胀后,冲击波迅速到达下气泡表面．冲击波在下气泡表面发生反射,并产生空化．随后下气泡开始产生向下的射流,上气泡开始向下凸起变化成“蛋”形．当t=51.26 ms时下气泡已经被射流击穿,上气泡向上的射流开始形成．t=70.09 ms时上气泡被射流击穿,并且上气泡的射流宽度相较于下气泡很小．在随后的运动中,两个气泡朝着相反的方向运动．

从图3中试验和数值计算的结果对比可以看出,在上气泡和下气泡的全周期过程中数值计算结果和试验结果都基本一致．本文建立的Euler有限元数值计算模型在和理论方法和试验结果的对比中都得到了验证．

2 结果与讨论

2.1 自由场单气泡和异相双气泡对比

异相双气泡之间的相互作用过程是十分复杂的,在分析相位等因素对异相爆炸气泡的影响之前有必要对自由场单气泡和异相双气泡进行对比分析．计算工况为100 kg TNT在50 m水深爆炸,气泡距离参数γbb=2Rm,相位差Δθ=π,压力测点为(-2Rm,1.5Rm)．

如图4(a)所示,在下气泡开始膨胀后上气泡等效半径开始迅速减小,减小的速度明显快于单气泡的情况．下气泡的膨胀速度也比单气泡的情况快,并且下气泡的最大体积也比单气泡时更大．结合图4(b)中气泡的迁移过程分析,在下气泡开始膨胀后上气泡被迅速推向上方,并且下气泡受上气泡的影响被推向相反的方向．Bjerknes力是超声场中声波作用下气泡受到的附加力．下气泡不仅运动方向发生了改变,气泡的射流方向也是向下的,这是因为气泡间的Bjerknes力大于下气泡的浮力作用．从图4(c)中可以看出,受到上气泡达到最大体积后收缩的影响,下气泡的形态更长并且气泡的整体轮廓都比单气泡时大．气泡的伸长率为气泡水深方向的长度和气泡宽度的比值,单气泡在自由场膨胀的过程中气泡伸长率基本维持在1.0,气泡保持球形．在异相爆炸中,延迟起爆的下气泡在膨胀过程中伸长率大于1,气泡在膨胀时是“蛋”形而不是球形．异相双气泡的流场压力曲线相较于单气泡要复杂得多,这主要来自气泡表面对冲击波的反射作用．从图4(d)中可以看出,下气泡起爆后不久冲击波到达上气泡表面发生反射,并产生空化效应．此外异相爆炸时上气泡的气泡载荷也要小于单气泡的气泡载荷,这与上气泡的能量转移有关．

(a) 等效半径(b) 气泡迁移(a) Equivalent radii(b) Bubble migration histories

(c) 气泡形态(d) 流场压力(c) Bubble shapes(d) Flow field pressures图4 异相爆炸双气泡和自由场单气泡对比Fig. 4 The comparison between out-of-phase explosion bubbles and the free-field single bubble

关于气泡能量的计算参考自Tian等[33,37]的论文,气泡的总能量可由下式计算:

(17)

如图5(a)所示,炸药在起爆后内能出现极大地下降,这部分减少的内能主要转化为冲击波的能量．自由场单气泡在起爆之后,冲击波向周围传播,不存在冲击波和气泡之间的相互作用．异相爆炸时,冲击波在气泡之间的反复作用使气泡的能量变化更为复杂．在异相爆炸气泡的相互作用过程中,气泡间的内能变化与单气泡相比有着很大的差异．对于先起爆的上气泡而言,气泡的内能在气泡坍缩时的峰值与单气泡相比要小．延迟起爆的下气泡的内能相对于单气泡而言峰值和脉宽都大幅提高．总体上看,异相爆炸中上气泡和下气泡的内能相较于单气泡都增加了．图5(b)中,在下气泡起爆后,上气泡和下气泡的势能都增加了,下气泡相较于上气泡增加的幅度更大．图5(c)所示为自由场单气泡流场动能和异相双气泡流场动能变化,异相爆炸气泡流场动能为双气泡的总流场动能．在上气泡达到最大体积时流场动能接近于零,此时下气泡起爆流场动能瞬间增大,但增加的动能大于单气泡时气泡起爆增加的动能．结合图5(d)分析气泡的总能量变化,在炸药起爆后总能量短时间内衰减的能量为冲击波的能量．只有下气泡存在时气泡系统损失的冲击波能量约为2.659×108J,异相爆炸时下气泡起爆后气泡系统损失的冲击波能量约为1.651×108J．异相爆炸时冲击波带走的能量明显小于单气泡时冲击波带走的能量,这说明下气泡起爆后冲击波对异相气泡系统做功,使异相气泡系统的动能和总能量增加,从而导致相应的气泡势能和内能增加．

(a) 内能(b) 势能(a) Internal energy(b) Potential energy

(c) 动能(d) 总能量(c) Kinetic energy (d) Total energy图5 异相爆炸双气泡和自由场单气泡的能量变化Fig. 5 Energy changes of out-of-phase explosion bubbles and the free-field single bubble

2.2 相位对异相双气泡耦合的影响

相位是影响异相双气泡相互作用的关键因素,本文针对垂向异相双气泡问题将相位差分为2.0π,1.5π,1.0π,0.5π,0,-0.5π,-1.0π,-1.5π以及-2.0π共9种情况来讨论．计算工况为100 kg TNT在50 m水深爆炸,气泡距离参数γbb=2Rm．

上气泡和下气泡在不同相位下的迁移曲线如图6所示,图例标签后面的箭头表示的是气泡在第一周期的射流方向．通过对比可以发现,上下气泡的迁移和相位差并不是线性关系．当Δθ的值为正时表示上气泡先起爆,为负时表示下气泡先起爆．通过之前的分析已经知道了冲击波对气泡做功是增加异相爆炸气泡系统能量的主要原因,当Δθ=π时上气泡体积最大,下气泡此时起爆产生的冲击波对上气泡的作用面积最大,因此图6(a)所示Δθ=π时上气泡运动得最远．当Δθ=-π时下气泡体积最大,因此图6(b)所示,Δθ=-π时下气泡运动得最远．对于射流方向,由于当|Δθ|=2π时先起爆的气泡已经产生很大位移,两个气泡之间的相互作用很弱,因此|Δθ|=2π不作为射流方向的参考．除|Δθ|=2π外,只有Δθ=0时上下两个气泡产生对向射流,其余参数下均产生相反方向的射流．这是因为Δθ=0时上下两个气泡近似镜像,气泡的运动和刚性壁面附近的气泡运动相似,产生朝向壁面的射流．后起爆的气泡在已经存在的气泡附近运动和自由面附近的气泡运动类似,产生反向射流．

(a) 上气泡(b) 下气泡(a) The upper bubble(b) The lower bubble图6 不同相位异相爆炸双气泡的迁移Fig. 6 Migration histories of out-of-phase explosion bubbles in different phases

在异相爆炸双气泡的耦合过程中,异相气泡的周期相对自由场单气泡会发生改变．如图7所示,如果下气泡先起爆的话,在上气泡起爆后受上气泡干扰下气泡会提前坍缩,并且随着相位差的增大下气泡坍缩得越来越早,周期也越来越短．反过来,上气泡受到下气泡收缩时的吸引力而被拉长,上气泡的周期也越来越长．如果上气泡先起爆, 在下气泡起爆后上气泡受到下气泡的作用会提前坍缩, 并且相位差越小上气泡的周期越小．反过来,下气泡受到上气泡收缩时的吸引力而被拉长,下气泡的周期随着相位差减小也越来越长．相位差为0是一种特殊的情况,上下气泡的周期均大于自由场单气泡的情况．以上规律可以总结为:先起爆的气泡周期受后起爆气泡影响会缩短,相位差绝对值越小周期越小．后起爆的气泡由于受到先起爆气泡收缩的影响周期变长,相位差绝对值越小周期越大．

图7 不同相位异相爆炸双气泡的周期图8 不同相位异相爆炸双气泡系统的总能量变化Fig. 7 Periods of out-of-phase explosion bubbles Fig. 8 Total energy changes of out-of-phase explosion in different phases bubbles system with different phases

相位对异相双气泡耦合过程的影响也包括对双气泡系统总能量的影响,如图8所示．在自由场水下爆炸中,爆炸冲击波向周围流场传播会带走相当大一部分能量．Δθ=0时,上下气泡同时起爆,在初始时刻上下气泡体积都很小,爆炸冲击波的作用面积有限,因此大部分冲击波穿过气泡散失到周围流场．Δθ=0时,初始时刻上下两气泡总能量为8.7×108J,冲击波过后系统总能量降至3.0×108J,损失的冲击波能量为5.7×108J,占比65.5%．当相位差不为零时,异相气泡系统总能量在先起爆的第一次冲击波过后从4.3×108J降至1.7×108J,损失的冲击波能量为2.6×108J,占比60.5%,参考同相爆炸理论上第二次爆炸将会损失3.1×108J．但实际上当|Δθ|=0.5π时,后起爆的爆炸冲击波带走的能量约为1.8×108J;当|Δθ|=π时,后起爆的爆炸冲击波带走的能量约为1.7×108J;当|Δθ|=1.5π时,后起爆的爆炸冲击波带走的能量约为1.6×108J．这部分冲击波少带走的能量通过冲击波和反射波对气泡做功留在异相气泡系统总能量内．并且后起爆的气泡越接近|Δθ|=π,气泡总能量损失得越少．

2.3 距离参数对异相双气泡耦合的影响

在异相爆炸气泡相互作用的过程中,除了相位的影响,距离参数也是影响耦合过程的重要因素．本小节计算工况为100 kg TNT在50 m水深爆炸,气泡相位差Δθ=-π,下气泡在上气泡体积最大时起爆．4种距离参数下的异相爆炸气泡动态特性如图9所示,每种距离参数下气泡的动态过程按照时间顺序从左到右排列．当γbb=1.25时,下气泡在起爆后很快就被坍缩的上气泡吸入,并且下气泡顶部在上气泡内部破裂,上下气泡内部连通,破裂之后的液滴和液膜撞击到上气泡顶部将上气泡撕裂．

(a) γbb=1.25

(b) γbb =1.50

(c) γbb =1.75

(d) γbb =2.00图9 不同距离参数异相爆炸双气泡动态特性Fig. 9 Dynamic characteristics of out-of-phase explosion bubbles with different distance parameters

图10 不同距离参数异相爆炸双气泡系统的总能量变化Fig. 10 Total energy changes of out-of-phase explosion bubbles with different distance parameters

图9(a)所示异相气泡耦合过程为穿透融合模式,与此作为区别的是图9(b)．图9(b)所示为γbb=1.50时异相双气泡的动态过程,当t=0.213 s时,上气泡刚要被射流击穿,此时两个气泡并未相连通,在随后的坍缩过程下气泡顶部突起的地方被环向射流分成上下两部分,下面的主体部分产生很细的向下的射流．图9(b)这种一个气泡进入另一个气泡内部但并未出现融合的情况定义为进入未融合模式．当距离参数继续增大时,如图9(c)和图9(d)所示,双气泡在脉动的过程中始终保持一定间隔,本文定义为未进入模式．除穿透融合模式外,距离参数越小后起爆的气泡被拉伸得越长,产生的射流也越细．

图10所示为不同距离参数时,异相爆炸气泡系统总能量的变化过程,在下气泡起爆之前4种距离参数下的气泡能量变化是完全相同的．当下气泡起爆后,在冲击波的传播过程中,异相双气泡系统总能量变化开始不同．距离参数越小,冲击波带走气泡系统的能量越少．γbb=1.25时,由于上气泡被穿透撕裂,撕裂之后会有很多小的气泡向流场辐射压力波,因此气泡能量随后相较于其他参数会出现较大下降．

3 结 论

本文基于EFEM和VOF方法建立了水下爆炸异相气泡动力学模型,对垂向异相双气泡非线性耦合问题进行了研究．本文通过将计算结果和气泡统一理论以及异相爆炸试验进行对比,验证了数值计算模型的有效性．异相爆炸气泡的脉动、坍缩和穿透融合等特性都得到了很好地模拟,并通过对比分析得到以下结论:

1) 与自由场单气泡相比,异相爆炸中先爆炸的气泡会提前坍缩周期会缩短．冲击波在气泡表面会发生反射并产生空化,受异相爆炸气泡干扰气泡的射流方向和自由场也不同．气泡起爆后冲击波对异相气泡系统做功,使异相气泡系统的动能和总能量增加,从而导致相应的气泡势能和内能增加．

2) 当相位差不为0时,异相双气泡会产生反向射流,先产生的气泡周期会缩短,相位差绝对值越小先产生的气泡周期越短,后产生的气泡周期越长．当相位差为0时,异相双气泡产生对向射流,上下气泡的周期都会变长．相位差越接近 ,冲击波带走的能量越少,气泡总能量损失得越少．

3) 根据不同距离参数下的异相爆炸气泡耦合现象,本文将异相爆炸气泡分为3种模式,即穿透融合模式、进入未融合模式和未进入模式．除穿透融合模式外,距离参数越小后起爆的气泡被拉伸得越长,产生的射流也越细．距离参数越小,冲击波带走气泡系统的能量越少．