# Modelling Magnetorheological Dampers in Preyield and Postyield Regions

Authors: E. Palomares; A. L. Morales; A. J. Nieto; J. M. Chicharro; P. Pintado
Journal: Shock and Vibration
 (2019)

---

## Abstract

The use of magnetorheological dampers has rapidly spread to many engineering applications, especially those related to transportation and civil engineering. The problem arises upon modelling their highly nonlinear behaviour: in spite of the huge number of apparently accurate models in the literature, most fail when considering the overall magnetomechanical behaviour. In this study, a brief but broad review of different magnetorheological damper models has been carried out, which includes characterisation, modelling, and comparison. Unlike many other studies, the analyses cover the behaviour from preyield to postyield regions of the MR fluid. The performance of the different models has been assessed by means of numerous experimental tests and by means of simulations in a simple and straightforward semiactive control case study. The results obtained prove that most models usually fail in predicting accurate low-velocity behaviour (before iron chains yield) and, as a result, may lead to bad estimations when used in control schemes due to modelling errors and chattering.


---

## Body

## 1. Introduction

Traditional materials and fluids are being replaced by smart materials and fluids which enhance features for cutting edge technological challenges. May be the most common example of smart fluid is magnetorheological (MR) fluids, which were discovered and developed by Rabinow at the U.S. National Bureau of Standards in the late 1940s [1]. These fluids modify their rheological behaviour as a response to applied magnetic fields [2–4]. They are noncolloidal suspensions where the continuous phase is mineral oil and the dispersed phase consists of high purity iron particles on the order of a few microns in volume fractions from 20% up to 50%. The magnetically-induced dipoles form chain-like structures parallel to the applied field, which restricts the motion of the fluid and increases viscosity. The higher the magnetic field, the higher the energy needed to break those chains and hence reduce viscosity. Several studies have been carried out not only to characterise pre- and postyield regions [5–7] but also many other aspects such as particle shape, which affects durability and abrasion [8], or the type of mineral oil, which influences behaviour as proven by Lim et al. [9] and Park et al. [10].Carlson and Jolly showed the huge variety of commercial applications that these fluids open up [11]. They are used in many fields and devices but one should emphasize their use as active dampers or shock absorbers [12] in applications which may range from civil engineering (buildings and bridges [13–16]) to transportation (vehicle seats [17, 18], automobiles [19–24], railway vehicles [25–27], and even landing gears [28]).The authors have found in the literature a huge number of studies dealing with modelling of MR dampers as a function of excitation and current supplied. Topical reviews, such as that carried out by Wang and Liao [29], provide support for this statement. The simplest and most widely known model is the Bingham model [30] where the MR fluid is treated as a Newtonian fluid only from a specific threshold force or yield stress [31]. Later, Wang and Gordaninejad [32] used a Herschel–Bulkley fluid model which considers a pseudoplastic behaviour in the postfluence region [33]. More intricate models also take into account the preyield region, such as the biviscous model [34, 35] or those developed by Soltane et al. [36] which model a soft transition between pre- and postyield regions based on a Papanastasiou fluid model [37].None of these models consider the inherent and characteristic hysteresis of the force vs. velocity cycles in MR dampers despite having been thoroughly studied. Among many other authors, one can mention Soltane et al., who used a nonhysteretic model as a base [36]; Guo et al., who modified the Bingham model and considered displacement along with velocity [38]; Pang et al., with a version with branches of the biviscous model [39]; Choi et al., with a brute-force polynomial fitting [40, 41]; or Spencer et al. who linked a phenomenological model with a general solution for a large class of hysteretic systems [42].The great number of scientific papers related to modelling MR dampers clearly indicates that this is a topic far from being solved. Although almost all the proposed models show accurate behaviour, none of them seem to hit the nail on the head. In general, characterisation, modelling, and validation tests are carried out in a relatively short range of medium velocities, so neither the lower nor the higher velocities can be successfully simulated. In addition, some significant aspects such as frictions, accumulator stiffness, and even hysteresis are usually neglected, which leads to over or underestimate the goodness of the proposed active or semiactive control technique. Nevertheless, one should not overcomplicate the models since that could hinder applicability in certain fields. A good example is the use of MR models in closed-loop control systems [43], which often demand concise and robust direct models (free from instabilities, chattering, and high computational costs), as well as effective and easily derivable inverse models.The authors present in this work a brief but broad review of different MR damper models as well as a thorough experimental comparison of all the models considered. Indeed, these comparisons are fair because not only are they quantitative but they are also based on the same commercial MR damper. The experiment design also permits a better understanding of the MR damper behaviour from quasistatic to wide-range dynamic conditions and serves to figure out advantages and drawbacks of the different models. A case study of semiactive control is also shown to obtain a fair and useful comparison between models.


## 2. Characterisation of a MR Damper

Experimental characterisation and the following model comparisons have been made on a commercial magnetorheological damper RD-8040-1 manufactured by the Lord Corporation. It is a monotube shock containing high-pressure nitrogen gas (300 psi). Experimental tests have been carried out in a Material Testing System MTS-810 where the damper, in the midstroke position, was jointed in order to avoid misalignment.MR fluid can accumulate sediment if it remains in the same static position for long periods. Thus, the specimen must be excited in order to obtain a homogeneous suspension of the iron particles prior to the tests. Moreover, temperature of the fluid may affect the obtained results. Figure1(a) shows the variation of temperature of the damper body while working (sinusoidal excitation of amplitude 8 mm at 0.5 Hz) at different currents applied (0, 0.5, and 1 A). The higher the current, the higher the steady temperature and the time required to achieve it. The significance of the working temperature is shown in Figure 1(b) where the damper force is shown to be temperature dependent. Differences are noticeable (up to 40%) and the behaviour of the MR damper becomes more symmetric and predictable once the steady working temperature has been achieved.Figure 1
(a) Time response of the MR damper body temperature for three different currents (0, 0.5, and 1 A) and (b) sample of a damper force cycle with a constant current of 1 A after different elapsed times (1, 25, and 90 minutes). In both cases, the excitation is a sinusoidal excitation of 8 mm amplitude and 2 Hz frequency.
(a)
(b)Quasistatic tests were carried out in order to characterise the stiffness due to the compressed nitrogen gas chamber along with friction forces. Figure2 shows the force vs. displacement curve after contracting and extending the shaft for 3 hours (aprox. 10 μm/s). One may observe that there is a slope which defines the stiffness of the nitrogen gas chamber as well as hysteresis due to friction forces. The friction is noticeable higher when reaching the top and bottom ends of the stroke, which suggests a more complicated source of friction than the usual Coulomb friction. Both stiffness and friction can explain the nonsymmetric force-velocity curves which appear when low velocities are tested.Figure 2
Quasistatic force-displacement cycle of the MR damper (3 hours per cycle).Dynamic tests are obtained by subjecting the damper to a sinusoidal excitation, where the amplitude is 4, 8, 12, or 16 mm, whereas the frequency is set to 0.1, 0.5, 1, 2, or 4 Hz. The resulting 20 combinations cover a wide range of velocities from 2.5 to 400 mm/s, which is considerably higher than those studied by other authors. These tests allows for studying different phenomenological situations:(a)
Low (L) velocity tests where preyield region dominates (up to 15 mm/s)(b)
Medium (M) velocity tests where pre- and postyield regions are balanced (up to 100 mm/s)(c)
High (H) velocity tests where postyield region dominates (over 100 mm/s)In addition, five different intensities have been considered (0, 0.25, 0.50, 0.75, and 1 A) which leads to a total of 100 different characterisation tests.The tests conducted permit the usual force-velocity hysteresis features (well described in [29]) to be characterised: it progresses along a counter-clockwise path, is strong only in the preyield domain, exhibits roll-off in the vicinity of zero velocity, and seems to depend on several parameters. Leaving aside the obvious intensity dependence, Wang and Liao argued that hysteresis shows displacement amplitude and frequency dependence for sinusoidal experimental tests where z=dsin2πft [29], but this means that hysteresis depends on maximum velocity: 2πdf. In order to determine dependences different from velocity, Figure 3 compares tests in which maximum velocity remains constant. Whether differences should be attributed to displacement (d), acceleration (4π2df2), or both of them is not clear, but not taking this behaviour into account would deteriorate model performance.Figure 3
Force-velocity hysteresis cycles for different experimental sinusoidal tests. (a)I=0A, (b) I=1A.
(a)
(b)


## 3. Review of MR Damper Models

In order to carry out a straightforward description of the state of the art related to modelling of MR dampers, all the parametric models referred to are expressed with the same (if possible) nomenclature. Figure4 shows the most recurrent variables in the parametric models. F is the actual damper force, but the models will describe the net force F−f0, that is, the force after subtracting the contribution of the nitrogen gas accumulator when the shaft remains in the midstroke position (f0). z refers to the shaft displacement from the midstroke position, where a positive value means a contraction of the rod. C0 and C∞ are usual damping constants which would linearly define the force-velocity curves in the preyield and postyield regions. In fact, these regions have a frontier in the velocity z˙y, whereas fy refers to the force that makes the iron chains yield (fy) in the Bingham model (that is, the origin ordinate for the C∞ line). For those models which predict hysteresis, z˙h means the zero force velocity intercept. Adopting such a nomenclature for all the models will shed some light on the currently fuzzy state of the art and will demonstrate that proposed models are not as different as they may initially seem.Figure 4
Main MR model parameters.In 1922, Eugene C. Bingham proposed the mathematical form of a viscoplastic model for those materials that behave as a rigid body at low stresses and as a viscous fluid at high stresses [31]. Perhaps from the works of Stanway et al. in the late 1980s [30], the so-called Bingham model has been used to model the rheological behaviour of MR fluids easily. The idealised mechanical model consists of a Coulomb friction element (related to the fluid yield stress) in parallel with a viscous damper:(1)F−f0=fysignz˙+C∞z˙.Wang and Gordaninejad in [32] argued that the postyield viscoplastic behaviour in MR fluids is not constant as assumed by the Bingham model but agrees with a non-Newtonian Herschel–Bulkley fluid model [33]. Although the mathematical form of this model was developed to describe MR fluid dynamics through pipes and parallel plates, this behaviour can be used to predict the MR damper force:(2)F−f0=fysignz˙+C∞z˙nsignz˙.This formulation provides a nonlinear behaviour in the postyield region which more closely resembles experimental data.Rather than assuming that the MR fluid is rigid in the preyield condition, as the Bingham model does, Stanway et al. [35] adopted the nonlinear biviscous model where the MR fluid is plastic in both the pre- and postyield conditions:(3)F−f0=C0z˙,ifz˙<z˙y,fysignz˙+C∞z˙−signz˙z˙y,ifz˙≥z˙y,where continuity in the piecewise function is given by the fact that fy=C0z˙y. In fact, the Bingham model can be regarded as a limiting case of this model in which z˙y⟶0 and hence C0=fy/z˙y⟶∞.Soltane et al. returned to the starting point very recently and proposed a new modification of the Bingham model which also includes a preyield region [36]. Unlike the biviscous model, the regularised Bingham model (RB in the following) proposes a continuous equation based on the Papanastasiou fluid model that smoothes the yield criterion [37], that is,(4)F−f0=fysignz˙1−e−signz˙z˙/z˙y+C∞z˙,where z˙y, although it can be still be considered as the transition velocity, is more properly an artificial regularisation parameter with the dimension of velocity which controls the exponential growth of the damping force. The smaller the parameter, the more similar the approximation to the initial Bingham model.Until now, the models described have shown an evolution which tends to model a biviscous behaviour with a soft transition between pre- and postyield regions, leaving aside the hysteretic behaviour of the MR damper. One of the first approaches to the hysteresis modelling was carried out by Gamota and Filisko with an extension of the Bingham model [44]. It is a viscoelastic-plastic model (viscoelastic in preyield region and plastic in postyield region) based on the Bingham model in series with a standard model of a linear solid [45]. Nevertheless, this model is usually discarded because its governing equations are extremely stiff.A simpler solution to include the hysteretic phenomenon is to duplicate any nonhysteretic model into two different branches: one for contraction and another for extension. May be the easiest approach (which includes pre- and postdomains) is the hysteretic biviscous model developed by Pang et al. in 1998 [39]. It is an extension of Stanway et al.’s biviscous model with an improved representation of the preyield hysteresis. This is accomplished by adding a fourth parameter, that is, the zero force velocity intercept z˙h which causes a translation in the axis of abscissas (z˙) as indicated by the following piecewise continuous function:(5)F−f0=C∞z˙−fy,ifz˙<−z˙ycandz¨>0,C0z˙−z˙h,if−z˙yc≤z˙≤z˙yeandz¨>0,C∞z˙+fy,ifz˙ye<z˙andz¨>0,C0z˙+fy,ifz˙yc<z˙andz¨<0,C∞z˙+z˙h,if−z˙ye≤z˙≤z˙ycandz¨<0,C0z˙+fy,ifz˙<−z˙yeandz¨<0,where the contraction yield velocity z˙yc and the extension yield velocity z˙ye are also introduced and given by the expressions:(6)z˙yc=fy−C0z˙hC0−C∞,z˙ye=fy+C0z˙hC0−C∞.In the same vein, Soltane et al. proposed the hysteretic regularised Bingham model (HRB) which modifies the nonhysteretic version by translating pre- and postyield regions to the left or right (a quantity defined by the zero force velocity interceptz˙h) depending on the compression or extension [36]:(7)F−f0=C∞z˙−z˙h+fy1−e−signz˙−z˙hz˙−z˙h/z˙ysignz˙−z˙h,ifz¨<0,C∞z˙+z˙h+fy1−e−signz˙+z˙hz˙+z˙h/z˙ysignz˙+z˙h,ifz¨>0.Guo et al. took a leap forward in 2006 when they developed a more concise and non-piecewise model which describes not only a smooth biviscous behaviour but also the hysteretic behaviour of MR dampers [38]:(8)F−f0=fytanhC0z˙+Ωz+C∞z˙+Ωz,Ω is defined as the absolute value of the hysteretic critical velocity divided by the absolute value of the hysteretic critical displacement, that is, zh/z˙h. The authors of this model emphasize its ability to predict the performance of MR dampers accurately with parameters which have a definite physical meaning while avoiding the inconveniences of other more complicated models.Some authors have opted to use nonparametric modelling methods based upon analytical expressions which avoid any physical interpretation. This is the case with the validated multifunction model [46], black-box model [47], query-based model [48], neural network model [49], fuzzy model [50], Ridgenet model [51], or the common polynomial model [40, 41]. This latter consists of a brute-force polynomial fitting process, the piecewise expression of which is as simple as(9)F−f0=∑i=0naiz˙i,ifz¨<0,∑i=0nbiz˙i,ifz¨>0.This model leads to2n+1 coefficients, which is a huge number considering that it has been proven that at least 6th order polynomials are required, and at least 12th order polynomials are recommended [41]. In addition, it is well known that these polynomials are extremely sensitive to velocities which exceed the fitting domain.Returning to parametric models, hysteresis operator-based dynamic models are also very common in the literature. They are conceived to represent a large class of hysteretic behaviour by means of a mathematical formulation which usually includes differential equations. One can find the Dahl hysteresis operator [52, 53], LuGre hysteresis operator [54, 55], the hyperbolic tangent function operator [56], and the well-known Bouc–Wen hysteresis operator (initially formulated by Bouc [57] and later generalised by Wen [58]). A modification of the Bouc–Wen model (the so-called modified Bouc–Wen model or Spencer model [42]) is widely used in the literature obtaining high accuracy in its force predictions. This phenomenological model is governed by the following equations:(10)F−f0=C0z˙mr+K0z,z˙mr=1C0+C∞αzbw+C∞z˙+K∞z−zmr,z˙bw=−γz˙−z˙mrzbwn−1zbw−βz˙−z˙mrzbwn+Az˙−z˙mr,where zmr is an intermediate variable for modelling pre- and postyield behaviour, K0 is the accumulator stiffness, K∞ controls stiffness at large velocities, and α, β, γ, A, and n, as well as the evolutionary variable zbw, define a general Bouc–Wen hysteresis [58]. Despite its proven accuracy, the extended number of model parameters (10 without considering the influence of the current applied) leads to difficulties in parameter identification.Intensity dependence has not yet been mentioned, but all the authors deal with this problem in the same way: they fit the proposed model for different intensities and the resulting parameters are subsequently fitted through polynomials with the power of applied current. The order of these polynomials varies with the author. Some of them use a simple first-order polynomial for all the current-dependent variables [36, 41, 42], whereas other authors carry out a more detailed study leading from second- to fourth-order polynomials depending on the variable considered [39]. In general, experimental force-velocity data for MR dampers suggest a nonlinear dependence of the force on the current applied.In any case, magnetic solenoid dynamics (usually modelled simply as an inductor in series with a resistor) should be included in order to obtain realistic results. This can be simply achieved by including a first-order differential equation in terms of voltage (or current) as Spencer et al. proposed [42]:(11)u˙=−ηu−v,where v is the voltage applied to the current driver, u the instantaneous voltage in the magnetic solenoid, and η a time constant.


## 4. Experimental Comparison of MR Damper Models

In order to carry out a fair comparison among MR models, all of them were fitted by using the same experimental data from a RD-8040-1 damper. The fitting procedure of the different models was conducted with a nonlinear least-squares method in Matlab software. Quantitative comparisons will be carried out by measuring the normalised root mean square error (NRMSE), that is, the usual root mean square error normalised by the range (maximum value minus minimum value).A first approach to the models confirms that they are all capable of predicting to a greater or lesser extent the experimental behaviour of a MR damper. That fact can be proven by fitting individually several experimental tests. Figure5 shows how the nonhysteretic models behave for four cases: a low velocity test at I=0A (Figure 5(a)), a low velocity test at I=1A (Figure 5(b)), a high velocity test at I=0A (Figure 5(c)), and a high velocity test at I=1A (Figure 5(d)). The same tests are later replicated with the hysteretic models in Figure 6.Figure 5
Nonhysteretic model predictions for four different experimental tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) high velocity without intensity; (d) high velocity with I=1A.
(a)
(b)
(c)
(d)Figure 6
Hysteretic model predictions for four different experimental tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) high velocity without intensity; (d) high velocity with I=1A.
(a)
(b)
(c)
(d)Differences in the parameters when different currents are applied (0 or 1 in these cases) are easily modelled by fitting a polynomial. The main problem arises when the parameters of the models, for a constant current, are different depending on the range of velocities that needs to be covered. That is what happens in Figures5 and 6: low velocity tests (where preyield dominates) and high velocity tests (where postyield dominates), for the same current, require very different values of the model parameters to provide an accurate performance.The influence of the range of velocities in the fitting procedure has been analysed by taking into account two different sets of experimental data:(a)
All the experimental results, which comprise low, medium, and high velocity tests (in the following LMH fitting)(b)
Only the habitual medium and high velocity tests, that is, discarding those where the postyield region is not fully developed and the preyield region dominates (in the following MH fitting).Predictions obtained from LMH fitting and LM fitting are later compared to six experimental results which cover three different maximum velocities (low, medium, and high) and two applied intensities (I=0A and I=1A). The sinusoidal excitation is set to A=4mm and f=0.5Hz for the low velocity test, A=4mm and f=4Hz for the medium velocity test, and A=16mm and f=2Hz for the high velocity test.The Bingham model (Figure7) stands out because of its simplicity and plainness, but returns high errors in general because it neglects both hysteresis and the preyield region. MH fitting obtains errors up to 24% in low velocity validations and up to 16% otherwise, whereas LMH fitting slightly improves those results (up to 19% and 16%, respectively).Figure 7
Bingham model (fitted with LMH and MH tests) predictions for six different tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) medium velocity without intensity; (d) medium velocity with I=1A; (e) high velocity without intensity; (f) high velocity with I=1A.
(a)
(b)
(c)
(d)
(e)
(f)Herschel–Bulkley model (Figure8) behaves reasonably well in all the velocity ranges since it smoothes the yield transition, but the lack of hysteresis also leads to medium errors. For MH fitting, NRMSEs are between 6% and 12%, whereas LMH fitting slightly improves predictions if yield transition has not totally developed yet but also slightly worsens predictions otherwise.Figure 8
Herschel–Bulkley model (fitted with LMH and MH tests) predictions for six different tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) medium velocity without intensity; (d) medium velocity with I=1A; (e) high velocity without intensity; (f) high velocity with I=1A.
(a)
(b)
(c)
(d)
(e)
(f)The biviscous model (Figure9) and regularised Bingham model (Figure 10) behave very similarly to the Herschel–Bulkley model. For MH fitting, they provide slightly lower errors except for low velocity predictions where these predictions clearly fail and NRMSEs rise to 23% and 28%, respectively. For LMH fitting, these two models fix their problems in the low velocity range (NRMSE around 10%) but the preyield region almost disappears if high velocities are reached. Thus, these models degenerate into the Bingham model when fitted this way.Figure 9
Biviscous model (fitted with LMH and MH tests) predictions for six different tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) medium velocity without intensity; (d) medium velocity with I=1A; (e) high velocity without intensity; (f) high velocity with I=1A.
(a)
(b)
(c)
(d)
(e)
(f)Figure 10
Regularised Bingham model (fitted with LMH and MH tests) predictions for six different tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) medium velocity without intensity; (d) medium velocity with I=1A; (e) high velocity without intensity; (f) high velocity with I=1A.
(a)
(b)
(c)
(d)
(e)
(f)Low NRMSEs do not appear until hysteretic models are used. On the one hand, for MH fitting, the hysteretic biviscous model (Figure11), the hysteretic regularised Bingham model (Figure 12), or the Shuqi-Guo model (Figure 13) significantly improves errors if the postyield domain is well developed (between 2.5% and 6.5%) but their behaviour is erratic if the preyield region dominates (NRMSE up to 40%). On the other hand, for LMH fitting, said erratic behaviour when preyield region dominates disappears but this is very harmful when simulating high velocities: NRMSEs increase (6–13%) because hysteresis almost vanishes and these models degenerate into the Bingham model.Figure 11
Hysteretic biviscous model (fitted with LMH and MH tests) predictions for six different tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) medium velocity without intensity; (d) medium velocity with I=1A; (e) high velocity without intensity; (f) high velocity with I=1A.
(a)
(b)
(c)
(d)
(e)
(f)Figure 12
Hysteretic regularised Bingham model (fitted with LMH and MH tests) predictions for six different tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) medium velocity without intensity; (d) medium velocity with I=1A; (e) high velocity without intensity; (f) high velocity with I=1A.
(a)
(b)
(c)
(d)
(e)
(f)Figure 13
Shuqi-Guo model (fitted with LMH and MH tests) predictions for six different tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) medium velocity without intensity; (d) medium velocity with I=1A; (e) high velocity without intensity; (f) high velocity with I=1A.
(a)
(b)
(c)
(d)
(e)
(f)The nonparametric 12th order polynomial model (Figure 14) is capable of obtaining low NRMSEs but only if the postyield domain is well developed. In addition, it presents a wavy shape and errors may be huge if velocities leave the fitting domain. NRMSEs are between 4% and 30% for MH fitting and 6% and 24% for LMH fitting. Lower polynomial order would compromise prediction accuracy, whereas higher orders would enhance its drawbacks.Figure 14
Polynomial model (fitted with LMH and MH tests) predictions for six different tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) medium velocity without intensity; (d) medium velocity with I=1A; (e) high velocity without intensity; (f) high velocity with I=1A.
(a)
(b)
(c)
(d)
(e)
(f)Finally, the Spencer model (Figure15) stands out from the others for providing the best qualitative and quantitative predictions when the postyield region dominates (NRMSEs between 1.7% and 4.4% for any fitting method). Predictions when the preyield region dominates are still accurate (between 4% and 8.1% for any fitting method), and they show no erratic behaviour in any case. Nevertheless, the fitting procedure is too complicated due to the high number of parameters: it requires a very precise seed and a trial-and-error procedure is usually worthwhile.Figure 15
Spencer model (fitted with LMH and MH tests) predictions for six different tests: (a) low velocity without intensity; (b) low velocity withI=1A; (c) medium velocity without intensity; (d) medium velocity with I=1A; (e) high velocity without intensity; (f) high velocity with I=1A.
(a)
(b)
(c)
(d)
(e)
(f)Tables1 and 2 show the fitted parameters which have a physical meaning and are shared by the different models. One can observe that not only are the values obtained similar in the different models but also close to those expected by inspection. Only one exception arises when analysing fy in the Herschel–Bulkley model: the least-squares fitting procedure leads to a low (and even zero if I=1A) value of fy because parameter n seems to be enough to model the pre- and postyield domains; nevertheless, if one sets an estimated fix value of fy, NRMSE and the qualitative behaviour of this model are equivalent.Table 1
Model parameters after fitting withI=0A.
Model
f
y (N)
C
∞ (Ns/m)
C
0 (Ns/m)
z
˙
y (mm/s)
z
˙
h (mm/s)
LMH
Bingham
75.26
961.87
—
—
—
Herschel–Bulkley
41.62
548.00
—
—
—
Biviscous
75.80
961.25
17124
4.7
—
Regularised Bingham
79.26
933.13
—
2.8
—
Hysteretic biviscous
76.44
955.97
20330
2.1
1.8
Hysteretic regularised Bingham
80.40
924.53
—
2.7
1.7
Shuqi-Guo
79.11
934.93
263
—
—
Spencer
—
961.25
20330
—
—
MH
Bingham
89.67
842.82
—
—
—
Herschel–Bulkley
19.69
549.58
—
—
—
Biviscous
86.28
869.84
5612
1.8
—
Regularised Bingham
89.80
851.50
—
9.9
—
Hysteretic biviscous
90.76
842.74
7579
3.1
9.2
Hysteretic regularised Bingham
92.09
829.93
—
7.0
7.9
Shuqi-Guo
90.60
844.22
85
—
—
Spencer
—
964.31
20156
—
—Table 2
Model parameters after fitting withI=1A.
Model
f
y (N)
C
∞ (Ns/m)
C
0 (Ns/m)
z
˙
y (mm/s)
z
˙
h (mm/s)
LMH
Bingham
963.1
3508.5
—
—
—
Herschel–Bulkley
0.0
2333.8
—
—
—
Biviscous
823.5
4503.7
147230
5.8
—
Regularised Bingham
845.5
4324.1
—
3.0
—
Hysteretic biviscous
837.4
4385.8
162000
3.3
2.0
Hysteretic regularised Bingham
863.7
4195.1
—
3.0
2.0
Shuqi-Guo
855.8
4293.8
222
—
—
Spencer
—
2401.3
162640
—
—
MH
Bingham
1046.7
2858.9
—
—
—
Herschel–Bulkley
0.0
2997.5
—
—
—
Biviscous
1005.0
3135.9
27392
4.1
—
Regularised Bingham
1137.0
2302.4
—
25.0
—
Hysteretic biviscous
1027.5
2974.6
40831
8.2
17.6
Hysteretic regularised Bingham
1099.9
2436.5
35
16.3
16.2
Shuqi-Guo
1084.7
2589.6
—
—
—
Spencer
—
2412.1
160680
—
—In short, one may come to the following main conclusions:(i)
If one uses MH fitting, the behaviour of the model when yield transition is not fully developed is clearly poor because predictions underestimate MR damper forces in the preyield region. This compromises the simulation results of control systems which, in fact, are specifically designed to achieve low vibration velocities (or accelerations).(ii)
If one uses LMH fitting, the models provide an accurate performance if the postyield domain is not fully developed but hysteresis almost vanishes and force is overestimated otherwise. This diminishes the claimed good performance of the complex and refined models, degenerating into the old Bingham model.(iii)
The Spencer model provides the best qualitative and quantitative behaviour, but the seed for the fitting procedure must be so accurate that least-squares methods are almost unnecessary.


## 5. Case Study: Skyhook Control of a Quarter-Car Suspension System

A quarter-car semiactive suspension system (shown in Figure16) with an on-off skyhook control strategy has been numerically analysed in order to determine the influence of the selected MR model. The chosen suspension parameters were M=500kg and K=78000N/m together with the characterised and modelled magnetorheological damper RD-8040-1.Figure 16
The quarter-car suspension system.Figure17 shows the excitement in the base which replicates a speed bump. The mathematical expression is(12)z0t=H2sin2πvLt−0.5−π2+1,0.5≤t≤Lv+0.5,0,otherwise,where H and L are the bump height and length, respectively, and v is the vehicle velocity (constant).Figure 17
Temporal excitation corresponding to a bump withH=0.02m and L=1m travelled at v=15km/h and example of the current required by the skyhook control system.The on-off skyhook control strategy can be defined as follows [59]:(13)Fsa=Fmax,z˙z˙−z˙0≥0,Fmin,z˙z˙−z˙0<0,where Fmax and Fmin are the MR damper forces when a current I=1A and I=0A, respectively, is supplied. An example of the current required by the skyhook control system, superimposed on the bump, is shown in Figure 17 (differences depending on the MR model are imperceptible). Note that, though not seen in Figure 17, the damper will switch off after the sprung mass has reached its equilibrium position again.Figures18 and 19 show the sprung mass acceleration (related to comfort in vehicles), both controlled and passive, for all the MR models considered after either MH fitting and LMH fitting. In all the analysed models, passive results for the two fitting cases are almost identical because the damper remains in the postyield region, where differences are limited. Nevertheless, controlled results highlight noticeable differences.Figure 18
Time response of the controlled (solid line) and passive (dotted line) sprung mass acceleration for the (a) Bingham model, (b) Herschel–Bulkley model, (c) biviscous model, and (d) regularised Bingham model, after LMH fitting (black) and MH fitting (grey).
(a)
(b)
(c)
(d)Figure 19
Time response of the controlled (solid line) and passive (dotted line) sprung mass acceleration for the (a) hysteretic biviscous model, (b) hysteretic regularised Bingham model, (c) Shuqi-Guo model, (d) polynomial model, and (e) Spencer model, after LMH fitting (black) and MH fitting (grey).
(a)
(b)
(c)
(d)
(e)Results in the control with the Bingham model (Figure18(a)) show chatter in the controlled sprung mass acceleration. This is due to the fact that this model does not consider a preyield region and there is a discontinuity between positive and negative postyield regions (see Figure 7). This occurs using both MH fitting and LMH fitting.Controls with nonhysteretic Herschel–Bulkley, biviscous, and regularised Bingham models show a similar behaviour (Figures18(b)–18(d)). Both the first peak (around 0.15 s), which corresponds to the ascent to the top of the bump, and the second (around 0.25 s), which is due to the drop from the top of the bump, are almost identical notwithstanding the tests which were used to fit the models. This happens because the relative velocity of the MR damper is still high and the postyield region (very similar in all the cases) dominates dynamics. Nevertheless, the third peak (around 0.40 s), which is due to a rebound, happens at lower velocities where yield transition is not fully developed and models have been proven to disagree. These three models, if MH fitted, predict a slower attenuation of the sprung mass acceleration because the dissipative force is underestimated (see Figures 8–10, subplots (a) and (b)). This fact is especially noticeable when considering the regularised Bingham model because it is the one which underestimates the force in a greater extent. Thus, the faster attenuation predicted by the model after LMH fitting would be more reliable.The use of hysteretic models such as the hysteretic biviscous, HRB, and 12th order polynomial models (Figures 19(a), 19(b), and 19(d), respectively) does not contribute perceptible improvements in comparison to the results obtained with the previous nonhysteretic models. Indeed, they are hysteresis models based on two different branches (depending on the sign of the relative acceleration of the damper), which may lead to low velocity chatter and discontinuities as shown in Figures 19(a), 19(b), and 19(d). In addition, since there is no branch which replicates the static state of the MR damper (that is, zero force when zero velocity), the sprung mass will achieve a different static position depending on either the supplied intensity and the active branch when it stops. As a consequence, these models provide acceleration time responses which are far from being acceptable.Results for the Shuqi-Guo model (Figure19(c)) avoid the control drawbacks inherent in the branched hysteresis models (chatter, static position, etc.,) but results are almost identical to those obtained with the nonhysteretic models: if MH fitting is used, the low velocity force is underestimated; if LMH fitting is used, hysteresis almost vanishes.Finally, results for the Spencer model (Figure19(e)) show a lower dissipation than the predicted by the rest of the models (rebounds reach slightly higher acceleration levels) and almost identical predictions in both fitting cases. Note that the Spencer model was proven to provide the most accurate behaviour in all the velocity ranges and for the two considered fittings.Tables3 and 4 show the rms value of the sprung mass acceleration. It is not a surprise that, even discarding chattering models, predicted rms values are in the range of 0.541–0.734 m/s2, which means variations up to 36% depend on the model one decides to use.Table 3
Root mean square (rms) value of the sprung mass acceleration for simulations with nonhysteretic models (units in m/s2).
MR model
Bingham
Herschel–Bulkley
Biviscous
Regularised Bingham
LMH
1.123
0.559
0.543
0.542
MH
1.178
0.582
0.627
0.596Table 4
Root mean square (rms) value of the sprung mass acceleration for simulations with hysteretic models (units in m/s2).
MR model
Hysterical biviscous
HRB
Shuqi-Guo
Polynomial
Spencer
LMH
0.674
0.531
0.541
0.944
0.734
MH
1.516
1.274
0.607
1.456
0.729We draw the following conclusions about the influence of the MR model on semiactive control simulations:(i)
Systems including discontinuous or piecewise models (such as the Bingham model, hysteretic biviscous model, HRB model, or polynomial model) are prone to chatter, discontinuities, and variations in the static deflection.(ii)
Nonhysteretic models provide similar predictions to hysteretic models because hysteresis actually affects a small fraction of an oscillation.(iii)
Differences between MH fitting and LMH fitting are mainly noticeable at low-velocity oscillations, where MH fitting underestimates the damping force and leads to a slower attenuation.(iv)
The Spencer model shows a lower dissipation (and hence worse rms accelerations) than the predicted by the other models, but it was shown to be the most accurate model.(v)
Nonnegligible influences such as friction are not included in the models, which mainly compromises low velocity predictions.


## 6. Conclusions

This article performs an updated review of the main MR damper models. Unlike other studies, this review includes a broad and thorough experimental comparison of nine models to predict the behaviour of the same Lord RD-8040-1 MR damper.The variety of experimental tests used for the fitting process has been proven to be a key decision: if tests where the preyield region dominates are discarded, low-velocity predictions are either poor or erratic; if these tests are also considered, hysteresis at higher velocities almost vanishes. The Spencer model avoids this problem and stands out for accurately predicting pre- and postyield regions with NRMSE lower than8% in all cases.These models have also been used in a straightforward semiactive control case study. Several problems have arisen: on the one hand, discontinuous or piecewise models produce chatter, discontinuities, and even variations in the static deflection; on the other hand, MH fitting underestimates damping forces and predicts slower attenuations. The a priori accurate Spencer model predicts a slower attenuation in comparison to the rest of the models.In the light of these results, there is still room for further research on MR models. NRMSE should be reduced in all the range of velocities which cover pre- and postyield domains, and proposed models must guarantee the lack of habitual problems such as chatter, discontinuities, or inaccuracies prior to claiming the goodness of any control scheme.



---

*Source: 5636053-2019-05-02.xml*