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.
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.
The parameters of the HPS are presented in Table 1.
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 aswhere 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.
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.
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.
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)
The calculation settings for the transient simulation are listed in Table 3.
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)
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 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)
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)
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)
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)
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) isthen 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.
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.
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.
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.
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.
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)
(a)
(b)
(c)
(d)
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)
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.
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)
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)).
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)
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.