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

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.

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

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 vector 𝐮, magnetic scalar potential 𝜙, and electric potential 𝜑 are chosen as basic variables. The basic equations are listed as follows.

Generalized kinematics equations are 𝜀_𝑖𝑗=12𝑢_𝑖,𝑗+𝑢_𝑗,𝑖,𝐻_𝑖=−𝜙_,𝑖,𝐸_𝑖=−𝜑_,𝑖, (1) where 𝜀_𝑖𝑗 is the strain, 𝐻_𝑖 is the magnetic field, and 𝐸_𝑖 is the electric field.

The generalized balance equations are 𝜎_{𝑖𝑗,𝑗}+𝑓_𝑖𝐵=0,_𝑖,𝑖𝐷=0,_𝑖,𝑖=𝑝^𝑒, (2)(3)(4) 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 𝜎_𝑖𝑗=𝑐_{𝑖𝑗𝑘𝑙}𝜀_𝑘𝑙−𝜀^∗_𝑘𝑙−𝑒_𝑛𝑖𝑗𝐸_𝑛,𝐵_𝑖=𝜇₀𝐻_𝑖+𝑀_𝑖,𝐷_𝑖=𝑒_𝑖𝑘𝑙𝜀_𝑘𝑙−𝜀^∗_𝑘𝑙+𝑘_𝑖𝑛𝐸_𝑛+𝐷^∗_𝑖,𝜀^∗_𝑘𝑙=𝜀^𝑚∗_𝑘𝑙+𝜀^𝑒∗_𝑘𝑙 (5)(6)(7)(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 𝜎_𝑖𝑗𝑛_𝑖=𝑇_𝑗. (9)

The given magnetic charge boundary on 𝑆^𝐵 and electric charge boundary on 𝑆^𝐷 are 𝐵_𝑖𝑛_𝑖=𝜎_𝑀,𝐷_𝑖𝑛_𝑖=𝜎_𝐸. (10)𝜎_𝑀 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 𝑢_𝑖=𝑢_𝑖,𝜙=𝜙,𝜑=𝜑. (11)

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: _𝑉𝜎_{𝑖𝑗,𝑖}+𝑓_𝑗𝛿𝑢_𝑗+𝐵_𝑖,𝑖𝛿𝜙+𝐷_𝑖,𝑖−𝑝^𝑒𝛿𝜑d𝑉−_𝑆𝜎𝜎_𝑖𝑗𝑛_𝑖−𝑇_𝑗𝛿𝑢_𝑗d𝑆−_𝑆𝐵𝐵_𝑖𝑛_𝑖−𝜎_𝑀−𝛿𝜙d𝑆_𝑆𝐷𝐷_𝑖𝑛_𝑖−𝜎_𝐸𝛿𝜑d𝑆=0. (12)

Using Gauss theorem and constitutive equations (5) and (6), the desired potential is Π=_𝑉12𝑐_{𝑖𝑗𝑘𝑙}𝜀_𝑖𝑗𝜀_𝑘𝑙−12𝜇₀𝐻_𝑖𝐻_𝑖−𝑒_𝑘𝑖𝑗𝐸_𝑘𝜀_𝑖𝑗−12𝑘_𝑖𝑗𝐸_𝑖𝐸_𝑗−d𝑉_𝑉𝑐_{𝑖𝑗𝑘𝑙}𝜀^∗_𝑖𝑗𝜀_𝑘𝑙+𝜇₀𝑀_𝑖𝐻_𝑖−𝑒_𝑖𝑗𝑘𝜀^∗_𝑗𝑘−𝐷^∗_𝑖−d𝑉_𝑉𝑓_𝑖𝑢_𝑖−𝑝^𝑒𝜑d𝑉−_𝑆𝜎𝑇_𝑖𝑢_𝑖𝑑𝑆−_𝑆𝐵𝜎_𝑀−𝜙d𝑆_𝑆𝐷𝜎_𝐸𝜑d𝑆. (13) 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: 𝐮=𝑢_𝑥𝑢_𝑦𝜑𝜙^𝑇,𝜺=𝜀_𝑥𝜀_𝑦𝛾_𝑥𝑦−𝐸_𝑥−𝐸_𝑦−𝐻_𝑥−𝐻_𝑦^𝑇,𝝈=𝜎_𝑥𝜎_𝑦𝜎_𝑥𝑦𝐷_𝑥𝐷_𝑦𝐵_𝑥𝐵_𝑦^𝑇. (14) The generalized stress and strain are linked by the constitutive relation as 𝝈=𝐃𝜺+𝐂^∗. (15) The matrix 𝐃 and matrix 𝐂^∗ are reduced from three-dimensional constitutive relations, which are given in Appendix A.

In an element, the generalized displacement can be interpolated as 𝐮=𝐍𝐚, (16) where 𝐚 is the generalized nodal displacement and 𝐍 is the shape function as follows: 𝐚=𝐚₁𝐚₂⋯𝐚_𝑚^𝑇,𝐚_𝑖=𝑢_𝑥^(𝑖)𝑢_𝑦^(𝑖)𝜑^(𝑖)𝜙^(𝑖),𝐍=𝐍₁𝐍₂⋯𝐍_𝑚,𝐍_𝑖=𝑁_𝑖0000𝑁_𝑖0000𝑁_𝑖0000𝑁_𝑖. (17) Then following the procedure of FEM, the stiff matrix and force matrix of an element are formed as 𝐊^𝑒=_𝑉12𝐁^𝑇𝐏𝐃𝐁d𝑥d𝑦,_𝑒=_𝑉𝐁^𝑇𝐂^∗d𝑥d𝑦+_𝑉𝐍^𝑇+𝐟d𝑥d𝑦_𝑆𝜎𝐍^𝑇𝐓d𝑆+_𝑆𝐵𝐍^𝑇𝐩^𝑚d𝑆+_𝑆𝐷𝐍^𝑇𝐩^𝑒d𝑆, (18) where 𝐁=𝐁₁𝐁₂⋯𝐁_𝑚,𝐁_𝑖=𝜕𝑁_𝑖0𝜕𝑥𝜕𝑁_𝑖0𝜕𝑦0000𝜕𝑁_𝑖𝜕𝑦𝜕𝑁_𝑖𝜕𝑥0000000𝜕𝑁_𝑖𝜕𝑥𝜕𝑁_𝑖𝜕𝑦0000000𝜕𝑁_𝑖𝜕𝑥𝜕𝑁_𝑖𝜕𝑦^𝑇,𝐟=𝑓_𝑥𝑓_𝑦−𝑝^𝑒0^𝑇,𝐓=𝑇_𝑥𝑇_𝑦00^𝑇,𝐩^𝑒=00𝜎_𝐸0^𝑇,𝐩^𝑚=000𝜎_𝑀^𝑇. (19)

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

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: 𝑀_𝑖^(𝑘)=𝑀_𝑠𝑚_𝑖,𝜀^(𝑘)_𝑖𝑗=32𝜆₁₀₀𝑚²_𝑖−13,𝑚_𝑖=𝑚_𝑗3𝜆₁₁₁𝑚_𝑖𝑚_𝑗,𝑚_𝑖≠𝑚_𝑗, (20) where 𝜆₁₀₀ and 𝜆₁₁₁ 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 𝑀_𝑖=⁹_𝑘=1𝑃_𝑘𝑀_𝑖^(𝑘),𝜀^𝑚∗_𝑖𝑗=⁹_𝑘=1𝑃_𝑘𝜀^(𝑘)_𝑖𝑗, (21) where 𝑃_𝑘 is the volume fraction of the 𝑘th 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] 𝑃_𝑘=𝑁exp−𝐸_𝑘𝜔, (22) where 𝐸_𝑘 and 𝑁 are domain energy and normalization factor, respectively, and 𝜔 is called energy distribution parameter, whose expression is 𝜎𝜔=𝐴+𝐵^eq𝜎₀, (23) where 𝐴, 𝐵, and 𝜎₀ are fitting parameters. 𝜎_eq is the Mises equivalent stress.

The total free energy of each domain is 𝐸=𝐸_𝐾+𝐸_el+𝐸_me+𝐸_𝐻+𝐸_𝜎, (24) 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 𝐾^′=𝐾₁+Δ𝐾₁≈𝐾₁+9/4𝜆²₁₀₀(𝐶₁₁−𝐶₁₂)−9/2𝜆²₁₁₁𝐶₄₄, where 𝐾₁ is the magnetocrystalline anisotropy constant and 𝐶₁₁, 𝐶₁₂, and 𝐶₄₄ are elastic constants, so that 𝐸_𝐾′=𝐸_𝐾+𝐸_el+𝐸_me=𝐾^′𝑚²₁𝑚²₂+𝑚²₂𝑚²₃+𝑚²₃𝑚²₁,𝐸_𝐻=−𝜇₀𝑀_𝑠𝐻_𝑖𝛼_𝑖,𝐸_𝜎3=−2𝜆₁₀₀𝑚²₁𝜎₁₁+𝑚²₂𝜎₂₂+𝑚²₃𝜎₃₃−3𝜆₁₁₁𝑚₁𝑚₂𝜎₁₂+𝑚₂𝑚₃𝜎₂₃+𝑚₁𝑚₃𝜎₃₁. (25)

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) in which the elastic moduli 𝑀^𝐸_{𝑖𝑗𝑘𝑛}, piezoelectric moduli 𝑑_𝑘𝑖𝑗, dielectric moduli 𝑘^𝜎_𝑗𝑘, remnant strain 𝜀^𝑒∗_𝑖𝑗, and electric displacement 𝐷^∗_𝑗 are calculated by averaging corresponding values of each domain: 𝑀^𝐸_{𝑖𝑗𝑘𝑛}=1𝑚^𝑚_𝑙=1𝑀^𝐸(𝑙)_{𝑖𝑗𝑘𝑛},𝑑_𝑖𝑗𝑘=1𝑚^𝑚_𝑙=1𝑑^(𝑙)_𝑖𝑗𝑘,𝑘^𝜎_𝑖𝑗=1𝑚^𝑚_𝑙=1𝑘^𝜎(𝑙)_𝑖𝑗,𝜀^𝑒∗_𝑖𝑗=1𝑚^𝑚_𝑙=1𝜀^∗(𝑙)_𝑖𝑗,𝐷^∗_𝑗=1𝑚^𝑚_𝑙=1𝐷_𝑗^∗(𝑙), (27) 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.

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 state 𝑆 is [30] 𝑔(𝝈,𝐄,𝑆)=−𝜺^{𝑒∗(𝑆)}∶𝝈+𝐃^∗(𝑆)1⋅𝐄+2𝝈∶𝐌^𝐄(𝑆)+1∶𝝈2𝐄⋅𝐤^𝜎(𝑆)⋅𝐄+𝐄⋅𝐝^(𝑆)∶𝝈. (28)

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.(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 Figure 1.

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 Tb_0.3Dy_0.7Fe_1.95 compound. Their material constants are listed in Tables 1 and 2, respectively.

Table 1 

Parameters of PLZT [27].

Table 2 

Parameters of Tb_0.3Dy_0.7Fe_1.95 [2].

5.1. Constitutive Curve

Our program is first validated by the experimental constitutive curves of PLZT and Tb_0.3Dy_0.7Fe_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.

(a)

(b)

(c)

(d)

Figure 2 

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

5.2. The Field Distribution in Three-Layer ME Composite

Figure 3 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 𝐿=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.

Figure 3 

Schematic illustration of the three-layer laminate ME composite.

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.

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 

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 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𝜎_𝑥^(𝑝)d𝑥,attheends. (29) 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) 𝐵=2𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑚2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝜀^∗_𝑥×2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑝𝑡_𝑚𝑡_𝑏𝐺^(𝑏)_𝑥𝑦+16𝑡_𝑝𝐺^(𝑝)_𝑥𝑦+13𝑡_𝑚𝐺^(𝑚)_𝑥𝑦⁻¹, (30) where 𝜀^∗_𝑥 is the magnetostriction. In the extreme case 𝐺^(𝑏)_𝑥𝑦/𝑡_𝑏→∞, we have 𝐵_perfect=2𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑚2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝜀^∗_𝑥×2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑝𝑡_𝑚16𝑡_𝑝𝐺^(𝑝)_𝑥𝑦+13𝑡_𝑚𝐺^(𝑚)_𝑥𝑦⁻¹ (31) which gives rise to the best interfacial bonding condition one can ever achieve. Then the interfacial transferring coefficient is defined as the normalized term 𝐵𝐵=𝐵_perfect=𝑡_𝑝/𝐺^(𝑝)_𝑥𝑦+2𝑡_𝑚/𝐺^(𝑚)_𝑥𝑦6𝑡_𝑏/𝐺^(𝑏)_𝑥𝑦+𝑡_𝑝/𝐺^(𝑝)_𝑥𝑦+2𝑡_𝑚/𝐺^(𝑚)_𝑥𝑦. (32)

Figure 6 

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

In the case of Section 5.2, (32) gives 𝐵=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 𝐵=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) 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 𝜏𝑥𝑦, (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 8 

Influence of interfacial bonding on ME coefficient. In the simulation, interfacial bonding is controlled by the shear module 𝐺𝑏 and thickness 𝑡𝑏 of the bonding layer. Large value of 𝐺𝑏/𝑡𝑏 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 (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.

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.

Figure 10 

Numerical results on the field distributions of a three-layer laminated composite with an interfacial crack in the middle. (a) Shear stress 𝜏𝑥𝑦, (b) normal stress 𝜎𝑥, (c) electric field 𝐸𝑦, (d) magnetization 𝑀𝑥. 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 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.

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

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 𝜏𝑥𝑦, (b) normal stress 𝜎𝑥, (c) electric field 𝐸𝑦, (d) magnetization 𝑀𝑥.

(a)

(b)

Figure 14 

(a) Interfacial transferring coefficient as function of 𝐺𝑏/𝑡𝑏, (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.

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.

Appendices

A. The Elastic Stiff Matrix 𝐷 and Matrix 𝐶^∗ in Planar Problems

In this appendix, we derive the matrix 𝐃 and 𝐂^∗ in (15) for planar strain and planar stress problems from three-dimensional constitutive relations. Considering the problem in 𝑥-𝑦 plane, the extended generalized stress and strain are related as 𝝈^′=𝐃^′𝜺^′−𝜺^′∗+𝝈^′∗=𝐃^′𝜺^′+𝐂^′∗, (A.1) where 𝜺^′=𝜀_𝑥𝜀_𝑦𝛾_𝑥𝑦−𝐸_𝑥−𝐸_𝑦−𝐻_𝑥−𝐻_𝑦𝜀_𝑧^𝑇,𝜺^′∗=𝜀^∗_𝑥𝜀^∗_𝑦𝛾^∗_𝑥𝑦0000𝜀^∗_𝑧^𝑇,𝝈^′=𝜎_𝑥𝜎_𝑦𝜎_𝑥𝑦𝐷_𝑥𝐷_𝑦𝐵_𝑥𝐵_𝑦𝜎_𝑧^𝑇,𝝈^′∗=000𝐷^∗_𝑥𝐷^∗_𝑦𝜇₀𝑀_𝑥𝜇₀𝑀_𝑦0^𝑇,𝐃^′=𝑐₁₁₁₁𝑐₁₁₂₂𝑐₁₁₁₂𝑒₁₁₁𝑒₂₁₁00𝑐₁₁₃₃𝑐₁₁₂₂𝑐₂₂₂₂𝑐₂₂₁₂𝑒₁₂₂𝑒₂₂₂00𝑐₂₂₃₃𝑐₁₁₁₂𝑐₂₂₁₂𝑐₁₂₁₂𝑒₁₂₁𝑒₂₂₁00𝑐₁₂₃₃𝑒₁₁₁𝑒₁₂₂𝑒₁₂₁−𝑘₁₁−𝑘₁₂00𝑒₁₃₃𝑒₂₁₁𝑒₂₂₂𝑒₂₂₁−𝑘₁₂−𝑘₂₂00𝑒₂₃₃00000−𝜇₀00000000−𝜇₀0𝑐₁₁₃₃𝑐₂₂₃₃𝑐₁₂₃₃𝑒₁₃₃𝑒₂₃₃00𝑐₃₃₃₃=𝐃_II(7×7)𝐃_III(7×1)𝐃_III(1×7)𝐷_IIII(1×1),𝐂^′∗=𝝈^′∗−𝐃^′𝜺^′∗=𝐂_I(7×1)𝐶_II(1×1). (A.2)(A.3)(A.4)(A.5)(A.6)(A.7) We rewrite (A.1) as 𝝈=𝐃_II𝜺+𝐃_III𝜀_𝑧+𝐂_I,𝜎_𝑧=𝐃_III𝜺+𝐷_IIII𝜀_𝑧+𝐶_II. (A.8)(A.9) In planar strain problem, 𝜀_𝑧=0, so we have 𝐃=𝐃_II,𝐂^∗=𝐂_I. (A.10) In planar stress problem, 𝜎_𝑧=0, and 𝜀_𝑧 is solved from (A.9) as 𝜀_𝑧1=−𝐷_IIII𝐃_III𝐶𝜺−^II𝐷_IIII. (A.11) Substitute (A.11) into (A.8), we have 𝝈=𝐃_II−1𝐷_IIII𝐃_III𝐃_III𝜺+𝐂_I−𝐶_II𝐷_IIII𝐃_III. (A.12) Therefore 𝐃=𝐃_II−1𝐷_IIII𝐃_III𝐃_III,𝐂^∗=𝐂_I−𝐶_II𝐷_IIII𝐃_III. (A.13)

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: 𝜀_𝑥^(𝑝)=1𝐸_𝑥^(𝑝)𝜎_𝑥^(𝑝),𝜀_𝑥^(𝑚)=1𝐸_𝑥^(𝑚)𝜎_𝑥^(𝑚)+𝜀^∗_𝑥,𝛾^(𝑝)_𝑥𝑦=1𝐺^(𝑝)_𝑥𝑦𝜏^(𝑝)_𝑥𝑦,𝛾^(𝑚)_𝑥𝑦=1𝐺^(𝑚)_𝑥𝑦𝜏^(𝑚)_𝑥𝑦,𝛾^(𝑏)_𝑥𝑦=1𝐺^(𝑏)_𝑥𝑦𝜏^(𝑏)_𝑥𝑦, (B.1)(B.2)(B.3)(B.4)(B.5) where the superscript (𝑝) stands for the piezoelectric layer, (𝑚) stands for the magnetostrictive layer, and (𝑏) stands for the bonding layer. 𝜎_𝑥^(𝑝), 𝜎_𝑥^(𝑚), 𝜏^(𝑝)_𝑥𝑦, 𝜏^(𝑚)_𝑥𝑦, and 𝜏^(𝑏)_𝑥𝑦 are stresses; 𝜀_𝑥^(𝑝), 𝜀_𝑥^(𝑚), 𝛾^(𝑝)_𝑥𝑦, 𝛾^(𝑚)_𝑥𝑦, and 𝛾^(𝑏)_𝑥𝑦 are strains; 𝐸_𝑥^(𝑝) and 𝐸_𝑥^(𝑚) are Young’s moduli; 𝐺^(𝑝)_𝑥𝑦, 𝐺^(𝑚)_𝑥𝑦, and 𝐺^(𝑏)_𝑥𝑦 are shear moduli; 𝜀^∗_𝑥 is the magnetostriction.

The displacements of the piezoelectric layer and the magnetostrictive layer are related by 𝑢^(𝑝)−𝑢^(𝑚)=𝑡_𝑏𝜏𝐺^(𝑏)_𝑥𝑦+16𝑡_𝑝𝜏𝐺^(𝑝)_𝑥𝑦+13𝑡_𝑚𝜏𝐺^(𝑚)_𝑥𝑦=𝑁𝜏, (B.6) where 𝜏 is the interfacial shear stress. In deriving (B.6), we use the shear-lag model [31]. The balance of magnetostrictive layer leads to d𝑡_𝑚𝜎_𝑥^(𝑚)d𝑥+𝜏=0 (B.7) and the overall balance equation is 2𝑡_𝑚𝜎_𝑥^(𝑚)+𝑡_𝑝𝜎_𝑥^(𝑝)=0. (B.8) Equations (B.1), (B.2), and (B.6)–(B.8) and the geometry relation 𝜀_𝑥^(𝑚)=d𝑢^(𝑚)d𝑥,𝜀_𝑥^(𝑝)=d𝑢^(𝑝)d𝑥 (B.9) give rise to the equation d²𝜎_𝑥^(𝑚)d𝑥²−2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝑁𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑝𝑡_𝑚𝜎_𝑥^(𝑚)=𝜀^∗_𝑥1𝑁𝑡_𝑚. (B.10) Under the traction free boundary condition at two ends, namely, 𝐿𝑥=±2,𝜎_𝑥^(𝑚)=0,𝜎_𝑥^(𝑝)=0, (B.11) (B.10) is solved as 𝜎_𝑥^(𝑚)𝐸=−^𝑥(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑝2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝜀^∗_𝑥×√1−sinh𝑀(𝐿/2+𝑥)sinh√−√𝑀𝐿sinh𝑀(𝐿/2−𝑥)sinh√,𝑀𝐿 (B.12) where 𝑀=2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝑁𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑝𝑡_𝑚=2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑝𝑡_𝑚×𝑡_𝑏𝐺^(𝑏)_𝑥𝑦+16𝑡_𝑝𝐺^(𝑝)_𝑥𝑦+13𝑡_𝑚𝐺^(𝑚)_𝑥𝑦⁻¹. (B.13) The stress in the piezoelectric layer is 𝜎_𝑥^(𝑝)=2𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑚2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝜀^∗_𝑥×√1−sinh𝑀(𝐿/2+𝑥)sinh√−√𝑀𝐿sinh𝑀(𝐿/2−𝑥)sinh√𝑀𝐿 (B.14) so d𝜎_𝑥^(𝑝)d𝑥_{𝑥=−𝐿/2}=2√𝑀𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑚2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝜀^∗_𝑥cosh√𝑀𝐿−1sinh√.𝑀𝐿 (B.15) When the length of the composite 𝐿 is sufficiently large cosh√𝑀𝐿−1sinh√𝑀𝐿⟶1 (B.16) the stress transferred by a unit length at the ends is 𝐵=2𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑚2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝜀^∗_𝑥×2𝐸_𝑥^(𝑚)𝑡_𝑚+𝐸_𝑥^(𝑝)𝑡_𝑝𝐸_𝑥^(𝑚)𝐸_𝑥^(𝑝)𝑡_𝑝𝑡_𝑚𝑡_𝑏𝐺^(𝑏)_𝑥𝑦+16𝑡_𝑝𝐺^(𝑝)_𝑥𝑦+13𝑡_𝑚𝐺^(𝑚)_𝑥𝑦⁻¹. (B.17)

Conflict of Interests

The authors do not have any conflict of interest with the content of the paper.

Acknowledgments

The authors acknowledge the support from National Natural Science Foundation of China (Grants nos. 11090334, 51232004, and 10820101048), National Basic Research Program of China (973 Program Grant nos. 2010CB832701), and Tsinghua University Initiative Scientific Research Program (no. 2011Z02173).