Nonlinear Vibration Control of a Piezoelectric–Elastic–Piezoelectric Sandwich Beam

1. Introduction

Structural vibrations are highly undesirable, as they can lead to issues such as structural fatigue, transmission of vibrations to other systems, and external or internal noise due to acoustic radiation, among others [1,2,3]. In many industrial and defense applications, noise and vibrations represent a major challenge. Conventional mitigation methods, which rely on passive damping techniques, often prove ineffective at low frequencies. In this context, active control methods appear to be more suitable [4,5,6,7,8,9]. The principle of these so-called active techniques is to generate a field that interferes with the disturbance field. The superimposed field must therefore match the disturbance in amplitude but be opposite in phase for each relevant frequency. While the principle is straightforward, its implementation is much more complex, as the disturbance is often unpredictable and composed of multiple frequencies. Moreover, disturbance minimization is often required over a wide spatial domain, further complicating the problem [10].

Although active control was conceived in the 1930s, it only truly advanced with the emergence of digital signal processors in the 1980s. While some applications of this technology have already been developed, many are still under research, particularly in the aerospace, avionics, and automotive sectors. Focusing specifically on active vibration control, advancements remain relatively recent. In fact, the additional size and mass introduced by the sensors and actuators required for active vibration control have long hindered the development of many applications. Only in the last few decades has the use of piezoelectric material-based transducers enabled significant progress. Due to their compactness, low weight, and electro-mechanical conversion capabilities, piezoelectric materials exhibit all the necessary qualities for use in active vibration control systems. Moreover, they can serve both as electromechanical actuators and vibration sensors in the system [11,12].

To reduce stresses in materials, extend their service life, enhance structural safety (e.g., in transportation), and improve user comfort, the control and damping of mechanical vibrations has been the subject of extensive scientific research over many decades. Furthermore, the recent proliferation of so-called “smart materials” which couple multiple physical fields such as mechanics and electricity has led to the development of reliable, robust, and efficient vibration control techniques that are also highly integrable. These techniques are therefore well-suited for embedded systems or structures with strict space constraints [11]. In this regard, Rechdaoui et al. [12,13,14] developed an active control method for nonlinear vibrations of a piezoelectric–elastic sandwich beam based on the method of multiple scales. Similarly, Belouettar et al. [15] proposed an active control approach for the same structure based on the harmonic balance method. Despite such contributions, vibration-related damage remains a persistent problem in our societies, and many avenues still remain to be explored hence the ongoing need for research efforts.

This work aims to contribute to the active control of nonlinear vibrations of a sandwich beam using the method of averaging [16,17,18,19,20].

1.1. Mathematical Modeling

1.1.1. Theoretical Formulation

The beam under investigation consists of an elastic core sandwiched between two piezoelectric face sheets polarized through their thickness, as shown in Figure 1. Euler–Bernoulli beam theory is applied to the face sheets, which are assumed to resist membrane and bending stresses. Timoshenko beam theory is adopted for the core, which is assumed to also resist transverse shear stress. The piezoelectric layers are fully covered on their top and bottom surfaces with electrodes. The elastic and piezoelectric materials are considered orthotropic, with their orthotropy axes aligned with those of the sandwich beam. All layers are assumed to be perfectly bonded. The transverse normal stress is considered negligible compared to the other stress components [21,22,23,24]. The length, width, and thickness of the beam are denoted by L, H, and h, respectively. The subscripts a, s, and c refer to the quantities associated with the bottom and top face sheets, and the core, respectively.

1.1.2. Kinematic Description of the Beam

According to the classical laminate theory based on the aforementioned assumptions, the displacement fields are described by [1,11,13,24].

$$\left\{\begin{array}{l}u(x,z,t)=u(x,t)-z{w}_{,x}(x,t)\\ v(x,z,t)=0\\ w(x,z,t)=w(x,t)\end{array}\right.$$

(1)

In the theory of beams undergoing large deformations, the Green-Lagrange strain tensor is considered without linearization and simplified according to the Von Kármán assumptions.

$u$ : longitudinal displacement;

$w$ : transverse displacement;

$z$ : coordinate along the beam thickness (thickness direction).

Thus, the strain in the x-direction is given by:

$$\left\{\begin{array}{l}\epsilon ={\epsilon }_0-z{w}_{,xx}\\ {\epsilon }_0=u_{,x}+{\displaystyle \frac{1}{2}}{w}_{,x}^2\end{array}\right.$$

(2)

1.1.3. Electromechanical Coupling

It is well known that in piezoelectric materials, the electric field and strain mutually influence each other. This property enables the use of piezoelectric materials as sensors and actuators for vibration control. More specifically, this relationship can be described by constitutive equations that characterize the coupling effects between mechanical and electrical properties [1,24].

$$\left\{\begin{array}{l}\sigma =c\epsilon -e^{t}E\\ D=e\epsilon +\in E\end{array}\right.$$

(3)

$\sigma$ : Stress tensor;

$\epsilon$ : Strain tensor;

D : Electric displacement vector;

E : Electric field vector;

$c$ : Elasticity matrix;

$e$ : Piezoelectric matrix;

ϵ : Dielectric permittivity.

The constitutive equations of piezoelectric materials can be written in the following expanded form:

$$\left\{\begin{array}{l}\left\{\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right\}\\ \left\{\begin{array}{l}{D}_1\\ D_2\\ D_3\end{array}\right\}\end{array}\right\}=\left[\begin{array}{cc}\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& 0& 0& 0\\ c_{21}& c_{22}& c_{23}& 0& 0& 0\\ c_{31}& c_{32}& c_{33}& 0& 0& 0\\ 0& 0& 0& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{55}& 0\\ 0& 0& 0& 0& 0& c_{66}\end{array}\right]& \left[\begin{array}{ccc}0& 0& -e_{13}\\ 0& 0& -e_{13}\\ 0& 0& -e_{13}\\ 0& -e_{13}& 0\\ -e_{13}& 0& 0\\ 0& 0& 0\end{array}\right]\\ \left[\begin{array}{cccccc}0& 0& 0& 0& e_{15}& 0\\ 0& 0& 0& e_{24}& 0& 0\\ e_{31}& e_{32}& e_{33}& 0& 0& 0\end{array}\right]& \left[\begin{array}{ccc}{\in }_{11}& 0& 0\\ 0& {\in }_{22}& 0\\ 0& 0& {\in }_{33}\end{array}\right]\end{array}\right]\left\{\begin{array}{l}\left\{\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\\ {\epsilon }_4\\ {\epsilon }_5\\ {\epsilon }_6\end{array}\right\}\\ \left\{\begin{array}{l}{E}_1\\ E_2\\ E_3\end{array}\right\}\end{array}\right\}$$

(4)

Displacements are considered independent of y and zero along the y-direction; the stress tensor is uniaxial, and the directions of the vectors D and E are parallel to the z-axis. Thus, the reduced constitutive relations are given by:

$$\left\{\begin{array}{l}{\epsilon }_2=0\\ {\sigma }_3=0\end{array}\right.\hspace{0.33em}\hspace{1em}\Rightarrow \hspace{0.33em}\hspace{0.33em}{\epsilon }_3={\displaystyle \frac{1}{c_{33}}}(e_{33}{E}_3-c_{13}\epsilon )$$

(4)

It follows that:

$$\left(\begin{array}{l}{\sigma }_1\\ D_3\end{array}\right)=\left[\begin{array}{cc}{c}_{11}^{*}& -e_{31}^{*}\\ e_{31}^{*}& {\in }_{33}^{*}\end{array}\right]\left(\begin{array}{l}\epsilon \\ E_3\end{array}\right)$$

(5)

With: ${\in }_{33}^{*}={\in }_{33}+{\displaystyle \frac{e_{33}^2}{c_{33}}}\hspace{0.17em}; e_{31}^{*}=e_{31}-{\displaystyle \frac{c_{13}}{c_{33}}}{e}_{33}\hspace{0.17em}; c_{11}^{*}=c_{11}-{\displaystyle \frac{c_{13}^2}{c_{33}}}.$

1.1.4. Feedback Control Law

Let us now consider an arbitrary piezoelectric layer, actuator or sensor, placed between $\hspace{0.17em}{z}_{-}, z_{+}\hspace{0.17em}(z_{-}<z_{+})$, with center$\hspace{0.17em}{z}_c=(z_{-}+z_{+})/2$ and thickness h. The electrostatic equilibrium equation, assuming no volume charge density, is given by:

$${\displaystyle \frac{\partial D_3}{\partial z}}=0$$

(6)

Using the boundary conditions,$\hspace{0.17em}D(z_{-})=0\hspace{0.17em}\hspace{0.17em}ou\hspace{0.17em}\hspace{0.17em}D(z_{+})=0$, we have:

$$D_3(z)=0$$

(7)

Thus, the electric field in the sensor, as a function of displacement, is given by:

$$E_3(z)=-{\displaystyle \frac{e_{31}^{*}}{{\in }_{33}^{*}}}\epsilon =-{\displaystyle \frac{e_{31}^{*}}{{\in }_{33}^{*}}}(u_{,x}+{\displaystyle \frac{1}{2}}{w}_{,x}^2-z{w}_{,xx})$$

(8)

Since the electric field derives from a potential, we have:

$$E_3=-{\displaystyle \frac{\partial \varphi }{\partial z}}$$

(9)

Consequently, the potential difference is given by:

$$\Delta \varphi =\varphi (z_{+})-\varphi (z_{-})=-{\displaystyle {\int }_{z_{-}}^{z_{+}}{E}_3(z)dz=\frac{e_{31}^{*}}{{\in }_{33}^{*}}{h}_i(u_{,x}+\frac{1}{2}{w}_{,x}^2-z_i{w}_{,xx})}$$

(10)

With: $\left\{\begin{array}{l}{z}_{-}=z_i-{\displaystyle \frac{h_i}{2}}\\ z_{+}=z_i+{\displaystyle \frac{h_i}{2}}\end{array}\right.\hspace{0.33em}\hspace{1em}i=a,\hspace{0.17em}s.$

From Equations (8) and (10), we have:

$$E_3(z)=-{\displaystyle \frac{\Delta \varphi }{h_i}}+{\displaystyle \frac{e_{31}^{*}}{{\in }_{33}^{*}}}(z-z_i)w_{,xx}$$

(11)

The core of the beam is assumed to be conductive with a uniform potential set to zero. The sensor potential, denoted$\hspace{0.17em}{\varphi }_s(x)$, is then given by:

$${\varphi }_s=\Delta {\varphi }_s={\displaystyle \frac{e_{31}^{*}}{{\in }_{33}^{*}}}{h}_s(u_{,x}+{\displaystyle \frac{1}{2}}{w}_{,x}^2-z_s{w}_{,xx})$$

(12)

The actuator potential$\hspace{0.17em}{\varphi }_a(x)$ depends on the sensor output potential$\hspace{0.17em}{\varphi }_s(x)$ through a proportional-derivative control law described by:

$${\varphi }_a=G_p{\varphi }_s+G_d{\dot{\varphi }}_s$$

(13)

Using Equations (11) and (12), the electric fields in the sensor and actuator are respectively given by:

$$E_3^s(z)=-{\displaystyle \frac{{\varphi }_s}{h_s}}+{\displaystyle \frac{e_{31}^{*}}{{\in }_{33}^{*}}}(z-z_s)w_{,xx}$$

(14)

$$E_3^a(z)={\displaystyle \frac{{\varphi }_a}{h_a}}+{\displaystyle \frac{e_{31}^{*}}{{\in }_{33}^{*}}}(z-z_a)w_{,xx}$$

(15)

The potentials $\hspace{0.17em}{\varphi }_a$ and$\hspace{0.17em}{\varphi }_s$ are independent of z and$\hspace{0.17em}{z}_s=(h_c+h_s)/2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}et\hspace{0.17em}\hspace{0.17em}\hspace{0.33em}{z}_a=-(h_c+h_s)/2$.

The direct and inverse piezoelectric coefficients were taken into account in these formulations and will be involved in the dynamic behavior of the beam [1,12,25,26].

1.1.5. Dynamic Equation

To determine the dynamic equation of the beam, we use the variational formulation, Hamilton’s principle, and Equations (5), (14), and (15). The beam is subjected to axial and transverse excitations F_x and F_z [1,24,27,28].

For the variational principle, we have:

$$\begin{array}{l}{\displaystyle \underset{V}{\int }{\sigma }_1\delta \epsilon }dv={\displaystyle \underset{V_s}{\int }{\sigma }_1\delta \epsilon d{v}_s}+{\displaystyle \underset{V_c}{\int }{\sigma }_1\delta \epsilon d{v}_c}+{\displaystyle \underset{V_a}{\int }{\sigma }_1\delta \epsilon d{v}_a}={\displaystyle {\int }_0^L(N\delta {\epsilon }_0+M\delta w_{,xx})dx}\\ ={\displaystyle {\int }_0^L(F_X\delta u+F_Z\delta w)dx}-{(\rho S)}_{*}{\displaystyle {\int }_0^L(\ddot{u}\delta u+\ddot{w}\delta w)dx}\end{array}$$

(16)

With: $\left\{\begin{array}{l}\hspace{0.33em}{\left(\rho \mathrm{S}\right)}_{*}={\rho }_{\mathrm{s}}{\mathrm{S}}_{\mathrm{s}}+{\rho }_c{\mathrm{S}}_{\mathrm{c}}+{\rho }_a{\mathrm{S}}_{\mathrm{a}}\\ \delta {\epsilon }_0=\delta u_{,x}+w_{,x}\delta w_{,x}\end{array}\right.$

According to Hamilton’s principle, we have:

$$\left\{\begin{array}{l}{N}_{xx}={\displaystyle \underset{s}{\int }{\sigma }_{xx}ds=ES\left[\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^2\right]}\\ M_{xx}={\displaystyle \underset{s}{\int }z{\sigma }_{xx}ds}=-EI{\displaystyle \frac{{\partial }^{2}w}{\partial x^2}}\end{array}\right.$$

(17)

If we integrate over the entire thickness and width, and assume that the piezoelectric layers are symmetrical$\hspace{0.17em}(h_a=h_s)$, the axial force N and the bending moment M are determined from the previous equation as follows:

$$N={\sigma }_1^a{S}_a+{\sigma }_1^c{S}_c+{\sigma }_1^s{S}_s=\left({\sigma }_1^a+{\sigma }_1^s\right)S_s+E_c{S}_c\epsilon$$

(18)

Using Equation (5), we obtain:

$${\sigma }_1^a+{\sigma }_1^s=2c_{11}^{*}\epsilon -e_{31}^{*}\left(E_3^a+E_3^s\right)$$

(19)

$$E_3^a+E_3^s={\displaystyle \frac{{\varphi }_a-{\varphi }_s}{h_s}}$$

(20)

Using Equations (2), (12)–(15), we obtain:

$$\begin{array}{l}N=\left(2c_{11}^{*}{S}_s+E_c{S}_c+\left(1-G_p\right)S_s{\displaystyle \frac{{(e_{31}^{*})}^2}{{\in }_{33}^{*}}}\right){\epsilon }_0-\left(1-G_p\right)S_s{\displaystyle \frac{{(e_{31}^{*})}^2}{{\in }_{33}^{*}}}{z}_s{w}_{,xx}\\ \hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}-G_d{S}_s{\displaystyle \frac{{(e_{31}^{*})}^2}{{\in }_{33}^{*}}}\left({\dot{u}}_{,x}+w_{,x}{\dot{w}}_x-z_s{w}_{,xx}\right)\end{array}$$

(21)

We also obtain the moment in the same manner; thus, we have:

$$\left\{\begin{array}{l}N={(ES)}_{*}{\epsilon }_0-B_N{w}_{,xx}-{(ES)}_{pe}{G}_d({\dot{u}}_{,x}+w_{,x}{\dot{w}}_{,x}-{\dot{w}}_{,xx}{z}_s)\\ M=-B_M{\epsilon }_0+{(EI)}_{*}{w}_{,xx}-{(ES)}_{pe}{z}_s{G}_d({\dot{u}}_{,x}+w_{,x}{\dot{w}}_{,x}-{\dot{w}}_{,xx}{z}_s)\end{array}\right.$$

(22)

With: $\left\{\begin{array}{l}{(ES)}_{*}=E_c{S}_c+2c_{11}^{*}{S}_s+{(ES)}_{pe}(1-G_p);\\ {(ES)}_{pe}=S_s{\displaystyle \frac{{(e_{31}^{*})}^2}{{\in }_{33}^{*}}};\\ B_N={(ES)}_{pe}(1-G_p)z_s;\\ B_M={(ES)}_{pe}(1+G_p)z_s;\\ {(EI)}_{*}=E_c{I}_c+2c_{11}^{*}(I_s+S_s{z}_s^2)+{\displaystyle \frac{{(ES)}_{pe}}{S_s}}(2I_s+(1+G_p)z_s^2{S}_s).\end{array}\right.$

By applying the variational principle to the displacements u and w, and integrating by parts once for the terms in$\hspace{0.17em}\delta u_{,x}$ and$\hspace{0.17em}\delta w_{,x}$, and twice for the terms in $\hspace{0.17em}\delta w_{,xx}$, we obtain the following partial differential equations [1,12,13,14,15].

$$\left\{\begin{array}{l}-N_{,x}+{(\rho S)}_{*}\ddot{u}=F_X\\ M_{,xx}-{(N{w}_{,x})}_{,x}+{(\rho S)}_{*}\ddot{w}=F_Z\end{array}\right.$$

(23)

Assuming that the axial force and the axial displacement inertia are negligible, system (23) becomes:

$$\left\{\begin{array}{l}-N_{,x}=0\\ M_{,xx}-{(N{w}_{,x})}_{,x}+{(\rho S)}_{*}\ddot{w}=F_Z\end{array}\right.$$

(24)

The axial force depending only on time$\left(N(x,t)=N(t)\right)$, system (24) becomes:

$$\left\{\begin{array}{l}N={\displaystyle \frac{1}{2}}{(ES)}_{*}{w}_{,x}^2-B_N{w}_{,xx}-{(ES)}_{pe}{G}_d(w_{,x}{\dot{w}}_{,x}-{\dot{w}}_{,xx}{z}_s)\\ M={\displaystyle \frac{1}{2}}{B}_M{w}_{,x}^2+{(EI)}_{*}{w}_{,xx}-{(ES)}_{pe}{G}_d{z}_s(w_{,x}{\dot{w}}_{,x}-{\dot{w}}_{,xx}{z}_s)\end{array}\right.$$

(25)

En intégrant la première équation du système (25) entre (O et L), on obtient :

$$N(t)={\displaystyle \frac{1}{2L}}{(ES)}_{*}{\displaystyle {\int }_0^L{w}_{,x}^{2}dx}-{\displaystyle \frac{B_N}{L}}{\displaystyle {\int }_0^L{w}_{,xx}dx}-{\displaystyle \frac{{(ES)}_{pe}}{L}}{G}_d{\displaystyle {\int }_0^L(w_{,x}{\dot{w}}_{,x}-{\dot{w}}_{,xx}{z}_s)dx}$$

(26)

By substituting (26) into the second equation of system (25), the following dynamic equation is obtained:

$$\begin{array}{l}{(\rho S)}_{*}\ddot{w}+{(EI)}_{*}{w}_{,xxxx}-N(t)w_{,xx}-B_M(w_{,xx}^2+w_{,x}{w}_{,xxx})\\ \hspace{1em}-{(\rho S)}_{pe}{G}_d{z}_s({\dot{w}}_{,x}{w}_{,xxx}+2{\dot{w}}_{,xx}{w}_{,xx}+w_{,x}{\dot{w}}_{,xxx}-z_s{\dot{w}}_{,xxxx})=F_Z\end{array}$$

(27)

This differential equation describes the transverse dynamic behavior of the piezoelectric-elastic-piezoelectric beam subjected to a transverse excitation and active control based on the feedback control law, when the axial force and axial displacement effects are neglected. The free and forced nonlinear vibrations, as well as the active control of the beam, can be analyzed by solving Equation (24) or (27) [12].

In this work, the axial effects are neglected; therefore, only Equation (27), which describes the dynamic behavior of the beam without axial effects, will be solved.

2. Solution Methodology

To solve Equation (27), and in order to simplify the calculations, the Galerkin approximation given by Equation (29) below is applied. The beam is transversely excited by an external uniformly distributed harmonic force of the form:

$$F_z(x,t)=f(x)\cos (\omega t)$$

(28)

The Galerkin approximation [12,28] is given by:

$$w(x,t)={\displaystyle \sum _{k=1}^n{q}_k(t){\phi }_k(x)}$$

(29)

$\hspace{0.33em}{\phi }_k(x)\hspace{0.17em}$ : are the vibration modes of the beam;

$q_k(t)$ : are the associated time-dependent amplitudes.

This mode superposition leads to a reduced-order approximate dynamic system model. To perform control in a simple manner, we consider a single mode. By substituting Equation (29) into Equation (27), integrating over the entire length, and omitting the indices since only one mode is considered, we obtain:

$$\ddot{q}(t)+2\mu \dot{q}(t)+{\omega }_L^{2}q(t)+{\alpha }_2{q}^2(t)+{\alpha }_3{q}^3(t)+{\alpha }_{4}q(t)\dot{q}(t)+{\alpha }_5{q}^2(t)\dot{q}(t)=F_1\cos (\omega t)$$

(30)

$$\begin{array}{l}2\mu ={\displaystyle \frac{{(ES)}_{pe}{G}_d{z}_s^2}{M}}{\displaystyle {\int }_0^L{\phi }_{xxxx}(x)}\phi (x)dx\hspace{0.33em};\\ {\omega }_L^2={\displaystyle \frac{{(EI)}_{*}}{M}}{\displaystyle {\int }_0^L{\phi }_{xxxx}(x)}\phi (x)dx;\\ {\alpha }_2={\displaystyle \frac{B_N}{ML}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)}\phi (x)dx{\displaystyle {\int }_0^L{\phi }_{xx}(x)}dx-{\displaystyle \frac{B_M}{M}}{\displaystyle {\int }_0^L\left\{{({\phi }_{xx}(x))}^2+{\phi }_x(x){\phi }_{xxx}(x)\right\}\phi (x)dx}\hspace{0.33em};\\ {\alpha }_3=-{\displaystyle \frac{{(ES)}_{*}}{2ML}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)}\phi (x)dx{\displaystyle {\int }_0^L({\phi }_x(x)}{)}^{2}dx\hspace{0.33em};\\ {\alpha }_4=-{\displaystyle \frac{{(ES)}_{pe}{G}_d{z}_s}{M}}\left({\displaystyle \frac{1}{L}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)}\phi (x)dx{\displaystyle {\int }_0^L{\phi }_{xx}(x)}dx+2\left({\displaystyle {\int }_0^L{({\phi }_{xx}(x))}^2+{\phi }_x(x){\phi }_{xxx}(x)\phi (x)dx}\right)\right)\hspace{0.33em};\\ {\alpha }_5={\displaystyle \frac{{(ES)}_{pe}{G}_d}{ML}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)}\phi (x)dx{\displaystyle {\int }_0^L({\phi }_x(x)}{)}^{2}dx\hspace{0.33em};\\ F_1={\displaystyle \frac{1}{M}}{\displaystyle {\int }_0^{L}f(x)}\phi (x)dx\hspace{0.33em};\\ M={(\rho S)}_{*}{\displaystyle {\int }_0^L({\phi }_x(x)}{)}^{2}dx\hspace{0.33em}.\end{array}$$

These coefficients depend on the control parameters G_p and G_d, and consequently, they can be significantly influenced by the control law considered. The resolution of Equation (30) will be carried out in the following using the method of averaging.

2.1. Primary Resonance

According to the principle of the method of averaging, Equation (30) can be rewritten by introducing the perturbation parameter, thus we have:

$$\ddot{q}+{\omega }_L^{2}q=-\epsilon \left[2\mu \dot{q}+{\alpha }_2{q}^2+{\alpha }_3{q}^3+{\alpha }_{4}q\dot{q}+{\alpha }_5{q}^2\dot{q}-F_1\cos (\omega t)\right]$$

(31)

If$\hspace{0.17em}\epsilon =0$, the general solution of Equation (31) is given by:

$$q=a\cos ({\omega }_{L}t+\beta )$$

(32)

Since$a$and$\beta$ are constants, the derivative of Equation (32) is:

$$\dot{q}=-a{\omega }_L\sin ({\omega }_{L}t+\beta )$$

(33)

If$\hspace{0.17em}\epsilon \ne 0$, the solution of Equation (31) takes the form of Equation (33), but with$\hspace{0.17em}a$and$\hspace{0.17em}\beta$now varying with time. Differentiating Equation (32) then yields:

$$\dot{q}=-a{\omega }_L\sin ({\omega }_{L}t+\beta )+\dot{a}\cos ({\omega }_{L}t+\beta )-a\dot{\beta }\sin ({\omega }_{L}t+\beta )$$

(34)

By comparing Equations (33) and (34), we obtain:

$$\dot{a}\cos ({\omega }_{L}t+\beta )-a\dot{\beta }\sin ({\omega }_{L}t+\beta )=0$$

(35)

Let us differentiate Equation (33) with respect to time:

$$\ddot{q}=-\dot{a}{\omega }_L\sin ({\omega }_{L}t+\beta )-a{\omega }_L^2\cos ({\omega }_{L}t+\beta )-a\dot{\beta }{\omega }_L\cos ({\omega }_{L}t+\beta )$$

(36)

By substituting (33), (34), and (36) into (31), we obtain:

$$\begin{array}{l}\dot{a}{\omega }_L\sin ({\omega }_{L}t+\beta )+a\dot{\beta }{\omega }_L\cos ({\omega }_{L}t+\beta )=-2\mu \epsilon \left(a{\omega }_L\sin ({\omega }_{L}t+\beta )\right)\\ +{\alpha }_2\epsilon \left(a^2{\cos }^2({\omega }_{L}t+\beta )\right)+{\alpha }_3\epsilon \left(a^3{\cos }^3({\omega }_{L}t+\beta )\right)\\ -{\alpha }_4\epsilon \left(a^2{\omega }_L\cos ({\omega }_{L}t+\beta )\sin ({\omega }_{L}t+\beta )\right)\\ -{\alpha }_5\epsilon \left(a^3{\omega }_L{\cos }^2({\omega }_{L}t+\beta )\sin ({\omega }_{L}t+\beta )\right)-\epsilon F_1\cos (\omega t)\end{array}$$

(37)

By using (35) and (37), we have:

$$\dot{a}={\displaystyle \frac{\epsilon }{{\omega }_L}}\left\{\begin{array}{l}-\mu a{\omega }_L+\mu a{\omega }_L\cos (2{\omega }_{L}t+2\beta )+{\alpha }_2\left({\displaystyle \frac{1}{4}}{a}^2\sin ({\omega }_{L}t+\beta )+{\displaystyle \frac{1}{4}}{a}^2\sin (3{\omega }_{L}t+3\beta )\right)\\ +{\alpha }_3\left({\displaystyle \frac{1}{4}}{a}^3\sin (2{\omega }_{L}t+2\beta )+{\displaystyle \frac{1}{8}}{a}^3\sin (4{\omega }_{L}t+4\beta )\right)-{\displaystyle \frac{1}{2}}{F}_1\sin \left(({\omega }_L-\omega )t+\beta \right)\\ -{\alpha }_4\left({\displaystyle \frac{1}{4}}{a}^2{\omega }_L\cos ({\omega }_{L}t+\beta )-{\displaystyle \frac{1}{4}}{a}^2{\omega }_L\cos (3{\omega }_{L}t+3\beta )\right)-{\displaystyle \frac{1}{2}}{F}_1\sin \left(({\omega }_L+\omega )t+\beta \right)\\ -{\alpha }_5\left({\displaystyle \frac{1}{8}}{a}^3{\omega }_L-{\displaystyle \frac{1}{8}}{a}^3{\omega }_L\cos (4{\omega }_{L}t+4\beta )\right)\end{array}\right\}$$

(38)

By substituting (38) into (35), we obtain:

$$a\dot{\beta }={\displaystyle \frac{\epsilon }{{\omega }_L}}\left\{\begin{array}{l}-\mu a{\omega }_L\sin (2{\omega }_{L}t+2\beta )+{\alpha }_2\left({\displaystyle \frac{3a^2}{4}}\cos ({\omega }_{L}t+\beta )+{\displaystyle \frac{a^2}{4}}\cos (3{\omega }_{L}t+3\beta )\right)\\ +{\alpha }_3\left({\displaystyle \frac{3a^3}{8}}+{\displaystyle \frac{a^3}{2}}\cos (2{\omega }_{L}t+2\beta )+{\displaystyle \frac{a^3}{8}}\cos (4{\omega }_{L}t+4\beta )\right)\\ -{\alpha }_4\left({\displaystyle \frac{a^2}{4}}{\omega }_L\sin ({\omega }_{L}t+\beta )+{\displaystyle \frac{a^2}{4}}{\omega }_L\sin (3{\omega }_{L}t+3\beta )\right)-{{\displaystyle \frac{F_1}{2}}}_1\cos \left(({\omega }_L-\omega )t+\beta \right)\\ -{\alpha }_5\left({\displaystyle \frac{a^3}{4}}{\omega }_L\sin (2{\omega }_{L}t+2\beta )+{\displaystyle \frac{a^3}{8}}{\omega }_L\sin (4{\omega }_{L}t+4\beta )\right)-{\displaystyle \frac{F_1}{2}}\cos \left(({\omega }_L+\omega )t+\beta \right)\end{array}\right\}$$

(39)

In primary resonance,$\hspace{0.17em}\omega \approx {\omega }_L$ and the expressions in$\hspace{0.33em}\sin \left(({\omega }_L-\omega )t+\beta \right)$ and $\cos \left(({\omega }_L-\omega )t+\beta \right)$ vary slowly with respect to time in Equations (38) and (39), respectively. We then have:

$$\left\{\begin{array}{l}\dot{a}=-\mu a\epsilon -{\alpha }_5\epsilon {\displaystyle \frac{a^3}{8}}-{\displaystyle \frac{\epsilon F_1}{2{\omega }_L}}\sin \left(({\omega }_L-\omega )t+\beta \right)\\ a\dot{\beta }={\alpha }_3\epsilon {\displaystyle \frac{3a^3}{8{\omega }_L}}-{\displaystyle \frac{\epsilon F_1}{2{\omega }_L}}\cos \left(({\omega }_L-\omega )t+\beta \right)\end{array}\right.$$

(40)

By setting $\gamma =({\omega }_L-\omega )t+\beta$ and$\hspace{0.17em}\omega ={\omega }_L+\epsilon \sigma$, the system (40) becomes:

$$\left\{\begin{array}{l}\dot{a}=-\mu a\epsilon -{\alpha }_5\epsilon {\displaystyle \frac{a^3}{8}}-{\displaystyle \frac{\epsilon F_1}{2{\omega }_L}}\sin (\gamma )\\ a\dot{\gamma }+\epsilon a\sigma ={\alpha }_3\epsilon {\displaystyle \frac{3a^3}{8{\omega }_L}}-{\displaystyle \frac{\epsilon F_1}{2{\omega }_L}}\cos (\gamma )\end{array}\right.$$

(41)

Initially, $a$ and $\gamma \hspace{0.17em}$oscillate, and as time increases, they become constants. Thus, for$\hspace{0.17em}\dot{a}=0$ and $\dot{\gamma }=0$, we have:

$$\left\{\begin{array}{l}{\displaystyle \frac{\mu }{{\omega }_L}}+{\displaystyle \frac{{\alpha }_5{a}^2}{8{\omega }_L}}={\displaystyle \frac{-F_1}{2a{\omega }_L^2}}\sin (\gamma )\\ \sigma -3{\displaystyle \frac{{\alpha }_3{a}^2}{8{\omega }_L}}={\displaystyle \frac{-F_1}{2a{\omega }_L^2}}\cos (\gamma )\end{array}\right.$$

(42)

From system (42), the following equation is obtained:

$${\left[{\displaystyle \frac{\omega }{{\omega }_L}}-\left(1+3{\displaystyle \frac{{\alpha }_3{a}^2}{8{\omega }_L^2}}\right)\right]}^2+{\left[{\displaystyle \frac{\mu }{{\omega }_L}}+{\displaystyle \frac{{\alpha }_5{a}^2}{8{\omega }_L}}\right]}^2={\left({\displaystyle \frac{F_1}{2a{\omega }_L^2}}\right)}^2$$

(43)

2.2. Secondary Resonance

In the case of secondary resonance, Equation (30) can be rewritten by introducing the perturbation parameter in the following form:

$$\ddot{q}+{\omega }_L^{2}q=-\epsilon \left[2\mu \dot{q}+{\alpha }_2{q}^2+{\alpha }_3{q}^3+{\alpha }_{4}q\dot{q}+{\alpha }_5{q}^2\dot{q}\right]+F_1\cos (\omega t)$$

(44)

If$\hspace{0.17em}\epsilon =0$ and using the principle of superposition, the general solution of Equation (44) is given by:

$$q=a\cos ({\omega }_{L}t+\beta )+2\Lambda \cos (\omega t)$$

(45)

$$2\Lambda ={\displaystyle \frac{F_1}{{\omega }_L^2-{\omega }^2}}$$

(46)

The derivative of (45) is given by:

$$\dot{q}=-a{\omega }_L\sin ({\omega }_{L}t+\beta )-2\Lambda \omega \sin (\omega t)$$

(47)

If$\hspace{0.17em}\epsilon \ne 0$, as in primary resonance, the derivative of (46), using the variation of constants, gives:

$$\dot{q}=-a{\omega }_L\sin ({\omega }_{L}t+\beta )+\dot{a}\cos ({\omega }_{L}t+\beta )-a\dot{\beta }\sin ({\omega }_{L}t+\beta )-2\Lambda \omega sin(\omega t)$$

(48)

By comparing (47) and (48), we obtain:

$$\dot{a}\cos ({\omega }_{L}t+\beta )-a\dot{\beta }sin({\omega }_{L}t+\beta )=0$$

(49)

The derivative of (47) gives:

$$\ddot{q}=-a{\omega }_L^2\cos ({\omega }_{L}t+\beta )-\dot{a}{\omega }_L\sin ({\omega }_{L}t+\beta )-a\dot{\beta }{\omega }_L\cos ({\omega }_{L}t+\beta )-2\Lambda {\omega }^2\cos (\omega t)$$

(50)

By substituting (45), (47), and (29) into (44), we obtain:

$$\begin{array}{l}\dot{a}{\omega }_L\sin ({\omega }_{L}t+\beta )+a\dot{\beta }\cos ({\omega }_{L}t+\beta )=-2\mu \epsilon \left[a{\omega }_L\sin ({\omega }_{L}t+\beta )+2\Lambda \omega sin(\omega t)\right]\\ +{\alpha }_2\epsilon {\left[a\cos ({\omega }_{L}t+\beta )+2\Lambda \cos (\omega t)\right]}^2+{\alpha }_3\epsilon {\left[a\cos ({\omega }_{L}t+\beta )+2\Lambda \cos (\omega t)\right]}^3\\ +{\alpha }_4\epsilon \left[\left(a\cos ({\omega }_{L}t+\beta )+2\Lambda \cos (\omega t)\right)\left(-a{\omega }_L\sin ({\omega }_{L}t+\beta )-2\Lambda \omega \sin (\omega t)\right)\right]\\ +{\alpha }_5\epsilon \left[{\left(a\cos ({\omega }_{L}t+\beta )+2\Lambda \cos (\omega t)\right)}^2\left(-a{\omega }_L\sin ({\omega }_{L}t+\beta )-2\Lambda \omega \sin (\omega t)\right)\right]\end{array}$$

(51)

Using (49) and (51), we obtain:

$$\begin{array}{l}\dot{a}={\displaystyle \frac{\epsilon }{{\omega }_L}}\left\{-\mu a{\omega }_L\left(1-\cos (2{\omega }_{L}t+2\beta )\right)+2\mu \Lambda \omega \cos \left(({\omega }_L+\omega )t+\beta \right)-2\mu \Lambda \omega \cos \left(({\omega }_L-\omega )t+\beta \right)\right\}\\ +{\displaystyle \frac{{\alpha }_2\epsilon }{{\omega }_L}}\left\{\begin{array}{l}\left(2{\Lambda }^2+{\displaystyle \frac{a^2}{4}}\right)\sin ({\omega }_{L}t+\beta )+{\displaystyle \frac{a^2}{4}}\sin (3{\omega }_{L}t+3\beta )+{\Lambda }^2\sin \left(({\omega }_L+2\omega )t+\beta \right)\\ +{\Lambda }^2\sin \left(({\omega }_L-2\omega )t+\beta \right)+a\Lambda \sin \left((2{\omega }_L+\omega )t+2\beta \right)+a\Lambda \sin \left((2{\omega }_L-\omega )t+2\beta \right)\end{array}\right\}\\ +{\displaystyle \frac{{\alpha }_3\epsilon }{{\omega }_L}}\left\{\begin{array}{l}\left(3a{\Lambda }^2+{\displaystyle \frac{a^3}{4}}\right)\sin (2{\omega }_{L}t+2\beta )+{\displaystyle \frac{a^3}{8}}\sin (4{\omega }_{L}t+4\beta )\\ \left(3{\Lambda }^3+{\displaystyle \frac{3a^2\Lambda }{4}}\right)\sin \left(({\omega }_L+\omega )t+\beta \right)+\left(3{\Lambda }^3+{\displaystyle \frac{3a^2\Lambda }{4}}\right)\sin \left(({\omega }_L-\omega )t+\beta \right)\\ +{\displaystyle \frac{3a^2\Lambda }{2}}\sin \left((2{\omega }_L+2\omega )t+2\beta \right)+{\displaystyle \frac{3a^2\Lambda }{2}}\sin \left((2{\omega }_L-2\omega )t+2\beta \right)+{\Lambda }^3\sin \left(({\omega }_L+3\omega )t+\beta \right)\\ +{\displaystyle \frac{3a^2\Lambda }{4}}\sin \left((3{\omega }_L+\omega )t+3\beta \right)+{\displaystyle \frac{3a^2\Lambda }{4}}\sin \left((3{\omega }_L-\omega )t+3\beta \right)+{\Lambda }^3\sin \left(({\omega }_L-3\omega )t+\beta \right)\end{array}\right\}\\ +{\displaystyle \frac{{\alpha }_4\epsilon }{{\omega }_L}}\left\{\begin{array}{l}-{\displaystyle \frac{a^2{\omega }_L}{4}}\cos ({\omega }_{L}t+\beta )+{\displaystyle \frac{a^2{\omega }_L}{4}}\cos (3{\omega }_{L}t+3\beta )-\Lambda a{\omega }_L\cos (\omega t)\\ +{\displaystyle \frac{a\Lambda }{2}}\left({\omega }_L+\omega \right)\cos \left((2{\omega }_L+\omega )t+2\beta \right)+{\displaystyle \frac{a\Lambda }{2}}\left({\omega }_L-\omega \right)\cos \left((2{\omega }_L-\omega )t+2\beta \right)\\ \omega {\Lambda }^2\cos \left((2{\omega }_L+\omega )t+2\beta \right)-\omega {\Lambda }^2\cos \left((2{\omega }_L-\omega )t+2\beta \right)\end{array}\right\}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}(52)\\ {\displaystyle \frac{{\alpha }_5\epsilon }{{\omega }_L}}\left\{\begin{array}{l}-\left(a{\omega }_L{\Lambda }^2+{\displaystyle \frac{a^3{\omega }_L}{8}}\right)-a{\omega }_L{\Lambda }^2\left(\cos (2\omega t)+\cos (2{\omega }_{L}t+2\beta )\right)+{\displaystyle \frac{a^3{\omega }_L}{8}}\cos (4{\omega }_{L}t+4\beta )\\ +\left(\omega {\Lambda }^3-{\displaystyle \frac{a^2{\omega }_L\Lambda }{2}}+{\displaystyle \frac{a^2\omega \Lambda }{4}}\right)\cos \left(({\omega }_L+\omega )t+\beta \right)-\left(\omega {\Lambda }^3+{\displaystyle \frac{a^2{\omega }_L\Lambda }{2}}+{\displaystyle \frac{a^2\omega \Lambda }{4}}\right)\cos \left(({\omega }_L-\omega )t+\beta \right)\\ +\left({\displaystyle \frac{a^2{\omega }_L\Lambda }{2}}+{\displaystyle \frac{a^2\omega \Lambda }{4}}\right)\cos \left((3{\omega }_L+\omega )t+3\beta \right)+\left({\displaystyle \frac{a^2{\omega }_L\Lambda }{2}}-{\displaystyle \frac{a^2\omega \Lambda }{4}}\right)\cos \left((3{\omega }_L-\omega )t+3\beta \right)\\ +\left(a\omega {\Lambda }^2+{\displaystyle \frac{a{\omega }_L{\Lambda }^2}{2}}\right)\cos \left((2{\omega }_L+2\omega )t+2\beta \right)+\left(-a\omega {\Lambda }^2+{\displaystyle \frac{a{\omega }_L{\Lambda }^2}{2}}\right)\cos \left((2{\omega }_L-2\omega )t+2\beta \right)\\ +\omega {\Lambda }^3\cos \left(({\omega }_L+3\omega )t+\beta \right)-\omega {\Lambda }^3\cos \left(({\omega }_L-3\omega )t+\beta \right)\end{array}\right\}\end{array}$$

By substituting (52) into (49), we obtain:

$$\begin{array}{l}a\dot{\beta }={\displaystyle \frac{\epsilon }{{\omega }_L}}\left\{-\mu a{\omega }_L\left(1-\sin (2{\omega }_{L}t+2\beta )\right)-2\mu \Lambda \omega \sin \left(({\omega }_L+\omega )t+\beta \right)+2\mu \Lambda \omega \sin \left(({\omega }_L-\omega )t+\beta \right)\right\}\\ +{\displaystyle \frac{{\alpha }_2\epsilon }{{\omega }_L}}\left\{\begin{array}{l}\left(2{\Lambda }^2+{\displaystyle \frac{3a^2}{4}}\right)cos({\omega }_{L}t+\beta )+2a\Lambda \cos (\omega t)+{\displaystyle \frac{a^2}{4}}\cos (3{\omega }_{L}t+3\beta )+{\Lambda }^2\cos \left(({\omega }_L+2\omega )t+\beta \right)\\ +{\Lambda }^2\cos \left(({\omega }_L-2\omega )t+\beta \right)+a\Lambda \cos \left((2{\omega }_L+\omega )t+2\beta \right)+a\Lambda \cos \left((2{\omega }_L-\omega )t+2\beta \right)\end{array}\right\}\\ +{\displaystyle \frac{{\alpha }_3\epsilon }{{\omega }_L}}\left\{\begin{array}{l}\left(3a{\Lambda }^2+{\displaystyle \frac{3a^3}{8}}\right)+\left(3a{\Lambda }^2+{\displaystyle \frac{a^3}{2}}\right)\cos (2{\omega }_{L}t+2\beta )+{\displaystyle \frac{a^3}{8}}\cos (4{\omega }_{L}t+4\beta )+3a{\Lambda }^2\cos (2\omega t)\\ \left(3{\Lambda }^3+{\displaystyle \frac{9a^2\Lambda }{4}}\right)\cos \left(({\omega }_L+\omega )t+\beta \right)+\left(3{\Lambda }^3+{\displaystyle \frac{9a^2\Lambda }{4}}\right)\cos \left(({\omega }_L-\omega )t+\beta \right)\\ +{\displaystyle \frac{3a^2\Lambda }{2}}\cos \left((2{\omega }_L+2\omega )t+2\beta \right)+{\displaystyle \frac{3a^2\Lambda }{2}}\cos \left((2{\omega }_L-2\omega )t+2\beta \right)+{\Lambda }^3\cos \left(({\omega }_L+3\omega )t+\beta \right)\\ +{\displaystyle \frac{3a^2\Lambda }{4}}\cos \left((3{\omega }_L+\omega )t+3\beta \right)+{\displaystyle \frac{3a^2\Lambda }{4}}\cos \left((3{\omega }_L-\omega )t+3\beta \right)+{\Lambda }^3\cos \left(({\omega }_L-3\omega )t+\beta \right)\end{array}\right\}\\ +{\displaystyle \frac{{\alpha }_4\epsilon }{{\omega }_L}}\left\{\begin{array}{l}-{\displaystyle \frac{a^2{\omega }_L}{4}}\sin ({\omega }_{L}t+\beta )-{\displaystyle \frac{a^2{\omega }_L}{4}}\sin (3{\omega }_{L}t+3\beta )-\Lambda a{\omega }_L\sin (\omega t)\\ -{\displaystyle \frac{a\Lambda }{2}}\left({\omega }_L+\omega \right)\sin \left((2{\omega }_L+\omega )t+2\beta \right)-{\displaystyle \frac{a\Lambda }{2}}\left({\omega }_L-\omega \right)\sin \left((2{\omega }_L-\omega )t+2\beta \right)\\ -\omega {\Lambda }^2\sin \left((2{\omega }_L+\omega )t+2\beta \right)+\omega {\Lambda }^2\sin \left((2{\omega }_L-\omega )t+2\beta \right)\end{array}\right\}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}(53)\\ {\displaystyle \frac{{\alpha }_5\epsilon }{{\omega }_L}}\left\{\begin{array}{l}-\left(a{\omega }_L{\Lambda }^2+{\displaystyle \frac{a^3{\omega }_L}{4}}\right)\sin (2{\omega }_{L}t+2\beta )-2a\omega {\Lambda }^2\sin (2\omega t)-{\displaystyle \frac{a^3{\omega }_L}{8}}\sin (4{\omega }_{L}t+4\beta )\\ -\left(\omega {\Lambda }^3+{\displaystyle \frac{a^2{\omega }_L\Lambda }{2}}+{\displaystyle \frac{3a^2\omega \Lambda }{4}}\right)\sin \left(({\omega }_L+\omega )t+\beta \right)+\left(\omega {\Lambda }^3-{\displaystyle \frac{a^2{\omega }_L\Lambda }{2}}+{\displaystyle \frac{3a^2\omega \Lambda }{4}}\right)\sin \left(({\omega }_L-\omega )t+\beta \right)\\ -\left({\displaystyle \frac{a^2{\omega }_L\Lambda }{2}}+{\displaystyle \frac{a^2\omega \Lambda }{4}}\right)\sin \left((3{\omega }_L+\omega )t+3\beta \right)-\left({\displaystyle \frac{a^2{\omega }_L\Lambda }{2}}-{\displaystyle \frac{a^2\omega \Lambda }{4}}\right)\sin \left((3{\omega }_L-\omega )t+3\beta \right)\\ -\left(a\omega {\Lambda }^2+{\displaystyle \frac{a{\omega }_L{\Lambda }^2}{2}}\right)\sin \left((2{\omega }_L+2\omega )t+2\beta \right)-\left(-a\omega {\Lambda }^2+{\displaystyle \frac{a{\omega }_L{\Lambda }^2}{2}}\right)\sin \left((2{\omega }_L-2\omega )t+2\beta \right)\\ -\omega {\Lambda }^3\sin \left(({\omega }_L+3\omega )t+\beta \right)+\omega {\Lambda }^3\sin \left(({\omega }_L-3\omega )t+\beta \right)\end{array}\right\}\end{array}$$

In secondary resonance, there are two cases: superharmonic resonance and subharmonic resonance.

2.2.1. Superharmonic Resonance

In this case, we have: $\omega \approx {\displaystyle \frac{1}{3}}{\omega }_L\hspace{0.33em}et\hspace{0.33em}\omega \approx {\displaystyle \frac{1}{2}}{\omega }_L$.

First case: $\hspace{0.33em}\omega \approx {\displaystyle \frac{1}{3}}{\omega }_L,\hspace{1em}\Lambda ={\displaystyle \frac{9F_1}{{\omega }_L^2}}$The expressions for$\hspace{0.33em}\sin \left(({\omega }_L-3\omega )t+\beta \right)$ and $\cos \left(({\omega }_L-3\omega )t+\beta \right)$ vary slowly with time. Equations (52) and (53) then become:

$$\left\{\begin{array}{l}\dot{a}+\mu a\epsilon +{\displaystyle \frac{{\alpha }_5\epsilon a^3}{8}}+{\alpha }_5\epsilon a{\Lambda }^2={\displaystyle \frac{{\alpha }_3\epsilon {\Lambda }^3}{{\omega }_L}}\sin \left(({\omega }_L-3\omega )t+\beta \right)-{\displaystyle \frac{{\alpha }_5\epsilon {\Lambda }^3}{3{\omega }_L}}\cos \left(({\omega }_L-3\omega )t+\beta \right)\\ a\dot{\beta }-{\displaystyle \frac{3{\alpha }_3\epsilon a^3}{8{\omega }_L}}-{\displaystyle \frac{3{\alpha }_3\epsilon a{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{{\alpha }_3\epsilon {\Lambda }^3}{{\omega }_L}}\cos \left(({\omega }_L-3\omega )t+\beta \right)+{\displaystyle \frac{{\alpha }_5\epsilon {\Lambda }^3}{3{\omega }_L}}\sin \left(({\omega }_L-3\omega )t+\beta \right)\end{array}\right.$$

(54)

By setting $\gamma =({\omega }_L-3\omega )t+\beta$ and $3\omega ={\omega }_L+\epsilon \sigma$, the system (54) becomes:

$$\left\{\begin{array}{l}\dot{a}+\mu a\epsilon +{\displaystyle \frac{{\alpha }_5\epsilon a^3}{8}}+{\alpha }_5\epsilon a{\Lambda }^2={\displaystyle \frac{{\alpha }_3\epsilon {\Lambda }^3}{{\omega }_L}}\sin (\gamma )-{\displaystyle \frac{{\alpha }_5\epsilon {\Lambda }^3}{3{\omega }_L}}\cos (\gamma )\\ a\dot{\gamma }+\epsilon a\sigma -{\displaystyle \frac{3{\alpha }_3\epsilon a^3}{8{\omega }_L}}-{\displaystyle \frac{3{\alpha }_3\epsilon a{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{{\alpha }_3\epsilon {\Lambda }^3}{{\omega }_L}}\cos (\gamma )+{\displaystyle \frac{{\alpha }_5\epsilon {\Lambda }^3}{3{\omega }_L}}\sin (\gamma )\end{array}\right.$$

(55)

If the system tends toward a steady state,$\dot{a}=0$ and $\dot{\gamma }=0,$ we have:

$$\left\{\begin{array}{l}\mu +{\displaystyle \frac{{\alpha }_5{a}^2}{8}}+{\alpha }_5{\Lambda }^2={\displaystyle \frac{{\alpha }_3{\Lambda }^3}{a{\omega }_L}}\sin (\gamma )-{\displaystyle \frac{{\alpha }_5{\Lambda }^3}{3a{\omega }_L}}\cos (\gamma )\\ \sigma -{\displaystyle \frac{3{\alpha }_3{a}^2}{8{\omega }_L}}-{\displaystyle \frac{3{\alpha }_3{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{{\alpha }_3{\Lambda }^3}{a{\omega }_L}}\cos (\gamma )+{\displaystyle \frac{{\alpha }_5{\Lambda }^3}{3a{\omega }_L}}\sin (\gamma )\end{array}\right.$$

(56)

From system (56), we obtain:

$${\left[{\displaystyle \frac{\omega }{{\omega }_L}}-\left({\displaystyle \frac{1}{3}}+{\displaystyle \frac{{\alpha }_3{a}^2}{8{\omega }_L^2}}+{\displaystyle \frac{{\alpha }_3{\Lambda }^2}{{\omega }_L^2}}\right)\right]}^2+{\left[{\displaystyle \frac{\mu }{3{\omega }_L}}+{\displaystyle \frac{{\alpha }_5{a}^2}{24{\omega }_L}}+{\displaystyle \frac{{\alpha }_5{\Lambda }^2}{3{\omega }_L}}\right]}^2={\left({\displaystyle \frac{{\alpha }_3{\Lambda }^3}{3a{\omega }_L^2}}\right)}^2+{\left({\displaystyle \frac{{\alpha }_5{\Lambda }^3}{9a{\omega }_L}}\right)}^2$$

(57)

Second case :$\hspace{0.33em}\omega \approx {\displaystyle \frac{1}{2}}{\omega }_L,\hspace{1em}\Lambda ={\displaystyle \frac{2F_1}{3{\omega }_L^2}}$The expressions for$\sin \left(({\omega }_L-2\omega )t+\beta \right)$and$\cos \left(({\omega }_L-2\omega )t+\beta \right)$vary slowly with time. Equations (52) and (53) then become:

$$\left\{\begin{array}{l}\dot{a}+\mu a\epsilon +{\displaystyle \frac{{\alpha }_5\epsilon a^3}{8}}+{\alpha }_5\epsilon a{\Lambda }^2={\displaystyle \frac{{\alpha }_2\epsilon {\Lambda }^2}{{\omega }_L}}\sin \left(({\omega }_L-2\omega )t+\beta \right)-{\displaystyle \frac{{\alpha }_4\epsilon {\Lambda }^2}{2{\omega }_L}}\cos \left(({\omega }_L-2\omega )t+\beta \right)\\ a\dot{\beta }-{\displaystyle \frac{3{\alpha }_3\epsilon a^3}{8{\omega }_L}}-{\displaystyle \frac{3{\alpha }_3\epsilon a{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{{\alpha }_2\epsilon {\Lambda }^2}{{\omega }_L}}\cos \left(({\omega }_L-2\omega )t+\beta \right)+{\displaystyle \frac{{\alpha }_4\epsilon {\Lambda }^2}{2{\omega }_L}}\cos \left(({\omega }_L-2\omega )t+\beta \right)\end{array}\right.$$

(58)

By setting $\hspace{0.17em}\gamma =({\omega }_L-2\omega )t+\beta \hspace{1em}et\hspace{1em}2\omega ={\omega }_L+\epsilon \sigma$, the system (58) becomes:

$$\left\{\begin{array}{l}\dot{a}+\mu a\epsilon +{\displaystyle \frac{{\alpha }_5\epsilon a^3}{8}}+{\alpha }_5\epsilon a{\Lambda }^2={\displaystyle \frac{{\alpha }_2\epsilon {\Lambda }^2}{{\omega }_L}}\sin (\gamma )-{\displaystyle \frac{{\alpha }_4\epsilon {\Lambda }^2}{2{\omega }_L}}\cos (\gamma )\\ a\dot{\gamma }+\epsilon a\sigma -{\displaystyle \frac{3{\alpha }_3\epsilon a^3}{8{\omega }_L}}-{\displaystyle \frac{3{\alpha }_3\epsilon a{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{{\alpha }_2\epsilon {\Lambda }^2}{{\omega }_L}}\cos (\gamma )+{\displaystyle \frac{{\alpha }_4\epsilon {\Lambda }^2}{2{\omega }_L}}\sin (\gamma )\end{array}\right.$$

(59)

If the system tends toward a steady state,$\dot{a}=0$ and $\dot{\gamma }=0,$ we have:

$$\left\{\begin{array}{l}\mu +{\displaystyle \frac{{\alpha }_5{a}^2}{8}}+{\alpha }_5{\Lambda }^2={\displaystyle \frac{{\alpha }_2{\Lambda }^2}{{\omega }_L}}\sin (\gamma )-{\displaystyle \frac{{\alpha }_4{\Lambda }^2}{2{\omega }_L}}\cos (\gamma )\\ \sigma -{\displaystyle \frac{3{\alpha }_3{a}^2}{8{\omega }_L}}-{\displaystyle \frac{3{\alpha }_3{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{{\alpha }_2{\Lambda }^2}{{\omega }_L}}\cos (\gamma )+{\displaystyle \frac{{\alpha }_4{\Lambda }^2}{2{\omega }_L}}\sin (\gamma )\end{array}\right.$$

(60)

From system (60), we obtain:

$${\left[{\displaystyle \frac{\omega }{{\omega }_L}}-\left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{3{\alpha }_3{a}^2}{16{\omega }_L^2}}+{\displaystyle \frac{3{\alpha }_3{\Lambda }^2}{2{\omega }_L^2}}\right)\right]}^2+{\left[{\displaystyle \frac{\mu }{2{\omega }_L}}+{\displaystyle \frac{{\alpha }_5{a}^2}{16{\omega }_L}}+{\displaystyle \frac{{\alpha }_5{\Lambda }^2}{2{\omega }_L}}\right]}^2={\left({\displaystyle \frac{{\alpha }_2{\Lambda }^2}{2a{\omega }_L}}\right)}^2+{\left({\displaystyle \frac{{\alpha }_4{\Lambda }^2}{4a{\omega }_L}}\right)}^2$$

(61)

2.2.2. Subharmonic Resonance

In this case, we have: $\hspace{0.33em}\omega \approx 3{\omega }_L\hspace{0.33em}et\hspace{0.33em}\omega \approx 2{\omega }_L$First case:$\hspace{0.33em}\omega \approx 3{\omega }_L,\hspace{1em}\Lambda ={\displaystyle \frac{-F_1}{16{\omega }_L^2}}$The expressions for$\sin \left((3{\omega }_L-\omega )t+\beta \right)$ and $\cos \left((3{\omega }_L-\omega )t+\beta \right)$vary slowly with time. Equations (52) and (53) then become:

$$\left\{\begin{array}{l}\dot{a}+\mu a\epsilon +{\displaystyle \frac{{\alpha }_5\epsilon a^3}{8}}+{\alpha }_5\epsilon a{\Lambda }^2={\displaystyle \frac{3{\alpha }_3\epsilon a^2\Lambda }{4{\omega }_L}}\sin \left((3{\omega }_L-\omega )t+3\beta \right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\displaystyle \frac{{\alpha }_5\epsilon a^2\Lambda }{2}}\cos \left((3{\omega }_L-\omega )t+3\beta \right)\\ a\dot{\beta }-{\displaystyle \frac{3{\alpha }_3\epsilon a^3}{8{\omega }_L}}-{\displaystyle \frac{3{\alpha }_3\epsilon a{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{3{\alpha }_3\epsilon a^2\Lambda }{4{\omega }_L}}\cos \left((3{\omega }_L-\omega )t+3\beta \right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \frac{{\alpha }_5\epsilon a^2\Lambda }{2}}\sin \left((3{\omega }_L-\omega )t+3\beta \right)\end{array}\right.$$

(62)

By setting$\hspace{0.17em}\gamma =({\omega }_L-\omega )t+3\beta \hspace{0.33em}et\hspace{0.33em}\omega =3{\omega }_L+\epsilon \sigma$, the system (62) becomes:

$$\left\{\begin{array}{l}\dot{a}+\mu a\epsilon +{\displaystyle \frac{{\alpha }_5\epsilon a^3}{8}}+{\alpha }_5\epsilon a{\Lambda }^2={\displaystyle \frac{3{\alpha }_3\epsilon a^2\Lambda }{4{\omega }_L}}\sin (\gamma )-{\displaystyle \frac{{\alpha }_5\epsilon a^2\Lambda }{4}}\cos (\gamma )\\ a\dot{\gamma }+\epsilon a\sigma -{\displaystyle \frac{9{\alpha }_3\epsilon a^3}{8{\omega }_L}}-{\displaystyle \frac{9{\alpha }_3\epsilon a{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{9{\alpha }_3\epsilon a^2\Lambda }{4{\omega }_L}}\cos (\gamma )+{\displaystyle \frac{3{\alpha }_5\epsilon a^2\Lambda }{4}}\sin (\gamma )\end{array}\right.$$

(63)

If the system tends toward a steady state,$\dot{a}=0$ and $\dot{\gamma }=0,$ we have:

$$\left\{\begin{array}{l}\mu +{\displaystyle \frac{{\alpha }_5{a}^2}{8}}+{\alpha }_5{\Lambda }^2={\displaystyle \frac{3{\alpha }_{3}a\Lambda }{4{\omega }_L}}\sin (\gamma )-{\displaystyle \frac{{\alpha }_{5}a\Lambda }{4}}\cos (\gamma )\\ \sigma -{\displaystyle \frac{9{\alpha }_3{a}^2}{8{\omega }_L}}-{\displaystyle \frac{9{\alpha }_3{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{9{\alpha }_{3}a\Lambda }{4{\omega }_L}}\cos (\gamma )+{\displaystyle \frac{3{\alpha }_{5}a\Lambda }{4}}\sin (\gamma )\end{array}\right.$$

(64)

From system (64), we obtain:

$${\left[{\displaystyle \frac{\omega }{{\omega }_L}}-\left(3+{\displaystyle \frac{9{\alpha }_3{a}^2}{8{\omega }_L^2}}+{\displaystyle \frac{9{\alpha }_3{\Lambda }^2}{{\omega }_L^2}}\right)\right]}^2+{\left[{\displaystyle \frac{3\mu }{{\omega }_L}}+{\displaystyle \frac{3{\alpha }_5{a}^2}{8{\omega }_L}}+{\displaystyle \frac{3{\alpha }_5{\Lambda }^2}{{\omega }_L}}\right]}^2={\left({\displaystyle \frac{9{\alpha }_{3}a\Lambda }{4{\omega }_L^2}}\right)}^2+{\left({\displaystyle \frac{3{\alpha }_{5}a\Lambda }{4{\omega }_L}}\right)}^2$$

(65)

Second case:$\hspace{0.33em}\omega \approx 2{\omega }_L,\hspace{1em}\Lambda ={\displaystyle \frac{-F_1}{6{\omega }_L^2}}$The expressions for$\sin \left((2{\omega }_L-\omega )t+2\beta \right)$ and $\cos \left((2{\omega }_L-\omega )t+2\beta \right)$vary slowly with time. Equations (52) and (53) then become:

$$\left\{\begin{array}{l}\dot{a}+\mu a\epsilon +{\displaystyle \frac{{\alpha }_5\epsilon a^3}{8}}+{\alpha }_5\epsilon a{\Lambda }^2={\displaystyle \frac{{\alpha }_2\epsilon a\Lambda }{{\omega }_L}}\sin \left((2{\omega }_L-\omega )t+2\beta \right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\displaystyle \frac{{\alpha }_4\epsilon a\Lambda }{2}}\cos \left((2{\omega }_L-\omega )t+2\beta \right)\\ a\dot{\beta }-{\displaystyle \frac{3{\alpha }_3\epsilon a^3}{8{\omega }_L}}-{\displaystyle \frac{3{\alpha }_3\epsilon a{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{{\alpha }_2\epsilon a\Lambda }{{\omega }_L}}\cos \left((2{\omega }_L-\omega )t+2\beta \right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \frac{{\alpha }_4\epsilon a\Lambda }{2}}\sin \left((2{\omega }_L-\omega )t+2\beta \right)\end{array}\right.$$

(66)

By setting$\hspace{0.17em}\gamma =(2{\omega }_L-\omega )t+2\beta ,\hspace{1em}\omega =2{\omega }_L+\epsilon \sigma$, the system (66) gives:

$$\left\{\begin{array}{l}\dot{a}+\mu a\epsilon +{\displaystyle \frac{{\alpha }_5\epsilon a^3}{8}}+{\alpha }_5\epsilon a{\Lambda }^2={\displaystyle \frac{{\alpha }_2\epsilon a\Lambda }{{\omega }_L}}\sin (\gamma )-{\displaystyle \frac{{\alpha }_4\epsilon a\Lambda }{2}}\cos (\gamma )\\ a\dot{\gamma }+\epsilon a\sigma -{\displaystyle \frac{3{\alpha }_3\epsilon a^3}{4{\omega }_L}}-{\displaystyle \frac{6{\alpha }_3\epsilon a{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{2{\alpha }_2\epsilon a\Lambda }{{\omega }_L}}\cos (\gamma )+{\alpha }_4\Lambda \sin (\gamma )\end{array}\right.$$

(67)

In a steady state, we obtain:

$$\left\{\begin{array}{l}\mu +{\displaystyle \frac{{\alpha }_5{a}^2}{8}}+{\alpha }_5{\Lambda }^2={\displaystyle \frac{{\alpha }_2\Lambda }{{\omega }_L}}\sin (\gamma )-{\displaystyle \frac{{\alpha }_4\Lambda }{2}}\cos (\gamma )\\ \sigma -{\displaystyle \frac{3{\alpha }_3{a}^2}{4{\omega }_L}}-{\displaystyle \frac{6{\alpha }_3{\Lambda }^2}{{\omega }_L}}={\displaystyle \frac{2{\alpha }_2\Lambda }{{\omega }_L}}\cos (\gamma )+{\alpha }_4\Lambda \sin (\gamma )\end{array}\right.$$

(68)

From system (68), we obtain:

$${\left[{\displaystyle \frac{\omega }{{\omega }_L}}-\left(2+{\displaystyle \frac{3{\alpha }_3{a}^2}{4{\omega }_L^2}}+{\displaystyle \frac{6{\alpha }_3{\Lambda }^2}{{\omega }_L^2}}\right)\right]}^2+{\left[{\displaystyle \frac{2\mu }{{\omega }_L}}+{\displaystyle \frac{{\alpha }_5{a}^2}{4{\omega }_L}}+{\displaystyle \frac{2{\alpha }_5{\Lambda }^2}{{\omega }_L}}\right]}^2={\left({\displaystyle \frac{2{\alpha }_2\Lambda }{{\omega }_L^2}}\right)}^2+{\left({\displaystyle \frac{{\alpha }_4\Lambda }{{\omega }_L}}\right)}^2$$

(69)

Briefly, the structural modeling allowed us to represent the nonlinear dynamic behavior of the beam in the form of a differential equation. Then, the application of the method of averaging to the dynamic equation enabled the calculation of the different resonances. In the following section, we present the results of the numerical simulations and their discussion.

3. Results and Discussion

3.1. Boundary Conditions and Beam Properties

In this study, three boundary conditions were used:

- S-S: simply supported beam;

- C-S: clamped-simply supported beam;

- C-C: clamped-clamped beam.

The geometric and material properties are presented in Table 1. The numerical values corresponding to the different boundary conditions are given in Table 2.

Using the feedback control parameters, the μᵢ and Cᵢ coefficients in Table 2 are given by:

$$\begin{array}{l}2\mu ={\mu }_1{G}_d\hspace{0.17em};\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{1em}\hspace{0.33em}{\alpha }_2=C_2+{\mu }_2{G}_P\hspace{0.17em};\hspace{1em}\hspace{1em}\hspace{0.33em}{\alpha }_3=C_3+{\mu }_3{G}_P\hspace{0.17em};\\ {\omega }_L^2=C_1+{\mu }_1{G}_P\hspace{0.17em};\hspace{1em}\hspace{0.33em}\hspace{0.17em}\hspace{0.33em}{\alpha }_4={\mu }_4{G}_d\hspace{0.17em};\hspace{0.17em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\alpha }_5={\mu }_5{G}_d.\end{array}$$

$$\begin{array}{l}{C}_1={\displaystyle \frac{1}{M}}\left(E_c{I}_c+2c_{11}^{*}\left(I_s+S_s{z}_s^2\right)+{\displaystyle \frac{{(ES)}_{pe}}{S_s}}\left(2I_s+S_s{z}_s^2\right)\right){\displaystyle {\int }_0^L{\phi }_{xxxx}(x)\phi (x)dx}\hspace{0.17em};\\ C_2={\displaystyle \frac{{(ES)}_{pe}{z}_s}{M}}\left({\displaystyle \frac{1}{L}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)\phi (x)dx{\displaystyle {\int }_0^L{\phi }_{xx}(x)dx}-{\displaystyle {\int }_0^L\left\{{\left({\phi }_{xx}(x)\right)}^2+{\phi }_x(x){\phi }_{xxx}(x)\right\}}}{\phi }_x(x)dx\right)\hspace{0.17em};\\ C_3=-{\displaystyle \frac{E_c{I}_c+2c_{11}^{*}{S}_s+{(ES)}_{pe}}{2ML}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)\phi (x)dx{\displaystyle {\int }_0^L{\left({\phi }_x(x)\right)}^2}}dx\hspace{0.17em};\\ {\mu }_1={\displaystyle \frac{{(ES)}_{pe}{z}_s^2}{M}}{\displaystyle {\int }_0^L{\phi }_{xxxx}(x)\phi (x)dx}\hspace{0.17em};\\ {\mu }_2=-{\displaystyle \frac{{(ES)}_{pe}{z}_s}{M}}\left({\displaystyle \frac{1}{L}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)\phi (x)dx{\displaystyle {\int }_0^L{\phi }_{xx}(x)dx}+{\displaystyle {\int }_0^L\left\{{\left({\phi }_{xx}(x)\right)}^2+{\phi }_x(x){\phi }_{xxx}(x)\right\}}}{\phi }_x(x)dx\right)\hspace{0.17em};\\ {\mu }_3={\displaystyle \frac{{(ES)}_{pe}}{2ML}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)\phi (x)dx{\displaystyle {\int }_0^L{\left({\phi }_x(x)\right)}^2}}dx\hspace{0.17em};\\ {\mu }_4=-{\displaystyle \frac{{(ES)}_{pe}{z}_s}{M}}\left({\displaystyle \frac{1}{L}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)\phi (x)dx{\displaystyle {\int }_0^L{\phi }_{xx}(x)dx}+2{\displaystyle {\int }_0^L\left\{{\left({\phi }_{xx}(x)\right)}^2+{\phi }_x(x){\phi }_{xxx}(x)\right\}}}{\phi }_x(x)dx\right)\hspace{0.17em};\\ {\mu }_5={\displaystyle \frac{{(ES)}_{pe}}{ML}}{\displaystyle {\int }_0^L{\phi }_{xx}(x)\phi (x)dx{\displaystyle {\int }_0^L{\left({\phi }_x(x)\right)}^2}}dx\hspace{0.17em};\\ F_1={\displaystyle \frac{1}{M}}{\displaystyle {\int }_0^{L}f(x)\phi (x)dx\hspace{0.17em};\hspace{1em}}M={(\rho S)}_{*}{\displaystyle {\int }_0^L{\left(\phi (x)\right)}^{2}dx\hspace{0.17em};\hspace{1em}{(ES)}_{pe}=S_s\frac{{(e_{31}^{*})}^2}{{\in }_{33}^{*}}}.\end{array}$$

3.2. Primary Resonance

In all the presented figures, the frequency is normalized with respect to the natural frequency of the studied beam, and the beam is uniformly excited.

Figure 2 compares three analytical methods: the method of averaging (MA), used in this work, the harmonic balance method (HBM) [16,17,18,19,20,29], and the multiple scales method (MSM) [12,30], which have been used in the literature. The curves obtained using the method of averaging and the multiple scales method overlap, as shown in the figure, whereas a slight deviation is observed with those obtained by the harmonic balance method. This is likely due to the lower accuracy of the latter method. The method of averaging can thus be considered a reliable technique for vibration control in structures with good accuracy.

Figure 3 and Figure 4 illustrate typical nonlinear phenomena. When G_p = 20, the stiffness increases, and the system is said to be hardening, which results in a frequency response curve leaning toward higher frequencies. When G_p = 35, the stiffness decreases, and the system is said to be softening, leading to a frequency response curve inclined toward lower frequencies. In both cases, multiple solutions may exist for the same excitation frequency. This gives rise to jump phenomena, depending on the direction of the frequency sweep.

Figure 5 presents the effects of boundary conditions on the amplitude-frequency response of the structure. The S-S beam is highly sensitive to changes in the G_p parameter compared to the C-S beam, and more sensitive than the C-C beam. The effects of geometric nonlinearities are significantly reduced when the value of this parameter is increased.

Figure 6 shows the free responses, i.e., when the external force is zero. However, depending on the initial shape of the structure at rest, some behaviors can be softening, while others are hardening. This behavior may depend on the static stiffness of the structure around its equilibrium point and the presence of nonlinearities.

In Figure 7, the proportional control gain G_p is fixed at 20 while the velocity control gain G_d is varied. As G_d increases, the vibration amplitudes decrease, and the system remains in the hardening regime. To observe the opposite effect, i.e., the softening behavior, one simply needs to increase the value of G_p, as shown in Figure 8. The hardening-softening transition is independent of G_d, since increasing or decreasing its value always results in a hardening behavior.

3.3. Secondary Resonance

3.3.1. Superharmonic Resonance

Figure 9 and Figure 10 show the effects of nonlinearities on beam S-S.

For $\omega \approx 1/3{\omega }_L$, on the stiffening side, the amplitudes decrease as G_p increases; however, on the softening side, the amplitudes increase with increasing G_p.

For $\omega \approx 1/2{\omega }_L$, the vibration amplitudes decrease on both sides as G_p increases. In this particular case, the transition can be controlled by the parameter G_p and the resonance shifts.

For C-S beam in Figure 11, nonlinear behaviors occur around the structure’s natural mode. In Figure 12, the excitation effect is negligible on C-C beam. These behaviors result from the presence of clamping in the structure.

3.3.2. Subharmonic Resonance

In Figure 13, the behavior of beam S-S exhibits hysteresis. This suggests that the nonlinear relationship between amplitude and frequency is due to mechanical losses within the structure. These losses, having become a form of energy, dissipate within the structure; therefore, control gains help reduce these losses by decreasing the area of the amplitude-frequency hysteresis loop. On the stiffening side for fixed G_d, this area can be reduced by increasing the gain G_p, whereas on the softening side, the area is instead reduced by decreasing the gain. The same phenomenon is observed in Figure 14 and Figure 15, with G_p fixed and G_d varied. For all values of G_d, the beam exhibits a stiffening behavior, and increasing its value reduces the hysteresis area. Thus, by minimizing this energy, the operational range of the beam can be optimized. These behaviors are significant for energy harvesting. Similarly, in Figure 16 and Figure 17, second-order nonlinear effects are observed, and one can also control the vibration amplitudes as well as the stiffening-softening transition. As for Figure 18, it shows that the influence of excitation is negligible on beam C-C.

4. Conclusions

This study aimed to explore the method of averaging for active control of nonlinear vibrations in a sandwich beam. Based on the results obtained and supported by findings in the literature, this method proves to be highly effective for active vibration control and may also be very relevant for energy harvesting. Three boundary conditions were highlighted: simply supported beam, clamped-simply supported beam, and clamped-clamped beam. In general, the dynamic behavior of the beam can be controlled through control gains for all these boundary conditions using a retractable control law. The stiffening-softening transition can be managed via the control gain G_p, while amplitude reduction is achieved through the parameter G_d.