液体内局部压强低于饱和蒸气压以下时，流体中会出现空化现象[1]．空化通常会经历气体空泡生成、发展、溃灭等过程．当空泡发生溃灭时，周围的流体会迅速向空泡内涌入，空泡内部及周围流场中压强和温度会快速升高[2]，若空泡溃灭在固体壁面附近，则会产生冲击波与微射流[3]，对固体壁面造成一定损害．高速运行的水下航行器表面会产生空泡．伴随着空泡的脱落、溃灭，航行器周围的流场的稳定被破坏，表面压力波动增加[4]，不利于航行器在水下平稳运行．对于空泡溃灭问题的研究，不仅可以降低由于空化效应带来的噪声、侵蚀及航行器表面压力波动大等不利影响，也可为空化效应更有效地被应用于生物、医学等领域提供参考依据．除理论与实验研究外，数值模拟也是研究空泡溃灭问题的重要手段．常见的数值模拟方法，如有限体积法[5]、有限元法[6]、边界元方法[7]等，往往须要求解压力泊松方程同时结合流体体积法[8]或水平集法[9]等界面追踪方法，增加额外的计算量．与前述几种传统的数值模拟方法不同，格子Boltzmann方法(lattice Boltzmann method，LBM) 通过介观分布函数的碰撞和迁移过程模拟流场时空演化过程[10]．作为一种新兴的数值方法，LBM具有物理意义清晰、边界处理简单、并行度高等优势，已在多相流模拟中[11-12]取得了大量的研究成果．LBM框架下有多种多相流模型，常见的有相场模型[13]、自由能模型[14]、颜色梯度模型[15]、伪势模型[16]等．其中基于匀场理论的伪势模型无须求解相界面追踪方程，可自动追踪相界面的运动并实现相的分离，计算效率高，广泛应用于多相流问题的模拟中．利用LBM伪势模型进行空泡溃灭模拟方面，文献[17]基于文献[18]改进伪势力模型研究了曲面边界处空泡溃灭过程．文献[19]研究了壁面疏水性对于近壁面单空泡溃灭的影响规律．文献[20]对壁面温度对于空泡溃灭的影响及其中表现出的热动理学行为进行了研究．然而，目前基于LBM的空泡溃灭问题的研究，多集中于二维情况，对于空泡溃灭的时间与空泡初始半径之间的定量关系、温度对于空泡溃灭过程影响规律的研究相对较少．基于此，本研究以LBM为核心求解器，利用Shan-Chen模型与C-S状态方程，开展定压边界三维空泡溃灭问题的模拟，通过对溃灭过程中半径与溃灭时间定量研究，获得了不同初始半径及温度条件下三维空泡的溃灭规律．1 数值方法1.1　格子Boltzmann方法本研究使用格子Boltzmann方法作为流场的核心求解器．与传统的模拟方法不同的是，格子Boltzmann方法是一种介观模拟方法，利用介观的分布函数的碰撞、迁移过程实现流场的演化，通过分布函数的线性化组合可以获取流场的宏观物理量，带外力项的分布函数演化方程为fi(x+eiΔt, t+Δt)=fi(x, t)-fi(x, t)-fieq(x, t)τ+Fi(x, t)Δt ,式中：fi(x,t)为t时刻，x处分布函数的第i个分量；Δt为最小时间步长；Fi(x,t)为外力项；fieq(x,t)为平衡态分布函数．又有fieq(x, t)=ωiρ1+ei∙ucs2+(ei∙u)22cs4-u22cs2，式中：ωi为权重系数、ω0=1/3，ω1～ω6均为1/18，ω7～ω18均为1/36；ρ为格点处的密度，另设参考密度ρ0取值参考[21]；ei为格子离散速度；u为格点处的速度；cs为格子声速．c=Δx/Δt为格子速度，且c与cs之间满足关系cs=c/3，τ为弛豫时间与运动黏性系数间满足v=cs2(τ-0.5)Δt．对格点处的分布函数作线性计算，可以获取：ρ=∑i=018fi；ρu=∑i=018fiei．1.2　伪势模型伪势模型最早是由文献[22]提出，该模型中引入粒子间相互作用力F(x)，其表达式为F(x)=-Gψ(x)∑i=08ωiψ(x+eiδt)ei，(1)式中：G为粒子间相互作用常数；ψ(x)为势函数，ψ(x)=2(Peos-ρcs2)/(Gcs2)，Peos为状态方程．有学者就状态方程的选取进行了研究[23]，结果表明C-S状态方程具有计算精度高、计算温度范围广、数值稳定性好、计算密度比大等优势．C-S状态方程为     Peos=ρRT∙1+bρ4+bρ42-bρ43 /1-bρ43-aρ2, (2)式中：a=0.496 3R2Tc2/pc；b=0.187 27RTc/pc；RT=1/3；Tc为参考温度；pc为参考压强．由式(1) 和(2)可以获取流场中粒子相互作用力．添加的力的构造形式为EDM格式[24]，力的作用效果体现在对于分布函数的修正上，有：Fi(x, t)=fieq(ρ, u+Δu)-fieq(ρ, u)；Δu=FΔt/ρ．此时流体的宏观速度定义为v=u+FΔt/(2ρ)．计算中C-S状态方程的参数选取为a=1，b=4，R=1，可得：Tc=0.094 33；参考密度ρc≈0.521 8/b，文中近似取值为0.130 45．流场中密度的初始化取值可以通过麦克斯韦等面积法曲线获取．选取不同的温度值T可以实现不同的密度比的模拟，T/Tc的值一般选取在0.55~1．当T/Tc=1时，ρ/ρc=1，此时气液相密度相同，不发生相分离的现象，当T/Tc越接近0.5时，两相密度比就越大，但计算的数值稳定性差，并且当温度比T/Tc越低时，产生的非物理的伪势速度会越大，因此为了同时兼顾密度比和计算的稳定性，如无特殊说明，T/Tc取值为0.6．流场中密度初始化采取如下方式，     ρ(x, y)=ρ1+ρv2+ρ1-ρv2×tanh2(x-x0)2+(y-y0)2+(z-z0)-R0W,式中：ρl为液相密度；ρv为气相密度；双曲正切函数tanh(x)=(ex-e-x)/(ex+e-x)；(x0，y0，z0)为空泡初始中心点的位置；R0为空泡初始半径；W为气液界面厚度，选取W=5．如果没有特殊说明，长度均采用Δx(Δx为最小网格尺寸)，进行无量纲化处理，时间采用Δt进行无量纲化处理，密度采用ρ0进行无量纲化处理，压力采用p0进行无量纲化处理，其中p0=ρ0(Δx/Δt)2．2 计算结果与分析2.1　算例验证为验证计算结果的准确性，首先在流场中设置两个立方形的空泡，由于表面张力的作用，因此立方形的空泡会逐渐向可以稳定存在的球形空泡转化，并且当空泡距离比较近时会出现双空泡融合的现象．图1展示的是在流场中设定了两个方形空泡随时间演化的过程，其中网格分辨率为723，初始流场各点速度为零，初始的分布函数取平衡态分布函数，边界条件使用的是周期性边界条件．在演化过程中，空泡会逐渐变成球形，并且随着时间的推进，双空泡会融合成一个大的空泡在流场中稳定存在．方形空泡的演化过程与文献[25]的结果相似，不同之处在于文献[25]模拟的是单空泡的演化过程，而这里模拟了双方形空泡的演化过程，并且由于两个空泡距离近而发生了空泡融合的现象．10.13245/j.hust.240630.F001图1双空泡演化与融合过程此外，以相分离过程为例，对于给定的随机密度场，在不同的初始条件下，会生成稳定的不同形状的两相界面．图2展示的是不同条件下的相分离的结果，流场初始密度选取为ρ+init(0.001-N/104)，N为区间[0，25]获得的随机整数，ρinit为流场初始的平均密度，计算时候选取ρinit为0.15和0.25，T/Tc取值为0.6和0.7，网格分辨率为723，边界处使用的是周期性边界方法．图2中：t*=t/Δt；ρ*=ρ/ρ0；ρinit*=ρinit/ρ0．三维的仿真结果生成的形状如图2(a)为圆柱状，图2(b)球形及图2(c)为方形层状．其中，图2(b)计算结果与文献[26]使用相场模型计算的结果相似．10.13245/j.hust.240630.F002图2不同温度和初始密度的相分离结果拉普拉斯定律是多相流的经典基准测试用例．该定律可以表述为静态气泡内外压差Δp*与空泡半径的倒数1/R*成线性关系，其中Δp*=Δp/p0，R*=R/Δx．为了更加定量地说明本文算法的准确性，这里使用不同半径的空泡内外压差与平衡状态下空泡半径之间的关系用于验证拉普拉斯定律，共选取了5组计算结果如图3所示，图中横纵坐标变量的单位均为格子单位．从计算结果中可以看出，LBM模拟出的结果与拉普拉斯定律相符．10.13245/j.hust.240630.F003图3空泡内外压差与半径倒数之间关系2.2　温度对空泡溃灭影响规律研究为研究温度对三维空泡溃灭过程的影响规律，这里系统模拟了边界压强和空泡初始半径恒定条件下，不同温度下空泡的溃灭过程．首先，对流场的计算域作网格无关性检验，在固定边界压力为0.001，空泡初始半径为15的情况下，统计了4组不同的网格分辨率下的空泡溃灭时间，如表1所示．其中相对误差为不同网格分辨下获取的空泡溃灭时间与最高网格分辨率下的空泡溃灭时间的相对误差．当空泡半径小于一个网格尺度时，已经无法被捕获，这里采用外插的方法得到了空泡溃灭的时间．从表中可以看出，当网格分辨率达到1003时，相对误差可以低于2%，因此本研究选用的计算域网格分辨率为1003．10.13245/j.hust.240630.T001表1不同网格分辨率下网格无关性检验50364310031283空泡溃灭时间156.2162164.1166.9相对误差/%6.42.91.70.0图4为当边界压强为0.001、初始半径R*=30、T/Tc=0.6时空泡溃灭过程中的流场密度等值面．其值为气液相密度的平均值，该等值面即为流场中空泡边界面．10.13245/j.hust.240630.F004图4T/Tc=0.6时空泡溃灭过程中密度等值面图图5(a)计算了无量纲空泡初始半径R0*为15，20，25，30，35的情况，其中R0*=R/Δx．可以看出：随着初始空泡半径的增大，空泡溃灭的时间也会逐渐增大，空泡溃灭过程中，溃灭速度|dR/dt|伴随着溃灭过程逐渐增大．观察溃灭总时间和初始半径之间的近似线性关系，对于所有的计算结果做线性拟合．图5(b)为温度T/Tc=0.6时空泡溃灭时间与空泡初始半径之间的变化关系．从图中可以发现:在恒定温度下，空泡初始半径越大，则空泡溃灭所需时间就越长，且二者服从线性关系t*∝14.68R*．10.13245/j.hust.240630.F005图5T/Tc=0.6时空泡半径与时间的关系使用空化初始半径R0和空泡溃灭总时间tco对溃灭过程的空泡半径和时间进行无量纲化处理，即R1*=R/R0，t1*=t/tco，如图6所示．从图中可以看出：相同条件下，初始半径不同的空泡具有相似的溃灭过程，溃灭速度随着时间不断增加．10.13245/j.hust.240630.F006图6无量纲后溃灭空泡半径随时间变化为进一步揭示温度对空泡溃灭过程的影响规律，这里又系统模拟了边界压强恒定条件下，不同温度时空泡的溃灭过程．图7(a)中T/Tc=0.65，计算的空泡初始半径R0*分别为15，20，25，30，35，可以看出：随着初始空泡半径的增大，空泡溃灭的时间也会逐渐增长，空泡溃灭过程中，溃灭速度伴随着溃灭过程逐渐增大．观察溃灭总时间和初始半径之间具有近似线性的关系，对于所有的计算结果作线性拟合．图7(b)显示的是温度T/Tc=0.65时空泡总溃灭时间与空泡初始半径之间的变化关系．从图中可以发现：在恒定温度下，空泡初始半径越大，则空泡溃灭所需时间就越长，且二者成线性关系t*∝13.68R*．10.13245/j.hust.240630.F007图7T/Tc=0.65时空泡半径与时间的关系图8中T/Tc=0.7，图8(a)空泡初始半径R0*分别为15，20，25，30，35，可以看出：随着初始空泡半径的增大，空泡溃灭的时间也会逐渐增长；空泡溃灭过程中，溃灭速度伴随着溃灭过程逐渐增大．图8(b)中所示为温度T/Tc=0.7时空泡溃灭时间与空泡初始半径之间的变化关系．从图中可以发现：在恒定温度下，空泡初始半径越大，则空泡溃灭所需时间就越长，且二者为线性关系t*∝15.12R*．10.13245/j.hust.240630.F008图8T/Tc=0.7时，空泡半径与时间的关系为了直观比较，不同温度和初始半径得到的三维空泡溃灭时间汇聚如图9所示．图9给出的是不同温度、空泡溃灭的时间与空泡初始半径的关系图．直线斜率[27]取自如下计算结果，即     tc=3ρ2P0∫0RmaxRmax3R3-1-(1/2)dR≈0.915Rmaxρ/P0. (3)10.13245/j.hust.240630.F009图9温度对三维空泡溃灭影响规律从图9中可以看出：不同温度情况下，空泡溃灭时间与初始半径均呈现线性关系，与式(3)计算结果是相符的，且斜率随着温度升高略有增大．图10为不同温度、初始半径下空泡溃灭速率随着时间的变化规律，图中：dR*=dR/Δx；dt*=dt/Δt．从整体趋势上来看，不同温度与初始半径的空泡溃灭速率都是随着时间不断增加．其中，相同温度条件下，空泡初始半径越大，空泡溃灭速率最大值出现的时间越晚，且该时间与空泡溃灭的时间是相近的．对于同一半径而言，T/Tc越大，空泡溃灭速率出现的时间越晚．不同温度及初始半径空泡的最大溃灭速率均在1附近．10.13245/j.hust.240630.F010图 10不同温度、初始半径对空泡溃灭速度的影响规律1—T/Tc=0.60，R=15；2—T/Tc=0.60，R=25；3—T/Tc=0.60，R=35；4—T/Tc=0.65，R=15；5—T/Tc=0.65，R=25；6—T/Tc=0.65，R=35；7—T/Tc=0.70，R=15；8—T/Tc=0.70，R=25；9—T/Tc=0.70，R=35．3 结语相同的温度条件下，空泡溃灭的总时间与空泡的初始半径之间呈线性相关关系，空泡溃灭的速度随着时间不断增加，且空泡溃灭的过程具有相似性；不同的温度条件下，空泡半径与空泡溃灭时间线性关系仍存在，但直线的斜率受到温度影响会有所改变．