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

Authors: He-Ling Wang; Bin Liu; Dai-Ning Fang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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.


---

## Body

## 1. 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–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–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 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 one-dimensional.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 Sections2–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.


## 2. 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 vectoru, magnetic scalar potential ϕ, and electric potential φ are chosen as basic variables. The basic equations are listed as follows.Generalized kinematics equations are(1)εij=12(ui,j+uj,i),Hi=-ϕ,i,Ei=-φ,i,
where εij is the strain, Hi is the magnetic field, and Ei is the electric field.The generalized balance equations are(2)σij,j+fi=0,(3)Bi,i=0,(4)Di,i=pe,
where σij, Bi, and Di are the stress, magnetic induction, and electric displacement. fi and pe 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
(5)σij=cijkl(εkl-εkl*)-enijEn,(6)Bi=μ0(Hi+Mi),(7)Di=eikl(εkl-εkl*)+kinEn+Di*,(8)εkl*=εklm*+εkle*
in which cijkl is the elastic stiff tensor, enij is the piezoelectric tensor, and kin is the dielectric tensor. εkl* is the sum of the strain caused by magnetostriction (εklm*) and electric domainswitching (εkle*), Mi is the magnetization caused by magnetic domain rotation, and Di* is the electric displacement caused by domainswitching.The boundary conditions are summarized below.The given force boundary onSσ is
(9)σijni=T-j.The given magnetic charge boundary onSB and electric charge boundary on SD are
(10)Bini=σM,Dini=σE.σM is the magnetic surface charge, an imaginary quantity, which is introduced in the formulation when a scalar magnetic potential is used [23]. σE is the electric surface charge, and ni represents the normal direction of the surface.The given displacement, magnetic potential, and electric potential boundarySu, Sϕ, and Sφ are
(11)ui=u-i,ϕ=ϕ-,φ=φ-.Next we derive a potential energy for FEM via the variation principle. Virtual displacementδui, 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:
(12)∫V[(σij,i+fj)δuj+Bi,iδϕ+(Di,i-pe)δφ]dV-∫Sσ(σijni-T-j)δujdS-∫SB(Bini-σM)δϕdS-∫SD(Dini-σE)δφdS=0.Using Gauss theorem and constitutive equations (5) and (6), the desired potential is
(13)Π=∫V(12cijklεijεkl-12μ0HiHi-ekijEkεij-12kijEiEj)dV-∫V[cijklεij*εkl+μ0MiHi-(eijkεjk*-Di*)]dV-∫V(fiui-peφ)dV-∫SσT-iuidS-∫SBσMϕdS-∫SDσEφdS.
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. Letu^, ε^, and σ^ denote the generalized displacement, strain, and stress:
(14)u^=[uxuyφϕ]T,ε^=[εxεyγxy-Ex-Ey-Hx-Hy]T,σ^=[σxσyσxyDxDyBxBy]T.
The generalized stress and strain are linked by the constitutive relation as
(15)σ^=D~ε^+C*.
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(16)u^=Na,
where a is the generalized nodal displacement and N is the shape function as follows:
(17)a=[a1a2⋯am]T,ai=[ux(i)uy(i)φ(i)ϕ(i)],N=[N1N2⋯Nm],Ni=[Ni0000Ni0000Ni0000Ni].
Then following the procedure of FEM, the stiff matrix and force matrix of an element are formed as
(18)Ke=∫V12BTD~Bdxdy,Pe=∫VBTC*dxdy+∫VNTfdxdy+∫SσNTT-dS+∫SBNTp-mdS+∫SDNTp-edS,
where
(19)B=[B1B2⋯Bm],Bi=[∂Ni∂x0∂Ni∂y00000∂Ni∂y∂Ni∂x0000000∂Ni∂x∂Ni∂y0000000∂Ni∂x∂Ni∂y]T,f=[fxfy-pe0]T,T-=[T-xT-y00]T,p-e=[00σE0]T,p-m=[000σM]T.


## 3. Coupling Magnetomechanical Field and Magnetic Domain Distribution

In Section2, the magnetostriction εklm* and magnetization Mi 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.


### 3.1. 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:(20)Mi(k)=Msmi,εij(k)={32λ100(mi2-13),mi=mj3λ111mimj,mi≠mj,
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
(21)Mi=∑k=19PkMi(k),εijm*=∑k=19Pkεij(k),
where Pk is the volume fraction of the kth domain, which is determined from magnetic field and stresses as discussed later in Section 3.2.


### 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]
(22)Pk=Nexp[-Ekω],
where Ek and N are domain energy and normalization factor, respectively, and ω is called energy distribution parameter, whose expression is
(23)ω=A+Bσeqσ0,
where A, B, and σ0 are fitting parameters. σeq is the Mises equivalent stress.The total free energy of each domain is(24)E=EK+Eel+Eme+EH+Eσ,
where EK, Eel, Eme, EH, and Eσ 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 K′=K1+ΔK1≈K1+9/4λ1002(C11-C12)-9/2λ1112C44, where K1 is the magnetocrystalline anisotropy constant and C11, C12, and C44 are elastic constants, so that
(25)EK′=EK+Eel+Eme=K′(m12m22+m22m32+m32m12),EH=-μ0MsHiαi,Eσ=-32λ100(m12σ11+m22σ22+m32σ33)-3λ111(m1m2σ12+m2m3σ23+m1m3σ31).


## 3.1. 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:(20)Mi(k)=Msmi,εij(k)={32λ100(mi2-13),mi=mj3λ111mimj,mi≠mj,
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
(21)Mi=∑k=19PkMi(k),εijm*=∑k=19Pkεij(k),
where Pk is the volume fraction of the kth domain, which is determined from magnetic field and stresses as discussed later in Section 3.2.


## 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]
(22)Pk=Nexp[-Ekω],
where Ek and N are domain energy and normalization factor, respectively, and ω is called energy distribution parameter, whose expression is
(23)ω=A+Bσeqσ0,
where A, B, and σ0 are fitting parameters. σeq is the Mises equivalent stress.The total free energy of each domain is(24)E=EK+Eel+Eme+EH+Eσ,
where EK, Eel, Eme, EH, and Eσ 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 K′=K1+ΔK1≈K1+9/4λ1002(C11-C12)-9/2λ1112C44, where K1 is the magnetocrystalline anisotropy constant and C11, C12, and C44 are elastic constants, so that
(25)EK′=EK+Eel+Eme=K′(m12m22+m22m32+m32m12),EH=-μ0MsHiαi,Eσ=-32λ100(m12σ11+m22σ22+m32σ33)-3λ111(m1m2σ12+m2m3σ23+m1m3σ31).


## 4. Coupling Electromechanical Filed and Electric Domainswitching

### 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(26)εij=εije*+MijknEσkn+dkijEk,Dj=Dj*+djknσkn+kjkσEk
in which the elastic moduli MijknE, piezoelectric moduli dkij, dielectric moduli kjkσ, remnant strain εije*, and electric displacement Dj* are calculated by averaging corresponding values of each domain:
(27)MijknE=1m∑l=1mMijknE(l),dijk=1m∑l=1mdijk(l),kijσ=1m∑l=1mkijσ(l),εije*=1m∑l=1mεij*(l),Dj*=1m∑l=1mDj*(l),
where m=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.


### 4.2. 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 stateS is [30]
(28)g(σ,E,S)=-(εe*(S):σ+D*(S)·E+12σ:ME(S):σhhhh+12E·kσ(S)·E+E·d(S):σ).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 Sections3 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 Section2.(2)
Update the magnetic domain distribution and the electric domain state according to the criteria in Sections3.2 and 4.2, respectively.(3)
Repeat steps(1)-(2) until both the magnetic and electric domains do not change.The details of this iterative solving procedure are further illustrated in Figure1.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.


## 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(26)εij=εije*+MijknEσkn+dkijEk,Dj=Dj*+djknσkn+kjkσEk
in which the elastic moduli MijknE, piezoelectric moduli dkij, dielectric moduli kjkσ, remnant strain εije*, and electric displacement Dj* are calculated by averaging corresponding values of each domain:
(27)MijknE=1m∑l=1mMijknE(l),dijk=1m∑l=1mdijk(l),kijσ=1m∑l=1mkijσ(l),εije*=1m∑l=1mεij*(l),Dj*=1m∑l=1mDj*(l),
where m=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.


## 4.2. 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 stateS is [30]
(28)g(σ,E,S)=-(εe*(S):σ+D*(S)·E+12σ:ME(S):σhhhh+12E·kσ(S)·E+E·d(S):σ).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 Sections3 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 Section2.(2)
Update the magnetic domain distribution and the electric domain state according to the criteria in Sections3.2 and 4.2, respectively.(3)
Repeat steps(1)-(2) until both the magnetic and electric domains do not change.The details of this iterative solving procedure are further illustrated in Figure1.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.


## 5. 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 Tb0.3Dy0.7Fe1.95 compound. Their material constants are listed in Tables 1 and 2, respectively.Table 1
Parameters of PLZT [27].
Name of parameters (unit)
Value
Young’s module (GPa)
68
Poisson ratio
0.3
Dielectric rate (farad/m)
5.625
×
10
-
8
Macro piezoelectric rated-333 (m/V)
1.188
×
10
-
9
Macro piezoelectric rated-311 (m/V)
-
5.5
×
10
-
10
Macro piezoelectric rated-113 (m/V)
1.856
×
10
-
9
Coercive fieldEc (MV/m)
0.36
Macroremnant strainε-11R, ε-22R
-
7.2
×
10
-
4
Macroremnant strainε-33R
1.44
×
10
-
3
Remnant polarizationPR (C/m2)
0.25Table 2
Parameters of Tb0.3Dy0.7Fe1.95 [2].
Name of parameters (unit)
Value
Elastic constantC11 (GPa)
141
Elastic constantC12 (GPa)
64.8
Magnetocrystalline  anisotropy constantK1 (J/m3)
60000
Magnetostriction  along[100]λ100
90
×
10
-
6
Magnetostriction   along[111]λ111
1640
×
10
-
6
Saturation  magnetizationMs (KA/m)
760
Fitting  parameterA (J/m3)
15000
Fitting  parameterB  (J/m3)
30000
Fitting  parameterσ0 (MPa)
75


### 5.1. Constitutive Curve

Our program is first validated by the experimental constitutive curves of PLZT and Tb0.3Dy0.7Fe1.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 D-E and ε-E 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.Predicted and experimental constitutive curves: (a) ferroelectric hysteresis loop, (b) butterfly hysteresis loop, (c) magnetization of Terfenol-D, and (d) magnetostriction of Terfenol-D under different prestresses. The experimental data are from [2, 3].
(a)
(b)
(c)
(d)


### 5.2. The Field Distribution in Three-Layer ME Composite

Figure3 schematically shows a typical three-layer laminate ME composite, in which one ferroelectric layer is sandwiched between two magnetostrictive layers. The geometry parameters are chosen as L=10.0 mm, tm=1.0 mm, tp=1.0 mm, and tb=0.05 mm. The spontaneous polarization of the ferroelectric material is alone y direction and a magnetic field is applied along x 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 x-y plane.Figure 3
Schematic illustration of the three-layer laminate ME composite.From the simulation results shown in Figure4, 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 σx 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.Figure 4
Numerical results on the field distributions of a three-layer laminated composite in perfect bonding condition. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx. The applied magnetic field is 630 Oe.Figure 5
The maximal normal stress in piezoelectric layer and the maximal shear stress increases as the applied magnetic field.


## 5.1. Constitutive Curve

Our program is first validated by the experimental constitutive curves of PLZT and Tb0.3Dy0.7Fe1.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 D-E and ε-E 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.Predicted and experimental constitutive curves: (a) ferroelectric hysteresis loop, (b) butterfly hysteresis loop, (c) magnetization of Terfenol-D, and (d) magnetostriction of Terfenol-D under different prestresses. The experimental data are from [2, 3].
(a)
(b)
(c)
(d)


## 5.2. The Field Distribution in Three-Layer ME Composite

Figure3 schematically shows a typical three-layer laminate ME composite, in which one ferroelectric layer is sandwiched between two magnetostrictive layers. The geometry parameters are chosen as L=10.0 mm, tm=1.0 mm, tp=1.0 mm, and tb=0.05 mm. The spontaneous polarization of the ferroelectric material is alone y direction and a magnetic field is applied along x 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 x-y plane.Figure 3
Schematic illustration of the three-layer laminate ME composite.From the simulation results shown in Figure4, 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 σx 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.Figure 4
Numerical results on the field distributions of a three-layer laminated composite in perfect bonding condition. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx. The applied magnetic field is 630 Oe.Figure 5
The maximal normal stress in piezoelectric layer and the maximal shear stress increases as the applied magnetic field.


## 6. The Influence of Interfacial Bonding on ME Effect

### 6.1. 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 Figure6, 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,
(29)B=dσx(p)dx,attheends.
If the bonding layer is tightly bonded to both the magnetostrictive and the ferroelectric layers without sliding, B is determined as (refer to Appendix B for the detailed derivation)
(30)B=2Ex(m)Ex(p)tm2Ex(m)tm+Ex(p)tpε-x*×2Ex(m)tm+Ex(p)tpEx(m)Ex(p)tptm[tbGxy(b)+16tpGxy(p)+13tmGxy(m)]-1,
where ε-x* is the magnetostriction. In the extreme case Gxy(b)/tb→∞, we have
(31)Bperfect=2Ex(m)Ex(p)tm2Ex(m)tm+Ex(p)tpε-x*×2Ex(m)tm+Ex(p)tpEx(m)Ex(p)tptm[16tpGxy(p)+13tmGxy(m)]-1
which gives rise to the best interfacial bonding condition one can ever achieve. Then the interfacial transferring coefficient is defined as the normalized term
(32)B^=BBperfect=tp/Gxy(p)+2(tm/Gxy(m))6(tb/Gxy(b))+tp/Gxy(p)+2(tm/Gxy(m)).Figure 6
The definition ofB, namely the stress transferred to the piezoelectric layer by a unit length at the ends of the composite.In the case of Section5.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
(33)αE=δEδH
under different bonding conditions is given in Figure 8, showing that a weak interfacial bonding significantly reduces the ME coefficient.Figure 7
Numerical results on the field distributions of a three-layer laminated composite in poor bonding condition. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx. The ineffective interfacial bonding results in a low stress and electric field the in ferroelectric material. The applied magnetic field is 630 Oe.Figure 8
Influence of interfacial bonding on ME coefficient. In the simulation, interfacial bonding is controlled by the shear moduleGb and thickness tb of the bonding layer. Large value of Gb/tb corresponds to good interfacial bonding.


### 6.2. 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 (Figure9) 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.Figure 9
Numerical results on the field distributions of a three-layer laminate composite with an interfacial crack at one end. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx. Stresses concentrate near the tip of the cracks ((a) and (b)). The applied magnetic field is 630 Oe.Figure 10
Numerical results on the field distributions of a three-layer laminated composite with an interfacial crack in the middle. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx. Stresses do not concentrate near the tip of the cracks. The applied magnetic field is 630 Oe.The effect of crack length on ME coefficient is also investigated. As shown in Figure11, 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.Figure 11
The ME coefficient as a function of the interfacial crack length.


### 6.3. Laminate ME Composite with Curved Interface

In order to overcome the influence of interfacial bonding condition on the mechanical coupling between adjacent layers, we propose a modified ME structure with curved interface as shown in Figure12, 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.Figure 12
The sketch of ME composite of curved interface.Figure 13
Numerical results of a three-layer laminated composite of a curved interface. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx.(a) Interfacial transferring coefficient as function ofGb/tb, (b) ME coefficient under different bonding conditions in the applied magnetic field of 630 Oe. Laminate ME composite of curved interface improves the ME coefficient when the interfacial bonding is poor but slightly decreases the ME coefficient in good interfacial bonding condition.
(a)
(b)


## 6.1. 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 Figure6, 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,
(29)B=dσx(p)dx,attheends.
If the bonding layer is tightly bonded to both the magnetostrictive and the ferroelectric layers without sliding, B is determined as (refer to Appendix B for the detailed derivation)
(30)B=2Ex(m)Ex(p)tm2Ex(m)tm+Ex(p)tpε-x*×2Ex(m)tm+Ex(p)tpEx(m)Ex(p)tptm[tbGxy(b)+16tpGxy(p)+13tmGxy(m)]-1,
where ε-x* is the magnetostriction. In the extreme case Gxy(b)/tb→∞, we have
(31)Bperfect=2Ex(m)Ex(p)tm2Ex(m)tm+Ex(p)tpε-x*×2Ex(m)tm+Ex(p)tpEx(m)Ex(p)tptm[16tpGxy(p)+13tmGxy(m)]-1
which gives rise to the best interfacial bonding condition one can ever achieve. Then the interfacial transferring coefficient is defined as the normalized term
(32)B^=BBperfect=tp/Gxy(p)+2(tm/Gxy(m))6(tb/Gxy(b))+tp/Gxy(p)+2(tm/Gxy(m)).Figure 6
The definition ofB, namely the stress transferred to the piezoelectric layer by a unit length at the ends of the composite.In the case of Section5.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
(33)αE=δEδH
under different bonding conditions is given in Figure 8, showing that a weak interfacial bonding significantly reduces the ME coefficient.Figure 7
Numerical results on the field distributions of a three-layer laminated composite in poor bonding condition. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx. The ineffective interfacial bonding results in a low stress and electric field the in ferroelectric material. The applied magnetic field is 630 Oe.Figure 8
Influence of interfacial bonding on ME coefficient. In the simulation, interfacial bonding is controlled by the shear moduleGb and thickness tb of the bonding layer. Large value of Gb/tb corresponds to good interfacial bonding.


## 6.2. 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 (Figure9) 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.Figure 9
Numerical results on the field distributions of a three-layer laminate composite with an interfacial crack at one end. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx. Stresses concentrate near the tip of the cracks ((a) and (b)). The applied magnetic field is 630 Oe.Figure 10
Numerical results on the field distributions of a three-layer laminated composite with an interfacial crack in the middle. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx. Stresses do not concentrate near the tip of the cracks. The applied magnetic field is 630 Oe.The effect of crack length on ME coefficient is also investigated. As shown in Figure11, 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.Figure 11
The ME coefficient as a function of the interfacial crack length.


## 6.3. Laminate ME Composite with Curved Interface

In order to overcome the influence of interfacial bonding condition on the mechanical coupling between adjacent layers, we propose a modified ME structure with curved interface as shown in Figure12, 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.Figure 12
The sketch of ME composite of curved interface.Figure 13
Numerical results of a three-layer laminated composite of a curved interface. (a) Shear stressτxy, (b) normal stress σx, (c) electric field Ey, (d) magnetization Mx.(a) Interfacial transferring coefficient as function ofGb/tb, (b) ME coefficient under different bonding conditions in the applied magnetic field of 630 Oe. Laminate ME composite of curved interface improves the ME coefficient when the interfacial bonding is poor but slightly decreases the ME coefficient in good interfacial bonding condition.
(a)
(b)


## 7. 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.



---

*Source: 197940-2013-10-30.xml*