Hybrid rayleigh–van der pol–duffing oscillator: Stability analysis and controller

Introduction

Nonlinear vibration provides an interesting potential example of the mathematical description of the nonlinear behavior of many phenomena in science, physics, and practical engineering; for example, the N/MEMS system vibrates nonlinearly.^1-7 The nonlinear wave equation of Kundu–Mukherjee–Naskar equation can be finally converted into a Duffing-like equation. ⁸ Fangzhou oscillator is a generalized Duffing equation (DE) with a singular term.^9-11 The nonlinear vibration systems in a porous medium can be converted to a fractal modification of the DE. ^12-14 The gecko-like vibration plays an important role in the accurate 3-D printing process. ¹⁵ The inherent pull-in instability of MEMS systems can be completely overcome by the fractal vibration theory.^16-18 The homotopy perturbation method (HPM) ¹⁹ and the Hamitonian approach ²⁰ are two main analytical tools for nonlinear vibration systems. The combination of the Laplace transforms, Lagrange multiplier, fractional complex transforms, and Mohand transform with HPM was employed to find approximate solutions for nonlinear partial differential equations.^21-25 Additionally, there are many other methods available in the literature.

Duffing equation is adopted in many problems such as optical stability, electrical circuit, oscillation of plasma, and the buckled beam. ²⁶ Nonlinear oscillations have been of paramount importance in practical engineering, physics, applied mathematics, and several real-world requirements for many years. Nayfeh ²⁷ and Mickens ²⁸ utilized different approaches to obtain approximate solutions for nonlinear oscillators. Moatimid ²⁹ combined the HPM and Laplace transform to obtain an approximate bounded solution for the parametric DE. This analysis resulted in an exact solution for the cubic DE and showed that the damped parameter and the cubic stiffness parameter have a destabilizing influence on the system. In the last decades, much research was devoted to examining the control of resonantly forced systems in several practical engineering areas. In the passive vibration absorbers, a physical practice is associated with primary organization, whereas in case of active absorbers, the device is replaced by a control system of sensors, actuators, and filters. The Van der Pol (VDP) oscillator is an essential nonlinear oscillator, which has been comprehensively investigated. However, the results of delayed universal VDP oscillators were relatively fewer. The combination of cubic–quintic Duffing–Van der Pol equation (DVdP) and two external periodic forcing terms was scrutinized by Ghaleb et al. ³⁰ In case of the autonomous system near the equilibrium points, the linearized stability was reached. Furthermore, in case of the nonautonomous system, by utilizing the multiple time scales, stability was analyzed. The bifurcation controller for a delayed extended DVdP oscillator was considered by Huang. ³¹ The impression of feedback advantage of the bifurcation point of the controlled oscillator was numerically demonstrated. Kimiaeifar et al. ³² employed the homotopy analysis method to analyze DVdP oscillator. A comparison between the analytical solutions and the numerical results has been proven effective and convenient. Therefore, the method is an effective implementation for resolving this kind of nonlinear problems.

Amer and Soleman ³³ examined the reaction of a dynamical system of a parametric excited pendulum. The delayed feedback control was applied to overcome the vibration of the system. The periodic solution of a nonlinear oscillator describing the generalized Rayleigh equation was investigated by Cveticanin et al. ³⁴ The obtained analytical solutions matched those calculated by the elliptic harmonic balance with generalized Fourier series and Jacobian elliptic functions. Miwadinou et al. ³⁵ examined the nonlinear dynamics of the hybrid Rayleigh–Van der Pol–Duffing oscillator (HRVD). They revealed that the system presents nine resonance states. They showed that the resonance phenomena were strong with respect to the nonlinear quadratic and cubic damping as well as the external force. Additionally, the numerical simulations were utilized to illustrate bifurcation. Kumar et al. ³⁶ proposed a single-degree-of-freedom oscillator in order to characterize the side force acting on a rigid ground owing to human walking. The proposed oscillator was a modification of the HRVD oscillator by presenting a supplementary nonlinear desensitization term. Stability analysis of the suggested oscillator has been performed by consuming the energy balance method and performance Poincaré–Lindstedt (P-L) technique. Amer et al. ³⁷ presented the nonlinear integrated positive position feedback (NIPPF) approach to combine the compensations of both integral resonant controller (IRC) and positive position feedback controller (PPF) of control nonlinear systems. They utilized the equation of the HRVD oscillator by combining NIPPF to control the vibrating system.

There are numerous ways to reduce the higher amplitude of vibrations that occur in nonlinear dynamical systems. In this regard, Amer et al. ³⁸ compared the evolution of time in three different types of control, namely, IRC, PPF, and NIPPF on their main model. They showed that the best one of them is the NIPPF controller, which reduces the vibration at larger or shorter times. Furthermore, as shown by Omidi and Mahmoodi, ³⁹ the NIPPF control has a great influence because it reduces the oscillation at the exact resonance frequency. Kandil and El-Gohary, ⁴⁰ regarding Lyapunov method, applied the proportional derivative control to the delay time effect of its performance. This effect was reduced with the reduction of the rotational beam oscillations at varying speed. The multiple scales method was utilized to obtain the analytical solution. The Lyapunov first technique was employed to carry out a stability examination for plotting the bifurcation profiles. The equation of the HRVD subjected to an external force was presented by Cveticanin et al. ³⁴ as follows

x ¨ + ω 0 2 x + 2 μ x ˙ + β 1 x x ˙ + β 2 x ˙ 2 + δ 1 x 2 x ˙ + δ 2 x ˙ 3 + λ x 3 + γ x 5 = F cos ϖ t

(1)

For simplicity, the coefficients that appear in equation (1) may be listed as follows:

OPEN IN VIEWER

No	Symbol	Definition
1	ω 0	Natural frequency
2	μ	Linear damping coefficient
3	β 1	Impure quadratic damping coefficient
4	β 2	Pure quadratic damping coefficient
5	δ 1	Impure cubic damping coefficient
6	δ 2	Pure cubic damping coefficient
7	λ	Cubic nonlinear Duffing coefficient
8	γ	Quintic nonlinear Duffing coefficient
9	F	External force
10	ϖ	External forcing frequency

Constructed with the aforementioned features, the current study aims to examine the NIPPF controller under an external force as the HRVD oscillator. The vibrating system is calculated analytically by using the frequency response equation. Furthermore, to determine the optimal conditions for the operation of the system, frequency response curves are plotted for different values ​of the controller and system parameters. The validations of time history figures and frequency response curves (FRCs) for the analytical and numerical results are satisfied. To crystallize the presentation of the current work, the rest of the article is structured as follows: Linearized Stability of the HRVD is devoted to adapting the Lindstedt–Poincaré technique in order to extract an analytic bounded solution of the HVPD oscillator. Controller Design via Multiple Scales Method examines the linearized stability of the considered problem. This section involves the phase portraits around the equilibrium points of the HRVD oscillator. The controller design that is proposed by using the multiple scales method and its subsections is introduced in Steady-State Solution. Finally, based on the fundamental findings, the concluding remarks are given throughout Conclusions.

Poincaré–Lindstedt technique

In perturbation theory, the Lindstedt–Poincaré technique (L-P) ²⁷ is a procedure for approximating consistent periodic solutions to differential equations which the classical perturbation approaches fail to provide. In this section, developing the L-P yields the bounded periodic solution of equation (1). In this approach, equation (1) may be written in light of the HPM sense as follows

x ¨ + ω 0 2 x + ρ ( 2 μ x ˙ + β 1 x x ˙ + β 2 x ˙ 2 + δ 1 x 2 x ˙ + δ 2 x ˙ 3 + λ x 3 + γ x 5 − F cos ϖ t ) = 0 ; ρ ∈ [ 0,1 ]

(2)

Along with the current methodology, the time-independent variable t will be transformed to another time-independent variable τ like τ = ω t , where ω is defined as an artificial frequency of the given oscillator.

It follows that equation (2) may be transferred to the following equation

ω 2 x ″ + ω 0 2 x + ρ ( 2 μ ω x ′ + β 1 ω x x ′ + β 2 ω 2 x ′ 2 + δ 1 ω x 2 x ′ + δ 2 ω 3 x ′ 3 + λ x 3 + γ x 5 − F cos ϖ τ − F ρ Ω 1 ⁡ sin ϖ τ ) = 0

(3)

In order to obtain acceptable secular terms, the initial conditions (I. C.) may be modified as follows

x ( 0 ) = ∑ j = 0 ∞ ρ j a j , and x ˙ ( 0 ) = 0

(4)

According to the methodology of the P-L, the function x and parameter ω are expanded in powers of ρ as follows

x ( t ) = ∑ j = 0 ∞ ρ j x j ( t ) = x 0 ( t ) + ρ x 1 ( t ) + ρ 2 x 2 ( t ) + .... , and, ω = 1 + ρ Ω 1 + ρ 2 Ω 2 + ........

(5)

Substituting equations (5) into (3) and identifying the coefficients of the same powers of ρ on both sides, one finds the following hierarchy equations

ρ 0 : x 0 ( τ ) = a 0 ⁡ cos ω 0 τ

(6)

ρ : x 1 ( τ ) = L T − 1 [ S a 1 S 2 + ω 0 2 ] − L T − 1 { 1 s 2 + ω 0 2 L T { 2 Ω 1 x 0 ″ + 2 μ x 0 ′ + β 1 x 0 x 0 ′ + β 2 x ′ 0 2 + δ 1 x 0 2 x 0 ′ + δ 2 x ′ 0 3 + λ x 0 3 + γ x 0 5 − F ⁡ cos ϖ τ } }

(7)

ρ 2 : x 2 ( τ ) = L T − 1 [ S a 2 S 2 + ω 0 2 ] − L T − 1 { 1 s 2 + ω 0 2 L T { 2 Ω 1 x 1 ″ + 2 Ω 1 μ x 0 ′ + 2 μ x 1 ′ + 2 Ω 1 β 2 x ′ 0 2 + 3 Ω 1 δ 2 x ′ 0 3 + 5 γ x 1 x 0 4 + Ω 1 β 1 x 0 x 0 ′ + β 1 x 0 x 1 ′ + 2 δ 1 x 0 x 1 x 0 ′ + β 1 x 1 x 0 ′ + 2 β 2 x 0 ′ x 1 ′ + 3 δ 2 x ′ 0 2 x 1 ′ + 3 λ x 1 x 0 2 + Ω 1 δ 1 x 0 ′ x 0 2 + δ 1 x 1 ′ x 0 2 + Ω 1 2 x 0 ″ + 2 Ω 2 x 0 ″ − F ϖ Ω 1 ⁡ sin ϖ τ } }

(8)

As identified by equations (7) and (8), the solution of first and second orders depends fundamentally on the zero-order solution.

On substituting equations (6) into (7), by using the Mathematica software version 12.2.0.0, the uniform valid expansion requires a cancellation of the substance of the secular terms. Therefore, the coefficients of the functions cos ω 0 τ and sin ω 0 τ must be canceled. Consequently, the parameters a 0 and Ω 1 are determined as follows

a 0 = 2 − 2 μ δ 1 + 3 δ 2 ω 0 2 and, Ω 1 = μ ( 20 γ μ − 3 λ δ 1 − 9 λ δ 2 ω 0 2 ) ω 0 2 ( δ 1 + 3 δ 2 ω 0 2 ) 2

(9)

It follows that the periodic solution x 1 ( τ ) then becomes

x 1 ( τ ) = 1 2 a 0 2 β 2 − F ω 0 2 − ϖ 2 cos ϖ τ + ( a 1 − 2 3 a 0 2 β 2 + 2 λ a 0 3 32 ω 0 2 + 2 γ a 0 5 24 ω 0 2 + F ω 0 2 − ϖ 2 ) cos ω 0 τ + a 0 2 ( − 32 β 1 − 9 δ 1 a 0 + 9 ω 0 2 δ 2 a 0 ) 96 ω 0 sin ω 0 τ + 1 6 a 0 2 β 2 ⁡ cos ( 2 ω 0 τ ) + β 1 a 0 2 6 ω 0 sin ( 2 ω 0 τ ) − a 0 3 ( 4 λ + 5 γ a 0 2 ) 128 ω 0 2 cos ( 3 ω 0 τ ) + a 0 3 ( δ 1 − δ 2 ω 0 2 ) 32 ω 0 sin ( 3 ω 0 τ ) − γ a 0 5 384 ω 0 2 cos ( 5 ω 0 τ )

(10)

Once more, substituting equations (6), (9), and (10) into (8), the invalidation of the secular terms makes the parameters a 1 and Ω 2 to be determined as follows

a 1 = 1 192 ω 0 2 ( ω 0 2 − ϖ 2 ) ( 3 a 0 2 ( 3 δ 2 ω 0 2 + δ 1 ) + 8 μ ) ( 1536 F μ ω 0 2 ( − 1 + Ω 1 a 0 ( ϖ 2 − ω 0 2 ) ) + 160 a 0 6 β 1 γ ( ϖ 2 − ω 0 2 ) − 8 a 0 5 ( ϖ 2 − ω 0 2 ) ( 8 γ μ + 9 δ 2 λ ω 0 2 + 6 δ 1 λ ) + 39 a 0 7 γ ( ϖ 2 − ω 0 2 ) ( 3 δ 2 ω 0 2 + δ 1 ) + 48 a 0 3 ( ϖ 2 − ω 0 2 ) ( ω 0 2 ( 4 β 1 β 2 + Ω 1 ( 39 δ 2 ω 0 2 + δ 1 ) ) + λ μ ) + 192 a 0 4 ( ϖ 2 − ω 0 2 ) ( β 1 λ − 2 β 2 ω 0 2 ( 3 δ 2 ω 0 2 + δ 1 ) ) − 64 a 0 2 ω 0 2 ( − 8 ϖ 2 ( ω 1 β 1 + 2 μ β 2 ) − 9 F δ 1 + ω 0 2 ( 8 β 1 Ω 1 + 16 μ β 2 − 27 F δ 2 ) )

(11)

Ω 2 = − 1 768 ω 0 4 ( 3 a 0 2 ( 3 δ 2 ω 0 2 + δ 1 ) + 8 μ ) ( 1920 a 0 7 β 1 γ λ + 1000 a 0 9 β 1 γ 2 + 165 a 0 10 γ 2 ( 3 δ 2 ω 0 2 + δ 1 ) + 5 a 0 8 γ ( 38 γ μ + 117 δ 2 λ ω 0 2 + 69 δ 1 λ ) + 6 a 0 6 ( 40 γ ( 5 β 1 β 2 ω 0 2 + λ μ ) + 40 γ ω 0 2 Ω 1 ( 48 δ 2 ω 0 2 + δ 1 ) + 3 ( 9 δ 1 λ 2 + ω 0 2 δ 1 3 ) + 27 δ 1 λ 2 − 81 δ 2 3 ω 0 8 + 45 δ 1 δ 2 2 ω 0 6 + 33 δ 1 2 δ 2 ω 0 4 ) − 3072 μ ω 0 4 Ω 1 2 + 1024 a 0 3 β 1 ω 0 2 ( 3 δ 2 μ ω 0 2 + δ 1 μ − 3 λ Ω 1 ) + 2048 a 0 β 1 ω 0 2 ( μ 2 + ω 0 2 Ω 1 2 ) − 64 a 0 2 ω 0 2 ( μ ( ω 0 2 ( 16 β 2 2 + 9 δ 2 μ ) + 4 β 1 2 − 9 δ 1 μ ) + 12 Ω 1 ( β 1 β 2 ω 0 2 − 9 λ μ ) + 3 ω 0 2 Ω 1 2 ( 21 δ 2 ω 0 2 − 5 δ 1 ) ) + 96 a 0 5 β 1 ( ω 0 2 ( ( 3 δ 2 ω 0 2 + δ 1 ) 2 − 40 γ Ω 1 ) + 9 λ 2 ) − 24 a 0 4 ( − 3 λ 2 μ + ω 0 2 ( − 36 β 1 β 2 λ − 11 δ 1 2 μ + 3 δ 2 ω 0 4 ( 16 β 2 2 + 15 δ 2 μ ) + 2 δ 1 ( ω 0 2 ( 8 β 2 2 − 17 δ 2 μ ) − 2 λ Ω 1 ) + 4 β 1 2 ( 3 δ 2 ω 0 2 + δ 1 ) − 4 Ω 1 ( 100 γ μ + 87 δ 2 λ ω 0 2 ) ) ) )

(12)

It follows that the uniform valid solution x 2 ( τ ) becomes

x 2 ( τ ) = A 0 + A 1 ⁡ cos ω 0 τ + A 2 sin ω 0 τ + A 3 ⁡ cos ( 2 ω 0 τ ) + A 4 ⁡ sin ( 2 ω 0 τ ) + A 5 ⁡ cos ( 3 ω 0 τ ) + A 6 ⁡ sin ( 3 ω 0 τ ) + A 7 ⁡ cos ( 4 ω 0 τ ) + A 8 ⁡ sin ( 4 ω 0 τ ) + A 9 ⁡ cos ( 5 ω 0 τ ) + A 10 ⁡ sin ( 5 ω 0 τ ) + A 11 ⁡ cos ( 6 ω 0 τ ) + A 12 ⁡ sin ( 6 ω 0 τ ) + A 13 ⁡ cos ( 7 ω 0 τ ) + A 14 ⁡ sin ( 7 ω 0 τ ) + A 15 ⁡ cos ( 9 ω 0 τ ) + A 16 ⁡ cos ( ϖ τ ) + A 17 ⁡ sin ( ϖ τ ) + A 18 ⁡ cos ( ϖ + ω 0 ) τ + A 19 ⁡ sin ( ϖ + ω 0 ) τ + A 20 ⁡ cos ( ϖ − ω 0 ) τ + A 21 ⁡ sin ( ϖ − ω 0 ) τ + A 22 ⁡ cos ( ϖ + 2 ω 0 ) τ + A 23 ⁡ sin ( ϖ + 2 ω 0 ) τ A 24 ⁡ cos ( ϖ − 2 ω 0 ) τ + A 25 ⁡ sin ( ϖ − 2 ω 0 ) τ + A 26 ⁡ cos ( ϖ + 4 ω 0 ) τ + A 27 ⁡ cos ( ϖ − 4 ω 0 ) τ

(13)

where the constants A i , i = 1 , 2 ,..., 27 . In order to follow the study more easily, they will be moved from the study.To obtain the constant a 2 , we need to eliminate the secular terms of the third-order solution. This procedure is very tedious, yet straightforward. For easier follow-up of the study, it will be removed.

Finally, the approximate bounded solution of the equation of motion that is given in equation (1) will be written as follows

x ( τ ) = lim ρ → 1 ( x 0 ( τ ) + ρ x 1 ( τ ) + ρ 2 x 2 ( τ ) )

(14)

To check the feasibility of the implication of the previous L-P, it is appropriate to compare this method with the numerical methodology known as RK4. Therefore, in what follows the analytic approximate solution as given by equation (14) is sketched in red. On the other hand, the RK4 of the considered system as given by equation (1) is plotted in blue. The following figure is graphed of a system having the following particulars:

ω 0 = 1 , μ = − 0.0005 , β 1 = 3.5 , β 2 = 0.125 , δ 1 = 1.05 , δ 2 = 0.85 , λ = − 0.6 , γ = 0.05 , F = 0.5 and ϖ = 200

According to these particulars, the first initial condition then becomes: x ( 0 ) = 0.0 3147 . This value is very essential as an initial condition for the implication of the software RK4.

As shown from this figure, the two curves are nearly coincident. This shows that L-P, as an analytic approximate solution, is a promising and powerful perturbed technique Figures 1–7.OPEN IN VIEWER

OPEN IN VIEWER

Figure 2.

Dynamical behavior (stable center) with given parameters listed in Table 1.

OPEN IN VIEWER

Figure 3.

Dynamical behavior (unstable saddle point) with given parameters listed in Table 1.

OPEN IN VIEWER

Figure 4.

Dynamical behavior (unstable saddle point) with given parameters listed in Table 1.

OPEN IN VIEWER

Figure 5.

Dynamical behavior (unstable saddle point) with given parameters listed in Table 1.

OPEN IN VIEWER

Figure 6.

Dynamical behavior (unstable saddle point) with given parameters listed in Table 1.

OPEN IN VIEWER

Figure 7.

Dynamical behavior (unstable saddle point) with given parameters listed in Table 1.

Linearized Stability of the HRVD

Throughout this section, the linearized technique is employed for the HRVD as given by equation (1). Unfortunately, the following procedure is valid only for the autonomous system. Therefore, an external force must be neglected, or F → 0 .

Along with this approach, the autonomous system of the HRVD then becomes

x ¨ + ω 0 2 x + 2 μ x ˙ + β 1 x x ˙ + β 2 x ˙ 2 + δ 1 x 2 x ˙ + δ 2 x ˙ 3 + λ x 3 + γ x 5 = 0

(15)

Considering the transformation: x ˙ = y , it follows that equation (15) may be converted to the following system

x ˙ = g ( x , y ) , y ˙ = h ( x , y )

(16)

g ( x , y ) = y , h ( x , y ) = − ω 0 2 x − 2 μ x ˙ − β 1 x x ˙ − β 2 x ˙ 2 − δ 1 x 2 x ˙ − δ 2 x ˙ 3 − λ x 3 − γ x 5

(17)

The fixed points (equilibrium points) happen at the points ( x 0 , y 0 ) , where

g ( x 0 , y 0 ) = 0 , h ( x 0 , y 0 ) = 0

(18)

It follows that

y 0 = 0

(19)

and

ω 0 2 x 0 + 2 μ y 0 + β 1 x 0 y 0 + β 2 y 0 2 + δ 1 x 0 2 y 0 + δ 2 y 0 3 + λ x 0 3 + γ x 0 5 = 0

(20)

It follows that there are five fixed points as follows

y 0 = 0

(21)

x 0 = ± − λ ± λ 2 − 4 γ ω 0 2 2 γ

(22)

In order to obtain real equilibrium points, the following conditions must be satisfied

( i ) λ 2 ≥ 4 γ ω 0 2

(23)

( ii ) − λ ± λ 2 − 4 γ ω 0 2 2 γ ≥ 0

(24)

Expanding the Taylor expansion up to the first order to develop the functions g ( x , y ) and h ( x , y ) round the critical points, one finds the following Jacobian matrix

J = ( 0 1 − ( ω 0 2 + β 1 y 0 + 2 δ 1 x 0 y 0 + 3 λ x 0 2 + 5 γ x 0 4 ) − ( 2 μ + β 1 x 0 + 2 β 2 y 0 + δ 1 x 0 2 + 3 δ 2 y 0 2 ) )

(25)

It follows that at the equilibrium points, the Jacobian matrix becomes

J = ( − Λ 1 − ( ω 0 2 + β 1 y 0 + 2 δ 1 x 0 y 0 + 3 λ x 0 2 + 5 γ x 0 4 ) − ( 2 μ + β 1 x 0 + 2 β 2 y 0 + δ 1 x 0 2 + 3 δ 2 y 0 2 ) − Λ )

(26)

From the above matrix, the eigenvalues are given as

Λ 1,2 = − 1 2 ( 2 μ + β 1 x 0 + 2 β 2 y 0 + δ 1 x 0 2 + 3 δ 2 y 0 2 ± ( 2 μ + β 1 x 0 + 2 β 2 y 0 + δ 1 x 0 2 + 3 δ 2 y 0 2 ) 2 − 4 ( ω 0 2 + β 1 y 0 + 2 δ 1 x 0 y 0 + 3 λ x 0 2 + 5 γ x 0 4 ) ) .

(27)

Characteristically, if the eigenvalue of the Jacobian, which is calculated at the equilibrium point, has a negative real part, the equilibrium point becomes a stable state. The equilibrium point, on the other hand, is unstable if at least one of the eigenvalues has a positive real part. It is more appropriate to consider a sample chosen system to indicate the stability/instability configuration in light of the equilibrium points. Consequently, the nature of the eigenvalues gives the criterion. This procedure may be done in the following Table 1.

OPEN IN VIEWER

Table 1.

Illustrates the different types of the eigenvalues and the corresponding stability/instability.

ω 0 = 1 , μ = 0 , β 1 = 0.1 , β 2 = 0.1 , δ 1 = 1 , δ 2 = 0.1 , λ = 5 , γ = 0.5	( 0.0 )	Λ 1,2 = ± i	A stable center. See Figure 2
ω 0 = 3 , μ = 1 , β 1 = 2 , β 2 = 0.1 , δ 1 = 1 , δ 2 = 4 , λ = − 5 , γ = 0.1	( − 1.37 , 0 )	Λ 1 = − 4.8 , Λ 2 = 3.6	Saddle point (unstable). See Figure 3
	( 1.37 , 0 )	Λ 1 = − 8.61 , Λ 2 = 2.01	Saddle point (unstable). See Figure 3
	( − 6.94 , 0 )	Λ 1,2 = − 18.13 ± i 10.80	Stable focus. See Figure 3
	( 6.94 , 0 )	Λ 1 = − 56.06 , Λ 2 = − 7.94	Stable node. See Figure 3
ω 0 = 0.1 , μ = 1 , β 1 = 0.2 , β 2 = 0.1 , δ 1 = 0.5 , δ 2 = 4 , λ = − 2 , γ = 0.5	( − 0.71 , 0 )	Λ 1 = − 1.999 , Λ 2 = 0.009	Saddle point (unstable). See Figure 4
	( 0.71 , 0 )	Λ 1 = − 2.036 , Λ 2 = 0.0098	Saddle point (unstable). See Figure 4
	( − 1.999 , 0 )	Λ 1,2 = − 1.79 ± i 3.56	Stable focus. See Figure 4
	( 1.999 , 0 )	Λ 1,2 = − 2.19 ± i 3.33	Stable focus. See Figure 4
ω 0 = 0.1 , μ = 2.794 , β 1 = 0.2 , β 2 = 0.1 , δ 1 = 0.5 , δ 2 = 4 , λ = − 2 , γ = 0.5	( − 0.71 , 0 )	Λ 1 = − 5.58 , Λ 2 = 0.0035	Saddle point (unstable) see Figure 5
	( 0.71 , 0 )	Λ 1 = − 5.58 , Λ 2 = 0.0035	Saddle point (unstable). See Figure 5
	( − 1.999 , 0 )	Λ 1,2 = − 3.59 ± i 1.74	Stable focus. See Figure 5
	( 1.999 , 0 )	Λ 1,2 = − 3.99	Stable node. See Figure 5
ω 0 = 0.1 , μ = − 5.1911 , β 1 = 0.2 , β 2 = 0.1 , δ 1 = 0.5 , δ 2 = 4 , λ = − 2 , γ = 0.5	( − 0.71 , 0 )	Λ 1 = − 0.0019 , Λ 2 = 10.39	Saddle point (unstable). See Figure 6
	( 0.71 , 0 )	Λ 1 = − 0.0019 , Λ 2 = 10.36	Saddle point (unstable). See Figure 6
	( − 1.999 , 0 )	Λ 1 = 2.56 , Λ 2 = 6.22	Unstable node. See Figure 6
	( 1.999 , 0 )	Λ 1 = 3.99 , Λ 2 = 3.99	Unstable node. See Figure 6
ω 0 = 1 , μ = − 1 , β 1 = 0.1 , β 2 = 0.1 , δ 1 = 1 , δ 2 = 0.1 , λ = 5 , γ = 0.5	( 0.0 )	Λ 1,2 = 1	Unstable node, see Figure 7

Controller Design via Multiple Scales Method

After adding the NIPPF controller, the equation of HRVD subjected to an external force as shown in Miwadinou et al. ³⁵ may be presented as follows

x ¨ + ω 0 2 x + 2 μ x ˙ + β 1 x x ˙ + β 2 x ˙ 2 + δ 1 x 2 x ˙ + δ 2 x ˙ 3 + λ x 3 + γ x 5 = F ⁡ cos ω t + α 1 y + α 2 z

(28)

y ¨ + ω 1 2 y + 2 μ 1 y ˙ = γ 1 x

(29)

z ˙ + ρ z = γ 2 x

(30)

As shown by Saeed and Kamel, ⁴¹ the previous coefficients may be chosen as follows

μ = ε 2 μ ˜ , μ 1 = ε 2 μ ˜ 1 , F = ε 2 F ˜ , α 1 = ε 2 α ˜ 1 , α 2 = ε 2 α ˜ 2 , γ 1 = ε 2 γ ˜ 1 and γ = ε − 2 γ ˜

Along with Amer et al. ³⁸ and Omidi and Mahmoodi, ³⁹ the NIPPF controller can be represented by using three degrees of freedom as shown in equations (33)–(35). It is worthy to notice that the NIPPF controller includes both types of controllers, PPF and IRC. As seen, PPF controller is represented by equation (34) in such a way that γ 1 is PPF gain and α 1 , which is given in equation (33), is the feedback gain. Additionally, the IRC controller is represented by equation (35) in such a way that γ 2 is the IRC gain, and α 2 sited in equation (33) is the feedback gain. The displacement of the HRVD is x ; meanwhile, y and z are the displacements of the controllers, μ 1 is the linear damping coefficient, ω 1 is the natural frequency of the controller, and ρ is the linear parameter.

Now, RK4 is applied to get the numerical results of the time history and phase portraits of HRVD. In what follows, a sample chosen system of particulars is given by the following parameters

F = 0.965 , ω 0 = 1 , β 1 = 0.5 , β 2 = 0.125 , μ = − 0.0005 , δ 1 = 1.05 , δ 2 = 0.85 , λ = 0.6 , γ = 0.05

The following analysis considers the primary resonance case at which ( ω = ω 0 ) and the controller parameters are

α 1 = 0.2 , α 2 = 0.5 , α 1 = 0.2 , γ 1 = 0.5 , γ 2 = 0.5 , μ 1 = 0.0005 , ρ = ω 0 at ω 1 = ω 0

To show the best control between the three controllers IRC, PPF, and NIPPF, the time history of the amplitude is shown in Figure 8. The suppression of the amplitude from its maximum value that is given in Figure 8(a) may be illustrated in the following Table 2OPEN IN VIEWER

OPEN IN VIEWER

Table 2.

The best control among three controllers of IRC, PPF and NIPPF.

	Figure 8	After time	Reduction amplitude percentage, %
PPF	(b)	200 s	87.99
IRC	(c)	10 s	2.25
NIPPF	(d)	100 s	97.36

On the other hand, the effectiveness of the previous controllers may be given in Table 2. Remember that the effectiveness ratio E a is calculated as: E a (= steady-state amplitude of the system before and divided by after controlling). By using Figure 8(a)–(d), the effectiveness of the various controllers may be listed as in Table 3.

OPEN IN VIEWER

Table 3.

Effectiveness of three controllers.

Controller	Effectiveness
PPF	8.33
IRC	1.02
NIPPF	37.88

This indicates that the effectiveness of the controller NIPPF is the maximum one. Therefore, from Tables 2 and 3, one can say that NIPPF is the best controller to suppress the vibration. Additionally, NIPPF takes a relatively short time to reduce the chaos when compared to both PPF and IRC.

The phase portraits for the main equation and the controllers are plotted in Figure 9. As seen from these figures, the velocity is pictured versus the amplitude. From Figure 9(a) and (c), it is clear that there is a little change between them. By contrast, as in Figure 9(b) and (c), there exists some uniform behavior away from Figure 9(a). Actually, Figure 9(d) gives the finest behavior. Again, this indicates that NIPPF controller is the optimal controller. The obtained results for the time history and phase portrait are in a good agreement with the previous works of Amer et al., ³⁸ Omidi and Mahmoodi, ³⁹ and El-Sayed and Bauomy. ⁴² OPEN IN VIEWER

Returning again to the multiple scales method, ²⁷ along with this approach, one gets the following expansions

x ( t ; ε ) = ε x 1 ( T 0 , T 1 , T 2 ) + ε 2 x 2 ( T 0 , T 1 , T 2 ) + ε 3 x 3 ( T 0 , T 1 , T 2 ) + O ( ε 4 )

(31)

y ( t ; ε ) = ε y 1 ( T 0 , T 1 , T 2 ) + ε 2 y 2 ( T 0 , T 1 , T 2 ) + ε 3 y 3 ( T 0 , T 1 , T 2 ) + O ( ε 4 )

(32)

z ( t ; ε ) = ε z 1 ( T 0 , T 1 , T 2 ) + ε 2 z 2 ( T 0 , T 1 , T 2 ) + ε 3 z 3 ( T 0 , T 1 , T 2 ) + O ( ε 4 )

(33)

The time derivatives can be established as follows

d d t ≡ D 0 + ε D 1 + ε 2 D 2

(34)

d 2 d t 2 ≡ D 0 2 + 2 ε D 0 D 1 + ε 2 ( D 1 2 + 2 D 0 D 2 )

(35)

Substituting equations (31)–(35) into equations (28)–(30), then one equates the same powers for ε in both sides to get the following equations:

Order ( ε )

( D 0 2 + ω 0 2 ) x 1 = 0

(36)

( D 0 2 + ω 1 2 ) y 1 = 0

(37)

( D 0 + ρ ) z 1 = γ 2 x 1

(38)

Order ( ε 2 )

( D 0 2 + ω 0 2 ) x 2 = F ˜ ( e i ω T 0 + e − i ω T 0 ) / 2 − β 1 x 1 D 0 x 1 − β 2 ( D 0 x 1 ) 2 − 2 D 0 D 1 x 1

(39)

( D 0 2 + ω 1 2 ) y 2 = − 2 D 0 D 1 y 1

(40)

( D 0 + ρ ) z 2 = γ 2 x 2 − D 1 z 1

(41)

Order ( ε 3 )

( D 0 2 + ω 0 2 ) x 3 = − λ x 1 3 − γ ˜ x 1 5 + α ˜ 1 y 1 + α ˜ 2 z 1 − β 1 ( x 1 D 1 x 1 + x 2 D 0 x 1 + x 1 D 0 x 2 ) − D 1 2 x 1 − 2 μ ˜ D 0 x 1 − δ 1 x 1 2 D 0 x 1 − 2 β 2 ( D 1 x 1 D 0 x 1 + D 0 x 1 D 0 x 2 ) − δ 2 ( D 0 x 1 ) 3 − 2 D 0 D 2 x 1 − 2 D 0 D 1 x 2

(42)

( D 0 2 + ω 1 2 ) y 3 = − γ ˜ 1 x 1 − D 1 2 y 1 − 2 μ ˜ 1 D 0 y 1 − 2 D 0 D 2 y 1 − 2 D 0 D 1 y 2

(43)

( D 0 + ρ ) z 3 = γ 2 x 3 − D 2 z 1 − D 1 z 2

(44)

The solutions of equations (36) and (37) may be expressed in the following forms

x 1 = A ( T 1 , T 2 ) e i ω 0 T 0 + A ¯ ( T 1 , T 2 ) e − i ω 0 T 0

(45)

y 1 = B ( T 1 , T 2 ) e i ω 1 T 0 + B ¯ ( T 1 , T 2 ) e − i ω 1 T 0

(46)

By using equation (45) and substituting equation (38), one finds

z 1 = γ 2 ( ( ρ − i ω 0 ) A e i ω 0 T 0 + ( ρ + i ω 0 ) A ¯ e − i ω 1 T 0 ) ρ 2 + ω 0 2 + C ( T 1 , T 2 ) e − ρ T 0

(47)

Substituting equations (45)–(47) into (39) and (40), one gets the following equations

( D 0 2 + ω 0 2 ) x 2 = F ˜ e i ω T 0 / 2 − β 2 ω 0 2 A A ¯ + ( − i β 1 ω 0 + β 2 ω 0 2 ) A 2 e 2 i ω 0 T 0 − 2 i ω 0 D 1 A e i ω 0 T 0 + c . c

(48)

( D 0 2 + ω 1 2 ) y 2 = − 2 i ω 1 D 1 B e i ω 1 T 0 + c . c

(49)

The aim of the current work is to suppress the vibration of the considered system. Therefore, the following analysis focuses only on the resonance case. Consequently, before we proceed further, we must determine the resonance cases of equations (48) and (49), which are the primary resonance ( ω = ω 0 ) and internal resonance ( ω 1 = ω 0 ) , respectively. Subsequently, the closeness of the considered resonance cases can be described by introducing the detaining parameters σ 1 and σ 2 as given by Saeed and Kamel ⁴¹

ω = ω 0 + σ 1 = ω 0 + ε σ ˜ 1 , ω 1 = ω 0 + σ 2 = ω 0 + ε σ ˜ 2

(50)

where σ 1 describes the nearness of ω 0 , and σ 2 represents the inherent difference between ω 0 and ω 1 . Substituting equations (50) into (48), one gets

( D 0 2 + ω 0 2 ) x 2 = − β 2 ω 0 2 A A ¯ + ( − i β 1 ω 0 + β 2 ω 0 2 ) A 2 e 2 i ω 0 T 0 + ( F ˜ e i ε σ ˜ 1 T 0 / 2 − 2 i ω 0 D 1 A ) e i ω 0 T 0 + c . c

(51)

The solvability conditions of equations (51) and (49) are

2 i ω 0 D 1 A = F ˜ e i σ ˜ 1 T 1 / 2 , 2 i ω 1 D 1 B = 0

(52)

Therefore, the particular solutions of equations (48) and (49) after eliminating the secular terms are

x 2 = − β 2 A A ¯ − ( − i β 1 + β 2 ω 0 ) A 2 3 ω 0 e 2 i ω 0 T 0 + c . c

(53)

y 2 = 0

(54)

By using equations (47), (52), and (53) then substituting equation (41), one gets

( D 0 + ρ ) z 2 = − β 2 γ 2 ω 0 2 A A ¯ + ( i β 1 − β 2 ω 0 ) γ 2 A 2 3 ω 0 e 2 i ω 0 T 0 + i F ˜ γ 2 ( ρ − i ω 0 ) 4 ω 0 ( ρ 2 + ω 0 2 ) e i σ ˜ 1 T 1 e i ω 0 T 0 − D 1 C ( T 1 , T 2 ) e − ρ T 0 + c . c

(55)

For the purpose of eliminating the secular term in equation (55), the function C ( T 1 , T 2 ) is then independent of T 1 . Therefore, one may write

C ( T 1 , T 2 ) = H ( T 2 )

(56)

The solution z 2 then becomes