Large-eddy simulation of plume dispersion within regular arrays of cubic buildings

Abstract. There is a potential problem that hazardous and flammable materials are accidentally or intentionally released within populated urban areas. For the assessment of human health hazard from toxic substances, the existence of high concentration peaks in a plume should be considered. For the safety analysis of flammable gas, certain critical threshold levels should be evaluated. Therefore, in such a situation, not only average levels but also instantaneous magnitudes of concentration should be accurately predicted. In this study, we perform Large-Eddy Simulation (LES) of plume dispersion within regular arrays of cubic buildings with large obstacle densities and investigate the influence of the building arrangement on the characteristics of mean and fluctuating concentrations.


Introduction
An accurate analysis of plume dispersion is important for emergency responses against accidental or intentional release of hazardous and flammable materials within populated urban areas.For the assessment of human health hazard or the safety analysis of the hazardous gas, not only mean but also fluctuating concentrations should be estimated, considering the effects of individual buildings.Therefore, the dispersion characteristics of a plume through obstacle arrays have been examined mainly by field and wind tunnel experiments.For example, Davidson et al. (1996) investigated the influence of building arrays on plume dispersion by wind tunnel experiments.Bezpalcova and Ohba (2008) conducted wind tunnel experiments of plume dispersion within various building arrays and investigated the effects of the building arrangement and obstacle density on the characteristics of mean and root mean square (RMS) concentrations.
In this study, we perform numerical simulations of plume dispersion within a regular array of cubic buildings as idealized urban canopy by Large-Eddy Simulation (LES) that can give detailed information on turbulent flow and concentration fields.The objective of this study is to perform LES of plume dispersion within building arrays with large obstacle densities, which corresponds to densely built-up urban areas Correspondence to: H. Nakayama (nakayama.hiromasa@jaea.go.jp) and investigate the distribution patterns of concentrations and the characteristics of peak concentration within the building array.

Numerical model
The basic equations for the LES model are the spatially filtered continuity equation, Navier-Stokes equation and the transport equation for concentration.The subgrid-scale (SGS) Reynolds stress is parameterized by using the standard Smagorinsky model (Smagorinsky, 1963), where the Smagorinsky constant is set to 0.1 for estimating the eddy viscosity (Murakami et al., 1987).The subgrid-scale scalar flux is also parameterized by an eddy viscosity model and the turbulent Schmidt number is set to 0.5.
The coupling algorithm of the velocity and pressure fields is based on the Marker and Cell (MAC) method (Harlow and Welch, 1965) with the second-order Adams-Bashforth scheme for time integration.The Poisson equation is solved by the Successive Over-Relaxation (SOR) method which is an iterative method for solving a Poisson equation for pressure.For the spatial discretization in the governing equation of the flow field, a second-order accurate central difference is used.For the dispersion field, Cubic Interpolated Pseudoparticle (CIP) method proposed by Takewaki et al. (1985) is used for the advection term.CIP is a very stable scheme that can solve generalized hyperbolic equations in space.For diffusion term, a second-order accurate central difference Published by Copernicus Publications.method is used.The time step interval ∆tU ∞ /H is 0.005 (∆t: time step).The maximum CFL (Courant-Friedrich-Levy) number is about 0.15.

Wind tunnel experiments for evaluating the model performance
The experiments were carried out by Bezpalcova and Ohba (2008) in the Boundary Layer Wind Tunnel at Wind Engineering Center of Tokyo Polytechnic University, Japan.The experimental set-up consists of buildings with dimensions: 70 mm (width), 70 mm (length), and 70 mm (height).In this paper, obstacle density λ f is defined as the ratio of the total floor projection area of buildings to the plan area of the study site.Buildings are arranged in the regularly square array with λ f = 0.25 and 0.33.There are 18 × 9 and 20 × 9 building arrays with λ f = 0.25 and 0.33, respectively.The ground-level point source is located at the center just behind a building of the 8th row and the 5th column, and the 9th row and the 5th column of the arrays in cases of λ f = 0.25 and 0.33, respectively.Here, the rows are numbered in increasing order in the streamwise direction from the leading edge of the array and the columns are numbered in increasing order in the spanwise direction.In their experiment, the lower part of the neutral atmospheric boundary layer is simulated by vortex generators set up at the wind tunnel section and roughness blocks as shown in Fig. 1.The scale of the modeled boundary layer is 1:400, i.e. the boundary layer height corresponds to 120 m in the full scale.The mean wind velocity vertical profile of approach flow can be approximated by a power law exponent of 0.25.Wind velocity was measured by Thermoanemometry using a split-fibre probe.The uncertainties of flow measurement were 5% for both mean and RMS quantities.Concentration is measured using a fast-response flame ionization detector.The uncertainties of concentration measurement were 9% and 17% for mean and RMS quantities, respectively.In this wind tunnel experiment, the building Reynolds numbers based on the cubical building height and wind speed at the building height is about 14 000.In this study, to evaluate the model performance, we compare our LES results with these wind tunnel experimental data.

Computational settings
Figure 2 shows a schematic illustration of the numerical model.Two computational domains are set up: The main region for a simulation of plume dispersion within a building array and the driver region for generating a spatiallydeveloping turbulent boundary layer flow.First, a thick turbulent boundary layer flow is generated by incorporating the inflow turbulence generation method of Kataoka and Mizuno (2002) into an upstream part of the driver region and, then, a  In the driver region, the Sommerfeld radiation condition (Gresho, 1992) is applied at the exit, a free-slip condition for streamwise and spanwise velocity components is imposed and vertical velocity component is 0 at the top.A periodic condition is imposed at the side and a non-slip condition for each velocity component is imposed at the ground surface.The size and the number of grid points for the driver region is 13.8δ L ×3.8δ L ×5.0δ L (δ L : the scale of the modeled boundary layer) and 460×250×100 in streamwise, spanwise and vertical directions, respectively.A tripping fence and each roughness block set up in the driver region are resolved by 7 × 250 × 24 and 3 × 6 × 12 grids in streamwise, spanwise and vertical directions, respectively.The Van Driest damping function (Van Driest, 1956) is incorporated to account for near-wall effects and the resolution of a grid above the ground surface is set to 2.0.Building effects are represented by the feedback forcing method proposed by Goldstein et al. (1993).The main idea of this method is to apply the external force inside the body.
In the main region, there are 25 × 8 and 28 × 9 obstacle arrays with λ f = 0.25 and 0.33, respectively.The ground-level point source is located just behind a building of the 8th row and the 4th column, and the 9th row and the 5th column of  the arrays in cases of λ f = 0.25 and 0.33, respectively.Each building of the array is resolved by 16 × 16 × 24 grids in the streamwise, spanwise and vertical directions, respectively.At the inlet of the main region, the inflow turbulence data obtained near the exit of the driver region is imposed.The other boundary conditions in a flow field are the same as those in the driver region but the damping function to account for near-wall effects is not incorporated.In a concentration field, zero gradient is imposed at all the boundaries (Shi et al., 2008).Assuming that the location of a plume source point in the wind tunnel experiment is x/H = 0.0, y/H = 0.0 and z/H = 0.0 (H: a building height), that of a plume source point in this LES model is x/H = 0.0, y/H = −0.03and z/H = 0.0.Because the number of grid points for individual cubic building is even number, the plume source position in the y/H coordinate is slightly different from that in the experimental condition.The size and the number of grid points for the main region are 18.0δ L × 3.8δ L × 5.0δ L and 1000 × 250 × 100 in streamwise, spanwise and vertical directions, respectively.The lengths of the domain in front of the first row and behind the last row are both 3.0δ L .The grid resolution above the ground surface is the same as the one in the driver region but the Van Driest damping function is not used in the main region.
The length of the simulation run to calculate the time averaged values of velocity and concentration T U ∞ /H (T : averaging time) is 500.The length of the simulation run before releasing the scalar is T U ∞ /H is 250.In the present LESs, the building Reynolds number is almost 5000.

Approach flow
Figure 3 compares the LES results of turbulence characteristics of approach flow with the wind tunnel experimental data of Bezpalcova and Ohba (2008) and the recommended data of Engineering Science Data Unit 85020 (ESDU 85020, 1985).ESDU 85020 provides comprehensive turbulence characteristics of neutrally stratified atmospheric boundary layer based on independent experimental measurements ranging from the ground surface to 300 m.ESDU 85020 recommends vertical profiles of turbulence intensities for each wind component and their relationship in dependence on surface roughness.δ L corresponds to 4.    Counihan (1975), it is shown that the average height of the constant shear stress layer is 100 m.The LES data lies within this range shown by Counihan (1975).Figure 4 shows the power spectrum of the approach flow of (a) the iment and (b) LES.f E( f ) and L ux indicate the frequency, the longitudinal velocity spectra and the integral length scale, respectively.The each power spectra obtained by wind tunnel experiment is consistent with the Karman type.Although LES power spectra rapidly decrease in higher frequency side f E( f )/σ u > 2, the LES data are found to show good agreement with the Karman type except the high frequency side.
The LES approach flow corresponds to a neutral atmospheric boundary layer based on comparison with the ESDU 85020 recommended data.Although some of the turbulence characteristics by LES are quantitatively different from those by the experiment, they both reasonable well model the neutral boundary layer above rough surface and can compared taking in account their differences.

Dispersion characteristics
Figure 5 shows instantaneous plume dispersion fields in case of λ f = 0.25 at times t * (= tU ∞ /H) = 15.0,54.0 and 117.0 after the plume release.The yellow areas on iso-surface indicate 0.01% of initial concentration.It shows that a portion of the plume is moved upwards by the rising airflow behind the upstream building at first, and then the plume is transported in the streamwise direction with being entrained into the street canyon normal to the wind direction.After enough time passing, the plume is found to be transported and dispersed within and above a building array.
Figures 6 and 7 compare the LES results with the wind tunnel experimental data (Bezpalcova and Ohba, 2008) of the spanwise profiles of mean (C ave ) and RMS (C RMS ) concentrations at a height of 0.29H at the 4th and 5th row behind the source in λ f = 0.25 and 0.33, respectively.The mean and RMS concentrations are normalized by wind velocity at the building height (U H ), the building height and the source strength (Q).The experimental data are shown with the error bars.In both cases of λ f = 0.25 and 0.33, the spanwise spread of the plume of the wind tunnel experiments by Bezpalcova and Ohba (2008) is enhanced by the influence of buildings and high concentration region is formed in the range −1.0 < y/H < 1.0.The mean concentration decreases towards the plume edge.Although LES data overpredict in both cases slightly, the tendency such as the formation of the high concentration region around y/H = 0.0 and the decrease the plume edge is the same as for the experimental data.The RMS concentration profile of the experiment shows the local minimum at y/H = 0.0 and the local maximum around y/H = −1.0 and 1.0.Although LES data overpredict around y/H = −1.0slightly, the shape of the RMS concentration profile of LES is the same as for the experiment.From these results, the mean and RMS concentrations of LES are found to be generally similar in magnitude to that of the experiment.Therefore, it is considered that our LES model for plume dispersion within a building array gives satisfactory results.Both concentration characteristic mean and RMS show an asymmetric pattern for LES data.The maximum value can be found at the left hand side (y/H < 0).This is due to even number of the computational cells and asymmetric placement of the source described in Sect.3.2.
Figure 8 shows the streamwise variation of mean and RMS concentrations at y/H = 0.0 at a height of 0.29H.At the shorter distances from the point source, x/H < 1.0, plume dispersion is enhanced by each building.Therefore, the mean concentration in λ f = 0.33 becomes smaller than that in λ f = 0.25.At the position located away from the point source, x/H > 1.0, the magnitude of the decrease with downwind distance becomes small and these data become quite similar due to the sheltering effect by the building array.Macdonald et al. (1997) investigated the influence of obstacle density on mean concentration of a plume by the field experiments.They mentioned that lateral dispersion is enhanced  in the denser arrays for short distances from a releasing point but is generally similar to the open-terrain case for larger distances due to the sheltering effect by the array.The tendency to decrease with downwind distance depending on obstacle density is similar the field experimental study of Macdonald et al. (1997).
The RMS concentration in λ f = 0.33 also becomes smaller than that in λ f = 0.25 at x/H < 1.0 due to the smoothing of concentration fluctuations by the smaller turbulent eddy with denser arrays.At x/H > 1.0, these data are quite similar due to the homogeneous mean concentration field in both cases.

Characteristics of the peak concentration
In case of accidental or intentional release of toxic or flammable gases into the atmosphere, it is important to estimate not only the mean but also the instantaneous high concentrations.In this section, we first investigate time series of concentration fluctuation and then discuss the characteristics of the peak concentrations.Figures 9 and 10 show time series of concentration fluctuation at the central street canyon and crossing section at the 4th and 5th row behind source location in cases of λ f = 0.25 and 0.33, respectively.In case of λ f = 0.25, concentrations fluctuate smoothly and continuously at the position of central street canyon, while instantaneous high concentrations which exceed the average level frequently occur at the crossing section.Mavroidis and Griffiths (2001) examined time series of concentration fluctuation under different conditions, such as in open-terrain, behind an isolated cube and within a building array by the wind tunnel experiments.It was shown that, particularly, higher concentration peaks in a gap between two cubes occur much larger than those behind a cube within a building array.These patterns of concentration fluctuation inside and outside of the cavity region of a building are similar to the wind tunnel experiment by Mavroidis and Griffiths.On the other hand, in case of λ f = 0.33, concentrations are found to fluctuate around the average level both at the central street canyon and crossing section.
Figure 11 shows probability distribution functions (1 − p(c)) of concentration fluctuation at the central street canyon and crossing section.The probability distribution functions in cases of λ f = 0.25 and 0.33 are found to be almost the same at the central street canyon among the cases.On the other hand, at the crossing section, instantaneous high concentrations in case of λ f = 0.25 occur much more frequently than those in case of λ f = 0.33.Furthermore, we evaluate the peak value c 99 defined as the values determined from 1 − p(c) = 0.99 in the wind tunnel experiment and LES.The peak concentration ratios (c 99 /C ave ) of LES at the central street canyon are 2.0 and 2.2 in λ f = 0.25 and 0.33, while those of the experiment are 2.0 and 1.9 in λ f = 0.25 and 0.33.At the crossing section, the peak ratios of LES are 4.5 and 2.7 in λ f = 0.25 and 0.33, while those of the experiment are 2.6 and 2.2 in λ f = 0.25 and 0.33.Although LES data at the central street canyon are in good agreement with the experimental data, those at the crossing position are overestimated.
Focusing on the obstacle density effect on the peak ratios, it is found from the experimental data that the peak at the crossing position becomes much larger than that at the central street canyon in λ f = 0.25 while the peak ratios are similar in both locations in λ f = 0.33.This tendency is similar to the results of our LES model.From these results, it is obvious that the peak concentration ratios show highly different values depending on the locations and obstacle density.This fact indicates that, for the assessment of human health hazard or the safety analysis of the hazardous gas within urban areas, it is important to locally evaluate the peak concentrations considering obstacle morphology.

Conclusions
In this study, we perform LES of plume dispersion within building arrays with large obstacle densities which corresponds to densely built-up urban areas and investigate mean and fluctuating concentrations.The obtained results are as follows: 1.The approach flow is generated by incorporating existing inflow turbulence generation method into an upstream small part of the driver region with a tripping fence and roughness blocks.By this turbulence generation method, the approach flow corresponding to a neutral atmospheric boundary layer is considered to be obtained.
2. When compared to the experimental results of Bezpalcova and Ohba (2008), the spanwise profiles of mean and RMS concentrations are generally similar in the magnitude to the experimental data.Therefore, it is considered that our LES model for plume dispersion within a building array gives satisfactory results and can be used for deeper and finer investigation of the flow and concentration field within the array.
3. The influence of obstacle density on mean concentration is very large for x/H < 1.0 due to the enhancement of plume spreads by each building.For x/H > 1.0, its influence becomes small and these data are quite similar due to the sheltering effect by the building array.RMS concentration becomes smaller with denser arrays for x/H < 1.0.However, these data are quite similar due to the homogeneous mean concentration field in both cases for x/H > 1.0.
4. Although LES peak concentration ratios at the central street canyon are in good agreement with the experimental data, those at the crossing position are overestimated.However, the patterns of the peak ratios depending on the locations and obstacle density are similar to the wind tunnel experiment.Focusing on the obstacle density effect on the peak ratios, it is found that those show different values depending on the locations and obstacle density.
These results imply that, for the assessment of human health hazard or the safety analysis of the hazardous gas within urban areas, it is important to locally evaluate the peak concentrations considering urban morphology.
Although the comparison of the LES data with the experimental data is not sufficiently discussed, we attempted to examine the influence of obstacle density on the spatial distribution of concentrations and the peak concentration characteristics.In order to quantitatively evaluate the obstacle density effects on the dispersion characteristics, the prediction accuracy of our LES model should be further investigated and improved.
Figure 1.The wind tunnel set-up.7

Figure 1 .
Figure 1.The wind tunnel set-up.

Figure 4 .
Figure3compares the LES results of turbulence characteristics of approach flow with the wind tunnel experimental data ofBezpalcova and Ohba (2008) and the recommended data of Engineering Science Data Unit 85020(ESDU 85020, 1985).ESDU 85020 provides comprehensive turbulence characteristics of neutrally stratified atmospheric boundary layer based on independent experimental measurements ranging from the ground surface to 300 m.ESDU 85020 recommends vertical profiles of turbulence intensities for each wind component and their relationship in dependence on surface roughness.δ L corresponds to 4.0H and is assumed to be 120 m in the full scale condition.The experimental data are shown with the error bars described in Sect.3.1.The profile of the mean wind velocity of LES is found to fit the experimental profile of 0.25 power law.LES

Figure 4 .
Figure 4. Power spectrum of approach flow obtained from (a) wind tunnel experiment and (b) LES.