油气两相动压密封动态特性的热流固耦合研究

李世聪，钱才富，李双喜，陈炼

（1 北京化工大学流体密封技术研究中心，北京100029； 2 中国船舶重工集团公司第七一二研究所，湖北武汉430064）

引 言

旋转设备的润滑与密封问题是制约设备性能的关键因素[1-2]，高速工况下的轴承润滑多采用喷射润滑[3]、环下润滑[4]、油气润滑[5]和油雾润滑[6]等微量高效润滑形式。国内外学者对高速轴承腔(如航空发动机轴承腔)内的润滑及油气分布状态做了大量的研究[7-8]，但随着旋转设备的转速和直径的增大，轴承腔内油气两相流体的密封问题也制约了轴承腔的发展。

油气两相动压密封是一种新型的用于高速轴承箱的非接触密封，正常工作时依靠螺旋槽的动压效应实现密封面的不接触，端面无摩擦磨损；同时依靠动压槽的回流泵送效应，实现润滑油的无泄漏。目前已有学者对油气两相动压密封的密封机理和稳态性能作了较为完整的阐述[9-11]，但油气两相密封处于高转速、变压差、变载荷、油气两相混合介质、操作条件多变等复杂工况[12-13]，对密封的动态特性要求非常高。良好的动态特性可防止密封在受到干扰时发生端面接触或密封间隙过大，从而避免密封环碰磨或密封泄漏量过大，既保证密封运转的稳定性，同时对维持高速轴承腔的稳定运行至关重要[14]。

研究密封动态特性的方法主要有步进法[15]、直接数值频率响应法[16]、摄动法[17]以及一些新的耦合方法[18-19]，然而当前大量的密封动态性能的研究建立在摄动法的基础上[20-22]。李双喜等[23]基于高阶函数的有限元方法，分析了高速螺旋槽端面密封的轴向微扰特性以及压缩因数和扰动频率因数对动态系数的影响。摄动法在动压型干气密封[24-25]、T 型槽干气密封[26-27]、静压型干气密封[28-29]、上游泵送密封[30-31]等非接触密封的动态性能研究中也得到了广泛应用。

前述的研究建立在密封环不变形、流体膜温度和黏度不变的假设上，然而随着密封转速、温度和压力的增加，密封环的热变形和力变形能达到微米级[32-33]，引起流体膜厚度的变化，进而影响到流体膜以及密封的性能[34-35]。

本文采用热流固耦合方法，在可压缩流体动态雷诺方程的基础上，考虑油气两相流体物性，通过摄动法得到油气两相动压密封的动态雷诺方程和动态性能参数，并采用有限元法进行求解。研究热流固耦合方法对油气两相动压密封的动态性能的影响，在此基础上分析操作参数(油气比、转速和压差)与动态性能的关系，为分析密封追随性提供理论依据，同时为油气两相动压密封的设计提供指导。

1 分析模型

1.1 几何模型

油气两相动压密封的结构示意图如图1 所示，润滑油从图中喷油嘴喷出，在高速旋转的轴承的作用下与空气混合形成油气两相介质。油气两相流体在轴承腔内的流动形式为油膜、油滴与空气共存，其中液相占密封腔空间体积的2%～20%，油滴颗粒直径大小为3～7 μm[36]。密封端面的外侧为压力较高的油气两相介质，内侧为空气。外侧的油气两相介质在压差的作用下进入密封端面间，形成油气两相润滑的动压密封。动环随轴旋转，动环与旋转轴之间的压紧力由传动轴套提供，石墨静环为补偿环，波形弹簧为补偿元件。采用O 形圈作为辅助密封，保证密封形成封闭空间。

图1 油气两相动压密封示意图Fig.1 Diagram of oil-gas miscible backflow pumping seal

1.2 数学模型

图2 热流固耦合方法和分析模型Fig.2 Thermal-fluid-solid coupling method and analysis model

1.2.1 热流固耦合分析模型 密封端面间的油气两相介质流体膜部分为密封环的变形提供流体膜压力，密封环的变形则会影响到流体膜的膜厚进而改变油气两相介质流体膜压力，密封环变形和油气两相介质流体膜压力互相之间的影响是一个耦合的过程。密封环受到热应力的影响产生热变形，同样会影响到油气两相介质流体膜压力，这使得密封环变形和油气两相介质流体膜压力之间互相影响，最终形成一个新的动态稳定状态，各分析模块之间的影响关系如图2(a)所示。

在MATLAB 中建立密封环与油气两相介质流体膜的有限元耦合模型，油气两相介质流体膜的分布同周期性螺旋槽一样，具有周期性，选取流体膜的单个周期作为分析区域建立有限元模型，如图2(b)所示。选取单个周期可以对流体膜有限元模型进行较为精细的网格划分，保证了网格精细度和计算准确度。

1.2.2 油气两相膜控制方程 由于轴承腔油气两相流体中的润滑油较少，因此在油气两相流体物理模型中，假设油滴均匀分布于空气中，油滴彼此间的相互作用忽略不计，且油滴颗粒互相之间不产生碰撞、破碎或聚合[37]，故而用均相流模型计算油气两相流体的物性，得到的油气两相流体的瞬态雷诺方程为(推导过程参考文献[38])：

式中，pm为油气两相流体膜的压力，μm为油气两相流体的等效黏度，ρm为油气两相流体的等效密度，ω为密封的转速，h为流体膜厚度，rm为油气两相流体计算域的半径，θ为角度，t为时间。

油气两相流体的等效黏度μm和等效密度ρm按式(2)计算：

式中，μgas为空气的黏度，c为油滴在密封腔气体中所占的体积率(即文中的油气比)，ρoil为润滑油密度，ρgas为空气密度。

由于油气两相动压密封的压力较低，因此可以忽略压力对流体黏度的影响，仅考虑温度和黏度的关系。在压力较低时，气体黏度与温度的关系用萨特兰公式计算：

式中，Bgas是与气体性质有关的常数，对于空气或氮气取110.4 K，μgas,0是温度为T0时对应的气体黏度。

1.2.3 温度分布与结构变形 油气两相动压密封的流体膜仅有几微米，因此端面变形和几何型槽结构对油气两相膜性能影响明显。由于端面变形受密封环几何结构形状和尺寸的影响，分析时采用有限元方法求解流体膜控制方程、传热方程及密封环变形[39]。

三维的固体传热分析可以离散为有限元方程式(4)进行描述：

式中，KT为动环和静环的热传导矩阵；T为节点的温度向量；PT为动环和静环的热通量向量。

三维的固体变形可以离散为有限元方程式(5)进行描述：

式中，Km为动环和静环的总刚度矩阵；δ为节点的位移向量；Fm为动环和静环的总载荷向量。

动环和静环均受到表面力(压力、弹簧力)和温差的作用，对于旋转的动环，还额外受到体积力(离心力)的作用。

则密封变形后的流体膜厚度为：

式中，h′0表示变形后的两相流体膜厚；h0表示非槽区油气两相流体膜厚；hg表示槽深；δj是槽深控制开关，δj=0 表示非槽区，δj=1 表示槽区；δp,R、δp,S分别为动环和静环轴向的力变形位移；δT,R、δT,S分别为动环和静环轴向的热变形位移。

1.2.4 密封动力学模型 图3 为油气两相动压密封动态特性分析的动力学模型。利用瞬态雷诺方程式(1)并运用小扰动法推导出微扰雷诺方程组(推导过程参考文献[23])，并引入无量纲变量如下：

式中，rin为流体膜内径；pin为内径处压力；j=r,i（r为密封微扰压力的实部,i为密封微扰压力的虚部）。

图3 油气两相动压密封动力学模型Fig.3 Dynamic model of oil-gas miscible backflow pumping seal

则油气两相流体的无量纲微扰雷诺方程组的表达式分别如式(8)～式(13)所示：

1.2.5 油气两相动压密封性能参数计算 开启力通过对密封端面流体膜的压力场积分得到，单位为N。

动态刚度系数按式(15)计算

动态阻尼系数按式(16)计算

1.3 边界条件

在迭代计算稳态、动态气膜压力时，需要给定边界条件，在迭代中边界条件分为两类。

(1)在密封环的外、内侧，即油气两相流体进、出口处满足

(2)密封端面周向(θ方向)的周期性边界条件

(3)动环和静环的热通量边界条件为：

1.4 油气两相动压密封的结构参数和工况参数

油气两相动压密封环的结构如图1 所示，其结构参数、操作参数和动静环材料属性参数见表1～表3，其中表2 为300 K 下润滑油以及空气的物性参数。

1.5 耦合方法和计算流程

密封端面间的流体膜部分为密封环的变形提供流体膜压力，密封环的变形则会影响到流体膜的膜厚进而改变流体膜压力，密封环变形和流体膜压力互相之间的影响是一个耦合的过程。密封环受到热应力的影响产生热变形，同样会影响到流体膜压力，这使得密封环变形和流体膜压力之间互相影响，最终形成一个新的动态稳定状态，各分析模块之间的影响关系如图4所示。

表1 油气两相动压密封结构参数Table 1 Structural parameters of sealing end face

表2 油气两相动压密封操作参数Table 2 Operating parameters of sealing end face

表3 密封环材料属性Table 3 Material properties of sealing ring

理论上，温度和变形这两个因素是互相影响的，但在一般的实际运用里，温度受到应力和应变的影响很微小，因此可以忽略变形对温度的影响。因此本研究中对密封环进行热分析的过程时，按照温度单向影响变形的计算方式进行耦合分析。

2 结果讨论与分析

2.1 算例验证

在润滑气体为理想气体、密封间隙进出口均采用强制压力出口边界的前提下，Ruan[40]采用有限元法给出了三自由度微扰下螺旋槽干气密封的动态刚度和阻尼数据。本节选用Ruan的计算参数，对本文的动态特性算法进行验证。

图4 密封热流固耦合方法和计算流程Fig.4 Thermal-fluid-solid coupling method and calculation process

案例计算参数：密封环内径rin= 30 mm；外径rout= 42 mm；螺旋槽槽根半径rg= 33.6 mm；螺旋角α= 20°；槽深hg= 5 μm；槽台比为1；气体黏度μ=1.79×10-5Pa·s；大气压力pa= 0.101 MPa；进口压力pout=0.303 MPa。

图5 文献计算结果对比验证Fig.5 Verification of calculated value with literature value

图5 所示为本文使用文献[40]的密封分析参数进行的计算结果对比，从图中可知，本文和文献[40]的计算结果有较好的一致性，这验证了本文中密封动态性能参数求解和程序编写的正确性。

2.2 油气两相动压密封的热流固耦合分析

图6是油气两相动压密封热流固耦合的流体膜性能。图6(a)和(b)是热流固耦合对流体膜厚度的影响，在当前计算条件下，端面的变形为0.3 μm，约为膜厚的11%，对密封性能有一定的影响，因此在分析油气两相动压密封性能时应采用热流固耦合的方法。油气两相动压密封的密封压力较低，但转速较高，密封面的变形主要由温差和转速带来的离心力引起，压差的影响很小。

稳态条件下采用热流固耦合的流体膜温度分布和压力分布分别如图6(c)和(d)所示，压力在螺旋槽根部出现极值，同时温度沿径向逐渐降低。不同的工作膜厚对应着不同的密封开启力，图6(e)是膜厚从1 ～4 μm变化时对应的开启力。当开启力与闭合力相等的时候，即对应着密封的工作膜厚。耦合后端面变形增加，最小流体膜厚度减小。

图6 热流固耦合对流体膜性能的影响Fig.6 Effect of thermal-fluid-solid coupling method on the properties of fluid film

2.3 油气两相动压密封动态性能计算

图7 为耦合前后，密封的主动态性能参数和非对角动态性能参数受到扰动频率Γ影响的规律。其中图7(a)为流体膜无量纲刚度系数受扰动频率Γ影响的规律，考虑热流固耦合后的动态性能参数加“[]”来区别(下同)。从图中可以看出，相比于不采用耦合方法的流体膜，采用耦合方法的主刚度系数(Kzz，Kxx，Kyy)较小，非对角刚度系数(Kxy，-Kyx)基本相同，热流固耦合方法对3 种动态刚度系数的影响程度分别达到了16.9%、16.3%和4.88%。

图7 流体膜动态性能参数与扰动频率的关系Fig.7 Relationship between dynamic parameters of fluid film and disturbance frequencies

在Γ≤1 时，刚度系数对扰动频率Γ的变化不太敏感。在Γ＞1 时，主刚度系数(Kzz，Kxx，Kyy)随扰动频率Γ的增大而线性增加，非对角刚度系数(Kxy，-Kyx)比主刚度系数小得多，并且随扰动频率的增大而略微减小。

图7(b)为流体膜无量纲阻尼系数受扰动频率Γ影响的规律。从图中可以看出，相比于不采用耦合方法的流体膜，采用耦合方法的密封主阻尼系数(Czz，Cxx，Cyy)较大，非对角阻尼系数(-Cxy，Cyx)基本相同，热流固耦合方法对3 种动态阻尼系数的影响程度分别达到了24.3%、31.2%和7.0%。

从图中可以看出，在Γ≤1 时，轴向主阻尼系数Czz比角向主阻尼系数(Cxx，Cyy)和非对角阻尼系数(-Cxy，Cyx)大得多，并且受到扰动频率影响的程度较大(随着扰动频率增大而增速递减)。在Γ＞1 时，阻尼系数皆随着扰动频率的增大而快速减小，并趋于稳定。

综合以上分析可以得知，虽然密封端面变形以及密封平衡膜厚的变化都不大，但对于密封动态性能的影响较大，主刚度系数和主阻尼系数的相对误差最大为16.9%和31.2%，因此在分析油气两相动压密封应采用热流固耦合方法。在Γ≤1 时，扰动频率对流体膜主动态性能参数的影响不大，但在Γ＞1 并逐渐增大时，流体膜主刚度系数略微增大，流体膜主阻尼系数快速减小。故而在后续计算中取Γ=1。

2.4 操作参数对动态性能参数的影响

2.4.1 油气比对密封动态性能参数的影响 油气比对密封动态性能参数的影响如图8 所示，由图可知，油气比的增大，油气两相动压密封的主刚度系数(Kzz，Kxx，Kyy)线性增大，而对于非对角刚度系数(Kxy，-Kyx)而言影响则不是很明显。油气比的增加导致油气两相介质的黏度增加，在同等操作参数下，密封的动压效应更强，因此密封的刚度系数变大。

图8 油气比对密封流体膜动态性能参数的影响Fig.8 Relationship between oil-gas ratio and dynamic parameters of fluid film

对于密封的阻尼系数而言，油气比的增大，主阻尼系数(Czz，Cxx，Cyy)和非对角阻尼系数(-Cxy，Cyx)均有不同程度的增加，其中轴向阻尼系数Czz的增速更快。说明较大流体黏度增加了密封的阻尼，油气两相动压密封能够更快速地使能量耗散。

油气两相流体黏度的增加不仅导致了动压效应的增加，还引起膜厚的增加。尽管膜厚的增加会降低密封的稳定性，但密封的动态刚度系数和动态阻尼系数均有增加，且阻尼系数的增速较快，因此油气比的增加有利于提高密封的动态性能。

2.4.2 转速对密封动态性能参数的影响 转速对密封动态性能参数的影响如图9所示。随着转速的增加，油气两相动压密封的主刚度系数(Kzz，Kxx，Kyy)都呈线性增加，而非对角刚度系数(Kxy，-Kyx)有略微的减小。可以认为这是由高转速引起的较好的动压效应，增强了密封的动态性能，进而提高了密封的动态刚度。

图9 转速与密封流体膜动态性能参数的关系Fig.9 Relationship between rotational speed and dynamic parameters of fluid film

对于阻尼系数而言，在转速较低(3000～6000 r·min-1)的情况下，动态阻尼系数快速增加，但随着转速的继续增大，阻尼系数的增速逐渐减缓，说明高转速时流体膜更稳定。

高转速虽然增强了密封的动压效应，但同时也增加了密封端面特别是动环端面的变形，反而降低了密封的动态性能。但综合考虑转速与变形的关系，高转速有利于油气两相动压密封的动态性能。

2.4.3 压差对密封动态性能参数的影响 压差对密封动态性能参数的影响如图10 所示，由图可知，压差的增加，油气两相动压密封的主刚度系数(Kzz，Kxx，Kyy)缓速增大，非对角刚度系数(Kxy，-Kyx)基本不变。

图10 压差与密封流体膜动态性能参数的关系Fig.10 Relationship between differential pressure and dynamic parameters of fluid film

从数值上看，压差的增加使得动态刚度系数的增幅远超油气比和转速对动态刚度系数的影响，而且压差的计算值最大为0.13 MPa，这是因为本文分析的油气两相动压密封的载荷系数较大，而转速为8000 r·min-1，当压差超过0.13 MPa 时，平衡膜厚计算值小于1.5 μm，较小的膜厚使得密封具有极大的刚度系数，因此密封的主刚度系数(Kzz，Kxx，Kyy)迅速增加。

对于阻尼系数而言，随着压差的增大，密封的主阻尼系数(Czz，Cxx，Cyy)迅速减小，当压差较大时，主阻尼系数(Czz，Cxx，Cyy)甚至会出现负值，而非对角阻尼系数(-Cxy，Cyx)略有增加。主阻尼系数减小意味着高压差下流体膜的稳定性变差，使得油气两相动压密封更容易发生失稳。

压差增加导致闭合力增加、膜厚减小、稳态刚度增加，导致动态刚度增加、动态阻尼减小，密封更易发生失稳。

3 结 论

本文建立了油气两相雷诺方程，采用热流固耦合方法并用有限元法求解，通过摄动法研究了油气两相动压密封的刚度系数和阻尼系数。不仅为分析密封追随性提供理论依据，同时也为油气两相动压密封的设计提供指导，得到的结论如下。

(1)采用热流固耦合方法计算的油气两相动压密封端面变形约为膜厚的11%，对密封的动态性能的影响较大，密封变形对刚度系数和阻尼系数的影响程度最大为16.9%和31.2%，因此在分析油气两相动压密封时应采用热流固耦合方法。耦合后油气两相膜的动态刚度降低，动态阻尼增加。从振动学的角度来说，密封端面变形提高了油气两相动压密封的动态性能。

(2)油气两相流体中的润滑油增加了工作介质的黏度，同等操作条件下，油气两相动压密封具有更好的动压效应、更大的动态刚度系数和动态阻尼系数，相较于干气密封，油气两相动压密封的动态性能更好。随着油气比的增加，密封的动态刚度系数和动态阻尼系数均有不同程度的增加，而动态阻尼系数的增速更快，有利于提高油气两相动压密封的动态性能。

(3)高转速使油气两相动压密封的主刚度系数和主阻尼系数增大，有利于密封的动态性能。而压差的增加使主刚度系数增大、主阻尼系数减小，不利于密封的动态性能。因此油气两相动压密封适用于高转速、低压差的油气润滑工况。

符 号 说 明

Bgas——与气体性质有关的常数，对于空气或氮气取110.4 K

Cxx,Cyy——角向主阻尼系数

Cxy,Cyx——非对角阻尼系数

Czz——轴向主阻尼系数

c——油滴在密封腔气体中所占的体积率

cp——比热容，J·kg-1·K-1

E——弹性模量，GPa

FC,FO——分别是密封的闭合力、开启力，N

Fm——动环和静环的总载荷向量

h,h′——分别为流体膜厚度、流体膜厚度瞬态变化量，μm

hg——槽深，μm

h0,h0′——分别为非槽区油气两相流体膜厚、变形后的两相流体膜厚，μm

KT,Km——分别是动环和静环的热传导矩阵、总刚度矩阵

Kxx,Kyy——角向主刚度系数

Kxy,Kyx——非对角刚度系数

Kzz——轴向主刚度系数

Ng——槽数

PT——动环和静环的热通量向量

p′——流体膜压力瞬态变化量，Pa

pin,pout——分别是轴承腔外压力、轴承腔内压力，MPa

pm——油气两相流体膜的压力，Pa

pzi,pzr——分别为气膜动态压力的实部、气膜动态压力的虚部，MPa

p0——油气两相流体膜的稳态压力，Pa

rb——平衡半径，mm

rg——螺旋槽根径，mm

rin,rout——分别是密封面内径、密封面外径，mm

rin,R,rin,S——分别是动环、静环密封面内径，mm

rm——油气两相流体膜计算域的半径，mm

rout,R,rout,S——分别是动环、静环密封面外径，mm

r1,r2,r3,r4——分别为静环O 形圈槽半径、静环外径、轴径、动环尾部外径，mm

T——温度，K

ΔT——温差，K

V——补偿环受到微扰的频率，Hz

u——泊松比

α——螺旋角，(°)

δj——槽深控制开关，δj=0：非槽区，δj=1：槽区

δp,R,δp,S——分别是动环轴向的力变形位移、静环轴向的力变形位移，μm

δR,δS——分别是动环厚度、静环厚度，μm

δT,R,δT,S——分别是动环轴向的热变形位移、静环轴向的热变形位移，μm

γ——槽堰比

λ——热导率，W·m-1·K-1

μm,μgas,μgas,0,μoil——分别是油气两相流体的等效黏度、空气的黏度、温度为T0时对应的空气黏度、润滑油黏度，Pa·s

θ——角度，(°)

ρ,ρm,ρgas,ρoil——分别是固体密度、油气两相流体的等效密度、空气密度、润滑油密度，kg·m-3

σ——动环尾部倾角，(°)

ω——密封的转速，r·min-1

下角标

in,out——分别是内径侧、外径侧

m,gas,oil——分别是油气两相流体、气体、润滑油

R,S——分别是动环、静环

r,i——分别是密封微扰压力的实部、密封微扰压力的虚部