Research on Rotordynamic Characteristics of Pump Annular Seals Based on a New Transient CFD Method

: Pump annular seals can cause ﬂuid reaction forces that have great e ﬀ ects on the vibration characteristic and stability of a pump system. For this reason, it is important to study rotordynamic characteristics of annular seals. In this paper, a new transient computational ﬂuid dynamics (CFD) method with dynamic mesh is proposed to investigate rotordynamic characteristics of the pump annular seal. The reliability of the transient CFD method is validated by comparison with the results from the experiment and the bulk-ﬂow method, and the relationship between the seal length and rotordynamic characteristics is investigated by the transient CFD method. The results indicate that direct sti ﬀ ness decreases sharply even turns to negative as the seal length increases, this phenomenon may change the direction of ﬂuid force on the rotor surface and a ﬀ ect supporting condition of the pump rotor. With the increasing seal length, the whirl frequency ratio gradually increases, which would weaken the stability of the pump rotor system.


Introduction
Liquid annular seals are necessary components of pumps to control leakage flow, improve efficiency and balance axial thrust. Annular seals mainly include impeller wear rings, inter-stage seals and balancing drum gaps. The length of an impeller wear ring is short and the aspect ratio (L/D) is generally around 0.1. The aspect ratios of inter-stage seals and balancing drum gaps can reach up to 3. Just like journal bearings, small clearances, large axial pressure drops and shearing effects of pump rotors mean that annular seals have an important effect on the vibration and stability characteristics of a rotor system.
In view of the important influence of annular seals on vibration characteristics of pump systems [1], the bulk-flow method is commonly used in industry to study the dynamic characteristics of annular seals. The bulk-flow method is based on the thin film assumption in Hirs' turbulent lubrication theory, which is actually a two-dimensional flow model as it ignores the radial gradients of flow variables. Childs [2] first achieved the finite length solutions of dynamic coefficients of annular seals based on the bulk-flow method. Then Dietzen et al. [3], Marquette et al. [4] and Sun et al. [5] used the bulk-flow method to investigate the static and dynamic characteristics of annular seals. However, the bulk-flow method lacks flow details and its accuracy largely depends on empirical correction coefficients such as friction factors that will change with seal geometry and operating conditions. With the development of computer technology, a CFD method [6] has been used as an effective method for predicting the leakage flow rates and rotordynamic characteristics of annular seals. Compared to the bulk-flow method, the CFD method is better at describing internal flow field, and is applicable for more complex sealing structures and operating conditions. Tam et al. [7] applied the quasi-steady CFD method to study fluid dynamic forces in seals and bearings. They assumed that the rotor performed the circle whirl motion around the seal center and solved the whirl flow field in the rotating frame of

Method for Determining Rotordynamic Coefficients
As the rotor location changes with time in a rotor-stator system, a transient CFD simulation is applied to investigate the unsteady flow within a full 360-degree model of the seal. Before starting the transient simulation, the rotor is at the seal center as shown in Figure 1a. During the simulation, the rotor position is perturbed by using a circle orbit whirl model, and the maximum clearance and the minimum clearance will appear in the seal as shown in Figure 1b. The circle whirl model consists of two types of rotor whirling motions, one is forward whirling where the rotor whirls and rotates in the same direction, and the other is backward whirling where the rotor whirls and rotates in opposite directions. In this circle orbit whirling model, the moving displacements of the rotor in the x and y directions are given in Equation (1), where the amplitude e denotes the whirl radius.
x d = ecos(Ωt), y d = esin (Ωt) (Forward whirl) x d = ecos(Ωt), y d = −esin(Ωt) (Backward whirl) (1) Processes 2020, 8, x FOR PEER REVIEW 3 of 15 opposite directions. In this circle orbit whirling model, the moving displacements of the rotor in the x and y directions are given in Equation (1), where the amplitude e denotes the whirl radius.
During the transient computations, the fluid forces (Fx, Fy) acting on the rotor at each time step are captured by integrating the pressure and shear stresses. In order to quantify the fluid forces for small perturbation around a concentric rotor position, the conventional rotordynamic model is expressed in Equation (2) where K and k are direct and cross stiffness coefficients, respectively; C and c are direct and cross damping coefficients, respectively; and M is a direct mass coefficient. According to the relationship between (Fx, Fy) and (Fr, Ft), the dynamic coefficients can be expressed as Equation (3).

Geometry Model and Grid
The annular seal used in the present study is shown in Figure 2a, and the main geometric parameters and operating conditions are shown in Table 1. The generated grid can be seen in Figure 2b. The whole flow field of the seal is meshed with structured hexahedral grids and the grid independence investigation under different operating pressures is carried out. Figure 3 shows the grid independence investigation result under the pressure of 1.38 MPa. It is found that the leakage remains almost constant when the mesh number reaches approximately 0.425 million. Therefore, the mesh used in this paper is composed of 424,960 nodes and 364,980 cells. During the transient computations, the fluid forces (Fx, Fy) acting on the rotor at each time step are captured by integrating the pressure and shear stresses. In order to quantify the fluid forces for small perturbation around a concentric rotor position, the conventional rotordynamic model is expressed in where K and k are direct and cross stiffness coefficients, respectively; C and c are direct and cross damping coefficients, respectively; and M is a direct mass coefficient. According to the relationship between (Fx, Fy) and (Fr, Ft), the dynamic coefficients can be expressed as Equation (3).

Geometry Model and Grid
The annular seal used in the present study is shown in Figure 2a, and the main geometric parameters and operating conditions are shown in Table 1. The generated grid can be seen in Figure 2b. The whole flow field of the seal is meshed with structured hexahedral grids and the grid independence investigation under different operating pressures is carried out. Figure 3 shows the grid independence investigation result under the pressure of 1.38 MPa. It is found that the leakage remains almost constant when the mesh number reaches approximately 0.425 million. Therefore, the mesh used in this paper is composed of 424,960 nodes and 364,980 cells.

Mesh Movement
Due to the circle orbit whirl model of the rotor as mentioned in Section 2, the shape of the fluid domain changes with time, which needs a transient CFD method to deal with it. To perform the transient CFD computations, the dynamic mesh problem should be solved. In this work, the dynamic mesh is achieved by applying a user-defined subroutine linked with FLUENT solver to strictly control the motion of every grid node in the seal model. The specific implementation procedures are as follows. Figure 4 shows a cross-section view of a simplified plain seal at an arbitrary axial position. As illustrated in the figure, the seal clearance Cr has been exaggerated, C1 represents the current eccentric position where the rotor moves from concentric position marked C0, d (xd, yd, 0) denotes the moving distance of the rotor in Cartesian coordinates. When the rotor is at the concentric position, Pi0 represents an arbitrary node in the clearance of the seal and is assumed to project along a mesh line that connects one node Pr0 on the rotor surface and another Ps0 on the stator surface. After the rotor moves by d, the movement displacements of the nodes on the rotor surface are also d. Therefore, the new coordinates (xr1, yr1) of Pr0 are defined as Equation (4).

Mesh Movement
Due to the circle orbit whirl model of the rotor as mentioned in Section 2, the shape of the fluid domain changes with time, which needs a transient CFD method to deal with it. To perform the transient CFD computations, the dynamic mesh problem should be solved. In this work, the dynamic mesh is achieved by applying a user-defined subroutine linked with FLUENT solver to strictly control the motion of every grid node in the seal model. The specific implementation procedures are as follows. Figure 4 shows a cross-section view of a simplified plain seal at an arbitrary axial position. As illustrated in the figure, the seal clearance C r has been exaggerated, C 1 represents the current eccentric position where the rotor moves from concentric position marked C 0 , d (x d , y d , 0) denotes the moving distance of the rotor in Cartesian coordinates. When the rotor is at the concentric position, P i0 represents an arbitrary node in the clearance of the seal and is assumed to project along a mesh line that connects one node P r0 on the rotor surface and another P s0 on the stator surface. After the rotor moves by d, the movement displacements of the nodes on the rotor surface are also d. Therefore, the new coordinates (x r1 , y r1 ) of P r0 are defined as Equation (4).
where x r0 , y r0 are x and y coordinates of P r0 when the rotor is at the concentric position. The nodes of the stator surface stay still, and the movement distances of the nodes between the rotor and the stator are determined by the interpolation algorithm [24] as shown in Equation (5a). Therefore, the new position coordinates (x i1 , y i1 ) of P i0 are transformed using Equation (5b) where x i0 , y i0 are x and y coordinates of P i0 when the rotor is at the concentric position. In Equation (5a), ra represents the ratio of the distance between the nodes in the clearance and the outer wall to the clearance. When the rotor is in the concentric and eccentric positions, ra can be expressed as Equation (6) where θ denotes the initial angular coordinate of the node P i0 as shown in Figure 4. By Equation (7), θ is expressed as Equation (8). Substituting Equation (8) into Equation (6) and then substituting Equation (6) into Equation (5b) gives new positions of the nodes in clearances after the rotor moves. x = x +x , y = y +y (4) where xr0, yr0 are x and y coordinates of Pr0 when the rotor is at the concentric position. The nodes of the stator surface stay still, and the movement distances of the nodes between the rotor and the stator are determined by the interpolation algorithm [24] as shown in Equation (5a). Therefore, the new position coordinates (xi1, yi1) of Pi0 are transformed using where xi0, yi0 are x and y coordinates of Pi0 when the rotor is at the concentric position. In Equation (5a), ra represents the ratio of the distance between the nodes in the clearance and the outer wall to the clearance. When the rotor is in the concentric and eccentric positions, ra can be expressed as Equation (6) where θ denotes the initial angular coordinate of the node Pi0 as shown in Figure 4. By Equation (7), θ is expressed as Equation (8). Substituting Equation (8) into Equation (6) and then substituting Equation (6)   By adopting the dynamic mesh algorithm, the displacement of each grid node is strictly calculated by the mathematical procedures to ensure the movement coordination of adjacent grid nodes. The algorithm has been tested and found that when the rotor whirled from the concentric position as shown in Figure 5a to the eccentric position with an amplitude of 90% seal clearance as shown in Figure 5b, the maximum skewness of the grids increases from 0.065 to 0.46. The grid skewness becomes larger, which means that when the rotor whirls with a big eccentricity, the grid distortion rate will increase (grid skewness ranges from 0 to 1, where the value is close to 1, which corresponds to low quality). However, the grids in the clearance still do not appear as highly distorted elements and negative volumes, which indicates that the dynamic mesh algorithm is applicable for transient flow simulations of the seal. During the rotor whirling, the locus of the rotor center changes with time so that the rotational speed of the rotor cannot be set by traditional methods. The DEFINE macro DEFINE_PROFILE is applied to read the coordinates of the rotor grids, which change according to the instantaneous position of the rotor. Tangential velocity in the surface of the rotor is divided into two component velocities as given in Equation (9), which are given to the surface nodes of the rotor by the DEFINE macro [27]. where v x , v y are velocities in the x and y directions of an arbitrary grid node of P 2 , and ω is the rotational speed of the rotor.
methods. The DEFINE macro DEFINE_PROFILE is applied to read the coordinates of the rotor grids, which change according to the instantaneous position of the rotor. Tangential velocity in the surface of the rotor is divided into two component velocities as given in Equation (9), which are given to the surface nodes of the rotor by the DEFINE macro [27].
where vx, vy are velocities in the x and y directions of an arbitrary grid node of P2, and ω is the rotational speed of the rotor.

Numerical Method
The commercial code FLUENT16.0 is employed to solve the three-dimensional incompressible Navier-Stokes equations. Figure 6 shows the flow chart for the calculation of unsteady CFD. A steady-state CFD solution is initially solved to obtain the initial conditions of the unsteady CFD calculation. In the unsteady CFD calculation, the rotational speed of the rotor is defined by Define macro Define_ profile, and the whirling speed of the rotor is loaded by user-defined function. The inlet boundary is set as total pressure and the outlet boundary is set as static pressure. Nonslip boundary conditions are used for the near-wall function. The realizable k-ε model and the enhanced wall function are chosen according to the wall Y+, which is located in the range of 20-40 [17]. The size of the time step is 1.634E-5s, i.e., the time spent by the rotor in one degree rotation. The whirl radius e is 0.0076 mm, which is selected as 10% of the radial clearance. Second-order, upwind discretization has been used for the convection and central difference schemes for diffusion terms. The SIMPLE algorithm is used to deal with the pressure-velocity coupling.

Numerical Method
The commercial code FLUENT16.0 is employed to solve the three-dimensional incompressible Navier-Stokes equations. Figure 6 shows the flow chart for the calculation of unsteady CFD. A steady-state CFD solution is initially solved to obtain the initial conditions of the unsteady CFD calculation. In the unsteady CFD calculation, the rotational speed of the rotor is defined by Define macro Define_ profile, and the whirling speed of the rotor is loaded by user-defined function. The inlet boundary is set as total pressure and the outlet boundary is set as static pressure. Nonslip boundary conditions are used for the near-wall function. The realizable k-ε model and the enhanced wall function are chosen according to the wall Y+, which is located in the range of 20-40 [17]. The size of the time step is 1.634E-5s, i.e., the time spent by the rotor in one degree rotation. The whirl radius e is 0.0076 mm, which is selected as 10% of the radial clearance. Second-order, upwind discretization has been used for the convection and central difference schemes for diffusion terms. The SIMPLE algorithm is used to deal with the pressure-velocity coupling.    Figure 7 shows the calculated reaction forces acting on the rotor respectively with forward and backward whirl motion models under three pressure differences. It can be seen from the figure that the fluid forces of forward and backward whirls show similar periodic behaviors. The fluid reaction forces increase as the pressure difference increases, and the force Fx in the X direction is slightly larger than the force Fy in the Y direction.

Theoretical and Experimental Verification
The calculation results based on the proposed method are compared to the results based on the bulk-flow method and the experimental results. Figure 8a shows direct stiffness changes with respect to the pressure differences. As illustrated in the figure, the maximum error of predicted direct stiffness K under three pressure differences based on the transient CFD method is less than 20%, while the maximum error of K from the bulk-flow method is almost up to 40%. The result of K from the transient CFD method shows better correlation with the experimental result in comparison to the result from the bulk-flow method. Figure 8b shows cross coupled stiffness k. The results of the transient CFD method and the bulk-flow method are all larger than the experimental one. The errors of the transient CFD method from 1.38 MPa to 3.45 MPa are approximately 29%, 26% and 32% respectively. However, the errors of bulk-flow analysis under corresponding pressure differences are about 38%, 42% and 29% respectively, which yield a worse prediction. Figure 8c shows direct damping C. The experimental result shows that the errors of the transient CFD method and the bulk-flow method for C are about 4~10% and 4~19%, respectively. The leakage of the seal is obtained by integrating the Z-directional flow velocity with respect to the cross section of the seal inlet and is shown in Figure 9. The results from the CFD simulation and bulk-flow method are respectively 16% and 22% of the experimental results. Therefore, the transient CFD analysis provides better improvements than the analysis taken from the bulk-flow method. The listed comparisons show that the transient CFD method proposed in this paper can significantly improve the prediction accuracy of leakage rates and dynamic characteristics of annular seals. by integrating the Z-directional flow velocity with respect to the cross section of the seal inlet and is shown in Figure 9. The results from the CFD simulation and bulk-flow method are respectively 16% and 22% of the experimental results. Therefore, the transient CFD analysis provides better improvements than the analysis taken from the bulk-flow method. The listed comparisons show that the transient CFD method proposed in this paper can significantly improve the prediction accuracy of leakage rates and dynamic characteristics of annular seals.

Results and Discussion
The length of the seal has a great effect on dynamic characteristics, especially direct stiffness. The length of the seal was changed from 13.13 mm (L/D = 0.17) to 91.44 mm (L/D = 1.2) with other boundary conditions maintained in order to investigate the relationship between the length and dynamic characteristics. The results of dynamic characteristics were obtained by using the transient CFD method, as shown in Table 2. Figure 10 shows the relationship between direct stiffness and the seal length. It can be seen that direct stiffness changes from positive to negative as the seal length increases and a critical length exists where direct stiffness is zero. In the present study, the critical points are at the length of 50.05 mm, 63.24 mm and 74.32 mm, when the pressure differences are 1.38, 2.41 and 3.45 MPa. This phenomenon is mainly caused by inlet loss and the change of the static pressure distribution in the flow channel of the clearances. Figure 11 shows that axial velocity in the eccentric side (minimum clearance) is smaller than that in the reverse side (maximum clearance). Luis [28] pointed out that a small clearance has a smaller axial velocity and lower inlet loss, which would bring higher static pressure to the eccentric side of the seal than that in the reverse side, thus creating a fluid reaction force to push the rotor to the seal center. This phenomenon is known as Lomakin effects, which can produce positive direct stiffness. Chen [29] explained that inlet loss is the main cause for

Results and Discussion
The length of the seal has a great effect on dynamic characteristics, especially direct stiffness. The length of the seal was changed from 13.13 mm (L/D = 0.17) to 91.44 mm (L/D = 1.2) with other boundary conditions maintained in order to investigate the relationship between the length and dynamic characteristics. The results of dynamic characteristics were obtained by using the transient CFD method, as shown in Table 2. Figure 10 shows the relationship between direct stiffness and the seal length. It can be seen that direct stiffness changes from positive to negative as the seal length increases and a critical length exists where direct stiffness is zero. In the present study, the critical points are at the length of 50.05 mm, 63.24 mm and 74.32 mm, when the pressure differences are 1.38, 2.41 and 3.45 MPa. This phenomenon is mainly caused by inlet loss and the change of the static pressure distribution in the flow channel of the clearances. Figure 11 shows that axial velocity in the eccentric side (minimum clearance) is smaller than that in the reverse side (maximum clearance). Luis [28] pointed out that a small clearance has a smaller axial velocity and lower inlet loss, which would bring higher static pressure to the eccentric side of the seal than that in the reverse side, thus creating a fluid reaction force to push the rotor to the seal center. This phenomenon is known as Lomakin effects, which can produce positive direct stiffness. Chen [29] explained that inlet loss is the main cause for Lomakin effects. So if the seal is not long enough to ignore the influence of inlet flow, a positive direct stiffness can be obtained. However, static pressure in the eccentric side of the seal is lower than that in the reverse side due to higher velocity in the eccentric side as shown in Figure 12 when the seal is long enough to ignore the influence of inlet flow. Therefore, direct stiffness of the long seal can be negative.

Results and Discussion
The length of the seal has a great effect on dynamic characteristics, especially direct stiffness. The length of the seal was changed from 13.13 mm (L/D = 0.17) to 91.44 mm (L/D = 1.2) with other boundary conditions maintained in order to investigate the relationship between the length and dynamic characteristics. The results of dynamic characteristics were obtained by using the transient CFD method, as shown in Table 2. Figure 10 shows the relationship between direct stiffness and the seal length. It can be seen that direct stiffness changes from positive to negative as the seal length increases and a critical length exists where direct stiffness is zero. In the present study, the critical points are at the length of 50.05 mm, 63.24 mm and 74.32 mm, when the pressure differences are 1.38, 2.41 and 3.45 MPa. This phenomenon is mainly caused by inlet loss and the change of the static pressure distribution in the flow channel of the clearances. Figure 11 shows that axial velocity in the eccentric side (minimum clearance) is smaller than that in the reverse side (maximum clearance). Luis [28] pointed out that a small clearance has a smaller axial velocity and lower inlet loss, which would bring higher static pressure to the eccentric side of the seal than that in the reverse side, thus creating a fluid reaction force to push the rotor to the seal center. This phenomenon is known as Lomakin effects, which can produce positive direct stiffness. Chen [29] explained that inlet loss is the main cause for Lomakin effects. So if the seal is not long enough to ignore the influence of inlet flow, a positive direct stiffness can be obtained. However, static pressure in the eccentric side of the seal is lower than that in the reverse side due to higher velocity in the eccentric side as shown in Figure 12 when the seal is long enough to ignore the influence of inlet flow. Therefore, direct stiffness of the long seal can be negative.    Figure 13 shows the distribution of static pressure on the rotor face with different seal lengths for Ω/ω = 0. Lx is the axial position and L is the seal length. Lx/L = 0 and Lx/L = 1 denote the position of seal inlet and outlet, respectively. The red circles in Figure 13 represent cross points, which indicate the axial position where the static pressure of the eccentric side equals that of the reverse side. Static pressure of the axial position before the cross point is consistent with Lomakin effects, while after the cross point is consistent with Bernoulli effects. For the short seals as shown in Figure 13a,b, static pressures of the eccentric side and the reverse side show an almost linear reduction across the seal length, and a cross point on the axial position does not exist in the short seal. Static pressure in the eccentric side is indeed larger than that in the reverse side, especially at a larger pressure difference, which is in line with Lomakin effects, and its direct stiffness is positive. For the long seals as shown in Figure 13c,d, static pressure on both sides also decreases in the direction of the leak and a cross point appears in the axial position at three pressure differences, where Bernoulli effects play a dominant role in the seal. This is due to inlet flow gradually stabilizing as the seal length increases so that the influence of inlet flow becomes weak and makes Bernoulli effects become noticeable, which corresponds to negative direct stiffness.  Figure 13 shows the distribution of static pressure on the rotor face with different seal lengths for Ω/ω = 0. Lx is the axial position and L is the seal length. Lx/L = 0 and Lx/L = 1 denote the position of seal inlet and outlet, respectively. The red circles in Figure 13 represent cross points, which indicate the axial position where the static pressure of the eccentric side equals that of the reverse side. Static pressure of the axial position before the cross point is consistent with Lomakin effects, while after the cross point is consistent with Bernoulli effects. For the short seals as shown in Figure 13a,b, static pressures of the eccentric side and the reverse side show an almost linear reduction across the seal length, and a cross point on the axial position does not exist in the short seal. Static pressure in the eccentric side is indeed larger than that in the reverse side, especially at a larger pressure difference, which is in line with Lomakin effects, and its direct stiffness is positive. For the long seals as shown in Figure 13c,d, static pressure on both sides also decreases in the direction of the leak and a cross point appears in the axial position at three pressure differences, where Bernoulli effects play a dominant role in the seal. This is due to inlet flow gradually stabilizing as the seal length increases so that the influence of inlet flow becomes weak and makes Bernoulli effects become noticeable, which corresponds to negative direct stiffness.
The whirl-frequency ratio f is introduced as Equation (10) [30] to weigh the effect of the seal length on the stability of the pump rotor system From Figure 14, the whirl-frequency ratio generally increases with the increasing seal length, due to the influence of cross-coupled stiffness. As the seal length increases to 91.44 mm, the whirl-frequency ratios are nearly up to 0.47 when the seal operates under three pressure differences. This phenomenon means that the long seal would reduce stability of the rotor system. The whirl-frequency ratio f is introduced as Equation (10) [30] to weigh the effect of the seal length on the stability of the pump rotor system k f C = Ω (10) From Figure 14, the whirl-frequency ratio generally increases with the increasing seal length, due to the influence of cross-coupled stiffness. As the seal length increases to 91.44 mm, the whirl-frequency ratios are nearly up to 0.47 when the seal operates under three pressure differences. This phenomenon means that the long seal would reduce stability of the rotor system.   The whirl-frequency ratio f is introduced as Equation (10) [30] to weigh the effect of the seal length on the stability of the pump rotor system k f C = Ω (10) From Figure 14, the whirl-frequency ratio generally increases with the increasing seal length, due to the influence of cross-coupled stiffness. As the seal length increases to 91.44 mm, the whirl-frequency ratios are nearly up to 0.47 when the seal operates under three pressure differences. This phenomenon means that the long seal would reduce stability of the rotor system.

Conclusions
In this paper, the transient CFD method with dynamic mesh was used to predict dynamic characteristics of the pump annular seal. The following conclusions can be made.
(1) The proposed transient CFD method can control the grid movement at each time step and achieve optimal grid quality of the displacement grid at any time step. The CFD results were compared with the results from the experiment and the bulk-flow method. It was shown that the transient CFD analysis can provide better improvements than the bulk-flow analysis. The dynamic characteristics of annular seals can be predicted accurately by the transient CFD method. (2) The relationship between the seal length and rotordynamic characteristics was also investigated by the transient CFD method. The results show that direct stiffness changes from positive to negative as the seal length increases. This phenomenon can change the direction of the fluid force on the rotor surface and reduce the supporting effect of the annular seal on the pump rotor. (3) For the short seal, the static pressure of the eccentric side is almost larger than that of the reverse side and direct stiffness is positive. For the long seal, a cross point exists in the axial position that makes Bernoulli effects predominant and direct stiffness negative. With the increasing seal length, the whirl-frequency ratio becomes larger, which decreases the stability of the pump rotor system. LGG18E060007).

Conflicts of Interest:
The authors declare no conflict of interest.