# Multiphysics Model of an MR Damper including Magnetic Hysteresis

Authors: M. Kubík; J. Goldasz
Journal: Shock and Vibration
 (2019)

---

## Abstract

Hysteresis is one of key factors influencing the output of magnetorheological (MR) actuators. The actuators reveal two primary sources of hysteresis. The hydro(mechanical) hysteresis can be related to flow dynamics mechanisms and is frequency- or rate-dependent. For comparison, the magnetic hysteresis is an inherent property of ferromagnetic materials forming the magnetic circuit of the actuators. The need for a good quality hysteresis model has been early recognized in studies on MR actuators; however, few studies have provided models which could be used in the design stage. In the paper we reveal a hybrid multiphysics model of a flow-mode MR actuator which could be used for that purpose. The model relies on the information which can be extracted primarily from material datasheets and engineering drawings. We reveal key details of the model and then verify it against measured data. Finally, we employ it in a parameter sensitivity study to examine the influence of magnetic hysteresis and other relevant factors on the output of the actuator.


---

## Body

## 1. Introduction

Magnetorheological (MR) dampers are fairly well-known devices utilizing MR fluids which, when subjected to magnetic stimuli of sufficient strength, generate yield stress [1]. So far, the unique technology has been commercialized in semiactive passenger vehicle suspensions, powertrain mounts [2], or optical finishing [3]. Low power consumption, fast and reversible responses, and high dynamic range have made the devices attractive for use in vibration control systems in particular [4]. As MR dampers are generally operated in real-time control systems, their dynamic performance is equally important as or more important than their steady-state characteristics. Steady-state characteristics only provide the evidence of an actuator or a damper meeting the required force/torque range (or force/torque) targets. Their dynamic behaviour needs to be quantified at the same time if it is used in a real-life control process. Thus, understanding the contributions of various factors complicating the force or torque build-up dynamic process is critical for the development of a realistic application. Briefly, with MR actuators, there is ample evidence of several factors complicating the force/torque generation process, namely, mechanical/hydraulic hysteresis, magnetic hysteresis, control circuit dynamics such as eddy currents, driver dynamics, temperature, flow losses, friction, and nonlinear relationship between the material’s yield stress and the induced flux [5–7]. These factors influence the device’s ability to generate the output force/torque and need to be accounted, for instance, for in the control algorithm development process.In this study, we pay particular attention to modeling the damper’s hysteretic behaviour. In general, MR devices reveal two primary sources of hysteresis. The (hydro-)mechanical hysteresis can be related to the damped dynamics of a heavy slug of MR fluid (MRF) bouncing against compliant columns of MRF in fluid chambers. The effect is rate- or frequency-dependent, and its magnitude varies with the current applied, too. It disappears as the mechanical excitation frequency approaches zero [8]. The magnetic hysteresis is different. It is present in all electromagnetic devices, e.g., electromagnetic solenoids [9], motors [10], and magnetorheological actuators [11]. First of all, it is the inherent property of ferromagnetic materials forming the magnetic circuit of the MR valve; the hysteresis of carbonyl iron- (CIP-) based MFRs is virtually nonexistent [12]. Next, it does not vanish as the inducing current frequency approaches DC limit. Also, temperature, load history, and mechanical stresses have a negative influence on the hysteresis and magnetization characteristics of ferromagnetic materials [13]. For instance, it is a common practice to subject machined ferromagnetic components to heat treatment for internal stress and hysteresis as well as coercive force reduction.The need for a good quality hysteresis model has been early recognized in studies on MR actuators, and the reader should refer, e.g., to Zheng et al. [14] for a review of suitable phenomenological models as well as to the well-known study of Spencer et al. [15]. In general, the posteriori parametric models were obtained by examining the force-position and force-velocity relationships by fitting by model response to the actuator’s output. Such models are suitable for control studies only. In the device’s development process, other approaches are required. In that aspect, many MR-related research studies neglected the particular contributor’s presence. The topic, however, has been well identified in the field of conventional solenoid actuators where various models were developed to copy the hysteretic behaviour of the actuators. For instance, Mayergoyz [16] applied the Preisach model to model the hysteretic behaviour of a solenoid actuator. In the Preisach model, the hysteresis is the sum of elementary hysteresis loops. Next, Coleman and Hodgdon [17] developed a first-order differential equation that links the field strength H and the flux density B. One model whose parameters can be related to physical properties of ferromagnetic materials is the Jiles–Atherton (J–A) model [18, 19]. The model was extended to include both the impact of eddy currents and temperature on hysteresis and magnetisation curves [20, 21]. All of the above models can be vectorized. Tellinen [22] proposed a simple scalar model for handling the hysteresis based on the limiting hysteresis loop from physical measurements of ferromagnetic materials. With MR actuators, however, although the significance of a good quality hysteretic model has been recognized early, the topic does not seem to have deserved enough attention. Significant contributions include Han et al. [23] who examined the field dependent hysteresis of ER fluids. The authors used the familiar Preisach approach. Moreover, Han et al. [7] used the Preisach model to identify the hysteresis of an MR fluid. Yadmellat and Kermani developed a model of an MR clutch in which the developed hysteresis model was assessed against the Preisach operator [24]. For comparison, in their early study using the Coleman-Hodgdon model, An and Kwon modelled the hysteretic behaviour of an MR clutch, and by examining the torque-current loops showed that the hysteresis is an important contributor to the device’s output [25]. The model parameters were identified from physical measurements (of magnetisation characteristics) of the materials forming the magnetic circuit of the clutch. Next, Jędryczka et al. presented a finite-element (FE) model of an MR clutch based on the J-A approach [26]. Moreover, Guo et al. presented a transient multidomain model of a flow-mode damper based on the J-A approach and then verified it against the novel FE vector hysteresis technique [27]. The inverse J-A model was recently examined by Zheng et al. [14] to copy the transient behaviour of an MR flow-mode damper. Goldasz et al. [28] proposed an extension of the Bouc–Wen model in an attempt to separate the magnetic hysteresis from the mechanical one. The authors proposed a simple lumped parameter model of the actuator including a hysteretic operator. The model was verified against selected sinusoidal AC excitation inputs and provided acceptable accuracy for practical purposes. Still, when compared to the vast number of research studies using parametric phenomenological hysteretic models, the topic does not seem intensively studied as already mentioned. That may be due to few existing comprehensive electromagnetic models of such actuators.Another aspect is dynamics. Clearly, the insight into the dynamics of MR actuators should be provided through transient models. Such a comprehensive model would attempt to copy not only the dynamics of the fluid flow through the valve and the flux dynamics but account for the physics outside the control valve as well. Suitable lumped parameter models usually utilize a network of elements representing physical domains of interest (hydraulic, thermal, electric, and magnetic) and connections (interfaces) between them [2]. For example, the electrical circuit of the actuator can be represented in the form of a resistor-nonlinear inductor network model [5]. Their main disadvantage is the necessity of using various simplifying assumptions, e.g., uniform yield stress/flux, fully developed flow, etc. On the contrary, continuum multiphysics (magnetics and flow dynamics) models utilize fewer assumptions and can be exercised on realistic geometries, however, at a significant computational expense [29].As such, in the paper, we propose a hybrid multiphysics model of the magnetorheological damper which separates the magnetic hysteresis of the magnetic circuit of the actuator from that of the mechanical hardware. Briefly, the electromagnetic domain is modelled using the vector hysteresis FE model (present in Ansys Maxwell) based on the extension of well-known Maxwell equations [30], and the hydraulic section is described through dimensionless biplastic Bingham approach [31].The paper is organized as follows. First, we present an MR damper geometry and key material properties. Then, in the following section, we reveal key details of the FE model of the actuator such as magnetic hysteresis and the coupled lumped parameter hydromechanical model. Next, we show measurements of magnetic hysteresis loops and a comparison of the measurements against the FE electromagnet model. Finally, we show results of a parametric study (also involving the hybrid model) in an attempt to examine the hysteresis influence on the output of the MR actuator and then draw conclusions.


## 2. Magnetorheological Damper

In the study, an MR flow-mode damper configuration having a single coil assembly in the electromagnet and one annular flow path in the control valve is of research interest. The damper is presented in Figure1. The hydraulic tube houses (1) the piston (2), the piston rod (3), the floating piston (4) and the rod guide assembly (5). The piston separates the MR fluid volume into rebound chamber volume and the compression chamber volume. The floating piston separates the fluid from the gas chamber. The MR valve located in the piston control controls the fluid flow between the rebound and compression chamber and vice versa. The MR valve is a conventional control valve by design. It consists of the piston core (6), the sleeve (7), the nonmagnetic flanges or plates (8), the coil assembly (9), and the connecting wires (10) for connecting to an external power supply. It is the most common single-tube MR damper configuration.Figure 1
Magnetorheological damper. (a) MR damper. (b) Piston (control valve).
(a)
(b)The MR valve’s magnetic circuit (6, 7) is manufactured out of annealed low-carbon steel 11SMn30 (see the components in blue in Figure1). The bronze (yellow) flanges (8) define the mutual position of the piston core and sleeve. The distance between the outer diameter of the annulus and the inner diameter of the sleeve defines the annular gap height. The coil assembly (9) incorporates N = 120 turns of copper (purple) wire (0.5 mm diameter).The connecting wires are routed through the thru-hole in the piston rod (3) made of steel 42CrMo4 (AISI 4140). The floating piston, the rod guide, and the remaining components are manufactured out of steel S235JR (green). The MR damper is filled with the fluid MRF132-DG by Lord Corp; see the magnetisation curve in Figure2(c). The damper key dimensions, rheological properties of MR fluid, and gas chamber details are shown in Table 1. The magnetisation curves of the annealed low-carbon steel 11SMn30 were determined using the measurement system Remagraph C-500. The obtained virgin & hysteresis data can be seen in Figures 2(a) and 2(b).Figure 2
Virgin and hysteresis magnetization curves. (a) Hysteresis curve of 11SMn30. (b) Static curve of 11SMn30. (c) Static curve of MRF 132-DG [32].
(a)
(b)
(c)Table 1
MR damper dimensions and material properties.
Name
Value
Symbol
Unit
Geometry and weight
Piston rod outer diameter
12
d
mm
Piston outer diameter
36
D
p
mm
Annular gap height
0.65
h
mm
Active zone length
16
L
a
mm
Core length
37
L
c
mm
Piston stroke
150
L
s
mm
Internal diameter of the gap
28
D
c
mm
Floating piston weight
0.07
m
f
kg
MR fluid properties (MRF-132DG [32])
MR fluid viscosity at 40°C
0.114
μ
Pa·s
MR fluid isothermal bulk modulus
1500
β
MPa
MR fluid density
3090
ρ
kg·m−3
Air content in the MR fluid
0.01
α
—
Gas chamber
Gas volume (at midstroke)
46000
V
g
0
mm3
Gas temperature
40
T
°C
Initial gas pressure
30
P
g
0
bar
Initial floating piston position
20
x
g
mm
Others
Initial rebound chamber (upper) MR fluid volume
63333
V
r0
mm3
Initial compression (lower) chamber MR fluid volume
71250
V
c0
mm3
Coil turns
120
N
—
Coil resistance
1.0
R
c
Ω
Nondimensional viscosity ratio (est.)
0.1
γ
—
Yield stress number (est.)
0.5
δ
—The coercitivity and remanence of the 11SMn30 material sample were determined from the measured hysteresis:Hc = 209 A/m and remanence Br = 1.09 T. The material’s bulk conductivity was set at 5.8 MS/m. Based on similar measurements of the rod material (42CrMo4), we set its coercivity to Hc = 1250 A/m and the bulk conductivity to 4.5 MS/m.


## 3. Modeling

In the section, we present modeling details. Specifically, we highlight the transient magnetic FE model of the MR valve including hysteresis followed by a description of a monotube damper lumped parameter model. The lumped parameter model is coupled with the transient FE model through the yield stress-flux density interface. We consider the integrated model as illustrated in Figure3. In the presented layout, the electromagnetic circuit (described in Section 3.1) is driven by the voltage u supplied by the current driver. The resulting output flux density Bg is then converted into the materials’ (fluid) field-induced yield stress τ0 (extracted from the material’s datasheet or rheological measurements). Given the input velocity vr or displacement and the yield stress, we then calculate the output force according to the equations in Section 3.2.Figure 3
Block diagram of the proposed model.


### 3.1. Transient Magnetic Model with Magnetic Hysteresis

To model the electromagnetic circuit of the MR valve, we applied the vector hysteresis modeling feature available in Ansys Maxwell R19. For isotropic material and 2D/3D problems, the vector play model was recognized to be more computationally efficient than a vector Preisach model [30, 33]. In general, the play model assumes a decomposition of the applied field H into the reversible component Hre and the irreversible one Hir. Then, the resulting flux density(1)B=Bre+Bir=μ0MHre+μ0Hir,where Bre is the reversible component of flux density and Bir is the irreversible component. The magnetization M varies with the reversible field component along an anhysteretic curve. The process can be visualized as in Figure 4. The parameters for the model can be identified from the major hysteresis loop. The major hysteresis loop incorporates two branches, the ascending branch and the descending branch, and they can be calculated from each other, the Maxwell model utilizes only one branch of particular B-H loops. The algorithm for constructing the major and symmetric minor hysteresis loops is given in [34].Figure 4
Magnetic field decomposed into reversible/irreversible components [35].To develop the FE model, we assumed the valve to be axially symmetrical around the centerline in a cylindrical coordinate system. The geometry of the MR damper piston (valve) was simplified for the transient simulations (Figure5(a)). The discretized model can be observed in Figure 5(b). As shown, the geometry was discretized using triangular elements. The element length-based refinement with the maximum length of 0.5 mm was applied to the coil core and the sleeve. Next, the value of 0.7 mm was applied to the piston rod and the hydraulic tube, and the criterion of 0.2 mm was applied to the MR fluid region in the active zone.Figure 5
Simplified geometry of the piston unit and mesh. (a) Geometry. (b) Mesh.
(a)
(b)To accomplish transient field simulations, the FE model was coupled to the external circuit revealed in Figure6(a). The lumped circuit documents an ideal current source (input) in series with a resistor (coil winding) and the FE valve object (as represented by the nonlinear inductor). The model is subjected to prescribed current waveforms, and the resulting flux density in the annulus Bg is extracted from the simulation results. The flux density Bg presented in the next sections was calculated by averaging flux density of the middle of the annulus. The information is required for coupling the FE model with the damper hydraulics.Figure 6
External circuit model and magnetisation curve: 11SMn30 (curve fit vs data). (a) External circuit. (b)B-H plot.
(a)
(b)The figure also shows the curve fit to the experimental data for the steel 11SMn30 (Figure6(b)). The agreement is satisfactory except for low magnetic field strength and on the initial magnetization curve only. The core and sleeve components were assigned the B-H properties of the 11SMn30 alloy and the MRF component that of MRF132-DG available from Lord Corp. Also, the rod component was assigned the material properties of the 42CrMo4 steel alloy as mentioned above. Finally, all data in subsequent simulations were obtained using the fixed step-size solver with the following settings: constant time step was 0.05 ms, nonlinear residual was 1e − 7, time integration method was Backward Euler.


### 3.2. Hydromechanical Model

To illustrate or reveal the effect of fluctuating (transient) magnetic field on the output of the actuator, a capable damper model is required. Modeling the behaviour of MR dampers has been clearly the subject of intensive research, to name only [36–38]. However, we chose to proceed further with the model of Goldasz and Sapinski in [2]. The approach is flexible, incorporates most key physical phenomena occurring in the MR valve and outside of it, and was successfully verified against several MR piston valve configurations (monotube damper, valve: single coil, single annular flow path, magnetic flux bypass feature). Therefore, in the sections that follow, we describe details of the lumped parameter model of the damper and the coupling method with the FE transient model.


#### 3.2.1. MR Damper Model: Theoretical Background

The schematic geometry of the damper is revealed in Figure7. In the presented illustration, the piston separates the upper (rebound) fluid chamber from the fluid below it (compression chamber). The pressure in the upper chamber is Pr, and its (initial) volume is Vr (Vr0). Accordingly, the pressure in the compression chamber is Pc, and its (initial) volume is Vc (Vc0). The gas pressure is Pg (Pg0), and the (initial) gas volume below the floating piston is referred to as Vg (Vg0). At static conditions, the pressure in each chamber is equal to Pg0. The cross-sectional area of the piston is Ap and that of the rod Ar. As the piston rod moves, it displaces the floating piston (separating the lower fluid chamber and the pressurised gas). The floating gas cup mass is mg, and its displacement is xg. The friction forces against the rod guide and the floating piston are Ffg and Ffr, respectively. We assume one annular flow path in the MR valve; dimensions: h is the gap height; w is the circumferential width (at perimeter); Ag = wh is the flow channel area. The flow rate through the annulus is referred to as Qa. Finally, we refer to the displacement of the piston rod as xr and to that of the cylinder tube as xt (not shown).Figure 7
Damper model schematic layout.The fluid’s behaviour is quantified with the viscosityμ, the density ρ, the compressibility β, and the field-induced yield stress τ0. The non-Newtonian rheology of the MR fluid is described using the biplastic Bingham model [31].We assume that the damper model would account for the following phenomena: MR effect (using the biplastic Bingham model mentioned above), compressibility of fluid, dynamics of the fluid element (“slug”) motion when forced through the annulus, entrance and exit losses in the annulus, floating piston mass inertia, and seal friction. Elasticity of the cylinder tube, various effects due to heating, and the dependency of seal friction on the damper internal pressure are not accounted for.First, the gas pressure in the volume below the floating piston can be modeled by assuming the adiabatic process (n = 1.4). The gas pressure Pg is then dependent on the position of the floating piston xg in the following manner:(2)Pg=Pg0Vg0Vg0+Apxgn.The pressure variation in the chambers below/above the piston is modeled assuming isothermal processes and the conservation of mass approach [39]. Next, the dynamics of the mass of fluid is considered by examining the motion of the fluid mass in the annular channel. The resulting system of ordinary differential equation in the state-space form which describes the mutual relationships between the rebound pressure chamber Pr, the compression chamber pressure Pc, the floating piston motion velocity vg, and the volumetric flow rate through the annulus Qa is shown below(3)P˙r=βAp−Arvp−QaVr,(4)P˙c=βvg−vpAp+QaVc,(5)v˙g=1mgApPc−Pg−Ffgsignvg,(6)Q˙a=AgρLPr−Pc−Δpa.As already mentioned, the behaviour of the energized MR fluid is described by incorporating the field-dependent losses into the pressure drop Δpa which is described in detail in Section 3.2.3; the reader should refer to [2, 39] for a more detailed derivation of the equations and the experimental verification method. The pressure term also incorporates the local flow losses Δpe. The local flow losses as the fluid enter/exits the annulus are accounted for using the semiempirical equation [40]:(7)Δpe=KρQa2A2Qa,(8)KSE=Kcor1−AgAp2,(9)KSC=Kcor1−AgAp0.75,where K/KSE is the pressure loss coefficient for the sudden enlargement (exit from the gap), K/KSC is the pressure loss coefficient for the sudden contraction (entrance to the annular gap), and Kcor is the correction factor. Finally, considering the forces acting on the piston yields the following relationship:(10)Fd=Ap−ArPr−PcAp+Ffr+Ffg.The friction force in MR damper is assumed to be the sum of Stribeck, Coulomb, and viscous components [40]. As already mentioned, the effects of viscosity change with temperature (heating) are not included.We solve the system of equations (2)–(10) using the multidomain modeling package Simscape which extends Simulink with tools for object oriented modeling and simulating multiphysics systems [40]. Our model as shown in Figure 8 consists of mechanical, hydraulic, and physical signals domains. Using that environment, the MR damper model was developed with isothermal hydraulic double-acting cylinder components, adiabatic gas blocks, and friction components. The MR fluid behaviour as copied specifically by equation (6) was defined by means of a custom component based on the biplastic Bingham model approach [31] in series with the local loss model.Figure 8
High-level Simscape model layout.To allow simulations of the transient performance of the damper, the model was coupled to the FE model in Ansys Maxwell through the magnetic flux density vs yield stress relationship,τ0 = τ0(Bg). The interface assumes zero delay between the electromagnetic response of the circuit and the MRF response. MRF measurements indicate the response time of the fluid to be below 0.6 ms; therefore, that particular contribution is omitted in the developed equation set.


#### 3.2.2. Magnetorheological Valve Model

In this section, we describe the biplastic Bingham computing scheme for determining the pressure drop across the magnetorheological valve. Specifically, the relationship between the flow rate through the annulusQa and the pressure drop Δpa is needed. The biplastic scheme is preferred rather than the conventional Bingham approach as it is more flexible and allows for a more effective modeling of low velocity features in the annulus, e.g., magnetic bypass [2]. Using the dimensionless representation of the scheme in terms of the pressure number G and the plasticity S, we express the relationship between the flow rate Qa and the pressure drop across the annulus Δpa as(11)Δpa=2τ2LahGS=2τ0Lah1−γ1−δGS,where(12)G=−hΔpa2Laτ2,S=12μQawh2τ2.The two additional parameters,γ and δ, are referred to as the artificial viscosity ratio and the (nondimensional) bypass number (yield stress ratio), respectively. As the biplastic model was well studied in prior research papers, the reader should refer there for in-depth details and the parameter estimation method. Briefly, δ controls the intercept force at the zero piston velocity, and γ influences the curve’s slope below the knee-point of the force-velocity characteristics [2]. The two parameters of the biplastic model are related to the valve’s geometry rather than material properties. The estimation procedure was highlighted, e.g., in [36]. For example, based on prior knowledge, the value of δ (0.5) was selected for a valve with no leakage flow path in the annulus. Using the model, we classify the valve’s behaviour into two flow regimes: preyield (G ≤ 1) and postyield (G > 1). The transition point coordinates at which the behaviour of the pseudomaterial changes from the preyield regime to the postyield regime are equal to G = 1, S0 = γ(2 − 3δ + δ3). Briefly, when in the postyield regime, the relationship between G and S can be expressed as(13)G=1631−γ1−δ+S2cos13atan2y,x+1,(14)y=12−81b2+12ba3,x=−108b+8a3,a=321−γ1−δ+12S,b=121−γ1−δ3,whereas in the preyield regime, the relationship between the pressure drop and the flow rate through the annulus is governed by the following formula:(15)G=δ16Sγδ+32cos13atan2y1,x1+1,where(16)y1=6327Sγδ+9Sγδ2+Sγδ3,x1=−27+27Sγδ+9Sγδ2+Sγδ3.The two model parameters (γ, δ) can be identified from real damper experimental data or CFD (computational flow dynamics) simulations. Finally, equation (11) can be modified to include the contribution of the nonenergized region above the coil of the length Lc − La through(17)Δpa=2τ2LahGS+12μLc−LaQawh3.


## 3.1. Transient Magnetic Model with Magnetic Hysteresis

To model the electromagnetic circuit of the MR valve, we applied the vector hysteresis modeling feature available in Ansys Maxwell R19. For isotropic material and 2D/3D problems, the vector play model was recognized to be more computationally efficient than a vector Preisach model [30, 33]. In general, the play model assumes a decomposition of the applied field H into the reversible component Hre and the irreversible one Hir. Then, the resulting flux density(1)B=Bre+Bir=μ0MHre+μ0Hir,where Bre is the reversible component of flux density and Bir is the irreversible component. The magnetization M varies with the reversible field component along an anhysteretic curve. The process can be visualized as in Figure 4. The parameters for the model can be identified from the major hysteresis loop. The major hysteresis loop incorporates two branches, the ascending branch and the descending branch, and they can be calculated from each other, the Maxwell model utilizes only one branch of particular B-H loops. The algorithm for constructing the major and symmetric minor hysteresis loops is given in [34].Figure 4
Magnetic field decomposed into reversible/irreversible components [35].To develop the FE model, we assumed the valve to be axially symmetrical around the centerline in a cylindrical coordinate system. The geometry of the MR damper piston (valve) was simplified for the transient simulations (Figure5(a)). The discretized model can be observed in Figure 5(b). As shown, the geometry was discretized using triangular elements. The element length-based refinement with the maximum length of 0.5 mm was applied to the coil core and the sleeve. Next, the value of 0.7 mm was applied to the piston rod and the hydraulic tube, and the criterion of 0.2 mm was applied to the MR fluid region in the active zone.Figure 5
Simplified geometry of the piston unit and mesh. (a) Geometry. (b) Mesh.
(a)
(b)To accomplish transient field simulations, the FE model was coupled to the external circuit revealed in Figure6(a). The lumped circuit documents an ideal current source (input) in series with a resistor (coil winding) and the FE valve object (as represented by the nonlinear inductor). The model is subjected to prescribed current waveforms, and the resulting flux density in the annulus Bg is extracted from the simulation results. The flux density Bg presented in the next sections was calculated by averaging flux density of the middle of the annulus. The information is required for coupling the FE model with the damper hydraulics.Figure 6
External circuit model and magnetisation curve: 11SMn30 (curve fit vs data). (a) External circuit. (b)B-H plot.
(a)
(b)The figure also shows the curve fit to the experimental data for the steel 11SMn30 (Figure6(b)). The agreement is satisfactory except for low magnetic field strength and on the initial magnetization curve only. The core and sleeve components were assigned the B-H properties of the 11SMn30 alloy and the MRF component that of MRF132-DG available from Lord Corp. Also, the rod component was assigned the material properties of the 42CrMo4 steel alloy as mentioned above. Finally, all data in subsequent simulations were obtained using the fixed step-size solver with the following settings: constant time step was 0.05 ms, nonlinear residual was 1e − 7, time integration method was Backward Euler.


## 3.2. Hydromechanical Model

To illustrate or reveal the effect of fluctuating (transient) magnetic field on the output of the actuator, a capable damper model is required. Modeling the behaviour of MR dampers has been clearly the subject of intensive research, to name only [36–38]. However, we chose to proceed further with the model of Goldasz and Sapinski in [2]. The approach is flexible, incorporates most key physical phenomena occurring in the MR valve and outside of it, and was successfully verified against several MR piston valve configurations (monotube damper, valve: single coil, single annular flow path, magnetic flux bypass feature). Therefore, in the sections that follow, we describe details of the lumped parameter model of the damper and the coupling method with the FE transient model.


### 3.2.1. MR Damper Model: Theoretical Background

The schematic geometry of the damper is revealed in Figure7. In the presented illustration, the piston separates the upper (rebound) fluid chamber from the fluid below it (compression chamber). The pressure in the upper chamber is Pr, and its (initial) volume is Vr (Vr0). Accordingly, the pressure in the compression chamber is Pc, and its (initial) volume is Vc (Vc0). The gas pressure is Pg (Pg0), and the (initial) gas volume below the floating piston is referred to as Vg (Vg0). At static conditions, the pressure in each chamber is equal to Pg0. The cross-sectional area of the piston is Ap and that of the rod Ar. As the piston rod moves, it displaces the floating piston (separating the lower fluid chamber and the pressurised gas). The floating gas cup mass is mg, and its displacement is xg. The friction forces against the rod guide and the floating piston are Ffg and Ffr, respectively. We assume one annular flow path in the MR valve; dimensions: h is the gap height; w is the circumferential width (at perimeter); Ag = wh is the flow channel area. The flow rate through the annulus is referred to as Qa. Finally, we refer to the displacement of the piston rod as xr and to that of the cylinder tube as xt (not shown).Figure 7
Damper model schematic layout.The fluid’s behaviour is quantified with the viscosityμ, the density ρ, the compressibility β, and the field-induced yield stress τ0. The non-Newtonian rheology of the MR fluid is described using the biplastic Bingham model [31].We assume that the damper model would account for the following phenomena: MR effect (using the biplastic Bingham model mentioned above), compressibility of fluid, dynamics of the fluid element (“slug”) motion when forced through the annulus, entrance and exit losses in the annulus, floating piston mass inertia, and seal friction. Elasticity of the cylinder tube, various effects due to heating, and the dependency of seal friction on the damper internal pressure are not accounted for.First, the gas pressure in the volume below the floating piston can be modeled by assuming the adiabatic process (n = 1.4). The gas pressure Pg is then dependent on the position of the floating piston xg in the following manner:(2)Pg=Pg0Vg0Vg0+Apxgn.The pressure variation in the chambers below/above the piston is modeled assuming isothermal processes and the conservation of mass approach [39]. Next, the dynamics of the mass of fluid is considered by examining the motion of the fluid mass in the annular channel. The resulting system of ordinary differential equation in the state-space form which describes the mutual relationships between the rebound pressure chamber Pr, the compression chamber pressure Pc, the floating piston motion velocity vg, and the volumetric flow rate through the annulus Qa is shown below(3)P˙r=βAp−Arvp−QaVr,(4)P˙c=βvg−vpAp+QaVc,(5)v˙g=1mgApPc−Pg−Ffgsignvg,(6)Q˙a=AgρLPr−Pc−Δpa.As already mentioned, the behaviour of the energized MR fluid is described by incorporating the field-dependent losses into the pressure drop Δpa which is described in detail in Section 3.2.3; the reader should refer to [2, 39] for a more detailed derivation of the equations and the experimental verification method. The pressure term also incorporates the local flow losses Δpe. The local flow losses as the fluid enter/exits the annulus are accounted for using the semiempirical equation [40]:(7)Δpe=KρQa2A2Qa,(8)KSE=Kcor1−AgAp2,(9)KSC=Kcor1−AgAp0.75,where K/KSE is the pressure loss coefficient for the sudden enlargement (exit from the gap), K/KSC is the pressure loss coefficient for the sudden contraction (entrance to the annular gap), and Kcor is the correction factor. Finally, considering the forces acting on the piston yields the following relationship:(10)Fd=Ap−ArPr−PcAp+Ffr+Ffg.The friction force in MR damper is assumed to be the sum of Stribeck, Coulomb, and viscous components [40]. As already mentioned, the effects of viscosity change with temperature (heating) are not included.We solve the system of equations (2)–(10) using the multidomain modeling package Simscape which extends Simulink with tools for object oriented modeling and simulating multiphysics systems [40]. Our model as shown in Figure 8 consists of mechanical, hydraulic, and physical signals domains. Using that environment, the MR damper model was developed with isothermal hydraulic double-acting cylinder components, adiabatic gas blocks, and friction components. The MR fluid behaviour as copied specifically by equation (6) was defined by means of a custom component based on the biplastic Bingham model approach [31] in series with the local loss model.Figure 8
High-level Simscape model layout.To allow simulations of the transient performance of the damper, the model was coupled to the FE model in Ansys Maxwell through the magnetic flux density vs yield stress relationship,τ0 = τ0(Bg). The interface assumes zero delay between the electromagnetic response of the circuit and the MRF response. MRF measurements indicate the response time of the fluid to be below 0.6 ms; therefore, that particular contribution is omitted in the developed equation set.


### 3.2.2. Magnetorheological Valve Model

In this section, we describe the biplastic Bingham computing scheme for determining the pressure drop across the magnetorheological valve. Specifically, the relationship between the flow rate through the annulusQa and the pressure drop Δpa is needed. The biplastic scheme is preferred rather than the conventional Bingham approach as it is more flexible and allows for a more effective modeling of low velocity features in the annulus, e.g., magnetic bypass [2]. Using the dimensionless representation of the scheme in terms of the pressure number G and the plasticity S, we express the relationship between the flow rate Qa and the pressure drop across the annulus Δpa as(11)Δpa=2τ2LahGS=2τ0Lah1−γ1−δGS,where(12)G=−hΔpa2Laτ2,S=12μQawh2τ2.The two additional parameters,γ and δ, are referred to as the artificial viscosity ratio and the (nondimensional) bypass number (yield stress ratio), respectively. As the biplastic model was well studied in prior research papers, the reader should refer there for in-depth details and the parameter estimation method. Briefly, δ controls the intercept force at the zero piston velocity, and γ influences the curve’s slope below the knee-point of the force-velocity characteristics [2]. The two parameters of the biplastic model are related to the valve’s geometry rather than material properties. The estimation procedure was highlighted, e.g., in [36]. For example, based on prior knowledge, the value of δ (0.5) was selected for a valve with no leakage flow path in the annulus. Using the model, we classify the valve’s behaviour into two flow regimes: preyield (G ≤ 1) and postyield (G > 1). The transition point coordinates at which the behaviour of the pseudomaterial changes from the preyield regime to the postyield regime are equal to G = 1, S0 = γ(2 − 3δ + δ3). Briefly, when in the postyield regime, the relationship between G and S can be expressed as(13)G=1631−γ1−δ+S2cos13atan2y,x+1,(14)y=12−81b2+12ba3,x=−108b+8a3,a=321−γ1−δ+12S,b=121−γ1−δ3,whereas in the preyield regime, the relationship between the pressure drop and the flow rate through the annulus is governed by the following formula:(15)G=δ16Sγδ+32cos13atan2y1,x1+1,where(16)y1=6327Sγδ+9Sγδ2+Sγδ3,x1=−27+27Sγδ+9Sγδ2+Sγδ3.The two model parameters (γ, δ) can be identified from real damper experimental data or CFD (computational flow dynamics) simulations. Finally, equation (11) can be modified to include the contribution of the nonenergized region above the coil of the length Lc − La through(17)Δpa=2τ2LahGS+12μLc−LaQawh3.


## 3.2.1. MR Damper Model: Theoretical Background

The schematic geometry of the damper is revealed in Figure7. In the presented illustration, the piston separates the upper (rebound) fluid chamber from the fluid below it (compression chamber). The pressure in the upper chamber is Pr, and its (initial) volume is Vr (Vr0). Accordingly, the pressure in the compression chamber is Pc, and its (initial) volume is Vc (Vc0). The gas pressure is Pg (Pg0), and the (initial) gas volume below the floating piston is referred to as Vg (Vg0). At static conditions, the pressure in each chamber is equal to Pg0. The cross-sectional area of the piston is Ap and that of the rod Ar. As the piston rod moves, it displaces the floating piston (separating the lower fluid chamber and the pressurised gas). The floating gas cup mass is mg, and its displacement is xg. The friction forces against the rod guide and the floating piston are Ffg and Ffr, respectively. We assume one annular flow path in the MR valve; dimensions: h is the gap height; w is the circumferential width (at perimeter); Ag = wh is the flow channel area. The flow rate through the annulus is referred to as Qa. Finally, we refer to the displacement of the piston rod as xr and to that of the cylinder tube as xt (not shown).Figure 7
Damper model schematic layout.The fluid’s behaviour is quantified with the viscosityμ, the density ρ, the compressibility β, and the field-induced yield stress τ0. The non-Newtonian rheology of the MR fluid is described using the biplastic Bingham model [31].We assume that the damper model would account for the following phenomena: MR effect (using the biplastic Bingham model mentioned above), compressibility of fluid, dynamics of the fluid element (“slug”) motion when forced through the annulus, entrance and exit losses in the annulus, floating piston mass inertia, and seal friction. Elasticity of the cylinder tube, various effects due to heating, and the dependency of seal friction on the damper internal pressure are not accounted for.First, the gas pressure in the volume below the floating piston can be modeled by assuming the adiabatic process (n = 1.4). The gas pressure Pg is then dependent on the position of the floating piston xg in the following manner:(2)Pg=Pg0Vg0Vg0+Apxgn.The pressure variation in the chambers below/above the piston is modeled assuming isothermal processes and the conservation of mass approach [39]. Next, the dynamics of the mass of fluid is considered by examining the motion of the fluid mass in the annular channel. The resulting system of ordinary differential equation in the state-space form which describes the mutual relationships between the rebound pressure chamber Pr, the compression chamber pressure Pc, the floating piston motion velocity vg, and the volumetric flow rate through the annulus Qa is shown below(3)P˙r=βAp−Arvp−QaVr,(4)P˙c=βvg−vpAp+QaVc,(5)v˙g=1mgApPc−Pg−Ffgsignvg,(6)Q˙a=AgρLPr−Pc−Δpa.As already mentioned, the behaviour of the energized MR fluid is described by incorporating the field-dependent losses into the pressure drop Δpa which is described in detail in Section 3.2.3; the reader should refer to [2, 39] for a more detailed derivation of the equations and the experimental verification method. The pressure term also incorporates the local flow losses Δpe. The local flow losses as the fluid enter/exits the annulus are accounted for using the semiempirical equation [40]:(7)Δpe=KρQa2A2Qa,(8)KSE=Kcor1−AgAp2,(9)KSC=Kcor1−AgAp0.75,where K/KSE is the pressure loss coefficient for the sudden enlargement (exit from the gap), K/KSC is the pressure loss coefficient for the sudden contraction (entrance to the annular gap), and Kcor is the correction factor. Finally, considering the forces acting on the piston yields the following relationship:(10)Fd=Ap−ArPr−PcAp+Ffr+Ffg.The friction force in MR damper is assumed to be the sum of Stribeck, Coulomb, and viscous components [40]. As already mentioned, the effects of viscosity change with temperature (heating) are not included.We solve the system of equations (2)–(10) using the multidomain modeling package Simscape which extends Simulink with tools for object oriented modeling and simulating multiphysics systems [40]. Our model as shown in Figure 8 consists of mechanical, hydraulic, and physical signals domains. Using that environment, the MR damper model was developed with isothermal hydraulic double-acting cylinder components, adiabatic gas blocks, and friction components. The MR fluid behaviour as copied specifically by equation (6) was defined by means of a custom component based on the biplastic Bingham model approach [31] in series with the local loss model.Figure 8
High-level Simscape model layout.To allow simulations of the transient performance of the damper, the model was coupled to the FE model in Ansys Maxwell through the magnetic flux density vs yield stress relationship,τ0 = τ0(Bg). The interface assumes zero delay between the electromagnetic response of the circuit and the MRF response. MRF measurements indicate the response time of the fluid to be below 0.6 ms; therefore, that particular contribution is omitted in the developed equation set.


## 3.2.2. Magnetorheological Valve Model

In this section, we describe the biplastic Bingham computing scheme for determining the pressure drop across the magnetorheological valve. Specifically, the relationship between the flow rate through the annulusQa and the pressure drop Δpa is needed. The biplastic scheme is preferred rather than the conventional Bingham approach as it is more flexible and allows for a more effective modeling of low velocity features in the annulus, e.g., magnetic bypass [2]. Using the dimensionless representation of the scheme in terms of the pressure number G and the plasticity S, we express the relationship between the flow rate Qa and the pressure drop across the annulus Δpa as(11)Δpa=2τ2LahGS=2τ0Lah1−γ1−δGS,where(12)G=−hΔpa2Laτ2,S=12μQawh2τ2.The two additional parameters,γ and δ, are referred to as the artificial viscosity ratio and the (nondimensional) bypass number (yield stress ratio), respectively. As the biplastic model was well studied in prior research papers, the reader should refer there for in-depth details and the parameter estimation method. Briefly, δ controls the intercept force at the zero piston velocity, and γ influences the curve’s slope below the knee-point of the force-velocity characteristics [2]. The two parameters of the biplastic model are related to the valve’s geometry rather than material properties. The estimation procedure was highlighted, e.g., in [36]. For example, based on prior knowledge, the value of δ (0.5) was selected for a valve with no leakage flow path in the annulus. Using the model, we classify the valve’s behaviour into two flow regimes: preyield (G ≤ 1) and postyield (G > 1). The transition point coordinates at which the behaviour of the pseudomaterial changes from the preyield regime to the postyield regime are equal to G = 1, S0 = γ(2 − 3δ + δ3). Briefly, when in the postyield regime, the relationship between G and S can be expressed as(13)G=1631−γ1−δ+S2cos13atan2y,x+1,(14)y=12−81b2+12ba3,x=−108b+8a3,a=321−γ1−δ+12S,b=121−γ1−δ3,whereas in the preyield regime, the relationship between the pressure drop and the flow rate through the annulus is governed by the following formula:(15)G=δ16Sγδ+32cos13atan2y1,x1+1,where(16)y1=6327Sγδ+9Sγδ2+Sγδ3,x1=−27+27Sγδ+9Sγδ2+Sγδ3.The two model parameters (γ, δ) can be identified from real damper experimental data or CFD (computational flow dynamics) simulations. Finally, equation (11) can be modified to include the contribution of the nonenergized region above the coil of the length Lc − La through(17)Δpa=2τ2LahGS+12μLc−LaQawh3.


## 4. Magnetic Flux Measurements

For the specific electromagnet geometry, we performed a series of measurements for extracting the flux density information with respect to the control (exciting) current input. The goal was to verify the FE model. Therefore, in this section, we reveal the experimental procedure for acquiring the magnetic flux relationship against the exciting current and present the obtained data.


### 4.1. Test Rig Configuration

The magnetic flux density was measured in the middle of the air gap with the ultrathin Hall transverse probe (STB1X-0201) and the magneto-meter F. W. Bell 5180 at the sampling frequency of 100 Hz. The coil current magnitude coil was simultaneously acquired by means of the Fluke i30s current clamp. The MR damper coil was excited using two laboratory power supplies: (1) Manson SDP2603 device for lower amplitude current excitations and (2) G. W. Instek PST-3202 power supply for higher current inputs. The two signals are recorded simultaneously using the front-end Dewetron USB-50-USB2-8 data acquisition module connected to the laptop (Figure9).Figure 9
Test rig configuration.The procedure was performed as follows: (1) current increase up to the maximum prescribed currentImax, which was followed by decreasing the current down to 0 A, (2) input voltage polarity change, (3) repeat Step 1, (4) repeat Step 2, and (5) repeat Step 1. Using the highlighted procedure, the magnetic flux density was measured for the maximum current levels Imax = {0.5, 1, 2, 3, 4, 5} A, respectively.The measurements of magnetic flux in the annular gap were performed without the MR fluid. Note that the relative permeability of the Hall sensor (μr = 1) placed in the thin annulus with MR fluid would distort the accuracy of the experiment as the flux flows around the probe as illustrated in Figure 10. In the simulations, we assume the presence of MR fluid would not degrade the accuracy of the model.Figure 10
Illustration of the measurement problem with Hall sensor in the annular gap in the presence of MR fluid [41].


### 4.2. Results

The obtained data are revealed in Figure11 as plots of flux density vs coil current.Figure 11
Flux density measurements: air gap, no MR fluid. (a)Imax = 0.5 A. (b) Imax = 1.0 A. (c) Imax = 2.0 A. (d) Imax = 3.0 A. (e) Imax = 4.0 A. (f) Imax = 5.0 A.
(a)
(b)
(c)
(d)
(e)
(f)Observations of the plots of flux density vs current reveal the presence of hysteresis and nonlinear behaviour with the actuator approaching the saturation at the highest current level (Imax = 5 A).


## 4.1. Test Rig Configuration

The magnetic flux density was measured in the middle of the air gap with the ultrathin Hall transverse probe (STB1X-0201) and the magneto-meter F. W. Bell 5180 at the sampling frequency of 100 Hz. The coil current magnitude coil was simultaneously acquired by means of the Fluke i30s current clamp. The MR damper coil was excited using two laboratory power supplies: (1) Manson SDP2603 device for lower amplitude current excitations and (2) G. W. Instek PST-3202 power supply for higher current inputs. The two signals are recorded simultaneously using the front-end Dewetron USB-50-USB2-8 data acquisition module connected to the laptop (Figure9).Figure 9
Test rig configuration.The procedure was performed as follows: (1) current increase up to the maximum prescribed currentImax, which was followed by decreasing the current down to 0 A, (2) input voltage polarity change, (3) repeat Step 1, (4) repeat Step 2, and (5) repeat Step 1. Using the highlighted procedure, the magnetic flux density was measured for the maximum current levels Imax = {0.5, 1, 2, 3, 4, 5} A, respectively.The measurements of magnetic flux in the annular gap were performed without the MR fluid. Note that the relative permeability of the Hall sensor (μr = 1) placed in the thin annulus with MR fluid would distort the accuracy of the experiment as the flux flows around the probe as illustrated in Figure 10. In the simulations, we assume the presence of MR fluid would not degrade the accuracy of the model.Figure 10
Illustration of the measurement problem with Hall sensor in the annular gap in the presence of MR fluid [41].


## 4.2. Results

The obtained data are revealed in Figure11 as plots of flux density vs coil current.Figure 11
Flux density measurements: air gap, no MR fluid. (a)Imax = 0.5 A. (b) Imax = 1.0 A. (c) Imax = 2.0 A. (d) Imax = 3.0 A. (e) Imax = 4.0 A. (f) Imax = 5.0 A.
(a)
(b)
(c)
(d)
(e)
(f)Observations of the plots of flux density vs current reveal the presence of hysteresis and nonlinear behaviour with the actuator approaching the saturation at the highest current level (Imax = 5 A).


## 5. Modelling Results

The series of modelling experiments was split into two stages. First, we validate the transient FE model against the experimental data, and then we study the behaviour of the hydraulic model.


### 5.1. FE Model Verification: Air Gap, No MR Fluid

Here, the FE model of the damper described in Section3.1 was verified against the obtained air gap flux density measurements. The comparison of the obtained data against the model output can be observed in Figure 12 as plots of flux density vs coil current. Due to the low current change rate, the eddy currents were neglected in the model, and only the hysteresis contribution was studied.Figure 12
Air gap flux density: measurements vs FE model. (a)Imax = 0.5 A. (b) Imax = 1.0 A. (c) Imax = 2.0 A. (d) Imax = 3.0 A. (e) Imax = 4.0 A. (f) Imax = 5.0 A.
(a)
(b)
(c)
(d)
(e)
(f)Again, observations of the plots reveal satisfactory agreement with the model anywhere except for the smallest exciting current. Overall, the plots prove the rationality of the proposed approach.


### 5.2. Hysteresis Assessment of the MR Valve

Due to reasons explained in Section4.1, direct assessment of the hysteretic behaviour of the MR valve with the fluid in the annulus was not possible with the available laboratory equipment. However, CIP- (carbonyl powder iron-) based MR fluids show virtually zero hysteresis [12]. Therefore, the presence of the fluid in the annulus only modifies the flux density-current relationship through its (nonlinear) magnetisation characteristics. Hence, it is reasonable to proceed further under the assumption that electromagnet model of the actuator was validated, and it would be accurate also in the scenarios in which the MR annulus would be filled with the fluid. Due to the magnetic circuit saturation above 2 A, we reveal the results for the exciting currents up to 2 A (Figure 13). The nonlinear contribution of the fluid is evident in the presented results.Figure 13
MR fluid contribution: calculated flux density vs. exciting current.Next, we examine the behaviour of the valve model for the two following variants:(i)
Hysteresis (core loss) ON, eddy currents ON (solid line)(ii)
Hysteresis (core loss) ON, eddy currents OFF (dashed line)The hysteresis model was applied to all 11SMn30 components (core, sleeve). As presented in Figure14, the calculated remanent flux density is relatively independent of the previous magnetic history within the examined coil current range from 0.5 A to 2 A, and the effect of eddy currents is rather insignificant in the examined case as already revealed in Figure 14.Figure 14
Time history of (gap) magnetic flux density: hysteresis ON (dashed line), hysteresis ON and eddy currents OFF (solid line), and electric current (dotted line).Furthermore, we repeated the flux density calculations for one selected electric current level (Imax = 2 A) for the following three model variants:(i)
Hysteresis switched OFF, eddy currents switched OFF (dashed line)(ii)
Hysteresis ON, eddy currents OFF (dotted line)(iii)
Hysteresis ON and eddy currents ON (solid line)The results are revealed in Figure15. It is now apparent that the hysteresis has the biggest impact on the initial output of the actuator (up to 100 mT), and the contribution of the eddy currents is negligible within the examined time scale. The eddy currents have little effect on the output flux density in the examined case (Figure 15). Clearly, the first and second scenarios are only hypothetical ones and unseen in real MR devices. However, they were included in the simulation for the purpose of isolating the contribution of a specific phenomenon.Figure 15
Time history: (gap) magnetic flux density,Imax = 2 A.Now, the obtained flux density output can be converted into the resulting yield stress and then used in all subsequent damper simulations (Figure16). As seen in the calculated data, the residual flux results in the field-induced yield stress of appr. 2.23 kPa.Figure 16
Time history: yield stress (full line) and electric current (dotted green line).


### 5.3. Damper Response to Sinusoidal Displacement Inputs

Using the set of equations (1)–(9) in Section 3.2.1 and (10)–(16) in Section 3.2.2 as well as the above yield stress map in Figure 16, we then modelled the response of the damper subjected to sinusoidal displacement inputs at fixed (constant) current levels. This was merely for demonstrating the force output of the damper corresponding to the prescribed current levels. The calculated results are highlighted as plots of force vs velocity and force vs displacement for the peak velocity Vr = 0.3 m/s in Figure 17. In the calculated example, the stroking amplitude Xr was 30 mm. Note that the data are shifted by the gas force proportional to the product of the gas pressure Pg and the piston rod area Ar (equal to 339 N). Hysteresis due to compressibility of the fluid (increasing with the current/flux level), force oscillations at piston motion reversal points, can be also observed in the generated plots, too.Figure 17
Fixed current: force vs. velocity and force vs. stroke,Vr=0.3m/s. (a) Fd vs xr. (b) Fd vs vr.
(a)
(b)Moreover, the influence of remnant flux can be seen in Figure18. In the figure, we present the results of off-state (zero current) simulations with and without demagnetization cycle. In the revealed example, the damper was subjected to a sinusoidal displacement input at the peak velocity V=0.3m/s. The blue loops reveal the damper response after demagnetization, and the red loops copy the off-state behaviour of the damper which was previously excited at Imax = 2 A. The presence of the remanent flux (appr. 54 mT) results in the contribution of 140 N.Figure 18
Off-state force: demagnetized (blue) and following 2 A current excitation, no demagnetization (red). (a)Fd vs xr. (b) Fd vs vr.
(a)
(b)


### 5.4. Damper Response to Control Current Step Inputs

The optimal performance of a S/A (semiactive) control system with an MR damper requires understanding its dynamic behaviour. In controller design studies as well as vehicle simulation, a need for modeling the MR damper system dynamics arises quite often; the time delay between the force and the control signal due to, e.g., eddy currents degrades the system performance. Therefore, in this section, we reveal the results of a series of simulations in which the influence of eddy currents, magnetic hysteresis, and hydromechanical hysteresis (due to compressibility of the fluid) is considered. In the simulations, we assumed a fast controller [6] capable of driving the coil to within 2 ms of the commanded values. The output is highlighted in Figure 19. As seen, the actuator’s off-state output is degraded by the presence of a static (residual) flux. The actuator is capable of reaching zero flux only in the hysteresis-free scenario.Figure 19
Transient FE model: time history of magnetic flux density and electric current (dashed green line).In Figure20, we further examined the force output varying the coercivity from 0 A/m to 200 A/m. It has no effect on the force output in the current (flux) rise stage; however, in the current drop case (Figure 20(b)), its contribution became evident as the actuator never reached the initial force of 240 N (prior to the current excitation).Figure 20
Force time history: damper response to current step input,Imax = 2 A (electric current course highlighted by the green dashed line). (a) Current rise. (b) Current drop.
(a)
(b)


### 5.5. Sensitivity to Gap Height: Parametric Study

Next, the following series of numerical experiments involved studying the contribution of the MR valve’s geometry on the residual force (and flux) output. The model setup incorporated the hysteresis loss model only and no eddy current mechanism. Throughout the experiment series, the coercivityHc was varied from 120 A/m to 600 A/m, and the gap size h from 0.65 mm to 1.0 mm. The generated B-H loop set (core material) is revealed in Figure 21.Figure 21
Coercivity: impact on the core material magnetisation characteristics (11SMn30).The calculated results are shown in Figure22 for the maximum electric current of 2 A. The exciting current input is identical with that in Figure 14. The coercivity influences the course of the magnetic flux density over time. The greater the coercive force applied, the slower the slope of magnetic induction observed. Moreover, the material coercivity also affects the maximum achievable magnetic flux induction. Note that the greater the coercive force, the lower the maximum magnetic flux density generated. However, the difference is minor on the order of several percent. The difference in maximum magnetic flux density can also be due to the different shape of the generated hysteresis B-H curve by Ansys Electronics. This effect will have to be further studied in follow-up research. However, one other conclusion can be drawn: the remanent flux density follows the coercivity change. The greater the coercivity, the greater the remanent flux density induced.Figure 22
Influence of gap height and coercivity on magnetic flux density. (a)h = 0.65 mm. (b) h = 0.8 mm. (c) h = 1.0 mm.
(a)
(b)
(c)The relationship of the remanent flux density vs coercivity within the range from 120 to 600 A/m is linear (Figure23). The remanent flux density decreases as the annular gap height is increased.Figure 23
Remanent flux density variation with coercivity and gap size.Furthermore, it was of research interest to compare the impact of the material coercivity on the ratioBmax/Brem, where Bmax is the gap maximum flux density, as there is a significant influence of the annular gap size on the calculated flux density ratio at low material coercivity levels (Figure 24). It is evident that the gap height does not affect the flux density ratio at higher coercivity levels. This ratio follows the turn-up ratio change of the MR damper (Figure 25). The increase of gap size and the coercivity degradation significantly increase the turn-up ratio Kf of MR damper. However, the maximum damping forces decrease at the same time as illustrated in Figure 26.Figure 24
Maximum-to-remanent flux density ratio vs material coercivity.Figure 25
Turn-up ratio variation with gap height and coercivity,Vr=0.3m/s.Figure 26
Gap height influence on damping force output,Imax = 2 A Vr=0.3m/s.


### 5.6. Parameter Sensitivity Study: Transient Response

In this section, we present the results of a parameter sensitivity study. Here, we examine the influence of the coercive forceHc, the electric conductivity σ, and the piston velocity vr on the damping force under constant velocity excitation (Vr=0.3m/s) and current step inputs.


#### 5.6.1. Influence of Material Coercitivity

To examine the contribution of the coercivityHc, we simulated the damper model response to current step inputs within the range from 120 A/m to 600 A/m. In each scenario, we assumed the material’s electric conductivity σ equal to 1 MS/m (which is typical to some silicon steel alloys or soft magnetic composites) and the presence of eddy currents. The FE results from the transient model are shown in Figure 27. The coil current step input is also shown in Figure 27 (green dashed line).Figure 27
CoercivityHc: (gap) averaged magnetic flux density time history.The corresponding time histories of the damping force following the current (flux) rise/decay are revealed in Figure28.Figure 28
Impact of the coercivityHc. Full line: damping force time history; green dashed line: current step input. (a) Force vs time: current rise (initial condition = demag.). (b) Force vs time: current drop. (c) Force vs time: current rise (no demag.).
(a)
(b)
(c)


#### 5.6.2. Influence of Electric Conductivity

It is well known that the ferromagnetic material’s electric conductivity has a rather dramatic effect on the eddy currents induced in solenoid structures. Here, we simulate the transient response of the damper subjected to 2 A current step inputs. The value of 5.8 MS/m is typical of low-carbon steel alloys, whereas 1 MS/m characterizes some soft magnetic composite materials. The resulting time history is revealed in Figure29. As shown, as the conductivity decreases, the eddy currents are reduced, and the response of the damper becomes faster.Figure 29
Core material’s electric conductivity: Full line: damping force time history; green dashed line: current step input.Hc = 200 A/m, Vr=0.3m/s.


#### 5.6.3. Influence of Piston Velocity

Setting the electrical conductivity to 1 MS/m and the coercitivity to 200 A/m, we then tested the influence of the piston velocity on the force. In the presented examples, the force output time histories were normalized for better comparison. The obtained results imply that the slower the piston velocity for control current rise, the slower the magnitude of the force change rate generated. After exceeding the piston velocity of 0.2 m/s, the actuator response in the current rise stage is independent of the piston velocity (Figure30). However, the exact velocity value will depend on the particular damper design. For comparison, the actuator response in the current drop stage is independent of the prescribed piston velocity.Figure 30
Influence of piston velocity for control current rise and drop (green dashed line) on the force output (full line). (a) Current rise. (b) Current drop.
(a)
(b)Apart from the eddy currents, the main source of slower force rise is the compressibility of the MR fluid itself. Three different values of MR fluid bulk modulus were tested to illustrate this effect (Figure31).Figure 31
Influence of MR fluid bulk modulus on the force output at velocity 0.1 m/s.The primary response time of force (63.3% of final force) was calculated from the simulated data (Figure30) (Figure 32). The primary response time for control current rise is influenced by the piston velocity. The lower the piston velocity, the lower the primary response time. However, the primary response time for control current drop is independent of piston velocity. It is noteworthy that similar trends were experimentally determined in other research studies [42, 43].Figure 32
Influence of piston velocity on the primary response time for the rise and drop control current.


## 5.1. FE Model Verification: Air Gap, No MR Fluid

Here, the FE model of the damper described in Section3.1 was verified against the obtained air gap flux density measurements. The comparison of the obtained data against the model output can be observed in Figure 12 as plots of flux density vs coil current. Due to the low current change rate, the eddy currents were neglected in the model, and only the hysteresis contribution was studied.Figure 12
Air gap flux density: measurements vs FE model. (a)Imax = 0.5 A. (b) Imax = 1.0 A. (c) Imax = 2.0 A. (d) Imax = 3.0 A. (e) Imax = 4.0 A. (f) Imax = 5.0 A.
(a)
(b)
(c)
(d)
(e)
(f)Again, observations of the plots reveal satisfactory agreement with the model anywhere except for the smallest exciting current. Overall, the plots prove the rationality of the proposed approach.


## 5.2. Hysteresis Assessment of the MR Valve

Due to reasons explained in Section4.1, direct assessment of the hysteretic behaviour of the MR valve with the fluid in the annulus was not possible with the available laboratory equipment. However, CIP- (carbonyl powder iron-) based MR fluids show virtually zero hysteresis [12]. Therefore, the presence of the fluid in the annulus only modifies the flux density-current relationship through its (nonlinear) magnetisation characteristics. Hence, it is reasonable to proceed further under the assumption that electromagnet model of the actuator was validated, and it would be accurate also in the scenarios in which the MR annulus would be filled with the fluid. Due to the magnetic circuit saturation above 2 A, we reveal the results for the exciting currents up to 2 A (Figure 13). The nonlinear contribution of the fluid is evident in the presented results.Figure 13
MR fluid contribution: calculated flux density vs. exciting current.Next, we examine the behaviour of the valve model for the two following variants:(i)
Hysteresis (core loss) ON, eddy currents ON (solid line)(ii)
Hysteresis (core loss) ON, eddy currents OFF (dashed line)The hysteresis model was applied to all 11SMn30 components (core, sleeve). As presented in Figure14, the calculated remanent flux density is relatively independent of the previous magnetic history within the examined coil current range from 0.5 A to 2 A, and the effect of eddy currents is rather insignificant in the examined case as already revealed in Figure 14.Figure 14
Time history of (gap) magnetic flux density: hysteresis ON (dashed line), hysteresis ON and eddy currents OFF (solid line), and electric current (dotted line).Furthermore, we repeated the flux density calculations for one selected electric current level (Imax = 2 A) for the following three model variants:(i)
Hysteresis switched OFF, eddy currents switched OFF (dashed line)(ii)
Hysteresis ON, eddy currents OFF (dotted line)(iii)
Hysteresis ON and eddy currents ON (solid line)The results are revealed in Figure15. It is now apparent that the hysteresis has the biggest impact on the initial output of the actuator (up to 100 mT), and the contribution of the eddy currents is negligible within the examined time scale. The eddy currents have little effect on the output flux density in the examined case (Figure 15). Clearly, the first and second scenarios are only hypothetical ones and unseen in real MR devices. However, they were included in the simulation for the purpose of isolating the contribution of a specific phenomenon.Figure 15
Time history: (gap) magnetic flux density,Imax = 2 A.Now, the obtained flux density output can be converted into the resulting yield stress and then used in all subsequent damper simulations (Figure16). As seen in the calculated data, the residual flux results in the field-induced yield stress of appr. 2.23 kPa.Figure 16
Time history: yield stress (full line) and electric current (dotted green line).


## 5.3. Damper Response to Sinusoidal Displacement Inputs

Using the set of equations (1)–(9) in Section 3.2.1 and (10)–(16) in Section 3.2.2 as well as the above yield stress map in Figure 16, we then modelled the response of the damper subjected to sinusoidal displacement inputs at fixed (constant) current levels. This was merely for demonstrating the force output of the damper corresponding to the prescribed current levels. The calculated results are highlighted as plots of force vs velocity and force vs displacement for the peak velocity Vr = 0.3 m/s in Figure 17. In the calculated example, the stroking amplitude Xr was 30 mm. Note that the data are shifted by the gas force proportional to the product of the gas pressure Pg and the piston rod area Ar (equal to 339 N). Hysteresis due to compressibility of the fluid (increasing with the current/flux level), force oscillations at piston motion reversal points, can be also observed in the generated plots, too.Figure 17
Fixed current: force vs. velocity and force vs. stroke,Vr=0.3m/s. (a) Fd vs xr. (b) Fd vs vr.
(a)
(b)Moreover, the influence of remnant flux can be seen in Figure18. In the figure, we present the results of off-state (zero current) simulations with and without demagnetization cycle. In the revealed example, the damper was subjected to a sinusoidal displacement input at the peak velocity V=0.3m/s. The blue loops reveal the damper response after demagnetization, and the red loops copy the off-state behaviour of the damper which was previously excited at Imax = 2 A. The presence of the remanent flux (appr. 54 mT) results in the contribution of 140 N.Figure 18
Off-state force: demagnetized (blue) and following 2 A current excitation, no demagnetization (red). (a)Fd vs xr. (b) Fd vs vr.
(a)
(b)


## 5.4. Damper Response to Control Current Step Inputs

The optimal performance of a S/A (semiactive) control system with an MR damper requires understanding its dynamic behaviour. In controller design studies as well as vehicle simulation, a need for modeling the MR damper system dynamics arises quite often; the time delay between the force and the control signal due to, e.g., eddy currents degrades the system performance. Therefore, in this section, we reveal the results of a series of simulations in which the influence of eddy currents, magnetic hysteresis, and hydromechanical hysteresis (due to compressibility of the fluid) is considered. In the simulations, we assumed a fast controller [6] capable of driving the coil to within 2 ms of the commanded values. The output is highlighted in Figure 19. As seen, the actuator’s off-state output is degraded by the presence of a static (residual) flux. The actuator is capable of reaching zero flux only in the hysteresis-free scenario.Figure 19
Transient FE model: time history of magnetic flux density and electric current (dashed green line).In Figure20, we further examined the force output varying the coercivity from 0 A/m to 200 A/m. It has no effect on the force output in the current (flux) rise stage; however, in the current drop case (Figure 20(b)), its contribution became evident as the actuator never reached the initial force of 240 N (prior to the current excitation).Figure 20
Force time history: damper response to current step input,Imax = 2 A (electric current course highlighted by the green dashed line). (a) Current rise. (b) Current drop.
(a)
(b)


## 5.5. Sensitivity to Gap Height: Parametric Study

Next, the following series of numerical experiments involved studying the contribution of the MR valve’s geometry on the residual force (and flux) output. The model setup incorporated the hysteresis loss model only and no eddy current mechanism. Throughout the experiment series, the coercivityHc was varied from 120 A/m to 600 A/m, and the gap size h from 0.65 mm to 1.0 mm. The generated B-H loop set (core material) is revealed in Figure 21.Figure 21
Coercivity: impact on the core material magnetisation characteristics (11SMn30).The calculated results are shown in Figure22 for the maximum electric current of 2 A. The exciting current input is identical with that in Figure 14. The coercivity influences the course of the magnetic flux density over time. The greater the coercive force applied, the slower the slope of magnetic induction observed. Moreover, the material coercivity also affects the maximum achievable magnetic flux induction. Note that the greater the coercive force, the lower the maximum magnetic flux density generated. However, the difference is minor on the order of several percent. The difference in maximum magnetic flux density can also be due to the different shape of the generated hysteresis B-H curve by Ansys Electronics. This effect will have to be further studied in follow-up research. However, one other conclusion can be drawn: the remanent flux density follows the coercivity change. The greater the coercivity, the greater the remanent flux density induced.Figure 22
Influence of gap height and coercivity on magnetic flux density. (a)h = 0.65 mm. (b) h = 0.8 mm. (c) h = 1.0 mm.
(a)
(b)
(c)The relationship of the remanent flux density vs coercivity within the range from 120 to 600 A/m is linear (Figure23). The remanent flux density decreases as the annular gap height is increased.Figure 23
Remanent flux density variation with coercivity and gap size.Furthermore, it was of research interest to compare the impact of the material coercivity on the ratioBmax/Brem, where Bmax is the gap maximum flux density, as there is a significant influence of the annular gap size on the calculated flux density ratio at low material coercivity levels (Figure 24). It is evident that the gap height does not affect the flux density ratio at higher coercivity levels. This ratio follows the turn-up ratio change of the MR damper (Figure 25). The increase of gap size and the coercivity degradation significantly increase the turn-up ratio Kf of MR damper. However, the maximum damping forces decrease at the same time as illustrated in Figure 26.Figure 24
Maximum-to-remanent flux density ratio vs material coercivity.Figure 25
Turn-up ratio variation with gap height and coercivity,Vr=0.3m/s.Figure 26
Gap height influence on damping force output,Imax = 2 A Vr=0.3m/s.


## 5.6. Parameter Sensitivity Study: Transient Response

In this section, we present the results of a parameter sensitivity study. Here, we examine the influence of the coercive forceHc, the electric conductivity σ, and the piston velocity vr on the damping force under constant velocity excitation (Vr=0.3m/s) and current step inputs.


### 5.6.1. Influence of Material Coercitivity

To examine the contribution of the coercivityHc, we simulated the damper model response to current step inputs within the range from 120 A/m to 600 A/m. In each scenario, we assumed the material’s electric conductivity σ equal to 1 MS/m (which is typical to some silicon steel alloys or soft magnetic composites) and the presence of eddy currents. The FE results from the transient model are shown in Figure 27. The coil current step input is also shown in Figure 27 (green dashed line).Figure 27
CoercivityHc: (gap) averaged magnetic flux density time history.The corresponding time histories of the damping force following the current (flux) rise/decay are revealed in Figure28.Figure 28
Impact of the coercivityHc. Full line: damping force time history; green dashed line: current step input. (a) Force vs time: current rise (initial condition = demag.). (b) Force vs time: current drop. (c) Force vs time: current rise (no demag.).
(a)
(b)
(c)


### 5.6.2. Influence of Electric Conductivity

It is well known that the ferromagnetic material’s electric conductivity has a rather dramatic effect on the eddy currents induced in solenoid structures. Here, we simulate the transient response of the damper subjected to 2 A current step inputs. The value of 5.8 MS/m is typical of low-carbon steel alloys, whereas 1 MS/m characterizes some soft magnetic composite materials. The resulting time history is revealed in Figure29. As shown, as the conductivity decreases, the eddy currents are reduced, and the response of the damper becomes faster.Figure 29
Core material’s electric conductivity: Full line: damping force time history; green dashed line: current step input.Hc = 200 A/m, Vr=0.3m/s.


### 5.6.3. Influence of Piston Velocity

Setting the electrical conductivity to 1 MS/m and the coercitivity to 200 A/m, we then tested the influence of the piston velocity on the force. In the presented examples, the force output time histories were normalized for better comparison. The obtained results imply that the slower the piston velocity for control current rise, the slower the magnitude of the force change rate generated. After exceeding the piston velocity of 0.2 m/s, the actuator response in the current rise stage is independent of the piston velocity (Figure30). However, the exact velocity value will depend on the particular damper design. For comparison, the actuator response in the current drop stage is independent of the prescribed piston velocity.Figure 30
Influence of piston velocity for control current rise and drop (green dashed line) on the force output (full line). (a) Current rise. (b) Current drop.
(a)
(b)Apart from the eddy currents, the main source of slower force rise is the compressibility of the MR fluid itself. Three different values of MR fluid bulk modulus were tested to illustrate this effect (Figure31).Figure 31
Influence of MR fluid bulk modulus on the force output at velocity 0.1 m/s.The primary response time of force (63.3% of final force) was calculated from the simulated data (Figure30) (Figure 32). The primary response time for control current rise is influenced by the piston velocity. The lower the piston velocity, the lower the primary response time. However, the primary response time for control current drop is independent of piston velocity. It is noteworthy that similar trends were experimentally determined in other research studies [42, 43].Figure 32
Influence of piston velocity on the primary response time for the rise and drop control current.


## 5.6.1. Influence of Material Coercitivity

To examine the contribution of the coercivityHc, we simulated the damper model response to current step inputs within the range from 120 A/m to 600 A/m. In each scenario, we assumed the material’s electric conductivity σ equal to 1 MS/m (which is typical to some silicon steel alloys or soft magnetic composites) and the presence of eddy currents. The FE results from the transient model are shown in Figure 27. The coil current step input is also shown in Figure 27 (green dashed line).Figure 27
CoercivityHc: (gap) averaged magnetic flux density time history.The corresponding time histories of the damping force following the current (flux) rise/decay are revealed in Figure28.Figure 28
Impact of the coercivityHc. Full line: damping force time history; green dashed line: current step input. (a) Force vs time: current rise (initial condition = demag.). (b) Force vs time: current drop. (c) Force vs time: current rise (no demag.).
(a)
(b)
(c)


## 5.6.2. Influence of Electric Conductivity

It is well known that the ferromagnetic material’s electric conductivity has a rather dramatic effect on the eddy currents induced in solenoid structures. Here, we simulate the transient response of the damper subjected to 2 A current step inputs. The value of 5.8 MS/m is typical of low-carbon steel alloys, whereas 1 MS/m characterizes some soft magnetic composite materials. The resulting time history is revealed in Figure29. As shown, as the conductivity decreases, the eddy currents are reduced, and the response of the damper becomes faster.Figure 29
Core material’s electric conductivity: Full line: damping force time history; green dashed line: current step input.Hc = 200 A/m, Vr=0.3m/s.


## 5.6.3. Influence of Piston Velocity

Setting the electrical conductivity to 1 MS/m and the coercitivity to 200 A/m, we then tested the influence of the piston velocity on the force. In the presented examples, the force output time histories were normalized for better comparison. The obtained results imply that the slower the piston velocity for control current rise, the slower the magnitude of the force change rate generated. After exceeding the piston velocity of 0.2 m/s, the actuator response in the current rise stage is independent of the piston velocity (Figure30). However, the exact velocity value will depend on the particular damper design. For comparison, the actuator response in the current drop stage is independent of the prescribed piston velocity.Figure 30
Influence of piston velocity for control current rise and drop (green dashed line) on the force output (full line). (a) Current rise. (b) Current drop.
(a)
(b)Apart from the eddy currents, the main source of slower force rise is the compressibility of the MR fluid itself. Three different values of MR fluid bulk modulus were tested to illustrate this effect (Figure31).Figure 31
Influence of MR fluid bulk modulus on the force output at velocity 0.1 m/s.The primary response time of force (63.3% of final force) was calculated from the simulated data (Figure30) (Figure 32). The primary response time for control current rise is influenced by the piston velocity. The lower the piston velocity, the lower the primary response time. However, the primary response time for control current drop is independent of piston velocity. It is noteworthy that similar trends were experimentally determined in other research studies [42, 43].Figure 32
Influence of piston velocity on the primary response time for the rise and drop control current.


## 6. Conclusions and Summary

In this paper, we present the results of a modeling study involving a multiphysics model of a flow-mode MR damper. The model allows integrating an FE electromagnet model of the device with a hydraulic lumped parameter model of the device. The modeling approach relies only on the information which can be extracted from engineering drawings (geometry), material data sheets (material properties), and therefore, it can be used in studies on the performance of real actuators or virtual prototypes.The electromagnet was verified experimentally. Based on the obtained data, we conclude that the model is capable of predicting the magnetic hysteretic behaviour of the MR valve. The results concerning the hydraulic model are simulated; however, it should be noted that the model is based on a well-established and experimentally verified theory [39]. To demonstrate the usefulness of the model, we applied it in a parametric study, in which the contribution of various geometric parameters and material properties to the output of the actuator was studied and analyzed. For instance, we can conclude the following.(i)
Residual magnetic flux is directly related to the current history and the annular gap height(ii)
Larger annular gaps induce lower remanent (residual) flux density and then less undesired force increase (Figures22 and 23)(iii)
Valves with larger annular gaps height reveal higher turn-up ratio (dynamic range), however, at the expense of maximum damping forces (Figures25 and 26)(iv)
Turn-up ratio (dynamic range) varies with the coercivity and gap height (Figure25)(v)
Demagnetizing current cycles are required to reduce/eliminate the residual flux (and the force increment due to the residual flux)(vi)
The electrical conductivity has a major influence on the dynamic behaviour of the actuator (Figure29)(vii)
The piston velocity influences the actuator’s response time. Low piston velocities degrade the response time in the current rise stage (Figures30 and 32). It is likely due to the compressibility of the fluid (Figure 31).The collected data enhance understanding the key mechanisms governing the flux/force output of MR actuators. They allow for a clear separation of various contributors to the static and dynamic behaviour of such devices. In our opinion, the proposed model can be a useful tool as it incorporates major key phenomena occurring in the actuator including magnetic hysteresis and remanence/coercivity, eddy currents, compressibility, and fluid inertia (hydraulic hysteresis), the MR effect. It allows separating the magnetic hysteresis from the hydromechanical one so that the two phenomena can be examined separately.



---

*Source: 3246915-2019-06-26.xml*