Nonlinear vibrations of a rotating thin-walled composite piezo-beam with circumferentially uniform stiffness (CUS)

A rotating system comprising a hub and a thin-walled laminate cantilever beam with embedded nonlinear piezoelectric layers is analysed in the paper. The reinforcing fibres set-up in composite material conforms to circumferentially uniform stiffness lamination scheme. This configuration exhibits the mutual bending couplings in two orthogonal planes. Nonlinear analytical model of a piezoelectric material embedded onto the beam walls is postulated by considering the higher-order constitutive relations with respect to electric field variable. Moreover, to properly model electromechanical structural behaviour, the full two-way coupling piezoelectric effect is considered. To this aim, the assumption of a spanwise electric field variation is postulated in the mathematical model of the structure. Based on previous authors’ research, the system of mutually coupled nonlinear equations of motion is formulated. In the numerical analysis the forced response of the system under zero and nonzero mean value harmonic torque excitation is considered. In particular, the influence of hub inertia, excitation amplitude and mean rotating speed on system dynamics is investigated. The results are presented in the form of appropriate frequency response plots and bifurcation diagrams.

rotating speed on system dynamics is investigated. The results are presented in the form of appropriate frequency response plots and bifurcation diagrams.

List of Symbolŝ b
Effective electrostrictive constants tensor D Electric displacement vector E Electric field vector e Tensor of piezoelectric coefficients J h Dimensionless inertia of the hub (calculated with respect to the beam inertia) a 22 Chordwise bending stiffness a 25 Coupled chordwise bending-flapwise transverse shear stiffness a 33 Flapwise bending stiffness a 34 Coupled flapwise bending-chordwise transverse shear stiffness a 44 Chordwise transverse shear stiffness a 55 Flapwise transverse shear stiffness b 1 Translational inertia of the cross section B 4 Cross-section moment of inertia about y axis B 5 Cross-section moment of inertia about z axis c Height of the cross section d Width of the cross section J h Hub inertia R 0 Hub radius T ext,z External driving torque applied to the hub v 0 Displacement of the cross-section origin O along the local y axis w 0 Displacement of the cross-section origin O along the local z axis X 0 Y 0 Z 0 Global (inertial) coordinate system x yz Local coordinate system attributed to beam cross-section α Reinforcing fibres orientation angle χ Third-order electric susceptibility tensor ξ Second-order permittivity tensor η Dimensionless abscissa along the beam span (η ∈ 0, 1 ) μ Dimensionless driving torque supplied to the hub μ 0 Constant component of the dimensionless driving torque ψ Angular position of the rotor in the global coordinate system X 0 Y 0 Z 0 ρ Dimensionless amplitude of the driving torque θ Beam presetting angle ϑ y Rotation of the beam cross section around the local y axis ϑ z Rotation of the beam cross section around the local z axis

Introduction
The continuous development of composite materials offers a great potential for modern structural systems that can take advantage of unique composite material properties as well as material tailoring technology. Historically, this later idea is exploited since late 1980s and relies on proper laminae stacking sequence and orientation of unidirectional (UD) reinforcing fibres. The main objective of this technology is to reduce weight of the structure, improve its stiffness and maintain static strength under compressive loads [10]. Other applications of the material tailoring technology in aerospace industry include the drag reduction [43], gust response [13] and optimum aeroelastic behaviour and flutter characteristics [1]. One of the most prominent demonstrators of this technology was the test programme of the forward swept wing experimental fighter Grummann X-29 [3]. Further progress within the material tailoring technology is the tow-steering technique, where the fibres orientations within a ply follow a predefined curvilinear path; thus the ply stiff-ness vary continuously through the plane of each ply. It has been demonstrated that these fibre paths can be optimized to increase the structural performance of a laminate beyond that of an equivalent regular UD laminate [26,32,38]. Further improvement in system performance and enhancement in its capabilities can be achieved by incorporating multifunctional materials into structural design. This technology offers, along with other concurrent abilities, advantages of active control, sensing, self-healing and thermal functionality of the structure. Typical applications of multifunctional designs range from civil structures, marine systems and automobiles, through machine tools and fine mechanics to aerospace and aeronautics systems and structures. Representative examples of active aerospace designs might be inflatable and deployable space structures such as solar panels or space antennas and mirrors that often undergo large controlled rotations and positioning when expanding from an initially compact configuration to a final geometry [25]. Referring to other aeronautical systems like fixed-wing aircraft, applications cover active control of flutter and handling qualities improvement. In particular, the potential of using piezoceramic transducers to suppress or delay the onset of panel aeroelastic instability has been studied. The conducted research have shown the possibility to increase the flutter speed up to 50% combined with reduction of the power-spectral density of response. Moreover, the potential savings in control effort resulting from optimal positioning of active elements on the wings and ailerons have been confirmed [34]. Further applications areas of piezoceramic materials cover structural health monitoring and interior/exterior noise reduction. The technology of smart material-based SHM sensors enhances the system flight safety and reliability, but also results in savings in operational and maintenance/inspection costs and extends the life cycle of an ageing aircraft [44]. For rotarywing aircrafts the embedded smart material sensors can be used for rotor tracking and head health monitoring, while active material actuators may induce changes in airfoil shape, which in turn help to control static and dynamic aeroelastic problems [6]. Numerous studies and comprehensive research programs run in recent years have demonstrated the feasibility of piezocomposite materials for a helicopter active rotor design and tilt rotor vehicles. Active twisting of blades with embedded piezoceramic transducers and shape mem-ory alloy actuated tabs and flaps systems have been examined in order to actively control vibration and blade-vortex-induced noise. Results have shown reductions of vibratory loads of about 80%, as well as reductions up to 6 dB in blade-vortex interaction and in-plane noise. Moreover, the impact of the active flap on rotor performance, rotor smoothing, and control power has been demonstrated [4,33]. Tests performed by Boeing on V-22 tilt rotor aircraft have revealed the potential of active flow control to minimize hover downloadlift during takeoff conditions and to improve the payload lift and angle-of-attack capability [4,11]. When compared to fixed-wing aircraft, helicopters appear to show much higher potential for a major payoff with the application of smart structures technology. The given above numerous examples of piezo-composite beamlike structures confirm the current relevance and timeliness of this subject matter and give a great stimulus and motivation to continue research in this area.
While currently there are many types of active materials available, like shape memory alloys, electrostrictives, and magnetorheological fluids etc. piezoelectric materials remain the most widely used "smart" ones. This is primarily due to their strong voltage-dependent actuation viability. Secondly, piezoceramics are capable of interacting with dynamic systems at high frequencies ranging up to even 1 MHz. Finally, the crucial feature of piezoceramic materials is their large energy density. Therefore, they can be effectively adapted to, e.g. energy-harvesting devices that can generate high electric potentials or may be used to supply power to wireless sensors and low-power electronics [20].
Composite materials are perfect candidates for implementing the concept of structural functionality because of their multiphase nature, suitable manufacturing technology, low density and inherent directional mechanical properties [6,28,39]. Studies on multifunctional structures started in late 1980s/early 1990s from modelling isotropic material systems followed next by the research on piezoelectric laminated composite materials. Smith et al. [31] proposed the constituent equations representing the behaviour of piezoelectric bimorphs for arbitrary mechanical boundary conditions. Weinberg [42] studied the piezoelectric multimorph for sensing and actuation by means of the Euler-Bernoulli beam model. Closed-form solutions for force, displacement and charge developed in piezoelectric beams were derived. His work was later extended by Tadmor and Kósa [36] by accounting for the effect of strain on the electric field in the active material domain. This was done by a simple correction factor to the moment of inertia of the piezoelectric layers. A review of different theories used for modelling and analysis of piezoelectric laminates at that time was published by Gopinathan et al. [9].
It is interesting to note the most papers studying structural behaviour of piezoelectric materials and their possible applications to multifunctional structures are focused on the classical d 33 or d 31 piezoelectric effects. This trend comes from the fact these phenomena are easily observed and can be directly exploited in axial and transverse deformations modes of structures. However, in some piezoceramics, the d 15 piezoelectric shear coefficient is much higher than the d 33 and d 31 ones. Therefore, transducers based on these materials should be operated in shear deformation modes since the electromechanical coupling coefficient of the piezoelectric material is one of the most significant parameters affecting energy conversion. Studies on this topic were done, e.g. by Dietl et al. [7]. Authors developed a mathematical model of a piezoelectric sensor based on the Timoshenko beam theory. Next, it was used to examine the frequency response of vibration-based energy harvesters; finally the obtained results were compared to ones coming from the classical Euler-Bernoulli model. Also Malakooti and Sodano [23] studied the shear deformation modes of piezoelectric materials. The proposed model was used to predict the electric power output from a cantilever piezoelectric sandwich beam under base excitations. It was shown that the energy harvester operating in the shear mode was able to generate up-to 50% more power if compared to the standard transverse mode operating one. Capabilities of shearable beam with embedded piezoelectric actuation for beam shape control were studied also by Voß and Scherpen [41]. Also mixed theories adopting the geometrically linear Bernoulli-Euler hypothesis for the mechanical components and a first-order theory for the electrical variable can be found in the recent literature [15,35].
Conclusions summarizing these papers indicate that the Euler-Bernoulli beam model severely over-predicts the structural behaviour of tested functional structures, especially at low frequencies and low length-towidth aspect ratios. This over-prediction is of particular importance in case of composite materials that exhibit relatively low shear stiffness when compared to classical isotropic materials.
Another important aspect in proper modelling the piezoelectric functional structures is mutual interaction of mechanical and electrical domains of the system. The initial studies did not consider the field coupling effects thus both mechanical and electrical domains were treated independently. Within this framework, the mechanical properties of the hybrid structure resulted solely from the combination of stiffnesses of two constituent materials. The functional nature of the structure was captured by introducing dynamic loading acting at the requested position on the structure. This simplified approach provides only approximate results, and therefore, in some cases like high frequency vibrations or thick piezoelectric material layers, significant errors arise [37].
The coupled electromechanical model of the smart composite materials was proposed by Mitchell et al. [24]. Authors presented a method of enhancing plate theory to account for the charge equations of electrostatics. Further studies on coupled domain models of piezoelectric composite plates were done by Li et al. [21] and Zhou et al. [45]. Chattopadhyay et al. [5] elaborated a higher-order displacement field model of a plate to investigate the behaviour of smart helicopter rotor blades. However, the proposed theory was based on a three-dimensional approach so the final equations were too complicated to be solved analytically and the finite element method was used to get the solution. Later, authors extended their model by addressing the nonlinear electromechanical coupling effect. To this aim, the polarization versus electric field hysteresis was taken into account and a new material constant was introduced to explicitly formulate the nonlinear constitutive governing equations of the structure. This paper is a continuation of the previous authors' research on dynamics of thin-walled composite beams. In particular, in reference [18] by Latalski et al. where the general partial differential equations of motion of a thin-walled composite beam clamped to the rigid hub were derived, the nonclassical effects like material anisotropy, rotary inertia, transverse shear, arbitrary pitch angle, and hub inertia were considered. The ordinary differential equations were formulated together with the associated orthogonality condition for the specific case of circumferential asymmetric stiffness (CAS) of the cross section.
In the following research [16], Latalski proposed a mathematical model of a beam with integrated piezoelectric layer accounting for both full electromechani-cal coupling effects as well as the higher-order constitutive relations. The detailed derivation procedure by means of Hamilton least action principle was presented. With respect to former research a new governing equation was formulated representing the electric field distribution within the piezoceramic domain. It has been shown the electromechanical coupling in the structure comes from the shear deformation in flapwise bending plane. This model has been adopted in later research for the analysis of a rotating thin-walled composite beam with embedded piezoelectric layer [17]. In studies, a specific case of circumferentially asymmetric stiffness lamination scheme that exhibits flapwise bending and twist mode elastic coupling was considered. Numerical results for system free vibrations were obtained to investigate the natural mode shapes and the electrical field spatial distribution depending on the system rotation speed and laminae fibre orientation angle.
In this research, the same fully coupled mathematical model of functional structure is adopted to study the dynamics of a rotating composite beam featuring bending to bending coupling. This dynamical behaviour is achieved by adopting the circumferentially uniform stiffness (CUS) along the profile cross section. In the formulation, the higher-order constitutive relations with respect to electric field variable are considered. In the numerical analysis, the forced responses of the system under zero and nonzero mean value harmonic excitation are investigated. Regarding the referenced papers and to the best of authors knowledge, these aspects of CUS laminated multifunctional beams with nonlinear piezoelectrics have not been studied in detail yet.

Structural model and problem formulation
Let us consider an elastic single-cell thin-walled beam with box-shaped cross section. The beam is clamped to the rigid hub at an arbitrary presetting angle θ -see Fig. 1a, b. The hub has radius R 0 and is driven by the external torque T ext,z resulting in rotation of a structure about fixed vertical axis C Z 0 .
The considered beam is made of unidirectional fibre-reinforced composite material obeying Hook law. Apart from regular UD material plies, two additional structural layers of transversely isotropic piezoceramic material are embedded onto the beam. These ones are placed on the outer faces of the two beam walls that are perpendicular to the flapwise bending plane (i.e. the plane of lower specimen stiffness-the xz one of the local coordinate system x yz , see Fig. 1c), and cover the full span of the beam. Moreover, it is assumed the active piezoelectric material is poled in the thickness direction. Thus, this configuration is a typical geometry corresponding to '3-1' mode operating transducers based on active piezoelectric materials. Moreover, in the performed analysis, it is assumed that the piezoceramic layers are not covered by electrodes and the poles of the piezoceramics are considered to be opencircuited.

Piezoelectric material model
The equations of motion of the structure are derived following the extended Hamilton principle of the least action where J is the action, T is the kinetic energy, U is the potential energy including mechanical (U m ) and electrical components (U e ), and the work done by the external loads (driving torque T ext,z ) is given by the W ext term.
Detailed derivation procedure of a six-DOF model of the beam with an arbitrary presetting angle θ and any circumferential lamination scheme has been presented in former authors papers [17] and [18]. These papers include also the list of adopted kinematic and structural assumptions to the mathematical model as well as the discussion on their significance. In particular, the second paper provides an efficient approach to properly model the full two-way electromechanical coupling interaction observed in piezoelectric structures. The relevance and suitability of this mathematical formulation is further enhanced by postulating the higherorder constitutive relations with respect to electric field In the above relations C stands for the second-order piezoceramic elasticity tensor at constant electric field, e is the tensor of piezoelectric coefficients, ξ is secondorder permittivity tensor,b is effective electrostrictive constants tensor, χ is third-order electric susceptibility tensor. Moreover, the variables σ , ε and D stand for stress and strain tensors and electric displacement vector, respectively. Finally, the electric field vector is denoted by E, and it contains just one component E 3 as results from the transducer topology and To capture the two-way electromechanical coupling effect, the electric field variable E 3 is assumed to be a spatially dependent one More details on the mathematical modelling of piezoceramic two-way coupling effect can be found in [37,45], as well as in previous study [16]. The considered constitutive equations postulate second-order nonlinearities expressed in electric field magnitude, rather than the electric field itself. Moreover, a third-order nonlinearities are not taken into account. In this way, the proposed model is capable to predict a softening/hardening curves with their peaks changing linearly with the excitation magnitude. This comes from the fact the linear change in peak response frequency with increased excitation level has been observed in experimental research on both stiff and soft piezoelectric materials-note papers by Usher and Sim [40], Mahmoodi et al. [22] and Leadenham and Erturk [19] among others.

CUS laminate configuration
One of the crucial features of structures made of fibre-reinforced laminates is the ability to tailor their mechanical properties to meet specific design requirements. This aim can be achieved by adjusting the orientation of the fibres in the composite material or/and through the use of the variable stiffness concept. This later idea implies the stiffness properties of the structure vary spatially throughout its volume. In case of thin-walled beam elements, this usually corresponds to changes in material thickness or laminate properties along the specimen span or along the cross-section midline (e.g. different orientation of reinforcing fibres in the subsequent segments of the cross section). When considering this second case, two designs are most typically adopted-they are referred to as circumferentially uniform stiffness (CUS) and circumferentially asymmetric stiffness (CAS) configurations [30]. For a thinwalled beam of arbitrary but rectangular cross section, the CUS lamination scheme can be generated by skewing the ply angles in the top and bottom walls (flanges) according to the formula α(z) = α(−z) and in the lateral walls (webs) as α(y) = α(−y), respectively. The reinforcing fibres orientation angle α is measured with respect to the in-plane axis (e.g. s) of a wall related coordinate system sxn ; note Fig. 2 for reference.
Accepting the discussed circumferentially uniform stiffness lamination scheme causes the general system of six (three translations and three rotations) fully coupled governing equations of the beam to be split into two independent subsystems. One of these subsystems describes the coupled in-plane and out-of-rotation plane shearable bendings of the blade. The second governing subsystem represents the coupled axial (specimen extension) and torsional dynamics of the beam. It should be noted that the circumferentially uniform stiffness ply angle configuration is preferred in the design of helicopter blades and tilt rotor aircraft [14].

Equations of motion
Presuming the discussed above circumferentially uniform stiffness of the profile as well as beam presetting angle θ = − π 2 the system of partial differential governing equations for the structure under consideration is as follows -the hub with boundary conditions with boundary conditions with boundary conditions where In the above equations, variables v 0 and w 0 represent displacements of the beam cross-section origin O along the y and z axes in the local coordinate system x yz -see Fig. 1b, while ϑ y , ϑ z are the corresponding rotations. Moreover, the coefficients b 1 , B 4 and B 5 are translational (b 1 ) and rotational inertias (B 4 , B 5 ), respectively. Terms a i j represent different stiffness coefficients of the specimen. In particular, the pair a 22 , a 44 represents chordwise (local x y plane) and pair a 33 , a 55 flapwise (local xz plane) bending and shear stiffnesses, respectively. Moreover, there are two additional stiffness coefficients a 25 and a 34 that lead to mutual coupling of beam bending deformations. The term a 25 is chordwise bending to flapwise transverse shear and a 34 is flapwise bending to chordwise transverse shear coupling stiffness. The last Eq. (8) of the system represents the distribution of the electric field along the span of the beam, and it is directly related to the lead-lag plane shear variable ϑ y Eq. (7).
Commenting on the given above equations of motion, it is worth to emphasize the mutual coupling within the system of beam governing equations is achieved by both shear deformation angles ϑ y and ϑ z . Therefore, within the framework of simplified unshearable beam theory (like Euler-Bernoulli one), not only the shear deformation Eqs. (5) and (7) vanish but the remaining displacements ones get decoupled. The relations are as follows: with boundary conditions v 0 This formulation is derived by extracting a 44 (v 0 + ϑ z ) and a 55 (w 0 − ϑ y ) from Eqs. (4) and (6), respectively, replacing these expressions in the remaining equations of motion Eqs. (5) and (7) and finally substituting ϑ y = −w 0 and ϑ z = −v 0 .
The definitions and ultimate formulas of the individual stiffness coefficients are given as follows In the above expressions, the additional superscripts (w) or (f) at K stiffnesses correspond to the part of the cross-section perimeter belonging to the web or to the flange, respectively. The coefficients K i j are expressed in terms of reduced two dimensional stiffnesses A, B, and D of the classical laminate theory It should be observed for the sections of the perimeter corresponding to lateral walls (webs) that the coefficient K 44 reduces just to bending stiffness D 22 , and two other ones, namely K 14 and K 24 , vanish. This is the result of material symmetry with respect to the wall midline that forces B 12 = B 22 = B 26 = 0. However, in case of flanges, the present outer piezoelectric layers break the wall cross-section symmetry and thus respective membrane-bending coupling stiffnesses B i j are different from zero.
The magnitudes of stiffness coefficients depend on reinforcing fibres orientation in UD laminae. In case of the bending stiffnesses a 22 and a 33 the relation is straightforward. They increase monotonically when the fibre orientation angle approaches the beam span direction. However, in case of the shear stiffnesses a 44 and a 55 as well as the coefficients a 25 and a 34 accounting for in-plane/out-of-plane shearable bendings coupling the relations are more complicated; see Fig. 3.
Both shear deformation coefficients a 44 and a 55 exhibit very similar behaviour with respect to reinforcing fibre orientation. Initially, their magnitude increases until α ≈ 54 • and ≈ 68 • , respectively. Next, the curves decrease rapidly reaching global minimum at α = π 2i.e. when fibres are oriented along the beam span.
Regarding the coefficient of chordwise bending to flapwise transverse shear coupling (a 25 ), one observes it reaches its global extrema at orientations α ≈ 74 • or α ≈ 106 • and stays antisymmetric with respect to beam span direction. For limit cases of fibres oriented along the cross-section circumference or along the beam spam, it reaches zero value. However, the qualitative behaviour of the second coupling coefficient a 34 (flapwise bending to chordwise transverse shear) is different. Apart from global extrema corresponding to angles 28 • and 152 • and zero values at the limits 0 • and 90 • , one can spot two additional roots corresponding to α ≈ 16 • and α ≈ 114 • . It means, at this specific orientations, the bending of the speci- Fig. 3 Governing equations stiffness coefficients with respect to reinforcing fibre orientation angle α (measured from circumferential axis direction s-see Fig. 2). Data used for this analysis correspond to specimen data used in numerical studies given in Sect. 4 men in the plane of rotation, namely lead-lag bending, do not contribute to shear deformation in the orthogonal plane. The qualitative difference of the a 34 curve if compared to a 25 one can be explained by the presence of the piezoceramics within the flanges of the cross section breaking the symmetry of material distribution with respect to these walls midline. Due to this fact, additional local membrane-bending couplings arise (B 12 , B 22 , B 26 = 0) contributing to the global specimen behaviour.

Solution procedure
To facilitate the solution and extend the generality of the problem formulation one starts from rewriting the orig-inal equations of the system (4)-(8) into their dimensionless form. To this aim, a dimensionless spanwise abscissa η = x /l (η ∈ 0, 1 ) is introduced as well as dimensionless time τ = ωt, where t is physical time and ω = √ a 33 (0) /b 1 l 4 . The parameter a 33 (0) corresponds to the flapwise bending stiffness presuming that the laminate fibres are oriented at α = 0 • with respect to profile circumferential axis s.
Next, the reformulated system of governing equations is transformed into ordinary differential form. To this aim the Extended Galerkin Method is used. At the first stage the simplified linear model of a nonrotating beam without piezoceramic layer is considered and the corresponding eigenvalue problem is solved. Within this procedure, the consistent admissible functions which have to fulfil all the geometric boundary conditions while not violating complementary boundary conditions [2,29] were used. Based on previous authors' studies [18], Duncan polynomials modified by Karanamoorthy have been adopted to approximate the displacements v 0 , w 0 and subsequent sin( 2i−1 2 πη), i = 1 . . . k, functions to represent the rotations ϑ y and ϑ z . As a solution of the eigenvalue problem, a series of complex mode shapes V (η), W (η), Y (η), Z (η) and corresponding natural frequencies have been determined. These mode shapes have been used to represent the original problem unknown variables by separating space and time dependent functions, respectively In the above relations, j stands for the mode order and the variables q j (τ ) are the corresponding generalized coordinates representing time dependent system behaviour, N is the number of mode shapes obtained from the eigenproblem.
The accuracy of this solution method, and the proposed analytical model has been verified in the previous authors' research [18] by comparing the analytically obtained natural frequencies to the outcomes of the FE analysis. For the discussed therein case of CAS lamination configuration, the margin of error was less than 5%.
Having found the individual components of beam modes shapes, the reduction in the system of partial differential equations to ordinary ones can be done. Combining all five equations representing the dynamics of the beam (subsystem of Eqs. (4)- (8)) and making use of the derived orthogonality condition one finally arrives at the system of two dimensionless equations where μ is dimensionless driving torque supplied to the hub and symbol J h denotes the relative mass moment of inertia of the hub calculated with respect to the beam inertia. In the general case, the torque is expressed as a sum of a constant component μ 0 and a periodic function of time t. Thus μ(t) = μ 0 + ρ cos ωt (16) where ρ and ω are the amplitude and the frequency of the excitation, respectively. Moreover, the factors ζ h and ζ i are hub and beam viscous damping coefficients and their magnitude has been estimated during test-stand experiments. The subscript i of the generalized coordinate q i indicates the mode order to be considered in the performed analysis.
The geometric characteristics of the rotating beam as well as composite and piezoceramic material data used in the subsequent numerical simulations are collected in Table 1. For the composite, the subscripts 1, 2 and 3 refer to parallel and transverse to the fibre directions respectively (the standard classical laminate theory assignment). The embedded active element is made of PZT-3203HD ceramic material and its properties are collected from the papers by Rao et al. [27] and by Kapuria et al. [12].

Numerical results
Numerical analysis of the rotating piezo-composite blades is performed around selected resonance zones of the rotating hub-beam system. By neglecting in Eqs. (15) the nonlinear terms, we can determine the first linear natural frequency of the structure. Therefore, the fundamental frequency of the linearized hub-blade system is approximated as Studying this relation, one may notice in the limit case of hub mass moment of inertia J h approaching infinity the natural frequency of the rotor tends to ω 0 = α 11 + α 13ψ 2 . In particular case when there is no rotation it simplifies to ω 0 = √ α 11 which corresponds to the natural frequency of a regular cantilever beam. We investigate dynamics of the system imposing periodic torque μ to the hub. Solutions of the reduced model are obtained numerically; a continuation method together with bifurcation analysis is performed for selected system parameters [8].
For the analysis two configurations of the laminate corresponding to reinforcing fibres orientations α = 40 • and α = 70 • ale selected; see Fig. 2 for notation of fibre direction. For each layout, the reduction to single mode is performed and dynamics of the system around the fundamental (mode 1) and the second (mode 2) natural frequency is investigated. The values of dimensionless coefficients obtained for both configurations after a system reduction are presented in Tables 2, 3, 4 and 5. Moreover, the shapes of investigated modes coming from the eigenvalue problem solution are presented in Fig. 4 with separation of individual components contribution.

Fundamental frequency
At first we consider a case when the structure is excited periodically by torque with constant component equal to zero, μ 0 = 0 (note Eq. (16) for reference). This scenario corresponds to the case of non-fully rotating structure but performing oscillations around its zero position due to zero mean value excitation μ(τ ) = ρ cos ωτ . The response of the blade q 1 and the angular velocity of the hubψ are computed for different torque amplitudes ρ and varied frequency ω in the neighbourhood of the resonance zone.
Results of the analysis for configuration corresponding to α = 40 • when the mode 1 is excited are presented in Fig. 6a, b. The values of coefficients used in numerical simulation are given in Table 2.
One can clearly observe the softening phenomenon near the resonance zone due to the nonlinear characteristics of the PZT layers. It should be noted the backbone curve reveals change in peak response frequency with increased excitation level. As commented in Section 2, this effect is related to the presence of second-order nonlinearities expressed in variables magnitude, rather than the variables itself. The softening effect occurs both for the beam Fig. 6a and for the hub Fig. 6b, despite the fact that active layers giving rise to system nonlinearity are embedded into the beam only. This layout of the composite reinforcing gives significant nonlinear effect observed even for relatively small excitation level ρ = 0.1 (black curve) and is increasing if torque amplitude is enlarged (ρ = 0.2-blue curve and ρ = 0.4-green curve).
Changing the composite configuration into α = 70 • layout (values of the corresponding coefficients are collected in Table 4) a similar set of fundamental resonance curves has been obtained-see Fig. 7. However, one notices the nonlinear softening effect is now much less pronounced comparing to the previous case α = 40 • . To observe the nonlinear phenomenon the excitation amplitude has to be increased up to ρ = 0.4 (green curve); for the record also ρ = 0.5 (red curve) ρ = 0.6 (magenta curve) cases are plotted. This observation can be easily explained by the significant difference in flapwise bending (in lead-lag direction, i.e. rotation plane) stiffness a 33 of the beam for these two specific casesrefer to Fig. 3. Moreover, at the orientation α = 70 • , there is a much stronger coupling of flapwise bending to the chordwise one (out-of-plane bending)-note coefficient a 34 when compared to α = 40 • case. This is accompanied by the similar strong reverse coupling (a 25 coefficient), but this modal component is barely present-see Figs. 4 and 5.
When analysing hub responses in both fibre configurations one can notice a zone where the hub vibrations  Fig. 6 and ω = 5.07 on the Fig. 7). This phenomenon is apparently occurring for all studied levels of excitation ρ. It corresponds to dynamic absorption of hub vibrations when the total oscillation energy is directed to the blade and keeping the hub at rest. Detailed analysis reveals the effect occurs at the excitation frequency corresponding to the natural frequency of the cantilever beam ω 0 = √ α 11 .
The direct comparison of both CUS 40 • and CUS 70 • layout variants is shown in Fig. 8. Changing the orientation of composite reinforcing fibres from 40 • to 70 • results in the shift of the resonance zone towards higher frequencies (i.e. to the right). Moreover, this stiffening effect is manifested in reduction of the resonance curve inclination and the magnitude of amplitudes as well.
In order to check the influence of the hub on the structure dynamics, we computed resonance curves for different magnitudes of hub mass moment of inertia J h . In Fig. 9, we present resonance curves for mode 1 and CUS 70 • configuration taking into account the relatively light hub (mass moment of inertia J h = 5 -black curve), next the cases of J h = 20-blue curve and J h = 50-green curve, and finally for the very heavy one J h = 100-red curve. One can observe the increase in the hub mass moment of inertia leads to shifting of the resonance zone into lower frequencies presuming the driving torque magnitude imposed to the hub is constant. The sharp decrease in the vibration amplitude is to be noted and-in consequence-the reduction in the nonlinearity effect. However, the already mentioned phenomenon of hub absorbtion is still present and occurs independently of the hub mass moment of inertia (Fig. 9b) and at the same frequency corresponding to natural vibrations of the cantilever beam. Now let us consider just for reference the black resonance curve (Fig. 9) computed for J h = 5 and the excitation amplitude ρ = 0.6. We intend to study the response of the system against J h parameter for frequencies before, inside and after the resonance zone corresponding to this curve. To this aim, the response plots for the cases ω = 5.2 corresponding to subresonance zone, ω = 5.4 and ω = 5.6 inside the reso-  Table 2) and CUS 70 • -blue curve (data in Table 4 nance zone as well as ω = 6.0, ω = 6.5, ω = 9.0 after the resonance have been computed. Appropriate plots are presented in Fig. 10. We note that the dynamics of the system changes qualitatively passing from the softening to hardening nature when the J h parameter is decreased. Moreover, the amplitudes of hub response (Fig. 10b) are higher for lower J h values and also if excitation frequencies are higher. In the case of beam response (Fig. 10a), the tendency is slightly different. Beam amplitudes increase as well, but for very low hub inertia, the tendency is opposite (orange curve in Fig. 10a).

Second natural frequency
The next mode of our interest is the second mode (mode 2). Coefficients of the discretised equations of motion corresponding to the α = 40 • layout are gathered in Table 3 and for 70 • configuration in Table 5.
We start the analysis from plotting the response curves for five different driving torque amplitudes, namely ρ = 0.1, ρ = 0.2, ρ = 0.4, ρ = 0.6 and also for the highest one ρ = 1.1. The curves are presented in Fig. 11 for layout CUS α = 40 • and in Fig. 12 for CUS 70 • case. It can be easily observed the amplitudes for the second mode are a few orders of magnitude smaller then those obtained for mode 1. Furthermore, the resonance curves demonstrate tendency to hardening effect which can be visible for very large, although non-realistic, external excitations which are not plotted here. Moreover, the dynamic absorption effect characteristic for fundamental frequency is not present any more since the resonate curves of the hub response decline monotonically. Comparison of the mode 2 resonance curves generated for CUS 40 • and CUS 70 • configurations is presented in Fig. 13. Since this mode can be activated in system response for higher excitation frequencies, therefore, in order to have comparison with fundamen-  Table 5. (Color figure online) tal mode 1, the analysis is performed for the same excitation level ρ = 0.6. In contrast to the similar comparison prepared for mode 1 in the resonance zone of mode 2, amplitudes are higher for CUS 70 • configuration. This behaviour can be explained by studying the mode plots in Figs. 4 and 5. For both configurations, the dominating deformation is related to out-of-plane bending (v(η) and ϑ z (η)), and thus the normalization of mode shapes to 1 has been done with respect to ϑ z in both cases. However, for the layout corresponding to α = 70 • , there is an additional component coming from lead-lag bending that is barely present for α = 40 • set-up. When comparing the response of the hub, both configurations reveal almost the same behaviour, and the curves overlap as presented in Fig. 13b.

Influence of rotating speed: full rotation case
In all the analysis performed so far the rotor has been excited by periodic torque with zero mean value. To capture a full rotation impact, we add the constant com-  Table 3) and CUS 70 • -blue curve (data in Table 5 Commenting these plots, one observes the evident centrifugal stiffening effect resulting in shift of resonance curves to the right. Moreover, due to the nonlinear properties of the PZT layers, the resonance curves still demonstrate a softening effect; however, it is less pronounced at higher rotating speeds (higher μ 0 value), again due to centrifugal stiffening of the blade. The softening, however, is much better demonstrated for mode 1 and CUS 40 • configuration and as expected less visible for CUS 70 • configuration. The change of nature of constant component μ 0 and for fixed amplitude ρ and excitation ρ are shown in Fig. 15a for beam response and in Fig. 15b for angular velocity of the hub. The increased value of constant component μ 0 makes a shift of the response; however, the nonlinear nature of the solution is clearly demonstrated.

Conclusions
The analysed model of the rotating hub with the composite thin-walled beam with embedded nonlinear piezoelectric layer demonstrates softening effect near its resonance zones. The nonlinear phenomenon is observed for the beam as well as for hub response due to the coupling of the hub to the blade by inertia terms. The intensity of the nonlinearity depends on the studied mode and the layout of reinforcing fibres within the discussed CUS configuration framework. Referring to the fundamental resonance zone the nonlinear nature of the response curve is much more pronounced for CUS 40 • layout than for the CUS 70 • case. This discrepancy comes from the significant difference in specimen bending/shear stiffness. In case of the second resonance zone, the oscillations are a few orders smaller if referred to fundamental resonance response. Furthermore, for the second mode, a small stiffening effect is observed with larger oscillation for CUS 70 • configuration, in contrast to the results obtained for first mode. The mass moment of hub inertia changes the dynamics of the rotating system. Depending on the excitation frequency response of the beam and hub as a function of hub inertia may vary from softening to hardening as presented in Fig. 10.
Another interesting phenomenon observed for the first mode and both configurations is the absorption of hub oscillations by beam motion. This dynamic absorption takes place for frequency equal to natural frequency of cantilever beam ω ≈ √ α 11 . In this singular zone oscillations of the hub are close to zero with vibration localized in the beam. For the second mode, the absorption is not present. The presented reduced model allows fast and efficient parametric analysis and optimal design of the rotating composite structure with active blades comparing to for example highcomputational cost and time consuming finite element studies.