Keywords: arrays simulation

Finally, I am

Titre : Analyse par simulation des grandes échelles de l'écoulement de couche limite au-dessus d'une canopée urbaine Mots clés : Canopée urbaine, Couche limite, Simulation des grandes échelles, Bilan d'énergie cinétique turbulente, OpenFOAM Résumé : L'urbanisation croissante fait émerger des enjeux sociétaux et environnementaux relatifs à la pollution atmosphérique et au microclimat urbain. La compréhension des phénomènes physiques de transport de quantité de mouvement, de chaleur et de masse entre la canopée urbaine et la couche limite atmosphérique est primordiale pour évaluer et anticiper les impacts négatifs de l'urbanisation. Les processus turbulents spécifiques à la couche limite urbaine sont étudiés par une approche de simulation des grandes échelles, dans une configuration urbaine représentée par un arrangement de cubes en quinconce. Le modèle de sous-maille de type Smagorinsky dynamique est implémenté pour mieux prendre en compte l'hétérogénéité de l'écoulement et les retours d'énergie des petites vers les grandes structures.

Le nombre de Reynolds basé sur la hauteur du domaine et la vitesse de l'écoulement libre est de 50000. L'écoulement est résolu dans les sous-couches visqueuses et le maillage est raffiné dans la canopée. Le domaine est composé de 28 millions de cellules. Les résultats sont comparés à la littérature et aux données récentes obtenues dans la soufflerie du LHEEA. Chaque contribution au bilan d'énergie cinétique turbulente est calculée directement en tout point. Cette information, rare dans la littérature, permet d'étudier les processus dans la sous couche rugueuse. Grâce à ces résultats 3D, l'organisation complexe de l'écoulement moyen (recirculations, vorticité, points singuliers) est analysée en relation avec la production de turbulence. Enfin, une simulation où les obstacles sont remplacés par une force de traînée équivalente est réalisée à des fins d'évaluation de cette approche.

Title : Analysis of the unsteady boundary-layer flow over urban-like canopy using large eddy simulation.

Keywords : Urban canopy, Boundary layer, Large-Eddy Simulation, Turbulent kinetic energy budget, OpenFOAM Abstract : The rapid development of urbanization raises social and environmental challenges related to air pollution and urban climate. Understanding the physical processes of momentum, heat, and mass exchanges between the urban canopy and the atmospheric boundary-layer is a key to assess, predict and prevent negative impacts of urbanization.

The turbulent processes occurring in the urban boundary-layer are investigated using computational fluid dynamics (CFD). The unsteady flow over an urban-like canopy modelled by a staggered arrangement of cubes is simulated using large eddy simulation (LES). Considering the high spatial and temporal inhomogeneity of the flow, a dynamic Smagorinsky subgrid-scale model is implemented in the code to allow energy backscatter from small to large scales.

The Reynolds number based on the domain height and free-stream velocity is 50000. The near-wall viscous sub-layers are resolved and the grid is refined in the canopy resulting in about 28 million grid cells. LES results are assessed by comparison with literature and data recently acquired in the wind tunnel of the LHEEA. The turbulent kinetic energy budget in which all contributions are independently computed is investigated. These rarely available data are used to analyse the turbulent processes in the urban canopy. By taking advantage of the three-dimensionality of the simulated flow, the complex 3D time-averaged organization of the flow (recirculation, vortices or singular points) is analyzed in relation with production of turbulence. Finally a drag approach where obstacles are replaced by an equivalent drag force is implemented in the same domain and results are compared to obstacle-resolved data.

"The two most important days in your life are the day you are born and the day you find out why." -Mark Twain. [START_REF] Castro | Are Urban-Canopy Velocity Profiles Exponential?[END_REF] The logarithmic profile for the mean velocity above the cube, before refining the mesh around the cube (White circle) and after refining the mesh around the cube (Red point), U + = z + is in Black line. . . . . . . . . . . . . . . . 2.14 Plan view of the refined mesh around the cubes. . . . . . . . . . . . . . . . 2. [START_REF] Castro | Measurements and Computations of Flow in an Urban Street System[END_REF] The contour plot of the two-point correlation of streamwise velocity Ruu, where the reference position is in the wake of the cube at a height of 1.5h. .

Introduction

With the constant development of urbanization, the ratio of population living in cities continues to grow. The United Nations said that 68% of the world's population projected to live in urban areas by 20501 . Followed from urban population growth, pollution problems in cities are more and more prominent: among them, air quality catches citizens special attention. Cities' micro-climates, including air quality, are the result of complex mechanical and thermal interactions between the atmospheric boundary layer (ABL) and the underlying buildings. Many different kinds of building layouts, heights and shapes exist in cities, such as the heterogeneous arrangement of high rise and low rise buildings in the downtown of Beijing (see Figure 1 The ground roughness plays a crucial role in the momentum, heat and mass transfer between the ABL and the urban canopy layer where populations live. From an aerodynamic point of view, the wind in the ABL is a high Reynolds number flow exhibiting high shears as well as high turbulence levels. Analyzing and understanding accurately the dynamics of the turbulent flow inside the urban boundary layer provides important information contributing to a better evaluation of the impact of buildings on the micro-scale climate. However, the high complexity of turbulent interactions at play remains an important scientific challenge. The present study focuses on the micro-scale analysis of the urban boundary layer at the scale of a district or streets. The urban canopy here is represented by a staggered arrangement of cubes of constant height. Such detailed studies, resolving explicitly the building geometry, are of high scientific interest to understand the dynamics of the urban canopy flow. This type of work also participates to the development of simplified models used at regional scale where buildings cannot reasonably be explicitly resolved. In this work, the performance of the drag-porosity approach where buildings are represented by an equivalent drag force is evaluated. The main steps of the work and organization of the manuscript are summarized in Figure 2.

Objectives and structures of the thesis

The first objective of this work is to set-up and run reliable numerical simulations of the flow field in an urban boundary layer. The bibliographic review proposed in Chapter 1 describes the up-to-date general knowledge on urban boundary layer and details some of the most salient work made on site, in a wind tunnel or in numerical simulations over the last decades on realistic and simplified urban geometries.

Several settings for numerical simulations using Large Eddy Simulation (LES) with OpenFOAM are evaluated in Chapter 2 in order to select the most appropriate parameters. This study is made on a [4h × 4h × 4h] domain for cost efficiency. A precise parametric study draws special attention to the choice of the mesh size and of the subgrid-scale model. Using the most appropriate numerical configuration found in Chapter 2, a simulation over a domain of [16h × 12h × 8h] using about 28 million cells is presented in Chapter 3. The assessment of the performance of the obstacle-resolved LES, the second objective of this work, is made throught a detailed comparison of the results against literature data. Mean flow and turbulent-related variables up to the third order momentum are carefully discussed in light of literature data from state-of-art wind tunnel experiments and numerical simulations. This comparison includes the latest experimental results obtained in the atmospheric wind tunnel of the LHEEA lab.

The third objective of the work is to contribute to the understanding of dominant turbulent processes at play in the urban canopy (Chapter 4). For that, all terms of the turbulent kinetic energy (TKE) budget are explicitly computed at every spatial locations of the domain. The spatial organization of the TKE budget terms are then depicted in 2D and 3D figures and carefully analyzed to extract the most salient features.

A drag-porosity approach is implemented in OpenFOAM (Chapter 5) with the objective to evaluate its performance compared to the simulation where obstacles are explicitly resolved. Its ability to reproduce the complex turbulent flow of an obstacle-resolved simulation is evaluated by comparing turbulent variables such as high order momentum and the TKE budget.

Finally, Chapter 6 presents conclusions and perspectives of the work.

Chapter 1

Atmospheric boundary layer: structure and characteristics

The scope of this chapter aims to introduce the structure of the atmospheric boundary layer and the morphological parameters used in the urban canopy flow study. The nature of the canopy flow and the method that is used to represent the urban canopy in numerical simulation are illustrated. Then, the statistical parameters that are used in the analysis procedure are described. Finally, the previous research background is introduced at the end of this chapter.

The structure of the atmospheric boundary layer

Atmospheric boundary layer (ABL) extends from the ground up to 1-2 kilometres, which is the primary living space for humans, animals, and plants. Unlike the overlying free atmosphere where the flow dynamics is insensible to the surface friction, the air inside the ABL is in direct contact with the surface and is significantly affected by the underlying roughness. Due to the complexity and diversity of the covers and roughness elements in the city areas, precisely defining the hierarchical structure of the atmospheric boundary layer will be very beneficial in understanding the critical turbulent mechanism of the boundary layer flow.

As quoted in [START_REF] Kaimal | Atmospheric boundary layer flows: their structure and measurement[END_REF], [START_REF] Sutton | Micrometeorology[END_REF] suggested dividing the ABL into two parts: the surface layer and the layer above. In the surface layer, which develops from the surface up to 50 -100m, the wind structure is controlled by surface friction and temperature gradients. The layer above the surface layer extends to a height of 500 -1000m.

The wind structure in this layer is affected by the surface friction, temperature gradients, and Earth's rotation.

With the in-depth understanding of the near-wall turbulence, researchers realized that the structure of the surface layer could be divided into the inertial sublayer (ISL) and the roughness sublayer (RSL). Figure 1.1 schematically depicts the atmospheric boundary layer structure including five different layers in the vertical direction. In boundary layer applications, identifying the atmospheric boundary layer structure is very helpful in demonstrating the critical features of each layer. For example, the momentum, turbulence, and mass exchange are very intense within the roughness sublayer, especially inside the urban canopy layer (UCL). It is worth mentioning that the comprehension of these vertical structural characteristics helps to understand the interaction between the atmosphere and underlying surface, as well as to identify the primary physical mechanisms of the flow processes that occur within the urban area. The following sections will briefly describe these layers separately.

Outer layer

The outer layer (also known as the mixed layer) is a layer in the upper part of the ABL that extends from the top of the surface layer up to the bottom of the free atmosphere. The potential temperature and specific humidity here are nearly constant with the height. The wind in this layer is the result of pressure gradient force and Coriolis force. The turbulence is generated by surface friction and heat transfer from the surface to the atmosphere. The thermal effect is very influential in the development of a well-mixed layer. By neglecting thermal effects, one assumes that the static stability of the atmosphere is neutral, and the turbulence is only generated by dynamic processes. In this case, this layer is called 'Ekman layer'. The thermal effects will not be considered or discussed in this thesis.

Inertial sublayer

The inertial sublayer (ISL) is the upper part of the sur f ace layer, wherein the flow characteristics are substantially homogeneous in the horizontal direction, and only depends on the height above the ground (z). In the neutral atmosphere, over a given period, the time-averaged streamwise velocity U can be characterized by a logarithmic profile, defined as:

U(z) = u τ κ ln( z -d z 0 ) (1.1)
In equation (1.1), u τ denotes the friction velocity (see section 1.4.3). The von Kármán constant κ can take the values ranging from 0.35 to 0.43, but κ = 0.4 is usually adopted in the boundary layer studies [START_REF] Kaimal | Atmospheric boundary layer flows: their structure and measurement[END_REF]. z 0 and d are aerodynamic parameters, which depend on the height and densities of the roughness elements. A zeroplane displacement height d is introduced in case of dense canopies to account for the effective height of the ground due to vertical flow displacement and drag force distribution through the canopy [START_REF] Raupach | Rough-wall turbulent boundary layers[END_REF]. Several parameterization methods have been proposed to determine the aerodynamic parameters d and z 0 as a function of morphological parameters such as canopy height, plan packing density, frontal area density (Macdonald et al. 1988, Kastern-Klein and[START_REF] Kastner-Klein | Mean Flow and Turbulence Characteristics in an Urban Roughness Sublayer[END_REF]). When the flow field is available (as in the present work), z 0 and d can be obtained from the mean velocity profile and friction velocity. Another way (used here) to calculate d in urban configuration is directly using the mean drag profile within the canopy, and assuming d is the mean height of the momentum absorption by the surface or the mean height at which the total drag acts [START_REF] Jackson | On the displacement height in the logarithmic velocity profile[END_REF][START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF]. It is then defined as:

d = h 0 zD(z) dz h 0 D(z) dz (1.2)
where D(z) is the drag profile within the canopy, which is calculated by the time-averaged pressure difference between the front and back face of an obstacle.

Roughness sublayer

The roughness sublayer (RSL) is the lower part of the surface layer located below the inertial sublayer, in which the flow characteristics depends explicitly on the properties of the roughness [START_REF] Raupach | Conditional statistics of Reynolds stress in rough-wall and smooth-wall turbulent boundary layers[END_REF]). The RSL extends from the ground surface up to a height of 2 ∼ 5 times the mean height of roughness elements (h). For instance, Rotach (1993a) observed that the level z/h = 1.55 is within the RSL from the experiment conducted in the center of Zürich. [START_REF] Roth | Turbulent transfer relationships over an urban surface. I. Spectral characteristics[END_REF] found evidence that the level z/h = 2.65 is still within the RSL over a suburban area. The airflow in the roughness sublayer above the urban canopy layer is very sensitive to the underlying canopy and atmospheric flow above [START_REF] Barlow | Progress in observing and modelling the urban boundary layer[END_REF]. Flow characteristics are significantly affected by obstacles and result in very complex three-dimensional flow structures.

Urban canopy layer

The urban canopy layer (UCL) is the bottom part of the RSL. The UCL extends from the ground up to the top of canopy [START_REF] Oke | The distinction between canopy and boundary-layer urban heat islands[END_REF][START_REF] Roth | Review of atmospheric turbulence over cities[END_REF]. The local climatic characteristics of urban areas are significantly affected by various factors, such as the frictional effects of rough surfaces, the exchange of heat and humidity between the surface and the atmosphere, and the anthropogenic heat sources. In addition, the presence of obstacles such as trees and buildings can strongly influence the mean flow and the characteristics of the turbulence. The behavior of the flow also depends on the shape, height and arrangement of the obstacles. These factors lead to a very complex flow dynamics within the urban canopy. For now, the understanding of the turbulence dynamics involved in the urban canopy is still under discussion and not fully understood even if recent studies have investigated the dynamic interactions between the large-scale structures of the boundary layer and the smaller structures associated with canopy flows [START_REF] Blackman | Non-linear interactions in a boundary layer developing over an array of cubes using stochastic estimation[END_REF][START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF].

Effect of the urban canopy roughness

The height and arrangement of obstacles in urban canopy play an important role in the characteristics and structure of turbulent flow, which significantly affects the transfer of momentum, heat and mass in the roughness sublayer. Therefore, it is necessary to have a detailed understanding about the type of urban canopy configuration and its impact on the flow.

Morphological parameters 1.2.1.1 Mean height of the buildings

The height of the building is an important parameter for canopy flow studies because it is closely related to flow patterns and frictional effects in urban areas. Xie et al. (2008) performed numerical simulations over a random height urban model and found that the relatively high pressure on the tallest building generates significant contributions to the total drag. However, the difference in the height of buildings throughout the city is very large, and there is no standard which can be used to fully understand the flow structure. Therefore, using the mean height of the buildings to simplify the problem can help to establish a unified research standard.

The mean height of the canopy (h) is defined as a weighted average of the buildings' height:

h = n ∑ i=1 h i × A pi n ∑ i=1 A pi (1.3)
where h i is the height of building i, n is the total number of buildings and, A pi is the plan area of the element i. More details about these parameters are shown in a schematic plot (Fig. 1.2).

Plan area density

The plan area density λ p is an important factor in classifying the flow and quantifying the impact of the canopy on the flow. It refers to the surface area of buildings in a particular urbanized area. λ p is defined as the ratio of the total plan area occupied by obstacles to the total plan area A Total (see Figure 1.2):

λ p = n ∑ i=1 A pi A Total
(1.4)

Frontal area density

The frontal area density λ f is the ratio of the total windward area of the buildings to total plan area (A Total ), and is defined as:

λ f = n ∑ i=1 A f i A Total (1.5)
where A f i is the cross-sectional area of the obstacle i normal to the incident flow (see Figure 1.2).

In the present work, the urban canopy will be represented by an array of cubic obstacles of constant height so that the canopy height and cube height are equal. The cubes being aligned with the flow, A pi and A f i have the same value and λ f is equal to λ p .

Volumic frontal density

The frontal density per unit volume α f (z) is a parameter which describes the frontal surface area per air volume. It is calculated according to the geometry of the obstacles present in the grid cell. When it is defined at different levels (z) of the canopy as in equation (1.6), it can be used in the drag-porosity approach (Sec. 1.3.3). It has a dimension of [m -1 ] and defined as:

α f (z) = n ∑ i=1 A f i (z) dxdydz - n ∑ i=1 V f i (z) (1.6)
where A f i (z) is the frontal area of the obstacle i normal to the incident flow at level z , and V f i (z) is the volume occupied by the obstacle i at level z (Fig. 1

.2).

Flow classification

The distinct architectures inside the urban area profoundly disturb the upcoming flow and make the study of flow characteristics within the boundary layer very difficult. To describe the flow regime inside the RSL, the ratio of the gap between obstacles (or street width) (W ) to the height of the buildings (H) is applied. [START_REF] Oke | Street design and urban canopy layer climate[END_REF] classified three typical flow patterns in 2D roughness geometry to state the canopy flow regimes. The roughness geometry was characterized by the aspect ratio H/W , which is the ratio of building height to building spacing. These regimes are referred as isolated roughness flow, wake interference flow, and skimming flow. The schematic plot of these flow regimes are shown in Figure 1.3. [START_REF] Perry | Rough wall turbulent boundary layers[END_REF] classified two types by adjusting the space between the neighboring cubes: d -type (W ≤ H) roughness and k -type (W ≥ 3 H) roughness. In k -type roughness, the friction effect plays an important role and is highly dependent on the friction Reynolds number. At the same time, a strong interaction between the flow and the obstacles occurs at the top of the roughness, causing the eddy shedding in the wake regions of the elements. In dtype roughness, the friction coefficient does not depend on the scale of the roughness, but on the outer scale such as the boundary layer height (δ ), and there is essentially little eddy shedding in the main flow. The behaviour in terms of turbulence between these two types are different [START_REF] Jimenez | TURBULENT FLOWS OVER ROUGH WALLS[END_REF][START_REF] Kanda | Large-Eddy Simulations on the Effects of Surface Geometry of Building Arrays on Turbulent Organized Structures[END_REF]). [START_REF] Leonardi | Properties of d-and k-type roughness in a turbulent channel flow[END_REF] demonstrated that the frictional drag is much larger than the pressure drag in dtype roughness, whereas the pressure drag dominates in ktype roughness.

The layout of the buildings

The layout of the buildings is not only a crucial factor in urban design but also in the urban turbulence surveys because the location and orientation of the turbulent structure are firmly related to the incoming ABL flow and local flow within the roughness sublayer. The relative position of the obstacles has a substantial influence on the direction and speed of the flow. Since complex architecture significantly limits the understanding of canopy flows, simplified cubic array models can be established to achieve the purpose of urban research. This section aims to describe three common types of the cubic array arrangements. As an example, we consider a plan area density λ p = 25% (Fig. 1.4).

In aligned arrays (Fig. 1.4 (a)), the cubes are placed in a single row in the wind direction, with a distance of one cube length between successive cubes in the streamwise direction and of three cubes length in the spanwise direction. In squared arrays (Fig. 1.4 (b)), the cubes are located along lines parallel and perpendicular to the wind direction with a distance of one cube length between cubes. In staggered arrays (Fig. 1.4 (c)), the distance between successive cubes in the the streamwise direction is three cubes length, and one cube length in the spanwise direction. It is interesting to distinguish the flow characteristics between these simplified cubic arrangements. [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF] found that the drag coefficient in the aligned and squared arrays are much smaller than in the staggered case. The consequence is a reduced spatialaveraged streamwise velocity within the staggered array compared to the other two types. The spatial-averaged spanwise and vertical turbulence intensities for the staggered array are found similar to the aligned array, but the streamwise turbulence intensity is much lower within the canopy due to the enhanced drag of the staggered cubic array.

Representation of the urban canopy in numerical simulations

In ABL models, there are different ways to represent the urban canopy, and the corresponding representation method can be selected in the model according to the size of the research domain. This section aims to introduce the common ways to take into account the urban canopy in computational fluid dynamics (CFD) applications.

Obstacle resolving method

The obstacle resolving method refers to a simulation in which the obstacles are explicitly described inside the computational domain. The benefit of this approach is that the presence of the obstacles produces a real fluid environment to simulate the flow features resulting from the interactions between the flow and the obstacles.

x/h z/h Although this method has apparent advantages, it also has drawbacks. For example, due to the existence of large velocity gradients near solid boundaries, high resolution is required in these regions. This results in a large number of cells near walls and makes prohibitive simulation cost (see Figure 1.5). These shortcomings prevent the use of obstacle resolving method combined with unsteady turbulence modeling for large-scale ABL studies, but it can be used in simplified academic configurations [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]Xie andCastro 2006). In the present study, obstacle resolving method will be used to study the turbulent flow over a staggered cube array.

Immersed boundary method

The application of obstacle resolving method with an arbitrary three-dimensional shape in space is still far from mature. The immersed boundary method (IBM) solves the problem of coupling between the complex obstacle shape and the flow field. The main benefit of this approach is that the force exerted by the obstacle on the flow is introduced in the momentum equation (equation 2.2), which allow to create complex and evolutive obstacle shapes. It was first used to represent a flexible and moving boundary exerting a force on the fluid [START_REF] Peskin | Flow patterns around heart valves: A numerical method[END_REF]. From the application perspective, there are two main ways: continuous forcing method and discrete forcing method. In the former, the force term is added to the continuous Navier-Stokes equation before discretization, whereas in the latter, the force is applied (explicitly or implicitly) to the discretized equations. The application of immersed boundary method to study the urban canopy has been performed by many researchers (Mittal and The results were fully validated by field measurements data under near-neutral conditions and confirmed the possibility of combining LES with this method to analyze the pollutant dispersion in urban areas.

Drag-porosity approach

Similar to the principle of the immersed boundary method, the drag-porosity method expresses the existence of obstacles and their influence on the turbulent flow through the action of an averaged drag force. This approach is often used in forest canopies [START_REF] Kanda | Organized structures in developing turbulent flow within and above a plant canopy, using a large eddy simulation[END_REF][START_REF] Su | Turbulent statistics of neutrally stratified flow within and above a sparse forest from large-eddy simulation and field observations[END_REF][START_REF] Dupont | Influence of foliar density profile on canopy flow: A large-eddy simulation study[END_REF][START_REF] Yue | Turbulent kinetic energy budgets in a model canopy: comparisons between LES and wind-tunnel experiments[END_REF]). The difference with IBM is that a set of obstacles can be considered as a porous medium. The configuration of the urban canopy is modeled by the profile of the drag coefficient C D (z) that is related to the drag force F D through the equation:

F Di = ρ 2 C D (z)α f (z)u i |u| (1.7)
where α f (z) is the volumic frontal density (Eq. 1.6), u i is the i th velocity component, and |u| = √ u j u j . This drag force constitutes a sink of momentum in the momentum conservation equations (See chapter 5, Eq. 5.3).

The advantage of the drag-porosity approach is to considerably reduce the number of computational cells, especially in the lower part of the simulation domain, thereby saving the calculation cost. This method has already shown a good agreement of the turbulent statistics in comparison with the literature [START_REF] Maché | Représentation multi-échelle des transferts entre couche de canopée urbaine et atmosphère à l'échelle de la ville[END_REF]). It can be used in applications at city scale or regional scale or to generate incoming flows in microscale applications. However, the local characteristics of the flow cannot be modeled using this method, so that some flow features which could be important in the flow dynamics and exchanges within the RSL can be wrongly predicted.

Aerodynamic roughness length

In the mesoscale atmospheric study, the flow can be considered as over an homogeneous surface since the scale is too large to explicitly describe the local effects of the roughness elements. Therefore, the aerodynamic parameter z 0 and d (see section 1.1.2) are prescribed to model the effect of the surface roughness on the wind field in the inertial sublayer. The vertical fluxes between canopy and atmosphere are usually estimated based on the Monin-Obukhov similarity theory. The typical values of aerodynamic roughness length z 0 (assuming d = 0) are summarized in Table 1

Statistical variables

In this section, the statistical variables and the scaling parameters (friction velocity) used in Chapter 3, Chapter 4 and Chapter 5 to validate the performance of the model and to make a quantitative analysis of the turbulent flow are described.

First-order statistics Time-averaged quantities

The time averaging is defined as the arithmetic mean of a series of values obtained over a given period. For instance, the time average of the variable ϕ is defined as:

ϕ = ∑ n k=1 ϕ k × ∆t k ∑ n k=1 ∆t k (1.8)
where the (..) symbol denotes a time-averaged value, ϕ k is the value of the variable ϕ at each time-step ∆t k .

Spatial-averaged quantities

The way to do the spatial averaging inside the canopy have been widely discussed recently [START_REF] Finnigan | Turbulence in plant canopies[END_REF][START_REF] Kono | Numerical Investigations of Mean Winds Within Canopies of Regularly Arrayed Cubical Buildings Under Neutral Stability Conditions[END_REF][START_REF] Böhm | Turbulence Structure Within and Above a Canopy of Bluff Elements[END_REF]and Castro 2017). The discussion is mainly about the quantities averaged over only the fluid area or over the entire domain including both fluid and obstacles.

In general, the spatial averaging of a flow property ϕ over a volume V can be expressed as:

⟨ϕ⟩ = 1 V V ϕ (x, y, z) dV (1.9)
where the ⟨..⟩ symbol denotes the spatial averaging value. When V corresponds to the entire domain including both fluid and obstacles, the spatial average is an extrinsic spatial average (ESA), defined as:

⟨ϕ⟩ ESA = 1 V f +V s V f +V s ϕ(x, y, z) dV = 1 V f +V s V f ϕ(x, y, z) dV (1.10)
where the V f is the fluid volume and V s is the solid volume.

When one considers only the fluid volume (V = V f ), then the spatial average is an intrinsic spatial average (ISA), defined as:

⟨ϕ⟩ ISA = 1 V f V f ϕ(x, y, z) dV (1.11)
There is no difference between these two averages above the canopy because V s = 0 there. The main difference occurs inside the canopy. For instance, for a cubic array with plan area density λ p = 0.25, the fluid volume V f is triple of the solid volume V s inside the canopy layer, and then:

⟨ϕ⟩ ESA = 3 4 × 1 V f V f ϕ(x, y, z) dV = 3 4 × ⟨ϕ⟩ ISA (1.12)
Kono et al. (2010) suggested that intrinsic spatial average could be reasonable when λ p is small because there is less effect by the obstacles. Xie and Fuka (2017) analyzed the differences between these two methods over aligned cuboid arrays and observed that extrinsic spatial average was better than intrinsic spatial average in avoiding the abrupt change of total shear stress at the canopy height. These two methods will be tested in chapter 5 in order to compare the drag approach to the obstacle resolving method.

Second-order statistics

Standard deviation

In statistics, the standard deviation (σ ) is a measure of the degree of the variation to the mean value of a dataset. It is defined as the arithmetic square root of the variance, which reflects the degree of dispersion between individuals within the dataset.

The standard deviation of turbulent variable (ϕ) is commonly used for statistical analysis of turbulence. The definition is as follows:

σ ϕ = 1 N N ∑ k=1 (ϕ k -ϕ k ) 2 = (ϕ ′ k ) 2 (1.13)
where ϕ k is the turbulent property of time-step k, ϕ k is time-averaged value and ϕ ′ k is the fluctuation with respect to ϕ k .

Velocity co-variances

Velocity co-variances are the mean statistical correlation between the components of the velocity fluctuations. It is defined as:

u ′ i u ′ j = 1 N N ∑ k=1 (u i -u i )(u j -u j ) (1.14)
where u i is the velocity components in the i th direction. The mean velocity u i determined by time averaging, and u ′ i is the fluctuation of the velocity.

In this study, velocity correlation is calculated using the following method:

u ′ i u ′ j = u i u j -u i u j (1.15)
The Reynolds shear stress (-u ′ w ′ ) is defined from the mean correlation of the streamwise and vertical components of the velocity fluctuations.

Friction velocity

Friction velocity is a reference velocity, which is widely used as a scaling parameter in the surface layer studies. It represents the effect of wind stress (τ) on the surface and varies with the nature of the roughness [START_REF] Kaimal | Atmospheric boundary layer flows: their structure and measurement[END_REF]). It can be derived from the average drag force (F D ) acting on the cubes through the relation:

u τ = τ ρ (1.16)
where τ = F D /A with A is total plan area (A Total ) and F D is the drag force. The contribution of the viscous force is neglected [START_REF] Cheng | Flow over cube arrays of different packing densities[END_REF]).

In some wind-tunnel experiments, the wind stress (τ) is not directly measured data. Therefore, assuming a constant flux layer, the friction velocity is estimated using the Reynolds shear stress data at the canopy top (e.g. [START_REF] Böhm | Turbulence Structure Within and Above a Canopy of Bluff Elements[END_REF][START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]). This scaling velocity is called here u * and defined as:

u * = -u ′ w ′ (1.17)

Skewness

Skewness (Sk) is a statistical parameter that evaluates the degree of asymmetry of a probability distribution. The value of the skewness can be a positive, negative or zero. The skewness for turbulent property ϕ is the third moment of fluctuation normalized using its standard deviation:

Sk ϕ = ϕ ′ 3 /(σ ϕ ) 3 (1.18)
The skewness of the velocity is computed to describe the high (or low) speed region, which links with the nonlinear process that occurs inside the canopy layer. It can provide insight into momentum transfer events such as sweep and ejection events [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]).

Two-point correlation coefficient

The two-point correlation of velocity is a crucial parameter in turbulent statistical studies, which states the relationship between the velocity fluctuations at two points in the flow field. The structure of turbulence and the spectral distribution of turbulent energy can be obtained through the analysis of two-point correlation. The two-point correlation coefficient of velocity is defined as:

R i j (x re f , y re f , z re f , x, y, z) = u ′ i (x re f , y re f , z re f ) × u ′ j (x, y, z) σ u i (x re f , y re f , z re f ) × σ u j (x, y, z) (1.19)
where (x re f , z re f ) is the coordinate of the reference location, and (x, z) is the coordinate of the random location.

The power spectral density characteristics

The power spectral density (E) shows the energy distribution with different frequencies (corresponding to different scales) of the turbulence. It is defined as :

E( f ) = 1 2π ∞ -∞ R(t)e -iωt dt (1.20)
where f is the frequency, t is time, R(t) is the autocorrelation of velocity signal.

Typically, the velocity fluctuations are used to reflect the distribution of turbulent kinetic energy in different sizes. As an example, Fig. 1.7 shows the temporal frequency spectrum obtained from the wind-tunnel experiment of [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF]. The most of the turbulence energy occurs at the low-frequency regions. In the range of 0.1 < kz ′ < 0.6, the curve has a slope of approximately -1, while in the range of 0.6 < kz ′ < 10, the curve has a slope of approximately -5/3. This slope implies the presence of an inertial subrange. Inertial subrange is between the large-scale energy-containing range at which energy is injected and the small-scale dissipation range at which energy is dissipated. Kolmogorov (Kolmogorov 1941a(Kolmogorov , 1941b) ) found that for high Reynolds numbers, the inertial subrange is universal and uniquely determined by the dissipation rate ε. From dimensional analysis, in the inertial subrange, the energy spectrum E(k) can be written as:

E(k) = α 1 ε 2 3 k -5 3 (1.21)
where α 1 is the Kolmogorov constant, with a value estimated between 0.5 and 0.6. That is the well known -5/3 power law in Kolmogorov's theory which characterizes the inertial subrange.

Previous researches for the urban canopy flow

In general, the urban canopy flow is studied using field measurement, wind-tunnel experiment or numerical simulations. In the past few decades, a great deal of research has been carried out by researchers using these approaches to explore the physical mechanisms of the turbulent flow over urban area. The purpose of this section is to briefly review the principal findings and conclusions that are helpful to evaluate and improve the numerical model and to analyze the turbulent characteristics within the urban canopy.

Field measurement

Field measurement means that the measuring devices are located outdoors and the measurement data collected from the survey area. In field measurement, all the scales of turbulence are obtained, up to the boundary layer scale.

The investigation of real flows in the urban environment began from the 1990s. Rotach (1993a[START_REF] Rotach | Turbulence close to a rough urban surface part II: variances and gradients[END_REF][START_REF] Rotach | Profiles of turbulence statistics in and above an urban street canyon[END_REF] studied the wind flow in the "Urban Climate Program" project of Switzerland over a W/H=0.83 street canyon. They found an inflection point in the mean wind profile at z = 1.2h, where h is the mean height of the canopy [START_REF] Rotach | Profiles of turbulence statistics in and above an urban street canyon[END_REF]. Later, a series of measurements were conducted in "Bubble" project of Switzerland [START_REF] Christen | Atmospheric turbulence and surface energy exchange in urban environments-Results from the Basel Urban Boundary Layer Experiment (BUBBLE)[END_REF][START_REF] Roth | Methodological Considerations Regarding the Measurement of Turbulent Fluxes in the Urban Roughness Sublayer: The Role of Scintillometery[END_REF][START_REF] Christen | Coherent structures in urban roughness sublayer turbulence[END_REF][START_REF] Christen | The Budget of Turbulent Kinetic Energy in the Urban Roughness Sublayer[END_REF]. They analyzed the flow over and within the urban canopy, as well as the surface energy balance. [START_REF] Christen | Atmospheric turbulence and surface energy exchange in urban environments-Results from the Basel Urban Boundary Layer Experiment (BUBBLE)[END_REF] did a lot of work about the TKE budget within the RSL. He pointed out that in the canopy, shear production is less significant, and turbulence is dominated by large coherent structures which is highly intermittent. The majority of TKE is transferred by turbulent and pressure transport from the roof layer where a large velocity gradient and an inflected wind profile are observed. The shear production dominates the TKE inside the shear layer. A large amount of TKE is exported from the shear layer by sweeps into the upper part of the street canyon and by ejections into the above-roof layer [START_REF] Christen | Coherent structures in urban roughness sublayer turbulence[END_REF]). As a consequence, the dissipation rate was found lower than local production. [START_REF] Christen | Coherent structures in urban roughness sublayer turbulence[END_REF] found that the inflection is above the top of the canopy between 1h and 1.5h. [START_REF] Davidson | Plume dispersion through large groups of obstacles-A field investigation[END_REF] carried out the first field survey of the flow characteristics around the staggered and aligned array of the three-dimensional cube in city area. The results of the comparison between the measured data and the wind tunnel experimental data were good.

In recent decades, a series of field measurement projects have been carried out in urban area. For example, the "COSMO" project in Japon (Kanda et However, only a few measurement points are generally available in real urban areas to conduct the analysis which cannot provide the complete three-dimensional flow field [START_REF] Giometto | Spatial Characteristics of Roughness Sublayer Mean Flow and Turbulence Over a Realistic Urban Surface[END_REF]). In addition, various infrastructures in urban areas may interfere with the understanding of measurement data. The inhomogeneous flow characteristics in the RSL raises questions about the use of point measurement. Recently, stochastic estimates are used to combine the wind tunnel data with the field measurement data [START_REF] Perret | Combining Wind-Tunnel and Field Measurements of Street-Canyon Flow via Stochastic Estimation[END_REF]) for the urban canopy flow studies, providing a direction for the widely use of field measurement data.

Wind-tunnel experiment

The wind-tunnel experiments are now a commonly used physical modeling method. In the field of aerodynamics research, wind tunnels are used to study the effects of air flow through the ad hoc objects. These tested objects are referred to as the wind-tunnel model and are instrumented with appropriate sensors to measure the velocity, aerodynamic forces, pressure distribution, or other aerodynamic-related characteristics. The measurement techniques can be Hot-Wire Anemometry (HWA), Particle Image Velocimetry (PIV) and Laser Doppler Velocimetry (LDV).

The wind-tunnel experiments conducted by the laboratory EnFlo and University of Southampton mainly focus on the flow over cubic array to study the development of turbulence inside the urban boundary layer. They have achieved many important results [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF][START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF][START_REF] Cheng | Flow over cube arrays of different packing densities[END_REF][START_REF] Reynolds | Measurements in an urban-type boundary layer[END_REF]. [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF] performed experiments over different height and arrangement cubic arrays with a same planning area density λ p = 25%. The results demonstrated the three-dimensionality of the turbulent flow inside the roughness sublayer (RSL). They also found that there is a greater drag to the flow in staggered array than in aligned array. [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF] studied the turbulence over staggered cube array under the neutral condition. The results confirmed that the dominant scale of turbulence inside the RSL is of the same order of magnitude as the roughness height. The eddy structure over the 3D cubic array is different from that over 2D roughness.The TKE budget results showed that the shear production acts as a principal energy source above the obstacles, and the turbulent transport is significant inside the RSL. [START_REF] Reynolds | Measurements in an urban-type boundary layer[END_REF] used PIV and LDV to measure the flow over staggered cube array. In the quadrant analysis, they found that sweep motions are dominant and contribute a large percent to the mean shear stress, whereas ejections are significant outside the RSL. The results of the two-point correlation analysis also show that the vertical velocity correlation is less dependent on the increased height, which is different from the streamwise velocity correlation.

Recently, the experiments performed in the wind tunnel of the laboratory LHEEA (Centrale Nantes) provide new dataset to analyze the turbulent flow and assess numerical simulations within the RSL (Blackman and Perret 2016; [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF][START_REF] Herpin | Investigation of the flow inside an urban canopy immersed into an atmospheric boundary layer using laser Doppler anemometry[END_REF]). Blackman and Perret (2016) carried out the wind-tunnel experiment in a boundary layer wind tunnel. The relationship between the large-scale structures in the flow and the small-scale structures inside the canopy was explored. In the another study of [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF], the energy transfer mechanism is analyzed through the TKE budget analysis, in which the dissipation is calculated from a LE-PIV method [START_REF] Sheng | A large eddy PIV method for turbulence dissipation rate estimation[END_REF]. The results confirmed that the shear layer induced by the underlying roughness produce and dissipate energy, and transfer the energy through other TKE budget terms. [START_REF] Herpin | Investigation of the flow inside an urban canopy immersed into an atmospheric boundary layer using laser Doppler anemometry[END_REF] analyzed dynamics of the flow in a staggered array with a 25% density. The power spectral density inside the canopy was analyzed using the LDV data. The spectral analysis provides more insight into the structure of turbulence near walls.

The advantage of the wind tunnel experiment is that as long as the fully developed turbulence is reached, the complex turbulence structure within the roughness sub-layer can be realistically simulated without relying on any turbulent closure model (e.g. in CFD). Additionally, complex urban morphology can be modeled in detail, including plants, street canyon and 3D urban models etc. Comparing to the field measurement, there is no limitation in the flow measurement position, and constant flow conditions can be maintained throughout the whole measurement. These experimental studies deepened the understanding of turbulent flow characteristics by elucidating the main turbulent processes that occur within the atmospheric boundary layer and provided a vastly useful dataset for validating the accuracy of numerical models (see section 1.5.3).

However, restrictions in experimental technique prevent from capturing the three dimensional flow field which often limits the analysis to one or two dimension and make researcher to rebuild missing information based on theoretical assumptions [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF][START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]). The assumptions used in the wind-tunnel experiments need to be carefully verified. Under the current situations, it is necessary to conduct numerical simulations to get more insight into the turbulent transport processes and verify the relevance of the assumption used in the analysis of experimental data.

Numerical simulations

Numerical simulation has been widely used in the study of turbulence and related processes in urban areas. It is expected to provide more insight into the turbulent structures and provide sufficient evidence to verify the relevancy of the assumptions employed in the experiments. The requirement for reliable prediction of unsteady flow dynamics has greatly spurred the development of numerical methods. [START_REF] Castro | Measurements and Computations of Flow in an Urban Street System[END_REF] pointed out the importance of unsteadiness of the turbulent flow in dispersion issues. In this section, numerical simulation studies of the unsteady flow including direct numerical simulation and large-eddy simulation over urban canopy are reviewed.

Direct numerical simulation

Direct numerical simulation (DNS) is a numerical method in which the Navier-Stokes equations are solved without using any turbulence model. This approach requires to resolve the flow down to the smallest turbulence length and time scale.

In the DNS studies over the 3D urban-like configuration, [START_REF] Yakhot | Turbulent flow around a wall-mounted cube: A direct numerical simulation[END_REF] used the immersed-boundary method to study the flow characteristics and coherent structure over one isolated cube in a fully developed turbulent channel with Reynolds number Re = 5610. The time-averaged data of turbulent mean flow, Reynolds shear stress, and turbulent kinetic energy are shown. The results confirmed that the vortex formed in front of, behind and on the both side of the cube generate the instability of the flow. In addition, negative production was observed in front of the cube where the main horseshoe vortex originates. [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]Coceal et al. ( , 2007a[START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF][START_REF] Coceal | Spatial Variability of Flow Statistics within Regular Building Arrays[END_REF]) published a series of papers that simulated turbulence features on aligned, square, and interlaced cubic arrays. The DNS results were well validated with the wind-tunnel measurement data from [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF] and [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF]. In [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF], three different types of cubic arrangements (aligned, squared and staggered) with a 25% packing density were used in the [4h × 4h × 4h], [4h × 4h × 6h] and [8h × 8h × 4h] computational domain with the Reynolds number based on the velocity at the top of the domain and the cube height is Re = 5000. Periodic boundary conditions were applied in the horizontal directions, and free-slip condition acted on the top, and no-slip boundary conditions were imposed in all solid surface. Good agreements were observed in the mean velocity and two-order statistics in comparison with the wind-tunnel data from [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF]. Besides, the mean flow pattern in the staggered array were more three dimensional than in the alined arrays and squared arrays. The section drag coefficient is constant with the vertical height in the squared arrays and alined arrays, while varies significantly in the staggered arrays. [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF] performed DNS study over a staggered cubic array in [16h × 12h × 8h] computational domain with a 25% packing density and Re = 5800 to investigate the turbulent structure of the turbulent flow. The results demonstrated the accuracy of DNS by comparing the mean velocity, turbulence intensities and spectra with the wind-tunnel experiment data (Cheng and[START_REF] Cheng | Near wall flow over urban-like roughness[END_REF][START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF]). The dynamic features and organized structures of the unsteady turbulent flow are analyzed in details. The presence of the low-momentum regions and ejection and sweep events are associated with hairpin-like vortices. Later, [START_REF] Leonardi | Channel flow over large cube roughness: a direct numerical simulation study[END_REF] performed DNS over staggered arrays of cubes with various density ranging from λ p = 4% to λ p = 25% in [8h × 6h × 8h] and [12h × 9h × 8h] computational domain with Re = 7000. The obstacles were represented using immersed boundary method. They examined the drag force coefficient difference in the different arrangement to provide a reference for the urban canopy model which treats the urban canopy as a porous medium modeled in terms of a drag force distributed in heights. (e.g. [START_REF] Coceal | A canopy model of mean winds through urban areas[END_REF].

The resolution in DNS requires a very narrow mesh to solve the smallest dissipative scales, the size of the mesh decreasing as the Reynolds number increases. Through the above introduction, we can find that although DNS is considered to be the most accurate numerical method, unfortunately, in most cases, conducting DNS in the high Reynolds number ABL flow is still an impossibility.

Large-eddy simulation

As described previous, DNS is computationally expensive and only limited to low Reynolds number flow applications. Large-eddy simulation (LES) explicitly resolve the large-scale flow and uses a subgrid-scale (SGS) model to parameterize the effect of the small-scale motions. It has a good advantage in reproducing accurately the unsteady features of the turbulent flow, better than the low-order method like Reynolds-averaged Navier-Stokes (RANS), which often fails to correctly capture the dynamic features of the flow field (Cheng et al. 2003, Xie andCastro 2006). [START_REF] Castro | Measurements and Computations of Flow in an Urban Street System[END_REF] demonstrated that LES is a powerful tool comparing to DNS and wind tunnel data, and pointed out that the careful use of LES can produce generally excellent agreement with wind tunnel and DNS results.

In the atmospheric boundary layer and canopy studies, LES approach is widely used. It has been applied to the atmospheric boundary layer for diagnosing and predicting turbulent flows since the early 1970s [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[END_REF][START_REF] Deardorff | Numerical investigation of neutral and unstable planetary boundary layers[END_REF]. For urban canopy studies, [START_REF] Kanda | Large-Eddy Simulations on the Effects of Surface Geometry of Building Arrays on Turbulent Organized Structures[END_REF] performed a set of LES simulations over staggered and squared arrangement of cubic arrays with different densities (from λ p = 0.03 to λ p = 0.44) and different height. They found that the drag coefficient was sensitive to the building area density in staggered array, while insensitive in square array. The turbulent organized structure was analyzed in the study through the flow visualization and quadrant analysis. [START_REF] Xie | LES and RANS for turbulent flow over arrays of wall-mounted obstacles[END_REF] carried out LES over staggered cubic arrangements with a 25% packing density in the [4h × 4h × 4h] computational domain. The simulations were performed under a range of Reynolds number from 5000 to 5000000. Good agreement of the mean velocity and turbulent statistics was observed in comparison with the DNS results [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]) and the wind-tunnel experiments [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF]. The result confirmed the weak influence of Reynolds number on the canopy flow. Later, Xie et al. (2008) performed LES over a random height roughness array. The mean velocity and turbulent intensities were validated against the wind-tunnel data from [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF]. This study found that the tallest buildings generate a disproportionate contribution to the surface drag, and the distinct blocks can significantly influence the canopy flow. Bou-Zeid et al. Turbulent kinetic energy budgets can provide a detailed turbulent energy information above and within the canopy. [START_REF] Giometto | Spatial Characteristics of Roughness Sublayer Mean Flow and Turbulence Over a Realistic Urban Surface[END_REF] performed LES in conjunction with the immersed boundary layer method on a real urban-like canopy model. The mean velocity data is in good agreement with the field measurement data. The novelty of this study is to quantify the impact of the non-measurable terms on TKE budget, and the pressure transport is found to be significant in the near-wall regions. [START_REF] Tomas | Stable stratification effects on flow and pollutant dispersion in boundary layers entering a generic urban environment[END_REF] performed LES over squared array in both stable and neutral conditions to study the stratification effects on the flow and pollutant dispersion. Good agreements were found in the profiles of mean streamwise velocity and mean Reynolds stress in comparison to the wind-tunnel data from [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF]. The results confirmed that the subgrid turbulent transport of TKE is negligible in the flow field, and subgrid dissipation account for 40% of the total dissipation of TKE. However, so far, no numerical simulation has been used to explore the TKE budget in a staggered array.

There are other works worth mentioning briefly, although they are not the main focus of current work, such as using LES simulations to study the dispersion issues over urban canopy (Boppana et They demonstrated that exponential mean velocity profiles could describe the mean velocity within the RSL under a wide range of arrangements. However, this exponential profile has been doubted by [START_REF] Castro | Are Urban-Canopy Velocity Profiles Exponential?[END_REF], who examined various datasets and mentioned that Yang's conclusions may be accurate when the grid resolution is sparse or the domain contains various building heights.

Finally, in the description of these urban canopy turbulence study, there are still many unanswered questions which need to be explored. Some unresolved issues such as using numerical simulation to analyze the TKE budget will be fully discussed in the later chapters.

Chapter 2 Numerical model and verification

The scope of this chapter is to introduce large-eddy simulation models, including governing equations, subgrid-scale (SGS) models, and numerical methods. Numerical selections such as mesh resolution and SGS models are tested in the preliminary simulation. The results using different SGS models are compared and discussed in regard to literature studies. The purpose of this section is to select the most appropriate simulation parameters to perform LES simulation.

Introduction

Numerical simulation has been widely used in turbulence research in urban areas. As described in the chapter 1, DNS is computationally expensive and is therefore limited to low Reynolds number applications. LES requires less computer resources than DNS and is considered to be the most beneficial and feasible numerical method for detailed analysis of unsteady turbulence in the atmospheric boundary layer.

However, the results obtained from the LES simulation depend on parameters such as mesh resolution, SGS model, and numerical scheme. Choosing the most appropriate simulation parameters is of utmost importance for properly modeling the turbulent flow. As far as the SGS model is concerned, the main drawback of the traditional eddy-viscosity SGS model is that it prescribes a constant model coefficient, which should vary in space and time. Additionally, these models prevent from simulating instantaneous energy backscatters through the cutoff length scale, which were evidenced both in the shear layer that develops at the top of the cube and in the overlying boundary layer [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]). The dynamic procedure proposed by [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF] overcomes this constraint and allows the model coefficient to change dynamically, which makes the dynamic SGS model be the best candidate to operate for the urban environmental simulation.

In this chapter, we first introduce the governing equations of incompressible turbulent flow, including the continuity equation and the momentum equation, as well as the common SGS model choices (see section 2.2). Simulation parameters used in OpenFOAM, such as the flow solver, boundary conditions, and numerical schemes, are then described in section 2.3. The comparison results of the different SGS models and of the different mesh resolution are discussed in regard to the literature study in Section 2.4. Section 2.5 summarizes the concluding remarks of the present chapter.

The LES equations

Resolved scale (filtered)

Un-resolved scale (subgrid scale) The filtered continuity and momentum conservation equations of the incompressible flow write as:

∂ u i ∂ x i = 0 (2.1) ∂ u i ∂t + ∂ u i u j ∂ x j = - 1 ρ ∂ p ∂ x i + ∂ ∂ x j (ν ∂ u i ∂ x j -τ i j ) (2.2)
where u i (for i = 1 ( u), 2 ( v), 3 ( w)) are the filtered velocity components in streamwise, spanwise and vertical directions, respectively, and p is the filtered pressure. ρ is the density and ν is the kinematic viscosity. τ i j is the subgrid-scale stress tensor defined as : τ i j = u i u ju i u j .

Subgrid-scale modeling

The LES application is sensitive to the choice of the subgrid-scale (SGS) model, especially in the vicinity of the walls. The proper selection of the SGS model will help to improve the accuracy and fidelity of the large eddy simulation. There are two main types of SGS models: one is the traditional eddy-viscosity SGS models such as Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment[END_REF]) and subgrid kinetic energy model [START_REF] Deardorff | Stratocumulus-capped mixed layers derived from a threedimensional model[END_REF]; the other type refers to dynamic models such as dynamic Smagorinsky model [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF]).

Smagorinsky SGS model

The well-known Smagorinsky SGS model is initially proposed in 1963 by Joseph Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment[END_REF]. It is an eddy-viscosity SGS model, and relies on the viscous analogy and the mixing length assumptions. The deviatoric (dev) part of the SGS stress tensor is defined as:

dev(τ i j ) = τ i j - 1 3 τ kk δ i j = -2 ν sgs S i j (2.3)
and the SGS kinetic viscosity ν sgs is modeled as:

ν sgs = (C S ∆) 2 | S i j | (2.4)
where C S is the Smagorinsky model coefficient, which is usually specified to be a constant value according to the flow type. ∆ = 3 √ ∆x∆y∆z is the mesh-filter width. | S i j | = 2 S i j S i j with the resolved strain rate tensor S i j defined as:

S i j = 1 2 ∂ u i ∂ x j + ∂ u j ∂ x i (2.5)
The Smagorinsky SGS model has been widely used in atmospheric flow field research (Xie et al. 2008;[START_REF] Boppana | Large-Eddy Simulation of Dispersion from Surface Sources in Arrays of Obstacles[END_REF]). The Smagorinsky SGS model is based on the assumption of a local equilibrium of turbulence at small scales, meaning the balance between shear production and dissipation [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment[END_REF][START_REF] Sullivan | A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows[END_REF]). The Smagorinsky model coefficient C S is prescribed to be a constant value. This prevents from the calculation of backscattering energy and results in an overestimated effect of model viscosity near walls. [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[END_REF] used C S = 0.1 in the turbulent channel flow after considering that the value of C S = 0.23 [START_REF] Lilly | On the application of the eddy viscosity concept in the inertial sub-range of turbulence[END_REF]) causes excessive damping of large-scale fluctuations. [START_REF] Mason | On the magnitude of the subgrid-scale eddy coefficient in large-eddy simulations of turbulent channel flow[END_REF] found that C S = 0.2 gave good results when the mesh resolution was sufficiently fine. In OpenFOAM, C S = 0.167 is the standard value which will be used in this study.

Turbulent kinetic energy SGS model

Turbulent kinetic energy (TKE) SGS model [START_REF] Deardorff | Stratocumulus-capped mixed layers derived from a threedimensional model[END_REF]) is also an eddy-viscosity model in which the subgrid-scale kinetic energy (k sgs ) is obtained by solving a transport equation (equation 2.9).

The subgrid-scale stress tensor τ i j is approximated as:

τ i j = 2 3 k sgs δ i j -2ν sgs S i j (2.6)
where the SGS kinetic energy k sgs is:

k sgs = 1 2 τ kk = 1 2 ( u k u k -u k u k ) (2.7)
and the SGS eddy viscosity ν sgs is:

ν sgs = C k k sgs ∆ (2.8)
The SGS kinetic energy k sgs transport equation writes:

∂ k sgs ∂t + ∂ ( u j k sgs ) ∂ x j = 2ν sgs S i j ∂ u i ∂ x j + ∂ ∂ x j (ν + ν sgs ) ∂ k sgs ∂ x j -C ε k 3/2 sgs ∆ (2.

Implementation of a dynamic Smagorinsky SGS model

The dynamic Smagorinsky SGS model has been implemented in OpenFOAM in the present work in order to improve the accuracy of the SGS model near walls and to capture the backscattering kinetic energy in LES simulation. In the dynamic approach a test-filter G = G G is defined, the characteristic width of which is ∆ = a × ∆. The momentum equation (2.2) filtered with the filter G involves a new SGS stress tensor T i j = u i u ju i u j which is unknown as well as τ i j . In the dynamic Smagorinsky model the deviatoric part of the SGS stress tensor at grid-and test-filter scales are modeled in a similar way as in equations (2.3) and (2.4) assuming that the model coefficient C is independent of the filter scale:

τ i j - 1 3 τ kk δ i j = -2C ∆ 2 | S i j | S i j (2.10) T i j - 1 3 T kk δ i j = -2C ∆ 2 | S i j | S i j (2.11)
Germano (1991)'s algebraic identify defines the tensor L i j as:

L i j = T i j -τ i j = u i u j -u i u j (2.12)
which can be calculated explicitly by applying the filter G to the resolved velocity field.

From equations (2.10), (2.11) and (2.12), one can write:

L i j - 1 3 L kk δ i j = CM i j (2.13)
where

M i j = -2 ∆ 2 (a 2 | S i j | S i j -| S i j | S i j
) with a = 2 in the present study.

As suggested by [START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF], the dynamic coefficient C is computed at each time step and each computational cell using least-squares method, through the relation: 2). However, the dynamic procedure may locally induce large negative values of the model coefficient, giving rise to the numerical instabilities. In order to avoid this problem, the total viscosity (ν + ν sgs ) is forced to be locally a non-negative value (the code of this procedure is given in Appendix B).

C = L i j M i j M i j M i j (2.14)

Numerical method details

In this work, LES simulations are carried out using OpenFOAM © 2.4.0 (Open source Field Operation And Manipulation) [START_REF]Openfoam user guide[END_REF], an open-source development environment based the C++ library aiming for solving various complex physical problems using finite volume method including fluid mechanics, thermodynamics, heat and mass transfer, etc. [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]). The code can be parallelized using an MPI interprocessor communication protocol. The great benefit of OpenFOAM is that it is released under the terms of the GNU general public license (GPL) and the object-oriented design makes it easy for users to freely redistribute and/or modify their own models and solvers for their specific purposes. Besides, OpenFOAM provides a variety of applications and simulation examples to help users better understand the principles and usage scenarios of the program.

Finite volume method

The finite volume method (FVM) is a method of representing and evaluating partial differential equations (PDEs) in the form of algebraic equations [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]). The filtered continuity equation and momentum equations of the incompressible flow (Eq. 2.1 and 2.2) write in finite volume method:

V ∂ u i ∂ x i dV = 0 (2.15) ∂ ∂t V u i dV + V ∂ u i u j ∂ x j dV = V - 1 ρ ∂ p ∂ x i dV + ν V ∆ u i dV - V ∂ τ i j ∂ x j dV (2.16)
The Gaussian divergence theorem is used to convert the volume integral of the divergence terms into a surface integral. The basic idea is to interpolate the volume center values to the values at face center, and then calculate the sum of these face center values. For example, in Figure 2.3 (in 2D for sake of clarity), the domain is divided into a finite number of small control volume (CV). The divergence term of variable ϕ in the computational nodes (such as at cell center P of the control volume) can be written as equation (2.17).

V ∇ • ϕ dV = S ϕ • - → n dS = ∑ f S f • ϕ f • - → n f = S w • ϕ w • -→ n w + S e • ϕ e • - → n e + S n • ϕ n • - → n n + S s • ϕ s • - → n s (2.17

Numerical schemes

Spatial schemes

In OpenFOAM, the user has several spatial interpolation scheme options (see Table 2.1). In the current work, interpolation involved in the gradient terms, laplacian terms and divergence terms are performed using the linear interpolation scheme (linear entry in Table 2.1). Linear interpolation is center difference and the weight calculated to each adjacent center is inversely proportional to its distance to the face. For example, see Figure 2.3, the face center value ϕ w can be calculated by interpolating the volume center value of two adjacent elements as:

ϕ w = a a + b (ϕ P ) + b a + b (ϕ W ) (2.18)

Time schemes

Time discretization is typically performed through discretizing the time integral equation. In OpenFOAM, the discretisation schemes for time can be selected from choices listed in Table 2.2.

In this work, the second order implicit backward scheme named by the keyword backward is used. For the quantity ϕ, the time derivative can be evaluated by:

∂ ϕ ∂t (t) = 3 ϕ n -4 ϕ n-1 + 2 ϕ n-2 2 ∆t (2.19)
where ∆t is the time step for the simulation. The superscripts denote time levels, ϕ n = ϕ(t), ϕ n-1 = ϕ(t -∆t) and ϕ n-2 = ϕ(t -2 × ∆t). In each time-step, the integration takes the value of 2 previous steps to obtain the solution.

OpenFOAM flow solver

The flow solver for the present LES study is modified from the PimpleFOAM solver. It is a combination of the PISO (Pressure Implicit with Splitting of Operator) algorithm [START_REF] Issa | Solution of the implicitly discretised fluid flow equations by operatorsplitting[END_REF]) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm [START_REF] Patankar | A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[END_REF]. It is well known that the time step required for the PISO algorithm to solve fast changing flow needs to be very small to keep the stability criterion of Courant number (Co) < 1, otherwise, the simulation will blow up. But the pimpleFOAM solver can speed up the simulation by using a larger time step and keep achieve good convergence. This is because SIMPLE (steady-state) processing is used to find the steady-state solution for each time step, especially when the nature of the solution is unstable. The code and flow chart of the PimpleFOAM algorithm is provided in Appendix A.

PimpleFOAM algorithm

The sequence of events can be summarized as follows in order to compute the solution at t n+1 = t n + ∆t 1. Use the initial guess solution at time t n for pressure, velocity and mass flow rate field ( p (n) , u (n) and ṁ(n) , respectively). 6. Using the latest available velocity and pressure field, calculate the coefficients of the momentum equation and solve it explicitly.

Construct momentum equation:

∂ u i ∂t + ∂ u i u j ∂ x j -∂ ∂ x j (ν ∂ u i ∂ x j -τ i j ) = -1 ρ ∂ p ∂ x i
7. Construct the pressure matrix using the momentum matrix with the updated velocities and obtain a new pressure correction field p ′′ .

8. Update the pressure, velocity, and mass flow rate field.

9. Judge whether the number of corrector steps exceeded (nCorrectors)

10. Judge whether the solution converged 11. Set the solution at time t n+1 to be equal to the converged solution and set the current time t n+1 to be t n .

12. Advance to the next time step.

13. Go back to step 1 and repeat until the last time step is reached.

Mean velocity driving the flow

In the present simulation, the flow is driven by balancing the flow speed to the desired velocity U bar through adjusting the pressure gradient at each time step. That is to say, the predicted velocity u * during the simulation is adjusted by comparing to the desired velocity U bar , and then one obtain the new velocity field u * * .

The sequence of events is summarized as follows:

1. Set desired mean velocity U bar .

2. Solve implicitly the momentum equation to obtain the velocity field u * and read the initial pressure gradient (∇ p (n) ).

3. Calculate the pressure gradient increment (∇ p ′ ) by adjusting the volume averaged velocity to the desired value U bar .

4. Then, add the pressure gradient source to the next loop ∇ p * = ∇ p (n) + ∇ p ′ .

5. Update the velocity field u * * .

Boundary conditions

Boundary conditions setting is a crucial step in numerical simulation. The role of boundary conditions in modeling is not only to construct the geometric objects, but also to affect the solutions associated with boundary conditions. In this section, the details about the boundary conditions will be described.

Periodic boundary conditions

Periodic boundary conditions (PBCs), also known as cyclic boundary conditions, mean that all variables are periodic at the boundaries of the computational domain. That is to say, the variable ϕ at the inlet and outlet has the equivalent value as equation (2.20).

ϕ(outlet) = ϕ(inlet) (2.20)
Although the computational power of computers are rapidly growing, it is still unrealistic to perform numerical simulations in very large computational domains. Using a finite field domain to reflect an infinite canopy can only be achieved by applying periodic boundary conditions. In this work, the periodic boundary conditions are imposed in span-wise and streamwise directions to simulate an infinite cube array.

No-slip boundary conditions

In fluid dynamics, the no-slip condition for viscous fluids assumes that the velocity of the fluid at the solid boundary is zero, which can also be an example of the Dirichlet boundary condition. A schematic plot of the no-slip boundary condition is shown in Figure 2.4.

In this work, the no-slip conditions are prescribed at the floor and on all obstacle surfaces for velocity. Free-slip boundary conditions in fluid dynamics assume that there is no friction between the fluid and the solid boundary.A schematic diagram of the free-slip boundary condition is shown in 2.5. Besides, free-slip condition is set at the top boundary of the domain.

Zero velocity at surface

In the present simulation, the free-slip condition is set for the boundary top where the fluid is not affected by the friction of the ground.

Zero gradient boundary conditions

A zero-gradient boundary condition is Neumann boundary condition, which prescribes the gradient of a quantity on a boundary to zero (see equation 2.21):

∂ ∂ - → n ϕ = 0 (2.21)
In the present simulation, the pressure is set to zero gradient boundary conditions at the floor and on all obstacle surfaces.

Preliminary simulation 2.4.1 Computational domain

Choice of the SGS model

Since LES only resolves the large-scale motions of the flow and the small-scale motions are modeled by an SGS model, it is necessary to choose the most suitable SGS model to perform the large eddy simulation. In this section, preliminary simulations are performed to evaluate three SGS models:

Standard deviation of velocity

Reynolds shear stress

Mesh generation and choice of the resolution

Basic mesh generation

The mesh used here is generated by OpenFOAM's meshing tool: blockMesh and Snappy-HexMesh. These tools generate hexahedral mesh cells.

Figure 2.11 shows an example of the mesh generated in OpenFOAM, where the blank square area is occupied by the cube. In the horizontal directions a uniform mesh h/16 is used. In the vertical direction, the mesh has a uniform resolution (h/16) from the ground up to z = 1.5h, and is stretched in the vertical direction. The maximum mesh expansion ratio reaches 4 at top of the simulation domain.

x / h z / h

Sensibility of the mesh

This section aims to compare the mean velocity in three mesh resolutions: ∆ = h/16, ∆ = h/32 and ∆ = h/64. The mesh settings gives about 0.16 million cells for the ∆ = h/16 simulation, 0.98 million cells for the ∆ = h/32 simulation and 9.50 million cells for the ∆ = h/64 simulation. The sensitivity of the time-averaged streamwise velocity to the mesh resolution at locations P0-P4 is shown in Figure 2.12. The computed velocities are almost identical at mesh resolution of ∆ = h/32 and ∆ = h/64, except for the lower part slightly above the ground. The minimum mesh resolution ∆ = h/16 gives a smaller velocity compared to other mesh resolutions near the ground (Figure 2.12 (a)) and in the middle of the canopy (Figure 2.12 (c)). 2006) have confirmed that simulations over cubic geometry are generally well resolved with a ∆ = h/32 mesh compared to a ∆ = h/64 mesh. [START_REF] Xie | LES and RANS for turbulent flow over arrays of wall-mounted obstacles[END_REF] also pointed out that a LES with a ∆ = h/32 mesh can successfully captured the peak of σ w around the top of cube. [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF] performed a DNS with ∆ = h/32 mesh in a [16h × 12h × 8h] computational domain and got good agreements in comparison with wind-tunnel data from [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF]. Therefore, the mesh resolution of ∆ = h/32 is applied for the present simulation to resolve the flow variation that develops within the roughness sublayer.

Coceal et al. (

Refined mesh around the cubes

For high Reynolds number flow, an inertial sublayer exists in the regions z + > 30, where the mean velocity over a smooth wall is logarithmic (see equation 2.22):

U + = 1 κ ln (z + ) (2.22)
Closer to the wall, the viscous sublayer is the region where z + < 5, and the mean velocity has the following relation:

U + = z + (2.23)
Here, U + is the dimensionless velocity (i.e. the streamwise velocity u parallel to the solid surface divided by the friction velocity u τ , U + = u/u τ ). κ is the Von Kármán constant. z + is the distance to the wall made dimensionless using the frictional velocity u τ and the kinematic viscosity coefficient ν, z + = z × u τ /ν.

In the current preliminary study, the shortage induced by the mesh resolution ∆ = h/32 have been pointed out (Figure 2.13). Using the present mesh resolution (∆ = h/32), the minimum value of z + is approximately 5, which means that the size of the first layer of cells is larger than the width of the viscous sublayer, so that, there are no sufficient cells above the top of the cube to resolve the flow variation inside the viscous sublayer. Usually, as no wall function is used here, at least 3 points have to be in the range z + < 5 to correctly resolve the viscous sublayer. Accurately simulating the flow characteristics within the near-wall region is a requirement for the wall-bounded turbulence simulation. Tomas et al. (2016) used a minimal vertical mesh spacing of ∆ = h/100 at the top of blocks, aiming to resolve the flow features inside the viscous sublayer.

P0

Figure 2. [START_REF] Castro | Are Urban-Canopy Velocity Profiles Exponential?[END_REF] The logarithmic profile for the mean velocity above the cube, before refining the mesh around the cube (White circle) and after refining the mesh around the cube (Red point), U + = z + is in Black line.

x / h z / h Therefore, to accurately deal with the turbulent flow inside the near-wall region, a refined mesh needs to be applied. However, it should be noted that as the mesh is refined, the number of cells will increase greatly and the computational cost will increase accordingly. Considering the computational cost, the vertical mesh spacing of ∆ = h/128 is chosen on the first two mesh layers above the top of the cube, and ∆ = h/64 mesh acted on the other near-wall regions as shown in Figure 2.14. The reason of applying a finer mesh on top of the cube is due to the larger velocity gradient and stronger shear stress in the shear layer lying above the cube. After refining the mesh in near-wall regions, the velocity profile above the cube is clearly modified well above the viscous sublayer (see Figure 2.13). Now, values of z + below 5 have three points and follow the linear law (equation 2.23). It indicates that the present ∆ = h/128 mesh truly help to capture the viscous sublayer flow.

Concluding remark

In this chapter, we first introduced the LES equations, followed by three commonly used SGS models. Then, the numerical method of the simulation such as finite volume meshed, numerical schemes, flow solver and boundary conditions are described.

To choose the most appropriate simulation parameters, preliminary simulations were performed. The turbulent statistics using these SGS models, including time-averaged streamwise velocity, Reynolds shear stress, and standard deviation of the velocity were compared to the literature data and the differences analyzed. The comparison showed that the dynamic Smagorinsky SGS model provides the closest results to the literature data, especially in the canopy. This demonstrated that the dynamic Smagorinsky model can significantly improve the prediction results of numerical simulations.

The meshing generation method in OpenFOAM and the sensitivity to the mesh resolution of the mesh are discussed in this chapter. The purpose is to select the most appropriate mesh resolution which can use less computer resource and well resolve the flow characteristics within the viscous sublayer. After analysis, the vertical mesh spacing of ∆ = h/128 is chosen on the first two mesh layers above the top of the cube, and ∆ = h/64 mesh acted on the other near-wall regions. , where the reference position is in the wake of the cube (P1) at the height of 1.5h. It clearly shows that the large-scale turbulent structure above the cube is blocked by the finite size of the computational domain, although applying periodic boundary conditions here. This figure shows that a larger computational domain is needed to provide complete turbulence structure information. In reference to [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF], the computational domain [16h × 12h × 8h] is chosen in the next chapter. The statistics will be validated against the literature data, and the 3D flow characteristics will be revealed in the analysis.

Chapter 3

Validation and discussion of LES model for simulating urban canopy flow

Introduction

The first purpose of this chapter is to evaluate the performance of present LES model. The second purpose is to describe the 3D features of the mean flow. To achieve these objectives, mean flow and second-order statistics are validated against the wind-tunnel experimental data and DNS data from the literature. The velocity skewness and the energy spectrum are given to explain the transport and the distribution of turbulent kinetic energy. Finally, a detailed analysis of the 3D structure of the mean flow is performed by flow visualization technique of the streamline, Q-criterion and vorticity analysis.

Numerical simulation setup 3.2.1 Simulation domain

Simulation grid set-up

Along with the size of the computational domain, the mesh size is another important parameter since it determines the size of the resolved eddies and the ability to account for large velocity gradients that occur near the walls. Comparing their numerical results obtained with a grid spacing ∆ = h/32 and ∆ = h/64, Coceal et al. (2006) concluded that simulations over cubic geometry are generally well resolved with ∆ = h/32. Xie and Castro (2006) also pointed out that a LES using ∆ = h/32 was able to capture the peak of vertical velocity fluctuations near the top of the canopy. In their DNS with ∆ = h/32 Coceal et al. (2007b) also obtained results in good agreement with wind-tunnel data [START_REF] Cheng | Near wall flow over urban-like roughness[END_REF]. However, our own sensitivity analysis (see Section 2.4.3) showed that a grid spacing ∆ = h/32 fails to resolve the flow gradients that occur within the viscous sublayer. The accurate description of near-wall region requires refined mesh as also suggested by [START_REF] Tomas | Stable stratification effects on flow and pollutant dispersion in boundary layers entering a generic urban environment[END_REF] who used a vertical grid spacing of ∆ = h/100 at the top of the obstacles. Hence, in the present study, a regular grid ∆ = h/32 is first selected from the floor up to z/h = 1.5 and then gradually expanded in the vertical direction throughout the rest of computational domain. Additionally, to accurately deal with the flow in the near-wall regions, a vertical grid spacing of ∆ = h/128 is applied at the first two mesh layers above the top of cubes while a grid spacing ∆ = h/64 is applied on the other near-wall regions (see Fig. 3.2). All of these mesh settings give about 28 million cells for the present simulation.

Simulation running time setting

The simulation is performed with a time step of 0.00076T over a total duration 500T where T = h/u τ can be interpreted as the eddy turnover time for the largest eddies shed by the cube [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]). After a first period of 300T , which is longer than the duration reported in previous studies [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]), the flow has reached a statistically steady state. The simulation is further continued during 200T in order to compute high-order statistics and turbulent kinetic energy budget. There is no standard rule for selecting the time over which the statistics have to be carried out.

Simulation parameters

The value of the Reynolds number Re based on the velocity at the top of the domain and the domain's height (8h) is about 50000, and the friction Reynolds number based on the friction velocity u τ (section 1.4.3) and h is Re τ = 481. Simulation details are summarized in Table 3.1.

Assessment of numerical approach

In order to evaluate the quality of our LES results, the mean vertical profiles of the longitudinal velocity and of the Reynolds shear stress are examined at four locations P0, P1, P2 and P3 indicated in Fig. 3.1. These locations are extracted from the simulation to be compared to various experimental and numerical literature data. In the following, the time averaging of the resolved fields is denoted by an overbar (..). The fluctuation with respect to this average is denoted by a prime symbol so that the resolved velocity component can be decomposed as u i = u i + u ′ i . The index i take the value 1, 2, or 3 referring to the streamwise, spanwise and vertical components, respectively. The turbulent Reynolds shear stress (-u ′ w ′ ) above the canopy and down to z/h = 0.2 at P1 and P2 corresponds fairly well to the DNS of [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF] and to wind-tunnel data (Figs. 3.4 (a), (c)). Upstream of the cube (P2), a local maximum is observed in our LES results at z/h = 0.1 (Fig. 3.4 (c)). [START_REF] Xie | LES and RANS for turbulent flow over arrays of wall-mounted obstacles[END_REF] made the same observation and suggested that it was due to the viscous sublayer that was well resolved close to the wall at this position. However, at P2 for z/h = 0.1 a reverse flow is observed in Fig 3 .3 (b). A peak is also observed in the standard deviation of the longitudinal velocity component (Fig 3 .6 (c)). Figure 3.5 shows that these peaks are more likely to be linked to the upstream recirculation area (P2) located upwind of the cube (more details about this area are given in Sec. 4.5, see also For the positions presented here, the maximum Reynolds shear stress is located at the cube height in both simulations and experiments. In the wake of the cube (at P1) the numerical simulations underestimate this maximum in comparison to experimental data of Castro et al. (2006) while the PIV results by [START_REF] Blackman | Non-linear interactions in a boundary layer developing over an array of cubes using stochastic estimation[END_REF] show the lowest maximum among all data (Fig. 3.4 (a)). Discrepancies are also found for the standard deviation of the vertical velocity component σ w : the peak of σ w obtained with numerical simulations falls between experimental results (Fig. 3.4 (b)). Similar observations were made by [START_REF] Reynolds | Measurements in an urban-type boundary layer[END_REF] who pointed out two main reasons explaining the differences in the peaks of -u ′ w ′ and σ w . The first one relies on the fully developed state of the boundary layer. They state that for h/δ < 10%, the boundary layer has reached a fully developed state so the Reynolds shear stresses are reaching a minimum asymptotic value. That means the higher h/δ , the higher the shear stress peak. This is also observed here in Figs. At point P2, experimental data from [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF] and [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF] are also available to be compared with numerical results of σ w (Fig. 3.4 (d)). Both DNS and LES values of σ w fall between the two experimental databases. As for point P3, vertical profiles of -u ′ w ′ and σ w are rarely reported in the literature. The LES shear stress profile is in rather good agreement with PIV data above the canopy and down to z/h = 0.6. Below, large discrepancies are observed and LES results are higher than PIV data (Fig. 3.4 (f)). σ w profiles are showing large discrepancies at all heights (Fig. 3.4 (f)).

Mean stream-wise velocity

Reynolds stress components

The standard deviation of streamwise (σ u ) and spanwise (σ v ) velocity components are presented in Fig. 3.6. LES results show that the accuracy in simulating σ u is close to that of the DNS [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF]) at P1 and P2. At P1, below z/h = 1, numerical data are in good agreement with PIV data (Figs. 3.6 (a)) but significantly lower than results from LDV technique. In such a wake region, the lack of space resolution is probably to blame to explain these differences. At P2 and P3 (Figs. 3.6 (c),(e)) numerical results match very well with the LDV data from [START_REF] Herpin | Investigation of the flow inside an urban canopy immersed into an atmospheric boundary layer using laser Doppler anemometry[END_REF]. A phenomenon similar to the shear stress (i.e high σ u ) is observed in the recirculation zone near the ground (Fig. 3.7).

The standard deviation of the spanwise velocity (σ v ) is certainly the component of the Reynolds stress the most difficult to measure and simulate. It is to be reminded that σ v extracted from PIV is reconstructed from stereo-PIV leading to a underestimation of the velocity fluctuations. Large discrepancies can be observed comparing the results obtained in the various experiments (Figs. 3.6 (b),(d),(f)). σ v is significantly smaller than in the experimental data of [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF], but larger than in the experimental data of [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]. This is mainly explained by the measurements and simulation resolution differences or the difference in h/δ ratio. However, the LES seems to provide the right tendency in the profile of σ v and it is in rather good agreement with the wind-tunnel data from Herpin et al. ( 2018) at P2 and P3 and at P1 above and inside the canopy top.

From this analysis it is clear that the assessment of the numerical method cannot be based on a unique experimental data set. The detailed reason why the different experimental approaches lead to different results is out of the scope of the present work. However the comparisons presented in this section demonstrate that the overall agreement between our LES results and available data including DNS is satisfactory for the first and second order moments.

Skewness of velocity component

As explained in section 1.4.4, skewness is a measure of the asymmetry of a statistical distribution. For a velocity component u i , it is the third order moment of velocity fluctuation normalized using the standard deviation:

Sk u i = u i ′ 3 /(σ u i ) 3 (3.1)
Skewness provides insight into momentum transfer events such as sweeps and ejections [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]). [START_REF] Brunet | A wind tunnel study of air flow in waving wheat: single-point velocity statistics[END_REF] demonstrated over a plant canopy that strong positive Sk u associated with negative Sk w are the evidence of energetic downward sweep events while strong negative Sk u and positive Sk w are linked to ejections.

Profiles P0-P3 presented in Fig. 3.8 are in the symmetry plane, so Sk v ≈ 0 for all positions. At P1, a negative peak in Sk w and a positive peak of Sk u near z/h = 1 may be the sign of important downward motions (sweeps) in the shear layer after the cube. Below z/h = 0.8, components of skewness have lower values indicating a more homogeneous mixing induced by the wake recirculation behind the cube. In contrast, at P2 and P3, high negative values of Sk w and low positive values of Sk u are found below z/h = 1 clearly indicating the presence of sweeps. For all positions studied, above z/h ≈ 1.25, Sk u < 0 and Sk w > 0 meaning that the flow is dominated by ejections, similarly to a rough boundary-layer flow. These results are consistent with the observation of [START_REF] Finnigan | Turbulence in plant canopies[END_REF] and [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF].

Energy spectrum

The energy spectrum E(kz ′ ) of the streamwise velocity component above the cubes (at z/h = 1.62) is shown at P2 in Fig. 3.10 and compared with the wind-tunnel data from [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF]). The wavenumber k is defined as k = 2π f /U, with f the frequency and U the local mean streamwise velocity. The variable z ′ is defined as z ′ = zd with d denoting the zero-plane displacement estimated as the height at which the total drag acts (see section 1.1.2). This method has been widely debated because the fitting of the logarithmic profile using the resulting value of d implied large discrepancies in the value of the von Kármán constant κ (see [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF]). However in our simulation, this method gives generally accepted values d = 0.74h and κ = 0.4 (κ is calculated from Fig. 3.9). In Fig. 3.10, the red doted-line denotes the lowest expected normalized wavenumber (k min z ′ ) computed based on domain size and mean velocity and the purple doted-line denotes the highest resolved wavenumber from the LES computation (cut-off) calculated based on the grid size and u τ . The amplitude of the spectrum is normalized through the integration of the spectrum ((σ u ) 2 ). The same procedure is applied to the dataset of [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF] to be comparable.

The simulated energy spectrum matches very well with the wind-tunnel data (Fig. 3.10) and the inertial subrange (-5/3 slope) is accurately captured. According to the data from [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF] the cut-off wavenumber of the LES (k max ) seems to be at the limit of the inertial sub-range. As expected, for a normalized wavenumbers higher than k max z ′ , the energy of the small scales starts to drop faster in the simulation than in the experiment. For low wave numbers, the limited size of the LES domain is expected to lead to a poor simulation of the large scale structures. However, the energy deficit for kz ′ lower than k min z ′ is rather limited compared to wind-tunnel data.

P2

Mean flow analysis

The focus of this section is on the mean flow field. Figure 3.11 shows the vertical crosssection with mean velocity vectors (u, w) in a vertical xz plane through the middle of the cube.

In the plot, a clockwise recirculation structure is shown at the bottom of the windward side of the cube (at P2). The flow is downward along the front face of the cube. The reverse flow is identified in the lower part close to the surface ahead of the cube. This circulation structure has a nearly identical position than in the previous studies of DNS [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]) and LES [START_REF] Xie | LES and RANS for turbulent flow over arrays of wall-mounted obstacles[END_REF]). A separation point at one cube length behind the cube leads to a reverse flow around location P1 in the lower part of the canopy. A strong updraft is observed on the leeward face of the cube, which merges into the uniform and unidirectional downstream flow slight above the top of the cube. We will show in Section 4.5 that these reversed flows have a strong connection with the production of the turbulent kinetic energy. Figure 3.12 shows the streamlines in an horizontal transect (x,y) near the ground (z = h/40). In particular, two small recirculation are located at the upwind corner of the cube. A big horseshoe vortex is observed in the vicinity of the cube, but the head of the horseshoe vortex is not clearly visible in the present plot. The upstream flow is deflected by the presence of the cube as a result of the flow separation occurring at the frontal face of the cube. The flow pattern presented here is different from the flow over one cube configuration which show a pair of counter-rotating flows behind the cube [START_REF] Yakhot | Turbulent flow around a wall-mounted cube: A direct numerical simulation[END_REF]). Besides, the reattachment point of the separated flow does not appear in the downstream area.

Vorticity analysis of mean velocity

The vorticity is used to measure the curl of the flow. Mathematically, the vorticity of a three-dimensional flow is usually defined as the curl (or rotational) of the flow velocity vector, denoted by -→ ω : The vorticity of the mean velocity is presented in Figure 3.14. The horseshoe vortex can be identified in front of the cube in an ellipsoid elongated in the cross-wind direction with two legs into the gaps of the cube. Besides, the vorticity of the mean velocity is observed on the edging sides of the cube and in the shear layer slight above the top of the cube. This indicates that when mean flow passing around an obstacle, the presence of the object causes the mean flow to change the flow direction and produce vorticity, while the laterally deflected flow moves upward, merging with the flow in the shear layer.

- → ω = ∇ × - → u (3.

Q-criterion analysis of mean velocity

The Q-criterion [START_REF] Haller | An objective definition of a vortex[END_REF]) allows to consider the vortex movement through separating the matrix of the gradient of velocity (equation 3.3), which can be used to present the vortex structure.

Q = 1 2 (Ω i j Ω i j -S i j S i j ) > 0 (3.3)
where Ω i j = 1 2 (∇u -(∇u) T ) is the vorticity tensor, and S i j = 1 2 (∇u + (∇u) T ) is the strain rate tensor. They are the symmetric and antisymmetric part of the velocity gradient ∇u.

Chapter summary

The investigations performed here are in the continuation of the efforts of Chapter 2 to progress toward higher fidelity levels of the LES modeling. In the context of assessing the performance of the model, large-eddy simulation (LES) with dynamic Smagorinsky SGS model are performed in this study to investigate the flow over urban-like cubic roughness.

The LES code is firstly assessed in comparison with the data from DNS [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]) and wind-tunnel [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF]; Blackman and Perret 2016; [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF][START_REF] Herpin | Investigation of the flow inside an urban canopy immersed into an atmospheric boundary layer using laser Doppler anemometry[END_REF]). Good agreement is observed in the vertical profile of mean flow, high-order turbulent statistics and energy spectrum. The investigation on the mean flow illustrates the inhomogeneous characteristic of the flow inside the urban boundary layer. The maximum simulated Reynolds shear stress is slightly above the top of the cube, but in magnitude smaller than the experimental data [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF]. Standard deviations of the streamwise and vertical velocity component match well with the wind-tunnel data above the canopy. The standard deviation of the span-wise velocity component is significantly smaller than the experimental data, but the general features are reproduced. The skewnesses of velocity show that sweep events mainly occur inside the canopy, and ejection events contribute more above the canopy. The energy spectrum confirms the existence of inertial subrange and the limitations due to the computational domain and the mesh size.

Besides, the mean flow is analyzed by presenting the vertical structure, 3D streamlines, and Q-criterion and vorticity analysis. The circulation in front of the cube is observed at an identical position as presented in the previous studies of [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF] and [START_REF] Xie | LES and RANS for turbulent flow over arrays of wall-mounted obstacles[END_REF]. Then, the horizontal streamlines confirm that the flow pattern is different from the flow over one isolated cube. Here, no vortex is observed at the top of the cubes, but the strong vortical structure is observed in front of the obstacles, which indicates that the horseshoe vortex generated in the recirculation zone is due to this vortex structure. The horseshoe vortex is clearly observed by the visualization means of Q-criterion and vorticity in front of the cube.

From the analysis in this chapter, the reliability of the large-eddy simulation model has been fully confirmed. At the same time, through the analysis of the flow statistics, we can conclude that the turbulent flow motion in the urban canopy is highly inhomogeneous and may lead to high turbulent energy exchanges. In order to understand the mechanism of turbulent energy transfer in the RSL, the characteristics of the turbulent kinetic energy will be described and discussed in detail in the next chapter.

Chapter 4 Turbulent kinetic energy budget over urban canopy

The scope of this chapter is to analyze the turbulent kinetic energy (TKE) budget over the urban-like canopy and to demonstrate the relevance of ignoring non-measurable terms in the wind-tunnel experiment when calculating turbulent transport. A comparison of the TKE budget against the experimental data is shown here to assess the fidelity level of the present LES simulation. The 3D visualization of the TKE budget terms in the vicinity of the cube is presented in order to identify the locations where they contribute more. Finally, the decompositions of the turbulent transport terms are compared to the wind-tunnel data.

Introduction

The growing interest in understanding the turbulence transport mechanisms that develop within the urban environment has promoted the development of the turbulent flow dynamic studies in the urban canopy. In the previous studies, mean flow, high-order statistics and coherent structures of the turbulent flow have been deeply analyzed by numerical modeling, such as in the work of [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]Coceal et al. ( , 2007a[START_REF] Coceal | Spatial Variability of Flow Statistics within Regular Building Arrays[END_REF] They demonstrated that the pressure transport is very significant in the near-wall area and the turbulence is transported into the urban canopy. [START_REF] Tomas | Stable stratification effects on flow and pollutant dispersion in boundary layers entering a generic urban environment[END_REF] used the LES method to study the spatial-averaged TKE budget over a squared array and highlighted the negligible effect of the subgrid-scale transport.

However, obtaining all the terms used to calculate the TKE budget in wind-tunnel experiments remains a challenging issue. Depending on the experimental facilities some terms in the TKE budget may be only partially measured. For instance, Blackman et al. (2017) used PIV technique to measure the three components of the velocity in the central plane of the cube: this makes it impossible to obtain the spanwise gradients involved in the transport term. These uncertainties raise the questions about the accuracy of turbulent transport data from the wind tunnel experiments [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]). To verify the relevance of the hypothesis used to rebuilt or neglect missing terms in wind tunnel data analysis, the LES modeling is performed in a [16h × 12h × 8h] computational domain with the dynamic Smagorinsky SGS model. The analysis carried out here is in the continuation of the effort and discussion of Chapter 3.

The objective of this chapter is to evaluate the performance of the LES model by validating and analyzing the TKE budget and to describe the main physical turbulent processes that occur within the urban canopy. The TKE budget terms are compared to the available experimental data from the wind-tunnel measurements [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF][START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF][START_REF] Herpin | Investigation of the flow inside an urban canopy immersed into an atmospheric boundary layer using laser Doppler anemometry[END_REF]) (see section 4.3 and section 4.4). 3D visualization of the TKE budget terms in the vicinity of the cube is used to present the spatial characteristics of the energy transfer (section 4.5). To the author's knowledge, this is the first attempt to present a 3D visualization of TKE budget using LES over the urban-like staggered cubic array. Finally, the effect of ignoring non-measurable terms when calculating TKE budget in wind-tunnel measurement is described and discussed.

The governing equations

From the LES equations (2.1), (2.2), (2.3) and (2.5), and assuming that the turbulent flow has reached a fully steady state, the turbulent kinetic energy (TKE) budget writes as:

0 = -u j ∂ ∂ x j 1 2 u ′ i u ′ i A -u ′ i u ′ j S i j P - 1 2 ∂ ∂ x j u ′ i u ′ i u ′ j T r - 1 ρ ∂ ∂ x i u ′ i p ′ T p + ∂ ∂ x j 2ν u ′ i S ′ i j D ν -2ν S ′ i j S ′ i j ε r - ∂ ∂ x j u ′ i τ ′ i j T sgs + τ ′ i j S ′ i j ε sgs (4.1)
where A represents advection by the mean flow, P is the production, T r is the turbulent transport by resolved velocity fluctuation, T p is the transport by pressure fluctuation, D ν is the viscous diffusion, ε r is the resolved dissipation, T sgs is the SGS transport, and ε sgs is the SGS dissipation, which represents the transfer of energy between resolved-scale and subgrid-scale through the cut-off filter. All contributions to the TKE budget are directly and individually computed in the simulation. Quantities are time averaged over 200T in each cell of the domain, and then combined to have the TKE budget terms everywhere in the domain.

TKE budget around a cube

Production and dissipation terms

At all locations, above the height z/h ≈ 1.5, the production (P) and dissipation terms (ε sgs + ε r ) are the major contributions to the TKE budget and tend to balance each other with the increasing height. However, the equilibrium is never exactly reached and the production remains stronger (|P| > |ε sgs + ε r |) (see Fig. 4.2).

The production P reaches a maximum just above the cubes at location P1 (Fig. 4.1 (a)), where the shear stress above the cube is the strongest. This sharp peak is then rapidly decreasing above and below the top of the cube (z/h = 1). With downwind distance (at P3 and P2), the production peak decreases in intensity and shifts slightly above z/h = 1. The region of high production becomes thicker (see Fig In the simulation, dissipation generally acts in a mirror of production, but with lower absolute values. The maximum ratio of |P|/|ε sgs + ε r | is found at about z/h = 1 (at P1), slightly above z/h = 1 (at P2 and P3) and at about z/h = 2 (at P0) (Fig. 4.2). The maximum ratio reaches 1.6, 6, 2.4, 2.25 for locations P0-P3, respectively. This phenomenon indicates a strong participation of other TKE budget terms in this area: advection, turbulent transport, and pressure transport. Above the roof level, the ratio decreases slowly at all locations. On the contrary, dissipation is higher than production inside the canopy in the regions, 0.3 < z/h < 0.9 at P1, 0.2 < z/h < 0.8 at P2, and 0.2 < z/h < 1 at P3. This is mainly due to the small production. For example, P < 0 is observed near the ground at P3 (Fig.

Turbulent transport terms

The turbulent transport (T r + T sgs ) is the major sink of TKE at P1. It transfers energy from the thin shear layer ((T r + T sgs ) < 0 at z/h = 1.05) downwards inside the canopy and upwards into the roughness sub-layer (T r + T sgs > 0 at z/h = 0.9 and z/h = 1.2). Going downwind, at P3 and P2, the turbulent energy sink still remains and broadens up and down, but the position of maximum energy sink becomes higher. Just upstream of the cube (at P2), the region where the turbulent transport is negative extends from z/h = 1.7 down to z/h = 0.9. The turbulent transport acts as a source of TKE inside the canopy. Similar behavior of turbulent transport is seen at P3. Besides, the turbulent transport, acting upward and downward at P1 is clearly unbalanced toward the canopy further downwind. Finally, the T sgs is found to be negligible throughout the whole flow field (not shown here). This phenomenon was also noted in the aligned cube arrays by [START_REF] Tomas | Stable stratification effects on flow and pollutant dispersion in boundary layers entering a generic urban environment[END_REF].

Pressure transport term

The transport by pressure fluctuations (T p ) acts in similar way as the turbulent transport in front of the cube (at P1) (see Figure 4.1 (a)) by transferring energy up and down. In the gap between the cubes (at P3), T p is very weak above z/h ≈ 0.8. Inside the canopy, it acts as a sink of TKE (see Figure 4.1 (c)). The contribution of T p is much more significant in front of a cube (at P2) (see Figure 4.1 (b)) than behind (at P1). At P2, it is a dominant sink with a broad extent between z/h = 0.4 up to z/h = 1.7, and even larger in magnitude than dissipation near the top of the cube (0.8 < z/h < 1.2) (see Figure 4.1 (b)). This peak of |T p | indicates that the flow in staggered configuration cubes impinges the windward side and edges of the cubes resulting in large pressure fluctuations in front of the cubes.

Advection term

Finally, the advection (A) is found to be significant in the shear layer at locations P1 and P2. According to the location in the wake of the cube or the upstream recirculation area, it may be negative (at P1) or positive (at P2). The large sink at P1 contributes to balance the strong production in the shear layer above the cube. Inside the canopy, the contribution of the advection is found to be significant at P2 and P3. Because of the strong link between the advection and mean flow, the negative and positive A observed in the lower part of the canopy (0.1 < z/h < 0.5) at location P2 reflects the existence of recirculation zone. This results from the flow impinging on the surface of the obstacle and creating a reverse flow. Besides, high advection area is observed at the lower part of the canopy (0.1 < z/h < 0.6) at P3 (Figure 4.1 (c)): it acts as a primary source, larger than the production. This phenomenon indicates that the turbulent energy present in the gap between the cubes is created elsewhere and advected. Indeed, as will be seen later (section 4.5.4), the flow is diverted laterally, and then converges in this area with speed up. The production (P) shows a good agreement with the experimental data inside the canopy and above z/h = 2 (at P2) and z/h = 1.5 (at P3). For both positions, the LES production peak is larger than PIV data [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]). Regarding dissipation (ε r + ε sgs ), the numerical data are larger in magnitude than in [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF], but slightly smaller than the data from [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF]. The same reasons as mentioned in Sec. 3.3.2 can be invoked here to explain the differences observed: (i) the lower boundary-layer height induces a higher shear resulting in higher production and higher dissipation in magnitude, (ii) the experimental PIV results may underestimate the peaks due to a too coarse resolution. The higher value of P and (ε r + ε sgs ) slightly above the canopy (at z ≈ 1.4h) and their rapid decrease with the height for both LES and wind-tunnel data from [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF] clearly show the effect by the small δ /h. With a four times bigger δ /h in Blackman et al. ( 2017), P and (ε r + ε sgs ) are decreasing very slowly with the height.

TKE budget comparison

As shown in [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF], through the decomposition of the production terms, the combination of the shear stress and the velocity gradient (-2 u ′ w ′ S 13 ) is the main source of production. This is confirmed in our LES result in the upper part of the canopy (0.8 < z/h < 1), and above the canopy (z/h > 1) by comparing (-2 u ′ w ′ S 13 ) to P (see Fig.

TKE budget in the vicinity of the canopy

The TKE budget terms presented in Figure 4.1 displayed a great inhomogeneity inside the canopy layer and within the roughness sub-layer. This section aims at describing and interpreting the three-dimensional results from the LES modeling by extracting 2D planes and displaying 3D iso-contours plots to understand deeper the spatial repartition of the turbulence.

Production

Turbulent transport

As discuss in section 4.3.2, the T sgs is found to be negligible throughout the whole flow field. Therefore, only the 3D iso-coutour of turbulent transport (T r ) is presented here.

A large area of negative T r are observed in the flow (see Fig. 4.8 and Fig. 4.9 (b)). As discussed in sec. 4.3.2, the turbulent transport is a major sink to balance the production inside the roughness sub-layer. Near the ground, the negative T r is observed at around h/2 upstream of the cube and in the gap of the cubes where T r < -8 (see Fig. 4.9 (b)). Moving upwards from the ground to the top of the cube, the negative value of T r are observed along the edging sides of the windward side and along the four edges of the top face and extend from the trailing edging to the next cube. Comparing Fig. 4.7 (a) and Fig. 4.8, it can be seen that the negative T r area is almost in the same location as the high production (P = 25) area. At the cube's foot, positive T r is forming on both side of the cube. Going up to z = h/2, three high turbulent transport (T r ) areas are observed: one is slightly downstream of the leading edge of the top face, another is an extended region in the wake area between two cubes (near location P3) and the third is on the vertical edges of the windward side of the cube (Fig. 4.9 (a) and Fig.

Pressure transport

The transport by pressure fluctuations (T p ) is significant (T p > 15) in front of the cube either on the windward side of the cube or in a region elongated in the spanwise direction which skirts the cube's foot (see Fig. 4.11 (a)). Going upwards with the increasing height, high pressure transport (T p > 15) passes the leading edge and on the cube's top face. The high T p area described here are clearly linked to the impinging flow on the surface near location N1 and N2 (repelling nodes).

Near the ground, two span-wise elongated negative T p zones (T p < -5) are observed: one is at around 2h/3 upstream of a cube above point S1, the other is along the foot edge of the cube above location N1 (Fig. 4.11 (b) and Fig. 4.12). Looking at the mean flow in this region (Fig. 4.12), the negative T p seems to be linked to the recirculation region upwind of the cube when w > 0 ( N2) and w < 0 (N1). Significant negative T p is observed at the leading edge of the windward side in a large width. Negative T p is also observed in the middle height of the canopy between the cubes. No significant negative T p is observed in the shear layer and in the wake of the cube. This phenomenon indicates that the T p does not contribute to transfer the high Production in the wake. But it certainly acts as a sink in the gap between the cubes (see Fig

Advection

Inside the canopy, the positive advection (A) region forms a span-wise tube at the cube's foot (see Fig. 4.13 (a)) and extends on the windward side up to the top of the cube, and on the cube's left and right sides, reaching the vertical elongated region situated in the gap of the cubes where A > 5 (near P3) (see Fig. 4.13 (a) and Fig. 4.14). Going up from z = h/2, high positive A is observed on leading edges in a large width, and downwind along the top face in the shear layer.

Near the ground, a negative cross-wind strip where A < -5 is observed located approximately z = h/3 in front of the cube (Fig. 4.13 (b)). This location is slightly downstream of point S1 where high P was observed (Fig. 4.14). Besides, stronger negative advection areas (A < -10) are observed on the leading edge and trailing edges of the top face. They extend downstream the second row cubes and merge into the canopy along the windward side vertical edges (see Fig. In the wake, two highly negative A regions are observed emanating from the trailing edge, forming two layers above and below the high P region and merging with the two vertical edges of windward side of the next cube (see Fig.

Dissipation

The main area to dissipate TKE is in the vicinity of the cube and near the solid boundaries (see Fig. 4.17). Inside the canopy, large values are observed on the two vertical edges of the windward face of the cube, and develops laterally along the downstream direction. Slightly above the cube, the peak of dissipation is clearly shown in the shear layer (Fig. 4.18) which extends downwind reaching the second row of cubes (see Fig. 4.17 and Fig. 4.18). The intensity of (ε r + ε sgs ) is larger in the gap (P3) than in the wake of the cube (P1) (Fig. 4.18).

Note that the little points in this graph are due to the subgrid-scale dissipation, in which the dynamic Smagorinsky model gives large SGS stress tensor values.

Analysis of the turbulent transport term

The commonly used experimental instrumentation (i.e., HWA, LDV or 2D PIV) are unable to measure gradients in the three directions restraining the calculation of some terms of the TKE budget. In the wind-tunnel experiment of [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF], span-wise gradients of triple correlations are ignored in the computation of turbulent transport (T r ). However, the influence of ignoring these terms has not been examined. In this section, the consequences of ignoring the spanwise terms of T r are examined by describing all terms of the turbulent transport (T r ) captured in CFD.

The decomposition of turbulence transport

The decomposition of turbulent transport (T r ) at location P3 is shown in Fig

∆T r = - 1 2 ∂ u ′ u ′ v ′ ∂ y - 1 2 ∂ w ′ w ′ v ′ ∂ y - 1 2 ∂ v ′ v ′ v ′ ∂ y . (4.2)
∆T r is small above the canopy and in the wake of the cubes. However, it is significantly large in the gap of the cubes. Ignoring these terms will result in a value of T r being overestimated or underestimated.

Chapter summary

Large-eddy simulations (LES) with dynamic Smagorinsky SGS model were performed to investigate the transport mechanisms of the turbulent kinetic energy within roughness sublayer.

The local TKE budget is demonstrates the complex behavior of turbulent flow. A peak of production is observed in the wake of the cube where the strongest shear layer induced by the presence of the cubic roughness develops. The peak decreases downward as the shear layer increases. The contribution of production and dissipation which balance each other are the major contributions to the TKE budget above z/h > 1.5. The turbulent transport acts as a sink in the shear layer, but as a source above and below. The pressure transport is significantly large in front of the cube. The advection acts to gain the energy in front of the cube but to lose energy in the wake. Production shows good agreement with the wind-tunnel results both above and within the canopy at location P3. The interesting suggestion from the comparison is that the overestimated absolute value of dissipation may be due to the lower boundary layer thickness. Besides, the 3D iso-surface visualization of the TKE budget terms reflects the high spatial variation. High production mainly concentrates on the edges of the cubes and in the shear layer. Dissipation is mainly in the vicinity of the cube. Advection and turbulent transport redistribute the energy in the flow and are the major sink in the wake of the cube. Pressure transport acts as a source in the recirculation area in front of the cube.

Finally, the relevance of ignoring the span-wise gradients when calculating turbulent transport in the experiment [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]) is verified with the present simulation results. The comparison of the results at location P3 is presented and highlight the contribution of span-wise terms inside the canopy. Neglecting these terms in the LES induces an underestimation but a better agreement with the wind-tunnel data. The 3D iso-coutour of the discrepancy provide more insight into the spatial distribution of this difference. The discrep-ancy is generally small above the canopy, but much more significant in the gap between the cubes than in the shear layer.

Chapter 5

Drag-porosity approach: assessment and suggestions

The scope of this chapter is to assess the drag-porosity approach in the application of unsteady urban canopy flow modeling. Spatial-averaged statistics such as mean flow, highorder statistics, TKE budget, and coherent structure are compared to the obstacle resolved LES simulation results and literature data. Suggestions are given for improving the prediction performance of the drag-porosity approach.

Introduction

In recent years, the urban environment has become a prime topic of concern, and increasing attention has been given to the various environmental problems which occur in urban areas. Based on current computing power, widely used building-resolved CFD models are only applicable to street or neighborhood scales and are not possible on a regional or city scale.

For the study of real urban areas, the detailed knowledge of the buildings organization is unavailable and a drag-porosity approach may be preferable. The drag-porosity approach, as commonly used for forest canopies [START_REF] Dupont | Influence of foliar density profile on canopy flow: A large-eddy simulation study[END_REF], models the presence of obstacles and their influence on the turbulent flow by a drag force that depends on averaged morphological characteristics of the canopy. It has the advantage of reducing the computation costs and can be applied for example to generate realistic inflow conditions for studies at the scale of some buildings. The model has been proved to reproduce the turbulence statistics characteristics within the atmospheric boundary layer [START_REF] Maché | Analysis of the dynamical interactions between atmosphere and urban canopies of different densities using a drag force approach[END_REF][START_REF] Rodrigues | Numerical and experimental tracer dispersion study for footprint assessment[END_REF]). However, still only few LES works exist on the use and efficiency of drag-porosity approach to represent the turbulent transfers between urban canopies and atmosphere [START_REF] Maché | Représentation multi-échelle des transferts entre couche de canopée urbaine et atmosphère à l'échelle de la ville[END_REF]). On the other hand, the flow within the roughness sublayer is very complex and unsteady with a vortex shedding in the shear layer. These vortical structures cannot be generated in the drag-porosity approach since obstacles are not explicitly represented. Identifying the drawbacks in order to improve the accuracy of the drag-porosity approach remain a challenge.

The aim of the present study is to get insights into the main advantages and drawbacks of the drag-porosity approach. To assess the model performance, turbulent statistics from the drag approach are compared with those obtained in the same configuration when obstacles are explicitly accounted for. They will be called DP-LES (drag-porosity LES) and OR-LES (obstacle resolving LES), respectively. The numerical simulation details including the governing equations and the simulation configurations are presented in section 5.2. The comparison results are discussed in section 5.3 and section 5.4. Section 5.5 gives the suggestions for improving the performance of the drag-porosity model.

Numerical simulation details

The governing equations

Considering the governing equations of OR-LES (see section 2.2.1), the filtered Navier-Stocks equations of the DP-LES are written as follow:

∂ u i ∂ x i = 0 (5.1) ∂ u i ∂t + ∂ u i u j ∂ x j = - 1 ρ ∂ p ∂ x i + ∂ ∂ x j (ν ∂ u i ∂ x j -τ i j ) - F i ρ (5.2)
Where the drag force F i is introduced here to model the influence of urban canopy in the drag-porosity approach using equation (5.3), see also section 2.2.1.

F i = ρ 2 C D (z) α f u i | u| (5.3)
Here, C D (z) is the drag force coefficient profile (see sec.5.2.2). α f is the volume frontal density, and in our configuration, α f = 1/3h (equation (1.6)). | u| = u j u j where u j is the velocity component in the j th direction. Hereafter, the ( ..) symbol is omitted during the comparison for clarity.

Drag force coefficient

The drag coefficient coefficient (C D ) used in the drag-porosity approach can be calculated from the OR-LES by the following equation:

C D (z) = 2 ∆p(z)/ρU 2 (z) (5.4)
where ∆p(z) is the mean of the laterally integrated pressure difference between the front and back faces of a cube and U(z) is the spatially-averaged streamwise mean velocity (ISA). In this approach, the viscous drag is omitted. In Leonardi and Castro 2010, the viscous drag was estimated as 7% of the pressure drag proving the large predominance of pressure drag.

The C D profile from the equation provided by [START_REF] Maché | Représentation multi-échelle des transferts entre couche de canopée urbaine et atmosphère à l'échelle de la ville[END_REF] was optimized in order to model as closely as possible the mean velocity profile inside the canopy. The equation for the staggered cube arrays with a 25% plan density is: From the ground up to the height z = 0.4h, rapid decrease C D with the height is observed in all these data. In the range 0.4h < z < 0.6h, constant values are observed in the literature. But different behaviors are observed in the "C D -resolved" and "C D -Mache", where "C D -resolved" decrease with the height and "C D -Mache" slightly increase with the height up to the height z = 0.8h. Then, all data are forced to zero at the top of the canopy. In general, the "C D -resolved" has a larger value at the same vertical altitude than the "C D -Mache".

C D (z) =

Simulation details

The present work is performed in the software OpenFOAM 2.4.0. Standard Smagorinsky SGS model (section 2.2.2) with C S = 0.167 is used in the drag-porosity simulation with "C D -resolved". Dynamic Smagorinsky SGS model (section 2.2.2) is used in the drag-porosity simulation using "C D -Mache". The drag-porosity LES performed with "C D -Mache" and "C Dresolved" are called "DP-LES M " and "DP-LES R ", respectively. Equations are solved using second-order implicit linear differencing for the spatial derivatives and second-order implicit linear backward differencing for the temporal integration as in OR-LES. The computational domain and boundary conditions are equivalent to the setup in the OR-LES modeling (see Sec. 3.2.1). The mesh used here is uniform (∆ = h/16) in all directions from ground to z = 1.5h, and the mesh is expended in the vertical direction (see Fig. 5.2). The roughness Reynolds number R τ based on the friction velocity and obstacle height is about 577. To ensure the flow temporal convergence, simulations were run for an initial duration about 200T . Statistics were collected and averaged over a further duration of 200T for computing flow statistics and turbulent kinetic energy budget. More simulation details are summarized in Table 5.1.

Assessment of the drag-porosity approach

To assess the performance of the drag-porosity approach inside the RSL, the spatial and temporal averaged statistics from the DP-LES R and DP-LES M are compared to OR-LES results. In OR-LES, the spatial averaging of the variable is denoted by an operator (⟨..⟩), which is spatial averaging over an horizontal plane. Both intrinsic and extrinsic spatial averaging (see section 1.4.1) are performed, they are called OR-LES ISA and OR-LES ESA , respectively. For all approaches, ⟨u⟩ h is the streamwise velocity averaged in horizontal direction and time at height z = h, and ⟨u⟩ 4h is the streamwise velocity averaged in horizontal direction and time at z = 4h. Friction velocity u * (see section 1.4.3) is calculated using equation (5.6).

u * = -⟨u ′ w ′ ⟩ h (5.6)

First-order statistic

The spatial and temporal averaged streamwise velocity (⟨u⟩) normalized by ⟨u⟩ 4h resulting from the drag-porosity LES performed with "C D -Mache" and "C D -resolved" are compared to the data of ⟨u⟩ using the intrinsic and extrinsic spatial averaging methods from OR-LES results (Figure 5.3 (a)). The mean velocity profiles of both porosity approaches show a satisfactory agreement to the obstacle-resolved simulation above the height z = h. Slight difference is observed between z = 1.5h and z = 3h, where porosity approaches underestimate the obstacle resolved data. Within the canopy, different flow behavior is observed in porosity approach. The DP-LES M seems to strongly increase the mean velocity and deviate from the obstacle-resolved results, while DP-LES R gives a slightly smaller value than the spatial averaging OR-LES data within the canopy. The reason why ⟨u⟩ from the DP-LES M is higher is due to the reduced drag coefficient that gives an enhanced mean velocity within the canopy. shows the ⟨u⟩ profile normalized by the mean velocity at the top of the cubes ⟨u⟩ h . Both drag-porosity approaches overestimate the ⟨u⟩ profile of OR-LES ESA in the range of 0.8 < z/h < 1. From the ground up to z = 0.8h, DP-LES M fits the OR-LES ISA profile better than the DP-LES R , but DP-LES R is closer to the OR-LES ESA profile than DP-LES M . However, the shape of the mean velocity profile (and in consequence the vertical velocity gradients) inside the canopy obtained with C D -Mache shows the best agreement with both OR-LES data.

In general, the drag-porosity approach models reasonably well the air-flow above the canopy, but the results inside are poor. The parametrized drag force coming from [START_REF] Maché | Représentation multi-échelle des transferts entre couche de canopée urbaine et atmosphère à l'échelle de la ville[END_REF] gives the better fit to the intrinsic spatial averaging data inside the canopy.

Second-order statistics

The spatial and temporal averaged σ u , σ u and σ w , Reynolds shear stress and TKE are shown in Figure 5.4. A discontinuity appears in the intrinsic spatial averaging method at z = h because the average is performed only on fluid regions. The drag-porosity simulations lead to an underestimation of ⟨σ u ⟩ (Figure 5.4 (a)), but agree well at the top of the canopy (0.9 < z/h < 1) and well above the canopy (z/h > 3). For ⟨σ v ⟩ ( In general, the second-order turbulence statistics from the drag-porosity approach shows a closer result to the extrinsic spatial averaging value than the intrinsic spatial averaging results.

Skewness of the velocity

The spatially averaged skewness of the streamwise (⟨Sk u ⟩), spanwise (⟨Sk v ⟩) and vertical (⟨Sk w ⟩) velocity components are shown in Figure 5.5.

The first point to note is that the spatially averaged skewness of the spanwise velocity component (⟨Sk v ⟩) is almost negligible in the domain (see Figure 5.5 (b)). Above the canopy, the streamwise velocity skewness (⟨Sk u ⟩) from DP-LES has a satisfactory agreement with the results from OR-LES (Figure 5.5 (a)). Inside the canopy, there is a good agreement of the (⟨Sk u ⟩) between DP-LES M and the OR-LES ISA data in the upper part of the canopy (0.7 < z/h < 1). But the discrepancy is observed at 0.3 < z/h < 0.7, where DP-LES M overestimates the ⟨Sk u ⟩ values. In terms of ⟨Sk w ⟩, good agreement between the DP-LES M and intrinsic spatial averaging data is observed from z/h = 0.7 to z/h = 1, but underestimated ⟨Sk w ⟩ appears in the drag porosity approach in the lower part of the canopy (0.1 < z/h < 0.7). The underestimated value indicates that the vertical fluctuations is not accurately considered in the DP-LES simulation. Since there are strong upward and downward flows around the cube (see Fig. 3.11), the lack of obstacles may reduce the vertical velocity fluctuation inside the canopy. ). After k max z ′ , the energy of the small scales starts to drop more rapidly in the DP-LES simulation comparing to the wind-tunnel data and OR-LES data. For low normalized wavenumbers, there is much lower energy in the DP-LES simulation for kz ′ < k min z ′ , but in general, the deficit of the energy spectrum is not too significant.

Energy spectrum

Turbulent structure analysis

To assess the performance of the drag-porosity approach, the turbulent structure is discussed in this section. The DP-LES analyzed in this section only refers to the simulation performed with "C D -Mache" (DP-LES M ).

Turbulent kinetic energy budget

The spatial averaged turbulent kinetic energy (TKE) budget of DP-LES writes here from Eqs. 5.1 and 5.2), which are modified from the filtered Navier-Stokes equations of OR-LES:

0 = -⟨ u j ∂ ∂ x j 1 2 u ′ i u ′ i A ⟩ -⟨ u ′ i u ′ j S i j P ⟩ -⟨ 1 2 ∂ ∂ x j u ′ i u ′ i u ′ j T r ⟩ -⟨ 1 ρ ∂ ∂ x i u ′ i p ′ T p ⟩ +⟨ ∂ ∂ x j 2ν u ′ i S ′ i j D ν ⟩ -⟨2ν S ′ i j S ′ i j ε r ⟩ -⟨ ∂ ∂ x j u ′ i τ ′ i j T sgs ⟩ +⟨ τ ′ i j S ′ i j ε sgs ⟩ -⟨ 1 ρ F ′ i u ′ i ε F ⟩ (5.7)
where A represents advection by the mean flow, P is the production by shear, T r is the transport by resolved velocity fluctuation, T p is the transport by pressure fluctuation, D ν is the viscous diffusion, ε r is the resolved dissipation, T sgs is the SGS transport, ε sgs is the SGS dissipation, which represents the transfer of energy from resolved to subgrid scales trough the cut-off, and ε F is the drag force dissipation. All these contributions to the TKE budget are directly computed from simulation results. Each component of the TKE budget, used the same proceduce as in section 4.2, and then are spatially averaged over horizontal planes.

Spatial-averaged TKE budget results

Since the drag-porosity approach give a better fit to the intrinsic spatial averaging method in the skewness analysis, shows that the production and dissipation balance above z = 3h. A rather sharp peak of production is observed at height z = h. The Advection term, of negative sign at z = h, shows a similar sign as at P1 (Fig. 4.15). The pressure transport is, with the space-average, nearly zero at z = h.

After comparison, the TKE budget terms from OR-LES ISA and DP-LES M present similar features: the balance between production and dissipation above z = 3h, and a production larger than the dissipation in the range h < z < 3h. Additionally, the productions from the two approaches are similar above z = 2h. Dissipation which includes ε F in the porosity approach, shows a good agreement with the OR-LES within the canopy. This indicates that the drag porosity approach reproduces well the TKE dissipation within the canopy and the TKE production well above the top of the canopy. However, the drag porosity approach completely underestimates the dissipation from z = h to z = 3h, and the peak value of dissipation and production at the roof level. These differences may be due to the poor reproduction of the shear layer developed on the top of the cube in the DP-LES simulation. According to [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF], the shear layer above the urban canopy can produce and dissipate TKE. Comparing to the OR-LES results, the pressure transport term in the DP-LES simulation is reversed but also negligible. The overall trend of turbulent transport is very similar, except for the underestimated peaks at the top of the canopy and the positive values within the canopy. It can also be noted that in the DP-LES, the advection term is almost non-existent. This agrees with the OR-LES results well above the top of the canopy, but lead to an underestimation of the negative value at the roof level.

Quadrant analysis

The quadrant analysis technique [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF]) considers the frequency of occurrence and contribution of the velocity fluctuations to the shear stress u ′ w ′ . Contributions to shear stress can be separated in four quadrants:

Q1 : u ′ > o, w ′ > 0 Q2 : u ′ < o, w ′ > 0 Q3 : u ′ < o, w ′ < 0 Q4 : u ′ > o, w ′ < 0 (5.8)
where Q2 events indicates the transport of low momentum fluid upwards through an ejection process, and Q4 (is linked to the transport of high-momentum fluid downwards by means of sweep events. as a function of height. Since the enormous computational cost in the simulation, there is no available data for quadrant analysis in the OR-LES simulation.

Q1 and Q3 events are less significant than the Q2 and Q4 events since the former two are rare and contribute little to the shear stress (Figure 5.8). Above the canopy, the contribution of Q3 and Q4 events are in satisfactory agreement with the DNS data (Figure 5.8 (b)). There is less contribution of the sweep events (Q4) to the Reynolds stress compared to the ejection events (Q2) above the canopy, but the opposite is observed inside the canopy. This agrees with the DNS data and the conclusions of [START_REF] Raupach | Conditional statistics of Reynolds stress in rough-wall and smooth-wall turbulent boundary layers[END_REF] that Q4 events contribute most of the shear stress near the rough surface. However, the DP-LES simulation underestimates the number and contribution of Q1 both above and within the canopy, Q2 events are overestimated in both number and contribution plots (Figure 5.8). Besides, the cross over point between the mainly sweep and mainly ejection contribution is at around z = 1.2h in the DP-LES simulation, which is slightly lower than in DNS (at z = 1.25h). It is worth noting that there is more sweep events than the ejection events well above the canopy in the DNS results, but it is not reproduced by the DP-LES simulation. These results are consistent with the finding that in the canopy, fewer but more energetic sweep events dominate the contribution to u

Low-momentum regions

Other large-scale structures in turbulence are the low-momentum regions (LMRs). Then are defined as the regions where the instantaneous flow velocity is below a certain percentage for the local averaged velocity [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF]). Figure 5.9 shows the contours of the streamwise velocity fluctuation u ′ at z = 1.5h and z = 2h above the cubes from DNS results [START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF]). The white regions correspond to areas where u ′ is greater than 20% of the local mean average and the black regions correspond to areas where u ′ is less than 20% of this averaged value. The low velocity streaks (the black regions) are elongated in the streamwise direction. They are more LMRs observed near the canopy (Fig. 5.9 (left)) than above (Fig. 5.9 (right)).

Two-point correlation

The two-point correlation (section 1.4.5) of the streamwise velocity component (R uu ) from the drag-porosity simulation at the reference location z = 1.5h and z = 3h, are presented and compared with the wind-tunnel experiment data from Rivet (2014) in Figure 5.12. In the plot, we set a correlation level equal to 0.7 in the DP-LES data to compare the size of the turbulent structures of the different studies.

Compared to the results of [START_REF] Rivet | Étude en soufflerie atmosphérique des interactions entre canopée urbaine et basse atmosphère par PIV stéréoscopique[END_REF], the R uu spatial correlation maps also highlight the presence of elongated structure, the streamwise dimension L x increases with altitude in agreement with the wind tunnel data.

At z = 1.5h, the turbulent structure in the DP-LES has the maximum streamwise dimension L x = 2h which is smaller than the maximum L x = 3h obtained in the wind tunnel. The same result is found at height z = 3h, where the maximum streamwise dimension L x = 4h in the wind tunnel is much larger than the maximum value in the DP-LES (L x = 3h). This is probably due to the fact that in the work of [START_REF] Rivet | Étude en soufflerie atmosphérique des interactions entre canopée urbaine et basse atmosphère par PIV stéréoscopique[END_REF], the wind tunnel develops a boundary layer δ = 19.5h much larger than the δ = 8h in the DP-LES. The major drawback of the drag-porosity approach is mainly related to the underestimation of the statistics moments inside the canopy. Indeed, it has been shown that even if the mean behavior is correctly predicted and a good agreement is found for the Reynolds shear stress inside the canopy, high-order statistics such as the standard deviation of velocity as well as TKE budget remains underestimated.

As the opinion in [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF], the flow in staggered configuration cubes is diverted laterally because of the staggered cubes in adjacent rows and speeds up in the gap, leading to greater pressure on the front face of the cubes. However, this enhanced pressure effect can not be captured in the drag-porosity approach. Therefore, the model should separately consider the near-wall flow features and this effect should be added into the model.

Spatial averaging methods

The local features is impossible to be captured in the drag-porosity approach, and the only way to compare the turbulence statistics is spatial averaging. The method for generating a spatial averaging value in obstacle resolved LES depends on whether one consider the solid volume. It has been shown that the difference between the two averaging methods is mainly inside the canopy. Although these two methods provide good results for comparing the DP-LES results, it still needs to choose a reasonable way to validate the drag-porosity statistics. Indeed, the intrinsic spatial averaging gives the better results in the mean flow and skewness, but the extrinsic spatial averaging in the second-order statistics is closer to the drag porosity data.

Add drag profile above the canopy

Currently, the C D profile is only used within the urban canopy. According to the results of the TKE budget, the underestimated dissipation above the canopy indicates that the lack of obstacles leads to the lack of a strong shear layer above the canopy, which in turn leads to the underestimation of dissipation and production of TKE. Therefore, the drag coefficient should be added to z = 1.5h above the canopy to simulate the effects by the strong shear layer.

Conclusions

The spatially and temporal averaged mean velocity profiles, second-order statistics profiles and skewness data from the drag-porosity approach are compared with the obstacle-resolving results in this chapter. Good agreement in these statistics is observed above the canopy. This good agreement achieved between the OR-LES and the DP-LES shows the possibility to use drag-porosity approach to generate mean flow characteristics in the urban canopy applications. The DP-LES using the C D profiles provided by [START_REF] Maché | Représentation multi-échelle des transferts entre couche de canopée urbaine et atmosphère à l'échelle de la ville[END_REF] gives a better fit to the intrinsic spatial averaging OR-LES data than using C D profiles from obstacleresolving method. Besides, the energy spectrum is presented in this chapter to demonstrate the performance of the present drag-porosity approach is able to capture the energy transfer in the inertial sub-range.

From the analysis of the spatially averaged TKE budget, the TKE terms differ greatly in obstacle resolved method compared to the drag-porosity approach. Although the TKE budget of the drag-porosity approach leads to a correct representation of production above z = h, it differs from the obstacle-resolved simulation is observed below the top of cubes. As for dissipation, good agreement is observed inside the canopy, but the sharp change at the roof the canopy and the value above are underestimated. Additionally, discrepancies concerning the lack of the reversed sign of pressure transport and the almost null advection need to be further considered.

In the coherent structure analysis, underestimated Q2 and overestimated Q1 and Q3 events are observed inside the canopy. The elongated streamwise LMRs is clearly shown in the drag porosity approach, but much thicker in the spanwise direction at the roof level. Finally, a satisfactory agreement of the streamwise velocity fluctuation correlation R uu from DP-LES was observed compared to the wind tunnel data.

Finally, there are still certain discrepancies observed inside the canopy. This discrepancies indicate that the understanding about the drag force approach is still not enough. Therefore, a more in-depth study of the characteristic of the drag porosity approach such as the proper orthogonal decomposition (POD) is needed to provide more suggestion for improving the performance of this model, especially within the canopy.

LES validation

The achievement of the present LES model in the field of the large-scale urban canopy opened the perspectives of improving the numerical applications. In the current work, LES using the dynamic Smagorinsky SGS model were performed with a refined mesh able to resolve the flow in the viscous sublayer. Such fine meshes are rarely reported in literature.

The LES code is firstly assessed in comparison with the data from DNS [START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]) and wind-tunnel [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF], Blackman and Perret 2016[START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF][START_REF] Herpin | Investigation of the flow inside an urban canopy immersed into an atmospheric boundary layer using laser Doppler anemometry[END_REF]). The investigation on the mean flow illustrates the inhomogeneous characteristic of the flow inside the urban boundary layer. The maximum simulated Reynolds shear stress is slightly above the top of the cube, but in magnitude smaller than the experimental data. Standard deviations of the streamwise and vertical velocity component match well with the wind-tunnel data above the canopy. The standard deviation of the span-wise velocity component is in rather good agreement with the wind-tunnel data from [START_REF] Herpin | Investigation of the flow inside an urban canopy immersed into an atmospheric boundary layer using laser Doppler anemometry[END_REF] above and inside the canopy top. However, it is significantly smaller than the experimental data of [START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF], but larger than the experimental data of [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]. This is mainly explained by the measurements and simulation resolution difference or the difference in h/δ ratio. The skewness shows that sweep events mainly occur inside the canopy, and ejection events contribute more above the canopy. The energy spectra confirm the existence of inertial subrange and the limitations by the computational domain and mesh resolution. On the other hand, the Q-criterion analysis confirmed that the present mesh resolution is fine enough to capture the vortex structures of the flow. This finding highlights in particular the need to consider the three-dimensionality of the flow.

TKE budgets and its 3D structure

Each term of the local TKE budget is computed in this work allowing to describe the complex behavior of turbulent flow. A peak of production is observed in the wake of the cube where the shear layer induced by the presence of the cubic roughness is the strongest. The peak decreases downwind as the shear layer increases. The production and dissipation are major contributions to the TKE budget above z/h > 2, and are nearly in equilibrium. The turbulence transport acts as a sink in the shear layer, but a source above and below. The pressure transport is significant in front the cube. Advection is a source of energy in front of the cube, but a sink in the wake. Production shows good agreement with the wind-tunnel results both above and within the canopy at location P3. The 3D iso-surface visualization of the TKE budget reflects the spatial variation of TKE terms. High production mainly concentrates on the edges of the cubes and in the shear layer. This high production is transferred by the turbulent transport upward and downward. A strong pressure transport acts as a source in front of the cube where the incoming flow is impinging on the cubes. Advection acts as a sink near the high production region. The negative advection in the wake can be explained by the positive streamwise gradient of TKE and positive streamwise velocity. Dissipation is mainly in the near-wall regions. High dissipation is also observed the shear layer.

Finally, the feasibility of ignoring the span-wise gradients when calculating turbulent transport in experiment [START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]) is verified with the present simulation results. The comparison of the results at location P3 highlights the contribution of span-wise terms inside the canopy. Underestimated turbulent transport is observed inside the canopy and shows a good agreement with the wind-tunnel data. 3D structures provide more insight into the spatial distribution of discrepancies, the largest are observed in the spanwise gap between two cubes cubes.

Drag porosity approach in LES

The drag-porosity approach is performed with LES to verify the capabilities of the proposed method. Satisfactory results were observed in the spatial and temporal averaged flow velocity. The good agreement with the obstacle resolved LES shows the ability of the drag-porosity approach to produce reasonable flow characteristics above the RSL that can be used to generate the unsteady flow. However, some discrepancies are observed inside the canopy. Besides, the TKE budgets differ greatly between the drag-porosity approach and obstacle resolving method. The TKE budget of the drag-porosity approach leads to a correct representation above z = 2h but differs from the obstacle-resolved simulation below. The discrepancies concern large lack of dissipation, the reversed sign of pressure transport and the almost null advection. Therefore, developments are still needed to improve the performance of the drag-porosity approach to capturing more detailed information, especially inside the canopy.

Perspectives

The remarkable performance of the present model made it interesting to continue the research over urban-like canopy. Further research will be conducted using the present LES model with the dynamic Smagorinsky SGS model to investigate the effects of configurations of buildings with more complex geometries, such as various building heights or different packing density.

Recently, the experimental data have shownthat large-scale coherent structures interact and influence the small-scale structures in the shear layer through a non-linear interaction [START_REF] Blackman | Non-linear interactions in a boundary layer developing over an array of cubes using stochastic estimation[END_REF]. Future work could combine with wind-tunnel data to study the nonlinear transfer relationships between different scale structures. This could help quantify the uncertainty of the present model. It is also worthwhile to improve the DP-LES simulation. A more advanced process that can help to eliminate the differences between the OR-LES simulations and the DP-LES simulations should be analyzed.

  High rise and low rise in the downtown Beijing, China; (b) Relatively uniform height of buildings in the downtown Paris, France (Photos are extracted from Google website). . . . . . . . . . . . . . . . . . . . . . . . 2 Schematic roadmap for numerical research work involving urban canopy flow studies. Blue plots: the work achieved, the green plots: data from literature and hatched yellow blocks: work that needs to be done in the future. 1.1 Schematic of the atmospheric boundary layer structure including five different layers in the vertical direction. . . . . . . . . . . . . . . . . . . . . . . 1.2 Schematic plot for calculating the morphological parameters in urban canopy, from Maché (2012). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Schematic plot of flow regimes classification over urban-like building arrays. (a) Isolated roughness flow W /H > 3.33, (b) Wake interference flow 1.53 < W /H < 3.33, (c) Skimming flow W /H < 1.53 (Figure from Oke et al. 2017). 1.4 Schematic plan view of the cubic arrays layouts with a plan area density λ p = 25% of there types: (a) aligned arrays, (b) squared arrays, (c) staggered arrays. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Example of the obstacle resolving method where the block within the grid is used to represent the urban canopy. . . . . . . . . . . . . . . . . . . . . . 1.6 Schematic plot of urban canopy created using immersed boundary method (Figure extracted from Giometto et al. (2016). . . . . . . . . . . . . . . . . 1.7 Spectrum of the axial turbulence component (figure extracted from Castro et al. 2006). Note that z ′ = zd and k = 2π f /u, with f is the frequency. . . . 2.1 Schematic of spectrum of turbulent kinetic energy and the energy cascade. The pink line marks the delimiting line between the resolved scale and the subgrid scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Instantaneous view of the C coefficient of the dynamic Smagorinsky model in a vertical symmetry plane. . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 2D Schematic plot of volume surrounding the node point on a mesh in finite volume method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Schematic plot of the no-slip condition. . . . . . . . . . . . . . . . . . . . 2.5 Schematic plot of the free-slip condition. . . . . . . . . . . . . . . . . . . . 2.6 Schematic 3D view of the computational domain [4h × 4h × 4h] occupied by the staggered cube array in the preliminary simulations, where h is the height of cubes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Schematic plan view of the computational domain in staggered cube arrays. 2.8 Vertical profiles of mean streamwise velocity at locations P1 (a) and P2 (b) and P3 (c) and P0 (d). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Vertical profiles of velocity standard deviation behind a cube at P1 (a, c) and in front of a cube at P2 (b, d). . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Vertical profiles of Reynolds shear stress behind a cube at P1 (a) and in front of a cube at P2 (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 View of the mesh generated by the OpenFOAM tools: blockMesh and Snap-pyHexMesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Comparison of the time-averaged streamwise velocity, normalized by the friction velocity (u τ ) at locations P1 (a) and P2 (b) and P3 (c) and P0 (d) for the three mesh resolutions. White circles: ∆ = h/64. Blue dashed line: ∆ = h/16. Red solid line: ∆ = h/32. . . . . . . . . . . . . . . . . . . . . . 2.

3. 1 xof figures 4 . 1

 141 Perspective view of the computational domain, where h is the cube height. . 3.2 Plan view of the refined mesh around the cubes in the present simulation. . 3.3 Vertical profiles of mean streamwise velocity at locations P1 (a), P2 (b), P3 (c), P0 (d) indicated in figure. 3.1. Blue solid line: LES computations. Red dashed line: DNS data from Coceal et al. (2007b). Circles: wind-tunnel data from Castro et al. (2006). Squares: wind-tunnel data from Herpin et al. (2018). Stars: PIV data from Blackman et al. (2017) . . . . . . . . . . . . 3.4 Vertical profile of normalized Reynolds shear stress and vertical velocity standard deviation at P1 ((a), (b)), P2 ((c), (d)), and P3 ((e), (f)) . Solid line: LES; Dashed line: DNS from Coceal et al. (2007b); Circles: wind-tunnel data from Castro et al. (2006); Triangles: wind-tunnel data from Blackman and Perret (2016); Stars: wind-tunnel data from Blackman et al. (2017). . . 3.5 Horizontal plane of Reynolds shear stress at z = 0.06h. . . . . . . . . . . . 3.6 Vertical profiles of standard deviation of streamwise and spanwise velocity components at P1 ((a), (b)), P2 ((c), (d)), and P3 ((e), (f)). Solid line: LES; Dashed line: DNS from Coceal et al. (2007b); Circles: wind-tunnel data from Castro et al. (2006). Squares: wind-tunnel data from Herpin et al. (2018); Triangles: wind-tunnel data from Blackman and Perret (2016); Stars: wind-tunnel data from Blackman et al. (2017). . . . . . . . . . . . . . . . . 3.7 Horizontal plane of the standard deviation of streamwise velocity components at z = 0.06h. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Vertical profiles of skewness at location P1 (a), P2 (b), P3 (c), P0 (d). Blue solid line: skewness of streamwise velocity component (Sk u ). Red-dashed line: skewness of spanwise velocity component (Sk v ). Green-dashed line: skewness of vertical velocity component (Sk w ). . . . . . . . . . . . . . . . 3.9 Vertical profiles of mean streamwise velocity plotted on a log scale, where ⟨u⟩ is the horizontally and temporally averaged streamwise velocity above the cubes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Turbulent spectra of the streamwise velocity component E(kz ′ ) at P2 from LES (Solid blue line) compared with wind-tunnel data (circles) from Castro et al. (2006). The black dash-line represents the -5/3 slope. The red dash-line and magenta dot-dash-line indicate the minimum and maximum wavenumber resolved in the LES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Mean flow structure in a vertical x-z plane through the middle of the cube in the staggered array showing (u, w) wind vectors. . . . . . . . . . . . . . . . 3.12 Top view of the horizontal transect (x,y) near the ground (z = h/40) with planar streamlines. Wind goes from left to right. . . . . . . . . . . . . . . . 3.13 3D view of the streamlines emanating from vertical lines at P1 (a), P2 (b) and P3 (c). Streamlines color denotes the mean velocity. Green lines indicated the streamlines sources. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 magVorticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Q-criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List Vertical profiles of TKE budget terms (a) behind cube (P1), (b) in front of cube (P2), (c) in the gap between two cubes (P3) and (d) on cube (P0), all terms are normalized by u 3 τ /h. . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Vertical profiles of |P|/|ε sgs + ε r | ratio (a) behind cube P1, (b) in front of cube P2, (c) in gap P3 and (d) on cube P0. . . . . . . . . . . . . . . . . . . 4.3 Production and dissipation terms of TKE budget compared to the wind-tunnel data from Castro et al. (2006) and Blackman et al. (2017) at location P2 (a) and P3 (b). All terms are normalized by u 3 τ /h. . . . . . . . . . . . . . . . . 4.4 (a) Decomposed production terms from LES at P3, all terms are normalized by u 3 τ /h. (b) Resolved strain rate tensor (S i j ) at P3, normalized by u τ /h. . . 4.5 Three-dimensional visualization of Production, normalized by u 3 τ /h. For clarity, only three cubes are shown. (a) Zoom inside the urban canopy with iso-contours of P = 10 (orange) and P = -10 (light blue). Iso-contours are shown only below z = 3h/4. (b) View of iso-contours of Production P = 25 (red) and P = -25 (dark blue). . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Top view of the horizontal transect (x,y) of Production normalized by u 3 τ /h. Wind goes from left to right. (a) Near the ground (z = h/40) with planar velocity streamlines where S, N and F locations stands for saddle point, node and focus node, respectively. (b) In the middle of the cube (z = h/2) with mean velocity vector field. Vectors are of equal length to better visualize flow structures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 (a) Vertical transect of Production normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures. (b) Production in the windward face of a cube with wall streamlines. . . . . . . . . . . . . . . . . . . . . . 4.8 Vertical transect of Turbulent transport normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures. The pink contour circles is Production (P = 25). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Three-dimensional visualization of Turbulent transport (T r ), normalized by u 3 τ /h. For clarity, only three cubes are shown. (a) View of T r = 5 (pink). (b) View of iso-contours of Turbulent transport T r = -8 (dark blue). . . . . . . 4.10 Top view of the horizontal transect (x,y) of Turbulent transport (T r ), normalized by u 3 τ /h, in the middle of the cube (z = h/2) with mean velocity vector field. Vectors are of equal length to better visualize flow structures. Wind goes from left to right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Three-dimensional visualization of Pressure transport, normalized by u 3 τ /h. For clarity, only three cubes are shown. (a) View of iso-contours of T p = 15 (pink). (b) View of iso-contours of Pressure transport T p = -5 (dark blue). . 4.12 Vertical transect of Pressure transport (T p ) normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures. . . . . . . . . . . . . . . . . 4.13 Three-dimensional visualization of Advection, normalized by u 3 τ /h. For clarity, only three cubes are shown. (a) View of iso-contours of Advection A = 5 (pink). (b) View of iso-contours of Production P = 25 (red) and of Advection A = -5 (dark blue). . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Vertical transect of Advection (A) normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures. The pink contour circles is Production (P = 25). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Vertical transect of TKE normalized by u 2 τ in the middle of the cube with a projection of the mean velocity vector field. Vectors are following the value of time-averaged velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Top view of the horizontal transect (x,y) of Advection normalized by u 3 τ /h, in the middle of the cube (z = h/2) with mean velocity vector field. Vectors are of equal length to better visualize flow structures. Wind goes from left to right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Three-dimensional visualization of iso-contours of Dissipation ε r + ε sgs < -5, normalized by u 3 τ /h. For clarity, only three cubes are shown. . . . . . . 4.18 Vertical transect of dissipation (ε r + ε sgs ) normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures. . . . . . . . . . . . . . . . . 4.19 (a) -(f) The decomposition of Turbulence transport (T r ) at location P3 into the contribution of each velocity gradient compared to the wind-tunnel data from Blackman et al. (2017). All terms normalized by u 3 τ /h. . . . . . . . . 4.20 (a) Turbulent transport (T r ) compared to the wind-tunnel data from Blackman et al. (2017) and T r without three span-wise gradients. (b) Three span-wise gradients of Turbulent transport (T r ). All terms normalized by u 3 τ /h. . . . . 4.21 3D visualization of the sum of the spanwise gradients of the turbulent transport terms discrepancy between the Turbulence transport and Turbulence transport without span-wise velocity fluctuation gradients. Red color means ∆T r × h/u 3 τ = 5 and blue color means ∆T r × h/u 3 τ = -5. . . . . . . . . . . . 5.1 Vertical profile of drag coefficient C D using the equation from Maché (2012) and extracted from obstacle resolving method, compared to the data from Cheng and Castro (2002) and Coceal et al. (2006). . . . . . . . . . . . . . . 105 5.2 Vertical view of the mesh generated by the OpenFOAM tool: blockMesh. . 107 5.3 Vertical profiles of (a) spatial and temporal streamwise mean velocity ⟨u⟩, normalized by ⟨u⟩ 4h : comparison between drag-porosity LES and obstacleresolved LES using intrinsic and extrinsic method. (b) spatial and temporal streamwise mean velocity ⟨u⟩ within the canopy, normalized by ⟨u⟩ h . Porosity C D -resolved (red solid line), porosity C D -Maché (blue dashed line), obstacle-resolved LES using intrinsic spatial averaging (white circle) and extrinsic spatial averaging (black dot). . . . . . . . . . . . . . . . . . . . . 109 5.4 Vertical profiles of the spatial averaged (a) standard deviation of the streamwise velocity, ⟨σ u ⟩; (b) standard deviation of the spanwise velocity ⟨σ v ⟩; (c) standard deviation of the vertical velocity ⟨σ w ⟩; (d) Reynolds shear stress, normalized by u * ; (e) turbulent kinetic energy (TKE), normalized by u 2 * . . . 110 5.5 Vertical profiles of spatial and temporal averaged (a) skewness of streamwise velocity component (⟨Sk u ⟩); (b) skewness of spanwise velocity component (⟨Sk v ⟩); (c) skewness of vertical velocity component (⟨Sk w ⟩). . . . . . . . . 112 5.6 Turbulent spectra of the streamwise velocity component E(kz ′ ) at P2 and z/h = 1.62 from OR-LES (Solid blue line) and from DP-LES (Solid cyan line) compared with wind-tunnel data (circles) from Castro et al. (2006). The black dash-line represents the -5/3 slope. The red dash-line and magenta dot-dash-line indicate the minimum and maximum wavenumber resolved in the LES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.7 Vertical profiles of spatial averaged TKE budget terms normalized by u 3 τ /h, (a) OR-LES ISA and (b) DP-LES M . . . . . . . . . . . . . . . . . . . . . . . 114 5.8 (a) Relative number of events in each quadrant as a function of height. (b) Relative contribution to u ′ w ′ of events in each quadrant as a function of height. Blue: Q1; Red: Q2; Green: Q3; Yellow: Q4. Points: DNS from Coceal et al. (2007b); dashed line: DP-LES data. . . . . . . . . . . . . . . 116 5.9 Instantaneous snapshot showing contour plots of streamwise velocity at two altitudes: z = 1.5h (left) and z = 2h (right) from Coceal et al. (2007b). Black regions: u < 0.8u. White regions: u > 1.2u. Here u is the local mean velocity.117 xiii 5.10 Instantaneous snapshot showing contour plots of streamwise velocity at different altitudes: a) z = 1.5h, b) z = 2h, c) z = 3h and d) z = 4h of the OR-LES result. Blue regions: u < 0.8u. Red regions: u > 1.2u. Here u is the local mean velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 Instantaneous snapshot showing contour plots of streamwise velocity at different altitudes: a) z = 1.5h, b) z = 2h, c) z = 3h and d) z = 4h of the DP-LES result. Blue regions: u < 0.8u. Red regions: u > 1.2u. Here u is the local mean velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 Two-point correlation plot in the XY plot from the wind-tunnel experiment from Rivet (2014). (a) and (c) from the DP-LES simulation; (b) and (d) from Rivet (2014). (a) and (b) at z = 1.5h; (c) and (d) at z = 3h. . . . . . . . . . A.1 Flow chart of the PIMPLE algorithm . . . . . . . . . . . . . . . . . . . . .

Figure 1

 1 Figure 1 (a) High rise and low rise in the downtown Beijing, China; (b) Relatively uniform height of buildings in the downtown Paris, France (Photos are extracted from Google website).

Figure 2

 2 Figure 2 Schematic roadmap for numerical research work involving urban canopy flow studies. Blue plots: the work achieved, the green plots: data from literature and hatched yellow blocks: work that needs to be done in the future.

Figure 1 . 1

 11 Figure 1.1 Schematic of the atmospheric boundary layer structure including five different layers in the vertical direction.

Figure 1 . 2

 12 Figure 1.2 Schematic plot for calculating the morphological parameters in urban canopy, from Maché (2012).

Figure 1 . 3

 13 Figure 1.3 Schematic plot of flow regimes classification over urban-like building arrays. (a) Isolated roughness flow W /H > 3.33, (b) Wake interference flow 1.53 < W /H < 3.33, (c) Skimming flow W /H < 1.53 (Figure from Oke et al. 2017).

Figure 1 . 4

 14 Figure 1.4 Schematic plan view of the cubic arrays layouts with a plan area density λ p = 25% of there types: (a) aligned arrays, (b) squared arrays, (c) staggered arrays.

Figure 1 . 5

 15 Figure 1.5 Example of the obstacle resolving method where the block within the grid is used to represent the urban canopy.

Figure 1 . 6

 16 Figure 1.6 Schematic plot of urban canopy created using immersed boundary method (Figure extracted from Giometto et al. (2016).

  Iaccarino 2005, Cheng and Porte-Agel 2015, Giometto et al. 2016). Giometto et al. (2016) used a discrete forcing immersed boundary method for representing the urban canopy (see Figure 1.6).

Figure 1 . 7

 17 Figure 1.7 Spectrum of the axial turbulence component (figure extracted from Castro et al. 2006). Note that z ′ = zd and k = 2π f /u, with f is the frequency.

  al. 2007; Inagaki and Kanda 2008); the "URBAN 2000" project (Allwine et al. 2002) and "Joint URBAN 2003" project (Brown et al. 2004; Klein and Clark 2007) in the United States of America, etc. In particular, environmental factors were considered in the study, such as Vachon et al. (2002),who investigated the impact of traffic on the development of flow field in street canopy.

  (2009) performed LES to understand the effect of the complexity of urban canopy representation on the flow. To close the equation, a Lagrangian[START_REF] Meneveau | A Lagrangian dynamic subgridscale model of turbulence[END_REF]) dynamic scale-dependent SGS model (Bou-Zeid et al. 2005) is applied. A large negative shear production was found inside the canopy at the transition between the surrounding area and the neighborhood of interest indicating that the energy is transferred from the turbulent motions to the mean flow.[START_REF] Kono | Numerical Investigations of Mean Winds Within Canopies of Regularly Arrayed Cubical Buildings Under Neutral Stability Conditions[END_REF] performed LES simulations with a dynamic SGS model over the aligned and staggered array in a range of packing densities. The results were validated against the wind tunnel data from[START_REF] Cheng | Near wall flow over urban-like roughness[END_REF]. They found that the drag coefficient is independent of the packing densities for the aligned array, but it was very sensitive to the packing density in the staggered array. Cheng and Porte-Agel (2015) performed LES with a modulated gradient SGS model (Cheng and Porte-Agel 2013) over staggered and aligned array in different packing densities. The study found that the displacement height (d) is insensitive to the canopy arrangement (staggered or aligned) and the roughness length (z 0 ) is larger in the staggered array than aligned array.[START_REF] Claus | Wind-Direction Effects on Urban-Type Flows[END_REF] using LES analyzed the effect of a set of upstream wind directions on staggered cubic array in a [8h × 8h × 4h] computational domain. The results demonstrated that the spatially averaged mean velocity in the staggered array is very dependent on the oncoming flow direction.

  2.2.1 Filtered Navier-Stokes equationsLES only resolves the large-scale structures in the flow by filtering the Navier-Stokes equations (see Figure2.1), and uses a subgrid-scale (SGS) model to simulate the effects of the unresolved scales on the filtered scales.

Figure 2 . 1

 21 Figure 2.1 Schematic of spectrum of turbulent kinetic energy and the energy cascade. The pink line marks the delimiting line between the resolved scale and the subgrid scale.

  9) where ν is the kinetic viscosity, C k and C ε are model coefficients. Schumann (1975) used C k = 0.094 and C ε = 1.02; Menon et al. (1996) used C k = 0.05 and C ε = 1.0. In OpenFOAM, the default values are those used in Yoshizawa (1986): C k = 0.094 and C ε = 1.048 .

  Germano et al. (1991) introduced a dynamic procedure for evaluating subgrid-scale model coefficients using information contained in the resolved velocity field. The aim is to adapt the coefficient value automatically in both space and time (Zang et al. 1993, Calmet and Magnaudet 1997, Giometto et al. 2016). It has two main advantages compared to the traditional eddy-viscosity SGS model. (a) Without using damping function or intermittency function, the model coefficient is neither a priori nor a fixed-constant value and varies according to the simulation process locally; (b) The calculation of the model coefficient could produce negative values, which eliminates the drawback of the conventional eddyviscosity SGS model that prevent energy backscatter from small scale to resolved scale (Piomelli et al. 1991, Yakhot et al. 2006, He et al. 2017).

Figure 2 . 2

 22 Figure 2.2 Instantaneous view of the C coefficient of the dynamic Smagorinsky model in a vertical symmetry plane.

  2

  .

Figure 2 . 3 2D

 23 Figure 2.3 2D Schematic plot of volume surrounding the node point on a mesh in finite volume method.

3 .

 3 Solve implicitly the momentum equation to obtain a new velocity u * 4. Compute the mass fluxes at faces ṁ( * ) f using the Rhie-Chow interpolation (Rhie and Chow 1983).

5 .

 5 Update the pressure and velocity fields at the cell centroids and the mass flow rate at the cell faces p ( * ) , u * * and ṁ( * * )f to obtain continuity-satisfying fields.

Figure 2 . 4 Figure 2 . 5

 2425 Figure 2.4 Schematic plot of the no-slip condition.

Figure 2 .Figure 2 . 6 Figure 2 . 7

 22627 Figure 2.6 shows the computational domain of the preliminary simulations with length (L) × width (W) × height (H) = [4h × 4h × 4h], where h = 0.02[m] is the cube height. The urban canopy is represented by a staggered array of cubic obstacles with constant height and packing density λ p = 25%. The size of the computational domain has been chosen in reference to Coceal et al. (2006) and Xie and Castro. (2006), who performed DNS and LES separately in the same domain size and obtained a good agreement of the turbulent statistics with the wind-tunnel experiment data from[START_REF] Cheng | Near wall flow over urban-like roughness[END_REF] and[START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF]. Therefore, this computational domain size is considered in the preliminary simulation for the selection of LES parameters.

Figure 2 .Figure 2 . 8

 228 Figure 2.8 presents the mean streamwise velocity normalized by the friction velocity u * at positions P0-P3, compared with DNS data from Coceal et al. (2007b) and wind-tunnel measurement data from Cheng and Castro (2002). The friction velocity u * is calculated by using the constant Reynolds shear stress just above the crest of the canopy. The comparison of the profiles show that the LES results are in good agreement with the literature data, which indicates that the present LES model gives a reliable result for the urban canopy study. It is found little difference between the velocity from various SGS models, even if the data from dynamic Smagorinsky model is closer to DNS data near the ground (Figure 2.8 (a)).

Figure 2 .

 2 Figure2.9 presents the comparisons of the standard deviation of velocity normalized by the friction velocity u * at two positions P1 (a,c) and P2 (b,d). Figure2.9 (a) shows that all the numerical simulation values underestimate the measured standard deviation of streamwise velocity σ u within the canopy, but the dynamic Smagorinsky SGS model results are in better agreement with DNS than the other non-dynamic SGS models. The differences between the SGS models at P2 (Figure2.9 (b)) are not very clear, but agreement between DNS and LES data is very good from 0.1 < z/h < 0.3. These similar performances can also be seen in standard deviation of vertical velocity component σ w : the dynamic Smagorinsky model provides the closest results to the measured data in Figure2.9 (c). Figure2.9 (d)) shows that the maximum difference of σ w between LES and wind-tunnel data is up to about 20 % at the

Figure 2 .

 2 Figure2.9 presents the comparisons of the standard deviation of velocity normalized by the friction velocity u * at two positions P1 (a,c) and P2 (b,d). Figure2.9 (a) shows that all the numerical simulation values underestimate the measured standard deviation of streamwise velocity σ u within the canopy, but the dynamic Smagorinsky SGS model results are in better agreement with DNS than the other non-dynamic SGS models. The differences between the SGS models at P2 (Figure2.9 (b)) are not very clear, but agreement between DNS and LES data is very good from 0.1 < z/h < 0.3. These similar performances can also be seen in standard deviation of vertical velocity component σ w : the dynamic Smagorinsky model provides the closest results to the measured data in Figure2.9 (c). Figure2.9 (d)) shows that the maximum difference of σ w between LES and wind-tunnel data is up to about 20 % at the

Figure 2 .Figure 2 . 9

 229 Figure2.9 presents the comparisons of the standard deviation of velocity normalized by the friction velocity u * at two positions P1 (a,c) and P2 (b,d). Figure2.9 (a) shows that all the numerical simulation values underestimate the measured standard deviation of streamwise velocity σ u within the canopy, but the dynamic Smagorinsky SGS model results are in better agreement with DNS than the other non-dynamic SGS models. The differences between the SGS models at P2 (Figure2.9 (b)) are not very clear, but agreement between DNS and LES data is very good from 0.1 < z/h < 0.3. These similar performances can also be seen in standard deviation of vertical velocity component σ w : the dynamic Smagorinsky model provides the closest results to the measured data in Figure2.9 (c). Figure2.9 (d)) shows that the maximum difference of σ w between LES and wind-tunnel data is up to about 20 % at the

Figures 2 .

 2 Figures 2.10(a, b) present the Reynolds shear stress (-u ′ w ′ ) profiles normalized by the friction velocity at two positions P1 (a) and P2 (b), using data from LES, DNS and measurements. The results from LES are fitting very well with the DNS data, especially in the wake of the cubes (Figures 2.10(a)). But it underestimates the measurement values at the top of canopy.Figure 2.10 (b)shows that there is a peak in DNS at the height of z = 0.1h just above the ground surface, which is overestimated by the LES data. Among these SGS models, dynamic SGS model gives the best simulation value in comparison with the DNS data.

  (a)). But it underestimates the measurement values at the top of canopy.

Figure 2 .P2Figure 2 .

 22 Figure 2.10 Vertical profiles of Reynolds shear stress behind a cube at P1 (a) and in front of a cube at P2 (b).

Figure 2 .

 2 Figure 2.11 View of the mesh generated by the OpenFOAM tools: blockMesh and Snappy-HexMesh.

Figure 2 .

 2 Figure 2.12 Comparison of the time-averaged streamwise velocity, normalized by the friction velocity (u τ ) at locations P1 (a) and P2 (b) and P3 (c) and P0 (d) for the three mesh resolutions. White circles: ∆ = h/64. Blue dashed line: ∆ = h/16. Red solid line: ∆ = h/32.

Figure 2 .

 2 Figure 2.14 Plan view of the refined mesh around the cubes.

P1Figure 2 .

 2 Figure 2.15 The contour plot of the two-point correlation of streamwise velocity Ruu, where the reference position is in the wake of the cube at a height of 1.5h.

Figure 2 .

 2 Figure 2.15 depicts a contour plot of the two-point correlation of streamwise velocity (section 1.4.5), where the reference position is in the wake of the cube (P1) at the height of 1.5h. It clearly shows that the large-scale turbulent structure above the cube is blocked by the finite size of the computational domain, although applying periodic boundary conditions here. This figure shows that a larger computational domain is needed to provide complete turbulence structure information. In reference to[START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF], the computational domain [16h × 12h × 8h] is chosen in the next chapter. The statistics will be validated against the literature data, and the 3D flow characteristics will be revealed in the analysis.

Figure 3 .Figure 3 . 1

 331 Figure 3.1 shows the computational domain of the present simulation with dimension [16h × 12h × 8h] in streamwise (x), span-wise (y) and vertical (z) directions, respectively. The size of the computational domain has been chosen in reference to Coceal et al. (2007b) who successfully performed DNS of the flow over the urban-like canopy. The floor is covered by staggered cubes array of uniform height (h) and packing density λ p = 25%. The four locations (P0-P3) will be used to compare the mean flow and high-order statistics with the literature. Besides, free-slip condition is set at the top boundary of the domain, and no-slip conditions are prescribed at the floor and on all obstacle surfaces. Periodic boundary

Figure 3 . 2

 32 Figure 3.2 Plan view of the refined mesh around the cubes in the present simulation.

  However,[START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF] observed that the dispersive stress was overestimated when the time used for statistics was only 50T and finally used 400T for the statistics collection.[START_REF] Xie | LES and RANS for turbulent flow over arrays of wall-mounted obstacles[END_REF] argued that 80T is enough for statistics gathering and even observed that the variation of the statistical data throughout the roughness sublayer was usually small after averaging duration longer than 20T . Coceal et al. (2007b and 2007c) chose 100T for each part of their simulation in the same computational domain as in the present work. In recent studies, Tomas et al. (2016) and Castro et al. (2017) used 800T and 710T for the statistics collection, respectively. However, in the present study the statistics obtained over a duration of 200T are converged at least up to the third order moments (see Section 3.3.3).

  The experimental work of[START_REF] Castro | Turbulence Over Urban-type Roughness: Deductions from Wind-tunnel Measurements[END_REF] is used for comparison and enriched by recent experimental datasets from the atmospheric wind tunnel of LHEEA (Nantes, France): the Laser Doppler Velocimetry (LDV) dataset from Herpin et al. (2018) and two datasets using Particle Image Velocimetry (PIV) by[START_REF] Blackman | Non-linear interactions in a boundary layer developing over an array of cubes using stochastic estimation[END_REF] and by[START_REF] Blackman | Turbulent kinetic energy budget in the boundary layer developing over an urban-like rough wall using PIV[END_REF]. Note that in order to be directly compared with numerical results, PIV datasets have two modifications compared to the one reported in the mentioned papers: first vertical profiles are here extracted at points P1, P2 and P3 without spatial averaging and, second, experimental data are here normalized by the friction velocity obtained from drag force measurements (see Section 1.4.3) as in[START_REF] Herpin | Investigation of the flow inside an urban canopy immersed into an atmospheric boundary layer using laser Doppler anemometry[END_REF], instead of from the constant shear layer as used in the original articles. Results are also compared to the direct numerical simulations (DNS) of[START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF].

Figure 3 .Figure 3 . 3

 333 Figure 3.3 shows vertical profiles of the mean streamwise velocity component normalized by the friction velocity (u τ ) at four locations around a cube. LES results (in blue lines) are in good agreement with DNS data from Coceal et al. (2007b) and from the wind-tunnel measurements of Castro et al. (2006) at the four locations presented. This observation is equally valid inside the canopy and up to z/h = 4.

Figures 3 .

 3 Figures 3.4 and 3.6 show turbulence characteristics of the flow normalized by u τ at locations P1, P2 and P3.

FigFigure 3 . 4

 34 Figure 3.4 Vertical profile of normalized Reynolds shear stress and vertical velocity standard deviation at P1 ((a), (b)), P2 ((c), (d)), and P3 ((e), (f)) . Solid line: LES; Dashed line: DNS from Coceal et al. (2007b); Circles: wind-tunnel data from Castro et al. (2006); Triangles: wind-tunnel data from Blackman and Perret (2016); Stars: wind-tunnel data from Blackman et al. (2017).

Figure 3 . 5

 35 Figure 3.5 Horizontal plane of Reynolds shear stress at z = 0.06h.

  3.4 where datasets have very different values of h/δ (see section 3.3.1). The second reason is linked to the vertical resolution of the datasets as a coarser resolution will have the tendency to smooth peak value. LDV measurements by Castro et al. (2017) have a vertical resolution of approximately 0.015h while the PIV resolution by Blackman et al. (2017) is approximately 0.038h. PIV is therefore expected to smooth local peaks of shear components. In LES, the resolution near z = h is approximately 0.031h, which is finer compared to PIV but still a factor of two bigger compared to LDV measurements. This explanation is correlated to observations depicted in Figs 3.4 (a) (b).

Figure 3 . 6

 36 Figure 3.6 Vertical profiles of standard deviation of streamwise and spanwise velocity components at P1 ((a), (b)), P2 ((c), (d)), and P3 ((e), (f)). Solid line: LES; Dashed line: DNS from Coceal et al. (2007b); Circles: wind-tunnel data from Castro et al. (2006). Squares: wind-tunnel data from Herpin et al. (2018); Triangles: wind-tunnel data from Blackman and Perret (2016); Stars: wind-tunnel data from Blackman et al. (2017).

Figure 3 . 7

 37 Figure 3.7 Horizontal plane of the standard deviation of streamwise velocity components at z = 0.06h.

Figure 3 . 8

 38 Figure 3.8 Vertical profiles of skewness at location P1 (a), P2 (b), P3 (c), P0 (d). Blue solid line: skewness of streamwise velocity component (Sk u ). Red-dashed line: skewness of spanwise velocity component (Sk v ). Green-dashed line: skewness of vertical velocity component (Sk w ).

Figure 3 . 9

 39 Figure 3.9 Vertical profiles of mean streamwise velocity plotted on a log scale, where ⟨u⟩ is the horizontally and temporally averaged streamwise velocity above the cubes.

Figure 3 .

 3 Figure 3.10 Turbulent spectra of the streamwise velocity component E(kz ′ ) at P2 from LES (Solid blue line) compared with wind-tunnel data (circles) from Castro et al. (2006). The black dash-line represents the -5/3 slope. The red dash-line and magenta dot-dash-line indicate the minimum and maximum wavenumber resolved in the LES.

Figure 3 .

 3 Figure 3.11 Mean flow structure in a vertical x-z plane through the middle of the cube in the staggered array showing (u, w) wind vectors.

Figure 3 .

 3 Figure 3.12 Top view of the horizontal transect (x,y) near the ground (z = h/40) with planar streamlines. Wind goes from left to right.

Figure 3 .

 3 Figure3.13 shows the streamlines pattern in a 3D view. Streamlines are generated upstream and downstream from three vertical lines starting at the ground and going up to z = h/2 at location P1 (Fig.3.13 a), P2 (Fig.3.13 b), and P3 (Fig.3.13 c). Streamlines generated at P1 show that the flow impacting the top of the wind faces of cubes "A" and "B" (Fig.3.13 a) is sinking along the cubes faces and converges in the wake of cube "C", where is goes backwards. The flow is then lifted up along the lee-side of cube "C" and goes finally back in the streamwise direction. Streamlines generated at P2 (Fig.3.13 b) show that the incoming flow on the lower part of a cube is "trapped" in the recirculation at the foot of the cube and slowly goes to the side of the cube. Streamlines generated at P3 (Fig.3.13 c) show that, below z = h/2, the flow between cubes comes from the two recirculation area of the two upstream cubes and then goes straight in the streamwise direction after P3.

Figure 3 .

 3 Figure 3.13 3D view of the streamlines emanating from vertical lines at P1 (a), P2 (b) and P3 (c). Streamlines color denotes the mean velocity. Green lines indicated the streamlines sources.

2 )Figure 3 .

 23 Figure 3.14 Visualization of the flow over the computational domain: iso-surface of the magnitude vorticity (| -→ ω | = 200 [s -1 ]).

Figure 3 .

 3 Figure 3.15 Flow visualization over the computational domain: instantaneous iso-surface of the Q-criterion (Q = 5000 [s -2 ]).

Fig. 3 .

 3 Fig. 3.15 shows an iso-contour of the mean flow Q-criterion (Q = 5000 [s -2 ]). This visualization plot provides a qualitative understanding of the flow complexity and dynamic features. The shear layer emanating from the sharp edges of the cube immediately separates and generates distinct boundary layers within the roughness sublayer. The high Reynolds number flow leads to strong interactions between different scales, and the individual shape and organization of the vortex are not clear to present. The vortex structure in the wake region may be induced by the recirculation on the edging of the leeward surface (Fig. 3.12).

  , Xie and Castro (2006), and Xie et al. (2008). However, few studies have been performed to investigate the turbulence kinetic energy (TKE) budget over the urban canopy (Tomas et al. 2016), although it has been well studied over plant canopy (Moeng 1984; Dwyer et al. 1997; Finnigan 2000; Yue et al. 2008). The TKE budget is critical for understanding the turbulent exchange in the atmospheric boundary layer, especially in the roughness sub-layer. Based on the conclusions of the wind-tunnel experiments of Brunet et al. (1994) and Castro et al. (2006), the production is approximately balanced by the dissipation in the surface layer, and the turbulent transport which is a major energy sink. Blackman et al. (2017) investigated in detail the TKE budget over staggered cube arrays in boundary layer wind tunnel using an LE-PIV model to estimate the dissipation rate. Giometto et al. (2016) performed LES over a realistic urban model.

Figure 4 .

 4 Figure 4.1 depicts all contributions to the TKE budget at P0-P3. Note that D ν is very small for high Reynolds number flows and not shown here for clarity.

  .

Figure 4 . 1

 41 Figure 4.1 Vertical profiles of TKE budget terms (a) behind cube (P1), (b) in front of cube (P2), (c) in the gap between two cubes (P3) and (d) on cube (P0), all terms are normalized by u 3 τ /h.

Figure 4 . 2

 42 Figure 4.2 Vertical profiles of |P|/|ε sgs + ε r | ratio (a) behind cube P1, (b) in front of cube P2, (c) in gap P3 and (d) on cube P0.

  TKE budget from LES at location P2 (Fig. 4.3 (a)) and P3 (Fig. 4.3 (b)) are compared to the experimental data from Blackman et al. (2017) and Castro et al. (2006). Only production and dissipation are shown here for clarity. Note that in the literature, only few experimental studies are available for comparison with the present work.

Figure 4 . 3

 43 Figure 4.3 Production and dissipation terms of TKE budget compared to the wind-tunnel data from Castro et al. (2006) and Blackman et al. (2017) at location P2 (a) and P3 (b). All terms are normalized by u 3 τ /h.

Figure 4 .

 4 Figure 4.4 (a) Decomposed production terms from LES at P3, all terms are normalized by u 3 τ /h. (b) Resolved strain rate tensor (S i j ) at P3, normalized by u τ /h.

4. 4

 4 (a)). Additionally, negative production P is observed in the lower part of the canopy at P3 (Fig. 4.3 (a)). It is the result of the span-wise velocity component fluctuation and the span-wise velocity gradient term (-v ′ v ′ S 22 ). This phenomenon is consistent with the observation of Blackman et al. (2017), who found that the span-wise term (-v ′ v ′ S 22 ) acts as a sink of energy instead of a source in the lower part of the canopy (below z/h = 0.5). Moreover, the contribution of S 12 and S 23 are negligible all over the vertical profile (Fig. 4.4 (b)) resulting in negligible contribution of (-2 u ′ v ′ S 12 ) and (-2 u ′ v ′ S 23 ) in the central plane of the cubes.

Figure 4 . 5

 45 Figure 4.5 Three-dimensional visualization of Production, normalized by u 3 τ /h. For clarity, only three cubes are shown. (a) Zoom inside the urban canopy with iso-contours of P = 10 (orange) and P = -10 (light blue). Iso-contours are shown only below z = 3h/4. (b) View of iso-contours of Production P = 25 (red) and P = -25 (dark blue).

Figure 4 . 6

 46 Figure 4.6 Top view of the horizontal transect (x,y) of Production normalized by u 3 τ /h. Wind goes from left to right. (a) Near the ground (z = h/40) with planar velocity streamlines where S, N and F locations stands for saddle point, node and focus node, respectively. (b) In the middle of the cube (z = h/2) with mean velocity vector field. Vectors are of equal length to better visualize flow structures.

Figure 4

 4 Figure 4.7 (a) Vertical transect of Production normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures. (b) Production in the windward face of a cube with wall streamlines.

Figure 4 . 8

 48 Figure 4.8 Vertical transect of Turbulent transport normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures. The pink contour circles is Production (P = 25).

Figure 4 . 9

 49 Figure 4.9 Three-dimensional visualization of Turbulent transport (T r ), normalized by u 3 τ /h. For clarity, only three cubes are shown. (a) View of T r = 5 (pink). (b) View of iso-contours of Turbulent transport T r = -8 (dark blue).

4

 4

  .10). The high positive T r moves downward along the vertical edges of the windward side to the foot of the cube (Fig.4.9 (a) and Fig. 4.8).

Figure 4 .

 4 Figure 4.10 Top view of the horizontal transect (x,y) of Turbulent transport (T r ), normalized by u 3 τ /h, in the middle of the cube (z = h/2) with mean velocity vector field. Vectors are of equal length to better visualize flow structures. Wind goes from left to right.

Figure 4 .

 4 Figure 4.11 Three-dimensional visualization of Pressure transport, normalized by u 3 τ /h. For clarity, only three cubes are shown. (a) View of iso-contours of T p = 15 (pink). (b) View of iso-contours of Pressure transport T p = -5 (dark blue).

  .

  4.11 (b)).

Figure 4 .

 4 Figure 4.12 Vertical transect of Pressure transport (T p ) normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures.

  4.13 (b) and Fig. 4.14).

Figure 4 .

 4 Figure 4.13 Three-dimensional visualization of Advection, normalized by u 3 τ /h. For clarity, only three cubes are shown. (a) View of iso-contours of Advection A = 5 (pink). (b) View of iso-contours of Production P = 25 (red) and of Advection A = -5 (dark blue).

  4.13 (b) and and Fig. 4.14). The negative A here can be explained by the positive streamwise gradient of TKE and positive streamwise velocity (see Fig. 4.15). Additionally, the negative A close to S1 (Fig. 4.14) can be explained by negative streamwise gradient of TKE and reversed streamwise velocity. It is also visible that the positive advection A follows immediately the negative A region in streamwise direction (see Fig. 4.16).

Figure 4 .

 4 Figure 4.14 Vertical transect of Advection (A) normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures. The pink contour circles is Production (P = 25).

Figure 4 .

 4 Figure 4.15 Vertical transect of TKE normalized by u 2 τ in the middle of the cube with a projection of the mean velocity vector field. Vectors are following the value of time-averaged velocity.

Figure 4 .

 4 Figure 4.16 Top view of the horizontal transect (x,y) of Advection normalized by u 3 τ /h, in the middle of the cube (z = h/2) with mean velocity vector field. Vectors are of equal length to better visualize flow structures. Wind goes from left to right.

Figure 4 .

 4 Figure 4.17 Three-dimensional visualization of iso-contours of Dissipation ε r + ε sgs < -5, normalized by u 3 τ /h. For clarity, only three cubes are shown.

Figure 4 .

 4 Figure 4.18 Vertical transect of dissipation (ε r + ε sgs ) normalized by u 3 τ /h in the middle of the cube with a projection of the mean velocity vector field. Vectors are of equal length to better visualize flow structures.

Figure 4 .

 4 Figure 4.19 (a) -(f) The decomposition of Turbulence transport (T r ) at location P3 into the contribution of each velocity gradient compared to the wind-tunnel data from Blackman et al. (2017). All terms normalized by u 3 τ /h.

  . 4.19 (a) -(f) in comparison with the wind-tunnel data from Blackman et al. (2017). A good agreement is observed both above and within the canopy. The term -source at the canopy height, while it acts as a sink in the canopy (Figs. 4.20 (a)). The main contribution is from is larger than in the wind-tunnel experiment, and the location of the peak is higher in LES (Fig. 4.19 (b)). Overestimated value are shown in the upper part of the canopy in Figs. 4.19 (c),(d),(f). This phenomenon may be due to the fact that the large size in wind-tunnel experiments weaken the intensity of velocity fluctuations and turbulent transport. However, very good agreement is observed in Fig. 4.19 (e).

P3

Figure 4 .

 4 Figure 4.20 (a) Turbulent transport (T r ) compared to the wind-tunnel data from Blackman et al. (2017) and T r without three span-wise gradients. (b) Three span-wise gradients of Turbulent transport (T r ). All terms normalized by u 3 τ /h.

4. 6 . 2

 62 The impact of ignoring non-measured terms

Figure 4 .∂

 4 Figure 4.20 (a) shows, at P3, the comparisons between T r completely computed from the LES, T r from LES data but without span-wise gradients and T r obtained from wind-tunnel experiment (Blackman et al. 2017). The agreement between the present T r and that from the wind-tunnel data is generally good above the canopy and close to the ground. A large discrepancy is observed in the middle of the canopy 0.4 < z/h < 1. It reflects that ignoring the three span-wise gradients results in an underestimated T r inside the canopy. Indeed, when the spanwise gradient are ignored in the present CFD results, the T r agrees much better with the experimental data. The three uncalculated spanwise velocity gradients in the wind-tunnel experiment are shown in Fig.4.20 (b). Their contributions are principally concentrated within the canopy. The two terms -1 2

Figure 4 .

 4 Figure 4.21 3D visualization of the sum of the spanwise gradients of the turbulent transport terms discrepancy between the Turbulence transport and Turbulence transport without spanwise velocity fluctuation gradients. Red color means ∆T r × h/u 3 τ = 5 and blue color means ∆T r × h/u 3 τ = -5.

Figure 4 .

 4 Figure 4.21 presents ∆T r , the contribution of the spanwise terms of turbulence transport (T r ) defined as

Figure 5 . 1

 51 Figure 5.1 Vertical profile of drag coefficient C D using the equation from Maché (2012) and extracted from obstacle resolving method, compared to the data from Cheng and Castro (2002) and Coceal et al. (2006).

Figure 5 . 2

 52 Figure 5.2 Vertical view of the mesh generated by the OpenFOAM tool: blockMesh.

Figure 5 .

 5 Figure 5.3 (b) shows the ⟨u⟩ profile normalized by the mean velocity at the top of the cubes ⟨u⟩ h . Both drag-porosity approaches overestimate the ⟨u⟩ profile of OR-LES ESA in the range of 0.8 < z/h < 1. From the ground up to z = 0.8h, DP-LES M fits the OR-LES ISA profile better than the DP-LES R , but DP-LES R is closer to the OR-LES ESA profile than DP-LES M . However, the shape of the mean velocity profile (and in consequence the vertical velocity gradients) inside the canopy obtained with C D -Mache shows the best agreement with both OR-LES data.

Figure 5 . 3

 53 Figure 5.3 Vertical profiles of (a) spatial and temporal streamwise mean velocity ⟨u⟩, normalized by ⟨u⟩ 4h : comparison between drag-porosity LES and obstacle-resolved LES using intrinsic and extrinsic method. (b) spatial and temporal streamwise mean velocity ⟨u⟩ within the canopy, normalized by ⟨u⟩ h . Porosity C D -resolved (red solid line), porosity C D -Maché (blue dashed line), obstacle-resolved LES using intrinsic spatial averaging (white circle) and extrinsic spatial averaging (black dot).

Figure 5 . 4

 54 Figure 5.4 Vertical profiles of the spatial averaged (a) standard deviation of the streamwise velocity, ⟨σ u ⟩; (b) standard deviation of the spanwise velocity ⟨σ v ⟩; (c) standard deviation of the vertical velocity ⟨σ w ⟩; (d) Reynolds shear stress, normalized by u * ; (e) turbulent kinetic energy (TKE), normalized by u 2 * .

  Figure 5.4 (b)) and ⟨σ w ⟩ (Figure 5.4 (c)), above the height z = h, the DP-LES M gives a better fit to the obstacle resolved data than DP-LES R . Similarly to ⟨σ u ⟩, lower values of ⟨σ v ⟩ (Figure 5.4 (b)) and ⟨σ w ⟩ (Figure 5.4 (c)) are observed in DP-LES inside the canopy. This may be due to the lack of a rigid shape of the obstacle resulting in a reduction in velocity fluctuations inside the canopy. The TKE from DP-LES M shows a good agreement in comparison with OR-LES data above z > 2h, and an underestimated value below (Figure 5.4 (e)).

Figure 5 . 5

 55 Figure 5.5 Vertical profiles of spatial and temporal averaged (a) skewness of streamwise velocity component (⟨Sk u ⟩); (b) skewness of spanwise velocity component (⟨Sk v ⟩); (c) skewness of vertical velocity component (⟨Sk w ⟩).

Figure 5 .

 5 Figure 5.6 shows the spectrum from the DP-LES M added into Fig.3.10. The simulated energy spectrum from the from DP-LES matches very well the wind-tunnel data and OR-LES data in the inertial subrange (-5/3 slope) is accurately captured (Fig.3.10). After k max z ′ , the energy of the small scales starts to drop more rapidly in the DP-LES simulation comparing to the wind-tunnel data and OR-LES data. For low normalized wavenumbers, there is much lower energy in the DP-LES simulation for kz ′ < k min z ′ , but in general, the deficit of the energy spectrum is not too significant.

Figure 5 . 6

 56 Figure 5.6 Turbulent spectra of the streamwise velocity component E(kz ′ ) at P2 and z/h = 1.62 from OR-LES (Solid blue line) and from DP-LES (Solid cyan line) compared with wind-tunnel data (circles) from Castro et al. (2006). The black dash-line represents the -5/3 slope. The red dash-line and magenta dot-dash-line indicate the minimum and maximum wavenumber resolved in the LES.

Figure 5 .Figure 5 . 7

 557 Figure 5.7 Vertical profiles of spatial averaged TKE budget terms normalized by u 3 τ /h, (a) OR-LES ISA and (b) DP-LES M .

Figure 5 .

 5 Figure 5.7 (a) shows that the production and dissipation balance above z = 3h. A rather sharp peak of production is observed at height z = h. The Advection term, of negative sign at z = h, shows a similar sign as at P1 (Fig.4.15). The pressure transport is, with the space-average, nearly zero at z = h.

Figure 5 .

 5 Figure 5.8 shows the relative number of each event and the relative contribution to the Reynolds stress u ′ w ′ for the DP-LES compared to the DNS data from Coceal et al. (2007b)

(

  [START_REF] Finnigan | Turbulence in plant canopies[END_REF][START_REF] Coceal | Structure of turbulent flow over regular arrays of cubical roughness[END_REF]).

Figure 5 .

 5 Figure 5.8 (a) Relative number of events in each quadrant as a function of height. (b) Relative contribution to u ′ w ′ of events in each quadrant as a function of height. Blue: Q1; Red: Q2; Green: Q3; Yellow: Q4. Points: DNS from Coceal et al. (2007b); dashed line: DP-LES data.

Figure 5 . 9

 59 Figure 5.9 Instantaneous snapshot showing contour plots of streamwise velocity at two altitudes: z = 1.5h (left) and z = 2h (right) from Coceal et al. (2007b). Black regions: u < 0.8u. White regions: u > 1.2u. Here u is the local mean velocity.

Figure 5 .

 5 Figure 5.10 Instantaneous snapshot showing contour plots of streamwise velocity at different altitudes: a) z = 1.5h, b) z = 2h, c) z = 3h and d) z = 4h of the OR-LES result. Blue regions: u < 0.8u. Red regions: u > 1.2u. Here u is the local mean velocity.

Figure 5 .

 5 Figure 5.10 shows the instantaneous visualization of LMRs on the horizontal plane at four different altitudes above the canopy from the OR-LES simulation. The red regions correspond to areas where u ′ is greater than 20% of the local average and the blue regions correspond to areas where u ′ is less than 20% of the average value. The low velocity streaks (the blue regions) are elongated in the streamwise direction in agreement with the DNS data. There are numerous LMRs observed near the canopy height, and the number of LMRs streaks is generally reducing with the higher altitude. The length of LMRs streaks in the spanwise direction L y increases with the higher altitude from L y = 2h (Figure 5.10 (a)) to L y = 4h (Figure 5.10 (c)). The high velocity streaks occupy more space and are much more numerous in the OR-LES simulation than the DNS results (Figure 5.9).

Figure 5 .

 5 Figure 5.11 Instantaneous snapshot showing contour plots of streamwise velocity at different altitudes: a) z = 1.5h, b) z = 2h, c) z = 3h and d) z = 4h of the DP-LES result. Blue regions: u < 0.8u. Red regions: u > 1.2u. Here u is the local mean velocity.

Figure 5 .

 5 Figure 5.11 shows the instantaneous visualization of LMRs at four different altitudes above the canopy from DP-LES.The LMRs streaks elongated along the streamwise direction are observed in the DP-LES simulation, but they are much thicker LMRs compared to the OR-LES and DNS data at the near roof level with L y = 3h (Figure 5.11 (a)). The LMRs are less thick in OR-LES simulation due to the staggered arrangement of the cubes changing the direction of the flow and reducing the magnitude of streamwise velocity, resulting in the separation of the LMRs streaks. In addition, the high velocity streaks are less presented in the DP-LES in agree with the data of the DNS (Coceal et al. 2007b) (see Figure 5.11 (a)). The high velocity streaks in the DP-LES starts to disappear from the height z = 3h (Figure 5.11 (c)), while in the OR-LES simulation is z > 4h (Figure 5.10 (d)).

Figure 5 .

 5 Figure 5.12 Two-point correlation plot in the XY plot from the wind-tunnel experiment from Rivet (2014). (a) and (c) from the DP-LES simulation; (b) and (d) from Rivet (2014). (a) and (b) at z = 1.5h; (c) and (d) at z = 3h.

5. 5

 5 Suggestion for improving the drag-porosity approach 5.5.1 Separately model the features in the near-surface region

Figure A. 1

 1 Figure A.1 Flow chart of the PIMPLE algorithm

  Effect of the urban canopy roughness . . . . . . . . . . . . . . . . . . . . . Representation of the urban canopy in numerical simulations . . . . . . . . 1.3.1 Obstacle resolving method . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Immersed boundary method . . . . . . . . . . . . . . . . . . . . . 1.3.3 Drag-porosity approach . . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Aerodynamic roughness length . . . . . . . . . . . . . . . . . . . . 1.4 Statistical variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 First-order statistics . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Second-order statistics . . . . . . . . . . . . . . . . . . . . . . . . 1.4.3 Friction velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . Pressure transport [m 2 s -3 ] T sgs Subgrid transport [m 2 s -3 ] ε r Viscous dissipation [m 2 s -3 ] ε sgs Subgrid dissipation [m 2 s -3 ]

	τ		Shear stress [kgm -1 s -2 ]
	ν		Kinematic viscosity[m 2 s -1 ]
	ν sgs Subgrid-scale viscosity [m 2 s -1 ]
	ρ		Reference density [kgm -3 ]
	U bar Desired mean velocity [ms -1 ]
	A		Advection [m 2 s -3 ]
	P		Production [m 2 s -3 ]
	T r		Turbulent transport [m 2 s -3 ]
	T p	Notations
		( )	Time-averaging operator
		⟨ ⟩	Horizontally spatial-averaging operator
		| |	Magnitude operator
	x y Abbreviations Streamwise coordinate [m] Spanwise coordinate [m] z ABL Atmosperic Boundary Layer Vertical coordinate [m] ISL Inertial Sub-Layer
	RSL Roughenss Sub-Layer
	UCL Urban Canopy Layer
	IBM Immersed Boundary Method
	CFD Computational Fluid Dynamics
	DNS Direct Numerical Simulation
	LES Large-Eddy Simulation
	RANS Reynolds-averaged Navier-Stokes equations
	OpenFOAM Open source Field Operation And Manipulation
	z 0 SGS Subgrid-Scale Aerodynamic roughness length [m] d TKE Turbulent Kinetic Energy Zero-plane displacement height [m] δ HWA Hot-Wire Anemometry Boundary layer height [m] κ PIV Particle Image Velocimetry Von Karman's constant t LDV Laser Doppler Velocimetry Time [s] h ISA Intrinsic Spatial Averaging Height of the cube [m] λ p ESA Extrinsic Spatial Averaging Plan area density λ f ECN Ecole Centrale de Nantes Frontal area density α f LHEEA Laboratoire de Recherche en Hydrodynamique, Energétique et Environnement Volume frontal density [m -1 ] F Drag force per unit volume and density [ms -2 ] Atmosphérique
		F D	Drag force [kgms -2 ]
		C D	Drag force coefficient

Introduction 1 Atmospheric boundary layer: structure and characteristics 1.1 The structure of the atmospheric boundary layer . . . . . . . . . . . . . . . 1.1.1 Outer layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Inertial sublayer . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Roughness sublayer . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.4 Urban canopy layer . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 1.2.1 Morphological parameters . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Flow classification . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 The layout of the buildings . . . . . . . . . . . . . . . . . . . . . . 1.3 1.4.4 Skewness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.5 Two-point correlation coefficient . . . . . . . . . . . . . . . . . . . u Instantaneous streamwise velocity component [ms -1 ] v Instantaneous spanwise velocity component [ms -1 ] w Instantaneous vertical velocity component [ms -1 ] u ′ Fluctuation of the streamwise velocity component [ms -1 ] v ′ Fluctuation of the spanwise velocity component [ms -1 ] w ′ Fluctuation of the vertical velocity component [ms -1 ] u τ Friction velocity [ms -1 ] u * Friction velocity [ms -1 ] u h Velocity at the top of the canopy [ms -1 ]

Table 1 .

 1 .1. 1 Terrain classification of aerodynamic roughness length z 0 fromWMO (2008)

	Terrain description	z 0 (m)
	Open sea, Fetch at least 5 km	0.0002
	Mud flats, snow, no vegetation, no obstacles	0.005
	Open flat terrain; grass, few isolated obstacle	0.03
	Low crops; occasional large obstacles, x/H > 20 0.10
	High crops; scattered obstacles, 15 < x/H < 20	0.25
	Land, bushes; numerous obstacles, x/H ≈ 10	0.5
	Regular large obstacle coverage (suburb, forest)	1.0
	City centre with high-and low-rise buildings	> 2

  al. 2010; Wong and Liu 2013), and to study the inflow conditions of the simulation (Xie and Castro 2008; Foti et al. 2017). Recently, Yang et al. (2016) performed LES within the roughness sublayer of boundary layer flow over rough surfaces.

Table 2 .

 2 1 Interpolation schemes from OpenFOAM user guide (2018)

	Interpolation schemes Numerical behaviour
	linear	Second order, unbounded
	skewLinear	Second order, (more) unbounded, skewness correction
	cubicCorrected	Fourth order, unbounded
	upwind	First order, bounded
	linearUpwind	First/second order, bounded
	QUICK	First/second order, bounded

Table 2 .

 2 2 Time schemes from OpenFOAM user guide (2018)

	Time schemes Numerical behaviour
	steadyState	sets time derivatives to zero.
	Euler	transient, first order implicit, bounded.
	backward	transient, second order implicit, potentially unbounded.
	CrankNicolson transient, second order implicit, bounded.
	localEuler	pseudo transient for accelerating a solution to steady-state
		using local-time stepping; first order implicit.

Table 2 .

 2 Smagorinsky model, TKE SGS model and dynamic Smagorinsky model. Obstacle-resolved method is employed to simulate the staggered cubic arrays. The mesh resolution is ∆ = h/32 for all cases (see section 2.4.3). The simulation is performed with a time step of 0.002T over a total duration 400T , including an initial duration of 200T , and a further duration of 200T used for the high-order statistics, where T = h/u τ can be interpreted as the eddy turnover time for the largest eddies shed by the cube[START_REF] Coceal | Mean flow and turbulence statistics over groups of urban-like cubical obstacles[END_REF]). The four locations (P0-P3) in Fig.2.7 will be used to compare the mean flow and high-order statistics with the literature data. More details of each simulation are listed in Table2.3.

	TKE SGS	Staggered 0.25 4h × 4h × 4h 0.02	0.357	452 ∆ = h/32
	Smagorinsky Staggered 0.25 4h × 4h × 4h 0.02	0.384	486 ∆ = h/32
	DS*	Staggered 0.25 4h × 4h × 4h 0.02	0.36	456 ∆ = h/32

3 Domain parameters information for SGS model comparison SGS type Array type λ p L x × L y × L z h [m] u τ [m/s] Re τ mesh *Dynamic Smagoinsky is noted as DS Time-averaged streamwise velocity

Table 3 .

 3

			1 Summary of parameters for urban-like arrays simulation
	Array type λ p	L x × L y × L z	h	u τ	Re τ	Re	U bar
	Staggered	0.25 16h × 12h × 8h 0.02 [m] 0.38 [m/s] 481 50000 3.5 [m/s]

Table 5 .

 5 1 Summary of parameters in the comparison

	Simulation C D	Array type λ p SGS C S	u *	Cell number
	OR-LES	--	Staggered 0.25 DS * --	0.379 28 Millions
	DP-LES M C D -Mache	Staggered 0.25 DS * --	0.456 2.7 Millions
	DP-LES R C D -Resolved Staggered 0.25 S * *	0.167 0.323 1.6 Millions

* Dynamic Smagoinsky is noted as DS ** Smagoinsky is noted as S, this simulation is performed before the Dynamic Smagoinsky model achieved.

https://www.un.org/development/desa/en/news/population/2018-revision-of-world-urbanizationprospects.html

Acknowledgements

First of all, I want to thank Dr. Laurent Perret. In the final year of the Master program, I had the opportunity to complete a six-month internship under his guidance. This work opened my vision of fluid mechanics research and inspired my interest in airflow modeling in urban environments.

I am also very gratefully acknowledge the financial support of the PhD scholarship from China Scholarship Council (CSC) under the grant CSC N • 20158070084, and the financial support of the French National Research Agency through the research grant URBANTURB N • ANR-14-CE22-0012-01. very grateful to my parents and my fiancee Dong Ru for understanding and encouraging in the process of completing this thesis. The completion of this paper is inseparable from their support. I love you always forever.

At z = h/2, only two main production regions are remaining. They are attached to the two vertical upwind edges of the cube and get stronger with the increasing height and, eventually, merge into the spanwise space between two cubes (Fig. 4.6 (b) and Fig. 4.5 (a)). These regions are corresponding to the side of a recirculation area identified by mean velocity vectors (Fig. 4.6 (b)). A secondary P peak is visible on the lower part of the windward face of the cube near point S4 (Fig. 4.7 and Fig. 4.5 (a)). It is an elongated span-wise structure that may be linked to a small counter-clock-wise secondary recirculation at the foot of the cube observed in Fig. 4.7.

Going up from z = h/2 (see Fig. 4.5 (b)), high production regions are observed on the two vertical edges of the windward side of the cube and the four edges of the top face. The strongest production area is seen starting from the trailing edge of a cube (downwind edge of the top face), and expanding downwind to reach the second row of cubes with high production (P > 25). The intensity of production in this area decreases with downwind distance from the cube and spreads in all direction but mainly upwards (Fig. 4.7). All high production regions around the cube are linked to the sharp corners with flow separation and to the high shear regions. Indeed, high P regions described are connected to the separation regions. Close to a surface, they are represented by saddle points S1, S2 and S4 where streamlines are converging/diverging with a great updraft. On cubes, they are situated on the edges, where flow separates. These regions of flow convergence may be linked to sources of vortex shedding.

The flow also presents areas where P < 0. The negative production reflects that energy is extracted from turbulence and returned to the mean flow [START_REF] Pope | Turbulent Flows[END_REF].

For z < h/2, a weak span-wise elongated negative production zone (P ≈ -15) is observed at the ground near point N1 (Fig. 4.6 (a) and Fig. 4.7 (a)), right upwind of a cube, parallel to the first high production region previously described (S1). Point N1 is a repelling node (here it is more a repelling line, see Helman and Hesselink 1989) representing flow impinging on a surface. Here it is very close to the re-attachment zone situated between the big clock-wise recirculation upwind of the cube and the secondary counter-clock-wise recirculation at the cube's foot.

The region with the most negative production (P < -25) starts at the cube's leading edge and extends downwards along the windward side of the cube with decreasing intensity to reach P = -10 at approximately z = h/2 (Fig. 4.7 (a) and Fig. 4.5 (b)). This region, denoted Chapter 6

Conclusions and perspectives

The high fidelity of numerical simulations will help to obtain better predictions results for the analysis of unsteady turbulence in urban environments. At the same time, the study of the dynamic characteristics of turbulence in the city area also helps to understand the formation and transport mechanism of turbulent kinetic energy, and it will be of great help to the research of pollutant dispersion and coherent structures.

Main results and achievements

Implementation of a dynamic Smagorinsky model in OpenFOAM

In LES, the very common Smagorinsky eddy-viscosity model used a prescribed value to model the subgrid-scale dissipation. This parameterization process does not allow obtaining model coefficients that can be locally changed like in a dynamic model. However, the existing dynamic Smagorinsky model in OpenFOAM used a horizontal average to compute the parameter. In this work, the dynamic Smagorinsky model implemented maintains the local features.

The method is not limited to the case of having a homogeneous flow and is easily applied to complex geometries and unsteady flows. In order to be able to implement the model with minimal computational instabilities, we considered Calmet and Magnaudet. (1997)'s suggestion that the local total viscosity is set to a non-negative value to keep the simulation stable.

Appendix A

PimpleFOAM Solver

To sum up, the PimpleFOAM algorithm has three layers of loops per time step: • The middle layer is PIMPLE corrector loop (pimple.correct()) similar to PISO loop , and is responsible for the pressure correction.

• Inner layer is non-orthogonal pressure corrector loop

• Then, the new pressure is used as the initial value to solve pressure equation

Dynamic Smagorinsky SGS model

In this work, the dynamic Smagorinsky SGS model is implemented in OpenFOAM © 2.4.0. As suggested by [START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF], the dynamic coefficient C is computed at each time step and each computational cell using least-squares method, the code is: However, the dynamic procedure may locally induce large negative values of the model coefficient C (see Figure 2.2), giving rise to numerical instabilities. In order to avoid this problem, the total viscosity (ν + ν sgs ) is forced to be locally a non-negative value, the code in OpenFOAM © 2.4.0 as follows, where the kinematic viscosity ν = 1.58 * 0.00001.