Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 56, 2, pp. 549-566, Warsaw 2018 DOI: 10.15632/jtam-pl.56.2.549 VIBRATION OF A MISTUNED THREE-BLADED ROTOR UNDER REGULAR AND CHAOTIC EXCITATIONS JERZY WARMINSKI, Jarosław Latalski, Zofia Szmit Lublin University of Technology, Department of Applied Mechanics, Lublin, Poland e-mail: j.warminski@pollub.pl This study considers forced vibrations of a rotating structure consisting of a rigid hub and three flexible beams.The blades are nominallymade of amultilayered laminatewith a speci- fic stacking sequence resulting in full isotropicmacroscopicmaterial behaviour. However, in the performedanalysis it is assumed that the rotor hasbeenmistunedbecause ofmanufactu- ring tolerances of the composite material. These inaccuracies are represented by deviations of reinforcing fibres orientations from their nominal values. The considered tolerances break the intendedmacroscopicmaterial isotropyandmake the laminate to exhibit the fully ortho- tropic behaviour. Based on previous authors research, the system of four mutually coupled dimensionless ordinary differential governing equations is adopted. Forced responses of the system under regular and chaotic excitations are investigated. Keywords: rotating beam, multi-bladed rotor, rotor mistuning, chaotic oscillations 1. Introduction Rotors are importantmachine componentswidely used in numerous industrial applications. The most common ones are helicopters, wind power turbines, fans, pumps, airplane propellers etc. Other advanced rotor design examples are rigid diskswith a series of beamelements combined in multi-stage assemblies. These are typically found in axial compressors, turbojet aircraft engines, steam and gas turbines etc. The considered systems are intensionally designed to be perfectly symmetric ones. However, multiple reasons may lead to symmetry break down in actual structures. The most common is rotor unbalance. This happens when mass centerline of the rotor does not coincide with the rotational axis. As a consequence, periodic inertia forces arise and large amplitude lateral vibra- tions of the structure can occur. The loss of rotor symmetrymight happen also due to operating wear of blades airfoil, possible crack development, randomvariations inmaterial properties and tolerances in manufacturing processes. Any of these factors leads to deviation in mechanical properties of the affected blade. This difference in blade to blade properties is referred to as rotor mistuning. The impact of systemmistuning on the free and forced response of multi-bladed rotors is of great scientific and practical interest. This stems from a few reasons. Primarily, the mistuning phenomenon has a fundamental influence on rotor dynamic characteristics. It turns out that mistuned multiple-bladed rotors can exhibit drastically larger forced response levels than the perfectly tuned designs (Xiao et al., 2004). For example, a 5% variation in the blade cantilever frequencies on a 92-bladed high pressure turbine disk can lead to one blade suffering a response magnitude of over 500% of the one observed on the perfectly tuned equivalent disk (Petrov and Ewins, 2003). This effect leads to large increases in stress and vibration amplitudes resulting in high cyclic fatigue and short lifetime of the rotor. An important observation regarding the increased vibration amplitudes is that they are not evenly distributed around all the rotor 550 J.Warminski et al. blades. The vibration energy is usually spatially concentrated making the modal motion to be limited to few blades only. This effect is known as mode localisation phenomenon and can be observed in any coupled periodic systems if a group of tuned structuremodes forms a complete set along the periodicity direction (Chen and Shen, 2015). The presence of irregularities in such structures restricts propagation of vibrations and confines the energy to a region close to the vibration source. However, as reported byVakakis et al. (1993), themistuning and the resulting localisation phenomena may originate also from structural nonlinearities and, thus, it may be observed even in perfectly periodic systems. Several concepts have been developed and discussed in the professional literature to reliably quantify mistuned rotor characteristics. These can be generally divided into either statistical or deterministic approach to the problem. Studying the dynamics of mistuned rotors within the frame of statistical analysis requires information about forced response probability density functions. These distributions are usual- ly determined within a series of Monte Carlo simulations. To this aim, a random population of test designs has to be given. This can be generated by assuming a priori the deviations in systemmass and stiffness matrices (Sinha, 1986). More representative designsmay be found by Taylor series expansions in terms of variables describing geometric variations of mistuned bla- des (Bhartiya and Sinha, 2013), from the experimental data tests (Li et al., 2006) or by means of a perturbation technique as proposed by Mignolet and Lin (1993). An alternative approach to approximate the mistuned rotor forced response probability density functions was proposed by Sinha (2006). The author tested the use of polynomial chaos for modal stiffness estima- tion and to compute the statistics of the forced response of a mistuned bladed disk assembly analytically. Although themistuning effect is a stochastic process (due to randomness of geometric and/or material perturbations) several deterministic approaches to the problem have been also deve- loped. The earliest studies considered multi-bladed rotors as a series of springs and lumped masses used to represent the blades and additional springs and dampers to model blade-to-hub and blade-to-blade interactions (Griffin andHoosac, 1984). Values of these systemvariables had to be determined through a parameter identification procedure. Despite its simplicity and ava- ilability of more advancedmodels these types of simplified formulations are still in use (Nikolic, 2006). Around the same time as analytical lumpedmasses models were being elaborated, the first finite element code software started showing up. This allowedmore strict investigation of actual rotor structures, taking into account complicated blades geometry, shrouds, aerodynamics and fluid-structure interactions etc. Over time, a reduced order treatment and sub-structuring tech- niques were developed to simplify a complete bladed disk finite element model to a smaller and more tractable problem. The most common ones used for mistuned rotors analysis are compo- nent mode synthesis (CMS), component mode based (CMB) method, subset of nominal modes (SNM)method andmodifiedmodal domain (MMD) analysis (Castanier and Pierre, 2006). For instance, the CMS method proposed by Castanier et al. (1997) is capable of reducing the ove- rall number of FE degrees of freedom up to three orders of magnitude if compared to the full structuremodel. Another deterministic approach to the analysis of mistuned rotors is based on a combination of the perturbation method and sensitivity analysis for natural frequencies and mode shapes (Shapiro, 1999). Alternatively, the analysis of rotors dynamics may be performed in the framework of an exact analytical formulation. The governing equations are usually de- rived by means of Hamilton’s principle. Next, these are solved analytically or numerically to test the system stability, individual blades motion synchronisation etc. (Crespo da Silva, 1998; Chandiramani et al., 2002; Sinha, 2013). Dynamics of rotating structures was also investigated by the authors of this publication. In particular, in paper (Warminski et al., 2014) a nonlinear system composed of two pendula Vibration of a mistuned three-bladed rotor... 551 attached to a hub and rotating in a horizontal plane was examined. The synchronisation phe- nomenon and transitions through resonances were analysed considering the influence of the hub inertia and order of nonliearity in the problem formulation. The existence of chaotic oscilla- tions of the system and paths leading to chaos were demonstrated as well. Next, Latalski et al. (2017) studied the dynamics of a hub-composite beam structure. Parametric studies regarding the laminate orientation angle and different regular driving torque scenarios were performed. The possibilities to control this kind of systems bymeans of the saturation control methodwere examined in a later research (Warminski and Latalski, 2017). In the present contribution, the former authors studies are extended to accommodate the three-bladed rotor case. In the performed analysis, it is assumed that the beams are made of a multilayered composite material with manufacturing tolerances of reinforcing fibres orienta- tions. The different magnitude of manufacturing inaccuracies in the individual blades results in their different stiffnesses followed by the rotor mistuning. The forced response of the system under regular and chaotic excitations is investigated as well as synchronisation of individual blades motions. 2. Statement of the problem Let us consider a rotor consisting of three slender and elastic composite beams clamped at the rigid hub of inertia Jh. The blades are fitted so that their flapwise bending plane coin- cides with the rotor plane. The system is driven by an external torque Text inducing rota- tion about a fixed frame vertical axis CZ0. The temporary angular position of the hub is denoted by an angle ψ(t) – see Fig. 1. The beams are made of an eighteen-layered laminate of an unidirectional graphite-epoxy pre-preg material. The adopted specific stacking sequence (0/−60/60/0/−60/603/−602/02/−60/02/602/−60) results in nominally full isotropic compo- site material behaviour (Vannucci andVerchery, 2002). However, in the performed analysis it is assumed that the rotor can bemistuned because of manufacturing inaccuracies in the laminate. These are represented by deviations of reinforcing fibres orientations from their nominal values. The discussed misalignments break the intended macroscopic material isotropy and make the laminate exhibit fully orthotropic behaviour. More detailed information regarding the way the rotor is mistuned and the considered mistuning magnitudes is given in the next Section and Table 3. Fig. 1. Model of the rotating hub with three elastic beams 552 J.Warminski et al. 2.1. Equations of motion The partial differential equations of motion of the considered structure have been derived according to the extended Hamilton’s principle δJ = t2 ∫ t1 (δT − δU + δWext) dt =0 (2.1) where J is the action, T is the kinetic energy, U is the potential energy and the work done by the external forces is given by the Wext term. The kinetic energy of the structure is defined as T = 1 2 Jhψ̇ 2(t)+ 1 2 3 ∑ i=1 ∫ Vi ρṘ T i Ṙi dVi (2.2) where the designation ρ refers to the average composite material density, Ṙi is the velocity vector of an representative infinitesimal element of volume Vi of the beam i. The total potential energy of the system U = ∑3 i=1Ui comes from elastic deformations of each beam. Posing the as- sumptions regarding a general shape of an open cross-section and its in-plane non-deformability, the energy Ui for an individual specimen is given by (Latalski et al., 2017) Ui = 1 2 ∫ Vi (σxxεxx+σxzγxz +σxsγxs) dVi (2.3) where σxx, σxn, σxs and εxx, γxz, γxs are stresses and strains in the axial direction and transver- se and lateral shear planes, respectively (referring to the individual blade). Although the posed mathematical model of the beam is limited to the linear case, the nonlinear axial strain εxx definition is adopted to represent properly the blade stiffening effect arising from system rota- tion ψ̇(t) (Latalski et al., 2017;Mayo et al., 2004). Therefore, the appropriate expressions are as follows εxx = ∂Dx ∂x + 1 2 [(∂Dx ∂x )2 + (∂Dy ∂x )2 + (∂Dz ∂x )2] γxz = ∂Dx ∂z + ∂Dz ∂x γxs = γxy +2zϕ ′ = ∂Dx ∂y + ∂Dy ∂x +2zϕ′ (2.4) where Dx, Dy and Dz are axial, lateral and transverse displacements of the cross-section re- presentative point written down in the local coordinates frame (x,y,z). The γxs strain includes an additional component coming from specimen torsion where ϕ is the profile twist angle. This is mandatory for the points located out of the beam mid-surface (Librescu and Song, 2006). Moreover, it is worth noting that the first nonlinear term present in εxx is skipped in further calculations due to the order of magnitude with respect to the other ones (Librescu and Song, 2006). Bearing inmind least action principle (2.1) and considering energy variations, after integra- tionwith respect to time, a set of general partial differential equations ofmotion of the structure is derived. By neglecting the deformations occurring out of the rotor plane, the problem is sim- plified.Therefore, in the final formulation, only the lead-lag plane displacement, transverse shear and profile twist for each individual beam are considered. These equations are supplemented by the additional one representing the rigid hub rotation dynamics. For the sake of brevity, the Vibration of a mistuned three-bladed rotor... 553 complex form of the partial differential equations and associated boundary conditions are omit- ted here. However, their full formulation as well as the detailed step-by-step derivation can be found in the previous authors paper (Latalski et al., 2017). The derived system of partial differential governing equations is transformed into ordinary differential ones taking into account the normalmodes projection and the associated orthogona- lity condition. To this aim, the Galerkin procedure for the first natural mode is applied. Next, the system is converted into the dimensionless notation. The coupled flexural-torsional mode projection results in the final set of nonlinear ODEs as follows ( Jh+ 3 ∑ i=1 Jbi+ 3 ∑ i=1 αhi2q 2 i ) ψ̈+ ζhψ̇+ 3 ∑ i=1 (αhi1q̈i+αhi3qiq̇iψ̇)= µ(τ) q̈1+ ζ1q̇1+α12ψ̈+(α11+α13ψ̇ 2)q1+α14q1q̇1ψ̇ =0 q̈2+ ζ2q̇2+α22ψ̈+(α21+α23ψ̇ 2)q2+α24q2q̇2ψ̇ =0 q̈3+ ζ3q̇3+α32ψ̈+(α31+α33ψ̇ 2)q3+α34q3q̇3ψ̇ =0 (2.5) where Jh, Jbi denote the mass moment of inertia of the hub and each subsequent beam, re- spectively. These are relative values calculated with respect to beam 1. The factors ζh and ζi are hub and beams viscous damping coefficients. They have been estimated during laboratory experiments. The approximate value for every ζi is 0.04 ratio of its corresponding beam natural frequency ω0i, thus ζi = 0.04 √ ai1. The hub damping ζh has been set at 0.1. The external di- mensionless torque imposed to the hub is denoted by µ(τ) where τ is dimensionless time. The parameters αhij (j = 1,2,3) present in (2.5)1 and αik (k = 1, . . . ,4) in (2.5)2-4 are coefficients obtained from themodal reduction procedure. Studying the system of governing equations (2.5) onemay observe that the individual beams equations are coupled by inertia terms present in the hub equation. If angular velocity of the structure is constant (ψ̈ = 0) then all equations become uncoupled. This happens despite the quadratic terms present in the first equation since these terms are of a higher order and can be neglected for a small oscillations case. By contrast, if the angular velocity is not constant (ψ̈ 6= 0) then all equations are mutually coupled, and the nonlinear quadratic terms as well as Coriolis forces are involved in the full structure dynamics. 2.2. Accounting for reinforcing fibres misalignment To take into consideration the composite fibres orientation tolerances let us assume that the fiber angle in each k-th ply (k = 1, . . . ,N) may be deviated from its nominal value αk by the maximum acceptable tolerance limit ∆αk. Hence the actual orientation angle stays within a range 〈αk − ∆αk;αk + ∆αk〉 – see Fig. 2a. Moreover, one assumes that the magnitude of this misalignment is set arbitrary and it does not depend on the nominal fiber orientations αk. Therefore, the tolerance stays equal for all layers and, thus, ∆αk = ∆α. Fig. 2. Modelling the tolerances of laminate reinforcing fibres orientations 554 J.Warminski et al. This approach to model the imperfect laminate where only the nominal values and peak deviations of fiber orientation angles are given renders their exact values unknown. However, the maximum possible impact of the accepted ∆α inaccuracies on the mechanical properties of thematerialmay be estimated bymeans of a sensitivity analysis. To this aim, let us consider the requested beam stiffnesses occurring in the partial differential equations of motion of the hub- composite beamrotor as formulated in (Latalski et al., 2017) Eqs. (24)-(28). These stiffnesses are calculated according to the Librescu composite beam theory (Georgiades et al., 2014; Librescu and Song, 2006) and expressed as a33 = ∫ c [( A22− A12A12 A11 ) z2−2 ( B22− A12B12 A11 )dy ds z + ( D22− B12B12 A11 )(dy ds )2] ds a37 =2 ∫ c [( B26− A12B16 A11 ) z− ( D26− B12B16 A11 )dy ds ] ds a55 = ∫ c [( A66− A16A16 A11 )(dz ds )2 +A44 (dy ds )2] ds a77 =2 ∫ c ( D66− B16B16 A11 ) ds (2.6) where the dummy variable s represents the beam profile coordinate measured along its width. Terms a33, a55, and a77 are flapwise bending, transverse shear and twist stiffnesses, respectively. The parameter a37 represents the laminate stiffness in coupled flapwise bending-twist deforma- tion. In the nominal design, this is equal to zero since the considered multilayered laminate is macroscopically isotropic. However, for a material with misaligned layers it is expected to be different from zero. This is confirmed by results given in Table 1. The above given stiffnesses are expressed in terms of individual elements ofA,D andB ten- sors representing stretching, bending and bending-stretching stiffnesses, respectively. Following the Classical Laminate Theory they are given as Aij = 18 ∑ k=1 Q (k) ij (zk −zk−1) Bij = 1 2 18 ∑ k=1 Q (k) ij (z 2 k −z2k−1) Dij = 1 3 18 ∑ k=1 Q (k) ij (z 3 k −z3k−1) (2.7) where zk, zk−1 are the distances from the reference plane to the two surfaces of the k-th ply (see Fig. 2b) and Q (k) ij are themembers of the reduced stiffnessmatrix of this ply. These recent ones depend on the fibres orientation angle αk Q (k) 11 = c 4Q11+s 4Q22+2c 2s2(Q12+2Q66) Q (k) 22 = s 4Q11+ c 4Q22+2c 2s2(Q12+2Q66) Q (k) 12 = c 2s2(Q11+Q22−4Q66)+(c4+s4)Q12 Q (k) 66 = c 2s2(Q11+Q22−2Q12)+(c2−s2)2Q66 Q (k) 16 = cs [ c2Q11−s2Q22− (c2−s2)(Q12+2Q66) ] Q (k) 26 = cs [ s2Q11− c2Q22+(c2−s2)(Q12+2Q66) ] (2.8) where c =cosαk and s =sinαk. Vibration of a mistuned three-bladed rotor... 555 Bearing the above relations in mind, the impact of the considered reinforcing fibers misa- lignment ∆αk on the beam mechanical properties given by Eqs. (2.6) may be evaluated by performing a sensitivity analysis. The change of any beam stiffness coefficient aij is ∆aij = 18 ∑ k=1 ∣ ∣ ∣ ∂aij ∂αk ∆αk ∣ ∣ ∣ (2.9) where ij pair is 33, 37, 55 or 77. The term ∂aij/∂αk represents the sensitivity of the specimen stiffness with respect to changes in fibre orientations in the k-th individual laminate layer. Since the signs of the derivative as well as the accepted inaccuracy ∆αk are arbitrary, an absolute value operator is used to consider the most unfavourable case. Finally, two limit values for any perturbed stiffness aij are possible, namely the lower one aij − ∆aij for a decreased stiffness and the upper one aij +∆aij for an increased value. The presented above treatment based on perturbation calculus is often encountered in reliability based structural design and represents the so called ‘a worst case scenario’ analysis (Gutkowski and Latalski, 2003). Calculations of ∂aij/∂αk sensitivities involve differentiation of the individual elements of stretching, bendingandbending-stretching tensors (2.7)with respect to thefibres angleα. These derivatives are easily found by substituting toEqs. (2.7) the expressions given by relations (2.8). The results of stiffnesses (2.6) calculations for the proposed 18 layers stacking sequence laminate beam are collected in Table 1. In the first section, the values for the nominal design are given. Next, the small misalignment ∆α = 1◦ is assumed and the perturbed values for increased and decreased stiffnesses are printed. Finally, the ∆α = 5◦ case is considered. The graphite-epoxy material data used for these calculations are given in Table 2. Table 1. Beam stiffnesses for the partial differential equations of motion; nominal design and two perturbed cases by ∆α =1◦ and ∆α =5◦ Nominal design values a33 =0.117568Nm 2 a37 =0.0Nm 2 a55 =67957.50N a77 =0.081397Nm 2 Fibres misalignment ∆α =1◦ perturbed values for increased stiffness a33 =0.121849Nm 2 a37 =0.000081Nm 2 a55 =68273.55N a77 =0.083996Nm 2 perturbed values for decreased stiffness a33 =0.113266Nm 2 a37 =−0.000138Nm2 a55 =67641.44N a77 =0.078523Nm2 Fibres misalignment ∆α =5◦ perturbed values for increased stiffness a33 =0.138754Nm 2 a37 =−0.000116Nm2 a55 =69537.77N a77 =0.091819Nm2 perturbed values for decreased stiffness a33 =0.095806Nm 2 a37 =−0.001301Nm2 a55 =66377.23N a77 =0.064011Nm2 Table 2.Rotor geometric data andmaterial properties used in numerical simulations Geometric properties l =0.350m d =0.034m h =0.0009m R0 =0.1× l Material properties of the laminate E1 =143.2GPa E2 =E3 =3.1GPa G23 =2.05GPa G12 = G13 =3.28GPa ν21 = ν31 =0.00758[–] ν32 =0.2439[–] ρc =1350.0kg·m3 556 J.Warminski et al. The performed numerical calculations confirm the already reported fully isotropic proper- ties of the assumed nominal design stacking sequence laminate (a37 = 0). However, as can be observed, the discussed stacking sequence configuration is sensitive to possible variations in fi- bres orientations as confirmed bymeaningful changes in the beam stiffnesses. In particular, the possible very small misalignment of the fibres angle leads to anisotropic material behaviour and induces mutual coupling of different components of specimen deformations (a37 6= 0). Finally, the values of the coefficients present in the ordinary differential equations of motion (system of Eqs. (2.5)) for the assumed mistuning magnitudes are listed in Table 3. The details regarding the three studied cases of mistuned rotor configurations are also given there. Table 3. The values of individual coefficients present in the ordinary differential equations of motion; three studied rotor configurations Case 1.Rotor with nominal design blades Beam 1 α11 =12.364453698 α12 =1.779913785 α13 =0.350955874 α14 =−1.551795958 Beam 2 α21 =12.364453698 α22 =1.779913785 α23 =0.350955874 α24 =−1.551795958 Beam 3 α31 =12.364453698 α32 =1.779913785 α33 =0.350955874 α34 =−1.551795958 Hub αh11 =−0.530660819 αh12 =−0.402771952 αh13 =−0.805543905 αh21 =−0.530660819 αh22 =−0.402771952 αh23 =−0.805543905 αh31 =−0.530660819 αh32 =−0.402771952 αh33 =−0.805543905 Case 2.Tolerance ∆α =1◦: beam 1 – nominal, beam 2 – decreased stiffness, beam 3 – increased stiffness Beam 1 α11 =12.364453698 α12 =−1.779913785 α13 =0.350955874 α14 =1.551795958 Beam 2 α21 =11.911908455 α22 =−1.779874982 α23 =0.350960254 α24 =1.551834362 Beam 3 α31 =12.814170650 α32 =−1.779950210 α33 =0.350951758 α34 =1.551759829 Hub αh11 =0.530660819 αh12 =−0.402771952 αh13 =−0.805543905 αh21 =0.530672090 αh22 =−0.402790584 αh23 =−0.805581169 αh31 =0.530650234 αh32 =−0.402754446 αh33 =−0.805508893 Case 3.Tolerance ∆α =5◦: beam 1 – nominal, beam 2 – decreased stiffness, beam 3 – increased stiffness Beam 1 α11 =12.364453698 α12 =−1.779913785 α13 =0.350955874 α14 =1.551795958 Beam 2 α21 =10.072877023 α22 =−1.779688026 α23 =0.350983037 α24 =1.552014103 Beam 3 α31 =14.591731254 α32 =−1.780075869 α33 =0.350937164 α34 =1.551634747 Hub αh11 =0.530660819 αh12 =−0.402771952 αh13 =−0.805543905 αh21 =0.530725703 αh22 =−0.402879065 αh23 =−0.805758131 αh31 =0.530613777 αh32 =−0.402693994 αh33 =−0.805387989 As the last remark to this Section, it should benoted that density of thematerial is not affec- ted by the deviations of composite reinforcing fibers from their nominal orientations. Therefore, the inertia coefficients for the nominal andmisaligned composite beams are similar. 3. Numerical results 3.1. Forced vibrations – regular oscillations To get inside into the dynamics of the structure the full nonlinear system of governing equations (2.5) is solved numerically and appropriate diagrams are prepared. To this end, the Auto-07p software (Doedel et al., 1998) adopting the multiparameter continuation method has been used. To start the analysis, let us consider the system response if the structure is excited by an external torque imposed to the hub. In general, we consider the forcing to be composed of two terms, namely the constant and the harmonic one. Thus, the driving torque is defined as µ(τ)= µ0+ρsin(ωτ) (3.1) Vibration of a mistuned three-bladed rotor... 557 where µ0 is a constant component, while ρ and ω correspond to the amplitude and frequency of the periodic term, respectively. For the purpose of an initial analysis, we assume the constant component to be equals to zero, µ0 = 0. The aim of this analysis is to establish the inherent dynamicproperties of thediscussed systemwithout the additional stiffening effect resulting from the full rotational motion. At the first stage, dynamics of the reference structure is examined (nominal design, case 1 in Table 3). In Fig. 3, we present resonance curves for all three beams as well as the angular velocity of the hub (i.e. maximummagnitudes of these four variables) if the constant component of the torque is neglected (µ0 =0). One can observe, for the fully symmetric rotor, just a single resonance zone is present with the peak response corresponding to the first natural frequency ω01 =3.0217. In this resonance zone, motions of all the blades are fully synchronised in magni- tude and phase and they stay in anti-phase with respect to the hub motion shown in Fig. 3b. It is worth to note that for very small values of excitation frequency ω angular speed of the system gets large values – Fig. 3b. This case corresponds to the zero natural frequency and the excitation of the systemwith the infinite period. However, this range of excitation is out of our interest in this paper, therefore in next figures we perform analysis around non zero resonance zones only. Fig. 3. Resonance curves of beams q1, q2, q3 (a) and angular velocity of the hub ψ̇ (b); variant 1 of the three reference (nominal) beams, µ0 =0, ρ =0.01, Jh =5 Fig. 4. Resonance curves for individual beams q1 – black, q2 – blue, q3 – green (a) and angular velocity of the hub ψ̇ (b). Rotor configuration in variant 2, ρ =0.01, Jh =5 558 J.Warminski et al. In the case of marginally mistuned beams (∆α = 1◦) represented by studied variant 2 the response amplitudes of individual blades are slightly different and are dominated by the most flexible beam 2 – presented in Fig. 4a. Similar to the above nominal design case, all the blades aremutually synchronised and oscillate together in anti-phasewith respect to the hub– see time histories plots shown in Fig. 5a. Apart from the main resonance, two other minor resonances occur due to the system mistuning. They are shifted to the right with respect to the main resonance peak – Fig. 4a. Time histories for those additional resonances are presented in plots shown in Figs. 5b and 5c. It is to be commented that the phase shift between individual beams and the hub is a constant value. Moreover, note that the response of the hub is not affected in any way by these two minor resonances – Fig. 4b. Fig. 5. Time histories of beams q1 – black, q2 – blue, q3 – green and the angle of hub rotation ψ –magenta for variant 2 and ω =3.01 (a), ω =3.50 (b), ω =3.58 (c); µ0 =0, Jh =5, ρ =0.01 The effect of bladesmistuning ismuchmore evident in variant 3 (see Table 3) corresponding to the fibre orientation tolerance limit ∆α = 5◦. The resonance curves and time histories are shown in Fig. 6 and Fig. 7, respectively. Fig. 6. Resonance curves for individual rotor beams q1 – black, q2 – blue, q3 – green (a) and angular velocity of the hub ψ̇ (b). Rotor configuration in variant 3, ρ =0.01, Jh =5 The direct comparison of system responses for the small and relatively large mistuning – Figs. 4 and 6 – reveals some qualitative differences. Primarily, in the case of a higher structural mistuning (variant 3) all three resonances are distinct and very well observed. This conclusion refers also to the response of the hub which increases in all three resonance zones – Fig. 6b – and is in contrast to the small mistuning case behaviour shown in Fig. 4b. Furthermore, a kind Vibration of a mistuned three-bladed rotor... 559 Fig. 7. Time histories of individual beams q1 – black, q2 – blue, q3 – green and the angle of hub rotation ψ –magenta for variant 3 and the frequency of excitation close to the subsequent resonance zones for ω =2.9 (a), ω =3.4 (b), ω =3.75 (c); Jh =5, ρ =0.01 of vibration localisation can be noticed. For the first resonance zone the motion is localised in the second beam (q2 coordinate in Fig. 6a), for the second resonance the motion is localised in beam number one (q1 coordinate), and finally for the third resonance the motion is localised in the third beam (q3 coordinate). This is confirmed by the time history plots and individual blades amplitudes as shown in Figs. 7a-7c. To examine the influence of rotational motion and anticipated stiffening effects, the relati- vely large mistuning (∆α = 5◦) variant is studied again and the bifurcation diagram for the µ0 component of the torque is presented in Fig. 8. It is to be noted here that any non-zeromean value of the driving torque represents the case when the system is accelerating from the zero initial velocity but only up to a certain moment where this torque is balanced by damping on the hub. Finally, the system is performing full rotation with a constant non-zero mean value angular velocity. Fig. 8. Bifurcation diagram for the constant torque component µ0. Rotor configuration in variant 3, ω =2.8 (a), ω =3.0 (b); ρ =0.01, Jh =5 To study the dynamics,we selected two frequencies around the first resonance peak as shown in Fig. 6, namely ω =2.8 and ω =3.0. These correspond to the situations before and after the resonance, respectively. As expected, full rotation of the structure changes the amplitudes of blades vibrations due to the centrifugal stiffening effect. However, comparing these two plots, one observes two different possible courses of the beams responsewhile the µ0 is varied.Thefirst 560 J.Warminski et al. one (Fig. 8a corresponding to the excitation frequency ω = 2.8) is the monotonic reduction in amplitudes.This occursparticularyup toµ0 ≈ 1andnextbecomesmuch less pronounced.Then the beam vibrations are strongly suppressed, the differences in blades oscillations are negligible and they behave almost like rigid bodies. The different scenario is observed for the case ω =3.0 (Fig. 8b).While changing the µ0 component, the amplitudes initially increase and then decrease rapidly to small values. The peak inFig. 8b is a direct result of the angular speedwhich shifts to the right the natural frequency of the blades. Therefore, the excitation ω =3.0 occurs right at the new resonance zone. Nevertheless, the further increase in µ0 results in a rapid reduction in beam amplitudes, and for µ0 > 1 they are suppressed similarly to the case presented in Fig. 8a. To study the influence of the torque amplitude ρ, the last analysis is repeated by setting µ0 =0.13 which is similar to the average torque value as assumed in computations for a chaotic signal presented later in Section 3.2. Two torque oscillations amplitudes are analysed, namely ρ = 0.27 and ρ = 1.0 as shown in Figs. 9 and 10, respectively. Studying these curves, one Fig. 9. Resonance curves for individual rotor beams q1, q2, q3 (a) and angular velocity of the hub ψ̇ (b). Rotor configuration in variant 3; µ0 =0.13, ρ =0.27, Jh =5 Fig. 10. Resonance curves for individual rotor beams q1, q2, q3 (a) and angular velocity of the hub ψ̇ (b). Rotor configuration in variant 3; µ0 =0.13, ρ =1.0, Jh =5 observes again the centrifugal stiffening effect that shifts the resonance peeks towards higher frequencies – check Fig. 6 for reference. Furthermore, direct comparison of Figs. 9 and 10 reveals the linear nature of the response curves for relatively small oscillations. By contrast, if the Vibration of a mistuned three-bladed rotor... 561 amplitude of excitation is increased then the resonance curves presented in Fig. 10 demonstrate a minor nonlinear softening phenomenon. However, these results are out of real dynamics of the structure. Therefore, it can be concluded that the higher order terms present in governing equations (2.5) do not have a significant impact on structural dynamic properties. When studying the system of governing equations (2.5), wemay expect that the hubmotion plays essential role in the entire rotor dynamics. To get more insight into this relation, let us consider the solution corresponding to the peak in the first resonance zone occurring at ω =2.9, as shown inFig. 6, andvary themassmomentof inertiaJh of thehub.Theperformedsimulations show the response peak to be reduced if the mass moment of inertia Jh = 5 is changed either towards lower or higher values. Theoutcomes of these tests are shown in the bifurcation diagram presented in Fig. 11. Fig. 11. Influence of hub inertia on the beams responses q1 – black, q2 – blue, q3 – green (a) and angular velocity of the hub ψ̇ (b). Calculations for the data corresponding to variant 3; ρ =0.1, ω =2.9 The next series of simulations demonstrates the influence of hub inertia if the excitation frequency is varied. The tests are performed for different hubs starting fromvery light ones with Jh =0.1 up to the heaviest with Jh =50. The corresponding curve plots are shown in Fig. 12. Evidently, larger hub inertias lead to vibration reduction. Moreover, we may also observe an interesting phenomenon of vibrations localisation clearly visible around frequencies ω =3.1 and ω =3.5, which is independent of the hub inertia (Fig. 12). 3.2. Chaotic oscillations Let us consider the casewhen the discussed three-bladed rotor is excited by a chaotic Duffing type oscillator. Thedimensionless driving torqueµ(τ) inEq. (2.5)1 is considered now to be given by the formula µ = αxx (αx =1), where the variable x is calculated fromDuffing’s equation ẍ+kẋ+x3 = β +ρcos(ωτ) (3.2) For numerical simulations, the constants k = 0.05, β = 0.03, ω = 1 and ρ = 0.16 are adopted as originally used by Ueda (1980). This results in the chaotic torque characteristic with the mean value of about 0.13, which is the one that has been used for harmonic excitation results presented in Figs. 9 and 10. The presented below Poincaré maps have been obtained by means of direct numerical simulations performed in Dynamics software. 562 J.Warminski et al. Fig. 12. Resonance curves of beams (a) q1, (b) q2, (c) q3, and (d) angular velocity of the hub ψ̇. Calculations for the data corresponding to variant 3 and various hub inertias: Jh =0.1 – black, Jh =0.5 – blue, Jh =5 – green, Jh =10 – red, Jh =50 –magenta; ρ =0.1 Fig. 13. Poincarémaps of beams responses (a) q1, (b) q2, (c) q3 subject to chaotic excitation given by Eq. (3.2). Rotor configuration in variant 3; Jh =5 The strange chaotic attractors are found for all the beams, both for slightly mistuned va- riant 2 and for highly mistuned configuration given by case 3. The results for this last one are shown in Fig. 13. Since the results for case 2 do not show evident differences between individual blades they are not presented here. Vibration of a mistuned three-bladed rotor... 563 As can be observed in Fig. 13, the strange attractor for the second beam is significantly bigger than the attractor for nominal beam 1 – subplots (a) and (b), respectively. Meanwhile, the strange attractor of beam 3 is smaller – Fig. 13c. This comes from the significant differences in stiffness of individual beams as already reported while discussing the regular excitation case. To study the possible synchronisation phenomenon, the time series plots are prepared. The time histories for a highly mistuned configuration given by case 3 are shown in Fig. 14. The motions of all the beams are well synchronised, although the amplitudes are slightly different. The amplitude of the second beam is the highest, while the amplitude of the third beam is the smallest one.Again, these disparities can be explainedby thedifferences in specimens stiffnesses. Fig. 14. Time histories of the beams responses, the full time domain (a) and zoomwindow (b); q1 – black, q2 – blue, q3 – green. Rotor configuration in variant 3, chaotic excitation; Jh =5 The final stage of numerical simulations comprises the studies of system dynamics if both the hub and the first beam are subjected to the same chaotic excitation signal µ = αxx, as given in the previous example. To this end, the additional forcing term µ is added to the right hand side of the governing equation of beam 1 – Eq. (2.5)2. The Poincaré maps are plotted for each beam individually as well as for the oscillator itself. Defined variant 3 of high structural mistuning is examined and corresponding plots are gathered in Fig. 15. It can be observed that the strange attractor for the nominal beam is the biggest one. This results from the excitation that has been applied to the hub and directly to this beam as well. Meanwhile, the Poincaré map for the most rigid beam (blade 3) is the smallest – Fig. 15c. In contrast to the previously studied case of chaotic excitation applied only to the hub now, the synchronisation scenario is different. Beam 1 oscillates in anti-phase with respect to complectly synchronised beams 2 and 3. 4. Final remarks and conclusions The presented paper considers forced vibrations of a rotating structure consisting of a rigid hub and three flexible beams made of a multilayered laminate. In the performed analysis, it is assumed that the rotor has beenmistuned due tomanufacturing tolerances of reinforcing fibres orientations in the composite material. Based on previous authors research, a system of four mutually coupled dimensionless ordinary differential governing equations has been formulated. The forced vibrations of the systemunder a regular excitation supplied to the hubhave been investigated. Thedirect comparison of system responses for small (∆α =1◦) and relatively large 564 J.Warminski et al. Fig. 15. Poincaré maps of the beams response (a) q1, (b) q2, (c) q3, (d) x chaotic excitation given by Eq. (3.2) and applied to the hub and nominal beam 1. (e) Time histories of individual beams responses q1 – black, q2 – blue, q3 – green. Rotor configuration in variant 3; Jh =5 (∆α = 5◦) mistuning cases has revealed some qualitative differences. This has concerned not only oscillation amplitudes but synchronisation of motions as well. Further studies have dealt with the influence of hub inertia on themistuned system dynamics. An interesting phenomenon of vibrations localisation clearly visible around two distinct frequencies has been observed. This effect has been shown to be independent of the hub inertia magnitude. Finally, the forced response of the structure under a chaotic excitation has been investigated. Two cases of forcing load have been examined, namely (a) forcing supplied to the hub only and (b) to the hub and to one of the rotor beams. The obtained results for both cases have revealed some disparities in individual beams strange attractors magnitudes due to differences in blades stiffnesses. Analysis of motion of the system with the excitation applied to the hub has showed full synchronisation of the three beams motions despite chaotic motion of the full structure. However, the chaotic excitation applied to the beam and to the hub has resulted in the oscillations of beam 1 in anti-phase with respect to complectly synchronised beams 2 and 3. Acknowledgments Thework is financially supported by the grantDEC-2015/19/N/ST8/03906 from thePolishNational Science Centre. References 1. BhartiyaY., SinhaA., 2013,Reducedordermodelingof abladed rotorwithgeometricmistuning via estimated deviations in mass and stiffness matrices, Journal of Engineering for Gas Turbines and Power, 135, 5, 052501, DOI: 10.1115/1.4007783 Vibration of a mistuned three-bladed rotor... 565 2. Castanier M.P., Óttarsson G., Pierre C., 1997, A reduced order modeling technique for mistuned bladed disks, Journal of Vibration and Acoustics, 119, 3, 439, DOI: 10.1115/1.2889743 3. Castanier M.P., Pierre C., 2006, Modeling and analysis of mistuned bladed disk vibration: current status and emerging directions, Journal of Propulsion and Power, 22, 2, 384-396, DOI: 10.2514/1.16345 4. Chandiramani N.K., Librescu L., SheteC.D., 2002,On the free-vibration of rotating compo- site beams using a higher-order shear formulation,Aerospace Science andTechnology, 6, 8, 545-561 5. Chen Y.F., Shen I.Y., 2015, Mathematical insights into linear mode localization in nearly cyc- lic symmetric rotors with mistune, Journal of Vibration and Acoustics, 137, 4, 041007, DOI: 10.1115/1.4029945 6. Crespo da Silva M.R.M., 1998, A comprehensive analysis of the dynamics of a helicopter ro- tor blade, International Journal of Solids and Structures, 35, 7-8, 619-635, DOI: 10.1016/S0020- 7683(97)00065-6 7. Doedel E., ChampneysA., FairgrieveT., KuznetsovY., SandstedeB.,WangX., 1998, Auto 97: Continuation and Bifurcation Software for Ordinary Differential Equations, Computatio- nal Mathematics and Visualization Laboratory, http://indy.cs.concordia.ca/aut 8. Georgiades F., Latalski J., Warminski J., 2014, Equations of motion of rotating compo- site beams with a nonconstant rotation speed and an arbitrary preset angle, Meccanica, 49, 8, 1833-1858, DOI: 10.1007/s11012-014-9926-9 9. Griffin J.H., Hoosac T.M., 1984, Model development and statistical investigation of turbine blademistuning, Journal of Vibration Acoustics Stress and Reliability in Design,106, 2, 204,DOI: 10.1115/1.3269170 10. GutkowskiW., Latalski J., 2003,Manufacturing tolerances of fiber orientations inoptimization of laminatedplates,EngineeringOptimization,35, 2, 201-213,DOI: 10.1080/0305215031000091550 11. Latalski J., Warminski J., Rega G., 2017, Bending-twisting vibrations of a rotating hub- thin-walled composite beam system,Mathematics andMechanics of Solids, 22, 6, 1303-1325,DOI: 10.1177/1081286516629768 12. LiJ.,CastanierM.P.,PierreC.,CeccioS., 2006,ExperimentalMonteCarlomistuningasses- sment of bladed disk vibration using forcing variations, [In:] 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, pp. AIAA 2006-1964 13. Librescu L., Song O., 2006,Thin-Walled Composite Beams: Theory and Application, Springer, Dordrecht and the Netherlands 14. Mayo J., Garćıa-Vallejo D., Doḿınguez J., 2004, Study of the geometric stiffening ef- fect: comparison of different formulations, Multibody System Dynamics, 11, 4, 321-341, DOI: 10.1023/B:MUBO.0000040799.63053.d9 15. Mignolet M.P., Lin C.-C., 1993, The combined closed form-perturbation approach to the ana- lysis of mistuned bladed disks, Journal of Turbomachinery, 115, 4, 771, DOI: 10.1115/1.2929315 16. Nikolic M., 2006,New Insights into the Blade Mistuning Problem, PhDThesis, Imperial College London, University of London, London 17. Petrov E.P., Ewins D.J., 2003, Analysis of the worstmistuning patterns in bladed disk assem- blies, Journal of Turbomachinery, 125, 4, 623-631, DOI: 10.1115/1.1622710 18. Shapiro B., 1999, Solving for mistuned forced response by symmetry, Journal of Propulsion and Power, 15, 2, 310-325, DOI: 10.2514/2.5429 19. Sinha A., 1986, Calculating the statistics of forced response of a mistuned bladed disk assembly, AIAA Journal, 24, 11, 1797-1801,DOI: 10.2514/3.9526 20. SinhaA., 2006,Computationof the statisticsof forcedresponseof amistunedbladeddiskassembly via polynomial chaos, Journal of Vibration and Acoustics, 128, 4, 449, DOI: 10.1115/1.2215620 566 J.Warminski et al. 21. Sinha S.K., 2013, Rotordynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contact-impact rub loads, Journal of Sound and Vibration, 332, 9, 2253-2283, DOI: 10.1016/j.jsv.2012.11.033 22. UedaY., 1980,Explosion of strange attractors exhibited byDuffing’s equation,Annals of the New York Academy of Sciences, 357, 1, 422-434, DOI: 10.1111/j.1749-6632.1980.tb29708.x 23. VakakisA., NayfehT., KingM., 1993,Amultiple-scales analysis of nonlinear, localizedmodes in a cyclic periodic system, Journal of Applied Mechanics, 60, 2, 388, DOI: 10.1115/1.2900806 24. Vannucci P., Verchery G., 2002, A newmethod for generating fully isotropic laminates,Com- posite Structures, 58, 1, 75-82, DOI: 10.1016/S0263-8223(02)00038-7 25. Warminski J., Latalski J., 2017, Nonlinear control of flexural-torsional vibrations of a rotating thin-walled composite beam, International Journal of Structural Stability and Dynamics, 17, 5, 1740003-1–1740003-17,DOI: 10.1142/S021945541740003X 26. Warmiński J., Szmit Z., Latalski J., 2014, Nonlinear dynamics and synchronisation of pendula attached to a rotating hub, European Physical Journal Special Topics, 223, 827847, DOI:10.1140/epjst/e2014-02143-9 27. Xiao B., Rivas-Guerra A., Mignolet M.P., 2004,Maximum amplification of blade response due to mistuning in multi-degree-of-freedom blade models,Proceedings of the ASME Turbo Expo: Power for Land, Sea, and Air, 6, 427-438 Manuscript received January 24, 2018; accepted for print March 14, 2018