超临界二氧化碳干气密封热-流-固耦合建模与变形特性分析

徐 洁 俞树荣 严如奇 丁俊华 丁雪兴

(兰州理工大学石油化工学院 甘肃兰州 730050)

以超临界二氧化碳(Supercritical carbon dioxide,S-CO2)为介质的布雷顿循环，相较于传统水蒸汽朗肯循环具有循环效率高、传热能力强、装置体积小等明显优势[1-2]，因此在舰船动力、航天航空、中小型核电领域极具应用潜力[3-6]。该循环主要包含透平主机、润滑系统、控制、密封系统等[7]。高性能的密封装置是宽负荷高强度工况条件下系统稳定运行的关键，由动、静2个相对非接触端面构成的干气密封已在各类高速旋转机械中广泛使用，因其使用寿命长、功率消耗小、工艺介质无泄漏等优点成为S-CO2循环系统的优选[8]。与传统的空气、氮气等干气密封相比，S-CO2干气密封介质物性呈强烈非线性变化，特别是在严苛运行环境中使得微间隙内流体流动复杂化，流场热力载荷及其他载荷的不稳定耦合作用引起密封环变形或者破裂，最终导致动静环碰磨失效，危害主机[9-10]。

现阶段，已有许多国内外学者发表了关于S-CO2干气密封的研究成果。2016年，FAIRUZ和JAHN等[11]利用CFD现代模拟仿真技术，研究了近临界点和远临界点两种S-CO2干气密封运行状况，结果表明近临界点处流体压力、温度、密度均变化较大，当选用较大密度的气体介质时不可忽略惯性力作用。2018年，许恒杰等[12]参考气体止推轴承的计算理论，在理想气体层流计算模型基础上，采用维里方程描述二氧化碳实际气体行为，通过考虑阻塞流效应和离心惯性效应对泵入式、泵出式2种干气密封进行研究，指出离心惯性效应对降低泵入式螺旋槽干气密封的泄漏率起积极作用，但对气膜开启力起消极作用，且当工况条件在二氧化碳临界点附近时，离心惯性效应愈加凸显。DU等[13-14]借助ANSYS CFX对比分析了空气和S-CO2为工作流体的干气密封，发现介质温度和压力更容易影响S-CO2模型的密封性能，同时表明应用S-CO2模型较空气模型更具优越性。2019年，沈伟等人[15]建立实际修正数学模型，对比讨论了高速工况下离心惯性效应、湍流扰动、实际气体效应对S-CO2干气密封关键密封性能参数的影响，指出各效应并非单一作用，而是强交互作用。李志刚等[10]综合分析了旋转机械动密封方向的国内外研究现状，表明在S-CO2的环境中应用时，泄漏量的控制、动静环碰磨导致的失效问题依旧是该领域面临的巨大困难。2020年初，DU等[16]针对S-CO2轴流式涡轮建立4种密封模型，指出与迷宫密封相比，干气密封的泄漏降低了99.38%，功率效率和等熵效率分别提高了0.88%和1.56%;同时，通过研究发现，密封端面会发生反向轴向变形，且相较于机械载荷对密封端面变形的影响，热载荷的影响更加显著。

综上所述，针对S-CO2干气密封流场热动力学的研究已取得很大进展，但流场与结构变形耦合机制的相关研究较少。本文作者选用经典螺旋槽干气密封结构，通过SolidWorks建立计算域模型，借助CFD和CSM计算机仿真技术，研究S-CO2干气密封动、静环在多重载荷共同作用下的变形规律。文中研究可加快S-CO2干气密封结构的研发进程，有利于降低试验成本，为S-CO2干气密封的设计与使用提供参考。

1 计算模型建立

通过激光打标技术，在密封动环表面雕刻螺旋槽，其结构如图1所示。气体在螺旋槽的动压效应下从外径高压侧泵送至内径低压侧，在摩擦副间形成一层2～6 μm且具有较大刚度的润滑气膜，从而实现零泄漏。

图1 动环结构示意

其中，槽型满足对数螺旋线方程，如公式(1)所示。

(1)

式中：R为密封动环端面任意处的半径，mm；Rg为螺旋槽槽根半径，mm；θ为展开角度，rad；α为螺旋角的余角，(°)。

选择如表1所示的端面结构尺寸参数，建立流体域和固体域的三维实体模型如图2所示，并将流体域计算模型导入ANSYS ICEM中进行网格划分。

表1 螺旋槽干气密封的几何参数

图2 三维实体模型

2 数值模拟过程

2.1 边界条件2.1.1 流体域边界

选取一个周期的流体域进行计算时，设置流体计算域的周期性边界为periodic，其中，与密封动环端面接触的面旋转，初始转速nr为15 000 r/min，与密封静环端面接触的面静止。在外径处设置入口压力pi为10 MPa，入口温度为Ti，其值为393.15 K。内径处设置出口压力p0为0.101 325 MPa。选择能量方程，开启黏性生热项。

2.1.2 密封环热力边界

腔体、推环、轴套、副密封O形环、弹簧、密封环组件等共同组成了S-CO2干气密封。动环围绕转轴中心轴线旋转，背部与动环座接触，因此，在动环内圆孔面建立柱面坐标系，设置动环旋转速度为nr，并选择周期性边界为Cyclic Symmetry。动环外径受介质压力pi。动环刻槽端面为传递流场压力载荷p和温度载荷T的流固耦合面。动环内径压力为标准大气压p0，周向施加约束，其余面视为绝热状态。密封动环边界条件设置如图3所示。

图3 密封动环边界条件

αm为动环外径壁面与周围介质的对流换热系数[17]，可由式(2)计算

(2)

式中:Nu为努塞尔数；Pr为密封介质的普朗特数；μ为介质动力黏度，Pa·s；Cp为介质定压比热容，J/(kg·K)；λ为介质导热系数，W/(m·℃)；Rec和Rea分别用于评价介质的绕流效应和横向绕流效应；ω为动环旋转角速度，rad/s；vm为密封动环外半径壁面处介质的平均轴向速度，m/s；ν为介质的运动黏度，m2/s。

施加在静环背部的弹簧力Fs，可看作在静环背部沿环形均匀分布的面载荷。设置周期性边界为Cyclic Symmetry。静环与气膜的接触端面受到气膜压力载荷p和温度载荷T的作用。尽管静环在弹簧的作用下可轴向移动，但在稳态运行时，静环的轴向位移近似为0，因此在轴向施加固定约束。此外，静环在静环座内不发生转动，故而对静环施加周向固定约束，密封动环边界条件设置如图4所示。

图4 密封静环边界条件

αs为静环外半径壁面与外界的对流换热系数，可由式(3)[17]计算

(3)

式中:δ为密封腔与静环的径向间隙，mm；ε为修正系数，取值一般为2；vs为密封静环外半径壁面处介质的平均轴向速度，m/s。

αl为密封环端面与气膜的对流换热系数，可由式(4)[18-19]计算

(4)

式中:Rm为密封环的中径，mm。

当气膜开启力等于闭合力时，系统达到平衡状态形成稳定的密封间隙，而闭合力Fc由弹簧力Fs和静环背侧所受气体压力共同产生。因此，静环背部所受弹簧力可由式(5)[20]计算

(5)

式中:Fs为弹簧力，N；B为平衡比，取值0.8；F为开启力，N。

2.2 流体域计算模型

以入口的工况条件确定介质物性，并针对表面光滑度较差且开设造型的端面密封结构，计算雷诺数Re判断流动状态。若雷诺数Re大于临界雷诺数Rec(Rec=500)时，采用RNGk-ε湍流模型计算。其中，雷诺数Re由公式(6)[21]估算得到。

(6)

式中:ρ为介质密度，kg/m3。

2.3 求解设置

密封介质选用CO2，在Fluent中使用NIST Real Gas材料模型，即可保证流场中每一点为真实气体，并选用压力和速度为Coupled Solver激活该模型。速度采用无滑移边界条件，空间离散方式全为二阶精度迎风格式，迭代精度均设为10-6。

2.4 固体域材料

S-CO2干气密封动环材料选用WC-Ni，静环材料选用SiC，具体材料物性参数如表2所示[22-23]。

表2 密封环物理及力学性能

2.5 热流传递方法

由于流场压力载荷、热载荷或离心力等多重物理场叠加导致密封环端面发生不同变形。在热流数据单向传递中，忽略密封环结构变形对密封间隙内流场的影响，只将流场计算结果单向传递给密封环进行结构分析计算，可在ANSYS Workbench中通过连接模块实现。

3 计算结果与分析

3.1 计算方法有效性验证

采用以上计算方法，选择文献[11]中的结构参数，设置介质压力8.5 MPa，出口压力0.101 3 MPa，进口温度和壁面温度为740 K，转速30 000 r/min。经过计算得到沿半径分布的平均压力如图5所示。对比发现，文中计算结果与文献计算结果一致性较好，但随着压力的增大，二者相对误差增大，最大误差发生在槽根径处，不超过9%，验证了文中计算方法的可行性。

图5 计算方法压力验证

3.2 流场计算结果分析

设置表1中的结构参数和前文工况参数，计算得到超临界二氧化碳干气密封摩擦副间流场的压力和温度分布，如图6所示。

图6 流场计算结果

由图6(a)可以看出，微间隙内流场压力从外径逐步向内径递减，槽根处压力只有略微上升，表明动压效果并不明显，出口处达到压力最小值。由图6(b)可以看出，流场温度沿半径减小方向而降低，在内径出口处出现最低温度，其值为310.88 K。同时发现，槽根至内径区域相较于外径至槽根区域的温度变化率大，即螺旋槽的存在使流场温度下降延缓，这主要是因为螺旋槽的存在产生了较多的黏性耗散热，抵消了一部分因CO2膨胀做功引起的温降。

3.3 密封环变形计算结果分析

导入流场数据并设置动环相关参数，经过求解获得不同影响因素下的动环轴向变形分布如图7所示，静环轴向变形分布如图8所示。

图7 动环轴向变形分布

图8 静环轴向变形分布

观察图7(a)和图8(a)发现，动、静环在只受到温度的作用时膨胀变形，最大轴向热变形均发生在密封环耦合面的外径处，最小轴向热变形均发生在密封环背面的外径处。由于CO2膨胀做功使得外径处入口温度大于内径处出口温度，且线速度在外径处最大产生较多的摩擦热[24]，因此密封环轴向热变形沿半径减小的方向减小。密封环轴向热变形沿厚度方向发生均匀阶梯递减变化，即密封环的耦合面发生最大轴向热变形，逐渐向密封环背侧递减，表明密封环的热变形主要受到微间隙流场的温度影响。同时可以看出，静环较动环发生的热变形小，这是由于动环选材为WC-Ni，它的热膨胀系数较静环材料SiC的大，因此动环对于热载荷较敏感，产生较大热变形。

由图7(b)可以看出，当只受到力载荷的作用时，动环被压缩发生反方向的变形，且变形量沿半径减小的方向减小，最大力变形出现在密封动环背侧。由图8(b)可以看出，静环的最大轴向力变形主要分布在2个地方，一处为静环内径壁面，另一处为静环背侧，二者均因防止静环周向旋转施加刚性约束所致。同时可以看出，静环较动环发生较大力变形，这是由于静环材料SiC的弹性模量较动环材料WC-Ni的小，因此静环抵抗变形的能力较差，受到力载荷的影响较大。

观察图7(c)和图8(c)发现，动、静密封环热-流-固耦合轴向变形的分布与轴向热变形的分布类似，且数量级与热变形为同一数量级，远大于力载荷引起的变形量，可见，密封环的热-流-固耦合变形中热载荷起主导作用。对比可知，动、静环的热-流-固耦合变形均大于热变形和力变形，表明力载荷引起的压缩变形在一定程度上抵消了热载荷引起的膨胀变形。

3.4 工况条件的影响3.4.1 转速的影响

设置动环旋转速度变化范围为5 000～15 000 r/min，保持其他参数不变，分析转速对密封环变形的影响情况，如图9所示。

图9 转速对密封环变形的影响

由图9(a)和图9(b)可知，随着转速的增大，动、静密封环的最大轴向热变形和力变形增大。转速增大，速度梯度在气膜厚度方向上增大，因而流体产生更多黏性剪切热，流场温度上升，致使热变形随之增大。转速增大，动压效应增强，造成螺旋槽区域应力增大，气膜开启力增大，因此力变形随之增大。由于静环只受到来自流场的热力作用，而动环因自身旋转还受到惯性力的作用，因此随转速的增大动环的力变形增幅明显，而静环增幅极小。由于动环的力变形增大明显，抵消热载荷膨胀变形的程度加强，最大轴向热-流-固耦合变形从5.12 μm减小至5.10 μm，而静环的力变形增幅极小，不足以减弱热变形的增大趋势，最大轴向热-流-固耦合变形从4.76 μm增大至4.81 μm。

3.4.2 介质压力的影响

设置介质压力变化范围为6～10 MPa，保持其他参数不变，分析介质压力对密封环变形的影响，如图10所示。

图10 介质压力对密封环变形的影响

由图10(a)和图10(b)可知，随着介质压力的增大，动、静密封环的最大轴向热变形均呈下降趋势，这是由于介质压力越大，进出口压差引起的膨胀温降越大[25]，流场温度降低，进而导致热变形减小。随着介质压力的增大，动、静密封环的最大轴向力变形近似线性增大，动环最大轴向力变形为0.116 μm，静环最大轴向力变形为0.138 μm，因为在转速不变的情况下，介质压力增大，气膜压力增大，介质压力和气膜压力二者共同作用于动、静密封环产生的压应力增大，因此动、静环的力变形随之增大。静环材料抗压能力小于动环材料，因而静环最大轴向力变形的增长率为66.25%，远大于动环最大轴向力变形的增长率6.18%。随着介质压力的增大，动、静密封环的热-流-固耦合变形均呈下降趋势，可见介质压力在该范围内增大引起的力变形不足以改变热变形的影响，热变形依旧起主导作用。

3.4.3 进口温度的影响

设置进口温度变化范围为393.15～433.15 K，保持其他参数不变，分析进口温度对密封环变形的影响，如图11所示。

图11 进口温度对密封环变形的影响

由图11(a)和图11(b)可以看出，动、静密封环的最大轴向热变形和最大轴向热-流-固耦合变形，均随进口温度的增大呈近似线性增大的方式变化。在其他条件不变的情况下，进口温度升高，流场温度上升，密封环热变形增大。当进口温度从393.15上升至433.15 K时，密封动环的最大轴向热变形增大了40.79%，密封静环的最大轴向热变形增大了34.90%，可见密封动环对热载荷的变化较为敏感。流场温度的上升使得CO2黏度增大，密度减小，气膜开启力减小，因此动、静密封环最大力变形减小。但由于文中所选进口温度变化范围较小，因而最大轴向力变形量减小幅度很小，动环最大轴向力变形保持在0.115 μm附近，静环最大轴向力变形保持在0.137 μm附近。结合表2对比可知，静环材料SiC的热膨胀系数和弹性模量均小于动环材料WC-Ni，所以静环在相同的热载荷下不容易发生正向热变形，在相同的力载荷下容易发生反向力变形，因此，随着进口温度的上升，静环的力变形抵消热变形的程度加大，热-流-固耦合变形增大趋势随之变缓。

4 结论

(1)动、静密封环在热载荷作用下发生膨胀变形，在力载荷作用下发生压缩变形；密封环在热-流-固耦合条件下，热变形起主导作用。

(2)密封环轴向热变形沿半径方向减小，沿厚度方向发生均匀阶梯递减变化；密封动环的轴向力变形沿半径的减小方向增大，最大力变形出现在密封动环背侧；静环最大轴向力变形出现在内径壁面和静环背侧。

(3)动环最大轴向热变形随着转速和进口温度的增大而增大，随着介质压力的增大而减小；动环最大轴向力变形随着转速和介质压力的增大而增大，随着进口温度的增大基本不发生变化；动环最大轴向热-流-固耦合变形随着转速和介质压力进口温度的增大而减小，随着进口温度的增大而增大。

(4)静环最大轴向热变形随着转速和进口温度的增大而增大，随着介质压力的增大而减小；静环最大轴向力变形随着转速和介质压力的增大而增大，随着进口温度的增大基本不发生变化；静环最大轴向热-流-固耦合变形随着介质压力的增大而减小，随着转速和进口温度的增大而增大。