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.
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 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 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, 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 is the piezomagnetic coefficient . Often it is tacitly assumed that 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 at magnetic saturation divided by the corresponding magnetic saturation field. However, this theoretical prediction does not agree well with the values of retrieved from the measurements of ; see [10, 14]. This discrepancy is overcome by introducing a correction factor for [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 explicit 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. Section 2 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
Figure 1 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.
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. Figure 2 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 is just the manifestation of the product ME effect in these specific structures.
In sensor applications, it is not that is measured, rather it is the change in the response to a modulated small-field signal which superimposes the bias magnetic field . 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: 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 field . It is the total field 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 (), normalized magnetic potential (), magnetization (), electric polarization (), 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 Figure 3. 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.
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 (Figure 3).
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, where is the magnetic flux density and 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 (); for example, for permanent magnets, we used in the range 300–4300 kA/m and fixed magnetic permeability . 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 , , and 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 dependence using a 99.95% Ni foil supplied by Alfa Aesar. The obtained 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, is plugged into simulations via the so-called curve [15] given by
It appears that the magnetostrictive strains obtained numerically for the ME laminate by plugging the 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 curve, reconstructed from versus measurements. Such a modified curve used in our simulation is shown in Figure 4.
Both the permanent magnets and the ME structure are surrounded by air for which . In the course of these simulations, COMSOL observes the well-known magnetostatics boundary conditions, that is, the continuity of normal component and tangential components of . In addition, magnetic insulation condition () 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, where is the th component of the relative magnetostrictive deformation along the direction , being the full magnetostrictive strain tensor, is the polycrystalline magnetostriction constant which is expressed through the magnetostriction constants of the corresponding crystal along two principal cubic axes [15, 20], , and being the th component of local magnetization in a grain/domain and the saturation magnetization, respectively; the bar means averaging over grains/domains. The expression for 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 () and perpendicular to the field () that is unique for polycrystalline material. Furthermore, and at magnetic saturation. Knowing the magnetization curve which is in fact , COMSOL allows us to calculate but, unfortunately, does not allow the computation of neither , nor . Therefore, to build a magnetostriction model in the COMSOL environment, we use the following approximation [19]: 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: kA/m and .
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: where is the strain tensor, is an initial strain one, is the stress tensor, is an initial stress, and is the elasticity matrix. The elements of the latter matrix are calculated for the Ni layers using Young’s modulus GPa and Poison’s ratio 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, , the expression for which is discussed previously, and . The coupled governing equations (6) for the laminate are solved with the traction-free boundary conditions () at outer surfaces and the perfect bonding at all the interfaces of the structure (continuity of and ). 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]: where is the electric field, is the electric displacement, is the remanent polarization (nonzero for a ferroelectric), is the compliance coefficients’ matrix at given electric field, 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 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 Figure 3) 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 () and the electric potential (): 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 with the surface charge density on the electrodes, that is, , and the continuity of the tangential components of 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 , is fully determined from the and 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 () 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 the 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 Section 3 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: mm, , 6, and 16 mm , 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 mm, 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 Figure 5. 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.
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 Section 3.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 (). 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 Figure 6.
(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
Figure 7 presents the computed magnetic flux density distribution and magnetic flux streamlines. Figure 8 presents the profiles of the -component of on the centerlines in middle plane inside both the Ni upper layer and PZT layer.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
The coupled ME simulations have been carried out at various values of external magnetic field meaning the magnetic field produced by the permanent magnets in the absence of the laminate. In what follows, the values of in the center of the laminate used in the simulations are referred to as simply “magnetic field .” To avoid misunderstanding, the field noted in connection with Figures 7 and 8 is induced by , the value of which is noted in the caption.
Examples of the electrostatics part of our simulations at values of between 10.2 and 54.7 kA/m are presented in Figures 9–11. Figure 9 shows the streamlines over the map on the central plane. Figure 10 shows the profiles on the -axis (the middle of the laminate, ) for different widths and total thickness of ME laminates at different magnetic fields. Figures 11(a) and 11(b) show the variation of as a function of and , respectively, on the centerlines of the PZT layer.
Figures 12(a) and 12(b) show the and displacement components (along the length and width of the laminate structure), correspondingly, taken in the center of the piezoelectric layer for various values of as shown in the figure. Figure 13 shows the strain component for both the piezoelectric and magnetostrictive layers; is calculated as the numerical derivative of with respect to , see Figure 12(a). Figure 14 shows the 3D distribution map of and the volume displacement field over it.
Figures 15(a) and 15(b) show the stress components and , respectively, in both the PZT layer and the top Ni layer, whereas Figures 16 and 17 show the stress map (middle plane ) and the stress map (middle plane ) correspondingly for several total thicknesses of ME laminates at various magnetic fields. Figure 18 presents the stress components 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 of 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 , as prescribed by (1), the output voltage was numerically differentiated by . 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, 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 Table 1. The Ni/PZT and Ni/PZT/Ni laminates’ samples were fabricated by bonding using Hernon 315 glue. For independent extraction of and , 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 , 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 () at the sample location. An ac magnetic field with constant amplitude of 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 and versus measured on the standalone Ni (Alfa Aesar 99.95% purity) plate and Ni/PZT/Ni laminate. Note that the laminate’s | | and || are overall smaller than the corresponding values for the Ni alone. The measured and simulated strains and as a function of are shown in Figure 23. The measured versus curve is shown in Figure 24.
5.3. Discussion of the Simulations Results
From Figures 7 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 and tangential components of , throughout the simulations.
It can be seen from Figures 9 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 . 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 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.
Figures 12(a), 12(b) and Figure 13 show nearly constant slopes of the components of with respect to the coordinates deep inside the structure. The and profiles clearly show a contraction of the laminate structure along the applied magnetic field (the -axis) and an elongation in the orthogonal direction (the -axis).
While for such a behavior continues to the ends of the laminate, shows deviation from the linearity near the edges. Figures 13, 14, and 15(a) show that the variation of and 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 , which is an interplay of both the geometry relation and the difference in the magnetostrictive strains . While in the PZT layer the curves describing both stresses are concave at all used, the curves in the Ni layer evolve from being convex in the region of the edges ( mm) to concave around the middle point (). In addition, the inhomogeneity of is such that it even can change sign at moderate . 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 (). 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 () and compressive stress () at moderate magnetic fields and finally compressive stress () 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 () in the direction of the magnetic field and tensile () 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 at the interlayer boundaries is seen in Figure 18 as it should occur for the other normal stress components ( and in the considered case). At the same time, and 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], is negligibly smaller than and . 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 Figures 17 and 22 shows a close semiquantitative agreement between the simulated and measured versus dependence. It is worth to note the agreement between the values of maximal , the measured 25 mV/Oe against the computed 24.8 mV/Oe, and the values of the magnetic field at which the maxima are attained, the experimental kA/m versus kA/m predicted by the COMSOL simulations. This agreement is quite expectable since the 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 that is close to the magnetic saturation field in the magnetosrictive material. However, for Ni samples, according to the literature data [15, 16, 20] and our measurements, ≈ 80–100 Oe (6.32– kA/m). At the same time, our experimental is notably larger. Note that such a shift of in the Ni-based ME laminates to the fields well above 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 to higher fields [15, 16, 20]. In particular, an shift comparable to that observed for 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 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 Section 1, the phenomenological piezomagnetic coefficient 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 is smaller than the former, ideal, [1–3, 5–8]. According to the qualitative estimation of noted in Section 1, the actual reduction of is due to both the above-noted decrease of the saturated magnetostriction (Figure 20) and the increase of (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 curve reconstructed from strain versus 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.
Appendix
A. Vector Notation and Linear Material Constants for Lead Zirconate Titanate (PZT-5H)
Vector notations of the strain and stress tensors are as follows: where is the engineering shear strain, The PZT-5H compliance coefficients in matrix form are those listed in [23]: The PZT-5H piezoelectric charge constant in matrix form is that listed in [23]: Relative permittivity [23] in matrix form is as follows:
Acknowledgments
The authors are grateful to M. Auslender (Ben-Gurion University) for supplying so much invaluable information and fruitful discussion of the results and I. Felner (Hebrew University) for helping us in the magnetization curves measurements for several magnetostrictive materials.