Material modeling of frequency, magnetic field and strain dependent response of magnetorheological elastomer

Accurate modeling of material behavior is very critical for the success of magnetorheological elastomer-based semi-active control device. The material property of magnetorheological elastomer is sensitive to the frequency, magnetic field and the input strain. Additionally, these properties are unique for a particular combination of matrix and the filler loading. An experimental-based characterization approach is costly and time consuming as it demands a large amount of experimental data. This process can be simplified by adopting material modeling approach. The material modeling of magnetorheological elastomer is an extension of conventional viscoelastic constitutive relations coupled with hysteresis and magnetic field sensitive attributes. In the present study, a mathematical relation to represent the frequency, magnetic field and strain dependent behavior of magnetorheological elastomer is presented. The viscoelastic behavior is represented by a fractional zener element and the magnetic field and strain dependent attributes incorporated in the model by a magnetic spring and linearized Bouc-–Wen element, respectively. The proposed model comprised of a total of eight parameters, which are identified by minimizing the least square error between the model predicted and the experimental response. The variations of each parameter with respect to the operating conditions are represented by a generalized expression. The parameters estimated from the generalized expression are used to assess the ability of the model in describing the dynamic response of magnetorheological elastomer. The proposed model effectively predicted the stiffness characteristics with an accuracy, more than 94.3% and the corresponding accuracy in predicting the damping properties is above 90.1%. This model is capable of fitting the experimental value with a fitness value of more than 93.22%.


Introduction
Decoupling the source of vibration is an important aspect in the design of a structure. An elastomer isolator is an ideal medium to interrupt the vibration path. These isolators are simple and cost-effective; however, its performance is severely effected as the excitation frequency is varied [1,2]. This limitation can overcome by replacing the elastomer resilient element with magnetorheological elastomer (MRE) [3]. MRE is a special category of smart material, exhibit field sensitive rheological properties. The properties of MRE are modified in response to the external magnetic field to meet the requirement of semi-active vibration mitigation devices [4][5][6][7]. The performance of these devices relies on the semi-active control strategies, which subsequently depends on the material models of MRE [8,9].
MRE exhibits the viscoelastic behavior coupled with the field sensitive characteristics. The fractional and integer order-based viscoelastic constitutive relations often employed to model the viscoelastic material behavior of MRE [10][11][12][13]. The magnetosensitive attributes encompassed in the models by incorporating the field-induced dipole interactions between the filers in a periodically spaced dipoles [14,15]. A comparative assessment on the fractional and integer order-based viscoelastic models to represent the frequency and magnetic field dependent behavior were discussed in the past work [16]. This study confirms that the fractional viscoelastic model with 3 point input is effective in portraying the frequency and magnetic field dependent behavior. However, the models based on viscoelastic relations effectively portray the frequency-dependent characteristics but are not fully adequate to describe the strain dependent properties of MRE.
The rheological properties of MRE vary significantly with respect to the input strain. The storage modulus of MRE decrease with an increase in strain [17][18][19][20] and the loss modulus exhibits a reverse trend. This typical phenomenon is referred to as Payne effect. The conventional viscoelastic models fail to concur these responses of MRE [21,22]. Consequently, the control strategies developed based on this framework fail to yield satisfactory results. This situation calls for an extension of the conventional viscoelastic model by incorporating the amplitude dependent effects while designing the MRE-based devices.
A general approach involves connecting a spring-Coulomb elements in parallel with the viscoelastic models. This element represents the stick-slip phenomena resulting from the interface slippage with an increase in amplitude [23,24]. The estimation of the parameters of a frictional element is a cumbersome process as it involves many variables and demands a large amount of experimental data. These limitations can be addressed by Berg's friction models comprising of lesser parameters. Although Berg's model predicts the amplitude dependent characteristics, its accuracy is significantly affected by the selected testing conditions [22,25,26]. This approach was simplified by the modifications proposed by Thaijaroen et al. [27] where the parameter estimation process is carried in two stages. The modeling approach based on an addition of friction element is simple, but its parameter estimation process requires low-frequency and high-amplitude tests, which are difficult to perform in a simple forced vibration test apparatus. Moreover, the stiffening effects induced by the geometrical constraints under high displacement amplitudes will reduce the accuracy of the measured results [27,28]. Alternatively, a new viscotribological model comprising of generalized Maxwell model and Dahl friction model is used to represent the amplitude dependent behavior. However, this model is having limited application due to a large number of parameters, which makes the parameter identification process more complex [29].
Alternatively, the hysteresis models based on Ramberg-Osgood model and Bouc-Wen are used to describe the nonlinear characteristics of MRE [30]. The Ramberg-Osgood model is successfully employed to tune the amplitude dependent characteristics of MRE. However, its effectiveness relies on the process of capturing the transition curve from elastic to the plastic response of a material [31]. This makes the parameter estimation process as path dependent and it has an influence on the identified parameters. On the contrary, Bouc-Wen models could overcome these limitations and are often used to represent the nonlinear characteristics of an isolator system. The advantage of this model is that it does not necessitate the initial path of the target hysteresis loop to estimate the parameters. Compared to other hysteresis models, the Bouc-Wen model comprises of more number of parameters [32], but it can effectively capture the strain stiffening and strain softening behavior often observed in the filled elastomers [33,34]. However, at intermediate strain levels, the distortion in the shape of the hysteresis curve is not significant. In such cases, a modified Bouc-Wen model with stochastic linearized random variable is sufficient to describe the amplitude dependent characteristics [35][36][37]. The simplified models based on linearized Bouc-Wen is reported in many literatures, and it is often used to model the dynamic characteristics of base isolation systems [38,39].
In the present study, the fractional Maxwell MRE model coupled with a stochastic linearized Bouc-Wen component is presented. Based on the model, a mathematical relation is formulated to represent the frequency, magnetic field and amplitude dependent behavior of MRE for the frequency range from 8 to 24 Hz. The model comprises of a total of eight parameters and these parameters are identified by minimizing the error between the experimental and the model predicted response using MATLAB optimization tool box. These parameters are represented in a generalized form to simulate the frequency, magnetic field and amplitude dependent effects simultaneously. The performance of this model is evaluated based on the ability to fit the experimental data, which is expressed in terms of fitness value.

Experimental MRE sample preparation
Cylindrical shaped MRE samples of dimension 10 mm thick and 19 mm in diameter were prepared [40]. Test samples comprised of 27% by volume CIP (diameter 5 lm; BASF; Type CN) and 73% by volume of two-component RTV silicone matrix (base to crosslinker ratio 100:5 from Dow Corning). This matrix-filler ratio represents the optimum level of reinforcement to attain maximum field-induced enhancements [41]. Synthesis process of isotropic MRE included two steps: mixing and curing. Mixing process involved the stirring of particle-matrix blend in a beaker to attain a homogenous mixture. The mixture was later poured into an aluminum mold and degassed in a vacuum chamber. Curing process is carried out at room temperature under constant pressure for 15 h.

Experimental setup
Dynamic compression behavior of MRE was evaluated from the dynamic blocked force measurement [42] for the harmonic input displacement. The schematic representation and the actual image of dynamic compression property measurement experimental setup are shown in Fig. 1a and b. Compressive properties of MRE were tested by the steady-state harmonic excitation from an electrodynamic shaker. The harmonic signals of different frequencies were generated by a NI PXI-5412 function generator and fed to the electrodynamic shaker through a power amplifier system. Force at the blocked end was measured by a force transducer (KISTLER, type 9712). To minimize the effect of rocking motion, two identical accelerometers (KISTLER, K-shear) were used [43]. The sensed signals from accelerometer and force transducer were acquired through NI PXI-4496 data-acquisition system.
The field dependent properties of MRE were tested by employing a pair of rare-earth neodymium permanent magnets (type: NdFeB 32, coercivity: 883,310 A/m, and relative permeability: 1.045). A special type of test rig was designed to alter distance between the magnets to vary the magnetic field. Magnetic flux at the opposite end of the permanent magnet was blocked by a ferromagnetic strip to enhance the magnetic field over the region occupied by MRE. This arrangement can induced a maximum flux density of 0.27 T for the minimum possible separation distance of 30 mm.
Tests were carried out at frequencies: 8 Hz, 10 Hz, 12 Hz, 14 Hz, 16 Hz, 18 Hz, 20 Hz and 24 Hz. Influence of magnetic field was evaluated by carrying out the viscoelastic property measurement studies at different magnetic flux densities (0 T, 0.1 T, 0.2 T and 0.27 T). The frequency and magnetic field dependent properties were measured at constant strain amplitude of 0.25%. The amplitude dependent properties of MRE were evaluated at 0.5%, 1.0% and 1.5% strain. To alleviate the effects of deformation history and stress softening on the measured viscoelastic properties, MRE samples were subjected to preconditioning cycles [44]. Each set of experiments was repeated 10 times, and the average value was considered for the analysis. The stable response of the system was ensured by measuring for 10 cycles in each loading case. To exclude the effect of deformation history on the measured viscoelastic properties, MRE samples were allowed a recovery period of 15 min between successive measurements [45].

Experimental results
Force-displacement hysteresis plots evaluated at different magnetic flux densities (0 T, 0.1 T, 0.2 T and 0. 27 T) corresponding to 8 Hz, 12 Hz, 16 Hz and 20 Hz are depicted in Fig. 2a-d. The slope and area of the hysteresis loop vary with the magnetic field and excitation frequency. Field-induced variation in the hysteresis behavior of MRE is associated with the interaction between the magnetically active fillers [14,46]. The frequency dependency in MRE is attributed to the viscoelastic behavior [42] which inherited from the elastomer matrix.
The amplitude dependent compression behavior of MRE is evaluated by varying the amplitude of harmonic excitation at constant frequency and magnetic field. Figure 3a-d represents the amplitude dependent hysteresis behavior of MRE corresponding to 8 Hz and 16 Hz. As evident from the graph, with an increase in amplitude, the slope of the hysteresis loop decreases and the area of the hysteresis loop increases. This indicates the phenomenon of Payne effect in MRE as reported in our past studies [47].
The dynamic test results revealed the existence of magnetic field, frequency and the strain dependent viscoelastic properties of MRE under compression mode. To describe these effects simultaneously, a mathematical relation to represent material behavior is essential. The material modeling approaches simplify the material property characterization studies by eliminating the costlier and time-consuming data collection processes.

Modelling magnetic field, frequency and amplitude dependent behavior of MRE
The material modelling of MRE based on viscoelastic constitutive relations is effective in portraying the magnetic field and frequency dependent properties. However, these models are not efficient as the input strain is varied. As evident force-displacement plots (Fig. 3), MRE exhibits a pronounced strain dependency like conventional filled rubbers. To address these effects, the fractional Maxwell MRE model developed in our previous work [16] coupled with frictional or hysteresis component is presented.
The proposed model of MRE comprises of a fractional Maxwell element connected in parallel with spring, magnetic element and a Bouc-Wen component (Fig. 4). The fractional Maxwell element describes the damping characteristics, the stiffness and field sensitive characteristics are addressed by elastic and magnetic springs, respectively. The Bouc-Wen element represents the amplitude dependent attributes of the model. The proposed MRE model is one dimensional, and its response depends on the contribution of elastic spring force (F e ), viscous force (F v ), magnetic force (F m ) and the force due to amplitude dependent effects (F A ). The generalized expression for the blocked force (F) experienced by MRE including the magnetic viscoelastic and amplitude dependent effect is given by [48][49][50][51], Elastic force in MRE is due to stiffness of the polymer matrix. The corresponding force associated with the stiffness K sf , for the input harmonic displacement x is given by [48], The coefficients of fractional Maxwell element are denoted by K, c and a . The expression for the force in a fractional Maxwell chain is given by [16,52,53], The relation between the magnetic force F m and the magnetic induction B is expressed as [16], where To address the amplitude dependent response, often nonlinear systems like Bouc-Wen hysteresis model [32] are used.
The expression for the force due to amplitude dependent effect is given by  The parameter K 2 represents the hysteretic component of the restoring force and z is a non-observable parameter described by a nonlinear equation, which is given by [49,54], where A,,, n are parameters that describe the hysteresis loop. Each parameter of the above governing equation influences the hysteresis loops in different ways. This could be in the form of vertical intercept, amplitude, thickness or area under the graph [54] Solving nonlinear Bouc-Wen system to obtain closed-form solutions for the differential equations is a complex task. So, it is a common practice to linearize the system to obtain a frequency response function. Additionally, the linearization in the present study justifiable as the force-displacement hysteresis loops (Fig. 3) revealed that MRE does not show the strain stiffening or softening phenomenon [54]. Moreover, for the tested strain levels, the distortion in the shape of the hysteresis curve is not significant. This further confirms that a modified Bouc-Wen model with stochastic linearized random variable is sufficient to describe the amplitude dependent characteristics.
The Bouc-Wen model can be modified in several ways to simplify the system and to make it computationally faster for parameter estimation. The system    can also be normalized to reduce the number of parameters without any compromise in rigorosity [32,35,54]. Deriving the frequency response function is vital to obtain the storage and the loss moduli of the elastomer. Although the parameter estimation can be accomplished using the nonlinear form of equation, it is not possible to obtain a transfer function for the differential equation. These shortcomings are addressed by a stochastic linearization, which replaces the nonlinear hysteresis component by an equivalent linear model [35,38,39,55,56]. In the linearization technique, the linear model is obtained by taking the reference of nonlinear equation, which portrays the hysteretic loop is assumed as [38], where C e and K e are the linearized parameters, which represent the stochastically equivalent damping and stiffness. Frequency domain transform of the linearized equation is given as follows: The expression for z() is simplified as follows: The overall expression for the blocked force F, with the inclusion of Bouc-Wen stiffness and damping parameters, is given by The expression for the complex form of dynamic stiffness K * (= K R ? iK i ), of proposed model is given by The dynamic stiffness comprises of a real part, denoted by KR and an imaginary part of the stiffness represented by Ki.
Under steady-state, the force-displacement response of MRE for the harmonic input displacement, x=x 0 e ixt is given by where d f is the phase difference between the force and displacement, and the magnitude of dynamic stiffness is given by

Parameter estimation
The proposed model uses the experimental results discussed in ''Experimental results'' Section as the input to estimate the parameters. Using these parameters, the real and imaginary part of the dynamic stiffness are computed. The parameter estimation process involves a nonlinear regression algorithm in order to minimize the error between the experimental and the model-predicted force for each set of experimental results. The expression for the objective function which optimizes the parameter of the model is given by [57,58] J Fði; jÞ À F e ði; jÞ ð Þ 2 ð15Þ where N and M are the number of different experimental data corresponding to the various loading conditions like, frequency and strain. The optimized parameters obtained for the different magnetic flux density (0 T, 0.1 T, 0.2 T and 0.27 T), frequency ( 8 Hz,16 Hz and 24 Hz) and the stain amplitudes (0.05 mm, 0.1 mm and 0.15 mm) are listed in Tables 1, 2 and 3. In total, 36 sets of compression test data were used to estimate the parameters of the model. From the parameters listed in Tables 1, 2 and 3, it is found that the parameters are sensitive to the change in magnetic field, frequency and the amplitude. These variations are represented in a mathematical form by using LINEST function (MS Excel). This function approximates a linear equation to simulate the simultaneous variations of the parameters with respect to frequency, magnetic field and strain. Assumed linear function is advantageous as it offers a simple mathematical formulation of MRE, which is easier to implement in MRE-based semi-active control strategies [59]. A generalized form to represent these variations is expressed as follows: where d is the parameter of the model. The coefficients of the expression are denoted by d 1 , d 2 , d 3 and d 4 . Input frequency, magnetic field and amplitude are denoted by f, B and A, respectively. The coefficients of the generalized expression are listed in Table 4. Using above relations, the real and the imaginary part of the stiffness is calculated to evaluate the ability of the model in reproducing the experimental results.

Comparison between experimental and model response
Using the parameters estimated from the optimization process, a comparison is made between the experimental, and the model predicted behavior. Figure 5 compares the force-displacement hysteresis loops at different magnetic field strengths. Correspondingly, a comparison between the force-displacement relations under different frequency and the input strain are represented in Figs. 6 and 7, respectively. From these graphs, it is observed that the proposed model is able to predict the frequency, magnetic field and strain dependent characteristics of MRE. The area of hysteresis loops is increased by increasing the magnetic field, which is very well portrayed by the model. Additionally, with an increase in frequency, the stiffness of the MRE is enhanced, which is consistent even in the presence of the magnetic field. Moreover, the slopes of the main axis and the area of hysteresis loops vary with the increase in the input strain. This indicates that the proposed model is effective in portraying the Payne effect in MRE under magnetized and non-magnetized state. The real and the imaginary part of the dynamic stiffness are computed using Eq. 13. A comparison between the K R and K i evaluated from the experiments and predicted by the model at different Table 1 Optimized parameters of the model corresponding to 0.05 mm input displacement amplitudes are presented in Figs. 8, 9 and 10. As evident from the graphs, by increasing the magnetic flux density, both K R and K i are increased. This behavior is observed for all the tested frequencies and the strain amplitudes. For example, at a strain of 1.5% (0.15 mm) and frequency of 16 Hz, increasing the magnetic field from 0 to 0.27 T, K R is increased from 168.06 to 196.836 N/mm. The corresponding K R estimated from the experiments is increased from 176.245 to 204.822 N/mm. From the graphs, it is observed that model predicted real part of the stiffness (K R ) is in good agreement with the experimental results. This proposed model is able to reproduce K R with an accuracy of more than 94.3%. However, the accuracy of predicted K i by the model is relatively lower than the K R . Moreover, this model is able to capture the experimentally determined K i with an accuracy greater than 90.1%. This disparity can be attributed to the variation in the damping characteristics of MRE,  Table 3 Optimized parameters of the model corresponding to 0.15 mm input displacement which is not exhibiting a linear variation with respect to the frequency [16]. Moreover, the parameter identification processes have a significant influence on the accuracy of the predicted stiffness values [16]. Apart from graphical analogy, the ability of the model in reproducing the experimental results is quantitatively investigated from the fitness value [57,58]. The fitness value represents the correlation between the force value predicted by the proposed model and the measured force from the experiments. This is calculated by the normalized root-meansquare function, and it is given by the equation: where F P is the predicted force and the F E is the measured force from the experiments. The resulting fitness values are listed in Table 5. In all the experimental cases, the values of the fitness are higher (more than 93.22%). This indicates that the proposed dynamic model of MRE can be effectively used to predict the frequency, magnetic field and strain dependent dynamic compressive characteristics for the tested operating conditions.  Comparison between the experimental data and model predicted K R and K i at different magnetic field corresponding to 0.05 mm input amplitude. transition region for the excitation frequency between 8 and 24 Hz. The model simplifies the overall material characterization process of MRE as it necessitates the experimental response corresponding to 8 Hz, 16 Hz and 24 Hz. This semi-physical modeling approach is beneficial as it reduces the amount of experimental tests substantially and saves time and cost. The proposed model effectively predicted the stiffness characteristics with accuracy more than 94.3% and the corresponding accuracy in predicting the damping characteristic is above 90.1%. Predicted stiffness and damping characteristics by the model depend upon the effectiveness of the parameter optimization process in reproducing the experimental results. The difference arising from the parameter optimization process carried over during the interpolation process, and it reflected on the predictive accuracies of the MRE model. This model is capable of fitting the experimental value with a fitness value of more than 93.22%. Moreover, the variation in the model parameters is represented by a planar function, which simplifies the control strategies to realize the MRE-based device.

Funding
Open access funding provided by Manipal Academy of Higher Education, Manipal.

Declarations
Conflict of interest The authors declare no conflict of interest in preparing this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licen ses/by/4.0/.