水中高压脉动气泡与浮体流固耦合特性研究1)

胡振宇 曹卓尔 李 帅 张阿漫

(哈尔滨工程大学船舶工程学院,哈尔滨 150001)

引言

脉动气泡广泛地存在于自然界中,并且被应用于水下爆炸[1-2]、空化空蚀[3-4]、超声清洗[5]、海底资源勘探[6-7]和医疗[8-9]等领域.迄今为止,前人对气泡动力学的研究主要集中于自由场[10]、壁面[11-12]和自由液面[13-14]等简单边界条件下的气泡;对复杂边界下的气泡尤其是自由液面−浮体−气泡三者非线性耦合作用的研究相对匮乏.气泡与近水面结构物的非线性耦合作用是海洋工程和军事研究领域的重要科学问题.浮体在气泡脉动压力的影响下运动并伴随着自由液面大变形,这种运动响应反过来又对气泡动力学特性产生影响.目前气泡−浮体耦合作用中所蕴含的力学机理仍未完全揭示,因此本文针对该问题进行数值研究,同时开展放电气泡实验验证数值模型,旨在增强对气泡与近水面结构物非线性流固耦合效应的认识,为相关领域的研究和应用提供基础性支撑.

Rayleigh[10]最早基于球状假设提出著名的Rayleigh-Plesset(RP) 方程,随后Herring[15]和Keller等[16]在RP 方程的基础上进一步考虑流体的黏性和可压缩性的影响,得到各种不同形式的改进RP 方程.近边界(自由液面、刚/弹性壁面、其他气泡界面)脉动气泡往往无法保持球状,这种非球型脉动极大地增加了理论研究的难度.在实验研究方面,研究者使用水下爆炸实验[17-19]、放电气泡实验[20-21]和激光气泡实验[22-23]来研究气泡与各种边界的耦合作用,高速摄影技术也被用于捕捉气泡的动力学行为.Lauterborn 和Bolle[24]使用激光聚焦产生单个气泡以研究气泡与刚性壁面的耦合作用,发现在气泡坍塌后期产生指向壁面的高速射流,其速度可达120 m/s.高速水射流与结构表面的砰击被认为是空蚀的重要成因[11,25].然而由于实验摄影中时空分辨率的限制,精细地捕捉射流内部形态仍然是一个挑战.通过数值模拟,能得到更多的流场信息以分析气泡与流场边界的非线性耦合作用.

Chahine 等[26]最早针对放电气泡与柱形浮体的流固耦合作用进行了实验和数值研究,分别使用刚固浮体与可动浮体进行研究发现浮体的运动响应对气泡坍塌模式、射流形成和气泡脉动周期具有重要影响.Klaseboer 等[27-28]对水下爆炸气泡与漂浮船体的耦合作用进行了数值研究,发现气泡射流方向受到自由液面、船体表面曲率和重力的影响.在数值计算中使用了镜像法来模拟自由液面,该方法对于γf>1.4(γf为无量纲气泡与自由液面距离,取气泡最大半径为特征长度) 的工况,精度较好,但对于近距离的强耦合工况(γf<1.4) 需要考虑全非线性自由液面条件.Zong 等[29]使用边界积分法(BIM)联合有限元法建立了气泡与刚性/可变形浮体耦合模型,考虑全非线性自由表面条件以及结构的运动和弹性变形,并采用传统的差分法求解结构表面压力(伯努利方程中∂ϕ/∂t项),然而这种方法可能不适用于结构运动很快的强耦合工况[30-31].Li 等[30]首次将辅助函数法应用于气泡与水中悬浮球体的非线性耦合效应研究中,相比传统差分法,计算更为准确,本文将该方法拓展至气泡与浮体的流固耦合特性研究中.近年来,有限体积法[32]、有限元法[33]等方法也被成功的应用于气泡动力学的数值研究中,相比于这些域离散方法,BIM 仅需对流场边界进行离散,这种降维处理使BIM 计算效率具有极大的优势.

本文联合BIM 与辅助函数法建立脉动气泡与浮体全非线性耦合模型.首先开展了一系列放电气泡与柱形浮体耦合作用机理实验,然后对比数值与实验结果,验证计算模型与方法的有效性.采用本文的数值模型系统地研究无量纲气泡与浮体距离参数γs对气泡动力学特性和结构响应的影响规律.本文的研究内容安排如下;第二章给出本文数值模型的基本理论;第三章首先验证数值模型然后分析讨论实验与数值结果,主要关注气泡射流模式、射流速度、气泡迁移以及浮体运动响应;第四章给出本文研究的结论.

1 基本理论和数值模型

1.1 边界元法

本文采用基于势流理论的轴对称边界元法,建立脉动气泡与浮体非线性流固耦合作用数值模型.基于流体无黏、无旋、不可压缩假设,流场(水域)可由拉普拉斯方程描述

上式可转化为积分形式的边界积分方程[34-35]

式中,p为控制点,q为源点,A为p点观察流场(水域)的立体角,S为包括气泡表面Sb、自由液面Sf、浮体湿表面Sw以及无穷大球面S∞的流体边界,n表示边界法向,以指向流场外为正.图1 所示为气泡与浮体非线性流固耦合作用数值模型示意图,引入柱坐标系O(r,θ,z),其坐标原点位于初始气泡中心,z轴指向重力的反方向且为对称轴,定义浮体的入水深度为d,初始气泡中心至初始自由液面的距离为h;本文考察圆柱型浮体与气泡的耦合作用,圆柱浮体半径为Rs,气泡最大等效半径为Rm.

图1 气泡与浮体非线性流固耦合作用数值模型示意图Fig.1 Sketch of the numerical model for investigating the nonlinear fluid-structure interaction between a bubble and a floating structure

离散气泡、浮体和自由液面边界,边界积分方程转化为一个线性方程组,求解此线性方程组后可获得气泡和自由液面节点的法向速度与浮体湿表面节点速度势,通过有限差分可获得气泡与自由液面的切向速度从而合成总速度[36].

对浮体附近脉动气泡动力学特性的时域模拟需显式更新边界表面节点位置和速度势.给出气泡表面与自由液面满足的运动学边界条件如下

其中x为节点位置矢量.

气泡表面和自由液面满足的动力学边界条件可由非定常伯努利方程给出,分别为

式中,p∞为z=0 平面内流场中无穷远处静水压力,pb为气泡内压,g表示重力加速度,ρ 为流体密度.气泡与自由液面满足的动力学边界条件用于更新节点速度势.由于气泡内部不可冷凝气体在脉动过程中满足理想气体绝热定律[37],忽略表面张力作用,pb可由下式计算

式中,pv为气体的饱和蒸汽压,本文取2337 Pa;p0为初始气泡内部不可冷凝气体压力,V0为初始气泡体积,V为瞬时气泡体积;κ 为气泡内部不可冷凝气体的比热比,本文取κ=1.4.

速度势ϕ 需满足浮体湿表面不可穿透条件如下

式中,n表示边界单位法向量,Vs表示浮体运动速度.

1.2 辅助函数法

通过牛顿第二定律

可求得浮体加速度a.式中m为浮体质量,G为浮体重力,空气压力Fa通过在浮体干表面上对压力(大气压patm) 积分即可获得,水动力Fh通过伯努利方程获得[29]

传统方法通过对相邻时间步的节点速度势采取向前差分近似求解∂ϕ/∂t项,但是时间步过小或浮体运动剧烈时,精度不足.理论上这种方法得到的是dϕ/dt,可通过物质导数与局部导数的关系得到修正的差分法即∂ϕ/∂t=dϕ/dt−Vs·∇ϕ.本文采用的辅助函数法可精确求解结构响应,后面将对比修正后的传统差分法、未修正的差分法和辅助函数法的计算精度.

本文中浮体仅存在z向运动,作用在浮体上的外力沿z轴,故可将式(8)改写为

上式左侧第一项代表惯性力项(其中积分项代表附加惯性力),左侧第二项代表重力,χ 与ξ 为两个辅助函数,Sd为浮体与大气接触的表面.

在气泡表面与自由液面,χ 满足边界条件[30]

在结构湿表面,χ 满足边界条件如下[38]

ξ 在气泡表面与自由液面满足的边界条件可分别表示为[30]

ξ 在浮体湿表面满足边界条件[31]

辅助函数χ 与ξ 满足拉普拉斯方程,并已知其在气泡和自由液面满足的狄利克莱边界条件与在浮体湿表面上满足的诺伊曼边界条件,采用类似于求解边界积分方程的方式可以得到在浮体湿表面上的值,代入式(10)中可求解出浮体加速度.

1.3 环状气泡模型

射流砰击时,原来的单连通气泡转变为双连通气泡,即变为环状气泡.此时气泡射流砰击处的速度势不连续.为保证计算的稳定性,需要解决砰击处节点的速度势突跳问题.本文采用Wang 等[35]提出的涡环模型,在气泡内部布置一个涡环,其环量为砰击处的速度势突跳.将气泡表面速度势分解为[34]

其中,ϕvor为涡环诱导速度势,ϕres为残余速度势.残余速度势ϕres在环状气泡表面连续,采用边界积分法进行计算,涡环诱导速度势采用一种半解析法求解.详细的计算过程可参考文献[39],本文不再赘述,这里仅对涡环模型的速度势分解思想进行介绍.

1.4 三相接触模型

对于水面漂浮结构存在气−液−固三相接触线,如图2 所示,接触线需要特殊处理以同时满足自由液面狄利克莱边界条件与浮体表面诺伊曼边界条件[40].本文采用双节点法追踪交界线,这种方法同样适用于气泡与结构接触问题.将接触点k分解为k1和k2两点,分别属于自由液面和柱体;使用下标fk1和sk2分别表示属于自由液面、浮体的交界点.下标b,f 和s分别表示除了交界点外气泡、自由液面和浮体上的节点.下标n表示速度势的法向导数.

图2 三相接触模型示意图Fig.2 Sketch of the solid-liquid-gas contact model

离散边界积分方程(2) 后可得到如下线性方程组

由于k1和k2两点在空间位置上完全相同,影响系数矩阵H和G第3 和第4 行完全相同,整理可得

求解线性方程组后可得到k1点的法向速度.

当气泡产生位置距浮体表面很近时,由于气泡膨胀或者气泡朝壁面迁移,导致气泡附着到浮体表面.若气泡节点与浮体节点距离小于单元尺度,传统BIM 产生严重的数值不稳定性;如图3(a) 所示,部分气泡节点穿过结构表面,产生了计算错误.因此本文将气泡与浮体间的薄水膜去除,连接气泡与浮体节点,进而保证数值计算的稳定性,如图3(b) 所示.Wang 等[41]的研究表明,这层薄水膜压力趋近气泡内压,将其去除未对气泡的动力学特性产生可见的影响.节点去除的范围和接触判定时刻由气泡最小单元长度∆s控制,当气泡节点与浮体节点最小距离小于∆s时,即认为气泡即将附着在浮体表面,并将距离小于2.5∆s[38]的气泡与浮体节点删除并连接气泡与浮体,采用双节点法处理产生的气−液−固三相交界线.气泡附着在浮体后,气泡内压直接作用在浮体表面上,此时浮体运动方程如下

式中,Sg为浮体与气泡接触表面.

图3 传统BIM 计算结果(a)与采用三相接触模型计算结果(b)对比Fig.3 Comparison between the traditional BIM(a)and the gas-liquid-solid contact model(b)

1.5 无量纲化

本文将所有物理量进行无量纲化,并使用无量纲量进行计算与讨论.取Rm,∆p=p∞−pv和ρ 分别作为长度、压力和密度的特征量,其他物理量如时间、速度、速度势和加速度的无量纲特征量分别为Rm(ρ/∆p)1/2,(∆p/ρ)1/2,Rm(∆p/ρ)1/2和∆p/(ρRm)[34-35,37].本文后续出现的物理量如无特殊说明均为无量纲形式.

定义浮力参数δ

强度参数ε

初始气泡中心位置与自由液面距离参数γf

初始气泡中心位置与浮体距离参数γs

尺度比ζ

由于放电气泡实验中很难得到气泡内压,并且Turangan 等[42]的研究发现强度参数ε 在100 ∼500 范围内变化,气泡射流速度变化很小,气泡的动力学特性变化也很小.经过计算发现,本文选取ε=300,实验与数值结果吻合良好.计算中气泡初始半径与强度参数相关[43],本文选取R0=0.113.计算中所用时间步长由下式计算得到[37]

其中C为单步最大的速度势改变量.

2 结果分析与讨论

2.1 收敛性分析

为了验证本文数值模型的正确性和有效性,首先进行收敛性分析,包括网格与时间步的收敛性分析.分别使用气泡节点数量N等于80,120,160,200和240 的模型进行计算.图4(a) 展示了不同网格密度下射流砰击瞬时气泡形态,节点数为80,120 和160 的结果中气泡形态具有一定的差异,当节点数N200 时气泡轮廓基本重合,证明当气泡网格节点取200 时能获得收敛的结果.

本文采用四阶Runge-Kutta 法对脉动气泡与浮体相互耦合的时域运动过程进行模拟,显式算法的收敛性依赖于时间步长.取C分别为0.04,0.02,0.01和0.005 分别进行计算.图4(b) 描述了浮体加速度随时间的变化,4 种方案取得了一致的结果,进一步证明了本文提出的数值模型的有效性.考虑程序计算的精度以及所需的计算成本,本文后续计算选取N=200,C=0.01.

图4 模型的收敛性分析(a)计算网格的收敛性(b)计算时间步的收敛性Fig.4 Convergence study with respect to(a)the mesh and(b)the time step

2.2 典型工况对比分析

本节将分别对3 个典型距离参数γs下气泡与浮体非线性耦合作用进行研究和讨论,展示实验中的物理现象,利用数值模拟分析现象背后的力学机理.本文实验方法详见已发表文献[31,44].

2.2.1 弱耦合工况

第一个实验中初始气泡深度h=74.7 mm,初始浮体入水深度d=45.0 mm,浮体半径Rs=16.0 mm,气泡最大等效半径Rm=16.8 mm.根据实验设定计算初始参数为:ε=300,R0=0.113,κ=1.4,γf=4.45,γs=1.77,ζ=0.95,δ=0.04.图5 所示为本工况数值与实验结果对比图,在此对本文的数值与实验对比图进行说明;图中左侧为高速摄影捕捉的气泡图像,右侧为对应时刻气泡运动模拟结果.每一张图片左上方显示了实验中气泡脉动时间,右侧数字为数值计算结果对应的气泡脉动时间,图5(a) 左下方给出了比例尺,图5(b) 中箭头指出了实验图像中自由液面位置,为了更好地分析浮体与气泡的非线性耦合效应,数值结果中还给出无量纲压力与速度场.

图5(a)中,气泡在内部高压的推动下猛烈膨胀.图5(b) 所示为气泡膨胀达到其最大体积时刻,此时气泡近似球状,自由液面和浮体未对气泡脉动产生明显的影响.由于过度膨胀气泡周围压力比静水压力小,在接下来的脉动过程中在流体惯性的作用下气泡开始坍塌.浮体对气泡界面收缩具有延时作用,在气泡坍塌前期,气泡顶部略微凸起,而底部具有较大坍塌速度因此表面较为平坦,如图5(c)所示.这种非球状坍塌导致气泡顶部在坍塌中期获得更大的曲率,如图5(d)所示.在气泡坍塌后期,气泡周围流体压力迅速增大,但气泡顶部压力梯度更大,使顶部界面向下(逆压力梯度方向)快速收缩,如图5(e)所示.第一个周期结束,并未形成射流,而在回弹阶段气泡顶部产生远离浮体的射流并穿透另一侧气泡表面,如图5(f) 所示.气泡顶部射流现象与Lauterborn[45]的预测一致,即气泡具有高曲率的部分更容易产生射流.此工况中气泡产生背离浮体的射流(称为“反射流”),说明浮体对气泡的吸引力较小,无法诱导指向壁面的射流,故将本工况定义为弱耦合工况.总体上,数值计算得到的气泡运动特性与实验吻合良好.

图5 弱耦合工况数值与实验结果对比Fig.5 Comparison between the numerical results and the experimental images for the weakly coupling case

2.2.2 中度耦合工况

第二个实验中气泡产生位置距离浮体更近,耦合作用更强,气泡初始深度h=60.0 mm,初始浮体入水深度d=42.6 mm,浮体半径Rs=16.0 mm,气泡最大等效半径Rm=16.4 mm.依据实验设定计算初始参数如下:ε=300,R0=0.113,κ=1.4,γf=3.66,γs=1.06,ζ=0.98,δ=0.04.

图6 给出了本工况中几个典型时刻数值与实验结果对比.铜丝燃烧生成一个高压气核,在内部高压气体作用下迅速膨胀,如图6(a) 所示.图6(b)表示气泡达到最大体积时刻,气泡上方流体运动受到浮体的阻碍,气泡上表面被压平,此时气泡上方流体继续向上运动推动浮体远离气泡.在气泡坍塌前期,远离浮体的气泡底部近似球形坍塌,如图6(c)所示.柱体阻碍了流体从顶部流向气泡中心,而对侧方的流体的运动影响较小,因此形成了“梨”形气泡.图6(d)∼图6(f)描述了气泡坍塌中后期的运动,气泡两侧以及底部收缩的更快,顶部运动速度很小,“梨”形气泡逐渐转变为“梭”形气泡,气泡在垂向被拉长.在坍塌后期气泡底部压力梯度变大,产生高压区,使底部气泡表面快速收缩形成射流.图6(f)的实验图像中气泡底部留下铜丝燃烧产物,表明形成了一束向上的射流.此工况中浮体对气泡的动力学行为影响更大,并能诱导气泡产生指向壁面的射流(称为“正射流”),因此将此工况命名为中度耦合工况,从实验中也能观察到浮体在膨胀阶段被气泡排斥而在收缩阶段被吸引的现象.

图6 中度耦合工况数值与实验结果对比Fig.6 Comparison between the numerical results and the experimental images for the intermediately coupling case

2.2.3 强耦合工况

为研究浮体与脉动气泡强烈的非线性耦合作用,第三个实验进一步减小了浮体与气泡间的距离.此时气泡在膨胀阶段附着到浮体底面,具有明显的流固耦合效应,固将此工况命名为“强耦合工况”.实验中参数如下;初始气泡深度h=42.6 mm,浮体入水深度d=36.0 mm,浮体半径Rs=14.0 mm,气泡最大等效半径Rm=13.1 mm.依据实验设定计算初始参数如下:ε=300,R0=0.113,κ=1.4,γf=3.25,γs=0.5,ζ=1.07,δ=0.04.

图7 给出了强耦合工况中几个典型时刻实验与数值计算结果对比.如图7(a)所示,在气泡膨胀初期,由于柱形浮体的阻碍作用,气泡顶部变平.在膨胀阶段气泡上表面逐渐靠近浮体底面,需启动三相接触模型以保证数值稳定性.图7(b)所示为气泡与浮体表面接触时刻,气泡即将附着到浮体表面.图7(c)表示气泡最大体积时刻,此时气泡顶部仍有一定的膨胀速度,而底部与侧面在收缩,进一步说明浮体对气泡运动具有延时的作用.侧面流体在气泡坍塌时能自由的流向气泡中心,顶部流体受到阻碍,如图7(d)所示.在接下来的坍塌过程中,气泡底部产生高压区,驱动射流直接砰击在浮体上.实验中还观察到射流前进阶段气泡撕裂为顶部大气泡和底部微小气泡,计算中并未产生撕裂现象,如图7(e)和图7(f)所示.图中最后两帧数值与实验结果具有较大偏差,这可能是由于气泡脉动过程中浮体底面附着的空化气泡(见图7(c))与放电气泡产生耦合作用所致.并且三相接触线附近气泡表面产生翻卷和破碎的现象,这可能与固体材料的润湿性能相关,而本文计算中未考虑接触角的影响.总体上气泡形态与周期与实验结果吻合良好.值得注意的是,柱形浮体在此强耦合工况中产生了显著的运动.

图7 强耦合工况数值与实验结果对比Fig.7 Comparison between the numerical results and the experimental images for the strongly coupling case

图8 给出了强耦合工况中浮体位移和浮体速度时历曲线,分别对比了本文辅助函数法模型、传统修正后的差分模型以及实验结果.实验结果中误差棒的中心代表4 次测量的平均值,误差棒的范围代表标准差;菱形标记表示气泡最大体积时刻.由图8(a)可知,浮体在气泡膨胀阶段具有向上的运动,由于惯性作用气泡达到最大体积后浮体仍继续向上运动.在坍塌阶段浮体向下加速运动再逐渐减速.可以看出,相比于传统方法本文采用的辅助函数方法与实验结果更加一致.在气泡膨胀早期,数值结果中浮体运动较快,这是由于数值计算中忽略了铜丝燃烧放热的过程;而在后期随着时间步减小,传统差分法结果逐渐偏离实验结果.从图8(b) 可以看出在气泡膨胀前期,两种方法求得的结构速度基本一致,而在后期传统差分法精度降低,在气泡附着于浮体表面后(对气泡拓扑结构进行“数值切割”可能进一步增大了不稳定性),速度时历曲线出现剧烈震荡(图中虚线表示气泡与浮体接触时刻),而本文方法能得到稳定且准确的结果.故对于研究近距离的脉动气泡与浮体非线性耦合作用,本文提出的数值模型更为有效.

图8 柱形浮体(a)位移与(b)速度时历曲线(菱形标记代表最大体积时刻,虚线表示气泡与浮体接触时刻)Fig.8 Comparison of(a)the displacement and(b)the velocity histories of the cylindrical floating structure(the diamonds denote the moment of reaching the maximum bubble volume,the dashed line denotes the moment of the bubble-structure attachment)

2.3 无量纲参数讨论

2.3.1 γs对气泡射流模式的影响

影响气泡与浮体非线性耦合效应的参数众多,主要有无量纲气泡与浮体距离参数γs、无量纲气泡与自由面距离参数γf与尺度比ζ.众多的影响参数导致了这一问题的复杂性.限于篇幅,本文仅研究γs的影响,其余参数与实验参数一致,本文提出的数值模型也可研究其余参数的影响.依据实验设定后续数值研究的初始参数:ε=300,R0=0.113,γ=1.4,d=2.6,ζ=1,δ=0.04,γs=0.2 ∼2.

图9 展示了不同距离参数γs下,射流砰击时刻的气泡形态,明显气泡射流模式分为5 种类型,分别为颈缩型环状射流(0.2γs0.3)、接触射流(0.4γs0.6)、非接触射流(0.7γs1)、对射流(1.1γs1.3)和反射流(1.4γs2),图中虚线代表浮体.

图9 不同无量纲气泡与浮体距离γs 下气泡射流模式Fig.9 Jetting patterns for different dimensionless distances γs

图10 给出了γs=0.2(颈缩型环状射流)工况中典型时刻的压力和速度场.气泡收缩时牵引浮体侧面流体向下运动,同时水平方向气泡能自由收缩,如图10(a)所示.水平方向流体与垂向流体在柱体底部最大半径处汇聚,在此处形成高压区(见图10(b) 和图10(c)).在压力梯度作用下流体向斜下方运动,形成颈缩现象.类似于Cui 等[44]的分析,浮体的向下运动可能进一步改变流体运动,使其偏离形成环状向外的射流导致气泡撕裂为一个大单连通气泡与小环状气泡.

图10 γs=0.2 工况中典型时刻压力与速度场(分别对应最大体积时刻t=0.910,气泡坍塌时刻t=1.132,环向射流砰击时刻t=1.403)Fig.10 Pressure and velocity fields at some typical moments for γs=0.2(corresponding times are t=0.910 at the moment of reaching the maximum bubble volume,t=1.132 during the bubble collapse stage and t=1.403 at the moment of the annular jet impact respectively)

图9 第3 ∼5 帧中气泡附着在浮体表面,此时产生的射流称为接触射流.图11 给出了γs=0.6(接触射流) 工况中典型时刻压力和速度场.在气泡坍塌后期,底部附近形成了局部高压区域,驱动气泡界面快速收缩形成高速细射流直接砰击在浮体上,如图11(a)和图11(b)所示.射流砰击到浮体底部后,在浮体底面产生很高的压力.流体受到壁面的作用变为径向流动,最终导致气泡脱离浮体表面,如图11(c)所示.在气泡脱离瞬时,浮体底面中心附近压力比射流砰击时刻更高,这是射流砰击产生的局部高压与气泡收缩产生的脉动压力叠加导致的结果.

图11 γs=0.6 工况中典型时刻压力与速度场(分别对应气泡坍塌时刻t=1.792，射流砰击时刻t=1.810 和气泡与浮体脱离时刻t=1.838)Fig.11 Pressure and velocity fields at some typical moments for γs=0.6(corresponding times are t=1.792 during the bubble collapse stage,t=1.810 at the moment of the jet impact and t=1.838 at the moment of the bubble-floating structure detachment respectively)

图9 第6 ∼9 帧气泡射流模式称为非接触射流.气泡底部形成的高速射流首先砰击在气泡上表面,穿过水膜后砰击在浮体底面.图12 给出了γs=0.9(非接触射流)工况中典型时刻压力和速度场.在坍塌阶段气泡下方压力梯度最大,导致正射流的形成.射流砰击在气泡上表面后,在气泡与浮体间形成局部高压区域,如图12(a)和图12(b)所示.环状气泡继续收缩同时向浮体运动,在顶部形成突出物最终与浮体底面接触.值得注意的是射流砰击后浮体底面中心附近压力比上文接触射流工况更大.环状气泡继续坍塌,壁面附近压力超过140 MPa,如图12(c)所示.

图12 γs=0.9 工况中典型时刻压力与速度场(分别对应气泡坍塌时刻t=1.806,射流砰击时刻t=1.807 和气泡与浮体接触时刻t=1.809)Fig.12 Pressure and velocity fields at some typical moments for γs=0.9(corresponding times are t=1.806 during the bubble collapse stage,t=1.807 at the moment of the jet impact and t=1.809 at the moment of the toroidal bubble-floating structure attachment respectively)

图9 第10 ∼12 帧展示了反射流与正射流相撞的现象(称为“对射流”).如图13(a)所示,在气泡坍塌末期,气泡顶部首先产生一个宽度较大的凹陷.气泡底部形成很高的曲率,在其下方产生一个高压区域.高压区使气泡底部迅速收缩形成一束指向壁面的正射流.最终正射流与宽而慢的反射流相撞(图13(b)和图13(c)).射流相互撞击后在撞击区域附近产生了很高的压力.反射流内部未与正射流砰击的部分继续向下运动,导致砰击区域环状凹陷的形成,如图13(d)所示.这种由射流相互撞击产生的环状凹陷无法发展成射流,在回弹阶段逐渐消失,如图13(e) 和图13(f) 所示.在坍塌和回弹阶段,气泡和高压区逐渐向浮体迁移.

图13 γs=1.3 工况中典型时刻压力与速度场(分别对应气泡坍塌时刻t=1.806,t=1.807,射流砰击时刻t=1.809,气泡坍塌时刻t=1.810 和气泡回弹时刻t=1.820,t=1.853)Fig.13 Pressure and velocity fields at some typical moments for γs=1.3(corresponding times are t=1.806,t=1.807 during the bubble collapse stage;t=1.809 at the moment of the jet impact;t=1.810 during the bubble collapse stage;t=1.820,t=1.853 during the bubble rebounding stage respectively)

图9 第13 帧中气泡产生了反向射流与径向射流,径向射流砰击将气泡撕裂为两个单连通的气泡.图9 第14 ∼15 帧中气泡只产生径向凹陷,无法穿透气泡.图14 给出了γs=1.5(反射流) 工况中典型时刻压力和速度场.在气泡坍塌后期,气泡顶部压力比周围压力更高导致气泡形成向下的凹陷.然而第一个周期内,向下的凹陷无法形成射流穿透气泡表面,如图14(a)所示.随后气泡两侧收缩较快形成径向的气泡凹陷.在气泡回弹阶段,上方凹陷发展为一束向下的反射流.由于气泡膨胀射流变得越来越细而径向凹陷开始消失,如图14(b)所示.在反射流尖端形成了非常奇特的压力分布,即在射流尖端具有较高的压力,而在两侧却形成较低的环形压力区域.反射流与下方气泡界面砰击后,在气泡底部形成了局部高压,而其余方向压力由于气泡过度膨胀降至较低水平.此时气泡呈现为“蝴蝶”形,如图14(c)所示.

图9 第16 ∼19 帧中气泡仅产生反射流.图15 给出了γs=1.9 工况中典型时刻的压力与速度场.与图14 工况相同,在第一个脉动周期内射流无法穿透气泡表面,而在回弹阶段射流变细并加速向下运动,如图15(a)和图15(b)所示.射流砰击产生后气泡下方形成局部高压,并且形成了向下的突出物,如图15(c)和图15(d)所示.

图14 γs=1.5 工况中典型时刻压力与速度场(分别对应气泡坍塌时刻t=1.803,气泡回弹时刻t=1.833,t=1.860)Fig.14 Pressure and velocity fields at some typical moments for γs=1.5(corresponding times are t=1.803 during the bubble collapse stage;t=1.833 and t=1.860 during the bubble rebounding stage respectively)

图15 γs=1.9 工况中典型时刻压力与速度场(分别对应气泡最小体积时刻t=1.816,气泡回弹时刻t=1.825,t=1.837,t=1.891)Fig.15 Pressure and velocity fields at some typical moments for γs=1.9(corresponding times are t=1.816 at the moment of reaching the minimum bubble volume;t=1.825,t=1.837,and t=1.891 during the bubble rebounding phase respectively)

图15 γs=1.9 工况中典型时刻压力与速度场(分别对应气泡最小体积时刻t=1.816,气泡回弹时刻t=1.825,t=1.837,t=1.891)(续)Fig.15 Pressure and velocity fields at some typical moments for γs=1.9(corresponding times are t=1.816 at the moment of reaching the minimum bubble volume;t=1.825,t=1.837,and t=1.891 during the bubble rebounding phase respectively)(continued)

当γs很小时(γs0.3),浮体底面两侧产生的局部高压使气泡产生颈缩,最后导致气泡撕裂.γs增大至(0.4γs1.3),气泡与浮体具有较强的流固耦合效应.由于Bjerknes 力[46-47]作用,在坍塌阶段气泡产生指向浮体的正射流,而当气泡与浮体较远时(1.1γs2.0),在坍塌阶段气泡顶部受到浮体的延时作用具有较高曲率并且在气泡上方压力梯度更大,使气泡上表面产生向下的反射流.自由液面对气泡具有排斥作用,进一步促进了反射流的形成[48-49].

2.3.2 γs对射流速度的影响

由于本文实验装置和条件的限制,很难精确测量射流速度,为验证本文模型能预测射流速度,首先对比本文模型计算结果与Cui 等[44]放电气泡与不同厚度浮冰耦合作用的实验.图16 所示为最大射流速度Vjet(正射流)与浮冰厚度的关系(γs=0.93).由图可知,当浮冰厚度大于28.1 mm 时,计算所得最大射流速度基本在实验测量范围内,而随着浮冰厚度减小或者射流速度增大,两者偏差增大,最大误差为16.8%(定义为实验与计算偏差/实验值).实验中由于气泡壁面作为一个“凹透镜”折射光线导致测量的射流速度可能偏小[43,50].并且当射流速度很大时,需要极高的拍摄帧率才能精确测量最大射流速度值,Cui等[44]的实验中拍摄帧率可能不足.这两种因素导致了数值计算结果相比实验结果偏大.浮冰厚度增大,射流速度减小,并且数值与实验结果吻合程度更高,一定程度上验证了本文计算模型能较准确的预测射流速度.本文后续进一步研究γs对脉动气泡与浮体非线性耦合效应的影响.

图16 最大射流速度数值与Cui 等[44] 实验结果对比Fig.16 Comparison of the maximum jet velocity between the present model and the experiment[44]

图17 给出了3.3.1 节中工况(γs=0.4 ∼2)最大射流速度Vjet和最大反射流速度Vcounter与γs的关系.正射流速度远高于反射流且随γs先增大后减小再增大,在γs=0.8 时达到最大值136.8(1368 m/s);反射流速度随γs增大而增大,最大值为27.2(272 m/s).图12 给出的γs=0.9 工况中,高速细射流砰击在另一侧气泡表面,导致浮体底面附近具有很高的压力(超过140 MPa).而当0.7γs0.9 时,射流速度均大于100,可能在浮体底部产生巨大的射流砰击载荷,造成结构局部损伤.与Lechner 等[51]对超近壁面附近脉动气泡动力学特性的数值研究结果相似,在气泡坍塌阶段与浮体底面平行的流动比垂直于底面的流动更快,导致气泡在远端形成很高的曲率,从而形成了速度约1000 m/s 的高速细射流(Lechner 等[51]发现的射流速度超过2000 m/s).本文数值模型中最大射流速度远大于Philipp 和Lauterborn[50]通过激光气泡与固定刚性板耦合作用实验测量的射流速度(150 m/s),可见浮体与脉动气泡的非线性流固耦合效应能显著影响气泡的动力学特性,在实际应用中需要考虑这种耦合作用的影响.但实际中如此高速的细射流可能存在不稳定性[51],并且流体可压缩性具有重要影响,这不在本文的讨论范围内.

图17 射流速度与无量纲气泡与浮体距离γs 的关系Fig.17 Maximum jet velocity versus the dimensionless bubble-floating structure distance γs

2.3.3 γs对气泡迁移和结构运动的影响

图18 展示了不同γs下气泡形心的位移时历曲线(计算至射流穿透前).当γs<1.5 时,在气泡膨胀阶段,由于浮体对气泡上方流体的运动具有很强的阻碍,导致气泡下表面相比上表面膨胀得更快,类似于气泡被浮体推开,形心远离结构.而在坍塌阶段气泡受到浮体Bjerknes 吸引[46]作用强于自由液面的Bjerknes 排斥[35,48]作用,气泡又被浮体吸回.在气泡坍塌阶段近浮体与远离浮体侧不平衡压力使气泡整体受到一个朝向壁面的力,从而导致了气泡朝浮体迁移,如图6(d)∼图6(f)所示.由图18 可知随着γs的增大,这种相互排斥、吸引作用减小.当γs1.5 时(此时γf4.1),气泡在膨胀阶段基本保持其形心位置不变,而在坍塌和回弹阶段气泡逐渐远离浮体.这可能是因为自由液面对气泡的排斥力比浮体的吸引力更强.

图18 不同γs 下气泡形心的位移时历曲线Fig.18 Time histories of the displacement of the bubble centroid for different γs

图19 浮体最大速度与无量纲气泡与浮体距离γs 的关系Fig.19 Maximum velocity of the floating body versus the dimensionless bubble-floating structure distance γs

3 结论

本文采用边界积分法联合辅助函数法,建立了高压脉动气泡与浮体非线性流固耦合数值模型,同时开展了柱形浮体与放电气泡流固耦合作用实验,对比数值与实验结果,二者吻合良好.采用该数值模型,本文对气泡与浮体无量纲距离参数γs进行了系统研究,得到的主要结论如下;

(1)将本文流固耦合模型、传统流固耦合模型(差分法计算速度势偏导数)与实验结果进行对比分析,发现传统方法在处理气泡与浮体接触工况时,精度明显不足,而本文的辅助函数法大大提高了计算稳定性与精度.

(3) 正射流速度随着γs增大先增大后减小再增大,而反射流速度随γs增大而增大.正射流速度能达到约1000 m/s 的量级,说明脉动气泡与浮体的流固耦合效应在实际应用中影响重大,其形成机理为气泡底部高曲率与浮体的Bjerknes 吸引力联合作用.该现象需要未来更精密的实验来验证.

(4) 总体而言,在气泡膨胀阶段,浮体排斥气泡,气泡形心向下运动;在坍塌阶段,当0.2γs1.4 时,浮体对气泡的Bjerknes 吸引力强于自由液面Bjerknes 排斥力,导致气泡向浮体迁移;当1.5γs2,浮体对气泡的Bjerknes 吸引力比自由液面Bjerknes 排斥力弱,导致气泡远离自由液面.