**Funding:** This research was funded by the National Natural Science Foundation of China (number 51804318 and number 51805188), and the National Key Research and Development Program of China (number 2018YFC0808401).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Numerical Investigation of Very-Large-Scale Motions in a Turbulent Boundary Layer for Di**ff**erent Roughness**

#### **Hehe Ren 1,2,3, Shujin Laima 1,2,3,\* and Hui Li 1,2,3**

Received: 8 January 2020; Accepted: 26 January 2020; Published: 4 February 2020

**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* (106–107). 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.

**Keywords:** WMLES; VLSMs; LSMs; turbulent boundary flow; roughness

#### **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δ − 3δ (where δ 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δ in turbulent boundary layer flows [7–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*−1 *x* region in pre-multiplied spectra. He further found the VLSMs in the log region to have a streamwise length of 10δ − 20δ in pipe flow, channel flow, and TBL flow. Del Alamo et al. [21] conducted a channel flow DNS up to friction Reynolds number Reτ = 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 Reθ = 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 Reτ = 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 Reτ = 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–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 Re9/5 τ , only a small savings over DNS, with its Re9/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 Re1/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 \mathbf{x}\_j} (\overline{u}\_j \overline{u}\_i) = -\frac{\partial}{\partial \mathbf{x}\_i} \overline{p} - \frac{\partial}{\partial \mathbf{x}\_j} (\tau\_{ij}^D) + f\_i^T \tag{1}$$

$$\frac{\partial \overline{u\_i}}{\partial \mathbf{x}\_i} = 0 \tag{2}$$

where *xi* are Cartesian coordinates, *ui* are corresponding velocity components; *p* is pressure divided by fluid density; τ*D ij* is the total stress; and *f T i* 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^{D}\_{ij,\text{sys}} = \frac{2}{3} k\_{\text{sys}} \delta\_{ij} - 2 \nu\_{\text{sys}} \overline{\mathfrak{S}}\_{ij} \tag{3}
$$

where τ*D ij*,*sgs* SGS stress part of the total stress τ*D ij* . *Sij* resolved-scale strain rate tensor, δ*ij* Kronecker delta function, *ksgs* SGS kinetic energy; υ*sgs* SGS stress eddy viscosity.

Accounting for the historical effect of *ksgs* because of diffusion, production, and dissipation, a transportation equation is derived:

$$\frac{\partial \mathbf{k}\_{\rm sgs}}{\partial t} + \frac{\partial \overline{\mathbf{u}}\_{j} \mathbf{k}\_{\rm sgs}}{\partial \mathbf{x}\_{j}} = -\boldsymbol{\pi}\_{ij, \rm sgs}^{D} \overline{\mathbf{S}}\_{ij} - \mathbf{C}\_{\varepsilon} \frac{\mathbf{k}\_{\rm sgs}^{3/2}}{\Delta} + \frac{\partial}{\partial \mathbf{x}\_{j}} (\boldsymbol{\upsilon}\_{\rm sgs} \frac{\partial \mathbf{k}\_{\rm sgs}}{\partial \mathbf{x}\_{j}}) \tag{4}$$

where Δ filter scale, *Ck* and *C*ε model constants.

*Ck* = 0.0673 and *C*ε = 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* (106), 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 = \begin{pmatrix}
0 & 0 & \tau\_{13}^{tot} \\
0 & 0 & \tau\_{23}^{tot} \\
\tau\_{13}^{tot} & \tau\_{23}^{tot} & 0
\end{pmatrix} \tag{5}
$$

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

$$
\pi\_{13}^{\text{tot}} = \left. -u\_\*^2 \frac{\left(\overline{u}\_{1/2} - \left<\overline{u}\_{1/2}\right>\right)}{\left(\left<\overline{u}\_{1/2}\right>^2 + \left<\overline{v}\_{1/2}\right>^2\right)^{1/2}}\right. \tag{6}
$$

$$\pi\_{23}^{\text{tot}} = -u\_\*^2 \frac{\left(\overline{v}\_{1/2} - \left\langle \overline{v}\_{1/2} \right\rangle\right)}{\left(\left\langle \overline{u}\_{1/2} \right\rangle^2 + \left\langle \overline{v}\_{1/2} \right\rangle^2\right)^{1/2}} \tag{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:

$$
\mu\_\*^2 = \left( \left< \tau\_{13}^{tot} \right>^2 + \left< \tau\_{23}^{tot} \right>^2 \right)^{1/2} \tag{8}
$$

which rough wall log law:

$$\frac{\left(\left<\overline{u}\_{1/2}\right> + \left<\overline{v}\_{1/2}\right>\right)^{1/2}}{u\_\*} = \frac{1}{k} \ln(\frac{z}{z\_0} + f(L))\tag{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, Reτ = 3.538 × 106), *z*0 = 0.002 m (case 2, Reτ = 4.894 × 106), *z*0 = 0.02 m (case 3, Reτ = 6.559 × 106), *z*0 = 0.2 m (case 4, Reτ = 9.927 × 106), and *z*0 = 2.0 m (case 5, Reτ = 1.199 × 107).

#### *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.

**Figure 1.** Schematic diagram of the computational domain.

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*+ = *zu*τ/*v* and *u*τ, and outer variables δ 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*τ 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 2.** Comparison of mean velocity profiles. The solid colored symbols indicate current simulation data; the black dashed line indicates *U*+ = ln(*z*+)/0.41 + 5.0 for a smooth-wall wind profile; the black open upward-facing triangles show Surface layer Turbulence and Environmental Science Test (SLTEST) field measured data, Reτ = 6.28 × 105; the solid black circles indicate the Melbourne wind tunnel experiment data, Reτ = 1.801 × 104; the solid black squares represent the William B. Morgan Large Cavitation Channel experiment data acquired by the US Navy, Reτ = 6.867 × 104; the solid black leftward-facing triangles show the Princeton Superpipe experiment data, Reτ = 9.819 × 104, the observational and experimental data obtained by Marusic et al. [33]; the open symbols present the QLOA field-measured data from Wang and Zheng [9]; the black open squares show Reτ = 1.63 × 106; the black open circles represent Reτ = 2.19 × 106; the black open leftward-facing triangles give Reτ = 2.95 × 106; the black open rightward-facing triangles provide Reτ = 3.14 × 106; and the black open stars show Reτ = 4.19 × 106.

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 Reτ = 7.7 × 105, and the QLOA field measured data with higher Reynolds number of Wang and Zheng [9].

**Figure 3.** Comparison of Reynolds stress with inner-scaled height: The solid colored symbols are current simulation data, the dashed line is Chauhan [34] for Reτ = 2.33 × 106, the solid line represents Reτ = 3.76 × 106, the solid black upward-facing triangles exhibit the wind tunnel results of Graaff and Eaton [37] at Reτ = 1.35 × 103, the solid black downward-facing triangles give the results of the field observations of Hutchins et al. [36] at Reτ = 7.7 × 105, the open symbols show the field-measured QLOA data of Wang and Zheng [9], the black open squares represent Reτ = 1.63 × 106, the black open circles give Reτ = 2.19 × 106, black open leftward-facing triangles present Reτ = 2.95 × 106, the black open rightward-facing triangles present Reτ = 3.14 × 106, and the black open stars are Reτ = 4.19 × 106.

In addition, *uu*/*u*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*/δ ≈ ∼ 0.036, was compared to the results of QLOA [9], SLTEST [36], and wind tunnel experiments [40,42–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 *uu*/*u*2 τ and the Reynolds number.

**Figure 4.** Comparison of *uu*/*u*2 τ in inner-scaled height simulation results with a similarity formulation from Marusic and Kunkel [41]: The solid colored circles are current simulation data at Reτ = 3.553 × 106, the open circles and squares are for the result of Kunkel and Marusic [38] at Reτ ≈ 3.1 × 106 and 3.8 × 106, respectively, the black pluses exhibit the result of Hutchins et al. [36] at Reτ = 7.7 × 105, the times and open stars represent Metzger et al. [39] at Reτ ≈ 8.3 × 105, the open upper triangles show the result of Hutchins et al. [40], Reτ = 2.8 × 103, 7.3 × 103, and 1.903 × 104, the lines represent the similarity formulations proposed by Marusic and Kunkel [41] calculated at Reτ = 7.7 × 105 (solid black line), Reτ = 3.2 × 106 (black dash line), Reτ = 3.8 × 106 (black dotted line).

**Figure 5.** Variation of normalized streamwise component of normal stress *uu*/*u*2 τ with the Reynolds number: Solid colored circles are the current simulation data; the open symbols indicate the wind tunnel results of Österlund et al. [42], Knobloch and Fernholz [43], Hutchins et al. [40], and Kulandaivelu and Marusic [44]; the solid black five-pointed stars represent the field observations of Hutchins et al. [36]; and the solid black squares are the field observations of Wang and Zheng [9].

#### *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 δ, *U*.

The Reynolds shear stress is shown in Figure 6 against the outer-scaled coordinates *z*/δ. 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*/δ = 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*/δ. 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.

**Figure 6.** Comparison of −*uw*/*u*2 τ in results of outer-scaled height simulation with the experimental data. The solid colored circles indicate current simulation data, the open symbols are experimental data, the squares show Reτ = 1.43 × 103, the upper triangles symbol give Reτ = 2.9 × 103, the stars indicate Reτ = 5.2 × 103, the circle symbols represent Reτ = 1.3 × 104, and the eight-spoked asterisks are Reτ = 3.1 × 104.

**Figure 7.** Comparison of *uu*/*u*2 τ in the results of outer-scaled height simulation with the experimental data. The symbols used are the same as those used in Figure 6.

**Figure 8.** Mean velocity-defect profiles, normalized by friction velocity. The symbols used are the same as those used in Figure 6.

#### **4. Results and Discussion**

#### *4.1. Velocity Spectra Characteristics under Di*ff*erent Surface Roughnesses*

#### 4.1.1. Streamwise Velocity Spectra

Figure 9 shows the wavenumber-normalized spectra of streamwise velocity *Su*/*u*2 τδ 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*δ < 2 (wavelengths larger than 3δ), 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*δ > 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 9.** *Cont*.

**Figure 9.** Power spectra of streamwise velocity fluctuations. (**a**) case 1, (**b**) case 2, (**c**) case 3, (**d**) case 4 and (**e**) case 5. The height indicators for each graph are the same as those shown in (**a**), that is, 7.5 m, 27.5 m, 47.5 m, 67.5 m, 87.5 m, 107.5 m, 127.5 m, and 147.5 m. Further, the boundary layer height is 250 m, 300 m, 350 m, 450 m, and 550 m for cases 1, 2, 3, 4, and 5, respectively.

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 *kSu* 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 λ/δ (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.

**Figure 10.** Pre-multiplied spectra for streamwise velocity fluctuation. (**a**) case 1, (**b**) case 2, (**c**) case 3, (**d**) case 4 and (**e**) case 5. The height indicators for each graph are the same as those shown in Figure 9a, that is, 7.5 m, 27.5 m, 47.5 m, 67.5 m, 87.5 m, 107.5 m, 127.5 m, and 147.5 m.

The *k*δ decreases at the first peak (meaning (*k*δ)0) with increasing height; when the height reaches to a certain extent, (*k*δ)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*δ)0 and the dimensionless height *z*/δ.

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.

**Figure 11.** Variation in the peak wavenumber of VLSMs with height. The solid colored circle is the current simulation data, the open symbols show the field-measured data of Wang and Zheng [9], and the solid black symbols represent the experimental data of Valliviki et al. [46].

#### 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.

**Figure 12.** *Cont*.

**Figure 12.** Power spectra of wall-normal velocity fluctuations. (**a**) case 1, (**b**) case 2, (**c**) case 3, (**d**) €e 4 and (**e**) case 5. The height indicators for each graph are the same as those shown in Figure 9a, that is, 7.5 m, 27.5 m, 47.5 m, 67.5 m, 87.5 m, 107.5 m, 127.5 m, and 147.5 m.

#### *4.2. E*ff*ects 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δ, as shown in Figure 13a,b, respectively. Similarly, band-pass filtering with cutoff lengths of 0.3δ − 3δ and δ − 3δ 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δ<λ< 3δ, Guala et al. [4]), the energy fraction is larger than the field-measured data; but in Figure 13d, the energy fraction of the LSMs (δ<λ< 3δ, 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δ − 3δ, 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.

**Figure 13.** Variations in the total energy (**a**) and energy fraction of VLSMs (**b**) and LSMs (**c**, 0.3δ<λ< 3δ; **d**, δ<λ< 3δ ) with height. Solid colored circles are current simulation data, the open symbols represent the field-measured data of Wang and Zheng [9], and the solid black symbols give the experimental data of Balakumar and Adrian [2].

#### *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*)=(2500*m*, 500*m*, 35*m*)) 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 λ/δ > 3 for the first three cases, indicating the existence of the VLSMs. For cases 4 and 5, the 0.3 < λ/δ < 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 (Reτ = 3.538 × 106), 14.04◦ for case 2 (Reτ = 4.894 × 106), 16.70◦ for case 3 (Reτ = 6.559 × 106), 14.04◦ for case 4 (Reτ = 9.927 × 106), and 19.29◦ for case 5 (Reτ = 1.199 × 107). 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 Reτ = 3.14 × 106, while for Reτ = 1.63 × 106, 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.

**Figure 14.** Two-point correlation contours of streamwise-fluctuating velocity.

In addition, the streamwise, wall normal, and spanwise extent of fluctuating-velocity two-point correlations are shown in Figure 15. The distance, *Lxuu*, 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 *Lxuu*/δ = 0.65 between *z*/δ = 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 Reτ = 3.14 × 106 measured data, almost same as in case 1 (Reτ = 3.553 × 106). 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, *Lxuu*/δ decreases, which agrees with the results of Krogstad and Antonia [56], who found *Lxuu*/δ was about 50% lower on their rough wall. The wall-normal extent of the two-point correlations of the fluctuating velocity, *Lzuu*, 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 *Lzuu* cannot be obtained for *z*/δ < 0.2. The conclusion is similar to the *Lxuu*, as shown in Figure 15b. The spanwise extent *Lyuu* shown in Figure 15c, with increasing roughness, the *Lyuu*/δ 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δ − 0.2δ in the spanwise direction. Comparing Figure 15a,b, the ratio of *Lxuu*/*Lzuu* 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].

**Figure 15.** The extent of two-point correlations of the fluctuating velocity, equal to the 0.5 contour as a function of *z*/δ. (**a**) Streamwise extent, (**b**) wall-normal extent, (**c**) spanwise extent, and (**d**) ratio between streamwise and wall-normal extent. The solid colored symbols are current simulation data, the solid black data are experimental data of Volino et al. [17], the circle symbols represent smooth wall condition, and squares represent roughness wall condition, the field measured data of Wang and Zheng [9], the anise star is Reτ = 3.14 × 106, and five star is Reτ = 1.63 × 106 which show in (**a**).

#### *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–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δ in streamwise direction, and the widths is around 0.1δ − 0.2δ. 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 16.** The normalized streamwise velocity fluctuations in *xy* plane. (**a**) case-1, (**b**) case-2, (**c**) case-3, (**d**) case-4, (**e**) case-5.

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.

**Figure 17.** Streamwise velocity fluctuations normalized by the friction velocity in *xz* plane of case-3. The lines (**A**), (**B**) and (**C**) in the Figure 16c are shown here.

#### *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 18.** Time evolution of the instantaneous fields at a height of 30 m in *xy* plane. The dashed line D1, D2 are added to show the growth of the VLSMs, dashed line E is added to show the extinction of the VLSMs, F1, F2 are added to show the growth first then extinction of the VLSMs and G1, G2 are added to show the extinction first then growth of the VLSMs. P-number represents the different part of VLSMs of Figure 18f. Subfigures (**a**–**j**) represent the time evolution of flow field.

Figures 19 and 20 show the associated flow features in the *xz* plane. The conclusion is consistent with the results show in Figure 18. Combining the Figures 18–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.

**Figure 19.** Evidence of the streamwise scale growth and extinction in *xz* plane. Subfigures (**a**–**j**) represent the time evolution of flow field.

**Figure 20.** Evidence of the streamwise scale growth in *xz* plane. Subfigures (**a**–**j**) represent the time evolution of flow field.

#### *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 λ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].

**Figure 21.** Instantaneous vortical structures visualized with an iso-surface of λ2 in *xy* plane.

**Figure 22.** Instantaneous vortical structures visualized with an iso-surface of λ2 in *xz* plane, the subfigures (**A**), (**B**) and (**C**) are the same with the Figure 18c, represent the cross section of vortices field.

#### **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* (106–107), 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*δ)0 decreases with height in an approximately log-linear manner. While, when height increases to a certain extent, (*k*δ)0 begins to reach a constant value or even to increase. In addition, *kSu* 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.

**Author Contributions:** Conceptualization, H.R. and S.L.; Methodology, H.R. and S.L.; Software, H.R.; Validation, H.R.; Formal Analysis, H.R.; Investigation, H.R.; Resources, H.L.; Data Curation, H.R.; Writing-Original Draft Preparation, HR.; Writing-Review & Editing, S.L.; Visualization, H.R.; Supervision, S.L. and H.L.; Project

Administration, H.L.; Funding Acquisition, S.L. and H.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study is financially suppborted by the National Key Research and Development Program of China under grant No. 2018YFC0705605, NSFC under grant No. 51878230, the China Postdoctoral Science Foundations (2015M580266 and 2016T90295), the Heilongjiang Postdoctoral Fund (LBH-Z15087), and the Fundamental Research Funds for the Central Universities (HIT. NSRIF. 2017064.)

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Surrogate Model with a Deep Neural Network to Evaluate Gas–Liquid Flow in a Horizontal Pipe**

**Yongho Seong 1, Changhyup Park 1,\*, Jinho Choi 2 and Ilsik Jang 3**

Received: 3 January 2020; Accepted: 20 February 2020; Published: 21 February 2020

**Abstract:** This study developed a data-driven surrogate model based on a deep neural network (DNN) to evaluate gas–liquid multiphase flow occurring in horizontal pipes. It estimated the liquid holdup and pressure gradient under a slip condition and different flow patterns, i.e., slug, annular, stratified flow, etc. The inputs of the surrogate modelling were related to the fluid properties and the dynamic data, e.g., superficial velocities at the inlet, while the outputs were the liquid holdup and pressure gradient observed at the outlet. The case study determined the optimal number of hidden neurons by considering the processing time and the validation error. A total of 350 experimental data were used: 279 for supervised training, 31 for validating the training performance, and 40 unknown data, not used in training and validation, were examined to forecast the liquid holdup and pressure gradient. The liquid holdups were estimated within less than 8.08% of the mean absolute percentage error, while the error of the pressure gradient was 23.76%. The R2 values confirmed the reliability of the developed model, showing 0.89 for liquid holdups and 0.98 for pressure gradients. The DNN-based surrogate model can be applicable to estimate liquid holdup and pressure gradients in a more realistic manner with a small amount of computating resources.

**Keywords:** surrogate model; deep neural network; multiphase flow; horizontal pipe; liquid holdup; pressure gradient

#### **1. Introduction**

The accurate evaluation of multiphase flow has been essential not only for optimum facility designs, but also for the estimation of transport features in pre-installed pipes. The oil and natural gas industries require time-consuming experimental analyses, but the flow characteristics, e.g., flow patterns, liquid holdup, superficial velocities, and pressure gradient, are uncertain. These uncertain parameters are important for demonstrating transport phenomena but they are linked nonlinearly. A few analytical models implement closure relationships that depend on experimental data and therefore they have the limitation of the applicable ranges used in the experimental works. The experiments validated the numerical simulations with the closure relationships, as well as the mechanistic modules, e.g., OLGA (Schlumberger, Houston, USA) and Ledaflow (Kongsberg, Norway). However, the detailed designs of these numerical simulations become complicated, e.g., computational fluid dynamics, and require a large amount of computing resources and skillful user interactions [1,2].

Many works have validated their effectiveness using the experimental database, e.g., TUFFPDB (University of Tulsa Fluid Flow Projects Database) and the experimental database of the University of Amsterdam. The empirical interpretations proposed different interpolation rules within the

experimental ranges [3–9]. Choi et al. [6] introduced a drift–flux closure relationship to interpolate the experimental results and estimate the liquid holdup. Their model estimated the liquid holdups of experimental data from 0.042 up to 0.156 of the mean absolute error. Lee et al. [7] examined the applicability of the minimum dissipated energy concept to explain the stratified gas–liquid flow in horizontal pipes and confirmed its usefulness by estimating liquid holdups and pressure gradients. Luo et al. [9] proposed an empirical power-law model not subject to the experimental ranges to predict the liquid holdup in vertical gas wells.

Machine learning methods have been actively implemented to predict the multiphase flow characteristics [10–15]: some studies implemented an artificial neural network [10–12] and Mask et al. [13] implemented dimensionless parameters to estimate the flow patterns. Kanin et al. [14] trained a machine learning algorithm to validate the mechanistic models. Mohammadi et al. [15] used a genetic algorithm to select closure relationships in multiphase flow models. These studies have focused on securing the most accurate ways to explain the nonlinearity more effectively between input and output data. Deep neural networks (DNNs), as a class of machine learning techniques, implement a few hierarchical hidden layers of non-linear transforms and improves the performances of a typical artificial neural network with one hidden layer [16–19]. The neural networks normalize the input data such that they can integrate different-scaled data. Added hidden layers reduce the number of calibrating parameters, thereby it can seek the non-linear relationships among the parameters more effectively, i.e., deep learning.

An objective of this paper was to develop a DNN framework, i.e., a data-driven surrogate (proxy) model, to estimate the liquid holdup and pressure gradient for a gas–liquid multiphase flow in horizontal pipes instead of using numerical simulations. DNN with multiple hidden layers was applied to secure an accurate empirical relationship between the experimental input factors and the estimating values. The experimental data were used to validate the reliability, to train the neural network, and to predict the flow-related parameters.

#### **2. Methodology**

The representative parameters used to explain gas–liquid flow characteristics are the liquid holdup and pressure gradient; the liquid holdup, i.e., the fraction of a part of pipe occupied by liquid, demonstrates the amount of liquid transported through pipes, while the pressure gradient, i.e., pressure drop divided by the given length of the horizontal pipe, is essential for designing the transport facilities. Various empirical correlations have been proposed to explain these parameters, e.g., Equations (1) and (2) for a horizontal pipe [20–22], but these correlations are challenging to use due to the complexity and the uncertainty of multiphase flow.

$$H\_L = \begin{cases} \lambda\_L = \frac{A\_l}{A\_p} = \frac{q\_L}{q\_L + q\_S} = \frac{v\_{SL}}{v\_{SL} + v\_{SG}}, & \text{no-slip condition} \\ \frac{a\_1 \lambda\_L^2}{\left(\lambda\_{FR}^2\right)^{0.5}}, & \text{slip condition} \end{cases} \tag{1}$$

$$\frac{dP}{dL} = \begin{cases} -\frac{f \rho\_v v\_m^2}{2d} - \rho\_n v\_m \frac{dv\_m}{dL}, & \text{homogeneous flow} \\ \qquad - \frac{f \rho\_v v\_m^2}{1 - E\_k}, & \text{non-homogeneous flow} \end{cases} \tag{2}$$

In Equations (1) and (2), the no-slip liquid holdup (λ*L*) is the ratio of the liquid volumetric flow rate to the total volumetric flow rate [22]. *AP* represents the cross-sectional area of pipe and *AL* is the cross-sectional area occupied by liquid in the pipe. *qL* is the flow rate of liquid and *qg* is the flow rate of gas. *vSL* = *qL AP* is the liquid superficial velocity and *vSG* = *qg AP* is the gas superficial velocity. *a*1, *a*2, and *a*3 are the coefficients for liquid holdup correlation, which are different for different flow patterns. *f* is the friction factor related to the pipe diameter (*d*) and roughness. *NFR* = *v*2 *m gd* is the mixture Froude number, where *vm* = *vSL* + *vSG* (mixture velocity). *Ek* is a dimensionless kinetic or acceleration energy term. ρ*n* is the no-slip density. The no-slip condition occurs where the input

phase-superficial velocities are constant in the course of the gas–liquid flow, e.g., a typical stratified flow pattern. A homogeneous two-phase flow is defined as the flow structure in which the two phases travel at the same in-situ velocity while a non-homogeneous flow means the flow structure in which the mixture physical properties vary across the pipe's cross-sectional area [22].

The experimental data were from Gokcal [23,24], which considered a viscous oil–air mixed flow in horizontal pipes, where the viscous oil was Citgo sentry 220 oil with a density range from 833.6 to 884.5 kg/m3 (Table 1). Gokcal [23] investigated various flow patterns while Gokcal [24] investigated slug flow. The experimental data showed a slip condition since the liquid holdups resulting from the multi-phase flow experiment were not equal to the no-slip liquid holdups (λ*L*). The total number of available data was 350: 279 for the supervised training of the DNN, 31 for validating the training performance, and the other 40 data were for testing the predictability. A reason why 310 training data were divided into the training and validation sets was to examine any overfitting or underfitting problems. A data-shuffling process was implemented to split the data. A total of 40 experimental data that were not used in either training or validation were randomly selected for testing the predictability (Table 2).

**Table 1.** Summary of experimental data for the input and the output parameters in the deep neural network (DNN).

1 *ID*: Inner diameter, *T*: Temperature; 2 (number) indicates the number of available experimental data with the specified flow pattern; The subscripts *L*, *G*, *SL*, *SG* stand for liquid, gas, liquid superficial velocity, and gas superficial velocity, respectively.

**Table 2.** Experimental data and their flow patterns used in the DNN model.

Figure 1 explains how to establish empirical correlations using a backpropagation algorithm between the input and output data in multi-hidden-layered neural networks (Equation (3); [18,25]).

$$\begin{array}{ll} a\_i^{(1)} = f\_i \Big( \sum\_{j=1}^n w\_{i,j}^{(1)} x\_j + b\_j^{(1)} \Big) \quad \forall i = 1, \ 2, \ 3, \ \cdots, m, \\\ a\_i^{(2)} = f\_i \Big( \sum\_{j=1}^n w\_{i,j}^{(2)} a\_j^{(1)} + b\_j^{(2)} \Big) \quad \forall i = 1, \ 2, \ 3, \ \cdots, l, \\\ g\_i = f\_i \Big( \sum\_{j=1}^h w\_{i,j}^{(3)} a\_j^{(2)} + b\_j^{(3)} \Big) \quad \forall i = 1, \ 2, \ 3, \ \cdots, l. \end{array} \tag{3}$$

**Figure 1.** DNN structure to establish the empirical relationship between the input and the output layer.

In Equation (3), *ak i* denotes the value of the *i*th neuron at the *k*th hidden layer. *fa* is an activation function. *w* is a weight and *b* is a bias. *x* is the value in the input layer and *y*ˆ is the estimated value at the output layer. Figure 2 depicts the DNN framework consisting of one input, one output, and two hidden layers in this study. The performance of estimating the output neurons depended on the design of hidden layers, e.g., the number of neurons in the hidden layers, the activation function, and the filters.

**Figure 2.** DNN design to estimate the liquid holdup and pressure gradient. One input layer (*X*), one output layer (*Y*ˆ), and two hidden layers (*HDNN1* and *HDNN2*) were implemented.

The input layer (*X*) was composed of the normalized properties, such as the liquid superficial velocity, gas superficial velocity, liquid density, gas density, and temperature (see Equation (4)). On the other hand, the output layer (*Y*ˆ) had the liquid holdup and pressure gradient. Two hidden layers were introduced and the effects of setting the number of neurons were examined as a result of changing the neurons in the hidden layers in terms of the computation cost, i.e., the processing time, and the validation loss, i.e., the mean square error between the estimated values and the true (experimental) data for cost function. The number of nodes in the second hidden layer (*HDNN2*) was assumed to be half of the node number in the first hidden layer (*HDNN1*).

$$\mathbf{x}\_{\text{norm}} = \frac{\mathbf{x} - \mathbf{x}\_{\text{MIN}}}{\mathbf{x}\_{\text{MAX}} - \mathbf{x}\_{\text{MIN}}} \tag{4}$$

In Equation (4), *xnorm* in *X* represents the normalized *x* value, and the subscripts *MIN* and *MAX* are the minimum value and the maximum value of *x*. The Nesterov accelerated gradient was introduced as a momentum optimizer (Equation (5); [26]) and the rectified linear unit (ReLU) was used for the activation function. The training continued until the optimizer converged before 40 epochs. Some methods, e.g., dropout, batch normalization, and weight regularization, was used to try to fix the overfitting matter.

$$\begin{aligned} \partial\_t &= \chi \partial\_{t-1} + \eta \nabla\_w J(w - \chi \partial\_{t-1}) \\ w &= w - \partial\_t \end{aligned} \tag{5}$$

In Equation (5), *w* is the parameter to update, i.e., weight. θ*t* denotes the update vector of the current time, *t*. γ is the momentum term that this work set to 0.9. *J* is a loss function, i.e., cost function, and ∇*wJ* is a gradient vector of the loss function. η represents the learning rate that was assigned 0.001 in this study. He et al.'s scheme initializes the weights applicable to ReLU (Equation (6); [27]). *W* is a set of weights. σ2(*W*) is the variance of weights and *nin* is the number of nodes in the previous layer.

$$\begin{array}{c} \mathcal{W} \sim \mathrm{N} \Big( 0, \,\sigma^2(\mathcal{W}) \Big) \\ \sigma^2 = \frac{2}{n\_{\mathrm{in}}} \end{array} \tag{6}$$

Three kinds of errors were applied: the root mean squared error (RMSE; Equation (7)), the percentage error (Equation (8)), and the mean absolute percentage error (MAPE; Equation (9)):

$$
\varepsilon\_{RMSE} = \sqrt{\frac{1}{n} \sum\_{i=1}^{n} \left(\boldsymbol{Y}\_i - \hat{\boldsymbol{Y}}\_i\right)^2} \tag{7}
$$

$$
\varepsilon\_{PE} = 100 \left| \frac{\mathbf{Y}\_i - \hat{\mathbf{Y}}\_i}{\mathbf{Y}\_i} \right|,
\tag{8}
$$

$$
\varepsilon\_{MAPE} = \frac{100}{n} \sum\_{i=1}^{n} \left| \frac{\mathbf{Y}\_i - \mathbf{Y}\_i}{\mathbf{Y}\_i} \right|. \tag{9}
$$

In Equations (7)–(9), ε*RMSE* denotes RMSE, *n* is the number of output data points, *Yi* represents the actual experimental data, and *Y*ˆ*i* are the values estimated using the DNN. The percentage error (ε*PE*) measures the difference between the predicted and the experimental value, divided by the experimental value, and expresses the error as a percentage. The MAPE (ε*MAPE*) is the arithmetic mean value of the percentage errors and is implemented for both the validation and testing sets. The optimum number of nodes in the hidden layers was determined using the results of the processing time, i.e., training time to construct the DNN model, and the MAPEs.

#### **3. Results and Discussion**

#### *3.1. Design of the Deep Neural Network: The Number of Hidden Neurons*

Both performances of training with validation loss and the predictability with the testing dataset were examined by changing the number of nodes in the hidden layers. Since the node number of the second hidden layer was set as half of that in the first hidden layer, the performances were observed in terms of the number of the first hidden layer from 10 to 1000. The processing time refers to the training time, i.e., time for constructing the DNN using supervised training. It was measured when the optimization scheme, i.e., the Nesterov accelerated gradient, converged until 40 epochs, thereby it can be different with the number of nodes in the hidden layers and the final epoch number. The measuring computer system was as follows: Intel Xeon CPU E5-1620v3 @3.50GHz, 64GB RAM, and NVIDIA GeForce GTX 1060 6GB.

Table 3 summarizes the range of the validation loss to evaluate the DNN training performances. The RMSEs of the validation dataset ranged from 0.0632 to 0.0742 (unitless) for the liquid holdup and from 615.87 to 903.84 Pa/m for the pressure gradient. The predictability of the pressure gradient was less than that of liquid holdup; the reasons for this could be the high nonlinearity and complexity of the pressure gradient. The large range of the pressure gradient, e.g., from 6000 to 8000 Pa/m, caused a high RMSE such that both the RMSE and MAPE should be considered.

**Table 3.** Summary of the processing time and the errors for validating the DNN training performances.

1 The number of neurons in the second hidden layer assumes the half of given number in the first hidden layer. RMSE: Root Mean Squared Error; MAPE: Mean Absolute Percentage Error. *HL*: liquid holdup (unitless); *dP*/*dL*: pressure gradient (Pa/m) in RMSE

Figure 3 depicts the processing time and the MAPE (ε*MAPE* defined in Equation (9)) of liquid holdups estimated using the validation data (31 data not used in training). Figure 4 demonstrates those of the pressure gradient. Both overall trends of processing time increased as the number of nodes in hidden layers increased. The processing time for less than 100 nodes in the first hidden layer was not high, but if applying over 100 nodes, the processing time increased dramatically in the case of predicting the liquid holdup (Figure 3). For the pressure gradient, the trend of processing time was similar to that of the liquid holdup but the errors fluctuated in the range of 30 to 100 for the number of nodes in the first hidden layer. The overall trend of errors, i.e., both for the liquid holdup and pressure gradient, converged after 70 nodes and decreased by up to 30% of the MAPE (Figures 3 and 4). The errors did not necessarily decrease even when the number of neurons in the hidden layer increased. This trend of processing time and error shows that the optimum conditions satisfying the acceptable training time and errors could be contentious because this study assumed several fixed constraints: the neuron number in the second hidden layer was half that in the first hidden layer and the converging criterion was limited until the given epoch. Figures 3 and 4 can be used to recommend 70 nodes in the first hidden layer and 35 nodes in the second hidden layer as the optimum DNN model with the consideration of the processing time and the training accuracy for both the results of validating the liquid holdup and pressure gradient.

**Figure 3.** Processing time and training accuracy of the liquid holdup (unitless) for the validation data (31 data points) in relation to the number of nodes in the first hidden layer (*HDNN1*) in the range of (**a**) 10 to 30, (**b**) 30 to 100, and (**c**) 100 to 1000.

**Figure 4.** Processing time and training accuracy of pressure gradients (Pa/m) for the validation data (31 data points) in relation to the number of nodes in the first hidden layer (*HDNN1*) in the range of (**a**) 10 to 30, (**b**) 30 to 100, and (**c**) 100 to 1000.

#### *3.2. Prediction Accuracy of the Liquid Holdups and Pressure Gradients*

Figure 5 compares the values predicted using the DNN with those of the experimental results; the proposed DNN implemented 70 nodes in the first hidden layer and 35 nodes in the second hidden layer. The MAPEs were 8.08% for the liquid holdup and 23.76% for the pressure gradient. The R2 showed a high correlation, i.e., 0.89 for the liquid holdup and 0.98 for the pressure gradient, which was enough to confirm the robustness of the proposed DNN (Table 4). The predictability of the pressure gradients was reliable given the high R2 value, notwithstanding that they had a wide range from 0 to near 8000 Pa/m in this study. The predictability of the liquid holdup was better or similar to the existing closure relationships [6,8,10,11], e.g., 3.7%–24.5% [6]. The errors of the pressure gradient were not less than the results of the previous works, e.g., around 4%–36% [5,7,8]. Since the available data and the methods in the previous works are different, a direct comparison of the results does not provide a reliable evaluation. Therefore, the errors and the R2 values, i.e., the correlation between the true and the estimated values, must examine the reliability of the developed model.

**Figure 5.** Prediction accuracy of the proposed DNN model: (**a**) liquid holdup and (**b**) pressure gradient. A total of 40 experimental data points, not from within the training and validation data sets, were used to evaluate the prediction performances.

Figure 6 depicts the percentage errors (Equation (8)) with the values estimated from the DNN. The largest percentage error was 67.23%, for the case of annular flow, with a 0.12 (relatively a small value) liquid holdup, while that of the pressure gradient was 167.0%, observed for slug flow, with 48 Pa/m (a small value). Given the definition of MAPE (Equation (9)), the smaller the value of the denominator, the more the error tends to increase. Five outliers, with percentage errors of 50% or more, influenced the MAPE; the number of outliers when predicting liquid holdups was 1 and that for pressure gradients was 4. All singular values were in the zone where the true values were small. Therefore, Figure 6 proves that the developed model could accurately estimate both the liquid holdup and the pressure gradient, except for the outliers.

**Figure 6.** Scatter plot of the percentage errors using the testing dataset: (**a**) liquid holdup and (**b**) pressure gradient.

A notable result is that the DNN-based model can be applicable to reliably estimate the flow characteristics without time-consuming numerical simulations and user interactions, even though the flow patterns are various in horizontal pipes. The factors influencing the liquid holdup or pressure gradient can be different and thereby the separate model composed of different input parameters can improve the prediction performances. This study discusses the characteristics of gas–liquid multiphase flows that are easily identified using only five simple input parameters. Obtaining enough experimental data to improve the DNN workflow is challenging; the time-consuming work and the complexity of the experimental designs limit the ability to obtain big data suitable for training neural networks and for improving the predictability under various experimental conditions. This paper is limited for the case of horizontal flow and a few flow patterns. If the training and testing data with different flow patterns is enough, a more quantitative analysis would be available that is related to the effects of flow patterns on the prediction accuracy. Constructing the standard forms of neural networks, e.g., determination of the optimum activation functions and elements of the neural networks, is still challenging to obtain to improve the robustness of the deep learning. Enhancement of the neural network requires the intensive examinations of various activation functions and details of designing neural networks to yield more reliable predictions.

#### **4. Conclusions**

This study developed a noble machine learning approach based on a deep neural network to evaluate the liquid holdup and pressure gradient of a gas–liquid two-phase flow in a horizontal pipe. The optimal design of two hidden layers in the neural network was obtained from the comparison of the processing time and training accuracy. Experimental results with different flow patterns in a horizontal pipe examined the reliability of the proposed DNN-based surrogate model: the mean absolute percentage errors were 8.08% for the liquid holdup and 23.76% for the pressure gradient, while the R2 values were 0.89 for the liquid holdup and 0.98 for the pressure gradient. This study developed the data-driven surrogate model to accurately estimate both the liquid holdup and pressure gradient without time-consuming numerical simulations.

**Author Contributions:** Conceptualization, C.P. and J.C.; formal analysis, C.P. and I.J.; investigation, Y.S. and I.J.; methodology, Y.S. and I.J.; resources, J.C.; software, Y.S.; supervision, C.P.; validation, C.P.; project administration, C.P.; funding acquisition, C.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2017R1D1A1B04033060), by the Korea Institute of Energy Technology Evaluation and Planning (KETEP) and the Ministry of Trade, Industry & Energy (MOTIE; Nos.<PHONE_NUMBER>2150;<PHONE_NUMBER>2160), and by 2016 Research Grant from Kangwon National University, Korea (No. 520160435).

**Acknowledgments:** The experimental data were from Gokcal's works [23,24] (The University of Tulsa, Tulsa, OK, USA).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Numerical Simulation of Coherent Structures in the Turbulent Boundary Layer under Di**ff**erent Stability Conditions**

#### **Shujin Laima 1,2,3, Hehe Ren 1,2,3,\*, Hui Li 1,2,3 and Jinping Ou 1,2,3**

Received: 7 February 2020; Accepted: 25 February 2020; Published: 1 March 2020

**Abstract:** Coherent structures in the turbulent boundary layer were investigated under different stability conditions. Qualitative analyses of the flow field, spatial correlation coefficient field and pre-multiplied wind velocity spectrum showed that the dominant turbulent eddy structure changed from small-scale motions to large- and very-large-scale motions and then to thermal plumes as the stability changed from strong stable to neutral and then to strong unstable. A quantitative analysis of the size characteristics of the three-dimensional turbulent eddy structure based on the spatial correlation coefficient field showed that under near-neutral stability, the streamwise, wall-normal and spanwise extents remained constant at approximately 0.3δ, 0.1δ and 0.2δ (δ, boundary layer height), respectively, while for other conditions, the extent in each direction varied in a log-linear manner with stability; only the spanwise extent under stable conditions was also independent of stability. The peak wavenumber of the pre-multiplied wind velocity spectrum moves towards small values from stable conditions to neutral condition and then to unstable conditions; thus, for the wind velocity spectrum, another form is needed that takes account the effects of the stability condition.

**Keywords:** coherent structures; turbulent boundary layer; stability; pre-multiplied wind velocity spectrum; spatial correlation coefficient field

#### **1. Introduction**

Over flat and homogeneous terrain, turbulence is a very important factor affecting turbulent motions near the ground. Consequently, the structure of the turbulent boundary layer depends mainly on vertical turbulence transport which has a great influence on mass, momentum, heat, water vapor and energy exchanges [1] with turbulent coherent structures playing an important role [2,3]. Field measurements, laboratory experiments and numerical simulations of the turbulent boundary layer have revealed the existence of a wide variety of coherent structures such as low-speed streaks, hairpin vortices, large-scale motions and very large-scale motions in canonical wall-bounded flows including pipe flows, channel flows and turbulent boundary layer flows [2,4–10]. Actually, coherent structures are characterized by a specific scale and lifetime, the recurrence of characteristic and flow-specific organization and regions of concentrated vorticity. Specifically, Theodorsen [4] proposed the existence of hairpin vortices, and such vortices, together with streamwise packets, have been suggested to be the primary elements forming coherent structures in the neutral turbulent boundary logarithmic layer under experimental conditions [6]. However, most studies have focused on neutral conditions. Under

neutral conditions, turbulent structures observed in the field are similar to those found in experiments. Furthermore, Marusic et al. [3] indicated that turbulence and even the pressure field in the turbulent boundary layer are also similar to those in canonical wall-bounded flows, and they suggested that the reason for this is the presence of coherent structures.

The prevalent view regarding the formation of coherent structures in the turbulent surface layer is the "bottom-up" mechanism, in which active structures emerge mainly at the surface and then grow up into the outer layer [8,11–13]. As is well known, the turbulent boundary layer is usually classified as being stable, neutral or unstable on the basis of the dominant mechanism of turbulence generation and stability [14,15]. Differences in stability will affect energy transport and fluxes between the near-surface layer and upper free layer, through the formation and development of coherent structures [16–20]. The stability depends on its ability to resist three-dimensional turbulent motions. In stable conditions, even slight three-dimensional turbulent motions are difficult or even impossible. By contrast, in unstable conditions, slight three-dimensional turbulent motions tend to become stronger, leading to turbulence and convective activity. Li and Bou-Zeid [18] concluded that coherent structures in the atmosphere must be different from those found in experiments, because in the turbulent boundary layer they can be affected by the surface buoyancy flux under stable (nocturnal) or unstable (daytime) conditions which are not taken into account in most experiments. Specifically, unstable turbulent boundary layers include three kinds of coherent structures: under slightly unstable conditions, large-scale horizontal roll vortices occur, owing to a combination of shear and buoyancy [21–24]; under moderately unstable conditions, individual hairpins arise from the ground, leading to an increase in the inclination angle of the vortex packets [25–28]; and under strongly unstable conditions, plume-shaped structures appear, together with streamwise elongated structures, and the magnitude of the vertical turbulent momentum flux increases [17,25,29,30]. Stable turbulent boundary layers can generally be divided into two kinds: weakly stable and strongly stable [31]. Under weakly stable conditions, because turbulence is nearly continuous and quasi-steady, Monin–Obukhov similarity theory applies [32]. Under strongly stable conditions, small eddy motions dominate in the boundary layer, and the layer is much shallower than those under neutral or unstable conditions; however, our understanding of this is still quite limited [31,33,34] partly because of the intrinsic complexity of fluid dynamical phenomena such as intermittency and Kelvin–Helmholtz instability [31,33,35].

In the present study, to obtain a better understanding of turbulent boundary layer flow under different stability conditions, a large eddy simulation (LES) approach was adopted. The LES technique was proposed for meteorological applications by Deardorff [36,37]; since then, it has become an indispensable tool to study the turbulent boundary layer [38–41]. In LES, large-scale structures are resolved explicitly, while the effect of small-scale structures to large-scale structures is not resolved but is parameterized using a subgrid-scale (SGS) model. Thus, because the unstable turbulent boundary layer is dominated by large-scale structures, LES is an appropriate approach [42–44]. Large eddy simulation has also been used successfully for neutral conditions [17,20,42,44,45]. However, simulation of the stable turbulent boundary layer faces a number of challenges. These arise because the way in which the radiative cooling of the ground leads to stable stratification modifying the flow turbulence is still not completely understood [46], and because the characteristic size of eddies decreases with increasing stability. The application of LES to stable conditions has, therefore, been mostly confined to weakly or moderately stable cases [47,48], and this approach has not been considered appropriate for strongly stable conditions [34,49]. However, Churchfield et al. [50] have suggested an approach that is capable of simulating complex terrain and stable conditions using wall model LES (WMLES), and it is this approach that is adopted in the present study.

Turbulent boundary layers have a complex structure and involve some phenomena that are still not completely understood, especially at high Reynolds numbers [9]. Therefore, the objective of this study was to investigate the characteristics of the turbulent boundary layer under different stability conditions and to quantify the variations in the shapes and sizes of eddy structures with stability. The

LES model and simulation setup are described in Section 2, and the model is validated in Section 3. Simulation results are discussed in Section 4. Finally, conclusions are given in Section 5.

#### **2. Simulation Details**

The governing equations of the LES model are:

$$\frac{\partial \overline{u}\_i}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (\overline{u}\_j \overline{u}\_i) = -\frac{\partial}{\partial \mathbf{x}\_i} \overline{p} - \frac{\partial}{\partial \mathbf{x}\_j} (\pi^D\_{ij}) + f\_i^T \tag{1}$$

$$\frac{\partial \overline{u}\_i}{\partial \mathbf{x}\_i} = 0 \tag{2}$$

For potential temperature transport, the governing equation is:

$$\frac{\partial \overline{\partial}}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (\overline{\mathbf{u}}\_j \overline{\theta}) = -\frac{\partial}{\partial \mathbf{x}\_j} (q\_j) \tag{3}$$

where *xi* represent Cartesian coordinates, *ui* represents the corresponding velocity components, *p* represents the pressure which, divided by the fluid density and τ*D ij* , represents the SGS stress, and *f T i* represents the density-normalized forces. The Coriolis force was not considered, the same as in Fang and Porté-Agel [20]. A bar on a variable indicates the mean grid-filtered quantity. θ is the resolved potential temperature, and *qj* represents the flux of temperature due to the viscous and sub-filter-scale effects. The details of how the SGS momentum flux τ*D ij* was modelled is shown in Ren et al. [51]. The heat flux *qj* at the surface was modelled using a similar method [38]:

$$q\_{\vec{\lambda}} = \begin{bmatrix} 0 \\ 0 \\ q\_3^{\text{tot}} \end{bmatrix} \tag{4}$$

where, in the simulation, the total average temperature flux is specified, and the fluctuating temperature flux *qtot* 3 is created by the wall model.

For numerical discretization, the finite volume method was used; for the viscous and convective terms, a second-order central difference scheme was adopted; and for the unsteady term, a second-order backward scheme was adopted.

Schumann's model [52] was used as the wall model to model the near wall shear stress, and Moeng's model [38] was adopted for the total temperature flux at the surface. The subgrid stress (SGS) model is the one-equation eddy viscosity model [53]; this is because the one-equation eddy viscosity model does not require the assumption in the algebraic eddy viscosity models, that is, local balance between the production and dissipation of SGS energy. The SGS stress can be modelled as:

$$
\pi^{D}\_{ij,\text{sys}} = \frac{2}{3} k\_{\text{sys}} \delta\_{ij} - 2 \nu\_{\text{sys}} \overline{\bf S}\_{ij} \tag{5}
$$

where τ*D ij*,*sgs* is the SGS stress part of the total stress, τ*D ij* . *Sij* is the resolved-scale strain rate tensor, δ*ij* is the Kronecker delta function, *ksgs* is the SGS kinetic energy; and υ*sgs* is the SGS stress eddy viscosity.

Accounting for the historic effect of *ksgs*, because of diffusion, production and dissipation, a transportation equation was derived:

$$\frac{\partial \mathbf{k}\_{s\xi s}}{\partial t} + \frac{\partial \overline{\mathbf{u}}\_{j} \mathbf{k}\_{s\xi s}}{\partial \mathbf{x}\_{j}} = -\tau\_{ij, s\xi s}^{D} \overline{\mathbf{S}}\_{ij} - \mathbf{C}\_{\varepsilon} \frac{\mathbf{k}\_{s\xi s}^{3/2}}{\Delta} + \frac{\partial}{\partial \mathbf{x}\_{j}} (\upsilon\_{s\xi s} \frac{\partial \mathbf{k}\_{s\xi s}}{\partial \mathbf{x}\_{j}}) \tag{6}$$

where Δ represents the filter scale, and *Ck* and *C*ε are the model constants.

In the present study, the size of the computational domain was 5 × 1 × 0.3 km for all cases. When the hexahedral mesh was applied, the mesh grid size was 5 m in each direction from the ground surface to 0.2 km, and above 0.2 km, the mesh grid size was 10 m. For the streamwise and spanwise directions, periodic boundary conditions were applied; this agrees well with Fang and Porté-Agel [20] of LES of turbulent boundary layer flow and with Lee et al. [54] of DNS of channel flow; for the upper boundary, the slip boundary condition was adopted; and for the solid wall, the wall model was used. For the initial conditions, the logarithmic mean wind speed profile was used in all cases (gradient wind height of 250 m and a gradient wind speed of 8 m/s).

The computation timestep was 0.5 s which ensures the stability and viscous stability by a given Courant number. First, to obtain a fully developed flow field and quasi-steady state, a spin-up simulation was performed. For the statistical analysis, to ensure statistical convergence, over 40,000 timesteps were simulated. For the perturbation component, the initial peak was set at the height of 0.015% of the computational domain, and for the horizontal direction, the maximum perturbation was 0.25 m/s; at the beginning of the simulation, for the streamwise and spanwise direction, the number of periods of the perturbations were 4.0 and 20.0, respectively.

The aim of this numerical simulation was to calculate the forms of coherent structures under specifically assumed stability conditions and to quantify the variations in the shapes and sizes of eddy structures with stability. The simulation was conducted within a pre-defined stability range. As is well known, stability refers to the air's tendency to either resist vertical movement (stability) or to rise and create storms (instability) which can be expressed in terms of a dimensionless stability parameter, −*z*/*L*, where *z* represents the height above the surface ground, and *L* represents the Obukhov length, given by *L* = −*u*3 ∗*Tv*/(*kgQv*0), where *k* represents the von Kármán constant, *u*∗ represents the friction velocity, *Qv*0 is the kinematic virtual temperature flux at the surface, *Tv* is a reference virtual temperature, and *g* is the gravitational acceleration. *L* expresses the relative roles of shear and buoyancy in the production/consumption of turbulence kinetic energy with positive and negative values indicating stable and unstable conditions, respectively, and it approaches zero in the limit of neutral stratification. For −*z*/*L*, positive and negative values indicate unstable and stable conditions. In this study, 17 different turbulent boundary layer cases were analysed, from case 1 to 17 which range from stable to unstable, with −*z*/*L* = −98.27 (case 1), −16.19 (case 2), −2.29 (case 3), −1.15 (case 4), −0.84 (case 5), −0.63 (case 6), −0.35 (case 7), −0.05 (case 8), 0 (case 9), 0.54 (case 10), 1.01 (case 11), 1.79 (case 12), 3.16 (case 13), 4.48 (case 14), 8.13 (case 15), 11.31 (case 16), and 31.46 (case 17). It should be noted that to represent the average stability of the turbulent boundary layer, the *z* used for calculating −*z*/*L* here is the boundary layer height δ; this is different from Figures 1 and 2 in which *z* represents the level height above the ground.

**Figure 1.** Ratio of horizontal and vertical heat fluxes *u* θ /*w* θ versus −*z*/*L*. In (**a**), for unstable conditions, the blue circles, red crosses, green plus signs, cyan stars, magenta squares, yellow asterisks, black diamonds and blue rightward-pointing triangles represent cases 17, 16, 15, 14, 13, 12, 11 and 10, respectively, from very to slightly unstable. The curve corresponds to Equation (7). In (**b**), for stable conditions, the blue circles, red crosses, green plus signs, cyan stars, magenta squares, yellow asterisks, black diamonds and blue rightward-pointing triangles represents cases 8, 7, 6, 5, 4, 3, 2, and 1, respectively, going from slightly to very stable.

**Figure 2.** Normalized standard deviation of *w*(σ*w*/*u*∗) versus −*z*/*L*: (**a**) unstable conditions; (**b**) neutral stability conditions; and (**c**) stable conditions. The symbols are the same as in Figure 1, except that the curve in (**a**) now corresponds to Equation (8).

#### **3. Validation**

The basic statistical properties of the velocity under neutral stability conditions were compared with those from previous studies, and good agreement was found when dimensioned both by inner and outer variables; the details are shown in previous studies by Ren et al. [51] and Ren et al. [55].

The ratio of the heat fluxes in the horizontal and vertical directions under unstable conditions is shown in Figure 2. Wyngaard and Cote [56] proposed the following formula for this ratio:

$$\frac{\overline{u'\theta'}}{\overline{w'\theta'}} = a\varphi\_{\mathfrak{m}}\varphi\_{\mathfrak{h}}\tag{7}$$

where ϕ*m* = (*kz*/*u*∗)(∂*u*/∂*z*) represents the non-dimensional wind shear function, *u*∗ is the surface friction velocity, ϕ*h* = (*kz*/*u*∗)(∂θ/∂*z*) is the non-dimensional potential temperature gradient function, and *a* is a constant. In the present study, ϕ*m* = (1 − 28*z*/*L*) −1/4, ϕ*h* = (1 − 20*z*/*L*) −1/2 [57,58], and *a* = 5 [57–59]. In Figure 1, the ratio of the horizontal and vertical heat fluxes decreases with increasing −*z*/*L* away from the near-surface region, and *u* θ /*w* θ tends towards 0 which indicates that at higher elevations the shear stress is very small compared with the thermal plume force. The simulation results show good agreement with theory under strongly unstable conditions, but under moderately or slightly unstable conditions, they deviate from the theoretical values. This deviation occurs because under these conditions, the shear stress driving the horizontal motions near the surface is greater than under strongly unstable conditions, while the vertical transport is weaker than under strongly unstable conditions. Therefore, Equation (7) is appropriate for strongly unstable conditions only, and perhaps, after the forms of *a*, ϕ*m* and ϕ*h* are changed in the present study, it will also prove to be appropriate for moderately or slightly unstable conditions. Another phenomenon revealed in Figure 1 is a shift in the peak of *w* θ /*w* θ towards larger −*z*/*L*; this is because as conditions become more unstable, thermal plumes reach higher. Similarly, the peak of *u* θ /*w* θ shifts to smaller −*z*/*L* as the stability increases under stable conditions. This suggests that as conditions become more unstable, −*z*/*L* becomes larger which is related to the upper location of thermal plumes.

Figure 2 shows the normalized standard deviation of the vertical velocity component, σ*w*/*u*∗, which is important for describing stability. The curves in Figure 2 represent the formula proposed by Panofsky et al. [59] in the case of the unstable condition:

$$
\sigma\_w/u\_\* = 1.3(1 - 3z/L)^{1/3} \tag{8}
$$

where *u*∗ is the friction velocity, *L* represents the Obukhov length, and *z* is the height above the ground. It can be seen that σ*w*/*u*∗ remains constant under slightly or moderately unstable conditions but increases with increasing height under strongly unstable conditions. The simulation results show good agreement with Equation (8) which means that σ*w*/*u*∗ obeys Monin–Obukhov similarity theory under unstable conditions, fitting well with the 1/3 power law. Under neutral conditions, σ*w*/*u*∗ is constant, independent of height and equal to 1.32 (shown in Figure 2b); indeed, a large number of observations show that under neutral conditions, σ*w*/*u*∗ should be constant [58,60–64]. In particular, Bowne and Ball [60] and Counihan [61] both concluded that this constant value is 1.3. Roth [64] summarized the values of different researchers, which ranged from 1.1 to 1.7. When the turbulent boundary layer is under stable conditions, the turbulent motions are weakly connected to the ground, and the height is no longer a similarity parameter of this motion, except in the region near the ground. As can be seen from Figure 2c, the normalized standard deviation of the vertical velocity component behaves as expected for stable conditions. Also, similar to the ratio of horizontal and vertical heat fluxes, the peak of σ*w*/*u*∗ shifts to smaller −*z*/*L* as stability increases, except for neutral conditions; this is related to the upper location of thermal plumes.

#### **4. Results**

The above verification analysis has illustrated the feasibility of the WMLES approach and the accuracy of simulation of turbulent boundary layer flow under different stability conditions. In the following, we investigate the characteristics of turbulent boundary layers under different stability conditions and quantify the variations in the shape and size of eddy structures with stability.

#### *4.1. Flow Field*

We now turn to an analysis of the flow field with the aim of providing a more intuitive view of its basic characteristics and those of the turbulent eddy scale. We consider a horizontal plane (whole domain) at the height of 0.2δ. Figure 3 shows the flow field for all 17 cases. On the whole, under stable conditions, the dominant turbulent eddy structure was small which is consistent with many previous results [31,33,34]. This is because, under stable conditions, the ground temperature is lower than that of the near-surface flow, and the vertical development of the turbulent eddy structure is limited which suppresses the development of three-dimensional turbulent eddy structure thus small-scale motion is prevalent. The dominant turbulent eddy structure scale increases as the stability decreases, and for the weakly stable condition in case 8, very-large-scale motion structures [8,10,65] even occur, just as under neutral conditions. It is obvious under neutral conditions that both large-scale and very-large-scale motion structures are visualized, based on the streamwise-elongated volumes associated with negative streamwise velocity, which is consistent with the results of Dennis and Nickels [66,67], Lee and Sung [65], Baltzer et al. [68] and Lee et al. [54]. Under slightly and moderately unstable conditions, obvious large-scale horizontal roll vortices form which are usually attributed to the combined effects of shear and buoyancy [21,22,24]. These can be seen as corresponding to the extended turbulent eddy structure that appears under neutral stability conditions. Under strongly unstable conditions, the effects of buoyancy dominate those of shear, so there is an obvious distinction between low-velocity and high-velocity regions which extends over the whole spanwise direction. This is consistent with the experimental results of Hommema and Adrian [25].

In Figure 3, going from the strongly stable condition to the neutral stable condition to the strongly unstable condition, some interesting phenomena can be seen. First, the form and shape of the dominant turbulent eddy structure changed from small-scale motions to large-scale and very-large-scale motions and then to regional motions as the turbulent boundary layer became unstable. Second, the dominant size of the turbulent eddy structure varied as the stability changed; a quantitative analysis of this is presented in Section 4.2 based on the three-dimensional spatial correlation coefficient field.

#### *4.2. Spatial Correlation Coe*ffi*cient Field*

To quantify the size characteristics of the three-dimensional turbulent eddy structure under different stability conditions, a flow field correlation analysis was performed with reference location (*x*, *y*, *z*) = (2500 m, 500 m, 35 m). Figure 4 presents contours of the spatial correlation coefficient for all 17 cases. The plots in the left column are for the *xz* plane, and those in the right column are for the *yz* plane; thus, it is possible to visualize the three-dimensional turbulent eddy structure. The conclusion is similar to that from the flow field analysis, namely, that with increasing stability, the dominant size of the turbulent eddy decreases. This clearly shows that under weakly stable, neutral and slightly and moderately unstable conditions, very-large-scale motion structures appeared, while strongly unstable conditions led to thermal plumes. Meanwhile, an interesting phenomenon was that the spanwise size of the dominant vortices remained almost constant, except under strongly unstable conditions which may be due to the fact that the effect of the shear force was mainly in the streamwise direction and that of the buoyancy force was mainly in the wall-normal direction; however, under strongly unstable conditions, the temperature of the ground was higher than that of the adjacent near-surface air, and so the buoyancy force had a positive effect on the upward motion of air near the ground in the whole spanwise direction. The spanwise and wall-normal velocity correlations were different from the streamwise velocity correlation. For the spanwise velocity correlation under stable conditions, the spanwise size of the correlation structure was smaller than the streamwise velocity correlation; moreover, the degree of tilt was smaller than for the streamwise velocity correlation. Under unstable conditions, even light unstable conditions, the plume-shaped structure occurred, meaning that this structure form occurred earlier than the streamwise velocity correlation. For the wall-normal velocity correlation, from stable to unstable conditions, the spanwise and wall-normal sizes of the correlation structure became wide and deep which is very different to the other two direction velocity correlations. These results are not shown here.

**Figure 3.** Streamwise velocity in the *xy* plane (whole domain) at a height of 0.2δ for each case. (**a**) to (**q**) represent cases 1 to 17, respectively, going from very stable, to neutral and, finally, very unstable conditions.

(**e1**) (**e2**)

**Figure 4.** *Cont.*

*\*

**Figure 4.** *Cont.*

**Figure 4.** *Cont.*

**Figure 4.** Contours of the two-point correlation of streamwise velocity fluctuations. Subfigures (**a1**) to (**q2**) means: each row represents case 1 to 17, from very stable, to neutral, and then very unstable conditions, with the left and right columns showing the contours in the *xz* and *yz* planes, respectively.

Figure 5 shows the size of the turbulent eddy structure under different stability conditions at a height of 0.2δ. The distances *Lxuu*, *Lyuu* and *Lzuu* are defined in Christensen and Wu [69] and Volino et al. [70] as twice the distances from the self-correlation peak to the most distant locations on the 0.5 contour of the two-point correlation of the fluctuating velocity in the downstream, spanwise and wall-normal directions, respectively. We can see that under strongly and moderately stable conditions, the streamwise extent *Lxuu* increased in a log-linear manner as −*z*/*L* increased, and that it behaved similarly under strongly unstable conditions. However, under near-neutral stability conditions, *Lxuu* was independent of −*z*/*L* and remained constant at approximately 0.3δ. Similar behaviour was found for the wall-normal extent *Lzuu* which remained constant at approximately 0.1δ. From strongly stable through neutral to slightly unstable conditions, the spanwise extent *Lyuu* remained constant at approximately 0.2δ, while under strongly unstable conditions, *Lyuu* also increased in a log-linear manner with increasing −*z*/*L*. This corresponds to the behaviour of the spatial correlation coefficient field.

**Figure 5.** Extents of the 0.5 contours of the two-point correlations of the fluctuating velocity as functions of −*z*/*L* at a height 0.2δ: (**a**) streamwise extent; (**b**) spanwise extent; and (**c**) wall-normal extent. The left and right columns are for stable and unstable conditions, respectively. For the neutral stability condition, *Lxuu*= 0.32δ, *Lyuu*= 0.20δ and *Lzuu*= 0.12δ.

#### *4.3. Pre-Multiplied Spectral Analysis*

Pre-multiplied spectral analysis is a popular approach for determining the distribution characteristics of vortices of different scales in canonical wall-bounded flows [8,10,71,72]. Figure 6 shows the multipoint averaged [73] pre-multiplied wind velocity spectra for all 17 cases. As we know, the most energetic point in the pre-multiplied wind velocity spectrum represents the dominant turbulent eddy structure. Therefore, under stable conditions, the peak wavenumber of the wind velocity spectrum is located in the high-wavenumber region which indicates that the dominant vortices are small. With decreasing stability, the peak wavenumber moved to smaller values which means that the size of the dominant turbulent eddy structure increased and very-large-scale motion structures occurred under the slightly stable conditions of case 8. Here, it should be noted that very-large-scale motion structures are those with turbulent eddy size larger than 3δ. Under near-neutral stability conditions (slightly stable, neutral and slightly unstable conditions), the stationary stage of the pre-multiplied wind velocity spectrum has a large width, because the stationary stage corresponds to the very-large-scale motion structure [6,8,10,12]. Under strongly unstable conditions, the peak wavenumber is located in the small wavenumber region which means that the dominant turbulent eddy structure is very large, and the shape of the wind velocity spectrum is similar to that of a typhoon wind velocity spectrum near the eye-wall region. This may be because the eye-wall region of a typhoon is also characterized by strongly convective conditions, accompanied by a strong updraft. For spanwise and wall-normal pre-multiplied wind velocity spectra, the maximum energy

points are consistent with the streamwise pre-multiplied wind velocity spectrum, and as the instability increases, it moves to a small wavenumber. The difference is that the steady stage is more obvious in the streamwise pre-multiplied wind velocity spectrum than the other two directions, and that the spanwise and wall-normal pre-multiplied wind velocity spectra have obvious cross-over phenomenon at different heights, which does not occur for streamwise pre-multiplied wind velocity spectrum. These results are not shown here.

**Figure 6.** *Cont.*

**Figure 6.** Pre-multiplied wind velocity spectra of the streamwise velocity fluctuation for each case at different heights: (**a**) to (**q**) are for cases 1 to 17, respectively. The dashed and solid black lines represent 0 and –2/3 slope, respectively.

It can be seen from Figure 6 that as the stability conditions changed from stable to neutral to unstable, the peak wavenumber of the pre-multiplied wind velocity spectrum moved towards smaller values which is consistent with the above analyses of the flow field and the spatial correlation coefficient field, and the stationary stage of the wind velocity spectrum existed only under near-neutral conditions. Under strongly stable or unstable conditions, the small-scale and large-scale turbulent eddy structures were the dominant structures, respectively, and only one dominant peak existed in the pre-multiplied wind velocity spectrum.

The canonical wind velocity spectral styles are based on the neutral stability condition, but in this paper, we found that the different stability conditions gave different wind velocity spectral forms. Therefore, in future studies, we should try to find a stability parameter that represents the differences

in the wind velocity spectra and allows us to obtain a universal wind velocity spectrum that is suitable for all stability conditions.

#### **5. Conclusions**

We investigated the coherent structures in the turbulent boundary layer under different stability conditions and performed qualitative and quantitative analyses of the shape and size of eddy structures.

First, the WMLES approach and the basic variables of the simulation results were validated by comparing these results with theoretical predictions, observed data, experimental data and other simulation results.

Then, based on a flow field analysis, it was found that the form and shape of the dominant turbulent eddy structure changed from small-scale motions to large-scale and very-large-scale motions and then to thermal plumes as the stability condition changed from strong stable conditions to neutral stability conditions and then, to strongly unstable conditions. A spatial correlation analysis was performed, and the results led to similar conclusions as those obtained from the flow field analysis. An interesting phenomenon was also found, namely, that the spanwise size of the dominant vortices remained almost constant, except under strongly unstable conditions which may be because the main effect of shear is in the streamwise direction while that of buoyancy is in the wall-normal direction. From an analysis of the pre-multiplied wind velocity spectrum, it was found that as the stability condition goes from stable to neutral to unstable, the peak wavenumber of the pre-multiplied wind velocity spectrum moves towards smaller values, and the stationary stage of wind velocity spectrum only exists under near-neutral stability, whereas under strongly stable or unstable conditions, the smalland large-scale turbulent eddy structures have an overwhelming advantage, so only one dominant peak exists in the pre-multiplied wind velocity spectrum. In future work, we need to consider a new style for the wind velocity spectrum that takes account of the effects of stability.

Finally, a quantitative analysis of the size characteristics of the three-dimensional turbulent eddy structure under different stability conditions based on the spatial correlation coefficient field showed that for near-neutral stability, the streamwise, wall-normal and spanwise extents remain constant at approximately 0.3δ, 0.1δ and 0.2δ, respectively. For other conditions, it was found that the streamwise, wall-normal and spanwise extents increase in a log-linear manner with increasing stability, except for the spanwise extent under stable conditions which is independent of stability.

**Author Contributions:** Conceptualization, H.R., S.L., H.L. and J.O.; Methodology, H.R. and S.L.; Software, H.R.; Validation, H.R.; Formal Analysis, H.R.; Investigation, H.R. and S.L.; Resources, S.L.; Data Curation, H.R.; Writing–Original Draft Preparation, H.R.; Writing–Review and Editing, S.L.; Visualization, H.R.; Supervision, H.L. and J.O.; Project Administration, H.L.; Funding Acquisition, H.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was financially supported by the National Key Research and Development Program of China under grant no. 2018YFC0705605, NSFC under grant no. 51878230, the China Postdoctoral Science Foundations (2015M580266 and 2016T90295), and the Heilongjiang Postdoctoral Fund (LBH-Z15087).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Numerical Study of the E**ff**ects of the Jet Fan Speed, Heat Release Rate and Aspect Ratio on Smoke Movement in Tunnel Fires**

#### **Ha Thien Khieu, Young Man Lee, Ji Tae Kim and Hong Sun Ryou \***

School of Mechanical Engineering, Chung-Ang University, 84, Heukseok-ro, Dongjak-gu, Seoul 06974, Korea<EMAIL_ADDRESS>(H.T.K<EMAIL_ADDRESS>(Y.M.L<EMAIL_ADDRESS>(J.T.K.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +82-2-820-5280

Received: 4 February 2020; Accepted: 3 March 2020; Published: 5 March 2020

**Abstract:** In this study, the effects of the jet fan speed, heat release rate and aspect ratio on smoke movement in tunnel fires have been investigated. The jet fan speed was changed from 6.25 (25%) to 12.5 m/s (50%), 18.75 m/s (75%), and 25 m/s (100%). The heat release rate was set up from 3.9 to 6 MW and 16 MW, the aspect ratio was changed from 0.6 to 1 and 1.5, respectively. The lower the jet fan speed is, the longer the smoke back-layering length is. With a higher velocity, the smoke tends to move out of the tunnel quickly; however, smoke stratification also occurs, and this reduces visibility. This could make it difficult for people to evacuate. With a higher heat release rate, the smoke tends to move far away from the fires quickly when compared with other cases. Additionally, the higher the heat release is, the longer the smoke back-layering is. Finally, with a higher aspect ratio, the smoke back layering length in the tunnel is also longer. The smoke layer thickness is also larger than in other cases. The correlation of velocity, heat release rate and aspect ratio has been investigated to avoid the smoke back layer length in tunnel fires.

**Keywords:** tunnel fires; jet fan speed; heat release rate; aspect ratio; smoke movement; visibility; smoke layer thickness; smoke stratification

#### **1. Introduction**

Recently, there have been many tunnels under construction worldwide. A tunnel is an underground passageway commonly used in transportation, and studying fires in tunnels is very critical for human safety. When a fire occurs, smoke and toxic gases will be accumulated in the tunnel due to the nature of a tunnel structure. Two parts of typical smoke are CO2 and CO, which are poisonous to humans. As smoke spreads across the tunnel, the smoke prohibits the evacuation of occupants and prevents firefighters from extinguishing the fire.

Many fire accidents have occurred in tunnels. These typically cause significant damage in both human and economic terms, causing many people to be killed and injured and many facilities destroyed; for these reasons, investigating tunnel fires is very important.

In order to prevent the smoke from spreading through an entire tunnel, mechanical ventilation systems are used to transport smoke from the tunnel to the outside. A jet fan is one of the main methods for ventilation in a tunnel; a jet fan can quickly move smoke from vehicles and fires outside of the tunnel. Correspondingly, investigating the jet fan speed on smoke movement in tunnels is necessary for decreasing human casualties.

 In that paper, the use of impulse ventilation for smoke control in underground car parks was performed by using Computational Fluid Dynamics (CFD) simulations. M.G.L.C et al. (2010) [2] performed modeling and simulation of a jet fan for controlling airflow in large enclosures. A

relatively simple CFD model of the free jet fan was proposed and validated. Moreover, the different series modes of jet fans in a longitudinal tunnel ventilation system were investigated by Guihong Pei and Jie Pan (2014) [3]. When most operational fans are located near the tunnel entrance, it is more conducive to increase the velocity and pressure inside the tunnel. Beyer M et al. (2016) [4] has considered the evaluation of the jet fan performance in tunnels. In a tunnel, due to the lack of space, the jet fan is usually installed close to the wall. However, this increases the wall shear stress around the jet, and thus, the effectiveness is significantly reduced.

An experimental study of the effect of the aspect ratio on the critical velocity in longitudinal ventilation tunnel fires was studied by Ryou et al. (2003) [5]. In that study, the aspect ratios of the tunnel were changed from 0.5, 0.667, 1.0, 1.5 and 2.0. The results showed that the aspect ratio affects the growth of smoke fires in the tunnel. The critical velocity increases with the aspect ratio; it was confirmed that high-temperature smoke moved along the ceiling and low-temperature air moved along the lower part of the tunnel.

Moreover, Haukur Ingason et al. (2008) [6] investigated the effect of air velocity on the heat release rate (HRR) and fire characteristics during a fire in the tunnel. The results showed that when increasing the ventilation rate, causing the maximum heat, the release rate increases. The effect of the cross-section on the maximum heat release rates in the tunnel fires was studied by Li et al. (2016) [7]. In that study, they made a model scale tunnel fire with heptane pool fires. For heptane pool fires, decreasing the tunnel height results in an increase of the heat release rate; the heat release rate is about three times that in a free burn test. However, in this research, there is a lack of information about the effect of the aspect ratio on smoke movement and stratification. On the other hand an experimental study on temperature evolution inside the compartment with fire growth and flame ejection through an opening under external wind was considered by Hu, Longhua, et al [8]. They also considered the characterization of buoyant flow stratification and flame extension lengths in the tunnel fires [9,10]. Wan, Huaxian, et al. [11] investigated a numerical study on the smoke back-layering length and inlet air velocity of fires in an inclined tunnel under natural ventilation with a vertical shaft. The tunnel inlet air velocity increases when increasing the tunnel slope and HRR. Besides, Ji, Jie, et al. (2019) [12] also performed a numerical investigation on the influence of the length–width ratio of the fire source on the smoke movement and temperature distribution in tunnel fires. Furthermore, many researches relating to fires or fans were carried out in literature [13–37].

Although there are many kinds of research regarding the jet fan operation and efficiency, there is no research relating to the effect of jet fan speed on smoke movement or stratification in tunnel fires. The purpose of this study is to analyze the effect of the jet fan speed, heat release rate and aspect ratio on smoke movement in a tunnel fire. For this purpose, a numerical study of the effect of the jet fan speed, heat release rate and aspect ratio in the tunnel fires has been carried out by using CFD simulation, as shown in Figure 1. The smoke movement, temperature, velocity and visibility are also investigated in this paper to clarify the impact of the jet fan location and velocity in the tunnel.

**Figure 1.** Numerical analysis objective.

#### **2. Theoretical Analysis**

#### *2.1. Jet Fan Speed*

The jet fan centerline velocity can be calculated as follows:

Awbi (2003) [38] and Baturin (1972) [39] provide the equation regarding the centerline velocity and distance for the experiment. The equation includes the centerline velocity, *um*, the supply velocity, *u0*, distance, *X*, supply diameter of the fan, *d0*, and constant, *a* (from 0.076 to 0.080 for cylindrical tubes).

With a square duct, Kümmel et al. (2007) [40] provided another equation of jet fan experimental data. It includes the length of the potential core, *x0*, the side length of the square duct, *h*, and constant, *m* (from 0.12 to 0.2).

From these equations, the jet fan centerline velocity can be calculated, as shown in the Figure 2, with the supply velocity at 25 m/s, and the side length of the square duct at 0.91 m.

**Figure 2.** Jet fan centerline velocity.

#### *2.2. Heat Release Rate*

In this study, the heat release rate was obtained from our experimental study; one car was assumed to be fired in the tunnel. The maximum heat release rate was 3896 kW at 1551 seconds as shown below Figure 3. Besides, the heat release rate was around 6 MW, and 16 MW for two vehicles [41] and truck fires, respectively.

**Figure 3.** The heat release rate from an experimental study (one vehicle) by Park et al. (2019) [41].

#### *2.3. Smoke Layer Height and Thickness*