CHEMICAL ENGINEERING TRANSACTIONS VOL. 57, 2017 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš, Laura Piazza, Serafim Bakalis Copyright © 2017, AIDIC Servizi S.r.l. ISBN 978-88-95608- 48-8; ISSN 2283-9216 Validation of a Numerical Model of Smoke Propagation in a Semi-Enclosed Space Aryslan Malik, Luis R. Rojas-Solórzano Dept. of Mechanical Engineering, School of Engineering, Nazarbayev University, Kazakhstan aryslan.malik@nu.edu.kz This paper aims to develop appropriate boundary conditions to be used in the numerical analysis (CFD) of smoke transport within a semi-enclosed space resembling a typical road tunnel. The study considers the following stages: (a) Model set up: the geometry of an existing tunnel, with available experimental data to validate numerical results, is recreated in 3D CAD model with physical properties and known boundary conditions; (b) Mesh verification: different aspects and techniques of mesh generation, such as local and global refinement and first/second-order discretization schemes were considered for the mesh sensitivity analysis. The aim of this stage is to demonstrate that numerical results are mesh-independent, i.e. use as large as possible element sizes to reduce the computational cost, yet refine grid in crucial locations to record important features of transport and thermal phenomena; (c) Model validation: This stage considers the inaccuracies associated to turbulence models, intensity of turbulence at inlet boundaries, buoyancy, turbulent diffusion and convection on the numerical results of the study. Different combinations of numerical aspects of the study were tested. Results of these combinations were compared with experimental data and the most accurate combination of variables was taken as benchmark. Obtained numerical model serves as a guideline tool capable of reproducing propagation of smoke in semi-enclosed spaces within minimal uncertainty than can be used in further studies as a predicting tool for engineering purposes. 1. Introduction In current era, computational techniques are burgeoning and arising as a complementary tool or under certain circumstances one of the main analysis tool used in investigating ventilation systems in cases of fire (Blanchard et al., 2012; Lin and Chuah, 2008; Lee and Ryou, 2006). Furthermore, due to high experimental cost of fire simulations, CFD (Computational Fluid Dynamics) is becoming more beneficial by virtue of its low cost and simplicity of implementation. This trend was developing on a par with related software and hardware in such a way that CFD analysis has become accessible as everyday tool for preliminary assessment and companion to final design stages with few key experiments for validation purposes. Particularly, fire safety analysis in road or underground tunnels is especially crucial because it is concerned with safety of human lives. This fact underscores and justifies the need for a detailed study of variables affecting the numerical simulation of these regretfully common hazards (Carvel and Beard, 2005), (Lacroix, 2001). An appropriate CFD analysis is supported on underlying principles and aspects of computational methods which must be checked and verified so that the results of such CFD simulations are accurate enough to predict actual flow rates, temperatures, velocities and other relevant information of the flow field within an acceptable level of error. In this paper computational investigation was carried out using the finite volume method supported on the CFD platform ANSYS CFX. Different turbulence models and meshing techniques were utilized and compared with experimental data and boundary conditions related to the ingestion of air at tunnel exits were scrutinized. The paper starts with a set of governing equations followed by computational approach and finalized with discussion and conclusions. 2. Governing Equations The numerical CFD analysis is performed by solving a set of governing equations which include fluid mass, DOI: 10.3303/CET1757246 Please cite this article as: Malik A., Rojas-Solorzano L., 2017, Validation of a numerical model of smoke propagation in a semi-enclosed space, Chemical Engineering Transactions, 57, 1471-1476 DOI: 10.3303/CET1757246 1471 momentum, energy and species conservation under turbulent regime. Therefore, essentially, Reynolds Averaged Navier-Stokes equations (RANS) and relevant accompanying Reynolds-averaged energy, species and two-equation turbulence model comprises the system of differential equations to be solved. Nowadays, several turbulence models are being used in order to approximate average effects of turbulence. Among these turbulence models, the k-ε model, which incorporates transport differential equations of the turbulent kinetic energy and dissipation, has proven to be effective in solving problems related to closed conducts (Calautit and Hughes, 2014; Rohdin and Moshfegh, 2007; Calautit et al., 2014; Hussain and Oosthuizen, 2012). Thus, the set of governing equations to be solved are presented as follows:  Continuity: ∇ ∙ U = 0 (1)  Momentum: 𝐷(𝜌𝑈) 𝑑𝑡 = 𝛻𝑝 + 𝜇𝛻2𝑈 + 𝜌𝛽(𝑇 − 𝑇𝑟𝑒𝑓 ) + 𝑆𝑀 (2)  Energy: 𝜕𝜌ℎ 𝜕𝑡 − 𝜕𝜌 𝜕𝑡 + 𝛻 ∙ (𝜌𝑈ℎ) = 𝛻 ∙ (𝜆𝛻𝑇) + 𝑆𝐸 (3)  Turbulent kinetic energy: 𝜕(𝜌𝑘) 𝜕𝑡 + 𝛻 ∙ (𝜌𝑈𝑘) = 𝛻 ∙ [(𝜇 + 𝜇𝑡 𝜎𝑘 ) ∇𝑘] + 𝑃𝑘 − 𝑝𝜀 (4)  Turbulent dissipation: 𝜕(𝜌𝜀) 𝜕𝑡 + 𝛻 ∙ (𝜌𝑈𝜀) = 𝛻 ∙ [(𝜇 + 𝜇𝑡 𝜎𝜀 ) ∇𝜀] + 𝜀 𝑘 (𝐶𝜀1𝑃𝑘 − 𝐶𝜀2𝜌𝜀) (5)  Turbulent viscosity relation: 𝜇𝑡 = 𝐶𝜇𝜌 𝑘2 𝜀 (6) Empirical coefficients correspond to the standard values (Versteeg and Malalasekera, 2007). Fire was simulated as an entry of smoke as a species diffusing in the background air and added heat flux. Despite, the combustion reaction will not be taken into consideration for the sake of decreasing computational cost, previous works using similar approach have proven to be enough representative and accurate when modeling diffusion and turbulent dispersion of smoke gases in closed spaces (Vittori et al., 2008; Hwang and Edwards, 2005; Lin and Chuah, 2008; Lee and Ryou, 2006). Smoke transport is considered by only one of its components, CO, which due to its homogeneity miscibility in air, can be simulated as species transport (Vittori et al., 2008; Kumar and Mazumder, 2010; Lin and Chuah, 2008). The equation for conservation of species concentration (φ) is given as follows:  Species transport: 𝜕𝜑 𝜕𝑡 + 𝛻 ∙ (𝑈𝜑) = 𝛻 ∙ [(𝜌𝐷𝜑 + 𝜇𝑡 𝑆𝑐𝑡 ) 𝛻 ∙ ( 𝜑 𝜌 )] + 𝑆𝜑 (7) The solution for these equations is nearly impossible to be obtained analytically. For practical reasons the approximate solution is found using numerical methods, subject to corresponding boundary and initial conditions. The descriptions of symbols shown in governing equations are listed in Table 1. 3. Computational Approach 3.1 Model set up The objective of this stage is to replicate the geometry of a tunnel which was used in experiments of fire simulation published by Vauquelin and Megret (2002). Authors of the study were concerned with effects of location, orientation and shape of exhaust ducts on the efficiency of the ventilation system, which is described as a ratio of the quantity of smoke extracted via the exhaust to the produced smoke: 𝜀 = �̇�𝐸𝑋𝑇 �̇�𝑆𝑀𝑂𝐾𝐸 (8) 1472 These experiments were conducted on a 1:20 scale tunnel which had original dimensions of 200m length, 10m width and 5m height. Thus, the source of smoke was also scaled accordingly. The CAD geometry in our numerical model was generated using original dimensions. Out of several experimental tests conducted by Vauquelin and Megret, one specific test was chosen for validation of our current model. Table 2 presents the boundary and source data corresponding to that experimental test which was later used for validation of our numerical model. It should be noted that for the purpose of saving computational resources, and given the steady state nature of our flow analysis, the symmetry of the tunnel across its vertical mid-plane was exploited so that half of the domain was cut as shown in Figure 1. Table 1: Nomenclature Symbol Description Units U velocity vector m/s ρ t fluid density time kg/m3 s p pressure Pa g gravity acceleration m/s2 μ μt viscosity turbulent viscosity kg/ms kg/ms μeffec effective viscosity kg/ms δ identity matrix SM SE source term of momentum source term of energy h total specific enthalpy J λ Thermal conductivity W/m-K T temperature K P Turbulence production due to viscous and buoyancy forces Cε1, Cε2, σε, σk, Cμ Constants of the turbulence model ṁEXT ṁSMOKE Mass flow rate of smoke extracted by the ventilation Mass flow rate of the smoke source kg/s kg/s 𝜑 Species concentration ε System efficiency, ratio of extracted smoke to the produced smoke Table 2: Boundary and source conditions Heat intensity Source diameter Smoke flow rate Source temperature Extraction velocity 1 MW 0.98 m 5.8 m3/s 267oC 0.36 m/s Figure 1: (a) Computational geometry of the tunnel in actual size generated in ANSYS CAD tool (b) boundary conditions incorporated into the domain a) b) Extraction Smoke inlet Opening Symmetry plane Extraction (magnified) Smoke inlet (magnified) 1473 3.2 Mesh Verification The purpose of this stage is to make mesh sensitivity analysis. Different meshing techniques were used. These meshing methods were utilized according to the model and foreseen needs of refinement near the walls and the source of fire. As the top of the domain (ceiling) right above the fire source is hit by inflow of the smoke, mesh refinement is applied to the top surface in order to capture the most important transport phenomena gradients. Moreover, mesh was refined close to the extraction and it was coarsened towards the center of the domain to reduce the overall computational cost. In order to prove that numerical results are mesh independent, global mesh parameters such as maximum tetrahedrons’ size and maximum face size were decreased each time the mesh was refined and the control variable (efficiency) was recorded and compared with previous result. Three meshes were generated by varying maximum tetrahedron and face sizes: 1.25m, 0.825m and 0.54m. Mesh refinements were applied to each wall, as well as near extractions. Figure 2: Depiction of the mesh with inflations and refined face sizing Table 3 shows that a 0.825m tetrahedral mesh produces accurate numerical results at relatively small number of nodes, thus requiring less time and computational cost. Therefore, this mesh will be used in further simulations. Table 3: Efficiency of ventilation system, effects of tetrahedron and face sizes Maximum Tetrahedron and Face sizes Number of nodes Number of elements Efficiency Relative Error 1.25m 0.825m 0.54m 113,541 151,504 250,234 105,130 140,467 233,500 0.1117 0.1139 0.1141 N/a 1.9% 0.2% 3.3 Model Validation This stage considers the effects of turbulence model, intensity of turbulence at inlet boundaries and inclusion of buoyancy effects, turbulent diffusion and convection on the numerical results of the study. Different combinations of numerical aspects of the study were tested. The results of these combinations are compared with experimental data. The most accurate combination of variables will be taken as benchmark and therefore, the obtained numerical tool, capable of reproducing propagation of smoke in semi-enclosed spaces within minimal uncertainty, will serve in further studies as a predicting tool for engineering purposes in similar hazard situations. First of all, as it was mentioned before that k-ε turbulence model is used as the standard along the whole process of modeling due to its proven effective characteristics for these type of flows (Calautit and Hughes, 2014; Rohdin and Moshfegh, 2007; Calautit et al., 2014; Hussain and Oosthuizen, 2012). Secondly, the turbulence intensity at inlet, extraction and outlets is varied from 1% to 10%. Preliminary results, not shown in this paper, demonstrate that these intensities have no tangible effects on the overall efficiency of the system. Thus, the turbulence intensity at all boundaries is set to be 5% as a default. At this stage the information regarding boundaries such as inlet and extraction is complete and known in advance from experimental setup. However, the opening boundaries at the extreme of the tunnel are somewhat ambiguous, which means that some conditions such as: whether air was entering the domain either contaminated or not, is completely unknown. In our numerical model, openings conditions, which prescribe uniform atmospheric pressure and fully developed flow, are used as boundary conditions, meaning that it is assumed that air is free to exit or enter the domain, as far as the prescribed mathematical and physical conditions are held. However, certain attributes such as the species concentration of air entering into the domain through these boundaries, if any, was not determined in the original experiment and thus, is unknown. In this case to reduce the scope of uncertainty the simulation is performed twice; during first simulation it is assumed that the air entering through 1474 opening boundaries is fresh, which means that the species concentration of air entering the domain through openings is zero. The result of this first simulation is analyzed and the maximum species concentration exiting at opening boundaries is recorded. Before performing second simulation, the recorded value is entered as an input of species concentration for any volume of air entering at opening boundaries. In order to check the validity of this method, three different case studies changing the power of the fire were assessed. Table 4 shows the recorded data for all three case studies with: 1 MW, 4 MW and 10 MW. After performing first simulation (fresh air entrainment) for each case study, the maximum species concentrations entering at opening boundaries are found to be 0.17, 0.34 and 0.5, respectively. Table 4: Mass flow rates and species concentrations* 1 MW simulation Fresh air at openings Air with species at openings (0.17 fraction) Extraction (kg/s) Inlet (kg/s) Opening 1 (kg/s) -0.85 (species: 0.45) 3.38 (species: 1.0) -1.378 -0.85 (species: 0.53) 3.38 (species: 1.0) -1.376 Opening 2 (kg/s) -1.15 -1.151 4 MW simulation Fresh air at openings Air with species at openings (0.34 fraction) Extraction (kg/s) Inlet (kg/s) Opening 1 (kg/s) -2.87 (species: 0.72) 11.3 (species: 1.0) -4.610 -2.87 (species: 0.80) 11.3 (species: 1.0) -4.605 Opening 2 (kg/s) -3.818 -3.824 10 MW simulation Fresh air at openings Air with species at openings (0.5 fraction) Extraction (kg/s) Inlet (kg/s) Opening 1 (kg/s) -6.73 (species: 0.92) 26.5 (species: 1.0) -10.78 -6.73 (species: 0.95) 26.5 (species: 1.0) -10.78 Opening 2 (kg/s) -8.94 -8.94 (*) Positive and negative values represent air flowing in and out of the domain, respectively. The efficiency is calculated using Eq(8), and the corresponding calculation is as follows: 𝜀 = 𝐸𝑥𝑡𝑟𝑎𝑐𝑡𝑖𝑜𝑛 ∗ 𝑆𝑝𝑒𝑐𝑖𝑒𝑠 𝑎𝑡 𝐸𝑥𝑡𝑟𝑎𝑐𝑡𝑖𝑜𝑛 𝐼𝑛𝑙𝑒𝑡 ∗ 𝑆𝑝𝑒𝑐𝑖𝑒𝑠 𝑎𝑡 𝐼𝑛𝑙𝑒𝑡 (9) Eq(9) yields efficiency value of 11.4% for 1 MW case study with fresh air entering through openings, whereas the experimental value of the efficiency obtained by Vauquelin and Megret (2002) is 13%. This implies that this method underestimates the efficiency, since a significant amount of contaminated air is entering the domain through opening boundaries. Thus, second method of obtaining efficiency will be checked, which will include contaminated air with species at openings as presented in Eq(10). 𝜀 = 𝐸𝑥𝑡𝑟𝑎𝑐𝑡𝑖𝑜𝑛 ∗ 𝑆𝑝𝑒𝑐𝑖𝑒𝑠 𝑎𝑡 𝐸𝑥𝑡𝑟𝑎𝑐𝑡𝑖𝑜𝑛 𝐼𝑛𝑙𝑒𝑡 ∗ 𝑆𝑝𝑒𝑐𝑖𝑒𝑠 𝑎𝑡 𝐼𝑛𝑙𝑒𝑡 + (𝑂𝑝𝑒𝑛𝑖𝑛𝑔𝑤𝑖𝑡ℎ 𝑠𝑝𝑒𝑐𝑖𝑒𝑠 − 𝑂𝑝𝑒𝑛𝑖𝑛𝑔𝑓𝑟𝑒𝑠ℎ ) (10) Eq(10) yields efficiency value of 13.5% for 1 MW study with air with species concentration of 0.17 entering the domain through openings. This method closely approximates the experimental results. However, it should be noted that this is maximum possible efficiency of the system, whereas the results of the first method show the minimal possible efficiency. Thus, both methods can be used in order to find the range within which the true effective efficiency lies. For 1 MW study the range is: 𝜀𝑙𝑜𝑤𝑒𝑟 𝑏𝑜𝑢𝑛𝑑 (1 𝑀𝑊) = 11.4% 𝜀𝑢𝑝𝑝𝑒𝑟 𝑏𝑜𝑢𝑛𝑑 (1 𝑀𝑊) = 13.5% As the experimental value of efficiency is 13%, the upper bound overestimates by 0.5%, and the lower bound underestimates by 1.6%. Similar procedure is carried out for 4 MW and 10 MW studies. The range for 4 MW is found to be: 𝜀𝑙𝑜𝑤𝑒𝑟 𝑏𝑜𝑢𝑛𝑑 (4 𝑀𝑊) = 18.4% 𝜀𝑢𝑝𝑝𝑒𝑟 𝑏𝑜𝑢𝑛𝑑 (4 𝑀𝑊) = 20.4% The experimental value for 4 MW is 20%, meaning that that the upper bound of the model overestimates the experimental value by 0.4%. The range for 10 MW is obtained in similar fashion and is found to be: 𝜀𝑙𝑜𝑤𝑒𝑟 𝑏𝑜𝑢𝑛𝑑 (4 𝑀𝑊) = 23.4% 𝜀𝑢𝑝𝑝𝑒𝑟 𝑏𝑜𝑢𝑛𝑑 (4 𝑀𝑊) = 24.2% The experimental value of efficiency in this case is 24%, the upper bound overestimates by 0.2%. 1475 4. Discussion and Conclusions Numerical simulation of smoke propagation from a fire in a semi-enclosed space using finite element method is presented. The results of the numerical simulation accurately approximate the experimental values of the extraction efficiency in a tunnel. The method of solving the simulation twice, i.e., with fresh entraining air and with contaminated air, allows the designer to obtain a range for the true value of efficiency. The proposed modified method of obtaining efficiency, using the contaminated entraining air instead of fresh air as usually is prescribed, slightly overestimates the experimental value, whereas the simulation with the fresh air entering the domain through openings underestimates it, due to neglecting the fact that there is a natural suction phenomenon happening near openings across part of these boundaries. This phenomenon of entrainment is expected to be smaller as the tunnel length/height ration increases, but its parametric effect on each given case should be the object of another study. It is observed that as the intensity of the fire source increases the deviation between the proposed contaminated air entrainment modeling and the experimental value is narrowed down. One possible explanation to this could be that when the intensity increases there is an increment in the mass flow rate of the smoke at inlet, while also the extraction experiences increase in the mass flow rate, which in turn makes the effects of the suction less significant. Aforementioned method of obtaining efficiency is able to reproduce the experimental value of the efficiency, and also together with the traditional evaluation of efficiency assuming only fresh air entrainment, provides the range for the efficiency, which takes into account uncertainties associated with the suction phenomena and species concentration in the openings of the tunnel. However, for safety purposes it is recommended to use the lower bound value for the efficiency, in order to design the extraction system for the worst case scenario, and to avoid overestimating the capability of the smoke extraction systems in tunnels. Acknowledgments This publication has been possible with the support and grant by the Shakhmardan Yessenov Foundation. References Beard, A. and Carvel, R. eds., 2005, The handbook of tunnel fire safety. Thomas Telford. Blanchard, E., Boulet P., Desanghere S., Cesmat E., Meyrand R., Garo J. P., and Vantelon J. P., 2012, Experimental and numerical study of fire in a midscale test tunnel. Fire Safety Journal 47: 18-31. Calautit, J. K. and Hughes, B.R., 2014, Wind tunnel and CFD study of the natural ventilation performance of a commercial multi-directional wind tower. Building and Environment 80: 71-83. Calautit, J.K., Chaudhry, H.N., Hughes, B.R. and Sim, L.F., 2014, A validated design methodology for a closed-loop subsonic wind tunnel. Journal of Wind Engineering and Industrial Aerodynamics 125: 180-194. Hussain, S. and Oosthuizen, P.H., 2012, Validation of numerical modeling of conditions in an atrium space with a hybrid ventilation system. Building and Environment 52: 152-161. Hwang, C.C. and Edwards, J.C., 2005, The critical ventilation velocity in tunnel fires—a computer simulation. Fire Safety Journal 40, no. 3: 213-244. Kumar, A. and Mazumder, S., 2010, Coupled solution of the species conservation equations using unstructured finite‐volume method. International Journal for Numerical Methods in Fluids 64, no. 4: 409- 442. Lacroix, D., 2001, The Mont Blanc tunnel fire. What happened and what has been learned. In proceedings of the fourth international conference on safety in road and rail tunnels, held Madrid, Spain, 2-6 april 2001. Lee, S.R. and Ryou, H.S., 2006, A numerical study on smoke movement in longitudinal ventilation tunnel fires for different aspect ratio. Building and Environment 41, no. 6: 719-725. Lin, C.J. and Chuah, Y.K., 2008, A study on long tunnel smoke extraction strategies by numerical simulation. Tunneling and underground space technology 23, no. 5: 522-530. Rohdin, P. and Moshfegh, B., 2007, Numerical predictions of indoor climate in large industrial premises. A comparison between different k–ε models supported by field measurements. Building and Environment 42, no. 11: 3872-3882. Vauquelin, O. & Megret, O., 2002, Smoke extraction experiment in case of fire in a tunnel. Fire Safety Journal, 37, 525-533 in Plastics, 2nd ed. vol. 3, J. Peters, Ed. New York: McGraw-Hill, 1964, pp. 15–64. Versteeg, H.K. and Malalasekera, W., 2007, An introduction to computational fluid dynamics: the finite volume method. Pearson Education. Vittori, F., Rojas-Solórzano, L., Blanco, A.J. and Urbina, R., 2008, Numerical Study of Smoke Propagation in a Simulated Fire in a Wagon Within a Subway Tunnel. In ASME 2008 Fluids Engineering Division Summer Meeting collocated with the Heat Transfer, Energy Sustainability, and 3rd Energy Nanotechnology Conferences, pp. 1017-1024. American Society of Mechanical Engineers. 1476