Abstract
The rod-fastened rotor is the core component of the gas turbine. It comprises several discs and tie rods. The flexural stiffness of the contact interface is the key factor for rotordynamic analysis. The contact interfaces of the discs are usually manufactured by grinding. The measured contour curve of the contact interfaces of an experimental rod-fastened rotor is analyzed by the structural function method, which shows that the contact interfaces can be well described by the double fractal model with fractal dimensions D1 and D2 and the fractal roughness parameters G1 and G2. The Hertz model is used to analyze the contact of the single asperity on the contact interface. On this basis, the flexural stiffness of the contact interface considering the pretightening force and the bending moment is derived. Modal frequencies of the experimental rod-fastened rotor under different pretightening forces and the bending moment (caused by gravity) are obtained by three-dimensional finite element analysis and experimental modal tests. It is observed that the modal frequencies increase with the nominal pressure of the contact interface, and the experimental results are consistent with the calculated results.
1. Introduction
The rod-fastened rotors are widely used in aeroengines and heavy-duty gas turbines (see Figure 1). In order to avoid the damage caused by flexural vibration, it is necessary to accurately calculate the flexural vibration modal frequencies of the gas turbine rotor at the design stage. Because the rotor is not a continuous whole, the contact effect between two discs has a great influence on the dynamic characteristics. There are some researches about the dynamic performance of such structures. Klompas [1] analyzed the influence of the bending moment on the dynamic characteristic of the rotor joints in the turbomachine. Lee and Lee [2] investigated the effect of the central tie rod on the rotordynamic performances of an auxiliary power unit gas turbine rotor. Lu et al. [3] calculated the modal frequencies of a rod-fastened rotor subjected to different pretightening forces by the three-dimensional finite element method. Liu et al. [4] studied the dynamic stability of a rod-fastened rotor bearing system considering the nonuniform pretightening forces. However, the effect of the contact stiffness on the dynamic characteristic of the rod-fastened rotor was not considered.
For the contact stiffness of the contact interface, some researchers have studied the contact stiffness in microscale. Rao et al. [5] calculated the contact stiffness using the Hertz elastic contact theory and Greenwood–Williamson statistical model. He et al. [6] studied the contact stiffness of rough contact interfaces using Greenwood–Williamson model and the elastic-plastic contact model. Zhang et al. [7] identified the contact stiffness by modal test and obtained the relationship between the contact stiffness and the pretightening force. There are also some researches to analyze the flexural stiffness considering the bending moment. Isa et al. [8] established a bilinear stiffness model to calculate the flexural stiffness of a rod-fastened rotor subjected to the bending moment. Gao et al. [9] studied the effect of the bending moment and the pretightening force on the flexural stiffness of the contact interface in the rod-fastened rotor. Several researchers calculated the equivalent stiffness of bolted structures by the finite element method. Lehnhoff et al. [10, 11] calculated the equivalent stiffness of the bolted joints using the axisymmetric finite element model. Yuan et al. [12] calculated the contact stiffness and flexural stiffness of the gas turbine rotor with curvic couplings using the three-dimensional finite element model.
In the above studies, the contact stiffness is calculated using Greenwood–Williamson statistical model. However, the contact interfaces of the discs in the actual gas turbine are manufactured by grinding. Due to the different scales of plastic deformation and microfracture during processing, such machined surfaces generally have double fractal characteristic [13]. Jiang et al. [14] calculated the normal contact stiffness of the joint surface based on the fractal model and compared the theoretical calculation results with the experimental data. Buczkowski et al. [15] used the Weierstrass–Mandelbrot function to express the fractal features of the rough surface and calculated the normal contact stiffness of the isotropic rough surface during extrusion deformation.
However, there are few investigations that focus on the analysis of the flexural stiffness of the rod-fastened rotor, considering the pretightening force and the bending moment based on double fractal model. In this paper, an experimental rod-fastened rotor is designed. The contact interfaces of the rotor are grinding. The contour curve of the contact interfaces is analyzed by the structural function method, which shows that the contact interfaces have double fractal characteristic. The normal contact stiffness of the contact interface is obtained by the Hertz contact theory and the double fractal model. The flexural stiffness of the contact interface considering the pretightening force and the bending moment is also derived. Finally, the modal frequencies of the rotor under different pretightening forces are obtained by the three-dimensional finite element method and the experimental modal method, respectively. Comparing the experimental test results with the calculated results, it is found that the double fractal model is reasonable to calculate the contact stiffness of the rod-fastened rotor.
2. Theoretical Analysis
2.1. The Experimental Rod-Fastened Rotor
The experimental rod-fastened rotor consists of discs, a front shaft end, a rear shaft end, and nine tie rods (see Figure 2). The rotor is clamped together by the circumferentially arranged tie rods. A rough contact interface is existed between two discs of the rotor. As shown in Figure 3, there are six rough contact interfaces of C1–C6 in the rotor. The structure of the rotor is bilaterally symmetrical. The material of the rotor is 40Cr.
(a)
(b)
The elastic modulus E is 206 GPa, the density ρ is 7870 kg/m3, and the Poisson ratio is 0.3.
2.2. The Mechanical Model of the Contact Segment in the Experimental Rod-Fastened Rotor
There are some contact asperities on the rough contact interface. The asperities come into contact under the pretightening force Fpre. According to the fractal model, a small plastic contact point first appears on the contact interface. As the pretightening force increases, the small plastic contact points merge into large contact points. When the contact area exceeds the critical value, the large contact points begin to make elastic contact. As shown in Figure 4, each asperity on the rough contact interface can be seen as a spring with bending stiffness, and all asperities on the contact interface can be characterized by a spring with equivalent flexural stiffness Kcc. Then, the equivalent flexural stiffness Keq of the contact segment is a series connection of the flexural stiffness Kcc of the contact interface and the flexural stiffness Ks of the continuous shaft with the length Ls. Keq is given by
For the continuous shaft with the length Ls, the flexural stiffness iswhere E is the elastic modulus, I is the moment of inertia, and Ls is the length of the contact segment.
Considering the effect Kcc on the Keq, the flexural stiffness correction factor of the contact segment η is introduced by
The flexural stiffness correction coefficient η represents the stiffness reduction of the contact interface to the continuous shaft. For the contact segment, η is mainly related to the pretightening force and the bending moment.
2.3. The Equivalent Flexural Stiffness of the Contact Interface Kcc
2.3.1. The Double Fractal Model of the Contact Interface in the Experimental Rod-Fastened Rotor
A rough surface profile with double fractal characteristic can be represented by Weierstrass–Mandelbrot (W-M) function [16]:where L is the sampling length, G1 and G2 are the fractal roughness parameters which determine the profile height, D1 and D2 are the fractal dimensions of the surface topography (1 < D1 < 2, 1 < D2 < 2), which determine the distribution ratio of high frequency and low frequency in the surface contour height, and γ is the contour spectral density.
The contact interfaces of the experimental rod-fastened rotor are machined by grinding, and the rough contact interfaces can be regarded as isotropic. The contour curves of the contact interfaces are measured by the SJ201P surface topographer. According to the contour curve, the fractal parameters can be obtained by the structural function method.
When τ ≤ τ12, the structural function of the W-M function is
When τ ≥ τ12, the structural function of the W-M function iswhere τ is the displacement along the x direction, C1 and C2 are the coefficients related to the fractal dimension, and < > indicates spatial average. The structural function diagram of the contour curve is shown in Figure 5. The structural function diagram is divided into two parts according to τ12 which represents the boundary displacement of the two fractal regions along the x direction. The region I corresponds to the fractal structure formed by the microfracture of the smaller scale, and the region II corresponds to the fractal structure formed by the plastic deformation of the larger scale. The contact of two rough surfaces can be equivalent to the contact between a rough face and a rigid face, and the structural function is the sum of the two rough face structure functions. According to the structural function diagram, the fractal parameters of the contact interface are obtained and shown in Table 1.
2.3.2. The Contact Model of the Double Fractal Surface
The contact of two rough surfaces can be equivalent to the contact between a rough face and a rigid face. Figure 6 shows the microcontact established between an asperity of the rough surface and opposing rigid plane. The contact asperity is assumed to be a sphere, when the small contact asperities are in the state of full plastic.
Deformation of the normal contact force is given bywhere H is the hardness value of the material, , Y is the yield stress, E is the elastic modulus, ν is the Poisson ratio, and is the truncated area of the contact asperity, is the radius of the truncated area.
When the pretightening force is increased, the small contact asperities are merged into a large contact asperity, and the large contact asperity begins to make elastic contact. When aL ≤ ac, the contact interface is in a fully plastic contact state. When aL > ac, the contact interface enters the elastic-plastic contact state, where aL is the truncated area of the maximum contact asperity and ac is the critical truncated area where the elastic contact occurs.
For the contact asperity, the normal deformation δ can be expressed as [17]
According to the Hertz contact theory, the normal elastic force of a single contact asperity can be expressed aswhere is the equivalent elastic modulus.where , and E1, E2 are Poisson’s ratio and elastic modulus of the contact interfaces, respectively. Because the same material is used for each part of the experimental rod-fastened rotor, E1 = E2 = E and = = .
In Figure 6, the relationship between and δ iswhere the radius of curvature Re is general much larger than δ. Then, equation (11) can be expressed as
Substituting equations (8) and (12) into (9) yields
The critical truncated area ac is given by [17]
For the double fractal contact interface, region I is corresponding to the critical truncated area ac1, region II is corresponding to the critical truncated area ac2, and the boundary displacement τ12 corresponds to the truncated area a12. The asperities deformation of the double fractal area can be divided into four cases according to the relationship among ac1, ac2, and a12: Case 1: when ac1 < ac2, the contact asperities satisfying a < ac1 in region I are plastic deformation, and the contact asperities satisfying ac1 < a < a12 are elastically deformed. Case 2: when ac1 > a12, all contact asperities in region I are plastically deformed. Case 3: when ac2 < a12, all contact asperities in region II are elastically deformed. Case 4: when ac2 < a12, the contact asperities satisfying a12 < a < ac2 in region II are plastically deformed, and the contact asperities satisfying a > ac2 are in the elastic deformation state.
According to the fractal parameter values of the double fractal contact interface listed in Table 1, the value of truncated areas ac1, ac2, and a12 can be obtained, respectively, and the relationship is a12 < ac2 < ac1. According to the relationship, the double fractal area of the contact interface satisfies the case 2 and case 4.
For the truncated area aL of the maximum contact asperity, when aL < a12, the frequency distribution function of the contact asperities is given by [18]
When aL > a12, the frequency distribution function of the contact asperities is expressed as
The total normal contact force of the contact interface FN is the sum of the plastic contact force and the elastic contact force of the two regions. FN can be expressed as
FN is determined by the pretightening force Fpre. Combining with equation (17), aL can be obtained. It will be applied to the following contact stiffness calculation.
2.3.3. The Flexural Stiffness of the Contact Interface of the Experimental Rod-Fastened Rotor
From the Hertz contact theory, the normal contact stiffness of each microcontact pair can be written as
Neglecting the contact stiffness of the asperities in the state of plastic deformation, the normal contact stiffness Kn of the entire contact interface is equal to the sum of the normal contact stiffness of all elastic contact pairs of the contact interface. The contact stiffness Kn can be expressed as
The value of Kn can be obtained by the method of numerical integration. For the annulus contact interface, the flexural stiffness of contact interface Kcc can be given bywhere A is the nominal contact area of the contact interface and y is the distance along the axis.
The value of Kcc is determined by the normal contact stiffness per unit area kn0 and the contact area. Since kn0 is a function of the contact pressure P, Kcc also can be expressed as
As shown in Figure 7, the contact interface of the rod-fastened rotor is a concentric annular surface and the outer and inner radius of the annulus are R1 and R2, respectively. When the contact interface is in the full contact state, the contact pressure is a linear function of y, thenwhere a and b are the constant coefficients.
The rod-fastened rotor is subjected to the pretightening force and the bending moment during the operation. When the pretightening force Fpre and the bending moment M are applied on the rod-fastened rotor, Fpre and M can be written aswhere Sx = Ayc is the static moment of the contact interface and yc is the coordinate of the centroid. Combining equations (23) and (24), the coefficients a and b can be obtained. Substituting them into equation (22) yieldswhere .
When the whole zone of the contact interface is in full contact, yc = 0. Equation (25) can be written aswhere ξ is the flexural dimensionless load coefficient, , and .
For the annulus contact interface,
Combining equations (21) and (27), the Kcc is expressed as
3. Application of Results and Validation
3.1. Equivalent Flexural Stiffness Keq of the Contact Segment
The flexural stiffness Kcc of the contact interface in the experimental rod-fastened rotor can be obtained by the above calculation method. Combining with equation (1), the equivalent flexural stiffness of the contact segment Keq can be obtained. Keq will be affected by the pretightening force, the bending moment, the contact area, and the length of the contact segment. To analyze the effect of above factors on Keq, the contact segment corresponding to the contact interface C2 is selected to conduct the analysis.
3.1.1. Effect of the Pretightening Force and the Bending Moment on Keq
The pretightening force and the bending moment applied on the experimental rod-fastened rotor will affect the contact pressure of the contact interface. In equation (26), the flexural dimensionless load coefficient ξ is introduced to represent the influence of the pretightening force and the bending moment. The effect of the contact interface on Keq is denoted by the flexural stiffness correction coefficient η (see equation (3)). For the contact interface C2, the value of R1 and R2 is 110 mm and 70 mm, respectively. The nominal contact pressure P0 can be obtained under different pretightening forces. Figure 8 shows the effect of the pretightening force and the bending moment on η. When ξ is constant, η increases with P0, but the increasing trend gradually slows down. When P0 is constant, η decreases with the increase of ξ, which indicates that the contact interface will appear to be separated with the increase of the bending moment.
3.1.2. Effect of the Contact Area on Keq
According to equation (28), the flexural stiffness Kcc is also influenced by the contact area. The outer radius R1 is given as 110 mm, and the ratio of the inner radius and outer radius α = R2/R1 is defined. The change of the contact area is obtained by changing α. Figure 9 shows the effect of the contact area on η. When ξ is constant, η increases with the decrease of α, indicating that the bending moment weakens against the bending stiffness with the increase of the contact area.
3.1.3. Effect of the Length of Contact Segment Ls on Keq
According to equations (1) and (3), η is influenced by the length Ls of the contact segment. As shown in Figure 10, when Ls is small, the length of the contact segment has a greater influence on η. When Ls reaches a certain value, η no longer increases with Ls. Therefore, it is necessary to select the appropriate length of the contact segment during the calculation. In this paper, Ls is equal to 50 mm.
3.2. Modal Analysis of the Experimental Rod-Fastened Rotor
The modal analysis of the rotor is conducted, considering the effect of the contact stiffness. In this paper, the effect of the contact interface on the flexural stiffness is simulated by modifying the elastic modulus of the contact segment. According to equation (3), the flexural stiffness correction coefficient η is obtained. Then,where Eeq is the equivalent elastic modulus of the contact segment.
The modal analysis is performed, considering the effect of the pretightening force and the gravity bending moment. There are six contact segments in the experimental rod-fastened rotor (see Figure 3). Due to the bilateral symmetrical structure, we only change the pretightening forces corresponding to tie rod set 2 and tie rod set 3, and the pretightening force of tie rod set 2 is equal to that of tie rod set 3. The pretightening forces corresponding to tie rod set 1 and tie rod set 4 are kept equal and constant. The different pretightening forces applied on the rotor are listed in Table 2.
The distribution of the gravity bending moment of the experimental rod-fastened rotor is shown in Figure 11. The distribution of the gravity bending moment is bilaterally symmetric. The gravity bending moments corresponding to the positions of C2 and C3 are selected as the bending moments of the contact interfaces. The length Ls of the two contact segments is both 50 mm. According to the pretightening force listed in Table 2, the flexural stiffness correction coefficients η of the two contact segments are listed.
As shown in Table 3, η increases with the pretightening forces.
The modal analysis of the experimental rod-fastened rotor is conducted using ANSYS Version 11.0. The block Lanczos method is applied to solve the modal analysis. The three-dimensional finite element model is shown in Figure 12. The first five flexural mode shapes of the experimental rod-fastened rotor are shown in Figure 13. It can be seen that the first five mode shapes are mainly the vibration of the front and rear shaft ends, and the mode shapes of the discs are not obvious. Table 4 shows the calculated values of the first five orders of modal frequencies of the experimental rod-fastened rotor under different pretightening forces. As the pretightening force increases, the modal frequency of the rotor gradually increases. Meanwhile, the influence of the pretightening force on the 3rd and 5th order modal frequencies is greater than that of other orders. This can be explained from the mode shapes of the rotor in Figure 13. Since the 3rd and 5th order mode shapes of the discs are larger than the other orders, the frequencies of the 3rd and 5th order are more sensitive to the change of the contact stiffness.
3.3. Experiment Validation
To validate the calculation method considering the pretightening force and the bending moment, a simplified experimental rod-fastened rotor is manufactured (see Figure 2(a)), and the modal experiment is conducted. The multipoint exciting method is used to measure the modal parameters of the rotor. The rotor is excited by the force hammer. An acceleration sensor is used to acquire the response signals. The range of the measured frequency is 0–3200 Hz with a frequency resolution of 1.25 Hz. The measurement system is shown in Figure 14.
In order to verify the calculated results, the pretightening forces applied on the experimental rotor are set according to Table 2. In order to reflect the effect of gravity bending moment, the experimental rotor is suspended by an elastic wire rope. The experimental results of the first five orders of modal frequencies are shown in Table 5. The comparison between the calculated and experimental results is shown in Figure 15. It can be seen that the modal frequencies gradually increase with the nominal contact pressure P0, and the change trend of the calculated results is consistent with that of the experimental results. For the 3rd order and the 5th order, the difference between the calculated and experimental results is large. Because the 3rd and 5th order mode shapes of the discs are larger than the other orders, the frequencies of the 3rd and 5th order are more sensitive to the change of the contact stiffness, which is determined by the pretightening force.
(a)
(b)
(c)
The relative errors between the calculated and experimental results are listed in Table 6. For the 3rd and 5th order, the relative errors are large when the nominal contact pressure P0 (pretightening force) is small. It could be explained by the relative errors of pretightening forces applied on the rotor being large when the pretightening forces are small. With the increase of P0, the relative errors decreased. According to the results of relative errors, it is found that the calculated results are correct.
4. Conclusions
In this paper, the contact stiffness of the rod-fastened rotor is analyzed considering the pretightening force and the bending moment. The results are summarized below:(1)The double fractal model and Hertz contact theory are used to calculate the normal contact stiffness of the contact interface. The effect of the pretightening force, the bending moment, the contact area, and the length of the contact segment on the equivalent flexural stiffness Keq are analyzed. When the flexural dimensionless load coefficient ξ is constant, the flexural stiffness correction coefficient η increases with the nominal contact pressure P0 (pretightening force), but the increasing trend gradually slows down. When P0 is constant, η decreases with the increase of ξ, which indicates that the contact interface will appear to be separated with the increase of the bending moment.(2)The calculated and experimental modal analysis is conducted considering the pretightening force and the bending moment (due to gravity). As P0 increases, the modal frequencies increase, and the experimental results are consistent with the calculated results. Therefore, the double fractal model is reasonable to calculate the contact stiffness of the rod-fastened rotor.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Nature Science Foundation of China (no. 11872289) and the Key Project of Natural Science Foundation of Xi’an Jiaotong University (no. zrzd2017025).