Assessment of air flow distribution and hazardous release dispersion around a single obstacle using Reynolds-averaged Navier-Stokes equations.

The flow around a cubical building, with a pollution source at the central point of the top of the cube, is studied. The Reynolds-averaged Navier-Stokes and species concentration equations are solved for Reynolds number, Re = 40,000, is based on the height of the cube. The predictions obtained with the standard, the Kato-Launder, and the low-Reynolds number k-epsilon models are examined with various wall functions for the near wall treatment of the flow. Results are compared against Martinuzzi and Tropea measurements (J. of Fluids Eng., 115, 85-92, 1993) for the flow field and against Li and Meroney (J. of Wind Eng. and Industrial Aerodynamics, 81, 333-345, 1983) experiments and Gaussian models for the concentration distribution. It is found that the present unstructured mesh model performs similarly to the structured mesh models. Results from the Kato-Launder model are closer to the experimental data for the flow patterns and contaminant distribution on the cube's roof. However, the Kato-Launder model has an over-prediction for the recirculation zone and the contaminant distribution windward of the cube. The standard k-epsilon and the low-Reynolds number k-epsilon models predict similar flow patterns and are closer to the experimental data of the cube's windward and side face.

Introduction

An accident that releases pollutant hazardous materials inside urban environment could lead to an important environmental pollution and harmful human situations. The airflow distribution and the mass transport phenomena around the urban buildings are defining the hazardous materials dispersion. The prediction of the dispersion of a hazardous pollutant is a difficult, because depends on complex physical phenomena and, contains lot of uncertainty on the determined results. The constrained open space of urban geometries obstructs the hazardous materials, favouring the mass trap inside the city's environment [1]. The accidental risk analysis techniques, define the systematic hazards identification, record the accidental causes and determine the protection measurements [2]. The identification of the pollutant concentration levels and the safety concentration limits is an important problem for the accidental risk analysis method [3].

In the complex terrain of urban areas, such as big cities, the air flow mechanisms and pollutant dispersion cannot be studied experimentally. This is the reason that simplified models should be used to help us understand the basic flow mechanisms. The flow around a cube is a widely-used experimental and numerical benchmark problem which could be used to study the airflow characteristics and the plume dispersion around isolated buildings.

Several experiments have been conducted so far regarding the flow around cubical geometries. Castro and Rodin [4] studied the dependence of the wake recirculation and the roof's vortex for uniform and turbulent approaching flows. The Reynolds number effect in the flow characteristics are examined by Lim and Castro [5]. Several experiments for fully turbulent flows have been examined [6, 7, 8]. The cost of field experiments and wind tunnel studies has led researchers to use Computational Fluid Dynamics (CFD) techniques more often in order to study the flow around similar geometrical shapes. Several studies examine the flow around isolated orthogonal geometries with different approaches [9, 10, 11, 12] and the pollutant dispersion around them [13]. Many studies exist for more complex orthogonal geometries such as a street canyon [14, 15, 16]. The usual technique is the numerical solution of the Reynolds-averaged Navier-Stokes equations (RANS or URANS for steady-state or the unsteady formulation, respectively). Lakehal and Rodi [17] studied the flow around a cube with several k-epsilon models and compared the results with Martinuzzi's experiments. Zhang, et al. [18] examined different types of flows around a building with a k-epsilon model and compared the results with Castro's experiments. Mochida, Tominaga et al. [19] examined different k-epsilon models for the flow around high-rise buildings with a focus on the basic boundary conditions. Yoshie, et al. [20] studied the differences between a modified k-epsilon model and a standard k-epsilon model to validate the flow characteristics around a cube and found that the modified k-epsilon model is more accurate for areas with high velocities than for regions with lower ones.

The study of the flow field around cubical geometries is significant to understand the pollutant dispersion mechanisms around buildings. In this direction, several wind tunnel experiments can be found that consider the pollutant dispersion from different source positions and different wind-induced angles [21, 22]. Robins and Castro [23] studied the plume dispersion around a cube. Different experiments have been realized for different source positions, various speeds of emission release and various wind directions. They clarify experimentally the importance of the turbulent shear stress and how it affects the pressure and the flow distribution. Thompson and Roger [24] studied the pollutant dispersion for different height and source location. Higson, et al. [25] compared field experiments against wind tunnel measurements and defined that the more massive plumes and their fluctuations are better defined in a real field area than in wind tunnel experiments. Mavroidis, et al. [26] studied the flow around cylindrical and rectangular obstacles, comparing field experiments against wind tunnel measurements. The mean concentration in the wind tunnel is found to be higher than the one in the real field, and this is mainly due to different turbulent scales.

A promising technique in turbulent calculations is the Large Eddy Simulation (LES). LES approaches to resolve the unsteady motion for the flow around a cube and define complex phenomena such as the vortex shedding, the dynamic loading, turbulence fluctuations, etc. [27]. The disadvantage of this method is that it is more expensive than the RANS methods. LES could be an option for the study of flow around isolated bodies such as the flow around a cube [28]. LES models the sub-grid (SGS) stress tensor and the scalar fluxes with a lot of detail, however they have the disadvantage of the computational cost because of the high-resolution requirement. Another much more uncomplicated and computationally more adequate approach is the Implicit Large Simulation (ILES) [29, 30, 31, 32, 33]. Direct Numerical Simulations method (DNS) is the most precise numerical approach but so computationally costly that cannot be adapted to geometries such as in an urban environment.

Similarly, numerical experiments have been conducted for the pollution distribution around cubical geometries. N. Meroney, et al. [34] observed that when the Reynolds-averaged turbulence models are used, the pollutant dispersion around bluff bodies is over-predicted. In addition, they demonstrated that the separation and the reattachment area, of the wind flow around a body, can be appropriately estimated by an adaptive grid without any extensive calculation time. Delaunay, et al. [35] compared numerical simulations against wind tunnel measurements around a rectangular building containing chimney outlets. The flow recirculation around the roof-top, the windward and the leeward sides of a building with a second-order turbulence model is reproduced. This model provided good estimations for the gas concentration on both sides of the building.

Tominaga and Stathopoulos [36] examined different types of k-epsilon models around a cube and the dispersion of a flush vent which is located on various points at the top of a cube with different air oblique angles. The study exhibited an underestimation of the concentration on the leeward and the lateral sides of the cube due to reduced diffusivity and concluded that the velocity field determines the concentration field. These results were successfully compared against the Li and Meroney [21] wind tunnel experiments. An essential problem of the urban aerodynamics is the need for large computational domains and thus the increase of mesh requirements. The flow characteristics define the grid resolution to achieve accurate and computationally efficient solutions. Different grid types could be applied to discretize an urban domain: Cartesian grids and unstructured grids. Generating an unstructured grid for complicated geometries is a significantly faster automated process compared to a structured one. On the other hand, all the numerical methods exhibit dependency on the quality of the grid employed, and since unstructured grids are mainly used for complicated geometrical definitions, they can consist of various element types and significant variations of grid quality can be noticed. This, in turn, can pose several challenges for the numerical methods and their formulations when employed in the RANS framework as identified by Antoniadis et al. [37] using unstructured meshes. However, it is not realistic to apply structured grids into a complicated geometry [38]. Non-orthogonal cells should be avoided near wall treatment cells and boundary conditions. Prismatic or hexagonal cells are preferred for the near wall turbulence models [39]. Blocken, et al. [40] have studied the grid resolution for a passage between two equal height buildings. A hybrid mesh which is the combination of a structural hexahedral grid and an unstructured tetrahedral grid is tested. The structured grid is applied at the building's roof, the building's sides and the passage. The unstructured grid is used for the connection of the building's roof and the domain's center. Van Hooff and Blocken [41] studied a combination of the urban wind flow to a building's natural ventilation approach. They described a body-fitted grid generation process for the complex internal and external environments. Their method defines the steps of the grid generation and required resolution to control the quality of the results. This process is in contrast to a semi-automatic unstructured grid generation procedure which can translate to insufficient control of the grid resolution, the volume skewness, the grid stretching and the aspect ratio. Gargallo-Peiró, et al. [42] are presenting a method which can introduce the methodology of meshing the geometry and the landscape of a city. They used 3D unstructured tetrahedral meshes for non-viscous urban simulations for flows around blocks of a city.

The focal point of the present work is to assess the performance of various RANS models using unstructured meshes, which is a combination that is at the fidelity level still used by the industry associated with insurance policy and the risk assessment of an urban environment, for the well-established test problem of air flow distribution around a cube. Simulation of the flow dispersion in the complex environment of an entire city would require the use of unstructured meshes and RANS models; therefore the subject work is relevant for assessing the performance of this combination and understanding which is the best combination to be employed for modeling an entire city.

In the present study, the turbulent model's capability to predict the contaminant distribution from a central vent release on the cube's roof with the most optimum manner is investigated. The predictions from different turbulence models of the flow field characteristics that affect the contaminant distribution around the cube are assessed and compared against a Gaussian dispersion model and validated against experimental data and numerical results. The detailed experiments of Martinuzzi and Tropea [6] are used for the comparison of the experimental flow characteristics and Li and Meroney [21] experiments are used for the comparison of the contaminant distribution around the cube. The model and simulation details are given in Section 2. Results and discussions of the flow patterns and contaminant distributions can be found in Section 3, and conclusions are listed in Section 4.

Theory/calculation

Flow field description

The applied computational field is presented in Fig. 1. The upstream computation length is 5H, the downstream computational length is 10H, the lateral width is 11H and the total height is Z = 2H, where H is the height of the cube in the base of which the axes origin is considered. The guidelines of the German Association of Engineers (VDI) [43] suggest that the maximum blockage should be below 10%. The size of the selected computational domain has a 5% blockage effect.

Fig. 1

Computational domain and boundary conditions. Letters indicates monitoring positions: A (X:-5H, Y:0, Z:0), B(0,0,H), C(0.5H,0,H), D(H,0,0), E(1.5 H,0,0), F(2.5H,0,0), G(3.5H,0,0).

The results of the simulations are validated against the Martinuzzi and Tropea [6] experimental results that were conducted in a 156 H × 24 H × 2H wind tunnel in fully turbulent flow with Reynolds's number equal to 40,000, based on the cube's height. For our simulations, we also kept the same Reynolds number.

The vent location is located at the central point of the cube's roof at point C of Fig. 1. Α point source with passive vent plumes and a low minimum dilution criterion is examined. The wind orientation is kept streamwise at . The plume dispersion results are compared against the experimental data of Li and Meroney's [21]. The experimental data for the concentration are collected at Re = 11,050, and since the critical Reynolds number for concentration variations is for Re = 11,000 [44], no changes on the concentration distributions are expected.

The study of this obstacle at a Reynolds number of 40,000 and a cube of 0.025 meters height corresponds to contaminant dispersion around a typical rectangular building in an urban landscape on a reasonably calm day, which however is challenging for turbulence models due to the low Re turbulent flow regime.

Turbulence models

The Reynolds-averaged Navier-Stokes (RANS) equations are used for the simulations of the present study, where three different turbulence models are tested: the standard k-epsilon model, the Kato-Launder standard k-epsilon model [45], and the low-Reynolds k-epsilon model. Gorji, et al. [46] summarized the model constants, the damping functions and near wall correction functions. The general form of the k-epsilon model with the Boussinesq hypothesis for the swirl turbulent viscosity can be written as:Where is the velocity in the three directions, is the density of the fluid, is the turbulence kinetic energy, is the dissipation rate, is the viscosity, is the turbulent viscosity, are models constants, is the turbulence kinetic energy production due to mean velocity fluctuations and described as , is the velocity fluctuations. D and E are near wall correction functions for k and epsilon equations, respectively.

For the standard k-epsilon model [47], , where the turbulent viscosity is and the rate-of-strain tensor, , is given by , , , D = E = 0, .

For the Kato-Launder model, , where the vorticity rate is and the vorticity tensor is , and , , , D = E = 0, .

In the Yang and Shih low-Reynolds k-epsilon model [46]: the turbulent viscosity is , , , , and , , where , and , . The damping function is, where the turbulent Reynolds number is , where .

Results from these three models were compared in the present work for the selection of the most suitable model in order to predict the flow around the cube and the hazardous release dispersion. The standard k-epsilon model over-predicts the turbulence kinetic energy due to an excessively generated term of at the stagnation points. In the Kato and Launder model the excessive turbulence energy production results from the critical shear strain rate S in the stagnation regions. For this reason an alternative form of the turbulence energy production, which introduces the term of the vorticity tensor rate is proposed. This model gives better results for flows around bluff bodies [45]. Finally, the low Reynolds model approach specifies ε and the turbulent viscosity in the near-wall cells. The domain is subdivided into two different regions, the fully-turbulent region and a viscosity-affected region, which are determined by the turbulent Reynolds number.

Species transport equation

The species transport equation of a contaminant concentration is expressed by:where c is the concentration and is any source term inside the flow field.

The species diffuse due to the turbulence is expressed as:Where is the diffusion coefficient for i species in the mixture and is the turbulent Schmidt number, an empirical number which may vary in general between 0.2-1.3 and plays an important role for the calculation [48]. In the present simulations, the turbulent Schmidt number is kept constant at 0.7.

The mean concentration is studied through the non-dimensional concentration coefficient K that is defined as [21, 49]:Where, is the measured tracer concentration, is the source tracer concentration, is the contaminant volumetric flow-rate and is the velocity at the building's height.

Boundary conditions

The velocity distribution at the inlet is defined with a logarithmic profile [50] obtained from the experimental data:where , is the bulk (average) velocity and H = 0.025m is the cube height. The inlet turbulence kinetic energy is calculated using experimental data from the Journal Engineering databank [6].

The inlet boundary conditions for the k and epsilon profiles are well described by Breuer, et al. [51]. The turbulence kinetic energy profile is expressed as where . The dissipation rate is described as , where is the turbulence length scale and set as 0.1 H.

No-slip boundary conditions are applied at the bottom, top and cube walls and two different wall function approaches are used. At the top of the computational domain, a wall boundary is applied. The standard wall functions that are based on the theory of Launder and Spalding [52], for the standard k-epsilon and standard k-epsilon with Kato-Launder models.

At the lateral sides, periodic boundary conditions are applied and usual outflow conditions are applied at the outlet, where the pressure is kept equal to zero and the streamwise derivatives of all other quantities are vanished. These boundary conditions are set far enough downstream from the cube location.

A passive scalar simulates the pollutant release from an orthogonal source on the centre and at the top of the cube with area. The exhaust velocity is kept equal to 1.54 m/sec in all simulations. This velocity is low enough to avoid the jet effect phenomenon [53].

Mesh type

Two types of meshes are used. Outside the boundary layer, an unstructured tetrahedral grid is used and, inside the boundary layer a prism mesh. The first cells at the walls are at m and the expansion ratio for the prism cells is 1.3 which corresponds to 1 < y+ < 5. This range of y+ satisfy the minimum values for the low-Reynolds model which is the most restricted. he log-law for the mean velocity near the walls is applied when and the laminar stress-strain relationship is applied for lower values. Furthermore, the enhanced wall treatment is applied for the standard k-epsilon with the low Reynolds model. This near-wall modeling method combines the two-layer model with enhanced wall functions.

The near-wall mesh is fine enough to resolve the laminar sublayer. As shown in Fig. 2, a grid independence test is conducted based on the standard k-epsilon model with the standard wall functions. It is found that the solution is grid-independent for 4,023,449 cells, where the maximum velocity difference is less than 0.5% from the finer case tested at the location X/H = -5. The coarse grid has 2,036,242 cells and the finer grid has 8,597,367 cells.

Fig. 2

A grid independence test based on the standard k-epsilon model with standard wall functions, a) Inlet velocity, location X/H = -5 b) location X/H = 0.5.

Numerical schemes

The CFD flow solver Ansys Fluent 17 is used for the flow calculation. The SIMPLE scheme is used for pressure and velocity coupling. The nonlinear terms are calculated with a second order upwind scheme, and second order schemes are used for the calculation of all the other terms. The converge criteria is kept less than based on the absolute error of all quantities.

Results & discussion

Flow field results

A 3D geometry with three different turbulence models is compared: the standard k-epsilon model with standard wall functions (St-ke-WF), the Standard k-epsilon with the Kato-Launder model (St-ke-KL) and the Standard k-epsilon with low Reynolds model (St-ke-low-Re). The results obtained from these turbulence models are compared against Martinuzzi and Tropea [6] experimental data in order to validate the flow field and Li and Meroney [21] experimental data for validation of the pollutant dispersion around a cubical building. The comparison of the flow field is made in terms of velocity and turbulence kinetic energy distributions and flow patterns. The flow fields at the symmetry plane, where the main recirculation zone exists and the major separation points formed are presented in the streamlines plots of Fig. 3a,b and c, for the St-ke-WF model, St-ke-KL model and St-ke-low-Re models, respectively.

Fig. 3

Streamlines of the flow at the symmetry plane for a) St-ke-WF model, b) St-ke-KL model, c) St-ke-low-Re model, d) Characteristic separation lengths.

The calculated lengths of the main separations points are: for the upstream locations of the cube, for the downstream, for the roof, is the roof's recirculation height. All are illustrated in Fig. 3d and summarized in Table 1. The centre of the vortex is also summarized in Table 1.

Table 1

Main separation point lengths.

Case	Model	X_f	X_r	X_b	Z_r	Center of the vortex
Standard k-epsilon with standard Wall functions	St-ke-WF	0.6 H	0.47 H	2.2 H	0.1 H	X = 1.55 HY = 0.59 H
Standard k-epsilon with Kato-Launder model	St-ke-KL	0.55 H	-	2.46 H	-	X = 1.5 HY = 0.76 H
Standard k-epsilon with low Reynolds model	St-ke-low-Re	0.4 H	0.42 H	2.3 H	0.1 H	X = 1.56 HY = 0.78 H
Martinuzzi and Tropea [6] experimental data		1.04 H	-	1.61 H	0.17 H	X = 1.5 HY = 0.93 H
Lakenhal and Rodi (1997)	k-epsilon with wall functions	0.651 H	0.432 H	2.182 H	-	X = 1.58 HY = 0.72 H

As the flow approaches the leeward surface of the cube, the main separation vortex appears. At the leeward cube's surface, the boundary layer detaches and binds the cube by forming the well-known shape of a horseshoe vortex [54]. All the k-epsilon models underestimate this recirculation zone which experimentally is found to extend to . The St-ke-WF model predicts the closest to the experimental separation point at . The St-ke-KL model predicts almost a similar separation point at H and the St-ke-low-Re model gives the worst prediction of the separation point at . Moreover, it is found that the St-ke-KL model calculates a long separation zone and does not predict a reattachment point on the cube's roof. This is in agreement with the experimental data and defines the pollutant dispersion. Results obtained using the St-ke-low-Re and St-ke-WF, models are quite similar i.e.: a small recirculation zone and a reattachment point are found at positions and for the St-ke-low-Re and the St-ke-WF models respectively. According to Table 1, the St-ke-KL model has a better approximation for the flow around the cube since it does not present a reattachment point at the roof of the cube.

The velocity profiles along Z at points C and D (see Fig. 1) are illustrated in Fig. 4a and b. It is found that the St-ke-WF and St-ke-low-Re models predict the velocity profile with a significant difference close to the wall in comparison to Martinuzzi and Tropea experiments. In contrast, similar behavior of the St-ke-WF and St-ke-low-Re models is found to the numerical results of Lakehal and Rodi [17] where the standard k-epsilon model is used in a structured mesh solver. Results obtained by the St-ke-KL model are closer to the experimental ones and the reverse flow is better predicted. At the windward edge of the cube's roof point, only the St-ke-KL model predicts the reverse flow at the edge compared to the other models where the reattachment flow is situated before the cube's centre point.

Fig. 4

Mean velocity (upper) and turbulence kinetic energy (bottom), along Z, normalized by and , respectively, for the locations X/H = 0.5 (a and c) and X/H = 1 (b and d).

Fig. 4c and d illustrate the turbulence kinetic energy distribution at the position C and D, respectively. Results at the C location from the St-ke-WF and St-ke-low-Re models are observed again to be similar; however, results of the St-ke-KL are found to be in better agreement with the experimental data (except at the peak value) due to the reduced calculated turbulence kinetic energy. The difference in the turbulence kinetic energy peak between models and experiments is more significant at location D. At this position all models underestimate the turbulence kinetic energy compared to the experimental data.

All the k-epsilon models over predict the length of the cavity recirculation zone as a result of the underestimation of the turbulence kinetic energy on all cases. Using the St-ke-WF model we obtain results closer to the experimental data with . The St-ke-KL model gives the worst overestimated recirculation length with , and finally, the St-ke-low-Re model predicts the separation at H. The windward area behind the cube could be separated into two different zones. The cavity zone (at location E) where the recirculation zone appears with low velocities and high turbulence [55] and the near-wake zone situated after the cavity. The velocity distribution from the unstructured grid calculation at the cavity zone and at the near-wake zone are in good agreement against the experimental data of Martinuzzi and Tropea [6] and the numerical results from Lakehal and Rodi [17] as shown in Fig. 5. This result is significant because it shows that the present unstructured solver has the potential to simulate flows in urban environments, despite their complex geometries, with the same accuracy as with a structured mesh solver. Fig. 5 a,b and c show the non-dimensional mean velocity distribution at positions E, F and G, respectively. At positions, E and F the mean velocity distribution is in good agreement with the experimental data. At G position the velocity distribution is in agreement at the part above the cube height's and shows differences at the lower part.

Fig. 5

Mean velocity at the symmetry plane for the locations a) X/H = 1.5, b) X/H = 2.5 c) X/H = 4.

The identification of the coherent structures and the vortices can be made with the iso-surfaces of the Q-criterion. The definition of the Q-criterion is [56]:where , is a constant for the impressions with a value , S is the Strain Rate, Ω is the vorticity rate.

Fig. 6 illustrates the iso-surface of the Q criteria, Q = 0.1 , for the three different models. It can be seen that the horseshoe vortex is formed on the leeward face of the cube. The horseshoe vortex has a considerable downstream extension for the St-ke-KL (Fig. 6 b) compared to the others models due to the vorticity based formulation of its production source term. The arc-shaped vortex on the leeward face of the cube is similar for all the three cases. The flow inside the recirculation zone is a strong mixing and turbulence generation region.

Fig. 6

Isosurface of the Q = 0.1 criteria for a) the St-ke-WF model, b) the St-ke-KL model, c) the St-ke-low-Re model.

Concentration

In order to define the advantages and disadvantages of each numeral model, the distribution of the non-dimensional concentration coefficient, K, of the pollutant release around the cube for all the different numerical model is compared. The K distribution is also compared against the Huber and Snyder model (Fig. 7d) and the Li and Meroney's experimental data (Fig. 7e).

Fig. 7

K distribution at the symmetry plane for a) St-ke-WF, b) St-ke-KL, c) St-ke-low-Re d) Huber and Snyder models and e) the experimental data Li and Meroney (1983).

Huber and Snyder [57] studied the wake effect for short emitted sources and developed a Gaussian equation model to predict the dispersion and transportation of emitted plumes from buildings. According to this model, that is suitable for flows of Re > 36,000, the non-dimensional coefficient distribution can be described from the equation:and the dispersion parameters can be calculated from the expression: is the source height that is situated in the cube, and in our case since the pollutant is released from the top of the cube.

Fig. 7, shows some qualitative characteristics regarding the dispersion of the non-dimensional concentration coefficient, K with isopleth graphs. In order to understand the mechanism of the concentration dispersion, the study of the mass diffusion is realized. According to Eqs. (3) and (4) the concentration is treated as species which is transferred by the advection-diffusion equation. The convective transfer of the mean concentration is defined as , the turbulent concentration flux is defined as , and the molecular diffusion flux is defined as .

Furthermore, in order to understand the mechanism of the concentration dispersion around the cube, the convective concentration fluxes (Fig. 8) and the turbulent concentration flux (Fig. 9) for the streamwise (x-direction) at the symmetry plane is realized for all the numerical cases. As shown in Fig. 8, the pollutant is primarily transferred by convection downstream of the cube and only a small quantity of the pollution reaches at the cube's cavity zone.

Fig. 8

Convective mass flux at the symmetry plane for a) St-ke-WF, b) St-ke-KL, c) St-ke-low-Re models.

Fig. 9

Turbulent mass flux at the symmetry plane for a) St-ke-WF, b) St-ke-KL, c) St-ke-low-Re models.

The St-ke-WF, St-ke-low-Re and Huber and Snyder models present similar distributions for the concentration coefficient, K, at the constant value of 1 (Fig. 7 a,c). On the other hand, the St-ke-KL model (Fig. 7b) presents a larger concentration length leeward the cube, at the constant value of 1. The Kato-Launder wall functions at the cube's surfaces have as result the negative values of the convective mass flux at its rooftop and present a larger concentration length due to the higher convective mass flux (Fig. 8b). Negligible differences at the convective mass flux distribution at the cube's roof between the St-ke-WF and St-ke-low-Re models have as result small differences at the concentration's length.

The St-ke-KL model presents important concentration values near the ground (Fig. 7b) because the turbulent mass flux area is extended in a higher limit at the cube's height and traps more pollutant into the cavity and wake area (Fig. 9b). Huber and Snyder models cannot predict the downwash effect that brings higher concentration near the ground (Fig. 7d). Li and Meroney's experimental data present a smaller extension for the pollutant concentration comparing to all the present numerical models and the Huber and Snyder model.

The k-epsilon models predict different K distributions also in the roof of the cube (Fig. 10) because of the velocity distribution and turbulence kinetic energy differences. The St-ke-WF and St-ke-low-Re models predict higher values of K windward of the release vent. Instead, the St-ke-KL model predicts a higher concentration area at the upwind vent area. The concentration distribution predicted by the St-ke-KL model is in better agreement with the Li and Meroney [21] experimental data. This agreement is due to the better calculation of the reverse flow at the top of the roof since the reattachment point is windward to the vent location. For this reason, the pollutant is trapped into the recirculation zone. Fig. 11 shows a 3-dimensional view of the non-dimensional concentration dispersion for the constant value of 1. The negative values of the convection mass flux (Fig. 8b) are moving the plume direction opposite of the wind direction.

Fig. 10

K distribution at the cube's roof a) St-ke-WF, b) St-ke-KL, c) St-ke-low-Re models and d) the experimental results from Li and Meroney (1983).

Fig. 11

Isosurface for K = 1 distribution at the cube's roof for a) St-ke-WF, b) St-ke-KL, c) St-ke-low-Re models.

Fig. 12 shows the contours of K on the downstream wall of the cube. The dominant experimental K values (Fig. 12d) on the windward face of the wall lie between 0.5 and 1.0 which is in a better agreement with the St-ke-KL model distribution (Fig. 12b). The St-ke-WF (Fig. 12a) and St-ke-low-Re (Fig. 12c) models calculate similar K distribution that differs from the experimental concentration one. Higher values of K are found mostly in the center area of the windward face of the cube and are presenting an expansion towards the cube's base.

Fig. 12

K distribution at the cube's windward face for a) St-ke-WF, b) St-ke-KL, c) St-ke-low-Re models and d) the experimental results from Li and Meroney (1983).

The K concentration at the side wall of the cube is illustrated in Fig. 13. The experimental data shows that the main concentration is transported from the roof to the upper part of the side wall (Fig. 13d). The present results (Fig. 13 a, b, c) under-estimate the concentration at the upper part of the side wall and found important concentration at the lower part of the wall which is transferred from the leeward face of the cube. The K distribution is quite similar for the St-ke-WF (Fig. 13a) and St-ke-low-Re (Fig. 13c) models, which present higher concentrations at a small area at the low part of the side wall cube. The St-ke-KL model (Fig. 13b) concentration values are at the lower part of the side wall of the cube but there is much higher dispersion at the side wall of the cube.

Fig. 13

Non-dimensional concentration coefficient distribution at the cube's left face a) St-ke-WF, b) St-ke-KL, c) St-ke-low-Re models, and d) the experimental results from Li and Meroney (1983).

The K level decreases as the flow passes through the cube edges. The plume concentration distribution follows the air flow behavior. The concentration rapidly decreases after passing the cube edge as is shown in Fig. 14. At the edge point of the cube, a significant underestimation of the non-dimensional concentration coefficient appears for all the k-epsilon models compared to Li and Meroney's experimental data. Predictions from St-ke-WF and St-ke-low-Re models are in better agreement with the experimental data while the one from the St-ke-KL model underestimates K.

Fig. 14

K distribution on the cube's roof and leeward face at the symmetry line (Distance line is marked in the right cube).

Different lateral isopleths planes are examined at the positions E, F and G in Fig. 15. The highest K value is observed at the centre of the isopleths at the source height. Moving downstream from the release source the values decrease and the K distribution is expanding laterally and longitudinally. Increasing the distance windward of the cube, all numerical results (Fig. 15 a,b,c,d,e,f,g,h,i) appear to have higher K values than in the Li and Meroney experimental data (Fig. 15 m, n, o). The K distribution for the Huber and Snyder model (Fig. 15 j, k, l) is denser around the plume centerline than the other results. Finally, the K distributions, as predicted by the numerical results, show a crucial vertical dispersion towards the lower part of the cube, contrary to the experimental results that show an important lateral dispersion.

Fig. 15

K profiles for the St-ke-WF model: (a) at X/H = 1.5, (b) at X/H = 2.5, (c) at X/H = 3.5), the St-ke-KL model: (d) at X/H = 1.5, (e) at X/H = 2.5, (f) at X/H = 3.5, and St-ke-low-Re model: at (g) at X/H = 1.5, (h) at X/H = 2.5, (i) at X/H = 3.5, the Huber and Snyder model: (j) at X/H = 1.5, (k) at X/H = 2.5, (l) at X/H = 3.5), and the experimental results from Li and Meroney (1983): (d) at X/H = 1.5, (e) at X/H = 2.5, (f) at X/H = 3.5.

Fig. 16 illustrates the K variation in the Z direction for different positions behind the cube. At point E, inside the recirculation zone, the prediction of the St-ke-KL model is in better agreement with the experimental data than the other two k-epsilon models. Moreover, K distributions from the St-ke-WF and St-ke-low-Re models predict similar behavior and an important dimensional concentration increase is found slightly above the cube height. At point F, near the limit where the recirculation zone ends, the St-ke-WF and St-ke-low-Re models have similar behavior and are in a better agreement with the experimental data. The St-ke-KL model overestimates the prediction of K compared to the other two models. At point G, all numerical models are in a good agreement with the experimental data. The St-ke-KL shows a small overprediction of the K values.

Fig. 16

Non-dimensional concentration coefficient distribution at a) X/H = 1.5, b) X/H = 2.5 and c) X/H = 3.5.

Conclusion

In this research, the flow around a cube with a contaminant source release at the roof is tested with different k-epsilon turbulent models. The obtained results are compared against the experimental data of Martinuzzi and Tropea [6] and Li and Meroney [21], the numerical results of Lakehal and Rodi [17], as well as the Gaussian model of Huber and Snyder [57]. All k-epsilon models are found to underestimate the flow characteristics around the cube something that affects the inflation of the pollutant dispersion around the cube.

At the cube's roof the St-ke-KL model predicts a long separation zone and does not have a reattachment point on the top of the cube's roof. The St-ke-WF and the St-ke-low-Re models predict a small recirculation zone and a reattachment point at a position which is situated before cube's centre. This estimation led to high concentrations windward of the release vent which is not confirmed from the experimental data.

All the k-epsilon models over-predict the length of the cavity recirculation zone. St-ke-KL gives the most overestimated recirculation length and the St-ke-WF gives results closer to the experimental data. St-ke-low-Re and St-ke-WF show similar recirculation zones. This long recirculation length results from the underestimation of the turbulence kinetic energy ( term).

The concentration level decreases as the flow passes through the cube's edges. At the edge point of the cube an important underestimation of the dimensional concentration appears for all the k-epsilon models compared to Meroney's experimental data. St-ke-WF and St-ke-low-Re models are in better agreement with the experimental data. The St-ke-KL plume dispersion is in better agreement with the Meroney's experimental data and shows a more diffusive main core than the St-ke-WF and St-ke-low-Re models.

According to the symmetry plane, the St-ke-WF and St-ke-low-Re models have similar dimensional concentration lengths to the Huber and Snyder model. St-ke-KL dimensional concentration length is slightly longer. The experimental data results give smaller lengths than RANS and Huber and Snyder model. But Huber and Snyder's model fail to predict the downwind shift of the dimensional concentration. On the other hand, RANS models are in better agreement with the experimental data. Comparing different lateral isopleths planes behind the cube, St-ke-WF is in better agreement with the experimental data compare to the other k-epsilon models near the wake cube. In contrast, moving away behind the cube the St-ke-low-Re is in a better agreement with the experimental data.

In order to define the hazardous release dispersion for safety approaches, it is important to study the advantages and disadvantages of each model. The St-ke-WF and St-ke-low-Re failed to accurately predict the central roof hazardous material release. St-ke-KL has a better approach and could be an option for this kind of problems. On the contrary, St-ke-KL overpredicts the hazardous zone compare to the other two models that are in a better agreement with the experimental data. None of the examined models were able to satisfactorily predict the lateral dispersion of the pollutant at the sides of the cube, as measured in the experiment of Li and Meroney [21]. St-ke-WF and St-ke-low-Re are in better agreement with the experimental data of the non-dimensional concentration variation with distance from the edge of the cube. St-ke-KL over predicts the non-dimensional concentration. It is found that the St-ke-WF and St-ke-low-Re models give a better approximation for the hazardous release dispersion windward of the cube, but the St-ke-KL model is better for the dispersion at the cube's roof.

Declarations

Author contribution statement

Konstantinos Vasilopoulos: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Ioannis E. Sarris, Panagiotis Tsoutsanis: Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data.

Funding statement

This work was supported by the Centre for Computational Engineering Sciences at Cranfield University (EEB6001R).

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.