基于多变量摄动的超临界CO2干气密封动态特性

2021-05-15江锦波滕黎明孟祥铠李纪云彭旭东

化工学报 2021年4期

江锦波，滕黎明，孟祥铠，李纪云，彭旭东

（浙江工业大学机械工程学院，浙江杭州310014）

引言

超临界二氧化碳（supercritical carbon dioxide，sCO2）布雷顿循环发电技术因具有效率高和体积小等优势而有望成为未来的主流发电技术[1]。离心压缩机和向心透平等旋转设备是sCO2布雷顿循环系统的心脏设备，而轴端密封的低泄漏和高稳定是保证其高运行效率和低发电成本的关键所在[2]。根据美国桑迪亚国家实验室的推荐[3]，当功率等级大于1MW时建议采用控漏效果好的干气密封。

近十年来，sCO2旋转设备中轴承和密封等关键零部件的性能分析和结构设计已引起国内外学者的广泛关注。Thatte 等[2]和李志刚等[4]综述了sCO2旋转机械动密封的研究进展及技术挑战，包括流体膜多相态转换、复杂湍流和阻塞流、动力失稳和物性突变影响等。Conboy[5]、Qin 等[6]和Guenat 等[7]关于sCO2气体轴承的研究指出高速运转下的sCO2轴承中存在高度湍流化的流体膜，且流体物性与传统空气介质差异很大，在膜压求解中有必要考虑湍流效应和实际气体效应。Fairuz 等[8]、沈伟等[9]、许恒杰等[10]、严如奇等[11]和袁韬等[12]分析了sCO2螺旋槽干气密封中的离心惯性效应、湍流效应、阻塞效应或实际气体效应等实际效应对其开启力和泄漏率等稳态性能的影响规律。值得注意的是，目前关于sCO2气体轴承和干气密封的研究还主要停留于稳态性能，对其受到外界扰动下的动力学特性分析还鲜见报道。

摄动法是求解气体轴承和密封动态特性的一种常用线性化方法[13]，其假设浮动环受到微小的位移或速度扰动，从而使气膜压力产生微小变化，密封运动状态是在稳态基础上的叠加。刘雨川等[14]、Chen 等[15]采用摄动法分析了三自由度微扰下干气密封的刚度和阻尼系数，指出轴向微扰和角向微扰之间的交叉作用很弱；张树强等[16]和孟祥铠等[17]则进一步将摄动法用于动静压型干气密封和液体润滑机械密封的动力学特性分析中。上述研究主要是针对忽略湍流、惯性作用和实际气体效应的常规介质密封，而近年来各种实际流体效应和多相态共存对气膜密封和轴承动态特性的影响也引起了关注。宋鹏云等[18-19]和沈伟等[20]基于摄动法分析了实际气体效应、阻塞效应和惯性效应对干气密封动态特性的影响，李世聪等[21-22]基于微扰法分析了油气两相密封的动态特性，温建全[23]基于摄动法分析了sCO2径向箔片轴承的动态特性，并进一步提出一种位移速度增量法。实际上，正如Bi 等[24]在关于sCO2径向轴承动态特性研究中所指出的，密封受到的位移或速度扰动不仅会引起气膜压力的变化，同时还会引起包括密度、黏度在内的物性参数和Reynolds数的变化，并进一步引起湍流系数和惯性系数的变化，特别是在近临界区物性突变区，忽略这些变量摄动可能会给轴承和密封动态特性预测造成较大误差，而这正是目前关于干气密封动态特性分析中所忽视和欠缺的。由此可见，有必要在近临界区sCO2干气密封动态特性的研究中综合考虑各种实际流体效应及变量摄动的影响，以提高其预测精度。

本文在综合考虑离心惯性、湍流、阻塞和实际气体效应的sCO2干气密封膜压控制方程基础上，基于摄动法推导了包括压力、密度、黏度、Reynolds数、湍流系数和惯性系数的多变量摄动模型。在近临界区对比分析了三自由度扰动下sCO2和N2干气密封的动态特性，探讨了各种实际流体效应和变量摄动对sCO2干气密封刚度和阻尼系数的影响规律，获得了不同扰动条件下实际流体效应和摄动变量的主要影响因素，为sCO2干气密封气膜动态特性预测和后续干气密封稳定性、追随性分析提供依据。

1 sCO2干气密封动态特性数值模型

1.1 sCO2干气密封几何模型

图1 超临界CO2旋转设备干气密封结构示意图Fig.1 Schematic diagram of dry gas seal used in supercritical CO2 rotating equipments

相较于传统迷宫密封，干气密封通过减小泄漏损失和电机鼓风损失，有望使10MW 级系统提高2%～3%的循环效率[25]，故成为高功率等级sCO2动力设备的首选轴端密封型式。图1所示为sCO2布雷顿循环发电系统旋转设备的轴端干气密封结构。sCO2干气密封的主要结构包括动环及其安装座、静环及置于其背部的弹簧、推环和辅助密封等。在动环端面上开设有沿周向均匀分布的流体动压浅槽，其产生的气膜承载力与作用于静环背部的流体介质静压力和弹簧力平衡，进而可在动静环端面之间形成一层微米级厚度的稳定流体膜。

干气密封端面流体动压结构的形式多样，其中对数螺旋槽因能产生显著的流体动压效应和均匀的压力分布[26]，故在高速干气密封中应用最广。开槽密封环端面上设有槽区、阻流密封坝和周向阻流密封堰，高压腔和低压腔介质压力分别为po和pi。动压槽与密封坝的交界半径为rg，密封环端面的内径、外径分别为ri和ro，对数螺旋槽的螺旋角为β，开槽深度为hg，非开槽区气膜厚度为hb。

1.2 考虑离心惯性的变黏度变密度湍流雷诺方程

对于sCO2干气密封而言，其介质物性和工况条件的特殊性使其表现出异于常规N2干气密封的特性。图2所示为sCO2干气密封膜压求解中需考虑的实际效应及其影响。第一，其“黏度似气体、密度似液体”的高密度、低黏度特性使其容易达到很高的Reynolds数，离心惯性作用突显，且流态可能由层流转变为湍流，此时需考虑离心惯性作用和实际流态以修正干气密封膜压控制方程，可称之为“离心惯性效应”和“湍流效应”。第二，其密度、黏度等物性参数在临界点附近对于温度和压力的变化非常敏感，故以往基于等黏度假设和采用理想气体状态方程描述密度与压力、温度之间的关系不再适用，而需根据各点的压力和温度值通过查询REFPROP 物性数据库或采用更准确的物性方程以获得实际物性，这种实际气体物性的影响可称为“实际气体效应”。第三，sCO2干气密封的压力一般高于7.4 MPa，高压差特性使其靠近密封出口处的径向流速增大，甚至达到声速，使得以往将密封出口压力设定为大气压的强制压力边界不再适用[27]，而需以保证密封出口流速不超过当地声速为准则以改变密封出口的压力边界条件，这种影响称为“阻塞效应”。

根据文献[9]，极坐标下考虑离心惯性的变黏度变密度湍流雷诺方程为：

图2 超临界CO2干气密封实际流体效应示意图Fig.2 The real fluid effects in supercticail CO2 dry gas seal

式中，p、h、ρ 和μ 分别为密封端面任意点(r,θ)处的瞬态压力、膜厚、介质密度和黏度，Ω 为角速度。Gθ和Gr分别为与Reynolds数有关的周向和径向湍流系数，λ 为离心惯性项中与Reynolds 数有关的惯性系数，表示为：

式中，Re为Reynolds数，Re=ρvh/μ，其中速度v为周向流速vθ和径向流速vr的合速度；Gθ和Gr中的常数由Ng-Pan 湍流模型确定[24]。当其处于湍流状态时，α1=0.0136，β1=0.900，α2=0.0043，β2=0.960，而在层流状态下，上述常数均为零。湍流状态下惯性系数λ 中α3=0.885，β3=0.367，而层流状态下惯性系数λ=0.34[9]。

1.3 多变量摄动模型

类似地，瞬态压力p 也可由平衡位置的稳态压力p0和微扰压力pd叠加而得。根据微扰压力pd的泰勒展开，并忽略高阶项，瞬态压力p可表示为：

式中，实部压力和虚部压力分别为：

对于理想气体，密度与压力之间的关系通过理想气体状态方程获得，压力摄动可间接引起密度摄动，不过因其黏度设为定值，故不会存在黏度摄动，称之为经典变量摄动模型。对于考虑惯性和湍流效应的sCO2干气密封，不但需要考虑位移扰动所引起的密度、黏度等物性参数扰动，同时还需要考虑局部Reynolds 数摄动，以及由此引起的湍流系数和惯性系数摄动，如图3所示。

引入y 表示物性参数，将物性参数分解为稳态量（零阶）和摄动量（一阶）的叠加：

sCO2的密度、黏度等物性参数是关于压力和温度的二元函数。根据文献[8]中关于sCO2干气密封温度场的数值计算结果可知，在密封气膜的大部分区域温度基本保持不变，只有在靠近内径很小的局部区域会出现温度降低，故暂可将介质在密封间隙内的流动过程视为等温过程，此时密度和黏度摄动仅仅是关于摄动压力的动态响应，而无须考虑摄动温度的影响，如图4所示。

图3 干气密封动态特性求解中的变量摄动关联关系Fig.3 The relationship among different perturbed variables on the calculation of dynamic characteristics in dry gas seal

图4 sCO2介质物性参数扰动示意图Fig.4 Schematic of physical parameters perturbation of supercritical CO2

考虑到Re 中黏度μ 位于分母，根据幂函数的泰勒展开，式（10）可化为：

将式（4）、式（9）和式（11）代入Reynolds数中，并忽略所有二阶和高阶项，则瞬态Reynolds数为：

式中，λ0为平衡位置的稳态惯性系数。

1.4 摄动雷诺方程及动特性系数

将摄动形式的膜厚、压力、密度、黏度、湍流系数和惯性系数代入式（1）中，保留一阶项，忽略高阶项。分离零阶项和一阶项，分别得到稳态和摄动雷诺方程，其中稳态雷诺方程（17）用于求解稳态压力p0，摄动雷诺方程（18）用于求解摄动压力pd：

式中，各系数Ai和Bi可表示为：

将摄动雷诺方程（18）分别对Δz、Δα 和Δβ 求偏导，再利用复数量的实部与虚部分别相等，即可获得关于z、x和y三个方向的摄动雷诺方程组。

干气密封的刚度系数和阻尼系数用于表征密封气膜抵抗或抑制外界扰动的能力。当受到轴向和两个角向共计三个方向的扰动时，理论上来说共计有9个刚度系数和9个阻尼系数，其表达式为：

求解稳态雷诺方程和摄动雷诺方程组的压力边界条件包括密封内、外径压力边界和周期性压力边界。在密封外径ro处，p=po，pkj=0，沿着密封周向满足p0(θ+2π)=p0(θ)，pkj(θ+2π)=pkj(θ)，其中k=z, α, β，j=r,i。忽略阻塞效应时，密封内径的强制性压力边界条件为：当r=ri时，p=pi，pkj=0；当考虑出口阻塞效应时，密封内径处的阻塞压力边界条件为：

图5 sCO2干气密封动特性系数计算流程Fig.5 Solution process of dynamic characteristics calculation for sCO2 dry gas seal

式中，Ma1和Maexit分别为出口上游节点位置和出口位置的Mach 数，p1和pexit分别为出口上游节点位置和出口位置的压力，γ为气体绝热指数。

1.5 计算流程

图5所示为计算sCO2干气密封动态刚度和阻尼系数的计算流程。在给定密封几何参数、工况参数、摄动参数、物性参数和压力边界条件的基础上，采用有限差分法和超松弛迭代法求解稳态雷诺方程以获得稳态气膜压力分布，其中稳态气膜压力的收敛判据为：

式中，m 和n 分别为周向和径向网格节点数，本文中取为m=1201 和n=101 时可满足计算精度；t 为计算次数，ε为收敛残差，取为10-6。

进一步判断密封出口处的Mach 数是否大于1，如大于1则根据式（21）给定密封内径处的阻塞压力边界，否则直接输出压力分布及稳态密度、黏度、Reynolds 数、湍流系数和惯性系数的分布。求解摄动雷诺方程组以获得一阶实部压力pzr、pαr、pβr和虚部压力pzi、pαi和pβi，并根据式（19）和式（20）计算干气密封的刚度系数和阻尼系数。

密度、黏度等物性参数根据各计算节点压力从REFPROP物性数据库获得，而当忽略实际气体效应时，密度根据理想气体状态方程获得。当忽略湍流效应时，湍流系数Gθ=Gr=1；当忽略离心惯性效应时，惯性系数λ=0。当忽略阻塞效应时，密封出口处设为强制性压力边界，即其压力恒为pi。

2 结果讨论与分析

sCO2压缩机入口端干气密封的压力和温度处于临界点附近，此时介质的物性参数对于压力和温度的变化最为敏感，本文将sCO2干气密封的入口压力和温度设在临界点附近以研究其物性非线性变化对动态特性的影响。表1所示为本文计算sCO2干气密封动态特性时所采用的初始计算参数，未作特别说明，下文数值计算时采用表1所示参数。

表1 本文干气密封动态特性数值计算初始参数Table 1 Initial parameters adopted in the calculation of dynamic characteristics of dry gas seal

2.1 程序正确性验证

图6 考虑离心惯性和湍流的计算程序正确性验证Fig.6 Verification of calculated program with consideration of turbulent and inertia

通过与文献[5,28]的数值计算结果对比以验证本文考虑惯性效应和湍流效应计算程序和数值方法的正确性。图6(a)所示为考虑和忽略惯性效应条件下气体推力轴承的数值计算结果和文献值，在不同的轴承数下，本文计算结果与文献值误差始终小于2%，验证了本文考虑惯性效应计算方法和程序的正确性。图6(b)给出了不考虑惯性效应时湍流流态下sCO2气体推力轴承承载力与文献[5]计算结果。在转速较小时，本文的承载力计算值较文献值略高；随着转速的增大，本文承载力计算值与文献值差距逐渐缩小，当转速增大至n=75000 r·min-1时，本文承载力的计算值较文献值略高。不过在给定转速范围内，本文承载力的计算值与文献值的相对误差始终在15%以内，验证了本文中考虑湍流效应时计算程序与数值方法的正确性。

2.2 sCO2和N2介质干气密封动特性对比

为研究临界点附近sCO2特殊物性及其实际流体效应所引起的干气密封动态特性与常用氮气介质干气密封的差异，基于综合考虑四种实际效应的修正数值模型计算了不同频率比条件下sCO2、N2介质干气密封动特性系数和基于不考虑上述四种实际效应的经典数值模型（不考虑惯性和阻塞，基于层流和等黏度假设，采用理想气体状态方程获得密度）计算了N2介质干气密封的动特性系数。

图7所示为基于两种不同数值模型计算的sCO2和N2干气密封动特性系数随频率比的变化规律。当频率比较小时，刚度系数基本不变，当达到某一临界频率比时，刚度系数随频率比的增大而迅速增大，并在频率比接近100时又逐渐趋于平缓；阻尼系数在频率比较小时也基本不变，当达到临界频率比时，随着频率比的增大而单调递减。由此可见，sCO2与N2干气密封动特性系数随频率比的变化规律基本一致，不过不同数值模型和介质条件下的临界频率比略有差别。对于N2干气密封，相较于基于经典数值模型的计算结果，基于修正数值模型计算所得的阻尼系数在低频下显著提高，而在高频下略有降低。其原因可能在于考虑实际气体效应后的氮气密度减小，气膜吸收外界干扰能量的能力增强，抑制振动发散的能力增强[18]，不过在不同扰动条件下其阻尼特性差异的原因还有待进一步研究。可见，忽略各种实际流体效应也会给高速高压N2干气密封的动特性预测带来一定的误差。

较大的主刚度系数和阻尼系数有助于提高干气密封抵抗和抑制外界扰动的能力。从主刚度系数来看，低频下的sCO2干气密封较N2干气密封差异要高出20%，而在高频下N2干气密封则显著更高，甚至达到sCO2干气密封的2倍；从主阻尼系数来看，相比于N2干气密封，低频下sCO2干气密封略高，而在高频下则显著更小。这说明在低频扰动下sCO2干气密封的动态特性更优，而在高频扰动下动态特性则明显表现不佳，抵抗和抑制扰动的能力显著更弱，应通过合理的结构设计和系统调控以提高其运行稳定性。

图7 sCO2和N2介质干气密封动态特性系数对比Fig.7 Comparison on dynamic characteristics of dry gas seals lubricated with sCO2 and N2

2.3 实际流体效应影响

实际气体效应、离心惯性效应、湍流效应和阻塞效应是sCO2干气密封膜压求解中需考虑的四种典型实际流体效应[18-20]。基于多变量摄动数值模型，通过对比考虑不同实际效应条件下sCO2干气密封刚度和阻尼系数，以探讨各实际流体效应对sCO2干气密封气膜动态特性的影响，进一步获得不同工况和扰动条件下动特性的关键影响因素。

图8 所示为不同实际流体效应下sCO2干气密封刚度和阻尼系数随频率比的变化规律。在低频条件下，相较于综合考虑四种实际效应的计算结果，忽略湍流效应会使刚度系数和阻尼系数都显著降低，如其轴向主刚度系数kzz和主阻尼系数czz只有前者的80%和40%，角向主刚度系数甚至会出现负值，这是因为考虑湍流效应后流体的等效黏度增加，动压效应增强，能形成具有更大承载能力和刚度的气膜；只有在高频时的阻尼系数略有提高。忽略实际气体效应会使低频下的刚度和阻尼系数都下降，而使高频下的刚度和阻尼系数显著提高。可见，以往忽略湍流效应或实际气体效应的数值计算模型会显著低估低频条件下的动特性系数，而会过高估计高频条件下阻尼系数。总体而言，低频下湍流效应的影响最大，实际气体效应次之；而高频下实际气体效应影响最大，湍流效应次之。

根据式（19）可知，轴向主刚度系数kzz为实部压力pzr在密封端面积分的相反数，通过对不同实际效应下实部压力分布的分析，可解释实际流体效应对轴向刚度系数影响的内在机制。定义某一实际效应对应实部压力差值Δpzr为忽略该实际效应时实部压力相反数与综合考虑四种实际效应时实部压力相反数的差值。图9 所示为频率比Γ=1 时，忽略不同实际效应时对应的sCO2干气密封实部压力差值分布。忽略实际气体效应或湍流效应时，实部压力差值总体为负值，且最小值都出现在槽根附近，说明忽略上述两种实际效应会显著降低槽根附近的实部压力峰值，进而使刚度系数显著降低；相对而言，忽略实际气体效应对应的实部压力差值更小。忽略离心惯性效应和阻塞效应时，实部压力差值恒为正值，这是其轴向主刚度系数略有增大的原因。具体而言，忽略离心惯性效应会使槽根附近的实部压力略有提高，忽略阻塞效应则主要使靠近内径的局部区域实部压力略有提高，而对密封端面其他区域的影响较小。

图8 实际效应对不同频率比下sCO2密封动特性系数影响Fig.8 The influence of real effects on dynamic characteristics of sCO2 dry gas seal at different frequency ratio

图9 不同实际效应下sCO2干气密封一阶实部压力差值分布Fig.9 First-order real pressure differential distribution of sCO2 dry gas seal at different real effects

为进一步分析不同实际效应对sCO2干气密封动特性系数的影响，获得了不同线速度下忽略单一实际效应时干气密封的轴向刚度和阻尼系数，如图10 所示。随着线速度的增大，各实际效应对刚度系数的影响程度增大，而对阻尼系数的影响则无显著差别。湍流效应对刚度和阻尼系数的影响最为显著，如当v=200 m·s-1时，忽略湍流效应时的刚度和阻尼系数降幅超过60%。这是因为随着线速度的增大，密封间隙内流体的湍流化程度提高，忽略湍流效应时的流体黏度差异进一步增大。离心惯性效应对刚度系数的影响在超高速条件下也逐渐显现，如当v=200 m·s-1时会给刚度系数计算带来20%的误差，这是因流体惯性作用随速度增大而增强所致。值得注意的是，当速度很高时，除了忽略湍流效应会给阻尼系数带来很大误差外，忽略其他效应对阻尼系数的影响很弱。

2.4 变量摄动影响

在关于摄动法对干气密封和机械密封气膜动态特性的以往研究中，往往只考虑膜厚和膜压摄动，密度摄动是由压力摄动通过气体状态方程的间接作用实现，这种摄动处理可称为“经典变量摄动”。在此基础上，本文进一步考虑了黏度和Reynolds 数的摄动，以及由此引起的湍流系数和惯性系数摄动。在多变量摄动且综合考虑四种实际效应模型的基础上，通过忽略某一变量摄动以考察其对动特性系数的影响，并定义动特性系数增量比R 为某一种变量摄动模型对应的动特性系数与多变量摄动模型对应的动特性系数的相对增量比。

图11 为不同变量摄动模型对应的sCO2干气密封刚度和阻尼系数增量比。从图中可看出，在低频条件下，经典变量摄动和忽略湍流系数摄动对应的刚度系数约为多变量摄动模型的2 倍，同时也使阻尼系数计算值偏小5%～10%；忽略惯性系数摄动也会给刚度系数计算带来5%～40%的误差。可见，低频扰动下sCO2干气密封动特性系数计算中湍流系数摄动和惯性系数摄动不可忽略，黏度摄动影响不大。随着频率比的增大，当Γ＞1 时，经典变量摄动和忽略黏度、惯性系数和湍流系数摄动的动特性系数增量比总体都迅速减小。当Γ=100 时，相较于多变量摄动模型，所述四种变量摄动模型对刚度系数计算误差不超过1%，除经典变量摄动模型外对阻尼系数计算误差不超过5%。

图10 实际效应对不同线速度下sCO2密封动特性系数影响Fig.10 The influence of real effect on dynamic characteristics of sCO2 dry gas seal at different seal velocity

图11 变量摄动对不同频率比下sCO2密封动特性系数影响Fig.11 The influence of variable perturbation on dynamic characteristics of sCO2 dry gas seal at different frequency ratio

图12 变量摄动对不同线速度下sCO2密封动特性系数影响Fig.12 The influence of variable perturbation on dynamic characteristics of sCO2 dry gas seal at different seal velocity

图12所示为Γ=1时，不同线速度下上述几种变量摄动模型对应的动特性系数增量比。在不同线速度下，经典变量摄动和忽略湍流系数摄动对刚度和阻尼系数计算的影响最大，而忽略黏度摄动基本无影响。随着线速度增大，经典变量摄动和忽略湍流系数摄动对刚度系数的影响逐渐减弱，而对阻尼系数的影响增强，如当v=200 m·s-1时，会使轴向主阻尼系数偏小15%。忽略惯性系数摄动的影响也随线速度的增加而逐渐变强，不过都控制在5%以内。

3 结论

（1）针对sCO2干气密封高Reynolds 数流动及近临界区强物性非线性特点，综合考虑sCO2干气密封的离心惯性、阻塞、湍流和实际气体效应，基于摄动法获得包括膜压、密度、黏度、Reynolds数、湍流系数和惯性系数在内的多变量摄动干气密封动特性数值模型，显著提高sCO2干气密封气膜动态特性的预测精度。

（2）相较于常规N2干气密封，高频扰动下sCO2干气密封的动态特性不佳，轴向和角向主刚度、阻尼系数低于前者的50%。

（3）湍流效应和实际气体效应对sCO2干气密封动态特性影响显著，在计算干气密封动特性时不可忽略，且实际流体效应和变量摄动模型对低频和高频下干气密封动特性系数影响差异很大，低频下采用经典变量摄动和忽略湍流系数摄动会带来很大计算偏差，而高频下经典变量摄动模型对于刚度系数的计算精度较高，而对阻尼系数的计算则存在较大偏差。

符号说明

c——阻尼系数，N·m·rad-1

Gθ,Gr——分别为周向和径向湍流系数

h,hb——分别为瞬态气膜厚度和非开槽区气膜厚度，μm

hg——螺旋槽深度，μm

k——刚度系数，N·m-1

m,n——分别为周向和径向网格节点数

p——气膜压力，MPa

pi,po——分别为低压腔压力和高压腔压力，MPa

R——动特性系数增量比

ri,rg,ro——分别为密封内半径、槽底半径和外半径，mm

T——气体温度，K