Numerical Modeling of the Mechanical Characteristics of the Magnetorheological Elastomers

. The paper presents the work results associated with the implementation of the developed mathematical model describing the magnetorheological elastomers mechanical characteristics under compressive stress. The employed mathematical instrument constitutes the classic viscoelastic model development, accounting in its structure for the effect of the strain amplitude together with the magnetic field induction on the measured stress values progressive characteristic. The necessity to introduce such changes was demonstrated as a result of earlier studies. Within the this work scope, the structure of the developed numerical model is presented. Further provided are the results of carried out simulations, which are finally juxtaposed and compared with the earlier experiments results. This gives ground to evaluating the proposed model accuracy.


Introduction
The term "Smart Materials" (adaptive materials, intelligent materials) is understood as a broad group of materials characterized by the change of their physical characteristics (e.g. viscosity) when exposed to another field (e.g. electromagnetic field) [1]. Their capabilities include, among others: detecting specific stimuli, processing the stimuli, typical reaction and return to initial state, all in the shortest possible time. The object of the current study are magnetorheological elastomers. These are classified as Smart Magnetic Materials (SMM). It is a relatively new material group undergoing dynamic development. Current studies by different centers focus on the manufacturing methods, material composition and structural formation of magnetorheological elastomers. Significant attention is directed to their practical applications. After year 2000 [1], there is an observed major increase in patent activity in regards to the material itself, but also its application. An analysis of solutions provided in subject literature allows to determine that the major area of use for magnetorheological elastomers focuses primarily on vibration dampers [2][3][4][5]. This application includes not only active control of linear vibration which use the discussed materials in compression [6] or shear [7,8] modes of action, but also active control of torsional vibration [9,10]. Many works were taken up involving the issue of mathematical representation of magnetorheological elastomers magnetic and mechanical characteristics. Despite this, no model to describe the behavior of magnetorheological elastomers was offered to date that would be deemed complete and finished. It is therefore justified to continue the research effort in this regard. The present paper demonstrates the ongoing research results into the material characteristics and the possible application of magnetorheological elastomers in the broadly understood machine building [11,2]..

Development of the mathematical model
As demonstrated earlier, the classic Kelvin-Voigt rheological model (viscoelastic material) cannot be deemed sufficient for the characteristics modeling of magnetorheological elastomers in a broad strain amplitude and excitation frequency range [11]. Considering the above, certain modifications were proposed to allow to account for the strain amplitude effect on the set of obtained results. One needs to emphasize that the variance of mechanical characteristics caused by excitation frequency, in the general case, cannot be omitted. The problems associated with the examination of this dependence are the subject matter of planned further studies. The developed mathematical model can be expressed as follows [11]: where: ( )the time course of a stress function, ( )the time course of a strain function, ( ) stands for strain rate, is the modified stiffness coefficient, represents the coefficient of viscosity. Considering that excitation shall be a known sinusoidal waveform of strain with respect to time, the formula (1) can be expressed as: where: and are the parameters for the hysteresis loop shape, whereas stands for the stiffness coefficient (classic model for viscoelastic body). The respective values in model (2) can be also characterized with the following dependence (3): where ′ stands for storage modulus, whereas ′′ represents the loss modulus. In the presented model, the viscosity coefficient is a linear function of magnetic field induction . The modified stiffness coefficient is the function of strain, the hysteresis loop shape parameters and as well as the magnetic field induction (indirectly via the stiffness coefficient ). So as to increase the representation accuracy of actual results, it was decided to utilize a higher number of hysteresis loop shape parameters than two ( and ). The results analysis of tests carried out to date [12][13][14] allows to determine that the highest variance in the registered hysteresis loops characterizes the material loading (upper part of the hysteresis loop) from the half of the value of strain amplitude. As regards the lower part of the characteristic, being responsible for unloading, one observes similar behavior of the graph lines; however they are occurring for different strain values. It was determined that the range in which the variety of the characteristic is the highest occurs from the maximum strain value 0 up to the value equal to 0.69 0 . It was therefore decided that this value shall constitute the threshold in the variability range of hysteresis loop shape parameters. On this basis, eight parameters of hysteresis loop shape were selected ( 11,12 , 21,22 oraz 11,12 , 21,22 ). A schematic representation of the graphical interpretation and hysteresis loop shape parameters assumption method is provided on Fig. 1. As a consequence of assuming a higher number of parameters, the expression (2)

Determination of model parameter values
Numerical values determination for the modified Kelvin-Voigt model parameters required carrying out a number of calculations. To this end, the authors proposed an algorithm of their own development allowing to calculate the values of parameters and as a magnetic field function, based on earlier measurements, as discussed in [12,13]. The diagram for hysteresis loop shape parameters determination is provided on Fig. 2 below. ). In the following stage, the parameter values for and are re-assumed with low scatter, close to the values which were defined in the previous step. This process aims to determine the sought parameters exact value. Afterwards, the data set for 1 and 2 , is determined once again and the algorithm searches for the lowest values to indicate the final values for 11,12 , 21,22 and 11,12 , 21,22 . In the earlier studies [12,13], the stress-strain curves (hysteresis loops) for magnetorheological elastomer samples were registered at different strain amplitude 0 values and magnetic induction . These were used as a basis for the parameters determination for the described mathematical model. First of all, the averaged moduli values were determined: the storage modulus ' and loss modulus " in the Kelvin-Voigt model (Fig. 3). Finally, the values for the hysteresis loop shape parameters 11,12 , 21,22 and 11,12 , 21,22 were determinedin accordance to the presented methodology ( Fig. 4 and  5).

Model implementation in the environment for numerical calculations
MATLab Simulink software package components were usedin the development of the model as defined in the previous chapter. The software enables carrying out simulations and calculations in real time. The implemented model view in this environment as well as its structure comprising three components defining the characteristics of the examined composite as appropriate, are provided in the schematic view on Fig. 6. The model component marked with 1 is responsible for modeling the viscous characteristics. These are dependent on the loss modulus '', which is described by the expression (3). This function is a linear dependence of the magnetic field induction , in accordance with the dependencies provided on Fig.3. The subsequent model's components, marked with numbers from 2 to 5, determine: the upper and lower part of the hysteresis loop. At the same time, functions for the storage modulus ' and hysteresis loop shape parameters 11,12 , 21,22 and 11,12 , 21,22 are used, with linear dependence on the magnetic field induction , in accordance with the dependencies provided on Fig. 4 and 5.    An analysis of the provided graphs shows that the biggest errors in stress value determination appear in the extreme strain values vicinity. A quality evaluation of the modified numerical model was therefore attempted. In order to determine the error value, it is necessary to utilize the measured stress value in a given time limit ( ), as well as the stress value in a given time limit from the model ( ). The deviation of the modified numerical model in the given time limit is defined as below (5): The numerical model allows to set any time of progression, allowing to determine the stress value at any time limit . This was used to determine the error value as a function of time for individual points of measurement in the entire strain range. The scheme representing the error value determination of the numerical model ( ) is provided on Fig. 9. The selected determination results of the model error value are provided on Fig. 10.  As follows from the analysis of the presented graphs, the agreement between the numerical model and the measured data is unsatisfactory for very low strain values. This stems in part from the assumptions regarding the linear influence of the magnetic induction on the model parameters, and in part from the imperfections in the Kelvin-Voigt model. The result of modifications introduced to the model (implemented hysteresis loop shape parameters: 11 , 12 , 21 , 22 oraz 11 , 12 , 21 , 22 ) brings about an improved consistency between the results of numerical analysis and actual measured results for higher strain values. The range and character of changes observed in the characteristics received from the numerical model are satisfactory in comparison to the observed behavior of the magnetorheological elastomer in the course of the experimental study. It should be pointed out that in a certain selected range (excluding low strain values), the deviation of the proposed model oscillates in a satisfactory error raterange. This is shown on Fig. 11 and 12.   Fig. 11. Error rate ( ) for the developed numerical model, excluding low strain values for frequency = 0.55 Hz, strain amplitude = 5% and magnetic field induction = 64 mT; A -loading, B -unloading.

Summary
A detailed variance analysis shown by the characteristics obtained in actual measurements and characteristics obtained in numerical calculations allow to ascertain that the model offers an adequate representation of actual values only above a specific value of strain amplitude. Higher discrepancy is observed, the more the strain value in the given time limit differs from its maximum value. The assumed model allows to describe the loop in the shape of an ellipsis. Therefore, certain differences between hysteresis loops determined numerically and actual hysteresis loops cannot be avoided. It should be concluded that the influence of the increasing cross-section of the examined samples is highly important. The variance range in characteristics of magnetorheological elastomer resulting from this phenomenon exceeds the variance range effected by the magnetic field. This is highly disadvantageous as the effects of these phenomena tend to overlap, increasing the difficulty of mathematical modeling. Taking into consideration the strain amplitude and frequency of excitation in a broader range than was proposed calls for performing additional studies focusing on the exact determination of the character of changes occurring as a result of significant crosswise strain (change in the cross-section of the examined samples of magnetorheological elastomers).