超临界二氧化碳涡轮升速过渡中密封-转子系统的动力学行为

2024-03-01司和勇王瑶俐曹丽华陈东超

航空学报 2024年2期

司和勇，王瑶俐，曹丽华，陈东超

东北电力大学能源与动力工程学院，吉林 132012

在提升超临界二氧化碳涡轮系统效率和容量的目标下，超临界二氧化碳的做功参数不断提高［1］，会不可避免地出现流体诱导振动的问题。尤其在动、静止部分的径向密封间隙内，受不平衡质量以及加工偏差的影响，密封间隙多处于不均匀分布状态，而高转速的旋流往往会形成较强的气体动压效应［2］。同时，转子微小涡动形成的流体涌动和密封齿节流形成的高低压涡流区会加剧密封气动作用，改变密封-转子系统的动力学特性，威胁涡轮运行的稳定性［3-5］。

对于密封所形成的气动性能研究多以动力特性分析为主，即通过密封实验或流场的数值求解获得流动参数，从而计算密封的静态力或动态特性。理论数值求解方面，密封腔室的控制体积模型和振荡流体力学模型均能有效反映密封内部的流动特征［6-7］。在实验方面，现有关于超临界二氧化碳的涡轮实验研究注重热力系统循环的分析较多，而在密封气动性能以及流体诱导的振动方面报道较少［8-9］。从传统涡轮机械的实验结果来看，流体形成的不平衡力多来源于转子偏心所形成的不均匀压力分布［10］。张万福等［11］针对密封的动静态特性开展相关研究，实验结果表明密封内会形成强烈的气动作用，不利于转子系统的稳定运行。目前，计算流体动力学以其高精度的数值求解和可视化功能广泛应用于密封流动研究。在密封数值仿真模型方面，学者们也不断追求建立与实际更为贴近的密封-转子运动模型，从而获得更加准确、具有实际参考性的结果。基于实验平台建立的静偏心模型虽然只考虑了转子的偏心和自转作用，但在密封静态特性分析中具有较好的准确性和可用性。丁学俊等［12］以600 MW 汽轮机隔板密封为例，分析了进出口压比对密封泄漏量和密封激振力的影响。由静偏心模型衍生的相对旋转模型可将动态涡动转化为静态流场的求解［13］，但从密封流场周向速度分布来看，其与实际情况相差较大，同时密封动力系数求解的准确性也尚未被实验证实。近十年来，国内西安交通大学李军等［14］率先提出了多频涡动模型，该模型基于小扰动理论和叠加扰动频率实现密封动力特性的精确求解，仿真结果与密封实验结果具有较好的一致性。在此基础上，李志刚等［15-16］对不同涡动轨迹、进口预旋等因素进行了分析。张万福等［17-18］在采用多频涡动模型的同时，建立了密封动力系数的差分求解方法，实现了转子大偏心涡动时密封的动力特性分析，为转子大幅涡动的密封动态仿真提供了有效手段。从密封内部流动和密封气动性能的形成机制来看，密封进出口压比、进口预旋以及密封结构等参数均会对内部流动产生影响，改变密封的动态特性。压比和正预旋增强均会导致密封动力系数增大［19］。通过改进密封结构从而削弱密封的气动作用是行之有效的，由此学者和工程师们设计出了蜂窝密封、螺旋篦齿密封、反向预旋板密封、刷式密封以及扇贝式等阻尼密封，并进行了相关研究。阻尼密封的抑振机理主要是通过降低工质周向旋流强度来减弱密封气动作用。Zhang 等［20］建立了篦齿-刷式密封模型进行分析，该密封不仅可以减少泄漏量，还可以利用刷丝的安装位置来控制密封的气动作用。鉴于超临界二氧化碳具有低黏度、高能量密度等特征，尹露等［21］开展了超临界二氧化碳高低齿密封的泄漏流动和动力特性的研究，并从结构角度进行了动力学特性的优化。孙丹等［22］对实际气体参数影响密封动力特性开展研究，结果表明密封的动力特性与工质的摩尔质量在低频时呈正相关变化，在高频时呈负相关变化。Si 等［23］利用多频平行涡动模型对比了蒸汽与超临界二氧化碳密封气动作用的显著程度，发现超临界二氧化碳的气动作用会在高频范围内形成较大的有效阻尼，在低频范围内形成较小的有效阻尼，而且密封的气动作用不可忽略。从转子动力学角度来看，转子裂纹引起的结构刚度变化、不平衡质量的存在以及流体激振均会使转子系统失稳。因此在判断转子运动稳定性时，应从系统的转子动力学分析入手。在转子动力学研究方面，转子系统的不确定因素分析及数值分析方法是当前研究的主要内容［24］。Jia 等［25］基于简单Jeffcott 模型建立了非概率凸模型来描述不确定参数的影响。通过嵌入切比雪夫展开函数提出评价方法，结果表明支承刚度和阻尼对转子系统的影响要大于质量和密度等材料参数。密封的气动作用与轴承类似，其对转子的影响最终也是以刚度和阻尼的形式改变了系统的动力特性，因此学者们将密封耦合到转子动力学分析中。在密封-转子动力学研究方面，Li 等［26］基于哈密顿原理提出了汽轮机系统的转子/轴承/密封系统非线性模型，采用Musz⁃nyska 模型和非稳态轴承油膜力模型来描述非线性蒸汽激振力和油膜力。国内学者瓮雷［27］也基于该模型分析了转子裂纹的动力响应。曹丽华等［28-29］将密封激振力拟合成关于频率和偏心的函数带入到转子运动方程中，实现密封与转子的耦合分析，但其拟合公式仍存在一定偏差。对于涡轮机械来说，密封转子系统是自成一体的耦合系统，密封气动性能会影响转子的运动，而转子运动变化后又导致密封内部流动改变，从而形成新的气动作用，二者互相干涉［30］。

通过不同频率叠加建立的多频涡动模型可以得到不同涡动频率的动力系数，但是转子弹性刚度、不平衡质量和密封气动作用对转子自由运动的影响却无法考虑，而且这种耦合影响极有可能使密封动力系数在设定频率点以外产生显著波动。在密封-转子系统的动力学分析方面，密封的气动作用多以拟合公式或经验参数公式嵌入转子运动方程中，主要考虑了密封作用对转子的影响，转子运动对密封气动特性的影响却被忽略或者存在一定偏差。因此，关于密封-转子系统的动力学特性研究需要综合考虑各种因素，有必要建立一个全面考虑转子结构刚度、不平衡质量以及密封气动作用的双向耦合动力学分析模型，开展更加深入的分析。从超临界二氧化碳涡轮实际运行的情况来看，不同的自转速度会使密封内部的周向流动强度变化，而在升速过程中，超临界二氧化碳涡轮常处于宽范围的变速过程，这也使得密封的气动性能处于大范围的变化，即动力特性大幅度改变，极易诱导转子失稳。通过对不同自转速度工况仿真可以获得每个自转速度下的密封动力系数，但升速过程中，转子由低速过渡到高速时的动态过程无法展现。

针对密封-转子系统非线性运动以及升速动态过程中密封动力特性模糊的问题，本文基于两层密封-转子运动控制方法分别推导转子自转运动与非线性涡动的耦合方程，建立了转子自由涡动的非线性运动模型。在此基础上，改进自转速度连续变化仿真模型，实现动力学参数随时间变化的动态求解。考虑转子弹性刚度、不平衡质量和密封气动作用与转子自由运动的双向耦合影响，实现升转速过程中密封-转子系统的动力学行为分析。

1 计算模型

1.1 密封物理模型及边界

超临界二氧化碳涡轮内多采用径向密封与轴向干气密封，其中径向密封是引起转子横向振动动力特性变化的主要原因，因此本文以文献［30-31］中1.5 级超临界二氧化碳涡轮的径向平齿密封为例，建立的三维物理模型如图1 所示。当转子发生涡动运动时，超临界二氧化碳会因转子偏心和自转作用在密封腔室以及密封间隙内形成周向流动，从而产生气体动压作用。为清晰展示密封气动作用与转子运动的耦合影响，该模型只考虑密封内流域特性变化以及密封-转子系统的动静边界。密封进出口流域与叶栅通流部分衔接，以充分考虑进口气流的实际流动特性。流体域采用结构网格划分，并在近转子壁面附近进行网格加密，以保证Y+值在合理的范围内，具体网格参数可见参考文献［23］。为获得与实际相符的流场特性，密封进出口采用压力边界，二氧化碳工质为可压的实际气体属性，从Fluent 的NIST REFPROP 数据库调用，边界参数根据额定工况数据设定。密封结构参数、仿真计算边界参数以及时间耦合参数如表1 所示。转子的升速过程是具有加速度的运动，计算时间步长取1×10−4s 以便准确计算转子加速运动时的轨迹，且该步长要与Runge-Kutta 的差分步长一致。在本研究中，时间序列只作为采样的依据，根据采样点将转子加速过程离散，因此转子的加速过程可以无限制延长，从而保证本文研究具有现实意义。

表1 密封结构参数与仿真边界参数Table 1 Seal structure parameters and simulation boundary parameters

图1 密封物理模型及网格Fig.1 Seal physical model and mesh

1.2 升转速的非线性运动方程

根据质量集中法建立超临界二氧化碳涡轮转子运动模型，考虑密封-转子系统中转子的弹性刚度、不平衡质量离心力和密封气动作用。其中轴承油膜力通过转子自身刚度来维持转子平衡，将其简化为铰支边界。该运动模型可近似为双自由度的弹簧-阻尼-质量系统，如图2 所示。

图2 密封-转子动力学模型Fig.2 Seal-rotor dynamic model

其动力学基本控制方程为

进一步扩展为

式中：Δc、Δk分别为密封气动作用所形成的附加阻尼和刚度。由于密封所形成的刚度和阻尼难以直接数值求解，可借助流体力学计算软件Flu⁃ent 求解密封流场，将获得的压力积分可得到密封作用力Fz、Fy。再以密封力的形式耦合到运动方程中，可得到考虑密封气动作用的转子运动方程。由密封气动作用形成的密封力与刚度和阻尼的关系为

将式（4）分别代入式（2）和式（3）可得

式中：M为模化质量；Fz、Fy为z和y方向的密封分力；e为不平衡质量偏心；ω为转子旋转速度；z、y分别为z、y方向上的位移；C为阻尼。

通过上述转换可避免直接求解密封刚度和阻尼，将密封气动作用以密封力Fz、Fy形式施加到转子上，而密封力Fz、Fy可通过UDF 数值传递在流场求解中直接获取。后续的密封动力特性分析则基于仿真结果所得到的密封力和转子位移，再通过式（4）进行求解。与以往研究不同，在转子变转速运动过程中，式（5）与式（6）中的加速度是未知的，而且整体运动呈非线性变化。因此，采用四阶Runge-Kutta 法求解运动微分方程，且差分步长为0.000 1。具体计算流程如图3所示。

图3 耦合计算流程Fig.3 Coupling calculation process

首先流场仿真可获得流场参数，通过UDF编程压力积分直接获得密封力，即式（5）与式（6）中的密封力Fz、Fy，然后在UDF 程序中求解式（5）和式（6），获得转子位移与速度，再将位移与速度施加到流场仿真中，通过动网格技术驱动转子运动，当前计算完成后进入下一步迭代计算，实现流场与动力学的双向耦合求解。在密封流场仿真中连续性方程收敛残差为10−4，动力学方程求解收敛精度为10−4。为保证求解数值的稳定性，不同阶段的转速工况需维持2 倍以上的周期运动后再进入下一阶段的升速过程。各阶段数据采样分布如图4 所示。

图4 转速变化过程分布Fig.4 Distribution of speed variation process

本文采用的数值仿真方法已在文献［23］中与实验结果进行了对比验证，结果表明该方法具有较好的准确性。研究采用的非线性转子涡动模型是基于多个包含工作频率时序叠加的方程，对于工作频率下非线性涡动模型的准确性已在文献［30］中得到了验证。动力系数的求解方法采用文献［30］中的微元差分法。

与前期研究不同，本文建立的运动方程能够使自转速度连续变化，实现动力学参数随时间变化的动态求解，即升转速的过程。而这一过程与稳定自转速度的转子运动具有截然不同的动力学特性，其动力学特性既是频率的函数也是时间的函数。因此，为进一步验证本文模型的准确性，将随时间变转速的非线性涡动模型与定转速同涡动中心的多频涡动模型对比进行分析，具体思路及结果如下。

提取非线性涡动模型中的工频动力系数可获得不同自转速度下的动力系数变化。从图5 可以发现，与同涡动中心的多频涡动模型对比，非线性涡动下的密封动力系数与多频涡动的结果高度相似，二者存在微小偏差的主要原因在于非线性涡动模型考虑了转子弹性刚度和不平衡离心作用，在验证中已尽可能保证二者涡动中心点一致，但这种偏差是实际中必然存在的一种扰动因素，这也是本文所建立模型的优势，即可以考虑转子弹性刚度和不平衡离心产生的影响。而多频涡动模型是基于原点对称的运动控制方程，无转子弹性刚度和不平衡离心作用，其自转速度与涡动频率的相关性较差。从实际情况来看，采用非线性涡动模型更具有实际意义。其中平均直接阻尼为z与y方向直接阻尼的平均值，表征流体对转子运动的阻滞作用。平均直接刚度为z与y方向直接刚度的平均值，表征密封力对转子运动的促进作用。有效阻尼则表征系统受交叉耦合作用和阻滞作用的综合结果，有效阻尼越高表示系统趋于稳定的能力越好。平均直接阻尼、平均直接刚度和有效阻尼的定义分别为

图5 动力系数验证Fig.5 Validation of dynamic coefficients

式中：kzz、kyy为直接刚度；kzy、kyz为交叉耦合刚度；czz、cyy为直接阻尼；czy、cyz交叉耦合阻尼。

2 结果分析

2.1 频域动力系数

计算所得转子升速过程中的轴心涡动轨迹如图6 所示。在T1 转速阶段，受转子弹性刚度作用，轴心从初始扰动位置逐渐沿竖直方向运动在平衡位置形成稳定涡动。随着自转速度的增加，转子涡动运动向进动方向偏移，涡动半径不断扩大。这是由于自转速度的升高导致不平衡质量的离心作用增强，从而引起涡动范围的扩大。随着速度的升高，转子涡动中心的位置在竖直方向呈现先下降后升高的变化。该系统中转子的重力主要依靠转子弹性刚度平衡（轴承传递到转子），在T1～T4 升速阶段，转子的涡动中心处于右斜下方，且涡动速度和涡动半径相对较小，密封形成的气动作用推动转子继续运动。在T5～T7 升速阶段，转子涡动速度和涡动半径相对较大，形成了类似气体轴承的支撑作用，促使转子涡动中心向上偏移。

图6 轴心涡动轨迹Fig.6 Whirling track of rotor axis

在升速过程中，转子自转速度随无量纲时间变化，属于时序的涡动频率叠加。采用传统的快速傅里叶变换无法获得全部转速的对应频率。本文采用增补的方式提取每个阶段的转子涡动数据，将各阶段内的涡动数据进行时均化，可获得密封动力特性。时间微元平均方法如下：

对数据点进行增补，形成周期性循环采样点，再以每个稳定运行阶段的时间为基础对其进行时间微元平均化：

基于式（10）和式（11）的时间平均方法，可得到采样周期内时间微元平均的密封力与位移，最后通过文献［30］中的微元差分理论求解动力系数。

图7 为时间平均下的密封直接刚度和交叉刚度。直接刚度的绝对值随涡动频率的增加而增大，表明转速升高时密封的气动作用逐渐增强。交叉刚度随频率的增大先增加后减小，这是由于转子速度升高过程中，在T1～T4（对应频率50～200 Hz）阶段内气动作用促使转子涡动中心向右斜下方偏移，竖直方向和横向的动静间隙偏差增加，所以竖直交叉作用先增强。在T5～T7（对应频率250～350 Hz）阶段内气动作用促使转子涡动中心向上方偏移，竖直方向的动静间隙偏差减小，所以竖直交叉作用减弱，但由于转子持续的横向偏移，引起交叉刚度kzy的方向改变，由正方向改为负方向。

图7 时间微元平均的直接刚度和交叉刚度Fig.7 Time-averaged direct and cross-coupling stiffness

图8 为时间平均下的密封直接阻尼和交叉阻尼。随着频率的增加，直接阻尼czz微增加后减小，直接阻尼cyy微减小后增加，交叉阻尼整体随频率的增加而增大。结合转子涡动轨迹可以看出，转子运动在横向上始终处于偏移状态，而在竖直方向先向下往复运动，后向上偏移。因此，y方向的运动处于反复的不稳定运动，但存在弹性约束力使得直接阻尼cyy和交叉阻尼cyz变化平缓。同时，横向的大幅度偏移使得流体对转子始终存在一个阻滞的作用，尤其当涡动频率较大时，转子涡动幅度增加，密封气动作用在横向显著。因此有效阻尼czz在50～200 Hz 波动较小，而在200～350 Hz 范围内变化显著。

图8 时间微元平均的直接阻尼和交叉阻尼Fig.8 Time-averaged direct and cross-coupling damping

图9 为时间平均的密封有效阻尼。从图中可以看出随着频率增加，有效阻尼先增大后减小，其变化趋势和转子涡动中心轨迹的变化趋势相反。这表明转子涡动所形成的密封气动效应与其稳定涡动的轨迹是相关的，而在升速过程中转子受密封气动作用会发生偏移，这也会进一步影响密封气动效应所形成的动力特性。

图9 时间微元平均的有效阻尼Fig.9 Time-averaged effective damping

综上分析可知，在升转速过程中密封动力系数会发生显著的变化，与同中心的多频涡动模型结果相比，密封的刚度和阻尼则失去了原有的对偶性。一方面，在非线性涡动模型中，转子涡动的稳定点由转子刚度、不平衡离心力和密封流体力决定，因此其涡动中心与系统坐标中心并不重合。另一方面，非线涡动模型中包含了转子刚度、不平衡离心力和密封流体力等影响，升速过程中转子中心偏移会包含更多影响频率，即密封动力特性系数还存在偏离工作频率的其他频率扰动。这也是密封时均有效阻尼发生显著波动的主要原因。

2.2 时频关联性分析

从上述分析可以看出，在连续转速变化的过程中，密封-转子系统包含了更多的扰动频率，且该过程中密封动力特性与时间和频率均有较强的关联性。为清晰辨识其他频率的影响，本文采用cmor2-2 小波信号处理方法对密封动力系数进行瞬时识别分析，其中小波识别尺度和采样周期均为2。经小波筛选后的密封力与位移数据可用来求解密封动力系数，具体过程参见文献［30］中密封动力系数求解部分。

获得的密封流体力的小波系数绝对值如图10 所示。结合图4 中时间与转速分配情况，从图10 可以看出，随着不同时间转速变化，不同频率下的密封流体力系数是显著增强的。而各个阶段和频率的流体力系数呈现连续变化表明了密封气流作用在相邻频率之间是互相干扰的。随着自转速度的增加，转子涡动会由于不平衡质量力而增强，由此也导致高频区域的密封流体力系数增大。以上分析表明本文的小波分析方法及其相关系数的确定较为合理，可以识别出本研究问题的主要影响频率。

图10 流体力小波系数绝对值Fig.10 Absolute value of fluid force wavelet coefficient

直接刚度和交叉刚度的小波系数如图11 所示。直接刚度kzz在运动初始阶段包含了大量的低频扰动，这是由于该阶段内转子在y方向的初始扰动和回复运动所形成的直流分量。通过对波动较大的系数统计发现，kzz小波系数的高幅值波动多集中在工作频率附近，但并未与工作频率相等。其中100 Hz、150 Hz 以及350 Hz 附近的波动较为明显。同理，直接刚度kyy在转子发生横向偏移时也出现的密频分布，即在0.025～0.15 s时间范围内频率波幅密集。其中在100 Hz、150 Hz 以及350 Hz 处同样具有明显的波动，而且在转速频率350 Hz 以外还存在更高的扰动频率，如379 Hz、385 Hz 处，结合时均直接刚度的变化（图8），这与kyy在高频范围内的变化较为显著是一致的。

图11 直接刚度和交叉刚度小波系数Fig.11 Wavelet coefficient of direct and cross-coupling stiffness

交叉刚度kzy的小波系数与kyy类似，在0.025～0.15 s 均存在密频，而且该密频范围多包含270～350 Hz 频率，表明在高频范围内，其变化较为明显。交叉刚度kyz所包含的扰动频率对应的系数波幅较为相近，除350 Hz 下的扰动较为明显外，其他频率下的系数波幅相近，因此其随频率的变化较为平缓，而在高频时达到较大值。

直接阻尼和交叉阻尼的小波系数如图12 所示。直接阻尼czz的小波系数同样包含个工作频率附近的扰动，其中低频和高频附近的扰动较多，集中在50～100 Hz 和250～300 Hz 范围内。其高频波动分布在升速的初始时刻，而低频分布在高转速时刻。表明转子升速运动过程中分频和倍频对直接阻尼czz的影响较为明显。直接阻尼cyy的扰动频率多集中在升速过程的中间时刻，其中50～200 Hz 区间内的扰动幅值多为负值。

图12 直接阻尼和交叉阻尼小波系数Fig.12 Wavelet coefficient of direct and cross-coupling damping

交叉阻尼的密频分布范围与交叉刚度相似，交叉阻尼czy的密频在升速中间阶段，交叉阻尼cyz的密频在初始阶段，在密频区域内同样存在显著波动的高频波幅，因此在高转速阶段的交叉阻尼变化显著。

通过上述结果可以看出，转子升速过程中的非线性运动包含工频以外的关联频率导致密封动力系数显著。值得注意的是，在低转速运行区间，直接刚度和交叉刚度存在较为显著的高频波幅。而在高转速区间，直接阻尼存在较为显著的低频波幅。这表明转子升速过程中，其非线性运动的分频和倍频影响较为明显。产生这些频率主要是由于转子非线性运动具有较高的自由性，即转子运动轨迹与系统受力相关。升速初始阶段，转子由扰动位置运动到稳定涡动位置的过程中主要受弹性回复力作用，沿竖直方向运动，因此密封流体的阻尼作用较为明显。而当转速较高时，转子轴心横向偏移和较大的涡动半径会导致流体产生较强周向旋转，产生较强的气动作用，初始转子中心抬升，因此密封流体的支撑刚度作用较显著。

对刚度和阻尼的小波系数差分可得到有效阻尼的小波系数，如图13 所示。从底部投影图可以看出，在初始阶段，有效阻尼小波系数的波动多集中在低频范围，随着转速的升高，有效阻尼小波系数的波动逐渐向高频区域迁移。而在升速的末尾阶段同时存在高频和低频的系数波动，表明密封动力特性更容易在高频区产生分频。从侧面投影图可以看出有效阻尼小波系数随频率的分布情况，在200 Hz 附近，有效阻尼的小波系数分布较为集中，且具有较高的幅值，而在350 Hz 附近有效阻尼幅值较小。这与图9 的时均有效阻尼变化相同。从时序分布来看，有效阻尼的小波系数在初始阶段的分布较为集中，且波幅较高，表明了弹性刚度的回复作用会形成较高的有效阻尼。