摘要：运行工况的瞬时变化严重影响密封性能。利用Matlab建立密封环端面间隙液膜三维模型，采用有限差分法离散基于JFO空化边界条件的雷诺方程，应用SOR迭代求解液膜压力分布，进一步耦合求解雷诺方程与瞬态动力学方程，分析工况连续变化及压力扰动对密封瞬态特性的影响。结果表明：相比于转速瞬时变化，压力瞬时变化过程中挤压效应对密封性能的影响更为显著，密封端面趋近速度越大，由液膜挤压产生的承载能力越高，端面流体被排出的速度越大；压力瞬时变化易引发静环轴向速度振荡，压差越大，振荡幅值越大；压力扰动情况下，空化率与泄漏量急剧突变后趋于稳定，压力突升相比于压力突降更易恢复稳定状态；摩擦扭矩在变工况过程中平稳变化，无较大幅度波动。
关键词：液膜密封；非定常工况；瞬态性能；动力学分析
引 言
全液膜润滑非接触式机械密封（简称液膜密封）有效改善了密封端面的润滑状况和操作稳定性，在石油、化工、冶金等行业用流体机械中不断得到推广应用[1-4]。尽管端面液膜有效降低了密封环表面的摩擦磨损及端面冷却消耗，但在实际运行过程中工况的瞬时变化将降低补偿环的追随性，使泄漏量变化加剧，严重影响密封系统稳定性。因此，有必要开展液膜密封在非定常工况下瞬态特性的研究，提高其在复杂工况下的稳定性及安全性。
目前，学者对非接触式机械密封的研究多集中于稳态及动态特性方面[5-7]，鲜有涉及瞬态启停行为及润滑状态转变[8-10]，非定常工况的瞬态性能研究多针对于干气密封[11-13]，尚未有液膜密封变工况瞬态性能的报道。Salant等[14-15]考虑接触和变形等因素，建立了预测干气密封和液膜密封瞬态启停的数学模型。Green等[16]确立了在定子和转子强迫失调时干气密封瞬态响应的计算方法，对于刚度较大的液膜密封其响应规律仍待进一步研究。刘录等[17]探讨了温度变化条件下密封端面发生闪蒸时对动环开启性能的影响，对转速与压力的瞬时扰动缺乏详细研究。Zhao等[18]考虑旋转密封蠕动特征，对转速冲击工况的密封动态特性进行了研究，未考虑压力冲击对密封行为的影响。陈汇龙等[19]研究了上游泵送机械密封变工况压力脉动特性对压力场的影响，对其进一步导致的密封系统稳定性尚未研究。Varney等[20]分析了非接触机械密封运行过程中受到的振动与冲击现象，未考虑微观气蚀在膜厚波动下对膜压的影响。胡松涛等[21]研究了干气密封内外压差波动工况下的瞬态特性，但其动力方程的求解过程仍不够完善。
总之，有关非定常工况下液膜密封瞬态行为仍缺乏精确的研究，压力及转速的瞬时变化对密封参数性能的影响也未有详细讨论，同时，端面液膜空化作为其重要影响因素，目前考虑JFO空化的含时间项雷诺方程求解算法也过于烦琐，一定程度降低了计算效率。鉴于上述问题，本文利用Matlab建立密封环端面间隙液膜三维模型，采用有限差分法离散基于质量守恒空化边界条件的雷诺方程，进一步耦合求解瞬态动力学方程，分析非定常工况下液膜密封的瞬态特性。
1 数学模型
1.1 密封系统模型
图1为液膜密封结构[22]，该结构由动环、静环、副密封（O形圈）、弹簧座等组成，弹簧座固定于腔体上，静环安装在弹簧座上，可通过弹簧做轴向运动，动环随转轴旋转，外径侧为高压密封介质，内径侧为缓冲液，静环刻有流体动压槽，动环旋转时可产生较强的动压使密封端面产生微米级间隙，外径侧高压密封介质在压差作用下向内径侧的泄漏被螺旋槽泵送效应反输回高压侧[23-24]，在减小泄漏的同时降低了端面摩擦磨损，改善了端面润滑状态，有效延长了使用寿命。
图1 液膜密封结构
Fig.1 Schematic diagram of liquid film seal structure
1.2 端面几何结构
选取工业应用中较为典型的螺旋槽液膜密封为研究对象，其密封环端面结构如图2所示，静环端面包括动压槽、密封堰及密封坝，ri、ro分别为螺旋槽的内径与外径，rg为螺旋槽的槽根半径，θ为转角，a为螺旋角，动压槽深度为hg，θa为台宽，θw为槽宽，ζ为螺旋槽周向槽宽与台宽之比，pi为内径处压力，po为外径处压力。
图2 密封环端面结构
Fig.2 Geometry schematic diagram of seal ring
槽型线为对数螺旋线，其在极坐标下表达式如式(1)所示，密封几何参数与工况参数如表1所示。
r=rgeqtana (1)
2 研究方法
2.1 基本假设
为便于方程求解，根据液膜在密封间隙中的流动特征，对密封间隙内薄膜流场做如下简化[25-26]：
（1）忽略惯性力的影响；
表1 几何结构与工况参数
Table 1 Geometrical and operational parameters
（2）密封端面间的产热较小，密封动环与静环为刚性体，为等温状态，不计表面变形；
（3）密封液膜中的流体做层流运动，液膜中没有旋涡和紊流，而且符合牛顿黏性定律；
（4）摩擦副曲率半径远大于液膜厚度，忽略液膜曲率影响；
（5）动环表面与黏附于表面的流体无相对滑动。
鉴于上述假设，推导出适用于全液膜润滑的液体润滑机械密封端面液膜流场的控制方程，如式(2)所示
(2)
式中，Φ为液膜密度比，Φ=ρ/ρc，ρ为端面间任意点液膜密度，ρc为密封介质的密度。
2.2 边界条件
强制性边界条件
内径处：p=pi，F(r=ri,q)=1
外径处：p=po，F(r=ro,q)=1
周期性边界条件
JFO空化边界条件[27]
2.3 控制方程离散
采用有限差分法对离散雷诺方程，采用向后欧拉法进行时间项离散，得到各系数表达式为
， 
，
，
，
，
，
，
式(2)在离散整理后得到如式(3)、式(4)所示的压力及液膜密度比表达式
(3)
(4)
利用超松弛迭代法求解端面液膜压力场及密度场分布
(5)
(6)
当迭代计算结果满足式(7)时计算结果满足精度要求，完成压力场求解。
(7)
2.4 运动微分方程求解
对图1所示密封结构图进行适当简化得到图3所示的液膜密封系统的动力学模型[28-29]，Kxt、Dxt分别为考虑O形圈及弹簧作用下系统的刚度和阻尼[30]。右侧为动环，左侧为浮动静环，在弹簧力作用下可做轴向移动。
图3 密封动力学模型
Fig.3 Dynamic model of mechanical seal
密封系统动力学方程如式(8)所示
(8)
式（8）中各参量计算式如下。
闭合力为
液膜承载能力为
运动微分方程可以写为
(9)
(10)
采用Newmark-β法对动力学方程积分[31]，尽量减小时间步长对计算结果影响。
(11)
(12)
其中，β=0.5；γ=0.25。
2.5 计算流程
液膜密封瞬态性能分析的计算流程如图4所示。取端面间隙流体为计算对象，因其呈周期性排布，为提高计算效率，选择整个密封端面的1/N周期作为计算区域。确定表1所示密封基本参数之后，式(2)与式(8)在工况控制函数命令下耦合求解，当压力场迭代满足相对误差的精度要求时，完成当前时间步下内循环的计算；当时间满足工况控制函数命令的时间要求时完成整个全局周期的计算并输出参量变化规律。网格无关性验证如图5所示，当单周期径向与周向网格数量大于60×60个时，液膜承载力及泄漏量的求解结果基本趋于稳定，在保证计算结果的准确性与稳定性的同时，为节省计算时间，提高计算效率，选取60×60的网格密度，选择时间步长为1×10-6 s。
图4 计算流程
Fig.4 Calculation flow chart
图5 网格数对密封参数的影响
Fig.5 Influence of mesh number on seal parameter
2.6 特性参数
在每个时间步雷诺方程迭代收敛之后，通过式(13)～式(15)可得到密封特性参数。
空化率
(13)
体积泄漏量
(14)
摩擦扭矩
(15)
3 结果分析
为模拟液膜密封装置实际运行的非稳态工况，本文分别研究了转速与外径侧（密封腔）压力在1 s内连续变化、密封腔压力瞬时扰动两种工况下的密封瞬态特性（表2）。初始时刻转速为2000 r·min-1，密封环内径侧与外径侧压力分别为0.1 MPa与0.4 MPa。压力扰动的变化规律为：在内径压力保持不变的情况下，外径侧压力由稳定状态下的0.4 MPa在0.002 s时分别瞬时变化至0.2、0.3、0.5、0.6和0.9 MPa，研究密封在多干扰工况下的密封瞬态特性。
表2 工况变化规律
Table 2 Changes of working conditions
3.1 算法正确性验证
采用电涡流传感器测量液膜密封浮动静环轴向位移，试验装置如图6(a)所示，其测量原理如图6(b)所示，铜片安装于动环，与动环端面水平，电涡流传感器探针固定于安装板，通过测量与铜片之间距离预测液膜厚度，探针测量精度0.1 μm，采样频率为500 Hz，0.05 s内选取25组数据并进行平均处理，得到一组数据点，通过变频器控制电机转速按预设时间匀速增加，得出膜厚随转速的变化规律。因腔体压力调节存在一定的滞后性，因此试验中仅研究转速连续变化对膜厚变化规律的影响。
图6 液膜密封试验装置
Fig.6 Liquid film seal testing device
图7 膜厚理论计算与试验对比
Fig.7 Comparison of film thickness between theoretical calculation and experimental values
图7示出了JFO空化边界、Sommerfeld边界两种瞬态计算方法与试验结果的比较，由图可得，随转速增加，螺旋槽流体动压效应不断增强，使得流体膜压不断升高。当端面液膜承载力大于闭合力时会使弹簧压缩，进一步使得液膜厚度变大，一定程度上削弱了动压承载力，密封系统又达到新的平稳位置。此外，因空化会提高液膜承载力，所以二者膜厚计算值均大于Sommerfeld条件的计算结果。转速越高，空化效应对液膜承载能力的影响越大，因此两种计算方法下膜厚差值不断扩大，因存在假设条件及试验环境限制，理论计算值与实测值存在一定程度误差，但整体变化趋势依然保持一致。另一方面也验证了该方法在求解液膜密封瞬态特性问题的准确性。
此外，参照文献[18]中研究对象的基本参数及运行工况，暂不考虑浮动静环滑动摩擦力及微凸体接触力的影响，将理论计算结果与其进行比较，结果如图8所示，文献计算结果与本文存在一定程度的差异，但两者最大误差在6%以内，且变化规律基本一致，证明了计算结果的准确性。
图8 膜厚变化规律对比
Fig.8 Comparison of film thickness changes
3.2 挤压效应对密封性能影响
建立稳态特性下JFO空化边界与不考虑空化的Sommerfeld边界条件的计算模型，并与瞬态JFO计算模型进行对比，研究工况变化过程中挤压效应密封性能的影响。图9(a)为3种计算方法下膜厚随转速的变化规律，由图可知，在内外侧压力保持不变前提下，考虑JFO空化边界的稳瞬态计算结果几乎重合。图9(b)为空化率与泄漏量在转速变化过程中计算结果比较，由图可知，稳瞬态两种计算方法的空化率和泄漏量基本一致，转速的瞬时变化过程中挤压效应对密封性能影响微弱。
图10(a)为外径侧压力瞬时变化时膜厚的变化曲线，与图9(a)相比较可以发现，挤压效应对液膜厚度变化有明显影响，虽然稳态计算下JFO与Sommerfeld计算差值保持平稳，但是两者与瞬态计算结果的差距逐渐增大。外径侧压力增加引起静环速度波动，端面液膜受到挤压作用，流体压力相互挤压迫使流体膜压强升高，静环速度波动越大，产生的动压承载能力越大，从而与稳态计算有较大差距。压力变化下空化率及泄漏量的稳态与瞬态计算结果比较如图10(b)所示，JFO与Sommerfeld两种边界条件下计算所得泄漏量基本一致，两者与瞬态计算结果差异较大，这是由于压差连续增加易引起静环扰动及速度波动，进一步使端面液膜发生间歇性挤压或扩张，将端面液体排出或吸入。稳态与瞬态计算所得空化率仍存在较大差异，变化趋势也不同，随压差的增大，稳态计算方法得到的空化率不断增加，而瞬态计算结果得到空化率在初期增长幅度较大，之后趋于平缓，这是膜厚的减小对空化的促进作用及外径侧压力增加对空化的抑制等综合原因造成的。
图9 参量在转速变化稳态与瞬态结果比较
Fig.9 Comparison between steady state and transient state of parameters during speed change
3.3 密封瞬态特性变化
3.3.1 膜厚变化 图11(a)为4种工况下膜厚随时间的变化规律，由图可知，随转速增加，膜厚呈线性增大趋势，压差增大，膜厚呈线性衰减趋势，转速及压差越大，膜厚间差值越大。相比之下，外径侧压力的变化会对闭合力产生较大影响。定转速下，随着压差逐渐增加，闭合力不断提高，因此液膜厚度不断减小以形成更高的承载能力。
图10 参量在压力变化稳态与瞬态结果比较
Fig.10 Comparison between steady state and transient state of parameters during pressure change
图11(b)为压差瞬时扰动过程中液膜厚度的变化情况，由图可知，在0.002 s受到不同程度的压力冲击后，膜厚经大幅变化后逐渐趋于稳定状态，压力突升相比于压力突降更易导致密封环大幅位移，但是需要更长的时间才能恢复至稳定状态。这是由于压差增大使液膜厚度减小，膜厚越小，液膜刚度越大，系统的稳定性越好。
3.3.2 速度变化 图12(a)为静环轴向速度随时间变化规律，由图可知，压差改变更易引起速度振荡，随着外径侧压力增加，波动幅值不断减小，波动频率下降；转速的瞬态变化对速度影响较微弱，基本保持平稳趋势。PA工况速度振荡幅值较小，振荡频率也少于PB工况，两种工况下静环轴向速度呈缓慢减小状态，转速改变导致速度有微小起伏，且转速越大，波幅与波数越多，外径侧压力升高导致液膜厚度减小，液膜刚度增大，浮动环移动速度较大进一步导致流体膜阻尼升高，降低了速度振荡频率，同时，PB工况空化率也在0.6 s之后趋于平缓，说明空化效应对液膜厚度及动压承载能力的影响保持稳定。因此，使整个速度波动逐渐衰减，最终达到稳定状态。
图11 膜厚变化规律
Fig.11 Variation of film thickness
图12(b)为受不同压差扰动下轴向速度变化情况，外径侧压力变化幅度越大，速度突变值越大，同等变化幅值下，压力突降下引起的速度瞬时增量基本一致，但是恢复至稳定的时间明显较长。说明在受压差扰动情况下，压力突升相比于压力突降轴向速度更易恢复稳定状态。
3.3.3 空化率变化 不同工况下空化率变化规律如图13所示，由图13(a)可知，WB工况下空化率随时间变化最为明显，转速越大，空化率越大，PA与WA的变化规律基本一致，PB工况下，空化率先增大后趋于平缓，说明当压差增至一定程度后，在转速不变的情况下，空化率基本保持平稳。图13(b)为压差突变下空化率的变化情况，由图可知，在受压差扰动之后，空化率发生急剧突变后逐渐趋于平稳，在0.003 s时多条曲线汇集于一点，该处空化率与初始扰动前空化率基本一致，然后各曲线变化幅度逐渐衰减至稳定状态。与其他变化曲线不同的是，外径侧压力突变为0.2 MPa时引起的空化率瞬态增值最大，压力突降首先引起空化率高于初始状态，经历一段时间变化之后逐渐减小，低于初始值后趋于稳定。这是由于虽然外径侧压力的减小一定程度上更易诱发空化，但是由图11(b)发现压力减小的同时也导致了液膜厚度的不断增加，因此空化区域呈减小趋势，说明膜厚在一定程度上为影响空化的主要因素。
图12 轴向速度变化规律
Fig.12 Variation of axial speed
3.3.4 泄漏量变化 图14(a)为4种工况下泄漏量随时间的变化规律，由图可知，随转速增加，膜厚呈逐渐线性增大的趋势，动静环端面间隙增大导致泄漏量也随之增加，同时转速升高一定程度上提升了螺旋槽泵送能力，使密封端面流体输送能力不断增强，压差增大，膜厚呈线性衰减趋势，端面间隙不断减小，与之对应的泄漏量逐渐减小；转速及压差越大，泄漏量间差值越大。
压差扰动下泄漏量变化规律如图14(b)所示，当压差减小时，泄漏量由负值缓慢增加至正值，泄漏方向由压差流占主导的外径向内径侧泄漏，被螺旋槽泵送效应引起的黏性剪切流反输回高压侧。在外径压力为0.2 MPa时，稳定后泄漏量近似为零，已基本达到零泄漏的理想状态。当外径侧压力瞬时增大时，泄漏量有较大的扰动量，压差增大相比于压差减小更易达到稳定状态，历时约为0.002 s。
图13 空化率变化规律
Fig.13 Variation of cavitation rate
3.3.5 摩擦扭矩变化 由图15(a)可知，摩擦扭矩随转速及压差瞬时增加均呈近似线性增大趋势，转速增值越大，摩擦扭矩变化越明显，由式(15)可知，摩擦扭矩随转速的变化主要是因受到转速及空化的影响，而压差变化影响液膜厚度，进一步使摩擦扭矩产生变化。PB工况摩擦扭矩增幅最大，约为43 %。相比之下，WA工况摩擦扭矩增幅最小，仅为16 %。图15(b)为摩擦扭矩随压差扰动的变化规律，其状态与膜厚变化规律相似，受扰动压力影响，摩擦扭矩呈抛物线递增或递减后逐渐至稳定状态，外径压力为0.3 MPa时变化最为微小，当外径侧密封介质压力为0.2 MPa时，内外径压差较小，流体膜的承载能力大部分由动压效应形成，膜厚较大，使得端面润滑状态为稳定的全液膜润滑，因此摩擦扭矩明显偏小。
图14 泄漏量变化规律
Fig.14 Variation of leakage rate
4 结 论
（1）转速瞬时变化工况下瞬态计算得到的密封性能参数与稳态结果基本一致。与之相比，压力瞬时变化过程中挤压项对密封参数求解的影响更显著。
（2）压力瞬时变化易引发静环轴向速度振荡，压差越大，振荡幅值越大；受压差扰动情况下，压力突升相比于压力突降更易恢复稳定状态。
（3）当外径侧压力瞬时增大时，空化率与泄漏量急剧突变后趋于稳定；摩擦扭矩在变工况过程中变化平稳，无较大幅度波动。
（4）本文提出考虑JFO空化边界液膜密封瞬态特性计算方法，主要计算了非稳态工况下密封性能的变化规律，为整个密封装置的现场运行提供有效指导，同时可进一步开展参数优化，提高密封系统运行稳定性。此外，考虑振动冲击及倾斜的瞬态动力学分析仍待进一步研究。
图15 摩擦扭矩变化规律
Fig.15 Variation of friction torque
符 号 说 明
