*2.1. 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.


**Table 1.** Parameters of plain annular seal.

**Figure 3.** Numerical model of concentric annular seal.

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.

**Figure 4.** Radial grid density study.



<sup>1</sup> Note: by comparing with Grid 3 (16 × 318 × 1448, radial × axial × tangential layers).

#### *2.2. 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* <sup>0</sup> (*x*<sup>0</sup> *<sup>f</sup>*, *y*<sup>0</sup> *<sup>f</sup>*) and *pb*<sup>0</sup> (*x*<sup>0</sup> *<sup>b</sup>*, *y*<sup>0</sup> *<sup>b</sup>*) represent the nodes of rotor surface and stator surface, respectively, when the rotor is at the concentric position. *pi*<sup>0</sup> (*x*<sup>0</sup> *i, y*<sup>0</sup> *<sup>i</sup>*) is an arbitrary node in the clearance domain of annular seal along the line between *pf* <sup>0</sup> 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*<sup>1</sup> *<sup>d</sup>*, *y*<sup>1</sup> *<sup>d</sup>*) represents the motion vector of moving rotor in Cartesian coordinates. *di*<sup>1</sup> denotes the vector from *pi*<sup>0</sup> to *pi*<sup>1</sup> (*x*<sup>1</sup> *<sup>i</sup>*, *y*<sup>1</sup> *i*). Then, the new coordinates (*x*<sup>1</sup> *<sup>f</sup>*, *y*<sup>1</sup> *<sup>f</sup>*) of *pf* <sup>0</sup> (current node *pf* 1) are defined as Equation (3).

$$\mathbf{x}\_f^1 = \mathbf{x}\_f^0 + \mathbf{x}\_{d'}^1 \ y\_f^1 = y\_f^0 + y\_{d'}^1 \tag{3}$$

**Figure 5.** Schematic diagram of moving grid node.

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*<sup>1</sup> *f, y*<sup>1</sup> *<sup>f</sup>*) of *pi*<sup>1</sup> could be expressed as Equation (4).

$$\mathbf{x}\_i^1 = \mathbf{x}\_i^0 + ra \times \mathbf{x}\_{d'}^1 \ y\_i^1 = y\_i^0 + ra \times y\_{d'}^1 \tag{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* <sup>0</sup> and the ratio of *ra* can be expressed by known parameters *R* and *Cr* and the coordinates of *pf* 0, *pi*<sup>0</sup> and *pb*<sup>0</sup> according to collinear geometric relations of *pf* 1, *pi*<sup>1</sup> and *pb*1. Then, the new position of *pi*<sup>1</sup> 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.

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-parallel-processes solver.

**Figure 6.** Diagram of cross section of the meshed rotor: (**a**) initial grids; and (**b**) moved grids.

### *2.3. 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).

$$\begin{cases} \Delta \mathbf{x} = \mathbf{x} - \mathbf{s}\mathbf{e} = \Delta \mathbf{c} \cos(\Omega t) \\ \Delta y = y - 0 = \Delta \mathbf{c} \sin(\Omega t) \end{cases} \tag{5}$$

$$\begin{cases} \Delta F\_{\chi} = F\_{\chi} - F\_{x0} \\ \Delta F\_{\mathcal{Y}} = F\_{\mathcal{Y}} - F\_{\mathcal{Y}0} \end{cases} \tag{6}$$

where *Fx0* and *Fy0* represent fluid forces at equilibrium position. Substituting Equations (5) and (6) into Equation (2), *Fx* and *Fy* can be expressed as the harmonic functions of time, as shown in Equation (7):

$$\begin{cases} F\_x = F\_{x0} + A\_1 \Delta \text{accos}(\Omega t) + B\_1 \Delta \text{csim}(\Omega t) \\\ F\_y = F\_{y0} + B\_2 \Delta \text{accos}(\Omega t) + A\_2 \Delta \text{csim}(\Omega t) \end{cases} \tag{7}$$

where

$$\begin{cases} \begin{aligned} A\_1 &= -k\_{xx} - c\_{xy}\Omega + m\_{xx}\Omega^2 \\ B\_1 &= -k\_{xy} + c\_{xx}\Omega + m\_{xy}\Omega^2 \\ A\_2 &= -k\_{yy} + c\_{yx}\Omega + m\_{yy}\Omega^2 \\ B\_2 &= -k\_{yx} - c\_{yy}\Omega + m\_{yx}\Omega^2 \end{aligned} \end{cases} \tag{8}$$

By simulating the transient flow field with rotor perturbation, the time histories of *Fx* and *Fy* can be recorded by integrating the fluid pressure at each time step. Then, they are used to evaluate *Fx*0, *Fy*<sup>0</sup> and the four constant coefficients (*A1*, *B1*, *A2* and *B2*) in Equation (7) by curve fittings. *A1*, *B1*, *A2* and *B2* are composed of a known whirling speed Ω and three unknown force coefficients, as shown in Equation (8). To obtain all the unknown force coefficients, *A1*, *B1*, *A2* and *B2* 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.
