key: cord-1040749-tep4h5l1 authors: Ding, L.; Lai, A.C.K. title: An efficient lattice Boltzmann model for indoor airflow and particle transport date: 2013-05-01 journal: J Aerosol Sci DOI: 10.1016/j.jaerosci.2013.04.004 sha: 7fda5d920eac31455a74e43611957a5aaac5fcac doc_id: 1040749 cord_uid: tep4h5l1 The three-dimensional multiple-relaxation-time LB (MRT-LB) and Lagrangian particle tracking methods were applied to simulate turbulent airflow and particle dispersion in a ventilated room with a partition. The turbulent airflow was simulated by large eddy simulation (LES) using the MRT-LB method with the Smagorinsky model. This method was verified by comparing it with experimental and other numerical results. Good agreement was observed between airflow simulation and experimental data. It is also demonstrated that the LES carried out by the MRT-LB method can produce airflow results very similar to the RNG LES and provide better prediction than the standard and RNG k–ε models. In order to further improve the efficiency of the MRT-LB method, the multi-block grid refinement (MBGR) technique was used. The accuracy and efficiency of the MBGR and the consistency of physical quantities across the interface were investigated. In simulation of particle dispersion in the model room, particles with diameters of 1 and 10 µm were considered. It is shown that this model can successfully capture dispersion characteristics of particles. After the outbreak of the Severe Acute Respiratory Syndrome (SARS) and Avian Influenza (H5N1), indoor air quality (IAQ) has been attracting increasing attention since it not only adversely impact human health (Risom et al., 2005; Weichenthal et al., 2007) and productivity (Bakó-Biró et al., 2004) , but also, if not managed properly, the public at large could be affected. Several experimental and numerical investigations have been carried out to gain in-depth understanding of IAQ (Lai, 2002; Zhang & Chen, 2006) . Experimental studies have provided direct evidence of airflow and particle transport phenomena and also benchmarks for verification of computational models. The complexity of indoor airflow makes high-quality experimental measurement not only difficult but also expensive, especially when the room is occupied (Jiang et al., 2009; Rey & Velasco, 2000) . Furthermore, if particles with various sizes are involved, expensive instruments have to be utilized. With advances in computer technology, CFD modeling has become another powerful tool widely adopted in IAQ research under various physical and geometric conditions (Tian & Ahmadi, 2007; Wang et al., 2012; Zhang & Chen, 2007; Zhao et al., 2004 Zhao et al., , 2008 . Compared with experiments, CFD modeling is better for conducting systematical parametric or sensitivity analysis. It can also overcome some deficiencies of experiments and extends the range of research. Since substantial computational resources are required when the CFD is applied to model the practical scenario, highly efficient numerical methods need to be developed. where subscripts i; j; k∈f0; 1; 2; :::; 18g represent the direction of discrete velocities. e i is the discrete velocity vector for the D3Q19 model, m k (x, t) are the moments of f i (x, t). M ij is the 19-by-19 transformation matrix and M −1 ij is the inverse of M ij .Ŝ jk is the collision matrix in moment space. It is a diagonal matrix with the form S ¼ diagðs 0 ; s 1 ; s 2 ; :::s 18 Þ ð2Þ where s i , for i¼0, 1, …, 18, specifies the relaxation frequencies for all moments toward their corresponding equilibria m eq k ðx; tÞ. Here, s 0 , s 3 , s 5 and s 7 , corresponding to conserved moments in the hydrodynamic limit, i.e. density and momentum, are set to zero; s 9 , s 11 , s 13 , s 14 and s 15 are related to kinematic viscosity ν, namely, The rest of the relaxation frequencies obtained through linear analysis are as recommended by d 'Humières et al. (2002) . In order to simulate the turbulent flow with LES, the Smagorinsky model is incorporated into the MRT-LB framework. According to Yu et al. (2006) , an additional viscosity called the turbulent eddy viscosity ν t is introduced to model the turbulence. Hence, the total viscosity ν total can be expressed as the sum of ν and ν t The turbulent eddy viscosity can be expressed as where C s denotes the Smagorinsky constant and is set to be 0.16, as suggested by Krafczyk et al. (2003) ; Δ represents the lattice spacing, S αβ (α; β∈fx; y; zg in Cartesian coordinates) is the strain rate tensor, i.e., In the MRT-LB method, S αβ can be computed directly from non-equilibrium moments (m neq k ðx; tÞ ¼ m k ðx; tÞ−m eq k ðx; tÞ). They can be expressed as Details of the derivation of Eq. (6) can be found in Yu et al. (2006) . Then, the corresponding relaxation frequencies were modified in terms of Eq. (3). As mentioned above, a multi-block grid refinement scheme, in accordance with Yu et al. (2002) , was employed. Only the basic theory of the method is demonstrated in the following. For simplicity, a two-block grid system in a plane is shown in Fig. 1 . The ratio of lattice spacing between the coarse and the fine block is defined as b r , i.e., b r ¼ Δ c =Δ f (superscripts c and f are used to represent variables on the coarse and fine blocks, respectively). In order to keep the entire flow field consistent, the Reynolds number Re should be the same in different blocks involving different lattice spacings, i.e., Re c ¼Re f . The convective refinement is chosen which means the characteristic velocities on the coarse and fine block are equal. Thus the relaxation times τ c and τ f must follow the following relationship: More attention should be paid to the way to deal with boundaries of coarse and fine blocks in order to make the variables and their derivatives continuous across the interface between the two blocks. For the MRT-LB method, the rule of data transfer from the fine block to the coarse block through the coarse block boundary is where subscripts l; n; p∈f0; 1; 2; :::; 18g.Ĉ kl is a diagonal matrix with the form Similarly, the rule on the fine block boundary is Details of the multi-block grid refinement strategy in the MRT-LB method are provided in Appendix A. Simulations of airflow in a model room with a partition, experimentally studied by Posner et al. (2003) , were performed in this paper. The schematic of the computational domain is shown in Fig. 2 . The dimensions of the model room in x, y and z directions are 0.914 m, 0.305 m and 0.457 m, respectively. A partition with half the room height, i.e., 0.1525 m, is located in the middle of the room. Ventilation air comes into the chamber vertically and exits through the outlet. The inlet and outlet are on the ceiling and dimensions of both are 0.101 m  0.101 m. The Reynolds number based on the vertical inlet velocity u in of 0.235 m/s and inlet size is 1628. The airflow properties were kept constant with values corresponding to a room temperature of 23 1C. The air density was 1.18 kg/m 3 and the kinematic viscosity was 1.46  10 −5 m 2 /s. It should be noted that the values in physical units were converted into lattice units for the computation (see Appendix B). A uniform cubic grid was used to discretize the computational domain. Based on the characteristics of the problem, a uniform velocity profile was assumed for airflow entering into the room, along with no-slip boundary condition for the walls and the partition. Implementation of such boundary conditions in the LB method was as follows. A velocity Dirichlet boundary condition, introduced by Skordos (1993) , was imposed on the inlet, walls and the partition while a constant pressure boundary condition was imposed on the outlet. The simulations were performed in parallel mode on two Dell Precision T7500 Workstations (Dual Intel Quad-Core Xeon E5520 2.26 GHz, 48 GB RAM each) with 14 cores. Since only airflow characteristics were investigated by Posner et al. (2003) , the corresponding particle results obtained through simulation by Tian et al. (2006) were utilized for verification. The particle motion was solved by the Lagrangian particle tracking method (see Appendix C). Particles with diameters d of 1 μm and 10 μm and a density of 800 kg/m 3 were considered. Their particle relaxation times were 3.01  10 −6 s and 2.62  10 −4 s, separately. At the beginning, only the airflow was simulated. After 70 s, 144 evenly distributed monodisperse particles were injected at the inlet with initial velocities identical to the ambient airflow. The particle injections were subsequently stopped when the simulation time reached 100 s. Eighty six thousand four hundred particles were introduced in the model room. As mentioned in Tian et al. (2006) , when the particles reached the wall surface, the process of bouncing was taken into consideration with a particle wall impact model where the normal and tangential coefficients, which defined the amount of momentum retained by the particles after collision in the direction normal and parallel to the wall, were assumed to be 0.9. The main parameters of airflow computation are listed in Table 1 . For the grid independence analysis, four different grid sizes, coarse (MRT-LB-LES1), medium (MRT-LB-LES2), fine (MRT-LB-LES3) and finer (MRT-LB-LES4), were employed. The time-averaged velocity vectors and their magnitude distributions in the mid-plane of the model room (z¼0.2285 m) are shown in Fig. 3 . It can be seen that the four meshes produce results similar to each other in spite of some differences between the coarse, medium and the other two finer grids in the vicinity of the wall regions. To further analyze the independence on these four grid sizes, the grid convergence index (GCI) (Roache, 1994) was used. Using the solution of MRT-LB-LES1 as the reference, the GCI for MRT-LB-LES2-MRT-LB-LES1, MRT-LB-LES3-MRT-LB-LES1 and MRT-LB-LES4-MRT-LB-LES1 are 1.62%, 4.68% and 4.10%, respectively, where they are calculated based on values at 180 selected points distributed uniformly in the computational domain. All of them are less than 5%, indicating that the accuracy of airflow results using the coarse grid is also acceptable, especially for regions away from the wall area. On the other hand, it is noticed that the GCI for MRT-LB-LES2-MRT-LB-LES1 is only 1.62%, indicating that the improvement of airflow field is not obvious when the grid resolution is refined from MRT-LB-LES1 to MRT-LB-LES2. When the grid is refined further, improvement of the results, especially for the near-wall area, can be observed. Since the near-wall grid resolution is of crucial importance for simulation of wall-bounded turbulent flow, which has a great effect on behavior of nearwall particles, especially for submicron particles (Lai & Chen, 2006) , the fine grid was used for the rest of analysis which represents a compromise between numerical accuracy and computational cost. It should also be noted that the near-wall resolution for the uniform grids of MRT-LB-LES1 and MRT-LB-LES2 may not meet the requirement for particle simulation, nevertheless, they may work well when the multi-block grid refinement is introduced where only the near-wall area is refined. The investigation of the applicability of multi-block grid refinement is presented in Section 4.1.3. In order to reduce the side effects of coarse block resolutions on the final results, the resolution of MRT-LB-LES2 is chosen as the resolution of coarse block in this paper. Figure 4 presents the comparison of time-averaged inlet jet centerline axial velocity predicted by the MRT-LB method and the experimental data (Posner et al., 2003) as well as other numerical results (Tian et al., 2006) . Good agreement is achieved between the MRT-LB predictions and measurements as well as the RNG LES results. The MRT-LB method performs better than the standard and RNG k-ε models, especially in the region of 0.15-0.25 m. Airflow enters the model room with an axial velocity of 0.235 m/s. The velocity increases slowly along the distance until about 0.17 m from the inlet according to the experimental results, 0.13 m and 0.14 m according to the standard and RNG k-ε models and 0.17 m according to the RNG Comparison of time-averaged vertical velocity along the horizontal line at mid-partition height between all the numerical models and measurements is shown in Fig. 5 . From x ¼0 m to the partition, all numerical models yield results that are almost similar to the experimental data, while the MRT-LB method produces slightly better predictions in the region from x ¼0.13 m to the partition. The first peak velocity near the partition is 0.109 m/s measured by the experiment, about 0.114 m/s by the standard k-ε model, 0.102 m/s by the RNG k-ε model, 0.138 m/s by the RNG LES and a slightly larger value of 0.171 m/s by the MRT-LB method. The absolute difference between the MRT-LB method and the measurement is 0.062 m/s. The vertical velocity then decreases to zero very rapidly and remains low over a short distance. After that, it increases very fast in the opposite direction because of the inlet jet. The peak value, found in the region right beneath the inlet, is −0.273 m/s from the experiment, about −0.19 m/s from the standard k-ε model, −0.205 m/s from the RNG k-ε model, −0.219 m/s from the RNG LES and −0.222 m/s from the MRT-LB method. It is further noted that the MRT-LB method successfully captures the velocity variation and yields a peak velocity closest to the measurements. The performance of the standard and RNG k-ε models are close to each other but they do not match the measurements well. The velocity profiles they predict are less steep on the sides and rounder at the bottom compared with the MRT-LB method and the RNG LES method. There is another peak near the right wall. The measured value is 0.115 m/s, while is about 0.089 m/s, 0.087 m/s, 0.115 m/s and 0.151 m/s found by the standard k-ε model, the RNG k-ε model, the RNG LES and the MRT-LB method, respectively. Although the MRT-LB result is still a bit larger than the measurement, the absolute difference decreases to 0.036 m/s. In general, although the two turbulence models (the standard and RNG k-ε models) capture the general trends of the experimental data, a better prediction of the airflow is achieved through the RNG LES and the MRT-LB methods. The time-averaged velocity fields predicted by the MRT-LB method and other numerical models in the mid-plane of the model room are given in Fig. 6 . It is shown that the MRT-LB method successfully produces results similar to the other three models, especially for the strong recirculation zone near the right wall. As mentioned above, in contrast to the standard and RNG k-ε models, which are inherently time-averaged approaches for turbulence simulation, one of the advantages of LES is its ability to capture the instantaneous turbulences. Hence, besides the time-averaged airflow conditions, the root-mean-square (RMS) velocity, which is a measure of the magnitude of time-dependent fluctuating velocity, is also employed to estimate the quality of LES using the MRT-LB method. A comparison of RMS velocities simulated by the MRT-LB method and the other three numerical models is shown in Fig. 7 . It can be clearly seen that the MRT-LB method produces results very similar to the RNG LES, especially in the left zone of the model room as well as the inlet jet flow region. They predict significantly higher (one order) RMS velocities than the standard and RNG k-ε models. This indicates that the LES using MRT-LB method can perform as well as the RNG LES and can provide better predictions of the turbulent airflow characteristics than the time-averaged approach. In terms of characteristics of wall-bounded turbulent flows, the near-wall area in the computational domain was refined. Two kinds of grids, i.e., fine and coarse, were utilized. The fine grid was employed in the region within 0.03 m from the wall boundary, while the coarse grid was applied in the remaining area of the domain. As mentioned above, the resolution of the coarse grid was chosen to be the same as the MRT-LB-LES2. The ratio of lattice spacing between the two grids b r ¼2. We refer to this computation as MRT-LB-LES-MBGR hereinafter. The schematic of grid arrangement is shown in Fig. 8 . For the sake of clarity, the grid shown is 10 times as large as the actual one in the mid-plane of the model room. The time-averaged velocity distributions in the mid-plane are presented in Fig. 9 . Compared with Fig. 3(b) , some improvement in accuracy can be observed near the fine grid region. For example, velocity distributions near the top of the partition become very close to the ones in Fig. 3(c) and (d) . Velocity distributions near the bottom of the left zone are also getting better than the ones in Fig. 3(b) . The dashed lines in the figure represent the interface between the fine and coarse grids. It is worth pointing out that the velocity contours are smooth across the interface, indicating that the consistency in the transfer of the physical quantity information across the interface is maintained. The RMS velocity distribution is given in Fig. 10 . A good agreement Tian et al. (2006).) with the corresponding ones in Fig. 7 can be seen and the smoothness of the RMS velocity contours is also ensured. After the accuracy of MBGR has been verified, its computational efficiency is also investigated. As in our recent paper (Ding et al., 2012) , the average time cost per iteration of each resolution was employed as the criterion and the results are shown in Table 2 . The refinement strategy results in nearly 3.0  10 7 elements, about 1.1 times as many as the number of MRT-LB-LES4. Meanwhile, the near-wall grid size is only 0.75 times as long as the MRT-LB-LES4. Although the average time cost is a little larger than that of MRT-LB-LES4, it should be mentioned that if the uniform grid with the same resolution as the fine grid of MBGR is used, the total number of elements will be nearly 6.5  10 7 , about 2.4 times MRT-LB-LES4. The time cost will also increase with respect to the total number of elements. Since the resolution near the wall is quite crucial to the simulation of particle-wall interaction (e.g., deposition), it is much more efficient to use the MBGR technique than the uniform grid in IAQ studies. Furthermore, our work in this paper is just an attempt to test the applicability of MBGR and, therefore, a simple refinement strategy was employed and the scope of refined area in terms of flow characteristics has not been systematically investigated. For occupied rooms, distributions of airflow intensity are not uniform. Different areas can be treated with different refinement strategies. Furthermore, the cut-off limit of LES on coarse grid resolution has not been studied. It is still possible to reduce the resolution of coarse grid, e.g., just as the resolution of MRT-LB-LES1. Obviously, substantial saving and better efficiency can be achieved using the MRT-LB with MBGR. The computation of particle dispersion in the model room was performed based on the MRT-LB simulated airflow conditions. Figure 11 shows the comparison of the number of tracked particles in the model room obtained from the current model and other numerical models for particle size d¼1 μm. It can be seen that the number of tracked particles increases with respect to time when the particles were injected into the model room continuously, while it decreases gradually after 100 s when particle injection was terminated. At 85 s, 15 s after the beginning of particle injection, all numerical results predict almost the same number of tracked particles. From 100 s to 130 s, more particles are predicted through the simulation using RNG LES than the standard and RNG k-ε models, while after 145 s, the number produced based on the RNG LES is less than the other two models. Since the three-dimensional transient turbulent velocity field, which strongly affects particle transport, can be adequately simulated by LES, it is anticipated that the simulation based on the LES can predict more accurate particle dispersions than the time-averaged models. This is particularly helpful for the simulation of transient events, e.g., sneezing and coughing . Generally, the present model provides results consistent with RNG LES, although smaller numbers are obtained after 100 s, it approaches the RNG LES based results after that. A further comparison was carried out. The particle concentration simulated by the current model and the other three models in the mid-plane at 130 s is shown in Fig. 12 . The particle source in cell (PSI-C) method (Zhang & Chen, 2007) was employed to transform the Lagrangian trajectory information into the form of concentration distributions. It can be seen that more particles have dispersed into the left part of the model room at this moment and very few particles can be found inside the inlet jet flow region. The present model provides results similar to the other three models. Furthermore, concentrations predicted by the current model and the simulation using RNG LES is higher than the standard and RNG k-ε models in the middle of the left zone, while a lower concentration is obtained through the present model in the lower left corner of the partition which may partially contribute to the smaller number of tracked particles shown in Fig. 11 . The time evolution of particle concentration from 70 s to 160 s is also presented as an animation. Comparisons of the number of tracked particles in the model room for particle size d ¼10 μm are presented in Fig. 13 . Before 100 s, the number of tracked particles also increases along with time and peaks at the end of the particle injection. This is very similar to the corresponding stage for particle size d ¼1 μm. All the models track nearly the same number of particles at 85 s. The number begins to decrease after 100 s, while simulations based on the time-averaged methods show a much more significant drop than the LES method and larger discrepancies can be observed between the two kinds of models. The present model presents a better agreement with the simulation using the RNG LES compared with the d¼1 μm case. It is also illustrated that more particles with a diameter of 10 μm are left in the model room than its counterpart at the decay stage. The dispersion characteristics of particles are analyzed as follows in order to find the possible reason that cause the differences. At the particle injection stage, particles injected into the model room are mainly concentrated in the right zone. Only a small portion disperses into the left zone and very few particles exit through the outlet. Therefore, the number of particles in the model room is mainly determined by the injection rate at the inlet. Since the injection rates for the two kinds of particles are the same, they show a very similar variation before 100 s. At the decay stage, as can be seen in Fig. 12 , more and more particles disperse into the left zone and then the number of exiting particles increases, which reduces the total number of particles in the model room. As the particle injection is stopped after 100 s, the ventilation rate then plays the dominant role in determining the number of particles in the model room. As discussed in our recent paper (Ding et al., 2012) , the effect of ventilation on accumulation mode particles is much more significant than on coarse mode particles, resulting in less particles remaining in the model room for particle size d ¼1 μm than for d ¼10 μm. The time evolution of particle concentration for particle size d ¼10 μm from 70 s to 160 s is also presented as an animation for comparison. Finally, one additional feature of LES should be addressed. Due to its nature, there is no need to introduce any discrete random walk (DRW) model to determine the fluctuating velocity components, in contrast to the time-averaged methods. Previous studies (MacInnes & Bracco, 1992) have indicated that the DRW model may give nonphysical results in strongly nonhomogeneous flows. The three-dimensional MRT-LB and Lagrangian particle tracking methods were applied to simulate the turbulent airflow and particle dispersion in indoor environments. The LES was used in this paper. It is different from the methods conventionally used (e.g., the standard and RNG k-ε models) in IAQ investigations, which modeled the turbulence from the time-averaged point of view. It was implemented by using the MRT-LB method with the Smagorinsky model. Experimental and numerical studies of airflow characteristics and particle dispersion in a ventilated model room with a partition were employed to verify the present model. Good agreement between the present results and the experimental data were observed. The MRT-LB method with the Smagorinsky model was able to capture the turbulent airflow characteristics. It was also demonstrated that LES performed by the MRT-LB method can produce airflow results very similar to the RNG LES. Furthermore, it can give better prediction than the standard and RNG k-ε models, showing the advantages of LES computation of turbulent airflow in indoor environments over the time-averaged methods. In order to enhance the efficiency of MRT-LB method, the multi-block grid refinement technique was used. The accuracy of MBGR was verified through comparison. The consistency of physical quantities across the interface between grids with different resolutions was also confirmed. Compared with a uniform grid, the MBGR can produce higher-resolution grids in the regions of interest (e.g., near-wall areas) at a lower computational cost. The MBGR thus provides further improvements of the efficiency of the MRT-LB method and will be very useful for further studies of practical problems in IAQ. In our simulation of particle dispersion in the model room, two sizes of particles, i.e., 1 μm and 10 μm, were considered. The computations were performed using the Lagrangian particle tracking method based on the airflow conditions provided by the LES with the MRT-LB method. It is shown that the current model can successfully capture dispersion characteristics of the particles. With the development of high performance computing, particle computation based on the LES should be carried out in IAQ studies and LES computations using the MRT-LB method is a better choice because of its higher accuracy and efficiency. by the number of meshes N pu , then the relationship of space step Δ between physical units and lattice units can be given as where Δ lu is often set to be unit in the LB method. Meanwhile, according to the definition of ν in the LB method, the relationship of ν between physical units and lattice units can be expressed as In order to make the flow simulation equivalent between the two unit systems, the Reynolds number should be the same, i.e. The relationship of characteristic velocity between the two unit systems can be determined when substituting Eqs. (B1) and (B2) into Eq. (B3) Then the relationship of time scale and acceleration between the two unit systems can be obtained through the following equations: Effects of pollution from personal computers on perceived air quality, SBS symptoms and productivity in offices Lattice Boltzmann method for fluid flows Multiple-relaxation-time lattice Boltzmann models in three dimensions Numerical study on particle dispersion and deposition in a scaled ventilated chamber using a lattice Boltzmann method Large eddy simulation of turbulent open duct flow using a lattice Boltzmann approach Grid refinement for lattice-BGK models A lattice subgrid model for high Reynolds number flows Investigation of particle dispersion and deposition in a channel with a square cylinder obstruction using the lattice Boltzmann method Experimental and numerical study of airflows in a full-scale room Large-eddy simulations with a multiple-relaxation-time LBE model Particle deposition indoors: a review Modeling particle deposition and distribution in a chamber with a two-equation Reynolds-averaged Navier-Stokes model Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability Dispersion and deposition of spherical particles from point sources in a turbulent channel flow Lattice Boltzmann method on composite grids Stochastic particle dispersion modeling and the tracer-particle limit Measurement and prediction of indoor air flow in a model room Generalized lattice Boltzmann equation with forcing term for computation of wall-bounded turbulent flows Lattice-Boltzmann approach for description of the structure of deposited particulate matter in fibrous filters Lattice BGK models for Navier-Stokes equation Experimental study of indoor air quality, energy saving and analysis of ventilation norms in climatised areas Oxidative stress-induced DNA damage by particulate air pollution. Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis Perspective: A method for uniform reporting of grid refinement studies Experimental and numerical investigation of interpersonal exposure of sneezing in a full-scale chamber Initial and boundary conditions for the lattice Boltzmann method Parallel simulation of particle suspensions with the lattice Boltzmann method The Lattice Boltzmann Equation for Fluid Dynamics and Beyond Particle deposition in turbulent duct flows-comparisons of different model predictions On the numerical study of contaminant particle concentration in indoor airflow Advanced turbulence models for predicting particle transport in enclosed environments Indoor ultrafine particles and childhood asthma: exploring a potential public health concern Viscous flow computations with the method of lattice Boltzmann equation A multi-block lattice Boltzmann method for viscous fluid flows LES of turbulent square jet flow using an MRT lattice Boltzmann model Experimental measurements and numerical simulations of particle transport and distribution in ventilated rooms Comparison of the Eulerian and Lagrangian methods for predicting particle transport in enclosed spaces Particle dispersion and deposition in ventilated rooms: Testing and evaluation of different Eulerian and Lagrangian models Comparison of indoor aerosol particle concentration and deposition in different ventilated rooms by numerical method The work described in this paper was partially supported by a grant from the Research Grants When the density distribution function is divided into two parts, i.e., the equilibrium part and the non-equilibrium part, Eq.(1) givesFor the coarse block, Eq. (A1) givesSimilarly, for the fine block,Since the velocity and density must be continuous across the interface between the two blocks, it is seen thatTo maintain continuity in deviatoric stress, it is required thatwhere I is the identity matrix. Then it can be obtained from Eq. (A5) thatwhereĈ kl can be expressed by Eq. (9) in Section 2.2. Substituting Eqs. (A4) and (A6) into Eq. (A2) one obtainsSimilarly, substituting Eqs. (A4) and (A6) into Eq. (A3) one obtainskl is the inverse ofĈ kl . Then Eqs. (8) and (10) where u pi is the particle velocity, u fi is the airflow velocity at the particle position, τ p is the particle relaxation time, ρ f and ρ p are the air and particle density, respectively, g i is the gravitational acceleration, x pi is the particle position. The first term on the right-hand side of Eq. (C1) represents the drag force caused by the relative motion between particles and flow. Because of the airflow velocities (less than 0.5 m/s) and small particle sizes involved in this work, such particle motion can be described by Stokes' law (Lai & Chen, 2006) . Then the particle relaxation time can be expressedwhere C c is the Cunningham correction factor. The second term represents the gravitational force. The third term represents the Brownian force. The Brownian force, very important for submicron particles, is modeled as a Gaussian white noise random process. The simulation procedure for Brownian excitation has been described in detail by Li & Ahmadi (1992) . It should be mentioned that due to the relatively small effect on indoor aerosol particles (Lai & Chen, 2006) , the lift force is not taken into consideration. The velocities and trajectories of particles can be obtained by solving Eqs. (C1) and (C2). For Eq. (C1), the Euler implicit discretization scheme is used as it is unconditional stable for all particle relaxation times while the trapezoidal discretization scheme is applied to solve Eq. (C2). Since the instantaneous particle locations may not coincide with air-phase grid points, the second-order Lagrangian polynomial interpolation scheme could be used, which yields the airflow velocity at particles' positions. Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.jaerosci. 2013.04.004.