Nonlinear Dynamics of an Electromagnetically Actuated Cantilever Beam Under Harmonic External Excitation

Abstract

The present work is devoted to the study of nonlinear vibrations of an electromagnetically actuated cantilever beam subject to harmonic external excitation. The soft actuator that controls the vibratory motion of such components of a robotic structure led to a strongly nonlinear governing differential equation, which was solved in this work by using a highly accurate technique, namely the Optimal Auxiliary Functions Method. Comparisons between the results obtained using our original approach with those of numerical integration show the efficiency and reliability of our procedure, which can be applied to give an explicit analytical approximate solution in two cases: the nonresonant case and the nearly primary resonance. Our technique is effective, simple, easy to use, and very accurate by means of only the first iteration. On the other hand, we present an analysis of the local stability of the model using Routh–Hurwitz criteria and the eigenvalues of the Jacobian matrix. Global stability is analyzed by means of Lyapunov’s direct method and LaSalle’s invariance principle. For the first time, the Lyapunov function depends on the approximate solution obtained using OAFM. Also, Pontryagin’s principle with respect to the control variable is applied in the construction of the Lyapunov function.

1. Introduction

Usually, a solenoid actuator contains housing, a moving beam, coils, and a spring. If voltage is applied through the coils, then a current flow results in an electromagnetic force, causing the beam to move [1]. An inherent nonlinear relationship between displacement and current input exists for solenoids. The study of the nonlinear vibration of electromagnetically actuated beams is one of the thrust areas of current research. An electromagnetic actuator is an important part of a classical control system and is extensively used in the automotive sector, biomedicine, electronics, transportation, military, and some other fields. Much research has been conducted on various nonlinear dynamic phenomena. By means of the Melnikov function, Poincare map, bifurcation diagram, and fuzzy control algorithm, Haghighi and Markazi [2] examined the effect of excitation amplitude on the vibration control of an electromagnetic actuator (EA). The influence of EA on the frequency response of a harmonically excited cantilever beam was explored analytically based on perturbation methods, both numerically and experimentally, by Belhaq et al. [3]. Chaotic behavior and the control of EA were employed by Tusset et al. [4] with the method of multiple scales and two control strategies: the state-dependent Riccati equation and optimal linear feedback control. The robustness of each strategy was tested considering the presence of parametric errors and the measurement of noise control.

Muscia [5] analyzed an engineering solution to manufacture EAs for moving rudders and find stabilizers for military ships. The finite element analysis method is used to identify the most stressed areas of some elements and to find the most convenient solution to join the pole pieces. Abba and Rachek [6] proposed multi-level coupling transient electromagnetic fields and mechanical structural dynamics based on the finite element method. The aim of this study was to predict with better accuracy the wave magnetic force density to obtain the mechanical deformation occurring in EA. The LuGre friction model was considered by Zuo et al. [7] to describe friction during the motion of a novel nonlinear lumped parameter model of EA for the low-speed condition.

Verma et al. [8] investigated a self-sensing EA for the vibration control of flexible structures that do not require any modifications to the design of the EA. It does not require an additional bridge structure to measure the voltage drop across the coil. Wei et al. [9] studied a novel EA that has been designed to produce energy via vibrations while suppressing the said vibrations, which can be installed in many kinds of mechanisms to improve machining accuracy and reduce energy loss. A control algorithm with a position controller and an acceleration controller was applied to the proposed actuator. Zhang et al. [10] considered a novel bistable nonlinear EA with an elastic boundary to enhance the actuation performance when controlled by a harmonic input signal. Phase portraits and Poincare maps were also presented to illustrate how the actuation response of the system evolves.

Seebacher et al. [11] analyzed pseudo-density topology optimization in nonlinear electromagnetism that followed the solid isotropic material with the penalization method. The linear characteristics of ferromagnetic material results were applied in a small or wrong EA dimension. A sliding mode control method was used by Al-Bakr et al. [12] to design the closed-loop coil voltage of a nonlinear EA. The reliability of the control law was approved under a range of steady-state actuator position dispersions and coil voltage disturbances using the Monte-Carlo simulation method. Mansour et al. [13] proposed a novel miniaturized and modular dual-axis EA. The modularity can be used in order to generate a smoke-like robot or to generate a quadruped robot.

Szmidt et al. [14] employed the stabilization of a cantilever pipe discharging fluid using the EA of the transformer type. They introduced the effect of negative stiffness and lateral vibrations in the magnetic circuits of the actuator. The equation of the pipe’s dynamics was discretized by the Galerkin procedure. A new concept of two-way EA with the help of magnetic damping was derived by Prajwal et al. [15], which can extend its arm on both sides to facilitate active suspension mechanisms in both humps and potholes. The neuro-fuzzy control theory was applied by Repinaldo et al. [16] to reduce the amplitude of vibration in the structure of two degrees of freedom using EA when submitted to an impulsive force.

Konig et al. [17] performed the design of a hysteresis-compensated self-sensing algorithm with low computational effort. The Integrator-Based Direct Inductance Measurement technique was applied for the resource-efficient estimation of the incremental inductance. Then, a modified Prandl–Ishlinski approach was used to model hysteretic behavior. Yang et al. [18] investigated the broadband active vibration isolation of a bistable nonlinear EA with an elastic boundary. The input-to-state stability of the control law for any non-negative control weights was proved, and thus, it is considered a model-free control method.

Luo et al. [19] evaluated the chaos control problem of MEMS resonators by using analog circuits. The dynamic analysis based on bifurcation diagrams, phase diagrams, and Lyapunov exponents illustrated how the transient chaos behaviors and chaotic behaviors strongly depend on system parameters and on the initial conditions of MEMS.

A hysteresis modeling method including a magnetic flux leakage effect was used by Zhang et al. [20] to compensate for the inherent hysteresis and magnetic flux leakage or large stroke Maxwell EA. The physical model of the actuator was established by means of Ampere’s loop law, the continuity principle of magnetic flux, and Ohm’s law of magnetic circuits. Lin et al. [21] examined the stochastic vibration responses of a bistable EA with an elastic boundary controlled by Gaussian white noise and low-pass-filtered stochastic noise. To realize the biostability and the elastic boundary, this actuator used an oblique spring and another linear spring.

The objective of this present work is to investigate the nonlinear, damped, and forced oscillations of an electromagnetically actuated cantilever beam in the neighborhood of primary resonance. The Optimal Auxiliary Functions Method was applied to obtain an explicit and very accurate analytical approximate solution in both cases for the governing nonlinear differential equation. The principal novelties of our procedure are the presence of so-called auxiliary functions and the optimal convergence-control parameters, which contribute to the accuracy of the solutions. Also, it is very important to have the construction of the initial and the first iteration. The solutions of each nonlinear differential equation can be reduced to only two components, which are determined from two linear differential equations. The local stability of the proposed model is studied using the Routh–Hurwitz criteria and global stability is analyzed through the Lyapunov theorem and LaSalle’s invariance principle. The Lyapunov function is dependent on the obtained approximate solution. Pontryagin’s principle using the control variable is applied in the construction of the Lyapunov function.

2. The Second Alternative of OAFM

In this section, we consider the general nonlinear differential equation [22,23,24,25,26,27,28,29]:

$$l[z(t)]+n[z(t)]=0.\hspace{1em}t\in D$$

(1)

with the following boundary/initial conditions:

$$b[z(t),{\displaystyle \frac{dz(t)}{dt}}]=0$$

(2)

where l is a linear operator, n is a nonlinear operator, t is the independent variable, z(t) is an unknown function, D is the domain of interest, and b is a boundary operator. We should mention that the linear operator l does not necessarily coincide entirely with the linear part of the governing equation. Henceforward, $\overline{z}(t)$ will be the approximate solution of Equations (1) and (2). Also, we suppose that $\overline{z}(t)$ can be expressed in the form with only the following two components:

$$\overline{z}(t)=z_0(t)+z_1(t,C_1,C_2,\dots ,C_n)$$

(3)

where C₁, C₂, …, C_n are unknown parameters at this moment and their number is arbitrary. The initial approximation z₀(t) and the first approximation z₁(t, C₁, C₂, …,C_n) can be determined as described below. Inserting Equation (3) into Equation (1) produces the following:

$$l[z_0(t)]+l[z_1(t,C_1,C_2,\dots ,C_n)]+n[z_0(t)+z_1(t,C_1,C_2,\dots ,C_n)]=0$$

(4)

It is natural that the initial approximation is determined from the following linear differential equation:

$$l[z_0(t)]=0$$

(5)

with the following boundary conditions:

$$b[z_0(t),{\displaystyle \frac{d{z}_0(t)}{dt}}]=0$$

(6)

It is clear that the linear operator l depends on the boundary/initial conditions (2). The initial approximation z₀(t) is well-determined from the linear differential Equation (5) with the boundary condition (6).

It follows that the first approximation z₁(t, C₁, C₂, …,C_n) is obtained from Equations (4) and (5):

$$l[z_1(t,C_1,C_2,\dots ,C_n)]+n[z_0(t)+z_1(t,C_1,C_2,\dots ,C_n)]=0$$

(7)

with the following boundary/initial conditions:

$$b[z_1(t,C_1,C_2,\dots ,C_n),{\displaystyle \frac{d{z}_1(t,C_1,C_2,\dots ,C_n)}{dt}}]=0$$

(8)

The nonlinear term of Equation (7) is developed in the following form:

$$n[z_0(t)+z_1(t,C_1,C_2,\dots ,C_n)]=n[z_0(t)]+{\displaystyle \sum _{k\ge 1}\frac{z_1^k(t,C_1,C_2,\dots ,C_n)}{k!}}{n}^{(k)}[z_0(t)]$$

(9)

in which k! = 1·2·…·k and n^(k) denotes differentiation of order _k of the nonlinear operator n.

Next, to avoid the difficulties that appear in solving the nonlinear differential Equation (7) and to accelerate the convergence of the approximate solution z₁(t, C₁, C₂, …,C_n) and implicitly of the approximate solution $\overline{z}$, we solve the following equation obtained from (7) and (9):

$$l[z_1(t,C_1,C_2,\dots ,C_n)]+n[z_0(t)]+{\displaystyle \sum _{k\ge 1}\frac{z_1^k(t,C_1,C_2,\dots ,C_n)}{k!}}{n}^{(k)}[z_0(t)]=0$$

(10)

We make the following crucial remarks: in general, the solution of linear differential Equations (5) and (6) can be expressed as follows:

$$z_0(t)={\displaystyle \sum _{i=1}^{n_1}{a}_i{f}_i(t)}$$

(11)

where the coefficients a_i, the functions f_i(t), and the positive integer n₁ are known by elementary calculations. On the other hand, the nonlinear operator n[z₀(t)] from Equation (9) calculated for z₀(t), given by the last equation, may be written as follows:

$$n[z_0(t)]={\displaystyle \sum _{j=1}^{n_2}{b}_j{g}_j(t)}$$

(12)

where the coefficients b_j, the function g_j(t), and the positive integer n₂ are evidently known and depend on the initial approximation z₀(t) and also on the nonlinear operator n.

In what follows, since Equation (10) is very difficult to solve, we did not solve it, but from the theory of differential equations [28], the Cauchy method, the method of influence functions, the operator method, and so on, it was more convenient to consider the unknown first approximation z₁(t, C₁, C₂, …,C_n) depending on z₀(t) and n(z₀(t)). More precisely, we had the freedom to choose the first approximation in the following form:

$$z_1(t,C_1,C_2,\dots ,C_n)={\displaystyle \sum _{i=1}^p{F}_i(C_j,f_k,g_r),}\begin{array}{c}B(z_1,{\displaystyle \frac{d{z}_1}{dt}})=0\end{array}$$

(13)

where F_i are p auxiliary functions depending on the n unknown convergence control parameters C_j and also on the n₁ known functions f_k defined in Equation (11) and on n₂ known functions g_j, which appear in the components of n(z₀(t)) from Equation (12).

In conclusion, the first approximation z₁(t, C₁, C₂, …, C_n) is determined from Equation (13). It follows that the approximate solution of Equations (1) and (2) can be determined from Equations (3), (5) and (13). Finally, the unknown parameters C_j, j = 1, 2, …, n can be optimally identified via rigorous mathematical procedures such as the Ritz method, the Galerkin method, the least squares method, the collocation method, the Kantorowich method, or by minimizing the square residual error. In this way, the optimal values of the convergence-control parameters C_j and the Optimal Auxiliary Functions F_i are known. Furthermore, with these optimal convergence-control parameters known, the approximate solution $\overline{z}(t)$ is well-determined. The accuracy of the results obtained through OAFM is growing along with the increasing number of convergence-control parameters C_i [25,26]. Let us note that the nonlinear differential Equations (1) and (2) are reduced to two linear differential equations that do not depend on all terms of the nonlinear operator n[z₀(t)]. Our original procedure led to very accurate results, was effective and explicit, and provides a rigorous way to control and adjust the convergence of the solutions using only the first iteration, without the presence of any small or large parameters in the governing equations or in the boundary/initial conditions.

3. Mathematical Model of Electrostatically Actuated Cantilever Beam in the Presence of Harmonic Excitation

Consider the single mode motion for a cantilever beam submitted to a harmonic external excitation and to a symmetric electromagnetic actuator presented in Figure 1 [2,3,4]:

$$m{x}^{″}+b{x}^{\prime }+K_{1}x+K_3{x}^3=F\cos \Omega t+F_{em},$$

(14)

where x is the vertical displacement of the cantilever beam, b is the damping coefficient, K₁ and K₃ are the stiffness, F and Ω are the amplitude and frequency of the external force, respectively, and F_em is the electromagnetic force:

$$F_{em}={\displaystyle \frac{1}{2}}{\displaystyle \frac{C_0{V}^2}{{(d-x)}^2}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{C_0{V}^2}{{(d+x)}^2}},$$

(15)

where C₀ is the capacitance of the parallel-plate actuator at rest, d is the initial gap width, and V is the bias voltage. The primes represent the time derivative. By means of the following dimensionless variables:

$$z={\displaystyle \frac{x}{d}};\begin{array}{c}{\omega }_0=\sqrt{{\displaystyle \frac{K_1}{m}}};\end{array}\begin{array}{c}c={\displaystyle \frac{b}{m{\omega }_0}};\end{array}\begin{array}{c}K={\displaystyle \frac{K_3}{m{\omega }_0^2{d}^3}};\end{array}\begin{array}{c}{a}_0={\displaystyle \frac{C_0{V}^2}{2m{\omega }_0^2{d}^3}};\end{array}\begin{array}{c}f={\displaystyle \frac{F}{md{\omega }_0^2}};\end{array}\begin{array}{c}\omega ={\displaystyle \frac{\Omega }{{\omega }_0}};\end{array}\begin{array}{c}\tau ={\omega }_{0}t\end{array}$$

(16)

Equation (14) can be rewritten as follows:

$$\ddot{z}+c\dot{z}+z+k{z}^3=f\cos \omega \tau +a_0\left[{\displaystyle \frac{1}{{(1-z)}^2}}-{\displaystyle \frac{1}{{(1+z)}^2}}\right],\begin{array}{c}{a}_0\ne 0\end{array},$$

(17)

where the dot represents the derivative with respect to τ.

The vibration of the cantilever beam is assumed to have small amplitudes around the equilibrium point z = 0, such that the electromagnetic force can be expressed in the following form:

$$G(z)={\displaystyle \frac{1}{{(1-z)}^2}}-{\displaystyle \frac{1}{{(1+z)}^2}}\approx 4z+8z^3+12z^5+16z^7=F(z),$$

(18)

It is remarkable that the maximum difference of $\left|G(z)-F(z)\right|$ is 1.075 × 10⁻⁵ and, therefore, the function F is a very good approximation of the function G(z), as can be seen in Figure 2.

In this way, the nonlinear differential Equation (17) can be rewritten as follows:

$$\ddot{z}+c\dot{z}+(1-4a_0)z+(k-8a_0)z^3-12a_0{z}^5-16a_0{z}^7=f\cos \omega \tau$$

(19)

The initial conditions for Equation (19) are the following:

$$z(0)=A;\begin{array}{c}\dot{z}(0)=0\end{array},$$

(20)

4. Geometric Nonlinearities and Potential Energy

The nonlinear restoring force can be obtained from Equation (19):

$$F_r={\Omega }_0^{2}z+(k-8a_0)z^3-12a_0{z}^5-16a_0{z}^7,\begin{array}{c}{\Omega }_0^2=1-4a_0,\end{array}\begin{array}{c}{a}_0\ne 0\end{array},$$

(21)

such that the potential energy function can be expressed as follows:

$$V_p={\displaystyle \frac{1}{2}}{\Omega }_0^2{z}^2+{\displaystyle \frac{1}{4}}(k-8a_0)z^4-2a_0{z}^6-2a_0{z}^8,$$

(22)

We considered that the constant of integration is zero.

If η = z², then Equation (22) can be rewritten in the following form:

$$V_p=-2a_0\eta \left({\eta }^3+{\eta }^2-{\displaystyle \frac{k-8a_0}{8a_0}}\eta -{\displaystyle \frac{{\omega }_0^2}{4a_0}}\right),$$

(23)

We considered the following parameters:

$$p={\displaystyle \frac{2}{9}}-{\displaystyle \frac{k}{24a_0}};\begin{array}{c}q={\displaystyle \frac{1}{27}}+{\displaystyle \frac{k-6}{48a_0}};\end{array}\begin{array}{c}D=q^2+p^3\end{array},$$

(24)

and as a result, the potential (23) has two or four equilibrium points, as follows.

If D > 0, then the potential (23) has two equilibrium points (Cardan formula):

$${\eta }_1=0,\begin{array}{c}{\eta }_2=\sqrt[3]{\sqrt{D}-q}\end{array}-\sqrt[3]{\sqrt{D}+q},$$

(25)

such that the potential (22) has the following three equilibrium points:

$$z_1=0,\begin{array}{c}{z}_2=\sqrt{{\eta }_2}-{\displaystyle \frac{1}{3}},\end{array}\begin{array}{c}{z}_3=-\sqrt{{\eta }_2}-{\displaystyle \frac{1}{3}}\end{array}$$

(26)

The first equilibrium point, z₁ = 0, is a stable point, and z₂ and z₃ are unstable points (saddle points).

If D < 0, then the potential (23) has four equilibrium points:

$$\begin{array}{l}{\eta }_1=0,\begin{array}{c}{\eta }_2=\sqrt[3]{\sqrt{D}-q}\end{array}-\sqrt[3]{\sqrt{D}+q},\begin{array}{c}{\eta }_3=\sqrt[3]{\sqrt{D}-q}{\displaystyle \frac{-1+i\sqrt{3}}{2}}-\sqrt[3]{\sqrt{D}+q}{\displaystyle \frac{-1-i\sqrt{3}}{2}}\end{array}\\ {\eta }_4=\sqrt[3]{\sqrt{D}-q}{\displaystyle \frac{-1-i\sqrt{3}}{2}}-\sqrt[3]{\sqrt{D}+q}{\displaystyle \frac{-1+i\sqrt{3}}{2}},\begin{array}{c}{i}^2=-1\end{array}\end{array}$$

(27)

Considering Equation (23) and equation,

$$V_p=0$$

(28)

which follows from Equations (23) and (28) where

$${\eta }_1+{\eta }_2+{\eta }_3=-1,\begin{array}{c}{\eta }_1{\eta }_2{\eta }_3={\displaystyle \frac{{\omega }_0^2}{4a_0}}\end{array}>0$$

(29)

and, therefore, Equation (28) has only one positive solution (η_i > 0) and two negative solutions. In this case, the potential (22) for D < 0 has only three equilibrium points:

$$z_1=0,\begin{array}{c}{z}_2=\sqrt{{\eta }_1},\end{array}\begin{array}{c}{z}_3=-\sqrt{{\eta }_2}\end{array}$$

(30)

In conclusion, the potential (22) always has only three equilibrium points:

$$z_1=0,\begin{array}{c}{z}_2=\sqrt{{\eta }_1}-{\displaystyle \frac{1}{3}},\end{array}\begin{array}{c}{z}_3=-\sqrt{{\eta }_2}-{\displaystyle \frac{1}{3}}\end{array}$$

(31)

where η₁ = η₂ if D > 0, z₁ is a stable point and z₃, z₃ are unstable points.

Figure 3 shows that for a₀ = 0.001, there exists one stable point and two saddle points for different values of parameter K.

Also, the potential energy function (22) increases with the increasing parameter K, and the abscissa of points of intersection between the potential and the horizontal axis increases with the increasing variable z for z > 0. For z < 0, this situation is reversed.

5. OAFM for Electromagnetically Actuated Cantilever Beam

With the help of the following transformation,

$$z(\tau )=A{e}^{-0.5c\tau }y(\tau )$$

(32)

where A is defined in Equation (20), Equations (19) and (20) and can be rewritten as follows:

$$\ddot{y}+p^{2}y+(k-8a_0)A^2{y}^3{e}^{-c\tau }-12a_0{A}^4{y}^5{e}^{-2c\tau }-16a_0{A}^6{e}^{-3c\tau }{y}^7=-{\displaystyle \frac{f}{A}}{e}^{0.5c\tau }\cos \omega \tau$$

(33)

$$y(0)=1;\begin{array}{c}\dot{y}(0)={\displaystyle \frac{1}{2}}c\end{array}$$

(34)

where $p=\sqrt{1-4a_0-0.25c^2}$. The linear and nonlinear operators corresponding to Equation (33) are the following, respectively:

$$l[y(\tau )]=\ddot{y}+p^{2}y;\begin{array}{c}n[\end{array}y(\tau )]=(k-8a_0)A^2{y}^3{e}^{-c\tau }-12a_0{A}^4{y}^5{e}^{-2c\tau }-16a_0{A}^6{e}^{-3c\tau }{y}^7+{\displaystyle \frac{f}{A}}{e}^{0.5c\tau }\cos \omega \tau$$

(35)

The approximate solution of Equations (33) and (34) takes the following form:

$$\overline{y}(\tau )=y_0(\tau )+y_1(\tau ,C_1,C_2,\dots ,C_n)$$

(36)

in which C_i, i = 1, 2, …, n are unknown parameters at this moment. The initial approximation y₀(τ) can be obtained from the linear differential Equations (5) and (6):

$${\ddot{y}}_0+p^2{y}_0=0;\begin{array}{c}{y}_0(0)=1;\end{array}\begin{array}{c}{\dot{y}}_0(0)={\displaystyle \frac{1}{2}}\end{array}c$$

(37)

the solution of which is the following (11):

$$y_0(\tau )=\cos p\tau +{\displaystyle \frac{c}{2p}}\sin p\tau$$

(38)

The nonlinear operator defined by Equation (35), calculated for the solution (38), becomes the following:

$$\begin{array}{l}n[y_0(\tau )]=(k-\gamma a_0)A^2{e}^{-c\tau }\left[\left({\displaystyle \frac{3}{4}}+{\displaystyle \frac{3c^2}{16p^2}}\right)\cos p\tau \right.+\left({\displaystyle \frac{3c}{8p}}+{\displaystyle \frac{3c^3}{32p^3}}\right)\sin p\tau +\left.\left({\displaystyle \frac{1}{4}}-{\displaystyle \frac{3c^2}{16p^2}}\right)\cos 3p+\left({\displaystyle \frac{3c}{8p}}+{\displaystyle \frac{c^3}{32p^3}}\right)\sin 3p\tau \right]-\\ -12a_0{A}^4{e}^{-2c\tau }\left[\left({\displaystyle \frac{5}{8}}+{\displaystyle \frac{5c^2}{10p^2}}+{\displaystyle \frac{5c^4}{128p^4}}\right)\cos p\tau \right.+\left({\displaystyle \frac{5c}{16p}}+{\displaystyle \frac{5c^3}{32p^3}}+{\displaystyle \frac{5c^5}{256p^5}}\right)\sin p\tau +\left({\displaystyle \frac{5}{6}}-{\displaystyle \frac{5c^2}{32p^2}}-{\displaystyle \frac{15c^4}{256p^4}}\right)\cos 3p\tau +\\ +\left({\displaystyle \frac{15c}{320p}}+{\displaystyle \frac{5c^3}{64p^3}}-{\displaystyle \frac{5c^5}{512p^5}}\right)\sin 3p\tau +\left({\displaystyle \frac{1}{16}}-{\displaystyle \frac{5c^2}{32p^2}}+{\displaystyle \frac{c^4}{128p^4}}\right)\cos 5p\tau +\left.\left({\displaystyle \frac{5c}{32p}}-{\displaystyle \frac{5c^3}{64p^3}}-{\displaystyle \frac{c^5}{512p^5}}\right)\sin 5p\tau \right]-\\ -16a_0{A}^6{e}^{-3c\tau }\left[\left({\displaystyle \frac{35}{64}}-{\displaystyle \frac{105c^2}{256p^2}}-{\displaystyle \frac{105c^4}{1024p^4}}+{\displaystyle \frac{35c^6}{4096p^6}}\right)\cos p\tau \right.+\left({\displaystyle \frac{35c}{128p}}+{\displaystyle \frac{105c^3}{512p^3}}+{\displaystyle \frac{105c^5}{2048p^5}}-{\displaystyle \frac{35c^7}{8192p^7}}\right)\sin p\tau +\\ +\left({\displaystyle \frac{21}{64}}-{\displaystyle \frac{105c^2}{256p^2}}-{\displaystyle \frac{105c^4}{1024p^4}}-{\displaystyle \frac{63c^6}{4096p^4}}\right)\cos 3p\tau +\left({\displaystyle \frac{63c}{128p}}-{\displaystyle \frac{105c^3}{512p^3}}+{\displaystyle \frac{21c^5}{2048p^5}}+{\displaystyle \frac{21c^7}{8192p^7}}\right)\sin 3p\tau +\\ +\left({\displaystyle \frac{7}{64}}-{\displaystyle \frac{63c^2}{256p^2}}-{\displaystyle \frac{35c^4}{1024p^4}}+{\displaystyle \frac{7c^6}{4096p^6}}\right)\cos 5p\tau +\left({\displaystyle \frac{35c}{128p}}+{\displaystyle \frac{35c^3}{512p^3}}-{\displaystyle \frac{63c^5}{2048p^5}}+{\displaystyle \frac{7c^7}{8192p^7}}\right)\sin 5p\tau +\\ +\left({\displaystyle \frac{1}{64}}-{\displaystyle \frac{21c^2}{256p^2}}+{\displaystyle \frac{35c^4}{1024p^4}}-{\displaystyle \frac{7c^6}{4096p^6}}\right)\cos 7p\tau +\left({\displaystyle \frac{7c}{128p}}-{\displaystyle \frac{35c^3}{512p^3}}+{\displaystyle \frac{21c^5}{2048p^5}}-{\displaystyle \frac{7c^7}{8192p^7}}\right)\sin 7p\tau \end{array}$$

(39)

From the expressions (38) and (39), it is clear that a lot of possibilities exist to choose the first approximation y₁(τ), taking into account Equation (13), such as the following:

$$y_1(\tau )=C_1{e}^{{\displaystyle \frac{c}{6}}\tau }(3\sin p\tau -\sin 3p\tau )+C_2(5\sin p\tau -\sin 5p\tau )$$

(40)

or

$$y_1(\tau )=e^{-{\displaystyle \frac{c}{2}}\tau }[C_1(\cos p\tau -\cos 3p\tau )+C_2(\cos 3p\tau -\cos 5p\tau )+C_3(\cos p\tau -\cos 5p\tau )+C_4(3\sin p\tau -\sin 3p\tau )]$$

(41)

or

$$\begin{array}{l}{y}_1(\tau )=e^{-c\tau }[C_1(\cos p\tau -\cos 3p\tau )+C_2(3\sin p\tau -\sin 3p\tau )+C_3(\cos p\tau -{\displaystyle \frac{2c}{p}}\sin p\tau -e^{-2c\tau }\cos 3p\tau )+\\ +C_4(\cos 3p\tau -\cos 5p\tau )+C_5(5\sin p\tau -\sin 5p\tau )+C_6({\displaystyle \frac{c}{p}}\sin p\tau -(1-e^{-c\tau })\cos p\tau )]\end{array}$$

(42)

and so on. It is known that there exists a value τ₀ for which the movement is damped in the domain [0, τ₀]. For τ > τ₀, the movement is undamped. Therefore, Equations (40)–(42) are valid on the domain [0, τ₀].

We remark that in the nonresonant case, ω is away from Ω₀ defined in Equation (21), but for the resonant case, we can write that ω is nearly Ω₀: ω = Ω₀ + δε where δ is a detuning parameter and ε << 1.

In the present paper, we consider only expression (40) for the first approximation. It follows that the approximate solutions of Equations (19) and (20) can be obtained from Equations (32), (36), (38), and (40):

$$\overline{z}=A{e}^{-0.5c\tau }[y_0(\tau )+y_1(\tau )]=A{e}^{-{\displaystyle \frac{c}{2}}\tau }(\cos p\tau +{\displaystyle \frac{c}{2p}}\sin p\tau )+C_1{e}^{-{\displaystyle \frac{c}{3}}\tau }(3\sin p\tau -\sin 3p\tau )+C_2{e}^{-{\displaystyle \frac{c}{2}}\tau }(5\sin p\tau -\sin 5p\tau )$$

(43)

Numerical Examples

We illustrate the applicability, accuracy, and effectiveness of OAFM by comparing the analytical approximate solution obtained in this research with the numerical integration results obtained using a fourth-order Runge–Kutta procedure.

For the resonant case ω = Ω₀ + δε, the optimal values of the convergence-control parameters for the domain τ$\in$[0, 90] for k = 2, A = 0.1, c = 0.04, f = 0.001, a₀ = 0.001, δ = 0.02, ε = 0.01 are the following:

$$C_1=0.0174168,C_2=-0.001988$$

(44)

In this domain, the movement is dampened. For τ > 90, the movement is undamped, and the approximate solution (33) becomes the following:

$$\overline{z}=K_3\cos p(\tau -90)+{\displaystyle \frac{K_4}{p}}\sin p(\tau -90)+C_3[3\sin p(\tau -90)-\sin 3p(\tau -90)]+C_4[5\sin p(\tau -90)-\sin 5p(\tau -90)]$$

(45)

where K₃ = $\overline{z}(90)$, K₄ = $\dot{\overline{z}}(90)$. In this case, the optimal values of the convergence-control parameters are C₃ = 0.0018647 and C₄ = −0.000439218.

Finally, the approximate solution (43) on [0, 90] and the approximate solution on [90, 200] can be written as follows, respectively:

$$\begin{array}{l}\overline{z}=0.1e^{-0.02\tau }(\cos p\tau +{\displaystyle \frac{0.02}{p}}\sin p\tau )+0.0174168e^{-0.015\tau }(3\sin p\tau -\sin 3p\tau )-0.00198861(5\sin p\tau -\sin 5p\tau )-\\ +e^{-0.0015\tau }(0.00136514107\cos p\tau -0.022967676\sin p\tau -2.557416\cdot {10}^{-7}\cos 3p\tau +2.98693\cdot {10}^{-5}\sin 3p\tau )\end{array}$$

(46)

$$\begin{array}{l}\overline{z}=K_3\cos p(\tau -90)+{\displaystyle \frac{K_4}{p}}\sin p(\tau -90)+0.0018647[3\sin p(\tau -90)-\sin 3p(\tau -90)]-\\ -0.000439218[5\sin p(\tau -90)-\sin 5p(\tau -90)]\end{array}$$

(47)

Figure 4, Figure 5 and Figure 6 depict the comparisons between the numerical solution of Equations (11) and (20) for the following resonant case and approximate solutions: (46) on [0, 90], (47) on [90, 200], and (46) and (47) on the domain [0, 200], respectively.

It can be seen that the solutions obtained by our procedure are nearly identical to those obtained from the numerical integration method, using only two parameters.

The behavior of the investigated model is very sensitive to the changes in some key parameters, as can be seen in the following figures.

The results in Figure 7 show that a small increase in the parameter a₀ related to the electromagnetic force produces a small decrease in the frequency of the vibration. The damping parameter c has an influence depicted in Figure 8, where one can see that a very small increase in its value produces a significant damping of the vibration. The influence of the parameter k can be observed in Figure 9, where it is shown that a relatively large increase in this parameter simultaneously produces a damping in the motion and an increase in the frequency of vibration.

6. Analysis of the Stability for the Primary Resonance

6.1. Study of Stability of Steady-State Motion

In this section, we reconsider the primary resonance ω = Ω₀ + δε,${\Omega }_0=\sqrt{1-4a_0}$, using the two-variable expansion method [29]. In what follows, we will use the notation ξ, which represents stretched time ωt, and η, which represents slow time εt:

$$\xi =\omega t,\eta =\epsilon t$$

(48)

In order to substitute these transformations into Equation (19), we needed the following expressions:

$$\begin{array}{l}{\displaystyle \frac{dz}{dt}}={\displaystyle \frac{\partial z}{\partial \xi }}{\displaystyle \frac{d\xi }{dt}}+{\displaystyle \frac{\partial z}{\partial \eta }}{\displaystyle \frac{d\eta }{dt}}=\omega {\displaystyle \frac{\partial z}{\partial \xi }}+\epsilon {\displaystyle \frac{\partial z}{\partial \eta }}\\ {\displaystyle \frac{d^{2}z}{d{t}^2}}={\omega }^2{\displaystyle \frac{{\partial }^{2}z}{\partial {\xi }^2}}+2\omega \epsilon {\displaystyle \frac{{\partial }^{2}z}{\partial \xi \partial \eta }}+{\epsilon }^2{\displaystyle \frac{{\partial }^{2}z}{\partial {\eta }^2}}\end{array}$$

(49)

Inserting Equation (49) into Equation (19) and then expanding z in the power series $z(\xi ,\eta )=z_0(\xi ,\eta )+\epsilon z_1(\xi ,\eta )$, we obtained the following partial differential equation:

$$\begin{array}{l}{({\Omega }_0+\delta \epsilon )}^2\left({\displaystyle \frac{{\partial }^2{z}_0}{\partial {\xi }^2}}+\epsilon {\displaystyle \frac{{\partial }^2{z}_1}{\partial {\xi }^2}}\right)+2({\Omega }_0+\delta \epsilon )\epsilon \left({\displaystyle \frac{{\partial }^2{z}_0}{\partial \xi \partial \eta }}+\epsilon {\displaystyle \frac{{\partial }^2{z}_1}{\partial \xi \partial \eta }}\right)+{\epsilon }^2\left({\displaystyle \frac{{\partial }^2{z}_0}{\partial {\xi }^2}}+\epsilon {\displaystyle \frac{{\partial }^2{z}_1}{\partial {\eta }^2}}\right)+\epsilon c\left[({\Omega }_0+\delta \epsilon )\left({\displaystyle \frac{\partial z_0}{\partial \xi }}+\epsilon {\displaystyle \frac{\partial z_1}{\partial \eta }}\right)\right.+\\ +\epsilon \left.\left({\displaystyle \frac{\partial z_0}{\partial \xi }}+\epsilon {\displaystyle \frac{\partial z_1}{\partial \eta }}\right)\right]+(z_0+\epsilon z_1){\Omega }_0^2+\epsilon (k-8a_0)(z_0^3+3z_0^2{z}_1\epsilon )-12a_0\epsilon (z_0^5+5z_0^4{z}_1\epsilon )-16a_0\epsilon (z_0^7+7z_0^6{z}_1\epsilon )-\epsilon f\cos \xi =0\end{array}$$

(50)

Neglecting the terms of 0(ε²), after collecting terms, we can obtain the following:

$${\Omega }_0^2\left({\displaystyle \frac{{\partial }^2{z}_0}{\partial {\xi }^2}}+z_0\right)=0$$

(51)

$${\Omega }_0^2\left[{\displaystyle \frac{{\partial }^2{z}_1}{\partial {\xi }^2}}\right.+z_1+{\displaystyle \frac{2\delta }{{\Omega }_0}}{\displaystyle \frac{{\partial }^2{z}_0}{\partial {\xi }^2}}+{\displaystyle \frac{2}{{\Omega }_0}}{\displaystyle \frac{{\partial }^2{z}_0}{\partial \xi \partial \eta }}+{\displaystyle \frac{c}{{\Omega }_0}}{\displaystyle \frac{\partial z_0}{\partial \xi }}+{\displaystyle \frac{k-\gamma a_0}{{\Omega }_0^2}}{z}_0^3-{\displaystyle \frac{12a_0}{{\Omega }_0^2}}{z}_0^5-{\displaystyle \frac{16a_0}{{\Omega }_0^2}}{z}_0^7-\left.{\displaystyle \frac{f}{{\Omega }_0^2}}\cos \xi \right]=0$$

(52)

The general solution of the linear differential Equation (51) is as follows:

$$z_0(\xi ,\eta )=A(\eta )\cos \xi +B(\eta )\sin \xi$$

(53)

For nonresonant terms in Equation (52), the coefficients of sinξ and cosξ are needed vanish, such that we have the following slow flow:

$$2{\displaystyle \frac{dA}{d\eta }}+cA+2\delta B-{\displaystyle \frac{3(k-8a_0)}{4{\Omega }_0}}B(A^2+B^2)+{\displaystyle \frac{15a_0}{2{\Omega }_0}}B{(A^2+B^2)}^2+{\displaystyle \frac{35a_0}{4{\Omega }_0}}B{(A^2+B^2)}^3=0$$

(54)

$$2{\displaystyle \frac{dB}{d\eta }}+cB-2\delta A+{\displaystyle \frac{3(k-8a_0)}{4{\Omega }_0}}A(A^2+B^2)-{\displaystyle \frac{15a_0}{2{\Omega }_0}}B{(A^2+B^2)}^2-{\displaystyle \frac{35a_0}{4{\Omega }_0}}A{(A^2+B^2)}^3={\displaystyle \frac{f}{{\Omega }_0^2}}$$

(55)

Equilibrium points of the slow flow (54) and (55) correspond to the periodic motion of Equation (19). To determine them, the following is set:

$${\displaystyle \frac{dA}{d\eta }}={\displaystyle \frac{dB}{d\eta }}=0$$

(56)

Using the notation B = Ay from Equations (54)–(56), it holds that

$$A={\displaystyle \frac{fy}{{\Omega }_{0}c(1+y^2)}},B={\displaystyle \frac{f{y}^2}{{\Omega }_{0}c(1+y^2)}}$$

(57)

where y is obtained from the following algebraic expression:

$${\left({\displaystyle \frac{y^2}{1+y^2}}\right)}^3+{\displaystyle \frac{6}{7}}{\left({\displaystyle \frac{{\Omega }_{0}c}{f}}\right)}^2{\left({\displaystyle \frac{y^2}{1+y^2}}\right)}^2-{\displaystyle \frac{3(k-8a_0)}{35{\Omega }_0}}{\left({\displaystyle \frac{{\Omega }_{0}c}{f}}\right)}^4{\displaystyle \frac{y^2}{1+y^2}}+{\displaystyle \frac{4c{\Omega }_0}{35a_{0}y}}{\left({\displaystyle \frac{{\Omega }_{0}c}{f}}\right)}^6+{\displaystyle \frac{8{\Omega }_0}{35a_0}}{\left({\displaystyle \frac{{\Omega }_{0}c}{f}}\right)}^6\delta =0$$

(58)

Figure 10, Figure 11 and Figure 12 plot the variation of y with respect to δ for δ$\in$[−0.4, 0)∪ (0, 0.04] and A and B with respect to δ for c = 0.08, f = 0.0008, a₀ = 0.001, and k = 2.

The stability of the steady-state motion is determined by the eigenvalues of the Jacobian matrix and obtained from Equations (54) and (55):

$$[J]=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]$$

(59)

where

$$a_{11}={\displaystyle \frac{\partial A^{\prime }}{\partial A}}\begin{array}{cc},& a_{12}=\end{array}{\displaystyle \frac{\partial A^{\prime }}{\partial B}}\begin{array}{cc},& a_{21}={\displaystyle \frac{\partial B^{\prime }}{\partial A}}\end{array}\begin{array}{cc},& a_{22}={\displaystyle \frac{\partial B^{\prime }}{\partial B}}\end{array}\begin{array}{cc},& A^{\prime }={\displaystyle \frac{dA}{d\eta }}\end{array}\begin{array}{cc},& B^{\prime }={\displaystyle \frac{dB}{d\eta }}\end{array}$$

(60)

After some manipulations from Equations (54) to (56) and (60), the coefficients a_ij take the following form:

$$\begin{array}{c}{a}_{11}=-{\displaystyle \frac{c}{2}}+{\displaystyle \frac{3(k-8a_0)}{4{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^2{\displaystyle \frac{y^3}{(1+y^2{)}^2}}-15{\displaystyle \frac{a_0}{{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^4{\displaystyle \frac{y^5}{(1+y^2{)}^3}}-{\displaystyle \frac{105}{4}}{\displaystyle \frac{a_0}{{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^6{\displaystyle \frac{y^7}{(1+y^2{)}^4}}\\ a_{12}=-\delta +{\displaystyle \frac{3(k-8a_0)}{4{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^2{\displaystyle \frac{y^2(1+3y^2)}{(1+y^2{)}^2}}-15{\displaystyle \frac{a_0}{{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^4{\displaystyle \frac{y^2(1+5y^2)}{(1+y^2{)}^3}}-{\displaystyle \frac{105}{4}}{\displaystyle \frac{a_0}{{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^6{\displaystyle \frac{y^2(1+7y^2)}{(1+y^2{)}^4}}\\ a_{21}=\delta -{\displaystyle \frac{3(k-8a_0)}{4{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^2{\displaystyle \frac{y^2(3+y^2)}{(1+y^2{)}^2}}+15{\displaystyle \frac{a_0}{{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^4{\displaystyle \frac{y^2(5+y^2)}{(1+y^2{)}^3}}+{\displaystyle \frac{105}{4}}{\displaystyle \frac{a_0}{{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^6{\displaystyle \frac{y^2(7+y^2)}{(1+y^2{)}^4}}\\ a_{22}=-{\displaystyle \frac{c}{2}}-{\displaystyle \frac{3(k-8a_0)}{4{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^2{\displaystyle \frac{y^3}{(1+y^2{)}^2}}+15{\displaystyle \frac{a_0}{{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^4{\displaystyle \frac{y^5}{(1+y^2{)}^3}}+{\displaystyle \frac{105}{4}}{\displaystyle \frac{a_0}{{\Omega }_0}}{\left({\displaystyle \frac{f}{{\Omega }_{0}c}}\right)}^6{\displaystyle \frac{y^7}{(1+y^2{)}^4}}\end{array}$$

(61)

The signs of real parts of the eigenvalues of the Jacobian matrix are obtained from the following characteristic equation:

$$\det \left([J]-\lambda [I_2]\right)=0$$

(62)

where [I₂] is the unity matrix of the second order and λ is the eigenvalue of the Jacobian matrix. Taking into account expression (6), the characteristic Equation (60) becomes the following:

$${\lambda }^2+tr[J]\lambda +\det [J]=0$$

(63)

where the trace of the Jacobian matrix is given by the following:

$$tr[J]=-(a_{11}+a_{22})=c$$

(64)

and the determinant J is the following:

$$\det [J]=a_{11}{a}_{22}-a_{12}{a}_{21}$$

(65)

The determinant of J is presented in Figure 13 for the above data.

The discriminant of Equation (63) is the following:

$$D={(tr[J])}^2-4\det [J]$$

(66)

The signs of the eigenvalues λ₁ and λ₂ are given by

$${\lambda }_1=-{\displaystyle \frac{1}{2}}tr[J]\pm {\displaystyle \frac{1}{2}}\sqrt{D}$$

(67)

to determine local stability. But c > 0, and from Figure 13, it is clear that det[J] > 0 such that the eigenvalues λ₁ and λ₂ are real, negative, or complex with the real part negative, and therefore the steady state is asymptotically stable for the above data.

Since det[J] ≠ 0, there does not exist a saddle-node point, and because tr[J] = c ≠ 0, Hopf bifurcation does not exist.

6.2. Global Stability by Lyapunov Function

The governing equation of an electromagnetically actuated cantilever beam (19) can be written by adding the control input U in the following form:

$$\begin{array}{l}{\dot{z}}_1=z_2\\ {\dot{z}}_2=-c{z}_2-{\Omega }_0^2{z}_1-(k-8a_0)z_1^3+12a_0{z}_1^5+16a_0{z}_1^7+f\cos \omega \tau +U\begin{array}{cc},& \omega ={\Omega }_0+\delta \epsilon \end{array}\end{array}$$

(68)

The tracking errors e₁ and e₂ are as follows:

$$e_1=z_1-\overline{z}\begin{array}{cc},& e_2=z_2-\dot{\overline{z}}+\theta e_1\end{array}$$

(69)

where $\dot{\overline{z}}$ is the approximate analytical solution of Equations (19) and (20) above that obtained by OAFM, θ is a positive parameter, and the control U will be defined later.

To define the errors of the parameters and the control U that appear in Equation (69), we introduced the following parameters:

$${\alpha }_1=k-8a_0\begin{array}{cc};& {\alpha }_2\end{array}=12a_0\begin{array}{cc};& {\alpha }_3\end{array}=16a_0$$

(70)

If $\overline{c}$, ${\overline{\Omega }}^2$, ${\overline{\alpha }}_1$, ${\overline{\alpha }}_2$, ${\overline{\alpha }}_3$ and $\overline{f}$ are so-called estimated parameters [19], then the estimated errors of parameters are defined as follows:

$$\tilde{c}=\overline{c}-c\begin{array}{cc};& {\tilde{\Omega }}_0^2\end{array}={\overline{\Omega }}_0^2-{\Omega }_0^2\begin{array}{cc};& {\tilde{\alpha }}_1={\overline{\alpha }}_1-{\alpha }_1\end{array}\begin{array}{cc};& {\tilde{\alpha }}_2={\overline{\alpha }}_2-{\alpha }_1\end{array}\begin{array}{cc};& {\tilde{\alpha }}_3={\overline{\alpha }}_3-{\alpha }_3\end{array}\begin{array}{cc};& \tilde{f}=\overline{f}-f\end{array}$$

(71)

The Lyapunov function is chosen using the Pontryagin principle in the following form:

$$V(e_1,e_2,\tilde{c},{\tilde{\Omega }}_0^2,{\tilde{\alpha }}_1,{\tilde{\alpha }}_2,{\tilde{\alpha }}_3,\tilde{f})={\displaystyle \frac{1}{2}}({\lambda }_1{e}_1^2+{\lambda }_2{e}_2^2+{\lambda }_3{\tilde{c}}^2+{\lambda }_4{\tilde{\Omega }}_0^2+{\lambda }_5{\tilde{\alpha }}_1^2+{\lambda }_6{\tilde{\alpha }}_2^2+{\lambda }_7{\tilde{\alpha }}_3^2+{\lambda }_8{\tilde{f}}^2$$

(72)

where λ_i, i = 1, 2, …, 8 are positive parameters. The time derivative of the Lyapunov function can be written in the following form, taking into account Equations (69)–(71):

$$\begin{array}{l}{\displaystyle \frac{dV}{d\tau }}={\lambda }_1{e}_1(e_2-\theta e_1)+{\lambda }_2{e}_2[-c{z}_2-{\Omega }_0^2{z}_1-{\alpha }_1{z}_1^3-{\alpha }_2{z}_1^5-{\alpha }_3{z}_1^7+f\cos \omega \tau -\ddot{\overline{z}}+\theta (e_2-\theta e_1)]+\\ +{\lambda }_3\tilde{c}\dot{\tilde{c}}+{\lambda }_4{\tilde{\Omega }}_0{\dot{\tilde{\Omega }}}_0+{\lambda }_5{\tilde{\alpha }}_1{\dot{\tilde{\alpha }}}_1+{\lambda }_6{\tilde{\alpha }}_2{\dot{\tilde{\alpha }}}_2+{\lambda }_7{\tilde{\alpha }}_3{\dot{\tilde{\alpha }}}_3+{\lambda }_8\tilde{f}\dot{\tilde{f}}\end{array}$$

(73)

We can define the input control U through Equation (73) as follows:

$$U=\overline{c}{z}_2+{\overline{\Omega }}_0^2{z}_1+{\overline{\alpha }}_1{z}_1^3+{\overline{\alpha }}_2{z}_1^5+{\overline{\alpha }}_3{z}_1^7-\overline{f}\cos \omega \tau +\ddot{\overline{z}}-\theta e_2$$

(74)

Then, Equation (73) can be rewritten through Equation (71) in the following form:

$$\begin{array}{l}{\displaystyle \frac{dV}{d\tau }}={\lambda }_1{e}_1{e}_2-{\lambda }_1\theta e_1^2+{\lambda }_2(\tilde{c}{z}_2+{\tilde{\Omega }}_0^2{z}_1+{\tilde{\alpha }}_1{z}_1^3+{\tilde{\alpha }}_2{z}_1^5+{\tilde{\alpha }}_3{z}_1^7-\tilde{f}\cos \omega \tau )e_2-{\lambda }_2{\theta }^2{e}_1{e}_2+{\lambda }_3\tilde{c}\dot{\tilde{c}}+{\lambda }_4{\tilde{\Omega }}_0{\dot{\tilde{\Omega }}}_0+\\ +{\lambda }_5{\tilde{\alpha }}_1{\dot{\tilde{\alpha }}}_1+{\lambda }_6{\tilde{\alpha }}_2{\dot{\tilde{\alpha }}}_2+{\lambda }_7{\tilde{\alpha }}_3{\dot{\tilde{\alpha }}}_3+{\lambda }_8\tilde{f}\dot{\tilde{f}}\end{array}$$

(75)

After simple manipulations, Equation (75) becomes the following:

$$\begin{array}{l}{\displaystyle \frac{dV}{d\tau }}=({\lambda }_1-{\displaystyle \frac{1}{2}}{\theta }^2)e_1{e}_2-{\lambda }_1\theta e_1^2+\tilde{c}({\lambda }_2{z}_2{e}_2+{\lambda }_3\dot{\tilde{c}})+{\tilde{\Omega }}_0({\lambda }_2{z}_1{e}_2+{\lambda }_4{\dot{\tilde{\Omega }}}_0)+{\tilde{\alpha }}_1({\lambda }_2{e}_2{z}_1^3+{\lambda }_5{\dot{\tilde{\alpha }}}_1)+\\ +{\tilde{\alpha }}_2({\lambda }_2{e}_2{z}_1^5+{\lambda }_6{\dot{\tilde{\alpha }}}_2)+{\tilde{\alpha }}_3({\lambda }_2{e}_2{z}_1^7+{\lambda }_7{\dot{\tilde{\alpha }}}_3)+\tilde{f}({\lambda }_8\dot{\tilde{f}}-{\lambda }_2{e}_2\cos \omega \tau )\end{array}$$

(76)

The estimated parameters from Equation (76) can be defined as follows:

$$\begin{array}{l}{\displaystyle \frac{d\tilde{c}}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_3}}{z}_2{e}_2\begin{array}{cc};& {\displaystyle \frac{d\tilde{\Omega }}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_4}}{z}_1{e}_2\end{array}\begin{array}{cc};& {\displaystyle \frac{d{\tilde{\alpha }}_1}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_5}}{z}_1^3{e}_2\end{array}\begin{array}{cc};& \begin{array}{cc};& {\displaystyle \frac{d{\tilde{\alpha }}_2}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_6}}{z}_1^5{e}_2\end{array}\end{array}\\ {\displaystyle \frac{d{\tilde{\alpha }}_3}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_7}}{z}_1^7{e}_2\begin{array}{cc};& {\displaystyle \frac{d\tilde{f}}{d\tau }}={\displaystyle \frac{{\lambda }_2}{{\lambda }_8}}{e}_2\cos \omega \tau \end{array}\end{array}$$

(77)

taking into account Equations (68) and (69):

$$\begin{array}{l}{\displaystyle \frac{d\tilde{c}}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_3}}\dot{z}[\dot{z}-\dot{\overline{z}}+\theta (z-\overrightarrow{z})]\begin{array}{cc};& {\displaystyle \frac{d{\tilde{\Omega }}_0}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_4}}z[\dot{z}-\dot{\overline{z}}+\theta (z-\overrightarrow{z})]\end{array}\begin{array}{cc};& {\displaystyle \frac{d{\tilde{\alpha }}_1}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_5}}{z}^3[\dot{z}-\dot{\overline{z}}+\theta (z-\overrightarrow{z})]\end{array}\\ {\displaystyle \frac{d{\tilde{\alpha }}_2}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_6}}{z}^3[\dot{z}-\dot{\overline{z}}+\theta (z-\overrightarrow{z})]\begin{array}{cc};& {\displaystyle \frac{d{\tilde{\alpha }}_3}{d\tau }}=-{\displaystyle \frac{{\lambda }_2}{{\lambda }_7}}{z}^7[\dot{z}-\dot{\overline{z}}+\theta (z-\overrightarrow{z})]\end{array}\begin{array}{cc};& {\displaystyle \frac{d\tilde{f}}{d\tau }}={\displaystyle \frac{{\lambda }_2}{{\lambda }_8}}[\dot{z}-\dot{\overline{z}}+\theta (z-\overrightarrow{z})]\end{array}\cos \omega \tau \end{array}$$

(78)

In this way, Equation (76) becomes the following:

$${\displaystyle \frac{dV}{d\tau }}=({\lambda }_1-{\lambda }_2{\theta }^2)e_1{e}_2-{\lambda }_1\theta e_1^2$$

(79)

The parameter θ defined in Equation (69) is chosen as follows:

$$\theta =\sqrt{{\lambda }_1{\lambda }_2^{-1}}$$

(80)

and, therefore, Equation (79) can be written in the final form as follows:

$${\displaystyle \frac{dV}{d\tau }}=-{\lambda }_1\sqrt{{\lambda }_1{\lambda }_2^{-1}}{e}_1^2$$

(81)

It is clear that dV/dτ < 0.

Using the Lyapunov function and LaSalle’s invariance principle, the system is globally asymptotically stable because V is a positive definite function and dV/dτ is a negative definite.

For the sake of simplicity, we consider λ_i = 1, i = 1, 2, …, 8. For the numerical case, A = 0.1, c = 0.04, k = 2, f = 0.01, a₀ = 0.001, δ = 0.02, ε = 0.01, and $\omega =\sqrt{1-4a_0}+\delta \epsilon$ in Figure 14 and Figure 15, which also depict the errors e₁ and e₂ given by Equation (69) with θ = 1.

Since e₁ and e₂ are known, from Equation (77), the estimation parameters $\tilde{c}$, $\tilde{\Omega }$, ${\tilde{\alpha }}_1$, ${\tilde{\alpha }}_2$, ${\tilde{\alpha }}_3$, and $\tilde{f}$ are known, and therefore, the input control U is known from Equation (74).

6.3. Circuit Model Verification

The term operational amplifier, or op-amp for short, denotes a special type of amplifier that, via the proper selection of its external components, can be configured for a variety of operations, such as amplification, addition, subtraction, differentiation, and integration. The ability to perform mathematical operations was the result of combining high gain with negative feedback.

Most operational amplifiers are powered by a double voltage source with opposite polarities and the usual values being +15 V and −15 V. A double source is obtained by connecting two simple DC sources, V_cc and V_ee, in series (Figure 16).

The inverting amplifier represents one of the most frequently used configurations and has the structure from Figure 17a.

The circuit is in a closed loop because the resistor R_f is connected between the output terminal and the one corresponding to the inverting input. Assuming linear and stable operation, the differential input voltage is forced to zero, and thus, the two input potentials are equal. But the non-inverting input is connected to the ground, so its potential is 0, and the inverting input will also have zero potential of the ground. The fact that the inverting input has a potential equal to that of the mass leads to the conclusion that the input voltage is found entirely at the terminals of resistor R₁. Since the inverting input is a virtual ground point, the output voltage ${\mathrm{V}}_{\mathrm{o}\mathrm{u}\mathrm{t}}$ is equal to the voltage drop across the resistor R_f but is opposite to the feedback voltage and can be written as follows: ${\mathrm{V}}_{\mathrm{o}\mathrm{u}\mathrm{t}}=-{\displaystyle \frac{{\mathrm{R}}_{\mathrm{f}}}{{\mathrm{R}}_1}}{\mathrm{V}}_{\mathrm{i}\mathrm{n}}$, where V_in is the input voltage.

As shown in Figure 17b, the summation of several input signals can be achieved if they are connected to the op-amp inputs, depending on the sign of the input voltage coefficients. By making all resistors equal, the sum (or difference) of the input voltages is obtained.

The derivation and integration circuits are circuits that perform the mathematical operations of integration and derivation. Both circuits change the shape of the processed signal in accordance with the associated mathematical operation. Depending on the placement of the capacitor at the inverting input (c) or in the reaction network (d), the derivative or integrator circuits can be obtained.

The derivation circuit from Figure 17c is the circuit where the input voltage V_in and the output voltage V_out establishes the following relationship: ${\mathrm{V}}_{\mathrm{o}\mathrm{u}\mathrm{t}}=-\mathrm{R}\mathrm{C}{\displaystyle \frac{\mathrm{d}\mathrm{V}\mathrm{i}\mathrm{n}\left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}}}$

The integrator configuration from Figure 17d is the circuit that implements the integration function, and the dependence of the output voltage on the input voltage is the following: ${\mathrm{V}}_{\mathrm{o}\mathrm{u}\mathrm{t}}=-{\displaystyle \frac{1}{\mathrm{R}\mathrm{C}}}{\int }_0^{\mathrm{t}}{\mathrm{V}}_{\mathrm{i}\mathrm{n}}\left(\mathrm{t}\right)\mathrm{d}\mathrm{t}+{\mathrm{V}}_{\mathrm{o}\mathrm{u}\mathrm{t}}\left(0\right)$, where V_out(0) represents the initial value of the output voltage (calculated at time t = 0).

In our approach, for the implementation of the ALPHA function, we designed an electronic circuit composed of AD741 op-amps and passive components, including resistors and capacitors, and an AD633 multiplier with output according to the following catalog specifications: $\mathrm{W}={\displaystyle \frac{({\mathrm{X}}_1-{\mathrm{X}}_2)({\mathrm{Y}}_1-{\mathrm{Y}}_2)}{10}}-\mathrm{Z}$, which is used to generate nonlinear items (where W is the output of the integrated circuit, X₁, X₂, Y₁, and Y₂ are the inputs of the integrated circuit, and the amplifier summing node Z allows the user to add two or more multiplier outputs, convert the output voltage to a current, and configure various analog computational functions; in our case, pin 6 is always connected to the ground). To avoid some confusion, it should be mentioned that Z (input of AD633) from the last relationship described has nothing to do with Z from the following, which represents the movement variable from Equation (17).

The first block of the circuit in Figure 18 has the Z function at the input and raises the third power of the input signal at the output. The first floor implements the square of the Z function divided by 10 and realizes the input and output offset compensation (POT1, R₈, and R₉, similar to the second stage). After the second stage, Z is raised to the third power divided by 100 (and this stage compensates for the input and output offset), and the third stage (non-inverting operational amplifier) divides the input by 100 because at the output Y is used to find Z³.

For the Z¹ case, the first block is completely removed.

The case with Z³ is described in this section.

For the case of Z⁵, an AD633 multiplier is interspersed between the second stage, and the last one, which has the following inputs: at X₁, the output of the second floor (Z³) and at Y₁, the output of the first floor (Z²) (the other inputs are at ground). Thus, after changing the value of the feedback resistor of the last stage (op-amp) to 10 MΩ, the output Y = Z⁵ is obtained.

A similar way is used for case Z⁷.

Figure 19 shows how the functions Z and $\overline{\mathrm{Z}}$ were derived. Each block also contains an inverting op-amp since the first stage has the derived function with negative polarity at its output. Signal summation is performed in the block from Figure 20.

The expression of the final ${\tilde{\mathsf{\alpha }}}_1$ function is implemented with the circuit in the last block (Figure 21). The first stage performs a multiplication and the second circuit in the block implements the integration function with an op-amp.

To validate an approximate solution of Equation (19) that describes an important technical phenomenon, an equivalent circuit model is presented in Figure 16, Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21. Determining it is an essential step for understanding a system or process. Building an electric circuit that implements this equation virtually allows it to be tested and verified without physically constructing a device. The model presented in Figure 16, Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21 can also be further employed in practical interface circuits. Force in the mechanical domain is equivalent to electrical voltage. The mass is represented by an induction, the flexibility is described by a capacitor, and the damping is represented as a resistor. The frequency responses of the vibration are then extracted from the time domain responses, and so on. Simulation can help identify potential problems or discrepancies between mathematical theory and practical reality. Once the equation has been validated by a simulation, the actual construction of the electric circuit allows for testing under real-world conditions. This final step confirms how the theoretical equation describes the behavior of the physical device. This will be performed in future work by the authors. This step-by-step methodology increases confidence in the accuracy of the equation and its potential for application in real systems.

7. Conclusions

The electromagnetic actuator is an important part of classical control systems and has a large domain of applicability in engineering practice. In the present paper, we studied the nonlinear vibrations of an electromagnetically actuated cantilever beam submitted to a harmonic external force. The governing equation contains nonlinear terms of seven orders; hence, this equation is very difficult to solve by exact or approximate methods. Our original approach, namely Optimal Auxiliary Functions Method—OAFM—was applied to solve the present problem. In comparison with any other procedures known in the literature, our technique is based on the presence of some auxiliary functions and optimal convergence-control parameters, which assure the high accuracy of the approximate solutions in comparison with numerical integration results. These parameters are determined using a rigorous mathematical procedure. The initial and the first iterations are determined by an original technique, and the initial nonlinear governing equation is reduced to two linear differential equations. OAFM is applied to obtain an approximate analytical solution in the neighborhood of the primary resonance. This method is easy to apply, producing very accurate results using only the first iteration. By applying OAFM, we reduced the complexity of the problem while maintaining accuracy, which is an advantage over traditional methods like perturbation techniques or the method of multiple scales. In future work, the investigation should be enlarged to other varieties of dynamic conditions, including, for example, Van der Waals and Casimir forces.

The analysis of the local stability of the present problem was made by means of the Routh–Hurwitz criteria. The eigenvalues of the Jacobian matrix depend on the trace of the Jacobian matrix and on the discriminant of the characteristic equation. The Hopf bifurcation and saddle-node points do not exist. Global stability is analyzed by means of Lyapunov’s direct method and LaSalle’s invariant principle. The Lyapunov function depends on the approximate solution obtained through our procedure with the help of a control variable and Pontryagin’s principle.

There is a limitation in this study that must be addressed in future papers: the experimental validation was not developed in this paper, and this is the main direction for future work. Moreover, in future work, the study of chaos related to these devices should be developed.