A study on the nonlinear vibration of the generator rotor based on the unbalanced electromagnetic force and the oil film force coupling model

2.1. Nonlinear motion equations

The rotor is supported by two same cylindrical sliding bearings. The nonlinear oil film force of the bearing is directly obtained by solving the Reynolds equation, and the effect of the temperature-viscosity relation is considered. The motion equation, applied to the rotor-bearing system, is built to take into account the effect of viscosity and elasticity of the rotor. The model of the symmetric viscoelastic rotor bearing system is shown in Figure 1, and a linear viscoelastic model and the lumped mass method are adopted to solve the motion equation of the system.

Fig. 1Symmetric viscoelasticity rotor-bearing system model

The position coordinates of the rotor and the rigid center of the disc are defined as $(x_1,y_1)$ and $(x_2,y_2)$, and the system satisfies the following motion equation:

where, $f_x$ and $f_y$ are the oil film force components in the horizontal and vertical directions. $g$ is the gravitational acceleration. $\omega$ is the rotor angular speed. $e$ is the eccentricity distance of the eccentric mass located on the rigid disc. ${\phi }_0$ is the phase angle. $K$ and $c_d$ represent the elasticity and damping coefficients of the rotor respectively, and $2k=48EI/l^3$, $2c_d=K_{0}48EI/l^3$ ($l$ is the center distance of a pair of bearings; $E$ and $K_0$ are the material constants; and I is the equatorial inertia moment of the rotor cross section). $m$ and $2M$ correspond to the concentrated mass of the rotor center and the disc center respectively, and $m=1/4m_1$, $2M=1/2m_1+m_2$ ($m_1$ and $m_2$ represent the total mass of the rotor and the disc respectively).

Applying the following dimensionless quantities to Eq. (1), we can get the dimensionless motion equation of the rotor system as shown in Eq. (2). The dimensionless quantities include: vibration component $X=x/c$, $Y=y/c$; oil film force components $F_X=f_x{\psi }^2/\left({\mu }_0\omega LD\right)\mathrm{}$, $F_Y=f_y{\psi }^2/{\mu }_0\omega LD$; mass $M_S=mc\omega {\psi }^2/{\mu }_{0}LD$, $M_B=Mc\omega {\psi }^2/{\mu }_{0}LD$; the unbalance eccentricity $\rho =e/c$; time $\tau =\mathrm{\omega }t$; and the shaft stiffness and the damping coefficient of the rotor $K=kc{\psi }^2/\left({\mu }_0\omega LD\right)$ and $C=c_{d}c{\psi }^2/\left({\mu }_{0}LD\right)$ ($L$, $D$, $c$, $\psi$ and ${\mu }_0$ represent bearing length, diameter, radius gap, the gap ratio, and lubricating oil viscosity respectively):

As is obvious in Eq. (2), for the balancing rotor, in case of no external excitation (i.e. $\rho =0$), when the rotor system is in the static equilibrium position, the resultant oil film forces in the $X$ direction and the $Y$ direction correspond to zero, and the total mass of the rotor and the disc respectively. The static equilibrium position, unaffected by the disk location, is contingent only on the total mass of the rotor. On the other hand, for the unbalancing rotor with external excitation (i.e. $\rho \ne 0$) such as the electromagnetic force generated by the inter-turn short circuit fault in the current study, the unbalance eccentricity $\rho$ can obtain. The electromagnetic force $F_e$ is expressed as:

where $n$ is the node number of the modeling rotor. $f_{ei}\left(i=1,2,3,\dots ,2n\right)$ are the components of the elecromagnetic force acting on the node, and its situation and coordinate direction are consistent with that of the node of the rotor.

2.2. The nonlinear oil film force model

For the sliding bearing supporting the rotor, the nonlinear oil film pressure distribution can be obtained by solving the dynamic Reynolds equation. Whereas the relevant bearing parameters: $\epsilon ,\theta ,{\epsilon }^{\text{'}},{\theta }^{\text{'}}$ influence the oil film pressure distribution, they are closely related to the rotor motion parameters $(X_1,Y_1,X_1^{\text{'}},Y_1^{\text{'}})$, as:

The dynamic Reynolds equation for obtaining the bearing dimensionless oil film pressure is expressed as follows [17]:

where, $H$ is the dimensionless film thickness; $P$ is the dimensionless film pressure; $\lambda$ is the axial dimensionless coordinates (the coordinate center is the center of the bearing, and the dimensionless unit is $L/2$); $\epsilon$ is the bearing eccentricity; $\theta$ is the bearing attitude angle. $\stackrel{-}{\mu }$ is the relative viscosity of the lubricating oil (22# turbine oil), and its relationship with the temperature of the lubricating oil is expressed in the following [18]:

The temperature distribution is consistent with the following energy equation [17]:

where, the relative temperature $\stackrel{-}{T}=T/T_0$$\left(T_0=\frac{2{\mu }_0\omega }{{\psi }^2}\frac{1}{\rho c_V}\right)$, and $\rho$ is oil density; $c_V$ is oil specific heat.

The finite element method is used to calculate the bearing oil film pressure distribution. The dimensionless finite element equation of the oil film pressure under small perturbation, derived from Galerkin Method based on Eq. (5), is expressed as:

where, $\mathrm{\Omega }$ is the oil film distribution area, $v$ is the shape function of the quadrilateral eight-node isoparametric element.

Substituting (4) into Eq. (8) and combining with Eq. (6) and (7), the dimensionless oil film pressure $P$ can be obtained.

The dimensionless oil film force components are derived as:

2.3. Calculation method for the periodic solution of the unbalanced rotor under the electromagnetic force

The vibration of the rotor system supported by sliding bearings under the unbalanced force is described as:

where $x$ consists of the displacement vectors and velocity vectors of the nodes of the modeling rotor. The unbalanced electromagnetic force generated by the inter-turn short circuit fault and the oil film force are included in$f\left(x,t\right)$. The coordinate center of displacement vectors is defined at the static equilibrium point of the rotor. Supposing that the period of the unbalance force is $T$, calculation of the periodic solution of the rotor system described above in Eq. (10) is equivalent to that of the solution of the following fixed boundary value problem with period $T$:

Eq. (10) can be reduced to one dimension using Poincare’s section because of its autonomous characteristics, and the Poincare’s section can be defined by setting the velocity component at 0. In this paper, the section where the horizontal displacement is equal to the vertical displacement and both are less than zero, is defined as the Poincare’s section, namely, $\sum =\{x_1=y_1<0\},\sum \in R^{n-1}$.

To help obtain the $x\left(0\right)$ in Eq. (11), we introduce the parameters vector $s=x\left(0\right)$. Once the value of $x\left(0\right)$ reaches certain level of precision, the periodic solutions of Eq. (11) can be obtained by directly applying the routines of numerical integration to solve the algebraic equations, as show in Eq. (12):

where, $x(T,s)$ is obtained by the following equation, as:

In this paper, a 4-stage Runge-Kutta method of order 4 is used to obtain the $x\left(T,s\right)$, and the integral step length is 1/50 of the period of the unbalanced exciting force. The convergence accuracy is 5 %.

The transformation from the Eq. (11) to Eq. (12) is beneficial to improving computation precision and speed, and a variety of methods to solve Eq. (12) are known as the shooting method [19]. In this paper, a quasi-Newton method for solving algebraic equations (12) is proposed, which has a wider range of convergence than otherwise. With the quasi-Newton method, the algebraic equations (12) can be solved even if the initial value is not a well-chosen one. In distinction from the general quasi-Newton shooting method, the method proposed in this paper is termed the improved quasi-Newton shooting method. We can rewrite the algebraic Eq. (12) with this shooting method with a variable parameter ${\beta }_k$, into the following iteration form:

where $k$ is the number of the iteration step; ${\beta }_k$ is a variable parameter in $(0,1]$ for the $k$th step; and $\Delta r\left(s\right)$ is the difference matrix of $r\left(s\right)$.

In $\Delta r\left(s\right)=\left({\Delta }_{1}r\right(s),{\Delta }_{2}r(s),\cdots ,{\Delta }_{s}r(s\left)\right)$, each component can be read as:

With the components $s_1=s_2,r_1\left(s\right)=r_2\left(s\right)$, the first row (column) of $\Delta r\left(s\right)$ is equal to the second row (column), and $\Delta r\left(s\right)$. will be solved by calculating the initial value problem (Eq. (13)) n times after is obtained. The variable parameter range ${\beta }_k$ is of $\left(0,1\right]$. In this paper, with ${\beta }_k=1$, $s^{\left(k+1\right)}$ is solved with Eq. (14). If $\Vert r\left(s^{(k+1)}\right)\Vert \ge \Vert r\left(s^{\left(k\right)}\right)\Vert$ (“$‖‖$” is the norm), then ${\beta }_k$ will be halved and the solution will be recalculated, until $\Vert r\left(s^{(k+1)}\right)\Vert <\Vert r\left(s^{\left(k\right)}\right)\Vert$.

The improved quasi-Newton shooting method in this paper has a wide range of convergence, and the number of iterations is less than that of the Newton shooting method. Nonetheless, it has to calculate the new difference quotient matrix $\Delta r\left(s\right)$ for each iteration so that when $n$ is quite large, the time cost will be very high. As a result, in this paper, we also present a new quasi-Newton method for solving the algebraic equation (12), which does not calculate the difference quotient matrix, with the Broyden rank being 1. With this method, although the time cost in each iteration is higher than that of the Newton method, it has obvious advantages, especially in case of a large $n$. This new quasi-Newton method with Broyden rank being 1 takes a similar form to Eq. (14), as in:

where ${\beta }_k$ is the parameter introduced in this paper, with the same purpose and usage as that in Eq. (14). Vector $\stackrel{-}{s}$ is computed by subtracting the first element from, and similarly vector $\stackrel{-}{r}\left(\stackrel{-}{s}\right)$ is computed by subtracting the first element from $r\left(s\right)$. Furthermore, iterative Eq. (14) can be rewritten as the iterative formula concerning $\stackrel{-}{s}$, which will facilitate the discussion on the stability of the system in the next section. $A_k$ is n-dimensional matrix introduced in reference [20]. When $k=0$, $A_k$ is a unit matrix; when $k>0$, it can be solved via the following formula:

To distinguish it from the general Broyden shooting method, the method introduced in this paper is called the improved Broyden shooting method. In brief, the convergence time cost of the two improved shooting methods is better than that of the general shooting method. In terms of the time cost, the improved quasi-Newton shooting method is preferred in the lower dimension system whereas the improved Broyden shooting method is favored in the higher dimension system.

When Eq. (12) is solved using the iterative Eq. (14) or (15), the convergence criteria is as follows:

where $\delta$ is a given accuracy. After $s$ and $r$, which meet the requirement of the aforementioned formula, are obtained, the periodic solution $g\left(t\right)$ (the period is $T$) will be calculated through Eq. (13).

2.4. Analysis method for the stability motion state of the rotor system

The Floquet theory is adopted to study the stability of the imbalance incentive periodic solution $g\left(t\right)$. For this purpose, it is necessary to solve the standard-based solution matrix $Y\left(T\right)$ of $\dot{y}=A\left(t\right)y$, $A(t+T)=A\left(T\right)$ with time $t$ being $T$; where $y$ is the disturbance vector $y=x-g\left(t\right)$, and $T$ is the period of the unbalanced exciting force,$A\left(t\right)={\frac{\partial f(x,t)}{\partial x}|}_{x=g\left(t\right)}$.

Solving $Y\left(t\right)$ of $\dot{Y}\left(t\right)=A\left(t\right)Y\left(t\right)$, $Y\left(0\right)=I$, is the initial value problem, and $Y\left(t\right)=I$ (unit matrix) under the $t=T$ condition. However, to obtain $Y\left(t\right)$, it is not necessary to directly solve the above equation. Instead, $Y\left(t\right)$ can be obtained by developing the relationship between $Y\left(t\right)$ and the difference quotient matrix of the aforementioned shooting method. The partial derivative of Eq. (13) is written as $\Delta r\left(s\right)$:

Then $Y\left(T\right)=\frac{\partial x(T,s)}{\partial s}$. From Eq. (12) we get $\frac{\partial x(T,s)}{\partial s}=I-\frac{\partial r}{\partial s}$. When replacing $\frac{\partial r}{\partial s}$ with the difference quotient matrix $\Delta r\left(s\right)$, we then obtain $Y\left(T\right)=I-\Delta r\left(s\right)$. Once a periodic solution that is precise enough is obtained, the difference quotient matrix can be recalculated. In doing so, $Y\left(T\right)$ can be obtained.

With the Floquet theory, the characteristic matrix $C=Y\left(T\right)$, the eigenvalues (characteristic coefficient) of the matrix $C$ can be determined by the following formula, as:

${\lambda }_i$ reflects the convergence characteristics of periodic solutions under the unbalanced incentive condition. The real part ${\lambda }_i$ is uniquely determined by formula (19), but the imaginary part is not uniquely determined. Assuming ${\lambda }_i=-u_i+j\left({\omega }_i+\frac{2\pi k}{T}\right)$, where $i=\mathrm{1,2},\cdots ,n$; $k=0,\pm 1,\pm 2,\cdots ,$ and $j=\sqrt{-1}$, the periodic solution stability margin indicator, namely, the characteristics damping, is presented, as:

Based on the Floquet theory, the stability criteria of the unbalanced incentive periodic solution, namely, $g\left(t\right)$ is described as follows: When $U_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is greater than zero, then it is in the stable state; when $U_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is less than zero, then it is in the unstable state; and when $U_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is equal to zero, it is in the critical state.

The whole calculation flow chart of the nonlinear dynamics coupling the nonlinear oil film force and the electromagnetic force induced by the inter-turn short circle fault is show in Fig. 2.

Fig. 2Nonlinear dynamics flowchart

The quantitative motion state of the instable rotor system and its convergence are determined by the time $T$ Poincaré mapping method. The procedures are as follows:

1) Define $Z_1=X_2-X_1$, $Z_2=Y_2-Y_1$, $Z_3=X_2^{\text{'}}$, $Z_4=Y_2^{\text{'}}$, where, ' represents $\frac{d}{d\tau }$; $\left(X_1,Y_1\right)$ is the initial bearing journal center position, as show in Eq. (2), also, the static equilibrium position.

2) Make $Z={(Z_1,Z_2,Z_3,Z_4)}^T$, and take $Z\left(0\right)={(\mathrm{0,0},\mathrm{0,0})}^T$ as the starting point. We can obtain $Z\left(kT\right),k=\mathrm{1,2},\cdots \cdots ;T=2\pi$ by solving Eq. (2). According to the $T$ Poincaré definition [19], the mapping point $g^k\left(Z\right(0\left)\right)=Z\left(kT\right)$, and the orderly mapping point column $g^k\left(Z\left(0\right)\right),(k=\mathrm{1,2},\cdots )$ are obtained, and $\Vert Z\Vert =\sum _{i=1}^4\left|Z_i\right|$ is taken as the flag value.

3) Check the convergence precision $R_i$ of the possible branch flag value $\Vert Z\Vert$ and determine if it has met the convergence requirement, as:

where $T$ is the period of the excitation force; $\epsilon$ is the convergence precision (5 % in this paper). When the mapping point is basically “fixed”, we can determine that the rotor system arrives at the steady motion state [20].

3.1. Calculation object and processes

A Dongfang 300 MW generator bearing-rotor system is taken as the research object. The main geometrical parameters and working parameters of the two same supporting bearings are shown in Table 1. The main parameters of the rotor are listed in Table 2. The model of the modeling rotor is consisted of 9 nodes and 8 segments, as shown in Figure 3. The generator rotor model and its fault positions (0° position, 90° position and 180° position) are described in Figure 4. In this paper, we will simulate the nonlinear vibration with the serious 90° position inter-turn short circuit fault (the affected turns are about 25 % of the total turns) when the rotor is at the rated speed of 3000 rpm. The electromagnetic force and the corresponding imbalance magnitude generated by the inter-turn short circuit fault are given in Table 3.

Fig. 3Dynamic model for the generator rotor

Fig. 4Generator model and diagram for the short circuit location

a) Turbine generator model

b) Diagram for the inter-turn short circuit location

4.1. Calculation process and convergence analysis

The characteristics damping $U_{min}$ of the balance rotor in its normal state (at the rated speed of 3000 rpm) is 0.0695; and the critical unbalanced excitation force of the rotor at its rated speed is approximately 8747.50 N ($U_{min}=0$). The range of the unbalanced electromagnetic force generated by the inter-turn short circuit fault is between 12213 and 127110 N. Therefore, it can be concluded that the rotor will fall into the unstable state when the fault occurs. Given that the motion state of the rotor in case of instability varies with the electromagnetic force value, we will investigate the nonlinear vibration in one specific fault mode (90° position short circuit serious fault). Figure 5 shows the convergence process. The flag $\Vert Z\Vert =\sum _{i=1}^4\left|Z_i\right|$ changes with the solution period. As shown in Fig. 5, when the 120 T period is reached, the convergence precision given in Eq. (21) $R=$4.37 % < 5 %, and thus the calculation is completed.

4.2. Discussion

The coordinates of the horizontal and vertical displacements $(X,Y)$ and oil film force $(F_x,F_y)$ are shown in Figure 6. In Fig. 6, $W$ is the bearing loads; $X$, $Y$, $F_x$ and $F_y$ are dimensionless parameters; the relative units of $(X,Y)$ are the side clearance of bearings $C$, the relative unit of $(F_x,F_y)$ is $W$. The results obtained in this paper include the displacement periodic solutions of nodes (2, 3, 4, 5 shown in Fig. 3) which are read as $X_1$, $Y_1$, $X_2$, $Y_2$, $X_3$, $Y_3$, $X_4$, $Y_4$, and the waveform of the oil film force $(F_x,F_y)$. The horizontal ordinate in the following figures represents the dimensionless time period $T$ in radians, and when the unit is second, $T=\frac{\pi \Omega t}{30}$. The maximum value of the horizontal ordinate is the cycle number of the stable periodic solutions. The vertical ordinate of the nonlinear vibration response and the rotor centerline orbit $(X,Y)$ are the relative vibration amplitude between the actual vibration response amplitude and the bearing side of the gap $C$. The vertical ordinate for the transient oil film force waveform $(F_x,F_y)$ is the ratio of the oil film force to the bearing loads.

Fig. 5Motion convergence process

Fig. 6Coordinate system (O is the center of bearing)

Figure 7 shows the rotor centerline orbit and the nonlinear vibration response generated by the electromagnetic force under the 90° position serious short circuit fault condition. Figure 8 is the waveform diagram of the transient process of the oil film force.

The FFT for the vibration response in the steady motion state is completed under the 90° position serious fault condition, and the maximum, minimum displacements and amplitude of the steady-state solution (for nodes 2, 3, 4, 5 and oil film force) are listed in Table 4. The harmonic components of the vibration response in nodes 2, 3, 4, 5 and the oil film force are shown in Table 5.

Fig. 7Rotor centerline orbit under 90° position serious fault

a) Node 2

b) Node 3

c) Node 4

d) Node 5

Fig. 8Oil film force under 90° position serious fault

a) Horizontal direction

b) Vertical direction

As is obvious in Figures 7 and 8, the steady motion state of the rotor system is a synchronized vortex. Figure 7 shows the rotor centerline orbit of the periodic motion. The results in Tables 4 and 5 show that the system motion is the constant periodic motion (the period is $T$), and the main motion state is a synchronized vortex. The vibration responses in Node 2 (the center of the bearing) with the $\Omega$ synchronizing component in the horizontal and vertical directions are 0.376 and 0.168 respectively; and the horizontal and vertical vibration sync components in the Node 5 (rotor symmetrical center) are 0.169 and 0.282, respectively. The $\Omega$ synchronizing components of the oil film forces in the horizontal and vertical directions are 0.346 and 0.546 respectively. Meanwhile, as can be seen from Tables 4 and 5, due to the electromagnetic force induced by the inter-turn short-circuit fault, the vibration amplitudes of the rotor system in Node 2, 3, 4 and 5 are excited. Compared with the normal case, the vibration amplitudes in the horizontal and vertical directions of Node 2, 3, 4 and 5 increase to 425 %, 168 %, 80 %, 28 % and 2702 %, 186 %, 164 %, 129 % respectively. The horizontal and vertical oil film forces effected by the electromagnetic force increase to 367 % and 728 %, respectively.

4.3. A comparison with the site test results

The above research object is from a power plant in Huangshi City, Hubei Province, China. The tested 300 MW generator unit is labeled as No. 2, which is used with the cooling fluid of water - hydrogen. The bearings supporting the generator rotor are labeled as the #7, #8 bearing respectively. Nodes 2 and 8 Figure 3 indicates the position, where the bearings are assembled. #7 and #8 bearings are at nodes 2 and 8 respectively. In 2005, an emergency stop occurred on the generator units, and the loss excitation protection action aroused. The duration of the asynchronous loss excitation conditions was about 14 s, and the current in the rotor was up to 16 kA. After re-booting, we found that the vibration in #7, #8 bearings was abnormal and the $\Omega$ (running frequency) component vibration was intensified. Therefore, we conducted a reactive power test. During the test, the active power was maintained at 107 MW, and the reactive power increased from 42 MVar, to 70 MVar, and finally to 75 MVar. The increase of the reactive power reflected that of the rotor current. The relation between vibration amplitude in the bearing (#7 bearing) and the reactive power is shown in Figure 9. From Fig. 9, we can see that the vibration in the horizontal and vertical directions increased with the reactive power.

Fig. 9Rotor vibration amplitude vs. reactive power

After the shutdown of the generator unit, the detected insulation of the excitation winding to ground was $0\Omega$. AC impedance value of the excitation winding was reduced to about 40 % of the normal unit. Voltage was applied to two slip rings of the rotor, and voltage values measured at two slip rings to the ground were significantly different. The implication was that the sudden failure of the generator unit was due to the inter-turn short circuit fault. For further examination, the generator unit was dismantled, and we found that the ground insulation was breakdown at approximately 90° (shown in Figure 4) position of the rotor symmetrical center. The fault led to the direct inter-turns short circuit from the upper incoming coil turn #30 slots to the upper outlet coil turn #29 slots, and the total number of the short circuit coil turns was up to 25 % of that of the whole field winding turns. Undoubtedly, it was a very serious inter-turn short circuit fault. The photo on the fault position of the dismantled generator rotor is shown in Figure 10.

Fig. 10Photo of field winding inter-turn fault in generator rotor

The results in section 4.2 show that vibration amplitude (Node 2, the bearing situation) in the horizontal direction (X direction) is 142 μm (0.4166×0.675/2 = 142 μm), and the vertical direction (Y direction) is 66 μm (0.1962×0.675/2 = 66 μm) (as shown in Table 4). The test results in Figure 9 show that the measured vibration amplitudes in the horizontal and vertical directions were about 138 μ m and 71 μm respectively. Because the error between theoretical prediction and test results is lower than 10 %, the current study is of high engineering value.