曲面固壁附近空化泡溃灭动力学的流体体积数值研究*

单鸣雷 杨 雨 邱顺俐 束方勇 韩庆邦

(1 河海大学常州市传感网与环境感知重点实验室并江苏省输配电装备技术重点实验室 常州 213022)

(2 江苏省“世界水谷”与水生态文明协同创新中心 南京 211100)

0 引言

空化是指液体中局部压强降低，其内部发生汽化并形成微小空穴(空化泡)的过程。流体中的空化泡在遇到周围液体压强增大的情况下，体积将急剧缩小或溃灭。空化泡溃灭时间很短(微秒级)，溃灭时会在局部产生极高的瞬时压强。当溃灭发生在固体表面附近时，流体中不断溃灭的空化泡所产生的极高压强反复作用于固体表面，对固壁面造成破坏，这种现象称为空蚀，也称气蚀[1-2]。一方面，空化效应会产生如水轮机的叶片损伤[3]、高速轴承表面材料呈点状脱落及螺旋桨的空蚀噪声[4-5]等不利影响。另一方面，人们发现合理利用空化效应会给生产生活带来极大的便利，如空化清洗技术用于裂缝表面的清洁[6-8]、使用空化泡作为药物的运输工具从而将药物送达病变组织表面[9-10]。因此，无论是避免还是利用空化效应，都有必要对固壁附近空化泡溃灭过程进行深入研究。

目前空化泡在固壁附近溃灭的动力学研究已相当丰富，模型结构也从原始的平直固壁附近空化泡溃灭模型进一步发展为任意几何形状固壁附近空化泡溃灭模型。Tomita 等[11]从理论和实验的角度给出了曲面边界影响下的气泡动力学特性，研究表明，如果边界稍微凹陷，最显著的迁移发生在第一次气泡溃灭的时候。液体射流冲击气泡表面远侧的速度随着边界曲率的增加而增加。Ma等[12]利用霍普金森杆(the Hopkinson bar,HPB)测量技术，对不同曲率边界和不同爆炸距离下的电火花气泡动力学和载荷特性进行了一系列实验研究，并通过实验结果总结了一些物理量的变化规律。Brujan等[13]对超声诱导产生的空化泡进行了实验研究，比较了两种尺寸的空化泡近壁溃灭时的射流速度与冲击压力的大小，发现尺寸较小的空化泡其射流速度与冲击压力都要小于大尺寸空化泡。Li 等[14]利用流体体积(Volume of fluid,VOF)法研究了空化泡在锥形固壁附近溃灭的过程，提出了锥形固壁角与空化泡溃灭形态、溃灭时间以及液体射流的关联，计算了近锥形固壁空化泡溃灭时间的理论值，并通过数值模拟验证了此理论值的有效性。魏梦举等[15]基于合理的简化模型，重点分析了空泡溃灭过程对球体的影响，研究表明，百微米级气泡能够通过射流方式驱动同尺度固体颗粒向前显著运动，而气泡受限的程度对射流驱动的效果具有重要影响，颗粒半径越小、气泡与颗粒间距越近，空泡溃灭时对颗粒的推动效果越明显。

以上研究都是关于非平直壁面附近或其他复杂边界附近空化泡溃灭动力学问题的相关研究工作，这些工作并没有涉及空化泡溃灭过程中的热力学特性的研究。然而在空化泡溃灭的最后阶段，热效应发挥了极其重要的作用，随着空化泡的溃灭收缩，泡内温度急剧上升，溃灭瞬间温度扩散到周围流体中，对泡外流体或固壁面都将产生重要影响。为更全面地了解空化泡溃灭过程，近年来，专家学者们对空化泡溃灭阶段的热力学问题开展了一系列实验。Vokurka 等[16]对电火花诱导生成的空化泡热力学特性进行了实验研究，空化泡表面的温度通过光学传感器和声传感器进行测量，实验结果显示，空化泡表面温度可以达到5900 K。但实验存在过多的不可控因素，且空化泡溃灭过程较为复杂，其温度测量存在较大的难度。因此，采用数值方法[17-21]研究空化泡热力学问题显得十分必要。Yang等[22]采用格子玻尔兹曼方法研究了空化泡溃灭的热力学效应，描述了气液相的温度演化规律。Liu 等[23-24]使用VOF 模型讨论了单壁与平行双壁附近空化泡溃灭对固壁表面的热影响，研究结果表明，在空化泡被射流穿透时，射流速度达到最大值，由于射流冲击，壁面温度明显降低。但以上研究的物理模型均为平直固壁，未探讨其他固壁模型对空化泡溃灭的影响。

本文采用VOF 数值模拟方法研究曲面固壁附近空化泡溃灭的过程，并对此模型下空化泡形态演化进行实验观测，验证数值模拟方法的有效性。首先，分析了空化泡近壁溃灭过程中相场、压力场、温度场的演化，讨论了空化泡近壁溃灭阶段的热动力学作用机理。其次，研究了曲面固壁尺寸、位置参数对空化泡溃灭时间、射流速度的影响。最后，详细探讨了不同条件下，曲面固壁温度的变化情况，深入分析了高温损伤对固壁表面的作用。

1 数值方法

数值计算固壁附近空化泡运动过程时，假设空化泡内为水蒸气，满足理想气体状态方程。整个空化泡溃灭过程不考虑重力因素的影响[25]，忽略气体与液体之间的质量传递。存在相界面的流动中，直接求解流场中的物理量会使其在界面处发生间断，导致计算的发散或计算结果的虚假，为保证界面处所要求解物理量值的合理性，对这类流动进行数值研究时，需要进行足够精确的界面捕捉。VOF方法可以较好地保持流体质量守恒并且可以处理界面拓扑结构的变化[26]，因此本文采用VOF 方法来计算空化泡运动过程。

1.1 控制方程

空化泡流场中的流动方程由连续性方程、动量方程及能量方程组成，表达式分别为[27]

其中，ρ表示流体的密度，μ表示流体的黏度，U表示流场速度，p表示流场压力，σ表示表面张力系数，κ表示表面曲率，n表示两相流界面中的单位法向量，T表示温度，E表示总能量，k表示流体的传热系数。

VOF 方法通过一个标量的输运方程来描述各相的体积分数，表达式为[28]

其中，α表示构成相的体积分数大小，当α= 0 代表该网格处于空化泡内部，当0＜α ＜1 代表该网格处于空化泡外即液体内部，当α= 1 代表该网格处于两相交界面上。

各相的密度以及黏度可分别表示为

其中，α1+α2= 1，α1代表液体的体积分数，α2代表气体的体积分数。

1.2 计算模型

本文采用VOF 方法研究空化泡邻近曲面固壁溃灭的过程，几何模型如图1 所示。曲面固壁为球形，半径设为r，球形空化泡的初始半径设为R0。为方便描述空化泡与固壁间的距离，引入无量纲位置参数λ=d/R0，其中d为空化泡中心到曲面固壁右侧一点的距离，在数值模拟过程中，改变曲面固壁的半径时始终保持曲面固壁右侧一点固定。为减少压力边界条件对溃灭过程的影响，将计算域大小设为30R0×30R0。计算域的四周均为压力出口边界，曲面固壁表面为无滑移固壁边界条件。

图1 计算模型Fig.1 Computational model

在数值模拟软件中导入网格文件，由于Gambit 中构建模型结构时默认的长度单位为m，因此要根据实际目标需求对网格的尺寸进行重新定义并检查网格质量。采用基于压力的瞬态求解器，选用多相流模型中的VOF 模型进行计算，开启能量方程。压力-速度耦合方式运用PISO 算法，空间离散方法中，梯度项使用Least Squares Cell Based 格式，压力项使用PRESTO!格式，体积分数项使用Geo-Reconstruct格式，密度项、能量项、动量项均采用Second Order Upwind 格式[29]。计算的时间步长设置为10-8s，开启残差曲线监控即可开始迭代计算，观察计算结果是否收敛。

1.3 网格无关性验证

不同的网格数量会产生不同的计算结果，为使计算结果误差较小，需要选取恰当的网格数。在比较不同网格数的计算结果时，采用如下初始计算条件：曲面固壁半径r= 0.2 mm，空化泡半径R0=0.2 mm，无量纲位置参数λ=1.1，液体中压强pl=3×106Pa，空化泡内压强pv=1.01×105Pa。

图2 为空化泡横轴上射流速度的变化情况，其中横轴表示空化泡轴上的位置坐标，纵轴表示轴上的射流速度。观察发现射流速度在某一位置达到最大值，且不同网格数下的空化泡轴上射流速度差异较小，几乎可以忽略不计。综合考虑计算时间、计算资源等的限制，最终选取网格数为200000的工况进行数值模拟。

图2 网格无关性验证Fig.2 Verification of the grid independence

1.4 实验验证

为验证数值模拟结果的有效性，本章节将曲面固壁附近空化泡演化的实验研究与数值仿真结果进行对比，如图3 所示。图3(a)为高压放电致泡实验下空化泡近曲面固壁溃灭的形态变化图，图3(b)为数值模拟仿真得到的空化泡近曲面固壁溃灭的相图。

图3 空化泡形态演化实验结果图和仿真结果图Fig.3 Comparison of cavitation bubble profiles between the experimental results and the simulation results

将实验结果与仿真结果进行对比，由于实验研究和数值仿真中空化泡尺寸的差异，其溃灭时间不能直接比较，需要进行归一化处理。比较两图可以发现空化泡的形态演化过程一致，且归一化后得到的无量纲时间相同。空化泡形态演化均为从初始的球状逐渐被压缩成类似椭球形状，可以明显地看出空化泡上下泡壁收缩速度大于左右泡壁收缩速度，左侧泡壁收缩的速度与右侧泡壁收缩的速度相比较慢。因此，在t*=0.5789时刻，左侧泡壁出现尖端现象，随后继续溃灭。在这过程中空化泡没有出现类似于近平直固壁溃灭时产生的凹陷或者被射流穿透现象，而是维持一个整体溃灭直至空化泡表面完全消失。由于空化泡溃灭形态演变的实验结果与数值仿真结果的吻合性较好，故进一步验证了数值模拟曲面固壁附近空化泡溃灭过程的正确性，为接下来对数值模拟结果的分析奠定基础。

2 结果分析

2.1 压力场、速度场、温度场演化分析

为了进一步探究空化泡邻近曲面固壁溃灭机理，对曲面固壁附近空化泡溃灭过程中压力场分布情况进行了分析。图4为曲面固壁半径r=0.3 mm，空化泡初始半径R0=0.2 mm，位置参数λ=1.4的工况下的计算域压力场演化图。

图4(a)～图4(d)为球形空化泡形变为长空化泡的压力场演化图。由于固壁的阻滞效应，空化泡纵向方向受到的挤压大于横向方向受到的挤压，空化泡呈现出被左右拉长的状态。受到液体回弹效应的影响，空化泡右侧形成锥形高压区域。图4(e)～图4(g)显示了空化泡溃灭末期压力场演化的过程，空化泡右侧的高压使得泡壁表面出现轻微凹陷现象，但是由于垂直于曲面固壁右侧的冲击压力没有达到能够击穿泡壁左侧表面的大小，因此空化泡壁没有被完全穿透，而是维持此状态继续溃灭直至泡壁表面完全消失，最后空化泡完全消失的压力场图如图4(g)所示。图4(h)和图4(i)为空化泡溃灭后一段时间内的压力场分布图，通过观察可以发现曲面固壁右侧表面形成圆形高压区，这是造成固壁损伤不可忽略的原因。

图4 曲面固壁附近空化泡溃灭阶段压力场演化Fig.4 Pressure field evolution of cavitation bubble collapse near curved wall

对曲面固壁半径r= 0.3 mm，空化泡初始半径R0= 0.2 mm，位置参数λ= 1.4 的工况进行空化泡溃灭时的速度场演化分析，如图5 所示，其中t=6.89 μs 时刻空化泡完全溃灭。计算域中速度最高的位置用红色箭头表示，速度最低的位置用蓝色箭头表示，箭头的方向即表示流场速度的方向。

观察图5(a)和图5(b)可以发现在空化泡溃灭初期，空化泡上下泡壁收缩的速度较快，而由于曲面固壁的阻滞作用使得靠近曲面固壁一侧的空化泡壁收缩的速度较慢，空化泡右侧泡壁收缩速度明显高于空化泡左侧收缩速度。图5(c)和图5(d)显示了空化泡溃灭末期流场速度的分布情况，图中指向曲面固壁的速度值较大，随着空化泡的溃灭，会产生对空化泡附近曲面固壁的射流冲击，对曲面固壁造成损伤。图5(e)和图5(f)为空化泡完全溃灭之后流场速度分布图，图中指向曲面固壁的速度仍然存在，然而随着时间的推移，指向曲面固壁的速度值越来越小，并且壁面附近的速度方向发生改变，沿着空化泡轴线处形成上下对称的涡流。

图5 曲面固壁附近空化泡溃灭阶段速度场分布Fig.5 Velocity field evolution of cavitation bubble collapse near curved wall

空化泡溃灭时会产生高温高压效应，其会带来空化泡周围温度的变化。图6 为空化泡溃灭过程中温度场的演化情况，其中图6(a)～图6(f)为空化泡溃灭阶段温度场的演化过程，图6(g)和图6(h)为空化泡完全溃灭后一段时间温度场的演化图。初始时刻，曲面固壁半径r= 0.3 mm，空化泡半径R0= 0.2 mm，位置参数λ= 1.4，空化泡内温度与泡外温度均为300 K。

从图6(a)～图6(c)中可以看出空化泡溃灭初期，泡内外温度差异较小，泡内温度上升缓慢。随着空化泡被压缩，泡内温度逐渐升高，泡内外温度差异逐渐增大，空化泡溃灭后期，泡内温度上升速度与空化泡溃灭初期相比较快，如图6(d)到图6(f)所示。然而在此过程中，泡外温度值变化较小，与泡内温度的变化形成鲜明对比。这是由于空化泡体积小，周围流体计算域范围大，当空化泡收缩产生高温效应时，周围流体对其冷却速度较快，导致空化泡周围流体总体温度上升较少。观察图6(g)和图6(h)可以发现空化泡壁完全消失后，随着流场中流体的运动，空化泡溃灭产生的高温区域将向曲面固壁一侧偏移。这将导致曲面固壁受到空化泡溃灭的影响，从而对曲面固壁表面的材料造成高温损伤或其他不可忽略的影响。

图6 曲面固壁附近空化泡溃灭阶段温度场图Fig.6 Temperature field evolution of cavitation bubble collapse near curved wall

2.2 溃灭时间分析

为了进一步考察曲面固壁下空化泡溃灭时间的变化情况，画出了空化泡初始半径R0=0.2 mm，曲面固壁半径r分别为0.1 mm、0.2 mm、0.3 mm时，空化泡溃灭时间随位置参数变化的折线图，如图7所示。图中横轴表示位置参数，纵轴表示溃灭时间，3 条折线分别表示不同曲面固壁半径与空化泡初始半径的比值下的空化泡溃灭时间。

从图7 中可以看出，曲面固壁半径与空化泡初始半径的比值相同时，空化泡溃灭时间随着位置参数的增大而减小，空化泡与固壁间的距离越远，溃灭速度越快。保持位置参数不变，曲面固壁半径与空化泡初始半径的比值越大，曲面固壁对空化泡溃灭过程的阻滞效应越强，空化泡溃灭时间越长。随着位置参数的增大，空化泡逐渐远离固壁面，大尺寸固壁下空化泡溃灭时间下降幅度较大，不同曲面固壁半径与空化泡初始半径的比值所带来的空化泡溃灭时间差异在逐渐减小。

图7 不同曲面固壁半径下溃灭时间随位置参数的变化Fig.7 The change of collapse time with position parameter under different curved wall radius

2.3 射流速度分析

空化泡近壁溃灭产生的射流冲击会对材料表面造成不可忽略的影响，选取空化泡初始半径R0= 0.2 mm，不同曲面固壁半径r，位置参数λ分别为1.1、1.4、1.7、2.0的工况进行空化泡轴上射流速度的讨论，如图8 所示。图中横轴表示空化泡轴线上的位置，纵轴表示空泡轴线上的射流速度。

图8 不同曲面固壁半径下射流速度随位置参数的变化Fig.8 The change of jet velocity with position parameter under different curved wall radius

从图8 中可以发现，射流速度大致变化情况为随着横轴值的增大，空化泡轴上射流速度逐渐增大再逐渐减小，在轴上某一位置会迎来射流速度的峰值，且相同位置参数下，射流速度峰值的位置随着曲面固壁半径与空化泡初始半径比值的增大而减小。出现以上现象的原因为，曲面固壁半径越大，空化泡溃灭点与固壁间的距离越小，射流速度峰值的位置与固壁间的距离也越小。随着位置参数的增大，空化泡与曲面固壁间的距离越远，固壁对空化泡溃灭过程中产生的液体射流的影响越小，空化泡轴上射流速度的峰值越小。另外，从λ= 1.4、λ= 1.7、λ= 2.0 三幅图中可以看出，同一位置参数下，曲面固壁半径与空化泡初始半径的比值越大，固壁对空化泡溃灭产生的阻滞效应越强，空化泡溃灭阶段其轴上射流速度的峰值越大。

2.4 固壁温度分析

本部分主要对空化泡溃灭过程中热力学效应进行讨论，由于实验研究中测量温度变化情况较为困难，且测得数据受较多因素的影响，所得结果与实际结果相比误差较大。而数值模拟可以便捷地调整参数大小，且数值仿真成本较低，故采用数值模拟的方式考察空化泡溃灭的热力学特性。根据仿真结果，重点研究了不同条件对壁面温度的影响。由于模型为轴对称结构，观察了曲面固壁上表面3 个位置点的温度随位置参数的变化情况来探究曲面固壁上的温度变化情况，其中A、B、C三个位置点分布如图9所示。

为探究不同位置参数下曲面固壁上的温度分布情况，画出了固定曲面固壁半径与空化泡初始半径的比值为1，位置参数λ分别为1.1、1.4、1.7、2.0时，曲面固壁上A、B、C三点的温度随时间变化的情况，如图10 所示。其中A、B、C三点在曲面固壁上的具体位置分布如图9 所示，图10 中横轴表示空化泡溃灭的时间，纵轴表示固壁上的温度值。

图9 A、B、C 点位置分布示意图Fig.9 Schematic diagram of location distribution of points A、B and C

图10 不同位置参数下曲面固壁温度分布Fig.10 Temperature distribution of curved wall under different position parameters

从图10中可以看出，固壁温度在空化泡溃灭前期基本保持不变，空化泡溃灭末期，固壁温度急剧上升，在溃灭时刻达到温度最大值，随后固壁温度迅速下降。这是因为溃灭前期流场温度变化较小，对固壁温度产生的影响几乎可以忽略不计。随着空化泡溃灭时刻的到来，空化泡内的高温在溃灭瞬间向外迅速扩散，且溃灭瞬间对壁面产生高压冲击，会对空化泡附近的固壁面造成较大的影响。固壁面在短暂的高温冲击之后，被周围流体迅速冷却，使温度急剧下降。当位置参数相同时，曲面固壁上距离空化泡越近的点，其温度峰值越大，空化泡的溃灭对其温度造成的影响越大。比较不同位置参数下壁面温度变化可得，位置参数越大，曲面固壁上温度峰值越小，即空化泡与曲面固壁间的距离越远，空化泡溃灭时刻到达固壁上的温度峰值越小，空化泡溃灭对固壁造成的高温损伤越小。

由于曲面固壁上C点受空化泡溃灭的影响最大，本文讨论了相同位置参数下曲面固壁右侧C点的温度变化情况，如图11 所示。图中横轴表示空化泡溃灭时间，纵轴表示固壁上的温度值。

图11 不同曲面固壁半径及泡内外压差下曲面固壁温度的变化Fig.11 The change of the temperature of curved wall under different curved wall radius and pressure difference between inside and outside the bubble

由图11可知，曲面固壁温度峰值随着曲面固壁半径与空化泡初始半径比值的增大而增大，且曲面固壁半径与空化泡初始半径的比值越大，固壁温度到达峰值所需的时间越长。固壁与空化泡间的距离相同时，曲面固壁的尺寸越大，空化泡轴上的射流速度越快，其温度到达固壁面时越高。而曲面固壁半径越大，空化泡溃灭时间越长，热量到达固壁所需的时间越长。泡内压强越大，泡内外压差越小，空化泡溃灭的时间越长，导致曲面固壁上温度到达最大值所需的时间越久，壁面温度峰值越小。

3 结论

本文采用VOF 方法对空化泡邻近曲面固壁溃灭的热动力学行为进行研究，将曲面固壁附近高压放电致泡实验结果与数值仿真结果进行对比，从压力场、速度场、温度场、空化泡溃灭时间、空化泡轴上射流速度等方面分析不同物理参数对空化泡溃灭过程的影响。实验所得结果与数值模拟所得空化泡溃灭形态变化相同，验证了数值模拟空化泡邻近曲面固壁溃灭过程的正确性。通过改变曲面固壁尺寸与空化泡初始半径的比值及空化泡与曲面固壁间的距离可以发现，空化泡溃灭时间随着位置参数的增大而减小，曲面固壁半径较小时可以加速空化泡的溃灭。空化泡轴上射流速度的峰值随位置参数的增大而减小，而曲面固壁半径与空化泡初始半径的比值越大，射流速度的峰值越大。

空化泡溃灭时产生的高温在压力和流体速度的影响下，热量沿着射流方向转移到壁面上。位置参数越大，壁面温度越小，壁面上距离空化泡越近的点，空化泡溃灭时对其温度造成的影响越大。曲面固壁的尺寸越大，固壁上温度峰值越大，到达温度峰值所需的时间越长。固壁上的温度峰值随着泡内压强的增大而降低，泡内压强增大，空化泡溃灭的时间变长，随之带来曲面固壁温度到达峰值所需的时间变长。

本文对空化泡近曲面固壁溃灭的热动力学研究，有利于人们深入了解空化泡近壁溃灭机理，并对工程中不同形状固壁附近空化效应的认知提供一些帮助。同时揭示了空化泡近壁溃灭的高温影响，对工程实际中空化在材料表面的应用或规避具有重要作用，为未来研究复杂几何固壁附近空化效应提供参考。