Numerical Study of Tire Hydroplaning Based on Power Spectrum of Asphalt Pavement and Kinetic Friction Coefficient

Hydroplaning is a driving phenomenon threating vehicle’s control stability and safety. It happens when tire rolls on wet pavement with high speed that hydrodynamic force uplifts the tire. Accurate numerical simulation to reveal the mechanism of hydroplaning and evaluate the function of relevant factors in this process is significant. In order to describe the friction behaviors of tirepavement interaction, kinetic friction coefficient curve of tire rubber and asphalt pavement was obtained by combining pavement surface power spectrum and complex modulus of tread rubber through Persson’s friction theory. Finite element model of tirefluid-pavement was established in ABAQUS, which was composed of a 225-40-R18 radial tire and three types of asphalt pavement covered with water film. Mechanical responses and physical behaviors of tire-pavement interaction were observed and compared with NASA equation to validate the applicability and accuracy of this model. Then contact force at tire-pavement interface and critical hydroplaning speed influenced by tire inflation pressure, water film thickness, and pavement types were investigated. The results show higher tire inflation pressure, thinner water film, and more abundant macrotexture enhancing hydroplaning speed. The results could be applied to predict hydroplaning speed on different asphalt pavement and improve pavement skid resistance design.


Introduction
Skid resistance of asphalt pavement on rainy days is an essential element for improving highway safety. When tire rolls on flooded asphalt pavement, contact force at tirepavement interface causes water flowing and hydrodynamic pressure, which would destroy previous vertical balance between loading and bearing capacity from pavement. Once a critical speed is obtained, tire floats entirely on the water and vehicle instability happens. The phenomenon is called hydroplaning and the critical speed is called hydroplaning speed. The loss of skid resistance on rainy days leads to the decrease of hydroplaning speed and friction coefficient, which concerns tire characteristics and asphalt pavement topography simultaneously.
In 1976, AASHTO published Guidelines for Skid-Resistant Pavement Design to improve asphalt pavement skid resistance through mixture design [1]. However, actual pavement friction is related to pavement-tire interaction, where tire material characteristics is significant. In 2001, Persson developed friction theory based on energy dissipation, which could be applied to depict tire-pavement friction behaviors [2]. Considering the effect of water on skid resistance, NASA conducted field test and obtained the famous empirical hydroplaning equation, where the critical hydroplaning speed was proportional to square root of tire inner pressure [3]. Wies investigated the impact of tire tread on hydroplaning speed by laboratory test, which contained the effect of tread area, direction, and position [4]. J.-Y. Jeong and H.-Y. Jeong analyzed the hydroplaning phenomenon when water film was less than 5 mm by finite element software LS-DYNA, and the power exponential function between fluid lifting force and traction force was regressed [5]. Sapragonas applied fluid-solid coupling method in PATRAN and DYTRAN to construct hydroplaning FE model, which could obtain critical hydroplaning speed [6]. Cheng et al. built hydroplaning model with complex tire tread and analyzed the relationship between tire-pavement interaction force and rolling velocity; hydroplaning speed prediction formula was then provided [7]. Fwa and Ong compared the hydroplaning 2 Advances in Materials Science and Engineering speed and skid number on pavement with longitudinal and transversal grooves and revealed the significance of using transversal grooves to reduce hydroplaning probability [8,9]. Srirangam et al. calculated friction coefficient of different pavement using FE hydroplaning model and verified it through comparison with field test data [10]. Macrotexture was considered important in improving skid resistance on flooded road. However, most FE models above did not realize the modelling of real asphalt pavement and failed in considering the fact that the kinetic friction coefficient varied with different vehicle speed. This could lead to the inaccuracy of numerical results. In this paper, pavement power spectrum of three typical types of asphalt pavement was calculated, and kinetic friction coefficient of rubber-pavement was deduced through Persson's friction theory. Pavement model with topography and tire model with tread pattern were established in ABAQUS; then water flow model was applied to simulate tire hydroplaning phenomenon. The FE hydroplaning model was verified by observing the contact force and water trace flowing in the tire pattern groove and comparing the hydroplaning speed with NASA hydroplaning equation. Then the influence of relevant factors in terms of tire inflation pressure, water film thickness, and asphalt pavement types on hydroplaning speed was analyzed.

Friction Model of Rubber and Asphalt Pavement
Previous FE model of hydroplaning mostly adopted constant friction coefficient derived from pavement surface characteristics. This value could not reflect the real contact mechanics of tire-pavement interaction and led to the inaccuracy of FE model. In order to overcome this weakness, pavement power spectrum and complex modulus of tire rubber were measured, and Persson's friction theory [2,11,12] considering energy dissipation was introduced to deduce the kinetic friction coefficient. Based on Chinese specification [13], slab specimens were made in the laboratory with length and width of 300 mm and height of 50 mm. According to the designed gradations shown in Table 1, the asphalt mixtures are asphalt concrete (AC), stone mastic asphalt (SMA), and open graded friction course (OGFC), which are widely used as surface layer in asphalt pavement.
Three-dimensional optical scanner (3DS) was used to acquire the three-dimensional topography of the slab specimens. In the process of measurement, several probes in 3DS worked together to collect the distance information. The coordinates of the surface points were obtained based on triangle measurement principle. Then the coordinate data were processed in MATLAB to display the 3D topography of different asphalt pavement in Figure 1.
Previous studies demonstrated the fractal features of asphalt pavement [14,15]. Pavement surface was self-affine when the wave length of asperities was between 1 and 0 , where 1 was the shortest distance cutoff and 0 was the upper distance cutoff which corresponded to the largest grain sizes in the mixture. Power spectrum of surface roughness ( ) was introduced to represent pavement surface characteristics, and = 2 / .
where ℎ( ) was the surface height measured from average surface plane. ⟨⋅ ⋅ ⋅ ⟩ indicates ensemble averaging. The power spectrums of asphalt pavement were calculated as shown in Figure 2.
Friction theory was introduced to describe the kinetic friction of three different types asphalt pavement. The surface power spectrum was divided into two regions approximately.
The dynamic elastic modulus at several frequencies and temperatures of a carbon and silica reinforced rubber compound used for the tire was measured using dynamic shear rheometer. Then a generalized Maxwell model with n terms was applied to describe the complex modulus of the rubber. Williams-Landel-Ferry (WLF) function was used and the master curve of viscoelastic modulus was obtained in Figure 3. The imaginary part of viscoelastic modulus depicted the energy dissipation when rubber is sliding on the oscillating surface asperities. The friction coefficient was calculated as

Advances in Materials Science and Engineering
where V is the sliding velocity, is the Poisson's ratio, and 0 is the nominal stress. ( 0 V cos ) is the viscoelastic modulus of rubber at the frequency = 0 V cos , where is the scaling factor in fractal and is the polar coordinate. The relation between sliding velocity and friction coefficient on three different asphalt pavement was illustrated in Figure 4. Compiling the FRIC subprogram in ABAQUS, the relationship was applied to monitor the real-time velocity and adopted the corresponding kinetic friction coefficient.

Asphalt Pavement
Model. Before the model construction work in ABAQUS, industrial computerized tomography and digital image processing technique were used to acquire the three-dimensional space structure of asphalt pavement.
In the first place, slab specimens were scanned by the X-ray tomography. 500 section images were acquired with a 4 Advances in Materials Science and Engineering scanning spacing of 0.1 mm, which met the accuracy requirements of macroscopic texture. After the section images were noise reduced, enhanced, and binary processed with MAT-LAB, voids in asphalt mixtures were completely removed in this process. Eventually, FE model of asphalt pavement with different macroscopic texture features was established in ABAQUS. These procedures are shown in Figure 5.

Tire Model.
In order to describe the characteristics of real tire, 225-40-R18 radial tire in this research was modelled as a composite structure. Rubber was considered as superelastic material and reinforcement was viewed as linear elastic material. The geometrical construction and meshing work were carried out as follows.
The procedures of constructing a 3D tire model with tread starts with a 2D half-section model. Then a 2D fullsection tire model with mesh was obtained by mapping the model of central symmetry (Figure 6(a)). The 3D tire model was obtained by rotating the 2D full-section tire model 360 degrees around the -axis ( Figure 6(b)). As a result, a slick tire with mesh was constructed. After the tire tread model was generated by rotating a single tread component (Figure 6(c)), these two models were combined by applying a tie constraint ( Figure 6(d)). Rubber material in the tire model was meshed as tetrahedron and hexahedron elements, the types of which were C3D8R and C3D6 in ABAQUS. The mesh size of tread part and sidewall part was 5 mm × 6 mm × 2 mm and 2 mm × 3 mm × 2 mm. Belt part was composed of SFM3D4R elements with the mesh size of 3 mm × 2 mm.
Belt was viewed as linear elastic metal material, and rubber materials were assumed to be homogeneous, isotropic, and incompressible. Considering the large deformation behavior of these materials during tire rolling on the pavement, superelastic model instead of linear elastic model was used to describe the mechanical properties of these materials. Excluding the effect of thermal expansion, the equation of unit volumetric energy was given by the principle of virtual work and fitted by polynomial, which could be expressed as follows: where and were the material constants, was the unit volumetric energy, was the volumetric variation, and 1 and 2 were the first and second invariant of the Cauchy-Green deformation tensor.
In general, there were two models widely applied for characterizing the mechanical properties of tire rubber materials with the assumption of = 1 or = 2, namely, Mooney-Rivlin model and Yeoh model [16].
Uniaxial tensile test was then carried out using universal testing machine on rubber specimens in the laboratory. Relationship of nominal stress and strain was acquired by transforming the axial force and deformation in the loading direction. Afterwards, the experimental data was fitted respectively in Mooney-Rivlin model and Yeoh model. The fitting curves of carcass and tread were shown in Figure 7.
Yeoh model in Figure 7 showed superiority in describing the mechanical properties of tire rubber materials, which were chosen in tire model in this paper. The parameters selected for Yeoh model in the research are shown in Table 2.

Application of Coupled Eulerian-Lagrangian Method.
Water film in the numerical hydroplaning model was approximated as Newtonian fluid, which obeyed the conservation of mass, momentum, and energy. Equation of state was utilized to represent all potential equilibrium states of fluid, which confirmed the relation among fluid pressure, density, and unit mass internal energy as follows [17].
Advances in Materials Science and Engineering 5  where was fluid pressure, was fluid density, was unit mass internal energy, and ( ) and ( ) were parameters in equation of state.
Mie-Grüneisen model [18] was applied in this paper, which is where 0 is initial fluid density; is nominal volumetric strain; 0 is propagation velocity of sound in pure water; is the model parameter, which is 1.979; Γ 0 is material constant, which is 0.11. Euler elements and Lagrange elements were used to model tire and fluid, respectively, which was combined by using the CEL method. In the model building process, the Lagrange meshes of tire were overlapped by the empty meshes in Eulerian element, which indicated that the tire immerged in part of the empty meshes. This step facilitated the loading of surface fluid pressure on the tire contacting region by tracking the motion of fluid free surface in calculating process. Simultaneously, the volume fraction and velocity boundary condition were updated according to the renewal of tire nodes displacement and velocities on the interaction surface using the volume of fluid (VOF) method. In this method, Eulerian volume fraction was defined on the fluid passing element meshes as follows, which was the function of space and time. where was the fluid velocity and was the time. When = 1, the mesh was filled with fluid; when = 0, it was empty; when 0 < < 1, free surface existed in the mesh. Figure 8 shows the mesh partition situation of the fluid model. The white and red areas represent, respectively, water region and air region. With the model size of 80 mm × 390 mm × 320 mm, the mesh density is increased in the potential tire-pavement interface, and the mesh sizes decrease from 17 mm × 10 mm × 10 mm to 1 mm × 1 mm × 1 mm. Total number of EC3D8R Euler elements in the fluid is 381420.

Description of the Numerical Model.
After the tire model with tread pattern was constructed and imported into ABAQUS, implicit analysis was carried out to apply tire inflation pressure and vehicle load. Then the numerical results of tire deformation and stress distribution were transferred in explicit analysis step to conduct hydroplaning simulation. The FE model for simulating the hydroplaning phenomenon X Y Z Figure 9: Tire hydroplaning FE model. was shown in Figure 9. There were three stages involved in this procedure.
A Uniformly distributed load was applied on the inner surface of tire model, which was viewed as initial tire pressure in the following discussions. Note that this pressure was completely incompatible with tire-road contact pressure in pavement mechanics.
B Then a small displacement in the negative direction of -axis was applied on the pavement model with the rim of the tire model fixed. After the tire model converged, quarter of a vehicle weight was applied on the pavement in the same direction as the vehicle load. The value in this research was set at 3822 N.
C The 3D tire model rolling on the asphalt pavement with water film was realized with an accelerated motion stage on dry pavement and a uniform motion stage transiting onto the water surface. In this process, an angular velocity ( ) around -axis and equivalent linear velocity (V = × ) in the negative direction of -axis were specified for tire model and asphalt pavement model with water film. In this way, tedious work for establishing fluid model could be reduced.

Verification of the Tire Hydroplaning Model.
In order to validate the applicability of the FE model in simulating the hydroplaning phenomenon, contact force in tire-pavement interface and water trace flowing in the tire pattern groove are observed and analyzed. Then the hydroplaning speeds with different tire inflation pressure obtained in the FE model are compared with NASA equation to verify the accuracy of the hydroplaning model. In this process, numerical simulations of tire hydroplaning phenomenon are carried out with constant load of 3822 N and the initial tire inflation pressure of 250 kpa.
The FE model provides a useful tool to monitor the mechanical responses and physical behaviors of tirepavement interaction. The contact force oscillates around a constant value on dry pavement; then it starts to decrease due to the beginning of tire-water film interaction ( Figure 10). The extrusion of water into the groove of tire tread pattern with velocity of 70 km/h is shown in Figure 11. Eventually, the normal contact force in tire-pavement interface declines to a constant value, which determines the occurrence of hydroplaning.
The hydroplaning speed obtained in the simulation and calculated in NASA equation is compared. In this process, the hydroplaning speed is determined with variable tire inflation pressure and fixed tire load and water film thickness of 4800 N and 7.62 mm. The results are shown in Figure 12.
With the tire inflation pressure of 100 kpa, the hydroplaning speed is 70 km/h in the simulation and 63.6 km/h in the NASA equation; with the tire inflation pressure of 250 kpa, the hydroplaning speed is 107.6 km/h in the simulation and 100.56 km/h in the NASA equation. The simulation deviations are 10.01% and 7%, respectively. The deviation could be attributed to tire tread pattern and pavement topography. However, the tendency confirms the accuracy of the tire hydroplaning model established in the research.

Results and Discussions
In this section, hydroplaning phenomenon under different conditions is simulated in the FE model described above, and different relevant factors are evaluated on the influences on hydroplaning characteristics. Contact force at tire-pavement interface and critical hydroplaning speed in these situations are calculated and compared. The loading in this part is set at 3822 N.

Effect of Tire Inflation Pressure on Hydroplaning.
The contact force and critical hydroplaning speed on AC pavement are illustrated in Figure 13. It is shown that contact force decreases faster for higher inflation tire in partly hydroplaning. However, once the tire entirely sinks into water film, tire with lower inflation pressure is matched with smaller contact force (Figure 13(a)). Figure 13(b) presents the variation of critical hydroplaning speed when tire inflation pressure changes from 100 kpa to 250 kpa. It is shown that the hydroplaning speed increases almost 50% linearly during the variation of tire inflation pressure for each water film thickness. According to the mechanism analysis in the FE model, the increase of moving speed is equivalent to the increase of tire inflation pressure in offering a higher lifting force at tire-pavement interface.

Effect of Water Film Thickness on Hydroplaning.
Vertical contact force at tire-pavement interface with different water film thickness is calculated. Figure 14(a) showed that the force decreases faster for thicker water film. When the tire entirely sinks into water film, ultimate contact force is higher for thinner water film. Figure 14 thickness increases from 0.5 mm to 10 mm. This could be attributed to higher lift force offered by thicker water film.

Effect of Asphalt Pavement Types on Hydroplaning.
In this section, the influence of different asphalt pavement on hydroplaning characteristics is evaluated. When tire rolls on dry pavement, the value of contact force at tire-pavement interface fluctuates around a fixed value, which is the loading value (3822 N). From rolling on the dry pavement to entirely sinking into water film, the contact force on OGFC pavement is the largest, followed by SMA and AC. Figure 15(a) reveals the superiority of OGFC pavement in increasing the contact force at tire-pavement interface.
Then the hydroplaning speeds are calculated with different asphalt pavement types. Figures 15(b) and 15(c) illustrate the similar variation of hydroplaning speed in three types of asphalt pavement. With thinner water film and higher tire inflation pressure, hydroplaning speed is increased. Compared with other two asphalt pavement, OGFC owns better performance in increasing hydroplaning speed, which could be attributed to higher roughness and more drainage channels.

Conclusions
In this paper, kinetic friction coefficient of rubber-pavement was deduced through Persson's friction theory combining pavement power spectrum and viscoelastic property of tire tread rubber material. Then tire-fluid-pavement hydroplaning model was established in ABAQUS. After the acquired hydroplaning speed was verified through comparison with NASA equation, the influences of relevant factors in hydroplaning were evaluated.
(1) Comparing hydroplaning speed acquired from the proposed model and NASA empirical equation, the applicability and accuracy of this numerical hydroplaning model were verified.
(2) Initial tire inflation pressure has a positive influence on improving skid resistance of flooded pavement.
The contact force at tire-pavement interface and critical hydroplaning speed are enhanced with the bigger tire inflation pressure.
(3) Thicker water film on the pavement leads to higher lift force on the tire, which would lead to smaller contact force and critical hydroplaning speed, which was unfavorable for vehicle safety.
(4) The types of asphalt pavement had an obvious impact on hydroplaning characteristics. OGFC owns better performance in increasing skid resistance, which could be attributed to higher roughness and more drainage channels.

Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.