Magnetorheological elastomer composites: Modeling and dynamic finite element analysis

Abstract Magnetorheological elastomers (MREs) are polymers reinforced by ferromagnetic particles that show magnetic dependent behavior. Mixing MREs with reinforcing fibers can create a new class of material so-called “MRE composites, MRECs” with additional functionalities and properties. Here, using a Generalized Maxwell model, we proposed a new magnetic-dependent rheological model by considering the hysteresis phenomenon for MREs to predict the dynamic damping responses of MREC plates reinforced by fibers in the frequency domain. We also investigated the influence of magnetic flux intensity, the volume fraction of the fiber, the orientation angle of the fibers, the number of layers, as well as the fiber-to-matrix stiffness ratio on the natural frequency, loss factor, and mode shapes of MRECs plates. Our results suggest that homogenously increasing the elastic properties of the MRECs through the spatial distribution of fibers and changing the fiber-to-matrix stiffness ratio can effectively tailor the dynamic properties of MRECs. Tailoring these properties can provide additional freedom for the fabrication of 4D-printed MRE-based composites.


Introduction
Composite materials with superior material properties and functionalities are the results of the precise placement of their constituents, namely fiber and matrix [1,2]. Examples are carbon fiber reinforced concretes with high tensile properties (i.e., low weight-to-strength ratio with remarkable stiffness properties, high flexural strength or toughness [3]), functionally graded composites with improved interfacial bending strength and their stability investigations [4] and high thermomechanical behavior [5], multi-layer composites including smart cores with excellent dynamic properties [6], and carbon fiber composites with damping enhancement [7].
The advanced composite materials can be fabricated from smart materials such as shape memory alloys (SMAs), shape memory polymers (SMPs), piezoelectrics, magneto-electro-elastics (MEE), electrorheological (ER) and magnetorheological (MR) materials. Using these materials can provide additional functionalities to the composites, such as shape memory effect, reversible cyclic behavior, and mag-netic or electric dependent behavior [8][9][10][11][12][13][14]67]. Among these smart materials, MR elastomers (MREs), which are polymers reinforced by ferromagnetic particles [15], have recently gained a great deal of attention due to their sensitivity to the magnetic fields [16]. MREs can be used as base materials in 3D printing processes [17], to create structures with piezoresistivity [18], positive piezoconductivity [19], and adaptive mounting [20] capabilities. Having ferromagnetic particles, MREs can be polarized in specific directions parallel to the longitudinal direction of the magnetic field, which can eventually lead to a hardening phenomenon [15]. Therefore, MREs can exhibit magnetoelastic interactions in the presence of damping behaviors. The dependency of MREs on the magnetic field makes them useful to control the properties of materials by merely changing the external stimuli (e.g., the amplitude or direction of the magnetic field). Being activated by an external stimulus makes MREs an appropriate candidate for 4D printing where a shapeshifting of a 3D printed structure under external stimuli overtime is needed [21][22][23][24]. Besides, tuning the elastic properties of MRE composites (MRECs) provides an excellent opportunity to create materials with func-tionally graded stiffnesses useful for several high-tech industries such as biomedical or automotive industries.
MREs were first used to test the dynamic characteristics of tunable adaptive vibration absorbers [25]. The damping properties of polymerbased composites in general and MRECs, in particular, can be influenced by dynamic losses. The fabrication process and consequently the damping response of MRECs can be optimized with respect to the influences of various parameters, including the magnetic flux intensity [26], the volume fraction of the magnetic particles, the orientation of the particles, and material properties of the matrix (e.g., using natural rubber as the matrix [27]). The damping phenomenon can be started from the matrix, fiber, or from the interphase between the matrix and fiber [28]. Therefore, the spatial distribution and direction of magnetic particles (i.e., isotropic or anisotropic orientations) in the MRECs can also tune the damping properties [29]. For example, it has been shown that anisotropic MREs exhibit higher modulus degradation than those with isotropic particle distribution [27,29].
In the past, several computational and experimental studies have been performed on the vibrational [30], damping [31,32], dynamic stability [33,34], non-linear static [35], and torsional dynamic [36] responses and viscoelastic behavior [37] of composites with an MRE core. In addition, different phenomenological models have been proposed to simulate the magneto-mechanical characteristics of MREs [38][39][40][41][42][43][44] where various theoretical models for predicting the magnetic field-dependent mechanical properties were used. Besides, MREs have applications in the fabrication of isolator devices where the effects of hysteresis are dominated. Several models have been proposed for such applications in the past [45][46][47] that are able to predict the hysteresis behaviors. These models, however, did not take into consideration the magnetic-dependent properties of the smart elastomers (i.e., MRE), and they did not present a general relationship of magneticdependent parameters such as the model presented in [43]. In addition, some available models in the literature are complicated and they cannot be calibrated, and implemented in commercial FEM software easily such as those presented in [48,49]. The development of such models requires the implementation of advanced constitutive models. This enables us to consider several multiphysics aspects into a model simultaneously. These advanced models also help to properly analyze and capture the response of the MREs under more complex loading scenarios. Furthermore, the dynamic analysis of MREs-based multilayer composite plates has not been investigated before.
Here, using a four-parameter viscoelastic model (i.e., Generalized Maxwell model with two branches and an equilibrium branch), a new constitutive equation for MREs was presented and its application in predicting the dynamic responses of MREC plates in the frequency domain was analyzed. Based on the proposed model, the storage and loss modulus of MREs in terms of magnetic flux intensity were derived.
We found the parameters of the proposed model in terms of the magnetic flux intensity obtained by a nonlinear regularization technique fitted on the experimental DMA reported in the literature [26]. Moreover, the hysteresis behavior of the MRE composites in different frequencies and magnetic flux intensities were investigated within the framework of the developed model. Then, based on the the firstorder shear deformation theory (FSDT), the formulation of the free vibrational analysis of MRE-based composite plates derived and solved using finite elemenet (FE) approach. Also, the verification of the problem was performed by comparison of FE code and those of commercial software (i.e., ABAQUS). We also parametrically analyzed the effect of the magnetic field magnitude, fiber orientation and volume fraction, and elastic properties of fibers on the free vibration behavior of the MREC plates.

Problem definition
In this study, we assumed thin MRE-based multi-layer composite plates. Therefore, to derive their constitutive equations, we used the equivalent single layer (ESL) approach combined with the FSDT of plates [50][51][52][53].

Governing equations of multi-layer MREC plates
We assumed that the multi-layer composite plate consisted of n layers (Fig. 1a). Also, because of the negligible thickness of the plate in comparison with its length and width, the plane stress assumption was assumed. To drive the equation of motions, the displacement field based on FSDT, with three degrees of translations (i.e., u 0 , v 0, and w 0 which are displacements of the mid-plane in the directions of x, y, and z) and two degrees of rotations (i.e., φ x and φ y which refer to the rotations with respect to y and x directions) were used. Therefore, the displacement field in the global coordinate system (x, y, z) can be presented as [54,55]: Here, we assumed infinitesimal strains to derive the governing equations. Therefore, the infinitesimal strain tensor can be expressed as [56]: By substituting the displacement field of the FSDT (i.e., Eq. (1)), in the strain tensor (Eq. (3)), we have: where, E S1 ,E S2, and E S3 are strain components where E S1 þ zE S2 represents the in-plane strains (i.e., ɛ xx ; ɛ yy ; γ xy ) and E S3 represents the transverse shear strains (i.e., γ xz ; γ yz ) and can be calculated as follows: ; The strain tensor for any element can be re-written in terms of the displacements once inserting Eq. (5) in Eq. (3). The strain tensor and displacement field vector are related to each other by means of the strain interpolation matrix of the elements of the media (known as B-matrix). Similar to strain components, the matrix B can be divided into three sub-matrices, namely B S1 ,B S2 and B S3 . The consistent submatrices can be expressed by: B S1 ¼ B S1i B S1j ::: ½ ¼ In the above equation, i, j subscripts denote the node numbers. N i corresponds to the shape function of the i th node of each of the elements. In this study, the bilinear quadrilateral Q4 elements with the Lagrange interpolation functions were used. Hence, the shape function of the elements mentioned above can be expressed in the following form [57]: in which ξ and η are the local coordinates of the Q4 elements. We assumed that the composite is made of MRE, which its matrix shows viscoelastic behaviors. The viscoelastic behavior was calculated using the complex modulus of the composite to include the loss storage moduli. The complex elastic and shear moduli of viscoelastic orthotropic composites can be expressed as [58]: where i denotes imaginary unit. Also, subscripts 1, 2, and 3 represent the longitudinal, transverse in-plane and transverse out-of-plane directions, respectively. The real part of each complex moduli is known as storage moduli, while its imaginary part shows loss moduli. Moreover, E and G indicate elastic modulus and shear modulus, respectively.
The rigidity matrix, D, which indicates the characteristics of the composite material, can be defined as [51]: in which, where, κ is the shear correction factor where depends on the geometry and in this paper becuause of assuming a homogeneous section, κ ¼ 5=6 is considered [2]. The equivalent stiffness tensor (i.e., ðQ À Ã ij Þ k ) of the k th layer in the fibers' direction which can be obtained as follows: in which, T is the transformation tensor and θ is the angle between the fiber and x-axis. On the other hand, using the principle of the minimum potential energy, the elementary mass matrix in the local coordinate can be written as [57]: where, ρ is the density matrix of layers, and det(J) or J (Jacobian) is the determinant of the Jacobi matrix that maps the local coordinate to the general one. For the Q4 elements, the Jacobian and parameter N can be defined as: The stiffness matrix of the elements of the composite plate can be defined by separation matrix operation and dividing the elementary stiffness matrix into five dependent sub-matrices as follows: where k S1S1 , k S1S2 or k S2S1 , k S2S2 and k S3S3 are respectively membrane, coupled membrane and bending, bending and shear stiffness matrices. The sub-matrices introduced in Eq. (14) can be defined as below [57]: Also, k S2S1 equals the transpose of the matrix k S1S1 (i.e., (k S2S1 ) ij = (k S1S2 ) ji ).
The integrations over the transformed area (i.e., Eq. (15)), were computed by means of the well-known numerical Gaussian integration technique. The integral, which corresponds to the shear stiffness of the composite, k S3S3 , was computed using the reduced integration to prevent shear self-locking in the elements, whereas other integrals were numerically computed.
Finally, the damped dynamic behavior of a free vibration problem of a viscoelastic composite plate can be expressed as follows: where, K R and K I are the real and imaginary parts of the total element stiffness matrix, respectively. Besides, M is the overall mass matrix of the system, and ω* is the complex eigenfrequency, while U* is the corresponding mode shape vector. Solving the above eigenvalue problem results in the natural frequency (ω) and loss factor (η) of the laminate in the following form [59]: where Re(.) and Im(.) operators show the real and imaginary parts of the desired complex argument, respectively.

Derivation of the mechanical behavior of the fiber and MRE
To obtain the constitutive equations of the MRECs, firstly, we showed the dependency of the material properties on the magnetic field. The novel model proposed here is developed based on the experimental data reported in [26]. We included the time-dependency of the viscoelastic properties using a Generalized Maxwell model with two branches and an equilibrium branch [58]. Finally, the equivalent constitutive equations were obtained by coupling the recently achieved relations of the MRE and fiber-reinforced composites (FRCs) using the micromechanical approaches available in the composites' literature [60][61][62].

Constitutive equation of MRE
The 1-D differential constitutive equation of the rheological model presented in Fig. 1b can be written as [58]: where, For a given small uniaxial sinusoidal strain input expressed as ɛ ¼ ɛ 0 e iωt the stress output can be obtained as σ t ð Þ ¼ σ Ã e iωt . The term σ Ã is the complex stress and can be assumed to be in the following form [58]: in which E* (or G*) is the complex modulus. Decomposition of the complex modulus, which leads to reaching two real and imaginary modules, can be written as [58]: Substituting strain input and stress output in Eq. (18) results in the following relations for the storage and loss moduli of the MRE Based on the Generalized Maxwell model used here, one can derive the simplified forms of storage and loss moduli of the MRE as follow: We aimed here to include the effects of the local magnetic field on the storage and loss moduli of the MRE. Thus, the following functions are developed to cover that effect: In order to determine the a i 's (i = 1,…,15), we fit our model to the experimental data reported in [26] (Fig. 2, Table 1) where was obtained from a DMA test in the constant frequency of 5 Hz for a natural rubber-based MRE with a Root Mean Square Error (RMSE) of 0.09 and 0.03 for storage modulus and loss factor, respectively. Toward this aim, an optimization technique (or a nonlinear regularization technique) based on the nonlinear least-squares of the objective function was used. Following this procedure, the best coefficients that can satisfy the identity between the recommended functions and those achieved from the experimental tests were developed. The following objective function based on the least-square principle with weight parameters of w 1 and w 2 was defined for the optimization algorithm: The subscripts "Model" and "Exp" show the result obtained from the present model and experimental data, respectively. Besides, the hysteric behavior of the present model was investigated. To this purpose, by substituting Eqs. (19) and (24) in Eq. (18), and considering a harmonic strain as an input (i.e., ɛ ¼ ɛ 0 sin 2πft ð Þ), the output stress was calculated. As a result, the hysteresis loop at different frequencies and magnetic flux intensities was measured (see Fig. 2c and d). It is noted that ɛ 0 was the amplitude of the applied strain, f was the loading frequency and t was time.

Material properties of the fiber
In this study, the glass fibers are implemented as reinforcing fibers in the structure of the MR composite. The linear elastic properties of the glass fibers are listed in Table 2.

Equivalent constitutive equations of the MREC
We used the modified rule of mixture presented by the Halpin-Tsai [61,62] to set a homogenization procedure based on the classical theories of the mechanics of material [60]. Based on this approach, the equivalent density, Poisson's ratio, and moduli of the MREC can be expressed by: In Eq. (26), the orthotropic storage and loss moduli can be defined as below [61,62]: It should be noted that for the elastic fibers, the relation E 00 F ¼ G 00 F ¼ 0 is available. In two above equations, ρ, ϕ, v, E, G and g are density, volume fraction, Poisson's ratio, elastic modulus, shear modulus, and fiber misalignment factor, respectively. The parameter g varies from 0.9 to 1, and in this study, g = 1 is assumed. The subscripts F, M, 1, 2 indicate the fiber, matrix, longitudinal and transverse   Table 2 The mechanical properties of glass fiber based on the values reported in [66]. directions, respectively. Afterward, the superscript 0 and ″ denote the storage and loss terms, respectively. In addition, η and ζ can be defined as following [60][61][62]:

Results and discussions
We used MATLAB (v. R2019a) software to develop our models. The MREC plates were constructed from multi-layer glass fibers and MRE with a dimension of 0.1 × 0.2 × 0.005 m 3 (see Fig. 1a). A clamped-free-free-free boundary condition was assumed for all analyses in this study. It is noteworthy that the whole simulations of the MRE plate in this study were performed at a fixed frequency of 5 Hz.
We validated the accuracy of our model by comparing the shear storage modulus and loss modulus obtained from our model with those of experimental data reported in [26] (Fig. 2). We also compared the natural frequencies and loss factors obtained from our rheological model with those of computationally simulated models. For the computational simulations, we used a commercial finite element code (ABAQUS v.6.17) modeled by 2D-linear shell elements (S4R) and used standard explicit solver (Fig. 3a). These comparisons showed a similar trend among our model, experimental values reported in the literature, and our finite element results.
We used our models to analyze the effects of the magnetic flux intensity, glass fiber volume fraction, fiber orientations, number of layers in composite, and stiffness ratio of fiber to the matrix on the dynamic response (i.e., natural frequencies, loss factors, and mode shapes) of the MRECs.
In what follows, the concentration will be on the free vibrational analysis of MREC plates via a new magnetic-dependent viscoelastic model in the presence of magnetic fields. However, it is worth mentioning that the developed methodology can also be efficiently implemented to extract the transient responses of systems as well as forced oscillation problems.

The effect of magnetic flux intensity
We observed that the storage modulus and natural frequency of all modes of a single-layer MRE plate with a glass fiber volume fractions of 20% and with θ = 0°non-linearly increased by increasing the magnetic flux intensity (Figs. 2 and 3a, b). These parameters reached a plateau after a magnetic flux intensity of 0.6 T (Figs. 2 and 3a, b). This can be due to the fact that increasing the magnetic field intensity increases the elastic properties (i.e., E 1 , E 2, and G 12 ) of the MRECs without affecting its mass properties, which can consequently increase the natural frequency of the composite (Fig. 3a). We observed a similar trend for storage modulus and natural frequency (see Figs. 2a and 3b). As a result, a comparable trend observed for tan (δ) and loss factor.
We also compared the first three-mode shapes of the MREC (Fig. 3c) when being subjected to 0.2 T magnetic flux intensity. The mode shapes were independent of the applied magnetic flux intensity. This can be explained by the fact that increasing the magnetic flux intensity simultaneously increased the elastic moduli of the composite in two orthogonal directions (i.e.,E 1 ; E 2 ), which, as a result, did not influence the overall mode shapes.
Considering the hysteresis effects in our models, we observed a counter-clockwise change of the hysteresis loop when increasing the values of frequency and magnetic flux intensity (Fig. 2c and d). These results suggested that under higher values of frequency and flux iten- Fig. 3. The effects of the magnetic flux intensity on the natural frequency (a) and loss factor (b) of a single-layer MREC with a constant glass fiber volume fraction equal to 20% with an orientation of θ = 0°at a fixed frequency of 5 Hz. The first three mode shapes of the MREC at a magnetic flux intensity of 0.2 T (c). sity, MRE compoistes could become stiffer. This is in agreement with the results of the other studies [39].

The effects of the volume fraction of glass fibers
In composites, the volume fraction of fibers has a significant effect on the behavior of the material. To analyze the impact of glass fiber volume fraction on the dynamic response of MRE composites, we changed the volume fraction of the glass fibers in a single-layer MREC with θ = 0°at a constant magnetic flux intensity (B = 0.5 T).
Increasing the volume fraction of glass fibers can stiffen the MRECs. This non-linearly increased the natural frequency of the plate (Fig. 4a). Contrary to the natural frequencies, the loss factor decreased by increasing the fiber volume fraction (Fig. 4b). This shows, changing the volume fraction of fibers can tune the properties of the MRECs inversely.
The glass fiber volume fraction did not have a significant influence on the first two-mode shapes of the MREC plates ( Fig. 4c and d). The third mode shape, however, changed when MRE plates with different fiber volume fractions (i.e., 0% and 40%) were compared ( Fig. 4c and  d).

The effect of the fiber orientations
We also studied the dynamic response of fiber orientations in a single-layer SMPC plate under a constant magnetic flux intensity of 0.5 T and a glass fiber volume fractions of 20%. When the fibers were aligned with the x-axis (i.e., the orientation angle of θ = 0°), the overall stiffness of the composite plate was maximum. Therefore, the natural frequency was highest in this configuration (Fig. 5a). Changing the orientation of the fibers to 45°, and 90°decreased the natural frequency of the plate (Fig. 5a). The fiber orientation had an inverse effect on the loss factor as fibers with the orientation of 90°showed the highest damping capacity.
The orientation of fibers with 0°, 90°did not change the first twomode shapes of MREC (Fig. 5c-e). However, it affected the first twomode shapes of composite with 45°fiber orientation (Fig. 5c-e, right, and middle subfigures). That is because the variation of the fibers' orientation could change the elastic stiffness of MREC in orthogonal directions. Interestingly, the third mode shape of the composite with different fiber orientation resulted in entirely different mode shapes (Fig. 5c-e, left subfigures).

a.
b.

The effects of the number of layers and fiber orientation in each layer
One way to homogeneously enhance the mechanical properties of the MRE is to make composites from layers with different fiber orientations. Here we considered a two-layer MREC plate with four different lay-up designs, namely; −45°+45°, 0°-90°, 0°-0°and 90°-90°fi ber orientations at a magnetic flux intensity of 0.5 T and a glass fiber volume fractions of 20%. Each layer had an equal thickness of 0.0025 m.
Such an implementation resulted in a broader range of natural frequencies and loss factors in different modes of deformations ( Fig. 6a  and b). For example, the natural frequency and loss factor of the third mode of deformation of the composite with 0°-90°layers reached its maximum values ( Fig. 6a and b). This can be used as an alternative design strategy to tailor the dynamic properties (i.e., damping capacity and frequency) of MRECs.
We also compared the mode shapes of bi-layer composites with fiber orientation of 0°0°and 45°45°with those bi-layer composites with + 45°, −45° (Fig. 6c-e). Changing the orientation of the fibers in each layer changed the amplitude of the deformation of the first mode ( Fig. 6c-e, left subfigures) and shape of deformations of the second and third mode of deformations (Fig. 6c-e, middle and right subfigures).

The influence of fiber to the matrix stiffness ratio
We varied the stiffness ratio of fiber to the matrix (i.e., E F =E M = 1-1000) while assuming a fixed volume fraction of fibers (5%) and magnetic flux intensity (0.5 T). Increasing the E F =E M non-linearly increased the natural frequency (Fig. 7a) and inversely decreased the loss factor (Fig. 7b) of the first and second deformation modes of MRECs. That is because increasing the stiffness ratio will enhance a. b. c.
d. the overall (or equivalent) stiffness of the MREC. The same trend was observed for the third mode of deformation until when 2.5 order of magnitude stiffer properties for the fiber was chosen (Fig. 7).

Conclusion
In the present study, we proposed a new magnetorheological model (i.e., a magnetic-dependent Generalized Maxwell rheological model) for MREs by considering their hysteresis behavior. We presented the material parameters in terms of magnetic flux intensity explicitly unlike, other models [43] with a straightforward approach in calibrating and implementing the present model in addition to what proposed in other literature [48,49]. Then, we used the proposed model for single-and multi-layer MREC plates and evaluated the effects of the composite properties and magnetic field on their dynamic responses in the frequency domain.
MRECs showed strength enhancement when being exposed to a magnetic field with higher magnetic flux intensity. Therefore, by merely changing the magnitude of the magnetic field, one can effectively tune the natural frequency, loss factors, and mode shapes of the composite. These properties can also be tailored by a rational distribution of the fiber orientations, fiber volume fraction, and fiber-tomatrix stiffness ratios. We summarized some of the highlights of the present study as follows: -Increasing the magnetic field could lead to an increase in the storage modulus of the MRE plates. This could consequently increase the values of natural frequencies similar to storage modulus and changing loss factors similar to tan (δ).

a. b.
c. d.
e. Fig. 6. The effects of glass fiber orientation on the natural frequency (a) and loss factor (b) of bi-layer MREC with a glass fiber volume fraction of 20% at a magnetic flux intensity of 0.5 T and at a fixed frequency of 5 Hz. The first three mode shapes of MREC plates for cases with a bi-layer MREC with fiber orientations of 0°0°(c), 45°45°(d) and + 45°−45°(e).
-The volume fraction of the glass fibers in the MRE composite could effectively tune its damping behavior. This could intensify the values of the natural frequency and lessen the loss factor of the MRE plate. -Increasing the fiber orientations from 0°to 90°resulted in a softer MREC, and eventually, decreased the MRE's natural frequencies.
Also, the loss factor increased because of its inverse relationship with the natural frequency. -The number of layers and the orientation of the fibers in each ply could significantly change the natural frequency, loss factor, and mode shapes of the MRE plate. -Increasing the fiber-to-matrix stiffness ratio caused a nonlinear increase in the natural frequency corresponding to the decrease of the loss factor of the MRE plate.
These properties can be used in the rational design of 4D printed structures where the specific shape of deformation overtime is expected. Such MRE materials can also provide more freedom for the designer to create materials with new functionalities and properties.
Also, it should be pointed out that the proposed model can be used once it is aimed to control the amplitude of any arbitrary system's fluctuation. Indeed, the magnetically responsive constitutive behaviors of such smart composites empower the designer to make the best gain from the tunable rheological features of such materials to control the vibration amplitude of the system without using any other controller, as reported in the open literature [63][64][65]. Fig. 7. The effects of the fiber to matrix stiffness ratio (i.e., E F =E M ) on the natural frequency (a) and loss factor (b) of a single-layer MREC with a glass fiber volume fraction of 5% at a magnetic flux intensity of 0.5 T and at a fixed frequency of 5 Hz.