On the Bending and Vibration Analysis of Functionally Graded Magneto-Electro-Elastic Timoshenko Microbeams

Abstract

In this paper, a new magneto-electro-elastic functionally graded Timoshenko microbeam model is developed by using the variational formulation. The new model incorporates the extended modified couple stress theory in order to describe the microstructure effect. The power-law variation through the thickness direction of the two-phase microbeams is considered. By the direct application of the derived general formulation, the static bending and free vibration behavior of the newly developed functionally graded material microbeams are analytically determined. Parametric studies qualitatively demonstrate the microstructural effect as well as the magneto-electro-elastic multi-field coupling effect. The proposed model and its classic counterpart produce significant differences for thin graded magneto-electro-elastic Timoshenko microbeams. The thinner the microbeam is, the larger the difference becomes.

1. Introduction

Currently, magneto-electro-elastic (MEE) materials have attracted more and more attention. MEE materials can realize the mutual conversion between magnetic, electrical, and mechanical energies. Such characteristics have found important applications in stability controlling, actuating, health monitoring, medical ultrasonic, and some smart structure technologies [1,2,3]. In addition, functionally graded materials (FGMs) are characterized by continuous changes in material properties [4,5,6]. The mechanical properties of MEE materials synthesized from functionally graded materials are of great significance in both research and industrial fields [7,8]. In recent years, the research on investigating magneto-electro-elastic functionally graded materials (MEE-FGMs) on thin beams and plates has become a major trend. Bhangale and Ganesan [9] studied the free vibration behavior of anisotropic and linear MEE-FGM plates. Sladek et al. [10] proposed a meshless method for the bend analysis of circular MEE-FGM plates. Vinyas et al. [11] studied the effectiveness of utilizing MEE-FGM plates in precise frequency responses control. Mahesh and Harursampath [12] and Mahesh [13] evaluated nonlinear deflections of MEE-FGM porous flat panels and shells subjected to mechanical, electrical, and magnetic loads, respectively. However, numerous experiments [14,15] have proved that thin beams and plates usually exhibit size effects, (i.e., the thinner, the stiffer). Such size effects arise from non-local interactions of material particles at a very small scale, which cannot be described by classical theories at the micron or nanometer level due to a lack of any material length scale parameters. Therefore, it is necessary to develop thin MEE-FGM structure models based on non-classical theories.

In order to predict the size effects, numerous theories have been proposed with additional material parameters, such as non-local theories [16], couple stress theories [17,18,19], strain gradient theories [20,21,22], and a series of simpler versions [23,24,25,26,27,28]. These theories were successfully applied to develop size-dependent structure models for very small scales. For example, based on nonlocal theories, a number of MEE/MEE-FGM beam and plate models have been developed to capture non-local size effects [29,30,31,32], in which a non-local medium, including long-range material interactions, is adopted. Lim et al. [33] proposed a non-local strain gradient theory to include both non-local and strain gradient effects, and the bending, buckling, and free variation problems of FGM beams have been solved [34,35]. In addition, the modified couple stress theory (MCST) [24,25] contains only one additional parameter for isotropic materials. This MCST and its extended versions only consider the symmetrical part of the curvature tensor, which leads to fewer material parameters than their classical counterparts. In view of the great difficulties for determining additional parameters and interpreting the relevant microstructures, these modified theories have been applied to build micro/nano-beam and periodic composite pipe models [36,37,38,39,40,41], from which a microstructure-dependent stiffness is revealed. Recently, three such models have been proposed for MEE Timoshenko homogeneous beams [39] and MEE homogeneous plates [42,43] based on the extended modified couple stress theory. However, to the best of our knowledge, the extended modified couple stress theory is not applicable to MEE-FGM microbeams, which are inhomogeneous and might be helpful for smart devices miniaturization [44,45,46,47,48]. This motivated the present work.

The present work uses the extended modified couple stress theory to develop a MEE-FGM Timoshenko microbeam model for the first time and analytically solves the static bending and free vibration problems of the new model.

2. Materials and Methods

Consider a two-phase FGM microbeam with length L, width b and thickness h under the combined electric, magnetic, and mechanical loadings, as shown in Figure 1. The effective material properties P(z) (i.e., elastic stiffness, couple stress stiffness, piezoelectric constant, piezomagnetic constant, dielectric constant, magnetic permeability constant, magneto-dielectric constant and density) of the current microbeam change continuously in the thickness direction based on a power-law distribution [36], where P₁ and P₂ are the material properties of material I and II, respectively. The functionally graded power-law index n determines the material distribution across the thickness.

$$P(z)=(P_1-P_2){\left(\frac{z}{h}+\frac{1}{2}\right)}^n+P_2,$$

(1)

Based on the extended modified couple stress theory [42,49], the constitutive equations for transversely isotropic magneto-electro-elastic materials are given by [39,42,50,51].

$$\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}\\ {\sigma }_{yz}\\ {\sigma }_{zx}\\ {\sigma }_{xy}\end{array}\right\}=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{12}& C_{11}& C_{13}& 0& 0& 0\\ C_{13}& C_{13}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& \frac{C_{11}-C_{12}}{2}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ 2{\epsilon }_{yz}\\ 2{\epsilon }_{zx}\\ 2{\epsilon }_{xy}\end{array}\right\}-\left[\begin{array}{ccc}0& 0& q_{31}\\ 0& 0& q_{31}\\ 0& 0& q_{33}\\ 0& q_{15}& 0\\ q_{15}& 0& 0\\ 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{H}_x\\ H_y\\ H_z\end{array}\right\}-\left[\begin{array}{ccc}0& 0& e_{31}\\ 0& 0& e_{31}\\ 0& 0& e_{33}\\ 0& e_{15}& 0\\ e_{15}& 0& 0\\ 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{E}_x\\ E_y\\ E_z\end{array}\right\},$$

(2)

$$\left\{\begin{array}{c}{m}_{xx}\\ m_{yy}\\ m_{zz}\\ m_{yz}\\ m_{zx}\\ m_{xy}\end{array}\right\}=\left[\begin{array}{cccccc}{A}_{11}& A_{12}& A_{13}& 0& 0& 0\\ A_{12}& A_{11}& A_{13}& 0& 0& 0\\ A_{13}& A_{13}& A_{33}& 0& 0& 0\\ 0& 0& 0& A_{44}& 0& 0\\ 0& 0& 0& 0& A_{44}& 0\\ 0& 0& 0& 0& 0& \frac{A_{11}-A_{12}}{2}\end{array}\right]\left\{\begin{array}{c}{\chi }_{xx}\\ {\chi }_{yy}\\ {\chi }_{zz}\\ 2{\chi }_{yz}\\ 2{\chi }_{zx}\\ 2{\chi }_{xy}\end{array}\right\},$$

(3)

$$\begin{array}{ll}\left\{\begin{array}{c}{D}_x\\ D_y\\ D_z\end{array}\right\}& =\left[\begin{array}{cccccc}0& 0& 0& 0& e_{15}& 0\\ 0& 0& 0& e_{15}& 0& 0\\ e_{31}& e_{31}& e_{33}& 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ 2{\epsilon }_{yz}\\ 2{\epsilon }_{zx}\\ 2{\epsilon }_{xy}\end{array}\right\}+\left[\begin{array}{ccc}{s}_{11}& 0& 0\\ 0& s_{11}& 0\\ 0& 0& s_{33}\end{array}\right]\left\{\begin{array}{c}{E}_x\\ E_y\\ E_z\end{array}\right\}\\ & +\left[\begin{array}{ccc}{d}_{11}& 0& 0\\ 0& d_{11}& 0\\ 0& 0& d_{33}\end{array}\right]\left\{\begin{array}{c}{H}_x\\ H_y\\ H_z\end{array}\right\}\end{array},$$

(4)

$$\begin{array}{ll}\left\{\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right\}& =\left[\begin{array}{cccccc}0& 0& 0& 0& q_{15}& 0\\ 0& 0& 0& q_{15}& 0& 0\\ q_{31}& q_{31}& q_{33}& 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ 2{\epsilon }_{yz}\\ 2{\epsilon }_{zx}\\ 2{\epsilon }_{xy}\end{array}\right\}+\left[\begin{array}{ccc}{\mu }_{11}& 0& 0\\ 0& {\mu }_{11}& 0\\ 0& 0& {\mu }_{33}\end{array}\right]\left\{\begin{array}{c}{H}_x\\ H_y\\ H_z\end{array}\right\}\\ & +\left[\begin{array}{ccc}{d}_{11}& 0& 0\\ 0& d_{11}& 0\\ 0& 0& d_{33}\end{array}\right]\left\{\begin{array}{c}{E}_x\\ E_y\\ E_z\end{array}\right\}\end{array}.$$

(5)

where σ_ij, m_ij, D_i, B_i are the Cauchy stress tensor, the deviatoric part of the couple stress tensor, the electric displacements, and the magnetic fluxes, respectively. C_αβ (α, β = 1, 2, …, 6) is the elastic stiffness tensor, A_αβ (α, β = 1, 2, …, 6) is the couple stress stiffness tensor, e_iα and q_iα are the piezoelectric and piezomagnetic tensors, s_ij and μ_ij are the dielectric and magnetic permeability tensors, d_ij is the magneto-dielectric tensor, and ε_ij and χ_ij are, respectively, the infinitesimal strain and the symmetric curvature tensors, which are defined by

$${\epsilon }_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right),$$

(6)

$${\chi }_{ij}=\frac{1}{4}\left({\epsilon }_{ipq}{u}_{q,pj}+{\epsilon }_{jpq}{u}_{q,pi}\right),$$

(7)

with u_i being the displacement, and ε_ijk is the Levi-Civita symbol. In addition, E_k and H_k are, respectively, the electric field intensity and magnetic field intensity read

$$E_k=-{\Phi }_{,k} ,H_k=-M_{,k},$$

(8)

where Φ and M are the electric and magnetic potentials.

For a MEE Timoshenko beam with a uniform cross-section shown in Figure 1, the displacement field and electric and magnetic potentials can be given by [52,53,54,55]

$$u_1=u\left(x,t\right)-z\phi \left(x,t\right),u_2=0,u_3=w\left(x,t\right),$$

(9)

$$\Phi =-\cos \left(\frac{\pi }{h}z\right)\gamma \left(x,t\right)+\frac{2z}{h}{\gamma }_0, M=-\cos \left(\frac{\pi }{h}z\right)\zeta \left(x,t\right)+\frac{2z}{h}{\zeta }_0,$$

(10)

where u and w are the beam extension and deflection, φ represents the rotation angle, γ and ζ are the spatial variations of the electric and magnetic potentials along the x-direction, respectively. γ₀ and ζ₀ are, respectively, the external electric and magnetic potentials.

By substituting Equations (9) and (10) into Equations (6)–(8) yields

$${\epsilon }_{xx}=\frac{\partial u}{\partial x}-z\frac{\partial \phi }{\partial x},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\epsilon }_{xz}=\frac{1}{2}\left(\frac{\partial w}{\partial x}-\phi \right), \mathrm{others}=0,$$

(11)

$${\chi }_{xy}=-\frac{1}{4}\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \phi }{\partial x}\right), \mathrm{others}=0,$$

(12)

$$E_x=\cos \left(\frac{\pi }{h}z\right)\frac{\partial \gamma }{\partial x}, E_z=-\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right)\gamma -\frac{2}{h}{\gamma }_0, E_y=0,$$

(13)

$$H_x=\cos \left(\frac{\pi }{h}z\right)\frac{\partial \zeta }{\partial x}, H_z=-\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right)\zeta -\frac{2}{h}{\zeta }_0, H_y=0.$$

(14)

Based on Equations (11)–(14), the constitutive equations in Equations (2)–(5) can be obtained as

$${\sigma }_{xx}=C_{11}{\epsilon }_{xx}-e_{31}{E}_z-q_{31}{H}_z, {\sigma }_{xz}=2C_{44}{\epsilon }_{xz}-e_{15}{E}_x-q_{15}{H}_x,$$

(15)

$$m_{xy}=\left(A_{11}-A_{12}\right){\chi }_{xy},$$

(16)

$$D_x=2e_{15}{\epsilon }_{xz}+s_{11}{E}_x+d_{11}{H}_x, D_z=e_{31}{\epsilon }_{xx}+s_{33}{E}_z+d_{33}{H}_z,$$

(17)

$$B_x=2q_{15}{\epsilon }_{xz}+{\mu }_{11}{H}_x+d_{11}{E}_x, B_z=q_{31}{\epsilon }_{xx}+{\mu }_{33}{H}_z+d_{33}{E}_z.$$

(18)

From Equations (11)–(18), the first variation of the total strain energy in the current beam satisfying the extended modified couple stress theory over the time span [0, T] takes the form [39,42]

$$\begin{array}{ll}\delta {\displaystyle {\int }_0^{T}Udt}=& {\displaystyle {\int }_0^T{\displaystyle {\int }_0^L{\displaystyle {\int }_A({\sigma }_{xx}\delta {\epsilon }_{xx}+2{\sigma }_{xz}\delta {\epsilon }_{xz}+2m_{xy}\delta {\chi }_{xy}}}}\\ & -D_x\delta E_x-D_z\delta E_z-B_x\delta H_x-B_z\delta H_z)dAdxdt\end{array},$$

(19)

where A is the cross-sectional area.

The first variation of the kinetic energy of the Timoshenko beam over the time interval [0, T] is given by [53]

$$\delta {\displaystyle {\int }_0^{T}K}dt={\displaystyle {\int }_0^T{\displaystyle {\int }_0^L{\displaystyle {\int }_A\rho }}}\left(\frac{\partial u_1}{\partial t}\frac{\partial \delta u_1}{\partial t}+\frac{\partial u_3}{\partial t}\frac{\partial \delta u_3}{\partial t}\right)dAdxdt,$$

(20)

where ρ is the mass density.

Furthermore, the virtual work performed by the applied forces acting on the current Timoshenko beam over the time span [0, T] can be written as [55,56]

$$\delta {\displaystyle {\int }_0^{T}Wdt}={\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left[f\delta u+q\delta w\right]dxdt}},$$

(21)

where f and q are, respectively, the x- and z-components of the body force per unit length along the x-axis.

According to Hamilton’s principle [53,56],

$$\delta {\displaystyle {\int }_0^T\left[K-(U-W)\right]}dt=0.$$

(22)

Substituting Equations (19)–(21) into Equation (22), applying the fundamental lemma of the calculus of variations [57], and considering the arbitrariness of δu, δw, and δφ yield

$$\frac{\partial N_{xx}}{\partial x}+f=m_0\frac{{\partial }^{2}u}{\partial t^2}-m_1\frac{{\partial }^2\phi }{\partial t^2},$$

(23)

$$-\frac{\partial M_{xx}}{\partial x}+N_{xz}-\frac{1}{2}\frac{\partial Y_{xy}}{\partial x}=m_2\frac{{\partial }^2\phi }{\partial t^2}-m_1\frac{{\partial }^{2}u}{\partial t^2},$$

(24)

$$\frac{\partial N_{xz}}{\partial x}+\frac{1}{2}\frac{{\partial }^2{Y}_{xy}}{\partial x^2}+q=m_0\frac{{\partial }^{2}w}{\partial t^2},$$

(25)

$$\frac{\partial {\mathrm{\Lambda }}_x}{\partial x}+{\mathrm{\Lambda }}_z=0,$$

(26)

$$\frac{\partial {\mathrm{\Sigma }}_x}{\partial x}+{\mathrm{\Sigma }}_z=0.$$

(27)

as the equation of motion, and

$$N_{xx}=0 \mathrm{or} u=\overline{u} \mathrm{at} x=0 \mathrm{and} x=L,$$

(28)

$$M_{xx}+\frac{1}{2}{Y}_{xy}=0 \mathrm{or} \phi =\overline{\phi } \mathrm{at} x=0 \mathrm{and} x=L,$$

(29)

$$-N_{xz}-\frac{1}{2}\frac{\partial Y_{xy}}{\partial x}=0 \mathrm{or} w=\overline{w} \mathrm{at} x=0 \mathrm{and} x=L,$$

(30)

$$\frac{1}{2}{Y}_{xy}=0 \mathrm{or} \frac{\partial w}{\partial x}=\frac{\partial \overline{w}}{\partial x} \mathrm{at} x =0 \mathrm{and} x=L,$$

(31)

$${\mathrm{\Lambda }}_x=0 \mathrm{or} \gamma =\overline{\gamma } \mathrm{at} x=0 \mathrm{and} x=L,$$

(32)

$${\mathrm{\Sigma }}_x=0 \mathrm{or} \zeta =\overline{\zeta } \mathrm{at} x=0 \mathrm{and} x=L$$

(33)

as boundary conditions, where the overbar denotes the prescribed value. Note that the stress, electric, magnetic resultants, and mass inertias can be expressed as

$$N_{xx}={\displaystyle {\int }_A{\sigma }_{xx}dA}=A_{xx}\frac{\partial u}{\partial x}-B_{xx}\frac{\partial \phi }{\partial x}+A_{31}^e\gamma +A_{31}^q\zeta +N_x^E+N_x^H,$$

(34)

$$M_{xx}={\displaystyle {\int }_{A}z{\sigma }_{xx}dA}=B_{xx}\frac{\partial u}{\partial x}-D_{xx}\frac{\partial \phi }{\partial x}+B_{31}^e\gamma +B_{31}^q\zeta +M_x^E+M_x^H,$$

(35)

$$N_{xz}={\displaystyle {\int }_A{k}_s{\sigma }_{xz}dA}=k_s^2{A}_{xz}\left(\frac{\partial w}{\partial x}-\phi \right)-k_s{A}_{15}^e\frac{\partial \gamma }{\partial x}-k_s{A}_{15}^q\frac{\partial \zeta }{\partial x},$$

(36)

$$Y_{xy}={\displaystyle {\int }_A{m}_{xy}dA=F_{xy}\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{\partial \phi }{\partial x}\right)},$$

(37)

$${\mathrm{\Lambda }}_x={\displaystyle {\int }_A{D}_x\cos \left(\frac{\pi }{h}z\right)dA=k_s{A}_{15}^e\left(\frac{\partial w}{\partial x}-\phi \right)}+A_{11}^s\frac{\partial \gamma }{\partial x}+A_{11}^d\frac{\partial \zeta }{\partial x},$$

(38)

$${\mathrm{\Lambda }}_z={\displaystyle {\int }_A{D}_z\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right)dA=A_{31}^e\frac{\partial u}{\partial x}-B_{31}^e}\frac{\partial \phi }{\partial x}-A_{33}^s\gamma -A_{33}^d\zeta -N_{33}^{Es}-N_{33}^{Hd},$$

(39)

$${\mathrm{\Sigma }}_x={\displaystyle {\int }_A{B}_x\cos \left(\frac{\pi }{h}z\right)dA=k_s{A}_{15}^q\left(\frac{\partial w}{\partial x}-\phi \right)}+A_{11}^{\mu }\frac{\partial \zeta }{\partial x}+A_{11}^d\frac{\partial \gamma }{\partial x},$$

(40)

$${\mathrm{\Sigma }}_z={\displaystyle {\int }_A{B}_z\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right)dA=A_{31}^q\frac{\partial u}{\partial x}-B_{31}^q}\frac{\partial \phi }{\partial x}-A_{33}^{\mu }\zeta -A_{33}^d\gamma -N_{33}^{H\mu }-N_{33}^{Ed},$$

(41)

$$\left(m_0,m_1,m_2\right)={\displaystyle {\int }_A\rho \left(z\right)\left(1,z,z^2\right)}dA$$

(42)

where k_s denotes the shape correction factor [58,59], and

$$(A_{xx}, B_{xx},{ \mathrm{D}}_{xx})={\displaystyle {\int }_A{C}_{11}(z)(1, z, z^2)dA},$$

(43)

$$A_{xz}={\displaystyle {\int }_A{C}_{44}\left(z\right)dA},$$

(44)

$$F_{xy}=\frac{1}{4}{\displaystyle {\int }_A\left(A_{12}\left(z\right)-A_{11}\left(z\right)\right)}dA,$$

(45)

$$\left(A_{31}^e, B_{31}^e\right)={\displaystyle {\int }_A{e}_{31}}\left(z\right)\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right)\left(1, z\right)dA,$$

(46)

$$\left(A_{31}^q, B_{31}^q\right)={\displaystyle {\int }_A{q}_{31}}\left(z\right)\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right)\left(1, z\right)dA,$$

(47)

$$\left(N_x^E, M_x^E\right)={\displaystyle {\int }_A{e}_{31}\left(z\right)\frac{2{\gamma }_0}{h}\left(1, z\right)dA},$$

(48)

$$\left(N_x^H, M_x^H\right)={\displaystyle {\int }_A{q}_{31}\left(z\right)\frac{2{\zeta }_0}{h}\left(1, z\right)dA},$$

(49)

$$\left(A_{15}^e, A_{15}^q\right)={\displaystyle {\int }_A\left\{e_{15}\left(z\right)\cos \left(\frac{\pi }{h}z\right), q_{15}\left(z\right)\cos \left(\frac{\pi }{h}z\right)\right\}dA},$$

(50)

$$\left(A_{11}^s, A_{33}^s\right)={\displaystyle {\int }_A\left\{s_{11}\left(z\right){\cos }^2\left(\frac{\pi }{h}z\right), s_{33}\left(z\right){\left(\frac{\pi }{h}\right)}^2{\sin }^2\left(\frac{\pi }{h}z\right)\right\}dA},$$

(51)

$$\left(A_{11}^d, A_{33}^d\right)={\displaystyle {\int }_A\left\{d_{11}\left(z\right){\cos }^2\left(\frac{\pi }{h}z\right), d_{33}\left(z\right){\left(\frac{\pi }{h}\right)}^2{\sin }^2\left(\frac{\pi }{h}z\right)\right\}dA},$$

(52)

$$\left(N_{33}^{Es}, N_{33}^{Hd}\right)={\displaystyle {\int }_A\left\{s_{33}\left(z\right)\frac{2{\gamma }_0}{h}\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right), d_{33}\left(z\right)\frac{2{\zeta }_0}{h}\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right)\right\}dA},$$

(53)

$$\left(A_{11}^{\mu }, A_{33}^{\mu }\right)={\displaystyle {\int }_A\left\{{\mu }_{11}\left(z\right){\cos }^2\left(\frac{\pi }{h}z\right), {\mu }_{33}\left(z\right){\left(\frac{\pi }{h}\right)}^2{\sin }^2\left(\frac{\pi }{h}z\right)\right\}dA},$$

(54)

$$\left(N_{33}^{H\mu }, N_{33}^{Ed}\right)={\displaystyle {\int }_A\left\{{\mu }_{33}\left(z\right)\frac{2{\zeta }_0}{h}\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right), d_{33}\left(z\right)\frac{2{\gamma }_0}{h}\frac{\pi }{h}\sin \left(\frac{\pi }{h}z\right)\right\}dA}.$$

(55)

Based on Equations (23)–(55), it is found that the current MEE-FGM beam model can additionally capture the effects of couple stress, piezomagnetism, piezoelectricity, and MEE coupling, when compared to the classical FGM Timoshenko beam model.

3. Analytical Solution

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

In order to illustrate the newly developed model in Section 2, the static bending and free vibration problems of the current beam are solved in this section.

According to Equations (28)–(33), the relevant boundary conditions of a simply supported beam can be identified as

$$N_{xx}=0,$$

(56)

$${w|}_{x=0}={w|}_{x=L}=0,$$

(57)

$${M_{xx}|}_{x=0}={M_{xx}|}_{x=L}=0,$$

(58)

$${Y_{xy}|}_{x=0}={Y_{xy}|}_{x=L}=0,$$

(59)

$${\gamma |}_{x=0}={\gamma |}_{x=L}=0,$$

(60)

$${\zeta |}_{x=0}={\zeta |}_{x=L}=0.$$

(61)

It should be noted that the boundary of the electric and magnetic conditions given in Equations (60) and (61) are for an open circuit.

3.1. Static Bending

Consider Fourier solutions for u(x), φ(x), w(x), γ(x), and ζ(x):

$$u\left(x\right)={\displaystyle \sum _{k=1}^{\infty }{U}_k\cos \left(\frac{k\pi x}{L}\right)},$$

(62)

$$\phi \left(x\right)={\displaystyle \sum _{k=1}^{\infty }{\Phi }_k\cos \left(\frac{k\pi x}{L}\right)},$$

(63)

$$w\left(x\right)={\displaystyle \sum _{k=1}^{\infty }{W}_k\sin \left(\frac{k\pi x}{L}\right)},$$

(64)

$$\gamma \left(x\right)={\displaystyle \sum _{k=1}^{\infty }{\Gamma }_k\sin \left(\frac{k\pi x}{L}\right)},$$

(65)

$$\zeta \left(x\right)={\displaystyle \sum _{k=1}^{\infty }{Z}_k\sin \left(\frac{k\pi x}{L}\right)}$$

(66)

where U_k, Φ_k, W_k, Г_k, and Z_k are the Fourier coefficients to be determined. It can be shown that the Fourier solutions in Equations (62)–(66) satisfy the boundary conditions in Equations (56)–(61). In addition, the body force f is equal to zero, and the uniform load q(x) can also be expanded in Fourier series as:

$$q(x)={\displaystyle \sum _{k=1}^{\infty }{Q}_k\sin \frac{k\pi x}{L}},$$

(67)

where Q_k is a Fourier coefficient calculated by q(x) = p₀ in the current case as

$$Q_k=\frac{2p_0}{k\pi }\left[1-\cos \left(k\pi \right)\right]$$

(68)

According to the Equations (23)–(27), (62)–(66) and (67), the equilibrium equations of static bending problems can be written as

$$A_{xx}\frac{{\partial }^{2}u}{\partial x^2}-B_{xx}\frac{{\partial }^2\phi }{\partial x^2}+A_{31}^e\frac{\partial \gamma }{\partial x}+A_{31}^q\frac{\partial \zeta }{\partial x}=0,$$

(69)

$$\begin{array}{l}-B_{xx}\frac{{\partial }^{2}u}{\partial x^2}+D_{xx}\frac{{\partial }^2\phi }{\partial x^2}-B_{31}^e\frac{\partial \gamma }{\partial x}-B_{31}^q\frac{\partial \zeta }{\partial x}+k_s^2{A}_{xz}\left(\frac{\partial w}{\partial x}-\phi \right)\\ -k_s{A}_{15}^e\frac{\partial \gamma }{\partial x}-k_s{A}_{15}^q\frac{\partial \zeta }{\partial x}-\frac{1}{2}{F}_{xy}\left(\frac{{\partial }^{3}w}{\partial x^3}+\frac{{\partial }^2\phi }{\partial x^2}\right)=0\end{array},$$

(70)

$$k_s^2{A}_{xz}\left(\frac{{\partial }^{2}w}{\partial x^2}-\frac{\partial \phi }{\partial x}\right)-k_s{A}_{15}^e\frac{{\partial }^2\gamma }{\partial x^2}-k_s{A}_{15}^q\frac{{\partial }^2\zeta }{\partial x^2}+\frac{1}{2}{F}_{xy}\left(\frac{{\partial }^{4}w}{\partial x^4}+\frac{{\partial }^3\phi }{\partial x^3}\right)=-q,$$

(71)

$$k_s{A}_{15}^e\left(\frac{{\partial }^{2}w}{\partial x^2}-\frac{\partial \phi }{\partial x}\right)+A_{11}^s\frac{{\partial }^2\gamma }{\partial x^2}+A_{11}^d\frac{{\partial }^2\zeta }{\partial x^2}+A_{31}^e\frac{\partial u}{\partial x}-B_{31}^e\frac{\partial \phi }{\partial x}-A_{33}^s\gamma -A_{33}^d\zeta =0,$$

(72)

$$k_s{A}_{15}^q\left(\frac{{\partial }^{2}w}{\partial x^2}-\frac{\partial \phi }{\partial x}\right)+A_{11}^{\mu }\frac{{\partial }^2\zeta }{\partial x^2}+A_{11}^d\frac{{\partial }^2\gamma }{\partial x^2}+A_{31}^q\frac{\partial u}{\partial x}-B_{31}^q\frac{\partial \phi }{\partial x}-A_{33}^{\mu }\zeta -A_{33}^d\gamma =0.$$

(73)

Substituting Equations (62)–(66) into Equations (69)–(73) results in

$$\left[\begin{array}{ccccc}{S}_{11}& S_{12}& 0& S_{14}& S_{15}\\ S_{12}& S_{22}& S_{23}& S_{24}& S_{25}\\ 0& S_{23}& S_{33}& S_{34}& S_{35}\\ S_{14}& S_{24}& S_{34}& S_{44}& S_{45}\\ S_{15}& S_{25}& S_{35}& S_{45}& S_{55}\end{array}\right]\left[\begin{array}{c}{U}_k\\ {\Phi }_k\\ W_k\\ {\Gamma }_k\\ Z_k\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ -Q_k\\ 0\\ 0\end{array}\right],$$

(74)

where

$$\begin{array}{l}{S}_{11}=-A_{xx}{\left(\frac{k\pi }{L}\right)}^2, S_{12}=B_{xx}{\left(\frac{k\pi }{L}\right)}^2, S_{13}=0, S_{14}=A_{31}^e\left(\frac{k\pi }{L}\right), S_{15}=A_{31}^q\left(\frac{k\pi }{L}\right), \\ S_{22}=-D_{xx}{\left(\frac{k\pi }{L}\right)}^2-k_s^2{A}_{xz}+\frac{1}{2}{F}_{xy}{\left(\frac{k\pi }{L}\right)}^2, S_{23}=k_s^2{A}_{xz}\left(\frac{k\pi }{L}\right)+\frac{1}{2}{F}_{xy}{\left(\frac{k\pi }{L}\right)}^3, \\ S_{24}=-A_{31}^e\left(\frac{k\pi }{L}\right)-k_s{A}_{15}^e\left(\frac{k\pi }{L}\right), S_{25}=-A_{31}^q\left(\frac{k\pi }{L}\right)-k_s{A}_{15}^q\left(\frac{k\pi }{L}\right),\\ S_{33}=-k_s^2{A}_{xz}{\left(\frac{k\pi }{L}\right)}^2+\frac{1}{2}{F}_{xy}{\left(\frac{k\pi }{L}\right)}^4, S_{34}=k_s{A}_{15}^e{\left(\frac{k\pi }{L}\right)}^2, S_{35}=k_s{A}_{15}^q{\left(\frac{k\pi }{L}\right)}^2, \\ S_{44}=A_{11}^s{\left(\frac{k\pi }{L}\right)}^2+A_{33}^s, S_{45}=A_{11}^d{\left(\frac{k\pi }{L}\right)}^2+A_{33}^d, S_{55}=A_{11}^{\mu }{\left(\frac{k\pi }{L}\right)}^2+A_{33}^{\mu }.\end{array}$$

(75)

According to Equation (74), the Fourier coefficients U_k, Φ_k, W_k, Г_k, and Z_k will be solved. The solutions of u(x), φ(x), w(x), γ(x), and ζ(x) for the current simple supported beam can also be given by inserting these results into Equations (62)–(66).

3.2. Free Vibration

In the free vibration problem of the current beam, both the external forces are vanished (i.e., f = q = 0). Consider the following Fourier series expansions for u(x, t), φ(x, t), w(x, t), γ(x, t), and ζ(x, t):

$$u\left(x,t\right)={\displaystyle \sum _{k=1}^{\infty }{U}_k^V\cos \left(\frac{k\pi x}{L}\right)}{e}^{i{\omega }_{k}t},$$

(76)

$$\phi \left(x,t\right)={\displaystyle \sum _{k=1}^{\infty }{\Phi }_k^V\cos \left(\frac{k\pi x}{L}\right)}{e}^{i{\omega }_{k}t},$$

(77)

$$w\left(x,t\right)={\displaystyle \sum _{k=1}^{\infty }{W}_k^V\sin \left(\frac{k\pi x}{L}\right)}{e}^{i{\omega }_{k}t},$$

(78)

$$\gamma \left(x,t\right)={\displaystyle \sum _{k=1}^{\infty }{\Gamma }_k^V\sin \left(\frac{k\pi x}{L}\right)}{e}^{i{\omega }_{k}t},$$

(79)

$$\zeta \left(x,t\right)={\displaystyle \sum _{k=1}^{\infty }{Z}_k^V\sin \left(\frac{k\pi x}{L}\right)}{e}^{i{\omega }_{k}t}$$

(80)

where ω_k is the kth vibration frequency, $U_k^V$, $W_k^V$ ${\Phi }_k^V$, ${Г}_k^V$$,$ and $Z_k^V$ are Fourier coefficients. It should be noted that the Fourier series expansions in Equations (76)–(80) satisfy the boundary conditions in Equations (56)–(61). Based on Equations (76)–(80) and Equations (23)–(27), the equations of motion can be expressed as

$$A_{xx}\frac{{\partial }^{2}u}{\partial x^2}-B_{xx}\frac{{\partial }^2\phi }{\partial x^2}+A_{31}^e\frac{\partial \gamma }{\partial x}+A_{31}^q\frac{\partial \zeta }{\partial x}=m_0\frac{{\partial }^{2}u}{\partial t^2}-m_1\frac{{\partial }^2\phi }{\partial t^2},$$

(81)

$$\begin{array}{l}-B_{xx}\frac{{\partial }^{2}u}{\partial x^2}+D_{xx}\frac{{\partial }^2\phi }{\partial x^2}-B_{31}^e\frac{\partial \gamma }{\partial x}-B_{31}^q\frac{\partial \zeta }{\partial x}+k_s^2{A}_{xz}\left(\frac{\partial w}{\partial x}-\phi \right)\\ -k_s{A}_{15}^e\frac{\partial \gamma }{\partial x}-k_s{A}_{15}^q\frac{\partial \zeta }{\partial x}-\frac{1}{2}{F}_{xy}\left(\frac{{\partial }^{3}w}{\partial x^3}+\frac{{\partial }^2\phi }{\partial x^2}\right)=m_2\frac{{\partial }^2\phi }{\partial t^2}-m_1\frac{{\partial }^{2}u}{\partial t^2}\end{array},$$

(82)

$$k_s^2{A}_{xz}\left(\frac{{\partial }^{2}w}{\partial x^2}-\frac{\partial \phi }{\partial x}\right)-k_s{A}_{15}^e\frac{{\partial }^2\gamma }{\partial x^2}-k_s{A}_{15}^q\frac{{\partial }^2\zeta }{\partial x^2}+\frac{1}{2}{F}_{xy}\left(\frac{{\partial }^{4}w}{\partial x^4}+\frac{{\partial }^3\phi }{\partial x^3}\right)=m_0\frac{{\partial }^{2}w}{\partial t^2},$$

(83)

$$k_s{A}_{15}^e\left(\frac{{\partial }^{2}w}{\partial x^2}-\frac{\partial \phi }{\partial x}\right)+A_{11}^s\frac{{\partial }^2\gamma }{\partial x^2}+A_{11}^d\frac{{\partial }^2\zeta }{\partial x^2}+A_{31}^e\frac{\partial u}{\partial x}-B_{31}^e\frac{\partial \phi }{\partial x}-A_{33}^s\gamma -A_{33}^d\zeta =0,$$

(84)

$$k_s{A}_{15}^q\left(\frac{{\partial }^{2}w}{\partial x^2}-\frac{\partial \phi }{\partial x}\right)+A_{11}^{\mu }\frac{{\partial }^2\zeta }{\partial x^2}+A_{11}^d\frac{{\partial }^2\gamma }{\partial x^2}+A_{31}^q\frac{\partial u}{\partial x}-B_{31}^q\frac{\partial \phi }{\partial x}-A_{33}^{\mu }\zeta -A_{33}^d\gamma =0.$$

(85)

Using Equations (76)–(80) in Equations (81)–(85), yields

$$\left[\begin{array}{ccccc}{S}_{11}& S_{12}& 0& S_{14}& S_{15}\\ S_{12}& S_{22}& S_{23}& S_{24}& S_{25}\\ 0& S_{23}& S_{33}& S_{34}& S_{35}\\ S_{14}& S_{24}& S_{34}& S_{44}& S_{45}\\ S_{15}& S_{25}& S_{35}& S_{45}& S_{55}\end{array}\right]\left[\begin{array}{c}{U}_k^V\\ {\Phi }_k^V\\ W_k^V\\ {\Gamma }_k^V\\ Z_k^V\end{array}\right]+\left[\begin{array}{ccccc}{m}_0{\omega }_k^2& -m_1{\omega }_k^2& 0& 0& 0\\ -m_1{\omega }_k^2& m_2{\omega }_k^2& 0& 0& 0\\ 0& 0& m_0{\omega }_k^2& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{U}_k^V\\ {\Phi }_k^V\\ W_k^V\\ {\Gamma }_k^V\\ Z_k^V\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right].$$

(86)

Therefore, the first natural frequency ω1 of the current beam can be solved from the smallest positive root of ω_k² (k = 1) of the Equation (86).

4. Numerical Results

The 50%-50% BaTiO₃-CoFe₂O₄ is adopted for material I [39,42,60,61,62], and the material II is taken to be epoxy [63], as listed in

Table 1. Material properties of the BaTiO₃-CoFe₂O₄ [39] and epoxy [63].

Physical Parameter	Material I	Material II
C11 (GPa)	226	4.889
C44 (GPa)	44.15	1.241
e15 (C/m2)	5.8	0
e31 (C/m2)	−2.2	0
e33 (C/m2)	9.3	0
s11 (10−9C2/(N·m2))	5.64	0
s33 (10−9C2/(N·m2))	6.35	0
q15 (N/(A·m))	275	0
q31 (N/(A·m))	290.15	0
q33 (N/(A·m))	349.85	0
d11 (10−12Ns/(V·C))	5.38	0
d33 (10−12Ns/(V·C))	2740	0
μ11 (10−6Ns2/C2)	297.5	0
μ33 (10−6Ns2/C2)	83.5	0
A11 (N)	11.7484	1.4014
A12 (N)	6.4980	0.6903
ρ (kg/m3)	5550	1180

. Note that the couple stress constants A₁₁ and A₁₂ are estimates based on the formula provided in [14,39]. The magnitude of the uniform load p₀ is equal to 1/2000h N/m, the shear correction factor k_s is 0.8^1/2, and the cross-sectional shape is kept at b = 2h and L = 20h.

In order to verify the correctness of the current model, a comparative study of the deflection of a simply supported microbeam subjected to uniform load between the current model (with Gradient index n = 0) and the model provided by Zhang et al. [39] are plotted in Figure 2. The beam parameters are adopted from Zhang et al. [39].

From Figure 2, it is obvious that the results of the classical and current model are the same as those in Zhang et al. [39]. In addition, this validates the current model and shows that the microstructure effect will always cause the deflection to decrease, as expected.

4.1. Static Bending

Figure 3 shows the distributions of the deformation, axial normal stress, and the electric and magnetic potentials of the current beam. In order to facilitate the observation of the deformation trend of the current beam, the x-component of the displacement vector u of a point (x, y, z) on the beam cross-section has been enlarged by 10 times. In addition, the thickness h is 20 μm, and the gradient index n is 5.

From Figure 3b, it can be observed that the axial normal stress in the middle of the current beam is relatively small, and the axial normal stress at the top of the beam is relatively large. From Figure 3c,d, it is clear that the distributions of electric and magnetic potentials in the current beam are center-symmetrical, and the maximum magnitudes both appear at the center of the beam.

Figure 4 and Figure 5 show the deflections and rotation angles with different thicknesses predicted by current and classical models. The gradient index n is 5. The numerical results for the current model (solid lines) incorporating the couple stress effect (with A₁₁ ≠ 0 and A₁₂ ≠ 0) are directly calculated from Equations (62), (63) and (74), while those for the classical model (dashed lines) are obtained using the same equations but with A₁₁ = A₁₂ = 0.

From Figure 4 and Figure 5, it can be found that the deflections and rotation angles of the current model are always smaller than those of the classical model in all cases. The difference between the results of the current and classical models is obvious when the beam thickness h is small, as expected.

Figure 6 shows the axial normal stress at the beam center (x = L/2) along the thickness direction of the current and classical models. From Figure 6, it is clear that the magnitude of the axial normal stress of the current model is always smaller than that of the classical model. The differences between the axial normal stress predicted by the two models also become smaller with the increase in the thickness h.

Figure 7 and Figure 8 display the electric and magnetic potentials of the FGM simply supported beam with different thickness of the current and classical models. From Figure 7 and Figure 8, it can be observed that the values of electric and magnetic potentials of the current model are always smaller than those of the classical model. When the beam thickness h is small, the differences between the two sets of results are very large. However, the differences become small when the beam thickness increases. This phenomenon also indicates that the microstructure effect is significant for very thin beams.

To illustrate the material inhomogeneity, Figure 9 shows the variation of the maximum deflections w_max (x = L/2) of the MEE-FGM beam with the different gradient index n for h = 20 μm and 20 mm, respectively. From Figure 9a, it can be seen that the maximum deflections w_max increases with the increase of the gradient index—for both current and classical models—and the deflection of the classical model is always larger than that of the current model. From Figure 9b, it is found that when the thickness of the beam is large enough, there is almost no difference in the prediction results of the maximum deflections predicted by the two models, which further confirms that the microstructure effect is only important for very thin beams. In addition, from Figure 9a,b, it is shown that the gradient index n does have a significant effect on the static bending response for all length scales.

Figure 10 shows the variation of the axial normal stress σ_xx(L/2, z) of the current MEE-FGM beam through the thickness for different gradient index n. From Figure 10, it can be found that the axial normal distribution of current MEE-FGM beam is different from that of a homogeneous beam for both h = 20 μm and 20 mm cases. In addition, the axial normal stress of homogenous beams on the geometric central axial (z = 0) is zero, but the zero-valued stresses positions of the current FGM beam are varying with the n. Furthermore, the axial normal stress of a homogeneous beam is linear, while those of current MEE-FGM beam are nonlinear at all length scales.

4.2. Free Vibration

Figure 11 shows the variation of natural frequency (with k = 1) of the MEE-FGM beam of the current and classical models with different beam thickness. From Figure 11, it is obvious that the natural frequencies of both the current and classical models decrease with the thickness increases. The results also show that the current model incorporating the couple stress effect always increases the value of the natural frequency (and thus increased the beam stiffness). When the beam thickness is small enough, the couple stress effect is significant.

Figure 12 shows the variation of the natural frequency ω₁ (k = 1) of the current and classical models with different gradient index n for h = 20 μm and 20 mm. From Figure 12a, it is clear that when the thickness of the beam is small (micro scale), the prediction results of the two models are very different. However, from Figure 12b, when the thickness is big enough (macro scale), the prediction results of the two models are almost the same. In addition, the effect of the gradient index is found to be important for all length scales.

5. Conclusions

Based on the extended modified couple stress theory, a new graded magneto-electro-elastic Timoshenko microbeam model is developed. The new model considers the effects of both three-field coupling and couple stress. The equations of motion and complete boundary conditions of the new microbeam model are determined through a variational approach.

As two direct applications of the new model, the static bending and free vibration properties of a simply supported microbeam subjected to uniformly distributed loads are analytically obtained. For the static bending problem, parametric studies demonstrate that the deflections, rotations, axial normal stresses, electric and magnetic potentials predicted by the current model are all smaller than those of the classical theory. The differences decrease with the increase in the microbeam thickness. For the problem of free vibration, the natural frequency obtained from the current model is found to be higher than that of the classical model. The difference increases as the thickness of the beam decreases. Such a behavior also indicates that the microstructure effect tends to make the graded magneto-electro-elastic microbeam stiffer, and the current model can predict the size effect for magneto-electro-elastic functionally graded microbeam. In addition, it was demonstrated that changing the gradient index significantly affects both the static and vibrational properties of the graded magneto-electro-elastic microbeam at all length scales. These findings are helpful in guiding the engineering design and optimization of graded magneto-electro-elastic materials in MEMS and NEMS devices.