# Modeling of Small DC Magnetic Field Response in Trilayer Magnetoelectric Laminate Composites

Authors: B. Zadov; A. Elmalem; E. Paperno; I. Gluzman; A. Nudelman; D. Levron; A. Grosz; S. Lineykin; E. Liverts
Journal: Advances in Condensed Matter Physics
 (2012)

---

## Abstract

We consider a magnetoelectric laminate which comprises two magnetostrictive (Ni) layers and an in-between piezoelectric layer (PZT). Using the finite-element method-based software COMSOL, we numerically calculate the induced voltage between the two faces of the PZT piezoelectric layer, by an external homogeneous small-signal magnetic field threading the three-layer Ni/PZT/Ni laminate structure. A bias magnetic field is simulated as being produced by two permanent magnets, as it is done in real experimental setups. For approaching the real materials’ properties, a measured magnetization curve of the Ni plate is used in the computations. The reported results take into account the finite-size effects of the structure, such as the fringing electric field effect and the demagnetization, as well as the effect of the finite conductivity of the Ni layers on the output voltage. The results of the simulations are compared with the experimental data and with a widely known analytical result for the induced magnetoelectric voltage.


---

## Body

## 1. Introduction

Magnetostrictive-piezoelectric laminates exhibiting magnetoelectric (ME) effect have drawn increasing interest due to their potential for many modern devices, such as sensors, gyrators, and energy harvesters [1–3]. The strong ME effect was recently observed [4] in artificially fabricated multiferroic composites, where the two different-phase materials, that is, the piezoelectric and magnetostrictive single-phase materials, are bonded together. The simplest type of the ME laminate comprises two parallel magnetostrictive and piezoelectric plates well bonded to each other. The analytical expressions for the ME coefficient αME of such a configuration are obtained so far [5–7] under the assumption of homogeneity of electric, magnetic, and elastic field and employing boundary conditions for mechanical stress in the integral meaning at the structure facets. As was argued [8], those expressions for αME are in fact leading asymptotic ones at the lateral-to-transverse dimensions ratios tending to infinity, which result from neglecting the true boundary conditions for all the previously mentioned physical fields. However, in this approach, αME depends only on the fractional thicknesses of the piezoelectric and magnetostrictive layers and not on the lateral dimensions of the laminate, which contradicts the new experimental data [9, 10]. Various numerical methods were proposed to study the ME effect in the multiferroic composite including the finite-element method (FEM) modeling [11–13]; the FEM modeling was also applied to the composite multiferroic device analysis [12].In addition, in most research on the ME effect in magnetostrictive/piezoelectric composites carried out so far, it was assumed that constitutive relations for both piezoelectric and magnetostrictive single phase are linear. In this framework, one of the principal parameters which enter the previousαME is the piezomagnetic coefficient q. Often it is tacitly assumed that q is defined by the derivative of the magnetostrictive strain with respect to the magnetic field of the bare magnetic material in the bias magnetic field and, therefore, can be estimated as the magnetostrictive strain λs at magnetic saturation divided by the corresponding magnetic saturation field. However, this theoretical prediction does not agree well with the values of q retrieved from the measurements of αME; see [10, 14]. This discrepancy is overcome by introducing a correction factor for q [10, 14].The idea behind our present study is that both the magnetostrictive and piezoelectric layers of the laminate are in a stressed state as a result of the bonding. The stressed state induces changes in many physical properties, especially, an appreciable change may occur in the magnetic and magnetostrictive characteristics of the magnetostrictive layer made of a ferromagnetic nickel (Ni) with cubic symmetry. Such changes are well studied for standalone samples of ferromagnetic Ni and Fe subjected to a homogeneous stress [15, 16]. This essential impact of the stress is in general nonlinear in magnetic field and was not taken into consideration for the ME laminates both in the framework of the analytical [1–3, 5–7] and numerical [11–13] approaches. Recently, a model for the magnetostriction nonlinearity both with respect to the stress and magnetization has been proposed for the Terfenol-D magnetostrictive material (see [17] and references therein). Very recently, this model has been utilized [18] in the numerical modeling of Terfenol-D-based ME laminates. Though this model does describe the previously noted effect of stress on the magnetization curve, a nonlinearity of the strain-stress-magnetization constitutive relation with respect to stress seems superfluous and introduces too much fitting constants. At the same time, adding there nonlinear terms of even orders in the magnetization, beyond the well-known second order, is acceptable [15, 16]. In the approach based on the strain-stress-magnetization constitutive relation linear in stress, the discussed effect emerges via the appearance of stress-induced excess magnetic anisotropy energy. Depending on the stress value, this excess energy may cause notable changes of the magnetization curve, which can be accounted for using full micromagnetic computations.In the present study, we overcome this limitation of the models employed thus far and avoid solving the accompanying nonlinear micromagnetic problem. For this purpose, we use an explicitH-B curve [15, 16] of the Ni plates bonded in the ME laminate, which we model in such a manner that the measured and computed strains of the laminate that develop under application of an external magnetic field are close to each other as much as possible.Currently, there is no commercial FEM package specifically designed to model multiferroic systems. On the other hand, general-purpose Multiphysics solvers can solve arbitrary sets of partial differential equations. In this paper, a case study of using COMSOL Multiphysics V4.1 FEM modeling software (COMSOL) to design and optimize a finite-size trilayer ME laminate structures will be described. Issues such as fringing electric field effect and the demagnetization as well as the effect of the finite conductivity of the Ni layers on the output voltage will be covered.Note that the fringing electric field and conductivity effect depend strongly on the ME laminate operation regime in a real circuit, for example, open- or closed-circuit, grounded or nongrounded electrode. Using the linear elastic model embedded into COMSOL, we have developed a magnetostrictive model which enables us to evaluate mechanical strain and stress induced by magnetic field in a magnetic material. A proper combination of all built-in physical models realized in COMSOL allowed us to investigate numerically the trilayer ME structures.This paper is organized as follows. Section2 describes a general physical concept of an ME response of magnetostrictive/piezoelectric laminate, its mostly used operation regime in a real circuit, and defines its principal figure of merit, namely, an ME voltage coefficient. Section 3 contains a detailed description of basic mathematical models, that is, governing equations, constitutive relations, and boundary conditions, used in the numerical simulations. In Section 4, the simulation setup is described in much detail. Section 5 presents the numerical simulations’ results and measurement’s apparatus and results as well as comparison between the numerical simulations and experiment. Finally, in Section 6 the conclusions are drawn.


## 2. Functional Principle of Magnetoelectric Action of a Trilayer-Laminate Structure

Figure1 shows how the ME effect appears in a trilayer magnetostrictive/piezoelectric laminate structure. The applied magnetic field that threads the structure causes an elastic deformation in both magnetostrictive layers which wrap the piezoelectric piece. In the particular case of Ni, this is a compressive strain in the direction of the applied magnetic field and a tensile strain in the directions transverse to the field, as shown in Figure 1. In turn, due to bonding conditions at the two layers’ interfaces, this deformation is transferred to the wrapped piezoelectric layer. The deformation of the piezoelectric layer causes appearance of an electric polarization directed perpendicular to the plane.Figure 1
The magnetoelectric effect in a trilayer magnetostrictive/piezoelectric/magnetostrictive laminate structure.One of the most used techniques to measure the electric response resulting from this polarization is to connect the ME laminate to an infinite-impedance measurement instrument like an ideal voltmeter or an electrometer. In this case, no net current is allowed to flow. Note that one should distinguish two types of magnetostrictive materials employed: highly conductive, for example, Ni, Terfenol-D, and poorly conductive, for example, ferrites, and two connection conditions—one of the electrodes grounded or not grounded. In the present study, we restrict ourselves to a laminate which comprises two (highly conductive) magnetostrictive Ni layers, serving also as the electrodes, and a PZT layer which we assume to be an insulator. Figure2 illustrates schematically an electrostatic diagram corresponding to the previously noted circuit connection with the bottom Ni plate being grounded. By virtue of the induced polarization and nonzero conductivity of the magnetostrictive material, there appears an inhomogeneous charge distribution on those surfaces of both Ni plates, which interface with the PZT layer. As a result of grounding the bottom and electrical isolating the top electrode, the grounded electrode acquires a nonzero net charge, while the isolated electrode will have zero net charge. The surface charge induced on the interfaces will spread over all the surfaces of the Ni plates, as shown in Figure 2. Because PZT is a functional, rather than plain, dielectric, the charges on the Ni/PZT interfaces do not follow the electric field alone but depend also on an elastic state of the PZT layer. At the same time, just the charges on the Ni surfaces in contact with air produce an electric field in the area surrounding the laminate structure. The volume-polarization and surface-charge distributions along with the measurement circuit conditions (see Figure 1) result in an induced voltage across the structure. This voltage VME(H) is just the manifestation of the product ME effect in these specific structures.Figure 2
Schematic electrostatic diagram of the ME laminate corresponding to the circuit connection with a grounded bottom electrode.In sensor applications, it is notVME(H) that is measured, rather it is the change δVME(H) in the response to a modulated small-field signal δH which superimposes the bias magnetic field H. Then, the essential parameter which defines the ac (up to a minute frequency) ME response of the laminates-based circuit is the ME voltage coefficient:(1)αME=δVME(H)δH.
It is worth to note a distinctive feature of the magnetostrictive part of the problem compared to the piezoelectric part. Namely, the magnetostrictive constitutive relation is linear in the strain components but nonlinear in magnetization distribution, which has been taken into account in [19] and is employed in the present study. At the same time, the piezoelectric constitutive relations are overall linear.In general, also the elastic constants of the magnetostrictive materials may depend on the magnetization, but this effect is not as prominent as the demagnetization. The latter effect is due to the fact that the magnetization distribution in finite-size samples induces internal excess magnetic fieldHm. It is the total field H+Hm that enters the constitutive relation for the magnetostrictive materials, the fact which is properly taken into consideration in the computational model presented in this paper, see the following.


## 3. Mathematical Modeling

To accurately analyze the response of the ME laminate structures qualitatively outlined previously, the coupling between the electric, magnetic, and mechanical fields in the magnetostrictive and piezoelectric layers should be taken into consideration. To solve the problem, the coupled fields of the three types mentioned above are to be computed in two steps. Firstly, the source magnetic field driven by the permanent magnets is to be computed, the key variable being the magnetic potential. Secondly, the coupled mechanical and magnetic and electric fields are to be computed in terms of the variables such as the components of the mechanical displacement vector (u), normalized magnetic potential (Ψ), magnetization (M), electric polarization (P), and electric potential (φ).To realize this program, we employed the previously noted COMSOL software. Different modules of COMSOL can be easily combined to simultaneously simulate any coupled physical fields in physical and engineering problems with various geometries, including the laminate structures under consideration. In our problem, to obtain the static deformations, three different modules were used, as shown in Figure3. Contrary to the well-known analytical model, see review [1–3, 5–7], which is in fact based on three assumptions: (i) plane stress, (ii) overall homogeneous physical fields, and (iii) zero average stress at the laminate’s facets, and the modification of this model by taking into account measurement circuit condition [8, 12] and introducing magnetometric demagnetizing factor [10], the approach using COMSOL allows, with proper combination of its modules and models, a solution of the problem numerically with fully correct boundary conditions both at the laminate’s interfaces and facets.Figure 3
Schematic description of the COMSOL Multiphysics V4.1 modules and models used in the simulation.Since there is no magnetoelectric module in COMSOL, we have used its built-in magnetic field (MF) and piezoelectric devices (PZD) modules to simulate ME effect (Figure3).


### 3.1. Modeling the Magnetic Fields

The MF module provides an interface for modeling the magnetic flux density and magnetic field when the magnetic components are part of a larger system that have an impact on other components of the system through other physical fields, such as in the considered case when Ni layers act on PZT layers through elastic displacement. MF module allows us to take into account the magnetic properties of the materials and generate external magnetic fields in the model.The MF module includes solving the magnetostatics equations,(2)∇×H=0,∇⋅B=0,H=-∇ψ,B=μ0(H+M),
where B is the magnetic flux density and μ0 the vacuum magnetic permeability. MF module allows us to model the permanent magnets which provide the ME structure with surrounding source magnetic field, by specifying the constant magnetization of the magnets (M0); for example, for permanent magnets, we used M0 in the range 300–4300 kA/m and fixed magnetic permeability μ=1.05. The model uses (2) of the MF module and physical interface from the AC/DC COMSOL module. The magnetic field is symmetric with respect to the xy,xz, and xz planes; these planes, therefore, may serve as exterior boundary to the geometry. If the air box is sufficiently large, the boundary condition used on its remaining exterior boundaries has little influence on the field in the vicinity of the magnets and between them. Although an infinite domain would give the best results, this model uses the magnetic insulation condition on a finite domain for convenience. All the exterior boundaries are magnetically insulated.In the case of magnetostrictive Ni, however, the magnetic permeability, and the magnetostrictive strains as well, becomes nonlinearly dependent on both the magnetic flux density and the mechanical stresses and strains that arise in the laminate. The simulation of magnetostatics in the magnetostrictive layers has been performed with the use of an experimentally measured nonlinear magnetization versus magnetic field dependenceM(H) using a 99.95% Ni foil supplied by Alfa Aesar. The obtained M(H) is in a good qualitative agreement with well-known experimental data [15]. We used this magnetostrictive material in our laminates for experimental validation of our modeling approach as well. As is customarily to COMSOL, M(H) is plugged into simulations via the so-called H-B curve [15] given by(3)H+M(H)=1μ0B.It appears that the magnetostrictive strains obtained numerically for the ME laminate by plugging theH-B curve, (3), for a standalone Ni foil deviate notably from those strains measured for the same foil bonded in the 3-layer structure described previously. Therefore, we used an H-B curve, reconstructed from εxx versus H measurements. Such a modified H-B curve used in our simulation is shown in Figure 4.Figure 4
TheH-B curve for a Ni layer bonded with PZT modified compared to bare Ni layer as described in the text.Both the permanent magnets and the ME structure are surrounded by air for whichμr=1. In the course of these simulations, COMSOL observes the well-known magnetostatics boundary conditions, that is, the continuity of normal component B and tangential components of H. In addition, magnetic insulation condition (Ψ=0) was applied to an artificial outer air box which wraps the entire magnets+ structure arrangement.


### 3.2. Building a Magnetostriction Model in COMSOL

In the present paper, the magnetostriction strains are modeled by a modification of the expression due to Becker and Döring [15, 20] for polycrystalline cubic ferromagnet,(4)εi(ms)=32λs(αi2¯-13),
where εi(ms)=εik(ms)nk is the ith component of the relative magnetostrictive deformation along the direction n̂, εik(ms) being the full magnetostrictive strain tensor, λs is the polycrystalline magnetostriction constant which is expressed through the magnetostriction constants of the corresponding crystal along two principal cubic axes [15, 20], αi=mi/Ms, mi and Ms being the ith component of local magnetization in a grain/domain and the saturation magnetization, respectively; the bar means averaging over grains/domains. The expression for εi(ms) vanishes in the absence of a magnetic field. When a magnetic field is applied and magnetization process starts, it is the deformation along the field (ε∥(ms)) and perpendicular to the field (ε⊥(ms)) that is unique for polycrystalline material. Furthermore, ε∥(ms)=λs and ε⊥(ms)=-λs/2 at magnetic saturation. Knowing the magnetization curve M(H) which is in fact M(H)=M∥(H)=m∥¯, COMSOL allows us to calculate M⊥(H)=m⊥¯ but, unfortunately, does not allow the computation of neither α∥2¯, nor α⊥2¯. Therefore, to build a magnetostriction model in the COMSOL environment, we use the following approximation [19]:(5)ε∥(ms)=λs(MMs)2,ε⊥(ms)=-λs2(MMs)2
which strictly preserves the two above-mentioned limiting cases of the nonmagnetized and magnetically saturated sample. In addition, this model turns out to be pretty much physically reasonable regarding nonlinear dependence of the magnetostriction on the magnetization, at least near saturation. The near-saturation range is especially important for our modeling because it is this range which produces maximum ME response. Similar approximation has recently been considered in [19]. The following parameters, validated by magnetometric and strain gauge measurements, were adopted in our simulation: Ms=480kA/m and λs=-3.4×10-5.


### 3.3. Elastic Model for Magnetostrictive Material

To simulate the elastic behavior of the magnetostrictive layers, we used a linear elastic model which is a part of PZD Module of COMSOL. This model is described by the well-known elastic constitutive relations:(6)ε=12[(∇u)T+∇u],σ-σ0=CE⋅(ε-ε0),∇σ=0,
where ε is the strain tensor, ε0 is an initial strain one, σ is the stress tensor, σ0 is an initial stress, and CE is the elasticity matrix. The elements of the latter matrix are calculated for the Ni layers using Young’s modulus Y=219 GPa and Poison’s ratio ν=0.3 adopted from [15]. Due to our approach, see also [19], the nonlinear magnetostriction is implemented through introduction of prestrain which is equal to the magnetostrictive strain, that is, ε0,ik=εik(ms), the expression for which is discussed previously, and σ0=0. The coupled governing equations (6) for the laminate are solved with the traction-free boundary conditions (σiknk=0) at outer surfaces and the perfect bonding at all the interfaces of the structure (continuity of u and σiknk). This model has been validated on a standalone Ni layer by simulation and measurements of the longitudinal and transverse strains induced by the magnetization process. There is a reasonable agreement between these measurements and the simulation.


### 3.4. Piezoelectric Material Modeling

For the PZT, we employ the linear model, for which the established constitutive relations, written in the strain-charge form, are as follows [21]:(7)ε-ε0=SE⋅(σ-σ0)+dT⋅E,D=Pr+d⋅(σ-σ0)+κσ⋅E,
where E is the electric field, D is the electric displacement, Pr is the remanent polarization (nonzero for a ferroelectric), SE is the compliance coefficients’ matrix at given electric field, d is the piezoelectric coefficients’ matrix, and κσ is the dielectric permittivity matrix at given stress. Here, the elasticity-theory quantities are denoted as in (6), but now they relate just to the PZT under consideration. The elements of the compliance coefficients’ matrix SE conform to the assumption that the polarized piezoelectric phase has the 6 mm symmetry.To build the full model of the PZT layer, we used the PZD module of COMSOL. The PZD Multiphysics interface combines solid mechanics and electrostatics (see Figure3) for modeling of piezoelectric devices, where all or some of the domains contain a piezoelectric material. The piezoelectric coupling can be in strain-charge form of (7). It allows one to take into account coupled piezoelectric, elastic, and electric properties of the piezoelectric layer in the considered structure. The PZD Multiphysics interface provides the following equations for modeling a piezoelectric layer, solving for the displacements (u) and the electric potential (φ):(8)∇⋅D=0,E=-∇φ,∇σ=0.
In the simulation, we used the parameters of a PZT-5H from COMSOL material library, which are listed in Appendix A.To obtain the electric potential (φ) in our model, we assume that we have metal electrodes (Ni plates) and a perfect insulator (PZT layer). Thus, the PZD interface solves the electrostatics within (8) with the proper account of charge balance. Solving the electrostatic part of the problem, the PZD interface observes all the electrostatic boundary conditions, that is, matching the normal component of D with the surface charge density on the electrodes, that is, Dn=ρs, and the continuity of the tangential components of E at the structure’s interfaces. In order to facilitate an adequate comparison between the simulations and measurements’ results, it is needed to further specify the boundary conditions to approach a realistic circuitry. As it was discussed in Section 2, the upper Ni plate is actually in contact with a piezoelectric medium but is electrically isolated from the external circuit. In this case, no net current is allowed to flow at the interface between this so-called electrically floating electrode (EFE) and the embedding media (PZT and air), so the zero net charge is held on this electrode. Just in the case of an EFE regime for this electrode, the constant potential φ on it, which is in fact the sought VME(H), is fully determined from the D and E distributions. It is instructive to note that COMSOL includes this EFE option for the electrodes conditions, which we effectively used in this study. In addition, to make the electrostatic problem well posed, the electric grounding (φ=0) is applied to the second electrode (bottom Ni layer) and to an artificial outer air box which wraps the entire structure, as well.As far as the elasticity part of the problem is concerned, the PZD module takes into consideration solid mechanics material model and the continuity of theu field and the normal components of the tensor σ field at the PZT/Ni interfaces, the traction-free condition on the facets of the PZT layer.


## 3.1. Modeling the Magnetic Fields

The MF module provides an interface for modeling the magnetic flux density and magnetic field when the magnetic components are part of a larger system that have an impact on other components of the system through other physical fields, such as in the considered case when Ni layers act on PZT layers through elastic displacement. MF module allows us to take into account the magnetic properties of the materials and generate external magnetic fields in the model.The MF module includes solving the magnetostatics equations,(2)∇×H=0,∇⋅B=0,H=-∇ψ,B=μ0(H+M),
where B is the magnetic flux density and μ0 the vacuum magnetic permeability. MF module allows us to model the permanent magnets which provide the ME structure with surrounding source magnetic field, by specifying the constant magnetization of the magnets (M0); for example, for permanent magnets, we used M0 in the range 300–4300 kA/m and fixed magnetic permeability μ=1.05. The model uses (2) of the MF module and physical interface from the AC/DC COMSOL module. The magnetic field is symmetric with respect to the xy,xz, and xz planes; these planes, therefore, may serve as exterior boundary to the geometry. If the air box is sufficiently large, the boundary condition used on its remaining exterior boundaries has little influence on the field in the vicinity of the magnets and between them. Although an infinite domain would give the best results, this model uses the magnetic insulation condition on a finite domain for convenience. All the exterior boundaries are magnetically insulated.In the case of magnetostrictive Ni, however, the magnetic permeability, and the magnetostrictive strains as well, becomes nonlinearly dependent on both the magnetic flux density and the mechanical stresses and strains that arise in the laminate. The simulation of magnetostatics in the magnetostrictive layers has been performed with the use of an experimentally measured nonlinear magnetization versus magnetic field dependenceM(H) using a 99.95% Ni foil supplied by Alfa Aesar. The obtained M(H) is in a good qualitative agreement with well-known experimental data [15]. We used this magnetostrictive material in our laminates for experimental validation of our modeling approach as well. As is customarily to COMSOL, M(H) is plugged into simulations via the so-called H-B curve [15] given by(3)H+M(H)=1μ0B.It appears that the magnetostrictive strains obtained numerically for the ME laminate by plugging theH-B curve, (3), for a standalone Ni foil deviate notably from those strains measured for the same foil bonded in the 3-layer structure described previously. Therefore, we used an H-B curve, reconstructed from εxx versus H measurements. Such a modified H-B curve used in our simulation is shown in Figure 4.Figure 4
TheH-B curve for a Ni layer bonded with PZT modified compared to bare Ni layer as described in the text.Both the permanent magnets and the ME structure are surrounded by air for whichμr=1. In the course of these simulations, COMSOL observes the well-known magnetostatics boundary conditions, that is, the continuity of normal component B and tangential components of H. In addition, magnetic insulation condition (Ψ=0) was applied to an artificial outer air box which wraps the entire magnets+ structure arrangement.


## 3.2. Building a Magnetostriction Model in COMSOL

In the present paper, the magnetostriction strains are modeled by a modification of the expression due to Becker and Döring [15, 20] for polycrystalline cubic ferromagnet,(4)εi(ms)=32λs(αi2¯-13),
where εi(ms)=εik(ms)nk is the ith component of the relative magnetostrictive deformation along the direction n̂, εik(ms) being the full magnetostrictive strain tensor, λs is the polycrystalline magnetostriction constant which is expressed through the magnetostriction constants of the corresponding crystal along two principal cubic axes [15, 20], αi=mi/Ms, mi and Ms being the ith component of local magnetization in a grain/domain and the saturation magnetization, respectively; the bar means averaging over grains/domains. The expression for εi(ms) vanishes in the absence of a magnetic field. When a magnetic field is applied and magnetization process starts, it is the deformation along the field (ε∥(ms)) and perpendicular to the field (ε⊥(ms)) that is unique for polycrystalline material. Furthermore, ε∥(ms)=λs and ε⊥(ms)=-λs/2 at magnetic saturation. Knowing the magnetization curve M(H) which is in fact M(H)=M∥(H)=m∥¯, COMSOL allows us to calculate M⊥(H)=m⊥¯ but, unfortunately, does not allow the computation of neither α∥2¯, nor α⊥2¯. Therefore, to build a magnetostriction model in the COMSOL environment, we use the following approximation [19]:(5)ε∥(ms)=λs(MMs)2,ε⊥(ms)=-λs2(MMs)2
which strictly preserves the two above-mentioned limiting cases of the nonmagnetized and magnetically saturated sample. In addition, this model turns out to be pretty much physically reasonable regarding nonlinear dependence of the magnetostriction on the magnetization, at least near saturation. The near-saturation range is especially important for our modeling because it is this range which produces maximum ME response. Similar approximation has recently been considered in [19]. The following parameters, validated by magnetometric and strain gauge measurements, were adopted in our simulation: Ms=480kA/m and λs=-3.4×10-5.


## 3.3. Elastic Model for Magnetostrictive Material

To simulate the elastic behavior of the magnetostrictive layers, we used a linear elastic model which is a part of PZD Module of COMSOL. This model is described by the well-known elastic constitutive relations:(6)ε=12[(∇u)T+∇u],σ-σ0=CE⋅(ε-ε0),∇σ=0,
where ε is the strain tensor, ε0 is an initial strain one, σ is the stress tensor, σ0 is an initial stress, and CE is the elasticity matrix. The elements of the latter matrix are calculated for the Ni layers using Young’s modulus Y=219 GPa and Poison’s ratio ν=0.3 adopted from [15]. Due to our approach, see also [19], the nonlinear magnetostriction is implemented through introduction of prestrain which is equal to the magnetostrictive strain, that is, ε0,ik=εik(ms), the expression for which is discussed previously, and σ0=0. The coupled governing equations (6) for the laminate are solved with the traction-free boundary conditions (σiknk=0) at outer surfaces and the perfect bonding at all the interfaces of the structure (continuity of u and σiknk). This model has been validated on a standalone Ni layer by simulation and measurements of the longitudinal and transverse strains induced by the magnetization process. There is a reasonable agreement between these measurements and the simulation.


## 3.4. Piezoelectric Material Modeling

For the PZT, we employ the linear model, for which the established constitutive relations, written in the strain-charge form, are as follows [21]:(7)ε-ε0=SE⋅(σ-σ0)+dT⋅E,D=Pr+d⋅(σ-σ0)+κσ⋅E,
where E is the electric field, D is the electric displacement, Pr is the remanent polarization (nonzero for a ferroelectric), SE is the compliance coefficients’ matrix at given electric field, d is the piezoelectric coefficients’ matrix, and κσ is the dielectric permittivity matrix at given stress. Here, the elasticity-theory quantities are denoted as in (6), but now they relate just to the PZT under consideration. The elements of the compliance coefficients’ matrix SE conform to the assumption that the polarized piezoelectric phase has the 6 mm symmetry.To build the full model of the PZT layer, we used the PZD module of COMSOL. The PZD Multiphysics interface combines solid mechanics and electrostatics (see Figure3) for modeling of piezoelectric devices, where all or some of the domains contain a piezoelectric material. The piezoelectric coupling can be in strain-charge form of (7). It allows one to take into account coupled piezoelectric, elastic, and electric properties of the piezoelectric layer in the considered structure. The PZD Multiphysics interface provides the following equations for modeling a piezoelectric layer, solving for the displacements (u) and the electric potential (φ):(8)∇⋅D=0,E=-∇φ,∇σ=0.
In the simulation, we used the parameters of a PZT-5H from COMSOL material library, which are listed in Appendix A.To obtain the electric potential (φ) in our model, we assume that we have metal electrodes (Ni plates) and a perfect insulator (PZT layer). Thus, the PZD interface solves the electrostatics within (8) with the proper account of charge balance. Solving the electrostatic part of the problem, the PZD interface observes all the electrostatic boundary conditions, that is, matching the normal component of D with the surface charge density on the electrodes, that is, Dn=ρs, and the continuity of the tangential components of E at the structure’s interfaces. In order to facilitate an adequate comparison between the simulations and measurements’ results, it is needed to further specify the boundary conditions to approach a realistic circuitry. As it was discussed in Section 2, the upper Ni plate is actually in contact with a piezoelectric medium but is electrically isolated from the external circuit. In this case, no net current is allowed to flow at the interface between this so-called electrically floating electrode (EFE) and the embedding media (PZT and air), so the zero net charge is held on this electrode. Just in the case of an EFE regime for this electrode, the constant potential φ on it, which is in fact the sought VME(H), is fully determined from the D and E distributions. It is instructive to note that COMSOL includes this EFE option for the electrodes conditions, which we effectively used in this study. In addition, to make the electrostatic problem well posed, the electric grounding (φ=0) is applied to the second electrode (bottom Ni layer) and to an artificial outer air box which wraps the entire structure, as well.As far as the elasticity part of the problem is concerned, the PZD module takes into consideration solid mechanics material model and the continuity of theu field and the normal components of the tensor σ field at the PZT/Ni interfaces, the traction-free condition on the facets of the PZT layer.


## 4. Realization of the Mathematical Model Using COMSOL

To realize the mathematics described in Section3 and maximally approach a real ME device, we performed numerical simulations by COMSOL and computed the ME response for the Ni/PZT-5H/Ni structure with the following specifications: L=16 mm, W=2, 6, and 16 mm tp=tm=0.5, 1, and 2 mm (for the notations, see Figure 1). The data for the PZT-5H material (similar to the data for PZT-854 [22] used in the fabrication) was taken from the COMSOL library [23]. The results of simulations and measurements are presented in Figures 8–22. All the presented simulation results refer to W=6 mm, tp=tm=1 mm, unless otherwise mentioned. The following computational structures for the simulations purposes were built within the COMSOL environment.


### 4.1. Geometrical Blocks

All geometrical blocks are shown in Figure5. They are the following.(a)
Air block represents the finite space area, where free-space electric and magnetic fields exist; it provides a natural environment to the model. The introduction of the air block allows COMSOL to solve finite-element problem spending much shorter CPU time than needed for time-consuming solution of the problem in infinite space.(b)
PZT block represents the piezoelectric layer that is located in the middle of the air block.(c)
Ni blocks represent magnetostrictive layers. The blocks are placed symmetrically above and below PZT block, which creates laminate-like form of the ME structure. At this stage of simulations, an ideal contact between Ni and PZT blocks is assumed.(d)
Permanent magnet block represents the permanent magnets that are symmetrically located at both sides of the ME laminate.Figure 5
The COMSOL simulations model blocks.The permanent magnets (neodymium magnets also known as NdFeB) are ones of cylindric form. These magnets are the sources of the external magnetic field acting on the structure. This field drives the ME effect in the laminate.


### 4.2. Boundary Conditions

In our model, we have hierarchically defined all blocks as ones to be solved within an ac/dc-magnetic fields’ solver. On the next step, we have defined the laminate blocks as ones to be solved by structural mechanics solver. In this solver, Ni has been defined as a linear elastic material with initial strain resulting from magnetostriction, as explained in Section3.3, whereas PZT has been defined as a linear piezoelectric material. To simulate electrostatic behavior of the model in addition to mechanical and magnetic behavior, we have defined the Ni blocks, air block, and permanent magnets block as electric materials within the PZD module.Magnetic insulation and electric grounding conditions were applied to the air block boundaries. The bottom Ni box was also electrically grounded (φ=0). In addition, symmetry relative to the geometrical center of the structure was effectively used.


### 4.3. Meshing Structure

One of the most important things of the FEM simulations is the mesh distribution. Overall finer meshing allows one to achieve optimal performance regarding accuracy and computation time. The latter is proportional to the mesh resolution and computation resources. The performance can be improved by selective meshing. The main idea of the selective meshing is to make mesh distribution more resolved in more critical areas of the computed structure such as edges and boundary regions and less resolved in less important areas as shown in Figure6.Mesh distribution model: (a) general meshing triangulation of the whole domain embedding the ME structure, (b) visualization of the meshing density in the different sub-domains, additional air blocks are added to the critical domains.
(a)(b)In our model, we defined five levels of mesh resolution in the order of decreasing fineness: the ME structure, permanent magnets, domain adjacent to the ME structure, domain adjacent to the permanent magnets, and outer air block.


## 4.1. Geometrical Blocks

All geometrical blocks are shown in Figure5. They are the following.(a)
Air block represents the finite space area, where free-space electric and magnetic fields exist; it provides a natural environment to the model. The introduction of the air block allows COMSOL to solve finite-element problem spending much shorter CPU time than needed for time-consuming solution of the problem in infinite space.(b)
PZT block represents the piezoelectric layer that is located in the middle of the air block.(c)
Ni blocks represent magnetostrictive layers. The blocks are placed symmetrically above and below PZT block, which creates laminate-like form of the ME structure. At this stage of simulations, an ideal contact between Ni and PZT blocks is assumed.(d)
Permanent magnet block represents the permanent magnets that are symmetrically located at both sides of the ME laminate.Figure 5
The COMSOL simulations model blocks.The permanent magnets (neodymium magnets also known as NdFeB) are ones of cylindric form. These magnets are the sources of the external magnetic field acting on the structure. This field drives the ME effect in the laminate.


## 4.2. Boundary Conditions

In our model, we have hierarchically defined all blocks as ones to be solved within an ac/dc-magnetic fields’ solver. On the next step, we have defined the laminate blocks as ones to be solved by structural mechanics solver. In this solver, Ni has been defined as a linear elastic material with initial strain resulting from magnetostriction, as explained in Section3.3, whereas PZT has been defined as a linear piezoelectric material. To simulate electrostatic behavior of the model in addition to mechanical and magnetic behavior, we have defined the Ni blocks, air block, and permanent magnets block as electric materials within the PZD module.Magnetic insulation and electric grounding conditions were applied to the air block boundaries. The bottom Ni box was also electrically grounded (φ=0). In addition, symmetry relative to the geometrical center of the structure was effectively used.


## 4.3. Meshing Structure

One of the most important things of the FEM simulations is the mesh distribution. Overall finer meshing allows one to achieve optimal performance regarding accuracy and computation time. The latter is proportional to the mesh resolution and computation resources. The performance can be improved by selective meshing. The main idea of the selective meshing is to make mesh distribution more resolved in more critical areas of the computed structure such as edges and boundary regions and less resolved in less important areas as shown in Figure6.Mesh distribution model: (a) general meshing triangulation of the whole domain embedding the ME structure, (b) visualization of the meshing density in the different sub-domains, additional air blocks are added to the critical domains.
(a)(b)In our model, we defined five levels of mesh resolution in the order of decreasing fineness: the ME structure, permanent magnets, domain adjacent to the ME structure, domain adjacent to the permanent magnets, and outer air block.


## 5. Simulation Results and Their Validation

### 5.1. Simulation Results

Figure7 presents the computed magnetic flux density distribution and magnetic flux streamlines. Figure 8 presents the profiles of the x-component of H on the centerlines in middle plane inside both the Ni upper layer and PZT layer.Figure 7
Magnetic flux density streamlines over magnetic flux map.Figure 8
Lateral variation of the magnetic field componentHx on the centerline in the middle planes of the Ni layers (the magenta line) and PZT layer (black line); H(ext)=54.7 kA/m.Figure 9
Electric field streamlines over the electric potential map in the middlex-z plane at H(ext)=54.7 kA/m.Variation of the electric potential withz-coordinate (x=y=0) and L=16 mm. (a): t=2 mm, w=2 mm; (b): t=2 mm, w=6 mm; (c): t=2 mm, w=16 mm; (d): t=1 mm, w=6 mm; (e): t=3 mm, w=6 mm; (f): t=4 mm, w=6 mm.
(a)(b)(c)(d)(e)(f)Variation of the electric field componentEz with (a) x-coordinate and (b) y-coordinate.
(a)(b)Variation of the elastic displacementux (a), uy (b) for various magnetic fields, as shown in the legend, along the PZT layer centerline of the middle plane.
(a)(b)Figure 13
Strainεxx along piezoelectric and magnetostrictive layers: the middle planes and centerline at H(ext)=42 kA/m.Figure 14
3D strainεxx map over displacement field obtained at H=3.95 kA/m.Stresses in the piezoelectric and magnetostrictive layers, middle planes and centerline: (a)σxx versus the length coordinate x; (b) σyy versus the width coordinate y. Curves # 1–4 correspond to the Ni layer and # 5–8 to the PZT layer; the increasing number of the curves correspond to the four increasing magnetic field values: 10.2, 22.9, 35.7, and 54.7 kA/m.
(a)(b)Stressσxx map in the middle plane (x-z) at (a) H=10.2 kA/m and (b) H=54.7 kA/m.
(a)(b)Stressσyy map in the middle plane y-z, L=16 mm, W=6 mm, (a, b): t=1 mm (c, d): t=2 mm; (e, f): t=4 mm.
(a)(b)(c)(d)(e)(f)Figure 18
Stressσzz through the piezoelectric and magnetostrictive layers versus the z-coordinate (x=y=0).Figure 19
SimulatedαME as a function of H.Figure 20
Schematics of the setup for measuring the ME response and elastic deformation.Figure 21
Measurement setup.Figure 22
The measured magnetostrictive strainsεxx and εyy for the Ni layer and Ni/PZT/Ni laminate as a function of H.The coupled ME simulations have been carried out at various values of external magnetic fieldH(ext) meaning the magnetic field produced by the permanent magnets in the absence of the laminate. In what follows, the values of Hx(ext) in the center of the laminate used in the simulations are referred to as simply “magnetic field H.” To avoid misunderstanding, the field H noted in connection with Figures 7 and 8 is induced by H(ext), the value of which is noted in the caption.Examples of the electrostatics part of our simulations at values ofH(ext) between 10.2 and 54.7 kA/m are presented in Figures 9–11. Figure 9 shows the E streamlines over the φ map on the central x-z plane. Figure 10 shows the φ profiles on the z-axis (the middle of the laminate, x=y=0) for different widths and total thickness of ME laminates at different magnetic fields. Figures 11(a) and 11(b) show the variation of Ez as a function of x and y, respectively, on the centerlines of the PZT layer.Figures12(a) and 12(b) show the ux and uy displacement components (along the length and width of the laminate structure), correspondingly, taken in the center of the piezoelectric layer for various values of H as shown in the figure. Figure 13 shows the strain component εxx for both the piezoelectric and magnetostrictive layers; εxx is calculated as the numerical derivative of ux with respect to x, see Figure 12(a). Figure 14 shows the 3D distribution map of εxx and the volume displacement field over it.Figures15(a) and 15(b) show the stress components σxx and σyy, respectively, in both the PZT layer and the top Ni layer, whereas Figures 16 and 17 show the stress σxx map (middle plane x-z) and the stress σyy map (middle plane y-z) correspondingly for several total thicknesses of ME laminates at various magnetic fields. Figure 18 presents the stress components σxx over the laminate.To estimate the modulation of the computed physical parameters, such as net strains and the voltage between the Ni layers, with respect to the magnetic field using the developed COMSOL approach, we swept parametrically the values ofM0 for the permanent magnets in the range 300–4300 kA/m. This procedure resulted in a parametric sweep of the magnetic field in the range 3.5–55 kA/m. Then, to obtain the ME voltage coefficient αME, as prescribed by (1), the output voltage was numerically differentiated by H. We especially focus here on the magnetostriction effect and the ME voltage, because we determined these parameters also experimentally. The quantity of most importance and final target of our simulations is, of course, αME which determines the ME response of the laminate and is shown in Figure 19.


### 5.2. Experimental Setup and Measurements

Our raw materials were polycrystalline ferromagnetic nickel (Ni) and piezoelectric ceramic lead zirconate titanate (PZT); 99.95% purity Ni and PZT-854 plates were obtained from Alfa Aesar and APC, respectively. The relevant data sheets of the materials are displayed in Table1. The Ni/PZT and Ni/PZT/Ni laminates’ samples were fabricated by bonding using Hernon 315 glue. For independent extraction of εxx and εyy, the strain measurements on the Ni layer and Ni/PZT bilayer were carried out using two-axis SR-4 strain gauges by Vishay Micro-Measurements. For the measurement of αME, we employed the setup developed in our lab by using a mechanical apparatus enabling accurate positioning of permanent magnets symmetrically located along the laminate axis, producing stable, variable, and reproducible magnetic bias fields (Hdc) at the sample location. An ac magnetic field with constant amplitude of hac=79,6 A/m at varied frequency was applied by a solenoid and buffer configuration amplifier AD549 with 0.11 fA current noise density. Data was collected with PXI measurement system driven by a LabVIEW program as shown in Figure 20. Figure 21 displays real view of that part of the setup, which corresponds to our computational environment (see Figure 5). Figure 22 shows the magnetostrictive εxx and εyy versus H measured on the standalone Ni (Alfa Aesar 99.95% purity) plate and Ni/PZT/Ni laminate. Note that the laminate’s |εxx | and |εyy| are overall smaller than the corresponding values for the Ni alone. The measured and simulated strains εxx and εyy as a function of H are shown in Figure 23. The measured αME versus H curve is shown in Figure 24.Table 1
Material parameters (compliance coefficientSE, piezoelectric charge constant d, and relative dielectric constant) for PZT [22] and Ni [15].
Properties/materials1110−12 (m2/N)s1210−12 (m2/N)d3110−12 (C/N)κσ/κ0Ni4.57− 1.37——PZT 85416.9−5.08−2403400Figure 23
The measured magnetostrictive strainsεxx and εyy as a function of H: the simulation (the green and blue lines) versus the measurements (the red and cyan markers).Figure 24
The measuredαME as a function of H.


### 5.3. Discussion of the Simulations Results

From Figures7 and 8, one can clearly see the demagnetization effect locally, rather than on average as in the treatment of our previous paper [10]. It is seen that, as is expected [15], the demagnetization effect is pretty much inhomogeneous and strongest near the edges. Figure 8 clearly demonstrates obeying the magnetostatics boundary conditions, namely, the continuity of the normal component of B and tangential components of H, throughout the simulations.It can be seen from Figures9 and 11 that the electric field is inhomogeneous across the structure, in particular one can see essential fringing field near the edges of the laminate, which is clear manifestation of the finite-size effects, whereas the field is nearly homogeneous around the center of the laminate. Figure 10 demonstrates that the COMSOL simulation obeys the grounding condition and confirms a charge accumulation on the electrodes’ surfaces (the φ curve has sharp kinks). Note that these kinks are determined by both (electric and elastic) parts of D. Also note that there exists an electric field in the laminate’s surroundings, though it is essentially smaller than the field inside. The external electric field appears to be unfavorable from the energetic viewpoint as it reduces the ME response defined by the internal electric field. The magnitude of that surrounding electric field depends on the ME laminate’s lateral dimensions and becomes very small if the total thickness is negligible relative to the lateral sizes of the ME laminate, see Figure 10(c). Thus, it can be argued that the boundary condition Dn=0 adopted in the homogeneous fields theory [5–7] does correspond to the approximation of an infinite plate, as it was claimed. Moreover, the simulations obtained in the present study show that the magnitude of the ME response predicted previously by the homogeneous fields theory [5–7] seems rather a lower than upper limit of the observed ME effect, in contrast with the conclusion in [12]. Further estimation, using Figures 9 and 10, of the electric field flux through any surface that surrounds the laminate, but inside the air box, shows existence of a nonzero net charge. It can be understood that this net charge appears on the bottom (grounded) electrode due to stresses induced by the magnetic field inside the PZT layer.Figures12(a), 12(b) and Figure 13 show nearly constant slopes of the components of u with respect to the coordinates deep inside the structure. The ux and uy profiles clearly show a contraction of the laminate structure along the applied magnetic field (the x-axis) and an elongation in the orthogonal direction (the y-axis).While foruy such a behavior continues to the ends of the laminate, ux shows deviation from the linearity near the edges. Figures 13, 14, and 15(a) show that the variation of εxx and σxx is rather slow inside the layers, and it is much more rapid approaching the edges. This behavior can be also detected for the polarization induced in the PZT layer which is of interest in this problem. So the complete electroelastic state of the PZT layer can be represented as a sum of two electroelastic states: the inner electroelastic state and that boundary (transition) layers. The inner electroelastic state varies relatively slowly along the coordinate lines of the middle plane. On the contrary, the boundary layers damp down quickly in the directions perpendicular to the edges and are described by truly 3D electroelasticity equations with appropriate boundary conditions. The FEM simulations in this study properly take account of the “edge effects” and show that the correction due to the boundary layers may introduce notable changes into the description of the internal electroelastic state of the ME laminate. Figures 15–17 are a true hallmark of the finite-size effects. Comparing Figures 15(a) and 15(b) shows that σxx≫σyy, which is an interplay of both the geometry relation L/W>2 and the difference in the magnetostrictive strains ε∥(ms)>ε⊥(ms). While in the PZT layer the curves describing both stresses are concave at all H used, the σyy curves in the Ni layer evolve from being convex in the region of the edges (y=±3 mm) to concave around the middle point (y=0). In addition, the inhomogeneity of σyy is such that it even can change sign at moderate H. From Figure 16 it is clearly seen that, despite the fact that the Ni layers suffer shortening in the direction of the magnetic field, the stress in the Ni layers is tensile (σxx>0). It is again a result of the functionality of the Ni material in which the resulting mechanical stress is influenced by both factors such as the elastic strain and magnetic field. On the contrary, in the direction of the sample width (perpendicular to the magnetic field), Ni layers suffer stretching but may have both the tensile stress (σyy>0) and compressive stress (σyy<0) at moderate magnetic fields and finally compressive stress (σyy<0) at strong magnetic fields, see Figure 17. It is clear that this inhomogeneity results from the demagnetization and can certainly affect the ME coupling. The stress induced in the PZT layer is of vital interest, since it determines the induced polarization in the PZT layer and, thus, the ME coupling. This stress is compressive (σxx<0) in the direction of the magnetic field and tensile (σyy>0) in the direction perpendicular to the magnetic field, and their values notably depend on the size of the sample, namely, on the ratio between length and width of the ME laminate and its thickness.The continuity ofσzz at the interlayer boundaries is seen in Figure 18 as it should occur for the other normal stress components (σxz and σyz in the considered case). At the same time, σxx and σyy are discontinuous at the above interfaces, as seen in Figures 15–17. Note that the stress does not simply follow the strain due to the hallmark of functional structure. Figures 15–18 show that, in accordance with an assumption of the homogeneous fields theory [1–3, 5–8], σzz is negligibly smaller than σxx and σyy. In addition, on the base of Figure 15(a), we infer on approximate obeying of the zero integral-traction condition employed in the analytic theory [1–3, 5–8]. On the other hand, Figures 8, 9, 13, and 15 demonstrate an essential nonhomogeneity of some important physical fields due to the finite sizes of the ME laminate structure. This nonhomogeneity, being a distinctive feature of the first-principles simulations, essentially contributes to the ME voltage coefficient.


### 5.4. Comparison of Simulations and Experimental Data

A comparison of Figures17 and 22 shows a close semiquantitative agreement between the simulated and measured αME versus H dependence. It is worth to note the agreement between the values of maximal αME, the measured 25 mV/Oe against the computed 24.8 mV/Oe, and the values of the magnetic field Hmaxat which the maxima are attained, the experimental Hmax=11 kA/m versus Hmax=12 kA/m predicted by the COMSOL simulations. This agreement is quite expectable since the H-B curve used in the simulations has been adjusted in order to fit the computed strains to the measured ones, see Figure 21.It is commonly accepted thatHmax is close to the magnetic saturation field Hsat in the magnetosrictive material. However, for Ni samples, according to the literature data [15, 16, 20] and our measurements, Hsat≈ 80–100 Oe (6.32–7.90 kA/m). At the same time, our experimental Hmax is notably larger. Note that such a shift of Hmax in the Ni-based ME laminates to the fields well above Hsat for Ni samples, and close to the values obtained in this work, was reported [7, 14] but was not explained before. We suggest that this shift proves that the Ni layer, bonded into the laminate structure embedded into a magnetic field, is under tension. It is known that the compressed bulk Ni samples exhibit a shift of Hsat to higher fields [15, 16, 20]. In particular, an Hsat shift comparable to that observed for Hmax in this work was produced by a prestress of ~1 MPa [15, 16, 20]. It can be seen from Figure 15 that the maximum stress in the Ni layer under H=54.7 kA/m is about 1.5 MPa, which shows that this mechanism is quite feasible.To the best of our knowledge, the changes in magnetic characteristics of the magnetostrictive layer embedded in an ME laminate as compared to those of the corresponding magnetostrictive material, which emerge as a result of inevitably stressed state of the layer, have not been considered previously. Thus, as it has been noted in Section1, the phenomenological piezomagnetic coefficient q cannot be simply referred to as being the magnetostrictive material property but rather a property of the ME laminate as a whole. Though we do not calculate it, the latter real q is smaller than the former, ideal, q [1–3, 5–8]. According to the qualitative estimation of q noted in Section 1, the actual reduction of q is due to both the above-noted decrease of the saturated magnetostriction (Figure 20) and the increase of Hmax (Figure 22) as compared to bare Ni material.


## 5.1. Simulation Results

Figure7 presents the computed magnetic flux density distribution and magnetic flux streamlines. Figure 8 presents the profiles of the x-component of H on the centerlines in middle plane inside both the Ni upper layer and PZT layer.Figure 7
Magnetic flux density streamlines over magnetic flux map.Figure 8
Lateral variation of the magnetic field componentHx on the centerline in the middle planes of the Ni layers (the magenta line) and PZT layer (black line); H(ext)=54.7 kA/m.Figure 9
Electric field streamlines over the electric potential map in the middlex-z plane at H(ext)=54.7 kA/m.Variation of the electric potential withz-coordinate (x=y=0) and L=16 mm. (a): t=2 mm, w=2 mm; (b): t=2 mm, w=6 mm; (c): t=2 mm, w=16 mm; (d): t=1 mm, w=6 mm; (e): t=3 mm, w=6 mm; (f): t=4 mm, w=6 mm.
(a)(b)(c)(d)(e)(f)Variation of the electric field componentEz with (a) x-coordinate and (b) y-coordinate.
(a)(b)Variation of the elastic displacementux (a), uy (b) for various magnetic fields, as shown in the legend, along the PZT layer centerline of the middle plane.
(a)(b)Figure 13
Strainεxx along piezoelectric and magnetostrictive layers: the middle planes and centerline at H(ext)=42 kA/m.Figure 14
3D strainεxx map over displacement field obtained at H=3.95 kA/m.Stresses in the piezoelectric and magnetostrictive layers, middle planes and centerline: (a)σxx versus the length coordinate x; (b) σyy versus the width coordinate y. Curves # 1–4 correspond to the Ni layer and # 5–8 to the PZT layer; the increasing number of the curves correspond to the four increasing magnetic field values: 10.2, 22.9, 35.7, and 54.7 kA/m.
(a)(b)Stressσxx map in the middle plane (x-z) at (a) H=10.2 kA/m and (b) H=54.7 kA/m.
(a)(b)Stressσyy map in the middle plane y-z, L=16 mm, W=6 mm, (a, b): t=1 mm (c, d): t=2 mm; (e, f): t=4 mm.
(a)(b)(c)(d)(e)(f)Figure 18
Stressσzz through the piezoelectric and magnetostrictive layers versus the z-coordinate (x=y=0).Figure 19
SimulatedαME as a function of H.Figure 20
Schematics of the setup for measuring the ME response and elastic deformation.Figure 21
Measurement setup.Figure 22
The measured magnetostrictive strainsεxx and εyy for the Ni layer and Ni/PZT/Ni laminate as a function of H.The coupled ME simulations have been carried out at various values of external magnetic fieldH(ext) meaning the magnetic field produced by the permanent magnets in the absence of the laminate. In what follows, the values of Hx(ext) in the center of the laminate used in the simulations are referred to as simply “magnetic field H.” To avoid misunderstanding, the field H noted in connection with Figures 7 and 8 is induced by H(ext), the value of which is noted in the caption.Examples of the electrostatics part of our simulations at values ofH(ext) between 10.2 and 54.7 kA/m are presented in Figures 9–11. Figure 9 shows the E streamlines over the φ map on the central x-z plane. Figure 10 shows the φ profiles on the z-axis (the middle of the laminate, x=y=0) for different widths and total thickness of ME laminates at different magnetic fields. Figures 11(a) and 11(b) show the variation of Ez as a function of x and y, respectively, on the centerlines of the PZT layer.Figures12(a) and 12(b) show the ux and uy displacement components (along the length and width of the laminate structure), correspondingly, taken in the center of the piezoelectric layer for various values of H as shown in the figure. Figure 13 shows the strain component εxx for both the piezoelectric and magnetostrictive layers; εxx is calculated as the numerical derivative of ux with respect to x, see Figure 12(a). Figure 14 shows the 3D distribution map of εxx and the volume displacement field over it.Figures15(a) and 15(b) show the stress components σxx and σyy, respectively, in both the PZT layer and the top Ni layer, whereas Figures 16 and 17 show the stress σxx map (middle plane x-z) and the stress σyy map (middle plane y-z) correspondingly for several total thicknesses of ME laminates at various magnetic fields. Figure 18 presents the stress components σxx over the laminate.To estimate the modulation of the computed physical parameters, such as net strains and the voltage between the Ni layers, with respect to the magnetic field using the developed COMSOL approach, we swept parametrically the values ofM0 for the permanent magnets in the range 300–4300 kA/m. This procedure resulted in a parametric sweep of the magnetic field in the range 3.5–55 kA/m. Then, to obtain the ME voltage coefficient αME, as prescribed by (1), the output voltage was numerically differentiated by H. We especially focus here on the magnetostriction effect and the ME voltage, because we determined these parameters also experimentally. The quantity of most importance and final target of our simulations is, of course, αME which determines the ME response of the laminate and is shown in Figure 19.


## 5.2. Experimental Setup and Measurements

Our raw materials were polycrystalline ferromagnetic nickel (Ni) and piezoelectric ceramic lead zirconate titanate (PZT); 99.95% purity Ni and PZT-854 plates were obtained from Alfa Aesar and APC, respectively. The relevant data sheets of the materials are displayed in Table1. The Ni/PZT and Ni/PZT/Ni laminates’ samples were fabricated by bonding using Hernon 315 glue. For independent extraction of εxx and εyy, the strain measurements on the Ni layer and Ni/PZT bilayer were carried out using two-axis SR-4 strain gauges by Vishay Micro-Measurements. For the measurement of αME, we employed the setup developed in our lab by using a mechanical apparatus enabling accurate positioning of permanent magnets symmetrically located along the laminate axis, producing stable, variable, and reproducible magnetic bias fields (Hdc) at the sample location. An ac magnetic field with constant amplitude of hac=79,6 A/m at varied frequency was applied by a solenoid and buffer configuration amplifier AD549 with 0.11 fA current noise density. Data was collected with PXI measurement system driven by a LabVIEW program as shown in Figure 20. Figure 21 displays real view of that part of the setup, which corresponds to our computational environment (see Figure 5). Figure 22 shows the magnetostrictive εxx and εyy versus H measured on the standalone Ni (Alfa Aesar 99.95% purity) plate and Ni/PZT/Ni laminate. Note that the laminate’s |εxx | and |εyy| are overall smaller than the corresponding values for the Ni alone. The measured and simulated strains εxx and εyy as a function of H are shown in Figure 23. The measured αME versus H curve is shown in Figure 24.Table 1
Material parameters (compliance coefficientSE, piezoelectric charge constant d, and relative dielectric constant) for PZT [22] and Ni [15].
Properties/materials1110−12 (m2/N)s1210−12 (m2/N)d3110−12 (C/N)κσ/κ0Ni4.57− 1.37——PZT 85416.9−5.08−2403400Figure 23
The measured magnetostrictive strainsεxx and εyy as a function of H: the simulation (the green and blue lines) versus the measurements (the red and cyan markers).Figure 24
The measuredαME as a function of H.


## 5.3. Discussion of the Simulations Results

From Figures7 and 8, one can clearly see the demagnetization effect locally, rather than on average as in the treatment of our previous paper [10]. It is seen that, as is expected [15], the demagnetization effect is pretty much inhomogeneous and strongest near the edges. Figure 8 clearly demonstrates obeying the magnetostatics boundary conditions, namely, the continuity of the normal component of B and tangential components of H, throughout the simulations.It can be seen from Figures9 and 11 that the electric field is inhomogeneous across the structure, in particular one can see essential fringing field near the edges of the laminate, which is clear manifestation of the finite-size effects, whereas the field is nearly homogeneous around the center of the laminate. Figure 10 demonstrates that the COMSOL simulation obeys the grounding condition and confirms a charge accumulation on the electrodes’ surfaces (the φ curve has sharp kinks). Note that these kinks are determined by both (electric and elastic) parts of D. Also note that there exists an electric field in the laminate’s surroundings, though it is essentially smaller than the field inside. The external electric field appears to be unfavorable from the energetic viewpoint as it reduces the ME response defined by the internal electric field. The magnitude of that surrounding electric field depends on the ME laminate’s lateral dimensions and becomes very small if the total thickness is negligible relative to the lateral sizes of the ME laminate, see Figure 10(c). Thus, it can be argued that the boundary condition Dn=0 adopted in the homogeneous fields theory [5–7] does correspond to the approximation of an infinite plate, as it was claimed. Moreover, the simulations obtained in the present study show that the magnitude of the ME response predicted previously by the homogeneous fields theory [5–7] seems rather a lower than upper limit of the observed ME effect, in contrast with the conclusion in [12]. Further estimation, using Figures 9 and 10, of the electric field flux through any surface that surrounds the laminate, but inside the air box, shows existence of a nonzero net charge. It can be understood that this net charge appears on the bottom (grounded) electrode due to stresses induced by the magnetic field inside the PZT layer.Figures12(a), 12(b) and Figure 13 show nearly constant slopes of the components of u with respect to the coordinates deep inside the structure. The ux and uy profiles clearly show a contraction of the laminate structure along the applied magnetic field (the x-axis) and an elongation in the orthogonal direction (the y-axis).While foruy such a behavior continues to the ends of the laminate, ux shows deviation from the linearity near the edges. Figures 13, 14, and 15(a) show that the variation of εxx and σxx is rather slow inside the layers, and it is much more rapid approaching the edges. This behavior can be also detected for the polarization induced in the PZT layer which is of interest in this problem. So the complete electroelastic state of the PZT layer can be represented as a sum of two electroelastic states: the inner electroelastic state and that boundary (transition) layers. The inner electroelastic state varies relatively slowly along the coordinate lines of the middle plane. On the contrary, the boundary layers damp down quickly in the directions perpendicular to the edges and are described by truly 3D electroelasticity equations with appropriate boundary conditions. The FEM simulations in this study properly take account of the “edge effects” and show that the correction due to the boundary layers may introduce notable changes into the description of the internal electroelastic state of the ME laminate. Figures 15–17 are a true hallmark of the finite-size effects. Comparing Figures 15(a) and 15(b) shows that σxx≫σyy, which is an interplay of both the geometry relation L/W>2 and the difference in the magnetostrictive strains ε∥(ms)>ε⊥(ms). While in the PZT layer the curves describing both stresses are concave at all H used, the σyy curves in the Ni layer evolve from being convex in the region of the edges (y=±3 mm) to concave around the middle point (y=0). In addition, the inhomogeneity of σyy is such that it even can change sign at moderate H. From Figure 16 it is clearly seen that, despite the fact that the Ni layers suffer shortening in the direction of the magnetic field, the stress in the Ni layers is tensile (σxx>0). It is again a result of the functionality of the Ni material in which the resulting mechanical stress is influenced by both factors such as the elastic strain and magnetic field. On the contrary, in the direction of the sample width (perpendicular to the magnetic field), Ni layers suffer stretching but may have both the tensile stress (σyy>0) and compressive stress (σyy<0) at moderate magnetic fields and finally compressive stress (σyy<0) at strong magnetic fields, see Figure 17. It is clear that this inhomogeneity results from the demagnetization and can certainly affect the ME coupling. The stress induced in the PZT layer is of vital interest, since it determines the induced polarization in the PZT layer and, thus, the ME coupling. This stress is compressive (σxx<0) in the direction of the magnetic field and tensile (σyy>0) in the direction perpendicular to the magnetic field, and their values notably depend on the size of the sample, namely, on the ratio between length and width of the ME laminate and its thickness.The continuity ofσzz at the interlayer boundaries is seen in Figure 18 as it should occur for the other normal stress components (σxz and σyz in the considered case). At the same time, σxx and σyy are discontinuous at the above interfaces, as seen in Figures 15–17. Note that the stress does not simply follow the strain due to the hallmark of functional structure. Figures 15–18 show that, in accordance with an assumption of the homogeneous fields theory [1–3, 5–8], σzz is negligibly smaller than σxx and σyy. In addition, on the base of Figure 15(a), we infer on approximate obeying of the zero integral-traction condition employed in the analytic theory [1–3, 5–8]. On the other hand, Figures 8, 9, 13, and 15 demonstrate an essential nonhomogeneity of some important physical fields due to the finite sizes of the ME laminate structure. This nonhomogeneity, being a distinctive feature of the first-principles simulations, essentially contributes to the ME voltage coefficient.


## 5.4. Comparison of Simulations and Experimental Data

A comparison of Figures17 and 22 shows a close semiquantitative agreement between the simulated and measured αME versus H dependence. It is worth to note the agreement between the values of maximal αME, the measured 25 mV/Oe against the computed 24.8 mV/Oe, and the values of the magnetic field Hmaxat which the maxima are attained, the experimental Hmax=11 kA/m versus Hmax=12 kA/m predicted by the COMSOL simulations. This agreement is quite expectable since the H-B curve used in the simulations has been adjusted in order to fit the computed strains to the measured ones, see Figure 21.It is commonly accepted thatHmax is close to the magnetic saturation field Hsat in the magnetosrictive material. However, for Ni samples, according to the literature data [15, 16, 20] and our measurements, Hsat≈ 80–100 Oe (6.32–7.90 kA/m). At the same time, our experimental Hmax is notably larger. Note that such a shift of Hmax in the Ni-based ME laminates to the fields well above Hsat for Ni samples, and close to the values obtained in this work, was reported [7, 14] but was not explained before. We suggest that this shift proves that the Ni layer, bonded into the laminate structure embedded into a magnetic field, is under tension. It is known that the compressed bulk Ni samples exhibit a shift of Hsat to higher fields [15, 16, 20]. In particular, an Hsat shift comparable to that observed for Hmax in this work was produced by a prestress of ~1 MPa [15, 16, 20]. It can be seen from Figure 15 that the maximum stress in the Ni layer under H=54.7 kA/m is about 1.5 MPa, which shows that this mechanism is quite feasible.To the best of our knowledge, the changes in magnetic characteristics of the magnetostrictive layer embedded in an ME laminate as compared to those of the corresponding magnetostrictive material, which emerge as a result of inevitably stressed state of the layer, have not been considered previously. Thus, as it has been noted in Section1, the phenomenological piezomagnetic coefficient q cannot be simply referred to as being the magnetostrictive material property but rather a property of the ME laminate as a whole. Though we do not calculate it, the latter real q is smaller than the former, ideal, q [1–3, 5–8]. According to the qualitative estimation of q noted in Section 1, the actual reduction of q is due to both the above-noted decrease of the saturated magnetostriction (Figure 20) and the increase of Hmax (Figure 22) as compared to bare Ni material.


## 6. Conclusions

The ME effect in the laminated composites of the magnetostrictive Ni plates stacked and perfectly bonded to the piezoelectric PZT plate, assumed insulating, has been simulated by employing a commercial FEM package COMSOL Multiphysics V4.1. The numerical results obtained in this study are in a good qualitative agreement with the numerical results reported by other groups for the ME laminates including different magnetostrictive layers [11–13].Unlike these works, we do not use the linear constitutive relations for the magnetostrictive layer (in our case, Ni). We rather perform the simulations using an initial prestrain of the magnetostrictive origin, which is modeled by a quadratic function of the magnetization. This method initially proposed in [19] for modeling the magnetostriction in Galfenol beams has been employed in the present work to model the ME Ni/PZT/Ni laminate. In addition, our modeling is partially based on experiment as, to compute the prestrain, we used an M(H) curve reconstructed from strain versus H measurements as compared to that of standalone Ni material to fit the measured magnetostriction.The simulation results for the ME response versus the magnetic field have been validated by thorough comparison with the experimental data obtained in this work and those reported by other groups [7, 14]. A close overall agreement between the simulations and experimental data has been stated. The present study has revealed the fact which seems of importance for designing the ME laminates. Namely, the simulations have shown the existence of an electric field surrounding the ME laminate, which may have undesirable effect on the magnetoelectric coupling.The FEM simulations reported in this work have been shown to give robust, physically reasonable, and accurate results. It can surely be argued that this approach is capable of solving a wide range of the ME laminates’ designs with different geometries, materials, and loading configurations.



---

*Source: 383728-2012-04-02.xml*