Research on Modeling of Hydropneumatic Suspension Based on Fractional Order

With such excellent performance as nonlinear stiffness, adjustable vehicle height, and good vibration resistance, hydropneumatic suspension (HS) has been more and more applied to heavy vehicle and engineering vehicle. Traditional modeling methods are still confined to simple models without taking many factors into consideration. A hydropneumatic suspension model based on fractional order (HSM-FO) is built with the advantage of fractional order (FO) in viscoelastic material modeling considering the mechanics property of multiphase medium of HS. Then, the detailed calculation method is proposed based on Oustaloup filtering approximation algorithm.The HSM-FO is implemented in Matlab/Simulink, and the results of comparison among the simulation curve of fractional order, integral order, and the curve of real experiment prove the feasibility and validity of HSM-FO.The damping force property of the suspension system under different fractional orders is also studied. In the end of this paper, several conclusions concerning HSM-FO are drawn according to analysis of simulation.


Introduction
FO calculus is a mathematical theory of investigating the differential and integral properties and application of any order.FO calculus is extension and promotion of integral order (IO) calculus, and they appeared almost simultaneously.The main difference between FO calculus and IO calculus lies in the global memory function and hereditary property.Not only is the order of calculus integers but also it can be real numbers and even complex numbers.In recent years, FO calculus has been widely applied to Non-Newtonian Fluid Mechanics, Multimedia Mechanics, Rheology, Electrical Conduction in Biological System, Robot Control, Electrochemical Analysis, Viscous Elastic Mechanics, Chaotic Dynamics and Molecular Spectroscopy, and so forth.Research on FO calculus has become one of the hot issues nowadays [1][2][3].
HS is a combination of multiphase medium, so its property can be affected by variation of air condition with the temperature, oil viscidity and compressibility, function time of external forces, and loading histories.At the same time, an excess of hypotheses and ignores and always focusing on variation of a narrow range of parameters or ensuring calculation accuracy by adding computational items make the traditional modeling methods difficult to describe the process precisely, quickly and globally.The property of HS that is related to variation of temperature, function time of external forces, and loading history perfectly accords with that of viscoelastic material.The viscoelastic materials have the mechanics characteristics of both elasticity and viscidity, which are between elastomer and Newtonian fluid.Its stress-strain response relies on time and strain rate and is related to external load and deformation history, which means that its stressstrain has memory.And viscoelastic model of FO calculus can describe the dynamic properties of many sophisticated viscoelastic materials in a wide range of frequency with only a few parameters [4][5][6].
In this paper, a HSM-FO is built for a certain multiaxle heavy vehicle, and the simulation results of IO model and FO model and the data of bench test are compared, which prove the HSM-FO feasible and demonstrate its advantage.

Structure Principle of HS and
Modeling Based on IO 2.1.Structure and Working Principle of Hydropneumatic Suspension.The one-wheel structure sketch of HS is shown in Figure 1, where 1 is vehicle body, 2 is single-action cylinder with oil in the upper chamber and connection to the air in the lower chamber, 3 is vehicle wheel, 4 is adjustable throttle, which serves as all-time through hole connected with other oil lines in parallel, 5 and 6 are compression and rebound relief valve, respectively, which can adjust the damping by changing the valve-opening pressure through spinning the tuning screw, and 7 is diaphragm accumulator with nitrogen of certain pressure in the upper chamber and connection to the oil line in the lower chamber.
With the input of road excitation, relative motion will come about between the wheel and body of vehicle while driving.When the piston goes upward, the HPS is in compression stroke, oil will be pushed towards the accumulator, and there will be two oil lines in parallel: throttle 4 and relief valve 5.These two lines work under different pressures (due to the piston speed).If the piston moves at a relatively low velocity, the system pressure is low, the relief valve remains shut, and oil can only flow into the accumulator through the throttle.If the piston moves at a relatively high velocity, the system pressure goes up till the opening pressure of relief valve 5.Then, oil will flow into the accumulator through a combination of oil line in parallel composed of the throttle line and the relief valve line.
When the piston goes downward, the HS is in rebound stroke and the oil in lower chamber of the accumulator will be pushed back to the tank through two parallel lines: throttle 4 and relief valve 6.These two lines work under different pressures (due to the piston speed).The working principle is the same as the compression stroke.
According to the analysis above, working states of each valve in compression and rebound strokes are shown in Table 1.

Components Models of HS.
The main components models including the accumulator, throttle, and relief valve models are first built with the structure principle analysis of the system.The elastic force model of single cylinder, single oil line, and single accumulator is as follows: These equations can be analyzed to develop the accumulator pressure at any time: where   is spring mass of single wheel;  is the lever rate;   is diameter of piston;  0 is precharged pressure of accumulator;  0 is original volume of accumulator;   is gas pressure at static balance;   is gas volume at static balance;  is gas pressure at random time;  is gas volume at random time;   is displacement of piston;  is gas polytropic index.The elastic force of the system can be obtained by converting the pressure variation of gas in the accumulator to the cylinder piston as shown in the following: As it can be seen in ( 3), the elastic force is a function on displacement of cylinder piston with all parameters of the system set.
The throttle is like an all-time through hole of thin wall, as shown in Figure 2. The variation of gravitational potential energy of oil while flowing in the cylinder is ignored and the cross-sectional area of pipe line is much larger than the alltime through hole of thin wall.
With the simplification above, the relationship between flow and pressure of all-time through hole is developed based on Bernoulli equation: where   is flow passing through the all-time through hole;   is flow coefficient;   is diameter of all-time through hole; Δ is pressure difference of two ends of the all-time through hole;  is density of oil.The model of relief valve can be built as a cone valve with a preloaded spring, as shown in Figure 3.Under ideal conditions, the radial distribution of liquid flow rate and pressure is symmetrical due to the symmetrical structure of the valve spool; thus, the external force on the valve spool can just be analyzed in the radial.According to the difference of the two ends of the valve port, the relief has three states: shut, open, and open to the maximum.
Thus, the mathematical model of flow rate and pressure difference of relief valve can be expressed as follows: where ℎ  = ((1/4) 2  Δ −   ℎ 0− )/  ;   is flow passing through the relief valve;   is flow coefficient;   is diameter of relief valve; Δ is pressure difference of two ends of relief valve;  is density of oil;  is stiffness of preloaded spring; ℎ 0 is the preload value of spring; ℎ  is opening value of relief valve; ℎ max is the maximum opening value of relief valve.

HS Model Based on Integral
Order.On the basis of components models of the system, the system model of HS can be built according to the connection relationship of the system.
In the compression stroke, if the vibrating velocity of the suspension is relatively low, the differential pressure produced by oil flowing through throttle 4 cannot open relief valve 5. Then the oil can only flow into the accumulator through throttle 4, and the volume of oil flowing through throttle 4 equals that of expellant oil due to piston movement of the cylinder, which is shown in the following: The above equations can be simplified to obtain the differential pressure of valves' ends in the following: The damping force is mainly from pressure drop produced by oil flowing through valves; thus, the damping force by oil flowing through valves is such that where   is the stiffness of preloaded spring of relief valve in compressor line; ℎ 0− is the original precompressed spring length of relief valve in compressor line;   is diameter of the hole of relief valve in compressor line.
In the compression stroke, if the vibrating velocity of the suspension is relatively high, the differential pressure produced by oil flowing through throttle 4 can open relief valve 5.The oil can flow into the accumulator through both throttle 4 and relief valve 5, and the volume of oil flowing through throttle 4 and relief valve equals that of expellant oil due to piston movement of the cylinder, as shown in the following: It can be simplified to obtain the relationship of differential pressure and the piston velocity, such that There is no explicit solution in (10) and the expression of differential pressure cannot be obtained, yet the differential pressure can be achieved in Matlab/Simulink, and the product of differential and the piston area is the damping force.Then the damping force is as follows: In the compression stroke, if the vibrating velocity of the suspension keeps increasing and the opening of the relief valve reaches the maximum, oil can flow into the accumulator through both throttle 4 and relief valve 5, and the volume of oil flowing through throttle 4 and relief valve equals that of expellant oil due to piston movement of the cylinder, as shown in the following: where ℎ max − is the maximum opening of relief valve in compressor line.
It can be simplified to obtain the differential pressure of valves as in the following: The damping force produced by oil flowing through valves is The mathematical model of the damping force can be obtained by analyzing the equations above, such that The same situation is in the rebound stroke; then the mathematical model of damping force in the rebound stroke is as follows: where   is the stiffness of preloaded spring of relief valve in rebound line; ℎ 0− is the original precompressed spring length of relief valve in rebound line;   is diameter of the hole of relief valve in rebound line; ℎ max − is the maximum opening of relief valve in rebound line.
The sign of the piston velocity indicates whether it is in the compression or the rebound stroke.The fractional calculus of Riemann-Liouville is defined as follows:

HS Modeling Based on Fractional Order
where     is differential operator; ,  are the lower and upper limit of differential operator respectively;  is the order of differential;  is an arbitrary integer; Γ(•) is gamma function, when  is a positive integer, Γ() = ( − 1)!;  is integration variable.
If the initial condition is 0, the FO calculus of Riemann-Liouville can be expressed as follows: (18)

HS Modeling Based on Fractional Order.
As in the analysis above, HS has the property of viscoelastic materials and can be described by a viscoelastic system.Then the FO model is introduced into the system to describe the elastic force and the damping item, and HSM-FO is built to analyze the output characteristics of the HS system.
On the basis of HSM-FO, the piston displacement in the elastic force model and the piston velocity in the damping force model are replaced.
The elastic force model based on FO is expressed in the following: where if  = 0,     =   .The damping force model based on FO in the compression stroke is shown in the following: The damping force model based on FO in the rebound stroke is shown in the following: where if  = 1,     = ẋ  .

The Calculation Methods for Fractional Order
At present, common calculation methods for FO include analytical method, numerical method, and filtering approximation method.The problem with the analytical method is the complicated formula transformation, huge computation, and too high dimension of the state space, which makes it inconvenient to be solved and further analyzed and only suitable for special equations.Thus, the numerical method and the filtering approximation method are more suited to engineering application.
Research on the traditional IO system is relatively mature, and a feasible research idea and method is to approximate the FO system with an IO system; then the issue is converted to approximating the FO operator with the standard IO operator.Research has shown that the method of rational number approximation is more effective with satisfactory accuracy and can operate differentiation and integration to displacement signal.The method of rational function approximating FO proposed by Oustaloup et al. is the most representative [11][12][13][14].
Oustaloup filtering approximation method is obtained based on the research on complex FO differentiator.The transfer function is as shown in the following: where  is positive time constant;  +  is order;  is real part of order;  is imaginary part of order; () is system input; () is system output.The expression of () in complex frequency domain is obtained with Laplace transform to (22).The unity-gain frequency (or inversion frequency)   = 1/ makes sense only if  = 0. Consider Thus, the transfer function is The excessively high and low frequency parts without research value can be removed, the chosen working frequency range is [  ,  ℎ ], and the   in ( 24) is replaced with the upper and low limits of the working frequency range; then the following is obtained: where   = (   ℎ ) 1/2 and  0 =   /  =   / ℎ .Equation ( 25) is rewritten in the form Equation ( 26) can be expressed as the recursive form of distribution of real poles and zeros through direct recursion based on FO theory, which is shown in the following: where where −   is the th zero; −  is the th pole.There are 2 + 1 poles and zeros, which means the order is  = 2 + 1. Apparently, the approximation accuracy will be higher with the increase of the IO order of the transfer function, yet the approximation accuracy will not go up proportionally to the order after it reaches a certain value.Instead, the complexity will be increased and computing speed decreased, which makes it meaningless to overall performance of system.Research shows that if  = 5, whether for amplitude or phase frequency, excellent approximation to the real frequency characteristics of FO filter can be obtained with different order  [15].
For Oustaloup filtering approximation method, the demands of accuracy and calculated amount can be met when  = 5; thus, rational function of fifth order is applied in this paper to approximate the FO function.The approximation frequency range [0.01 rad/s, 100 rad/s] is chosen considering the effective range of vehicle working frequency, where   = 0.01 rad/s and  ℎ = 100 rad/s.The order of FO items,  and , to make the FO model closer to reality is not confirmed, so the trial and error method is applied to obtain the optimal values of  and .The values of  and  and the corresponding rational function of fifth order are shown in Table 2.

Simulation Analysis and Experimental Verification
5.1.Bench Test.In order to verify the validity of the mathematical model of HS elastic and damping force, a bench test on the HS system of single cylinder, single oil line, and single accumulator that loaded with throttle and relief valve is conducted in this paper.The testing equipment is the FCS testing system of vehicle suspension.The cylinder, valves, and the accumulator are connected to each other and the cylinder is installed on the testing system upright according to the system diagram, as shown in Figure 4.
The test conditions are listed in Table 3.
The loading spectrum is shown in Table 4.
Under the loading spectrum excitation of 0.01 Hz/30 mm, the damping force generated by oil flowing through the damper valve is almost zero; thus, the output force of the piston rod can be approximately seen as elastic force.The damping force of the system is obtained by removing the elastic force from the system output force under the loading spectrum excitation of 1.00 Hz/30 mm.

Simulation of FO Model and Experimental Verification.
The elastic and damping force model of the system with single cylinder, single oil line, and single accumulator are built in Matlab/Simulink; then the elastic force and damping force are simulated and analyzed by inputting different excitation signal of piston displacement.
The simulation parameters are shown in Table 5.
With the excitation of 0.01 Hz/30 mm, comparison between the simulation curve of the elastic force based on IO and FO and the experimental curve is shown in Figure 5.
Figure 5 is the comparison between simulation curves of the elastic force based on IO and FO and the experimental curve of the elastic force.The results of comparison show that the simulation curve of the elastic force is closer to the experimental curve when the order of fractional operator of the elastic force model  = −0.1.In Figure 5, the blue solid line is the experimental curve of the elastic force.The black dot line and the red dashed line are the simulation curves of the elastic force model based on IO and FO, respectively.It can be seen from the figure that the simulation curve closer to the experimental curve can be obtained by adjusting the order of fractional operator  in the elastic force model, which proves the feasibility and validity of the elastic force model based on FO and demonstrates the advantage of the order of fractional operator being arbitrarily and continuously adjustable.
With the excitation of 1.0 Hz/30 mm, comparison between the simulation curve of the damping force based on IO and FO and the experimental curve is shown in Figure 6.Utilizing the same method for the damping force model based on FO, the comparison shows that, under the excitation of 1.0 Hz/30 mm, the simulation curve of the damping force is closer to the experimental curve when the order of fractional operator of the damping force model  = 1.1.In Figure 6, the blue solid line is the experimental curve of the damping force.The black dot line and the red dashed line are the simulation curves of the damping force model based on IO and FO, respectively.It can be seen from the figure that the simulation curve closer to the experimental curve can be obtained by adjusting the order of fractional operator  in the elastic force model, which proves the feasibility and validity of the elastic force model based on FO and reflects the contribution of the order of fractional operator being arbitrarily and continuously adjustable to describing a physical process precisely [16,17].model and indicates the effect of the order of fractional operator on the amplitude and phase of the model.The following will be the research on the influence of the order of fractional operator on the velocity of the cylinder piston, the dampingtime curve, the damping-displacement curve, and dampingvelocity curve.The excitation is a sine wave of 1.0 Hz/30 mm.

The Influence of the
Displayed in Figure 7 is the velocity variation curve of cylinder piston with different fractional operator's order.In the figure, the black lateral dashed line is the line of zero velocity, the longitudinal dashed line is the time line, and the purple solid line is the velocity curve of cylinder piston in condition of  = 1.0.The other six curves are the velocity

Figure 1 :
Figure 1: The one-wheel structure sketch of HS.

Figure 4 :
Figure 4: The testing system of vehicle suspension and the experimental components.
Order of Fractional Operator on the Damping Force.The comparison of simulation of FO model and the experiments proves the feasibility and validity of FO 0.00 0.01 0.02 0.03 0

Figure 5 :
Figure 5: Comparison between simulation and the experimental curve of the elastic force.

Figure 6 :
Figure 6: Comparison between simulation and the experimental curve of the damping force.

Table 1 :
Working states of each valve in compression and rebound strokes.

Table 2 :
Values of  and  and the corresponding Oustaloup transfer function of fifth order.

Table 4 :
Loading spectrum of sine wave excitation.

Table 5 :
Settings of the simulation parameters.