Dynamic Model of MR Dampers Based on a Hysteretic Magnetic Circuit

As a key to understand dynamic performances ofMRdampers, a comprehensive dynamicmagnetic circuitmodel is proposed in this work on the basis of Ampere’s and Gauss’s laws. It takes into account not only the magnetic saturation, which many existing studies have focused on, but also the magnetic hysteresis and eddy currents in a MR damper. The hysteresis of steel parts of MR dampers is described by Jiles-Atherton (J-A) models, and the eddy current is included based on the field separation. Compared with the FEM results, the proposedmodel is validated in lowand high-frequency studies for the predictions of the magnetic saturation, the hysteresis, and the effect of eddy currents. A simple multiphysics model is developed to demonstrate how to combine the proposed magnetic circuit model with the commonly used Bingham fluid model. The damping force in the high-frequency case obviously lags behind the coil current, which exhibits a hysteresis loop in the current-force plane. The lag of damping force even exists in a low-frequency varying magnetic field and becomes more severe in the presence of eddy currents.


Introduction
The dynamic magnetic behavior of MR dampers cannot be fully understood yet due to the lack of a comprehensive electromagnetic model, in which the magnetic saturation, the ferromagnetic hysteresis, and the eddy current in the steel parts of MR dampers should be fully considered under a time-varying magnetic field.
The finite element method (FEM) has been the main tool for the design and analysis of the magnetic circuit of MR dampers, with many previous studies focusing on the magnetic saturation analysis [1][2][3][4][5].In contrast with FEM, the electric circuit analogy method based on Ampere's and Gauss's laws is much easier to implement, flexible, and less time consuming and proved to be particularly effective in the optimization of magnetic circuits of MR devices [6,7].
Due to the hysteresis of ferromagnetic materials, magnetic induction follows different paths with increasing and decreasing magnetic field, and this ambiguous characteristic further complicates the predicting of the performances of MR devices.So far, there have been limited studies on the magnetic hysteresis modeling of MR devices.And not surprisingly, few commercial FEM software tools support the hysteresis modeling except that ANSYS Maxwell5 starts to include such capability until its recent release [8].Several research efforts have been made to integrate numerous hysteresis models such as Preisach model [9] and J-A model [10] into FEM for a general electromagnetic modeling [11,12].In this way, Je ¸dryczka et al. [13] showed that the residual magnetic flux density had a significant influence on the torque of a MR clutch, resulting in an increase of the friction torque in the disengagement state.An and Kwon [14] developed a nonlinear model of MR actuator using the Hodgdon hysteresis model [15] for the magnetic circuit and Bingham model for the MR fluid.They successfully simulated the electric current-torque hysteresis.Instead of starting from the hysteresis of materials, Deur et al. [16] described the electric current-torque hysteresis in a phenomenological way based on the experimental data and showed that the modeling accuracy was significantly improved after including the effect of hysteresis.
From an energy point of view, the effect of eddy current can be taken into account by incorporating in a physical model the classical (macroscopic) and excess (microscopic) eddy current losses, with the total loss () decomposed into hysteresis loss ( ℎ ), classical loss (  ), and excess loss (  ); that is,  =  ℎ +   +   . (1) The hysteresis loss ( ℎ ) increases linearly with frequency while the eddy current loss (  ) increases with the frequency squared.There is usually a discrepancy between the measured loss and the loss expected from the sum of hysteresis and eddy current losses, and this is usually referred to as the anomalous/excess loss (  ) [17].Thus, under a high-frequency changing magnetic field, the eddy current in an electrically conducting material, particularly in a solid electrically conducting material (e.g., the steel piston in a MR damper), is significant and delays the magnetic responses of devices.Replacing the high conductivity aluminum alloy by an insulating material, Takesue et al. [18] sharply improved the magnetic response of a magnetic actuator by reducing the eddy current.The existence of eddy current changed the inductance dynamics of a 200 kN largescale MR damper, as shown by Jiang and Christenson [19], and such effect should be considered in order to fully capture the time-varying effective current that magnetizes the MR fluid.Nam and Park [20] qualitatively analyzed the magnetic response characteristic of a MR damper and concluded that the magnetic reluctance should be increased as much as possible in order to improve the magnetic response.Under the assumptions of ignorance of the hysteresis, Farjoud and Bagherpour [21] recently performed a more detailed analysis on the eddy current in MR devices.The flux leakage and magnetic saturation are taken into account by integrating the static FEM into the whole Simulink5 model as an function.However, it is worth noting that the hysteresis loss in steel sheet accounts for more than 70% of the total loss even under a 50 Hz high-frequency varying magnetic field [22], so the assumption of negligible hysteresis effect should be used carefully.
In this paper, a general dynamic magnetic circuit model of MR dampers will be proposed.It is expected to have the following advantages compared to the available studies.
(1) The proposed model is based on Ampere's and Gauss's laws and is essentially a simple lumped parameter model, with high computational efficiency achieved.All the model parameters have clear physical meanings, and the results are easy to interpret.Additionally, the proposed magnetic circuit model can be easily extended by coupling it with mechanical models (e.g., Bingham model for MR fluid) to realize a multiphysics analysis.
(2) By using J-A models to describe the ferromagnetic materials (e.g., the piston and house cylinder of a MR damper), the proposed model is a fully dynamic magnetic circuit model with the magnetic saturation, hysteresis, and eddy current all included.This makes it possible to conduct more comprehensive studies on the dynamic performances of a MR damper.Moreover, the proposed model can be also applied to the permanent magnets based hybrid-magneticcircuit [23], self-sensing [24], and energy-harvesting [25] MR dampers which receive more and more research interests, because permanent magnets [26,27] can also be accurately characterized by J-A models.

J-A Model for Ferromagnetic Materials
The J-A model [10] is based on physical premises concerning domain wall movement in soft magnetic materials.It is interesting from both the theoretical and practical points of view.It offers insights into the details of the magnetization processes.The form of J-A model has evolved over the last few decades, and more and more physical phenomena have been taken into account, such as the effects from mechanical stress, anisotropy, temperature, eddy current, and dependency of frequency.
J-A model consists of algebra-differential equations in which the magnetization () is either calculated from the magnetic induction () or from the magnetic field ().These two forms, () and (), are equivalent and the model parameters share the same set of values.The () form is more convenient to characterize a material, while the () form, referred to as the inverse J-A model, is preferred in a magnetic circuit analysis in which magnetic induction () is often used as the independent variable.Hereafter the J-A model refers to the () form unless otherwise stated.
As shown by Zirka et al. [28] the dynamic J-A model can be expressed in a convenient form after some algebraic manipulations: where and the anhysteretic magnetization is given by the Langevin function the effective field the dynamic field and the directional parameter  = sign(/).The parameter   is introduced to eliminate the unphysical negative susceptibility and calculated by Model parameters,   , , , , ,   ,   , are obtained from measurements, with their physical meanings listed in Table 1.The corresponding inverse J-A model showed a similar form [28], except that and /, rather than /, is used in parameters  and   . 0 is the permeability of free space.In addition, efforts also have been devoted to improving the accuracy and convergence of J-A model.For example, Zirka et al. [28] suggested that   =  1 (1 +  2  2 ) in order to refine the accuracy of J-A model under highfrequency magnetic fields.Wilson et al. [29] proposed that the convergence can be improved by making parameter  dependent on the applied field in a form of Gaussian function: where  0 is the default value of the  parameter in the original J-A model,  is the applied field, and  is the standard deviation of the Gaussian function.Before using J-A model in the magnetic circuit of MR dampers, it is safe to validate its effectiveness of characterizing magnetic hysteresis.Kis and Iványi [30] tested the magnetic hysteresis of C19 structural steel by using a twocoil measurement system, in which a prescribed magnetic field was generated by the primary coil and the magnetic flux density was calculated from the induced voltage on the secondary coil.More details of the measurement can be found in the references therein.Here, their experimental data is used for investigating the effectiveness of J-A model, and a simple least-squares method is used to determine the model parameters, with the block scheme of identification shown in Figure 1.The ranges of parameters proposed by Chwastek and Szczygłowski [31] in Table 1 are used here, and the initial parameters (x 0 in Figure 1) are chosen as the medians of the ranges.
The following criterion function is used to generate an error by comparing the experimental and the simulated hysteresis: where  is the number of measurement points,  exp, and  sim, are the th measured and simulated magnetic induction, and  max is the maximum of the measured magnetic induction.
The detailed identification process of the model parameters is as below: (a) Input the initial guesses of model parameters.
(b) Using the adopted parameter values, solve the J-A model and obtain the simulated magnetic flux densities ( sim, ) at different magnetic fields points specified by measurements.
(c) Compare the simulated results ( sim, ) with the measured ones ( exp, ), and calculate the error using (10).
(d) If the error is larger than the tolerance (), the nonlinear least-squares method automatically modifies the input values of model parameters.If the error is below the tolerance, the corresponding parameter values are accepted, and the identification process is finished.Figure 2: The simulated and measured [30] hysteresis curves for different frequencies.
After about 50 iterations, the parameters are identified as As shown in Figure 2, the model results agree well with the experimental data, so J-A model is effective and will be used in the following sections to describe the magnetic hysteresis of the steel parts in MR dampers.

A Dynamic Magnetic Circuit Model of MR Dampers
As shown in Figure 3, the magnetic circuit of a MR damper is composed of a gap and three steel parts (denoted by A, B, and C).An independent current source is considered here because current drivers are preferred for faster magnetic responses [32].
For simplicity, a static J-A model in which the dynamic field (  ) vanishes will be used to describe the steel parts of MR dampers, and the effect of eddy current will be included later in the magnetic circuit modeling.Another main reason for choosing a static J-A model is that the hysteresis algorithm in Maxwell software, which will be compared with the magnetic circuit model proposed in this paper, is still a black box with little information available.However, such a temporary simplification will not undermine the generality of the proposed model, because the dynamic hysteresis model would be recovered by simply restoring the dynamic field (  ).
The steels are described by the following three inverse J-A models which constitute three differential algebraic equations (DAEs) in terms of the magnetizations (  ): g Similarly, the magnetic saturation/hysteresis of MR fluid can be described by introducing another J-A model.However, the magnetic permeability of MR fluid is very low in comparison with that of steels, so it is neglected here and regarded as air.
Without loss of generality, the piston and house cylinder are assumed to be made of the same steel, so the three J-A models share the same set of parameters.
If the effects of magnetic flux fringing and leakage are ignored, according to Gauss's law, the magnetic inductions in the different parts of the circuit are related by where   and   are the average magnetic inductions in the air gap and steels, Φ is the magnetic flux in the circuit loop as shown in Figure 3,  1 and  2 are the cross-sectional areas of the piston core and the house cylinder,  3 is the average lateral area of the piston head (i.e., the lateral area at position III in Figure 3), and   is the average lateral area of the annular gap.These four areas are given by with the meanings of all the variables shown in Figure 3.The flux leakage, fringing flux, and the nonuniformity of magnetic field/induction, as shown in Figure 4, are still challenging to be accurately considered, and their combined effect is reflected here by introducing parameter   for the gap and  , for the steels; that is, The governing equation of the whole magnetic circuit can be obtained by Ampere' law as where  tot is the total current in the source,   is the magnetic circuit length of the steel parts,  1 =  2 = (  +   )/2,  3 =   −   ,  is the gap width.The magnetic field in different steel parts (  ) are coupled to the corresponding J-A model by And the magnetic field in the gap is simply calculated by In a magnetic circuit analysis, the effect of either classical or excess eddy current can be taken into account based on the following field separation [28,33,34]: where the total magnetic field ( tot ) is the sum of the hysteretic magnetic field (), the dynamic field due to classical eddy current (  ), and the dynamic field concerning excess eddy current (  ).In the cases of low frequencies and/or magnetic circuits with small cross-sections, the magnetic field and induction can be assumed to be uniform.Consequently, the field separation is equivalent to the aforementioned loss separation (1) and the dynamic fields can be estimated by [35,36] where  is the resistivity,  is the cross-sectional dimension (thickness for laminations, diameter for cylinders and spheres), and  is a geometrical factor which varies form  = 6 in laminations to  = 16 in cylinders and  = 20 in spheres.
In order to keep the generality of the proposed model, the general form ( 6) is used for the dynamic magnetic field; that is, Only the dynamic field due to classical eddy current will be considered for simplicity.Accordingly, ( 15) is extended to include the effect of classical eddy current: Finally, ( 11) and ( 21), expressed in terms of four independent variables ( 1 ,  2 ,  3 , Φ), constitute a complete magnetic circuit model of MR dampers.This model is essentially a set of DAEs and can be numerically solved by the Matlab built-in ODE/DAE solver, ode15i, with the relative error tolerance and absolute tolerance, respectively, set to 1e − 5 and 1e − 6.

Model Validation and Analysis
The effectiveness of the proposed magnetic circuit model is investigated for describing the low-and high-frequency magnetic performances of a MR damper with the main structural parameters chosen as   = 34.6 mm,   = 30 mm,   = 20 mm,  = 0.6 mm,   = 25 mm, and   = 30 mm.The coil has 284 turns.The current driver supplies a sinusoidal current with an amplitude of 2 A and a frequency of 0.2 Hz in the low-frequency study and 10 Hz in the high-frequency study.
For the validation of the proposed magnetic circuit model, it should be noted that, although comparison with experimental data is always the most direct way, it is often (e.g., in the predesign stage of a MR damper) not a preferred way for validation because the magnetic response of a MR damper, indispensable for identifying the parameters of a lumped parameter model, is not available until the damper is finally fabricated and then tested.In this case, a FEM simulation, although coming with drawbacks, such as being time consuming and hard-to-interpret, is often a practical and reliable alternative.
More importantly, although validation by an indirect test of damping forces is always much easier, it should not be a preferred way for validating a magnetic circuit model because the damping force itself is greatly dependent on many factors such as the compressibility of MR fluid and the compliance of instrument system.Moreover, the effect of stress on magnetic hysteresis is often unavoidable, which further makes the damping force an unsuitable indicator to validate a magnetic circuit model.
Since the FEM simulation will be used to validate the proposed model, the FEM results will act as the "experimental data" in (10), and the similar error will be calculated as where  fem, is the th magnetic induction obtained by FEM,  lum, is the th magnetic induction obtained by the proposed model,  FEM,max is the maximum FEM magnetic induction, and  is the number FEM results.It should be noted that the summation index, , implicitly includes all the average FEM results at three parts (lines I, II, and III in Figure 3).ANSYS Maxwell employs a vector hysteresis model to predict the hysteresis loops and losses for soft and hard magnetic materials and permanent magnets.Hysteresis is automatically activated by properly setting the material data, and the detailed material settings in ANSYS Maxwell are given below: (a) Specify the steel parts of the damper as materials with nonlinear relative permeabilities and edit the  curves.
(b) Input the descending branch of the saturation hysteresis loop.It should be noted that the first point and the last point must mirror each other; otherwise, the curve is considered invalid.
(c) After finishing the  curve input, a coercivity value is automatically calculated and listed.
(d) Set all three components of the "Unit vector" to zero; otherwise, it acts like a permanent magnet.

Low-Frequency Study: Magnetic Saturation and Hysteresis.
Under a low-frequency (0.2 Hz) current excitation, the effect of eddy current is insignificant, and the magnetic hysteresis/saturation predicting capability of the proposed model is examined in this case.Accordingly, the three parameters,  ,1 ,  ,2 ,  ,3 , related to the dynamic field vanish.
The magnetic responses in the three steel parts are simulated by the proposed model with the parameters identified as  = 1.1440 × 10 −3 ,   = 1.1527 × 10 6 ,  = 6.3750 × 10 2 ,  = 2.2428 × 10 −1 ,  = 6.3750 × 10 2 ,   = 1.3284, and  ,1 =  ,2 =  ,3 = 1.0000.The three average cross-sectional magnetic responses at positions I, II, and III (Figure 3) are also obtained in Maxwell (v16.0)for comparison.As shown in Figure 5, the results predicted by the proposed model are in a satisfactory agreement with the FEM results, so the proposed magnetic circuit model is reliable for describing the static magnetic hysteretic behavior of MR dampers.
It can be seen from values of the parameters ( 1 =  2 =  3 = 1.0000) that the combined effects of magnetic fringing, leakage, and nonuniformity of magnetic field/induction are negligible in all the three steel parts and will be ignored in the following high-frequency study.

High-Frequency Study: Effect of Eddy Current.
The highfrequency study is performed under a sinusoidal current of 10 Hz, and the maximum of the eddy current density is one order of magnitude larger than that under a current of 0.2 Hz, as shown in Figure 6.
Because the five J-A material parameters have been already identified in the low-frequency study, only four parameters ( ,1 ,  ,2 ,  ,3 ,   ) are left to be determined.Using  an identification method similar to that shown in the above section, they are obtained as  ,1 = 1.6650 × 10 2 ,  ,2 = 2.7785 × 10 1 ,  ,2 = 4.8578 × 10 1 , and   = 1.0593.As shown in Figure 7, the predicted results in the steel parts show an overall satisfactory agreement with the FEM results, although a larger discrepancy observed as compared to the results in the low-frequency study.The numerical errors in low-and high-frequency studies, calculated by (22), are compared in Figure 8.The agreement between FEM and the proposed model results in the low-frequency study is much better than that in the highfrequency study, with the largest error occurring in steel part A.
The larger discrepancy or error, especially that in steel part A, is believed to be the consequence of the larger nonuniformity of magnetic induction due to the effect of eddy current.The hysteresis curves along lines I, II, and III (in Figure 3) are shown in Figure 9 with obvious nonuniformities of magnetic fields observed under high frequency.Compared to the nonuniformities of magnetic fields in parts B and C, the nonuniformity in part A is much larger, so using average values in the lumped parameter model produces larger error.

Shock and Vibration
Moreover, the comparison between the two groups of hysteresis curves in Figures 8 and 9 shows that the magnetic induction, which directly determines the yields stress of MR fluid, is apparently reduced under high-frequency changing fields due to the effect of eddy currents.

A Demonstration of How to Use the Proposed Magnetic Circuit Model
For the purpose of the demonstration of how to use the proposed magnetic circuit model, a simple multiphysics model will be developed in this section.The flow field in a MR damper can be controlled by the magnetic field due to changes of the shear yield-strength of MR fluid.In contrast, the magnetic field is generally assumed to be independent of the fluid flow field, so the analyses of these two fields can be conducted separately in series.The experimental yield stresses (  ) of the MR fluid are shown in Figure 10 under different magnetic inductions (), and their relationship, which is required for the damping force modeling, can be obtained by a least square curve fitting method as As shown in Figure 10, the fitted curve agrees well with the experimental data.As shown in Figure 11, after validating the proposed magnetic circuit model, the damping force () and the flow field () of a MR damper can be analyzed by substituting the yieldstrength history (obtained from the magnetic circuit model) into any extensively studied fluid model such as the Bingham model, the Herschel-Bulkley model, and the biviscous model.The simple Bingham fluid model is adopted here for the demonstration purpose.Since the proposed magnetic circuit model has been validated in previous sections, the results from the combinations of these two reliable models should be also reliable.
Describing the flow of MR fluid by the Bingham flow equations (i.e., combination of the Navier-Stokes equations and the Bingham constitutive equation) in a rectangular duct, the damping force () of a MR damper can be related to the piston velocity (V  ) by the following mechanical model [37]: with the dimensionless pressure gradient (P) and yield stress (T) defined by where  is the mean circumference of the annular flow path, Δ is the pressure drop responding to the volumetric flow rate  = V    ,   is the effective area of the piston head, the damping force  = Δ  , and  is the postyield plastic viscosity and taken as 1.3 Pa⋅s.
Using the time history of the shear yield stress (23) as the bridge connecting the proposed magnetic circuit model and the Bingham fluid based mechanical model (24), the damping forces corresponding to the low-and high-frequency studies in the above section are examined, with the piston velocity chosen as 10 mm/s.
As shown in Figure 12, compared to the results of the low-frequency study, an obvious damping force lag with respect to the coil current is observed in the high-frequency study due to the effect of eddy currents.Such lag exhibits a hysteresis loop in the current-force plane, as shown in Figure 13, which even exists in a low-frequency varying magnetic field because of the magnetic hysteresis.This hysteresis behavior produces a variation in the damping force between an increasing and decreasing coil current.Thus, in a dynamic environment such as oscillatory current excitation, uncertainty is introduced in any calculated torque.Worse still, the eddy currents make such hysteresis more severe.As expected, an apparent decrease in the maximum damping force is also observed in the high-frequency study because of the decreasing maximum magnetic induction as shown in Figures 5 and 7.
For better understanding of dynamic performances of MR dampers, it is often desirable to gain an insight into the flow analysis of MR fluid in the working gap, and this  becomes straightforward once the pressure drop (Δ) has been obtained by solving (24).The flow velocity () can be directly calculated by [ where  is the gap coordinate and  is the plug thickness, as shown in Figure 14.
The velocity profiles corresponding to the low-and highfrequency studies are shown in Figure 15, in which the magnitudes of the coil currents are indicated by the surface colors.The parabolic velocity profiles occurring at the instants of zero current and flat velocity profiles at other instants appear alternately.Typically for a flow of MR fluids, the larger the currents are, the flatter the flow velocity profiles are.

Conclusion
In this research work, a novel dynamic magnetic circuit model of MR dampers is proposed on the basis of Ampere's

Figure 1 :
Figure 1: Block diagram of the identifications of J-A model parameters.

Figure 4 :
Figure 4: The fringing flux, leakage flux, and nonuniform magnetic field near the working gap of a MR damper.

Figure 8 :
Figure 8: Comparison of numerical errors under low-and highfrequency studies.

Figure 9 :
Figure 9: Hysteresis curves at points of lines (in Figure 3) under low-and high-frequency current excitation (nonlinear residual of FEM solver: 1e − 5).

Table 1 :
Physical meanings and estimated value ranges of J-A model parameters.