Guidance to Predict and Mitigate Dynamic Hydroplaning on Roadways (2021)

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

27 3. INTEGRATED HYDROPLANING MODEL Models for tire dynamics, vehicle dynamics, fluid dynamics, water accumulation, and driver behavior have traditionally been developed as separate entities for different applications. For this project, we integrated the models specifically for hydroplaning prediction. Many features that may have been considered unimportant in other applications became critical to predicting hydroplaning. This research developed a research-grade Integrated Hydroplaning Model composed of water accumulation, tire-water-pavement, and vehicle dynamics submodels, developed as separate entities for different applications. The most appropriate modeling approaches and tools specifically for hydroplaning prediction were integrated. Although the original intent was to develop a fully coupled model that could run the water accumulation, tire, and vehicle response simultaneously, this proved impractical and thus, the final model includes a series of submodels that run one after the other. The maximum WFT is first calculated along the pavement and is then taken as a constant during the vehicle performance simulation. Similarly, the interaction between the tire, pavement surface, and water accumulated on the pavement is modeled and the resulting forces are used as an input to the vehicle dynamics model. To simplify the use of the model without the need of commercial software, the research also developed simplified relationships from the results of the research-grade hydroplaning model and implemented these relationships in a Hydroplaning Potential Assessment Tool that is discussed in the next chapter. 3.1. MODEL ARCHITECTURE The first step in the integration process included understanding the required fidelity of the models and choosing existing model architectures that possess the appropriate level of fidelity. The interaction between the models is shown schematically in Figure 3. Note that for highway engineering applications, the simplified model provides a direct relationship between the road characteristics and PM that is used as a measure of HP. The rainfall intensities are given for a particular location, and the other characteristics can be fixed for a particular analysis (e.g., a mid- size car with a standard tire with new tread cruising along the road segment at a constant speed). Since WFT develops at a significantly different time scale than the vehicle dynamics, it can be treated as a constant when considering vehicle response. That is, the WFT is simulated at different depths but taken as a constant during the vehicle performance simulation. For example, the rate at which rainfall occurs is cited in inches per hour. Presumably, the time constant for the dynamics of the WFT is on the order of minutes, while the vehicle suspension dynamics occur on the order of 5 ms to 2 s; thus, the dynamics of the WFT are well separated from the dynamics of the vehicle suspension. Similarly, the movement of the tire is considered when determining the forces created by the tire, whereas the inertia of the tire is lumped with the unsprung mass. In this way, the vertical (normal) force generated by the tire is considered a function of the spindle position and velocity along with the location of the road surface, and the dynamics of the tire occur more quickly than those of the vehicle suspension.

28 Figure 3. Flow of information between the vehicle response submodels. Consider a vehicle dynamics simulation in which the spindle position and velocity are calculated at a 5 ms time step. The road surface at this moment in time is known, and the spindle and road information are passed to the tire model. The tire model then calculates the force at that moment in time and returns that information to the vehicle model. Implicit in this modeling configuration is the fact that the dynamics of the tire occur much more quickly than 5 ms. That is, the local deformation at the tire surface due to the water on the road occurs on time scales less than 5 ms, the vehicle suspension responds on time scales of 5 ms to 2 sec, and the WFT occurs on the order of minutes. Finally, the dynamics of the water at the tire/pavement interface are even faster than the tire dynamics. That is, the tire position at any moment in the tire-water-pavement (fluid dynamics) model is considered a constant, and the water responds dynamically to the tire configuration. The tire shape and location will change at the next time step in the fluid dynamics simulation (perhaps 100 µs later), but the tire configuration is considered a constant at any particular moment. The speed of the vehicle is also considered as an input for the vehicle response model as it is controlled by the driver, who often adjusts speed based on the maneuver and weather conditions. 3.2. WATER ACCUMULATION SUBMODEL To get away from simplistic “rules” and increase our understanding of the physics of the water flow, as determined by rainfall rate, road geometry, and pavement texture, CFD can be used to simulate the formation of the water film. A 3-D CFD model was developed to predict the accumulation of water on multilane highways using ANSYS Fluent software. The RANS equations representing the conservation of mass and momentum are solved in the CFD model.

29 The two-fluid Eulerian model was employed, and the water and air phases are treated mathematically as interpenetrating continua, where the concept of phasic volume fraction is introduced. The Eulerian model is the most complex of the multiphase models and solves a set of momentum and continuity equations for each phase. The continuity equations governing the mass balance of each phase and the phase-averaged Reynolds-averaged momentum equation, respectively are: 𝜙𝜌 + ∇⃑ ∙ 𝜙𝜌𝑢 = 0 (11) 𝜙𝜌𝑢 + 𝑢 ∙ ∇⃑ 𝜙𝜌𝑢 = −𝜙∇⃑𝑃 + ∇⃑ ∙ 𝜏̿ + 𝜌 − 𝜌 ?⃑? (12) where: ρ = fluid density, ϕ = volume fraction scalar, 𝑢 = field velocity, 𝜏̿ = deviatoric fluid stress tensor, ?⃑? = gravity field, and P = pressure. The terms on the right-hand side in equation 12 represent, from left to right, the pressure gradient, the shear stress, the gravitational force, and the surface tension body force. This is an improvement with respect to previous work with 1-D analytical models or kinematic wave models, which assumed that gravity was the only driving force (Reed and Kibler 1983; Reed et al. 1989; Huebner et al. 1997; Anderson et al. 1998; Jeong 2008; Jeong and Charbeneau 2010). Volume fraction ϕ represents the space occupied by each phase. Mass and momentum conservations are specified by each phase individually, where the water and air volume fractions sum to one: 𝜙 + 𝜙 = 1 (13) where φw = volume fraction of water, and φa = volume fraction of air. Equations 11 to 13 for water and air phases were solved in ANSYS Fluent 16.0 using a segregated double precision cell-centered solver. Details for discretization schemes and the combinations for Eulerian-Eulerian two-fluids model can be found in ANSYS (2011). The segregated pressure- based Navier-Stokes (PBNS) solver is used, which applies the projection method, where conservation of mass is rewritten in terms of the primitive variable pressure (ANSYS 2011). The PBNS solver is an effective way to solve incompressible flow, which is applicable for the flow in current study since local speed is subsonic and the pressure change in the system is small. The semi-implicit method for pressure linked equations (Simple) is applied for pressure-velocity coupling. For the gradient and spatial discretization, least-squares cell-based and the QUICK scheme, which is formally higher than second-order accuracy, are used for gradients and momentum terms. The modified implicit “high resolution interface capturing” scheme is applied to capture and interpolate the air-water interface. The convergence criteria for the continuity, momentum, and volume fraction of each phase ensure that the normalized residuals reduce to 10-

30 6. Time step size Δt has been decided based on the Courant-Friedrichs-Levy (CFL) number of 0.1, which can be expressed as: U tCFL x Δ= Δ (14) where: U = maximum fluid velocity and ∆x = average minimum cell size. The vertical position and thickness of the water film layer are determined from the simulated density profiles and fluid pressure acting on the pavement. The water film depth is calculated as: 𝑝 = 𝜌 − 𝜌 𝑔ℎ + 𝜌 𝑔𝐻 (15) where: ρw = density of water, ρa = density of air, p = the relative pressure, H = total height of the domain, and h = water film depth. Under constant rainfall conditions, the WFT distribution on a straight and real road can reach a quasi-equilibrium state when the mass from the rainfall is balanced with the water discharge out of the system and the water accumulates on the pavement. The quasi-equilibrium state is known when the instantaneous results of WFT fluctuate in a stable range and mean (time-averaged) and root mean squared (RMS) results become constant values. Reed and Kibler (1983) and Reed et al. (1989) defined the equilibrium flow time, teq, on an impervious surface under constant rainfall excess rate q: 𝑡 = / (16) where: L = length along flow path and width of runway with the unit of ft, and α and m = parameters of the hydraulic friction relationship. However, such equilibrium flow time only applies on a 1-D flow path with one longitudinal slope. In the current work, the equilibrium flow time and quasi-steady behavior are established using Equation 16 and obtained by dynamically monitoring the mass balance and the mean and RMS error of the WFT results. 3.2.1. Rainfall Intensity The major source of precipitation data is the National Weather Service (NWS), an agency of the National Oceanic and Atmospheric Administration (NOAA). Historically, NWS publications have been the primary source for precipitation depth-duration-frequency data across the United States. The NOAA Atlas 14 – Precipitation Frequency Data Server (PFDS) (NOAA 2004) provides the

31 most up-to-date information for precipitation frequencies. The NOAA’s National Climatic Data Center contains the raw rainfall data. As an example, Figure 4 show the rainfall rates for different durations at a weather station in Tampa Florida. This location will experience a rain event of approximately 50 mm in 1 hr on average very year. According to past records, an event of 125 mm over 1 hr will only occur every 1000 years. Each agency will need to select an appropriate rainfall rate depending on the level of risk that it is willing to take. In the example the analysis can select a value of 50 mm knowing that it will be exceeded once a year, a higher value up to 125 mm for a much higher certainty that it will not be exceeded, or something in between. Figure 4. Rainfall rates versus duration for Tampa City, FL (from NOAA 2014). In addition, the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) site developed by NASA can be used for reference (NASA 2020). This database provides climatic information with enhancements as it assimilates space-based observations of aerosols and represent their interactions with other physical processes in the climate system (Collow et al. 2017). 3.2.2. Initial Mesh Definition and Verification Work In the validation work, the study of 1-D flow with a longitudinal slope by Reed et al. (1989) was adopted, which includes both analytical and experimental work of water runoff. The experimental work measured the WFT on a 9.144-m by 4.572-m concrete slab with constant longitudinal slope Sx of 1.5% and an average texture depth of 0.66 mm under rain intensities of 29.21 mm/h. The experimental work reported the WFT measurements as a set of transversal and longitudinal sampling points. The schematic of the experimental work (Reed et al. 1989) is shown in Figure 5, where a slight slope was added on the pavement and a uniform rainfall was specified on the top using a rain 0.1 1.0 10.0 100.0 1000.0 0.0 0.1 1.0 10.0 100.0 1000.0 R ai nf al l R at e ( m m /h r) Duration (hr) 1 year 10 years 50 years 100 years 1000 years

32 simulator. The computational domain has been marked as a dashed line, based on the geometric characteristics shown in Figure 5 and the 1-D flow behavior reported by Reed et al. (Reed et al. 1989). The rain rate at the top surface of the domain is specified with a water velocity of 29.21 mm/h. The boundary condition of atmospheric pressure is specified along the side surfaces to allow water to exit the system. The pavement represented as the bottom of the domain is modeled as a non-slip boundary. To model the effects of mean texture depth (MTD), the pavement surface is modeled with a non-zero but evenly distributed macrotexture profile using the roughness height Ks in the wall function. Figure 5. Schematic of experimental work by Reed et al. (1989). The schematic shown in Figure 6 depicts the pavement parameters. The roughness effects are represented by two roughness parameters; namely, roughness height Ks and roughness constant Cs need to be specified. Ks equals to zero corresponds to smooth walls. As shown in Figure 6, an “equivalent” macrotexture parameter can be specified by the roughness height, Ks. The roughness constant Cs is designed to specify the uniformity of the rough surface. However, trials by the authors with different Cs values have revealed that it is not sensitive to the water film value, as long as the roughness height is the same. The main influential factor of the pavement roughness is the roughness height, which is twice the MTD, as shown in Figure 6. The water film value is extracted from the average water depth level for the CFD simulations by subtracting the MTD or half the Ks value. Figure 6. A schematic of macrotexture and CFD parameters For this initial analysis, the roughness height of the pavement, which represents the amplitude of the macrotexture, has been specified with the value of 1.32 mm, which is the average texture profile of the pavement. A roughness height Ks of zero corresponds to a smooth wall. The roughness constant, Cs, has been specified with the value of 0.5, indicating a uniformly distributed texture.

33 Since the WFT is usually a very small value, the grid resolution plays a very crucial role in capturing the water film and details of the flow field. A grid resolution study was conducted to optimize the predictive performance and simulation computational time. The current study tested three grid resolutions with the cell size and numbers listed in Table 5. In the Z coordinate direction, the mesh spacing increases according to a geometric progression with a modest expansion ratio of 5% so that the grid lines are concentrated in the region near the pavement. Initial conditions were set to pure air with air volume fractions of 1, indicating a dry condition on the pavement before the rain comes. Several locations were sampled in the computational domain to monitor when a quasi-steady-state of water is achieved. Table 5. Details of the Grids Used in the Comparisons Mesh ΔX (mm) # cells ΔZmin (mm) ΔZmax (mm) # cells # Total Coarse 10 1000 1.5 3.8 20 20,000 Medium 6.7 1500 1 3.4 26 39,000 Fine 4.4 2250 0.67 3 32 72,000 The CFD predictions are compared with the experimental work and 1-D models developed by previous researchers, shown in Figure 7. The solid circles with error bars show the data measured by Reed et al. (1989). The CFD results are shown as hollow symbols and represent the predictions using the coarse, medium, and fine meshes. Figure 7. Validation work with Reed at al. (1989) and 1-D correlations. Various lines are used to include the 1-D equation by Gallaway et al. (1971; 1979), Reed and Kibler (1983), and Andersen et al. (1998). The WFT measured by Reed et al. (1989) at the end of the road length has a very large error bar and was hardly captured by the CFD or 1-D equations, which indicates that the last measured point (close to the collecting channel) may have a large bias due to the experimental operating conditions. The CFD predictions using the medium and fine

34 meshes show good agreement with the first three measured values by Reed et al. (1989). The CFD results (medium and fine meshes) also show very similar WFT profiles with Gallaway et al. (1971) and Gallaway et al. (1979), where the slopes of WFT increase along road length. 3.2.3. Model Validation Through Sensitivity Analysis The model was further validated using more complex geometries similar to those used by Charbeneau et al. (2008). This section presents the results of a thorough sensitivity analysis of WFT on superelevation transitions using the developed CFD model, including the impact of different weather conditions, road characteristics, and pavement surface properties. Specifically, factors like rainfall rate, longitudinal and cross slopes, road lanes, pavement texture depth, and road types have been tested to study the effects on WFT distributions. Superelevation is a roadway design feature used to control centrifugal forces as vehicles move through curves. Superelevation transitions are used to change the road transverse section from a normal crown to a superelevation and then back to a normal crown. The schematic of such transitions is shown as sections AB and DE in Figure 8 (upper left). Figure 8. Schematic of two superelevated transitions (upper left A-B and D-E) and the change of longitudinal and cross slopes in transition section AB (lower left and right). In transitions leading to curved alignments (shown as section AB in Figure 8), the inward travel lane edge is lowered in elevation while the outer travel lane edge is raised. A 0% cross slope occurs at the point where the outward side of the road passes from a negative cross slope to a positive cross slope. The change of the cross section slopes is shown in Figure 8 (right). Since the road also has a longitudinal slope (Sx %), shown as the lower left schematic in Figure 8, the combination of the cross and longitudinal slopes causes the road elevation to change along both the X and Y directions, which requires a 3-D analysis. The section DE represents the transition

35 after the curved alignment, where the cross slope changes from superelevation back to normal crown. Since the changes of longitudinal and cross slopes are opposite to section AB, the two sections are symmetric about point C. In this section, a thorough sensitivity analysis of WFT on superelevation transitions using the developed CFD model is conducted, including the impact of different weather conditions, road characteristics and pavement surface properties. Table 6 shows the cases with various combinations, including longitudinal slope, cross slope, macrotexture, road type, and rainfall rate. Table 6. Operating Conditions for 3-D Simulation of WFT Cases Longitudinal Slope Type Rainfall Rate (mm/h) Lanes Texture Depth (mm) Purpose 1 1% Type I 60 2 0.5 Base case 2 1% Type I 20 2 0.5 Rain fall rate 3 1% Type I 100 2 0.5 4 0% Type I 60 2 0.5 Slopes 5 2% Type I 60 2 0.5 6 1% Type I 60 2 0 Macrotexture 7 1% Type I 60 2 1 8 1% Type II 60 2 0.5 Type 9 1% Type I 60 4 0.5 Lanes 10 1% Type I 60 8 0.5 Case 1 is the base case to which other cases are compared. For example, cases 2 and 3 test the effects of rainfall rate from normal rain to heavy rain, cases 4 and 5 test the effects of longitudinal slope, and cases 6 and 7 are testing macrotexture effects. To compare different road transitions, Case 1 is the transition before superelevation and Case 8 tests the transition after superelevation. Cases 1, 9, and 10 compare different number of lanes. The simulations investigate the WFT distribution due to the small changes of reference height on a 50-meter long pavement surface with different road parameters and rainfall rates. As elaborated in the introduction and shown in Figure 8, the numerical simulation focuses on the superelevation transition, where various slopes are designed and distributed laterally and longitudinally. Figure 9 presents the geometry of the computational model based on the transition section AB of a real two-lane road. Based on the grid resolution study in the previous section, the cell size in the vertical direction has been kept at the medium resolution for all of the following simulations of real road conditions.

36 Figure 9. Schematic of 3-D computational domain for real road conditions. Effect of Rainfall Rate The effect of rain intensity was studied by changing it from 20 mm/h to 60 mm/h to 100 mm/h. Figure 10 shows the WFT distribution with different rainfall rates. It should be noted that the shoulders on each side of lanes have been excluded from the figures due to the specific interest of WFT on the road rather than the shoulder. When the outside lane changes its cross slope from a negative value to a positive value, water tends to move toward the inner edge downstream. (a) Case 2 with 20 mm/h (b) Case 1 with 60 mm/h (c) Case 3 with 100 mm/h Figure 10. WFT distribution on pavement with different rainfall rate (ZCS = zero cross slope; RR = rainfall rate).

37 It is clear in Figure 10 that the WFT results all show similar distribution patterns with changing rainfall rate, while the magnitude of WFT has greatly increased with increasing rainfall rate. Thus, WFT distribution patterns depend mostly on the road characteristics rather than on rain intensity. However, as expected, the magnitude of the WFT increases and the water ponding around the zero- slope section increases with increasing rain intensity. Effect of Longitudinal Slope Figure 11 shows the WFT distribution with different longitudinal slopes, Sx. The longitudinal slope changes the pavement reference height along the x-direction and this change is represented in the 3-D CFD geometric model. Figure 11(a) has zero longitudinal slope (Sx = 0%) and the water accumulates around the zero cross slope (ZCS) and the start of the inner side road (normal crown). With the longitudinal slope of 1%, shown as Figure 11(b), water accumulates further downstream on the inner side. Increasing longitudinal slope to 2% produces a smaller WFT or improved water drainage for both the ZCS region and downstream along the inner side. (a) Case 4 with Sx = 0% (b) Case 1 with Sx = 1% (c) Case 5 with Sx = 2% Figure 11. WFT distribution on pavement with different longitudinal slope. Effects of Macrotexture Road surface properties play a critical role in several key elements affecting hydroplaning. For example, the surface texture affects tire traction and the road evenness affects the vehicle and tire dynamics. The effect of surface macrotexture on WFT has received little attention in the published literature on hydroplaning. In this section, macrotexture is modeled and assessed relative to hydroplaning using WFT.

38 Based on Ks and Cs, the effects of macrotexture are investigated for 0 mm, 0.5 mm, and 1mm and a constant Cs=0.5. Figure 12 shows that with increasing macrotexture depth, the WFT increases near ZCS and along the inner travel edge downstream. The changes in WFT can be attributed to the increased roughness and therefore higher shear stress along the pavement. Therefore, the water film values increase slightly with the texture depth. However, these effects are small compared to the rainfall rate and longitudinal slope, road cross slope, and other road geometry changes. (a) Case 6 with texture depth of 0 mm (b) Base case with texture depth of 0.5 mm (c) Case 7 with texture depth of 1 mm Figure 12. WFT distribution on pavement with macrotextures. Type I and Type II Road Geometries (Two-lane Road; New Geometry) Finally, the two types of road segments were investigated to assess the transitions before and after superelevation. In the before superelevation transition (Type I road), the cross slope changes from normal crown to superelevated. In the after superelevation transition (Type II road), the cross slope changes from superelevated back to normal crown status. Figure 13 compares the WFT results of the two road transitions, and due to the symmetrically distributed slopes, the WFT distribution is relatively symmetric as well. The symmetrical distribution of the WFT is caused mainly by the symmetry of the slopes for Type I and II road according to the point C, as shown in Figure 8. Figure 13(c) and (d) also depict the velocity vectors of the moving water. In the normal crown section, where the cross slopes are −2% on both sides, the velocity vectors indicate that the water moves towards each edge of the road. However, at the ZCS section, the velocity vectors change direction toward the inner side of the road and show acceleration.

39 The vectors also show a change of flow direction near the ZCS section, which causes more water accumulation and explains the higher WFT value near ZCS. Meanwhile, the majority of the water flows downstream to the inner side road, which leads to relatively higher WFT values. (a) WFT distribution on Case 1 Type I road (b) WFT distribution on Case 8 Type II road (c) Velocity vector on Type I road (d) Velocity vector on Type II road Figure 13. WFT distribution and velocity vectors before and after superelevation. Summary A summary of the 3-D simulation results for the cases tested in this investigation is shown in Figure 14. The figure presents the data using different symbols to identify different locations along the road. Figure 14(a) depicts the relationship between WFT and rainfall rate. It shows that WFT is highest at the road–shoulder interface, and increasing rainfall rate significantly increases WFT at all locations. Figure 14(b) shows the effects of macrotexture depth. Although WFT does increase with increasing MTD, the effect of macrotexture depth on WFT clearly is not dramatic. Figure 14(c) shows that increasing longitudinal slope can both increase and decrease WFT near the shoulder, but will tend to decrease mean WFT for the entire road surface. Figure 14(d) shows the WFT results with 2, 4, and 8 lanes. The mean roadway WFT and the maximum WFT at 30 cm from the shoulders increase very slightly compared to the increase at the road–shoulder interface. The sensitivity analysis confirmed that one of the most important impact factors is the pavement geometry, which can be changed by cross slope, longitudinal slope, and the number of lanes. The sensitivity analysis of real road scenarios has shown several emerging patterns: the change in cross slope causes more water to accumulate near the ZCS position and the inward travel lane edge; and the increase in longitudinal slope leads to a higher WFT value near the end of the drainage path

40 around the end of the inner edge of the pavement but helps to improve the overall and mean water drainage in the normal crown and higher region of the superelevated section. (a) WFT changes with rainfall rates (b) WFT changes with MTD (c) WFT changes with longitudinal slopes (d) WFT changes with number of lanes Figure 14. Summary of the WFT Sensitivity Analysis. Further, as expected, the magnitude of the WFT increases with increasing rain intensity. Increasing texture depth increases the shear stress of the fluid near the pavement boundary, which increases the friction from fluid acting parallel to the surface and therefore slightly increases the water film values. However, such influence is relatively small compared to vertical height changes of the pavement due to road geometric characteristics.

41 Compared with the texture effects, the symmetrical water runoff patterns due to the symmetrical design of longitudinal and cross slopes have also proved the significant importance of geometric factors between the before and after superelevation transitions. The plots of fluid flow vectors have shown the reason for water ponding at ZCS and inward travel edge at downstream locations. The thicker WFT contributes to a higher HP. 3.2.4. Discussion The validation exercise proved that the CFD model produce reasonable results even in the complex geometries evaluated. In addition, the results of the sensitivity analysis were used to fine-tune the results of the simplified procedure presented in the next chapter. Although the model produced appropriate WFT predictions in the area of interest, the boundary condition on the sides of the modeled areas (outside the shoulders) created artificially high WFT on the edges of the simulation area, as illustrated in Figure 15. (a) Normal crown section (b) ZCS section Figure 15. WFT lateral profile at specified cross sections with different longitudinal slopes. Figure 15(a) shows the WFT profile at x = 5 m, which corresponds to the normal crown section and produces a symmetrical WFT distribution on either side of centerline for each case. Figure 15(b) shows the WFT profile at x =16.41 m, where the cross slope at the outside edge changes from −2% to 0%, which increases WFT in the ZCS area (0 to 3 m). Since the inner side maintains a slope of −2%, the inner edge becomes the lowest height in the cross section and higher water depths are produced due to large accumulation areas and gravity/momentum effects. 3.3. VEHICLE RESPONSE SUBMODEL The Integrated Hydroplaning Model models the vehicle response to water accumulation on the pavement by combining vehicle dynamics, tire, and hydraulic submodels (as presented in Figure 3). The tire-water-road interaction was modeled by combining detailed (Abaqus) tire and hydraulic submodels into an integrated FSI submodel that interacts with the vehicle response model, as follows:

42 1. The complex FSI simulation results are resolved at the spindle and include detailed tire- water-pavement interactions and estimates of the pressures and deformations along the surface of the tire (which is computationally intense). a. The FSI model predicts the movement that the tire induces on the water captured between the tire and the pavement and the induced tire deformation and longitudinal and lateral forces. The motion of the fluids (water and air) that interact with the tire is solved using the 3-D CFD solver Star-CCM+, adopting RANS equations, with conservation of mass, momentum, and energy. b. The motion of the tire, on the other hand, is solved by the Abaqus FE model. c. The coupling between the CFD model and FE model allows data transfer between the different solvers to represent the tire deformation. Star-CCM+ is built with commands to directly link to Abaqus, enabling communication and data transfer between the two solvers. Using direct co-simulation coupling will result in more efficient data transfer compared to the use of an external “middleware” software. 2. The complex underlying FSI simulation results are used to define simplified equations that depict the relationship between the road (including the geometry and WFT) and the spindle (including loads, positions, and velocities). These equations are used to replace the default tire model in the dynamics simulation (CarSim). 3. The integrated vehicle response model is then used to determine the HP for different driving maneuvers and WFT on different road geometries and surfaces. 3.3.1. Assessing Vehicle Handling Capabilities A key question for this research was estimating the effects of wet pavement on the handling capabilities of a vehicle. The limits of these capabilities must be defined in a clear way so that HP can be assessed for the various operating conditions. Specifically, we must assess the ability of the vehicle to brake and maneuver given the WFT and the other conditions (vehicle type, grade, superelevation, tire condition, etc.). To simplify the modeling of the vehicle response to the different road conditions and WFTs, we used two standard vehicle dynamics scenarios to determine the vehicle limit-handling capabilities: • Constant radius test – During the constant radius test, the vehicle operates at increasing speeds around a constant radius circle to determine the maximum lateral force that can be generated. This is (almost always) written in terms of maximum lateral acceleration relative to gravity (i.e., 9.81 m/s2 = 1 g). This is a convenient scaling factor when defining vehicle handling capabilities. Specific information from each of the four tire submodels (i.e., the vertical force, Fz, and the lateral force, Fy, for every operating condition) are needed for this vehicle simulation scenario. Again, these tire forces are derived from the complex underlying FSI simulation results. • Straight-line braking – The straight-line braking test involves a vehicle beginning at some prescribed initial speed, exerting the maximum braking effort, and then determining the maximum longitudinal deceleration. Again, this deceleration is described in units of g. Specific information from each of the four tire submodels (i.e., the vertical force, Fz, and the longitudinal force, Fx, for every operating condition) is needed for this vehicle simulation scenario.

43 The usefulness in defining the handling performance in terms of g-forces is illustrated by the estimation of friction. This concept simplifies the description of HP discussed in a later section of this report. The maximum lateral and longitudinal g-forces are often written as a vector combination in a “G-G” diagram (sometimes called the friction circle or friction ellipse). The G- G diagram can be used to describe the characteristics of the tire, as shown in Figure 16(a), or can be used to describe the characteristics of the vehicle’s limit-handling capabilities, as shown in Figure 16(b). (a) Tire Adapted from: Milliken and Milliken (1995) (b) Vehicle Figure 16. Examples of G-G diagrams. There are differences in the G-G diagram for a tire and the G-G diagram for a whole vehicle due to the vehicle dynamics. Specifically, the G-G diagram for the vehicle varies with vehicle parameters (e.g., due to brake proportioning, understeer gradient, etc.) and due to the vehicle dynamics (e.g., changes in vertical loads). For this reason, three different vehicle types were used in the simulations. It should be clear that the vehicle G-G diagram is the operating envelope for the vehicle in a given situation. The operating envelopes can be compared to determine the HP. This process is described in a subsequent section of this report. Briefly put, the HP is defined in terms of the degradation of the vehicle Performance Envelope due to the accumulation of water on the pavement, which is defined as the PM.

44 In order to simulate the constant radius test and the straight-line braking test, the tire model used in the vehicle simulation transforms geometric constraints to forces. Specifically, at each time step the vehicle prescribes the pose of the wheel, including spindle position (the spindle location with respect to the pavement) and orientation (slip angle, α). Given these geometric constraints on each of the four tires, the FSI tire submodel, which combines the tire model and tire-water-pavement interaction models, provides the forces (Fx, Fy, and Fz) to the vehicle model. Constant Radius Test To run the vehicle simulations for a constant radius test, the tire submodel must provide the relationship between the lateral force, Fy and the slip angle, α. These are typically represented as a family of curves, with each curve representing a constant vertical load, Fz. The FSI model was used to generate the relationships for the lateral force versus slip angle for different spindle positions. These simple relationships, derived from the complex FSIs at the tire-pavement interface, were then used as a simplified tire submodel for vehicle simulations. Straight-line Braking Test The straight-line braking test is harder to address, as we must relate the torque on the tire (braking) and the slip ratio, and this requires some approximations. For example, the braking effects can be considered by changing or increasing friction factors and shear stress in the FSI model or by directly adding extra forces in the Abaqus model. 3.3.2. FSI Model The interaction between the tire and the pavement when the surface is wet was modeled using Computational Domain We optimized the computational domain by balancing the accuracy of the representation of the hydrodynamics and tire dynamics with a reasonable computational demand. Figure 17 illustrates the domain adopted to capture the fluid dynamics, aerodynamics, and interactions with the moving tire. The software analyzes the fluid properties and movement, and the pressure and lift forces on the tire within a control volume. The spindle position of the tire is used as the reference point such that the tire rotates with angular velocity ωY, which can be calculated as the ratio between the vehicle speed, Vtire, and tire radius, Rtire. The air and water were defined with the relative velocity magnitude, which is the same as the vehicle speed but in the opposite direction of the moving tire. The pavement was also specified with the same relative velocity magnitude in the z-direction. The geometry of the tire was imported from the FE model of the tire (described in the following section) and then meshed in Star-CCM+. The mesh of the tire is shown in Figure 18 and the tire and the fluid domain are shown in Figure 19.

45 Figure 17. Schematic of the problem setup. Figure 18. Tire mesh profile. Figure 19. Mesh profile of the whole fluid model.

46 In Figure 18, the mesh is refined near the tire tread to capture the tread detail to predict the tire characteristics. To save computational time, the pavement mesh is fine near the tire and gets coarser away from the tire as shown in Figure 19. Two enlarged views of the tire mesh are shown in Figure 18. The yellow part shows the contact patch between the tire and pavement, which was dynamically coupled using Abaqus and Star-CCM+. Modeling Method The properties of the fluids are solved using the RANS equations for mass, momentum, and energy conservation in the generated control volume, similar to the one used for the water accumulation but using a much shorter time scale. The water surface is determined using the VOF model that solves the momentum equations and tracks the volume fraction of each of the fluids throughout the domain. The volume fraction of water, α, varies from 1, when the cell is full of water (the primary fluid), to 0, when it is full of air (the secondary fluid). A value between 0 and 1 indicates an interface between water and air. Figure 20 shows a schematic representing the VOF method to capture the free surface. Figure 20. Schematic of using the VOF method to track free surface. The k-ω turbulence model is employed to include two extra transport equations that determine the turbulent properties of the flow: the turbulent kinetic energy and the dissipation rate. The k-ω model is one of the most commonly used turbulence models since it can represent effects like convection and diffusion of turbulent energy. Boundary Conditions The boundary conditions have been specified to represent, as closely as possible, the actual scenario. The front surface numbered ① in Figure 17 is specified with uniform velocity corresponding to the vehicle speed (i.e., 70 km/h). The inflow of fluids from the front surface includes water and air. The pavement has been specified as the no-slip wall boundary condition with the same relative velocity. The WFT is added above the pavement as the initial condition with the same velocity of the pavement. The back surface numbered ③ is defined with atmospheric pressure for the outflow. The tire is modeled as a no-slip surface with no axial velocity but a constant angular velocity of 66.54 rad/s. To reduce the CFD computational domain and the effects of the no-slip boundary conditions on the flow field, the two side surfaces labeled ② and ③ and the top surface labeled ④ are specified as a free-slip boundary.

47 Preliminary Verification Figure 21 presents the velocity vectors for the pavement and the tire. The tire is rotating and shows a velocity of 70 km/h at the edge of the tire, but relatively smaller velocities with decreasing tire radius. An enlarged view of the velocity vector shows the road and tire moving in opposite directions. Figure 21. Velocity vector distribution on the pavement and the tire surfaces. Development of the Tire Model A Michelin premier A/S series tire (P205/60R15 91H) was modeled in Abaqus using a novel reverse engineering process developed at CenTiRe. The process first obtains the 2-D cross- sectional geometry data (Figure 22) and builds a 2-D geometric model in Abaqus. Figure 22. Reverse engineering process used to build a 2-D tire model in Abaqus.

48 The designed tread pattern of choice is added in a second step. Figure 23(a) shows the Michelin A/S Premier series tread pattern created using the commercial computer-aided design package Siemens NX. Figure 23(b) shows the completely meshed tread pattern in Abaqus. Figure 23. Tread design (a) created in NX and (b) meshed tread pattern created in Abaqus. The results of attaching the meshed tread pattern to the 2-D cross section are shown in Figure 24(a). The complete 3-D meshed tire model is obtained using the symmetric model generation feature in Abaqus (Figure 24[b]). Figure 24. 3-D FE model development in Abaqus: (a) section and (b) full 3-D model. Fluid-Structure Interaction Model This section presents the results of the FSI model using bold and treaded tires and variable slip angles. The integrated FSI tire model has been developed and coupled in two numerical models, the CFD model and the FEA model. In vehicle dynamics, slip angle is the angle between the rolling wheel and the pointing direction, which can significantly increase the hydroplaning risk with increased speed or braking, as shown in Figure 25.

49 Figure 25. Schematic of tire with slip angle. As previously mentioned, the fluid was modeled using 3-D RANS equations, including the conservation of mass, momentum, and energy. Initially, the 3-D geometry of the bald tire was imported from the Abaqus/CAE model where the vertical force of 3300 N was added to create the static tire deflection caused by the vehicle weight. Figure 26 shows the general overview of the mesh with details of tire, contact patch, and pavement, etc. The figure also shows the pressure distribution on the bald tire surface caused by water. Similar to the treaded tire, the region where water hits the tire’s surface has the largest pressure magnitude. Figure 26. Bald tire mesh profile and pressure distribution on bald tire surface.

50 The 3-D geometry of a treaded tire model was then imported from the Abaqus/CAE model, where the weight load of 4000 N was added to create the static tire deflection. Starting from the treaded tire model with zero slip angle, the operating conditions in the Abaqus/CAE model were changed to the 2-degree slip angle, which changed the profile of the contact patch. With the updated operating conditions in both the CFD and Abaqus models, the exchanged nodal information was updated for the mesh. Figure 27 shows a comparison of mesh profile and contact patch (highlighted) between the treaded tire with and without slip angle. As shown, although the slip angle is only 2 degrees, the contact patch changes significantly, and the profile is no longer symmetric. Figure 28 compares the lift forces for the two cases over time. The red line shows the lift force with 0 degrees of slip angle for the treaded tire simulation and the green line shows the lift force with 2 degrees of slip angle. (a) 0° slip angle (b) 2° slip angle Figure 27. Comparison of treaded tire mesh and contact patch Figure 28. Lift force profiles of 0° slip angle and 2° slip angle over time.

51 Figure 29 presents water distribution and splash for the bald tire, treaded tire with 0 degrees of slip angle, and treaded tire with 2 degrees of slip angle. The water film and the tire are presented from a 30-degree side view to depict the water splash near the front and the side of the tire. Comparing the cases of the bald and treaded tire, both of which do not have a slip angle, water splashes more because of the tread of the tire. Near the front region of the tire, water stays in the tread groove and moves rotationally due to the tire’s rotating motion. More water splash is also observed on the sides of the treaded tire. (a) Bald tire with 0° slip angle (b) Treaded tire with 0° slip angle (c) Treaded tire with 2° slip angle Figure 29. Water film distributions. A comparison of Figure 29(b) to (c) illustrates the effect of slip angle; the splash at the front of the tire is similar but there is significantly more side splash in Figure 29(c). The slip angle of the tire increases the water splash on the side, which increases the forces along the lateral direction. Finally, a model of a half-tread tire was developed to add an intermediate point between the bald and new tires. Experimental Verification of the Effect of Water on the Tire Contact Patch The tire contact patch modeled in Abaqus was validated using experimental data. Figure 30 shows the effect of WFT and slip angle on the tire contact patch. (a) 0° slip angle and dry pavement (b) 2° slip angle and 5-mm WFT (c) 0° slip angle and 5-mm WFT Figure 30. Effect of slip angle on the contact patch at 40 mph.

52 For the zero-degree slip angle, the difference in the tire contact patch between a dry pavement and a pavement with 5 mm WFT is illustrated by comparing parts (a) and (c) of the figure. The number of treads in contact with the pavement for dry conditions is shown in part (a). This number is higher than the number of treads in contact with the pavement for a 5-mm WFT shown in part (c) of the figure. A snapshot of the volume fraction of water on the pavement as the tire is moving is shown in Figure 31. The tire is rotating so that the vehicle moves from the right to left. The dark red color represents 100% volume fraction (water), and the dark blue color represents 0% volume fraction (air). We can see the effect of the tire in removing the water and channeling the water along the openings between the treads. Figure 31. Volume fraction of the water flowing in the tire pattern groove with 5-mm WFT at 40 mph. An instrumented VW Jetta (Figure 32) equipped with intelligent tires (sensors inside the tires) was used to collect data on the Virginia Smart Road under dry and wet conditions. The Smart Road’s weather-making towers were used to generate rain at a relatively high rate. Table 7 shows the recorded rainfall at different locations. Table 7. Approximate Rainfall Intensities Tower No. Total Rain (mm) Rain Duration (min) 36 79 59 41 70.8 59 46 129 59 Testing was performed closer to towers 36 and 41, where the rainfall was approximately 75 mm/h. Approximate WFT measurements of 2 to 3 mm were obtained with a ruler (accuracy in the order of 0.5 mm). A triaxial accelerometer placed on the tire inner liner was used to monitor the interaction between the tire and the road at a data collection rate of 5000 Hz. The triaxial accelerometer was used to determine the contact patch length. When the tire section where the accelerometer is placed is in contact with the pavement, the vertical acceleration is zero. This occurs when the accelerometer is in the contact patch area. The length of time per tire rotation that the vertical acceleration is zero can be combined with the speed of the vehicle to determine the contact patch length.

53 Figure 32. The instrumented vehicle used to collect the data. Figure 33 shows the contact patch length from the testing conducted on the Smart Road. In wet conditions, the contact patch is reduced due to the reduction in the normal load on the tire. The length of the contact patch on the dry surface at 50 mph is on average 128 mm and the length of the contact patch on the wet surface at the same speed is on average 118 mm. Figure 34 shows contact patch length calculated with the FE model. For wet conditions, the contact patch lengths from the measurement and the FE model are very similar. For the dry condition, the contact patch length from the FE model was about 9% larger. Figure 33. Measurements of contact patch length for dry and wet conditions performed on the Smart Road.

54 Figure 34. FE calculations of contact patch length for dry and wet conditions. Results of the FSI Simulations Once the model was completed, the research team ran a matrix of cases that included different speeds, tire pressures, slip angles, and WFTs. The results of these runs are summarized in Table 8 and were subsequently used to define a simplified tire model as discussed in the following section. These runs illustrate the effect of WFT on the lateral and vertical (lift) forces on the tire for all cases simulated. For example, the first four rows show how the lift force increases and the lateral force decreases with a thicker WFT, as depicted in Figure 35. Table 8. FSI Simulation Results for Lift and Lateral Forces Under Different Conditions Vertical Force (N) Speed (mph) Tire Pressure (psi) Slip Angle (degrees) Tire Tread(2) WFT (mm) Lift Force by water (N) Lateral Force (N) 4100 40 32 10 New tread 0 X 3566 4100 40 32 2 New tread 2 1050 1650 4100 (1) 40 32 2 New tread 5 1676 1490 4100 40 32 2 New tread 10 2205 1250 4100 60 32 2 New tread 5 2876 720 4100 40 32 5 New tread 0 X 2870 4100 40 32 5 New tread 5 1550 2210 4100 60 32 5 New tread 0 X 2980 5300 40 25 5 Bald 0 X 2376 2100 80 32 5 Bald 2 1250 1270 2100 60 32 10 Bald 10 815 2596 2100 40 25 5 Half Tread 0 X 1260 (1) Baseline state used for referencing the tire model parameters (represents a typical state of the system). (2) The tire tread description is converted to “wear” in the following section (to prepare the model summarized in Table 9): a tire with new tread has no wear (wear = 0) and a bald tire is completely worn (wear = 1). This way the baseline state has a wear value of zero.

55 Figure 35. Illustration of the effect of WFT on vertical and lateral forces on the tire. Simple Tire Model The values obtained with the FSI simulations were used to define a “simple” tire model, which follows the “Magic Formula” in Equation 1, which can be implemented in CarSim. The formula was developed by Hans Bastiaan Pacejka at Delft University of Technology and is extensively used in vehicle dynamic simulations (Bakker et al., 1987). 𝐹 𝛼, 𝐵, 𝐶, 𝐷,𝐹 = 𝐷 × 𝐹 × sin 𝐶 × atan (𝑎𝑡𝑎𝑛(𝐵𝛼 (17) where, 𝐵 ≈ 𝐵 + 𝐵 × ∆𝑠𝑝𝑒𝑒𝑑 + 𝐵 × ∆𝑝𝑟𝑒𝑠𝑠𝑢𝑟𝑒 + ⋯ 𝐶 ≈ 𝐶 + 𝐶 × ∆𝑠𝑝𝑒𝑒𝑑 + 𝐶 × ∆𝑝𝑟𝑒𝑠𝑠𝑢𝑟𝑒 + ⋯ 𝐷 ≈ 𝐷 + 𝐷 × ∆𝑠𝑝𝑒𝑒𝑑 + 𝐷 × ∆𝑝𝑟𝑒𝑠𝑠𝑢𝑟𝑒 + ⋯ 𝐵 = 𝑆 𝐵 = 𝑆 𝐶 = 𝑆 ... - 500 1,000 1,500 2,000 2,500 3,000 3,500 0 2 4 6 8 10 Fo rc e (N ) WFT (mm) Lift (vertical) Force Lateral Force

56 The parameters were optimized to predict the simulation results summarized in Table 9. The optimization used a cost function with two parts, the RMS error in the prediction of the lateral force and the variation in the parameters from the initial set. The RMS error in the lateral force estimate is straightforward. The quadratic cost function for the parameters was established using the following assumptions: • B should not vary more than about ± 1. • C should not vary more than about ± 0.2. • D should not vary more than about ± 0.1. • Speed varies by about ± 25 mph. • Tire pressure varies by about ± 7 psi. • Tire wear varies by about ± 0.5. • WFT varies by about ± 5 mm. • Vertical force varies by about 2000 N. Then, the normalized variations about the baseline state are the following: • varB = [1, 25, 7, 0.5, 5, 2000] ./ 1 • varC = [1, 25, 7, 0.5, 5, 2000] ./ 0.2 • varD = [1, 25, 7, 0.5, 5, 2000] ./ 0.1 The resulting weighting matrix is simply a diagonal matrix with the following entries: weight = [varB.2, varC. 2, varD. 2] The final optimized parameters of the “simple” tire model, which follows the “Magic Formula,” are summarized in Table 9. Each row represents a coefficient and each column a subscript (e.g., B0 =10.2950). Table 9. FSI Simulation Results for Lift and Lateral Forces under Different Conditions Variable Baseline Variation with Speed (mph) Variation with Tire Pressure (psi) Variation with Tire Tread (1) WFT (mm) Vertical Force (N) Coefficient / Subscript 0 1 2 3 4 5 B 10.2950 4.578E-04 -1.588E-04 -5.220E-05 -1.226E-03 -4.954E-04 C 1.9921 7.261E-04 9.116E-04 1.928E-03 1.217E-03 2.762E-05 D 0.5326 -3.527E-03 2.559E-03 -1.902E-01 -3.586E-02 4.981E-02 (1) (1 = new, 0 = bald) The fit of the selected model is illustrated in Figure 36

57 Figure 36. Fit of the selected simple tire model. 3.3.3. Performance Margin The PM is redefined as the additional performance capability that can be drawn upon beyond that which is demanded by the current operating condition. Clearly, this definition requires knowledge of both the current performance requirements and the limits of performance capability. The Performance Envelope describes the limit of performance capability that can be drawn upon during any maneuver (at any operating condition) and is developed in terms of the Effective Friction, which varies due to vehicle dynamics and road conditions (including roughness and WFT). A rigorous definition of the Performance Envelope is the first step in redefining the PM. The Performance Envelope is defined by the locus of points for which the vehicle traction equals the available vehicle traction (both acting in the road plane), or equivalently the locus of points for which the Required Acceleration equals the available acceleration (both acting in the ground plane). The components of the Required Acceleration for a given operating condition are written as (𝐴 ,𝐴 , indicated by a cross in Figure 37, and the components of the available acceleration are written as (𝐴∗ ,𝐴∗ ; shown as a solid curved line, where the asterisk notation denotes the boundary of the Performance Envelope. The intercepts with the axles are defined by the Effective Friction in the longitudinal direction, 𝜇 , effective coefficient of friction in the lateral direction, 𝜇 , longitudinal slope or grade angle, θs, and cross slope angle, θb, as explained later in this section. It is noted that the longitudinal and transversal slopes of the road are represented as percentages in most of the report but are analyzed as angles for the vehicle dynamic calculations.

58 Figure 37. Performance Envelope and PM (θb = cross slope angle; θs =grade angle) The PM for this operating condition is then the minimum difference between the Required Acceleration for a given operating condition and the locus of points that define the available acceleration. Some advantages of the redefined PM are that vehicles with different operating weights can be more easily compared. The units are the intuitive units of gravity, and the measure is convenient when including the dimensionless coefficient of friction. Specifically, the PM typically has a value between 0 and 1; the value of zero indicates that the limiting performance capabilities have been reached and there is no additional traction remaining. Values below zero indicate that the limit of traction has been exceeded. The upper bound of 1 may not be achievable except on dry pavement with specialized tires, but values of 0.8 are typical on clean dry pavement. However, for wet pavements the Pm will be much lower as the magnitude of the available Performance Envelope is reduced significantly. The PM is defined in terms of accelerations in the ground plane in physically meaningful and intuitive units of acceleration relative to the acceleration due to earth’s gravity. The choice of representing the measure in terms of acceleration allows comparison between vehicles of different sizes. The choice of defining the PM in the ground plane rather than the road plane enables a more concise representation of the concept and a more intuitive tool for geometric road design The Performance Envelope is a function of the road geometry and the Effective Friction. The Effective Friction is defined as the maximum fraction of the normal force that can be used to generate tractive force by the vehicle at the current operating condition. This definition of the Effective Friction is not the traditional relationship between a single tire and the road surface, but accounts for both the limitations imposed by the vehicle dynamics and road surface properties. For example, vehicles may be designed such that the front tires reach their performance limits before the rear tires in order to maintain directional stability (but not directional control). This implies that the front and rear tires will not reach their performance limits simultaneously and correspondingly the available vehicle traction is less than the theoretical tire force. The Effective Friction may also be influenced by road roughness and by water or contaminates on the surface. Consider two specific cases from which a generalization is made: the effective coefficient of friction in the longitudinal direction, 𝜇 , and effective coefficient of friction in the lateral direction, 𝜇 . Specifically, consider that the point on the Performance Envelope for pure braking is simply 𝐶𝑜𝑟𝑛𝑒𝑟𝑖𝑛𝑔 𝐵𝑟𝑎𝑘𝑖𝑛𝑔(𝐴 ,𝐴 (𝐴 ∗ ,𝐴∗ 𝜇 + tan𝜃 𝜇 + tan𝜃 𝑷𝒆𝒓𝒇𝒐𝒓𝒎𝒂𝒏𝒄𝒆𝑴𝒂𝒓𝒈𝒊𝒏 Performance Envelopes

59 the sum of the effective coefficient of friction, 𝜇 , and the grade; this point is indicated in Figure 37 by a small grey dot on the braking axis. Similarly, the point on the Performance Envelope for pure cornering is the sum of the effective coefficient of friction, 𝜇 , and the cross slope; this point is indicated in Figure 37 by a small grey dot on the cornering axis. The Effective Friction is developed for these limiting cases from which a generalization to the Performance Envelope is made. One issue to note is that the coefficient of friction acts in the road plane, which requires the Performance Envelope to be developed in both the road plane and ground plane. This road plane is identified with the addition of the subscript “𝑉” (for vehicle), which is consistent with SAE notation. Effective Friction Estimation To estimate the coefficient of friction, either through simulation or experimentation, two standard vehicle dynamics scenarios, performed on a flat road surface (i.e., the road plane is coincident with the ground plane), were used: a constant radius test and a straight-line braking test. During the constant radius test, the vehicle operates at increasing speeds around a constant radius circle to determine the maximum lateral force that can be generated. This condition is written in terms of maximum lateral acceleration relative to gravity. The straight-line braking test involves a vehicle beginning at some prescribed initial speed, then the maximum braking effort is exerted and the maximum longitudinal deceleration is determined. Again, this deceleration is described in acceleration units relative to gravity. Consider a particle model of a vehicle traveling in a circle prescribed in the ground plane. The effective coefficient of friction is estimated by increasing the vehicle speed until it slips free of the circle, at which point that maximum force in the horizontal plane is achieved 𝜇 =max(𝑚𝑎 )/𝑚𝑔 = max(𝑎 /𝑔). Accelerations defined in terms of units of gravity, or 𝑔’s, simplify this relationship, so that 𝜇 = max(𝐴 ) where 𝐴 is in units of 𝑔’s and the description of the friction using the vehicle’s acceleration information is straightforward. Similarly, 𝜇 =max(𝐴 ) for straight-line braking in the ground plane in a steady-state (trim) operating condition. For the remainder of this work, it is assumed that the vehicle is operating in steady-state (trim) operating condition, unless explicitly indicated otherwise. It should be clear that estimating the lateral and longitudinal coefficients of friction are critical in establishing the Performance Envelope and that this definition of Effective Friction accounts for losses due to factors including excess understeer, brake proportioning, load transfer, and powertrain limitations. Lateral Friction, 𝝁𝒚, Estimation Consider the simplest model of a vehicle traversing a properly banked curve (positive cross slope) without any grade (slope), as shown in Figure 38. The centripetal force, 𝐹 , which acts in the ground plane, is generated by forces acting on the vehicle in the road plane, 𝐹 and 𝐹 . Similarly, the vehicle operating weight, 𝑚𝑔, is balanced by forces acting on the vehicle in the road plane, 𝐹 and 𝐹 .

60 Figure 38. Lateral forces acting on a vehicle system. The centripetal force, vehicle lateral force, and vehicle normal force are given by Equations 18 to 20. 𝐹 = (18) 𝐹 = 𝑚𝑔 cos 𝜃 + 𝐹 sin𝜃 (19) 𝐹 = 𝐹 cos 𝜃 − 𝑚𝑔 sin𝜃 (20) The centripetal acceleration, 𝐴 , is defined in terms of longitudinal velocity and the radius of the turn in Equation 21. 𝐴 = and 𝐴∗ = ∗ (21) The required vehicle lateral acceleration is therefore given in Equation 22. 𝐴 = = 𝐴 cos 𝜃 − sin𝜃 (22) Consider the limiting condition of a banked, steady-state cornering maneuver, without any longitudinal force, in which the Performance Envelope is reached (as denoted with an asterisk). The available vehicle lateral force, 𝐹∗ is equal to the available traction and, using the definition of the coefficient of friction, is written as Equation 23. Note that the critical velocity at which this available vehicle lateral force is reached, 𝑣∗ , must be determined. 𝐹∗ = 𝜇 𝐹∗ = 𝜇 𝑚𝑔 cos 𝜃 + ∗ sin𝜃 (23) This available vehicle lateral force is normalized by the vehicle operating weight to form the available vehicle lateral acceleration, as shown in Equation 24. 𝐴∗ = ∗ = 𝜇 cos 𝜃 + ∗ sin𝜃 = 𝜇 cos 𝜃 + 𝐴∗ sin𝜃 (24)

61 Longitudinal Friction, 𝝁𝒙, Estimation Similarly, consider the particle model of a vehicle moving on an inclined road, as shown in Figure 39. The vehicle operating weight, 𝑚𝑔, is balanced by forces acting on the vehicle in the road plane, 𝐹 and 𝐹 . Figure 39. Longitudinal forces acting on a vehicle system. The tire braking force, and vehicle normal force are given by Equations 25 and 26. 𝐹 = −𝑚𝑎 −𝑚𝑔 sin𝜃 (25) 𝐹 = 𝑚𝑔 cos 𝜃 (26) The vehicle longitudinal acceleration is written as Equation 27. 𝐴 = − − sin𝜃 (27) Consider the limiting condition without any lateral force, in which the limit on the force generation capability of the vehicle is reached. For this case, the available vehicle longitudinal force, 𝐹∗ , is equal to the available traction. The force is written as Equation 28. 𝐹∗ = 𝜇 𝐹∗ = 𝜇 𝑚𝑔 cos 𝜃 (28) This available vehicle longitudinal force is normalized by the vehicle operating weight to form the available normalized vehicle longitudinal force, as shown in Equation 29. 𝐴∗ = − ∗ − sin𝜃 = −𝜇 cos 𝜃 − sin𝜃 (29) When the cross slope (banking) is zero and tan𝜃 ≪ 1 𝐴∗ = 𝐴∗ = −𝜇 − tan𝜃 (30) Since this study only deals with braking cases, not acceleration, the sign of the available vehicle longitudinal acceleration is changed from negative to positive for the convenience of analysis.

62 Summary of PM Results The maximum deceleration that can be accomplished in the absence of cornering forces is 𝜇 +tan𝜃 and the maximum lateral acceleration that can be achieved in the absence of longitudinal force is 𝜇 + tan𝜃 . The equation defining this elliptical locus of points is then: ∗( ) + ∗( ) = 1 (31) This is the equation that defines the Performance Envelope. That is, all points on the envelope (𝐴∗ ,𝐴∗ ) must satisfy this equation, which is a function of the Effective Friction (𝜇 and 𝜇 ) and the road geometry (slope, tan𝜃 , and cross slope, tan𝜃 ). The minimum difference between this envelope and the current operating condition is the PM, represented in the ground plane, in terms of units of gravity. Note that the distinction between the road plane and the ground plane must be accounted for when the grade or cross slope become large, which is outside the scope of the current work. Based on recommendation of the panel, the notation was changed for 𝐴∗ to 𝐴∗ to denote braking and the equation for calculating the PM was rewritten as follows: ∗ + ∗ = 1 (32) This revised equation represents an offset ellipse, which was used to compute the PM and incorporated in the Hydroplaning Potential Assessment Tool discussed in the next chapter. 𝑃𝑀 = min[ (𝐴∗ − 𝐴 ) + (𝐴∗ − 𝐴 ) ] ,∀𝐴∗ ,𝐴∗ (33) The PM is formally defined as the additional performance capability by the current operating condition. Specifically, the PM for a given operating condition is then the minimum difference between the Required Acceleration and the locus of points that define the available acceleration. The PM as defined in the equation is represented in terms of units of gravity. These concepts are exploited when defining the HP in the next section. With the new axis system defined, the ellipse is not center in the origin and the maximum deceleration is μx + tan θs, and the maximum lateral acceleration is μy + tan θb. Note that in the case of wet pavement, the coefficient of friction is also a function of the WFT, so that an effective coefficient of friction must be estimated. Specifically, μx and μy in the wet pavement condition refer to this Effective Friction. These two points on the Performance Envelope are each indicated by a bold “x” in Figure 37. 3.3.4. Vehicle Dynamics Simulation The CarSim vehicle simulation uses the relationships developed based on the highly detailed FSI tire simulations on wet and dry pavement, with different WFTs and tire conditions, as the tire submodel. This tire model is a functional relationship between the inputs, nominal tire deflection, slip angle, slip ratio, and vehicle speed (δ, α, κ, and v, respectively), and the force outputs (Fx, Fy, Fz). For each vehicle simulation, the WFT and tire condition are prescribed; the inputs for each of

63 the four wheels at each time step will be (δ, α, κ, v) and the outputs for each of the four tires are the resulting spindle loads (Fx, Fy, Fz). Hydroplaning Vehicle Simulator The vehicle response to different road conditions and WFTs is assessed using the two standard vehicle dynamics scenarios described in the previous sections for determining the vehicle limit- handling capabilities and estimating HP. In order to analyze a large number of conditions, we have developed a Hydroplaning Vehicle Simulator that can, by changing the CarSim simulation factors, automatically run a batch simulation written in MATLAB and CarSim’s own code language (.par file). The developed code changes CarSim’s simulation factors, such as vehicle type, maneuver, grade, cross slope, roughness, and tire model with respect to WFT. Synthetic road roughness data are also generated based on ISO Standard 8608. The user can execute the various batch simulations by selecting the simulation factors through the GUI shown in Figure 40. For each case, the simulation estimates the maximum acceleration for each case and produces a vector with all maximum g- forces. Figure 40. GUI of Hydroplaning Vehicle Simulator. The software, written in MATLAB, selects case numbers based on the simulation factors selected in the GUI, produces the CarSim files, and imports them to the CarSim automatically. Simulink then conducts all the simulations and generates the final outputs, a matrix of maximum g-forces, and a text file with the data for all runs. Figure 41 shows a schematic of the system configuration and Figure 42 shows the outputs of the simulator.

64 Figure 41. Vehicle simulation flow diagram. Example of Simulation Results The procedure discussed in the previous section was implemented in a simple tool for computing the PM. This is not the final product, but rather an intermediate program that is being used to model vehicle response and generate the PM in the background as an intermediate step to estimate the HP in the Hydroplaning Potential Assessment Tool. For illustration purposes only, Figure 43 show the effect of WFT and tire condition on the PM and HP. It must be noted that hydroplaning occurs when there is accumulation of water on the road; however, the vehicle can lose control even on a dry pavement when the PM approaches zero. Figure 42. Examples of G vector and output text file.

65 (a) Base Condition – New Tire (b) 5-mm WFT – New Tire (c) Base Condition – Bald Tire (d) 5-mm WFT – Bald Tire Figure 43. Performance Envelopes and PMs with four different conditions. Similarly, Figure 44 illustrate the effect of changing the vehicle type on the available friction and the speed at which the vehicle will reach the Performance Envelope (PM =0) for moderate braking on the same curve with a 5-mm WFT. In this case, the SUV will reach the critical condition at a lower speed (148 km/h) than the sedan (159 km/h) and the hatchback (182 km/h). Thus, the SUV is the most critical vehicle of the three considered. 3.4. HYDROPLANING POTENTIAL AND RISK The final step after estimating the PM in the process is to estimate the HP and hydroplaning risk on the roadway section under investigation. The approach proposed to implement these calculations is discussed in this section.

66 (a) Sedan (b) SUV (c) Hatchback Figure 44. Effect of vehicle type on the Performance Envelopes and hydroplaning speed. However, they are not implemented in the tool as, based on discussions with the project panel, it was decided that more validation data would be needed to support the calibration for the models proposed. This validation would require an extensive database of hydroplaning crashes linked with the roadway characteristics at the locations of the crashes. Since the research team only had access to data from Virginia, the panel felt that it was preferable not to implement the calculation of the hydroplaning risk. Therefore, the research team proposes a simplified approach in which the agency defines a critical PM and uses this to evaluate new or existing alignments. Once the agency gains experience, it can start fine-tuning the criteria and develop a more systematic approach based on crash analysis. This approach is consistent with the ones used for roughness and friction over the last couple of decades. 3.4.1. Hydroplaning Potential for a Specific Roadway Segment The PM can then be used to estimate the HP using, for example, a sigmoidal function, as follows:

67 𝐻 = 1 + (34) where: α and β are calibration parameters. The HP, thus, is developed as an estimated conditional probability of hydroplaning occurring for a specific vehicle under specific conditions. First consider several sets to describe the vehicle/tire combinations (V), pavement surfaces (S), and WFT (W) being analyzed. The elements of each set are formed by partitioning the set (the elements are mutually exclusive and collectively exhaustive). It should be clear that the partitioning of each set can be as course or as fine as is needed. For example, the set of all vehicle/tire combinations could be partitioned as simply as V = {car/bald, car/tread, noncar/bald, noncar/tread}, or the partition could be defined to specify every make, model, and year for each vehicle and make, model, and tread depth for each tire. The important concept is that whenever there is a valid partition of sets, the proposed approach is equally applicable. The results from applying these developments using different partitions will of course vary depending on the quality and refinement of the partition chosen. Similarly, the pavement surface partition could include many discrete ranges of grade, cross slope, roughness, texture, etc. Again, the important concept is that every combination of grade, cross slope, etc., is indexed separately so that the union of all the combinations accounts for all the pavement surfaces being considered without intersection. Next an event, H, is defined as hydroplaning occurring. The probability that hydroplaning occurs given a specific vehicle/tire combination, WFT, and pavement surface is: 𝑃(H | 𝑉 𝑆 𝑊) (35) Clearly, the true probability distribution is unknowable; therefore, it is proposed that the HP, be defined as an appropriate estimate of this conditional probability. Specifically, it is proposed that the HP be estimated from the PM, which is a function of the vehicle/tire combination, WFT, and pavement surface, as follows: 𝐻 = 𝑃(𝐻 | 𝑉 𝑆 𝑊) = 𝑓(𝑃𝑀(𝑉 𝑊 𝑆) (36) An admissible function can be developed to estimate this conditional probability. An admissible function of the PM should have the following required properties and desired characteristics: • Obey the axioms of probability. The function of the PM is an estimate of a conditional probability and, as such, must obey the probability axioms, including being nonnegative and less than or equal to one for all values of the PM. • Be a monotonically decreasing function of all positive finite values of the PM. The PM is a measure of any additional performance capability that is available. As the additional performance capability increases, the probability of hydroplaning must decrease. • Asymptote to one as the PM approaches zero. As the limit of vehicle performance is exceeded, the vehicle hydroplanes with a probability of 1. Furthermore, operating the vehicle very near the limit does not significantly change the HP.

68 • Asymptote to zero as the PM approaches infinity. There is always some very small probability of hydroplaning, even if the PM is large. The limiting case is that the probability of hydroplaning only approaches zero as the PM approaches infinity. The admissible function proposed as the HP is a logistic function of the PM presented in the previous section (which, in turn, is a function of the pavement surface, vehicle/tire combination, and WFT). 𝐻 = 𝑃(𝐻 | 𝑉 𝑆 𝑊) = 1 + (37) where: V = vehicle/tire combinations (e.g., sedan with new tires); S = pavement surfaces, including many discrete levels of grade, cross slope, roughness, and texture; W = WFTs; α = the value of the PM for which there is a 50% probability that hydroplaning will occur; β = the slope of the curve at PM = α. This parameter affects the sharpness of the transition. The parameter of the curve could be calibrated using a comparison of network-level hydroplaning crashes with data from fleets and roadways along the network. 3.4.2. Hydroplaning Risk It is further proposed that the hydroplaning risk, HR, is defined as an appropriate estimate of the conditional probability of hydroplaning occurring given the pavement surface being considered, specifically: HR = 𝑃(𝐻|𝑆) (38) Implicit in this statement is that hydroplaning risk is defined in the aggregate, for all vehicles at all times traveling over the road surface in question. This is accomplished by establishing the HP for sets of vehicle/tire combinations (V) and WFT (W) and then establishing the probabilities of these sets occurring for the road in question. First recall that V and W are partitioned so that: 𝑃(𝐻|𝑆) = ( )( ) = ∑ ∑ ( )( ) (39) 𝑃(𝐻𝑉𝑊𝑆) = 𝑃(𝐻|𝑉𝑊𝑆)𝑃(𝑉|𝑊𝑆)𝑃(𝑊|𝑆)𝑃(𝑆) (40) which is simplified by assuming that the vehicle driven is independent of the WFT and surface 𝑃(𝑉|𝑊𝑆) = 𝑃(V) (41) Combining

69 𝑃(𝐻|𝑆) = ∑ ∑ 𝑃(𝐻|VW𝑆)𝑃(V)𝑃(𝑊|𝑆) (42) So that 𝐻 = ∑ ∑ 𝐻 𝑃(V)𝑃(𝑊|𝑆) (43) The three terms in the double summation include: • An estimate for the first term is developed in this work as the HP. • An estimated probability distribution for each vehicle/tire combination on the road, which is outside the scope of this work. • The estimated probability of a particular WFT given the pavement surface (including the geographical location and rainfall), estimated by the WFT model.