Nomenclature

Af,	frontal area ratio;	ke,	turbulent kinetic energy in the testing case;
Ap,	projection area ratio;	Lx, Ly, Lz,	domain sizes in the streamwise, lateral, and vertical directions;
Dk0,	diffusion term in the turbulent kinetic energy budget equation in the baseline case;	l0 (z),	characteristic length scale in the baseline case;
Dke,	diffusion term in the turbulent kinetic energy budget equation in the testing case;	le (z),	characteristic length scale in the testing case;
E,	(ke/k0)1/2;	Pk0,	production term in the turbulent kinetic energy budget equation in the baseline case;
Euu,	wave spectrum of the streamwise velocity;	Pke,	production term in the turbulent kinetic energy budget equation in the testing case;
h,	side length of the cubical roughness;	δ,	boundary layer depth;
k0,	turbulent kinetic energy in the baseline case;	ρ,	air density;
,	filtered pressure;	x, y, z,	coordinate in the streamwise, lateral, and vertical directions;
Rij0 (z),	Reynolds stress in the baseline case;	νt,	subgrid-scale eddy-viscosity;
Rije(z),	Reynolds stress in the testing case;	Φij,	redistribution term in the budget equation of ;
,	deformation tensor of the resolved velocity;	τij,	subgrid-scale stress;
P,	averaged mean pressure;	εk0,	dissipation term in the turbulent kinetic energy budget equation in the baseline case;
ν,	kinematic viscosity coefficient;	εke,	dissipation term in the turbulent kinetic energy budget equation in the testing case.
U, V, W,	averaged mean velocities in the streamwise, lateral, and vertical directions, respectively;
ue,	characteristic turbulent velocity in the testing case;
u0,	characteristic turbulent velocity in the baseline case;
,	filtered velocities in the streamwise, lateral, and vertical directions;

1 Introduction

With the ever-increasing computational resources, the eddy-resolving simulation technics, e.g., direct numerical simulation (DNS) and large-eddy simulation (LES), have been widely used in the turbulence predictions of microscale urban atmospheric environment (MUAE) flows^[1-3]. For such simulations, the underlying surface may change rapidly, and it is more appropriate to use the inflow-outflow boundary condition to simulate a spatially developing turbulent field rather than to use the periodic condition in the horizontal direction. Therefore, it is of great importance to add proper inflow turbulence, which is crucial for the contaminant dispersion^[4].

So far, many inflow turbulence generation methods have been proposed. There are basically three kinds of turbulence generation methods, i.e., synthesis methods such as the digital filtering method^[5], the inverse Fourier method^[6], the synthetic eddy method^[7], and the proper orthogonal decomposition (POD) method^[8-9], forcing methods^[10], and precursor simulation methods^[11].

Synthesis methods generate the turbulent field with some of the prescribed turbulent statistics, e.g., temporal correlation, spatial correlation, energy spectra, integral length scale, integral time scale, and Reynolds stress tensor. Alternatively, random body force can be added to accelerate the turbulence-development in a transition region^[10]. Similarly, the amplitude of the random force is based on the Reynolds stress. In the simulations of simple geometry distribution, one can run a precursor simulation, and the inlet data can then be provided by transforming the turbulent field of the fully developed precursor simulation^[11]. Schlüter et al.^[12] transformed a databank to generate turbulent field for the coupling boundary of the hybrid approach of RANS and LES, in which the precursor time series are transformed based on the ratio between the normal stress of the target problem and that of the databank. To sum up, second order statistics of turbulent field are needed for most of inflow turbulence generation methods.

It is important to know how each of the statistics affects the turbulence-development. Keating et al.^[13] tested the spectrum effect of the inflow turbulence in a plane channel flow, and found it important to capture the eddies with dimensions equal to or larger than the integral length scale of the flow. However, how the turbulent statistics affect the short-distance development of turbulence in MUAE flows, e.g., within 2δ, especially in the roughness sublayer or the canopy layer, which is of great importance to human activities, remains unsolved. Schlüter et al.^[12] found that, for confined strong swirling jet, in which the high shear regions were prominent, the test case with only quasi-laminar inlet flow yielded a flow field that agreed reasonably well with the experimental data. The results suggested that, for different types of flows or different parts of a single flow field, the effect of the inlet turbulence on the development could be quite different.

Besides, it is not easy to get such turbulent statistics for a genuine MUAE flow. Particularly, in the roughness sublayer, the turbulent statistics change substantially for different packing densities^[14], and the turbulent field around the buildings with the same height is quite different if the layout of the surrounding buildings is different^[15]. In the simulations of MUAE flows with the multi-scale coupling method, in which the boundary velocity fields are provided by meso-scale meteorological models^[1], the turbulent kinetic energy derived from the mesoscale models may lead to large errors^[16]. Moreover, the turbulent kinetic energy in the canopy layer cannot be provided by meso-scale models because of the coarse grid in the near-surface region and the employed urban canopy models.

The target of the present work is to investigate the effect of inflow turbulence statistics on the development of the turbulent field in MUAE flows. Among all the turbulent statistics aforementioned, the Reynolds stress and spectra of the inflow turbulent are focused on. How different parts, e.g., the inner-layer which denotes the roughness sublayer hereafter, i.e., z < 2h in the present work, and the outer-layer which denotes the region above the roughness sublayer hereafter, i.e., z > 2h, affect the turbulence-development is also investigated.

The paper is organized as follows. In Section 2, the numerical procedure and physical model are introduced. In Section 3, a theoretical analysis and the numerical results are given. Conclusions are made in Section 4.

2 Numerical procedure and physical model 2.1 Numerical procedure

Neutral MUAE flows are considered, which are governed by incompressible Navier-Stokes equations. The LES, which has been widely used in simulating the turbulent statistics and flow structures^{[2, 15, 17-20]} in MUAE flows, is applied. The filtered continuity equations and the momentum equations are

(1)

(2)

where i, j=1, 2, 3, and is the filtered pressure. dP/dx₁ is the mean pressure gradient in the streamwise direction. It balances the total drag and the subgrid-scale stress τ_ij defined by

which is closed by the Lagrangian dynamical model as follows:

(3)

where is the subgrid-scale eddy-viscosity, and Δ is the filter scale. The model coefficient C_s is determined dynamically by averaging over the flow pathlines^[21]. is the deformation tensor of the resolved velocity, i.e.,

The air density ρ is 1.208 kg m^-3, and the kinematic viscosity coefficient υ is 1.5 ×10 ^-5m²s^-1.

The equations are discretized with the finite volume method on non-staggered grids. The filtering is conducted in the finite volume implicitly, and thus the filter scale is the grid resolution. A third-order explicit Runge-Kutta scheme is applied in time, and the time step Δt is 2 × 10 ^-3T, where T = h/u_∗ is the turnover time for the largest eddy shed by the buildings, and u_∗ is the friction velocity. Such a time step ensures the numerical stability and that the Courant-Friedrichs-Lewy (CFL) number is smaller than 0.4. A quadratic upstream interpolation for convective kinematics (QUICK) scheme is applied in space, and the momentum interpolation is used to avoid non-physical pressure oscillation. The initial duration of the simulation is no less than 100T when the turbulent kinetic energy of the motoring point is statistically steady. Another 600T is simulated to obtain sufficient samples for time averaging.

Henceforth, x, y, and z are the same as x₁, x₂, and x₃, respectively, and u, v, and w are the same as u₁, u₂, and u₃, respectively.

2.2 Physical model and case description

The flow over a cluster of cubical roughness distributed regularly on the wall is used to simulate the MUAE (see Fig. 1). The cubes with a side length h are staggered in the streamwise direction, and are aligned in the lateral direction. The projection area ratio A_p of the project area to the total floor area and the frontal area ratio A_f of the frontal area to the total floor area are both 0.25.

The shear free boundary condition is employed on the upper bound of the domain, i.e., , and w = 0. On the solid wall including ground and building surfaces, the non-slip condition is used. A wall model is used because the atmospheric flow is a high Reynolds number flow and the first grid point can be beyond the viscous sublayer. The wall model is from Grötzbach et al.^[22] and the friction velocity u_∗ is calculated with the assumption that the instantaneous velocity of the first grid point satisfies a logarithmic or linear distribution as follows:

(4)

where is the velocity of the first grid point, z⁺ = zu_∗/ν, κ = 0.41, and B = 5.6.

The periodic condition is used in the lateral direction, because the underlying surface is supposed to extend infinitely in the lateral direction.

The inflow data is specified with the turbulent field generated by different methods, and the non-reflecting boundary condition is used at the outlet, i.e., , where ϕ is a dependent variable, and U is the averaged streamwise velocity at the outlet.

The sketches of the computational domain are shown in Figs. 1(a) and 1(b). The size of the domain L_x × L_y × L_z is 22h × 8h × 8h. The effects of the computational domain size on the flow have been tested in Ref. ^[23], where the domain sizes of 4h×4h×4h, 8h×8h×4h, 4h×4h×6h, and 12h × 12h × 8h are tested, and it is concluded that the smallest domain 4h × 4h × 4h is not able to capture the largest scale turbulence. Besides, the differences in the mean velocity and turbulence statistics are negligible except near the top boundary in that paper. Therefore, the domain height 8h is supposed to be adequate. The domain width is chosen as 8h, which is based on the following two factors:

(ⅰ) The lateral domain size should be several times larger than the lateral integral length scale, which is less than h^[24].

(ⅱ) The largest structure should be contained in the lateral direction, of which the size is supposed to be smaller or approximately the same as the domain height 8h.

Therefore, 8h is a proper choice for the domain width. The domain length should be large enough to capture the streamwise development of turbulence, and is thus set to be 22h, which is also far larger than the streamwise integral length scale.

A precursor simulation P1, which uses the periodic condition in the streamwise direction, is carried out, and the velocity time series of a yz-plane are saved. The sketches of the computational domain of the precursor simulation are shown in Figs. 1(c) and 1(d), and the size of the computational domain is 12h×8h×8h. The original precursor time series or transformed time series are used as the inlet flow for the testing cases to investigate the development of the inflow turbulence in the MUAE flows.

The computational grids are shown in Fig. 2. There are 16 grid points in all the three directions of a cube, and the girds are finer near the edge of the cubes. The first grid point from the cube surface is at a distance of 0.01h. The grids are stretched in the vertical direction above z=2h. For turbulent flows over the bluff bodies, the shape drag dominates and the vortices with the length-scale larger than h/10 contribute more significantly for the mass and momentum transfer in the roughness sublayer^[15]. Xie et al.^[15] have given appropriate results for the average flow and low order statistics with the resolutions of h/16 and h/13. Therefore, the present resolution is supposed to be fine enough.

The present numerical procedure is the same as that of Wang et al.^[20], in which the LESs of urban-like cube arrays were carried out and the results showed good agreement with those of Cheng and Castro^[25]. In Ref. [25], the roughness Reynolds number ranged from 361 to 868, and it was concluded that the dependence of the drag coefficient on the Reynolds number was weak. In the present work and Wang et al.^[20], the roughness Reynolds number is about 1 000, which indicates that the comparison is convincing and the code and numerical procedure are appropriate to investigate the low-order turbulent statistics in MUAE flows.

The parameters of the testing cases are summarized in Table 1. The mean velocity fields at the inlet of all cases are the same. In Case I0, the precursor time series is specified as the inflow turbulence. By the transformation of the precursor time series, the errors of the inlet turbulence are brought in and the development of such errors will be investigated.

The purpose of the present work are as follows:

(ⅰ) To test the effects of the turbulent kinetic energy amplitude and the characteristic length scale on the turbulence-development. The characteristic length scale is kept constant when the effects of the turbulent kinetic energy amplitude is tested, e.g., for Case IA1, all the precursor data are transformed by amplitude scaling and then imposed on the inlet of the successor domain. When the effects of the characteristic length scale is tested, the turbulent kinetic energy amplitude is kept constant by scaling the fluctuation amplitude after low-pass filtering, e.g., for Case IS1, all the precursor data are processed by the low-pass filtering followed by amplitude scaling. The low-pass filtering is conducted in the lateral direction and time. The choice of the cut-off wavenumber depends on the lateral integral length scale of the turbulence, which is approximately h in the flow over a cluster of uniform cubes^[24]. The grid resolution of the present work is Δ = h/16. Therefore, 64Δ is four times of the integral length scale, and the energy-containing structures are maintained after the filtering. The cut-off frequency is correspondingly chosen as 2π/(64dt), which ensures that the turbulence is low-pass filtered with a similar scale in the streamwise direction under Taylor's frozen turbulence hypothesis. The results of Cases IA1 and IS1 are discussed in Subsection 3.3.

(ⅱ) To test the effects of the inlet inner-and outer-layer turbulence on the turbulence-development. For Cases IA2 and IS2, the precursor data are transformed in the outer-layer to test the effects of the inlet outer-layer turbulence, while for Cases IA3 and IS3, the data are transformed in the inner-layer to test the effects of the inner-layer data. In practical processes, the precursor time series is firstly transformed by amplitude scaling or low-pass filtering, which is denoted by u_T'. Then, the combination of the original data u_o' and the transformed data u_T' are conducted by a weighted average of the two data set. For the outer-layer transformation,

(5)

For the inner-layer transformation,

(6)

The weight function H(η) is transformed from the following hyperbolic tangent function:

(7)

where η=z/δ, b is the transition height between the outer-layer and the inner-layer defined by b = 2h/δ, and α is the transition width defined by α=10. The weighting function is 0 at η=0, 0.5 at η=b, and 1 at η=1.

The results of Cases IA2, IA3, IS2, and IS3 are given in Subsection 3.4.

3 Results and discussion 3.1 Turbulent statistics of Case I0

In Case I0, the precursor time series are used as the inflow turbulence. The mean velocity and Reynolds stress of Case I0 and the precursor results are shown in Fig. 3. It is observed that the results of Case I0 are approximately the same as those of the precursor simulations. Therefore, in Case I0, no redevelopment distance is needed, and the turbulent field of Case I0 is considered to be fully developed. This will be treated as the baseline for other cases.

3.2 Mean velocity field

Figure 4 shows the vertical profiles of the streamwise velocity U and the vertical gradient of the streamwise velocity dU/dz for Cases I0, IA1, and IS1. The values of U and dU/dz of the two cases show little adjustment in the streamwise direction, and are almost identical with those for Case I0. Overall, relative to the variation in the turbulent kinetic energy, e.g., 50% at the inlet for Case IA1, the effect of the inflow turbulence on the mean steamwise velocity is small.

The lateral profiles of the streamwise and vertical mean velocity within the canopy layer at x=10h are shown in Fig. 5. The streamwise velocities of Cases IA1 and IS1 are larger than those of Case I0. The magnitude of the vertical velocity is increased for Cases IA1 and IS1, indicating that the vertical convection and meanwhile the canopy vortex are strengthened, which is consistent with the results of Kim et al.^[4]. The momentum flux downward into the canopy increases when the turbulent intensity increases, and thus the streamwise velocity in the upper level of the canopy is larger than that of Case IA0. For other locations of the inlet, the data yield the same conclusions. Therefore, only the profiles at x=10h are shown.

3.3 Turbulent statistics 3.3.1 Theoretical analysis for the development of the turbulent kinetic energy in the outer-layer

For Case I0, the turbulent flow field is homogeneous in the streamwise direction. Therefore, the mean velocity field is supposed to be

(8)

where U, V, and W are averaged temporally and in the spanwise direction. Therefore, V (z) = 0. For z > 1.5h, Eq. (8) can be further simplified as W(z) ≈ 0. The characteristic turbulent velocity, characteristic length scale, Reynolds stress, and turbulent kinetic energy are denoted by u⁰, l⁰(z), R_ij⁰ (z), and k⁰ = R_ii⁰ /2, respectively.

For other cases, the perturbation of the mean velocity is negligible (see Fig. 4), and thus Eq. (8) is still valid. The characteristic turbulent velocity, characteristic length scale, and Reynolds stress are u^e, l^e(x, z), and Re_ij^e(x, z), respectively. The relationship between u^e and u⁰ is u^e = E(x, z)u⁰. Therefore, R_ij^e (x, z) = E²R_ij⁰ (z), and the turbulent kinetic energy is

(9)

The development of E(x, z) is then governed by

(10)

which is derived in Appendix A. λ⁰ is defined by −〈u'w'〉⁰ = λ⁰k⁰. In the derivation of the governing equation of E, the spectrum effect is simplified as the characteristic length scale.

When l⁰/l^e = 1, Eq. (10) can be simplified as follows:

(11)

Integrating Eq. (11) from 0 to x gives

(12)

where , and it is the length scale that characterizes the decaying rate of the inlet flow error. In MUAE flows, λ⁰ > 0, and it is constant for most parts of the boundary depth (see Fig. 10 in Subsection 3.3.2). Therefore, D(z) increases with the distance from the ground. It is apparent that E|_x=∞ = 1.0. Therefore, some qualitative conclusions can be made as follows:

(ⅰ) If the turbulent kinetic energy amplitude of the inlet turbulence is larger than that of the fully developed level, i.e., E₁ > 1, the turbulent kinetic energy and E will decay while the changing rate will weaken when z increases.

(ⅱ) If the turbulent kinetic energy amplitude of the inlet turbulence is smaller than that of the fully developed level, i.e., E₁ < 1, the turbulent kinetic energy and E will increase while the changing rate will weaken when z increases.

When E₁=1, i.e., the amplitude of the turbulent kinetic energy is the same as that of the fully developed level while the spectra are different from those of the fully developed level, l^e adjusts to l⁰ because of the rebuild of the turbulent cascade process. For simplicity, the short distance is considered so that the adjustment of l^e is negligible. Thus, l^e ≈ l₁^e(z). Then, integrating Eq. (10) gives

(13)

where , and it is the length scale that characterizes the changing rate of the flow field. Some qualitative conclusions can be made in the region of the inlet immediately as follows:

(ⅰ) If the inflow integral length scale is larger than that of the fully developed level, i.e., l^e/l⁰ > 1 (e.g., the precursor time series are low-pass filtered), the turbulent kinetic energy and E will increase while the changing rate will weaken when z increases.

(ⅱ) If the inflow characteristic length scale is smaller than that of the fully developed level, i.e., l^e/l⁰ < 1 (e.g., the time series are high-pass filtered), the turbulent kinetic energy and E will decrease while the changing rate will weaken when z increases.

The length scale is partly determined by the convection velocity U, which carries the error of the inlet data and partly determined by the gradient of the mean velocity determining the shear production level. Thus, for the flows in which the high shear regions are prominent, the distance needed for the turbulence to develop fully is quite short. Therefore, the nature of the inlet data makes little difference.

3.3.2 Effects of the inlet turbulent statistics

In Case IA1, the amplitude of the precursor time series is scaled such that the Reynolds stress of the transformed data is 1.5 times of the original data and the transformed data are used as the inlet turbulence in Case IA1. The development of the turbulent kinetic energy ratio E²(x, z) in Case IA1 is shown in Fig. 6(a). E²(x, z) decreases when x increases, and the decaying rate weakens when z increases, which is consistent with the conclusions in Subsection 3.3.1. In Fig. 6(b), x is normalized by D(z), and the lines, except for the ones denoting the results at z=1.5h and z=2.0h, show a good collapse. The results confirm that, D(z) is the appropriate length scale which characterizes the development distance in the outer-layer when it comes to the effect of the inlet turbulent kinetic energy amplitude E₁. The line for z=2.0h decays with the same rate as the collapse in Fig. 6(b), and then departs from x/D(z) ≈ 0.2. The line for z = 1.5h decays faster when x/D(z) < 0.2, indicating that, in the inner-layer where turbulent mixing is strong, the inlet error decays faster than the outer-layer. The line for z = 1.5h decays more slowly than the other lines when x/D(z) > 0.2. The departures of the lines for z=1.5h and z=2.0h show that the convection term and diffusion term become more important in the energy budget in the inner-layer and the far field region than those in the outer-layer. Thus, Eq. (12) is not appropriate to some extent.

For Case IS1, the precursor time series is low-pass filtered, and then adjusted to match the Reynolds stress profiles. Thus, the inlet turbulent kinetic energy amplitude is the same as that of the precursor time series, i.e., E₁=1, while the inlet turbulent characteristic length scale is larger than that of the precursor time series, i.e., l^e > l0. The development of E²(x, z) for Case IS1 is shown in Fig. 7(a). E² increases immediately (see Eq. (13)), because the inflow data are in a lack of small-scale turbulent energy, which dominates the dissipation process of turbulent cascade. For the lines denoting z ≤ 4.0h, E² starts to decrease after a peak, and the locations of the peaks in terms of h, i.e., x/h, vary with z. For the two lines, when z > 4.0h, E² is expected to decrease further. In Fig. 7(b), x is normalized by D(z), and the lines collapse to a single line except for the one when z=1.5h. The locations of the peaks in terms of D(z), i.e., x/D(z), are approximately 0.28. Similarly, the results confirm D(z) as the appropriate characteristic length scale of the turbulence-development in the outer-layer when it comes to the effect of the inlet turbulent characteristic length scale l₁^e.

The developments of the lateral energy spectra of u' for Case IS1 are shown in Fig. 8. At the inlet, the energy of the large wavenumber is zero because of the low-pass filter operation and turbulent kinetic energy increases when the small-scale turbulence rebuilds through the turbulent cascade. It is evident that, in Fig. 8(a) (the results of z=0.5h) and Fig. 8(b) (the results of z=1.5h), the cascade reconstruction is fast, and the energy spectra at x=4h are identical with the fully developed turbulence (the results of Case I0). However, for z=4.0h (see Fig. 8(c)), the energy spectra are not fully adjusted until x=12h. These results demonstrate that the turbulent spectra in the inner-layer rebuild faster than those in the outer-layer.

Figure 9 shows the streamwise development of the components of the Reynolds stress for different heights. In Case IA1, both 〈 u'u'〉 and 〈 w'w'〉 decrease when it flows, which is consistent with the turbulent kinetic energy development. For Case IS1, 〈 u'u'〉 starts to increase, and then decreases at a further distance, which is also consistent with the turbulent kinetic energy development. However, 〈 w'w'〉 of Case IS1, especially for the line of z=1.5h, experiences a short-range decrease immediately (see Fig. 9(b)), and then 〈 w'w'〉 starts to increase after a short-range decrease.

In the budget equation of 〈 w'w'〉, the redistribution term, which is the correlation between p' and , is the dominant production term. Since small scale turbulence contributes more to the term than large scale turbulence, the redistribution term is decreased when the inlet turbulent data are low-pass filtered. Thus, 〈 w'w'〉 will experience a decrease.

Figure 10 shows the ratio of each component of the Reynolds stress to the turbulent kinetic energy. For Cases IA1 and IS1, the ratio in the region z < h is almost identical with that of Case I0. However, for Case IS1 in the region h < z < 3.0h, 2〈 u'u'〉/(〈 u'u'〉+〈 v'v'〉+〈 w'w'〉) is larger than those of the other cases, and the remaining terms 2〈 v'v'〉/(〈 u'u'〉+〈 v'v'〉+〈 w'w'〉) and 2〈 w'w'〉/(〈 u'u'〉+〈 v'v'〉+〈 w'w'〉) are smaller than those of the other cases, which is caused by the decreased redistribution term. The departure of the ratio of Case IS1 disappears at x=8.0h. The ratios -2〈 u'w'〉/(〈 u'u'〉+〈 v'v'〉+〈 w'w'〉) of all the three cases are identical with each other, and are almost constant with the height above z=1.5h, indicating that the assumption of Eq. (A10) in Appendix A is reasonable.

The results suggest that, in the simulations of the spatially developing problems, a buffer region is necessary to let the error of the inlet turbulence decay to a low level, and the length of the buffer region needed can be characterized by D(x).

Besides, in the applications of synthetic turbulent generation methods, the turbulent kinetic energy always decays immediately, and a redevelopment process is in need^{[6, 13]}. The decaying can perhaps be offset by increasing the amplitude of the Reynolds stress or the characteristic length scale of the inflow turbulence.

3.4 Effects of the inlet inner-and outer-layer turbulence

The precursor time series is transformed in the inner-layer and outer-layer separately to investigate the effect of the inlet inner-and outer-layer turbulence on the flow field, especially in the inner-layer, which is of great importance in human activities and short-range dispersion.

In Case IA2, the amplitude of the time series is scaled so that the Reynolds stress of the transformed data is 1.5 times of the original data in the outer-layer. In Case IA3, the same procedure is implemented in the inner-layer. The development of E²(x, z) and the normalized turbulent kinetic energy in the inner-layer of Cases IA2 and IA3 are shown in Fig. 11. For the case with the transformed outer-layer data (see Case IA2 in Fig. 11(a)), the error of the turbulent kinetic energy field is small, e.g., E²(x, z) < 110%, in the near field. However, E²(x, z) increases further, and E²(x, z) is larger than 110% in the far field, e.g., x > 5.0h, indicating the region where the inlet outer-layer turbulence dominates. For the case with the transformed inner-layer data (see Case IA3 in Fig. 11(a)), the inner-layer turbulence with high turbulent intensity, e.g., E²(x, z) > 110%, is restricted within x < 7.0h because the extra turbulent kinetic energy decays rapidly in the inner-layer.

The developments of the normalized turbulent kinetic energy profiles are shown in Figs. 11(b) and 11(c). In Fig. 11(b), the vertical profiles are averaged in the spanwise direction. In Fig. 11(c), the data are averaged at the corresponding positions relative to a cube in the lateral direction to obtain the spanwise profiles at z=0.5h. In the vertical profiles, the peak values locate slightly above the top of the building. In the spanwise profiles, the peak values locate at the cube sides. For the case with the transformed inlet inner-layer data, large turbulent kinetic energy errors appear within x < 4.0h. For the case with the transformed inlet outer-layer data, large turbulent kinetic energy errors appear after x > 8.0h. For different inlet turbulence configurations, the effect of the inlet inner-layer turbulence disappears at x=12.0h, and thus the turbulent field is controlled by the inlet outer-layer turbulence.

In Cases IS2 and IS3, the inlet data are obtained by the low-pass filtering precursor time series in the outer-layer and inner-layer, respectively. The development of E²(x, z) and the turbulent kinetic energy in the inner-layer for Cases IS2 and IS3 are shown in Fig. 12. For Case IS2 with the transformed outer-layer data displayed by the contours in Fig. 12(a), the turbulent kinetic energy ratio E² is small in the near field, and increases further when the turbulent kinetic energy ratio E² is larger than 110% in the canopy layer after x=4.0h. For the case with the transformed inner-layer data of Case IS3, the turbulent kinetic energy increases rapidly, and then decays rapidly, such that the region affected by the inlet inner-layer flow (E²(x, z) > 110%) is within the region x < 17.0h, indicating again that the inlet inner-layer turbulence dominates only the near field.

The development of the turbulent kinetic energy profiles is shown in Figs. 12(b) and 12(c). The data are averaged in the same way as that in Figs. 11(b) and 11(c). For different inflow turbulence configurations, at x=16.0h, the effect of the inlet inner-layer turbulence is negligible, and thus the turbulent field is controlled by the inlet outer-layer data.

The lateral spectra of the streamwise velocity u' in the inner-layer of Case IS2 are shown in Fig. 13. For all curves, the spectral energy at low wavenumber is larger than that of Case I0, which can be induced by the high-energy large scale structures in the outer-layer. Accordingly, the turbulent kinetic energy in both the inner-layer and the outer-layer increases when it flows. Since the adjustment of the turbulent spectra in the outer-layer is slow, the large scale turbulence affects the spectral energy in the inner-layer and consequently the turbulent kinetic energy in the inner-layer.

In the rough boundary layer, the turbulent production is active near the top of the canopy, and the turbulent kinetic energy will be partly transported to the outer-layer. If the turbulent kinetic energy increases in the outer-layer somehow as it does in Cases IA2 and IS2 because of the errors of the inflow turbulence, the turbulent kinetic energy diffusion from the inner-layer to the outer-layer will decrease. For a fixed height z, the turbulent flux (〈 u'u'w'〉+〈 vv'w'〉+〈 w'w'w'〉) is the dominant term of the turbulent diffusion flux of the turbulent kinetic energy, and the normalized turbulent flux defined by

is shown in Fig. 14. For Case IA0, the turbulent flux is approximately constant when x increases, while the turbulent flux of Case IA2 changes of the inlet and is smaller than that of Case IA0, i.e., the fully developed case. For Case IS2, the normalized vertical turbulent diffusion flux is smaller than that of Case I0 at the early development stage. Therefore, the turbulent kinetic energy accumulates in the inner-layer, and increases correspondingly.

In conclusion, the error of the inlet inner-layer turbulence dominates the near field because it decays rapidly (see Subsection 3.3). The error of the inlet outer-layer turbulence decays slowly and affects the inner-layer turbulence in the far field through diffusion and convection.

The turbulent statistics in the roughness sublayer are hard to obtain (see Section 1). Fortunately, such uncertainty affects only the near field region in the inner-layer. However, to limit the distance of redevelopment and to ensure the accuracy of the turbulence in the near field, the inlet inner-and outer-layer turbulence should be specified properly.

4 Conclusions

The development of the inlet turbulence in MUAE flows is investigated by transforming the Reynolds stress amplitude and spectra of the precursor time series in the LES of turbulent flow over a cluster of model buildings and with a theoretical analysis.

With the analysis of the turbulent statistics, it is concluded that, (ⅰ) if the characteristic length scale of the inlet data is the same as that of the fully developed level, the error of the turbulent kinetic energy decays until the fully developed level is achieved, and (ⅱ) if the inlet turbulent kinetic energy is the same as that of the fully developed level while the characteristic length scale is larger than the fully developed level, the turbulent kinetic energy will increase immediately and then decrease further because of the adjustment of the turbulent length scale. A length scale, D(z), is deduced by the theoretical analysis, which indicates that the changing rate of the errors of the inlet data in the streamwise direction weakens when the distance from the ground increases.

The detailed numerical testing of the effects of the inlet inner-and outer-layer data on the development shows that the inlet inner-layer data dominate the near turbulent field and the outer-layer inlet data dominate the far field.

Appendix A

The budget of k⁰ for Case I0 is

(A1)

(A2)

(A3)

(A4)

(A5)

Since Case I0 is statistically steady, . Because Case I0 is homogeneous in the horizontal direction (see Subsection 3.1) and W=0, . Therefore, . Since V=0 and W=0, we obtain

Outside the roughness sublayer, the diffusion term D_k is much smaller than the other terms, and thus the dominant balance is between the shear production term and the dissipation, i.e.,

(A6)

Although the diffusion term is negligible in the turbulent kinetic energy budget, it is of significant importance in transporting the turbulent kinetic energy vertically when it comes to the effect of the outer-layer turbulence on the inner-layer turbulence further (see Subsection 3.4).

The budget of k^e for other inlet problems is

(A7)

(A8)

(A9)

(A10)

(A11)

Since the testing cases are statistically steady, . Because V=0 and W=0, we have

The diffusion term is again neglected, and thus the balance is

(A12)

where

(A13)

On the dimensional ground,

Inserting Eq. (A16) into Eq. (A15) yields

(A14)

Define

(A15)

Thus, according to Fig. 10(d),

(A16)

Inserting Eq. (A2) and Eq. (A18) into Eq. (A17), we obtain

(A17)