Joint PDF modelling of turbulent flow and dispersion in an urban street canyon

Joint PDF modelling of turbulent flow and dispersion in an urban street canyon

J. Bakosi Los Alamos National Laboratory, Los Alamos, NM 87545 P. Franzese and Z. Boybeyi College of Science, George Mason University, Fairfax, VA 22030

The joint probability density function (PDF) of turbulent velocity and concentration of a passive scalar in an urban street canyon is computed using a newly developed particle-in-cell Monte Carlo method. Compared to moment closures, the PDF methodology provides the full one-point one-time PDF of the underlying fields containing all higher moments and correlations. The small-scale mixing of the scalar released from a concentrated source at the street level is modelled by the interaction by exchange with the conditional mean (IECM) model, with a micro-mixing time scale designed for geometrically complex settings. The boundary layer along no-slip walls (building sides and tops) is fully resolved using an elliptic relaxation technique, which captures the high anisotropy and inhomogeneity of the Reynolds stress tensor in these regions. A less computationally intensive technique based on wall functions to represent boundary layers and its effect on the solution are also explored. The calculated statistics are compared to experimental data and large-eddy simulation. The present work can be considered as the first example of computation of the full joint PDF of velocity and a transported passive scalar in an urban setting. The methodology proves successful in providing high level statistical information on the turbulence and pollutant concentration fields in complex urban scenarios.

Langevin equation; Monte-Carlo method; Probability density function method; Scalar dispersion; Urban-scale turbulence
journal: Boundary-Layer MeteorologyCorresponding author.

1 Introduction

Regulatory bodies, architects and town planners increasingly use computer models to assess ventilation and occurrences of hazardous pollutant concentrations in cities. These models are mostly based on the Reynolds-averaged Navier-Stokes (RANS) equations or, more recently, large-eddy simulation (LES) techniques. Both of these approaches require a series of modelling assumptions, including most commonly the eddy-viscosity and gradient-diffusion hypotheses. The inherent limitations of these approximations, even in the simplest engineering flows, are well known, and detailed for example by Pope (2000). Although the effects of the assumptions are more pronounced in RANS than in LES models, where they are confined to the smaller (modelled) scales, there is clearly a need to develop higher-order models. In pollutant dispersion modelling it is also desirable to predict extreme events such as peak values or probabilities that concentrations will exceed a certain threshold. In other words, a fuller statistical description of the concentration field is required (Chatwin and Sullivan, 1993; Kristensen, 1994; Wilson, 1995; Pavageau and Schatzmann, 1999). These issues have been explored in grid turbulence and in the unobstructed atmosphere, and models capable of predicting higher-order statistics have also appeared (Yee et al., 1994; Yee and Wilson, 2000; Luhar et al., 2000; Cassiani and Giostra, 2002; Franzese, 2003; Cassiani et al., 2005a, b), but more research is necessary to extend these capabilities to cases of built-up areas.

Probability density function (PDF) methods have been developed mainly within the combustion engineering community as an alternative to moment closure techniques to simulate chemically reactive turbulent flows (Lundgren, 1969; Pope, 1985; Dopazo, 1994). Because many-species chemistry is high-dimensional and highly non-linear, the biggest challenge in reactive flows is to adequately model the chemical source term. In PDF methods, the closure problem is raised to a statistically higher level by solving for the full PDF of the turbulent flow variables instead of its moments. This has several benefits: convection, the effect of mean pressure, viscous diffusion and chemical reactions appear in closed form in the PDF transport equation. Therefore these processes are treated mathematically exactly without closure assumptions, eliminating the need for gradient-transfer approximations. The effects of fluctuating pressure, dissipation of turbulent kinetic energy and small-scale mixing of scalars still have to be modelled. The rationale is that, since the most important physical processes are treated exactly, the errors introduced by modelling assumptions for less important processes amount to a smaller departure from reality. Moreover, the higher level description provides more information that can be used in the construction of closure models.

The PDF transport equation is a high-dimensional scalar equation. Although techniques of solution based on stochastic Eulerian methods have been developed (Valiño, 1998; Mobus et al., 2001; Soulard and Sabel’nikov, 2006), Lagrangian Monte Carlo methods are a more natural choice because their computational cost increases only linearly with the dimensionality of the problem, favourably comparing to the more traditional finite difference, finite volume or finite element methods.

The numerical development of Lagrangian PDF methods has mainly centred around three distinctive approaches, all of them representing a finite ensemble of fluid particles with Lagrangian particles. A common approach is the stand-alone Lagrangian method, where the flow is represented by particles whereas the Eulerian statistics are obtained using kernel estimation (Pope, 2000; Fox, 2003). Another technique is the hybrid methodology, which builds on existing Eulerian computational fluid dynamics (CFD) codes based on moment closures (Muradoglu et al., 1999, 2001; Jenny et al., 2001; Rembold and Jenny, 2006). Hybrid methods use particles to solve for certain quantities and provide closures for the Eulerian moment equations using the particle/PDF methodology. A more recent approach is the self-consistent non-hybrid method (Bakosi et al., 2007, 2008), which also employs particles to represent the flow, and uses the Eulerian grid only to solve for inherently Eulerian quantities (such as the mean pressure) and for efficient particle tracking. Since the latter two approaches extensively employ Eulerian grids, they are particle-in-cell methods (Grigoryev et al., 2002).

The current study presents an application of the non-hybrid method to a simplified urban-scale case where pollution released from a concentrated line source between idealized buildings is simulated and results are compared to data from wind-tunnel experiments and LES.

PDF methods in atmospheric modelling mostly focus on the simulation of passive pollutants, wherein the velocity field (mean and turbulence) is either assumed or obtained from experiments (Sawford, 2004, 2006; Cassiani et al., 2005a, b, 2007). In contrast, the current model directly computes the joint PDF of the turbulent velocity, characteristic turbulent frequency and scalar concentration, extending the use of PDF methods in atmospheric modelling to represent more physics at a higher statistical level. A computed full joint PDF also has the advantage of providing information on the uncertainty originating from turbulence on a physically sound basis.

In this study the turbulent boundary layers developing along solid walls are treated in two different ways: either fully resolved or via the application of wall functions (i.e. the logarithmic “law of the wall”). The full resolution is obtained using Durbin’s elliptic relaxation technique (Durbin, 1993), which was incorporated into the PDF methodology by Dreeben and Pope (1997a, 1998). This technique allows for an adequate representation of the near-wall low-Reynolds-number effects, such as the high inhomogeneity and anisotropy of the Reynolds stress tensor and wall-blocking. Wall conditions for particles based on the logarithmic “law of the wall” in the PDF framework have also been developed by Dreeben and Pope (1997b). These two types of wall treatments are examined in terms of trade-off between computational cost and performance, addressing the question of how important it is to adequately resolve the boundary layers along solid walls in order to obtain reasonable scalar statistics.

At the urban scale the simplest settings to study turbulent flow and dispersion patterns are street canyons. Due to increasing concerns for environmental issues and air quality standards in cities, a wide variety of canyon configurations and release scenarios have been studied both experimentally (Hoydysh et al., 1974; Wedding et al., 1977; Rafailids and Schatzmann, 1995; Meroney et al., 1996; Pavageau and Schatzmann, 1999) and numerically (Lee and Park, 1994; Johnson and Hunter, 1995; Baik and Kim, 1999; Huang et al., 2000; Liu and Barth, 2002). Street canyons have a simple flow geometry, they can be studied in two dimensions, and a wealth of experimental and modelling data are available for different street-width to building-height ratios. This makes them ideal candidates for testing a new urban pollution dispersion model. We validate the computed velocity and scalar statistics with the LES results of Liu and Barth (2002) and the wind-tunnel measurements of Meroney et al. (1996); Pavageau (1996); Pavageau and Schatzmann (1999). The experiments have been performed in the atmospheric wind tunnel of the University of Hamburg, where the statistics of the pollutant concentration field have been measured in an unusually high number of locations in order to provide fine details inside the street canyon.

The rest of the paper is organized as follows. In Section 2, the exact and modelled governing equations are presented. Several statistics are compared to experimental data and large-eddy simulations in Section 3. Finally, Section 4 draws some conclusions and elaborates on possible future directions.

2 Governing equations

We write the Eulerian governing equations for a passive scalar released in a viscous, Newtonian, incompressible flow as


where , , , , and are the Eulerian velocity, pressure, constant density, kinematic viscosity, scalar concentration and scalar diffusivity, respectively. Based on this system of equations the exact transport equation that governs the one-point, one-time Eulerian joint PDF of velocity and concentration can be written as (Pope, 1985, 2000),


where and denote the sample space variables of the stochastic velocity and concentration fields, respectively, and the pressure is decomposed into its mean and fluctuation part . In Equation (4) the physical processes of advection (second term on the left), viscous diffusion (first term on the right) and transport of in velocity space by the mean pressure gradient (second term on the right) are represented mathematically exactly. The last three terms in the form of conditional expectations have to be modelled: these are respectively the effect of dissipation of turbulent kinetic energy, the effect of fluctuating pressure and the small-scale diffusion of the scalar. After appropriate modelling of the unclosed terms, Equation (4) can, in principle, be solved with a traditional numerical method. However, the high-dimensionality makes Monte Carlo methods more appealing. In particular, because Equation (4) is a Fokker-Planck equation, it can be written as an equivalent system of stochastic differential equations (SDEs) (van Kampen, 2004). We use the generalized Langevin model (GLM) of Haworth and Pope (1986) for the velocity increment, and the interaction by exchange with the conditional mean (IECM) model for the scalar. The physics and characteristics of the IECM model are discussed in detail by Fox (1996); Pope (1998) and Sawford (2004). Thus our system of SDEs that solve Equation (4) is written as


Equation (5) governs the Lagrangian particle position , and it consists of advection by the instantaneous particle velocity and molecular diffusion represented by the isotropic Wiener process , which is a given Gaussian process with zero mean and variance . The particle velocity is governed by Equation (6). The second and third terms on the right-hand side are a direct consequence of the particular Lagrangian representation of the viscous diffusion in Equation (5) (for a derivation, see Dreeben and Pope, 1997a), thus the first three terms on the right-hand side govern the particle velocity in a laminar flow mathematically exactly. The last two terms involving the tensor and jointly model the effect of pressure redistribution and anisotropic dissipation of turbulent kinetic energy. is a second-order tensor function of the mean velocity gradients , the Reynolds stresses and the rate of dissipation of turbulent kinetic energy , while is a positive constant. Note that the Wiener process appearing in both Equations (5) and (6) is the same process, i.e. the same exact series of random numbers, and is independent of the other process . The concentration of the passive scalar is governed by Equation (7), which represents the physical process of diffusion by relaxation of the particle scalar towards the velocity-conditioned scalar mean with a timescale . Note that the Eulerian statistics denoted by angled brackets are to be evaluated at the fixed particle locations , and in particle-in-cell methods, this is usually achieved using an Eulerian grid and computing ensemble averages in grid elements and/or around vertices. The energy dissipation rate is defined as


where is a model constant. We adopt the model of van Slooten et al. (1998) for the stochastic characteristic turbulent frequency


where is a source/sink term for the mean turbulent frequency


where is the production of turbulent kinetic energy, is a scalar valued Wiener process, while and are model constants. To define the micro-mixing time scale for a scalar released from a concentrated source in a geometrically complex flow domain we follow Bakosi et al. (2007, 2008) and specify the inhomogeneous as


where denotes the radius of the source, is a characteristic velocity at that we take as the Euclidean norm of the mean velocity vector, is the location of the source, and and are model constants.

2.1 Elliptic relaxation modelling of near-wall turbulence

The turbulence model GLM (represented by the last two terms of Equation (6) for the velocity increments) is in fact a family of models. A specific definition of and corresponds to a particular closure from a wealth of models, many of which can be made equivalent (at the level of second moments) to popular Reynolds-stress closures (Pope, 1994). To be able to capture the near-wall low-Reynolds-number effects on the Reynolds stresses in fully resolved boundary layers, we follow Durbin (1993) and Dreeben and Pope (1998) and specify and through the tensor


where denotes the turbulent kinetic energy and is obtained by solving the elliptic equation


where the fourth-order tensor is given by




the Reynolds stress anisotropy


and are model constants. The characteristic length scale in Equation (14) is defined by the maximum of the turbulent and Kolmogorov length scales


with , where and are model constants, and is the unit wall-normal of the closest wall element pointing outward of the flow domain. For the applied wall-boundary conditions in this fully-resolved case the reader is referred to Dreeben and Pope (1998) and Bakosi et al. (2008). More details on the inflow and outflow conditions for the mean pressure and on the wall conditions for the tensor are given in (Bakosi et al., 2008). Equation (14) was developed in conjunction with turbulent channel flow (Durbin, 1993; Dreeben and Pope, 1998). Modifications and different forms of this idea have also been proposed (Whizman et al., 1996; Dreeben and Pope, 1997a, 1998; Wacławczyk et al., 2004), and an application in channel flow with the current non-hybrid model is presented by Bakosi et al. (2007). The model to compute the joint PDF of velocity, scalar, and characteristic turbulent frequency is now complete.

2.2 Wall functions modelling of near-wall turbulence

When the wall region has sufficient resolution, and as defined by Equations (12) and (13) enable the model to adequately capture the near-wall effects on the higher-order turbulence statistics. Full resolution at the wall is strictly required in certain cases such as, for example, computations of heat transfer at walls embedded in a flow or detaching boundary layers with high adverse pressure gradients. In many other realistic simulations, however, full resolution of high-Reynolds-number boundary layers is not always possible and may not be necessary when the flow details close to walls are not important because the analysis focuses on the boundary-layer effects at farther distances. In these cases an alternative option is to model the near-wall turbulence using wall functions instead of the elliptic relaxation technique. Employing wall functions for no-slip walls provides a trade-off between the accuracy of fully resolved boundary layers and computational speed. Wall functions are widely applied in atmospheric simulations, where full wall resolution is usually prohibitively expensive even at the microscale or urban scale (Bacon et al., 2000; Lien et al., 2004). It is worth emphasizing that one of the main assumptions used in the development of wall functions is that the boundary layer remains attached, which is not always the case in simulations of complex flows. However, since wall functions are the only choice for realistic atmospheric simulations, they are still routinely employed with reasonable success.

To investigate the gain in performance and the effect on the results, we implemented the wall treatment for complex flow geometries that has been developed for the PDF method by Dreeben and Pope (1997b). In this case, the tensor is defined by the simplified Langevin model (SLM) (Haworth and Pope, 1986) and is simply a constant:


with . In line with the purpose of wall functions, boundary conditions have to be imposed on particles that hit the wall so that their combined effect on the statistics at the first grid point from the wall will be consistent with the universal logarithmic wall function in equilibrium flows, i.e. in boundary layers with no significant adverse pressure gradients. The development of boundary conditions based on wall functions rely on the self-similarity of attached boundary layers close to walls. These conditions are applied usually at the first grid point from the wall based on the assumption of constant or linear stress distribution. This results in the well-known self-similar logarithmic profile for the mean velocity. For the sake of completeness the conditions on the particles developed by Dreeben and Pope (1997b) are repeated here. The condition for the wall-normal component of the particle velocity reads


where the subscripts and denote reflected and incident particle properties, respectively. The reflected streamwise particle velocity is given by


where the coefficient is determined by imposing consistency with the logarithmic law at the distance of the first grid point from the wall, :


where is a characteristic velocity scale of the turbulence intensity in the vicinity of , defined as


with . , and are, respectively, the mean streamwise velocity, the wall-normal component of the Reynolds stress tensor and the turbulent kinetic energy at . In Equation (22) is the magnitude of the equilibrium value of the mean velocity at and is specified by the logarithmic law


where is the Kármán constant and the surface roughness parameter for a smooth wall. The friction velocity is computed from local statistics as




In Equations (20-25) the streamwise and wall-normal coordinate directions are defined according to the local tangential and normal coordinate directions of the particular wall element in question. In other words, if the wall is not aligned with the flow coordinate system then the vectors and , and the Reynolds stress tensor , need to be appropriately transformed into the wall-element coordinate system before being employed in the above equations.

The condition on the turbulent frequency is given by




which completes the description of the wall function approach.

In summary, the flow is represented by a large number of Lagrangian particles whose position , velocity , scalar concentration and characteristic turbulent frequency are governed by Equations (5), (6), (7) and (9), respectively. These equations are discretised and advanced in time by the explicit forward Euler-Maruyama method (Kloeden and Platen, 1999). The mean pressure, required in Equation (6), is obtained via a pressure projection scheme (Bakosi et al., 2008). Full wall resolution is obtained through Equations (12-18), while wall functions are applied through Equations (19-28). The pressure-Poisson and elliptic relaxation (14) equations are solved using an unstructured Eulerian grid with the finite element method. The grid is also used to track particles throughout the domain and to estimate Eulerian statistics using ensemble averaging. In practical simulations using PDF methods a few hundred particles per element is usually employed. Adequate stability can already be achieved using as little as 50–100 particles, however, 300–500 particles per elements are recommended to exploit the bin structure to compute (Bakosi et al., 2008) and to decrease the statistical error. The numerical algorithm and performance issues are detailed in Bakosi et al. (2008).

3 Modelling the street canyon

Figure 1: Geometry and Eulerian mesh (consisting of approximately 12,000 cells) for the computation of turbulent street canyon with full resolution of the wall boundary layers using elliptic relaxation. The grid is generated by the general purpose mesh generator Gmsh (Geuzaine and Remacle, 2009). The positions labelled by bold numbers indicate the sampling locations for the passive scalar, equivalent with the combined set of measurement tapping holes of Meroney et al. (1996), Pavageau (1996) and Pavageau and Schatzmann (1999). In the zoomed area the refinement is depicted, which ensures an adequate resolution of the boundary layer and the vortices forming in the corner.
# 1 2 3 4 5 6 7 8 9 10 11 12
0.5 1 1.5 2 2 2 2 2 2 2.5 2 2
2 2 2 2 1.93 1.5 1.33 1 0.67 0.5 0.33 0.17
# 13 14 15 16 17 18 19 20 21 22 23 24
4 4 4 4 4 4 4 4 4 4.5 5 5.5
0.17 0.33 0.5 0.67 1 1.33 1.5 1.93 2 2 2 2
Table 1: Concentration sampling locations at building walls and tops according to the experimental measurement holes of Meroney et al. (1996), Pavageau and Schatzmann (1999) and Pavageau (1996). See also Figure 1.
1.85 0.63 5 0.25 6 0.134 72 1.4 0.1 0.5 0.73 0.02 0.7
Table 2: Constants for modelling the joint PDF of velocity, characteristic turbulent frequency and transported passive scalar.
Figure 2: Geometry and Eulerian mesh (consisting of approximately 500 cells) for the computation of flow in a turbulent street canyon with wall functions at . The domain is stripped at no-slip walls so that it does not include the close vicinity of the wall at . The positions for sampling the scalar concentrations are the same as in Figure 1.

Street canyons are often used to study flow and pollutant dispersion patterns in urban areas. A wealth of experimental data for this simplified urban-scale setting is available from wind-tunnel and LES data (Meroney et al., 1996; Pavageau and Schatzmann, 1999; Liu and Barth, 2002), making it a natural choice to validate the current, newly developed method. We will simulate the “urban roughness” case of Meroney et al. (1996), which is a model for a series of street canyons in the streamwise direction. The simulations are performed for statistically two-dimensional flow geometry, with periodic inflow and outflow boundary conditions in the free stream above the buildings (i.e. the particles crossing the outflow boundary are re-injected at the inflow boundary). The Reynolds number based on the maximum free stream velocity and the building height was . This corresponds to based on the friction velocity and the free stream height, , if the free stream above the buildings is considered as the lower part of an approximate fully developed turbulent channel flow. The velocity-conditioned scalar mean required in Equation (7) has been computed using the general method described by Bakosi et al. (2007) using a bin structure of . After the flow has reached a statistically stationary state, time averaging is used to collect velocity statistics and a continuous scalar is released from a street level line source at the centre of the canyon (corresponding to a point source in two dimensions). Particles travelling through the source of diameter are assigned a unit source strength. The scalar field is also time averaged after it has reached a stationary state.

The simulations with the full resolution model have been made with the constants given in Table 2, using 300 particles per element. The Eulerian mesh used for this simulation is displayed in Figure 1, which shows the considerable refinement along the building walls and tops necessary to resolve the boundary layers. In this case, the high anisotropy and inhomogeneity of the Reynolds stress tensor in the vicinity of walls are captured by the elliptic relaxation technique, using Equation (14).

The simulations using wall functions were performed on the Eulerian mesh displayed in Figure 2, also using 300 particles per element. The particle-boundary conditions described in Sec. 2.2 were implemented for arbitrary geometry. Note that the first grid point where the boundary conditions based on wall functions are to be applied should not be closer to the wall than , where is the non-dimensional distance from the wall in wall units, but should still be sufficiently close to the wall to lie in the inertial sublayer (Dreeben and Pope, 1997b). Accordingly, the grid in Figure 2 only contains the domain stripped from the wall region at .

Turbulence and scalar statistics are obtained entirely from the particles that represent both the flow itself and the scalar concentration field. The Eulerian meshes displayed in Figure 1 for the full resolution and in Figure 2 for the wall function cases are used to extract the statistics, to track the particles throughout the domain and to solve the Eulerian equations, namely Equation (14) and the mean-pressure-Poisson equation in the fully resolved case and only the latter in the wall function case.

Figure 3: Velocity vectors (first row) and iso-contours of turbulent kinetic energy (second row) of the fully developed turbulent street canyon at based on the maximum free stream velocity and the building height . Left – full resolution with elliptic relaxation, right – coarse simulation with wall functions.

In Figure 3, the mean velocity vector field and the iso-contours of the turbulent kinetic energy are displayed for both fully resolved and wall function simulations. It is apparent that the full resolution captures even the smaller counter-rotating eddies at the internal corners of the canyon, while the coarse grid resolution with wall functions only captures the overall flow pattern characteristic of the flow, such as the large steadily rotating eddy inside the canyon. The turbulent kinetic energy field is captured in a similar manner. Both methods reproduce the highest turbulence activity at the building height above the canyon, with a maximum at the windward building corner. The full resolution simulation shows a more detailed spatial distribution of energy, whereas the coarse resolution of the wall-function simulation still allows tha capture of the overall pattern.

Figure 4: Dimensionless turbulent intensities (first column) and (second column) computed using full wall resolution (first row) and using wall functions (second row) at compared with the LES results (third row) of Liu and Barth (2002).

In Figure 4 the normalized turbulent intensities and are displayed for both simulation cases and compared with the LES results of Liu and Barth (2002). In the large eddy simulations the filtered momentum equations are solved by the Galerkin finite element method using brick three-dimensional elements, while the residual stresses are modelled using the Smagorinsky closure.

The full resolution simulation shows very good agreement with LES. The contour plots of correctly display two local maxima, at the windward external and at the leeward internal corners. The contour plots of show distributed high values at the building level above the canyon, along the windward internal corner and wall, and at the street level downstream of the source. By contrast, the wall-function contour plots are in general less detailed, failing to reproduce the internal maximum of , and showing a more uniform representation of .

Several wind-tunnel measurements have been carried out for a scalar released from a street level continuous line source at the centre of the canyon, providing concentration statistics above the buildings, at the walls, and inside the canyon (Meroney et al., 1996; Pavageau, 1996; Pavageau and Schatzmann, 1999). To examine the concentration values along the building walls and tops, we sampled the computed mean concentration field at the locations depicted in Figure 1 and listed in Table 1.

The excellent agreement of the results using both full resolution and wall functions with a number of experiments is shown in Figure 5. The concentration peak is precisely captured at the internal leeward corner and the model accurately reproduces the pattern of concentration along both walls including the higher values along the leeward wall.

Figure 5: Distribution of mean concentrations at the boundary of the street canyon. The experimental data are in terms of the ratio , where is the actual measured mean concentration (ppm), is the free-stream mean velocity () taken at the reference height and is the line source strength () in which denotes the scalar flow rate and is the source length. The calculation results are scaled to the concentration range of the experiments. References for experimental data: Meroney et al. (1996); , , Pavageau and Schatzmann (1999); Pavageau (1996). See also Figure 1 and Table 1 for the measurement locations.
Figure 6: Comparison of the spatial distribution of the normalized mean (left column) and variance (right column) of the scalar released at the centre of the street level. The normalization and the scaling of the calculated results are the same as in Figure 5. First row – PDF calculations with full wall resolution, second row – PDF calculations with wall functions, third row – experimental data of Pavageau and Schatzmann (1999) and fourth row – LES calculations of Liu and Barth (2002).

In Figure 6, the first two statistical moments of the concentration inside the canyon are compared with experimental data and LES. The agreement with observations indicates that both the fluid dynamics and the micro-mixing components of the model provide a good representation of the real field.

Because the one-point one-time joint PDF contains all higher-order statistics and correlations of the velocity and scalar fields resulting from a close, low-level interaction between the two fields, a great wealth of statistical information is available for atmospheric transport and dispersion calculations. As an example, in Figure 7 the cumulative distribution functions (CDF) of scalar concentration fluctuations


are depicted after time averaging at selected locations of the domain for the full resolution case. No experimental data are available to assess the distributions in Figure 7, although the irregular shape of the CDF at , , which corresponds to a multimodal PDF, is likely to be an artifact of the micro-mixing model.

The performance gain obtained by applying wall functions as opposed to full resolution was about two orders of magnitude already at our moderate Reynolds number. The gain for higher Reynolds numbers is expected to increase faster than linearly.

Figure 7: Cumulative distribution functions of scalar concentration fluctuations (left) at and (right) at using full resolution at walls.

4 Discussion

We have used an Eulerian unstructured grid, consisting of triangular element types, to estimate Eulerian statistics, to track particles throughout the domain, and to solve for inherently Eulerian quantities in conjunction with a PDF method. The boundary layers developing close to solid walls are fully captured with an elliptic relaxation technique, but can also be represented by wall functions, which use a coarser grid resolution and require significantly less particles, resulting in substantial savings in computational cost. We found that the one-point statistics of the joint PDF of velocity and scalar are well-captured by the wall function approximation. In view of its affordable computational load and reasonable accuracy, this approximation appears to hold a realistic potential for application of the PDF method in atmospheric simulations, where the natural extension of the work is the implementation of the model in three spatial dimensions.

In hybrid PDF models developed for complex chemically reacting flows, numerical treatments for boundary conditions have been included for symmetric, inflow, outflow and free-slip walls employing the ghost-cell approach common in finite volume methods (Rembold and Jenny, 2006). The representation of no-slip boundaries adds a significant challenge to the above cases. This is partly due to the increased computational expense entailed by the higher Eulerian grid resolution that is required to fully resolve the boundary layers. In addition, there is an increased complexity in specifying the no-slip particle conditions for both fully resolved and wall function representations. We presented an implementation of both approaches to treat no-slip boundaries with unstructured grids in conjunction with the finite element method. This also obviates further complications with ghost cells.

In the case of full wall resolution we employed the Lagrangian equivalent of a modified isotropisation of production (IP) model as originally suggested by Dreeben and Pope (1998). The elliptic relaxation technique, however, allows for the application of any turbulence model developed for high-Reynolds-number turbulence (Durbin, 1993; Whizman et al., 1996). The standard test case for developing near-wall models is the fully developed turbulent channel flow. In this case, we explored the simpler Rotta (1951) model, which is the Eulerian equivalent of the simplified Langevin model (SLM) in the Lagrangian framework (Pope, 1994). This is simply achieved by eliminating the term involving the fourth-order tensor from the right-hand side of Equation (14). While the SLM makes no attempt to represent the effect of rapid pressure (Pope, 2000) (in fact it is strictly correct only in decaying homogeneous turbulence), it is widely applied due to its is simplicity and robustness. Our experience showed a slight degradation of the computed velocity statistics (as compared to direct numerical simulation) using SLM for the case of channel flow. Since we experienced no significant increase in computational expense or decrease in numerical stability, we retained the original IP model.

Similarly, in the case of wall functions, several choices are available regarding the employed turbulence model. The methodology developed by Dreeben and Pope (1997b) uses the SLM, but it is general enough to include other more complex closures, such as the Haworth and Pope (1986, 1987) models (HP1 and HP2), the different variants of the IP models (IPMa, IPMb, LIPM) (Pope, 1994) or the Lagrangian version of the SSG model of Speziale et al. (1991). All these closures can be collected under the umbrella of the generalized Langevin model, by specifying its constants as described by Pope (1994). These models have been all developed for high-Reynolds-number turbulence and need to be modified in the vicinity of no-slip walls. Including them in the wall-function formulation is possible by specifying the reflected particle frequency at the wall as


instead of Equation (28). This involves the additional computation of the statistics and at , which does not increase the computational cost significantly, but may result in a numerically less stable condition since the (originally constant) parameter that appears using the SLM has been changed to a variable that fluctuates during simulation. We implemented and tested all the above turbulence models using the wall function technique. Without any modification of the model constants we found the IPMa and SLM to be the most stable, providing very similar results. Thus we retained the original (and simplest) SLM along with Equation (28).

The most widely employed closure to model the small scale mixing of the passive scalar in the Lagrangian framework is the interaction by exchange with the mean (IEM) model of Villermaux and Devillon (1972) and Dopazo and O’Brien (1974). This simple and efficient model, however, fails to comply with several physical constraints and desirable properties of an ideal mixing model (Fox, 2003). The interaction by exchange with the conditional mean (IECM) model overcomes some of the difficulties inherent in the IEM model. In this study we justify the use of the IECM model by its being more physical and more accurate, but we acknowledge that it markedly increases the computational cost (Bakosi et al., 2008).

5 Acknowledgment

It is a pleasure to ackowledge the fruitful discussions during the course of the present work with Dr. Guido Cervone at George Mason University.


  • Bacon et al. (2000) Bacon, D. P., Ahmad, N. N., Boybeyi, Z., Dunn, T. J., Hall, M. S., Lee, P. C. S., Sarma, R. A., Turner, M. D., III., K. T. W., Young, S. H., Zack, J. W., 2000. A dynamically adapting weather and dispersion model: The Operational Multiscale Environment Model with Grid Adaptivity (OMEGA). Mon. Weather Rev. 128, 2044–2076.
  • Baik and Kim (1999) Baik, J.-J., Kim, J.-J., 1999. A numerical study of flow and pollutation dispersion characteristic in urban street canyons. J. Appl. Meteorol. 38, 1576–1589.
  • Bakosi et al. (2007) Bakosi, J., Franzese, P., Boybeyi, Z., 2007. Probability density function modeling of scalar mixing from concentrated sources in turbulent channel flow. Phys. Fluids 19 (11), 115106.
  • Bakosi et al. (2008) Bakosi, J., Franzese, P., Boybeyi, Z., 2008. A non-hybrid method for the PDF equations of turbulent flows on unstructured grids. J. Comput. Phys. 227 (11), 5896–5935.
  • Cassiani et al. (2005a) Cassiani, M., Franzese, P., Giostra, U., 2005a. A PDF micromixing model of dispersion for atmospheric flow. Part I: development of the model, application to homogeneous turbulence and to neutral boundary layer. Atmos. Environ. 39 (8), 1457–1469.
  • Cassiani et al. (2005b) Cassiani, M., Franzese, P., Giostra, U., 2005b. A PDF micromixing model of dispersion for atmospheric flow. Part II: application to convective boundary layer. Atmos. Environ. 39 (8), 1471–1479.
  • Cassiani and Giostra (2002) Cassiani, M., Giostra, U., 2002. A simple and fast model to compute concentration moments in a convective boundary layer. Atmos. Environ. 36 (30), 4717–4724.
  • Cassiani et al. (2007) Cassiani, M., Radicchi, A., Albertson, J. D., 2007. Modelling of concentration fluctuations in canopy turbulence. Boundary-Layer Meteorol. 122 (3), 655–681.
  • Chatwin and Sullivan (1993) Chatwin, P. C., Sullivan, P. J., 1993. The structure and magnitude of concentration fluctuations. Boundary-Layer Meteorol. 62, 269–280.
  • Dopazo (1994) Dopazo, C., 1994. Recent developments in pdf methods. In: Libby, P. A. (Ed.), Turbulent reactive flows. Academic, New York, pp. 375–474.
  • Dopazo and O’Brien (1974) Dopazo, C., O’Brien, E. E., 1974. An approach to the autoignition of a turbulent mixture. Acta Astronaut. 1, 1239–1266.
  • Dreeben and Pope (1997a) Dreeben, T. D., Pope, S. B., 1997a. Probability density function and Reynolds-stress modeling of near-wall turbulent flows. Phys. Fluids 9 (1), 154–163.
  • Dreeben and Pope (1997b) Dreeben, T. D., Pope, S. B., 1997b. Wall-function treatment in pdf methods for turbulent flows. Phys. Fluids 9 (9), 2692–2703.
  • Dreeben and Pope (1998) Dreeben, T. D., Pope, S. B., 1998. Probability density function/Monte Carlo simulation of near-wall turbulent flows. J. Fluid Mech. 357, 141–166.
  • Durbin (1993) Durbin, P. A., 1993. A Reynolds stress model for near-wall turbulence. J. Fluid Mech. 249, 465–498.
  • Fox (1996) Fox, R. O., 1996. On velocity-conditioned scalar mixing in homogeneous turbulence. Phys. Fluids 8 (10), 2678–2691.
  • Fox (2003) Fox, R. O., 2003. Computational models for turbulent reacting flows. Cambridge University Press.
  • Franzese (2003) Franzese, P., 2003. Lagrangian stochastic modeling of a fluctuating plume in the convective boundary layer. Atmos. Environ. 37, 1691–1701.
  • Geuzaine and Remacle (2009) Geuzaine, C., Remacle, J.-F., 2009. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities (accepted for publication). Int. J. Numer. Meth. Eng.
  • Grigoryev et al. (2002) Grigoryev, Y. N., Vshivkov, V. A., Fedoruk, M. P., 2002. Numerical “particle-in-cell” methods: theory and applications. Utrecht, Boston.
  • Haworth and Pope (1986) Haworth, D. C., Pope, S. B., 1986. A generalized Langevin model for turbulent flows. Phys. Fluids 29 (2), 387–405.
  • Haworth and Pope (1987) Haworth, D. C., Pope, S. B., 1987. A pdf modeling study of self-similar turbulent free shear flows. Phys. Fluids 30 (4), 1026–1044.
  • Hoydysh et al. (1974) Hoydysh, W. G., Griffiths, R. A., Ogawa, Y., 1974. A scale model study of the dispersion of pollution in street canyons. In: 67th Annual Meeting of the Air Pollution Control Association, Denver, Colorado, APCA Paper No. 74-157.
  • Huang et al. (2000) Huang, H., Akutsu, Y., Arai, M., Tamura, M., 2000. A two-dimensional air quality model in an urban street canyon: Elevation and sensitivity analysis. Atmos. Environ. 34, 689–698.
  • Jenny et al. (2001) Jenny, P., Pope, S. B., Muradoglu, M., Caughey, D. A., 2001. A hybrid algorithm for the joint PDF equation of turbulent reactive flows. J. Comput. Phys. 166, 218–252.
  • Johnson and Hunter (1995) Johnson, G. T., Hunter, L. J., 1995. A numerical study of dispersion of passive scalar in city canyons. Boundary-Layer Meteorol. 75, 235–262.
  • Kloeden and Platen (1999) Kloeden, P. E., Platen, E., 1999. Numerical Solution of Stochastic Differential Equations. Springer, Berlin.
  • Kristensen (1994) Kristensen, L., 1994. Recurrence of extreme concnetrations. In: S.-V. Grining, M. M. Millán (Eds.), Air Pollution Modelling and Its Application. X. Plenum Press, New York.
  • Lee and Park (1994) Lee, I. Y., Park, H. M., 1994. Parameterization of the pollutant transport and dispersion in urban street canyons. Atmos. Environ. 28, 2343–2349.
  • Lien et al. (2004) Lien, F. S., Yee, E., Cheng, Y., 2004. Simulation of mean flow and turbulence over a 2d building array using high-resolution CFD and a distributed drag force approach. J. Wind Eng. Ind. Aerod. 92 (2), 117–158.
  • Liu and Barth (2002) Liu, C.-H., Barth, M. C., 2002. Large-eddy simulation of flow and scalar transport in a modeled street canyon. J. Appl. Meteorol. 41 (6), 660–673.
  • Luhar et al. (2000) Luhar, A. K., Hibberd, M. F., Borgas, M. S., 2000. A skewed meandering-plume model for concentration statistics in the convective boundary layer. Atmos. Environ. 34, 3599–3616.
  • Lundgren (1969) Lundgren, T. S., 1969. Model equation for nonhomogeneous turbulence. Phys. Fluids 12 (3), 485–497.
  • Meroney et al. (1996) Meroney, R. N., Pavageau, M., Rafailidis, S., Schatzmann, M., 1996. Study of line source characteristics for 2-d physical modelling of pollutant dispersion in street canyons. J. Wind Eng. Ind. Aerod. 62 (1), 37–56.
  • Mobus et al. (2001) Mobus, H., Gerlinger, P., Bruggemann, D., 2001. Comparison of eulerian and lagrangian monte carlo pdf methods for turbulent diffusion flames. Combust. Flame 124, 519–534.
  • Muradoglu et al. (1999) Muradoglu, M., Jenny, P., Pope, S. B., Caughey, D. A., 1999. A consistent hybrid finite-volume/particle method for the PDF equations of turbulent reactive flows. J. Comput. Phys. 154, 342–371.
  • Muradoglu et al. (2001) Muradoglu, M., Pope, S. B., Caughey, D. A., 2001. The hybrid method for the PDF equations of turbulent reactive flows: consistency conditions and correction algorithms. J. Comput. Phys. 172, 841–878.
  • Pavageau (1996) Pavageau, M., 1996. Concentration fluctuations in urban street canyons – groundwork for future studies. Tech. rep., Meteorological Institute of the University of Hamburg.
  • Pavageau and Schatzmann (1999) Pavageau, M., Schatzmann, M., 1999. Wind tunnel measurements of concentration fluctuations in an urban street canyon. Atmos. Environ. 33, 3961–3971.
  • Pope (1985) Pope, S. B., 1985. PDF methods for turbulent reactive flows. Prog. Energ. Combust. 11, 119–192.
  • Pope (1994) Pope, S. B., 1994. On the relationship between stochastic Lagrangian models of turbulence and second-moment closures. Phys. Fluids 6 (2), 973–985.
  • Pope (1998) Pope, S. B., 1998. The vanishing effect of molecular diffusivity on turbulent dispersion: implications for turbulent mixing and the scalar flux. J. Fluid Mech. 359, 299–312.
  • Pope (2000) Pope, S. B., 2000. Turbulent flows. Cambridge University Press, Cambridge.
  • Rafailids and Schatzmann (1995) Rafailids, S., Schatzmann, M., 1995. Physical modelling of car exhaust disperision in urban street canyons. In: Proc. 21st Int. Meeting on Air Pollution Modelling and Its Applications, Baltimore, Nov. 6-10.
  • Rembold and Jenny (2006) Rembold, B., Jenny, P., 2006. A multiblock joint pdf finite-volume hybrid algorithm for the computation of turbulent flows in complex geometries. J. Comput. Phys. 220, 59–87.
  • Rotta (1951) Rotta, J. C., 1951. Statistiche theorie nichthomogener turbulenz. Z. Phys. 129, 547.
  • Sawford (2004) Sawford, B. L., 2004. Micro-mixing modelling of scalar fluctuations for plumes in homogeneous turbulence. Flow Turbul. Combust. 72, 133–160.
  • Sawford (2006) Sawford, B. L., 2006. Lagrangian modeling of scalar statistics in a double scalar mixing layer. Phys. Fluids 18 (8), 085108.
  • Soulard and Sabel’nikov (2006) Soulard, O., Sabel’nikov, V. A., 2006. Eulerian Monte Carlo method for the joint velocity and mass-fraction probability function in turbulent reactive gas flows. Combust. Explo. Shock. 42, 753–762.
  • Speziale et al. (1991) Speziale, C. G., Sarkar, S., Gatski, T. B., 1991. Modelling the pressure-strain correlation of turbulence: an invariant dynamical systems approach. J. Fluid Mech. 227, 245–272.
  • Valiño (1998) Valiño, L., 1998. A field Monte Carlo formulation for calculating the probability density function of a single scalar in a turbulent flow. Flow Turbul. Combust. 60, 157–172.
  • van Kampen (2004) van Kampen, N. G., 2004. Stochastic processes in physics and chemistry, 2nd Edition. North Holland, Elsevier B.V., Amsterdam, The Netherlands.
  • van Slooten et al. (1998) van Slooten, P. R., Jayesh, Pope, S. B., 1998. Advances in PDF modeling for inhomogeneous turbulent flows. Phys. Fluids 10 (1), 246–265.
  • Villermaux and Devillon (1972) Villermaux, J., Devillon, J. C., 1972. Représentation de la coalescence et de la redispersion des domaines de ségrégation dans un fluide par un modèle d’interaction phénoménologique. In: Proceedings of the 2nd Intl. Symp. on Chemical Reaction Engineering. Elsevier, New York, pp. 1–13.
  • Wacławczyk et al. (2004) Wacławczyk, M., Pozorski, J., Minier, J.-P., 2004. Probability density function computation of turbulent flows with a new near-wall model. Phys. Fluids 16 (5), 1410–1422.
  • Wedding et al. (1977) Wedding, J. B., Lambert, D. J., Cermak, J. E., 1977. A wind tunnel study of gaseous pollutants in city street canyons. Air Pollut. Control Assoc. J. 27, 557–566.
  • Whizman et al. (1996) Whizman, V., Laurence, D., Kanniche, M., Durbin, P. A., Demuren, A., 1996. Modeling near-wall effects in second-moment closures by elliptic relaxation. Int. J. Heat Fluid Fl. 17, 255–266.
  • Wilson (1995) Wilson, D. J., 1995. Concentration fluctuations and averaging time in vapor clouds. Center for Chemical Process Safety. American Institute of Chemical Engineers, New York.
  • Yee et al. (1994) Yee, E., Chan, R., Kosteniuk, P. R., Chandler, G. M., Biltoft, C. A., Bowers, J. F., 1994. Incorporation of internal fluctuations in a meandering plume model of concentration fluctuations. Boundary-Layer Meteorol. 67 (1–2), 11–39.
  • Yee and Wilson (2000) Yee, E., Wilson, D. J., 2000. A comparison of the detailed structure in dispersin tracer plumes measured in grid-generated turbulence with a meandering plume model incorporating internal fluctuations. Boundary-Layer Meteorol. 94 (2), 253–296.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description