Shock and Vibration
Volume 2020 (2020), Article ID 8886631, 21 pages
https://doi.org/10.1155/2020/8886631

Research Article

High-Precision Dynamics Characteristic Modeling Method Research considering the Influence Factors of Hydropneumatic Suspension

Wenguang Wu,1,2 Hongliang Tang,1,2 Sha Zhang,³ Lin Hu,1,2 and Fanhao Zhang1,2

¹School of Automotive & Mechanical Engineering, Changsha University of Science & Technology, Hunan 410000, China
²Hunan Key Laboratory of Smart Roadway and Cooperative Vehicle-infrastructure Systems, Changsha University of Science and Technology, Hunan 410000, China
³CRRC Zhuzhou Institute, Hunan 412000, China

Correspondence should be addressed to Lin Hu; hulin@csust.edu.cn

Received 15 June 2020; Revised 1 July 2020; Accepted 9 July 2020; Published 29 July 2020

Academic Editor: Hong-jun Zhu

Copyright © 2020 Wenguang Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In recent years, hydropneumatic suspension (HPS) has come into widespread use for improving the ride comfort and handling of mining dump trucks and off-road vehicles. Therefore, it is critical to improve the mathematical modeling accuracy to enhance the design and control efficiency and accuracy of HPS. This paper aims to propose a model for improving the modeling precision by considering the effect of different factors on HPS characteristics. A computational fluid dynamic (CFD) model of a HPS was developed, and the volume of fluid (VOF) method was used for the transient calculations in order to simulate the fluid dynamic characteristics and determine the damping and stiffness forces of HPS. The effect of temperature, oil viscosity, nitrogen dissolution rate, and suspension vibration speed on the nonlinear characteristics of HPS was investigated. A limited number of simulation sample points were designed based on the variation ranges of the above factors using the design of experiment (DOE) method. The corresponding damping and stiffness force of each sample point were calculated by CFD simulation. The obtained simulation data were utilized for the fitting of a Kriging model. The results demonstrated that the Kriging model can provide high accuracy, with a prediction error lower than 5%. The proposed modeling method of the HPS nonlinear characteristics is highly efficient, accurate, and faster than traditional methods.

1. Introduction

Due to their nonlinear stiffness and damping characteristics, hydropneumatic suspensions (HPSs) are able to attenuate the body vibration and impact caused by road unevenness in a transitory period, improving handling stability and ride comfort [1]. Mining dump trucks and off-road vehicles are usually equipped with HPSs. In order to improve the control efficiency and accuracy of HPSs under instantaneous working conditions, it is necessary to build a mathematical model that considers the effect of several factors on HPS nonlinear characteristics. Some studies have suggested the oil temperature, oil viscosity, gas and hydraulic oil dissolution, and suspension vibration speed as the key factors that affect the nonlinear characteristics of HPSs [2–4]. Therefore, to accurately model the HPS nonlinear characteristics, these factors should be taken into consideration.

Most of the vibration energy absorbed by the suspension system is dissipated in the form of thermal energy, which leads to the increase of temperature in the HPS. Thermodynamic models have always been used to evaluate the temperature increment [3, 5]. A lumped parameter thermal model considering the heat capacity of the cylinder and piston rod was developed to investigate the temperature increase in a mining dump truck [6]. Moreover, another study demonstrated that the hydraulic oil in the HPS can be divided into several regions, since the oil temperature is gradient. Their experiments indicated that the multiregion thermodynamic model method could improve the accuracy of temperature dynamics [7].

At present, the ideal and real gas law approaches have been employed to calculate the effect of HPS gas state change [8, 9]. The Benedict–Webb–Rubin model, which includes the energy equation, has exhibited strong accuracy; however, it does not take into consideration the effect of gas-oil emulsion [10]. In [2], an analytical model of the HPS strut was proposed considering polytropic change in the gas state and gas-oil emulsion in the HPS. The acquired data exhibited a decreasing trend in the HPS damping and stiffness characteristics, which was related to the fact that the fluid density and effective bulk modulus underwent a decrease due to gas-oil emulsion. Moreover, an analytical model was developed considering the gas-oil emulsion flows through damping valves. It was reported that, due to high fluid pressure, the single orifice configuration provides greater gas entrapment within the oil and thus significantly higher compressibility of the gas-oil emulsion [9].

Furthermore, hydraulic oil viscosity changes with the variation of temperature, which significantly affects the HPS damping characteristic [11]. It has been shown that oil viscosity decreases gradually with the increase of temperature in the oil shock absorber, having a great effect on its damping performance [12]. In [13], the GenShock active suspension system prototype was used to investigate the effects of fluid viscosity variation on the power transmission of the shock absorber and the energy recovery efficiency of the suspension system. The optimal fluid viscosity that enables the system to achieve the best performance was determined.

Geometrical parameters and operating conditions have been proven to have notable effect on the HPS nonlinear characteristics. For example, the diameter of the damping valve has active effect on damping force [14]. The stiffness is affected by the gas volume and the area of the piston and piston rod [15]. The HPS damping characteristics are mainly affected by the structural parameters, such as the rebound chamber and damping valve area and the check valve effective flow area [16].

Based on approaches reported in recent studies, the research cost, time, and complexity can be reduced and the modeling accuracy can be improved by combining statistical theory with computer simulation [17]. According to [18–22], the optimization of the research object can be achieved by combining the design of experiment (DOE) method with the approximate model technology, which can significantly reduce the number of experiments and design costs. Furthermore, it has been proven that computation fluid dynamics (CFD) is an effective and accurate way to simulate the performance of the shock absorber. In [23], the CFD method was employed to develop a dynamic shock absorber model and simulate its performance. The results exhibited good agreement between simulation and experiment. In [24], the pressure distribution acting on the components of a deformable damper valve was numerically investigated using CFD to predict HPS performance. Moreover, since a single-chamber HPS operates through the interaction between oil and gas, the volume of fluid (VOF) method can be appropriately applied, which can be used to solve the free surface problem [25, 26] and simulate the gas-oil interaction in the HPS.

Previous studies on the mathematical modeling of the HPS nonlinear characteristics have been mostly based on mathematical derivation and experimental verification [27, 28], taking only few influence factors of the HPS nonlinear characteristics into consideration, while ignoring the fact that the HPS internal flow field varies over time. Consequently, the modeling accuracy is poor. In order to accurately model the nonlinear characteristics of HPSs, experiments must be performed under real working conditions involving many factors. Due to that experimental equipment is usually associated with high costs and limitations, it is unrealistic to conduct experiments to study the HPS nonlinear characteristics considering multiple factors. Therefore, it is of particular importance to propose a high-precision and fast method for the mathematical modeling of HPSs.

In this study, a novel method for HPS modeling is proposed. In particular, the principal motivation for this study is to propose a novel modeling method for HPS, which combines CFD simulations with the approximate modeling, and can serve as guidance in designing and controlling HPSs. A flowchart of the overall modeling method is illustrated in Figure 1. The present work aims to provide insights into which parameters are the most sensitive to the HPS nonlinear characteristics and how to build a functional or interactional relationship between these parameters and the nonlinear characteristics. Moreover, the transient flow field characteristics in the HPS are analyzed, including the variations of fluid volume fraction, velocity, and pressure distribution.

Figure 1: Flowchart of the overall mathematical modeling method of HPS.

This paper is assembled as follows: the structure, working principle, and mathematical modeling method of the HPS under investigation are introduced in Section 2. An assessment of the CFD simulation process is addressed in Section 3. In Section 4, the simulation results of the HPS under different working conditions are first demonstrated and then approximated by the Kriging model. The fitting and prediction precision of the Kriging model are also validated. Subsequently, the accuracy of the proposed model is verified in a real car experiment. Finally, the conclusions and some future recommendations are drawn in Section 5.

2. Dynamic Characteristics of HPS

2.1. Structure and Working Principle

The structure diagram of an HPS with single chamber is illustrated in Figure 2, where three main components, namely, cylinder, piston, and piston rod, can be observed. Damping valves and check valves are symmetrically distributed in pairs on the piston rod. Two cavities are formed in the inner volume of the HPS: a compression cavity (I) and a rebound cavity (II). The hydraulic oil and the nitrogen contained in cavity (I) are not separated by a diaphragm. During the compression stroke, the check valves are opened, and the hydraulic oil flows through the damping and check valves from cavity (I) to cavity (II). In this way, the stiffness force is larger than the damping force, the vibration of the vehicle is attenuated by the HPS, and massive heat is generated around the damping valve area. During the extension stroke, the check valves are closed, and the oil flows from cavity (II) to cavity (I) only through the damping valves. As the volume of cavity (I) expands, the volume of the contained nitrogen increases, leading to a decrease in pressure. The pressure in cavity (II) increases rapidly due to the compression of the hydraulic oil. At this stroke, the damping force is larger than the stiffness force, which is helpful to attenuate the vibration and enable vehicle stabilization.

Figure 2: Structure diagram of the HPS.

The parameters of the HPS are presented in Table 1.

Table 1: Nomenclature of the HPS parameters.

Parameter
Description
Value

Area of the suspension piston

Area of the suspension piston rod

Density of the hydraulic oil

Diameter of the check valve

Diameter of the damping valve

Area of the check valve

Area of the damping valve

Flow rate coefficient

2.2. Mathematical Model of HPS

2.2.1. Stiffness Mathematical Model

In general, nitrogen is regarded as an ideal gas in the calculations, and its state change process is described by the polytropic state equation:

where and are the initial pressure and volume of the HPS, respectively, and are the instantaneous gas pressure and volume during HPS operation, respectively, and is the polytropic constant. In general, when the change of the gas state is regarded as an isothermal process, , while when it is regarded as an adiabatic process, .

The compression force in the cylinder can be calculated as follows:

where is the inner effective bearing area of HPS.

According to the polytropic state equation, the instantaneous volume and pressure of the nitrogen at any time can be calculated by the following equation:

Then, the stiffness of the hydropneumatic suspension can be calculated as follows:

where is the relative displacement between the cylinder and piston rod. As it can be seen in Figure 2, the effective bearing area in a single-chamber HPS does not change with the stroke; therefore,

Consequently, the stiffness of the single-chamber HPS can be obtained as follows:

2.2.2. Damping Mathematical Model

When the suspension cylinder and the piston rod move up and down relative to each other, the oil flow rate between cavity (I) and cavity (II) can be calculated as follows:

where is the cross-sectional area of cavity (II).

When the hydraulic oil flows through the damping valve, the flow is purely turbulent. Thus, it can be assumed to satisfy the flow formula of thin-walled orifice, which can be written as follows:

where is the density of the hydraulic oil, and are the flow coefficients of the damping and check valves, respectively, and are the effective overcurrent areas of the damping and check valves, respectively, and is the sign function. When the suspension is in an extension stroke, ; when the suspension is in a compression stroke, .

The pressure drop at both sides of the damping valve can be calculated as follows:

where is the average oil flow velocity, is the length of the damping valve, and is the diameter of the damping valve.

When the Reynolds number , , the pressure drop can be written as

where is the drag coefficient, which can be calculated by the following equation:

where is the roughness of the damping valve wall.

According to Pascal's law, the damping force of the HPS can be calculated by the following equation:

3. Two-Phase Flow Simulation Method for HPS

The traditional empirical formula does not take into consideration the effect of the HPS structure and working conditions on stiffness and damping force. However, flow characteristics of the fluid in the HPS, such as pressure and velocity, at different time points can be obtained through CFD simulation, which can be used to calculate the HPS damping and stiffness forces and analyze the nonlinear characteristics. Considering that the HPS under investigation is a single-chamber suspension with gas-oil emulsion and that the hydraulic oil flows freely and the nitrogen state changes with stroke, the changes at the interface between nitrogen and oil need to be tracked. Therefore, the VOF method was employed for the simulations in this study.

3.1. VOF Fluid Phase Interface Tracking Method

The VOF method can simulate the flow and interaction of two or more immiscible fluids by solving a single set of momentum equations and tracking the volume fraction of each fluid throughout the domain. The VOF formulation relies on the fact that the different fluids (or phases) are not interpenetrating. For each additional phase that is added to the model, its volume fraction is introduced in the computational cell. If the volume fraction of the qth fluid in the cell is denoted as αq, then the following three conditions are possible:

 : the cell is empty (of that fluid) : the cell is full (of that fluid) : the cell contains one or more interfaces between that fluid and one or more other fluids

The governing equations of the VOF method include the volume fraction equation, the additional scalar equations, and the continuum surface force model.

3.1.1. Volume Fraction Equation

The tracking of the interface between the phases is accomplished by solving a continuity equation for the volume fraction of one (or more) of the phases. For the fluid, this equation takes the following form:

where is the mass migrated from phase to phase , is the mass migrated from phase to phase , is the source term, is the density of phase , and is the velocity of phase .

The volume fraction equation is not solved for the primary phase; the primary-phase volume fraction is computed based on the following constraint:

3.1.2. Additional Scalar Equations

The density and viscosity of the mixed fluid in the unit are calculated as follows:

where is the density of the gas, is the density of the liquid, is the viscosity of the gas, and is the viscosity of the liquid.

3.1.3. Continuum Surface Force Model

The surface tension force is generated by the gravitational force between the molecules in the fluid. It acts only on the surface and minimizes the surface free energy by reducing the interface area. In CFD analysis, the continuum surface force (CSF) model proposed by Brackbill has been used to implement the effect of surface tension force. In this model, the pressure drop across the surface depends upon the surface tension coefficient and the surface curvature as measured by two radii in orthogonal directions, and :

where and are the pressures of the two fluids on either side of the interface.

The surface tension force can be expressed in terms of the pressure jump across the surface. Using the divergence theorem, the force at the surface can be expressed as a volume force as follows:

where is the surface tension coefficient between phase i and phase j, and are the volume fraction and density of phase i (the primary phase is usually liquid), respectively, and are the volume fraction and density of phase j (the secondary phase is usually gas), respectively, and and are the surface curvatures of phase i and phase j, respectively.

If only two phases are present in a cell, the above equation can be simplified as follows:

3.2. Nitrogen Dissolution in Oil

The volume and pressure of cavity (I) and cavity (II) change dynamically during the compression and extension strokes. In this process, nitrogen dissolves partially in the oil. Due to its low density, the intermolecular distance between nitrogen molecules is relatively larger than that of hydraulic oil molecules, and hence, the nitrogen dissolved in the oil will lead to a decrease in hydraulic oil density and a slight increase in volume.

The dissolution rate of nitrogen in hydraulic oil is defined as the volume of dissolved nitrogen as a percentage of the initial nitrogen volume in HPS. It can be shown as follows:

where is the initial volume of nitrogen in the HPS and is the volume of nitrogen dissolved in hydraulic oil.

The density of the gas-oil mixture can be calculated as follows:

where is the total mass of oil in the HPS and and are the initial volume and density of the oil, respectively. Based on the volume of the dissolved gas, the volume of nitrogen and the density of hydraulic oil were adjusted in the Fluent software to simulate the transient characteristics of gas-oil emulsion.

3.3. Design and Meshing of the HPS Model

3.3.1. HPS Geometry Design

Based on the actual HPS structure of a specific mining dump truck, the geometry model was designed using UG software. The internal flow domain of the HPS was extracted for the CFD analysis (as shown in Figure 3). Thereby, the outer surface of the isolated flow domain was the boundary of the fluid motion in the calculations.

Figure 3: 3D model and isolated flow domain of the HPS.

3.3.2. Meshing

The unstructured grid method was used to mesh the volume of the flow domain. Since the size of the damping and check valves was smaller than that of the whole model, the mesh at these regions was finer. Finally, the total number of elements in the model was more than 1.6 million. Figure 4 shows a cross-section of the generated mesh of the model.

Figure 4: Cross-section of the generated mesh of the model.

3.4. Numerical Calculation

3.4.1. Model Setup

The ANSYS Fluent software was used to perform the transient calculation of the flow in the HPS, and the parameters set to the designed VOF model are listed in Table 2. The oil used in the HPS was L-HM46, which can be regarded as Newtonian fluid. Thus, the oil was modeled as incompressible Newtonian fluid.

Table 2: VOF model parameters.

VOF model properties
Settings

Primary/secondary phase
Nitrogen gas/oil

Viscous model
RNG k-ε model

Gas state
Ideal gas

Oil compressibility
Incompressible Newtonian fluid

Precharge nitrogen pressure ()
6.276

Precharge nitrogen height ()
100

Temperature ()
300

Oil dynamic viscosity ()
0.1059

Oil density ()
994

The dynamic mesh method was used to simulate the relative movement between cylinder and piston rod. The moving boundaries were the cavity (I) upper wall and the cavity (II) bottom wall, as shown in Figure 5(a). The boundary movement was prescribed using a profile file in the Fluent software. The motion functions of the moving boundary in compression and extension stroke are presented in Figures 5(b) and 5(c), respectively. Since during the compression and extension strokes, the boundary moves in opposite directions; in compression stroke, the boundary velocity was defined by a negative value, and in extension stroke by a positive value.

(a)

(b)

(c)

Figure 5: (a) The defined moving boundaries; (b) plot of boundary velocity versus time in compression stroke; (c) plot of boundary velocity versus time in extension stroke.

The calculation settings for the transient simulation are listed in Table 3.

Table 3: Calculation settings for the transient simulation.

Solving method
Settings

Solver type
Pressure-based

Time
Transient

Near-wall treatment
Standard wall function

Dynamic mesh method
Smoothing, remeshing

Pressure-velocity coupling scheme
PISO

Pressure spatial discretization
Presto!

Density spatial discretization
First-order upwind

Volume fraction spatial discretization
Georeconstruct

Transient formulation
First-order implicit

3.4.2. Simulation Results and Flow Field Analysis

T = 0.2, 0.4, 0.8, and 1.0 s were chosen as representative time points to present the dynamic volume ratio of oil and nitrogen in compression stroke as shown in Figure 6. During this process, the nitrogen in cavity (I) is compressed, and its volume is decreased. The oil in cavity (I) flows to cavity (II) through the check and damping valves, and the volume of cavity (II) becomes increased.

(a)

(b)

(c)

(d)

Figure 6: Volume changes of oil and nitrogen during the compression stroke at (a) t = 0.02 s, (b) t = 0.04 s, (c) t = 0.08 s, and (d) t = 0.1 s. The blue area is the hydraulic oil, the red area is the nitrogen, and the green area is the interface between the two phases.

The internal flow field distribution characteristics of the HPS can be obtained based on the simulation data. The contours of velocity distribution at a cross-section of the damping valve at different time points during the compression stroke are illustrated in Figures 7(a) and 7(b). It can be seen that the velocity distribution of the flow was symmetrical, and the velocity of the oil increased sharply after flowing from cavity (I) through the damping and check valves and gradually decreased after entering cavity (II). When the vibration speed of the suspension was 0.5  at t = 0.05 s, the maximum oil flow velocity was 42 . However, at t = 0.1 s, the maximum fluid velocity reached 80 , while the vibration speed of the suspension reached 1 , which was twice as fast as the velocity of the oil at t = 0.05 s. Despite that the flow area of the damping valves was smaller than that of the check valves, the maximum fluid velocity occurred at the junction area between cavity (II) and the check valves. This is due to that the flow rate through the check valves was greater than that through the damping valves.

(a)

(b)

Figure 7: Velocity distribution contours at a cross-section of the damping valve in compression stroke at (a) t = 0.05 s and (b) t = 0.1 s.

Figure 8(a) and 8(b) display the pressure distribution contours at a cross-section of the damping valve at t = 0.05 s and t = 0.1 s, respectively. It can be seen that the pressure distribution in the suspension was also symmetrical, and the pressure in cavity (I) was higher than that in cavity (II) during the compression stroke. This can be attributed to the energy dissipation when the oil flows from cavity (I) to cavity (II) and passes through the damping and check valves. At t = 0.05 s, the maximum pressure was 7.2 , and the pressure difference between cavity (I) and cavity (II) was 1.2 . At t = 0.1 s, the maximum pressure reached 12.5 , and the pressure difference between the two cavities was 4 . It can be seen that the pressure difference between the two cavities increased with the increase of the suspension vibration speed.

(a)

(b)

Figure 8: Pressure distribution contours at a cross-section of the damping valve in compression stroke at (a) t = 0.05 s and (b) t = 0.1 s.

The volume ratios of oil and nitrogen in the HPS at representative time points (t = 0.2, 0.4, 0.8, and 1.0 s) during the extension stroke are presented in Figure 9. During this working condition, the check valves are closed, and the oil flows from cavity (II) to cavity (I) through the damping valve. Therefore, the volume of the oil in cavity (II) decreases and the pressure increases, while the volume of the oil in cavity (I) increases significantly. At the same time, the nitrogen in cavity (I) expands due to the pressure decrease in cavity (I).

(a)

(b)

(c)

(d)

Figure 9: Volume changes of oil and nitrogen during the extension stroke at (a) t = 0.02 s, (b) t = 0.04 s, (c) t = 0.08 s, and (d) t = 0.1 s. The blue area is the hydraulic oil, the red area is the nitrogen, and the green area is the interface between the two phases.

The velocity distribution contours of flow at a cross-section of the damping valve at t = 0.05 s and t = 0.1 s during the extension stroke are shown in Figures 10(a) and 10(b), respectively. During this process, the hydraulic oil flows from cavity (II) to cavity (I) only through the damping valves, while the check valves are closed; thus, the flow velocity at the check valve cross-sectional area was equal to zero. The maximum velocity at t = 0.05 s and t = 0.1 s was 95  and 190 , respectively, which was higher than that of the compression stroke. Furthermore, it can be observed that the oil flow in cavity (I) was very complex and irregular, with symmetrically generated flow vortices. As the suspension vibration speed increased from 0.5  to 1 , the vortices became more apparent. In general, vortices convert the vibration energy of the suspension into heat loss, which is helpful to attenuate the vibration of the vehicle.

(a)

(b)

Figure 10: Velocity distribution contours at a cross-section of the damping valve in extension stroke at (a) t = 0.05 s and (b) t = 0.1 s.

The pressure contours at a cross-section of the damping valve in extension stroke at t = 0.05 s and t = 0.1 s are shown in Figures 11(a) and 11(b), respectively. It is apparent that the pressure of cavity (II) was higher than that of cavity (I), indicating the hydraulic oil in cavity (II) was compressed. At t = 0.05 s, the suspension vibration speed was 0.5  and the maximum pressure was 16 , which was much higher than that in the compression stroke (7.2 ). When the suspension vibration speed increased to 1  at t = 0.1 s, the maximum pressure reached 50 . As the suspension vibration speed increased from 0.5  (t = 0.05 s) to 1  (t = 0.1 s), the pressure difference between cavity (I) and cavity (II) increased from 10  to 45 , respectively.

(a)

(b)

Figure 11: Pressure distribution contours at a cross-section of the damping valve in extension stroke at (a) t = 0.05 s and (b) t = 0.1 s.

3.5. Stiffness and Damping Calculation

Through the above simulation analysis, the transient pressure distribution in the HPS under specific operating conditions can be determined. Consequently, the transient pressures in cavity (I) and cavity (II) can be obtained to calculate the stiffness and damping of the HPS.

If the friction between cylinder and piston rod can be ignored, the output force of the HPS under external excitation can be calculated as follows:

where and are the average total pressures of cavity (I) and cavity (II), respectively, and and are the cross-sectional areas of cavity (I) and cavity (II), respectively.

Since the pressure difference between cavity (I) and cavity (II) is

then equation (21) can be written as

The output force of the HPS is divided into two parts, stiffness force and damping force , as follows:

Consequently, as long as the transient pressures of cavity (I) and cavity (II) are obtained, the stiffness and damping force of the HPS can be accurately calculated.

4. Surrogate Model Construction and Verification

4.1. Surrogate Modeling

The specific implementation scheme of the mathematical modeling process for HPS is presented in Figure 12.

Figure 12: Flowchart of the HPS modeling process.

The surrogate modeling method is an analysis method used in the field of multidisciplinary design optimization (MDO) founded in recent years. It is a highly efficient method for determining the functional or interactional relationship between the input and output of a studied system. In the case of HPS design, the DOE method can be employed to design a large number of simulation sample points of several different factors that affect the HPS nonlinear characteristics. Subsequently, according to the working conditions of these sample points, several groups of CFD simulations can be performed to obtain the HPS damping and stiffness forces under various working conditions. Finally, the simulation data can be fitted by surrogate models, which can build an approximate function relationship between HPS nonlinear characteristics and their influence factors.

4.2. Kriging Model

The Kriging model, which was employed in this study, is an interpolation method developed in the field of spatial statistics and geostatistics. It predicts the distribution of function values at unknown points instead of the function values themselves. From the distribution of the function values, both function values and their uncertainty at new points can be estimated. Simpson et al. [29] indicated that the Kriging model is the best choice for deterministic and highly nonlinear problems with a moderate number of variables (<50). The Kriging model can be given as follows [17]:

where is the sample point with m variables, is an approximate function fitted to the n sample points, is a linear or nonlinear function of , is the regression coefficient to be estimated, and is a stochastic function with mean zero and variance .

4.3. Optimal Latin Hypercube Experiment Design

The optimal Latin hypercube experiment design was used to provide a uniform distribution of the experimental sample points in the experiment space. In this study, the effects of oil temperature (T), oil viscosity (U), gas dissolution rate (D), and suspension vibration speed () on the HPS nonlinear characteristics were taken into consideration. The hydraulic oil used in the HPS was L-HM46. As it can be seen in Figure 13, the viscosity of the L-HM46 hydraulic oil varies considerably with temperature. Thus, it is necessary to investigate the effect of oil viscosity variation with temperature on the HPS damping and stiffness characteristics.

Figure 13: L-HM46 hydraulic oil viscosity variation with temperature.

The variation range of each factor was determined empirically (Table 4). The temperature range is between 25 and 65. The hydraulic oil viscosity, which varies over the temperature, ranges from 0.0159 to 0.1252 . In [30], it has been indicated that the amount of nitrogen dissolved in hydraulic oil can be calculated as 78% of the amount of air that can be dissolved. Therefore, in this study, the variation range of gas dissolution rate was determined to be 0 to 10%. The vibration speed, which ranges from to , takes negative values during compression stroke and positive during extension stroke. The Latin hypercube sampling method was used for sampling within the design space of the four factors. 30 groups of sample points with different combinations of values for the four factors were suggested. The sample points were then utilized to perform simulation experiments using the two-phase flow simulation method mentioned in Section 2. The sample points and the corresponding simulation results are listed in Table 5.

Table 4: Variation ranges of the four different factors.

Factors
Variation range

Temperature ()
25∼65

Oil viscosity ()
0.0159∼0.1252

Gas dissolution rate (%)
0∼10

Vibration speed ()
−1.0∼1.0

Table 5: Simulation results.

1.0
8.95
52
0.0258
5.91
29.3
−363043
536437

1.0
4.95
33
0.0686
5.92
30.9
−387672
537618

0.9
3.68
54
0.0239
5.96
25.4
−300640
541068

0.9
5.37
50
0.0288
5.96
25.3
−300514
540714

0.8
5.26
38
0.0507
5.99
21.3
−237626
543828

0.8
7.89
63
0.0179
5.98
21.6
−242167
543193

0.7
0.53
42
0.0408
6.03
17.9
−183556
547769

0.7
8.32
42
0.0408
6.01
17.7
−181170
546126

0.6
2.11
25
0.1252
6.07
15.8
−150268
550820

0.5
6.32
65
0.0159
6.09
12.2
−94817
552954

0.4
8.42
29
0.0914
6.12
10.6
−68643
556096

0.3
7.89
46
0.0338
6.16
8.57
−37284
559609

0.2
1.05
63
0.0179
6.21
7.31
−17129
563423

0.1
3.16
48
0.0308
6.24
6.59
−5377
566556

0.1
4.53
46
0.0338
6.24
6.59
−5371
566519

−0.1
4.21
31
0.0795
6.31
6.24
1093
573284

−0.2
0.00
36
0.0616
6.35
6.15
3046
576462

−0.2
6.63
61
0.0189
6.35
6.17
2798
576961

−0.3
9.47
59
0.0199
6.40
6.04
5530
580675

−0.3
7.47
44
0.0378
6.39
6.03
5684
580557

−0.4
4.74
61
0.0189
6.43
5.83
9322
583790

−0.4
3.68
25
0.1193
6.43
5.75
10473
583662

−0.5
10.0
40
0.0457
6.48
5.57
14164
588575

−0.5
3.26
59
0.0199
6.47
5.56
14113
587176

−0.6
5.79
44
0.0378
6.51
5.22
19931
590963

−0.7
6.84
27
0.1054
6.56
4.77
27643
595303

−0.8
1.58
50
0.0288
6.58
4.36
34340
597219

−0.8
8.74
57
0.0219
6.60
4.43
33657
599407

−0.9
2.63
33
0.0696
6.63
3.81
43575
601623

−1.0
7.37
57
0.0219
6.68
3.34
51739
606771

4.4. Modeling of HPS Nonlinear Characteristics

4.4.1. Model Fitting Evaluation

The simulation results were fitted by the Kriging model and response surface model, respectively. Error analysis on the two models was performed, and the statistical results are shown in Table 6.

Table 6: Evaluation indicators of the Kriging model.

Evaluation indicator
Error acceptance level
Error value

Kriging model
Response surface model

Maximum
0.3
0.095
0.098

Mean
0.2
0.045
0.052

Root mean square
0.2
0.057
0.048

R-square
0.9
0.97
0.95

The root-mean-square error can be calculated as follows:

where and are the observation and real value of the model, respectively, and n is the number of samples.

R-square of the model is given by the following equation:

where is the sample estimates and is the mean value.

From the comparisons between the evaluation index values and the acceptance levels in Table 6, it can be observed that both the Kriging model and response surface model fit the CFD sample points very well. The R-square value of the Kriging model is higher than that of the response surface model. Consequently, the Kriging model fits the simulation results better.

4.4.2. Model Prediction Accuracy Verification

In order to compare the prediction accuracy of the Kriging model and response surface model on the stiffness and damping characteristics of HPS, additional simulations are required. Therefore, based on the variation ranges of the four factors, 10 groups of sample points were randomly generated, and the validity of the model was verified by comparing the simulated values with the predicted values. Tables 7 and 8 present a comparison between simulated and predicted values. It can be seen that the maximum and mean relative error between the simulated value and predicted value of the Kriging model are 2.49% and 0.94%, respectively, while the maximum and mean relative error between the simulated value and predicted value of response surface model are 14.87% and 1.93%, respectively. The prediction accuracy of the Kriging model on damping and stiffness force is higher than that of the response surface model. Therefore, it is more appropriate to model the nonlinear characteristics of HPS by the Kriging model. Meanwhile, in all cases, the relative error of the Kriging model was below 2.56%. Consequently, it can be concluded that the Kriging model developed in this study can accurately predict the damping and stiffness characteristics of HPS under specific working conditions within a certain error range. Therefore, the method proposed in this paper to model the HPS nonlinear characteristics considering multiple influencing factors is validated.

Table 7: Damping force validation.

No.
V ()
D
T (°C)
U ()
Simulation ()
Prediction
Relative error (%)

Kriging model ()
Response surface model ()
Kriging model
Response surface model

1
−1
0.027
50
0.0288
52.44
51.13
53.98
2.49
−2.93

2
−1
0.05
46
0.0338
52.24
51.44
51.19
1.54
2.02

3
−0.6
0.09
32
0.0408
19.79
20.30
19.06
−2.54
3.71

4
−0.6
0.03
50
0.0288
19.93
19.42
19.32
2.56
3.08

5
−0.5
0.05
40
0.0457
14.33
14.41
14.15
−0.54
1.31

6
−0.3
0.036
52
0.0258
5.65
5.54
4.81
1.9
14.87

7
0.6
0.05
60
0.0209
−136.20
−135.70
−138.19
0.37
−1.46

8
0.8
0.07
45
0.0358
−240.95
−237.86
−237.71
1.29
1.35

9
0.9
0.05
46
0.0338
−301.90
−304.48
−298.34
−0.85
1.18

10
1
0.056
46
0.0338
−366.77
−371.06
−372.66
−1.17
−1.61

Table 8: Stiffness force validation.

No.
V ()
D
T (°C)
U ()
Simulation ()
Prediction
Relative error (%)

Kriging model ()
Response surface model ()
Kriging model
Response surface model

1
−1
0.027
50
0.0288
605.32
604.85
605.96
0.08
−0.11

2
−1
0.05
46
0.0338
587.35
589.83
609.08
−0.42
−0.55

3
−0.6
0.09
32
0.0408
591.80
591.98
590.87
−0.03
0.16

4
−0.6
0.03
50
0.0288
590.73
592.89
593.69
−0.36
−0.50

5
−0.5
0.05
40
0.0457
587.35
589.83
589.96
−0.42
−0.44

6
−0.3
0.036
52
0.0258
580.15
580.54
580.96
−0.07
−0.14

7
0.6
0.05
60
0.0209
550.15
547.65
556.30
0.46
−1.12

8
0.8
0.07
45
0.0358
543.20
539.25
538.64
0.73
0.84

9
0.9
0.05
46
0.0338
540.71
537.05
536.05
0.68
0.86

10
1
0.056
46
0.0338
537.61
535.46
539.87
0.4
−0.42

4.4.3. Surrogate Model of HPS Nonlinear Characteristics

The Kriging model cannot be shown as an explicit formulation. The surrogate models of HPS nonlinear damping and stiffness are shown in Figures 14 and 15, respectively, reflecting the relationship between the HPS nonlinear characteristics and the temperature (T), oil viscosity (U), gas dissolution rate (D), and vibration velocity () of the suspension.

(a)

(b)

(c)

(d)

Figure 14: Surrogate models of the relationship between (a) damping force, oil viscosity, and gas dissolution rate; (b) damping force, temperature, and oil viscosity; (c) damping force, suspension vibration speed, and temperature; and (d) damping force, suspension vibration speed, and oil viscosity.

(a)

(b)

(c)

(d)

Figure 15: Surrogate models of the relationship between (a) stiffness force, gas dissolution rate, and oil viscosity; (b) stiffness force, suspension vibration speed, and dissolution rate; (c) stiffness force, suspension vibration speed, and temperature; and (d) stiffness force, suspension vibration speed, and oil viscosity.

Figure 14(a) reflects the effect of the interaction between oil viscosity (U) and gas dissolution rate (D) on damping force. It can be seen that the effect of gas dissolution rate on damping force is different between the extension and compression strokes, exhibiting strong nonlinearity. During the compression stroke, as the nitrogen dissolution rate increases and the damping force increases first and then presents a decreasing trend. However, during the extension stroke, as the nitrogen dissolution rate increases and the damping force decreases first and then presents an increasing trend. This difference is attributed to the changes in the physical properties of the hydraulic oil caused by gas-oil emulsion. Figure 14(b) reflects the effect of the interaction between U and temperature (T) on damping force. It can be noticed that the effects of temperature and oil viscosity on damping force exhibit a consistent trend, since oil viscosity varies with temperature. At the same time, it can be observed that within a certain range, the decrease of oil viscosity leads to an increase in damping force, and then, as the oil viscosity continues to increase and the damping force decreases. The effect of the interaction between suspension vibration speed () and T on damping force is illustrated in Figure 14(c). In this case, compared to T, plays a dominant role in affecting the damping force of the HPS, and the function curve between and Fc on the curved surface is nearly parallel, independent of the value of T. This means that the interaction between and T is weak. Similarly, Figure 14(d) presents the effect of the interaction between V and U on damping force. It can be seen that, compared to U, the effect of on damping force is more apparent, while the interaction between and U is not obvious.

Figure 15(a) reflects the highly nonlinear stiffness characteristics of HPS under the effect of the interaction between oil viscosity (U) and gas dissolution rate (D). It can be observed that the stiffness force reached the minimum value under a low gas dissolution rate and oil viscosity. In addition, the stiffness force reached its peak value when the gas dissolution rate and oil viscosity ranges from 4% to 6% and 0.06 to 0.08, respectively. Figure 15(b) illustrates the effect of the interaction between and D on stiffness force. It can be observed that has a significant impact on stiffness force, which, in compression stroke, increases with the increase of , while in extension stroke, the increase of leads to a decrease in stiffness force. The effect of the interaction between and T on stiffness force is presented in Figure 15(c). It can be clearly observed that, compared to , the variation of T plays a secondary role in affecting the stiffness force. In terms of the degree of parallelism of the Fk- and Fk-T function curves on the surface, no interaction between and T can be identified. Figure 15(d) reflects the effect of the interaction between and U on stiffness force. Through analysis, it can be observed that the effect of on stiffness force is more significant than that of U. The interaction between U and is not strong. Furthermore, in Figures 15(b)–15(d), it can be seen that the stiffness force of the HPS in compression stroke is larger than that in extension stroke.

4.5. Response Sensitivity Analysis

By postprocessing the surrogate model, the response sensitivity to factors can be analyzed. The Pareto plots for the damping and stiffness of the HPS are shown in Figures 16. It can be seen that the factor that has the maximum influence on the damping force is the oil viscosity; next is temperature. Therefore, it can be concluded that the viscosity and temperature effects of hydraulic oil have a significant impact on the damping force of HPS. The variations in suspension vibration speed have the greatest effect on stiffness force; next are the temperature and oil viscosity. The influence of gas dissolution rate on the stiffness force is not apparent.

(a)

(b)

Figure 16: Pareto plots of HPS: (a) damping and (b) stiffness characteristics.

4.6. Experiment and Verification

In order to verify the accuracy of the Kriging model, an experiment using a mining dump truck and a 11-DOF mathematical model are proposed. The accuracy of the traditional HPS model and the proposed HPS model was compared by analyzing the RMS acceleration of the vehicle wheel center and measuring the point of driver’s feet.

4.6.1. Experimental Procedure

The driving experiment using the mining dump truck was performed in an open pit mine in western China (Figure 17). The experimental devices are listed in Table 9. Accelerometers were used to measure the acceleration that separated all of truck, especially the engine and suspensions. Some measurement points are shown in Figure 18. The truck speed ranged between 10 and 30 km/h, and the input was random road.

Figure 17: Mining dump truck and measurement instrumentations.

Table 9: Experimental devices.

Experimental device
Assignment

Battery and invertor
Supply energy for the device

Accelerometer, ,
Measuring the acceleration

Amplifier
Amplify the measured signal

Recorder, MR30 C
Record the signal

Data analysis device and software, SD380
Analyze the signal

Figure 18: Measurement points of experiment.

The experimental results concerning the suspension upper point at a speed of are shown in Figure 19. It can be seen that the amplitude of the rear suspension vibration was larger than that of the front suspension. The root-mean-square values of the acceleration of the front left, front right, rear left, and rear right suspensions were , , , and , respectively.

(a)

(b)

(c)

(d)

Figure 19: Plots of the measured suspension upper point acceleration versus time (at 30 km/h) for the (a) front left suspension; (b) front right suspension; (c) rear left suspension; and (d) rear right suspension.

4.6.2. Simulation and Verification

In order to verify the proposed HPS model, an 11-DOF mathematical model of the mining dump truck was built (Figure 20), and the proposed HPS Kriging model was employed (equations (25) and (26)).

Figure 20: Diagram of the 11-DOF truck model, where is the road excitation; is the unsprung mass of the front and rear suspension; is the deflection of the unsprung mass; is the suspension force, which consists of the stiffness force and the damping force; is the deflection of the sprung mass; is the pitch angle; and is the roll angle.

According to Newton’s Law, the dynamic mathematical equation of vertical axles, the pitch, and roll can be written as follows:

where is the body mass, is the vertical displacement of body mass center, is the pitch moment of inertia, is the roll moment of inertia, and are the height of pitch center and rolling center, respectively, a and b are the distances from the front and rear axles to the body mass center, respectively, and B is the wheel base.

Considering the kinetic coordination of the four suspensions, the suspension deflection equation can be calculated as follows:

where is the suspension deflection.

The plots of RMS acceleration of the four suspensions versus vehicle velocity are illustrated in Figure 21. It can be seen that the results of the mathematical model were very close to the experimental results. About 2–10% of the accuracy was evaluated by the proposed model comparing the traditional model. The accuracy of the proposed model was validated by the indirect verification method.

(a)

(b)

(c)

(d)

Figure 21: Plots of RMS acceleration versus vehicle velocity for the four suspensions. Comparison between the experimental, traditional models, and proposed model results concerning the upper point of the (a) front left suspension; (b) front right suspension; (c) rear left suspension; and (d) rear right suspension.

5. Conclusions

The aim of this study was to propose a high-precision modeling method of the nonlinear characteristics of HPS in which the effect of multiple influence factors can be considered. The oil temperature, oil viscosity, gas dissolution rate, and suspension vibration speed were investigated as the influencing factors of the nonlinear characteristics of HPS. The findings of this study are summarized as follows:

(1)The mathematical modeling method for HPS proposed in this paper, which combines CFD simulation with the approximate model method, provided a highly efficient and accurate way to develop a simple and direct mathematical model describing and predicting the relationship between HPS nonlinear characteristics and the influencing factors regardless of structural changes in HPSs, which will substantially reduce the experimental costs and speed up the development process. This has very important guiding significance for the analysis and development of HPSs.(2)The Kriging model demonstrated excellent performance for approximating the mathematical relationship between HPS nonlinear characteristics and influencing factors. The results suggested that it can provide a good prediction accuracy for the damping and stiffness forces of HPS with the prediction errors not exceeding 5%.

However, this modeling method cannot be expressed explicitly, and the effect of using different gas state equations on the accuracy of the transient calculations is not taken into consideration. In future research work, where more factors will be considered, the applicability of this method needs to be verified. In addition, further investigation is worth to be done to determine how such a method can be applied to enhance active suspension control and autonomous driving, in order to improve driving comfort.

Data Availability

All data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

Special thanks to the National Natural Science Foundation of China (51705035) and the Natural Science Foundation of Hunan Province (2017JJ3336). Contribution was also supported by XEMC.

References

D. Wang, D. Zhao, M. Gong, B. Yang, and C. Trigona, “Nonlinear predictive sliding mode control for active suspension system,” Shock and Vibration, vol. 2018, Article ID 8194305, 10 pages, 2018.

Y. Yin, S. Rakheja, J. Yang, and P. E. Boileau, “Effects of entrapped gas within the fluid on the stiffness and damping characteristics of a hydro-pneumatic suspension strut,” SAE International Journal of Commercial Vehicles, vol. 10, pp. 204–209, 2017.

C. Yi-jie, G. U. Liang, and G. Ji-fu, “Research on the different throttle methods thermodynamic model of hydro-pneumatic spring,” Mechanical Design and Manufacturing, no. 6, pp. 54–56, 2008.

Z. Chang-cheng, Z. Zhi-yun, and G. U. Liang, “Simulation and test of oil viscosity's influence on hydro-pneumatic suspension resistance characteristics,” Computer Simulation, vol. 24, no. 3, pp. 257–261, 2007.

C. Yi-jie, G. U. Liang, and G. Ji-fu, “Research on natural convection thermodynamic model of hydro-pneumatic spring,” Journal of System Simulation, vol. 20, no. 15, pp. 4198–4201, 2008.

X. Huang, Y. Shen, J. Yang, and W. Zhang, “Thermal analysis of hydro-pneumatic suspension system for dumper based on a lumped-parameter thermal model,” Journal of Agricultural Engineering, vol. 10, pp. 64–70, 2013.

H. Xia-xu, Y. Jue, S. Yan-hua, and Z. Wen-ming, “Characteristic analysis of oil temperature in a hydropneumatic suspension based on multiple regions,” Journal of University of Science and Technology Beijing, vol. 37, no. 7, pp. 950–954, 2015.

Y. Yong-zeng and Z. Hui-li, “Analysis based on benedict-webb-rubin equation of state for temperature change characteristics of hydro-pneumatic suspension,” Chinese Hydraulics & Pneumatics, no. 8, pp. 44–53, 2018.

Y. Yin, S. Rakheja, J. Yang, and P.-E. Boileau, “Characterization of a hydro-pneumatic suspension strut with gas-oil emulsion,” Mechanical Systems and Signal Processing, vol. 106, pp. 319–333, 2018.

S. F. van der Westhuizen and P. Schalk Els, “Comparison of different gas models to calculate the spring force of a hydropneumatic suspension,” Journal of Terramechanics, vol. 57, pp. 41–59, 2015.

Y. Shen, L. Chen, Y. Liu, X. Zhang, and M. Lallart, “Modeling and optimization of vehicle suspension employing a nonlinear fluid inerter,” Shock and Vibration, vol. 2016, p. 2623017, 2016.

W. U. Guo-xiang, “Study of viscosity & temperature characteristics of oil pressure damper of rail vehicle,” Journal of Lanzhou Jiaotong University, vol. 29, no. 1, pp. 82–86, 2010.

S. N. Avadhany, Analysis of Hydraulic Power Transduction in Regenerative Rotary Shock Absorbers as Function of Working Fluid Kinematic Viscosity, Massachusetts Institute of Technology, Cambridge, MA, USA, 2009.

Yanhua, J. Shen, Zhao et al., “Research on test and simulation of hydro-pneumatic suspension,” in Proceedings of the International Conference on Consumer Electronics, Communications and Networks (CECNet 2011), pp. 678–681, Xianning, China, 2011.

C. Linshan, “Analysis on influence factors of nonlinear characteristics of hydro-pneumatic suspension based on simulink,” Machine Tool & Hydraulics, vol. 45, no. 19, pp. 179–184, 2017.

Z. Sha, G. Zhengqi, Z. Jingkai, X. Ya, and W. Wenguang, “Gas-liquid two-phase flow numerical simulation of a hydro-pneumatic suspension by VOF model and dynamic mesh method,” China Mechanical Engineering, vol. 27, no. 15, pp. 2091–2099, 2016.

Y. Gao and X. Wang, “An effective warpage optimization method in injection molding based on the Kriging model,” The International Journal of Advanced Manufacturing Technology, vol. 37, no. 9-10, pp. 953–960, 2008.

K. Kasza, Ł. Malinowski, and I. Królikowski, “Optimization of pin fin heat sink by application of CFD simulations and DOE methodology with neural network approximation,” International Journal of Applied Mechanics and Engineering, vol. 18, no. 2, pp. 365–381, 2013.

A. Todoroki and T. Ishikawa, “Design of experiments for stacking sequence optimizations with genetic algorithm using response surface approximation,” Compos Struct, vol. 64, no. 3-4, pp. 349–357, 2004.

R. Wang, X. Su, and W. Wang, “Optimization and research of wheel alignment parameters based on approximates model,” Machinery Design & Manufacture, vol. 11, p. 22, 2015.

L. Hu, J. Ou, J. Huang, Y. Chen, and D. Cao, “A review of research on traffic conflicts based on intelligent vehicles,” Ieee Access, vol. 8, pp. 24471–24483, 2020.

L. Hu, X. Hu, J. Wan, M. Lin, and J. Huang, “The injury epidemiology of adult riders in vehicle-two-wheeler crashes in China, Ningbo, 2011-2015,” Journal of Safety Research, vol. 72, pp. 21–28, 2020.

F. G. Guzzomi, P. O'Neill, and A. Tavner, in Proceedings of the 16th Australasian Fluid Mechanics Conference, P. Jacobs, T. McIntyre, M. Cleary et al., Eds., pp. 1123–1130, School of Engineering, The University of Queensland, Gold Coast, Australia, 2017.

F. Herr, T. Mallin, J. Lane, and S. Roth, “A shock absorber model using CFD analysis and Easy5,” SAE Transactions, pp. 2421–2434, 1999.

L. Xu, H. Liu, N. Si, and E. W. M. Lee, “Numerical simulation of a two-phase flow in the electrospinning process,” International Journal of Numerical Methods for Heat & Fluid Flow, vol. 24, no. 8, 2014.

S. Zhang, J. Yanke, D. R. Johnson, and M. J. Krane, “Modeling defects in castings using a volume of fluid method,” International Journal of Numerical Methods for Heat & Fluid Flow, vol. 24, no. 2, 2014.

A. Konieczny and M. Ragulskis, “Analysis of simplifications applied in vibration damping modelling for a passive car shock absorber,” Shock and Vibration, vol. 2016, Article ID 6182847, 9 pages, 2016.

P. Wang, Q. Wang, X. Xu, N. Chen, and E. Petrov, “Fractional critical damping theory and its application in active suspension control,” Shock and Vibration, vol. 2017, Article ID 2738976, 8 pages, 2017.

T. W. Simpson, J. D. Poplinski, P. N. Koch, and J. K. Allen, “Metamodels for computer-based engineering design: survey and recommendations,” Engineering with Computers, vol. 17, no. 2, pp. 129–150, 2001.

M. Feinstein, B. LaFrence, K. Safford, and C. Timmons, Hydraulic Fluid Dissolved Gas Tester, 2009.