Page 6 - 摩擦学学报2025年第8期

1104 摩擦学学报（中英文） 第 45 卷

式中：p为液膜压力；h为液膜厚度；U 和U 分别为动环 常明显的对流占优现象，流线迎风(Streamline-Upwind/
x
y
沿x和y方向的速度；F 为流体膜单位体积内的液相分 Petrov-Galerkin，SUPG)有限元方法(FEM)对于处理对
L
数；ρ和μ分别为均相流体的密度和动力黏度；其表达 流占优问题具有明显的优越性，也可有效地改善流体
式分别为 膜相变引起的数值稳定性问题. 故采用流线迎风有限
元法可有效解决此现象. 由于式(1)及式(2)均为非线性
(3)
ρ = F L ρ L + F G ρ G
方程，因此采用Newton-Raphson方法求解，求解步骤
(4)
µ = F L µ L + F G µ G
如图2所示.
式中：F 为气相体积分数，F =1−F ；ρ 和ρ 为液氧和
L
G
L
G
G

气氧的密度；μ 和μ 为液氧和气氧的动力黏度；S为源 Start
L G
项，表示相变过程中汽液相之间的传质速率 [25-26] ；当流
Initialize the flow domain
体由液相转变为汽相时，气泡变化引起的相变速率可 (k=0)
with F L
表示为
Solve Reynolds equation to get

2 ˙ (k)
S = −C 1 Nρ G 4πR R b (5) the initial values, p
b
Update the mixture viscosity
当流体由汽相转变为液相时，气泡变化引起的相 and density with F L (k=k+1)
Solve equation (2) to get liquid
变速率可估计为 volume fraction, F L (k)

(6)
2 ˙
b
S = C 2 Nρ G 4πR R b
Check convergence No
式中：C 和C 为经验系数，取0.05；R 为气泡直径，取 of p and F L
b
1
2
[27]
1 μm ；N为气泡的实际数量，其表达式为 Yes
 Output the results
F L
 ,
 p < p sat
 4

 π R 3 b

3
N = (7) End
 F G
 ,
 p > p sat
 4

 π R 3 Fig. 2 Flow chart of numerical solution for the
3 b hydrodynamic-hydrostatic mechanical seal
基于气泡动力学理论Rayleigh-Plesset方程(RPE)， phase change lubrication model
图 2 动静压机械密封相变润滑模型数值求解流程图
忽略表面张力与黏性项，RPE取一阶近似则有

˙
R b = 2 |p sat − p| (8) (1) 对如图1(a)中所示的周期性计算域进行三角
3 ρ L
形网格划分，单元数共8 700个.
式中：p 为流体介质温度所对应的饱和蒸汽压. (2) 假定计算域液膜液相体积分数F 的初值.
sat
L
为求解上述方程以获得液膜压力和气液相比，需 (3) 基于有限单元法并采用牛顿迭代法对雷诺方
给定密封环内外径压力边界条件和液相体积分数F L 程(1)进行求解，获得液膜压力场分布.
的值，其表达式为 (4) 基于有限单元并结合牛顿迭代法对传质控制
 ®
1, p i > p sat 方程(2)进行求解，获得端面液相体积分数F .
 r = r i , p = p i , F L = L


 0, p i ⩽ p sat (5) 基于方程(3)和(4)，更新流体介质的动力黏度
(9)
®
 1,
 p o ⩾ p sat μ及密度ρ.
r = r o , p = p o , F L =

0, p o < p sat (6) 计算前后2次液膜压力和液相体积分数值的
式中：p 为端面外径侧压力；p 为端面内径侧压力. 由 误差，判定是否满足收敛准则，满足则程序结束，计算
o
i
[24]
于均压槽节流孔处压力脉动较小 ，故可认为孔区压 密封性能参数；否则返回步骤3进行重复计算.
力为一定值p ，此外针对如图1(a)和1(b)所示的周期性 本文中基于C语言编制了有限元求解程序，求得
o
计算区域，需赋予周期性边界条件. 密封端面液膜压力后，可以求出密封的稳态性能参
在高速工况下，式(1)及式(2)所示Reynolds方程和 数，其中开启力F open 、摩擦系数f和泄漏率Q的计算公
传质控制方程等式中流项数值非常大，方程表现出异 式分别如式(10~12)所示. 其中，Ω为二维液膜计算域；