- Email: [email protected]
Prediction of four concentration moments of an airborne material released from a point source in an urban environment
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
Contents lists available at ScienceDirect
Journal of Wind Engineering & Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia
Prediction of four concentration moments of an airborne material released from a point source in an urban environment George C. Efthimiou Environmental Research Laboratory, INRASTES, NCSR Demokritos, Patriarchou Grigoriou & Neapoleos Str., 15310, Aghia Paraskevi, Greece
A R T I C L E I N F O
A B S T R A C T
Keywords: RANS Beta distribution MUST Wind tunnel Skewness Kurtosis
Іn case of the dispersion of an airborne material from a point source in an urban environment the reliable prediction of the concentration statistical distribution by a numerical dispersion model presupposes the capability of the model to predict at least four statistical moments (mean, variance, skewness and kurtosis). In the present study, the beta distribution, the selection of which is justified based on a previous study, is incorporated in the Reynolds Averaged Navier Stokes (RANS) methodology. The shape parameters of the beta distribution are calculated using the numerical results of the mean, variance and maximum concentration. The latest is calculated through a deterministic model which uses also the numerical results of the mean and variance concentration as well as a hydrodynamic time scale. The validation of the new hybrid model “RANS-beta” is performed using the experimental dataset of the MUST wind tunnel experiment. The performance of the “RANS-beta” model for the skewness is very good (FAC2 ¼ 0.811) while for the kurtosis it is acceptable (FAC2 ¼ 0.557). The discrepancies are observed mainly at the edges of the plume. Future research will be focused on the optimization of the mean flow field, turbulent quantities (Reynolds stresses, turbulent kinetic energy) and on new parameterizations for the turbulent diffusion term.
1. Introduction Increasing urbanization and the resulting localization of human living and industrial production within relatively small areas lead to very demanding standards for human safety and security in urban environments. In this context, there is a necessity to be able to provide reliable decision support to minimize or even eliminate the potential health (or even life) hazard from accidental or deliberate atmospheric releases of harmful toxic, flammable, radioactive or bacteriological substances that might happen within complex urban and industrial areas. Accidental or deliberate airborne hazardous releases might occur in urban or industrial environments (e.g. COST Action ES1006, http://www.elizas.eu/) and they are threatening for human health and life, as people can be exposed to released substances through, e.g., inhalation, skin or eye contact. Additionally, everyday pollution sources like cars and chimneys, cause a permanent health risk for city dwellers. The chaotic nature of atmospheric turbulence increases the uncertainties for a realistic prediction of human exposure. For this reason, the use of statistics is imperative for this research field. Thus, the exposure is inherently stochastic due to turbulence variability. The problem becomes even more complex when dispersion occurs over complex
terrain such as the urban environment. If the concentration is measured for relatively short time periods at a specific location downwind of the source, then it can be observed that the range of the concentration values is wide. If a human is at the same location for this time period, then it is crucial to predict the probability – or the percentage of time - that the concentration will be higher than or equal to a threshold that could cause adverse health effects. In other words, the ability to predict the upper tail of the concentration probability distribution downstream of the source would be an important asset contributing to an efficient emergency management. Many researchers in the literature have proposed various statistical distributions to describe concentration fluctuations in the atmospheric surface layer such as exponential (Lewis and Chatwin, 1997), gamma (Nironi et al., 2015), beta (Bartzis et al., 2015), lognormal (Gailis et al., 2007), generalized Pareto (Munro et al., 2003) etc. In the last years, special interest is given by researchers about the gamma and the beta distribution. In Efthimiou et al. (2017b) we proved that the beta distribution predicts better the concentration fluctuations of an airborne material released from a point source in an urban-like area (MUST field experiment). The question that arises is that if this distribution can be incorporated in numerical dispersion models. The advantage of such an
E-mail address: [email protected]. https://doi.org/10.1016/j.jweia.2018.11.032 Received 29 May 2018; Received in revised form 3 October 2018; Accepted 24 November 2018 Available online 3 December 2018 0167-6105/© 2018 Elsevier Ltd. All rights reserved.
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
turbulence closure, Launder and Spalding, 1974) are calculated in the computational grid. Also, standard wall functions have been used in the present study for rough walls.
effort is that simple models that can predict two concentration moments (mean and variance) and a turbulent timescale will have the capability to predict the peak concentration. Then, the beta distribution can be used to predict all the concentration percentiles, necessary for probabilistic health impact assessments. The two higher moments (skewness and kurtosis) are used to validate the assumption of a beta distribution. Furthermore, in the literature, many researchers have selected the Computational Fluid Dynamics for the study of the flow field and the pollutant dispersion in urban areas. The modeling of the atmospheric turbulence is the main difficulty and the basic element of differentiation of research works. The most advanced method for this type of studies is the Direct Numerical Simulation (e.g. Goulart et al., 2018; Coceal et al., 2014). This method can be used e.g. for the understanding of the physical mechanisms of turbulence, the study of coherent structures etc. There is also the Large Eddy Simulation (e.g. Tolias et al., 2018; Jiang and Yoshie, 2018; Hertwig et al., 2018; Li et al., 2018; Aristodemou et al., 2018; Cheng and Port e-Agel, 2016; Nakayama et al., 2014; Liu et al., 2012) which has many applications e.g. it can be used also for basic understanding of turbulence and it can be used for emergency response. Concerning the latest, reliable conclusions are needed (e.g. Koutsourakis et al., 2012) in short times (e.g. less than 10s) and the selection of LES and the simple Reynolds Averaged Navier Stokes (RANS) model can be considered attractive. In this case, “smart” algorithms could extract the information from an LES or RANS database. The present author works this period on emergency preparedness and response tools in case of source identification (e.g. Efthimiou et al., 2018). This is the main reason why RANS is examined in this study. However, RANS can predict the first (mean) and the second (variance) moments through numerical solution of the corresponding transport equations. The incorporation of the beta distribution in the RANS methodology will allow the calculation of two additional moments (skewness and kurtosis). The objective of the present paper is to perform the incorporation of the beta distribution in the RANS methodology and to examine the validity of the hybrid “RANS-beta” model. A similar work has been performed in Efthimiou et al. (2016) where the hybrid “RANS-gamma" model was validated with wind tunnel measurements of the Michelstadt and CUTE experiments. In the present paper, the validation is performed with the MUST wind tunnel experiment. In this work the in-depth investigation of sophisticated numerical schemes, turbulence models etc. is not performed but the interest is focused on the operationality of the method in order to be used for emergency response in the future. Optimization of the methodology will be performed in a future study. In Chapter 2 the methodology is described i.e. the incorporation of the beta distribution in the RANS methodology. In Chapter 3 the experiment is described, in Chapter 4 the numerical simulations and in Chapter 5 the results are presented. The conclusions are presented in Chapter 6.
2.2. Numerical solution of the concentration mean and variance equations In this step, the partial differential equations (pde) of concentration mean, C; and variance, C'2 , are solved numerically in the computational grid using the velocities and the turbulent parameters of step 1. The pde of the mean concentration in Cartesian coordinates is usually expressed as follows:
∂ ∂ ∂ ∂C ρC þ ρui C ¼ ρD ρu0 i C0 ∂xi ∂t ∂xi ∂xi
þ SC
(1)
In Eq. (1), C is the mean concentration, t is the time (s), ρ is the density of the air-pollutant mixture (kg m3), xi is the coordinate (m) along the i0 direction,ui is the mean velocity (m s1) in the i-direction, u i is the ve1 locity fluctuation (m s ) in the i-direction, C’ is the concentration fluctuation, D is the molecular mass diffusivity (m2 s1) and SC is the pollutant sources/sinks. The overbar denotes Reynolds average, while repeated indices in a term imply summation over all possible values of the index, according to the Einstein summation convention. To close Eq. (1) the turbulent diffusion term u0 i C0 needs to be modeled. This is done in the present study through the gradient transfer assumption, due to its wide use in the modeling of engineering flows and its computational simplicity. So, following Csanady (1967) the turbulent diffusion term is modeled using the eddy mass diffusivityK C : u0 i C 0 ¼ K C
∂C K m ∂C ¼ ∂xi σ C ∂xi
(2)
In Eq. (2) the eddy mass diffusivity is calculated by dividing the eddy viscosity K m by the turbulent Schmidt number σ C. The pde for the mean square fluctuation of concentration is expressed as follows:
0
1
∂ ∂ ∂C ∂ @ ∂C 0 2 2 ρC0 2 þ ρui C0 2 ¼ 2ρu0 i C0 þ ρD ρu0 i C 0 A ∂xi ∂xi ∂xi ∂t ∂xi 2ρD
∂C 0 ∂C 0 ∂xi ∂xi
(3)
where C0 2 is the concentration variance. Eq. (3) accounts for all the physical phenomena related to the advection, production, diffusion and dissipation of this quantity. To close Eq. (3), two terms need to be modeled: the turbulent diffusion term and the dissipation term.
2. Methodology
The turbulent diffusion term u0 C0 2 is modeled by using again the
The present methodology consists of four steps:
0
gradient transport hypothesis with an eddy diffusion coefficient K C that is assumed to be equal to the eddy mass diffusivity K C (Csanady, 1967; Donaldson, 1973; Sykes et al., 1984):
1. Numerical simulation of the hydrodynamic problem. 2. Numerical solution of the concentration mean and variance equations using the results of step 1. 3. Calculation of maximum concentration using the results of steps 1 and 2. 4. Calculation of concentration skewness and kurtosis using the results of steps 2 and 3.
2
u0 i C 0 ¼ K C
0
∂C 0 2 ∂C 0 2 K m ∂C 0 2 ¼ K C ¼ ∂xi ∂xi σ C ∂xi
(4)
The dissipation term (last term) of Eq. (3) is modeled as the ratio of the mean square fluctuations of the mass fraction to the dissipation time
Analysis of each step follows.
0
2
scale T C of C0 . Following Fackrell and Robins (1982), the dissipation time scale is expressed as the ratio of a length scale LC to the velocity scale k1/2, where k is the turbulent kinetic energy:
2.1. Numerical simulation of the hydrodynamic problem In this step, the 3-dimensional, steady-state fields of the wind velocity components u, v, w, the turbulence kinetic energy, k, and the turbulence kinetic energy dissipation rate, ε, (in case of using the standard k-ε
D
248
1=2 0 2 ∂C 0 ∂C 0 C 0 2 C 0 k ¼ 0 ¼CC ∂xi ∂xi T C LC
(5)
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
0
a¼
where C C is a constant equal to 0.5 and LC is an effective length scale of turbulence which for the selected turbulence model (standard kε) is 1:5 given by LC ¼ C0:75 μ k =ðκ εÞ, where κ ¼ 0.4 is the von Karman constant and Cμ ¼ 0.09. In the present study the eddy viscosity K m (appearing in Eqs. (2) and (4)) is calculated by the standard k-ε model, which solves transport equations for the turbulent kinetic energy k and the dissipation rate ε.
1 η 1 1þη I
ξ ¼ ηa
η¼
Cmax ðΔτÞ C C
(13) (14)
(15)
It is reminded that the variable Cmax(Δτ) is calculated in step 3 while C in step 2.
2.3. Calculation of maximum concentration 3. The MUST wind tunnel experiment
In this step, the maximum concentration at an interval time Δτ is calculated by the deterministic model Bartzis et al. (2015). The equation is: " Cmax ðΔτÞ ¼ 1 þ 9I
Δτ TC
The methodology has been validated against data of the MUST wind tunnel experiment (Bezpalcova and Harms, 2005) which have been scaled up for the conditions of the corresponding field experiment (Yee and Biltoft, 2004). In MUST experiment the obstacles were arranged in 12 rows, each consisting of 10 obstacles. The obstacles were nearly identical and had average length, width and height 12.2 m 2.42 m 2.54 m respectively. The contaminant's concentration has been measured by a 256-detectors array arranged along obstacle rows in the part of the domain covered by the plume (Fig. 1). All detectors were placed at the same height equal to 1.275 m. Non-zero concentrations have been measured by nearly all (244 out of 256) detectors. Here it is noted that only the measurement data that were available to the authors (through data base of COST Action 732) are taken into account. The link to the database of COST Action 732 is: https://mi-pub.cen.uni-hamburg.de/index.php?id¼484. The wind flow was characterized by neutral stratification, wind speed at the roof level Uref ¼ 8 m/s, (field scale) and wind direction 45 (Fig. 1) in the experimental coordinate system. The contaminant originated from a point source located at the ground level. The volume flow rate of gas at the source was: 3.3 106 m3s1. Wind tunnel measurements of u and v (and w for intercomparison purposes with Tominaga et al., 2013) were used also to validate the numerical results of the hydrodynamic problem. Two groups of sensors have been used which correspond to measurements in a coarse network and in a fine network (Fig. 2). The height of the velocity and concentration sensors is the same and equal to 1.275 m.
0:3 # C
(6)
where TC is the turbulence integral time scale, which under fully mixed conditions it is given by the relationship: TC ¼
Cμ k 3=2 ; Cμ ¼ 0:09 Uε
(7)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and Uis the mean wind speed: U ¼ u2 þ v2 þ w2 (ms1). The fluctuation intensity I is defined as: I¼
C0
2
C
2
(8)
It is reminded that the variables u, v, w, k and ε are calculated in step 1 while C and C'2 in step 2. Δτ is the exposure time considered. The relevant duration of the exposure time Δτ is selected either on health impact considerations or on the hazardous substance release duration or on real conditions (e.g., the time an individual is expected to stay at the particular location). Especially near the emission source there might be cases where the concentrations are very high and related time scales are as low as breathing time (i.e., a few seconds). On the other hand very short releases in the order of a few seconds up to a few minutes are common scenarios for the above mentioned accidental or deliberate release events. In this work the averaging time interval Δτ is considered equal to the measurement time interval of the instrument.
2.4. Calculation of concentration skewness and kurtosis The general definitions of the skewness (Sk) and kurtosis (Ku) are: Sk ¼
Ku ¼
C0 3 C0 2 C0 C0
1:5
(9)
4
2
2
(10)
In case of the beta distribution Sk and Ku are calculated based on the following theoretical equations (Efthimiou et al., 2017b): Sk ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ðξ aÞ a þ ξ þ 1 pffiffiffiffiffi ða þ ξ þ 2Þ aξ
(11)
Ku ¼
6 ða ξÞ2 ða þ ξ þ 1Þ aξða þ ξ þ 2Þ þ3 aξða þ ξ þ 2Þða þ ξ þ 3Þ
(12)
Fig. 1. The computational domain of the MUST case. The sensors are presented with yellow circles. The source is depicted with a star. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
where 249
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
Fig. 2. The velocity sensors (yellow circles) of the MUST wind tunnel experiment placed at 1.275 m above the ground (left: coarse network, right: fine network). The white rectangles are the obstacles. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
“More details about the experiment can be found in the website of COST Action 732 (https://mi-pub.cen.uni-hamburg.de/index.php? id¼484) and personal communication with the experimentalists.”
should be noticed that in the experiment not the entire bottom was rough but only a part of it. Concerning CFD guidance on the treatment of walls, detailed information can be found in Franke et al. (2007). The fixed vertical profiles at the inlet planes (-x and þy) for u, v, k and ε have been constructed based on the experimental corresponding profiles, while for w a zero value has been set. The outlet boundary conditions at planes þ x and –y are justified, as these planes are located at a sufficient distance (see above) downwind of the last buildings. For turbulence modeling, the standard k–ε model (Launder and Spalding, 1974) was used. The problem was treated as a pseudo-transient case i.e. time integrated towards steady state. The unsteady Reynolds-Averaged NavierStokes equations for total mass, the three components of momentum, the turbulence kinetic energy and its dissipation rate were solved. At the end of the solution, steady-state values were produced i.e. constant (in time) values of the variables at the sensors. The total calculation time (1200s) was set about equal to the travel time from the source point to the outlet of the computational domain of a hypothetical particle moving with the lowest speed observed in the experiment. This time was sufficient to reach a steady-state situation. The time step was automatically adapted according to the desired Courant–Friedrichs–Lewy (CFL) number whose maximum value was set to 2 based on the author's experience for numerical stability (better convergence behaviour of the solver). Nevertheless, as an implicit scheme is adopted and a steady-state solution is sought for, the magnitude of the final time step does not play here an important role. A first order scheme was used for the time derivative. For the calculation of pressure, the linear system that arises from the discretized equations was solved using the preconditioned BiCGstab method. A block-ILU[0] type of preconditioner was used (additive Swartz preconditioner, Saad, 2003).
4. The numerical simulations All computations (wind flow and dispersion) were performed with the CFD code ADREA-HF (Bartzis, 1991; Andronopoulos et al. 1994, 2002; Venetsanos et al., 2010). This is a Eulerian model, solving the unsteady RANS equations. The code uses finite volumes with rectangular cells for the discretization of the Eulerian transport equations. To describe the complex geometry, the volume porosity concept is used with solid surfaces of any orientation allowed to cross the computational cells. The ADREA-HF code has been setup for the simulation in the field scale. 4.1. The wind field simulations The spatial discretization of the computational domain is presented in Table 1. The domain extends horizontally by 17.55 m (6.9 H, H ¼ 2.54 m being the buildings height) upwind of the first buildings and by 52.65 m (20.7 H) downwind of the last buildings. The vertical dimension of the domain is about 8.3 H. The above dimensions conform to the recommendations of COST Action 732 (Franke et al., 2007). The grid is Cartesian. Cubic cells have been selected inside the urban area in order to decrease the discretization errors. Outside the urban area, the grid increases logarithmically by a factor of 1.1. It is well known by CFD practices for microscale modeling (Franke et al., 2007) that the building height should be described by at least ten cells. Thus, in case of MUST, a minimum cell size should be equal to dxmin ¼ dymin ¼ dzmin ¼ 0.25 m. The total number of cells of the present simulation is 21,807,800 which is finer than the one used in Efthimiou et al. (2017c). The boundary conditions for the hydrodynamic variables at each boundary plane or surface of the domain are presented in Table 2. It
Table 2 Boundary conditions for the hydrodynamic variables (u, v, w: velocity components in the x-, y-, z-axis respectively, k: turbulence kinetic energy, ε: turbulence kinetic energy dissipation rate).
Table 1 Dimensions and spatial discretization of the computational domain. Domain dimensions x/ y/z (m)
Total number of cells
Number of cells in each axis
Minimum/maximum cell sizes (m)
x
y
z
dx
dy
dz
243/268/21
21,807,800
740
842
35
0.25/ 5.1
0.25/ 5.1
0.25/ 1.84
250
Plane
Boundary condition
-x þx -y þy Ground Top Building walls
Inlet: fixed vertical profiles for u, v, k, ε, w ¼ 0 Outlet: ∂ϕ=∂x ¼ 0, φ ¼ u, v, w, k, ε Outlet: ∂ϕ=∂y ¼ 0, φ ¼ u, v, w, k, ε Inlet: fixed vertical profiles for u, v, k, ε, w ¼ 0 Standard wall functions, roughness length ¼ 0.0165 m Fixed values for u, v, k, ε, w is calculated from the cell mass balance Standard wall functions, roughness length ¼ 1 105 m
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
experiment. For example, in Tominaga et al. (2013) the HRs of 8 models, for 498 points (xz plane), for the velocity component u was ranged from 0.63 to 0.91 while for the velocity component w was ranged from 0.12 to 0.26. The HR of the same velocity components of the present simulation is inside these ranges (0.75 for u/Uref and 0.18 for w/Uref) indicating the consistency with the results obtained by Tominaga et al. (2013) and by COST 732 group (Schatzmann et al., 2010). The HR of each velocity component at each sensor of Fig. 2 is presented in Fig. 3. The black circle indicates HR ¼ 1 and a good prediction while the red circle indicates HR ¼ 0 and a poor prediction. For the velocity u/Uref, there are a lot of hits everywhere in the urban area. On the other hand, the HR of v/Uref is good only upwind and in the first rows of buildings. It should be noticed that at the selected elevation the influence of the prescribed wall boundary conditions is strong and most likely causing the observed deviations from the experiment. Furthermore, the numerical data were linearly interpolated in order to match the spatial positions of the experimental measurements. This procedure adds further uncertainties to the predictions which are hard to distinguish from errors that originated from the simulation itself. However, the error/uncertainty introduced by interpolation should be viewed part of the
The 1st order upwind scheme has been used for approximation of convective terms. The second order central-difference scheme has been used for approximation of diffusion terms. It should be noticed that, since the grid is finer than the one used in Efthimiou et al. (2017c), where a grid independent numerical solution was attained, the selection of the numerical schemes for the spatial derivatives does not have significant effects on the results. The computations were run in parallel using 20 cores of 5 processors (i.e. 100 processors). ADREA-HF is parallelized for distributed memory architectures using the MPI protocol. 4.2. The pollutant dispersion simulations The transport equations for the concentration mean and variance were solved on the grid simultaneously with the hydrodynamic problem from which the flow and turbulence related variables are used for the advection and diffusion terms. The boundary conditions of the concentration mean (C) and variance (C'2 ) for each boundary plane or surface of the domain are presented in Table 3. The total calculation time was set equal to the one used for the hydrodynamic computations, i.e., the 1200s, which was adequate for obtaining a steady-state solution. For the discretization of the convective term in the transport equations, the upwind scheme was used, while for the diffusion term the central scheme was used. The time step was automatically adapted according to the desired Courant–Friedrichs–Lewy (CFL) number whose maximum value was set to 2. 5. Results and discussion 5.1. Validation of the hydrodynamic problem solution Author has checked in a past paper the profile development across the domain of the same experiment (Efthimiou et al., 2017c). It was noticed that the model presented a very good performance above the buildings, where the Hit Rate (HR) (VDI Guideline 3783/9, 2005) had almost the ideal value (¼1.0) for the velocities u/Uref and v/Uref. Below the building roofs the HR was gradually decreased for all velocity components. u/Uref presented in total better agreement with the experimental data than v/Uref and this is normally the case for the components being aligned with the mean flow direction in x-direction. The reason for this is that the magnitudes of the lateral component, i.e. in y-direction, are very much smaller than the magnitudes of the component in the mean flow direction. In this case the coordinate system corresponds to the wind tunnel. Very close to the ground the velocity v/Uref presented the lowest HRs. The HR of the velocities u/Uref and v/Uref for the measurements below the half height of the buildings was found equal to 0.48 and 0.14 respectively in Efthimiou et al. (2017c). Similar values of the HR were found also in the present simulation which for the specific sensors at the height of 1.275 m (Fig. 2) are equal to 0.55 and 0.16 respectively. It should be noticed that in the literature there are exercises where various models have been intercompared using the MUST wind tunnel Table 3 Boundary conditions for the concentration mean (C) and variance (C'2 ). Plane
Boundary condition
-x
∂C=∂x ¼ 0, ∂C0 2 =∂x ¼ 0
þx
∂C=∂x ¼ 0, ∂C0 2 =∂x ¼ 0
-y
∂C=∂y ¼ 0, ∂C0 2 =∂y ¼ 0
þy
∂C=∂y ¼ 0, ∂C0 2 =∂y ¼ 0
Ground
∂C=∂z ¼ 0, ∂C0 2 =∂z ¼ 0
Top
∂C=∂z ¼ 0, ∂C0 2 =∂z ¼ 0
Building walls
∂C=∂n ¼ 0, ∂C0 2 =∂n ¼ 0, n ¼ x, y, z
Source (mass fraction)
C ¼ 1, ∂C=∂n ¼ 0, ∂C0 2 =∂n ¼ 0, n ¼ x, y, z
Fig. 3. The Hit Rate (HR) per sensor of the velocity components u/Uref (upper) and v/Uref (lower). The black circle indicates HR ¼ 1 and a good prediction while the red circle indicates HR ¼ 0 and a poor prediction. The results are presented for height 1.275 m above the ground for the coarse and fine networks. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) 251
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
5.2. Validation of the concentration statistical moments
simulation (Oberkampf and Barone., 2006). Finally, as has been discussed in Hertwig et al. (2012), unsteady flow effects triggered by the presence of buildings are a further influencing factor for the model performance. Another way to compare the results of the model and of the experiment is the vector plots. Fig. 4 presents horizontal vector plots of the experiment (upper) and the model (lower) at the area of the fine network of sensors (Fig. 1). In this case the coordinate system corresponds to the model and it is not the same as the one used for the calculation of the HRs. The direction of the flow is almost the same but an underprediction of the velocity magnitudes by the model is clear. This underprediction was observed also by author in a past work (Hertwig at al. 2012).
In the present study, the validation of the four concentration statistical moments (mean, variance, skewness and kurtosis) was performed using scatter plots, validation metrics and contour point-by-point plots. The experimental database provided normalized results of the concentration mean and root mean square (RMS) and as a result, the numerical results were normalized in the same way based on the following equation: C* ¼
CH 2 Uref Q
Fig. 4. Horizontal vector plots of the experiment (upper) and the model (lower) for the area covered by the fine network of sensors (Fig. 1). 252
(16)
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
intermittency). The model predicted very well the values near the centerline as many of the points lie on the 1 to 1 line. The poor prediction of the highest skewness values (at the edges of the plume) is an indication that the concentration mean and RMS are inversely proportional to the skewness. The optimization of the turbulent diffusion term will optimize also the performance of the highest skewness values. The performance is also good for the kurtosis. The positive sign reveals that the distributions are leptokurtic. The sensors at the edges of the plume have much fatter tails than the sensors at the plume centerline. Also, the scatter of the points around the ideal line is higher than skewness. Table 4 presents the Normalized Mean Square Error (NMSE), the Fractional Bias (FB) and the Factor of 2 (FAC2) for the four concentration moments. The explanation of these metrics is given in Chang and Hanna (2005). The ideal values are NMSE ¼ 0, FB ¼ 0 and FAC2 ¼ 1. Based on the available information, the experimental confidence is provided only for the concentration mean and RMS and is equal to 0.003. However, this value was not considered in the present study for the calculation of the FAC2. The FAC2 of the skewness (0.811) is very good and better than the other moments supporting the robustness of the methodology. However, the validation metrics can lead to erroneous conclusions about the performance of a model. Even if in Fig. 5 the scatter of the points around the ideal line (1–1) for the kurtosis is not so high in
where C is the predicted concentration (either mean or RMS), H is height of the buildings, Uref is the reference velocity and Q is the source rate. Fig. 5 presents the scatter plots with logarithmic axes of the concentration mean, RMS, skewness and kurtosis. For the concentration mean we observe overprediction for the higher and the medium values and underprediction for the lower values. For the RMS the prediction is better than the mean as more values lie on the 1 to 1 line and we observe also the underestimation for the lower values. There are mainly two parameters that influence the performance of the model: a) the mean flow field and turbulent quantities (Reynolds stresses, turbulent kinetic energy) which are the driving force of the dispersion. The optimization in this case is performed with sophisticated turbulence models (e.g. Bartzis, 2005) and boundary conditions (e.g. Franke et al., 2007) and b) the turbulent diffusion term. It is reminded that this term is based on the Schmidt number the value of which is constant everywhere and this raises questions about the validity of this assumption. Furthermore, the turbulent diffusion term is based on the full mixing hypothesis which means that we assume that the concentration and the wind flow have the same length/time turbulent scales everywhere. Yee et al. (2009) have separated these scales, using also the pollutant travel time, for the calculation of the dissipation time scale of the concentration variance and their results were very good. Beyond concentration variance, Efthimiou and Bartzis (2011) have used the pollutant travel time for the calculation of the concentration autocorrelation time scale. The usefulness of this parameter is high as it can be used in combination with the eddy viscosity and the hydrodynamic autocorrelation time scale for the calculation of the turbulent diffusion term. This is a task that will be performed by author in the future. Concerning the concentration skewness, the positive sign reveals that the right tails are longer. Close to the plume centerline (lower intermittency) the skewness is lower than at the plume edges (higher
Table 4 Validation metrics (NMSE, FAC2 and FB) for the concentration moments.
Mean RMS Skewness Kurtosis
NMSE
FAC2
FB
2.85 4.42 47.00 927.21
0.496 0.627 0.811 0.557
0.267 0.036 0.765 1.701
Fig. 5. Scatter plots of the concentration mean, RMS, skewness and kurtosis. The mean and RMS are normalized values to be in accordance with the experiment. 253
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
comparison with the other moments, its NMSE is very high and this is due to the very high discrepancies of the model from the experiment at the edges of the plume. Similar behaviour is observed also for the skewness. The higher overprediction is observed for the kurtosis and then for the skewness and this is again due to the highest values. The overprediction of the concentration mean (FB ¼ 0.267) is based on many parameters such as the turbulent diffusion term, the underprediction of the velocity field by the model, the turbulence model, the numerical schemes etc. Overall the results are acceptable. The FAC2 of the mean concentration (0.5) is close to the lower values reported by Tominaga et al. (2013) (from 0.55 to 0.88), for 256 points (xy plane). However, in the present calculation of the FAC2 all the experimental concentration values were considered. Finally, in Fig. 6 and Fig. 7, point by point contour plots are presented for the skewness and kurtosis respectively for the experiment (upper) and the model (lower). The range of values for each color is common for both moments. The similarities between the experiment and the model are very good for both moments especially at the main body of the plume
supporting the robustness of the “RANS-beta” model. The discrepancies are observed at the edges of the plume supporting the previous thoughts.
Fig. 6. Contour point-by-point plots of the concentration skewness for the experiment (upper) and the model (lower).
Fig. 7. Contour point-by-point plots of the concentration kurtosis for the experiment (upper) and the model (lower).
6. Conclusions In this work the operationality of the new, hybrid “RANS-beta” model was examined in order to predict four concentration moments (mean, variance, skewness and kurtosis) in case of the dispersion of an airborne material from a point source in the atmospheric surface layer. The computational cost concerns only the numerical simulation of the partial differential equations of concentration mean and variance as the concentration skewness and kurtosis are calculated during the postprocessing of the results. This process, overall, can be considered as an additional advantage for the RANS methodology during an emergency response. Concluding, it seems that, both RANS and LES are viable methods to generate data bases for operational emergency response modeling.
254
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
The validation was focused mainly on the skewness and the kurtosis. The agreement of the model with the MUST wind tunnel experiment for the concentration skewness is high (FAC2 ¼ 0.811) and acceptable for the kurtosis (FAC2 ¼ 0.557). The errors in the prediction of the concentration mean and variance were transferred in the prediction of concentration skewness and kurtosis and this was obvious at the edges of the plumes. Initially, optimization of the mean flow field is necessary as it is the driving force of the dispersion. Also, more sophisticated parameterizations for the turbulent diffusion term are necessary. These could include the relation of this coefficient with the pollutant travel time and the autocorrelation time scales of the concentration and the wind flow. These issues will be examined more carefully by author in a future publication. With better results for the prediction of the flow field and lower concentration moments one might obtain better results for skewness and kurtosis. Optimization of the methodology using e.g. higher order numerical schemes, more sophisticated turbulence models will be performed in a future study. It would be interesting also to incorporate the beta distribution in Large Eddy Simulation and this will be a future task.
Efthimiou, G.C., Bartzis, J.G., 2011. Atmospheric dispersion and individual exposure of hazardous materials. J. Hazard Mater. 188, 375–383. Efthimiou, G.C., Kovalets, I.V., Argyropoulos, C.D., Venetsanos, A., Andronopoulos, S., Kakosimos, K., 2018. Evaluation of an inverse modelling methodology for the prediction of a stationary point pollutant source in complex urban environments. Build. Environ. 143, 107–119. Efthimiou, G.C., Kovalets, I.V., Venetsanos, A., Andronopoulos, S., Argyropoulos, C.D., Kakosimos, K., 2017c. An optimized inverse modelling method for determining the location and strength of a point source releasing airborne material in urban environment. Atmos. Environ. 170, 118–129. Fackrell, J.E., Robins, A.G., 1982. Concentration fluctuations and fluxes in plumes from point sources in a turbulent boundary layer. J. Fluid Mech. 117, 1–26. Franke, J., Hellsten, A., Schlünzen, H., Carissimo, B. (Eds.), 2007. Best Practice Guideline for the CFD Simulation of Flows in the Urban Environment. COST Action 732 Quality Assurance and Improvement of Microscale Meteorological Models. University of Hamburg. ISBN: 3-00-018312-4. Gailis, R.M., Hill, A., Yee, E., et al., 2007. Extension of a fuctuating plume model of tracer dispersion to a sheared boundary layer and to a large array of obstacles. Bound LayerMeteorol 122, 577–607. Goulart, E.V., Coceal, O., Belcher, S.E., 2018. Dispersion of a passive ScalarWithin and above an urban street network. Boundary-Layer Meteorol. 166, 351–366. Hertwig, D., Efthimiou, G.C., Bartzis, J.G., Leitl, B., 2012. CFD-RANS model validation of turbulent flow in a semi-idealized urban canopy. J. Wind Eng. Ind. Aerod. 111, 61–72. Hertwig, D., Soulhac, L., Fuka, V., Auerswald, T., Carpentieri, M., Hayden, P., Robins, A., Xie, Z.T., Coceal, O., 2018. Evaluation of fast atmospheric dispersion models in a regular street network. Environ. Fluid Mech. 18, 1007–1044. Jiang, G., Yoshie, R., 2018. Large-eddy simulation of flow and pollutant dispersion in a 3D urban street model located in an unstable boundary layer. Build. Environ. 142, 47–57. Koutsourakis, N., Bartzis, J.G., Markatos, N.C., 2012. Evaluation of Reynolds stress, k-ε and RNG k-ε turbulence models in street canyon flows using various experimental datasets. Environ. Fluid Mech. 12, 379–403. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 3, 269–289. Lewis, D.M., Chatwin, P.C., 1997. A three-parameter pdf for the concentration of an atmospheric pollutant. J. Appl. Meteorol. 36, 1064–1075. Li, H., Cui, G., Zhang, Z., 2018. A new scheme for the simulation of microscale flow and dispersion in urban areas by coupling large-eddy simulation with mesoscale models. Boundary-Layer Meteorol. 167, 145–170. Liu, Y.S., Miao, S.G., Zhang, C.L., Cui, G.X., Zhang, Z.S., 2012. Study on microatmospheric environment by coupling large eddy simulation with mesoscale model. J. Wind Eng. Ind. Aerod. 107–108, 106–117. Munro, R.J., Chatwin, P.C., Mole, N., 2003. A concentration pdf for the relative dispersion of a contaminant plume in the atmosphere. Bound Layer Meteor 106, 411–436, 2003. Nakayama, H., Leitl, B., Harms, F., Nagai, H., 2014. Development of local-scale highresolution atmospheric dispersion model using large-eddy simulation. Part 4: turbulent flows and plume dispersion in an actual urban area. J. Nucl. Sci. Technol. 51 (5), 626–638. Nironi, C., Salizzoni, P., Marro, M., et al., 2015. Dispersion of a passive scalar fuctuating plume in a turbulent boundary layer. Part I: velocity and concentration measurements. Bound Layer Meteorol 156, 415–446. Oberkampf, William L., Barone, Matthew F., 2006. Measures of agreement between computation and experiment: validation metrics. J. Comput. Phys. 217 (1), 5–36. Saad, Y., 2003. Iterative Methods for Sparse Linear Systems, second ed. SIAM, Philadelphia. Schatzmann, M., Olesen, H., Franke, J. (Eds.), 2010. Model Evaluation Case Studies: Approaches and Results. COST Action 732 COST Action 732 Quality Assurance and Improvement of Microscale Meteorological Models. University of Hamburg. ISBN 300-018312-4. Sykes, R.I., Lewellen, W.S., Parker, S.F., 1984. A turbulent-transport model for concentration fluctuations and fluxes. J. Fluid Mech. 139, 193–218. Tolias, I.C., Koutsourakis, N., Hertwig, D., Efthimiou, G.G., Venetsanos, A.G., Bartzis, J.G., 2018. Large Eddy Simulation study on structure of turbulent flow in a complex city. J. Wind Eng. Ind. Aerod. 177, 101–116. Tominaga, Y., et al., 2013. Cross comparisons of CFD results of wind and dispersion fields for MUST experiment: evaluation exercises by AIJ. J. Asian Architect. Build Eng. 12 (1), 117–124. Venetsanos, A., Papanikolaou, E., Bartzis, J., 2010. The ADREA-HF CFD code for consequence assessment of hydrogen applications. Int. J. Hydrogen Energy 35, 3908–3918. VDI Guideline 3783/9, 2005. Environmental Meteorology—prognostic Microscale Windfield Models—evaluation for Flow Around Buildings and Obstacles. Beuth Verlag, Berlin. Yee, E., Biltoft, C., 2004. Concentration fluctuation measurements in a plume dispersing through a regular array of obstacles. Bound.-Layer Meteorol. 111, 363–415. Yee, E., Wang, B.C., Lien, F.S., 2009. Probabilistic model for concentration fluctuations in compact source plumes in an urban environment. Boundary-Layer Meteorol. 130, 169–208.
Acknowledgments The simulation was supported by computational time granted by the Greek Research & Technology Network (GRNET) in the National HPC facility ARIS (http://hpc.grnet.gr) under project CFD-URB (pr004009). Present author would like to thank Dr Tolias Ilias from NCSR Demokritos for giving permission to use the MPI version of ADREA-HF. Present author would like to thank all members of COST Action 732 (https://mi-pub.cen.uni-hamburg.de/index.php?id¼484). References Andronopoulos, S., Bartzis, J.G., Würtz, J., Asimakopoulos, D., 1994. Modelling the effects of obstacles on the dispersion of denser-than-air gases. J. Hazard Mater. 37, 327–352. Andronopoulos, S., Grigoriadis, G., Robins, A., Venetsanos, A., Rafailidis, S., Bartzis, J.G., 2002. Three-dimensional modelling of concentration fluctuations in complicated geometry. Environ. Fluid Mech. 1, 415–440. Aristodemou, E., Boganegra, L.M., Mottet, L., Pavlidis, D., Constantinou, A., Pain, C., Robins, A., ApSimon, H., 2018. Environ. Pollut. 233, 782–796. Bartzis, J.G., 1991. ADREA-HF: a Three Dimensional Finite Volume Code for Vapour Cloud Dispersion in Complex Terrain. EC Joint Research Centre. CEC JRC Ispra Report EUR 13580 EN, Ispra. Bartzis, J.G., 2005. New approaches in two-equation turbulence modelling for atmospheric applications. Boundary-Layer Meteorol. 116 (3), 445–459. Bartzis, J.G., Efthimiou, G.C., Andronopoulos, S., 2015. Modelling short-term individual exposure from airborne hazardous releases in urban environments. J. Hazard Mater. 300, 182–188. Bezpalcova, K., Harms, F., 2005. EWTL Data Report/Part I: Summarized Test Description Mock Urban Setting Test, Technical Report, Environmental Wind Tunnel Lab. Center for Marine and Atmospheric Research, University of Hamburg. Chang, J.C., Hanna, S.R., 2005. Technical Descriptions and User's Guide for the BOOT Statistical Model Evaluation Software Package. www.harmo.org/kit. Cheng, W.C., Porte-Agel, F., 2016. Large-eddy simulation of flow and scalar dispersion in rural-to-urban transition regions. Int. J. Heat Fluid Flow 60, 47–60. Coceal, O., Goulart, E.V., Branford, S., Thomas, T.G., Belcher, S.E., 2014. Flow structure and near-field dispersion in arrays of building-like obstacles. J. Wind Eng. Ind. Aerod. 125, 52–68. Csanady, G.T., 1967. Concentration fluctuations in turbulent diffusion. J. Atmos. Sci. 24, 21–28. Donaldson, C. duP., 1973. Construction of a dynamic model for the production of atmospheric turbulence and the dispersal of atmospheric pollutants. In: Haugen, D.A. (Ed.), Workshop on Micrometeorology, 313–392. American Meteorological Society, Boston. Efthimiou, G.C., Andronopoulos, S., Tolias, I., Venetsanos, A., 2016. Prediction of the upper tail of concentration probability distributions of a continuous point source release in urban environments. Environ. Fluid Mech. 16 (5), 899–921. Efthimiou, G.C., Andronopoulos, S., Bartzis, J.G., 2017b. Evaluation of probability distributions for concentration fluctuations in a building array. Physica A 484, 104–116.
255
Contents lists available at ScienceDirect
Journal of Wind Engineering & Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia
Prediction of four concentration moments of an airborne material released from a point source in an urban environment George C. Efthimiou Environmental Research Laboratory, INRASTES, NCSR Demokritos, Patriarchou Grigoriou & Neapoleos Str., 15310, Aghia Paraskevi, Greece
A R T I C L E I N F O
A B S T R A C T
Keywords: RANS Beta distribution MUST Wind tunnel Skewness Kurtosis
Іn case of the dispersion of an airborne material from a point source in an urban environment the reliable prediction of the concentration statistical distribution by a numerical dispersion model presupposes the capability of the model to predict at least four statistical moments (mean, variance, skewness and kurtosis). In the present study, the beta distribution, the selection of which is justified based on a previous study, is incorporated in the Reynolds Averaged Navier Stokes (RANS) methodology. The shape parameters of the beta distribution are calculated using the numerical results of the mean, variance and maximum concentration. The latest is calculated through a deterministic model which uses also the numerical results of the mean and variance concentration as well as a hydrodynamic time scale. The validation of the new hybrid model “RANS-beta” is performed using the experimental dataset of the MUST wind tunnel experiment. The performance of the “RANS-beta” model for the skewness is very good (FAC2 ¼ 0.811) while for the kurtosis it is acceptable (FAC2 ¼ 0.557). The discrepancies are observed mainly at the edges of the plume. Future research will be focused on the optimization of the mean flow field, turbulent quantities (Reynolds stresses, turbulent kinetic energy) and on new parameterizations for the turbulent diffusion term.
1. Introduction Increasing urbanization and the resulting localization of human living and industrial production within relatively small areas lead to very demanding standards for human safety and security in urban environments. In this context, there is a necessity to be able to provide reliable decision support to minimize or even eliminate the potential health (or even life) hazard from accidental or deliberate atmospheric releases of harmful toxic, flammable, radioactive or bacteriological substances that might happen within complex urban and industrial areas. Accidental or deliberate airborne hazardous releases might occur in urban or industrial environments (e.g. COST Action ES1006, http://www.elizas.eu/) and they are threatening for human health and life, as people can be exposed to released substances through, e.g., inhalation, skin or eye contact. Additionally, everyday pollution sources like cars and chimneys, cause a permanent health risk for city dwellers. The chaotic nature of atmospheric turbulence increases the uncertainties for a realistic prediction of human exposure. For this reason, the use of statistics is imperative for this research field. Thus, the exposure is inherently stochastic due to turbulence variability. The problem becomes even more complex when dispersion occurs over complex
terrain such as the urban environment. If the concentration is measured for relatively short time periods at a specific location downwind of the source, then it can be observed that the range of the concentration values is wide. If a human is at the same location for this time period, then it is crucial to predict the probability – or the percentage of time - that the concentration will be higher than or equal to a threshold that could cause adverse health effects. In other words, the ability to predict the upper tail of the concentration probability distribution downstream of the source would be an important asset contributing to an efficient emergency management. Many researchers in the literature have proposed various statistical distributions to describe concentration fluctuations in the atmospheric surface layer such as exponential (Lewis and Chatwin, 1997), gamma (Nironi et al., 2015), beta (Bartzis et al., 2015), lognormal (Gailis et al., 2007), generalized Pareto (Munro et al., 2003) etc. In the last years, special interest is given by researchers about the gamma and the beta distribution. In Efthimiou et al. (2017b) we proved that the beta distribution predicts better the concentration fluctuations of an airborne material released from a point source in an urban-like area (MUST field experiment). The question that arises is that if this distribution can be incorporated in numerical dispersion models. The advantage of such an
E-mail address: [email protected]. https://doi.org/10.1016/j.jweia.2018.11.032 Received 29 May 2018; Received in revised form 3 October 2018; Accepted 24 November 2018 Available online 3 December 2018 0167-6105/© 2018 Elsevier Ltd. All rights reserved.
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
turbulence closure, Launder and Spalding, 1974) are calculated in the computational grid. Also, standard wall functions have been used in the present study for rough walls.
effort is that simple models that can predict two concentration moments (mean and variance) and a turbulent timescale will have the capability to predict the peak concentration. Then, the beta distribution can be used to predict all the concentration percentiles, necessary for probabilistic health impact assessments. The two higher moments (skewness and kurtosis) are used to validate the assumption of a beta distribution. Furthermore, in the literature, many researchers have selected the Computational Fluid Dynamics for the study of the flow field and the pollutant dispersion in urban areas. The modeling of the atmospheric turbulence is the main difficulty and the basic element of differentiation of research works. The most advanced method for this type of studies is the Direct Numerical Simulation (e.g. Goulart et al., 2018; Coceal et al., 2014). This method can be used e.g. for the understanding of the physical mechanisms of turbulence, the study of coherent structures etc. There is also the Large Eddy Simulation (e.g. Tolias et al., 2018; Jiang and Yoshie, 2018; Hertwig et al., 2018; Li et al., 2018; Aristodemou et al., 2018; Cheng and Port e-Agel, 2016; Nakayama et al., 2014; Liu et al., 2012) which has many applications e.g. it can be used also for basic understanding of turbulence and it can be used for emergency response. Concerning the latest, reliable conclusions are needed (e.g. Koutsourakis et al., 2012) in short times (e.g. less than 10s) and the selection of LES and the simple Reynolds Averaged Navier Stokes (RANS) model can be considered attractive. In this case, “smart” algorithms could extract the information from an LES or RANS database. The present author works this period on emergency preparedness and response tools in case of source identification (e.g. Efthimiou et al., 2018). This is the main reason why RANS is examined in this study. However, RANS can predict the first (mean) and the second (variance) moments through numerical solution of the corresponding transport equations. The incorporation of the beta distribution in the RANS methodology will allow the calculation of two additional moments (skewness and kurtosis). The objective of the present paper is to perform the incorporation of the beta distribution in the RANS methodology and to examine the validity of the hybrid “RANS-beta” model. A similar work has been performed in Efthimiou et al. (2016) where the hybrid “RANS-gamma" model was validated with wind tunnel measurements of the Michelstadt and CUTE experiments. In the present paper, the validation is performed with the MUST wind tunnel experiment. In this work the in-depth investigation of sophisticated numerical schemes, turbulence models etc. is not performed but the interest is focused on the operationality of the method in order to be used for emergency response in the future. Optimization of the methodology will be performed in a future study. In Chapter 2 the methodology is described i.e. the incorporation of the beta distribution in the RANS methodology. In Chapter 3 the experiment is described, in Chapter 4 the numerical simulations and in Chapter 5 the results are presented. The conclusions are presented in Chapter 6.
2.2. Numerical solution of the concentration mean and variance equations In this step, the partial differential equations (pde) of concentration mean, C; and variance, C'2 , are solved numerically in the computational grid using the velocities and the turbulent parameters of step 1. The pde of the mean concentration in Cartesian coordinates is usually expressed as follows:
∂ ∂ ∂ ∂C ρC þ ρui C ¼ ρD ρu0 i C0 ∂xi ∂t ∂xi ∂xi
þ SC
(1)
In Eq. (1), C is the mean concentration, t is the time (s), ρ is the density of the air-pollutant mixture (kg m3), xi is the coordinate (m) along the i0 direction,ui is the mean velocity (m s1) in the i-direction, u i is the ve1 locity fluctuation (m s ) in the i-direction, C’ is the concentration fluctuation, D is the molecular mass diffusivity (m2 s1) and SC is the pollutant sources/sinks. The overbar denotes Reynolds average, while repeated indices in a term imply summation over all possible values of the index, according to the Einstein summation convention. To close Eq. (1) the turbulent diffusion term u0 i C0 needs to be modeled. This is done in the present study through the gradient transfer assumption, due to its wide use in the modeling of engineering flows and its computational simplicity. So, following Csanady (1967) the turbulent diffusion term is modeled using the eddy mass diffusivityK C : u0 i C 0 ¼ K C
∂C K m ∂C ¼ ∂xi σ C ∂xi
(2)
In Eq. (2) the eddy mass diffusivity is calculated by dividing the eddy viscosity K m by the turbulent Schmidt number σ C. The pde for the mean square fluctuation of concentration is expressed as follows:
0
1
∂ ∂ ∂C ∂ @ ∂C 0 2 2 ρC0 2 þ ρui C0 2 ¼ 2ρu0 i C0 þ ρD ρu0 i C 0 A ∂xi ∂xi ∂xi ∂t ∂xi 2ρD
∂C 0 ∂C 0 ∂xi ∂xi
(3)
where C0 2 is the concentration variance. Eq. (3) accounts for all the physical phenomena related to the advection, production, diffusion and dissipation of this quantity. To close Eq. (3), two terms need to be modeled: the turbulent diffusion term and the dissipation term.
2. Methodology
The turbulent diffusion term u0 C0 2 is modeled by using again the
The present methodology consists of four steps:
0
gradient transport hypothesis with an eddy diffusion coefficient K C that is assumed to be equal to the eddy mass diffusivity K C (Csanady, 1967; Donaldson, 1973; Sykes et al., 1984):
1. Numerical simulation of the hydrodynamic problem. 2. Numerical solution of the concentration mean and variance equations using the results of step 1. 3. Calculation of maximum concentration using the results of steps 1 and 2. 4. Calculation of concentration skewness and kurtosis using the results of steps 2 and 3.
2
u0 i C 0 ¼ K C
0
∂C 0 2 ∂C 0 2 K m ∂C 0 2 ¼ K C ¼ ∂xi ∂xi σ C ∂xi
(4)
The dissipation term (last term) of Eq. (3) is modeled as the ratio of the mean square fluctuations of the mass fraction to the dissipation time
Analysis of each step follows.
0
2
scale T C of C0 . Following Fackrell and Robins (1982), the dissipation time scale is expressed as the ratio of a length scale LC to the velocity scale k1/2, where k is the turbulent kinetic energy:
2.1. Numerical simulation of the hydrodynamic problem In this step, the 3-dimensional, steady-state fields of the wind velocity components u, v, w, the turbulence kinetic energy, k, and the turbulence kinetic energy dissipation rate, ε, (in case of using the standard k-ε
D
248
1=2 0 2 ∂C 0 ∂C 0 C 0 2 C 0 k ¼ 0 ¼CC ∂xi ∂xi T C LC
(5)
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
0
a¼
where C C is a constant equal to 0.5 and LC is an effective length scale of turbulence which for the selected turbulence model (standard kε) is 1:5 given by LC ¼ C0:75 μ k =ðκ εÞ, where κ ¼ 0.4 is the von Karman constant and Cμ ¼ 0.09. In the present study the eddy viscosity K m (appearing in Eqs. (2) and (4)) is calculated by the standard k-ε model, which solves transport equations for the turbulent kinetic energy k and the dissipation rate ε.
1 η 1 1þη I
ξ ¼ ηa
η¼
Cmax ðΔτÞ C C
(13) (14)
(15)
It is reminded that the variable Cmax(Δτ) is calculated in step 3 while C in step 2.
2.3. Calculation of maximum concentration 3. The MUST wind tunnel experiment
In this step, the maximum concentration at an interval time Δτ is calculated by the deterministic model Bartzis et al. (2015). The equation is: " Cmax ðΔτÞ ¼ 1 þ 9I
Δτ TC
The methodology has been validated against data of the MUST wind tunnel experiment (Bezpalcova and Harms, 2005) which have been scaled up for the conditions of the corresponding field experiment (Yee and Biltoft, 2004). In MUST experiment the obstacles were arranged in 12 rows, each consisting of 10 obstacles. The obstacles were nearly identical and had average length, width and height 12.2 m 2.42 m 2.54 m respectively. The contaminant's concentration has been measured by a 256-detectors array arranged along obstacle rows in the part of the domain covered by the plume (Fig. 1). All detectors were placed at the same height equal to 1.275 m. Non-zero concentrations have been measured by nearly all (244 out of 256) detectors. Here it is noted that only the measurement data that were available to the authors (through data base of COST Action 732) are taken into account. The link to the database of COST Action 732 is: https://mi-pub.cen.uni-hamburg.de/index.php?id¼484. The wind flow was characterized by neutral stratification, wind speed at the roof level Uref ¼ 8 m/s, (field scale) and wind direction 45 (Fig. 1) in the experimental coordinate system. The contaminant originated from a point source located at the ground level. The volume flow rate of gas at the source was: 3.3 106 m3s1. Wind tunnel measurements of u and v (and w for intercomparison purposes with Tominaga et al., 2013) were used also to validate the numerical results of the hydrodynamic problem. Two groups of sensors have been used which correspond to measurements in a coarse network and in a fine network (Fig. 2). The height of the velocity and concentration sensors is the same and equal to 1.275 m.
0:3 # C
(6)
where TC is the turbulence integral time scale, which under fully mixed conditions it is given by the relationship: TC ¼
Cμ k 3=2 ; Cμ ¼ 0:09 Uε
(7)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and Uis the mean wind speed: U ¼ u2 þ v2 þ w2 (ms1). The fluctuation intensity I is defined as: I¼
C0
2
C
2
(8)
It is reminded that the variables u, v, w, k and ε are calculated in step 1 while C and C'2 in step 2. Δτ is the exposure time considered. The relevant duration of the exposure time Δτ is selected either on health impact considerations or on the hazardous substance release duration or on real conditions (e.g., the time an individual is expected to stay at the particular location). Especially near the emission source there might be cases where the concentrations are very high and related time scales are as low as breathing time (i.e., a few seconds). On the other hand very short releases in the order of a few seconds up to a few minutes are common scenarios for the above mentioned accidental or deliberate release events. In this work the averaging time interval Δτ is considered equal to the measurement time interval of the instrument.
2.4. Calculation of concentration skewness and kurtosis The general definitions of the skewness (Sk) and kurtosis (Ku) are: Sk ¼
Ku ¼
C0 3 C0 2 C0 C0
1:5
(9)
4
2
2
(10)
In case of the beta distribution Sk and Ku are calculated based on the following theoretical equations (Efthimiou et al., 2017b): Sk ¼
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ðξ aÞ a þ ξ þ 1 pffiffiffiffiffi ða þ ξ þ 2Þ aξ
(11)
Ku ¼
6 ða ξÞ2 ða þ ξ þ 1Þ aξða þ ξ þ 2Þ þ3 aξða þ ξ þ 2Þða þ ξ þ 3Þ
(12)
Fig. 1. The computational domain of the MUST case. The sensors are presented with yellow circles. The source is depicted with a star. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
where 249
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
Fig. 2. The velocity sensors (yellow circles) of the MUST wind tunnel experiment placed at 1.275 m above the ground (left: coarse network, right: fine network). The white rectangles are the obstacles. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)
“More details about the experiment can be found in the website of COST Action 732 (https://mi-pub.cen.uni-hamburg.de/index.php? id¼484) and personal communication with the experimentalists.”
should be noticed that in the experiment not the entire bottom was rough but only a part of it. Concerning CFD guidance on the treatment of walls, detailed information can be found in Franke et al. (2007). The fixed vertical profiles at the inlet planes (-x and þy) for u, v, k and ε have been constructed based on the experimental corresponding profiles, while for w a zero value has been set. The outlet boundary conditions at planes þ x and –y are justified, as these planes are located at a sufficient distance (see above) downwind of the last buildings. For turbulence modeling, the standard k–ε model (Launder and Spalding, 1974) was used. The problem was treated as a pseudo-transient case i.e. time integrated towards steady state. The unsteady Reynolds-Averaged NavierStokes equations for total mass, the three components of momentum, the turbulence kinetic energy and its dissipation rate were solved. At the end of the solution, steady-state values were produced i.e. constant (in time) values of the variables at the sensors. The total calculation time (1200s) was set about equal to the travel time from the source point to the outlet of the computational domain of a hypothetical particle moving with the lowest speed observed in the experiment. This time was sufficient to reach a steady-state situation. The time step was automatically adapted according to the desired Courant–Friedrichs–Lewy (CFL) number whose maximum value was set to 2 based on the author's experience for numerical stability (better convergence behaviour of the solver). Nevertheless, as an implicit scheme is adopted and a steady-state solution is sought for, the magnitude of the final time step does not play here an important role. A first order scheme was used for the time derivative. For the calculation of pressure, the linear system that arises from the discretized equations was solved using the preconditioned BiCGstab method. A block-ILU[0] type of preconditioner was used (additive Swartz preconditioner, Saad, 2003).
4. The numerical simulations All computations (wind flow and dispersion) were performed with the CFD code ADREA-HF (Bartzis, 1991; Andronopoulos et al. 1994, 2002; Venetsanos et al., 2010). This is a Eulerian model, solving the unsteady RANS equations. The code uses finite volumes with rectangular cells for the discretization of the Eulerian transport equations. To describe the complex geometry, the volume porosity concept is used with solid surfaces of any orientation allowed to cross the computational cells. The ADREA-HF code has been setup for the simulation in the field scale. 4.1. The wind field simulations The spatial discretization of the computational domain is presented in Table 1. The domain extends horizontally by 17.55 m (6.9 H, H ¼ 2.54 m being the buildings height) upwind of the first buildings and by 52.65 m (20.7 H) downwind of the last buildings. The vertical dimension of the domain is about 8.3 H. The above dimensions conform to the recommendations of COST Action 732 (Franke et al., 2007). The grid is Cartesian. Cubic cells have been selected inside the urban area in order to decrease the discretization errors. Outside the urban area, the grid increases logarithmically by a factor of 1.1. It is well known by CFD practices for microscale modeling (Franke et al., 2007) that the building height should be described by at least ten cells. Thus, in case of MUST, a minimum cell size should be equal to dxmin ¼ dymin ¼ dzmin ¼ 0.25 m. The total number of cells of the present simulation is 21,807,800 which is finer than the one used in Efthimiou et al. (2017c). The boundary conditions for the hydrodynamic variables at each boundary plane or surface of the domain are presented in Table 2. It
Table 2 Boundary conditions for the hydrodynamic variables (u, v, w: velocity components in the x-, y-, z-axis respectively, k: turbulence kinetic energy, ε: turbulence kinetic energy dissipation rate).
Table 1 Dimensions and spatial discretization of the computational domain. Domain dimensions x/ y/z (m)
Total number of cells
Number of cells in each axis
Minimum/maximum cell sizes (m)
x
y
z
dx
dy
dz
243/268/21
21,807,800
740
842
35
0.25/ 5.1
0.25/ 5.1
0.25/ 1.84
250
Plane
Boundary condition
-x þx -y þy Ground Top Building walls
Inlet: fixed vertical profiles for u, v, k, ε, w ¼ 0 Outlet: ∂ϕ=∂x ¼ 0, φ ¼ u, v, w, k, ε Outlet: ∂ϕ=∂y ¼ 0, φ ¼ u, v, w, k, ε Inlet: fixed vertical profiles for u, v, k, ε, w ¼ 0 Standard wall functions, roughness length ¼ 0.0165 m Fixed values for u, v, k, ε, w is calculated from the cell mass balance Standard wall functions, roughness length ¼ 1 105 m
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
experiment. For example, in Tominaga et al. (2013) the HRs of 8 models, for 498 points (xz plane), for the velocity component u was ranged from 0.63 to 0.91 while for the velocity component w was ranged from 0.12 to 0.26. The HR of the same velocity components of the present simulation is inside these ranges (0.75 for u/Uref and 0.18 for w/Uref) indicating the consistency with the results obtained by Tominaga et al. (2013) and by COST 732 group (Schatzmann et al., 2010). The HR of each velocity component at each sensor of Fig. 2 is presented in Fig. 3. The black circle indicates HR ¼ 1 and a good prediction while the red circle indicates HR ¼ 0 and a poor prediction. For the velocity u/Uref, there are a lot of hits everywhere in the urban area. On the other hand, the HR of v/Uref is good only upwind and in the first rows of buildings. It should be noticed that at the selected elevation the influence of the prescribed wall boundary conditions is strong and most likely causing the observed deviations from the experiment. Furthermore, the numerical data were linearly interpolated in order to match the spatial positions of the experimental measurements. This procedure adds further uncertainties to the predictions which are hard to distinguish from errors that originated from the simulation itself. However, the error/uncertainty introduced by interpolation should be viewed part of the
The 1st order upwind scheme has been used for approximation of convective terms. The second order central-difference scheme has been used for approximation of diffusion terms. It should be noticed that, since the grid is finer than the one used in Efthimiou et al. (2017c), where a grid independent numerical solution was attained, the selection of the numerical schemes for the spatial derivatives does not have significant effects on the results. The computations were run in parallel using 20 cores of 5 processors (i.e. 100 processors). ADREA-HF is parallelized for distributed memory architectures using the MPI protocol. 4.2. The pollutant dispersion simulations The transport equations for the concentration mean and variance were solved on the grid simultaneously with the hydrodynamic problem from which the flow and turbulence related variables are used for the advection and diffusion terms. The boundary conditions of the concentration mean (C) and variance (C'2 ) for each boundary plane or surface of the domain are presented in Table 3. The total calculation time was set equal to the one used for the hydrodynamic computations, i.e., the 1200s, which was adequate for obtaining a steady-state solution. For the discretization of the convective term in the transport equations, the upwind scheme was used, while for the diffusion term the central scheme was used. The time step was automatically adapted according to the desired Courant–Friedrichs–Lewy (CFL) number whose maximum value was set to 2. 5. Results and discussion 5.1. Validation of the hydrodynamic problem solution Author has checked in a past paper the profile development across the domain of the same experiment (Efthimiou et al., 2017c). It was noticed that the model presented a very good performance above the buildings, where the Hit Rate (HR) (VDI Guideline 3783/9, 2005) had almost the ideal value (¼1.0) for the velocities u/Uref and v/Uref. Below the building roofs the HR was gradually decreased for all velocity components. u/Uref presented in total better agreement with the experimental data than v/Uref and this is normally the case for the components being aligned with the mean flow direction in x-direction. The reason for this is that the magnitudes of the lateral component, i.e. in y-direction, are very much smaller than the magnitudes of the component in the mean flow direction. In this case the coordinate system corresponds to the wind tunnel. Very close to the ground the velocity v/Uref presented the lowest HRs. The HR of the velocities u/Uref and v/Uref for the measurements below the half height of the buildings was found equal to 0.48 and 0.14 respectively in Efthimiou et al. (2017c). Similar values of the HR were found also in the present simulation which for the specific sensors at the height of 1.275 m (Fig. 2) are equal to 0.55 and 0.16 respectively. It should be noticed that in the literature there are exercises where various models have been intercompared using the MUST wind tunnel Table 3 Boundary conditions for the concentration mean (C) and variance (C'2 ). Plane
Boundary condition
-x
∂C=∂x ¼ 0, ∂C0 2 =∂x ¼ 0
þx
∂C=∂x ¼ 0, ∂C0 2 =∂x ¼ 0
-y
∂C=∂y ¼ 0, ∂C0 2 =∂y ¼ 0
þy
∂C=∂y ¼ 0, ∂C0 2 =∂y ¼ 0
Ground
∂C=∂z ¼ 0, ∂C0 2 =∂z ¼ 0
Top
∂C=∂z ¼ 0, ∂C0 2 =∂z ¼ 0
Building walls
∂C=∂n ¼ 0, ∂C0 2 =∂n ¼ 0, n ¼ x, y, z
Source (mass fraction)
C ¼ 1, ∂C=∂n ¼ 0, ∂C0 2 =∂n ¼ 0, n ¼ x, y, z
Fig. 3. The Hit Rate (HR) per sensor of the velocity components u/Uref (upper) and v/Uref (lower). The black circle indicates HR ¼ 1 and a good prediction while the red circle indicates HR ¼ 0 and a poor prediction. The results are presented for height 1.275 m above the ground for the coarse and fine networks. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) 251
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
5.2. Validation of the concentration statistical moments
simulation (Oberkampf and Barone., 2006). Finally, as has been discussed in Hertwig et al. (2012), unsteady flow effects triggered by the presence of buildings are a further influencing factor for the model performance. Another way to compare the results of the model and of the experiment is the vector plots. Fig. 4 presents horizontal vector plots of the experiment (upper) and the model (lower) at the area of the fine network of sensors (Fig. 1). In this case the coordinate system corresponds to the model and it is not the same as the one used for the calculation of the HRs. The direction of the flow is almost the same but an underprediction of the velocity magnitudes by the model is clear. This underprediction was observed also by author in a past work (Hertwig at al. 2012).
In the present study, the validation of the four concentration statistical moments (mean, variance, skewness and kurtosis) was performed using scatter plots, validation metrics and contour point-by-point plots. The experimental database provided normalized results of the concentration mean and root mean square (RMS) and as a result, the numerical results were normalized in the same way based on the following equation: C* ¼
CH 2 Uref Q
Fig. 4. Horizontal vector plots of the experiment (upper) and the model (lower) for the area covered by the fine network of sensors (Fig. 1). 252
(16)
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
intermittency). The model predicted very well the values near the centerline as many of the points lie on the 1 to 1 line. The poor prediction of the highest skewness values (at the edges of the plume) is an indication that the concentration mean and RMS are inversely proportional to the skewness. The optimization of the turbulent diffusion term will optimize also the performance of the highest skewness values. The performance is also good for the kurtosis. The positive sign reveals that the distributions are leptokurtic. The sensors at the edges of the plume have much fatter tails than the sensors at the plume centerline. Also, the scatter of the points around the ideal line is higher than skewness. Table 4 presents the Normalized Mean Square Error (NMSE), the Fractional Bias (FB) and the Factor of 2 (FAC2) for the four concentration moments. The explanation of these metrics is given in Chang and Hanna (2005). The ideal values are NMSE ¼ 0, FB ¼ 0 and FAC2 ¼ 1. Based on the available information, the experimental confidence is provided only for the concentration mean and RMS and is equal to 0.003. However, this value was not considered in the present study for the calculation of the FAC2. The FAC2 of the skewness (0.811) is very good and better than the other moments supporting the robustness of the methodology. However, the validation metrics can lead to erroneous conclusions about the performance of a model. Even if in Fig. 5 the scatter of the points around the ideal line (1–1) for the kurtosis is not so high in
where C is the predicted concentration (either mean or RMS), H is height of the buildings, Uref is the reference velocity and Q is the source rate. Fig. 5 presents the scatter plots with logarithmic axes of the concentration mean, RMS, skewness and kurtosis. For the concentration mean we observe overprediction for the higher and the medium values and underprediction for the lower values. For the RMS the prediction is better than the mean as more values lie on the 1 to 1 line and we observe also the underestimation for the lower values. There are mainly two parameters that influence the performance of the model: a) the mean flow field and turbulent quantities (Reynolds stresses, turbulent kinetic energy) which are the driving force of the dispersion. The optimization in this case is performed with sophisticated turbulence models (e.g. Bartzis, 2005) and boundary conditions (e.g. Franke et al., 2007) and b) the turbulent diffusion term. It is reminded that this term is based on the Schmidt number the value of which is constant everywhere and this raises questions about the validity of this assumption. Furthermore, the turbulent diffusion term is based on the full mixing hypothesis which means that we assume that the concentration and the wind flow have the same length/time turbulent scales everywhere. Yee et al. (2009) have separated these scales, using also the pollutant travel time, for the calculation of the dissipation time scale of the concentration variance and their results were very good. Beyond concentration variance, Efthimiou and Bartzis (2011) have used the pollutant travel time for the calculation of the concentration autocorrelation time scale. The usefulness of this parameter is high as it can be used in combination with the eddy viscosity and the hydrodynamic autocorrelation time scale for the calculation of the turbulent diffusion term. This is a task that will be performed by author in the future. Concerning the concentration skewness, the positive sign reveals that the right tails are longer. Close to the plume centerline (lower intermittency) the skewness is lower than at the plume edges (higher
Table 4 Validation metrics (NMSE, FAC2 and FB) for the concentration moments.
Mean RMS Skewness Kurtosis
NMSE
FAC2
FB
2.85 4.42 47.00 927.21
0.496 0.627 0.811 0.557
0.267 0.036 0.765 1.701
Fig. 5. Scatter plots of the concentration mean, RMS, skewness and kurtosis. The mean and RMS are normalized values to be in accordance with the experiment. 253
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
comparison with the other moments, its NMSE is very high and this is due to the very high discrepancies of the model from the experiment at the edges of the plume. Similar behaviour is observed also for the skewness. The higher overprediction is observed for the kurtosis and then for the skewness and this is again due to the highest values. The overprediction of the concentration mean (FB ¼ 0.267) is based on many parameters such as the turbulent diffusion term, the underprediction of the velocity field by the model, the turbulence model, the numerical schemes etc. Overall the results are acceptable. The FAC2 of the mean concentration (0.5) is close to the lower values reported by Tominaga et al. (2013) (from 0.55 to 0.88), for 256 points (xy plane). However, in the present calculation of the FAC2 all the experimental concentration values were considered. Finally, in Fig. 6 and Fig. 7, point by point contour plots are presented for the skewness and kurtosis respectively for the experiment (upper) and the model (lower). The range of values for each color is common for both moments. The similarities between the experiment and the model are very good for both moments especially at the main body of the plume
supporting the robustness of the “RANS-beta” model. The discrepancies are observed at the edges of the plume supporting the previous thoughts.
Fig. 6. Contour point-by-point plots of the concentration skewness for the experiment (upper) and the model (lower).
Fig. 7. Contour point-by-point plots of the concentration kurtosis for the experiment (upper) and the model (lower).
6. Conclusions In this work the operationality of the new, hybrid “RANS-beta” model was examined in order to predict four concentration moments (mean, variance, skewness and kurtosis) in case of the dispersion of an airborne material from a point source in the atmospheric surface layer. The computational cost concerns only the numerical simulation of the partial differential equations of concentration mean and variance as the concentration skewness and kurtosis are calculated during the postprocessing of the results. This process, overall, can be considered as an additional advantage for the RANS methodology during an emergency response. Concluding, it seems that, both RANS and LES are viable methods to generate data bases for operational emergency response modeling.
254
G.C. Efthimiou
Journal of Wind Engineering & Industrial Aerodynamics 184 (2019) 247–255
The validation was focused mainly on the skewness and the kurtosis. The agreement of the model with the MUST wind tunnel experiment for the concentration skewness is high (FAC2 ¼ 0.811) and acceptable for the kurtosis (FAC2 ¼ 0.557). The errors in the prediction of the concentration mean and variance were transferred in the prediction of concentration skewness and kurtosis and this was obvious at the edges of the plumes. Initially, optimization of the mean flow field is necessary as it is the driving force of the dispersion. Also, more sophisticated parameterizations for the turbulent diffusion term are necessary. These could include the relation of this coefficient with the pollutant travel time and the autocorrelation time scales of the concentration and the wind flow. These issues will be examined more carefully by author in a future publication. With better results for the prediction of the flow field and lower concentration moments one might obtain better results for skewness and kurtosis. Optimization of the methodology using e.g. higher order numerical schemes, more sophisticated turbulence models will be performed in a future study. It would be interesting also to incorporate the beta distribution in Large Eddy Simulation and this will be a future task.
Efthimiou, G.C., Bartzis, J.G., 2011. Atmospheric dispersion and individual exposure of hazardous materials. J. Hazard Mater. 188, 375–383. Efthimiou, G.C., Kovalets, I.V., Argyropoulos, C.D., Venetsanos, A., Andronopoulos, S., Kakosimos, K., 2018. Evaluation of an inverse modelling methodology for the prediction of a stationary point pollutant source in complex urban environments. Build. Environ. 143, 107–119. Efthimiou, G.C., Kovalets, I.V., Venetsanos, A., Andronopoulos, S., Argyropoulos, C.D., Kakosimos, K., 2017c. An optimized inverse modelling method for determining the location and strength of a point source releasing airborne material in urban environment. Atmos. Environ. 170, 118–129. Fackrell, J.E., Robins, A.G., 1982. Concentration fluctuations and fluxes in plumes from point sources in a turbulent boundary layer. J. Fluid Mech. 117, 1–26. Franke, J., Hellsten, A., Schlünzen, H., Carissimo, B. (Eds.), 2007. Best Practice Guideline for the CFD Simulation of Flows in the Urban Environment. COST Action 732 Quality Assurance and Improvement of Microscale Meteorological Models. University of Hamburg. ISBN: 3-00-018312-4. Gailis, R.M., Hill, A., Yee, E., et al., 2007. Extension of a fuctuating plume model of tracer dispersion to a sheared boundary layer and to a large array of obstacles. Bound LayerMeteorol 122, 577–607. Goulart, E.V., Coceal, O., Belcher, S.E., 2018. Dispersion of a passive ScalarWithin and above an urban street network. Boundary-Layer Meteorol. 166, 351–366. Hertwig, D., Efthimiou, G.C., Bartzis, J.G., Leitl, B., 2012. CFD-RANS model validation of turbulent flow in a semi-idealized urban canopy. J. Wind Eng. Ind. Aerod. 111, 61–72. Hertwig, D., Soulhac, L., Fuka, V., Auerswald, T., Carpentieri, M., Hayden, P., Robins, A., Xie, Z.T., Coceal, O., 2018. Evaluation of fast atmospheric dispersion models in a regular street network. Environ. Fluid Mech. 18, 1007–1044. Jiang, G., Yoshie, R., 2018. Large-eddy simulation of flow and pollutant dispersion in a 3D urban street model located in an unstable boundary layer. Build. Environ. 142, 47–57. Koutsourakis, N., Bartzis, J.G., Markatos, N.C., 2012. Evaluation of Reynolds stress, k-ε and RNG k-ε turbulence models in street canyon flows using various experimental datasets. Environ. Fluid Mech. 12, 379–403. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 3, 269–289. Lewis, D.M., Chatwin, P.C., 1997. A three-parameter pdf for the concentration of an atmospheric pollutant. J. Appl. Meteorol. 36, 1064–1075. Li, H., Cui, G., Zhang, Z., 2018. A new scheme for the simulation of microscale flow and dispersion in urban areas by coupling large-eddy simulation with mesoscale models. Boundary-Layer Meteorol. 167, 145–170. Liu, Y.S., Miao, S.G., Zhang, C.L., Cui, G.X., Zhang, Z.S., 2012. Study on microatmospheric environment by coupling large eddy simulation with mesoscale model. J. Wind Eng. Ind. Aerod. 107–108, 106–117. Munro, R.J., Chatwin, P.C., Mole, N., 2003. A concentration pdf for the relative dispersion of a contaminant plume in the atmosphere. Bound Layer Meteor 106, 411–436, 2003. Nakayama, H., Leitl, B., Harms, F., Nagai, H., 2014. Development of local-scale highresolution atmospheric dispersion model using large-eddy simulation. Part 4: turbulent flows and plume dispersion in an actual urban area. J. Nucl. Sci. Technol. 51 (5), 626–638. Nironi, C., Salizzoni, P., Marro, M., et al., 2015. Dispersion of a passive scalar fuctuating plume in a turbulent boundary layer. Part I: velocity and concentration measurements. Bound Layer Meteorol 156, 415–446. Oberkampf, William L., Barone, Matthew F., 2006. Measures of agreement between computation and experiment: validation metrics. J. Comput. Phys. 217 (1), 5–36. Saad, Y., 2003. Iterative Methods for Sparse Linear Systems, second ed. SIAM, Philadelphia. Schatzmann, M., Olesen, H., Franke, J. (Eds.), 2010. Model Evaluation Case Studies: Approaches and Results. COST Action 732 COST Action 732 Quality Assurance and Improvement of Microscale Meteorological Models. University of Hamburg. ISBN 300-018312-4. Sykes, R.I., Lewellen, W.S., Parker, S.F., 1984. A turbulent-transport model for concentration fluctuations and fluxes. J. Fluid Mech. 139, 193–218. Tolias, I.C., Koutsourakis, N., Hertwig, D., Efthimiou, G.G., Venetsanos, A.G., Bartzis, J.G., 2018. Large Eddy Simulation study on structure of turbulent flow in a complex city. J. Wind Eng. Ind. Aerod. 177, 101–116. Tominaga, Y., et al., 2013. Cross comparisons of CFD results of wind and dispersion fields for MUST experiment: evaluation exercises by AIJ. J. Asian Architect. Build Eng. 12 (1), 117–124. Venetsanos, A., Papanikolaou, E., Bartzis, J., 2010. The ADREA-HF CFD code for consequence assessment of hydrogen applications. Int. J. Hydrogen Energy 35, 3908–3918. VDI Guideline 3783/9, 2005. Environmental Meteorology—prognostic Microscale Windfield Models—evaluation for Flow Around Buildings and Obstacles. Beuth Verlag, Berlin. Yee, E., Biltoft, C., 2004. Concentration fluctuation measurements in a plume dispersing through a regular array of obstacles. Bound.-Layer Meteorol. 111, 363–415. Yee, E., Wang, B.C., Lien, F.S., 2009. Probabilistic model for concentration fluctuations in compact source plumes in an urban environment. Boundary-Layer Meteorol. 130, 169–208.
Acknowledgments The simulation was supported by computational time granted by the Greek Research & Technology Network (GRNET) in the National HPC facility ARIS (http://hpc.grnet.gr) under project CFD-URB (pr004009). Present author would like to thank Dr Tolias Ilias from NCSR Demokritos for giving permission to use the MPI version of ADREA-HF. Present author would like to thank all members of COST Action 732 (https://mi-pub.cen.uni-hamburg.de/index.php?id¼484). References Andronopoulos, S., Bartzis, J.G., Würtz, J., Asimakopoulos, D., 1994. Modelling the effects of obstacles on the dispersion of denser-than-air gases. J. Hazard Mater. 37, 327–352. Andronopoulos, S., Grigoriadis, G., Robins, A., Venetsanos, A., Rafailidis, S., Bartzis, J.G., 2002. Three-dimensional modelling of concentration fluctuations in complicated geometry. Environ. Fluid Mech. 1, 415–440. Aristodemou, E., Boganegra, L.M., Mottet, L., Pavlidis, D., Constantinou, A., Pain, C., Robins, A., ApSimon, H., 2018. Environ. Pollut. 233, 782–796. Bartzis, J.G., 1991. ADREA-HF: a Three Dimensional Finite Volume Code for Vapour Cloud Dispersion in Complex Terrain. EC Joint Research Centre. CEC JRC Ispra Report EUR 13580 EN, Ispra. Bartzis, J.G., 2005. New approaches in two-equation turbulence modelling for atmospheric applications. Boundary-Layer Meteorol. 116 (3), 445–459. Bartzis, J.G., Efthimiou, G.C., Andronopoulos, S., 2015. Modelling short-term individual exposure from airborne hazardous releases in urban environments. J. Hazard Mater. 300, 182–188. Bezpalcova, K., Harms, F., 2005. EWTL Data Report/Part I: Summarized Test Description Mock Urban Setting Test, Technical Report, Environmental Wind Tunnel Lab. Center for Marine and Atmospheric Research, University of Hamburg. Chang, J.C., Hanna, S.R., 2005. Technical Descriptions and User's Guide for the BOOT Statistical Model Evaluation Software Package. www.harmo.org/kit. Cheng, W.C., Porte-Agel, F., 2016. Large-eddy simulation of flow and scalar dispersion in rural-to-urban transition regions. Int. J. Heat Fluid Flow 60, 47–60. Coceal, O., Goulart, E.V., Branford, S., Thomas, T.G., Belcher, S.E., 2014. Flow structure and near-field dispersion in arrays of building-like obstacles. J. Wind Eng. Ind. Aerod. 125, 52–68. Csanady, G.T., 1967. Concentration fluctuations in turbulent diffusion. J. Atmos. Sci. 24, 21–28. Donaldson, C. duP., 1973. Construction of a dynamic model for the production of atmospheric turbulence and the dispersal of atmospheric pollutants. In: Haugen, D.A. (Ed.), Workshop on Micrometeorology, 313–392. American Meteorological Society, Boston. Efthimiou, G.C., Andronopoulos, S., Tolias, I., Venetsanos, A., 2016. Prediction of the upper tail of concentration probability distributions of a continuous point source release in urban environments. Environ. Fluid Mech. 16 (5), 899–921. Efthimiou, G.C., Andronopoulos, S., Bartzis, J.G., 2017b. Evaluation of probability distributions for concentration fluctuations in a building array. Physica A 484, 104–116.
255