Page 52 - 《应用声学》2022年第4期

550 2022 年 7 月

接求解流场中的物理量会使其在界面处发生间断， 半径时始终保持曲面固壁右侧一点固定。为减少压
导致计算的发散或计算结果的虚假，为保证界面处 力边界条件对溃灭过程的影响，将计算域大小设为
所要求解物理量值的合理性，对这类流动进行数值 30R 0 × 30R 0 。计算域的四周均为压力出口边界，曲
研究时，需要进行足够精确的界面捕捉。VOF方法 面固壁表面为无滑移固壁边界条件。
可以较好地保持流体质量守恒并且可以处理界面 在数值模拟软件中导入网格文件，由于 Gam-
拓扑结构的变化 [26] ，因此本文采用 VOF 方法来计 bit 中构建模型结构时默认的长度单位为 m，因此
算空化泡运动过程。 要根据实际目标需求对网格的尺寸进行重新定义
并检查网格质量。采用基于压力的瞬态求解器，选
1.1 控制方程
用多相流模型中的 VOF 模型进行计算，开启能量
空化泡流场中的流动方程由连续性方程、动量
方程。压力 -速度耦合方式运用 PISO 算法，空间离
方程及能量方程组成，表达式分别为 [27]
散方法中，梯度项使用 Least Squares Cell Based 格
∂ρ
+ ∇ · (ρU) = 0, (1) 式，压力项使用 PRESTO! 格式，体积分数项使用
∂t
∂ (ρU) Geo-Reconstruct格式，密度项、能量项、动量项均采
+ ∇ · (ρUU)
∂t 用 Second Order Upwind 格式 [29] 。计算的时间步
= − ∇p + ∇ · (µ∇U) + σκn, (2) 长设置为 10 −8 s，开启残差曲线监控即可开始迭代
∂ (ρE) 计算，观察计算结果是否收敛。
+ ∇ · [U (ρE + p)] = ∇ · (k∇T) , (3)
∂t
其中，ρ表示流体的密度，µ表示流体的黏度，U 表示 ԍҧѣ԰᣸ႍ
๯ʹ
流场速度，p 表示流场压力，σ 表示表面张力系数，κ
表示表面曲率，n 表示两相流界面中的单位法向量， ԍ y ԍ
ҧ R  ҧ
T 表示温度，E 表示总能量，k 表示流体的传热系数。 ѣ ѣ
԰ r x ԰
VOF 方法通过一个标量的输运方程来描述各 ᣸ d ቇӑจ ᣸
ႍ ႍ
相的体积分数，表达式为 [28] ڍܞ᣸ႍ
d
λ=
Dα ∂α R 
= + (U · ∇) α = 0, (4)
Dt ∂t
ԍҧѣ԰᣸ႍ
其中，α 表示构成相的体积分数大小，当α = 0 代表
该网格处于空化泡内部，当 0 < α < 1 代表该网格 图 1 计算模型
处于空化泡外即液体内部，当 α = 1 代表该网格处 Fig. 1 Computational model

于两相交界面上。
1.3 网格无关性验证
各相的密度以及黏度可分别表示为
不同的网格数量会产生不同的计算结果，为
ρ = α 1 ρ 1 + α 2 ρ 2 , (5) 使计算结果误差较小，需要选取恰当的网格数。

µ = α 1 µ 1 + α 2 µ 2 , (6) 在比较不同网格数的计算结果时，采用如下初始
计算条件：曲面固壁半径 r = 0.2 mm，空化泡半径
其中，α 1 + α 2 = 1，α 1 代表液体的体积分数，α 2 代
R 0 = 0.2 mm，无量纲位置参数λ = 1.1，液体中压强
表气体的体积分数。
6
5
p l = 3 × 10 Pa，空化泡内压强p v = 1.01 × 10 Pa。
1.2 计算模型 图 2 为空化泡横轴上射流速度的变化情况，其
本文采用 VOF 方法研究空化泡邻近曲面固壁 中横轴表示空化泡轴上的位置坐标，纵轴表示轴上
溃灭的过程，几何模型如图 1 所示。曲面固壁为球 的射流速度。观察发现射流速度在某一位置达到最
形，半径设为 r，球形空化泡的初始半径设为R 0 。为 大值，且不同网格数下的空化泡轴上射流速度差异
方便描述空化泡与固壁间的距离，引入无量纲位置 较小，几乎可以忽略不计。综合考虑计算时间、计算
参数 λ = d/R 0 ，其中 d 为空化泡中心到曲面固壁右 资源等的限制，最终选取网格数为200000的工况进
侧一点的距离，在数值模拟过程中，改变曲面固壁的 行数值模拟。