Effective nonlinear responses of three-phase magnetoelectric composites

A computational method, dubbed simplified unit-cell micromechanics model, is generalized and applied to establish the effective nonlinear responses of three-phase magnetoelectric composites that are composed of two distinct magnetostrictive and piezoelectric phases embedded in elastic polymer matrices. The nature of nonlinear constitutive behavior of each constituent is expected to significantly influence the overall responses of the composites. To obtain the effective nonlinear responses, a mathematical linearization is first introduced to perform the constitutive linearization for the nonlinear materials, and the resulting constitutive equations are then unified and nested into the micromechanics model followed by iterations in order to minimize errors from the linearization process. For the purpose of comparison, we also reformulate the well-established Mori–Tanaka micromechanics model insofar as its mathematical structure is aligned with that of the simplified unit-cell model. Numerical results are first validated against limited experimental measurements available in literature. Parametric studies are then conducted in order to reveal the effect of phase constitutive laws, volume fractions, and geometries on the overall nonlinear responses of there-phase magnetoelectric composites. The contributions of this work complement those of earlier studies that prevalently devoted to two-phase magnetoelectric composites and linear magneto-electro-elastic coupled responses only.


Basic equations of nonlinear magneto-electro-elasticity
A three-phase composite concerned in this study is constituted by piezoelectric, magnetostrictive, and polymeric phases. Herein, we make a simplifying assumption, for analytical purposes, that the composite is under small deformations and within the scopes of electrostatics and magnetostatics (i.e., we disregard the transient response of a composite). The constitutive relations of these three phases are summarized here, respectively. First, for piezoelectric phase, this study is restricted in a polarized piezoelectric material subject to a large driving electric field but without exceeding the coercive field. In such case, it is reasonable to assume that the effect of depolarization and polarization switching are negligible during the operation for actuation or sensing. As a result, the nonlinear constitutive laws derived by Tiersten 24 is utilized in this study as shown below (1)  www.nature.com/scientificreports/ where the field variables are stress σ ij , strain ε ij , electric displacement D i , and electric field E i . The material properties are the elastic stiffness C ijkl , the third-order and fourth-order piezoelectric stress constants e ijk and b ijkl respectively, and the second-order and third-order permittivities κ ij and χ ijk respectively. Second, for magnetostrictive phase, this study, similarly, assumes a magnetostrictive material subject to a large driving magnetic field but without exceeding the coercive field. Under the circumstances, it is reasonable to presume that the effect of magnetization/demagnetization and domain reorientations are negligible during the operation for actuation or sensing. Consequently, the nonlinear constitutive laws proposed by Carman and Mitrovic 25 is employed here given as where B i and H i are magnetic flux density and magnetic field, respectively. q ijk and x ijkl are the third-order and fourth-order piezomagnetic stress constants, respectively, while µ ij and ω ijk are the second-order and thirdorder magnetic permeabilities, respectively. Third, for polymeric phase, it is assumed as a conventional solid which is linear behavior among all elastic, electric, and magnetic fields and is uncoupled among any of them. It is noted that, although C ijkl in Eqs. (1) and (3) are the same symbol, the elastic stiffness are in general different for different phases.
Here, a first-moment secant linearization is applied to linearize the nonlinear constitutive equations. For example, Eqs. (1) and (2)  where With the linearized constitutive equations, Eqs. (5), (6), (9) and (10), in hand, a set of constitutive equations that directs the nonlinear interaction of magnetic, electric and elastic fields in a magnetoelectric medium can be derived through a mathematical superposition of those equations that share common strain and stress field variables. Subsequently, magnetoelectric coupling tensors a ij and ij are instinctively derived from the resultant equations, given as Magnetoelectric coefficients are caused by mechanically coupled magnetostrictive and piezoelectric phases in a three-phase composite: it is present in neither phases individually.
In addition, the kinematic equations that relate elastic displacement u i , electric potential φ i , and magnetic potential ϕ i to strain, electric field, and magnetic field, respectively, are given by: www.nature.com/scientificreports/ where the comma notation indicates a derivative. Moreover, the equilibrium equations that satisfy mechanical equilibrium and the conservation of magnetic and electric fluxes in a magneto-electro-elastic solid, in absence of body forces and free charges, are defined by: Equations (13)-(21) summarizes the governing equations of a quasi-static magneto-electro-elastic problem for a magnetoelectric composite containing nonlinear constituents.
For coupled multi-physics problems, it will be convenient in the sequel to define the vectors and Z as follows: such that a set of the constitutive equations in Eqs. (13)- (15), consequently, can be unified into a single equation by using Voigt and Nye's contracted notations, i.e., where L(E, H) is a 12 × 12 matrix. It contains field-dependent (or arithmetic) material properties that are the outcomes of the linearization process. The parenthesis used in Eq. (24) explicitly indicates that the arithmetic matrix L is indeed a function of fields E and H , and should not be confused with a multiplication operation. As convention, a boldface letter is used to denote a vector, matrix or tensor of any order, while a lightface letter is utilized to represent a scalar or a component of a tensor. The superscript T in Eqs. (22) and (23) stands for transpose operation. Equation (24) will be utilized throughout the subsequent micromechanics analysis that deals with the nonlinear responses of three-phase magnetoelectric composites.

Micromechanics formulations for three-phase magnetoelectric composites
The simplified unit-cell micromechanics model that simulates three-phase magnetoelectric composites with 1-3, 0-3 and 2-2 connectivities is first formulated in this section. For the purpose of comparison, we further reformulate the Mori-Tanaka micromechanics model insofar as its mathematical structure is aligned with that of the simplified unit-cell model.
Simplified unit-cell model. The simplified unit-cell model has been recently employed by Lin and Lin 26 for the prediction of the nonlinear behavior of two-phase magnetostrictive-piezoelectric composites, and by Lin and Liu 27 for the simulation of three-phase polymer matrix smart composites. Lin and Lin 28 39 . In these papers, extensive comparisons with experimental data or other computational approaches, i.e., finite element, Mori-Tanaka etc., have been performed in order to verify the reliability of the two-phase unit-cell model. In the present article, the simplified unit-cell model is further generalized and applied to establish the effective nonlinear responses of three-phase magnetoelectric composites. A three-phase composite microstructure is first idealized as a periodically distributed array that is constructed by elementary cubes shown in Fig. 1. A unit cell is then defined by a cubic array to the extent that it is able to generate the overall microstructure by means of repeating itself. For a three-phase composite, the simplest unit cell is therefore identified as a collection of 64 subcells. The unit cell here can also be viewed as an idealized representative volume element whose effective behavior is representative of that of a composite material as a whole. In the case of a 0-3 magnetoelectric composite, as shown in Fig. 1, the magnetostrictive and piezoelectric particles (blue and red ones) occupy the subcell 1, 3, 9, 11, 33, 35, 41 and 43 in a spatially equivalent manner, while the polymer matrix (white one) engages the rest of the subcells. The field of each subcell is regarded as an average field; as a result, the expression of the average field of a magnetoelectric composite can be derived by volume-weighted average over the entire unit cell, that is www.nature.com/scientificreports/ An overbar denotes that a quantity is in an average sense. The superscript (α) represents the identity of a subcell and is named as a set of consecutive integers, 1-64. c (α) is the volume fraction of the subcell α and is determined by is the volume of the subcell α . Depending on each constituent that occupies specific subcells, the volume fraction of each phase in a three-phase magnetoelectric composite can be calculated. The constitutive equation of a unit cell is defined as follows: where L is regarded as the effective property of a three-phase magnetoelectric composite. The constitutive equation of a subcell is of course given by The volume average of the local field in the subcell α is related to the applied overall field by a concentration factor tensor A (α)40 below Substituting Eq. (29) into Eq. (28) to replace Z (α) , we obtain Substituting for (α) from Eq. (30) into Eq. (26), we arrive at the expression for the effective response: Comparing Eq. (27) with Eq. (31), we arrive at the result for the effective property: www.nature.com/scientificreports/ From Eqs. (31) and (32), it is obvious that once the concentration factor tensor A (α) is determined, the overall response and average property L of a three-phase magnetoelectric composite can be obtained. In order to evaluate the concentration factor tensor in the simplified unit-cell model, micromechanical relations between the subcells and constitutive models for all subcells are imposed. The micromechanical relations are derived based on the continuity conditions of displacements, tractions, electric potentials, normal electric displacements, magnetic potentials, and normal magnetic flux densities at the interfaces among the subcells. For the simplified unit-cell model having 64 subcells as shown in Fig. 1 and employing three physical domains (magnetic, electric and elastic), Eq. (29) results in 768 (64 × 12) independent variables Z (α) . Thus, we need to define 768 equations which are obtained from the micromechanical relations as listed in Supplementary 1 that also provides additional pictograms to clearly illustrate how the three-phase model simulates 1-3 and 2-2 connectivity types. The novelty of this work is that the present model can unify three typical connectivity types, i.e., 0-3, 1-3 and 2-2 with the least number of subcells, i.e., 64 subcells, into a unit cell for modeling of three-phase magnetoelectric composites. The equations can be written in a matrix form as: From Eq. (33), the concentration factor matrix A (α) in linearized relations are therefore determined by where the superscript −1 represents inverse operation. The linearized micromechanical relations are exactly satisfied only when all subcells exhibit linear constitutive responses; as a result, Eq. (34) directly provides the required concentration factor matrix to estimate the effective linear response of a three-phase magnetoelectric composite. In contrast, due to the nonlinear responses in the magnetostrictive and piezoelectric subcells, the linearized micromechanical relations will usually violate the constitutive equations. This inconsistency is defined by a residual written as Numerical iterations using computers are typically required to minimize the residual. Once the residual has been minimized, the concentration factor matrix is again determined by Eq. (34). Finally, the effective nonlinear response of a three-phase magnetoelectric composite is calculated via using Eq. (31).
At this point, one might assume the present three-phase unit-cell model is essentially the same as the traditional two-phase unit-cell model, as both of them utilize periodic parallelepipeds (subcells) to discretize a composite material. The fact is the three-phase model contains many unique features not found in the two-phase model. In Table 1 these differences are specifically identified in terms of 0-3 connectivity type. The aim of this table is to draw a clear distinction between the two-and three-phase models. For the purpose of comparison, Table 1. Comparison between the present model and some micromechanics models in terms of 0-3 connectivity type. a In present study, 0-3 connectivity means that cubic particles are embedded in a matrix; 1-3 connectivity stands for that continuous fibers having square cross sections are embedded in a matrix; 2-2 connectivity represents that a bilaminated composite having rectangular cross sections for its two phases.

Aspects
Three-phase unit-cell model (Present study) Two-phase unit-cell model 31 Mori-Tanaka model 41 www.nature.com/scientificreports/ we also include the Mori-Tanaka model, which will be presented shortly, in the Table 1. The intention here is not to undermine the well-established two-phase model and Mori-Tanaka approach or claim that they should be replaced, but rather to establish the present three-phase model as a method in its own right.
Mori-Tanaka mean field approach. The Mori-Tanaka approach is commonly employed for the predictions of the effective moduli of two-phase smart composites comprising linear constituents only. Here, we reformulate this method for magnetoelectric composites insofar as nonlinear magnetostrictive and piezoelectric constituents are concerned and three-phase inclusion-matrix composites are addressed. In a similar manner as for the average field of a unit cell in the simplified unit-cell micromechanics model, the average field of a representative volume element of a three-phase magnetoelectric composite is written as The superscript (r) indicates the identity of a phase. For example, r = 1 is polymer matrix; r = 2 is magnetostrictive inclusion; r = 3 is piezoelectric reinforcement. c (r) is the volume fraction of the rth phase and satisfies the condition, 3 r=1 c (r) = 1 . The composite constitutive equation is later defined as (in an average sense) The phase constitutive equation is given by Apply the concentration factor tensor A (r)40 to relate the field between the composite and its individual phase: Substitute Eq. (40) into Eq. (39) to replace Z (r) ; then the resulting equation is used in Eq. (37) to replace (r) ; as a result, the expression for the effective response is arrived below Compare Eq. (41) with Eq. (38); consequently the effective arithmetic stiffness of a three-phase magnetoelectric composite is obtained below Up to this point, the concentration factor tensor A (r) is the only unknown to calculate the effective response and effective property L in Eqs. (41) and (42), respectively. The Mori-Tanaka theory 41 provides a unique way to determine the concentration factor tensor. It is the primary use of the elegant Eshelby 42 formalism, based on eigenstrain concept, which is used to determine the solution to the problem of an ellipsoidal inclusion embedded in an infinite elastic matrix of material under uniform exterior mechanical loading. For a magnetoelectric material modeled by the Mori-Tanaka method, an inclusion is surrounded by the matrix with uniform strain and electric and magnetic fields the same as the matrix's averaged strain and averaged electric and magnetic fields. The inclusion's averaged fields are calculated from the solution for one ellipsoidal inclusion embedded in an infinite magneto-electro-elastic matrix. The extension from elastic problems to magneto-electro-elastic problems for three-phase composites makes the concentration factor tensor of the Mori-Tanaka approach having the following form: where A dil, (r) is the concentration factor tensor of the so-called dilute scheme (or Eshelby approach) for the rth phase of a magnetoelectric composite, given as I in Eqs. (43) and (44) is a 12 × 12 identity matrix. S (r) is the magneto-electro-elastic Eshelby tensor, expressed as      www.nature.com/scientificreports/ where A shorthand notation that treats the elastic, electric, and magnetic variables on an equal footing has adopted here. It is similar to conventional indicial notation with the exception that lowercase subscripts take on the range 1-3, while uppercase subscripts take on the range 1-5 and repeated uppercase subscripts are summed over 1-5. When an ellipsoidal inclusion is embedded in a general anisotropic magneto-electro-elastic material, S (r) is given in terms of a double quadrature as shown in Eq. (45) that has to be carried out numerically. In the present study, however, the closed-form expressions of magneto-electro-elastic Eshelby tensor for the aligned elliptic-cylindrical inclusion (1-3 connectivity), spherical inclusion (0-3 connectivity), and thin-disc inclusion (2-2 connectivity) in an isotropic elastic medium are obtained, and are shown in Supplementary 2. In the case of linear constitutive laws applied for all phases in a magnetoelectric composite, Eq. (43) immediately gives a necessitate concentration factor tensor for Mori-Tanaka's predictions. In contrast, when nonlinear constitutive laws are employed for magnetostrictive and piezoelectric phases in a multifunctional composite, the linearized constitutive relation in Eq. (24) would lead to an error estimation of the strain and electric and magnetic fields in each phase. This error can be defined by a residual as: After the residual is minimized, the Mori-Tanaka concentration factor tensor is evaluated by Eq. (43). Eventually, the effective nonlinear response of a three-phase magnetoelectric composite is obtained by the use of Eq. (41).

Numerical implementation
This section is dedicated to the study of effective nonlinear responses in three-phase 1-3, 0-3, and 2-2 magnetoelectric composites under a prescribed set of magnetic fields. We first validate the simplified unit-cell model by comparing micromechanics predictions with limited experimental data available in literature. Next, we present parametric studies in order to reveal the effect of phase constitutive laws, phase volume fractions, and composite connectivities on the overall responses of magnetoelectric composites.

Comparison with experimental data. Most of experimental investigations for magnetoelectric compos-
ites are focused on two-phase composites. Three-phase setting for such composites is relatively limited in literature. However, few smart composites reported in literature can be still treated as three-phase magnetoelectric composites. For example, Petrov et al. 43 experimentally examined the magnetoelectric effect in a porous magnetostrictive/piezoelectric bulk composite. The composite was made by PZT particles surrounded by a nickel ferrite matrix which had been modified to include porous microstructure. This composite can be regarded as a three-phase composite having PZT particles and voids in a nickel ferrite solid. The volume fraction of the PZT particles and the nickel ferrite matrix are 40% and 60%, respectively, where the 60% nickel ferrite matrix further shares some of spatial portion with the pores. The material properties required for micromechanics simulations are reported by Petrov et al. 43 as well with the exception that instead of including nonlinear material parameters, the linear material coefficients were only used by them. The two micromechanics models with nonlinear responses developed in this study should be capable in predicting the overall linear magnetoelectric properties of magnetoelectric composites. In Fig. 2, it is seen that the unit-cell predictions (red solid line and blue dashed line) of the magnetoelectric voltage coefficients, a E33 and a E31 in longitudinal and transverse directions, respectively, are in good agreement with the experimental results (red dots and blue circles) for the pore volume fraction 0-40%. Although the Mori-Tanaka estimations (red solid-cross line and blue dashed-cross line) have the same trend as the unit-cell predictions but the Mori-Tanaka approximations show a relatively large deviation from the experimental results than the unit-cell model does. It may be owing to the various morphology used (50) www.nature.com/scientificreports/ to approximate the PZT particles and pores. A cubic shape by the unit-cell model while a spherical geometry by the Mori-Tanaka scheme are utilized to approximate a particulate reinforcement. It also probably due to various number of fields employed to approximate the variation of the nickel ferrite matrix fields. Fifty-six matrix subcells enable the unit-cell model to capture some variation of the matrix fields while the Mori-Tanaka model discerns only the average matrix fields, i.e., a single value applicable to the entire matrix. It is however noted that the intention here is not to undermine the well-established Mori-Tanaka theory or claim that it should be replaced, but rather to establish the simplified unit-cell model as a method in its own right.
Parametric studies. We conduct parametric studies here in order to investigate the overall nonlinear magnetoelectric responses particularly for a composite composed of PZT-G1195 and Terfenol-D that show significant material nonlinearity. Figure 3 depicts Crawley and Anderson's 16 experimental measurements (red dots) on the in-plate strain response ε 11 in a PZT-G1195 plate due to an applied electric field E 3 through the plate thickness aligned with the poling direction. It is obvious that the nonlinear constitutive law (solid line) can capture the entire range of the response while the linear constitutive equation (dashed line) works well only for low applied field, i.e., E 3 ≈ 0.1 MV/m. As for engineering applications, it is desirable to utilize the greatest possible strain available in piezoelectric materials; as a result, a large electric driving field is needed, which results in significantly nonlinear piezoelectric responses. Similarly, Fig. 4 illustrates Jiles and Thoelke's 17 experimental results (dots) on the strain ε 33 and magnetic flux density B 3 responses in a Terfenol-D rod owing to an applied magnetic field H 3 . It is evident that the nonlinear constitutive relation (solid line) can predict the whole range of the responses while the linear constitutive law (dashed line) is useful for low magnetic field, i.e., H 3 ≈ 30 kA/m. As for practical applications, it is desirable to utilize the greatest possible operating range for applied magnetic field available in magnetostrictive materials. In such case, pronouncedly nonlinear piezomagnetic responses are observed.
The nonlinear material parameters used to fit the experimental data of the PZT-G1195 and Terfenol-D specimens are listed in Table 2. For the subsequently parametric studies, an exemplary three-phase composite is taken as the nonlinear PZT-G1195 and Terfenol-D reinforcements in a linear elastic polymer matrix Araldite D. The material property of the Araldite D used for computation is shown in Table 2 as well.  www.nature.com/scientificreports/ First of all, we investigate the nonlinear magnetoelectric responses in a three-phase composite having 1-3, 0-3 and 2-2 connectivities, respectively. The volume fractions of the PZT-G1195 and Terfenol-D inhomogeneities and the Araldite D matrix are taken of 0.25, 0.25 and 0.5, respectively. Both magnetization and electric polarization are aligned with the 3-direction that is the same as the fiber direction for 1-3 connectivity and is along with the layering direction for 2-2 type composite. Figure 5a shows the micromechanics predictions of the effective electric displacement D 3 due to an applied magnetic field H 3 for the 1-3 composite under a stress-free condition. Significant differences between nonlinear (solid lines) and linear (dashed lines) micromechanics estimations are observed. Besides, the unit-cell (without cross symbols) and Mori-Tanaka (with cross symbols) give close prediction regardless of linear or nonlinear cases in the 1-3 composite. The higher magnetic field is applied, the more  www.nature.com/scientificreports/ pronounced nonlinearity is exhibited. It implies that it is necessary to consider the material nonlinearity when a three-phase magnetoelectric composite is subjected to a large applied field and is formed of the 1-3 connectivity, but it is not the case for 0-3 and 2-2 composites as shown in Fig. 5b,c. The magnetoelectric responses for the 0-3 and 2-2 composites show not only fewer overall responses but also less mismatch between linear and nonlinear estimations than 1-3 composite does. This is due to the fact that in the 0-3 and 2-2 magnetoelectric composites, the Araldite D polymer matrix, which is a passive constituent, dominates the overall responses; as a result, only a small fraction of the applied magnetic field reaches the magnetostrictive Terfenol-D reinforcements. It is also noted that the magnetic field H 3 is applied along the fiber direction (3-direction), and the microstructural arrangement in the fiber composite allows for conducting this magnetic field directly through the Terfenol-D fibers, while in the case of the 0-3 and 2-2 composites, the inactive matrix limits the magnetic field in reaching the Terfenol-D particles and laminas, respectively. In addition, Fig. 5b depicts a relatively large discrepancy between the unit-cell and Mori-Tanaka results, regardless of linear or nonlinear cases, in the 0-3 composite in comparison with those in the 1-3 or 2-2 composites (Fig. 5a,c). It is because of the different geometries utilized in the unit-cell and Mori-Tanaka models to simulate the particle reinforcements in the 0-3 composite as discussed previously. This combination of findings supplies some support for the conceptual premise that nonlinear active constituents lead to significant overall nonlinear responses of three-phase magnetoelectric composites when the composites are subjected to large applied fields. Moreover, a three-phase magnetoelectric composite with continuous active fibers provides the best performance. Second, we examine the effect of the volume fraction of the active constituents, i.e., PZT-G1195 and Terfenol-D, on the overall responses of a three-phase composite. We keep the volume fraction of the Araldite D matrix as a constant value of 0.5 and then vary the volume fraction of the PZT-G1195 reinforcement from 0 to 0.5 against the one of the Terfenol-D inhomogeneity from 0.5 to 0. All settings of the connectivity and the direction of magnetization, electric polarization, fiber, and layering are the same as the previous parametric study except for the active phase volume fraction. Figure 6 illustrates the effective electric displacement D 3 as a function of the volume fraction of the PZT-G1195 for a stress-free magnetoelectric 1-3, 0-3 and 2-2 composite, respectively, subjected to a constant magnetic field, i.e., H 3 = 100 kA/m. It is of course that the magnetoelectric coupling vanishes when only one active phase (PZT-G1195 or Terfenol-D) is present; that is shown by the zero response of D 3 at the PZT-G1195 volume fraction of 0 or 0.5. In contrast, a relatively large magnetoelectric coupling can be obtained by evening the proportion of two active phase. In the present example, it can be seen at the PZT-G1195 volume fraction of around 0.25. This is due to that the magnetoelectric coupling is engineered by the two active phases, PZT-G1195 and Terfenol-D, in a three-phase composite. Evened the proportion of the two active phases maximizes the overall magnetoelectric coupling. As the previous finding in the first parametric study, the unit-cell and Mori-Tanaka predictions are close to each other for the 1-3 composite (Fig. 6a) regardless of linear or nonlinear cases. On the other hand, for the 0-3 and 2-2 types (Fig. 6b,c) the difference between the linear and nonlinear predictions is insignificant regardless of the unit-cell or Mori-Tanka predictions, although www.nature.com/scientificreports/ the two micromechanics estimations are deviated from each other. It is probably due to, as the foregoing discussion, only a small amount of the magnetic field reaches the magnetostrictive Terfenol-D particles and layers, respectively. These findings deduce that it is necessary to employ the nonlinear constitutive laws for the design of a three-phase magnetoelectric composite having 1-3 connectivity. These findings also raise an intriguing question regarding the prediction accuracy of the two micromechanics models applying to the cases of the 0-3 and 2-2 connectivities. The accuracy is required to be evaluated by a large amount of experimental data, which is outside the scope of this study. Finally, we explore the effect of the proportion of the passive phase, i.e., Araldite D, on the effective responses of a three-phase composite. All prescribed conditions for this simulation are the same as the last parametric study. At this moment, however, we retain a constant ratio of the volume fractions of the two active phases to be 1:1, and then modulate the volume fraction of the Araldite D matrix from 0 to 1. Figure 7 shows the effective response of electric displacement D 3 as a function of the volume fraction of the passive Araldite D matrix for the three types of connectivity, respectively. It is expected that the magnetoelectric coupling reaches the extrema at the Araldite D volume fraction of either 0 or 1, respectively. Between them, the magnetoelectric coupling descends as the proportion of the Araldite D ascends. It suggests that the maximum magnetoelectric effect can be achieved when the Araldite D matrix vanishes. However, the exclusion of the polymer matrix from a three-phase composite, it may cause a noncompliant composite that is in general not preferable for practical applications. In other words, a flexible polymer matrix makes a three-phase magnetoelectric composite more compliant at the expense of reducing overall magnetoelectric coupling. It further deduces that the present micromechanics models provide a means to reveal the overall nonlinear response, which is quite essential for tailoring a three-phase magnetoelectric composite for desirable performance. An important physical mechanism is also unveiled by all the above numerical results. The magnetoelectric effect in a three-phase composite is caused by mechanically coupled magnetostrictive and piezoelectric phases that are both surrounded by a polymer matrix: it is present in none of the individual phase separately. Under magnetic field owing to the piezomagnetism of the magnetostrictive constituent, there are mechanical strains which are elastically transmitted through the polymer intermediator onto the piezoelectric phase resulting in polarization changes via piezoelectricity.

Conclusion
Three-phase nonlinear magnetoelectric composites have been studied using the simplified unit-cell model. In development of the mathematical foundation, we utilize concepts from conventional micromechanics discipline such as representative volume elements, volume-weighted averages, and concentration factor tensors. The significance of this formulation is that a first-moment secant linearization is employed and is then nested www.nature.com/scientificreports/ into the micromechanics model followed by an iterative scheme in order to minimize the error resulted from the linearization process. The outcome of the simplified unit-cell model is the predations of effective nonlinear responses and its results are comparable to those from the well-established Mori-Tanaka model. We elucidate the similarity between the two micromechanics models in terms of not only their mathematical structures but also their numerical results. The developed unit-cell model can serve as a constitutive equation of a three-phase magnetoelectric composite material, that can be further integrated into finite element processes, via integration points, to modeling of composite structures such as beams, plates, and shells, which is a great merit in a variety of practical engineering applications. The main points found from the numerical implementation are summarized here. First, the developed simplified unit-cell model with 64 subcells has validated against experimental data. Second, the overall responses of a three-phase magnetoelectric composite is significantly influenced by nonlinear constituents particularly when the composite is subjected to large applied fields. Third, 1-3 connectivity provides the best figure of merit than the magnetoelectric composites have 0-3 or 2-2 connectivities. Fourth, including a polymer matrix lets a magnetoelectric composite to be easily shaped but would also deteriorate the overall magnetoelectric coupled performance. Fifth, the micromechanics predictions from the simplified unit-cell and Mori-Tanaka models are in general close to each other especially for 1-3 magnetoelectric composites. Taken together, the contributions of this work complement those of earlier studies that prevailingly focused on two-phase composites and linear constitutive responses, in general.

Data availability
The datasets generated during or analysed during the current study are available from the corresponding author on reasonable request.