Abstract
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.
1. 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–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 stress-strain 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–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.
2. 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.
2.2. 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; is precharged pressure of accumulator; 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 all-time 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 ; 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; is the preload value of spring; is opening value of relief valve; is the maximum opening value of relief valve.
2.3. 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; 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 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; is the original precompressed spring length of relief valve in rebound line; is diameter of the hole of relief valve in rebound line; 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.
3. HS Modeling Based on Fractional Order
3.1. The Definition of Fractional Order
The three common definition modes of fractional order are the following: Riemann-Liouville, Grunwald-Letnikov, and Caputo. The most applied one in engineering is Riemann-Liouville [7–10].
The fractional calculus of Riemann-Liouville is defined as follows: 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, ; is integration variable.
If the initial condition is 0, the FO calculus of Riemann-Liouville can be expressed as follows:
3.2. 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 , .
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 , .
4. 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–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) makes sense only if . 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 and .
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 poles and zeros, which means the order is . 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 , 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 ; 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 rad/s and 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.
5. 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.
5.2. 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 . 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 . 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].
5.3. The Influence of the 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 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 damping-time curve, the damping-displacement curve, and damping-velocity curve. The excitation is a sine wave of 1.0 Hz/30 mm.
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 . The other six curves are the velocity curves of cylinder piston with , 0.8, 0.9, 1.1, 1.2, 1.3, respectively.
As can be seen in the figure, when , the piston velocity is 0 at 0.75 s, reaches the maximum velocity at 1.00 s, and returns to 0 at 1.25 s with the amplitude of 0.188 m/s. The piston velocity under is synchronous in phase and has the same amplitude as the sine excitation signal 1.0 Hz/30 mm. When , the piston velocity has delays in phase and decrease in amplitude compared to the sine excitation signal, and the delays will be more and the amplitude smaller with the decrease of the order (as the curves of piston velocity when , 0.8, 0.9 shown in Figure 7). When , the piston velocity is relatively advanced to the sine excitation signal in phase and has larger amplitude; also, the phase advance will be more and the amplitude larger with the increase of order (as the curves of piston velocity when , 1.2, 1.3).
Figure 8 indicates the damping force-time curves under different orders of FO, and the similar conclusions can be drawn. When , the piston velocity and the sine excitation signal are synchronous in phase and the amplitude is in the middle of the 7 curves. When , delays of the damping force in phase compared to the sine excitation signal will be more and the amplitude smaller with the decrease of the order. When , the phase advance of the damping force will be more and the amplitude larger with the increase of order. The wave crests and troughs in the figure are different in amplitude and waveforms due to the difference of the parameters settings of the relief valve in the compression and rebound strokes.
Figure 9 shows the damping force-displacement curves under different orders of FO; the parts above zero are the damping force curves in the compression stroke and the others are that of the rebound stroke. It can be seen in the figure that when the left and the right parts of the damping force curves based on the zero line of displacement are symmetric. When , at the displacement of the same absolute value (e.g., the displacement of 0.02 m and −0.02 m), the left half curves of the compression stroke are smaller in amplitude than the right half and the left half curves of the rebound stroke are larger in amplitude than the right half. When , at the displacement of the same absolute value, the left half curves of the compression stroke are larger in amplitude than the right half and the left half curves of the rebound stroke are smaller in amplitude than the right half, which is just opposite to the situation of . And this asymmetry will be more distinct with the decrease or increase of the order. The overall variation of the amplitude in the 7 curves indicates that the amplitude of the damping force increases along with the order.
Figure 10 displays the damping force-velocity curves under different orders of FO; the parts above zero are the damping force curves in the compression stroke and the others are that of the rebound stroke. When , the damping force-velocity curves overlap as one curve in a complete compression and rebound stroke. When , there is not an overlap of curves but an envelope area in a complete compression and rebound stroke, and the envelope area becomes larger with the decrease of the order. When , no overlap occurs and the envelope area increases along with the order. The overall variation in the 7 curves indicates that the damping force increases along with the order, and the farther the order is from 1.0, the bigger the misalignment will be between the curves, that is, the larger the envelope area will be.
6. Conclusions
A HSM-FO is built with the advantage of FO in viscoelastic material modeling considering the mechanics property of multiphase medium of HS, and the detailed calculation method is given utilizing the Oustaloup filtering algorithm. Then the bench test is performed to demonstrate the validity of the model built.
The conclusions drawn from the simulation analysis of the model are as follows.(1)Comparison between simulation curves of the HS models based on IO and FO and the experimental curve of the HS proves the feasibility and validity of the HSM-FO and demonstrates the advantage of the order of fractional operator being arbitrarily and continuously adjustable.(2)Variation of the order in HSM-FO will give rise to delays or advance in phase and increase or decrease in amplitude of the piston velocity and damping force of the suspension.(3)Variation of the order in HSM-FO will cause no overlap of curves but an envelope area in a complete compression and rebound stroke, and the farther the order is from the integral order with actual physical meaning, the bigger the misalignment will be between the curves, that is, the larger the envelope area will be.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The project is supported by the National Natural Science Foundation of China (Grant no. 51375046 and Grant no. 51205021).