A Piecewise Hysteresis Model for a Damper of HIS System

A damper of the hydraulically interconnected suspension (HIS) system, as a quarter HIS, is prototyped and its damping characteristic is tested to characterize the damping property.The force-velocity characteristic of the prototype is analyzed based on a set of testing results and accordingly a piecewise hysteresis model for the damper is proposed.The proposed equivalent parametric model consists of two parts: hysteresis model in low speed region and saturation model in high speed region which are used to describe the hysteresis phenomenon in low speed and nonhysteresis phenomenon in high speed, respectively. The parameters of themodel are identified based on genetic algorithm by setting the constraints of parameters according to their physical significances and the corresponding testing results. The advantages of the model are highlighted by comparing to the nonhysteresis model and the permanent hysteresis model. The numerical simulation results are compared with the testing results to validate the accuracy and effectiveness of the proposedmodel. Finally, to further verify the proposedmodel’s wide applicability under different excitation conditions, its results are compared to the testing results in three-dimensional space. The research in this paper is significant for the dynamic analysis of the HIS vehicle.


Introduction
As the suspension is the key part influencing the vehicle posture, ride comfort, road holding, and vehicle agility [1], the accuracy of its model is important for precisely predicting the properties of the vehicle.In past years, with the development of advanced suspensions, such as semiactive suspension and active suspension, the corresponding shock absorber models are developed [2,3].HIS is a novel kind of passive suspension whose chambers of different dampers are interconnected mutually by hydraulic pipes so that the damping coefficients in different motions are able to be independently tuned [4,5].Smith et al. developed the mathematical model of HIS, in which the characteristic of damper valves, however, is assumed to be linear [5,6].To improve the accuracy of the model of HIS, a more precise and complex model for the damper in HIS is required.
Models of passive dampers have been developed over the last decades, which can be classified into parametric and nonparametric types [7].Physical parametric models indicate these models which consider the oil flow rate and the valves deformation inside the damper so that the fluid equations are coupled to the mechanical equations.Accordingly, this kind of model can describe the mechanism of the system in detail, for example, the nonlinear physical parametric model containing 82 parameters developed by Lang [8] and the nonlinear hysteretic parametric models developed by Reybrouck et al. [9,10].These models are able to be used as subroutines in the full-car simulation but its complexity also sets a barrier for the parameter identification and limits its practicability.Differently, the equivalent parametric models are relatively much simpler and they are the equations derived from fitting the testing results and built by ignoring some practical structures including operating conditions inside the damper.The simplest equivalent parametric models are the linear model and speed piecewise linear model which are too simple to precisely predict the damping force.An equivalent parametric model consisting of a dashpot, a spring, a backlash, and friction elements was developed by Karadayi and Masada, which can fit the testing results of low frequency well and describe the hysteresis phenomenon by considering the damping force as a function of both displacement and velocity [11].A damper model with hysteresis loop based on Matlab/Simulink was developed by Ko et al. and used in full-car simulations in their research; however, the mathematical model was not shown explicitly and it is observed that the hysteresis permanently occurs according to the displayed figures in [12].
The nonparametric model based on an alternative restoring force method gets more attention in the community of testing modal analysis [7].The results of [13] show that the nonparametric model as a function of both velocity and acceleration predicts more precisely than the nonparametric model based on both displacement and velocity.
Recently, with the development of the semiactive shock absorber and for the excellent properties of magnetorheological (MR) damper such as rapid dynamic response, low power consumption, and high-fidelity control, the research to modeling MR damper is increasing [14,15].Bouc-Wen model is the most popular model for MR dampers which has been developed by many other researches [16,17] and applied in various kinds of systems [18].
The structure of this paper is as follows.In Section 2, the HIS system and the prototype of a quarter HIS system are described and the test setup is introduced.The testing results are discussed to characterize the damping property in Section 3 and accordingly the model is built and the parameters of the proposed model are identified in Section 4. In Section 5, the effectiveness of the proposed model is validated in various ways.The conclusion is drawn in Section 6.

Prototype and Test Setup
2.1.Description of the Prototype.HIS system is the kind of suspension in which the dampers of different wheels are interconnected by pipes as shown in Figure 1 [5].However, as the test rig is only able to characterize one damper, we have to separate this system and build a quarter HIS system to be tested as shown in Figure 2. Different from the conventional damper, there is no damping valve on the piston of HIS damper but, as with the fluid cylinder, there are ports where damping valves are installed.The quarter HIS prototype employs a single damper of HIS where the ports of two chambers are connected by pipes.To accommodate to the variation of the inside volume of the system while the piston rod moves in or out, an accumulator is connected to the pipe via a tee joint.
As shown in Figure 2, when the damper is at compression stroke, the piston moves downward and the fluid flows out from compression chamber into both the rebound chamber and the accumulator by compressing the gas in it.Compression spring valves can be easily opened at compression stroke but the compression disc valve can only be opened while the downward speed of piston is sufficiently high, which results in decreased damping which prevents the force from being too great.At rebound stroke, the piston moves upward.The fluid flows back into compression chamber via only the orifice and rebound check valves, until the upward speed of piston is high enough to open the rebound spring valve as well.

Test Setup.
As shown in Figure 3, a computer-controlled MTS is employed to drive the damper and test the characteristic of the damper.The piston rod is fixed on the upside of the hydraulic actuator while the damper body is fixed on and driven by the bottom of the equipment.
To eliminate the influence of the temperature, both the prototype and the MTS test rig were warmed up for 1 minute with the excitation peak speed of 0.5 m⋅s −1 and the amplitude of 45 mm before measuring the data.It is also worth mentioning that the volume of gas in the accumulator changes with the piston rod moving in or out, which can influence the pressure in the system.In addition, because of the different area of two sides of the piston, the pressure in the system can result in an additional rebound force.Considering the purpose of this research is only to model the property of the damper, the effects of the variation of the gas volume inside the accumulator, that is, the variation of the mean pressure caused by the different relative displacements between the piston and the body of the damper, should be eliminated.To this end, a test under slow excitation with a peak velocity of 0.005 m⋅s −1 is implemented to investigate the effects of relative displacement and the testing results are shown in Figure 4.As the excitation speed is relatively slow compared to that in other tests, it is reasonable to consider that the force in this test is completely produced by the relative displacement and not influenced by the damping valves.To eliminate such effects, the forces are subtracted from other testing results according to the corresponding displacement data in the same direction of velocity.The data shown in the following sections has been processed and for brevity it is not going to be recurrently mentioned.
In the tests for modeling, the imposed signals are sinusoidal displacements with 45 mm amplitude and their peak velocities for different tests are varied from 0.1 m⋅s −1 to 1.0 m⋅s −1 with an increment of 0.1 m⋅s −1 .To validate the suitability of the proposed model, the tests in other excitation conditions are also implemented.More detail about them is shown in Section 5. A force sensor is used to measure the damping force.A displacement sensor is used to obtain the displacement.Based on the difference derivation, the corresponding velocity and acceleration can be obtained according to sampling force data and time data.

Discussion of Testing Results
The peak force-peak velocity values of a series of the testing results for modeling are noted in Table 1 and, in order to clearly observe the hysteresis, the curve drawn by connecting the points of peak force-peak velocity is drawn together with sets of testing results in Figure 5.As the maximal force or the minimal force is not always at the maximal velocity point or the minimal velocity point, respectively, it is noted that peak force here refers to the force at the peak velocity point.The positive velocity refers to the rebound while the negative velocity refers to the compression.
In Figure 5, the testing results in three different tests are compared with the peak force-peak velocity relation.In low frequency range as shown in Figure 5(a), when the peak velocity of excitation is 0.1 m⋅s −1 , the speed of piston Figure 1: Roll-plane HIS system [5]. is not high enough to open the high-pressure valves (the compression disc valve and the rebound spring valve) and the fluid flows via only the orifice and low-pressure valves (the compression spring valve and the rebound check valve).The phenomenon of hysteresis is negligible in this case.However, when the peak velocity reaches 0.5 m⋅s −1 as shown in Figure 5(b), the hysteresis phenomenon becomes obvious in the region of low speed.The authors suppose that the main reason for the hysteresis is that the opening or closing of the high-pressure valves will take a period of time during which the piston speed has changed significantly.However, in the region of high speed, the high-pressure valves keep ultimately open and accordingly the damping coefficient keeps constant and thus the hysteresis is hard to be observed.With the increase of frequency, as shown in Figure 5(c), when the peak velocity of excitation reaches 1.0 m⋅s −1 , due to, for example, cavitation and abnormal inside leakage, the fluid cannot return quickly enough to the low-pressure chamber, which results in the interaction force between piston and the fluid to be around zero for a short period of time.During the restoration period from around zero to normal state, the characteristic is distorted and the force sharply fluctuates, which is mainly caused by the imperfect sinusoidal excitation under the great force.

Modeling and Parameter Identification
To detect the common characteristics, the testing results under excitation frequencies of 1.0610 Hz, 1.7684 Hz, 2.8294 Hz, and 3.5368 Hz and the curve of peak forcepeak velocity are drawn together in Figure 6.For easy identification, the curve of sine 3.5368 Hz-45 mm is marked by letters from A to J in chronological order.
The model is divided into two parts: the boundary of saturation and the hysteresis in low speed region.The saturation and time delay are employed to describe the two parts, respectively, which is described in detail as follows.

Saturation Model.
As shown in Figure 6, we can also find that the damping force seems to be limited in a diagonal band and never goes beyond the white range unless distortion occurs, which is named as "saturation" in this paper.The boundaries of the saturation can be generally divided into two parts: maximum boundary A-B-C  -D ( = 1, 2, 3, 4) (A-D) and minimum boundary F-G-H  -I ( = 1, 2, 3, 4) (F-I) which are modeled separately.
If we use  max and  min to describe the maximal boundary and the minimum boundary, respectively, at the moment , the saturation model can be formulated as follows: where  max and  min are modeled in the following content.As shown in Figure 6, it is obvious that the maximum boundary is the part A-D.For this part, the slope of curve A-B almost remains constant, which represents the damping coefficient while all rebound valves are open.The slope of curve B-C  ( = 1, 2, 3, 4) shows a downward trend with the increase of the excitation frequency, which is caused by tendency of the rebound spring valve to close.Also of note is that, irrespective of the peak velocity, the horizontal line C-D around zero force can keep the force nonnegative, which represents the piston travel in practice caused by the cavitation.Thus, the maximum boundary of saturation can be divided into three regions: part A-B with constant slope, part B-C with variable slope, and the horizontal part C-D.
Part A-B with constant slope is represented by the constant damping coefficient  1 .Considering the produced force always passes the position of point B, the force of part A-B can be represented by  1 ( ẋ () − ẋ B ) +  B while ẋ () is greater than ẋ B .ẋ () is the relative velocity between the body and the piston of the damper, and its positive value represents the rebound stroke and negative value represents the compression stroke.The horizontal part C-D denoting the piston travel is represented by a constant  C-D which is around zero.It is necessary to point out that the horizontal line C-D does not work until the force generated by hysteresis model is hysteretic enough to touch it in the situation of high frequency excitation.The force generated by hysteresis model will be discussed later.
To sum up, the maximum boundary can be represented as follows: Similarly, for the minimum boundary F-I, it can be divided into three regions: part F-G with constant slope, part G-H with variable slope, and the horizontal part H-I.As the characteristic at the compression stroke is different from that  at rebound stroke, it is required to be modeled separately.The original value of  min , the force of part F-H, can be calculated according to the excitation velocity and acceleration firstly.If the value-calculated  min original is greater than  H-I , it is set as  H-I , as follows: min original ,  min original <  H-I . (3)

Hysteresis Model.
The hysteresis refers to misalignment between part D-E and part J-I in Figure 6.As shown in this figure, the slopes of the two curves D-E and J-I remain stable until the excitation reaches such a high frequency that the distortion of damping characteristic occurs.The degree of hysteresis shows an upward tendency with the increase of the excitation frequency.The authors suppose that the time taken to open and close the high-pressure valves is constant, so the degree of hysteresis is related to the variation of velocity over the period of opening or closing high-pressure valves.Accordingly, it is acceptable to use a linear model with a time delay to describe the hysteresis in low speed region unless the distortion happens.Constant  7 is employed to denote the damping coefficient for the condition that only the compression spring valve and the rebound check valve are opened in low speed region.Constant  is used to denote the time interval from the force reaching the threshold of being able to open or close the high-pressure valve to the valve being ultimately opened or closed.The hysteresis model can be formulated as follows: where  is the excursion comparing with zero point.

Integrative Model.
The whole model integrates the saturation model and the hysteresis model by calculating the force of the hysteresis model firstly and then limiting it by setting the boundary via the saturation model.The integrative damping force  is calculated as shown in (5). max ,  min , and  delay are calculated by ( 2), (3), and ( 4), respectively: delay ,  min ≤  delay <  max ,  max ,  max ≤  delay . (5)

Identification Procedure.
As the functions are already given, the values of the parameters in them need to be identified based on the commonly used genetic algorithm.Optimized variables and objective function are the two necessary factors for the optimization while reasonable constraints can markedly enhance the computational efficiency.
The interested reader can be referred to [19] for further information about the algorithm.Optimization variables in this paper refer to all coefficients in the integrative model: where test(  ) refers to the testing results,  refers to different tests, and  refers to the sample point in each test.
One of the advantages of the proposed model is that the applicable constraints can be generally obtained according to testing results in Figure 5 to boost the optimization efficiency.The rough position of point B (corresponding to ẋ B and  B ) can be observed, so the constraints for ẋ B and  B can be estimated and, similarly, the general position of point G (corresponding to ẋ G and  G ) can be observed, so the constraints for ẋ G and  G can be estimated. 1 relates to the slope of A-B.The constraints of  2 and  3 can be estimated according to the slope of B-C. 4 relates to the slope of F-G while the constraints of  5 and  6 can be estimated according to G-H;  7 can be assessed according to the slopes of both D-E and I-J in the situation of low excitation frequency without distortion.The value of  approximates the force corresponding to the velocity of zero in Figure 5(a). is the time of hysteresis and its maximal constraint is set to 0.2 s.The constraints and the achieved optimization results are shown in Tables 2 and 3, respectively.

Comparison between the Proposed Model and Other Models.
According to the testing results, the authors of this paper believe that, irrespective of the rebound or the compression stroke, the damping coefficient can remain constant while all the valves are completely open, and accordingly there is no hysteresis at this stage.However the model in [12] is always hysteretic during the operation which, for brevity, is named as the permanent hysteresis model in this paper.The model in ADAMS can be nonlinear but without hysteresis.Accordingly, the models both in [12] and in ADAMS are not perfectly suitable for the HIS damper, while the proposed model can be both nonlinear and hysteretic in the low speed region.In order to demonstrate the advantage of the proposed model, a qualitative comparison between the proposed model, permanent hysteresis model, nonhysteresis model, and the testing results is carried out.The employed nonhysteresis model is similar to the proposed piecewise hysteresis model, while the delay time in (4) is set to zero.The nonhysteresis model can be written according to (5) as follows: where the time of  delay is calculated  in advance so that the time delay in proposed model is canceled.The damping force calculated by permanent hysteresis model is   and its value can be derived according to   in (7) by setting the delay time , at all times, as follows: All three models and the corresponding testing results are demonstrated in Figure 7.This figure shows that while the permanent hysteresis model and the nonhysteresis model can only accurately predict the force in low speed and high speed region, respectively, the proposed model can make a relatively accurate prediction in both low speed and high speed regions.
In order to evaluate the degree of the discrepancy, both the simulation and testing results are sampled every 0.01 s to calculate the relative errors.The root mean square (RMS) of relative errors is used to evaluate the accuracy of models: where   is the sampling time and   (  ) and   (  ) are the sampled force in simulations and tests, respectively.The relative errors in the tests with different excitation frequency are demonstrated in Figure 8.As shown, the three kinds of models are of similar accuracy in low excitation frequency because both the permanent hysteresis model and the proposed model do not apparently show a characteristic of hysteresis in the low frequency excitation.With the increase of the excitation frequency, the proposed model can evidently predict the force more accurately.For all models, the relative RMS errors show an obvious upward trend with the increase of the frequency when the frequency is over 3.1831 Hz.This is caused by the undesired distortion of the damping characteristic while the force is too great.It is noted that all of the three models do not reproduce the damping force in low frequency as accurately as that in intermediate frequency, which is caused by the nonlinear damping  characteristic in the nearby speed region and needs to be improved in the future.Section 5.2 details more information about the comparisons between the damping characteristics of the proposed model and the testing results.

Comparison between Simulations and Tests.
The simulation results are compared to the testing results in Figure 9 where the amplitude of excitation displacement is fixed at 45 mm and the frequency ranges from 0.74 Hz to 3.54 Hz.As shown in the figure, most of the simulation results agree well with the testing results until the excitation frequency rises to reach 3.5368 Hz where the distortion is evident.It is pointed out that Figure 9 shows a discrepancy at low frequency of 0.74 Hz in rebound, which is caused by the undesired transition between the nonhysteretic part and hysteretic part.This needs to be improved in the future model.
In order to demonstrate the suitability of the model in a wide range of circumstances, the amplitude of excitation displacement of 34 mm and 17 mm was used.Figure 10 shows that the simulation results agree well with the testing results.

Validation of the Proposed Model
Based on the Nonparametric Model.In Section 5.2, the accuracy of the model is validated in the conventional way by giving examples of comparisons between the simulation and the test, which is limited because it is impossible to demonstrate all the damping characteristics under all excitation conditions.In order to enhance the validation of the suitability of the proposed model, the proposed model is compared to the three-dimensional model which is directly built according to testing results.
Duym found that the damping characteristic can be described by a curved surface in three-dimensional space and verified that the force-velocity-acceleration model is of higher accuracy than the force-velocity-displacement model [13].Therefore, the force-velocity-acceleration model is employed to validate the proposed model.
In the force-velocity-acceleration model, the curved surface of damping characteristic formed by testing results is displayed in both Figures 11 and 12. Figure 11 displays the simulation results under the same sinusoidal excitation conditions while Figure 12 shows the simulation results under different excitation conditions.As shown in Figure 11, the simulation results under the same excitation conditions agree well with testing results until the hysteretic characteristic is aggravated sharply caused by the distortion.Figure 12 demonstrates the comparison between testing results and simulation results under different excitation conditions.For clarity, (a) and (b) show it in two different visual angles, respectively.According to Figure 12(b), the accuracy of the proposed model deteriorates while the distortion of the damping characteristic occurs.It is worth noting that the distortion of the damping characteristic is undesired in practice and needs to be eliminated by optimizing the prototype in the future.

Conclusion
This paper developed a piecewise hysteresis model which consists of the saturation model and the hysteresis model for a prototype of the damper of a quarter HIS system and its parameters have been identified.The proposed model can accurately predict the damping force of a quarter HIS until the distortion of the damping characteristic of the prototype appears.This is significant for improving the accuracy of the full-car model with HIS system to predict the practical dynamic characteristics better.
The time delay is set as a constant which actually should be more complex and the transition between the nonhysteretic part and the hysteretic part should be nonlinear and smooth.This needs to be optimized in the future.Additionally, the undesired distortion of the damping characteristic also needs to be eliminated by optimizing the prototype.to gratefully acknowledge the help of Mr. David Burt in the final language editing of this paper.

Figure 3 :Figure 4 :
Figure 3: Test setup: (a) photo of the test rig; (b) damper sensors and the hydraulic actuator.

Figure 7 :
Figure 7: Different simulations and the corresponding testing results.

Figure 8 :
Figure 8: Relative RMS errors in various exciting frequencies with amplitude of 45 mm.

Figure 10 :l o c i t y o f t h e p i s t o n ( m • s − 1 )
Figure 10: Comparisons while the amplitude of displacement is (a) 34 mm and (b) 17 mm.

Figure 11 :Figure 12 :
Figure 11: Three-dimensional comparisons under the same excitation conditions.

Table 2 :
Constraints of parameters.

Table 3 :
The results of the optimization algorithm.