2018, Vol. 39
转子与静子间的密封结构是航空发动机、燃气轮机和压缩机等透平机械的重要组成部件[1]，起着降低工作介质泄漏和提高工作效率的关键作用[2, 3]。近年来，透平机械工作介质不断向高压比、高转速和高温度等高参数方向发展，由密封引起的气流激振力成为转子失稳的重要因素[4, 5]。而密封动力特性系数是评价透平机械转子稳定性的重要参数，准确预测密封动力特性系数至关重要[6, 7]。
国内外研究者最初采用控制体方法建立密封动力特性求解模型，如图 1所示，控制体模型包括单控体模型、双控体模型、三控体模型和流线型模型。Iwatsubo [8]提出将一个密封腔作为一个控制体单元的单控体密封动力特性求解模型。Childs等[9, 10]通过研究发现单控体模型未考虑气流轴向速度的影响且认为气流周向速度沿径向不变，这与实际情况不符，并提出将密封间隙作为喷射区和密封腔室作为核心区的两个控制体的双控体模型。Dan等[11]、张万福等[12]和赵三星等[13]分别应用双控体模型分析了转子偏心对密封动力特性的影响。Arqhir等[14]在双控体模型的前提下，将密封齿顶间隙作为第三个控制体，提出三控体模型，该求解模型更细节地描述了密封流场的变化情况。Komeoka [15]认为核心区与喷射区应形成倾斜的交界线并提出了流线型模型，但交界线的倾斜角很难确定。研究表明，控制体方法的优点是计算速度快，但缺少对流动细节特征的描述，且只适用于简单的密封结构，当运行工况复杂时，求解误差较大[16]。
|
Fig. 1 Numerical methods of seal rotordynamic coefficients |
随着计算机功能的提高和CFD技术的不断发展，研究者们开始应用CFD方法建立密封动力特性求解模型。如图 1所示，CFD方法分为两类：稳态旋转坐标系方法和瞬态动网格方法。Hirano等[17]应用稳态旋转坐标系方法求解了迷宫密封动力特性，并将所得结果与双控体模型的计算结果进行了对比。孙婷梅等[18]、刘晓峰[19]和陈陆淼等[20]分别应用稳态旋转坐标系方法研究了高低齿、压比、转速、偏心率和间隙变化等因素对密封动力特性的影响。瞬态动网格方法包括平动法和单/多频法。瞬态平动法[21]是通过给转子施加微小位移扰动求得刚度系数，施加微小速度扰动求得阻尼系数。Thorat等[22]和Giuseppe等[23]通过实验发现密封动力特性系数具有转子涡动频率相关性。Xin等[24]、Chocgua等[25]、Nielsen等[26]和Sreedharan等[27]应用瞬态单频法分析了转子正弦直线、圆、椭圆三种涡动轨迹下的密封动力特性系数随转子涡动频率变化的规律。为了获得多个涡动频率下的密封动力特性系数，瞬态单频法需进行多次非定常计算，耗时长，Jun等和Zhigang等[28~31]提出了瞬态多频法，只需两次非定常求解，便可得到多个转子涡动频率下的密封动力特性系数。研究表明，CFD方法能够获得密封内流场细节，计算结果更加接近密封真实工作情况，且能计算复杂的密封结构和运行工况[32]。
现有文献关于密封动力特性研究方法比较单一，综合分析密封动力特性求解模型的文献较少。本文分别建立了基于控制体方法的双控体求解模型和基于CFD方法的稳态旋转坐标系求解模型与转子平动、转子单/多频椭圆涡动求解模型。综合比较分析了双控体法，稳态旋转坐标系法和瞬态平动法、单/多频法的求解结果，研究了不同方法的优缺点及适用范围，并研究了转子涡动频率对密封动力特性的影响。
2 密封动力特性理论求解模型 2.1 控制体模型图 2给出了控制体模型的示意图。单控体模型是把一个密封腔室作为一个控制体单元，如图 2(a)所示；双控体模型将密封流场分为喷射区和核心区两部分，如图 2(b)所示；三控体模型将双控体模型的齿顶间隙作为第三个控制体，如图 2(c)所示；流线型模型的核心区与喷射区交界线为倾斜的，如图 2(d)所示。
|
Fig. 2 Control volume model |
双控体模型是目前控制体方法常用的密封动力特性求解模型，如图 3所示。双控体模型的密封间隙和密封腔室分别作为控制体Ⅰ(CVⅠ)和控制体Ⅱ (CVⅡ)，其中，mI为流入控制体Ⅰ的轴向流量，H为密封间隙，mr为两控制之间的流量，mI+1为流出控制体Ⅰ的轴向流量，R1为密封顶部半径，R2为转子半径，L为密封腔室宽度，B为齿高。
|
Fig. 3 Geometry and coordinate system of the seal |
双控体模型第i腔室的控制方程表示为
| $ \left\{ \begin{array}{l} \frac{{\partial \rho {A_1}}}{{\partial t}} + \frac{{\partial \rho {w_1}{A_1}}}{{{R_1}\partial \theta }} + {m_{{\rm{I}} + 1}} - {m_{\rm{I}}} + \frac{{\partial \rho {A_2}}}{{\partial t}} + \frac{{\partial \rho {w_2}{A_2}}}{{{R_2}\partial \theta }} = 0\\ \rho {A_1}\frac{{\partial {w_1}}}{{\partial t}} + \frac{{\rho {w_1}{A_1}}}{{{R_1}}}\frac{{\partial {w_1}}}{{\partial \theta }} + {m_{\rm{I}}}\left( {{w_{1i}} - {w_{1i - 1}}} \right) + \\ \left( {\frac{{\partial \rho {A_2}}}{{\partial t}} + \frac{{\partial \rho {w_2}{A_2}}}{{{R_2}\partial \theta }}} \right)\left( {{w_0} - {w_1}} \right) = - \frac{{{A_1}}}{{{R_1}}}\frac{{\partial p}}{{\partial \theta }} + \tau L - {\tau _{\rm{s}}}{A_{\rm{s}}}L\\ \rho {A_2}\frac{{\partial {w_2}}}{{\partial t}} + \frac{{\rho {w_2}{A_2}}}{{{R_2}}}\frac{{\partial {w_2}}}{{\partial \theta }} + \left( {\frac{{\partial \rho {A_2}}}{{\partial t}} + \frac{{\partial \rho {w_2}{A_2}}}{{{R_2}\partial \theta }}} \right)\cdot\\ \left( {{w_2} - {w_0}} \right) = - \frac{{{A_2}}}{{{R_2}}}\frac{{\partial p}}{{\partial \theta }} - {\tau _{\rm{r}}}{A_{\rm{r}}}L + \tau L\\ p = \rho RT \end{array} \right. $ | (1) |
式中A1，A2分别为两控制体的周向通流面积，Ar为两控制体之间的周向通流面积，w1，w2为控制体Ⅰ，Ⅱ的平均周向速度，w0为两控制体之间的平均周向速度，τr为转子表面切应力，τs为静子表面切应力，τ为自由切应力，p为气体压力，ρ为流体密度，R为气体常数，T为气体温度。
当对位于平衡位置的转子施加微小扰动时，相应的压力和流量等参数可以泰勒展开为
| $ \varphi = {\varphi _0} + \varepsilon {\varphi _1} $ | (2) |
式中φ为压力和温度等流动参数，φ0为零阶摄动参数，φ1为一阶摄动参数，ε为转子偏心率。
将式(2)代入控制方程(1)中，将得到的方程分解为零阶和一阶摄动方程，通过求解零阶摄动方程可得到密封的定常流场，通过求解一阶方程可得到压力分布，对压力进行积分可得到密封动力特性系数，如式(3)所示
| $ \left\{ \begin{array}{l} {K_{xx}} = \sum\limits_{z{\rm{ = 1}}}^{{N_{\rm{T}}}-1} {} {\rm{ \mathsf{ π} }}{R_2}L\Re \left( {\tilde P_{1, z}^{cx}} \right), {K_{xy}} = - \sum\limits_{z{\rm{ = 1}}}^{{N_{\rm{T}}}-1} {} {\rm{ \mathsf{ π} }}{R_2}L\Im \left( {\tilde P_{1, z}^{cy}} \right)\\ {K_{yx}} = \sum\limits_{z{\rm{ = 1}}}^{{N_{\rm{T}}}-1} {} {\rm{ \mathsf{ π} }}{R_2}L\Re \left( {\tilde P_{1, z}^{sx}} \right), {K_{xy}} = - \sum\limits_{z{\rm{ = 1}}}^{{N_{\rm{T}}}-1} {} {\rm{ \mathsf{ π} }}{R_2}L\Im \left( {\tilde P_{1, z}^{sy}} \right)\\ {C_{xx}} = \sum\limits_{z{\rm{ = 1}}}^{{N_{\rm{T}}}-1} {} {\rm{ \mathsf{ π} }}{R_2}L\Im \left( {\tilde P_{1, z}^{cx}} \right)/\mathit{\Omega} , {C_{xy}} = \sum\limits_{z{\rm{ = 1}}}^{{N_{\rm{T}}}-1} {} {\rm{ \mathsf{ π} }}{R_2}L\Re \left( {\tilde P_{1, z}^{cy}} \right)/\mathit{\Omega} \\ {C_{yx}} = \sum\limits_{z{\rm{ = 1}}}^{{N_{\rm{T}}}-1} {} {\rm{ \mathsf{ π} }}{R_2}L\Im \left( {\tilde P_{1, z}^{sx}} \right)/\mathit{\Omega} , {C_{xy}} = \sum\limits_{z{\rm{ = 1}}}^{{N_{\rm{T}}}-1} {} {\rm{ \mathsf{ π} }}{R_2}\Re L\left( {\tilde P_{1, z}^{sy}} \right)/\mathit{\Omega} \end{array} \right. $ | (3) |
式中NT为齿腔数，Kxx，Kxy，Kyx和Kyy为刚度系数，Cxx，Cxy，Cyx和Cyy为阻尼系数，s和c分别为正弦项和余弦项，
当转子受到微小位移和速度扰动时，可将转子受到的气流力与扰动位移和扰动速度的关系线性化表示为[33]
| $ - \left[ {\begin{array}{*{20}{l}} {{\rm{\Delta }}{F_x}}\\ {{\rm{\Delta }}{F_y}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{K_{xx}}{\rm{}}{K_{xy}}}\\ {{K_{yx}}{\rm{}}{K_{yy}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{\rm{\Delta }}x}\\ {{\rm{\Delta }}y} \end{array}} \right] + \left[ {\begin{array}{*{20}{l}} {{C_{xx}}{\rm{}}{C_{xy}}}\\ {{C_{yx}}{\rm{}}{C_{yy}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{\rm{\Delta }}\dot x}\\ {{\rm{\Delta }}\dot y} \end{array}} \right] $ | (4) |
在稳态旋转坐标系求解模型中，转子既绕其轴心自转，又绕静子中心做同心圆轨迹涡动。转子在做同心圆轨迹涡动时，密封齿腔内静态流动参数沿周向均匀分布，所以密封动力特性系数直接项相等，交叉项大小相等，符号相反，即Kxx = Kyy = K，Kxy = -Kyx = k，Cxx = Cyy = C，Cxy = -Cyx = c。
如图 4所示，将坐标系移到转子几何中心并使其以转子涡动转速Ω旋转，在旋转坐标系下，转子绕其几何中心旋转的角速度为ω-Ω，静子绕原点旋转的角速度为-Ω。
为了简化计算，取转子中心旋转到原点正下方的位置，此时Δx = 0，Δy = -δ，Δ
| $ \left\{ {\begin{array}{*{20}{l}} {\frac{{{F_\tau }}}{\delta } = k - C \mathit{\Omega} }\\ {\frac{{{F_r}}}{\delta } = - K - c\mathit{\Omega} } \end{array}} \right. $ | (5) |
|
Fig. 4 Frame transfer from stationary to rotating |
由式(5)可以看出，当转子做同心圆轨迹涡动时，密封气流切向力和径向力与转子涡动角速度Ω成线性关系。通过拟合转子涡动角速度Ω与其对应的密封气流力得到线性方程，该方程的斜率即为阻尼系数，截距即为刚度系数。
2.3 瞬态CFD动网格模型 2.3.1 转子平动模型图 5给出了平动法的密封动力特性求解模型，如图 5所示，在转子x方向上施加一个微小的位移扰动Δx，分别计算转子移动前后转子x和y方向向上受到的气流力，代入到式(6)可得到Kxx和Kyx。同理，在转子的y方向向上施加位移扰动可获得Kxy和Kyy。
| $ \left\{ {\begin{array}{*{20}{l}} {{K_{xx}} = \frac{{{\rm{\Delta }}{F_x}}}{{{\rm{\Delta }}x}}}\\ {{K_{yx}} = \frac{{{\rm{\Delta }}{F_y}}}{{{\rm{\Delta }}x}}} \end{array}} \right. $ | (6) |
| $ \left\{ {\begin{array}{*{20}{l}} {{C_{xx}} = \frac{{{\rm{\Delta }}{F_x}}}{{{\rm{\Delta }}\dot x}}}\\ {{C_{yx}} = \frac{{{\rm{\Delta }}{F_y}}}{{{\rm{\Delta }}\dot x}}} \end{array}} \right. $ | (7) |
|
Fig. 5 Translation model |
在转子的水平方向上施加一个微小的扰动速度Δ
图 6给出了考虑转子涡动频率时，转子做直线、圆、椭圆三种涡动的轨迹。
|
Fig. 6 Sketch diagram of three kinds of precession orbits |
(1) 单频法
转子的涡动方程为
| $ \left\{ {\begin{array}{*{20}{l}} {x = b{\rm{sin}}\mathit{\Omega} t}\\ {y = a{\rm{cos}}\mathit{\Omega} t} \end{array}} \right. $ | (8) |
| $ \left\{ {\begin{array}{*{20}{l}} {x = a{\rm{cos}}\mathit{\Omega} t}\\ {y = b{\rm{sin}}\mathit{\Omega} t} \end{array}} \right. $ | (9) |
当a=b≠0时为圆轨迹；当a≠b≠0时为椭圆轨迹；a≠b=0或b≠a=0时为正弦直线轨迹。本文以椭圆涡动轨迹为例，为了获得密封的8个动力特性系数需要一个有8个方程的方程组，任选转子的两个运动状态代入式(8)，为了简化计算，取Ωt=2πn和Ωt=π/2+2πn(n为周期个数)两个时间下转子运动状态代入式(8)，将得到的扰动位移和速度代入到式(4)，化简可得到式(10)和式(11)。同理，将这两个时间下转子运动状态代入式(9)，化简后为式(12)和(13)。通过联立求解式(10)~(13)便可得到密封的8个动力特性系数。
Ωt = 2πn时
| $ \left\{ {\begin{array}{*{20}{l}} { - {\rm{\Delta }}{F_{x1}} = {K_{xy}}a + {C_{xx}}b\mathit{\Omega} }\\ { - {\rm{\Delta }}{F_{y1}} = {K_{yy}}a + {C_{yx}}b\mathit{\Omega} } \end{array}} \right. $ | (10) |
Ωt = π/2 + 2πn时
| $ \left\{ {\begin{array}{*{20}{l}} { - {\rm{\Delta }}{F_{x2}} = {K_{xx}}b - {C_{xy}}a\mathit{\Omega} }\\ { - {\rm{\Delta }}{F_{y2}} = {K_{yx}}b - {C_{yy}}a\mathit{\Omega} } \end{array}} \right. $ | (11) |
Ωt = 2πn时
| $ \left\{ {\begin{array}{*{20}{l}} { - {\rm{\Delta }}{F_{x3}} = {K_{xx}}a + {C_{xy}}b\mathit{\Omega} }\\ { - {\rm{\Delta }}{F_{y3}} = {K_{yx}}a + {C_{yy}}b\mathit{\Omega} } \end{array}} \right. $ | (12) |
Ωt = π/2 + 2πn时
| $ \left\{ {\begin{array}{*{20}{l}} { - {\rm{\Delta }}{F_{x4}} = {K_{xy}}b - {C_{xx}}a\mathit{\Omega} }\\ { - {\rm{\Delta }}{F_{y4}} = {K_{yy}}b - {C_{yx}}a\mathit{\Omega} } \end{array}} \right. $ | (13) |
(2) 多频法
转子的涡动方程为
| $ \left\{ {\begin{array}{*{20}{l}} {x = a \cdot \sum\limits_{z = 1}^N {} {\rm{cos}}\left( {{\mathit{\Omega} _z}t} \right)}\\ {y = b \cdot \sum\limits_{z = 1}^N {} {\rm{sin}}\left( {{\mathit{\Omega} _z}t} \right)} \end{array}} \right. $ | (14) |
| $ \left\{ {\begin{array}{*{20}{l}} {x = b \cdot \sum\limits_{z = 1}^N {} {\rm{cos}}\left( {{\mathit{\Omega} _z}t} \right)}\\ {y = a \cdot \sum\limits_{z = 1}^N {} {\rm{sin}}\left( {{\mathit{\Omega} _z}t} \right)} \end{array}} \right. $ | (15) |
对式(4)进行快速傅里叶变换(FFT)可得到频域内密封气流力变化量与密封动力特性系数和微小扰动位移的关系式为
| $ \left\{ {\begin{array}{*{20}{l}} { - {\rm{\Delta }}{F_x} = \left( {{K_{xx}} + {\rm{j}}\mathit{\Omega} {C_{xx}}} \right) \cdot {D_x} + \left( {{K_{xy}} + {\rm{j}}\mathit{\Omega} {C_{xy}}} \right) \cdot {D_y}}\\ { - {\rm{\Delta }}{F_y} = \left( {{K_{yy}} + {\rm{j}}\mathit{\Omega} {C_{yy}}} \right) \cdot {D_y} + \left( {{K_{yx}} + {\rm{j}}\mathit{\Omega} {C_{yx}}} \right) \cdot {D_x}} \end{array}} \right. $ | (16) |
式中
| $ \left\{ {\begin{array}{*{20}{l}} { - {\rm{\Delta }}{F_{xx}} = \left( {{K_{xx}} + {\rm{j}}\mathit{\Omega} {C_{xx}}} \right) \cdot {D_{xx}} + \left( {{K_{xy}} + {\rm{j}}\mathit{\Omega} {C_{xy}}} \right) \cdot {D_{xy}}}\\ { - {\rm{\Delta }}{F_{xy}} = \left( {{K_{yy}} + {\rm{j}}\mathit{\Omega} {C_{yy}}} \right) \cdot {D_{xy}} + \left( {{K_{yx}} + {\rm{j}}\mathit{\Omega} {C_{yx}}} \right) \cdot {D_{xx}}} \end{array}} \right. $ | (17) |
| $ \left\{ {\begin{array}{*{20}{l}} { - {\rm{\Delta }}{F_{yy}} = \left( {{K_{yy}} + {\rm{j}}\mathit{\Omega} {C_{yy}}} \right) \cdot {D_{yy}} + \left( {{K_{yx}} + {\rm{j}}\mathit{\Omega} {C_{yx}}} \right) \cdot {D_{yx}}}\\ { - {\rm{\Delta }}{F_{yx}} = \left( {{K_{xx}} + {\rm{j}}\mathit{\Omega} {C_{xx}}} \right) \cdot {D_{yx}} + \left( {{K_{xy}} + {\rm{j}}\mathit{\Omega} {C_{xy}}} \right) \cdot {D_{yy}}} \end{array}} \right. $ | (18) |
式中Duv和ΔFuv的第一个下角标u为扰动位移激励的方向，第二个下角标v为在u方向的扰动位移激励下，v方向的转子涡动位移和密封气流力的变化量。定义密封阻抗系数为
| $ \left\{ {\begin{array}{*{20}{l}} {{H_{xx}} = {K_{xx}} + {\rm{j}}\mathit{\Omega} {C_{xx}}}\\ {{H_{xy}} = {K_{xy}} + {\rm{j}}\mathit{\Omega} {C_{xy}}}\\ {{H_{yx}} = {K_{yx}} + {\rm{j}}\mathit{\Omega} {C_{yx}}}\\ {{H_{yy}} = {K_{yy}} + {\rm{j}}\mathit{\Omega} {C_{yy}}} \end{array}} \right. $ | (19) |
将式(19)代入到式(17)和(18)后阻抗系数变为
| $ \left\{ \begin{array}{l} {H_{xx}} = \frac{{\left( { - {F_{xx}}} \right) \cdot {D_{yy}} - \left( { - {F_{yx}}} \right){D_{xy}}}}{{{D_{xx}} \cdot {D_{yy}} - {D_{yx}} \cdot {D_{xy}}}}\\ {H_{xy}} = \frac{{\left( { - {F_{xx}}} \right) \cdot {D_{yx}} - \left( { - {F_{yx}}} \right) \cdot {D_{xx}}}}{{{D_{xy}} \cdot {D_{yx}} - {D_{yy}} \cdot {D_{xx}}}}\\ {H_{yy}} = \frac{{\left( { - {F_{yy}}} \right) \cdot {D_{xx}} - \left( { - {F_{xy}}} \right) \cdot {D_{yx}}}}{{{D_{xx}} \cdot {D_{yy}} - {D_{yx}} \cdot {D_{xy}}}}\\ {H_{yx}} = \frac{{\left( { - {F_{yy}}} \right) \cdot {D_{xy}} - \left( { - {F_{xy}}} \right) \cdot {D_{yy}}}}{{{D_{xy}} \cdot {D_{yx}} - {D_{yy}} \cdot {D_{xx}}}} \end{array} \right. $ | (20) |
密封动力特性系数为
| $ \left\{ {\begin{array}{*{20}{l}} {{K_{xx}} = {\rm{Re}}\left( {{H_{xx}}} \right)}\\ {{C_{xx}} = \frac{{{\rm{Im}}\left( {{H_{xx}}} \right)}}{\mathit{\Omega }}}\\ {{K_{xy}} = {\rm{Re}}\left( {{H_{xy}}} \right)}\\ {{C_{xy}} = \frac{{{\rm{Im}}\left( {{H_{xy}}} \right)}}{\mathit{\Omega }}}\\ {{K_{yy}} = {\rm{Re}}\left( {{H_{yy}}} \right)}\\ {{C_{yy}} = \frac{{{\rm{Im}}\left( {{H_{yy}}} \right)}}{\mathit{\Omega }}}\\ {{K_{yx}} = {\rm{Re}}\left( {{H_{yx}}} \right)}\\ {{C_{yx}} = \frac{{{\rm{Im}}\left( {{H_{yx}}} \right)}}{\mathit{\Omega }}} \end{array}} \right. $ | (21) |
本文求解模型为文献[17]的迷宫密封模型，该模型的转子表面光滑且静子面带有5个锥形直齿，几何参数如图 7所示，工况参数如表 1所示。
|
Fig. 7 Size of the seal(mm) |
|
Table 1 Seal model work condition parameters |
本文采用双控制体模型求解密封动力特性系数，根据密封的结构参数、边界条件和经验参数建立控制方程组(1)。将方程组分解为零阶方程和一阶摄动方程，求解一阶摄动方程得到密封内的压力分布，由式(3)对压力进行积分可得到密封的8个动力特性系数。
3.3 CFD方法 3.3.1 网格划分密封齿腔部位的轴向剖面网格划分如图 8所示。为划分高质量网格来增加求解准确性，密封整体流场的网格结构均选用六面体。为了细化流动特性变化较大的近壁面区域，相邻两节点的间距比定为1.1，增加近壁面处的网格密度。为了选择合适的网格密度，本文首先分析了网格密度对泄漏量的影响，当泄漏量随网格密度变化率小于5%时，再增加网格密度，泄漏量基本保持不变，最终将密封齿顶间隙的轴向和径向节点数分别设置为10和15，密封齿腔的轴向和径向节点数分别为50和25，密封周向节点数设置为96。
|
Fig. 8 Grid of seal axil section |
本文分别求解了转子涡动角速度在0，290.5，581，871.5，1162rad/s下转子所受到的径向和切向密封气流力。当计算残差达到控制标准10-6，密封进出口流量差小于0.1%，密封的径向和切向气流力相邻两次迭代的变化量均小于5%时认为计算已经收敛。
图 9给出了转子受到的径向和切向密封气流力与转子涡动角速度的关系，由图 9可以看出转子在做同心圆轨迹涡动时，转子受到的径向和切向密封气流力与转子涡动角速度近似成线性化关系，将密封的径向气流力和切向气流力分别与转子涡动角速度进行线性拟合，由式(5)可得出密封的动力特性系数。
|
Fig. 9 Seal force versus whirl speed |
(1) 瞬态平动法
瞬态平动法在求解刚度系数时只需建立转子偏移前后的CFD模型，然后分别求解作用在转子x和y方向的力。在求解密封的阻尼系数时，使用转子偏移后的CFD模型，应用CFD软件的瞬态动网格技术来实现转子的匀速直线运动的条件。本文应用CFX软件的自定义函数(UDF)给转子施加一个匀速直线的动边界条件。同时考虑非定常计算下各个时间步的收敛性和求解时间，最终将时间步长设置为0.1ms，并且根据转子偏移距离和设置的微小扰动速度决定时间步数，且需保证每一时间步下计算残差达到控制标准10-6。
(2) 瞬态单/多频法
转子做单频椭圆涡动时将转子涡动位移振幅a，b分别设置为：a=0.05Cr，b=0.1Cr，Cr为密封间隙。转子做多频椭圆涡动时将转子涡动位移振幅a，b分别设置为：a=0.01Cr，b=0.005Cr。本文应用CFD瞬态动网格技术实现转子的单频和多频椭圆涡动，求解了转子涡动频率在40，80，120，160，200，240，280，320Hz共8个涡动频率下的密封动力特性系数。首先对求解模型进行定常计算，此时只需设置转子的自转转速，当残差系数达到控制标准10-6时认为该定常计算已收敛。将定常计算得到的结果作为非定常计算的初值，当对求解模型进行非定常计算时，设置转子的涡动位移方程。当非定常计算下密封气流力Fx和Fy呈周期性变化，且相邻两周期转子涡动到同一位置时密封气流力的变化幅度小于0.2%时，可认为该非定常计算已收敛。
4 求解结果分析 4.1 密封流场特性分析图 10给出了本文数值模型求解得到的密封流场轴向剖面速度矢量图，从图中可以清晰看出迷宫密封流场分为齿腔内的旋流区和齿顶处的射流区两部分。
|
Fig. 10 Scattergram of speed on axial direction of seal |
图 11给出了迷宫密封流场压力沿轴向的变化情况，可以看出流场压力呈阶梯式降低且压力降主要发生在齿顶处，齿腔内的压力值基本相等。
|
Fig. 11 Scattergram of pressure on axial direction of seal |
通过对迷宫密封流场流线和压力分布的分析可得到其封严机理：迷宫密封上依次排列着许多环形密封齿，当气流经过每一个密封齿时，由于密封齿顶处流动截面小，气流的流动速度会迅速增大，将气流的压力能转变为动能，当气流进入密封齿腔后，由于齿腔内流动截面比齿顶流动截面大，气流的流动速度会突然降低，同时气流会迅速膨胀而产生一个强烈的漩涡，将气流的大部分动能转化为热量而散失掉。正是由于密封齿顶处的节流作用和密封齿腔内的动能耗散作用使气流的压力沿密封轴向逐渐降低，从而达到封严的效果。为了提高密封的封严效果，可通过增强密封齿顶处的节流作用和齿腔内的动能耗散作用。增大密封齿的尖锐度，减小密封齿的厚度可增强节流作用；适当增加密封腔室的高度，使得涡流刚好充满腔室，此时动能耗散作用最强；减小密封间隙可减小直通泄漏面积，进而可减小泄漏量。
4.3 转子涡动位移对动力特性求解精度的影响本文以稳态旋转坐标系法为例，分析了转子涡动半径大小对密封动力特性系数求解精度的影响。图 12给出了不同转子涡动半径下，旋转坐标系法与文献[17]的平均相对误差。由图 12可以看出，随着转子涡动半径的增加，平均相对误差随之增大。这由于随着涡动半径的增加，气流径向力和切向力与转子涡动频率的线性化变差，从而导致求解误差增大，故稳态旋转坐标系法只适用于小涡动半径密封动力特性求解模型。
|
Fig. 12 Relative error versus whirling radius |
图 13给出了本文CFX模型与文献[17]的双控体模型和Tascflow模型沿密封轴向压力分布的对比结果，表 3给出了泄漏量的对比结果。
|
Fig. 13 Scattergram of pressure on axial direction of seal |
|
|
Table 3 Leakage comparison |
由图 13和表 3的对比结果得出，本文CFX求解模型的流场特性计算结果与Tascflow模型计算结果[17]相近，与双控体模型相差较大。这是由于双控体模型将迷宫密封流场分为喷射区和核心区两部分，并将每部分看成一个整体，未考虑齿顶间隙处流场变化，缺少对密封流场细节的描述。
4.4.2 动力特性计算结果对比分析表 4给出了本文双控体法、稳态旋转坐标系法和瞬态平动法求得的密封动力特性系数。表 5给出了本文数值方法求得结果与文献[17]的相对误差。由表 5可以看出，与文献[17]相比，本文稳态旋转坐标系法求得各个动力特性系数的相对误差较小，这是由于两者采用的数值方法相同；本文瞬态平动法求得各个动力特性系数的相对误差也较小，但各个动力特性系数的相对误差与稳态旋转坐标系法不同，这是由于瞬态平动法是应用动网格技术求解密封动力特性系数，两者方法不同；双控体法求得各个动力特性系数的相对误差较大，这主要是由于双控体法采用的控制体方程是N-S方程的简化，采用半经验公式，动能系数和流动系数等经验参数的选取对求解结果影响很大。
|
|
Table 4 Rotordynamic characteristic coefficients |
|
|
Table 5 Relative error of rotordynamic characteristic coefficients |
本文所使用的计算机配置：至强E5-2650V3处理器，64G内存。表 6给出了不同密封动力特性数值方法计算所需时间。由表 6可以看出，未考虑转子涡动频率时，控制体方法计算所需时间最少，但该方法的求解精度较低；稳态旋转坐标系法计算所需时间少于瞬态平动法，通过两种数值方法求得的密封动力特性数对比分析已知，两种数值方法的相对偏差较小，当求解结构为周向对称密封的动力特性系数时应采用稳态旋转坐标系法，节约求解时间。考虑转子涡动频率后，在单次非定常求解下，瞬态多频法计算时间大于瞬态单频法，这是由于在相同的时间步长下，为了保证求解结果的精度，瞬态多频法所需时间步数较多，但瞬态多频法可通过两次非定常求解得到多个转子涡动频率下的密封动力特性系数，所需总时间少于瞬态单频法。
|
|
Table 6 Computation time |
(1) 转子单频涡动模型
在给转子施加单频椭圆激励时，图 14给出了转子涡动位移的时域图与频域图，图 15给出了气流力的时域图与频域图。由图 14和图 15的时域图可以看出，转子在x方向和y方向做正弦涡动，对应的气流力在x方向和y方向也呈正弦规律波动。由图 14和图 15的频域图可以看出气流力的波动频率与转子涡动频率相同。图 15气流力时域图中x方向气流力第三个峰值为第235步，大于图 14涡动位移时域图中x方向涡动位移第三个峰值对应的第200步，由此可以得出，气流力的变化相位滞后于转子涡动位移相位。图 15气流力频域图中涡动频率为0Hz时的气流力峰值为转子在偏心处自转所形成的。
|
Fig. 14 Diagram of rotor whirling displacement(single-frequency) |
|
Fig. 15 Time domain and frequency domain of air force(single-frequency) |
(2) 转子多频涡动模型
在给转子施加多频椭圆激励时，图 16给出了转子涡动位移时域图和经过快速傅里叶变换后得到的频域图。由图 16转子涡动位移频域图可以看出，本文瞬态多频法可以准确捕捉到设定在转子上的各个涡动频率成分及相应的涡动位移幅值。
|
Fig. 16 Diagram of rotor whirling displacement(multi-frequency) |
图 17给出了密封气流力时域图和经过快速傅里叶变换后得到的频域图。由图 17气流力频域图可以明显看出，8个密封气流力幅值对应的涡动频率与施加在转子上的涡动位移频率成分相同，且随着转子涡动频率的增加，在同一涡动幅值下，x方向和y方向的气流力均随涡动频率的增加而增大。
|
Fig. 17 Time domain and frequency domain of air force(multi-frequency) |
图 18给出了转子涡动频率对密封动力特性系数的影响。由图 18还可以看出，Kxx和Kyy随着转子涡动频率增加而增大，Kxy，Kyx和Cyx的绝对值随着转子涡动频率增加而减小，Cxx，Cxy和Cyy随着转子涡动频率增加而减小。与未考虑转子涡动频率的文献[17]相比，密封动力特性系数相对偏差达到6.01%~91.54%，故转子涡动频率对密封动力特性系数影响很大。在密封实际工况下，转子不仅做自转，同时还带有涡动，所以瞬态单/多频法更符合密封的实际工作情况。
|
Fig. 18 Effects of whirling frequencies of rotor on rotordynamic characteristics |
通过对现有透平机械动密封动力特性的数值方法，即控制体法、CFD稳态旋转坐标系法和CFD瞬态平动法、单/多频法的研究，得出以下结论：
(1) 迷宫密封齿顶处的节流效应和齿腔内的动能耗散效应使迷宫密封起到封严的作用。可通过增加密封齿的尖锐度和适当增加腔室的高度，减小密封齿的厚度和密封间隙，达到降低泄漏量的效果。
(2) 控制体法求解速度快，但求解精度较低；稳态旋转坐标系法只适用于转子做小轨迹同心涡动的求解模型；瞬态平动法求解速度较慢；瞬态单/多频法考虑了转子的涡动频率，更符合密封的实际工作情况。
(3) 密封气流力随着转子涡动频率的变化而变化，且气流力的频率成分与转子涡动位移的频率成分相同，但相位滞后于转子涡动位移的相位。
(4) 密封动力特性系数随着转子涡动频率的变化而变化，主刚度系数随转子的涡动频率增大而增大，交叉刚度系数和阻尼系数的绝对值随转子的涡动频率增大而减小。
| [1] |
曹树谦, 陈予恕. 现代密封转子动力学研究综述[J]. 工程力学, 2009, 26(增刊Ⅱ): 68-79. (  0) |
| [2] |
李伟, 乔渭阳, 许开富, 等. 涡轮叶尖迷宫式密封对泄漏流场影响的数值模拟[J]. 推进技术, 2009, 30(1): 88-94. (LI Wei, QIAO Wei-yang, XU Kai-fu, et al. Numerical Simulation of Labyrinth Seal on Tip Leakage Flow in Partially and Full Shrouded Axial Turbine[J]. Journal of Propulsion Technology, 2009, 30(1): 88-94.)
( 0) |
| [3] |
李钰洁, 刘永葆, 贺星. 间隙变化对新型涡轮密封气动性能影响的数值分析[J]. 推进技术, 2015, 36(8): 1179-1185. (LI Yu-jie, LIU Yong-bao, HE Xing. Numerical Simulation for Effects of Clearance Change of a Novel Seal on Aerodynamic Performance of Gas Turbine[J]. Journal of Propulsion Technology, 2015, 36(8): 1179-1185.)
( 0) |
| [4] |
钟一谔, 何衍宗, 王正, 等. 转子动力学[M]. 北京: 清华大学出版社, 1987.
( 0) |
| [5] |
贾兴运, 徐国印, 张海, 等. 转子振动对T型交错式迷宫密封性能影响[J]. 推进技术, 2017, 38(6): 1370-1378. (JIA Xing-yun, XU Guo-yin, ZHANG Hai, et al. Effects of Rotor Vibration on T Type Labyrinth Seal Performance[J]. Journal of Propulsion Technology, 2017, 38(6): 1370-1378.)
( 0) |
| [6] |
田爱梅, 朱梓根. 涡轮泵转子稳定性计算[J]. 推进技术, 2000, 21(3): 43-45. (TIAN Ai-mei, ZHU Zi-gen. Calculation for Rotor Stability in Turbopumps[J]. Journal of Propulsion Technology, 2000, 21(3): 43-45.)
( 0) |
| [7] |
高庆, 李军. 涡轮蜂窝面径向轮缘密封封严性能的数值研究[J]. 推进技术, 2016, 37(5): 937-944. (GAO Qing, LI Jun. Numerical Investigations on Sealing Performance of Turbine Honeycomb Radial Rim Seal[J]. Journal of Propulsion Technology, 2016, 37(5): 937-944.)
( 0) |
| [8] |
Iwalsubo T. Evaluation of Instability Force of Labyrinth Seals in Turbine or Compressors[R]. NASA CP2133, 1980.
( 0) |
| [9] |
Childs D W, Scharrer J K. Iwatsubo-Based Solution for Labyrinth Seals: Comparison to Experimental Results[J]. Journal of Engineering for Gas Turbines and Power, 1986, 108(4): 325-331.
( 0) |
| [10] |
Childs D W, Scharrer J K. Theory Versus Experiment for the Rotordynamic Coefficients of Labyrinth Gas Seals: Part Ⅱ A Comparison to Experiment[J]. Journal of Vibration, Acoustics, Stress, and Reliability in Design, 1988, 110(7): 281-287.
( 0) |
| [11] |
Dan Sun, Jiangang Yang, Rui Guo, et al. A Trigonometric series Expansion based Method for the Research of Static and Dynamic Characteristics of Eccentric Seals[J]. Journal of Mechanical Science and Technology, 2014, 28(6): 2111-2120. DOI:10.1007/s12206-014-0408-8
( 0) |
| [12] |
张万福, 杨建刚, 曹浩, 等. 偏心密封内切向气流力及其对稳定性影响的理论与试验研究[J]. 机械工程学报, 2011, 47(17): 92-98. ( 0) |
| [13] |
赵三星, 孟光, 冷小磊, 等. 用双控制体模型计算直通式篦齿密封的动特性系数[J]. 机械工程学报, 2006, 42(5): 103-109. ( 0) |
| [14] |
Arqhir Mihai, Frene Jean. A Bulk-Flow Analysis of Static and Dynamic Characteristics of Eccentric Circumferentially: Grooved Liquid Annular Seals[J]. Journal of Tribology, 2004, 126(4): 317-324.
( 0) |
| [15] |
Kameoka T. Theoretical Approach for Labyrinth Seal Forcescross Coupled Stiffness of a Straight-Through Labyrinth Seal[R]. NASA CP2338, 1984.
( 0) |
| [16] |
晏鑫, 蒋玉娥, 李军. 迷宫密封转子动力学特性的数值模拟[J]. 热能动力工程, 2009, 24(5): 566-570. ( 0) |
| [17] |
Hirano T, Zenglin Guo, Kirk R G. Application of Computational Fluid Dynamics Analysis for Rotating Machinery-Part Ⅱ: Labyrinth Seal Analysis[J]. Journal of Engineering for Gas Turbines and Power, 2005, 127(4): 820-826. DOI:10.1115/1.1808426
( 0) |
| [18] |
孙婷梅, 郑水英. 基于FLUENT迷宫密封动力特性分析[J]. 流体机械, 2008, 36(8): 24-27. ( 0) |
| [19] |
刘晓峰. 汽轮机气流激振CFD数值分析与转子动力稳定性研究[D]. 南京: 东南大学, 2007. http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y1234604
( 0) |
| [20] |
陈陆淼, 秦朝烨, 褚福磊. 航空发动机转子篦齿封严诱导的气流激振力分析[J]. 推进技术, 2016, 37(3): 527-533. (CHEN Lu-miao, QIN Zhao-ye, CHU Fu-lei. Analysis on Airflow Induced Forces of Aero-Engine Rotor with Labyrinth Seal[J]. Journal of Propulsion Technology, 2016, 37(3): 527-533.)
( 0) |
| [21] |
孙丹, 王双, 艾延廷, 等. 负错位密封静力与动力特性的数值研究[J]. 航空动力学报, 2017, 32(4): 791-799. ( 0) |
| [22] |
Thorat M, Childs D. Predicted Rotordynamic Behavior of a Labyrinth Seal as Rotor Surface Velocity Approaches Mach 1[R]. ASME GT2009-59256.
( 0) |
| [23] |
Giuseppe V, Thorat M, Childs D W, et al. Impact of Frequency Dependence of Gas Labyrinth Seal Rotordynamic Coefficients on Centrifugal Compressor Stability[R]. ASME GT2010-22039.
( 0) |
| [24] |
Xin Yan, Li J, Jun Li, et al. Investigations on the Rotordynamic Characteristics of a Hole-Pattern Seal Using Transient CFD and Periodic Circular Orbit Model[J]. ASME Journal of Vibration and Acoustics, 2011, 133(4).
( 0) |
| [25] |
Chochua G, Soulas T A. Numerical Modeling of Rotordynamic Coefficients for Deliberately Roughened Stator Gas Annular Seals[J]. Journal of Tribology, 2007, 129(2): 424-429. DOI:10.1115/1.2647531
( 0) |
| [26] |
Nielsen K K, Janck K, Underbakke H. Hole-Pattern and Honeycomb Seal Rotordynamic Forces: Validation of CFD-Based Prediction Techniques[J]. Journal of Engineering for Gas Turbines and Power, 2012, 134(12).
( 0) |
| [27] |
Sreedharan S S, Vannini G. CFD Assessment of Rotordynamic Coefficients in Labyrinth Seals[R]. ASME GT 2014-26999.
( 0) |
| [28] |
Jun Li, Zhigang Li, Zhenping Feng. Investigations on the Rotordynamic Coefficient of Pocket Damper Seals Using the Multi-Frequency, One-Dimensional, Whirling Orbit Model and RANS Solutions[J]. Journal of Engineering for Gas Turbines and Power, 2012, 134(10).
( 0) |
| [29] |
Zhigang Li, Jun Li, Zhenping Feng. Numerical Investigations on the Leakage and Rotordynamic Characteristics of Pocket Damper Seals, Part Ⅰ: Effects of Pressure Ratio, Rotational Speed, and Inlet Preswirl[J]. Journal of Engineering for Gas Turbines and Power, 2014, 137(3).
( 0) |
| [30] |
Zhigang Li, Jun Li, Zhenping Feng. Numerical Investigations on the Leakage and Rotordynamic Characteristics of Pocket Damper Seals, Part Ⅱ: Effects of Partition Wall Type, Partition Wall Number, and Cavity Depth[J]. ASME Journal of Engineering for Gas Turbines and Power, 2014, 137(3).
( 0) |
| [31] |
Zhigang Li, Jun Li, Xin Yan. Multiple Frequencies Elliptical Whirling Orbit Model and Transient RANS Solution Approach to Rotordynamic Coefficients of Annual Gas Seals Prediction[J]. ASME Journal of Vibration, Acoustics, Stress, and Reliability in Design, 2013, 135(3).
( 0) |
| [32] |
Zhigang Li, Jun Li, Zhenping Feng. Comparisons of Rotordynamic Characteristics Predictions for Annular Gas Seals Using the Transient CFD Method Based on Different Single-Frequency and Multi-Frequency Rotor Whirling Models[J]. ASME Journal of Tribology, 2016, 138(1).
( 0) |
| [33] |
孙丹, 王双, 艾延廷, 等. 阻旋栅对密封静力与动力特性影响的数值与实验研究[J]. 航空学报, 2015, 36(9): 3002-3011. ( 0) |