- Email: [email protected]

Contents lists available at SciVerse ScienceDirect

Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia

Inﬂow turbulence generation for large eddy simulation in non-isothermal boundary layers Guoyi Jiang a,n, Ryuichiro Yoshie b, Taichi Shirasawa c, Xinyang Jin d a

Tokyo Polytechnic University, Atsugi, Kanagawa, Japan Tokyo Polytechnic University, Atsugi, Kanagawa, Japan c Otsuma Women’s University, Tama, Tokyo, Japan d China Academy of Building Research, Beijing, China b

a r t i c l e i n f o

abstract

Available online 21 March 2012

Many studies have presented the importance of turbulent inﬂow data for large eddy simulation, and several techniques have been proposed for generating inﬂow ﬂuctuations for LES in a neutral boundary layer, but few studies have considered non-isothermal cases. In this study, ﬁrst wind tunnel experiments were conducted to generate both unstable and stable turbulent thermal boundary layers in a thermally stratiﬁed wind tunnel at Tokyo Polytechnic University. The decrease of turbulence due to the effect of buoyancy was clearly observed in a stable boundary layer. Then a precursor method using a long domain and a recycling procedure using a short domain to generate the turbulence statistics corresponding to the experimental data were investigated for both unstable and stable conditions. The characteristics of the generated ﬂow (mean proﬁles and ﬂuctuation proﬁles) by both methods agreed relatively well with those of the wind tunnel experimental data. When a recycling procedure was adopted to generate turbulence statistics in a stable boundary layer, large damping was necessary to decrease the turbulence level. For the temperature ﬁeld, giving only a mean proﬁle (no temperature ﬂuctuation) as an inﬂow boundary condition of the driver region was the simplest way. Even if the inﬂow temperature ﬂuctuation was not given, temperature ﬂuctuation gradually developed in the downstream direction with velocity ﬂuctuation in the turbulent boundary layer. Crown Copyright & 2012 Published by Elsevier Ltd. All rights reserved.

Keywords: Non-isothermal Unstable Stable LES Inﬂow turbulence

1. Introduction Turbulent ﬂows are complex states of ﬂuid motion, and are characterized by eddies with a wide range of length and time scales. When simulating turbulent atmospheric boundary layers using Large Eddy Simulation (LES), a crucial issue is how to impose physically correct ﬂuctuating inﬂow data. The incoming ﬂow should have these spatial and temporal characteristics. The inﬂuence of inﬂow turbulence on large eddy simulation has been presented by Tominaga et al. (2008). It has been conﬁrmed that the inﬂow turbulence for LES is extremely important. Several techniques have been proposed for generating inﬂow turbulence for LES in a neutral boundary layer. Most of them can be classiﬁed into four basic categories: random noise (white noise), synthetic method, precursor simulation, and recycling method. The simplest way is to superimpose random ﬂuctuations on the mean velocity proﬁle with amplitude determined by the n

Corresponding author. Tel./fax: þ 81 46 242 9545. E-mail addresses: [email protected] (G. Jiang), [email protected] (R. Yoshie), [email protected] (T. Shirasawa), [email protected] (X. Jin).

turbulence intensity level. But due to the lack of the required characteristics of turbulent ﬂow, random noise is not an appropriate inlet condition. One commonly used way is the artiﬁcial generation method (synthetic method) in which velocity ﬂuctuations are obtained from an inverse Fourier transform of a prescribed spectrum that satisﬁes the characteristics of power spectral density and spatial correlations of the turbulent boundary layer. A method of generating the inﬂow turbulence for CWE applications based on this approach was developed by Kondo et al. (1997), Maruyama et al. (1999), and Iizuka et al. (1999). A digital-ﬁlter-based generation of turbulent inﬂow conditions was presented by Xie and Castro (2008), and the method was validated by simulating plane channel ﬂows with smooth walls and ﬂows over arrays of staggered cubes (a generic urban-type ﬂow). Another synthetic method for generating realistic inﬂow conditions was presented by Druault et al. (2004) based on proper orthogonal decomposition (POD) and linear stochastic estimation (LSE). But its spatial resolution is not sufﬁcient due to the limited number of hot wire probes that simultaneously measure velocity. Perret et al. (2006) took the same approach based on the use of POD, but they coupled it with a database obtained by stereoscopic particle image velocimetry (SPIV).

0167-6105/$ - see front matter Crown Copyright & 2012 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jweia.2012.02.030

370

G. Jiang et al. / J. Wind Eng. Ind. Aerodyn. 104–106 (2012) 369–378

Nomenclature f /fS x, y, z ui L

y d

instantaneous value of a quantity time-averaged value of f three components of space coordinates (stream-wise, span-wise, vertical) (m) three components of velocity vector (m/s) distance between aluminum plates (0.2 m) temperature (1C) boundary layer thickness (unstable 0.25 (m); stable 0.2 (m))

Another method that has been used to generate inﬂow conditions involves running a separate precursor calculation of an equilibrium ﬂow to generate a library of turbulent data that can be introduced into the main computation at the inlet (Tabor and Baha-Ahmadi, 2010). This has the advantage that the inﬂow conditions for the main computation are taken from a genuine simulation of turbulence, and thus should possess many of the required characteristics, including temporal and spatial ﬂuctuation with correlation and a correct energy spectrum. Instead of simulating the entire upstream region using precursor simulation, the most commonly used way is to employ a recycling method in which the velocity in a downstream plane is used (recycled) for the inﬂow boundary. In this method, the region used to generate turbulence is usually referred to as ‘‘driver region’’, and the second region that we are interested in is usually referred to as ‘‘main region’’. Lund et al. (1998) ﬁrst proposed this rescaling recycling method to generate developing turbulent inﬂow data for LES. Kataoka and Mizuno (2002) simpliﬁed Lund’s method by assuming that the boundary layer thickness is constant within the driver section. Nozawa and Tamura (2002, 2005) extended Lund’s method to a rough-wall boundary layer ﬂow using a roughness block arrangement. Now these methods are widely used to generate inﬂow ﬂuctuations for LES in neutral boundary layers. The non-isothermal boundary layer (unstable or stable) is a very common atmospheric phenomenon, but there have been few studies on stratiﬁed atmospheric ﬂows. When LES is applied to a non-isothermal ﬁeld, not only inﬂow velocity ﬂuctuation but also temperature ﬂuctuation is necessary. Kong et al. (2000) proposed a method for generating inﬂow ﬂuctuation for temperature with reference to Lund’s method according to the similarity between temperature and stream-wise velocity. To the knowledge of the author, this is the ﬁrst time inﬂow turbulence in the nonisothermal condition has been dealt with. In this method, the velocity ﬂuctuation was generated using Lund’s method, and the inﬂow temperature ﬂuctuation was generated using the same rescaling and recycling law as that used for stream-wise velocity. Hattori et al. (2007) conducted an unstable and stable boundary layer simulation, and the turbulent inﬂow data were generated by this method (Lund et al., 1998; Kong et al., 2000). In their calculation, the bulk Richardson number was 0.01 for the unstable condition and 0.01 for the stable condition, which corresponds to very weak thermal stratiﬁcation. In the driver region, the neutral boundary layer was simulated, and the temperature was treated as a passive scalar. Tamura et al. (2003) proposed a method for dealing with the thermally stratiﬁed effect. In the driver region, velocity ﬂuctuation was generated using the quasi-periodic boundary condition for a rough wall, while temperature was treated as a passive scalar, and a mean temperature proﬁle was given to the inﬂow condition of the driver region. Generated inﬂow data for temperature as well as velocity were introduced into the main

Ud

Yd Yf Y0 DY k Ri

mean stream-wise velocity at boundary layer height mean temperature at boundary layer height ﬂoor temperature (unstable 45 (1C); stable 18 (1C)) vertically averaged air temperature along boundary layer thickness temperature difference, ¼9Yf Yd9 (unstable 36 (1C); stable 34 (D)) turbulent kinetic energy (m2/s2) local gradient Richardson number

computation domain, where the solution of physical quantities took into account buoyancy effects. But they also pointed out that ﬂuctuation of the passive scalar was too large and not appropriate for the inﬂow condition of a stable turbulent boundary layer. Therefore, for stable conditions, the generated inﬂow velocity data at the recycling station were introduced into the main computational region, but a mean temperature proﬁle without ﬂuctuation was imposed at the inﬂow boundary of the main region (Tamura, 2008). Tamura’s treatments for dealing with the thermally stratiﬁed effect are shown in Fig. 1. In the simulation of the urban heat island phenomena of the downtown region of Tokyo, as the ﬂow coming from the sea was cold air, they used two driver regions to generate inﬂow turbulence for LES (Tamura et al., 2006b). The driver region in domain 1 generated a neutral turbulent boundary layer by using the rescaling technique, and domain 2 thermally stabilized the turbulent boundary layer based on the sea breeze characteristics. Then the generated data were used for domain 3 to simulate the urban heat island phenomena of Tokyo. Abe et al. (2008) conducted LES calculation of gas dispersion in a convective boundary layer (CBL), and investigated the characteristics of turbulence structures and gas dispersion behavior. The inﬂow velocity and temperature ﬂuctuations were generated using the method described by Tamura et al. (2003) for unstable boundary layers. Brillant et al. (2008) also developed a thermal turbulent inﬂow condition based on parallel ﬂows in order to simulate a turbulent thermal boundary layer. Then they tested this thermal turbulent inﬂow condition through a turbulent plane channel simulation. In

Fig. 1. Tamura’s treatments for dealing with thermally stratiﬁed effect. (a) Unstable boundary layer; (b) stable boundary layer.

G. Jiang et al. / J. Wind Eng. Ind. Aerodyn. 104–106 (2012) 369–378

371

temperature ﬂuctuation was not given, temperature ﬂuctuation gradually developed toward the downstream direction with velocity ﬂuctuation in the turbulent boundary layer.

Pre-simulation domain

u, v, w (y, z, t) (y, z, t) Data base

Main domain Fig. 2. Schematic view of two domains for LES.

this simulation, when velocity and temperature inﬂow ﬂuctuations were given simultaneously, the ﬂuctuating temperature proﬁles were well maintained as the ﬂow proceeded downstream. If velocity inﬂow ﬂuctuation and only a mean temperature proﬁle (no ﬂuctuation) were given as inﬂow condition, the temperature ﬂuctuation gradually developed as the ﬂow proceeded downstream. But the turbulent velocity ﬁeld did not quickly generate thermal ﬂuctuations. The present authors used a precursor method to generate velocity and temperature ﬂuctuations simultaneously in a nonisothermal boundary layer (Yoshie et al., 2011). In this method, the total simulations were composed of two domains: a presimulation domain to generate inﬂow velocity and temperature ﬂuctuations for LES, and a main domain to simulate gas dispersion around a high-rise building in a non-isothermal boundary layer (shown in Fig. 2). In the pre-simulation, the whole wind tunnel and all the aluminum plates (roughness elements) were reproduced by large eddy simulation using a buoyant solver. However, these kinds of researches on large eddy simulation of non-isothermal turbulent boundary layers are still quite rare, and methods of generating temperature ﬂuctuation have not been sufﬁciently examined yet. As a continuation of our previous study, which only focused on an unstable case using a precursor method (Yoshie et al., 2011), this study investigated two inﬂow turbulence generation methods (precursor method and recycling method) in both unstable and stable boundary layers. Wind tunnel experiments were conducted to generate both unstable and stable turbulent thermal boundary layers in a thermally stratiﬁed wind tunnel at Tokyo Polytechnic University. The enhancement of turbulence in the unstable condition and suppression of turbulence in the stable condition due to the effect of buoyancy was clearly observed in the wind tunnel experiment. The measured local gradient Richardson number near the boundary layer height was 0.3 for the unstable condition and 0.23 for the stable condition, which corresponded to weak thermal stratiﬁcation. Then a precursor simulation method using a long domain and a recycling procedure (Kataoka’s method) using a short domain to generate turbulence statistics corresponding to these experimental data were investigated for both unstable and stable thermal stratiﬁcations. In this precursor simulation, the whole wind tunnel and all the roughness elements were reproduced by large eddy simulation using a buoyant solver. The boundary layer proﬁles generated by both precursor simulation and the recycling method showed relatively good agreement with the wind tunnel experimental data. When the recycling procedure was adopted to generate turbulence statistics in a stable boundary layer, a large damping effect by a damping function was necessary for Kataoka’s method to decrease the turbulence level accordingly. For the temperature ﬁeld, only giving a mean proﬁle (no temperature ﬂuctuation) as inﬂow boundary condition of the driver region was the simplest way. Even if the inﬂow

2. Wind tunnel experiment Fig. 3 shows the experimental setup. The experiments were conducted in a thermally stratiﬁed wind tunnel at Tokyo Polytechnic University. The unstable and stable turbulent boundary layers were generated by 26 very thin aluminum plates 9 mm high, placed at L¼200 mm pitch on the wind tunnel ﬂoor. In order to form the thermal boundary layers, the oncoming air at the wind tunnel inlet was pre-cooled (9 1C) for the unstable condition and pre-heated (52 1C) for the stable condition. The ground temperature was set at 45 1C for the unstable condition, and at 18 1C for stable condition. The measuring station was located 5.6 m downstream from the inlet (0.1 m behind the last aluminum plate), where the turbulence and boundary layer were already fully developed. A split-ﬁlm and a cold-wire were used to measure wind velocity and temperature simultaneously (Yoshie et al., 2007). The characteristics of the generated boundary layers at the measuring station are summarized in Table 1. Both the mean velocity and mean temperature proﬁles closely followed a logarithmic law. d is the boundary layer thickness: d ¼0.25 m for the unstable condition and d ¼0.2 m for the stable condition. Ud and Yd are the mean streamwise velocity and mean temperature at the boundary layer height, and Ud ¼1.53 m/s, Yd ¼9 1C for the unstable condition, Ud ¼1.48 m/ s, Yd ¼52 1C for the stable condition. Y0 is the average temperature from the ground to the boundary layer height, and Y0 ¼16 1C for the unstable condition and 42 1C for the stable condition. Local gradient Richardson number Ri is the most widely used indicator of stability in the atmospheric boundary layer, and is deﬁned as:

ðg=Y0 Þ @ y [email protected] Ri ¼ ð1Þ 2 @hu1 [email protected] where g is acceleration due to gravity (9.8 m/s2), /u1S is mean stream-wise velocity proﬁle, and /yS is mean temperature proﬁle.

Fig. 3. Experimental setup and roughness elements arrangement.

Table 1 Characteristics of generated boundary layers.

d (m) Ud (m/s) Yd (D) Yf (D) DY (D) Y0 (D)

Unstable

Stable

0.25 1.53 9 45 36 16

0.2 1.48 52 18 34 42

372

G. Jiang et al. / J. Wind Eng. Ind. Aerodyn. 104–106 (2012) 369–378

The magnitude of Ri gradually increased along the vertical direction, and reached 0.3 near the boundary layer height for the unstable condition and 0.23 for the stable condition, which corresponded to weak thermal stratiﬁcation.

3. Computational method Two techniques for generating inﬂow ﬂuctuations for LES in a non-isothermal boundary layer (both unstable and stable conditions) were investigated in this study: precursor simulation using a long domain and a recycling procedure (Kataoka’s method) using a short domain. The generated turbulence characteristics were compared with the wind tunnel experimental data. 3.1. Precursor simulation The method for generating both velocity and temperature ﬂuctuations simultaneously in non-isothermal boundary layers by precursor simulation is discussed here. In the previous study by Ohya and Uchida (2008), they simulated a very long atmospheric boundary layer, and investigated the inﬂuence of thermal stability on turbulence. We also simulated a long spatially developing turbulent boundary layer (including the transition from laminar ﬂow to turbulent ﬂow), but focused on the purpose of generating inﬂow turbulence for LES. This method is referred to as ‘‘precursor method’’ hereafter. The whole wind tunnel (6.5 m long) and all the aluminum plates (shown in Fig. 3) were reproduced by LES using a buoyant solver. The plates were treated as zero-thickness plates in the simulation. The wind velocity and temperature distributions at the wind tunnel inlet were spatially uniform and turbulence intensity was very small (less than 1%), so a uniform velocity U¼ 1.42 m/s (unstable), 1.4 m/s (stable) and a uniform temperature Y ¼9 1C (unstable), 52 1C (stable) without turbulence were given to the inﬂow boundary of the precursor simulation. A zero gradient condition was used for the outlet boundary condition. The depth of the ﬁrst ﬂuid cells on the ﬂoor was 0.1 mm. A no-slip boundary condition was applied to the wall shear stress on the ﬂoor. The non-dimensional distances from the surfaces to the ﬁrst ﬂuid cells were less than 1.0 for most regions. As thermal boundary conditions, the surface temperature was 45 1C for the unstable condition and 18 1C for the stable condition, and a Fourier’s law of heat conduction was applied for the heat ﬂux on the ﬂoor surface. A wall function of Spading law (1961) was used for the two lateral boundaries (side walls of the wind tunnel) and the top boundary (ceiling of the wind tunnel). Three mesh systems (here referred to as M-1, M-2 and M-3) were tested for the precursor simulations to examine the inﬂuence of mesh distributions (stream-wise and vertical direction). They are summarized in Table 2. The mesh distribution for M-1 and M-2 is the same in the span-wise direction (Dy¼10 mm). But M-1 has more mesh numbers in the stream-wise direction (Dx¼6 mm), while the mesh density of M-2 is relatively coarse in the stream-wise direction (Dx ¼10 mm). Although the vertical Table 2 Mesh systems used for precursor simulation. Mesh system

Dx, Dy

M-1

6, 10

M-2

10, 10

M-3

6, 10

Nx Ny Nz

Vertical mesh distribution

(mm) 1042 120 80 50 points (zo 0.3 m)þ 30 points (z 40.3 m) 651 120 80 61 points (zo 0.3 m) þ19 points (z 40.3 m) 1042 120 80 61 points (zo 0.3 m) þ19 points (z 40.3 m)

mesh numbers of M-1 and M-2 are the same (80 points), their distributuions are different. For M-2, ﬁner grids were given below the height z¼0.3 m. M-3 has the same mesh distributions as M-2 in the cross-section (y and z directions), and has the same mesh distribution in the stream-wise direction as M-1. For the unstable condition, the standard Smagorinsky sub-grid scale model (SM model) with Smagorinsky constant 0.12 was applied, and the Van Driest-type damping function was used to account for the near wall effect. In addition, the SGS Prandtl number was set to 0.7. For stable boundary layer simulation, both the standard Smagorinsky model (SM model) and the dynamic Smagorinsky model (DSM model) were examined. For the DSM model, the model coefﬁcient CS was determined using a dynamic procedure, but the value was locally averaged over cell faces in order to avoid numerical instabilities, and the negative values of the effective eddy viscosity were removed by clipping it to zero (neff ¼ nSGS þ n ¼0, in which nSGS is the sub-grid scale (SGS) eddy viscosity and n is the molecular viscosity). A second-order backward differencing scheme (using the current and two previous time-step values) was adopted for time integration. The second order central difference scheme was adopted for the convection terms. The PISO algorithm (Pressure Implicit with Splitting of Operators) was used for the pressure– velocity calculation procedure. PISO involves one predictor step and two corrector steps. At each time step, velocity and temperature were predicted, and then pressure and velocity were corrected. The PISO algorithm was adopted because it has better performance in the simulation of unsteady ﬂow (Issa, 1986). For the unstable condition, SM model was examined using M-1 and M-2 meshes. For the stable condition, SM model was examined using M-1 and M-2 meshes, and DSM model was examined for all the three mesh systems. A unique time step Dt ¼0.001 s for mesh M-1 and M-3 and Dt¼0.002 s for mesh M-2 was used to make sure that the maximum courant number was less than 1 in all positions of the domain. The sampling station for obtaining ﬂuctuating velocity and temperature data was put 0.1 m downstream of the last aluminum plate (the same position as the experimental measuring station), where the turbulence and the boundary layer were already fully developed. The data at this section can be stored and used for main simulation. These calculations were conducted by a high-performance PC with 2 CPU 6 cores (12 cores of parallel simulation). Fig. 4 shows the vertical distributions of turbulent kinetic energy and the r.m.s. value of temperature ﬂuctuation at the sampling station from the results using different mesh systems and different SGS models. For the unstable condition (Fig. 4(a) and (b)), the difference between the two LES results is relatively small, and the results are not sensitive to the mesh system. Hereafter, only the results using mesh system M-1 and SM model are shown in a later section for the unstable condition. For stable boundary layer simulation (Fig. 4(c) and (d)), inaccurate results were obtained for both mesh systems when the SM model was applied. Of the ﬁve results, those obtained from the DSM model using the M-2 and M-3 mesh systems were reasonable, and the difference between the results of M-2 and M-3 was small. The results were sensitive to the SGS model and the mesh system for the stable condition, and the inﬂuence of mesh distributions in vertical direction was larger than that in stream-wise direction. More mesh numbers were needed inside the boundary layer for stable boundary layer simulation. Although large eddy simulation has gained a victory in studying unstable and neutrally stratiﬁed boundary layers, some researches (Beare and Macvean, 2004; Beare et al., 2006; Tamura and Mori, 2006a) have focused on the difﬁculties of LES with stably stratiﬁed turbulent boundary layers. Comparing unstable and neutral conditions, the dominant eddies are usually smaller in stable boundary

G. Jiang et al. / J. Wind Eng. Ind. Aerodyn. 104–106 (2012) 369–378

1.2

1.2 Exp M-1-SM M-2-SM

Exp M-1-SM M-2-SM

0.9

z/δ

z/δ

0.9

0.6

0.6

0.3

0.3

0.0 0.00

0.01

0.02

0.0 0.00

0.03

0.05

0.10

k/Uδ2 1.2 EXP M-1-SM M-2-SM M-1-DSM M-2-DSM M-3-DSM

0.6

Exp M-1-SM M-2-SM M-1-DSM M-2-DSM M-3-DSM

0.9

z/δ

0.9

z/δ

0.15

σθ/ΔΘ

1.2

0.3

0.0 0.00

373

0.6

0.3

0.01

0.02

0.03

0.0 0.00

0.05

0.10

k/Uδ2

0.15

σθ/ΔΘ

Fig. 4. Turbulence statistics obtained using different mesh systems and different SGS models. (a) Turbulent kinetic energy (unstable); (b) r.m.s. value of temperature ﬂuctuation (unstable); (c) turbulent kinetic energy (stable); (d) R.m.s. value of temperature ﬂuctuation (stable).

layers and the ﬂow is complicated due to the gravity waves. In this study, the precursor simulation results using the DSM model and the M-2 mesh system are shown for the stable condition in a later section. 3.2. Recycling method For this case, the velocity ﬂuctuation was generated using Kataoka’s method (Kataoka and Mizuno, 2002) with the roughness ground arrangement described by Nozawa and Tamura (2002). The roughness elements were exactly the same as those used in the precursor simulation, but a short domain (about 6d) was adopted here. A mean velocity proﬁle that came from the experimental measurement was prescribed for the inﬂow condition of the driver region, and only the ﬂuctuating part was recycled between outlet station and inlet station. The velocity components at the inﬂow boundary of the driver section are given as follows: uinlt ðy,z,tÞ ¼ huiinlt ðzÞ þ jðZÞ uðy,z,tÞhuiðy,zÞ recy vinlt ðy,z,tÞ ¼ jðZÞ vðy,z,tÞhviðy,zÞ recy winlt ðy,z,tÞ ¼ jðZÞ wðy,z,tÞhwiðy,zÞ recy ð2Þ The following damping function (Kataoka, 2008) was used to restrain the velocity ﬂuctuation from developing for the unstable

condition:

fðZÞ ¼ 12 1tanh

8:0ðZ1:0Þ =tanhð8:0Þ 0:4Z þ0:82

ð3Þ

where, Z ¼z/d. Due to suppression effect of turbulence in stable condition, the turbulence in the stable boundary layer becomes small. Therefore, when generating inﬂow turbulence using the recycling procedure (Kataoka’s method), this should be taken into account. Here, a modiﬁed damping function with a larger damping effect than Eq. (3) was adopted, and different damping functions were used for different velocity components for the stable condition: 8:0ðZ1:0Þ =tanhð8:0Þ ð4Þ fu,v ðZÞ ¼ 0:8 12 1tanh 0:4Z þ 0:82

fw ðZÞ ¼ 0:2 12 1tanh

8:0ðZ0:2Þ =tanhð8:0Þ 0:4Z þ0:82

ð5Þ

A no-slip condition was used for the ground surface, and a wall function was used for the two sides and the top boundaries. For the recycling method, the standard Smagorinsky sub-grid scale model (SM model) was used. The temperature was treated as a passive scalar, and no buoyancy term was included in the momentum equation. A mean temperature proﬁle for the experiment was prescribed at the inﬂow boundary of the driver region,

374

G. Jiang et al. / J. Wind Eng. Ind. Aerodyn. 104–106 (2012) 369–378

and we tried to use the ﬂuctuating velocity ﬁeld to generate a ﬂuctuating temperature ﬁeld. Total mesh number 150 120 55¼990,000 was used, which was much less than that used for the precursor simulation. The arrangement for the recycling procedure is shown in Fig. 5

4. Results and discussion ! Vorticity o is one of the most widely used functions for identifying coherent vortex structures and representing vortex cores, and is deﬁned as: !

o ¼r! u

ð6Þ

Recycling fluctuating part of velocity

Main region Buoyancy effect

Outlet

Driver region

Inlet

<θ>

θ

Sampling Plane

u

Neutral boundary layer

Bottom heating or cooling

Fig. 5. Arrangement for recycling procedure.

Fig. 6 shows the instantaneous vortex structures (9oy9¼1, 5, 10) from the results of precursor simulations. The turbulence and the boundary layer were gradually generated by these 26 aluminum plates (roughness elements) as the ﬂow proceeded downstream. Compared with the unstable condition (Fig. 6(a)), decrease of turbulence due to the effect of buoyancy was clearly observed in the stable boundary layer simulation (Fig. 6(b)), which is in accordance with the experimental results. Fig. 7 shows the distributions of mean stream-wise velocity, turbulent kinetic energy and the r.m.s. value of temperature ﬂuctuation at the sampling station for both experimental data and LES results (black: unstable; green: stable). It took a total of 15 s to perform time-averaging for the LES results, and we conﬁrmed that the deference of the averaged values for 10 s and 15 s was extremely small. All the variables were averaged along the span-wise direction (the regions near the side boundaries were excluded, and the averaging did not include these regions). The generated (calculated) mean stream-wise velocity proﬁles by precursor simulation and recycling method match very well with those of the experiment for both unstable and stable conditions (Fig. 7(a)). The experimental results showed that, compared with the unstable condition, the turbulent kinetic energy is smaller in the stable boundary layer (Fig. 7(b)). This is because the turbulence was enhanced in unstable condition and was suppressed in stable condition due to the effect of buoyancy. This feature was well captured by precursor simulations. Ohya

Fig. 6. Instantaneous vorticity magnitude of 9oy9¼ 1 (green), 5 (yellow), 10 (red). (precursor simulation, y¼ 0): (a) unstable; (b) stable. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

1.2

0.6

0.3

0.3

0.0 0.0

0.4

0.8

〈u1〉/Uδ

1.2

0.0 0.00

Exp_u Precursor_u Recycling_u Exp_s Precursor_s Recycling_s

0.9

z/δ

0.6

Exp_u Precursor_u Recycling_u Exp_s Precursor_s Recycling_s

0.9

z/δ

0.9

z/δ

1.2

1.2

Exp_u Precursor_u Recycling_u Exp_s Precursor_s Recycling_s

0.6

0.3

0.01

0.02

k/Uδ

2

0.03

0.0 0.00

0.05

0.10

0.15

σθ/ΔΘ

Fig. 7. Vertical distributions of mean stream-wise velocity, turbulent kinetic energy and r.m.s. value of temperature ﬂuctuation. In the legends, ‘‘u’’ means unstable, and ‘‘s’’ means stable. (a) Mean stream-wise velocity; (b) turbulent kinetic energy; (c) r.m.s. value of temperature ﬂuctuation.

G. Jiang et al. / J. Wind Eng. Ind. Aerodyn. 104–106 (2012) 369–378

and Uchida (2008) investigated the phenomena of stable boundary layer for a wide range of stabilities using wind tunnel experiment and numerical simulations. They found that stable stratiﬁcation rapidly suppresses the ﬂuctuations of velocity and temperature, and these turbulence statistics become nearly zero over the whole boundary depth with very strong stability. For the recycling method, a large damping effect (Eqs. (4) and (5)) was adopted here to take into account this effect when generating inﬂow turbulence for the stable boundary layer simulation. These damping functions were determined by trial and error. They were adjusted several times in the simulation until a reasonable turbulent kinetic energy proﬁle was obtained (dashed green line in Fig. 7(b)). It has to be pointed out that, for stronger unstable thermal condition we need to use smaller damping effect to increase the ﬂuctuations. For stronger stable thermal stratiﬁcation, we need larger damping effect to decrease the ﬂuctuations accordingly. Due to this low turbulence in the stable boundary layer, the temperature ﬂuctuations also became smaller in the experimental results (Fig. 7(c)). The calculated r.m.s. value of temperature ﬂuctuations by both precursor simulation and the recycling method agreed fairly well with the experimental data for the unstable condition, although these values were somewhat overestimated for the stable condition. The difference between sy values for the unstable and stable conditions is relatively small in the results of both precursor simulation and recycling method. Fig. 8 shows vertical distributions of the r.m.s. value of velocity ﬂuctuations (su, sv, sw). The decrease of these values due to the effect of buoyancy was clearly observed in the stable boundary layer. Overall, these variables generated by both precursor simulation and

0.9

z/δ

z/δ

0.6

1.2 Exp_u Precursor_u Recycling_u Exp_s Precursor_s Recycling_s

0.6

0.3

0.3

0.08

0.16

0.0 0.00

0.24

Exp_u Precursor_u Recycling_u Exp_s Precursor_s Recycling_s

0.9

z/δ

Exp_u Precursor_u Recycling_u Exp_s Precursor_s Recycling_s

0.9

0.0 0.00

the recycling method agreed relatively well with those of the wind tunnel experimental data, except that sw was somewhat overestimated by both of the LES results. Figs. 9 and 10 show the distributions of mean temperature proﬁle and turbulent heat ﬂuxes in the unstable and stable boundary layer, respectively. The calculated mean temperature by both methods agreed well with those of the experimental data, but the vertical turbulent heat ﬂux was overestimated by LES in both the unstable and stable boundary layers. The decrease of turbulent heat ﬂux due to the effect of buoyancy was also observed in both experiment and precursor simulation in the stable boundary layer from comparison of Figs. 9 and 10. Fig. 11 shows the instantaneous stream-wise velocity and temperature distributions from the results of the recycling method in the unstable boundary layer (symmetry plane, y¼0). L is the distance between the aluminum plates (200 mm), and x/L¼0 corresponds to the inlet position of the driver section. The black triangle shows the sampling station. As the recycling procedure was adopted for the velocity ﬁeld, the velocity was ﬂuctuating over the whole driver region (Fig. 11(a)). But only a mean temperature proﬁle was given as an inﬂow condition, so the temperature ﬂuctuation was very small near the inlet position. As the ﬂow proceeded downstream, the temperature ﬁeld was caused to ﬂuctuate gradually by turbulence (Fig. 11(b)). Fig. 12 shows the vertical distributions of the r.m.s. values of temperature ﬂuctuation and turbulent heat ﬂuxes at six different stations along the stream-wise direction of the driver region (recycling method) for the unstable (Fig. 12(a)–(c)) and stable (Fig. 12(d)– (f)) conditions. The temperature ﬂuctuation is very small near the

1.2

1.2

375

0.6

0.3

0.08

0.16

0.0 0.00

0.24

0.08

0.16

0.24

Fig. 8. Vertical distributions of r.m.s. values of velocity ﬂuctuations. (a) Stream-wise component; (b) span-wise component; (c) vertical component.

1.2 Exp Precursor Recycling

0.9

z/δ

z/δ

0.9

0.6

0.3

0.0 -1.2

1.2 Exp Precursor Recycling

0.6

0.3

-0.8

-0.4

0.0

0.0 -0.009

Exp Precursor Recycling

0.9

z/δ

1.2

0.6

0.3

-0.006

-0.003

0.000

0.0 0.000

0.001

0.002

0.003

Fig. 9. Vertical distributions of mean temperature and turbulent heat ﬂuxes (unstable). (a) Mean temperature; (b) stream-wise turbulent heat ﬂux; (c) vertical turbulent heat ﬂux.

376

G. Jiang et al. / J. Wind Eng. Ind. Aerodyn. 104–106 (2012) 369–378

1.2

1.2 Exp Precursor Recycling

Exp Precursor Recycling

0.6

0.3

0.9

z/δ

0.9

z/δ

z/δ

0.9

1.2

0.6

0.3

0.0 0.0

0.4

0.8

1.2

Exp Precursor Recycling

0.6

0.3

0.0 0.000

0.003

0.006

0.009

0.0 -0.003

-0.002

-0.001

0.000

4

4

3

3 z/δ

z/δ

Fig. 10. Vertical distributions of mean temperature and turbulent heat ﬂuxes (stable). (a) Mean temperature; (b) stream-wise turbulent heat ﬂux; (c) vertical turbulent heat ﬂux.

2 1 0

2 1

0 0.5

1.5

2.5

3.5 4.5 x/L

5.5

6.5

7.5

0

0 0.5

1.5

2.5

3.5 4.5 x/L

5.5

6.5

7.5

Fig. 11. Instantaneous variables in symmetry plane (unstable, recycling method, y¼ 0). (a) Instantaneous contour of stream-wise velocity; (b) instantaneous contour of temperature.

inﬂow boundary (black line) because no temperature ﬂuctuation was given to the inﬂow boundary, but as the ﬂow proceeded downstream, the temperature ﬂuctuation became larger and larger until position x/L¼4.5. Downstream of this position, a proﬁle that matched the experimental data was formed and nearly maintained. According to the results of Fig. 12, the method proposed by Tamura et al. (2003) (see Fig. 1(a)) is considered to be valid if sufﬁcient fetch distance is ensured before the target region. If the recycling procedure is adopted not only for velocity but also for temperature (like Kong’s method), the ﬂuctuation of temperature may develop faster. According to our result for the recycling method for the stable condition (Fig. 7(c)), the method for dealing with the unstable condition in the driver region (Fig. 1(a) left) proposed by Tamura et al. (2003) may also succeed in stable condition by combining a large damping function.

5. Conclusions In this study, we focused on the generation of inﬂow turbulence for large eddy simulation in non-isothermal boundary layers. First, wind tunnel experiments were carried out to generate both unstable and stable turbulent boundary layers in a thermally stratiﬁed wind tunnel at Tokyo Polytechnic University, using the simultaneous measuring techniques developed by the present authors. The decrease of turbulence due to the effect of buoyancy was clearly observed in the stable boundary layer. Then large eddy simulation was conducted to generate turbulence statistics corresponding to these experimental conditions. For LES, two methods were investigated here, including a precursor

simulation using a long domain and a recycling method using a short domain. For the precursor method, the whole wind tunnel and all the roughness elements were reproduced by LES using a buoyant solver. We simulated this long spatially developing turbulent thermal boundary layer (including transition from laminar ﬂow to turbulent ﬂow), but focused on the purpose of generating inﬂow turbulence for LES. As a recycling method, Kataoka’s method was examined here. The issue was discussed when it was adopted for the stable condition. A method for dealing with temperature ﬂuctuation was examined. The main conclusions of this study are as follows: (1) The characteristics of the generated ﬂow (mean proﬁles and ﬂuctuation proﬁles) by both precursor simulation and the recycling method agreed relatively well with those of the wind tunnel experimental data for both the unstable and stable conditions, and both of these methods can be used to generate inﬂow ﬂuctuations for LES in a non-isothermal boundary layer. (2) Because a long spatially developing boundary layer (including transition from laminar ﬂow to turbulent ﬂow) needs to be simulated in the precursor simulation, a mesh sensitivity check should be done ﬁrst. As the dominant eddies are usually small and the ﬂow ﬁeld is complicated in a stable boundary layer, a more elegant SGS model is necessary. Because no inﬂow turbulence is required for precursor simulation (only a uniform velocity and temperature value are given as inﬂow condition), this method is the simplest, but it requires a long computational domain and long simulation time. (3) For the recycling method, a short domain can be adopted as the driver region for generating turbulence. The requirement of

G. Jiang et al. / J. Wind Eng. Ind. Aerodyn. 104–106 (2012) 369–378

1.2

0.6

0.9

0.3

0.0 0.00

0.03

0.06

0.0 -0.009

0.09

z/δ

-0.006

-0.003

0.0 0.000

0.000

0.6

0.3

x/L=0.5 x/L=1.5 x/L=2.5 x/L=3.5 x/L=4.5 x/L=5.5

0.9

0.6

0.9

0.06

0.09

0.0 0.000

0.002

0.003

0.6

-0.001

0.000

x/L=0.5 x/L=1.5 x/L=2.5 x/L=3.5 x/L=4.5 x/L=5.5

0.3

0.3

0.03

0.001

1.2

z/δ

0.9

0.6

0.3

1.2 x/L=0.5 x/L=1.5 x/L=2.5 x/L=3.5 x/L=4.5 x/L=5.5

x/L=0.5 x/L=1.5 x/L=2.5 x/L=3.5 x/L=4.5 x/L=5.5

0.9

0.3

1.2

0.0 0.00

0.6

1.2 x/L=0.5 x/L=1.5 x/L=2.5 x/L=3.5 x/L=4.5 x/L=5.5

z/δ

z/δ

0.9

z/δ

x/L=0.5 x/L=1.5 x/L=2.5 x/L=3.5 x/L=4.5 x/L=5.5

z/δ

1.2

377

0.003

0.006

0.009

0.0 -0.003

-0.002

Fig. 12. Development of temperature ﬂuctuation in the driver region (recycling method). (a) R.m.s. value of temperature ﬂuctuation (unstable); (b) stream-wise turbulent heat ﬂux (unstable); (c) vertical turbulent heat ﬂux (unstable); (d) r.m.s. value of temperature ﬂuctuation (stable); (e) stream-wise turbulent heat ﬂux (stable); (f) vertical turbulent heat ﬂux (stable).

computational resources is relatively small. When Kataoka’s method is applied for inﬂow turbulence generation in a nonisothermal ﬁeld, a large damping effect is necessary to decrease the turbulence level accordingly in the stable boundary layer. For the temperature ﬁeld, giving a mean proﬁle as the inﬂow condition of the driver region is the simplest way (in both unstable and stable conditions), because the temperature ﬂuctuation develops as the ﬂow proceeds downstream in a turbulent boundary layer. Since the turbulent velocity ﬁeld does not generate thermal ﬂuctuations quickly, a sufﬁcient distance should be ensured for it to occur.

Acknowledgments This study was partially supported by the Ministry of Education, Culture, Sports, Science and Technology, Japan, through the Global Center of Excellence Program of Tokyo Polytechnic University. We would also like to express our gratitude to the Japan Society for the Promotion of Science (JSPS) for Grant-in-Aid for Scientiﬁc Research (B), (No. 21360283).

References Abe, S., Tamura, T., Nakayama, H., 2008. LES on turbulence structures and gaseous dispersion behavior in convective boundary layer with the capping inversion. Journal of Wind Engineering 33. Beare, R.J., Macvean, M.K., 2004. Resolution sensitivity and scaling of large-eddy simulations of the stable boundary layer. Boundary-Layer Meteorology 112, 257–281.

Beare, R.J., Macvean, M.K., Holtslag, A.M., Cuxart, J., Esau, I., Golaz, J.C., Jimenez, M.A., Khairoutdinov, M., Kosovic, B., Lewellen, D., Lund, T.S., Lundquist, J.K., Mccabe, A., Moene, A.F., Noh, Y., Raasch, S., Sullivan, P., 2006. An intercomparison of large eddy simulations of the stable boundary layer. Boundary-Layer Meteorology 118, 247–272. Brillant, G., Husson, S., Bataille, F., Ducros, F., 2008. Study of the blowing impact on a hot turbulent boundary layer using thermal large eddy simulation. International Journal of Heat and Fluid Flow 29, 1670–1678. Druault, P., Lardeau, S., Bonnet, J.P., Coiffet, F., Delville, J., Lamballais, E., Largeau, J.F., Perret, L., 2004. Generation of three-dimensional turbulent inlet conditions for large eddy simulation. AIAA Journal 42 (3), 447–456. Hattori, H., Houra, T., Nagano, Y., 2007. Direct numerical simulation of stable and unstable turbulent thermal boundary layers. International Journal of Heat and Fluid Flow 28, 1262–1271. Issa, R.I., 1986. Solution of the implicit discretized ﬂuid ﬂow equations by operator-splitting. Journal of Computational Physics 62, 40–65. Iizuka, S., Murakami, S., Tsuchiya, N., Mochida, A., 1999. LES of ﬂow past 2D cylinder with imposed inﬂow turbulence. Proceedings of 10th International Conference on Wind Engineering 2, pp. 1291–1298. Kataoka, H., Mizuno, M., 2002. Numerical ﬂow computation around aeroelastic 3D square cylinder using inﬂow turbulence. Wind and Structures 5, 379–392. Kataoka, H., 2008. Numerical simulations of a wind-induced vibrating square cylinder within turbulent boundary layer. Journal of Wind Engineering and Industrial Aerodynamics 96, 1985–1997. Kondo, K., Murakami, S., Mochida, A., 1997. Generation of velocity ﬂuctuations for inﬂow boundary condition of LES. Journal of Wind Engineering and Industrial Aerodynamics 67/&68, 51–64. Kong, H., Choi, H., Lee, J.S., 2000. Direct numerical simulation of turbulent thermal boundary layers. Physics of Fluids 12, 2555–2568. Lund, T.S., Wu, X.H., Squires, K.D., 1998. Generation of turbulent inﬂow data for spatially-developing boundary layer simulations. Journal of Computational Physics 140, 233–258. Maruyama, T., Rodi, W., Maruyama, Y., Hiraoka, H., 1999. Large eddy simulation of the turbulent boundary layer behind roughness elements using an artiﬁcially generated inﬂow. Journal of Wind Engineering and Industrial Aerodynamics 83, 381–392. Nozawa, K., Tamura, T., 2002. Large eddy simulation of the ﬂow around a low-rise building immersed in a rough-wall turbulent boundary layer. Journal of Wind Engineering and Industrial Aerodynamics 90, 1151–1162.

378

G. Jiang et al. / J. Wind Eng. Ind. Aerodyn. 104–106 (2012) 369–378

Nozawa, K., Tamura, T., 2005. Large eddy simulation of wind ﬂows over large roughness elements. Proceedings of EACWE4, 1–6. Ohya, Y., Uchida, T., 2008. Laboratory and numerical studies of the atmospheric stable boundary layers. Journal of Wind Engineering and Industrial Aerodynamics 96, 2150–2160. Perret, L., Delville, J., Manceau, R., Bonnet, J.P., 2006. Generation of turbulent inﬂow conditions for large eddy simulation from stereoscopic PIV measurements. International Journal of Heat and Fluid Flow 27, 576–584. Spading, D.B, 1961. A single formula for the law of the wall. Journal of Applied Mechanics, Transaction ASME, Series E 28, 455–458. Tabor, G.R., Baha-Ahmadi, M.H., 2010. Inlet conditions for large eddy simulation: a review. Computer & Fluids 39, 553–567. Tamura, T., Tsubokura, M., Cao, S., Furusawa, T., 2003. LES of spatially-developing stable/unstable stratiﬁed turbulent boundary layers. Direct and Large-Eddy Simulation 5, 65–66. Tamura, T., Mori, K., 2006a. LES of spatially-developing stably stratiﬁed turbulent boundary layers. Direct and Large-Eddy Simulation 6, 583–590. Tamura, T., Nakayama, J., Ohta, K., Takemi, T., Okuda, Y., 2006b. LES estimation of environmental degradation at the urban heat island due to densely-arrayed

tall buildings. In: 17th Symposium on Boundary Layers and Turbulence, AMS, pp. 22–25. Tamura, T., 2008. Towards practical use of LES in wind engineering. Journal of Wind Engineering and Industrial Aerodynamics 96, 1451–1471. Tominaga, Y., Mochida, A., Murakami, S., Sawaki, S., 2008. Comparison of various revised k–e models and LES applied to ﬂow around a high-rise building model with 1:1:2 shape placed within the surface boundary layer. Journal of Wind Engineering and Industrial Aerodynamics 96, 389–411. Xie, Z.T., Castro, I.P., 2008. Efﬁcient generation of inﬂow conditions for large eddy simulation of street-scale ﬂows. Flow, Turbulence and Combustion 81, 449–470. Yoshie, R., Jiang, G.Y., Shirasawa, T., Chung, J., 2011. CFD simulations of gas dispersion around high-rise building in non-isothermal boundary layer. Journal of Wind Engineering and Industrial Aerodynamics 99, 279–288. Yoshie, R., Tanaka, H., Shirasawa, T., 2007. Technique for simultaneously measuring ﬂuctuating velocity, temperature and concentration in non-isothermal ﬂow. In: Proceedings of the 12th International Conference on Wind Engineering, pp. 1399–1406.