#### Abstract

To overcome the overstiffness and imprecise magnetoelectroelastic coupling effects of finite element method (FEM), we present an inhomogeneous cell-based smoothed FEM (ICS-FEM) of functionally graded magnetoelectroelastic (FGMEE) structures. Then the ICS-FEM formulations for free vibration calculation of FGMEE structures were deduced. In FGMEE structures, the true parameters at the Gaussian integration point were adopted directly to replace the homogenization in an element. The ICS-FEM provides a continuous system with a close-to-exact stiffness, which could be automatically and more easily generated for complicated domains, thus significantly decreasing the numerical error. To verify the accuracy and trustworthiness of ICS-FEM, we investigated several numerical examples and found that ICS-FEM simulated more accurately than the standard FEM. Also the effects of various equivalent stiffness matrices and the gradient function on the inherent frequency of FGMEE beams were studied.

#### 1. Introduction

Functionally graded magnetoelectroelastic (FGMEE) materials are generally multiphase composites with continuously varying mechanical properties. FGMEE materials can convert magnetic, electric, and mechanical energy from one type into another and have received wide attention recently due to their electroelastic, magnetoelastic, and electromagnetic coupling effects [1, 2]. Therefore, FGMEE materials have been adopted in various smart structures, such as magnetic field probes, smart vibration sensors, optoelectronic devices, and medical ultrasonic transducers [3, 4]. The smart FGMEE structures are commonly fabricated in the beam pattern. However, to better apply FGMEE beams, researchers must analyse the statics and free vibration, and to predict the coupled response of FGMEE beams for practical applications, they should accurately calculate the properties of free vibrations.

Several computational techniques were proposed to investigate the electroelastic, magnetoelastic, and electromagnetic coupling effects of smart structures, such as finite element method (FEM), mesh-free method, and scaled boundary FEM [5–10]. Bhangale and Ganesan analyzed the static behaviors of linear anisotropic FGMEE plates using semianalytical FEM and investigated the free vibration of FGMEE plates and cylindrical shells [11, 12]. A layerwise partial mixed FEM was proposed to model MEE plates [13]. Phoenix et al. analysed the static and dynamic behaviors of coupled MEE plates using FEM with the Reissner mixed variational theorem [14]. Buchanan used FEM to study the free vibrations of infinite magnetoelectroelastic cylinders [15]. However, these FEMs overestimated the stiffness of solid structures and were limited by low accuracy. Therefore, Sladek et al. proposed a mesh-free method to more accurately study the static behavior of a circular FGMEE plate [16], but this method reduced the computational effectiveness. Recently, Liu et al. solved the deformations of a nonuniform MEE plate using scaled boundary FEM [17]. By incorporating the nonlocal theory into scaled boundary FEM, Ke and Wang more accurately and effectively studied the free vibrations of MEE beams [18]. However, the effectiveness of scaled boundary FEM is still low and should be improved.

FEM, as a powerful computational tool to investigate MEE coupling behaviors, yet overestimates the stiffness of FGMEE structures, which may result in locking behavior and inaccurate eigenvalue solutions [19]. To overcome these limitations, a series of cell-based smoothed FEM (CS-FEM) [20–26] and node-based, edge-based, or face-based smoothed FEMs [27–32] were proposed. In recent years, many CS-FEM-based formulations were proposed. Moreover, CS-FEM does not require the shape function derivatives or high generosity of program and is insensitive to mesh distortion because of the absence of isoparametric mapping.

CS-FEM has been successfully extended into dynamical control of piezoelectric sensors and actuators, topological optimization of linear piezoelectric micromotor, and analysis of static behaviors, frequency, and defects of piezoelectric structures [33–43]. Due to its versatility, CS-FEM becomes a simple and effective numerical tool to solve numerous electric and mechanical physical problems. However, the application of CS-FEM to investigate MEE properties is still a challenge.

In this work, free vibrations of FGMEE structures were studied. Inhomogeneous CS-FEM (ICS-FEM) for FGMEE materials was formulated by incorporating gradient smoothing into the standard FEM for the multi-physics field of FGMEE. Then the equations of free vibration computation were deduced under the multi-physics coupling field for FGMEE materials. Finally, FGMEE beams in the functional gradient exponential form or power law form were calculated under different boundary conditions. ICS-FEM outperformed FEM when compared with the reference solution.

This paper is organized as follows: Section 2 introduces the basic formulations for FGMEE materials. Section 3 briefly describes the properties of the FGMEE materials. Section 4 contains the detailed formulation of ICS-FEM. In Section 5, two numerical examples and the model of typical MEMS-based FGMEE energy harvester are investigated in detail. Final conclusions from the numerical results are drawn in Section 6.

#### 2. Basic Formulations

The material properties of a functionally graded material (FGM) plate vary continuously, which is an advantage over the discontinuity across adjoining layers in a laminated plate. The wide range of engineering applications of FGM has attracted many scientists to investigate the behaviors of FGM.

Considering the transverse isotropy of the FGMEE medium [9, 44] and for the plane stress problem, we set stress components = = = 0, electric displacement component = 0, and magnetic induction component = 0. The geometric parameters and the chosen Cartesian coordinate system are illustrated in Figure 1.

The basic formulations for MEE materials include equilibrium equations, geometric equations, and constitutive equations. The equilibrium equations arewhere , and denote stress components; and are electric displacement components; and are magnetic induction components.

The geometric equations arewhere , , and denote strain components; and are displacement components; and are electric field components; is electrical potential; and are magnetic field components; is magnetic potential.

The constitutive equations arewhere , , and are the elastic, dielectric, and magnetic permeability coefficients, respectively; , , and are piezo-electric, piezomagnetic, and magnetoelectric coefficients, respectively. For FGMEE materials, we havewhere is an arbitrary function; and , , , , and are the values on the plane = 0.

#### 3. FGMEE Materials

An FGMEE structure is characterized by the high heterogeneity of material properties with a distribution prescribing the volume fractions of constituent phases. For particular analysis, it is functional to idealize them as continua with smooth gradual variation of material properties in the spatial coordinates. Hence, the proper micromechanical model should be able to characterize the material property distribution of a system in accurate sense.

Previous literatures focus on two types of gradation methods widely applied to solve many problems. Among various methods for composites, some are also used for FGMEE materials, including the exponential and Voigt rule of mixture scheme.

For FGMEE materials with exponential variation in the thickness direction (-direction), (4) can be rewritten aswhere is the exponential factor governing the degree of -direction gradient, is the thickness, the superscript 0 indicates the -independent coefficients, and = 0 in homogeneous MEE materials.

The volume fraction of an FGMEE structure across the thickness direction is assumed as a simple power law type as follows:where and is the power law index. The bottom surface of the material () is whereas the top surface () is . The total volume of the constituents should be

Based on (6) and (7), the effective material property is defined as follows: where “MC” is general notation for material property. With (3), the effective coefficients can be written aswhere “eff” stands for effective properties corresponding to a specific value of . Material coefficients of piezoelectric BaTiO_{3}, magnetostrictive CoFe_{2}O_{4}, and MEE BatiO_{3}-CoFe_{2}O_{4} are given in Table 1. Figure 2 depicts the through-the-thickness distribution of the volume fraction changing with different values of . For = 1.0, the variation of effective material property is linear.

#### 4. ICS-FEM

The solution domain is discretized into elements containing nodes, the approximation displacement , the approximation electrical potential , and the approximation magnetic potential . For FGMEE materials, we havewhere** u**, , and are the vectors of node displacement, node electrical potential, and node magnetic potential, respectively; , , and are displacement shape, electrical potential shape, and magnetic potential shape functions of ICS-FEM, respectively. , , and were expressed in similar shape functions. Four-node element divided into four smoothing subdomains [27], field nodes, edge smoothing nodes, center smoothing nodes, edge Gaussian point, outer normal vector distribution, and shape function values are shown in Figure 3.

At any point in the smoothing subdomain , the smoothed strain , smoothed electric field , and smoothed magnetic field arewhere , , and are the strain, electric field, and magnetic field in FEM, respectively; is the constant function:where

Substituting (12) into (11), we getwhere is the boundary of ; , , and are the outer normal vector matrices of the boundary:

Eqs. (14) can be rewritten aswhere is the number of smoothing elements

At the Gaussian point of the smoothing boundary , (17) are rewritten aswhere is the length of the smoothing boundary; is the total number of boundaries for each smoothing subdomain.

As for the essential difference, FEM has to derive the shape function matrix of the element, but ICS-FEM avoids this step and simply uses the shape function at , which reduces the requirement for continuity of the shape function and improves the accuracy and convergence by the use of gradient smoothing.

The thermodynamic potential of a 2D FGMEE problem is given aswhere** S**,** E,** and** H** are independent variables of strain, electric field, and magnetic field, respectively.

By applying (1) into (19), we get the variational expression of MEE plane:

By minimizing (20) for nodal variables of shape functions for strain-displacement, electric field–electric potential, and magnetic field–magnetic potential, we get the ICS-FEM equations for MEE plane:where is the eigenvalues.

Different elemental stiffness matrices used for FGMEE beams are expressed as follows:where ; is the mass of smoothing element ; is the smoothing element thickness; is density of Gaussian integration point in smoothing subdomain ; , , , , , and are the matrices of elastic constant, dielectric coefficient, magnetic permeability, piezomagnetic coefficient, piezomagnetic coefficient, and magnetoelectric coefficient, respectively.

The inhomogeneous smoothing element was adopted to calculate its stiffness matrix. Because the parameters of four smoothing subdomains differed in element , the actual parameters at the Gaussian integration point were taken directly to simulate the changes of material property in each element.

By eliminating the terms of electric and magnetic potentials using a condensation technique, we get the equivalent stiffness matrix :where

The eigenvectors corresponding to and are given as

To study the effect of magnetoelectric constant on system frequencies, we derived by neglecting the magnetoelectric coupling effect.

By making , we get the reduced cell-based finite element equations:

The reduced stiffness matrix is

To evaluate the effect of PE phase on beam frequency, we derived the stiffness matrix by setting the magnetic potential = 0:

To study the magnetic effect of PM phase on system frequency, we obtained by plugging electric potential to zero in (26):

#### 5. Results and Discussion

##### 5.1. C-F Beam

The free vibrations on FGMEE beams were calculated by changing the exponential factor (Figure 4). The material properties of FGMEE beams were governed by the -direction exponential variation. The following geometrical parameters were considered: length = 0.3 m and width = 0.02 m with the assumption of plane stress. Boundary conditions were = = = = 0 at the clamped end. Table 2 gives the various structural frequencies in the study.

Firstly, the convergence of ICS-FEM was verified by using BaTiO_{3}-CoFe_{2}O_{4} FGMEE beams, with properties listed in Table 1. The natural frequencies of these beams were calculated using ICS-FEM with different meshes (30 × 2, 60 × 4, 120 × 8, 150 × 10) (Figure 5). The simulation results with different meshes agree well, which prove the good convergence of ICS-FEM. The first- to third-order modes of clamp-free BaTiO_{3}–CoFe_{2}O_{4} FGMEE beams using different elements with exponential factor = 1.0 and 5.0 were calculated by ICS-FEM, and the results were summarized in Figures 6 and 7, respectively. It was found that the first- to third-order modes of BaTiO_{3}–CoFe_{2}O_{4} FGMEE beams in the same gradient distribution were basically not affected by equivalent stiffness matrix or mesh number. the first- to third-order modes of BaTiO_{3}–CoFe_{2}O_{4} FGMEE beams were basically not different between = 1.0 and 5.0.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

Secondly, the free vibration frequencies of CoFe_{2}O_{4} FGMEE beams were studied by both ICS-FEM and FEM using extremely irregular elements, with domain discretization shown in Figure 8. The frequencies of clamp-free CoFe_{2}O_{4} FGMEE beams with different values of exponential factor are shown in Figure 9; the first eleven natural frequencies calculated by ICS-FEM are smaller than those calculated by FEM. The validity of ICS-FEM is verified by the agreements between the calculations and the reference solutions. The shape of quadrilateral element in FEM cannot be severely distorted but was eliminated in ICS-FEM. ICS-FEM abstains from calculating the derivative of the shape functions of an element, and the area integral of the solution domain is converted to the boundary integral. The stiffness of FGMEE structures is improved because ICS-FEM does not require continuity of the shape function. The ICS-FEM provides a continuous system with a close-to-exact stiffness, which could be automatically and more easily generated for complicated domains, thus significantly decreasing the numerical error. The free vibration of CoFe_{2}O_{4} FGMEE beam, a pure CoFe_{2}O_{4} material without piezoelectric or magnetoelectric material coefficients, influences structural frequency because the magnetic effect is marginally higher compared with . coincides with since piezoelectric phase is absent in CoFe_{2}O_{4}. Similarly, coincides with of the CoFe_{2}O_{4} FGMEE beam as the magnetoelectric effect is absent in pure CoFe_{2}O_{4} FGMEE beam. The natural frequencies of CoFe_{2}O_{4} FGMEE beams increase with the rise of the exponential factor.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

Finally, the comparison of calculation time between ICS-FEM and FEM at Intel (R) Xeon (R) CPU E3-1220 v3 @ 3.10 GHz, 16 G RAM is shown in Figure 10, with the element number of 60, 240, 960, and 1500. As showed in Figure 10, the time required to solve algebraic equations by ICS-FEM is similar to that of FEM. Because the stiffness construction of ICS-FEM is based on smoothing cells inside each element, no coupling occurs between nodal degrees-of-freedom that are the distance of up to two elements. In other words, the bandwidth of ICS-FEM stiffness matrix is the same as that of FEM. Nevertheless, ICS-FEM is more effective in terms of generalized displacement (including displacement, electrical potential and magnetic potential) and computational efficiency (computation time for the same accuracy).

##### 5.2. S-S Beam

The free vibrations on BaTiO_{3}–CoFe_{2}O_{4} FGMEE beams were studied by changing the gradient function form in Figure 4. The geometrical parameters were the same as the C-F beam. The simply supported boundary conditions were used: at (, ) and at (, ).

Firstly, the free vibration frequencies for BaTiO_{3}–CoFe_{2}O_{4} FGMEE beams were calculated by both ICS-FEM with 60 × 4 meshes and FEM with 240 × 16 meshes (Figure 11). Results show the first eleven natural frequencies calculated by ICS-FEM are closer to the reference solutions than those calculated by FEM, indicating ICS-FEM is more efficient than FEM due to the reduced number of meshes. The differences in natural frequencies , , and are marginal, so the magnetic effect only slightly impacts the natural frequencies of FGMEE beams. The natural frequencies increase with the rise of the exponential factor.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

Secondly, as for BaTiO_{3}–CoFe_{2}O_{4} FGMEE beams with power law type (, 1.0, 5.0, 10.0, 20.0), the bottom surface () is BaTiO_{3}–CoFe_{2}O_{4} whereas the top surface () is CoFe_{2}O_{4}. The free vibration frequencies for FGMEE beams calculated by ICS-FEM with 60 × 4 meshes and FEM with 240 × 16 meshes are shown in Figure 12. Results show the first eleven natural frequencies calculated by ICS-FEM are closer to the reference solutions than those calculated by FEM. Also the differences in , , and are marginal, so the magnetic effect does not largely impact the natural frequencies of FGMEE beams. Meanwhile, the natural frequency of BaTiO_{3}–CoFe_{2}O_{4} FGMEE beams is between those of BaTiO_{3}–CoFe_{2}O_{4} MEE beams and CoFe_{2}O_{4} MEE beams.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

##### 5.3. Typical MEMS-Based FGMEE Energy Harvester

The model of FGMEE energy harvester developed by ICS-FEM is shown in Figure 13. The free vibrations on the FGMEE energy harvester were studied by changing the exponential factor. The geometrical parameters were mm, mm, and its structure was fabricated with BaTiO_{3} FGMEE.

**(a)**

**(b)**

The free vibration frequencies for the FGMEE energy harvester calculated by ICS-FEM with 60 × 4 meshes and FEM with 240 × 16 meshes are shown in Figure 14. The first eleven natural frequencies calculated by ICS-FEM are closer to the reference solutions than those calculated by FEM, indicating ICS-FEM is more efficient than FEM owing to the reduced number of meshes. The ICS-FEM does not take the derivative of the shape functions of the element and can be much easily generated automatically for complicated domains, thus significantly decreasing the numerical errors. The natural frequencies and agree well with each other since the piezomagnetic phase is absent from the BaTiO_{3} FGMEE energy harvester. Moreover, coincides with of the BaTiO_{3} FGMEE energy harvester as the magnetoelectric effect is absent in pure BaTiO_{3} materials. The and are very close, so the piezoelectric effect only slightly affects the natural frequencies of the energy harvester, which increase with the rise of the exponential factor.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

The Wilson-*θ* method and the equivalent stiffness matrix [] were employed to solve the dynamic response of the FGMEE energy harvester. The parameters were set as time step = 0.005 s, ; without damping; sine-wave transient load with a time period of 2 s; 4 cycles of loading (Figure 15). The dynamic behaviors for the harvester calculated by ICS-FEM with 60 × 4 meshes and FEM with 240 × 16 meshes are shown in Figures 16 and 17, respectively. The temporal variations of displacement and electric potential calculated by ICS-FEM are closer to the reference solutions than those by FEM, which validate the accuracy of ICS-FEM. The temporal variations of and Φ of the FGMEE energy harvester decrease with the increase of the exponential factor when the material properties of the harvester are in exponential distribution.

#### 6. Conclusions

The free vibrations on FGMEE structures were studied. Firstly, ICS-FEM for FGMEE materials was formulated by incorporating gradient smoothing into the FEM-based computation for the FGMEE multi-physics field. Then the equations of free vibration computation were deduced for the multi-physics coupling field of FGMEE materials. Finally, the FGMEE beams were calculated with functional gradient exponential form or power law form under different boundary conditions.

(i) ICS-FEM reduced the systematic stiffness of the finite element, which improved the computational accuracy compared with FEM under the same element number. ICS-FEM was more efficient than FEM in terms of computation time for the same accuracy.

(ii) Due to the material property changes in each smoothing element, the true parameters at the Gaussian integration point were adopted directly. ICS-FEM avoided the derivative of the shape functions but only transformed the area integral to the boundary integral in the solution domain, which omitted the requirement of continuity of the shape function.

(iii) The magnetic effect slightly influenced the natural frequencies of FGMEE beams, which increased with the increase of the exponential factor when the material properties of FGMEE beams were under exponential distribution. The natural frequency of FGMEE beams lied in between those of the FGMEE beams using the bottom surface of materials and the FGMEE beams using the upper surface, when the material properties of FGMEE beams were under the power law distribution.

(iv) The natural frequencies and general displacements of the FGMEE energy harvester developed by ICS-FEM were more accurate compared with FEM, owing to the reduced number of meshes.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

Special thanks are due to Professor Guirong Liu for the S-FEM Source Code in http://www.ase.uc.edu/~liugr/software.html. This work was financially supported by the National Natural Science Foundation of China (Grant no. 11502092), Jilin Provincial Department of Science and Technology Fund Project (Grant nos. 20160520064JH, 20170101043JC), and the Fundamental Research Funds for the Central Universities.