# Dynamic Behaviour Analysis of Turbocharger Rotor-Shaft System in Thermal Environment Based on Finite Element Method

Authors: Zhihao Liu; Renren Wang; Fang Cao; Pidong Shi
Journal: Shock and Vibration
 (2020)

---

## Abstract

The 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. The temperature field will increase the tangential force generated by the internal damping of the composite material. The 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. The analysis was done according to the modal damping coefficient, stability limit speed, and unbalance response. The results show that accurate prediction of internal damping energy dissipation in a temperature field is crucial for accurate prediction of rotor dynamic performance. This is an important step to understand dynamic rotor stress and rotor dynamic design.


---

## Body

## 1. Introduction

Turbochargers are mechanical devices that can improve fuel efficiency and reduce greenhouse gas emissions. The core component of a turbocharger is a rotor composed mainly of a turbine and a compressor. The 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. This introduces high requirements for the stable rotation of the rotor. Therefore, 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. They 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. They 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. They 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. The 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. The 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].The 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 rotor-bearing finite element model. The 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. The 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 high-temperature 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. Their 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. The 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. The study found that, as the temperature increased, the modal frequency of the beam structure decreased accordingly, and the damping ratio increased accordingly. That 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.Therefore, 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 rotor-shaft system. The temperature field was coupled with the finite element model of the rotor, compared with the instance without considering the temperature field. The rotor finite element is verified by experiment.


## 2. Numerical Model and Research Methods

### 2.1. Fluid and Heat Transfer Sections

Before the thermal modal analysis, the aerodynamic thermal analysis was performed. The current industry standard modelling approaches assume the turbine and compressor operate under adiabatic conditions [15]. The CHT simulation method [16–22] is used to obtain the temperature distribution of the rotor, and the node temperature is provided for the modal simulation part.


#### 2.1.1. Turbocharge and Compressor Geometry

The turbocharger had an impeller with 10 blades, and the compressor had an impeller with 6 blades and 6 splitters. The turbocharger and compressor parameters are shown in Table1. To improve the accuracy of the calculations, we modelled both the turbocharger and the compressor using full passages. The solid parts and the air flow are shown in Figure 1.Table 1
Turbocharger and compressor parameters.
Turbine sideCompressor sideParametersValue and unitsParametersValue and unitsBlades number10Blades number6 + 6 (−)Impeller inlet diameter55.05 mmImpeller outlet diameter56.5 mmTip clearance0.41 mmTip clearance0.26 mmFigure 1
Centrifugal turbocharger and compressor with a solid casing.


#### 2.1.2. 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. This allows the noninterpolated exchange of heat flux between adjacent grids [23]. The calculation accuracy of CHT is very sensitive to the resolution of the fluid boundary layer grid. Therefore, 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. The 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. The impeller grid and the fluid grid were connected in a general grid interface mode in Ansys CFX software.Figure 2
Global grid for conjugate heat transfer calculation.Figure 3
Computational grids for conjugate heat transfer calculation.


#### 2.1.3. 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 Table2.Table 2
The boundary conditions for aerothermal analysis.
Turbine sideCompressor sideMedium (intensity = 5%)Medium (intensity = 5%)Inlet mass flow0.074 kg/sInlet total pressure99935.9 PaInlet total temperature872.97 KInlet total temperature297.201 KOutlet static pressure96892.7 PaOutlet static pressure127578 PaThe 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].


#### 2.1.4. 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. The 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(1)∂ρf∂t+∇⋅ρfuf=0,(2)∂ufρf∂t+∇⋅ρfufuf=−∇p+∇⋅τ,(3)∂htotρf∂t+∇⋅ρfufhtot=∇⋅λf∇T+∇⋅uf⋅τ,where uf, ρf, p, τ, λf, T, htot, 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. The energy equation is(4)∂htotρs∂t+∇⋅ρsushtot=∇⋅λs∇T+SE,where us is the velocity vector of the solid, us=uf; SE is the optional volume heat source, SE=∇⋅us⋅τ; ρs is the density of the solid; and λs is the thermal conductivity of the solid.This study did not directly solve (1)–(4). Instead, they were converted to the steady-state Reynolds average Navier–Stokes method to calculate turbulence. The 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]. That 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].


### 2.2. Modal Part

#### 2.2.1. Rotor-Shaft System Geometry

The model of the rotor-shaft system was provided by the turbocharger company and agreed with the real one. The 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 Figure4 shows. The turbine wheel and shaft were connected by friction welding, whereas the compressor was axially and circumferentially fixed by a left-hand nut.Figure 4
Composition of the rotor-shaft system.The materials of the main parts of the rotor-shaft system are shown in Table3, based on the actual situation.Table 3
Material of the main part of the rotor-shaft system.
Part nameMaterial nameTurbo impellerK418Compressor impeller and nutZL105Shaft42CrMoThrust bearing and seal sleeve316 stainless steelThis study investigated the effect of temperature and found that the material properties changed as the temperature changed. The function curves of the material properties of each part as a function of temperature are shown in Figure5 [25, 26].Figure 5
Material properties of rotor-shaft system as a function of temperature (Young’s modulus, Poisson ratio, speciﬁc heat capacity, and thermal conductivity). (a) K418. (b) ZL105. (c) 42CrMo. (d) 316 stainless steel [25, 26].
(a)(b)(c)(d)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.Table4 shows what we assumed was the material damping coefficient of each part of the rotor.Table 4
Damping coefficient of each part.
Part nameDamping coefficientTurbo impeller0.002Compressor impeller and nut0.001Shaft0.005Thrust bearing and seal sleeve0.005


#### 2.2.2. Floating Ring Bearing Parameters

The rotor shaft was supported by a floating ring bearing whose detailed parameters are shown in Table5.Table 5
Parameters of floating ring bearing that supported the rotor shaft.
Parameter nameParameter valueFloating ring mass5.17 gInner length5.14 mmOuter length7.96 mmInner diameter7.50 mmOuter diameter12.9 mmTo 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.Figure 6
Total impedance of floating ring bearing (equivalent stiffness and damping coefficients). (a) Turbine side. (b) Compressor side.
(a)(b)


#### 2.2.3. Finite Element Model of Rotor

This study used a large-scale 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 Figure7. The blue lines (solid and dotted) in the figure represent the division positions.Figure 7
Rotor shaft divided into four parts and 1/4 cylinder.The solid model was discretised, and an Ansys Solid186 element was used. The Solid186 element is a high-order, 3D 20-node element with quadratic displacement characteristics. The element is defined by 20 nodes, each node with three degrees of freedom (translation of each node in thex-, y-, and z-directions). The 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 “ith” element after segmentation.Figure 8
Solid186 element generated by the sweep method.We transformed (deformed) the original hexahedron element to an approximately one-quarter cylinder element. In Figure8, the red dotted line represents the actual edge of the element, and the red solid line represents the edge of the element. The 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.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. Figure9 shows the (i+1), (i+2), and (i+3) elements generated by three methods.Figure 9
Solid186 elements generated by tetrahedral, pyramid, and prism methods.The 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. The rotor model was discretised intoNe elements, and the total number of nodes became a. The node displacement vector of the ith element is(5)qe=xiI,yiI,ziI,xiJ,yiJ,ziJ,…,xiμ,yiμ,ziμ,…,xiA,yiA,ziA,xiB,yiB,ziBT,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. These 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(6)Mq¨t+Cq˙t+Kqt=ft,where M=Mtrs+Mrot; C=−ΩG+Ccon+Cbrg; and K=Kbrg+Ksh+Ktem+Ktemσ.Here,M is a mass matrix that includes mass components caused by translation along the axis (Mtrs) and rotation along the axis (Mrot). C is a damping matrix (asymmetric matrix coefficient of the velocity vector) that includes a skew-symmetric gyro matrix (G), a constant structural damping matrix (Ccon=βKcon) with hysteresis characteristics caused by internal friction of the material and damping caused by the floating ring bearing support matrix (Cbrg). K is the stiffness matrix (matrix coefficient of the displacement vector) that includes the stiffness matrix (Kbrg) caused by the support of the floating ring bearing, the spin-softening bending matrix (Ksh) generated by the rotor rotation, the stiffness matrix (Ktem) of the elastic modulus, the stiffness matrix (Ktemσ) caused by the nonuniform thermal stress, and thermal deformation caused by temperature changes. The size of all matrices is 3a×3a. Ω is the rotor speed in rpm. β is the material damping (Rayleigh damping), stiffness damping coefficient.The total displacement and global force vectorsqt and ft are expressed as(7)qt=xta×1yta×1zta×13a×1=x1tx2tx3t,…,xa−1txatTy1ty2ty3t,…,ya−1tyatTz1tz2tz3t,…,za−1tzatT,ft=fxta×1fyta×1fzta×13a×1=fx1tfx2tfx3t,…,fxa−1tfxatTfy1tfy2tfy3t,…,fya−1tfyatTfz1tfz2tfz3t,…,fza−1tfzatT.Writing (6) in state-space form yields(8)Aw˙t=Bwt+Ft,where matrices A and B; state vector wt; and force vector Ft are(9)A=0MMC6a×6aB=M00−K6a×6awt=q˙tqt6a×1Ft=0ft6a×1,where 0 is an empty matrix of size 3a×3a and an empty vector of size 3a×1.Figure10 shows a simplified finite element model of the rotor-shaft system. The grey-blue points in the figure are used to observe the modal shapes, line trajectories, and unbalanced harmonic response analysis.Figure 10
Finite element model of rotor-shaft system.Figure11 shows the rotor finite element global grid. The 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.Figure 11
Global grid of rotor-shaft system.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.


## 2.1. Fluid and Heat Transfer Sections

Before the thermal modal analysis, the aerodynamic thermal analysis was performed. The current industry standard modelling approaches assume the turbine and compressor operate under adiabatic conditions [15]. The CHT simulation method [16–22] is used to obtain the temperature distribution of the rotor, and the node temperature is provided for the modal simulation part.


### 2.1.1. Turbocharge and Compressor Geometry

The turbocharger had an impeller with 10 blades, and the compressor had an impeller with 6 blades and 6 splitters. The turbocharger and compressor parameters are shown in Table1. To improve the accuracy of the calculations, we modelled both the turbocharger and the compressor using full passages. The solid parts and the air flow are shown in Figure 1.Table 1
Turbocharger and compressor parameters.
Turbine sideCompressor sideParametersValue and unitsParametersValue and unitsBlades number10Blades number6 + 6 (−)Impeller inlet diameter55.05 mmImpeller outlet diameter56.5 mmTip clearance0.41 mmTip clearance0.26 mmFigure 1
Centrifugal turbocharger and compressor with a solid casing.


### 2.1.2. 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. This allows the noninterpolated exchange of heat flux between adjacent grids [23]. The calculation accuracy of CHT is very sensitive to the resolution of the fluid boundary layer grid. Therefore, 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. The 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. The impeller grid and the fluid grid were connected in a general grid interface mode in Ansys CFX software.Figure 2
Global grid for conjugate heat transfer calculation.Figure 3
Computational grids for conjugate heat transfer calculation.


### 2.1.3. 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 Table2.Table 2
The boundary conditions for aerothermal analysis.
Turbine sideCompressor sideMedium (intensity = 5%)Medium (intensity = 5%)Inlet mass flow0.074 kg/sInlet total pressure99935.9 PaInlet total temperature872.97 KInlet total temperature297.201 KOutlet static pressure96892.7 PaOutlet static pressure127578 PaThe 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].


### 2.1.4. 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. The 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(1)∂ρf∂t+∇⋅ρfuf=0,(2)∂ufρf∂t+∇⋅ρfufuf=−∇p+∇⋅τ,(3)∂htotρf∂t+∇⋅ρfufhtot=∇⋅λf∇T+∇⋅uf⋅τ,where uf, ρf, p, τ, λf, T, htot, 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. The energy equation is(4)∂htotρs∂t+∇⋅ρsushtot=∇⋅λs∇T+SE,where us is the velocity vector of the solid, us=uf; SE is the optional volume heat source, SE=∇⋅us⋅τ; ρs is the density of the solid; and λs is the thermal conductivity of the solid.This study did not directly solve (1)–(4). Instead, they were converted to the steady-state Reynolds average Navier–Stokes method to calculate turbulence. The 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]. That 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].


## 2.1.1. Turbocharge and Compressor Geometry

The turbocharger had an impeller with 10 blades, and the compressor had an impeller with 6 blades and 6 splitters. The turbocharger and compressor parameters are shown in Table1. To improve the accuracy of the calculations, we modelled both the turbocharger and the compressor using full passages. The solid parts and the air flow are shown in Figure 1.Table 1
Turbocharger and compressor parameters.
Turbine sideCompressor sideParametersValue and unitsParametersValue and unitsBlades number10Blades number6 + 6 (−)Impeller inlet diameter55.05 mmImpeller outlet diameter56.5 mmTip clearance0.41 mmTip clearance0.26 mmFigure 1
Centrifugal turbocharger and compressor with a solid casing.


## 2.1.2. 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. This allows the noninterpolated exchange of heat flux between adjacent grids [23]. The calculation accuracy of CHT is very sensitive to the resolution of the fluid boundary layer grid. Therefore, 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. The 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. The impeller grid and the fluid grid were connected in a general grid interface mode in Ansys CFX software.Figure 2
Global grid for conjugate heat transfer calculation.Figure 3
Computational grids for conjugate heat transfer calculation.


## 2.1.3. 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 Table2.Table 2
The boundary conditions for aerothermal analysis.
Turbine sideCompressor sideMedium (intensity = 5%)Medium (intensity = 5%)Inlet mass flow0.074 kg/sInlet total pressure99935.9 PaInlet total temperature872.97 KInlet total temperature297.201 KOutlet static pressure96892.7 PaOutlet static pressure127578 PaThe 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].


## 2.1.4. 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. The 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(1)∂ρf∂t+∇⋅ρfuf=0,(2)∂ufρf∂t+∇⋅ρfufuf=−∇p+∇⋅τ,(3)∂htotρf∂t+∇⋅ρfufhtot=∇⋅λf∇T+∇⋅uf⋅τ,where uf, ρf, p, τ, λf, T, htot, 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. The energy equation is(4)∂htotρs∂t+∇⋅ρsushtot=∇⋅λs∇T+SE,where us is the velocity vector of the solid, us=uf; SE is the optional volume heat source, SE=∇⋅us⋅τ; ρs is the density of the solid; and λs is the thermal conductivity of the solid.This study did not directly solve (1)–(4). Instead, they were converted to the steady-state Reynolds average Navier–Stokes method to calculate turbulence. The 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]. That 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].


## 2.2. Modal Part

### 2.2.1. Rotor-Shaft System Geometry

The model of the rotor-shaft system was provided by the turbocharger company and agreed with the real one. The 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 Figure4 shows. The turbine wheel and shaft were connected by friction welding, whereas the compressor was axially and circumferentially fixed by a left-hand nut.Figure 4
Composition of the rotor-shaft system.The materials of the main parts of the rotor-shaft system are shown in Table3, based on the actual situation.Table 3
Material of the main part of the rotor-shaft system.
Part nameMaterial nameTurbo impellerK418Compressor impeller and nutZL105Shaft42CrMoThrust bearing and seal sleeve316 stainless steelThis study investigated the effect of temperature and found that the material properties changed as the temperature changed. The function curves of the material properties of each part as a function of temperature are shown in Figure5 [25, 26].Figure 5
Material properties of rotor-shaft system as a function of temperature (Young’s modulus, Poisson ratio, speciﬁc heat capacity, and thermal conductivity). (a) K418. (b) ZL105. (c) 42CrMo. (d) 316 stainless steel [25, 26].
(a)(b)(c)(d)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.Table4 shows what we assumed was the material damping coefficient of each part of the rotor.Table 4
Damping coefficient of each part.
Part nameDamping coefficientTurbo impeller0.002Compressor impeller and nut0.001Shaft0.005Thrust bearing and seal sleeve0.005


### 2.2.2. Floating Ring Bearing Parameters

The rotor shaft was supported by a floating ring bearing whose detailed parameters are shown in Table5.Table 5
Parameters of floating ring bearing that supported the rotor shaft.
Parameter nameParameter valueFloating ring mass5.17 gInner length5.14 mmOuter length7.96 mmInner diameter7.50 mmOuter diameter12.9 mmTo 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.Figure 6
Total impedance of floating ring bearing (equivalent stiffness and damping coefficients). (a) Turbine side. (b) Compressor side.
(a)(b)


### 2.2.3. Finite Element Model of Rotor

This study used a large-scale 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 Figure7. The blue lines (solid and dotted) in the figure represent the division positions.Figure 7
Rotor shaft divided into four parts and 1/4 cylinder.The solid model was discretised, and an Ansys Solid186 element was used. The Solid186 element is a high-order, 3D 20-node element with quadratic displacement characteristics. The element is defined by 20 nodes, each node with three degrees of freedom (translation of each node in thex-, y-, and z-directions). The 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 “ith” element after segmentation.Figure 8
Solid186 element generated by the sweep method.We transformed (deformed) the original hexahedron element to an approximately one-quarter cylinder element. In Figure8, the red dotted line represents the actual edge of the element, and the red solid line represents the edge of the element. The 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.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. Figure9 shows the (i+1), (i+2), and (i+3) elements generated by three methods.Figure 9
Solid186 elements generated by tetrahedral, pyramid, and prism methods.The 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. The rotor model was discretised intoNe elements, and the total number of nodes became a. The node displacement vector of the ith element is(5)qe=xiI,yiI,ziI,xiJ,yiJ,ziJ,…,xiμ,yiμ,ziμ,…,xiA,yiA,ziA,xiB,yiB,ziBT,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. These 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(6)Mq¨t+Cq˙t+Kqt=ft,where M=Mtrs+Mrot; C=−ΩG+Ccon+Cbrg; and K=Kbrg+Ksh+Ktem+Ktemσ.Here,M is a mass matrix that includes mass components caused by translation along the axis (Mtrs) and rotation along the axis (Mrot). C is a damping matrix (asymmetric matrix coefficient of the velocity vector) that includes a skew-symmetric gyro matrix (G), a constant structural damping matrix (Ccon=βKcon) with hysteresis characteristics caused by internal friction of the material and damping caused by the floating ring bearing support matrix (Cbrg). K is the stiffness matrix (matrix coefficient of the displacement vector) that includes the stiffness matrix (Kbrg) caused by the support of the floating ring bearing, the spin-softening bending matrix (Ksh) generated by the rotor rotation, the stiffness matrix (Ktem) of the elastic modulus, the stiffness matrix (Ktemσ) caused by the nonuniform thermal stress, and thermal deformation caused by temperature changes. The size of all matrices is 3a×3a. Ω is the rotor speed in rpm. β is the material damping (Rayleigh damping), stiffness damping coefficient.The total displacement and global force vectorsqt and ft are expressed as(7)qt=xta×1yta×1zta×13a×1=x1tx2tx3t,…,xa−1txatTy1ty2ty3t,…,ya−1tyatTz1tz2tz3t,…,za−1tzatT,ft=fxta×1fyta×1fzta×13a×1=fx1tfx2tfx3t,…,fxa−1tfxatTfy1tfy2tfy3t,…,fya−1tfyatTfz1tfz2tfz3t,…,fza−1tfzatT.Writing (6) in state-space form yields(8)Aw˙t=Bwt+Ft,where matrices A and B; state vector wt; and force vector Ft are(9)A=0MMC6a×6aB=M00−K6a×6awt=q˙tqt6a×1Ft=0ft6a×1,where 0 is an empty matrix of size 3a×3a and an empty vector of size 3a×1.Figure10 shows a simplified finite element model of the rotor-shaft system. The grey-blue points in the figure are used to observe the modal shapes, line trajectories, and unbalanced harmonic response analysis.Figure 10
Finite element model of rotor-shaft system.Figure11 shows the rotor finite element global grid. The 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.Figure 11
Global grid of rotor-shaft system.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.


## 2.2.1. Rotor-Shaft System Geometry

The model of the rotor-shaft system was provided by the turbocharger company and agreed with the real one. The 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 Figure4 shows. The turbine wheel and shaft were connected by friction welding, whereas the compressor was axially and circumferentially fixed by a left-hand nut.Figure 4
Composition of the rotor-shaft system.The materials of the main parts of the rotor-shaft system are shown in Table3, based on the actual situation.Table 3
Material of the main part of the rotor-shaft system.
Part nameMaterial nameTurbo impellerK418Compressor impeller and nutZL105Shaft42CrMoThrust bearing and seal sleeve316 stainless steelThis study investigated the effect of temperature and found that the material properties changed as the temperature changed. The function curves of the material properties of each part as a function of temperature are shown in Figure5 [25, 26].Figure 5
Material properties of rotor-shaft system as a function of temperature (Young’s modulus, Poisson ratio, speciﬁc heat capacity, and thermal conductivity). (a) K418. (b) ZL105. (c) 42CrMo. (d) 316 stainless steel [25, 26].
(a)(b)(c)(d)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.Table4 shows what we assumed was the material damping coefficient of each part of the rotor.Table 4
Damping coefficient of each part.
Part nameDamping coefficientTurbo impeller0.002Compressor impeller and nut0.001Shaft0.005Thrust bearing and seal sleeve0.005


## 2.2.2. Floating Ring Bearing Parameters

The rotor shaft was supported by a floating ring bearing whose detailed parameters are shown in Table5.Table 5
Parameters of floating ring bearing that supported the rotor shaft.
Parameter nameParameter valueFloating ring mass5.17 gInner length5.14 mmOuter length7.96 mmInner diameter7.50 mmOuter diameter12.9 mmTo 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.Figure 6
Total impedance of floating ring bearing (equivalent stiffness and damping coefficients). (a) Turbine side. (b) Compressor side.
(a)(b)


## 2.2.3. Finite Element Model of Rotor

This study used a large-scale 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 Figure7. The blue lines (solid and dotted) in the figure represent the division positions.Figure 7
Rotor shaft divided into four parts and 1/4 cylinder.The solid model was discretised, and an Ansys Solid186 element was used. The Solid186 element is a high-order, 3D 20-node element with quadratic displacement characteristics. The element is defined by 20 nodes, each node with three degrees of freedom (translation of each node in thex-, y-, and z-directions). The 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 “ith” element after segmentation.Figure 8
Solid186 element generated by the sweep method.We transformed (deformed) the original hexahedron element to an approximately one-quarter cylinder element. In Figure8, the red dotted line represents the actual edge of the element, and the red solid line represents the edge of the element. The 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.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. Figure9 shows the (i+1), (i+2), and (i+3) elements generated by three methods.Figure 9
Solid186 elements generated by tetrahedral, pyramid, and prism methods.The 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. The rotor model was discretised intoNe elements, and the total number of nodes became a. The node displacement vector of the ith element is(5)qe=xiI,yiI,ziI,xiJ,yiJ,ziJ,…,xiμ,yiμ,ziμ,…,xiA,yiA,ziA,xiB,yiB,ziBT,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. These 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(6)Mq¨t+Cq˙t+Kqt=ft,where M=Mtrs+Mrot; C=−ΩG+Ccon+Cbrg; and K=Kbrg+Ksh+Ktem+Ktemσ.Here,M is a mass matrix that includes mass components caused by translation along the axis (Mtrs) and rotation along the axis (Mrot). C is a damping matrix (asymmetric matrix coefficient of the velocity vector) that includes a skew-symmetric gyro matrix (G), a constant structural damping matrix (Ccon=βKcon) with hysteresis characteristics caused by internal friction of the material and damping caused by the floating ring bearing support matrix (Cbrg). K is the stiffness matrix (matrix coefficient of the displacement vector) that includes the stiffness matrix (Kbrg) caused by the support of the floating ring bearing, the spin-softening bending matrix (Ksh) generated by the rotor rotation, the stiffness matrix (Ktem) of the elastic modulus, the stiffness matrix (Ktemσ) caused by the nonuniform thermal stress, and thermal deformation caused by temperature changes. The size of all matrices is 3a×3a. Ω is the rotor speed in rpm. β is the material damping (Rayleigh damping), stiffness damping coefficient.The total displacement and global force vectorsqt and ft are expressed as(7)qt=xta×1yta×1zta×13a×1=x1tx2tx3t,…,xa−1txatTy1ty2ty3t,…,ya−1tyatTz1tz2tz3t,…,za−1tzatT,ft=fxta×1fyta×1fzta×13a×1=fx1tfx2tfx3t,…,fxa−1tfxatTfy1tfy2tfy3t,…,fya−1tfyatTfz1tfz2tfz3t,…,fza−1tfzatT.Writing (6) in state-space form yields(8)Aw˙t=Bwt+Ft,where matrices A and B; state vector wt; and force vector Ft are(9)A=0MMC6a×6aB=M00−K6a×6awt=q˙tqt6a×1Ft=0ft6a×1,where 0 is an empty matrix of size 3a×3a and an empty vector of size 3a×1.Figure10 shows a simplified finite element model of the rotor-shaft system. The grey-blue points in the figure are used to observe the modal shapes, line trajectories, and unbalanced harmonic response analysis.Figure 10
Finite element model of rotor-shaft system.Figure11 shows the rotor finite element global grid. The 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.Figure 11
Global grid of rotor-shaft system.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.


## 3. Results and Discussion

### 3.1. Heat Transfer Results and Experimental Verification

After the calculation, we verified the numerical results of the mass flow and the supercharger ratio (expansion ratio). The numerical results agreed with the experimental data, and the error control was approximately 7%, as Figure12 shows.Figure 12
Numerical simulation verified by experiments. (a) Turbine side. (b) Compressor side.
(a)(b)The temperature distribution of the rotor was obtained through postprocessing, as Figure13 shows; the temperature distribution law was the same as that in the Bohn [19] study result and is shown in Figure 14.Figure 13
Postprocessing rotor temperature distribution.Figure 14
Rotor temperature distribution of Bohn [19].Figures13 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.


### 3.2. Modal Results

#### 3.2.1. Campbell Diagram and Stability Diagram of Normal Temperature Environment

Figure15 shows the Campbell diagram of the rotor-shaft system in case 1, that is, without considering the temperature. The figure was obtained from the whirl frequency (obtained from the imaginary part of the eigenvalue) and the rotational speed. To present a clear figure and express its basic characteristics, only the first four whirling directions are described. In the figure, there are two forward whirls (the whirl direction of the rotor agrees with the rotation direction of the speed), which are marked as 1F and 2F in order of increased frequency. There are also two backward whirls (the whirl direction of the rotor is opposite to the rotation direction of the speed) that are marked as 3B and 4B, respectively. The synchronous whirl line is called the synchronous excitation line. In the figure, the velocities corresponding to the intersection of the synchronous whirl line and the frequency line of the forward whirl are marked as Ωcr1 and Ωcr2, respectively, and were calculated as follows: Ωcr1=5,011.5rpm and Ωcr2=43,622rpm. Because the backward whirl could not be excited by the unbalanced force, it is not marked or considered.Figure 15
Campbell diagram showing whirl frequency and rotor-spin speed for the first four modes.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. That 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.Figure16 shows that the curve of the damping ratio of the four modes in case 1 changed with the rotation speed. The 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. The point where the damping ratio curve intersects with the zero line is the stable limit speed (SLS). The figure shows that the rotor-shaft system became unstable when SLS1=56,550rpm in 1F mode and then became unstable when SLS2=97,946rpm in 2F mode.Figure 16
Stability diagram of rotor-shaft system without considering temperature.Figure17 shows the four modes close to the SLS 1 stability limit speed at 55,000 rpm. The 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. The 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. The 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 17
Mode shapes at 55,000 rpm (SLS 1). (a) Mode 1F. (b) Mode 2F. (c) Mode 3B. (d) Mode 4B.
(a)(b)(c)(d)Figure18 shows the four modes close to the SLS 2 stability limit speed at 97,946 rpm. The rotor was close to the unstable state at this time. The 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.Figure 18
Mode shapes at 97,000 rpm (SLS 2). (a) Mode 1F. (b) Mode 2F. (c) Mode 3B. (d) Mode 4B.
(a)(b)(c)(d)


#### 3.2.2. Effect of Temperature on the Campbell Diagram and Stability Diagram

We exported the temperature data of each node in Figure12 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.Figure19 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. This decrease was also caused by the intrinsic frequency drop.Figure 19
Effect of temperature on the Campbell diagram.Figure20 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. The 2F mode change was more obvious, and the same SLS 2 decreased from the original 97,946 rpm to 95,969 rpm. The reason for this can be explained by the following formula: ξr=α/2ωhr+βωhr/2. Here, ξr is the damping ratio coefficient of the rth mode, α is the mass damping coefficient, β is the stiffness damping coefficient, and ωhr is the natural frequency of the rth 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. Therefore, 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.Figure 20
Effect of temperature on stability diagram.


#### 3.2.3. Unbalance Response Analysis

We applied a dummy unbalance mass of 0.01kg⋅mm at node 2, and the frequency range of the steady-state synchronous response analysis was 0to2500Hz (corresponding speed range 0to150,000rpm). 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. The 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. The 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.Figure 21
Response comparison between case 1 and case 2 of the turbine side floating ring bearing center position (node 5).Figure 22
Partial enlargement of Campbell diagram (case 1).Figure 23
Partial enlargement of Campbell diagram (case 2).


#### 3.2.4. 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. Figure24 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. Table 6 compares the natural frequency experimental value with the simulated value. The experimental and simulated values of the first natural frequency were 1386 Hz and 1453 Hz, respectively. The second-order natural frequencies were 4100 Hz and 4383 Hz, respectively, and the error rate was within 7%. But, the higher mass of the accelerometer is the main cause of error value between the experimental value and the simulated value.Figure 24
Hammering experiment for a rotor.Table 6
Experimental and calculated values of natural frequencies.
Order numberExperimental valueSimulated valueError rate11386 Hz1453 Hz0.04624100 Hz4383 Hz0.064A more favorable verification was the oil film data shown in Figure6. The 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. This also illustrates the reliability of the simulation method.


## 3.1. Heat Transfer Results and Experimental Verification

After the calculation, we verified the numerical results of the mass flow and the supercharger ratio (expansion ratio). The numerical results agreed with the experimental data, and the error control was approximately 7%, as Figure12 shows.Figure 12
Numerical simulation verified by experiments. (a) Turbine side. (b) Compressor side.
(a)(b)The temperature distribution of the rotor was obtained through postprocessing, as Figure13 shows; the temperature distribution law was the same as that in the Bohn [19] study result and is shown in Figure 14.Figure 13
Postprocessing rotor temperature distribution.Figure 14
Rotor temperature distribution of Bohn [19].Figures13 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.


## 3.2. Modal Results

### 3.2.1. Campbell Diagram and Stability Diagram of Normal Temperature Environment

Figure15 shows the Campbell diagram of the rotor-shaft system in case 1, that is, without considering the temperature. The figure was obtained from the whirl frequency (obtained from the imaginary part of the eigenvalue) and the rotational speed. To present a clear figure and express its basic characteristics, only the first four whirling directions are described. In the figure, there are two forward whirls (the whirl direction of the rotor agrees with the rotation direction of the speed), which are marked as 1F and 2F in order of increased frequency. There are also two backward whirls (the whirl direction of the rotor is opposite to the rotation direction of the speed) that are marked as 3B and 4B, respectively. The synchronous whirl line is called the synchronous excitation line. In the figure, the velocities corresponding to the intersection of the synchronous whirl line and the frequency line of the forward whirl are marked as Ωcr1 and Ωcr2, respectively, and were calculated as follows: Ωcr1=5,011.5rpm and Ωcr2=43,622rpm. Because the backward whirl could not be excited by the unbalanced force, it is not marked or considered.Figure 15
Campbell diagram showing whirl frequency and rotor-spin speed for the first four modes.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. That 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.Figure16 shows that the curve of the damping ratio of the four modes in case 1 changed with the rotation speed. The 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. The point where the damping ratio curve intersects with the zero line is the stable limit speed (SLS). The figure shows that the rotor-shaft system became unstable when SLS1=56,550rpm in 1F mode and then became unstable when SLS2=97,946rpm in 2F mode.Figure 16
Stability diagram of rotor-shaft system without considering temperature.Figure17 shows the four modes close to the SLS 1 stability limit speed at 55,000 rpm. The 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. The 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. The 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 17
Mode shapes at 55,000 rpm (SLS 1). (a) Mode 1F. (b) Mode 2F. (c) Mode 3B. (d) Mode 4B.
(a)(b)(c)(d)Figure18 shows the four modes close to the SLS 2 stability limit speed at 97,946 rpm. The rotor was close to the unstable state at this time. The 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.Figure 18
Mode shapes at 97,000 rpm (SLS 2). (a) Mode 1F. (b) Mode 2F. (c) Mode 3B. (d) Mode 4B.
(a)(b)(c)(d)


### 3.2.2. Effect of Temperature on the Campbell Diagram and Stability Diagram

We exported the temperature data of each node in Figure12 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.Figure19 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. This decrease was also caused by the intrinsic frequency drop.Figure 19
Effect of temperature on the Campbell diagram.Figure20 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. The 2F mode change was more obvious, and the same SLS 2 decreased from the original 97,946 rpm to 95,969 rpm. The reason for this can be explained by the following formula: ξr=α/2ωhr+βωhr/2. Here, ξr is the damping ratio coefficient of the rth mode, α is the mass damping coefficient, β is the stiffness damping coefficient, and ωhr is the natural frequency of the rth 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. Therefore, 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.Figure 20
Effect of temperature on stability diagram.


### 3.2.3. Unbalance Response Analysis

We applied a dummy unbalance mass of 0.01kg⋅mm at node 2, and the frequency range of the steady-state synchronous response analysis was 0to2500Hz (corresponding speed range 0to150,000rpm). 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. The 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. The 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.Figure 21
Response comparison between case 1 and case 2 of the turbine side floating ring bearing center position (node 5).Figure 22
Partial enlargement of Campbell diagram (case 1).Figure 23
Partial enlargement of Campbell diagram (case 2).


### 3.2.4. 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. Figure24 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. Table 6 compares the natural frequency experimental value with the simulated value. The experimental and simulated values of the first natural frequency were 1386 Hz and 1453 Hz, respectively. The second-order natural frequencies were 4100 Hz and 4383 Hz, respectively, and the error rate was within 7%. But, the higher mass of the accelerometer is the main cause of error value between the experimental value and the simulated value.Figure 24
Hammering experiment for a rotor.Table 6
Experimental and calculated values of natural frequencies.
Order numberExperimental valueSimulated valueError rate11386 Hz1453 Hz0.04624100 Hz4383 Hz0.064A more favorable verification was the oil film data shown in Figure6. The 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. This also illustrates the reliability of the simulation method.


## 3.2.1. Campbell Diagram and Stability Diagram of Normal Temperature Environment

Figure15 shows the Campbell diagram of the rotor-shaft system in case 1, that is, without considering the temperature. The figure was obtained from the whirl frequency (obtained from the imaginary part of the eigenvalue) and the rotational speed. To present a clear figure and express its basic characteristics, only the first four whirling directions are described. In the figure, there are two forward whirls (the whirl direction of the rotor agrees with the rotation direction of the speed), which are marked as 1F and 2F in order of increased frequency. There are also two backward whirls (the whirl direction of the rotor is opposite to the rotation direction of the speed) that are marked as 3B and 4B, respectively. The synchronous whirl line is called the synchronous excitation line. In the figure, the velocities corresponding to the intersection of the synchronous whirl line and the frequency line of the forward whirl are marked as Ωcr1 and Ωcr2, respectively, and were calculated as follows: Ωcr1=5,011.5rpm and Ωcr2=43,622rpm. Because the backward whirl could not be excited by the unbalanced force, it is not marked or considered.Figure 15
Campbell diagram showing whirl frequency and rotor-spin speed for the first four modes.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. That 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.Figure16 shows that the curve of the damping ratio of the four modes in case 1 changed with the rotation speed. The 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. The point where the damping ratio curve intersects with the zero line is the stable limit speed (SLS). The figure shows that the rotor-shaft system became unstable when SLS1=56,550rpm in 1F mode and then became unstable when SLS2=97,946rpm in 2F mode.Figure 16
Stability diagram of rotor-shaft system without considering temperature.Figure17 shows the four modes close to the SLS 1 stability limit speed at 55,000 rpm. The 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. The 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. The 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 17
Mode shapes at 55,000 rpm (SLS 1). (a) Mode 1F. (b) Mode 2F. (c) Mode 3B. (d) Mode 4B.
(a)(b)(c)(d)Figure18 shows the four modes close to the SLS 2 stability limit speed at 97,946 rpm. The rotor was close to the unstable state at this time. The 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.Figure 18
Mode shapes at 97,000 rpm (SLS 2). (a) Mode 1F. (b) Mode 2F. (c) Mode 3B. (d) Mode 4B.
(a)(b)(c)(d)


## 3.2.2. Effect of Temperature on the Campbell Diagram and Stability Diagram

We exported the temperature data of each node in Figure12 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.Figure19 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. This decrease was also caused by the intrinsic frequency drop.Figure 19
Effect of temperature on the Campbell diagram.Figure20 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. The 2F mode change was more obvious, and the same SLS 2 decreased from the original 97,946 rpm to 95,969 rpm. The reason for this can be explained by the following formula: ξr=α/2ωhr+βωhr/2. Here, ξr is the damping ratio coefficient of the rth mode, α is the mass damping coefficient, β is the stiffness damping coefficient, and ωhr is the natural frequency of the rth 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. Therefore, 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.Figure 20
Effect of temperature on stability diagram.


## 3.2.3. Unbalance Response Analysis

We applied a dummy unbalance mass of 0.01kg⋅mm at node 2, and the frequency range of the steady-state synchronous response analysis was 0to2500Hz (corresponding speed range 0to150,000rpm). 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. The 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. The 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.Figure 21
Response comparison between case 1 and case 2 of the turbine side floating ring bearing center position (node 5).Figure 22
Partial enlargement of Campbell diagram (case 1).Figure 23
Partial enlargement of Campbell diagram (case 2).


## 3.2.4. 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. Figure24 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. Table 6 compares the natural frequency experimental value with the simulated value. The experimental and simulated values of the first natural frequency were 1386 Hz and 1453 Hz, respectively. The second-order natural frequencies were 4100 Hz and 4383 Hz, respectively, and the error rate was within 7%. But, the higher mass of the accelerometer is the main cause of error value between the experimental value and the simulated value.Figure 24
Hammering experiment for a rotor.Table 6
Experimental and calculated values of natural frequencies.
Order numberExperimental valueSimulated valueError rate11386 Hz1453 Hz0.04624100 Hz4383 Hz0.064A more favorable verification was the oil film data shown in Figure6. The 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. This also illustrates the reliability of the simulation method.


## 4. Conclusions

This 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 higher-order frequencies. The 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. Therefore, 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 Figure21 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.



---

*Source: 8888504-2020-08-14.xml*