Nonlinear Vibration Control of Piezoelectric-Elastic-Piezoelectric Sandwich Beam

Barthelemy Zra Mha, Guy Edgar Ntamack

Mechanics, Materials and Acoustics Group, Department of Physics, Faculty of Science, University of Ngaoundéré, Ngaoundéré, Cameroon.

DOI: 10.4236/wjm.2025.155005                

Abstract

This work presents the active control of nonlinear vibrations of a piezoelectric-elastic-piezoelectric sandwich beam, subjected to transverse excitation while neglecting axial displacement effects. By using a structure with piezoelectric actuators and sensors, and taking into account geometric nonlinearities, a nonlinear vibration control model was obtained through a feedback control law, which is a proportional-derivative controller. The dynamic equation of the structure is derived by applying the variational principle and Hamilton’s principle. This equation is solved under primary and secondary resonance by adopting the method of averaging as a perturbation scheme and Galerkin’s approximation. The simulation results of amplitude-frequency responses are presented and discussed for different values of control gains and for three boundary conditions. Our results are in good agreement with those obtained by other methods.

Keywords

Piezoelectric, Sandwich Beam, Nonlinear Vibrations, Active Control, Method of Averaging

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]-[3]. In many industrial and defense applications, noise and vibrations represent a major challenge. Conventional vibration attenuation methods based on passive devices, such as tuned mass dampers or viscoelastic treatments, exhibit limited efficiency at low frequencies. At such frequencies, the kinetic energy involved is relatively small, which requires very large auxiliary masses to achieve noticeable effects, often making the solution impractical. In addition, the intrinsic dissipation of damping materials decreases when the strain rate is low, thereby reducing their ability to absorb vibratory energy. Finally, passive devices have a narrow bandwidth and are highly sensitive to mistuning: even slight shifts in the natural frequency of the structure, caused, for instance, by temperature or load variations, can significantly degrade their performance. These limitations explain why passive control strategies are generally ineffective at low frequencies and highlight the relevance of active control methods [4]-[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 electromechanical 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 have 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]-[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]-[20].

2. Mathematical Modeling

2.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]-[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.

Figure 1. Piezoelectric-elastic-piezoelectric sandwich beam [14].

2.2. Kinematic Description of the Beam

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

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

In the theory of beams undergoing large deformations, the Green-Lagrange strain tensor is considered without linearization and is 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:

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

2.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].

$\{\begin{array}{l}\sigma =c\epsilon -e^{t}E\\ D=e\epsilon +ϵE\end{array}$ (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}{ϵ}_{11}& 0& 0\\ 0& {ϵ}_{22}& 0\\ 0& 0& {ϵ}_{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:

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

It follows that:

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

With: ${ϵ}_{33}^{*}={ϵ}_{33}+\frac{e_{33}^2}{c_{33}};e_{31}^{*}=e_{31}-\frac{c_{13}}{c_{33}}{e}_{33};c_{11}^{*}=c_{11}-\frac{c_{13}^2}{c_{33}}$ .

2.4. Feedback Control Law

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

$\frac{\partial D_3}{\partial z}=0$ (7)

Using the boundary conditions, $D\left(z_{-}\right)=0$ or $D\left(z_{+}\right)=0$ , we have:

$D_3\left(z\right)=0$ (8)

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

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

Since the electric field is derived from a potential, we have:

$E_3=-\frac{\partial \varphi }{\partial z}$ (10)

Consequently, the potential difference is given by:

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

With: $z_{-}=z_i-\frac{h_i}{2}$ and ${z_{+}=z_i+\frac{h_i}{2}|}_{i=a,s}$ .

From Equations (9) and (11), we have:

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

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

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

The actuator potential ${\varphi }_a\left(x\right)$ depends on the sensor output potential ${\varphi }_s\left(x\right)$ through a proportional-derivative control law described as follows:

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

In this formulation, G_p represents the proportional gain, which acts directly on the sensor displacement and thus modifies the effective stiffness perceived by the structure. By adjusting G_p, the stiffening or softening behavior of the system can be influenced, and the effective natural frequency can be shifted. G_d, on the other hand, is the derivative (or velocity) gain, which acts on the sensor velocity and contributes to the effective damping of the system, reducing vibration amplitudes and the hysteresis associated with nonlinear responses.

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

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

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

The potentials ${\varphi }_a$ and ${\varphi }_s$ are independent of z and $z_s=\left(h_c+h_s\right)/2$ et $z_a=-\left(h_c+h_s\right)/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].

2.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}{c}{\displaystyle \underset{V}{\int }{\sigma }_1\delta \epsilon \text{d}v}={\displaystyle \underset{V_s}{\int }{\sigma }_1\delta \epsilon \text{d}{v}_s}+{\displaystyle \underset{V_c}{\int }{\sigma }_1\delta \epsilon \text{d}{v}_c}+{\displaystyle \underset{V_a}{\int }{\sigma }_1\delta \epsilon \text{d}{v}_a}\\ ={\displaystyle {\int }_0^L\left(N\delta {\epsilon }_0+M\delta w_{,xx}\right)\text{d}x}\\ ={\displaystyle {\int }_0^L\left(F_X\delta u+F_Z\delta w\right)\text{d}x}-{\left(\rho S\right)}_{*}{\displaystyle {\int }_0^L\left(\ddot{u}\delta u+\ddot{w}\delta w\right)\text{d}x}\end{array}$ (17)

With: $\{\begin{array}{l}{\left(\rho S\right)}_{*}={\rho }_s{S}_s+{\rho }_c{S}_c+{\rho }_a{S}_a\\ \delta {\epsilon }_0=\delta u_{,x}+w_{,x}\delta w_{,x}\end{array}$

According to Hamilton’s principle, we have:

$\{\begin{array}{l}{N}_{xx}={\displaystyle \underset{s}{\int }{\sigma }_{xx}\text{d}s}=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}\text{d}s}=-EI\frac{{\partial }^{2}w}{\partial x^2}\end{array}$ (18)

If we integrate over the entire thickness and width, and assume that the piezoelectric layers are symmetrical ($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$ (19)

Using Equation (5), we obtain:

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

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

Using Equations (2), (13)-(16), we obtain:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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 (25) or (28) [12].

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

3. Solution Methodology

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

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

The Galerkin approximation [12] [28] is given as follows:

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

${\phi }_k\left(x\right)$ : are the vibration modes of the beam;

$q_k\left(t\right)$ : 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 (30) into Equation (28), integrating over the entire length, and omitting the indices since only one mode is considered, we obtain:

$\begin{array}{l}\ddot{q}\left(t\right)+2\mu \dot{q}\left(t\right)+{\omega }_L^{2}q\left(t\right)+{\alpha }_2{q}^2\left(t\right)+{\alpha }_3{q}^3\left(t\right)+{\alpha }_{4}q\left(t\right)\dot{q}\left(t\right)+{\alpha }_5{q}^2\left(t\right)\dot{q}\left(t\right)\\ =F_1\cos \left(\omega t\right)\end{array}$ (31)

With:

$\begin{array}{l}2\mu =\frac{{\left(ES\right)}_{pe}{G}_d{z}_s^2}{M}{\displaystyle {\int }_0^L{\phi }_{xxxx}\left(x\right)\phi \left(x\right)\text{d}x};\\ {\omega }_L^2=\frac{{\left(EI\right)}_{*}}{M}{\displaystyle {\int }_0^L{\phi }_{xxxx}\left(x\right)\phi \left(x\right)\text{d}x};\\ {\alpha }_2=\frac{B_N}{ML}{\displaystyle {\int }_0^L{\phi }_{xx}\left(x\right)\phi \left(x\right)\text{d}x}{\displaystyle {\int }_0^L{\phi }_{xx}\left(x\right)\text{d}x}-\frac{B_M}{M}{\displaystyle {\int }_0^L\left\{{\left({\phi }_{xx}\left(x\right)\right)}^2+{\phi }_x\left(x\right){\phi }_{xxx}\left(x\right)\right\}\phi \left(x\right)\text{d}x};\\ {\alpha }_3=-\frac{{\left(ES\right)}_{*}}{2ML}{\displaystyle {\int }_0^L{\phi }_{xx}\left(x\right)\phi \left(x\right)\text{d}x}{\displaystyle {\int }_0^L{\left({\phi }_x\left(x\right)\right)}^2\text{d}x};\\ {\alpha }_4=-\frac{{\left(ES\right)}_{pe}{G}_d{z}_s}{M}\left(\frac{1}{L}{\displaystyle {\int }_0^L{\phi }_{xx}\left(x\right)\phi \left(x\right)\text{d}x}{\displaystyle {\int }_0^L{\phi }_{xx}\left(x\right)\text{d}x}+2\left({\displaystyle {\int }_0^L{\left({\phi }_{xx}\left(x\right)\right)}^2+{\phi }_x\left(x\right){\phi }_{xxx}\left(x\right)\phi \left(x\right)\text{d}x}\right)\right);\\ {\alpha }_5=\frac{{\left(ES\right)}_{pe}{G}_d}{ML}{\displaystyle {\int }_0^L{\phi }_{xx}\left(x\right)\phi \left(x\right)\text{d}x}{\displaystyle {\int }_0^L{\left({\phi }_x\left(x\right)\right)}^2\text{d}x};\\ F_1=\frac{1}{M}{\displaystyle {\int }_0^{L}f\left(x\right)\phi \left(x\right)\text{d}x};\\ M={\left(\rho S\right)}_{*}{\displaystyle {\int }_0^L{\left({\phi }_x\left(x\right)\right)}^2\text{d}x}.\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 (31) will be carried out in the following, using the method of averaging.

3.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 \left(\omega t\right)\right]$ (32)

If $\epsilon =0$ , the general solution of Equation (31) is given by:

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

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

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

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

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

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

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

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

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

By substituting (34), (35), and (37) into (32), we obtain:

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

By using (36) and (38), we have:

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

By substituting (39) into (36), we obtain:

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

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

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

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

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

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

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

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

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

3.2. Secondary Resonance

In the case of secondary resonance, Equation (31) 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 \left(\omega t\right)$ (45)

If $\epsilon =0$ and using the principle of superposition, the general solution of Equation (45) is given by:

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

$2\Lambda =\frac{F_1}{{\omega }_L^2-{\omega }^2}$ (47)

The derivative of (46) is given as follows:

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

If $\epsilon \ne 0$ , as in primary resonance, the derivative of (47), using the variation of constants, is given by:

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

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

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

The derivative of (47) is given by:

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

By substituting (46), (48), and (30) into (45), we obtain:

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

Using (50) and (52), we obtain:

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

$\begin{array}{c}\stackrel{}{}+\omega {\Lambda }^2\cos \left(\left(2{\omega }_L+\omega \right)t+2\beta \right)-\omega {\Lambda }^2\cos \left(\left(2{\omega }_L-\omega \right)t+2\beta \right)\}\\ +\frac{{\alpha }_5\epsilon }{{\omega }_L}\{-\left(a{\omega }_L{\Lambda }^2+\frac{a^3{\omega }_L}{8}\right)-a{\omega }_L{\Lambda }^2\left(\cos \left(2\omega t\right)+\cos \left(2{\omega }_{L}t+2\beta \right)\right)+\frac{a^3{\omega }_L}{8}\cos \left(4{\omega }_{L}t+4\beta \right)\\ +\left(\omega {\Lambda }^3-\frac{a^2{\omega }_L\Lambda }{2}+\frac{a^2\omega \Lambda }{4}\right)\cos \left(\left({\omega }_L+\omega \right)t+\beta \right)-\left(\omega {\Lambda }^3+\frac{a^2{\omega }_L\Lambda }{2}+\frac{a^2\omega \Lambda }{4}\right)\cos \left(\left({\omega }_L-\omega \right)t+\beta \right)\\ +\left(\frac{a^2{\omega }_L\Lambda }{2}+\frac{a^2\omega \Lambda }{4}\right)\cos \left(\left(3{\omega }_L+\omega \right)t+3\beta \right)+\left(\frac{a^2{\omega }_L\Lambda }{2}-\frac{a^2\omega \Lambda }{4}\right)\cos \left(\left(3{\omega }_L-\omega \right)t+3\beta \right)\\ +\left(a\omega {\Lambda }^2+\frac{a{\omega }_L{\Lambda }^2}{2}\right)\cos \left(\left(2{\omega }_L+2\omega \right)t+2\beta \right)+\left(-a\omega {\Lambda }^2+\frac{a{\omega }_L{\Lambda }^2}{2}\right)\cos \left(\left(2{\omega }_L-2\omega \right)t+2\beta \right)\\ \stackrel{}{}+\omega {\Lambda }^3\cos \left(\left({\omega }_L+3\omega \right)t+\beta \right)-\omega {\Lambda }^3\cos \left(\left({\omega }_L-3\omega \right)t+\beta \right)\}\end{array}$ (53)

By substituting (53) into (50), we obtain:

(54)

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

3.2.1. Superharmonic Resonance

In this case, we have $\omega \approx \frac{1}{3}{\omega }_L$ et $\omega \approx \frac{1}{2}{\omega }_L$ .

First case: $\omega \approx \frac{1}{3}{\omega }_L,\Lambda =\frac{9F_1}{{\omega }_L^2}$

The expressions for $\sin \left(\left({\omega }_L-3\omega \right)t+\beta \right)$ and $\cos \left(\left({\omega }_L-3\omega \right)t+\beta \right)$ vary slowly with time. Equations (53) and (54) then become:

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

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

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

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

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

From system (57), we obtain:

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

Second case: $\omega \approx \frac{1}{2}{\omega }_L,\Lambda =\frac{2F_1}{3{\omega }_L^2}$

The expressions for $\sin \left(\left({\omega }_L-2\omega \right)t+\beta \right)$ and $\cos \left(\left({\omega }_L-2\omega \right)t+\beta \right)$ vary slow-ly with time. Equations (53) and (54) then become:

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

By setting $\gamma =\left({\omega }_L-2\omega \right)t+\beta$ et $2\omega ={\omega }_L+\epsilon \sigma$ , the system (59) becomes:

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

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

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

From system (61), we obtain:

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

3.2.2. Subharmonic Resonance

First case: $\omega \approx 3{\omega }_L,\Lambda =\frac{-F_1}{16{\omega }_L^2}$

The expressions for $\sin \left(\left(3{\omega }_L-\omega \right)t+\beta \right)$ and $\cos \left(\left(3{\omega }_L-\omega \right)t+\beta \right)$ vary slow-ly with time. Equations (53) and (54) then become:

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

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

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

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

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

From system (65), we obtain:

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

Second case: $\omega \approx 2{\omega }_L,\Lambda =\frac{-F_1}{6{\omega }_L^2}$

The expressions for $\sin \left(\left(2{\omega }_L-\omega \right)t+2\beta \right)$ and $\cos \left(\left(2{\omega }_L-\omega \right)t+2\beta \right)$ vary slowly with time. Equations (53) and (54) then become:

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

By setting $\gamma =\left(2{\omega }_L-\omega \right)t+2\beta$ , $\omega =2{\omega }_L+\epsilon \sigma$ , the system (67) gives:

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

In a steady state, we obtain:

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

From system (69), we obtain:

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

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.

4. Results and Discussion

4.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.

Table 1. Geometric and material properties of the beam [12].

	Elastic layer	Piezoelectric layer
L: Length (m)	1	1
H: Width (m)	$H=5h$	$H=5h$
h:Total thickness (h = 0.01)	$h_c=5h÷6$	$h_a=h_s=H÷12$
Young’s modulus (Pa)	$E_c=6.9\times {10}^{10}$	---
Density (Kg·m−3)	${\rho }_c=2766$	${\rho }_s=7500$
$c_{11}^{*}$ (Pa)	---	6.98 × 1010
$e_{11}^{*}$ (C·m−2)	---	−23.2
${ϵ}_{11}^{*}$ (F·m−1)	---	1.73 × 10−8

Table 2. Coefficient values corresponding to the three boundary conditions [12].

	S-S	C-S	C-C
${\mu }_1$	1.4924 × 103	3.6420 × 103	7.6689 × 103
${\mu }_2$	−3.4548 × 105	−2.9565 × 105	0
${\mu }_3$	−1.7760 × 107	−2.3151 × 107	−2.1882 × 107
${\mu }_4$	−4.8367 × 105	−4.3922 × 105	0
${\mu }_5$	−3.5520 × 107	−4.6302 × 107	−4.3765 × 107
$C_1$	1.7333 × 104	4.2300 × 104	8.9070 × 104
$C_2$	1.5075 × 107	0.18569 × 107	0
$C_3$	4.9133 × 108	6.4047 × 108	6.0537 × 108

Using the feedback control parameters, the ${\mu }_i$ and $C_i$ coefficients in Table 2 are given by:

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

$\begin{array}{l}{C}_1=\frac{1}{M}\left(E_c{I}_c+2c_{11}^{*}\left(I_s+S_s{z}_s^2\right)+\frac{{\left(ES\right)}_{pe}}{S_s}\left(2I_s+S_s{z}_s^2\right)\right){\displaystyle {\int }_0^L{\phi }_{xxxx}\left(x\right)\phi \left(x\right)\text{d}x};\\ C_2=\frac{{\left(ES\right)}_{pe}{z}_s}{M}\left(\frac{1}{L}{\displaystyle {\int }_0^L{\phi }_{xx}\left(x\right)\phi \left(x\right)\text{d}x}{\displaystyle {\int }_0^L{\phi }_{xx}\left(x\right)\text{d}x}-{\displaystyle {\int }_0^L\left\{{\left({\phi }_{xx}\left(x\right)\right)}^2+{\phi }_x\left(x\right){\phi }_{xxx}\left(x\right)\right\}{\phi }_x\left(x\right)\text{d}x}\right);\end{array}$

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

4.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]-[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.

Figure 2. Comparison of the three analytical methods on the amplitude-frequency responses of the S-S beam.

Figure 3. Jumps in the increasing direction of frequencies.

Figure 4. Jumps in the decreasing direction of frequencies.

Figure 5. Amplitude-frequency responses for the three boundary conditions.

Figure 6. Free nonlinear amplitude-frequency responses for the three boundary conditions.

Figure 7. Amplitude-frequency responses of the three boundary conditions.

Figure 8. Amplitude-frequency responses of the three boundary conditions.

4.3. Secondary Resonance

4.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 the C-S beam in Figure 11, nonlinear behaviors occur around the structure’s natural mode. In Figure 12, the excitation effect is negligible on the C-C beam. These behaviors result from the presence of clamping in the structure.

Figure 9. Amplitude-frequency superharmonic responses for $\omega \approx 1/3{\omega }_L$ of beam S-S.

Figure 10. Amplitude-frequency superharmonic responses for $\omega \approx 1/2{\omega }_L$ of beam S-S.

Figure 11. Amplitude-frequency superharmonic responses for $\omega \approx 1/2{\omega }_L$ of beam C-S.

Figure 12. Super-harmonic amplitude-frequency response for $\omega \approx 1/2{\omega }_L$ of the C-C beam.

4.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 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.

Figure 13. Subharmonic amplitude-frequency responses for $\omega \approx 3{\omega }_L$ of the S-S beam.

Figure 14. Subharmonic amplitude-frequency responses for $\omega \approx 3{\omega }_L$ of S-S beam.

Figure 15. Subharmonic amplitude-frequency responses for $\omega \approx 3{\omega }_L$ of the C-S beam.

Figure 16. Subharmonic amplitude-frequency responses for $\omega \approx 2{\omega }_L$ of the S-S beam.

The control of vibrational behavior can significantly influence the performance of energy harvesting devices. Maintaining the excitation amplitude and frequency close to the resonance point of a piezoelectric energy harvesting device allows for an increase in the harvested power. In this context, the proportional-derivative control applied in this study could optimize energy harvesting by adjusting the effective stiffness and damping, thereby modulating the vibrational response to maximize energy conversion efficiency.

Figure 17. Subharmonic amplitude-frequency responses for $\omega \approx 2{\omega }_L$ of S-S beam.

Figure 18. Subharmonic amplitude-frequency responses for $\omega \approx 2{\omega }_L$ of the C-C beam.

5. Conclusion

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 has proven 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.

Author Contributions

Barthelemy Zra Mha: Conceptualization, Methodology, Software, Writing, original draft. Guy Edgar Ntamack: Supervision, Methodology, Validation, Writing.

Data Availability Statement

The data supporting the findings of this study were generated via numerical simulation and are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Trindade, M.A. (2000) Contrôle hybride actif-passif des vibrations de structures par des matériaux piézoélectriques et viscoélastiques: Poutres sandwich-multicouches intelligentes. Thèse de doctorat, Conservatoire national des arts et métiers-CNAM.
[2]	Boudaoud, H. (2007) Modélisation de l’amortissement actif-passif des structures sandwichs. Thèse de doctorat, Université de Paul Verlaine-Metz.
[3]	Ekeocha, R.J.O. (2018) Vibration in Systems. Journal of Mechanical Engineering Research, 10, 1-6.[CrossRef]
[4]	Nelson, P.A. and Elliott, S.J. (1991) Active Control of Sound. Academic Press.
[5]	Fuller, C.R. and Von Flotow, A.H. (2002) Active Control of Sound and Vibration. IEEE Control Systems Magazine, 15, 9-19.
[6]	Eriksson, L.J. (1996) Active Sound and Vibration Control: A Technology in Transition. Noise Control Engineering Journal, 44, 1-9.[CrossRef]
[7]	Fuller, C.C., Elliott, S. and Nelson, P.A. (1996) Active Control of Vibration. Academic Press.
[8]	Hansen, C.H., Snyder, S.D., Qiu, X., Brooks, L.A. and Morceau, D.J. (1997) Active Control of Noise and Vibration. Routledge.
[9]	Olson, H.F. (1956) Electronic Control of Noise, Vibration, and Reverberation. The Journal of the Acoustical Society of America, 28, 966-972.[CrossRef]
[10]	Rizet, N. (1999) Contrôle actif de vibration utilisant des matériaux piézoélectriques. Thèse de doctorat, Institut National des Sciences Appliquées de Lyon.
[11]	Yan, L. (2013) Contrôle de vibration large bande à l’aide d’éléments piézoélectriques utilisant une technique non-linéaire. Thèse de doctorat, Institut National des Sciences Appliquées de Lyon.
[12]	Rechdaoui, M.S. (2010) Stabilité et contrôle actif des vibrations non linéaires des poutres. Thèse de doctorat, Université Abdelmalek Essaadi.
[13]	Rechdaoui, M.S., Azrar, L., Belouettar, S. and Potier-Ferry, M. (2009) Active Vibration Control of Piezoelectric Sandwich Beams at Large Amplitudes. Mechanics of Advanced Materials and Structures, 16, 98-109.[CrossRef]
[14]	Rechdaoui, M.S. and Azrar, L. (2013) Stability and Nonlinear Dynamic Analyses of Beam with Piezoelectric Actuator and Sensor Based on Higher-Order Multiple Scales Methods. International Journal of Structural Stability and Dynamics, 13, Article ID: 1350042.
[15]	Belouettar, S., Azrar, L., Daya, E.M., Laptev, V. and Potier-Ferry, M. (2008) Active Control of Nonlinear Vibration of Sandwich Piezoelectric Beams: A Simplified Approach. Computers & Structures, 86, 386-397.[CrossRef]
[16]	Nayfeh, A.H. and Mook, D.T. (1979) Nonlinear Oscillations. Wiley.
[17]	Richard, H.R. (1987) Perturbation Methods, Bifurcation Theory, and Computer Algebra. Springer.
[18]	Nayfeh, A.H. (1981) Introduction to Perturbation Techniques. Wiley.
[19]	Nayfeh, A.H. (1973) Perturbation Methods. Wiley.
[20]	Verhulst, F. (1996) Nonlinear Differential Equations and Dynamical Systems. Springer.
[21]	Rechdaoui, M.S., Azrar, L., Daya, E.M. and Belouettar, S. (2009) Mathematical Modeling of Subharmonic Sandwich Beams. Revue de Mécanique Appliquée et Théorique, 2, 3-14.
[22]	Nazemnezhad, R. and Hosseini-Hashemi, S. (2017) Exact Solution for Large Amplitude Flexural Vibration of Nanobeams Using Nonlocal Euler-Bernoulli Theory. Journal of Theoretical and Applied Mechanics, 55, 649-658.[CrossRef]
[23]	Awrejcewicz, J., Saltykova, O.A., Chebotyrevskiy, Y.B. and Krysko, V.A. (2011) Nonlinear Vibrations of the Euler-Bernoulli Beam Subjected to Transversal Load and Impact Actions. Nonlinear Studies, 18, 329-364.
[24]	Voß, T. and Scherpen, J.M.A. (2014) Port-Hamiltonian Modeling of a Nonlinear Timoshenko Beam with Piezo Actuation. SIAM Journal on Control and Optimization, 52, 493-519.[CrossRef]
[25]	Piollet, E. (2014) Amortissement non-linéaire des structures sandwichs à matériaux d’âme en fibres enchevêtrées. Thèse de doctorat, Université de Toulouse.
[26]	Fessal, K. (2016) Formulations et modélisation des vibrations par éléments finis de type solide-coque: Application aux structures sandwichs viscoélastiques et piézoélectriques. Thèse de doctorat, Université de Lorraine.
[27]	Daya, E.M., Azrar, L. and Potier-Ferry, M. (2004) An Amplitude Equation for the Non-Linear Vibration of Viscoelastically Damped Sandwich Beams. Journal of Sound and Vibration, 271, 789-813.[CrossRef]
[28]	Azrar, L., Benamar, R. and White, R.G. (2002) A Semi-Analytical Approach to the Non-Linear Dynamic Response Problem of Beams at Large Vibration Amplitudes, Part II: Multimode Approach to Steady State Forced Periodic Response. Journal of Sound and Vibration, 255, 1-41.
[29]	Ghadimi, M. and Kaliji, H.D. (2013) Application of the Harmonic Balance Method on Nonlinear Equations. World Applied Sciences Journal, 22, 532-537.
[30]	Adoukatl, C., Ntamack, G.E. and Azrar, L. (2022) High Order Analysis of a Nonlinear Piezoelectric Energy Harvesting of a Piezo Patched Cantilever Beam under Parametric and Direct Excitations. Mechanics of Advanced Materials and Structures, 30, 4835-4861.[CrossRef]