# Research on MRD Parametric Model Based on Magic Formula

Authors: Wei Pan; Zhengtao Yan; Jingjun Lou; Shijian Zhu
Journal: Shock and Vibration
 (2018)

---

## Abstract

In order to get a better description to the nonlinear characteristics of magnetorheological dampers, the magic formula is introduced into the general method of parametric modelling of magnetorheological dampers to propose a new parameterized model called magic formula-hysteresis loop model (MFM). The new model is simple in structure, the physical meaning of each parameter is clear, and the parameter identification is convenient. The fitting and experimental data of MFM and the phenomenon model under different conditions are applied for error analysis and comparison. The results show that the errors of MFM are more accurate and have better fitting and experimental data under different working conditions, which also have better adaptability and versatility.


---

## Body

## 1. Introduction

Magnetorheological Fluid (MRF) is a new intelligent material composed of nonmagnetic liquid and tiny soft magnetic particles. The particles are featured by high magnetic flux and low hysteresis [1]. This material exhibits low viscosity Newtonian fluid properties at zero magnetic field, but under high magnetic field, it shows high viscosity and low fluidity Bingham fluid characteristics. Under applied magnetic field, MR fluid’s rheological properties such as viscoelasticity, yield strength, and so on can change rapidly (milliseconds) and reversibly [2]. Magnetorheological damper (MRD), which takes magnetorheological fluid as the actuating medium, has many advantages such as high yield stress, high temperature adaptability, short response time, strong stability, low requirement voltage, and convenient semiactive control, etc [3]. It has drawn more and more attention in the fields of automobile [4], architecture [5], aviation [6], and other fields that require vibration control.By adjusting the damper input current, MRD could change the strength of the magnetic field to control the size of the damping force. Therefore, MRD is often applied in vibration semiactive control. The premise of accurate semiactive control is to establish accurate and reliable MRD model. But the output damping force of MRD, which is influenced by the coupling of external excitation, applied magnetic field, and internal flow field, shows a strong nonlinear character. This results in difficulties of precise MRD modelling. In order to accurately describe the properties of MRD, many MRD models have been established, which are divided into two types: physical model and phenomenological model [7]. Among them, the physical model refers to the use of magnetorheological fluid constitutive model and damper geometry to do the modelling and determine the damping force [8]. Because there are still many assumptions about the mechanism of magnetorheological fluid, the established physical model is not accurate enough. Therefore, the physical model can often be used to guide the damper design but seldom be applied to the control. Phenomenological model is divided into nonparametric model and parametric model. In the nonparametric model, each parameter has no explicit physical meaning and is only used to identify the hysteresis loop that controls damping force, such as fuzzy logic model [9], neural network model [10], and so on. Besides, the parameters in parametric models have comparatively clear physical meaning. The commonly used parametric models are Bingham model [11], nonlinear double-viscosity model [12], Dahl model [13], Bouc–Wen model [14], and phenomenon model (PM) etc [15].Among them, PM, which is an improved version based on the Bouc–Wen Model, is a parameterized model that conforms well to the experimental phenomena. It can accurately describe the MRD low-velocity attenuation, yield characteristic, and hysteresis characteristic, which owns higher precision and incentive adaptability, that is, the model parameters identified by a certain excitation can be applied to a wider range of excitation. Therefore, PM is the most widely used MRD parameter model, and its model expression is listed in the following equation:(1)FPM=c1y˙+k1x−x0,y˙=1c0+c1αz+c0x+k0x−y,z˙=−γx˙−y˙zzn−1−βx˙−y˙zn+Ax˙−y˙.However, PM is composed of two differential equations and some absolute and power function terms, which is difficult to establish and to solve by general numerical method [16]. The model needs to identify up to 12 parameters, of which the existing intermediate variable has ambiguous physical meaning; the relationship between other parameters and external incentives is not clear, so it is difficult to determine the initial value or a fixed range of values, and repeated trial and error and iterative calculations to complete the model parameter identification [17–19] are needed. All of these have brought great difficulties to the establishment and use of MRD model. Therefore, in order to describe the dynamic characteristics of MR fluid damper better, it is necessary to establish an MRD parameterization model which is convenient for modelling and solving and easy to identify parameters with high accuracy and strong universality.The magic formula [20] can get a good approximation to the nonlinear curve through the combination of several trigonometric functions, therefore, is usually used to establish tire mechanics model. The general form of the magic formula is shown in the following equation:(2)F=DsinCarctanB1−Eα+EarctanBα.Among them,B, C, D, and E are parameters which control the curve shape; α is independent variable, and based on the actual situation it can select speed, displacement, deflection angle, and so on. The magic formula can well describe the mechanical characteristics of the tire within the scope of the experimental data. At the same time, the mechanical properties of the tire can still be accurately expressed in a certain range outside the experimental data. That is, the magic formula can predict from the limited working conditions and ensure a high degree of confidence. At the same time, we can see that, as shown in Figure 1, the curve of tire force-deflection angle is very similar to the plastic fluid curve of MR damper; however, the plastic fluid curve cannot describe the hysteresis characteristics of the magnetorheological damper. Considering that the tire mechanics model has similar nonlinear characteristics with the MRD, the Magic Formula is introduced to the MRD modelling process, and a magic formula-hysteresis loop model (MFM) is established. The new model is of simple structure, definite physical parameters and convenient parameters identification. According to the fitting with the experimental data, MFM is compared with PM and error analysis is done. The results prove the accuracy and universality of the new model, which also provide a reference for the application of MRD in vibration control.Figure 1
Tire force-deflection angle.


## 2. Magnetorheological Damper Experiment

The general structure of MRD is shown in Figure2. The magnetorheological fluid is filled in the cavity of MRD. The piston has some fixed damper channels inside, and the coil is wound inside the piston. When an external force acts axially on the piston rod, the piston moves within the MRD cavity, and the magnetorheological fluid flows along the fixed damping channels. When the coil is energized, a magnetic field can be generated, and the properties of the magnetorheological fluid in the damping channels change as the magnetic field changes. The gas compensation chamber at the end of the damper can provide an appropriate compensation force.Figure 2
Magnetorheological damper structure.An Instron8802 performance testing machine (Figure3) is used for a type of MRD data test. The main structure parameters of MRD are as follows: the maximum input current of the damper is 3 A, the maximum displacement of the piston is 0.055 m, the length of the damper cylinder is 0.212 m, the outer cylinder diameter is 0.042 m, piston rod’s diameter is 0.01 m, and the weight is 2.3 kg. The sinusoidal signal is used as the excitation to do the test. By changing the amplitude of the sinusoidal signal Amp (0.005 m, 0.01 m, 0.015 m, and 0.02 m), the frequency f (0.5 Hz, 1 Hz, 1.5 Hz, and 2 Hz) and the control current I (0 A, 0.5 A, 1 A, 1.5 A, 2 A, 2.5 A, and 3 A), the MRD damping force data under different working conditions were obtained. There are too many working conditions, therefore, the working condition Amp = 0.01 m, f = 1 Hz in different currents is selected, and the experimental data are collected as shown in Figure 4. It can be seen from Figure 3 that MRD shows obvious yield characteristic and hysteresis characteristic during its working process.Figure 3
Instron8802 Testing Machine.Figure 4
Experimental date ofAmp = 0.01 m, f = 1 Hz under different currents. (a) Piston speed at different currents-damping force curve. (b) Piston displacement at different currents-damping force curve.
(a)
(b)


## 3. The Establishment of the Magic Formula Model

By studying the existing MRD parametric model (such as Bouc–Wen Model), the basic principle of parametric modelling is shown in Figure5. The MRD moves at a speed of x˙ under the action of external force F at the moment the damping force generated by the damper mainly consisting of three parts: hysteresis element, damping element, and stiffness element. The hysteresis element mainly describes nonlinear characteristics such as yield characteristics and hysteresis characteristics; the damping element describes the viscous damping after yielding; the stiffness element describes the stiffness characteristics of the damper. Part of the damper also needs to add a offset element to represent the fixed compensation force of the damper due to structural problem.Figure 5
MRD parametric model principle diagram.In order to describe the hysteresis characteristic, an absolute value term is introduced, and a complete MRD parametric model is established based on the principle of MRD parametric model modelling to generate the damping forceFMFM as is shown in the following equation:(3)FMFM=DsinCarctanB1−Ex˙+Ax+EarctanBx˙+Ax+kcx˙+kx+f0.In Equation (3), x indicates the displacement of damper piston, x˙ indicates the speed of damper piston. The other parameters’ definition and the range of their values are in Table 1.Table 1
Parameters’ definition and the range of their values.
A
The hysteresis factor, determines the half width of MRD hysteresis, generally positive
B
The stiffness factor, generally positive
C
The shape factor, determines the “S” curve shape and change around 1
D
The peak factor, determines the MRD damping force saturation, generally positive
E
The curvature factor, determines the yield characteristics, the general value is less than 1
k
c
The viscous damping coefficient
k
MRD stiffness coefficient
f
0
The bias force


## 4. Parameter Identification

For the convenience of calculation,Amp = 0.01 m, f = 1 Hz is selected as the reference condition for fitting. Curve fitting was performed by a combination of genetic algorithm and nonlinear least square method.The fitting steps are as follows. Firstly, the range of each parameter obtained from the model establishment process is optimized by genetic algorithm, and approximated global optimum can be obtained quickly by genetic algorithm in this interval, and then the approximate global optimal is taken as the nonlinear least square’s initial value of the input and the accurate global optimum of each parameter can be obtained. The data of fitting are as listed in Table2.Table 2
The fitting data ofAmp = 0.01 m, f = 1 Hz in different current conditions.
I (A)
A
B
C
D
E
k
c
k
f
0
0
0.0046
283.68
1.07
189.45
−1.34
1993.45
1272.37
−2.36
0.5
0.0049
251.30
1.05
349.28
−1.63
2534.51
1416.91
−2.37
1
0.0051
230.31
1.04
532.18
−1.89
3211.59
1600.30
−2.36
1.5
0.0052
216.04
1.02
717.44
−2.15
4065.40
1796.11
−2.37
2
0.0052
207.28
1.01
885.38
−2.38
5152.79
1981.94
−2.35
2.5
0.0052
202.71
1.00
1020.85
−2.57
6544.65
2139.31
−2.35
3
0.0051
201.68
0.99
1113.97
−2.68
8329.57
2259.40
−2.36According to the data in Table1, to identify each parameter:(1)
f0 represents the bias, as a fixed value, taking the average of the current, we can obtain(4)f0=−2.36.(2)
C, D, E, kc, and k in this condition can all be fitted as linear function of the current I, as shown in Figure 6.Figure 6
ParameterC, D, E, kc, and k and current fit graph. The points represent the experimental values, and the line represents the fitted value. (a) Parameter C and current fit graph. (b) Parameter D and current fit graph. (c) Parameter E and current fit graph. (d) Parameter kc and current fit graph. (e) Parameter k and current fit graph.
(a)
(b)
(c)
(d)
(e)The fitting results are shown in the following equations:(5)C=−0.023I+1.061,(6)D=318.463I+208.694,(7)E=−0.445I−1.469,(8)kc=296.411I+1241.836,(9)k=143.369I+1180.632.When the magnitude of the external excitation frequency changes, the values ofC, D, E, kc, k show only minor fluctuations. It can be considered that these parameters do not change with the amplitude of the external excitation frequency.(3)
A, B in this condition can be fit into a quadratic function of the current I, shown in Figure 7.Figure 7
ParameterA, B and current fit graph. The points represent the experimental values, and the line represents the fitted value. (a) Parameter A and current fit graph. (b) Parameter B and current fit graph.
(a)
(b)By changing the amplitude or frequency, the relationship ofA and amplitude is indicated in Figure 8(a), the relationship of A and frequency is indicated in Figure 8(b).Figure 8
The fitting map of ParameterA and amplitude and frequency. The points represent the experimental values, and the line represents the fitted value. (a) Parameter A and amplitude fit graph. (b) Parameter A and frequency fit graph.
(a)
(b)By changing the amplitude or frequency, the relationship of B and amplitude is indicated in Figure9(a), and the relationship of B and frequency is indicated in Figure 9(b).Figure 9
The fitting map of ParameterB and amplitude and frequency. The points represent the experimental values, and the line represents the fitted value. (a) Parameter B and amplitude fit graph. (b) Parameter B and frequency fit graph.
(a)
(b)Therefore, the final fittedA, B is as follows:(10)A=−0.205I2+0.914I+6.57849.042Amp+0.4571.042f−0.058×10−3,(11)B=8.450I2−43.804I+201.1402976.1Amp2−120.920Amp+20.804f2−2.978f+3.310.Take Equations (4)–(11) of the recognized parameters in Equation (3) to get a complete model of the MRD magic formula. Numerical simulation is taken when the complete MFM in MATLAB is under different conditions for different current, and it is compared with the experimental data. Because there are so many different working conditions, only a few of them are selected, as shown in Figure 10.Figure 10
Comparison of magic formula model and experimental data under different operating conditions. (a) WhenAmp = 0.015 m, f = 2 Hz, the piston speed-damping force curve at different currents. (b) When Amp = 0.02 m, f = 1 Hz, the piston speed-damping force curve at different currents.
(a)
(b)It can be seen from Figure9 that the parameters of MFM under different conditions are in good agreement with the experimental data.


## 5. Error Analysis and Comparison

### 5.1. PM Parameter Identification

To further investigate the accuracy and applicability of MFM, MFM is compared with PM. Taking into account the differential equation in PM is more difficult to directly model, therefore, it is chosen to establish the PM in MATLAB/Simulink, as shown in Figure11.Figure 11
Phenomenon model in MATLAB/Simulink.From previous studies, the parametersc0, c1, α, and A can be expressed as a polynomial function of the current I in PM, and other parameters are constants. The genetic algorithm is used to further identify the parameters in PM based on experimental data. Taking into account the unknown value range of each parameter of PM, the method of continuously narrowing the parameter value interval gradually improves the recognition accuracy of genetic algorithm. The specific method is as follows: set the value range of each parameter as [−1×108, 1×108], and identify each parameter by genetic algorithm, respectively, under basic condition (Amp = 0.01 m, f = 1 Hz) and different currents (0∼3 A). The minimum value of each parameter obtained by different current identification is used as the lower limit, and the maximum value is taken as the upper limit, which is reset as the parameter value range and identified by genetic algorithm. After narrowing the parameter interval several times, the final parameter identification result can be obtained (Table 3).Table 3
Phenomenon model parameters.
c
0
1.870
+
2.240
I
×
10
3
c
1
1.727
+
2.815
I
×
10
5
k
0
42.83
k
1
561.2
x
0
−
6.678
×
10
−
2
α
71.44
+
94.57
I
β
6.079
×
10
3
γ
1.083
×
10
4
A
4.731
+
0.29
I
×
10
4
n
0.9929Taking into account many more working conditions, the basic conditionAmp = 0.01 m, f = 1 Hz is selected to compare their own fitting accuracy; select the condition Amp = 0.02 m, f = 2 Hz, which is far away from the basic condition, to do the damping force prediction accuracy comparison. At the two fixed conditions, compare the error of FMFM with the error of EPM and the errors of EMFM and EPM are defined as follows:(12)EMFM=∑i=1NFMFM−FP2∑i=1NFP2,EPM=∑i=1NFPM−FP2∑i=1NFP2.


### 5.2. Self-Fitting Accuracy Comparison

MFM and PM are identified under the basic conditionAmp = 0.01 m, f = 1 Hz, and the average error for different currents under this condition is listed in Table 4. When I = 3 A, the difference in accuracy is the greatest, as shown in Figure 12.Table 4
The average current error whenAmp = 0.01 m, f = 1 Hz.
I (A)
0
0.5
1
1.5
2
2.5
3
E
MFM (%)
1.44
1.90
2.09
2.29
2.43
2.54
2.47
E
PM (%)
4.39
4.61
3.98
4.18
3.585
2.96
5.77Figure 12
The comparison between MFM, PM, and experimental data whenAmp = 0.01 m, f = 1 Hz, and I = 3 A. (a) Piston speed-damping force curve. (b) Piston displacement-damping force curve.
(a)
(b)It can be seen from the error data and Figure12 that under the condition Amp = 0.01 m, f = 1 Hz, the average error of each equation of MFM is less than 2.5% under different currents, and the average error of MFM under each current is less than the average error of PM at each current, that is, the fitting equation of MFM to the experimental data is higher.


### 5.3. Prediction Accuracy Comparison

Select the conditionAmp = 0.02 m, f = 2 Hz, which is the farthest from the basic condition, and compare the prediction accuracy of damping force. Calculate the average error under this condition as listed in Table 5. When I = 1 A, the difference in accuracy is the greatest, as shown in Figure 13.Table 5
The average current error whenAmp = 0.02 m, f = 2 Hz.
I (A)
0
0.5
1
1.5
2
2.5
3
E
MFM (%)
2.04
2.49
2.46
2.60
2.62
2.54
2.43
E
PM (%)
3.65
9.51
11.31
10.13
8.12
5.87
3.72Figure 13
The comparison between MFM, PM, and experimental data whenAmp = 0.02 m, f = 2 Hz, and I = 1 A (a) Piston speed-damping force curve. (b) Piston displacement-damping force curve.
(a)
(b)As can be seen from the error data and Figure13, under the condition Amp = 0.02 m, f = 2 Hz, the average error of the equation formula of MFM at each current is increased compared with the condition Amp = 0.01 m, f = 1 Hz, and the largest error is only 2.6%, which is still able to meet the accuracy requirements. At this time, the average error of MFM at each current is still less than the average error of PM at each current, indicating that the prediction accuracy of MFM is higher than that of PM.Based on the above error analysis, it can be seen that compared with PM, MFM has higher accuracy; when the amplitude of the external excitation frequency changes, the magnitude and standardization tendency of the damping force can still be predicted accurately, and the prediction accuracy is higher than the phenomenon model with better applicability.


## 5.1. PM Parameter Identification

To further investigate the accuracy and applicability of MFM, MFM is compared with PM. Taking into account the differential equation in PM is more difficult to directly model, therefore, it is chosen to establish the PM in MATLAB/Simulink, as shown in Figure11.Figure 11
Phenomenon model in MATLAB/Simulink.From previous studies, the parametersc0, c1, α, and A can be expressed as a polynomial function of the current I in PM, and other parameters are constants. The genetic algorithm is used to further identify the parameters in PM based on experimental data. Taking into account the unknown value range of each parameter of PM, the method of continuously narrowing the parameter value interval gradually improves the recognition accuracy of genetic algorithm. The specific method is as follows: set the value range of each parameter as [−1×108, 1×108], and identify each parameter by genetic algorithm, respectively, under basic condition (Amp = 0.01 m, f = 1 Hz) and different currents (0∼3 A). The minimum value of each parameter obtained by different current identification is used as the lower limit, and the maximum value is taken as the upper limit, which is reset as the parameter value range and identified by genetic algorithm. After narrowing the parameter interval several times, the final parameter identification result can be obtained (Table 3).Table 3
Phenomenon model parameters.
c
0
1.870
+
2.240
I
×
10
3
c
1
1.727
+
2.815
I
×
10
5
k
0
42.83
k
1
561.2
x
0
−
6.678
×
10
−
2
α
71.44
+
94.57
I
β
6.079
×
10
3
γ
1.083
×
10
4
A
4.731
+
0.29
I
×
10
4
n
0.9929Taking into account many more working conditions, the basic conditionAmp = 0.01 m, f = 1 Hz is selected to compare their own fitting accuracy; select the condition Amp = 0.02 m, f = 2 Hz, which is far away from the basic condition, to do the damping force prediction accuracy comparison. At the two fixed conditions, compare the error of FMFM with the error of EPM and the errors of EMFM and EPM are defined as follows:(12)EMFM=∑i=1NFMFM−FP2∑i=1NFP2,EPM=∑i=1NFPM−FP2∑i=1NFP2.


## 5.2. Self-Fitting Accuracy Comparison

MFM and PM are identified under the basic conditionAmp = 0.01 m, f = 1 Hz, and the average error for different currents under this condition is listed in Table 4. When I = 3 A, the difference in accuracy is the greatest, as shown in Figure 12.Table 4
The average current error whenAmp = 0.01 m, f = 1 Hz.
I (A)
0
0.5
1
1.5
2
2.5
3
E
MFM (%)
1.44
1.90
2.09
2.29
2.43
2.54
2.47
E
PM (%)
4.39
4.61
3.98
4.18
3.585
2.96
5.77Figure 12
The comparison between MFM, PM, and experimental data whenAmp = 0.01 m, f = 1 Hz, and I = 3 A. (a) Piston speed-damping force curve. (b) Piston displacement-damping force curve.
(a)
(b)It can be seen from the error data and Figure12 that under the condition Amp = 0.01 m, f = 1 Hz, the average error of each equation of MFM is less than 2.5% under different currents, and the average error of MFM under each current is less than the average error of PM at each current, that is, the fitting equation of MFM to the experimental data is higher.


## 5.3. Prediction Accuracy Comparison

Select the conditionAmp = 0.02 m, f = 2 Hz, which is the farthest from the basic condition, and compare the prediction accuracy of damping force. Calculate the average error under this condition as listed in Table 5. When I = 1 A, the difference in accuracy is the greatest, as shown in Figure 13.Table 5
The average current error whenAmp = 0.02 m, f = 2 Hz.
I (A)
0
0.5
1
1.5
2
2.5
3
E
MFM (%)
2.04
2.49
2.46
2.60
2.62
2.54
2.43
E
PM (%)
3.65
9.51
11.31
10.13
8.12
5.87
3.72Figure 13
The comparison between MFM, PM, and experimental data whenAmp = 0.02 m, f = 2 Hz, and I = 1 A (a) Piston speed-damping force curve. (b) Piston displacement-damping force curve.
(a)
(b)As can be seen from the error data and Figure13, under the condition Amp = 0.02 m, f = 2 Hz, the average error of the equation formula of MFM at each current is increased compared with the condition Amp = 0.01 m, f = 1 Hz, and the largest error is only 2.6%, which is still able to meet the accuracy requirements. At this time, the average error of MFM at each current is still less than the average error of PM at each current, indicating that the prediction accuracy of MFM is higher than that of PM.Based on the above error analysis, it can be seen that compared with PM, MFM has higher accuracy; when the amplitude of the external excitation frequency changes, the magnitude and standardization tendency of the damping force can still be predicted accurately, and the prediction accuracy is higher than the phenomenon model with better applicability.


## 6. Conclusions

In the modelling process of MRD, MFM is established after magic formula is introduced. Based on the experimental data of MRD, the parameters of the new model are identified, and the accuracy and the error analysis of the fitted MFM are compared with the commonly used PM. Under the basic conditionAmp = 0.01 m, f = 1 Hz, MFM is more accurate than PM at different currents, and the maximum precision error of MFM is only 2.5%; under the condition Amp = 0.01 m, f = 1 Hz, the fitting accuracy of MFM at different currents is also higher than that of PM, and the maximum accuracy error of MFM is only 2.6%. The results show that MFM is simple in form and does not include differential equations and intermediate variables. The physical meaning and range of values of each parameter are clear, and parameter identification and simulation calculations can be performed quickly and easily. The simulation results under different conditions are in good agreement with the experimental data. With high precision and good adaptability, it has great potential in the field of MRD control.



---

*Source: 1871846-2018-10-10.xml*