Abstract
Field studies have shown that local climate is strongly influenced by urban structures. This influences both energy consumption and the pedestrian comfort. It is thus useful to be able to simulate the urban environment to take these effects into account in building and urban design. But for computational reasons, conventional computational fluid dynamics (CFD) codes cannot be used directly on a grid fine enough to resolve all scales found in a city. For this, we use mesoscale models, variants of CFD codes in which the 3D conservation equations are solved on grids having a resolution of a few kilometers. At this resolution, the effects of subgrid scales need implicit representations. In other words, phenomena such as momentum and energy exchanges averaged over the mesoscale grid contribute necessary sources/sinks to the corresponding equations. Such spatial averaging results in additional terms called dispersive fluxes. Until now these fluxes have been ignored. To better understand these fluxes, we have conducted large eddy simulations (LESs) over an array of cubes for different inter-cube spacings. The study shows that these fluxes are as important as the turbulent fluxes and exhibit trends which are related to the eddy formations inside the canopies.
1. Introduction
The large and continuous variety of scales present in the atmospheric flow over a city generates an intrinsic difficulty in the numerical treatment of the atmospheric conservation equations. From scaling considerations, the ratio between the smallest flow scale and the characteristic length scale is approximately proportional to [1], where is the Reynolds number. This means that a 3D representation of the planetary boundary layer resolving all scales will require about grid cells. This number is far from being handled conveniently nowadays or in the near future by any computing device. Therefore, the transport phenomena over larger distances (in our case covering a city and its bounding context) must be handled by mesoscale atmospheric codes with spatial resolutions of a few kilometers either by subgrid urban parameterizations [2, 3] or by a coupling with a higher resoltion microscale model [4, 5]. In the former approach, the mesoscale codes cannot “see” buildings explicitly. Yet buildings and urban land use significantly impact the micro- and mesoscale flow, altering the wind, temperature, and turbulence fields and radiation exchanges [6, 7]. Since mesoscale numerical models do not have the spatial resolution to directly simulate the fluid dynamics and thermodynamics in and around urban structures, urban parameterizations are needed to approximate the drag, heating, and enhanced turbulent kinetic energy (tke) produced and dissipated by the subgrid-scale urban elements. The drag forces offered by the buildings as well as the heat transfer characteristics are a function of the local velocity field. Local turbulent fluxes, dispersive fluxes (generally ignored in mesoscale models), and drag coefficients can significantly impact the exchange of mass, momentum, and energy. However, a mesoscale code as described earlier, does not have the spatial resolution to generate the profiles of these quantities. Several attempts have been made in the recent past to estimate the velocity profiles in the urban canopy. These so-called urban canopy models are based on either single-layer [8] or multilayer considerations of the canopy [9]. In almost all these models, the dispersive fluxes which result from the averaging of the governing equation in the horizontal plane are either ignored or assumed to behave in the same way as the turbulent fluxes. However, recent work [10, 11] tends to confirm that these stresses can be significant and sometime comparable to the turbulent stresses themselves and may behave differently. However, most of these findings were based on simulations conducted using Reynolds Averaged Navier Stokes (RANS) codes whose validity for this kind of application is somewhat questionable [12]. Moreover, in the study we wanted to compare the magnitude of the dispersive fluxes and the turbulent fluxes. A RANS model only gives the modelled turbulent fluxes against which the dispersive fluxes cannot be compared. In an attempt to resolve this issue, we therefore employ an LES code which is capable of resolving the turbulent fluxes.
An LES “resolves” the large-scale fluid motions and “models” the subgrid-scale motions by filtering the Navier Stokes equations. When unsteady RANS methods are used, it is implicitly assumed that there is a fair degree of scale separation between the large timescale of the unsteady flow features and the time scale of the genuine turbulence. However, in reality it is hard to find an evident time scale gap for many turbulent flows. Furthermore, RANS generally does not intend to capture most of the genuinely turbulent fluctuation information. A RANS approach thus has obvious weaknesses and poses serious uncertainties in flows for which large-scale organized features dominate, such as flows around building like obstacles. Against this background, it may be argued that the use of a sound although computationally expensive LES approach is fully justified. Although a city might not be well represented by a regular array of buildings, this is nevertheless a sound pragmatic starting point because these shapes are the basic building blocks of any city and also because there is good availability of data for validation purposes. Thus, the study of the flow over a matrix of cubes (resembling an array of buildings) can provide some fundamental understandings of the various physical phenomena involved in the flow through an urban area. Various characteristics like the vortex shedding, flow separation and velocity profiles have been experimentally investigated for such problems in the past [13, 14].
In this study, all the simulations have been conducted with the Smagorinsky model because of its numerical simplicity and stability. It has also provided excellent results for the case we are interested in, when compared against the experimental data [13]. Since, we are interested in general behaviors rather than results for a particular point inside the domain, we present the spatially averaged profiles of the velocity, turbulent fluxes, and dispersive fluxes and explain their nature especially that of the dispersive fluxes.
2. Governing Equations for LES
2.1. Transport Equations
In LES, only the large-scales are explicitly resolved by the numerical grid while the smaller ones are represented by a subgrid-scale model. The motivation for this approach is that the large-scale vortices are dominated by geometrical constraints and boundary conditions. Due to turbulent transport phenomena, these vortices pass their kinetic energy onto smaller scales while the orientation of the initial vortices gets lost during this energy cascade. Therefore, the small-scale turbulence is expected to be isotropic without any preferred orientation and should consequently be much easier to model than the whole spectrum of turbulence. Starting with the governing equations for an incompressible three-dimensional (3D) unsteady flow field, we apply a top-hat filter function to separate large- and small-scale motions leading to the filtered equation set where are the filtered velocity components, is the filtered pressure, denotes the filtered strain-rate tensor, and denotes the molecular viscosity. The correlation within the convective term () is a priori unknown and has to be modeled.
The most common way is to rewrite this term into , where is the unresolved stress resulting from the subgrid-scale contribution and needs to be modelled by an appropriate subgrid-scale (SGS) model. The additional stresses are split into an anisotropic part (where is the eddy viscosity) and an isotropic part which is added to the pressure , leading to the LES equation set which forms the basis of this investigation:
2.2. Numerics
For conducting large eddy simulations, we used the Transat code [15] which is based on a finite volume discretization. It solves for mass, momentum, and heat transport in both single- and two-phase flow conditions and provides the option of using the reynolds averaged or unsteady turbulence modelling (LES or direct numerical simulation). For the present simulation, we used LES because of the accuracy we needed for a better understanding of the flow. A Runge Kutta 3rd-order scheme is used for time integration while the convective schemes used for density, velocity and temperature were HYBRID, CDS (central differencing scheme), and HLPA (hybrid linear parabolic approximation), respectively. A preconditioned (multigrid) GMRES (generalized minimal residual method) pressure solver is used for solving the acoustic equation. The Standard Smagorinsky model is used to simulate the effects of the subgrid scales on the flow. Although there are more accurate models available, the established accuracy of prediction of this particular model has proved to be sufficient for our purpose.
2.3. Subgrid-Scale Modelling
An eddy-viscosity-based model has been employed in the computations presented in this paper, where the anisotropic part of the SGS stress is modelled using The eddy viscosity is determined using Smagorinsky’s expression [16] with and , determined using an explicit box filter of width twice the mesh size in wall-parallel planes, together with averaging in the spanwise direction and relaxation in time with a factor of . The model coefficient is taken to be equal to . The dynamic Smagorinsky model (DSM) of [17], with the modification of [18], could also be applied here, but the required averaging of the model coefficient (which will now depend on the resolved invariant ) is made difficult by the absence of a clear homogeneous averaging direction. The near-wall behavior of the model is such that it yields an eddy viscosity that reduces as the wall is approached, using an explicit damping following the van Driest relationship [19]: where is the distance from the wall in viscous wall units for which is the friction velocity.
3. Governing Equations for Mesoscale Model
Since we intend to interpret the results from the LES simulations in a mesoscale modelling context, it is worth mentioning the standard Reynolds Averaged Navier Stokes equations that form the basis of a large variant of existing Mesoscale Codes. The mass and momentum equations can be mathematically represented as: where is the mean part of the velocity, are the fluctuation of velocities in time, the mean pressure, the Reynolds stresses and the source term. Equation (6) when averaged over space takes the following form: where the angular bracket is the space averaging operator. Such space averaging as we will show later in the paper results in an additional term called the dispersive flux.
4. Geometric Description and Test Cases
For the purpose of determining the validity of the LES model for our study, we used the experimental results of [13, 14] who carried out detailed measurements of the mean flow and turbulence characteristics using Doppler anemometry throughout an array of cubes. The experimental setup consisted of 250 cubes, each of height placed at a distance of from their neighbouring cubes in an aligned configuration (Figure 1(a)) consisting of 25 rows of 10 columns. The depth of the plane channel was . Due to the high computation cost associated with LES, our numerical simulations were based on a domain of (streamwise × spanwise × height) with the cube located at the center (consistent with the experimental setup) and with periodic boundary conditions in the streamwise direction. Since a symmetric boundary condition is inappropriate for the instantaneous velocity field, a periodic boundary condition was also applied across the pair of vertical boundary planes of the flow domain in the spanwise direction. For the top and bottom walls of the channel, as well as for the surfaces of the cube, no-slip and impermeability conditions were specified (for the tangential and wall normal velocity components, resp.). In accordance with the experiment, the Reynolds number for the simulation was , based on the mean bulk velocity and the height of the cube. The domain was discretized into nodes in each direction, with the grids being preferentially fine near the wall surfaces. Although no grid independence test was attempted, it should be pointed out that the resolution we used in this simulation is finer than that used in a similar simulation reported in another study [12]. We refer to this particular case () which has been validated against the experimental data, as Case .
(a) Top view of the domain
(b) Locations of investigation
Four more simulations were conducted for the domain sizes of , , , and corresponding to ( is the inter-cube spacing and is the cube height) ratios of , , , and . These we refer to as Cases , , , and . The number of nodes used in all these cases was out of which nodes have been used to represent the cube. It should be noted that Case corresponds to , which is the same as that for Case . Indeed the only difference between these two cases is the type and placement of the upper boundary: in Case a free slip boundary condition is applied at a height of . Although, we have not been able to directly validate Cases because of the unavailability of experimental data, we found that the profiles of velocities and stresses inside the canopy predicted in Case are almost identical to those obtained in Case A, so that we may have a reasonably good degree of confidence in our own results. The volume between to in Case was meshed with a uniform grid of dimensions similar to those of the topmost level in the mesh of Case A. For Cases , , and , the domain size was reduced but the same number of nodes was used, so resulting in finer resolutions in Cases , , and as compared to Case . The profiles of mean velocity and turbulence statistics were obtained using a time-averaging procedure. After carrying out the simulation for several large-eddy turnover times to ensure that the final time-averaged results were independent of the initial conditions, we averaged the instantaneous quantities over 20,000 time steps. The corresponding averaged quantities were compared with those that were similarly obtained for 40,000 time steps. Very little difference was observed between the two cases implying statistical convergence. However, for the statistical convergence of the dispersive fluxes, 200,000 time steps were required.
5. Results
5.1. Validation of LES Model
5.1.1. Velocity Profiles
For the validation of the simulation results we used the experimental data of [13]. We present results, on the mid-plane at different positions (), Figure 1(b). In Figure 2, one can see that the normalized velocity component remains positive over the cube, however, inside the spanwise canyons (between two cubes) the streamwise velocity component is negative implying a reverse flow or a strong vortex formation. Figure 3 presents the horizontal profiles of the mean streamwise velocity on the plane at at the same locations. Again one can see that the normalized streamwise velocity is positive in the streamwise canopy. However, the profiles inside the spanwise canopy is negative again due to the vortex formed between the cubes. Figure 4 displays predictions of the horizontal profiles of the mean spanwise velocity on the plane at at the same locations. A nonzero value of spanwise velocity represents recirculation, and as the figure shows this is more dominant in the vicinity of the cube. The agreement between LES results and the experimentation is generally very good. Indeed the prediction of the velocity profile () in general is very much acceptable.
5.1.2. Stress Profiles
Figures 5(a) and 5(b) present the Reynolds normal stress , in the plane at and on the plane at at five selected locations of , , , , and . From these results, it is clear that the normal stress is always at a maximum near the walls (top or side). These peaks correspond to the generation and development of thin intense vertical and horizontal shear layers along the roof and side walls of the cube, respectively. As is very clear from Figures 5(a) and 5(b), the values of the streamwise Reynolds normal stress predicted by LES are in very good agreement with the experimental data. However, the results for the spanwise Reynolds stresses (Figures 6(a) and 6(b)) are not very encouraging. While the predicted values compare well for , inside the canopy the magnitude is underpredicted. At the moment, this is an unresolved issue. Furthermore, profiles of Reynolds shear stress on the horizontal plane at half cube height (Figure 7) are in good agreement at most locations but for the agreement is not good specially inside the canopy. A more detailed physical explanation of the results can be found in the paper by [12]. In general, one can say that the LES results are in good agreement with the experimental data, and its use for further investigation of similar cases is reliable.
(a)
(b)
(a)
(b)
5.2. Effects of the Change in Street Width to Building Height Ratio on the Spatially Averaged Quantities
Since our numerical simulations (LESs) were conducted at a very high-resolution detailed profiles of the velocity field, pressure, stresses, and turbulent kinetic energy have been obtained. However, we are interested in spatially averaged quantities that can be representative of a mesoscale grid. Assuming that the values computed by the LES at every grid point are also representative of the volume average of the corresponding grid volume, we apply (7) to the domain on which LES was conducted. The equation for the streamwise velocity component can be expanded to the following form:
Since we are looking for the steady-state equation, we can neglect the first term on the left-hand side. The second term on the LHS, using the flux divergence theorem, can be written as: Here, is the volume of air over which the average is performed. is the surface delimiting the volume over which the average is performed. is the component of the normal entering in the surface ( in this case because the derivative is respect to ). For the horizontal surfaces, the value of is zero. There are two types of vertical surfaces: those at the boundaries of the domain and those delimiting the obstacle. For the surfaces at the boundary of the domain, since we have periodic boundary conditions the contribution is zero. Over the surfaces of the obstacle, the velocity is zero. So, the second term on the right-hand side of (8) is zero. Similarly, the third terms is also zero. Similarly, the first and second term on the RHS can be neglected. Because the flow is turbulent, the viscous terms (fifth, sixth, and seventh) can be neglected as well. Thus, we are left with the simplified equation: Split and , where , are spatially averaged velocity components in the direction of flow and vertical directions and and are fluctuation in space. Using these expressions, one gets Since (11) reduces to . Introducing this in (10), one gets the following equation: where the first and second terms on LHS are the gradient of turbulent and dispersive fluxes, respectively, in the vertical direction and the third term is the gradient of pressure in the flow direction. To study the vertical profiles of the streamwise velocity, turbulent stress, and dispersive stress, we evaluate these quantities from the result obtained form the simulation using (13) through (16): where for the blocked regions and , , are the indexes in the streamwise, spanwise, and vertical directions. As is clear from (13) through (16) the averaging is performed over horizontal planes at different heights.
5.2.1. Mean Velocity
Figure 8 gives the spatially averaged streamwise velocity profile. This appears to be logarithmic above the canopy but inside the canopy it can be strongly affected by adjacent cubes, as can be seen comparing the cases corresponding to different ratios. In the case of , the profile inside the canopy also takes a logarithmic profile because the flow has sufficient time and space to redevelop within the wide canopy. However, as ratio decreases the profile starts to deviate from the normal logarithmic profile, becoming linear for .
5.2.2. Turbulent Stresses
The mesoscale unresolved fluxes can be split into two components: the turbulent part and the dispersive part. The vertical profile as shown in Figure 9 of the turbulent stresses is negative throughout the profile (there is a downward transfer of momentum) implying that the flux is a downgradient. Within the canopy, the turbulent fluxes decrease with height until the top of the cube where a minima is seen. Above the canopy the magnitude increases in height to a height of . These fluxes are then absent above this height, so that there is very little sign of turbulence within this region. Also noteworthy are the linear profiles of the turbulent stress both inside and above the canopy with negative and positive slopes, respectively. These turbulent stresses, which are actually an indication of the transport property of turbulence, decrease with ratio, indicating that an area of widely spaced cubes can experience more turbulence because of greater penetration of eddies within these streets. Conversely, very narrow streets will experience little turbulence because of low eddy penetration. However, there is a need for more exhaustive data analysis to support such a generalization.
5.2.3. Dispersive Stresses
In most mesoscale models, the dispersive stresses are neglected. In order to study and understand the behavior, these stresses are plotted in Figure 10. From this it is clear that these stresses can be significant and in some cases comparable to turbulent stresses (Figure 9). These dispersive stresses are absent near the bottom wall but increase or decrease (depending upon the ) to attain a maxima or minima at half the cube height. Above the canopy these stresses vanish. As with turbulent stresses they reduce to zero above . Of particular interest is that these dispersive stresses do not exhibit a regular trend when expressed as a function of ratio. When the cubes are wide apart (), the dispersive fluxes inside the canopy are negative but upon decreasing the inter-cube spacing beyond a certain point a switch to a positive flux is experienced. The physical explanation of the dispersive stresses lies in the coherent vortex formed in the canyon.
In order to better understand the behavior of these dispersive fluxes contours of vertical velocity (at the top of the canopy), horizontal velocity and spatially averaged dispersive flux for different are plotted side by side, as shown in Figure 12. Since dispersive flux is defined as , its sign will also depend on and . Here, is always positive (Figure 8) and inside the canopies is mostly negative implying that will mostly be positive. Thus, the sign of the dispersive flux will depend on the sign of . Also, since is negligible one can conclude that the sign of will be opposite to that of . Now let us consider each of our cases in turn.
Case 1 . In Figure 11(a), one can see that a jet of fluid impinges the top of the windward side of the cube and is deflected downward, resulting in two clockwise rotating vortices one above the other with the stronger one at the top. These vortices are eccentric with their “eyes” shifted toward the leeward side of the cube. The formation of such eccentric vortices results in negative and in most of the regions inside the canopy (Figure 12(a)). However, there are regions near the bottom and top of the leeward side of the cube where there is a positive vertical velocity. At the same time, there are regions near the mid-plane inside the canopy where the vertical velocity is mostly negative. The existence of a mixture of regions of positive and negative vertical velocity fields results in the cancelation of the dispersive fluxes during the averaging over horizontal planes near the top and bottom planes. Near the mid-plane there is no such canceling, and hence a maxima of dispersive flux is obtained.
| (a) |
| (b) |
| (c) |
| (d) |
| (a) |
| (b) |
| (c) |
| (d) |
Case 2 . An increase of the ratio from 1 to 1.5 results in the shift of the eye of the primary vortex towards the windward side of the cube. The primary vortex in this case is concentric (Figure 11(b)). The formation of such a nearly concentric vortex results in positive vertical velocities near the leeward side of the cube and negative velocities near the windward side as show in Figure 12(b). On any horizontal plane below there is greater flux injection (negative vertical velocity) and less ejection (positive vertical velocity). However, above the situation is exactly the opposite. At the mid-plane itself both the ejection and injection balance each other resulting in a cancelation of dispersive fluxes. As explained earlier, the profile of the dispersive flux depends on the sign of the vertical velocity, so that in this particular case the dispersive fluxes are positive below and negative above .
Cases 3 and 4 . For both and , there is a clear demarcation between the zones of positive and negative vertical velocities. Injection from the top takes place near the windward side and ejection near the leeward side of the cube. This is the result of a strong large vortex formed inside the canopy as a consequence of the strong shear forces at the top. However, this behavior changes significantly when is increased to 2 and 3. Contrary to the observations in Cases 1 and 2, in these cases the injection occurs in the middle of the top plane just above the canopy (Figures 11(c) and 11(d)) and the jet impinges not on the top of the cube but on the lower half of the cube. This situation leads to the formation of several tilted vortices resulting in strong ejection at both the leeward and the windward sides of the cube. Also, the tilt in the vortices, which is a result of the injection in the middle, results in positive vertical velocities in most regions and hence negative dispersive fluxes (Figures 12(c) and 12(d)).
In all the cases , the vertical velocity above the cube and the canopy is positive and hence it always resulted in a negative dispersive flux. Another thing to be observed in Figures 8 and 9 is that the effect of the cube can be felt up to a height of to , which is consistent with the observations from field experiments [20].
6. Conclusion
LES with standard Smagorinsky’ model was used to compute a fully developed turbulent flow over a matrix of cubes. Detailed comparsions between the numerical predictions obtained with the LES and the corresponding experimental data of [13] were conducted. Later on the numerical data generated was used to study the spatially averaged profiles of the velocity, turbulent flux, and dispersive flux. The results of this investigation allow the following conclusions to be drawn. (i) Validation of the numerical results: qualitatively, the profiles of mean velocities and Reynolds stresses, the latter including , , and , are generally well represented by the LES model. The greatest discrepancy between the predictions and observations was for within the street canyon of the obstacle array. The underestimation of will lead to an underestimation of the turbulent kinetic energy and hence must be kept in mind. Overall the numerical predictions were very good for the job we are interested in. (ii)Spatially averaged quantities: we then carried out the tests for an array of cubes with different inter-cube spacings. The results that were obtained from LES were spatially averaged to derive information useful for urban mesoscale simulations. It was evident that the profiles of turbulent flux and dispersive flux go to zero at nearly three times the cube height, a fact which has been observed in various field experiments [20]. Of particular significance has been the use of the results to explain the behavior of dispersive fluxes. It was observed that for widely spaced array of cubes these fluxes were negative (same sign as the turbulent fluxes). This implies that these fluxes can be modeled in the same way as the turbulent fluxes. However, for high packing density, these profiles started assuming a negative profile which may result in the canceling of other sources in the equation and therefore need to be modelled differently.
However, it must be stressed here that the conclusions drawn from this study are only valid for a regular array of cubes in a neutral atmosphere. In order to generalize such conclusion, more numerical and wind tunnel experiments are required. On a related note, this work leads to pose some interesting scientific questions, such as how much complexity must be added to produce a configuration that gives spatially averaged values similar to those of a real city? Or in other words, which is the simplest configuration that represents a real city? Which combination of parameters (building heights, building shapes, building width, street widths, etc.) is sufficient to characterize city morphology?