# Study on Angular Free Vibration Stability and Parameters Influence of Dry Gas Seal Based on the Characteristic Equation

Authors: Liming Teng; Jinbo Jiang; Xudong Peng; Jiyun Li; Simin Zheng
Journal: Shock and Vibration
 (2022)

---

## Abstract

Based on Laplace transform and Hurwitz stability criterion, the system characteristic equation and algebraic stability criterion of the two degree-of-freedom angular free vibration of dry gas seal are derived. The effects of rotational inertia, stiffness and damping parameters on the angular free vibration stability of dry gas seal are analyzed by using the root locus method and the concept of the closed-loop dominant pole. The results show that the constraint condition of rotational inertia is the most demanding of all stability conditions. Both the angular main damping and main stiffness are not simply the larger the better, but there are preferred values in the interval greater than the stability threshold. The approach of the absolute value of the cross coefficient to zero is beneficial to the suppression of angular free vibration. There is a threshold of rotational inertia, which makes the cross stiffness change from having only the lower limit value to both the upper and lower critical values at the same time.


---

## Body

## 1. Introduction

Dry gas seal (DGS) is a kind of noncontacting mechanical seal with dry and clean gas as the sealing medium. Because of its low leakage, low power consumption, and long service life, DGS has been successfully applied in rotating machineries such as compressor, fan, centrifugal pump, and stirred kettle, and has become the preferred type of shaft-end seal for numerous rotating machineries. Several researches have been carried out on the steady-state performance of DGS, including the optimization of groove geometric parameters [1, 2], the design of the surface texture [3–5], the thermoelastic deformation analysis of sealing rings [6, 7], and the influence of real gas effects [8, 9]. With the development of rotating machinery to high-parameter working conditions, especially under the condition of high-speed or even ultra-high-speed, the poor dynamic characteristics of DGS become increasingly prominent because of the stiffness and damping characteristics of sealing gas film, vibration or swing of the rotating shaft, and impact disturbance. The sealing pair is prone to wear caused by direct contact of the end face or excessive leakage caused by an excessive opening, which will cause fragmentation and failure of the sealing ring in severe cases. Thus, it can be seen that, in order to ensure the operational reliability of high-parameter DGS, it is necessary to conduct an in-depth study on the dynamic performance of DGS.Due to the inevitable installation deviation and the possible bending vibration of the rotating shaft under high-speed conditions, in the actual operation of DGS, the sealing ring will swing angularly along the two orthogonal axes under the action of angular excitation. The results of the current research [10–14] show that compared to the axial motion of the sealing ring caused by axial perturbation, the angular deflection of the sealing ring is more likely to cause seal failure. Therefore, it is imperative to study the angular motion law of DGS to ensure stable operation of DGS. In addition to obtaining the angular response of the disturbed sealing ring [15–17], one of the more important works is to discuss the stability characteristics of the angular motion of DGS. The condition of angular vibration stability of DGS is of great significance for predicting and evaluating the angular dynamic behavior of DGS.The stability condition of the angular motion of the sealing ring has attracted wide attention of scholars as early as the end of the last century. Etsion and Dan [18, 19] obtains the condition for angular stability of mechanical seal by using Hurwitz criterion for the first time and gives the expression of the critical speed above which the seal becomes dynamically unstable. Green and Etsion [20–22] systematically analyze the dynamic behavior of a non-contacting coned-face seal and provides the expression for the critical rotor runout and stator misalignment causing face contact. Subsequently, in the research on the steady-state response of a grooved face seal, Person et al. [23] tests the dynamic stability using the method developed by Etsion and Green for the case of narrow, smooth face seals, which proves the correctness of the previous theory about the stability threshold. Wileman and Green [24] and Kundera [25] investigate the stability of a mechanical seal with two flexibly mounted rotors. The result shows that the stability threshold strongly depends upon the ratio of the transverse to the polar moment of inertia. Green and Barnsby [26] and Miller and Green [27] find that the inertia product I·Ω2 matters for the stability threshold, and this term with units of energy should be the critical indicator. On this basis, Hu et al. [28] and Hao et al. [29] respectively derive the critical moment of inertia for angular stability by using the perturbation method and carry out the parametric research. Thus, it can be seen that current research on the stability condition of angular motion is mainly focused on smooth face seals with liquid lubrication. Due to the grooved face structure and the nonlinearity of the gas film pressure, the closed-form solution of gas film stiffness and damping properties cannot be obtained. Therefore, there is little research on the angular stability of grooved face DGS. In addition, most of the stability thresholds obtained at present are macroscopic physical quantities, such as moment of inertia, rotational speed, and installation deflection angle. However, there have been few reports on dynamic characteristics of the gas film. It is urgent to study the stability conditions of gas film dynamic characteristics of spiral groove DGS to guide seal performance prediction and structure design, and to provide a guarantee for the stable operation of high-speed rotating machinery.Only the angular excitation caused by the installation deviation of the nonfloating ring is considered. In this case, the angular motion response of the floating ring can be regarded as the response under harmonic excitation, and the angular frequency of the angular excitation is always equal to the angular velocity of the shaft. According to the theory of mechanical vibration and linear ordinary differential equations, the response of a mechanical system under harmonic excitation consists of two parts: the transient response and steady-state response. Steady-state response, i.e. the forced vibration, is a continuous constant-amplitude vibration with the frequency equal to the excitation frequency and the phase lagging behind the excitation phase. Transient response is the free vibration, which determines the characteristics of the system in transition stage and the stable state of the system vibration. Therefore, it is necessary to study the stability of two degree-of-freedom (DOF) angular free vibration of DGS and its influencing factors.Without loss of generality, this paper takes the flexibly mounted stator DGS, the most widely used under high-speed working conditions, as the research object. On the basis of the geometric and physical model used to solve the angular vibration characteristics of DGS, the system characteristic equation and algebraic stability criterion of two DOF angular free vibration of DGS are derived based on Laplace transform and Hurwitz stability criterion. Based on the dynamic response and simplified analysis methods of the high-order system, the maximum value of the real part of the characteristic root is proposed as the evaluation index of angular free vibration stability of DGS. The root locus method is used to study the single factor influence law and the interaction relationship of rotational inertia, stiffness, and damping parameters on the angular free vibration stability of DGS. The physical meaning of the cross coefficient is clarified, and the reasonable value range is pointed out. It provides a new method for angular vibration stability analysis and sealing characteristic parameter design of DGS.


## 2. Calculation Model

### 2.1. Physical and Geometric Model

Figure1(a) shows the geometric structure model of the flexibly mounted stator DGS. At the sealing face of the rotor near the outer radius ro, there are circumferential uniformly distributed logarithmic spiral grooves with a quantity of Ng and a depth of hg, and the spiral angle is α; an ungrooved sealing dam is arranged adjacent to the inner radius ri to achieve sealing during shutdown; the boundary radius between the sealing dam and the sealing groove is rg; the equilibrium film thickness at the sealing dam is hb. Due to installation deviation or shafting vibration or swing, the rotor fixed on the rotating shaft and rotating at an angular velocity Ω has a certain axial and angular motion, and applies continuous axial and angular excitation to the floating stator through the sealing gas film, in which the axial excitation amplitude and frequency are Arz and ω respectively, and the fixed installation angle deviation of the rotor is Ar.Figure 1
Schematic diagram of the geometric structure and kinematic model of the flexibly mounted stator DGS: (a) geometric structure model; (b) kinematic model and (c) positional relationship.
(a)(b)(c)The sealing gas film between rotor and stator can be regarded as a “spring-damping” structure with certain stiffness and damping. The floating stator can respond to the three DOF excitation of the rotor and keep following the disturbance of the rotor under the joint action of the thin film of gas, spring, and auxiliary sealing ring. Because the axial and angular motion can be decoupled [30, 31], only the angular disturbance is considered in this analysis.During seal transient, the transient gas film thickness is determined by the relative position of the seal rings. Assuming only rigid body motion, the position of each face can be completely described by their angular orientation. As shown in Figure1(c), the inertial coordinate system Oxyz can be considered as fixed in space. z is the axis about which the shaft rotates. Osxsyszs is the stator reference coordinate system, where axis xs is parallel to the Oxy plane. Axis ys points to the point in the stator plane of maximum distance from the Oxy plane, and zs is the central axis of the stator.


### 2.2. Mathematic Model

#### 2.2.1. Angular Free Vibration Model

Considering the stiffness of the supporting spring and the damping of the auxiliary seal, and ignoring the influence of the axial disturbance on the angular motion response, according to Newton’s law of motion, the two DOF angular free vibration motion equation of the floating ring can be obtained as follows:(1)Ixφ¨x=Mx′−Mxsp′−Mxse′,Iyφ¨y=My′−Mysp′−Myse′,where Ix and Iy are the rotational inertia of the floating ring around the x-axis and y-axis, respectively; φ¨x and φ¨y are the angular acceleration of the response of the floating ring around the x-axis and y-axis, respectively; Mx′ and My′ are the x-direction and y-direction deflection moment increment of the gas film under perturbation respectively, which can be obtained from the Taylor expansion and the definition of the dynamic characteristic coefficient of the gas film:(2)Mx′=−kxxφx−kxyφy−cxxφ˙x−cxyφ˙y,My′=−kyxφx−kyyφy−cyxφ˙x−cyyφ˙y,where φx, φy, φ˙x and φ˙y are the angular displacement and angular velocity of the angular response of the floating ring, respectively, and the subscript indicates the direction of the corresponding physical quantities; kij and cij (i, j = x, y) are the angular dynamic stiffness and damping coefficients of the gas film, respectively. The first subscript indicates the direction of the moment. The second subscript indicates the direction of the disturbing displacement or velocity that causes the moment change. The negative sign indicates that the increment of the moment always prevents the displacement and velocity of the disturbance.Mxsp′ and Mysp′ are the deflection recovery moment of the supporting spring to the floating ring caused by the displacement perturbation in the x and y directions, respectively. If ksx and ksy are used to express the equivalent recovery moment stiffness of the spring system in the corresponding direction under the angular displacement disturbance φx and φy respectively, then(3)Mxsp′=ksxφx,Mysp′=ksyφy.Mxse′ and Myse′ are the deflection recovery moment of the auxiliary seal to the floating ring caused by the velocity perturbation in the x and y directions, respectively. If csx and csy are used to express the equivalent recovery moment damping of the auxiliary seal in the corresponding direction under the angular velocity disturbance φ˙x and φ˙y respectively, then(4)Mxse′=csxφ˙x,Myse′=csyφ˙y.By substituting equations (2)–(4) into equation (1), the equation of motion of the angular free vibration of the floating ring can be transformed into(5)Ixφ¨xt+cxx+csxφ˙xt+cxyφ˙yt+kxx+ksxφxt+kxyφyt=0,Iyφ¨yt+cyxφ˙xt+cyy+csyφ˙yt+kyxφxt+kyy+ksyφyt=0.


#### 2.2.2. Film Dynamic Coefficients

In this paper, the perturbation method is used to calculate the linearized angular stiffness and damping coefficients appearing in equation (2). The compressible perturbed Reynolds equations are as follows (the theoretical derivation refers to references [31, 32]):(6a)1r∂∂r3r2h02sinθ2μ∂p02∂r+rh03μ∂p0pαr∂r+1r2∂∂θ3rh02sinθ2μ∂p02∂θ+h03μ∂p0pαr∂θ−6Ω∂pαrh0+p0rsinθ∂θ+12ωpαih0=0,1r∂∂rrh03μ∂p0pαi∂r+1r2∂∂θh03μ∂p0pαi∂θ−6Ω∂pαih0∂θ−12ωpαrh0+p0rsinθ=0,(6b)1r∂∂rrh03μ∂p0pβr∂r−3r2h02cosθ2μ∂p02∂r+1r2∂∂θh03μ∂p0pβr∂θ−3rh02cosθ2μ∂p02∂θ−6Ω∂pβrh0−p0rcosθ∂θ+12ωpβih0=0,1r∂∂rrh03μ∂p0pβi∂r+1r2∂∂θh03μ∂p0pβi∂θ−6Ω∂pβih0∂θ−12ωpβrh0−p0rcosθ=0,where p0 and h0 are the steady-state gas film pressure and thickness, respectively. pkj (k = α,β; j = r,i) is the first order perturbed pressure. The first subscript represents the direction of the perturbed pressure, and the second subscript represents the real part or imaginary part. ω is the frequency of angular excitation motion. In general cases (such as the angular excitation introduced by installation deviations of seal rings), the excitation frequency ω is equal to the rotational angular velocity Ω.Solving equation (6) to obtain gas film perturbed pressure distribution, and eight dynamic stiffness and damping coefficients can be obtained by integrating on the whole seal face according to the following formula:(7a)kxxkxykyxkyy=−∫02π∫riropαrrsinθpβrrsinθ−pαrrcosθ−pβrrcosθrdrdθ,(7b)cxxcxycyxcyy=−∫02π∫riropαirsinθωpβirsinθω−pαircosθω−pβircosθωrdrdθ.


#### 2.2.3. System Characteristic Equation

According to equation (5), the acceleration of the angular response of the floating ring can be expressed as(8)φ¨xtφ¨yt=−cxx+csxIx−cxyIx−cyxIy−cyy+csyIyφ˙xtφ˙yt+−kxx+ksxIx−kxyIx−kyxIy−kyy+ksyIyφxtφyt.Applying Laplace transform to both ends of equation (8) at the same time and using the second differential property, we can get(9)s2Φ⟶s−sφ⟶0−φ˙⟶0=CsΦ⟶s−φ⟶0+KΦ⟶s,where the state variables and coefficient matrices are defined as follows:(10a)φ¨⟶t=φ¨xt,φ¨ytT,φ˙⟶t=φ˙xt,φ˙ytT,φ⟶t=φxt,φytT,Φ⟶s=Φxs,ΦysT,(10b)C=−cxx+csxIx−cxyIx−cyxIy−cyy+csyIy,K=−kxx+ksxIx−kxyIx−kyxIy−kyy+ksyIy,where s is a variable in the complex frequency domain. Φx(s) and Φy(s) are image functions of the angular displacements φx(t) and φy(t), respectively.From equation (9), the transfer function Φ⟶s can be obtained as(11)Φ⟶s=s2E2−sC−K−1sE2−Cφ⟶0+φ˙⟶0,where E2 represents the second-order identity matrix, and the matrices (s2E2 − sC − K) and (sE2 − C) are as follows:(12)s2E2−sC−K=s2+cxx+csxIxs+kxx+ksxIxcxyIxs+kxyIxcyxIys+kyxIys2+cyy+csyIys+kyy+ksyIy,(13)sE2−C=s+cxx+csxIxcxyIxcyxIys+cyy+csyIy.Due to the symmetry of the floating ring and the angular dynamic characteristic coefficient of the gas film, there are(14)Ix=Iy=I,kxx=kyy,cxx=cyy,kxy=−kyx,cxy=−cyx.According to reference [2], there are(15)ksx=ksy=0.5ksrs12,csx=csy=0.5csrs22,where ks is the spring stiffness; cs is the auxiliary seal damping; rs1 and rs2 are the radial installation positions of the spring and the auxiliary sealing ring respectively, and rs1 = 0.5(ro + rb), rs2 = rb; rb is the equilibrium radius, which can be determined by the following formula according to the equilibrium ratio B:(16)rb=ro2−Bro2−ri2.Substituting equations (14) and (15) into equations (12) and (13), we can get(17)s2E2−sC−K=s2+cxx+csxIs+kxx+ksxI−cyxIs−kyxIcyxIs+kyxIs2+cxx+csxIs+kxx+ksxI,(18)sE2−C=s+cxx+csxI−cyxIcyxIs+cxx+csxI.Calculating the inverse matrix of equation (17), we can get(19)s2E2−sC−K−1=1s2+cxx+csx/Is+kxx+ksx/I2+cyx/Is+kyx/I2⋅s2+cxx+csxIs+kxx+ksxIcyxIs+kyxI−cyxIs−kyxIs2+cxx+csxIs+kxx+ksxI.By substituting equations (18), and (19) into (11), the transfer function Φ⟶s, can be transformed into(20)ΦxsΦys=1s2+cxx+csx/Is+kxx+ksx/I2+cyx/Is+kyx/I2⋅s2+cxx+csxIs+kxx+ksxIcyxIs+kyxI−cyxIs−kyxIs2+cxx+csxIs+kxx+ksxI⋅s+cxx+csxIφx0−cyxIφy0+φ˙x0cyxIφx0+s+cxx+csxIφy0+φ˙y0.The two DOF angular vibration response of the floating ring can be obtained by applying the inverse Laplace transform to equation (20). It is worth noting that the more convenient way to judge the stability of the angular free vibration of the system is to use the system’s characteristic equation instead of solving the vibration differential equation completely to obtain the angular vibration mode of the floating ring. In automatic control theory, the root of the characteristic equation is called the closed-loop pole, and each closed-loop pole λi corresponds to a factor eλit in the time response, that is, the closed-loop pole determines the mode of the system response. According to the different positions of the poles on the complex plane, it can be known that the corresponding transient components present as decreasing amplitude oscillation, exponential attenuation, increasing amplitude oscillation, exponential growth, or constant amplitude oscillation. The closed-loop pole distribution determines the motion form of the transient component and reflects the system’s stability. For a stable high-order system, the closed-loop poles are all on the left half-open plane of the complex plane, that is, the roots of the characteristic equation of the system all have negative real parts, and only negative exponents exist in the time response relation of the system. Over time, the transient components will attenuate to zero.Figure2 shows the closed-loop pole distribution and the corresponding transient component motion in the complex plane. It can be seen that the distribution of poles on the left or right side of the imaginary axis determines the final value of the response: when λi is on the left side of the imaginary axis, the transient component finally attenuates to zero; when λi is on the right side of the imaginary axis, the transient component must diverge and tend to infinity; when λi is on the imaginary axis (except the origin), the transient component oscillates with constant amplitude. The imaginary and real distribution determines the mode shapes: when λi is on the real axis, the transient component is aperiodic motion, and when λi is on the complex plane, the transient component is periodic motion. The distance between the pole and the imaginary axis determines the attenuation rate. When λi is on the left side of the imaginary axis, the farther λi is from the imaginary axis, the faster the transition process attenuates.Figure 2
Closed-loop poles distribution and motion form of transient components.Set the denominator of equation (20) to zero, and the characteristic equation of the system can be obtained as follows:(21)s4+2cxx+csxIs3+cxx+csxI2+2kxx+ksxI+cyxI2s2+2cxx+csxkxx+ksxI2+2cyxkyxI2s+kxx+ksxI2+kyxI2=0.


#### 2.2.4. Root Locus Analysis and Algebraic Stability Criterion

The system characteristic equation (21) of two DOF angular free vibration of DGS is a univariate quartic equation. The characteristic root can be solved directly by the Ferrari method, and then the root locus diagram under the influence of a single factor can be drawn. Root locus analysis is a common method for qualitative analysis of system performance, and its advantage is that the root locus diagram can directly reflect dynamic performance such as stability, rapidity, damping coefficient, vibration mode, etc. However, it is difficult to directly solve all roots of the univariate quartic equation in engineering because the Ferrari method is complicated in solving processes and involves hard operations such as extracting the square root of complex numbers. Using the relationship between the root of the system characteristic equation and its coefficient, Hurwitz obtains the method of judging the real part symbol of the characteristic root according to the coefficient of the characteristic equation, thus determining the system’s stability, which is called Hurwitz stability criterion. Hurwitz stability criterion is a simple algebraic criterion method to obtain the necessary and sufficient conditions for system stability without solving higher-order characteristic equations. Hurwitz stability criterion can judge whether the system is absolutely stable, but cannot reflect the degree of system stability and the dynamic performance such as oscillation frequency, amplitude, and adjusting time. It is challenging to analyze the influence of changes in parameters on performance. However, it is an excellent method to quantitatively determine the parameters in the critical instability state of the system.Figure3 shows the method for qualitatively analyzing the performance of a second-order underdamped system through the position of closed-loop poles, in which two closed-loop poles are conjugate, respectively expressed as(22)λ1,2=−ζωn±iωn1−ζ2,where ζ is the damping ratio; ωn is the undamped natural oscillation frequency (the system’s natural frequency).Figure 3
Qualitative analysis of system performance: (a) definition of transient performance parameters and (b) relationship between closed-loop pole position and performance index.
(a)(b)Figure3(a) defines two transient performance parameters, adjusting time ts and overshoot σ%, which can be expressed as follows:(23)ts≥1ζωnln1Δ1−ζ2,σ%=exp−πζ/1−ζ2×100%,where Δ is the allowable fluctuation range.As can be seen from Figure3, the distance between the closed-loop pole and the imaginary axis is ζωn, that is, the real part reflects the adjusting time ts of the system. The larger the absolute value of the real part, the shorter the adjusting time, and the faster the system’s response speed. The value of the imaginary part of the closed-loop pole reflects the oscillation frequency ωd of the output response of the system. The distance between the closed-loop pole and the coordinate origin reflects the undamped natural oscillation frequency ωn of the system. The angle β (cosβ = ζ) between the closed-loop pole and the negative real axis reflects the overshoot σ% of the system. The larger β is, the larger the overshoot is, and the worse the stability of the system response is.According to Hurwitz stability criterion, for the characteristic equation of the system:(24)ansn+an−1sn−1+⋯+a1s+a0=0.The sufficient and necessary condition for the stability of the system is that the coefficients of the characteristic equation of the system are all positive, that is,ai > 0 (i = 0,1,2, …, n), and the leading principal minors of each order Δ1, Δ2, …, Δn of the n-order Hurwitz determinant H composed of various coefficients are all greater than zero. The construction method of the Hurwitz determinant (as shown in equation (25)) is: the items on the principal diagonal are the second coefficient an − 1 to the last coefficient a0 of the characteristic equation; the subscripts of the coefficients increase gradually in the rows below the principal diagonal, and decrease progressively in the rows above the principal diagonal; when the subscript is greater than n or less than 0, the item in the determinant is taken as 0.(25)H=an−1an−3an−5an−7⋯0anan−2an−4an−6⋯00an−1an−3an−5⋯00anan−2an−4⋯0⋮⋮⋮⋮⋮⋮0000⋯a0.In particular, for the fourth-order characteristic equation, according to the above algebraic criterion method, the necessary and sufficient conditions for the stability of the system are as follows:(26)ai>0a3a2a1>a4a12+a32a0,i=0,1,2,3,4.By substituting equation (21) into equation (26), the algebraic criterion for judging the stability of two DOF angular free vibration of DGS can be obtained:(27a)cxx+csx>0,(27b)kxx+ksx>−cxx+csx2+cyx22I,(27c)cyxkyx>−cxx+csxkxx+ksx,(27d)I<cxx+csx4cxx+csx2+cyx2kxx+ksxkyx2+cxx+csx2cyx2cxx+csx2+cyx2kxx+ksxkyx2+cxx+csx3cyxcxx+csx2+cyx21kyx+cxx+csxcyx3cxx+csx2+cyx21kyx.It can be seen from equations (27a)–(27d) that, different from the axial free vibration stability of DGS, the stability of two DOF angular free vibration of DGS is not only related to the angular main damping cxx + csx, but jointly determined by five factors: main damping cxx + csx, main stiffness kxx + ksx, cross damping cyx, cross stiffness kyx and rotational inertia I. If any of the four conditions in equations (27a)–(27d) is not satisfied, the floating ring of DGS will lose its stability under any instantaneous angular excitation. Furthermore, equations (27a)–(27d) are used to determine whether the main damping, the main stiffness, the cross coefficient, and the rotational inertia satisfy the stability conditions of angular vibration, respectively.In previous studies on the film dynamic characteristic coefficient of DGS, only the influence of angular main damping and main stiffness is often concerned. Still, the method, which judges the angular gas film stability only by whether the main damping is greater than zero, ignores the influence of cross coefficient. According to equation (27c), the cross coefficient will also affect the critical value of the angular main damping and the main stiffness. It is not true that the system must be stable when the main damping is greater than zero. Due to symmetry, there will be positive and negative values of the cross stiffness and the cross damping, and few scholars have clearly proposed a reasonable range of cross coefficient. However, the uncertainty of the physical meaning and influence of the cross coefficient will lead to errors in judging the angular vibration stability of DGS. In addition, some scholars have put forward the concept of critical rotational inertia in the study of angular free vibration stability of DGS. However, the influence law of rotational inertia on stability and its interaction with gas film dynamic characteristic parameters have not been discussed in depth.


## 2.1. Physical and Geometric Model

Figure1(a) shows the geometric structure model of the flexibly mounted stator DGS. At the sealing face of the rotor near the outer radius ro, there are circumferential uniformly distributed logarithmic spiral grooves with a quantity of Ng and a depth of hg, and the spiral angle is α; an ungrooved sealing dam is arranged adjacent to the inner radius ri to achieve sealing during shutdown; the boundary radius between the sealing dam and the sealing groove is rg; the equilibrium film thickness at the sealing dam is hb. Due to installation deviation or shafting vibration or swing, the rotor fixed on the rotating shaft and rotating at an angular velocity Ω has a certain axial and angular motion, and applies continuous axial and angular excitation to the floating stator through the sealing gas film, in which the axial excitation amplitude and frequency are Arz and ω respectively, and the fixed installation angle deviation of the rotor is Ar.Figure 1
Schematic diagram of the geometric structure and kinematic model of the flexibly mounted stator DGS: (a) geometric structure model; (b) kinematic model and (c) positional relationship.
(a)(b)(c)The sealing gas film between rotor and stator can be regarded as a “spring-damping” structure with certain stiffness and damping. The floating stator can respond to the three DOF excitation of the rotor and keep following the disturbance of the rotor under the joint action of the thin film of gas, spring, and auxiliary sealing ring. Because the axial and angular motion can be decoupled [30, 31], only the angular disturbance is considered in this analysis.During seal transient, the transient gas film thickness is determined by the relative position of the seal rings. Assuming only rigid body motion, the position of each face can be completely described by their angular orientation. As shown in Figure1(c), the inertial coordinate system Oxyz can be considered as fixed in space. z is the axis about which the shaft rotates. Osxsyszs is the stator reference coordinate system, where axis xs is parallel to the Oxy plane. Axis ys points to the point in the stator plane of maximum distance from the Oxy plane, and zs is the central axis of the stator.


## 2.2. Mathematic Model

### 2.2.1. Angular Free Vibration Model

Considering the stiffness of the supporting spring and the damping of the auxiliary seal, and ignoring the influence of the axial disturbance on the angular motion response, according to Newton’s law of motion, the two DOF angular free vibration motion equation of the floating ring can be obtained as follows:(1)Ixφ¨x=Mx′−Mxsp′−Mxse′,Iyφ¨y=My′−Mysp′−Myse′,where Ix and Iy are the rotational inertia of the floating ring around the x-axis and y-axis, respectively; φ¨x and φ¨y are the angular acceleration of the response of the floating ring around the x-axis and y-axis, respectively; Mx′ and My′ are the x-direction and y-direction deflection moment increment of the gas film under perturbation respectively, which can be obtained from the Taylor expansion and the definition of the dynamic characteristic coefficient of the gas film:(2)Mx′=−kxxφx−kxyφy−cxxφ˙x−cxyφ˙y,My′=−kyxφx−kyyφy−cyxφ˙x−cyyφ˙y,where φx, φy, φ˙x and φ˙y are the angular displacement and angular velocity of the angular response of the floating ring, respectively, and the subscript indicates the direction of the corresponding physical quantities; kij and cij (i, j = x, y) are the angular dynamic stiffness and damping coefficients of the gas film, respectively. The first subscript indicates the direction of the moment. The second subscript indicates the direction of the disturbing displacement or velocity that causes the moment change. The negative sign indicates that the increment of the moment always prevents the displacement and velocity of the disturbance.Mxsp′ and Mysp′ are the deflection recovery moment of the supporting spring to the floating ring caused by the displacement perturbation in the x and y directions, respectively. If ksx and ksy are used to express the equivalent recovery moment stiffness of the spring system in the corresponding direction under the angular displacement disturbance φx and φy respectively, then(3)Mxsp′=ksxφx,Mysp′=ksyφy.Mxse′ and Myse′ are the deflection recovery moment of the auxiliary seal to the floating ring caused by the velocity perturbation in the x and y directions, respectively. If csx and csy are used to express the equivalent recovery moment damping of the auxiliary seal in the corresponding direction under the angular velocity disturbance φ˙x and φ˙y respectively, then(4)Mxse′=csxφ˙x,Myse′=csyφ˙y.By substituting equations (2)–(4) into equation (1), the equation of motion of the angular free vibration of the floating ring can be transformed into(5)Ixφ¨xt+cxx+csxφ˙xt+cxyφ˙yt+kxx+ksxφxt+kxyφyt=0,Iyφ¨yt+cyxφ˙xt+cyy+csyφ˙yt+kyxφxt+kyy+ksyφyt=0.


### 2.2.2. Film Dynamic Coefficients

In this paper, the perturbation method is used to calculate the linearized angular stiffness and damping coefficients appearing in equation (2). The compressible perturbed Reynolds equations are as follows (the theoretical derivation refers to references [31, 32]):(6a)1r∂∂r3r2h02sinθ2μ∂p02∂r+rh03μ∂p0pαr∂r+1r2∂∂θ3rh02sinθ2μ∂p02∂θ+h03μ∂p0pαr∂θ−6Ω∂pαrh0+p0rsinθ∂θ+12ωpαih0=0,1r∂∂rrh03μ∂p0pαi∂r+1r2∂∂θh03μ∂p0pαi∂θ−6Ω∂pαih0∂θ−12ωpαrh0+p0rsinθ=0,(6b)1r∂∂rrh03μ∂p0pβr∂r−3r2h02cosθ2μ∂p02∂r+1r2∂∂θh03μ∂p0pβr∂θ−3rh02cosθ2μ∂p02∂θ−6Ω∂pβrh0−p0rcosθ∂θ+12ωpβih0=0,1r∂∂rrh03μ∂p0pβi∂r+1r2∂∂θh03μ∂p0pβi∂θ−6Ω∂pβih0∂θ−12ωpβrh0−p0rcosθ=0,where p0 and h0 are the steady-state gas film pressure and thickness, respectively. pkj (k = α,β; j = r,i) is the first order perturbed pressure. The first subscript represents the direction of the perturbed pressure, and the second subscript represents the real part or imaginary part. ω is the frequency of angular excitation motion. In general cases (such as the angular excitation introduced by installation deviations of seal rings), the excitation frequency ω is equal to the rotational angular velocity Ω.Solving equation (6) to obtain gas film perturbed pressure distribution, and eight dynamic stiffness and damping coefficients can be obtained by integrating on the whole seal face according to the following formula:(7a)kxxkxykyxkyy=−∫02π∫riropαrrsinθpβrrsinθ−pαrrcosθ−pβrrcosθrdrdθ,(7b)cxxcxycyxcyy=−∫02π∫riropαirsinθωpβirsinθω−pαircosθω−pβircosθωrdrdθ.


### 2.2.3. System Characteristic Equation

According to equation (5), the acceleration of the angular response of the floating ring can be expressed as(8)φ¨xtφ¨yt=−cxx+csxIx−cxyIx−cyxIy−cyy+csyIyφ˙xtφ˙yt+−kxx+ksxIx−kxyIx−kyxIy−kyy+ksyIyφxtφyt.Applying Laplace transform to both ends of equation (8) at the same time and using the second differential property, we can get(9)s2Φ⟶s−sφ⟶0−φ˙⟶0=CsΦ⟶s−φ⟶0+KΦ⟶s,where the state variables and coefficient matrices are defined as follows:(10a)φ¨⟶t=φ¨xt,φ¨ytT,φ˙⟶t=φ˙xt,φ˙ytT,φ⟶t=φxt,φytT,Φ⟶s=Φxs,ΦysT,(10b)C=−cxx+csxIx−cxyIx−cyxIy−cyy+csyIy,K=−kxx+ksxIx−kxyIx−kyxIy−kyy+ksyIy,where s is a variable in the complex frequency domain. Φx(s) and Φy(s) are image functions of the angular displacements φx(t) and φy(t), respectively.From equation (9), the transfer function Φ⟶s can be obtained as(11)Φ⟶s=s2E2−sC−K−1sE2−Cφ⟶0+φ˙⟶0,where E2 represents the second-order identity matrix, and the matrices (s2E2 − sC − K) and (sE2 − C) are as follows:(12)s2E2−sC−K=s2+cxx+csxIxs+kxx+ksxIxcxyIxs+kxyIxcyxIys+kyxIys2+cyy+csyIys+kyy+ksyIy,(13)sE2−C=s+cxx+csxIxcxyIxcyxIys+cyy+csyIy.Due to the symmetry of the floating ring and the angular dynamic characteristic coefficient of the gas film, there are(14)Ix=Iy=I,kxx=kyy,cxx=cyy,kxy=−kyx,cxy=−cyx.According to reference [2], there are(15)ksx=ksy=0.5ksrs12,csx=csy=0.5csrs22,where ks is the spring stiffness; cs is the auxiliary seal damping; rs1 and rs2 are the radial installation positions of the spring and the auxiliary sealing ring respectively, and rs1 = 0.5(ro + rb), rs2 = rb; rb is the equilibrium radius, which can be determined by the following formula according to the equilibrium ratio B:(16)rb=ro2−Bro2−ri2.Substituting equations (14) and (15) into equations (12) and (13), we can get(17)s2E2−sC−K=s2+cxx+csxIs+kxx+ksxI−cyxIs−kyxIcyxIs+kyxIs2+cxx+csxIs+kxx+ksxI,(18)sE2−C=s+cxx+csxI−cyxIcyxIs+cxx+csxI.Calculating the inverse matrix of equation (17), we can get(19)s2E2−sC−K−1=1s2+cxx+csx/Is+kxx+ksx/I2+cyx/Is+kyx/I2⋅s2+cxx+csxIs+kxx+ksxIcyxIs+kyxI−cyxIs−kyxIs2+cxx+csxIs+kxx+ksxI.By substituting equations (18), and (19) into (11), the transfer function Φ⟶s, can be transformed into(20)ΦxsΦys=1s2+cxx+csx/Is+kxx+ksx/I2+cyx/Is+kyx/I2⋅s2+cxx+csxIs+kxx+ksxIcyxIs+kyxI−cyxIs−kyxIs2+cxx+csxIs+kxx+ksxI⋅s+cxx+csxIφx0−cyxIφy0+φ˙x0cyxIφx0+s+cxx+csxIφy0+φ˙y0.The two DOF angular vibration response of the floating ring can be obtained by applying the inverse Laplace transform to equation (20). It is worth noting that the more convenient way to judge the stability of the angular free vibration of the system is to use the system’s characteristic equation instead of solving the vibration differential equation completely to obtain the angular vibration mode of the floating ring. In automatic control theory, the root of the characteristic equation is called the closed-loop pole, and each closed-loop pole λi corresponds to a factor eλit in the time response, that is, the closed-loop pole determines the mode of the system response. According to the different positions of the poles on the complex plane, it can be known that the corresponding transient components present as decreasing amplitude oscillation, exponential attenuation, increasing amplitude oscillation, exponential growth, or constant amplitude oscillation. The closed-loop pole distribution determines the motion form of the transient component and reflects the system’s stability. For a stable high-order system, the closed-loop poles are all on the left half-open plane of the complex plane, that is, the roots of the characteristic equation of the system all have negative real parts, and only negative exponents exist in the time response relation of the system. Over time, the transient components will attenuate to zero.Figure2 shows the closed-loop pole distribution and the corresponding transient component motion in the complex plane. It can be seen that the distribution of poles on the left or right side of the imaginary axis determines the final value of the response: when λi is on the left side of the imaginary axis, the transient component finally attenuates to zero; when λi is on the right side of the imaginary axis, the transient component must diverge and tend to infinity; when λi is on the imaginary axis (except the origin), the transient component oscillates with constant amplitude. The imaginary and real distribution determines the mode shapes: when λi is on the real axis, the transient component is aperiodic motion, and when λi is on the complex plane, the transient component is periodic motion. The distance between the pole and the imaginary axis determines the attenuation rate. When λi is on the left side of the imaginary axis, the farther λi is from the imaginary axis, the faster the transition process attenuates.Figure 2
Closed-loop poles distribution and motion form of transient components.Set the denominator of equation (20) to zero, and the characteristic equation of the system can be obtained as follows:(21)s4+2cxx+csxIs3+cxx+csxI2+2kxx+ksxI+cyxI2s2+2cxx+csxkxx+ksxI2+2cyxkyxI2s+kxx+ksxI2+kyxI2=0.


### 2.2.4. Root Locus Analysis and Algebraic Stability Criterion

The system characteristic equation (21) of two DOF angular free vibration of DGS is a univariate quartic equation. The characteristic root can be solved directly by the Ferrari method, and then the root locus diagram under the influence of a single factor can be drawn. Root locus analysis is a common method for qualitative analysis of system performance, and its advantage is that the root locus diagram can directly reflect dynamic performance such as stability, rapidity, damping coefficient, vibration mode, etc. However, it is difficult to directly solve all roots of the univariate quartic equation in engineering because the Ferrari method is complicated in solving processes and involves hard operations such as extracting the square root of complex numbers. Using the relationship between the root of the system characteristic equation and its coefficient, Hurwitz obtains the method of judging the real part symbol of the characteristic root according to the coefficient of the characteristic equation, thus determining the system’s stability, which is called Hurwitz stability criterion. Hurwitz stability criterion is a simple algebraic criterion method to obtain the necessary and sufficient conditions for system stability without solving higher-order characteristic equations. Hurwitz stability criterion can judge whether the system is absolutely stable, but cannot reflect the degree of system stability and the dynamic performance such as oscillation frequency, amplitude, and adjusting time. It is challenging to analyze the influence of changes in parameters on performance. However, it is an excellent method to quantitatively determine the parameters in the critical instability state of the system.Figure3 shows the method for qualitatively analyzing the performance of a second-order underdamped system through the position of closed-loop poles, in which two closed-loop poles are conjugate, respectively expressed as(22)λ1,2=−ζωn±iωn1−ζ2,where ζ is the damping ratio; ωn is the undamped natural oscillation frequency (the system’s natural frequency).Figure 3
Qualitative analysis of system performance: (a) definition of transient performance parameters and (b) relationship between closed-loop pole position and performance index.
(a)(b)Figure3(a) defines two transient performance parameters, adjusting time ts and overshoot σ%, which can be expressed as follows:(23)ts≥1ζωnln1Δ1−ζ2,σ%=exp−πζ/1−ζ2×100%,where Δ is the allowable fluctuation range.As can be seen from Figure3, the distance between the closed-loop pole and the imaginary axis is ζωn, that is, the real part reflects the adjusting time ts of the system. The larger the absolute value of the real part, the shorter the adjusting time, and the faster the system’s response speed. The value of the imaginary part of the closed-loop pole reflects the oscillation frequency ωd of the output response of the system. The distance between the closed-loop pole and the coordinate origin reflects the undamped natural oscillation frequency ωn of the system. The angle β (cosβ = ζ) between the closed-loop pole and the negative real axis reflects the overshoot σ% of the system. The larger β is, the larger the overshoot is, and the worse the stability of the system response is.According to Hurwitz stability criterion, for the characteristic equation of the system:(24)ansn+an−1sn−1+⋯+a1s+a0=0.The sufficient and necessary condition for the stability of the system is that the coefficients of the characteristic equation of the system are all positive, that is,ai > 0 (i = 0,1,2, …, n), and the leading principal minors of each order Δ1, Δ2, …, Δn of the n-order Hurwitz determinant H composed of various coefficients are all greater than zero. The construction method of the Hurwitz determinant (as shown in equation (25)) is: the items on the principal diagonal are the second coefficient an − 1 to the last coefficient a0 of the characteristic equation; the subscripts of the coefficients increase gradually in the rows below the principal diagonal, and decrease progressively in the rows above the principal diagonal; when the subscript is greater than n or less than 0, the item in the determinant is taken as 0.(25)H=an−1an−3an−5an−7⋯0anan−2an−4an−6⋯00an−1an−3an−5⋯00anan−2an−4⋯0⋮⋮⋮⋮⋮⋮0000⋯a0.In particular, for the fourth-order characteristic equation, according to the above algebraic criterion method, the necessary and sufficient conditions for the stability of the system are as follows:(26)ai>0a3a2a1>a4a12+a32a0,i=0,1,2,3,4.By substituting equation (21) into equation (26), the algebraic criterion for judging the stability of two DOF angular free vibration of DGS can be obtained:(27a)cxx+csx>0,(27b)kxx+ksx>−cxx+csx2+cyx22I,(27c)cyxkyx>−cxx+csxkxx+ksx,(27d)I<cxx+csx4cxx+csx2+cyx2kxx+ksxkyx2+cxx+csx2cyx2cxx+csx2+cyx2kxx+ksxkyx2+cxx+csx3cyxcxx+csx2+cyx21kyx+cxx+csxcyx3cxx+csx2+cyx21kyx.It can be seen from equations (27a)–(27d) that, different from the axial free vibration stability of DGS, the stability of two DOF angular free vibration of DGS is not only related to the angular main damping cxx + csx, but jointly determined by five factors: main damping cxx + csx, main stiffness kxx + ksx, cross damping cyx, cross stiffness kyx and rotational inertia I. If any of the four conditions in equations (27a)–(27d) is not satisfied, the floating ring of DGS will lose its stability under any instantaneous angular excitation. Furthermore, equations (27a)–(27d) are used to determine whether the main damping, the main stiffness, the cross coefficient, and the rotational inertia satisfy the stability conditions of angular vibration, respectively.In previous studies on the film dynamic characteristic coefficient of DGS, only the influence of angular main damping and main stiffness is often concerned. Still, the method, which judges the angular gas film stability only by whether the main damping is greater than zero, ignores the influence of cross coefficient. According to equation (27c), the cross coefficient will also affect the critical value of the angular main damping and the main stiffness. It is not true that the system must be stable when the main damping is greater than zero. Due to symmetry, there will be positive and negative values of the cross stiffness and the cross damping, and few scholars have clearly proposed a reasonable range of cross coefficient. However, the uncertainty of the physical meaning and influence of the cross coefficient will lead to errors in judging the angular vibration stability of DGS. In addition, some scholars have put forward the concept of critical rotational inertia in the study of angular free vibration stability of DGS. However, the influence law of rotational inertia on stability and its interaction with gas film dynamic characteristic parameters have not been discussed in depth.


## 2.2.1. Angular Free Vibration Model

Considering the stiffness of the supporting spring and the damping of the auxiliary seal, and ignoring the influence of the axial disturbance on the angular motion response, according to Newton’s law of motion, the two DOF angular free vibration motion equation of the floating ring can be obtained as follows:(1)Ixφ¨x=Mx′−Mxsp′−Mxse′,Iyφ¨y=My′−Mysp′−Myse′,where Ix and Iy are the rotational inertia of the floating ring around the x-axis and y-axis, respectively; φ¨x and φ¨y are the angular acceleration of the response of the floating ring around the x-axis and y-axis, respectively; Mx′ and My′ are the x-direction and y-direction deflection moment increment of the gas film under perturbation respectively, which can be obtained from the Taylor expansion and the definition of the dynamic characteristic coefficient of the gas film:(2)Mx′=−kxxφx−kxyφy−cxxφ˙x−cxyφ˙y,My′=−kyxφx−kyyφy−cyxφ˙x−cyyφ˙y,where φx, φy, φ˙x and φ˙y are the angular displacement and angular velocity of the angular response of the floating ring, respectively, and the subscript indicates the direction of the corresponding physical quantities; kij and cij (i, j = x, y) are the angular dynamic stiffness and damping coefficients of the gas film, respectively. The first subscript indicates the direction of the moment. The second subscript indicates the direction of the disturbing displacement or velocity that causes the moment change. The negative sign indicates that the increment of the moment always prevents the displacement and velocity of the disturbance.Mxsp′ and Mysp′ are the deflection recovery moment of the supporting spring to the floating ring caused by the displacement perturbation in the x and y directions, respectively. If ksx and ksy are used to express the equivalent recovery moment stiffness of the spring system in the corresponding direction under the angular displacement disturbance φx and φy respectively, then(3)Mxsp′=ksxφx,Mysp′=ksyφy.Mxse′ and Myse′ are the deflection recovery moment of the auxiliary seal to the floating ring caused by the velocity perturbation in the x and y directions, respectively. If csx and csy are used to express the equivalent recovery moment damping of the auxiliary seal in the corresponding direction under the angular velocity disturbance φ˙x and φ˙y respectively, then(4)Mxse′=csxφ˙x,Myse′=csyφ˙y.By substituting equations (2)–(4) into equation (1), the equation of motion of the angular free vibration of the floating ring can be transformed into(5)Ixφ¨xt+cxx+csxφ˙xt+cxyφ˙yt+kxx+ksxφxt+kxyφyt=0,Iyφ¨yt+cyxφ˙xt+cyy+csyφ˙yt+kyxφxt+kyy+ksyφyt=0.


## 2.2.2. Film Dynamic Coefficients

In this paper, the perturbation method is used to calculate the linearized angular stiffness and damping coefficients appearing in equation (2). The compressible perturbed Reynolds equations are as follows (the theoretical derivation refers to references [31, 32]):(6a)1r∂∂r3r2h02sinθ2μ∂p02∂r+rh03μ∂p0pαr∂r+1r2∂∂θ3rh02sinθ2μ∂p02∂θ+h03μ∂p0pαr∂θ−6Ω∂pαrh0+p0rsinθ∂θ+12ωpαih0=0,1r∂∂rrh03μ∂p0pαi∂r+1r2∂∂θh03μ∂p0pαi∂θ−6Ω∂pαih0∂θ−12ωpαrh0+p0rsinθ=0,(6b)1r∂∂rrh03μ∂p0pβr∂r−3r2h02cosθ2μ∂p02∂r+1r2∂∂θh03μ∂p0pβr∂θ−3rh02cosθ2μ∂p02∂θ−6Ω∂pβrh0−p0rcosθ∂θ+12ωpβih0=0,1r∂∂rrh03μ∂p0pβi∂r+1r2∂∂θh03μ∂p0pβi∂θ−6Ω∂pβih0∂θ−12ωpβrh0−p0rcosθ=0,where p0 and h0 are the steady-state gas film pressure and thickness, respectively. pkj (k = α,β; j = r,i) is the first order perturbed pressure. The first subscript represents the direction of the perturbed pressure, and the second subscript represents the real part or imaginary part. ω is the frequency of angular excitation motion. In general cases (such as the angular excitation introduced by installation deviations of seal rings), the excitation frequency ω is equal to the rotational angular velocity Ω.Solving equation (6) to obtain gas film perturbed pressure distribution, and eight dynamic stiffness and damping coefficients can be obtained by integrating on the whole seal face according to the following formula:(7a)kxxkxykyxkyy=−∫02π∫riropαrrsinθpβrrsinθ−pαrrcosθ−pβrrcosθrdrdθ,(7b)cxxcxycyxcyy=−∫02π∫riropαirsinθωpβirsinθω−pαircosθω−pβircosθωrdrdθ.


## 2.2.3. System Characteristic Equation

According to equation (5), the acceleration of the angular response of the floating ring can be expressed as(8)φ¨xtφ¨yt=−cxx+csxIx−cxyIx−cyxIy−cyy+csyIyφ˙xtφ˙yt+−kxx+ksxIx−kxyIx−kyxIy−kyy+ksyIyφxtφyt.Applying Laplace transform to both ends of equation (8) at the same time and using the second differential property, we can get(9)s2Φ⟶s−sφ⟶0−φ˙⟶0=CsΦ⟶s−φ⟶0+KΦ⟶s,where the state variables and coefficient matrices are defined as follows:(10a)φ¨⟶t=φ¨xt,φ¨ytT,φ˙⟶t=φ˙xt,φ˙ytT,φ⟶t=φxt,φytT,Φ⟶s=Φxs,ΦysT,(10b)C=−cxx+csxIx−cxyIx−cyxIy−cyy+csyIy,K=−kxx+ksxIx−kxyIx−kyxIy−kyy+ksyIy,where s is a variable in the complex frequency domain. Φx(s) and Φy(s) are image functions of the angular displacements φx(t) and φy(t), respectively.From equation (9), the transfer function Φ⟶s can be obtained as(11)Φ⟶s=s2E2−sC−K−1sE2−Cφ⟶0+φ˙⟶0,where E2 represents the second-order identity matrix, and the matrices (s2E2 − sC − K) and (sE2 − C) are as follows:(12)s2E2−sC−K=s2+cxx+csxIxs+kxx+ksxIxcxyIxs+kxyIxcyxIys+kyxIys2+cyy+csyIys+kyy+ksyIy,(13)sE2−C=s+cxx+csxIxcxyIxcyxIys+cyy+csyIy.Due to the symmetry of the floating ring and the angular dynamic characteristic coefficient of the gas film, there are(14)Ix=Iy=I,kxx=kyy,cxx=cyy,kxy=−kyx,cxy=−cyx.According to reference [2], there are(15)ksx=ksy=0.5ksrs12,csx=csy=0.5csrs22,where ks is the spring stiffness; cs is the auxiliary seal damping; rs1 and rs2 are the radial installation positions of the spring and the auxiliary sealing ring respectively, and rs1 = 0.5(ro + rb), rs2 = rb; rb is the equilibrium radius, which can be determined by the following formula according to the equilibrium ratio B:(16)rb=ro2−Bro2−ri2.Substituting equations (14) and (15) into equations (12) and (13), we can get(17)s2E2−sC−K=s2+cxx+csxIs+kxx+ksxI−cyxIs−kyxIcyxIs+kyxIs2+cxx+csxIs+kxx+ksxI,(18)sE2−C=s+cxx+csxI−cyxIcyxIs+cxx+csxI.Calculating the inverse matrix of equation (17), we can get(19)s2E2−sC−K−1=1s2+cxx+csx/Is+kxx+ksx/I2+cyx/Is+kyx/I2⋅s2+cxx+csxIs+kxx+ksxIcyxIs+kyxI−cyxIs−kyxIs2+cxx+csxIs+kxx+ksxI.By substituting equations (18), and (19) into (11), the transfer function Φ⟶s, can be transformed into(20)ΦxsΦys=1s2+cxx+csx/Is+kxx+ksx/I2+cyx/Is+kyx/I2⋅s2+cxx+csxIs+kxx+ksxIcyxIs+kyxI−cyxIs−kyxIs2+cxx+csxIs+kxx+ksxI⋅s+cxx+csxIφx0−cyxIφy0+φ˙x0cyxIφx0+s+cxx+csxIφy0+φ˙y0.The two DOF angular vibration response of the floating ring can be obtained by applying the inverse Laplace transform to equation (20). It is worth noting that the more convenient way to judge the stability of the angular free vibration of the system is to use the system’s characteristic equation instead of solving the vibration differential equation completely to obtain the angular vibration mode of the floating ring. In automatic control theory, the root of the characteristic equation is called the closed-loop pole, and each closed-loop pole λi corresponds to a factor eλit in the time response, that is, the closed-loop pole determines the mode of the system response. According to the different positions of the poles on the complex plane, it can be known that the corresponding transient components present as decreasing amplitude oscillation, exponential attenuation, increasing amplitude oscillation, exponential growth, or constant amplitude oscillation. The closed-loop pole distribution determines the motion form of the transient component and reflects the system’s stability. For a stable high-order system, the closed-loop poles are all on the left half-open plane of the complex plane, that is, the roots of the characteristic equation of the system all have negative real parts, and only negative exponents exist in the time response relation of the system. Over time, the transient components will attenuate to zero.Figure2 shows the closed-loop pole distribution and the corresponding transient component motion in the complex plane. It can be seen that the distribution of poles on the left or right side of the imaginary axis determines the final value of the response: when λi is on the left side of the imaginary axis, the transient component finally attenuates to zero; when λi is on the right side of the imaginary axis, the transient component must diverge and tend to infinity; when λi is on the imaginary axis (except the origin), the transient component oscillates with constant amplitude. The imaginary and real distribution determines the mode shapes: when λi is on the real axis, the transient component is aperiodic motion, and when λi is on the complex plane, the transient component is periodic motion. The distance between the pole and the imaginary axis determines the attenuation rate. When λi is on the left side of the imaginary axis, the farther λi is from the imaginary axis, the faster the transition process attenuates.Figure 2
Closed-loop poles distribution and motion form of transient components.Set the denominator of equation (20) to zero, and the characteristic equation of the system can be obtained as follows:(21)s4+2cxx+csxIs3+cxx+csxI2+2kxx+ksxI+cyxI2s2+2cxx+csxkxx+ksxI2+2cyxkyxI2s+kxx+ksxI2+kyxI2=0.


## 2.2.4. Root Locus Analysis and Algebraic Stability Criterion

The system characteristic equation (21) of two DOF angular free vibration of DGS is a univariate quartic equation. The characteristic root can be solved directly by the Ferrari method, and then the root locus diagram under the influence of a single factor can be drawn. Root locus analysis is a common method for qualitative analysis of system performance, and its advantage is that the root locus diagram can directly reflect dynamic performance such as stability, rapidity, damping coefficient, vibration mode, etc. However, it is difficult to directly solve all roots of the univariate quartic equation in engineering because the Ferrari method is complicated in solving processes and involves hard operations such as extracting the square root of complex numbers. Using the relationship between the root of the system characteristic equation and its coefficient, Hurwitz obtains the method of judging the real part symbol of the characteristic root according to the coefficient of the characteristic equation, thus determining the system’s stability, which is called Hurwitz stability criterion. Hurwitz stability criterion is a simple algebraic criterion method to obtain the necessary and sufficient conditions for system stability without solving higher-order characteristic equations. Hurwitz stability criterion can judge whether the system is absolutely stable, but cannot reflect the degree of system stability and the dynamic performance such as oscillation frequency, amplitude, and adjusting time. It is challenging to analyze the influence of changes in parameters on performance. However, it is an excellent method to quantitatively determine the parameters in the critical instability state of the system.Figure3 shows the method for qualitatively analyzing the performance of a second-order underdamped system through the position of closed-loop poles, in which two closed-loop poles are conjugate, respectively expressed as(22)λ1,2=−ζωn±iωn1−ζ2,where ζ is the damping ratio; ωn is the undamped natural oscillation frequency (the system’s natural frequency).Figure 3
Qualitative analysis of system performance: (a) definition of transient performance parameters and (b) relationship between closed-loop pole position and performance index.
(a)(b)Figure3(a) defines two transient performance parameters, adjusting time ts and overshoot σ%, which can be expressed as follows:(23)ts≥1ζωnln1Δ1−ζ2,σ%=exp−πζ/1−ζ2×100%,where Δ is the allowable fluctuation range.As can be seen from Figure3, the distance between the closed-loop pole and the imaginary axis is ζωn, that is, the real part reflects the adjusting time ts of the system. The larger the absolute value of the real part, the shorter the adjusting time, and the faster the system’s response speed. The value of the imaginary part of the closed-loop pole reflects the oscillation frequency ωd of the output response of the system. The distance between the closed-loop pole and the coordinate origin reflects the undamped natural oscillation frequency ωn of the system. The angle β (cosβ = ζ) between the closed-loop pole and the negative real axis reflects the overshoot σ% of the system. The larger β is, the larger the overshoot is, and the worse the stability of the system response is.According to Hurwitz stability criterion, for the characteristic equation of the system:(24)ansn+an−1sn−1+⋯+a1s+a0=0.The sufficient and necessary condition for the stability of the system is that the coefficients of the characteristic equation of the system are all positive, that is,ai > 0 (i = 0,1,2, …, n), and the leading principal minors of each order Δ1, Δ2, …, Δn of the n-order Hurwitz determinant H composed of various coefficients are all greater than zero. The construction method of the Hurwitz determinant (as shown in equation (25)) is: the items on the principal diagonal are the second coefficient an − 1 to the last coefficient a0 of the characteristic equation; the subscripts of the coefficients increase gradually in the rows below the principal diagonal, and decrease progressively in the rows above the principal diagonal; when the subscript is greater than n or less than 0, the item in the determinant is taken as 0.(25)H=an−1an−3an−5an−7⋯0anan−2an−4an−6⋯00an−1an−3an−5⋯00anan−2an−4⋯0⋮⋮⋮⋮⋮⋮0000⋯a0.In particular, for the fourth-order characteristic equation, according to the above algebraic criterion method, the necessary and sufficient conditions for the stability of the system are as follows:(26)ai>0a3a2a1>a4a12+a32a0,i=0,1,2,3,4.By substituting equation (21) into equation (26), the algebraic criterion for judging the stability of two DOF angular free vibration of DGS can be obtained:(27a)cxx+csx>0,(27b)kxx+ksx>−cxx+csx2+cyx22I,(27c)cyxkyx>−cxx+csxkxx+ksx,(27d)I<cxx+csx4cxx+csx2+cyx2kxx+ksxkyx2+cxx+csx2cyx2cxx+csx2+cyx2kxx+ksxkyx2+cxx+csx3cyxcxx+csx2+cyx21kyx+cxx+csxcyx3cxx+csx2+cyx21kyx.It can be seen from equations (27a)–(27d) that, different from the axial free vibration stability of DGS, the stability of two DOF angular free vibration of DGS is not only related to the angular main damping cxx + csx, but jointly determined by five factors: main damping cxx + csx, main stiffness kxx + ksx, cross damping cyx, cross stiffness kyx and rotational inertia I. If any of the four conditions in equations (27a)–(27d) is not satisfied, the floating ring of DGS will lose its stability under any instantaneous angular excitation. Furthermore, equations (27a)–(27d) are used to determine whether the main damping, the main stiffness, the cross coefficient, and the rotational inertia satisfy the stability conditions of angular vibration, respectively.In previous studies on the film dynamic characteristic coefficient of DGS, only the influence of angular main damping and main stiffness is often concerned. Still, the method, which judges the angular gas film stability only by whether the main damping is greater than zero, ignores the influence of cross coefficient. According to equation (27c), the cross coefficient will also affect the critical value of the angular main damping and the main stiffness. It is not true that the system must be stable when the main damping is greater than zero. Due to symmetry, there will be positive and negative values of the cross stiffness and the cross damping, and few scholars have clearly proposed a reasonable range of cross coefficient. However, the uncertainty of the physical meaning and influence of the cross coefficient will lead to errors in judging the angular vibration stability of DGS. In addition, some scholars have put forward the concept of critical rotational inertia in the study of angular free vibration stability of DGS. However, the influence law of rotational inertia on stability and its interaction with gas film dynamic characteristic parameters have not been discussed in depth.


## 3. Results Analysis and Discussion

Based on the system characteristic equation of two DOF angular free vibration of DGS derived above and Hurwitz stability criterion, the root locus method and the concept of the closed-loop dominant pole are used to analyze the influence of rotational inertia, stiffness, and damping parameters on the angular free vibration stability of DGS. In the calculation examples in this paper, unless otherwise specified, the following parameters are used: transverse rotational inertia of the flexibly mounted statorI = 2.57 × 10−3 kg·m2, spring stiffness ks = 1 × 105 N/m, auxiliary seal damping cs = 200 N·s/m, equilibrium ratio B = 0.6, film angular main damping cxx = 165.2039 N·m·s/rad, film angular main stiffness kxx = 1.7755 × 105 N·m/rad, film angular cross damping cyx = 33.4574 N·m·s/rad, film angular cross stiffness kyx = −1.1824 × 105 N·m/rad.


### 3.1. Influence of Parameters on Angular Free Vibration Stability of DGS

According to the automatic control theory, for a stable high-order system, the distance between the closed-loop pole and the imaginary axis determines the attenuation speed of the corresponding transient component, and the farther the distance, the faster the attenuation. The closed-loop zero-point (the point satisfying the molecular of the transfer function to be zero) affects the relative intensity of each transient component. If there is a zero-point near a certain pole, the intensity of its corresponding transient component will become smaller, and the transient response component of a pair of zero-point and pole that are close to each other can be ignored. Suppose among all the closed-loop poles, there is no closed-loop zero-point around the pole closest to the imaginary axis, and other closed-loop poles are far away from the imaginary axis. In that case, the corresponding transient component of the closed-loop pole closest to the imaginary axis in the unit step response not only has the largest initial value whent = 0, but also attenuates the slowest in all responses. Therefore, it plays an important role in the time response of the system. Such a closed-loop pole is called the closed-loop dominant pole, and other poles are called nondominant poles. The dominant poles are generally a pair of complex conjugate poles, which are located on the left side of the imaginary axis, in the sector with an angle of 60° to the real axis, and have a certain distance from the imaginary axis.If the distribution pattern of the closed-loop poles of the system agrees with the existing conditions of the dominant poles, the concept of closed-loop dominant poles is often used in engineering to simplify the high-order system to a second-order underdamped system in order to approximately predict the dynamic performance of the high-order system without bringing significant errors to the analysis results. Based on the dynamic response and the simplified analysis method of the high-order system, the maximum value of the real part of the characteristic root max(Re) is proposed as the evaluation index of the angular free vibration stability of DGS. If max(Re) ≥ 0, the floating ring of DGS will not reach a stable state after being instantaneously excited in any angular direction. On the contrary, the transient component of the angular motion time response of the floating ring of DGS will eventually attenuate to zero, and the system is stable. The poles corresponding to max(Re) are generally closed-loop dominant poles, which can indirectly reflect the system’s dynamic performance. The smaller max(Re) is, that is, the farther the dominant poles are from the imaginary axis, the better the angular free vibration stability of DGS is. According to the four stability criteria given in equations (27a)–(27d), the values of the five influencing factors corresponding to the critical instability state of the DGS system can be determined quantitatively, and the intersection of the value ranges of a certain influencing factor required by the four stability criteria is taken as the final value of the influencing factor, as shown in Table 1. It can be seen that only the value of rotational inertia is determined solely by equation (27d). On the contrary, the values of the other four influencing factors are jointly determined by multiple stability criteria. At the same time, under the given parameters, the stability criterion (27d) is the most demanding, and the final value of each characteristic parameter is determined by this stability criterion.Table 1
Characteristic parameter values of DGS system under critical instability state according to Hurwitz stability criterion.
Influencing factorsUnitEquation (27a)Equation (27b)Equation (27c)Equation (27d)Intersection of valuescxx + csxN·m·s·rad−1>0Permanent establishment>22.2474<−6.9260 or >29.1703>29.1703kxx + ksxN·m·rad−1No influence>−5.5563 × 106>2.3881 × 104>2.5190 × 104>2.5190 × 104cyxN·m·s·rad−1No influencePermanent establishment<249.1186<247.2842<247.2842kyxN·m·rad−1No influenceNo influence>−8.8036 × 105>−6.7149 × 105 and <2.8280 × 106>−6.7149 × 105 and <2.8280 × 106Ikg·m2No influencePermanent establishmentNo influence<0.3021<0.3021For the system characteristic equation of the two DOF angular free vibration of DGS shown in equation (21), its four characteristic roots are symmetrical about the real axis on the complex plane. In order to avoid data repetition, the root locus diagram only shows the part above the real axis in the later analysis. In order to study the influence of each influencing factor on the characteristic root, the moving trajectories of two system characteristic roots with nonnegative imaginary parts in three-dimensional space with the change of a certain parameter are made first, and the distance of each characteristic root relative to the imaginary axis is compared to determine the dominant pole. Then the root locus of the dominant pole is enlarged and projected to the complex plane of the characteristic root and the “real part-influencing factor” plane, respectively. The influence of parameters on dynamic performance can be further analyzed qualitatively.Figure4 shows the moving trajectories of the characteristic roots of the system characteristic equation (21) in the complex plane with various influencing factors. Generally, in the study of the angular vibration stability of DGS, the angular main damping cxx of the gas film is taken as the optimization objective, and it is considered that the larger cxx is, the better the stability of the system is. As seen from the line “Re(λ) − cxx + csx” in Figure 4(a), the absolute value of the real part of a pair of conjugate complex roots increases monotonously with the increase of cxx + csx, and the real part is always less than zero. The corresponding component remains a stable decreasing amplitude oscillation, and the attenuation rate of the oscillation is accelerated. However, the real part value of the other pair of characteristic roots is more considerable, which determines the dynamic performance of the angular vibration of DGS. The variation of the real part value with cxx + csx is not monotonous and remains near the cxx + csx axis all the time. When cxx + csx is less than its critical value, max(Re) is positive, and the angular vibration of DGS is unstable. It can be seen from Table 1 that the critical value of cxx + csx will be restricted by the rotational inertia criterion (27d), which is a positive number with a small absolute value. It is not true that the angular vibration of DGS will be stable when the angular main damping cxx + csx of the system is greater than zero, and cxx + csx > 0 is a necessary and insufficient condition. When cxx + csx is further increased beyond the critical value, the distance between the characteristic root and the imaginary axis on the complex plane increases at first and then decreases, while the damping angle β decreases continuously. When the damping angle β is about 45°, the system has both good damping characteristics and a good response speed. When cxx + csx approaches infinity, the characteristic root approaches the origin infinitely in the second quadrant. Although the system’s time response is stable, the transient performance at this time is abysmal.Figure 4
Root locus diagram under the influence of a single factor: (a)cxx + csx; (b) kxx + ksx; (c) cyx; (d) kyx and (e) I.
(a)(b)(c)(d)(e)In the “Re(λ) − kxx + ksx” diagram (Figure 4(b)), it can be seen that when kxx + ksx tends to positive infinity, both variation curves have horizontal asymptotes below the transverse axis. With the increase of kxx + ksx, max(Re) decreases monotonously. There is a critical value of kxx + ksx in the interval of kxx + ksx > 0, and the system is stable when kxx + ksx is greater than its critical value. Similarly, the critical value is mainly determined by the rotational inertia criterion (27d). Under the calculation parameters, when kxx + ksx is large enough (>4 × 106 N·m/rad), if kxx + ksx continues to increase, the characteristic root’s real part value only slightly decreases, but the imaginary part value increases rapidly. On the complex plane, the characteristic root moves to the upper left in the second quadrant, and there is a vertical asymptote, indicating that the adjustment time of the system does not change much, but the amplitude and frequency of the time response will increase. The floating ring of DGS will be in a high-frequency and large-amplitude swing state in a short time after being instantaneously excited in any angular direction. If the installation structure of DGS limits the maximum deflection angle of the floating ring, the excessive kxx + ksx will make the angular vibration amplitude of DGS exceed the acceptable range and cause failure.With the increase ofcyx, max(Re) first decreases, then increases and exceeds zero as shown in Figure 4(c), indicating that the angular free vibration of DGS will be unstable, and when it reaches the maximum value, max(Re) decreases slightly. As shown in Table 1, the system is stable when cyx is less than its critical value under the calculated parameters. Still, it is not true that the smaller the cyx, the better the system’s dynamic performance. The root locus on the complex plane has an intersection point with the negative half axis of the real axis, reflecting that the dominant poles are not a pair of conjugate complex roots at this time, but multiple second-order roots. The corresponding response function is monotonous exponential attenuation, and its response speed is not as fast as that of the second-order underdamped system. If cyx is further reduced, the damping angle β will increase while the absolute value of the real part will decrease, and the rapidity and stationarity of the system response will be weakened.As can be seen from Figure4(d), with the increase of kyx, max(Re) first decreases and then increases. When the absolute value of kyx is large, the angular vibration response of DGS will be unstable and divergent. The moving trajectory of the dominant pole on the complex plane passes through the imaginary axis twice. In conjunction with Table 1, it is known that under the calculation parameters, the angular free vibration of DGS is stable under certain conditions, and kyx has both upper and lower critical values at the same time. With increasing kyx, the characteristic root moves to the lower left of the complex plane through the imaginary axis and has an intersection point with the negative half-axis of the real axis. With the further increase of kyx, the absolute value of the real part and the damping angle β also increase. Under the premise of ensuring the stability of the system, appropriately increasing the damping angle β reduces the damping ratio ζ, which can make the system have a faster response speed. At this time, kyx is greater than zero but less than its upper limit, that is, there is a preferred value of kyx in a small range of kyx > 0 so that the system has the best transient performance.For the equation of motion (5) of the angular free vibration of the floating ring, the main stiffness and the main damping are separated from the cross stiffness and the cross damping, and the following equation can be obtained:(28)Ixφ¨xt+cxx+csxφ˙xt+kxx+ksxφxt=−cxyφ˙yt−kxyφyt,Iyφ¨yt+cyy+csyφ˙yt+kyy+ksyφyt=−cyxφ˙xt−kyxφxt.It can be seen from the form of equation (28) that the angular deflection φx around the x-axis and the angular deflection φy around the y-axis, which are related by the angular cross stiffness and the angular cross damping, are the inputs of the system to each other, and coordinate and influence each other. Therefore, the approach of the absolute value of the cross coefficient to zero is beneficial to the suppression of angular free vibration, which is consistent with the laws shown in Figures 4(c) and 4(d).With the increase ofI, the real part values of the two pairs of characteristic roots of the system all increase rapidly to the maximum at first and then decrease slightly as shown in Figure 4(e). There is a unique critical value for I, and the system is stable only when I is less than its critical value. In the complex plane, with the increase of I, the characteristic root first moves towards the imaginary axis in the second quadrant. In this process, the distance relative to the real axis increases first and then decreases. But the damping angle β keeps the trend of increasing, indicating that the damping degree and the response speed of the system decrease simultaneously, that is, the dynamic performance index deteriorates, and the free vibration stability decreases. After passing through the imaginary axis, the characteristic root first moves away from the imaginary axis, then returns, and gradually moves to the coordinate origin. Under the calculated parameters, the smaller the rotational inertia, the better the angular free vibration stability of the system.Based on the above analysis, the effects of rotational inertia, stiffness and damping parameters on the angular free vibration stability of DGS are different. When analyzing and designing the geometric parameters of DGS, it is necessary to reasonably match the relevant parameters in combination with the system characteristic equation and the stability criterion to achieve the optimal transient performance. It is worth noting that the constraint condition of rotational inertia (27d) is the most demanding of all stability conditions, which determines the dynamic characteristic coefficients of the gas film in the critical instability state of the system. Therefore, the following discussion will focus on the interaction between rotational inertia and other parameters on the angular vibration stability of DGS.


### 3.2. Interaction Influence of Rotational Inertia and Stiffness or Damping Parameters

Figure5 shows the interaction between rotational inertia and stiffness or damping parameters. From the perspective of parameter matching, there is an obvious interaction between rotational inertia and stiffness or damping parameters. Under different parameter matching, the influence law of each parameter on the stability of angular free vibration may change. Specifically, as shown in Figure (5a), when the absolute value of kyx is large, max(Re) first increases rapidly with the increase of I and then changes little after reaching the maximum value. When the absolute value is the same, the change rate of max(Re) in the case of positive cross stiffness is significantly greater than that in the case of negative cross stiffness. This phenomenon results in the existence of a rotational inertia threshold, which changes the variation law of max(Re) with kyx. When I is tiny, max(Re) decreases monotonously with the increase of kyx, and the greater the kyx, the better the angular free vibration stability of the system. After exceeding the threshold of rotational inertia, max(Re) first decreases, then increases, and exceeds zero with the increase of kyx. At this time, the angular vibration response is unstable and divergent. When I is greater than a certain threshold, kyx will have both upper and lower critical values at the same time, and when its absolute value is small, the angular stability of the system is better.Figure 5
Interaction between rotational inertia and stiffness or damping parameters: (a)I − kyx; (b) I − cyx; (c) I − (kxx + ksx); and (d) I − (cxx + csx).
(a)(b)(c)(d)As seen from Figure5(b), the value of I does not change the variation law of max(Re) with cyx, but only the critical cross damping value and the position of the local minimum point. With the increase of I, the critical cross damping value and the local minimum point will all shift to the left. The local minimum value of max(Re) gradually increases and approaches zero. When I is large enough, the critical cross damping will be negative, and when cyx is less than its critical value, the system is stable. Under different cyx, the variation law of max(Re) with I is also different. When the absolute value of cyx is large, max(Re) increases monotonously with the increase of I. When the absolute value of cyx is close to zero, max(Re) decreases with the increase of I at first, then increases gradually and will exceed zero. That is, the stability of the angular free vibration becomes better at first and then becomes unstable gradually. The system has better angular free vibration stability when the absolute values of I and cyx are small.With the increase ofkxx + ksx, max(Re) decreases monotonously as shown in Figure 5(c), and the stability of angular free vibration becomes better. Therefore, when the optimization goal is to improve the angular free vibration stability of the system, the angular main stiffness of the gas film should be moderately enhanced by designing the end face structure. With the increase of I, the decreasing rate of the “max(Re) − kxx + ksx” line decreases. The zero-point where the line crosses the transverse axis shifts to the right, that is, the critical main stiffness increases with the increase of I. When kxx + ksx is greater than its critical value, the system is stable, so the increase of I will reduce the stability range of the system. The effect of I on max(Re) is significantly different under different main stiffness values. When kxx + ksx is large enough, with the increase of I, max(Re) first decreases rapidly, then increases after reaching the local minimum value, and the rate of increase decreases, and finally max(Re) gradually exceeds zero. When I is small, and kxx + ksx is large enough, the angular free vibration stability of the system is better.It can be seen in Figure5(d) that with the increase of cxx + csx, max(Re) first decreases rapidly, then increases after reaching the local minimum value, and the rate of increase decreases. When cxx + csx tends to infinity, max(Re) will approach infinitely without exceeding zero. This trend means that for the stability of angular free vibration, cxx + csx is not the larger the better, but there is an optimal range of values. On the premise of ensuring the system’s stability, appropriately reducing the main damping value is beneficial to improve the response speed of the angular free vibration of the system and reduce the adjustment time. With the increase of I, the critical main damping value increases, and the system is stable when the main damping is greater than its critical value. When cxx + csx is large enough, the change of I has little effect on max(Re). Under the condition that cxx + csx is close to zero (the angular vibration is unstable at this time), the resistance of the sealing ring to the change of its motion state is weak when I is small, and it is more prone to deflection due to angular excitation. Therefore, there is a larger max(Re) at this time, and the amplitude of the transient response oscillation of the angular motion of the floating ring increases faster. In severe cases, maintenance personnel on the front line of work may not have time to respond at all, and direct contact friction will occur on sealing faces. That is, when the angular motion is unstable, the gas seal system with the more negligible rotational inertia of the floating ring has a higher risk of sudden failure.


### 3.3. Application Case and Validity Proof

In order to verify the validity of the method for judging the angular free vibration stability of DGS in this paper, the seal stability under typical working conditions is discussed in this section. Firstly, based on the frequency perturbation method mentioned in Section2.2.2, the angular dynamic characteristic coefficients of gas film under given geometric and operating parameters are solved, and the algebraic stability criterion derived in this paper is used to calculate the critical rotational inertia of DGS, so as to evaluate the stability range of angular vibration of DGS. Furthermore, the angular motion response of the flexibly mounted stator under instantaneous excitation is accurately and completely obtained by using the nonlinear direct numerical simulation method (DNSM). Figure 6 shows the specific flow of using DNSM to trace the dynamic sealing behavior after disturbance. DNSM involves solving the gas film lubrication equation and the dynamic equations of the seal ring simultaneously at each moment to obtain the complete time histories of the gas film pressure distribution and the motion state of the seal ring. Specifically, the successive over relaxation method is used to solve the transient gas film pressure distribution, and the sealing performance parameters are calculated. The overturning moments on the seal face are substituted into the dynamic equations, and the angular response displacement and velocity of the stator at the next moment are calculated recursively by using the fourth-order classical Runge-Kutta scheme. The film thickness distribution between the seal faces is updated according to the relative position of the seal rings (as shown in Figure 1(c)). The coupling process of the above fluid governing equation, dynamic equations and time-varying function of film thickness distribution is repeated at each time step until all time nodes are calculated.Figure 6
Calculation process of the direct numerical simulation method.Before comparing the prediction results of the angular vibration stability range of DGS between the algebraic stability criterion proposed in this paper and the traditional nonlinear numerical solution method, the calculated values of the film dynamic coefficients in this paper are compared with the numerical results in reference [32]. The results are shown in Figure 7(a). It can be seen that with the increase of spiral angle α, the calculated results of dynamic characteristic coefficients of gas film in this paper are close to those in the literature, and the trend is consistent. In addition, the steady-state angular response of the floating ring under the continuous excitation of the rotor runout is calculated in reference [30]. By changing the initial values of the kinematic variables of the floating ring and the excitation form of the nonfloating ring in the film thickness equation, DNSM can be used to solve both the free vibration response and the forced vibration response of DGS. Figure 7(b) shows the comparison between the angular response displacement obtained in this paper and the calculated result in reference [30]. It can be seen that the angular motion of the floating ring calculated by the nonlinear simulation method in this paper is in good agreement with the result in the literature, which verifies the correctness of the relevant numerical calculation methods in this paper. In the following examples, unless otherwise specified, the default values of geometric and operating parameters are shown in Table 2.Figure 7
Verification of angular dynamic characteristics of DGS: (a) film dynamic coefficients and (b) angular response displacement of the floating ring.
(a)(b)Table 2
The default values of geometric and operating parameters.
ParametersValueSeal ring inner radiusri, mm30Seal ring outer radiusro, mm42Spiral groove inner dam radiusrg, mm33.6Groove-to-land ratioκ1Spiral angleα, °20Number of groovesNg12Equilibrium ratioB0.8Gas viscosityμ, Pa·s1.79 × 10−5Ambient pressurepa, MPa0.101Sealed gas pressurepo, MPa0.303Rotational angular velocity Ω, rad·s−1376Equilibrium film thicknesshb, μm8Spring stiffnessks, N·m−11 × 105Auxiliary seal dampingcs, N·s·m−1200Transverse rotational inertia of the flexibly mounted statorI, kg·m21.33 × 10−4Figure8(a) shows the variation of the gas film dynamic characteristics with the spiral groove depth. It can be seen that with the increase of groove depth, the absolute values of cross stiffness kyx and cross damping cyx decrease at first and then tend to smooth. The influence laws of groove depth on main stiffness kxx and main damping cxx are opposite. With the increase of groove depth, the main stiffness increases nearly linearly, while the main damping decreases gradually. However, in the research range of groove depth, the main damping always remains positive. Therefore, according to the previous viewpoint of taking damping as the evaluation index, it will be mistaken that the angular vibration of DGS is stable, and the increase of groove depth will weaken the angular vibration stability.Figure 8
Variation of dynamic characteristics of DGS with groove depth: (a) film dynamic coefficients; (b) critical rotational inertia; and (c) root locus.
(a)(b)(c)According to the algebraic stability criterion (27) provided in this paper, the value of rotational inertia corresponding to the critical instability state of the DGS system is calculated. Figure8(b) shows the variation of the difference between the critical rotational inertia Icr and the transverse rotational inertia I of the flexibly mounted stator with the groove depth. When Icr − I < 0, the angular vibration of DGS is easy to lose stability and produce self-excited vibration. As can be seen from Figure 8(b), Icr − I shows an increasing trend with the increase of groove depth, but the rate of increase decreases. Since the main stiffness kxx is negative when the groove depth is small, there is an unstable region of angular vibration (the pink area in Figure 8(b)). According to the method for judging the angular free vibration stability of DGS in this paper, when hg < 4.3 μm, the angular vibration of DGS is unstable. At this time, DGS is prone to failure accidents due to the instability of gas film.The influence of spiral groove depth on the angular dynamic performance of DGS can be qualitatively analyzed by observing the moving trajectory of the dominant pole in the complex plane with the change of groove depth in Figure8(c). It can be found that with the increase of groove depth, the dominant pole will pass through the imaginary axis in the complex plane, and then continue to move to the upper left in the second quadrant. During the movement of the dominant pole, the damping angle β increases slightly. This means that when the groove depth is large, there is a large saltation peak after the floating ring is disturbed. However, since the distances from the dominant pole to both the imaginary axis and the real axis increase, the system has better damping degree and response speed. The angular vibration of the floating ring will attenuate rapidly, and the amplitude of the subsequent oscillation is smaller. The gas seal system can return to steady state more quickly.Based on the direct numerical simulation method, the complex nonlinear flow field between the seal faces and the two DOF angular motion of the seal rings are coupled to obtain the complete time history of angular vibration of the floating ring, and the angular vibration stability of DGS under different groove depths is analyzed. An initial swing angular velocity of 1.88 × 104 μrad/s is applied to the stator of DGS in the y direction. Thus, the stator is excited and begins to vibrate in both x and y directions. Figure 9 shows the time domain response of the stator in the y direction to the initial swing angular velocity under different groove depths. When hg = 3 μm, the stator response under the initial velocity disturbance is a divergent response process (as shown in Figure 9(a)). The angular response displacement of the stator cannot reach a steady state and gradually diverges with the passage of time. Within just three rotation cycles, the angular displacement φy has reached 80 μrad. As can be seen from Figure 9(b), when hg > 5 μm, the gas seal system is in a stable state, and the angular response of the stator is an attenuated oscillation. By comparing the vibration peaks and the time to restore steady state of each response curve, it can also be found that with the increase of groove depth, the first peak value of the angular response displacement of the stator increases, but both the intensity and response time of the subsequent oscillation obviously decrease. In other words, the angular vibration stability of the system is improved. The comparative analysis with Figure 8 reveals that the stability range predicted based on the algebraic stability criterion and the dynamic performance obtained by the root locus method are basically consistent with the variation law of angular vibration characteristics of DGS with groove depth shown in Figure 9, which proves the feasibility, effectiveness and accuracy of using the methods of control theory to study the angular free vibration stability of DGS.Figure 9
Angular response displacement of the floating ring solved by DNSM under different groove depth: (a)hg = 3 μm and (b)hg > 5 μm.
(a)(b)In summary, combined with the direct numerical simulation results of the stator transient response, the validity of the method for judging the angular free vibration stability of DGS is verified. The results show that the rules and conclusions obtained by the two research methods are consistent, and the method presented in this paper can effectively predict the angular vibration stability of the DGS system and its dynamic performance in the transition stage. There is no doubt that the most accurate research method is to solve the complete angular vibration behavior of the floating ring. However, DNSM has a large amount of calculation and storage, and each calculation is only for the specified parameters, which is not suitable for seal design and structure optimization. In this paper, the linearization method is used to represent the features of gas film between the seal faces, and the angular dynamic performance of DGS are directly analyzed according to the system characteristic equation, without real-time solution of the flow field. Therefore, the calculation efficiency is much higher, and it is suitable for large-scale parameter discussion. In addition, the definition and calculation formula of critical rotational inertiaIcr are given in this paper. It is generally known that the larger the film angular main damping cxx is, the slower the response speed of the system is. Thus, Icr can more fully reflect the angular free vibration characteristics of DGS (including damping degree and response speed). As a comprehensive performance indicator, Icr can be used as an optimization objective of the angular dynamic performance of DGS in future research. On the premise that the main damping is greater than zero, the larger the Icr − I is, the more excellent the system transient performance in the transition stage is.


## 3.1. Influence of Parameters on Angular Free Vibration Stability of DGS

According to the automatic control theory, for a stable high-order system, the distance between the closed-loop pole and the imaginary axis determines the attenuation speed of the corresponding transient component, and the farther the distance, the faster the attenuation. The closed-loop zero-point (the point satisfying the molecular of the transfer function to be zero) affects the relative intensity of each transient component. If there is a zero-point near a certain pole, the intensity of its corresponding transient component will become smaller, and the transient response component of a pair of zero-point and pole that are close to each other can be ignored. Suppose among all the closed-loop poles, there is no closed-loop zero-point around the pole closest to the imaginary axis, and other closed-loop poles are far away from the imaginary axis. In that case, the corresponding transient component of the closed-loop pole closest to the imaginary axis in the unit step response not only has the largest initial value whent = 0, but also attenuates the slowest in all responses. Therefore, it plays an important role in the time response of the system. Such a closed-loop pole is called the closed-loop dominant pole, and other poles are called nondominant poles. The dominant poles are generally a pair of complex conjugate poles, which are located on the left side of the imaginary axis, in the sector with an angle of 60° to the real axis, and have a certain distance from the imaginary axis.If the distribution pattern of the closed-loop poles of the system agrees with the existing conditions of the dominant poles, the concept of closed-loop dominant poles is often used in engineering to simplify the high-order system to a second-order underdamped system in order to approximately predict the dynamic performance of the high-order system without bringing significant errors to the analysis results. Based on the dynamic response and the simplified analysis method of the high-order system, the maximum value of the real part of the characteristic root max(Re) is proposed as the evaluation index of the angular free vibration stability of DGS. If max(Re) ≥ 0, the floating ring of DGS will not reach a stable state after being instantaneously excited in any angular direction. On the contrary, the transient component of the angular motion time response of the floating ring of DGS will eventually attenuate to zero, and the system is stable. The poles corresponding to max(Re) are generally closed-loop dominant poles, which can indirectly reflect the system’s dynamic performance. The smaller max(Re) is, that is, the farther the dominant poles are from the imaginary axis, the better the angular free vibration stability of DGS is. According to the four stability criteria given in equations (27a)–(27d), the values of the five influencing factors corresponding to the critical instability state of the DGS system can be determined quantitatively, and the intersection of the value ranges of a certain influencing factor required by the four stability criteria is taken as the final value of the influencing factor, as shown in Table 1. It can be seen that only the value of rotational inertia is determined solely by equation (27d). On the contrary, the values of the other four influencing factors are jointly determined by multiple stability criteria. At the same time, under the given parameters, the stability criterion (27d) is the most demanding, and the final value of each characteristic parameter is determined by this stability criterion.Table 1
Characteristic parameter values of DGS system under critical instability state according to Hurwitz stability criterion.
Influencing factorsUnitEquation (27a)Equation (27b)Equation (27c)Equation (27d)Intersection of valuescxx + csxN·m·s·rad−1>0Permanent establishment>22.2474<−6.9260 or >29.1703>29.1703kxx + ksxN·m·rad−1No influence>−5.5563 × 106>2.3881 × 104>2.5190 × 104>2.5190 × 104cyxN·m·s·rad−1No influencePermanent establishment<249.1186<247.2842<247.2842kyxN·m·rad−1No influenceNo influence>−8.8036 × 105>−6.7149 × 105 and <2.8280 × 106>−6.7149 × 105 and <2.8280 × 106Ikg·m2No influencePermanent establishmentNo influence<0.3021<0.3021For the system characteristic equation of the two DOF angular free vibration of DGS shown in equation (21), its four characteristic roots are symmetrical about the real axis on the complex plane. In order to avoid data repetition, the root locus diagram only shows the part above the real axis in the later analysis. In order to study the influence of each influencing factor on the characteristic root, the moving trajectories of two system characteristic roots with nonnegative imaginary parts in three-dimensional space with the change of a certain parameter are made first, and the distance of each characteristic root relative to the imaginary axis is compared to determine the dominant pole. Then the root locus of the dominant pole is enlarged and projected to the complex plane of the characteristic root and the “real part-influencing factor” plane, respectively. The influence of parameters on dynamic performance can be further analyzed qualitatively.Figure4 shows the moving trajectories of the characteristic roots of the system characteristic equation (21) in the complex plane with various influencing factors. Generally, in the study of the angular vibration stability of DGS, the angular main damping cxx of the gas film is taken as the optimization objective, and it is considered that the larger cxx is, the better the stability of the system is. As seen from the line “Re(λ) − cxx + csx” in Figure 4(a), the absolute value of the real part of a pair of conjugate complex roots increases monotonously with the increase of cxx + csx, and the real part is always less than zero. The corresponding component remains a stable decreasing amplitude oscillation, and the attenuation rate of the oscillation is accelerated. However, the real part value of the other pair of characteristic roots is more considerable, which determines the dynamic performance of the angular vibration of DGS. The variation of the real part value with cxx + csx is not monotonous and remains near the cxx + csx axis all the time. When cxx + csx is less than its critical value, max(Re) is positive, and the angular vibration of DGS is unstable. It can be seen from Table 1 that the critical value of cxx + csx will be restricted by the rotational inertia criterion (27d), which is a positive number with a small absolute value. It is not true that the angular vibration of DGS will be stable when the angular main damping cxx + csx of the system is greater than zero, and cxx + csx > 0 is a necessary and insufficient condition. When cxx + csx is further increased beyond the critical value, the distance between the characteristic root and the imaginary axis on the complex plane increases at first and then decreases, while the damping angle β decreases continuously. When the damping angle β is about 45°, the system has both good damping characteristics and a good response speed. When cxx + csx approaches infinity, the characteristic root approaches the origin infinitely in the second quadrant. Although the system’s time response is stable, the transient performance at this time is abysmal.Figure 4
Root locus diagram under the influence of a single factor: (a)cxx + csx; (b) kxx + ksx; (c) cyx; (d) kyx and (e) I.
(a)(b)(c)(d)(e)In the “Re(λ) − kxx + ksx” diagram (Figure 4(b)), it can be seen that when kxx + ksx tends to positive infinity, both variation curves have horizontal asymptotes below the transverse axis. With the increase of kxx + ksx, max(Re) decreases monotonously. There is a critical value of kxx + ksx in the interval of kxx + ksx > 0, and the system is stable when kxx + ksx is greater than its critical value. Similarly, the critical value is mainly determined by the rotational inertia criterion (27d). Under the calculation parameters, when kxx + ksx is large enough (>4 × 106 N·m/rad), if kxx + ksx continues to increase, the characteristic root’s real part value only slightly decreases, but the imaginary part value increases rapidly. On the complex plane, the characteristic root moves to the upper left in the second quadrant, and there is a vertical asymptote, indicating that the adjustment time of the system does not change much, but the amplitude and frequency of the time response will increase. The floating ring of DGS will be in a high-frequency and large-amplitude swing state in a short time after being instantaneously excited in any angular direction. If the installation structure of DGS limits the maximum deflection angle of the floating ring, the excessive kxx + ksx will make the angular vibration amplitude of DGS exceed the acceptable range and cause failure.With the increase ofcyx, max(Re) first decreases, then increases and exceeds zero as shown in Figure 4(c), indicating that the angular free vibration of DGS will be unstable, and when it reaches the maximum value, max(Re) decreases slightly. As shown in Table 1, the system is stable when cyx is less than its critical value under the calculated parameters. Still, it is not true that the smaller the cyx, the better the system’s dynamic performance. The root locus on the complex plane has an intersection point with the negative half axis of the real axis, reflecting that the dominant poles are not a pair of conjugate complex roots at this time, but multiple second-order roots. The corresponding response function is monotonous exponential attenuation, and its response speed is not as fast as that of the second-order underdamped system. If cyx is further reduced, the damping angle β will increase while the absolute value of the real part will decrease, and the rapidity and stationarity of the system response will be weakened.As can be seen from Figure4(d), with the increase of kyx, max(Re) first decreases and then increases. When the absolute value of kyx is large, the angular vibration response of DGS will be unstable and divergent. The moving trajectory of the dominant pole on the complex plane passes through the imaginary axis twice. In conjunction with Table 1, it is known that under the calculation parameters, the angular free vibration of DGS is stable under certain conditions, and kyx has both upper and lower critical values at the same time. With increasing kyx, the characteristic root moves to the lower left of the complex plane through the imaginary axis and has an intersection point with the negative half-axis of the real axis. With the further increase of kyx, the absolute value of the real part and the damping angle β also increase. Under the premise of ensuring the stability of the system, appropriately increasing the damping angle β reduces the damping ratio ζ, which can make the system have a faster response speed. At this time, kyx is greater than zero but less than its upper limit, that is, there is a preferred value of kyx in a small range of kyx > 0 so that the system has the best transient performance.For the equation of motion (5) of the angular free vibration of the floating ring, the main stiffness and the main damping are separated from the cross stiffness and the cross damping, and the following equation can be obtained:(28)Ixφ¨xt+cxx+csxφ˙xt+kxx+ksxφxt=−cxyφ˙yt−kxyφyt,Iyφ¨yt+cyy+csyφ˙yt+kyy+ksyφyt=−cyxφ˙xt−kyxφxt.It can be seen from the form of equation (28) that the angular deflection φx around the x-axis and the angular deflection φy around the y-axis, which are related by the angular cross stiffness and the angular cross damping, are the inputs of the system to each other, and coordinate and influence each other. Therefore, the approach of the absolute value of the cross coefficient to zero is beneficial to the suppression of angular free vibration, which is consistent with the laws shown in Figures 4(c) and 4(d).With the increase ofI, the real part values of the two pairs of characteristic roots of the system all increase rapidly to the maximum at first and then decrease slightly as shown in Figure 4(e). There is a unique critical value for I, and the system is stable only when I is less than its critical value. In the complex plane, with the increase of I, the characteristic root first moves towards the imaginary axis in the second quadrant. In this process, the distance relative to the real axis increases first and then decreases. But the damping angle β keeps the trend of increasing, indicating that the damping degree and the response speed of the system decrease simultaneously, that is, the dynamic performance index deteriorates, and the free vibration stability decreases. After passing through the imaginary axis, the characteristic root first moves away from the imaginary axis, then returns, and gradually moves to the coordinate origin. Under the calculated parameters, the smaller the rotational inertia, the better the angular free vibration stability of the system.Based on the above analysis, the effects of rotational inertia, stiffness and damping parameters on the angular free vibration stability of DGS are different. When analyzing and designing the geometric parameters of DGS, it is necessary to reasonably match the relevant parameters in combination with the system characteristic equation and the stability criterion to achieve the optimal transient performance. It is worth noting that the constraint condition of rotational inertia (27d) is the most demanding of all stability conditions, which determines the dynamic characteristic coefficients of the gas film in the critical instability state of the system. Therefore, the following discussion will focus on the interaction between rotational inertia and other parameters on the angular vibration stability of DGS.


## 3.2. Interaction Influence of Rotational Inertia and Stiffness or Damping Parameters

Figure5 shows the interaction between rotational inertia and stiffness or damping parameters. From the perspective of parameter matching, there is an obvious interaction between rotational inertia and stiffness or damping parameters. Under different parameter matching, the influence law of each parameter on the stability of angular free vibration may change. Specifically, as shown in Figure (5a), when the absolute value of kyx is large, max(Re) first increases rapidly with the increase of I and then changes little after reaching the maximum value. When the absolute value is the same, the change rate of max(Re) in the case of positive cross stiffness is significantly greater than that in the case of negative cross stiffness. This phenomenon results in the existence of a rotational inertia threshold, which changes the variation law of max(Re) with kyx. When I is tiny, max(Re) decreases monotonously with the increase of kyx, and the greater the kyx, the better the angular free vibration stability of the system. After exceeding the threshold of rotational inertia, max(Re) first decreases, then increases, and exceeds zero with the increase of kyx. At this time, the angular vibration response is unstable and divergent. When I is greater than a certain threshold, kyx will have both upper and lower critical values at the same time, and when its absolute value is small, the angular stability of the system is better.Figure 5
Interaction between rotational inertia and stiffness or damping parameters: (a)I − kyx; (b) I − cyx; (c) I − (kxx + ksx); and (d) I − (cxx + csx).
(a)(b)(c)(d)As seen from Figure5(b), the value of I does not change the variation law of max(Re) with cyx, but only the critical cross damping value and the position of the local minimum point. With the increase of I, the critical cross damping value and the local minimum point will all shift to the left. The local minimum value of max(Re) gradually increases and approaches zero. When I is large enough, the critical cross damping will be negative, and when cyx is less than its critical value, the system is stable. Under different cyx, the variation law of max(Re) with I is also different. When the absolute value of cyx is large, max(Re) increases monotonously with the increase of I. When the absolute value of cyx is close to zero, max(Re) decreases with the increase of I at first, then increases gradually and will exceed zero. That is, the stability of the angular free vibration becomes better at first and then becomes unstable gradually. The system has better angular free vibration stability when the absolute values of I and cyx are small.With the increase ofkxx + ksx, max(Re) decreases monotonously as shown in Figure 5(c), and the stability of angular free vibration becomes better. Therefore, when the optimization goal is to improve the angular free vibration stability of the system, the angular main stiffness of the gas film should be moderately enhanced by designing the end face structure. With the increase of I, the decreasing rate of the “max(Re) − kxx + ksx” line decreases. The zero-point where the line crosses the transverse axis shifts to the right, that is, the critical main stiffness increases with the increase of I. When kxx + ksx is greater than its critical value, the system is stable, so the increase of I will reduce the stability range of the system. The effect of I on max(Re) is significantly different under different main stiffness values. When kxx + ksx is large enough, with the increase of I, max(Re) first decreases rapidly, then increases after reaching the local minimum value, and the rate of increase decreases, and finally max(Re) gradually exceeds zero. When I is small, and kxx + ksx is large enough, the angular free vibration stability of the system is better.It can be seen in Figure5(d) that with the increase of cxx + csx, max(Re) first decreases rapidly, then increases after reaching the local minimum value, and the rate of increase decreases. When cxx + csx tends to infinity, max(Re) will approach infinitely without exceeding zero. This trend means that for the stability of angular free vibration, cxx + csx is not the larger the better, but there is an optimal range of values. On the premise of ensuring the system’s stability, appropriately reducing the main damping value is beneficial to improve the response speed of the angular free vibration of the system and reduce the adjustment time. With the increase of I, the critical main damping value increases, and the system is stable when the main damping is greater than its critical value. When cxx + csx is large enough, the change of I has little effect on max(Re). Under the condition that cxx + csx is close to zero (the angular vibration is unstable at this time), the resistance of the sealing ring to the change of its motion state is weak when I is small, and it is more prone to deflection due to angular excitation. Therefore, there is a larger max(Re) at this time, and the amplitude of the transient response oscillation of the angular motion of the floating ring increases faster. In severe cases, maintenance personnel on the front line of work may not have time to respond at all, and direct contact friction will occur on sealing faces. That is, when the angular motion is unstable, the gas seal system with the more negligible rotational inertia of the floating ring has a higher risk of sudden failure.


## 3.3. Application Case and Validity Proof

In order to verify the validity of the method for judging the angular free vibration stability of DGS in this paper, the seal stability under typical working conditions is discussed in this section. Firstly, based on the frequency perturbation method mentioned in Section2.2.2, the angular dynamic characteristic coefficients of gas film under given geometric and operating parameters are solved, and the algebraic stability criterion derived in this paper is used to calculate the critical rotational inertia of DGS, so as to evaluate the stability range of angular vibration of DGS. Furthermore, the angular motion response of the flexibly mounted stator under instantaneous excitation is accurately and completely obtained by using the nonlinear direct numerical simulation method (DNSM). Figure 6 shows the specific flow of using DNSM to trace the dynamic sealing behavior after disturbance. DNSM involves solving the gas film lubrication equation and the dynamic equations of the seal ring simultaneously at each moment to obtain the complete time histories of the gas film pressure distribution and the motion state of the seal ring. Specifically, the successive over relaxation method is used to solve the transient gas film pressure distribution, and the sealing performance parameters are calculated. The overturning moments on the seal face are substituted into the dynamic equations, and the angular response displacement and velocity of the stator at the next moment are calculated recursively by using the fourth-order classical Runge-Kutta scheme. The film thickness distribution between the seal faces is updated according to the relative position of the seal rings (as shown in Figure 1(c)). The coupling process of the above fluid governing equation, dynamic equations and time-varying function of film thickness distribution is repeated at each time step until all time nodes are calculated.Figure 6
Calculation process of the direct numerical simulation method.Before comparing the prediction results of the angular vibration stability range of DGS between the algebraic stability criterion proposed in this paper and the traditional nonlinear numerical solution method, the calculated values of the film dynamic coefficients in this paper are compared with the numerical results in reference [32]. The results are shown in Figure 7(a). It can be seen that with the increase of spiral angle α, the calculated results of dynamic characteristic coefficients of gas film in this paper are close to those in the literature, and the trend is consistent. In addition, the steady-state angular response of the floating ring under the continuous excitation of the rotor runout is calculated in reference [30]. By changing the initial values of the kinematic variables of the floating ring and the excitation form of the nonfloating ring in the film thickness equation, DNSM can be used to solve both the free vibration response and the forced vibration response of DGS. Figure 7(b) shows the comparison between the angular response displacement obtained in this paper and the calculated result in reference [30]. It can be seen that the angular motion of the floating ring calculated by the nonlinear simulation method in this paper is in good agreement with the result in the literature, which verifies the correctness of the relevant numerical calculation methods in this paper. In the following examples, unless otherwise specified, the default values of geometric and operating parameters are shown in Table 2.Figure 7
Verification of angular dynamic characteristics of DGS: (a) film dynamic coefficients and (b) angular response displacement of the floating ring.
(a)(b)Table 2
The default values of geometric and operating parameters.
ParametersValueSeal ring inner radiusri, mm30Seal ring outer radiusro, mm42Spiral groove inner dam radiusrg, mm33.6Groove-to-land ratioκ1Spiral angleα, °20Number of groovesNg12Equilibrium ratioB0.8Gas viscosityμ, Pa·s1.79 × 10−5Ambient pressurepa, MPa0.101Sealed gas pressurepo, MPa0.303Rotational angular velocity Ω, rad·s−1376Equilibrium film thicknesshb, μm8Spring stiffnessks, N·m−11 × 105Auxiliary seal dampingcs, N·s·m−1200Transverse rotational inertia of the flexibly mounted statorI, kg·m21.33 × 10−4Figure8(a) shows the variation of the gas film dynamic characteristics with the spiral groove depth. It can be seen that with the increase of groove depth, the absolute values of cross stiffness kyx and cross damping cyx decrease at first and then tend to smooth. The influence laws of groove depth on main stiffness kxx and main damping cxx are opposite. With the increase of groove depth, the main stiffness increases nearly linearly, while the main damping decreases gradually. However, in the research range of groove depth, the main damping always remains positive. Therefore, according to the previous viewpoint of taking damping as the evaluation index, it will be mistaken that the angular vibration of DGS is stable, and the increase of groove depth will weaken the angular vibration stability.Figure 8
Variation of dynamic characteristics of DGS with groove depth: (a) film dynamic coefficients; (b) critical rotational inertia; and (c) root locus.
(a)(b)(c)According to the algebraic stability criterion (27) provided in this paper, the value of rotational inertia corresponding to the critical instability state of the DGS system is calculated. Figure8(b) shows the variation of the difference between the critical rotational inertia Icr and the transverse rotational inertia I of the flexibly mounted stator with the groove depth. When Icr − I < 0, the angular vibration of DGS is easy to lose stability and produce self-excited vibration. As can be seen from Figure 8(b), Icr − I shows an increasing trend with the increase of groove depth, but the rate of increase decreases. Since the main stiffness kxx is negative when the groove depth is small, there is an unstable region of angular vibration (the pink area in Figure 8(b)). According to the method for judging the angular free vibration stability of DGS in this paper, when hg < 4.3 μm, the angular vibration of DGS is unstable. At this time, DGS is prone to failure accidents due to the instability of gas film.The influence of spiral groove depth on the angular dynamic performance of DGS can be qualitatively analyzed by observing the moving trajectory of the dominant pole in the complex plane with the change of groove depth in Figure8(c). It can be found that with the increase of groove depth, the dominant pole will pass through the imaginary axis in the complex plane, and then continue to move to the upper left in the second quadrant. During the movement of the dominant pole, the damping angle β increases slightly. This means that when the groove depth is large, there is a large saltation peak after the floating ring is disturbed. However, since the distances from the dominant pole to both the imaginary axis and the real axis increase, the system has better damping degree and response speed. The angular vibration of the floating ring will attenuate rapidly, and the amplitude of the subsequent oscillation is smaller. The gas seal system can return to steady state more quickly.Based on the direct numerical simulation method, the complex nonlinear flow field between the seal faces and the two DOF angular motion of the seal rings are coupled to obtain the complete time history of angular vibration of the floating ring, and the angular vibration stability of DGS under different groove depths is analyzed. An initial swing angular velocity of 1.88 × 104 μrad/s is applied to the stator of DGS in the y direction. Thus, the stator is excited and begins to vibrate in both x and y directions. Figure 9 shows the time domain response of the stator in the y direction to the initial swing angular velocity under different groove depths. When hg = 3 μm, the stator response under the initial velocity disturbance is a divergent response process (as shown in Figure 9(a)). The angular response displacement of the stator cannot reach a steady state and gradually diverges with the passage of time. Within just three rotation cycles, the angular displacement φy has reached 80 μrad. As can be seen from Figure 9(b), when hg > 5 μm, the gas seal system is in a stable state, and the angular response of the stator is an attenuated oscillation. By comparing the vibration peaks and the time to restore steady state of each response curve, it can also be found that with the increase of groove depth, the first peak value of the angular response displacement of the stator increases, but both the intensity and response time of the subsequent oscillation obviously decrease. In other words, the angular vibration stability of the system is improved. The comparative analysis with Figure 8 reveals that the stability range predicted based on the algebraic stability criterion and the dynamic performance obtained by the root locus method are basically consistent with the variation law of angular vibration characteristics of DGS with groove depth shown in Figure 9, which proves the feasibility, effectiveness and accuracy of using the methods of control theory to study the angular free vibration stability of DGS.Figure 9
Angular response displacement of the floating ring solved by DNSM under different groove depth: (a)hg = 3 μm and (b)hg > 5 μm.
(a)(b)In summary, combined with the direct numerical simulation results of the stator transient response, the validity of the method for judging the angular free vibration stability of DGS is verified. The results show that the rules and conclusions obtained by the two research methods are consistent, and the method presented in this paper can effectively predict the angular vibration stability of the DGS system and its dynamic performance in the transition stage. There is no doubt that the most accurate research method is to solve the complete angular vibration behavior of the floating ring. However, DNSM has a large amount of calculation and storage, and each calculation is only for the specified parameters, which is not suitable for seal design and structure optimization. In this paper, the linearization method is used to represent the features of gas film between the seal faces, and the angular dynamic performance of DGS are directly analyzed according to the system characteristic equation, without real-time solution of the flow field. Therefore, the calculation efficiency is much higher, and it is suitable for large-scale parameter discussion. In addition, the definition and calculation formula of critical rotational inertiaIcr are given in this paper. It is generally known that the larger the film angular main damping cxx is, the slower the response speed of the system is. Thus, Icr can more fully reflect the angular free vibration characteristics of DGS (including damping degree and response speed). As a comprehensive performance indicator, Icr can be used as an optimization objective of the angular dynamic performance of DGS in future research. On the premise that the main damping is greater than zero, the larger the Icr − I is, the more excellent the system transient performance in the transition stage is.


## 4. Conclusions

(1)
The constraint condition of rotational inertia is the most demanding of all stability conditions, which determines the dynamic characteristic coefficient of gas film in the critical instability state of the gas seal system.(2)
If the angular main stiffness is considerable, although the absolute stability of the system can be guaranteed, the floating ring of DGS will be in a high-frequency and large-amplitude swing state in a short time after being instantaneously excited in any angular direction.(3)
The approach of the absolute value of the cross coefficient to zero is beneficial to the suppression of angular free vibration.(4)
There is an obvious interaction between rotational inertia and stiffness or damping parameters. The four stiffness and damping parameters will all change the influence law of rotational inertia on the stability of angular free vibration. Except for the cross stiffness, the rotational inertia does not change the variation law of max(Re) with the other three parameters. There is a threshold of rotational inertia, which makes the cross stiffness change from having only the lower limit value to both the upper and lower critical values at the same time.



---

*Source: 7378999-2022-04-01.xml*