Numerical Investigations on the Blade Tip Clearance Excitation Forces in an Unshrouded Turbine

: The purpose of this study was to investigate the characteristics of the blade tip excitation forces represented as the rotordynamic coe ﬃ cients (sti ﬀ ness and damping coe ﬃ cients) in an unshrouded turbine using the three-dimensional computational ﬂuid dynamic (CFD) numerical method. The blade geometrical parameters were based on a SNECMA transonic experimental rig. The simulations were performed by solving the compressible Reynolds-averaged Navier–Stokes (RANS) equations. The multi-frequency elliptical whirling orbit model and an improved mesh deformation method based on the transient analysis were utilized. The e ﬀ ects of operating conditions on the rotordynamic coe ﬃ cients and the unsteady ﬂow were also found. The results show that the positive direct sti ﬀ ness, which conﬁrmed the direct force contribution in the tip excitation forces and the cross-coupling sti ﬀ ness, were dependent on the whirling frequencies. Damping e ﬀ ects were shown to be negligible. The rotational speed, inlet ﬂow angle, eccentric ratio ( ER ), and mean tip clearance had impacts on the sti ﬀ ness, and some e ﬀ ects of these variables on the rotordynamic coe ﬃ cients were found to be frequency dependent. Additionally, increasing the rotor eccentricity and the mean tip clearance led to the nonuniformity of the circumferential pressure distributions.


Introduction
In an unshrouded turbine with an eccentric rotor, tip clearance excitation forces are generated owing to the variation of blade efficiency with the tip clearance in the circumferential direction. These forces, also known as the "Thomas-Alford forces," could result in self-excited vibration and severe rotordynamic instabilities since they increase with the turbine load. Therefore, to ensure the stability of a rotor-bearing system, it is of great significance to investigate the characteristics of the excitation forces in the unshrouded turbine.
As mentioned above, these forces were first founded by Thomas [1] and Alford [2]. They both theoretically analyzed the mechanism of these forces and proposed respective models. Alford indicated that these forces are related to the stage torque T and mean tip clearance t m as described in Equation (1): where R m is the mean radius, h is the blade height, and β is a correction factor, which is related to the blade efficiency. Based on their theory, several experimental and theoretical analyses were performed to investigate these forces. Urlichs [3] and Wohlrab [4] measured the excitation forces in Figure 1 shows the full 360 • , three-dimensional computational model with an unshrouded concentric rotor geometry. The blade geometrical parameters were based on a SNECMA transonic experimental rig [19] and the detailed information is illustrated in Table 1. The outlet Reynolds number based on the outlet velocity and rotor blade chord was 2.5 × 10 5 . In addition, it can be seen that the inlet extension (200% axial chord length) and outlet extension (500% axial chord length) were constructed to prevent the boundary influence. Different eccentricities and rotor motions were achieved by using the improved mesh deformation method, which is introduced in detail below. Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 15

Computational Method and Boundary Conditions
The commercial software ANSYS-CFX 16.1 was applied to solve the compressible Reynoldsaveraged Navier-Stokes (RANS) equations. To obtain the full set of the frequency-dependent rotordynamic coefficients, transient analysis along with the multi-frequency elliptical whirling model were performed. The corresponding steady solutions were also found using the initial flow fields of the transient solutions. The convective terms and turbulence were found using high-order discretization. The desired converged target was the root mean square (RMS) residuals of each variable reaching a value of 10 −4 or lower, the imbalance of the mass and momentum being less than 0.1%, the response forces caused by the rotor motion approaching periodic vibration, and the difference of results between two adjacent vibration periods being lower than 0.1%. The RNG (Re-Normalization Group) k-ε turbulence model, as well as the scalable wall function, were used to solve the turbulence characteristics. On one hand, RNG k-ε model has a higher precision and a faster convergence speed compared to the standard k-ε model, and it has been validated in many studies [20,21]. On the other hand, it demands a relatively high y+ value compared to the SST (Shear Stress Transfer) model, which is more suitable for mesh movement.
As for the boundary conditions, the absolute total pressure and total temperature were given at the inlet and the static pressure was applied at the outlet. The fixed temperature and no-slip wall conditions were applied. It should be noted that Arts et al. [19] did not take the rotational effects into consideration in their experiment. Therefore, the inlet total pressure (16.613 kPa) with no rotational speed was applied to simulate the experiment for the validation part. For the rest of the analysis, the inlet total pressure (243.5 kPa) and rotational speed (10000 rpm) was applied to keep the mass flow rate same as that in the experiment. Other detailed information is listed in Table 2.

Computational Method and Boundary Conditions
The commercial software ANSYS-CFX 16.1 was applied to solve the compressible Reynoldsaveraged Navier-Stokes (RANS) equations. To obtain the full set of the frequency-dependent rotordynamic coefficients, transient analysis along with the multi-frequency elliptical whirling model were performed. The corresponding steady solutions were also found using the initial flow fields of the transient solutions. The convective terms and turbulence were found using high-order discretization. The desired converged target was the root mean square (RMS) residuals of each variable reaching a value of 10 −4 or lower, the imbalance of the mass and momentum being less than 0.1%, the response forces caused by the rotor motion approaching periodic vibration, and the difference of results between two adjacent vibration periods being lower than 0.1%. The RNG (Re-Normalization Group) k-ε turbulence model, as well as the scalable wall function, were used to solve the turbulence characteristics. On one hand, RNG k-ε model has a higher precision and a faster convergence speed compared to the standard k-ε model, and it has been validated in many studies [20,21]. On the other hand, it demands a relatively high y+ value compared to the SST (Shear Stress Transfer) model, which is more suitable for mesh movement.
As for the boundary conditions, the absolute total pressure and total temperature were given at the inlet and the static pressure was applied at the outlet. The fixed temperature and no-slip wall conditions were applied. It should be noted that Arts et al. [19] did not take the rotational effects into consideration in their experiment. Therefore, the inlet total pressure (16.613 kPa) with no rotational speed was applied to simulate the experiment for the validation part. For the rest of the analysis, the inlet total pressure (243.5 kPa) and rotational speed (10000 rpm) was applied to keep the mass flow rate same as that in the experiment. Other detailed information is listed in Table 2.  Figure 2 presents the mesh details of the rotor stage. A multi-block structural grid was generated for the computational domain using the commercial software ICEM CFD 16.1. The grids near the wall were refined to meet the requirements of the RNG k-ε turbulence model. Three kinds of mesh density in all were used to carry out the mesh sensitivity study and confirm that the results were not dependent on the mesh size. Figure 3 shows the variation of the mass flow rate and stage torque with the element number. It can be seen that the results with 3.67 million elements coincided with the results with 7.30 million elements. Therefore, the rest of the computations were conducted with 3.67 million elements. Furthermore, Figure 4 shows the comparison of the isentropic Mach number distribution between the numerical results using the RNG k-ε turbulence model and the experimental data. The isentropic Mach number Ma is is defined as:

Mesh Verification and Validation
where P s is the static pressure at each location and γ is the ratio of specific heats (γ = 1.4).   Figure 2 presents the mesh details of the rotor stage. A multi-block structural grid was generated for the computational domain using the commercial software ICEM CFD 16.1. The grids near the wall were refined to meet the requirements of the RNG k-ε turbulence model. Three kinds of mesh density in all were used to carry out the mesh sensitivity study and confirm that the results were not dependent on the mesh size. Figure 3 shows the variation of the mass flow rate and stage torque with the element number. It can be seen that the results with 3.67 million elements coincided with the results with 7.30 million elements. Therefore, the rest of the computations were conducted with 3.67 million elements. Furthermore, Figure 4 shows the comparison of the isentropic Mach number distribution between the numerical results using the RNG k-ε turbulence model and the experimental data. The isentropic Mach number Mais is defined as:

Mesh Verification and Validation
where Ps is the static pressure at each location and γ is the ratio of specific heats (γ = 1.4). It can be seen that the numerical results agree well with the experimental data at most streamwise locations. Additionally, the results here are consistent with that in Kim et al. [22], which is also based on the same turbine rig. Therefore, the numerical method with the RNG k-ε model is validated as reliable.

Rotor Motion and Improved Mesh Deformation Method
The multi-frequency elliptical orbit whirling model was used in this paper to obtain the entire set of the frequency-dependent rotordynamic coefficients. This model has been validated in Li et al. [23]. The rotor motion Equations (3) and (4) are defined as harmonic functions with specific vibration amplitudes. Table 3 presents the details of the rotor motion parameters. The vibration amplitudes were obtained based on the small motion theory [24]. For x-direction excitation: and for the y-direction excitation: As mentioned above, the model we used was one of a concentric unshrouded rotor. Therefore, to realize the specific rotor motion, a mesh deformation method [25], which can dynamically control

Rotor Motion and Improved Mesh Deformation Method
The multi-frequency elliptical orbit whirling model was used in this paper to obtain the entire set of the frequency-dependent rotordynamic coefficients. This model has been validated in Li et al. [23]. The rotor motion Equations (3) and (4) are defined as harmonic functions with specific vibration amplitudes. Table 3 presents the details of the rotor motion parameters. The vibration amplitudes were obtained based on the small motion theory [24]. For x-direction excitation: and for the y-direction excitation: As mentioned above, the model we used was one of a concentric unshrouded rotor. Therefore, to realize the specific rotor motion, a mesh deformation method [25], which can dynamically control It can be seen that the numerical results agree well with the experimental data at most streamwise locations. Additionally, the results here are consistent with that in Kim et al. [22], which is also based on the same turbine rig. Therefore, the numerical method with the RNG k-ε model is validated as reliable.

Rotor Motion and Improved Mesh Deformation Method
The multi-frequency elliptical orbit whirling model was used in this paper to obtain the entire set of the frequency-dependent rotordynamic coefficients. This model has been validated in Li et al. [23]. The rotor motion Equations (3) and (4) are defined as harmonic functions with specific vibration amplitudes. Table 3 presents the details of the rotor motion parameters. The vibration amplitudes were obtained based on the small motion theory [24]. For x-direction excitation: and for the y-direction excitation: (Ω i t). As mentioned above, the model we used was one of a concentric unshrouded rotor. Therefore, to realize the specific rotor motion, a mesh deformation method [25], which can dynamically control the movement of the mesh nodes, was needed. It should be noted that the mesh quality, especially for grids of boundary layers, is likely to become worse during the mesh deformation procedure. Moreover, compared to the seal model in Li et al. [23], the unshrouded turbine model was more complicated. Hence, to guarantee the mesh quality, an improved mesh deformation method that could precisely set specific locations of every node at every moment was proposed. Figure 5 presents the two partitions of the cross-section and the illustration of the node movement.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 15 the movement of the mesh nodes, was needed. It should be noted that the mesh quality, especially for grids of boundary layers, is likely to become worse during the mesh deformation procedure. Moreover, compared to the seal model in Li et al. [23], the unshrouded turbine model was more complicated. Hence, to guarantee the mesh quality, an improved mesh deformation method that could precisely set specific locations of every node at every moment was proposed. Figure 5 presents the two partitions of the cross-section and the illustration of the node movement. For the nodes in part A (green area, Rh < R0 < Rt), the specific locations of the nodes at the time t (Xt, Yt) could be described using: and for the nodes in part B (yellow area, Rt < R0 < Rc), the specific locations could be described using: where (X0, Y0) depict the locations of nodes at the initial time, R0 can be calculated using = + , and (x, y) can be obtained via Equations (3) and (4). Additionally, the mesh deformation method could also be used to alter the mean tip clearance without building a new model, which dramatically improved the computing efficiency. For an unshrouded turbine with a new mean tip clearance , the specific locations of the nodes in part A (green area, Rh < R0 < Rt) could be described using: where ⁄ is the correction factor of the vibration amplitude. The specific locations of the nodes in part B (yellow area, Rt < R0 < Rc) could be described using: For the nodes in part A (green area, R h < R 0 < R t ), the specific locations of the nodes at the time t (X t , Y t ) could be described using: and for the nodes in part B (yellow area, R t < R 0 < R c ), the specific locations could be described using: where (X 0 , Y 0 ) depict the locations of nodes at the initial time, R 0 can be calculated using R 0 = X 0 2 + Y 0 2 , and (x, y) can be obtained via Equations (3) and (4). Additionally, the mesh deformation method could also be used to alter the mean tip clearance without building a new model, which dramatically improved the computing efficiency. For an unshrouded turbine with a new mean tip clearance t m , the specific locations of the nodes in part A (green area, R h < R 0 < R t ) could be described using: Appl. Sci. 2020, 10, 1532 where t m /t m is the correction factor of the vibration amplitude. The specific locations of the nodes in part B (yellow area, R t < R 0 < R c ) could be described using: where the first items to the right of the equal sign are the new coordinates of nodes in this part caused by the mean tip clearance change.

Rotordynamic Coefficients Solution
Based on the small motion theory, Childs [24] presented a response-force model, as defined in Equation (9): The rotordynamic coefficients K ij and C ij , including the direct stiffness (K xx , K yy ), the cross coupling stiffness (K xy , K yx ), the direct damping (C xx , C yy ), and the cross coupling damping (C xy , C yx ), could be calculated using the response forces and the corresponding rotor displacements and velocities. For K ij and C ij , the first subscript (i) represents the direction of the response force and the second one (j) represents the direction of the rotor motion. The response forces in this paper are the excitation forces on the rotor blades caused by the rotor eccentricity, which are also known as the "Thomas-Alford" forces. By using the Fast Fourier Transform (FFT), Equation (9) can be transformed into the frequency domain expression as follows, where D ij and G ij are the frequency-domain components of the rotor displacement and response forces, respectively. For D ij and G ij , the first subscript (i) represents the excitation direction and the second one (j) represents the direction of response forces. Then, the corresponding impedances (H ij ) could be solved using Equation (11): Finally, the stiffness and damping coefficients could be solved using Equation (12): Appl. Sci. 2020, 10, 1532 8 of 15 Figure 6 presents the monitored rotor motion data of two cycles in the time domain during the transient calculation process and the corresponding transformed results in the frequency domain with the x-direction excitation. According to the results in the frequency domain, the motion amplitude in the x-direction of each frequency was 5 µm and that in the y-direction was 2.5 µm, which were exactly the target values in Table 3. This indicates that the rotor movement followed Equation (3) Figure 7 shows the response forces in the time and frequency domains with the x-direction excitation. It is shown that the direct response force (Gx) increased first and decreased afterward with the frequency in a small range, where the frequency of the maximum direct force amplitude was about 100-120 Hz. The cross-response force (Gy) increased with frequency in a relatively larger range, where the maximum amplitude appeared at a frequency of 260 Hz. The variations of the direct and cross-response forces with the whirling frequency indicates that the rotordynamic coefficients were dependent on the whirling frequency. Using the monitored data above, the rotordynamic coefficients were solved. Figure 8 presents the variation of the rotordynamic coefficients with frequency. It was shown that the rotordynamic coefficients were dependent on frequency. Moreover, the direct stiffness and damping coefficients were the same (Kxx = Kyy, Cxx = Cyy) and the cross-coefficients were equal in magnitude with opposite sign (Kxy = −Kyx, Cxy = −Cyx). The direct stiffness coefficients Kxx and Kyy were positive at all selected frequencies, which indicates the direct force contributions in the tip clearance excitation forces.   Figure 7 shows the response forces in the time and frequency domains with the x-direction excitation. It is shown that the direct response force (Gx) increased first and decreased afterward with the frequency in a small range, where the frequency of the maximum direct force amplitude was about 100-120 Hz. The cross-response force (Gy) increased with frequency in a relatively larger range, where the maximum amplitude appeared at a frequency of 260 Hz. The variations of the direct and cross-response forces with the whirling frequency indicates that the rotordynamic coefficients were dependent on the whirling frequency. Using the monitored data above, the rotordynamic coefficients were solved. Figure 8 presents the variation of the rotordynamic coefficients with frequency. It was shown that the rotordynamic coefficients were dependent on frequency. Moreover, the direct stiffness and damping coefficients were the same (Kxx = Kyy, Cxx = Cyy) and the cross-coefficients were equal in magnitude with opposite sign (Kxy = −Kyx, Cxy = −Cyx). The direct stiffness coefficients Kxx and Kyy were positive at all selected frequencies, which indicates the direct force contributions in the tip clearance excitation forces. Using the monitored data above, the rotordynamic coefficients were solved. Figure 8 presents the variation of the rotordynamic coefficients with frequency. It was shown that the rotordynamic coefficients were dependent on frequency. Moreover, the direct stiffness and damping coefficients were the same (K xx = K yy , C xx = C yy ) and the cross-coefficients were equal in magnitude with opposite sign (K xy = −K yx , C xy = −C yx ). The direct stiffness coefficients K xx and K yy were positive at all selected frequencies, which indicates the direct force contributions in the tip clearance excitation forces. Thomas and Alford did not predict the direct forces in their model; instead, these forces were confirmed in later experiments. It also should be noted that since the direct stiffness coefficients decreased with the frequency, the direct stiffness might be zero at a frequency around 270 Hz. Moreover, the cross-stiffness coefficient K xy was positive at all frequencies and it increased with the frequency, which indicates that the tip clearance excitation forces could promote forward rotor whirl. Compared to the stiffness coefficients, the damping coefficients were quite small, which indicates that there was no damping effect. Hence, the damping coefficients are not discussed in the following sections.

Effects of the Operating Conditions
In this section, the effects of the operating conditions on the stiffness coefficients are presented. Figure 9 shows the variation of stiffness coefficients with frequency at different rotational speeds. Both the direct and cross-coupling stiffness coefficients decreased with the rotational speed, which indicates that the increase of rotational speed caused a decrease of the tip clearance excitation effect. This was because the torque on the rotor decreased with increasing rotational speed and hence caused a decrease in the response forces in accordance with Equation (1). Moreover, the direct stiffness coefficients at lower frequencies were more easily affected compared to the results at high frequencies, while the decreases of the cross-coupling stiffness coefficients appeared independent of the frequency.

Effects of the Operating Conditions
In this section, the effects of the operating conditions on the stiffness coefficients are presented. Figure 9 shows the variation of stiffness coefficients with frequency at different rotational speeds. Both the direct and cross-coupling stiffness coefficients decreased with the rotational speed, which indicates that the increase of rotational speed caused a decrease of the tip clearance excitation effect. This was because the torque on the rotor decreased with increasing rotational speed and hence caused a decrease in the response forces in accordance with Equation (1). Moreover, the direct stiffness coefficients at lower frequencies were more easily affected compared to the results at high frequencies, while the decreases of the cross-coupling stiffness coefficients appeared independent of the frequency. Figure 10 shows the effects of the inlet flow angle on the stiffness coefficients. It was shown that a higher inlet flow angle led to lower direct stiffness coefficients and higher cross-coupling stiffness coefficients. This was because changing the inlet flow angle only caused a slight change in the blade efficiency and had little impact on the rotor torque. Hence, the effect of the inlet flow angle was much smaller than that of the rotational speed, which indicates that the tip clearance excitation could not be reduced by altering the inlet flow angle. In addition, the results at lower frequencies showed less dependency on the inlet flow angle.
indicates that the increase of rotational speed caused a decrease of the tip clearance excitation effect. This was because the torque on the rotor decreased with increasing rotational speed and hence caused a decrease in the response forces in accordance with Equation (1). Moreover, the direct stiffness coefficients at lower frequencies were more easily affected compared to the results at high frequencies, while the decreases of the cross-coupling stiffness coefficients appeared independent of the frequency.  Figure 10 shows the effects of the inlet flow angle on the stiffness coefficients. It was shown that a higher inlet flow angle led to lower direct stiffness coefficients and higher cross-coupling stiffness coefficients. This was because changing the inlet flow angle only caused a slight change in the blade efficiency and had little impact on the rotor torque. Hence, the effect of the inlet flow angle was much smaller than that of the rotational speed, which indicates that the tip clearance excitation could not be reduced by altering the inlet flow angle. In addition, the results at lower frequencies showed less dependency on the inlet flow angle.  Figure 11 presents the variation of the stiffness coefficients with frequency with different rotor eccentricities in the y-direction. The eccentricity could be introduced using a mesh deformation. For an unshrouded turbine with a specific rotor eccentric ratio (ER) in the y-direction, the locations of nodes in part A were: The locations of the nodes in part B were: According to the results, the eccentricity had little impact on the direct stiffness coefficient, especially for the result at low frequencies. When the rotor eccentric ratio (ER) increased, the direct stiffness coefficient showed an accelerated decreasing trend from ER = 0−0.3 and then a decelerated decreasing trend from ER = 0.3−0.7. For the cross-coupling stiffness coefficients, the results at low frequencies increased with the eccentricity from ER = 0−0.3 and decreased rapidly from ER = 0.3−0.7. However, the results at high frequencies decreased continually with the eccentricity from ER = 0−0.7. These results indicate that the tip clearance excitation effect reduced when the unshrouded turbine rotor was at a high eccentricity.  Figure 11 presents the variation of the stiffness coefficients with frequency with different rotor eccentricities in the y-direction. The eccentricity could be introduced using a mesh deformation. For an unshrouded turbine with a specific rotor eccentric ratio (ER) in the y-direction, the locations of nodes in part A were: The locations of the nodes in part B were: According to the results, the eccentricity had little impact on the direct stiffness coefficient, especially for the result at low frequencies. When the rotor eccentric ratio (ER) increased, the direct stiffness coefficient showed an accelerated decreasing trend from ER = 0−0.3 and then a decelerated decreasing trend from ER = 0.3−0.7. For the cross-coupling stiffness coefficients, the results at low frequencies increased with the eccentricity from ER = 0−0.3 and decreased rapidly from ER = 0.3−0.7. However, the results at high frequencies decreased continually with the eccentricity from ER = 0−0.7. These results indicate that the tip clearance excitation effect reduced when the unshrouded turbine rotor was at a high eccentricity.  Figure 12 shows the rotordynamic coefficient with different mean tip clearances. It was shown that the direct stiffness at selected frequencies decreased dramatically with the mean tip clearance, especially for the results with tm/h = 2%-3%. In addition, the results at low frequencies were more easily affected by the mean tip clearance than at high frequencies. Moreover, it should be noted that the direct stiffness coefficients became negative at frequencies over 220 Hz for tm/h = 3%−4%. The variation of the cross-coupling stiffness showed a different pattern compared to that of the direct stiffness. When tm/h increased from 1% to 2%, the cross-coupling stiffness coefficients increased for all frequencies, which indicates increased instability forces on the rotor. This was because too small of a tip clearance leads to a change of tip leakage and hence causes a change of torque on the blades. When tm/h continued to increase from 2%-3%, the cross-coupling stiffness at low frequency (Ω < 160 Hz) increased and that at high frequency (Ω > 160 Hz) decreased. Meanwhile, the dependency of the cross-coupling stiffness with the whirling frequency reduced a lot. As tm/h increased further from 3% to 4%, the cross-coupling stiffness decreased for all frequencies. Therefore, when the mean tip clearance was too small (tm/h = 1%), reducing the tip clearance excitation forces by increasing the mean tip clearance may not be feasible. However, when the mean tip clearance was relatively large (tm/h = 3%), increasing the mean tip clearance would be a useful way to reduce the instability forces caused by a fluid.

Static Pressure Distributions
To better understand the flow characteristics in this turbine and the cause of the tip clearance excitation forces, the cross-sections are illustrated in Figure 13a, while Figure 13b-d shows the  Figure 12 shows the rotordynamic coefficient with different mean tip clearances. It was shown that the direct stiffness at selected frequencies decreased dramatically with the mean tip clearance, especially for the results with t m /h = 2%-3%. In addition, the results at low frequencies were more easily affected by the mean tip clearance than at high frequencies. Moreover, it should be noted that the direct stiffness coefficients became negative at frequencies over 220 Hz for t m /h = 3%−4%. The variation of the cross-coupling stiffness showed a different pattern compared to that of the direct stiffness. When t m /h increased from 1% to 2%, the cross-coupling stiffness coefficients increased for all frequencies, which indicates increased instability forces on the rotor. This was because too small of a tip clearance leads to a change of tip leakage and hence causes a change of torque on the blades. When t m /h continued to increase from 2%-3%, the cross-coupling stiffness at low frequency (Ω < 160 Hz) increased and that at high frequency (Ω > 160 Hz) decreased. Meanwhile, the dependency of the cross-coupling stiffness with the whirling frequency reduced a lot. As t m /h increased further from 3% to 4%, the cross-coupling stiffness decreased for all frequencies. Therefore, when the mean tip clearance was too small (t m /h = 1%), reducing the tip clearance excitation forces by increasing the mean tip clearance may not be feasible. However, when the mean tip clearance was relatively large (t m /h = 3%), increasing the mean tip clearance would be a useful way to reduce the instability forces caused by a fluid.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 15 (a) (b) Figure 11. Variation of the rotordynamic coefficient with rotor eccentricity: (a) direct stiffness coefficient and (b) cross-coupling stiffness coefficient. Figure 12 shows the rotordynamic coefficient with different mean tip clearances. It was shown that the direct stiffness at selected frequencies decreased dramatically with the mean tip clearance, especially for the results with tm/h = 2%-3%. In addition, the results at low frequencies were more easily affected by the mean tip clearance than at high frequencies. Moreover, it should be noted that the direct stiffness coefficients became negative at frequencies over 220 Hz for tm/h = 3%−4%. The variation of the cross-coupling stiffness showed a different pattern compared to that of the direct stiffness. When tm/h increased from 1% to 2%, the cross-coupling stiffness coefficients increased for all frequencies, which indicates increased instability forces on the rotor. This was because too small of a tip clearance leads to a change of tip leakage and hence causes a change of torque on the blades. When tm/h continued to increase from 2%-3%, the cross-coupling stiffness at low frequency (Ω < 160 Hz) increased and that at high frequency (Ω > 160 Hz) decreased. Meanwhile, the dependency of the cross-coupling stiffness with the whirling frequency reduced a lot. As tm/h increased further from 3% to 4%, the cross-coupling stiffness decreased for all frequencies. Therefore, when the mean tip clearance was too small (tm/h = 1%), reducing the tip clearance excitation forces by increasing the mean tip clearance may not be feasible. However, when the mean tip clearance was relatively large (tm/h = 3%), increasing the mean tip clearance would be a useful way to reduce the instability forces caused by a fluid.

Static Pressure Distributions
To better understand the flow characteristics in this turbine and the cause of the tip clearance excitation forces, the cross-sections are illustrated in Figure 13a, while Figure 13b-d shows the

Static Pressure Distributions
To better understand the flow characteristics in this turbine and the cause of the tip clearance excitation forces, the cross-sections are illustrated in Figure 13a, while Figure 13b-d shows the normalized static pressure contours on the selected cross-sections. At this moment, the rotor displacement in the x-direction caused by the excitation was at the maximum (x = 65 µm, ER = 0.13, shown in Figure 6a) and that in the y-direction was about 350 µm (ER = 0.7) due to the rotor eccentricity. Meanwhile, the rotor was moving in the positive y-direction.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 15 normalized static pressure contours on the selected cross-sections. At this moment, the rotor displacement in the x-direction caused by the excitation was at the maximum (x = 65 μm, ER = 0.13, shown in Figure 6a) and that in the y-direction was about 350 μm (ER = 0.7) due to the rotor eccentricity. Meanwhile, the rotor was moving in the positive y-direction. According to Figure 13b, the eccentricity had an impact on the pressure distributions of upstream flow, which validates the hypothesis in Song and Martinez-Sanchez [8]. The static pressure at the smaller tip gap was relatively higher than that at the larger tip gap when comparing the result on the upper and right surfaces with that on the lower and left surfaces. According to Figure 13c, a larger tip gap led to a larger tip leakage flow, and therefore caused a decrease of static pressure at the tip gap. Due to the circumferentially uneven upstream flow and tip leakage flow, the torque on the blade was also circumferentially non-uniform, and hence caused the tip clearance excitation forces. According to Figure 13d, at the larger tip gap region, a low-pressure area caused by a larger tip leakage flow can be seen. Figure 14 shows the distributions of selected variables of numerical results (left) and experimental results in Song [9] (right). In the right figure, ρ is the density, m is the mass flow rate, U is the linear velocity at the average radius, and the subscript av indicates circumferential averaged value. Comparing these two figures showing the phase differences between the local torque and the local tip clearance, the local pressure and the local tip clearance were quite similar. According to the numerical results, the static pressure on the tip basically decreased with the local tip clearance. However, it should be noted that the phase difference between the local tip clearance and the local According to Figure 13b, the eccentricity had an impact on the pressure distributions of upstream flow, which validates the hypothesis in Song and Martinez-Sanchez [8]. The static pressure at the smaller tip gap was relatively higher than that at the larger tip gap when comparing the result on the upper and right surfaces with that on the lower and left surfaces. According to Figure 13c, a larger tip gap led to a larger tip leakage flow, and therefore caused a decrease of static pressure at the tip gap. Due to the circumferentially uneven upstream flow and tip leakage flow, the torque on the blade was also circumferentially non-uniform, and hence caused the tip clearance excitation forces. According to Figure 13d, at the larger tip gap region, a low-pressure area caused by a larger tip leakage flow can be seen. Figure 14 shows the distributions of selected variables of numerical results (left) and experimental results in Song [9] (right). In the right figure, ρ is the density, m is the mass flow rate, U is the linear velocity at the average radius, and the subscript av indicates circumferential averaged value. Comparing these two figures showing the phase differences between the local torque and the local tip clearance, the local pressure and the local tip clearance were quite similar. According to the numerical results, the static pressure on the tip basically decreased with the local tip clearance. However, it should be noted that the phase difference between the local tip clearance and the local torque on the blade was not 180 • , but instead approximately 120 • , which indicates that the lowest local tip clearance could not lead to the highest local torque on the blade. Based on this, a positive direct stiffness, which indicates the direct force contribution in the tip excitation forces, was generated.

Conclusions
A three-dimensional CFD numerical method, along with an improved mesh deformation technique based on the transient analysis and multi-frequency elliptical orbit whirling model, was proposed to investigate the tip clearance excitation forces in an unshrouded turbine. The tip clearance excitation forces were represented as the rotordynamic coefficients (stiffness and damping coefficients) in this paper. The effects of several variables, including the whirling frequency, rotational speed, inlet flow angle, eccentricity, and the mean tip clearance, on the rotordynamic coefficients and aerodynamic performance were analyzed in detail. The conclusions are as follows: (1) Positive direct stiffness coefficients confirmed the existence of a direct force contribution in the tip clearance excitation forces. The damping effects were found to be quite small and could be ignored. Both the stiffness and damping showed a dependency on the rotor whirling frequency. Specifically, the direct stiffness increased with the frequency, while the cross-coupling stiffness decreased with the whirling frequency.
(2) The increase of the rotational speed caused a considerable decrease of the direct stiffness and cross-coupling stiffness due that the torque on the rotor decreasing with the rotational speed, which indicates a decreasing tip clearance excitation effect. The increase of the mean tip clearance led to a decreasing direct stiffness and a complicated change of the cross-coupling stiffness. For turbines with a small mean tip clearance (tm/h = 1%), increasing the mean tip clearance caused the increase of cross-coupling stiffness, which led to rotor instability. For a turbine with a higher mean tip clearance (tm/h = 3%), the cross-coupling stiffness decreased with the increase of the mean tip clearance, which enhanced the rotor stability. The effects of the inlet flow angle and eccentricity on the rotordynamic coefficients were relatively smaller.
(3) The rotor eccentricity caused a pressure redistribution of both upstream and downstream flow. A higher eccentricity and higher mean tip clearance resulted in a greater nonuniformity of the circumferential static pressure distribution. The characteristics of leakage flow changed with the local tip clearance. The phase difference between the local tip clearance and the local torque on the blade was about 120°, which led to the direct force contribution.

Conclusions
A three-dimensional CFD numerical method, along with an improved mesh deformation technique based on the transient analysis and multi-frequency elliptical orbit whirling model, was proposed to investigate the tip clearance excitation forces in an unshrouded turbine. The tip clearance excitation forces were represented as the rotordynamic coefficients (stiffness and damping coefficients) in this paper. The effects of several variables, including the whirling frequency, rotational speed, inlet flow angle, eccentricity, and the mean tip clearance, on the rotordynamic coefficients and aerodynamic performance were analyzed in detail. The conclusions are as follows: (1) Positive direct stiffness coefficients confirmed the existence of a direct force contribution in the tip clearance excitation forces. The damping effects were found to be quite small and could be ignored. Both the stiffness and damping showed a dependency on the rotor whirling frequency. Specifically, the direct stiffness increased with the frequency, while the cross-coupling stiffness decreased with the whirling frequency.
(2) The increase of the rotational speed caused a considerable decrease of the direct stiffness and cross-coupling stiffness due that the torque on the rotor decreasing with the rotational speed, which indicates a decreasing tip clearance excitation effect. The increase of the mean tip clearance led to a decreasing direct stiffness and a complicated change of the cross-coupling stiffness. For turbines with a small mean tip clearance (t m /h = 1%), increasing the mean tip clearance caused the increase of cross-coupling stiffness, which led to rotor instability. For a turbine with a higher mean tip clearance (t m /h = 3%), the cross-coupling stiffness decreased with the increase of the mean tip clearance, which enhanced the rotor stability. The effects of the inlet flow angle and eccentricity on the rotordynamic coefficients were relatively smaller.
(3) The rotor eccentricity caused a pressure redistribution of both upstream and downstream flow. A higher eccentricity and higher mean tip clearance resulted in a greater nonuniformity of the circumferential static pressure distribution. The characteristics of leakage flow changed with the local