- Email: [email protected]
Modelling and simulations of turbulent flows in urban areas with vegetation
J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
Contents lists available at ScienceDirect
Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia
Modelling and simulations of turbulent flows in urban areas with vegetation Saša Kenjereš a,n, Benjamin ter Kuile b a Transport Phenomena Section, Department of Chemical Engineering, Faculty of Applied Sciences and J.M. Burgers Centre for Fluid Dynamics, Delft University of Technology, Julianalaan 136, 2628 BL Delft, The Netherlands b UFlows, Ezelsveldlaan 153, 2611 DM Delft, The Netherlands
art ic l e i nf o
a b s t r a c t
Article history: Received 19 June 2012 Received in revised form 3 September 2013 Accepted 10 September 2013 Available online 30 October 2013
A comparative assessment of the Reynolds-Averaged Navier–Stokes (RANS) two-equation k ɛ model, with different proposals for a proper treatment of the vegetation (trees) effects in urban areas, is reported here. The most common variants proposed in the literature based on a two-equation k ɛ eddy viscosity model are revisited. An extensive model validation is performed on generic cases with homogeneous and non-homogeneous vegetation and compared with available experimental data. It is shown that the minimal level of the proper modeling of the non-homogeneous (i.e. with the vertically variable leaf area density) vegetation requires a model with additional four terms in the equations for the turbulence kinetic energy and its dissipation rate. The proposed version of the model is then applied to a computational study of flow over the Delft University of Technology site with different scenarios of the vegetation islands. By performing series of simulations, some practical guidelines are suggested to reduce intensity of the side-wind gust within the campus. In conclusion, it is demonstrated that the adopted turbulence model for non-homogeneous vegetation, in combination with the passive element approach for buildings, proved to be a physically accurate and numerically robust method. The model is recommended for future use in simulating turbulent flows in complex urban areas with vegetation. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Environmental flows Urban canopy Trees Vegetation Turbulence modeling RANS-CFD simulations
1. Introduction Severe air pollution in urban areas is usually caused by the combined effects of intensive traffic/industrial emissions and critical meteorological conditions (windless days, strong stratification, temperature inversion capping, etc.). Since the origin of all environmental pollution is in complex spatial and temporal flow, turbulence and dispersion interactions, computer simulations can be powerful tools for estimating the impacts of different parameters on pollution. Air flow and pollutant dispersion at environmental scales are characterized by high values of Reynolds numbers (indicating highly turbulent flows) and complex geometries (topography of terrain simulated, arrays of buildings and streets). Obviously, direct numerical resolving of all flow and turbulence scales is beyond the reach of present computational power, so mathematical models for scales not resolved by underlying numerical grid (subgrid/subscale contributions) must be introduced. In our previous work, we addressed modeling and simulations of some important aspects of the environmental flows, including the turbulent dispersion (passive scalars) in street canyons or over
n
Corresponding author. Tel.: þ 31 152783649. E-mail addresses: [email protected], [email protected] (S. Kenjereš).
0167-6105/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jweia.2013.09.007
wavy-walls in the laboratory conditions and some real-scale applications dealing with the time-dependent (diurnal cycles) simulations of flow, turbulence and dispersion over a complex terrain orography, e.g. Kenjereš and Hanjalić (2002, 2009), Kenjereš et al. (2005), Hanjalić and Kenjereš (2005, 2006, 2008), Kuhn et al. (2010), and Wagner et al. (2011). This is a part of the development of an integral in-house computational code, which is able to simulate the laboratory scale and real-scale environmental flows with advanced numerical methods and turbulence models. The latter also include a recently developed seamless hybrid RANS/LES approach (Kenjereš and Hanjalić, 2006). In this paper, special focus will be placed on modeling and validation of effects of vegetation (trees) on flow and turbulence in urban areas, in the framework of the two-equation eddy-viscosity models. Such models can be used for simulations of urban areas covered by significant portions of green areas (parks) or where trees are distributed along streets. These models can also be used to assess effects of the urban vegetation on dispersion of the traffic exhaust gases in the street canyons, which is an important prerequisite to achieve a sustainable development of urban areas (Gromke et al., 2008; Balczo et al., 2009). The proper modeling of vegetation effects on flow, turbulence and turbulent dispersion represents a challenging task. The flow reorganization inside and above an area covered by vegetation involves complex transfer of the turbulence
44
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
kinetic energy (Raupach and Thom, 1981; Raupach and Shaw, 1982; Wilson, 1988; Wilson et al., 1998; Ayotte et al., 1999; Dupont and Brunet, 2009; Yi, 2008; Belcher et al., 2012). Due to a blocking effect (the flow is slowing down inside the vegetation area due to the drag force) strong velocity gradients can be created producing a significant source of the turbulence kinetic energy. This blocking can also lead to the formation of long recirculative regions (wakes) similar to the flows over bluff bodies. On the other hand, inside the vegetation area, leaves act as turbulence suppressors – enhancing its dissipation rate. This is done by generating small scale eddies, without undergoing the entire process of cascade interactions (a so-called shortcut in the turbulence cascade process, Wilson, 1988). To account for these different mechanisms, we tested different variants of the RANS eddy-viscosity based models by performing a series of validation studies for flows over pockets of homogeneously or non-homogeneous distributed vegetation (trees), for a range of Reynolds numbers and the leaf area density parameters.
2. Mathematical formulation The simplest form of the vegetation model is based on the inclusion of the drag force term into the momentum equation. Since the fluctuating drag force also directly affects the turbulent kinetic energy, the more elaborate models include additional source/sink terms in both turbulence kinetic energy ðSD k Þ and its dissipation rate (SD ɛ ) equations as follows (Green, 1992): ∂U i ∂U ∂ ∂U i ∂U j 1 ∂P þ Uj i ¼ ν þ C D ajUjU i ð1Þ ui uj ffl{zfflfflfflfflfflffl} ∂t ∂xj ∂xj ∂xj ∂xi ρ ∂xi |fflfflfflfflffl FD i
∂k ∂k ∂ þ Uj ¼ ∂t ∂xj ∂xj
νþ
νt ∂k þP k ɛ sk ∂xj
þ C D aðβ p jUj3 βd jUjkÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
ð2Þ
ð3Þ
SD ɛ
qffiffiffiffiffiffi where CD, a and jUj ¼ U 2i are the drag coefficient, the leaf area density (LAD) and the speed, respectively. The standard set of model coefficients in equations of the turbulence kinetic energy and its dissipation rate is applied (C ɛ1 ; C ɛ2 ; sk ; sɛ ; C μ Þ ¼ ð1:44; 1:92; 1:; 1:3; 0:09Þ (Launder and Spalding, 1974). The production of turbulence kinetic energy Pk is calculated as P k ¼ 2νt jSj2 , where S is the rate of strain production defined as Sij ¼ 1=2ð∂U i =∂xj þ ∂U j =∂xi Þ.p Tffiffiffiis the Durbin time-scale limiter evaluated as T ¼ minðk=ɛ; 0:6= 6C μ jSjÞ (Durbin, 1991, 1995). Its main function is to suppress an excessive rate of the turbulent kinetic energy production in a stagnant region or in regions with a strong curvature of the flow. The turbulent stress and turbulent viscosity are calculated as follows: 2 ∂U i ∂U j ; νt ¼ C μ kT: ui uj ¼ kδij νt þ ð4Þ 3 ∂xj ∂xi For flows at realistic environmental scales, two important elements should be included in the treatment of the wall-boundary: the surface roughness and non-equilibrium wall-functions. The velocity wall-functions with specified surface roughness imply the calculation of the wall-friction as follows: U tP ER ynP ln en
1=2 τW ¼ ρC 1=4 μ kP κ
Model
βp
βd
C ɛ4
C ɛ5
Momentum only Svensson and Häggkvist (1990) Katul et al. (2004) Sanz (2003) Liu et al. (1998)
0.0 1.0 1.0 1.0 1.0
0.0 0.0 5.1 5.1 4.0
0.0 1.95 0.9 C ɛ4 ðαÞ 1.5
0.0 0.0 0.9 C ɛ5 ðαÞ 0.6
where the subscript ‘P’ indicates that the variable is evaluated in the cell center of the first numerical grid next to the wall, en is the non-dimensional surface roughness parameter, ER ¼30 is a constant, U tP is the wall tangential velocity and κ ¼ 0:41 is the von Kármán constant. Note that to normalize the distance to the wall 1=2 (ynP ), the definition U nτ ¼ C 1=4 μ kP is used to avoid singularity when τW tends to zero. To be able to deal at least partially with nonequilibrium wall-effects, a generalized form of the wall-friction is introduced as ! C yn τGW ¼ τW 1 tU Pn ð6Þ UP κUτ where the correction parameter CU is evaluated from the combination of the convective and pressure-gradient terms evaluated from the previous iteration (Kim and Choudhury, 1995). For example, for situations where it is assumed that the flow is mainly oriented in the direction of the streamwise U velocity component, the correction parameter can be evaluated as follows: CU ¼ U
SD k
∂ɛ ∂ɛ ∂ ν ∂ɛ þ Uj ¼ νþ t ∂t ∂xj ∂xj sɛ ∂xj C ɛ1 P k C ɛ2 ɛ ɛ þC D a C ɛ4 βp jUj3 C ɛ5 βd jUjɛ þ T k |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
Table 1 Vegetation model coefficients.
ð5Þ
∂U ∂U ∂U 1∂P þV þ W þ ∂x ∂y ∂z ρ ∂x
ð7Þ
The vegetation effects are modeled with additional terms in the momentum equation (FiD) and in equations for turbulence variables D (SD k , Sɛ ) (Green, 1992). The drag coefficient CD and the leaf area density (a) are the given quantities. The typical values of the drag coefficient for most vegetation are 0:1 r C D r 0:3 (Katul et al., 2004). The different types of vegetation can be easily mimicked by the nonuniform vertical distribution of the leaf area density parameter (a). While the vegetation acts as a sink term in the momentum equation, it can be the sink or the source term in equations for turbulence variables. As such, it displays a complex exchange of energy redistribution among the different velocity components inside and in the proximity of the area covered by vegetation. We compared different modeling approaches for the vegetation – starting from a simple inclusion of the drag force in the momentum equations only, to more elaborate models, which contain four additional terms in the equations for turbulence kinetic energy and its dissipation rate (Sanz, 2003; Katul et al., 2004; Liang et al., 2006; Mochida and Lun, 2008; Mochida et al., 2008). An intermediate modeling level consisting of only two additional terms is also validated (Svensson and Häggkvist, 1990). Finally, a model that included a spatial variation of the model 2=3 coefficients C ɛ4 ¼ C ɛ5 ¼ sk ½2=sɛ ðC 1=2 ðC ɛ2 C ɛ1 Þ, as proμ =6Þð2=αÞ posed in Sanz (2003), is analyzed. Here, the parameter α is evaluated from the characteristic mixing length and the Kolmogorov relation, 3=2 3=2 α ¼ lm C D a, lm ¼ C 3=4 =ɛ, i.e. we have α ¼ C D aC 3=4 =ɛ (Sanz, μ k μ k 2003). Now we will perform a comparative assessment of the above-listed models (Table 1) for vegetation (trees) canopy effects.
3. Numerical method The equations set (Eqs. (1)–(3)) is solved numerically with an in-house code based on a finite volume Navier–Stokes solver in structured non-orthogonal geometries (Kenjereš and Hanjalić, 2002, 2006). A collocated variable arrangement is employed for
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
all variables. The second-order central difference scheme (CDS) is used for all diffusive terms of conservative variables. The secondorder quadratic upwind scheme (QUICK) is applied for convective
7
7
6
6
5
5
4
4
45
terms. The coupling between velocity and pressure fields is provided through the SIMPLE algorithm. An important aspect of simulating flows in complex urban areas is how to efficiently generate a numerical mesh with required refinements in the proximity of the buildings or vegetation areas. We adopted the passive element approach presented in Kenjereš et al. (2001). In contrast to the conventional multi-block approach in structured geometries, where building-blocks occupying only the fluid domain are used, in the passive elements approach, the buildings (blocked regions) are excluded from the basic underlying mesh (mesh covering the empty-fluid domain) and treated as new solid boundaries by imposing proper wall-function boundary conditions. This method proved to be very efficient both in terms of the mesh generation and in terms of computational costs of performed simulations.
z / h [-]
z / h [-]
4. Results
3
We start our investigation by performing comparative assessments of different model extensions within the RANS two-equation k ɛ model on generic configurations. These include a turbulent boundary layer over and within a homogeneous (constant leaf area density) and a non-homogeneous (vertically variable leaf area density) vegetation region. Then, we focus on a complex real-scale urban area case with green islands (site of the TU Delft university campus).
3
2
2
1
1
0
0 0
2
2
4
6
8
0
10
0.5
1
4.1. Turbulent boundary layer over homogeneous vegetation
σu / U [-]
U [m / s]
Fig. 1. The simulated turbulent boundary layer development over a homogeneous vegetation (canopy), Reh ¼1.39 106, a¼ 2.1 m2/m3, CD ¼0.15. Geometrical specifications (top), the vertical profiles of the streamwise velocity (U) and of the variance of the horizontal velocity (sU =U) at the specified location (x ¼ 192.5 m) (bottom). The undisturbed approach flow (dashed lines), present simulations (solid lines), results of the CFD-RANS of Svensson and Häggkvist (1990) (symbols).
The development of the turbulent boundary layer over a homogeneously distributed vegetation is considered first (Fig. 1-top). This case was numerically solved and presented in Svensson and Häggkvist (1990) and is used in the present study for the very first verification of the numerical implementation of the additional drag
1
z / h [-]
0.8 0.6 0.4 0.2 0 0
10
20 2
30
40
3
a [m /m ] Fig. 2. The geometric settings for simulations over a non-homogeneously distributed vegetation (top) and measured characteristic leaf area density in the vertical direction, experiments by Green (1992) and CFD-RANS by Liang et al. (2006) (bottom).
46
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
force in the momentum equation and of two additional terms for turbulence variables. The region with vegetation is represented by a volume of 370 370 2.5 m3 in the streamwise, spanwise and vertical directions. Note that the vegetation expands over the entire spanwise (y) direction. The region between the inlet and the vegetation is 20 m. The segment behind the vegetation region has a length of 1480 m. The side and top boundaries are treated as a symmetry boundary condition. At the inlet, identical flow conditions are specified as in Svensson and Häggkvist (1990). The streamwise velocity profile is given as UðzÞ ¼ U 0 ðz=δ0 Þ1=7 , where U 0 ¼ 10 m=s and δ0 ¼ 110 m. The turbulence kinetic energy is specified as kðzÞ ¼ k0 ðLz zÞ=Lz , with k0 ¼1.5 m2/s2 and where Lz is the vertical extension of the computational domain. The dissipation profiles is defined as ɛðzÞ ¼ 0:41kðzÞ3=2 =z. The Reynolds number based on the vegetation height is Reh ¼ 1:39 106 . The numerical mesh consisted of 255 78 61 control volumes (CVs) and the vegetation region was represented by 65 78 13 CVs, in the x-, y- and z-directions. Numerical cell sizes range from 1.5 m 1.5 m 0.15 m to 10 m 6 m 4 m, with a grid expansion factor of 1.1 and refinement near the edges of the vegetation domain. The vegetation region is a rectangular block with a constant leaf area density a¼2.1 m2/m3 and a drag coefficient of CD ¼0.15. It can be seen that the horizontal velocity component and turbulence kinetic energy are significantly affected by the vegetation layer (Fig. 1, below, left). This reduction of the horizontal velocity is
present not only within the vegetation layer, but also extends vertically up to z/h¼5 vegetation heights. The non-dimensional variance of the horizontal velocity (sU =U) also shows a very strong increase in the vegetation region, but it rapidly decays toward values for the undisturbed flow in the upper part of the domain (Fig. 1-below, right). For both profiles, a good agreement between the present simulation and the results of Svensson's numerical study are obtained, confirming the proper implementation of the simplified model for the homogeneous vegetation. 4.2. Turbulent boundary layer flow over non-homogeneous vegetation The next case is a ‘vegetation island’ shown in Fig. 2-top, characterized by a non-uniform vertical distribution of the leaf area density parameter ‘a’ (Fig. 2-bottom). The total simulation domain extends for 10 m in the front and 35 m in the back of the vegetation region (with a total length of 5 m, a width of 5 m and a height of 1 m) in the streamwise (x) direction, and for 10 m from each side in the spanwise (y) direction, as shown in Fig. 2-top. Along the side- and top-boundaries, the symmetry boundary condition is applied, whereas the zero-gradient outlet boundary condition is applied at the end of the simulation domain. The approaching flow conditions are specified to be identical to those used by Liang et al. (2006). The inlet horizontal velocity
Fig. 3. Results of simulations with parameters specified in Fig. 2. Contours of the pressure (in Pa), streamwise and spanwise velocity (in [m/s]) in the characteristic horizontal (x y) plane z/h ¼ 0.25 (x- and y- are in [m]). The vegetation region (trees) is also indicated with a white box. For present plots Katul's model is used.
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
profiles is defined as UðzÞ ¼ 1:25 U 0 ðz=δ0 Þ1=6 , where U0 ¼6 m/s and δ0 ¼ 10 m. The turbulence kinetic energy profile is given as kðzÞ ¼ 3:33 U 2f ð1 z=δ1 Þ2 , where δ1 ¼ 5 m is the gradient height and Uf ¼0.15 m/s is the friction velocity. Finally, the dissipation rate is calculated as ɛðzÞ ¼ U 3f ð1 z=δ1 Þ2 =ðκ ðz þ z0 ÞÞ, where z0 ¼ 0.0025 m is the surface roughness. The characteristic Reynolds number based on the vegetation height is Reh ¼ U 0 h=ν ¼ 3:3 105 . The leaf area density parameter is now a function of the height (0:5 r aðzÞ r 36 m2 =m3 ), to mimic a realistic situation of the small scale trees, as was used in the wind tunnel experiments of Green (1992) and CFD RANS studies of Liang et al. (2006) (Fig. 2-bottom). It starts from small values in the lower part of the vegetation island, reaches its peak value at z/h¼0.6, and then decreases gradually toward the vegetation top. This is numerically mimicked by placing 10 homogeneous layers on top of each other, where each of the layers has the same drag coefficient of CD ¼ 0.1 (Fig. 2-bottom). The numerical mesh of 145 82 56 CVs is used from which 20 20 20 CVs are used within the vegetation region. Numerical cell sizes range from 0.25 m 0.25 m 0.05 m within the vegetation region to 0.4 m 0.4 m 0.2 m in the regions farther away. These regions are matched by a gradually increasing mesh size with a grid expansion factor of 1.05 in all directions. Different variants of the vegetation models (presented in Table 1) are tested. Changes in the mean flow caused by the presence of the vegetation island in the
47
horizontal plane (z/h¼0.25) are shown in Fig. 3. Three main physical phenomena can be observed: the pressure variations within the vegetation region (porous region), Fig. 3-top; generation of the wake region, Fig. 3-middle; and appearance of edge effects introducing the spanwise velocity component (blockage effect), Fig. 3-bottom. Contours of the turbulence kinetic energy (TKE) in the central vertical plane for different models are shown in Fig. 4. It can be seen that the Svensson model (Fig. 4-middle) strongly suppresses the TKE within the vegetation region while both the momentum-only (Fig. 4-top) and Katul's model (Fig. 4-bottom) show significant levels of TKE at the same location. To conclude which of the listed models gives the best physical representation of the vegetation, horizontal velocity and TKE profiles are compared with the available wind tunnel experimental data of Green (1992) along the three horizontal lines at heights z/h¼0.25, 0.75 and 1.25 (the first two are within the vegetation region) (Fig. 5). The horizontal velocity along the z/h¼0.25 (Fig. 5-top) line shows a slight decrease from the fully developed state up to the edge of the vegetation at x/h¼10 location. The air flow then accelerates, and after reaching its second peak at x/h¼12.5, strong suppression of the velocity takes place. Flow starts to recover after leaving the vegetation region after reaching the x/h¼20 position. All models reasonably predict appearances of the first peak for x/h¼0.25 and 1.25 locations, and of the flow recovery zone for all locations. The
Fig. 4. Contours of the turbulence kinetic energy TKE (in [m2/s2]) in the central vertical (x z) plane for a flow over vegetation island for three models for vegetation effects: momentum only (above), Svensson model (middle), Katul model (below) (x and z are in [m]).
48
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
5
z / h = 0.25 Experiment Momentum only Svensson model Katul model Sanz model
TREES
4
1.5
1
k [m2/s2]
3
U [m/s]
z / h = 0.25 Experiment Momentum only Svensson model Katul model Sanz model
TREES
2
0.5
1 0 -1
0 0
10
20
30
40
50
0
10
20
x / h [-]
5
z / h = 0.75 Experiment Momentum only Svensson model Katul model Sanz model
TREES
4
1.5
40
50
z / h = 0.75 Experiment Momentum only Svensson model Katul model Sanz model
TREES
1
k [m2/s2]
3
U [m/s]
30
x / h [-]
2
0.5
1 0 -1
0 0
10
20
30
40
50
0
10
20
x / h [-]
6
1.5
z / h = 1.25 Experiment Momentum only Svensson model Katul model Sanz model
TREES
5
30
40
50
x/h z / h = 1.25 Experiment Momentum only Svensson model Katul model Sanz model
TREES
1
k [m2/s2]
U [m/s]
4 3
0.5 2 1 0 0
0
10
20
30
40
50
x / h [-] Fig. 5. Results of simulations with parameters specified in Fig. 2. Horizontal profiles of the streamwise velocity component (in m/s) at three characteristic heights: z/h ¼0.25, 0.75, 1.25.
horizontal velocity exhibits different behaviors along the z/h¼0.75 line (Fig. 5-middle). Here, in contrast to the previous location, the second peak inside the forest is not present anymore and velocity is strongly suppressed along the entire vegetation region. Also, along this line, the flow recovery starts more rapidly, i.e. immediately after leaving the vegetation region. It is interesting to see that Svensson's model is outperformed along the z/h¼0.75 and 1.25 locations even by the model that includes only the drag force into the momentum
0
10
20
30
40
50
x / h [-] Fig. 6. Horizontal profiles of the turbulence kinetic energy (in m2/s2) along three characteristic locations: z/h ¼0.25, 0.75, 1.25.
equations. This indicates the importance of the inclusion of both mechanisms (production and dissipation) due to the presence of vegetation in equations with turbulence parameters. The profiles of TKE along identical lines show significant differences. As already observed in Fig. 4, Svensson's model suppresses the TKE too strongly within the vegetation region along all the three locations (Fig. 6). The momentum-only model entirely fails to capture the characteristic first peak of the TKE immediately within the vegetation region at z/h¼0.75 location (Fig. 6-middle). Here, both the Katul and the Sanz model are properly capturing characteristic peak
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
49
behavior. While Katul's model predicts the peak value closer to the experimental value at z/h¼0.75 location, it significantly overpredicts the peak value at z/h¼1.25 location (Fig. 6-bottom). It can be concluded that Katul's and Sanz's models provide the best overall agreement of the TKE profiles when compared to the experimental data. It should be noted that Katul's model proved to be more numerically stable than Sanz's model. To be able to get fully convergent results with Sanz's model, it was necessary to start simulation from the fully convergent solutions of Katul's model. 4.3. Case study: Delft University of Technology campus site Because of its overall best performances (physical accuracy and numerical robustness), the model of Katul is selected and applied to a representative urban area case (Delft University of Technology campus site) where changes in flow patterns and turbulence levels, due to the presence of vegetation (trees) (to mimic current situation in the summer period of 2012), are analyzed. The Reynolds number based on the characteristic height of the tallest building is Re¼3 107. The averaged values of the wind intensity and direction are collected at an undisturbed monitoring location in the proximity of the Rotterdam airport (located approximately 10 km from the campus site) and for a monitoring location 10 m above ground for a period of 50 years and used in generating
Fig. 8. Contours of the streamwise velocity (V) (top) and of the turbulence kinetic energy (TKE) in three characteristic vertical (y–z) planes, x ¼480, 660 and 840 m (in the coordinate system shown in Fig. 7).
Fig. 7. Top: the wind-rose map for the TU Delft campus site. Bottom: geometry of buildings and vegetation islands used for simulations of turbulent flow over Delft University of Technology site (a zoom-in). The entire geometry (including vegetation) covers a domain of 0 r x r 1200 m, 0 r y r 2100 m, 0 r z r 400 m and is represented by a numerical mesh consisting of 367 328 53 CVs. Buildings are represented by 86 passive (flow blockage) elements and in total 8 vegetation islands (locations marked by green elements) represented by porous volumes with a¼ 3 m2/m3 and CD ¼ 0.10. The arrow indicates approaching wind direction (the wind is blowing from the south-west).
inlet conditions based on simulated approaching fully developed neutral atmospheric boundary layer. The pre-simulations are performed to obtain a fully developed profile for all variables. This is done by simulating an empty domain with simple inflow and outflow boundary conditions, which were continuously updated (cyclic exchange of outflow profiles which are projected on the inlet plane after each cycle). The wind-rose map1 used for the present simulations is shown in Fig. 7-top. In total, 86 passive elements are used to represent all main features of characteristic buildings with an accuracy of 1 m (Fig. 7-bottom). The averaged values of the non-dimensional wall distances from all buildings are in the 30 rynþ r 100 range. A numerical mesh of 367 328 53 CVs is used to represent the entire simulation domain. The total simulation domain extends over 0 r x r 1200 m, 0 r y r 2100 m, 0 r z r 400 m. This is sufficient to allow for development of the wake region behind buildings along the y-coordinate direction. The pockets of trees along the main campus street are also taken into account for the present simulation. In total, eight vegetation islands are taken into account in simulations, as shown in Fig. 7-bottom.
1 The wind-rose map is created by using the raw wind-data from the Royal Dutch Meteorological Institute (KNMI) site; http://www.knmi.nl.
50
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
Fig. 9. Contours of the streamwise velocity component at z¼ 4 m (top) and z¼ 40 m (bottom).
The two-equation k ɛ model with the Durbin time-scale limiter proved to be numerically very efficient and robust in obtaining fully convergent solutions for a reference case, i.e. without vegetation. Then, starting from this solution, regions covered with trees are activated. The ‘green pockets’ shown in Fig. 7 mark the locations occupied with trees with a typical height of 15 m (where the vertical canopy is distributed between 5 and 15 m) and with a characteristic value of a ¼3 m2/m3 and CD ¼0.10. The flow details are revealed by plotting contours of the streamwise velocity component and of the turbulence kinetic energy in characteristic vertical and horizontal planes (Figs. 8 and 9). The flow reversals (recirculative regions) are present in the wakes of the buildings (Figs. 8-top and 9-top). Note that a significantly large wake is created for the tallest building, with a horizonal extension of up to six times the height of a building, as shown in Fig. 9-bottom. For the lower buildings, this horizontal extension of the recirculative region is considerably smaller (Fig. 9-top). The contours of the turbulence kinetic energy show significantly enhanced levels in the wake region at larger distances from the ground (Fig. 8-bottom). The next step is to activate vegetation effects. This is done for two characteristic scenarios of vegetation islands within the campus. The first corresponds to the current distribution of vegetation (situation of summer 2012), as shown in Figs. 7-bottom and 10-top. The current
Fig. 10. The locations of the vegetation islands within the TU Delft campus: the existing situation (summer 2012) (top) and a new proposal (bottom). In both cases, the vegetation islands are with the leaf area density of a ¼3 m2/m3. The goal is to reduce the intensive side-wind occurring in the region indicated by the red-dashed line in the upper figure. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)
situation generates strong side-winds around the ‘EWI’-building (the tallest building on the campus). This region is indicated by red colored ellipse in Fig. 10-top. The alternative new proposal is shown in Fig. 10-bottom. The goal is to perform a comparative assessment of three situations: campus without vegetation, present case with vegetation (summer 2012), and finally, the alternative (new) proposal. The contours of the modulus of velocity and of the kinetic energy of turbulence in the characteristic horizontal plane (z¼ 3.5 m, beneath the vegetation canopy, which starts at z¼5 m) for analyzed scenarios are shown in Fig. 11. Building cross-sections are shown in white. It can be seen that the presence of vegetation reduces significantly the intensity of wind at desired location (500 rx r600 m and 1400 r y r1550 m) (Fig. 11-left). Here, the vegetation cross-sections are marked by black colored rectangles. The characteristic side jet is practically not visible anymore for the case of new vegetation, indicating strong suppression of the wind intensity (from jV j ¼ 4 m=s to less than 1.5 m/s) (Fig. 11-left/bottom). This suppression is also clearly visible in the contours of the turbulence kinetic energy (Fig. 11-right). The reorganization of flow and turbulence is even more visible
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
51
Fig. 11. Contours of the modulus of velocity (left) and superimposed velocity vectors and contours of the turbulence kinetic energy (right) in the characteristic horizontal plane beneath the vegetation islands (z ¼ 3.5 m) for three considered cases: without vegetation (top), the current situation (middle) and a new alternative proposal (bottom). The approaching wind direction is from top to bottom (indicated by the arrow).
52
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
Fig. 12. Same as in Fig. 11, only now at z ¼10 m (the plane crosses vegetation islands that are marked by black rectangles).
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
53
Fig. 13. The contours of the turbulence kinetic energy in the characteristic horizontal plane (z ¼10 m) and indications of the characteristic lines (L1–L2–L3– L4) along which the profiles are extracted. The presented solution is for the new arrangement of the vegetation islands with a ¼3 m2/m3. Note that (L5) coincides with (L4) but is extracted at z ¼3.5 m.
Fig. 15. Same as in Fig. 14, only now for contours of the turbulence kinetic energy.
Fig. 14. Contours of the modulus of velocity in the characteristic vertical plane along the L1 line (y¼ 1530 m) (a zoom-in) for: the reference state (no vegetation) (top), current situation (middle) and proposed new alternative (bottom). The direction of the approaching wind is perpendicular to the plot.
in the horizontal plane that crosses the vegetation islands at z¼ 10 m shown in Fig. 12. For this location, even more dramatic changes are observed. Again, the EWI building side jet is completely gone for the new vegetation (Fig. 12-left/bottom). This is achieved by adding an additional vegetation island to affect the wind conditions (this new vegetation island is located between 425 rx r525 m and 1600r y r 1700 m). It can be clearly observed that the velocity vectors within the vegetation canopy are significantly reduced (from 5 m/s to 1 m/s) (Fig. 12-right). To provide more detailed insights into intensity of the wind and its energy of fluctuations for different vegetation distributions, the profiles along the selected locations, shown in Fig. 13, are analyzed next. The characteristic lines (L1–L2–L3) are perpendicular whereas (L4–L5) are aligned with the mean wind direction. The first four lines are extracted at middle of the vegetation canopy height (z ¼10 m), while the last one (L5) coincides with (L4) and is extracted at height of z ¼3.5 m (just beneath vegetation which starts at z ¼5 m). The contours of the velocity and of the turbulence kinetic energy in the vertical plane across the L1 line (a zoom-in) are shown in Figs. 14 and 15. The direction of the approaching wind is perpendicular to this plane. It can be seen that the flow is significantly suppressed within the vegetation regions (Fig. 14). The contours shown in Fig. 15 show that the TKE is locally enhanced along the upper edges of the vegetation islands, while a very low turbulence intensity regions are observed within the vegetation. The profiles of modulus of velocity along the L1–L3 lines are shown in Fig. 16. The typical cross-section with vegetation islands are marked by green colored rectangles. Along the L3 line, there is practically no difference between the a¼0 and 3
54
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
. |V| [m/s]
6 4
2
3
2
3
a=0 m /m 2 3 a=3 m /m (new)
-L3-
a=3 m /m (old)
2 0
300
400
500
600
700
x [m]
|V| [m/s]
6 4
2
3
2
3
a=0 m /m 2 3 a=3 m /m (new)
-L1-
a=3 m /m (old)
2 0
300
400
500
600
700
x [m]
|V| [m/s]
6 2
3
2
3
a=0 m /m 2 3 a=3 m /m (new)
4
(old) m2/m3 vegetation cases. The situation is significantly different when additional vegetation is placed here (425r x r 525 m), which dramatically reduced the intensity of the wind (less than 1 m/s within the vegetation canopy). Note that the zero values correspond to segments that cross the buildings. Similarly, along the L1 and L2 lines, the new vegetation arrangement produces a dramatic reduction in the wind intensity. Almost a complete suppression is now obtained along the L2 (Fig. 16-bottom). This is the result of the combined effects of the extended vegetation islands located here, but also from the initial suppression of the incoming wind by placing additional vegetation island at the L3. It stresses importance of the taking into account combined effects of the different vegetation islands within the complex urban areas (Fig. 17). Finally, the profiles along the L4 and L5 are shown in Fig. 18. Now, the approaching wind is blowing from right to left, as indicated by the upper arrow. At z¼10 an 3.5 m, peaks in jVj are observed at the start of the first vegetation cross-section. It is followed by rapid (L4) or gradual (L5) decrease within the vegetation region. This is the consequence of the absence of the vegetation resistance along the L4 (just beneath the vegetation canopy). Interesting behavior is
-L2-
a=3 m /m (old)
2 0
300
400
500
600
700
x [m] Fig. 16. The horizontal profiles of the modulus of velocity along the characteristic locations indicated in Fig. 13 (i.e. in the z ¼ 10 m plane). The locations of the crosssections with vegetation islands are marked by green rectangles. The direction of the approaching wind is perpendicular to the plot. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)
Fig. 17. Same as in Fig. 16 only now for the turbulence kinetic energy profiles.
Fig. 18. The characteristic profiles of the modulus of velocity and of turbulence kinetic energy along lines L4 (z ¼10 m) and L5 (z ¼3.5 m) (shown in Fig. 13) for different leaf area densities (a¼ 0, 1 and 3 m2/m3). Note that the direction of the approaching wind is from right to left. The location of the line L5 is just beneath the line L4 (Fig. 13).
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
observed to the TKE profile along the L5 line (Fig. 18-bottom). Here, the TKE is even enhanced beneath the vegetation at the first location (starting at y¼1700 m) and then is significantly suppressed (more than fourfold reduction) in the second segment (starting at y¼1475 m). The results obtained with a reduced leaf area density of a¼ 1 m2/m3 for the new proposal are also shown for locations L4 and L5. It can be seen that just some marginal differences are obtained – apart from the TKE profile at z¼ 3.5 m, along the first green segment (1600 r yr 1700 m). 5. Summary and conclusions In this study, we analyzed different variants of the RANS twoequation k ɛ turbulence model for proper treatment of the vegetation in urban areas. Starting with a homogeneously distributed vegetation in a developing turbulent boundary layer, the additional terms in momentum, and equations for turbulence variables, are analyzed. Then a turbulent flow over a non-homogeneous (in the vertical direction) forest island is studied. Based on comparative assessment of the different models with the available experimental data, it is shown that the overall best agreement was obtained with the model that includes four additional terms in the turbulence equations. Simplified versions of the model were not able to properly mimic the changes in the velocity fields, and consequently in the turbulence kinetic energy, for cases with vertically non-homogeneously distributed vegetation. Then, the selected version of the model was applied for a real case of flow in a complex urban area (TU Delft campus). Simulations for three different scenarios are conducted. These included: no-vegetation case, present (old) vegetation situation (summer 2012) and a new proposal. It is shown that the velocity components are suppressed inside the vegetation regions. In contrast to that, the turbulence kinetic energy is firstly enhanced at the edges of these regions, and then significantly suppressed. This gives some practical guidelines for possible suppressions of the very strong side-winds along street-canyons in the complex urban areas. It is concluded that the vegetation distribution mimicking the situation of TU Delft campus in summer of 2012 partially affected flow and turbulence distributions within the campus site. In contrast to this, the intensity of the notorious side-wind gust around the tallest building at the campus site was significantly reduced with a new distribution of the vegetation islands. Finally, it is concluded that the variant proposed here for the non-homogeneous vegetation model proved to be numerically robust and physically accurate, making it a potentially useful tool for analysis and assessments of sustainable solutions in urban areas. References Ayotte, K.W., Finnigan, J.J., Raupach, M.R., 1999. A second-order closure for neutrally stratified vegetative canopy flows. Boundary-Layer Meteorology 90 (2), 186–216. Balczo, M., Gromke, C., Ruck, B., 2009. Numerical modeling of flow and pollutant dispersion in street canyons with tree planting. Meteorologische Zeitschrift 18 (2), 197–206. Belcher, S.E., Harman, I.N., Finnigan, J.J., 2012. The wind in the willows: flows in forest canopies in complex terrain. Annual Review of Fluid Mechanics 44, 479–504.
55
Dupont, S., Brunet, Y., 2009. Coherent structures in canopy edge flow: a large-eddy simulation study. Journal of Fluid Mechanics 630, 93–128. Durbin, P.A., 1991. Near-wall turbulence closure modeling without damping functions. Theoretical and Computational Fluid Dynamics 3, 1–13. Durbin, P.A., 1995. Separated flow computations with the k ɛ v2 model. AIAA Journal 33, 659–664. Green, S.R., 1992. Modeling turbulent air flow in stand of widely spaced trees. Phoenics Journal 23, 294–312. Gromke, C., Buccolieri, R., Di Sabatino, S., Ruck, B., 2008. Dispersion study in a street canyon with tree planting by means of wind tunnel and numerical investigations – evaluation of CFD data with experimental data. Atmospheric Environment 42, 8640–8650. Hanjalić, K., Kenjereš, S., 2005. Dynamic simulation of pollutant dispersion over complex urban terrains: a tool for sustainable development, control and management. Energy 30 (8), 1481–1497. Hanjalić, K., Kenjereš, S., 2006. RANS-based VLES of thermal and magnetic convection at extreme conditions. Journal of Applied Mechanics – Transactions of the ASME 73 (3), 430–440. Hanjalić, K., Kenjereš, S., 2008. Some developments in turbulence modelling for wind and environmental engineering. Journal of Wind Engineering & Industrial Aerodynamics 96 (10–11), 1537–1570. Katul, G.G., Mahrt, L., Poggi, D., Sanz, C., 2004. One- and two equation models for canopy turbulence. Boundary-Layer Meteorology 113, 81–109. Kenjereš, S., Hanjalić, K., Gunarjo, S.B., 2001. A T-RANS/LES approach to indoor climate simulations. In: Proceedings of ASME Fluids Engineering Division, Paper No. FEDSM20002-31400, Montreal, Canada. Kenjereš, S., Hanjalić, K., 2002. Combined effects of terrain orography and thermal stratification on pollutant dispersion in a town valley: a T-RANS simulation. Journal of Turbulence 3 (026), 1–25. Kenjereš, S., Gunarjo, S.B., Hanjalić, K., 2005. Contribution to elliptic relaxation modelling of turbulent natural and mixed convection. International Journal of Heat and Fluid Flow 26 (4), 569–586. Kenjereš, S., Hanjalić, K., 2006. LES, T-RANS and hybrid simulations of thermal convection at high Ra numbers. International Journal of Heat and Fluid Flow 27 (5), 800–810. Kenjereš, S., Hanjalić, K., 2009. Tackling complex turbulent flows with transient RANS. Fluid Dynamics Research 41 (1), 1–32. (article number: 012201). Kim, S.E., Choudhury, D., 1995. A near-wall treatment using wall functions sensitized to pressure gradient. In: Separated and Complex Flows, ASME FED vol. 217, ASME. Kuhn, S., Rudolf von Rohr, P., Kenjereš, S., 2010. Large eddy simulations of wall heat transfer and coherent structures in mixed convection over a wavy wall. International Journal of Thermal Sciences 49 (7), 1209–1226. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering 3, 269–289. Liang, L., Xiaofeng, L., Borong, L., Yingxin, Z., 2006. Improved k ɛ two-equation turbulence model for canopy flow. Atmospheric Environment 40, 762–770. Liu, J., Chen, J.M., Black, T.A., Novak, M.D., 1998. E–ɛ modelling of turbulent air flow down wind of a model forest edge. Boundary-Layer Meteorology 72, 21–44. Mochida, A., Lun, I.Y.F., 2008. Prediction of wind environment and thermal comfort at pedestrian level in urban area. Journal of Wind Engineering and Industrial Aerodynamics 96, 1498–1527. Mochida, A., Tabata, Y., Iwata, T., Yoshino, H., 2008. Examining tree canopy models for CFD prediction of wind environment at pedestrian level. Journal of Wind Engineering and Industrial Aerodynamics 96, 1667–1677. Raupach, M.R., Thom, S.M., 1981. Turbulence in and above plant canopies. Annual Review of Fluid Mechanics 13, 97–129. Raupach, M.R., Shaw, R.H., 1982. Averaging procedures for flow within vegetation canopies. Boundary-Layer Meteorology 22, 79–90. Sanz, C., 2003. A note on k ɛ modelling of vegetation canopy air-flows. BoundaryLayer Meteorology 108, 191–197. Svensson, U., Häggkvist, 1990. A two-equation turbulence model for canopy flows. Journal of Wind Engineering and Industrial Aerodynamics 35, 201–211. Wagner, C., Kenjereš, S., Rudolf von Rohr, P., 2011. Dynamic large eddy simulations of momentum and wall heat transfer in forced convection over wavy surfaces. Journal of Turbulence 12 (7), 1–27. Wilson, J.D., 1988. A second-order closure model for flow through vegetation. Boundary-Layer Meteorology 42, 371–392. Wilson, J.D., Finnigan, J.J., Raupach, M.R., 1998. A first-order closure for disturbed plant-canopy flows, and its application to winds in a canopy on a ridge. Quarterly Journal of the Royal Meteorological Society 124, 705–732. Yi, C., 2008. Momentum transfer within canopies. Journal of Applied Meteorology and Climatology 47, 262–275.
Contents lists available at ScienceDirect
Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia
Modelling and simulations of turbulent flows in urban areas with vegetation Saša Kenjereš a,n, Benjamin ter Kuile b a Transport Phenomena Section, Department of Chemical Engineering, Faculty of Applied Sciences and J.M. Burgers Centre for Fluid Dynamics, Delft University of Technology, Julianalaan 136, 2628 BL Delft, The Netherlands b UFlows, Ezelsveldlaan 153, 2611 DM Delft, The Netherlands
art ic l e i nf o
a b s t r a c t
Article history: Received 19 June 2012 Received in revised form 3 September 2013 Accepted 10 September 2013 Available online 30 October 2013
A comparative assessment of the Reynolds-Averaged Navier–Stokes (RANS) two-equation k ɛ model, with different proposals for a proper treatment of the vegetation (trees) effects in urban areas, is reported here. The most common variants proposed in the literature based on a two-equation k ɛ eddy viscosity model are revisited. An extensive model validation is performed on generic cases with homogeneous and non-homogeneous vegetation and compared with available experimental data. It is shown that the minimal level of the proper modeling of the non-homogeneous (i.e. with the vertically variable leaf area density) vegetation requires a model with additional four terms in the equations for the turbulence kinetic energy and its dissipation rate. The proposed version of the model is then applied to a computational study of flow over the Delft University of Technology site with different scenarios of the vegetation islands. By performing series of simulations, some practical guidelines are suggested to reduce intensity of the side-wind gust within the campus. In conclusion, it is demonstrated that the adopted turbulence model for non-homogeneous vegetation, in combination with the passive element approach for buildings, proved to be a physically accurate and numerically robust method. The model is recommended for future use in simulating turbulent flows in complex urban areas with vegetation. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Environmental flows Urban canopy Trees Vegetation Turbulence modeling RANS-CFD simulations
1. Introduction Severe air pollution in urban areas is usually caused by the combined effects of intensive traffic/industrial emissions and critical meteorological conditions (windless days, strong stratification, temperature inversion capping, etc.). Since the origin of all environmental pollution is in complex spatial and temporal flow, turbulence and dispersion interactions, computer simulations can be powerful tools for estimating the impacts of different parameters on pollution. Air flow and pollutant dispersion at environmental scales are characterized by high values of Reynolds numbers (indicating highly turbulent flows) and complex geometries (topography of terrain simulated, arrays of buildings and streets). Obviously, direct numerical resolving of all flow and turbulence scales is beyond the reach of present computational power, so mathematical models for scales not resolved by underlying numerical grid (subgrid/subscale contributions) must be introduced. In our previous work, we addressed modeling and simulations of some important aspects of the environmental flows, including the turbulent dispersion (passive scalars) in street canyons or over
n
Corresponding author. Tel.: þ 31 152783649. E-mail addresses: [email protected], [email protected] (S. Kenjereš).
0167-6105/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jweia.2013.09.007
wavy-walls in the laboratory conditions and some real-scale applications dealing with the time-dependent (diurnal cycles) simulations of flow, turbulence and dispersion over a complex terrain orography, e.g. Kenjereš and Hanjalić (2002, 2009), Kenjereš et al. (2005), Hanjalić and Kenjereš (2005, 2006, 2008), Kuhn et al. (2010), and Wagner et al. (2011). This is a part of the development of an integral in-house computational code, which is able to simulate the laboratory scale and real-scale environmental flows with advanced numerical methods and turbulence models. The latter also include a recently developed seamless hybrid RANS/LES approach (Kenjereš and Hanjalić, 2006). In this paper, special focus will be placed on modeling and validation of effects of vegetation (trees) on flow and turbulence in urban areas, in the framework of the two-equation eddy-viscosity models. Such models can be used for simulations of urban areas covered by significant portions of green areas (parks) or where trees are distributed along streets. These models can also be used to assess effects of the urban vegetation on dispersion of the traffic exhaust gases in the street canyons, which is an important prerequisite to achieve a sustainable development of urban areas (Gromke et al., 2008; Balczo et al., 2009). The proper modeling of vegetation effects on flow, turbulence and turbulent dispersion represents a challenging task. The flow reorganization inside and above an area covered by vegetation involves complex transfer of the turbulence
44
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
kinetic energy (Raupach and Thom, 1981; Raupach and Shaw, 1982; Wilson, 1988; Wilson et al., 1998; Ayotte et al., 1999; Dupont and Brunet, 2009; Yi, 2008; Belcher et al., 2012). Due to a blocking effect (the flow is slowing down inside the vegetation area due to the drag force) strong velocity gradients can be created producing a significant source of the turbulence kinetic energy. This blocking can also lead to the formation of long recirculative regions (wakes) similar to the flows over bluff bodies. On the other hand, inside the vegetation area, leaves act as turbulence suppressors – enhancing its dissipation rate. This is done by generating small scale eddies, without undergoing the entire process of cascade interactions (a so-called shortcut in the turbulence cascade process, Wilson, 1988). To account for these different mechanisms, we tested different variants of the RANS eddy-viscosity based models by performing a series of validation studies for flows over pockets of homogeneously or non-homogeneous distributed vegetation (trees), for a range of Reynolds numbers and the leaf area density parameters.
2. Mathematical formulation The simplest form of the vegetation model is based on the inclusion of the drag force term into the momentum equation. Since the fluctuating drag force also directly affects the turbulent kinetic energy, the more elaborate models include additional source/sink terms in both turbulence kinetic energy ðSD k Þ and its dissipation rate (SD ɛ ) equations as follows (Green, 1992): ∂U i ∂U ∂ ∂U i ∂U j 1 ∂P þ Uj i ¼ ν þ C D ajUjU i ð1Þ ui uj ffl{zfflfflfflfflfflffl} ∂t ∂xj ∂xj ∂xj ∂xi ρ ∂xi |fflfflfflfflffl FD i
∂k ∂k ∂ þ Uj ¼ ∂t ∂xj ∂xj
νþ
νt ∂k þP k ɛ sk ∂xj
þ C D aðβ p jUj3 βd jUjkÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
ð2Þ
ð3Þ
SD ɛ
qffiffiffiffiffiffi where CD, a and jUj ¼ U 2i are the drag coefficient, the leaf area density (LAD) and the speed, respectively. The standard set of model coefficients in equations of the turbulence kinetic energy and its dissipation rate is applied (C ɛ1 ; C ɛ2 ; sk ; sɛ ; C μ Þ ¼ ð1:44; 1:92; 1:; 1:3; 0:09Þ (Launder and Spalding, 1974). The production of turbulence kinetic energy Pk is calculated as P k ¼ 2νt jSj2 , where S is the rate of strain production defined as Sij ¼ 1=2ð∂U i =∂xj þ ∂U j =∂xi Þ.p Tffiffiffiis the Durbin time-scale limiter evaluated as T ¼ minðk=ɛ; 0:6= 6C μ jSjÞ (Durbin, 1991, 1995). Its main function is to suppress an excessive rate of the turbulent kinetic energy production in a stagnant region or in regions with a strong curvature of the flow. The turbulent stress and turbulent viscosity are calculated as follows: 2 ∂U i ∂U j ; νt ¼ C μ kT: ui uj ¼ kδij νt þ ð4Þ 3 ∂xj ∂xi For flows at realistic environmental scales, two important elements should be included in the treatment of the wall-boundary: the surface roughness and non-equilibrium wall-functions. The velocity wall-functions with specified surface roughness imply the calculation of the wall-friction as follows: U tP ER ynP ln en
1=2 τW ¼ ρC 1=4 μ kP κ
Model
βp
βd
C ɛ4
C ɛ5
Momentum only Svensson and Häggkvist (1990) Katul et al. (2004) Sanz (2003) Liu et al. (1998)
0.0 1.0 1.0 1.0 1.0
0.0 0.0 5.1 5.1 4.0
0.0 1.95 0.9 C ɛ4 ðαÞ 1.5
0.0 0.0 0.9 C ɛ5 ðαÞ 0.6
where the subscript ‘P’ indicates that the variable is evaluated in the cell center of the first numerical grid next to the wall, en is the non-dimensional surface roughness parameter, ER ¼30 is a constant, U tP is the wall tangential velocity and κ ¼ 0:41 is the von Kármán constant. Note that to normalize the distance to the wall 1=2 (ynP ), the definition U nτ ¼ C 1=4 μ kP is used to avoid singularity when τW tends to zero. To be able to deal at least partially with nonequilibrium wall-effects, a generalized form of the wall-friction is introduced as ! C yn τGW ¼ τW 1 tU Pn ð6Þ UP κUτ where the correction parameter CU is evaluated from the combination of the convective and pressure-gradient terms evaluated from the previous iteration (Kim and Choudhury, 1995). For example, for situations where it is assumed that the flow is mainly oriented in the direction of the streamwise U velocity component, the correction parameter can be evaluated as follows: CU ¼ U
SD k
∂ɛ ∂ɛ ∂ ν ∂ɛ þ Uj ¼ νþ t ∂t ∂xj ∂xj sɛ ∂xj C ɛ1 P k C ɛ2 ɛ ɛ þC D a C ɛ4 βp jUj3 C ɛ5 βd jUjɛ þ T k |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}
Table 1 Vegetation model coefficients.
ð5Þ
∂U ∂U ∂U 1∂P þV þ W þ ∂x ∂y ∂z ρ ∂x
ð7Þ
The vegetation effects are modeled with additional terms in the momentum equation (FiD) and in equations for turbulence variables D (SD k , Sɛ ) (Green, 1992). The drag coefficient CD and the leaf area density (a) are the given quantities. The typical values of the drag coefficient for most vegetation are 0:1 r C D r 0:3 (Katul et al., 2004). The different types of vegetation can be easily mimicked by the nonuniform vertical distribution of the leaf area density parameter (a). While the vegetation acts as a sink term in the momentum equation, it can be the sink or the source term in equations for turbulence variables. As such, it displays a complex exchange of energy redistribution among the different velocity components inside and in the proximity of the area covered by vegetation. We compared different modeling approaches for the vegetation – starting from a simple inclusion of the drag force in the momentum equations only, to more elaborate models, which contain four additional terms in the equations for turbulence kinetic energy and its dissipation rate (Sanz, 2003; Katul et al., 2004; Liang et al., 2006; Mochida and Lun, 2008; Mochida et al., 2008). An intermediate modeling level consisting of only two additional terms is also validated (Svensson and Häggkvist, 1990). Finally, a model that included a spatial variation of the model 2=3 coefficients C ɛ4 ¼ C ɛ5 ¼ sk ½2=sɛ ðC 1=2 ðC ɛ2 C ɛ1 Þ, as proμ =6Þð2=αÞ posed in Sanz (2003), is analyzed. Here, the parameter α is evaluated from the characteristic mixing length and the Kolmogorov relation, 3=2 3=2 α ¼ lm C D a, lm ¼ C 3=4 =ɛ, i.e. we have α ¼ C D aC 3=4 =ɛ (Sanz, μ k μ k 2003). Now we will perform a comparative assessment of the above-listed models (Table 1) for vegetation (trees) canopy effects.
3. Numerical method The equations set (Eqs. (1)–(3)) is solved numerically with an in-house code based on a finite volume Navier–Stokes solver in structured non-orthogonal geometries (Kenjereš and Hanjalić, 2002, 2006). A collocated variable arrangement is employed for
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
all variables. The second-order central difference scheme (CDS) is used for all diffusive terms of conservative variables. The secondorder quadratic upwind scheme (QUICK) is applied for convective
7
7
6
6
5
5
4
4
45
terms. The coupling between velocity and pressure fields is provided through the SIMPLE algorithm. An important aspect of simulating flows in complex urban areas is how to efficiently generate a numerical mesh with required refinements in the proximity of the buildings or vegetation areas. We adopted the passive element approach presented in Kenjereš et al. (2001). In contrast to the conventional multi-block approach in structured geometries, where building-blocks occupying only the fluid domain are used, in the passive elements approach, the buildings (blocked regions) are excluded from the basic underlying mesh (mesh covering the empty-fluid domain) and treated as new solid boundaries by imposing proper wall-function boundary conditions. This method proved to be very efficient both in terms of the mesh generation and in terms of computational costs of performed simulations.
z / h [-]
z / h [-]
4. Results
3
We start our investigation by performing comparative assessments of different model extensions within the RANS two-equation k ɛ model on generic configurations. These include a turbulent boundary layer over and within a homogeneous (constant leaf area density) and a non-homogeneous (vertically variable leaf area density) vegetation region. Then, we focus on a complex real-scale urban area case with green islands (site of the TU Delft university campus).
3
2
2
1
1
0
0 0
2
2
4
6
8
0
10
0.5
1
4.1. Turbulent boundary layer over homogeneous vegetation
σu / U [-]
U [m / s]
Fig. 1. The simulated turbulent boundary layer development over a homogeneous vegetation (canopy), Reh ¼1.39 106, a¼ 2.1 m2/m3, CD ¼0.15. Geometrical specifications (top), the vertical profiles of the streamwise velocity (U) and of the variance of the horizontal velocity (sU =U) at the specified location (x ¼ 192.5 m) (bottom). The undisturbed approach flow (dashed lines), present simulations (solid lines), results of the CFD-RANS of Svensson and Häggkvist (1990) (symbols).
The development of the turbulent boundary layer over a homogeneously distributed vegetation is considered first (Fig. 1-top). This case was numerically solved and presented in Svensson and Häggkvist (1990) and is used in the present study for the very first verification of the numerical implementation of the additional drag
1
z / h [-]
0.8 0.6 0.4 0.2 0 0
10
20 2
30
40
3
a [m /m ] Fig. 2. The geometric settings for simulations over a non-homogeneously distributed vegetation (top) and measured characteristic leaf area density in the vertical direction, experiments by Green (1992) and CFD-RANS by Liang et al. (2006) (bottom).
46
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
force in the momentum equation and of two additional terms for turbulence variables. The region with vegetation is represented by a volume of 370 370 2.5 m3 in the streamwise, spanwise and vertical directions. Note that the vegetation expands over the entire spanwise (y) direction. The region between the inlet and the vegetation is 20 m. The segment behind the vegetation region has a length of 1480 m. The side and top boundaries are treated as a symmetry boundary condition. At the inlet, identical flow conditions are specified as in Svensson and Häggkvist (1990). The streamwise velocity profile is given as UðzÞ ¼ U 0 ðz=δ0 Þ1=7 , where U 0 ¼ 10 m=s and δ0 ¼ 110 m. The turbulence kinetic energy is specified as kðzÞ ¼ k0 ðLz zÞ=Lz , with k0 ¼1.5 m2/s2 and where Lz is the vertical extension of the computational domain. The dissipation profiles is defined as ɛðzÞ ¼ 0:41kðzÞ3=2 =z. The Reynolds number based on the vegetation height is Reh ¼ 1:39 106 . The numerical mesh consisted of 255 78 61 control volumes (CVs) and the vegetation region was represented by 65 78 13 CVs, in the x-, y- and z-directions. Numerical cell sizes range from 1.5 m 1.5 m 0.15 m to 10 m 6 m 4 m, with a grid expansion factor of 1.1 and refinement near the edges of the vegetation domain. The vegetation region is a rectangular block with a constant leaf area density a¼2.1 m2/m3 and a drag coefficient of CD ¼0.15. It can be seen that the horizontal velocity component and turbulence kinetic energy are significantly affected by the vegetation layer (Fig. 1, below, left). This reduction of the horizontal velocity is
present not only within the vegetation layer, but also extends vertically up to z/h¼5 vegetation heights. The non-dimensional variance of the horizontal velocity (sU =U) also shows a very strong increase in the vegetation region, but it rapidly decays toward values for the undisturbed flow in the upper part of the domain (Fig. 1-below, right). For both profiles, a good agreement between the present simulation and the results of Svensson's numerical study are obtained, confirming the proper implementation of the simplified model for the homogeneous vegetation. 4.2. Turbulent boundary layer flow over non-homogeneous vegetation The next case is a ‘vegetation island’ shown in Fig. 2-top, characterized by a non-uniform vertical distribution of the leaf area density parameter ‘a’ (Fig. 2-bottom). The total simulation domain extends for 10 m in the front and 35 m in the back of the vegetation region (with a total length of 5 m, a width of 5 m and a height of 1 m) in the streamwise (x) direction, and for 10 m from each side in the spanwise (y) direction, as shown in Fig. 2-top. Along the side- and top-boundaries, the symmetry boundary condition is applied, whereas the zero-gradient outlet boundary condition is applied at the end of the simulation domain. The approaching flow conditions are specified to be identical to those used by Liang et al. (2006). The inlet horizontal velocity
Fig. 3. Results of simulations with parameters specified in Fig. 2. Contours of the pressure (in Pa), streamwise and spanwise velocity (in [m/s]) in the characteristic horizontal (x y) plane z/h ¼ 0.25 (x- and y- are in [m]). The vegetation region (trees) is also indicated with a white box. For present plots Katul's model is used.
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
profiles is defined as UðzÞ ¼ 1:25 U 0 ðz=δ0 Þ1=6 , where U0 ¼6 m/s and δ0 ¼ 10 m. The turbulence kinetic energy profile is given as kðzÞ ¼ 3:33 U 2f ð1 z=δ1 Þ2 , where δ1 ¼ 5 m is the gradient height and Uf ¼0.15 m/s is the friction velocity. Finally, the dissipation rate is calculated as ɛðzÞ ¼ U 3f ð1 z=δ1 Þ2 =ðκ ðz þ z0 ÞÞ, where z0 ¼ 0.0025 m is the surface roughness. The characteristic Reynolds number based on the vegetation height is Reh ¼ U 0 h=ν ¼ 3:3 105 . The leaf area density parameter is now a function of the height (0:5 r aðzÞ r 36 m2 =m3 ), to mimic a realistic situation of the small scale trees, as was used in the wind tunnel experiments of Green (1992) and CFD RANS studies of Liang et al. (2006) (Fig. 2-bottom). It starts from small values in the lower part of the vegetation island, reaches its peak value at z/h¼0.6, and then decreases gradually toward the vegetation top. This is numerically mimicked by placing 10 homogeneous layers on top of each other, where each of the layers has the same drag coefficient of CD ¼ 0.1 (Fig. 2-bottom). The numerical mesh of 145 82 56 CVs is used from which 20 20 20 CVs are used within the vegetation region. Numerical cell sizes range from 0.25 m 0.25 m 0.05 m within the vegetation region to 0.4 m 0.4 m 0.2 m in the regions farther away. These regions are matched by a gradually increasing mesh size with a grid expansion factor of 1.05 in all directions. Different variants of the vegetation models (presented in Table 1) are tested. Changes in the mean flow caused by the presence of the vegetation island in the
47
horizontal plane (z/h¼0.25) are shown in Fig. 3. Three main physical phenomena can be observed: the pressure variations within the vegetation region (porous region), Fig. 3-top; generation of the wake region, Fig. 3-middle; and appearance of edge effects introducing the spanwise velocity component (blockage effect), Fig. 3-bottom. Contours of the turbulence kinetic energy (TKE) in the central vertical plane for different models are shown in Fig. 4. It can be seen that the Svensson model (Fig. 4-middle) strongly suppresses the TKE within the vegetation region while both the momentum-only (Fig. 4-top) and Katul's model (Fig. 4-bottom) show significant levels of TKE at the same location. To conclude which of the listed models gives the best physical representation of the vegetation, horizontal velocity and TKE profiles are compared with the available wind tunnel experimental data of Green (1992) along the three horizontal lines at heights z/h¼0.25, 0.75 and 1.25 (the first two are within the vegetation region) (Fig. 5). The horizontal velocity along the z/h¼0.25 (Fig. 5-top) line shows a slight decrease from the fully developed state up to the edge of the vegetation at x/h¼10 location. The air flow then accelerates, and after reaching its second peak at x/h¼12.5, strong suppression of the velocity takes place. Flow starts to recover after leaving the vegetation region after reaching the x/h¼20 position. All models reasonably predict appearances of the first peak for x/h¼0.25 and 1.25 locations, and of the flow recovery zone for all locations. The
Fig. 4. Contours of the turbulence kinetic energy TKE (in [m2/s2]) in the central vertical (x z) plane for a flow over vegetation island for three models for vegetation effects: momentum only (above), Svensson model (middle), Katul model (below) (x and z are in [m]).
48
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
5
z / h = 0.25 Experiment Momentum only Svensson model Katul model Sanz model
TREES
4
1.5
1
k [m2/s2]
3
U [m/s]
z / h = 0.25 Experiment Momentum only Svensson model Katul model Sanz model
TREES
2
0.5
1 0 -1
0 0
10
20
30
40
50
0
10
20
x / h [-]
5
z / h = 0.75 Experiment Momentum only Svensson model Katul model Sanz model
TREES
4
1.5
40
50
z / h = 0.75 Experiment Momentum only Svensson model Katul model Sanz model
TREES
1
k [m2/s2]
3
U [m/s]
30
x / h [-]
2
0.5
1 0 -1
0 0
10
20
30
40
50
0
10
20
x / h [-]
6
1.5
z / h = 1.25 Experiment Momentum only Svensson model Katul model Sanz model
TREES
5
30
40
50
x/h z / h = 1.25 Experiment Momentum only Svensson model Katul model Sanz model
TREES
1
k [m2/s2]
U [m/s]
4 3
0.5 2 1 0 0
0
10
20
30
40
50
x / h [-] Fig. 5. Results of simulations with parameters specified in Fig. 2. Horizontal profiles of the streamwise velocity component (in m/s) at three characteristic heights: z/h ¼0.25, 0.75, 1.25.
horizontal velocity exhibits different behaviors along the z/h¼0.75 line (Fig. 5-middle). Here, in contrast to the previous location, the second peak inside the forest is not present anymore and velocity is strongly suppressed along the entire vegetation region. Also, along this line, the flow recovery starts more rapidly, i.e. immediately after leaving the vegetation region. It is interesting to see that Svensson's model is outperformed along the z/h¼0.75 and 1.25 locations even by the model that includes only the drag force into the momentum
0
10
20
30
40
50
x / h [-] Fig. 6. Horizontal profiles of the turbulence kinetic energy (in m2/s2) along three characteristic locations: z/h ¼0.25, 0.75, 1.25.
equations. This indicates the importance of the inclusion of both mechanisms (production and dissipation) due to the presence of vegetation in equations with turbulence parameters. The profiles of TKE along identical lines show significant differences. As already observed in Fig. 4, Svensson's model suppresses the TKE too strongly within the vegetation region along all the three locations (Fig. 6). The momentum-only model entirely fails to capture the characteristic first peak of the TKE immediately within the vegetation region at z/h¼0.75 location (Fig. 6-middle). Here, both the Katul and the Sanz model are properly capturing characteristic peak
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
49
behavior. While Katul's model predicts the peak value closer to the experimental value at z/h¼0.75 location, it significantly overpredicts the peak value at z/h¼1.25 location (Fig. 6-bottom). It can be concluded that Katul's and Sanz's models provide the best overall agreement of the TKE profiles when compared to the experimental data. It should be noted that Katul's model proved to be more numerically stable than Sanz's model. To be able to get fully convergent results with Sanz's model, it was necessary to start simulation from the fully convergent solutions of Katul's model. 4.3. Case study: Delft University of Technology campus site Because of its overall best performances (physical accuracy and numerical robustness), the model of Katul is selected and applied to a representative urban area case (Delft University of Technology campus site) where changes in flow patterns and turbulence levels, due to the presence of vegetation (trees) (to mimic current situation in the summer period of 2012), are analyzed. The Reynolds number based on the characteristic height of the tallest building is Re¼3 107. The averaged values of the wind intensity and direction are collected at an undisturbed monitoring location in the proximity of the Rotterdam airport (located approximately 10 km from the campus site) and for a monitoring location 10 m above ground for a period of 50 years and used in generating
Fig. 8. Contours of the streamwise velocity (V) (top) and of the turbulence kinetic energy (TKE) in three characteristic vertical (y–z) planes, x ¼480, 660 and 840 m (in the coordinate system shown in Fig. 7).
Fig. 7. Top: the wind-rose map for the TU Delft campus site. Bottom: geometry of buildings and vegetation islands used for simulations of turbulent flow over Delft University of Technology site (a zoom-in). The entire geometry (including vegetation) covers a domain of 0 r x r 1200 m, 0 r y r 2100 m, 0 r z r 400 m and is represented by a numerical mesh consisting of 367 328 53 CVs. Buildings are represented by 86 passive (flow blockage) elements and in total 8 vegetation islands (locations marked by green elements) represented by porous volumes with a¼ 3 m2/m3 and CD ¼ 0.10. The arrow indicates approaching wind direction (the wind is blowing from the south-west).
inlet conditions based on simulated approaching fully developed neutral atmospheric boundary layer. The pre-simulations are performed to obtain a fully developed profile for all variables. This is done by simulating an empty domain with simple inflow and outflow boundary conditions, which were continuously updated (cyclic exchange of outflow profiles which are projected on the inlet plane after each cycle). The wind-rose map1 used for the present simulations is shown in Fig. 7-top. In total, 86 passive elements are used to represent all main features of characteristic buildings with an accuracy of 1 m (Fig. 7-bottom). The averaged values of the non-dimensional wall distances from all buildings are in the 30 rynþ r 100 range. A numerical mesh of 367 328 53 CVs is used to represent the entire simulation domain. The total simulation domain extends over 0 r x r 1200 m, 0 r y r 2100 m, 0 r z r 400 m. This is sufficient to allow for development of the wake region behind buildings along the y-coordinate direction. The pockets of trees along the main campus street are also taken into account for the present simulation. In total, eight vegetation islands are taken into account in simulations, as shown in Fig. 7-bottom.
1 The wind-rose map is created by using the raw wind-data from the Royal Dutch Meteorological Institute (KNMI) site; http://www.knmi.nl.
50
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
Fig. 9. Contours of the streamwise velocity component at z¼ 4 m (top) and z¼ 40 m (bottom).
The two-equation k ɛ model with the Durbin time-scale limiter proved to be numerically very efficient and robust in obtaining fully convergent solutions for a reference case, i.e. without vegetation. Then, starting from this solution, regions covered with trees are activated. The ‘green pockets’ shown in Fig. 7 mark the locations occupied with trees with a typical height of 15 m (where the vertical canopy is distributed between 5 and 15 m) and with a characteristic value of a ¼3 m2/m3 and CD ¼0.10. The flow details are revealed by plotting contours of the streamwise velocity component and of the turbulence kinetic energy in characteristic vertical and horizontal planes (Figs. 8 and 9). The flow reversals (recirculative regions) are present in the wakes of the buildings (Figs. 8-top and 9-top). Note that a significantly large wake is created for the tallest building, with a horizonal extension of up to six times the height of a building, as shown in Fig. 9-bottom. For the lower buildings, this horizontal extension of the recirculative region is considerably smaller (Fig. 9-top). The contours of the turbulence kinetic energy show significantly enhanced levels in the wake region at larger distances from the ground (Fig. 8-bottom). The next step is to activate vegetation effects. This is done for two characteristic scenarios of vegetation islands within the campus. The first corresponds to the current distribution of vegetation (situation of summer 2012), as shown in Figs. 7-bottom and 10-top. The current
Fig. 10. The locations of the vegetation islands within the TU Delft campus: the existing situation (summer 2012) (top) and a new proposal (bottom). In both cases, the vegetation islands are with the leaf area density of a ¼3 m2/m3. The goal is to reduce the intensive side-wind occurring in the region indicated by the red-dashed line in the upper figure. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)
situation generates strong side-winds around the ‘EWI’-building (the tallest building on the campus). This region is indicated by red colored ellipse in Fig. 10-top. The alternative new proposal is shown in Fig. 10-bottom. The goal is to perform a comparative assessment of three situations: campus without vegetation, present case with vegetation (summer 2012), and finally, the alternative (new) proposal. The contours of the modulus of velocity and of the kinetic energy of turbulence in the characteristic horizontal plane (z¼ 3.5 m, beneath the vegetation canopy, which starts at z¼5 m) for analyzed scenarios are shown in Fig. 11. Building cross-sections are shown in white. It can be seen that the presence of vegetation reduces significantly the intensity of wind at desired location (500 rx r600 m and 1400 r y r1550 m) (Fig. 11-left). Here, the vegetation cross-sections are marked by black colored rectangles. The characteristic side jet is practically not visible anymore for the case of new vegetation, indicating strong suppression of the wind intensity (from jV j ¼ 4 m=s to less than 1.5 m/s) (Fig. 11-left/bottom). This suppression is also clearly visible in the contours of the turbulence kinetic energy (Fig. 11-right). The reorganization of flow and turbulence is even more visible
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
51
Fig. 11. Contours of the modulus of velocity (left) and superimposed velocity vectors and contours of the turbulence kinetic energy (right) in the characteristic horizontal plane beneath the vegetation islands (z ¼ 3.5 m) for three considered cases: without vegetation (top), the current situation (middle) and a new alternative proposal (bottom). The approaching wind direction is from top to bottom (indicated by the arrow).
52
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
Fig. 12. Same as in Fig. 11, only now at z ¼10 m (the plane crosses vegetation islands that are marked by black rectangles).
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
53
Fig. 13. The contours of the turbulence kinetic energy in the characteristic horizontal plane (z ¼10 m) and indications of the characteristic lines (L1–L2–L3– L4) along which the profiles are extracted. The presented solution is for the new arrangement of the vegetation islands with a ¼3 m2/m3. Note that (L5) coincides with (L4) but is extracted at z ¼3.5 m.
Fig. 15. Same as in Fig. 14, only now for contours of the turbulence kinetic energy.
Fig. 14. Contours of the modulus of velocity in the characteristic vertical plane along the L1 line (y¼ 1530 m) (a zoom-in) for: the reference state (no vegetation) (top), current situation (middle) and proposed new alternative (bottom). The direction of the approaching wind is perpendicular to the plot.
in the horizontal plane that crosses the vegetation islands at z¼ 10 m shown in Fig. 12. For this location, even more dramatic changes are observed. Again, the EWI building side jet is completely gone for the new vegetation (Fig. 12-left/bottom). This is achieved by adding an additional vegetation island to affect the wind conditions (this new vegetation island is located between 425 rx r525 m and 1600r y r 1700 m). It can be clearly observed that the velocity vectors within the vegetation canopy are significantly reduced (from 5 m/s to 1 m/s) (Fig. 12-right). To provide more detailed insights into intensity of the wind and its energy of fluctuations for different vegetation distributions, the profiles along the selected locations, shown in Fig. 13, are analyzed next. The characteristic lines (L1–L2–L3) are perpendicular whereas (L4–L5) are aligned with the mean wind direction. The first four lines are extracted at middle of the vegetation canopy height (z ¼10 m), while the last one (L5) coincides with (L4) and is extracted at height of z ¼3.5 m (just beneath vegetation which starts at z ¼5 m). The contours of the velocity and of the turbulence kinetic energy in the vertical plane across the L1 line (a zoom-in) are shown in Figs. 14 and 15. The direction of the approaching wind is perpendicular to this plane. It can be seen that the flow is significantly suppressed within the vegetation regions (Fig. 14). The contours shown in Fig. 15 show that the TKE is locally enhanced along the upper edges of the vegetation islands, while a very low turbulence intensity regions are observed within the vegetation. The profiles of modulus of velocity along the L1–L3 lines are shown in Fig. 16. The typical cross-section with vegetation islands are marked by green colored rectangles. Along the L3 line, there is practically no difference between the a¼0 and 3
54
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
. |V| [m/s]
6 4
2
3
2
3
a=0 m /m 2 3 a=3 m /m (new)
-L3-
a=3 m /m (old)
2 0
300
400
500
600
700
x [m]
|V| [m/s]
6 4
2
3
2
3
a=0 m /m 2 3 a=3 m /m (new)
-L1-
a=3 m /m (old)
2 0
300
400
500
600
700
x [m]
|V| [m/s]
6 2
3
2
3
a=0 m /m 2 3 a=3 m /m (new)
4
(old) m2/m3 vegetation cases. The situation is significantly different when additional vegetation is placed here (425r x r 525 m), which dramatically reduced the intensity of the wind (less than 1 m/s within the vegetation canopy). Note that the zero values correspond to segments that cross the buildings. Similarly, along the L1 and L2 lines, the new vegetation arrangement produces a dramatic reduction in the wind intensity. Almost a complete suppression is now obtained along the L2 (Fig. 16-bottom). This is the result of the combined effects of the extended vegetation islands located here, but also from the initial suppression of the incoming wind by placing additional vegetation island at the L3. It stresses importance of the taking into account combined effects of the different vegetation islands within the complex urban areas (Fig. 17). Finally, the profiles along the L4 and L5 are shown in Fig. 18. Now, the approaching wind is blowing from right to left, as indicated by the upper arrow. At z¼10 an 3.5 m, peaks in jVj are observed at the start of the first vegetation cross-section. It is followed by rapid (L4) or gradual (L5) decrease within the vegetation region. This is the consequence of the absence of the vegetation resistance along the L4 (just beneath the vegetation canopy). Interesting behavior is
-L2-
a=3 m /m (old)
2 0
300
400
500
600
700
x [m] Fig. 16. The horizontal profiles of the modulus of velocity along the characteristic locations indicated in Fig. 13 (i.e. in the z ¼ 10 m plane). The locations of the crosssections with vegetation islands are marked by green rectangles. The direction of the approaching wind is perpendicular to the plot. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)
Fig. 17. Same as in Fig. 16 only now for the turbulence kinetic energy profiles.
Fig. 18. The characteristic profiles of the modulus of velocity and of turbulence kinetic energy along lines L4 (z ¼10 m) and L5 (z ¼3.5 m) (shown in Fig. 13) for different leaf area densities (a¼ 0, 1 and 3 m2/m3). Note that the direction of the approaching wind is from right to left. The location of the line L5 is just beneath the line L4 (Fig. 13).
S. Kenjereš, B. ter Kuile / J. Wind Eng. Ind. Aerodyn. 123 (2013) 43–55
observed to the TKE profile along the L5 line (Fig. 18-bottom). Here, the TKE is even enhanced beneath the vegetation at the first location (starting at y¼1700 m) and then is significantly suppressed (more than fourfold reduction) in the second segment (starting at y¼1475 m). The results obtained with a reduced leaf area density of a¼ 1 m2/m3 for the new proposal are also shown for locations L4 and L5. It can be seen that just some marginal differences are obtained – apart from the TKE profile at z¼ 3.5 m, along the first green segment (1600 r yr 1700 m). 5. Summary and conclusions In this study, we analyzed different variants of the RANS twoequation k ɛ turbulence model for proper treatment of the vegetation in urban areas. Starting with a homogeneously distributed vegetation in a developing turbulent boundary layer, the additional terms in momentum, and equations for turbulence variables, are analyzed. Then a turbulent flow over a non-homogeneous (in the vertical direction) forest island is studied. Based on comparative assessment of the different models with the available experimental data, it is shown that the overall best agreement was obtained with the model that includes four additional terms in the turbulence equations. Simplified versions of the model were not able to properly mimic the changes in the velocity fields, and consequently in the turbulence kinetic energy, for cases with vertically non-homogeneously distributed vegetation. Then, the selected version of the model was applied for a real case of flow in a complex urban area (TU Delft campus). Simulations for three different scenarios are conducted. These included: no-vegetation case, present (old) vegetation situation (summer 2012) and a new proposal. It is shown that the velocity components are suppressed inside the vegetation regions. In contrast to that, the turbulence kinetic energy is firstly enhanced at the edges of these regions, and then significantly suppressed. This gives some practical guidelines for possible suppressions of the very strong side-winds along street-canyons in the complex urban areas. It is concluded that the vegetation distribution mimicking the situation of TU Delft campus in summer of 2012 partially affected flow and turbulence distributions within the campus site. In contrast to this, the intensity of the notorious side-wind gust around the tallest building at the campus site was significantly reduced with a new distribution of the vegetation islands. Finally, it is concluded that the variant proposed here for the non-homogeneous vegetation model proved to be numerically robust and physically accurate, making it a potentially useful tool for analysis and assessments of sustainable solutions in urban areas. References Ayotte, K.W., Finnigan, J.J., Raupach, M.R., 1999. A second-order closure for neutrally stratified vegetative canopy flows. Boundary-Layer Meteorology 90 (2), 186–216. Balczo, M., Gromke, C., Ruck, B., 2009. Numerical modeling of flow and pollutant dispersion in street canyons with tree planting. Meteorologische Zeitschrift 18 (2), 197–206. Belcher, S.E., Harman, I.N., Finnigan, J.J., 2012. The wind in the willows: flows in forest canopies in complex terrain. Annual Review of Fluid Mechanics 44, 479–504.
55
Dupont, S., Brunet, Y., 2009. Coherent structures in canopy edge flow: a large-eddy simulation study. Journal of Fluid Mechanics 630, 93–128. Durbin, P.A., 1991. Near-wall turbulence closure modeling without damping functions. Theoretical and Computational Fluid Dynamics 3, 1–13. Durbin, P.A., 1995. Separated flow computations with the k ɛ v2 model. AIAA Journal 33, 659–664. Green, S.R., 1992. Modeling turbulent air flow in stand of widely spaced trees. Phoenics Journal 23, 294–312. Gromke, C., Buccolieri, R., Di Sabatino, S., Ruck, B., 2008. Dispersion study in a street canyon with tree planting by means of wind tunnel and numerical investigations – evaluation of CFD data with experimental data. Atmospheric Environment 42, 8640–8650. Hanjalić, K., Kenjereš, S., 2005. Dynamic simulation of pollutant dispersion over complex urban terrains: a tool for sustainable development, control and management. Energy 30 (8), 1481–1497. Hanjalić, K., Kenjereš, S., 2006. RANS-based VLES of thermal and magnetic convection at extreme conditions. Journal of Applied Mechanics – Transactions of the ASME 73 (3), 430–440. Hanjalić, K., Kenjereš, S., 2008. Some developments in turbulence modelling for wind and environmental engineering. Journal of Wind Engineering & Industrial Aerodynamics 96 (10–11), 1537–1570. Katul, G.G., Mahrt, L., Poggi, D., Sanz, C., 2004. One- and two equation models for canopy turbulence. Boundary-Layer Meteorology 113, 81–109. Kenjereš, S., Hanjalić, K., Gunarjo, S.B., 2001. A T-RANS/LES approach to indoor climate simulations. In: Proceedings of ASME Fluids Engineering Division, Paper No. FEDSM20002-31400, Montreal, Canada. Kenjereš, S., Hanjalić, K., 2002. Combined effects of terrain orography and thermal stratification on pollutant dispersion in a town valley: a T-RANS simulation. Journal of Turbulence 3 (026), 1–25. Kenjereš, S., Gunarjo, S.B., Hanjalić, K., 2005. Contribution to elliptic relaxation modelling of turbulent natural and mixed convection. International Journal of Heat and Fluid Flow 26 (4), 569–586. Kenjereš, S., Hanjalić, K., 2006. LES, T-RANS and hybrid simulations of thermal convection at high Ra numbers. International Journal of Heat and Fluid Flow 27 (5), 800–810. Kenjereš, S., Hanjalić, K., 2009. Tackling complex turbulent flows with transient RANS. Fluid Dynamics Research 41 (1), 1–32. (article number: 012201). Kim, S.E., Choudhury, D., 1995. A near-wall treatment using wall functions sensitized to pressure gradient. In: Separated and Complex Flows, ASME FED vol. 217, ASME. Kuhn, S., Rudolf von Rohr, P., Kenjereš, S., 2010. Large eddy simulations of wall heat transfer and coherent structures in mixed convection over a wavy wall. International Journal of Thermal Sciences 49 (7), 1209–1226. Launder, B.E., Spalding, D.B., 1974. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering 3, 269–289. Liang, L., Xiaofeng, L., Borong, L., Yingxin, Z., 2006. Improved k ɛ two-equation turbulence model for canopy flow. Atmospheric Environment 40, 762–770. Liu, J., Chen, J.M., Black, T.A., Novak, M.D., 1998. E–ɛ modelling of turbulent air flow down wind of a model forest edge. Boundary-Layer Meteorology 72, 21–44. Mochida, A., Lun, I.Y.F., 2008. Prediction of wind environment and thermal comfort at pedestrian level in urban area. Journal of Wind Engineering and Industrial Aerodynamics 96, 1498–1527. Mochida, A., Tabata, Y., Iwata, T., Yoshino, H., 2008. Examining tree canopy models for CFD prediction of wind environment at pedestrian level. Journal of Wind Engineering and Industrial Aerodynamics 96, 1667–1677. Raupach, M.R., Thom, S.M., 1981. Turbulence in and above plant canopies. Annual Review of Fluid Mechanics 13, 97–129. Raupach, M.R., Shaw, R.H., 1982. Averaging procedures for flow within vegetation canopies. Boundary-Layer Meteorology 22, 79–90. Sanz, C., 2003. A note on k ɛ modelling of vegetation canopy air-flows. BoundaryLayer Meteorology 108, 191–197. Svensson, U., Häggkvist, 1990. A two-equation turbulence model for canopy flows. Journal of Wind Engineering and Industrial Aerodynamics 35, 201–211. Wagner, C., Kenjereš, S., Rudolf von Rohr, P., 2011. Dynamic large eddy simulations of momentum and wall heat transfer in forced convection over wavy surfaces. Journal of Turbulence 12 (7), 1–27. Wilson, J.D., 1988. A second-order closure model for flow through vegetation. Boundary-Layer Meteorology 42, 371–392. Wilson, J.D., Finnigan, J.J., Raupach, M.R., 1998. A first-order closure for disturbed plant-canopy flows, and its application to winds in a canopy on a ridge. Quarterly Journal of the Royal Meteorological Society 124, 705–732. Yi, C., 2008. Momentum transfer within canopies. Journal of Applied Meteorology and Climatology 47, 262–275.