基于伪势模型的表面张力对空化泡溃灭的影响

邹建锋，何小泷，袁 浩，汪凯迪，周昔东

(1.中国电建集团贵阳勘测设计研究院有限公司，贵阳 550081； 2.重庆交通大学 重庆西南水运工程科学研究所，重庆 400074； 3.中国三峡建工(集团)有限公司，成都 610000； 4.重庆交通大学 河海学院，重庆 400074)

1 研究背景

水利枢纽中，高水头和高流速往往伴随着负压的产生，导致空化空蚀破坏，影响水利枢纽运行的安全，因此空化空蚀在水利研究中一直倍受重视。前人研究表明，空化泡外流体与泡内气体压差在短时间内迅速变化，导致空化泡溃灭，溃灭过程中产生的微射流、冲击波、高压、高温是空蚀破坏的主要原因[1-3]。

受限于复杂的边界和有限的观测手段，基于理论分析和试验的空化泡溃灭过程研究多着眼于空化泡形态的演化，难以获得空化泡及周边流体流场的连续变化过程。与之相比，计算流体力学可以直观地反映出空化泡溃灭过程中流场、压力场和温度场的连续变化过程，有助于对空化泡溃灭中丰富的力学机理展开系统研究。传统宏观空化泡溃灭模型主要基于Navier-Stokes(N-S)方程和Rayleigh-Plesset(R-P)方程，此类方法需要引入界面捕捉模型(Volume of fluid，Level-set)对气液交界面进行捕捉，通常需要较小的步长和较密的网格尺度，往往会带来极大的计算量[4-5]。近30 a来，格子玻尔兹曼方法(Lattice Boltzmann Method，LBM)迅速发展为多相流研究的重要方法之一，该方法基于玻尔兹曼动理学理论，求解介观粒子分布函数，统计获得宏观尺度的密度场与流场。与传统宏观数值模拟方法相比，LBM无需求解宏观方程中的对流项和泊松方程，适用于复杂边界，且具有高并行性，已被用于模拟沸腾、蒸发、空化等复杂多相流现象[6-8]。

基于LBM的空化模型主要依托于伪势模型，该类模型利用粒子间的作用势，自动形成气液交界面[8]。Sukop等[9]首先模拟了各向异性条件下空化泡的生长过程，结果与R-P理论解一致，但其液气密度比仅为6.5。Chen等[10]引入精确差分法(exact-difference-method，EDM)外力格式后，将原始伪势空化模型密度比提高到60。Shan等[11]引入多松弛(multi-relaxation-time，MRT)碰撞算子，进一步将LBM伪势空化模型密度比提高到750，液气密度比接近于真实水与水蒸气密度比。随着热力学方程被引入，LBM空化模型对空化现象的研究进一步拓展到空化溃灭时温度演化规律的研究[12-13]。

近壁区空化泡溃灭过程中，常伴随着气液交界面的剧烈变形，局部气液交界面具有极大的曲率，在这些地方，表面张力往往对空化泡演化形态、溃灭强度、微射流分布具有重要的影响。而目前基于LBM伪势模型的空化泡溃灭过程的模拟极少考虑表面张力对空化泡演化的影响。因此本研究引入表面张力调节外力项，建立可调节表面张力的LBM伪势空化模型，分析不同表面张力对空化泡溃灭过程的影响，揭示表面张力对空化泡微射流流速分布、最大压力变化、空化泡形态演化的影响。

2 数值模拟方法

2.1 LBM伪势模型

本研究中采用多松弛系数碰撞算子、Li修正外力格式[14]、Carnahan-Starling(C-S)状态方程等手段获得大密度比、大黏滞系数比条件下，数值稳定性更好的LBM伪势模型，其对应粒子分布方程为[14]

(1)

(2)

式中：ωi为平衡态函数对应权重；ρ为密度；ueq为平衡态速度；cs为格子声速。对于D2Q9格子模型，ω0=4/9，ω1～4=1/9，ω5～8=1/36。离散速度ei为

(3)

式中：c为格子常数，c=Δx/Δt；松弛矩阵Λ由不同松弛参数构成，即

(4)

(5)

宏观密度ρ和流场实际速度u可通过以下计算方式求得：

(6)

外力项Si采用Li等[14]提出的改进外力格式，满足大密度比条件下热力学一致性，即

meq=Mfeq=

(7)

式中：ε为调节热力学一致性的参数；ψ为粒子间的伪势；τe为能量松驰系数；F为作用在粒子上的合力；Fm为粒子间的相互作用力，可以通过伪势强度进行求解，即

(8)

式中：G为粒子间的相互作用强度；w为不同方向的权重，w0=0，w1～4=1/3，w5～8=1/12[15]。引入非理想流体状态方程后，可获得求解伪势ψ[16]的表达式为

(9)

式中P为压力，本研究可用C-S状态方程求解[16]，即

其中：

式中：Tc为临界温度；Pc为临界压力。

本模型采用线性插值获得不同密度条件下的松弛系数，以获得不同黏滞系数比[8]，即

(11)

式中：τυ(x)为混合流体黏滞系数；τg为气相黏滞系数；τl为液相黏滞系数；ρg为初始气相密度；ρl为初始液相密度。

为探究表面张力对空化泡溃灭过程的影响，本模型引入源项C，获得可单独调节表面张力的伪势模型[17]，即

(12)

式中Q可以通过下列方程求解，即

(13)

式中κ为表面张力调节参数。

2.2 边界条件及初始参数

近壁区空化泡溃灭计算示意图如图1所示，计算区域为200 lu×200 lu矩形区域，r0为空化泡初始半径，d为空化泡中心与壁面之间的间距，λ=d/r0为空化泡与壁面之间的无量纲距离。计算域顶部边界为Zou-He压力边界[8]，底部边界为无滑移边界，左右边界为非平衡外推边界。

图1 近壁区空化泡溃灭计算域示意图Fig.1 Computational domain for cavitation bubble collapse in near-wall region

为促使模拟过程尽快达到平衡态，不影响空化泡演化过程，空化泡及周边计算域流场按照式(14)进行初始化。

(14)

式中：(x0，y0)为液滴中心坐标；w为界面厚度，本研究中选取界面厚度w=5 lu。

3 模型验证

本研究采用Laplace定律和无限域空化泡溃灭过程验证空化模型的可靠性。

3.1 Laplace定律

首先利用Laplace定律验证并获得不同κ值条件下的表面张力。Laplace定义为，当空化泡或液滴达到平衡状态时，空化泡内外压差与半径的倒数成正比，即其斜率即为该条件下的表面张力，Δp=σ/r。式中Δp为空化泡内外压差，σ为表面张力，req为空化泡平衡态条件下半径。

本文选取5种初始半径用于Laplace定律的研究，分别为r0=35、40、45、50、55 lu。不同值条件下空化泡内外压差Δp和平衡态半径r之间关系如图2所示。本模型获得的空化泡内外压差与其平衡时刻的半径的倒数具有较好的线性关系。随着κ值的减小，表面张力由0.001 9增大到0.018 3，表明当前模型可以有效调节表面张力大小。

图2 不同κ值下空化泡内外压差和平衡态半径的关系Fig.2 Relations of pressure difference inside and outside of the bubble against bubble radius under different κ values

3.2 无限域空化泡溃灭过程

为进一步验证模型的准确性，本研究还对无限域内空化泡溃灭过程进行模拟，并将LBM模拟结果与Rayleigh-Plesset方程数值解进行对比验证。考虑表面张力，对应的Rayleigh-plesset空化泡演化方程为[18]

(15)

模拟中设置计算域四周均为Zou-He压力边界，初始液相压力p∞=0.006 79 mu/(lu·tu2)，空化泡内初始压力pv=7×10-7mu/(lu·tu2)，空化泡半径r=50 lu，表面张力σ=0.010 2 mu/tu2，其他参数与2.2节一致。无限域内空化泡溃灭中半径演化过程的LBM模拟结果与Rayleigh-Plesset方程对比如图3所示，图中r*为无量纲空化泡半径，t*为无量纲空化泡溃灭时间，LBM模拟结果与R-P方程理论解一致，尽管在空化泡溃灭阶段初期存在部分误差，其可能是由于LBM模拟中对初始化流场调整而引入的。

图3 R-P方程理论解与LBM模拟结果对比Fig.3 Comparison between the result of Rayleigh-Plesset equation and LBM simulation

4 表面张力对空化过程的影响

本节选取空化泡内外初始压差Δp=0.003 8、0.006 8、0.010 2 mu/(lu·tu2)三种压力，无量纲距离λ=1.6条件下，不同表面张力系数对空化泡溃灭过程中溃灭形态、流场、最大微射流速度和最大溃灭压力的影响。

图4展示了Δp=0.006 8 mu/(lu·tu2)不同表面张力条件下空化泡溃灭最终形态及流场。对比不同表面张力条件下压力分布，随着表面张力减小，空化泡上方高压区面积增大，同时空化泡与壁面之间的低压区面积减小。在较小表面张力条件下，空化泡更易发生变形，界面曲率半径更小，空化泡溃灭最终时刻形变更大，导致溃灭形成的微射流流速更为集中，这与Zhang等[19]的研究结论一致。

图4 不同表面张力条件下空化泡溃灭最终时刻密度、速度和压力分布密度Fig.4 Distribution of density，velocity，and pressure at the final collapse stage under different surface tensions

图5展示了Δp=0.006 8 mu/(lu·tu2)时不同表面张力条件下空化泡溃灭过程中最大流速和最大压力演化过程。溃灭过程中更大的表面张力条件下空化泡积蓄了更多的表面能，并在溃灭最后阶段迅速释放，导致了更大的空化泡溃灭速度。随着表面张力增加，空化泡溃灭时间增加，同时空化泡溃灭最大流速也相应增加。当表面张力由0.001 9 mu/tu2增加到0.018 3 mu/tu2时，空化泡溃灭流速由0.404 lu/tu增加到0.512 lu/tu，速度增加了26.8%。由图5(b)可知，在空化泡溃灭最后阶段，由于空化泡内水蒸气在短时间内发生相变，溃灭压力迅速增加。且更大的表面张力阻碍了空化泡的形变和溃灭，导致空化泡溃灭压力减小。在Δp=0.006 8 mu/(lu·tu2)压力条件下，当表面张力由0.001 9 mu/tu2增加到0.018 3 mu/tu2时，空化泡溃灭过程中最大压力由0.049 mu/(lu·tu2)增加到0.064 mu/(lu·tu2)。

图5 不同表面张力条件下空化泡溃灭过程中最大速度和最大压力演化过程Fig.5 Evolution of maximum velocity and maximum pressure in the collapse process under different surface tensions

不同表面张力和初始空化泡内外压差条件下空化泡最大微射流流速、压力和溃灭时间变化规律如图6所示。随着初始内外压差增大，不同表面张力条件下空化泡溃灭产生的最大微射流流速、最大溃灭压力均随之增大，但溃灭时间则随之减小。初始空化泡内外压差为0.003 8 mu/(lu·tu2)时，当表面张力由0.001 9 mu/tu2增加到0.018 3 mu/tu2，最大微射流流速增加了55%，最大溃灭压力增加了74%。初始空化泡内外压差增加到0.010 2 mu/(lu·tu2)时，空化泡溃灭最大微射流流速仅增加了22%，溃灭压力仅增加了6%，最大微射流流速和最大溃灭压力增加幅度均随着空化泡初始内外压差的增加而减小，说明空化泡初始内外压差的增加会减小表面张力对溃灭过程的影响。

图6 不同压力条件下空化泡溃灭过程中各参数随表面张力的变化规律Fig.6 Vaiation law of parameters with surface tension under different pressures during the process of civitation bubble collapse

而在相同压力条件下，空化泡溃灭速度、溃灭压力和溃灭时间均随着表面张力的增加而增加。根据Bjerknes力的定义[20]，有

(16)

空化泡溃灭时间随着表面张力的减小而减小，而时间质量变化率则随之增大，导致Bjerknes力增大，说明较小的表面张力会促使空化泡朝向壁面发生溃灭。

依托于试验研究和宏观数值模拟方法，研究者们针对表面张力对空化的影响开展了系统性研究[21-22]。LBM模拟结果表明表面张力对空化泡溃灭强度具有重要影响，其Bjerknes力随表面张力减小而增大，但溃灭强度却随之减小，模拟结果定性上与前人试验结果和数值模拟结果一致。掺气减蚀是有效减小空化空蚀的重要手段，研究表明掺气浓度将减小流体表面张力，进而减小空蚀强度。

5 结 论

本文基于可调节表面张力的LBM伪势空化模型，研究了表面张力对近壁区空化泡溃灭的影响，得到结论如下：

(1)表面张力直接影响空化泡附近区域压力演化，并对空化强度产生影响。较小的表面张力在增加空化泡上部高压范围，同时减小壁面与空化泡之间的低压区范围。空化泡水汽交界面抗形变能力随之减弱，导致微射流范围更集中。

(2)空化泡溃灭时产生最大微射流和溃灭压力均随表面张力增加而增加。较大的表面张力条件下空化泡演化过程中积蓄了更多的表面能，并在溃灭阶段释放，加剧溃灭强度。

(3)空化泡溃灭时间随表面张力减小而减短，导致Bjerknes力增强，加剧了空化泡朝壁面溃灭的趋势。