Ynhong MA, Yongfeng WANG, Cun WANG, Jie HONG,*

a School of Energy and Power Engineering, Beihang University, Beijing 100083, China

b Collaborative Innovation Center of Advanced Aero-Engine, Beijing 100083, China

KEYWORDS

Abstract The dynamic influence of joints in aero-engine rotor systems is investigated in this paper.Firstly, the tangential stiffness and loss factor are obtained from an isolated lap joint setup with dynamic excitation experiments. Also, the influence of the normal contact pressure and the excitation level are examined, which revel the uncertainty in joints. Then, the updated Thin Layer Elements (TLEs) method with fitted parameters based on the experiments is established to simulate the dynamic properties of joints on the interface. The response of the rotor subjected to unbalance excitation is calculated, and the results illustrate the effectiveness of the proposed method. Meanwhile,using the Chebyshev inclusion function and a direct iteration algorithm,a nonlinear interval analysis method is established to consider the uncertainty of parameters in joints. The accuracy is proved by comparison with results obtained using the Monte-Carlo method. Combined with the updated TLEs,the nonlinear Chebyshev method is successfully applied on a finite model of a rotor.The study shows that substantial attention should be paid to the dynamical design for the joint in rotor systems,the dynamic properties of joints under complex loading and the corresponding interval analysis method need to be intensively studied.

1. Introduction

The high-level vibration of aero-engines caused by the vibration of a rotor system is a major concern of researchers.During the past decades, the rotor dynamic analysis is usually applied to a rotor system without joints and interfaces. However, the rotor shaft is designed with thin shells to reduce the weight of aero-engine, for which a stiffness decrease is unavoidable. The dynamic influence of a joint cannot be ignored, especially in the low-pressure rotors of high bypass turbine fan engines, which are a typically flexible rotor.

Many studies have been performed on the joint dynamics in astronautics,and a large number of phenomenological models have been established.1,2The rocket or pipe structures are usually modeled as beams; the joint is modeled by a spring and a sliding block, which can be implanted into the beam model easily.3Due to the relatively complicated construction of aero-engine rotors,three-dimensional finite element solid modeling is widely used to keep features of the structure. The implantation of a phenomenological model to corresponding nodes and the calculation will be very complicated and costly.

Therefore, an approach of modeling the joint dynamics called Thin Layer Elements(TLE)is proposed,which is firstly used in soil-structure interactions and rock joints.4,5The joint interface between two contacting bodies is modeled with continuous elements of very small but finite thickness.A special constitutive model is used and various deformation modes are incorporated. Ahmadian et al.6incorporated linear TLE with isotropic material properties into existing commercial finite element codes to establish a parameterized model of a pipeline with a bolted flange joint; the stiffness are identified by experiments,while the damping is ignored.Bograd et al.7–9uses frequency independent damping and linear joint stiffness parameters to predict the vibration response of structures,and the simulations are based on the experimentally determined joint parameters acquired from an isolated joint.A fast prediction of the eigenfrequencies and modal damping factors of the structures can be performed before there exists a physical prototype for experimental testing.Zhai et al.10establishes a finite element model of an aero-engine stator with a joint modeled by TLE. The complexity of the Finite Element (FE) model of an aero-engine stator with many joints is reduced. However, the linear assumption used in these analyses is inaccurate; the different deformation modes,such as stick,micro and macro slip,will lead to different localized shear stiffnesses and damping.Ahmadian and Iranzad11use elastoplastic materials that can describe the entire process from stick to macro slip to simulate the joint behavior at different states and achieves good agreement with the nonlinear response of experiments. Above all,the TLE can describe the dynamic of joints accurately and is compatible with solid elements for the prediction of the dynamic response of a system with high efficiency compared to the contact model.However,the damping of joints is usually ignored in an aero-engine, and the accuracy of nonlinear TLE and the efficiency of linear TLE could not be achieved simultaneously.

Additionally, the stiffness and damping of a joint or the parameters of equivalent TLE are the functions of the geometry, material, surface condition, load, etc. The sparse experimental data available on the very few joint geometries studied showed huge variability.1These data were collected from a small number of specimens, and testing on larger data sets might well have shown even more variability. Clearly,much more investigation by laboratory experiment must be performed in the future to obtain both the nominal properties and statistical statements of the variability of joints. Before that research is performed, the probabilistic methods12are invalid as the statistical distribution of parameters cannot be acquired. The interval analysis method, which only needs the bound of parameters, seems to be more effective.

The interval analysis method of linear systems is mature.Moens and Vande Pitte applied the interval method to the frequency response envelope analysis of uncertain structures.13Qiu et al.14,15developed the parameter perturbation method and Taylor’s expansion method, while the second order Taylor’s expansion16and updated second-order Taylor series expansion17are used to achieve higher accuracy. Ma et al.18established the interval analysis method on rotor dynamics,studied the modal characteristics, the steady response and the transient response using perturbation method, the interval modal superposition method and the Taylor expansion,respectively.To cope with the overestimation problem in Taylor series-based algorithms, Qiu and Qi19further proposed a collocation interval analysis method based on the Chebyshev polynomials. This method has excellent control in the overestimation compared with a Taylor series and has been used in the torsional vibrations of wind turbine geared transmission systems20and the start-up transient process of an overhung rotor.21Furthermore,it was used in the analysis of parametrically excited systems.22,23

The application of the interval analysis method on nonlinear systems is relatively less common. Qiu et al.24performed studies on the nonlinear structural dynamics with interval uncertain parameters using Taylor expansions. Wu et al.25–28introduced the Chebyshev inclusion function into nonlinear multibody mechanical problems. All the analyses are applied on a simple example with few degrees of freedom; the Taylor expansion method, especially, needs to be applied on the partial differential form of the original equation, which makes it approximately impossible to implant into commercial finite element software. However, the non-intrusive Chebyshev method has the potential to be used on complicated models.Until now, a nonlinear interval rotor dynamic analysis method, which is applicable to systems with huge degrees of freedom, has not been established; the uncertainty caused by the joint has not been considered in rotor either.

In this paper, an updated TLE method will be given based on the experimental data to consider both the accuracy and efficiency and to take the dynamic influence of joint into account in aero-engines. The nonlinear interval analysis method based on the Chebyshev expansion and commercial finite element software will be introduced to obtain the response on frequency domain of a rotor with uncertain joints.

2. Experimental investigation of dynamic characteristics on lap joints

2.1. Test setup and parameter estimation

The dumbbell joints experiments carried out by Gaul29and Sandia laboratory1are widely used to acquire the stiffness and damping of joints. The similar experimental layout was built, which consists of two small flat patches joined in the middle with a bolt, as shown in Fig.1. To excite the joint in the tangential direction, the parts of the joint are attached to additional masses one of which is excited with a shaker. The mass is excited by the shaker constructed with a leaf spring to amplify the response when the system vibrates in resonance.This design allows the creation of larger tangential forces in the lap joint, and provides a better signal to noise ratio for measurements in the axial direction. The parameters of the joint are sampled at the first order axial vibration frequency,which is 360 Hz according to the hammer test, with different contact pressures on the interface and excitation levels. This system behaves like a 1-DOF (degree of freedom) system with two masses connected by a spring with a friction element.

Fig.1 Test setup for determination of joint patch dynamic parameters.

The normal force in the joint is regulated using a bolt and is loaded using a torque wrench. As shown in Fig.2, three piezoelectric accelerometers 1–3# with the sensitivities of 2.828 m/s2,2.655 m/s2and 1.806 m/s2,respectively,are located on both sides of the joint. Integrating the acceleration signals twice with respect to time, the relative displacement in the joint can be calculated by the indefinite integration with same initial conditions,7shown as Eq.(1).During the measurement,the setup is supported with two thin wires at the centers of gravity of the masses to reduce the influence of the support stiffness.

In fact,the geometry of the setup is not completely symmetrical, and the excitation point could not be perfectly on the centerline of the mass, which means that the response will include some other harmonics despite the axial vibration mode.Furthermore,there will be a larger drift of the displacement signals due to the potential drift of the accelerometers,and the high-frequency noise will be magnified if the signals are integrated on the time domain directly. Therefore, the displacement is calculated on the frequency domain, and the signals are filtered near the first order axial vibration frequency.

The tangential force in the joint is approximated as the product of the free hanging massm2and its acceleration.

Fig.2 Location of accelerometers.

Fig.3 Schematic of the hysteresis curve with the dissipated and elastic energies.

Using the results from Eqs. (1) and (2), a hysteresis curve that depends on the displacement and the transmitted tangential force is evaluated in Fig.3.The damping in the joint could be represented by the loss factor

where, ΔEis the energy dissipated per period of vibration,which equals the area of the hysteresis curve, whileUmaxis the maximum potential energy.The structural damping coefficient ξ=η/2.

The slope of the hysteresis curve is the tangential stiffnesskTof the joint

Based on the existing research, the contact status of the joint interface between two contacting bodies could be stick,micro slip, or macro slip for different loads. While the energy dissipation is caused by the slip, the area of hysteresis loops will decrease with the diminishing of the excitation level. The macro slip of joints in an aero-engine is not allowed to occur.Once the joint parameters are determined experimentally,they can be implemented in the FE model,and the eigenfrequencies and modal damping factors of structures can be predicted.

2.2. Results for structural damping and shear stiffness

The experimental data are collected on the specimen with different torques and excitation levels.The damping of the joint is usually assumed to be frequency-independent.7Therefore, the excitation frequency is set as 360 Hz, and the influence of the frequency is ignored. Fig.4 shows the hysteresis curve over one period for different torques or excitation levels. As we can see,the torques in the joint have influence on both the area and the slope of the hysteresis curves, when the excitation levels are fixed.As the torques increase,the area of the hysteresis curves decreases, which means that there is less energy dissipated per period.But the slope will increase for the tangential stiffness is getting larger when the torque increase.The curve is similar to that in Fig.3,and the loss factor and tangential stiffness can be calculated using Eqs. (3) and (4).

The loss factor and tangential stiffness of joints for different torques and excitation levels are given in Fig.5. The contact area is assumed to be the same for the different torques under a quasi-static state.30The normal contact stiffness obviously increases with the increasing normal load but tends to be constant when the normal load is high. Simultaneously, the tangential stiffness directly proportional to the normal stiffness based on the Hertz contact theory. Therefore, the tangential stiffness increases with the increasing torque,but once the torque exceeds 5 N·m, the influence of the torque is insignificant in this specimen. In contrast, the loss factor decreases due to the reduction of the sliding distance caused by the increasing of the normal pressure. The loss factor increases with the increasing displacement, since the larger sliding distance will lead to higher energy dissipation, while the stiffness declines with a relatively small slope. Notably, the phase difference due to the damping is small, and the amplitude of the relative slip could be approximately equal to the difference between the amplitudes of masses. This conclusion will be useful for the subsequent dynamic analysis of structures with TLE.

The next step is to link the experimental results with an FE model. The examination of the data and the established joint models suggest a nonlinear relation between the relative displacement and the loss factor, but the linear regression shown as Eq. (5) is chosen in the paper for simplicity. When the torque is 2 N·m,the linear regression results of fitting the slopec0of loss factor η to the experimental data is 22.46 /mm. The intercept of the tangential stiffnessk0is 1.31×105N/mm,and the slope coefficient is 4.72×106N/mm2, as shown in Fig.6. The following section will go into more detail about how the linear regression result is used in the updated TLE model.

Fig.4 Hysteresis curve over one period for different torques or excitation levels.

Fig.5 Parameters with different torques and excitation levels.

Fig.6 Linear fit of data with 2 N·m torque.

Fig.7 Hysteresis curve for different assembly (Torque 2 N·m and excitation 50 N).

The experimental data collected in different load periods have shown nonnegligible variability and uncertainty, even though the specimen and the load condition are all the same.For the same specimen, multiple assembly will increase the uncertainty as well,as shown in Fig.7.The probabilistic methods are needed to take the uncertainty into consideration, and the interval analysis method, which only needs the bound of parameters,seems to be more effective,when the statistical distribution of parameters are difficult to be acquired.

3. Thin layer element model for the rotor system with joints

3.1. Model description and FE simulation

3.1.1. Parameters of a thin-layer element

The stiffness of the joint from the experiment should be transferred into the FE model as parameters of the TLEs.As shown in Fig.8, the force acting on both sides of the joint will produce shear stress τ in the TLEs.

whereGis the shear modulus, γ is the shear angle,uis the shear displacement calculated by the difference value between the top and bottom nodes andhis the thickness of the thin layer. The shear stress τ can also be estimated by the ratio of the tangential forceFand the nominal area of contactAarea

By combining Eqs. (6) and (7).

wherekTcan be obtained from the lap joint experiment. The shear modulus could be estimated as

This fact means the shear modulus is directly proportional to the thickness. As previously mentioned, the TLEs are actually solid elements with a special length to thickness ratio,and the constitutive equation of which can be expressed as a diagonal matrix, as in Eq. (10).

E11andE22are zero,for the interface has no principal stresses in the direction parallel to the joint’s surface.Meanwhile,there is no transversal contraction invoked by the contact interface,which meansE12,E13andE23will disappear. Since the joint exhibits no stiffness for in-plane shearing,E44is also zero.7E33define the normal stiffness of the joint, while theE55andE66represents the tangential stiffnessGthat estimated in Eq. (9).

Fig.8 Lap joint with thin layer elements.

The normal elastic modulusE, which depends on the contact stiffness in the normal direction,is assumed to be constant with a higher value. Despite the stiffness, the damping should be modeled in the TLE also.The following equation of motion for rotor systems is used in the calculation of the vibrational response

whereMis a mass matrix,Kis a real-valued stiffness matrix,Cis a damping matrix,Gis a gyroscopic matrix,uis the displacement vector(similar to the Δdin Eq.(10)),andpis the applied load vector. The mass unbalance will produce a specific synchronous force with the imposed circular frequency Ω=2πf,wherefis rotational frequency of rotor.All points in the structure are moving at the same known frequency but not necessarily in phase. The presence of damping will cause phase shifts. The excitation vector and displacements vector can be defined as follows:

where ψ is the force phase shift, φ is the displacement phase shift,umaxis the maximum displacement,andPmaxis the maximum force amplitude.Note that they may be different at each DOF.

The structural damping, which is independent of the frequency, is part ofCand can be written as

where ηjis the loss factor,which is twice the value of the structural damping coefficient,Kjis the stiffness matrix of the TLE with the constitutive equation of stress-strain shown as Eq.(10), andNmis the number of materials for TLEs with different ηj, which are acquired in Section 2 at different load condition (as shown in Figs. 5 and 6).

Thus, the Eq. (13) can be rewritten as

whereKj(1+iηj)is the complex stiffness matrix,andK0is the stiffness matrix exclusive of the TLEs. This equation can be solved for the response with some commercial FE packages,and in this case, the equation was solved using ANSYS. The experimentally determined stiffness and dissipation parameters are used as material properties of the TLEs. When the parameters are constant, the equation could be solved using harmonic analysis. According to the experimental results in Section 2,the joint in the support of a rotor will have different dynamic parameters with different excitation levels, which are proportional to the square of the rotating speed Ω. Thus, the constant value will lead to a slightly incorrect result,especially at the critical speed when the vibration amplitude is high.However,the nonlinear material will lead to a lower computational efficiency. Therefore, the paper will provide an updated TLE model based on harmonic analysis in ANSYS17.0 and iterative algorithm.

3.1.2. Updated TLE model

The transfer function is

whereKjand ηjare the functions of the displacements of the nodes, which decides the relative displacement in the TLE.

The responses of the system can be obtained when Eq.(15)is solved using a direct iterative algorithm.The method can be described by

whereu（r-1） andu（r）,respectively,designate the iteration results of the vectoruat the(r-1)th andrth iteration steps;H-1represents the inverse of the matrix ofH(Ω,u) that is estimated at the (r-1)th iteration step; ε is the predefined error limit and the iteration terminal rule is

The initial values ofGj(0)and ηj(0)are given in ANSYS to calculate responseu（r）; once the relative displacement is acquired, the values ofGj(1)and ηj(1)are obtained using Eq. (5). The TLEs and FE model are updated, and the calculation should be repeated until the results converge. The calculation of the relative displacement is achieved using a subroutine.As previously mentioned,the amplitude of the relative displacements between the top and bottom nodes can be estimated by the difference value of the displacement amplitude if the phase difference is small.

3.2. Application in the rotor system

The bolted flange joint can be modeled in use of the TLEs as shown in Fig.9. The bolts were removed so that the model could be simplified, and the TLEs were inserted between the flanges. The extent of TLE thickness will depend on various factors such as the interface roughness,asperities size and their distribution, particle size and influence of neighboring elements.5Desai et al.4carried out a limited parametric study and suggested that satisfactory simulation of interface behavior can be obtained for thickness to width ratio (h/B) in the range from 0.01 to 0.1, while the ratio is set as 0.05 in this study.

The dynamic effect of joints in the support of the rotor is discussed in this section. As shown in Fig.10, the support and rotor are established by SOLID185, and the elastic modulus is 200 GPa, Poisson’s ratio is 0.3, and the density is 7.8×103kg/m3.The total mass of the rotor is 33.01 kg;during the analysis, the unbalance mass 200 g·mm is applied on the center of the left disk. The MPC184 rigid beam elements are connected to the solid elements by MPC184 weld elements,the beams welded on the support and shaft are connected through MPC184 spherical hinge elements to simulate the ball bearings, while the rear support is simulated by COMBIN14 directly.The model includes 28,285 elements and 34,268 nodes.

Fig.9 Modeling of bolted flange joint using thin layer elements (TLEs).

Fig.10 Schematic of an FE model of a rotor system with TLEs.

The joint in the front support of the rotor system is modeled by TLEs, the parameters of which are chosen based on Fig.6. The thicknesshof the TLEs is about 0.5 mm, when the width ratio (t/B) is set as 0.05. And it is assumed that the connect area has the same parameters as the test piece.The shear modulus in Eq. (10),G=kT·h/Aarea=(134.81–4857.26/mm·u)MPa,whereAarea=485.87 mm2is the contact area of the test piece andkTis the tangential stiffness acquired in Section 2.2. The normal modulusE=200 MPa, and the loss factor η acquired by experiments is (22.46/mm·u). The influence of joints in the rotor is ignored in this example for simplicity, with the assumption that the rotor is in the ideal synchronous forward precession state,and there will be no relative motion in the joint of the rotor.but in practice,the rotor is usually very difficult to operate in the ideal state due to the influence of gravity, unsymmetrical support stiffness and so on, which still need to be investigated in future work.

Notably, for the unbalance excitation with components in thexandydirections,the relative motion between the contact nodes on the interface can be an ellipse, a circle, or a line,depending on the amplitudes and phases, which is totally different from the one-dimensional motion in the test set-up.Usually, a time domain integration method will be used to calculate the steady response to obtain the actual relative motion in the interface.As previously mentioned,if the damping is small enough, the phase difference between the top and bottom nodes of the TLEs can be ignored. Additionally, the elements on the same radial position will have the same shape of relative motion, leading to the same energy dissipation and loss factor (the ellipse of elements at the radial position ofRandrare different,as shown in Fig.11).According to the calculation, the amplitude on one direction is much larger than that of another. Therefore, the contact is simplified to onedimension motion.

The front support stiffness is 3.8×107N/m if the influence of joints is ignored and the material of the TLEs is same as that of the rotor,while the stiffness of the rear support is 4×106N/m. The first critical frequency is 96 Hz corresponding to the rotor mode, but at the second critical speed (257.0 Hz), the vibration of the front support is obvious,as shown in Fig.12.

Fig.11 Relative motion between the top and bottom nodes of the TLEs.

Fig.12 Critical speeds and corresponding modes.

In Fig.13, the dynamic responses of the rotor versus the excitation frequency around the critical speeds are shown,and the vibration amplitude ignoring the influence of joints(without TLE), with linear TLEs and with updated nonlinear TLEs are given. The difference near the first critical speeds is negligible, which could be explained from the mode shape in Fig.12, for there is almost no vibration in front support. In contrast, the influence of joints on the response near the second critical speeds could not be ignored anymore, while the peak frequency reduces to 255.4 Hz, and the amplitude decreases from 9.0 mm to 6.0 mm as the updated nonlinear TLEs are used in the front support joints. Compared to the response of rotor systems with linear TLEs, the rotor using updated TLEs has lower critical speeds, and the amplitude decreases, due to the more accurate stiffness and damping set in updated TLEs for different excitation levels.

The corresponding shear modulus and loss factor when the excitation frequency is near the second peak frequency are shown in Fig.14. The parameters of TLEs can be concluded to have a large variation when the excitation frequency approaches critical speeds of the system. The updated TLEs must be used near the critical speeds to improve the accuracy,while the linear TLEs can be used when the excitation frequency is far away from the critical speeds.

In the actual structure of a bolt joint in an aero-engine,the tangential relative motion is usually limited by the contact.The clearance or the interference is graded in 10-2mm, while the relative displacement near the joint is usually graded in 10-3mm or lower. Therefore, the slip could occur locally, and the influence should be considered.

4. Interval analysis for nonlinear rotor systems

4.1.Nonlinear Chebyshev method for solving the steady response with interval parameters

As previously mentioned, the parameters of the joint show some variability and uncertainty,which can be described using the interval method. The interval parameterxIis used,

where the superscript I denotes the interval character,idenotes the serial number of the uncertain parameter,xicrepresents the mid-value of the parameter, and βiis the uncertain coefficient to define the degree of uncertainty.In this analysis,c0andk0in Eq. (5) are the interval parameters.

The direct iterative algorithm should be revised as an interval form

Fig.13 Diagram of the amplitude versus the excitation frequency.

Fig.14 Diagram of the TLEs’ parameters versus the excitation frequency.

The transfer functionHis a function of the interval parameters. If the interval arithmetic is directly applied in the iterative algorithm, the results will be easily overestimated because of the interval dependency in the interval arithmetic.The existing interval analysis method for nonlinear dynamic systems is usually based on Taylor or Chebyshev expansions and focus on the time domain response of a system. By comparison with the Taylor inclusion function, the Chebyshev inclusion function can effectively control the overestimation.The interval solution is obtained by the combination of the deterministic solutionsuat the interpolation points.

Without loss of generality,the number of uncertain parameters is set asp,and the interpretation point ofm-order Chebyshev inclusion function is calculated using

The responsesuat all the interpolation points are calculated using Eq. (17) until the iteration is stopped. Then, the Chebyshev coefficients are calculated using the Gaussian-Chebyshev integration formula and are given as

The approximation of the responseuwith respect to the uncertain parameters can be given as

whereldenotes the total number of zeros that occur in the subscriptsi1,...,ip. The inclusion function is given as

The upper and lower bounds can be easily obtained by experiments or simulation methods. It can be found that the above analysis process is a non-intrusive interval method that can be used on the FE model with many degrees of freedom in commercial FE software. The parametermis usually not less thann+1. The basic knowledge of Chebyshev polynomials and the derivation of Chebyshev coefficients, which can be found in literatures25, are not given in this paper. In the analysis, the number of uncertain parameterspis 2.

4.2. Example of a Jeffcott rotor

As shown in Fig.13,the dynamic influence of joints is obvious in the rotor system, which also cause the uncertainty of rotor parameters. This section will address a Jeffcott rotor as an example to validate the effectiveness of the interval method proposed in Section 4.1 for the uncertain rotor system.

In Fig.15(a),mdis the mass of the disk,kdis the stiffness of the shaft at the disk,cdis the viscous damping at the disk,mbis the mass of the support,kx=kyare the stiffnesses of the supports, and the viscous damping of the support iscx=cy. The dynamic equation is shown as

wherez=x+iy,zb=xb+iyb, ξ is the eccentricity of the rotor.

By using complex exponents, Eq. (27) can be solved as

The viscous damping in the supports is ignored and the structural damping is considered in Fig.15(b). Eqs. (27) and(28) are revised as

Fig.15 Jeffcott rotor with elastic supports.

The linear response is calculated whenmd=1,mb=0.1,cd=0.05,kd=10,kx=1, ξ=1, and η=0.02. The nonlinear result is calculated when η andkxare valued based on Eq.(5), wherekx=kT,k0=1, coff=0.01,c0=0.50, and Δdis assumed to be the vibration amplitude of the support Δd=abs(zb).

The responses are calculated using an iteration method at different rotating speeds and are shown in Fig.16. The influence rule of the parameters is similar to that presented in Section 3.2. Note that the values of the system have no practical significances.

The dynamic responses of uncertain rotor are calculated by nonlinear Chebyshev interval method, as well as Monte Carlo method, as shown in Figs. 17 and 18. The interval parameters are [c0-βc0,c0+βc0] and [k0-βk0,k0+βk0], where β=0.02.The results show that,the proposed nonlinear interval method can be used for both time domain and frequency domain responses, and the results agree well with those obtained using the Monte Carlo method, which demonstrates the validity of the proposed method.

The uncertain coefficients β has a great influence on the variation range of dynamic response. As β increases from 0.02 to 0.05, the ranges of responses are also increase in both time domain and frequency domain, which mainly caused by the greater variation of stiffness and damping.Note that there is a decline at 0.84 in the lower bound of frequency domain response. As the uncertainty increases, the overestimation could be unavoidable. Next, the nonlinear interval method based on the frequency domain will be used on a complicated FE model.

4.3. Application in the FE model with thin layer elements

Fig.16 Diagram of amplitude and phase of disk versus excitation frequency.

Fig.17 Results of interval analysis and Monte Carlo method when β=0.02.

Fig.18 Results of interval analysis and Monte Carlo method when β=0.05.

The existing interval analysis method is mostly applied on dynamic systems with few degrees of freedom. In the Refs.18,22,the models have no more than a few hundred DOFs.The analysis is performed using MATLAB or some other computing languages rather than commercial FE software due to the intrusive interval analysis algorithm. Once the number of DOFs is large,the efficiency will be incredibly inefficient.Therefore,the method for nonlinear systems proposed in Section 4.1 will have a strong applied value. All the calculations at the interpolating points are finished in ANSYS, while the results are combined to form interval bounds by a post process.In theory,the method could be used on a model of any size.

The double disk rotor system same as that in Section 3.2 is built,as shown in Fig.19,and the same parameters are set.The uncertain coefficients β is set as 0.02,0.05 and 0.20,to consider the variability of joints’parameters.Therefore,the shear modulus of TLEsG=(134.81(1±β)-4857.26/mm·u)MPa,and the loss factor η=(22.46 /mm·(1±β)·u), which are used to model the bolted flange joint in front support.

The Campbell diagram,as shown in Fig.20,shows the variation of the eigenfrequencies that correspond to the mode shapes. Due to the variation of joints’ parameters, the front support stiffness is interval value, and the critical speeds are no longer fixed. The influence of the dynamic parameters on the first order critical mode is relatively small, as shown in Fig.13(a), and the influence of the uncertainty is also small for the same reason. But the 2nd-order critical speed is more sensitive to the variation of parameters, and the 2nd critical speed changes from 255.4 Hz to [250.7, 261.8] Hz when the uncertain coefficients β is set as 0.2.

The interval response of rotor is shown in Fig.21, whenGand η of TLEs are uncertain parameter with uncertain coefficient β=0.2. The Monte Carlo method is employed as well, and the numerical results are in good agreement, which demonstrates the validity of the Chebyshev interval method.The response of rotor is not fixed but an interval value due to the variation of parameters, especially at the frequency closed to the critical speeds when the response is huge. The variation of response near the first critical speeds is negligible for the front support has little influence, but the vibration amplitude closed to the 2nd critical varies in a wide range.Thus, only the dynamic response around the 2nd critical speeds are given out and analyzed below for the conciseness.

The influence on the second order response for different uncertain coefficients is given in Fig.22. When β=0.02, the upper and lower bounds are similar to the results when the parameters are set as the mid-value. The interval grows with the increase of β. The max amplitude of the upper bound increases from 6.06 mm to 7.45 mm, while the peak frequency changes slightly. The max amplitude of the lower bound decreases from 5.83 mm to 3.75 mm. In conclusion, the proposed method, namely, the nonlinear Chebyshev method,can be conveniently applied to the complex system with uncertainty. The joint should be fastened and manufactured with high quality to diminish the uncertainty.

Fig.19 Finite-element model of double disk rotor system considering joints in front support.

Fig.20 Campbell diagram for rotor with interval support stiffness (β=0.2).

5. Conclusions

The paper takes the rotors in aero-engines as the research object, first studies the dynamic characteristics of the joint experimentally and then establishes an updated TLE model,which is applied on the rotor system to simulate joint dynamics.Finally,a nonlinear interval analysis method is established to consider the uncertainty of the parameters and is used in an FE rotor model. It is possible to draw the following conclusions:

Fig.21 Dynamic response of rotor at various rotation speeds (β=0.2).

Fig.22 Bounds of amplitude in different uncertain coefficients.

(1) Experiments on a lap joint are performed to find the damping and stiffness parameters of the joint, and the variation characteristics is verified. The tangential stiffness increases with the increasing torque, while the loss factor decreases due to the reduction of the sliding distance. The loss factor increases with the increasing displacement, while the stiffness declines with a small slope. The results can be fitted using a linear regression.The experimental data collected in different load periods or different surface condition of joints show nonnegligible variability, which need to be analyzed by the probabilistic methods.

(2) Based on the experiments data from a lap joint, the updated TLE method is established to simulate the tangential stiffness and damping of joints for different excitation levels. The application in rotor systems indicated that the proposed method could represent the dynamic properties of joints.The degradation of the stiffness will lead to the left drift of the peak frequency, while the damping decreases the amplitude. The updated model is convenient compared to the contact model and has a better accuracy compared to the linear TLE method.The accurate and rapid response prediction of the structure with joints is achieved.

(3) A nonlinear interval analysis method is established based on Chebyshev polynomials and a direct iteration algorithm to take the uncertainty of joints into consideration.To validate the accuracy of the proposed method,a Jeffcott rotor with elastic supports subjected to an unbalance excitation is demonstrated as an example,and the stiffness and damping of the support are set as uncertain parameters.The presented approach is proved to be effective and accurate by comparison with results from the Monte-Carlo method. The calculation of the interval equations is transferred to the solution of a set of equations with deterministic parameters.The application of the approach on the rotor with joints modeled by updated TLEs indicates that it could be used in a model with the huge amount of DOFs. The results show that,the critical speeds, as well as the dynamic responses of rotor,are no longer fixed,due to the variation of parameters in joints,and the interval grows with the increase of uncertain coefficient β.

In this paper, the tangential stiffness and damping are obtained from the dynamic excitation experiments of a dumbbell test set-up. The joint in an actual structure may suffer complex loads including torsion or bending. Therefore, the parameters of TLEs can be different on the localized position of the interface.The difference should be considered,and some other experiments should be performed in future work. Additionally, the corresponding method to control the overestimation in the interval analysis should be proposed. The method to identify the uncertainty of systems from the uncertain response is another important topic.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (Nos. 51575022, 11772022 and 51475021).