MR Damper Modeling Performance Comparison Including Hysteresis and Damper Optimization

This research paper represents the analysis and simulation of semi-active suspension using non-linear modeling of the Magneto-Rheological (MR) suspension with consideration of the hysteresis behavior for a quarter car model. The research is based on the assumption that each wheel experiences the same disturbance excitation. Hysteresis is analyzed using Bingham, Dahl’s, and Bouc-Wen models. This research focuses on simulation of passive, Bingham, Dahl, and Bouc-Wen models and analysis for the five road profiles. The desired damping force determines the optimum working conditions based on optimized critical design parameters. An integrated approach towards the numeric design optimization by computational methods has been used to find the optimum working conditions. The critical parameters of MR damper are determined, and multi-objective optimization is performed considering the pole length, piston radius, gap thickness, piston internal radius, piston velocity and coil current. Sensitivity analysis also is performed to identify the sensitive parameter in the MR damper geometry towards the damping force. Analysis shows that gap thickness is the most sensitive geometric parameter of the MR damper. Furthermore, the comparative study of the models for the highest comfort with less overshoot and settling time are executed. The Bouc-Wen model is 36.91% more accurate than passive suspension in terms of damping force requirements, has a 26.16% less overshoot, and 88.31% less settling time. The simulation of the Bouc-Wen model yields the damping force requirement to be 2150N which is 91% of the analytical results.


I. INTRODUCTION
Automotive technology has been continuously incorporating developments over the past few decades to provide end users with better riding comfort.One such subsystem undergoing rigorous changes is the automotive suspension system, and the invention of electromagnetic capabilities of certain materials suspended in viscous oils has led to the development of the Magnetorheological (MR) damper technology.MR suspension has been an over an overlooked technology in the past, when it could substitute the conventional damper system with a smart controlled damper to reduce sprung mass acceleration based on road profile conditions.MR fluid is responsible for the operational behavior and the significant performance changes in the MR damper.Magnetorheological fluids have recently been gaining popularity in automotive component applications such as engine mounts, clutches, brakes and dampers for suspension systems [1], [2].
The associate editor coordinating the review of this manuscript and approving it for publication was Feiqi Deng .
Controlling the motion of the vehicle body undergoing road profile variations has been a challenge to the engineering world owing to the complexities arising in the multiple degree of freedom vibration response of the vehicle.Although electrically controlled suspension systems have been used to improve the dynamic performance of vehicles, such systems are still limited due to discontinuous damping forces, structural complexity, and high cost [3], [4].MR dampers eliminate the requirement of a bulky reservoir in the counterpart of pneumatic shock absorbers and come with the additional advantage of simple construction.MR dampers are essentially dampers with variable effective viscosity due to yield shear stress changes induced by excitation current controlled magnetic field strength in the damper coil [5], [6].
MR damper is a non-linear device with hysteretic characteristics, thus output of the MR damper is dependent on both the previous outputs and on the instantaneous values [6].The hysteretic study of the MR damper is carried out in this research considering four approaches: passive model without hysteresis, Bingham model, Dahl's model, and Bouc-Wen model.The mathematical formulation for all of these models is well-established [7]- [9].The comparative analysis illustrates the accurate model for design by comparing the force generated by the damper for the road disturbance characteristics at any given point of time.The analytical calculations comprise the most critical phase of the design process as we define the MR damper parameters in terms of mathematical equations.This process involved the determination of spring force from values of quarter car mass, motion ratio of the vehicle, desired spring frequency and the damping ratio range [10].The critical parameters of MR damper are determined using multi-objective optimization considering the critical parameters: pole length, piston radius, gap thickness, piston internal radius, piston velocity and coil current.The objective function constitute of damping force optimization for the desired damping force range with constraint on the critical design parameters.

II. MATHEMATICAL MODELING OF QUARTER CAR SUSPENSION
A quarter car model with passive and semi-active suspension is examined for the performance analysis of MR suspension.The results are compared to ascertain the performance in terms of vertical displacement (overshoot) and the settling time of the sprung mass due to road disturbances.The governing equations explained in the following sections are used to model the quarter car suspension followed by its analysis for various road profiles.

A. PASSIVE SUSPENSION
The passive suspension system features a fixed spring and damper.The spring and damper characteristics are determined according to the goals and intended application.Governing equations for the passive suspension are determined from the free body diagram of quarter car suspension physical model [11].
Sprung mass dynamics for the passive suspension can be written from the free body diagram shown in figure 1.
Rearranging the equation ( 1), we can find the acceleration of the sprung mass.
Unsprung mass dynamics for the passive suspension can be written from the free body diagram shown in figure 1.
Rearranging the equation (3), we can find the acceleration of the unsprung mass.

B. SEMI-ACTIVE SUSPENSION
A semi-active suspension system consists of a conventional spring and a controllable shock absorber or damper.The viscous damping coefficient of the damper can be controlled in real time, which gives an advantage over the passive suspension system.The present study focuses on MR damper.
Governing equations for the semi-active suspension are determined from the free body diagram of semi-active quarter car suspension physical model [11].Sprung mass dynamics for the semi-active suspension can be written from the free body diagram shown in figure 2.
Uc is the variable and controllable control force that can be generated according to the particular fluid properties used as a damper fluid.Rearranging the equation ( 5), we can find the acceleration of the sprung mass.
m s (6) Unsprung mass dynamics for the semi-active suspension can be written from the free body diagram shown in figure 1.
Rearranging the equation ( 7), we can find the acceleration of the unsprung mass.

C. HYSTERESIS MODELS
Hysteresis is a phenomenon common to a broad spectrum of physical systems.As such, it is often present in plants for which controllers are being designed, where it introduces a nonlinear multi-valued behavior [12].The hysteresis in the MR suspension is an important phenomenon to simulate and study since the non-linear effects of hysteresis result in discontinuous relationship between the current in electromagnet and damping provided by the MR damper [13].This research takes into consideration the hysteresis existed in the MR models.Bingham, Dahl's, and Bouc-Wen models are modeled and simulated using the MATLAB/Simulink software package for the analysis of quarter car model with hysteresis loop in the simulation to identify the accurate response of the system for the various road profiles.The hysteresis parameters were obtained based on magnetic fluid properties for the given road profile.

D. BINGHAM MODEL
Bingham model is one of the initial models of the MR damper.Bingham plastic model behaves as solid until the minimum yield stress is reached.After reaching the minimum yield stress point, it then follows the linear relationship between stress and the deformation.Bingham plastic model was proposed in 1985, and it can be formulated as follows [7]: Signum (sgn) function takes care of direction of friction force Fc relative to the relative velocity of hysteresis variable z.

E. DAHL MODEL
Dahl model of MR damper [8] considers the quasi-static bonds in the origin of friction [9].The dahl model is formulated as follows: F. BOUC-WEN MODEL Bouc-Wen model takes into consideration the spring stiffness element, conventional damper, and the Bouc-Wen hysteresis loop elements [9], [15].Bouc-Wen model is formulated as follows [14]- [16]: y is the evolutionary variable, which is dependent on ϒ, β, and A.
The nature of y changes from the sinusoidal to the quasirectangular function.The resulting force generated by the MR damper is calculated using the following equation: Coefficients C s (u) and α (u) are determined by the following equations:

III. COMPARATIVE ANALYSIS AND DAMPING FORCE DETERMINATION
Suspension models are developed based on the governing equations mentioned in section II, and then the comparative analysis is performed.The models are simulated using five input road profiles: step, sine wave, white noise, random number and mixed inputs.The step response simulates the sudden disturbance in the road profile.Sine wave simulates the wavy road disturbance that is simultaneous crest and trough profile of the road.The band limited white noise and random uniform number road profile simulates the unpredicted/random excitation of road profile.Combined sine wave and random uniform number simulates the off-roading road profile.

A. STEP INPUT
A road step input of 75 mm is used as road disturbance to the wheel [8].The assumption is that all wheels experience the same road excitation.Figure 3 and figure 4 represent the response of each model in terms of sprung mass displacement, and magnitude and phase change respectively for the given road profile.other models.The step input simulates the sudden bump condition, which is a rare phenomenon for real operating conditions.Thus, consideration of settling time as well as overshoot is very important for this road profile.Bouc-Wen model has the lowest overshoot and lowest settling time for such conditions.Also, from figure 4 it is observed that Bouc-Wen model goes through the least phase changes until the vibrations are damped, which results into least hysteresis loss.

B. SINE WAVE INPUT
Sine wave input of amplitude 75 mm and frequency of 20.8 rad/sec is given as road input [8].The sine wave simulates the wavy road disturbance -simultaneous crest and trough profile of the road.Figures 5 and 6 represent the responses of each model in terms of sprung mass displacement, magnitude and phase change respectively for the sine road profile input.From figure 5 it is observed that the Bouc-Wen model has least overshoot as well as least sprung mass vertical displacement superposition, and hence results in the maximum comfort of passengers compared to other models.

C. WHITE NOISE INPUT
White noise is a random vibration signal with uniform intensity over the time with varying frequency which represents the road roughness variations for simulation.From Figures 7 and 8, it is observed that the Bouc-Wen model has least overshoot and settling time compared to the remaining models.Figure 7 shows the vehicle sprung mass displacement

D. UNIFORM RANDOM NUMBER INPUT
There is no single definite method to analyze or synthesize the vibrations generated by a vehicle travelling over an irregular terrain.The apt method thus assumes that the vibrations can be approximated by a zero mean, normally distributed random (Gaussian) signal [17].Models are simulated using the uniform random number road profile with minimum and maximum bound −75 mm and 75 mm respectively.Figure 9 represents that the Bouc-Wen is the most efficient model under such conditions and, in turn, provides the maximum ride comfort.Figure 10 shows the sprung mass displacement amplitude and phase change for the given road profile.

E. MIXED SINE WAVE AND UNIFORM RANDOM NUMBER INPUT
This road profile is a combination of the wavy and rough road profile, which resembles off-roading conditions.Sine wave input of amplitude 75 mm, frequency 20.8 rad/sec and random number road profile with minimum and maximum bound −75 and 75 respectively and 0.1 as a sampling time are coupled as the excitation from road to wheel [8]. Figure 11 shows that the Bouc-Wen model has both the least overshoot and settling time compared to the other models, thus it provides the maximum ride comfort to the passengers.Figure 12 further illustrates the sprung mass displacement in vertical direction due to the road profile roughness, and the phase changes due to change in the road profile.The analysis of the models has been performed for the various road profiles as discussed in previous section.The Bouc-Wen model shows the most efficient and optimal results with least overshoot and less settling time for the sine wave input, white noise input, uniform random number input, and the mixed sine wave and uniform random number input.The road profiles considered for the simulation and analysis reflect real time operating conditions.For the above stated reason Bouc-Wen model is considered for the MR damper design.The damping force requirement analysis is performed for the Bouc-Wen with the maximum disturbance of the road profile.Figure 13 shows the damping force requirement for the maximum disturbance of 75 mm.The peak finder analysis shows the maximum damping force requirement for the damper is 2006 N as mentioned in the table 1. Considering the factor of safety in unexpected conditions, damper is designed and geometrically optimized for the force of 2150 N.

IV. ANALYTICAL CALCULATIONS OF MR DAMPER
The damper is designed to operate within the damping force range of 790 to 2150N, which is derived from the damping force analysis shown in figure 13 and table 1.The analytical design of the MR damper involves identification of important geometric parameters, using mathematical calculations and appropriate fluid flow equations through the damper body to define the operating range of the damping force.A Quasistatic quarter car model of the damper formulation involves determination of the spring force from basics of vehicle dynamics for desirable ride frequencies of vibration of the sprung mass.This is based on various factors such as the type of vehicle for which the suspension is being designed, terrains for which the suspension is being designed i.e. highways or off-roading.The frequency of vibration for a vehicle being driven mostly on off-road terrains should be around 0.7 Hz and that for a regular car being driven on smooth roads is 1.6 Hz [16].To facilitate a smooth operation under both conditions, the frequency adapted in this study is 1.2 Hz [16].
The equations for the design of basic spring characteristics can be obtained using equations 17 and 18 [10].
For a passenger car, a feasible spring length is 38.10 cm (15 inches) and the allowable deflection for smooth operation is 11.43 cm (4.5 inches) [18].Based on the requirements, the maximum allowable force acting on the spring can be determined using equations 19 & 20.
Using equations 17 to 20 for a passenger car the spring force is 3306N.This value is in the acceptable range of 3000−4500 N [18].
Further, an automotive suspension system is required to be underdamped for the smooth transition of vibrations with minimal shock to the passenger [10].A critically and overdamped system will quickly reduce the sprung mass displacement magnitude but that will result in a jerk to the passenger.Considering an underdamped system, the damping force range is calculated for various operating conditions from the spring force using the minimum and maximum damping conditions or damping ratios.The limits on damping ratios (ζ min and ζ max ) used is between 0.25 and 0.65, where 0.25 is for an off-road terrain where minimum damping is desirable for a greater force transmissibility from ground to the sprung mass for a smoother ride [18].Thus, the desired damping forces are 826 N to 2150 N using equations 21 and 22 [10].
Damping Force min = ζ min * F max (21) Damping Force max = ζ max * F max (22) Thus, the damper needs to generate a minimum force of 826 N, when there is no magnetic excitation of the MR fluid and the damping force is a function of the viscosity of the MR fluid and the piston velocity.Conversely, the damper needs to generate a maximum force of 2150N when the MR fluid is excited using electromagnetic induction.
The primary factors affecting the damping force are damper geometry, inductive current and the MR fluid characteristics [1].The Pole Length (L), Piston Radius (R), and gap thickness (g) are the dominant geometrical parameters whereas the current through the coil and piston velocity are the dominant non-geometric factors.Figure 14 shows a typical damper cross sectional view and the magnetic links of the circuit in the MR damper.From figure 14 and basic concepts of fluid mechanics, each parameter and its effect on the damping force can be evaluated as described further.
The viscous damping force is directly proportional to the piston radius which is a function of fluid flow resistance.The current dependent damping force is also directly proportional to the piston radius as the increase in piston radius increases the coil width which eventually means more current passing through the coil.
Gap thickness also affects both the viscous and current dependent damping force.It is the only parameter, which varies inversely with respect to the damping force, i.e. a wider gap lessens the damping force, as the magnetic flux density is less in a wider gap because fluid provides a high reluctance in the magnetic circuit.Wide gap also allows for easier flow of the fluid from the top of the piston to its bottom which reduces the viscous damping force due to the decrease in resistance to the flow.
Pole length is a parameter which only affects the current dependent damping force.Pole length is normal to the axial flow direction, and shear resistance due to geometry is negligible.However, the greater the pole length, the greater the area available for magnetic flux lines to pass through the circuit.
Coil current is the most important variable for the MR damper design since the current dependent damping force varies exponentially with a change in coil current.This is due to the nonlinear increase in shear stress with a linear increase in the magnetic field strength.
Piston velocity is another critical parameter in damper design.The emulation of viscous damping with MR dampers needs to increase the coil current more or less in half-sine shaped mode every half period: See reference [19], [20] The criticality of geometrical and non-geometrical parameters was identified using Design of Experiments (DOE) and the critical parameters were optimized using pattern search approach which is explained in section 6. Results of the experimentation show that the coil current is the primary variable contributing to the desired variation in force.From the design perspective, the damper is designed for a specific operating range of the current between 0 A to 2 A [1], [21], [22].
The MR fluid considered for this study is MRF-132EG [1], [6] which is the common fluid used for automobile applications due to its low apparent viscosity when the fluid is not electromagnetically excited.The maximum piston velocity is derived from the desired spring stiffness, the road input, and ride frequency.
Road disturbances can be of any form, including sinusoidal, triangular, rectangular periodic waves, or a random combination of all of the disturbances.However, all these signals can be modelled as a homogenous sinusoidal signal using Fourier and Laplace approximations [10].The maximum disturbance producing sine wave input for the vehicle is considered for this model and the maximum piston displacement is set equal to the maximum spring deflection [10].
To determine the velocity, the road bump amplitude is modelled to be between 5 cm to 10 cm [8].Using the equations 23-25 for harmonic motion for a road amplitude of 7 cm and the frequency of 1.2 Hz, the maximum piston velocity is determined to be 0.5 m/s.
The lengths of magnetic links (L1-8) are initially calculated from the geometry shown in figure 14 and subsequently the cross-sectional areas (A1-8) of the links are calculated.Using equation 26, the total magnetic reluctance of the magnetic circuit is calculated as 0.3967 A.T/Wb [1].
Using equation 28, the magnetic flux density (B) in the circuit is calculated as 1.3039 Tesla for a current of 2 A [1].
The term µ0 refers to the permeability of vacuum which is 4π × 10 −7 H/m.

Magnetic Flux Density
Since the magnetic flux density in the flow gap primarily contributes to damping force changes, only the flux generated in the gap is considered for damping force calculation.The flux density in the steel components does not contribute towards shear stress for development of the damping force and hence can be neglected.
Using equations 29-30, the flow rate is calculated from the piston velocity and the piston shaft area to be 2.89 × 10 −4 m 3 /s [1].
The shear stress in the damper is determined as 83.8 kPa, using the equation 31 and MRF-122EG fluid datasheet [6].
The viscous and current dependent damping force components (F v ) and (F τ ) respectively, are determined using equation 32-33 [1].Finally, the total damping force is determined as 1.998 kN using equation 34.

V. DAMPER DESIGN OPTIMIZATION AND SENSITIVITY ANALYSIS
A. OPTIMIZATION APPROACH MR fluid constants of hysteresis equations were obtained from the fluid properties and other hysteresis constants were determined through simulation iterations for the designed spring mass damper system.To optimize output damping force, an integrated computation approach is used for each critical design parameter.The total damping force is defined as the objective function with damping force maximization criterion subjected to bounds on critical design parameters.
The initial values of these parameters were adopted from the calculations and literature review [1], [6], [12].The parameters in the constraint function are assigned minimum and maximum bounds for the generation of desired total damping force in the range of 826-2150 N which is calculated in the analytical model.The optimization process incorporates the use of the pattern search methodology to determine the best possible values of parameters which satisfies the damping force conditions for minimum and maximum values of current through the coil.A pattern search method of optimization is useful for optimizing discrete functions, which is desirable for the MR damper analysis.
The damping force is a function of piston radius (R), gap thickness (g), pole length (L), cylinder thickness (t) and piston internal radius (R c ).The optimization is performed to maximize the damping force given by equation 35-36 constrained by inequalities given by 37 [12].
The total damping force (Fmr), shear stress (τ ) and dynamic force range (D) are constrained with respect to the flux density (B), current (I) and each geometric parameter given by G. G is a set of critical geometrical parameters under consideration which are piston radius, piston core radius, pole length, gap thickness and cylinder thickness.
The optimization process yields a combination of a set of values for the critical parameters using applied voltage of 5 Volts as shown in tables 2 and 3.The iterations are performed simultaneously on each variable, starting with the lower bound on variables, and the process is stopped when the desired damping force in both inactive and active states is achieved.
Furthermore, a sensitivity analysis for parameters is performed at the minimum and maximum values of currents passing through the electromagnetic coil to find the extent of effect of optimum parameters on the damping force.

B. SENSITIVITY ANALYSIS
During the damper's inactive state, the current passing through the electromagnetic coil is zero and the damper should provide a viscous damping force of 826 N. A tolerance of ±5% is considered for the sensitivity analysis.Figures 15 and 16 represent the sensitivity of all the parameters, at a coil current of 0A and 2A respectively.Gap thickness, as presented in the graphs, has a higher sensitivity since the slope is steeper as compared to other parameters during both inactive and active operating conditions.This phenomena is observed in the optimization process as shown in tables 2 and 3.

C. SENSITIVITY RESULTS
The optimization process successfully evaluated the optimal values of the five critical geometric design parameters for both inactive and active states: Piston Radius (R) = 21mm Pole Length (L) = 14.5mmGap Thickness (g) = 1.3mmCylinder Thickness (t) = 9mm Internal Piston Radius(R c ) = 10mm   These optimized geometric parameters satisfy the total damping force requirements in compliance with the constraints.
The parameters are validated by modelling the damping force as a function of current range of 0-2 A. The minimum and maximum damping force obtained at this range including a 5% sensitivity tolerance are 798N and 2205N respectively for the desired operating range of 826 N to 2100N.The damping force increases exponentially as a function of the input current as shown in figure 17.This is due to the nonlinear relation between the magnetic flux density and shear stress of the MR fluid.
The models are analyzed using five different road profiles to determine the most efficient and accurate model with minimum overshoot, least settling time, and maximum damping factor for the quarter car semi-active suspension design.Tables 4 and 5 show the comparison of all the models with passive suspension.The statistics show that the Bouc-Wen model has a highest improvement and is the superior model compared to all other models considering sprung mass acceleration, overshoot, settling time, logarithmic decrement of spring mass vibration, and damping factor.

CONCLUSION AND FUTURE SCOPE
Present research successfully compares the analytical design of the passive suspension and semi-active suspension using the Bingham hysteresis, Dahl hysteresis, and the Bouc-Wen hysteresis models in conjunction with the MATLAB/Simulink package.
The results show that the Bouc-Wen model is most efficient and the most apt model for the design of the semi-active suspension system.Further analysis is performed to determine the required maximum damping force under the utmost road disturbance to wheel.The MR damper is successfully designed and analytical results are then verified with the MATLAB simulation results, which show 91% agreement between two.The geometric parameters of the damper are then optimized to develop the maximum required damping force and multi-objective optimization is successfully done with the pattern search approach.Sensitivity analysis is performed to determine the sensitivity of all five critical parameters and vindicated with the simulation results that gap thickness is the most sensitive parameter for the damper design considering the maximum damping force requirement.Future scope for this research is implementation of the conventional controller such as PID controller as well as optimal controllers such as LQG and LQR controller for the developed models.

FIGURE 3 .
FIGURE 3.Comparative response of all the models for step input.

Figure 3
Figure 3 shows that Bouc-Wen model has minimum overshoot and the lowest settling time as compared to the

FIGURE 4 .
FIGURE 4. Magnitude and phase analysis of all the models for step input.

FIGURE 5 .
FIGURE 5.Comparative response of all the models for sine input.

Figure 6
represents the magnitude of the displacement and signal phase change.The Bouc-Wen model follows the road profile with the least vertical displacement of vehicle sprung mass.

FIGURE 6 .
FIGURE 6. Magnitude and phase analysis of all the models for sine input.

FIGURE 7 .
FIGURE 7.Comparative response of all the models for white noise input.

FIGURE 8 .
FIGURE 8. Magnitude and phase analysis of all the models for white noise input.

FIGURE 9 .
FIGURE 9. Comparative response of all the models for uniform random number input.

FIGURE 10 .
FIGURE 10.Magnitude and phase analysis of all the models for uniform random number input.

FIGURE 11 .
FIGURE 11.Comparative response of all the models for mixed sine wave and uniform random number input.

FIGURE 12 .
FIGURE 12.Magnitude and phase analysis of all the models for mixed sine wave and uniform random number input.

FIGURE 13 .
FIGURE 13.Damping force requirement vs time for the MR damper design.

TABLE 4 .TABLE 5 .
Damping force statistic with simultaneous changes in the parameters at passive off (I = 0A).Damping force statistic with simultaneous changes in the parameters at passive on (I = 2A).

TABLE 2 .
Damping force statistic with simultaneous changes in the parameters at passive off (I = 0A).

TABLE 3 .
Damping force statistic with simultaneous changes in the parameters at passive on (I = 2A).