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

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. )erefore, it is critical to improve the mathematical modeling accuracy to enhance the design and control efficiency and accuracy of HPS. )is 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.)e 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.)e corresponding damping and stiffness force of each sample point were calculated by CFD simulation.)e obtained simulation data were utilized for the fitting of a Kriging model. )e results demonstrated that the Kriging model can provide high accuracy, with a prediction error lower than 5%. )e proposed modeling method of the HPS nonlinear characteristics is highly efficient, accurate, and faster than traditional methods.


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][3][4]. erefore, 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. ermodynamic 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. eir 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]. e 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. e 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. e 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]. e stiffness is affected by the gas volume and the area of the piston and piston rod [15]. e 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][19][20][21][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. e 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 gasoil 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. erefore, it is of particular importance to propose a highprecision 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.
e 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.
is 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. e 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.

Structure and Working Principle.
e 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). e 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. e 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. e parameters of the HPS are presented in Table 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 P 0 and V 0 are the initial pressure and volume of the HPS, respectively, P and V are the instantaneous gas pressure and volume during HPS operation, respectively, and r is the polytropic constant. In general, when the change of the gas state is regarded as an isothermal process, r � 1, while when it is regarded as an adiabatic process, r � 1.4. e compression force in the cylinder can be calculated as follows: where A 0 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: en, the stiffness of the hydropneumatic suspension can be calculated as follows: where x 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, dA 0 dX � 1, Consequently, the stiffness of the single-chamber HPS can be obtained as follows: where A is the cross-sectional area of cavity (II). When the hydraulic oil flows through the damping valve, the flow is purely turbulent.
us, 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, C z and C d are the flow coefficients of the damping and check valves, respectively, S z and S d are the effective overcurrent areas of the damping and check valves, respectively, and sign( _ x) is the sign function. When the suspension is in an extension stroke, sign( _ x) � −1; when the suspension is in a compression stroke, sign( _ x) � 1. e pressure drop Δp at both sides of the damping valve can be calculated as follows: where v is the average oil flow velocity, l is the length of the damping valve, and d is the diameter of the damping valve. When the Reynolds number Re � vd/υ, q � πd 2 v/4, 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:

Two-Phase Flow Simulation Method for HPS
e 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. erefore, the VOF method was employed for the simulations in this study.

VOF Fluid Phase Interface Tracking Method.
e 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. e 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 q th fluid in the cell is denoted as α q , then the following three conditions are possible: α q � 0: the cell is empty (of that fluid) α q � 1: the cell is full (of that fluid) 0 < α q < 1: the cell contains one or more interfaces between that fluid and one or more other fluids e governing equations of the VOF method include the volume fraction equation, the additional scalar equations, and the continuum surface force model.

Volume Fraction Equation.
e 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 q th fluid, this equation takes the following form: where _ m qp is the mass migrated from phase q to phase p, _ m pq is the mass migrated from phase p to phase q, S α q is the source term, ρ q is the density of phase q, and v → q is the velocity of phase q. e volume fraction equation is not solved for the primary phase; the primary-phase volume fraction is computed based on the following constraint:

Additional Scalar Equations.
e density and viscosity of the mixed fluid in the unit are calculated as follows: where ρ g is the density of the gas, ρ l is the density of the liquid, μ g is the viscosity of the gas, and μ l is the viscosity of the liquid.

Continuum Surface Force Model.
e 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, R 1 and R 2 : where p f1 and p f2 are the pressures of the two fluids on either side of the interface. e 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 σ ij is the surface tension coefficient between phase i and phase j, α i and ρ i are the volume fraction and density of phase i (the primary phase is usually liquid), respectively, α j and ρ j are the volume fraction and density of phase j (the secondary phase is usually gas), respectively, and κ i and κ j 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:

Nitrogen Dissolution in Oil.
e 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.
e 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 V gi is the initial volume of nitrogen in the HPS and V g is the volume of nitrogen dissolved in hydraulic oil. e density of the gas-oil mixture ρ mix can be calculated as follows: where M o is the total mass of oil in the HPS and V oi and ρ oi 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.

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

3.3.2.
Meshing. e 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.

Numerical Calculation
e 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. e oil used in the HPS was L-HM46, which can be regarded as Newtonian fluid. us, the oil was modeled as incompressible Newtonian fluid. e dynamic mesh method was used to simulate the relative movement between cylinder and piston rod. e moving boundaries were the cavity (I) upper wall and the cavity (II) bottom wall, as shown in Figure 5(a). e boundary movement was prescribed using a profile file in the Fluent software. e 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. e calculation settings for the transient simulation are listed in Table 3. . 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 m/s at t � 0.05 s, the maximum oil flow velocity was 42 m/s. However, at t � 0.1 s, the maximum fluid velocity reached 80 m/s, while the vibration speed of the suspension reached 1 m/s, 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. is is due to that the flow rate through the    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. is 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 MPa, and the pressure difference between cavity (I) and cavity (II) was 1.2 MPa. At t � 0.1 s, the maximum pressure reached 12.5 MPa, and the pressure difference between the two cavities was 4 MPa. It can be seen that the pressure difference between the two cavities increased with the increase of the suspension vibration speed. e 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. erefore, 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).  First-order implicit e velocity distribution contours of flow at a crosssection 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. e maximum velocity at t � 0.05 s and t � 0.1 s was 95 m/s and 190 m/s, 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 m/s to 1 m/s, 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.

Simulation Results and Flow Field
e 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 m/s and the maximum pressure was 16 MPa, which was much higher than that in the compression stroke (7.2 MPa). When the suspension vibration speed increased to 1 m/s at t � 0.1 s, the maximum pressure reached 50 MPa. As the suspension vibration speed increased from 0.5 m/s (t � 0.05 s) to 1 m/s (t � 0.1 s), the pressure difference between cavity (I) and cavity (II) increased from 10 MPa to 45 MPa, respectively. 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 P 1 and P 2 are the average total pressures of cavity (I) and cavity (II), respectively, and A 1 and A 2 are the crosssectional areas of cavity (I) and cavity (II), respectively.
Since the pressure difference ΔP between cavity (I) and cavity (II) is then equation (21) can be written as e output force of the HPS is divided into two parts, stiffness force F k and damping force F c , 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.

Surrogate
Modeling. e specific implementation scheme of the mathematical modeling process for HPS is presented in Figure 12. e 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.

Kriging Model.
e 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). e Kriging model can be given as follows [17]:

Optimal Latin Hypercube Experiment Design.
e 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 (V) on the HPS nonlinear characteristics were taken into consideration. e 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. us, it is necessary to investigate the effect of oil viscosity variation with temperature on the HPS damping and stiffness characteristics. e variation range of each factor was determined empirically (Table 4). e temperature range is between 25 and 65°C. e hydraulic oil viscosity, which varies over the temperature, ranges from 0.0159 to 0.1252 Pa · s. 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. erefore, in this study, the variation Shock and Vibration 11 range of gas dissolution rate was determined to be 0% to 10%. e vibration speed, which ranges from −1 m/s to 1 m/s, takes negative values during compression stroke and positive during extension stroke. e 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. e sample points were then utilized to perform simulation experiments using the two-phase flow simulation method mentioned in Section 2. e sample points and the corresponding simulation results are listed in Table 5.

Model Fitting Evaluation.
e 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. e root-mean-square error can be calculated as follows: where y obs,i and y model,i 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 y i is the sample estimates and y 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. e 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.

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. erefore, 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. e prediction accuracy of the Kriging model on damping and stiffness force is higher than that of the response surface model. erefore, 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. erefore, the method proposed in this paper to model the HPS nonlinear characteristics considering multiple influencing factors is validated.

Surrogate Model of HPS Nonlinear Characteristics.
e Kriging model cannot be shown as an explicit formulation. e 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 (V) of the suspension. 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. is difference is attributed to the changes in the physical properties of the hydraulic oil caused by gas-oil emulsion. Figure 14

Factors
Variation range Temperature (°C) 25∼65 Oil viscosity (Pa · s) 0.0159∼0.1252 Gas dissolution rate (%) 0∼10 Vibration speed (m/s) −1.0∼1.0    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. e effect of the interaction between suspension vibration speed (V) and T on damping force is illustrated in Figure 14(c). In this case, compared to T, V plays a dominant role in affecting the damping force of the HPS, and the function curve  between V and F c on the curved surface is nearly parallel, independent of the value of T. is means that the interaction between V 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 V on damping force is more apparent, while the interaction between V 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 V and D on stiffness force. It can be observed that V has a significant impact on stiffness force, which, in compression stroke, increases with the increase of V, while in extension stroke, the increase of V leads to a decrease in stiffness force. e effect of the interaction between V and T on stiffness force is presented in Figure 15(c). It can be clearly observed that, compared to V, the variation of T plays a secondary role in affecting the stiffness force. In terms of the degree of parallelism of the F k -V and F k -T function curves on the surface, no interaction between V and T can be identified. Figure 15(d) reflects the effect of the interaction between V and U on stiffness force. rough analysis, it can be observed that the effect of V on stiffness force is more significant than that of U. e interaction between U and V 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.  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. 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. erefore, it can be concluded that the viscosity and temperature effects of hydraulic oil have a significant impact on the damping force of HPS. e variations in suspension vibration speed have the greatest effect on stiffness force; next are the temperature and oil viscosity. e influence of gas dissolution rate on the stiffness force is not apparent.

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. e 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.

Experimental Procedure.
e driving experiment using the mining dump truck was performed in an open pit mine in western China (Figure 17). e 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. e truck speed ranged between 10 and 30 km/h, and the input was random road. e experimental results concerning the suspension upper point at a speed of 30 km/h 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. e root-mean-square values of the acceleration of the front left, front right, rear left, and rear right suspensions were 2.2 m/s 2 , 2.2 m/s 2 , 3.1 m/s 2 , and 3.3 m/s 2 , respectively.

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:     € θ � m s gh 2 θ − aF fl − aF fr + bF rl + F rr , where m s is the body mass, z s is the vertical displacement of body mass center, I xxs is the pitch moment of inertia, I xxy is the roll moment of inertia, h 2 and h 1 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 △z sij is the suspension deflection.
e 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. e accuracy of the proposed model was validated by the indirect verification method.

Conclusions
e 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. e oil temperature, oil viscosity, gas dissolution rate, and suspension vibration speed were investigated as the influencing factors of the nonlinear characteristics of HPS. e findings of this study are summarized as follows: (1) e 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. is has very important guiding significance for the analysis and development of HPSs.
(2) e Kriging model demonstrated excellent performance for approximating the mathematical relationship between HPS nonlinear characteristics and influencing factors. e 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
e authors declare that they have no conflicts of interest.