In this paper, particle dispersion and spatial distribution in a full scale (5.5 m x 2.4 m x 3.7 m) forced ventilated room are investigated using four different multiphase flow models, including passive scalar model, discrete particle phase model, mixture model and Eulerian model. The main differences between these four models lie in how the particles are modeled. A two layer k-ε turbulence model is used to calculate airflows. Simulated airflow characteristics and particle concentration are compared with corresponding experimental data. The results show that only discrete particle phase model could predict particle concentration distribution close to experimental values and satisfy the published validation criteria (ASTM D5157-97). The reasons for the failure and success of these models in the present case are discussed. Furthermore, the effects of turbulence models of airflows and treatment of boundary conditions on the particle concentration are also investigated.

## Introduction

Airborne particulate spatial distribution (APSD) has a direct and profound effect on indoor environment quality and occupational health. Many experimental studies have been conducted to investigate the particle behaviors and concentration in the past years in either prototype or scale model rooms. While valuable information has been provided, measurements are always costly and time-consuming. In recent years, CFD (computational fluid dynamics) becomes the routine tool in the field of indoor air quality due to its ability to quickly provide detailed information on airflow, particle concentration, deposition and movement in various ventilated spaces with relatively low cost.

The approaches for the numerical simulations of particle matters in the airflows can be divided into two categories: Euler method and Lagrangian method.^{1,2} Passive scalar model (PSM), mixture model (MIX) and Eulerian model (EUL) belong to the first approach while discrete phase model (DPM) falls into the second method.

In earlier times, PSM model was adopted as the main tool to simulate particle transport.^{3} In this model, particulate matters are assumed as passive scalars, that is, particles move in the same manner as airflow and there is no relative movement between particle and airflow. The limitation of this model is obvious: only small particles could be assumed as passive scalar and thus their movement could be determined by airflow. In fact, many factors could make PM behave differently from airflow, such as coagulation, sedimentation and deposition. Among those factors, sedimentation caused by particle gravity plays an important role in determining the properties of particles in indoor airflow. Murakami et al^{4} proposed the modified passive scalar model (MPSM), which takes into account the settling velocities (due to gravity of PM) of particulate matters. Compared to the PSM, MPSM gives much better predictions in some convection dominant indoor flow fields.^{4,5} However, the underlying assumption of MPSM that particles follow the airflows instantly is not always valid in indoor airflow. Furthermore, the assumption that gravity is imposed on particle instantly is another limitation of the MPSM.

Both of the above models neglect the inertial effects of particles. Discrete phase model (DPM) is a good candidate for overcoming the above restriction. In this model, all possible forces on the particulate matters are considered and the particle concentration could be derived from the individual particle trajectories. In recent years, the DPM model has been popular in indoor airflow field.^{6–789} The major limitation of DPM is that it requires much more computation time than PSM and MPSM. Furthermore, using computational particles to represent a packet of real particles implies that DPM model could only be applicable to very diluted particle flows.^{10}

MIXand EUL models are two widely used multiphase models for airflows with suspended particles. In these two models, both phases (air and particulate matters) are treated as interpenetrating continua and interaction between the two phases are taken into account. In the MIX model, momentum equations for the mixture of air and particles are solved and relative velocity between particles and air is prescribed on the assumption that a local equilibrium between the phases is reached over a short spatial length scale. In EUL model, momentum and continuity equations for both phases are solved simultaneously. The coupling between air and particles is taken into account through pressure and interphase exchange coefficients. Both MIX and EUL models have been widely used in the two-phase flow motion of particles.^{1,2,11}

It is well known that no single turbulence model is universally accepted as being superior for all turbulent flows. The same case seems to hold for multiphase flows: Multiphase models should be selected properly according to the simulation purpose and fluid property of the problem. Thus it is necessary and useful to know the performence of each multiphase model for indoor air quality simulation. This investigation will conduct a detailed analysis of these four models’ performance in predicting the partricle dispersion and spatial concentration in a full scale ventialted room by comparising with corresponding experimental data.^{12}

There have been many papers devoted to investigate the performance comparisons of mathematical models,^{13} but criteria for validating different models adopted by different authors are various. It is accepted that predicting trends may be more important than achieving high accuracy of measurement.^{14} However basic criteria should be established to make the comparisons meaningful and useful. The ASTM guide^{15} provides such criteria for assessing performance of mathematical models in predicting particle concentrations. In this paper, all the predictions based on the four multiphase flow models are compared to experimental results according to the ASTM guide.

## Experiments

The dimension of ventilation room is 5.5 m x 2.4 m x 3.7 m (length x height x width). The room has a continuous slot inlet along one side wall and a slot outlet on the opposite wall. The inlet height was 0.05 m and the outlet height is 0.2 m. In this room, the inlet is 2.05 m from floor on the left side and outlet is 1.3 m from the floor on the right side. The midplane of the room is shown in Figure 1. For more details about the room, see Zhang et al.^{16}

A multi-sampler, which used an array of critical venturi orifices for controlling the airflow rate to each sampling point, was applied to measure the particle matter spatial distribution.^{12,17} Particle concentrations at 25 points within the testing room were measured and measurements points were uniformly distributed in the central cross-section of the test room (Fig. 2) based on the assumption that room airflows was two dimensional.^{18,19} Particles were uniformly emitted downwards using a dust emission system with 25 emitting ports evenly distributed over the entire floor area.^{12} The initial velocity of particles was estimated in the range of (-0.1, -1) m/s, where “-” indicated that particles were emitted towards floors (negative y direction). The emission ports with 1.6 mm diameter were located on the bottom of tubes that were 12 mm above the floor. A rotating-table dust generator with precisely controlled rotating speed was used to feed the particle to the particle emission system.

The measured particle size distribution is shown in Figure 3. The mean aerodynamic diameter of dust is 1.63 µm with density of 2.65 g/cm^{3}. The maximal uncertain is estimated to be 2.3%. The ventilation rates, particle generation rates and related quantities are listed in Table 1. The experimental data^{20} are used in this paper to validate the four multiphase models.

## Mathematical Models and Numerical Methods

### Turbulence model

In the present study, the air in the room is regarded as continuous and incompressible fluid since the air velocity is much less than the speed of sound and the density of the air does not change with pressure. Reynolds Averaged Navier-Stokes (RANS) equations with turbulence models are used to predict turbulent airflows in the room. Appropriate treatment of wall boundary conditions is necessary for a successful implementation of turbulent models. It is well known that standard k-∊ turbulence model is only applicable to high Reynolds number flows and fails to predict satisfactorily the constants in the wall law.^{21} While wall functions are commonly used to overcome this difficulty, it should be mentioned that they are not always applicable in indoor airflow predictions due to the existence of separation and reattachment phenomena in indoor airflows. To overcome this difficulty, a low-Reynolds number k-∊ turbulence model^{22} that could be integrated down to the walls is used in the present study.

In the two layers near wall k-∊ model,^{22} the Reynolds stress tensor $\overline{-{u}_{i}{u}_{j}}$ is closured based on eddy-viscosity assumption:

_{ij}is the Kronecker delta,

*k*is turbulent kinetic energy,

*v*

_{t}is turbulent eddy viscosity,

*U*

_{i}(or

*U*

_{j}) is mean velocity compent, and

*x*

_{i}(or

*x*

_{j}) is spatial coordinate component. In the fully developed turbulent region where turbulent Reynolds number ${\mathrm{Re}}_{y}=y\sqrt{k}/v\ge 200$ (

*v*is molecular viscosity and

*y*is the normal distance from the wall at the cell centers), the standard k-∊ model model

^{23}is applied while in viscosity-affected region, the one-equation turbulence model

^{24}is used. Thus turbulent eddy viscosity would adopt different expressions in the two regions: where ∊ is turbulent dissipation rate and the length scale

*l*

_{µ}is specified by

Similar to turbulent eddy visocity *v*_{t}, turbulent dissipation rate is also treated separately in the two regions: in the fully developed region it is obtained by solving the transport equation while in the viscosity-effected region it is calculated algebraically by the expression:

The transport equations for turbulent kinetic energy and dissipation rates are written as

## Table 1

Summary of experimental parameters.

where ρ_{k}, and ρ

_{∊}are the turbulent Prandtl numbers for turbulent kinetic energy

*k*and turbulent dissipation rate ∊, respectively.

*S*is the modulus of the mean rate-of-strain tensor, defined as

The model constants are listed in Table 2.

### Boundary conditions

Three kinds of boundaries are used: wall boundary, inlet boundary and outlet boundary. All the walls are assumed to be smooth and non-slip. Uniform mean velocity across the inlet was widely used^{25–2627} however, in some cases this method gave poor agreement with measurements while detailed profiled velocities got from measurements led to better predictions.^{28,29} Thus in this paper mean velocity and turbulent intensity (*I*) at the inlet are taken from previous hotwire measurements at the same Reynolds number.^{16} The hydrodynamic diameter for inlet is estimated about 0.025 m. Thus turbulence length scale *l* is approximated as 0.00175 m. Turbulent kinetic energy is determined from the following relation:

## Table 2

Constant.

While turbulent dissipation is determined by### Multiphase Flow model

#### Passive scalar model (PSM)

In the passive scalar model, the particles are assumed to behave like gases and have no effects on ariflows. The space-time evolution for the concentration *C* of the passive scalar is govern end by the following advection-diffusion equation:

*D*is the molecular diffusivity of the scalar, and

*S*is source term.

Zero-gradient boundary condition of concentration is applied at the outlet. At walls and inlet, zero particle concentration is used.

#### Discrete particle model (DPM)

As discussed in the previous section, DPM calculates the individual particles trajectories by considering the effects of all forces on particles in Lagrangian frame. The governing equation for each particle is as follows:

where*u*

_{p},

*u*are particle and airflow velocities respectively,

*F*

_{D}(

*u*-

*u*

_{p}) is drag force per unit particle mass,

*F*is any additional acceleration term, ρ

_{p}is particle density, ρ is fluid density, and

*g*ρ

_{p}- ρ/ρ

_{p}is gravity force.

The drag coefficient is calculated as:

where µ is fluid viscosity,*d*

_{p}is particle diameter and

*C*

_{c}(Cunningham correction factor) is calculated using where » is the molecular mean free path.

*F* represents the additional forces exerted on particle, including pressure gradient force, Basset force, virtual mass force, Brownian force and Saffman lift force etc. Considering the particle size distribution in the experiments (Fig. 2), Brownian force and Saffman lift force were included in the calculation.

The amplitudes of Brownian force take the following form:

where ζ is normally distributed random numbers, Δ*t*is time step and

*S*

_{0}is the spectral intensity.

^{30}

The Saffman lift force is calculated as:^{30}

*K*= 2.594 and

*d*

_{mn}is the deformation tensor.

Since the grids adopted in this study are sufficiently fine near the wall regions, trap-type boundary conditions are set for all the walls except ceiling, that is, the walls will trap particles when they reach the walls. For the ceiling, the resititution coefficients are set to 0.7, that means, particles reaching the ceiling can escape from the ceiling once enough momentum is gained. When particle reach the inlet and outlet, they will escape from the room and the trajectories terminate.

### Mixture Model(Mix)

In the mixture model, the air and particles are treated as interpenetrating continua. Momentum equations of air and particle mixture are solved and relative velocity between air and particles are considered.

The governing equations are the continuity equation for the mixture, the momentum equation for the mixture, the volume fraction equation for the particles, and the algebraic expression for the relative velocity.

where ρ^{m}is mixture density, µ

_{m}is mixture viscosity,

*u*

_{m}is mass-averaged velocity and

*u*

_{Dk, i}is drift velocity, which is related to relative velocity

*v*

_{Dk}according to the following relation:

The algebraic relationship for the relative velocity is obtained under the assumption that local equilibrium between particle and airflow is attained in a very short spatial length scale:

where α is particle acceleration and given by gravity in the present study.Volume fraction equation for the particles is given by

The following assumptions are adopted in the current study: (1) mass transfer between air and particles are neglected; (2) there are no particle coagulation.

### Eulerian Model(Eul)

In the Eulerian model, the momentum and continuity equations for air and particles are solved simultaneously and coupling between gas and solid is considered through the pressure and interphase exchange coefficients.

The continuity equations for both airflow and particles could be taken in the following form:

where α*f*, ρ

*f*and

*u*

_{fj}are fluid phase (air) volume fraction, density and velocity respectively. Correspondingly, α

_{s}, ρ

_{s}and

*u*

_{s, i}. are solid phase (particle) volume fraction, density and velocity respectively.

The conservation of momentum for the fluid phase (air) is

and for the solid phase (particles) where*p*is the pressure shared by both air and particle phases, τ is the phase stress tensor,

*F*is the extrenal body force,

*F*

_{lift}is the lift force,

^{31}

*F*

_{vm}is the virtual mass force,

^{31}

*K*

_{fs}=

*K*

_{sf}is the interphase momentrum exhange coefficients

^{32}and subsrcipt

*f*and

*s*denote the fluid (air) and solid (particle) phases, respectively.

The phase stress tensors for fluid and solid phases take the following forms

where µ_{f}and µ

_{s}are viscosity for fluid (air) and solid (particle) respectively.

The following boundary conditions are applied for Mixture and Eulerian model: zero volume fraction at the inlet and zero-grdient of volume fraction of particles are set at the inlet and outlet respectively.

The intial emission velocities of particles are needed for the simulations. As exact values are not available from the correspoinding experiments, nine velocities in the range (-0.1, -1) m/s with the increase of 0.1 m/s are tested. It is found that the -0.5 m/s generated closest particle concentration with experimental vlaues for DPM simulation. For the other three models, it is hard to find a preferred value which performed better in the predicition than the rest values when compared with corresponding experimental data. In fact, the differences in the predictions with values of initial velocities in the range (-0.4, -0.6) were relatively small. Thus all the results shown later are based on a value of -0.5 m/s initial velocity. Finally particle size distribution (Fig. 2) is implemented through user-defined functions (UDF).

### Numerical methods

The numerical solutions of the governing equations for airflow and particles were carried out with a commercial finite-volume based program, ie, CFD package Fluent6.2 (Lebanon, NH, USA). The second order upwind discretization scheme was used for the momentum, turbulent kinetic energy, turbulent energy dissipation rate, species (for passive scalar model), and volume fraction (for Mixture model and Eulerian model). The SIMPLE algorithm (semi-implicit method for pressure-linked equations) was employed to evaluate pressure-velocity coupling. The solution was iterated until convergence was achieved, such that the residual for each equation fell below 10^{-4}. Improving the convergence criterion to 10^{-5} had a negligible effect on the simulation results. In the DPM model, the airflow field was first simulated, and then the trajectories of individual particles were evaluated using particle equation of motion given in (12). The Euler implicit method, a stable scheme regardless of the value of integration time step, was used to solve the particle equation (12). Note that the Euler implicit scheme is second order method, thus the schemes adopted in the present study are second-order accurate in time and space.

All the grids listed in Table 3 are generated with Gambit2.3. Non-uniform grids are used with finer grids in the near wall regions. The wall unit values (*y*^{+} = *u*_{τ}*y*_{p}/_{v}, where *u*_{τ} is the wall friction velocity) are less than 1. The maximal difference of mean velocity between Medium and Fine meshes (along 8 axial positions) are about 5% simulation respectively. Thus only results from Medium meshes would be presented in the following sections. The converged solutions are achieved when scaled residuals of all the variables were less than 10^{-6}. Integrated quantities such as skin friction coefficients are also used to check the convergence. It is found that scaled residuals less than 10^{-6} are enough to guarantee the converged solutions.

## Numerical Simulations and Model Assessment

### Airflow Patterns and Mean velocity

Figure 4 shows the airflow patterns at the symmetry plane for both flow visualization.^{16} particleimage velocimetry measurement,^{19} and numerical simulations. The smoke visualization was obtained by injecting titanium tetrachloride smoke at different locations of the room.^{16} It can be seen that numerical prediction correctly imitates the flow behavior observed in the experiments: incoming air first attached to the celing due to Coanda effect after enterning the room, then it travelled along the ceiling for a certain distance before it separated from the ceiling and air below the jet was entrained so that a recirculation flow pattern (vortex) formed below the jet. Coanda effect resulted from the pressure difference between the two sides of the airjet and numerical pressure distribution given in Figure 5 clearly shows smaller pressure in the top left corner than outside air. Comparising the numerical results in Figure 4 and Figure 5 shows that vortex center (near zero velocities) from flow pattern is almost identical to that from pressure distribution (local minimum pressure), which is consistent with.^{33}

Plots of the mean velocity determined from simulation and hotwire measurements^{16} at seven streamwise (x) positions on the symmetry plane of the room are shown in Figure 6. The general agreement is quite good and simulation results capture the general trends of the experiemntal data. In fact, excellent agreement is achieved at first two positions (x/L = 0.125, 0.25). At the other positions, discrepancy between experimental data and simulation results mainly happen near the floor regions.

### Comparisons of Contour plots

Figure 7A shows the particle spatial concentration from experimental measurements. The contour plot is interpolated from the particle concentrations at 25 measurement points. It can be easily seen that the particle spatial distribution is strongly non-uniform: the highest concentration is near the lower left corner while the lowest concentration is near the upper right corner of the room.

Figure 7B gives the spatial distribution of particle concentration predicted by PSM. In some parts of the room, the pattern predicted by PSM is similar to that of the experiment. PSM correctly predicts the the highest and lowest concentration regions but much smaller values than measurements. In the other regions, the discrepancy between prediction and experimental data is nonnegligible.

The prediction by discrete particle model (DPM) is shown in Figure 7C. Simulation imitats experimental values with fairly good accurary. In particular, the highest and lowest concentration regions predicted by DPM agree with those of the experiments. Compared to PSM, DPM predicts more accurate particle concentration in almost all the regions at the symmetry plane.

Figure 7D and E depict the predictions by the Mixture model and Eulerian model respectively. It can be noticed that the former predicts lower values than experimental data while the latter overpredicts the particle concentration in most regions. Morevoer, the largest concentration region predicted by Mixture model is located in the middle region near the floor (2.5 m < X < 3 m) instead of left low coner.

### Model Assessment and validation

In this paper, ASTM D5157-97 is adopted as the standard guide for statistical evaluation of numerical models. Six parameters are chosen to assess the performance of models, ie, correlation coefficient (CC), regression slope (RS), regression intercept (RI), normalized mean square error (NMSE), fractional bias (FB), and fractional variance bias (FVB). The first three parameters reflect the general agreement between simulated and experimental concentrations while the other two parameters indicate the bias in the mean or variance of simulated concentrations relative to experimental concentrations. It should be noted that CC could only detect linear dependencies and zero value cannot prevent the existence of dependence in other types. Linear relationship between simulated and experimental values could also be measured with regression slope (RS) and regression intercept (RI): best fit linear relationship means RS of one and RI of zero. In case that systematical differences exist between simulated and experimental data, CC would take value near one but RS and RI would refelect the systematic differences. The differences between experimental and simulated values also could be measured by NMSE, which would have a value of zero for a perfect agreement and tend to toward higher values when the differences increase. Fractional bias (FB) and fractional variance bias (FVB) are used as measure prerformance tools by assess the mean (FB) or vairiance (FVB) bias between experimental and simulated values. These two parameters are symmetrical and bounded: perfect agreement would have zero values of FB and FVB while larger differences would make FB and FVB towards 2 (extreme overprediction) or -2 (extreme underprediction). Values of ±0.67 for FB and FVB are equivalent to underprediction (positive value) or overprediction (negative values) by a factor of two.

Table 4 presents the statistical analysis of the four models performance. According to the standard guide of ASTM, a model performance will be considered as adequate when the six criteria listed in the last row of table 1 are satisfied. It can be seen that, at the ventilation rate 19.5 ACH, none of the four models completely satisfies the ASTM criteria. The worst case is the performance of PSM: it fails to satisfy all six criteria. On the contrary, DMP performs the best among the four models: only one parameter (RI) exceeds the value given by ASTM criteria. The performance of Eulerian and Mixture models are similar: two criteria are satified while the other four ones are not.

The failure of PSM model has been expected, as noted in the previous sections, the underlying assumption in PSM is that the particle behaviors are totally determined by airflow and inerial effects of particles are negelcted. As a consequence, PSM cannot calculate the turbulent diffusion of the particle near the sources and information of velocity fluctuation correlation is lost. Note that the diameters of the particles in the experiments^{20} span a wide range (Fig. 2) and larger particles would not follow the airflows instantly. Stoke number, which presents an imporant criteria to investigate whether particles are in kinetic equilibrium with surrounding airflows, is defined as the ratio between particle relaxation time and airflow relaxation time:

*U*

_{f}and

*L*

_{f}are airflow characteristical velocity and length scales. The Stokes numbers of particles used in experiments

^{12}rang from 0.0002 to 0.2.

Larger Stokes number indicate particles can not follow the airflows closely due to their inertial effects. Lemaire^{34} showed the neglect of inertial effects is only valid for time scales larger than the turbulent integral time scales.^{35} In contrast, DPM does not have this limitation and thus predictions based on DMP agree much closer with experimental data. It should be noticed that RI of DPM slightly exceeded the ASTM criteria.

The failure of Mixture and Eulerian models in the present predictions of particle concnetrations may be due to the low volume fraction of particles and existence of larger particles (larger Stokes number). Although the volume fraction of particles near the emission source was relatively high^{20} the overall voluem fractions were much less than 1%. When the voluem fractions of particles are low, predictions based on Eulerian and Mixture models are always unreliable.^{11} Furthermore, Boffetta et al^{36} showed that Eulerian model would break down at larger Stokes number (≥0.1, the atual number depending on the flow charateristics of carrier phase).

## Table 3

Nx, Ny, Nz denote the number of cells in stream-wise, wall-normal and spanwise direction respectively.

To further validate the four models in the prediction of particle concentration, two higher ventilation rates (28ACH and 66ACH) in the same room are also investigated (Table 4). Mean velocity simulated by the same turbulence model^{22} agree well with correspoinding measurements from Hotwire^{16} and PIV.^{37} Statistical performance of DPM with experimental data from^{12} is listed in Table 4 and results show that DPM predictions achieved excellent agreement with measurement. Wang et al^{17} developed a 2-D mathematical model (zonal model) based on the mass balance of particulate matters in the airspace which is divided into small element zones. The central assumption of this model is that the particles are well mixed in each element zones. The zonal model results at the two higher ventilation rates are listed in Table 4 and it can be found that the agreement with measurement data is not good especially at highest ventilation rate. In fact, significant differences in the over mean particle concentration pattern were observed between measured data and zonal model prediction.^{17} The assumption of well mixed element zone in zonal model is the most possible reason for the disagreement.

## Table 4

Statistical evaluation of the four multiphase flow models.

Table 5 shows the comparison of average particle concentrations in the whole room between the measurements and numerical simulations for three ventilation rates. As expected, only DPM predict values very close to experimental ones.

## Discussion

This section is devoted to the investigation of several factors which may affect the accurary of numerical predictions of particle concentration.

Effects of turbulence models of airflows. Undoubtly particle behavior and spatial concentration are closely related to surrounding airflows, thus accurately predicting airflow characteristics is necessary and essential for particle predictions. It is known one drawback of k-∊ turbulence models is that they cannot predict the anisotropy effects in the three dimensioanl wall jet flows^{38} damping of turblence flucuations perperndicular to the wall is absent. As a consequence, the growth rate of wall jet in the direction of parallel surface is predicted much smaller than measurements.^{39} On the other hand, it has been shown that Reynold stress transport models (RSTM) with wall reflection terms^{40} could predict the growth rate of wall jet correctly. Thus it is useful to investigate the differences between RSM and k-∊ models in the predictions of flow characteristics and particle concentration. Mean velocity from RSM simulation is included Figure 6 and it is obivious that the differences between two simulated results are very small. Same conclusion has been obained for the particle concentration predictions for the three ventilation rates. In fact, RSM do not show much improvement in the predictions of both mean velocity and particle concentration.^{41} However, it should be emphasized that both mean velocity and particle concentration were measured and calcuated in the middle plane of the room. In the regions away from the middle plane, nonneglible differences for both mean velocity and particle concentration existed between two models predictions.

## Table 5

comparison of mean particle concentrations.

Effects of boundary conditions. Proper settings of boundary condition in the CFD simulations are necessary for successful predictions. For example, Lee et al^{29} found that profiled inlet profile from experiments showed better agreement qualitatively and quantitaviely with measured particle concentration than uniform inlet velocity. However, in the present study, simualtions based on uniform inlet velocity (estimated from profiled inlet velocity of experiments) were also conducted and the differences in both mean velocity and particle concentration were very small compared with uniform velocity profile except the inlet regions (x/L < 0.1, y/H < 0.25). The underlying reason should be ascribed to the relative small inlet height compared to the room height (h:H = 0.021). On the other hand, the treatement of particle reaching the wall is very critical for the present case. Prelimary simulations based on standard k-∊ turbulence model along with wall functions^{23} showed the trappring of particles near the wall are not suitable when near wall regions are not resolved (20 < y + <100 in that simulation). Zhang and Chen^{8} suggested to set restitution coefficients to a small values when using high Reynolds number turbulence models. By doing so, the particles are immediatedly stopped but not trapped when reaching the walls. They would escape from the walls and resupend in the airflows once they obtain enough normal velocity. However, preliminary simulations showed somewhat strong dependence of particle concentration on the values of restitution coefficients. Furthermore, it was found that predictions based on different restitution coefficients for ceiling, floor and side walls sometimes agreed better with experimental data than same coefficients. It seems that there should exist certain relationship between these coefficients and further study is needed. For the reasons list above, the present study adopts the two layer k-∊ turbulence model without using the wall functions by fully resolving the wall regions. But the treatment of celing is different from other walls, as hinted from simulations with high Reynolds k-∊ model, restitutioin coefficients are not set to zero. Several simulations were tested and it was found a value of 0.7 generated particle concentration closest to experimental data.

Effects of particle size distribution. The simulated particle concentrations in the room was found to be highly dependent on particle size distribution. Preliminary simulations with DPM model showed predictions based on single size (1.63 µm, mean particle size) gerenrated poor agreement with expeirmental data. As seen from Figure 2, the particle size distribution used in the experiment covered a wide range and particles with different sizes would behave very differently. Thus this study shows that using the real particle size distribution in the calculation can provide more reliable estimates.

## Conclusion

This investigation studied the particle dispersion and spatial distribution in a full scale ventilated room by numerical simulation. The numerical study was achieved by applying two turbulence models (two layer k-∊ model, Reynolds stress model), four multiphase models (PSM, DPM, MIX, EUL). Corresponding experimental results were used to evaluate the performance of the four models and ASTM guide^{15} was used as validation criteria. The following conclusions can be drawn from this study:

1. Among the four multiphase models, only results from DPM agreed very with corresponding experimental data and satified the validation criteria.

^{15}For this model, it is imporant to include all particle forces which are expected to affect the particle trajectories;2. PSM is not suitable for modeling larger particles in the room due to the its neglect of particle settling velocity and inerial effects. The application of MIX and EUL models in indoor air study should be careful especailly when the volume fractions of particles are less than 10%;

3. The effects of velocity distribution at the inlet on the airflows and particle concentration inside the room depend on the the relative size of inlet. For the present case with a small inlet, uniform velocity and profiled velocity from experiments generated almost the same flow patterns and particle spatial distribution. On the other hand, detailed information for the particle size distribution and particle behaviors at the walls are critical for the accurary of simulations.

## Disclosures

Author(s) have provided signed confirmations to the publisher of their compliance with all applicable legal and ethical obligations in respect to declaration of conflicts of interest, funding, authorship and contributorship, and compliance with ethical requirements in respect to treatment of human and animal test subjects. If this article contains identifiable human subject(s) author(s) were required to supply signed patient consent prior to publication. Author(s) have confirmed that the published article is unique and not under consideration nor published by any other publication and that they have consent to reproduce any copyrighted material. The peer reviewers declared no conflicts of interest.