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

应用数学和力学编委会，：水下爆炸异相气泡动力学特性的 有限元数值模拟研究唐 皓，刘云龙，冯集团，鞠欣洋，张阿漫，（哈尔滨工程大学 船舶工程学院，哈尔滨；哈尔滨工程大学南海研究院，海南 三亚）（我刊编委张阿漫来稿）摘要：该研究基于 有限元方法建立了水下爆炸异相气泡动力学模型，将计算结果和气泡统一理论以及异相爆炸试验进行对比验证计算模型的有效性 通过和自由场单气泡进行对比分析，发现异相爆炸冲击波对气泡做功是异相爆炸气泡总能量增加的原因，相位差绝对值越接近、距离参数越小气泡总能量损失越少，后产生的气泡会使先产生的气泡提前坍缩 气泡的射流方向受相位差影响，相位差为零时产生对向射流，其他相位差产生反向射流 关 键 词：水下爆炸；有限元；异相气泡中图分类号：文献标志码：，（，；，）（，）：，应用数学和力学 卷 期 年 月 ，收稿日期：；修订日期：基金项目：海南省财政科技计划（）；国防基础科研项目（）作者简介：唐皓（），男，博士（：）；刘云龙（），男，教授，博士，博士生导师（通讯作者：）引用格式：唐皓，刘云龙，冯集团，鞠欣洋，张阿漫 水下爆炸异相气泡动力学特性的 有限元数值模拟研究 应用数学和力学，（）：；引 言气泡是自然界中常见的物理现象，气泡动力学在人类的生产生活中有着广泛的应用 超声波气泡可以用来清洗，气枪产生的高压气泡可以用来勘探海洋资源 除此之外，螺旋桨产生的空泡会对结构产生剥蚀，水下爆炸气泡会对舰艇造成严重毁伤 水下爆炸不仅可以通过冲击波毁伤舰艇，还可以通过气泡脉动和气泡射流对舰艇造成进一步的毁伤 关于气泡动力学的研究不仅可以提高生产生活的效率，而且对国防安全和海军装备具有重要的意义 气泡往往不是以单一气泡的形式出现，而是以多气泡的形式出现 在海战中，舰艇受到水下和水面武器攻击时，通常会面临多发武器的威胁，由于鱼雷和水雷等武器触发时间不同，所产生的气泡往往是不同相位的 单气泡的动力学过程本身就包含射流和气泡撕裂等复杂拓扑结构变化，异相多气泡之间耦合过程更为复杂，气泡之间的穿透、撕裂和融合存在很强的几何非线性，异相双气泡问题是多气泡问题中较为典型的一种 关于多气泡的耦合过程很难用单一理论来进行探究，研究人员主要通过实验和数值方法对多气泡的耦合进行探究 利用激光气泡实验，等发现了同相双气泡的运动和气泡间距的关系，并且研究了刚性壁面附近的双气泡运动特性 等通过实验发现在一定参数下，气泡附近的小气泡可以使气泡射流的冲量达到原来的两倍 等研究了在不同时间产生的两个激光诱导气泡的射流行为，并发现通过控制控制参数可以实现强化射流冲击 等利用电火花气泡探究了异相气泡初始距离和相位差的影响，并将异相双气泡的耦合类型分为融合、弹射、对向射流和反向射流四种 等进一步利用电火花实验对异相双气泡进行探究，增加了未弹射的类别 等利用电火花气泡探究了双气泡破冰问题，开拓了双气泡应用的新方向 由于实验研究难以探究双气泡耦合的详细过程，气泡撕裂融合的细节以及气泡内部的情况很难通过摄像机记录到，气泡的能量等物理量也很难检测 因此许多研究学者利用数值计算方法去完善双气泡的研究，边界积分方法（）是最为典型的计算方法 边界积分方法具有精度高、计算速度快的优势，但是其在计算气泡撕裂和融合等复杂拓扑结构变化时需要人工干预 采用 方法和 方法的多相流求解器是求解气泡动力学的另外一种途径，在处理界面大变形问题时具有很大的优势 方法和 方法都是 界面捕捉方法，通过对示踪函数进行处理从而确定界面的位置和形状 除了 界面捕捉方法之外还有 追踪方法（界面追踪方法），例如 方法等 等采用任意 方法（）针对水下爆炸在无限介质中产生的气泡动力学问题进行了建模和研究 等利用 结合 方法对同相双气泡耦合问题进行了模拟计算，对双气泡的耦合方式进行了模式分类 和 等提出了 有限元方法（），将 方程拆分为 计算步和 计算步进行求解，在模拟强间断冲击波的同时，也能很好地处理水下爆炸气泡动力学问题中多项流界面变化问题 在此基础上，等基于 结合 方法实现了对多相界面的追踪，建立了同相水下爆炸双气泡耦合的轴对称力学模型，模拟了强间断冲击波传播与双气泡相互作用的连续变化过程，探究了距离参数和浮力参数对同相双气泡的影响 本文基于 建立了水下爆炸异相气泡动力学模型，实现了异相水下爆炸的全过程数值模拟 首先本文介绍了 的基本理论，然后将数值计算结果和水下爆炸实验以及气泡统一理论进行对比，验证了数值计算模型的正确性 在此基础上本文对比了异相气泡和单气泡动力学的区别，探究了相位和距离参数对异相双气泡耦合过程以及流程压力载荷的影响 有限元数值计算模型 问题描述和变量声明本文用轴对称模型计算异相爆炸问题，如图 所示 以上气泡的初始位置作为坐标原点，并建立坐标应 用 数 学 和 力 学 年 第 卷系 选取上气泡作为参考气泡，上气泡在自由场的最大半径为，第一脉动周期为，起爆时刻为 下气泡的起爆时刻为，第一脉动周期为，上气泡和下气泡的间距为，无量纲距离参数定义为 ，异相双气泡的相位差定义为 （），如果上气泡先起爆，则 ，反之为 计算域的大小为 ，边界为无反射边界条件，无反射边界条件可以有效地减弱压力波在边界的反射 图 异相双气泡问题示意图 水下爆炸气泡的尺度较大，运动过程属于高 数问题，计算域流体的黏性和表面张力可以忽略 假设水下爆炸问题涉及流体是无黏性且可压缩的，因此可以用 方程来描述流场的运动过程，体积分数、质量、动量和能量方程如下：（），（）（），（）（），（）（），（）式中 为流体类别；为体积分数，当值为 时表示完全被该流体占有，值为 时表示不含该流体；为时间；为流场速度矢量；为流场密度；为压力；为重力加速度；代表单位质量的内能；为流体体积模量，表示单元声速 当两种流体在一个单元混合时，平均密度为 ，（）平均模量可以表示为 ，（）混合单元内的平均声速为 （）第 期 唐皓，等：水下爆炸异相气泡动力学特性的 有限元数值模拟研究对于可压缩流体而言，需要引入流体状态方程来封闭 方程组 本文数值计算水的状态方程用的是 方程：（），（）式中，取值为，取值为 ，声速可以表示为 （）本文涉及高压气泡的计算，气泡内部气体采用理想气体状态方程，具体如下：（），（）式中 为气体的热容比，是一个无量纲数，取值为 理想气体状态方程也可以写成 方程 针对水下爆炸气泡，爆轰气体采用 状态方程，具体如下：（）本文数值计算时炸药种类为 炸药，为，为，为，为，为，初始值为 ，为 在 步，计算域网格节点跟随流体运动 本文采用显式有限元方法求解，通过分部积分和 公式得到不包含对流项的半离散形式的动量方程：（），（）式中 为计算域，为计算域的边界，是计算域边界的法向量，和 为节点 和 的形函数 节点在 方向的加速度向量 可由以下矩阵计算得到：，（），（）式中 为质量矩阵，是 对时间的导数，表示节点速度对应的列向量，表示节点力 依据数值分析的 条件，的稳定时间步长为 ，（）式中 本文取值，上标“”表示计算域内所有单元的最小值，为单元的最小尺寸 在 步，需要将网格移动回原始位置，流体与网格之间发生相对运动 通过变形之后的网格和原始网格之间的守恒关系，计算出相邻单元之间的物理量输运 具体为利用 算法重建多介质单元中的介质界面，之后应用 方法计算输运体积，使用 解单元中心变量输运的质量和能量 使用半指数位移法输运节点中心动量，密度以及内能更新之后利用流体的状态方程求解单元中心压力，最后根据单元声速计算时间步增量 将水下爆炸过程简化为轴对称问题，可以极大地提高计算效率 在柱坐标系的 方程中，散度算子和梯度算子如下：，（）（）数值计算模型验证及收敛性分析本文首先将数值计算模型和气泡统一理论进行对比并进行收敛性分析 气泡统一理论是由（张阿漫）等推导出的，其考虑气泡迁移特性的脉动方程为 ，（）式中，为声速，为时间，为气泡半径 方程左侧相当于对应问题的驱动力，可以根据具体的物理问应 用 数 学 和 力 学 年 第 卷题确定 气泡统一理论可以根据不同的物理问题进行扩展，同时保持统一、简洁和优雅的数学形式 其展开形式为 ，（）式中，是气泡的迁移速度，为气泡表面的焓差，和 分别表示气泡脉动、迁移和环境流场引起的等效力 方程（）是气泡统一理论的核心，可以扩展到各种场景，如多气泡的相互作用研究和不同边界的气泡动力学研究 对比工况为两个高压异相气泡，相位差为，上气泡和下气泡的能量比为 不考虑重力，环境压力为 上气泡先开始膨胀，初始内压为 倍，初始半径为 下气泡位于上气泡下方 的位置，初始内压为 倍，初始半径为 ，相比于上气泡延迟 开始膨胀 本文除了对比两种方法计算的气泡的等效半径变化，也对比了流场的压力，流场中压力测点为（，）本文利用 有限元数值计算模型进行了 种网格精度的计算 当下气泡在第一次收缩达到最小体积时，气泡统一理论计算的气泡半径为 ，和气泡统一理论的计算结果相比，网格精度为，和时的 计算结果相对误差分别为，和 从图（）中可以看出当网格精度从 提高到时，下气泡等效半径变化很小，说明用的网格计算满足计算精度要求（）气泡等效半径变化对比（）流场压力对比（）（）图 计算结果与气泡统一理论计算结果对比 注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同 的计算结果和气泡统一理论在上气泡和下气泡的第一个周期吻合较好，从第二个计算周期相差开始明显，并且 计算的气泡等效半径略小于气泡统一理论，这主要是由于 计算时的能量耗散较大 在图（）中上气泡冲击波阶段和下气泡冲击波阶段两种计算方法的结果几乎完全一致，包括下气泡的冲击波在上气泡表面的反射波扰动细节都很一致 气泡载荷除在时间上两种方法有差别，在载荷峰值和脉宽上相差很小 除了与理论方法进行对比，本文还将数值计算结果和试验进行对比，从气泡形态上验证数值计算模型 试验工况选自 等关于异相爆炸的论文，试验是在一个 的水箱中进行的 试验中用的炸药当量为 ，转换成 需要乘以转换系数，等效为 试验中上气泡爆炸水深为，相位差 当 ，下气泡膨胀到最大体积时上气泡产生 如图 中 所示，在上气泡开始膨胀后，冲击波迅速到达下气泡表面 冲击波在下气泡表面发生反射，并产生空化 随后下气泡开始产生向下的射流，上气泡开始向下凸起变化成“蛋”形 当 时下气泡已经被射流击穿，上气泡向上的射流开始形成 时上气泡被射流击穿，并且上气泡的射流宽度相较于下气泡很小 在随后的运动中，两个气泡朝着相反的方向运动 第 期 唐皓，等：水下爆炸异相气泡动力学特性的 有限元数值模拟研究图 计算结果与试验结果的气泡形态对比 应 用 数 学 和 力 学 年 第 卷从图 中试验和数值计算的结果对比可以看出，在上气泡和下气泡的全周期过程中数值计算结果和试验结果都基本一致 本文建立的 有限元数值计算模型在和理论方法和试验结果的对比中都得到了验证 结果与讨论 自由场单气泡和异相双气泡对比异相双气泡之间的相互作用过程是十分复杂的，在分析相位等因素对异相爆炸气泡的影响之前有必要对自由场单气泡和异相双气泡进行对比分析 计算工况为 在 水深爆炸，气泡距离参数，相位差 ，压力测点为（，）如图（）所示，在下气泡开始膨胀后上气泡等效半径开始迅速减小，减小的速度明显快于单气泡的情况 下气泡的膨胀速度也比单气泡的情况快，并且下气泡的最大体积也比单气泡时更大 结合图（）中气泡的迁移过程分析，在下气泡开始膨胀后上气泡被迅速推向上方，并且下气泡受上气泡的影响被推向相反的方向 力是超声场中声波作用下气泡受到的附加力 下气泡不仅运动方向发生了改变，气泡的射流方向也是向下的，这是因为气泡间的 力大于下气泡的浮力作用 从图（）中可以看出，受到上气泡达到最大体积后收缩的影响，下气泡的形态更长并且气泡的整体轮廓都比单气泡时大 气泡的伸长率为气泡水深方向的长度和气泡宽度的比值，单气泡在自由场膨胀的过程中气泡伸长率基本维持在，气泡保持球形 在异相爆炸中，延迟起爆的下气泡在膨胀过程中伸长率大于，气泡在膨胀时是“蛋”形而不是球形 异相双气泡的流场压力曲线相较于单气泡要复杂得多，这主要来自气泡表面对冲击波的反射作用 从图（）中可以看出，下气泡起爆后不久冲击波到达上气泡表面发生反射，并产生空化效应 此外异相爆炸时上气泡的气泡载荷也要小于单气泡的气泡载荷，这与上气泡的能量转移有关（）等效半径（）气泡迁移（）（）（）气泡形态（）流场压力（）（）图 异相爆炸双气泡和自由场单气泡对比 第 期 唐皓，等：水下爆炸异相气泡动力学特性的 有限元数值模拟研究关于气泡能量的计算参考自 等，的论文，气泡的总能量可由下式计算：（），（）式中 为流场动能，为气泡形心的垂向迁移，为气泡内部的平均内能 式（）中等号右侧括号内前两项为周围流体势能，括号内第三项为气体的内能 如图（）所示，炸药在起爆后内能出现极大地下降，这部分减少的内能主要转化为冲击波的能量 自由场单气泡在起爆之后，冲击波向周围传播，不存在冲击波和气泡之间的相互作用 异相爆炸时，冲击波在气泡之间的反复作用使气泡的能量变化更为复杂 在异相爆炸气泡的相互作用过程中，气泡间的内能变化与单气泡相比有着很大的差异 对于先起爆的上气泡而言，气泡的内能在气泡坍缩时的峰值与单气泡相比要小 延迟起爆的下气泡的内能相对于单气泡而言峰值和脉宽都大幅提高 总体上看，异相爆炸中上气泡和下气泡的内能相较于单气泡都增加了 图（）中，在下气泡起爆后，上气泡和下气泡的势能都增加了，下气泡相较于上气泡增加的幅度更大 图（）所示为自由场单气泡流场动能和异相双气泡流场动能变化，异相爆炸气泡流场动能为双气泡的总流场动能 在上气泡达到最大体积时流场动能接近于零，此时下气泡起爆流场动能瞬间增大，但增加的动能大于单气泡时气泡起爆增加的动能 结合图（）分析气泡的总能量变化，在炸药起爆后总能量短时间内衰减的能量为冲击波的能量 只有下气泡存在时气泡系统损失的冲击波能量约为，异相爆炸时下气泡起爆后气泡系统损失的冲击波能量约为 异相爆炸时冲击波带走的能量明显小于单气泡时冲击波带走的能量，这说明下气泡起爆后冲击波对异相气泡系统做功，使异相气泡系统的动能和总能量增加，从而导致相应的气泡势能和内能增加（）内能（）势能（）（）（）动能（）总能量（）（）图 异相爆炸双气泡和自由场单气泡的能量变化 相位对异相双气泡耦合的影响相位是影响异相双气泡相互作用的关键因素，本文针对垂向异相双气泡问题将相位差分为，应 用 数 学 和 力 学 年 第 卷，以及 共 种情况来讨论 计算工况为 在 水深爆炸，气泡距离参数 上气泡和下气泡在不同相位下的迁移曲线如图 所示，图例标签后面的箭头表示的是气泡在第一周期的射流方向 通过对比可以发现，上下气泡的迁移和相位差并不是线性关系 当 的值为正时表示上气泡先起爆，为负时表示下气泡先起爆 通过之前的分析已经知道了冲击波对气泡做功是增加异相爆炸气泡系统能量的主要原因，当 时上气泡体积最大，下气泡此时起爆产生的冲击波对上气泡的作用面积最大，因此图（）所示 时上气泡运动得最远 当 时下气泡体积最大，因此图（）所示，时下气泡运动得最远 对于射流方向，由于当 时先起爆的气泡已经产生很大位移，两个气泡之间的相互作用很弱，因此 不作为射流方向的参考 除 外，只有 时上下两个气泡产生对向射流，其余参数下均产生相反方向的射流 这是因为 时上下两个气泡近似镜像，气泡的运动和刚性壁面附近的气泡运动相似，产生朝向壁面的射流 后起爆的气泡在已经存在的气泡附近运动和自由面附近的气泡运动类似，产生反向射流（）上气泡（）下气泡（）（）图 不同相位异相爆炸双气泡的迁移 在异相爆炸双气泡的耦合过程中，异相气泡的周期相对自由场单气泡会发生改变 如图 所示，如果下气泡先起爆的话，在上气泡起爆后受上气泡干扰下气泡会提前坍缩，并且随着相位差的增大下气泡坍缩得越来越早，周期也越来越短 反过来，上气泡受到下气泡收缩时的吸引力而被拉长，上气泡的周期也越来越长 如果上气泡先起爆，在下气泡起爆后上气泡受到下气泡的作用会提前坍缩，并且相位差越小上气泡的周期越小 反过来，下气泡受到上气泡收缩时的吸引力而被拉长，下气泡的周期随着相位差减小也越来越长 相位差为 是一种特殊的情况，上下气泡的周期均大于自由场单气泡的情况 以上规律可以总结为：先起爆的气泡周期受后起爆气泡影响会缩短，相位差绝对值越小周期越小 后起爆的气泡由于受到先起爆气泡收缩的影响周期变长，相位差绝对值越小周期越大 相位对异相双气泡耦合过程的影响也包括对双气泡系统总能量的影响，如图 所示 在自由场水下爆炸中，爆炸冲击波向周围流场传播会带走相当大一部分能量 时，上下气泡同时起爆，在初始时刻上下气泡体积都很小，爆炸冲击波的作用面积有限，因此大部分冲击波穿过气泡散失到周围流场 时，初始时刻上下两气泡总能量为，冲击波过后系统总能量降至，损失的冲击波能量为，占比 当相位差不为零时，异相气泡系统总能量在先起爆的第一次冲击波过后从 降至，损失的冲击波能量为，占比，参考同相爆炸理论上第二次爆炸将会损失 但实际上当 时，后起爆的爆炸冲击波带走的能量约为；当 时，后起爆的爆炸冲击波带走的能量约为；当 时，后起爆的爆炸冲击波带走的能量约为 这部分冲击波少带走的能量通过冲击波和反射波对气泡做功留在异相气泡系统总能量内 并且后起爆的气泡越接近 ，气泡总能量损失得越少 第 期 唐皓，等：水下爆炸异相气泡动力学特性的 有限元数值模拟研究图 不同相位异相爆炸双气泡的周期图 不同相位异相爆炸双气泡系统的总能量变化 距离参数对异相双气泡耦合的影响在异相爆炸气泡相互作用的过程中，除了相位的影响，距离参数也是影响耦合过程的重要因素 本小节计算工况为 在 水深爆炸，气泡相位差 ，下气泡在上气泡体积最大时起爆 种距离参数下的异相爆炸气泡动态特性如图 所示，每种距离参数下气泡的动态过程按照时间顺序从左到右排列 当 时，下气泡在起爆后很快就被坍缩的上气泡吸入，并且下气泡顶部在上气泡内部破裂，上下气泡内部连通，破裂之后的液滴和液膜撞击到上气泡顶部将上气泡撕裂（）（）应 用 数 学 和 力 学 年 第 卷（）（）图 不同距离参数异相爆炸双气泡动态特性 图 不同距离参数异相爆炸双气泡系统的总能量变化 图（）所示异相气泡耦合过程为穿透融合模式，与此作为区别的是图（）图（）所示为 时异相双气泡的动态过程，当 时，上气泡刚要被射流击穿，此时两个气泡并未相连通，在随后的坍缩过程下气泡顶部突起的地方被环向射流分成上下两部分，下面的主体部分产生很细的向下的射流 图（）这种一个气泡进入另一个气泡内部但并未出现融合的情况定义为进入未融合模式 当距离参数继续增大时，如图（）和图（）所示，双气泡在脉动的过程中始终保持一定间隔，本文定义为未进入模式 除穿透融第 期 唐皓，等：水下爆炸异相气泡动力学特性的 有限元数值模拟研究合模式外，距离参数越小后起爆的气泡被拉伸得越长，产生的射流也越细 图 所示为不同距离参数时，异相爆炸气泡系统总能量的变化过程，在下气泡起爆之前 种距离参数下的气泡能量变化是完全相同的 当下气泡起爆后，在冲击波的传播过程中，异相双气泡系统总能量变化开始不同 距离参数越小，冲击波带走气泡系统的能量越少 时，由于上气泡被穿透撕裂，撕裂之后会有很多小的气泡向流场辐射压力波，因此气泡能量随后相较于其他参数会出现较大下降 结 论本文基于 和 方法建立了水下爆炸异相气泡动力学模型，对垂向异相双气泡非线性耦合问题进行了研究 本文通过将计算结果和气泡统一理论以及异相爆炸试验进行对比，验证了数值计算模型的有效性 异相爆炸气泡的脉动、坍缩和穿透融合等特性都得到了很好地模拟，并通过对比分析得到以下结论：）与自由场单气泡相比，异相爆炸中先爆炸的气泡会提前坍缩周期会缩短 冲击波在气泡表面会发生反射并产生空化，受异相爆炸气泡干扰气泡的射流方向和自由场也不同 气泡起爆后冲击波对异相气泡系统做功，使异相气泡系统的动能和总能量增加，从而导致相应的气泡势能和内能增加）当相位差不为 时，异相双气泡会产生反向射流，先产生的气泡周期会缩短，相位差绝对值越小先产生的气泡周期越短，后产生的气泡周期越长 当相位差为 时，异相双气泡产生对向射流，上下气泡的周期都会变长 相位差越接近，冲击波带走的能量越少，气泡总能量损失得越少）根据不同距离参数下的异相爆炸气泡耦合现象，本文将异相爆炸气泡分为 种模式，即穿透融合模式、进入未融合模式和未进入模式 除穿透融合模式外，距离参数越小后起爆的气泡被拉伸得越长，产生的射流也越细 距离参数越小，冲击波带走气泡系统的能量越少 参考文献（）：贺铭，张阿漫，刘云龙 近场水下爆炸气泡与双层破口结构的相互作用 爆炸与冲击，（）：（，（）：（）张阿漫，王诗平，彭玉祥，等 水下爆炸与舰船毁伤研究进展 中国舰船研究，（）：（，（）：（）贺啸秋，熊永亮，徐顺，等 底部加热肥皂泡上准二维湍流的数值模拟 应用数学和力学，（）：（，（）：（）孙涛，庞明军，费洋 气泡间距对受污染球形气泡界面性质和尾流的影响 应用数学和力学，（）：（，（）：（），：，：，（）：杜志鹏，张磊，谌勇，等 泡沫覆盖层对水下爆炸气泡射流防护机理缩比试验研究 应用数学和力学，（）：（，（）：（）彭玉祥，张阿漫，薛冰，等 强冲击作用下舰船结构毁伤的三维无网格 方法数值模拟 中国科应 用 数 学 和 力 学 年 第 卷学：物理学 力学 天文学，（）：（，：，（）：（），（）：李帅 水下爆炸气泡射流载荷特性研究 博士学位论文 哈尔滨：哈尔滨工程大学，（：，（）张阿漫，曾令玉，王诗平，等 水下爆炸气泡融合动态特性研究 应用数学和力学，（）：（，（）：（），（）：，：，：，：，：，（）：，：，（）：，：，：，（），（）：，：，：，：，（）：，（）：，（）：，（）：，（）：，：，：，（）：，第 期 唐皓，等：水下爆炸异相气泡动力学特性的 有限元数值模拟研究 ，：，：，：，：，：，（）：，（）：田昭丽 改进的欧拉有限元方法及水下爆炸近场载荷特性研究 博士学位论文 哈尔滨：哈尔滨工程大学，（：，（），：，（）：，（）：应 用 数 学 和 力 学 年 第 卷