Nonlinear vibrations produced by unbalanced motors

Nonlinear Vibrations Produced by Unbalanced Motors Doctoral Thesis Javier González Carbajal Supervised by Prof. Jaime Domínguez and Prof. Daniel García Department of Mechanical and Manufacturing Engineering Faculty of Engineering University of Seville March 2017 i Contents 1 Introduction 1 1.1 State of the Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 Motivation and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Organization of the Document . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Mathematical Methods 13 2.1 First Order Averaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Second Order Averaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Singular Perturbation Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4 Hopf Bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5 The Poincaré-Béndixson Theorem . . . . . . . . . . . . . . . . . . . . . 25 3 The Case of Large Slope of the Motor Characteristic: Analytical Approach 27 3.1 Problem Statement and Assumptions . . . . . . . . . . . . . . . . . . . . 28 3.2 Alternative First Order Averaging . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Perturbation Approach: Derivation of the Reduced System . . . 36 3.4 Analysis of the Reduced System . . . . . . . . . . . . . . . . . . . . . . . 45 3.5 Classification of the Hopf bifurcations . . . . . . . . . . . . . . . . . . 55 ii 3.6 Conditions under which all System Trajectories are Attracted towards a Limit Cycle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4 The Case of Large Slope of the Motor Characteristic: Numerical Simulations 71 4.1 Global Bifurcations of the Limit Cycles . . . . . . . . . . . . . . . . . 72 4.2 Numerical Validation of Analytical Results . . . . . . . . . . . . . . . 80 5 The Case of Small Slope of the Motor Characteristic: Analytical Approach 89 5.1 Outer Region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.2 Inner Region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 6 The Case of Small Slope of the Motor Characteristic: Numerical Simulations 113 7 Torque-Speed Curves for the Whole Frequency Range 127 7.1 Computation of the Torque-Speed Curves . . . . . . . . . . . . . . . 128 7.2 A Global Perspective for the Cases of Large and Small Slope 136 8 Modelling and Simulation of the Vibrocompaction Process 141 8.1 Some Notes on the Real Process . . . . . . . . . . . . . . . . . . . . . . 142 8.2 Full Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 8.3 Simplified Model and Torque-Speed Curves . . . . . . . . . . . . . 156 8.4 Analytical Investigation of a Quasistatic Vibrocompaction Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 8.5 Numerical Results and Discussion . . . . . . . . . . . . . . . . . . . . . 171 9 Summary and Conclusions 203 Appendix 211 iii Acknowledgements Tras estos cuatro años de tesis, estoy en deuda con muchas personas y me gustaría poder incluirlas a todas en este agradecimiento. Confío en que sepan disculparme las que, por falta de espacio o por mi despiste existencial, no aparezcan expresamente nombradas. En primer lugar, quiero agradecer a mis directores, Jaime y Dani, todos sus buenos consejos y el clima de libertad que han propiciado desde el principio de este trabajo. Estoy convencido de que no podría haber tenido mejores maestros. Gracias también al profesor Emilio Freire, por orientarme sabiamente cuando me adentraba en terrenos matemáticos pantanosos. I would like to express my gratitude to Professor Gaëtan Kerschen for my stay at the University of Liège. Working with his group was a wonderful research experience, in an environment which could not have been friendlier. Thanks to Peppe, who kindly contributed to correct typos in the thesis, Lionel, Thibaut, Chiara, Edouard, Vincent, Iván, Mirco and many more good friends that I hope to keep for many years and made my stay in Liège warm and enjoyable. Thanks to Professor Jan A. Sanders, who was very willing to advise me about the attraction properties of averaged systems. His keen e-mail responses were truly helpful for the development of Chapter 5 in the thesis. Gracias a Joselu, compañero de fatigas desde los míticos tiempos en que inventamos la rosca cuadrada, a Aceituno, Jorge Julio, Sergio, Merche, Juan, David y al resto de compañeros del departamento, porque ha sido una alegría poder trabajar rodeado de buenos amigos. Extiendo con gusto este agradecimiento al gran Guido, pese a su obsesión posesiva con el espacio de trabajo y con las papas fritas. A Charly, Merchante, Diego, Ortega, Bill, Angelito, Andrés, Miguel, Yoshua, Koctel, Nacho, Derey, Valera y muchos más, por los buenos ratos y cervezas que una tesis necesita, aunque sigan pensando que la mía trata de tornillos. iv A Marta, mi equipo alpha particular, que ha sabido animarme como nadie y hacerlo todo muy llevadero estos últimos meses. Muchas gracias a mis padres, que son el mejor ejemplo que tengo, y a Pablo y Ana, por su apoyo y su cariño diarios. También al resto de mi gran familia, desde mis abuelas y abuelos hasta Julia y Rafilla. Este trabajo habría sido mucho más duro si no hubiera estado tan bien rodeado. Y principalmente gracias a Dios, a quien le debo todo. 1 INTRODUCTION The linear theory of vibrations is one of the most useful tools in the handbook of a mechanical engineer. It furnishes a solid and well-established mathematical framework for concepts such us resonance and linear modes of vibration, which are central for the modal analysis of structures and mechanical systems. A most appealing quality of this linear theory is its relative simplicity. Linear models exhibit special attributes which make them particularly useful from a practical point of view, such us the principle of superposition –whereby the response to a linear combination of excitations can be obtained as the linear combination of individual responses– or the fact that several attractors never coexist –the stationary motion of a forced, damped, linear system, is independent of the initial conditions–. It is thanks to these properties, together with the fact that many real structures are well represented by linear models in their ranges of operation, 2 1 Introduction that linear vibration analysis has become a chief tool in the design of mechanical systems and structures in industry. However, linearity is clearly an idealization. Nature is not linear in general, although it can be approximately represented by linear models in some cases. For instance, the linear theory of elasticity is known to give good results for structures undergoing small displacements and small strains. However, as soon as displacements become significant, nonlinear effects need to be considered (Luongo, Rega, & Vestroni, 1986). Some other possible sources of nonlinearity in mechanical systems are (Thomsen, 2003) - Material nonlinearities due to a nonlinear relation between stresses and strains in some materials. - Nonlinear body forces, such as magnetic or aerodynamic interactions. - Nonlinearities due to the physical configuration, such as those associated to discontinuous couplings, clearances or stops. Nonlinearity plays a key role in the dynamic behaviour of numerous real-world applications. Some examples are the motion of large wind turbines, the crashworthiness of vehicles or the vibrocompaction of granular materials. The aeronautical industry, with great interest in minimizing weight, is increasingly producing very light and slender structures, with the subsequent activation of geometrical nonlinearities. The use of materials such as carbon-fiber composites in aerospace applications can also produce a significant deviation from linearity. It is evident that, in all these situations, a linearized model would not be able to capture the real system dynamics. Thus, in order to have reliable predictions, nonlinearity would need to be included in the model, which in turn implies entering the complex field of nonlinear dynamics. Nonlinear dynamics is the branch of mathematics which intends to uncover the temporal evolution of nonlinear dynamical systems. Unlike its linear counterpart, nonlinear dynamics is not a closed subject –not even a mature one–. The main reason is that, while all linear systems are essentially the same, each nonlinear 1 Introduction 3 system is nonlinear in its own way. This means that, in general, conclusions about the behaviour of a particular nonlinear system cannot be generalized to any other. Besides, the dynamics of nonlinear systems is extremely rich, exhibiting a wide range of phenomena which cannot occur in linear systems, such as multistability, chaos or limit cycles. All this complexity renders it extremely hard to obtain analytical solutions to the nonlinear differential equations governing the system dynamics. In general, it is necessary to resort to numerical computations in order to get some insight into the system behaviour. From the above considerations, it is clear that a significant research effort needs to be oriented to a better comprehension of the dynamic behaviour of nonlinear mechanical systems. This will hopefully lead to efficient predictions about the performance of systems where nonlinearity has a significant effect, and will also motivate the design and development of new nonlinear components which are able to outperform their linear counterparts. It is the aim of this thesis to contribute to this general objective, by analysing in detail a particular class of nonlinear systems, namely those excited by unbalanced motors. The motion of unbalanced rotors constitutes one of the most common vibration sources in mechanical engineering (Boyaci, Lu, & Schweizer, 2015; Yang et al., 2016). Vibrations due to unbalance may occur in any kind of rotating systems, such as turbines, flywheels, blowers or fans (Shabana, 1996). Actually, in practice, rotors can never be completely balanced because of manufacturing errors such as porosity in casting, non-uniform density of the material, manufacturing tolerances, etc. (Xu & Marangoni, 1994). Even a subsequent balancing process will never be perfect due to the tolerances of the balancing machines. Moreover, some amount of unbalance generally appears during the operation of the machine, as a consequence of uneven wear, corrosion or unequal build-up of deposits (dirt, lime, etc.). Usually, rotor unbalance has a harmful effect on rotating machinery, since vibration may damage critical parts of the machine, such as bearings, seals, gears and couplings (Xu & Marangoni, 1994). However, there are applications where 10 1 Introduction fixed points of both averaged systems are obtained and their stability is analysed. Admittedly, some of the developments of this Chapter are not totally new, but a reformulation of already published treatments of the problem. However, the stability of the stationary motions of the system near resonance was not totally solved hitherto. By using attraction arguments, a detailed stability analysis of these solutions is conducted. Chapter 6 presents numerical results which validate the results and conclusions of Chapter 5. In Chapter 7, an alternative approximate method is used to obtain the stationary solutions of the system studied in previous chapters. This approach is based on the Method of Direct Separation of Motions, proposed by Blekhman (Blekhman, 2000) and has the advantage of providing a graphical representation of the stationary motions which is highly convenient with a view to the vibrocompaction analysis of the next chapter. Furthermore, this alternative procedure provides a very clear graphical comparison between the cases of large and small slope considered in preceding chapters. A new model is introduced in Chapter 8, more complex than the one analysed in previous chapters. This system is a first attempt to model the vibrocompaction process, which has not been done before, to the author’s knowledge. In addition to the nonlinearity intrinsically associated to nonideality, the model includes contact and impacts between the mixture and the platform supporting the unbalanced motor, and also between the mixture and the mould where it is contained. Furthermore, a nonlinear constitutive law for the mixture, which allows modelling the compaction itself, is proposed. It is shown that, under some conditions, this model can be transformed into the simpler system analysed in Chapters 3-7, which makes the preceding analytical developments useful for the vibrocompaction analysis. Several numerical simulations illustrate how the proposed model can be used to investigate the effect of different parameters on the final level of compaction achieved. 1.3 Organization of the Document 11 Finally, Chapter 9 summarises the conclusions of this work and proposes possible further investigations. 2 MATHEMATICAL METHODS Before the analysis of the problem under study, a brief description of the mathematical methods used within the thesis is presented. 2.1 First Order Averaging Perturbation methods constitute a broad class of mathematical techniques aimed at finding approximate solutions to a problem, based on the solution of a simpler related problem. In particular, averaging procedures are among the most widely used perturbation methods, presenting two relevant strengths (Sanders et al., 2007): - They are supported by rigorously proved theorems. - They can be systematically extended to any order of accuracy. In order to understand the basic idea behind averaging, consider a system of the form 󰇗 ( ), where is a small parameter and is -periodic in . The 14 2 Mathematical Methods dynamics of such systems contains two different time scales: a fast scale, associated to the fact that depends on , and a slow scale, associated to the fact that is a slow variable ( 󰇗 ( )). Then, it can be shown that the essential features of the system are maintained when it is replaced by its corresponding averaged system 󰇗 ∫ ( ) ( ). The idea is to take into account the mean effect of the fast oscillatory dynamics through averaging, so that the second system retains the long term behaviour of the first one. This transformation is useful because the averaged system is autonomous and, therefore, considerably easier to analyse than the non-autonomous original system. Although the above description has been given for averaging over time, it is sometimes useful to average over fast rotating angles, which play essentially the same role as time. Once the intuitive idea has been explained, some more rigorous results are now given, which will be used throughout the thesis. Then, consider an initial value problem of the form { 󰇗 ( ) 󰇗 ( ) ( )} ( ) ( ) (2.1) where is a vector of slow variables, is a fast rotating phase and is a sufficiently small, positive, dimensionless parameter, . Note that represents the circumference. Hence to say that merely means that , with functions and being -periodic in . Let be the solution of 󰇗 ( ) ( ) (2.2) where 2.1 First Order Averaging 15 ( ) ∫ ( ) (2.3) Then, if ( ) remains in , ( ) ( ) ( ) ( ⁄) (2.4) The proof can be found in (Sanders et al., 2007). This is the standard version of the theorem for first order averaging over angles. Note that, thanks to this theorem, an asymptotic approximation to the solution ( ) can be obtained by replacing the original system (2.1) with the approximate system (2.2). In other words, it is possible to reduce the system dimension, from to , upon averaging over the fast rotating phase. Note also that, as explicitly stated in (2.4), the approximation obtained through averaging is only valid during a limited time scale. This is a general feature of perturbation approaches, which can be overcome in some situations (e.g. when attraction exists, as described later on in this section). The limitation in the time scale where the approximation is valid will prove to be crucial in Section 5.2. Generalization to multi-frequency systems The generalization of this theorem to multi-frequency systems is straightforward (Sanders et al., 2007). The only caveat that needs to be taken into account is that function must be written as a sum of functions, each of one depending on only one of the fast angles. Thus, consider the multi-frequency system given by { 󰇗 ( ) 󰇗 ( ) ( )} ( ) ( ) (2.5) where represents the -torus. Then, assuming that function can be written as 16 2 Mathematical Methods ( ) ∑ ( ) (2.6) the first order averaged system can be obtained as 󰇗 ∑  ( ) ( ) (2.7) with  defined as in (2.3). Then, the error estimate is given by ( ) ( ) ( ) ( ⁄) (2.8) This multi-frequency version of the theorem will be used in Section 5.1. The concept of ‘Resonance Manifold’ A crucial point of this averaging process is that, for the approximation to be valid, each frequency ( ) must be bounded away from zero. The reason is that, when any ( ) approaches zero, a vanishing denominator appears in the higher order terms of eqs. (2.4) and (2.8). (This will be seen clearly in equation (2.12) of Section 2.2.) The sets of for which one or more functions ( ) vanish are called ‘resonance manifolds’. The failure of averaging in the vicinity of resonances can be easily understood by noticing that, near a resonance manifold, one or more angles are not fast and, consequently, we cannot average over them. The main consequence is the following: away from resonances, it is possible to average over all angles. However, in the vicinity of a resonance manifold, we can only average over particular angles, namely those whose frequencies remain bounded away from zero. This distinction will be necessary in Chapter 5. 2.1 First Order Averaging 17 Averaging with Attraction Asymptotic approximations obtained by averaging are, in general, valid on a time scale ⁄. However, this time scale can be enlarged if there exists attraction in the averaged system (Sanders et al., 2007). Consider again a system of the form (2.1), with averaged system (2.2). Suppose is an asymptotically stable fixed point of system (2.2), with domain of attraction . Then, for all , it can be shown that ( ) ( ) ( ) , ) (2.9) Thus, the approximation is uniformly valid –i.e. valid for all time– for solutions of the averaged system which are attracted by an asymptotically stable fixed point. Two different proofs of this theorem have been given by Sánchez-Palencia (Sanchez-Palencia, 1975) and Eckhaus (Wiktor Eckhaus, 1975). A similar result holds for trajectories of the averaged system which are attracted by an asymptotically stable limit cycle. In this case, the approximate solution is valid on , ) for all variables except the angular one, i.e. the variable which measures the flow on the closed orbit (Sanders et al., 2007). This is equivalent to say that the closeness to the limit cycle can be uniformly approximated, yet not the position on it. The reason is that any small deviation on the frequency is accumulated over the cycles, giving rise to large errors after a sufficient number of periods. 2.2 Second Order Averaging In some situations, more accurate approximations than those provided by a first order averaging (Section 2.1) are required, as will be the case in Chapter 5 of this 18 2 Mathematical Methods thesis. Fortunately, there exist general results for order averaging which allow increasing the precision of the approximation as much as needed. Here we only show some results for order 2, which is enough for the purpose of this thesis. Hence consider a system of the form 8 󰇗 ( ) ( ) , -( ) 󰇗 ( )9 ( ) ( ) (2.10) where, following the notation in (Sanders et al., 2007), brackets are used in , - to stress that this term is a remainder of an expansion in powers of . This is also noted by the fact that , - depends on , while and do not. Let * ( ) ( )+ be the solution of 8 󰇗  ( )  ( ) 󰇗 ( )9 ( ) ( ) ( ) (2.11) where  is defined as in (2.3) and ( ) ∫[ ( )  ( )] (2.12)  ( ) ∫[ ( ) ( ) ( )] (2.13) Constant in (2.12) is chosen in such a way that ∫ ( ) (2.14) 2.2 Second Order Averaging 19 while symbol in (2.13) represents differentiation with respect to . Then, the error estimate is ( ) ( ) ( ( ) ( )) ( ) ( ⁄) (2.15) This result will be used in Chapter 5. The proof can be found in (Sanders et al., 2007). Note the presence of in the denominator of expression (2.12), in connection with the concept of ‘resonance manifold’ explained in Section 2.1. A vanishing frequency would make unbounded and, according to (2.15), the approximation would not be valid anymore. In other words, the averaging transformation becomes singular when approaches zero. 2.3 Singular Perturbation Theory The Singular Perturbation Theory (SPT) explains the behaviour of a particular class of fast-slow systems (Hunter, 2004; Lesne, 2006; Verhulst & Bakri, 2006). As has been done with averaging, an intuitive description of the theory is first given, followed by the exposition of some more rigorous results which will be used later on in the thesis. Systems where the SPT is applicable are those exhibiting a singular limit. This means that the system depends on a parameter in such a way that, when the parameter approaches some limiting value, the general solution of the problem is qualitatively different to the solution of the limiting problem. Typically, systems of this type exhibit two separate time scales –or space scales– with very different behaviours. A well-known example is the boundary layer theory in fluid mechanics, where the effects of viscosity are relevant in a very thin layer of fluid, close to a bounding surface, while being negligible for the rest of the domain. 26 2 Mathematical Methods is oriented inwards, as depicted in Fig. 2.3. Then, all trajectories which start inside are restricted to remain inside. Fig. 2.3 Phase portrait showing a closed, bounded region with the flow directed inwards at the boundaries. If contains no fixed points, then it contains at least one stable limit cycle. R 3 THE CASE OF LARGE SLOPE OF THE MOTOR CHARACTERISTIC: ANALYTICAL APPROACH This Chapter investigates the dynamics of a 2-DOF system consisting in an unbalanced motor attached to the fixed frame by a nonlinear spring and a linear damper. As commented in the introduction, two different scenarios need to be considered separately, depending on the order of magnitude of the slope of the motor characteristic. The Case of a Large Slope is considered in this Chapter and the following, while the alternative situation will be studied in Chapters 5 and 6. 28 3 The Case of Large Slope: Analytical Approach 3.1 Problem Statement and Assumptions Consider the system depicted in Fig. 3.1. It consists in an unbalanced motor attached to a fixed frame by a nonlinear spring –whose force has linear and cubic components– and a linear damper. The cubic component of the spring gives the possibility to model a nonlinear behavior for the structure supporting the motor (Mettler, 1962). The effect of gravity can be shown to have no relevance (Dimentberg et al., 1997) and, therefore, it will not be included in the model. Fig. 3.1 Model Variable stands for the linear motion, is the angle of the rotor, is the unbalanced mass with eccentricity , is the rest of the vibrating mass, is the rotor inertia (without including the unbalance), is the viscous damping coefficient and and are, respectively, the linear and cubic coefficients of the spring. The equations of motion for the coupled 2-DOF system are (El-Badawy, 2007) 󰇘 󰇗 ( 󰇗 󰇘 ) 󰇘 ( 󰇗) 󰇘 (3.1) where , and an overdot represents differentiation with respect to time, . Function ( 󰇗) is the driving torque produced by the motor –given by its torquespeed curve, also known as static characteristic– minus the losses torque due to 0 I 0 m 1 m r b x  , k  3.1 Problem Statement and Assumptions 29 friction at the bearings, windage, etc. We assume this net torque to be a linear function of the rotor speed: ( 󰇗) 󰇗 (3.2) Although ( 󰇗) includes the damping of rotational motion, we will usually refer to it shortly as ‘the motor characteristic’. As will be seen later, it is convenient for the purpose of this chapter to write the driving torque in an alternative way. Then, denoting by the linear natural frequency of the oscillator, given by √ ⁄, the motor torque can be written as ( 󰇗) ( 󰇗 ) (3.3) where represents the driving torque at resonance ( ( ) ). From equations (3.2) and (3.3), the relation between constants and can be directly deduced: (3.4) Along the whole thesis, the motor characteristic will be written as (3.2) or (3.3), depending on the situation. It should be kept in mind that these two expressions are totally equivalent. The important point is that the driving torque is assumed to follow a linear relation with the rotor speed. It is further assumed that –the driving torque decreases with the rotor speed–, as is usual for most kinds of motor. This assumption will prove to be of major importance. 30 3 The Case of Large Slope: Analytical Approach Fig. 3.2 Typical static characteristic for an asynchronous motor Fig. 3.3 Static characteristic corresponding to equation (3.3) As an example, the static characteristic of an induction motor is depicted in Fig. 3.2. Note that such a motor is usually designed to work on the region 󰇗 , where the curve could be reasonably approximated by a straight line with negative slope. The simplified motor characteristic given at (3.3) is represented in Fig. 3.3. m L  m L n  C 1 D peak   3.1 Problem Statement and Assumptions 31 In the second of equations (3.1), which imposes the equilibrium of the rotor, the last term is of great significance, since it accounts for the torque on the rotor caused by linear motion of the system. Its physical interpretation can be readily understood with the aid of Fig. 3.4. Due to displacement ( ), a horizontal inertial force acts on the unbalanced mass and generates a torque with respect to the rotor axis. This particular term of the equations of motion is what makes the excitation nonideal, for it takes into account how vibration influences rotation. If this torque due to vibration did not exist –or if it was negligible–, the rotor equilibrium equation would reduce to 󰇘 ( 󰇗), and it could be solved for ( ) regardless of the linear motion. Then, this solution ( ) could be introduced in the first of equations (3.1) as a prescribed excitation. Fig. 3.4 Torque on the rotor due to vibration By defining ⁄ ⁄ √ ( ) (3.5) the equations of motion can be written in a more convenient dimensionless form 󰇘 󰇗 ( 󰇗 󰇘 ) 󰇘 ( 󰇗 ) 󰇘 (3.6) 1 mx sin r   r x 1 m 32 3 The Case of Large Slope: Analytical Approach where a dot now represents differentiation with respect to dimensionless time, . In order to apply perturbation techniques to system (3.6), some assumptions on the order of magnitude of the system parameters have to be made. Thus, we assume the damping, the unbalance and the nonlinearity to be small. This is expressed by making the corresponding coefficients proportional to a sufficiently small, positive and dimensionless parameter : (3.7) where parameters with subscript are -independent. It is also assumed that the torque generated by the motor at resonance ( 󰇗 ) is sufficiently small: (3.8) Finally, the slope of the motor characteristic is assumed to be of the order of unity, i.e. independent of : (3.9) This assumption corresponds to what we have called ‘large slope characteristic’. The case of small slope, with proportional to , is treated in Chapters 5 and 6. Taking the proposed scaling (3.7)-(3.9) into account and dropping the subscript ‘ ’ for convenience, system (3.1) takes the form 󰇘 [ 󰇗 ( 󰇗 󰇘 )] 󰇘 ( 󰇗 ) , 󰇘 - (3.10) 3.2 Alternative First Order Averaging 33 3.2 Alternative First Order Averaging Before turning to the treatment of system (3.10) through some perturbation techniques, an alternative averaging procedure is developed in this section, which will be useful in what follows. In order to make the procedure as general as possible, consider a system of the form { 󰇗 , ( )- ( ) 󰇗 ( ) ( ) 󰇗 ( )} (3.11) where and are matrices of constant coefficients and is a scalar constant, bounded away from zero. It will be shown in the next section that system (3.10) can be written in the form (3.11). Suppose that we try to perform a first order averaging on system (3.11) over angle . According to the results explained in Section 2.1, such a technique is not applicable in this case, because the set of variables is not slow. In order to use the theorem of Section 2.1 we would need to have 󰇗 ( ) and 󰇗 ( ), which is not the case. This justifies the introduction of the modified technique presented below. First, the averaged variables are defined as ( ) ∫ ( ) ⁄ ⁄ ( ) ∫ ( ) ⁄ ⁄ (3.12) where ⁄. As illustrated in Fig. 3.5, the effect of the operator defined in (3.12) is to smooth out the short-term fluctuations of each variable, while retaining the long-term behavior. 34 3 The Case of Large Slope: Analytical Approach Fig. 3.5 Definition of the averaged variables Suppose we are interested in the evolution of the averaged variables ( ) and ( ). Then, we can average the first two equations in (3.11), which yields { 󰇗 [  ∫ ( ( ) ( )) ⁄ ⁄] ( ) 󰇗  ∫ ( ( ) ( )) ⁄ ⁄ ( ) } (3.13) where it has been used that the average, as defined in (3.12), is a linear operator (the average of the sum is the sum of the averages). The next step consists in transforming the integrals in (3.13). Since the process is exactly the same for both integrals, we only focus on the first of them. First, we can write 6.4 6.6 6.8 7 7.2 0 0.2 0.4 0.6 0.8 1 x x t 3.2 Alternative First Order Averaging 35 ∫ ( ( ) ( )) ⁄ ⁄ ∫ ( ( ) ( )) ⁄ ⁄ ( ) (3.14) where it has been used the property that, in one period , ( ) can only change by ( ), according to (3.11). Thus, we can write ( ) ( ) ( ). Changing the integration variable from to yields ∫ ( ( ) ( )) ⁄ ⁄ ( ) ∫ ( ( ) ) ( ⁄) ( ⁄) ( ) (3.15) where the last of relations (3.11) has been used ( ( )). The integration limits can also be transformed by using again 󰇗 ( ): ∫ ( ( ) ) ( ⁄) ( ⁄) ( ) ∫ ( ( ) ) ( ) ( ) ( ) (3.16) Finally, as function is -periodic in , we can write ∫ ( ( ) ) ( ) ( ) ( ) ∫ ( ( ) ) ( ) (3.17) By comparing this last expression to definition (2.3), system (3.13) can be rewritten as 8 󰇗 ,  ( )- ( ) 󰇗  ( ) ( )9 (3.18) 42 3 The Case of Large Slope: Analytical Approach can be written as ( ) (3.41) which is the expression of the Slow Manifold: a 2D surface in the 3D phase space. Thus, the first condition is satisfied. 2. For fixed values of and , it is found that ( ) is a globally asymptotically stable fixed point of the 1D system 󰇗 (3.42) provided that assumption holds. Therefore, the second condition is also fulfilled. Once both requirements have been verified, it can be stated that system (3.38) displays two qualitatively different behaviors at two sequential time scales –see Section 2.3–, which correspond to the second and third stages of the original system (3.26). Using the results of Section 2.3, we have that, at the first of these stages – second stage of (3.26)–, the system can be written as ( ) ( ) 󰇗 ( ) (3.43) where it has been taken into account that, at the beginning of stage 2, ( ) and ( ). Then, at this stage, the system is attracted towards the Slow Manifold, with the slow variables remaining nearly constant: ( ) (3.44) 3.3 Perturbation Approach: Derivation of the Reduced System 43 Summing up, the second stage corresponds to a time length ( ), just as the first one. It ends once variable has reached an ( )–distance to ( ). During this phase of the motion, and do not change significantly. Third stage The third stage of the original system (3.26) –which is the second stage of the averaged system (3.38)– occurs at a time scale ( ⁄). This can be easily understood by noticing that, once the system is near the slow manifold, variable becomes slow (introducing (3.41) in (3.38) leads to 󰇗 ( )). Therefore, near the slow manifold, all variables are slow and, as a consequence, the system natural time scale is ( ⁄). By introducing the expression of the slow manifold in (3.38), the equations corresponding to the third phase of the motion are obtained: { 󰇗 ( ) ( ) 󰇗 ( ( ) ) ( ) ( ) ( ) } (3.45) As usual, higher order terms in (3.45) can be eliminated, giving rise to an ( ) approximation for a time length ( ⁄): { 󰇗 ( ) 󰇗 ( ( ) ) ( ) } (3.46) It is convenient to observe that, although (3.46) contains three equations, only two of them are differential equations. Thus, (3.46) represents a 2D autonomous dynamical system. The evolution of and no longer depends on , once is written as a function of and . The last equation is written with the only purpose of tracking the evolution of variable . 44 3 The Case of Large Slope: Analytical Approach In summary, the third stage corresponds to a time length ( ⁄). At this phase of the motion, the averaged system evolves along the slow manifold given by (3.41). Variables , and obey equations (3.46), with ( ) precision. Fig. 3.6 shows a schematic representation of the three different stages of the system dynamics, summing up the results obtained in the present section. Note that, in Fig. 3.6, the use of overbars for the averaged variables is recovered. The most relevant result is that, once the initial transient corresponding to the first two stages has finished, the evolution of variables and is governed by equations (3.46) –within an ( ) error–. From Fig. 3.6, it is clear that suitable initial conditions for system (3.46) are * +. Recalling definition (3.31), this can be written as * ( ) ⁄+, where * + is the set of initial conditions for system (3.26) Fig. 3.6 Overview of the system dynamics, with * + being the solution of system (3.46) with appropriate initial conditions. Stage 3 Stage 2 Stage 1 0 0 0 a   0 0 1 () () () * aa O O O            00 0 0 () () ( ) ( ) * * *, aa a O O O               0   0 0 1 () * aa O       00 0 0 () () () * **, aa a O O             () () ( ) ( ) *, RR R R aa a O O O                (1)O   (1)O   () 1 O   3.3 Perturbation Approach: Derivation of the Reduced System 45 However, we may be interested in a particular set of initial conditions for system (3.10), given as { 󰇗 󰇗 }. It is, then, convenient, to express the initial conditions for (3.46) as functions of the initial conditions for (3.10): √ 󰇗 ( 󰇗 ) 󰇗 (3.47) as can be readily deduced from relations (3.20), (3.25) and (3.31). Recapitulating, we have been able to eliminate from the formulation variable by Averaging, and variable by applying the Singular Perturbation Theory. 3.4 Analysis of the Reduced System This section focuses on the behaviour of system (3.46), once it has been shown to capture, with ( ) precision, the dynamics of the original system (3.10) during the third stage of the motion. Firstly, it is useful to make a comparison between the system under study and its ideal counterpart, where the rotor speed is constant. Clearly, for this ideal case, the equation of motion of the system shown in Fig. 3.1 is given by 󰇘 󰇗 󰇗 (3.48) with 󰇗 fixed. Equation (3.48) describes a Duffing oscillator, subjected to harmonic excitation. This is a very well-known problem, which has been widely studied in the literature (Brennan, Kovacic, Carrella, & Waters, 2008; Fidlin, 2006; Nayfeh & Mook, 1995; Thomsen, 2003). Under the assumptions of small damping, small 46 3 The Case of Large Slope: Analytical Approach nonlinearity, small unbalance and near-resonant excitation ( 󰇗 ), the Averaging Method can be applied to system (3.48), leading to { 󰇗 ( ) 󰇗 ( )} (3.49) where all the parameters and variables are defined as in Sections 3.1 and 3.3. It is easy to verify that system (3.49) is exactly the same as (3.46), with the only difference of replacing ( ) by the constant value . This is a clear illustration of the concept of nonideal excitation. In the ideal case, the rotor speed appears in equations (3.49) as a constant value , externally imposed by the motor. However, in the nonideal case, the rotor speed enters equations (3.46) as a function of the system vibratory motion, ( ). It is also important to observe that an ideal motor displays a vertical static characteristic, corresponding to the limit case . The motor is, then, able to generate any torque for the same rotor speed. This suggests the idea that a real motor with a static characteristic of very large slope (in absolute value) is more likely to behave in an ideal manner than another one with a smaller slope. Fixed points Going back to the objective of analyzing system (3.46), it is first convenient to look for its fixed points, { }: ( ) ( ) (3.50) From the first of equations (3.50), we have 3.4 Analysis of the Reduced System 47 √ (3.51) Combining (3.41), (3.50) and (3.51) yields √ (3.52) Solutions of (3.52), for both values of , give for all the fixed points of (3.46). This can be done analytically, but the expressions become cumbersome and difficult to interpret. An alternative procedure is proposed, which leads to the fixed points of (3.46) in a graphical way. To this end, the last of equations (3.46) can be rewritten as (3.53) where definition (3.41) has been used. Now, recall the last of equations (3.38), which governs the evolution of the rotor speed for the averaged system: 󰇗 ( ) (3.54) In the light of (3.54), (3.53) can be interpreted as an equilibrium between two torques on the rotor. The left hand term in (3.53) represents the driving torque produced by the motor, while the right hand term represents the resisting torque due to vibration. Thus, the fact that the averaged system is on the slow manifold –which is expressed in equation (3.53)– can be understood as a torque equilibrium condition. Equation (3.53), particularized for the fixed point { }, takes the form 48 3 The Case of Large Slope: Analytical Approach (3.55) where (3.50) has been used. We now define the following functions: ( ) ( ) (3.56) Clearly, according to the comments below equation (3.54), represents the driving torque produced by the motor, while corresponds to the resisting torque due to vibration. Then, (3.55) can be rewritten as ( ) ( ) (3.57) which is the torque equilibrium condition, particularized for the fixed point. In order to solve (3.57) in a graphical way, it would be desirable to write both torques explicitly in terms of . However, this would in turn need explicitly writing in terms of , which produces long and complicated expressions. Thus, an implicit procedure for the graphical representation is proposed. Combining (3.50) and (3.51) results in ( ) (3.58) where function ( ) is defined as ( ) √ (3.59) The proposed representation can be constructed as follows: first, graph versus according to (3.56). Then, graph on the same plot the parametric curve given by 3.4 Analysis of the Reduced System 49 * ( ) ( )+, for and ( -. The fact that is strictly positive comes from the definition of as the radius of a polar coordinate transformation – see (3.20)–. On the other hand cannot be greater than 1, according to the first of equations (3.50). The above procedure gives rise to a plot like that shown in Fig. 3.7. Considering equation (3.57), the fixed points can be found as the intersections of the two torque curves. In the particular case displayed in Fig. 3.7, there are three equilibrium points, marked with circles. Note that the curve associated to the vibration torque is composed of two branches, which collide at the maximum of the curve. They correspond to the two possible values of parameter , as specified in Fig. 3.7. Fig. 3.7 Fixed points of system (3.46) We note that the ‘Sommerfed effect’, which was described in the introduction, can be readily explained by using Fig. 3.7. For such an explanation, the interested reader can refer to (Blekhman, 2000; Dimentberg et al., 1997; Kononenko, 1969; Nayfeh & Mook, 1995). -2 0 2 4 6 8 10 12 0 0.5 1 1.5 2 2.5 1z  T v T 1z m T 50 3 The Case of Large Slope: Analytical Approach Stability Analysis Once the fixed points of the reduced system have been obtained, it is convenient to investigate their stability. For a 2D system, this reduces to calculating the trace and determinant of the jacobian matrix, evaluated at the equilibrium point of interest: [ ( ) ] (3.60) where stands for √ . The conditions for a fixed point to be asymptotically stable are ( ) (3.61) ( ) (3.62) After some algebra, these conditions can be expressed as (3.63) { } (3.64) where denotes the slope of the curve at the considered equilibrium point (see Fig. 3.8 and Fig. 3.9), and has the expression 3.4 Analysis of the Reduced System 51 (3.65) as can be deduced from (3.56), (3.59). Conditions (3.63) and (3.64) are now applied to evaluate stability regions in different scenarios. The procedure is as follows. Consider parameters fixed, so that the curve –see Fig. 3.7– is fixed too. Consider a pair of values ( ) which gives a particular curve ( ). The intersections between the two curves represent the equilibrium points of the system. Select one of them –if there are more than one– and let parameters ( ) vary in such a way that the selected equilibrium point remains an equilibrium point. In other words, let parameters ( ) vary so as to make the curve ( ) rotate around the selected equilibrium point, satisfying restriction . Finally, use conditions (3.63) and (3.64) to analyze how the stability of the fixed point is affected by the slope of the motor characteristic. Fig. 3.8 displays the outcome of applying the above procedure for a fixed point located at the left branch of the vibration torque curve ( ). Two scenarios are considered, depending on the sign of slope , evaluated at the fixed point under consideration. It is observed that a change of stability occurs when both torque curves become tangent ( ). This can be shown to correspond to a transcritical bifurcation. Note that, in Fig. 3.8, the motor curve corresponding to has been directly labeled as , instead of ( ). This shortened notation will be widely used in the figures of the document. Fig. 3.9 shows analogous results for a fixed point located at the right branch of the vibration torque curve ( ). The system behavior is richer in this case, since stability may change in two different ways, depending on the comparison where is defined below. 58 3 The Case of Large Slope: Analytical Approach Transformation to the real eigenbasis of matrix A new change of variables, using the real eigenbasis of matrix , is defined: 0 1 0 1 (3.76) where the columns of matrix are the real and imaginary parts of the complex conjugate eigenvectors of , denoted by : 0 1 0 1 0 1 (3.77) with √4 5 (3.78) System (3.73), written in terms of the new variables, takes the form [ 󰇗 󰇗] ([ ]0 1 [ ( ) ( )]) (3.79) where functions and , containing the nonlinear terms of the system, can be written as Taylor series: 3.5 Classification of the Hopf Bifurcations 59 ( ) ∑ ( ) ∑ (3.80) Coefficients and are specified in the Appendix. Note that the system is finally written in the form (2.21). Thus, the result explained in Section 2.4 can be directly applied. Transformation to Normal Form System (3.79) can be transformed to its Normal Form by a standard procedure (Guckenheimer & Holmes, 1983; Kuznetsov, 1998), as described in Section 2.4: 󰇗 (3.81) where parameter is obtained as { , ( ) ( ) -} (3.82) In summary, it can be said that, after a large number of variable transformations, system (3.46) can be written as (3.81), from which it is concluded that the bifurcation is supercritical (subcritical) if ( ). Despite the fact that coefficients and are of rather complicated form, we find –with the aid of software for symbolic computation (Matlab)– that the condition for supercriticality or subcriticality can be expressed in a surprisingly simple manner: 60 3 The Case of Large Slope: Analytical Approach (3.83) Fig. 3.11 Definition of slope From (3.83), it is clear that a nonlinearity of the softening type ( ) is needed to have a supercritical bifurcation. It is also worth noting that conditions (3.83) admit a very clear graphical interpretation. Consider a curve which intersects at the equilibrium point under consideration and also at the peak of curve . Let denote the slope of this particular motor characteristic, as depicted in Fig. 3.11. In order to obtain , the coordinates of the two points defining the straight line are defined below. First, the highest peak of curve can be shown to correspond to . Substituting this condition in (3.56) and (3.59) yields (3.84) T v T P d  3.5 Classification of the Hopf Bifurcations 61 On the other hand, the ( ) coordinates of the equilibrium point under study are directly given in (3.56) and (3.59): (3.85) Then, from (3.84) and (3.85), the expression of can be readily obtained: (3.86) By comparing (3.86) and (3.66), conditions (3.83) can be expressed as (3.87) This last manner of characterizing the bifurcation is certainly appealing from a graphical point of view, since the basic information about the bifurcation can be directly observed from the torque–speed curves, as shown in Fig. 3.12 for two particular examples. 62 3 The Case of Large Slope: Analytical Approach Fig. 3.12 Examples of (a) subcritical and (b) supercritical bifurcations. (a) (b) -4 -3 -2 -1 01234 0 0.1 0.2 0.3 0.4 0.5  T -10 -5 0 5 0 0.1 0.2 0.3 0.4 0.5  T (a) (b) v T P d H d v T P d H d 3.6 Conditions for the System to be Always Attracted by a Limit Cycle 63 3.6 Conditions under which all System Trajectories are Attracted towards a Limit Cycle In Section 3.5, a simple condition has been obtained to ascertain whether the Hopf bifurcation under study is subcritical or supercritical, which in turn allows predicting the kind of limit cycle generated by the bifurcation (see Fig. 3.10). Although this distinction is relevant, it is based on a local analysis and, consequently, it only gives local information about the system behaviour. This is so in two senses: the analysis of Section 3.5 provides insight into the system dynamics - for values of close enough to (results are local in the parameter space) and - for trajectories close enough to the investigated fixed point (results are local in the phase plane). In view of the aforementioned limitations, this section addresses a new global result that complements those of Section 3.5. It will be shown that, under certain circumstances, the Poincaré-Bendixson (P-B) theorem can be used to prove that all trajectories of the system under study are attracted towards a limit cycle. For a brief explanation of the theorem, see Section 2.5. First, it can be easily deduced from (3.46) that 󰇗 (3.88) Let and represent polar coordinates on the phase plane, according to (3.70), and let denote a circle centred at the origin of the phase plane with a radius slightly greater than , say . From (3.88), it can be said that every trajectory starting outside region will enter and remain inside for all subsequent time. Obviously, trajectories starting inside will also remain inside forever. This kind of behavior would present as a suitable candidate for the role of region in the P-B theorem –see Section 2.5 –, if it were not for the presence of fixed points inside . 64 3 The Case of Large Slope: Analytical Approach Consider now the following particular situation: { } (3.89) whose torque curves are depicted in Fig. 3.13. We suppose that the only fixed point of the system is on the right branch of curve and undergoes a Hopf bifurcation. It is also assumed that the actual slope of the motor characteristic is and, therefore, the equilibrium is unstable. First, let us prove that the fixed point is a repeller. Since the equilibrium is already known to be unstable, we only need to prove that it is not a saddle. Let be the jacobian matrix of system (3.46), evaluated at the equilibrium point. Taking into account that a saddle point has two real eigenvalues with different signs, we can state ( ) (3.90) With some simple algebra, it can be shown that, for , condition ( ) can be written as . Then, it is clear that, for a fixed point satisfying (3.89), we have ( ) . Thus, the equilibrium is a repeller. A new region is now defined as minus a circle of infinitesimal radius around the equilibrium point. From the above considerations –all trajectories enter and the fixed point is a repeller–, it is clear that the flow on the boundary of is directed inwards, as depicted in Fig. 3.14. In summary, a closed, bounded region of the phase plane has been obtained, which contains no fixed points and such that all trajectories of the system enter and remain inside forever. Then, all conditions of the P-B theorem are fulfilled, and it can be assured that any trajectory of the system is attracted towards a closed orbit as , if it is not a closed orbit itself. 3.6 Conditions for the System to be Always Attracted by a Limit Cycle 65 Fig. 3.13 Schematic view of the torque curves corresponding to conditions (3.89) Finally, it should be noted that, although the P-B theorem does not guarantee that all trajectories tend to the same closed orbit, all the numerical experiments conducted within this thesis show the presence of only one stable limit cycle, namely that created by the Hopf bifurcation. This suggests that, for a system verifying (3.89), all the system dynamics is attracted towards a unique limit cycle. -8 -6 -4 -2 0 2 4 0 0.1 0.2 0.3 0.4 0.5 T v T d H d  66 3 The Case of Large Slope: Analytical Approach Fig. 3.14 Flow on the boundary of region (dashed), under conditions (3.89) 3.7 Discussion Time Validity A crucial point in any perturbation analysis is the time scale for which the obtained approximate solution is valid. It has been shown in Section 3.3 that the solution given by the reduced system is valid, at least, for a time scale ( ⁄) –see Fig. 3.6–. However, the situation is even better than that. As described in Section 2.1 (Averaging with Attraction), the asymptotic approximations attained through averaging are valid for all time, whenever they are attracted by a stable fixed point or a stable limit cycle. In the latter case, the uniform validity holds for all variables except the angular one, i.e. the variable which measures the flow on the limit cycle. As will be seen later, all the numerical solutions obtained in Chapter 4 fulfill the above condition of attraction. Region 𝑄 Fixed Point 101 . a x y 3.7 Discussion 67 Comparison with other authors’ results In this subsection, the presented approach and results are compared to some proposed by other authors. First of all, as far as the authors know, there has been no attempt in the literature to use the SPT for the analysis of nonideally excited systems. Thus, the analytical procedure addressed in this Chapter appears to be a novel approach to the problem. On the other hand, the possibility of a Hopf bifurcation on the right branch of the vibration torque curve (Fig. 3.9b.2) has been addressed. An important implication of this result is that the stability of the stationary solutions near resonance does not only depends on the comparison between the slopes of the two torque curves ( ), as commonly stated in the literature (Blekhman, 2000; Dimentberg et al., 1997; Kononenko, 1969; Nayfeh & Mook, 1995). Let us try to explain this divergence in the results. Kononenko’s book (Kononenko, 1969) is one of the most relevant references in the subject. He considered several linear and nonlinear systems excited by nonideal motors. By using the averaging method, he was able to analytically investigate the stationary motions of the motor and their stability. His approach was as follows. Considering the rotor speed to be in the vicinity of resonance, he expanded it as 󰇗 (3.91) Thus, he found equations of motion of the form { 󰇗 ( ) 󰇗 ( ) 󰇗 ( ) 󰇗 } (3.92) 74 4 The Case of Large Slope: Numerical Simulations Fig. 4.2 Phase portraits corresponding to parameters (4.1). The fixed points are marked with dots. The dashed loop represents the unstable limit cycle (a) , (b) The dynamical mechanism whereby the limit cycle is destroyed, which turns out to be a homoclinic bifurcation (Kuznetsov, 1998), is shown in Fig. 4.2 and Fig. 4.3. Let us follow the evolution of the phase portrait. From Fig. 4.2 (a) to Fig. 4.2 (b), the Hopf bifurcation takes place: the focus becomes stable, while an unstable limit cycle is born around it. In Fig. 4.3(a), the cycle has swelled considerably and passes close to saddle point . The homoclinic bifurcation occurs when the cycle touches -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -0.8 -0.6 -0.4 -0.2 0 0.2 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -0.8 -0.6 -0.4 -0.2 0 0.2 x y x y (a) (b) 4.1 Global Bifurcations of the Limit Cycles 75 the saddle point ( ), becoming a homoclinic orbit. In Fig. 4.3(b), we have and the loop has been destroyed. Fig. 4.3 Phase portraits corresponding to parameters (4.1). The fixed points are marked with dots. The dashed loop represents the unstable limit cycle (a) , (b) -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -0.8 -0.6 -0.4 -0.2 0 0.2 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -0.8 -0.6 -0.4 -0.2 0 0.2 x y x y (a) (b) S 76 4 The Case of Large Slope: Numerical Simulations It is worth noting that, when the unstable limit cycle exists –namely, for –, it acts as a frontier between the domains of attraction of the two stable equilibrium points of the system –see Fig. 4.2(b) and Fig. 4.3(a)–. Many other cases exhibiting a subcritical bifurcation, which are not shown here, have also been numerically solved. In all of them, the unstable limit cycle has been found to disappear through a homoclinic bifurcation. The Supercritical Case Consider the following set of dimensionless parameters: (4.5) which might be associated to dimensional parameters { ⁄ ⁄ ⁄ } (4.6) with . Equations (3.66) and (3.86) yield the values of slopes and , depicted in Fig. 4.4. (4.7) Criterion (3.87) allows characterizing the bifurcation as supercritical. Then, as represented in Fig. 3.10, it can be assured that a stable limit cycle encircles the unstable equilibrium for , within a certain neighborhood of . As a matter of fact, the results of Section 3.6 can be used here to investigate the range of slopes for which the limit cycle exists. 4.1 Global Bifurcations of the Limit Cycles 77 Consider the curve which intersects at the fixed point under study and is tangent to curve at another point. Let stand for the slope of that particular torque curve, as displayed in Fig. 4.4. Then, it is straightforward to show that, for , conditions (3.89) are fulfilled and, consequently, it can be assured that all system trajectories tend to a periodic orbit. In the case under analysis, we have (4.8) Fig. 4.4 Torque curves corresponding to parameters (4.5) Note that the Poincaré-Bendixson Theorem gives sufficient, but not necessary, conditions for the existence of a stable periodic orbit. Thus, it cannot be deduced from the Theorem whether the limit cycle survives or not when . To the end of answering this question, we resort again to a numerical resolution of system (3.46), for increasing values of . The results are displayed in Fig. 4.5 and Fig. 4.6. -6 -4 -2 0 2 0 0.2 0.4 0.6 0.8 1  T v T H d T d P d 78 4 The Case of Large Slope: Numerical Simulations Fig. 4.5 Phase portraits corresponding to parameters (4.5). The fixed points are marked with dots. The solid loop represents the stable limit cycle ( ) , (b) Let us track the evolution of the phase portrait. In Fig. 4.5(a) we have and all system trajectories are attracted towards the only fixed point of the system. It may seem from Fig. 4.5(a) that trajectories are actually attracted towards a limit cycle surrounding the fixed point. The reason for this false impression is that the attraction of the fixed point is very weak, as it is close to becoming unstable ( is -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 x y x y (a) (b) 4.1 Global Bifurcations of the Limit Cycles 79 close to ). Hence the required time for trajectories to approach the equilibrium is extremely long. Fig. 4.6 Phase portraits corresponding to parameters (4.5), for . The fixed points are marked with dots. The solid loop represents the stable limit cycle Fig. 4.5b corresponds to . The Hopf bifurcation has occurred and, therefore, the focus has lost its stability at the same time that a stable limit cycle has appeared around it. Note that, in Fig. 4.5(b), conditions (3.89) hold. Consequently, all system trajectories are attracted towards a periodic orbit. Actually, Fig. 4.5b can be observed as a particular example of the general picture shown in Fig. 3.14. The numerical results mentioned above are only useful to confirm the analytical developments of previous sections. By contrast, Fig. 4.6 does provide new information about the global dynamics of the system. It shows that the stable limit cycle is destroyed through a saddle-node homoclinic bifurcation (Kuznetsov, 1998), which occurs at . This means that the cycle disappears exactly when conditions (3.89) are not fulfilled anymore. The mechanism is as follows. At a new fixed point, which immediately splits into a saddle and a node, is created through a saddle-node bifurcation. This new equilibrium appears precisely on the limit cycle, transforming it into a homoclinic orbit. What is found at , -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 x y 80 4 The Case of Large Slope: Numerical Simulations as observed in Fig. 4.6, is that the limit cycle has been replaced by a couple of heteroclinic orbits connecting the saddle and the node. It has been shown that, for the particular set of parameters (4.5), conditions (3.89) are necessary and sufficient for the existence of a stable limit cycle. Thus, the periodic orbit never coexists with any other attractor of the system. Nevertheless, it should be stressed that this is not always the case. In fact, cases have also been found where the stable limit cycle is destroyed through a homoclinic bifurcation, just like in the subcritical case. In these situations, the global bifurcation occurs at certain slope and, therefore, the limit cycle coexists with a stable equilibrium for . As an example, consider a case with satisfying . Clearly, according to (3.87), the Hopf bifurcation is supercritical. However, it is not possible for the limit cycle to be destroyed through a saddle-node homoclinic bifurcation, because the saddle and the node are created before the limit cycle. In fact, in these cases, the closed orbit has been found to die in the same way as the unstable limit cycle shown in Fig. 4.3, i.e. through a homoclinic bifurcation due to the presence of a saddle point. In summary, the simulations carried out suggest that, while unstable limit cycles are destroyed by homoclinic bifurcations, the stable ones can disappear either through homoclinic bifurcations or saddle-node homoclinic bifurcations. 4.2 Numerical Validation of Analytical Results A Subcritical Case Consider again the set of parameters given at (4.1), which gives rise to a subcritical Hopf bifurcation, as depicted in Fig. 4.2 and Fig. 4.3. Two different scenarios are studied, corresponding to the following slopes of the motor characteristic: 4.2 Numerical Validation of Analytical Results 81 (4.9) By comparing (4.9) with Fig. 4.2 and Fig. 4.3, it can be verified that, for , the system has a stable focus surrounded by an unstable limit cycle, while, at , the focus has become unstable through a Hopf bifurcation. As pointed out in Section 4.1, the unstable limit cycle for is the boundary which separates the basins of attraction of the two attracting fixed points present in the system–see Fig. 4.3(a)–. For , two sets of initial conditions, I.C. (1) and I.C. (2), are selected, outside and inside the limit cycle, respectively: ( ){ } ( ){ } (4.10) Then, by using relations (3.47), corresponding initial conditions for the original system can be computed: ( ) { 󰇗 󰇗 } ( ) { 󰇗 󰇗 } (4.11) Note that this step has not a unique solution, because different sets of original initial conditions can produce the same reduced initial conditions. The obtained numerical solutions are shown in Fig. 4.7, for . A good agreement between solutions of both systems is observed. Clearly, the two considered sets of initial conditions lead the system to different attractors. 82 4 The Case of Large Slope: Numerical Simulations Fig. 4.7 Comparison of numerical solutions of the original (solid line) and reduced (dashed line) systems for parameters (4.1), and (a) Displacements (b) Rotor Speed It is convenient to make here an observation about the size of parameter . The procedure used in Chapter 3 to transform the original system into a simpler reduced system is based on perturbation methods. These techniques are useful for dynamical systems which contain a small parameter , and they explain how such systems behave for a sufficiently small . This means that the smaller is, the more accurate 0 0.5 1 1.5 2 2.5 3 3.5 4 x 104 -1 -0.5 0 0.5 1  u, a 0 1 2 3 4 x 104 0.998 1 1.002 1.004   I.C. (2) I.C. (1) I.C. (2) I.C. (1) (a) (b) 4.2 Numerical Validation of Analytical Results 83 perturbation predictions are. Fig. 4.7 shows that, for the case under consideration, a value of gives a remarkable accordance between solutions of the original and reduced system. As an illustrative example, the same numerical computation is done, for initial conditions I.C. (2) and . This larger gives rise to a less accurate prediction, as displayed in Fig. 4.8. The required to have an accurate result depends on the case under study. For instance, in the following simulation (Fig. 4.9), it was necessary to take for a good matching between solutions of the exact and approximate systems. However, in the majority of simulations conducted within this work, proved to be small enough. Consider now the case where, according to Fig. 4.2(a), the focus is unstable and there is a unique attracting fixed point in the system. Initial conditions ( ){ } (4.12) are selected for the reduced system, from which corresponding initial conditions for the original system can be obtained: ( ) { 󰇗 󰇗 } (4.13) The original and reduced systems are numerically solved with and initial conditions (4.13) and (4.12), respectively. The results are displayed in Fig. 4.9, where it is clearly observed how the system moves away from the unstable focus, as the oscillation amplitude increases, until it is attracted to the stable node. I.C. (1) 90 5 The Case of Small Slope: Analytical Approach torque at resonance–, and replacing assumption (3.9) –large slope of the motor characteristic– with (5.1) The assumption is kept within this chapter. Moreover, we assume , with defined in (3.5). With these new assumptions, the dimensionless equations of motion are 󰇘 [ 󰇗 ( 󰇗 󰇘 )] 󰇘 [ ( 󰇗 ) 󰇘 ] (5.2) where subscript ‘0’ has been dropped for convenience. It is useful to transform system (5.2), according to change of variables { ( ) 󰇗 ( )} (5.3) and define a new variable for the rotor speed: 󰇗 (5.4) Notice that the procedure followed in Chapter 3 is being repeated: a change to polar coordinates is performed by replacing the pair of variables * ( ) 󰇗( )+ with the pair of amplitude–phase variables * ( ) ( )+. Thus, the intermediate steps can be skipped, since they are exactly the same as in Chapter 3. The system, written in the new variables, takes the form 5 The Case of Small Slope: Analytical Approach 91 { 󰇗 ( ) ( ) ( ) 󰇗 , ( ) ( )- ( ) 󰇗 ( ) ( ) ( ) 󰇗 } (5.5) where ( ) ( ) (5.6) Equations (5.5) and (5.6) are analogous to (3.26) and (3.27). A direct inspection of system (5.5) reveals that it contains two non-angular real variables * + which are slow –they evolve with rate ( )– and two angular variables * + which are, in principle, fast –they evolve with rate ( ) unless or –. Hence this is a suitable scenario for averaging over the fast angles. However, in order to average over several angles, the system needs to be written in the form (2.5), (2.6), as explained in Section 2.1. To this end, new angular variables are defined: (5.7) Then, by expanding the products of sines and cosines in (5.5), the system can be written as { 󰇗 [ ( ) ( ( ) ( ))] ( ) 󰇗 0 ( ) 1 ( ) 󰇗 ( ) 󰇗 ( ) 󰇗 ( ) } (5.8) 92 5 The Case of Small Slope: Analytical Approach Now, assume a positive rotor speed, . Then, the only resonance manifold present in system (5.8) is given by condition (5.9) As explained in Section 2.1, it is necessary to distinguish between two scenarios, depending on whether or not the system is close to the resonance manifold. 5.1 Outer Region Suppose the rotor speed is away from . Then, we can average system (5.8) over the three fast angles , and . The resulting averaged system is { 󰇗 󰇗 ( )} (5.10) where ( ) ( ) (5.11) According to the averaging theorem stated in Section 2.1, system (5.10) is valid on a time scale ( ⁄), with ( ) precision. A straightforward analysis of system (5.10) yields the conclusion that it has one only fixed point, given by { ( ) } 8  9 (5.12) which is globally asymptotically stable as long as . 5.1 Outer Region 93 Note that, according to assumptions , the equilibrium point (5.12) corresponds to a post-resonant regime, . This solution has a very clear physical interpretation. First, note by comparing (5.11) to (3.3) that function ( ) is simply a dimensionless version of the motor characteristic ( 󰇗): ( ) ( ) (5.13) Clearly,  is the only zero of ( ), as represented in Fig. 5.1. Then, the outer fixed point (5.14) corresponds to a post-resonant motion where the oscillation amplitude is zero and the rotor speed takes the value which makes the motor torque vanish. Note that this holds for the averaged system (5.10), whose solutions are at an ( ) distance to those of the original system (5.2). Then, regarding the original system, it can be said that the outer fixed point (5.14) represents a post-resonant motion with small oscillation amplitudes ( ) and with the rotor speed close to the zero of the motor characteristic  ( ). Fig. 5.1 Dimensionless motor characteristic ( ) In physical terms, it is consistent that a non-resonant excitation produces a small oscillation of the vibrating system. Note that, for such small oscillations, the torque   () m H 1 c 1 d 94 5 The Case of Small Slope: Analytical Approach on the rotor due to vibration is very small ( 󰇘 ( )). Then, during this post-resonant motion, the motor does not have to provide any significant torque to maintain the system vibration and, consequently, the rotor speed takes such a value that the driving torque is virtually zero:  ( ) , ( )- ( ). Two different scenarios can be considered: - If ( ) , system (5.10) is exponentially attracted towards equilibrium (5.12) without approaching the resonance manifold. Then, according to Section 2.1 (Averaging with Attraction), the outer averaged system is valid for all . - If ( ) , system (5.10) is also exponentially attracted by equilibrium (5.12). However, in its way towards the equilibrium, the system will necessarily reach the vicinity of the resonance manifold, making system (5.10) no longer valid. There are, in principle, two options: o The system remains close to the resonance manifold for all subsequent time (resonant capture). o The system stays near the resonance manifold for some finite time, after which it continues its evolution towards fixed point (5.12) (passage through resonance). The next Section investigates the dynamics of (5.5) close to the resonance manifold. 5.2 Inner Region In order to study the system behavior in the vicinity of the resonance manifold, the rotor speed is expanded as √ (5.14) 5.2 Inner Region 95 Note that the definition of the detuning variable is not the same as in the case of large slope –compare (5.14) to (3.35)–. The reason is that, under the assumption of small slope made in this chapter, the problem requires a different perturbation approach, which in turn requires a different scaling of the rotor speed. It can be checked that a scaling such as (3.35) would not yield any relevant result in the present case. Replacing (5.14) in (5.5), (5.6) yields { 󰇗 ( ) ( ) ( √ ) 󰇗 √ ( ) ( ) ( √ ) 󰇗 √ , ( )- ( √ ) 󰇗 √ } (5.15) with ( ) ( ) (5.16) Clearly, system (5.15) contains three slow variables * + and a fast rotating phase . It is, then, suitable for a second order averaging procedure. Following the procedure described in Section 2.2, we arrive at averaged system { 󰇗 (  ) 󰇗 √  4    5 󰇗 √ 0  1  󰇗 √  } (5.17) where, with an appropriate relation between the initial conditions for the original and averaged systems, the error estimates are 96 5 The Case of Small Slope: Analytical Approach {  ( )  ( )  √  (  ) ( )  (√ ) } ( √ ) (5.18) Note that the (√ ) terms in the first two of relations (5.18) turn out to be zero in this particular case. Despite the fact that the evolution {   } is independent of  –as is evident, since this is precisely the purpose of averaging–, system (5.17) includes  as a state variable. The reason is that variable ( ) is necessary to construct the error estimates in (5.18). However, in order to investigate the dynamics of the averaged system, it is convenient to rewrite it without the fast angle: { 󰇗 (  ) 󰇗 √  4    5 󰇗 √ 0  1  } (5.19) A direct analysis of (5.19) allows deducing that, if ⁄, there are no fixed points in the inner region. On the other hand, if ⁄, system (5.19) exhibits two fixed points, given by { √ ( ) √ ( ) √ ( )} (5.20) with 5.2 Inner Region 97 { √ ( ) } (5.21) where √ . The two possible values of correspond to the two different equilibrium points. Note that, unlike in the case of large slope –see (3.52)–, the analytical expressions for the fixed points are very simple in the present scenario. However, the torquespeed plots used in Section 3.4 may also be illustrative here and will provide an interesting comparison with the case of large slope. Thus, consider the equilibrium condition applied to variable . From the third of equations (5.19) we have, at first order, (5.22) On the other hand, condition 󰇗 yields (5.23) which allows writing (5.22) as (5.24) Equation (5.24) can be clearly interpreted as a torque equilibrium condition: ( ) (5.25) with 98 5 The Case of Small Slope: Analytical Approach ( ) (5.26) represents the driving torque produced by the motor and corresponds to the resisting torque due to vibration. In order to obtain the usual torque-speed plot, we would need to write and in terms of . Nevertheless, this would in turn require writing in terms of , and then substitute in (5.25). Since this yields very long and cumbersome expressions, we resort to an alternative implicit procedure for the graphical representation. From condition 󰇗 , we have √ (5.27) Then, we can write ( ) (5.28) where function ( ) is defined as ( ) √ (5.29) The proposed representation can be constructed as follows. First, graph versus (in this case, as , a constant function is obtained). Then, represent the parametric curve given by * ( ) ( )+ for and ( -. The fact that is strictly positive comes from its definition as the radius of a polar coordinate transformation –see (5.3)–, while condition (5.23) forbids to be greater than . This procedure gives rise to a plot like that shown in Fig. 5.2. 5.2 Inner Region 99 Fig. 5.2 Fixed points of system (5.19) A direct comparison between Fig. 5.2 and Fig. 3.7 is quite illustrative as for the difference between the two scenarios considered in this thesis. In the case of large (small) slope, the driving torque curve exhibits a slope which is comparable (negligible) with respect to that of the vibration torque curve. Fig. 5.2 shows the existence of two fixed points, as long as ⁄. The stability of these equilibria is now investigated. To this end, the jacobian matrix of system (5.19) needs to be obtained and evaluated at the equilibrium point of interest: √ ( √ ) (5.30) with [ ] [ ] (5.31) The eigenvalues of matrix are now computed. After some algebra, we find, for the right branch of the torque curve ( ), -2 0 2 4 6 8 10 12 0 0.5 1 1.5 2 2.5 1z  T v T 1z m T c 2  106 5 The Case of Small Slope: Analytical Approach where is a real matrix whose columns contain the real eigenvectors of the jacobian of (5.45), evaluated at . In the case of complex conjugate eigenvectors, the real and imaginary parts are stored in different columns of . Clearly, change of variables is composed of two transformations. First, the coordinate system is translated so that the origin coincides with the fixed point. Then, a transformation to the real eigenbasis of the jacobian is performed, by means of matrix . Vector contains the state variables of the system, expressed with respect to the basis formed by the eigenvectors of the jacobian. Then, it is clear that these eigenvectors are necessarily orthogonal in the space of coordinates . Consequently, for two solutions starting in the Poincaré-Lyapunov domain, we have ‖ ( ) ( )‖ ‖ ( ) ( )‖ (5.53) By considering a time increment ⁄ in (5.53), we have ‖ ‖ ‖ ( ) ( )‖ (5.54) where ( ) (5.55) Thus, it has been shown that ( ) exhibits exponential contraction on the time scale ⁄, even though this contraction is weak. Now, variable ( ) –solution of the original system (5.35)– can also be transformed according to (5.52): (5.56) 5.2 Inner Region 107 Consider the following partition of time in intervals of ( ⁄): [ ] [ ] [ ] ( ) (5.57) where constants are chosen such that 4 ( ) ( )5 (5.58) Note that this choice of constants is always possible, because vanishes at least once per period, according to (5.44). We define ( ) ( ) ( ( ) ) (5.59) ( ) ( ) ( ( ) ) (5.60) ( ) ( ) ( ) (5.61) Thus, ( ) represents the value of at the end of an interval when, as an initial condition, is imposed to be equal to at the beginning of the interval. First, we have that ‖ ( )‖ (5.62) as is clear from (5.37),(5.49) and (5.56), using that at the considered instants. 108 5 The Case of Small Slope: Analytical Approach On the other hand, by virtue of relation (5.54), we can write ‖ ( )‖ ‖ ‖ (5.63) Combining (5.62) and (5.63), we have ‖ ‖ ‖ ‖ (5.64) By using (5.64) recursively, we arrive at ‖ ‖ ( ) ‖ ‖ (5.65) Note that, according to (5.37), (5.49) and (5.56), we can write ‖ ‖ (5.66) thanks to the fact that ( ( ⁄) ( ⁄)) . Finally, introducing (5.66) into (5.65) and taking the limit for yields ‖ ‖ (5.67) where (5.55) has been used. Although (5.67) only holds, in principle, for the particular instants ⁄, it can be readily generalized for any . Note that, as stated in (5.44), any is at ( )– distance from an instant where . Clearly, could be taken as ⁄ and, therefore, (5.67) holds at . On the other hand, and can only undergo ( ) variations in an ( ) time increment, which justifies the generalization of (5.67) to any . Then, for , we can write 5.2 Inner Region 109 ‖ ( ) ( )‖ (5.68) Recovering the original variables, we have ‖ ( ) ( )‖ ‖ ‖‖( ( ) ( ))‖ , ) (5.69) Finally, introducing (5.68) into (5.69) yields ‖ ( ) ( )‖ , ) (5.70) where ‖ ‖ (5.71) The conclusion is that, for initial conditions close enough to the considered equilibrium, the solution of the averaged system is at an ( ) distance from the solution of the original system, for all . Then, if the equilibrium is asymptotically stable in the averaged system, it is asymptotically stable as well for the original system. Now, the obtained result is particularized for the case of the motor with small slope characteristic. Equation (5.70) takes the form {  (√ )  (√ )  (√ )} , ) (5.72) 110 5 The Case of Small Slope: Analytical Approach for solutions starting close enough to the fixed point (5.20), (5.21), with . By comparing (5.70) with (5.18), it is clear that the time validity of the approximation has been extended, paying the price of a less accurate solution. Final Remarks In summary, the system has been found to exhibit two equilibrium points in the resonance region as long as ⁄. Fig. 5.2 represents both equilibria on a torque-speed plot, where the fixed point on the right branch is unstable and the one on the left branch is stable. Note that the existence of a stable fixed point in the inner region justifies the possibility of ‘resonance capture’. Recall that the system reaches the resonance manifold whenever ( ) . For some sets of initial conditions, the trajectory will enter the Poincaré-Lyapunov domain of the stable fixed point and, therefore, it will remain near resonance for all subsequent time –resonant capture–. Clearly, there may also be sets of initial conditions such that the trajectory does not reach the Poincaré-Lyapunov domain of the stable fixed point. In these cases, the system will probably leave the inner region and evolve towards fixed point (5.12) in the outer region –passing through resonance–. In principle, it would also be possible that the system was attracted by a different object in the inner region, such as a stable limit cycle or a chaotic attractor. This would represent another kind of resonance capture, not due to the presence of the stable fixed point analysed in this section. However, the numerical simulations carried out have not revealed the existence in the inner region of any attractor other than the analysed fixed point. Note also that this chapter has coped with the inner and outer approximations separately. We have not tried to construct a ‘composite expansion’ by matching the inner and outer solutions, which is a rather intricate and complex subject, treated, for example, in (W. Eckhaus, 1979; Sanders et al., 2007). 5.2 Inner Region 111 Before showing numerical results to confirm the analytical developments of this Chapter, it is convenient to comment some other works on the subject. Sanders, Verhulst and Murdock considered the system studied in this Section as an illustrative example in Chapters 7 and 8 of their book (Sanders et al., 2007). Regarding the outer region of the phase space, they conducted the same analysis as in this Chapter, averaging over the three fast angles in (5.8) and obtaining equilibrium (5.12). However, in the inner region they only carried out a first order averaging in contrast to the second order averaging addressed in this Chapter. This procedure did not allow them to analyse the stability of the fixed points near resonance, since a first order averaging is not accurate enough for this purpose. In addition, Alexander Fidlin devoted Chapter 5 of his book (Fidlin, 2006) to the study of nonideal excitations, taking system (5.2) as a relevant example. He focused only on the resonance region, arriving at a system analogous to (5.19) after a second order averaging. However, there are two main differences between his results and those presented in the present Chapter: - According to the analysis proposed in this thesis, the equilibrium point for is stable as long as . However, Fidlin came to the conclusion that the condition for stability is . The reason for this difference is a small erratum in the eigenvalues computation in (Fidlin, 2006), namely in the step from equation (5.21) to (5.22) of the mentioned reference: where it reads , it should read . - Fidlin addressed the short time scale where the inner averaged system is valid, ( √ ⁄), as an important limitation of the analysis. He proposed a hierarchic averaging procedure as a way to enlarge the time of validity of the approximation (Pechenev, 1992). However, this method fails precisely in the vicinity of the equilibrium point of interest, because the required variable transformation becomes singular at that point. Therefore, the hierarchic averaging scheme cannot be used to justify the asymptotic stability of the fixed point. 112 5 The Case of Small Slope: Analytical Approach In conclusion, the chief contribution of this Chapter with respect to previous published works is the rigorous justification of the asymptotic stability of one of the stationary motions of the system near resonance, which in turn gives a solid explanation of the possibility of resonant capture, or ‘locking into resonance’. 6 THE CASE OF SMALL SLOPE OF THE MOTOR CHARACTERISTIC: NUMERICAL SIMULATIONS This Chapter is intended to verify the results of Chapter 5 by means of numerical simulation. Following an analogous scheme to that in Chapter 4, particular values are assigned to the system parameters and both the original and approximate systems are numerically solved in order to compare the obtained solutions. Thus, consider the following parameters 114 6 The Case of Small Slope: Numerical Simulations (6.1) which might be associated to dimensional parameters { ⁄ ⁄ ⁄ } (6.2) with . This set of parameters gives rise to the torque-speed curves depicted in Fig. 6.1. Note that condition ⁄ is fulfilled, which implies that there exist two fixed points in the inner region of the phase space, corresponding to the two intersections between and in Fig. 6.1. Fig. 6.1 Torque-speed curves corresponding to parameters (6.1). S and U label the stable and unstable fixed points, respectively. -6 -4 -2 0 2 4 6 8 0.5 1 1.5 2  T m T v T S U 6 The Case of Small Slope: Numerical Simulations 115 As discussed in Chapter 5, the equilibrium on the left branch of curve is stable, while the one on the right branch is unstable. The fixed points can be readily computed by introducing (6.1) into (5.21): { } (6.3) { } (6.4) The third equilibrium, in the outer region of the phase space, can be obtained by introducing (6.1) into (5.12): 2 3 (6.5) The simulations have been carried out as follows. A set of initial conditions for the original system (5.2) is chosen with ( ) , i.e. in the pre-resonant region of the phase space. Then, the original system of equations is numerically solved for a time interval [ ] which is long enough to ascertain whether the system is captured or passes through resonance. Suppose the system passes through resonance. Looking at the numerical solution of the original equations, two particular instants, and are defined, at which the system enters and leaves the resonance region, respectively. Although the choice of these two values is somewhat arbitrary, they give an approximation to the limits between the inner and outer solutions. Then, the outer approximate system (5.10) is solved for , - and [ ], with the initial conditions obtained as the solution of the original system evaluated at and , respectively. The inner approximate system (5.17) is solved for , -, with the initial conditions corresponding to the solution of the original equations particularized at . Finally, the solutions of the original, outer and inner systems are represented 122 6 The Case of Small Slope: Numerical Simulations approximate equations. Once the system has overcome the resonance region, it evolves towards the outer stable equilibrium given by (6.5). It is also illustrative to consider a different situation. Suppose parameter is increased until there are no fixed points in the averaged system near resonance. The condition to meet is ⁄ . Then, consider the following set of parameter values: (6.8) which is exactly the same as (6.1) except for a larger driving torque at resonance. The torque-speed graph for this scenario, obtained through relations (5.26) and (5.29), is depicted in Fig. 6.6, exhibiting no intersections between the curves. Fig. 6.6 Torque-speed curves corresponding to parameters (6.8) Clearly, resonance capture cannot occur in this situation, unless an attractor other than a fixed point existed in the inner region. As stated before, no numerical evidence of such an attractor has been found. The conclusion is that the system will pass through resonance for any pre-resonant initial condition and will lead towards the outer stable equilibrium, which is now given by -6 -4 -2 0 2 4 6 8 0 0.5 1 1.5 2  T m T v T 6 The Case of Small Slope: Numerical Simulations 123 2 3 (6.9) as can be obtained by introducing (6.8) into (5.12): Despite the resonance being not active –i.e. there are no attractors in the resonance region–, it can be expected that trajectories are somehow distorted when passing through the resonance manifold. To the end of observing this effect, two different simulations have been conducted with different initial conditions: { 󰇗 󰇗 } (6.10) { 󰇗 󰇗 } (6.11) The results are displayed in Fig. 6.7-Fig. 6.9. Note that, for initial conditions (6.10), the evolution of the rotor speed is nearly unaffected by resonance, while there is a significant effect on the vibration amplitude. It is interesting that exactly the opposite case is encountered for initial conditions (6.11): whereas the structure vibration is almost unaltered by resonance, the rotor speed undergoes significant oscillations when the system passes through the resonance manifold. Thus, it is clear that the influence of resonance on the system behaviour depends on the initial conditions. In general, we can state that some transient resonant effects can be expected in the system, even when there are no attractors in the resonance region. 124 6 The Case of Small Slope: Numerical Simulations Fig. 6.7 Numerical solutions for parameters (6.1), initial conditions (6.10) and . Solid, dashed and dotted lines correspond to the original, inner and outer systems, respectively a) Displacement b) Rotor speed 0500 1000 1500 2000 2500 -0.1 -0.05 0 0.05 0.1  u, a 0500 1000 1500 2000 2500 0 1 2 3   (a) (b) 1  2  1  2  6 The Case of Small Slope: Numerical Simulations 125 Fig. 6.8 Numerical solutions for parameters (6.1), initial conditions (6.11) and . Solid, dashed and dotted lines correspond to the original, inner and outer systems, respectively a) Displacement b) Rotor speed 0500 1000 1500 2000 2500 -2 0 2  u, a 0500 1000 1500 2000 2500 0 1 2 3   1  2  1  2  (a) (b) 126 6 The Case of Small Slope: Numerical Simulations Fig. 6.9 Close-up around resonance of numerical solutions for parameters (6.1), initial conditions (6.11) and . Solid, dashed and dotted lines correspond to the original, inner and outer systems, respectively a) Displacement b) Rotor speed 150 250 350 450 0.5 1 1.5   150 250 350 450 -3 -2 -1 0 1 2 3  u, a 1  2  1  2  (a) (b) 7 TORQUE-SPEED CURVES FOR THE WHOLE FREQUENCY RANGE This Chapter serves as a connector between the analyses presented in Chapters 3-6 and the study of the vibrocompaction process expounded in Chapter 8. Recall that torque-speed curves have already been successfully used to obtain the stationary motions of the vibrating unbalanced motor in Chapters 3-6. However, in these previous approaches, the torque-speed plot is only represented for the resonance region, like in Fig. 3.7 and Fig. 5.2, or for the non-resonant region, like in Fig. 5.1. This necessary distinction between resonant and non-resonant regions of the phase 128 7 Torque-Speed Curves for the Whole Frequency Range space is a direct consequence of the perturbation approaches that have been utilized in previous chapters. The objective is now to generalize the use of these curves, so that the vibration torque and the motor torque can be plotted together in a single graph for the whole frequency range, thereby representing the resonant and non-resonant stationary motions of the motor. This will turn out to be very useful in the next chapter, as will be seen, and will also provide a clear global perspective of the problem studied in Chapters 3-6. 7.1 Computation of the Torque-Speed Curves Consider again the mechanical system shown in Fig. 3.1, whose equations of motion (3.1) are rewritten below. 󰇘 󰇗 ( 󰇗 󰇘 ) 󰇘 ( 󰇗) 󰇘 (7.1) Note that the cubic nonlinearity is now assumed to be zero for simplicity. As usual, the motor characteristic is assumed to be a linear function of the rotor speed: ( 󰇗) 󰇗 (7.2) with . We note that the analysis shown in what follows is based on Blekhman’s approach of direct separation of motions (Blekhman, 2000). In order to approximately obtain the stationary motions of the system, it is reasonable to look for solutions where the rotor speed has the form 7.1 Computation of the Torque-Speed Curves 129 󰇗( ) ( ) (7.3) with constant and ( ) a periodic function of time with zero average. It is also assumed that the solution satisfies { ( ) 󰇗 ( ) } (7.4) Conditions (7.4) will be verified afterwards. Introducing (7.3) into the first of equations (7.1) yields 󰇘 󰇗 [( ) 󰇗 ] (7.5) By taking (7.4) into account, equation (7.5) can be approximated as 󰇘 󰇗 ( ) (7.6) where ( ) has been assumed for simplicity. Note that, to first approximation, the small oscillation of the rotor speed ( ) does not affect the system vibration. It may seem from (7.6) that the problem has been rendered linear with the proposed approximation. In fact, equation (7.6) represents a harmonically forced linear oscillator. However, the system as a whole is still nonlinear, due to the nonideal interaction with the exciter. This can be seen by noticing that constant in (7.6) is not known a priori. Hence the linear motion ( ) needs to be solved as a function of . Then, the torque produced by this vibration will be introduced in the rotor equilibrium equation –second of equations (7.1)–, which will allow obtaining . Therefore, there is still a two-way coupling between vibration and rotation. The stationary solution of (7.6) is very well-known from linear vibration theory: 130 7 Torque-Speed Curves for the Whole Frequency Range ( ) ( ) (7.7) with . / √6. / 7 0 1 ( . / ) (7.8) where √ ⁄, ( )⁄ . The vibration amplitude is represented against the average rotor speed in Fig. 7.1, according to (7.8). Fig. 7.1 Amplitude of the stationary vibration versus averaged rotor speed, corresponding to equation (7.8) Once the linear motion has been obtained, it can be introduced in the rotor equilibrium equation in order to compute the stationary rotor speed. First, the proposed solution for the rotor speed (7.3) is replaced in the second of equations (7.1): 0 xmax n  7.1 Computation of the Torque-Speed Curves 131 󰇗 ( ) 󰇘 ( ) (7.9) where assumption (7.4) has been used. Then, introducing solution (7.7) into (7.9) yields 󰇗 ( ) ( ) ( ) (7.10) It is convenient to rewrite the last term in (7.10) as the sum of its mean value and an oscillating component: 󰇗 ( ) ( ) (7.11) Note that (7.11) contains constant and oscillating terms with zero average. Clearly, if equation (7.11) is averaged, only the constant terms remain: ( ) (7.12) Substracting (7.12) to (7.11) yields 󰇗 ( ) (7.13) Equation (7.12) can be interpreted as an equilibrium conditions between the average torques acting on the rotor during the stationary motion. Actually, by inserting (7.8) into (7.12), the following relation is obtained: ( ) ( ) (7.14) 138 7 Torque-Speed Curves for the Whole Frequency Range Fig. 7.4 Torque-speed curves corresponding to a case of small slope of the motor characteristic (a) Resonant capture can occur (b) Resonant capture does not occur Torque Torque m L v L Torque n    n  m L v L (a) (b) 7. 2 A Global Perspective for the Cases of Large and Small Slope 139 The near-resonant stationary motions in Fig. 7.4(a) correspond to the fixed points shown in Fig. 5.2. As was widely discussed in Section 5.2, the first of these two fixed points is asymptotically stable, whereas the second one is unstable. The third stationary motion represented in Fig. 7.4(a), which is outside the resonance region, corresponds to the fixed point of the outer averaged system, depicted in Fig. 5.1. This solution was shown to be stable in Section 5.1. Then, for the case of small slope –assuming the motor torque at resonance to be smaller than the resonance peak of the vibration torque curve–, two stable stationary behaviours exist, one of them in the resonance region, and the other being far away from resonance. For any pre-resonant initial state, the system can be attracted by either the near-resonant or the post-resonant stable stationary motions. These two scenarios are referred to as resonant capture and passage through resonance, respectively. In the simpler case represented in Fig. 7.4(b), where no stationary motions close to resonance exist, the system always passes through resonance, and evolves towards its only attractor, away from the resonance region. Hopefully, it has been shown that most of the conclusions about the system behaviour obtained in previous chapters can be summed up and easily retained by using the torque-speed curves described in this chapter. However, the rigorous perturbation approaches of Chapters 3 and 5 are necessary to assess the stability of the stationary motions. 8 MODELLING AND SIMULATION OF THE VIBROCOMPACTION PROCESS This Chapter focuses on the vibrocompaction process which has motivated the whole thesis. After describing the real industrial procedure, a 4-DOF numerical model of the vibrocompaction system is presented. Although this model is suitable for numerical investigation of the process, it is still too complex for an analytical treatment which may reveal more general information about the system dynamics. Then, a second model with 2 DOFs is derived, through some reasonable simplifications, which turns out to be very useful in order to analyse the process and even tune the parameters of the compacting machine. Finally, numerical simulations on the first model (full model) and comparison with the predictions of 142 8 Modelling and Simulation of the Vibrocompaction Process the second (simplified model) illustrate their ability to predict the effect of different parameters on the final level of compaction achieved. 8.1 Some Notes on the Real Process Quartz agglomerates, made of granulated quartz mixed with a polyester resin, are widely used as an artificial stone for countertops in kitchens or bathrooms. The manufacturing process of a slab of this material starts with the filling of a mould with the mixture of quartz and resin. Once the mould is full, a conveyor belt carries it to the vibrocompaction zone, where the thickness of the slab is reduced to nearly half of its initial value, by eliminating the air out of the material. Then, the mixture is cured in a kiln, during a specified time interval, at a suitable temperature for the polymerization of the resin. After the resin is polymerized, an air stream is used to cool the slab before it enters the mechanical finishing stage. During this process the edges are cut, producing a slab of prescribed dimensions, and the surfaces are polished. Then, the product is ready for the quality control stage. It is worth giving some more insight into the vibrocompaction stage of the process, which is the one of interest for the purpose of this study. Before the mixture has been compacted, it is composed of three different phases: solid (the quartz grains), liquid (the resin) and gas (air). The air is present in the material in two different ways: as bubbles within the resin or as gaps between grains of quartz that the resin has not been able to fill. The aim of the compaction process is to eliminate the air out of the mixture, since the presence of pores at the surface of the final countertop is clearly detrimental from a practical point of view: the pores tend to accumulate dirt and are rather difficult to clean. The compaction is conducted by means of several unbalanced electric motors, mounted on a piston with the dimensions of the slab surface. At the beginning of the vibrocompaction process, the piston descends onto the mixture and exerts a static pressure, due to its weight and to an air pressure applied on it. Then, the air 8.1 Some Notes on the Real Process 143 pressure inside the mould is reduced by using a vacuum system, after which the motors are switched on. The vibration produced by the unbalanced motors is the main responsible for the compaction. During the motion of the system, there can be separations and impacts between the piston and the slab, which are generally beneficial for the compaction, as they produce very high peaks of compression forces. In order to reduce vibrations in the vicinity of the compaction machine, elastic elements are placed between the foundation of the machine and the ground, acting as a vibration absorber and thus protecting nearby equipment. Fig. 8.1 shows a pilot plant used for testing purposes, which preserves the main features of the actual industrial machine. It is interesting to note that there are two motors mounted on the piston, which rotate in opposite directions in order to cancel the horizontal components of the centrifugal forces on the unbalanced masses. Hence the net effect of the rotation of both motors is an oscillating vertical force. Fig. 8.1 Pilot plant for the analysis of the vibrocompaction process 144 8 Modelling and Simulation of the Vibrocompaction Process From the above comments, it is clear that the vibrocompaction process is extremely complex from a physical point of view. A large number of factors –some of them being intrinsically nonlinear– influence the final result of the compaction: - The quartz granulometry, the rheological properties of the resin and the mass ratio between quartz and resin affect the mechanical behaviour of the compacting mixture. This behaviour is necessarily nonlinear, since the mixture suffers irreversible deformation during compaction. Moreover, an accurate description of this constitutive law would require modelling the motion of the bubbles through the mixture, the friction between quartz particles, the interaction between quartz and resin, etc. Some investigations about the behaviour these types of three-phase mixtures can be found in (Alonso, Gens, & Josa, 1990; Pietruszczak & Pande, 1996; Stickel & Powell, 2005). - The dynamic properties of the different elements of the machine –the piston, the conveyor belt supporting the mould, the elastomer between the foundation and the ground, etc. – may influence the vibrocompaction as well. - The speed of the motors, their available power and the amount of unbalance are key parameters of the process. - The final result of the compaction may also depend on the duration of the process. - The spatial distribution of the vacuum channels influences the extraction of the air out of the mixture, thereby affecting the compaction. 8.2 Full Model Building a reliable model of such a complex manufacturing process, able to accurately predict the result of the compaction depending on the system parameters, is an extremely hard task, which clearly exceeds the scope of this thesis. It should be noted that, as far as the author know, such a model is not available yet. 8.2 Full Model 145 The aim of this Section is to present an approximate model which, without intending to give accurate quantitative predictions, provides useful qualitative results regarding the vibrocompaction process. This may be seen as a first step towards the ambitious goal of achieving a more complex model which reliably captures the dynamics of the real system. Note that the name full model is used here only for distinction from the simplified model presented in the next section. The simplification carried out can be observed in Fig. 8.2 and Fig. 8.3. The former shows a schematic picture of the real machine, while the later displays the approximate 4-DOF model. The quartz-resin mixture is represented in the model by a couple of masses attached to each other by a linear damper and a nonlinear spring, which models the compaction itself by allowing for permanent deformation when the spring is compressed. Then, the distance between both masses would represent the thickness of the compacting mixture. The mould is modelled as a rigid base, while the piston with the unbalanced motors is represented by a mass with a single unbalanced motor. The mixture is in contact –with separations and impacts allowed– with the mould at the bottom and with the piston at the top. The vacuum system is not included in the model. It should be noted that the model assumes the horizontal motion of the piston to be completely restrained, which makes unnecessary to include a couple of motors rotating in opposite directions. As represented in Fig. 8.3, the model has 4 DOFs: , , and , which correspond, respectively, to position of the bottom of the mixture, position of the top of the mixture, position of the piston and rotation of the motor. The parameters represented in Fig. 8.3 are as follows: stands for the mass of the mixture, is the unbalanced mass, is the mass of the piston and the motor, is the eccentricity of the unbalance, is the rotor inertia, is the damping