Development of a Numerical Investigation Framework for Ground Vehicle Platooning

This paper presents a study on the flow dynamics involving vehicle interactions. In order to do so, this study first explores aerodynamic prediction capabilities of popular turbulence models used in computational fluid dynamics simulations involving tandem objects and thus, ultimately presents a framework for CFD simulations of ground vehicle platooning using a realistic vehicle model, DrivAer. Considering the availability of experimental data, the simulation methodology is first developed using a tandem arrangement of surface-mounted cubes which requires an understanding on the role of turbulence models and the impacts of the associated turbulence model closure coefficients on the prediction veracity. It was observed that the prediction accuracy of the SST k−ω turbulence model can be significantly improved through the use of a combination of modified values for the closure coefficients. Additionally, the initial validation studies reveal the inability of the Unsteady Reynolds-Averaged Navier-Stokes (URANS) approach to resolve the far wake, and its frailty in simulating tandem body interactions. The Improved Delayed Detached Eddy Simulations (IDDES) approach can resolve the wakes with a reasonable accuracy. The validated simulation methodology is then applied to the fastback DrivAer model at different longitudinal spacing. The results show that, as the longitudinal spacing is reduced, the trailing car’s drag is increased while the leading car’s drag is decreased which supports prior explanations of vortex impingement as the reason for drag changes. Additionally, unlike the case of platooning involving Ahmed bodies, the trailing model drag does not return to an isolated state value at a two car-length separation. However, the impact of the resolution of the far wake of a detailed DrivAer model, and its implication on the CFD characterization of vehicle interaction aerodynamics need further investigations.


Introduction
In recent years, the influence of climate change and a push toward sustainable environmental practices have begun to change the face of the automotive industry. Governing bodies, such as the EU, have initiated actions and policies to reduce emissions by 60% in the coming years [1] which has led to a push in research for more sustainable methods of vehicle traffic, such as platooning. Platooning, also known as convoying, is the act of two or more cars traveling in the aerodynamic wake of a lead car in order to improve gas mileage or gain a competitive edge in motorsports. This practice, also known as drafting or slipstreaming, is seen quite heavily in automotive racing where racers travel in extremely close proximity, this is done in order to conserve fuel and increase top speed through a reduction in drag. While this practice is less seen in consumer vehicles, among heavy transport vehicles, such as freight liners, the practice is quite popular and can produce great returns in fuel savings [2][3][4], to cite a few. With the rapid advances in self-driving technology it is believed by many that soon even consumer vehicles will begin platooning in a connected transportation system due to the decreasing follow distances that self-driving cars can achieve with their low reaction times and vehicle to vehicle communication [5,6].
The increased interest in platooning has produced much research in the optimization of platoon formation with Sun et al. [7] developing a method of optimal platoon formation based on the relative benefit to the platoon member's gas consumption and speed. Deng [8] developed a framework which could be used to simulate mixed traffic flows that included heavy vehicle platoons (HVP) and found ways to optimize the formation times of platoons. While CFD research has been done on platooning bodies, the vast majority has been limited to simplified bodies such as the Ahmed body or non-detailed car geometries like in the work of Gan et al. [9]. McNabola et al. [10] performed experiments of the intake of pollutants info a car's ventilation system at low speeds and found that at low speeds increasing follow distance drastically decreased the intake of pollutants into the system. They then continued the work into CFD simulation of heavily simplified automotive body, finding that the k − model was able accurately predict the intake of pollutants. Ebrahim et al. [11] used the Realizable k − model to compare CFD results of electric passenger vehicle platoons to data gathered through track testing and found that the Realizable k − was not able to predict the drag on the vehicles with great accuracy, leading to a large discrepancy between the calculated power usage of the vehicles and the actual. A notable CFD work in the study of drafting is Pagliarella [12] who performed wind tunnel testing on a variety of slant angles, determining that the vortices created from the rear of the lead Ahmed body were impinged upon by the trailing body. Rajamnai [13] also performed CFD on up to six Ahmed bodies and saw that at close spacings the lead body decreased in drag while the trailing body increased in drag. He then noted that at large spacings between x/L = 1 and 2 that the drag on the bodies become constant as they resumed a flow like that of the isolated body. Papers from Kim and Le Good [14,15] both used simplified bodies to investigate the aerodynamic impact on the geometry of the platoon members and its effects on the drag of the platoon. Also notable is the work that has been done on HVPs by the likes of Watts et al. [2] who used CFD studies to analyze the interactions of tractor trailer platooning as well as the flow around surface mounted cubes in a tandem orientation. This work was followed by works such as Mokhtar [3] who investigated ways of lowering drag through transverse offset between vehicle and Tornell [4] who found that even at a follow distances of 20L the trailing vehicle still saw a large reduction in cooling flow.
In recent years a fully detailed open source car model has been introduced called the DrivAer as a joint project from Audi AG, BMW, and Institute of Aerodynamics and Fluid Mechanics of the Technische Universitat Munchen (TUM) [16][17][18]. This geometry presents three forms a fastback, notchback and square-back. Along with the publishing of this geometry, the same group released experimental data containing load data, surface pressure, and the effect of ground simulation. CFD was also released using the k − ω SST which can be used for simulation validation. Soares et al. [19] added to this data by performing RANS simulations on the far wake, determining that for the fastback geometry side mirror vortices persisted up to 1L behind the car and that other vortex cores were able to persist up to 3L behind the car. Ashton et al. [20] compared the DES methods to the RANS simulation methods and found that the IDDES SST k − ω improved the drag prediction but still overestimated the recirculation length due to an over-prediction of TKE in the initial separation. With all this considered, the DrivAer geometry provides a great platform for a fully detailed platooning simulation which will be performed in this paper. However, the geometry has one drawback, because of its relatively recent release, large amounts of data are not available as they are for some older testing geometries like a surface mounted cube or Ahmed body.
To address this problem, this paper seeks to create a high fidelity CFD simulation framework which can be used for the testing of detailed car geometries in a platooning scenario. Because of the abundance of flow feature data, this methodology starts by using wall mounted cubes to tune and test the accuracy of turbulence models using URANS and IDDES. Once the turbulence models have been tuned and proven, they will then be applied to a tandem wall mounted cube case of various separations in order to test the validity of the tuned turbulence models in a tandem case. Note that works presented in this paper are extensions of the master's thesis work of Sudhan Rajasekar [21]. While a wall mounted cube is not an exact analogy to a detailed car model, it will allow the tuning and validation of turbulence models and meshing schemes to be carried out in a less computationally demanding simulation which also has large amount of literature investigations. Once the methodology has been established, the same procedure will then be applied to the DrivAer model simulations.

Governing Equations
The Navier-Stokes (N-S) equations are the basic governing equations for viscous fluid flow. These equations express conservation of mass (also known as the continuity equation) and conservation of momentum (derived by applying Newton's second law of motion to a fluid element) for a Newtonian fluid. These are shown in Equations (1) and (2), respectively, for an incompressible fluid flow using Einstein summation notation where repeating index variables (i) imply summation over all possible values (i = 1, 2, 3).
In Equations (1) and (2), t represents time and the variables u i , p, ρ and T ij represent the instantaneous values (a function of time) of the velocity in x i direction, pressure, fluid density, and fluid viscous stress tensor, respectively. The viscous stress tensor, T ij , is defined as: where µ is the fluid kinematic viscosity and s ij represents instantaneous rate of strain tensor defined as: The Navier-Stokes equations describe the complex fluid flow accurately and in complete detail. The computational cost of solving for the flow-field around a full-scale DrivAer model (the object of our interest) at a speed of 80 mph using a Direct Numerical simulation (DNS) approach would require a near body mesh size of about 264 trillion cells and would require 30 thousand time-steps to simulate 1 s of physical time [22]. As such, the DNS is normally used as a powerful research tool to understand turbulent flows at moderate Reynolds numbers. Furthermore, as the N-S equations describe every detail of the turbulent flow field from the largest to the smallest scales, the volume of information is overwhelming and requires statistical methods to get a useful solution for real flows of engineering interest.

Reynolds Averaged Navier-Stokes (RANS) Approach
Reynolds Averaged Navier-Stokes equation are achieved by first decomposing the time-dependent velocity field u i into mean (U i ) and fluctuating (u i ) components using the Reynolds decomposition, mathematically expressed as u i = U i + u i , and then ensembleaveraging the original N-S equations. The RANS equations in the conservative form are expressed in Equations (5) and (6). Note that, we described here the turbulent flow statistically in terms of the mean velocity field U i (x, t) and mean rate of strain S ij (x, t) instead of the instantaneous velocity field u i (x, t) and instantaneously rate of strain field s ij (x, t), respectively, and these equations are commonly referred to as the Unsteady Reynolds Averaged Navier-Stokes (URANS) equations.
By comparing Equations (2) and (6), it is clear that the Reynolds averaging process introduced six additional terms, −ρu i u j , known as the Reynolds stresses (τ ij ). The Reynolds stress tensor, τ ij , is symmetric tensor, and has six independent terms which gives rise to the issue known as the turbulence closure problem, as we now have 10 variables to determine using 4 equations. One way of solving this closure issue is using the turbulent-viscosity hypothesis introduced by Boussinesq in 1877 (see Equation (7)) which related the turbulent stresses with the mean rate of strain via a flow variable turbulent eddy viscosity, ν t , the determination of which is the central element of turbulence modeling approach.
where k is turbulence kinetic energy per unit mass, k ≡ (1/2) u i u i , and δ ij is Kronecker delta. Subsequently, eddy viscosity based turbulence models differ (1) in the way they calculate eddy-viscosity, (2) how the eddy viscosity coefficient, ν t , is used to relate the mean turbulence field to the mean velocity field, and (3) what additional transport equations are solved in order to determine ν t . The current study restricts itself to three different two equation turbulence modeling approaches: Realizable k − (RKE) [23], SST Menter k − ω (SST) [24,25], and v 2 − f (V2F) [26,27] models. For brevity, we are omitting the transport equations and constitutive relationships for the V2F and RKE models, however, an interested reader is directed to Zhang et al. [28] for the details. As will be seen later that, although the SST model appears to be superior to the other two models, for improved flow predictions, some of the turbulence model closure coefficients of the SST required fine tuning. As such, for the benefit of the reader, we are including below the details of SST model which uses two additional transports equations: one for the turbulence kinetic energy, k, and the other for the specific turbulent dissipation rate, ω.
The well-known sensitivity to free-stream/inlet conditions flaw of the standard k − ω model of Wilcox [29,30] is addressed in the SST model of [24] through blending the k − ε model of Launder and his coworkers [31,32] in the far-field with a k − ω model in the near wall region. The relevant transport equations for the steady incompressable formulation of the SST model are given below.
where α, β, β * , σ k , σ ω and σ ω2 are turbulence model closure coefficients. F 1 and F 2 are blending functions; F 1 blends the near-wall value or the k − ω model value of a model closure coefficient with its k − ε or far away from the wall value. That means, a closure coefficient φ is determined from its k − ω value φ 1 and k − ε value φ 2 , as The default values for the coefficients corresponding to the k − ω region are: β 1 = 0.0750, σ k1 = 0.85, and σ ω1 = 0.5. The default values for the coefficients corresponding to the k − ε region are: β 2 = 0.0828, σ k2 = 1.0, and σ ω2 = 0.856. Default values for the other model closure constants are: von Kármán constant κ = 0.41, β * = 0.09, and α = 1.

Detached Eddy Simulation (DES)
One of the major flaws of the RANS approach is that it models all scales, which often causes a poor prediction in separated flows dominated by multiple-scales of motions. Thus, scale-resolving simulation (SRS) approaches, such as, the Large Eddy Simulation (LES), which resolves the large scales of turbulence and models only the small-scale motions, seem to be more attractive alternatives. Large Eddy Simulation (LES) was proposed based on the argument that small-scale turbulence contribution to Reynolds stresses is much weaker when compared to largest eddies, and can be modeled without much loss of accuracy. Thus, LES segregates the wide turbulent scales into the resolved larger scales and the modeled sub-grid scales (SGS). This means that the grids need not be as fine as required for a DNS study, and it is, thus, possible to simulate higher Reynolds number flows with the LES than at an order of magnitude less computational cost less than the DNS. However, for most of the automotive external aerodynamics flows, LES is still too computationally costly as it requires billions of cells and millions of CPU hours.
This issue is somewhat mitigated in the Detached Eddy Simulation (DES) approach proposed by Spalart [33,34]. DES is a hybrid modeling approach which combines the best of both RANS and LES: a RANS simulation in the boundary layers (which would require very fine mesh to be resolved using an LES approach) and a LES in the unsteady separated flow regions. This greatly reduced computational requirement as RANS, unlike LES, does not require very fine isotropic cells in the boundary layer. In each cell, the effective turbulent length scale, T , given in Equation (18) and the local finite-difference cell size, LES , are computed.
where C DES is a constant and ∆ DES is the filter width based on local mesh resolution, and wherever in the domain T < LES , RANS is implemented. A number previous studies showed that the DES methodology is prone to producing an artificial separation, known as the Grid Induced Separation (GIS), when refining the maximum cell edge length inside the boundary layer below a critical value (see Menter et al. [35]). This fictitious separation is an artifact of extra fine-mesh used in resolving a wallbounded flow and can be avoided by shielding the RANS model from the DES formulation for wall boundary layers. Spalart and his coworkers had an extensive investigations of this phenomena and suggested a shielding function in terms of the eddy-viscosity and the wall distance [36]. This DES scheme, which incorporates a delay factor that enhances model's ability to differentiate between the LES and RANS regions on meshes where spatial refinement could give rise to ambiguous behavior, is popularly known a the Delayed Detached Eddy Simulation (DDES). The Improved Delayed Detached Eddy Simulation (IDDES) model is an improvement over the DES and DDES, and consists of a combination of various new and existing techniques aimed at making the DDES approach versatile for high Reynolds number flows, see [37,38]. In addition to providing a shielding against GIS, and, thus, preventing the overestimation of dissipation due to premature switch from RANS to LES, the IDDES approach also allows the model to run in Wall-Modeled LES (WMLES) mode, for a more efficient resolution of wall bounded flows. However, the model formulations are relatively complex for both DDES and IDDES. Again, for brevity, we are including the equations relevant to the calculation of the specific dissipation rate and C DES only; an interested reader is referred to the original publications of Shur et al. [39] and and Gritskevich et al. [38] for additional details.
The specific dissipation rate and the blending function in the RANS-SST is modified for IDDES-SST asω where, f β * is the free-shear modification factor and Hybrid is the hybrid length scale determined based on the effective length scale and the local grid size using a blending function and an elevating function.
where C DES,k−ε = 0.61 and C DES,k−ω = 0.78, are the values of C DES in the k − ε and k − ω regions of the SST-RANS model, respectively.

Computational Setup and Boundary Conditions
All simulations performed in this study were done using a commercial finite volume solver CFD code, STAR-CCM+. Simulations were performed with a segregated incompressible unstructured solver using the SIMPLE algorithm. In order to solve the momentum equations and turbulent quantities, an Algebraic Multi-grid (AMG) solver was used with V cycle and flex cycle. A two-layer y + wall-function was used in order to gain accurate results even in areas where the mesh was unable to remain fine because of sharp or complex model geometry. A CFL number independence study was performed and a CFL number of approximately 4 was chosen for all URANS simulations and a CFL number of approximately 1 for all IDDES simulations. Table 1 shows the time-step used for each of the simulation cases ran in this study. In order to create a free flowing virtual wind tunnel with minimal blockage, a blockage ratio, based on the frontal area, less than 0.3% was maintained for all cases. For the surface mounted cube, the domain was separated into two distinct sections, a slip boundary and no-slip boundary. The slip boundary was located from the front of the channel to 2L before the front face of the cube, where the no-slip boundary began and continued to the end of the tunnel. This slip boundary was used in order to ensure that the flow was fully developed and that the boundary layer thickness was comparable to that of Martinuzzi and Havel [40], whose work will be the main comparison for the current study. The inlet condition velocity was set for 8.8 m/s which when combined with a cube height (h) of 40 mm was used, producing a Reynolds number of 22,500. For the DrivAer cases, a 40% scale model was placed in a domain with dimensions of 46L × 26W × 12H in the stream-wise, cross stream, and normal directions, respectively, where L is the wheelbase of the vehicle, W is the trackwidth of the vehicle, and H is the height of the vehicle. The DrivAer model was first simplified and surface meshed using a commercial pre-processing software, ANSA. The wheel rotation was simulated using a moving reference frame method to better simulate the pumping effect of the wheels. A no-slip patch was placed under the car and extended 0.5L in front of the car, 10.5L behind the car, and 2.5W on either side of the car. To simulate the ground effect, the patch was given a tangential velocity of 40 m/s, same as the inlet velocity. This inlet velocity of 40 m/s corresponds to a Reynolds number of 4.87 × 10 6 .
An overset mesh technique was chosen for its cost effectiveness between cases, ease of use in transient simulations, and ability to be modified for future works. Figures 1 and 2 show the generated mesh for the cube in its respective virtual wind tunnel. The base size of 2 mm was selected by performing a mesh independence study varying the base size from 4 mm down to 1 mm. The mesh size around the cube was gradually increased using three volume sources which surround the cube. The difference in drag prediction when moving from 2 mm to 1 mm was less than 0.1%, so a 2 mm (0.05 h) base size was selected for computational efficiency. The final meshed cube domain contained 5.7 million unstructured trimmer cells for the isolated case and 7 million for the tandem case.  The DrivAer geometry was then meshed in a similar fashion to the cube. The water tight geometry was imported from ANSA where volume sources were then added. The volume sources were made in the shape of the car. Four main volume sources were used with the first composed of 20 layers of 6 mm base size cells which extended 0.5 length behind the car. The second volume source was made of 15 layers before transitioning to the final volume source that extended 3L behind the car. A final volume source extended 8L behind the car in order to act as overlap between the overset and background regions. Using this meshing strategy yielded a final isolated model mesh of 27.7 million cells, shown in Figure 3 and 54.2 million for the tandem case. In both meshing strategies prism layer cells were used over all external surfaces in order to capture the high gradient small-scale motions in viscous regions the flow. Table 2 shows the settings used to generate the prism layers for each geometry, where TPLT represents the total prism layer thickness. While wall functions were used in this study, all y + values were kept below 1 for the wall mounted cube and below 5 for the drivAer geometry. Figure 4 shows scalar scenes of the y + values achieved for this study both for the wall mounted cube and the DrivAer model. These y + values were achieved using anisotropic prism layer cells near the no-slip boundaries of the domain. The slight increases in the DrivAer y + values are due to the complexity of the model but the authors believe this to still be in an acceptable range. Table 2. Prism layer meshing details for both the surface mounted cube and DrivAer models.

Results and Discussion
The first step in testing methodology was to analyze three different common RANS turbulence models to see which model gave the most accurate predictions and showed the most room for prediction improvement. These testing simulations were performed using the isolated cube and three different RANS models. The three turbulence models chosen for this investigation included Realizable k − (RKE) [23], SST Menter k − ω (SST) [25], and v 2 − f (V2F) [27] models. It was seen that over the cube the SST model showed the best correlated prediction of the overall flow field. Furthermore, this model was also found to have the most literature available on the improvement of model predictions. The SST's closure coefficients were then tuned through a testing regime which involved varying the closure coefficients β * , σ ω1 , and σ ω1 which have been shown by Bounds et al. [41] to be the most sensitive coefficients. The results of each variation are shown in Table 3, where the length of the recirculation bubble is the main flow feature monitored, and was used as the criterion to determine the veracity of the chosen coefficient values. The final closure coefficients were β * = 0.07, σ ω1 = 1.1, and σ ω2 = 1.7, for details see [21]. These coefficients showed the best flow feature agreement to the experimental recirculation of X R = 1.7h with a predicted length of X R = 2.8h. This improved set of coefficients created a vast improvement over the prediction of X R = 5.3h by the default set of coefficients. The SST model with the new coefficient values was then tested in an isolated cube setup in order to validate the turbulence model's predictions in both URANS and IDDES schemes. To investigate the improvements made to the IDDES and URANS models, Figure 5 shows the streamlines on the symmetry plane of the current CFD compared to experimental works of Martinuzzi (see Yakhot et al. [42] and Rodi et al. [43]). It can be seen in the figure that the URANS case slightly out performs the IDDES in the prediction of the separation bubble on the top of the cube and the vortex core placement of the rear recirculation bubble. The IDDES case predicted the vortex core slightly higher and further downstream of the experimental case thus generating a wake which is slightly larger than the experimental case and modified URANS case. The URANS case on the other hand, while it predicts the location of the vortex core in good agreement with the experiment, under predicts the size of the vortex core. When viewing the top streamlines, the IDDES model predicts the flow field with better accuracy than URANS. Both models failed to achieve agreement with experimental recirculation length, with both models predicting a length around 2.5h while the experimental value was 1.7h. This is believed to mainly come from the difference in blockage ratio used in the experimental work versus the CFD work. This can be further supported by works, such as Castro and Robins [20] where they showed that the blockage ratio in experimental tests can constrain the wake, an effect which increases with blockage ratio. While the URANS slightly outperforms the IDDES model in the near flow field, the IDDES model is able to predict with great accuracy the far field wake of the cube, an area where the URANS model completely fails. Figures 6-9 show the spanwise velocity profile at four different locations compared to the experimental work of Castro [20]. Both models show good agreement for the predictions at x/h = 2, however, at x/h = 6.5 only the IDDES model is able to predict the peaks which occur along the block's axis. This also agrees with the past work of Watts [2] in which he explored the effect of two platooning trucks at a separation distance of 1000 ft. Watts reported that the URANS model failed to reproduce the correct wake and even caused residual effects which influenced the trailing truck. While the improvement of the URANS model should not be overlooked, this has shown that the use of the URANS model for close follow distances is ill-advised because of the failure to produce correct far field predictions. Note that uncertainties in the CFD predicted velocity field are less than 4% of the reference mean velocity, which is comparable to the measurement uncertainty estimates of 1% and 4% in free-stream and in the shear layer, respectively, as reported by Martinuzzi and Havel [40]. These considerations suggest that the CFD predicted velocity field, using the IDDES methodology, is actually very well correlated with the experimental measurements. With this in mind, the IDDES methodology will be used for the remaining simulations.
When comparing the force coefficient predictions, the Reynolds number dependence of the geometry must be taken into account. As such, a Reynolds number Re = 5610 was chosen in order to be in agreement with the DNS study performed by Curley and Uddin [43]. The results of the modified IDDES case compared to the DNS work of Curly and Uddin [43] are presented in Table 4. This comparison shows that all variation between force coefficients is kept below 7% which is believed to be acceptable for this work.     A second block was then placed into the simulation and meshed with the same schemes as described previously. This tandem cube simulation was designed to recreate the work of Martinuzzi and Havel [40,45]. This validation with the tandem cubes against Martinuzzi and Havel's work will allow for the validation of the simulation methodology which can then be applied to the DrivAer model simulations. To compare data with this work, four positions where chosen so that they were placed into the three regimes identified by Martinuzzi. The first is x/h < 1.5 in which vortex shedding is interrupted and random fluctuations are picked up on the spectra frequency. The lock-in regime is between 1.5 < x/h < 2.5 in which vortex shedding is continuous and scales linearly with s. The final regime is the bi-stable regime and is between 4 < x/h < 6 where two non-harmonic shedding frequencies appear. These four positions are shown in Figure 10 which displays the differences between the current simulations and the experimental work [45]. While some slight disagreements can be seen in the coefficient of drag (C D ) values, overall a good agreement can be seen between the two works, with only a slight over estimate in the lead block's drag being the most noteworthy difference. Note that the drag coefficient, and U ∞ represent the total drag force, frontal reference area, fluid density and reference free-stream velocity, respectively. The differences in drag are again explain by the differences in blockage ratios between the two works with the current having almost no effect from blockage. While Martinuzzi and Havel [40,45] didn't explicitly mentioned the uncertainty in drag measurements, based off of the given uncertainties in velocity measurements, we anticipate approximately ±5% uncertainty in experimental drag measurements; the uncertainty in the CFD predicted drag value is less than 1%. While these agreements help validate the current methodology, to completely validate the work, the mean and instantaneous flow fields around the two cubes were compared with the experimental work of [40]. The mean velocity vector fields on the planes x/h = 4 and x/h = 2 are shown in Figure 11. These show that two main recirculation bubbles occur, the first, and larger of the two, occurs behind the leading cube and the smaller occurs behind the trailing cube. When comparing the two cases, good agreement is seen between the large-scale flow features, such as the vortex core positions, with the only main disagreements occurring in the amount of separation over the top face of the leading cube which was already noted in the isolated case to be over predicted by the IDDES model. When analyzing the flow features, the rear cube's vortex is believed to be smaller than the leading cube due to the leading cube's wake curving streamlines onto the top of the rear cube, creating a flow analogous to a backward facing step. This caused the flow over the top of the rear cube to remain attached, unlike the front cube, and so reduced the size of the recirculation bubble. Since both the force coefficients and main flow field show good agreement in the cube validation cases, the DrivAer geometry was simulated using the same testing methodology. Again, the isolated DrivAer model was first be validated using the modified IDDES-SST model against the experimental data of Heft et al. [16]. It is important to note that a 40% scaled DrivAer model was used in line with Heft et al. in order to reduce computational costs. During the wind tunnel testing of Heft et al., restraining devices were used on the car and a correction to the force data must be made in order to compensate for these restraints. The testing procedure for the isolated case was identical to the platooning case but with the trailing car removed, both time step and Reynolds number remained the same through both cases. Table 5 presents the overall C D values presented in both the current work and the past work of Heft with a correction, as suggested by Collin et al. [46], to allow for the vehicle restraints used in the wind tunnel during the experimental work. The results of this comparison shows that the two differ by less than 3% which is well within the acceptable range or measurement and numerical prediction uncertainty. Figure 12 shows the C p over the top surface of the vehicle as compared to same experimental work of Heft et al. [16].
While the C p data shows good agreement, at x = 500 mm and x = 1200 mm which correspond to the lead and trailing edges of the roof respectively. These differences are most likely explained by the absences of a stinger in the CFD cases while one was present during the wind tunnel testing. This conjecture is also supported by Heft's observations that the sting significantly affected the flow on the rear glass for all three car configurations. To further validate the CFD modeling approach, the complex wake of the geometry is compared to an experimental work by Strangfeld et al. [47] in Figure 13. This cut plane shows two counter rotating vorticies from the bottom and top of the vehicle, compared using the orange and yellow lines. The current work shows a slight difference in position of the vortex core, with the CFD predicting a slightly higher position for the under-body vortex core and more rearward prediction for the upper-body vortex core. However, the saddle point of the vorticies are in good agreement with both occurring at approximately x/L = 0.17 for the rear of the vehicle.
The wake of the body was then compared to the numerical works of Soares et al. [19] and Ashton et al. [48,49]. Figure 14 compares the wake flow field of Soares et al. [19] at x/L = 1 to the current CFD work. This comparison of the total pressure loss shows that the current study differs from the previous work significantly only in the side flow near the wheels. The difference in the wake near the wheels can be attributed to how wheel motion was simulated. In the current study a moving reference frame was used while in Soares et al. [19] study only a tangential velocity was used, creating less of a wake effect from the wheels. This is further supported when comparing the areas of separation in the current work to that of Ashton et al. [48] in Figure 15 who also only used a tangential velocity to simulate wheel movement. From this it is seen, that much more separation is seen in the wheel areas of the current work than in the other works, thus creating the difference in pressure loss between the current work and Soares et al. [19]. Further comparing the separation regions Figure 15, it is seen that the current work slightly overestimates the length of the recirulation bubble as compared to Ashton's highly resolved IDDES SST simulation. This, however, does not present a large enough disagreement to discredit either simulation as the bulk trend of the flows are in agreement. So, it was noted that this high prediction of separation may influence the results of close drafting separations. Figure 13. Wake Flow of Isolated DrivAer Fastback Simulation versus experimental [47] with key circulation zones compared using dashed lines.  The second drivAer model was then added to the simulation and set up under the same meshing principals used for the tandem block cases. The IDDES model with modified coefficients was once again used for the simulation of flow. Figure 16 shows that at the closest vehicle spacings, both the lead and trailing cars experienced a drag reduction when compared to the isolated model. As the spacing increased, the drag increased for the leading car until the isolated value was achieved at around 2L. However, the trailing car's drag never achieved the isolated value but instead still experienced a drag reduction. The reduction in drag at closer spacing is hypothesized by Pagliarella [12] and Watkins [50] to be from the impingement of the lead car's vorticies by the trailing car. The vortices normally produce downwash and low pressure on the base of the car, when the trailing car approaches the leading car, it impinges on these vorticies and creates a higher pressure region in between the two vehicles. Figure 17 shows this vortex impingement occurring during the closest separation of the DrivAer models.
To better understand the why the drag of the trailing vehicle did not return to the isolated state like the lead vehicle, the normalized velocity was taken along the z = 0.15 plane and displayed in Figure 18. From this plane it is clear to see that while the leading vehicle has returned to a isolated flow case, the trailing vehicle is still being influence greatly by the wake of the leading vehicle. While this influence is slight compared to the smaller spacing, it does show that dirty air produced by the leading vehicle will have an effect on the trialing vehicle until the flow is able to relaminarize. For more complex geometries, the relaminarization process may take far longer than a simple body like the DrivAer due to the large amount of energy imparted into the wake which will cause the effects of the leading car to be seen well past the even the 8 x/L seen in the current study.   This impingement also causes the pressure on the front of the trailing car to decrease compared to its isolated state. Figures 19 and 20 show this decrease in pressure which occurs on the front of the vehicle. The decrease is most drastic from the x/L = 1 and decreases its effects until at x/L = 8 almost no difference is seen. This decrease in pressure poses interesting implications for the use of drafting in motorsports situations and long term heavy vehicle platooning. In motorsports, this decreased pressure may also lead to a loss in downforce for vehicles equipped with front downforce devices such as front splitters, and so create an adverse effect on the car's ability to brake and corner due to the reduced front downforce. In both HVP and motorsports this decreased pressure on the rear vehicle may also decrease cooling efficiency due to the decrease pressure on the nose of the vehicle where the majority of the cooling openings are located. This may raise the need for a switching of the leading and trailing vehicles in long term platoons, like those used in heavy transport vehicles, in order to not overheat the engine of the trailing vehicle. Optimal switching configurations may need to be developed in addition to such platoon formation algorithms proposed by He et al. [51] and others, see [5,7].

Conclusions
During this investigation, several main goals were achieved. The first was to investigate the prediction capabilities and potential improvement of the default URANS and IDDES models in highly turbulent shear flows, such as the wakes of surface mounted cubes and vehicles, through the use of closure coefficient tuning. It was found that the modified URANS was able to slightly outperform the IDDES model in the near wake of the surface mounted cube but failed completely to capture far field wake features, where the IDDES showed great agreement with the experimental works. This failure of the URANS model was believed to be a failure of the model to predict the relamarization of the flow. The greater accuracy in the far field then led the IDDES model to be chosen for the tandem case studies of both the block and DrivAer models. During the tandem block studies, it was seen that good agreement between the experimental work and this simulation methodology occurred both in the force coefficients and in prediction of the mean flow field. So, this methodology was continued to the DrivAer models where it was able to reproduce the same success by predicting the drag coefficient value to within 3% of the work performed by Heft [16]. The DrivAer study revealed several key insights into the drafting behavior of real-world automotive vehicles.
The specific conclusions drawn from this study are: 1.
The wake of the fully detailed models extends much father down stream than the simplified geometries such as the Ahmed body. This causes the effects lead vehicles wake to be felt by the trailing vehicle even past 4L, unlike in the Ahmed body studies.

2.
Unlike simplified geometries, the DrivAer trailing vehicle saw a reduction in drag at the closest separations instead of an increase. This is believed to be because the vortices of the leading vehicle were impinged upon by the front area of the trailing vehicle, decreasing the pressure on the front of the vehicle, reducing the drag. 3.
The pressure on the front of the trailing vehicle was seen to decrease significantly throughout all separations tested. This implies that for vehicles with downforce producing devices, such as splitters, a significant decreasing in front downforce may be experienced which following in the wake of another vehicle. In addition to a loss in downforce a reduction in cooling efficiency may also be seen due to the reduction in pressure on the front of the radiator.
This study has shown that is possible to predict, with reasonable accuracy, the complex interactions that occur in a fully detailed vehicle platooning situation using the setup and validation methodology presented in this paper. Future studies seek to use this methodology to investigate the flow of fully detailed vehicles during an overtake maneuver, in which the trailing vehicle has a velocity relative to the leading vehicle.

Data Availability Statement:
The data presented in this paper is available from the authors upon request.