Stability of a nonideally excited Duffing oscillator

Nonlinear Dyn (2022) 110:3075–3095 https://doi.org/10.1007/s11071-022-07791-y ORIGINAL PAPER Stability of a nonideally excited Duffing oscillator Javier González-Carbajal ·Daniel García-Vallejo ·Jaime Domínguez Received: 9 June 2022 / Accepted: 9 August 2022 / Published online: 3 September 2022 © The Author(s) 2022 Abstract This paper investigates the dynamics of a Duffing oscillator excited by an unbalanced motor. The interaction between motor and vibrating system is considered as nonideal, which means that the excitation provided by the motor can be influenced by the vibrating response, as is the case in general for real systems. This constitutes an important difference with respect to the classical (ideally excited) Duffing oscillator, where the amplitude and frequency of the external forcing are assumed to be known a priori. Starting from preresonant initial conditions, we investigate the phenomena of passage through resonance (the system evolves towards a post-resonant state after some transient nearresonant oscillations) and resonant capture (the system gets locked into a near-resonant stationary oscillation). The stability of stationary solutions is analytically studied in detail through averaging procedures, and the obtained results are confirmed by numerical simulations. Keywords Duffing oscillator ·Sommerfeld effect · Nonideal excitation ·Unbalanced motor ·Resonant capture ·Passage through resonance J. González-Carbajal (B)·D. García-Vallejo · J. Domínguez Department of Mechanical Engineering and Manufacturing, University of Seville, Seville, Spain e-mail: [email protected] 1 Introduction The analysis of forced vibrations in mechanical systems and structures is generally based on the hypothesis of ideal excitation, which means that the forcing amplitude and frequency are assumed to be known a priori. However, the ideality assumption fails to explain and predict certain kinds of behaviour found in real systems [25,43]. In these situations, the excitation needs to be considered as nonideal, meaning that there is a reciprocal interaction between excitation and vibrating response. The study of nonideally excited oscillations started with the pioneering work of Sommerfeld in 1904 [43]. The experimental setup consisted in an unbalanced electric motor mounted on an elastically supported table. Sommerfeld measured the rotor speed and the oscillation amplitude, while gradually increasing the input power in order to make the rotor speed pass through the resonance frequency of the vibrating system [30]. As represented in Fig. 1, it was discovered that, with a steady rate of input power increment, the increase in rotor speed was notably slowed down when the system reached the vicinity of the resonance region, while the oscillation amplitude grew considerably. This smooth evolution of the system state was found to be interrupted by an abrupt jump of the rotor speed towards a post-resonant value, together with a sudden decrease in the amplitude of the oscillations. A less pronounced jumped phenomenon occurred when passing through resonance with decreasing input power (see Fig. 1). 123 3076 J. González-Carbajal et al. Actual Evolution Expected Evolution Natural Frequency of the System Fig. 1 Numerical simulation of Sommerfeld’s experiment This atypical behaviour is often referred to as the Sommerfeld effect. The Sommerfeld effect is closely related to the more general notion of resonant capture [32,35,38]. This phenomenon occurs when a nonlinear system becomes locked in the vicinity of a particular surface of the phase space, known as a resonance manifold, which usually corresponds to a critical value of one of the frequencies present in the system. The alternative, known as passage through resonance, occurs when the system undergoes a transient motion near the resonance manifold and, after some time, it sets out towards the nonresonant region of the phase space. In 1992, Rand et al. [22,35] reported an interesting real case of nonideal excitation. They investigated a specific type of failure found in dual-spin spacecrafts when first placed in orbit, known as precession phase lock. Dual-spin spacecrafts consist of two bodies (platform and rotor) connected by a bearing assembly that allows for relative rotation. An internal motor is used to generate a relative spin motion between platform and rotor. When precession phase lock occurs, the energy provided by the motor is found to be channelled into motions other than the intended rotation. Rand et al. concluded that the phenomenon of precession phase lock is mathematically analogous to the resonant capture of an unbalanced motor attached to an elastic support and driven by a constant torque. Quinn et al. investigated approximate analytical methods to distinguish initial conditions leading to capture or passage through resonance [32]. Quinn also extended the analysis to the case of an unbalanced motor supported by orthogonal, linearly elastic supports, constrained to have a planar motion [33]. Fidlin [14] analysed the near-resonant oscillations produced by an unbalanced motor using two different approaches: a standard second-order averaging and a hierarchic averaging procedure, based on a coordinate transformation suggested by Pechenev [31]. Sanders [38] investigated the same system as Fidlin through first-order averaging, looking for conditions leading to resonant capture or passage through resonance. The literature mentioned up to this point is particularly connected to the present work, due to the similarity in the considered system and assumptions. (The relation between these references and the results of our investigation will be highlighted in Sect. 6.) Some additional works of relevance are briefly described in the next lines. Cveticanin et al. [7,8] used averaging techniques to investigate different configurations of nonideally excited systems, including cases of variable mass. The behaviour of elastically supported unbalanced motors with fractional damping was numerically analysed in [2,47]. Bharti et al. [3] addressed the appearance of the Sommerfeld effect in the torsional vibrations of a double-Cardan joint driveline. Drozdetskaya and Fidlin used a particular type of averaging method, namely the averaging technique for partially strongly damped systems, to study the dynamics of an unbalanced motor with a passive self-balancing system [11]. As an alternative to averaging, many researchers have exploited the method of Direct Separation of Motions [4] for the analysis of nonideally excited oscillators [9,21,24,39,41,42,49,50]. Of course, together with the analytical and numerical research, experimental studies are critical to understand nonideally excited oscillations. In this regard, it is worth mentioning some recent works from Varanis et al. [45,46] where different setups—a cantilever beam, a portal frame and a 3-DOF shear-building structure—are used for experimental validation. Also relevant is the experimental demonstration given by Kossoski et al. [26] of the use of shape memory alloys to mitigate the Sommerfeld effect. From the above paragraphs, it is clear that nonideally excited oscillations have attracted wide attention from researchers in the last decades. On the other hand, the Duffing oscillator is surely one of the most studied systems within the field of nonlinear dynamics [44]. The fact that such an apparently simple equation—a linear oscillator with an extra cubic term in the restoring force—can give rise to jumps, chaotic behaviour, 123 Stability of a nonideally excited Duffing oscillator 3077 bursting oscillations, etc., makes the Duffing oscillator a very useful system to understand the nature of these complex dynamic phenomena [27,34]. In addition, the Duffing equation has great practical importance because it is able to capture the nonlinear behaviour of many physical systems: pendulums [37], cables [48], elastic beams undergoing large displacements [19,29,36], nonlinear electrical circuits [27], etc. The present work is intended to analyse how the Duffing oscillator behaves when, instead of being subjected to the usual harmonic forcing with known frequency and amplitude, it is excited by a nonideal energy source. The organization of this article is as follows. Section 2defines the problem under consideration and the underlying assumptions. The analytical treatment of the equations through averaging techniques is presented in Sect. 3. Section 4shows that an extension of the time scale of validity of the averaged system is required to prove stability of near-resonant solutions. Section 5presents some numerical results to confirm and illustrate the analytical developments. A discussion and interpretation of the obtained analytical and numerical results is offered in Sect. 6, which also highlights the relation between the present paper and some relevant literature. Section 7summarizes the work and presents its conclusions. 2 Problem statement and assumptions We investigate the dynamics of an unbalanced electric motor attached to a fixed frame through a nonlinear spring—whose force has linear and cubic components—and a linear damper, as represented in Fig. 2. Variable xstands for the linear motion, φis the angle of the rotor, m1is the unbalanced mass with eccentricity r,m0is the rest of the vibrating mass, I0 is the rotor inertia (without including the unbalance), b is the viscous damping coefficient, and kand λare, respectively, the linear and cubic coefficients of the spring. The equations of motion for this 2DOF system are [13] m¨x+b˙x+kx +λx3=m1r(˙ φ2cos φ+¨ φsin φ) I¨ φ=Lm(˙ φ) +m1r¨xsin φ(1) where m=m0+m1,I=I0+m1r2and a dot represents differentiation with respect to time t. Gravity can Fig. 2 Representation of the system under analysis: an unbalanced motor with known torque–speed curve connected to a fixed frame through a Duffing-like spring and a linear damper be shown to have no relevance for the purposes of this work [10], since the corresponding term in the equations becomes zero after averaging. For this reason, no gravity term has been included in Eq. (1). Function Lm(˙ φ) represents the motor characteristic. We assume this driving torque to be linearly related to the rotor speed: Lm(˙ φ) =C+D(˙ φ−ωn), (2) with C>0, D<0 and ωn=√k/m. It should be pointed out that writing the motor torque as a function of the rotor speed, and not of some internal electrical variables (voltages, currents or magnetic fluxes), corresponds to a quasi-steady state hypothesis.This means that the dynamics of the electrical variables are assumed to be much faster than the dynamics of the mechanical variables. For a thorough discussion on the validity and limitations of this approach, see [6]. It is useful to define the following dimensionless parameters Rm=m1 m,RI=m1r2 I,ξ=b 2√km α=RIRm 2ξ,ρ=λr2 kRm 2ξ2 c=C Iω2 n ,d=D Iωn (3) Let us also define dimensionless versions of displacement xand time t: u=x r·2ξ Rm ,τ=ωnt(4) Nowsystem(1) can be rewritten in dimensionless form: ¨u+u=−2ξ˙u−ρu3+2ξ(˙ φ2cos φ+¨ φsin φ) ¨ φ=c+d(˙ φ−1)+α¨usin φ(5) where a dot now represents differentiation with respect to dimensionless time τ. 123 3078 J. González-Carbajal et al. In order to investigate system (5) using perturbation techniques, some assumptions on the order of magnitude of the parameters are needed. We assume the damping, the unbalance, the spring nonlinearity and the driving torque at resonance to be small. This is expressed by making the corresponding parameters proportional to a sufficiently small, positive, dimensionless parameter ε: ξ=εξ0,α=εα0,ρ=ερ0,c=εc0(6) where parameters with subscript 0 are ε-independent. It can be shown that the system dynamics and the required perturbation approach are very different depending on the order of magnitude of the motor characteristic [16,17]. We distinguish between the cases of large slope (d=O#(1)) and small slope (d= O#(ε))—symbol O#represents a sharp estimate,see Definition 1.4.4 in [38]orDefinition 1.1.4 in [12]. While the case of large slope was investigated by the authors in [16–18], the case of small slope is the focus of the present work: d=εd0(7) A global comparison between the cases of large and small slope can be found in Chapter 7.2 of [15]. 3 Perturbation approach By introducing Eqs. (6) and (7) into Eq. (5), the system can be written as ¨u+u=ε−2ξ˙u−ρu3+2ξ(˙ φ2cos φ+¨ φsin φ) ¨ φ=εc+d(˙ φ−1)+α¨usin φ(8) where subscript 0 has been dropped for convenience. It is useful to transform system (8) according to change of variables u=acos (φ +β) ˙u=−asin (φ +β). (9) Let us also define a new variable for the rotor speed: ≡˙ φ. (10) It can be shown that system (8), written in terms of the new variables, takes the form (see [17] for details) ˙a=−εsin (φ +β)F1(a,β,φ,)+O(ε2) ˙ =ε[c+d( −1)−αasin φcos (φ +β)]+O(ε2) ˙ β=1−−εcos (φ +β) aF1(a,β,φ,)+O(ε2) ˙ φ= (11) where F1=2ξasin (φ +β) −ρa3cos3(φ +β) +2ξ2cos φ. (12) A direct inspection of system (11) reveals that it contains two nonangular real variables {a,}that are slow—they evolve with rate O(ε)—and two angular variables {β,φ}that are, in principle, fast—they evolve with rate O(1)unless ≈0or≈1. Hence, this is a suitable scenario for averaging over the fast angles (see Chapter 7 in [38] for a thorough treatment of this kind of techniques). However, in order to average over several angles, the system needs to be rewritten with the r.h.s. being a sum of functions, each of them depending on only one of the fast angles (see Remark 7.7.1 in [38]). To this end, new angular variables are defined: ϕ1≡β, ϕ2≡φ+β, ϕ3≡2φ+β. (13) Then, by expanding the products of sines and cosines in (11), the system can be written as ˙a=ε−ξa−ξsin ϕ1+F2(a,ϕ 2)−ξ2sin ϕ3+O(ε2) ˙ =εc+d( −1)+α 2asin ϕ1−α 2asin ϕ3+O(ε2) ˙ϕ1=1−+O(ε) (14) ˙ϕ2=1+O(ε) ˙ϕ3=1++O(ε) where F2(a,ϕ 2)=ρa31 4sin (2ϕ2)+1 8sin (4ϕ2) +ξacos (2ϕ2) (15) System (14) is in the appropriate form for averaging over the 3 fast angles {ϕ1,ϕ 2,ϕ 3}. However, this averaging procedure fails in the vicinity of =−1 and 123 Stability of a nonideally excited Duffing oscillator 3079 =1—note that ϕ1or ϕ3is no longer a fast variable when =1or=−1, respectively. Since we assume the rotor speed to be always positive, the only relevant singularity occurs at =1. This can also be expressed by stating that the system has a resonance manifold at =1. Thus, two different regions of the phase space need to be distinguished: far away from the resonance manifold (outer region), system (14) can be averaged over the 3 fast angles, while in the vicinity of the resonance manifold (inner region) a different averaging approach is required. 3.1 Outer region If the rotor speed is away from 1, we can average system (14) over the three fast angles {ϕ1,ϕ 2,ϕ 3}.Note that system (14)isoftheform ˙ x=ε3 i=1Xi(x,ϕ i)+ O(ε2), with x=[a]T. Its corresponding first-order averaged system can be obtained as ˙ y=ε3 i=1Xi(y), with Xi(·)=1 2π2π 0Xi(·,ϕ i)dϕi. This general procedure is described in detail in Section 7.7 of [38]. The resulting averaged system is ˙a=−εξa ˙ =εHm() (16) where Hm() ≡c+d( −1)(17) System (16) approximates system (14) with O(ε)- precision on a time scale 1/ε,asisknownfromaveraging theory [38]. Note that, in the outer region of the phase space, the vibration amplitude and the rotor speed evolve independently of each other. A direct analysis of system (16) yields the conclusion that it has one only fixed point, given by a=0 Hm() =0⇒a=0 =˜ ≡1−c d(18) which is globally asymptotically stable as long as d<0. Note that, according to the assumption that c>0 and d<0, the equilibrium point given in Eq. (18) corresponds to a post-resonant regime (>1), as represented in Fig. 3. Two different scenarios can be considered: Fig. 3 Dimensionless motor characteristic Hm() •If (0)>1, system (16) is exponentially attracted towards equilibrium (18) without approaching the resonance manifold. It can also be proved that, in this case, the outer averaged system is valid for all τ∈[0,∞)(see Sect. 5.5 in [38]). •If (0)<1, system (16) is also exponentially attracted towards fixed point (18). However, in its way towards the equilibrium, the considered trajectory will necessarily reach the neighbourhood of the resonance manifold, making system (16)no longer valid. There are, in principle, two options: – The system remains close to the resonance manifold for all subsequent time (resonant capture or locking into resonance). – The system stays near the resonance manifold for some finite time, after which it continues its evolution towards fixed point (18)(passage through resonance). 3.2 Inner region In order to study the system behaviour in the neighbourhood of the resonance manifold, the rotor speed is expanded as =1+√εσ. (19) Introducing Eq. (19) into system (11) yields ˙a=−εsin (φ +β)F3(a,β,φ)+O(ε√ε) ˙ β=− √εσ −εcos (φ +β) aF3(a,β,φ)+O(ε√ε) ˙σ=√ε[c−αasin φcos (φ +β)]+εdσ+O(ε√ε) ˙ φ=1+√εσ (20) 123 3080 J. González-Carbajal et al. with F3=2ξasin(φ +β) −ρa3cos3(φ +β) +2ξcos φ (21) System (20) contains three slow variables {a,β,σ}and a fast angle φ. It is useful in this case to apply a secondorder averaging procedure, which provides more accurate results than the usual first-order averaging. Observe that system (20)isoftheform˙ x=√εG1(x,φ) + εG2(x,φ)+..., with x=[aβσ]T. Its corresponding second-order averaged system can be obtained as ˙ y= √εG1(y)+εG2 ∗(y), with G1(·)=1 2π2π 0G1(·,φ)dφ and G2 ∗(·)=1 2π2π 0[G2(·,φ)+H(·,φ)]dφ, where H is a term that depends on G1and G1. The details of this general methodology can be found in Section 7.9 of [38]. The obtained approximate system is ˙ ¯a=−εξ ¯a+sin ¯ β ˙ ¯ β=− √ε¯σ+ε3 8ρ¯a2−ξcos ¯ β ¯a ˙ ¯σ=√εc+α 2¯asin ¯ β+εd¯σ ˙ ¯ φ=1+√ε¯σ. (22) Note that an overbar has been used to emphasize the difference between the original variables {a,β,σ,φ} (solutions of system (20)) and the averaged variables {¯a,¯ β, ¯σ, ¯ φ}(solutions of system (22)), whose relation is given by ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a=¯a+O(ε) β=¯ β+O(ε) σ=¯σ+√εα 4¯acos(2¯ φ+¯ β) +O(ε) φ=¯ φ+O(√ε) ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ for 0 ≤τ≤L √ε (23) with Lan ε-independent constant. The general procedure that allows obtaining the error estimates (23)is described in Section 7.9 of [38]. Despite the fact that the evolution of {¯a,¯ β, ¯σ}is independent of ¯ φ—which is precisely the purpose of averaging—system (22) includes ¯ φas a state variable. The reason is that variable ¯ φ(τ) is required to construct the error estimates in (23). However, in order to investigate the dynamics of the averaged system, it is convenient to rewrite it without the fast angle: ˙ ¯a=−εξ ¯a+sin ¯ β ˙ ¯ β=− √ε¯σ+ε3 8ρ¯a2−ξcos ¯ β ¯a ˙ ¯σ=√εc+α 2¯asin ¯ β+εd¯σ. (24) A direct analysis of system (24) reveals that, if c>α/2, there are no fixed points. On the other hand, if c<α/2, the system exhibits two fixed points, given by aeq =a0+√εa1+O(ε) βeq =β0+√εβ1+O(ε) σeq =σ0+√εσ1+O(ε) (25) with a0=2c α,a1=0 β0=tan−1−a0 −zR0,β 1=0 σ0=0,σ 1=zξR0 a0+3 8ρa2 0 (26) where R0≡1−a2 0and z=±1. The two possible values of zcorrespond to the two different equilibrium points. It is now illustrative to represent the obtained fixed points on a torque–speed plot. This will provide a clear comparison between the cases of large and small slope of the motor characteristic. To this end, consider the condition ˙ ¯σ=0 applied to Eq. (24), which yields c=−α 2a0sin β0.(27) This can clearly be interpreted as a torque balance condition, since the terms on the r.h.s. of the third equation in (24) represent the different torques acting on the rotor. On the other hand, condition ˙ ¯a=0 yields sin β0=−a0,(28) which allows writing Eq. (27)as c=α 2a2 0(29) If we now define functions Tm(dimensionless motor torque) and Tv(dimensionless vibration torque) as Tm≡c Tv(a0)≡α 2a2 0,(30) the fixed points are defined by Tm=Tv(a0). (31) 123 Stability of a nonideally excited Duffing oscillator 3081 0 1 0 Dimensionless Torque Fig. 4 Torque–speed curves in the resonance region for ρ>0 (for ρ<0 the curve is inclined to the left). The fixed points of system (24) correspond to the intersections between Tmand Tv In addition, condition ˙ ¯ β=0 yields σ0=0,σ 1=σv(z,a0)≡3 8ρa2 0+ zξ1−a2 0 a0 . (32) The torque–speed curves can be represented as follows. First, graph Tmversus σ1(in this case, since Tm≡c,a constant function is obtained). Then, plot the parametric curve given by {σv(z,a0), Tv(a0)}for z=±1 and a0∈(0,1]. The fact that ais strictly positive comes from its definition as the radius of a polar coordinate transformation—see Eq. (9)—while condition (28)forbids a0from being greater than 1. This procedure gives rise to a torque–speed representation such as the one shown in Fig. 4. The reader interested in the relation between the cases of large and small slope of the motor characteristic may find illustrative to compare Fig. 4with Fig. 4 in [17]. Figure 4shows the existence of two fixed points, as long as c<α/2. The stability of these equilibria is now investigated. To this end, the Jacobian matrix of system (24) needs to be obtained and evaluated at the equilibrium point of interest: Jeq =√εJ1+εJ2+O(ε√ε) (33) with J1=⎡ ⎣ 000 00−1 −α 2a0−zα 2a0R00⎤ ⎦,(34) J2=⎡ ⎢ ⎣ −ξzξR00 −zξR0 a2 0+3 4ρa0−ξ0 00d ⎤ ⎥ ⎦.(35) Frequency Stable Unstable Fig. 5 Stability of near-resonant solutions for a classical (ideally excited) Duffing oscillator, with ρrepresenting the cubic term in the restoring force Now we compute the eigenvalues of matrix Jeq.After some algebra, the following result is obtained for the fixed point at the right branch of the torque curve (z=1): z=1→⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ λ1=−2εξ +O(ε√ε) λ2=√εαa0R0 2+O(ε) λ3=− √εαa0R0 2+O(ε) ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (36) which clearly corresponds to an unstable equilibrium of the saddle type. For the left branch (z=−1)we obtain z=−1→⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ λ1=−2εξ +O(ε√ε) λ2=√εiαa0R0 2+εd 2+O(ε√ε) λ3=− √εiαa0R0 2+εd 2+O(ε√ε) ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ , (37) where irepresents the imaginary unit. According to Eq. (37), the fixed point on the left branch is stable as long as d<0, since all eigenvalues have negative real parts. Note a significant feature of Eqs. (36) and (37): the stability of both fixed points is independent of the cubic nonlinearity. If condition d<0 is met, then the left equilibrium is stable and the right one is unstable, regardless of the value of parameter ρ.Thisisan important difference with respect to the classical Duffing oscillator, for which it is well known that the stability of near-resonant solutions is highly dependent on the cubic nonlinearity [44], as represented in Fig. 5. 123 3082 J. González-Carbajal et al. 4 Analysis of the time scales This section focuses on a very relevant detail about the stability results obtained in Sect. 3. Note that not only give the eigenvalues information about stability, but also about the time needed by trajectories to be attracted to or repelled by a fixed point. For z=1 (right branch of the vibration torque curve), we have Re(λ2)=O(√ε), where λ2is the only eigenvalue with a positive real part—see Eq. (36). Thus, repulsion of trajectories from the fixed point occurs on the time scale 1/√ε. Note that this is precisely the time scale on which the inner approximation is valid, as specified in Eq. (23). Consider now a trajectory of the averaged system (24) starting sufficiently close to the considered fixed point. If this trajectory is tracked for the time length on which the averaged system is valid (O(1/√ε)), it is found to be repelled out of the equilibrium. Therefore, it can be stated that this fixed point, which is unstable in the averaged system (24), is also unstable in the exact system (20). For z=−1 (left branch of the vibration torque curve), the situation is more intricate. From Eq. (37), we have Re(λ1,λ 2,λ 3)=O(ε), with all three eigenvalues having negative real parts as long as d<0. Hence, the attraction of trajectories towards the fixed point takes place on the time scale 1/ε. On the other hand, from Eq. (23), the averaged system is known to be valid on the time scale 1/√ε. Since 0 <ε1, we can write τvalidity τattraction.(38) With a similar reasoning to the one used for z=1, we can consider an arbitrary trajectory of the averaged system (24) starting in the neighbourhood of the fixed point. If we follow the trajectory for the time length on which the averaged system is valid (O(1/√ε)), it is found to be neither attracted nor repelled by the equilibrium, since the time scale is not long enough. On the other hand, if a longer time scale is considered, solutions of the averaged system may not be good approximations to those of the exact system anymore. From the above considerations, it is clear that a question remains open at this point of the analysis. The fixed point on the left branch of the vibration torque curve has been shown to be asymptotically stable in the averaged system (24)ifd<0. However, is it also stable in the exact system (20)? We answer this question based upon the following result. Result 4.1 Assuming d <0, and for initial conditions sufficiently close to the fixed point corresponding to z=−1, solutions of the averaged system (24) and solutions of the exact system (20) satisfy ⎧ ⎪ ⎨ ⎪ ⎩ a=¯a+O(√ε) β=¯ β+O(√ε) σ=¯σ+O(√ε) ⎫ ⎪ ⎬ ⎪ ⎭ for τ∈[0,∞). (39) Note by comparing Eq. (39)toEq.(23) that the time of validity of the averaged approximation has been extended (in the vicinity of the stable fixed point), paying the price of losing an order of magnitude in √εin the accuracy of the estimate. The proof of this result is based on the attraction properties of system (24) in the vicinity of an asymptotically stable fixed point. We note that the use of attraction properties to extend the time of validity of asymptotic approximations is standard in averaging theory (see Chapter 5 in [38]). However, the problem under study presents a particular difficulty in this regard, due to the weakness of the attraction—as has been stressed before in this section, the real part of the eigenvalues in Eq. (37)isnotO(√ε),butO(ε). This makes necessary a modified version of the general theorem presented in Section 5.5 of [38]. In order to avoid making the main body of the paper too long and cumbersome, the details of the proof can be found in Appendix I. In summary, Result 4.1 allows stating that, if the fixed point corresponding to z=−1 is asymptotically stable in the averaged system (24), the stability is preserved in the exact system (20). 5 Numerical simulations This section presents some numerical results intended to confirm and illustrate the analytical developments of previous sections. The simulations are carried out as follows. After assigning specific values to the system parameters, a set of initial conditions for the original system (8)is chosen with (0)<1, i.e. in the pre-resonant region of the phase space. Then, differential equations (8)are numerically solved for a time interval [0,τf]that is long enough to ascertain whether the system is captured or passes through resonance. Suppose the system passes through resonance. Looking at the obtained numerical solution, two par- 123 Stability of a nonideally excited Duffing oscillator 3083 ticular instants τ1and τ2are 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 averaged system (16) is solved for τ∈[0,τ 1]and τ∈ [τ2,τf], with the initial conditions obtained as the solution of the original system evaluated at and τ=0 and τ=τ2, respectively. The inner averaged system (22)is solved for τ∈[τ1,τ 2], with the initial conditions corresponding to the solution of the original system particularized at τ=τ1. Finally, the solutions of the original, outer and inner systems are represented together, in order to confirm the accuracy of the averaged solutions. In the event of resonant capture, instant τ2does not exist. Then, the outer and inner averaged systems are solved for τ∈[0,τ 1]and τ∈[τ1,τf], respectively. It is worth noting that, regarding the inner approximation, we do not represent the solution {¯a,¯ β, ¯σ, ¯ φ} of system (22), but the more accurate solution given by {¯a,¯ β, ¯σ+√εα 4¯acos(2¯ φ+¯ β), ¯ φ}, according to Eq. (23). The selected set of dimensionless parameters is ξ=2,α=4,ρ=−20,c=1,d=−2(40) which might be associated with dimensional parameters m=100 kg k=5·104N/m b=8.94 Ns/m m1=0.1kg λ=−1.6·106N/m3 r=0.1m I=6.25 ·10−2kgm2 C=3.1·10−2Nm D=−2.8·10−3Nms (41) with ε=0.001. This set of parameters gives rise to the torque–speed curves in the resonance region depicted in Fig. 6. Note that condition c<α/2 is fulfilled, which implies that there exist two fixed points in the inner region of the phase space (marked as Aand Bin Fig. 6). As discussed in Sects. 3and 4, the equilibrium on left branch of the vibration torque curve (point A)is -10 -5 0 5 10 1 0 0.5 1 1.5 2 Dimensionless Torque Fig. 6 Torque–speed curves near resonance for parameter values (40) stable—since we have d<0—while the one on the right branch (point B) is unstable. These fixed points can be readily computed by introducing parameter values (40) into expressions (26): Fixed points in the inner region a0=0.71 β0=5.50 σ0=0 a1=0β1=0σ1=−5.75Point A Stable a0=0.71 β0=3.93 σ0=0 a1=0β1=0σ1=−1.75Point B Unstable (42) The third equilibrium, in the outer region of the phase space, can be obtained through Eq. (18): Fixed point in the outer region a=0 =1.5 Stable (43) This scenario is particularly useful to highlight the difference between the classical (ideally excited) Duffing oscillator and its nonideal counterpart. In order to see this in a clear manner, it is convenient to resort to the amplitude–frequency plot, as usually done when analysing the primary resonance of the classical Duffing oscillator [44]. This can be easily done by representing the last of relations (26) for the parameter values given in Eq. (40), as shown in Fig. 7. Note that the two fixed points Aand Bare also represented on this curve. There is a straightforward relation between the torque– speed plot in Fig. 6and the amplitude–frequency (or amplitude–speed) plot of Fig. 7, since the vibration 123 3090 J. González-Carbajal et al. In order to apply a contraction argument, we will need to guarantee that, after a time interval O(1/ε), the norm of the difference between y1and y2has decreased. By introducing t=L/ε into Eq. (63), we have y1(L/ε) −y2(L/ε)≤ky1(0)−y2(0), with 0 <k<1(64) where k=e−εμ0L=1−εμ0L+O(ε2). (65) Let x(t),φ(t)denote the solution of system (51)for initial conditions x(0)=x0,φ(0)=φ0. Likewise, let y(t)denote the solution of system (58) for initial conditions y(0)=x0−εu1(x0,φ 0), with y(0)being within the Poincaré–Lyapunov domain of the fixed point y= yeq. Now, consider the following partition of time: !0,L ε"&!L ε,2L ε"&... &!mL ε,(m+1)L ε"&... m=1,2, ... (66) On each segment Im=mL ε,(m+1)L εwe define y(m)(t)as the solution of system (58) with initial condition y(m)mL ε =xmL ε−εu1xmL ε,φmL ε. (67) By defining ˜ x(t)=x(t)−εu1(x(t), φ(t)), we can rewrite Eq. (67)as y(m)mL ε=˜ xmL ε.(68) For all finite mwe have from Eq. (57) ˜ x(t)−y(m)(t)≤cε2t∈Im,Im−1(69) with can ε-independent constant and where y(m)(t)has been continued on Im−1(the existence properties of the solutions permit this). On the other hand, we can use Eq. (64) to write y(t)−y(m)(t)Im≤ky(t)−y(m)(t)Im−1(70) with 0 <k<1 (observe the notation · Im≡ supt∈Im·). Combining Eqs. (69) and (70) and using the triangle inequality yields ˜ x(t)−y(t)Im≤cε2+ky(t)−y(m)(t)Im−1.(71) We can now combine Eqs. (71) and (69) and apply again the triangle inequality: ˜ x(t)−y(t)Im≤cε2+k˜ x(t)−y(t)Im−1 +k˜ x(t)−y(m)(t)Im−1 ≤cε2(1+k)+k˜ x(t)−y(t)Im−1. (72) Using relation (72) recursively, we obtain ˜ x(t)−y(t)Im≤cε2(1+k)(1+k+k2+...+km). (73) Taking the limit for m→∞yields the following result for t→∞: ˜ x(t)−y(t)≤cε21+k 1−k.(74) Let us now introduce Eq. (65) into Eq. (74): ˜ x(t)−y(t)≤ cε22+εμ0L+O(ε2) εμ0L+O(ε2)=2c μ0Lε+O(ε2), (75) which can be rewritten as ˜ x(t)−y(t)≤c1ε, (76) with c1an ε-independent constant such that c1>2c μ0L. Note by comparing Eqs. (76) and (74) that we have lost an order of magnitude in εin the estimate due to the fact that constant kis close to 1, which is in turn associated to the very weak attraction of the fixed point of interest. Nevertheless, this attraction proves to be sufficient to guarantee the closeness between ˜ x(t)and y(t)for all t. Finally, by recalling definition ˜ x(t)=x(t)− εu1(x(t), φ(t)) and taking into account that function 123 Stability of a nonideally excited Duffing oscillator 3091 u1is bounded for any x∈Dand φ∈S1, we can write x(t)−y(t)≤c2εfor all t∈[0,∞)(77) which completes the proof. Note that Eq. (39)isthe particularization of Eq. (77) to the problem of the nonideally excited Duffing oscillator. Appendix II: Transformation to real modal form This appendix is intended to justify the assumption C=1inEq.(61). Let us first summarize the situation of interest: we investigate a system of the form ˙ y=εX1(y)+ε2X2 ∗(y), (78) whichisassumedtohaveanasymptoticallystablefixed point y=yeq with all eigenvalues of the Jacobian matrix, Jeq =εJ1+ε2J2+O(ε3), being distinct. For any pair of solutions y1(t)and y2(t)starting within the Poincaré–Lyapunov domain of the fixed point, the Poincaré–Lyapunov theorem states [38] y1(t)−y2(t)≤Cy1(0)−y2(0)e−μt,(79) with μ>0, C≥1. There is, however, one specific scenario where relation (79) holds with C=1, as proved in Appendix III: the case in which Jeq is in real modal form [5]. A real matrix Awith distinct eigenvalues is said to be in real modal form if its nonzero elements coincide with its eigenvalues in the following manner. The real eigenvalues are located on the diagonal of A, while the complex conjugate eigenvalues appear as 2 x2 blocks, with the real part on the diagonal terms and the imaginary part, with alternate signs, on the off-diagonal terms. The rest of the elements of Aare zero. For example, a 3 x3 matrix Ain real modal form with one real eigenvalue λ1 and two complex conjugate eigenvalues λ2,3=α±iβ would be A=!λ100 0αβ 0−βα". Since the proof presented in Appendix I relies upon the assumption that relation (79) holds with C=1, it is clear that the obtained result—Eq. (77)—is valid if Jeq is in real modal form. This appendix will show that the conclusions of Appendix I are also valid in the general case. Consider again the situation portrayed in the first paragraph of this appendix, with the Jacobian matrix Jeq assumed to not be in real modal form. Let us perform a set of coordinate transformations. First, we define a new vector of variables zas z≡y−yeq.(80) In the vicinity of the fixed point of interest, system (78) can be written as ˙z=Jeq z+O(z2). (81) We now define Tas a matrix whose columns constitute a real eigenbasis [28]ofJeq. This means that Tis constructed using the eigenvectors of Jeq as follows. The real eigenvectors are directly stored as columns of T, while, for the complex conjugate eigenvectors, the real and imaginary parts are stored in different columns of T. For example, if Jeq had one real eigenvector v1and a couple of complex conjugate eigenvectors v2,3=u±iw, we would have T=[v1|u|w].Now we define a new vector of variables qas z=Tq.(82) Clearly, system (81), written in terms of q, takes the form ˙ q=Beqq+O(q2), (83) with Beq =T−1JeqT. Due to the way in which matrix Beq is defined, it has some relevant properties, which are well known from the theory of linear dynamical systems [5]: •Beq has the same eigenvalues as Jeq. •Beq is in real modal form. Then, according to the result presented in Appendix III, we can write q1(t)−q2(t)≤q1(0)−q2(0)e−μt,(84) with μ>0 and q1(t),q2(t)being any pair of solutions of system (83) starting within the Poincaré–Lyapunov domain of the fixed point q=0. The conclusion is that, using variable qinstead of y, the Poincaré–Lyapunov theorem can be written with C=1. 123 3092 J. González-Carbajal et al. It has been shown that system (78), subjected to the coordinate transformation y→ q, can be written as system (83). Analogously, system (51) can also be transformed according to change of variables {x,φ} → {p,φ}, with x=yeq +Tp.(85) Then, using variables pand qinstead of xand y, we can legitimately follow the steps from Eq. (66)toEq.(77), which clearly yields p(t)−q(t)≤c3εfor all t∈[0,∞)(86) with c3an ε-independent constant. Now the original variables can be recovered: x(t)−y(t)=T(p(t)−q(t)) ≤Tp(t)−q(t).(87) Hence, we have that x(t)−y(t)≤c4εfor all t∈[0,∞)(88) with c4=Tc3. It has been shown that the conclusion of Appendix I is valid whether or not matrix Jeq is in real modal form. The only difference between both cases is the need for an appropriate coordinate transformation. In Appendix I, the Jacobian matrix was assumed to be in real modal form in order to avoid making the proof unnecessarily cumbersome. Appendix III: The Poincaré–Lyapunov Theorem with the Jacobian matrix being in real modal form The Poincaré–Lyapunov Theorem is a classic result on the stability of dynamical systems, whose general formulation and proof can be found in any textbook on the subject, such as [38]. In this appendix, we are concerned with the specific form of the theorem when the Jacobian matrix of the system, at the fixed point of interest, is in real modal form: Theorem A.3.1 (Poincaré–Lyapunov) Consider the system ˙ x=Ax +g(x), x(0)=a.(89) where x,a∈Rn;Ais a constant n xn matrix in real modal form, with all eigenvalues having strictly negative real part. The vector field g(x)is continuously differentiable with respect to xin a neighbourhood of x=0, with g(x)=o(x)as x→0.(90) Then, there exist constants δ, μ>0such that if a<δ x(t)≤ae−μtfor all t ∈[0,∞). (91) Remark A.3.2 The domain a<δwhere the attraction is exponential is called the Poincaré–Lyapunov domain of the system at x=0. Proof Consider first the linear system associated with system (89): ˙ y=Ay,y(0)=a.(92) The solution of Eq. (92) can be written as y(t)=Φ(t)a,(93) where Φ(t)is the fundamental matrix, i.e. the solution of ˙ Φ=AΦ,Φ(0)=I,(94) where Iis the nxn-identity matrix. Recall that matrix Ais assumed to be in real modal form, which, by definition, means that Ais of the form A= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ λ10... 0 0λ2 ... . . .λ m α1β1 −β1α1 α2β2 −β2α2 ... αrβr 0−βrαr ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (95) where λi(i=1,2,...,m) are the real eigenvalues of Aand λj,j+1=αj±iβj(j=1,2,...,r)arethe pairs of complex conjugate eigenvalues of A. 123 Stability of a nonideally excited Duffing oscillator 3093 The solution of Eq. (94)is Φ(t)=eAt,(96) which, for a matrix Asuch as shown in Eq. (95), takes the form Φ(t)= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ eλ1t0... 0 0eλ2t ... . . .eλmt T1(t) T2(t) ... 0Tr(t) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ,(97) where Tj(t)=!eαjtcos(βjt)−eαjtsin(βjt) eαjtsin(β jt)eαjtcos(βjt)", j=1,2,...,r. (98) Now, let us obtain the norm of Φ(t): Φ(t)=sup{Φ(t)χ:χ=1}(99) Introducing Eqs. (97) and (98) into Eq. (99) yields, after some algebra, Φ(t)2=supχ=1{e2λ1tχ2 1 +e2λ2tχ2 2+···+e2λmtχ2 m+e2α1t(χ2 m+1+χ2 m+2) +e2α2t(χ2 m+3+χ2 m+4)+···+e2αrt(χ2 n−1+χ2 n)} (100) An upper bound for the squared norm can be readily provided: Φ(t)2≤e2Ht(χ2 1+χ2 2+···+χ2 n)=e2Ht (101) with H=max{λ1,λ 2,...,λ m,α 1,α 2,...,α r}. On the other hand, a useful lower bound for Φ(t)2 can be found by choosing an appropriate vector χas follows: if H=λi(i≤m), then χ=ei;ifH=αj, then χ=em+2j(vectors ekrepresent the standard basis, with the k-th component being 1 and the rest being zero). With such a choice of vector χ,wehave that Φ(t)2≥Φ(t)χ2=e2Ht (102) Since it has been shown that Φ(t)2≤e2Ht and Φ(t)2≥e2Ht, the plain conclusion is Φ(t)=eHt.(103) Recall that all eigenvalues are assumed to have negative real parts, which means that H<0. Hence, we can write Φ(t)=e−μt,(104) with μ=−H>0. From Eq. (93), an upper bound for the norm of y(t)can be obtained as y(t)=Φ(t)a≤Φ(t)a.(105) Introducing Eq. (104) into Eq. (105) yields y(t)≤ae−μtfor all t∈[0,∞). (106) It should be stressed that, for a system of the form of Eq. (94) with Abeing a generic real matrix with all eigenvalues having negative real parts, the standard estimate is y(t)≤Cae−μtfor some constants μ>0, C≥1 (see Section 5.2 of [38]). What we have justified is that, when Ais in real modal form, this inequality holds with C=1. The next step would be to show that relation (106), which applies to the linear system (94) for any initial condition a, is also valid for the nonlinear system (89) as long as a<δ,forsomeδ>0. Since this second part of the proof is standard, meaning that it does not depend on whether or not matrix Ais in real modal form, it will not be displayed here. The interested reader can consult it, for example, in Section 5.2 of [38].  Lemma A.3.3 Consider two solutions x1(t)and x2(t) for the equation ˙ x=Ax +g(x)(107) for which the conditions of the Poincaré–Lyapunov Theorem A.3.1 are satisfied. Then, starting in the Poincaré–Lyapunov domain, we have x1(t)−x2(t)≤x1(0)−x2(0)e−μtfor all t ∈[0,∞) (108) 123 3094 J. González-Carbajal et al. for some constant μ>0. Proof See Section 5.2 of [38].  References 1. Avanço, R.H., Tusset, A.M., Balthazar, J.M., Nabarrete, A., Navarro, H.A.: On nonlinear dynamics behavior of an electro-mechanical pendulum excited by a nonideal motor and a chaos control taking into account parametric errors. J. Braz. Soc. Mech. Sci. Eng. 40(1), 1–17 (2018) 2. Balthazar, J., Brasil, R., Felix, J., Tusset, A., Picirillo, V., Iluik, I., Rocha, R., Nabarrete, A., Oliveira, C.: Dynamics behaviour of an elastic non-ideal (nis) portal frame, including fractional nonlinearities. In: journal of physics: conference series, vol. 721, p. 012004. IOP Publishing (2016) 3. Bharti, S.K., Samantaray, A.K.: Resonant capture and Sommerfeld effect due to torsional vibrations in a double cardan joint driveline. Commun. Nonlinear Sci. Num. Simul. 97, 105728 (2021) 4. Blekhman, I.I.: Vibrational mechanics: nonlinear dynamic effects, general approach, applications. World Scientific (2000) 5. Chen, C.T.: Linear system theory and design. Oxford University Press (1999) 6. Cunha, A., Pereira, M., Avanço, R., Tusset, A.M., Balthazar, J.M.: On the reduction of nonlinear electromechanical systems. Meccanica pp. 1–19 (2022) 7. Cveticanin, L., Zukovic, M., Balthazar, J.M.: Dynamics of mechanical systems with non-ideal excitation. Springer (2018) 8. Cveticanin, L., Zukovic, M., Cveticanin, D.: Oscillator with variable mass excited with non-ideal source. Nonlinear Dyn. 92(2), 673–682 (2018) 9. Dasgupta, S.S.: Instability attenuation and bifurcation studies of a non-ideal rotor involving time delayed feedback (2021) 10. Dimentberg, M., McGovern, L., Norton, R., Chapdelaine, J., Harrison, R.: Dynamics of an unbalanced shaft interacting with a limited power supply. Nonlinear Dyn. 13(2), 171–187 (1997) 11. Drozdetskaya, O., Fidlin, A.: Passing through resonance of the unbalanced rotor with self-balancing device. Nonlinear Dyn. 106(3), 1647–1657 (2021) 12. Eckhaus, W.: Asymptotic analysis of singular perturbations. Elsevier (1979) 13. El-Badawy, A.A.: Behavioral investigation of a nonlinear nonideal vibrating system. J. Vib. Control 13(2), 203–217 (2007) 14. Fidlin, A.: Nonlinear oscillations in mechanical engineering. Springer (2005) 15. González-Carbajal, J.: Nonlinear vibrations produced by unbalanced motors. Ph.D. thesis, Universidad de Sevilla (2017) 16. González-Carbajal, J., Domínguez, J.: Limit cycles in nonlinear vibrating systems excited by a nonideal energy source with a large slope characteristic. Nonlinear Dyn. 87(2), 1377–1391 (2017) 17. González-Carbajal, J., Domínguez, J.: Non-linear vibrating systems excited by a nonideal energy source with a large slope characteristic. Mech. Syst. Signal Process. 96, 366–384 (2017) 18. González-Carbajal, J., García-Vallejo, D., Domínguez, J.: On the stability and long-term behaviour of structural systems excited by nonideal power sources. In: nonlinear vibrations excited by limited power sources, pp. 191–237. Springer (2022) 19. González-Carbajal, J., Rincón-Casado, A., García-Vallejo, D., Domínguez, J.: Nonlinear solutions for the steady state oscillations of a clamped-free rotating beam. Eur. J. Mech. Solids 91, 104413 (2022) 20. Habib, G.: Dynamical integrity assessment of stable equilibria: a new rapid iterative procedure. Nonlinear Dyn. 106(3), 2073–2096 (2021) 21. Jiang, J., Kong, X., Chen, C., Zhang, Z.: Dynamic and stability analysis of a cantilever beam system excited by a nonideal induction motor. Meccanica 56(7), 1675–1691 (2021) 22. Kinsey, R., Mingori, D., Rand, R.: Nonlinear controller to reduce resonance effects during despin of a dual-spin spacecraft through precession phase lock. In: [1992] Proceedings of the 31st IEEE Conference on decision and control, pp. 3025–3030. IEEE (1992) 23. Kinsey, R., Mingori, D., Rand, R.: Limited torque spinup of an unbalanced rotor on an elastic support. In: Proceedings of 1995 American Control Conference-ACC’95, vol. 6, pp. 4368–4373. IEEE (1995) 24. Kong, X., Li, W., Jiang, J., Dong, Z., Wang, Z.: Dynamic characteristics of a simply supported elastic beam with three induction motors. J. Sound Vib. 520, 116603 (2022) 25. Kononenko, V.O.: Vibrating systems with a limited power supply. Iliffe, London (1969) 26. Kossoski, A., Tusset, A.M., Janzen, F.C., Rocha, R.T., Balthazar, J.M., Brasil, R.M., Nabarrete, A.: Jump attenuation in a non-ideal system using shape memory element. In: MATEC Web of Conferences, vol. 148, p. 03003. EDP Sciences (2018) 27. Kovacic, I., Brennan, M.J.: The Duffing equation: nonlinear oscillators and their behaviour. Wiley (2011) 28. Kuznetsov, Y.A.: Elements of applied bifurcation theory, vol. 112. Springer (2013) 29. Lacarbonara, W., Yabuno, H.: Refined models of elastic beams undergoing large in-plane motions: theory and experiment. Int. J. Solids Struct. 43(17), 5066–5084 (2006) 30. Nayfeh, A.H., Mook, D.T.: Nonlinear oscillations. Wiley (1979) 31. Pechenev, A., Fidlin, A.: Hierarchy of resonant motions excited in a vibroimpact system with contact zones by an inertial source of limited power. Izw. AN SSSR Mechanika Tverdogo Tela 27, 46–53 (1992) 32. Quinn, D., Rand, R., Bridge, J.: The dynamics of resonant capture. In: advances in nonlinear dynamics: methods and applications, pp. 1–20. Springer (1995) 33. Quinn, D.D.: Resonance capture in a three degree-of- freedom mechanical system. Nonlinear Dyn. 14(4), 309–333 (1997) 34. Rakaric, Z., Kovacic, I.: Mechanical manifestations of bursting oscillations in slowly rotating systems. Mech. Syst. Signal Process. 81, 35–42 (2016) 123 Stability of a nonideally excited Duffing oscillator 3095 35. Rand, R.H., Kinsey, R., Mingori, D.L.: Dynamics of Spinup through resonance. Int. J. Non-Linear Mech. 27(3), 489–502 (1992) 36. Rincón-Casado, A., González-Carbajal, J., García-Vallejo, D., Domínguez, J.: Analytical and numerical study of the influence of different support types in the nonlinear vibrations of beams. Eur. J. Mech. Solids 85, 104113 (2021) 37. Salas, A.H., et al.: Exact solution to duffing equation and the pendulum equation. Applied Math. Sci. 8(176), 8781–8789 (2014) 38. Sanders, J.A., Verhulst, F., Murdock, J.: Averaging methods in nonlinear dynamical systems (Second Edition), vol. 59. Springer (2007) 39. Shahlaei-Far, S., Balthazar, J.M.: Method of direct separation of motions applied to a non-ideal electromechanical pendulum system. In: dynamical systems theory and applications, pp. 325–333. Springer (2017) 40. Silveira, M., Gonçalves, P.J., Balthazar, J.M.: Multi-core and many-core SPMD parallel algorithms for construction of basins of attraction. J. Theor. Appl. Mech. (2019). https:// doi.org/10.15632/jtam-pl/112463 41. Sinha, A., Bharti, S.K., Samantaray, A., Chakraborty, G., Bhattacharyya, R.: Sommerfeld effect in an oscillator with a reciprocating mass. Nonlinear Dyn. 93(3), 1719–1739 (2018) 42. Sinha, A., Bharti, S.K., Samantaray, A.K., Bhattacharyya, R.: Sommerfeld effect in a single-dof system with base excitation from motor driven mechanism. Mech. Mach. Theor. 148, 103808 (2020) 43. Sommerfeld, A.: Naturwissenchftliche ergebnisse der neuren technischen mechanik. Verein Deutscher Ingenieur Zeitscchrift 18, 631–636 (1904) 44. Thomsen, J.J.: Vibrations and stability, 3rd edn. Springer (2021) 45. Varanis, M., Balthazar, J., Silva, A., Mereles, A., Pederiva, R.: Remarks on the Sommerfeld effect characterization in the wavelet domain. J. Vib. Control 25(1), 98–108 (2019) 46. Varanis, M., Mereles, A., Silva, A.L., Barghouthi, M.R., Balthazar, J.M., Lopes, E.M., Bavastri, C.A.: Numerical and experimental investigation of the dynamic behavior of a cantilever beam driven by two non-ideal sources. J. Braz. Soc. Mech. Sci. Eng. 42(10), 1–18 (2020) 47. Varanis, M.V., Tusset, A.M., Balthazar, J.M., Litak, G., Oliveira, C., Rocha, R.T., Nabarrete, A., Piccirillo, V.: Dynamics and control of periodic and non-periodic behavior of duffing vibrating system with fractional damping and excited by a non-ideal motor. J. Franklin Inst. 357(4), 2067–2082 (2020) 48. Vassilopoulou, I., Gantes, C.: Nonlinear dynamic phenomena in a SDOF model of cable net. Arch. Appl. Mech. 82(10), 1689–1703 (2012) 49. Yaroshevich, N., Yaroshevych, O., Lyshuk, V.: Drive dynamics of vibratory machines with inertia excitation. In: vibration engineering and technology of machinery, pp. 37–47. Springer (2021) 50. Zhang, X., Li, Z., Li, M., Wen, B.: Stability and sommerfeld effect of a vibrating system with two vibrators driven separately by induction motors. IEEE/ASME Trans. Mech. 26(2), 807–817 (2020) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123