A Computationally Efficient Approach to Resolving Vehicle-Induced Turbulence for Near-Road Air Quality

It is estimated that 30–45% of people in large North American cities live within 300 to 500 m from a highway or a major road, and exposure to traffic-related air pollution (TRAP) is associated with a range of health effects including reduced lung function, cardiovascular and respiratory problems, and exacerbated asthma [1]. Mitigating TRAP is essential to promote healthier and sustainable cities [2–4]. Computational fluid dynamics (CFD) tools have been employed in near-road air quality research, resolving physical and chemical processes [5], explaining experimental results [6], and assisting in developing mitigation strategies for traffic-related pollution through road and roadside designs [7–10]. Vehicle-induced turbulence (VIT) has a major impact on the initial dispersion of pollutants [11,12]. Earlier studies found that the contribution of VIT to the overall turbulence on-road varies from 45% to 95% depending on atmospheric, traffic, and meteorological conditions [13] and that neglecting VIT in dispersion modeling could lead to overestimating pollutant concentration by more than a factor of two in areas right next to the roadway [14]. A regional modeling study found that the lack of VIT contributes to over-predicting pollutant concentration in grid sections containing highways [15]. The incorporation of VIT in dispersion models has been shown to lead to a significant improvement in predicting pollution concentration [13,14]. Recent European studies show that with stricter regulations on traffic pollution there is a need for more precise modeling of VIT to accurately predict pollutant dispersion [16]. Similarly, exploring the effectiveness of road design options (such as road configurations, solid and vegetation barriers) to mitigate near-road air pollution also needs a better representation of VIT to assess the potential effectiveness of those designs on pollutant dispersion [9].

Several field measurements were carried out to quantify VIT, and the data collected were then used to parameterize VIT as a function of traffic volume and vehicle types [16–18]. However, since those parameterizations are developed based on specific roadways and are meant to be used in dispersion models, they are not generally applicable for CFD applications. Multiple CFD-based simulation techniques have been proposed to simulate VIT, notably, the wind tunnel and dynamic mesh methods. Both methods can resolve the realistic vehicle shapes. In the wind tunnel method, vehicles are held stationary in the computational domain and are subjected to flow coming from the inlet at a speed representing the vector sum of vehicle speed and ambient wind speed. This method produced reasonable results as demonstrated in earlier studies [19,20]. However, it is very challenging to model two-way traffic using the wind tunnel method. In addition, it is only applicable to the on-road domain because the inlet wind becomes unphysical for the near-road domain. To overcome this limitation, a multi-scale approach was introduced, where on-road and near-road simulations are tackled separately with the wind tunnel method being employed to estimate VIT [19]. While the wind tunnel method simulates the relative vehicle motion, the dynamic mesh method captures the absolute vehicle motion. For example, vehicles can be modeled as moving boundaries and a moving zone is created around each vehicle to add or remove layers of computational cells adjacent to a moving boundary [21]. In theory, the dynamic mesh method is the most realistic representation of VIT, but it is also the most computationally intensive. The computational costs are usually prohibitive to resolving complex on-road and near-road environments.

Therefore, there is a need to develop computationally efficient methods to simulate VIT with sufficient accuracy. The computational efficiency of modeling VIT can be improved if instead of explicitly modeling vehicles in the domain, a force or a turbulent source is injected into the domain to generate turbulence. By doing so, vehicles are not explicitly included in the domain and hence there is a significant reduction in mesh cell count and enhanced computational performance as the refined mesh required to resolve the fluid flow and turbulence around those vehicles is eliminated. This approach of injecting a force or turbulence source instead of having a physical geometry in the domain is implemented in other fields to reduce mesh size and computational cost. In the wind energy field, it is important to study the wake of turbine blades and its interaction with both turbines located downwind and the turbulent atmospheric boundary layer. To do so, instead of including the geometry of wind turbine blades in the computational domain, lines carrying body forces corresponding to the loading of the blades are used [22,23]. The reduced number of mesh cells needed to resolve the flow close to the blades allows for simulations of entire wind farms. Another example is modeling turbulent flows passing through vegetation [24,25]. Resolving the flow around individual leaves of vegetation in CFD would be computationally prohibitive. Therefore, equivalent drag force and turbulent kinetic energy (TKE) sinks are implemented in regions of the computational domain where vegetation is present. To implicitly model VIT using a similar technique, a few CFD studies inject a turbulent kinetic energy source in the traffic domain [26,27], referred to as the TKE method in this paper. However, the accuracy of the TKE method has never been evaluated.

This paper presents and assesses two novel methods, namely the force and moving force methods, introducing a force source to the momentum equations to account for VIT. We compare the accuracy of the two methods, together with the TKE method described above, against field measurements. The objective is to provide a much-needed assessment of robust implicit methods to model VIT using CFD and help future studies incorporate VIT to capture its impact on fluid flow and pollutant dispersion. The research products will assist researchers with better assessing TRAP concentrations in urban environments and developing mitigation strategies for TRAP to create sustainable and healthier communities.

This paper is organized as follows. Section 2 describes the three different methods to model VIT, the field measurement used for evaluation, the computational domain, setup, governing equations, and the modeling approach for each method. Section 3 discusses the results evaluating the three methods against the field measurements, and their impact on fluid flow and turbulence production. Section 4 concludes the paper and provides recommendations for future studies.

Figure 1 illustrates the previously adopted TKE method and two new methods introduced in this paper, i.e., the force method and the moving force method, in simulating VIT at a section of a highway with northbound (NB) and southbound (SB) traffic.

The TKE method entails creating a zone in the computational domain where a constant turbulent kinetic energy volumetric source is injected. The TKE injected in the zone emulates that produced by vehicle motion on the highway.

The force method injects a constant momentum force in the direction of motion of the vehicles. The constant force is the average drag force from the different vehicles in each traffic lane. This captures more of the physical phenomena that occur since the force from vehicles induces both fluid motion and generates turbulence. Additionally, VIT is anisotropic, with more fluctuations occurring in the direction of motion of the vehicles. By injecting a momentum source in the driving direction, the turbulence produced by the force method was also anisotropic. The force method also accounts for traffic flow direction. In our simulations, traffic flow in the northbound direction had a positive force while traffic in the southbound direction had a negative force.

The moving force method takes into account the convection of turbulence by fluid motion, as the surrounding fluid of a moving vehicle is being displaced around the vehicle and some of it moves in the vertical direction. Unlike the force method which applies a constant force throughout the traffic lane, the moving force method applies a constant force in a box that has similar dimensions to vehicles on road. That box force moves in the computational domain with the same velocity as the vehicles in each respective lane also accounting for traffic flow direction. The motion of the box in the domain produces a non-constant force in the driving direction that enhances vertical fluid flow.

To evaluate the performance of each of the three computational approaches to modeling VIT, the computational results were compared to field measurements that took place near I-15 in Las Vegas, NV. Those field measurements were part of a study to assess how complex urban roadway configurations affect local-scale air quality [28,29]. The section of the I-15 studied had ten lanes (five northbound, five southbound) and supported over 200,000 vehicles per day. Heavy-duty diesel trucks constituted approximately 10% of the total traffic volume, and the remaining were passenger cars and medium-duty vehicles. Figure 2(a) shows a top view of the layout of the field measurements. Four 3D Sonic anemometers were deployed at different locations to provide local meteorological and wind variability information. Two anemometers were located at the highway shoulder at heights 3 m and 6 m (sonics 1 and 2), one located 20 m downwind from the highway at a height of 10 m (sonic 3), and one located 100 m upwind from the highway at a height of 10 m (sonic 4). Figure 2(b) displays the locations of the anemometers. All anemometers collected 3D wind speed data at a frequency of 10 Hz.

Data from Sonic anemometers were collected from November 18, 2009 to February 1, 2010. However, lane-by-lane traffic flow information of vehicle speed, vehicle type, and traffic volume, necessary to evaluate the drag force, was only available from January 22, 2010 to February 1, 2010. In addition, we focused on the hours when the airport-based meteorological measurements were available and reported low wind speed (≈1 m/s) to ensure that the turbulence produced was mainly dominated by VIT and not due to other factors like structure-induced turbulence which can become significant at higher wind velocities. Furthermore, we selected hours when the wind was directed to the east at all sonics to ensure that sonics 1, 2, and 3 were located downwind of the highway and that VIT produced on the highway was convected to that sonics. As a result, we identified only 3 hours that met all the filtering criteria, i.e., availability of lane-by-lane traffic data, on-site anemometer measurements, and airport-based meteorological measurements, as well as low wind speed conditions. Two of those hours had similar traffic conditions resulting in similar TKE. Therefore, we presented two cases in this paper one with high traffic volume and the other with low traffic volume. Details of the two cases are highlighted in Sec. 2.7.

Figure 2(c) displays the computational domain used, where major geometric features were included: the bridge, the depression on both sides of the highway, and the trailer located 20 m downwind from the highway. The dimensions of the domain were 150 m (the driving direction) × 120 m (the span-wise direction) × 60 m (the vertical direction). Sonic 3 was located on top of the trailer at a height of 8.25 m. The height of the domain, 60 m, was about five times the height of the tallest object in the domain, to ensure there was no unphysical flow acceleration and prevent blocking effects [30]. To reduce the computational cost, the computational domain did not include 100 m upwind from the highway, rather only 2 m were included before the depression. As shown in Fig. 2(a), the upwind anemometer was located along an isolated trial way with no building and few barriers between it and the depression, and hence the TKE and wind velocity will not change substantially in that region.

We employed the Reynold’s stress model (RSM), which is a Reynold’s averaged Navier–Stokes (RANS) model that solves six transport equations for the individual Reynolds stress terms [32]. While the RSM model is one of the most computationally expensive RANS models, it can capture the anisotropy of turbulence as it does not use the isotropic eddy viscosity hypothesis used in other RANS models. VIT is anisotropic, and to evaluate how each of the three models captures it, the Reynolds stresses in different directions are needed and are accessible from the RSM model. Second-order schemes were used to solve for pressure, velocity, and Reynold’s stress terms. The relevant governing equations are provided in Sec. 2.5.

We conducted a mesh independence study to ensure that our results are insensitive to the mesh size. Figure 3 displays the TKE versus height at the highway shoulder for three different mesh sizes for both the force and moving force methods with a total of 13, 20, and 27 million cells, respectively. For both methods, the predicted TKE for both the 20 and 27 million cells mesh were similar to each other. Therefore, a mesh of 20 million cells was used for all the simulations in this study.

As mentioned in Sec. 2.3, two different traffic cases were chosen from field measurements to evaluate the computational methods proposed. For the high traffic case, a total of 6215 vehicles passed through both bounds, and the lane speed varied between 22.6 m/s and 30.3 m/s across the ten respective lanes. For the low traffic case, a total of 2335 vehicles passed through both bounds, and the lane speeds varied from 22.8 m/s to 28.5 m/s. Hourly average velocity was used to calculate the TKE at all sonic locations. Finally, to evaluate each of the three methods the percentage difference between the field measurements and the value obtained from the simulations was used as indicated in Eq. (10). A positive percentage indicates that the simulation overestimated the results while a negative percentage indicates that the simulation underestimated that value.

The height of the traffic zone in the computational domain affects the magnitude of either the TKE or force implemented in the zone. As demonstrated in both Eqs. (8) and (9), as the height of the zone L_z increases, the magnitude of the TKE or force implemented in the zone decreases. This in turn affects the TKE being produced in the domain. Three different zone heights, i.e., 1.8, 3, and 4 m, were evaluated. A height of 1.8 m was chosen because it corresponds to a passenger car height, a height of 4 m corresponds to a heavy-duty vehicle height, and 3 m is the average of both.

Table 2 displays the relative differences in TKE obtained at the three sonic locations with the three different zone heights for the high traffic case. From Table 2, a zone height of 1.8 m overestimated the TKE at sonics 1 and 2 at the highway shoulder while it came close to the field measurements value at sonic 3. A height of 4 m underestimated the TKE at all three sonic locations, which indicates that a zone height of 4 m resulted in a dilution of the force applied in the zone. A zone height of 3 m predicted the highway shoulder TKE well; however, it underestimated the TKE at sonic 3, i.e., 20 m away from the highway. Finally, for the high traffic and low traffic cases the heavy-duty vehicles represented 9.6% and 10.4% of the traffic flow but accounted for 48.9% and 50.4% of the total force, respectively. Therefore, it is reasonable that a force zone height of 3 m, which is the average height of a heavy-duty vehicle and a passenger car, produced optimal results. A zone height of 3 m was used for the rest of the discussion.

For the moving force simulation, it was important to explore the impact that the length of the box has on the turbulence produced in the domain. Similar to the zone height discussion for the force method, as the length of the moving box, L_box, increases, the value of the constant force implemented in the box decreases (Eq. (9)). Three different box lengths, i.e., 5, 11.5, and 18 m, were evaluated, corresponding to a passenger car length, a heavy-duty vehicle length, and the average of the two.

Table 2 displays the relative differences in TKE obtained at the three sonic locations with the three different box lengths for the high traffic case, which shows that the three different lengths overestimated the TKE at the highway shoulder at sonics 1 and 2 while they estimated the TKE at sonic 3 well within a reasonable range. In particular, a box of length 5 m produced the highest TKE at the highway shoulder location. This could be attributed to the fact that with a box length of 5 m the magnitude of the force in the box was too high generating more turbulence. By increasing the length to 11.5 m, the magnitude of the force applied in the box decreased and the TKE generated at the highway shoulder was closer to the values obtained from the field measurements. Further increasing the length to 18 m while it resulted in a decreased force being applied in the box, this reduced force was applied for a longer box length and resulted in an increase in the TKE at the highway shoulder. This demonstrates that the TKE produced was sensitive to both the overall force applied within the box and the length over which it is applied.

From the table, the TKE produced for moving boxes of length 18 m at sonic 3 was greater than the TKE produced by moving boxes of length 5 m, even though the TKE at the highway shoulder for a length of 18 m was smaller than that of the 5 m. This indicates that by increasing the length of the zone the overall vertical convection of turbulence was enhanced. Therefore, the TKE values for the three different lengths were within agreement with the field measurements at sonic 3. The moving boxes with a length of 11.5 m were chosen for the rest of the discussions since they performed better in comparison with field measurements at sonics 1 and 2.

VIT is anisotropic indicating that most of the turbulence fluctuations occur in the driving direction. U, V, and W refer to the crosswind velocity, velocity along the driving direction, and vertical velocity, respectively, while u′, v′, w′ are the fluctuations around the mean velocity in the crosswind, driving, and vertical directions respectively. The field measurements also showed that indeed the Reynolds stress term along the driving direction, v′v′, is two to four orders of magnitude higher than that of the crosswind, u′u′, and vertical direction, w′w′, at the highway shoulder. A computational method needs to account for that to properly assess VIT. Figure 5 displays the ratio of v′v′/ w′w′ obtained for the different methods at the three sonic locations for both the high and low traffic cases, which indicated that the turbulence is anisotropic for all three stations for both traffic cases since the ratio of the Reynolds stress terms is not one. Both the force and moving force methods were capable of capturing the anisotropy of VIT since the obtained magnitude of v′v′ is higher than that of u′u′ and w′w′ and the ratio of v′v′/w′w′ for both methods are close to that obtained from the measured values. Both methods captured the anisotropy of turbulence better at the highway shoulder at sonics 1 and 2 compared to sonic 3, which could be attributed to the fact that the decay of turbulence is modeled in addition to the rough estimate of the ABL used in the simulation that has a stronger influence over the TKE further away from the highway. Both the force and moving force methods inject a force in the driving direction, which resulted in most of the turbulence being produced in that direction. This is advantageous in comparison to the TKE method, which does not account for the anisotropy of VIT. The u′u′, w′w′, and v′v′ were nearly the same for the TKE method. As shown in Fig. 5, the ratio of v′v′/w′w′ derived from the TKE method is close to one for all three sonic locations, indicating that turbulence is almost isotropic. By directly injecting a TKE source, there is no directionality associated with it and therefore the TKE is nearly isotropic. Figure S2 available in the Supplemental Materials on the ASME Digital Collection displays the ratio of v′v′/u′u′ for the three methods at the sonic locations. Figure S2 available in the Supplemental Materials on the ASME Digital Collection also demonstrates that both the force and moving force methods captured the anisotropy of turbulence at the highway shoulder while the TKE method was nearly isotropic.

It is important to capture most of the physical fluid interaction that occurs between the vehicles and surrounding air to properly evaluate VIT. For example, by accounting for the flow around the vehicles, the moving force method captured some of the vertical convection of turbulence. Another important flow feature is the directionality of fluid flow due to traffic from the northbound to southbound. Both the force and moving force methods accounted for the bi-directionality of traffic. Traffic directionality affects both turbulence production and interaction. Figure 6 depicts the velocity V along the driving direction at a height of 2 m across the highway. The highway is centered around distance of zero and values greater than zero lie in the northbound region while values less than zero lie in the southbound region. Figure 6 shows that both the force and moving force methods have negative velocity in the southbound region and positive velocity in the northbound region indicating that both methods accounted for the directionality of traffic. Additionally, in each bound, there were fluctuations in the mean velocity accounting for the different velocities in each lane. The velocity was not zero at distance zero since the wind is directed to the east causing the velocity to shift slightly to the right. By contrast, the TKE method did not account for fluid flow. There are only minor fluctuations in the velocity due to the injected turbulence; however, the mean velocity was nearly zero. This demonstrates that for both the force and moving force methods the injected force resulted in a fluid motion that consumed some of the force whereas the TKE method converted all the drag force into turbulence which resulted in overestimating the TKE and minimal fluid flow.

Neglecting or improperly modeling VIT can result in significant under or overestimating of near-road pollutant concentrations. Proper incorporation of VIT in CFD modeling is critical for simulating transport and dispersion of traffic-related air pollution in on-road and near-road environments. We have presented two novel methods (force and moving force methods) that can model turbulence produced in complex traffic scenarios and compared then to the currently used TKE method. Using implicit modeling, i.e., instead of including vehicles in the computational domain, an equivalent force or turbulence source is injected into the domain. Those methods reduced the number of mesh cells that otherwise would have been needed to resolve the flow around vehicles. This implies that the methods proposed are computationally robust and can model intricate traffic scenarios that would have been otherwise computationally prohibitive using explicit methods.

The TKE method, which is currently used in CFD studies to account for VIT, resulted in an overestimation of VIT produced on the highway due to converting all the equivalent drag force into turbulence. Additionally, injecting a TKE source does not have a direction associated with it and hence it did not account for turbulence anisotropy and fluid convection due to vehicle motion. Both the force and moving force methods involved injecting a force into the domain and both methods captured turbulence anisotropy, traffic directionality, and fluid convection due to vehicle motion.

The force method matched the field measurements at the highway shoulder for both low and high traffic scenarios within a reasonable range. The force method is useful when steady and averaged traffic scenarios are to be modeled. The moving force overestimated VIT at the highway shoulder; however, it captured more of the flow phenomena around vehicles like vertical fluid convection. A simple model was used for the moving force method where the force was applied in a moving box. Future improvements for the moving force method include modifying the shape over which the force is applied and potentially having an uneven distribution of the applied force. With further improvements, the moving force method could better capture flow structures from individual vehicles in addition to varying structures from different vehicle types.

The Cornell team would like to acknowledge the funding support from the National Science Foundation (NSF) under Grant No. 1605407. We especially want to thank Dr. Kirk L Clawson from the NOAA ARL Field Research Division for providing meteorological data from Sonic anemometers as a part of the Las Vegas Roadway Vehicle Turbulence Study. This paper has been subjected to the U.S. Environmental Protection Agency review and approved for publication. Approval does not signify that the contents reflect the views of the Agency nor does mention of trade names or commercial products constitute endorsement or recommendation for use. The views expressed in this article are those of the author(s) and do not necessarily represent the views or policies of the U.S. Environmental Protection Agency.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.