氦气介质干气密封热-流固耦合建模及性能分析*

黄伟峰 王伟达 刘 莹 刘向锋 何 强

(清华大学摩擦学国家重点实验室 北京 100084)

氦气透平在高温气冷堆核电站[1]、航天等领域[2-3]具有重要应用。对于氦气透平，轴封技术是其关键技术和难点之一，而干气密封由于其低泄漏、高可靠性特性，是一种非常理想的氦气透平轴密封型式。干气密封是一种非接触式气体润滑密封[4]，其工作性能与基于流场、传热、变形的综合仿真分析设计能力密切相关。

自干气密封于20世纪60年代出现以来，大批专家学者在干气密封的理论建模及仿真计算方面开展了研究。相关研究大致可分为两类。其中之一是基于等温假设的纯流场分析计算[5-8]，由于流场分析是干气密封性能分析的基础，因此相关研究较多。采用的方法包括用于螺旋槽密封的解析法[5]，基于有限元法[6]、有限差分法[7]等的全数值法，以及结合了解析法和离散方法特征的半解析法[8]等。另一类是考虑了传热过程、固体变形两者之一或者全部的多场耦合分析，这类方法考虑的因素完整，对干气密封性能预测更加准确，是干气密封分析方法的发展方向。LEBECK[9]在其著作中从求解机械密封环端面位移的圆环理论出发，构建了基于线性位移分布的多场耦合模型。THOMAS等[10]通过数值求解流场、传热、变形场，建立了高压干气密封的多场耦合模型。BAI等[11]对高压及高转速工况下的干气密封进行了耦合分析。FAIRUZ等[12]面向工程中的密封结构，在较为准确的几何模型和边界条件下进行了单向热-流固分析。约翰克兰等代表性国际密封企业，开发了自己的机械密封多场耦合分析软件，但这些软件仅在企业内部使用，软件中采用的算法和模型的具体内容未发表在公开文献中。总体而言，可以看到目前已有一些关于干气密封的多物理场模型的研究，为干气密封耦合分析方法的发展奠定了很好的基础。但将这些模型用于干气密封设计尚存在一些问题：大部分耦合计算的研究采用简化边界条件，如补偿环完全固定；部分研究未纳入干气密封轴向开启力和闭合力的平衡，在给定膜厚条件下进行计算；部分研究的能量方程求解尚存在一些不足之处等。相关简化处理方法往往可以大大降低耦合分析的收敛难度，但与干气密封实际受力状态并不一致，使这些模型在工程设计中的应用受到限制。

本文作者以高温气冷堆用氦气干气密封为研究对象，在氦气介质条件下，建立流场、传热、变形及耦合过程模型，以及更接近干气密封实际力热状态的边界条件，发展干气密封的多场耦合分析方法，并对氦气干气密封在不同工况参数和槽型参数下的性能变化规律进行分析，为氦气干气密封的机制研究和结构设计提供一定的依据和帮助。

1 数学模型和数值方法

1.1 几何模型

文中以螺旋槽干气密封为研究对象，图1(a)为螺旋槽干气密封的端面结构示意图，各几何参数分别为：密封环带内径ri=96.5 mm，外径ro=120.3 mm，槽根半径rg=108.4 mm，平衡半径rb=100.5 mm，螺旋角α=15°，槽深hg=10 μm，槽数Ng=30，堰槽宽度比γ=1。密封介质为氦气。从干气密封结构特征出发，取旋转环与静止环密封环进行轴对称几何建模，如图1(b)所示，其他结构件的作用以边界条件体现。

图1 螺旋槽干气密封的端面结构示意图及固体域分析模型

1.2 流体域方程

流场分析采用螺旋槽机械密封的半解析式方法[8]。由MUIJDERMAN解析模型[13]，螺旋槽密封的压力分布如下：

(1)

式中：p(r)为气膜压力在密封环带上的径向分布；ωa为角速度；η为流体的动力黏度；ρ为流体密度；st为质量泄漏率；h1为槽区膜厚；h2为非槽区膜厚；C为积分常数。

干气密封运行时，密封环间隙中气体的密度、黏度及气膜厚度均会沿径向发生变化，MUIJDERMAN模型无法直接使用。为此，文中采用了基于该模型发展出的半解析式模型[8]。将密封端面流场在径向进行离散化处理，可用于介质黏度、介质密度以及膜厚沿径向位置发生变化等条件下的机械密封流场分析。考虑氦气介质的可压缩性，引入真实气体状态方程[14]：

(2)

式中：ρ为氦气的密度；p为气体的绝对压力；T为气体的绝对温度。

氦气的动力黏度可由黏温公式[14]确定：

η=3.674×10-7T0.7

(3)

在忽略气膜内沿径向和周向的热传导、忽略气体体积力条件下，流体域能量方程[15]如下：

(4)

式中：e为比内能；v为速度；p为压力；τ为黏性应力张量；k为热传导系数；S为体积源项。

当密封运行达到稳定状态时，补偿环的轴向工作位置可由力平衡方程确定：

Fo=Fc

(5)

式中：Fo为开启力；Fc为闭合力。

开启力可通过将气膜压力在密封环带上积分进行求解，其表达式为

(6)

闭合力由密封入口压力、出口压力与弹簧力共同组成，其表达式为

(7)

式中：po为密封入口压力；pi为密封出口压力；psp为弹簧比压。

1.3 固体域方程

旋转环和静止环中，均采用轴对称稳态热传导方程[16]：

(8)

用于求解热弹性固体域变形问题的控制方程[17]如下：

平衡方程为

(9)

式中：σrr、σθθ与σzz分别为径向、周向与轴向正应力；τrz为圆柱面上的切应力；br与bz分别为径向与轴向的单位体积力。

几何方程为

(10)

式中：εrr、εθθ与εzz分别为径向、周向与轴向正应变；γrz为径向与轴向间的切应变；ur与w分别为径向与轴向的位移分量。

物理方程为

(11)

式中：E为材料的弹性模量；ν为材料的泊松比；G=E/[2(1+ν)]为材料的切变模量；αT为材料的热膨胀系数；Tref为热应变的参考温度,文中取Tref=293.15 K。

1.4 边界条件

对于流体域，压力边界条件为

p(r=ri)=pi;p(r=ro)=po

(12)

温度边界条件为

T(r=ro)=To

(13)

式中：To为密封入口温度。

在流体域与固体域界面处，满足热流密度和温度连续条件：

Tg=Ts;qg=qs

(14)

式中：Tg和Ts分别为界面处气体与固体的温度；qg和qs分别为界面处气体与固体的法向热流密度。

对于固体域，图2(a)所示为固体域的受力与约束边界条件。其中密封环带部分施加由流场分析求解得出的密封气体压力分布，其余部分施加入口压力或出口压力，动环与静环背部靠近内侧处施加轴向位移约束。动环以一定的速度随轴旋转，因此对其施加离心力载荷。

图2(b)所示为固体域热边界条件。其中密封环带部分施加由流体域求解得出的热流密度分布；两密封环背部和底部与周围部件形成的间隙较小，气体流动速度慢，导热能力差，因此采用绝热边界条件；动环和静环外径处部分取对流换热边界条件(由于动环旋转引发的外径处气体流动会作用于静环，因此图2(b)中静环外径处与动环外径处所取对流换热系数相同)，即

(15)

图2 密封环边界条件示意

式中：h*为对流换热系数;ks为固体的热传导系数。

目前常用于机械密封对流换热系数求解的经验公式有BECKER[18]公式、DOANE[19]公式以及TACHIBANA[20]公式等。文中计算采用经典的BECKER公式，其表达式为

(16)

式中：kf为流体的热传导系数；D为旋转环相应边界处的直径；Re=ρuD/η为以D为特征长度的雷诺数，其中ρ为流体密度，η为流体的动力黏度，u为旋转环相应边界处的线速度；Pr=cpη/kf为普朗特数，其中cp为流体的定压比热容。

1.5 耦合过程和数值方法

干气密封运行中，流体场的压力分布与黏性剪切生热会导致密封环产生变形，进而改变气膜厚度分布；而气膜厚度分布的改变又会对流场与温度场产生影响。为准确反映干气密封的这种耦合效应，文中利用上述理论对干气密封的流场、传热、变形进行计算，并通过界面间的物理量传递对三者之间进行反复迭代。其中主要包含3个迭代过程：开启力闭合力平衡迭代，流体域固体域传热迭代，流体压力与密封环变形迭代。整体计算流程和3个迭代的嵌套关系如图3所示。

图3 耦合计算流程

在文中研究中，流体域方程采用半解析法求解；固体域传热方程和变形方程均采用有限元法进行数值离散和求解。

2 结果与讨论

基于文中建立的热-流固耦合模型，首先对氦气干气密封在典型工况下的基本工作特性进行分析；进而分别以转速和槽深作为工况参数和几何参数的代表，研究热-流固耦合条件下工况参数和几何参数变化对密封产生的影响；对所建立的模型的能力进行检验，并通过分析仿真结果深化对密封工作机制和特性的认识。

2.1 基本工况分析

文中计算设定的基本工况为：入口压力po为7 MPa；出口压力pi为0.2 MPa；弹簧比压psp为0.01 MPa；入口温度To为320.15 K，转速ω为4 000 r/min。旋转环采用碳化硅材料，静环采用石墨材料，其性能参数如表1所示。

表1 密封环材料性能参数

为了方便地评估密封环变形，以密封环带内径为参考点，定义密封环带上任意一点的轴向相对变形量为

Δz=z-zri

(17)

式中：z为密封环带上任意一点处的轴向位移分量；zri为密封环带内径处的轴向位移分量。

图4所示为密封环端面的轴向相对变形量分布情况。可以看出，动环端面最大轴向变形量为0.222 μm，其中热致变形沿径向分布使密封间隙趋于收敛形(这里指由外径向内密封间隙由大变小)，最大变形量为0.217 μm；力致变形沿径向呈曲线分布，与热致变形量综合作用后的总变形沿径向呈近似线性收敛分布。静环端面总变形总体上使密封间隙趋于收敛形，最大轴向变形量为1.607 μm，其中热致变形沿径向呈近似线性分布，最大变形量为0.748 μm；力致变形为曲线分布，导致总变形在密封环带内径附近有轻微的弯曲，其余部分为近似线性分布。由此可见，动环的力致变形量小于静环，这是因为碳化硅材料的弹性模量较大，因此动环抵抗压力变形的能力更强。热致变形会形成有利于密封稳定性的收敛锥度，并且在密封环端面的轴向变形中占据重要地位。

图4 密封环端面力致轴向相对变形量、热致轴向

图5所示为热-流固耦合模型与纯流场模型计算得到的气膜厚度沿径向分布情况。热-流固耦合模型中考虑了密封环变形对膜厚的影响，最小膜厚出现在内径处，大小为3.238 μm，最大膜厚出现在外径处，大小为5.066 μm；而纯流场模型由于无法考虑密封环变形的影响，因此气膜厚度沿径向保持不变，大小为2.960 μm。通过热-流固耦合模型计算得到的密封泄漏率为5.9 m3/h(标况下体积泄漏率，下同)；通过纯流场模型计算得到的密封泄漏率为4.1 m3/h。由2种模型的计算结果对比可知，密封环产生的变形会使得平均膜厚增加，从而增大密封气体的泄漏。

图5 2种模型的气膜厚度分布对比

图6所示为热-流固耦合模型计算得到的温度分布。可以看出，密封环端面温度自气体入口至出口先升高后降低，于内侧达到极大值。这是因为在气体强制对流换热的作用下，密封环外径处具有较好的散热性能，端面产生的热量可以及时传出，因此温度较低；距离入口越远的位置散热性能越差，因此温度不断升高。气体出口附近压力梯度较大，压力骤降使得气体产生明显的膨胀吸热现象，导致该区域温度下降。密封环端面远离入口的部分平均温度相对较高，热膨胀效应更加明显，因此密封环的热致变形呈收敛分布。碳化硅材料与石墨材料相比热传导系数较高，因此密封运转中产生的剪切热量主要由动环导出。

图6 密封环温度场分布

2.2 转速对密封性能的影响

以转速为代表量，研究了工况参数变化时密封性能和其他参数变化的情况。采用文中模型分别对密封转速为3 000、4 000、5 000、6 000、7 000、8 000 r/min时的工况进行仿真分析，计算过程中采用的其他参数同基本工况。

图7所示为不同转速下的气膜沿径向的温度分布。可以看出，转速为3 000 r/min时气膜的最大温升为2.68 K，而转速为8 000 r/min时最大温升达到了10.27 K。虽然转速升高后密封环外侧的对流换热条件会得到一定改善，但是气体黏性剪切效应更加明显，生热量急剧增大，导致气膜整体温度显著升高，可见黏性剪切生热是影响干气密封温升的主要因素。

图7 不同转速下气膜温度的径向分布

图8所示为不同转速下的动环与静环变形分布。

图8 不同转速下密封环端面轴向相对变形量

从图8可以看出，转速为3 000 r/min时动环的最大变形量为0.146 μm，静环的最大变形量为1.480 μm；而转速为8 000 r/min时动环的最大变形量为0.606 μm，静环的最大变形量为2.159 μm。这是因为在转速较高时密封温升较大，会造成更大的热致变形，从而导致动环与静环形成较大的收敛变形。

图9所示为不同转速下的气膜厚度分布。转速升高主要可能带来几个方面的效应：(1)螺旋槽产生的流体动压效应变大，导致气膜厚度增加；(2)螺旋槽对气体的泵送能力增强，也会使泄漏增大；(3)气膜剪切生热增大(不考虑膜厚变化)，使得温度和热变形发生变化；(4)压力分布变化，力致变形随之发生变化。从图9可以看出，气膜厚度整体随转速升高而增大，变形形成的锥度也同步增大，这是以上4种效应综合产生的结果。随之带来的密封泄漏和刚度的变化如图10所示，可见，泄漏率随转速升高显著增加，但气膜刚度有所下降，后者主要是因为膜厚增大导致的结果。

图9 不同转速下的气膜厚度分布

图10 不同转速下的泄漏率与气膜刚度

2.3 槽深对密封性能的影响

以槽深为代表量，研究了槽型参数变化时密封性能和其他参数变化的情况。采用文中模型在槽深为7～12 μm时，如2.1节所述的基本工况条件下进行了仿真分析。

图11所示为不同槽深下的气膜厚度分布。可见气膜厚度随槽深单调增加，这是因为槽深的增加会使流体动压效应、泵送能力、平均膜厚均增大(在文中选取参数范围内)。泄漏率对槽深变化十分敏感，随着槽深增加，密封泄漏率显著增大，如图12所示，但气膜刚度随转速增高而变小，其原因也是气膜厚度增加。

图11 不同槽深下的气膜厚度分布

图12 不同槽深下的泄漏率与气膜刚度

图13所示为不同槽深下的气膜沿径向的温度分布。可以看出，槽深为7 μm时气膜的最大温升为5.66 K，而槽深为12 μm时最大温升仅2.91 K。这是因为增加槽深会导致气膜整体厚度的增大，气体黏性剪切生热量随之降低，使得气膜温度下降；同时泄漏量的增加也会带走更多热量。

图13 不同槽深下气膜温度的径向分布

图14所示为不同槽深下的动环与静环变形分布。可以看出，动环变形随槽深增大逐渐减小，静环变形随槽深增大不断增大，但动、静环的整体变形量总体变化较小，对膜厚形状影响不大。

图14 不同槽深下密封环端面轴向相对变形量

3 结论

(1)考虑温度分布、密封环变形与气体性质、膜厚分布等因素之间的相互作用关系，建立了螺旋槽干气密封热-流固耦合模型，可用于泄漏率、气膜刚度、膜厚分布以及温度分布等密封性能参数的预测，可为干气密封结构设计提供有力的技术支撑。

(2)利用热-流固耦合模型与纯流场模型分别对典型工况下的密封性能进行了计算分析，结果表明，在文中所给参数范围内，密封环的热、力综合变形会明显改变间隙形状，从而对泄漏率等密封性能造成显著影响；而纯流场模型无法考虑密封环变形及温度梯度的影响，会导致计算结果产生较大误差。

(3)对不同转速及槽深下的密封性能进行了计算分析，结果表明，在文中所给参数范围内，转速及槽深的变化会直接对流场动压及传热性能产生影响，并通过热、流、固三场之间的相互耦合作用影响密封性能。