Flow-thermo-elastic coupled vibration characteristics of aeroengine high pressure rotor at the typical start and stop conditions

2.1. Vibration equation of thermal-elastic coupling

According to the balance theory of elastic mechanics, thermodynamic elastic equation can be obtained from the expressions of thermal stress and thermal strain. It is expressed as:

where $G$ is shear modulus of the elasticity, $F_x$, $F_y$, $F_z$ is the component of unit volume force at the $x$, $y$ and $z$ directions, ${\nabla }^2$ is the Laplace operator, ${\nabla }^2=\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}$, $e$ is the total strain, $e={\epsilon }_x+{\epsilon }_y+{\epsilon }_z$, $\beta$ is the thermal coefficient.

Based on the interpolation principle of 2D elasticity plane, the strain component is obtained:

where $\left[B\right]$ is geometric matrix, $\{\delta {\}}^e$ is the component of node displacement.

The strain of thermoelastic body is the sum of the strain generated by external load and the strain generated by temperature change. It is:

where $\left[D\right]$ is elastic matrix, $\left\{{\epsilon }_0\right\}$ is initial strain, and the expression is:

Substituting Eq. (3) into Eq. (2), it is obtained:

Based on the principle of finite element, the structure is divided into a finite number of unit bodies. The strain energy generated by thermal stress on the element $e$ can be obtained:

where $\stackrel{-}{t}$ is the thickness of element. It is obtained by further processing:

where ${\left[K\right]}^e$ is stiffness matrix of the element, ${\left[K\right]}^e={\underset{e}{\iint }\left\{\left.B\right\}\right.}^T\left[D\right]\left\{\left.B\right\}\right.\stackrel{-}{t}dxdy$.

Let ${\left\{\left.L\right\}\right.}^e={\underset{e}{\iint }\left[B\right]}^T\left[D\right]{\left\{\left.{\epsilon }_0\right\}\right.}^e\stackrel{-}{t}dxdy$, $C_e=\frac{1}{2}{\underset{e}{\iint }\left\{\left.{\epsilon }_0\right\}\right.}^T\left[D\right]\left\{\left.{\epsilon }_0\right\}\right.\stackrel{-}{t}dxdy$, Eq. (7) can be processed as:

The work of external load is:

In equation, $\{R{\}}_{pe}$ is load vector of the unit volume force $\left\{p\right\}$. $\{R{\}}_{qe}$ is load vector of the surface force $\left\{q\right\}$ on element boundary. $\{R{\}}_{Me}$ is load vector of the concentrated force on node $M$.

Total potential energy $U$ of the element $e$ is equal to the sum of strain energy and external force energy:

The symbol ${\left\{\left.F\right\}\right.}^e$ is junction force, it is the superposition of various loads acting on the element $e$:

Then Eq. (10) can be simplified as:

Potential energy of the entire elastomer is obtained by adding the $U_e$ in all $E_o$ cells. It is expressed as:

According to the principle of minimum potential energy:

So, it can be obtained:

Compared with the elastic formula, load vector $\left\{\left.L\right\}\right.$ is only added to thermoelastic formula. From the significance of numerical calculation, coupling effect of temperature field is caused by the “equal” load array due to the change of temperature. So, the calculation of transient temperature and the calculation of internal force and deformation can be performed independently. It is necessary to calculate the history of transient temperature firstly. Then the temperature of time-discrete point is selected as the load value, stress and strain can be obtained by the corresponding relationship.

2.2. Equation of natural convection heat transfer

In the high-pressure rotor cavity, natural convection satisfies the conservation of mass, energy and momentum. So, they can be described by mathematical expressions.

The law of mass conservation denotes that mass increase in the fluid micro-element per unit time is equal to net mass of the micro-element at the same time interval. The continuum of mass conservation can be described in tensor form:

For natural convection with small temperature differences and low speeds, the assumption of incompressibility is appropriate. So, the continuous Eq. (16) can be simplified as follows:

State equation is:

The law of momentum conservation denotes that the change rate of fluid momentum in a micro-body is equal to the sum of various forces acting on the micro-body. The momentum equation can be expressed based on Newton’s second law:

$D{u}_i/Dt$ is derivative operator with the object, it is simplified as:

In Eq. (20), volume force $F_i$ is only the buoyancy in $y$ direction.

The law of energy conservation denotes that the rate of energy increasing in the micro-element is equal to the sum of net heat entering the micro-element and the work done by physical forces on the micro-element. This law is actually the first law of thermodynamics. Energy equation can be described by the variable of temperature $T$:

Because the velocity of natural convection is generally small in a cavity, viscous energy dissipation and pressure work can be negligible based on incompressible assumptions. Energy Eq. (21) is simplified as follows:

In order to analyze the various control equations, general form of each basic equation is established firstly. Assuming that $\phi$ is generic variable. The continuous, momentum and energy equations can be described as a unified transport equation:

where $\partial \phi /\partial t$ is the transient item, $\mathrm{d}\mathrm{i}\mathrm{v}\left(\overrightarrow{v}\phi \right)$ is the convective term, $\mathrm{d}\mathrm{i}\mathrm{v}\left(\mathrm{\Gamma }\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}\phi \right)$ is the diffusion term, and $S$ is the source item. The definition of transporting equation is shown in Table 1.

3.1. The establishment of physical model

Fig. 1 shows the 3D model of high pressure rotor. Symbol 1 indicates the section of front journal, symbol 2 indicates the first stage disk, symbol 3 indicates the second stage disk, symbol 4 indicates the third stage disk, symbol 5 indicates the fourth stage disk, symbol 6 indicates the fifth stage disk, symbol 7 indicates the sixth stage disk, symbol 8 indicates the seventh stage disk, symbol 9 indicates the eighth stage disk, symbol 10 indicates the ninth stage disk, symbol 11 indicates the high pressure turbine. It can be found that the model is symmetric about the central axis, so the structure is simplified by a 2D model.

In this research, thermo-elastic coupling element Plane77 is applied to mesh the high pressure rotor. For improving the calculation accuracy, here the number of elements is 2203, and the number of nodes is 7602. Fig. 2 shows physical model of the high pressure rotor in natural convection heat transfer.

Fig. 1Cutaway view of 3D solid model

Fig. 2Finite element model of high pressure rotor

As the model is very complicated, the calculation requires plenty of time and is also difficult to converge. So, applying the symmetry of geometric model, the 2D cavity model is adopted to research natural convection heat transfer in the high-pressure rotor, as is shown in Fig. 3(a).

The cavity of high pressure rotor is meshed by flow and thermal coupling element Fluid141, and the CFD model is obtained as shown in Fig. 3(b). In the meshing process, solid wall and the region where the flow changes significantly are meshed finely. Here the number of elements is 11310, and the number of nodes is 12442.

Fig. 3Physical model of nature convection for high pressure rotor

a) The two-dimensional geometric model

b) The cavity CFD model of high pressure rotor

3.2. The definition of boundary conditions

The heat exchanges between turbine disk and boundary airflow mainly include heat exchange between turbine disk surface and cooling air, heat exchange between outer edge surface and blade root, heat exchange between the outside of rotor drum and cooling air, heat transfer between the passage and cooling air. The formulas of heat transfer applied in each part are as follows:

(1) Heat exchange between turbine disk surface and cooling air.

According to the principle of convective heat transfer [31], total heat on the surface of turbine disk is expressed as:

Total heat on the surface is transmitted by conduction, so it is also expressed as:

Further processing, it can be obtained:

For making the equation dimensionless, the radius of local wheel $R_b$ is multiplied on both sides:

The right side of Eq. (27) is the ratio of surface temperature gradient to the reference temperature gradient. The left side is the ratio of conductive thermal resistance to fluid convection thermal resistance, and it is defined as Nusselt number $Nu$. So, it is obtained:

By the empirical formula [32], the value of Nusselt number $Nu$ can be determined by Reynolds number $R_e$. When $R_e$ < 2.4×10⁵, $Nu=$ 0.675×$R_e$^0.5, and when $R_e$ > 2.4×10⁵, $Nu=$ 0.0217×$R_e$^0.8:

where ${\lambda }_c$ is thermal conductivity of cooling air, $R_b$ is the radius of local wheel, $u$ is circumferential speed at wheel radius $R_b$, $v$ is kinematic viscosity of cooling air.

(2) Heat exchange between outer edge surface and blade root.

The heat transfer coefficient between outer edge and blade root is converted into a third type of boundary condition between outer edge and the gas. Through the convective heat transfer between the gas and blade, the heat is transferred from the blade to blade root, and then to the outer edge of turbine disk. So, the direct equivalent relationship of heat transfer is established between the gas and outer edge of turbine disk. The equivalent third type of boundary condition is:

where $h_e$ is the equivalent coefficient of convective heat transfer and is a function of the gas and blade parameters, $T_g$ is the temperature of gas, other parameters have the same meaning as before.

Fig. 4 shows a typical operating condition curve for aeroengine high-pressure rotor. The high pressure rotor accelerates from 0 rpm to the idle speed 10500 rpm. After a short constant speed platform, it continues to accelerate to 15000 rpm, and runs for 690 s at a constant speed of 15000 rpm. Then the rotor slows down from 15000 rpm to 0 rpm, the whole process has gone through 915 s.

Boundary conditions are changed with the engine state, while heat transfer empirical formula, calculation parameters and methods are the same in thermal boundary conditions. So, the heat transfer boundary for typical working condition can be obtained by the steady state condition, and other heat transfer boundary condition can be obtained by linear interpolation. In this research, the calculation of temperature is a pseudo-transient temperature field.

Temperature distribution of rotor outer wall at the last stop moment is taken as the initial boundary condition, as is shown in Fig. 5. The temperature of casing inner wall is constant at 20 ℃, the other solid walls are insulated. The velocity of all solid walls is assumed to be no slip, and that is 0 m/s. In the fluid region, physical properties are adopted the ideal gas parameters at 20 ℃.

4.1. Thermal-elastic coupling characteristics at typical start and stop conditions

(1) The distribution of transient temperature at typical start and stop conditions.

Temperature curves at typical operating conditions are obtained on different nodes, as are shown in Fig. 6.

Fig. 6(a) shows the temperature curves of different nodes in the whole process. During the process of start, the highest temperature is at outer edge of high pressure turbine. This zone exchanges heat with high temperature gas, so the temperature is higher, and is basically synchronized with the boundary temperature.

Fig. 4Speed curve during star and stop times

Fig. 5Boundary condition of initial temperature

When start process is completed, the highest temperature of zone reaches a stable value. The time required for temperature field to stabilize is longer than the stabilization time of heat exchange boundary condition for 150 s. The stabilization time of highest temperature is exceeded by 190 s. It also shows that the temperature rise rate and warm-up are significantly reduced in the start process. Due to the vibration of thermal bending, it usually makes the amplitude increase significantly in the process of slow speed start. So, the temperature change is given in this start period, as is shown in Fig. 6(b).

Fig. 6(c) shows the temperature curves of different nodes at the stop process. As the stopping time is short, the whole rotor has a higher temperature when the stop speed is reduced to 0 rpm. For the temperature difference of high pressure rotor, natural convective heat transfer occurs at the process of natural cooling.

Fig. 6Temperature curves of high pressure rotor at start and stop process

a) At the whole process

b) At the slow speed start process

c) At the stop process

Fig. 7 shows the temperature distribution of high pressure rotor at different moments. It can be found that isoline density of the temperature at the third stage disk is larger in rotor start process, and indicates that the change of temperature is more obvious in this region. Rotor temperature increases gradually from the lower stage to the higher stage, and the temperature of outer edge in each disk is higher than the temperature of inner edge.

When start process reaches a steady state, the lowest temperature appears in the region of front section, temperature value is 357.6℃. The highest temperature appears in the outer edge at the first stage disk, the temperature value is 601.6℃, and it is shown in Fig. 7(c). Because the time of stop process is shot, temperature change isn’t shown a constant speed (10500 rpm) platform for 50 s, and cooling rate is always large. Fig. 7(f) shows the final temperature distribution at the end of whole process. It is adopted as the temperature boundary condition for flow-thermal coupling calculation in natural convective heat transfer.

Fig. 7Temperature distributions of high pressure rotor at different start moments

a)$t=$ 30 s

b)$t=$ 60 s

c)$t=$ 130 s

d)$t=$ 400 s

e)$t=$ 850 s

f)$t=$ 915 s

(2) Thermal-elastic coupling characteristics at start and stop process.

The calculation of stress field is taken a sequential coupling approach. Based on the distribution of temperature field, node temperature is loaded to the model of stress calculation. Then the degree of freedom is constrained, blade centrifugal force and fluid pressure of external flow field are ignored.

Transient displacements of high pressure rotor at different nodes and moments are shown in Fig. 8. It can be found that radial displacement of each node is basically in the same law, the change is significant at start and stop process. The radial displacement curves of high pressure rotor are shown in Fig. 9.

In Fig. 9(a), radial displacement of each node gradually increases with the increase of start time. When rotational speed is 10500 rpm, the value of radial displacement is decreased suddenly. It is because the coefficient of convective heat transfer at this moment is kept constant. The change of each node temperature is decreased, and the difference of temperature is gradually reduced, so the displacement is gradually decreased. Then the temperature is gradually increased with rotational speed, the difference of temperature is increased, and radial displacement is increased.

Fig. 9(b) shows the radial displacement curves of high pressure rotor at stop process. As the time of stop process is shot, radial displacement isn’t shown a constant speed 10500 rpm platform for 50 s, and is increased with stop time. This illustrates that temperature difference, thermal stress and thermal deformation are much large in stop process.

Fig. 8Transient displacement of high pressure rotor at different nodes and moments

Fig. 9Radial displacement distributions of high pressure rotor at start and stop process

a) The process of start

b) The process of stop

Fig. 10 shows response curves of the first stage disk on high pressure rotor. It can be found in Fig. 10(a) that the radial and axial stresses have the same law, and the larger thermal stresses are generated at the start and stop process. The maximum stress occurs at the moment when the speed is accelerated to 15000 rpm (starting time 150 s). The value of rotor stress becomes small with the stabilization of rotor temperature. Fig. 10(b) exhibits that radial displacement of the rotor is much larger than axial displacement. Total displacement is basically the same as radial displacement, and axial displacement has little effect on the total displacement. Hence, the effect of axial thermal deformation on high pressure rotor is ignored in the research of thermal vibration.

Fig. 10Response curves of the first stage disk on high pressure rotor

a) The curves of thermal stress

b) The curves of displacement

Fig. 11 shows displacement distributions of high pressure rotor at different moments. It reflects the change of displacement at different moments.

Fig. 11Displacement distributions of high pressure rotor at different moments

a)$t=$ 30 s

b)$t=$ 160 s

c)$t=$ 600 s

d)$t=$ 850 s

e)$t=$ 910 s

4.2. Flow-thermal coupling characteristics of natural convection conditions after stopping

(1) Flow law of rotor cavity at stop process.

Fig. 12 shows the distribution of steady velocity contour in high pressure rotor cavity. Due to the complexity of calculating geometric model and temperature boundary condition, velocity distribution of the whole computational domain is also very complex. In particular, velocity contour of each disk cavity are distributed densely, and it changes more obviously.

Fig. 13 shows the distribution of velocity vector in rotor cavity. In Fig. 13(a), it is macroscopically exhibited the velocity vector distribution of natural convective air. As physical model of the high pressure rotor is much complex, the density of discretized grid is large. In order to more clearly understand the velocity distribution of air flow, velocity vector distributions of different typical regions are shown in Fig. 13(b).

Fig. 12Velocity contour distribution in high pressure rotor cavity

Fig. 13Velocity vector distribution in high pressure rotor cavity

a) The diagram of macroscopic velocity vector

b) Sketch map of different typical regions

Fig. 14 shows the local magnification diagrams of velocity vector at different regions. In Fig. 14(c), a large typical ellipse vortex is formed between upper and lower surfaces. Its flow characteristic is similar to natural convection in a square cavity. But the position of large ellipse vortex is taken a certain offset, and it is located at bottom right position of the cavity.

Fig. 14(a) and (d) exhibit the flow characteristics of disk cavities in the upper and lower regions. In the adjacent upper cavities formed by the 7th, 8th and 9th disks, the characteristics of air flow are basically similar, and a typical ellipse vortex is formed in the cavity. It is different in the two adjacent lower cavities, three ellipse vortexes are formed in the high temperature disk cavities. While only one ellipse vortex is formed in the low temperature region, and the flow velocity also has a significant difference.

Fig. 14(b) shows the velocity vector distributions in the cavities of 3rd disk and 4th disk. When the air flows from narrow passage, most of it flows into the cavities of 3rd disk and 4th disk, and immediately forms two vortexes in left and right sides. Only a little amount of air flows into the small gap formed by high pressure rotor and the casing.

In Fig. 14(e), a large ellipse vortex is formed in the E region, and the velocity of flow is large. Flow characteristics of other areas are not described in detail. In short, due to geometrical complexity of high pressure rotor and nonlinear variation of temperature boundary along the axial direction, natural convection of the air is very complicated in high pressure rotor cavity. Many annular vortexes are formed in the whole cavity, which is completely different from natural convection of simple square cavities. It shows the diversity and complexity of flow.

Fig. 14Local magnification diagrams of velocity vector at different regions

a) Region A

b) Region B

c) Region C

d) Region D

e) Region E

f) Region F

g) Region G

h) Region H

(2) Flow-thermal coupling characteristics at natural convection condition.

Different coefficients are selected as convective boundary conditions in each rotor part, which is obtained by the calculation of natural convection flow-thermal coupling. The distribution of thermal-elastic coupling field is analyzed at different stop moments. Fig. 15 shows the temperature distribution in rotor upper region.

At the beginning of stop, the temperature of rotor is much higher. The density of temperature contour is larger, especially at the 3rd and 4th disks. The temperature of each node gradually decreases with the increasing of stop time. The density of temperature contour decreases in the same time, it illustrates that temperature distribution of high pressure rotor gradually becomes uniform. By comparison, it is found that the reduction of temperature is relatively smaller at the 3rd and 4th disks, and the density of contour is much larger. This is mainly because the velocity of air flow is relatively small at the 3rd and 4th disks, and the coefficient of convective heat transfer is small. Hence, it makes drop amplitude of the temperature in a lower value.

Fig. 15Temperature distributions of high pressure rotor at different stopping moments

a)$t=$ 0 s

b)$t=$ 100 s

c)$t=$ 300 s

d)$t=$ 1000 s

e)$t=$ 3000 s

Fig. 16 shows temperature curves of different locations after stopping. The temperature of each node gradually decreases with the increase of stopping time. Cooling speed of the temperature is the slowest at the 3rd disk, while cooling speed is the fastest on the 8th disk. After stopping time reaches 90 minutes, temperature change of each node becomes much slower. Eventually the temperature of each node is equal to the temperature of room, and the heat exchange tends to be balanced.

Fig. 17 shows the curves of temperature difference between the upper and lower surfaces after stopping. The difference of temperature isn’t gradually decreased with the increase of stopping time, it presents a phenomenon of increasing firstly and then decreasing. Besides that, the maximum temperature difference of different positions appears at different moments. Here the maximum temperature difference is about 23 ℃ at the 3rd disk, and appears near 60 minutes after stopping. The minimum temperature difference is about 4.8 ℃ at the 8th disk, and appears near 15 minutes after stopping. It also can be found that temperature difference at the 3rd disk is larger in the stopping time of 40-90 min, and is more than 20 ℃. This is basically consistent with the operating experience, so the correctness of calculation is verified.

Fig. 16Temperature curves at different moments

Fig. 17Difference curves of the temperature

Fig. 18 shows thermal stress curves of high pressure turbine in $x$ and $y$ directions. After stopping, the cold air and hot air form a convection process in the high pressure turbine cavity. It makes the temperature difference of high pressure turbine change alternately, so thermal stress fluctuates greatly at the beginning of stopping period. The convection of cold air and hot air becomes slower with the time consumption, and temperature difference becomes small and stable. When stopping time reaches 1500 s, stress value remains basically unchanged. It illustrates that heat exchange has been completed at this moment, and temperature difference of the rotor changes little.

Fig. 18Thermal stress curves of high pressure turbine in x and y directions

a) In the $x$ direction

b) In the $y$ direction

Fig. 19 shows temperature curves of idle speed in restart process after stopping. The temperature of each node decreases firstly and then gradually increases in the restart process after parking 300 s, as is shown in Fig. 19(a). The time corresponding to the lowest temperature is also different in each position. Here it is the longest in the 3rd disk, and is about 35 s. While it is the shortest in the 8th disk, and is about 12 s. Except that the temperatures of the 3rd disk and 4th disk decrease firstly and then increase gradually when the parking time is increased to 1200 s in Fig. 19(b), the other nodes are gradually increased with starting time. When the stopping time is increased to 5000 s in Fig. 19(c), temperature change of each node has become balanced on high pressure turbine, and it is equal to room temperature. So, the start at this moment is equivalent to a cold start.

Fig. 20 shows difference curves of temperature in restart process at different stop moments and positions. Fig. 20(a) exhibits temperature difference of the 3rd disk gradually decreases with the time of restart after different stop moments. When start time increases to 50 s, temperature difference between the upper and lower surfaces is 0 ℃ at the 3rd disk. Fig. 20(b) exhibits temperature difference of different position in restart process after stopping 3600 s. As the start time increases, temperature difference between the upper and lower surfaces gradually decreases. Here temperature difference of the 3rd disk is the largest, and temperature difference of the 6th disk is the smallest. When the starting time is increased to 50 s, temperature difference is the smallest in each node.

Fig. 19Temperature curves of idle speed in restart process

a) After stopping 300 s

b) After stopping 1200 s

c) After stopping 5000 s

Fig. 20Difference curves of temperature in restart process

a) At the 3rd disk

b) At different positions