Simulation of indoor harmful gas dispersion and airflow using three-dimensional lattice Boltzmann method based large-eddy simulation

The lattice Boltzmann method (LBM) and large-eddy simulation (LES) are combined with a scalar subgrid-scale model to simulate the indoor air velocity field and harmful gas dispersion. The LBM-LES model is validated by comparing its results with published experimental and numerical simulation results. Taking a simplified chemical building as the scenario, the relative ventilation efficiency is evaluated based on the maximum harmful gas concentration, and configurations with centralized and distributed harmful gas sources with both mixing ventilation (MV) and displacement ventilation (DV) systems are considered. According to the results, if the density of the harmful gas is less than the air density, the DV system is more efficient than the MV system. The DV system is more stable than the MV system under fluctuating relative ventilation efficiency due to changes in the distance between the ventilation vents and in the distance between the centralized gas sources and the exhaust air vent. © 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0045120


I. INTRODUCTION
Indoor harmful gases seriously threaten human safety and health. Harmful gases emit from the indoor facilities into the air such as formaldehyde from the furniture in residential buildings and methane from chemical products in chemical buildings, and methane may even cause an explosion under certain conditions. Designing an effective indoor ventilation system is critical to protecting people from harmful gases and gas explosion. A numerical study of the indoor harmful gas dispersion and airflow can guide the design of the indoor ventilation system.
Nowadays, computational fluid dynamics (CFD) approaches have become powerful tools for simulating the airflow and designing the appropriate ventilation system for indoor environments.
Many researchers have studied indoor ventilation systems, harmful gas behavior, and air quality based on CFD methods. For instance, Chang 1 used a CFD technique to evaluate the effect of traffic pollution on the indoor air quality of a naturally ventilated building with various ventilation control strategies. Su et al. 2 tested three large eddy simulation (LES) models [Smagorinsky subgrid-scale (SGS) model, dynamic SGS model, and stimulated SGS model] to simulate the indoor airflow with different ventilation systems and found that the dynamic and stimulated models perform slightly better than the Smagorinsky model. Steeman et al. 3 combined CFD with an effective penetration depth model to simulate the indoor air and hydrothermal wall interaction; the CFD results agreed well with those of the traditional methods. Min and Xu 4 studied the structure of displacement ventilation (DV) convection with the Reynolds averaged Navier-Stokes (RANS) equation and k-epsilon turbulence model. They also investigated the effect of the Grashof number (Gr) and squared Reynolds number Re ratio on the indoor air convection and ventilation efficiency. D'Agostino et al. 5 used CFD tools to model the effect of ventilation strategies on the microclimate of the "Crypt of Lecce" Cathedral (South Italy) and validated the CFD results with the experimental microclimate data of the crypt. Besides, Zhuang et al. 6 built a CFD model to study the effect of the furniture layout on the indoor air quality with typical office ventilation. They discovered that adjusting the furniture layout significantly influenced indoor air quality. Yang et al. 7 simulated the wind velocity, temperature, and air age fields in a bedroom with wall-mounted air conditioners. Kong et al. 8 adopted a commercial CFD code (STAR-CCM) to simulate the airflow field in an office cubicle. Besides, they evaluated the environment quality in an office cubicle with and without personalized ventilation (PV) based on their model. Wang et al. 9 proposed the index "normalized concentration in the target zone (NC-TZ)" to evaluate the control of a local ventilation system on the indoor air quality and assessed its capability based on the CFD method. Pesi et al. 10 employed the Fire Dynamics Simulator and LES model to investigate the effects of wind flow on plume dispersion generated by indoor fires. Kempe and Hantsch 11 presented a numerical scheme for the LES model of indoor airflow and heat transfer and validated their method with numerical and experimental reference data. Yuce and Pulat 12 conducted the forced, natural, and mixed convection benchmark studies for indoor thermal environments using ANSYS-Fluent. The CFD studies of indoor environments have been mainly concentrated on the airflow in city buildings, and CFD simulations of the harmful gas behavior in chemical buildings are rare.
In recent years, the lattice Boltzmann method (LBM) has been used in the simulations of many sophisticated physical phenomena, in particular, multi-component flow 13 and reactive flow. 14 The LBM has several advantages compared with the traditional CFD methods; 15 these include the ability to handle complex geometries and capture non-equilibrium effects. 16 Furthermore, the LBM can be easily implemented in large-scale parallel computing environments such as graphics processing units (GPUs) 17 because it is an explicit numerical method and comprises of only local operations. 18 The computational speed of the GPU architecture is many orders of magnitude higher than that of the traditional CPU simulation.
Some researchers have investigated indoor airflow with the LBM. For instance, Hasan et al. 19 conducted the LBM simulation of airflow and heat transfer in a model ward of a hospital with no turbulence model. Bin et al. 20 modeled the performance of air filters for cleanrooms using the LBM with no turbulence model. Zhang and Lin 21 applied the Bhatnagar-Gross-Krook (BGK) model with no turbulence model to simulate indoor airflow. They reported that the LBM results, experimental data, and CFD data were in good agreement. The high Reynolds number of indoor airflow makes it a typical turbulent system. A larger Reynolds number usually requires much more lattices to get a certain value of relaxation time, which plays an important role in the rate of convergence and stability. However, the simulation time takes a toll if the length of a lattice cell is too small, 21 and it requires tremendous computing power to obtain the acceptable simulation time. In fact, the LBM results are equivalent to direct numerical simulation (DNS) results if no turbulence model is included, 22 but it is also expensive for engineering applications. Some researchers introduced turbulence models from LES 23 and RANS 24 simulation into the LBM to solve engineering problems. The RANS simulation is widely used in engineering problems because of its low economic cost. However, it cannot provide detailed fluid velocity fluctuations because of its time-averaged character. 25 The LES model requires fewer computations than the DNS 26 and provides detailed fluid velocity fluctuations. Therefore, many researchers have adopted the LBM-LES model to solve engineering problems such as the gust index in an urban area, 27 gust structure in the atmospheric boundary layer, 28 and wind comfort in a full-scale city area; 29 the model has also been employed in nuclear engineering studies. 30 Khan et al. 31 incorporated the subgrid effect into the LBM to simulate indoor airflow and temperature and validated the model results with a traditional CFD method. Sajjadi et al. 32,33 coupled the LBM and LES models to simulate indoor airflow and particle behavior; in their studies, the LBM-LES results and experimental data agreed well. Mengtao et al. 34  The indoor harmful gas dispersion and airflow have two characters: turbulent and multi-component. Therefore, to simulate the indoor harmful gas dispersion and airflow, the LBM model must have two abilities: the ability to simulate the turbulent flow and describe the multi-component flow. The simulation of the turbulent flow can be realized by the LBM-LES model. To describe the multicomponent flow, there are mainly four kinds of multi-component models based on the LBM framework: the scalar model, 40,41 the single-collision model, 42,43 the split-collision model, 44,45 and the forcing-term model. 46 In recent years, based on the framework of above models, some researchers presented new models to simulate the multi-component flow. For instance, Hosseini et al. 47 presented a mass-conserving advection-diffusion lattice Boltzmann model for multi-species reacting flows. Chai et al. 13 established an MRT LBM model for diffusion in multi-component mixtures. Lin et al. 48 proposed an MRT discrete Boltzmann model to simulate the multi-component mixture with nonequilibrium effects. Among above models, the scalar model is the most simple and stable one, and it requires the least computing resources than other models. Considering that the large-scale simulation needs huge computing resources and the high Re turbulent flow requires high numerical stability, the LBM-LES model is combined with the scalar model to simulate the indoor harmful gas dispersion and airflow in this research.
The aim of this study is to validate the applicability of the LBM-LES model combined with the scalar model in simulating the indoor airflow and harmful gas behavior first and then, taking ARTICLE scitation.org/journal/adv a simplified chemical building as the scenario, evaluate the ventilation efficiency of the MV and DV systems in removing the indoor harmful gas. Section II introduces the theories of the LBM model, LES model, and scalar model and shows how these models combine. Section III introduces the validation of the presented model and the mesh sensitivity analysis. Section IV introduces the test cases and the simulation results. Section V gives the main conclusions of this study.

II. METHODOLOGY
In this study, the lattice BGK (LBGK) model 49,50 is coupled with the LES 51 and scalar model to simulate the indoor airflow velocity fields and harmful gas dispersion.

A. LBGK model
The evolution equation of the LBGK model can be written, as in Ref. 50, as where f α is the velocity distribution along the α direction, ⃗ r is the space position, Δt is the lattice time interval (equal to 1), τα is the relaxation time along the α direction (equal to ν l0 /c 2 s f Δt + 1/2 in normal LBGK models), and ν l0 is the air molecular viscosity in lattice units.
The eα in Eq. (1) represents the discrete fluid particle velocity vector. The D3Q15 model is selected as the discrete fluid particle velocity vector to simulate the airflow velocity field (Fig. 1). It can be expressed, as in Ref. 52, as where c is the lattice speed (equal to Δx/Δt) with Δx being the lattice space step (equal to 1). The f eq α in Eq. (1) represents the equilibrium velocity distribution along the α direction and can be expressed as follows: where ⃗ u is the macro-particle velocity, ρ is the macro-density, c sf is the lattice sound speed (equal to c/ √ 3), and ωα represents the weight coefficient in the α direction with ω 0 = 2/9, ω 1-6 = 1/9, and ω 7-14 = 1/72.
The relationship between ρ, ⃗ u, and the velocity distribution f α can be expressed as follows:

B. LES model
The LES Smagorinsky model is applied in this study. In the Smagorinsky model, 53 the turbulence viscosity νt is related to the strain rate Sij = 1 2 (∂iuj + ∂jui) and filter length scale Δx, where CS is the Smagorinsky constant. According to Ref. 51, the turbulence viscosity νt can be written as follows: where τ 0 is the laminar relation time (equal to v l0 /c 2 s f Δt + 1/2) and v l0 is the lattice laminar viscosity. In this study, CS and Δx are set to 0.16 and 1, respectively; Q ij can be calculated as follows: The turbulence relaxation time τt equals 3νt. The relationship of τα in Eq. (1), τ 0 , and τt can be expressed as τα = τ 0 + τt.
An explicit power-law-based wall model is used in the LBM-LES model. The details about the explicit power-law-based wall model could be found in Ref. 54.

C. Scalar model
By assuming that the harmful gases distribute slowly inside buildings, the effect of gas emission on the airflow velocity can be ignored. The scalar model is applied to simulate harmful gas dispersion. Therefore, the concentration distribution g must be introduced. The evolution equation of g can be written, as in Ref. 41, as where τ β is the relaxation time of the harmful gas concentration distribution along the β direction (equal to D/c 2 sg Δt + 1/2), c 2 sg equals 1/4 for the D3Q7 space, and D is the diffusion coefficient of the harmful gas in air. In Eq. (8), the e β represents the discrete harmful gas-particle velocity vector. Furthermore, the D3Q7 model is selected to simulate the harmful gas behavior (Fig. 2), which can be expressed as In Eq. (8), the g eq β represents the equilibrium concentration distribution along the β direction and can be expressed as follows: where Cg is the concentration of the harmful gas and ω β represents the weight coefficient in the β direction (ω 0 = 1/4, ω 1-6 = 1/8).
In Eq. (8), the source term F β can be written as 32,33 where F b is the buoyancy of the harmful gas. The relationship between the concentration distribution g and the concentration of harmful gas Cg can be written as follows:

D. Unit conversion
The principle of unit conversion between the lattice units and the physical units is to keep the Reynolds number Re equal in both the lattice system and the physical system. Assuming that the characteristic length in lattice units is L l0 , the characteristic velocity in lattice units is U l0 , the characteristic length in physical units is L p0 , the characteristic velocity in physical units is U p0 , and the air molecular viscosity in physical units is ν p0 , the Re in both systems can be written as

III. VALIDATION
To validate the LBM-LES model, two rooms that have been experimentally and numerically studied by previous researchers are  velocity profile is used at the room inlet and a zero gradient condition for the air velocity is used at the room outlet. In this study, a computational mesh of 90 × 30 × 30 is adopted in the LBM simulation, which is slightly coarser than the 130 × 34 × 34 mesh used by Su et al., and the boundary conditions of the LBM model at the inlet and outlet of room 2 are the same as those in room 1. Figure 7 shows the mid-plane velocity field of room 2. Based on a comparison of the results and those in Ref. 2, it is evident that the LBM-LES model and the previous pure LES provide similar velocity fields. The time-averaged X-velocity at the mid-plane along the Z-direction at the lines X = H and X = 2H is compared with the experimental and Dynamic Sub-Grid Scale (DSGS) LES results in Fig. 8. It can be seen that the results of the LBM-LES model show good agreement with the experimental and DSGS model results.
The characteristic of the LES model compared with the timeaveraged model is its ability to capture the transient velocity fluctuation. Figure 9 shows the transient X-velocity at the position whose

IV. TEST CASES
Mixing ventilation (MV) and displacement ventilation (DV) systems are two typical indoor ventilation systems. Many researchers have claimed that DV systems generally have better ventilation efficiency than MV systems. 58,59 However, other researchers have reported opposing results for certain scenarios. 60,61 Furthermore, the relative positions of the supply air vent, exhaust air vent, and contamination sources can significantly affect the ventilation efficiency. Therefore, buildings with centralized and distributed harmful gas sources with both mixing ventilation (MV) and displacement ventilation (DV) systems are considered herein.

A. Test scenarios
The test scenarios are shown in Fig. 10. Methane is a typical chemical industrial material, which explodes when its concentration in the air is 5%-16% and is exposed to an ignition source. 62 For convenience, methane is selected as the target harmful gas in this research. Two typical methane tank layouts (centralized and distributed) under two typical ventilation systems (MV and DV) are considered. Figure 10(a) shows the geometries and dimensions of the building. It is simplified to a cube with dimensions of 9 × 4.5 × 3 m 3 , and the ceiling and sidewall of the room each have a ventilation vent. The dimensions of the ventilation vents are 1 × 1 m 2 , and the methane tanks are cylindrical, each with a height of 1 m and a diameter of 0.6 m. Besides, the centralized methane tanks are simplified to a single tank, as shown in Fig. 10(a). The positions of the distributed methane tanks are shown in Fig. 10(c). The distance between the two ventilation vents is defined as L 1 , and the distance between the centralized methane tank and the sidewall ventilation vent is L 2 . The ventilation vents indicated by blue arrows in Fig. 10   the supply air vents; the vents indicated by red arrows are the exhaust air vents. Methane enters the room from the upper face of the methane tanks. The supply air rate was set to 0.005 m 3 /s, and the methane emission rate of each tank was 0.00005 m 3 /s. The mesh sensitivity analysis was conducted based on the study scenario in Fig. 10(a) for L 1 = 7.5 m and L 2 = 4.5 m. The two mesh schemes are listed in Table I.
To evaluate the ventilation efficiency of the ventilation system, the relative ventilation efficiency based on the maximum methane concentration Emax was introduced based on 63 where Ce is the methane concentration near the exhaust air vent, Csp is the methane concentration near the supply air vent, and Cmax is the maximum methane concentration in the room. The Ce and Csp can be defined as follows: where node i is the internal node that is adjacent to the air outlet, node j is the internal node that is adjacent to the air inlet, m and n are the total numbers of nodes i and j, and Ci and Cj are the methane concentrations at nodes i and j. Figure 11 shows a comparison of the velocity and methane concentration along the horizontal line at 2/3 room height for the two mesh schemes, which both provide similar results with minor differences. Considering the statistical and fluctuating character of the turbulent flow, the slight differences are reasonable and are negligible for engineering applications. Because mesh scheme 1 requires fewer computations than scheme 2, scheme 1 is chosen for further investigation. Figure 12 shows the mid-plane velocity and methane concentration fields in different study scenarios. The different ventilation systems lead to different air velocity fields and methane concentration fields. The minimal methane concentration occurs near the supply air vents in all study scenarios. The methane buoyancy makes the methane concentration near the roof of the buildings higher than that in other positions. For MV systems with both centralized and distributed methane tanks, methane accumulates at the top corner near the air exhaust vent, as shown in Figs. 12 Figure 12 shows the characteristic of the velocity and methane concentration fields of different test scenarios. However, it cannot give the comparison of different study scenarios in ventilation efficiency because the relative positions of the supply air vents, exhaust air vents, and methane tanks significantly affect the airflow field and methane concentration field, even for equal ventilation systems. The effect of the relative position of the supply air vents, exhaust air vents, and methane tanks is discussed below.

B. Results and discussion
To generalize the discussion and compare the ventilation system efficiency quantitatively, the relative ventilation efficiency Emax of different ventilation systems and L 1 and L 2 are calculated based on the simulation results. The definition of the relative ventilation efficiency Emax is shown in Eq. (12). Figure 13 shows a comparison of the relative ventilation efficiencies of the scenarios with centralized methane tanks under different ventilation systems and different L 2 (L 1 = 7.5 m). With increasing L 2 , Emax of the DV system increases; however, Emax of the MV system decreases. It should be noted that the supply and exhaust air vents of DV and MV systems have opposite relative positions. Thus, the decreasing Emax with increasing L 2 is due to the decreasing distance between the centralized methane tank and the exhaust air vent. It implies that the distance between the centralized methane tank and the exhaust air vent should be as small as possible to increase the ventilation efficiency. Its reason is that the decreasing distance between the centralized methane tank and the exhaust air vent can shorten the path of the methane flow. If the centralized methane tank is placed at the position close to the supply air vent, the emitted methane needs to flow through the entire path from the supply air vent to the exhaust air vent, and it will increase the risk of gas accumulation. Besides, it can also be seen from Fig. 13 that the DV system is better than the MV system in most cases, which agrees well with the conclusions presented in previous studies, 58,59 that the fluctuations in the Emax values of the DV system are within 5%, whereas the fluctuations in the Emax value of the MV system can reach 15%, and that the DV system is more stable than the MV system regarding the fluctuations of the relative ventilation efficiency caused by location changes of the methane sources.
Based on the previously presented analysis, a smaller distance between the centralized methane tank and the exhaust air vent results in improved ventilation. However, the exact effect of the AIP Advances distance between the supply and exhaust air vents on ventilation efficiency is still unknown. The distance between the centralized methane tank and the exhaust air vent should be fixed in both ventilation systems when exploring the effect of the distance between the supply and exhaust air vents on ventilation efficiency. In the DV system, the distance between the centralized methane tank and the exhaust air vent remains 2 m for L 1 = L 2 according to the dimensions in Fig. 10. In the MV system, L 2 is set to 2 m to maintain the distance between the centralized methane tank and the exhaust air vent equal to that of the DV system. Figure 14 shows a comparison of the relative ventilation efficiencies of scenarios with centralized methane tanks under different ventilation systems and different L 1 when the centralized methane tank is positioned 2 m away from the exhaust air vent. With increasing L 1 , Emax of both ventilation systems increases. The DV system exhibits a higher Emax and better stability than the MV system. For buildings with centralized methane sources, considering the ventilation efficiency only, the methane sources should be placed as close as possible to the exhaust air vent and the distance between the supply and exhaust air vents should be as large as possible. Besides, the DV system is always better than the MV system with different L 1 . Figure 15 shows a comparison of the relative ventilation efficiencies of the scenarios with differently distributed methane tanks under different ventilation systems and different L 1 . Accordingly, the DV system is always better than the MV system in these scenarios and also exhibits better stability than the MV system with different L 1 , and Emax of both ventilation systems increases with increasing L 1 .
It can be seen from Figs. 13-15 that, for the scenarios with both centralized methane sources and distributed methane sources, the DV system is always better than the MV system, Emax of both ventilation systems increases with increasing L 1 , and the DV system exhibits better stability than the MV system with different L 1 and L 2 . According to the mass conservation, the air velocity near the exhaust air vent of both the MV and DV systems should be equal, and this can also be seen in the velocity fields shown in Fig. 12. However, for the DV system, the exhaust air vent is located at the top of the buildings and the methane buoyancy makes the methane concentration near the building top higher than that at other positions, which causes more methane to be removed from the building even with the same air velocity near the exhaust air vent. For the MV system, the methane buoyancy makes the methane concentration near the exhaust air vent lower than that at the higher positions, which results in less methane being removed from the building than the DV system. Therefore, it can be concluded that the buoyancy of methane in the air is the main reason that the DV system is always better than the MV system in these scenarios.
It can be seen from Fig. 10 that L 1 represents the distance between the supply air vent and the exhaust air vent. The distance between the vents can affect the area that fresh air can reach in the building. If the distance between the vents is too small, the air short circuit will happen, which means most of the fresh air from the supply air vent flows out of the building through the exhaust air vent quickly and the fresh air reaches only a small area between the vents. Methane tends to accumulate in areas where the fresh air cannot reach. Therefore, the reason that Emax of both ventilation systems increases with increasing L 1 is that the increasing L 1 can expand the area where fresh air can reach. The reason why the DV system exhibits better stability than the MV system with different L 1 mainly lies in the building geometry and the methane buoyancy. The flow of inlet air can be regarded as a jet, which means there is still some fresh air that can get rid of the air short circuit and reach other positions of the building, counteracting some of the effects of the air short circuit. Figure 12 shows that the jet of the DV system is limited by the two sidewalls along the length direction and that of the MV system is limited by the building roof and floor. Figure 10 shows that the length of the building is 9 m and is three times as much as the height of the building. The jet direction of the DV system is along the length direction, and that of the MV system is along the height direction. It means that there is more fresh air to get rid of the air short circuit in DV systems than MV systems. Besides, the methane buoyancy can also counteract some of the effects of the air short circuit in the DV system. It is the reason why the DV system exhibits better stability than the MV system with different L 1 . The reason why the DV system exhibits better stability than the MV system with different L 2 also mainly lies in the methane buoyancy.
In summary, the ventilation efficiency of DV systems is higher than that of MV systems in most cases and the DV system exhibits better stability than the MV system under fluctuating relative ventilation efficiency due to changes in L 1 and L 2 if the harmful gas is methane. The reason lies in the methane buoyancy. When designing the ventilation systems of chemical buildings with harmful gas sources, the comparison between the density of the harmful gas and the air density should be conducted first; if the density of the harmful gas is smaller than the air density, the DV system should be selected; otherwise, more studies should be conducted. For buildings with distributed gas sources or centralized gas sources, the distance between the supply and exhaust air vents should be as large as possible to avoid the air short circuit. For buildings with centralized gas sources, the centralized gas sources should be placed as close to the air exhaust vent as possible to shorten the flow path of the harmful gas.

V. CONCLUSION
An LBM-LES modeling framework for turbulent flow simulation has been extended to indoor airflow and pollutant gas dispersion in buildings. The results of the LBM-LES model match those of published experimental and numerical studies and can correctly simulate the indoor airflow and harmful gas behavior. The presented model can capture the transient airflow velocity fluctuations.
If the density of the harmful gas is smaller than the air density, the ventilation efficiency of DV systems is higher than that of MV systems in most cases and the DV system exhibits better stability than the MV system under fluctuating relative ventilation efficiency due to changes in the distance between the supply air vent and the exhaust air vent and in the distance between the centralized gas sources and the exhaust air vent. When designing the ventilation systems of the buildings with harmful gas sources, the comparison between the density of the harmful gas and the air density should be made first.
For buildings with distributed gas sources or centralized gas sources and with MV systems or DV systems, the distance between the supply and exhaust air vents should be as large as possible to avoid the air short circuit. For buildings with centralized gas sources and with MV systems or DV systems, the centralized gas sources should be placed as close to the air exhaust vent as possible to shorten the flow path of the harmful gas.

ACKNOWLEDGMENTS
The support from the UK Engineering and Physical Sciences Research Council under the project "UK Consortium on Mesoscale Engineering Sciences (UKCOMES)" (Grant No. EP/R029598/1) is acknowledged. This research was also supported by the National Natural Science Foundation of China (Grant Nos. 51274206 and 51404277).
No conflict of interest exits in the submission of this manuscript, and the manuscript is approved by all the authors for publication. I would like to declare on behalf of my co-authors that the work described was original research that has not been published previously, and not under consideration for publication elsewhere, in whole or in part.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.