Dynamic Analysis of a Uniform Microbeam Resting on a Nonlinear Foundation Considering Its Curvature Subjected to a Mechanical Impact and Electromagnetic Actuation

Abstract

This study proposes an investigation into the nonlinear vibration of a simply supported, flexible, uniform microbeam associated with its curvature considering the mechanical impact, the electromagnetic actuation, the nonlinear Winkler–Pasternak foundation, and the longitudinal magnetic field. The governing differential equations and the boundary conditions are modeled within the framework of a Euler–Bernoulli beam considering an element of the length of the beam at rest and using the second-order approximation of the deflected beam and the Galerkin–Bubnov procedure. In this work, we present a novel characterization of the microbeam and a novel method to solve the nonlinear vibration of the microactuator. The resulting equation of this complex problem is studied using the Optimal Homotopy Asymptotic Method, employing some auxiliary functions derived from the terms that appear in the equation of motion. An explicit closed-form analytical solution is proposed, proving that our procedure is a powerful tool for solving a nonlinear problem without the presence of small or large parameters. The presence of some convergence-control parameters assures the rapid convergence of the solutions. These parameters are evaluated using some rigorous mathematical procedures. The present approach is very accurate and easy to implement, even for complicated nonlinear problems. The local stability near the primary resonance is studied.

1. Introduction

The flexural vibration of a microbeam under mechanical impact, electromagnetic actuation, a moving body, or a nonlinear foundation is interesting for basic research on vibration problems that has in view their practical interest. The main component of an electromagnetic actuator is the solenoid actuator. This is comprised of housing, a moving rod, coils, and a spring [1]. The moving rod is made of a ferroelectric material that is surrounded by coils. When a voltage is applied through the coils, a current flow appears through the coils, which results in an electromagnetic force exerted on the moving rod, causing the rod to move. An inherent nonlinear relationship between displacement and current input exists for solenoids.

Some previous studies have been the works of Abiala [2], which used the finite element method and Newmark’s integration to obtain the dynamic response of a beam under uniformly distributed moving loads. The effect of electromagnetic actuation on the dynamics of a single-sided Hertzian contact-forced oscillator is explored near primary resonance and secondary resonance by Bichri et al. [3].

Muscia [4] obtained an engineering solution to manufacture electromagnetic linear actuators for moving rudders and fin stabilizers on military ships. The finite element method analysis helped them to find the most convenient solution for joining the pole pieces to the respective bars. Bichri et al. [5] considered the effects of combined DC and fast AC electromagnetic actuations on the dynamic behavior of an excited cantilever beam considering a single-mode approximation. Also, the influence of the fast AC actuation and the air gap on the nonlinear behavior of the system was analyzed. The time-stepping solution of the electromagnetic field model resulting in the magnetic force densities, computed from the Lorentz force and based on the space–time variation of the induced eddy current, was employed by Abba and Rachek [6]. In this way, this framework was dedicated to predicting the wave magnetic force density for the mechanical deformations occurring in electromagnetic actuators.

A novel nonlinear lumped parameter model of the electromagnetic actuator is proposed by Zuo et al. [7] for the engine’s idling or low-speed condition. The LuGre friction model is applied to describe the friction during the mover’s reciprocating motion. A specific experimental scheme is proposed for validation in which a suspended ring-shaped iron is added as the actuator’s load. Another novel electromagnetic actuator that has been designed to produce energy via vibrations while suppressing said vibrations, which can be installed in many kinds of mechanisms to improve machining accuracy and reduce energy loss, is employed by Wei et al. [8]. Based on this prototype, the characteristics of the electromagnetic actuator are employed by means of measuring experiments. On the other hand, a novel, bistable, nonlinear electromagnetic actuator with an elastic boundary to enhance the actuation performance when controlled by a harmonic input signal is attempted by Zhang et al. [9]. This model has an inclined spring, one end of which is supported by an elastic boundary, which produces bistable nonlinearity to realize the large-amplitude inter-well actuation responses. The governing equations are controlled by a harmonic input signal.

A self-sensing electromagnetic actuator for vibration control of flexible structures is studied by Verma et al. [10]. This scheme does not require any modifications in the design of the electromagnetic actuator and has the advantage that the voltage drop across the coil can be measured in relation to motion. Seebacher et al. [11] introduced a pseudo-density topology optimization method in nonlinear electromagnetism that follows the Solid Isotropic Material with Penalization procedure. The adjoint technique is used to determine the sensitivity of the objective functions. An applied linear characteristic of the ferromagnetic material results in a small or wrongly actuated electromagnetic dimension. The mathematical model of lateral vibration of a cantilever pipe discharging fluid using an electromagnetic actuator of the transformer type in the form of a partial differential equation coupled with two ordinary differential equations is presented by Szmidt et al. [12]. The system is discretized in a set of ordinary differential equations and is solved numerically. The critical velocity and the amplitude of post-critical vibrations are determined and validated experimentally.

Al Bakri et al. [13] considered a sliding mode control method to design the closed-loop coil voltage such that the nonlinear system of an electromagnetic actuator slides along a certain sliding surface. The reliability of the system was approved under a wide range of steady-state actuator position dispersions and coil voltage disturbances using the Monte Carlo simulation method.

A novel miniaturized and modular dual-axis electromagnetic actuator for several applications is proposed by Mansour et al. [14]. The modularity of this actuator gives an extra privilege to connect several modules together, providing more potential applications, and can be used to generate a quadrupled robot. The robot shows a high level of manipulability and has a light weight and a small size. A new concept of two-way electromagnetic actuation with the help of magnetic damping, which can extend its arm on both sides to facilitate active suspension mechanisms in both humps and potholes, has been derived by Prajwal et al. [15]. Catia V5 software was used to design and simulate the actuator model. Repinaldo et al. [16] investigated the parameter identification and application of the active vibration control technique within a structure with two degrees of freedom using an electromagnetic actuator. The neuro-fuzzy control theory was applied to reduce the amplitude of vibration in the structure when subjected to an impulsive force.

Konig et al. [17] explored the design of a hysteresis-compensated self-sensing algorithm with low computational effort and integrator-based direct inductance. The measurement technique was applied for the resource-efficient estimation of the incremental inductance. Also, a modified Prandtl–Ishlinski approach was used for modeling this hysteretic behavior. Yang et al. [18] examined the broadband active vibration of a bistable, nonlinear electromagnetic actuator with an elastic boundary. A new feedback control law was proposed to control the actuator’s input curves to significantly attenuate the broadband vibration transmissibility from the base excitation to the actuator’s movements. The input-to-state stability of the control law for any non-negative control weights was proven.

Zhang et al. [19] proposed the physical model of the actuator via Ampere’s loop law, the continuity principle of magnetic flux, and Ohm’s law of magnetic circuits. To enhance hysteresis model accuracy, the magnetic flux leakage coefficients of the actuator with different positions were obtained. The experimental results demonstrate that the trajectory and tracking errors compared with the feedback control and the unoptimized hysteresis method feedforward control of the optimized hysteresis were substantially reduced. Lim et al. [20] investigated the stochastic vibration responses of a bistable electromagnetic actuator with an elastic boundary controlled by Gaussian white noise and low-pass-filtered stochastic noise. This actuator used an oblique spring and another linear spring to realize the bistability and the elastic boundary, respectively.

Jorshari et al. [21] investigated the vibration of single-walled boron nitride nanotube induced by a moving nanoparticle based on the Rayleigh beam model and used the Pasternak substrate to model the elastic medium. Baroudi et al. [22] predicted quantitatively the motion of a soft-beam clamped Euler’s Elastica by adapting the Hencky technique in a regime of large displacement and deformations. A nonlinear planar beam with stretch and shear deformation is described by Ren et al. [23]. The static equations of equilibria are discretized by the finite difference method. The nonlinear vibration of an elastic string with large amplitude and large curvature has been investigated by Zhao et al. [24]. The modified compeer novel form method and the finite difference scheme are used to calculate the critical parameters, the time history diagram, configuration, total length, or fundamental frequency.

The present study is devoted to the nonlinear forced vibration of simply supported, flexible, uniform microbeams. The nonlinearities of the equations are caused by the curvature of the beam, mechanical impact, the electromagnetic actuator, the nonlinear Winkler–Pasternak elastic foundation, and the longitudinal magnetic field. The governing nonlinear differential equations and the corresponding boundary conditions are established for a Euler–Bernoulli beam considering an element of the length of the beam at rest and applying the second-order approximation of the deflected beam and the Galerkin–Bubnov procedure. The two differential equations are reduced to only one nonlinear differential equation, and this equation is solved by means of the Optimal Homotopy Asymptotic Method (OHAM) with the help of so-called auxiliary functions obtained from the terms that appear in the governing equation. A very accurate solution is obtained using a moderate number of convergence-control parameters.

Our original approach is a very powerful tool for solving a nonlinear problem without the presence of any small parameters in the governing equation or in the boundary conditions. The presence of some convergence-control parameters assures the rapid convergence of the approximate solutions after the first iteration. The convergence-control parameters are evaluated using some rigorous mathematical procedures. In general, we have great freedom to select both auxiliary functions and the number of convergence-control parameters. This approach has proven to be very accurate, simple, and easy to implement for any complicated nonlinear problems. The local stability near the primary resonance is studied using the Homotopy Perturbation Method, appealing to the Jacobian and eigenvalues of the system.

2. Formulation of the Problem

In Figure 1 is represented the physical model of a simply supported microbeam of length L subjected to a mechanical impact by the transversal constant force, an electromagnetic load V_DC, the nonlinear Winkler–Pasternak elastic foundation F_WP, and longitudinal magnetic field F_m.

The transversal and longitudinal displacements are W(x,t) and u(x,t), respectively. The Young’s modulus E, the mass density ρ, and the cross-sectional area A of the microbeam are supposed to be constant.

Considering an element of length dx of the microbeam at rest, the corresponding deflection can be written as

$$ds={\left[{\left(1+{\displaystyle \frac{\partial u}{\partial x}}\right)}^2+{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2\right]}^{1/2}dx$$

(1)

The unit vector parallel to the deflected element is of the form

$$\overline{k}=\left[\left(1+{\displaystyle \frac{\partial u}{\partial x}}\right)\overline{i}+{\displaystyle \frac{\partial W}{\partial x}}\overline{j}\right]{\displaystyle \frac{dx}{ds}}$$

(2)

The length of the microbeam is variable during the motion, such that the tension becomes

$$T=EA{\displaystyle \frac{ds-dx}{dx}}$$

(3)

The governing equations are

$$\mathsf{\rho }A{\displaystyle \frac{{\partial }^{2}u}{\partial t^2}}=-{\displaystyle \frac{\partial }{\partial x}}(T\overline{k})\overline{i}$$

(4)

$$\mathsf{\rho }A{\displaystyle \frac{{\partial }^{2}W}{\partial t^2}}=-{\displaystyle \frac{\partial }{\partial x}}(T\overline{k})\overline{j}+EI{\displaystyle \frac{{\partial }^2}{\partial t^2}}\left({{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}/\sqrt{1+{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2}}^3\right)+F_{mi}+F_{em}+F_{WP}+F_m$$

(5)

where I is the moment of inertia of the microbeam cross-section, F_mi is the impact force, and F_em is the electromagnetic actuation.

Neglecting the terms of the forms ${\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^2$, ${\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}$, ${\displaystyle \frac{{\partial }^{3}u}{\partial x^3}}$, ${\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}$, ${\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}$, and ${\displaystyle \frac{{\partial }^{4}u}{\partial x^4}}$, from Equation (1), one can obtain

$${\displaystyle \frac{dx}{ds}}={\left[{(1+u_x)}^2+W_x^2\right]}^{-1/2}\approx 1-u_x-{\displaystyle \frac{1}{2}}{W}_x^2+{\displaystyle \frac{3}{8}}{W}_x^4$$

(6)

When

$$e={\displaystyle \frac{ds-dx}{ds}}=1-{\displaystyle \frac{dx}{ds}},$$

(7)

then, from Equations (6) and (7), it follows that

$$e={\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2-{\displaystyle \frac{3}{8}}{\left({\displaystyle \frac{{\partial }^{4}W}{\partial x^4}}\right)}^4$$

(8)

Equation (4) can be rewritten, taking into account Equations (1)–(3) and (7), as

$$\rho {\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}+E{\displaystyle \frac{\partial e}{\partial x}}=0$$

(9)

Usually, the longitudinal inertial term ${\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}$ can be neglected, such that, from Equation (9), it is clear that e is a constant, e = C, and, therefore, from Equation (8), it holds that

$${\displaystyle \frac{\partial u}{\partial x}}=C-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2+{\displaystyle \frac{3}{8}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^4$$

(10)

By integrating the last equation with the boundary conditions

$$u(0,t)=u(L,t)=0$$

(11)

the constant C is given by

$$C=e={\displaystyle \frac{1}{2L}}{{\displaystyle {\int }_0^L\left(\frac{\partial W}{\partial x}\right)}}^{2}dx-{\displaystyle \frac{3}{8L}}{{\displaystyle {\int }_0^L\left(\frac{\partial W}{\partial x}\right)}}^{4}dx$$

(12)

Substituting Equation (12) into Equation (5), we obtain

$$\mathsf{\rho }A{\displaystyle \frac{{\partial }^{2}W}{\partial t^2}}+EI{\displaystyle \frac{{\partial }^2}{\partial t^2}}\left({{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}/\sqrt{1+{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2}}^3\right)={\displaystyle \frac{EA}{2L}}{\displaystyle \frac{{\partial }^{2}W}{\partial t^2}}\left[{\displaystyle {\int }_0^L{\left(\frac{\partial W}{\partial x}\right)}^{2}dx-\frac{3}{4}}\right.\left.{\displaystyle {\int }_0^L{\left(\frac{\partial W}{\partial x}\right)}^4}dx\right]+F_{mi}+F_{em}+F_{WP}+F_m$$

(13)

The term F_mi represents the mechanical impact of load F₀:

$$F_{mi}=F_0\delta (x-vt)$$

(14)

where v is the speed of the load, and δ is the Dirac function. The vibration of the microbeam depends on the action of the DC voltage source, such that electromagnetic actuation can be considered as

$$F_{em}={\displaystyle \frac{1}{2}}{\displaystyle \frac{C_0{V}_{DC}^2}{{\left(d-W(x,t)\right)}^2}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{C_0{V}_{DC}^2}{{\left(d+W(x,t)\right)}^2}}$$

(15)

in which C₀ is the capacitance of the actuator, d is the gap width, and V_DC is the voltage.

Expression (15) can be simplified in the following form:

$$\begin{array}{l}{\displaystyle \frac{1}{2}}{\displaystyle \frac{C_0{V}_{DC}^2}{{\left(d-W\right)}^2}}-{\displaystyle \frac{1}{2}}{\displaystyle \frac{C_0{V}_{DC}^2}{{\left(d+W\right)}^2}}={\displaystyle \frac{C_0{V}_{DC}^2}{2d^2}}\left[1/{\left(1-{\displaystyle \frac{W}{d}}\right)}^2-1/{\left(1+{\displaystyle \frac{W}{d}}\right)}^2\right]\approx \\ \approx {\displaystyle \frac{2C_0{V}_{DC}^2}{d^2}}\left[{\displaystyle \frac{W}{d}}+2{\left({\displaystyle \frac{W}{d}}\right)}^3+3.0125{\left({\displaystyle \frac{W}{d}}\right)}^5\right]\end{array}$$

(16)

The nonlinear elastic Winkler–Pasternak and shear foundation of the beam is represented by

$$F_{WP}=K_{1}W+K_3{W}^3-K_s{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}$$

(17)

where K₁ and K₃ are linear and nonlinear coefficients, respectively, and K_s is the shear coefficient of the elastic foundation.

The longitudinal magnetic field is given by magnetic permeability, area A, and the component H_x in the x-direction.

$$F_m=-\mathsf{\eta }A{H}_x^2{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}$$

(18)

According to the governing Equation (13), the term that defines the curvature of the microbeam can be rewritten in the following form:

$${{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}/\sqrt{1+{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2}}^3\approx {\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}\left[1-{\displaystyle \frac{3}{2}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2+{\displaystyle \frac{15}{8}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^4\right]$$

(19)

Inserting Equations (14)–(19) into Equation (13), one retrieves

$$\begin{array}{l}\mathsf{\rho }A{\displaystyle \frac{{\partial }^{2}W}{\partial t^2}}+EI\left[{\displaystyle \frac{{\partial }^{4}W}{\partial x^4}}\right.-{\displaystyle \frac{3}{2}}{\displaystyle \frac{{\partial }^{4}W}{\partial x^4}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2-9{\displaystyle \frac{\partial W}{\partial x}}{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}{\displaystyle \frac{{\partial }^{3}W}{\partial x^3}}-3{\left({\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}\right)}^3+{\displaystyle \frac{45}{2}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2{\left({\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}\right)}^3+\\ +{\displaystyle \frac{45}{2}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^3{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}{\displaystyle \frac{{\partial }^{3}W}{\partial x^3}}+{\displaystyle \frac{15}{8}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^4\left.{\displaystyle \frac{{\partial }^{4}W}{\partial x^4}}\right]-{\displaystyle \frac{EA}{2L}}{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}\left[{\displaystyle {\int }_0^L{\left(\frac{\partial W}{\partial x}\right)}^{2}dx-\frac{3}{4}}\right.\left.{\displaystyle {\int }_0^L{\left(\frac{\partial W}{\partial x}\right)}^4}dx\right]-\\ -{\displaystyle \frac{2C_0{V}_{DC}^2}{d^2}}\left[{\displaystyle \frac{W}{d}}+2{\left({\displaystyle \frac{W}{d}}\right)}^3+3.0125{\left({\displaystyle \frac{W}{d}}\right)}^5\right]-K_{1}W-K_3{W}^3-(K_s+\mathsf{\eta }{H}_x^{2}A){\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}=F_0\delta (x-vt)\end{array}$$

(20)

The above equation is obtained for a simply supported microbeam such that the boundary conditions are

$$W(0,t)=W(L,t)={\displaystyle \frac{{\partial }^{2}W(0,t)}{\partial x^2}}={\displaystyle \frac{{\partial }^{2}W(L,t)}{\partial x^2}}=0$$

(21)

To express the governing Equations (20) and (21) in nondimensional form, the following parameters are defined:

$$\begin{array}{l}\overline{W}={\displaystyle \frac{W}{d}}\begin{array}{cc},& \overline{x}={\displaystyle \frac{x}{L}}\end{array}\begin{array}{cc},& \overline{t}={\displaystyle \frac{t}{L^2}}\sqrt{{\displaystyle \frac{EI}{\mathsf{\rho }A}}}\end{array}\begin{array}{cc},& \overline{v}=vL\sqrt{{\displaystyle \frac{\mathsf{\rho }A}{EI}}}\end{array}\begin{array}{cc},& {\overline{F}}_0={\displaystyle \frac{F_0{L}^5}{dEI}}\end{array},\\ \overline{\alpha }={\displaystyle \frac{d^{2}A}{2I}}\begin{array}{cc},& \overline{\beta }={\displaystyle \frac{d^2}{L^2}}\end{array}\begin{array}{cc},& {\overline{\gamma }}_1={\displaystyle \frac{2C_0{L}^4{V}_{DC}^2}{d^{3}EI}}\end{array}\begin{array}{cc},& {\overline{\gamma }}_2={\displaystyle \frac{L^2}{EI}}(K_s+\mathsf{\eta }A{H}_z^2)\end{array}\begin{array}{cc},& {\overrightarrow{K}}_1={\displaystyle \frac{K_1{L}^4}{EI}}\end{array}\begin{array}{cc},& {\overrightarrow{K}}_3={\displaystyle \frac{d^2{L}^3{K}_3}{EI}}\end{array}\end{array}$$

(22)

Omitting the bar, Equation (20) can be rewritten in nondimensional form as

$$\begin{gathered} {\displaystyle \frac{{\partial }^{2}W}{\partial t^2}}+{\displaystyle \frac{{\partial }^{4}W}{\partial x^4}}-3\beta \left[{\displaystyle \frac{1}{2}}\right.{\displaystyle \frac{{\partial }^{4}W}{\partial x^4}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2+3{\displaystyle \frac{\partial W}{\partial x}}{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}{\displaystyle \frac{{\partial }^{3}W}{\partial x^3}}+\left.{\left({\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}\right)}^3\right]+ \\ {\displaystyle \frac{15}{2}}{\beta }^2\left[3{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^2\right.{\left({\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}\right)}^3+3{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^3{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}{\displaystyle \frac{{\partial }^{3}W}{\partial x^3}}+{\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^4{\displaystyle \frac{{\partial }^{4}W}{\partial x^4}}-\alpha {\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}\left[{\int }_0^1{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^{2}dx-{\displaystyle \frac{3}{4}}\beta {\int }_0^1{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^{4}dx\right] \\ -{\gamma }_1(W+2W^3+3.0925W^5)-K_{1}W-K_3{W}^3+{\gamma }_2{\displaystyle \frac{{\partial }^{2}W}{\partial x^2}}=F_0\delta (x-vt) \end{gathered}$$

(23)

The solution of Equation (23) can be assumed to be of the following form:

$$W(x,t)=X(x)T(t)$$

(24)

where, taking account of the boundary conditions (21), the eigenfunction X(x) can be expressed as

$$X(x)=\sqrt{2}\sin \pi x$$

(25)

which satisfies the orthonormality condition

$${\displaystyle {\int }_0^1{X}^2(x)dx=1}$$

(26)

Substituting Equation (24) into (23), multiplying with X(x), integrating on the domain [0, 1], using the expression

$${\displaystyle {\int }_0^{1}X(x)\delta (x-vz)}dx=X(vt)$$

(27)

and then using the Galerkin–Bubnov procedure, Equation (23) will be written as follows:

$$\ddot{T}+{\mathsf{\omega }}^{2}T+a{T}^3+b{T}^5=f\sin \pi vt$$

(28)

where the dot denotes the derivative with respect to time and

$$\begin{array}{l}{\mathsf{\omega }}^2={\displaystyle {\int }_0^1\frac{d^{4}X}{d{x}^4}Xdx}-(K_1+{\gamma }_1){\displaystyle {\int }_0^1{X}^{2}dx}\begin{array}{cc};& a=-\alpha \left({\displaystyle {\int }_0^1\frac{d^{2}X}{d{x}^2}Xdx}\right.\left({\displaystyle {\int }_0^1{\left(\frac{dX}{dx}\right)}^{2}Xdx}\right)-\end{array}\\ -{\displaystyle \frac{3}{4}}\beta \left.{\displaystyle {\int }_0^1{\left(\frac{dX}{dx}\right)}^{4}Xdx}\right)-3\beta \left[{\displaystyle \frac{1}{2}}\right.{\displaystyle {\int }_0^1\frac{d^{4}X}{d{x}^4}{\left(\frac{dX}{dx}\right)}^{2}Xdx+3}{\displaystyle {\int }_0^1\frac{dX}{dx}\frac{d^{2}X}{d{x}^2}\frac{d^{3}X}{d{x}^3}Xdx}+\\ +\left.{\displaystyle {\int }_0^1{\left(\frac{dX}{dx}\right)}^{3}Xdx}\right]-(K_3+2{\gamma }_1){\displaystyle {\int }_0^1\frac{d^{2}X}{d{x}^2}{X}^{4}dx}\begin{array}{cc};& b={\displaystyle \frac{3}{4}}\alpha \beta \left[{\displaystyle {\int }_0^1\left.\frac{d^{2}X}{d{x}^2}Xdx\right]}\right.\left[{\displaystyle {\int }_0^1{\left(\frac{dX}{dx}\right)}^{4}Xdx}\right]-\end{array}\\ -3.0925{\gamma }_1{\displaystyle {\int }_0^1{X}^{6}dx}+{\displaystyle \frac{15}{2}}{\beta }^2\left[3{\displaystyle {\int }_0^1{\left(\frac{dX}{dx}\right)}^2{\left(\frac{d^{2}X}{d{x}^2}\right)}^{3}Xdx}\right.+\\ +3{\displaystyle {\int }_0^1{\left(\frac{dX}{dx}\right)}^3\frac{d^{2}X}{d{x}^2}\frac{d^{3}X}{d{x}^3}Xdx}+{\displaystyle \frac{1}{4}}{\displaystyle {\int }_0^1{\left(\frac{dX}{dx}\right)}^4\left.\frac{d^{4}X}{d{x}^4}Xdx\right]}\begin{array}{cc};& f=\sqrt{2}{F}_0\end{array}\end{array}$$

(29)

The initial conditions for Equation (28) are

$$T(0)=A\begin{array}{cc},& \dot{T}(0)=0\end{array}$$

(30)

Equation (28) with the initial conditions (30) is a second-order nonlinear forced differential equation of the fifth order and thus is very difficult to analytically solve by an exact solution. In what follows, for Equations (28) and (29), the OHAM is applied to study the nonlinear forced vibration.

It should be mentioned that the coefficients ω, a, b, and f include the effect of the curvature, mechanical impact, electromagnetic actuation, and nonlinear elastic foundation.

3. The Optimal Homotopy Asymptotic Method (OHAM)

The OHAM was proposed for the first time in 2008 by Marinca and Herisanu [25,26,27,28]. To apply the OHAM, we consider the following nonlinear differential equation:

$$L\left(T(t)\right)+N\left(T(t)\right)=0\begin{array}{cc},& t\in D\end{array}\begin{array}{cc},& B\left(T(t),{\displaystyle \frac{dT(t)}{dt}}\right)=0\end{array}$$

(31)

where L[T(t)] is a linear operator, t denotes independent variable, T(t) is an unknown function, N[T(t)] is a nonlinear operator, D is the domain of interest, and B is a boundary operator.

By means of the OHAM, one first constructs a family of equations:

$$(1-p)L[\overline{T}(t,p,C_i)]=H(t,p,C_i)\left[L(t,p,C_i)+N(t,p,C_i)\right]$$

(32)

where p $\in$ [0, 1] is an embedding parameter, H is a nonzero auxiliary function for p $\ne$ 0 and H(t,0,C_i) = 0, $\overline{T}$ is the unknown approximate solution, and C_i, i = 1,2,…,n are n unknown parameters which will be determined later.

Obviously, when p = 0, it holds that

$$L\left(T_0(t)\right)=0\begin{array}{cc},& B\left(T_0(t),{\displaystyle \frac{d{T}_0(t)}{dt}}\right)=0\end{array}$$

(33)

from which can be determined the initial approximation T₀(t). If $\overline{T}\left(t\right)=T_0\left(t\right)+p{T}_1(t,C_1,\dots ,C_n)$, then, from Equation (32), the first approximation T₁(t, C₁, C₂, …, C_n) can be determined for p = 1:

$$L\left(T_1(t,C_1,C_2,\dots ,C_n)\right)=H(t,C_1,C_2,\dots ,C_n)N[T_0(t)]\begin{array}{cc},& B\left(T_1(t),{\displaystyle \frac{d{T}_1(t)}{dt}}\right)\end{array}=0$$

(34)

where H(t, C₁, C₂, …, C_n) is an arbitrary auxiliary function which can be determined such that the products HN[T₀(t)] and N[T₀(t)] are of the same shape. The parameters C₁, C₂,…., C_n which appear in the first-order approximate solution T₁ can be determined in many ways, using the least square method, the Galerkin method, the collocation method, the Ritz method, or by minimizing the square residual error as

$$J(C_1,C_2,\dots ,C_n)={\displaystyle {\int }_0^1{R}^2}(t,C_1,C_2,\dots ,C_n)dt$$

(35)

where R is given by Equation (31):

$$R(t,C_1,C_2,\dots ,C_n)=L\left[\overline{T}(t)\right]+N\left[\overline{T}(t)\right]$$

(36)

The unknown parameter C_i can be identified optimally from the condition

$${\displaystyle \frac{\partial J}{\partial C_1}}={\displaystyle \frac{\partial J}{\partial C_2}}=\dots ={\displaystyle \frac{\partial J}{\partial C_n}}=0$$

(37)

With these parameters known (namely, the convergence-control parameters), the approximate analytical solution T(t) is well determined.

It should be emphasized that the auxiliary function and the number n of the convergence-control parameters are not unique.

The OHAM is a very efficient procedure since any nonlinear differential equation is reduced to only two linear differential equations, and, on the other hand, the convergence-control parameters are optimally determined.

4. Application of Optimal Homotopy Asymptotic Method to Study the Nonlinear Vibration of the Microbeam

In what follows, we apply the OHAM to obtain an approximate solution for Equations (28) and (30). Making the transformations

$$\tau =\mathsf{\Omega }t\begin{array}{cc},& T(t)=A\mathsf{\theta }(\tau )\end{array}$$

(38)

where Ω is the unknown natural frequency of the beam, Equation (28) becomes

$${\mathsf{\theta }}^{″}+{\displaystyle \frac{{\mathsf{\omega }}^2}{{\mathsf{\Omega }}^2}}\mathsf{\theta }+{\displaystyle \frac{a{A}^2}{{\mathsf{\Omega }}^2}}{\mathsf{\theta }}^3+{\displaystyle \frac{b{A}^4}{{\mathsf{\Omega }}^2}}{\mathsf{\theta }}^5-{\displaystyle \frac{f}{A{\mathsf{\Omega }}^2}}\sin \pi {\displaystyle \frac{v\tau }{\mathsf{\Omega }}}=0$$

(39)

and corresponding initial conditions are

$$\begin{array}{cc}& \mathsf{\theta }(0)=1\end{array}\begin{array}{cc},& {\mathsf{\theta }}^{\prime }(0)=0\end{array}$$

(40)

where the prime denotes the derivative with respect to τ.

The linear and nonlinear operators are, respectively,

$$L[\mathsf{\theta }(\tau )]={\mathsf{\theta }}^{″}+\mathsf{\theta }\begin{array}{cc};& N[\mathsf{\theta }(\tau )]=\left({\displaystyle \frac{{\mathsf{\omega }}^2}{{\mathsf{\Omega }}^2}}-1\right)\mathsf{\theta }+{\displaystyle \frac{a{A}^2}{{\mathsf{\Omega }}^2}}{\mathsf{\theta }}^3+{\displaystyle \frac{b{A}^4}{{\mathsf{\Omega }}^2}}{\mathsf{\theta }}^5-{\displaystyle \frac{f}{A{\mathsf{\Omega }}^2}}\sin \pi {\displaystyle \frac{v\tau }{\mathsf{\Omega }}}\end{array}$$

(41)

The initial approximation is determined from Equation (33), which is rewritten as

$${{\mathsf{\theta }}^{″}}_0+{\mathsf{\theta }}_0=0\begin{array}{cc},& {\mathsf{\theta }}_0(0)=1\end{array}\begin{array}{cc},& {{\mathsf{\theta }}^{\prime }}_0(0)=0\end{array}$$

(42)

whose solution is

$${\mathsf{\theta }}_0(\tau )=\cos \tau$$

(43)

The first approximation θ₁(τ) is obtained from Equation (34) as

$${{\mathsf{\theta }}^{″}}_1+{\mathsf{\theta }}_1=H(\tau ,C_1,C_2,\dots ,C_n)N[{\mathsf{\theta }}_0(\tau )]\begin{array}{cc},& {\mathsf{\theta }}_1(0)=\end{array}{{\mathsf{\theta }}^{\prime }}_1(0)=0$$

(44)

where N[θ₀(τ)] is determined by inserting Equation (43) into the nonlinear operator (41):

$$N[{\mathsf{\theta }}_0(\tau )]=A_1\cos \tau +A_3\cos 3\tau +A_5\cos 5\tau -{\displaystyle \frac{f}{A{\mathsf{\Omega }}^2}}\sin {\displaystyle \frac{\pi v\tau }{\mathsf{\Omega }}}$$

(45)

The coefficients from the above equations are

$$A_1={\displaystyle \frac{{\mathsf{\omega }}^2}{{\mathsf{\Omega }}^2}}-1+{\displaystyle \frac{3a{A}^2}{4{\mathsf{\Omega }}^2}}+{\displaystyle \frac{5b{A}^4}{8{\mathsf{\Omega }}^2}}\begin{array}{cc},& A_3={\displaystyle \frac{a{A}^2}{4{\mathsf{\Omega }}^2}}+\end{array}{\displaystyle \frac{5b{A}^4}{16{\mathsf{\Omega }}^2}}\begin{array}{cc},& A_5={\displaystyle \frac{b{A}^4}{16{\mathsf{\Omega }}^2}}\end{array}$$

(46)

The auxiliary function H from Equation (44) is chosen such that the product between HN[θ₀(τ)] and N[θ₀(τ)] is of the same form. For example, we can choose these auxiliary functions in the following expressions:

$$H_1(t,C_1,C_2,C_3,C_4)=C_1+2C_2\cos 2\tau +2C_3\cos 4\tau +2C_4\cos 6\tau$$

(47)

$$H_2(t,C_1,C_2,C_3)=C_1+2C_2\cos 2\tau +2C_3\cos 4\tau$$

(48)

$$H_3(t,C_1,C_2,C_3,C_4,C_5)=C_1+2C_2\cos 2\tau +2C_3\cos 4\tau +2C_4\cos 6\tau +2C_5\cos 8\tau$$

(49)

$$H_4(t,C_1,C_2)=C_1+2C_2\cos 4\tau$$

(50)

Using only expression (47), Equation (44) can be rewritten after some manipulations as follows:

$$\begin{array}{l}{{\mathsf{\theta }}^{″}}_1+{\mathsf{\theta }}_1=[A_1(C_1+C_2)+A_3(C_2+C_3)+A_5(C_3+C_4)]\cos \tau +[A_3(C_1+C_4)+\\ +A_1(C_2+C_3)+A_5{C}_2]\cos 3\tau +[A_5{C}_1+A_3{C}_2+A_1(C_3+C_4)]\cos 5\tau +(A_5{C}_2+\\ +A_3{C}_3+A_1{C}_4)\cos 7\tau +(A_5{C}_3+A_3{C}_4)\cos 9\tau +A_5{C}_4\cos 11\tau \begin{array}{cc},& {\mathsf{\theta }}_1(0)=\end{array}{{\mathsf{\theta }}^{\prime }}_1(0)=0\end{array}$$

(51)

No secular term in Equation (51) requires that the coefficient of cosτ be zero. From this condition, one can determine the natural frequency:

$${\mathsf{\Omega }}^2={\mathsf{\omega }}^2+{\displaystyle \frac{3a{A}^2}{4}}+{\displaystyle \frac{5b{A}^4}{8}}+\left({\displaystyle \frac{a{A}^2}{4}}+{\displaystyle \frac{5b{A}^4}{16}}\right){\displaystyle \frac{C_2+C_3}{C_1+C_2}}+{\displaystyle \frac{b{A}^4}{16}}{\displaystyle \frac{C_3+C_4}{C_1+C_2}}$$

(52)

The solution of linear Equation (51) is

$$\begin{array}{l}{\mathsf{\theta }}_1(\tau )={\displaystyle \frac{A_3(C_1+C_4)+A_1(C_2+C_3)+A_5{C}_2}{8}}(\cos \tau -\cos 3\tau )+\\ +{\displaystyle \frac{A_5{C}_1+A_3{C}_2+A_1(C_3+C_4)}{24}}(\cos \tau -\cos 5\tau )+{\displaystyle \frac{A_5{C}_2+A_3{C}_3+A_1{C}_4}{48}}(\cos \tau -\\ \cos 7\tau )+{\displaystyle \frac{A_5{C}_3+A_3{C}_4}{80}}(\cos \tau -\cos 9\tau )+{\displaystyle \frac{A_5{C}_4}{120}}(\cos \tau -\cos 11\tau )\end{array}$$

(53)

The approximate solution $\overline{T}$ is obtained from Equations (38), (41) and (55):

$$\begin{array}{l}\overline{T}(t)=A\cos \mathsf{\Omega }t+{\displaystyle \frac{A[A_3(C_1+C_4)+A_1(C_2+C_3)+A_5{C}_2]}{8}}(\cos \mathsf{\Omega }t-\cos 3\mathsf{\Omega }t)+\\ +{\displaystyle \frac{A[A_5{C}_1+A_3{C}_2+A_1(C_3+C_4)]}{24}}(\cos \mathsf{\Omega }t-\cos 5\mathsf{\Omega }t)+\\ {\displaystyle \frac{A[A_5{C}_2+A_3{C}_3+A_1{C}_4]}{48}}(\cos \mathsf{\Omega }t-\cos 7\mathsf{\Omega }t)+{\displaystyle \frac{A[A_5{C}_3+A_3{C}_4]}{80}}(\cos \mathsf{\Omega }t-\cos 9\mathsf{\Omega }t)+\\ +{\displaystyle \frac{A{A}_5{C}_4}{120}}(\cos \mathsf{\Omega }t-\cos 11\mathsf{\Omega }t)\end{array}$$

(54)

where Ω is given by Equation (32).

The OHAM can be applied to any nonlinear differential equations without a supplementary hypothesis, such as, for example, the existence of a small parameter in the differential equation or in the boundary/initial conditions.

5. Numerical Results

The efficiency of the OHAM can be proved through numerical investigations for the following case considering α = 0.007, β = 0.001, γ₁ = 0.3, γ₂ = 0.5, k₁ = 1, and k₃ = 0.4. It follows that ω = 1.1, a = 0.023, b = 0.03, f = 0.03, and v = 0.95.

The optimal values of the convergence-control parameters are obtained by minimizing the residual error: C₁ = 0.430116872481, C₂ = 0.000287558193, C₃ = 0.006874588249, C₄ = 0.049982573619, and Ω = 1.116.

To validate the obtained analytical solution, the analytical results are compared with numerical integration results obtained by means of a fourth-order Runge–Kutta approach. The explicit analytical solution benefits from the optimal values identified for the convergence-control parameters and is plotted in Figure 2, demonstrating the high accuracy of the employed method. In Figure 2, we present the approximate solution (54) in comparison with the numerical integration results. It can be observed that the error between the approximate analytical solution obtained by the OHAM and the solution obtained by numerical integration is very small, these solutions being in an excellent agreement. In this way, the results obtained by the OHAM are validated through this comparison, and an excellent accuracy of the analytical solution is proved.

6. Analysis of the Stability of the Steady-State Motion near the Primary Resonance

The procedure applied in this section proposes to distinguish between two time-scales by associating a separate independent variable with each one and using the Homotopy Perturbation Method [29]. We reconsider the primary resonance for Equation (28):

$$\mathsf{\omega }=\overline{\mathsf{\omega }}+\overline{\lambda }p\begin{array}{cc},& \overline{\mathsf{\omega }}=\pi v\end{array}$$

(55)

with λ being a detuning parameter and p $\in$ [0, 1] the embedding parameter. As in Section 3,

$$\xi =\mathsf{\omega }t\begin{array}{cc},& \mathsf{\eta }=pt\end{array}$$

(56)

In accordance with the Homotopy Perturbation Method, for p = 0, any nonlinear differential equation reduces to a linear differential equation, and, for p = 1, one obtains the original equation, such that the original Equation (28) is rewritten as

$$\ddot{T}+{\mathsf{\omega }}^{2}T+p(a{T}^3+b{T}^5-f\sin \overline{\mathsf{\omega }}t)=0$$

(57)

The variable T can be written as a power series of p:

$$T=S_0+p{S}_1+p^2{S}_2+\dots$$

(58)

where S_i are the solutions of the i-th order homotopy equations.

To substitute Equation (56) into (57), we need the following expressions:

$${\displaystyle \frac{dT}{dt}}={\displaystyle \frac{\partial T}{\partial \xi }}{\displaystyle \frac{\partial \xi }{\partial t}}+{\displaystyle \frac{\partial T}{\partial \mathsf{\eta }}}{\displaystyle \frac{\partial \mathsf{\eta }}{\partial t}}=\mathsf{\omega }{\displaystyle \frac{\partial T}{\partial \xi }}+p{\displaystyle \frac{\partial T}{\partial \mathsf{\eta }}}$$

(59)

$$\ddot{T}={\displaystyle \frac{d^{2}T}{d{t}^2}}={\mathsf{\omega }}^2{\displaystyle \frac{d^{2}T}{d{\xi }^2}}+2\mathsf{\omega }p{\displaystyle \frac{d^{2}T}{d\xi \partial \mathsf{\eta }}}+p^2{\displaystyle \frac{d^{2}T}{d{\mathsf{\eta }}^2}}$$

(60)

Substituting Equations (57) and (58) into Equation (60), one can obtain

$$\ddot{T}={\overline{\mathsf{\omega }}}^2{\displaystyle \frac{{\partial }^2{S}_0}{\partial {\xi }^2}}+p\left[{\overline{\mathsf{\omega }}}^2{\displaystyle \frac{{\partial }^2{S}_1}{\partial {\xi }^2}}+2\overline{\mathsf{\omega }}\overline{\lambda }{\displaystyle \frac{{\partial }^2{S}_0}{\partial {\xi }^2}}+2\overline{\mathsf{\omega }}{\displaystyle \frac{{\partial }^2{S}_0}{\partial \xi \partial \mathsf{\eta }}}\right]+\dots$$

(61)

From Equations (57), (58) and (61), the homotopy equations of order pⁱ, i = 0,1 are, respectively,

$$p^0\begin{array}{cc}:& {\overline{\mathsf{\omega }}}^2\left({\displaystyle \frac{{\partial }^2{S}_0}{\partial {\xi }^2}}+S_0\right)\end{array}=0$$

(62)

$$p^1\begin{array}{cc}:& {\overline{\mathsf{\omega }}}^2\left({\displaystyle \frac{{\partial }^2{S}_1}{\partial {\xi }^2}}+S_1\right)\end{array}+2\overline{\lambda }\overline{\mathsf{\omega }}{\displaystyle \frac{{\partial }^2{S}_0}{\partial {\xi }^2}}+2\overline{\mathsf{\omega }}{\displaystyle \frac{d^2{T}_0}{d\xi \partial \mathsf{\eta }}}+2\overline{\lambda }\overline{\mathsf{\omega }}{S}_0+a{S}_0^3+b{S}_0^5-f\sin \xi =0$$

(63)

The solution of the linear Equation (62) with partial derivatives is

$$S_0(\xi ,\mathsf{\eta })=A(\mathsf{\eta })\cos \xi +B(\mathsf{\eta })\sin \xi$$

(64)

where A and B are functions depending on η.

Now, substituting Equation (64) into Equation (63) and avoiding the secular terms yields

$$2\overline{\mathsf{\omega }}{\displaystyle \frac{dB}{d\mathsf{\eta }}}+{\displaystyle \frac{3aA}{4}}(A^2+B^2)+{\displaystyle \frac{5bA}{8}}{(A^2+B^2)}^2=0$$

(65)

$$-2\overline{\mathsf{\omega }}{\displaystyle \frac{dA}{d\mathsf{\eta }}}+{\displaystyle \frac{3aB}{4}}(A^2+B^2)+{\displaystyle \frac{5bB}{8}}{(A^2+B^2)}^2-f=0$$

(66)

The equilibrium points correspond to the periodic solution of Equation (28). To determine them, set the following:

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

(67)

In this way, the equilibrium points can be obtained from the algebraic equations

$$A_e\left[6a+5b(A_e^2+B_e^2)=0\right]$$

(68)

$$\left[6a+5b(A_e^2+B_e^2)\right]B_e(A_e^2+B_e^2)=8f$$

(69)

It is possible only in the case where

$$A_e=0\begin{array}{cc},& 5b{B}_e^5+6a{B}_e^3-8f=0\end{array}$$

(70)

In Figure 3, Figure 4 and Figure 5 is depicted the equilibrium point B_e in respect to the parameter f for different values of the parameters a and b: a = 0.25, b = 0.3 (Figure 3); a = 0.5, b = 1 (Figure 4); a = −1, b = 0.25 (Figure 5).

It can be seen that there is a single equilibrium point B_e in every case, except in situations where there are two equilibrium points B_e in a small neighborhood of f = 0 for a = 0.5 and b = 1 (Figure 4).

The stability of the steady-state motion is determined by the eigenvalues of the Jacobian matrix obtained from Equations (66) and (67) (the Routh–Hurwitz criteria):

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

(71)

where

$$a_{11}={\left.{\displaystyle \frac{\partial A\prime }{\partial A}}\right|}_{(A_e,B_e)}\begin{array}{cc};& a_{12}=\end{array}{\left.{\displaystyle \frac{\partial A\prime }{\partial B}}\right|}_{(A_e,B_e)}\begin{array}{cc};& a_{21}={\left.{\displaystyle \frac{\partial B\prime }{\partial A}}\right|}_{(A_e,B_e)}\end{array}\begin{array}{cc};& a_{22}={\left.{\displaystyle \frac{\partial B\prime }{\partial A}}\right|}_{(A_e,B_e)}\end{array}\begin{array}{cc};& A\prime ={\displaystyle \frac{\partial A}{\partial \mathsf{\eta }}}\end{array},B\prime ={\displaystyle \frac{\partial B}{\partial \mathsf{\eta }}}$$

(72)

After some manipulations, from Equations (66), (67) and (69), the coefficients a_ij have the form

$$a_{11}=a_{22}=0\begin{array}{cc},& a_{12}={\displaystyle \frac{(18a+25b{B}_e^2)B_e^2}{16\overline{\mathsf{\omega }}}}\end{array}\begin{array}{cc},& a_{21}=-{\displaystyle \frac{(6a+5b{B}_e^2)B_e^2}{16\overline{\mathsf{\omega }}}}\end{array}$$

(73)

The eigenvalues of the Jacobian matrix are obtained from the characteristic equation

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

(74)

where I is the unity matrix of the second order, and λ is the eigenvalue of the Jacobian matrix. The characteristic equation can be rewritten as

$${\lambda }^2+{\displaystyle \frac{(18a+25b{B}_e^2)(6a+5b{B}_e^2)B_e^4}{256{\overline{\mathsf{\omega }}}^2}}=0$$

(75)

We have the following possible cases:

Case 6.1 If ab $\ge$ 0 and a² + b² $\ne$ 0, then the eigenvalues are purely imaginary with opposite signs of their imaginary parts. This is a center, such that there is no net motion towards or away from the equilibrium point.

Case 6.2 If ab < 0 and $B_e^2=-{\displaystyle \frac{18a}{25b}}$ or $B_e^2=-{\displaystyle \frac{6a}{5b}}$, then λ ₁= λ ₂ = 0. The situation is similar to the above.

Case 6.3 If ab < 0 and $-{\displaystyle \frac{18a}{25b}}<B_e^2<-{\displaystyle \frac{6a}{5b}}$, then λ₁ and λ₂ are real, and λ₁ = −λ₂. This case corresponds to a saddle point.

Case 6.4 If ab < 0 and $B_e^2<-{\displaystyle \frac{18a}{25b}}$ or $B_e^2>-{\displaystyle \frac{6a}{5b}}$, then the eigenvalues are purely imaginary, as in the situation of case 6.1.

7. Results and Discussion

In Figure 6, Figure 7, Figure 8 and Figure 9, the effects of different parameters on the vertical displacements of the Euler–Bernoulli beam are shown.

From Figure 6, we deduce that, if the curvature parameter β increases, then the abscissa of the intersection point of the vertical displacement and the horizontal axis decreases. Different values of β are taken as 0.01, 0.02, and 0.03. The effect of the nonlinear Winkler–Pasternak parameter K₃ on the vertical displacement is shown in Figure 7. This figure shows that, if the parameter increases from K₃ = 0.4 to K₃ = 0.6, then the abscissa of the intersection point of the vertical displacement and the horizontal axis increases. The effect of the gap parameter α is illustrated in Figure 8; the abscissa of the intersection point of the vertical displacement and the horizontal axis decreases with the growth of the gap parameter. From Figure 9, it follows that, if the electromagnetic parameter γ₁ increases, then the above-mentioned abscissa increases even more substantially for the time t > 14.

It can be observed that the longitudinal inertia term u(x,t) can be neglected only if E/ρ satisfies the condition (83).

8. Conclusions

According to the present results, the oscillatory behavior of a simply supported flexible uniform microbeam is studied. The Euler–Bernoulli microbeam is subjected to a mechanical impact using the Dirac delta function, an electromagnetic actuation, a nonlinear Winkler–Pasternak foundation, and a longitudinal magnetic field. In the governing differential equations, the curvature of the beam, the second-order approximation of the deflected beam, and the boundary conditions are taken into consideration. The OHAM is applied to solve the complex nonlinear differential equation, introducing so-called auxiliary functions and some convergence-control parameters without supplementary hypotheses. These parameters are optimally determined by rigorous mathematical procedures. The proper procedure leads to a very accurate solution using only the first iteration. It should be emphasized that the nonlinear dynamical system is reduced to only two linear differential equations.

The main aspects of novelty presented in this work, which are not present in previous related research, are summarized as follows:

(a)

The combination of the effects of the curvatures, mechanical impact, electromagnetic actuation, nonlinear elastic Winkler–Pasternak foundation, and longitudinal magnetic field. This combination is ensured by the new dynamical model which is considered in the study.

(b)

The Optimal Homotopy Asymptotic Method is used to solve a very complex differential equation by means of some auxiliary functions, providing novel solutions.

(c)

The values of the convergence-control parameters are optimally determined, ensuring a high degree of accuracy.

(d)

It is proved that the longitudinal inertia term can be neglected, but only if E/ρ satisfies the condition (83).

(e)

Local stability for the primary resonance is studied by means of the Homotopy Perturbation Method.

The effect of some physical parameters is highlighted. The local stability depends on the nature of the solutions to the characteristic equation, and different cases are discussed. The local stability near the primary resonance is achieved by using the Homotopy Perturbation Method, the variable expansion method, equilibrium points, and the Jacobian matrix. The effects of different parameters on the vertical displacements of the beam are considered.