A continuum magneto-mechanical model for magnetorheological elastomers

Magnetorheological elastomers (MREs) consist of micron-sized magnetizable particles embedded in a rubber matrix. Properties of these magneto-sensitive materials are changed reversibly upon application of external magnetic fields. They exhibit highly non-linear magneto-mechanical response which allows developing new devices and applications. However, the coupled magneto-mechanical behavior makes mathematical modeling of MREs quite complicated. So development of a reliable constitutive framework is essential for further understanding of this coupled behavior as well as simulation of the systems that utilize MREs. In this paper, a finite strain continuum model is developed for MREs where the effect of magnetization on material stiffness is directly introduced in the material shear modulus. It is shown that this approach simplifies the constitutive models and also perceives the magnetic saturation of MREs. Moreover, the coupled effects of magnetization, deformation and particle-chains orientation on the mechanical response are also taken into account in the introduced parameter. This reduces the number of material parameters, the required experimental tests for parameters identification and also simplifies the mathematical formulation of the developed constitutive equations which is beneficial for numerical formulations. A systematic two-step method is then introduced for material parameters identification which assures uniqueness of the parameters set. The predictive capabilities of the proposed model are examined via available mechanical and magneto-mechanical experimental data on both isotropic and anisotropic MRE samples at different configurations of magnetic field and loading with respect to the preferred direction of the samples. It is shown that the model can well predict the magneto-mechanical response of MREs at different deformation modes and magnetic fields.


Introduction
Magnetorheological elastomers (MREs) as a kind of magnetosensitive materials consist of micron-sized magnetizable particles (e.g. iron/iron oxides) embedded in a soft rubber matrix [1][2][3]. Properties of these magneto-sensitive materials are changed reversibly and almost instantaneously upon application of a magnetic field [4][5][6]. They exhibit a highly non-linear magneto-mechanical response which allows developing new devices such as vibration absorbers [3,7] and isolators [8], sensing devices [9][10][11][12][13], and engine mounts [14,15]. Presence of a magnetic field through the curing makes the magnetic particles form column-like arrangements which result in both mechanical and magnetic anisotropy while curing without external field leads to isotropic MREs, as shown in figure 1 [16,17].
The magnetizable particles interact mutually in effect of the magnetic field raising the overall mechanical stiffness of material which is more evident for anisotropic samples. In addition to stiffness, a magnetic field can also change natural frequency and damping properties of MREs.
The controlled deformation of magneto-active materials can be utilized in programmable matter and soft actuators [18][19][20][21][22], magnetically activated pumps [23,24], and active scaffolds for cell and drug delivery [25,26]. However, the highly coupled magneto-mechanical behavior makes mathematical modeling of MREs quite complicated. In the large strain regime, the effect of magnetic field on mechanical response of magnetic elastomers is not well recognized yet. So development of a theoretical framework is essential for further realizing of this behavior. Some works have been devoted to providing such a theoretical framework [27][28][29][30][31] and constitutive models based on those theories [32][33][34][35][36]. The adopted methods are generally classified in two categories: (a) direct methods that use continuum mechanics conservation laws, and (b) energy methods which use variation calculus to minimize a suitable potential energy. Although these approaches have been developed independently and according to their inherent assumptions they result in different representations of Maxwell stress and interface tractions, they have been shown to be equivalent [33].
Dorfmann and Ogden established a continuum theory for magneto-sensitive elastomers which uses magnetic induction/field as an independent variable rather than the magnetization [28][29][30]. They proposed a strain energy function (SEF) which results in a relatively simple constitutive equation. The theory has been applied only for isotropic magneto-elastic solids while the formulation is quite general and serves as a good starting point for developing anisotropic theories for MREs. Nedjar extended the theory by splitting the contribution of matrix material and particles chains to account for the anisotropic mechanical response [37]. Saxena et al [38]  followed a similar approach to model magneto-viscoelastic behavior of anisotropic MREs.
Kankanala and Triantafyllidis [33] developed a magnetomechanical variational formulation for the finite strain magneto-elastic problems. Their formulation obtains mechanical and magnetic equations as well as boundary/interface conditions by minimizing a potential energy.
A variety of micromechanical models have also been proposed. In such models, the microstructure and distribution of particles are also noticed in order to consider particles interactions [39][40][41][42][43].
In this paper, a continuum model is proposed for quasistatic magneto-mechanical behavior of MREs where the effect of magnetization on stiffness is introduced in the shear modulus of the material. Moreover, the coupled effects of magnetization, deformation and particle-chains orientation on the mechanical behavior of material are reflected in the same introduced parameter. This approach is followed to reduce the number of invariants and material parameters, required experimental tests for parameters identification and also to simplify the mathematical formulation of the developed constitutive equations which is of great concern for numerical formulations. A systematic two-step method is also introduced for material parameters identification and then, the predictive capabilities of the proposed model are examined via available experimental data in the literature. The model is examined for mechanical and magneto-mechanical loading of isotropic/anisotropic MRE samples at different configurations of magnetic field and loading with respect to preferred direction of the samples.

Basics of magneto-mechanics
Consider a magneto-elastic material in the reference configuration B 0 where the material points are defined by a position vector X as shown in figure 2. Upon application of a mechanical deformation combined with a magnetic field, the material deforms into a new configuration B where the point X occupies a new position x = x (X) The body deformation is described by the deformation gradient tensor F [44,45]: The magnetic field and induction vectors in the deformed configuration are shown by H and B, respectively. In vacuum or non-magnetizable materials they are related by [46]: where µ 0 stands for the magnetic permeability in vacuum. For a magnetizable material, (2) takes the following form: with the relative permeability µ r defined for that particular material. The magnetic susceptibility χ m which is closely related to the relative permeability is defined as: Using the susceptibility in (4), the magnetic induction can also be defined as: where M is the vector of magnetization. According to the Maxwell's equations, if there is no distributed currents and dependency on time, B and H should satisfy the following field equations [46]: Equations (6) apply both within and outside the material, where div and curl stand for the divergence and curl operators with respect to x, respectively. For detailed background on the electromagnetics theory, one can refer to the Jackson's classic text [47]. If there is no surface currents the continuity conditions are as follows:  (7) it can be shown that [46]: In the absence of mechanical body forces, one can treat magnetic body forces equivalently as stresses which can be added to the Cauchy stress to define a total stress, which is a symmetric tensor denoted by σ. In this case, the equilibrium equations can simply be written as [29]: According to the non-linear elasticity theory [48], a total nominal stress tensor T can be defined in relation to σ: where J = det F.

Preliminaries
In a MRE boundary-value problem, in order to derive the total stress in terms of F and an independent magnetic variable, a constitutive equation is required. According to [33], for the total stress tensor σ, we have the following expression: Here, I is the second-order identity tensor and W is the SEF per unit undeformed volume of material. Though the first term vanishes outside the material and the magnetization is zero within the free space, the total stress even outside the material would be a non-zero value called the Maxwell stress: Determination of an appropriate SEF that fits the experimental data in different modes of deformation is of concern in modeling the magneto-mechanical behavior of MREs. For an isotropic magneto-active material, the general form of W is given by: where the right Cauchy-Green deformation tensor c associated with F is defined as: For a transversely isotropic material which is the case for MREs cured in the presence of a magnetic field, W should also depend on the preferred direction which is shown by a unit vector N in the reference configuration [49]. Thus W would depend on c, N and the magnetization M: In accordance with the general theory of transversely isotropic functions [50][51][52] and also to satisfy objectivity and integrity conditions, it is preferred to define the energy function in terms of the following invariants instead of the strain components, magnetization and the preferred direction [53,54]: where the invariants are defined as: It is concluded that the characterization of MREs includes determination of an appropriate SEF as well as identification of the introduced material parameters to define the constitutive response of the material.

The constitutive magneto-mechanical model
It is common to simplify the constitutive equations as much as possible regarding the existent material symmetries and other appropriate considerations. For ordinary isotropic materials only the invariants of c are enough W = W (I 1 , I 2 , I 3 ). For incompressible materials which is a good assumption for elastomers (I 3 = 1) the SEF is independent of I 3 , which implies that the current density ρ is equal to the initial value ρ 0 . Furthermore, due to the lack of a clear physical description for I 2 , most isotropic hyperelastic models for elastomers are defined only in terms of I 1 [55][56][57]. So, for a general anisotropic magneto sensitive elastomer we can write: Exponential hyperelastic strain energies possess remarkable features in predicting the rubbery behavior and satisfying the necessary characteristics conditions [58,59]. Here we use the exponential SEF proposed by Khajehsaeid et al [60] for the isotropic part W 1 of the model: where A, a and b are material parameters. As thoroughly discussed in [59][60][61], this model reduces to the neo-Hookean model at small strains (I 1 → 3). This implies that A = G (= nkT) /2 where G, n, k and T are the small strain shear modulus, molecular chain density in unit undeformed volume, the Boltzman's constant and absolute temperature, respectively [62]. Behavior of the model at large deformations (as λ approaches the limiting chain extensibility) is dominantly governed by the parameter 'a'. So, this parameter is related to the limiting chain extensibility √ N where N is the average number of monomers in each chain segment between two adjacent cross-links [63][64][65]. On the other hand, the parameter 'b' which appears in the logarithmic term plays a significant effect at moderate strains. By increasing 'b', the material behavior approaches the mechanical behavior of elastomeric foams. Thus the logarithmic term may represent the effect of material porosity, free volume, etc [60].
In terms of the pseudo-invariants related to the material anisotropy, we recall (18) As the preferred direction N has a crucial effect on the material behavior, we need to preserve I 4 but we may be able to neglect the higher order terms for the sake of simplicity [31]: One can use the modified 'Fiber-Reinforced Elastic' model proposed in [66] for W 3 due to its simplicity. Moreover, as I 6 accounts for effect of the magnetization on the behavior of material, we will introduce a new parameter A 1 to account for the change of the material stiffness due to the magnetic field rather than including I 6 explicitly in W 4 We will show later this not only simplifies the constitutive model but also works satisfactorily in magneto-mechanical loadings as it allows capturing of magnetic saturation of the material as well. The interaction of magnetic field with mechanical deformation is comprised in the invariants I 7 and I 8 but here for simplicity we only retain I 7 and neglect the higher order terms carried in I 8 [31]: where the parameter 'r' is related to the resistance of the particle chains against deformation. This parameter depends on particles shape and surface condition as well as the volume fraction of the particles. The parameter 'k' is related to the resistance of material against deformation in the magnetization direction. It depends on magnetic permeability of the material and intensity of the applied field. One may conclude that as the parameter A is equal to the shear modulus of the material which governs the behavior at small strains [67] the new parameter A 1 reflects the change in the shear modulus due to the magnetic field intensity. The following definition is proposed for A 1 which includes the saturation effect of the magnetic field on the material stiffness: where s is a magnetic material parameter, h and h s are the applied and saturation magnetic flux densities, respectively. The saturation magnetic flux of MREs are obtained either by conventional measurement methods (e.g. Hall effect [68], induction coil [34]) or the following equation [69,70]: where ϕ is the carbonyl iron particles (CIP) volume fraction and h CIP s is the saturation magnetic flux of bulk CIP [71]. For a transversely isotropic material, invariants I 9 and I 10 reflect combined effects of magnetization, deformation and particle-chains orientation on the material behavior. For such materials, the largest effect of magnetic field is observed when the field is aligned in the particles direction. We reflect this effect by introducing the component of the magnetization vector in the preferred direction of the material (e.g. N. M), in the parameter A 1 which results in a new parameter A * corresponding to the effect of magnetic field on the stiffness of the material. This approach not only reduces the number of invariants, material parameters and the required experimental tests for material parameters identification but also leads to a simpler mathematical formulation of the developed constitutive equation which is of great concern in numerical formulation of the constitutive equations: where A * is defined as: It follows from (11)-(28) that the total nominal stress can be written as follows: where the hydrostatic pressure p should be determined from boundary conditions of the problem. It is to be noted that, for isotropic MREs there is no preferred direction but upon application of a magnetic field the particles temporarily tend to the field direction which means that the induced temporary N would be in the direction of M.

Identification of the parameters and the results
In this section, the developed constitutive model will be examined via comparisons with associated experimental data in literature. To assure uniqueness of the identified material parameters for the model, the parameters identification procedure should be done in two steps: (a) The material parameters corresponding to the nonmagnetic terms of the constitutive model are identified by means of purely mechanical tests data. (b) Keeping the parameters found in the previous step, the parameters associated with the magnetic response of material are identified by means of magneto-mechanical loading data.
Though experimental data on MREs are too limited, however some relevant data can be found in [68,72]. To the best of the authors' knowledge, these are the most comprehensive experimental studies in the field and that is why they are used here to evaluate the proposed model.

Experiments reported by Schubert and Harrison
Schubert and Harrison conducted remarkable mechanical and magneto-mechanical experimental investigations on siliconerubber based MREs in different deformation modes [68]. They tested both isotropic and anisotropic samples with different CIP contents h CIP s = 2.15 T in large uniaxial compression, tension and pure shear modes.

Isotropic samples.
We first perform material parameters identification process for the isotropic samples. In the first step we only use the compression data to determine the mechanical parameters. The least-squares optimization method (see [73]) is used to determine the parameters which are reported in table 1. The model results and the loading configuration are shown in figure 3. The identified parameters are then used to reproduce the uniaxial tension and pure shear data as also shown in figure 3. It should be noted that the tensile and shear data have not been used in the material parameters identification procedure. It is noted that, the term 'stretch' has been used as the ratio of deformed length to the initial length of material line elements. As the results show, the model is capable of predicting the mechanical behavior of the isotropic samples in different modes of deformation. Also it is observed that, the MR effect increases non-linearly with the volume fraction of particles [74]. The small discrepancy seen in the tensile regime is likely to be due to the errors in the experiments as mentioned in [68] which will be discussed later in details.
In the second step, the magneto-mechanical data will be utilized to determine the parameters associated with the magnetic response of the material. It should be noted that, for this particular data set each deformation mode has been performed at a different magnetic field intensity, i.e. compression, tension and shear tests have been conducted at 450, 289.2 and 290 mT magnetic fields, respectively. We use the compression data carried at 450 mT to identify the magnetic parameters and then to predict behavior of the material in tension and pure shear at 289.2 mT and 290 mT fields. This way the capability of the proposed model can be examined in both deformation mode change and field intensity variation. Moreover, the compression test includes a wider range of data which makes the parameters identification process more efficient. Figure 4 shows the model results for magneto-mechanical loading in comparison with the experimental data. Comparing the results with figure 3 it is observed that, application of the magnetic field has led to hardening of the material where this effect is more evident for samples with higher particle content. It is also concluded that the model well predicts the magneto-mechanical behavior of the material at different deformation modes and magnetic fields. The model exhibits an outstanding accuracy in the shear tests and quite acceptable results in the tension mode whilst the authors of [68] had doubts about the accuracy of their data in tensile mode because of relatively large distances covered in tensile tests which complicates maintaining a uniform field during the tests.
The effect of magnetic field on the (small strain) stiffness of the samples is shown in figure 5. It is observed that the effect of magnetic field is more pronounced for samples with higher particle content and also a higher magnetic field is required for saturation of these samples.
In terms of the lateral deformations, since the magnetic field is applied parallel to the mechanical loading, the transverse stretches would be equal λ 1 = λ 2 = λ −1/ 2 3 It should be noted that, if the field was perpendicular to the loading direction, the transverse stretches could have been different which will be shown later.  Unlike the isotropic samples, the data for anisotropic samples are only reported for 10% and 20% samples but mechanical tests have been conducted in two loading directions; (1) in the preferred direction, and (2) perpendicular to the preferred direction of the samples. Due to the wider range of the data in the compression test, we again use the compression data for material parameters identification. In this case, as a common practice for anisotropic materials, the compression data in both loading directions are used simultaneously to determine the material parameters. Using the algorithm proposed in [73] assures uniqueness of the identified parameters set. The material parameters are reported in table 2 and results of the model are shown in figures 6(a) and (c) in comparison with the compression test data. The small discrepancies are due to the fact that the curves of both directions have been fitted simultaneously rather than independently. Predictions of the model for tensile and shear modes are reported in figures 6(b) and (d) and figure 7, respectively. It is observed that though the parameters identification has been done only by means of the compression data, the model predictions for both tension and shear are good.
Magneto-mechanical loading data are also provided in [68] for two different configurations; (1) magnetic field parallel to the preferred direction of the material ( N ∥ H), and (2) magnetic field perpendicular to the preferred direction ( N ⊥ H). In both cases the mechanical load is parallel to the magnetic field.
Keeping the same parameters identification procedure previously conducted for isotropic samples, we use the compression magneto-mechanical data to determine the magnetic parameters given that we have already determined the mechanical parameters from pure mechanical tests. Figure 8 shows the reproduced compression curves at 450 mT as well as predicted  tensile curves at 289.2 mT in comparison with the experimental data. It is observed that for (H ⊥ N) the model predictions are in good agreement with the experiments while for the case (H ∥ N) the initial slope of the experimental tension curves are surprisingly high which is underestimated by the model. It is thought to be due to the non-uniform magnetic field during the tensile tests as the authors of [68] have reported that the distribution of the field across the samples was in a range of 96.5%-160.7% of the nominal field intensity. This means that the samples might have experienced up to 60% stronger field than what has been reported and accordingly considered in the simulations. It should also be noted that, the theoretical stress induced by magnetic interaction is generally defined based on an ideal structure which requires perfect alignment of the particles within the material, and maintaining a constant mean distance between adjacent particles while the MRE is being deformed. Obviously this may not be true in practice, though. Furthermore, the magnetic interaction between particle chains are often neglected in simulations for the sake of simplicity. This simplifications can always result in some discrepancies between the modelling results and real test data.
In terms of the lateral deformations, for the case (H ∥ N) as the loading direction is in the symmetry axis of the sample, the transverse stretches are expected to be equal as happened for the isotropic samples. However, it is not the case for (H ⊥ N). Figure 9 shows the transverse stretches induced in the samples in the compression test for the case (H ⊥ N). As the material is loaded in its plane of isotropy, the transverse stretches are expected to be different while the incompressibility constraint is still fulfilled. The results show that, the transverse expansion is lower in the direction of chains which can be explained by  the mutual attraction of the adjacent particles under the effect of magnetic field. The issue is more pronounced for samples with higher particle content which causes a stronger mutual attraction between the particles. Figure 10 shows the predictions of the model for shear tests at h = 290 mT. As shown, the results are in good agreement with the experimental data for (H ⊥ N) but the initial slopes are again underestimated for (H ∥ N).
Verifying the effectiveness of the model in different configurations of magnetic field and loading with respect to the preferred direction of the samples, the model is now employed to predict the magneto-mechanical behavior of the anisotropic samples for the case the magnetic field is applied perpendicular to the loading direction. The results are shown in figure 11 both for (H ⊥ N) and (H ∥ N) where it can be observed that the material behaves stiffer in the latter case. This can be  explained by the more pronounced effect of magnetic field on aligning the chains and mutual attraction of particles when it is applied in the chains direction. Comparing figure 8 with figure 11 one concludes that the hardening effect of magnetic field is more pronounced when it is applied in the preferred direction of MREs and the loading is also in the same direction.
As mentioned before, the transverse stretches for the configuration shown in figure 11(a) would almost be equal but for the configuration mentioned in figure 11(b) the associated transverse stretches are shown in figure 12. In this case the transverse expansion in the direction of magnetic field is significantly lower than the other direction which proves the stiffening effect of the magnetic field in the chains direction. To have  As shown in the above figures and also depicted in [68], it is expected to have the largest effect of magnetic field in case the anisotropic samples are loaded as their particles aligned in the direction of the magnetic field. This issue is well reflected in the introduced parameter A * which takes into account the angle θ between N and M. Figure 14 shows the magnitude of A * for different respective configurations of the material and the applied field. It is shown that the effect of magnetic field is maximum for (H ∥ N) while it gradually decreases as θ tends to the right angle. Even for the case (H ⊥ N) there might be a small effect in practice as the particles are not often aligned ideally in the samples. Schubert and Harrison reported that the anisotropic samples with particle chains perpendicular to the loading direction exhibited very little MR effect [68] which is in correlation with the predictions of the model according to the above explanations.

Experiments reported by Khanouki et al
To further assess the model, we also examine the data reported by Khanouki et al [72] on some isotropic MRE samples. The samples are consisted of different volume fractions of CIPs (BASF SQ ® ) embedded in Ecoflex TM silicone rubber. The investigations includes mechanical as well as magnetomechanical simple shear tests conducted in several magnetic field intensities.
In order to calibrate the model for this particular MREs, the mechanical test data are used to identify the mechanical parameters in the first step. Then, an arbitrary magneto-mechanical data (e.g., at 0.2 T) is used for identification of the magnetic parameters. The other curves are predicted by the determined parameters which are reported in table 3. Figure 15 shows the model predictions for samples with 5%, 15% and 25% particle volume fractions. It is observed that, the model provides a fairly good representation of the shear behavior under the application of magnetic field. The samples show an almost linear behavior at zero magnetic field     within the entire deformation range but presence of the magnetic field results in a non-linear material response exhibited as a hardening which is more pronounced for samples with higher particle volume fraction. This hardening can be explained by the stronger alignment as well as increased resistance of the magnetized particles against movement in the rubbery matrix.
Khanouki et al [72] also proposed a model based on dipoles interaction called 'Dipole model' to predict shear moduli of the tested samples. Results of the proposed model are compared with the Dipole model in figure 16. For low and moderate particle-filled samples (5% and 15%) both models provide good predictions in the whole range but for the 25% sample the dipole model is no longer reliable at h >0.8 T while the present model is still accurate in the full range. It is concluded that, for high particle-filled samples the dipole model fails at higher fields while the proposed model works satisfactorily in almost the entire test range.

Summary and conclusions
In the present work, a finite strain continuum model was proposed for MREs where the effect of magnetization on stiffness was directly introduced in the material shear modulus. It was shown that this approach simplifies the constitutive models and also takes into account the magnetic saturation of MREs. The coupled effects of magnetization, deformation and particle-chains orientation on the mechanical behavior were also reflected in the same new parameter. It reduces the number of material parameters, required experimental tests for parameters identification and also simplifies the mathematical form and numerical formulation of the developed constitutive equations. A two-step method was also introduced for material parameters identification. The predictive capabilities of the proposed model were examined through mechanical and magneto-mechanical experimental data on both isotropic and anisotropic MREs at different configurations of magnetic field and loading with respect to preferred direction of the samples. It was observed that the model satisfactorily predicts the characteristic magneto-mechanical behavior of both isotropic and anisotropic MREs at different deformation modes and loading conditions by using a small number of invariants and material parameters.