A Nonlinear Finite Element Method for Magnetoelectric Composite and the Study on the Influence of Interfacial Bonding

Magnetoelectric composite material is effective in transferring magnetic field into electric signal. In this paper, a nonlinear finite element method is present to model the magnetoelectric composite of ferroelectric and magnetostrictive material. In the method, the nonlinear and coupling behavior of magnetostrictive material such as Terfenol-D is considered. The nonuniform magnetic, electric, and mechanical field distributions are present. An interfacial transferring coefficient is defined to investigate the performance of interfacial mechanical coupling quantitatively, and the influence of the properties of interfacial bonding material and interfacial cracks on magnetoelectric coefficient is discussed. A new laminate ME composite of curved interface is proposed to overcome weak interfacial bonding.


Introduction
Magnetoelectric (ME) composite materials, composed of magnetostrictive and ferroelectric materials, are effective in converting magnetic field into electric signal.In particular, the magnetic field deforms the magnetostrictive material due to magnetostriction effect, which then stresses the ferroelectric material through interfacial mechanical coupling, and finally the ferroelectric material generates electric field due to its piezoelectric property.ME composite is able to generate a magnetoelectric coefficient of several V/cm-Oe, much higher than single-phase material, which makes it a good candidate for potential applications in sensors, transformers, energy harvesters, and so on [1].
The studies on ME composites can be dated back to 1972 when van Suchtelen proposed the idea of ME conversion through mechanical coupling between two-phase materials [4].Since then, several researchers fabricated such composite and the obtained ME coefficient has been increasing [5][6][7][8][9][10][11][12].At the same time, much attention is also focused on the theoretical and numerical investigations.Harshe et al. developed a theoretical model dealing with the linear behavior of ME composites [13].Nan et al. predicted a giant magnetoelectric effect of Terfenol-D and P (VDF-TrFE) composites using Green's function technique [14,15].The shear-lag and demagnetization effects in laminated ME composites were considered by Chang and Carman [16].Nan et al. discussed the influence of interfacial bonding on the magnetoelectric coefficient [17].The resonance behavior of ME composite was investigated in [18][19][20].Numerically, Liu et al. [21] and Zhou et al. [22] calculated the magnetoelectric effect using finite element method (FEM).In their models, the nonlinear behavior of magnetostrictive material is considered.Linnemann et al. [23] proposed a constitutive model for magnetostrictive and piezoelectric materials and gave FEM examples.Nguyen et al. [24] modeled the nonlinear behavior of magnetic sensor using FEM.In recent years, there are several review articles on ME composites, which give a comprehensive knowledge of the development and current status about these materials [1,25,26].
Compared to analytical models, the finite element method (FEM) has two advantages.First, FEM is able to simulate the magnetic, electric, and mechanical field distribution in the composite regardless of the geometry and configuration, while analytical models simplify the field distribution under certain assumptions which are only valid in certain cases.Second, FEM can adopt nonlinear constitutive relations more conveniently, which is crucial for modeling Mathematical Problems in Engineering ME composite because the magnetostriction of the most widely used magnetostrictive material (e.g., Terfenol-D) is highly nonlinear and magnetomechanical coupling; that is, the magnetostriction depends both on the magnetic field and stresses.
In the existing works, FEM is usually used to predict the ME coefficient.In our opinion, it is also desirable to get knowledge of the nonuniform magnetic, electric, and mechanical field distribution in this composite material, especially for the laminate ME composite with cracks.Because the ferroelectric ceramics usually used in ME composites are quite brittle, the maximum stress in this material should be monitored to assure their reliability.To achieve this goal, the nonlinear and coupling constitutive relation of the magnetostrictive material should be used in FEM.The existing FEM models usually did not consider the coupling behaviors, that is, the influence of stresses on magnetostriction.Exceptionally, in [22] this coupling behavior was fully considered using a constitutive relation which fit well with experimental data.However their FEM example was onedimensional.
In this paper, we propose a finite element model for ME composites, considering the nonlinear behaviors of both magnetostrictive and ferroelectric materials.The nonlinear and coupling constitutive models of magnetostrictive materials have been proposed by several scholars.We adopt the model based on magnetic domain rotation [2,28] in our FEM because it predicts the experimental observed curve well and has profound physical background.For the ferroelectric material, the nonlinearity originates from domain switching under electric field and stresses.In our FEM, this ferroelectric behavior is modeled using a FEM developed by Liu et al. [27,29] which reveals the fully coupled behaviors of electromechanical field and electric domainswitching, and the computation is quite efficient.
The interfacial bonding condition significantly influences the mechanical coupling between layers of ME composites.In most widely used ME composites, the mechanical coupling is achieved through interfacial shearing, which strongly depends on the properties of bonding material [17].The bonding layer of small shear modulus and large thickness leads to weak interfacial coupling and reduced ME coefficient.ME coefficient is also affected by the imperfections of the interfacial bonding, for example, the interfacial cracks.In this paper, we define an interfacial transferring coefficient to describe the effect of interfacial mechanical coupling quantitatively and investigate the influence of stiffness and thickness of bonding material on ME effect using our proposed FEM.The influence of the interfacial cracks on ME coefficient is also obtained.Afterwards we propose a modified ME structure with curved interface and the FEM results show that this structure achieves effective ME conversion even if the interfacial bonding is weak.
We organize this paper as follows.The nonlinear FEM is developed in Sections 2-4.In Section 2, the basic governing equations are given and the FEM formula is derived.The nonlinear behaviors of magnetostrictive material and ferroelectric material are simulated in Sections 3 and 4 separately.
Several numerical examples are present in Sections 5 and 6.
In Section 5, the magnetic, electric, and mechanical field distribution in the typical three-layer ME composite is simulated.Section 6 focuses on the influence of interfacial bonding on the properties of laminate ME composites.Finally, the major conclusions are summarized in Section 7.

Basic Governing Equations and FEM Formula
In our FEM model the material is assumed to have the magnetic, electric, mechanic, magnetomechanical coupling and electromechanical coupling properties.If the actual material lacks a certain property, we set the corresponding parameters and (if necessary) the field variables to be zero.In this paper, only static magnetic and electric field is considered, so that the displacement vector u, magnetic scalar potential , and electric potential  are chosen as basic variables.The basic equations are listed as follows.
Generalized kinematics equations are where   is the strain,   is the magnetic field, and   is the electric field.
The generalized balance equations are where   ,   , and   are the stress, magnetic induction, and electric displacement.  and   are the mechanical body force and electric volume charge.In (2), we neglect the influence of magnetic body force as in [23], and the symmetry of stress tensor is guaranteed.Constitutive equations are *  =   *  +   *  (8) in which   is the elastic stiff tensor,   is the piezoelectric tensor, and   is the dielectric tensor. *  is the sum of the strain caused by magnetostriction (  *  ) and electric domainswitching (  *  ),   is the magnetization caused by magnetic domain rotation, and  *  is the electric displacement caused by domainswitching.
The boundary conditions are summarized below.The given force boundary on   is The given magnetic charge boundary on   and electric charge boundary on   are is the magnetic surface charge, an imaginary quantity, which is introduced in the formulation when a scalar magnetic potential is used [23].  is the electric surface charge, and   represents the normal direction of the surface.
The given displacement, magnetic potential, and electric potential boundary   ,   , and   are Next we derive a potential energy for FEM via the variation principle.Virtual displacement   , magnetic potential , and electric potential  are assumed to satisfy the given displacement and magnetic and electric potential boundary conditions.The generalized balance equations and given force, magnetic, electric charge boundaries are equivalent to the following variation equations: Using Gauss theorem and constitutive equations ( 5) and (6), the desired potential is The variation of (13) Π = 0 is equal to (12).
For simplicity, we only study two dimensional cases in this paper, although it is not difficult to extend to three dimensional problems.Let û, ε, and σ denote the generalized displacement, strain, and stress: The generalized stress and strain are linked by the constitutive relation as The matrix D and matrix C * are reduced from three-dimensional constitutive relations, which are given in Appendix A.
In an element, the generalized displacement can be interpolated as where a is the generalized nodal displacement and N is the shape function as follows: Then following the procedure of FEM, the stiff matrix and force matrix of an element are formed as B  DB d d, where

Coupling Magnetomechanical Field and Magnetic Domain Distribution
In Section 2, the magnetostriction   *  and magnetization   are determined from magnetic domain distribution, which evolve in the existence of magnetic field and stresses.This corresponding coupling constitutive model is proposed by Armstrong [28] and Pei and Fang [2], so we only give a brief introduction here.

Determine Magnetostriction and Magnetization from Domain Distribution.
For each material point, nine representative magnetic domain directions are chosen, which are the eight easy magnetization axes plus the direction of external field.The magnetization and magnetostriction of each magnetic domain are: where  100 and  111 are the saturated magnetostriction measured in [100] and [111] directions, respectively.The magnetization and magnetostriction of a material point are obtained by averaging ( 20) over all representative domain directions so that where   is the volume fraction of the th domain, which is determined from magnetic field and stresses as discussed later in Section 3.2.

Determine Magnetic Domain Distribution from Magnetomechanical Field.
It is assumed in [2,28] that the domain having lower free energy occupies more volume fraction and a quantitative expression is given as [2] where   and  are domain energy and normalization factor, respectively, and  is called energy distribution parameter, whose expression is where , , and  0 are fitting parameters. eq is the Mises equivalent stress.The total free energy of each domain is where   ,  el ,  me ,   , and   are the magnetocrystalline anisotropy energy, elastic energy, magnetoelastic energy, magnetic field energy, and stress energy, respectively.The first three energy terms can be combined into one term with the introduction of where  = 42 is the number of electric domains at each material point.In deriving (27), the Reuss assumption is adopted so that the stress and the electric field are the same for all domains at each material point.

Determine Domainswitching from Electromechanical
Filed.We determine the domainswitching behavior based on double Gibbs free energy criterion.The Gibbs free energy of an electric domain in state  is [30]  (, E, ) = − (  * () : The Gibbs free energy criterion assumes that if the difference of Gibbs energy between two domain states exceeds a critical value, electric domain can switch to the state with lower Gibbs energy.In our program, we first determine the domain state with the lowest Gibbs energy, and then the difference of free energy between the present state and the new state is calculated.If the difference exceeds the critical value, the domain switches to the new state.Next the electric and mechanical fields are recalculated using the newly switched domains and then the Gibbs free energy criterion is used again to check whether this domain rotation is allowed or not.
Summarizing Sections 3 and 4, the solving procedures of our FEM are as follows.
(1) Determine the magnetic, electric, and mechanical fields from an initial guessed domain states by solving the FEM equations in Section 2.
(2) Update the magnetic domain distribution and the electric domain state according to the criteria in Sections 3.2 and 4.2, respectively.
The details of this iterative solving procedure are further illustrated in Figure 1.

The Field Distribution in Laminate ME Composite
Several numerical examples of the proposed FEM are given in this and the following section.The simulated ferroelectric material is PLZT (La/Zr/Ti = 8/65/35) and the magnetostrictive material is Tb 0.3 Dy 0.7 Fe 1.95 compound.Their material constants are listed in Tables 1 and 2, respectively.

Constitutive Curve.
Our program is first validated by the experimental constitutive curves of PLZT and Tb 0.3 Dy 0.7 Fe 1.95 from Hwang et al. [3] and Pei and Fang [2].
From Figure 2, we can see that our numerical model can show the irreversible - and - curves of ferroelectric material and captures the key nature of ferroelectric hysteresis loop and butterfly hysteresis loop.The nonlinear dependence of magnetostriction and magnetization on magnetic field under different stresses are also revealed.All of the four predicted constitutive curves reasonably agree with the experimental results.

The Field Distribution in Three-Layer ME Composite.
Figure 3 schematically shows a typical three-layer laminate Saturation magnetization   (KA/m) 760 Fitting parameter  (J/m 3 ) 15000 Fitting parameter  (J/m 3 ) 30000 Fitting parameter  0 (MPa) 75 ME composite, in which one ferroelectric layer is sandwiched between two magnetostrictive layers.The geometry parameters are chosen as  = 10.0 mm,   = 1.0 mm,   = 1.0 mm, and   = 0.05 mm.The spontaneous polarization of the ferroelectric material is alone  direction and a magnetic field is applied along  direction.The magnetostrictive and ferroelectric materials are bonded through a thin elastic bonding layer whose elastic constants and thickness can be adjusted to model different bonding conditions.In this example, we assume the materials are in perfect bonding condition so we set the shear module of the bonding layer to be 50 GPa.The simulation is carried out in the - plane.
From the simulation results shown in Figure 4, it is obvious that the magnetic, electric, and mechanical fields are not uniform in each layer.Shear stress concentrates at two ends of the interface (Figure 4(a)), while the normal stress   reaches its maximum in the middle of the material (Figure 4(b)).As shown in Figure 5, the maximal shear stress and normal stress increase as the applied magnetic field increases, but not in a proportional way since our FEM is a nonlinear analysis.Analytical models usually assume the stress is uniform along the thickness direction in each layer.Our FEM results show that this assumption is valid only in the middle of the material.

Yes
Next loading Figure 1: FEM program flowchart.The electric domain status is recorded by an integer array ED, while the magnetic domain distribution, that is, the volume fraction of each domain, is represented by a real array MD.

The Effect of the Modulus and the Thickness of Bonding
Layer on ME Coefficient.In order to investigate the influence of the modules and thickness of bonding layer quantitatively, we define an interfacial transferring coefficient to describe the interfacial bonding performance.As shown in Figure 6, at the ends of the composite nominal stress increases from zero to the peak value.In an effective interfacial bonding, the stress should increase as quickly as possible, so the performance of interfacial coupling can be measured by the stress transferred to the piezoelectric layer by a unit length at the ends of a sufficient long composite, namely, d , at the ends.
If the bonding layer is tightly bonded to both the magnetostrictive and the ferroelectric layers without sliding,  is determined as (refer to Appendix B for the detailed derivation) where  *  is the magnetostriction.In the extreme case G ()   /  → ∞, we have which gives rise to the best interfacial bonding condition one can ever achieve.Then the interfacial transferring coefficient is defined as the normalized term In the case of Section 5.2, (32) gives B = 0.95, which is very close to the best interfacial bonding condition.In this section we first reduce the shear module of the bonding layer from 50 GPa to 50 MPa, where B = 0.10 to model a weak interfacial bonding condition.The simulated electric, magnetic and mechanical fields are shown in Figure 7. Comparing with Figure 4, we can see that the stress and generated electric field are much lower.The magnetoelectric coefficient under different bonding conditions is given in Figure 8, showing that a weak interfacial bonding significantly reduces the ME coefficient.

The Effect of Interfacial Cracks on the Properties of ME Composite.
In this example, we investigate the case when the interfacial cracks exist in the structure.Two typical cases are studied.The first one is cracked at one end (Figure 9) while the other is cracked in the middle (Figure 10).Simulation results show that cracks at the ends of the material cause stress concentration near the crack tip (Figures 9(a) and 9(b)) and a decrease in ME coefficient due to the decreased effective interface length to transfer mechanical coupling.However, when the crack is in the middle, its influence is insignificant (Figure 10) and the ME coefficient is not affected.This can be explained by the nature of mechanical coupling through laminated layers.The shearing transfer of stress starts from the ends and the shear stress decreases towards the middle.In fact, in the middle of the interface, the shear stress is almost zero, indicating that the coupling is almost accomplished.Since the middle part of the interface plays little role in mechanical coupling, a crack there has no influence.The effect of crack length on ME coefficient is also investigated.As shown in Figure 11, when cracks locate at the end, ME coefficient decreases almost linearly as the crack grows, while middle cracks do not affect ME coefficient before they grow to the ends.a modified ME structure with curved interface as shown in Figure 12, where the interfacial normal stress plays an important role in the mechanical coupling.The simulation results of the curved structure are shown in Figure 13, and the applied magnetic field and interfacial bonding condition are the same as in Section 6.1.Comparing Figures 13(a) and 13(b) with Figures 7(a) and 7(b), we can conclude that the modified structure achieves mechanical coupling effectively even in poor bonding condition.As shown in Figure 14(a), the interfacial transferring coefficient of the structure with curved interface is higher than the flat interfacial.The ME coefficient of this curved structure in different applied magnetic field is shown in Figure 14(b), suggesting that this structure improves the ME coefficient when interfacial bonding is weak.However, the curved interfacial slightly decreases the ME coefficient in good bonding conditions, because the volume of magnetostrictive layer is decreased in this structure.

Conclusion
A finite element model is proposed for simulating the nonlinear magnetomechanical and electromechanical behavior based on domain evolvement.This model is used to simulate the laminated magnetoelectric composite of magnetostrictive and ferroelectric material.The magnetic, electric, and mechanical fields as well as the magnetoelectric coefficient of a three-layer composite are calculated and discussed.The influence of interfacial bonding is discussed based on the simulation results.The effect of interfacial mechanical coupling is investigated quantitatively by defining an interfacial transferring coefficient, which depends on the shear modulus and thickness of bonding layer.A quantitative relation between the interfacial transferring coefficient and magnetoelectric coefficient is obtained.It is found that the magnetoelectric coefficient strongly depends on the interfacial transferring coefficient, where a poor interfacial bonding results in ineffective mechanical coupling and low magnetoelectric coefficient.To improve the performance of laminated magnetoelectric composite, a modified structure with curved interface is proposed.The simulated result proves that the new structure can achieve good mechanical coupling even when the interfacial bonding is weak.We also investigate the influence of interfacial cracks on magnetoelectric effect.Particularly, interfacial cracks located at the ends of the material decrease the ME effect, while cracks in the middle have little influence.The examples show that our FEM is useful in the designing of magnetoelectric structures.

B. The Interfacial Bonding Performance of Three-Layer ME Composite
Considering the three-layer ME composite in Figure 3, we assume that   ≪   in both the magnetostrictive and the piezoelectric material, so one-dimensional constitutive relations are used: where the superscript () stands for the piezoelectric layer, () stands for the magnetostrictive layer, and () stands for

Figure 4 :
Figure 4: Numerical results on the field distributions of a three-layer laminated composite in perfect bonding condition.(a) Shear stress   , (b) normal stress   , (c) electric field   , (d) magnetization   .The applied magnetic field is 630 Oe.

Figure 5 :
Figure 5: The maximal normal stress in piezoelectric layer and the maximal shear stress increases as the applied magnetic field.

Figure 6 :
Figure 6: The definition of , namely the stress transferred to the piezoelectric layer by a unit length at the ends of the composite.

Figure 7 :
Figure 7: Numerical results on the field distributions of a three-layer laminated composite in poor bonding condition.(a) Shear stress   , (b) normal stress   , (c) electric field   , (d) magnetization   .The ineffective interfacial bonding results in a low stress and electric field the in ferroelectric material.The applied magnetic field is 630 Oe.

Figure 9 :Figure 10 :Figure 11 :Figure 12 :Figure 13 :Figure 14 :
Figure 9: Numerical results on the field distributions of a three-layer laminate composite with an interfacial crack at one end.(a) Shear stress   , (b) normal stress   , (c) electric field   , (d) magnetization   .Stresses concentrate near the tip of the cracks ((a) and (b)).The applied magnetic field is 630 Oe.
4.1.Determine Electric and Mechanical Properties from Electric Domains.In our FEM model, at each material point 42 representative electric domains with the same volume fraction are considered.At each material point, the constitutive equations are written as   =   *  +      +     ,   =  *  +     +