Nonlinear dynamics of self-, parametric, and externally excited oscillator with time delay: van der Pol versus Rayleigh models

The regular and chaotic vibrations of a nonlinear structure subjected to self-, parametric, and external excitations acting simultaneously are analysed in this study. Moreover, a time delay input is added to the model to control the system response. The frequency-locking phenomenon and transition to quasi-periodic oscillations via Hopf bifurcation of the second kind (Neimark–Sacker bifurcation) are determined analytically by the multiple time scales method up to the second-order perturbation. Approximate solutions of the quasi-periodic motion are determined by a second application of the multiple time scales method for the slow flow, and then, slow–slow motion is obtained. The similarities and differences between the van der Pol and Rayleigh models are demonstrated for regular, periodic, and quasi-periodic oscillations, as well as for chaotic oscillations. The control of the structural response, and modifications of the resonance curves and bifurcation points by the time delay signal are presented for selected cases.


Introduction
Interactions among different vibration types can occur in various engineering systems. Classical examples are vortex-induced vibrations occurring owing to fluid flow [1], which may interact with other vibrations generated by different mechanisms, such as parametric excitation or external loading. Under specific conditions, different vibration interactions may lead to dangerous resonant states. Interacting self-and parametric vibrations have led to large riser motion being observed in the offshore industry [2,3]. The huge riser vibration was caused by the so-called frequency-locking phenomenon. Self-, parametric, and external oscillations have also been studied for a tower under turbulent flow [4,5]. An interesting example of fluid structure interactions is rotating composite helicopter blades [6], which performed complex bending-twisting oscillations owing to structural couplings. By adding aerodynamic forces and considering the rotation of blades as well as possible flight manoeuvres, interactions may occur between self-excitation and other possible vibrations, such as parametric or external excitations.
An in-depth theoretical study of self-and parametric systems has been presented by Tondl and co-authors. The effect of cancelling self-excitation by parametric vibrations was presented in [7,8]. Detailed bifurcation analysis through Hopf or Neimark-Sacker bifurcations leading to the quenching of self-excited systems was presented by Verhulst [9] and for self-excited autoparametric models by Abadi in a Ph.D. thesis [10].
In the papers of [11,12], new phenomena were detected for nonlinear systems with an external force added to a vibrating self-parametric oscillator. If the frequency of the added force was in the ratio of 1:2 with respect to the frequency of the parametric excitation, the resonance curve changed shape from classical to new, with an internal loop. It has been demonstrated that, in the selected frequency domain, up to five steady states may exist, but stability analysis has indicated that only two of these are stable. This fact has also been confirmed for a two-degree-of-freedom (2-DOF) model for the principal parametric resonances around the first and second natural frequencies [13].
To control the structural response, a model with a time delay was considered with application to the MEMS oscillator [14]. The delay signals (the time delay displacement and velocity) were considered in the study to control the system dynamics. It has been demonstrated that the existence of the internal loop inside the resonance curve of the MEMS device can be controlled by a time delay. Appropriate selection of the gain and time delay may completely reduce or create the loop.
The above-mentioned studies have mainly dealt with the analytical determination of the periodic solutions, bifurcation points, and scenarios of transition from periodic to quasi-periodic oscillations through Hopf (or Neimark-Sacker) bifurcation. However, in recent years, the interactions among various vibration types have been of interest in terms of energy harvesting [15,16]. Authors have considered van der Pol systems with a time-varying delay amplitude coupled with an electromagnetic harvester induced by self-excitations and oscillating near the main parametric resonance. A delayed Mathieu-van der Pol-Duffing oscillator coupled with a piezoelectric mechanism was presented in [17]. To optimise the harvested energy, it was necessary to determine not only the periodic solutions, but also the quasi-periodic oscillations of the harvester. Quasiperiodic solutions have been determined by applying the multiple time scales method twice, firstly to obtain modulation equations representing the slow flow and secondly to obtain periodic solutions of the slow flow. The periodic solutions of the slow flow represent the quasi-periodic motion of the original system. This method has been applied to an externally excited Mathieu oscillator with quadratic and cubic nonlinearities [18] and subsequently to systems with time delay [19], which has previously been studied for periodic response in the paper of [20].
In the above-mentioned literature, self-excited oscillations were studied based on a phenomenological model represented by van der Pol or equivalent Rayleigh nonlinear damping. In both self-excitation models, for the appropriate selection of parameters, the equilibrium position of the system becomes unstable and the solution tends towards a limit cycle (LC), which is stable in such a case. This type of self-excitation is known as soft, in contrast to the so-called hard selfexcitation with an unstable LC. In the latter case, the solution tends towards a stable equilibrium or, under certain initial conditions, approaches infinity. Thus, hard self-excitation is often referred to as catastrophic. Another method of self-excitation modelling is to apply the time delay signal, which may also generate an LC. Therefore, self-excited vibrations are also modelled by differential equations with a time delay, for example in the studies of [21][22][23].
In this study, we investigate a 1-DOF Duffing-type oscillator, which includes all possible excitation types, self-excitation, parametric terms, external forcing, and a time delay signal, which can be treated as a control input. Furthermore, we compare two different selfexcitation types given by van der Pol or Rayleigh terms and demonstrate their differences. The solution of such a general model is determined analytically by the multiple time scales method, taking into account the periodic motion (slow flow) and quasi-periodic oscillations (slow-slow flow). The studied general model can be simplified by equalling certain parameters to zero, and then, the simpler cases can easily be compared with results already published in the literature. This study is a continuation of the paper [14]; however, in this case, the model considers van der Pol or Rayleigh damping and nonlinear parametric excitation producing higherorder terms. The analytical solutions are determined for periodic as well as quasi-periodic oscillations.

Model
Let us consider a model of a structure reduced to an oscillator with 1-DOF. Despite its structural simplicity, we assume that vibrations can be induced by three different mechanisms: self-excitation arising from nonlinear damping, parametric excitation caused by periodically varying stiffness, and external periodic force imposed on the oscillator. All three of these vibration types may interact, in which case various important phenomena may occur. Furthermore, we assume that there exists an additional time delay signal, which is considered as an input of the closed-loop control ( Fig. 1). If the gain g 1 of the input signal is equal to zero, the dynamics of the 1-DOF system (the plant) is not controlled. However, by changing the gain g 1 and time delay τ we may influence the plant dynamics, reduce the oscillations, or induce complex vibrations.
The dynamics of the nonlinear oscillator is governed by a nonlinear ordinary differential equation given in the dimensionless form We assume nonlinear, periodically varied stiffness of the structure, defined by the parameter μ and frequency 2 , which, owing to the coupling with linear and nonlinear components of the stiffness characteristic, generates linear and nonlinear parametric excitation terms. The natural linear frequency of the oscillator is defined by ω 0 and normalised to the dimensionless value ω 0 = 1. However, in the further analytical computations, we maintain the notation ω 0 to observe its influence in the analytical solutions. Cubic nonlinearity is represented by the dimensionless parameter γ . As mentioned previously, the oscillator is also driven by an external periodic force with amplitude f and frequency ω and additionally controlled by a time delay input with gain g 1 and delay τ . The damping of the system is considered as the nonlinear function f d (x,ẋ), producing self-excited vibrations. In the literature, there are two commonly used phenomenological models of self-excited oscillations, namely van der Pol or Rayleigh models, which are defined by the nonlinear functions f v = (−α v +β v x 2 )ẋ and f R = (−α R + β Rẋ 2 )ẋ, with subscripts v and R, respectively. In the literature, both models are often treated as equivalent. In fact, starting with the Rayleigh model and differentiating it with respect to time, we obtain ...

Analytical approach: multiple time scales method
We assume that the differential equation of motion is not strongly nonlinear. Therefore, the nonlinear terms, damping, external, and parametric excitations, as well as time delay, can be expressed by the formal small parameter ε: wherẽ for the van der Pol model and The coefficients in Eq. (2) are expressed by the formal small parameter α = εα, β = εβ, γ = εγ , μ = εμ, f = εf , and g 1 = εg 1 . To simplify the notation, we omit the 'tilde' in further analytical computations.
Owing to nonlinear parametric excitation, the system is excited parametrically in both the first and second perturbation orders; therefore, an analytical solution should be sought, at least up to the second approximation order.
The system response depends on the excitation frequency, and in the case of parametric and external excitations acting simultaneously, we may have an infinite number of combinations of parametric and external excitations leading to parametric and/or external resonances [24]. The parametric resonances occur for m ≈ nω 0 , where n and m are natural numbers.
In this study, we concentrate on the principal parametric resonance, which is the most essential from a practical point of view. Therefore, we assume that ≈ ω 0 . Furthermore, we assume excitation of the external force ω = . This means that parametric and external terms excite the system in the ratio of 2:1 around its linear natural frequency ω 0 . Thus, for the accepted assumptions, we can write where σ 1 is the frequency detuning parameter.
To determine the analytical solution, we apply the multiple time scales method. According to Nayfeh [25,26], the solution is expressed in a series of the small parameter where x j (T 0 , T 1 , T 2 ) and x jτ (T 0 , T 1 , T 2 ) denote the solution and delayed solution, respectively, in the zeroth, first, and second perturbation orders ( j = 0, 1, 2). According to this method, the time is also expressed by scales of the small parameter where T 0 is the fast time scale (fast flow) and T 1 and T 2 are the slow time scales (slow flow). The time definition (5) yields new definitions of the first-and second-order time derivatives where D k j = ∂ k ∂ T j k denotes the kth-order partial derivative with respect to the T j time scale.
Solutions are sought near the principal parametric resonance (3). Thus, substituting (4) into (2) and considering the definitions of the time scales (5) and time derivatives (6), we obtain a set of successive ordinary differential equations in the zeroth, first, and second perturbation orders: The functions F d0 and F d1 are defined in "Appendices A and B" for the van der Pol (40) and Rayleigh (48) models, respectively. The solution to Eq. (7) takes the form where i is the imaginary unit, A is the complex amplitude, andĀ is its complex conjugate. Following substitution of (10) into (8) and grouping the exponential terms, we obtain where cc denotes the complex conjugate functions for the right-side terms and ST 1 represents the secular generating terms of the first order, which must vanish to eliminate the secular terms in the solution. Thus, sub- The functions Q 1 and Q 3 are defined by (41) and (49) for the van der Pol and Rayleigh models, respectively. Next, for ST 1 = 0, we determine the particular solutions to Eq. (11): Substituting solutions (13) into (9), we obtain where N ST 2 are nonsecular generating terms of the second order and, for the sake of brevity, are not reported in this paper, while ST 2 are secular generating terms which we consider as equal to zero: with P 1 defined by (42) and (50) for both damping models.
The complex amplitude A is expressed in the polar form: with a and φ as the amplitude and phase, respectively. Substituting (17) into (16) and separating the real and imaginary parts, we obtain the modulation equations for the amplitude a and phase φ, the so-called slow flow, in the form where W 1 , W 2 , and W 3 are defined by (44) and (52). The approximate solution is constructed based on Eq. (4). Applying (17) and thereafter transforming the exponential into trigonometric functions, the final form of the solution is obtained: with C 1 provided in (45) and (53). The amplitude a and phase φ in the above solution are defined by the modulation equations (18) and (19). In the steady state, da dt = 0 and dφ dt = 0, and thus, by equalling the righthand sides of Eqs. (18) and (19) to zero, we obtain the amplitude and phase for the steady-state solution, which correspond to the periodic oscillations of the original system.

Steady-state periodic oscillations
The system dynamics can be analysed based on the modulation equations (18) and (19). However, we begin from the modulation equations, taking into account the first-order perturbation, and then, we move to the second-order approximation. Therefore, in the first step, we neglect terms of the ε 2 order in the modulation equations (18) and (19). For a steady state, we obtain The above equations allow for preliminary analysis of the system and understanding the meaning of the terms responsible for the self-excitation (parameters α and β included in the variable W 1 ), parametric excitation (coefficient μ), external force (coefficient f ), and time delay (gain g 1 and delay τ ). The above set of equations can be solved analytically for the case without external force, namely f = 0. By eliminating the functions sin 2φ and cos 2φ, we can obtain the resonance curve equation with the coefficients U 1 , U 2 , and U 3 defined in (46) for the van der Pol model and in (54) for the Rayleigh model. From Eq. (22), we can determine the amplitude of oscillations in a steady state for a self-parametric system with time delay. Setting a = 0 in Eq. (22), we can also determine the bifurcation points of the trivial (a = 0) into nontrivial (a > 0) solutions, defined as The bifurcation points of the trivial into nontrivial solutions exist if > 0; thus, we obtain We note that condition (25) does not depend on the β parameters, which means that it does not depend on a model of self-excitation, and the conditions of existence are the same for the van der Pol and Rayleigh models.
In the case without time delay, namely g 1 = 0, we obtain and assuming that α is small, which results in α 4 ≈ 0, we obtain the approximate condition for the existence of the bifurcation points μ α ≥ 2.
For the complete model with self-, parametric, and external force, we must solve the set of Eq. (21). Then, the solution depends on all parameters of the system and has been determined numerically.

Slow-slow motion of self-excited oscillator with time delay, and parametric and external excitations
Modulation equations allow for determining oscillations in a steady state by equalling them to zero. The obtained fixed points correspond to the periodic solutions of the original system (1). Nevertheless, the modulation equations (18) and (19) also provide information on the transient behaviour or periodic motion of the slow motion, which corresponds to the quasiperiodic oscillations of the original system. Therefore, using Eqs. (18) and (19), we may attempt to determine the slow motion of the slow flow, the so-called slow-slow motion. Thus, we expect to obtain analytical solutions for the quasi-periodic oscillations, which play an important role in the studied problem. For this purpose, we apply an extension of the multiple time scales method proposed by Belhaq and Houssni in [18] for parametric and externally excited systems and later also applied to the oscillator with time delay [19] or a more advanced model of a self-excited MEMS device [16,17].
To determine a periodic solution of the modulation equations, we transform them into the Cartesian form by substituting new coordinates Then, by introducing a new small parameter as a bookkeeping device, we obtain Solutions to the equation in (30) are sought in the form of Substituting (31) into (30), taking into account the derivative definitions (6), and collecting the ε 0 and ε 1 terms, we obtain where λ = S 2 2 − S 2 9 represents the LC frequency and depends on the parameters S 2 and S 9 .
The solution of Eq. (32) takes the form Substituting the solutions in (34) into the first-order perturbation, expanding the nonlinear terms into the sine and cosine functions of multiple angles, and then grouping the harmonics, we can determine the secular generating terms, which must be equal to zero. Thus, we obtain modulation equations of the second kind, which represent the amplitude and phase of the quasiperiodic motion dR dt where the coefficients c and d depend on the parameters S 1 . . . S 15 . By setting Eq. (35) equal to zero, we determine the steady state, which is defined by the amplitude The amplitude of the quasi-periodic vibrations a = √ u 2 + v 2 can be computed from (34) in the first approximation order and takes the following form: The envelope of the modulated amplitude is defined by the maximal a max and minimal a min values of Eq. (37). The frequency of a LC corresponding to the quasiperiodic oscillations of the original system is defined by the parameter λ. Considering the definitions of S 2 and S 9 , we obtain λ = 1 64 3 8σ 2 1 + 32σ 1 2 − 3μ 2 + 8α 2 2 + 16g 1 σ 1 + 2 2 cos τ + 8g 2 1 cos 2τ 2 − 256μ 2 4 .
The derived equations facilitate determining the amplitude of slow-slow motion and its frequency λ, which correspond to a quasi-periodic limit cycle of the original system. When analysing Eq. (38), we note that the period of the modulated amplitude T λ = 2π/λ depends on the detuning parameter σ 1 , damping coefficient α, and parametric excitation μ, as well as the delay parameters g 1 and τ . It is interesting to note that it is independent of the nonlinear component of the nonlinear damping β and external force f , taking into account the first-order approximation of the slow-slow motion. This result is in agreement with the results presented in the paper of [15].

Analysis of self-and parametrically excited system with time delay
The analytical solutions obtained in the previous sections enable the analysis of a very general Duffing-type oscillator, which consists of self-excitation, a parametric term, external loading, and time delay input. The approximate solutions are also determined for periodic and quasi-periodic oscillations. By equalling the selected parameters to zero, we can study several interesting dynamic cases. Let us consider the self-and parametric excitation acting on the model. Therefore, we set the external force equal to zero, namely f = 0. For numerical examples, we take the following parameters: The values of the parameters provided in (39) are treated as basic, but several can be varied in specific analyses. The parameter of nonlinear damping β is maintained in an appropriate ratio to compare the van der Pol and Rayleigh models. van   The resonance curves for the principal parametric resonance are obtained based on the second-order approximation by equalling Eqs. (18) and (19) to zero and then determining the amplitude a as a function of the frequency . The comparison of the resonance curves for the van der Pol (black) and Rayleigh (magenta) models is presented in Fig. 2.

(b) (a)
The two considered models exhibit similar qualitative behaviour. The curves represent periodic oscillations, which occur following Hopf bifurcation of the second kind (Neimark-Sacker bifurcation), as indicated by H B v or H B R for the van der Pol or Rayleigh models, respectively. The Hopf bifurcation point occurs earlier, and the resonance zone is wider for the van der Pol model (see Fig. 2b). This fact is also confirmed in Fig. 3

for various parametric excitation parameters μ.
For smaller values of parametric excitation, the resonance zone decreases and the amplitude is smaller. The decreasing tendency of the resonance zones is stronger for the Rayleigh function of self-excitation.
The bifurcation points of the trivial into nontrivial solutions (points at which the resonance curves arise) can be determined from Eq. (26) and depend on the parameters μ and α, but not on β. Thus, considering the first-order approximation, they are located in the same place for the van der Pol and Rayleigh models. Furthermore, they may exist if the condition (27) is satisfied. In fact, for μ = 0.02, the resonance curve arises almost from a single point, as illustrated in Fig. 3a, b by the blue resonance curves. For values of μ smaller than the determined critical μ = 0.02, the resonance curves do not bifurcate from the trivial solutions, but form the shape of a separate island located above the axis (not presented here).
To determine the influence of the parametric excitation amplitude on the system dynamics, we fix the excitation frequency , while maintaining the remaining parameters as provided in (39), and then compute the amplitude versus bifurcation parameter μ. In Fig. 4a, b, we present the bifurcation diagrams for the van der Pol and Rayleigh functions, respectively. The bifurcation curves computed for = 0.95, = 1.0, and = 1.2 begin from unstable periodic solutions, which, while increasing μ, become stable periodic oscillation. Two possible scenarios exist: (1) after supercritical bifurcation of the trivial solution, the curve becomes stable above the Hopf bifurcation point, which corresponds to bifurcation of quasi-periodic to periodic motion, and (2) the curve starts from subcritical bifurcation of the trivial solution and then becomes stable after the limit point (turning point). The first scenario is presented in Fig. 4a for = 0.95 and = 1.0 (black and blue curves, respectively) and Fig. 4b for = 0.95 (black curve). The second scenario is observed in Fig. 4a for = 1.2 (green line) and Fig. 4b for = 1.0 (blue curve) and = 1.2 (green curve). The behaviour of the two models is completely different in the studied cases. The periodic solutions in the Rayleigh model are stable for large values of the parameter μ, while for the van der Pol model, the second turning point occurs and the solutions become unstable.
The periodic solutions presented above exist near the principal parametric resonance zone. When moving beyond this zone, the system response is no longer periodic, but becomes quasi-periodic. To obtain the analytical solutions of quasi-periodic motion, the second type of multiple time scales is applied, as presented in Chapter 5. The maximal and minimal values of the modulated amplitude, as well as the periodic resonance curves, are presented in Fig. 5a. The black lines represent the analytical solutions the for van der Pol model, while the solutions for the Rayleigh model are plotted in magenta. The dynamics of both models is qualitatively similar. However, for the Rayleigh model, the modulations are smaller either before or after the resonance zone, and the transition from quasi-periodic to periodic oscillations begins slightly earlier and ends later for the van der Pol model. This is in agreement with the result presented earlier in Fig. 2b, where a shift in the Hopf bifurcation point for both self-excitation models is illustrated. The numerical verification of the analytical results by means of direct numerical simulation of the original model for the van der Pol self-excitation is presented in Fig. 5b. The blue vertical lines indicate the maximal and minimal values of the modulated amplitude, while the blue dots represent the amplitudes of the periodic motion. The analytical prediction of the periodic solutions exhibits very strong agreement with the numerical simulations, while the analytical results of the quasi-periodic motion exhibit larger differences. This can be explained by the fact that the slow-slow flow is solved in the first-order approximation, while the periodic motion is determined up to the secondorder perturbation terms.
The period of slow-slow motion T λ computed from (38) is independent of the self-excitation type, at least up to the first-order approximation, but depends on the parametric excitation amplitude μ. When approaching the resonance zone, the modulated amplitude period tends towards infinity (Fig. 6). For higher values of parametric excitation μ, the increasing tendency of the period T λ is higher, as indicated by the blue curve in Fig. 6. The period of slow-slow motion for μ = 0.2 obtained from the analytical formula is presented in Fig. 6 by a black solid line, while the results of the corresponding direct numerical simulation are marked by black dots. The gain of the input signal is equal to zero in the mentioned solutions, namely g 1 = 0. The solution is in strong agreement with a slight shift of the analytical curves into the lower frequency direction. The solutions for the period T λ do not exist in the zone near = 1, as in this domain, only periodic motion may exist. The analytical prediction of the amplitude and period of the quasi-periodic motion enables identification of the most important features of the system dynamics. The analytical results are further verified by a bifurcation diagram (Fig. 7) plotted based on the direct numerical simulation of Eq. (1). The solution is computed starting from various basins of attraction, and the transient response is neglected by the rejection of 200 periods. The stroboscopic effect based on the excitation frequency is used to create the diagram. The periodic response versus excitation frequency is represented by the double solid line, which indicates the subharmonic solution in the principal parametric resonance zone. The quasi-periodic solutions with modulated amplitudes are represented by the dark areas. The transition from quasi-periodic motion to the socalled frequency-locking zone with periodic oscillations for the van der Pol oscillator is demonstrated by Poincaré sections in Fig. 8. Starting from low frequencies and moving forward in the bifurcation diagrams in Figs. 5b and 7, we approach the resonance zone. The Poincaré portraits in Fig. 8a-d demonstrate the changes in the system response. For = 0.95, we observe a quasi-periodic LC with an unstable focus (UF) point in the centre. For = 0.962, two unstable focuses arise inside the LC, with a saddle point in the centre, while for = 0.967, the UFs bifurcate into two stable focuses (SFs) inside the LC. A small increase in up to = 0.968 leads to destruction of the LC, and only periodic oscillations represented by double SF may exist. Moving backwards, starting from  Fig. 8e, depending on the initial conditions, the response may be quasi-periodic, represented by a LC with UF in the centre, or periodic, with two SF points located outside of the LC. A small variation in the bifurcation parameter up to = 1.105 destroys the LC, and in Fig. 8f, only the periodic solution exists.
Analysing Eq. (38), we notice that the control parameters g 1 and τ influence the period of slow-slow motion. Thus, we consider the quasi-periodic system response before and after the resonance zone. We fix the excitation frequency as = 0.8 and = 1.2 and vary the control parameters g 1 and τ . The three-dimensional (3D) plots of the period T λ computed based on the analytical solution (38) are presented in Fig. 9. We may effectively change the period and amplitude of the modulated oscillations; however, the tendency when selecting the control parameters for the slow-slow motion before and after the resonance zone is the opposite. As observed in Fig. 9a ( = 0.8), when varying the time delay τ from zero we decrease the period T λ , while in Fig. 9b ( = 1.2), the period is increased.
The time delay signal enables control of the periodic and quasi-periodic response in terms of the frequency and amplitudes, as well as a shift in the resonance zone. This means we can design the system response according to the required criteria, based on the analytical solution, which may differ, for example for energy harvesting or vibration absorption. In this study, we do not present a detailed analysis of any specific dedicated control strategy, but in Fig. 10 we provide an exam- ple of the possible influence of delay parameters for the van der Pol model. By introducing a small gain g 1 = 0.05 and delay τ = 0.1, we can reduce the peri-odic oscillations slightly and move the resonance curve towards the lower frequencies, in addition to reducing the amplitudes of the quasi-periodic response. This is demonstrated by the magenta curves in Fig. 10. The change in the response for the quasi-periodic oscillations is confirmed by the time histories computed for selected excitation frequencies = 0.9 and = 1.15 by direct numerical integration of Eq. (1). The time histories in black in Fig. 11a-c correspond to a system without control, while those in magenta, (b) and (d), correspond to a system with control. The results demonstrate very strong agreement with the analytical solution presented in Fig. 10.
The obtained analytical solutions provide very deep insight into the system dynamics and its possible control; however, they are limited to mainly weakly nonlinear systems. To determine the difference between the van der Pol and Rayleigh self-excitation for a strongly nonlinear model and large oscillations, we perform numerical analysis of Eq. (1). The two bifurcation diagrams in Fig. 12 illustrate the difference in response of the self-and periodically excited oscillator when the parametric excitation amplitude is increased with a frequency = 1 and the remaining parameters are fixed, as given in (39). For very small values of μ, we obtain the quasiperiodic oscillations represented by a black zone on the left-hand side of the bifurcation diagrams in Fig. 12a, b, which, by means of the second type of Hopf bifurcation, bifurcates into periodic, subharmonic motion, represented by the double line. For large values of the bifurcation parameter μ ≈ 1.4, period-doubling bifurcation occurs. However, above this point, the bifurcation scenarios differ for varying self-excitation models. For the van der Pol model, after period doubling, a cascade of period doubling takes place and the system enters chaotic motion (Fig. 12a). Then, after the boundary crisis bifurcation, the system returns to regular subharmonic oscillations. The transition from chaos to regular motion is illustrated in Fig. 13. For μ = 1.55 (Fig. 13a), a regular attractor, represented by the double point, coexists with a strange chaotic attractor located in the middle of the phase plane, which almost touches the boundaries of its attraction basins. A small change in the bifurcation parameter μ = 1.60 destroys the chaotic attractor (owing to the boundary crisis), and the periodic oscillation is maintained, with a fractal basin of attractions in the middle (Fig. 13b).
The analysed 1-DOF nonlinear oscillator with Rayleigh self-excitation does not demonstrate a chaotic response for the range of studied parameters. In contrast, chaotic motion is detected for the van der Pol oscillator with equivalent parameters.

Influence of external force on self-and parametrically excited system
We now consider the self-, parametric, and externally excited system. We assume that the system is excited by three different mechanisms simultaneously. There- fore, an external force occurs on the right-hand side of Eq. (1). We assume that the parametric term and external force excite the system in a ratio of 2:1. This means that we analyse the oscillations near the principal resonance zone, with the additional external force having a frequency corresponding to the subharmonic system response. Such a situation is common in numerous mechanical systems [14]. The influence of the external force is illustrated in Fig. 14 for the van der Pol model. By imposing a very small force f = 0.01 on the system, we obtain a resonance curve with an internal loop and two Hopf bifurca-tion points, located very close to one another (Fig. 14a). This is a qualitatively new response with five periodic solutions; however, only the two upper branches are stable, while the three lower ones are unstable (dashed line). The periodic solutions beyond the resonance zone are unstable. The phenomenon of the resonance curve with the loop and double Hopf bifurcation point is observed for very small values of the external force f . If the force f is increased, the loop becomes smaller and the Hopf bifurcation point located on the left side of the loop disappears (Fig. 14b). A further increase in the parameter f eliminates the loop and stabilises part of the lower branch of the resonance curve, with a new Hopf bifurcation point occurring on the right side of the resonance curve. This situation of the loop almost disappearing is illustrated in Fig. 14c and a zoom presenting the very small loop in Fig. 14d. The external force influences both self-excited oscillators in a similar manner. The resonance curves for f = 0.05 are presented in Fig. 15 for the (a) van der Pol and (b) Rayleigh models. In both cases, we can observe the existence of five solutions, with ∼ (1.05, 1.32) for the van der Pol and ∼ (1.07, 1.23) for the Rayleigh models. The stability analysis based on the modulation equations (18), (19) and then confirmed numerically demonstrates that only the two upper branches of the curve are stable (solid line), while the three lower ones are unstable (dashed line). A direct comparison of the resonance curves for f = 0.1 for both models is presented in Fig. 16a. We note that, for the van der Pol model, the Hopf bifurcation point is located further beyond the resonance zone, and the loop is larger compared to that of the Rayleigh model (see black and magenta curves in Fig. 16a). The loop in the resonance zone for the Rayleigh model is smaller and disappears earlier if the amplitude of excitation f is increased, as illustrated in Fig. 16b. The accuracy of the analytical solutions is verified by means of direct numerical simulation of Eq. (1). Bifurcation diagrams of the solution x versus frequency for f = 0.2 are presented in Fig. 17 for the (a) van der Pol and (b) Rayleigh models. The dark zones correspond to quasi-periodic oscillations, while the solid lines indicate stable periodic solutions. The unstable solutions are not presented in the figures. The solutions corresponding to the internal loop are located in the negative parts of the diagrams, which means that they are shifted in phase compared to the excitation phase. The dynamics of the response can be changed by adding a time delay signal. We can modify the solution stability, shift the resonance zone, move the Hopf bifurcation points, and eliminate the additional resonance loop. The influence of the delay parameters on the periodic and quasi-periodic dynamics of the overall system will be elaborated in detail in a separate paper.

Conclusions
The interactions among self-, parametric, and external oscillations have been studied for two different self-excitation models, represented by the van der Pol or Rayleigh functions. For a weakly nonlinear sys-tem and relatively small oscillations, both considered models exhibit quantitatively similar behaviours. The frequency-locking phenomenon is observed for the system with self-and parametric excitations near the principal resonance zone. The system transits from quasi-periodic to periodic oscillations via Hopf bifurcation of the second kind. The periodic and quasi-periodic solutions are obtained by applying the multiple time scales method twice to determine the amplitudes of the periodic motion (slow flow) and next amplitude of the quasi-periodic motion (slow-slow flow). The solutions are in strong agreement with the direct numerical simulations.
In general, the vibration amplitudes and resonance zones are smaller for the Rayleigh model compared to its van der Pol counterpart. Moreover, the Hopf bifurcation points are slightly shifted. The amplitude of the modulated quasi-periodic response demonstrates similar tendencies. The time delay signal enables control of the dynamics of self-parametric systems, with either periodic or quasi-periodic motions. The analytical solutions obtained provide deep insight into the system dynamics and may aid in the design of a dedicated control strategy.
For strongly nonlinear vibrations, the system response differs quantitatively and qualitatively for varying self-excitation functions. The increase in the parametric excitation transits the van der Pol model to chaotic motion by means of period-doubling bifurcation and then escapes from chaotic motion to a periodic response by means of boundary crisis bifurcation. For the Rayleigh model and equivalent parameters, only period doubling is observed. The van der Pol model demonstrates a higher tendency of moving to chaotic oscillations. For the Rayleigh 1-DOF model, chaotic motion is not observed in the studied parameters domain.
The added small external force changes the dynamics essentially. Near the principal parametric resonance zone, additional solutions occur with an internal loop shape. The stability analysis demonstrates that only the upper parts of the resonance curve are stable. If the amplitude of the external force increases, the loop decreases or totally disappears. The effect of the loop exists for both studied models and is quantitatively similar; however, it exhibits a tendency of decreasing amplitudes for the Rayleigh model.