空化是一种由于局部压力震荡引起的空泡初生、发展及溃灭的物理现象，其往往与水力机械性能下降、振动噪声、空蚀破坏密切相关，是制约水力机械大型化、高速化发展的主要原因．然而如今，工业界也开始发展一些基于空化现象的新兴技术，如污水降解[1]、药物输送[2]、增材制造[3]、化工混合[4]、清洗除锈[5]等．以上领域无论是空化防治或是空化应用，都涉及到空泡在复杂边界环境下的演化特性，因此开展空泡边界相互作用研究对于空化在工业界中的应用具有非常重要的指导意义．近年来基于实验与数值模拟手段对复杂边界环境下的空泡动力学特性已开展大量研究，其中对于不同界面环境下的单空泡演化特性的研究最为广泛．为掌握空蚀产生机理，研究者开展了固体平板壁面环境下的空泡动力学研究，重点关注不同近壁距离下空泡的瞬态溃灭特征，微射流与压力波的产生机制，壁面切应力的演化特性等[6]．出于对水下爆破的研究，自由表面环境下的空泡动力学特性也得到了较多的关注．另外，由于空化在医药领域的应用，研究者开始关注非简单几何环境下的单空泡演化机制[7]．例如Tomita等[8]研究了不同凹凸曲率下的曲面附近空泡溃灭时微射流的动力特性；Tagawa等[9]实验研究了不同角度的角形壁面附近空泡的动力特性；Brujan等[10]通过高速摄影观察了直角壁面附近不同距离下空泡的瞬态行为；Zhang等[11]研究了球形颗粒附近的空泡行为，并分析了不同球形颗粒及空泡尺寸下两者的相互作用；Andrews等[12]实验研究了多孔介质表面特性(孔隙率、孔径和孔形状)对空泡位移及回弹尺寸的影响机制；Sun等[13]通过实验发现固壁表面含气时将对空泡的动力特性产生大的影响；Shi等[14]数值研究了三种不同仿生壁面形状对微射流及空泡溃灭行为的影响；Dualr等[15]通过实验与数值方法研究了当微空泡在与其尺寸相近的固体纤维附近产生时的生长、溃灭等动力学行为．此外，由于空化大多是以多空泡或者空泡团的形式存在，因此多空泡共存环境下的空泡动力学特性也得到了较多关注，但多空泡相互作用的研究大多是在无限域环境下开展．例如Chew等[16]通过实验得到了当两个不同尺寸的空泡共存时，其相互作用所呈现的四种不同的动力学行为；Liang等[17]则通过实验研究了两个相似尺寸的空泡间的相互作用；Luo等[18]实验研究了当双空泡共存时的微射流和压力波演化机制．此外，Yin等[19]通过数值方法研究了无限域内两个空泡的相互作用，重点分析了不同泡间距离下空泡溃灭时的热力学特性；同样地，Liu等[20]数值研究了无限域下双空泡配置下，泡间距离和浮力参数对泡间相互作用及产生的压力载荷的影响，Li等[21]通过数值方法与实验发现了不同泡间距离下空泡融合的三种特征．除以上针对无限域环境下的空泡相互作用的研究之外，Robinson等[22]通过实验和数值研究了当自由液面下存在两个不同尺寸的空泡时，空泡-空泡-自由表面间的相互作用．以上研究大多关注不同界面环境下的单空泡运动特性，或者重点关注无限计算域下当两个空泡共存时的相互影响机制．然而针对近壁环境下的多空泡运动特性研究较少，文献[23-24]仅关注了当近壁区两个或多个空泡横向排列时空泡间的相互作用，文献[25-26]重点分析当多个空泡近壁区溃灭时的瞬态特征．基于此，本研究着重分析当两个空泡尺寸相近时，在不同近壁距离上下排列的配置下，由于近壁面及空泡彼此共存所引起的溃灭瞬态特征的改变及对空泡生命周期的影响．该研究与双空泡横向排列及多空泡近壁溃灭研究相互补充，有助于进一步掌握和理解复杂环境下的空泡相互作用机制及空泡瞬态演化特性．1 数值方法1.1　数值模型本研究采用可压缩流体体积(VOF)方法对空泡运动进行数值模拟．研究介质为水、气两相，且由于空泡膨胀和溃灭过程极快，故与文献[6,27]的假设一致，忽略了两相之间的传质传热．因此，当前数值计算控制方程为∂ρ/∂t+∇∙(ρu)=0;(1)    ∂(ρu)/∂t+∇∙(ρuu)=-∇p+∇∙[μ(∇u+(∇u)T-(2/3)(∇∙u)I)]+Fσ；(2)∂(αiρi)/∂t+∇∙(αiρiu)=0,(3)式中：ρ为密度；t为时间；u和p分别为流速和压力；μ为动力黏度；I为单位张量；Fσ为表面张力，其由连续体表面力模型计算[28]得出；α为气液交界面处的体积分数；下标i表示不同相，i=l表示液体，i=g表示气体．μ和ρ定义为：μ=αgμg+αlμl;(4)ρ=αgρg+αlρl,(5)式中µg和µl分别为气体和水的动力黏度，分别设置为µg=1 mPa∙s和µl=17.1 µPa∙s．此外，在当前计算中假定气液两相均可压缩，根据文献[29-31]，气液两相属性定义如下p=(ρ0c02/nl)[(ρ/ρ0)nl-1]+p0;(6)p/ρk=pref/ρrefk,(7)式中：ρ0为水的密度，ρ0=998.2 kg/m3；c0为声速，c0=1 450 m/s；nl=7.15；p0=2.34 kPa；pref为参考压力，pref=101.325 kPa；ρref为参考密度，ρref=1.225 kg/m3；k=1.4为空气的绝热指数．相应地，水和气体的声速定义为c=(dp/dρ)．(8)根据上述假设和定义，采用ANSYS Fluent软件，基于有限体积法对控制方程进行离散．水相和气相的属性通过Fluent UDF定义．采用文献[29]中的求解方法，基于PISO算法处理压力-速度耦合，采用二阶迎风格式对密度和动量进行离散化，并采用压缩界面捕捉格式对VOF交界面形状计算．瞬态项采用二阶隐式时间格式，采用自适应方法进行时间推进，最大CFL(Courant-Friedrichs-Lewy)数为0.2，初始时间步长设置为1×10-9，每个时间步长迭代50次，以确保在每个计算时间步达到收敛标准(1×10-6)．1.2　计算域和网格细化当前设置的计算域如图1所示，其中计算域顶部及右侧边界均设置为无反射压力出口(恒定值为101.325 kPa)条件，底部边界设为无滑移壁面；同时，考虑到空泡在第一个溃灭周期内具有很好的轴对称特性，将计算域左侧边设为轴对称条件，从而实现空泡运动的三维动力学仿真．此外，如图1(a)所示，初始时初始态为整个计算域充满水，而在靠近壁面的对称轴上设置为两个压力非常高的半圆形气体域，气体域的初始半径及其内压对空泡动力学影响很大，包括影响空泡半径的演变及空泡寿命；因此，在本研究中通过调整空泡的初始半径和内部压力，直至达到期望的空泡大小及其形态演变．10.13245/j.hust.241210.F001图12D计算域与边界坐标系在网格方面，采用了ICEM-CFD软件来创建整个计算域的结构化网格，如图1(b)所示，在两个空泡所在区域内，网格进行细化，而在该区域之外至外边界，网格逐渐稀疏过渡．须要注意的是，本研究的空泡的最大半径大致为400 μm，因此网格细化区域的大小设置为Lx×Ly=5Rmax×11Rmax，其中Ly=11Rmax的原因是上方空泡的圆心到壁面的最远距离为10Rmax．为了保证当前的计算域大小及网格的疏密程度不会对空泡模拟结果产生影响，以单个空泡为研究对象开展了计算域与网格无关性研究．在模拟中对于不同的计算域及网格采用相同的初始条件，即设置相同的空泡初始半径及其内部压力，以空泡无量纲半径(R/Rmax)随时间的演化作为基准，其中R是由空泡体积Ω定义的等效空泡半径，R=(4Ω/(3π))1/3．为研究计算域的收敛性，构建了边长从10 mm(25Rmax)到40 mm(100Rmax)的四个计算域，空泡区域的网格尺寸保持一致，即3 μm×3 μm (0.007 5Rmax)，计算域信息见表1．10.13245/j.hust.241210.T001表1计算域与网格信息计算域计算域尺寸/mm网格数/106网格尺寸/μm110×100.938 403×3220×201.059 503×3330×301.131 523×3440×401.188 073×3图2显示了不同计算域下空泡半径随时间t的变化，可以清楚地看到：随着计算域尺寸的增加，空泡半径演化呈现出收敛的趋势；计算域1预测的空泡半径演化与其他三个计算域的预测结果有较大的差异，而计算域2，3和4的预测结果之间只有微小的差别，甚至计算域3和4给出了完全相同的结果．因此，为了降低计算域大小对于空泡模拟的影响，选取边长为40 mm(100Rmax)的计算域进行后续研究．10.13245/j.hust.241210.F002图2计算域无关性分析随后主要考虑近空泡区域中网格尺寸的大小进行网格无关性研究，网格信息见表2，网格无关性分析如图3所示．网格1得到的空泡振荡周期比其他三个网格的结果都要长，而网格2，3和4预测的结果几乎相同，只有非常微小的差别．因此，选择网格3开展后续研究．须要说明的是，上述研究表明当网格计算域尺寸达到100Rmax，空泡区网格分辨率为0.007 5Rmax时，可满足空泡计算域及网格无关性要求．10.13245/j.hust.241210.T002表2空泡域网格信息网格网格数/105空泡域网格数/105网格尺寸/μm网格15.324 403.525×5 (0.012 5Rmax)网格27.561 255.504×4 (0.010 0Rmax)网格311.880 709.783×3 (0.007 5Rmax)网格423.520 0022.002×2 (0.005 0Rmax)10.13245/j.hust.241210.F003图3网格无关性分析此外，为了验证当前数值方法的准确性，先针对两种不同近壁距离(γ=1.05，γ=0.56，γ=d/Rmax，其中d为初始空泡中心到壁面的距离)下的空泡瞬态演化过程进行了计算．为此，先定义无量纲时间系数τ=t/Tc，其中溃灭时间Tc为微射流到达气泡下表面或气泡达到其最小体积时的时间．图4和图5为实验与数值结果的比较，其中高速摄影中左上角出现的数字对应于一个无量纲的溃灭时间周期，0表示空泡产生，1表示空泡溃灭，更多实验结果详情可参阅文献[32]．可以看到：在这两个近壁距离下，当前计算均很好地捕捉了空泡形状的演变，包括空泡的膨胀、收缩和溃灭过程中的各种瞬时形态，因此可以证明本研究采用的数值方法能够准确还原实验中空泡的瞬态动力学特征．10.13245/j.hust.241210.F004图4实验与数值模拟结果在γ=1.05距离下的比较(Tc=87 μs，Rmax=420 μm)10.13245/j.hust.241210.F005图5实验与数值模拟结果在γ=0.56距离下的比较(Tc=85 μs，Rmax=390 μm)1.3　双空泡配置本文主要研究两个大小相近的空泡在壁面附近的相互作用，旨在分析上方空泡的存在对下方近壁空泡溃灭时瞬态行为的影响．为此在计算域中设置两空泡，下方的空泡称为“空泡1”，上方的空泡称为“空泡2”(如图1所示)．此外，为了使这两个空泡具有相似的最大尺寸，并考虑到本数值方法已被证明能够准确再现最大尺寸约为400 μm的空泡动力学特征，因此本研究主要聚焦于最大半径约为400 μm的空泡相互作用．为此先调整了初始空泡半径及其内部压力，直至不同距离壁面的单个空泡的最大半径与400 μm的最大空泡半径相近．但须要注意的是，由于与壁面的距离不同，导致最大空泡半径不可能完全相同．图6显示了不同近壁距离下空泡达到最大半径时的空泡形状，可以看出空泡最大直径为393.5～412.0 μm，均在400 μm左右．为简便起见，选择400 μm作为空泡特征半径(Rchar=400 μm)，则初始空泡中心的无量纲位置定义为10.13245/j.hust.241210.F006图6不同初始位置的空泡最大半径γβ=dβ/Rchar；(9)Δγ=Δd/Rchar，(10)式中：dβ(β=1，2)为初始空泡中心到壁面的距离；Δd为两个初始空泡中心之间的距离．为了分析不同距离差Δγ下两个近似等尺寸空泡的行为，并比较有无空泡2存在对于空泡1的影响，表3汇总了两个空泡的数值模拟配置，表中以“配置γ1-γ2”的格式为每个配置命名，以便快速索引空泡1的无量纲位置．γ1设置为1～5，这样配置的原因如下：当γ1＜1时，空泡1的最大半径受到明显影响，难以形成与空泡2类似的最大尺寸；而当γ1＞5时，壁面对空泡1的动力学几乎没有影响(这将在后续的讨论中反映出来)．此外，值得注意的是，空泡2的无量纲位置γ2设置为3～10，且最小Δγ设置为2，这样可以保证当两空泡在生长到最大尺寸时，最大半径相似且都在400 μm左右．10.13245/j.hust.241210.T003表3双空泡配置γ1γ23456789101配置1-3配置1-4配置1-5配置1-6配置1-7配置1-8配置1-9配置1-102配置2-4配置2-5配置2-6配置2-7配置2-8配置2-9配置2-103配置3-5配置3-6配置3-7配置3-8配置3-9配置3-104配置4-6配置4-7配置4-8配置4-9配置4-105配置5-7配置5-8配置5-9配置5-102 结果和讨论2.1　不同初始距离下的单个空泡动力学分析为了分析空泡之间的相互作用，须掌握不同近壁距离下单个气泡的溃灭过程．表4显示了溃灭时间Tc时的气泡形状，从表4的溃灭时间和空泡溃灭形态可以看出，越靠近壁面空泡在溃灭过程中受壁面的影响越大．当γ=1时，微射流在87.5 μs时刺穿空泡下表面，且此时空泡仍维持较大的体积；而在2≤γ≤5的范围内，壁面对于空泡的影响明显逐渐减弱，空泡的溃灭时间逐渐减少；但由于壁面的影响依然存在，导致在空泡体积缩小的过程中空泡上下表面压力分布不均，微射流仍得到发展并能够刺穿空泡下表面．当γ≥6时，壁面对空泡第一次溃灭瞬态行为的影响已可忽略不计，此时由于空泡距离壁面较远，溃灭过程中上下表面压力梯度较小，已不能形成稳定的微射流，只有在空穴体积反弹的过程中微射流得以发展．此外，值得注意的是，随着空泡与壁面距离的增大，空泡的溃灭时间单调减少，且当γ=9，10时，溃灭时间的1/2(39.8 μs)已接近于空泡在无限空间内溃灭的理论解TRay=38 μs．通过上述分析发现：只有当空泡与壁面的距离在一定范围内(1≤γ≤5)时，壁面才对空泡第一次溃灭的动力学行为产生较为明显的影响．因此，本研究设置了如表3所示的双空泡配置，近壁空泡距离壁面的范围设置为1≤γ≤5，并且本研究主要关注当在近壁空泡上部存在一个近似尺寸的空泡时，对近壁空泡第一次溃灭的动力学行为的影响．10.13245/j.hust.241210.T004表4不同初始位置空泡在溃灭时刻的形状γ12345678910Tc/μs87.584.982.581.881.180.380.179.879.679.6形态2.2　空泡2对近壁空泡溃灭的影响为了探究空泡2的存在对空泡1溃灭动力特性的影响，本研究主要选取了近壁空泡即将溃灭即空泡1即将被微射流击穿的时刻进行分析．图7显示当γ1=1时，不同空泡间距离下，空泡即将被微射流刺穿时的速度云图；同时，在图7最后显示了当只有单个空泡存在时，空泡在溃灭前期的速度云图，从该图的对比中很容易发现空泡2的存在对空泡1溃灭时微射流的速度存在明显影响．当空泡间距离较小时配置1-4，空泡1溃灭时的微射流速度最大，且随着空泡间距离的增大，空泡1中微射流的速度逐渐减小．但即使空泡间距离加大至配置1-10，空泡1中的微射流速度仍明显大于单个空泡配置下微射流的速度．10.13245/j.hust.241210.F007图7γ1=1时空泡溃灭前微射流云图为了更好地分析微射流速度受空泡距离的影响，图8中显示了该条线上的速度分布，可以看出与单空泡状态下对比，双空泡状态下的近壁空泡微射流速度明显更大，且微射流的最大速度明显随着空泡距离的增大而减小，但减小的趋势并非线性．对于配置1-4，微射流的速度达到160 m/s，大约是单空泡状态下微射流速度的两倍，但是对于配置1-5，该最大速度锐减至122 m/s，此后随着空泡间距离的增大，减小的趋势也逐渐减弱．10.13245/j.hust.241210.F008图8γ1=1时中心轴线速度分布为了进一步分析影响微射流速度的原因，类似于图7和图8，图9和图10分别显示了双空泡及单空泡状态下，空泡溃灭前的压力分布．可以看出：驱动压力随空泡间距离的变化趋势与微射流速度相一致，即随着空泡间距离的增大，驱动压力的最大值逐渐减小．对于双空泡配置1-4及配置1-5，最大压力达到3.05 MPa和2.1 MPa，然后逐渐减小，在配置1-10下，最大压力减小至1.25 MPa，但仍然明显大于单空泡配置下的最大压力(0.82 MPa)．图11~14分析了当γ1=2时，双空泡及单空泡状态下，近壁空泡即将溃灭时的速度与压力分布．从图11和图12中可以看出：与当γ1=1时得到的结论类似，双空泡配置下近壁空泡溃灭前的微射流最大速度明显更大，且该最大速度随着空泡间距离增大而逐渐减小．同样地，从图13和图14中可以看出驱动微射流发展的最大压力也随着空泡间距离的增大而减小．10.13245/j.hust.241210.F009图9γ1=1时空泡溃灭前驱动压力云图10.13245/j.hust.241210.F010图10γ1=1时中心轴线压力分布10.13245/j.hust.241210.F011图11γ1=2时空泡溃灭前微射流云图10.13245/j.hust.241210.F012图12γ1=2时中心轴线速度分布10.13245/j.hust.241210.F013图13γ1=2时空泡溃灭前压力云图10.13245/j.hust.241210.F014图14γ1=2时中心轴线压力分布通过以上对比可以发现：相比于单空泡状态，当在近壁空泡上部存在等尺寸空泡时，近壁空泡溃灭前的微射流速度增大，驱动压力更强，且增幅随空泡间距离的增大而减小．2.3　空泡2对近壁空泡振荡周期的影响在此主要分析空泡2的存在对于近壁空泡振荡周期的影响．为了对比双空泡配置下，不同空泡间距离对近壁空泡振荡周期的影响，图15分别显示当γ1=1，2，3，4，5时，τ=1时刻下不同空泡间距下的双空泡形态．10.13245/j.hust.241210.F015图15不同近壁距离和不同空泡间距下的双空泡形态对比图15与表4中空泡形态可发现，双空泡配置明显改变了空泡溃灭过程中的瞬时形态；从图15中可以看出，空泡1溃灭时的瞬时形态与空泡2及壁面的距离密切相关．当γ1=1时，在配置1-3下，由于壁面及空泡2下表面与空泡1的距离都相近，因此两种壁面均强烈影响着空泡1的溃灭进程，致使其被明显拉伸，上下表面都呈现出扁平状，瞬时形态类似为长方形，这与实验中观察到的现象相同[33]．而随着空泡2与空泡1距离的增大，这种拉伸效应逐渐减小．当γ1=2，3，4，5时也可以发现类似的结论．此外，对比当Δγ=2时的双空泡瞬时形态，即在配置1-3、配置2-4、配置3-5、配置4-6和配置5-7下可以发现：随着空泡1与壁面的距离越来越大，空泡2对空泡1的瞬时形态影响变得更大，以至于在γ1≥2后空泡1的上表面在溃灭瞬态都呈现出扁平状．此外，对比表4和图15，更重要的发现是在双空泡配置下，空泡1和空泡2的振荡周期都被明显延长．对应于表4中γ=1，2，3，4，5状态下的溃灭时刻，近壁空泡均被微射流刺穿，而在图15中同一时刻，可以看到在5种不同距离下的近壁空泡均尚未形成明显的微射流，且当空泡2与空泡1距离更接近时，在同一个时刻空泡1的体积更大，意味着需要更长的溃灭时间．此外，观察不同γ1状态下空泡2的瞬时形态，例如单空泡状态下，当γ=3和γ=4时，空泡分别在t=82.5 μs和t=81.8 μs时刻被微射流刺穿，而在双空泡配置1-3和配置1-4，在t=87.5 μs时刻微射流仍未到达空泡下表面，且当γ=9和γ=10时，空泡在t=79.6 μs时刻达到其第一次溃灭的最小体积，而在双空泡配置3-9和配置3-10下，在t=82.5 μs时刻空泡2仍未达到其最小体积．为了进一步量化双空泡配置下空泡2的存在对空泡1生命周期的影响，统计了不同近壁距离下空泡1被微射流刺穿的时间，即空泡1的生命周期，得到如图16所示的不同空泡间距下空泡生命周期比率，图中：τlife为无量纲的时间系数；Δγ为无量纲的空泡间距离．从图16可以直观地看出：近壁空泡的生命周期随着空泡间距离的增大而不断减小，且减小的趋势越来越弱．此外，比较当γ1=1与γ1=2，3，4，5时的生命周期比率可以发现：在同样的泡间距离Δγ下，当γ1=1时，近壁空泡的生命周期比率明显更大；而当γ1≥2时，上空泡对近壁空泡生命周期的影响规律趋于一致，这意味着尺寸相近的两个空泡在同等的泡间距离下，当空泡1与壁面的距离非常近时，其生命周期不仅与泡间距离相关，而且也受到近壁面的影响．但在空泡1与壁面的距离逐渐变远后，近壁空泡的生命周期则主要与泡间距离相关．10.13245/j.hust.241210.F016图16双空泡不同间距下近壁空泡生命周期比率3 结论本研究基于可压缩的VOF多相流数值方法分析了近壁表面两个尺寸相近的空泡上下排列下的相互作用，重点关注空泡2的存在对近壁空泡1溃灭瞬态特征的影响，得到结论如下．a．通过开展计算域及网格无关性研究，并通过与实验空泡瞬时形态的对比，确定了当前的数值方法能够准确还原近壁区空泡的瞬时特性．b．与单空泡状态相比，双空泡配置下近壁空泡溃灭前的微射流最大速度及溃灭压力明显增大，两空泡的瞬时形态及溃灭时间也由于彼此的存在而发生明显改变．c．在双空泡配置下，近壁空泡溃灭的微射流速度及驱动压力场均随着泡间距离的增大而逐渐缩小，且缩小的趋势也随着泡间距离的增大而逐渐减弱．d．在双空泡配置下，上下两空泡的生命周期均被不同程度延长，且当γ1≥2时，在同样的泡间距离下，上空泡对近壁空泡生命周期的影响规律趋于一致．此外，值得注意的是高压激波的释放也是空泡溃灭的重要特征，关于双空泡配置对激波释放及空泡形态改变的影响机制将在后续的研究中重点分析．