Dynamic Behaviour Analysis of Turbocharger Rotor-Shaft System in Thermal Environment Based on Finite ElementMethod

'e stable operation of a high-speed rotating rotor-bearing system is dependent on the internal damping of its materials. In this study, the dynamic behaviours of a rotor-shaft system with internal damping composite materials under the action of a temperature field are analysed. 'e temperature field will increase the tangential force generated by the internal damping of the composite material. 'e tangential force will also increase with the rotor speed, which can destabilise the rotor-shaft system. To better understand the dynamic behaviours of the system, we introduced a finite element calculation model of a rotor-shaft system based on a 3D high-order element (Solid186) to study the turbocharger rotor-bearing system in a temperature field. 'e analysis was done according to the modal damping coefficient, stability limit speed, and unbalance response.'e results show that accurate prediction of internal damping energy dissipation in a temperature field is crucial for accurate prediction of rotor dynamic performance. 'is is an important step to understand dynamic rotor stress and rotor dynamic design.


Introduction
Turbochargers are mechanical devices that can improve fuel efficiency and reduce greenhouse gas emissions. e core component of a turbocharger is a rotor composed mainly of a turbine and a compressor. e turbine is crucial because it can recover energy from the exhaust gas and increase the intake air volume by driving the compressor [1]. It is characterized by light weight and high speed. To meet requirements for higher speeds, greater power density, and the ability to adapt to harsh operating environments for a long time, structural designers use various materials to manufacture the rotor parts. is introduces high requirements for the stable rotation of the rotor. erefore, to analyse the stability of a rotor and its dynamic characteristics, it is necessary to accurately predict the damping effect. Damping is divided into external damping, such as bearing damping, and internal damping, such as material damping, which in turn is mainly modelled by viscous damping and hysteretic damping [2]. In a dynamic composite rotor, internal damping is meaningful only when the matrix can damp [3]. Also, most single materials, such as metals, have vibration damping characteristics similar to those of hysteretic internal damping but not viscous internal damping [4].
Many researchers have studied the effect of material damping on rotor dynamics and stability behaviour. Sin [5] studied the instability phenomenon of a composite shaft with internal damping and calculated the natural frequency and instability threshold of the shaft by using the finite element model for beams. ey found that the stacking sequence of layers of composite materials, the arrangement directions of material fibers, and the transverse shear forces of materials affected the natural frequency and instability threshold of the shaft. Wittgren [6] studied the stability of flexible rotors with flexible supports. ey found that the correct combination of the asymmetry of the internal damping of the shaft material and the anisotropy of bearings could significantly reduce the amplitude at critical speed. Montagnier [7] used Euler-Bernoulli shaft model to analyse the critical speed of a flexible internal resistance rotor supported by elastic bearings to obtain the stable working range of the rotor. ey found that most of the rotor instability was caused by the internal damping of the material. Increasing the bearing damping could make the rotor run more stably, and they established a rotor instability analysis model that included the internal resistance of the material. Vitta [8] used linear and nonlinear evaluation methods to analyse the effect of internal damping on the dynamic behaviour of a shaft. e results show that the critical, subcritical, and stable limit speeds of a rotor can be obtained by linear evaluation. Newkirk [9] studied internal damping and found that a rotor may vibrate violently at a speed higher than the first critical speed. Genta [10] explained that the hysteretic damping of rotating structures is stable in the subcritical range but unstable in the supercritical range. All those researchers proved that the instability speed is higher than the first-order critical speed when considering the hysteretic damping of rotating structures. Also, many researchers have studied the combined effects of material damping and bearing damping. e results show that increasing the bearing damping can improve the stability of the rotor, but increasing the material damping can reduce the instability threshold [6]. e effect of the surrounding temperature on material damping of a turbocharger rotor cannot be disregarded in high-temperature environments. It is very important to understand the dynamic behaviours of the structure in various temperatures when designing a rotor to operate in thermal extremes. San [11] established a nonlinear rotorbearing finite element model. e model considered the thermal effects of lubricating oil and the thermal expansions of the rotor shaft and bearing, and San Andres analysed the natural frequency, stability, and unbalanced response of the rotor.
e results show that the natural frequency and synchronous speed amplitude obtained by the simulation were completely consistent with the experimental values, which verified the feasibility of the finite element model. Zych [12] used a finite element method to calculate the thermal stress of a radial axial microturbine in a hightemperature environment. In the calculation, they considered the mass of the disc, the rotation speed of the rotor, and the complex shape at the rear of the disc. eir work showed that numerical calculation helps to choose the best optimization method, and they reduced the turbine's von Mises stress by approximately 45%. Jeyaraj [13] uniformly heated a thin plate that incorporated various internal damping composite materials and did free vibration and forced vibration analyses.
e study found that, with increasing temperature, the vibration amplitude (response) of the thin plate structure decreased, but the modal loss factor increased significantly with increasing temperature, thereby reducing the vibration amplitude. Its response frequency also decreased with increasing temperature. Guo [14] did thermal vibration analysis of a rotating beam structure with constrained layer damping. e study found that, as the temperature increased, the modal frequency of the beam structure decreased accordingly, and the damping ratio increased accordingly. at study provides a basis for dynamic analysis of high-speed rotating blades in various thermal environments. However, the number of articles on this subject is limited. erefore, this research studied the dynamic characteristics of internal material damping and an oil film force turbocharger rotor-bearing system under thermal environments (temperature fields). We used the conjugate heat transfer (CHT) method to simulate the temperature field of the solid part of a rotorshaft system. e temperature field was coupled with the finite element model of the rotor, compared with the instance without considering the temperature field. e rotor finite element is verified by experiment.

Fluid and Heat Transfer Sections.
Before the thermal modal analysis, the aerodynamic thermal analysis was performed.
e current industry standard modelling approaches assume the turbine and compressor operate under adiabatic conditions [15].
e CHT simulation method [16][17][18][19][20][21][22] is used to obtain the temperature distribution of the rotor, and the node temperature is provided for the modal simulation part.

Turbocharge and Compressor Geometry.
e turbocharger had an impeller with 10 blades, and the compressor had an impeller with 6 blades and 6 splitters. e turbocharger and compressor parameters are shown in Table 1. To improve the accuracy of the calculations, we modelled both the turbocharger and the compressor using full passages. e solid parts and the air flow are shown in Figure 1.

Grid Generation.
CHT involves the direct coupling of fluids and solids. ICEM grid discretisation software (Ansys, Inc., USA) uses the same numerical principles and grid discretisation for both regions.
is allows the noninterpolated exchange of heat flux between adjacent grids [23]. e calculation accuracy of CHT is very sensitive to the resolution of the fluid boundary layer grid. erefore, the dimensionless distance y + � 1 or less (in this paper y + � 1) of the wall distance of the first layer of the grid can determine the local heat flux with enough accuracy. As shown in Figures 2 and 3, the fluid domain, solid domain (rotor impeller), and boundary layer grid were generated for CHT calculations. e total number of global grids was 12,329,631, with 2,067,903 global grid nodes. Following the shape of the outer contour of the rotor, we used a triangular tetrahedral mesh with good adaptability to the outer contour, and we used smoothing to optimize the mesh. e impeller grid and the fluid grid were connected in a general grid interface mode in Ansys CFX software.

Boundary Conditions.
We determined the boundary conditions for aerodynamic thermal analysis using experimental data provided by the turbocharger company; the boundary conditions are shown in Table 2. Shock and Vibration     e heat transfer of the surrounding air was disregarded, and the out-wall of the turbocharger was assumed to be adiabatic. We applied "no-slip" boundary conditions to all inner walls. An interface was added between the rotating domain and the fixed domain, and the interface was connected by a "frozen rotor" [24].

Numerical
Methodology. CHT refers to a coupled heat transfer phenomenon in which the thermal properties of two materials occur through a medium or in direct contact. e CHT method can calculate the heat transfer between fluid and solid and calculate the temperatures of fluids and solids at the same time. In this study, we used the commercial software Ansys CFX for numerical simulation. CFX is a computational fluid dynamics software package based on the control volume method to solve Navier-Stokes equations. In the fluid domain, the mass conservation, momentum, and energy transport equations are described as where u f , ρ f , p, τ, λ f , T, h tot , and t represent the velocity vector, density, pressure, stress tensor, thermal conductivity, temperature, total enthalpy, and time of the fluid, respectively. In a solid domain, the conservation of energy equation can explain the heat transfer caused by solid motion, conduction, and a volume heat source. e energy equation is where u s is the velocity vector of the solid, u s � u f ; S E is the optional volume heat source, S E � ∇ · (u s · τ); ρ s is the density of the solid; and λ s is the thermal conductivity of the solid. is study did not directly solve (1)-(4). Instead, they were converted to the steady-state Reynolds average Navier-Stokes method to calculate turbulence. e fluid medium (exhaust gas and air) is an ideal gas, and the shear stress transmission (SST) turbulence model was used because the SST turbulence model has good accuracy for CHT calculations [20]. at model has both the accuracy of the k − ω model in high-pressure gradient flow boundary layer prediction and the stability of the k − e model in mainstream prediction [18].

Rotor-Shaft System Geometry.
e model of the rotorshaft system was provided by the turbocharger company and agreed with the real one. e rotor-shaft system comprised a turbine wheel, a compressor wheel, a rotating shaft, a floating ring bearing, a thrust bearing, a seal sleeve, and a nut, as Figure 4 shows. e turbine wheel and shaft were connected by friction welding, whereas the compressor was axially and circumferentially fixed by a left-hand nut. e materials of the main parts of the rotor-shaft system are shown in Table 3, based on the actual situation.
is study investigated the effect of temperature and found that the material properties changed as the temperature changed. e function curves of the material properties of each part as a function of temperature are shown in Figure 5 [25,26].
Publication [29] proposes a set of damping coefficients, and we verified the influence of damping coefficients on the stability limit displacement of the rotor through theoretical analysis and numerical simulation. We selected this set of damping coefficients based on the theory of the publication. Table 4 shows what we assumed was the material damping coefficient of each part of the rotor.

Floating Ring Bearing Parameters.
e rotor shaft was supported by a floating ring bearing whose detailed parameters are shown in Table 5.
To aid the simulation, we combined the stiffness and damping coefficients of the inner and outer oil film into a total impedance (equivalent stiffness and damping coefficients) according to [27,28]. Combining the stiffness and damping coefficients of the inner and outer oil film provided by the turbocharger company with the parameters of the floating ring bearing, we calculated the total impedance using Dyrobes rotor dynamics software (Dyrobes, USA). Figure 6 shows the total impedance (equivalent stiffness and damping coefficients) of the floating ring bearing at various speeds.

Finite Element Model of Rotor.
is study used a largescale general software Ansys Workbench for simulation. To obtain the motion trace on the axis, the rotor shaft was equally divided into four parts along the axis direction, as shown in Figure 7. e blue lines (solid and dotted) in the figure represent the division positions. e solid model was discretised, and an Ansys Solid186 element was used. e Solid186 element is a high-order, 3D 20-node element with quadratic displacement characteristics. e element is defined by 20 nodes, each node with three degrees of freedom (translation of each node in the x-, y-, and z-directions).
e element supports plasticity, superelasticity, creep, stress hardening, large deflection, and large strain capacity. In the shaft part, the grid was generated by a sweep, as shown in Figure 8, depicting the "i th " element after segmentation.
We transformed (deformed) the original hexahedron element to an approximately one-quarter cylinder element. In Figure 8, the red dotted line represents the actual edge of the element, and the red solid line represents the edge of the element. e advantage of a solid model is that it retains the properties of the original element and describes the outer contour of the shaft model well. 4 Shock and Vibration To adapt to the more complicated outer contour of the blade, tetrahedral, pyramid, and prism methods were used to generate the unit in the blade and the remaining part. Figure 9 shows the (i + 1), (i + 2), and (i + 3) elements generated by three methods. e effects of rotational inertia, translational inertia, gyro moment, support stiffness, support damping, material damping, rotational damping, and thermal stress stiffness were considered in the model to establish the equation of motion. e rotor model was discretised into N e elements, and the total number of nodes became a.
e node displacement vector of the i th element is q e � x iI , y iI , z iI , x iJ , y iJ , z iJ , . . . , x iμ , y iμ , z iμ , . . . , x iA , where μ � I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y, Z, A, B is the node number and the superscript T is the matrix/vector transpose. ese writing formats are also mentioned later in the text.
After combining the governing equations of all elements and combining the boundary conditions, the equation of motion of the rotor-bearing system is where Here, M is a mass matrix that includes mass components caused by translation along the axis (M trs ) and rotation along the axis (M rot ). C is a damping matrix (asymmetric matrix coefficient of the velocity vector) that includes a skewsymmetric gyro matrix (G), a constant structural damping matrix (C con � βK con ) with hysteresis characteristics caused by internal friction of the material and damping caused by the floating ring bearing support matrix (C brg ). K is the stiffness matrix (matrix coefficient of the displacement vector) that includes the stiffness matrix (K brg ) caused by the support of the floating ring bearing, the spin-softening bending matrix (K sh ) generated by the rotor rotation, the stiffness matrix (K tem ) of the elastic modulus, the stiffness matrix (K temσ ) caused by the nonuniform thermal stress, and thermal deformation caused by temperature changes. e size of all matrices is 3a × 3a. Ω is the rotor speed in rpm. β is the material damping (Rayleigh damping), stiffness damping coefficient. e total displacement and global force vectors q(t) and f(t) are expressed as      Shock and Vibration Writing (6) in state-space form yields where matrices A and B; state vector w(t); and force vector F(t) are where 0 is an empty matrix of size 3a × 3a and an empty vector of size 3a × 1. Figure 10 shows a simplified finite element model of the rotor-shaft system. e grey-blue points in the figure are used to observe the modal shapes, line trajectories, and unbalanced harmonic response analysis. Figure 11 shows the rotor finite element global grid. e total number of grids was 3,176,485, with 546,067 grid nodes. In the geometry software, each portion of the rotor was treated as a part to avoid extra contact stiffness in the finite element simulation and prevent affecting the simulation accuracy.
To aid the comparison, we considered two cases in this study: case 1, the rotor-shaft system without the influence of temperature, and case 2, the rotor-shaft system affected by temperature. In the following sections, r-th forward and r-th backward modes are represented as "r-F" and "r-B" modes, respectively.

Heat Transfer Results and Experimental Verification.
After the calculation, we verified the numerical results of the mass flow and the supercharger ratio (expansion ratio). e numerical results agreed with the experimental data, and the error control was approximately 7%, as Figure 12 shows. e temperature distribution of the rotor was obtained through postprocessing, as Figure 13 shows; the temperature distribution law was the same as that in the Bohn [19] study result and is shown in Figure 14. Figures 13 and 14 show that the highest temperature was at the outer edge of the turbine wheel blades and that the lowest temperature was at the compressor wheel. Overall, the rotor had a large temperature gradient.

Campbell Diagram and Stability Diagram of Normal
Temperature Environment. Figure 15 shows the Campbell diagram of the rotor-shaft system in case 1, that is, without considering the temperature. e figure was         Figure 14: Rotor temperature distribution of Bohn [19].  marked as Ω cr1 and Ω cr2 , respectively, and were calculated as follows: Ω cr1 � 5, 011.5 rpm and Ω cr2 � 43, 622 rpm. Because the backward whirl could not be excited by the unbalanced force, it is not marked or considered.
When the speed was in the critical speed range between Ω cr1 and Ω cr2 in the 1F mode, the whirling frequency had a clear increasing trend. at trend may have been caused by floating ring bearings, because in the critical speed range of Ω cr1 and Ω cr2 , the stiffness and damping of the oil film changed sharply. When the critical speed of Ω cr2 was exceeded, the stiffness and damping of the oil film changed smoothly, and the eddy frequency in the corresponding 1F mode also increased at a constant rate. Figure 16 shows that the curve of the damping ratio of the four modes in case 1 changed with the rotation speed. e modal damping ratio is the ratio of actual damping to the critical damping of a certain mode, and the damping ratio is the decisive parameter to describe the stability. A positive damping ratio means that the rotor-shaft system was in a stable region when the vibration energy was dissipated. A negative damping ratio means that the rotation energy supported the rotor rotation by adding vibration energy, resulting in the rotor-shaft system becoming unstable. e point where the damping ratio curve intersects with the zero line is the stable limit speed (SLS). e figure shows that the rotor-shaft system became unstable when SLS 1 � 56, 550 rpm in 1F mode and then became unstable when SLS 2 � 97, 946 rpm in 2F mode. Figure 17 shows the four modes close to the SLS 1 stability limit speed at 55,000 rpm. e black dotted line represents the axis, the combination of the orange line and the orange sphere represents the starting position of the trajectory, and the incomplete track curve represents the whirl direction. e figure shows that the 1F mode was cylindrical, thus combining rigid body motion with rotor bending. Furthermore, the rotor was close to the unstable state at this time. In contrast, the 2F mode was a conical, bending mode. e figure shows that the 3B mode was conical, with both rigid body motion and rotor bending. Lastly, the 4B mode was cylindrical, with both rigid body motion and rotor bending. Figure 18 shows the four modes close to the SLS 2 stability limit speed at 97,946 rpm. e rotor was close to the unstable state at this time. e 1F and 3B modes were cylindrical and conical, combining rigid body motion and rotor bending. In contrast, the 2F mode was a conical, bending mode.

Effect of Temperature on the Campbell Diagram and
Stability Diagram. We exported the temperature data of each node in Figure 12 and inserted the temperature value into the corresponding rotor-shaft system node through the Ansys software commands function. We analysed the critical speed and stability of the rotor-shaft system in the thermal environment. Figure 19 shows that temperature affected the intrinsic frequency (whirl frequency). Compared with case 1, the whirl-frequency curves in the 1F, 2F, and 3B modes were slightly lower, whereas for the 4B mode, the whirl-frequency curve was much lower. In the figure, the starting point at the left end (when the rotating speed was 0 rpm) was taken as the observation object. In case 1, the frequency was 2,300.9 Hz, whereas in case 2, the frequency decreased to 2,182.5 Hz. Because the temperature changed the material's elastic modulus and Poisson's ratio, it affected its structural stiffness. For higher-order frequencies, the effect on frequency became greater as the order increased. However, the effect on the critical speed, Ω cr1 , did not change much, whereas Ω cr2 decreased from 43,622 rpm to 41,815 rpm. is decrease was also caused by the intrinsic frequency drop. Figure 20 shows the stability diagram in case 2 and that the temperature had an important effect on the damping ratio and SLS. Compared with case 1, SLS 1 of the 1F mode decreased from the original 56,550 rpm to 55,078 rpm and made the 1F mode enter the unstable region ahead of time.
e 2F mode change was more obvious, and the same SLS 2 decreased from the original 97,946 rpm to 95,969 rpm. e reason for this can be explained by the following formula: ξ r � (α/2ω hr ) + (βω hr /2). Here, ξ r is the damping ratio coefficient of the r th mode, α is the mass damping coefficient, β is the stiffness damping coefficient, and ω hr is the natural frequency of the r th mode. In most cases, α (α � 0) can be disregarded; that is, ξ r � βω hr /2. β did not change, whereas ω hr decreased because of the effect of temperature, resulting in a decrease in ξ r . In other words, SLS 1 and SLS 2 in the 1F and 2F modes were reduced. erefore, when predicting the stability limit reset and dynamic performance of a turbocharger rotor, the effect of temperature should have a great influence on the dynamic design of the rotor.

Unbalance Response Analysis.
We applied a dummy unbalance mass of 0.01 kg · mm at node 2, and the frequency range of the steady-state synchronous response analysis was 0 to 2500 Hz (corresponding speed range 0 to 150, 000 rpm). We then measured the response amplitude of the turbine end floating ring bearing center position node 5. In Figure 21, case 1 is represented by a black curve (normal temperature), and case 2 by a red curve (high temperature). In case 1 (case 2), the curves represent 6F, 6B, and 20F modes. e frequency at the wave crest corresponding to the frequency corresponding to the fast reverse rotation (the frequency of the rotation speed within one minute) in the case 1 (case 2) Campbell diagram, as shown in Figures 22 and 23), is similar to the black (red) number in Figure 20. Because the steering mode is dominant in the response diagram, the critical speeds Ω cr1 and Ω cr2 in case 1 (case 2) 1F and 2F modes are not shown in the diagram and followed a phenomenon similar to that reported by Chouksey and Roy [29,30]. Because of the effect of temperature, the steering frequency of case 2 shifted to the left, making the change of response amplitude more obvious. e response value at the peak of the 6F and 6B modes decreased from 0.147 mm to 0.014 mm, and that at the peak of the 20F and 20B modes decreased from 0.036 mm to 0.0037 mm.

Verification of Rotor FEM.
Natural frequency analysis is very important as the basis of a dynamic analysis. We obtained the first two natural frequencies of the rotor by hammering experiments and verified them by numerical simulation. Figure 24 shows the measured parts (rotor system), force hammer, and piezoelectric acceleration sensor. To reduce the measurement error rate, we adopted a sensor with a nonlinearity of less than 1.2% and a sensitivity of 10 mv/g and placed the sensor in the middle of the rotor shaft.  A more favorable verification was the oil film data shown in Figure 6. e oil film data were experimental data provided by the turbocharger company. In the simulation, the oil film was the key factor in supporting or restraining the rotor. is also illustrates the reliability of the simulation method.

Conclusions
is study used the CHT numerical simulation method to obtain the temperature of a rotor-shaft system and found that there was a large temperature gradient when the system was in use.
Splitting the rotation shaft solved the problem that a solid model cannot obtain an axis orbit on the Ansys Workbench software platform.
Analysing the Campbell diagram shows that the temperature reduced the stiffness of the rotor structure, causing the intrinsic frequency to decrease, especially for higherorder frequencies. e decrease in stiffness was greater and the critical speed (Ωcr1) in the 1F mode did not change much, whereas the critical speed (Ωcr2) in the 2F mode decreased a large amount. erefore, dynamic design and stability prediction of a turbocharger rotor must consider the temperature during its operation.
For the analysis of an unbalanced response, the response corresponding to the critical speed in the 1F and 2F modes does not appear in Figure 21 because these two modes have higher damping and because the two steering modes were dominant. However, the response diagram can well reflect the rotor steering mode. Also, the temperature effect on the response was mainly reflected in the response amplitude, and the response frequency was slightly shifted to the left (reduced), which was similar to the decreased intrinsic frequency relation.