An Inhomogeneous Cell-Based Smoothed Finite Element Method for Free Vibration Calculation of Functionally Graded Magnetoelectroelastic Structures

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 gradedmagnetoelectroelastic (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.


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][6][7][8][9][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][21][22][23][24][25][26] and node-based, edge-based, or facebased smoothed FEMs [27][28][29][30][31][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][34][35][36][37][38][39][40][41][42][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 MEMSbased FGMEE energy harvester are investigated in detail.Final conclusions from the numerical results are drawn in Section 6.

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 are where   ,   , and   denote stress components;   and   are electric displacement components;   and   are magnetic induction components.The geometric equations are where   ,   , 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 are where   ,   , and   are the elastic, dielectric, and magnetic permeability coefficients, respectively;   ,   , and   are piezo-electric, piezomagnetic, and magnetoelectric coefficients, respectively.For FGMEE materials, we have where () is an arbitrary function;  0  and  0  ,  0  ,  0  ,  0  , and  0  are the values on the plane  = 0.

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 as where  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 −ℎ/2 ≤  ≤ ℎ/2 and  is the power law index.The bottom surface of the material ( = −ℎ/2) is   whereas the top surface ( = ℎ/2) 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 as where "eff" stands for effective properties corresponding to a specific value of .Material coefficients of piezoelectric BaTiO

ICS-FEM
The solution domain Ω is discretized into   elements containing   nodes, the approximation displacement u, the approximation electrical potential Φ, and the approximation magnetic potential Ψ.For FGMEE materials, we have where u, Φ, and Ψ are the vectors of node displacement, node electrical potential, and node magnetic potential, respectively; N  , N Φ , and N Ψ are displacement shape, electrical potential shape, and magnetic potential shape functions of ICS-FEM, respectively.N  , N Φ , and N Ψ 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 x  in the smoothing subdomain Ω   , the smoothed strain (x  ), smoothed electric field E(x  ), and smoothed magnetic field H(x  ) are where S(x), E(x), and H(x) are the strain, electric field, and magnetic field in FEM, respectively; (x − x  ) is the constant function: where Substituting ( 12) into (11), we get where Γ   is the boundary of Ω   ; n   , n  Φ , and n  Ψ are the outer normal vector matrices of the boundary:

Gaussian points Outward normal vectors
Eqs. ( 14) can be rewritten as where   is the number of smoothing elements At the Gaussian point of the smoothing boundary x   , ( 17) are rewritten as where    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 x   , 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 as where 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   =   ×   ;   =       ( = 1, 2, 3, 4) is the mass of smoothing element ;  is the smoothing element thickness;   is density of Gaussian integration point in smoothing subdomain ; [C], [], [], [e], [q], and [m] 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    ( = 1, 2, 3, 4) 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 [K eq ]: where The eigenvectors corresponding to Φ and Ψ are given as To study the effect of magnetoelectric constant on system frequencies, we derived [K eq reduced ] by neglecting the magnetoelectric coupling effect.By making [K ΦΨ ] = 0, we get the reduced cell-based finite element equations: The reduced stiffness matrix [K eq reduced ] is To evaluate the effect of PE phase on beam frequency, we derived the stiffness matrix [K eq ΦΦ ] by setting the magnetic potential = 0:

Structural frequency
Matrix used to compute the structural frequency To study the magnetic effect of PM phase on system frequency, we obtained [K eq ΨΨ ] by plugging electric potential to zero in (26):

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.
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  eq because the magnetic effect is marginally higher compared with  uu . ΦΦ coincides with  uu since piezoelectric phase is absent in CoFe 2 O 4 .Similarly,  eq re coincides with  eq 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.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).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  eq ,  eq re , and  eq ΦΦ 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.

S-S Beam. The free vibrations on BaTiO
Secondly, as for BaTiO  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  eq ,  eq re , and  eq ΦΦ are marginal, so the magnetic effect does not largely impact the natural frequencies of FGMEE beams.Meanwhile, the natural frequency of BaTiO

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  = 30 mm,  = 2 mm, and its structure was fabricated with BaTiO 3 FGMEE.
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  eq and  eq ΦΦ agree well with each other since the piezomagnetic phase is absent from the BaTiO 3 FGMEE energy harvester.Moreover,  eq re coincides with  eq of the BaTiO 3 FGMEE energy harvester as the magnetoelectric effect is absent in pure BaTiO 3 materials.The  uu and  eq 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.
The Wilson- method and the equivalent stiffness matrix [K eq ] were employed to solve the dynamic response of the FGMEE energy harvester.The parameters were set as time step = 0.005 s,  = 1.4; without damping; sinewave 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.

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.

Figure 2 :
Figure 2: Variation of the volume fraction function versus the nondimensional thickness /ℎ with varying .

Figure 3 :
Figure 3: Smoothing subdomains and the values of shape functions.

Figure 4 :
Figure 4: Geometry of an FGMEE beam and the coordinates.

Figure 6 :
Figure 6: 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 by ICS-FEM.

Figure 7 :
Figure 7: The first-to third-order modes of clamp-free BaTiO 3 -CoFe 2 O 4 FGMEE beams using different elements with exponential factor  = 5.0 by ICS-FEM.

3 -
CoFe 2 O 4 FGMEE beams were studied by changing the gradient function form in Figure 4.The geometrical parameters were
3 , magnetostrictive CoFe 2 O 4 , and MEE BatiO 3 -CoFe 2 O 4 are given in Table1.Figure2depicts the throughthe-thickness distribution of the volume fraction changing with different values of .For  = 1.0, the variation of effective material property is linear.

Table 2 :
List of f used in the study. sol.