Numerical Investigation of Very-Large-Scale Motions in a Turbulent Boundary Layer for Different Roughness

Abstract

Wall-model large eddy simulations (WMLES) are conducted to investigate the spatial features of large-scale and very-large-scale motions (LSMs and VLSMs) in turbulent boundary flow in different surface roughnesses at a very high Reynolds number, O (10⁶–10⁷). The results of the simulation of nearly smooth cases display good agreement with field observations and experimental data, both dimensioned using inner and outer variables. Using pre-multiplied spectral analysis, the size of VLSMs can be reduced or even disappear with increasing roughness, which indirectly supports the concept that the bottom-up mechanism is one of the origins of VLSMs. With increases in height, the power of pre-multiplied spectra at both high and low wavenumber regions decreases, which is consistent with most observational and experimental results. Furthermore, we find that the change in the spectrum scaling law from −1 to −5/3 is a gradual process. Due to the limitations of the computational domain and coarse grid that were adopted, some VLSMs and small-scale turbulence are truncated. However, the size of LSMs is fully accounted for. From the perspective of the spatial correlation of the flow field, the structural characteristics of VLSMs under various surface roughnesses, including three-dimensional length scales and inclination angles, are obtained intuitively, and the conclusions are found to be in good agreement with the velocity spectra. Finally, the generation, development and extinction of three-dimensional VLSMs are analyzed by instantaneous flow and vorticity field, and it shows that the instantaneous flow field gives evidence of low-speed streamwise-elongated flow structures with negative streamwise velocity fluctuation component, and which are flanked on each side by similarly high-speed streamwise-elongated flow structures. Moreover, each of the low-speed streamwise-elongated flow structure lies beneath many vortices.

1. Introduction

Large-scale and very-large-scale motions (LSMs and VLSMs) exist in pipe flow, channel flow, and turbulent boundary layer (TBL) flow. LSMs with lengths of $0.3\delta -3\delta$ (where $\delta$ is the boundary layer thickness, pipe radius, or channel half-height) in the streamwise direction has been observed during experiments conducted in turbulent boundary layer flows [1,2], pipe flows [3,4], and channel flows [5,6]. LSMs are packets of hairpin vortices moving upward in the near-wall region. In addition, VLSMs, much longer meandering structures with streamwise lengths larger than $3\delta$ in turbulent boundary layer flows [7,8,9], pipe flows [3,4,10,11], and channel flows [6] were also observed. Both LSMs and VLSMs carry a considerable amount of the turbulent kinetic energy (TKE). In contrast with the case of the LSMs, the mechanism of the generation and development of VLSMs remains an open question.

One widespread proposed explanation is the bottom-up mechanism suggested by Kim and Adrian [3], which implies that VLSMs are formed from coherently aligned bottom hairpin packets. This mechanism is approved by Balakumar and Adrian [2] and Lee and Sung [8] for boundary layer flows, by Guala et al. [4] for pipe flows, and by Balakumar and Adrian [2] for channel flows. However, Hutchins and Marusic [7,12], Mathis et al. [13] and Chung and McKeon [14] observed that the interaction of LSMs and VLSMs was not only a mere superposition on near-wall fluctuations, but also that near-wall small-scale structures were subject to an inhibiting modulation effect from the LSMs and VLSMs in the log region. Therefore, Hutchins and Marusic [7] and Hunt and Morrison [15] proposed a new mechanism for the origin of VLSMs, the top-down mechanism, according to which VLSMs may originate in the upper turbulent boundary layer before moving downward and breaking into small-scale structures in the lower boundary layer as a result of shear forces. As we know that the top-down mechanism is at this point, only conjecture for some experimental observations. VLSMs cannot be distinguished at the top of TBL [2]. In addition, Wang and Zheng [9] found that the product of the wavenumber and the power spectra decreases with height on the right side of the maximum value in the pre-multiplied spectra, which can be explained using the bottom-up mechanism. While, on the left side of the peak, the product of the power spectra and the wavenumber increases with height, indicating that the evolution of the VLSMs may not be entirely attributed to the bottom-up mechanism, and other candidate mechanisms (the top-down mechanism) may, therefore, coexist with the bottom-up mechanism.

The simulation study of the characteristics of VLSMs in TBL flow on a rough surface has not yet been performed. A recent review of studies turbulent flows over rough walls was conducted by Jiménez [16]. Turbulence structures on rough-wall boundary layers and spectrum characteristics in rough-wall pipe flow were investigated by Volino et al. [17] and Rosenberg et al. [18], respectively. In particular, although this was not mentioned by Volino et al. [17], Figure 1 in this paper clearly shows the effect of the roughness wall on pre-multiplied spectra, and this specific phenomenon is consistent with results of this study.

In this study, different surface roughness is selected to indirectly study the origin and the mechanism of development of VLSMs. This is because roughness alters near-wall flow structures. Typically, a spacing of about 100 wall units spanwise and extend from the wall about the same distance for near-wall streaks on smooth boundary layers, as was documented by Kline et al. [19]. This size or larger roughness elements undoubtedly disrupt the streaks.

Studies using direct numerical simulation (DNS) have been recently conducted in an attempt to resolve the temporal evolution of spatial structures fully. Jiménez [20] summed up previous works on the subject before identifying the beginning and end of the $k_x^{-1}$ region in pre-multiplied spectra. He further found the VLSMs in the log region to have a streamwise length of $10\delta -20\delta$ in pipe flow, channel flow, and TBL flow. Del Alamo et al. [21] conducted a channel flow DNS up to friction Reynolds number ${\mathrm{Re}}_{\tau }=1901$ and found evidence for the variant behavior of VLSMs and LSMs in pre-multiplied spectra. Lee and Sung [8] conducted DNS simulations of a TBL flow, with a momentum height Reynolds number of up to ${\mathrm{Re}}_{\theta }=2560$. They showed that the development process of adjacent packet-type structures combines to form VLSMs, and they employed a modified feature-extraction algorithm to identify the properties of the VLSMs. Lee et al. [6] investigated the spatial features of LSMs and VLSMs in a turbulent channel flow with a DNS of ${\mathrm{Re}}_{\tau }=930$, the results indicating that the streamwise length of the VLSMs linearly determines the number of outer LSMs and that the formation of VLSMs possibly comes from the alignment of the positive and negative streamwise-fluctuation structures.

Although DNS is an ideal approach to flow research at low Reynolds numbers, the cost of computation is prohibitively expensive as Reynolds numbers increase. Nowadays, a DNS investigation of these VLSMs at higher Reynolds numbers and with a larger streamwise length domain is not possible. To date, the DNSs of wall-bounded turbulence have been limited to ${\mathrm{Re}}_{\tau }=5000$ [22]. Thus, large-eddy simulation (LES) is very useful for investigating turbulence phenomena in the TBL. This is because LES enables simulations with high Reynolds numbers, and it is computationally more efficient. The horizontally averaged and time-averaged vertical profiles of relevant flow parameters were studied in most LES studies, and hence, they use rather short sizes for the streamwise computational domain [23,24,25,26]. Fang and Porté-Agel [27] first used LES to study the VLSMs in the neutral TBL, followed by Salesky and Anderson [28]. As we know, the LES grid is typically constrained by the near-wall region, where turbulent and viscous momentum transport occurs on small scales, especially where the atmospheric surface is covered with roughness elements that could be more expensive to resolve with the grid. Piomelli and Balaras [29] and Smits and Marusic [22] indicated that for wall-bounded flows, the LES computational costs are estimated as ${\mathrm{Re}}_{\tau }^{9/5}$, only a small savings over DNS, with its ${\mathrm{Re}}_{\tau }^{9/4}$. If, instead, just empirically model the flow that at the near-wall region, the grid requirements would be greatly reduced, which would also significantly reduce computational costs, currently reaching about ${\mathrm{Re}}_{\tau }^{1/5}$. Most important is that wall-model LES (WMLES) tends to be more accurate at high Reynolds numbers than at moderate or low ones [29]. The reason for this is that the WMLES approach depends only on large eddies. Consequently, it can effectively simulate high Reynolds-number flows. Chung and McKeon [14] used the WMLES approach to study large-scale structures in a long channel flow, focusing on their amplitude modulation in near-wall small-scale fluctuation. In the present study, LES coupled with a wall model [30] was indirectly used to investigate the origin of VLSMs in TBL, here the word “indirectly” means we cannot consider the top-down mechanism in these idealized simulations. The present WMLES can produce a promising approach to the study of the characteristics of VLSMs.

In this paper, the VLSMs in the TBL under different surface roughnesses were studied progressively, from a single, two-point correlation field analysis to flow field analysis. Our study is arranged as follows. Section 2 introduces our numerical method, wall model, and subgrid-scale (SGS) model. Section 3 shows the validation of the WMLES approach using DNS and field observation data. The simulation results for VLSMs are presented and analyzed in Section 4. Section 5 provides the conclusion.

2. Numerical Method

2.1. Wall Model LES (WMLES)

For incompressible flow, the time-dependent filtered N-S equations are written as follows:

$$\frac{\partial {\overline{u}}_i}{\partial t}+\frac{\partial }{\partial x_j}({\overline{u}}_j{\overline{u}}_i)=-\frac{\partial }{\partial x_i}\overline{p}-\frac{\partial }{\partial x_j}({\tau }_{ij}^D)+f_i^T$$

(1)

$$\frac{\partial \overline{u_i}}{\partial x_i}=0$$

(2)

where $x_i$ are Cartesian coordinates, $u_i$ are corresponding velocity components; $p$ is pressure divided by fluid density; ${\tau }_{ij}^D$ is the total stress; and $f_i^T$ is the other density-normalized forces. The bar on top of the variables represents a filtering operator.

In Equation (1), on the left hand side, the first term represents the time change rate, the second term represents convection, on the right hand side, the first term gives the density normalized pressure gradient, the second term represents the stresses, and the third term is the other density-normalized forces.

SGS model used to represent the effect of unresolved motions on the resolved motions is critical for accurate simulations, especially for the specific purpose of capturing the dynamics of large-scale structures. In the present study, the one-equation of the eddy viscosity model [31] was adopted. The reason for this being that the crucial feature of the Smagorinsky model lies in the assumption of the balance between SGS energy production and dissipation. On using the ensemble-averaging procedure in place of the filtering in LES, the balance of energy production and dissipation holds only in some restricted flow situations. In fact, this balance breaks down in the vicinity of walls in channel and boundary-layer flows, as in jets and wakes. Moreover, for the Smagorinsky model, the eddy-viscosity representation for the Reynolds stress based on the ensemble mean cannot reproduce the anisotropy of turbulence observed in some typical flows like channel flows. However, the one-equation eddy viscosity model is not based on the assumption of SGS energy balance and uses a more accurate representation for SGS Reynolds stress, the deficiency of local balance assumption adopted in algebraic eddy viscosity models can be overcome. At high Reynolds number flows and for coarse grid resolution such a phenomenon may occur, which fits perfectly with the simulation considered in this study.

The SGS stress can be modeled as:

$${\tau }_{ij,sgs}^D=\frac{2}{3}{k}_{sgs}{\delta }_{ij}-2{\upsilon }_{sgs}{\overline{S}}_{ij}$$

(3)

where ${\tau }_{ij,sgs}^D$ SGS stress part of the total stress ${\tau }_{ij}^D$. ${\overline{S}}_{ij}$ resolved-scale strain rate tensor, ${\delta }_{ij}$ Kronecker delta function, $k_{sgs}$ SGS kinetic energy; ${\upsilon }_{sgs}$ SGS stress eddy viscosity.

Accounting for the historical effect of $k_{sgs}$ because of diffusion, production, and dissipation, a transportation equation is derived:

$$\frac{\partial k_{sgs}}{\partial t}+\frac{\partial {\overline{u}}_j{k}_{sgs}}{\partial x_j}=-{\tau }_{ij,sgs}^D{\overline{S}}_{ij}-C_{\epsilon }\frac{k_{sgs}^{3/2}}{\Delta }+\frac{\partial }{\partial x_j}({\upsilon }_{sgs}\frac{\partial k_{sgs}}{\partial x_j})$$

(4)

where $\Delta$ filter scale, $C_k$ and $C_{\epsilon }$ model constants.

$C_k=0.0673$ and $C_{\epsilon }=0.93$ were adopted in this study.

The computational requirement of high Reynolds number wall-bounded flow of fully-resolved LES dependents strongly on Reynolds number. Piomelli and Balaras [29] pointed out that at high Reynolds number o (10⁶), for smooth boundary layers flow, 99% of the computational mesh grids are used to resolve the inner layer whose thickness is only 10% of the boundary layer height; furthermore, more than 99% of the computational mesh grids is required for the rough boundary layers flow. Consequently, the only economical way to perform LES of high Reynolds number turbulent boundary layer flow is just to compute the outer layer. The outer-flow eddies can determine the grid size in this case, while for the near wall eddies, the grid is too coarse, so must model the wall layer. In particular, near the wall, discrete differentiation cannot evaluate the momentum flux, since the coarse grid cannot resolve the sharp velocity gradient and the quasi-streamwise in this region. Moreover, Piomelli and Balaras concluded that the WMLES tend to be more accurate at high Reynolds number than at moderate or low Reynolds number.

Moreover, for the fully-resolved LES, the SGS model and numerical errors very large in the near-wall region usually; thus, it is conceivable that comparison with fully-resolved LES, using “perfect” wall model can get a more accurate predictions.

Another important point we need to consider involves the case of an atmospheric surface covered with high roughness elements (e.g., sands, rocks, vegetation) that would greatly increase the computation cost when used in conjunction with the fully-resolved LES.

Thus, in the present study, the surface stress model is the following matrix:

$${\tau }_{ij}^D=\left(\begin{array}{ccc}0& 0& {\tau }_{13}^{tot}\\ 0& 0& {\tau }_{23}^{tot}\\ {\tau }_{13}^{tot}& {\tau }_{23}^{tot}& 0\end{array}\right)$$

(5)

and, the Schumann’s model [30] was adopted:

$${\tau }_{13}^{tot}=-u_{*}^2\frac{({\overline{u}}_{1/2}-\langle {\overline{u}}_{1/2}\rangle )}{{({\langle {\overline{u}}_{1/2}\rangle }^2+{\langle {\overline{v}}_{1/2}\rangle }^2)}^{1/2}}$$

(6)

$${\tau }_{23}^{tot}=-u_{*}^2\frac{({\overline{v}}_{1/2}-\langle {\overline{v}}_{1/2}\rangle )}{{({\langle {\overline{u}}_{1/2}\rangle }^2+{\langle {\overline{v}}_{1/2}\rangle }^2)}^{1/2}}$$

(7)

where 1/2 is at first level cell centers near the surface; the angle brackets represent horizontal average at a height.

The friction velocity:

$$u_{*}^2=({\langle {\tau }_{13}^{tot}\rangle }^2+{\langle {\tau }_{23}^{tot}\rangle }^2{)}^{1/2}$$

(8)

which rough wall log law:

$$\frac{{(\langle {\overline{u}}_{1/2}\rangle +\langle {\overline{v}}_{1/2}\rangle )}^{1/2}}{u_{*}}=\frac{1}{k}\ln (\frac{z}{z_0}+f(L))$$

(9)

where $f(L)$ atmospheric stability-related function, $L$ Obuhkov length, $z_0$ aerodynamics roughness height.

In this study, $f(L)$ is zero for neutral stability, $z_0$ is related to roughness elements shape, height, and distribution.

In this paper, the five different cases of surface roughness are investigated: $z_0$ = 0.0002 m (case 1, ${\mathrm{Re}}_{\tau }$ = 3.538 × 10⁶), $z_0$ = 0.002 m (case 2, ${\mathrm{Re}}_{\tau }$ = 4.894 × 10⁶), $z_0$ = 0.02 m (case 3, ${\mathrm{Re}}_{\tau }$ = 6.559 × 10⁶), $z_0$ = 0.2 m (case 4, ${\mathrm{Re}}_{\tau }$ = 9.927 × 10⁶), and $z_0$ = 2.0 m (case 5, ${\mathrm{Re}}_{\tau }$ = 1.199 × 10⁷).

2.2. Simulation Setup

Figure 1 gives the schematic diagram of the computational domain (5000 × 1000 × 300 m for case 1, 5000 × 1000 × 350 m for case 2, 5000 × 1000 × 400 m for case 3, 5000 × 1000 × 500 m for case 4, 5000 × 1000 × 600 m for case 5) adopted in this study. The streamwise, spanwise, and wall-normal coordinates are x, y, and z. At the solid wall, wall model was adopted, and at the upper surface, slip boundary conditions were adopted. For the streamwise and spanwise directions, periodic boundary conditions were applied, same with the LES TBL flow simulation by Fang and Porté-Agel [27], and DNS channel flow simulation conducted by Lee et al. [6]. The logarithmic mean wind-speed profile was used for the inflow, specifically for five cases: Case 1, 250 m, 8 m/s; case 2, 300 m, 8 m/s; case 3, 350 m, 8 m/s; case 4, 450 m, 8 m/s; and case 5, 550 m, 8 m/s.

The hexahedral meshes type was adopted in the simulation. The computational domain at the height below 200 m, the mesh grid size is 5 m in each direction. Except for this region, the mesh grid size is 10 m in each direction.

The time step is set to 0.5 s in order to satisfy the Courant number for the stability condition and viscous stability. To guarantee a fully developed flow field and quasi-steady conditions, a spin-up simulation is first made for enough time. To ensure statistical convergence, statistical calculations are then performed over 40,000 timesteps. At 0.015% of domain’s height, the initial peak in perturbation magnitude occurs, with both the maximum perturbations of streamwise/spanwise flow near the surface running at 0.25 m/s, a 4.0/20.0 ratio was adopted at the beginning of the simulation for total periods of streamwise/spanwise perturbations.

The simulations were performed using the Simulator for Wind Farm Application (SOWFA) [32], which is developed using OpenFOAM C++ library. The numerical discretization was performed by a finite volume method. A second-order central difference scheme was used for the convection and viscosity terms, and a second-order backward scheme was employed for the unsteady term.

3. Validations

Because the wall layer is modeled in WMLES, the accuracy of the simulation results must be validated before the study can proceed. Here, the simulation results based on non-dimensional inner variables $z^{+}=z{u}_{\tau }/v$ and $u_{\tau }$, and outer variables $\delta$ and $U$, are compared with the results of DNS, experiments, and field measurements.

3.1. Validation by Non-Dimensional Inner Variables

To confirm that the WMLES approach is suitable for TBL flows at high Reynolds numbers, variations in mean velocity, Reynolds stress, and turbulent variation with height and Reynolds number are compared with theoretical results and experimental, field-measured data for the TBL flow.

The non-dimensional mean streamwise velocity profiles, shown in Figure 2. The simulated non-dimensional mean velocities $U^{+}=U/u_{\tau }$ are displayed at different heights from case 1. The mean velocity profile agrees well with the logarithmic law and Qingtu Lake Observation Array (QLOA) field-measured data, indicating that the WMLES is highly accurate for TBL flow, despite the down shift in ΔU⁺ due to ground roughness.

Figure 3 shows the non-dimensional Reynolds shear-stress profiles for case 1. These agree well with theoretical predictions [34,35] for canonical TBL flow, the field-observation data of Hutchins et al. [36] at ${\mathrm{Re}}_{\tau }$ = 7.7 × 10⁵, and the QLOA field measured data with higher Reynolds number of Wang and Zheng [9].

In addition, $\overline{uu}/u_{\tau }^2$ is plotted against the inner-scaled height, as shown in Figure 4, including the field measurement data from the SLTEST site [36,38], laboratory data [39,40], and the similarity formulation presented in Marusic and Kunkel [41]. This reasonable agreement suggests that the simulated data agree well with the theoretical formulation [34] and previous measurements at the SLTEST. The magnitude of the normalized streamwise component of normal stress is strongly dependent on the Reynolds number, except for the sublayer. Specifically, the normalized streamwise component of normal stress, simulated at $z/\delta \approx ~0.036$, was compared to the results of QLOA [9], SLTEST [36], and wind tunnel experiments [40,42,43,44] in Figure 5. The current simulation result falls within the QLOA results because the Reynolds number for case 1 is located within the five cases of Wang and Zheng [9], indicating an approximate log-linear relationship between $\overline{uu}/u_{\tau }^2$ and the Reynolds number.

3.2. Validation by Non-Dimensional Outer Variables

In the outer layer of TBL, the flow is assumed to be independent of viscosity and to depend on the global characteristics of the flow, as represented by $\delta$, $U$.

The Reynolds shear stress is shown in Figure 6 against the outer-scaled coordinates $z/\delta$. It is clear that the simulated result of case 1 is in a good agreement with the experimental results of Graaff and Eaton [37]; most importantly, around $z/\delta =0.02$, a trend in the shift of Reynolds shear stress is also captured by this simulation, which indicates the accuracy and high precision of the WMLES adopted in this study.

Normalized streamwise component of normal stress is shown in Figure 7 against outer-scaled coordinates $z/\delta$. The simulated result of case 1 agrees well with the experimental results of Graaff and Eaton [37]; like the Reynolds shear stress, a normalized streamwise component of normal stress also shows a shift in the trend with increasing height.

Figure 8 shows the mean velocity defect profiles according to the outer coordinates. The profile of case 1 is in good agreement with the experiment results of Graaff and Eaton [37] and the DNS results of Lee and Sung [8] throughout the outer region of the boundary layer, albeit at different Reynolds numbers.

4. Results and Discussion

4.1. Velocity Spectra Characteristics under Different Surface Roughnesses

4.1.1. Streamwise Velocity Spectra

Figure 9 shows the wavenumber-normalized spectra of streamwise velocity $S_u/u_{\tau }^2\delta$ for all five cases at different heights. Taylor’s hypothesis of spatial–temporal transformations was adopted to calculate the wavenumber spectrum, to smooth the spectra, a multi-point averaging method was used, following Balasubramaniam [45]. As Figure 9a–c shows, for $k\delta <2$ (wavelengths larger than $3\delta$), each spectral curve has an individual inflection point where the slope of the curves conspicuously changes. The inflection point represents a spectral peak in low-wavenumber region of the pre-multiplied spectra (as shown in Figure 10) associated with the VLSMs [3,9,46]. Thus, it can also be demonstrated that the VLSMs exist in the first three cases for the current TBL study. While, for cases 4 and 5, individual inflection points exist for $k\delta >2$, which means that VLSMs does not exist in these two cases. Above all, with where roughness increases, the scale of the VLSMs (the width of −1 region) decreases or even disappears. We predict that the reason for this may be associated with roughness elements of similar or larger sizes, which undoubtedly disrupt near-wall quasi-streamwise streaks. This phenomenon is identical to what is shown in Figure 1 from Volino et al. [17], although it is not mentioned as such in the report of that study.

Figure 10 shows the pre-multiplied spectra for all five roughness cases. There is a distinct peak in the spectra for all the heights in Figure 10a–c in the low-wavenumber region that corresponds to VLSMs; another peak appears in high-wavenumber region that corresponds to the LSMs in the spectra for all heights in Figure 10a–e.

In addition, the product of wavenumber and power spectra $k{S}_u$ decreases with height throughout all of the spectra under various ground roughness, which is consistent with the result of the turbulent-boundary layer [2,46], and the pipe-flow result [4]. This phenomenon can be explained by “bottom-up” mechanism, and VLSMs weaken in the wall-normal direction, due to the TKE produced in near-wall region gradually dissipates.

The slope between the first and second peaks of the streamwise velocity spectra is also analyzed. This is governed by the transformation scaling law of VLSMs and LSMs (cases 1, 2, and 3) or the transformation scaling law of LSMs and the small-scale turbulence vortex (cases 4 and 5); instead of observing jumps from −1 to −5/3, we observe a gradual change as height increases. With increasing roughness, the scaling law grows, and the transformation rate from −1 to −5/3 increases, with the roughness increasing as well, while the VLSMs go extinct. Another interesting phenomenon observed here is the scaling law in cases 4 and 5, showing height increases of −5/3. For other cases, the scaling law is located between the values −1 and −5/3, with the smoother terrain yielding the larger scaling law.

Based on this finding, a new concept could be proposed, according to which the presence or absence of VLSMs cannot simply be defined by the value −1 or −5/3; instead, the scale factor $\lambda /\delta$ (the length of VLSMs and LSMs to the boundary-layer height) and the scaling law should be combined to illustrate the location of the vortex structure in VLSMs, LSMs, or small-scale turbulence vortex regions.

The $k\delta$ decreases at the first peak (meaning ${(k\delta )}_0$) with increasing height; when the height reaches to a certain extent, ${(k\delta )}_0$ begins to attain a constant value or even to increase, as determined by the characteristics of the VLSMs and is consistent with the literature [18,46]. These two phenomena can be seen more clearly in Figure 11, which plots the relationship between ${(k\delta )}_0$ and the dimensionless height $z/\delta$.

In Figure 11, the result of case 1 agrees with the field-measured data, as a result of the similarity of terrain between this and the QLOA site, as shown in Figure 2. With increases in roughness, cases 2 and 3 reach a certain degree of agreement with field data. While, for cases 4 and 5, the results diverge from the field data, implying that the surface roughness has a great influence on the structure of VLSMs. As we know, there is also a certain deviation between the experimental value and the field-measured value, which may be due to the difference in Reynolds number between them.

4.1.2. Wall-Normal Velocity Spectra Characteristics

Wall-normal velocity fluctuation spectra are shown in Figure 12. The trend for wall-normal location is similar to the streamwise velocity fluctuations, with the exception of the crossover phenomenon, which is consistent with the finding of Guala et al. [4]. For the left part of the crossover point, energy density increases from the wall to a wall-normal location, and then it decreases again; for the right part of the crossover point, the energy density always decreases. In addition, it is clear that as the roughness increases, the crossover point moves to a large wavenumber, which indicates that the rougher the terrain, the less the vortex interacts at different heights. Comparison with streamwise velocity spectra, the wall-normal velocity spectra has lower values, which reflects the difference between the two components in the total energy, and the ultimate difference is close to one order of magnitude.

4.2. Effects of VLSMs and LSMs on Turbulent Kinetic Energy

One of the principal questions addressed in this study is how far VLSMs and LSMs contribute to TKE. We will examine the distribution of TKE as a function of height.

The TKE of VLSMs and the corresponding energy fraction are plotted against height using a low-pass filter with a cutoff length of $3\delta$, as shown in Figure 13a,b, respectively. Similarly, band-pass filtering with cutoff lengths of $0.3\delta -3\delta$ and $\delta -3\delta$ are used to demonstrate the energy fraction of the LSMs, as shown in Figure 13c,d, respectively. The total kinetic energy of VLSMs of both case 1 in the present study and the experimental data of Balakumar and Adrian [2] slightly decreases with increases in height. However, the total kinetic energy of VLSMs increases with height increases as in Wang and Zheng [9], which may relate to fact that the case investigated by Wang and Zheng [9] found energy increases with height in the low-wavenumber region of the wind spectra, unlike what was found in other studies [2,46]. Case 1 shows a close agreement with the experimental data of Balakumar and Adrian [2], given the increasing ground roughness and nondimensional decreases in TKE. The reason for this is found in the fact that VLSMs decrease or even disappear as roughness increases. It is worth noticing that Figure 13 also shows a result found in cases 4 and 5 that does not exist in VLSMs. It simply is the mathematical expression of the region of VLSMs and LSMs, as shown in Figure 10.

In addition, the energy fraction of the VLSMs in the present study increases with height, which tends to close agreement with the wind-tunnel measurements of Balakumar and Adrian [2], as shown in Figure 13b; however, little deviation from the field-measured data of Wang and Zheng [9] appears, due to the spectrum characteristics mentioned above. The phenomenon, shown in Figure 13 indicates, that contrary to most other studies [4,8,47], which have generally found that VLSMs contribute 40–50% of the total TKE, the present study shows that VLSMs contribute up to 40% only at the highest height, shown in Figure 13b. In Figure 13c, it is seen that in LSMs ($0.3\delta <\lambda <3\delta$, Guala et al. [4]), the energy fraction is larger than the field-measured data; but in Figure 13d, the energy fraction of the LSMs ($\delta <\lambda <3\delta$, Lee et al. [48]) is similar to the field measured data. It is clear that the energy fraction of the LSMs is independent of ground roughness, possibly owing to the fact that the size of the element roughness was adopted in the present study only to disrupt the VLSMs, yielding no effect on LSMs, at least from the viewpoint of energy fraction.

The reason why the energy fraction of LSMs is found to be large in our study can be explained as follows. Because the WMLES approach is used, the smallest scale falls in the inertial subrange, while some or even most of the inertial subrange remains unresolved. It can be seen from the wind spectra shown in Figure 10 that it quickly falls into the dissipative region. Therefore, the total energy in this study is small. Due to the relative limitations of the computational domain, the structure of VLSMs is not completely contained, so that the simulation of the energy of VLSMs is less than that of the real atmosphere or the sufficient computational-domain case. This might also explain why the experimental data of Balakumar and Adrian [2] are lower than the field-measured of from Wang and Zheng [9], as shown in Figure 13b, the experimental domain being also limited to some extent. The scale of LSMs, ranging between $0.3\delta -3\delta$, can be fully resolved. In summary, the energy of VLSMs is relatively smaller than the field-measured data, while the energy of LSMs is larger.

4.3. Flow Field Correlation Analysis

To intuitively visualize the structure of VLSMs, the coefficients of spatial correlation of the streamwise fluctuating velocity between the spatial and reference locations ($(x,y,z)=(2500m,500m,35m)$) were calculated. Figure 14 shows the contours of the spatial correlation coefficient for five types of roughness cases. The contour level of the correlation coefficient indicating the coherent structure edge was set to be 0.05, following Hutchins and Marusic [7]. As the figure shows, the length scale of the streamwise coherent structure satisfies $\lambda /\delta >3$ for the first three cases, indicating the existence of the VLSMs. For cases 4 and 5, the $0.3<\lambda /\delta <3$ range only bears on the existence of the LSMs. For increasing roughness, the scale of VLSMs and LSMs decreases, indicating that the generation and development of VLSMs and LSMs are affected by the bottom terrain. It also indirectly illustrates that the bottom-up mechanism controls the generation and development of VLSMs.

Along the streamwise direction, the correlation contour lines are tilted, which indicates the inclination of turbulent structures that observed by Marusic and Heuer [49] and Wang and Zheng [9] is also present in the current TBL simulation at higher Reynolds numbers. The inclination angle of the two-point correlation of the fluctuating velocity is related to the average extent of the inclination of the structures of the VLSMs or LSMs. For the present simulation, the inclination angles in Figure 14 are 11.31° for case 1 (${\mathrm{Re}}_{\tau }$ = 3.538 × 10⁶), 14.04° for case 2 (${\mathrm{Re}}_{\tau }$ = 4.894 × 10⁶), 16.70° for case 3 (${\mathrm{Re}}_{\tau }$ = 6.559 × 10⁶), 14.04° for case 4 (${\mathrm{Re}}_{\tau }$ = 9.927 × 10⁶), and 19.29° for case 5 (${\mathrm{Re}}_{\tau }$ = 1.199 × 10⁷). The results of case 1 agree closely with those found by Wang and Zheng [9], who indicated that the inclination angles of field-measured data were 12.14° for ${\mathrm{Re}}_{\tau }$ = 3.14 × 10⁶, while for ${\mathrm{Re}}_{\tau }$ = 1.63 × 10⁶, the angle is 20.35°. The difference among all five cases is comparable to the range reported in the literature for smooth- and rough-wall boundary layers. Christensen and Wu [50] found an inclination angle of 11° for their smooth-wall channel flow. Head and Bandyopadhyay [51] observed inclination angles between 15–20°. Christensen and Adrian [52] reported 12–13°. Adrian et al. [53] found an inclination angle of 12°. Tomkins and Adrian [54] found inclination angles between 10–20°. Nakagawa and Hanratty [55] found an angle of 9° on the rough wall in a channel flow, compared to smooth-wall values of 6–8°. In general, the inclination angle of large-scale coherent structures increases with roughness.

In addition, the streamwise, wall normal, and spanwise extent of fluctuating-velocity two-point correlations are shown in Figure 15. The distance, $L{x}_{uu}$, twice the distance from the most downstream location to the self-correlation peak on two-point correlations of the fluctuating streamwise velocity equal to 0.5 contour, as defined in Christensen and Wu [50]. Volino et al. [17] found that the results for rough and smooth walls agree well with a value of about $L{x}_{uu}/\delta =0.65$ between $z/\delta =0.1$ and 0.6 with their experimental data. Close to the wall, some differences are visible. In Figure 15a, the result of case 1 is much smaller than that found by Volino et al. [17]; however, the field observation data of Wang and Zheng [9] with solid symbols are shown in Figure 15a, for ${\mathrm{Re}}_{\tau }$ = 3.14 × 10⁶ measured data, almost same as in case 1 (${\mathrm{Re}}_{\tau }$ = 3.553 × 10⁶). This indicates that the results of simulation are correct. This can be explained by the fact that field-measured and simulation data have a similarity ratio comparable to the experimental data. In Figure 15a, as roughness increases, $L{x}_{uu}/\delta$ decreases, which agrees with the results of Krogstad and Antonia [56], who found $L{x}_{uu}/\delta$ was about 50% lower on their rough wall. The wall-normal extent of the two-point correlations of the fluctuating velocity, $L{z}_{uu}$, is determined based on the wall-normal distance between the points farthest and closest from the surface at a 0.5 contour. As the contours merge with the wall, reliable estimates for $L{z}_{uu}$ cannot be obtained for $z/\delta <0.2$. The conclusion is similar to the $L{x}_{uu}$, as shown in Figure 15b. The spanwise extent $L{y}_{uu}$ shown in Figure 15c, with increasing roughness, the $L{y}_{uu}/\delta$ range from 0.2 to 0.1, which agrees with the DNS result of Lee and Sung [8], who concluded that the streamwise negative motions meander in the downstream direction with a characteristic width of approximately $0.1\delta -0.2\delta$ in the spanwise direction. Comparing Figure 15a,b, the ratio of $L{x}_{uu}/L{z}_{uu}$ obtained in Figure 15d, the value is roughly 2.5 for both the experimental data of Volino et al. [17] and the present simulation cases, also consistent with Nakagawa and Hanratty [55], and the smooth-wall results of Krogstad and Antonia [56].

4.4. Instantaneous Flow Field Analysis

Figure 16 shows visualizations of the VLSMs or LSMs structures in the instantaneous flow fields for $z$ = 30 m in the $xy$ plane, both the VLSMs and LSMs were visualized based on the streamwise-elongated negative velocity, which is consistent with the results of previous studies [6,8,57,58,59]. The direction of flow is from left to right. Figure 16a–c show several apparent very long negative (dark color) motions meandering in the streamwise direction, which in the spanwise direction are flanked by positive streamwise fluctuation (white color), and these motions often extend more than $10\delta$ in streamwise direction, and the widths is around $0.1\delta -0.2\delta$. This phenomenon is in agreement with TBL and turbulent channel DNS studies [7,8]. For Figure 16d,e, there are some intermediate long negative streamwise structures, and the scale characteristics are similar to the LSMs structures.

Figure 17 shows visualizations of the VLSMs or LSMs structures in the instantaneous flow fields for line A, B and C of Figure 16c in the $xz$ plane, and the dashed line in Figure 17 corresponds to the horizontal plane height of Figure 16. In Figure 16c, line A locates at the very long streamwise structure, line B locates at the intermediate long streamwise structure in order to see the inclination angle of VLSMs, and line C locates at the high-velocity region. Comparing these three lines streamwise fluctuation in Figure 17, there are inclined at angles of 18.26° of line B, very close to the spatial correlation result, 16.70°, which is calculated by the reference point of streamwise direction. Note that the inclination angle is 10°–15° in an almost linear type [8], it also can be seen in the present study, the line A has nearly same inclination angle. Moreover, the flow patterns of these aligned packets are consistent with the observations [53]. For line C, which locates at the high-velocity region, the inclination angles at A and B do not exist at line C, this indicates that LSMs and VLSMs structures located at low momentum region, and comes from the near wall aligned packets, as described by “bottom-up” mechanism.

4.5. The Generation and Development of VLSMs

The results of Section 4.4 illustrate the presence of VLSMs and their associate with low-momentum regions. In this section, we address the creation and extinction procession of VLSMs in the flow field. As we know, the hairpin packet model could give a reasonable explanation for logarithmic layer features, for example, quasi-streamwise vortices and bursting processes. An idealized conceptual model Based on the hairpin packet paradigm and in the outer layer, most features of the coherent structures, an idealized conceptual model is built [53]. In contrast, Kim and Adrian [3] suggested that VLSMs is not a new type of turbulent motion. They conjectured that the alignments of hairpin packets could form VLSMs; however, the created and extincted of VLSMs in TBL is little to know. Thus, to obtain spatial information about the formation and disappearance of VLSMs is our object.

The mechanism of Kim and Adrian [3] means that in the streamwise direction continually spawning new hairpins to create packets. After growing, the packets become larger and longer, and then probably merge with adjacent packets or broken up into small structures. In this study, a representative case of merging and extincting of vortex packets is discussed in detail, which VLSMs showed in the dashed red rectangle of Figure 16c, The time evolution of flow field in the $xy$ plane with a time interval of 20 s and the height is 30 m, shown in Figure 18. The dashed lines D, E, F and G were added to show the development of the VLSMs. It is obvious that the VLSMs in Figure 18f originates in several packet motions P1, P2, P3, P4, P5, P6, P7, and P8 present in Figure 18a. For dashed line D1 in Figure 18a, two adjacent low-momentum regions of P1 and P2 are travelling downstream, then they merge together from Figure 18f, the similar phenomenon is observed at dashed line D2. For dashed line E, the P2 and P3 are merge together, as they move downstream, and the low-momentum regions separate in the streamwise direction, as shown in Figure 18j. For dashed line F1, two adjacent low-momentum regions of P3 and P4 are travelling downstream, and the first stage is that they merge together from Figure 18a–e, then separate from each other again in the Figure 18j, the similar phenomenon is observed at dashed line F2. For dashed line G1, as part P5 and P6 move downstream, and the first stage is that they separate from Figure 18c and then merge together again in the Figure 18f, the similar phenomenon is observed at dashed line G2. The four cases of D1, D2, E, F1, F2 and G, represent three typical development process of adjacent hairpin vortices. Due to the low-momentum regions merge in streamwise direction, resulting in a longer structure at the streamwise scale of the flow field in Figure 18f,g. Another phenomenon is that although the streamwise growth of the VLSMs occurs; however, the spanwise length of VLSMs becomes thinner slightly.

Figure 19 and Figure 20 show the associated flow features in the $xz$ plane. The conclusion is consistent with the results show in Figure 18. Combining the Figure 18, Figure 19 and Figure 20, which indicates that the development and extinction of VLSMs structure is a three-dimensional process. As can be seen in these figures, in the downstream direction, these packets are growing continuously and further away from the wall, part of the upstream packets go through the higher mean flow velocity. The lower parts slower than these part flows, which leads lift vortex away from the wall by a higher velocity, and in the higher part of the upstream packet causes in greater stretching. Consequently, in downstream the adjacent packets are merged to form VLSMs with a shallow angle. These processes are same with hairpin structure dynamics [5] and DNS results [8].

For the formation of VLSMs, present simulations describe one possible mechanism, however, in the future, the analysis of the general mechanism of VLSMs formation is needed, because the trigger of the formation of new packets still no definitive evidence.

4.6. Vortices Field Analysis

The packet merging process is presented schematically in Figure 17 of Tomkins and Adrian [54]. The similar conclusion was obtained in figure 25 of Adrian, Meinhart and Tomkins [53]. This is the reason for the vorticity field analysis, which could show the VLSMs structure is formed by hairpin vortex or hairpin vortex package development.

Figure 21 shows three-dimensional vortical structures that contoured by ${\lambda }_2$ [60] in the horizontal plane. It is clear that several highly elongated low-speed regions are shown in the flow field, and many hairpin-type vortices lie above each of them. Corresponding to the Figure 17, here Figure 22 is shown as the vertical plane to illustrate the negative velocity region is covered by vortices, while positive velocity is not. By comparing A, B, and C, it is clearly shown that vortices are more density for A than B than C, which locates in high-speed region—this is consistent with the preceding results. These flow patterns, together with three-dimensional vortical structures to form longer structures that erupt from the wall and grow towards the outer region, are consist with Lee and Sung [8].

5. Conclusions

In this paper, the spatial features for LSMs and VLSMs in the TBL flow under different surface roughness at very high Reynolds number, O (10⁶–10⁷), were investigated using the WMLES approach. The VLSMs in TBL were studied progressively from a single point, in correlation with a field analysis. The statistical properties of the mean and fluctuating velocity for the nearly smooth case were compared with those found in the literature, including field observations, laboratory studies, and DNSs, and a close agreement was found using the inner and outer variables. These results, thus, provide a useful database for the turbulence statistics of TBL.

Our analysis of the streamwise velocity spectra demonstrates the tendency is that with increasing roughness, the scale of VLSMs decreases or even disappears. The ${(k\delta )}_0$ decreases with height in an approximately log-linear manner. While, when height increases to a certain extent, ${(k\delta )}_0$ begins to reach a constant value or even to increase. In addition, $k{S}_u$ decreases with height throughout all spectra.

For wall-normal velocity spectra, the crossover phenomenon occurs, and as roughness increases, the crossover point moves to larger wavenumber.

The VLSMs and LSMs carry a considerable portion of the TKE. While, the energy fraction of VLSMs in the present study is a little smaller than the field-measured data, but it is similar to experimental data, which can be explained by spectral characteristics.

The spatial correlation coefficient pertaining to five different roughness cases were calculated. The results have reached a close agreement with the wind spectra. For the ratio of the streamwise extent and wall-normal extent of the two-point correlations of the fluctuating velocity, the near-constant is 2.5. Moreover, the inclination angle of nearly smooth cases is consistent with field observation data at almost the same Reynolds number. The tendency of the inclination angle of different cases of roughness is similar to what has been shown in other studies.

Then, we conducted the instantaneous flow field analysis to illustrate the generation and development of VLSMs. It is clearly shown that both the VLSMs and LSMs were visualized based on the streamwise-elongated volumes associated with the negative streamwise velocity, and which are flanked on either side by similarly elongated positive streamwise velocity. In addition, the generation and extinction process of VLSMs is seen by three-dimensional flow field analysis.

In the end, Vortices field analysis is conducted to explain further the VLSMs locate in negative streamwise velocity region and each of them lies beneath many hairpin-type vortices.

Above all, in this study, the results at least indirectly support that the bottom-up mechanism is one of the origins of VLSMs. It is confirmed here that the WMLES method has the ability to investigate the spatial features of LSMs and VLSMs in the TBL flow at high Reynolds numbers at neutrally stratified TBL, and the applicability of these results to a real atmospheric boundary layer flow that stability variation will be investigated in the future study.