Journal of Heat Island Institute International Vol. 12-2 (2017)

Academic Article

Large-Eddy Simulation of Non-Isothermal Flow around a Building
Using Artificially Generated Inflow Turbulent Fluctuations ofWind Velocity and Air Temperature
Tsubasa OKAZE *1

Akashi MOCHIDA *2*1 School of Environment and Society, Tokyo Institute of Technology*2 Graduate School of Engineering, Tohoku University
Corresponding author: Tsubasa Okaze, okaze.t.aa@m.titech.ac.jp
ABSTRACTThis study aims to apply artificially generated turbulent fluctuations of wind velocity and air temperature as
inflow boundary conditions for a large-eddy simulation of a non-isothermal flow around a building. The building
was an isolated high-rise structure with a length ratio of 1:1:2. The stratification condition of the flow field was
weakly unstable. The artificial turbulent fluctuations generated satisfied the prescribed profiles for the turbulent
fluxes of momentum and temperature. Moreover, they fulfilled the prescribed spatial and time correlations
expressed using exponential functions. The mean wind velocity and mean air temperature predicted from thelarge-eddy simulation agreed well with the experimental results. The turbulent kinetic energy and the fluctuation of
the air temperature also corresponded fairly well with the experiment except for the region immediately behind the
inflow boundary.
Key Words : Inflow turbulence generation, Temperature fluctuation, Large-eddy simulation,
Unstable stratification, Synthetic method, Cholesky decomposition

1. Introduction

region was established consisting of two domains in the
preliminary calculation for generating inflow turbulence and air

In recent years, large-eddy simulation (LES) computations thattemperature fluctuations. Matsuda et al.(7) investigated theconsider the atmospheric stability in built-up environments havegreenery effect in the Tokyo Bay area using LES. They set a

been reported. Generation of fluctuations of wind velocity andlong approaching section over the Tokyo Bay to generate the

air temperature fluctuation used as an inflow boundary condition

wind velocity and air temperature fluctuations. However, thehas become an important issue in implementing LES for

methods for generating the velocity and temperature fluctuations

non-isothermal flow fields. Kong et al.(1) proposed a method forbased on a preliminary or additional LES need additional

generating the inflow fluctuations of wind velocity and air

computational costs in the preliminary or additional calculations.

temperature. This method is based on the rescaling and recycling

Furthermore, in these methods, the distributions of mean wind(2)

laws for the flow field proposed by Lund et al.

and is extended

velocity and air temperature and its fluctuations are obtained

to the temperature field. Hattori et al.(3) conducted direct

after the preliminary computation. This means it is not easy tonumerical simulations for very weak unstable and stable

impose the prescribed profiles of wind velocity, air temperatureboundary layers using the method proposed by Kong et al.(1) toand the turbulent fluxes of the momentum and air temperature.

generate inflow turbulent fluctuations of wind velocity and air

Thus, the process of the trial and error is usually needed.

temperature. Jiang et al. (4) performed LESs for boundary layers,

Recently, the present authors proposed a method for generating

both weakly stable and weakly unstable stratifications, using

turbulent fluctuations of wind velocity and various scalar

Kataoka’s recycling method(5) to generate inflow turbulence andquantities,

air temperature fluctuations. Tamura et al.

(6)

such

as

air

temperatureandcontaminantconducted an LES

concentration, based on Cholesky decomposition of thefor a non-isothermal flow field in the Tokyo metropolitan areatime-averaged turbulent flux tensors of momentum and scalarwith inflow turbulence and air temperature fluctuations. A driver

quantities(8). The artificial turbulent fluctuations generated using- 29 -

this method satisfy both the prescribed profiles for the turbulent

Reynolds stress tensors as a 3  3 matrix. Here, it is extended tofluxes of momentum and the considered scalar quantities. The

account for the turbulent fluxes of a scalar within a 4  4 matrix.

authors then confirmed the flow field and dispersion of a passive

To impose time and space correlations for each component of

scalar in a boundary layer were well reproduced with artificiallythe fluctuation, a two-dimensional digital-filter method proposed

generated inflow boundary conditions, including wind velocity

by Xie and Castro(10) is employed. The prescribed time andand passive scalar fluctuations. Nevertheless, the performances

space correlations are assumed using exponential functions withof artificially generated fluctuations of wind velocity and air

an integral time scale, T, and a length scale, L: j t  j t   temperature in a non-isothermal flow field have not been  
 exp   ,
 2T 

(5)  
 exp   .
 2L 

(6) j t  j t 

investigated.
The objective of this study was to therefore apply artificiallyandgenerated turbulent fluctuations of wind velocity and air

 j r  j r   temperature as inflow boundary conditions for a large-eddy j r  j r 

simulation of a non-isothermal flow field around an isolated

The update equation for the artificially generated fluctuations on

high-rise building. The LES results were compared with those ofa grid point (m, n) is expressed asa wind tunnel experiment.

 t  j t  t , m, n    j t , m, n  exp  
 2T 

2. Overview of inflow turbulence generation and1 t  2
  j t  t , m, n1  exp   ,
 T 
scalar fluctuation
where

We developed an artificial generation method of the fluctuations

NyNz j t , m, n  of wind velocity and scalar based on Cholesky decomposition ofb b rm N y n N z(8)

the Reynolds stress tensor and turbulent scalar flux . In thism n m m,nn

.

(8)Here, r is a random number satisfying rj = 0 and ri rj = ij, Ny

study, the proposed method was employed to generate the wind

and Nz relate to the length scale of the filter in each direction,

velocity and air temperature fluctuations for LES inflow

and bk is a digital-filter coefficient for the integral length scale inboundary conditions.

the generated plane in each direction, as described below:

The proposed method is summarized as follows. The wind

1

 N
 2
~ 
~ 
bk  bk / 
bj2  ,
 j N

velocities and scalars are expressed as fi, the time-averaged
values of fi are denoted as fi, and the deviations from thetime-averaged value are expressed as fi'
f i  f i  f i ,

(7)



(9)

 x k ~.
bk  exp  

L 


(10)

where
(1)

where the subscripts, i = 1, 2, 3, represent the wind velocity
components in the streamwise, lateral, and vertical directions (u,By substituting j (as obtained using the new dataset of randomv, w), respectively, and i = 4 represents the scalar value, . A

numbers from Eqs. (8) to (10) into Eq. (7) for each time step, jregular matrix of the turbulent fluxes of the momentum andfor the next time step is obtained. Then, the fluctuations, fi, arescalar, Rij, is defined asgiven by substituting j into Eq. (4).

 u u  v u 
Rij  w u   u u v v v 
w v  v 

u w
v w
w w  w

u   v    .
w      

3. LES of non-isothermal flow around a building(2)

3.1. Calculation conditions

The Cholesky decomposition of Rij produces a lower triangular

For LES validation data, we selected a wind tunnel experimentmatrix, aij:

for a non-isothermal flow field around an isolated high-rise

 R11

 R21
a11

aij   R31

a11
R
 41
a11


0

0

2
R22  a210

R32  a21a31
R42  a21a41

a22
a22

building model with dimensions DD2D (H = 2D) in the02
2
R33  a31 a32

R43  a31a41  a32 a42





.
0


2
2
2
R44  a41
 a42
 a43


0

a33

streamwise (x), lateral (y), and vertical (z) direction, respectively.

(3)

It was conducted in a thermally stratified wind tunnel at Tokyo
Polytechnic University(11). The stratification of the approaching

With aij and a variable j that satisfies j = 0 and ij = ij,

flow was weakly unstable. The total measurement points were

the fluctuations, fi, can be rewritten as
f i  f i  f i  f i  aij  j .

260. Of these points, 152 were in the vertical center section, y =0.0 H. Another 108 points were located on the horizontal plane

(4)
(9)

at z = 1/4 H.

The transformation originally proposed by Le and Moin uses

- 30 -

The bulk Richardson number, Ri, is 0.10 and is defined as Ri =the first-order upwind scheme was automatically mixed to 40%(g /0)  (/ ) / (U / ) . Here,  is the boundary layer height,

when the convection boundedness criterion was violated toand U and  denote the mean streamwise wind velocity and air

avoid the appearance of unphysical oscillations in the result.

temperature at the boundary layer height, respectively. fThis scheme was implemented as “filteredLinear2V” in

represents the mean temperature of the wind tunnel floor, isOpenFOAM(12).

2Fluctuations of windvelocity andair

the temperature difference between and f:f , 0temperature for the LES inflow boundary condition were

denotes the mean air temperature in the boundary layer, and g is

artificially generated at intervals of 1/16 H on the inflow

the gravity acceleration.

boundary (7.5 H (y) × 6.25 H(z)) based on the method described

Fig. 1 illustrates the computational domain. The computationalin Section 2. The characteristic turbulent length, L, in eachdomain was set to 12.5H(x)6.0H(y)5.0H(z), where H is thecomponent of wind velocity was assumed to be a constant value,building height. The cross section of the computational domain

0.15 where  is the boundary layer height. The characteristic

corresponds to that of the wind tunnel. An orthogonal grid

time scale, T, was calculated from L and the mean streamwise

system was applied. The calculation conditions are summarized

wind velocity at the boundary layer height U based on Taylor’s

(12)

.

hypothesis of frozen turbulence, T = L / U. The turbulent length

The standard Smagorinsky model was used for the sub-gridand time scales of the air temperature were assumed to be same

scale turbulence model. The Smagorinsky constant was set tovalues of these of wind velocity in this study.

0.12. The building width and height were uniformly discretized

A run-up calculation was performed over t* = 200 with

into 20 and 40 with computational grids, respectively. The grid

dimensionless time, t*, which is defined as t* = t  uH / H.

size was expanded at a distance from the building with an

Here, t is the time calculation, and uH denotes the mean windexpansion ratio of 1.08. A second-order linear-interpolation

velocity at the building height of the inflow boundary. From that

scheme was used for the diffusion terms. In convection terms,

point, the turbulent statistics were collected and averaged over t*

in Table 1. LES was performed with OpenFOAM ver.2.4.0= 400 as the main simulation.
3.2. Results and discussion
The wind velocity and turbulent statistics were normalized withthe mean wind velocity at the building height of the inflow
boundary uH and building height H. The air temperature was

Figure 1: Computational domain
Table 1: Calculated conditions

Outlet boundary
Upper and sideboundaries

Ground andbuildingboundaries

Reynolds number

generatedtarget0147(x)  114(y)  90(z)Standard Smagorinsky model (Cs = 0.12)
First-order backward
Artificially generated fluctuations based on
the method proposed in Section 2
Advective outflow condition
Flow field: Slip wall
Temperature field: Adiabatic condition
Flow field: a two-layer model based on thelinear-power law expression proposed by
Werner and Wengle(13) was used toestimate wall shear stress.
Temperature field: The wall temperatures
were fixed with the experimental result. A
two-layer model based on thelinear-power law expression proposed by
Tuchiya et al. (14) was used to estimate
heat flux on the walls.

0.5
1u / uH []1.5(1) u / uHgeneratedtarget-1-0.8
(θ  θf ) /  []-0.6(2) (  f ) / 

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

generatedtarget

0

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
-1.2

z / H []Inlet boundary

12.5 H(x)  7.25 H(y)  6.5 H(z)z / H []Computational
domainGrid arrangement
SGS model
Time advancement

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

z / H []z / H []normalized with the temperature difference between the mean

0.01
0.02
u' 2 / uH2 [](3) u' 2 / uH2

0.03

2
1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0

generatedtarget

0

0.001
0.002
θ' 2 / θ2 []

0.003(4) ' 2 / 2

Figure 2: Comparison of targeted values and artificially

1.5  104 (uH  H )generated values for inflow boundary condition- 31 -

from the wind tunnel experiment shown in Fig. 3. The mean
wind velocity and air temperature predicted from LES agreed
very well with those of the experiment. The degree of agreement
between the experimental and LES results was evaluated by the
hit rate(15) and a factor of two of observations (FAC2) (16).
Evaluated values of these two validation metrics are comparedin Table 2. The hit rate, q, is used as a metric by the German
VDI Guideline(15) for the comparison of velocities, and it iscalculated fromq(1) Scalar wind velocity

1 N
 ni
N i 1
Pi  OiDOi,Pi  Oi  W

1 for

with

ni 

1 for
0

(11)

else

while FAC2 is evaluated using
N

1
FAC2 ＝ ∑ ni
N
i=1

Pi
2
Oi
.
Oi  W and Pi  W
0
else
0.5 

1 for

with

ni 

1 for

(2) Air temperature
Figure 3: Distributions of instantaneous scalar wind velocity(12)Here, N is the total number of measurement points, O is theand air temperature at x = 2.5 H (inflow boundary),observation result, and P is the predicted result. In this study, D

y = 0 H (center section) and z = 0.25 H

and W were set to 0.25 and 0.05, respectively.
Table 2 shows the metrics evaluated for the streamwise

temperature of the floor and that of the building height at the
inflow boundary (f − H).

component of the mean wind velocity and turbulent kinetic

Fig. 2 shows a comparison between the targeted values and theq  0.66 could be used as the quality acceptance criterion. Theenergy (TKE). According to the VDI guideline (15), the limit

artificially generated values for the inflow boundary condition

metrics are high for the streamwise component of the mean windbased on the method in Section 2. The generated mean windvelocity. However, the TKE has relatively low metrics.

velocity and air temperature completely agree with the targetedvalues. The variance of the fluctuation of streamwise component

1.5

that the generated turbulent fluctuations, as obtained by this
proposed method, satisfy the prescribed time-averaged values.
Fig. 4 illustrates the distributions of scalar wind velocity and airtemperature at x = 2.5 H (inflow boundary), y = 0 H (center-0.4
y = 1.0012x
R² = 0.9667

1.0Normalized Θ (LES)

in good agreement with the targeted values. Thus, we confirmed

Normalized u (LES)of wind velocity and the fluctuation of air temperature are also

0.5
0.0
-0.5

velocity and air temperature were observed at the planes y = 0 H

0.0
0.5
1.0
Normalized u (Exp.)

1.5(1) Streamwise component ofand z = 0.25 H. The fluctuations generated on the inflow

-0.8
-1.0
-1.2

-0.5

section) and z = 0.25 H. The spatial heterogeneities of windy = 0.9829x
R² = 0.9277

-0.6

-1.2

-1.0
-0.8
-0.6
Normalized Θ (Exp.)-0.4(2) Mean temperaturemean wind velocity

boundary were advected to the downstream region. On the

Figure 4: Comparison of experimental results with LES results

leeward areas of strong wind regions on the building side, theflow near the floor was elevated. Then, the high temperature air

Table 2: Validation metrics with q and FAC2

was transported to the upper region and the high temperatureareas were scattered on the leeward of the strong wind region.
The streamwise component of the mean wind velocity and mean
air temperature obtained from LES were compared with those- 32 -q

FAC2

Streamwise component ofmean wind velocity

0.87

0.98

TKE

0.44

0.77

fact that the artificially generated velocity fluctuations do not

Furthermore, Fig. 5 shows the distributions of the streamwise
component of the mean wind velocity, mean air temperature,satisfy the continuity and momentum equations, as shown by

TKE, and the variance of temperature fluctuation in vertical

Xie and Castro(10). The variance of air temperature fluctuationlines at x = −1.0 H, 0.0 H, 1.0 H, and 2.0 H, and y = 0. The mean

agrees well with the experiment, except for behind the inflow

wind velocity and air temperature predicted from LES agree

boundary. The distribution of the variance of air temperaturevery well with those of the experiment. TKE corresponds

fluctuation shows a similar trend as the distribution of TKE

relatively well with the experiment, except for just behind thebehind the inflow boundary. The variance rapidly decreased,

inflow boundary. TKE decreased by 40% of the imposed value

similar to the decrease of TKE.

as the inflow boundary condition immediately behind the inflow
4. Conclusions

boundary.
The cause of the above decrease of TKE likely relates to the
0

2.0

Exp.

u / uH [-]In this study, a time series of wind velocity and air temperature

LES

2

fluctuationsz / H [-]

generated

based

on

the

Choleskydecomposition of the Reynolds stress tensor and turbulent heat

1.5

flux. The artificially generated fluctuations were used as the1

inflow boundary condition of LES of flow around a building
model under a weakly unstable condition.

0.5

The mean wind velocity and air temperature predicted by

0
-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

LES using the generated inflow turbulence were generally inx / H [-]

good agreement with the experimental results.

(1) Streamwise component of mean wind velocity
-1.0

-0.5

2

Exp.

(f /  [-]Although TKE decreased by 40% of the imposed value as

LES

the inflow boundary condition immediately behind the inflow
boundary, TKE and the variance of air temperature fluctuation

1.5

z / H [-]

wasagreed reasonably well with the experiment, except for the1

region immediately behind the inflow boundary.
0.5

Acknowledgement

0

-1.5

-1

-0.5

0

0.5
x / H [-]

1

1.5

2

2.5

The authors are grateful for the support provided by the(2) Mean air temperature
0

0.1

Exp.

k / uH2 [-]Grant-in-Aid for Scientific Research (C) (Grant No. 26420578).

LES

2

Referencesz / H [-]

1.5

1(1) Kong H., Choi H., Lee J.S., Direct numerical simulation of
turbulent thermal boundary layers. Phys. Fluids, 12 (2000).

0.5(2) Lund T.S., Wu X., Squires K.D., Generation of Turbulent Inflow

0

-1.5

-1

-0.5

0

0.5

1

1.5

2

2.5

Data for Spatially-Developing Boundary Layer Simulations. J.

x / H [-]

Comput. Phys., 258 (1998), pp. 233-258.

(3) Turbulent kinetic energy
0

0.005

Exp.

'2 / 2 [-](3) Hattori H., Houra T., Nagano Y., Direct numerical simulation ofLES

2stable and unstable turbulent thermal boundary layers. Int. J. Heat
Fluid Flow, 28 (2007), pp.1262-1271.

z / H [-]

1.5

(4) Jiang G., Yoshie R., Shirasawa T., Jin X., Inflow turbulence1generation for large eddy simulation in non-isothermal boundary
layers. J. Wind Eng. Ind. Aerodyn., 104 (2012), pp.369-378.

0.5(5) Kataoka H., Mizuno M., Numerical flow computation around
0

-1.5

-1

-0.5

0

0.5

1

1.5

2

aeroelastic 3D square cylinder using inflow turbulence. Wind

2.5

x / H [-]

Struct., 5 (2002), pp.379-392.

(4) Variance of air temperature fluctuation(6) Tamura T., Nozu T., Okuda Y., Ohashi M., Umakawa H. Hybrid

Figure 5: Distributions of statistics in vertical lines

method of meteorological and LES models for heat environment.

at x = 1.0 H, 0.0 H, 1.0 H and 2.0 H

Proc. 8th Int. Conf. Urban Clim., (2012), Dublin, Ireland.

- 33 -

(7) Matsuda K, Onishi R, Takahashi K. High-resolution large-eddy(12) OpenFOAM User Guide (2014), http://openfoam.org/.

simulations for analyzing greenery effect in Tokyo Bay Zone. 8th(13) Werner H., Wengle H.. Large-eddy simulation of turbulent flowproc. Japanese-German Meet. Urban Climatol., (2017), pp. 69-72,

over and around a cube in plate channel. Proc. 8th Symp. Turbul.
Shear Flows, (1991), pp.155-158.

Osaka Japan.
(8) Okaze T., Mochida A. Generation of artificial inflow turbulence(14) Tsuchiya N., Iizuka S., Ooka R., Murakami S., Kato S. Analysis ofincluding scalar fluctuation for LES based on Cholesky

turbulence structure of a non-isothermal room airflow throughdecomposition. Proc. 9th Int. Conf. urban Clim., (2015), Toulouse,budgets of Reynold stresses and turbulent heat fluxes based on LES

France.

database. J. Archit. Plann. Environ. Eng., AIJ, 550 (2001),
pp.47-54.

(9) Le H, Moin P. Direct numerical simulation of turbulent flow over a
back-facing step. Rep TF-58 (1994), Thermosciences Div., Stanford(15) VDI. Environmental meteorology - Prognostic microscale wind

Univ.

field models evaluation for flow around buildings and obstacles,(10) Xie Z.T., Castro I.P. Efficient Generation of Inflow Conditions for

VDI Guideline 3783, Part 9, (2005), pp.53.

Large Eddy Simulation of Street-Scale Flows. Flow, Turbul.

(16) Schatzmann M., Olesen H., Franke J.. COST 732 model evaluation

Combust., 81 (2008),pp. 449-470.

case studies: approach and results, (2010), COST Action.

(11) Yoshie R., Jiang G., Shirasawa T., Chung J., CFD simulations of(Received Apr. 29, 2017, Accepted Jul. 7, 2017)

gas dispersion around high-rise building in non-isothermalboundary layer. J. Wind Eng. Ind. Aerodyn., 99 (2011), pp.279-288.

- 34 -