Wakes in stratified atmospheric boundary layer flows: an LES investigation

TL;DR

This study uses large eddy simulations and proper orthogonal decomposition to analyze wakes in a stratified atmospheric boundary layer, revealing how thermal stratification influences wake dynamics and turbulence.

Contribution

It introduces a LES-based investigation of bluff body wakes in stratified flows, highlighting the impact of thermal stratification on turbulent wake structures.

Findings

Higher energy POD modes dominate in more turbulent wakes.

Thermal stratification slightly affects turbulent statistics.

Stratification significantly influences the number of dominant POD modes.

Abstract

Large eddy simulations and three-dimensional proper orthogonal decomposition were used to study the interaction between a large stationary and moving bluff body and a high Reynolds number stably-stratified turbulent boundary layer. An immersed boundary method is utilized to take into account the effect of the bluff body on the boundary layer flow and turbulent inflow conditions upstream of the bluff body are imposed by employing a concurrent precursor simulation method both because a pseudo-spectral method is utilized in the horizontal directions. The dominant POD mode type was observed to occur at higher energy levels in the more turbulent wakes. Varying the thermal stratification showed only a slight effect on the turbulent statistics, but a significant effect on the number of dominant POD modes at the highest energy levels.

Figures

Tables

Table 1. Table 1: POD mode number ( k 𝑘 k ) where λ k λ 1 < 0.01 superscript 𝜆 𝑘 superscript 𝜆 1 0.01 \frac{\lambda^{k}}{\lambda^{1}}<0.01 .

	$U_g=8$ m/s	$U_g=18$ m/s
$\Delta \theta =0.2$ K	$504$	$953$
$\Delta \theta =0.4$ K	$673$	$1807$

Table 2. Table 3: Box length Strouhal number for the dominate helical modes.

$S{t}_{l_{box}}$	$U_g=8$ m/s	$U_g=18$ m/s
$\Delta \theta =0.2$ K	$0.26$	$0.31$
$\Delta \theta =0.4$ K	$0.23$	$0.25$

Equations

$$\frac{\partial {\tilde{u}}_i}{\partial t}+{\tilde{u}}_j\frac{\partial {\tilde{u}}_i}{\partial x_j}=-\frac{\partial {\tilde{p}}^{\ast }}{\partial x_i}-\frac{\partial {\tilde{\tau }}_{ij}}{\partial x_j}+{\delta }_{i3}g\frac{\tilde{\theta }-⟨\tilde{\theta }⟩}{{\theta }_0}+F_i$$

$$\frac{\partial \tilde{\theta }}{\partial t}+{\tilde{u}}_j\frac{\partial \tilde{\theta }}{\partial x_j}=-\frac{\partial {\pi }_j}{\partial x_i},\frac{\partial {\tilde{u}}_i}{\partial x_i}=0$$

$$\frac{u_i^{n+1}-u_i^n}{\Delta t}=RH{S}_i+f_i$$

$$f_i=\{\begin{array}{ll}\frac{(u_i{)}_{box}-u_i^n}{\Delta t}-RH{S}_i& ;boxpoints\\ 0& ;fluidpoints\end{array}$$

$${\tau }_{i3}{|}_{z=0}=-u_{\ast }^2\frac{\widetilde{u_i}}{V_f}=-{\left(\frac{\kappa V_f}{\ln \left(\frac{z}{z_0}\right)-{\Psi }_M}\right)}^2\frac{\widetilde{u_i}}{V_f};$$

$$i=1,2$$

$$⟨w^{\prime }{\theta }^{\prime }{⟩}_{z=0}=\frac{u_{\ast }\kappa \left({\theta }_s-\tilde{\theta }\right)}{\ln \left(\frac{z}{z_{0s}}\right)-{\Psi }_H}$$

$$\displaystyle\phi^{blend}_{i}(x,y,z,t)=\$$

$$w(x)=\{\begin{array}{ll}\frac{1}{2}\left[1-cos\left(\pi \frac{x-L_s}{L_{pl}-L_s}\right)\right]& ;L_s\le x\le L_{pl}\\ 1& ;L_{pl}<x\le L_x\end{array}$$

$$u_i(x,y,z,t)={\overline{u}}_i(x,y,z)+\sum _{k=1}^N{a}^k(t){\psi }_i^k(x,y,z),$$

$$u_i(x,y,z,t_n)=u_i(x,y,z,n\tau )$$

$$u_i^{\prime }(x,y,z,t_n)=u_i(x,y,z,t_n)-⟨u_i(x,y,z,t_n){⟩}_n$$

$$C_{nm}=\frac{1}{N}⟨u_i^{\prime }(x,y,z,t_n)u_i^{\prime }(x,y,z,t_m){⟩}_{xyz}$$

$$a^k(t_n)=\sqrt{{\lambda }^{k}N}{e}^k(t_n)$$

$${\psi }_i^k(x,y,z)=\sum _{n=1}^N\frac{a^k(t_n)}{{\lambda }^{k}N}{u}_i^{\prime }(x,y,z,t_n)$$

$$ϵ=⟨⟨u_i^{\prime }(x,y,z,t_n)u_i^{\prime }(x,y,z,t_n){⟩}_{xyz}{⟩}_n=\sum _{k=1}^N{\lambda }^k$$

$$\frac{{\sum }_{k=1}^n{\lambda }^k}{{\sum }_{k=1}^{\infty }{\lambda }^k}>0.99$$

Full text

1 Introduction

Dynamically significant flow structures that are developing around bluff bodies interacting with laminar or turbulent boundary layers can be captured by the use of large eddy simulation (LES), which offers a great potential over other lower order models, such as Reynolds Averaged Navier-Stokes Equations, and a reduced computational cost versus higher order models, such as direct numerical simulations. LES is based on the assumption that the large eddies in the flow are anisotropic and depend on the mean flow and the configuration geometry, while smaller eddies are isotropic and homogeneous and have a universal character. In LES, the instantaneous fluctuating flow field is split into resolved and subgrid scale (SGS) terms using a spatial filter of characteristic size which must be larger (or at least equal) to the numerical grid spacing. The equations are spatially filtered and flow structures that are larger than the filter size are explicitly predicted by the numerical solution to the Navier-Stokes equations, while the influence of smaller scales is accounted for via subgrid scale models.

Multiple investigations of different configurations of bluff bodies interacting with boundary layer flows have been performed over the years. In the early years, Murakami et al. (1987) demonstrated the potential of LES to predict a boundary layer flow around a wall-mounted cube, using a very coarse grid resolution given the computational resources of that time. Later on, Shah and Ferziger (1996, 1997) performed new simulations of a boundary layer flow around bluff bodies and showed that LES is accurate enough in predicting various quantities of interest. In a review article, Rodi (2002) compared LES and RANS predictions of channel flows around a cube and demonstrated the potential of LES to predict such complex flows accurately.

In Nozawa and Tamura (2002), an LES study was performed to quantify the mean and fluctuating surface pressures on a wall-mounted cube, with a precursor simulation used to generate appropriate boundary layer inlet conditions. They found good agreement with experimental data, which was measured in a boundary layer flow over a smooth surface. Ono et al. (2006) conducted an LES study of flows generated by a building that faced the incoming flow at 45 degrees. While, Zheng et al. (2012) used a combination of a RANS simulation with the Fluent inherent random flow generation to generate a turbulent boundary layer flow in order to compare the wind loads on a single tall building calculated from a LES study to wind tunnel experiments. Thomas and Williams (1999) used LES of a surface-mounted cube in a turbulent boundary layer to study the effect of varying cube orientation on the surrounding flow. Sohankar et al. (2000) conducted an LES study of a square cylinder in uniform flow to compare different subgrid scale models. In order to study the effects of varying anisotropy in turbulent flow on the mean flow around a square cylinder, Haque et al. (2014) compiled a library of grid-generated isotropic turbulence flow slices from a wind tunnel that were rescaled and used as inflow for the large eddy simulations. Focusing on the flow in the immediate vicinity of a bluff body, Yakhot et al. (2006) carried out a direct numerical simulation (DNS) of a surface-mounted cube in a fully-developed turbulent channel. With the goal of comparing to wind tunnel experiments, Lim et al. (2009) performed a wall-modeled LES study of a surface-mounted cube interacting with a turbulent boundary layer. De Stefano et al. (2016) used wall-resolved adaptive wavelet-collocation LES of a square cylinder in uniform flow to show that a variable threshold method is feasible for turbulent flows around bluff bodies. A grid resolution study was performed by Tseng et al. (2006) using a wall-modeled LES framework of both a square cylinder and a periodic matrix of cubes in uniform flow. These results were then applied to particulate transport through a urban environment, which was modeled as a configuration of different sized boxes to resemble a cluster of buildings. Bou-Zeid et al. (2009) then conducted a building resolution study using wall-modeled LES to investigate the effect the level of building detail had on the surrounding neutral atmospheric boundary layer (ABL).

An LES approach and the fully three-dimensional proper orthogonal decomposition (POD) technique are utilized to investigate how the dynamics of a stably-stratified turbulent boundary layer are affected by the presence of a large bluff body. The study is focused on the wakes that are the result of the interaction between the bluff body and the boundary layer flow. The numerical tool, consisting of a pseudo-spectral LES algorithm, solves the incompressible Navier-Stokes equations with the Boussinesq approximation in the momentum equations introduced to take into account the buoyancy effect. The effect of the body on the boundary layer flow is taken into account via an immersed boundary method. One of the main challenges of the study is the imposition of a realistic turbulent conditions in the upstream of the body, and this is achieved by utilizing a concurrent precursor simulation performed in a separate flow domain; a blending region is introduced at the end of the main simulation to transfer the data from the precursor simulation. This study is one of the first attempts to characterize the wake generated by the bluff body in high Reynolds number turbulent boundary layer by a fully three-dimensional POD technique. Whereas the body height Reynolds number is on the order of $10^{4}$ to $10^{5}$ in the above studies, this work considers a Reynolds number on the order of $10^{7}$. Results consisting of instantaneous velocity, potential temperature, and POD modes for streamwise and wall-normal velocity components contour plots along with vertical profiles of vertical turbulent momentum and heat fluxes and turbulence intensity are reported and discussed.

2 LES framework

The governing equations employed are the filtered momentum, potential temperature transport, and continuity equations for incompressible flow. Included in the momentum equations is the Boussinesq approximation term.

[TABLE]

[TABLE]

$\tilde{u}_{i},i=1,2,3$ represents the velocity components along the axial, spanwise, and vertical ($x_{1}$-,$x_{2}$-,$x_{3}$-) directions, respectively. The potential temperature is $\tilde{\theta}$, while $\theta_{0}$ is a characteristic potential temperature. The bracketed angles represent a horizontal averaging, $g$ represents the gravitational acceleration, and the Kronecker delta is given by $\delta_{ij}$. The effective pressure, which is divided by density, is represented by $\tilde{p}^{*}$. $F_{i}$ is a forcing term. A Lagrangian scale-dependent model, developed by Bou-Zeid et al. (2005) and later extended by Porté-Agel et al. (2000) to scalar transport, is used to model the subgrid scale (SGS) stress, $\tau_{ij}$, and the SGS heat flux, $\pi_{j}$.

Because the Reynolds number associated with the ABL is very large, the flow in the proximity to the ground is modeled using Monin-Obukhov similarity theory (Monin and Obukhov (1954)), and the molecular viscous diffusion term in the momentum equation is neglected because it is assumed to be very small. The presence of a bluff body in the boundary layer flow is modeled using the direct forcing immersed boundary method (IBM) first introduced by Mohd-Yusof (1997). In the direct forcing approach, the velocity boundary conditions are imposed by a force added to the discretized momentum equations.

[TABLE]

The imposed force can then be calculated by specifying the velocity at the boundary $(u_{i})_{box}$.

[TABLE]

Through substituting the imposed force into the discretized momentum equations, it is shown that the boundary condition,$u^{n+1}_{i}=(u_{i})_{box}$, is satisfied at each discrete point.

Under the assumption that the turbulence in a ABL is horizontally-homogeneous, periodic boundary conditions are used in the horizontal directions. However, here a concurrent precursor simulation provides inflow boundary conditions that are blended at the end of the domain, while keeping the periodicity condition in the streamwise direction. The vertical gradients of velocity and scalars and the vertical component of velocity are set to zero at the top boundary. The top of the domain is located well above the height of the boundary layer. The wall stress is related to the horizontal velocity components at the first grid point via the Monin-Obukhov similarity theory at the bottom boundary (Businger (1971)). While the similarity theory was developed for averaged quantities, Moeng (1984) imposed the wall stress in a strictly local sense. Bou-Zeid et al. (2004, 2005) then showed that by using local filtered velocities, the local similarity formulation produces an average stress that is very close to the averaged similarity formulation.

[TABLE]

The wall stress components are $\tau_{13}|_{z=0}$ and $\tau_{23}|_{z=0}$ and the friction velocity is $u_{*}$. $\kappa$ is the von Kármán constant (taken to be $\kappa=0.4$) and $z_{0}$ is the effective roughness length. $\Psi_{M}$ provides a stability correction and $V_{f}=\left[\tilde{u_{1}}(\Delta z/2)^{2}+\tilde{u_{2}}(\Delta z/2)^{2}\right]^{1/2}$ is the filtered horizontal velocity at the first grid point. Within the same theory, the heat flux at the surface is computed as

[TABLE]

where $\theta_{s}$ is the imposed potential temperature at the surface, $\tilde{\theta}$ is the resolved potential temperature in the first grid point, $z_{0s}$ is the roughness length for scalars, $\Psi_{H}(\zeta)=\int_{0}^{\zeta}[1-\phi_{H}(\zeta^{\prime})d\zeta^{\prime}/\zeta^{\prime}]$ provides a stability correction where $\zeta=z/L$ and $\phi_{H}\left(\zeta\right)=Pr_{t}+\beta\zeta,$ for stable conditions (here, the Prandtl number is $Pr_{t}=0.74$, and $\beta\approx 5$). An artificial heat source or sink, developed by Sescu and Meneveau Sescu, A. and Meneveau (2013), is added above the top of the ABL within the precursor simulation to enable the desired atmospheric stability.

The numerical method implemented in a pseudo-spectral LES code uses a pseudo-spectral horizontal discretization and a centered finite difference vertical discretization (Calaf et al. (2010), Sescu and Meneveau (2015)). The aliasing errors introduced through the handling of nonlinear terms by the pseudo-spectral method are corrected according to the $3/2$ rule (Canuto et al. (1988)). The continuity equation is satisfied through the solution of the Poisson equation resulting from taking the divergence of the momentum equation. In the vertical direction, a staggered grid is used to avoid decoupling the velocity and pressure. This is implemented by storing the vertical velocity component at a location halfway between the location of the other dependent quantities. The solution is marched in time using a fully-explicit second-order Adams-Bashforth scheme (Butcher (1971)).

3 Concurrent precursor simulation

Turbulent inflow conditions at the inflow boundary can be imposed using a precursor simulation, wherein a library of turbulent data is generated by a separate simulation performed prior to the running of the main simulation (Tabor and Baba-Ahmadi (2010)). The library of turbulent data is created in a periodic domain by saving a slice of the flow variables normal to the streamwise direction at every time step to disk. These slices are imposed as inflow boundary conditions for the main flow simulation. In this way the main simulation is provided with a realistic turbulent velocity field and scalars at the inlet (for example, Park et al. (2014), Wu and Porté-Agel (2012, 2013)). The memory that is necessary to store all of the turbulent data and the time used to read and impose the turbulent slices at the inlet of the main simulations increases as the number of grid points increases (Dhamankar et al. (2015)).

Another approach is to run the precursor and main simulations simultaneously; this method was first introduced by Lund et al. (1998) in the context of turbulent flat plate boundary layers and then later improved upon by Ferrante and Elghobashi (2004). In the method of Lund et al., a slice of flow data taken at a streamwise location of the precursor simulation is rescaled and imposed at the inlet of the main simulation. Mayor et al. (2002) used a recycling region upstream of the domain where perturbations at the end of the recycling region are added to a mean flow from a previous simulation. The outlet of the recycling region serves as the inflow conditions for the domain and the combination of the perturbations and the previous mean flow serve as the inflow for the recycling region.

Consisted with the approach proposed by Stevens et al. (2014), the precursor and main flow domains considered here are identical, except the bluff body is added to the main domain. To avoid time interpolation the two simulations are synchronized in time. After each time step, a region of the flow data near the outlet of the precursor domain is blended on to the flow data in a region located in the proximity of the outflow boundary of the main domain. The blending region ensures that the flow data is smoothly transferred from the precursor simulation to the main simulation. Since the precursor domain and the main domain are identical, except for the presence of the body in the main domain, the Stevens method is capable of providing accurate turbulent inflow conditions. However, since the two domain have similar sizes, there is an increase in computational time and memory compared to other techniques, such as those based on recycling and rescaling. Assuming that the length of the blending region is $L_{blend}$ and ranges from $x=L_{s}$ to $x=L_{x}$, where $L_{x}$ is the length of entire domain, a generic flow variable (velocity or temperature) in the blending region can be calculated according to

[TABLE]

where $\phi^{pre}_{i}$ and $\phi^{main}_{i}$ are flow variables from the precursor and main domains, respectively. The blending function used is

[TABLE]

where $L_{pl}=L_{x}-\frac{1}{4}L_{blend}$. Figure 1 shows an example involving an LES of a turbulent boundary layer with inflow conditions obtained from a concurrent precursor simulation; contours of the streamwise velocity are plotted for the precursor and main simulations. The transition inside the blending region can be seen at the end of the main simulation in figure 1.

4 Proper orthogonal decomposition

Proper orthogonal decomposition (POD) is a modal decomposition technique used to determine the highest energy associated with coherent structures in a given flow. In a typical POD analysis, the flow is decomposed into orthogonal modes, with the corresponding eigenvalues representing contributions of the modes to the flow turbulent kinetic energy. Lumley (1967) was the first to apply POD analysis to turbulent flows. Later, the snapshot POD analysis was introduced by Sirovich (1987) to deal with large flow data resulting from three-dimensional simulations. As examples of previous studies, Huang et al. (2009) used a two-dimensional POD analysis to study the coherent structures in ABL, and Liu et al. (2012) analyzed the vertical velocity profiles of the ABL using one-dimensional snapshot POD method to investigate mesoscale thermal structures. Shah and Bou-Zeid (2014) simulated the ABL over a range of stability conditions and used a two-dimensional snapshot POD method to look at the effect buoyancy had on the very large scale flow structures. Bastine et al. (2015) used the two-dimensional snapshot POD analysis of a single wind turbine wake to work towards a simplified dynamic wake model. Recently, VerHulst and Meneveau (2014) employed a three-dimensional snapshot POD method to study the flow within large wind farms; this framework is utilized in the current work. For an in-depth overview of modern modal decomposition methods the reader is advised to read the recent review article by Taira et al. (2017) (where various tools are discussed and presented, and strengths and weaknesses of each modal decomposition method are discussed).

In the snapshot POD analysis, the decomposition of the velocity field is defined as

[TABLE]

where $\overline{u}_{i}$ is the time averaged velocity, the number of snapshots is represented by $N$, the time-coefficients are $a^{k}$, and $\psi^{k}_{i}$ are the POD modes. The POD modes in the decomposition are orthonormal to each other, that is, $\langle\psi^{k}_{i}\psi^{l}_{i}\rangle_{xyz}=\delta_{kl}$, and the time-coefficients satisfy the relationship $\overline{a^{k}a^{l}}=\lambda^{k}\delta_{kl}$, where $\lambda^{k}$ represents the strength of the POD mode $k$. The contribution of the POD mode $k$ to the time averaged turbulent kinetic energy is represented by $\lambda^{k}$.

Snapshots of the velocity field are saved every $\tau$ which is a time scale of the same order as the correlation time. The definition of the $n^{th}$ snapshot is

[TABLE]

and it is assumed that none of the snapshots are correlated with each other. The fluctuating velocity field associated with each snapshot can be determined by subtracting out the mean velocity.

[TABLE]

Invoking the ergodicity assumption and using the fluctuating velocities, the two point correlation matrix is determined as

[TABLE]

where $n,m=1,2,3,\dots,N$. The eigenvalue problem is solved to calculate the POD mode energies, $Ce=\lambda e$, which together with the correlation matrix eigenvectors $e$ can be used to calculate the time coefficients as

[TABLE]

Finally, the POD modes can be calculated as

[TABLE]

Some important properties regarding the POD method can be highlighted. First, the POD modes are arranged according to decreasing mode energy. This implies that a small number of POD modes can characterize the most energetic coherent structures of a turbulent flow. Secondly, POD modes can be used to reconstruct any velocity field snapshot. Finally, the resolved turbulent kinetic energy can be calculated according to

[TABLE]

Because the ergodicity assumption, a long enough time span is necessary for the collection of the velocity field snapshots. To this end, Sirovich (1987) proposed a nominal criterion that implies that the number of required snapshots should resolve at least $99\%$ of the turbulent kinetic energy.

[TABLE]

5 Results and discussions

5.1 Comparison to experimental results

For the purpose of validating the numerical methodology, a boundary layer flow developing over a wall-mounted bluff body is considered, at a Reynolds number $Re=3.0\times 10^{4}$ based on the height of the cube and the freestream velocity. In this test case, the viscous term in the momentum equation (1) is included because the viscous effects may be non-negligible. The objective of this numerical simulation is to validate the LES results with the experimental measurements taken in a wind tunnel data by Castro and Robins (1977). A $0.2$ m wall-mounted cube in a domain with streamwise, spanwise, and vertical dimensions of $1.92\times 0.8\times 2.7$ m is considered, where the height of both the domain and the cube are equal to the height of the wind tunnel test section and the experimental cube. The (streamwise, spanwise, vertical) direction names correspond to the ($x$,$y$,$z$) axes with ($u$,$v$,$w$) velocities in those directions. A representation of the physical domain is shown in figure 2. The main domain is discretized uniformly using $192\times 64\times 324$ grid points in the streamwise, spanwise, and vertical directions, respectively, which gives $20\times 16\times 24$ grid points that are associated with the cube. The dimensions and discretization of the precursor domain are identical to the main domain. A $2.0$ m turbulent boundary layer with a freestream velocity of $2.02$ m/s is imposed in the precursor simulation using time averaged streamwise velocity, turbulence intensity, and vertical turbulent momentum flux profiles given by Castro and Robins (1977) and the controlled forcing scheme developed by Spille-Kohoff and Kaltenbach (2001). The nondimensional time step is $\frac{\Delta tU_{g}}{h_{box}}=0.005$ and the nondimensional grid spacings are $\frac{\Delta x}{h_{box}}\times\frac{\Delta y}{h_{box}}\times\frac{\Delta z}{h_{box}}=0.05\times 0.0625\times 0.0417$. ($l_{box}$, $w_{box}$, $h_{box}$) are the (length, width, height) dimensions of the box and correspond to the (streamwise, spanwise, vertical) directions or ($x$,$y$,$z$) axes.

The time averaged streamwise velocity and turbulence intensity profiles along the vertical direction are plotted in figures 3 and 5, respectively, alongside the experimental data from Castro and Robins at three streamwise locations in the wake, $x/h_{box}=\{0,1,2\}$. Figure 4 shows a spanwise profile of time averaged streamwise velocity at one streamwise location in the wake, $x/h_{box}=0.75$. The agreement was found to be good for all profiles.

5.2 Flow around large bluff bodies

A large rectangular body with streamwise, spanwise, and vertical dimensions of $196\times 79\times 70$m, respectively, is considered in a stratified boundary layer flow. The body is first considered at rest, and then in motion at velocity $10$ m/s in a wind in the opposite direction of velocity of $8$ m/s. This results in a total wind velocity of $18$ m/s for the moving body case. The Reynolds numbers based on the height of the box for the stationary and moving box are both on the order of $10^{7}$. Two small thermal stratifications are considered: $\Delta\theta=0.2$ K and $\Delta\theta=0.4$ K. The precursor and the main flow domains have the same streamwise, spanwise, and vertical dimensions of $3770\times 628\times 600$ m. The nondimensional time steps for the stationary and moving boxes are $\frac{\Delta tU_{g}}{h_{box}}=\{0.0023,0.0051\}$, respectively. Figure 6 shows a model of the physical domain. Uniform grid spacing is used in the horizontal directions and also along the vertical direction because Monin-Obukhov similarity theory is used to model the flow at the wall. The nondimensional grid spacings for the $384\times 64\times 128$ point grid are $\frac{\Delta x}{h_{box}}\times\frac{\Delta y}{h_{box}}\times\frac{\Delta z}{h_{box}}=0.14\times 0.14\times 0.06$. The difference in potential temperature between the ground, $\theta_{s}$, and the top of the ABL, $\theta_{t}$, is used to set the thermal stratification, $\Delta\theta=\theta_{t}-\theta_{s}$.

5.2.1 Contour plots and vertical profiles

Several contours of instantaneous streamwise velocity and potential temperature are presented in figures 7 and 8. In particular, figures 7(a) and 7(b) illustrate the streamwise component of velocity using contours through $xz$ and $xy$ planes, respectively, sectioning the box through the center. It shows that the wake is more intense in the proximity to the box, about two to three box lengths in the downstream. Since the thermal stratification is small, large turbulent structures, in the same order of magnitude as the size of the box, are moving with the flow in the downstream direction (they interact with the box and and the wake). In figures 7(c) and 7(d), corresponding to a wind velocity of $18$ m/s, one can notice that the wake becomes more intense and persists for a longer distance. They also show that the length of the wake along the spanwise direction is larger than the length of the wake generated by the box moving against a $8$ m/s wind. In figures 8(a) and 8(b), contour plots of the potential temperature through an $xz$ plane and an $xy$ plane, respectively, sectioning the box in the center, are depicted. The thermal wake seems to be much longer than the momentum wake (light red and yellow are associated with larger temperatures). Similar contour plots of potential temperature are shown in figures 8(c) and 8(d) for the wind velocity $18$ m/s. One can notice that the thermal wake is more compact in the case of the larger incoming velocity.

The vertical profiles shown in the next figures 9-14 represent time averaged flow data that was plotted in the downstream of the box at every two box lengths from the center of the box. In all figures, the vertical spatial coordinate is scaled by the height of the box. The streamwise time averaged turbulence intensity is revealed in figures 9 and 10 for various cases. In figure 9, one can notice that by increasing the wind velocity, the turbulence intensity in the proximity to the box is enhanced starting from the bottom boundary to approximately one and a half box heights. As the streamwise distance is increased, the turbulence intensity difference between the two cases diminishes. The increase in the thermal stratification slightly impacts the vertical distribution of the turbulence intensity as shown in figure 10. Only very small differences can be noticed in the vicinity of the box (this is because the thermal stratification is very small).

Figures 11 and 12 present profiles of the vertical turbulent momentum flux at four locations downstream of the body. In figure 11, the thermal stratification is kept constant at $0.2$, and averaged profiles from simulation cases with wind velocities $8$ m/s and $18$ m/s are compared to each other. One can notice that by increasing the freestream velocity the vertical momentum flux magnitude decreases, with a maximum attained in the vicinity of the top of the box. A small increase in the momentum flux can be noted in the proximity to the ground. There are clear differences between the profiles for the two wind velocities beyond seven box lengths downstream from the box. In figure 12, the wind velocity is constant at $18$ m/s and averaged profiles for the thermal stratifications of $0.2$ and $0.4$ are compared to each other. Only very small differences are noted in the vicinity of the box.

In the next two figures, comparisons in terms of the vertical turbulent heat flux are presented. Profiles of vertical turbulent heat flux are shown in figure 13 for a stratification of $\Delta\theta=0.2$ K in order to compare the both wind velocity cases. The magnitude of the vertical turbulent heat flux behind the box is larger in the case of the larger freestream velocity, up to approximately $x/l_{box}=3$; at this location, there is a decrease in the turbulent flux above the box. By seven box lengths in the downstream, there is no significant differences between the two profiles. Figure 14 shows that there are no noticeable differences in the heat flux profiles between the two thermal stratifications.

5.2.2 POD results

The snapshot POD analysis was performed utilizing the same tool used in VerHulst and Meneveau (2014). The three-dimensional velocity field was saved every $0.2$ seconds for a total of 8000 snapshots. When first inspecting contours of the streamwise velocity POD modes, three main mode types were identified in the region of the wake: helical, symmetric, and combined helical-symmetric. Figures 15(a), 15(b), and 15(c) show examples of the three types of POD modes, while figures 15(d) and 15(e) show the helical and symmetric modes, respectively, using contours cross-sectioning the wake.

In figure 16(a), the energy levels of first twenty POD modes normalized by the highest energy containing mode are shown for all cases. The relative mode energy levels decrease rapidly over the first twenty POD modes and by the twentieth mode, the corresponding energy is less than half of the energy contained in the strongest mode. So, the behavior of these first modes is important in representing the most energetic flow structures. Figure 16(b) plots the relative energy levels of all the modes that contain at least one percent of energy of the strongest POD mode on a log-log scale. Table 1 highlights the mode number at the one percent threshold. When comparing at constant thermal stratification, the initial modes of the lower wind velocity cases contain more of the relative energy than the higher wind velocity cases until modes $k=90$ and $k=35$ for the $\Delta\theta=0.2$ K and $\Delta\theta=0.4$ K thermal stratification levels, respectively. After those crossover modes, the modes of the higher wind velocity cases are at higher relative energy levels. Increasing the wind velocity spreads the turbulent kinetic energy over a larger number of POD modes. This can be seen in the mode number of the one percent threshold; there are $1.9$ times, for $\Delta\theta=0.2$ K, and $2.7$ times, for $\Delta\theta=0.4$ K, more modes with energy levels above the one percent threshold. Holding wind velocity constant, increasing the thermal stratification raised the relative energy levels of all the modes.

Notice that in figure 16(a), certain consecutive modes occur at nearly equal relative energy levels: $k=1$ and $k=2$ for $\Delta\theta=0.4$ K and $U_{g}=8$ m/s and $k=9$ and $k=10$ for $\Delta\theta=0.4$ K and $U_{g}=8$ m/s. These consecutive modes are a pair of similar helical modes differentiated by a wavelength shifting. In order to preserve the orthogonality of the two modes, the wavelength of one of the modes is shifted by a quarter of the wavelength of the other mode. Figures 17(a) and 17(b) show two examples, at different wavenumbers, of the helical mode pairing. For each pair, the wavenumbers are the same, but the modes are only shifted by a quarter wavelength. By looking at figure 16(a), the number of paired helical modes at the highest relative energy levels increases with increasing wind velocity and thermal stratification. Through this analysis of the significant energy contributing POD modes, the helical mode was found to be dominant.

Next, contour plots of several streamwise velocity POD modes are compared at constant relative energy levels. The horizontal lines in figure 16(a) highlight the mode numbers across all four cases, within the first twenty modes, that have constant relative energy levels, $\frac{\lambda^{k}}{\lambda^{1}}$, which are shown in Table 2. The contours of the streamwise velocity POD modes shown in figure 18 correspond to the relative energy levels in table 2. For example, the second contour from the top in figures 18(a), 18(b), 18(c), and 18(d) corresponds to the second relative energy level in table 2. For the thermal stratification $\Delta\theta=0.2$ K and wind velocity $U_{g}=8$ m/s contours shown in figure 18(a), the modes of the first four energy levels are of the weak combined helical-symmetric type that shows no significant effect from the box. The fifth energy level mode shows the behavior of the dominant helical type, although in a weaker form. The contours of the modes at a thermal stratification $\Delta\theta=0.2$ K and wind velocity $U_{g}=18$ m/s in figure 18(b) show a clear dominant helical mode at the third energy level, while the modes of the first two energy levels are largely unaffected and the modes of the last two energy levels show a stronger combined helical-symmetric type behavior. Figure 18(c) shows the mode contours for the thermal stratification $\Delta\theta=0.4$ K and wind velocity $U_{g}=8$ m/s. Clear dominant helical modes are seen at the highest two energy levels. The next two energy levels display modes with a weak combined helical-symmetric behavior and the mode at last energy level exhibits weaker helical behavior at a larger streamwise wavenumber than the modes at the first two energy levels. The contours of the modes at a thermal stratification $\Delta\theta=0.4$ K and wind velocity $U_{g}=18$ m/s are seen in figure 18(d). The dominant helical modes are still seen at the first two highest energy levels, except stronger, while the last three energy levels all display clearer dominant helical type mode behavior.

Increasing the wind velocity causes the dominant helical modes to appear at higher energy levels. This is seen in the $\Delta\theta=0.2$ K cases where there is a weaker helical mode at the fifth energy level when $U_{g}=8$ m/s and a stronger helical mode at the third energy level when $U_{g}=18$ m/s and also in $\Delta\theta=0.4$ K cases where helical modes appear at the last three energy levels when $U_{g}=18$ m/s, but only at the last energy level when $U_{g}=18$ m/s. The $\Delta\theta=0.4$ K cases are also interesting because clear helical modes are present at the highest two energy levels for both wind velocities, with the modes associated with the higher velocity appearing stronger than the modes associated with the lower velocity at a constant energy level. The appearance of helical modes at higher energy levels with increasing wind velocity and the strength of helical modes at constant energy levels increasing with increasing wind velocity can be related back to the increase of turbulence intensity and vertical turbulent momentum flux in the wake found in figures 9 and 11. Since the helical modes are the dominant modes, as the wake becomes more turbulent these helical modes contain more of the turbulent energy. Increasing the thermal stratification causes a similar effect as increasing the the wind velocity; the dominant helical modes appear at higher energy levels at higher thermal stratifications. Looking back at figures 10 and 12, the profiles in the wake of the turbulence intensity and the vertical turbulent momentum flux were only slightly affected by varying the thermal stratification. This can be interpreted as the increase in thermal stratification concentrating the turbulent energy in the higher energy levels, which manifests as more helical modes at higher energy levels, while leaving the overall statistics largely unaffected.

While the previous comparisons of the POD modes were made at constant relative energy levels, next we focus on comparisons of modes at constant streamwise wavenumber. Contours of the spanwise velocity POD modes are presented in figure 19. These contours show the first occurrence of a dominant helical mode at a constant streamwise wavenumber for each of the four cases. For a constant streamwise wavenumber and comparing at a constant thermal stratification, figures 19(a) and 19(b) for $\Delta\theta=0.2$ K and figures 19(c) and 19(d) for $\Delta\theta=0.4$ K, it is seen that the strength of the helical mode increases with increasing wind velocity. This agrees with the trends seen above. Now, comparing at a constant wind velocity, figures 19(a) and 19(c) for $U_{g}=8$ m/s and figures 19(b) and 19(d) for $U_{g}=18$ m/s, it is seen that the mode inclination angle (with respect to the bottom boundary) decreases a small amount as thermal stratification increases. The increased thermal stratification forces the modes to spread more in the horizontal directions and less in the vertical direction. This is consistent with the small decrease in vertical turbulent momentum flux around the box height seen in figure 12; increasing the thermal stratification “contains” the wake closer to the wall.

The Strouhal number based on the length of the box for the dominate helical modes for each case is shown in table 3. It can be observed in figure 18 that an increase in wind velocity causes the wavelength of the of the helical mode to decrease, which in turn causes an increase in Strouhal number. The decrease in wavelength occurs because the helical modes contain more higher energy structures at the higher wind velocity. Also from figure 18, the wavelength of the helical mode increases as the thermal stratification increases, which yields a decrease in Strouhal number. This relates back to the increased thermal stratification limiting the vertical spread of the modes. The smaller thermal stratification allows for a larger increase in the Strouhal number as the wind velocity increases because the modes are allowed to spread more in the vertical direction.

6 Conclusions

Large eddy simulations were used to study the dynamics of the wake created by a large bluff body in a high Reynolds number stably thermally-stratified boundary layer. The boundary conditions at the surface of the body were modeled using a direct forcing immersed boundary method. Since the LES tool uses a pseudo-spectral discretization, requiring the use of periodic boundary conditions in the horizontal directions, the inflow conditions were handled via a concurrent precursor simulation method aimed at generating realistic turbulent inflow conditions. This allowed the simulation of a single row of boxes (the periodic condition in the spanwise direction was applied). The three-dimensional snapshot POD method was used to study the structure of the wake generated by the body. Two different thermal stratifications and two wind velocities were considered in the analysis. The computational method was validated using experimental results from a wind tunnel flow around a wall-mounted cube.

A large bluff body in a stable boundary layer was then considered, and results consisting of vertical profiles of vertical turbulent momentum and heat fluxes and turbulence intensity were reported and discussed. The turbulent statistics were seen to be significantly affected by the presence of the bluff body. The thermal stratification variation was found to have just a slight effect on all statistical quantities, but the larger wind velocity significantly affected all the vertical profiles when compared to the smaller wind velocity. From the POD analysis, two distinct POD modes were observed in the wake. These two classifications of modes were labeled helical and symmetric and had different streamwise wavenumbers. The helical mode was found to be dominant; the more turbulent the wake became, the more helical modes were found at the highest energy levels. Contrary to the slight effect seen in the statistical quantities, the distribution of helical modes at the highest energy levels was significantly affected by varying the thermal stratification.

Future work will focus on an improvement of the current concurrent precursor method in order to provide a more natural transition from the main to the precursor flow at the outlet of the main domain. In addition, the concurrent precursor method will also be compared with a synthetic eddy method, also in development, in an attempt to improve the computational speed while still generating realistic turbulent inflow conditions because the current concurrent precursor method is computationally expensive.

Acknowledgement

The authors would like to thank Dr. Claire VerHulst for Johns Hopkins University for her help in regard to the POD analysis.