CFD Study on the Influence of Atmospheric Stability on Near-field Pollutant Dispersion from Rooftop Emissions

2. 1 Numerical Methods

The Fluent CFD software package (Fluent, 2012) was used to simulate wind flow and pollutant dispersion around a cubic building. Reynolds-averaged conservation equations for mass and momentum were used to simulate the processes of interest. Mass and momentum conservation equations are written as follows;

where u_j is the velocity of j component, t is the time, x_j is the j coordinate, ρ is the air density, μ is the dynamic viscosity; g_i is the gravitational body force.

The SST k-ω model is selected for closure of the system of equations composed of the continuity equation and the Reynolds-Averaged Navier-Stokes (RANS) equations. The SST k-ω model of the Fluent software uses two different transport equations to express the turbulence kinetic energy (TKE) k and the specific dissipation rate ω. The SST k-ω accounts for the principal turbulent shear stress and uses across-diffusion term in the ω equation to blend both the k-ω and k-ε models and to ensure that the model equations behave appropriately in both the near-wall and far-field zones. Thus, the SST k-ω model offers superior simulation performance compared with the individual k-ω and k-ε models (Menter et al., 2003).

The SST k-ω model has the following two transport equations (Fluent ver. 14, 2012):

and

In these equations, $$ {\tilde{G}}_k$$ represents the generation of turbulence kinetic energy due to mean velocity gradients, G_ω represents the generation of ω, and Γ_k and Γ_ω represent the effective diffusivity of k and ω, respectively. Y_k and Y_ω represent the dissipation of k and ω due to turbulence. D_ω represents the cross-diffusion term. S_k and S_ω are user-defined source terms.

The pollutant dispersion patterns were analyzed after solving the species transport equation in conjunction with the turbulence model equations. The advection–diffusion (AD) module was applied to study the species transport process by analyzing the mass fraction of pollutants in the mixture. Fluent analyzes the mass diffusion process based on the following equations:

Where J_i is the diffusion flux of the species i (kgㆍm^-2ㆍs^-1);

In these equations, ρ is the density of the mixture (kgㆍm^-3), D_i,m is the mass diffusion coefficient of the pollutant in the mixture (m²ㆍs^-1), D_T,i is the thermal diffusion coefficient (m²ㆍs^-1), Y_i is the mass fraction of the species i (kgㆍkg^-1), μ_t is the turbulent viscosity (kg·m^-1ㆍs^-1). Turbulent Schmidt number Sc_t was specified as 0.1 in this study.

2. 2 Computational Domain and Boundary Conditions

The computational domain and boundary conditions are shown in Fig. 1. Dimensions of computational domains are equal to 25H×11H×6H (H=10 m is building height) in the streamwise (x), lateral (z) and vertical (y) directions, respectively. These values were selected by following the COST Action 732 (Franke et al., 2007) and AIJ (Tominaga et al., 2008) guidelines to limit the deterioration of the prescribed inflow profiles along the empty fetch upstream of the building.

Fig. 1.

Computational domain and boundary conditions

The inlet velocity profile of the horizontal wind velocity is defined (Pieterse and Harms, 2013) as:

where $$ x={\left(1-\frac{16z}{L_{MO}}\right)}^{1/4}$$, z(m) is height from the surface. L_MO is Monin-Obukhov length. The friction velocity u_∗ = 0.3 mㆍs^-1 and roughness length z₀ = 0.13 m were used in this study.

As illustrated by Lin et al. (2007) the vertical profile of temperature can be defined as follows:

where, T(z) is the air temperature (K) at $$ z,x_0={\left(1-\frac{16z}{L_{MO}}\right)}^{1/2},x_s={\left(1-\frac{16z_s}{L_{MO}}\right)}^{1/2}$$, z_s (1.35 m) is a height above the earth’s surface, ϒ_d (0.01 Kㆍm^-1) is the dry adiabatic lapse rate of T_s, (K) is the air temperature at z_s and g (9.8 mㆍs^-2) is the gravitational acceleration constant.

The profile of k for thermally stratified atmospheric boundary layers (Pieterse and Harms, 2013) is:

The turbulent dissipation of inlet (Pieterse and Harms, 2013) is defined as follows:

When using a commercial CFD code, Blocken et al. (2007) provide the wall boundary condition by deriving a relationship which brings the rough wall functions in equilibrium with the inlet profiles. For Fluent this relation is given by:

Where k_s is the roughness height and C_s is a roughness constant required for the wall function. Roughness height should be smaller than the height of the center point of the wall adjacent cell (=0.25 m) and was consequently set to 0.19 m in this study. The resulting value for C_s is 6.7 is defined through a User Defined Function. At the ground surface, the non-slip condition and zero heat flux are assumed.

2. 3 Simulation Plan

Twelve simulations were designed to test the effect of atmospheric stability on plume dispersion. Two emission velocities are used in this study. To compare the profile of velocity and streamline, emission velocity 0 ms^-1 was used. Hydrogen sulfide was selected as a pollutant. Emission velocity 0.619 ms^-1 was used to simulate the concentration field of wind tunnel results. As shown in Table 1, six atmospheric stabilities are selected using the Monin-Obukhov length scale L_MO. The simulated flow and concentration results of neutral stability condition (D) were compared with the wind tunnel experiments.

Table 1.

Simulation plan to test the effect of atmospheric stability conditions.

3. 1 Grid Sensitivity Study

To test grid sensitivity, three grids with cubical cells of 1.0, 0.5 and 0.4 m around building surfaces were used in this study. Franke et al. (2004) recommend the use of hexahedral cells over tetrahedral cells, as hexahedral cells yield smaller truncation errors and better iterative convergence. Therefore, this study selected hexahedral cells. Franke et al. (2004) also recommend keeping the stretching ratio below 1.3 in regions of high gradients, to limit the truncation error. Thus, outside the building block array an expansion ratio of 1.1 was applied. Fig. 2 shows the results of the grid sensitivity study that presents the comparison of the mean horizontal wind velocity U along vertical profiles at two locations with different flow regimes. The first profile (Fig. 2(a)) was located upstream (X/H=-1.0) of cubic building. The second velocity profile (Fig. 2(b)) was located in the leeward wake region (X/H=3.0). As shown in Fig. 2, the profiles show some grid dependency, which becomes smaller on the finer grids. For both locations, very little quantitative differences are present between the 0.4 m-grid and the 0.5 m-grid, whereas for the 1.0 m-grid, even qualitative differences are found in comparison to the 0.5 m-grid. This suggests the 0.5 m-grid with a cell is appropriate for reliably predicting flows around an isolated cubic building. The details of grid sensitivity study can also be found in Jeong (2017). The 0.5 m-grid was modeled around H₂S outlet in this study.

Fig. 2.

Mean wind velocities U (m/s) on different grids at X/H=-1.0 (a) and X/H=3.0 (b), Δ is mesh size, H (=10 m) is building height.

3. 2 Comparison with Wind Tunnel Results

Table 2 compares the reattachment lengths on the roof (X_R) and behind the building (X_F). The X_R values obtained by computations show good agreement with the experimental values. However, X_F is greatly overestimated in the simulation result. These results are similar to those shown in a previous study using Realizable k-ε model (Tominaga and Stathopoulos, 2010).

Table 2.

Comparison of reattachment lengths on roof and behind cube of stability D.

Fig. 3 shows profiles of streamwise velocities on the roof and behind the cube at the centerline obtained by the SST k-ω model and the experiment. In this figure, the experimental data for the flows around the cube without vent emission provided by Castro and Robins (1977) is compared with the CFD result for the validation. The differences between the velocity distributions of the experiment and present CFD were rather small.

Fig. 3.

Comparison of vertical distribution of streamwise velocity on (a) roof center (X/H=0.5) and (b) behind cube (X/H=1.5) at centerline.

Fig. 4 shows the contours of the time-averaged dimensionless concentration, K, on the roof and wall surfaces obtained from the SST k-ω model and the experiment by Li and Meroney (1983). In this study, the experimental and the numerical results presented as dimensionless concentration were as follows:

Fig. 4.

Distribution of dimensionless concentration K on roof and wall surfaces.

where C is the mass fraction of the tracer gas, L the size of the cube (m), U_H is mean wind speed at roof height (ms^-1) and Q_e the emission rate of the pollutant (m³ s^-1).

The concentration of the high concentration region (K>50) upwind of the vent in the experiment was larger than those in SST k-ω on the roof surface. The results of K in SST k-ω model showed a spread smaller than those of the experiment in the downstream direction. However, the concentrations are widely expanded in the horizontal direction in SST k-ω. The general distribution of K given by SST k-ω at the side and leeward wall surfaces is similar to that of the experiment, although the SST k-ω result tends to be a little diffusive. The concentrations at the side wall in the SST k-ω are mainly coming from the roof, as in the experiment.

3. 3 Effect of Atmospheric Stability on Flow Field

In order to evaluate the effect of the atmospheric stability on the flow field behind the building, the x-velocity and streamline of centerline under various atmospheric stabilities are shown in Fig. 5. Generally, reattachment length behind the building was increased in stable conditions compared to in unstable conditions. This means that the mixing effect near the cube under unstable conditions was stronger than in stable conditions.

Fig. 5.

Velocity and streamline at Z=0 m for six stability conditions

To compare the rooftop cavity, Fig. 6 shows velocity and streamline at Z=0 m to illustrate the reattachment distance for various stability conditions. The result showed that the unstable condition (stability A, B, C) decreased the length of the rooftop cavity. This is because the vertical velocity is increased and the longitudinal velocity decreased around the rooftop of the building. The stable condition (Stability E, F) increased the horizontal velocity in the lower part around the building and decreased the vertical velocity around the building.

Fig. 6.

Velocity and streamline around building top of centerline (Z=0 m).

In order to compare the temperature distribution around a building for six stability conditions, temperature contours around a building are illustrated in Fig. 7. This comparative view of results illustrates typical temperature profiles in the vertical direction for all boundary layers. However, in strong stable and unstable conditions, the near surface temperature along x direction was slightly changed from inlet conditions. The reason for such results is deemed to be partly because of the use of zero heat flux at the ground surface.

Fig. 7.

Temperature distribution and streamlines around building.

Fig. 8 illustrates the velocity and streamline at Y=1 m for six stability conditions. Stratification increases the horizontal cavity length and width near the surface. However, unstable stratification decreases horizontal cavity length and width near the surface. This means a stable stability condition will increase horizontal pollutant dispersion compared to an unstable stability condition.

Fig. 8.

Velocity and streamline at Y=1 m for six stability conditions.

3. 4 Effect of Atmospheric Stability on Concentration Field

Fig. 9 show mass fraction contours for various stability conditions to compare the plume length. As shown in Fig. 9, stratification increases the plume length and results in a larger pollution impact region at leeward facet. However, more pollutants move to the side facet of the building in an unstable stability condition. Although stable stratification suppresses vertical motion and reduces vertical spreading, it also increases lateral motion. As a result, stable stability increases the lateral spread of the plume on the leeward surface. The concentration levels close to the ground surface were higher under stable conditions than under unstable and neutral conditions. This is because the lateral flow carries the contaminant away from the building.

Fig. 9.

Mass fraction contour of surface for various stability conditions.

Fig. 10 shows the plume contours at Z=0 m for six stability conditions. Unstable stratification decreases horizontal plume length but increases vertical plume height. However, stratification increases the horizontal plume length and decreases vertical plume height. This means a stable stability condition will increase horizontal pollutant dispersion compared to an unstable stability condition.

Fig. 10.

Centerline plume contours for various stability conditions.