Magneto-electro-elastic modelling and nonlinear vibration analysis of bi-directional functionally graded beams

In the paper, a novel magneto-electro-elastic model of bi-directional (2D) functionally graded materials (FGMs) beams is developed for investigating the nonlinear dynamics. It is shown that the asymmetric modes induced by the 2D FGMs may significantly affect the nonlinear dynamic responses, which is tremendously different from previous studies. Taking into account the geometric nonlinearity, the nonlinear equation of motion and associated boundary conditions for the beams are derived according to the Hamilton’s principle. The natural frequencies and numerical modes of the beams are calculated by the generalized differential quadrature method. The frequency responses of the nonlinear forced vibration are constructed based on the Galerkin technique incorporating with the incremental harmonic balance approach. The influences of the material distributions, length–thickness ratio, electric voltage, magnetic potential as well as boundary condition on the nonlinear resonant frequency and response amplitude are discussed in details. It is notable that increasing in the axial and thickness FG indexes, negative electric potential and positive magnetic potential can lead to decline the nonlinear resonance frequency and amplitude peak, which is usually applied to accurately design the multi-ferroic composite structures. Furthermore, the nonlinear characteristics of motion can be regulated by tuning/tailoring the 2D FG materials.


Introduction
Many engineering devices can be regarded as a certain nonlinear dynamic system, which is described by the relevant nonlinear differential equations of motion. Due to its complexity, the investigation on the nonlinear dynamics of the engineering structures is a challenging subject. The development of the theories of nonlinear dynamics, the nonlinear free vibration [1], nonlinear primary resonance [2], nonlinear parametric resonance [3], nonlinear internal resonance [4] and nonlinear isolation [5] has been widely studied by using some developed approaches such as the method of multiple scales [6], harmonic balance method [7], the Galerkin truncation method [8], the differential quadrature method [9], and the finite difference method [10]. The predictions and understanding become possible for complicated nonlinear phenomena in the dynamic system, such as jumping, saturation, bifurcations and chaotic behaviors [10,11] to be beneficial for the application of the engineering devices.
Magneto-electro-elastic (MEE) material is one type of smart composites constructed by the interactions of the mechanical, electric and magnetic fields, which is extensively applied in various engineering areas, including health monitoring, nano-sensors and actuators applications, medical ultrasonic imaging, vibration control and energy harvesting. The MEE structure has a series of striking superiorities such as high strength-to-weight ratios, strong capabilities of loading-carry, and converting energy due to the magnetoelectro-elastic energy in the multi-ferroic composites can be exchanged from one to another form [12]. Therefore, the ability of the desired coupling effect in multi-physical fields has been paid attention to by numerous scholars. It is not surprising that the mechanical and physical behaviors in the MEE composites like wave propagation [13], dislocation [14], contact [15], crack [16], buckling and postbuckling [17], and vibration control problems [18] have been widely investigated for providing a theoretical basis to design and manufacture the intelligent composites.
Functional graded material (FGM) is a nonhomogeneous composite with the smooth and continuous gradation procedures, which is realized by changing the volume fractions of the constituents [19]. Combination of the MEE materials and FGMs can remarkably increase the energy conversion efficiency and has plenty of incorporated benefits compared with each of them. Thus, a number of engineering devices fabricated by functional graded magneto-electro-elastic (FGMEE) materials such as beam [20], plate [21], shell [22] and other micro-and nano-structures [23,24] have been concerned by some scientists to study the performances. For example, Li et al. [25] adopted the Rayleigh-Ritz method to examine vibration and acoustic behaviors of FGMEE porous plates with the assumption of the modified power-law distribution and the Voigt model. With the aid of the first-order shear deformation theory in conjunction with the nonlocal strain gradient theory, Liu and Lyu [26] put forward the dynamic analysis on the FGMEE nanofilm integrated with graphene face sheets. However, most of the above studies are limited to onedirectional (1D) FGMEE nanostructures. More generally, in this work the theoretical framework of a bidirectional (2D) FGMEE is proposed to satisfy the requirement of flexibility and high reliability design in lot of intelligent materials.
Up to date, the concept of 2D FGMs [27] has recently been regarded as an improved method and applied for designing some advanced structures. The mechanical properties of 2D FGMs can be tuned/tailed along two directions for manufacturing enhanced composites in comparison with one-directional FGM [28]. Some mechanical and physical behaviors of composites made by 2D FGMs are also deeply analyzed. For instance, Shafiei et al. [29] concentrated on the porous and 2D FGMs micro-and nano-beams and investigated the frequency characteristics influenced by the FGMs indexes, nonlocal and small-scale parameters and the porous types based on the generalized differential quadrature method (GDQM). Tang et al. [30] presented a novel model of 2D FGMs beams undergoing the large deformation motion and used the Galerkin technique combined with the homotopy analysis method to investigate the nonlinear frequency and free vibration response. Consider the hygrothermal effect, Tang and Ding [31] employed the GDQM to analyze the nonlinear dynamics of a 2D FGMs beam incorporating with coupled longitudinal and transverse motions. Chen et al. [32] established the mathematical model of 2D FGMs microbeams based on the Hamilton's principle and utilized the GDQM and pseudo-arc length continuation technique to study the static bifurcation and postbuckling configuration. Esmaeilzadeh and Kadkhodayan [33] proposed a porous 2D FGMs plate enhanced by outside stiffeners model and employed the dynamic relaxation technique method combined with the Newmark method to conduct the dynamic analysis.
Some intelligent structures, such as the MEE composites, have been widely applied in some advanced devices. However, they are constantly suffered by multi-directional loads. Here we herein show that the concept of 2D FGMs can be used to reconstruct the MEE composites, which can resist the loads in two directions synchronously, therefore effectively improving the structural stability and mechanical behaviors. But there is no published literature on the nonlinear dynamic analysis of 2D FGMEE mechanical system like beam structures. The mechanic mechanisms existing in the coupling structure between the MEE material and 2D FGMs are unclear and an interesting issue to unveil. Moreover, the nonlinearity is unavoidable for some advanced structures with large deformation. To address those problems, we firstly established the nonlinear governing equations of 2D FGMEE beams, then used the Galerkin technique in conjunction with the IHB method to trace the resonance solutions, and then explored the novel dynamic phenomena, which may provide an insight reference to the potential application of MEE composites with the action of multidirectional loads.

Theory and formulation
A 2D FGMEE beam made of CoFe 2 O 4 and BaTiO 3 materials subjected to an electric potential C x; z; t ð Þ and magnetic potential c x; z; t ð Þ, with thickness h, width b B and length L is taken into account, as shown in Fig. 1. Due to the heterogeneous distribution of the composites, we define that z and z m , respectively, represent the locations of the physical and geometric neutral surface, as well as the symbol x denotes the neutral axial coordinate. The effective physical and mechanical properties including the 2D FGMEE beams' elastic modulus c 11 , piezoelectric constant e 31 , dielectric constantk 33 , piezomagnetic constantf 31 , magnetoelectric constantd 33 , magnetic constantk 33 and mass density q are regarded as continuously varying based on the arbitrary functions along both thickness and axial orientations. Thus, the mixed rule can be stated as follows where the superscripts c and b, respectively, refer to the properties of the CoFe 2 O 4 and BaTiO 3 materials, ce represents the distance between the physical and geometric centers, the calculation of which is obtained in the previous literature [31]. In addition, n denotes the thickness FG index, which estimates the volume fraction exponent along the beam thickness orientation. The axial strain e x;z ð Þ of 2D FGMEE beams in the framework of the Euler-Bernoulli beam model is given by [28] whereû andŵ, respectively, denote the transverse and axial displacements at the central cylindrical plane. According to the Maxwell's equation, the electric and magnetic fields along the thickness orientation of 2D FGMEE beams can be written as [20] According to the linear elastic model [20], we have where r xx , D z and B z , respectively, are the axial stress, electric displacement and magnetic induction. Then, the virtual strain energy is in which V and A, respectively, denote the volume and cross-sectional area of 2D FGMEE beams. The normal result force N x and bending moment M x , respectively, are The variation of the kinetic energy for 2D FGMEE beams is expressed as The inertia components m pm1 and m pm2 in Eq. (15) are, respectively, calculated by The boundary conditions of the magnetic and electric fields for the 2D FGMEE beams are assumed as [23] The external virtual work can be expressed as in which YBp YBm where N p and N m , respectively, are the electric and magnetic forces along the longitudinal orientation. Introducing Hamilton's principle in the following form Substituting Eqs. (13) and (15) into Eq. (21), and then applying the transformation of integration by parts, we obtain the governing equations and boundary conditions In order to study the coupling effect of the MEE materials on the mechanical behaviors of 2D FG beams, we derive the magneto-and electro-profiles along the beam thickness orientation based on the magnetic and electric equations in Eqs. (24) and (25). Therefore, substituting Eqs. (11) and (12) into Eqs. (24) and (25), we have According to the Gramer's rule, one obtains wherê h pm1 ¼k 33ê31 Àd 33f31 k 33k33 Àd 2

33
: Based on Eqs. (17) and (28), we can yield the gradients of the magnetic and electric potentials for the 2D FGMEE beams as follows Consider Eqs. (10)- (12) and the gradients of the magnetic and electric potentials (Eq. (30)), one obtains the following expressions about the normal result force and bending moment as In order to describe the large deformation effect induced by the geometrical nonlinearity, the essential boundary conditions satisfied by two supported at both ends for the 2D FGMEE beams are considered as [28] By only taking into account the transverse motion, one yields By conducting the integration of Eq. (37) twice with respect to x, one obtainŝ where C 1 and C 2 , respectively, denote the integral constants calculated by the essential boundary conditions. Substituting Eq. (38) into Eq. (36), we have where Substituting Eqs. (39) and (40) into Eq. (38), one obtainŝ Substituting Eqs. (41) and (43) into Eq. (23) in conjunction with some mathematical simplifications, one obtains Equation (44) is the nonlinear equation of motion for the 2D FGMEE beams.
In the paper, the material properties along the axial orientation according to an exponential function are taken as follows where b are the axial FG index of 2D FGMEE beams. Substituting Eq. (45) into Eq. (44), we have Based on Eqs. (26) and (36), the boundary conditions reduce tô Introducing the following dimensionless quantities gets to , we obtain the dimensionless form of nonlinear dynamic model for the 2D FGMEE beams as follows b Fig. 3 The first two order modes of 2D FGMEE beams against the initial magnetic potential for different types of boundary Substituting Eq. (49) into Eqs. (47) and (48), the boundary conditions of 2D FGMEE beams are expressed as in the following form w 0; for a clamped-clamped (C-C) beam, w 0; for a hinged-hinged (H-H) beam, and w 0; for a clamped-hinged (C-H) beam.

Numerical modes
In view of the dynamic system of 2D FGMEE beams which is the weak nonlinearity, the dynamic behavior predicted by the nonlinear modes is very close to that predicted by the corresponding linear modes. Therefore, we apply the numerical modes of corresponding linear system to study the nonlinear vibration characteristic of 2D FGMEE beams. By neglecting the Use the GDQM to solve Eq. (54); we obtain the natural frequencies and numerical modes. According to the GDQM combined with the Chebyshev-Gauss-Lobatto technique [34], one leads to the distribution of mesh points in the following form The mth derivative of the beam displacement w is where N denotes the number of mesh points; the symbol C m ð Þ kj m ¼ 1 ð Þcan be calculated as follows Applying the derived procedure of the weighting coefficient in the GDQM leads to the higher-order derivatives Consider Eq. (56) and refer the literature [34], we can derive the discrete formats associated with the terms about kinematic variables in Eq. (54) as follows where the over dot is the differentiation with respect to t. Substituting Eq. (60) into Eq. (54), one obtains a series of ordinary differential equations to the nonboundary mesh points Using the GDQM approximation to discretize the boundary conditions, we have for a C-C 2D FGMEE beam, for a H-H 2D FGMEE beam, and b Fig. 7 Influence of the axial FG index on frequency responses of 2D FGMEE beams for various boundary conditions (L=h ¼ 20; c ¼ 0:01; P ¼ 0:03; n ¼ 1) ð64Þ for a C-H 2D FGMEE beam. By assembling the matrices associated with boundary conditions and Eq. (61), one leads to where d and j, respectively, are the domain and the boundary points. The matrixes corresponding to nonboundary elements are written as and the associated terms in the boundary points are given as b Fig. 8 Influence of the thickness FG index on frequency responses of 2D FGMEE beams for various boundary conditions (L=h ¼ 20; c ¼ 0:01; P ¼ 0: Eq. (65) are considered as b Fig. 9 Influence of the length-thickness ratio on frequency responses of 2D FGMEE beams for various boundary conditions ( c ¼ 0:01; P ¼ 0:03; n ¼ 1; Ã T and x L , respectively, are the numerical modes and natural frequency of 2D FGMEE beams. Substituting Eq. (71) into Eq. (65), and then eliminating the boundary points yields to a standard eigenvalue problem in the following form where jj K jd . By solving Eq. (72), we can obtain the frequencies and numerical modes of 2D FGMEE beams. During the simulation procedure, we consider the beams comprised of CoFe 2 O 4 and BaTiO 3 composites with the mechanical properties, as listed in Table 1. In addition, the parameters of the beam geometry are taken as: beam thickness h ¼ 2 m, width b B ¼ 1 m and L ¼ 40 m. Figures 2 and 3, respectively, depict the effects of the initial electric potential and initial magnetic potential on the first two modes of the 2D MEE beams for different types of boundary conditions. It is observed that the 2D FG structures apparently provide the asymmetric modes, which is different from the classical case. It is also revealed that the initial electric potential effect is inconspicuous. In addition, the negative values of magnetic potential have significantly impact on the first two modes, especially for the 2D FGMEE beams with the H-H and C-C end conditions. It is concluded that the MEE materials essentially affect the asymmetric modes of 2D FG structures, which can further change the nonlinear mechanical behaviors.

Nonlinear vibration analysis of the 2D FGMEE beams
In the section, the nonlinear vibration of the 2D FGMEE beams is studied. According to the secondorder Galerkin technique, the transverse response is defined as b Fig. 10 Influence of the damping coefficient on frequency responses of 2D FGMEE beams for various boundary conditions (L=h ¼ 20; Þis the ith numerical mode of corresponding degenerate system obtained from Sect. 3 and X ¼ X 1 X 2 ½ T represent the truncation mode matric. Considering the transverse forcing excitation of the beam as P cos x t ð Þ (x denotes the excitation frequency and P is the forcing amplitude), and Substituting Eq. (73) into Eq. (50), then conducting the Galerkin procedure, we havẽ whereM,K,K 2 ð Þ andK 3 ð Þ , respectively, represent the matrices of the mass, linear stiffness, quadratic and cubic nonlinear stiffness. These matrices can be written as Fig. 11 Influence of the forcing amplitude on frequency responses of 2D FGMEE beams for various boundary conditions (L=h ¼ 20; c ¼ 0:01; n ¼ 1; Besides,F is the vector of generalized force, which is written asF ¼ R 1 0 XPdn. In the study, we assume the viscous damping in the numerical simulation. Therefore, the proportional damping matrixC can be defined as [36,37] C ¼ 2c where c represents the damping coefficient. In order to accurately predict the nonlinear responses of 2D FGMEE beams, we use the incremental harmonic balance (IHB) method [38] For the IHB technique, the first step is performing the incremental procedure. Considering the excitation and response parameters at a state as F j0 ; -0 ; f j0 , one leads to the neighboring state Inserting Eq. (81) into (80) and neglecting higherorder terms, we have the incremental equation in the following form Fig. 12 Influence of the initial electric potential on frequency responses of 2D FGMEE beams for various boundary conditions (L=h ¼ 20; c ¼ 0:01; P ¼ 0:03; n ¼ 1; Besides,Q 0 ¼ f 10 f 20 ½ T , DQ ¼ Df 1 Df 2 ½ T and RE represents the residual vector, which is considered as tending to zero when the responses are gradually closed to the accurate value. The second step is taken as conducting the process of Galerkin averaging. Thus, the response and its incremental variables are written as follows Then,Q 0 and DQ can be rewritten as the following forms Substituting Eq. (87) into Eq. (82) and applying the harmonic balance method in one cycle to conduct the Galerkin procedure, one gets b Fig. 13 Influence of the initial magnetic potential on frequency responses of 2D FGMEE beams for various boundary conditions (L=h ¼ 20; c ¼ 0:01; P ¼ 0:03; n ¼ 1; By using the integration to Eq. (88), we merge into a series of algebraic equations with DB _ , Dand DF as the unknown variables, where It is noticed that Dand DF existed in Eq. (89) make the number of equations smaller than that of incremental variables. In the view of this, we only focus on the solutions of 2D FGMEE beams at a given excitation in this study. Consequently, Eq. (89) reduces to The solving procedure starts from an initial conjecture solution, and then the incremental and iterative procedures [38] are applied to construct the frequency response curves according to point by point. The iterative process stops until the 2-norm of R is smaller than the error (1:0 Â 10 À5 ). Finally, one can obtain the frequency responses of 2D FGMEE structures with desired accuracy.
Based on Floquet's theory [39], the stability of a steady-state responseQ 0 obtained by the IHB technique is determined by the following form of linearized perturbation equation For convenience, Eq. (81) can be expressed in the following state equation It is mentioned that the steady-state responseQ 0 has the same period (T ¼ 2p) as the state matrixÃs ð Þ. Therefore, the stability ofQ 0 is dependent on the eigenvalues of the matrixÃ generated by transforming Ys 0 ð Þ toỸs 0 þ T ð Þ. The stability criteria are thatQ 0 is unstable if the moduli of all eigenvalues of the transition matrixÃ are bigger than or equal to 1, otherwiseQ 0 is stable.
In order to illustrate the accuracy of the analytical solutions based on the IHB method, we compare the frequency response curves of uniform and 2D FGMEE beams calculated by the present method with those obtained by Runge-Kutta method, as shown in Fig. 4. It is seen that the numerical method can only provide the stable solutions, but the unstable solutions are also obtained by the analytical method. Besides, the stable solutions of IHB technique exhibit very good agreement with those of Runge-Kutta method, which validates the effectiveness of the IHB solutions. It is also noticed that there are one apparent superharmonic response existing in Fig. 2b and some typical nonlinear behaviors including multi-values, bifurcation and jump phenomena observed in the MEE system. Compared to the uniform MEE beams, the 2D FG structures exhibit more bending behaviors in the resonance region. This is because the presence of the bend-extension coupling influence in the FGMs beams produces a quadratic nonlinear term. To further interpret this bending phenomenon, we give the gap between the frequency responses of beams predicted by the original model and those obtained in the neglecting quadratic nonlinearity case, as shown in Fig. 5. It is indicated that the quadratic nonlinearity caused by the FGMs cannot be neglected to predict the mechanical behaviors of 2D FGMEE beams, especially for the responses in the resonance region.
To further demonstrate the essentiality of frequency responses obtained by the numerical modes in the 2D FGMEE beams, Figure 6 displays a comparison between the frequency responses obtained using the numerical modes and the analytical modes (sin kpn ð Þ; k ¼ 1; 2) for different values of axial and thickness FG indexes. It is indicated that the response predicted by the numerical modes approaches the classical case for the uniform MEE beams, but the difference between the two responses enlarges as increasing of the values of b and n. It can be concluded that the analytical modes do not accurately predict the solutions of the 2D FGMEE beams with higher values  This is because the growth of b and n will increase the contribution of the BaTiO 3 material on the 2D FGMEE structure which reduces its stiffness. It is also mentioned that with the ends become stiffer, the linear and resonance frequencies raise, the response peak descends and the hardening-type behavior of 2D FGMEE beams is weakened.
The variations in response amplitude of 2D FGMEE beams with changes as different values of length-thickness ratio for various end supports are shown in Fig. 9. For all types of end support, it can be seen that the increasing in the amplitude peak, the decreasing in the linear and resonance frequencies are obtained as well as the typical hardening nonlinearity reduces when an increasing in L/h. Furthermore, the influence of L/h in the case of 2D FGMEE beams with the H-H condition is more remarkable than that of other cases.   11, respectively, display the response amplitude of 2D FGMEE beams versus the frequency ratio associated with different values of the damping coefficient and forcing amplitude for different types of edge supports. Regardless of all types of end condition, it is revealed that the ascent of damping coefficient leads to the lessening of amplitude peak, while the amplitude peak increases with the growth of forcing amplitude. In addition, the decreasing of damping coefficient or the increasing of forcing amplitude will broaden the width of the frequency domain of multi-valued solutions.
The effects of initial electric potential on the frequency response curves of 2D FGMEE beams with H-H, C-H and C-C boundary conditions are presented, as shown in Fig. 12. It is demonstrated that initial electric potential significantly affects the response curves, especially for H-H beams. For the positive values of electric potential, the resonance frequency, response peak and hardening-type behavior are altogether intensified by increasing of initial electric potential. However, the opposite trends are obtained for the negative values of electric potential. The reason is that imposing positive/negative electric potential means to produce the axial compressive/ tensile forces in the 2D FGMEE beams. Figure 13 explores the influences of initial magnetic potential on the frequency response curves corresponding to various end supports. It is found that the positive magnetic potential reduces the resonance frequency, response peak and hardening-type behavior, whereas the negative magnetic potential has an  enhancing influence. Furthermore, the influence of initial magnetic potential is more prominent for the H-H edge support compared with that of C-H and C-C boundary conditions. In order to further evaluate the 2D FG indexes effects on the nonlinear dynamics of the MEE beams, the midpoint of the beam is set as the example of calculation; we use the numerical methods of nonlinear dynamics such as bifurcation, phase-plane, time history, and Poincaré section diagrams to explore the complicated dynamic phenomena. The excitation amplitude is chosen as the bifurcation parameter; the 2D FG indexes are set to (n ¼ 0; b ¼ 0), (n ¼ 0; b ¼ 1), (n ¼ 1; b ¼ 0) and (n ¼ 1; b ¼ 1). which, respectively, represents the uniform, axially FG, thickness FG, and 2D FG cases, as shown in Figs. 14, 15, 16, and 17. It is seen that these systems all occur the jump and bifurcation phenomena. With the increase in the excitation amplitude, the uniform beam system displays multiple forms of motion such as period-1, period-2, quasiperiodic oscillations, those complicated characteristics of motion are, respectively, illustrated in Figs. 18, 19 and 20. As the excitation amplitude is increased, the axially FG beam shows period-1, period-2 motions, Figs. 21, 22 and 23 confirm the motion behaviors. To the thickness FG and 2D FG cases, the two systems exhibit only the period-1 oscillation. The type motion characteristics are, respectively, shown in Figs   A novel model of 2D FG Euler-Bernoulli beams subjected to magneto-electro-elastic effects based on the geometrical nonlinearity is proposed to study the nonlinear dynamics. According to the proposed model, the nonlinear equation of motion is established using a variational formulation in Hamilton's principle. The GDQM is employed to predict the natural frequencies and numerical modes. We find that the MEE materials significantly affect the asymmetric modes, which can further vary the mechanics of 2D FG beams. Then, the IHB method incorporating with the Galerkin technique is applied to trace the solution trajectory for different types of end support. The stability of solution is discerned using the Floquet's theory. A comparison between the frequency responses of 2D FG MEE beams obtained by the IHB and Runge-Kutta methods confirms the accuracy of present solutions. Numerical examples exhibit the effects of different parameters, such as the electric potential and magnetic potential, on the nonlinear behaviors of the 2D FGMEE beams. Bifurcation, phase-plane, time history, and Poincaré section diagrams further examine the complicated nonlinear phenomena of the system. The main findings are listed as follows: (1) The nonlinear resonance response of 2D FGMEE beams influenced by the asymmetric modes in the 2D FG materials is tremendously different from that of the uniform and 1D case.