Dual Hysteresis Model of MR Dampers

This study concerns the modeling of the hysteretic behavior of magnetorheological (MR) dampers. In general, hysteresis is one of key factors influencing the output of such actuators. So far, more attention has been paid to studying the combined hysteretic behavior of MR actuators by observing the relationships between the output (force/torque) and the inputs (current, velocity, and position). However, these devices feature two distinct hysteretic mechanisms: mechanical/hydraulic and magnetic. The mechanical hysteresis is of different nature than the magnetic hysteresis due to the properties of ferromagnetic materials forming the actuator's electromagnet circuit, and these should be split in the modeling process. In the present study, we separate the magnetic hysteresis from the mechanical/hydraulic one by investigating the magnetic flux vs. exciting current relationship of a commercial flow-mode MR damper subjected to sinusoidal current loading and independently of the mechanical excitations. The resulting behavior of the electromagnetic circuit is then examined using the non-linear inductor approach with hysteresis. Total hysteresis is then modeled using a non-linear inductor model in combination with a phenomenological parametric Maxwell type model of the damper.


INTRODUCTION
Magnetorheological (MR) fluids are well-known representatives of so-called smart materials. The material develops a yield stress when exposed to magnetic field (Rabinow, 1948), and it has been successfully used in commercial applications, i.e., semi-active vehicle dampers or powertrain mounts (Jolly et al., 1999).
Fundamentally, a typical MR flow-mode/shear-mode damper features an internal control valve in the form of an electromagnet with a fixed height planar/annular flow channel. The electromagnet's core contains a coil assembly. Supplying the electrical current to the coil results in inducing magnetic field in the flow channel, thus activating the fluid. The effect is a resistanceto-flow build-up manifested by changes in the output force or torque. The effect is reversible and fast. At the same time, various factors make the conversion process complicated (Gołdasz et al., 2018a), namely, temperature, friction, material's liability to sedimentation, non-linear magnetization characteristics of the materials of which the damper's magnetic circuit is built, magnetic hysteresis, mechanical hysteresis, current driver dynamics, control coil dynamics, nonlinear relationship between the magnetic flux and the field-dependent yield stress, etc. (Kubik et al., 2017). Generally, the contributors make the control algorithm development process a challenge and the optimal control algorithm difficult or impossible to obtain. Of the above contributors, the magnetic hysteresis is present virtually in any ferromagnetic material, though on a different scale. Although the electromagnetic circuits of MR valves are usually developed with soft magnetic alloys, the flux's hysteresis should be accounted for in engineering an MR damper-based control system. As MR dampers are by principle solenoid actuators (in which the topic has been subject to extensive research, , the significance of hysteresis is similarly important. In latching two-position solenoids, for instance, the hysteresis and the residual magnetism affect the control circuit's capacity to maintain the solenoids at a given position. In MR dampers the hysteresis complicates controlling the output of these devices, and residual forces have a usually negative impact on the damper's output at off-state forces in particular. In literature, magnetic hysteresis models are classified into two groups (Mazgaj, 2010). The first one includes energy-based models of which the Stoner-Wohlfart (S-W) model or Jiles-Atherton (J-A) model are the most wellknown representatives (Jiles and Atherton, 1984). The other class incorporates the so-called phenomenological models of Preisach (Mayergoyz and Friedman, 1988;, Duhem/Hodgdon-Coleman (Macki et al., 1993) or Chua (Chua and Stromsmoe, 1970). For this study, we chose one representative of the phenomenological group-the Duhem model.
The contribution of hysteresis to the force output of MR dampers has been recognized rather early (Dyke et al., 1996). Following on from Dyke et al., numerous phenomenological or lumped parameter models have been developed and applied for the purpose of predicting the hysteretic output of MR dampers (Wang and Liao, 2011). Typical approach involved studying the relationship between force/torque vs. input current and displacement/velocity and then fitting a particular model to damper data from experiments. With such a methodology, only total hysteresis can be studied. Few authors attempted to analyze the magnetic hysteresis independently of other hysteretic contributors (Szczygłowski, 2001;An and Kwon, 2003;Zheng et al., 2015;Guo et al., 2018;Li et al., 2019). The mechanical/hydraulic hysteresis has its origin in flow mechanics due to the fluid's compressibility and the lumped fluid mass (inertia) traveling back and forth through the annulus (Gołdasz and Alexandridis, 2012). For comparison, the magnetic hysteresis is an inherent property of ferromagnetic materials and does not disappear in the DC limit, which is what the hydraulic hysteresis does. Including the magnetic hysteresis operator, therefore, brings numerous benefits. First, it facilitates control through either sensor-based or sensorless approach as proven with other electrical actuators and drives (Erol et al., 2012). The implementation of the magnetic hysteresis operator will allow for a model-based control, too. Second, it allows taking into account the effects already in the design process. For instance, Kubik and Gołdasz (2019), in the parametric study on the MR damper dynamic behavior with the so-called vector hysteresis approach, showed that the use of such models may yield improvements in understanding the MR damper physics and deliver a better product. Simply analyzing prior art, e.g., the study of hysteresis of an MR brake as in An and Kwon (2003) accomplished with the Hodgdon model highlights significant advantages in predicting the output of the brake on including the magnetic hysteretic operator.
In general, the purpose of the study is to present means for distinguishing the contribution of particular hysteretic operators to the output of MR dampers. To accomplish the goal, we reconstruct the magnetic behavior of a longstroke flow-mode MR damper with one magnetic hysteretic operator to obtain the induced flux variation against the exciting current. The sensorless flux reconstruction technique we used allows inspection of the damper's magnetic hysteretic behavior without any knowledge of its internals. Finally, a first-order phenomenological damping force model based on a Maxwell model is employed to explore the contribution of the remaining mechanisms.
The paper is organized as follows. Section 2 contains the description of the dual hysteresis model concept of an MR damper. In the section we describe both the magnetic hysteretic operator and the mechanical model. Furthermore, in section 3, we explain the test setup and the hardware used. Next, section 4 presents model identification results and a comparison of the hysteretic models' performance against selected experimental data sets. Finally, conclusions are presented in section 5.

DAMPER MODELING
To predict the hysteretic behavior of a typical single-tube flow-mode MR damper as in Figure 1, we developed a phenomenological model of the device. The cylinder tube houses a piston assembly that separates the fluid volumes into the upper (rebound) volume and the lower (compression) volume. The floating gas cup separates the compression chamber from pressurized gas for volume compensation. The piston features a solenoid with a fixed-height annular gap for magnetizing the MR fluid flowing through it as in Figure 1. There is a thru-rod electrical connection between the current driver (not shown) and the coil. The induced flux travels across the core, the annular gap, the ring and the cylinder tube to return to the core through the annular gap again. As illustrated, there is also a flux leakage path from the tube through the MR fluid and the rod. Specifically, we model viscous losses, friction, and offset force due to the accumulator. Moreover, we assume the model would capture the effects of fluid compressibility, field-dependent pressure drop, as well as basic dynamics of the MR actuator's electromagnet although in a simplified manner. The model incorporates two fundamental components reflecting the force build-up process in the actuator: the electromagnet model and the MR effect model.
First, the dynamics of the electromagnet is captured with a non-linear hysteretic inductor concept. The model can be analyzed by looking up the relationship between the exciting (magnetizing) current and the resulting induced flux (linkage). Second, the behavior of the fluid is examined further with a first-order phenomenological model of the damping force. The inductor model and the force model are related by a non-linear flux-to-force coupling.

Non-linear Inductor Modeling With Hysteresis
Let us consider the simplest non-linear lumped parameter model of an MR actuator as shown in Figure 2. The model includes the input voltage source u(t), the coil resistance R c , and the non-linear inductance L(i) in series with the resistor. The model equation is then as follows: where u -supply voltage, i -coil current, and λ = L(i)imagnetic flux linkage, and L(i) -coil inductance. Equation (1) represents the dynamics of a simple non-linear inductor. The relationship λ(i) can be further expanded to include other phenomena.
In the presented form the model copies the average flux induced in the electromagnet, and its parameters can be estimated independently of the hysteretic force component. Equation (1) may represent a simple non-linear inductor model with no hysteresis or assume more complex forms. In the present study we chose to proceed with the Duhem hysteresis model.

Duhem Model
The relationship between the flux linkage λ and the coil current i of the inductor in Figure 2 can be described using the following Duhem hysteresis operator We select the shape functions f (i) and g(i) as follows:  The model incorporates a set of four tuning parameters (a, b 1 , b 2 , d), which can be identified from experimental data. In the study, we assume the parameters are current dependent-the examined object is non-linear. Typical identification procedures rely on fitting the model to measurements as demonstrated by Chwastek and Szczygłowski (2008) or Gołdasz et al. (2019), usually by means of least-squares quality metrics. The Duhem model parameters are not directly related to the material's FIGURE 5 | Block diagram of the hysteretic model (gas force not shown); physical properties, e.g., remanence, coercivity, contrary to the famous Jiles-Atherton model, for instance. In the anhysteretic case (a = 0) the inductor model reduces to that of a non-linear one without hysteresis, i.e., λ = d(i)i. Then, an initial estimate of the parameter d at a given current level can be provided from coil inductance calculations using lumped parameter models or more advanced finite-element techniques, for example. Once the parameter d is fixed, the parameters a and b 1 can be played with to arrive at the correct hysteresis width. Finally, the rate tuning  parameter b 2 of the hyperbolic tangent function can be adjusted to match the λ − i curve's slope. The particular form of the shape function was preferred for consistency and ease of use. The model parameter set can then be used for varying the hysteresis width and shape. In Figure 3, the authors reveal the impact that each parameter of the model has on the λ − i normalized loop shape; the arrow in each figure shows the direction of particular parameter increase. In the normalized plots, the flux and current variation range is from −1 to 1.
In this context, the Duhem operator is a non-linear system model driven by the coil current i. Its output is the flux linkage λ(t). In the analyzed form, the Duhem model has the advantage over other hysteretic models, e.g., the Bouc-Wen model, for being less complex, as it accepts half of the parameters (that the standard Bouc-Wen model requires) for successful operation as shown in Gołdasz et al. (2018b).

Phenomenological Model of an MR Damper
Next, we examine the simplest first-order Maxwell model incorporating a non-linear dashpot in series with a spring. The model is assumed to ignore higher-order dynamics due to mass effects. Only the compressibility of fluid chambers is taken into account as main contributor to the damper hysteresis. The nonlinear Maxwell model is functionally identical to the model examined by Gołdasz and Alexandridis (2012) as well as Simms et al. (2004). As a partial proof, let us then consider the cylinder housing with an MR valve located in the piston as in Figure 4. The piston is driven by the prescribed displacement x(t), thus forcing the fluid through the valve at the volumetric rate Q. The pressure drop across the piston is p = P 1 − P 2 . Considering the variation of pressures in each fluid chamber, we yield where β-bulk modulus. Assuming the volumes V 1 = V 2 = V and transforming the above system of equations gives where A is the piston cross-section area. In general, the pressure drop across the MR valve and the flow rate are related by a generic non-linear function Q = f ( p, λ). By way of simplification, we further assume Q = p/R(λ) and R = R(λ) is field-dependent hydraulic resistance (which further reduces to a (non)linear hydraulic resistance in the absence of flux). Substitution then gives or multiplying by the piston cross-section area to get the output force F yields This can be rewritten into a more generic form as follows: where ζ is the field-dependent time constant. As an example, we represent the force shape function F d to be where the force component F v represents viscous losses proportional to the velocityẋ, R v denotes hydraulic resistance due to viscous losses, F λ copies the field-dependent force due to the MR effect. R λ = R(λ)A 2 is a link between the damping force model and the inductor model. The model may incorporate friction force as well as force offset due to the presence of the pressurized accumulator. Moreover, c defines the force increase rate against the velocityẋ. The model's simple form, as shown in Figure 5, is particularly convenient in model/parameter identification problems; the parameter k d represents all compliances present in the damper. Note that if fixed flux input is assumed, then the effect of force evolution against the velocity may be easily analyzed independently of the other excitations. Under time-varying current excitations, the flux build-up is copied by the inductor component, including a relevant hysteretic operator (see Equation 1), whereas the force change against the flux (as well as the velocityẋ and the displacement x) is taken care of by Equation (8).
As shown, the model parameters can be deduced from the damper's geometry and material properties. However, due to unknown internals of the tested commercial MR damper, we estimate them from damping force measurements.

TEST INPUTS AND HARDWARE
The model requires measurements of quantities from which the investigated models can be identified and their parameters estimated, namely, flux and force. Although both could be estimated simultaneously, and all measurements were carried out at the same time, we decided to rely on the sequential approach. First, to acquire the λ − i relationship for the damper, magnetic flux was extracted from voltage and current data via integration (and high-pass filtering for drift removal). The measurements were accomplished over a 20 s time span and by sampling the data at the frequency of 1 kHz. The data acquisition system was under the control of an AD/DA I/O board. The tested unit was a flow-mode linear commercial damper by BWI Group for a passenger car featuring an MR valve with one annular gap and one coil assembly. The long-stroke damper's coil resistance is appr. 1.1 incl. electrical connections, and the device can be operated up to 5 A. More detailed information on the tested damper incl. transient response, frequency response and other characteristics can be found in Jastrzębski and Sapiński (2017), and the reader should refer to these for further details.
The damper was tested on a mechanical MTS810 shaker, and the displacement was acquired from the stroker's internal sensor, see Figure 6. In our experiments, the damper was subjected to constant velocity (triangular displacement) inputs using the above 1 Hz sine wave current input superimposed on the displacement profile x(t) as in Figure 7A. The peakto-peak displacement was 150 mm, and the triangular displacement wave frequency was set to be 0.5 Hz. During the experiments, the input voltage was adjusted to result in peak coil currents I = {1, 1.5, 2, 2.5, 3, and 3.5} A. The magnetic flux time histories were then reconstructed as already mentioned.

RESULTS AND DISCUSSION
Below, we show the obtained experimental results and then apply the proposed model for predicting the dampers characteristics.
Damper measurement results are highlighted in Figures 7-13. In particular, Figure 7A contains an exemplary time history of the excitation inputs (displacement, current) plotted over a selected 4 s time span. In Figure 7B, we plot time histories of the damping force corresponding to the abovementioned current peak levels. The results are complemented by the plots of the damping force vs. the input displacement in Figure 7C and the damping force vs. the measured coil current in Figure 7D. In the plots, the damping force is offset by the gas force due to the pressurized accumulator present in the monotube damper; the commercial damper was gas-charged at P 0 = 2.6 MPa. The gas pressure was measured at mid-stroke. The gas volume was estimated by pushing the piston rod downwards from a fully extended position (full rebound) to the fully collapsed one (full compression). The other bias force which can be observed in the data is friction (determined directly from the experimental data-F r ≈ 70 N). The viscous damping coefficient was estimated by plotting the force vs. velocity at zero current condition. Furthermore, the flux linkage information was acquired by post-processing voltage and current as already mentioned; the flux integration procedure was applied for predicting the magnetic hysteresis of an MR damper and explained in detail in Gołdasz et al. (2019). There are several unavoidable issues with flux integration. First, the flux's initial value was unknown; that issue was solved by demagnetizing the damper after each measurement. Second, as the damper is stroked and/or the current is applied to the coil, the internal temperature increases. As a result, the coil resistance varies with temperature; the issue could be solved by, for example, using a sensory coil wrapped around the core, though this was not possible due to a lack of access to the internals of the damper. Third, measurement noise and integration errors are further accumulated in the process resulting in the flux linkage signal drift. That issue was at least partially solved by filtering to remove the signal bias and trend/drift. These issues may have influenced the outcome of the study.
When the peak current increases from 1 to 3.5 A, on the inductor part, the parameter a was determined to decrease from the maximum value of 3 down to 1.4. On the other hand, the parameter b 1 increased from 0.03 to 0.08 (3.5 A) with the increasing current. At the same time, the rate tuning parameter b 2 varied from 0.7 to 0.25, and the parameter d was found to vary from 0.015 to 0.004. The time constant ζ was determined to be nearly constant -0.015 to 0.018 s, and the off-state damping coefficient R v = 120 Ns/m, and c = 15 s/m. The parameter R(λ) was found to be equal to 18 kN/Wb (1 A), increasing with the current up to 28.8 kN/Wb (3.5 A). Based on the observations of plots in Figures 8-13, it seems that the evolution of damping force against the current and/or magnetic flux can be well-studied with the proposed dual hysteretic approach. At all current levels above 1 A, the model is capable of providing good quality performance in predicting both the F − x loops, the F − i relationship, as well as the F − λ plots. The force behavior at transition points as well as the hysteresis width and shape have been wellcaptured with the proposed approach. At the lowest current case (1 A), only the F − λ predictions are not acceptable. The flux as well as force levels are well-predicted; however, the F − λ loop's width is poorly captured, and the transition points are in the wrong quadrants of the F − λ system of coordinates which can be observed also in the F − x loop of the same figure. The model's poor prediction at the lowest peak current level is likely to be due to the above mentioned flux estimation problems.

SUMMARY AND CONCLUSIONS
Hysteresis has been a well-known phenomenon influencing the force/torque output of MR dampers. Prior art on this topic usually included studying the total or combined hysteresis by considering the relationship of the damping force or torque vs. displacement/velocity and current, thus ignoring the distinct nature of magnetic hysteresis and the hydraulic/mechanical hysteresis. The first one is an inherent property of ferromagnetic materials forming the damper's electromagnetic circuit. It is of a different nature to the mechanical (hydraulic) hysteresis of the devices since it does not vanish in the DC limit (as the excitation frequency approaches zero). Understanding the contribution of the two mechanisms is then vital in developing a good quality model. In MR dampers, magnetic hysteresis not only degrades any current-based control scheme performance but also reduces the actuator's dynamic range.
To distinguish between the two mechanisms, the present study extends the concept of a non-linear inductor with magnetic hysteresis (which is then linked to the damping force model). The inductor concept (based on the Duhem hysteretic operator) is employed for copying the λ − i characteristics of the solenoid. In the presented form, it simply captures the average flux variation in the structure. It requires estimating the model parameter values from the damper measurements. As such, it may be only used during a control algorithm development stage and is not suitable, for example, for damper sizing studies. The lumped parameter form may, however, be convenient for modelbased control studies. Replacing the Duhem operator with the Jiles-Atherton (whose parameters are linked to magnetization characteristics features) model or the vector hysteresis modeling approach would, however, make it suitable for solving such problems at early stages of the development process, though often at the expense of higher computing cost.
The mechanical mechanisms contribution to the force output are then captured with a phenomenological Maxwell-based model of the damper featuring a spring in series with a non-linear dashpot. As demonstrated, the presented procedure allows independent interpretations of the contribution of each analyzed mechanism.
The phenomenological approach undertaken by the authors employs a posteriori models whose parameters can be computed from experimental data. Only control studies can thus be supported with it. That is in contrast with the multiphysics technique presented by Kubik and Gołdasz (2019) where the authors showed a hybrid finite-element/lumped parameter approach that can be employed at a design stage.
To further summarize and clarify, in the present study we highlighted an approach targeted toward separation of the hysteretic output of an MR flow-mode damper into two distinct hysteretic operators, namely, mechanical/hydraulic hysteresis and magnetic flux hysteresis. With the proposed approach each mechanism can be studied independently. The presented technique relies on the flux sensorless estimation technique which is particularly convenient if no access to the internal components of a damper can be gained. The damper model is phenomenological and comprises parameters that can be extracted from measurements. Moreover, it allows for the independent analyses of each hysteresis mechanisms with any existing hysteretic model. Both the Duhem model and the derived Maxwell type model were used here for illustration purposes.
Finally, application of the model in more complex transient studies requires adopting a more advanced inductor concept so that the effects of eddy currents and the hysteresis dependence on the excitation input frequency are well-captured. The advanced inductor concept is a subject of on-going study. Moreover, work on a state estimation technique based on the non-linear Kalman filter is in progress. Implementing it is crucial given the flux integration problems mentioned in the previous section.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JG conceived the model and carried out simulations. BS and ŁJ planned and carried out experiments. MK processed data, participated in model verification and contributed to the interpretation of the results. All authors provided critical feedback and helped shape the research, modeling results and manuscript.