Use of the POD and the Coherent Structure Detection Criteria to Study the Flow Dynamics Downstream a Confined Square Obstacle

Studying the concept of the vortex is an essential tool in the comprehension of fluid dynamics. Despite this, it is still very difficult to find a universal definition of a vortex. Several methods of detection and characterization of vortex structures has been developed and performed. Specifically tailored for PIV data, they present an important topic in modern experimental fluid mechanics. In fact, the performance of the method is related to its effectiveness in the velocity data analysis and to successfully detect and locate the vortices as well as calculate the characteristic vortex parameters. In the present study, we explore the efficiency of different vortex detection algorithms such as the vorticity  ,  2 and Q criteria by studying the vortices structures in the wake of a confined square obstacle.


INTRODUCTION
The confined flow around bluff bodies has received considerable attention in recent years, due to its fundamental relevance. This type of geometry have a variety of engineering applications in aerodynamics, wind engineering, electronics cooling and flow metering devices, etc. Bluff body cross sections that are often employed are circular and rectangular especially, square. The flow details behind these geometries depend mainly on Reynolds number and blockage ratio. Flow past a square cylinder resembles flow past a circular cylinder as far as instabilities are concerned. However, the vortex formation region is significantly broader and longer for a square cylinder compared to the circular. At low Reynolds number, aspect ratio and end conditions play a significant role in determining the flow properties. The study on the effect of aspect ratio in the literature has been largely focused to a circular cylinder (Rehimi et al. [1], Paranthoën et al. [2]). These experimental works have studied the wake of the cylinder: the effect of the confinement and gap parameter on the different flow regimes and the reorgonisation of the cylinder wake. Sushanta et al. [3] studied the sensitivity of flow patterns in the wake of square cylinder cross section to aspect ratio and orientation with respect to the mean flow (the incidence angle). Numerical simulations are documented in the works of Breuer et al. and Turki et al. [4,6]. These works, especially based on 2D laminar flows simulations, have studied the effect of the Von Karman vortex in the wake behind the obstacle. In this paper, we have restricted our attention to the range of Reynolds were exist a laminar Von Karman vortex street. The appears of this coherent structure is conditioned by the These coherent structures play an essential role in heat and mass transfer. Their interaction and the mecanism of generation/dissipation constitues a fundamental topic in fluid mechanics research. Many definitions of the concept of coherent structures have been proposed (Lugt et al., 1983;Chassaing et al., 2000) [7][8]. However, some definitions are still debated, and the most essential question with coherent structure includes: how do we detected the center of vortex and how do we define its limits? This search of the coherent structure definition has led to the development of the detection criterion. The present study reports experimental measurements of flow patterns in the wake downstream of a square cylinder localazed in a rectangular cross-section and centred in a channel, at low Reynold numbers. The main objective of this work is to study the influence of confinement on different regimes and the reorganisation of the flow in the cylinder wake. For this purpose, in the first step, the vorticity , the Q and  2 criteria have been applied to the PIV velocity fields. It is also essential to add that such configuartion is suitable for the validation of numerical simulations codes. Thus, the same criteria are applied, in the second step to a numerical velocity fields obtained by the Lattice Boltzmann Method simulation [22]. Proper orthogonal decomposition (POD) of the flow was used for filtring and extracting the energetic contribution of different modes. For the experimental investigation, a square obstacle of 10 mm height was placed in a channel with a rectangular cross section 300x30 mm². The confinement ratio was fixed in this study to r=1/3 and the gap to =1.

GENERALIZED LATTICE BOLTZMANN METHOD
Lattice Boltzmann method is a recent numerical method alternative to the classical techniques such as finite-volume, finite-element and finite-difference for solving Navier-Stokes equations.
This method treats the fluid on a statistical level, simulating the movement and interaction of single particles by solving a discrete velocity Boltzmann equation on a particle distribution function noted f(x,c,t).
With the BGK approximation and single relaxation time, this equation can be written as: The discretized form is given by: The right term of Eq. (2) represents the collision process through a relaxation time , while the left one is the streaming process. Many authors studied different approaches to write the collision process and the simplest one in order to obtain a linearized form with a single relaxation time. This form, was introduced by Bhatnagar, Gross and Krook (BGK) [18]. The common notation used for Lattice pattern is refered to the dimension of the problem and the number of speed DnQm where n represents the dimension of the problem and m the speed models.
In 2D problems, there are generally two types of Lattice patterns: square lattice and hexagonal lattice suggested in the literature. Generally, for LBM computation, an equidistant orthogonal Lattice is chosen and on each lattice node, the distribution function f is stored and updated on each new time step. As usual, the 9-velocity square Lattice is the most widely used. In fact, numerical studies proved that both D2Q9 and 7-velocity hexagonal lattice D2Q7 have sufficient Lattice symmetry, which is a dominant requirement for recovery of the correct flow equations [19] and recent studies have shown that D2Q9. It usually gives more accurate results than that based on hexagonal Lattice [20]. Hence, we have used in this study the 9-velocity square lattice D2Q9, on which each particle moves one Lattice unit only along the eight links indicated with 1-8. The 0 indicates the rest particles with zero speed and the discrete velocity vector of particles is defined by: where is the Lattice speed, eq f is the equilibrium distribution functions (EDFs) that plays an essential role in the Lattice Boltzmann Method and recovering the Navier-Stokes equations [18].
This equilibrium distribution function can be considered as a Gaussian quadrature of the Maxwell-Boltzmann continuous distribution function, which is expanded as a Taylor series in macroscopic velocity to its second order to recover incompressible Navier-Stokes equation. According to He et al. [21], it is given by the following form: This local EDF is computed at each time step for every node from the local flow velocity, the fluid density  , the Lattice geometry weighting factor i w and the sound velocity c s .
In this equation, i w represent the weight coefficients and c s is called the sound velocity in LBM, which is related to c and i w .
For D2Q9 model, we have [22]: 6,7,8 and 0 w = 4/9. At each simulation time step, the density  and the velocity field u are respectively defined in terms of the particle distribution function as follows: The kinematic viscosity of the fluid is given by: The great progress of the LB Method during the last decades is due to its many advantages for fluid dynamics simulation.
In fact, only a loop of two steps is implemented during programming: the collision and streaming step which can recover the nonlinear macroscopic advection terms. Also, LBM is a linear equation, where nonlinearity is implicitly imbedded in the collision operator. For incompressible flow, the pressure satisfies a simple state equation unlike the common CFD method, solving Poisson equation becomes unnecessary [12]. The pressure is then defined by the following state equation as:

EXPERIMENTAL SET UP AND MEASUREMENT TECHNIQUES
The experimental device is a hydraulic channel made up of a transparent Plexiglas (figure. 1). The length of the testsection is L=3 m, the width is l=0.3 m and H=0.03 m in height. A square obstacle (of dimensions 0.01x0.01 m²) is placed in the test-section perpendicularly to the mean flow direction at a distance l o =1.6m from the entry. The gap is γ=Δ/D=1, were Δ is the distance between the square and the nearest wall. The confinement ratio r=D/H is fixed at r=1/3. The mean flow velocity flow in the channel section upstream the obstacle was in the range 0.2 to 1.8 cm/s.
In the exit of the test-section, two Brooks flowmeters (range from 0 to 0.530 m 3 /h and from 0 to 1.17 m 3 /h) are assembled in parallel in order to measure the liquid flowrate for a range of Reynolds numbers included between 20 < Re < 400. The liquid is then recovered in a storage tank of 0.5 m 3 of capacity. A pump allows filling a tank of 0.15 m 3 capacity, where the water level is maintained constant by an overflow.
A PIV system has been used to determine instantaneous velocity fields downstream the square. It was composed of a CCD camera of 1600x1184 pixels resolution (Dantec Dynamics Flow sense), a pulsed ND-Yag laser (New wave solo). The whole system is driven by the "DANTEC" software "Flow Manager". The flow seeding was performed using polyamide particles with 20 m of diameter.

. Physical model
A 45° mirror was used to reflect the horizontal laser plane vertically at z = 0. The physical field of view of camera is 81.73x 61 mm². A mask was applied to retain only 30 mm in the y-direction. For the measurements, the interrogation's zones were taken equal to 32x32 pixels. An adaptive cross correlation was used with an overlap of 50%. The sampling frequency was fe=15Hz for a duration of 20 to 40s.

COHERENT STRUCTURES DETECTION
Many criteria used in practice are mainly deduced from the velocity gradient tensor. The objective of these criteria is to characterise the full coherent structure.
The vorticity tensor, as an invariant Galilean, appears as the most widely choice to identify the coherent structures in the flow. Hussain et al. [9] proposed that the coherent structure is a turbulent zone in which, instantaneous vorticity is correlated with the phase of the flow. The vector   is used to compute the rotation rate in the flow. It is defined as Strawn et al. [10] defined the center of the vortex as a local maximum of the vorticity module written for a 2D flow as follow: Jeong and Hussain [11] studied the use of vorticity  and show that the shear and the rotation of the vortex are not distinguished by this criterion. In our case, the presence of shear at the upper and lower walls of the channel deforms the vorticity pattern of the coherent structures. Based an this observation, another criteria is proposed. Known under the name of the Q criterion, it was proposed by Hunt et al. [12]. It can be interpretd as a balance between the rotation rate and the strain rate and can not be affected by the local shear. For a 2D flow, this Q creterion is defined as follow: The vortex structures are identified by the positive iso-values of Q, while their centres are identified by the maximum values of this criterion.
Graftieaux et al. [13] developed a kinematic criterion noted  1 for the vortex center identification considering only the topology of the flow. Although it is not a Galilean invariant, this criterion has presented a simple and robust method for the vortex structures identification. Thereafter this criteria was improved to became a Galilean invariant. This new criterion named  2 has the advantage to consider the local vortex convection. The extremum is used for the vortex identification. This criterion gives accurate results in comparison with the vorticity magnitude or the Q criterion and it is defined as: , e z is a unit vector perpendicular to the velocity map, P is the point where the criterion is calculated and S is the surface arrounding the point P. S should be as small as possible since for a large S, a small scale is filtred. For our case the size of S is: S=8.x  y.

PROPER ORTHOGONAL DECOMPOSITON OF THE FLOW (POD)
The velocity fields obtained from PIV measurements are generally noised. By applying vorticity  or velocity gradient tensor based on the velocity derivatives, high noise frequencies are amplified.
Using filter is not adequate in many cases. For these reasons, the POD applied to the flow velocity fields seems as a good solution. It takes into consideration the energetic distribution of the flow for the acquisition duration.
Proper Orthogonal Decomposition (POD) brings new opportunities because this post-processing extracts spatial or temporal structures using a rigorous mathematical decomposition basis, which breaks free from all arbitrary steps. It became popular for aerodynamics research in the 2000s, starting with Tang et al. [14] works, although it was first proposed in 1960 by Lumley [15] in the context of coherent structure's investigations.
The basic purpose of POD is to extract the most probabilistic structures in an energetic point of view from a statistical series of signals. Its principle is the creation of a mathematical model that decouples the spatial from the temporal variations in the flow.
This modal base needs only few modes to describe all principal aspects of the signal (majority of the large scale behaviour).
Sirovich [16] has introduced the method of snapshots as a way to reduce the computational requirements of the POD calculations, especially when the spatial data size is much larger than the number of images like PIV measurements. For these reasons, we used snapshot in this work.
Considering the measurements obtained by PIV, it yields, generally, N snapshots of a 2D section of the flow field taken at the sequence time t 1 ,…,t N . These snapshots will usually give feature information on the velocity vectors [U(x,y,t 1 ),…, U(x,y,t N )]. A 2D flow field described by this velocity can thus be expressed as: where the temporal coefficients ) ( ) ( t a n and the spatial modes To filter the intantaneous PIV velocity fields, the equation (13) is truncated at a reduced order model M << N, which is chosen as a compromise between model simplicity and model accuracy.
The filterd velocity is written as followed:

Flow modes behind the obstacle
An example of a velocity profile obtained by 2D PIV measurements upstream of the square obstacle is shown in figure 3. In order to verify that the studied flow is developed and symmetrical in the channel, we compare the obtained profile with the theoretically solution of Lundgren and Sparrow [17]. The measurements are in good agreement with the theoretical results.
To characterize the flow and the various coherent structures present in the flow, we used the PIV measurements obtained by varying the Reynolds number. In fact, for a fixed gap and confinment ratio, only the Reynolds number controls the flow regime. In our case, it is defined as where m U is the average axial velocity in the inlet channel flow and  the kinematic viscosity of the flow. The regimes downstream of the square obstacle, similar to those of an obstacle placed in an infinite medium appears by varying the Reynolds number. These schemes are also comparable to those obtained numerically by various authors cited in the literature.  In the second part, the velocity profiles were compared with numerical ones obtained by the Lattice boltzmann Method simulation. An example of experimental velocity profile, superposed to the numerical one obtained by LBM simulation, is illustrated by the figure 5.
For the Reynolds number Re=80 (above the critical regime), and for a longitudinal position close to the obstacle x/D=2, the axial and radial velocity profiles are in accordance with the numerical results. The velocity profile of U x according to the y/D position, presents two peaks on the higher and the lower part of the obstacle (y=±D), and then converges towards zeros on the two wall sides. These peaks show well the presence of two vortices on the same level of the x value.

Coherent structure detection Criteria
In this section, different criteria were applied at first to the PIV velocity fields, then to the numerical velocity fields obtained by the LBM simulation. For this reason, we fixed a specific Reynolds number Re=90 for a periodic regime, which is characterized by an instable wake and the appearance of the von Karman vortices. The streamlines and the axial velocity magnitude contours ( figure 6) show clearly the oscillatory form of the wake. However, the von Karman vortices are not well highlighted as they are embedded in the mean flow. The vorticity fields brings up two families vortex well detectable (figure 7). Using the vorticity  as a creterion for the structure's detection, is not effective in the presence of shearing in the channel near the wall. Hence, we adopted Q and Γ 2 criteria in order to correctly identify the vortices. Figure 8 shows the Q criterion applied to the PIV velocity fields, at first, then to the numerical (LBM) velocity fields. This criterion brings up the structures that we can define their limits by the nearest values around the vortex centre, which is represented by the maximum of Q values. The comparaison between the two results shows a good agreement. The Γ 2 criterion applied to PIV and numerical velocity fields (synchronized at the same time), shows the presence of two types of vortices (figure 9). A first family (P 1 , P 2 ) mainly composed of two counter-rotating vortices similar to those of von Kármán and alternatively detached behind the obstacle. Containment, created by the channel walls, leads to a significant change in the vortex trajectory as well as in the structure of the wake.
For example, the passage of the detached vortex on the upper side (P 1 ), creates a change in the space properties near the top wall. This latter, opposes it by creating a vortex with inverse direction, (P 1 ') which is ejected by forcing (P 1 ) through the channel along the y-axis towards the bottom wall. The same phenomenon occurs for detached vortex (P 2 ) on the lower side and the vortex (P 2 ') is ejected towards the higher wall. This proves that a vortex tends to roll and not slide on a wall.

Figure 8. Q criterion for Re=90: (a) PIV velocity field; (b) LBM simulation velocity field
It should be noted that the distribution of these vortices is different from the infinite medium case. In fact, the vortex (P 1 ) detached from the upper wall does not reach directly the low wall of the channel, but it encounters, first, a vortex (P 2 ), which creates a coalescence between the vortices (at x=6D for this case). 8 It is to be noted also that the confinement reduces the expansion of the vortices in the transverse direction of the flow.
An example of  2 profile at a specific zone of the channel obtained from the PIV measurements, is shown in figure 10. It is in good accordance with result obtained when the  2 criterion is applied to velocity fields obtained by the LBM simulation.

POD DECOMPOSITION OF THE FLOW
The POD was applied to the experimental PIV measurements obtained for a number of snapshots N=300. It was verified that number of snapshots N=300, is largely sufficient for the statistcal convergence of the data for a given Reynolds number. The four first modes are directly related to the von Karman structures and are the most energetic ( figure  11). The reconstruction of filtered data on the basis of this modes permits to recuperate more than 90% of the fluctuating flow energy as % 90 It should be noted that the POD was applied to two positions separetly. For this reason, the issued spatial modes are interpreted independently. The spatial modes ) (n  are paired for Re=90 and present a spatial periodicity ( figure 12). The right columns of figure 12 show the spatial modes for the measurements taken in the position 2 far from of the obstacle. The spatial modes relative to the zone near the obstacle (position 1) are shown in the left column of figure 12. The coherent mode can be clearly observed mainly in the position 1 where the von Karman vortex street is more intense.

CONCLUSIONS
In this work, we presented a study of the wake behind a square obstacle in a confined environment. Experimental results of velocity fields issued from PIV measurements were used and compared with numerical ones obtained by Lattice Boltzmann simulations. The comparison between results shows a good accordance. The presence of the square obstacle in channel creates, from the Re > Re c the von Kármán vortices interacting with the walls to form secondary vortices advected downstream the obstacle. These results are important to show the fundamental difference with von Kármán vortices generated in unconfined medium. Thus, we applied, firstly, the structure criteria detection (, Q and Γ 2 ) to the instantaneous experimental and numerical velocity fields. The objective is to properly extract and describe the dynamics of vortex structures present in the flow, and which play a significant role in transfer and transport phenomena. The obtained results exhibit the interactions between a vortex and a wall, which are characterized by the creation of a secondary one at the wall with opposite vorticity signs. This latter is advected to the downstream then coalesce with the opposite vortices. The decomposition of the velocity field that has been implemented in this research has shown that the relevant structures of the field can be described with a limited number of modes (four) though capturing more than 90% of the kinetic energy of the complete field (mean and fluctuating fields). Subdomain 1 (near the obstacle) has revealed the presence of spatial modes paired where the von Karaman vortex street is mainly intense. The LBM simulation results are validated with our experimental ones which highlights the accuracy of this LBGK method for simulation of flow Dynamics behind bluff bodies.