Frequency Dependent Spencer Modeling of Magnetorheological Damper Using Hybrid Optimization Approach

Magnetorheological dampers have been widely used in civil and automotive industries. The nonlinear behavior of MR fluid makes MR damper modeling a challenging problem. In this paper, a frequency dependent MR damper model is proposed based on Spencer MR damper model. The parameters of the model are identified using an experimental data based hybrid optimization approach which is a combination of Genetic Algorithm and Sequential Quadratic Programming approach. The frequency in the proposed model is calculated using measured relative velocity and relative displacement between MR damper ends. Therefore, the MR damper model will be function of frequency. The mathematical model is validated using the experimental results which confirm the improvement in the accuracy of the model and consistency in the variation damping with the frequency.


Introduction
The semiactive control system provides both features of passive and active devices in terms of reliability and adaptability.Using semiactive system, the rate of energy dissipation becomes controllable, while, in active control devices, the energy can be added to the system to control the dynamic response.Magnetorheological (MR) damper is a semiactive control device which is commonly used in vehicle industries and structural applications.The MR damper contains MR fluid instead of regular oil.The MR fluid is a smart material which contains micron sized magnetic polarized metal particles which provide variable viscous damping with changing magnetic field [1].The application of MR fluid is dependent on three different operational modes: flow, squeeze-flow, and shear [2].For instance, MR dampers and servovalves are designed based on flow mode of MR fluid [3].The squeezeflow mode of MR fluid is utilized in the application of impact control dampers for large forces [3].The shear mode can be used for brakes, clutches, and damping layer of sandwich structures [3].In order to describe the dynamic behavior of the MR damper, different mathematical models have been proposed in both discrete and continuous time domain.
Modeling of MR damper by using black box nonlinear models is carried out in discrete time domain [4] where MR damper hysteresis, function of displacement, and velocity are modeled by Neural Network (NN).In this model, three parameters are indicated based on the collected experimental input and output data.However, the nonmodel based parameter identification is only valid for the system operation range during which the experimental data are collected and used for training the NN model.Therefore the accuracy of model cannot be guaranteed for extrapolating the range of operation.
The Bingham viscoplastic MR damper model [5] is built in continuous time domain and it describes the dynamic behavior of MR damper based on measuring the shear stress and the shear strain rate.The Bingham model consists of a Coulomb friction element parallel with viscous damping.Using Bingham model, the storage energy in the MR damper cannot be modeled.Moreover, the difference between the simulated force and the real force increases when the velocity is near zero.
The modified Bingham MR damper model proposed by Gamota and Filisko [6] is a viscoelastic-plastic model.The socalled Gamota and Filisko model is a Bingham model which is in series with a parallel set of spring and viscous damper.The Gamota and Filisko MR damper model improves the accuracy of the model in describing the hysteresis loop and storage energy of the MR damper.However, the simulation of this modified Bingham model needs step size in the order of 10 −6 which is the main drawback of this model [1].
The Bouc-Wen model [7] is a continuous, viscoelasticplastic model which can describe a wide range of hysteresis behavior [1].The hysteresis behavior of Bouc-Wen model is described by an evolutionary variable with three coefficients of velocity which results in smoothness of transition from the preyield to the postyield region.The roll-off effect cannot be simulated using Bouc-Wen model in the region of the small magnitude of the velocity where the velocity and acceleration have opposite directions [1].
In the Spencer MR damper model [1], a spring and a viscous damping element are added to the Bouc-Wen model to simulate the roll-off effect at small velocities.Therefore, the other pair of damping and stiffness elements can be adjusted for small velocities or high frequency region.The Spencer model is capable of simulating the roll-off effect in all velocity and acceleration regions.In the Spencer model, the assigned damping coefficients only depend on the changing current.However, the MR damper viscosity depends on the frequency of excitation [8] and temperature of MR fluid [9].And in the literature, MR damper models cannot describe such frequency dependent behavior.
The present study deals with the variation of MR fluid viscosity with the frequency of the excitation in Spencer MR damper model.The viscosity of the MR damper is modeled using two viscous damping elements for large and small velocities.The MR model is identified by minimizing the error between the experimental data and simulated data of the proposed model.The hybrid optimization approach is used for the identification which is a combination of Genetic Algorithms (GA) and Sequential Quadratic Program (SQP).The excitation frequency in real application can be calculated by measuring the velocity and the displacement of the MR damper.Therefore, the viscosity of MR fluid in Spencer model is described by exponential and Gaussian equations which are the functions of velocity and displacement for small and large velocity regions, respectively.
The rest of the paper is organized as follows.The modeling of MR damper is presented in Section 2. The experimental set-up and procedure are explained in Section 3. The optimization approach and characterization of MR damper model are discussed in Section 4. A comparison between the proposed model and experimental data is presented in Section 5. Finally the conclusions and future work are presented in Section 6.

Spencer Magnetorheological Damper Model
Due to the nonlinearity in dynamic behavior of MR damper, the accuracy and validity of the Spencer MR damper model over wide range of frequencies are not consistent.The schematic of the Spencer MR damper model is shown in where  mr1 and  mr1 are the MR damper accumulator stiffness and viscous damping coefficients for small velocities, respectively. 0 is the initial displacement associated with spring  mr1 .Further,  mr0 and  mr0 are accumulator stiffness and viscous damping coefficients for large velocities, respectively.Equation ( 2) represents the internal state  which is used to define the roll-off due to the damping coefficient  mr1 .The MR damping force in (1) involves the nondimensional auxiliary variable  to define the hysteresis.The constants  and  are the nondimensional values to present the hysteretic loop in the negative and positive slopes in (3). describes the hysteresis loop size with respect to velocity.The scalar value  is used to represent the smoothness of transition of MR fluid from elastic to plastic [1,10].
The  variable in ( 2) is a polynomial function of voltage to describe the MR fluid yield stress which is defined as a first order polynomial in (6).The variations of  mr1 and  mr0 with respect to the frequency of excitationare proposed by  mr1 (Δ, Δ Ż ) and  mr0 (Δ, Δ Ż ) in ( 4) and ( 5), respectively.The damping coefficients  mr1 and  mr0 are also functions of voltage as a first order polynomial in ( 4) and ( 5), respectively.Equation ( 7) presents a filter to reach rheological equilibrium.Accordingly, where  is constant for changing rate of magnetic field and  is the applied voltage.The peak value of relative velocity and peak to peak value of relative displacement of MR damper ends are Δ Ż and Δ, respectively.

Experimental Setup
In  2 and 3.The experimental data are gathered in two categories: displacement, velocity, and force for variation of frequency in range of 1.5 to 5 Hz in 0.5 (Hz) intervals.
The operating range of MR damper RD 8041 is 0 to 2 A and maximum temperature is 71 ∘ C. By increasing the frequency of excitation, the temperature of the MR fluid inside MR damper increases and is measured by a digital thermometer attached to the body of the MR damper as shown in Figure 3. Considering this constraint, the applicable maximum frequency is 5 Hz.The amplitude of the harmonic excitation is 0.00635 m.

Characterization of Magnetorheological Damper Model Using Hybrid Optimization Approach
In order to identify the parameters of Spencer model and to study the trend of viscosity versus frequency, the objective function is defined as the error between experimental data and mathematical data calculated using Spencer model.The objective function is defined in (8).The model variables defined in (1) to (7)  steps: (i) identifying ten parameters defined in (1) to (7) using experimental data at frequency of 2.5 Hz and zero volt and (ii) identifying  mr1 (Δ, Δ Ż ) and  mr0 (Δ, Δ Ż ) for frequency from 1.5 to 5 (Hz) in 0.5 (Hz) interval at zero current.Therefore, in the first step, the model has ten design variables and ten lower bounds and upper bounds as linear constraints which are presented in the following inequalities: The lower bound and upper bound of all design variables are given in Table 1 based on identified parameters for the MR damper in the same range of the force as in [11] considering ±50% variation.A nonlinear constraint is defined in (11) to limit the feasible search region and to guarantee that the error between experimental data and mathematical model is less than 4%: In order to identify the design variables  mr1 (Δ, Δ Ż ) and  mr0 (Δ, Δ Ż ), eight hybrid optimization algorithms are formulated based on eight experimental data sets in the range of  The objective function is defined as the error between the experimental and mathematical model as presented in (8).Linear constraints are the lower bound and upper bound and nonlinear constraint is as defined in (11).The quantities of lower bound and upper bound are assigned based on the sensitivity of the objective function to the design variables.
The contour plot of an optimization algorithm in Figure 4 presents the sensitivity of objective function to design variables.
The range of design variables in Figure 4 is considered to be ±40% from the identified parameters based on 2.5 Hz (first optimization step), in order to identify viscosity parameters of MR damper at frequency of 1.5 Hz.The lower bound and upper bound for data sets in other frequencies (in the range of 1.5 Hz to 5 Hz) are considered with ±20% variation with  respect to the assigned viscosities from previous data set.For example, the lower bound and upper bound of optimized algorithm for 2 Hz are defined based on ±20% from the identified  mr1 (Δ, Δ Ż ) and  mr0 (Δ, Δ Ż ) for 1.5 Hz.
Figure 4 shows that the variation of the objective function with respect to the design variable is convex, which guarantees that a global optimum point will be reached in the process of identification.It should be noted that the sensitivity of objective function with respect to viscosity parameters in 20% variation for all data sets is studied and convexity in the feasible region is checked.

Sequential Quadratic Programming Technique.
The Sequential Quadratic Programming is a methodology for nonlinear optimization problems considering nonlinear equality and inequality constraints.The SQP is an iterative procedure based on Quadratic Programming (QP) to solve QP subproblems and to define new iterations [12].The QP should be solved to satisfy feasibility of problem considering local properties of current iteration.Objective functions and linear and nonlinear constraints of QP subproblems are defined in the following [12]: In nonlinear optimization problems, the SQP method may find one of the extremum points which is near the selected initial point.Therefore, the sensitivity of the initial point with respect to the optimized parameters should be studied.The presented objective function is nonlinear and it is not robust in variation of initial point.Therefore, the GA is employed to assign the initial point in SQP method which makes up the hybrid GA-SQP optimization algorithm.

Genetic Algorithm Optimization Technique.
The GA is a multivariable optimization method for solving linear, nonlinear, discrete, continuous, differentiable, and nondifferentiable constrained problems by using stochastic search algorithms [12][13][14][15].The Genetic Algorithms are a random search based method which can avoid stopping in local optimum point.However, using GA the globality of the optimum point cannot be proved or guaranteed.The optimization problem for identifying the MR damper parameters is highly dependent on the initial point.The combination of GA and SQP can solve the problem by using initial point obtained through GA based on random search in SQP algorithm [16,17].The hybrid of GA-SQP approaches can reach global optimum point.But, the globality cannot be proved mathematically.
The GA operates based on random point selection in random generation [14].Therefore, running the same algorithm each time may give different optimal points [14].As a result, by increasing the number of populations and repeating the same algorithm, the relaxation of the obtained initial points is studied.The GA is designed based on the presented objective function (8) and constraints (10).

Results and Discussions
Two optimization approaches are employed to investigate the variation of viscosity of MR fluid with frequency.In the first algorithm, the parameters of MR damper are identified for 2.5 Hz and zero current.For the second part, the variations  0 and  1 are identified for the frequencies in the range of 1.5 Hz to 5 Hz at step of 0.5 Hz in eight suboptimization algorithms.

Parameter Identification for Constant
Frequency.The parameters of Spencer MR damper model are identified using hybrid optimization of GA-SQP.The identified parameters are presented in Table 2 based on experimental data at 2.5 Hz and zero current.The initial values for SQP are assigned using GA.The GA algorithms are used for different number of populations from 80 to 360 in the interval of 80.Moreover, to study the result of relaxation, the algorithms are repeated ten times for each population.The presented results in Table 2 are obtained based on 240 numbers of populations resulting in minimum error.The identified parameters in this table are used for second step optimization in order to find the viscosity variation of MR damper at different frequencies.Figures 5 and 6 show the force velocity and force versus time, respectively.The objective function is defined as the absolute value of the error between the experimental measured force and the generated force in mathematical model.Figure 6 demonstrates that the mathematical model is fitted by experimental data with 2.85% error.However, due to the neglected effect of the velocity in objective function (7), the experimental and mathematical force-velocity curves have some discrepancies in Figure 5.

Parameter Identification for Variable
Frequency.The variation of  0 and  1 for the frequencies in the range of 1.5 Hz to 5 Hz in interval 0.5 Hz is calculated by solving eight suboptimization problems.The optimization objective is to minimize the error between the experimental and mathematical data.The constraints are defined in the range of ±40% variation from the optimized values of  0 and  1 for each frequency from the assigned viscosities in 2.5 Hz.The other linear constraints are assigned considering ±20% variation with respect to the assigned viscosities in the previous step.The sensitivity of the objective function to  0 and  1 in the feasible search region is examined to guarantee the convexity of the function in feasible range defined by lower bound and upper bound.The suboptimization problems are solved using hybrid optimization technique.
Figures 7 and 8 show the variation of  0 and  1 for different frequencies by blue points.The variation of  0 is modeled by an exponential function which can be used for extrapolation.Therefore, the variation of  0 is defined as a function of frequency in (13).The frequency can be measured by displacement sensor, and the relation between velocity, displacement, and frequency is presented in (15).The variation of  1 with frequency in Figure 8 is modeled using the second order Gaussian function in (14): 1 () = 7087 ((4.282−)/1.197) 2 + 18520 ((1.783−)/3.25) 2 , ( 14) Figure 7 shows that parameter  0 is sensitive in the low frequency region and becomes constant at high frequency, whereas parameter  1 in Figure 8 is constant at low frequency and becomes variable in the high frequency region.In order to compare the effect of frequency dependent model and frequency independent model (existing model), the error between the experimental and analytical models is presented in Table 3 for seven frequencies.The parameters of the frequency independent model are identified at frequency of 2.5 Hz.Therefore, the errors of both frequency independent and frequency dependent models for data samples near 2.5 Hz are close to each other.However, at high frequencies, the difference becomes significant.For example, for 5 Hz data sample, the error is decreased by 1.1%.The extrapolation of viscosity variation in Figure 8 shows that the difference at high frequencies would be significant.
Figure 9 shows that the existing independent frequency model has overestimation in modeling for hysteresis of MR damper in comparison with the proposed frequency dependent model.The main objective in modelling the MR damper is to simulate MR damper force at different frequencies of excitation and current.It shows that the proposed frequency dependent model improves the accuracy of force simulation.Moreover, the proposed model is in agreement with the force simulation in variable frequency and operating conditions of the MR damper.

Conclusions
In this study, the variation of MR damper viscosity at different frequencies has been studied using analytical and experimental approaches.

Figure 1 .
Figure 1.The governing equations of Spencer MR damper model are presented in the following [1]:

Figure 4 :
Figure 4: Sensitivity of objective function to variation of viscosity parameters of MR damper in 1.5 Hz.

Figure 5 :
Figure 5: Force-velocity hysteresis in 2.5 Hz and zero current.

Figure 6 :
Figure 6: Experimental and mathematical MR damper force in 2.5 Hz and zero current.
=  mr0 ,  mr1 ,  0 ,  1 ,  0 ,  1 ,   ,   , , , , ,  0 .(9)Theconstraints are defined to limit the search region.The linear constrains are defined as lower bounds and upper bounds in both GA and SQP optimization algorithms.The identification of the model parameters is carried out in two

Table 2 :
Identified parameters of MR damper in 2.5 Hz and zero current.

Table 3 :
Comparison between frequency dependent and original Spencer model.
The hybrid GA-SQP optimization techniques are used to identify the parameters of the MR damper based on experimental data.The variation of MR damper viscosity with frequency is studied for frequency range from 1.5 HZ to 5 Hz at 0.5 Hz interval using hybrid GA-SQP technique.A mathematical MR frequency dependent model is proposed based on Spencer MR damper model where the viscosity of MR damper is modeled by exponential and Gaussian functions.The results confirm that the frequency dependent MR damper model improves the accuracy of the model in force simulation at high frequency region and shows consistency in force simulation as well.The proposed MR damper model can be used in application of impact control and aircraft landing gear in which these operational conditions need high precision in high frequency region.