Improved delayed detached-eddy simulation and proper orthogonal decomposition analysis of turbulent wake behind a wall-mounted square cylinder

The turbulent flow past a wall-mounted square cylinder with an aspect ratio of 4 was investigated with the aid of Spalart–Allmaras improved delayed detached-eddy simulation and proper orthogonal decomposition (POD) analysis. The Reynolds number was equal to 12 000 (based on the free-stream velocity and obstacle width). The boundary layer thickness was ∼0.18 of the obstacle height. This study focused on analyzing the vortical structure of the wake and vortex shedding process along the obstacle height. A quantitative comparison of the first- and second-order flow statistics with the available experimental and direct numerical simulation data was used to validate the numerical results. The numerical model coupled with the vortex method of generating the turbulent inflow conditions could successfully reproduce the flow field around and behind the obstacle with commendable accuracy. The flow structure and vortex shedding characteristics near the wake formation region were discussed in detail using time-averaged and instantaneous flow parameters obtained from the simulation. Dipole type mean streamwise vortices and half-loop hairpin instantaneous vortices with energetic motions were identified. A coherent shedding structure was reported along the obstacle using two-point correlations. Two types of vortex shedding intervals were identified, namely, low amplitude fluctuations (LAFs) and high amplitude fluctuations (HAFs) [P. Sattari et al., Exp. Fluids 52, 1149–1167 (2012)]. The HAFs’ interval exhibits von Karman-like behavior with a phase difference of ∼180°, while the LAFs’ interval shows less periodic behavior. It was observed that the effect of the LAFs’ interval tends to weaken the alternating shedding along the obstacle height. The POD analysis of the wake showed that for the elevations between 0.25 and 0.5 of the obstacle height, the first two POD modes represent the alternating shedding and contribute to 66.6%–57.6% of the total turbulent kinetic energy. However, at the free end of the obstacle, the first two modes have a symmetrical shedding nature and share 36.5% of the kinetic energy, while the rest of the energy is distributed between the alternating and the random shedding processes. A simple low-order model based on the vortex-shedding phase angle and the spectrum of the time coefficients obtained from POD was developed to predict the wake dynamics at the range of elevations where the alternating shedding is dominated.


I. INTRODUCTION
The flow around wall-mounted bluff bodies has long been a subject of interest in fluid dynamics, owing to its wide range of applications around engineering structures. For instance, understanding the dynamic behavior of the vortices that are induced by building structures is a significant factor in building design. This helps in avoiding any overlap of the vortex shedding frequency and natural frequency of the structures, which may result in a potential structural resonance. Additionally, understanding the behavior of the wake dynamics can lead to a better comprehension of the air pollution propagation behind buildings, chimneys, and so on. However, this type of flow has a complex behavior as it is highly three-dimensional and contains a wide range of flow regimes (Wang et al., 2006). These flow characteristics require appropriate choice of a suitable experimental setup and numerical model, which can accurately capture most of the flow features, which is a challenging topic in the fluid dynamics and wind engineering community. Numerous experimental and numerical studies have investigated the flow around finite cylinders with circular, square, and rectangular cross sections aiming to obtain a better understanding of the wake behavior behind the obstacle with various Reynolds numbers (Re), aspect ratios (ARs) (the ratio between the height and the diameter or width of the cylinder), and thicknesses of the boundary layer.
Several previous studies have demonstrated that the dynamic behavior of the flow behind the obstacle is strongly affected by a specific range of aspect ratios and Reynolds numbers for both developing and fully developed boundary layers. Hence, quasi-periodical, von Kármán-like vortex shedding dominates over a certain range of Reynolds numbers and aspect ratios. If the aspect ratio is less than the prescribed limit, this vortex vanishes due to the effect of tip-induced downwash (Sumner et al., 2004;Wang and Zhou, 2009;Williamson, 1996;and Zhang et al., 2017). Wang et al. (2006) conducted an experiment to study the near wake of a finite-length square cylinder using a closed-loop low-speed wind tunnel, hotwire anemometry, laser Doppler anemometry, and particle image velocimetry (PIV). In their study, the aspect ratio of the test cylinder ranged from three to seven. The range of Reynolds number was 221-9300 (based on the free stream velocity and the obstacle width). They recorded two types of primary symmetric and antisymmetric (arch type) vortices. Additionally, they observed that the flow behavior depends strongly on the aspect ratio, while the overall flow structure does not. Okajima (1982) studied the flow around a rectangular cylinder in a wide range of Reynolds numbers between 70 and 2 × 10 4 using both the wind tunnel and water tank. The author observed that within a specific range of Reynolds numbers and aspect ratios, the flow vortical structure and the reattachment points abruptly changed with a sudden discontinuity in the Strouhal number. Lim et al. (2007) investigated the effect of the Reynolds number on the turbulent flow over a cubic obstacle based on wind-tunnel experiments. They recorded that within the range of Reynolds numbers in which the flow did not contain strong concentrated-vortex motions, the mean flow quantities were independent of the Reynolds number, while the fluctuating quantities such as the root mean square (rms) and fluctuating surface pressures were dependent on the Reynolds number. However, with the presence of strong vortex motions, the Reynolds number significantly affected both the mean and fluctuating flow quantities. Martinuzzi et al. (2007) experimentally investigated the turbulent flow around square-based, surfacemounted pyramids mounted in thin and thick boundary layers based on hot-wire measurements. They presented that for both boundary layers, the normalized pressure of the ground plane within the wake region can be scaled with respect to the length measured from the upstream origin of the separated shear layer to the near wake attachment point. Wang and Zhou (2009) investigated the flow around a square cylinder with AR = 5 at three boundary layer thicknesses. They observed that the probability of anti-symmetrical vortex shedding at z/D = 1 (where z is the distance normal to the ground and D is the cylinder width) was reduced from 84% to 34.5% when the boundary layer thickness increased from 0.07D to 0.245D. However, at z/D = 4, an opposite behavior was observed with the probability of anti-symmetrical shedding varying from 19.5% to 46.5%. These results indicate that by increasing the boundary layer thickness, the base vortex is enhanced, inducing a stronger upwash flow from the cylinder base, which leads to a weakening of the downwash free-end shear layer and the tip vortex. Wang (2012) experimentally studied the near wake of a wall-mounted finite-length square cylinder with AR = 7 and Re = 9300. The proper orthogonal decomposition (POD) was used to analyze the experimental data. The results showed that under the effects of free-end downwash flow, the near wake is highly three-dimensional and drastically different from that of a 2D cylinder and both symmetrical and anti-symmetrical vortices' structures occur in the finite-length cylinder wake. Bourgeois et al. (2011) and Sattari et al. (2012) experimentally investigated the flow characteristics around a wall-mounted square cylinder with aspect ratio = 4 and Re = 12 000 in terms of large scale vortical structures and quasi-periodical shedding flow pattern. They recorded two von Kármán-like alternating vortices and two co-existing vortices in the lee side of the obstacle throughout the shedding cycle observed within low-amplitude pressure fluctuation intervals.
In addition to the progress in experimental studies, a significant improvement has been achieved in the numerical simulations of complex flows over the last three decades, owing to rapid development in computer capabilities and numerical models. Rodi (1997) investigated the ability of large-eddy simulation (LES) and Reynolds averaged Navier-Stokes (RANS) approaches of turbulent flows over a wall-mounted square cylinder at Re = 22 000. The author reported that the turbulence fluctuations were highly underpredicted with the RANS method, while LES could capture considerable details of the unsteady turbulent flow. Fröhlich et al. (1998) performed LES of flow around a circular cylinder at Reynolds numbers 3900 and 140 000 (based on the cylinder diameter and free-stream mean velocity). Their results revealed excellent compliance with the experimental data. They inferred that LES is better suited for simulating this type of flow than the RANS model. Dousset and Pothérat (2010) performed a direct numerical simulation (DNS) study of the laminar shedding of hairpin vortices in the wake of a square cylinder placed in a rectangular duct at Reynolds numbers varying from 10 to 400 and AR = 4. They discussed the development of hairpin vortices as well as the different flow structures that developed with an increase in Reynolds number. Zhang et al. (2017) performed a DNS study of flow past a square cylinder with AR = 4 and low Reynolds numbers from 250 to 1000. In their study, they found that the streamwise vortex structure is significantly affected by the Reynolds number apart from the aspect ratio. They identified three types of mean streamwise vortices, namely, "quadrupole type" at Re = 50 and Re = 100, "six vortices' type" at Re = 150 and Re = 250, and "dipole type" at Re = 500 and Re = 1000. Furthermore, they observed three types of spanwise vortex-shedding models, namely, "hairpin vortex model" at Re = 150, "C, reverse-C and hairpin vortex model (symmetric shedding)" at Re = 250, and "C, reverse-C and hairpin vortex model (symmetric/antisymmetric shedding)" at Re = 500 and Re = 1000.
Motivated by the rapid increase in computational capability, Saeedi et al. (2014) performed a DNS study of the flow around a square cylinder based on the wind-tunnel experiment of Bourgeois et al. (2011) and Sattari et al. (2012), which was considered as a challenging case at the computational fluid dynamics (CFD) Society of Canada Annual Conference in 2012. They achieved an excellent agreement between their numerical results and experimental measurements, and the wake was spread wide behind the cylinder and exhibited complex and energetic vortex motions. Furthermore, the vortex shedding patterns and other flow characteristics such as the coherent structures of the flow downstream of the cylinder, the turbulent kinetic energy, and the instantaneous pressure distribution in the wake region were discussed in their study. Although their study demonstrates the possibility of performing DNS to simulate the flow over obstacles at moderate and relatively high Reynolds numbers, it is still impractical and sometimes impossible to apply DNS for high Reynolds numbers and flows over complex structures.
In the aforementioned studies, LES can be considered as a good CFD model that can simulate the flow over wall-mounted cylinders. Nevertheless, the requirement of an extremely fine mesh near the ground and the walls of the obstacles deem it impractical for high Reynolds number applications. This issue motivated Wang et al. (2019) to use delayed detached-eddy simulation (DDES) in their study on the effect of wall proximity on the flow around a cube with Re = 50 000. Their model validation demonstrated a commendable agreement with the available experimental data. Chen (2018) used the shear stress transport model-improved delayed detached-eddy simulation (k-ω SST-IDDES) to model the flow around a wall-mounted square cylinder with Re = 12 000 and AR varying from 1 to 4. The time-averaged statistical results were in compliance with the previous experimental and LES data. The author reported that the influence of the downwash flow becomes less dominant with increasing AR as accompanied by a near-bed upwash flow at the rear of the obstacle.
In this study, the Spalart-Allmaras improved delayed detachededdy simulation (S-A IDDES) was used to model the flow past a finite square cylinder with AR = 4 at Re = 12 000 in a thin developing boundary layer, and a proper orthogonal decomposition (POD) analysis was applied to the simulation results with the aim of investigating the vortical structure and the shedding process of the turbulent wake behind the obstacle by using the time-averaged and instantaneous results of the simulation as well as the POD modes and the corresponding time coefficients. This paper is organized as follows: In Sec. II, the case parameters, governing equations, turbulent inflow conditions, and the numerical method are presented, which are used to simulate the flow. A brief explanation of proper orthogonal decomposition is presented in Sec. III. The numerical model is validated using a quantitative comparison with the available experimental measurements and DNS data in Sec. IV. The results of this study are thoroughly discussed in Sec. V followed by conclusions in Sec. VI.

A. Case geometry and boundary conditions
The schematic diagram of the numerical domain for the wallmounted square cylinder is shown in Fig. 1. The square cylinder has a height of H = 4D, where D is its width. The Reynolds number was 12 000 (based on the cylinder width and free-stream mean velocity). The overall size of the domain was equal to 29D × 18D × 15D. Therefore, the blockage ratio was ∼1.5%, which was lower than the maximum value of 3% as recommended by Franke et al. (2004). The obstacle was located 8D downstream from the inlet and 20D upstream from the outlet. The distance of the obstacle from each side of the domain was set to 8.5D, while the distance from the end of the obstacle to the top surface of the domain was set to 11D. In large-eddy simulation and detachededdy simulation, the inlet boundary conditions significantly affect the simulation reliability. Hence, the inflow boundary conditions were carefully treated and the vortex method (VM) (Sergent, 2002) was implemented in OpenFOAM-5.x to generate unsteady turbulent inflow conditions. A brief explanation of the method can be found in Sec. II B of this paper. For the outlet plane of the domain, the pressure outlet boundary condition was assigned and the periodic boundary condition was used for the side boundaries. In addition, the no-slip wall boundary condition was applied to the obstacle walls and the ground, while the symmetry plane was applied for the top surface of the computational domain.

B. Turbulent inflow conditions
One of the most important factors that affect the accuracy of the simulation is the inlet conditions of the turbulent flow. The velocity components should contain correlated fluctuating parts with prescribed time and length scales. In this section, the vortex method (VM) of generating turbulent inflow conditions is briefly explained. This method has been further elaborated by Sergent (2002) and Mathey et al. (2006). The vortex method is essentially derived from the Lagrangian form of the two-dimensional evolution equation of vorticity. The fluctuating velocity field in the direction normal to the streamwise of the flow is resolved by using a fluctuating twodimensional vorticity field in the inlet plane, as can be seen in Fig. 2. The resulting discretization for the tangential velocity field is obtained by using the Biot-Savart law as where nx is the unit vector in the streamwise direction, |p − p i | is the distance between the center of the vortex p i and the position p, N represents the total number of vortices injected on the inlet plane, and Γi represents the circulation of the vortex that can be estimated according to the turbulent kinetic energy k of the inlet plane with area S as . (2) The symbol σ in Eq. (1) is a parameter used to control the size of the vortex and can be estimated by using the turbulent mixing length hypotheses, where L is the turbulent length scale, which can be approximated as where Cμ is a model constant equal to 0.09 and ε is the turbulence dissipation rate. The mean velocity profile, turbulent kinetic energy, and turbulent dissipation rate were specified by performing the RANS standard k-ε model simulation of channel flow with the same domain size of the main simulation and with a turbulent intensity of 0.8%. The inlet mean velocity profile of the main simulation was carefully chosen to ensure that the boundary layer thickness is ∼0.18H at the obstacle location, where H is the height of the obstacle. The number of vortices injected on the inlet plane N was set to 100. The perturbation of the velocity in the streamwise direction was estimated using a simplified linear kinematic model (LKM) (Mathey et al., 2006).

C. Governing equations and turbulence modeling
The momentum equation for an incompressible viscous fluid in the absence of external forces is written as Here, u is the velocity of the fluid, ρ is the density, p is the pressure, and v is the kinematic viscosity. The incompressibility of the fluid is expressed by the continuity equation Different numerical approaches have been used to solve and model the flow field parameters in Eqs. (5) and (6). DNS resolves all  scales of turbulent flow structures from the large energy-containing scales to the smallest Kolmogorov length scale in the dissipation range (Kolmogorov, 1941). To ensure this, an extremely fine grid size is required and results in a high computational cost, thus making DNS impractical for several applications, especially those with high Reynolds numbers and complex geometries (Moin and Mahesh, 1998). The RANS equations solve only the time-averaged flow field, i.e., the effects of the fluctuating components are considered through the modeled Reynolds stresses. Thus, the RANS model is not suitable to describe the details of the flow field such as the instantaneous effect of multi-scale eddies (Schmitt, 2007;Wilcox, 1993). In LES, Kolmogorov's first hypothesis is considered, i.e., the small eddies in the turbulent flow tend to exhibit a universal behavior. This property allows the possibility to model small eddies and directly solve large eddies (Smagorinsky, 1963;Stolz et al., 2001). This process can be conducted using spatial filters to separate the large scales from the small ones. The small-scale structures (small eddies) are modeled and introduced as a subgrid-scale (SGS) term in the filtered Navier-Stokes equations. Thus, LES is considerably more accurate for describing the flow characteristics compared with RANS. Nevertheless, it still requires high computational effort, which makes it impractical to simulate the flow around engineering structures at high Reynolds numbers (Wang et al., 2019). Hybrid RANS-LES models, such as the detached-eddy simulation (DES), have been introduced by Spalart and co-authors (Spalart et al., 1997;Spalart, 2000;and Travin et al., 2000) to bridge the gap between LES and RANS. DES has been designed to simulate the flow at high Reynolds numbers and flow around bodies with massive flow separation. By applying RANS near the solid walls within the attached boundary layer and using LES outside the boundary layer in the separated flow region, the grid size for numerical simulations can be reduced significantly in the near-wall region, compared with that required by LES. DES is a promising turbulence modeling tool that can be applied for simulations including wall-bounded flows at high Reynolds numbers. However, some undesirable phenomena such as modeled stress depletion and grid induced separation may occur and affect the model accuracy (Spalart, 2009). This could be attributed to the mesh resolution inside the boundary layer that may result in DES attempting LES rather than RANS in the near-wall region. To overcome these limitations inherited with DES, variant versions such as DDES (Spalart et al., 2006) and IDDES (Shur et al., 2008) have been introduced to overcome these model limitations such that the adapted models can be suitable for different grid spacings inside the boundary layer regardless of the boundary layer thickness. Therefore, considering the relatively high Reynolds number used in this study and the complex flow that is expected to be obtained from the simulation, the S-A IDDES has been chosen to simulate the flow. The transport equation of the S-A IDDES can be defined as whereṽ is the modified eddy viscosity. The turbulent eddy viscosity is defined as vt = f v1ṽ , and the functions f v1 and f w are used as corrections in the near-wall region.S is the strain rate tensor, andr is a non-dimensional term defined as vt/(Sk 2 dw 2 ), where k here is the von Kármán constant and dw is the distance from the wall. σ, c b1 , c b2 , and cw are the model constants of the Spalart-Allmaras model (Spalart and Allmaras, 1992). lIDDES represents the modified length scale used to trigger the scale-resolving mode. The primary idea of IDDES is to switch the transition from URANS to a scaleresolving mode, i.e., LES, depending on a criterion based on the turbulent length scale. The lIDDES term is introduced into the destruction term of Eq. (7) to decrease the turbulent eddy viscosity with increasing wall-normal distance from the wall. This leads to a gradual switch to the scale resolving mode (Spalart, 2009). The lIDDES term is defined as where lLES is defined as CDESψΔ, CDES is a constant equal to 0.65, and ψ is a correction for low Reynolds numbers (Spalart et al., 2006). In addition, Δ is the characteristic cut-off length scale used to calculate lLES, which can be defined as where Cw is an empirical constant equal to 0.15, hmax is the maximum of the local cell size in the streamwise, wall-normal, and lateral directions, and hwn is the wall-normal grid spacing. The functionf d in Eq. (8) is defined asf The switching from URANS accrues inside an intermediate region known as the gray region. Thus, 0 <f d < 1 in this region. Here, f d is called the delaying function and is defined as where r d is inherited from the Spalart-Allmaras model (Spalart and Allmaras, 1992). The blending function f B in Eq. (10) is purely grid dependent and is based on the distance from the wall and the local maximum cell edge length, The grid-dependent parameter α is calculated as The elevating function f e included in Eq. (8) prevents excessive reduction in the Reynolds stresses in the near-wall region and can be defined as where f e1 is the grid-dependent function and f e2 is the function of flow field parameters.

D. Numerical framework
The finite volume open source CFD code OpenFOAM-5.0x was used to perform the simulation. A structured hexahedral type mesh was used for spatial discretization. Local mesh refinement was applied using the stretching mesh technique in regions where high velocity gradients were expected. As shown in Fig. 3, the grid refinement was primarily introduced near the ground and around the obstacle walls. Grid-independence tests were performed using three different grid sizes, namely, 132 ×122 × 90, 152 ×142 × 110, and 172 × 162 × 130. The overall time-averaged drag coefficient and Strouhal number were chosen as the parameters for the test. As shown in Table I, less than 1% variation was observed in the drag coefficient, while there was no variation in the Strouhal number with the progressive refinement of the mesh, ultimately converging the grid-independence test. Considering the maximum dimensionless first-cell spacing y + , the finest mesh size of 172 × 162 × 130 was chosen for this study.  Figure 4 shows the distribution of y + along the ground surfaces and around the obstacle surfaces in the central x-z plane (y/D = 0). The maximum value of y + is less than 0.5 on the ground, while it is ∼15 on the obstacle. It has been observed that the maximum value of y + on the sides does not exceed 2. These values of y + are within the applicable range for IDDES (Spalart, 2009). Advancement was done through a time step of size 5 × 10 −6 s while maintaining the mean Courant number (CFL) at 0.52 and the maximum Courant number less than 0.95 to ensure a stable simulation. To ensure that the flow reached a statistically stationary state, the simulation was first extended over 30 shedding cycles, and then the flow statistics were collected through ∼60 shedding cycles, where the shedding cycle is the time required to complete one full alternative shedding cycle. The PIMPLE algorithm, which is a merge of the PISO algorithm (pressure implicit split operator) and SIMPLE algorithm (semi-implicit method for pressure-linked equations), was employed to solve the coupled pressure momentum system. The convective fluxes were discretized with a second-orderaccurate linear upwind scheme, and all other discretization schemes used in the simulation had second-order accuracy. In order to perform parallel computing, the computational domain was divided into 80 sub-domains and the Scotch method was used to decompose the domains. Correspondingly, 80 processors were used to solve the flow field at each time step. All the computations were performed using a computer cluster with an Intel Xeon Skylake (Gold 6148) 2.4 GHz processor.

III. PROPER ORTHOGONAL DECOMPOSITION
Proper orthogonal decomposition (POD) (Lumley, 1967;Berkooz et al., 1993) was applied to different x-y planes along the obstacle height to investigate the flow behind the obstacle in terms of the coherent structure and the dynamics of the turbulent wake. POD is a technique that decomposes the fluctuating parts of the velocity components into spatial orthogonal modes ϕ m and corresponding time coefficients am(t). The extracted modes represent the most coherent structure of the field, where u ′ here represents the fluctuating part of the streamwise velocity component and M is the number of POD modes. It is worth mentioning here that the direct POD was used in this study, which is slightly different from the snapshots of POD (Sirovich, 1987). The modes are the eigenvectors obtained from solving the eigenvalue problem of the correlation matrix A T A/M − 1, where A is a matrix that contains the fluctuating part of the velocity component, where N represents the position in each of the snapshot (instantaneous) data of the flow and v ′ is the fluctuating part of the spanwise velocity component. The corresponding eigenvalues λm represent the turbulent kinetic energy captured by the respective POD modes. The eigenvalues are arranged from the largest to the smallest (in decreasing order) such that

IV. VALIDATION OF THE NUMERICAL RESULTS
In order to validate the numerical results, the first-and second-order flow statistics obtained from the simulation were compared with the available experimental measurements of Bourgeois et al. (2011) and the DNS data of Saeedi et al. (2014). Figure 5 shows the recovery of the time-averaged nondimensionalized streamwise velocity U/U∞ in the central x-z plane (y/D = 0) for two different elevations behind the obstacle (z/D = 1 and 3). For both elevations, the velocity has negative values near the obstacle wall due to the reverse flow as a result of recirculation. For x/D ≈ 3, the velocity value becomes positive and the flow recovers to the free stream velocity with a higher rate at elevation z/D = 3. Both Figs. 5(a) and 5(b) show a commendable agreement of the results with the experimental measurements and the DNS results. Figure 6 shows the streamwise profiles of the turbulence intensity (urms/U∞) obtained from the simulation, experimental measurements, and the DNS data in the central x-z plane (y/D = 0) for two different elevations (z/D = 1 and 3). The values of turbulent intensity have been well reproduced at the near-wall peak region for both elevations. However, under-prediction of the values in the far downstream of the obstacle at elevation z/D = 3 can be observed in Fig. 6(b). This is attributed to the turbulent level attenuation in this region, which demands more computational effort.
The time-averaged non-dimensionalized streamwise velocity along the lateral direction at elevation z/D = 3 and streamwise location x/D = 2 is shown in Fig. 7(a). The streamwise velocity values become negative at y/D between −0.5 and 0.5, indicating the effect of recirculation inside the wake. The velocity profile obtained from the numerical simulation in general complies with the experimental measurements and DNS results. However, the streamwise velocity values are slightly over-predicted in the recirculation region. Figure 7(b) compares the numerical results with the experimental measurements and the DNS results of the non-dimensionalized Reynolds stress component u ′ v ′ /U 2 ∞ at elevation z/D = 3 along the lateral direction and at the streamwise location x/D = 2. As shown in Fig. 7(b), the Reynolds stress profiles have maximum negative and positive values at y/D = ±1, which decrease rapidly farther in the lateral direction, indicating that the shear stress accrues due to the combination of the vortices induced behind the obstacle. A favorable agreement with the experimental and DNS results can be observed from Fig. 7(b). Here, the S-A IDDES model successfully reproduces the flow field around and behind the obstacle with commendable accuracy, and this could be attributed to the turbulent inflow conditions obtained from the VM and the small blockage ratio used in this study.

A. Flow structure
In order to obtain statistically stationary mean flow results, the turbulence statistics have been collected with a sampling time equal to 60 vortex-shedding cycles. Figure 8 shows the time-averaged streamlines and the non-dimensionalized time-averaged streamwise velocity contours in the x-y planes located at three different elevations. The flow as expected is symmetrical about the central line, i.e., y/D = 0. Figure 8(a) shows that when the flow is very close to the ground (z/D is less than 0.08D), the horseshoe vortex A and the base vortex B are clearly observable and the streamlines assume a bell shape around the obstacle with two circulation regions right behind the obstacle wall. As shown in Figs. 8(b) and 8(c), for higher elevations, the flow behavior is completely different from the near-ground case. The horseshoe vortex does not exist in this range of elevations, while two large counter-rotating vortices E and F are formed behind the obstacle and two small vortices C and D are located near the obstacle sides. Moreover, the boundary of these vortices decreases and becomes closer to the obstacle walls with increasing elevation. Figure 9 shows the time-averaged non-dimensionalized streamwise velocity contours and streamlines in the central x-z plane located at y/D = 0. A big recirculation bubble G can be observed clearly, and as observed from Fig. 8, the boundary of the recirculation bubble decreases toward the obstacle free-end due to the downwash effect of the flow. Figure 10 shows the non-dimensionalized Reynolds stress components in the x-y planes located at three different elevations along  Fig. 10(a), two distinct symmetric spots of normal Reynolds stress u ′ u ′ /U 2 ∞ of high intensity are recognized on the sides of the flow central line and tend to shift toward the obstacle walls with increasing elevation. This behavior is directly related to the variation in the strength of the vortex shedding along the height of the obstacle. Figure 10(b) shows the spanwise Reynolds stress component v ′ v ′ /U 2 ∞ . The maximum value accrues at z/D = 1 (0.25H), which represents the end of the mean recirculation bubble at plane y/D = 0 in Fig. 9. Here, the flow is considered as quasi two-dimensional, i.e., it is primarily affected by the von Kármán process. In Fig. 10(c), as expected, the Reynolds stress u ′ v ′ /U 2 ∞ has two peaks with negative and positive signs with low magnitude compared with the other Reynolds stress components. As u ′ v ′ /U 2 ∞ represents the correlation between the streamwise and the spanwise velocity components, the flow, in general, is dominated by the twodimensional vortex shedding along the three quarters of the obstacle height. However, the flow behavior varies for the last quarter as a result of the tip effect, and this has been discussed in detail in Secs. V B and V C.
To obtain the main vorical structure of the turbulent wake, the time-averaged normalized vorticity component contours at different locations of the x-y and y-z planes are demonstrated in Figs. 11 and 12. As observed in Fig. 11(a), the z-component of the mean vorticity near the ground surface is decomposed into two small antisymmetric vortices that represent together the horseshoe vortex and two big antisymmetric vortices located on the sides and behind the obstacle. It is evident from Fig. 11(a) that the vortices in front of the obstacle switch signs with the main vortices, owing to the adverse pressure gradient in front of the obstacle. As shown in Figs. 11(b) and 11(c), for higher elevations, the horseshoe vortices no longer exist and the two big vortices approach the obstacle walls as z/D increases. Figure 12 shows the streamwise normalized mean vorticity in the y-z planes located at four different streamwise locations. In the region between the rear wall of the obstacle and the center of the recirculation bubble, a pair of counter-rotating tip vortices can be recognized. Additionally, two more counter-rotating vortices originating from the two streamwise-parallel leading corners of the obstacle can be observed in Fig. 12(a). The streamwise spatial evolution of the mean streamwise vorticity beyond the center of the recirculation bubble is shown in Figs. 12(b)-12(d). It can be observed that the two tip vortices disappear directly beyond the center of the recirculation bubble, i.e., x/D = 2, while the horseshoe and base vortices gradually decrease with increasing distance from the obstacle. The leading corner vortices continually descend toward the ground surface until the upwash effect dominates, thus causing a bend in the flow and a remarkable reduction in the vortex pair strength as the streamwise distance increases. Notably, in contrast with the far downstream region, the leading corner vortices at distances inside the recirculation bubble have a vortex tail that formed along the obstacle height, as observed in Figs. 12(b) and 12(c).
The λ 2 criterion (Jeong and Hussain, 1995) is used in this study to visualize the mean and the instantaneous vortical structure of the

ARTICLE
scitation.org/journal/adv wake. λ 2 considers the symmetric Sij and antisymmetric Ωij parts of the velocity gradient tensor ∇Uij, Sij = 1/2(Ui,j + Uj,i), Ωij = 1/2(Ui,j − Uj,i), In Eq. (18), λ 2 is the second eigenvalue of Mij, where i, j, and k are the three components of the Cartesian coordinates. The condition for a pressure minimum associated with the vortex formation is λ 2 < 0, which represents the vortex core at this particular location. Figure 13 shows the isosurface of the mean flow vortex cores identified by λ 2 = −0.3. A single pair of dipole type vortices can be recognized downstream from the obstacle, which is consistent with Fig. 12. This behavior is expected for flow with ARs between 3 and 5. At higher ARs, the quadrupole type vortices are expected due to the downwash being less effective and the increase in the upwash effect (Wang and Zhou, 2009). Figure 14 shows the isosurface of the instantaneous vortical structure identified by λ 2 = −24. The alternating half-loop shedding structure is shown clearly with the hairpin vortex type. The vortex structure has a main core, namely, the "principal core" (Bourgeois et al., 2011) that connects the vortex with the ground surface where the vortex is divided into sub-vortices that overlap with the boundary layer vortices.

B. Shedding process analysis
In this section, the shedding process is investigated using the temporal evolution of the instantaneous velocity, vorticity, and pressure results obtained from the numerical simulation.
As the flow in the vortex formation region is considered inhomogeneous, the two-point correlation function along the obstacle height can be defined as where the angle brackets represent the mean value, α = v ′ or ω Z , and Δz is the distance lagging along the obstacle height. Figure 15 shows the two-point correlations of the spanwise velocity fluctuations R v ′ v ′ and the z-component of the vorticity Rω Z ω Z along the obstacle height for different streamwise locations. As observed in Fig. 15, the correlation coefficient for all the streamwise locations decreases to zero just before the distance lagging equals the obstacle height. This gives an integral length scale ≈1D, indicating a strong spatial correlation along the obstacle height due to alternating vortex shedding.

ARTICLE scitation.org/journal/adv
In order to examine the frequency associated with the shedding process, a power spectral density (PSD) analysis was applied to the spanwise velocity fluctuations at y/D = 1.2 for different heights along the obstacle and different streamwise locations, as shown in Fig. 16. The frequency associated with the periodic shedding is located at the peak of the PSD with a value of St = 0.1 ±0.008, where St is the Strouhal number = f D/U∞ with f being the frequency in Hz. The distribution of the PSD amplitude is slightly broader at z/D = 4 (for all streamwise locations) comparing with the other elevations, indicating a change in the vortex shedding behavior near the end of the obstacle where the flow is expected to be more three-dimensional. However, the obtained PSD peak is maintained coherent along the obstacle height. These results are consistent with the results obtained from the hot wire measurements of Sattari et al. (2012). Figure 17 shows the simultaneous temporal evolution of the pressure coefficient, Cp = (p − p∞)/(1/2ρU 2 ∞ ), on the sides of the obstacle, i.e., at y/D = ± 0.5 and x/D = 0 for different heights along the obstacle. Two intervals can be obtained from Fig. 17, namely, the HAFs' interval and LAFs' interval (Sattari et al., 2012). These two periods correspond to the occurrence of two types of vortex shedding processes. During HAFs, primary counter-rotating vortices with a von Kármán-like shedding behavior are formed and shed alternately from both sides of the obstacle. The phase difference between the vortices on the opposing sides is ∼180 ○ in this interval. The formation of co-existing secondary vortices with opposite signs can be obtained from the LAFs' interval. It is evident from Fig. 17 that the LAFs' interval shows less periodic behavior compared to the HAFs' interval. Notably, even though the interval size of the LAFs for the same elevation changes randomly, the secondary vortices become considerably active as the elevation increases. This is due to the increase in the LAFs' interval size until reaching the free end of the obstacle, where the HAFs have almost the same amplitude as the secondary vortices except that the alternating shedding is kept coherent.
The normalized probability density function (PDF) of the phase difference between the pressure coefficients on the opposite sides of the obstacle for different elevations is shown in Fig. 18. The phase difference broadens as the height increases, which further demonstrates the increasing effect of the secondary vortices along the obstacle height. Nevertheless, the results obtained from

C. POD analysis
POD was applied to the x-y planes located at different elevations along the obstacle height. A total of 500 snapshots spanning by 6 ×10 −4 s were collected for each elevation and subtracted from the mean value. The sample size of the data was equivalent to 36 shedding cycles, which makes it sufficient to obtain converged POD analysis. Figure 19 shows the contours of the first three POD modes of the streamwise and spanwise velocity components at elevation z/D = 1. It can be seen clearly that the first and the second modes accrue alternately, indicating that the first two modes represent the alternative shedding of the wake. However, the third mode exhibits a random and complex behavior comparing with the first two modes. For elevation z/D = 4, the three modes show a completely different behavior comparing with z/D = 1. As shown in Fig. 20, the first and second modes represent the symmetrical shedding of the wake, while the third mode shows here a weak alternative behavior. This change in the behavior of the modes is attributed to the downwash effect from the free-end of obstacle, which tends to weaken the alternative shedding and increase the effect of the symmetrical shedding. Figure 21 shows the normalized eigenvalues of the first 20 POD modes of the x-y planes located at different elevations along the obstacle height. Since the normalized eigenvalues represent the fractional contribution to the total turbulent kinetic energy, as can be seen that for elevation z/D = 1 (0.25H), the first two modes account for 66.6% of the total turbulent kinetic energy, while the fractional energy from the third mode is 2.6%, indicating that the turbulent wake is dominated by the alternating shedding. However, for z/D = 4 (1H), the first two modes account for 36.5% of the total turbulent kinetic energy, suggesting that the turbulent wake at this elevation is highly affected by the symmetrical shedding, which tends to weaken the alternating shedding as a result of the free-end downwash effect.
The simultaneous temporal evolution of the first three normalized time coefficients at z/D = 1 and 4 is shown in Figs  a complex behavior. It is observed from Fig. 22(b) that the clear periodic pattern no longer exists for the first two coefficients. Figure 23 shows the distribution of the normalized time coefficients of the first and second modes for different elevations along the obstacle height. As can be seen from Fig. 23(a), at z/D = 1, the points scatter close to the boundary of the unit circle, indicating that there is a phase difference angle between a 1 and a 2 , which is similar to the correlation between the first two time coefficients of the wake behind an infinite cylinder with a similar Reynolds number (Oudheusden et al., 2005). This observation could lead to a simple approach of building a low-order model for the flow at the elevations where this correlation exists using the periodic nature of the time coefficients of the first two modes as in Eqs. (20) and (21). However, since the contribution to the total turbulent kinetic energy of the corresponding modes is less than the first two modes of the infinite cylinder, the contribution of the other modes should be taken into account, as expressed in Eq. (22), a 2 (t) = √ 2λ 2 cos(2π f t + θ),

ARTICLE scitation.org/journal/adv
where f is the shedding frequency and θ is the corresponding phase shift, where uLOM and vLOM represent the low-order model of the streamwise and spanwise flow velocity, respectively. U and V are the corresponding time-averaged velocities, while the rest of the time coefficients are represented by a ′ m and modeled by using time-series gm(t) of random numbers that have a Gaussian distribution with a zero mean and standard deviation equal to one. In order to enforce a realistic spectral representation of the coefficients, a simple spectral transfer function is applied aŝ whereâ ′ m( f ) andĝm( f ) are the Fourier transform of a ′ m (t) and gm(t), respectively. Ea( f ) is the energy spectrum of the time coefficients obtained from POD, while the energy spectrum of gm(t) is represented by Eg( f ). Since the sets of the time coefficients are orthogonal to each other, i.e., uncorrelated in time, every time-series is modeled separately and normalized to satisfy the condition acceptable agreement with the instantaneous results of the numerical simulation. Despite the simplicity of the low-order model, it could successfully reproduce the main fetchers of the turbulent wake behind the obstacle for the range of z/D = 1-2, i.e., from 0.25 H to 0.5 H. However, for higher elevations, and as can be seen from Figs. 23(c) and 23(d), the clear periodic nature of a 1 and a 2 does not exist and the flow becomes more complex and unorganized, so the low-order model obtained from Eq. (22) is not applicable at this range of elevations.

VI. CONCLUSIONS
The turbulent wake behind a wall-mounted square cylinder of aspect ratio 4 placed in a thin, developing boundary layer at Reynolds number 12 000 was investigated in terms of vortical structure and shedding process characteristics using S-A IDDES and POD. The quantitative comparisons of the first-and second-order flow statistics with the available experimental and direct numerical simulation data demonstrate a commendable agreement, indicating that the S-A IDDES model and the vortex method of generating the turbulent inflow conditions are suitable for this study. An inspection of the time-averaged results revealed that the horseshoe vortex no longer exists for elevations higher than 0.08D. Instead, two counterrotating vortices are observed behind the obstacle and two small antisymmetric vortices are located near the obstacle sides. It is found that the size of the recirculation bubble decreases toward the free end of the obstacle.
The components of the Reynolds stress show a strong correlation with the change in the vortex structure along the obstacle height. The maximum spanwise velocity fluctuations accrue at an elevation equal to 1D (0.25H), where the shedding process is considered twodimensional. The analysis of the mean streamwise vorticity revealed that it originates from the two tip vortices and another pair of vortices originates from the streamwise-parallel leading corners of the obstacle. Furthermore, the λ 2 visualization showed that the dipole type vortex exists in the mean-field in the downstream direction to the obstacle, while the half-loop hairpin vortex structure is identified from the visualization of the instantaneous flow field. The two-point correlation analysis of the flow near the vortex formation region showed a coherent vortical structure along the obstacle height. Two intervals of fluctuations are obtained from the simultaneous temporal evolution of the surface pressure coefficient on the obstacle sides, i.e., the HAFs, corresponding to the formation of the primary vortex shedding with a phase difference of ∼180 ○ , and the LAFs, with less periodic behavior. The size of the LAFs' interval increases along the ARTICLE scitation.org/journal/adv obstacle height, indicating that the probability of alternating shedding decreases toward the free end of the obstacle. Even though the HAFs' interval has an amplitude similar to that of the LAFs' interval near the free end of the obstacle, the alternating shedding is maintained coherent. The POD analysis of the wake at different elevations revealed that at z/D = 1 to z/D = 2 (0.25H-0.5H), the first two POD modes exhibit a clear periodic behavior and the turbulent kinetic energy contribution is between 66.6% % and 57.6%. For z/D = 4, the first two modes represent the symmetric shedding of the wake and contribute to 36.5% of the total turbulent kinetic energy. The periodic nature of a 1 and a 2 for the elevations where the shedding process is dominated by the alternative shedding, i.e., z/D = 1 to 2, leads to a possibility of developing a simple low-order model based on the vortex-shedding phase angle and the spectrum of the time coefficients obtained from POD. For the elevations higher than z/D = 2, a 1 and a 2 do not show a clear periodic nature, so the simple low-order model is no longer applicable and more computational efforts such as the Galerkin projection or deep learning approach are needed to predict the temporal evolution of the time coefficients.