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

Authors: Hongxing Gao; Maoru Chi; Liangcheng Dai; Jungang Yang; Xiaozhi Zhou
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The objective of this paper is to establish an accurate nonlinear mathematical model of the hydraulic damper during the orifice-working 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 orifice-working 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. The mathematical model proposed in this paper is more applicable to the railway vehicle system dynamics and individual system description of the hydraulic damper.


---

## Body

## 1. 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–5], air springs [6–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. Therefore, obtaining the optimal performance parameters of the hydraulic damper has been the permanent goals of vehicle designers. The 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. The common models are the Maxwell model [14–16], improved Maxwell model [17], experimental model [18–20], and parametric model [21–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.). The 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. The 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. The 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. Therefore, 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.Figure 1
Configuration of a typical hydraulic damper: (a) two-dimensional assembly graph, (b) piston unit, and (c) base valve unit.
(a)
(b)
(c)Figure 2
(a) Relative displacement and (b) the corresponding frequency spectrum of a yaw damper in service.
(a)
(b)With the rapid development of railway in our country, the performance of hydraulic dampers is demanded higher. The 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?The 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. Therefore, 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. The simulation results obtained from the above model are compared with the experimental results in order to determine the model’s accuracy.


## 2. Hydraulic Damper Description

A typical hydraulic damper for railway vehicles is shown in Figure1. The devices generating the damping force are mainly the piston unit (Figure 1(b)) and base valve unit (Figure 1(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.The 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.The fluid flow directions of the compression stroke and rebound stroke are shown in the dotted arrow and the solid arrow in Figure1(a), respectively.


## 3. Mathematical Modelling of the Hydraulic Damper

### 3.1. Orifice

The small orifice (not necessarily circular) which plays a key role during the service of railway vehicles is the restricting device of a hydraulic damper. The 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.Figure3 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. The volumetric discharge of an incompressible noncavitating fluid through an orifice is given by the Bernoulli equation [28]:(1)Q=CcCvAo1−CcAo/Au22Pup−Pcρ,where Q is the volumetric flow rate of the fluid, Cc is the contraction coefficient, Cv is the velocity coefficient defined as the ratio of the mean velocity at the vena contracta over the ideal speed, Ao is the cross-sectional area of the orifice, Au is the upstream cross-sectional area, Pup is the upstream pressure, Pc is the pressure on the contraction section, and ρ is the density of the fluid, treated as a constant.Figure 3
Orifice schematic diagram of the flow characteristic.As areaAu is much larger than Ao, we can ignore the effect of Au on the volumetric flow rate Q, given by(2)Q=CdAo2Pup−Pcρ,where Cd is the discharge coefficient, Cd = CcCv.The discharge coefficientCd varies with the orifice length l, hydraulic diameter dh (equal to the diameter d for the circular orifice), aspect ratio l/dh, Reynolds number Re, and so on [29], which means that the coefficient is not a constant value. In order to obtain the discharge coefficient Cd, we must measure the volumetric flow rate Q, upstream pressure Pup, and pressure on the contraction section Pc. Unfortunately, it is almost impossible to measure the Pc, 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 pressurePc with the downstream pressure Pdown, so does the discharge coefficient Cd with a new coefficient Cq called the flow coefficient. Ultimately, equation (2) can be written as follows [28, 29]:(3)Q=CqAo2Pup−Pdownρ.As with the discharge coefficient, the flow coefficientCq 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 Pup − Pdown, 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 Cq is certainly not a constant.Before determining the flow coefficientCq, the Reynolds number must be determined because it has a great influence on the flow coefficient Cq. The conventional definition of the Reynolds number Re is in terms of the mean velocity U and orifice hydraulic diameter dh and is given as follows [28]:(4)Re=Udhν=Udhμ/ρ=Q/Aodhμ/ρ=ρQdhμAo,dh=4AoXwet,where ν is the coefficient of the kinematic viscosity, μ is the coefficient of the dynamic viscosity, and Xwet is wetted perimeter of the orifice.The formula above means that the Reynolds number is a function of the volumetric flow rateQ related to the flow coefficient Cq. In other words, the Reynolds number is also a function of the flow coefficient Cq, which is not suitable for damper modelling because we would prefer an explicit relationship between Cq 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]:(5)λ=dhν2Pup−Pdownρ.From the modelling point of view,λ contains the quantities we easily know, which means λ can be calculated without knowing the flow coefficient Cq, which is more easier to obtain the measurement Cq=Cqλ than Cq=CqRe, so the flow number λ has many advantages.Actually, the flow coefficientCq relates not only to the flow number λ but also the orifice geometry (i.e., diameter ratio dh/Du and length/diameter l/dh) [30], and Cq may conveniently be expressed by(6)Cq=fλ,ldh,dhDu,where f is the function of the unrestricted form and Du 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 dh/Du 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. Therefore, the flow coefficient Cq is only related to the flow number λ and geometric ratio l/dh.The influence of the flow numberλ on the flow coefficient Cq is not the same for different l/dh. In our study, according to the different length/diameter l/dh, the orifice is divided into three categories: thin-bladed orifice (l/dh < 2), short orifice (2 ≤ l/dh ≤ 10), and long orifice (l/dh > 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/dh < 2), the experimental data has been taken by numerous workers and the results indicate that the substantially approximate constant value of Cq will be obtained when the flow number λ reaches a high value λlt where the turbulent occurs [31]. This 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.Figure 4
Variation of the flow coefficient for the thin-bladed orifice.To make numerical integration easy to run, a smooth curve between zero andλlt is expected. Therefore, 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:(7)Cq=CqubtanhCb⋅λλlt,where Cqub is the ultimate constant value of Cq at high flow number, λlt is the critical flow number where laminar is converted to turbulent, and Cb is a constant.In order to ensure the flow coefficient to be at least 95% of the ultimate valueCqub when λ equals to the critical number λlt, the following equation is obtained:(8)Cq=CqubtanhCb×1.0≥0.95×Cqub.Solving equation (8), we can obtain Cb≥1.8318, and the value of 1.9 is selected.Next, we will concentrate on the ultimate constant valueCqub. 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/dh was already performed by numerous workers [28, 30–32]. Test results show that Cqub rises from 0.61 to about 0.78 for l/dh between zero and unity, then continues to rise to a value of 0.81 at l/dh = 2. In view of this, in our modelling process, a piecewise linear function is used to determine Cqub for different l/dh as follows:(9)Cqub=0.61+0.17ldh,if0<ldh≤1,0.03ldh+0.75,if1<ldh<2.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]. Therefore, 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 coefficientCq for the thin-bladed orifice are identified explicitly. The variation of the flow coefficient Cq with the flow number λ for different l/dh is shown in Figure 5(a).Figure 5
Variation of the flow coefficient with the flow number for different l/dh. (a) Thin-bladed orifice. (b) Short orifice.
(a)
(b)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/diameterl/dh. The flow coefficient is expressed as the function of the flow number λ and the ratio length/diameter l/dh. An experimental fitting formula describing Cq with the flow number λ above 10 and the ratio length/diameter l/dh is given by [30](10)Cq=11/Cqus+20/λ1+2.25l/dh−0.005l/dh/1+7.5log100.00015λ2,where Cqus is the limiting value of the flow coefficient for the short orifice.Referring the study in [30, 31], the limiting value Cqus as the function of the ratio length/diameter l/dh is expressed as follows:(11)Cqus=0.827−0.0085ldh.Equation (10) is only applicable for the value of the flow number above 10 [30]. The flow number less than 10 means that the flow through the short orifice is laminar. This 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 Cq is directly proportional to the square root of the flow number λ at the flow number less than 10 and can be formulated as follows:(12)Cq=θλ,where θ is a constant value related to the ratio length/diameter l/dh.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:(13)θ=0.3162Cqus1+2+4.5l/dhCqus.Figure5(b) shows its evolution of the flow coefficient Cq with the flow number λ for different values of the length/diameter l/dh.


### 3.2. 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.


#### 3.2.1. Density and Bulk Modulus

The density is the mass of a substance per unit volume, which is the function of the pressure and temperature. The 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. The 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]:(14)Bfluid=ρfluid∂P∂ρfluidT,where Bfluid 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 Bfluid are just the function of the pressure, and we can integrate equation (14) from the atmospheric pressure Patm to the current instantaneous pressure P yields(15)ρfluidP=ρfluidPatmexp∫PatmPdPBfluidP.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. The entrained air/gas is in the form of air/gas bubbles which are dispersed in the fluid. The existence of the entrained air/gas significantly reduces the bulk modulus of the fluid. The 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].The air/gas content is described by the volumetric fraction at the atmospheric pressurePatm and the temperature 273.15 K, which is(16)Xg0=Vg0Vg0+Vl0,where Xg0 is the volumetric fraction of air/gas at Patm and 273.15 K, Vl0 is the volume of the liquid at Patm and 273.15 K, and Vg0 is the volume of air/gas at Patm and 273.15 K.If we consider a unit volume of the liquid (i.e.,Vl0 = V0 = 1 m3), the corresponding volume of air/gas is given by(17)Vg0=Xg01−Xg0V0.The volumetric fractionXg0 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. The air/gas release starts when the pressure drops below a pressure called saturation pressurePsat, where all air/gas are dissolved, and finishes when the pressure reaches a pressure called vapor pressure Pvap. The cavitation occurs when the pressure drops below the vapor pressure Pvap [35, 36]. These abovementioned phenomena are shown in Figure 6. Both Psat and Pvap are important to the modelling of the hydraulic damper, which can be obtained from the technical document of oil supplied.Figure 6
Air/gas release and cavitation process.ForP > Psat, 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(18)Vfluid_P_T=Vliq_P_T,where Vfluid_P_T is the volume of the fluid at the current pressure P and current temperature T and Vliq_P_T is the volume of liquid at the current pressure P and current temperature T.According to equation (17), the mass of fluid mfluid is given by(19)mfluid=mliq+mgas=V0ρliq_Patm_273+Xg01−Xg0⋅ρair_Patm_273,where mliq is the mass of the liquid, mgas is the mass of air/gas, ρliq_Patm_273, and ρair_Patm_273 are the density of the liquid and air/gas at Patm 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(20)Vliq_P_T=mliqρliq_P_T=V0ρliq_Patm_273ρliq_P_T,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(21)ρliq_P_T=ρliq_Patm_Texp∫PatmPdPBliq_P_T,where Bliq_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 Bliq_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(22)ρliq_P_T=ρliq_Patm_273expP−PatmBliq_P_T.Using equations (18)–(20) and (22), the density of the fluid at P and T can be written as(23)ρfluid_P_T=mfluidVfluid_P_T=V0ρliq_Patm_273+Xg01−Xg0⋅ρair_Patm_273⋅expP−PatmBliq_P_T.Using equations (14) and (23), the bulk modulus of the fluid at P and T is given by(24)Bfluid_P_T=ρfluid_P_T1dρfluid_P_T/dP=Bliq_P_T=constant.In equations (23) and (24), the bulk modulus of the liquid at the current pressure P and current temperature T Bliq_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(25)Bliq_P_T=1.3+0.15log10ν_Patm_20deg104+20−Tdeg/435,where ν_Patm_20deg is the kinematic viscosity of the liquid at Patm and 20°C, which can be determined in equation (36), and Tdeg is the current operating temperature in °C.ForPvap < P < Psat, 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(26)ϕ=Vair_entrainedVair_all=Vair_entrainedVair_dissolved+Vair_entrained=1−P−PvapPsat−Pvap,where Vair_entrained is the volume of the entrained air/gas at P and T, Vair_dissolved is the volume of the dissolved air/gas at P and T, and Vair_all is the volume of the total air/gas at P and T.Equation (17) is the total volume of air/gas at Patm and 273.15 K. The 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 Patm:(27)Vair_all_T=V0Xg01−Xg0T273.15,where Vair_all_T is the volume of the total air/gas at Patm and T.The polytropic equationPVγ=constant gives the volume from Patm to the current pressure P:(28)Vair_all=V0Xg01−Xg0T273.15PatmP1/γ,where γ 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(29)Vair_entrained=ϕV0Xg01−Xg0T273.15PatmP1/γ.Applying equations (20) and (22), the volume of the liquid (with the dissolved air/gas) Vliq_with_dissolved at P and T is given by(30)Vliq_with_dissolved=Vliq_P_T=exp−P−PatmBliq_P_T.According to the definition of the density, equations (29) and (30), the fluid density at P and T can be obtained as follows:(31)ρfluid_P_T=mfluidVair_entrained+Vliq_with_dissolved=ρliq_Patm_273+Xg0/1−Xg0⋅ρair_Patm_273ϕXg0/1−Xg0T/273.15Patm/P1/γ+exp−P−Patm/Bliq_P_T.Substituting equations (26) and (31) into (14), the bulk modulus of the fluid at P and T is given by(32)Bfluid_P_T=ρfluid_P_T1dρfluid_P_T/dP=ϕXg0/1−Xg0T/273.15Patm/P1/γ+exp−P−Patm/Bliq_P_TXg0/1−Xg0T/273.15Patm/P1/γϕ/γP+1/Psat−Pvap+1/Bliq_P_Texp−P−Patm/Bliq_P_T.ForP < Pvap, there are the vapor and entrained air/gas but no liquid. The dynamics of this process is complex and the physical is not fully understood. The 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 atPvap, substituting P = Pvap into equation (31) yields(33)ρfluid_Pvap_T=ρliq_Patm_273+Xg0/1−Xg0⋅ρair_Patm_273ϕXg0/1−Xg0T/273.15Patm/Pvap1/γ+exp−Pvap−Patm/Bliq_P_T,where ρfluid_Pvap_T is the density of the fluid at Pvap and T.The polytropic equation change fromPvap to the current pressure P gives(34)ρfluid_P_T=ρliq_Patm_273+Xg0/1−Xg0⋅ρair_Patm_273ϕXg0/1−Xg0T/273.15Patm/Pvap1/γ+exp−Pvap−Patm/Bliq_P_T⋅P/Pvap1/γ.Substituting equation (34) into equation (14), the bulk modulus of the fluid at P and T is given by(35)Bfluid_P_T=ρfluid_P_T1dρfluid_P_T/dP=γP.


#### 3.2.2. Viscosity

Viscosity is a measure of the resistance of the fluid to flow. This characteristic has both positive and negative effects on the hydraulic damper. A low viscosity leads to leakages between the hood and spool, resulting in an unstable work diagram, but a high one will lead to the orifice blocked. There are two kinds of viscosity: dynamic viscosityμand kinematic viscosity ν. In general, suppliers always provide the pure liquid (without air/gas) kinematic viscosity at different temperatures (−20, 40, 100°C, etc.) under the same atmospheric pressure Patm. If the kinematic viscosity of a specified temperature is not provided, it can be calculated from these provided values. Or more accurately, if the viscosity at T1 and T2 in °C and at atmospheric pressure νliq_T1_Patm and νliq_T2_Patm in cSt are provided, the viscosity at a specified temperature T in °C and atmospheric pressure νliq_T_Patm in cSt is calculated as follows [38]:(36)νliq_T_Patm=1010c1⋅T+273.15c2−0.7.With(37)c2=log10log10νliq_T1_Patm+0.7−log10log10νliq_T2_Patm+0.7log10T2+273.15−log10T1+273.15,c1=log10log10νliq_T1_Patm+0.7+c2⋅log10T1+273.15.ForP > Psat, 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(38)μfluid_P_T=μliq_T_Patm=νliq_T_Patm⋅ρliq_Patm_T=νliq_T_Patm⋅ρliq_Patm_273,where μfluid_P_T is the absolute viscosity of the fluid at a specified pressure P and T and μliq_T_Patm is the absolute viscosity of the liquid at Patm and T.ForPvap < P < Psat, the fluid is a mixture of the liquid with a constant absolute viscosity μliq_T_Patm 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:(39)μfluid_P_T=Vair_entrainedμair+Vliq_with_dissolvedμliq_T_PatmVair_entrained+Vliq_with_dissolved.Substituting equations (29) and (30) into (39) yields(40)μfluid_P_T=ϕXg0/1−Xg0T/273.15Patm/P1/γμair+exp−P−Patm/Bliq_P_Tμliq_T_PatmϕXg0/1−Xg0T/273.15Patm/P1/γ+exp−P−Patm/Bliq_P_T.ForP < Pvap, 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(41)μfluid_P_T=μair.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.


### 3.3. Rebound Chamber and Compression Chamber

For a fluid volume, the rate of change of mass is(42)dmdt=dρVdt=ρdVdt+Vdρdt=∑min•−∑mout•,where m is the mass, ρ is the instantaneous density, V is the instantaneous volume, t is the time, min is the input mass, and mout is the output mass.Assuming the instantaneous density is distributed uniformly, equation (42) can be written as follows:(43)∑Qin−∑Qout−dVdt=VρdP/dρ⋅dPdt,where Qin and Qout are the input volume flow rate and the output volume flow rate, respectively.Substituting equation (14) into equation (43) yields(44)dPdt=Bfluid⋅∑Qin−∑Qout−dV/dtV.It should be noted that the flow rate ofQin and Qout 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. The 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(45)dPrebdt=Bfluid_Preb_T⋅qorifice2Preb−qorifice1Preb+Apiston−Arodx•Apiston−Arod⋅Lreb0−x,dPcomdt=Bfluid_Pcom_T⋅qorifice1Pcom−qorifice2Pcom−qorifice3Pcom+qcheckPcom−Apistonx•Apiston⋅Lcom0+x,where Preb and Pcom are, respectively, the pressure of rebound and compression chamber; qorifice1, qorifice2, and qorifice3 are, respectively, the flow rate through the orifice1, orifice2, and orifice3; qcheck is the flow rate through the check valve; Apiston is the crossing area of the piston; Arod 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 Lreb0 and Lcom0 are, respectively, the initial length of the rebound and compression chamber, as shown in Figure 7(a). Bfluid_Preb_T is the bulk modulus of the fluid at Preb and T and Bfluid_Pcom_T is the bulk modulus of the fluid at Pcom and T.Figure 7
The 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.
(a)
(b)These flow rates ofqorifice1, qorifice2, qorifice3, and qcheck can be calculated through equation (3). Significantly, the density and kinematic viscosity used for those above flow rates are evaluated at the mean pressurePup+Pdown/2. These 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(46)qorifice1Pcom=qorifice1⋅ρPup+Pdown/2ρPcom.


### 3.4. Reservoir Chamber

The 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. Figure7(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 bagVgas0 is(47)Vgas0=∑i=1kVk−π8Apiston2−Arod2⋅Lc0−Lp,where k is the number of the air bag, Vk is the precharge volume of the air bag, Lc0 is the active length of the pressure tube, and Lp is the width of the piston, as shown in Figure 7(a).According to the isothermal lawP1V1 = P2V2, the initial pressure Pgas0 is calculated as follows:(48)Pgas0=Patm⋅∑1kVkVgas0.The air in the reservoir is assumed to obey the polytropic gas law of the formP⋅Vγ=constant. Differentiating the polytropic equation, the air pressure derivative is calculated as follows:(49)dPgasdt=γ⋅PgasVgas⋅dVgasdt,where Pgas is the air instantaneous pressure and Vgas the total air instantaneous volume.Assuming the compressibility of the fluid is ignored, the air volume derivative is(50)dVgasdt=qcheckPgas−qorifice3Pgas.Using equations (47)–(50), the states of the reservoir chamber can be described.


### 3.5. 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. The integration method is the SODASRT2. The new hydraulic damper model contains 56 inputs with definite physical meanings, 6 force states, and 47 outputs. Figure8 shows the flow chart of the characteristic calculation of the hydraulic damper. Main simulation parameters of the proposed model are described in Table 1.Figure 8
Flow chart of the characteristic calculation of the hydraulic damper.Table 1
Main simulation parameters of the proposed model.
Parameters
Values
Parameters
Values
Diameter of orifice1
0.5 mm
Length of orifice1
1.75 mm
Diameter of orifice2
0.5 mm
Length of orifice2
1.75 mm
Diameter of orifice3
0.4, 0.5, 0.6 mm
Length of orifice3
0.75 mm
Diameter of pressure tube
60 mm
Diameter of piston
59.98 mm
Diameter of rod
30 mm
Width of piston
35 mm
X
g0
0.1 %
P
sat
75 bar
P
vap
−0.6 bar
ρ
liq
_
P
atm
_
273
854 kg/m3
ρ
air
_
P
atm
_
273
1.2 kg/m3
B
liq
_
P
_
T
15281 bar
T
25°C
γ
1.4
T
1
−20°C
ν
liq
_
T
1
_
P
atm
1430 cSt
T
2
100°C
ν
liq
_
T
2
_
P
atm
6.1 cSt
μ
air
0.02 cP
k
4
L
reb0
168 mm
L
com0
169 mm
L
c0
372 mm
V
k
110 ml


## 3.1. Orifice

The small orifice (not necessarily circular) which plays a key role during the service of railway vehicles is the restricting device of a hydraulic damper. The 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.Figure3 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. The volumetric discharge of an incompressible noncavitating fluid through an orifice is given by the Bernoulli equation [28]:(1)Q=CcCvAo1−CcAo/Au22Pup−Pcρ,where Q is the volumetric flow rate of the fluid, Cc is the contraction coefficient, Cv is the velocity coefficient defined as the ratio of the mean velocity at the vena contracta over the ideal speed, Ao is the cross-sectional area of the orifice, Au is the upstream cross-sectional area, Pup is the upstream pressure, Pc is the pressure on the contraction section, and ρ is the density of the fluid, treated as a constant.Figure 3
Orifice schematic diagram of the flow characteristic.As areaAu is much larger than Ao, we can ignore the effect of Au on the volumetric flow rate Q, given by(2)Q=CdAo2Pup−Pcρ,where Cd is the discharge coefficient, Cd = CcCv.The discharge coefficientCd varies with the orifice length l, hydraulic diameter dh (equal to the diameter d for the circular orifice), aspect ratio l/dh, Reynolds number Re, and so on [29], which means that the coefficient is not a constant value. In order to obtain the discharge coefficient Cd, we must measure the volumetric flow rate Q, upstream pressure Pup, and pressure on the contraction section Pc. Unfortunately, it is almost impossible to measure the Pc, 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 pressurePc with the downstream pressure Pdown, so does the discharge coefficient Cd with a new coefficient Cq called the flow coefficient. Ultimately, equation (2) can be written as follows [28, 29]:(3)Q=CqAo2Pup−Pdownρ.As with the discharge coefficient, the flow coefficientCq 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 Pup − Pdown, 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 Cq is certainly not a constant.Before determining the flow coefficientCq, the Reynolds number must be determined because it has a great influence on the flow coefficient Cq. The conventional definition of the Reynolds number Re is in terms of the mean velocity U and orifice hydraulic diameter dh and is given as follows [28]:(4)Re=Udhν=Udhμ/ρ=Q/Aodhμ/ρ=ρQdhμAo,dh=4AoXwet,where ν is the coefficient of the kinematic viscosity, μ is the coefficient of the dynamic viscosity, and Xwet is wetted perimeter of the orifice.The formula above means that the Reynolds number is a function of the volumetric flow rateQ related to the flow coefficient Cq. In other words, the Reynolds number is also a function of the flow coefficient Cq, which is not suitable for damper modelling because we would prefer an explicit relationship between Cq 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]:(5)λ=dhν2Pup−Pdownρ.From the modelling point of view,λ contains the quantities we easily know, which means λ can be calculated without knowing the flow coefficient Cq, which is more easier to obtain the measurement Cq=Cqλ than Cq=CqRe, so the flow number λ has many advantages.Actually, the flow coefficientCq relates not only to the flow number λ but also the orifice geometry (i.e., diameter ratio dh/Du and length/diameter l/dh) [30], and Cq may conveniently be expressed by(6)Cq=fλ,ldh,dhDu,where f is the function of the unrestricted form and Du 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 dh/Du 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. Therefore, the flow coefficient Cq is only related to the flow number λ and geometric ratio l/dh.The influence of the flow numberλ on the flow coefficient Cq is not the same for different l/dh. In our study, according to the different length/diameter l/dh, the orifice is divided into three categories: thin-bladed orifice (l/dh < 2), short orifice (2 ≤ l/dh ≤ 10), and long orifice (l/dh > 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/dh < 2), the experimental data has been taken by numerous workers and the results indicate that the substantially approximate constant value of Cq will be obtained when the flow number λ reaches a high value λlt where the turbulent occurs [31]. This 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.Figure 4
Variation of the flow coefficient for the thin-bladed orifice.To make numerical integration easy to run, a smooth curve between zero andλlt is expected. Therefore, 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:(7)Cq=CqubtanhCb⋅λλlt,where Cqub is the ultimate constant value of Cq at high flow number, λlt is the critical flow number where laminar is converted to turbulent, and Cb is a constant.In order to ensure the flow coefficient to be at least 95% of the ultimate valueCqub when λ equals to the critical number λlt, the following equation is obtained:(8)Cq=CqubtanhCb×1.0≥0.95×Cqub.Solving equation (8), we can obtain Cb≥1.8318, and the value of 1.9 is selected.Next, we will concentrate on the ultimate constant valueCqub. 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/dh was already performed by numerous workers [28, 30–32]. Test results show that Cqub rises from 0.61 to about 0.78 for l/dh between zero and unity, then continues to rise to a value of 0.81 at l/dh = 2. In view of this, in our modelling process, a piecewise linear function is used to determine Cqub for different l/dh as follows:(9)Cqub=0.61+0.17ldh,if0<ldh≤1,0.03ldh+0.75,if1<ldh<2.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]. Therefore, 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 coefficientCq for the thin-bladed orifice are identified explicitly. The variation of the flow coefficient Cq with the flow number λ for different l/dh is shown in Figure 5(a).Figure 5
Variation of the flow coefficient with the flow number for different l/dh. (a) Thin-bladed orifice. (b) Short orifice.
(a)
(b)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/diameterl/dh. The flow coefficient is expressed as the function of the flow number λ and the ratio length/diameter l/dh. An experimental fitting formula describing Cq with the flow number λ above 10 and the ratio length/diameter l/dh is given by [30](10)Cq=11/Cqus+20/λ1+2.25l/dh−0.005l/dh/1+7.5log100.00015λ2,where Cqus is the limiting value of the flow coefficient for the short orifice.Referring the study in [30, 31], the limiting value Cqus as the function of the ratio length/diameter l/dh is expressed as follows:(11)Cqus=0.827−0.0085ldh.Equation (10) is only applicable for the value of the flow number above 10 [30]. The flow number less than 10 means that the flow through the short orifice is laminar. This 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 Cq is directly proportional to the square root of the flow number λ at the flow number less than 10 and can be formulated as follows:(12)Cq=θλ,where θ is a constant value related to the ratio length/diameter l/dh.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:(13)θ=0.3162Cqus1+2+4.5l/dhCqus.Figure5(b) shows its evolution of the flow coefficient Cq with the flow number λ for different values of the length/diameter l/dh.


## 3.2. 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.


### 3.2.1. Density and Bulk Modulus

The density is the mass of a substance per unit volume, which is the function of the pressure and temperature. The 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. The 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]:(14)Bfluid=ρfluid∂P∂ρfluidT,where Bfluid 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 Bfluid are just the function of the pressure, and we can integrate equation (14) from the atmospheric pressure Patm to the current instantaneous pressure P yields(15)ρfluidP=ρfluidPatmexp∫PatmPdPBfluidP.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. The entrained air/gas is in the form of air/gas bubbles which are dispersed in the fluid. The existence of the entrained air/gas significantly reduces the bulk modulus of the fluid. The 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].The air/gas content is described by the volumetric fraction at the atmospheric pressurePatm and the temperature 273.15 K, which is(16)Xg0=Vg0Vg0+Vl0,where Xg0 is the volumetric fraction of air/gas at Patm and 273.15 K, Vl0 is the volume of the liquid at Patm and 273.15 K, and Vg0 is the volume of air/gas at Patm and 273.15 K.If we consider a unit volume of the liquid (i.e.,Vl0 = V0 = 1 m3), the corresponding volume of air/gas is given by(17)Vg0=Xg01−Xg0V0.The volumetric fractionXg0 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. The air/gas release starts when the pressure drops below a pressure called saturation pressurePsat, where all air/gas are dissolved, and finishes when the pressure reaches a pressure called vapor pressure Pvap. The cavitation occurs when the pressure drops below the vapor pressure Pvap [35, 36]. These abovementioned phenomena are shown in Figure 6. Both Psat and Pvap are important to the modelling of the hydraulic damper, which can be obtained from the technical document of oil supplied.Figure 6
Air/gas release and cavitation process.ForP > Psat, 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(18)Vfluid_P_T=Vliq_P_T,where Vfluid_P_T is the volume of the fluid at the current pressure P and current temperature T and Vliq_P_T is the volume of liquid at the current pressure P and current temperature T.According to equation (17), the mass of fluid mfluid is given by(19)mfluid=mliq+mgas=V0ρliq_Patm_273+Xg01−Xg0⋅ρair_Patm_273,where mliq is the mass of the liquid, mgas is the mass of air/gas, ρliq_Patm_273, and ρair_Patm_273 are the density of the liquid and air/gas at Patm 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(20)Vliq_P_T=mliqρliq_P_T=V0ρliq_Patm_273ρliq_P_T,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(21)ρliq_P_T=ρliq_Patm_Texp∫PatmPdPBliq_P_T,where Bliq_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 Bliq_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(22)ρliq_P_T=ρliq_Patm_273expP−PatmBliq_P_T.Using equations (18)–(20) and (22), the density of the fluid at P and T can be written as(23)ρfluid_P_T=mfluidVfluid_P_T=V0ρliq_Patm_273+Xg01−Xg0⋅ρair_Patm_273⋅expP−PatmBliq_P_T.Using equations (14) and (23), the bulk modulus of the fluid at P and T is given by(24)Bfluid_P_T=ρfluid_P_T1dρfluid_P_T/dP=Bliq_P_T=constant.In equations (23) and (24), the bulk modulus of the liquid at the current pressure P and current temperature T Bliq_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(25)Bliq_P_T=1.3+0.15log10ν_Patm_20deg104+20−Tdeg/435,where ν_Patm_20deg is the kinematic viscosity of the liquid at Patm and 20°C, which can be determined in equation (36), and Tdeg is the current operating temperature in °C.ForPvap < P < Psat, 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(26)ϕ=Vair_entrainedVair_all=Vair_entrainedVair_dissolved+Vair_entrained=1−P−PvapPsat−Pvap,where Vair_entrained is the volume of the entrained air/gas at P and T, Vair_dissolved is the volume of the dissolved air/gas at P and T, and Vair_all is the volume of the total air/gas at P and T.Equation (17) is the total volume of air/gas at Patm and 273.15 K. The 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 Patm:(27)Vair_all_T=V0Xg01−Xg0T273.15,where Vair_all_T is the volume of the total air/gas at Patm and T.The polytropic equationPVγ=constant gives the volume from Patm to the current pressure P:(28)Vair_all=V0Xg01−Xg0T273.15PatmP1/γ,where γ 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(29)Vair_entrained=ϕV0Xg01−Xg0T273.15PatmP1/γ.Applying equations (20) and (22), the volume of the liquid (with the dissolved air/gas) Vliq_with_dissolved at P and T is given by(30)Vliq_with_dissolved=Vliq_P_T=exp−P−PatmBliq_P_T.According to the definition of the density, equations (29) and (30), the fluid density at P and T can be obtained as follows:(31)ρfluid_P_T=mfluidVair_entrained+Vliq_with_dissolved=ρliq_Patm_273+Xg0/1−Xg0⋅ρair_Patm_273ϕXg0/1−Xg0T/273.15Patm/P1/γ+exp−P−Patm/Bliq_P_T.Substituting equations (26) and (31) into (14), the bulk modulus of the fluid at P and T is given by(32)Bfluid_P_T=ρfluid_P_T1dρfluid_P_T/dP=ϕXg0/1−Xg0T/273.15Patm/P1/γ+exp−P−Patm/Bliq_P_TXg0/1−Xg0T/273.15Patm/P1/γϕ/γP+1/Psat−Pvap+1/Bliq_P_Texp−P−Patm/Bliq_P_T.ForP < Pvap, there are the vapor and entrained air/gas but no liquid. The dynamics of this process is complex and the physical is not fully understood. The 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 atPvap, substituting P = Pvap into equation (31) yields(33)ρfluid_Pvap_T=ρliq_Patm_273+Xg0/1−Xg0⋅ρair_Patm_273ϕXg0/1−Xg0T/273.15Patm/Pvap1/γ+exp−Pvap−Patm/Bliq_P_T,where ρfluid_Pvap_T is the density of the fluid at Pvap and T.The polytropic equation change fromPvap to the current pressure P gives(34)ρfluid_P_T=ρliq_Patm_273+Xg0/1−Xg0⋅ρair_Patm_273ϕXg0/1−Xg0T/273.15Patm/Pvap1/γ+exp−Pvap−Patm/Bliq_P_T⋅P/Pvap1/γ.Substituting equation (34) into equation (14), the bulk modulus of the fluid at P and T is given by(35)Bfluid_P_T=ρfluid_P_T1dρfluid_P_T/dP=γP.


### 3.2.2. Viscosity

Viscosity is a measure of the resistance of the fluid to flow. This characteristic has both positive and negative effects on the hydraulic damper. A low viscosity leads to leakages between the hood and spool, resulting in an unstable work diagram, but a high one will lead to the orifice blocked. There are two kinds of viscosity: dynamic viscosityμand kinematic viscosity ν. In general, suppliers always provide the pure liquid (without air/gas) kinematic viscosity at different temperatures (−20, 40, 100°C, etc.) under the same atmospheric pressure Patm. If the kinematic viscosity of a specified temperature is not provided, it can be calculated from these provided values. Or more accurately, if the viscosity at T1 and T2 in °C and at atmospheric pressure νliq_T1_Patm and νliq_T2_Patm in cSt are provided, the viscosity at a specified temperature T in °C and atmospheric pressure νliq_T_Patm in cSt is calculated as follows [38]:(36)νliq_T_Patm=1010c1⋅T+273.15c2−0.7.With(37)c2=log10log10νliq_T1_Patm+0.7−log10log10νliq_T2_Patm+0.7log10T2+273.15−log10T1+273.15,c1=log10log10νliq_T1_Patm+0.7+c2⋅log10T1+273.15.ForP > Psat, 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(38)μfluid_P_T=μliq_T_Patm=νliq_T_Patm⋅ρliq_Patm_T=νliq_T_Patm⋅ρliq_Patm_273,where μfluid_P_T is the absolute viscosity of the fluid at a specified pressure P and T and μliq_T_Patm is the absolute viscosity of the liquid at Patm and T.ForPvap < P < Psat, the fluid is a mixture of the liquid with a constant absolute viscosity μliq_T_Patm 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:(39)μfluid_P_T=Vair_entrainedμair+Vliq_with_dissolvedμliq_T_PatmVair_entrained+Vliq_with_dissolved.Substituting equations (29) and (30) into (39) yields(40)μfluid_P_T=ϕXg0/1−Xg0T/273.15Patm/P1/γμair+exp−P−Patm/Bliq_P_Tμliq_T_PatmϕXg0/1−Xg0T/273.15Patm/P1/γ+exp−P−Patm/Bliq_P_T.ForP < Pvap, 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(41)μfluid_P_T=μair.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.


## 3.2.1. Density and Bulk Modulus

The density is the mass of a substance per unit volume, which is the function of the pressure and temperature. The 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. The 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]:(14)Bfluid=ρfluid∂P∂ρfluidT,where Bfluid 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 Bfluid are just the function of the pressure, and we can integrate equation (14) from the atmospheric pressure Patm to the current instantaneous pressure P yields(15)ρfluidP=ρfluidPatmexp∫PatmPdPBfluidP.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. The entrained air/gas is in the form of air/gas bubbles which are dispersed in the fluid. The existence of the entrained air/gas significantly reduces the bulk modulus of the fluid. The 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].The air/gas content is described by the volumetric fraction at the atmospheric pressurePatm and the temperature 273.15 K, which is(16)Xg0=Vg0Vg0+Vl0,where Xg0 is the volumetric fraction of air/gas at Patm and 273.15 K, Vl0 is the volume of the liquid at Patm and 273.15 K, and Vg0 is the volume of air/gas at Patm and 273.15 K.If we consider a unit volume of the liquid (i.e.,Vl0 = V0 = 1 m3), the corresponding volume of air/gas is given by(17)Vg0=Xg01−Xg0V0.The volumetric fractionXg0 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. The air/gas release starts when the pressure drops below a pressure called saturation pressurePsat, where all air/gas are dissolved, and finishes when the pressure reaches a pressure called vapor pressure Pvap. The cavitation occurs when the pressure drops below the vapor pressure Pvap [35, 36]. These abovementioned phenomena are shown in Figure 6. Both Psat and Pvap are important to the modelling of the hydraulic damper, which can be obtained from the technical document of oil supplied.Figure 6
Air/gas release and cavitation process.ForP > Psat, 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(18)Vfluid_P_T=Vliq_P_T,where Vfluid_P_T is the volume of the fluid at the current pressure P and current temperature T and Vliq_P_T is the volume of liquid at the current pressure P and current temperature T.According to equation (17), the mass of fluid mfluid is given by(19)mfluid=mliq+mgas=V0ρliq_Patm_273+Xg01−Xg0⋅ρair_Patm_273,where mliq is the mass of the liquid, mgas is the mass of air/gas, ρliq_Patm_273, and ρair_Patm_273 are the density of the liquid and air/gas at Patm 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(20)Vliq_P_T=mliqρliq_P_T=V0ρliq_Patm_273ρliq_P_T,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(21)ρliq_P_T=ρliq_Patm_Texp∫PatmPdPBliq_P_T,where Bliq_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 Bliq_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(22)ρliq_P_T=ρliq_Patm_273expP−PatmBliq_P_T.Using equations (18)–(20) and (22), the density of the fluid at P and T can be written as(23)ρfluid_P_T=mfluidVfluid_P_T=V0ρliq_Patm_273+Xg01−Xg0⋅ρair_Patm_273⋅expP−PatmBliq_P_T.Using equations (14) and (23), the bulk modulus of the fluid at P and T is given by(24)Bfluid_P_T=ρfluid_P_T1dρfluid_P_T/dP=Bliq_P_T=constant.In equations (23) and (24), the bulk modulus of the liquid at the current pressure P and current temperature T Bliq_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(25)Bliq_P_T=1.3+0.15log10ν_Patm_20deg104+20−Tdeg/435,where ν_Patm_20deg is the kinematic viscosity of the liquid at Patm and 20°C, which can be determined in equation (36), and Tdeg is the current operating temperature in °C.ForPvap < P < Psat, 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(26)ϕ=Vair_entrainedVair_all=Vair_entrainedVair_dissolved+Vair_entrained=1−P−PvapPsat−Pvap,where Vair_entrained is the volume of the entrained air/gas at P and T, Vair_dissolved is the volume of the dissolved air/gas at P and T, and Vair_all is the volume of the total air/gas at P and T.Equation (17) is the total volume of air/gas at Patm and 273.15 K. The 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 Patm:(27)Vair_all_T=V0Xg01−Xg0T273.15,where Vair_all_T is the volume of the total air/gas at Patm and T.The polytropic equationPVγ=constant gives the volume from Patm to the current pressure P:(28)Vair_all=V0Xg01−Xg0T273.15PatmP1/γ,where γ 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(29)Vair_entrained=ϕV0Xg01−Xg0T273.15PatmP1/γ.Applying equations (20) and (22), the volume of the liquid (with the dissolved air/gas) Vliq_with_dissolved at P and T is given by(30)Vliq_with_dissolved=Vliq_P_T=exp−P−PatmBliq_P_T.According to the definition of the density, equations (29) and (30), the fluid density at P and T can be obtained as follows:(31)ρfluid_P_T=mfluidVair_entrained+Vliq_with_dissolved=ρliq_Patm_273+Xg0/1−Xg0⋅ρair_Patm_273ϕXg0/1−Xg0T/273.15Patm/P1/γ+exp−P−Patm/Bliq_P_T.Substituting equations (26) and (31) into (14), the bulk modulus of the fluid at P and T is given by(32)Bfluid_P_T=ρfluid_P_T1dρfluid_P_T/dP=ϕXg0/1−Xg0T/273.15Patm/P1/γ+exp−P−Patm/Bliq_P_TXg0/1−Xg0T/273.15Patm/P1/γϕ/γP+1/Psat−Pvap+1/Bliq_P_Texp−P−Patm/Bliq_P_T.ForP < Pvap, there are the vapor and entrained air/gas but no liquid. The dynamics of this process is complex and the physical is not fully understood. The 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 atPvap, substituting P = Pvap into equation (31) yields(33)ρfluid_Pvap_T=ρliq_Patm_273+Xg0/1−Xg0⋅ρair_Patm_273ϕXg0/1−Xg0T/273.15Patm/Pvap1/γ+exp−Pvap−Patm/Bliq_P_T,where ρfluid_Pvap_T is the density of the fluid at Pvap and T.The polytropic equation change fromPvap to the current pressure P gives(34)ρfluid_P_T=ρliq_Patm_273+Xg0/1−Xg0⋅ρair_Patm_273ϕXg0/1−Xg0T/273.15Patm/Pvap1/γ+exp−Pvap−Patm/Bliq_P_T⋅P/Pvap1/γ.Substituting equation (34) into equation (14), the bulk modulus of the fluid at P and T is given by(35)Bfluid_P_T=ρfluid_P_T1dρfluid_P_T/dP=γP.


## 3.2.2. Viscosity

Viscosity is a measure of the resistance of the fluid to flow. This characteristic has both positive and negative effects on the hydraulic damper. A low viscosity leads to leakages between the hood and spool, resulting in an unstable work diagram, but a high one will lead to the orifice blocked. There are two kinds of viscosity: dynamic viscosityμand kinematic viscosity ν. In general, suppliers always provide the pure liquid (without air/gas) kinematic viscosity at different temperatures (−20, 40, 100°C, etc.) under the same atmospheric pressure Patm. If the kinematic viscosity of a specified temperature is not provided, it can be calculated from these provided values. Or more accurately, if the viscosity at T1 and T2 in °C and at atmospheric pressure νliq_T1_Patm and νliq_T2_Patm in cSt are provided, the viscosity at a specified temperature T in °C and atmospheric pressure νliq_T_Patm in cSt is calculated as follows [38]:(36)νliq_T_Patm=1010c1⋅T+273.15c2−0.7.With(37)c2=log10log10νliq_T1_Patm+0.7−log10log10νliq_T2_Patm+0.7log10T2+273.15−log10T1+273.15,c1=log10log10νliq_T1_Patm+0.7+c2⋅log10T1+273.15.ForP > Psat, 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(38)μfluid_P_T=μliq_T_Patm=νliq_T_Patm⋅ρliq_Patm_T=νliq_T_Patm⋅ρliq_Patm_273,where μfluid_P_T is the absolute viscosity of the fluid at a specified pressure P and T and μliq_T_Patm is the absolute viscosity of the liquid at Patm and T.ForPvap < P < Psat, the fluid is a mixture of the liquid with a constant absolute viscosity μliq_T_Patm 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:(39)μfluid_P_T=Vair_entrainedμair+Vliq_with_dissolvedμliq_T_PatmVair_entrained+Vliq_with_dissolved.Substituting equations (29) and (30) into (39) yields(40)μfluid_P_T=ϕXg0/1−Xg0T/273.15Patm/P1/γμair+exp−P−Patm/Bliq_P_Tμliq_T_PatmϕXg0/1−Xg0T/273.15Patm/P1/γ+exp−P−Patm/Bliq_P_T.ForP < Pvap, 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(41)μfluid_P_T=μair.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.


## 3.3. Rebound Chamber and Compression Chamber

For a fluid volume, the rate of change of mass is(42)dmdt=dρVdt=ρdVdt+Vdρdt=∑min•−∑mout•,where m is the mass, ρ is the instantaneous density, V is the instantaneous volume, t is the time, min is the input mass, and mout is the output mass.Assuming the instantaneous density is distributed uniformly, equation (42) can be written as follows:(43)∑Qin−∑Qout−dVdt=VρdP/dρ⋅dPdt,where Qin and Qout are the input volume flow rate and the output volume flow rate, respectively.Substituting equation (14) into equation (43) yields(44)dPdt=Bfluid⋅∑Qin−∑Qout−dV/dtV.It should be noted that the flow rate ofQin and Qout 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. The 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(45)dPrebdt=Bfluid_Preb_T⋅qorifice2Preb−qorifice1Preb+Apiston−Arodx•Apiston−Arod⋅Lreb0−x,dPcomdt=Bfluid_Pcom_T⋅qorifice1Pcom−qorifice2Pcom−qorifice3Pcom+qcheckPcom−Apistonx•Apiston⋅Lcom0+x,where Preb and Pcom are, respectively, the pressure of rebound and compression chamber; qorifice1, qorifice2, and qorifice3 are, respectively, the flow rate through the orifice1, orifice2, and orifice3; qcheck is the flow rate through the check valve; Apiston is the crossing area of the piston; Arod 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 Lreb0 and Lcom0 are, respectively, the initial length of the rebound and compression chamber, as shown in Figure 7(a). Bfluid_Preb_T is the bulk modulus of the fluid at Preb and T and Bfluid_Pcom_T is the bulk modulus of the fluid at Pcom and T.Figure 7
The 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.
(a)
(b)These flow rates ofqorifice1, qorifice2, qorifice3, and qcheck can be calculated through equation (3). Significantly, the density and kinematic viscosity used for those above flow rates are evaluated at the mean pressurePup+Pdown/2. These 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(46)qorifice1Pcom=qorifice1⋅ρPup+Pdown/2ρPcom.


## 3.4. Reservoir Chamber

The 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. Figure7(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 bagVgas0 is(47)Vgas0=∑i=1kVk−π8Apiston2−Arod2⋅Lc0−Lp,where k is the number of the air bag, Vk is the precharge volume of the air bag, Lc0 is the active length of the pressure tube, and Lp is the width of the piston, as shown in Figure 7(a).According to the isothermal lawP1V1 = P2V2, the initial pressure Pgas0 is calculated as follows:(48)Pgas0=Patm⋅∑1kVkVgas0.The air in the reservoir is assumed to obey the polytropic gas law of the formP⋅Vγ=constant. Differentiating the polytropic equation, the air pressure derivative is calculated as follows:(49)dPgasdt=γ⋅PgasVgas⋅dVgasdt,where Pgas is the air instantaneous pressure and Vgas the total air instantaneous volume.Assuming the compressibility of the fluid is ignored, the air volume derivative is(50)dVgasdt=qcheckPgas−qorifice3Pgas.Using equations (47)–(50), the states of the reservoir chamber can be described.


## 3.5. 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. The integration method is the SODASRT2. The new hydraulic damper model contains 56 inputs with definite physical meanings, 6 force states, and 47 outputs. Figure8 shows the flow chart of the characteristic calculation of the hydraulic damper. Main simulation parameters of the proposed model are described in Table 1.Figure 8
Flow chart of the characteristic calculation of the hydraulic damper.Table 1
Main simulation parameters of the proposed model.
Parameters
Values
Parameters
Values
Diameter of orifice1
0.5 mm
Length of orifice1
1.75 mm
Diameter of orifice2
0.5 mm
Length of orifice2
1.75 mm
Diameter of orifice3
0.4, 0.5, 0.6 mm
Length of orifice3
0.75 mm
Diameter of pressure tube
60 mm
Diameter of piston
59.98 mm
Diameter of rod
30 mm
Width of piston
35 mm
X
g0
0.1 %
P
sat
75 bar
P
vap
−0.6 bar
ρ
liq
_
P
atm
_
273
854 kg/m3
ρ
air
_
P
atm
_
273
1.2 kg/m3
B
liq
_
P
_
T
15281 bar
T
25°C
γ
1.4
T
1
−20°C
ν
liq
_
T
1
_
P
atm
1430 cSt
T
2
100°C
ν
liq
_
T
2
_
P
atm
6.1 cSt
μ
air
0.02 cP
k
4
L
reb0
168 mm
L
com0
169 mm
L
c0
372 mm
V
k
110 ml


## 4. 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 Figure9(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. The 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”. The geometry dimensions of the modified damper are shown in Table 1.Figure 9
Illustration of the test apparatuses and the orifice. (a) Test bench of SPTR100. (b) Hydraulic fixture. (c) Test orifices.
(a)
(b)
(c)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 Figure9(c) are installed into the base valve unit so as to obtain as much experimental data as possible during the orifice-working stage.


### 4.1. Static Characteristics

The 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. Figure10 shows these comparisons of the force-displacement diagram under the sinusoidal amplitude of ±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 force-velocity 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. The 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.Figure 10
Comparisons of the force-displacement diagram with the amplitude ±5 mm at different excitation frequencies (diameter of orifce3 is 0.4 mm). (a) 0.2 Hz. (b) 0.3 Hz. (c) 0.4 Hz. (d) 0.5 Hz.
(a)
(b)
(c)
(d)Figure 11
Comparisons of the force-velocity diagram with the amplitude ±5 mm at different excitation frequencies (diameter of orifce3 is 0.4 mm). (a) 0.2 Hz. (b) 0.3 Hz. (c) 0.4 Hz. (d) 0.5 Hz.
(a)
(b)
(c)
(d)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 Figures12 and 13 under the excitation 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. The 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. The 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 orifice-working stage, the effect of relief valves are beyond the scope of this article and there is no point in comparison. The simulation reappearance of the pressure limited effect and the fluid shortage indicates that the fluid modelling is accurate.Figure 12
Comparisons of the force-displacement diagram with the amplitude ±5 mm at different excitation frequencies (diameter of orifce3 is 0.5 mm). (a) 0.2 Hz. (b) 0.3 Hz. (c) 0.4 Hz. (d) 0.5 Hz.
(a)
(b)
(c)
(d)Figure 13
Comparisons of the force-displacement diagram with the amplitude ±10 mm at different excitation frequencies (diameter of orifce3 is 0.6 mm). (a) 0.2 Hz. (b) 0.3 Hz. (c) 0.4 Hz.
(a)
(b)
(c)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.


### 4.2. Dynamic Characteristics

The 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. Figures14 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 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.Figure 14
Comparisons of the force-displacement diagram with the amplitude ±0.5 mm at different excitation frequencies (diameter of orifce3 is 0.4 mm). (a) 2 Hz. (b) 5 Hz. (c) 10 Hz. (d) 15 Hz.
(a)
(b)
(c)
(d)Figure 15
Comparisons of the force-displacement diagram with the amplitude ±1 mm at different excitation frequencies (diameter of orifce3 is 0.4 mm). (a) 2 Hz. (b) 4 Hz. (c) 6 Hz.
(a)
(b)
(c)Like the static characteristic, different diameters of the base orifice are used for the purpose of the adequate verification of the established mathematical model. Figure16 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. The 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.Figure 16
Comparisons of the force-displacement diagram with the amplitude ±0.5 mm at different excitation frequencies (diameter of orifce3 is 0.5 mm). (a) 2 Hz. (b) 5 Hz. (c) 10 Hz. (d) 15 Hz.
(a)
(b)
(c)
(d)Figure 17
Comparisons of the force-displacement diagram with the amplitude ±1 mm at different excitation frequencies (diameter of orifce3 is 0.6 mm). (a) 3 Hz. (b) 5 Hz. (c) 7 Hz.
(a)
(b)
(c)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.


### 4.3. 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 Figure18. 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 orifice-working stage.Figure 18
The flow coefficient of different orifices with the amplitude ±0.5 mm at 4 Hz: (a) orifice1 and (b) orifice3.
(a)
(b)As stated previously, fluid properties are crucial to the modelling of the hydraulic damper. Similarly, as an example, Figure19 illustrates the fluid properties of the rebound 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. The 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 volumetric content also has a great effect on the density related to the flow rate of the orifice. The 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.Figure 19
The fluid characteristics of the rebound chamber with the amplitude ±10 mm at 0.4 Hz (diameter of orifce3 is 0.6 mm). (a) Ddensity. (b) Bulk modulus. (c) Absolute viscosity. (d) Air/gas volumetric content.
(a)
(b)
(c)
(d)


## 4.1. Static Characteristics

The 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. Figure10 shows these comparisons of the force-displacement diagram under the sinusoidal amplitude of ±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 force-velocity 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. The 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.Figure 10
Comparisons of the force-displacement diagram with the amplitude ±5 mm at different excitation frequencies (diameter of orifce3 is 0.4 mm). (a) 0.2 Hz. (b) 0.3 Hz. (c) 0.4 Hz. (d) 0.5 Hz.
(a)
(b)
(c)
(d)Figure 11
Comparisons of the force-velocity diagram with the amplitude ±5 mm at different excitation frequencies (diameter of orifce3 is 0.4 mm). (a) 0.2 Hz. (b) 0.3 Hz. (c) 0.4 Hz. (d) 0.5 Hz.
(a)
(b)
(c)
(d)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 Figures12 and 13 under the excitation 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. The 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. The 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 orifice-working stage, the effect of relief valves are beyond the scope of this article and there is no point in comparison. The simulation reappearance of the pressure limited effect and the fluid shortage indicates that the fluid modelling is accurate.Figure 12
Comparisons of the force-displacement diagram with the amplitude ±5 mm at different excitation frequencies (diameter of orifce3 is 0.5 mm). (a) 0.2 Hz. (b) 0.3 Hz. (c) 0.4 Hz. (d) 0.5 Hz.
(a)
(b)
(c)
(d)Figure 13
Comparisons of the force-displacement diagram with the amplitude ±10 mm at different excitation frequencies (diameter of orifce3 is 0.6 mm). (a) 0.2 Hz. (b) 0.3 Hz. (c) 0.4 Hz.
(a)
(b)
(c)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.


## 4.2. Dynamic Characteristics

The 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. Figures14 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 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.Figure 14
Comparisons of the force-displacement diagram with the amplitude ±0.5 mm at different excitation frequencies (diameter of orifce3 is 0.4 mm). (a) 2 Hz. (b) 5 Hz. (c) 10 Hz. (d) 15 Hz.
(a)
(b)
(c)
(d)Figure 15
Comparisons of the force-displacement diagram with the amplitude ±1 mm at different excitation frequencies (diameter of orifce3 is 0.4 mm). (a) 2 Hz. (b) 4 Hz. (c) 6 Hz.
(a)
(b)
(c)Like the static characteristic, different diameters of the base orifice are used for the purpose of the adequate verification of the established mathematical model. Figure16 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. The 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.Figure 16
Comparisons of the force-displacement diagram with the amplitude ±0.5 mm at different excitation frequencies (diameter of orifce3 is 0.5 mm). (a) 2 Hz. (b) 5 Hz. (c) 10 Hz. (d) 15 Hz.
(a)
(b)
(c)
(d)Figure 17
Comparisons of the force-displacement diagram with the amplitude ±1 mm at different excitation frequencies (diameter of orifce3 is 0.6 mm). (a) 3 Hz. (b) 5 Hz. (c) 7 Hz.
(a)
(b)
(c)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.


## 4.3. 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 Figure18. 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 orifice-working stage.Figure 18
The flow coefficient of different orifices with the amplitude ±0.5 mm at 4 Hz: (a) orifice1 and (b) orifice3.
(a)
(b)As stated previously, fluid properties are crucial to the modelling of the hydraulic damper. Similarly, as an example, Figure19 illustrates the fluid properties of the rebound 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. The 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 volumetric content also has a great effect on the density related to the flow rate of the orifice. The 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.Figure 19
The fluid characteristics of the rebound chamber with the amplitude ±10 mm at 0.4 Hz (diameter of orifce3 is 0.6 mm). (a) Ddensity. (b) Bulk modulus. (c) Absolute viscosity. (d) Air/gas volumetric content.
(a)
(b)
(c)
(d)


## 5. 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. The 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.(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)
The Maxwell model cannot accurately reproduce the static and dynamic characteristics.(5)
The modelling of leakages is suggested to be taken into consideration under the high frequency.(6)
The short orifice and thin-bladed orifice both have obvious nonlinear characteristics, and the nonlinear flow characteristics of them should be taken seriously rather than a simple constant value of the flow coefficient.(7)
The 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)
The characteristics of the fluid are not steady but dynamic. The 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.



---

*Source: 1830150-2020-01-06.xml*