An efficient lattice Boltzmann model for indoor airflow and particle transport.

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.

Introduction

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 and 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 and 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 and Ahmadi, 2007, Wang et al., 2012, Zhang and 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.

In recent years, the lattice Boltzmann (LB) method has emerged as a promising alternative computational technique for simulating fluid flows with complex physics (Chen and Doolen, 1998, Succi, 2001). Unlike conventional numerical methods, which solve the discrete macroscopic Navier–Stokes (N–S) equations, the LB method is based on the solution of discrete Boltzmann equation on a lattice with discrete velocity fields (Qian et al., 1992). It has several advantages over conventional CFD. As the pressure in the LB method is calculated by using the equation of the state, there is no need to solve the Poisson equation as the conventional CFD methods do, which often costs significant computational efforts. This makes the LB method computationally more efficient and preferred in large-scale simulations (Stratford & Pagonabarraga, 2008). It can deal with an arbitrary complex curved boundary readily; its intrinsic parallel nature also makes it computationally efficient in that high performance computers can be fully utilized.

The potential applicability of the LB method in IAQ has not been explored fully. Only a few related works have been reported, e.g., Przekop et al. (2003) applied 2D lattice Bhatnagar–Gross–Krook (LBGK) model to study the local structure of deposited particles in fibrous filters. Jafari et al. (2010) investigated particle dispersion and deposition over a square cylinder in a 2D channel flow using LBGK. The current authors also employed multiple-relaxation-time LB (MRT-LB) method to study particle dispersion and deposition in a 3D scaled ventilated chamber (Ding et al., 2012). Previous work was mainly focused on low Reynolds number conditions. As indoor airflow is mainly turbulent, there is a great demand for the application of LB for turbulent flows.

Large eddy simulation (LES) is often a good choice for capturing the turbulence as computational costs are significantly lower than direct numerical simulations (DNS) while its resolution is higher as compared with Reynolds-averaged Navier–Stokes (RANS) equations. Furthermore, LES directly resolves the large eddies mainly attributed to the transport of momentum and energy of turbulence, and only the subgrid-scale eddies are modeled. This can result in transient airflow conditions which play important roles in particle transport. Conventionally, the IAQ problems are almost resolved by RANS approach, and recently the LES based on the N–S equations has become more popular (Tian et al., 2006).

LES can also be performed using the LB method. The simplest way is to incorporate the LB method with the Smagorinsky model. The computational procedure of the Smagorinsky model in LB is entirely local, which makes it compatible with the parallel nature of the LB method. The pioneering work of Hou et al. (1996) paved the way to implement the Smagorinsky model in LBGK. They applied it to study dynamics and the Reynolds number dependence of flow structures of a 2D-lid-driven cavity flow. Fernandino et al. (2009) performed LES of 3D turbulent open duct flow using LBGK. In addition, it has been well demonstrated that significant improvements in numerical stability can be achieved by using MRT-LB (Lallemand & Luo, 2000), which directly results in a drastic gain in computational efficiency. It is therefore natural to incorporate MRT-LB with the Smagorinsky model to improve turbulence simulations using the LB method. Krafczyk et al. (2003) introduced the Smagorinsky model into MRT-LB and carried out simulations of the flow over a surface-mounted cube in a channel. Yu et al. (2006) applied this method to simulate turbulent square jet flow. Based on these previous studies, there is a potential to extend the application of LES using the LB method, especially MRT-LB, to indoor environments.

One basic criterion about the application of MRT-LB in practical indoor environments is to make the number of grids as small as possible. In wall-bounded turbulent flows, lengths are very small when the flow is close to a wall, requiring relatively fine grids to resolve turbulent structures adequately. When it is farther from the wall, coarser grids can be introduced as the turbulence length scale is larger. Thus, there is no need to employ a uniform grid fine enough to resolve the wall regions through the domain which would entail significant computational cost. Moreover, the scales of turbulent flow vary locally in general situations (Premnath et al., 2009). Therefore, the local grid refinement strategy in conjunction with the LB method (Yu et al., 2003) can be an attractive choice. Filippova and Hänel (1998) proposed a local refinement method in which a coarse grid covered the whole domain, and a finer grid in critical sub-regions was superimposed on the basic coarser grid. This allowed two-way information exchange at the grid interface at the post-collision stage after rescaling the density distribution functions on both grid levels. Based on Filippova and Hänel (1998), Lin and Lai (2000) proposed a method which allowed one-way information exchange at the post-streaming stage without rescaling the density distribution function. Yu et al. (2002) introduced the multi-block grid to the LB method where the computational domain was divided into different blocks which did not overlap each other. The grid size was uniform with the desired resolution in each block. The blocks were connected only through interface information exchange which followed the scheme of Filippova and Hänel (1998). Performances of these methods, including accuracy and efficiency aspects with emphasis on the conservative properties across the grid interface, were evaluated systematically by Yu et al. (2003).

In this paper, LES computation of turbulent airflow in indoor environments was carried out using the three-dimensional MRT-LB method with the Smagorinsky model. Experimental and numerical data for a model room (Posner et al., 2003, Tian et al., 2006) were employed to verify the present simulation. Moreover, we incorporated variable resolutions in the MRT-LB method by introducing a local grid refinement approach within the multi-block grid (Yu et al., 2002). Computational accuracy and efficiency of the MRT-LB method with grid refinement were evaluated by comparing with the corresponding results of the uniform grid case. Then the particle dispersion in the model room was simulated based on the airflow conditions obtained above.

Flow simulation

Multiple-relaxation-time lattice Boltzmann method with the Smagorinsky model

In this section, the MRT-LB method with 19 discrete velocities in three dimensions (i.e. D3Q19 model), employed in our recent work (Ding et al., 2012), were also used for flow simulation. The evolution equation for density distribution function f (x, t) can be expressed aswhere subscripts represent the direction of discrete velocities. e is the discrete velocity vector for the D3Q19 model, m (x, t) are the moments of f (x, t). M is the 19-by-19 transformation matrix and is the inverse of M . is the collision matrix in moment space. It is a diagonal matrix with the formwhere s , for i=0, 1, …, 18, specifies the relaxation frequencies for all moments toward their corresponding equilibria . 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 is introduced to model the turbulence. Hence, the total viscosity can be expressed as the sum of and

The turbulent eddy viscosity can be expressed aswhere C denotes the Smagorinsky constant and is set to be 0.16, as suggested by Krafczyk et al. (2003); represents the lattice spacing, ( in Cartesian coordinates) is the strain rate tensor, i.e., and. In the MRT-LB method, can be computed directly from non-equilibrium moments (). 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).

Multi-block grid refinement

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 , i.e., (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 =Re . The convective refinement is chosen which means the characteristic velocities on the coarse and fine block are equal. Thus the relaxation times and must follow the following relationship:

Fig. 1

Schematic of grid refinement.

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 iswhere subscripts . 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.

Computational domain and boundary conditions

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 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 °C. The air density was 1.18 kg/m3 and the kinematic viscosity was 1.46×10−5  m2/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.

Fig. 2

Schematic of the computational domain.

Particle dynamic solutions

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/m3 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.

Results and discussion

Flow simulation results

Grid independence test

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 near-wall 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.

Table 1

Computational parameters.

No.	Reynolds number (Re)	Inlet velocity (uin, m/s)	Meshes in each direction	Grid points in the partition
MRT-LB-LES1	1628	0.235	271×91×136	4×46×137
MRT-LB-LES2			362×121×181	5×61×182
MRT-LB-LES3			453×151×226	6×76×227
MRT-LB-LES4			543×181×271	7×92×272

Fig. 3

Time-averaged velocity vectors and their magnitude distributions in the mid-plane of the model room (z=0.2285 m) ((a) MRT-LB-LES1; (b) MRT-LB-LES2; (c) MRT-LB-LES3 and (d) MRT-LB-LES4).

Comparison with experimental and other numerical results

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 LES as well as the MRT-LB method. The measured peak velocity is 0.249 m/s, while the computational peak values are 0.244 m/s, 0.246 m/s, 0.242 m/s and 0.242 m/s by the standard k–ε model, the RNG k–ε model, the RNG LES and the MRT-LB method respectively. There is a difference of only 0.007 m/s between the MRT-LB method and the experimental results. The axial velocity then decreases as the airflow approaches the floor and the computational velocities drop off faster than the experimental values. The largest discrepancy occurs between the standard k–ε model and the measurements. However, it can be seen that compared with the RNG k–ε model, significant improvement has been made when the RNG LES and MRT-LB methods are employed.

Fig. 4

Comparison of the time-averaged vertical velocity along the vertical inlet jet axis.

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.

Fig. 5

Comparison of the time-averaged vertical velocity along the horizontal line at mid-partition height.

Fig. 6

Comparison of computational time-averaged velocity field in the mid-plane of the model room. (Results of the standard and RNG k–ε model as well as the RNG LES are reproduced with permission from Fig. 6 of Tian et al. (2006).)

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.

Fig. 7

Comparison of computational RMS velocity distribution in the mid-plane of the model room. (Results of the standard and RNG k–ε model as well as the RNG LES are reproduced with permission from Fig. 5 of Tian et al. (2006).)

Computation with multi-block grid refinement

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 =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 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×107 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×107, 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.

Fig. 8

Schematic of grid arrangement in the mid-plane of the model room.

Fig. 9

Time-averaged velocity distributions in the mid-plane of the model room (z=0.2285 m) (MRT-LB-LES-MBGR, dashed line: interface between the fine and coarse grids).

Fig. 10

Computational RMS velocity distribution in the mid-plane of the model room (MRT-LB-LES-MBGR, dashed line: interface between the fine and coarse grids).

Table 2

Time cost of computation.

No.	Total number of elements	Average time cost per iteration (s)
MRT-LB-LES1	3,353,896	0.17
MRT-LB-LES2	7,928,162	0.35
MRT-LB-LES3	15,459,078	0.63
MRT-LB-LES4	26,634,693	1.03
MRT-LB-LES-MBGR	29,316,002	1.28

Particle dispersion in the model room

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 (Seepana & Lai, 2012). 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.

Fig. 11

Comparison of the number of tracked particles in the model room (d=1 μm).

Fig. 12

Comparison of computational particle concentration in the mid-plane of the model room at 130 s (d=1 μm). (Results of the standard and RNG k–ε model as well as the RNG LES are reproduced with permission from Fig. 9 of Tian et al. (2006).)

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.

Fig. 13

Comparison of the number of tracked particles in the model room (d=10 μm).

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.

Conclusions

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.