Transient CFD Simulation on Dynamic Characteristics of Annular Seal under Large Eccentricities and Disturbances

: Annular seals of turbomachinery usually su ﬀ er from various degrees of eccentricities and disturbances due to the rotor–stator misalignment and radial loads, while the discussion of annular seal under both large static eccentricities and dynamic disturbances is relatively limited. In this paper, the applicability of linear assumption and reliability of nonlinear dynamic model for eccentric annular seals under large eccentricities and disturbances is discussed based on the investigation of seals with various rotor motions through computational ﬂuid dynamics (CFD). After the validation of transient CFD methods by comparison with experimental and bulk theory results, the dynamic behaviors of annular seal are analyzed by adopting both direct transient simulations and the nonlinear Muszynska model. The results show that the nonlinear dynamic model based on rotor circular whirls around seal center can predict the ﬂuid excitations of di ﬀ erent types of rotor motions well under small static eccentricities, while it is limited severely with large static eccentricities, which indicates that the dynamic characteristics of annular seal under large eccentricities are related with the rotor’s motion ways. The paper provides a reference for studies of rotor–seal system with complex rotor motions considering radial loads or running across the resonance region.


Introduction
Hydraulic machinery such as pumps and turbines is widely applied in various energy fields, playing a significant role in energy development, utilization and transformation. The vibration caused by the fluid forces generated in gap seals of hydraulic machinery tend to have important effects on the efficiency and vibration of rotor system [1]. Due to the rise of safety and efficiency concerns, dynamic characteristics of various annular seals have been studied by researchers [2][3][4][5]. Almost all of these studies are based on the assumption of small perturbation, hence linear dynamic characteristics of annular seals can be investigated. Generally, the annular seal is not the supporting element in design. Under the condition of static equilibrium, the rotor is normally concentric with the annular seal. Due to the axial-symmetry of seal geometry, as shown in Figure 1, the force coefficients of concentric seals show symmetric or skew symmetric features, as shown in Equation (1), where F x , F y are the X and Y components of fluid forces respectively; K and k denote direct and cross stiffness coefficients, Energies 2020, 13, x FOR PEER REVIEW 2 of 21 respectively; and M is direct mass coefficient. These five coefficients can be numerically computed by using the bulk flow model [6], CFD simulations by introducing moving reference frame [7,8] or transient method [9] and measured by perturbing the rotor or the stator [10].
However, under actual condition, the static eccentricity of the rotor may exist in annular seal due to the misalignment during assembly process or the effects of various side loads (e.g., impeller weight). The dynamic characteristics of eccentric annular seals, as shown in Figure 2, were also investigated based on the bulk flow model by Nelson and Nguyen [11]. The fluid force increments (ΔFx and ΔFy) induced by the small perturbation around static eccentricity position are similarly expressed in linearized rotordynamic form, as shown in Equation (2) [12].
where Δx and Δy define the rotor motion relative to the equilibrium position. Unlike concentric seals, the force coefficients of eccentric seals are no longer symmetric or skew symmetric due to rotor misalignment. This brings difficulties to the numerical solutions of force coefficients. Arghir and Frene [13] compared the bulk flow model of concentric seals and eccentric seals, the results showing that the terms of circumferential partial derivatives emerge in all bulk flow equations due to the static eccentricity of flow field. This can result in the coupling effect between circumferential momentum equation and continuity equation and make the solutions of both bulk flow equations and their perturbation equations very complex. As to the CFD method, the seal flow field disturbed by rotor circular whirl is not axisymmetric, as shown in Figure 2, and the steady-state simplified treatment by introducing moving reference frame is no longer applicable [8]. This means that transient simulations are necessary for evaluating force coefficients of eccentric seals.
To overcome numerical difficulties in eccentric seal research, Venkataraman and Palazzolo [14] determined the circumferential derivatives through a cubic spline interpolation method and However, under actual condition, the static eccentricity of the rotor may exist in annular seal due to the misalignment during assembly process or the effects of various side loads (e.g., impeller weight). The dynamic characteristics of eccentric annular seals, as shown in Figure 2, were also investigated based on the bulk flow model by Nelson and Nguyen [11]. The fluid force increments (∆F x and ∆F y ) induced by the small perturbation around static eccentricity position are similarly expressed in linearized rotordynamic form, as shown in Equation (2) [12].
− ∆F x ∆F y = k xx k xy k yx k yy ∆x ∆y + c xx c xy c yx c yy where ∆x and ∆y define the rotor motion relative to the equilibrium position. Unlike concentric seals, the force coefficients of eccentric seals are no longer symmetric or skew symmetric due to rotor misalignment. This brings difficulties to the numerical solutions of force coefficients. Arghir and Frene [13] compared the bulk flow model of concentric seals and eccentric seals, the results showing that the terms of circumferential partial derivatives emerge in all bulk flow equations due to the static eccentricity of flow field. This can result in the coupling effect between circumferential momentum equation and continuity equation and make the solutions of both bulk flow equations and their perturbation equations very complex. As to the CFD method, the seal flow field disturbed by rotor circular whirl is not axisymmetric, as shown in Figure 2, and the steady-state simplified treatment by introducing moving reference frame is no longer applicable [8]. This means that transient simulations are necessary for evaluating force coefficients of eccentric seals.
Energies 2020, 13, x FOR PEER REVIEW 2 of 21 respectively; and M is direct mass coefficient. These five coefficients can be numerically computed by using the bulk flow model [6], CFD simulations by introducing moving reference frame [7,8] or transient method [9] and measured by perturbing the rotor or the stator [10].
However, under actual condition, the static eccentricity of the rotor may exist in annular seal due to the misalignment during assembly process or the effects of various side loads (e.g., impeller weight). The dynamic characteristics of eccentric annular seals, as shown in Figure 2, were also investigated based on the bulk flow model by Nelson and Nguyen [11]. The fluid force increments (ΔFx and ΔFy) induced by the small perturbation around static eccentricity position are similarly expressed in linearized rotordynamic form, as shown in Equation (2) [12].
where Δx and Δy define the rotor motion relative to the equilibrium position. Unlike concentric seals, the force coefficients of eccentric seals are no longer symmetric or skew symmetric due to rotor  To overcome numerical difficulties in eccentric seal research, Venkataraman and Palazzolo [14] determined the circumferential derivatives through a cubic spline interpolation method and simplified the bulk flow equations of eccentric seals. Athavale and Hendricks [15] presented a small perturbation CFD method for calculation of rotordynamic coefficients of concentric and eccentric seals, and the SCISEAL code along with a modified SIMPLEC algorithm was adopted. Wu et al. [16] developed a new transient CFD method, which is based on rotor's variable-speed whirl; the results show that this new method can keep good accuracy of traditional transient method and save much computational time and cost in the meantime.
The research for fluid force presented above has mainly focused on linear fluid force analysis, and it was performed under the strict restriction and assumption that the whirl amplitude is relatively very small compared to the seal clearance (within 0.1 C r ; C r denotes the seal clearance). While large amplitude vibration often occurs during the passage of the critical speed of actual turbomachinery, the linear bulk flow analysis may not be applicable for the accurate fluid force characteristics in such situations with large amplitude. To describe the fluid forces of annular seal induced by large disturbances, the nonlinear dynamic model should be established. San Andres and Jeung [17] presented an orbit analysis method based on extended Reynolds equation to investigate force coefficients valid over a wide frequency range of a squeeze film damper bearing with large amplitude and static eccentricity. Ikemoto et al. [6] investigated the nonlinear fluid forces for the concentric seal with large whirl amplitude up to about a half of the clearance by using extended perturbation analysis of the bulk flow theory. Currently, the Muszynska's model proposed by Bently and Muszynska [18] is commonly used by researchers as a nonlinear dynamic model. Li and Chen [19] adopted the Muszynska's seal force model with the empirical parameters to investigate the 1:2 subharmonic resonance of labyrinth seal-rotor system. These empirical parameters obtained by employing the CFD analysis are used in the subsequent nonlinear analysis, regardless of whether the whirl amplitude is around the concentric position or not. He and Jing [20] indicated that Muszynska's model will not describe the dynamic characteristics of the rotor-seal system well when the rotor-seal system has larger eccentricity ratio. However, the present paper is devoted to develop nonlinear dynamic models of concentric seal with large whirl amplitude or eccentric seal with large static eccentricity and rather small whirl amplitude. The applicability of linear assumption and reliability of nonlinear model for seals under large static eccentricities and disturbance amplitude is rarely discussed in the literature. Thus, an investigation on the applicability of nonlinear Muszynska's model under large eccentricities and disturbances is wished for, particularly in nonlinear rotor-seal system research considering radial loads.
In experimental studies of eccentric seals, Marquette, Childs and Andres [21] measured the force coefficients of a plain liquid annular seal under different static eccentricities, and the results show that the force coefficients were more sensitive to the changes of static eccentricity than theoretically predicted. Childs, Arthur and Mehta [22] measured the net reaction forces of gas annular seals as the eccentricity ratios increased; negative stiffness created by unanticipated eccentricities may lead to over prediction of critical speeds, which illustrates the importance of concentric assembly of annular seals.
In this paper, three-dimensional (3D) transient CFD simulations based on dynamic mesh method are performed to evaluate the static and dynamic characteristics of eccentric annular seals. The obtained force coefficients and leakage rates are compared with Marquette's experiment [21] for validating the reliability of the transient CFD method. The effects of rotor disturbance amplitude on the dynamic characteristics of eccentric annular seals are analyzed to investigate the linear ranges of seal dynamic characteristics. In addition, transient CFD simulations and a nonlinear dynamic model are adopted to study the fluid excitations of annular seals induced by different rotor large motions. The nonlinear dynamic model is based on the famous Muszynska's model [18,23,24] and is obtained by fitting the "nominal" force coefficients of concentric annular seal under different whirl amplitude, as shown in eccentricities and disturbances are investigated in detail, which provides a solid basis for the research of seal-rotor system analysis by using Muszynska' model as nonlinear seal force.

Geometry Model and Grid
The plain annular seal adopted to perform the studies in this paper is applied in high speed hydrostatic journal bearings, which is tested in the apparatus and facility in Marquette's experiment. The work medium is water at 20°C. The geometric and operating parameters of the seal are listed in Table 1. As shown in Figure 3, the structured grids are generated in the concentric annular fluid domain by the CFD Preprocessor Gambit, which is geometry and mesh generation commercial software for computational fluid dynamics (CFD) analysis.

Geometry Model and Grid
The plain annular seal adopted to perform the studies in this paper is applied in high speed hydrostatic journal bearings, which is tested in the apparatus and facility in Marquette's experiment. The work medium is water at 20 ℃. The geometric and operating parameters of the seal are listed in Table 1. As shown in Figure 3, the structured grids are generated in the concentric annular fluid domain by the CFD Preprocessor Gambit, which is geometry and mesh generation commercial software for computational fluid dynamics (CFD) analysis.  The grid independence is checked by comparing the several grids with different radial grid densities. Under 80% eccentricity ratio, the radial and tangential components of fluid force are evaluated according to different grid models, as shown in Figure 4. The curves of "Fr refined" and "Ft refined" represent the radial and tangential fluid force of refined grid model, which has 36 radial layers with more than 10 layers near the both walls to keep y+ less than 5. The grid model of 16 radial layers is adopted considering the accuracy and computational time. With respect to the tangential and axial density, it can be seen in Table 2 that the results of fluid force show good convergence at Grid 3 (16 × 318 × 1448, i.e., there are 16 layers of grids generated along seal clearance in radial direction, 318 layers in axial direction and 1448 layers in circumferential direction) as the grid density The grid independence is checked by comparing the several grids with different radial grid densities. Under 80% eccentricity ratio, the radial and tangential components of fluid force are evaluated according to different grid models, as shown in Figure 4. The curves of "Fr refined" and "Ft refined" represent the radial and tangential fluid force of refined grid model, which has 36 radial layers with more than 10 layers near the both walls to keep y+ less than 5. The grid model of 16 radial layers is adopted considering the accuracy and computational time. With respect to the tangential and axial density, it can be seen in Table 2 that the results of fluid force show good convergence at Grid 3 (16 × 318 × 1448, i.e., there are 16 layers of grids generated along seal clearance in radial direction, 318 layers in axial direction and 1448 layers in circumferential direction) as the grid density changes to 1.25 or 1.5 times. This indicates that the present grid resolution (16 × 318 × 1448, 7,358,770 grid cells) is suitable for this research considering about the accuracy and efficiency of simulations.

3D Transient Solutions
Under various rotor disturbances, the static and dynamic characteristics of plain annular seal are investigated by simulating the transient flow in seal clearance. In this paper, the commercial CFD solver, ANSYS Fluent, is chosen to solve the 3D Reynolds-averaged Navier-Stokes equations. To achieve transient simulations, dynamic mesh problem should be firstly settled. As shown in Figure 2, the motion of rotor (i.e., rotating wall) can change the shape of fluid domain, and grids will change accordingly. However, due to high aspect ratio of grid cells in the clearance, the three types of dynamic methods in Fluent-spring-based smoothing, local remeshing and dynamic layering methods-tend to cause bad orthogonality or negative volume of grids.
To ensure good grid quality, the dynamic mesh model based on interpolation method [9,25] is adopted in this paper, which can effectively control the movement of the girds. First, nodes on rotating wall (i.e., rotor surface) are controlled to move according to the motion equation of the rotor and nodes on static wall keep stationary. Then, the ratio of nodes in the clearance is deduced according to the geometric relations of position of nodes in the clearance and movement of rotor. After that, the motions of grid nodes in the domain are determined by using the interpolation method based on the distances of the nodes from rotor and stator walls. Finally, the positions and velocities of grid nodes in the domain are obtained after the movement of rotor. Figure 5 shows the grid nodes moving in the clearance of annular seal. As illustrated in the figure, pf 0 (x 0 f, y 0 f) and pb 0 (x 0 b, y 0 b) represent the nodes of rotor surface and stator surface, respectively,

3D Transient Solutions
Under various rotor disturbances, the static and dynamic characteristics of plain annular seal are investigated by simulating the transient flow in seal clearance. In this paper, the commercial CFD solver, ANSYS Fluent, is chosen to solve the 3D Reynolds-averaged Navier-Stokes equations. To achieve transient simulations, dynamic mesh problem should be firstly settled. As shown in Figure 2, the motion of rotor (i.e., rotating wall) can change the shape of fluid domain, and grids will change accordingly. However, due to high aspect ratio of grid cells in the clearance, the three types of dynamic methods in Fluent-spring-based smoothing, local remeshing and dynamic layering methods-tend to cause bad orthogonality or negative volume of grids.
To ensure good grid quality, the dynamic mesh model based on interpolation method [9,25] is adopted in this paper, which can effectively control the movement of the girds. First, nodes on rotating wall (i.e., rotor surface) are controlled to move according to the motion equation of the rotor and nodes on static wall keep stationary. Then, the ratio of nodes in the clearance is deduced according to the geometric relations of position of nodes in the clearance and movement of rotor. After that, the motions of grid nodes in the domain are determined by using the interpolation method based on the distances of the nodes from rotor and stator walls. Finally, the positions and velocities of grid nodes in the domain are obtained after the movement of rotor. Figure 5 shows the grid nodes moving in the clearance of annular seal. As illustrated in the figure, pf 0 (x 0 f , y 0 f ) and pb 0 (x 0 b , y 0 b ) represent the nodes of rotor surface and stator surface, respectively, Energies 2020, 13, 4056 6 of 21 when the rotor is at the concentric position. pi 0 (x 0 i , y 0 i ) is an arbitrary node in the clearance domain of annular seal along the line between pf 0 and pb 0 . θ denotes the initial angular coordinate of node pf 0 . The superscript denotes the moving step of nodes and the subscript denotes the position of nodes.
represents the motion vector of moving rotor in Cartesian coordinates. di 1 denotes the vector from pi 0 to pi 1 (x 1 i , y 1 i ). Then, the new coordinates (x 1 f , y 1 f ) of pf 0 (current node pf 1 ) are defined as Equation (3).
Energies 2020, 13, x FOR PEER REVIEW 6 of 21 when the rotor is at the concentric position. pi 0 (x 0 i, y 0 i) is an arbitrary node in the clearance domain of annular seal along the line between pf 0 and pb 0 . θ denotes the initial angular coordinate of node pf 0 . The superscript denotes the moving step of nodes and the subscript denotes the position of nodes. d 1 (x 1 d, y 1 d) represents the motion vector of moving rotor in Cartesian coordinates. di 1 denotes the vector from pi 0 to pi 1 (x 1 i, y 1 i). Then, the new coordinates (x 1 f, y 1 f) of pf 0 (current node pf 1 ) are defined as Equation (3). The node of stator surface is assumed to stay still, the movement distance between rotor and stator is determined by the interpolation algorithm. Then, the new position coordinates (x 1 f, y 1 f) of pi 1 could be expressed as Equation (4).
where ra denotes the ratio of the distance between the nodes in the clearance domain and the static outer wall to the clearance. When the rotor is in concentric and eccentric position, the initial angular coordinate θ of pf 0 and the ratio of ra can be expressed by known parameters R and Cr and the coordinates of pf 0 , pi 0 and pb 0 according to collinear geometric relations of pf 1 , pi 1 and pb 1 . Then, the new position of pi 1 in the clearance after the movement of rotor can be obtained by substituting ra to Equation (3). The displacement of each node is restricted and calculated by mathematical procedures, which strictly ensures the movement coordination of adjacent grid nodes. The whole dynamic mesh process is implemented by adopting a subroutine linked with the CFD solver. This algorithm has been tested and the results show that, when the rotor whirled from the concentric position (with exaggerated seal clearance Cr), as shown in Figure 6a, to the eccentric position, as shown in Figure 6b, the grid distortion rate will increase but there is no negative volumes and highly distorted elements. The maximum grid aspect ratio will not exceed 200 even with eccentricity ratios (e/Cr, e denotes the rotor eccentricity) of 80%, which indicates that this dynamic mesh algorithm is suitable for the transient simulation with large eccentricity. The node of stator surface is assumed to stay still, the movement distance between rotor and stator is determined by the interpolation algorithm. Then, the new position coordinates (x 1 f , y 1 f ) of pi 1 could be expressed as Equation (4).
where ra denotes the ratio of the distance between the nodes in the clearance domain and the static outer wall to the clearance. When the rotor is in concentric and eccentric position, the initial angular coordinate θ of pf 0 and the ratio of ra can be expressed by known parameters R and C r and the coordinates of pf 0 , pi 0 and pb 0 according to collinear geometric relations of pf 1 , pi 1 and pb 1 . Then, the new position of pi 1 in the clearance after the movement of rotor can be obtained by substituting ra to Equation (3). The displacement of each node is restricted and calculated by mathematical procedures, which strictly ensures the movement coordination of adjacent grid nodes. The whole dynamic mesh process is implemented by adopting a subroutine linked with the CFD solver. This algorithm has been tested and the results show that, when the rotor whirled from the concentric position (with exaggerated seal clearance C r ), as shown in Figure 6a, to the eccentric position, as shown in Figure 6b, the grid distortion rate will increase but there is no negative volumes and highly distorted elements. The maximum grid aspect ratio will not exceed 200 even with eccentricity ratios (e/C r , e denotes the rotor eccentricity) of 80%, which indicates that this dynamic mesh algorithm is suitable for the transient simulation with large eccentricity.
Due to rotor eccentricity, one side of the grids is compressed and the maximum aspect ratio of the grids increases on basis of the initial grid model in Figure 6a. Considering the extreme thin grid layers, numerical computations are performed under double precision to ensure the stability and reliability of the result. The boundary conditions of 5.52 MPa total pressure and 0 Pa static pressure are, respectively, adopted at inlet and outlet. Both walls are set as no-slip walls and the rotating wall possesses a rotation speed which equals to r/min. The wall y+ of flow field under various disturbances is generally located in the range of 20-40 and the Realizable k-ε model with enhanced wall function is suitable to handle the Energies 2020, 13, 4056 7 of 21 situation [9,16]. The first-order implicit scheme is used for the discretization of time term. The chosen time step is equal to the wall rotation for 1 degree so that the courant number in most regions can be confined within 5 for stability. More than 360 steps are performed to ensure the stability of transient simulation according to different rotor motions. The second-order up-wind scheme with numerical under-relaxation is adopted to the convection term in the equations. The central-differencing scheme is employed to discretize the diffusion. The velocity-pressure coupling is solved by using the well-known SIMPLE strategy. Each simulation case for one revolution costs about 50 h on the platform of CPU is Intel ® Xeon ® Gold 6240 @ 2.60GHz with an average 16-parallel-processes solver. Due to rotor eccentricity, one side of the grids is compressed and the maximum aspect ratio of the grids increases on basis of the initial grid model in Figure 6a. Considering the extreme thin grid layers, numerical computations are performed under double precision to ensure the stability and reliability of the result. The boundary conditions of 5.52 MPa total pressure and 0 Pa static pressure are, respectively, adopted at inlet and outlet. Both walls are set as no-slip walls and the rotating wall possesses a rotation speed which equals to r/min. The wall y+ of flow field under various disturbances is generally located in the range of 20-40 and the Realizable k-ε model with enhanced wall function is suitable to handle the situation [9,16]. The first-order implicit scheme is used for the discretization of time term. The chosen time step is equal to the wall rotation for 1 degree so that the courant number in most regions can be confined within 5 for stability. More than 360 steps are performed to ensure the stability of transient simulation according to different rotor motions. The second-order up-wind scheme with numerical under-relaxation is adopted to the convection term in the equations. The central-differencing scheme is employed to discretize the diffusion. The velocity-pressure coupling is solved by using the well-known SIMPLE strategy. Each simulation case for one revolution costs about 50 h on the platform of CPU is Intel ® Xeon ® Gold 6240 @ 2.60GHz with an average 16-parallelprocesses solver.

Computing Static and Dynamic Characteristics of Eccentric Annular Seals
Transient CFD simulations are used to compute the static and dynamic characteristics of eccentric seals. The six cases with different static eccentricity ratios (se/Cr, where se denotes the static eccentricity of rotor) are investigated, respectively, 0%, 10%, 20%, 30%, 40% and 50%. The eccentric direction in +X direction is shown in Figure 2. The leakage rates of eccentric annular seals can be obtained by simulating the steady-state flow fields without rotor disturbances. The dynamic characteristics of eccentric seals can be analyzed by considering small rotor perturbations. The adopted perturbation is the circular whirl with a small whirl amplitude Δe (termed as dynamic eccentricity), as shown in Figure 2. The whirling speed Ω is constant. Given the suitability of small perturbation assumption, dynamic eccentricity ratio (Δe/Cr) should be very small (1% in the study). The small whirls are described by Equation (5). The fluid force increments (ΔFx and ΔFy) induced by perturbations are expressed by Equation (6).

Computing Static and Dynamic Characteristics of Eccentric Annular Seals
Transient CFD simulations are used to compute the static and dynamic characteristics of eccentric seals. The six cases with different static eccentricity ratios (se/C r , where se denotes the static eccentricity of rotor) are investigated, respectively, 0%, 10%, 20%, 30%, 40% and 50%. The eccentric direction in +X direction is shown in Figure 2. The leakage rates of eccentric annular seals can be obtained by simulating the steady-state flow fields without rotor disturbances. The dynamic characteristics of eccentric seals can be analyzed by considering small rotor perturbations. The adopted perturbation is the circular whirl with a small whirl amplitude ∆e (termed as dynamic eccentricity), as shown in Figure 2. The whirling speed Ω is constant. Given the suitability of small perturbation assumption, dynamic eccentricity ratio (∆e/C r ) should be very small (1% in the study). The small whirls are described by Equation (5). The fluid force increments (∆F x and ∆F y ) induced by perturbations are expressed by Equation (6).
where F x0 and F y0 represent fluid forces at equilibrium position. Substituting Equations (5) and (6) into Equation (2), F x and F y can be expressed as the harmonic functions of time, as shown in Equation (7): where By simulating the transient flow field with rotor perturbation, the time histories of F x and F y can be recorded by integrating the fluid pressure at each time step. Then, they are used to evaluate F x0 , Energies 2020, 13, 4056 8 of 21 F y0 and the four constant coefficients (A 1 , B 1 , A 2 and B 2 ) in Equation (7) by curve fittings. A 1 , B 1 , A 2 and B 2 are composed of a known whirling speed Ω and three unknown force coefficients, as shown in Equation (8). To obtain all the unknown force coefficients, A 1 , B 1 , A 2 and B 2 under at least three (generally five is desired considering the fitting error) whirling speeds should be determined. Hence, at least three transient CFD simulations should be performed.

Fitting Nonlinear Dynamic Model
Force coefficients of annular seals can only be used to describe fluid forces induced by rotor small perturbations. The Muszynska's model is adopted to describe the fluid forces of annular seal induced by large disturbances. It is derived based on a serial of experiments and adopts nonlinear dynamic parameters similar with force coefficients to associate fluid forces with rotor motion, as shown in Equation (9): for ε = x 2 + y 2 /C r n, τ 0 and b are empirical factors for certain seal structure; S 0 , D 0 and m f can be computed using Black-Childs formulas [26]. When the seal is under steady working condition (constant ω and ∆P), the S 0 , D 0 , m f and empirical factors in Muszynska's model become constant values. Namely, nonlinear dynamic parameters in matrices are only related with eccentricity ratio ε. Thus, Equation (9) can be expressed in a simplified form, as shown in Equation (10): The Muszynska's model under constant working condition is very similar to the linear dynamic model of concentric annular seal in Equation (1). The only difference is that dynamic parameters in Equation (10) are nonlinear functions of ε and can be used to describe fluid forces induced by rotor large disturbances. The expressions of nonlinear dynamic parameters can be determined based on the formulas in Equation (9), but proper empirical factors need to be chosen. In addition, the nonlinear expressions can be fitted based on the "nominal" force coefficients of concentric seal under different eccentricities (as shown in Figure 2). These "nominal" force coefficients can be computed by using CFD methods [8,27]. Rotor perturbation is the circular whirl around seal center. Usually, whirl amplitude (i.e., rotor eccentricity) is controlled within 0.1C r for satisfying the linear assumption. To obtain "nominal" force coefficients under different eccentricities, the limitation is broken here, and the adopted whirl amplitudes are located in the range of 0.01C r -0.8C r (i.e., ε in 1%-80%).
Transient CFD simulations are conducted to solve the flow field disturbed by constant-speed circular whirls, and fluid-induced forces can be obtained. Based on these fluid forces, the "nominal" force coefficients of concentric annular seal can be evaluated [8] and used to generate the nonlinear dynamic model. With respect to the Muszynska's model, the new nonlinear model does not need any empirical factors. Theoretically, it can describe fluid forces (seal forces) induced by various rotor motions within eccentricity ratio 80%. Its reliability and suitability are discussed in Section 3.

Dynamic Characteristics of Different Static Eccentric Seals and Comparisons
The

Dynamic Characteristics of Different Static Eccentric Seals and Comparisons
The leakage rates and force coefficients of annular seal under different static eccentricity positions and same whirl amplitude ratio (Δe/Cr = 1%) are computed using the numerical scheme based on transient CFD simulations (see Section 2.3). They are, respectively, shown in              In Figures 7 and 8, the measured values of direct and cross stiffness coefficients are larger than numerical values from transient CFD simulations. The test apparatus of Marquette's experiment does not include any device to control or measure the inlet swirl, and the non-uniformity of incoming flow is not considered during the testing procedure, which may lead to variations of stiffness coefficients with eccentric directions according to Wu's research [28]. The boundary condition of seal model illustrated in Marquette's research is only pressure condition for inlet and outlet. There is no geometry information or measurement of inlet, which is also mentioned by the authors as a drawback. The reason of the lower accuracy is that it is much difficult to ensure the boundary condition of CFD analysis, especially for the inlet, consistent with the Marquette's experiment. Tae Woong Ha's study [27] shows the similar difference between the stiffness results of CFD analysis and the experimental results. Despite the differences of values in Figure 7, it shows a high level of consistency between results of transient simulations and measured values as static eccentricity ratio increases.
In Figures 9 and 10, direct damping coefficients from transient simulations and bulk flow method are both a little higher than measured values, and CFD results are closer to experimental values. Cross damping coefficients from the three approaches are all close in size. As shown in Figure 11, although numerical results of mass coefficients do not coincide with measured values in variation trend, they are close in size. Figure 12 shows the variations of seal leakage rate with the eccentricity. As rotor eccentricity increases, seal leakage rate just slightly rises. In Figure 12, leakage rates computed by transient simulations are very close to measured values, which further indicate the reliability of CFD method. This also shows that leakage rates of annular seals mainly depend on sealed differential pressures and are not sensitive to flow status at seal inlet, unlike force coefficients [28].
On the whole, rotor eccentricities change the static and dynamic characteristics of annular seals to some extent; the behaviors of annular seals should be specially considered when the rotor is far away from seal center. By comparing with experimental data, transient CFD simulations are effective in computing the static and dynamic behaviors of eccentric seals.

Effects of Disturbance Amplitude on Force Coefficients of Eccentric Annular Seal
The force coefficients of eccentric seals in Figures 7-12 are computed using a very small dynamic eccentricity (eccentricity ratio 1%). To study the effects of disturbance amplitude, two more dynamic eccentricity ratios, 5% and 10%, are separately adopted to evaluate force coefficients of eccentric seals. The results are presented in Figures 13-17  concentric annular seal are still linear when dynamic eccentricity ratio reaches 10%, which has been widely recognized by researchers [7,8,29]. However, as to the eccentric annular seal, its force coefficients tend to be sensitive to dynamic eccentricities (i.e., disturbance amplitude). As rotor static eccentricity increases, the sensitivity to rotor disturbances strengthens gradually and the linear range of seal dynamic characteristics narrows. With respect to the concentric seal, the annular seal under large eccentricity is more likely to show nonlinear characteristics.       As shown in Figures 13-17, force coefficients based on three dynamic eccentricities are not the same. Their differences grow gradually as rotor static eccentricity increases, especially the coefficients k xx , k yy , k xy , c xx , m xx and m yy . For the concentric seal with static eccentricity zero, its force coefficients are generally unchanged with varying dynamic eccentricities. Namely, the dynamic characteristics of concentric annular seal are still linear when dynamic eccentricity ratio reaches 10%, which has been widely recognized by researchers [7,8,29]. However, as to the eccentric annular seal, its force coefficients tend to be sensitive to dynamic eccentricities (i.e., disturbance amplitude). As rotor static eccentricity increases, the sensitivity to rotor disturbances strengthens gradually and the linear range of seal dynamic characteristics narrows. With respect to the concentric seal, the annular seal under large eccentricity is more likely to show nonlinear characteristics.

Nonlinear Dynamic Model
As discussed above, force coefficients of eccentric seals can only be used to describe fluid forces induced by very small perturbations due to their strong sensitivity to disturbance amplitude. To express seal forces under large eccentricities or large disturbances, the nonlinear dynamic model is determined based on the thought in Section 2.3. The "nominal" force coefficients of concentric seal are computed based on different whirl amplitudes. They are presented in Figure 18 along with polynomial fitting curves. Piecewise fittings are used for direct stiffness K and cross damping c. As shown in Figure 18, the fitting effects of five force coefficients are satisfactory. The nonlinear expressions of these coefficients are listed as follows: M(ε) = −5.28ε 4 + 6.63ε 3 − 3.58ε 2 + 0.36ε + 3.50, ε ≤ 0.8, R 2 = 0.9905 (13) k(ε) ×10 7 = 14.5ε 5 − 24.2ε 4 + 15.5ε 3 − 3.22ε 2 + 0.250ε + 0.517, ε ≤ 0.8 R 2 = 0.9999 (14) C(ε) ×10 4 = 10.7ε 4 − 11.6ε 3 + 5.57ε 2 − 0.682ε + 3.01, ε ≤ 0.8, R 2 = 0.9994 (15) Energies 2020, 13 Substituting these nonlinear expressions into Equation (10), the nonlinear dynamic model is obtained. It will be used to evaluate fluid forces induced by rotor large disturbances along with transient CFD simulations. Substituting these nonlinear expressions into Equation (10), the nonlinear dynamic model is obtained. It will be used to evaluate fluid forces induced by rotor large disturbances along with transient CFD simulations.

Fluid Excitations under Large Disturbances
The fluid forces of annular seal under rotor large motions are computed by nonlinear dynamic model as well as transient CFD simulations. The suitability of nonlinear dynamic model for various rotor disturbances is investigated through comparisons with direct transient simulations. Several typical rotor motions are chosen for the investigation.

Constant-Speed Circular Whirl Around Seal Center
The rotor performs circular whirl around seal center with speed 10,200 r/min ( Figure 1) and the whirl magnitude is 55% C r . The motion equation is shown as below. Substituting Equations (11)- (16) into Equation (10), the induced fluid forces (F x and F y ) are obtained by the nonlinear dynamic model. They are shown in Figure 19 along with the results from transient CFD simulations.
Because the prescribed initial solution is not absolutely accurate, the transient CFD computation needs passing a period of time to eliminate its effects. Fluid forces computed at initial some time steps are not true and can be ignored. In Figure 19, fluid force curves from nonlinear dynamic model and transient CFD simulations are in good agreement. Maximum differences are only 16.5 N for both F x and F y , and they are mainly caused by fitting and computing errors. This indicates that the present nonlinear dynamic model is adequate to describe seal fluid excitations induced by circular whirls around seal center.
Because the prescribed initial solution is not absolutely accurate, the transient CFD computation needs passing a period of time to eliminate its effects. Fluid forces computed at initial some time steps are not true and can be ignored. In Figure 19, fluid force curves from nonlinear dynamic model and transient CFD simulations are in good agreement. Maximum differences are only 16.5 N for both Fx and Fy, and they are mainly caused by fitting and computing errors. This indicates that the present nonlinear dynamic model is adequate to describe seal fluid excitations induced by circular whirls around seal center.

Constant-Speed Circular Whirl Around Static Position
The rotor is assumed to perform circular whirl around static eccentricity position, as shown in Figure 2. The dynamic eccentricity ratio is 10%. Two whirl centers correspond to static eccentricity ratio 20% and 50%, respectively. The whirling speed is same with the rotating speed, 10,200 r/min. Rotor movements (x, y) can be expressed by Equation (3). Substituting Equations (3) and (11)-(15) into Equation (10)

Constant-Speed Circular Whirl Around Static Position
The rotor is assumed to perform circular whirl around static eccentricity position, as shown in Figure 2. The dynamic eccentricity ratio is 10%. Two whirl centers correspond to static eccentricity ratio 20% and 50%, respectively. The whirling speed is same with the rotating speed, 10,200 r/min. Rotor movements (x, y) can be expressed by Equation (3). Substituting Equations (3) and (11)  As shown in Figure 20, for the circular whirl around the static position with eccentricity ratio 20%, fluid forces from nonlinear dynamic model and transient CFD simulations agree well.
Maximum differences are 9.42 N for Fx and 4.65 N for Fy. However, when the static eccentricity ratio increases to 50%, fluid force curves obtained by these two approaches are no longer consistent, as shown in Figure 21. The maximum differences are 111 N for Fx and 36.7 N for Fy. This indicates that the nonlinear dynamic model has low reliability for rotor disturbances under large eccentricity. However, the circular whirl around concentric position is an exception (see Figure 19). The nonlinear dynamic model can deal with it well. The two rotor motions corresponding to Figures 19 and 21 are  As shown in Figure 20, for the circular whirl around the static position with eccentricity ratio 20%, fluid forces from nonlinear dynamic model and transient CFD simulations agree well.
Maximum differences are 9.42 N for Fx and 4.65 N for Fy. However, when the static eccentricity ratio increases to 50%, fluid force curves obtained by these two approaches are no longer consistent, as shown in Figure 21. The maximum differences are 111 N for Fx and 36.7 N for Fy. This indicates that the nonlinear dynamic model has low reliability for rotor disturbances under large eccentricity. However, the circular whirl around concentric position is an exception (see Figure 19). The nonlinear dynamic model can deal with it well. The two rotor motions corresponding to Figures 19 and 21 are As shown in Figure 20, for the circular whirl around the static position with eccentricity ratio 20%, fluid forces from nonlinear dynamic model and transient CFD simulations agree well.
Maximum differences are 9.42 N for F x and 4.65 N for F y . However, when the static eccentricity ratio increases to 50%, fluid force curves obtained by these two approaches are no longer consistent, as shown in Figure 21. The maximum differences are 111 N for F x and 36.7 N for F y . This indicates that the nonlinear dynamic model has low reliability for rotor disturbances under large eccentricity. However, the circular whirl around concentric position is an exception (see Figure 19). The nonlinear dynamic model can deal with it well. The two rotor motions corresponding to Figures 19 and 21 are both under large eccentricities. The only difference is their motion ways. The rotor motion corresponding to Figure 19 is the whirl around seal center, and "nominal" force coefficients used for generating the nonlinear dynamic model are based on this motion way. Therefore, it is understandable that the nonlinear model applies well. From this point of view, the reason that nonlinear dynamic model does not suitable for circular whirls around large eccentricity position (as shown in Figure 21) can be assumed to be that the force coefficients (or dynamic characteristics) of annular seal under large eccentricity are related to rotor motion ways.

1D Harmonic Shaking Motions
To validate the assumption proposed above, the circular whirl around the static position with eccentricity ratio 50% is divided into two separate 1D shaking motions, as shown in Figure 22. One is the shaking motion in X direction; the other is in Y direction. The Y direction is also the tangential direction of the concentric whirl. Namely, the Y-directional shaking is somewhat similar to the circular whirl around seal center. The expressions of two shaking motions are, respectively, presented in Equations (17) and (18). The amplitude A is 10% C r and the harmonic frequency Ω corresponds to the speed 10,200 r/min.
x = se + Acos(Ωt) y = 0 (17) Energies 2020, 13, x FOR PEER REVIEW 17 of 21 in Equations (17) and (18). The amplitude A is 10% Cr and the harmonic frequency Ω corresponds to the speed 10,200 r/min.  Figure 24. This can be attributed to the slight likeness of Y-directional shaking with the circular whirl around seal center. In Figure 20, maximum differences of two approaches are 18.5 N for Fx and 17.7 N for Fy, and they are much smaller than those in Figure 23. Namely, the nonlinear dynamic model is more applicable to rotor disturbances similar to the circular whirl around concentric position. The assumption proposed in Section 4.2 is confirmed and can explain the low reliability of nonlinear dynamic model. Under large The comparison with transient CFD simulations shows that the reliability of nonlinear dynamic model is low in predicting fluid forces induced by X-directional shaking. Maximum differences of two approaches are 119 N for F x and 30.8 N for F y . However, as to the Y-directional shaking, the reliability of the nonlinear model is obviously improved, as shown in Figure 24. This can be attributed to the slight likeness of Y-directional shaking with the circular whirl around seal center. In Figure 20, maximum differences of two approaches are 18.5 N for F x and 17.7 N for F y , and they are much smaller than those in Figure 23. Namely, the nonlinear dynamic model is more applicable to rotor disturbances similar to the circular whirl around concentric position. The assumption proposed in Section 4.2 is confirmed and can explain the low reliability of nonlinear dynamic model. Under large eccentricity, the dynamic characteristics of seals are varied with the motion ways of the rotor, and the nonlinear dynamic model based on a specific motion way is incompetent in dealing with all kinds of rotor disturbances. two approaches are 119 N for Fx and 30.8 N for Fy. However, as to the Y-directional shaking, the reliability of the nonlinear model is obviously improved, as shown in Figure 24. This can be attributed to the slight likeness of Y-directional shaking with the circular whirl around seal center. In Figure 20, maximum differences of two approaches are 18.5 N for Fx and 17.7 N for Fy, and they are much smaller than those in Figure 23. Namely, the nonlinear dynamic model is more applicable to rotor disturbances similar to the circular whirl around concentric position. The assumption proposed in Section 4.2 is confirmed and can explain the low reliability of nonlinear dynamic model. Under large eccentricity, the dynamic characteristics of seals are varied with the motion ways of the rotor, and the nonlinear dynamic model based on a specific motion way is incompetent in dealing with all kinds of rotor disturbances.

Quasi-Circular (Spiral) Whirl
In this section, the rotor is assumed to perform the circular whirl around seal center with growing whirl radius, i.e., the spiral whirl, in order to validate the suitability of nonlinear dynamic model for quasi-circular whirls no matter eccentricity magnitudes. The whirling speed is r/min and the rotor eccentricity ratio rises linearly to 60% within six whirl periods. The whirl equation is presented in Equation (17) and the whirl orbit is shown in Figure 25.
where T is the whirl period and f = 60% Cr/(6 T), indicating the eccentricity speed of the rotor. In actual applications, the spiral whirl in Figure 25 represents the destabilizing process of the rotor. Fluid forces induced by the destabilizing whirl are obtained, respectively, by nonlinear dynamic model and transient CFD simulations. They are shown in Figure 26. The rise of rotor eccentricity with time leads to the increasing fluid force (i.e., the resultant force of Fx and Fy). The Fx and Fy from nonlinear model are in good agreement with those from transient CFD simulations. Maximum differences of two approaches are 42.7 N (i.e., relative error 2.3%) for Fx and 35.7 N (relative error 2.3%) for Fy. Namely, nonlinear dynamic model is reliable in evaluating fluid forces induced by the quasi-circular whirl around seal center without special limitations on eccentricity magnitudes.

Quasi-Circular (Spiral) Whirl
In this section, the rotor is assumed to perform the circular whirl around seal center with growing whirl radius, i.e., the spiral whirl, in order to validate the suitability of nonlinear dynamic model for quasi-circular whirls no matter eccentricity magnitudes. The whirling speed is r/min and the rotor eccentricity ratio rises linearly to 60% within six whirl periods. The whirl equation is presented in Equation (17) and the whirl orbit is shown in Figure 25.
where T is the whirl period and f = 60% C r /(6 T), indicating the eccentricity speed of the rotor. In actual applications, the spiral whirl in Figure 25 represents the destabilizing process of the rotor. Fluid forces induced by the destabilizing whirl are obtained, respectively, by nonlinear dynamic model and transient CFD simulations. They are shown in Figure 26. The rise of rotor eccentricity with time leads to the increasing fluid force (i.e., the resultant force of F x and F y ). The F x and F y from nonlinear model are in good agreement with those from transient CFD simulations. Maximum differences of two approaches are 42.7 N (i.e., relative error 2.3%) for F x and 35.7 N (relative error 2.3%) for F y . Namely, nonlinear dynamic model is reliable in evaluating fluid forces induced by the quasi-circular whirl around seal center without special limitations on eccentricity magnitudes.
rotor. Fluid forces induced by the destabilizing whirl are obtained, respectively, by nonlinear dynamic model and transient CFD simulations. They are shown in Figure 26. The rise of rotor eccentricity with time leads to the increasing fluid force (i.e., the resultant force of Fx and Fy). The Fx and Fy from nonlinear model are in good agreement with those from transient CFD simulations. Maximum differences of two approaches are 42.7 N (i.e., relative error 2.3%) for Fx and 35.7 N (relative error 2.3%) for Fy. Namely, nonlinear dynamic model is reliable in evaluating fluid forces induced by the quasi-circular whirl around seal center without special limitations on eccentricity magnitudes.

Conclusions
In the paper, dynamic characteristics of the annular plain liquid seal under various large rotor disturbance motions are studied using the transient CFD method based on dynamic mesh technique and nonlinear Muszynska' model.
Force coefficients and leakage rates of annular seal under different static eccentricities are evaluated. The reliability of transient CFD simulation is validated by comparing the force coefficients and leakage rates with those from the Marquette's experiment and bulk flow method. With increasing static eccentricity, these force coefficients show clearly asymmetric behavior and obvious changes. The force coefficients from transient CFD simulations show a high consistency with experimental values despite the different values of stiffness. The error sources are mainly form the influence of upstream and inlet boundary condition due to the drawback of the experimental apparatus for absent inlet control. Leakage rates computed by the CFD method fit better to measured values than those from the bulk flow method, which indicates that leakage rates are insensitive to static eccentricity.
As to the concentric annular seal, its dynamic characteristics are usually supposed to be linear (namely, constant force coefficients) when the rotor disturbance is within 10% Cr. However, this conclusion is not suitable for the eccentric annular seal, especially the seal under large static eccentricity. As rotor static eccentricity increases, the force characteristics of annular seal become more sensitive to whirl amplitude, in other words, the linear range of dynamic characteristics narrows gradually. With respect to the concentric seal, the annular seal with large eccentricity is easier to show nonlinear characteristics.
According to the Muszynska's model, a nonlinear dynamic model is presented in the paper for describing nonlinear seal forces induced by rotor large disturbances. The suitability of the nonlinear model for all kinds of rotor disturbances is studied through four forms of rotor motions. The nonlinear dynamic model is suitable for various rotor disturbances when the rotor is under small static eccentricity (e.g., eccentricity ratio under 20%). However, when rotor static eccentricity is large (e.g., eccentricity ratio 50%), the nonlinear dynamic model based on circular whirls around eccentric center becomes incompetent and unsatisfactory. It shows high reliability only for circular or quasicircular whirls around concentric center. This means that dynamic characteristics of annular seal under large disturbance are related to rotor motion ways. For the annular seals under large dynamic eccentricity (whirl amplitude) and rather small static eccentricity (e.g., static eccentricity ratio under 20% in this case), the nonlinear Muszynska's model performs well when dealing with large rotor

Conclusions
In the paper, dynamic characteristics of the annular plain liquid seal under various large rotor disturbance motions are studied using the transient CFD method based on dynamic mesh technique and nonlinear Muszynska' model.
Force coefficients and leakage rates of annular seal under different static eccentricities are evaluated. The reliability of transient CFD simulation is validated by comparing the force coefficients and leakage rates with those from the Marquette's experiment and bulk flow method. With increasing static eccentricity, these force coefficients show clearly asymmetric behavior and obvious changes. The force coefficients from transient CFD simulations show a high consistency with experimental values despite the different values of stiffness. The error sources are mainly form the influence of upstream and inlet boundary condition due to the drawback of the experimental apparatus for absent inlet control. Leakage rates computed by the CFD method fit better to measured values than those from the bulk flow method, which indicates that leakage rates are insensitive to static eccentricity.
As to the concentric annular seal, its dynamic characteristics are usually supposed to be linear (namely, constant force coefficients) when the rotor disturbance is within 10% C r . However, this conclusion is not suitable for the eccentric annular seal, especially the seal under large static eccentricity. As rotor static eccentricity increases, the force characteristics of annular seal become more sensitive to whirl amplitude, in other words, the linear range of dynamic characteristics narrows gradually. With respect to the concentric seal, the annular seal with large eccentricity is easier to show nonlinear characteristics.
According to the Muszynska's model, a nonlinear dynamic model is presented in the paper for describing nonlinear seal forces induced by rotor large disturbances. The suitability of the nonlinear model for all kinds of rotor disturbances is studied through four forms of rotor motions. The nonlinear dynamic model is suitable for various rotor disturbances when the rotor is under small static eccentricity (e.g., eccentricity ratio under 20%). However, when rotor static eccentricity is large (e.g., eccentricity ratio 50%), the nonlinear dynamic model based on circular whirls around eccentric center becomes incompetent and unsatisfactory. It shows high reliability only for circular or quasi-circular whirls around concentric center. This means that dynamic characteristics of annular seal under large disturbance are related to rotor motion ways. For the annular seals under large dynamic eccentricity (whirl amplitude) and rather small static eccentricity (e.g., static eccentricity ratio under 20% in this case), the nonlinear Muszynska's model performs well when dealing with large rotor disturbances. The range of capability of this nonlinear model depends on the typical parameters of annular seals. It can also explain why Muszynska's model is out of action when rotor-seal system has a large eccentricity ratio in He's research.
On the whole, dynamic characteristics of annular seals under large disturbance are very complex. They are very sensitive to various rotor motion ways including whirl amplitude and static eccentricity. For the seal with large disturbances motion of a small static eccentricity, the nonlinear Muszynska's model performs reliably, which provides a solid basis for the seal-rotor system analysis using nonlinear seal force model. The capability and limitation of nonlinear dynamic model under large disturbances needs further investigation.  Rotor eccentricity or whirl radius (mm] f Eccentricity speed of the rotor F x , F y Fluid forces in X and Y directions (N) F x0 , F y0 Fluid forces at equilibrium position (N) ∆F x , ∆F y The increments of fluid forces relative to F x0 and F y0 (N) K, k Direct and cross stiffness coefficients of concentric seal (N/m) k xx , k yy , k xy , k yx Stiffness coefficients of eccentric annular seal (N/m) M Direct mass coefficient of concentric annular seal (kg) m xx , m yy Direct mass coefficients of eccentric annular seal (kg) ra Ratio of the distance between the nodes in the clearance domain and the outer static wall to the clearance se Rotor static eccentricity (mm) se/C r Static eccentricity ratio t Time (s) x, y The displacements of rotor center (mm) ∆x, ∆y Rotor displacements relative to equilibrium position (mm) ∆e Rotor dynamic eccentricity (mm) ε Eccentricity ratio (e/C r ) ω Rotating speed of the rotor (rpm) Ω Whirling speed of the rotor (rpm)