Mathematical Modelling and Computational Simulation of the Hydraulic Damper during the Orifice-Working Stage for Railway Vehicles

,e objective of this paper is to establish an accurate nonlinear mathematical model of the hydraulic damper during the orificeworking stage. A new mathematical model including the submodels of the orifices, hydraulic fluids, pressure chambers, and reservoir chambers is established based on theories of the fluid mechanics, hydropneumatics, and mechanics. Subsequently, a force element based on the established model of the hydraulic damper which contains 56 inputs, 6 force states, and 47 outputs is developed with the FORTRAN language in the secondary development environment of the multibody dynamics software SIMPACK. Using the force element, the damping characteristics of the modified yaw damper with different diameters of the base orifice are calculated under different amplitudes and frequencies of the sine excitation, and then the simulation results are compared with the experimental results which are obtained under the same conditions. Results show that during the orificeworking stage, the new established mathematical model can accurately reproduce the nonlinear static and dynamic characteristics of hydraulic dampers such as the force-displacement characteristic, force-velocity characteristic, fluid shortage, hysteresis effect, and pressure limited effect. Furthermore, it also shows that the nonlinear characteristics of the orifice, air release, cavitation, leakage for high frequencies, and dynamic characteristics of fluid (i.e., the density, bulk modulus, and air/gas content) should be taken seriously during the modelling of the hydraulic damper at the orifice-working stage. ,e mathematical model proposed in this paper is more applicable to the railway vehicle system dynamics and individual system description of the hydraulic damper.


Introduction
For railway vehicles, multibody vehicle dynamics simulation has become a useful design instrument, which can predict vehicle performances and reduce the manpower and funding.Although this approach can help to understand the basic phenomena of the vehicle dynamics such as the ride quality, stability, and curving behaviour, a number of complex issues cannot be solved because the vehicles are directly built with some simple linear elements which are far from the representative of the real one [1].Besides, too simple models may introduce misleading results in terms of the wheel-rail force and stability.In order to obtain precise results when simulating the dynamic behaviour of a railway vehicle, it is necessary to develop accurate models of all the elements that play an important part in the dynamics of the vehicle (wheel-rail contact [2], coil springs [3][4][5], air springs [6][7][8][9][10][11], rubber springs [12,13], dampers, etc.).
One of the elements that most affects the railway vehicle dynamics is the hydraulic damper.Many practical engineering problems such as the hunting instability and abnormal vibration are all correlated to the performance of the hydraulic damper.
erefore, obtaining the optimal performance parameters of the hydraulic damper has been the permanent goals of vehicle designers.e optimal design of parameters can be good or bad depending mainly on the accuracy of the model of the hydraulic damper.In recent years, many researchers have dedicated to the model of the hydraulic damper and many models have been established.
e common models are the Maxwell model [14][15][16], improved Maxwell model [17], experimental model [18][19][20], and parametric model [21][22][23][24][25][26][27].For railway vehicles, the most commonly used model is the Maxwell model, which consists a linear spring in series with a linear dashpot, and this model has been implemented into the simulation software (SIMPACK, ADAMS/RAIL, etc.).e linear or nonlinear force-velocity characteristic which is obtained by putting the maximum damping forces versus the maximum velocities under different sinusoidal signals is an input for the Maxwell model, where two main problems exist.Firstly, the vehicle performance relates closely to the shape of the damping curve, not only to some absolute values of the damping force as a function of maximum velocities.Secondly, this damping force generated by the hydraulic damper is the function not only of the velocity, but also the displacement, temperature, frequency, and amplitude.As a whole, the Maxwell model can hardly satisfy the needing of the vehicle dynamics.e improved Maxwell model and experimental model are based on experimental results, which means that both models are not suitable for the design process.Whether for the parameter optimization of the hydraulic damper or the performance prediction of railway vehicles, the parametric model is the best choice because it creates a direct link between the physical properties of hydraulic dampers and the vehicle performances.A typical hydraulic damper for railway vehicles is shown in Figure 1, which will be described in detail in the Section 2. As can be seen, the orifices and valves are the mainly damping elements.e existing parametric models mostly focus on the modelling of different valves but weaken the role of the orifice.Actually, the idea of these models is putting the cart before the horse for railway vehicles.Figure 2(a) shows the relative displacement of both ends of a yaw damper when the EMU was accelerated from 80 km/h to 320 km/h on the Beijing-Shanghai line, and the corresponding frequency spectrum is shown in Figure 2(b).As shown, the relative displacement is mostly less than 2 mm and the frequency is mainly concentrated within the frequency of 15 Hz, which means that the relative velocity between the piston and pressure tube is not too high.erefore, the orifices mainly provide the damping force and the valves which are distributed in the piston unit and base unit generally do not work in service.For this reason, the valves inside the hydraulic damper are usually called relief valves for railway vehicles.
From the above analysis we can see that the performance of the hydraulic damper during the orifice-working stage should be a focal research for railway vehicles.
With the rapid development of railway in our country, the performance of hydraulic dampers is demanded higher.
e key challenges that have been brought out in the vehicle design process can be summarised as follows: (i) How can we ensure that the stiffness and damping of the hydraulic damper are reasonable during the early stage of the vehicle design process?
(ii) Can the input parameters of the hydraulic damper be obtained without experiment when the whole vehicle is simulated?
(iii) How to design internal parameters (orifice length, rod diameter, piston diameter, etc.) for damper manufacturers?(iv) How can an accurate model of the hydraulic damper be convenient for vehicle dynamics simulation? e first three challenges can be solved by developing an accurate hydraulic damper model, and the last challenge can be solved by developing a new damper force element in the multibody dynamics software.
erefore, our study focuses on the mathematical model of the hydraulic damper during the orifice-working stage.A new mathematical model of the hydraulic damper is established including the submodels of the orifice, fluid, and chambers.Subsequently, a new force element of the hydraulic damper with the FORTRAN language in the secondary development environment of the multibody dynamics software SIMPACK is developed based on the established nonlinear model.A modified yaw damper with different diameters of the base orifice is selected as the test object.e simulation results obtained from the above model are compared with the experimental results in order to determine the model's accuracy.

Hydraulic Damper Description
A typical hydraulic damper for railway vehicles is shown in Figure 1.
e devices generating the damping force are mainly the piston unit (Figure 1 e orifice, commonly arranged in the spool of relief valves, mainly provides the damping force during the hydraulic damper operation because the relief valves do not work normally; in other words, the hydraulic damper generates the damping force mainly due to the resistance of a fluid passing through the orifice, particularly during the low velocity zone.For the common hydraulic damper, there are three orifices: one is placed in one of the spool of the two rebound relief valves distributed in the piston unit, other is placed in one of the spool of the two compression relief valves distributed in the piston unit, and the third one is placed in one of the spool of the three compression relief valves distributed in the base valve unit.Two orifices distributed in the piston unit work all the time during the rebound and compression stroke.
In order to establish a mathematical model of the hydraulic damper during the orifice-working stage, a modified hydraulic damper is selected for research, the only modification of which is that the unloading force of all the relief valves inside the piston unit and base valve unit is adjusted to 25 KN in order to obtain as much experimental data as possible, removing the disturbance of relief valves and producing the damping force from orifices.
e fluid flow directions of the compression stroke and rebound stroke are shown in the dotted arrow and the solid arrow in Figure 1(a), respectively.Mathematical Problems in Engineering

Mathematical Modelling of the Hydraulic Damper
3.1.Orifice.e small orifice (not necessarily circular) which plays a key role during the service of railway vehicles is the restricting device of a hydraulic damper.e characteristics of the flow through the orifice, generally accurately obtained by the experiment, must be precise for the performance prediction of a hydraulic damper.For the limited number of the reliable experiment and the inevitable deviation between test and actual working conditions, the experimental method is considered as an inefficient way for the orifice modelling.Nevertheless, in some certain circumstances, the analytical method for the orifice's nonlinear modelling is feasible and desirable.
Figure 3 shows a basic flow through an orifice, from upstream to downstream.A vena contracta, the neck in the flow where the cross-sectional area of the flow is the minimum and the velocity is the maximum, is observed.e volumetric discharge of an incompressible noncavitating fluid through an orifice is given by the Bernoulli equation [28]: where Q is the volumetric flow rate of the fluid, C c is the contraction coefficient, C v is the velocity coefficient defined as the ratio of the mean velocity at the vena contracta over the ideal speed, A o is the cross-sectional area of the orifice, A u is the upstream cross-sectional area, P up is the upstream pressure, P c is the pressure on the contraction section, and ρ is the density of the fluid, treated as a constant.
As area A u is much larger than A o , we can ignore the effect of A u on the volumetric flow rate Q, given by where e discharge coefficient C d varies with the orifice length l, hydraulic diameter d h (equal to the diameter d for the circular orifice), aspect ratio l/d h , Reynolds number Re, and so on [29], which means that the coefficient is not a constant value.In order to obtain the discharge coefficient C d , we must measure the volumetric flow rate Q, upstream pressure P up , and pressure on the contraction section P c .Unfortunately, it is almost impossible to measure the P c , so equation (2) cannot be directly used for modelling the hydraulic damper.
As the pressure at the vena contracta is difficult to measure, we can replace the pressure P c with the downstream pressure P down , so does the discharge coefficient C d with a new coefficient C q called the flow coefficient.Ultimately, equation (2) can be written as follows [28,29]: As with the discharge coefficient, the flow coefficient C q is also a function of the orifice geometry and Reynolds number.Traditional orifice modelling regards the coefficient as a constant value, typically 0.62 or 0.81, which is obviously unreasonable.If we do take a constant value, we are forced to have the gradient which will be infinity at the origin if Q against P up − P down , which cannot be and it is a numerical  disaster if you try to implement it.Clearly, the flow is laminar for sufficiently small pressure drops which means that the coefficient C q is certainly not a constant.Before determining the flow coefficient C q , the Reynolds number must be determined because it has a great influence on the flow coefficient C q .e conventional definition of the Reynolds number Re is in terms of the mean velocity U and orifice hydraulic diameter d h and is given as follows [28]: where ] is the coefficient of the kinematic viscosity, μ is the coefficient of the dynamic viscosity, and X wet is wetted perimeter of the orifice.e formula above means that the Reynolds number is a function of the volumetric flow rate Q related to the flow coefficient C q .In other words, the Reynolds number is also a function of the flow coefficient C q , which is not suitable for damper modelling because we would prefer an explicit relationship between C q and Re.In order to solve this problem, another new dimensionless number λ called the flow number instead of the Reynolds number Re is introduced, which is defined as follows [30]: From the modelling point of view, λ contains the quantities we easily know, which means λ can be calculated without knowing the flow coefficient C q , which is more easier to obtain the measurement C q � C q (λ) than C q � C q (Re), so the flow number λ has many advantages.
Actually, the flow coefficient C q relates not only to the flow number λ but also the orifice geometry (i.e., diameter ratio d h /D u and length/diameter l/d h ) [30], and C q may conveniently be expressed by where f is the function of the unrestricted form and D u is the upstream inlet diameter (i.e., the diameter of a pressure tube).Lichtarowicz et al. [30] pointed that the effect of the diameter ratio d h /D u on the flow coefficient can be ignored when the value is less than 0.25.In the real hydraulic damper system, the diameter ratio d/D is far less than 0.25.erefore, the flow coefficient C q is only related to the flow number λ and geometric ratio l/d h .
e influence of the flow number λ on the flow coefficient C q is not the same for different l/d h .In our study, according to the different length/diameter l/d h , the orifice is divided into three categories: thin-bladed orifice (l/d h < 2), short orifice (2 ≤ l/d h ≤ 10), and long orifice (l/d h > 10).As the long orifice is seldom used in the hydraulic damper for railway vehicles, it is beyond the scope of this article.
For the thin-bladed orifice (i.e., l/d h < 2), the experimental data has been taken by numerous workers and the results indicate that the substantially approximate constant value of C q will be obtained when the flow number λ reaches a high value λ lt where the turbulent occurs [31]. is flow characteristic of the thin-bladed orifice is described in Figure 4.For the flow number λ between zero and λ lt , its value is not specified and an arbitrary curve is plotted in the figure .To make numerical integration easy to run, a smooth curve between zero and λ lt is expected.erefore, in our study, the hyperbolic tangent function with a good smoothness at the point of zero or λ lt is introduced to describe the flow characteristic for different flow numbers between zero and λ lt .According to this rule, this flow coefficient for the bladed orifice is written as follows: where C qub is the ultimate constant value of C q at high flow number, λ lt is the critical flow number where laminar is converted to turbulent, and C b is a constant.In order to ensure the flow coefficient to be at least 95% of the ultimate value C qub when λ equals to the critical number λ lt , the following equation is obtained: Solving equation ( 8), we can obtain C b ≥ 1.8318, and the value of 1.9 is selected.
Next, we will concentrate on the ultimate constant value C qub .It is not difficult to understand that this value is obtained when the fluid flow is turbulent.Determining this value by the experiment for different l/d h was already performed by numerous workers [28,[30][31][32].Test results show that C qub rises from 0.61 to about 0.78 for l/d h between zero and unity, then continues to rise to a value of 0.81 at l/d h � 2. In view of this, in our modelling process, a piecewise linear function is used to determine C qub for different l/d h as follows: As for the critical flow number λ lt , the most accurate method used to determine λ lt is the experiment.Fortunately, as for the thin-bladed orifice, the value of λ lt is almost constant around 1700 in actual hydraulic dampers [29].erefore, the critical number can be regarded as the constant during the modelling of the thin-bladed orifice.
Up to now, all the parameters used to determine the flow coefficient C q for the thin-bladed orifice are identified explicitly.e variation of the flow coefficient C q with the flow number λ for different l/d h is shown in Figure 5(a).
For the short orifice, it has the characteristic of laminar or turbulent and the switch is made at a flow number by the ratio length/diameter l/d h .e flow coefficient is expressed Mathematical Problems in Engineering as the function of the flow number λ and the ratio length/ diameter l/d h .An experimental fitting formula describing C q with the flow number λ above 10 and the ratio length/diameter l/d h is given by [30] where C qus is the limiting value of the flow coefficient for the short orifice.
Referring the study in [30,31], the limiting value C qus as the function of the ratio length/diameter l/d h is expressed as follows: Equation ( 10) is only applicable for the value of the flow number above 10 [30].e flow number less than 10 means that the flow through the short orifice is laminar.
is situation is possible for the low velocity in the hydraulic damper.Attempting to extend equation (10) to low flow number has been made by some researchers [28,29,33].Results show that the flow coefficient C q is directly proportional to the square root of the flow number λ at the flow number less than 10 and can be formulated as follows: where θ is a constant value related to the ratio length/diameter l/d h .
Inserting λ of 10 into equations ( 10) and ( 12), for the sake of continuity at the flow number of 10, the constant value θ can be obtained as follows: Figure 5(b) shows its evolution of the flow coefficient C q with the flow number λ for different values of the length/ diameter l/d h .

Fluid Properties.
In the practical hydraulic damper, the fluid is a mixture of the basic fluid, dissolved air/gas, air/gas bubbles, and sometimes vapor [34].In our study, it is worth noting that the term "fluid" means the homogeneous mixture of the liquid and air/gas.For the air/gas free fluid, the term "liquid" will be used.Fluid properties are very important input parameters for the modelling of the hydraulic damper.From the modelling point of view, we will concentrate on three fluid properties: density, bulk modulus, and viscosity.

Density and Bulk Modulus.
e density is the mass of a substance per unit volume, which is the function of the pressure and temperature.e fluid bulk modulus, the reciprocal of the compressibility, represents the resistance of a liquid to compress, leading to the damper hydraulic stiffness [34,35].During the modelling of the hydraulic damper, many researchers treat the density and bulk modulus as a constant value, which leads to the abnormal evolution of the pressure, especially when the dissolved air in the fluid appears in the form of bubbles [35].In our modelling, we assume that all transformations are isothermal, which means that the temperature of the fluid remains constant.e isothermal bulk modulus or, for simplicity, the bulk modulus of the hydraulic fluid is defined in terms of the instantaneous density of the fluid as follows [34]: where B fluid is the instantaneous bulk modulus of the fluid (with air/vapor bubbles), ρ fluid is the instantaneous density of the fluid, P is the instantaneous pressure, and T is the current temperature of the fluid.
As the temperature is a constant value, ρ fluid and B fluid are just the function of the pressure, and we can integrate equation ( 14) from the atmospheric pressure P atm to the current instantaneous pressure P yields ρ fluid (P) � ρ fluid P atm exp Once one of the density and bulk modulus is determined, the other can be calculated from equation (14) or (15).During the modelling, the values for the density and bulk modulus must be consistent with these formulas or there will be the loss of mass conservation, resulting in unexplained phenomena.
In the real hydraulic damper, the hydraulic fluid always contains some air/gas which is known to have a substantial effect on the compressibility of the fluid so that the bulk modulus value varies as well.Air/gas is known to exist in the hydraulic damper in two forms: entrained air/gas and dissolved air/gas.It is possible for air/gas to change from one to the other depending on the conditions to which the fluid is subjected.
e entrained air/gas is in the form of air/gas bubbles which are dispersed in the fluid.e existence of the entrained air/gas significantly reduces the bulk modulus of the fluid.e dissolved air/gas is also in the form of bubbles which are invisible but stored in the space between the bigger molecules of the liquid.Test results indicate that the dissolved air/gas has no effect on the bulk modulus of the fluid [35].
e air/gas content is described by the volumetric fraction at the atmospheric pressure P atm and the temperature 273.15 K, which is where X g0 is the volumetric fraction of air/gas at P atm and 273.15 K, V l0 is the volume of the liquid at P atm and 273.15 K, and V g0 is the volume of air/gas at P atm and 273.15 K.
If we consider a unit volume of the liquid (i.e., V l0 � V 0 � 1 m 3 ), the corresponding volume of air/gas is given by e volumetric fraction X g0 is often supplied by oil suppliers, and its value is basically between 0.05%-0.2%.
Two kinds of air phenomena should be pointed out: one is the air/gas release and the other is the cavitation.e air/ gas release starts when the pressure drops below a pressure called saturation pressure P sat , where all air/gas are dissolved, and finishes when the pressure reaches a pressure called vapor pressure P vap .e cavitation occurs when the pressure drops below the vapor pressure P vap [35,36].ese abovementioned phenomena are shown in Figure 6.Both P sat and P vap are important to the modelling of the hydraulic damper, which can be obtained from the technical document of oil supplied.
For P > P sat , all the small air/gas molecules hidden between the bigger molecules of the liquid, in other words, the dissolved air/gas do not increase the volume but increase the mass, which is where V fluid_P_T is the volume of the fluid at the current pressure P and current temperature T and V liq_P_T is the volume of liquid at the current pressure P and current temperature T.
According to equation (17), the mass of fluid m fluid is given by where m liq is the mass of the liquid, m gas is the mass of air/ gas, ρ liq P atm 273 , and ρ air P atm 273 are the density of the Mathematical Problems in Engineering liquid and air/gas at P atm and 273.15 K, respectively, which are obtained easily from the technical document of oil supplied.
According to the definition of the density, we obtain where ρ liq P T is the density of the liquid at the current pressure P and current temperature T.
Referring to equation ( 15), we obtain where B liq P T is the bulk modulus of the liquid at the current pressure P and current temperature T.
Because the working pressure of the hydraulic damper is less than 10 MPa, the effect of the pressure on the liquid bulk modulus is neglected [35].Besides, as abovementioned, all transformations are isothermal, and then the bulk modulus of liquid B liq P T is a constant value.
Assuming the density of the liquid is not dependent on the temperature [29], equation ( 21) can be written as Using equations ( 18)-( 20) and ( 22), the density of the fluid at P and T can be written as Using equations ( 14) and ( 23), the bulk modulus of the fluid at P and T is given by In equations ( 23) and ( 24), the bulk modulus of the liquid at the current pressure P and current temperature T B liq P T is determined by Hayward [38], declaring that the bulk modulus of any hydraulic oil can be predicted about 5% deviation with only knowing the kinematic viscosity of oil at the atmospheric pressure and 20 °C, and we obtain where ] P atm 20 deg is the kinematic viscosity of the liquid at P atm and 20 °C, which can be determined in equation ( 36), and T deg is the current operating temperature in °C.
For P vap < P < P sat , part of the air is dissolved and the rest is free (or entrained).According to Herry's law [35], the volume fraction of undissolved air/gas is given by where V air_entrained is the volume of the entrained air/gas at P and T, V air_dissolved is the volume of the dissolved air/gas at P and T, and V air_all is the volume of the total air/gas at P and T. Equation ( 17) is the total volume of air/gas at P atm and 273.15 K. e equation of the state for an ideal gas PV � mRT gives the total volume of air/gas from 273.15 K to T at the constant Pressure P atm : where V air_all_T is the volume of the total air/gas at P atm and T. e polytropic equation PV c � constant gives the volume from P atm to the current pressure P: T 273.15 where c is the polytropic index for the air/gas/vapor content.
Replacing equation (28) in equation ( 26), the volume of the entrained air/gas at P and T will be T 273.15 Applying equations ( 20) and ( 22), the volume of the liquid (with the dissolved air/gas) V liq_with_dissolved at P and T is given by According to the definition of the density, equations ( 29) and ( 30), the fluid density at P and T can be obtained as follows: Substituting equations ( 26) and ( 31) into ( 14), the bulk modulus of the fluid at P and T is given by ( For P < P vap , there are the vapor and entrained air/gas but no liquid.e dynamics of this process is complex and the physical is not fully understood.e primary objective is to ensure the bulk modulus of the fluid is reduced by the presence of air/gas and vapor. In order to maintain the continuity at P vap , substituting P � P vap into equation ( 31) yields Substituting equation (34) into equation ( 14), the bulk modulus of the fluid at P and T is given by

Mathematical Problems in Engineering
For P > P sat , all air is dissolved in the fluid and the absolute viscosity of the fluid can be regarded as a constant value of the liquid.Besides, the density of the liquid is independent to the temperature, so μ fluid P T � μ liq T P atm � ] liq T P atm • ρ liq P atm T � ] liq T P atm • ρ liq P atm 273 , (38) where μ fluid P T is the absolute viscosity of the fluid at a specified pressure P and T and μ liq T P atm is the absolute viscosity of the liquid at P atm and T.
For P vap < P < P sat , the fluid is a mixture of the liquid with a constant absolute viscosity μ liq T P atm and free air/gas with a constant absolute viscosity μ air (usually 0.02 cP).In our study, the fluid absolute viscosity is the mean value based on the volumes of the liquid and air/gas as follows: Substituting equations ( 29) and ( 30) into (39) yields For P < P vap , the fluid is a mixture of the air/gas and vapor.Because of the same magnitude of the air/gas viscosity and vapor viscosity, the fluid absolute viscosity can be the same as the air/gas absolute viscosity, which is Once the absolute viscosity is obtained, the kinematic viscosity at a specified pressure can be calculated by dividing the fluid density at the same pressure, which can be used for calculating the flow number and flow rate.

Rebound Chamber and Compression Chamber. For a fluid volume, the rate of change of mass is
where m is the mass, ρ is the instantaneous density, V is the instantaneous volume, t is the time, m in is the input mass, and m out is the output mass.
Assuming the instantaneous density is distributed uniformly, equation (42) can be written as follows: where Q in and Q out are the input volume flow rate and the output volume flow rate, respectively.Substituting equation ( 14) into equation (43) yields It should be noted that the flow rate of Q in and Q out are the flow rate at P.
For the convenience of the following description, the orifice placed in one of the two rebound relief valves is designed as orifice1, the other orifice distributed in the piston unit is designed as orifice2, and the orifice distributed in the base valve unit is designed as orifice3.
We assume that the direction of the flow from the rebound chamber to the compression chamber through the orifice1 is positive and the direction of the flow from the compression chamber to the rebound chamber through the orifice2 is positive.
e direction of the flow from the compression chamber to the reservoir chamber through the orifice3 is positive, and the direction of the flow from the reservoir chamber to the compression chamber through the check valve in the base valve unit is positive.
According to equation (44), the time derivatives of the pressure in the rebound chamber and compression chamber are where P reb and P com are, respectively, the pressure of rebound and compression chamber; q orifice1 , q orifice2 , and q orifice3 are, respectively, the flow rate through the orifice1, orifice2, and orifice3; q check is the flow rate through the check valve; A piston is the crossing area of the piston; A rod is the crossing area of the rod; x and x • are the position and velocity of the piston with the positive direction of the rebound stroke; and L reb0 and L com0 are, respectively, the initial length of the rebound and compression chamber, as shown in Figure 7(a).B fluid_P reb _T is the bulk modulus of the fluid at 10 Mathematical Problems in Engineering P reb and T and B fluid_P com _T is the bulk modulus of the fluid at P com and T. ese flow rates of q orifice1 , q orifice2 , q orifice3 , and q check can be calculated through equation ( 3).Significantly, the density and kinematic viscosity used for those above flow rates are evaluated at the mean pressure(P up + P down )/2.ese values of flow rates which cannot be used directly must be corrected to the corresponding pressure, for instance, the flow rate through the orifice1 for the compression chamber is 3.4.Reservoir Chamber.e compressibility of the fluid in the reservoir is insignificant compared with the air bag and hence the hydraulic fluid in the reservoir chamber is assumed to have the same pressure as the air bag, so does their pressure derivative.
We have been conditioned to think that the initial pressure of the air bag is the atmospheric pressure, which is actually wrong if you know the assembly process of the hydraulic damper.Figure 7(a) illustrates the end state of the assembly and operating state.When the hydraulic damper is assembled, the piston is at the state of the top of the pressure tube, at which the pressure of the air bag is equal to the atmospheric pressure, that is, these air bags are in the initial state, as shown in Figure 7(b).Subsequently, the hydraulic damper is compressed to the operating state where the piston is usually in the middle of the pressure tube.It is obvious that the pressure of the air bag is higher than the atmospheric pressure.
In the operating state, the initial total volume of the air bag V gas0 is where k is the number of the air bag, V k is the precharge volume of the air bag, L c0 is the active length of the pressure tube, and L p is the width of the piston, as shown in Figure 7(a).According to the isothermal law P 1 V 1 � P 2 V 2 , the initial pressure P gas0 is calculated as follows: e air in the reservoir is assumed to obey the polytropic gas law of the form P • V c � constant.Differentiating the polytropic equation, the air pressure derivative is calculated as follows: where P gas is the air instantaneous pressure and V gas the total air instantaneous volume.
Assuming the compressibility of the fluid is ignored, the air volume derivative is Using equations ( 47)-( 50), the states of the reservoir chamber can be described.

Computational Simulation of the Hydraulic Damper.
Once mathematical models of main elements are established, the lumped nonlinear model of the hydraulic damper can be constructed according to the law of the flow rate equilibrium.Using these mathematical equations established above, a new hydraulic damper force element is developed with the FORTRAN language in the secondary development environment of the multibody dynamics software SIMPACK.e integration method is the SODASRT2.e new hydraulic damper model contains 56 inputs with definite physical meanings, 6 force states, and 47 outputs.Figure 8 shows the flow chart of the characteristic calculation of the hydraulic damper.Main simulation parameters of the proposed model are described in Table 1.

Model Verification and Simulation Results
In order to determine the accuracy of the established mathematical model during the orifice-working stage, large amounts of experimental data have been gathered under different amplitudes and frequencies of sine excitation, using the self-designed test bench SPTR100 with adjustable angle, as shown in Figure 9(a).Besides, in order to eliminate the influence of the rubber joints every hydraulic damper has in their ends, a special designed hydraulic fixture shown in Figure 9(b) is used to clamp the tested damper.e tested hydraulic damper is a modification of a yaw damper used for China Railway High-speed, which is described in the above section of "hydraulic damper description".e geometry dimensions of the modified damper are shown in Table 1.
During the orifice-working stage, the orifice is the primary damping element inside the hydraulic damper, thus modelling the orifice is critical to the mathematical model.For this reason, orifices with different diameters (i.e., 0.4, 0.5, and 0.6 mm) shown in Figure 9(c) are installed into the base valve unit so as to obtain as much experimental data as possible during the orifice-working stage.

Static Characteristics.
e static characteristics of the hydraulic damper are based on large excitation amplitudes and low excitation frequencies.In order to eliminate these errors caused by the signal processing, experimental signals of the force and the displacement from corresponding sensors (i.e., force sensor and displacement sensor) without processing are compared with the simulation results.In addition, for the sake of description of hysteresis of the force-velocity diagram, a calculated differential velocity signal is used appropriately.For the purpose of the comparison, computer simulation results together with the experimental results for each case are depicted in the same figure.
New model Test Maxwell model - ±5 mm at 0.2, 0.3, 0.4, and 0.5 Hz when the diameter of the base orifice (i.e., orifice3) is 0.4 mm. Figure 11 shows these comparisons of the force-velocity diagram under the same conditions.As can be seen, in general, the experimental results are very close to the simulation results obtained by the established mathematical physical model.For low frequencies (for example, 0.2 and 0.3 Hz), the values of the damping force are about a few hundred Newton less than the experimental results, the reason of which is that the Coulomb friction which can be observed at maximum or minimum displacement experimentally in Figure 10(a) is not included in the established mathematical model.Figure 11 indicates that the damping force is not the simple single-value function of velocity which is considered traditionally, while in reality, the damping force developed by the hydraulic damper is a nonlinear function of velocity with the hysteresis characteristic, which means the hysteresis loop occurs in the forcevelocity diagram.Moreover, the characteristic of the damping force versus the velocity is frequency dependent; that is, the hysteresis effect becomes more obvious as the frequency increases.e deviation between the simulation results of the Maxwell model and experimental results indicate that the Maxwell model cannot accurately simulate the static characteristics of the hydraulic damper.
In order to fully validate the established mathematical model during the orifice-working stage, the diameter of the orifice3 of 0.4 mm is substituted to the 0.5 and 0.6 mm.Due to the space limitation, some of the comparison results of simulation and test for the 0.5 and 0.6 mm are merely shown in Figures 12 and 13 under the excitation 14 Mathematical Problems in Engineering amplitude of 5 mm and 10 mm, respectively.As can be seen, the simulation results obtained by the established the mathematical physical model are basically consistent with the experimental results.In Figure 12, the damping force increases very slowly at the beginning of the rebound stroke because the bulk modulus of the fluid is reduced by the presence of the undissolved air or vapor.In this case, although the obvious fluid shortage is not observed, the pressure rise is limited.As the diameter of the base orifice increases to 0.6 mm, as shown in Figure 12, the obvious phenomena of the fluid shortage appears in the rebound chamber so that the rebound damping force will not be formed until the air hollow layer fades away.e fluid shortage occurs because of the unreasonable configuration, orifice blocking, leakage, and so on.In our study, the large base orifice makes the fluid easier to flow through and a small amount of fluid is forced into the rebound chamber, resulting in the fluid shortage.e reappearance of the fluid shortage owes to the accurate modelling of the air release and cavitation.Negative absolute pressure, obviously unrealistic, will occur if the air release and cavitation are ignored.In particular, greater discrepancies are found between simulation results and actual measurements in Figure 12(c), the main reason of which is that the relief valves do work when the damping force is around 25 kN but not considered in the process of simulation.As our study focuses on the orificeworking stage, the effect of relief valves are beyond the scope of this article and there is no point in comparison.e simulation reappearance of the pressure limited effect and the fluid shortage indicates that the fluid modelling is accurate.
As a whole, the established mathematical model can effectively reproduce the static characteristics of the hydraulic damper, for instance, the fluid shortage, hysteresis effect, and pressure limited effect, but the Maxwell cannot.

Dynamic Characteristics.
e dynamic characteristics of the hydraulic damper are based on small excitation amplitudes and high excitation frequencies.In our study, the experimental amplitudes are selected as 0.5 and 1 mm with the frequency ranged from 2 to 15 Hz.Similarly, these direct signals without processing from different sensors are used for comparison.Figures 14 and 15 show the comparisons of the force-displacement diagram at different frequencies for the amplitudes of 0.5 and 1 mm, respectively, when the diameter of the base orifice is 0.4 mm.As the relief valves do work when the excitation frequency reaches above 6 Hz for the amplitude of 1 mm, the comparisons above 6 Hz are not carried out.As can be seen, the simulation results of the new model are basically consistent with the experimental results.Of course, some inevitable deviations exist (for example, 0.5 mm at 15 Hz) because the coulomb friction and leakage are not

16
Mathematical Problems in Engineering considered in the established mathematical model of the hydraulic damper.Besides, the higher the excitation frequency is, the greater the deviation of simulation and experiment will be, which can be observed for both 0.5 and 1 mm.So, in other words, the leakage should be considered in the mathematical model for high frequencies.However, the results of the Maxwell model show that it cannot reproduce the dynamic characteristics especially for the small amplitude.Like the static characteristic, different diameters of the base orifice are used for the purpose of the adequate verification of the established mathematical model.Figure 16 shows the comparisons of the force-displacement diagram for the diameter of 0.5 mm, and Figure 17 is used for the diameter of 0.6 mm.In addition, different amplitudes of 0.5 and 1 mm are used in order to obtain as much experimental data as possible.All of these results show the validity and accuracy of the established mathematical model.As can be seen, a more obvious pressure limited effect (for example, the dashed region in Figure 16(a) where the pressure rises slowly) is observed during the rebound stroke when the diameter of the base orifice is 0.5 mm by comparing Figures 14 and 16. e main reason for this is that the larger base orifice makes the fluid easier to flow through resulting in the air release and cavitation in the rebound chamber.It is not difficult to infer that the pressure limited effect for the diameter of 0.6 mm is more obvious, indeed confirmed in Figure 17.Moreover, the point where the pressure limited effect starts moves to a higher displacement as the frequency increases, which is shown in Figure 17 as an example.
In generally, the established mathematical model of the hydraulic damper can simulate the nonlinear dynamic characteristics efficiently during the orifice-working stage; however, all these characteristics cannot be described by the Maxwell model.

Simulation Results of the Orifice and Fluid.
As an example, the flow characteristics of the orifice under the amplitude of 5 mm at 4 Hz are shown in Figure 18.In our study, the orifice1 is the short orifice and the orifice3 is the thin-bladed orifice.As can be seen, the flow coefficients of both orifices change from about 0 to 0.7 rather than a constant value.Besides, for the orifice1, the time of flow coefficient changes from about 0 to 0.35 is quite short, indicating that the nonlinearity of the hydraulic damper is significant strong when the piston changes the direction.In summary, the nonlinear flow characteristics of the short orifice and thin-bladed orifice should be taken seriously rather than a simple constant value of the flow coefficient for the modelling of the hydraulic damper during the orificeworking stage.
As stated previously, fluid properties are crucial to the modelling of the hydraulic damper.Similarly, as an example, Figure 19 illustrates the fluid properties of the rebound  Mathematical Problems in Engineering chamber under the amplitude of 10 mm at 0.4 Hz when the diameter of the base orifice is 0.6 mm.As can be seen, the air/gas volumetric content increases resulting from the air release and cavitation.As the value of the air/gas volumetric content reaches to about 1%, the fluid bulk modulus is near to zero, which we do not expect but often face with, especially for the KONI damper which has no air bags in the reservoir chamber.e fluid bulk modulus is greatly influenced by the air/gas volumetric content, and for our current case, the value of the bulk modulus changes from 0 to 11000 bar, which means that the fluid bulk modulus cannot be considered as a constant value.Moreover, the air/gas content also has a great effect on the density related to the flow rate of the orifice.
e density decreases rapidly with the increasing air/gas content in the fluid.In our preliminary research process, the bulk modulus and density are considered as approximate constant values, which means the air/gas release and cavitation are not introduced in the fluid, and as a result, the numerical integration is terminated in the process of the simulation.Fortunately, the absolute viscosity of the fluid changes relatively little and can be regarded as a constant in the modelling of the hydraulic damper.In a word, the characteristics of the fluid are not steady but dynamic, so the instantaneous density, bulk modulus, and air/gas content should be calculated seriously, but the absolute viscosity can be regarded as a constant value under the isothermal condition.Mathematical Problems in Engineering

Conclusions
In conclusion, below is a list of the main results obtained from this work: (1) A new mathematical model of the hydraulic damper during the orifice-working stage has been proposed.e lumped model, including the submodels of the orifice, fluid, and chambers, has been implemented into a new force element of the hydraulic damper so as to be used for the performance analysis of hydraulic dampers and full vehicle simulation, developed by the FORTRAN language in the multibody dynamics software SIMPACK.e density, bulk modulus, and air/gas content of the fluid cannot be regarded as the constant value but the absolute viscosity be under the isothermal condition.

Data Availability
All data included in this study are available upon request by contacting the corresponding author.
(b)) and base valve unit (Figure1(c)).Each of the relief valves has a preloaded coil spring and will not open until the pressure differential across the valve reaches the specified value.

Figure 1 :
Figure 1: Configuration of a typical hydraulic damper: (a) two-dimensional assembly graph, (b) piston unit, and (c) base valve unit.

Figure 2 :
Figure 2: (a) Relative displacement and (b) the corresponding frequency spectrum of a yaw damper in service.

Figure 3 :
Figure 3: Orifice schematic diagram of the flow characteristic.

Figure 4 :
Figure 4: Variation of the flow coefficient for the thin-bladed orifice.

FlowFigure 5 :
Figure 5: Variation of the flow coefficient with the flow number for different l/d h .(a) in-bladed orifice.(b) Short orifice.

Figure 7 :
Figure 7: e states of the hydraulic damper and configuration of the air bag.(a) Comparison of the end state of assembly and the operating state.(b) Air bag.

Figure 8 :
Figure 8: Flow chart of the characteristic calculation of the hydraulic damper.

Figure 9 :
Figure 9: Illustration of the test apparatuses and the orifice.(a) Test bench of SPTR100.(b) Hydraulic fixture.(c) Test orifices.

( 2 )
Large numbers of experimental results of different damper configurations prove that the above mathematical model can accurately reproduce the nonlinear static characteristics including the fluid shortage, hysteresis effect, and pressure limited effect.(3) Compared with test results, it indicates that the model can simulate the nonlinear dynamic characteristics during the orifice-working stage.(4) e Maxwell model cannot accurately reproduce the static and dynamic characteristics.(5) e modelling of leakages is suggested to be taken into consideration under the high frequency.(6) e short orifice and thin-bladed orifice both have obvious nonlinear characteristics, and the nonlinear flow characteristics of them should be taken seriously rather than simple constant value of the flow coefficient.(7) e phenomena of the air release and cavitation are crucial to the modelling of the hydraulic damper during the orifice-working stage.If ignored, unrealistic results or the integrator interruption will occur.(8) e characteristics of the fluid are not steady but dynamic.
the viscosity at T 1 and T 2 in °C and at atmospheric pressure ] liq T 1 P atm and ] liq T 2 P atm in cSt are provided, the viscosity at a specified temperature T in With c 2 � log 10 log 10 ] liq T 1 P atm + 0.7     − log 10 log 10 ] liq T 2 P atm + 0.7     log 10 T 2 + 273.15  − log 10 T 1 + 273.15  , c 1 � log 10 log 10 ] liq T 1 P atm + 0.7     + c 2 • log 10 T 1 + 273.15 .

Table 1 :
Main simulation parameters of the proposed model.