Indoor Air Quality (IAQ) issues have been receiving much research attention from various disciplines over the last couple of decades [1] [2] [3] [4] [5] [6] . This has been mainly triggered by the need of promoting more comfortable and healthy indoor environment on one hand, and by the necessity of protecting indoor environment against the intentional release of biological and/or chemical agents on the other hand.
After the epidemic outburst of severe acute respiratory syndrome (SARS) and avian influenza in East and Southeast Asia, there has been a growing research interest in studying the transport and control of airborne bacteria and viruses indoors and in confined environments such as aircraft cabin [7, 8] . Dispersion of microorganism-laden aerosols exhaled from infected patients was recognised as a potential airborne transmission pathway. Thus, proper understanding of aerosol transport is required to improve exposure assessment tools and models and adopt better ventilation strategies that can substantially reduce indoor particle concentrations and improve the indoor air quality.
Depending on the original and final size, droplet nucleus can remain suspended in air for several hours and thus distribute widely throughout indoors. This also depends on the ventilation scheme used. Indeed, the ventilation system determines the airflow pattern in the room which in turn decides the droplet nucleus fate. Displacement ventilation has been acknowledged to provide better indoor air quality than mixing ventilation. Some existing studies that used passive gaseous and small inertial particles as contaminant sources stated that in the presence of heat sources, displacement ventilation is more contributive to pollutants removal without much mixing to the whole indoor environments and this compared to other types of ventilation [9, 10] . These studies also showed that in order to design an effective ventilation system, it is crucial to have a reliable tool that is capable to predict airflow pattern and particle distribution and dispersion indoors. This can be achieved through the use of Computational Fluid Dynamics (CFD) which has the capacity to provide microscopic information on the indoor air environment like the air velocity, pressure, temperature, and pollutant's concentration distribution which are useful to obtain pertinent macroscopic parameters for engineering goals.
Among these CFD tools, Reynolds averaged Navier-Stokes (RANS) simulation has been extensively used for simulation of airflow indoors to predict the averaged velocity or temperature. Large Eddy Simulation (LES) solves directly for the transient behavior of the large-scale turbulent motion which tends to have the greatest influence on the turbulent transport but less adjustable constants are required [11] . LES approach to model turbulent flows has been recognised as a powerful tool that is able to satisfy a continuing desire for higher fidelity of predictive capabilities. Its use has been invigorated by the extraordinary development of computational facilities. Since LES solves time-dependant turbulent flows, it can provide detailed description of the turbulence phenomenon, such as three-dimensional instantaneous velocity field. Thus, LES properly accounts for the history and transport effects of turbulence on aerosol dispersion and deposition. The very informative picture of the turbulent transport provided by LES has usually been associated with a large computational expense that has prevented its use in the past. However, many recent works have shown that LES computations carried out on current computer workstations have been successfully applied to several indoor airflow problems with reasonable computation overheads [3] . Among these studies, we can mention the work of Emmerich and McGrattan [12] and Musser and McGrattan [13] who applied LES to investigate isothermal room airflow. Jiang and Chen [14] who explored the potential of LES to study natural ventilation in buildings. Zhang et al. [15] evaluated various turbulence models, including LES, in predicting airflow and turbulence in enclosed environments and Zhang and Chen [11] proposed a filtered dynamic subgrid scale model to be used in conjunction with LES of indoor airflow. All these investigations have shown the great potential of LES in predicting turbulent airflow in enclosed environments but none of them have considered thermal effects generated by human.
For LES of two-phase turbulent flows, the numerical simulation of the particulate trajectories and dispersion pattern in airflows are treated along two streams namely the Lagrangian and the Eulerian methods. The Eulerian methods consider the particle phase as another continuum. Transport equation for the mass concentration is derived from the mass (species) conservation conditions and solved to give details of the particle concentration field. When the particle phase is considered as scalar species, the gravitational settling and the deposition rate should be taken properly into account to reflect the aerosol's inertial character [1, 6] . For the Lagrangian methodology, the motion of the discrete matters are described by the force balance including those induced by the interactions with the carrier phase. In indoor particle studies, each method gains its own reputation depending on the research goals. The Eulerian method is widely used to predict particle concentration distributions in rooms. Generally these simulations agree well with experimental data although remarkable discrepancies exist in some studies [4] . If the particle motion and dispersion history is of interest, the Lagrangian method is more suitable. Indeed, the Lagrangian approach appears to provide a more natural description of the actual physical phenomena compared to the Eulerian approach. Investigations using two-phase LES to study particle transport in enclosed environments are rare in literature. We can mention the works of Tian et al. [16] , Bouilly et al. [17] and Chang et al. [18] .
Almost all of the previous studies did not consider thermal effects. For conventional ceiling schemes the influence of high momentum airflow overwhelms the buoyancy effect and the pollutant transport is mainly governed by the bulk airflow [5] . Nevertheless the working principle of displacement ventilations relies completely on the existence of internal heat sources, i.e. occupancy. Hence it is very important to study aerosol transport by displacement ventilation with internal heat sources. In the present work, measurements and numerical simulations of contaminant particle concentration evolution in a mechanically ventilated room with heated sources have been carried out. This investigation aims to provide a realistic simulation of the time-dependant flow field in a chamber populated with two heated manikins using the wellresolved LES approach. The potential of both the Lagrangian and the Eulerian approaches coupled to LES of non-isothermal airflow to study dispersion characteristics of expiratory aerosols has been explored.
It is important to mention that by conducting this study, there was no intention to treat the experimental and numerical findings as portable results and hence draw practical conclusions about the effectiveness of the displacement ventilation scheme for the contaminant removal in full-size model. This is because the scaling rules were not respected, in particular for the thermal aspects due to the physical restrictions imposed by the type of chamber used for the experiment. However, the aim of this work is twofold: (i) to test the capabilities of the in-house code Code_Saturne to investigate aerosol transport under displacement ventilation with two heated human models using LES, (ii) to compare the aerosol concentration fields computed by both the Lagrangian and the drift-flux Eulerian approaches in the framework of LES, and (iii) to compare the experimental results with model predictions.
Both measurements and numerical simulations of airflow field and particle concentration were carried out in a downscaled chamber with two identical model occupants as depicted in Fig. 1 . The chamber is mechanically ventilated using the displacement ventilation scheme. Table 1 shows the details of the room and occupants geometries and the boundary conditions. The center point of each object is indicated in the table. The geometry of the human occupant used in this study is the one originally proposed by Brohus and Nielsen [9] . An opening (0.004 m Â 0.004 m), measured at 0.33 m above floor, located at the centerline of the head was added to simulate the mouth of the occupant. Table 2 shows the details of the geometry of the model occupant used for this study. Two planes are defined in the geometry; a plane (X-Y) crossing through the human model near at the inlet and a midplane (Y-Z) at X ¼ 0.25 m. The wall temperature is set to 297 K while the occupant temperature is set to 317 K. One of the occupant is emitting droplets (source) and faces directly the second occupant (receptor). The source emits water spherical droplets with initial velocity of 10 m/s lasting for 0.1 s. Similar duration period has been adopted in the literature [19] .
In this work, some well-justified assumptions were made. First, the evaporation period of the emitted droplets was ignored and only the droplet nuclei is directly modelled. Two recent reviews have demonstrated the time scale of evaporation. Nicaset et al. [20] and Morawska [21] estimated that the shrinkage time from the original droplet to droplet nuclei is rapid and is in the order of 0.5 s. This time scale is at least an order of magnitude shorter than the residence time of the droplet nuclei suspended in the room. Also, the droplets are assumed trapped once they touch any surfaces and do not re-suspend or break up. These assumptions are valid for the present low air velocity environment. Coagulation effect has been examined by applying a simple estimation [22] . The results revealed that the coagulation effect can be neglected.
The filtered spatial and temporal evolution of non-isothermal incompressible Newtonian fluid flow is governed by the following equations:
where u i is the component of the filtered fluid velocity in the x i direction, p i is the fluid pressure, T is the filtered temperature, r is the fluid density, n is the fluid kinematic viscosity, g j is the gravitational acceleration, b is the thermal expansion coefficient, Pr is the molecular Prandtl number and d ij is the Kronecker symbol. s ij and Q j are the sub-grid scale (SGS) stress tensor and heat flux, respectively. The sub-grid scale (SGS) stress tensor s ij is modelled using the algebraic eddy-viscosity model proposed by Smagorinsky [23] : Table 1 Geometry and boundary conditions. 
where n SGS is the sub-grid scale viscosity:
Here C s is the Smagorinsky constant, its value is taken equal to 0.12.
of-strain tensor. The length scale D is equal to the grid spacing,
is the Van Driest damping function that accounts for the reduction of the sub-grid length near solid walls [24] . It is based on the dimensionless distance from the wall y þ ¼ yu s /n. u s is the friction velocity and A þ is taken equal to 25.
In the same way the heat flux Q j is modelled:
The subgrid scale Prandtl number Pr SGS is taken equal to 0.5 [25] . An unstructured grid (non-conforming embedded refinement) consisting of 1,360,000 cells was used for the computational domain discretisation. The first grid point near the chamber walls and manikins at which the velocity is computed is located at y þ z 1. Two grid points are placed within the viscous sublayer, the depth of which equals 5 wall units. A non-uniform grid is employed in the normal-to-the-wall and normal-to-manikins directions. This is done in order to locate more gridpoints where they are most needed.
Care should be taken when non-conforming embedded refinement is used. It is well known that the near-solid region is characterised by steep gradients and very small energy containing eddies that should be well captured. These near-solid coherent structures contain most of the turbulence and are responsible for the correct distribution of the turbulent energy from the streamwise into the other directions. Moreover, the near-solid turbulence has a significant impact on the deposition of inertial particles and therefore it should be properly resolved. The different non-conforming layers were designed such that the aspect ratios between the different directions for all the cells do not exceed 5. This is done to avoid numerical instabilities that arise when using very flat cells in the LES context. Thus the use of non-conforming layers should be efficient since it allows a cell distribution that responds to the requirement of the flow dynamics without introducing further numerical errors.
Speziale et al. [26] pointed out that a reliable LES is the one that becomes a DNS when the grid resolution is as small as the Kolmogorov scales. In this numerical study, the Reynolds number is 1300 based on the length of the inlet opening and the expected Kolmogorov microscale is 70 microns. Consequently, one cannot seek a grid independent LES, as we usually do for RANS. This is because a grid independent LES is essentially DNS, and the philosophy of LES that is based on grid dependency, looses its meaning. Celik et al. [27] developed a method to assess the quality of LES results. It consists of estimating an index of quality which is a measure of the percentage of the resolved turbulent kinetic energy. They stated that if more than 75% of the kinetic energy is captured, then the LES is considered adequate. Pope [28] has also stated that the amount of turbulent kinetic energy that is carried by the SGS scales should not exceed 25% for the LES to be well resolved.
A simple way of estimating the SGS kinetic energy can be made based on the assumption of equilibrium at the cutoff. Then the dissipation rate and the SGS kinetic energy can be evaluated as the following:
Typically C e z 1 [29] . For non-equilibrium turbulent flows characterised for instance by zones of strong recirculation and/or boundary layers detachment, the assumption of equilibrium at the cutoff is not valid and it is necessary to solve a transport equation for the SGS kinetic energy to get an accurate estimation of the filtered out kinetic energy [30] . In this work, Eqs. (9) and (10) were used to give such an estimation. A time step, Dt ¼ 0.01t * was used to advance calculations. t * is the integral time scale defined as the ratio of the inlet height to the inlet velocity. The size of the time step was dictated by the numerical stability. The LES computations are initiated from a randomly generated instantaneous inlet velocity with mean velocity and turbulent kinetic energy profiles fitted to analytical formulae [31] . The time advancement was carried out until t ¼ 150t * to achieve a flow field independent of the initial conditions. At t ¼ 150t * , residuals of Eqs. (1)-(3) became smaller than the set convergence tests indicating that the computations had reached a nearly statistically steady state. From t ¼ 150t * , the calculations were continued until t ¼ 200t * . In this interval, the final statistical data was accumulated.
A flow solver from the R&D section of Electricite de France named Code_Saturne (http://www.code-saturne.org) was used as starting point of the present work. The discretisation in Code_Saturne is based on the collocated finite-volume approach. It allows solving Navier-Stokes and scalar equations on hybrid and nonconform unstructured grids. Velocity and pressure coupling is ensured by the prediction/correction method with a SIMPLEC algorithm. The collocated discretisation requires the Rhie and Chow interpolation in the correction step to avoid oscillatory solutions. A second order centred scheme (in space and time) is used. The flow solver has been extensively tested for LES of single-phase flows [32] .
Aerosols are released and tracked in the turbulent flow that is described in the previous section. The physical properties of these inertial particles are summarized in Table 3 .
As a result of the high density ratio between particle and fluid densities, the equation describing particle motion is reasonably simple and only the drag and gravity forces will be retained since other forces are in this case negligible [33] . Since the dispersion of very small particles is investigated in this work, the Brownian force should be taken into account. Thus, the tracking of the inertial particles within the turbulent flow obeys the following system of equations: dx p;i ¼ u p;i dt; du p;i ¼ u p;i À u s;i s p dt þ n i ðtÞdt À g i dt; Here x p and u p are the particle position and velocity, u s is the fluid velocity seen by an inertial particle along its trajectory, g is the gravity force per unit of mass, d p and r p are the diameter and the material density of inertial particles, s p is the particle response time, C D is the drag coefficient and Re p is the particle Reynolds number, Re p ¼ d p ju s À u p j=v with n is the kinematic fluid viscosity. C n is Cunningham slip correction factor. It is considered herein to correct the drag coefficient in order to take into account the free-slip boundary conditions that occur at the surface of the particles.
Here l is the molecular mean free path. n i (t) is the Brownian force per unit of mass.
Here G i is zero-mean, unit variance-independent Gaussian random number, Dt is the time step and S 0 is spectral intensity which is computed using:
The number of particles injected was set to 500 per time step for a period of time lasting 0.1 s which corresponds to a number of iteration equal to 500 iterations. Both particle and fluid seen velocities were set equal to the fluid velocity at the secondary inlet (occupant's mouth); i.e. 10 m/s. Calculations were repeated using 1000 particles and only slight differences in the concentration fields were noticed (around 2%). Because of the high velocity at the secondary inlet, the time step of the simulation was decreased to Dt ¼ 0.001t * to keep the Courant-Friedrichs-Lewy number around 1.
For the Eulerian approach, the concept of a particulate phase consisting of individual, distinguishable droplets is abandoned. The particulate matter phase is considered as continuum that can be described using a set of generalized equations similar to the equations used to solve the gas phase. This approach is known in the literature as the two-fluid approach. If only the instantaneous mass concentration of the particulate phase C is of interest as it is the case of this study, a transport equation by turbulent motion for this property can be derived:
where G is the particulate matter mass diffusivity and S c is the rate of creation or destruction of the mass concentration per unit volume. In the context of LES, Eq. (16) should be spatially filtered giving rise to the filtered transport equation for particulate phase mass concentration:
The term vF j =vx j is considered as a supplementary turbulent diffusional process occurring at the sub-grid scales. This additional concentration flux is approximated following:
where Sc SGS is taken equal to 1. Eq. (18) The sink term S c is computed as the mass wall flux per unit of volume. It is computed for all the cells that have at least one wall face. This sink terms account for a decrease in mass concentration due to loss of mass to the walls. This sink term corresponds to particle deposition in the Lagrangian approach.
where A is the area of the wall face linked to its corresponding cell. v d is the deposition velocity estimated for the different wall orientation using the three-layer model developed by Lai and Nazaroff [34] . This was done by integrating the particle loss across the boundary layer. This afformentioned approach is not a new approach and its use in the framework of RANS has been discussed by Holmberg and Li [1] and Chen et al. [6] . residual charge. The airflow was induced by means of a DC fan and regulated by a power supply. The inlet duct length was sized at least 60 times of hydraulic diameter of the size of the duct to ensure that the flow is fully developed at the chamber inlet.
Monodisperse particles of 10 mm were generated by atomization of diluted standard polystyrene microsphere suspensions (Thermo Fisher Scientific). The diluted suspension was pre-filled to the cup container attached to a spray gun. The spray gun was fixed at the centreline (X ¼ 0.25 m, Y ¼ 0.33 m, Z ¼ 0.4 m) of the head of the source. A flow regulating valve was used to adjust the emission velocity to match with the simulation boundary conditions. The expiratory process was mimicked by a short release of particles through a spray gun connected to a compressor. The spraying duration was controlled by a timer circuit and it can be adjusted by a LabVIEW program. In this work 0.1 s was selected.
HEPA filters were installed at inlet and outlet to minimize the background particle count inside the chamber and to prevent crosscontamination. Type T thermocouples were selected as it has better accuracy than the other types of thermocouples. Prior to the temperature measurement, the thermocouples were all calibrated in-situ by 5-points measurement.
The two manikins that are made of aluminuim were heated by wrapping heating wire around the body. Thirteen thermocouples were used to measure inlet, manikins' surface temperature and different pre-fixed locations inside the chamber. The air inlet temperature is the room temperature and during the entire experimental period the inlet air temperature was maintained at approximate 24 C. The temperature difference between the inlet air and the surface temperature of the manikins was kept at 20 C. All the temperatures were monitored and controlled through the same LabVIEW program.
The emission velocity was measured by a thermal anemometer (TSI, 9555). The particle concentration was measured by an optical particle counter (TSI, 3775). Conductive sampling tube was used to sample the particles to minimize the electrostatic loss. Since there was one counter available, only a single point measurement made at a time. Background concentration was measured 5 min prior to the emission starting. The measured concentration was subtracted from the background.
Figs. 2 and 3 show the predicted steady state averaged velocity contour at the mid-plane of section Y-Z and at the plane crossing the manikin closer to the air inlet respectively. A cooled-jet can be easily identified near the low-level of the floor and a strong vertical buoyancy-plume is formed above the heated manikin. One key feature of displacement ventilations is the low inlet velocity. As the cold air flows around the manikin and picks up the heat, the velocity increases rapidly. It is noticed that the low momentum cooled air (0.2 m/s) near the floor level absorbs heat from the two occupants and creates a dominant vertical thermal plume in the boundary layer around each occupant. The airflow velocity is fairly weak in all regions except inside the buoyancy plume [35] and the bulk velocity is less than 0.1 m/s. This low velocity will influence significantly the convective transport of aerosol. Fig. 4 shows the predicted averaged temperature contour along the mid-plane of Y-Z direction. This shows you a clearer picture on how the temperature varied from the cooled-inlet to the exhaust. Fig. 5 shows an instantaneous velocity at the same plane. Two similar but unequal strengths of vertical plumes can be observed. Inferring carefully from both figures, it can be seem that the buoyancy plume of the heated manikin near the inlet is warmer and stronger than another one. It is due to the coolest air contacting that heated manikin while the air contacting the other manikin has been warmed.
It is not straightforward to compare the experimental data with those modeling results directly. The particle counter output is number (or mass) concentration expressed in particle number (mass) per cm 3 . The native output for the Eulerian and Lagrangian models are concentration or dimensionless concentration. To compare with both Eulerian and Lagrangian modeling results, the (concentration) data obtained by the counter was transformed to a dimensionless concentration and is defined as
where C(x, t) is the temporal concentration at the measuring point x and C Vol ðt ¼ 0Þ is the initial volume-averaged concentration for the entire chamber. It was evaluated by calculating the particle injected during the emission period (0.1 s) and divided by the volume of the chamber. Figs. 6-8 show the experimental measured steady state airflow pattern and temperature distribution prior to the emission of aerosols. The results were compared with those predicted by the LES. All the results were measured at mid-plane along Y-Z. The velocity profile along the centerline of the inlet duct (i.e. Y ¼ 0.075 m) is shown in Fig. 6 . Excellent agreement between the modeling and experiment can be observed. Due to the very low velocity and the sensitivity of the velocity probe used, there was no measurement data for Z greater than 0.6 m. It is interesting to note that a small recirculation zone is existed at Z ¼ 0.7-0.75 m. Fig. 7 shows the temperature measurement along a central vertical line (i.e. Y from 0 to 0.5 m). The computational prediction experimental data consistently overestimates the measured temperature for most of the locations. It is interesting to notice that the discrepancy decreases with the Y-direction. The simulated results were fluctuated temperature while the experimental results were averaged value. Besides, the discrepancy may also be attributed to the uncertainty of the thermocouples. Fig. 8 depicts the temperature profile along the Z-direction at the vertical height of the spray gun. Since the separation of the two manikins was less than 0.2 m, seven thermocouples were placed. It can be seen that the measurements and modeling results agree well as it shows the increasing trend when it approaches the heated manikin.
Concentration evolutions were measured in three streamwise locations for two different transverse points. Figs. 9-12 show the experimental results for four locations and the results are compared with the two modeling predictions. The x-axis is the elapsed time after the aerosol emission and the y-axis is the dimensionless concentration. It should be noted that unlike Reynolds averaging methods LES gives instantaneous concentrations of particles that are convected by a turbulent flow. All the reported results were measured downstream of the emission point. Hence it can be anticipated that the concentration will decay with elapsed time. Regardless of the modeling approaches or experimental results, first three figures show clear decay profiles. For Fig. 12 where the location is at the exhaust, the result is more complex as both Lagrangian and experimental results depict non-decay profiles. It can be seen that the variation of the experimental measurement at point F is very small which is ranging from 0.01% to 0.07%. Further investigation is required to explain the observation. It is also interested to compare the ''decay'' rate for different locations as the decay magnitude should be correlated to the ''local'' air velocity. This issue will be addressed in future work. The experimental results were best fitted by an exponential function using the least square method and the corresponding decay rates are shown (Figs. 9 and 10). Inferring from the results Figs. 9 and 10, the decay rates are comparable to each. By taking into account of the uncertainly of the particle counter used which is AE12%, the agreement between the experimental results and the computational modeling predictions is acceptable [36] . In this work observable discrepancy between the modeling results by Lagrangian and Eulerian approaches is seen. More experimental results are required to validate the accuracy of the model.
In this work particle dispersion in a scaled chamber was studied numerically and experimentally. A well-resolved LES was performed to study aerosol dispersion in a turbulent scaled chamber populated by two manikins. Numerical predictions were compared to the experimental observations. Due to the complexity of the turbulent flow in such an enclosed environment, a good deal of care has been taken to ensure that the simulation of the carrier phase is reasonably accurate. Efforts were made to adapt the mesh to the dynamical features of the flow. The numerical predictions showed a good quantitative agreement with the experimental findings for both phases though some discrepancies were noticed in some regions of the chamber.
Both steady state and transient parameters were measured experimentally. The measured steady state temperature and velocity agreed well with the simulation. Particle concentration was measured a CPC. The measured concentration was lower than that predicted by Lagrangian and Eulerian models. By taking into account of the measurement uncertainty of the instrument and some limitations of the set-up, the result comparison is acceptable. 