Computational Fluid Dynamic (CFD) Simulation of Flow around Tall Buildings

: Significant improvements of computer resources in recent past years allow to use computational fluid dynamics (CFD) methods as an alternative test method for wind tunnel tests in various wind engineering aspects. However the accuracy of CFD simulation mainly depends on careful setup of three main components of a CFD simulation, which are domain size, adopted mesh, and boundary conditions. This paper presents basic theoretical background of use of CFD for atmospheric boundary layer simulations and proper methods recommended for creating domains and meshes for CFD models. It also demonstrates several empirical methods that can be used as boundary conditions in the absence of more accurate data for simulation. The CFD simulation results of pressure distribution of 112m tall buildings is compared with the wind tunnel test results and found that performance of those empirical methods is satisfactory. The use of CFD simulation for flow visualization around a tall building and evaluating pedestrian level wind velocities are also demonstrated in this paper.


Introduction
Computational Fluid Dynamics (CFD) simulations are being widely used by engineers for various wind engineering studies such as determine wind loads on buildings, evaluating wind flow patterns in built areas, predicting pollutants depression patterns in urban areas, evaluate pedestrian level wind comfort, etc. Numerical simulations are more flexible and robust in terms of simulating wind flow conditions together with detailed surrounding features such as buildings, mountains, trees, etc. CFD technique is more economical and widely available compared to wind tunnel facilities. Rapid evolution of computer facilities enables to solve complex flow situations by using vigorous mathematical equations. However, the lack of knowledge in underlying theories and embedded features of CFD packages create this valuable tool as a black box for most of the users. According to the authors' point of the view wind engineering is not as much as developed compared to the other branches of civil engineering in Sri Lanka. Nowadays, the concern of wind engineering is increasing among Sri Lankan engineers and policy makers due to increase of damages due to frequent occurrence of high wind events and construction of many tall buildings in city centres, which are more susceptible to wind loads. Currently, there is no wind tunnel facility available in Sri Lanka to conduct detailed wind analysis. The only available method is use of wind loading standards, which are in fact, adopted from other countries with some modifications to comply with prevailing conditions in Sri Lanka. Therefore, CFD simulations would be a smart tool for practicing engineers to analyse wind conditions in vicinity of a tall building. According to the author's knowledge there is no any published paper in Sri Lanka on this area with a critical review. Thus, this paper is to cover basic knowledge of CFD simulation for simple modelling case, such as flow around a tall building within atmospheric boundary layer. However, the author would like to remind that, this paper only presents a brief explanation of some basic features of CFD simulation and for more details direct to read vigorous and critical studies [1,2,3].
The arrangement of this paper is as follows. In the section 2, it presents the basic concepts in atmospheric boundary layer flow and CFD calculations with necessary equations. The procedure of creating a CFD simulation with main components such as the domain, meshing, and assigning boundary conditions are discussed in section 3. Section 4 presents the results of CFD simulations and those results are compared with wind tunnel test results. Finally, conclusion and recommendations are given in section 5 for more detailed and precise simulations.

2.
Atmospheric boundary layer and CFD simulations 2.1. Governing equations for fluid flow within atmospheric boundary layer Wind engineers study more about the lower part of the atmosphere though entire earth atmosphere extends few kilometres above the earth surface. This lower part of the atmosphere is called as atmospheric boundary layer (ABL), which is directly under influence of earth surface itself such as shape, friction, thermal with time scales of less than a day and turbulent motion length scales of the order of boundary layer depth [4]. This boundary layer depth can be varying from several hundred meters to more than a kilometre aloft. Thus most of manmade structures are well within the atmospheric boundary layer, governing flow equations can apply in this layer easily. However, both time and length scales of atmospheric flows have large variations. Thus, numerical simulations are divided in to microscale and meso-scale based on time and length scale for easiness of study. In this study, term CFD simulation is strictly used to describe the simulation of smaller length (~10 cm -100 m ) and time scales (~1 minute -1 hour), thus probably in the category of micro-scale modelling.
Most of governing equations in fluid dynamics can be applied to the atmospheric flows. The main governing equations are about conservation of mass (Eq 1) and momentum (Eq 2). The latter is also known as Navier-Stokes equation for motion of the fluid Where, ui, uj are velocity components, ρ is air density, P is air pressure, µ is dynamic viscosity, and t is time.
In this equation Coriolis force and buoyance force are not considered as their effect is negligible in smaller length and time scale, which is valid for micro-scale CFD simulations.

Reynolds Average Navier -Stokes equation (RANS Equation)
The original Navier -Stokes equation was developed for the laminar flow but atmospheric flow is turbulent with three dimensional random unsteady motions. To combine the turbulence effect it is suggested a statistical approach where the value of velocity, uj , at any time can be split into a mean component and a fluctuating value, with Uj as the mean velocity, and uj' representing the turbulent variation.
The extra source terms on the right-hand side of the momentum equation, known as Reynolds number, defined as Where, τij is the turbulent stress tensor.
The rearranged equation is also known as Reynolds average Navier-Stokes (RANS) equation and it cannot be solvable directly. Therefore, turbulence models were introduced for the solution of the flow problems.

2.3.
Standard k -ε model The widely used 'industry standard' version of the turbulence model is standard k -ε model developed by Hanjalic and Launder [5]. This two equation model characteristics small scale turbulence by using two numbers kinetic energy (k) , and the energy dissipation rate (ε). This helps to calculate the turbulence transport as well as an empirical length scale. The kinetic energy per unit mass (k) is given by where, u, v, and w are wind speed components along, lateral and vertical directions.
Where k and ε are subject of the transport equations The model cannot be integrated through the near wall region as a singularity occurs at the wall surfaces. Thus wall functions must be introduced [6]. Verteesg and Malalasekara [7] highlighted that standard k -ε model is really successful in flow where the normal Reynolds stresses are less important. Even though this condition is rare in wind engineering, still researchers use this model due to its low computational cost, high numerical stability, and availability of verification data for wide variety of flows.

Main components of a CFD simulation
A CFD simulation includes three basic components such as; 1. Computational domain 2. Meshing 3. Boundary conditions for inlet, outlet, and walls.
Above components are either modelled as an objective or as a numerical function. The degree of accuracy of these components would directly affect the accuracy of results of the simulation.

Computational domain
The size of computational domain depends on the region that shall be represented by the simulation [8]. Domain should be large enough to avoid reflecting of fluid streams, which may cause abnormal pressure fields around the model. The blockage ratio should be below 3% as suggested by Baetke and Werner [9] is also a consideration for select size of the computation domain. For the single-building model, lateral and top boundaries should be set 5H or more away from the building, where H is the height of the target building [10].
The distance between the inlet boundary and the building should be set to correspond to the upwind area with proper roughness height. The outflow boundary should be set at least 10H behind the building. Where the building surroundings are included, the height of the computational domain should be set to correspond to the boundary layer height determined by the terrain category of the surroundings [11]. The lateral size of the computational domain should extend about 5H from outer edges of the target building. However, recommendations of Franke [8] about lateral boundaries are 2 to 3 times W, (where W is the width of built area) and the outflow boundaries 15Hmax, (where Hmax is the height of the tallest building) is a conservative approach. It should be noted that there is a possibility of unrealistic results, if the computational region is expanded without proper representation of surrounding of the building [12]. However, it is also true that unnecessary large fluid domain may demand higher computational cost. Therefore, researchers are being used arbitrary domain sizes, which are large enough to give satisfactory level of accurate results with reasonable computational cost. Figure 1

Eq (8)
Eq (9) Eq (10) Eq (11) Eq (12) Ground ENGINEER 47 ENGINEER 5 In Table 1, yp + is the dimensionless distance from the wall to centre of the first mesh element can be calculated as Equation (13). It is used to check the location of the first node away from a wall.
Where ρ is air density, Δyp is the distance from the wall to the first node, µ is the dynamic viscosity and uτ is frictional velocity at the wall, which equals to Where τw is the wall shear stress. This is more important when standard wall functions were used, yp + values should not be smaller than approximately 20 [8]. With the scalable wall functions and the automatic wall treatment, these values are only provided for information on the near wall resolution. In order to have a better simulation of flow at walls, these two criterions are followed, (1) minimum distance between nodes in the boundary layer; (2) minimum number of nodes within boundary layer. The goal is to determine the required near wall mesh spacing, (Δy), in terms of Reynolds number (Re in Equation (12) is the grid Reynolds number), length scale, and a suitable Δyp + target value. This target value depends on flow type and turbulence model used in, or in other words the near wall treatment in use. As a general guideline, a boundary layer should be resolved with at least:

Nnormal,min
Where Nnormal, min is the minimum number of nodes that should be placed in the boundary layer in the direction normal to the wall. The constructed mesh for all simulation by using meshing tool ICEM CFD® is shown in Figure 3. This structured mesh has over 273000 hexahedral cells. Structured mesh is more effective in cubic type models and gives undistorted flow conditions easily.

3.3.
Boundary conditions Boundary conditions represent the influence of surroundings that have been cut off by the computational domain. As they determine to a large extent the solution inside the computational domain, their proper selection is very important in a CFD simulation [8]. The boundary conditions for inlet, outlet and outer walls should be provided. These boundary conditions may be of the values for velocity, pressure or mass flow rates etc. In turbulent flow computations, additional boundary conditions for turbulence parameters need to be specified at inlet and outlet locations [12]. This information can be supplied in the form of convenient derived quantities such as turbulent intensity, length scale, viscosity ratio, hydraulic diameter, etc. In this study, velocity inlet used as inlet boundary with wind and turbulence profiles, outlet is modelled as outflow boundary and side walls and top wall have symmetrical boundary conditions, and ground (bottom wall) has no-slip boundary condition with suitable roughness length to represent upstream terrain roughness characteristics. In CFD simulations, roughness length is expressed as equivalent sand roughness length Ks. Blocken et al. ( [12], [13]) derived the relationship between the equivalent sand-grain roughness height for the ABL, KS,ABL and the roughness length z0, as KS,ABL≈30z0. This relationship generally could not be fully satisfied in the numerical simulation practice for the calculated KS, ABL usually exceeding the height of the first grid layer, yp. Thus, a compromised formula, KS = min (30z0, 1/2yp), is adopted to determine the value of roughness height KS.

Wind profiles
Vertical wind profile within the atmospheric boundary layer extends from ground level to few hundred meters (occasionally more than a kilometre). The shape, speed and direction of the wind profile are changing with various factors such as terrain type, topography features, atmospheric stability, etc. Therefore, real wind profile is changing in both time and space. Thus, defining an actual wind profile is extremely costly, tedious and time consuming procedure. Therefore, wind engineers use simple empirical wind profiles to define wind speed variation with height at a location. Two popular empirical methods are power-law method and logarithmic wind profiles. Both methods are combining surface roughness and wind speed to encounter the influence of surface friction on the wind velocity. Logarithmic wind profile has a solid theoretical background compared to power law wind profile. The standard logarithmic wind profile, used to describe fluid flow over rough surfaces, and requires the use of two unknown scaling 10 for wall function 15 for low Re model normal,
The friction velocity (u*) is a scaling velocity of the surface shear stress and defined by the relationship which depends on nature of the surface and mean velocity value.
Power law wind profile describes the mean profile by a simple power function of height [16] as Equation (17)   10 10 Where, α is the power law exponent, which changes with terrain roughness. This exponent has a relation to roughness length as Where zref is the reference height Figure 04 shows targeted wind profiles in wind tunnel test as defined by log-law and power law together with used wind profile for wind tunnel test. All Wind speeds are normalized with reference wind speed at test building height (0.28 m above ground).
The 1.2 m is the gradient height of wind profile used for wind tunnel test which represents 480 m in height in field condition. The wind tunnel wind profile and the power law wind profile are more identical. It is reasonable as wind tunnel test conducted according to the AIJ guidelines which also in power-law format. The power law exponent is 0.27 for AIJ wind profile. However, the deviation of log-law wind profile from power law wind profile indicates that, the conversion between two methods by using α and zo is not much accurate. The corresponding values are 0.885 m and 0.007 m for u* and zo respectively.

Turbulence model
Modelling of equilibrium boundary layer is an important precondition for the numerical simulation of flows around buildings. The horizontal inhomogeneity of the simulated ABL will result in additional errors to numerical results, and sometimes the influences of the additional errors are also quite significant [8].
In recent years, many researches ( [8], [11], [17], [18]) have emphasized the requirements of modelling of equilibrium ABL for the numerical investigation of flow around buildings. Blocken et al. ([12], [13]) has highlighted that CFD simulation of a horizontally homogeneous atmospheric boundary layer flow was difficult and it was vital for the successful application of CFD in wind engineering studies. Many researchers have devoted to improve the horizontal inhomogeneity of ABL ( [12], [19], [20]). The problem of constructing equilibrium boundary layer was further investigated from the viewpoint of the turbulence model itself ( [18], [21]). New sets of inflow turbulence boundary conditions for modelling equilibrium ABL was proposed by researchers ([19], [22]) to comply with near wall treatment. Richards and Hoxey [19] suggested a set of inflow turbulence boundary conditions, which were based on assumptions such as the vertical velocity, is zero, the pressure is constant and the shear stress is constant, etc. In their model, the kinetic energy 'k' and energy dissipation rate 'ε' can be expressed as Equation (19) and (20) respectively.

Figure 4 -Normalized Wind profiles from power-law, logarithmic law and wind tunnel test
Where, z -height at where energy dissipation rate measure from the wall z0-roughness length Further investigation of the turbulence model by Yang et al [22] proposed modifications for original k-ε equation proposed by Richards and Hoxey [19] .These modified equations are used for this study together with logarithmic wind profile model. The Equations (21) and (22) show the method to determine new k and ε values respectively as proposed by Yang et al [22] 2 0 Where, C1 and C2 are constant, which are found by nonlinear curve fitting techniques.
Cμ is usually smaller than the standard value of 0.09 in the sub-layer of ABL, in which the turbulence level is high. In order to adapt the characteristics of high turbulent flow, the value of Cμ was changed to 0.04 based on the statistics data for these equations [22]. During a wind tunnel test kinetic energy cannot be measured directly. Thus, turbulence intensity as shown in Equation (23) is first calculated from a wind tunnel test. Then it can be used to derive turbulence kinetic energy (k) as shown in Equation (24). The recommendations of Architectural Institute of Japan (AIJ) [23] to calculate turbulent kinetic energy (k) and turbulent energy dissipation rate (ε) are as shown in Equation (25) and (26).
Where σi is fluctuating component of velocity

Where,
Where, zG is gradient height of the atmospheric boundary layer, zs is the reference height, and α is power law exponent.

Setting up of turbulence models
In this study there are three models were created with three different wind profiles and k -ε turbulence models. One wind profile is based on wind tunnel test data and other two wind profiles are based on power-law method and log-law method respective. Those input boundary conditions are inserted to the commercial CFD software FLUENT® as user defined functions (udf).

4.1.
Turbulence models Turbulence kinetic energy (k) and energy dissipation rates (ε) used for three models are shown in Figure 5 and Figure 6 respectively. In wind tunnel test, kinetic energy was calculated by multiplying turbulent intensity data and wind speed data. However, energy dissipation was calculated by using AIJ guidelines as shown in Equation (26). AIJ method uses empirical formulae to calculate turbulence intensity (I) and turbulent kinetic energy (k). Energy dissipation rate was calculated according to Equation (26). These AIJ profiles were used together with power-law wind profile. For third model, logarithmic wind profile used with Equations (22) and (23) to calculate turbulent kinetic energy and energy dissipation rate respectively.
While wind tunnel kinetic energy profile has a complex shape, both empirical models have simple exponential curve shapes due to their simplifications. However, it is difficult to decide the best method to use for future simulations as both methods have their own advantages and disadvantages. Wind tunnel Eq (21) Eq (22) Eq (23)

Eq (25)
Eq (26) captures actual turbulence characteristics at different height but questionable whether it complies with the applied near wall treatment. On the other hand simplified empirical models cannot represent sharp turbulence variations at different heights but they more comply with the CFD properties. However, this kind of a difference might be reflected on CFD simulation results in a larger scale.

Figure 6 -Energy dissipation rate for different wind profiles
There is minor variation in empirical energy dissipation profiles from wind tunnel profile as shown in Figure 6. Neither power-law (AIJ method) nor log-law (Yang's method) could exactly replicate the wind tunnel energy dissipation profile accurately. However, two empirical profiles are same except at lower heights. This differences lead to some variations in final results of CFD simulations from wind tunnel test results.

Mean pressure coefficient
The Cp values are compared to determine the degree of successful of simulations. Cp is a normalized pressure coefficient defined as below Pr    [25]). Maximum negative value is -0.5 at the top of the building and it is -0.3 from middle of the building to the base. The maximum normalised pressure coefficient on windward side lies within a range of 0.85-0.89 for the CFD simulation of actual wind tunnel data. Two empirical k-ε models yield higher normalised pressure coefficient values compared to both wind tunnel test and CFD simulation of wind tunnel test data. This may be primarily due to inaccuracies in used k -ε models. However, these differences are not significant in magnitude and simulation results are acceptable.
Another indirect use of pressure distribution patterns observed from CFD simulation is designing wind tunnel test. With the knowledge of pressure distribution on sides of a building faces it is easy to decide pressure tap distribution on a building model. For an example, wind ward side needs more pressure taps at higher level due to the sharp pressure gradient compared to less number of pressure taps need for leeward side of the building.

Flow around a tall building
Flow simulation is another advantage of CFD simulation. Figures 10 to 12 show the stream patterns at ground level, mid-height level and top level of the building. It is clearly see that stream line deviation at front corners of the building and wake structure at leeward side of the building. However the strength of the wake is reduce at higher level as inflow wind speed is increasing. An important consideration is under estimation of wake size, a common short come of standard k -ε model highlighted by many researchers.
The stream line pattern as see with side view of the building shows the downwash and deviated flow over the building clearly in Figure 13. The contour map of wind speed at pedestrian level (2 m height above ground level) can be used to determine pedestrian level wind comfort. Figure 14 shows high wind speed areas extending from front corners of the building to sideways. Main wind shedding area is in leeward side of the building due to building wake. In front of the building there is a small wind speed up area in upstream from the building due to the down wash of the building. However, due to the under predicting nature of flow separation and wake structure of standard k -ε model, it is advisable to use this figure as a qualitative measure rather than as a quantitative measure.

Conclusion
Computational Fluid Dynamic (CFD) simulations are getting popular with the improvement of performance of computer resources. However, it is important to use this tool properly to obtain accurate results. The proper control over main components of a CFD simulation, fluid domain, meshing, and boundary condition would lead to an accurate simulation. Dimensions of fluid domain should be selected as its boundaries are not inducing false flow conditions on the model. Correct size of cells should create to have better performance of adopted wall function. Finer or coarser mesh selection depends on the degree of detail need from model area. In the absence of field data and/or wind tunnel data, modeller can use empirical boundary conditions for simulations. These empirical methods give satisfactory results, especially for simple static wind analysis, compare with wind tunnel test. Standard k -ε model has some advantages over other turbulence model such as simplicity, computational cost and numerical stability, but tend to under estimate flow separation, wake structure, etc. flow visualization is another advantage of CFD simulation, which can be used to understand flow around a tall building and determine pedestrian level wind comfort. However, it is recommended to use improved two equations type turbulence model such as realizable k -ε model or even