Computational simulations of unsteady flow field and spray impingement on a simplified automotive geometry

Abstract Accurately predicting vehicle soiling is important for maintaining a clear view for the driver and on board camera and sensor systems. In this work we study the soiling process on a scale model of generic SUV body, which is a vehicle type particularly susceptible to base contamination. The Spalart-Allmaras formulation of the IDDES model is used to compute the continuous phase and the dispersed phase is computed using Lagrangian particle tracking, both concurrently with the flow-field, and also as a post-processing approach using time averaged statistical information of turbulence in a stochastic dispersion model. The results are compared against experimental data and the discrepancies discussed with regard to the predicted and measured flow field and base pressure distribution. Good agreement with experiment is shown for the contamination pattern using the fully unsteady method, but the more economic stochastic model does not recover some important details. This is attributed to the role of spatially correlated flow structures around the wheel in entraining particles into the wake that the stochastic model cannot accurately represent. This leads to the conclusion that base soiling is a function of unsteady modes, elimination of which may potentially reduce spray deposition.


Introduction
Management of surface contamination is an increasingly important development objective for vehicle manufacturers as it can reduce vehicle visibility, driver's vision, sensor and camera systems performance; as well as compromising aesthetics. With the development of autonomous vehicles, which rely on sensor systems that must be kept clean, the need to understand and predict vehicle soiling processes is likely to increase in importance (Kutila et al., 2015;Shearman et al., 1998). In addition, contaminants can be transferred to hands and clothing on contact with the vehicle exterior, which can be a source of customer dissatisfaction. Three useful reviews of this topic have been published to date. Kuthada and Cyr (2006) provide a short general review, before focusing on calibrating a tyre spray model and simulating body-side soiling. Hagemeier et al. (2011) explored side-window water management, providing a comprehensive review of numerical techniques. More recently, Gaylard et al. (2017b) have provided a comprehensive general overview. The three sources of vehicle surface contamination are: rain (primary contamination); spray generated by other vehicles (third-party contamination); spray generated by the rotation of the vehicle's wheels (self-contamination) (Kuthada and Cyr, 2006). The work reported here focuses on the latter phenomenon and is particularly relevant for vehicles with blunt rear geometries, such as Sports Utility Vehicles (SUVs), estates (Station wagons), hatchbacks (Maycock, 1966), as well as bus bodies (Lajos et al., 1986). This is due to strong large-scale recirculating vortices generated by these vehicles that advect the spray and deposit it onto the base (Gaylard and Duncan, 2011).
Concerns over vehicle rear soiling had become evident by the mid-1960s, when Dawley (1965) had tested deflecting vanes at the rear corners of the vehicle body to provide clean air into the wake and thus minimize dirt impinging on the rear surface. It was later realized that applying these to the roof trailing edge, as proposed by Goetz (1971), was most effective for reducing base soiling. Nevertheless, this was achieved at the expense of increased overall drag. A deployable roof trailing edge deflecting vane was also studied by Costelli (1984), who linked the mechanism of rear-surface soiling with the pressure variation in the wake. He postulated that soiling tends to accumulate in regions of relatively high surface static pressure. However, recent simulations of rear face soiling for a relatively simple bluff body showed this to be a necessary though not sufficient condition and that local contaminant availability is also required (Jilesen et al., 2017). The work of Gaylard and Duncan (2011), Gaylard et al. (2014) and Jilesen et al. (2013) presented a comprehensive analysis of experimental and numerical base soiling results for a number of real vehicles. In the numerical approach, particles were emitted from the rotating wheels using the technique proposed by Kuthada and Cyr (2006). These studies showed that the main source of spray that gets deposited on the rear surface are rear wheels. In addition, the wakes of rear wheels contribute to the advection of droplets into the wake, which then transports them onto the base. It was also found that base contamination was better predicted using the fully unsteady particle tracking approach in Jilesen et al. (2013) rather than the one in which particles were post-processed for a number of frozen transient data frames and then superimposed to obtain the overall deposition, as was used in Gaylard and Duncan (2011)).
Most reported studies used detailed geometries of specific vehicles along with complex computational models, mainly focusing on the final outcome with regard to surface contamination. However, such an approach may not always be useful if one aims to understand the fundamentals of this phenomenon. This is due to the inherent complexity of the flow associated with geometric features and general level of detail of the physical model used. According to Jilesen et al. (2013), even small changes in the details of the vehicles' geometry can strongly influence contamination mechanisms. Although lacking the details and styling of specific vehicles, simplified bodies allow investigation of the relevant flow features in well characterised and repeatable settings, giving better insight into the studied phenomenon. In addition, they can be used to reveal trends due to design changes such as rear body shaping. The use of simplified, generic, vehicle bodies also allows soiling studies to be undertaken for cases for which detailed flow field data such as wake Particle Image Velocimetry (PIV) or base pressures is also known and publicly available. This, along with control of factors such as details of the contaminant source, allows the accuracy of simulation methods and the importance of modelling choices to be correctly assessed. One of the few examples of using a simple body in the studies of vehicle contamination is the study by Kabanovs et al. (2016), which used Computational Fluid Dynamics (CFD) and experiments to investigate rear surface contamination for a well known simple bluff body. This work showed that neither RANS nor URANS turbulence modelling is able to capture the relevant flow structures in the region of a strong pressure gradient, failing to accurately predict the flow field. This, therefore, resulted in an inaccurate prediction of spray dynamics and rear-surface contamination. High-fidelity eddy resolving methods have been used to produce encouraging predictions of the rear soiling of simple bluff body (Gaylard et al., 2017a;Kabanovs et al., 2017), the latter correctly predicting the trend in soiling pattern when the upper taper angle of the rear of a generic SUV geometry is changed.
In this work we apply an eddy resolving simulation technique, namely the Detached Eddy Simulation (DES) method also used in Kabanovs et al. (2017)), to simulate the rear soiling of two configurations of the generic SUV and compare with detailed PIV measurements and rear soiling experiments. It has been shown in previous work (Gaylard et al., 2017a;Kabanovs et al., 2016) that a high fidelity eddy resolving method is required to correctly predict the mean flow field and that this is a prerequisite to accurate soiling predictions. However what is less clear is whether, given an accurate mean flow field, is this sufficient to allow accurate soiling predictions? i.e. Can soiling be predicted as a post-processing step if a suitable set of data is available? Hence, in this paper we also investigate the level of detail required in the simulation of the dispersed phase by comparing: A fully unsteady, time-resolved, spray simulation with post-processed spray tracking simulations using only the mean flow field. A fully unsteady, time-resolved, spray simulation with post-processed spray tracking simulations using the mean flow field with stochastic dispersion due to the time averaged turbulence field. The time averaged statistical information of the turbulence is available from the DES simulation and is used by a stochastic dispersion algorithm that models velocity perturbation seen by a particle advected through a time averaged flow field.
These post-processing methods represent a significant saving of computational time. The simulations are also used to investigate the mechanisms by which contaminant is entrained from wheel spray into the wake and then dispersed and distributed onto the rear surfaces of the vehicle. The influence of vehicle geometry on these mechanisms is investigated, both from the point of view of how this affects contamination, and how they should be modelled. Discrepancies between experiment and simulation are discussed together with suggestions for further improvements to the current methodology.

Experimental method and test geometries
A quarter scale generic SUV body was used in this study. The model is representative of a typical SUV and was designed within the department of Aeronautical and Automotive engineering at Loughborough University. The aft body has sharp edges while the leading edge radii are large to prevent separation. The model has configurable elements such as ride heights and the aft body. Four pins, attached to stationary wheels connect the model to the balance in the wind tunnel. The model gives a blockage ratio of 5.58% in a 2.5 m 2 working section. More information can be found in Wood et al. (2015). In this work, two configurations of this model were used. Both configurations had the ride height fixed at 0.065 m and the roof taper angle set to 0 . The angle of the diffuser, however, was set to 30 (configuration 1) and 0 (configuration 2). These two configurations were chosen on purpose, as they are recognized to produce very different wake structures and hence were expected to give unique soiling distributions. On the other side, geometries with bluff and square aft bodies such as configuration 2 are known to be very challenging to simulate, which is associated with the unsteady, sometimes bi-stable, wake of these models. For example, a number of experimental studies reported the presence of large scale (both in time and in size) lateral oscillations of the entire wake (flapping mode), accompanied by a "breathing" mode of the mean recirculation region of the generic square-back SUV (Al-Garni et al., 2004), and other boxy geometries (Perry et al., 2016;Khalighi et al., 2001). Both configurations are illustrated in Fig. 1. It also shows the location of spray injector used in soiling experiments, discussed below.
Aerodynamics and soiling experiments were performed in the wind tunnel at Loughborough University. The wind tunnel has an open loop, circuit, closed working section configuration (see Fig. 2), full details of which can be found in Johl et al. (2004). All tests were carried out at a free-stream velocity of 40 m/s, giving a Reynolds number of 2.77 million based on model length. Full details of the experimental method, especially that used in soiling tests, can be found in Kabanovs et al. (2017) but important details are included in the two following sub-sections.

Aerodynamics tests
The experimental force and base pressure data used in this paper were collected for each configuration as part of two separate studies, reported by Wood et al. (2015) (configuration 1) and Varney et al. (2017) (configuration 2). The data taken for configurations 1 and 2 were sampled for 30 s and 300 s at a frequency of 300 Hz, respectively. The base pressure tappings were limited to one half of the base in the experiment that used configuration 1. The entire base of configuration 2 was populated with the pressure tappings due to a highly unsteady nature of the wake of this configuration. The pressure measurements were made with two 64 channel pressure scanners with samples triggered at 260 Hz. The pressure data and force coefficients presented in this paper account for blockage by applying a continuity based correction method shown in Eqs. (1) and (2), respectively (Littlewood and Passmore, 2010). In the equations, TA and MA correspond to the wind tunnel cross section area and the model frontal area, respectively. C p and C f are the pressure coefficient and force coefficient, respectively.
The PIV data for configuration 1 was collected and presented by Wood et al. (2015). Although the data for configuration 2 was also obtained in the same study, it was not included in any previous publications. For details of the setup readers are referred to Wood et al. (2015). In summary, the measurements were collected with a 4 mega-pixel, 14 bit camera equipped with a 50 mm lens in combination with a frequency doubled Nd:YAG Litron LASER with 200 mJ pulse energy. Data collection and subsequent analysis where performed using the LaVision commercial software. The image spatial resolution was 0.2 mm/pixel based on 400Â400 image area. A total of 1000 image pairs were captured for each plane at a frequency of 7.26 Hz.

Soiling tests
In soiling tests, the seeding setup featured a single nozzle installed behind the left rear wheel and directed 45 downstream. Such arrangement was chosen to reasonably represent tyre spray. A deliberate choice was made to use a single source rather than one located behind each rear wheel because such an approach ensures that a clear analysis of the source and outcome is possible. Tests with the first configuration used a seeder placed 0.086 m downstream of the rear wheel centreline, while a distance of 0.16 m was used in the tests of the seconds configuration. The position of the seeder is shown in Fig. 1. The inconsistent position of the nozzle between two tests is associated with a large cone angle and hence the necessity to move it further downstream to stop the contaminant impinging on the under-body of the second configuration and subsequently being stripped back into the wake. These processes are not considered in the CFD model and they are unlikely to be reproduced on a full scale automotive geometry. Hence it was decided to minimize the impingement of spray on the under-body by moving the contaminant further downstream.
The contaminant was a mixture of de-ionised water and UV dye (Uvitex at 0.03% concentration), injected at 110 bar, giving a mass flow rate of _ m ¼ 3:48g=s. The model was run in the wind tunnel for some time to settle prior to injecting the contaminant and each soiling test was 30 s in duration. Following every test, the UV lamp and a DSLR camera were placed in the wind tunnel behind the model and an image of the rearsurface captured. The photographs were then processed based on the fluorescence to produce soiling intensity plots, used later to validate CFD results.
To ensure simulations used accurate spray parameters, a separate experiment was performed employing a Phase Doppler Anemometry (PDA) system to provide statistical information with regard to droplet population and their velocities in the spray. More information about the PDA system used can be found in Jiang et al. (2016) and Wigley et al. (1999). More details of the PDA testing of the spray used in this work can be found in Kabanovs et al. (2017).

Modelling of the continuous phase
Numerical simulations were performed using OpenFOAM ® (version  , an open source CFD tool. The continuous phase was calculated using the Improved Delayed DES (IDDES) (Shur et al., 2008) method, that is designed to treat the entire boundary layer using a RANS method and apply LES to separated regions. The advantage of the IDDES formulation over the original DES (Spalart et al., 1997) and DDES (Spalart et al., 2006) formulations is that the transition to resolved turbulence near the wall is quicker and it has an enhanced capability to identify which mode (RANS or LES) it should be operating in near the wall (Shur et al., 2008). This is achieved by using a universal sub-grid length scale formulation that satisfies the nature of flows in the vicinity of walls and further away from wall surfaces: In the equation, d w is the distance to the wall, Δ wn is the grid space in the wall-normal direction, Δ max is the maximum local grid spacing, and C w is a constant equal to 0.15 based on LES of a developed channel flow. In this work, the RANS branch of the IDDES model used the Spalart-Allmaras turbulence model; an eddy viscosity model that follows the Boussinesq relationship to close the Reynolds stresses and provides a single equation to derive eddy viscosity (Spalart and Allmaras, 1992).
The near-wall treatment uses a "universal" velocity profile suggested by Spalding (1961). The advantage of this formulation is that it allows a wider range of wall-normal computational grid spacing as it provides the velocity profile that includes the viscous sub-layer, buffer layer and log-layer. It is usually referred to as an all-y þ wall treatment and has been applied in other studies (Joubert et al., 2015).
The PISO (Pressure Implicit with Splitting of Operators) algorithm was used to solve the Navier-Stokes equations. A constant global timestep Δt ¼ 2:5 Â 10 À5 s was used throughout the simulation, which kept the mean Courant Friedrichs Lewy (CFL) number (CFL ¼ UΔt=Δx) at 0.035, allowing stable simulations. The maximum CFL number was 8 for a few cells located in a 4 mm gap between the floor and the base of the wheels. This is believed to be due to a significant acceleration of the flow in that region. The CFL number was below 1 in the detached regions, ensuring that turbulent eddies were resolved in space and time. Second order numerical schemes were generally used.

Modelling of the dispersed phase
The dispersed phase was computed using the Lagrangian approach (Migdal and Agosta, 1967). This technique has been widely used to study various multi-phase problems, including the problems associated with atmospheric pollutants (Ahmadi and Li, 2000), sand (Paz et al., 2015), wind-driven rain (Wu et al., 2017;Persoon et al., 2008), and pesticide spray (Xu et al., 1998). A total of 17.5 million parcels per second were released, each parcel representing a number of particles of specific size and mass. The forces of main significance were considered to be the drag (F D ), gravity (F G ) and shear lift (F L ) (also called Saffman-Mei lift (Saffman, 1965), see Eq. (3)). The shear lift is particularly important in the regions of strong shear (Barton, 1995).
The drag C D was modelled using Liu correlation, which is based on the assumption of solid spheres (Liu et al., 1993) and is shown in Eq. (4).
where C LS is calculated as follows:
The influence of the presence of particles on the flow structures can be estimated using a momentum coupling parameter Π, given in Eq. (10). The momentum coupling parameter Π represents the ratio of particle drag to the carrier fluid momentum flux and is greater than unity in regions where the back-coupling becomes important (Paschkewitz, 2006a). In the equation, the parameter C corresponds to the ratio of the mass flow rates of the dispersed and continuous phases and St (given in Eq. (11)) quantifies the relative importance of particle inertia.
Similarly, the likelihood of liquid particles to break-up can be estimated by calculating the ratio of inertia forces to the surface tension of a particle (Weber number). The definition of Weber number is given in Eq. (12), where U slip , ρ f and σ p are the slip velocity, the density of the fluid and the surface tension of the liquid (water), respectively. According to van Basshuysen and Schafer (2005), the secondary break-up mechanism should be considered when the Weber number exceeds 10. We To assess the effects of the particles on the flow structures and the possibility of particle break-up, two images that show instantaneous results for configuration 1 are presented in Fig. 3. The spray is coloured by the coupling momentum parameter Π (a) and the Weber number (b). It can be seen that there are very few regions with the parameter Π exceeding 1. In addition, the chance of particle break-up is very unlikely due to the small Weber number. As a result, neither back-coupling nor the break-up of particles was included in this study. Weber number and coupling parameter were also found to be similarly low for configuration 2. Although the collision of particles and particle evaporation may play an important role, these were not considered in the current study but may be of interest in future research.

Injection model
The particles were introduced into the computational domain at the same location as in experiments using ten cone injectors assembled together to account for droplet size distribution which changes with lateral distance from the centre of the spray. This was based on experimental results obtained using the PDA system mentioned previously. The ensemble-averaged velocity measured at each lateral position was used as the initial velocity of computational particles leaving each injector. The total cone angle was 90 . The spray consists predominantly of lowmomentum particles. The full details of particle size distribution, as well as a detailed description of the injection approach used in the numerical model can be found in Kabanovs et al. (2017), as it was essentially identical.

Analysis methods
In order to assess the level of detail needed to capture the unsteady transport of the dispersed phase three methods of particle tracking were used. In the first the particle tracking is performed each time step as part of the IDDES simulation. This method allows full details of the turbulent transport to be included. The second and third methods use time averaged data from the IDDES simulation and tracks particles as a postprocessing step. This would allow a saving of computational time over the fully unsteady method, particularly in a vehicle development process where the data will have been produced as part of the aerodynamic design work. The three methods are further described below.
3.2.2.1. Fully unsteady simulations. The first method used a fully unsteady simulation, along with the Lagrangian particles tracked concurrently with the flow solver. In addition to the numerical settings explained previously, a stochastic dispersion model was applied that is commonly used to account for unresolved turbulent scales and their influence on the velocity of computational particles. Even though the fraction of the unresolved content in DES computations tends to be relatively small, these unresolved structures may play an important role in a particle dispersion mechanism and therefore should be considered (Pesmazoglou et al., 2013). In this work, the formulation of the Gosman and Loannides (1983) stochastic dispersion model was used, which is essentially a discrete random walk (DRW) model. This approach adds a Gaussian-distributed random perturbation u 0 (Eq. (13)) onto the resolved air velocity, seen by the particle and used in the drag model (Eq. (3)) to mimic the sub-grid turbulent motion. The characteristic eddy lifetime t eddy % k sgs =ε sgs determines the period over which the random value of u 0 is kept constant, as it is not necessarily updated each time step. This allows both velocity and timescale information to be included in the dispersion of the particles. In the model, ζ is a normally distributed random number, d represents spatial randomness of turbulence and σ ¼ ffiffiffiffiffiffiffiffiffiffi ffi 2k=3 p is the standard deviation, which is related to the turbulent kinetic energy.
here, the stochastic dispersion model uses sub-grid turbulent kinetic energy (k ¼ k sgs ) and the sub-grid turbulent dissipation rate to account for unresolved scales. The k sgs is supplied by the one-equation sub-grid scale LES model (Ni ceno et al., 2008;Davidson, 1997) via Eq. (14), while the dissipation rate ε sgs is computed using Eq. (15). Here, ν t is the sub-grid turbulent viscosity, Δ f is the filter size, C k is an empirical constant that is equal to 0.07, ν eff is the effective viscosity (ν eff ¼ ν þ ν t ) and S ij corresponds to the strain rate tensor.
The effect of the stochastic dispersion model was studied by running simulations with and without this model activated. The use of the stochastic dispersion model was found to increase soiling rate by 5%, without affecting the pattern of contamination.
3.2.2.2. Particle tracking in a mean flow-field. To better analyse the role of flow unsteadiness on the dynamics of vehicle soiling, particles have been tracked in a mean flow-field without any unsteadiness taken into account. These simulations used the accurate mean velocity field computed using the fully unsteady computations.
3.2.2.3. Lagrangian particle tracking in a mean flow field with an isotropic dispersion model used to replicate flow unsteadiness. Due to the capability of high fidelity CFD methods, such as IDDES, to resolve a large percentage of turbulent scales, good statistical information of turbulence can be obtained. This statistical information can be used to enhance the particle tracking through the time-mean velocity field. The variance of velocity perturbations in all three directions, recorded from a fully unsteady simulation, is used here to construct the mean resolved turbulent kinetic energy field (see Eq. (16)).
This field, along with the mean sub-grid turbulent kinetic energy field, is then fed to a stochastic dispersion model (Eq. (13)) to compute a perturbation velocity vector for each parcel, based on the total mean turbulent kinetic energy (k ¼ k sgs þ k res ). To accurately estimate the characteristic time scale of eddies, a model-based estimation of the sub- grid turbulent dissipation rate, proposed in Moeng and Wyngaard (1988) Eq. (17), is used.
In the equation, Δ is the grid spacing and C ε ¼ 1:048 is a constant used in the Smagorinsky sub-grid model. This estimation is based on the assumption that eddies dissipate at smallest scales. This method tracks the dispersed phase advecting through a time averaged flow field (calculated beforehand using a fully unsteady simulation) with the turbulent stochastic dispersion switched on to impose the artificial unsteadiness of particles. Results from this method are described as being from a stochastic modelling approach in subsequent sections of this paper.

Computational domain and boundary conditions
The computational meshes were generated using snappyHexMesh. This meshing algorithm uses hexahedra and split-hexahedra cells to construct grids and is available within the OpenFOAM ® software suit. The particle tracking approach available within OpenFOAM ® makes use of a tetrahedral decomposition of mesh cells that is compatible with meshes generated with snappyHexMesh.    4.7 m ahead of the model. This length is based on the distance required for the boundary layer to reach the thickness determined in empty-tunnel tests. The middle sub-domain has divergent walls in order to match the physical wind tunnel arrangement that is designed to ensure that the longitudinal pressure gradient is zero when the domain is empty. Table 1 provides further dimensional information of the computational domain. There are prism layers around the model and refinement zones, particularly in the rear of the model, to capture strong gradients in the flow. The thickness of the first cell in computational meshes used in this study was such that the mean yþ was around 30 (high yþ). This complies with the all-y þ near-wall treatment (Spalding, 1961), used in this work. All meshes used the same topology and contained approximately 63 million cells. A very similar computational procedure had been used by Forbes et al. (2017), who studied the aerodynamics of the generic SUV and achieved good agreement with experiment in the same wind tunnel arrangement.

Results and discussion
Computations have been run on 2.80 GHz (20 CPU per node) and 64 GB RAM nodes. An initial single-phase simulation was run for each configuration to establish the flow field, requiring 32016.8 CPU hours of computational effort to compute 1 s. The unsteady soiling simulations were then run for 2 s of simulated time each, requiring additional 90912 CPU hours per simulation. The time averaged flow field data and turbulence statistics were also obtained in these simulations. These were later employed in computations that used stochastic modelling of velocity perturbation and required 7047.34 CPU hours to compute 2 s worth of particle tracking. Figs. 6 and 7 present the time averaged velocity streamlines in the wake of the same SUV configuration, projected on vertical and horizontal planes, respectively. The PIV data is also shown. The correlation between the experimental PIV data and numerical results is very good, with a slight mismatch in the centre of the upper vortex on the vertical centreplane. This is believed to be associated with the flow upwash, which is stronger in the experimental results, thus pushing the upper vortex closer to the base. Fig. 8 shows the time averaged pressure distribution on the base of configuration 2. Although the base pressure distribution is matched well, the values predicted by the computational model are somewhat high. Nevertheless, Figs. 9 and 10 show a good agreement of flow structures between CFD and PIV, but with some difference in details. The CFD predicts a stronger upper vortex on the centre plane ( Fig. 9 (a) and (d)) leading to a greater upwash towards the upper part of the centre of the base, consistent with the higher pressure on this part of the base.

Aerodynamics results
Clear asymmetry in the experimental data is seen in (a) of Fig. 10. It should be noted that, due to its boxy geometry, this configuration of the generic SUV appears to be very sensitive to yaw angle. It is possible that the model had been at a slight yaw during the experiments leading to a mismatch between the location of the centre of the bottom vortex seen in (a) and (d) of Fig. 9. It is clear that there is also asymmetry in the CFD results, particularly visible in Fig. 8, although not to the extent seen in the PIV. It was found that by laterally moving the model geometry in the CFD fractionally ($ 0:05% of vehicle width) within the tunnel that an  asymmetric wake biased the other way could be predicted. This emphasizes the sensitivity of the wake of this flow to any asymmetry, and also demonstrates the challenging nature of this flow for CFD. This may explain why the prediction is less accurate for this case than for configuration 1, in which the wake is dominated by the presence of the diffuser, making it much more laterally stable. Table 2 presents experimental and CFD lift and drag coefficients. Although the base pressures match very well, CFD predicts a higher drag for configuration 1. However, the excellent agreement with base pressure and PIV gives us confidence that the wake is being predicted accurately and that the discrepancy in drag is due to other sources such as the modelling of the attached boundary layer. In addition, the predicted lift is negative, in disagreement with the experiment. This has been also reported in Forbes et al. (2017), who suggests that the source of this is possibly due to inaccurate modelling of the interaction of the under-body flow and the boundary layer on the floor. In addition, an annular gap around the balance pins is present in the experimental setup. It is fair to assume that the leakage flow through the annular gap affects the flow between the truncated wheel and floor, thus affecting pressure and lift. In the CFD setup the annular gap around pins is not present. The drag predicted for configuration 2 is in good agreement with the experiment. However, similar to configuration 1, there is a mismatch in the lift coefficient.
Finally, Figs. 11 and 12 show the turbulent intensity (I ¼ ffiffiffiffiffiffiffiffiffiffi ffi 2=3k p =U ref ) and the resolved turbulent content on the symmetry plane for both configurations, respectively. In Fig. 11, the total intensity is calculated based on the sum of sub-grid and resolved turbulent kinetic energies. Sub-plots (c) and (d) of Fig. 11 indicate the total level of unsteadiness, while (a) and (b) show the modelled content. According to Pope (2000), the filter and grid are sufficiently fine if the resolved turbulent content is greater than 80%. The resolved turbulent content in the wake, shown in Fig. 12 is high and is in the range of 70% À 95%, giving confidence that the unsteady flow is being modelled correctly. There appears to be a discontinuity in the flow field in Figs. 11 and 12 at the interface between the two different mesh resolution zones. This can only be found in the results that include sub-grid turbulence statistics and is associated with the way the sub-grid turbulence is estimated. As can be seen from Eq. (14), the k sgs is estimated using the filter size Δ f . The computational cells located in the transition region between the two different mesh resolution zones are irregular in shape, which may have caused discontinuities in the sub-grid results. This is believed to have no effect on the wake structure as the Spalart-Allmaras model does not actually consider k sgs . There may be an insignificant effect on the results obtained using the stochastic modelling approach.

Distribution pattern
The base soiling obtained with particle tracking using only the mean velocity field was found to be negligible and is not shown here, but will be discussed later. This immediately indicates that the flow unsteadiness must be included in some way even if an accurate mean velocity field is available. Figs. 13 and 14 present soiling results obtained experimentally and numerically for configurations 1 and 2, respectively. The numerical results shown have been obtained using fully unsteady simulations, as well as using computations in which the random motion of numerical particles was estimated with the stochastic model, explained previously. The numerical time-dependant point data was averaged across the base surface using cells of 12 mm to make numerical results comparable to the experimental ones, obtained over a longer time period. The experimental results do not give a quantitative measurement of the mass deposited on the rear surface, instead they give the soiling pattern and a relative soiling intensity for each case. As such, soiling results are presented coloured by the local soiling intensity. For the experimental results this is A. Kabanovs et al. Journal of Wind Engineering & Industrial Aerodynamics 171 (2017) 178-195 normalised by maximum UV intensity whereas for CFD it is normalised by maximum predicted deposited mass. This allows identification of the most contaminated areas on the base; however, the quantitative information is lost, since the results are normalised by the local maximum. Therefore, soiling patterns are compared here and soiling rates are discussed later. Experimental results obtained for configuration 1 and presented in (a) of Fig. 13 show a small contamination region slightly inwards of the source, in an approximately semi-circular pattern. This is replicated very well with the fully unsteady CFD (see (b) of Fig. 13). Moreover, the most contaminated regions, such as the upper edge of the base and the circular area below, are well predicted. What is not captured is a small contamination region at the top of the vehicle base in the central region. Possibly, particles should be injected over a longer time period to make that region become visible in the CFD results. Also in fully unsteady CFD results, an area of high soiling intensity can be seen stretching along the left edge of the base in (b) of Fig. 13. This is due to a large number of small, low-inertia, particles that are predicted to be deposited along this edge. These particles generally have a small response time to the changes in the flow, following the streamlines and being unable to exit flow recirculation. In the experiment it is possible that these small particles fail to reach the base either because they have evaporated or coalesced into larger droplets. A further study of the effect of including models for these processes may provide further insight. Fig. 13 (c) shows that the stochastic model predicts the general location of particle impingement reasonably well, however, the areas of highest intensities are not matched and the pattern is shifted towards the left edge of the base. In addition, this method predicts a much stronger impingement of particles on the bottom left part of the base than in the experiment. Fig. 14 shows the experimental and predicted soiling patterns on the base of configuration 2. The change from configuration 1 to 2 leads to the contamination to be located more centrally on the base and dispersed over a wider region. This is broadly captured by both the fully unsteady and the stochastic modelling approaches. In the experimental results the centre of the contaminated area for configuration 2 is slightly to the opposite side of the source, as can be seen in (a) of Fig. 14. It covers a significant portion of the vehicle base with the highest intensity contamination near mid-height. This contradicts Costelli (1984), who suggests that particles tend to deposit in areas of higher base pressure. However it should be noted that the base pressure for this case is   . 11. DES-performance analysis; Top to bottom: Mean sub-grid intensity, mean total intensity; Configuration 1 on the left, configuration 2 on the right; Symmetry plane shown.  relatively uniform compared to configuration 1. While the CFD results correctly predict a wider distribution of soiling, both methods put the highest intensity soiling too far up the base and more centrally than in experiment. Two explanations for this are suggested. Firstly, CFD is unable to capture the unsteady wake of this configuration to the extent required to provide accurate soiling distribution pattern. For example, there is a stronger upwash towards the top centre of the base in the CFD than in the experiment (see Fig. 9) and this carries contaminant higher up the base. In addition, large scale flapping may be present in the dynamics of the wake that is not captured by CFD. This would explain the contaminant being deposited on the other side of the source. Secondly, evaporation and coalescence may change the size of the particles in the wake. Although the location of highest soiling intensity does not agree particularly well between the experiment and the unsteady simulation, the match in the spread of liquid on the base is apparent. It is pertinent to note that the turbulence levels seen in Fig. 11 are broadly similar for configuration 1 and 2. Therefore the increased dispersion of contamination seen in configuration 2 is not due to higher turbulence in the wake but rather due to the increased residence time of the spray in the wake of configuration 2. This is explored further in the next section. The numerical results obtained with the stochastic model shown in (c) of Fig. 14, although less dispersed, agree reasonably well with the first numerical approach. The broad trend of changed soiling from configuration 1 to 2 is also captured with the stochastic modelling post-processing approach. This suggests that this method could be used to gain some understanding about likely soiling trends using pre-existing simulations run to generate aerodynamic data or potentially using experimental data if both mean and fluctuating velocities were available. However to produce accurate simulations of soiling patterns it is clear that the full details of the unsteady particle transport must be included through the use of concurrent particle tracking as part of an unsteady eddy-resolving simulation. The reasons for this, together with the mechanisms which transport the contamination are discussed in the following sections. Fig. 15 shows instantaneous spray tracked in a mean flow-field for both configurations with no unsteadiness modelled. It can be seen that the particles released behind the rear wheel of configuration 1 ((a) of Fig. 15) instantly lose their momentum due to a steady downwash from the wheel arch. As a result, they are spread across the floor rather than entrained into the wake. A small fraction of particles released behind configuration 2 are picked up by the wake and deposited on the base in a non-physical manner, as can be seen in (b) of Fig. 15.

Deposition mechanism
In contrast, Fig. 16 shows iso-surfaces of the time-averaged spray volume fraction together with the mean velocity streamlines for two configurations obtained using the fully unsteady and stochastic methods. The spray is coloured by the mean particle size. A large fraction of spray  A. Kabanovs et al. Journal of Wind Engineering & Industrial Aerodynamics 171 (2017) 178-195 inside the wake region compared to that seen with the mean field shows that flow unsteadiness plays a crucial role in transporting particles into the wake. It can be seen in (b) and (d) of Fig. 16 that the stochastic model is able to simulate some particle entrainment into the wake, though the fraction of spray that is entrained is smaller than that computed with a fully unsteady simulation. Due to the upwash-dominated smaller wake of configuration 1, particles are drawn earlier into the wake than the particles behind configuration 2. This also explains the greater dispersion of the soiling on the base of configuration 2. The particles are entrained into the rear of the longer wake and are dispersed more widely as they are transported upstream towards the base, is clearly illustrated in Figs. 17 and 18 that show the paths of a sub-set of the parcels that end up on the base for configurations 1 and 2, respectively. The paths are coloured by particle size and the effect of the longer wake length can be seen for configuration 2. It is evident that the particles released behind configuration 1 are more disturbed in the vicinity of the wheel wake than those released behind configuration 2. This is most likely due to a complex flow interaction of the wheel wake and the flow from the wheel arch near the injector. It can be seen that the entrainment of particles into the wake is different for the two numerical approaches, particularly for configuration 1. However, once entrained into the wake, they are transported towards the base in a similar manner. This suggests that large scale unsteady flow structures are responsible for entrainment of particles into the wake, while their dynamics in the wake can be adequately described by the mean flow-field and a stochastic description of the turbulence.
Differences can also be seen in the size of the particles predicted to be deposited on the base with the two numerical methods. Figs. 17 and 18 show that the stochastic modelling approach predicts a very large fraction of small particles in the wake region. These low-momentum particles have a small response time and are easily influenced by any velocity perturbation. In addition, they may become trapped in a recirculating flow, which is possibly the reason for high contamination in the lower region of the base, seen in (c) of Fig. 13, and in the left region of the base, seen in (b) of Fig. 13. It is possible that coalescence and evaporation models may be important, as these models are expected to reduce the amount of small particles in the spray. Fig. 19 presents the mass-weighted size distribution of particles that have deposited on the rear-surfaces of configurations 1 and 2, modelled with two numerical approaches. The mass-weighted particle size distribution in the spray is also shown. This parameter provides statistical information about the fraction of mass that is carried by particles of specific sizes. It can be seen that the massweighted particle size distribution on the base of configuration 2 modelled with the stochastic approach is very close to the mass-weighted particles size distribution in the injected spray. Although less apparent, the stochastic modelling approach used for configuration 1 also predicts a larger fraction of small particles on the base than the fully unsteady method. This can also be seen in Figs. 17 and 18, which show that the fully unsteady method predicts a higher proportion of large particles in the wake and a greater vertical spread on the base. These large particles also seem to have a higher concentration towards the top of the base. This has been also reported in Paschkewitz (2006aPaschkewitz ( , 2006b, who studied Fig. 16. Three iso-surfaces of mean water volume ratio (3:0 Â 10 À7 ; 1:0 Â 10 À6 ; 1:0 Â 10 À5 ) coloured by the mean particle diameter, shown with streamlines on the Y ¼ 0.09 m plane for configuration 1 (top), and on the Y ¼ 0 m plane for configuration 2 (bottom); Fully unsteady simulation on the left, stochastic modelling approach on the right. particle dispersion through the wake of commercial vehicles and demonstrated that the vertical distribution of the spray is significantly different for steady-state and unsteady flow fields. As the large particles have a greater mass this also has an impact on the total soiling. Fig. 20 shows the predicted evolution of total contamination on the base of the two configurations. The rate of soiling modelled with the stochastic method is approximately three times smaller, this will be in part due to the greater number of large particles entrained in the fully unsteady method. In addition, the stochastic modelling approach predicts a very steady contamination rate, which is not the case in the fully unsteady simulations as they show the presence of multiple dynamics within the process of contamination. This is particularly true for configuration 1, while the contamination for configuration 2 is somewhat more linear. A similar transiency in the results was reported in Kabanovs et al. (2017), where similar experimental and numerical approaches were applied and same physical model was used, yet configured differently. The level of unsteadiness in the numerical results was stronger, possibly due to a different model arrangement. Interestingly, the total mass computed by the unsteady simulations is very similar for both configurations. However, as has been pointed out in Kabanovs et al. (2017), the apparent transiency in the contamination process makes the comparisons of soiling rates between experimental and numerical tests difficult, unless averaged over a sufficient amount of time. This contradicts Gaylard et al. (2017a) and Jilesen et al. (2017), who found that accumulation of contaminant on the rear surface is linear with time for a different geometry.

Rate of deposition and unsteady behaviour of spray
The unsteady process of contamination seen in the results obtained for configuration 1 and shown previously can be explained using Fig. 21. This figure presents the instantaneous spray behind the model, created with a volume fraction value of 3 Â 10 À6 and coloured by the average droplet size in the cell. The spray is captured in the figure at a frequency of 125 Hz. According to the figure, a fraction of the spray is pulled towards the diffuser region as a body. Then, the high momentum underbody flow pushes it into the wake region, where the cluster of particles is dispersed and a portion of it is carried by the recirculating motion of the wake and deposited on the base. To account for such an unsteady mechanism, consideration of the full unsteadiness of the wake is crucial as this large scale correlated flow structure cannot be fully captured by a stochastic description of turbulence. It is believed that these spatially correlated flow structures are responsible for entraining large particles  into the wake region, leading to a higher concentration of larger particles on the base in the concurrent eddy-resolving simulations, as seen in Figs. 17 and 18. The velocity fluctuations modelled by the stochastic approach are random and not spatially correlated, thus influencing the motion (drag forces) of large particles to a smaller extent than that of smaller particles. Some unsteadiness in the dynamics of the spray can be also seen behind configuration 2, as shown in Fig. 22. These flow structures, which the stochastic model can not simulate correctly, are believed to be responsible for entrainment of particles into the wake. However, the effect of these structures on contamination is less pronounced than in the case of configuration 1, probably due to a longer wake which smooths out deposition rate. Although this needs to be further studied, it appears that the unsteady behaviour in the contamination mechanism seen in the results is governed by the complex flow interaction of the wheel counterrotating vortices and the downwash from the wheel arch. This effect is likely to be stronger with the injector closer to the wheel. This also suggests that unsteady flow structures around the wheel wake can be responsible for increasing rear soiling compared to a steady turbulent shear layer behind the wheel. In this case the wheel is not rotating but a true rotating wheel is likely to have greater unsteadiness. Methods of controlling the flow around the wheel may provide a means of reducing

Conclusion
The importance of modelling the full unsteadiness of the flow in vehicle soiling predictions has been assessed through the comparison of results for two configurations of a generic SUV model. These were obtained in fully unsteady simulations and using just the mean flow and a simple statistical approach to model velocity perturbations seen by the dispersed phase. The two chosen configurations produce very different wake structures, giving unique base contamination patterns and therefore representing two good test cases for this study. The wake of the first configuration is dominated by a strong upwash associated with the diffuser, while the square aft body of the second configuration produces a wake that is seen to be sensitive to even very small asymmetries in the setup. The flow-field was computed using the Improved Delayed DES method, together with the Spalart-Allmaras RANS model. The dispersed phase was tracked using the Lagrangian particle tracking model. Experimental data were used to validate CFD results.
It has been shown that the overall qualitative agreement between the computational results using the fully unsteady method and experimental results is good. This statement is particularly true for the SUV configuration with a 30 diffuser installed. For this case, the high fidelity eddy resolving CFD method was able to accurately predict both the flow structure and the pattern of base contamination. In addition, the area of highest soiling intensity was also matched. The prediction of soiling on the base of the SUV configuration with no diffuser installed is less accurate, however, the unsteady flow field is also less well captured. Although the trend from one configuration to the other is predicted, the predicted pattern is located higher on the base than in the experiments and is also not as widely spread across the surface. This is consistent with differences seen in the predicted wake structure compared to PIV and base pressure distribution. This shows the importance of accurate unsteady flow field prediction as prerequisite to accurate soiling prediction.
The soiling pattern with the stochastic modelling approach shows some agreement with the fully unsteady results. The general part of the base contaminated is the same for both methods, but details of the pattern and the soiling rate do not match. It has been shown that the mechanism of particle entrainment into the wake is very unsteady, revealing large scale, spatially correlated, flow structures in the process. The whole process incorporates high frequency events that are superimposed on lower frequency events. The energetic low frequency events also seem to be responsible for transporting large particles into the wake. This suggests that an accurate representation of flow unsteadiness, especially in the injection region, is essential to accurately model particle advection into the recirculation zone and thus obtain the correct rate of liquid deposition. In order to achieve this, full unsteadiness of the flow must be considered. Control of these flow structures is also a potential method of reducing base contamination. Once in the wake, the dynamics of particles can be adequately described by the mean flow-field and a statistical description of turbulence. This suggests that the statistical model may still be used to get a crude estimate of the most contaminated regions on the surface, considering their small computational price, particularly if flow simulations have already been carried out for other purposes. For example, such an approach may be used to get a crude idea of the places on a vehicle's exterior where sensors or cameras could be installed.
A possible reason for the discrepancies between experiment and fully unsteady results is that a limited number of sub-models were used. For example, a large number of small particles are seen to be predicted to impinge on the very edge of the base in a position that is not seen in experiment. It is possible that additional processes such as particle coalescence and evaporation, which would remove small particles, may be important and are proposed as topics for further study. In addition, work is planned to study the dynamics of spray released from the contact patch of a rotating wheel, which should give a more realistic picture of the phenomenon.