*2.4. 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):

$$
\begin{Bmatrix} F\_x \\ F\_y \end{Bmatrix} = -\begin{bmatrix} S - m\_f \tau\_1^2 \omega^2 & \tau\_1 \omega D \\ -\tau\_1 \omega D & S - m\_f \tau\_1^2 \omega^2 \end{bmatrix} \begin{Bmatrix} \mathbf{x} \\ \mathbf{y} \end{Bmatrix} - \begin{bmatrix} D & 2\tau\_1 \omega m\_f \\ -2\tau\_1 \omega m\_f & D \end{bmatrix} \begin{Bmatrix} \dot{\mathbf{x}} \\ \dot{\mathbf{y}} \end{Bmatrix} - \begin{bmatrix} m\_f & 0 \\ 0 & m\_f \end{bmatrix} \begin{Bmatrix} \ddot{\mathbf{x}} \\ \ddot{\mathbf{y}} \end{Bmatrix}, \tag{9}
$$

where *S* = *S*<sup>0</sup> <sup>1</sup> <sup>−</sup> <sup>ε</sup><sup>2</sup> <sup>−</sup>*<sup>n</sup>* , τ<sup>1</sup> = τ0(1 − ε) *b* , *D* = *D*<sup>0</sup> <sup>1</sup> <sup>−</sup> <sup>ε</sup><sup>2</sup> <sup>−</sup>*<sup>n</sup>* , for ε =  *x*<sup>2</sup> + *y*2/*Cr*

*n*, τ<sup>0</sup> and *b* are empirical factors for certain seal structure; *S*0, *D*<sup>0</sup> and *mf* can be computed using Black–Childs formulas [26]. When the seal is under steady working condition (constant ω and Δ*P*), the *S0*, *D0*, *mf* 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):

$$
\begin{Bmatrix} \mathbf{F\_x} \\ \mathbf{F\_y} \end{Bmatrix} = - \begin{bmatrix} \mathbf{K}(\varepsilon) & \mathbf{k}(\varepsilon) \\ -\mathbf{k}(\varepsilon) & \mathbf{K}(\varepsilon) \end{bmatrix} \begin{Bmatrix} \mathbf{x} \\ \mathbf{y} \end{Bmatrix} - \begin{bmatrix} \mathbf{C}(\varepsilon) & \mathbf{c}(\varepsilon) \\ -\mathbf{c}(\varepsilon) & \mathbf{C}(\varepsilon) \end{bmatrix} \begin{Bmatrix} \dot{\mathbf{x}} \\ \dot{\mathbf{y}} \end{Bmatrix} - \begin{bmatrix} \mathbf{M}(\varepsilon) & \mathbf{0} \\ \mathbf{0} & \mathbf{M}(\varepsilon) \end{bmatrix} \begin{Bmatrix} \ddot{\mathbf{x}} \\ \ddot{\mathbf{y}} \end{Bmatrix} \tag{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.1*Cr* 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.01*Cr*–0.8*Cr* (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.

#### **3. Results and Discussions**

#### *3.1. Dynamic Characteristics of Di*ff*erent 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 Figures 7–12 along with the results from San Andres' bulk flow method and Marquette's experiments.

**Figure 7.** Direct Stiffness coefficients of eccentric annular seals.

**Figure 8.** Cross Stiffness coefficients of eccentric annular seals.

**Figure 9.** Direct damping coefficients of eccentric annular seals.

**Figure 10.** Cross damping coefficients of eccentric annular seals.

**Figure 11.** Mass coefficients of eccentric annular seals.

**Figure 12.** Leakage rates of eccentric annular seals.

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.

#### *3.2. E*ff*ects of Disturbance Amplitude on Force Coe*ffi*cients 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 for comparisons.

**Figure 13.** Direct stiffness coefficients under different dynamic eccentricities.

**Figure 14.** Cross stiffness coefficients under different dynamic eccentricities.

**Figure 15.** Direct damping coefficients under different dynamic eccentricities.

**Figure 16.** Cross damping coefficients under different dynamic eccentricities.

**Figure 17.** Mass coefficients under different dynamic eccentricities.

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 *kxx*, *kyy, kxy, cxx, mxx* and *myy*. 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.

#### *3.3. 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:

$$K(\varepsilon) \Big( \times 10^6 \Big) = \begin{cases} 2.71 \varepsilon^2 - 0.703 \varepsilon + 16.9, & \varepsilon \le 0.1, R^2 = 1 \\ 1.70 \varepsilon^3 - 8.67 \varepsilon^2 + 4.21 \varepsilon + 16.5, & 0.1 < \varepsilon \le 0.8, R^2 = 0.9911 \end{cases} \tag{11}$$

$$\mathbf{c}(\varepsilon) \begin{pmatrix} \varepsilon \, \mathbf{1}0^3 \end{pmatrix} = \begin{pmatrix} -16.2\varepsilon^4 + 10.8\varepsilon^3 - 2.20\varepsilon^2 + 0.0767\varepsilon + 3.29, \ \varepsilon \le 0.4, \ R^2 = 0.9995\\ -2.28 + 6.15\varepsilon^2 - 5.43\varepsilon + 4.58, \ 0.4 < \varepsilon \le 0.8, \ R^2 = 0.9997 \end{pmatrix} \tag{12}$$

$$M(\varepsilon) = -5.28\varepsilon^4 + 6.63\varepsilon^3 - 3.58\varepsilon^2 + 0.36\varepsilon + 3.50, \; \varepsilon \le 0.8, R^2 = 0.9905 \tag{13}$$

$$\begin{aligned} k(\varepsilon) \{ \times 10^7 \} &= 14.5\varepsilon^5 - 24.2\varepsilon^4 + 15.5\varepsilon^3 - 3.22\varepsilon^2 + 0.250\varepsilon + 0.517, \; \varepsilon \le 0.8\\ R^2 &= 0.9999 \end{aligned} \tag{14}$$

$$\mathbf{C}(\varepsilon) \Big( \times 10^4 \Big) = 10.7 \varepsilon^4 - 11.6 \varepsilon^3 + 5.57 \varepsilon^2 - 0.682 \varepsilon + 3.01, \ \varepsilon \le 0.8, R^2 = 0.9994 \tag{15}$$

**Figure 18.** "Nominal" force coefficients of concentric annular seal: (**a**) direct stiffness; (**b**) cross damping; (**c**) direct mass; (**d**) cross stiffness; and (**e**) direct damping.

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.

## **4. 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.

#### *4.1. 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% *Cr*. The motion equation is shown as below. Substituting Equations (11)–(16) into Equation (10), the induced fluid forces (*Fx* and *Fy*) are obtained by the nonlinear dynamic model. They are shown in Figure 19 along with the results from transient CFD simulations.

$$\begin{cases} \begin{aligned} \mathbf{x} &= \operatorname{ecos}(\Omega t) \\ \mathbf{y} &= \operatorname{esin}(\Omega t) \end{aligned} \tag{16}$$

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.

**Figure 19.** Fluid forces induced by the circular whirl around seal center.

#### *4.2. 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), fluid-induced forces (*Fx* and *Fy*) are obtained by nonlinear dynamic model. They are compared with the forces from transient CFD simulations, as shown in Figures 20 and 21

**Figure 20.** Fluid forces induced by the circular whirl around seal center (se/*Cr* 20%).

**Figure 21.** Fluid forces induced by the circular whirl around seal center (se/*Cr* 50%).

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