Nonlinear and dual-parameter chaotic vibrations of lumped parameter model in blisk under combined aerodynamic force and varying rotating speed

In this paper, the nonlinear and dual-parameter chaotic vibrations are investigated for the blisk structure with the lumped parameter model under combined the aerodynamic force and varying rotating speed. The varying rotating speed and aerodynamic force are, respectively, simplified to the parametric and external excitations. The nonlinear governing equations of motion for the rotating blisk are established by using Hamilton’s principle. The free vibration and mode localization phenomena are studied for the tuning and mistuning blisks. Due to the mistuning, the periodic characteristics of the blisk structure are destroyed and uniform distribution of the energy is broken. It is found that there is a positive correlation between the mistuning variable and mode localization factor to exhibit the large vibration of the blisk in a certain region. The method of multiple scales is applied to derive four-dimensional averaged equations of the blisk under 1:1 internal and principal parametric resonances. The amplitude–frequency response curves of the blisk are studied, which illustrate the influence of different parameters on the bandwidth and vibration amplitudes of the blisk. Lyapunov exponent, bifurcation diagrams, phase portraits, waveforms and Poincare maps are depicted. The dual-parameter Lyapunov exponents and bifurcation diagrams of the blisk reveal the paths leading to the chaos. The influences of different parameters on the bifurcation and chaotic vibrations are analyzed. The numerical results demonstrate that the parametric and external excitations, rotating speed and damping determine the occurrence of the chaotic vibrations and paths leading to the chaotic vibrations in the blisk.


Introduction
With the development of the aviation industry, the aero-engines are developing toward the light weight, high thrust/weight ratio and high reliability. The blisk is a novel structure of the aero-engine. Due to the coupling effect between the blade and disk sector, we call the blade-disk structure as the blisk. The conventional tenon connections are deleted in the blisk structure, which incorporates the blade and disk as an integrated structure. The blisk can decrease the compressor weight and enhance the thrust/weight ratio of the aero-engines. Comparing with the traditional mortise-tenon joint blade-disk structures, the advantages of the integral blisk structures lead to the widely application in the aero-engines. The blisk structures of the fan and compressors are the effective technologies to improve the performance of the durability and reliability for the aero-engines [1]. However, due to the light and thin structural characteristics and high rotating speed of the blisk, the vibration amplitudes increase rapidly with the change of the aerodynamic force. The thinner the blisk is, the larger the vibration amplitudes of the blisk become, which may even result in the chaotic vibrations. The chaotic vibrations are harmful to the stability and safety of the aero-engines. Therefore, in order to avoid the chaotic vibrations of the aero-engines, we must investigate the chaotic dynamics of the blisk thoroughly. The novelty of this paper is that the nonlinear and dual-parameter chaotic vibrations are investigated for the blisk with the lumped parameter model under combined the varying rotating speed and aerodynamic force. We focus on the dual-parameter chaotic dynamic responses of the blisk. The influences of several parameters on the nonlinear dynamics are mainly analyzed for the blisk structure.
Based on the research results about the free vibration [2,3] and stability analysis [4,5] of the rotary blade, scholars pay more attentions to study the influences of a single parameter on the chaotic vibrations [6,7] of the blade. Especially, when considering the complexity of the operating condition for the blisk structure, the current researches are far from the requirements for analyzing the complex nonlinear dynamics of the blisk. If we want to investigate the chaotic vibrations of the blisk systematically, we must consider the combination of different parameters as much as possible, which may lead to the chaotic vibrations of the blisk. The dual parameters are referred to any combination of the parametric and external excitation, damping and rotating speed in the blisk structure.
The chaos, a seemingly random motion, highly depends on the initial conditions of the nonlinear dynamical systems [8]. When the systems are subjected to the complex excitations, the chaotic vibrations more likely happen in the rotary machines [9]. Researchers devoted more attentions to the nonlinear vibrations of the blades [10][11][12][13][14][15]. However, these studies ignored the influences caused by the interactions among different parameters on nonlinear vibrations of the blade structure. In order to accurately predict and investigate the nonlinear dynamic responses of the blisk structures, it is necessary for us to discuss modal interactions of the blisk structure. The free vibration of the blisk structure is more complex than that of a single blade [16]. Thus, as the basis of the nonlinear coupling vibrations for the blisk structures, the modal analysis is of great importance. Generally, three different methods were used to study the natural vibration characteristics of the blisk structures, namely the reduced-order modeling (ROM) techniques [17][18][19], experiment methods [20,21] and numerical simulations [22,23]. The harmonic balance method [24], eigenmode method [25,26] and parameter identification method [27,28] were used to establish the reduced-order model of the blisk structures. Tang et al. [29,30] presented an extension of the reduced-order model technique based on the component mode mistuning. Sun et al. [31,32] utilized the energy-based finite element method (FEM) to investigate the free vibration of the hardcoating blisk. Zhang et al. [33] studied the vibrations of the rotating pre-twisted composite tapered blade with the graphene coating layers. Zhang et al. [34] investigated the multi-pulse jumping double-parameter chaotic dynamics of the eccentric rotating ring truss antenna under combined the parametric and external excitations. Unfortunately, up till now, almost all ROM methods, including system identification and eigenmode method, are the data-driven empirical models. When the parameters change slightly, the accuracy of the ROM will be greatly reduced. Since the dimension of the equations constructed by the reduced-order model is still too high for the analysis of the chaotic vibration for the blisk structures, few researchers carried out the qualitative investigation about the chaotic vibrations of the blisk structures based on the ROM technique.
In fact, the vibrations of the blisk structures are mostly nonlinear. To understand and solve engineering problems, researchers begin to study the nonlinear vibrations of the blisk structures. Salas et al. [35] used two methods to comparatively study the validity of the noncyclic reduced lumped parameter models and the forced vibrations of the mistuning blisk structures. Considering the effects of the vibration reduction for the damping coating, Yan et al. [36] discussed the forced vibrations of the coated aero-engine blisk. Han and Mignolet [37] predicted the forced dynamic responses of the damped mistuned disk. Willeks et al. [38] studied the response reduction of the bladed disks with the aerodynamic damping by using the harmonic mistuning. Gao et al. [39] gave the application of the hard-coating damper to the mistuned blisk for the passive vibration reduction. Chen [40] investigated the vibration reduction of the blisk by the damping hard coating and gave its mistuning design. Feng and Jing [41] found that at the low frequencies, the vibration isolation of the structures can be remarkably improved by using the nonlinear properties.
Usually, the contact state among the components of the mechanical structures is not consummate. Therefore, the impact of the contact and friction cannot be ignored when studying the nonlinear vibrations of the blisk structures. Laxalde et al. [42] presented the qualitative analysis of the forced dynamic responses for the blisk with the friction ring dampers. Salas et al. [43] gave a mistuned analysis of the forced dynamic responses for an embedded compressor blisk using a reduced-order model. Beirow et al. [44] utilized the equivalent model to explore the influence of the aerodynamic load on the dynamic responses of a blisk. Sarrouy et al. [45] simplified the blades to six identical beams clamped at the disk and explored the complete bifurcation phenomena of the blisk. Hoskoti et al. [46] studied the instability of the vortex induced vibrations for a blade. Gu et al. [47] investigated the free vibration of the rotating cantilever pre-twisted panel with the geometric imperfection of the initial exponential function type. Najafi [48] researched the stability and nonlinear vibrations of a rotating bladed disk at the critical speed.
Although many researches are focused on the nonlinear dynamics of the blades, the studies related to the chaotic dynamics of the lumped parameter model for the blisk structures are less. Researches seem to be blank on the dual-parameter chaotic vibrations of the blisk structures. Considering the coupling influences of various parameters on the nonlinear vibrations of the blisk, the dual-parameter nonlinear vibrations of the lumped parameter model are investigated for the blisk. The governing equations of motion for the lumped parameter model of the blisk are established by applying Hamilton's principle. The influences of different parameters on the bifurcations are analyzed by using numerical simulations. The free vibration and mode localization phenomenon [49,50] of the mistuning blisk are investigated. When the structure is the mistuning, the energy of the system is no longer uniformly distributed. The larger the mode localization factor is, the larger the vibration amplitude is in a certain region. The amplitude-frequency response curves of the blisk are investigated based on the averaged equations in the polar coordinate under 1:1 internal resonance and principal parametric resonance [51]. The dual-parameter Lyapunov exponents and bifurcation diagrams of the blisk are depicted. The effects of the damping, perturbation rotating speed, steady-state rotating speed and external excitation on the nonlinear dynamics of the blisk structure are numerically investigated. The phase portraits, waveforms and Poincare maps are given to demonstrate the nonlinear dynamic responses of the blisk structure.

Free vibration and frequency analysis of blisk
The lumped parameter model of the blisk structure is shown in Fig. 1a. A single sector of two-degree-offreedom lumped parameter blisk structure is simplified as a dual-mass model. This model is developed by introducing the nonlinear part of the stiffness for the beam model of the blade and disk [52]. The introduction of the nonlinear terms for the stiffness makes the nonlinear dynamic analysis of the blisk to be possible. The blade and disk sectors are simplified as the lumped masses and massless nonlinear elastic beams. There are the linear springs between two mass blocks of the simplified blisk structure, which represent the coupling effect between two blades. The equivalent masses of the blade and disk are m b and m d . The stiffness terms of the blade include the linear part k b and nonlinear part e 1 k b . The stiffness terms of the disk include the linear part k d and nonlinear part e 2 k d . The elastic coefficient of the linear coupling spring is k l . The displacements of the blade and disk sectors in the circular direction are assumed to be x i ðtÞ and y i ðtÞ. The rotating speed for the blisk is considered as where X 0 is the steady-state rotating speed and f cos X 1 t is the periodic perturbation rotating speed. Cartesian coordinate system of the blisk rotates around the point o. The unit vectors of the coordinate system are i and j, as shown in Fig. 1b.
The greater the damping is, the faster the dissipation is. Otherwise, the energy dissipates slowly in the blisk structure. In order to keep the free vibration of the blisk structure not attenuating, the damping is not considered in the blisk structure. Without considering the damping, the energy dissipation of the blisk structure is small, and the energy dissipation depends entirely on the characteristics of the springs in the blisk structure.
When the free vibration of the blisk structure is studied without the damping, the mode localization phenomenon of the blades with the stiffness mistuning or mass mistuning in four cases can be studied more directly. When the mode localization factor is large, the vibration energy will concentrate on the specific blades. It is necessary to analyze the nonlinear vibrations of the blisk structure. Thus, we focus on the undamped free vibration for the blisk structure. Therefore, the damping is neglected in Eq. (1).
Without considering the damping, the energy dissipation of the system is small, and the energy dissipation depends entirely on the characteristics of the springs in the system. When the free vibration of the system is studied without damping, the mode localization phenomenon the blades with stiffness mistuning or mass mistuning can be studied more directly. Thus, we focus on the free vibration of the blisk without damping.
It is known from Fig. 1 that the equation of the undamped free vibration for the blisk structure is obtained as follows where M ½ and K ½ are the blisk mass and stiffness matrices, and X ½ is the displacement vector. The matrices M ½ , K ½ and X ½ are expressed as follows ; ð2aÞ ð2bÞ where u is the phase angle between two adjacent blades.
The disk vibrates freely with the frequency x. The following formula is obtained as Equation (2) is expressed by the eigenvalues and eigenvectors Fig. 1 The lumped parameter model of the rotating blisk is given by considering the phase angle, a lumped parameter model of blisk, b rotating Cartesian coordinate system of blisk where A ½ ¼ M ½ À1 K ½ , and I ½ is the identity matrix. The geometric parameters and material properties in Eq. (1)(2)(3)(4) are determined by references [53,54], which are about the continuous modeling of the blade and disk. References [53,54] give the values of the density, Young's modulus and geometry of the blade and disk. For example, according to the relationship among the mass, density and volume, we can calculate the mass of a single blade based on the data in Table 1 of reference [53], namely mb ¼ 7850 kg/m 3 Â 0:1 m Â 0.03 mÂ 0.002 m ¼ 0.0471 kg. In this paper, we select the mass of a single blade to approximately equal to m b ¼ 0.05 kg and disk m d ¼ 0:1 kg.
The stiffness of the structure is related to Young's modulus. Besides Young's modulus of the materials, the stiffness of the structure is also related to the geometric shape of the structure, boundary conditions and external forces. Therefore, we cannot deduce the exact stiffness value based on the data given by in Table 1 in reference [53]. However, we can roughly give the stiffness value of the blade and disk when considering the geometric parameters and Young's modulus of the structure. In material mechanics, the product of Young's modulus and moment of the inertia about the cross section is defined as the flexural stiffness. The unit of the flexural stiffness is Nm 2 . For example, considering b Â h 3 =12 as the moment of the inertia about the cross section of the blade, we derive the flexural stiffness of the blade E Â I z ¼ 200 GPaÂ b Â h 3 =12 ¼ 4 Nm 2 . The unit of the stiffness in our manuscript is N/m. In order to convert Nm 2 to N/m, considering three geometric dimensions of the length, width and height in Table 1 of reference [53], we define the stiffness of the blade as The effects of the randomness on the frequencies of the blades are determined by how the samples are selected. The simulation results may be different due to different selection criteria of the samples. In this paper, we mainly focus on the samples having the normal distributions. It is assumed that the mistuned mass and mistuned stiffness of the blades follow the random normal distributions with the mean value 0. We consider that a and b are defined based on the variance. a is the variance of the random mistuned mass of the blades. b is the variance of the random mistuned stiffness of the blades. The percent sign can be found in the captions of Figs. 4 and 5. The captions   are the variance of the random mistuned mass and variance of the random mistuned stiffness for the blades. The natural frequencies and mode shapes of the blades are investigated. Tables 1 and 2 give the first fifth natural frequencies of the blades and disk sectors without the mistuning, which means that a and b are taken as 0 in this case. Figure 2 demonstrates the first fifteen mode shapes of the blades without the mass mistuning or stiffness mistuning, where the parameter n represents the serial number of the blades. The vibration modes of the blades present the standard sinusoidal manner, which indicates that the distribution of the energy on the blades is completely uniform. Figure 3 exhibits the first fifteen vibration modes of the blades with the mass mistuning, where the variance of the random mistuned mass for the blades is a ¼ 2%. The mass mistuning values Dm b of the blades are indicated in Table 2. In Fig. 3, it is seen that vibration modes of the blades are no longer sinusoidal. However, the mass mistuning destroys the periodic properties of the blisk, which leads to the transfer of the energy between the blades, and finally leads to the mode localization. The aforementioned results illustrate that the stiffness mistuning also leads to the mode localization of the blades.
In this paper, the numbers of the blades and disk sectors remain unchanged. Figure 4 discusses the first twenty natural frequencies of the blades rather than the characteristics of the blisk structure. The horizontal axis of Fig. 4 is the frequency order of the blades, and vertical axis is the frequency values.
When the variance of the random mistuned mass is, respectively, chosen as a ¼ 0, a ¼ 1%, a ¼ 3% and a ¼ 5%, the first twenty natural frequencies of the blades without and with the mass mistuning in four cases are given, as shown in Fig. 4a. It is indicated that the curves of the first eleven natural frequencies for the blades with the mass mistuning a 6 ¼ 0 are smaller than the curve of the frequencies without the mass mistuning a ¼ 0 for the blades. The mass mistuning leads to that the curves of the first eleven natural frequencies for the blades with the mass mistuning a 6 ¼ 0 are lower than the curve of the frequencies without the mass mistuning a ¼ 0 for the blades. However, the mass mistuning of the blades almost has no the effect on higher than the eleventh-order natural frequencies of the blades. When the variance of the random mistuned mass changes from a ¼ 1% to a ¼ 3% or a ¼ 5%, the natural frequencies for lower than the eleventh-order of the blades almost keep consistent in three cases. Figure 4b depicts the first twenty natural frequencies of the stiffness tuning and stiffness mistuning blades when the variance of the random mistuned stiffness are taken as Figure 4b illustrates that the stiffness mistuning b 6 ¼ 0 in the blades results in lower of the first sixteen natural frequencies than one without the mass mistuning b ¼ 0 for the blades. The curve of the natural frequency corresponding to b ¼ 0 is above curves when b 6 ¼ 0 in the blades. However, the natural frequencies of the blades between the sixteenth-order to twentieth-order almost keep consistent. There is no difference in the natural frequencies for higher than the sixteenth-order when b ¼ 0 and b 6 ¼ 0.
Using Euclidean norm to define the mode localization factors, we have Fig. 4 The natural frequencies of the blades are obtained without the mistuning and with mistuning, a influence of masses without mistuning and with mistuning when a ¼ 0, a ¼ 1%, a ¼ 3% and a ¼ 5%, b influence of stiffness terms without mistuning and with mistuning when where A max is the largest vibration amplitude of the blades, n is the quantity of the blades, j is the number of the blade with the largest amplitude, and A i is the amplitude of the ith blade. Figure 5a demonstrates the mode localization factors of the blades without and with the mass mistuning in four cases when the variance of the random mistuned mass is chosen as a ¼ 0, a ¼ 1% and a ¼ 3%, respectively. Figure 5b illustrates the mode localization factors of the blades without and with the stiffness mistuning when the variance of the random mistuned stiffness, respectively, is The horizontal axis and vertical axis in Fig. 5 are the serial number of the blades and mode localization factor, respectively. Based on Fig. 5, the mode localization factors are less than or equal to 1 when the blades vibrate without mistuning. By contrary, when the blades vibrate with the mass mistuning or stiffness mistuning, the mode localization factors are greater than or equal to 1. When the mode shapes of the blades are sinusoidal or the mode localization factors are less than 1, the mode localization phenomenon will not appear. Correspondingly, when the mode shapes of the blades do not distribute sinusoidally or the mode localization factors are greater than 1, the mode localization phenomenon occurs. The larger a or b are, the larger the mode localization factors of the blades become.
We use the modal localization factor to study the modal localization phenomenon of the blisk structure. The mode localization factor is defined by Euclidean norm, as shown in Eq. (5). When the free vibration of the blisk structure is studied without the damping, the mode localization phenomenon of the blades with the stiffness mistuning or mass mistuning in four cases can be studied more directly. When the mode localization factor is large, the vibration energy will concentrate on the specific blades. It is necessary to analyze the nonlinear vibrations of the blisk structure. Thus, we focus on the undamped free vibration of the blisk structure. Therefore, the damping is neglected in Eq. (1).
As a combination, the small damping of the blisk structure will lead to the large amplitude nonlinear vibrations. Thus, it is urgent for us to study the nonlinear dynamical behaviors of the blisk.

Nonlinear governing equations of motion
Based on the lumped parameter model of the blisk in Fig. 1, we establish the dynamic model of the blisk. The displacements of the blades m b and disks m d are expressed as The velocity and acceleration of the blade and disk are obtained as follows and Fig. 5 The mode localization factors of the blisk are obtained without mistuning and with mistuning, a influence of masses without mistuning and with mistuning when a ¼ 0, a ¼ 1%, a ¼ 3% and a ¼ 5%, b influence of stiffness terms without mistuning and with mistuning when The variation of the kinetic energy is given by The potential energy U consists of the elastic potential energy U 1 of the blades, elastic potential energy U 2 of the disks and elastic potential energy U 3 of the linear springs. These potential energy expressions are given as follows ð10cÞ The variation of the potential energy is obtained as The mass blocks of the blade are subjected to the aerodynamic force. Therefore, we have The variation of the virtual work done by the first type of the aerodynamic force is denoted as follows The differential governing equation of motion for the rotating blisk is established by adopting Hamilton where t is time and d is the variation.
Substituting dK, dU and dW into Eq. (14), the differential governing equations of motion for the rotating blisk are derived in terms of the displacements. dx i : where l 1 and l 2 are the structural damping terms of the blisk. The displacement relationships between two adjacent blades and two adjacent disk sectors are described as follows where u is the phase angle between two adjacent blades, and v is the phase angle between two adjacent disk sectors. We introduce the following dimensionless transformation Substituting Eqs. (16) and (17) into Eq. (15), for convenience, we remove the '-' above the parameters, two-degree-of-freedom dimensionless differential governing equations of motion for the blisk are obtained as where m is the total number of the blades. Substituting X, _ X and e ix ¼ cos x þ i sin x into Eq. (18) and separating the real and imaginary parts, the nonlinear governing equations of motion for describing the transverse vibrations of the rotating blisk are obtained as follows where the coefficients of Eq. (19) are calculated based on simplifying equations obtained by separating the real part and imaginary part of Eq. (18). We introduce the following transformations Equation (19) is transformed into the first-order ordinary differential equations. Therefore, we obtain the dimensionless first-order ordinary differential governing equations of motion for the blisk 4 Perturbation analysis and amplitude-frequency curves The perturbation analysis is performed by using the method of multiple scales [53]. The dimensionless transformations are introduced as follows where x 2 1 ¼ a 13 À a 12 X 2 0 and x 2 2 ¼ a 23 À a 22 X 2 0 . For a 12 , a 12 , a 12 and a 23 , the parameter transformations are not introduced, which can be verified from Eq. (23). The method of multiple scales is adopted to obtain averaged equations. We let that the linear solution of Eq. (23) be in the following form where T 0 ¼ t and T 1 ¼ et.
The resonant phenomena easily happen for the blade and disk. Due to the energy transfer between different modes, the relationship between two natural frequencies leads to the large amplitude nonlinear vibrations. It is found from Tables 1 and 2 that comparing the second-order and first-order natural frequencies of the blade and disk sectors, we obtain the conclusion in which the blades and disk sectors can exhibit 1:1 internal resonance. Considering 1:1 internal resonance, principal parametric resonance and primary resonance, the resonant relations are presented as follows. (25). where x 1 and x 2 , respectively, are the second-order and first-order natural frequencies of the blade and disk sector, r 1 and r 2 are two detuning parameters, X 1 is the frequency of the periodic perturbation rotating speed, X 2 is the frequency of the aerodynamic force.
Thus, the derivatives with respect to t become.
Assuming X 2 ¼ 1, substituting Eqs. (24)-(26) into Eq. (23) and balancing coefficients of the similar power of e yield. e 0 order: e 1 order: The general solution of Eq. (27) is expressed as where A 1 and A 2 are the complex conjugate terms of A 1 and A 2 . The coefficients of x 10 and x 20 in Eq. (27) are the natural frequencies of the rotating blisk. Tables 1 and 2 give the first 15 natural frequencies of the blades and disk sectors in Hertz. It is seen from Tables 1 and 2 that the blades and disk sectors exhibit 1:1 internal resonance. Substituting Eq. (29) into Eq. (28) yields :: where cc denotes the complex conjugate of the right terms in Eq. (30), and NST is the non-secular part.
It is no doubt that there must exist the terms related to e iT 0 , e i2T 0 and so on in Eq. (30). This leads to that r 1 appears in the coefficients of e iT 0 in Eq. (30a). r 2 appears in the coefficients of e iT 0 in Eq. (30b). The coefficients of e iT 0 are called as the secular term. The coefficients of e i2T 0 and so on are called as the nonsecular term. Based on the method of multiple scales, these secular terms must equal to zero. In the right of Eq. (30), all terms except cc and NST are related to e iT 0 . Obviously, Eq. (30) is a non-autonomous dynamical system. Letting the secular terms in Eq. (30) to be zero and eliminating the secular terms, the following equations are derived as To obtain the averaged equations in the polar coordinate, A 1 and A 2 are expressed as follows To obtain the averaged equations in the polar coordinate system, A 1 and A 2 in Eq. (31) are presented by the amplitudes and phases, as shown in Eq. (32). Substituting Eq. (32) into Eq. (31) and separating the real part and imaginary part, the average equations in the polar coordinate system are derived as follows Based on Eq. (33), the amplitude-frequency response curves are investigated for the rotating blisk. Some parameters obviously affect the nonlinear dynamic behaviors of the lumped parameter model for the blisk. We mainly study influences of the rotating speed, aerodynamic force, damping, coupling stiffness and cubic nonlinear terms on the amplitudefrequency response curves of the blisk (Table 3).
Since the dimensionless parameters are introduced in Eq. (22), Eq. (23) are dimensionless. For simplicity, we take X 2 ¼ 1. In other references [6,10], we also take that a frequency is 1. This is the classical research method for the nonlinear dynamical systems. In Eq. (34), the frequencies are represented by r 1 and r 2 . System (34) is a nonlinear algebraic equation on the frequencies r 1 and r 2 . In Eq. (34), the natural frequencies are fixed to 1.
After Eq. (21), numerical simulations on the amplitude-frequency curves and dual-parameter chaos of the lumped parameter model for the rotating blisk are performed. After the internal resonance and multiscale perturbation analyses, due to the complexity of the derivation process, we cannot obtain the corresponding mathematical relationship between two groups of parameters. Therefore, we give the independent parameters in numerical simulations of the amplitude-frequency response curves and dual-parameter chaos. We study the influences of different parameters on the nonlinear vibrations of the rotating blisk and focus on revealing the qualitative influences of the relevant parameters. Even if two sets of parameters are completely independent, our work is still of the certain significance for further understanding the nonlinear vibrations of the lumped parameter model for the blisk structure. The equations of the steady-state vibrations for the blisk are sixth-order nonlinear algebraic equations with four variables. Therefore, it is impossible for us to analytically solve sixth-order nonlinear algebraic equations with four variables. We utilize some assumptions to analyze this equation. The amplitude-frequency response curves are related to the amplitudes and frequencies. Numerical simulations of the amplitude-frequency curves are carried out based on the dimensionless Eq. (34). Theoretically, if the relationship between the amplitudes a 1 , a 2 and frequencies r 1 , r 2 can be solved separately, the amplitude-frequency response curves can be drawn directly. However, both amplitudes a 1 and a 2 appear in Eqs. (34a) and (34b), which mean a 1 and a 2 are highly coupling. This makes it impossible for us to express the amplitudes a 1 and a 2 with the frequencies r 1 , r 2 directly. Therefore, we focus on a special case of solving the amplitude-frequency response curves. Considering the weak coupling condition, namely, letting a 2 ¼ 1 in Eq. (34a), the relationship between the amplitude a 1 and frequency r 1 can be solved. Similarly, letting a 1 ¼ 1 in Eq. (34b), the amplitude a 2 can be expressed by the frequency r 2 . In this case, the amplitudes change when the frequencies vary. It is noticed that there are two possible cases, namely a 1 is nonzero and a 2 is zero, or a 1 and a 2 are nonzero. In reference [55,56], the authors discussed the amplitude-frequency response curves of the system when a 1 is nonzero and a 2 is zero. We analyze the special case, namely a 1 is nonzero and a 2 ¼ 1 in Eq. (34a), together with a 1 ¼ 1 and a 2 is nonzero in Eq. (34b).
Letting a 2 ¼ 1 in Eq. (34a) and a 1 ¼ 1 in Eq. (34b), the amplitude-frequency response curves of the blades and disk subjected to the aerodynamic force are obtained, as shown in Figs. 6 and 7. In Figs. 6 and 7, other parameters remain unchanged except for the controlling parameters. The typical Fig. 6 The amplitude-frequency response curves are obtained for the blades, a influence of damping, b influence of amplitude of perturbation rotating speed, c influence of coupling stiffness parameter, d influence of aerodynamic force, e influence of cubic term Fig. 7 The amplitude-frequency response curves are obtained for the disk, a influence of damping, b influence of perturbation rotating speed, c influence of coupling stiffness parameter multi-valued phenomena occur, which results in the jumping of the solution in Eq. (34). Figure 6a demonstrates that the amplitudes and bandwidths of the blades decrease with the increase of the damping l 1 . The increase in the perturbation rotating speed f and coupling stiffness k 12 result in the increase in the amplitude and bandwidth for the blades, as shown in Fig. 6b, c. The right-ward shift of the amplitude-frequency curve occurs for the blades, as shown in Fig. 6b. Figure 6d illustrates that as the aerodynamic force p 1 increases, the amplitudes of the amplitude-frequency response curves for the blades increase, but the bandwidth deceases. From Fig. 6a-d, it is shown that the blades have the hard spring properties. As the frequency r 1 changes, the multivalued phenomena occur. When the cubic term a 15 j j increases, Fig. 6e depicts that the amplitude-frequency response curves bend and are away from zero solution, and the amplitudes of the dynamic responses decreases. When the nonlinear parameter a 15 is positive, the blades exhibit the hard spring properties. However, when the nonlinear parameter a 15 is negative, the blades have the soft spring properties. When value of a 15 is opposite, the amplitude-frequency response curves are symmetric. Figure 7 illustrates the influences of several parameters on the amplitude-frequency response curves of the disk, including the damping l 2 , perturbation rotating speed f and coupling stiffness k 21 . It is illustrated that the disk has the hard spring properties. The results show that damping l 2 and coupling stiffness k 21 have almost the same effect on the amplitude-frequency curves for the disk, as shown in Fig. 7a, c. In Fig. 7a, the damping l 2 has no difference to the response bandwidth of the disk. Figure 7b demonstrates that the right-ward shift or left-ward shift of the amplitude-frequency response curves for the disk does not occur with the change of the perturbation rotating speed f . It is pointed out that a smaller f results in the remarkable narrowing in the response bandwidth and the enhancement of the amplitude for the disk.
The internal resonance relationship is determined by the frequency relationship. Assuming one amplitude being a constant value and solving another amplitude, this has no effect on the internal resonance relationship. The internal resonance causes the decrease or increase of the amplitudes. When a 2 ¼ 1, we obtain Fig. 6, in which a 1 varies with the detuning parameter r 1 . When a 1 ¼ 1, we obtain Fig. 7, in which a 2 varies with the detuning parameter r 2 .

Numerical simulation of dual-parameter chaos
In recent decades, a large variety of numerical methods have been developed to detect the periodic and chaotic vibrations of the nonlinear dynamical systems. The bifurcation diagrams and Lyapunov exponents have been the common choices for investigating the chaotic vibrations of the structures with the geometric nonlinearity and subjected to various kinds of excitations. Bifurcation diagrams are the effective tools detecting the periodic and chaotic vibrations with a single controlling parameter. However, only depending on two-dimensional bifurcation diagrams is insufficient to illustrate the existence of the chaotic dynamics for the dual-parameter nonlinear system in three-dimensional (3D) space. Lyapunov exponents are used to distinguish the periodic and chaotic vibrations. However, it is difficult to identify which type of the periodic vibration appears in the nonlinear dynamical system according to Lyapunov exponents. Therefore, to exceed the limitations of the traditional bifurcation diagrams and Lyapunov exponent, the dual-parameter bifurcation diagrams are considered to study the periodic and chaotic vibrations of the structure.  Figure 8 shows 3D dual-parameter Lyapunov exponents and bifurcation diagrams of the blisk versus the perturbation rotating speed f and aerodynamic force p 0 . Obviously, there is an excellent correspondence between Lyapunov exponents and bifurcation diagrams. Based on 3D dual-parameter Lyapunov exponents, we give the influences of one parameter on the nonlinear dynamics of the blisk structure when another parameter is regarded as a controlling parameter, as shown in Figs. 9. These figures reveal clearly how the system evolves between the periodic and chaotic vibrations.
The dimensionless analyses of the nonlinear dynamic behaviors are carried out to avoid being confused in the physical parameters for the blisk structure. In order to analyze the chaotic vibrations of the blisk structure, the selections of the dimensionless initial conditions and parameters are not related to the physical parameters. Thus, the abundant nonlinear dynamic phenomena can be revealed for the blisk structure. The any changes of the parameters will affect the nonlinear vibration behaviors of the blisk structure. The relations of these parameters to the mechanical properties are not presented in the blisk structure. We pay attention to the qualitative analyses on the influences of the dimensionless parameter changes on the nonlinear vibration characteristics of the blisk structure. These provide the theoretical reference for the design of the mechanical structures.

Effect of Perturbation Rotating Speed
Considering the effect of the perturbation rotating speed f and aerodynamic force p 0 , the dual-parameter chaotic vibrations are studied for the blisk, and the paths to chaos are explored. Based on 3D dualparameter bifurcation diagrams and Lyapunov Fig. 12 The chaotic vibrations of the lumped parameter model are studied for the blisk when f ¼ 72, a phase portrait on plane ðx 1 ; x 2 Þ, b waveform on plane ðt ; x 1 Þ, c phase portrait on plane ðx 3 ; x 4 Þ; d waveform on plane ðt ; x 3 Þ, e 3D phase portrait in space ðx 1 ; x 2 ; x 3 Þ; f Poincare map exponents, single parameter bifurcation phenomena of the blisk are obtained by using the variable controlling method. Selecting the perturbation rotating speed f and aerodynamic force p 0 as two controlling parameters, the dual-parameter Lyapunov exponents exhibit the regions where the dual-parameter chaotic vibrations occur, as shown in Figs. 8a, b, where damping terms l 1 and l 2 , and the steady-state rotating speed X 0 are given as l 1 ¼ 0:46, l 2 ¼ 0:46, X 0 ¼ 200. In order to show the evolution process of the periodic and chaotic vibrations, 3D dual-parameter bifurcation diagram is given in Fig. 8c.
In Fig. 8a, two horizontal axes, respectively, are the perturbation rotating speed f and aerodynamic force p 0 , and the vertical axis is the dual-parameter Lyapunov exponents. We calculate four Lyapunov exponents of each point. Therefore, there are four surfaces along the vertical axis in Fig. 8a. According to the legend on the right side of Fig. 8a, Lyapunov exponents in red are greater than zero, which indicates that the blisk performs the dual-parameter chaotic vibrations under combined the perturbation rotating speed f and aerodynamic force p 0 when the parameters are located at the regions in red. The striking red regions in Fig. 8a, b are called as the chaotic regions. To achieve the transformation from one chaotic region to next chaotic region, the periodic vibration between two chaotic regions must happen, where Lyapunov exponents are less than or equal to 0. Meanwhile, 3D dual-parameter bifurcation diagram in Fig. 8c is given to further verify that this phenomenon occurs in Fig. 8a, b. The results demonstrate that these figures correspond to each other perfectly.
The evolution paths of the dual-parameter chaotic vibrations are studied in Fig. 9 when p 0 is located in the interval 0; 100 ð Þ. Choosing the aerodynamic force p 0 as the controlling parameter, the effect of the perturbation rotating speed f on the dynamics of the blisk is analyzed. Figure 9a presents the dual-parameter Lyapunov exponents.
perturbation rotating speeds f increase, the periodic vibration of the blisk turns to the chaotic vibrations gradually, as shown in Fig. 9b-d. When f ¼ 2, the points of the bifurcation diagram in Fig. 9b are sparse, which means that the blisk only has the periodic vibration. When f ¼ 10, the intensive points appear in Fig. 9c, which implies that the chaotic vibrations arise in the blisk. When f increases to f ¼ 25, the intensive points occupy the larger range, which means that the chaotic vibrations occur in a wider interval of the aerodynamic force p 0 , as shown in Fig. 9d. The amplitudes of the chaotic vibrations increase. Moreover, the larger the rotating speed is, the larger the amplitude of the chaotic vibrations is. The results indicate that the perturbation rotating speed f decides whether the chaotic vibrations of the blisk occur under combined the perturbation rotating speed f and aerodynamic force p 0 when p 0 is in the interval 0; 100 ð Þ. Three-dimensional bifurcation diagrams of the perturbation rotating speed f versus the blade displacement x 1 and disk displacement x 3 are presented in the interval 0; 100 ð Þ, as shown in Fig. 10. Where the aerodynamic force p 0 is p 0 ¼ 50, and other parameters remain the same as these given in Fig. 8. Projecting 3D bifurcation diagram, namely, the black part in Fig. 10a, to the planes f ; x 1 ð Þand x 1 ; x 2 ð Þ, the twodimensional bifurcation diagram and Poincare map are obtained. Two-dimensional bifurcation diagram and corresponding Poincare map for x 3 are given, as shown in Fig. 10b. Due to figures obtained by numerical simulations, an accurate judgment on the nonlinear vibrations of the blisk is realized.
The enlarged version of the bifurcation diagrams for the perturbation rotating speed f versus the blade displacement x 1 and disk displacement x 3 are, respectively, displayed in Fig. 11a, b. Figure 11c depicts Lyapunov exponents of the rotating blisk when considering the perturbation rotating speed f as a variable. Lyapunov exponents and bifurcation diagrams have high consistency. Lyapunov exponents generated by the system can be positive, negative or zero. The positive Lyapunov exponents imply that the system must be the chaotic or hyperchaotic vibrations. If a chaotic system has the multiple positive exponents in high-dimensional phase space, it is hyperchaotic. Fig. 15 The dual-parameter Lyapunov exponents and bifurcation diagrams of X 0 and p 0 are obtained, a dualparameter Lyapunov exponents on space ðX 0 ; p 0 ; LEs Þ, b projection of dual-parameter Lyapunov exponents on plane ðX 0 ; p 0 Þ, c dual-parameter bifurcation diagram on space ðX 0 ; p 0 ; x 1 Þ It is seen from Fig. 11 that the periodic and chaotic vibrations of the rotating blisk appear repeatedly. The numerical results demonstrate that the perturbation rotating speed f has a major effect on the nonlinear dynamics of the blisk. Both negative Lyapunov exponents in Fig. 11c and sparse points in Fig. 11a, b illustrate the periodic vibrations occur in the blisk. The high density points in the interval f ¼ 28 À 65 indicate that the chaotic vibrations appear, this is proved by Lyapunov exponents in the corresponding interval in Fig. 11c. When the perturbation rotating speed f is f ¼ 72, the max Lyapunov exponent of the blisk is positive, as shown in Fig. 11c. The nonlinear vibrations of the large amplitude appear when f ¼ 72, as shown in Fig. 12 Figure 12e is 3D phase portrait in space Þ . The blue and black parts in Fig. 12f Fig . 16 The dual-parameter Lyapunov exponents and five bifurcation diagrams are given with X 0 as the control parameter, a dual-parameter Lyapunov exponents on plane ðX 0 ; p 0 Þ, b bifurcation diagram on plane ðp 0 ; In the aforementioned analysis, the phase portraits, waveforms and Poincare map are highly consistent.

Effects of external excitation and steady-state rotating speed
The dynamic behaviors of the rotating blisk under the aerodynamic force p 0 are presented in Figs. 13 and 14 when the perturbation rotating speed f is f ¼ 70. Other parameters keep unchanged. It is seen obviously from Fig. 14 that the chaotic vibrations mainly occur in two intervals, in which the plenty of points are presented, and there are positive of Lyapunov exponents. It is hard to decide whether there are the periodic vibrations in these two intervals depending on Fig. 14a, b. However, Lyapunov exponents in Fig. 14c illustrates that there are exactly the period-1, multi-periodic and chaotic vibrations with the increase of the aerodynamic force p 0 . In the same way, the dual-parameter chaotic vibrations are investigated for the blisk under combined the steady-state rotating speed X 0 and aerodynamic force p 0 when the perturbation rotating speed f is f ¼ 70, as shown in Fig. 15. The initial condition and other parameters keep the same as Fig. 8. Figure 15a is 3D dual-parameter Lyapunov exponents of the blisk. Projecting Fig. 15a on the horizontal plane, we obtain Fig. 15b. Figure 15 is the corresponding 3D dual-parameter bifurcation diagram for the blisk structure. These diagrams correspond to each other completely. Figure 15b shows that the left side of the dualparameter chaotic region extends along 45 degrees and there is an obvious boundary between the chaotic and periodic regions. To study the relation between the boundary and dual-parameter chaotic vibrations of the Fig. 18 The bifurcation diagrams and Lyapunov exponents of the steady-state rotating speed for the blisk are given in the interval X 0 ¼ 160 $ 260, a bifurcation diagram on plane ðx 1 ; X 0 Þ, b bifurcation diagram on plane ðx 3 ; X 0 Þ, c Lyapunov exponents Fig. 17 3D bifurcation diagrams of the steady-state rotating speed for the blisk are given in the interval X 0 ¼ 160 $ 260, a 3D bifurcation diagram on space ðx 1 ; x 2 ; X 0 Þ, b 3D bifurcation diagram on space ðx 3 ; x 4 ; X 0 Þ blisk under combined the steady-state rotating speed X 0 and aerodynamic force p 0 , we obtain Fig. 16. Figure 16b-f is five bifurcation diagrams about the aerodynamic force when the steady-state rotating speeds, respectively, are X 0 ¼ 180, 190, 195, 210 and 220. According to the motion of the vertical red line in these diagrams, we see clearly that the chaotic region moves toward right along the horizontal axis p 0 with almost the same maximum vibration amplitude.
Single period vibration only appears on the right side in Fig. 16b, d. Figure 16e, f illustrates when X 0 increases to a certain value, the vibration of the blisk becomes more complex and the periodic vibrations alternately appear with the chaotic vibrations, which indicates the aerodynamic load p 0 has an important effect on the dual-parameter chaotic dynamics for the blisk. Fig. 19 The chaotic vibrations of the lumped parameter model for the rotating blisk are studied when X 0 ¼ 225, a phase portrait on plane ðx 1 ; _ x 1 Þ, b waveform on plane ðt ; x 1 Þ, c phase portrait on plane ðx 3 ; _ x 3 Þ, d waveform on plane ðt ; x 3 Þ, e 3D phase portrait in space The effect of the steady-state rotating speed on the nonlinear dynamics of the blisk is studied when the aerodynamic force is p 0 ¼ 50 and other parameters are the same as Fig. 15. 3D bifurcation diagrams of the steady-state rotating speed X 0 for the blisk in the interval X 0 ¼ 160 $ 260 are given, as shown in Fig. 17. Two-dimensional bifurcation and Poincare maps are obtained by projecting 3D bifurcation diagram to the specific planes. Two-dimensional bifurcation diagrams are shown using the blue and red. Poincare maps are described in the purple. The bifurcation diagrams of X 0 versus the x 1 and x 3 are represented, as shown in Figs. 18a, b. Figure 18c gives Lyapunov exponents of the blisk. The numerical results illustrate that the periodic and chaotic vibrations appear. When the steady-state rotating speed of the blisk is taken as X 0 ¼ 225, the chaotic vibrations are obtained, as shown in Fig. 19.

Effect of damping on nonlinear dynamics
In order to explain the effect of the damping on the dual-parameter vibrations of the blisk, three-dimensional dual-parameter Lyapunov exponents, two-dimensional dual-parameter Lyapunov exponents and three-dimensional dual-parameter bifurcation diagram of the blisk with two damping terms l 1 and l 2 are, respectively, provided, as shown in Fig. 20. In this case, the aerodynamic force p 0 , the steady-state rotating speed X 0 and the perturbation rotating speed f are, respectively, selected as p 0 ¼ 50, X 0 ¼ 200 and f ¼ 20. Four layers in Fig. 20a represent four groups of Lyapunov exponents for the blisk when l 1 and l 2 are considered as the controlling parameters. The colors in Fig. 20a indicate that the blisk presents complex nonlinear dynamic behaviors. The four quadrants in Fig. 20, c are called 1, 2, 3 and 4, respectively. It is observed that Fig. 20b has good agreement with Fig. 20c.
When the damping terms l 1 and l 2 are limited to quadrant 1, the blisk only performs the periodic vibrations. Otherwise, the chaotic vibrations happen. The blisk is most likely to perform the chaotic vibrations with a small damping. This fully explains why in quadrant 3 has the largest proportion of the red region among four quadrants. The dense points in Fig. 20c show that the blisk performs the chaotic vibrations. Fig. 20 The dual-parameter Lyapunov exponents and bifurcation diagrams of l 1 and l 2 are given for the blisk, a dual-parameter Lyapunov exponents on space ðl 1 ; l 2 ; LEsÞ, b projection of dual-parameter Lyapunov exponents on plane ðl 1 ; l 2 Þ, c dualparameter bifurcation diagram on space ðl 1 ; l 2 ; x 1 Þ In Fig. 20, five bifurcation diagrams of the blisk are obtained when the damping terms are l 2 ¼ 0:1, 0:15, 0:2, 0:3 and 1:0. To reveal the variation rule between the chaotic phenomena and bifurcation, Fig. 21 is given. Four straight lines with the arrows are marked along the boundary of the chaotic regions in Fig. 21a. Figure 21(b-d demonstrates that the chaotic regions move to the left as l 2 increases in a certain range. The nonlinear vibrations of the blisk give the regular changes: period-1 ? multi-periodic ? chaotic vibrations when l 1 is selected in the interval between the red and black lines. There are the period-2, multiperiodic and chaotic vibrations in the region between the green and black lines in Fig. 21d. With l 2 becomes larger in the quadrants 3 and 4, the nonlinear vibrations of the blisk become more and more diversification. The periodic vibrations occupy a wider interval, which means that the blisk is less likely to perform the large amplitude nonlinear vibrations when l 2 is large enough.
The influences of the damping l 1 on the nonlinear vibrations of the blisk are presented, as shown in Fig. 21 The dual-parameter Lyapunov exponents and five bifurcation diagrams are obtained for the blisk with the damping l 2 as the control parameter, a dualparameter Lyapunov exponents on plane ðl 1 ; l 2 Þ, b bifurcation diagram on plane ðl 1 ; x 1 Þ when l 2 ¼ 0:1, c bifurcation diagram on plane ðl 1 ; x 1 Þ when l 2 ¼ 0:15, d bifurcation diagram on plane ðl 1 ; x 1 Þ when l 2 ¼ 0:2, e bifurcation diagram on plane ðl 1 ; x 1 Þ when l 2 ¼ 0:3, f bifurcation diagram on plane ðl 1 ; x 1 Þ when l 2 ¼ 1:0 Fig. 22. The bifurcation diagrams and Lyapunov exponents of the blisk are depicted. The damping l 2 is l 2 ¼ 0:46. Other parameters are as same as Fig. 20. It is found that the dynamic responses of the blisk exhibit the following evolution law: chaotic ? periodic vibrations, with the increase of the damping l 1 . It is seen that the periodic vibration happens when the damping l 1 increases to l 1 ¼ 0:44. The damping l 1 has significantly effects on the nonlinear dynamics of the blisk.

Conclusions
The dual-parameter chaotic vibrations of the lumped parameter model for the rotating blisk are studied under combined the varying rotating speed and aerodynamic force. The free vibration and mode localization phenomenon of the lumped parameter model are investigated for the tuning and mistuning blisk structure. The nonlinear dynamic model of the rotating blisk subjected to the aerodynamic force is presented by using Hamilton's principle. The amplitude-frequency response curves are obtained by applying the method of multiple scales under 1:1 internal resonance and principal parametric resonance. The evolution process to the chaotic vibrations for the rotating blisk is described by the dual-parameter Lyapunov exponents and bifurcation diagrams. The influences of different parameters on the nonlinear responses of the blisk are analyzed based on the amplitude-frequency response curves of the first two modes. The main conclusions are obtained.
1. Comparing with the mode shapes of the tuning blades, the mode shapes of the mass mistuning or stiffness mistuning blades are not sinusoidal.
There is a positive correlation between the mistuning variable and mode localization factor. The larger the mode localization factor is, the more irregular the mode shape becomes. The large amplitude vibrations appear in a certain part and lead to the chaotic vibrations of the blisk. 2. Some parameters seriously affect the vibration amplitudes, response bandwidths and stiffness properties of the blisk. The spring property of the blisk remains unchanged when the perturbation rotating speed, coupling stiffness and aerodynamic force are changed. The nonlinear terms determine the spring property of the blisk. If it is positive, the blisk shows the hard spring property. When it becomes negative, the blisk presents the soft spring property. The corresponding amplitude-frequency curves are perfectly symmetrical about the longitudinal axis as the values of the nonlinear term are opposite. 3. The dual-parameter Lyapunov exponents and bifurcation diagrams of the blisk illustrate that the evolution paths leading to the chaotic vibrations and varying pattern of the chaotic regions are determined three different parameters of the blisk, including the rotating speed, aerodynamic force and damping. 4. The phase portraits, waveforms and Poincare maps are depicted to verify the nonlinear dynamic behaviors of the blisk according to the dualparameter Lyapunov exponents and bifurcation The aforementioned conclusions indicate that chaotic vibrations are closely related to the comprehensive effect of the various factors for the blisk. Due to the small damping and approximate stiffness of the disk and blade, the large amplitude coupling vibrations of the blisk structure are more likely to occur, and even cause a disaster. Thus, our study has practical significance for the nonlinear vibration control of the blisk structure.