Nonlinear vibrations and time delay control of an extensible slowly rotating beam

Nonlinear dynamics of a rotating flexible slender beam with embedded active elements is studied in the paper. Mathematical model of the structure considers possible moderate oscillations thus the motion is governed by the extended Euler–Bernoulli model that incorporates a nonlinear curvature and coupled transversal–longitudinal deformations. The Hamilton’s principle of least action is applied to derive a system of nonlinear coupled partial differential equations (PDEs) of motion. The embedded active elements are used to control or reduce beam oscillations for various dynamical conditions and rotational speed range. The control inputs generated by active elements are represented in boundary conditions as non-homogenous terms. Classical linear proportional (P) control and nonlinear cubic (C) control as well as mixed (P-C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P-C$$\end{document}) control strategies with time delay are analyzed for vibration reduction. Dynamics of the complete system with time delay is determined analytically solving directly the PDEs by the multiple timescale method. Natural and forced vibrations around the first and the second mode resonances demonstrating hardening and softening phenomena are studied. An impact of time delay linear and nonlinear control methods on vibration reduction for different angular speeds is presented.


Introduction
Slender beam-like elements play an important role in engineering and structural design. Typical applications might be cranes, aircraft wings, diving boards at swimming pools, overhang structural elements of civil engineering like masts or roof supports. When rotary motion of the structure is considered further examples might be wind turbines blades, helicopter blades, aircraft propellers, etc.
Advances over the years in composite materials technology as well as increasing demands particularly from aeronautics and off-shore engineering have stimulated the extensive use of light and flexible elements. Predominantly, they are subjected to large elastic deformations that affect the precision and stability of the structure motion. Moreover, large deformations com-bined with low structural damping may lead to fatigue damage and shorten the lifespan of the design. Therefore, the attenuation of large vibrations observed in highly flexible beam-like elements is a problem of primary importance.
Driven by practical needs as well as theoretical challenges, efficient and accurate modeling of flexible beams dynamics and their control have received great attention in the literature. Depending on complexity, current compact beam theories can be divided into three main groups: (a) un-shearable theory including the classical Euler-Bernoulli one, (b) shear deformable models-e.g., Timoshenko theory, the third-order shear theory etc., and (c) three-dimensional beam theories capturing different phenomena neglected by the former two approaches. Typically, each of commonly accepted theories can be formulated within linear or nonlinear framework. Linear models are very useful in a case of relatively stiff systems, performing small oscillations. For thin slender flexible elements that undergo large deformations, use of nonlinear theories is much more favorable [32].
The dynamics of a rotating beam system within nonlinear framework was examined by several researchers. Weidenhammer [44] studied rotating beams by adopting a (not-complete) nonlinear theory of Bernoulli-Euler beams; the governing equations of motion were derived applying Hamilton's principle. The influential work presenting a comprehensive nonlinear beam analysis was published by Crespo da Silva and Glynn in [8,9]. The nonlinear-order three differential governing equations were derived by Hamilton's method accounting for contributions resulting from nonlinear curvature as well as nonlinear inertia. It was shown that both these effects may had a significant influence on the nonplanar moderately large oscillations of the system. This initial research was later enhanced to consider complex deformations involving flexure along two principal directions as well as torsion [7]. Later, author extended his analytical model to a rotating beams case with the main application to helicopter rotor blades [10].
Hamdan and El-Sinawi in [16] studied the inextensible nonlinear Euler-Bernoulli beam model accounting for relatively large planar deformations and exact expression for the beam curvature. Influence of a setting angle and other selected structural parameters on rotor response characteristics for a prescribed hub torque scenarios was discussed. It was shown that for a soft base and a low preset angle unstable vibrations of the rotor might have occur.
Fenili et al. [13] presented a nonlinear mathematical model of a flexible beam-like structure in slewing motion. Authors used the method of multiple timescales to find analytical solutions in selected primary and secondary resonance states. Recently, the importance of nonlinear effects coming from the large displacement oscillations observed in rotating inextensible beams leading to rich dynamic behavior has been presented by Thomas et al. in [38]. The performed analysis demonstrated both softening or hardening frequency response characteristics dependent on the resonance order. Interestingly, the observed softening effect originating from geometric nonlinearities prevailed the centrifugal stiffening phenomenon. Also Tian et al. [39] presented a general formulation for nonlinear vibration analysis of rotating beams. The numerical solutions demonstrated that Coriolis effect could essentially change dynamics of the hub-beam system in the case of small hub radius, large beam slenderness, and high angular velocity.
When considering structures made of composite materials the kinematic relations become much more complicated since geometrically originated couplings are accompanied by strong directional properties of the constituent material. Studies accounting for coupling between flexural and longitudinal vibrations were published recently by Lenci and Rega [27] and next by Babilio and Lenci [3,4]. Besides the mentioned bending-axial couplings, authors took into account the shear effect and imperfect boundary conditions. The importance of the adopted definition of geometric curvature was discussed by Lenci and his group in [24,26]. The proposed strict nonlinear shearable beam model was solved analytically by attacking directly the nonlinear partial differential equations and the results were verified by the Abaqus/CAE finite element analysis [18,19]. The results showed softening versus hardening dichotomy in the resonance curves and also strong interactions between flexural and longitudinal (axial) vibrations leading to internal resonances occurring for specific combinations of beam to axial end spring stiffness ratios. Analytical and numerical predictions were accompanied by experimental studies on a slender beam-spring system subject to kinematic excitation in [20].
Widely reported large amplitudes observed in oscillations of highly flexible beam-like elements motivated practicing engineers and scholars to extensively study the methods of oscillations suppression and to improve structural positioning accuracy. The fundamental idea was the application of any control methods, that is passive or active ones. Possibly the most promising results can be obtained by closed-loop control strategy achieved by the feedback of the system state recorded by sensors to drive actuation forces/moments generated by, e.g., piezoelectric devices [14].
With reference to large oscillations of the cantilever beams the problem of structural control by means of piezo sensor-actuator layers was studied by Rechdaoui et al. in [36]. To mitigate vibrations the proportional and time derivative potential feedback control was formulated. It was shown by tuning the control parameters and gains the nonlinear dynamic behavior of the structure could be actively suppressed. The problem of large amplitude oscillations of beams and their control was studied also by Nguyen et al. in [33] with reference to marine risers represented by long tensioned Euler-Bernoulli beam undergoing bending in two orthogonal planes. Based on the set of equations, the boundary controller applied at the top end of the riser was designed using Lyapunov's direct method. The proposed algorithm was used to effectively stabilize the riser at its equilibrium position. Similar problem was studied by Do in [12]. The structure was modeled analytically as extensible and shearable slender beam undergoing large 3D translational and cross-section rotary motions. The control design and stability analysis were based on two Lyapunov-type theorems developed for a class of evolution systems in Hilbert space.
The control of complex flexural-torsional vibrations of a rotating composite beam was studied also by coauthors of this paper in [43]. It was proposed to use the saturation control method to reduce vibrations of flexible beam clamped to the hub that was excited by the prescribed harmonic torque. Conducted tests showed narrowing down the zones of effective suppression of beam vibrations while increasing the rotating speed of the structure. Other control methods exploited by piezo actuator devices for beam oscillations suppressing include LQR technique [5], LQG [30], P control [28], PID control [35]. Further reading on beam vibrations and various control laws can be found in a comprehensive book by He and Liu [17].
A significant effort and research works are also focused on accounting for the effect of time delays as present in all actual dynamical systems with system control. Daqaq et al. [11] examined the effect of feedback delays on the nonlinear vibrations of a piezoelectric actuated cantilever beam and analyzed the effect of feedback delays on a blade when subjected to harmonic base excitations. Alhazza et al. [1] investigated the effect of time delays on the stability, amplitude and frequency-response characteristics of a beam. The authors found that even the small time delays could completely altered the behavior and stability of the parametrically excited beam and might lead to unexpected structural phenomena. Liu et al [29] studied the piezoelectric based optimal delayed feedback control method when applied to large amplitude nonlinear vibrations of a beam. The time-delayed feedback control to reduce the nonlinear resonant vibration of a piezoelectric elastic beam was studied also by Peng et al. in [34]. Authors tested three different single-input linear time-delayed feedback control methodologies, namely displacement, velocity and acceleration timedelayed feedback. It was shown the time-delayed feedback control could act as a vibration absorber at specific values of time delay magnitude.
It is interesting to note the vast majority of papers proposing the use of piezoelectric transducers for feedback structural control adopt the simplified mechanical model of the piezoelectric active elements and their mutual interaction with hosting structure. Most of them capture just the phenomenological behavior of the combined system involving functional material and the host member. These models consider solely the converse piezoelectric effect where the commanded actuation strain is transferred as a shear force to the master structure along the actuator-subsystem interface. This force-since located off beam mid-line-generates the control moment.
The discussed above new results obtained for nonlinear stationary beams, followed by the observed interactions between transversal and longitudinal vibrations and the dichotomy of softening and hardening response behavior as well as the potential of feedback control methods motivated authors of the present contribution to a more detailed study of the problem under discussion. In particular, the nonlinear phenomena observed for the rotating inextensible beams [38,39] and provided conclusions suggested to analyze a more general case of an extensible beam. As reported in the literature, the effect of longitudinal vibrations can be even more important if the rotating beam carries a heavy tip mass. A preliminary study of a rotating extensible beam model which takes into account transversal and also longitudinal vibrations was presented in [41] and then extended in [42]. In the current work, we consider the configuration of the system corresponding to an arbitrary preset angle of the beam and added tip mass. We also modify definition of the beam curvature adopted in the mathematical model, according to the comments presented by Lenci in [25]. The model studied in this paper considers an active structure with embedded piezo-layers and a boundary control method with time delay. On the basis of the proposed complete active beam model we determine approximate solutions on the basis of the multiple timescales method applied directly to the partial differential equations with time delay and associated dynamic boundary conditions.
Following the given above comments the rest of this manuscript is organized as follows: In Sect. 2, the dynamic nonlinear model of a highly flexible beam is derived using the Hamilton's principle and then solved by the multiple timescales method in Sect. 2.2. Next, identification of model parameters based on laboratory experiment is presented (Sect. 3). Finally, results of numerical simulations are presented in Sect. 4. The paper is concluded by results discussion and final remarks. "Appendix" contains listing of adopted shape functions of the first-and second-order perturbation solution.

Derivation of governing equations
We consider a rotating hub-beam structure which is carrying a tip mass m t as shown in Fig. 1. The highly elastic isotropic and slender beam with spanwise embedded piezo-layers is attached to the hub at an arbitrary preset angle θ (Figs. 1b and 2b). In order to describe motion of the structure, we introduce a fixed (X 0 , Y 0 , Z 0 ) coordinates system originated at the hub centre C, a set of coordinates (X, Y, Z ≡ Z 0 ) rotating with the hub and with the same origin, and a set of local coordinates (x, y, z) fixed to the beam and oriented along the symmetry axes of the undeformed blade with the origin at clamping point 0. We assume the rotation is performed about Z ≡ Z 0 axis and it is described by a temporary position angle ψ(T ). The longitudinal and transversal displacements u(x, T ) and v(x, T ) of an elementary segment attributed to point A, which is next moved to A position, are defined in the rotating coordinate frame as presented in Fig. 2a. Due to its high slenderness and low stiffness, the beam can be deformed with large amplitudes. The kinematics of the beam deformations at any arbitrary time instant is shown in Fig. 2b where dx denotes the initial length of an infinitesimal beam element, while ds is its extended length. Moreover, we assume that the beam is allowed to deform only in (x, y) plane because of its high stiffness in the orthogonal direction z. For the specific preset angle θ = 0, this plane would overlap with the plane of rotation. Moreover, in the presented analysis, twist of the beam is neglected leaving these more general cases for future developments. In its reference configuration, the beam has length L, cross section A, second moment of inertia I , and mass per unit length ρ 1 ; finally the hub radius is denoted R h and its mass moment of inertia J h .
The kinematic relations in the extensible-unshearable (Euler-Bernoulli) beam element are given as where (...) denote derivative with respect to spatial coordinate ∂... ∂ x , in the following notation(...) describes a time derivative ∂... ∂ T . Assuming the material to be uniform (the Young modulus E is constant), the axial force N along the deformed beam element is related to the axial stiffness E A, and the bending moment M to the mechanical curvature and bending stiffness E I as where elongation due to extension of the beam's midplaneê and mechanical curvature κ m are given asê Assuming the linear elastic behavior of the material and given above kinematic relations, the potential energy of the beam element is The kinetic energy of the beam-hub system has three components, which arise from motion of the hub T h , beam T b and tip mass T t : where The equations of motion are derived using the Hamilton's principle of the least action where D and W are works done by dissipative and nonconservative forces in the system, respectively Terms c u , c v are two independent linear viscous damping coefficients in longitudinal (x) and transverse (y) directions of the beam, respectively. Note that the ydirection damping (second component of Eq. (8) 1 ) contains terms responsible for beams oscillation about xz plane and slewing motion of the rotor. External excitation applied to the structure can be represented by any time-dependent force functions f v (x, T ) and f u (x, T ).
Variations of kinetic and potential energies, dissipation, and external works are given by Substituting Eqs. (1)-(6) and (8)-(9) into Eq. (7), then integrating by parts, we get a set of partial differential equations of motion and corresponding boundary conditions. Due to computational difficulties, we limit ourselves to the specific case of constant hub speeḋ ψ(T ) = const. Moreover, since the problem will be solved by the multiple timescales method up to the third order of approximation, the formulas (1) and (3) are expanded in Taylor series up to fourth order of geometric nonlinearities. After some mathematical manipulations, one arrives at the system of two nonlinear partial differential equations governing the longitudinal and transverse motion: as well as the associated boundary conditions The underlined terms Q up , Q uc , Q vp and Q vc present on the right-hand side of conditions (13) 1 and (13) 3 are artificiality introduced control generalized loads resulting from the action of piezoelectric actuators. It can be shown if the active layer is distributed along the whole span of hosting specimen (e.g., beam or plate), the action of the piezoelectric actuators is mathematically represented as non-homogenous boundary conditions at the free end of the structure. These nonzero terms represent the induced dynamic moment/shear force and are effective and robust computational ways toward the implementation of a feedback control in mathematical models of active structures. Since the proposal by Lagnese [21] this boundary control moment/force methodology has been successfully adopted by many investigators to study the behavior of the plates, shells, and beam structures with feedback control [6,15,22,31].
The subsequent terms Q up , Q vp , and Q uc , Q vc represent proportional (subscript 'p') and cubic (subscript 'c') functions of longitudinal and transverse displacements, respectively. Considering the time delay τ in the system, these are expressed by formulas

Solution by perturbation method
The analysis is applied for the pure nth flexural mode without internal resonance interactions. To solve the problem, we introduce three timescales where t 0 , t 1 , t 2 are fast and slow timescales, respectively. Parameter ε is a formal small parameter and serves as a book keeping device grouping small terms in a proper perturbation order. Thus using the chain rule the first-and second-order time derivatives arė Next, we seek the solutions in a series of small parameter at the subsequent three orders of perturbation: where v(x, T ) denotes the solution without time delay, and v τ (x, T, τ ) is its counterpart when time delay τ is considered. The variables u corresponding to the longitudinal direction are expanded in the similar way. For the sake of text brevity, they are deliberately omitted here.
To decouple linear response of the structure, we assume the angular speed to be of ε order (relatively small quantity) and writė By this assignment, the small book-keeping parameter ε shifts the angular velocity term from first to the second and higher orders of the problem. Furthermore, eliminating from the first-order solution will result in a classical cantilever beam equation (linear part of the solution). Without this assumption, static predeformations due to constant angular speed could be studied.
The external excitations f u , f v present on RHS in (10) and (11) as well as damping coefficients and control loads are shifted to the third and second perturbation order of the problem, respectively. Following this assumption, they are declared by where ξ v (ξ u ) describe the amplitudes of uniformly distributed external excitations at frequencies ω v (ω u ). The later ones are expressed by corresponding natural frequencies ω 0v (ω 0u ) at the fast timescale and detun-ing parameters σ v (σ u ) attributed to the second slow timescale t 2 .
Inserting equations (15)-(18) into expressions (10)-(13) and collecting ε power-like terms, we obtain the following sets of equations and corresponding boundary conditions • O(ε 1 ): boundary conditions at x = 0 and boundary conditions at x = 0 and

First-order solution
The solutions to the first-order problem Eqs. (19)- (22) are given as products of time dependent functions and linear mode shapes where are complex amplitudes dependent on both slow timescales and the over bar symbol denotes complex conjugate.
Parameters r 1 -r 4 are trigonometric coefficients which satisfy the boundary value problem (20)- (21). When considering the time delay the sought delayed solutions to the first-order problem u τ 0 , v τ 0 are again defined as products of time dependent amplitude functions and linear mode shapes. However, the former ones are shifted by appropriate time delays τ u and τ v whereas the deformation modes stay similar to their counterparts in regular solution as given by Eq. (28) This treatment of time delay effect within the multiple timescales method framework was successfully adopted by Rusinek et al. in [37].
The natural frequencies ω 0u , ω 0v in axial and transverse directions are given by respectively, while φ u and φ have to be found by solving the transcendental equations Subsequent roots determine the mode shapes Eq. (28) 3,4 and corresponding natural frequencies Eq. (30).

Second-order solution
Considering the second-order problem, we solve Eqs. (22)- (24). Namely, in the first attempt, we are focused on secular generating terms (the ones containing e ±iω 0u t 0 and e ±iω 0v t 0 ) to vanish. Following this, the solvability conditions are After integration by parts and substituting the boundary conditions Eqs. (23)-(24) and performing simple algebraic manipulations, we get Then analyzing Eqs. (22)-(24), we decompose harmonics, calculate time derivatives, and then reduce partial differential equations to ordinary differential ones. Finally, the second-order solutions are given by where r 5 -r 12 are real valued functions dependent on A, E, I, ρ 1 , ω 0u , ω 0v that represent mode shapes, while using parameters R h , and θ we can manipulate their magnitudes. To keep the length of the paper concise these functions are reported just graphically only for a set of fixed parameters as given in "Appendix A". Studying the Eq. (22) and associated boundary conditions (23)-(24) one observes the time delays to be absent at this order of perturbation. Therefore, the second-order solution of time-delayed problem is identically zero.

Third-order solution
In the third-order problem Eqs. (25)-(27), we consider only terms proportional to e ±iω 0v t 0 and e ±iω 0u t 0 , to which solvability conditions are applied (multiply functions by linear solutions and next integrate by parts in 0 to L limits). After cumbersome computations, taking into account boundary conditions (26)-(27), one obtains the following set of four modulation equations in the slow timescale t 2 : where q 1 -q 7 and p 1p 9 are real valued coefficients calculated during integration process. We preferentially convert the above equations expressing complex ampli-tudes in their polar form by introducing new variables and obtain a modified set of four modulation equations of amplitudes a, b and phase angles β a , β b Coefficients q 1 -q 7 and p 1p 9 contain terms of first-and second-order solutions and can be calculated numerically for fixed properties of the beam. Material data and geometry of the structure used to derive approximate solution were measured within laboratory experiment described in Sect. 3. Further parametric analysis of the system is presented in Sect. 4.
Finally, the approximate solutions of Eq. (16) take the form where (38)

Stability analysis
Stability analysis of analytical solutions is performed on the basis of modulation equations (37) which are rewritten in a shorter forṁ 1 (a, β a , b, β b ) , , β a , b, β b ) , , β a , b, β b ) , where "dot" means time derivative with respect to timescale t 2 . In the steady state the right-hand sides of Eq.
If any of the eigenvalue has a positive real part, the system becomes unstable.

Laboratory experiment
To identify the coefficients of the analytical model the laboratory experiment was conducted. The test stand comprised the laser vibrometer PSV-500 camera system-see Fig. 3 marks (a,b), and a highly flexible specimen (d) mounted to a dedicated grip fixed to an anti-vibrational slip table TIRA TGT MO 48XL (c). A series of reflective markers were sticked on the beam surface to be traced by the scanning laser. Moreover, piezoelectric macro fiber composite (MFC) patches M8528-P1 by Smart-Material Corp. were bonded on the tested blade as shown in Fig. 4. These are the elements operating in d 33 mode, so the piezoceramic poling direction (being in-plane of the transducer) when oriented along the beam span axis x excites the specimen bending (note also Fig. 1). To amplify the signal from the vibrometer generator, a dedicated MFC highvoltage amplifier was used.
The above data of the tested prototype have been adopted for the numerical investigations. In the further analysis, the governing equations and BCs are transformed to dimensionless forms. The space coordinates are normalized to the beam length:x = x/L, u = u/L,v = v/L,R h = R h /L and angles remain as they areθ = θ,ψ = ψ. Dimensionless time is defined ast = t/ω , where ω = E I ρ 1 L 4 , and thus corresponding dimensionless natural frequencies and dimensionless angular velocity are now given byω u = ω u /ω ,  The adopted transformation rules result in the dimensionless bending and axial stiffnesses normalized asĒ I = 1,Ē A = E AL 2 E I , dimensionless beam mass per unit lengthρ 1 = 1 and dimensionless tip massm t = m t ρ 1 L . Dimensionless values of the physical parameters given in Table 1 arē L = 1,Ē I = 1,Ē A = 5.24481 × 10 6 ,ρ 1 = 1, (42) and the tip mass is equal to zero for the studied case, m t = 0.
Coefficients q 1 -q 7 and p 1p 9 present in Eq. (35) were calculated for two basic cases:

Parametric studies of the system
The derived analytical solutions of the governing equations expanded up to ε 3 approximation order as obtained in Sect. 2.2 enable comprehensive parametric study of the system. To this aim, bifurcation diagrams are constructed to present the effects of various structural parameters on system response characteristics. In particular, parameters , σ v , ξ , R h , θ are kept in Eq. (37) and thus can be used in these analyses.
The modulation equations obtained from MTS method involve also control signals C u , C v . They can be modified according to the adopted control law to meet the required beam behavior. In this paper, we consider two control strategies, namely a proportional P control and cubic control C. Proportional strategy considers the input signal multiplied by a gain coefficient g p and supplied with some time delay τ to the actuators. In the case of nonlinear cubic C control, the signal is raised to power 3, gained by g c factor and than supplied with delay τ .
Moreover, by calculating coefficients p i (i = 1, . . . , 9) and q i (i = 1, . . . , 7) as present in (35) for individual longitudinal and transverse vibration modes the multiple different dynamic states of the system can be considered.
We start the analysis from solving the modulation equations analytically. To this aim, we introduce new variables γ like , and next substitute them into Eq. (37). In a steady state therefore, the modulation formulas simplify to the set of nonlinear but algebraic equations. Solution to this system involve response amplitude as a function of selected bifurcation parameters including the detuning parameters σ u and σ v . Moreover, these relations can be used to tune control parameters for effective vibration reduction. We note that for a general case all four modulation Eq. (37) are involved in the solution. However, if we neglect axial excitation and control in this direction, one can set a (t 2 ) = 0 followed by a(t 2 ) = 0 relation due to damping. Then, the dynamics of the structure is governed only by the third and the forth equations of the system. It should be emphasized this simplification can be done if we confine the transverse excitation frequency only to the vicinity of the corresponding natural frequency and exclude the possibility of internal resonance. Study of the structural dynamics involving the internal longitudinal-transversal resonances for simply supported beams can be found in [23]. This problem for rotating beams with a tip mass will be discussed in a separate contribution.

Natural vibrations
We start the analysis from analytical solutions to determine backbone curves which represent nonlinear natural vibrations around the first and second flexural frequencies. Dimensionless values of the linear natural frequencies corresponding to non-rotating system con- It can be observed for the adopted structural parameters the first longitudinal frequency ω 0u is located far away from the flexural ones. Thus, as it has been mentioned above, the internal resonance case does not occur in the given problem formulation.
To obtain the expressions for backbone curves, we substitute zero values on the excitation, damping, and control gains in Eq. (37). Then, we get the final expression for the first and second flexural modes, respectively. When analyzing the relation (45), one observes the presence of terms involving the rotor angular velocity in power 2 and 4. The plotted shape of the backbone curve corresponding to fundamental vibration mode and zero angular velocity (non-rotating structure) is presented in Fig. 6a. The curve reveals the evident hardening nature. If the angular speed is increased, the curve gets shifted right toward higher frequencies, and its hardening-like property is maintained (Fig. 6b). The linear natural frequency increases as a parabolic function of rotor angular speed which is presented in Fig. 6c.
In contrast to the first mode, the backbone curve for second flexural mode exhibits a softening (Fig. 7a) behavior for non-rotating structure. This nature of the curve is preserved also for rotating system-cases of = 2, = 3, = 5 presented in Fig. 7b. It is interesting to note that the inherent softening behavior is strong enough to prevail the stiffening effect due to centrifugal loadings.
The proposed analytical model of the structure and obtained approximate solutions show that both modes exhibit opposite characteristics. This fact is in an agreement with results published by Thomas [38] consider- ing relatively small rotations, but in contrast to Arvin et al. [2] where the introduced so-called equivalent nonlinearity coefficient depends on the rotor angular speed as well. Based on the analytical results, we may conclude the angular velocity does not affect neither hardening nor softening features of the backbone curves but just shifts them toward higher frequencies. This means the linear part of the stiffness matrix gets larger if rotation increases. But the nonlinear effect resulting from large amplitude oscillations maintains the same magnitude revealing stiffening/softening features depending on mode order. The first and the second modes having opposite nature will be studied for excited vibrations and structural control.

Forced vibrations
Forced vibrations of the rotating beam structure are studied assuming the excitation is imposed in the transverse direction only. Presuming steady-state vibrations and neglecting excitation in axial direction in Eq. (37), the last two formulas yield Eliminating phase angle γ b (t 2 ) and after mathematical transformations, we get where α 1 = 1 32 p 2 5 + g 2 vc p 2 6 + 2g vc p 5 p 6 cos τ v α 2 = 1 4 p 5 2 p 7 + p 8 R h + p 9 cos 2 θ − σ v +g vp p 3 (g vc p 6 + p 5 cos τ v ) +C v p 2 p 6 g vc sin τ v + p 6 g vc 2 p 7 + p 8 R h In the above equation, we can distinguish terms dependent on angular velocity , hub radius R h , external excitation ξ v , and related to linear g vp and nonlinear g vc control laws and time delay τ v . Equating to zero the control parameters, we get equations of the resonance curves around the first and the second resonance zones as a function of the angular velocity and the preset angle θ . For the following analysis, we assume that the preset angle is equal zero (θ = 0). This correspond for the beam oscillating in the rotor plane (see Fig. 5). Modal damping is approximated on the base of experimental tests as 2% of first natural frequency value: C v = 0.07032 and 1% for the second mode C v = 0.22035. Amplitude of excitation ξ v is varied as reported in figure captions.
Resonance curves around the first flexural frequency are presented in Fig. 8. Despite that large oscillations in the system just a very faint stiffening effect is observed for this mode (Fig. 8a), much too small for the jump phenomenon to occur. This feature is observed for zero angular velocity and also for rotating structure-cases = 2, = 3 and = 5 corresponding to blue, red and green curves, respectively (Fig. 8b). Again the results are similar to the backbone curves, the so-called stiffening effect due to rotation leads solely to the right shift of the curves without a change of the stiffening effect caused by the nonlinear beam oscillations.
For the second resonance zone, nonlinear softening is rather strong. In Fig. 9a, computed for = 0, the softening behavior is observed for amplitudes below of 8% of the beam length. As it is demonstrated in Fig. 9b the increased angular speed (cases = 2, = 3, = 5) does not change the intensity of softening phenomenon but just shifts the resonance zone. Stability of frequency response curves is studied by computing eigenvalues of Jacobian matrix (41); dashed line represents unstable part of the solution and solid line defines stable solution in the following Figures.
The results presented in this section shall be used in the next steps when the proportional and cubic control strategies are engaged in order to reduce vibrations around the discussed resonance zones. The outcomes obtained in this section are consistent with authors previous work in [42].

Controlled forced vibrations
The derived analytical model of the active blade and approximate solutions to governing equations enable analysis of the system under assumed control law. This involves additional terms present in the boundary conditions Eq. (35) and later modulation equations-see terms with g up , g uc , and g vp , g vc in formulas (37). They represent loads generated by control unit according to the boundary control method approach adopted in this research. In a real set-up structural control can be achieved by the controlled voltage supplied to active elements. The mathematical model of the combined blade-control unit structure considers also a time delay τ v which may occur due to inherent properties of actual devices or can be added intensionally to enhance system's performance. It was decided to study two control strategies, namely the linear control (P) and cubic control (C). In the linear control method, the input signal is multiplied by a gain coefficient g p and supplied to the actuators. In the case of nonlinear C control the signal is raised to power 3 and gained by g c factor.
Remembering that dynamics of the beam close to its first natural frequency is almost linear, we apply P control with time delay. Based on the Eq. (47), substituting g vc = 0 we get Equation (48) depends on controller gain, time delay τ v , angular velocity , and amplitude and frequency of excitation, ξ v and σ v , respectively. Therefore, we can evaluate the effectiveness of the controller at various combinations of excitation conditions. The response of the system for fixed amplitude of excitation ξ v = 0.1 and activated linear controller is shown in Fig. 10. The figures present an overview of the response for fixed gain value g vp = 0.2 with respect to detuning frequency σ v and time delay τ v varied in −2π , 2π interval. The computed 3D surface shows periodically repeated peaks where the applied control does not work properly. However, there are also valleys corresponding to small amplitudes and effective vibrations suppression. Comparing graphs for = 0 (Fig. 10a) and = 2 (Fig. 10b), we observe the controller works in a very similar way. The only difference is the shift of the response surface toward higher frequencies. This is visible by localization of peaks along σ v axis (note that scales for σ v are different on both plots). Studying these plots, one may also conclude about the influence of the time delay. It is evident the proper tuning of the τ v parameter with respect to time delay the effective suppression of vibrations may be achieved.
The selected cross sections of these 3D surfaces are presented in Fig. 11a and b for the stationary and rotating configuration, respectively. In Fig. 11a the resonance curve of the system without control is sketched by the grey line. If the P controller is engaged and the input signal is supplied with zero time delay then the curve is just shifted toward lower frequencies without any amplitude change (black line). However, if the time delay is introduced, then we may simultaneously shift the resonance curve and reduce vibrations. As it is seen in Fig. 11a the amplitude can be decreased even over three times for τ v = 2.0 (orange curve).
In the case of the rotating beam (Fig. 11b), the resonance curve without control due to rotation is shifted toward higher frequencies (magenta line with respect to grey curve as copied from (a)). Then if the P controller is operational without time delay (τ v = 0) the curve is shifted back toward lower frequencies (black continuous line). Similar to the stationary configuration case, the proper selection of time delay τ v may significantly reduce vibration amplitudes-see, e.g., orange curve corresponding to τ v = 2.0. However, the optimal delays are different than for non-rotating system. It should be also noted the response amplitude can be not only reduced, but also magnified to large values if the time delay is tuned improperly. Meanwhile, the characteristics are shifted toward higher or lower frequencies as it is presented in Fig. 11c for τ v = 3 (blue) and τ v = 6.5 (red) in Fig. 11d, respectively.
Direct comparison of P control effectiveness around the first resonance zone for fixed excitation amplitude and frequency is presented in Fig. 12. The detuning of excitation frequency σ v is selected close to the first resonance and time delay is varied within −2π, 2π interval. As reported above, the proper selection of time delay results in a large reduction of vibration amplitudes both for stationary (Fig. 12a) and rotating beam (Fig. 12b). Furthermore, one observes the zones of effective vibrations suppression are significantly wider for the rotating beam case.
As the second option, we propose to apply the nonlinear cubic C controller. Thus the P control is switched off by setting g vp = 0. Substituting the condition to equations (47) we get Solutions of Eq. (49) enable to determine impact of C control on the system response. The idea of this controller is to influence mainly the fundamental property of the structure characteristics that is represented by the slope of the backbone curve. In fact, setting time delay τ v = 0 and varying gain g vc we may modify the properties of the system from hardening to softening (Fig. 13a) behavior. For the gain g vc = 0.12, we may obtain linear-like response (red curve in Fig. 13a).
Results of the rotating system analysis are presented in Fig. 13b. The introduced rotation does not affect the controller performance apart from the fact that the resonance curves are shifted and their slopes are modified around the shifted linear natural frequency value.
The effectiveness of the C controller as well as combined P − C control we demonstrate for the second natural frequency. Due to its strong nonlinear softening effect, the comparison of nonlinear versus linear control strategy is more appealing.
The analysis of the second-mode dynamics is started from P control governed by Eq. Since vibrations amplitudes around this resonance zone are much smaller then for the first resonance the gain g vp need to be adjusted to higher values. Based on experimental data, we accept g vp = 1 as a reference value. The second-mode response surface plots are prepared keeping the excitation amplitude fixed ξ v = 0.4 and varying detuning parameter σ v (excitation frequency) and time delay τ v . The amplitudes of oscillations for the stationary and the rotating beam configuration = 5 are presented in Fig. 14a and b, respectively.
The darker blue surfaces represent the desired controller performance with just a few hills corresponding to slightly increased amplitudes. Studying the plots it may be concluded for the given excitation frequency one may adjust the time delay to maximize the effec- Fig. 11 Resonance curves of the system with P control with various time delay; a for = 0, τ v = 0 (black), τ v = 0.5 (blue), τ v = 1.0 (green), τ v = 1.5 (pink), tiveness of the controller and retain possibly lowest vibration amplitudes. The summits of pale blue hills are located right at the second natural frequency (i.e., for σ v = 0) in the case of the non-rotating beam or shifted toward higher frequencies to about σ v = 4 for the = 5 rotating beam. This offset results purely from the centrifugal stiffening effect. Nevertheless, the inherent strong nonlinear behavior (softening type see Fig. 9) of the second mode brings additional large amplitudes solutions represented by orange surfaces in Fig. 14 even if P controller is on.
The individual response curves for selected time delays are clearly presented in Fig. 15. The resonance curve of the system without control is marked by grey line. When the control is triggered, the overall performance of the system is quite similar to the first mode. Increase in time delay first shifts the curves slightly toward lower frequencies (Fig. 15a), but next the trend is reversed and peaks are shifted a bit to the right (Fig. 15b). By proper tuning of the time delay parameter τ v , we can decrease amplitudes of the response while keeping the natural frequency almost unchanged. Meanwhile, the softening type of the characteristics is maintained. On the other hand the poor choice of time delay causes the dramatic increase in vibrations amplitudes and enhances nonlinearity of the system as presented in Fig. 15a for τ v = 0.5 and Fig. 15b for τ v = 3.0. These cases correspond to orange areas in Fig. 14. Possible consequences of the poorly tuned P control are presented also in Fig. 16. Adjusting time delay, we may get low amplitude oscillations but meanwhile above these low amplitude sections there are isolated closed loops representing large amplitude oscillations with unstable bottom part of the "isola." For specific time delays where the "isolas" exist the basins of attractions of all possible solutions have to be carefully checked.
To eliminate the presented above drawbacks of Pcontrol we propose to apply C strategy which-as it has been demonstrated earlier-allows to modify the slope of the resonance curve. The coefficients of the govern- First it should be noted the value of the cubic control gain g vc may get even two orders higher values than its counterpart applied for the first mode. This is related to the fact that the second mode dimensionless amplitudes are much smaller than the first mode ones. Since these values are put to the power 3 the difference must be "compensated" by higher magnitudes of gain at second mode. Accordingly, the second mode cubic gain may approach values of about 1500 but we decided to take a smaller value setting g vc = ± 500 as a reference level.
In contrast to the first-mode analysis now, we show that the C-control strategy applied to mode 2 can modify the slope of the response characteristics by tuning the time delay. The original resonance curve (grey line in Fig. 17a) can be modified in terms of amplitudes and also its slope can be changed from softening to hardening. Moreover, the further increase in time delay results in reduced amplitudes combined with amplified softening behavior evidenced by increased slope (Fig. 17b). However, for negative time delays, the reduction of response is not so evident. Furthermore, the double nonlinearities, namely the structural and control ones result in presence of additional solutions observed by isolated curves (Fig. 17c) of high amplitude. The spotted additional solutions are adverse events in terms of vibration reduction but they may be interesting while designing the energy harvesters.
Having above results now we apply mixed P − C control for the nonlinear second mode resonance. To benefit from advantages of both controllers, we apply cubic controller to get a pseudo-linear system and then the linear controller to reduce vibrations. Therefore, we set g vc = −500 and τ v = 1 for the analysis which correspond to linear-like response in Fig. 17 (red curve). Next, we are to adjust the settings of the additional P controller. The response surface in Fig. 18 presents the characteristics of the beam against detuning parameter σ v and gain g vp . Obviously, for g vp = 0 we get purely  cubic control strategy. To get more reduction, we have to avoid the peak and select value out of this zone. The resonance curve with activated P −C controller at time delay τ v = 1 is presented in Fig. 19a by the black line. The vibrations amplitude is reduced and the only one solution exists. For comparison we study control setting τ v = −1.5 (red line) where the greater reduction of the amplitude can be achieved. However, in this case, additional isolated solutions exist (red isola), and they are of large amplitude. Such situation requires a careful detection of their basins of attraction to guarantee safe system dynamics. The parameters of the combined P − C controller can be applied for rotating beam as well. The only difference is that the resonance curve is shifted toward higher frequencies, as presented in Fig. 19b.

Conclusions
The presented mathematical model of the rotating cantilever beam structure with active elements enables analysis of free and forced oscillations of the structure corresponding to moderate-large amplitudes. The partial differential equations and associated boundary conditions are derived considering longitudinal and transversal vibrations of the beam rotating with fixed angular velocity and preset to the hub at any arbitrary angle. The performance of active elements has been represented by non-homogenous terms in BCs. These involved linear and nonlinear (cubic) control signals with time delay. The set of partial differential equations and associated BCs have been solved directly by the multiple timescales method up to the third-order perturbation. The physical parameters of the model including electromechanical parameters of the active elements have been determined within the laboratory tests.
The analytically derived modulation equations of the beam dynamics include most important structural parameters and control gains that may be considered as bifurcation parameters. The performed analysis of the solutions revealed hardening for the first mode and softening of the second vibration bending mode. These effects are maintained also if angular velocity is increased. It has been shown the rotation does not change the curve slope (hardening or softening) but just shifts the resonance zones toward higher frequencies which is in an agreement with results presented in paper [38] if angular velocity is not high.
To reduce vibrations of the system the linear proportional P controller, nonlinear cubic C controller and finally mixed P − C controller is proposed. It has been shown for the weakly nonlinear system around the first resonance the P controller is effective both for the stationary as well as rotating beam. For the moderate or strongly nonlinear system behavior corresponding to second natural frequency multiple solutions occur and thus more advanced control techniques are applied. The suggested strategy is based on the cubic or mixed linear-cubic control results. It has been shown these two approaches are capable to significantly reduce nonlinearities and to eliminate multiple solutions. By a proper tuning of control gains and time delay we can determine analytically domains of system parameters corresponding to unique (single) solution and zones with safe controller performance. Therefore, it is possible to adjust controller parameters for safe operation of the system at large oscillations and various angular speeds of the beam.
The additional nonlinearity coming from controller combined with beam's second mode structural nonlinearity may lead to additional isolated solutions corresponding to large amplitudes. That specific operation of the system can be exploited in terms of energy harvesting when nonlinear oscillations are much desired effect. The obtained analytical results will be tested in the laboratory for actual implementations.