Proper orthogonal decomposition analysis of coherent motions in a turbulent annular jet

A three-dimensional incompressible annular jet is simulated by the large eddy simulation (LES) method at a Reynolds number Re = 8 500. The time-averaged velocity field shows an asymmetric wake behind the central bluff-body although the flow geometry is symmetric. The proper orthogonal decomposition (POD) analysis of the velocity fluctuation vectors is conducted to study the flow dynamics of the wake flow. The distribution of turbulent kinetic energy across the three-dimensional POD modes shows that the first four eigenmodes each capture more than 1% of the turbulent kinetic energy, and hence their impact on the wake dynamics is studied. The results demonstrate that the asymmetric mean flow in the near-field of the annular jet is related to the first two POD modes which correspond to a radial shift of the stagnation point. The modes 3 and 4 involve the stretching or squeezing effects of the recirculation region in the radial direction. In addition, the spatial structure of these four POD eigenmodes also shows the counter-rotating vortices in the streamwise direction downstream of the flow reversal region.


Introduction
An annular jet, which is regarded as the limiting case of a coaxial jet [1] , is utilized widely in combustion devices [2][3] . Early studies showed that annular jets share the similarity properties (such as expansion, decay, and mean axial velocity profiles) of those found in turbulent round jets in the far-field region [4][5][6] . Apart from the jet characteristics, the flow field of an annular jet also features a wake flow near the nozzle. Due to the existence of a central obstacle, flow separation occurs, and a recirculation region is developed behind it, which is served as a way to enhancing a fuel-air mixing rate and reducing pollutant emission in combustion systems. Also, the central bluff-body can have a function as an injector, from which fuel can be injected axially or radially [7] . In addition, two shear layers originating from the annular orifice exist in the near-field. The outer shear layer is between an ambient medium and a jet, and the inner shear layer is between a recirculation zone and a jet. These two shear layers meet and interact with each other downstream of the reversal flow region [8] .
Although the flow has both a symmetric geometry and symmetric boundary conditions, the investigations on turbulent annular jets have shown that spontaneous break of symmetry of the time-averaged flow field could appear [9][10][11][12][13][14] . This type of asymmetry is observed in laminar annular jets as well, if the Reynolds number is above a certain threshold value, which is around 150 for the geometry used in this study [15] . In turbulent regimes, this asymmetry is illustrated by a low-frequency oscillating vortex shedding from the buff-body that is convected to the end of the recirculation zone, creating a shift of the stagnation point in the radial direction [8,9,11,16] . On the other hand, in the laminar regime, Del Taglia et al. [13] found that the occurrence of this flow asymmetry is influenced by the jet blockage ratio and the Reynolds number, which can be integrated as a state parameter. Similar kinds of asymmetric flow patterns were also observed in other flow configurations, such as in turbulent wakes past blunt bodies [17][18] .
Several investigations have been conducted to address the asymmetric mean wake in turbulent annular jets. Patte-Rouland et al. [8] applied the proper orthogonal decomposition (POD) to the particle image velocimetry (PIV) velocity fields and showed that the asymmetric flow is related to the first two POD modes. It should also be noted that the POD analysis is not applied on the fluctuating velocity part but on the instantaneous fields, which means that the first POD mode denotes the mean flow. However, their PIV velocity fields are two-dimensional, and no temporal dynamics is obtained owing to the relative low sampling frequency. Vanierschot and Van den Bulck [11] investigated the precession of the central toroidal vortex behind the central bluff-body using the phase averaging method, which treats the coherent fluctuations in the flow field as a whole. The symmetry breaking and vortex precession in low-swirling annular jets were studied by Vanierschot et al. [12] with PIV and POD techniques, showing that the motion of the asymmetric wake is random in time. The drawback of their study is that the measured PIV velocity vectors were two-dimensional as well. The first experimental study to reveal the three-dimensional structure of the wake was done by Vanierschot et al. [16] . Using tomographic PIV measurements, they found several asymmetric vortical structures which were shed from the inner shear layer towards the stagnation point.
This study aims to reveal the connection between the asymmetric mean flow field and the turbulent coherent motions in an annular jet. First, the flow is simulated by the large eddy simulation (LES) method. Then, the snapshot POD technique is adopted to capture the most energetic coherent structures in the fluctuating flow. The main concerns here are the spatial structures of the low-frequency coherent motions and their relationship with the flow dynamics, especially on the wake or recirculation region behind the central bluff-body. To the authors' knowledge, a three-dimensional POD analysis of detailed LES data has not been conducted yet. The remainder of this paper is organized as follows. The numerical method is briefly introduced in Section 2. Section 3 presents the asymmetric time-averaged flow field. The POD analysis of the velocity fields is given in Section 4. Finally, we summarize our findings and present the conclusions in Section 5.

Numerical method
The governing equations for the concerned flow are the filtered three-dimensional incompressible Navier-Stokes equations for the Newtonian fluids, Here, u i are the resolved velocities, p is the filtered kinematic pressure, t is time, x i are spatial coordinates, and ν is the kinematic viscosity of the fluid. The subgrid-scale stress tensor τ ij = u i u j − u i u j in the above formulation is modeled by the localized dynamic k model introduced by Kim and Menon [19][20] . Recently, we have used the same numerical approach to simulate the turbulent annular swirling jet flow in a similar configuration to identify the helical vortex cores in the flow field [21] . It was found that the numerical results were in line with experimental measurements, validating the approach. The only difference in this work is that, for a zero swirl intensity annular jet, there is no tangential injection of fluid in the swirl generator, as depicted in Fig. 1(a). The annular jet is generated by a fluid (water with a viscosity ν = 1.005×10 −6 m 2 /s) flowing through a swirl generator with twelve fixed and movable blocks, and then guided into the passage between two coaxial cylinders. The configuration of the annular tube near the jet orifice is shown in Fig. 1   The governing equations are solved by the open-source C++ library OpenFOAM based on the finite volume method. The temporal discretization scheme is second-order implicit backward. For the spatial derivative terms, the convective term is discretized by a second-order upwind scheme, whereas the diffusive term is approximated by a second-order central difference scheme. The computational grid is hexahedrally structured and nonuniform and has a total number of around 10.5 million cells. There are 384 control volumes in the azimuthal direction and 356 cells in the main flow direction from the inlet to the outlet. The finest grid spacing is situated at the jet orifice with ∆x = D o /600 in the axial direction and ∆r = D o /1 200 in the radial direction. The time step is fixed at 8 × 10 −6 s to meet the Courant-Friedrichs-Lewy condition, and the maximum Courant number is less than 0.55 during the simulation. The detailed information about the numerical method and its validation can be found in our previous work [21] . can see that the recirculation zone behind the inner cylinder is asymmetric although the flow configuration has an axisymmetric geometry. It is worth mentioning that the backflow region is precessing gradually around the jet central axis as the averaging time interval increases, which means that it contains very low frequency dynamics. Depending on the shape of the recirculation zone, both an "asymmetric" plane and a "symmetric" plane can be defined, similar to the study of Vanierschot et al. [16] , and both planes are perpendicular to each other, as indicated by the slices in Fig. 2(b). Here, the "asymmetric" plane can be obtained by the cross-section cutting through the axis of symmetry and the stagnation point, whereas the "symmetric" plane is orthogonal to it [10,16] . Figures 3(a) and 3(b) show the corresponding time-averaged axial velocity field u /u 0 and velocity vectors in these two planes, where the radius r is calculated as r = y 2 + z 2 . The stagnation point that indicates the end of the recirculation region is located at x/D o = 0.55. As shown in the figures, the reversal flow region is not parallel to the jet central axis in the asymmetric plane, and it is almost parallel to the jet center axis in the symmetric plane. These results agree well with the findings reported by Ryzhenkov et al. [10] and Vanierschot et al. [16] . Furthermore, the normalized normal Reynolds stress u u /u 2 0 in these two planes is expressed in Figs. 3(c) and 3(d). The figures show that large velocity fluctuations exist in the inner and outer shear layer regions owing to turbulent mixing. The largest normal Reynolds stress is located in the inner shear layer in the asymmetric plane, which is caused by the highly dynamic properties of the reversal region or oscillations of the stagnation point [8] . This phenomenon was also observed in a laminar annular jet at Re = 180 by Ogus et al. [14] .

POD modes
POD looks for spatial base functions that capture the most energetic components of a dynamical system [22] . Suppose that N snapshots of the three-dimensional velocity field u(x, t) are recorded during the numerical simulation, and POD decomposes the fluctuating velocity u (x, t) into a set of spatial modes ψ i (x) with temporal coefficients a i (t),  by separating the independent variables x and t. In this work, the snapshot POD approach of Sirovich was adopted [23][24] . The mode amplitude a i (t) is obtained from the eigenvectors of the auto-correlation matrix Here, the eigenvalues λ i stand for the energy contained by each spatial mode. The eigenfunctions ψ i (x) are calculated by projection of the snapshots onto the temporal coefficients a i (t), The asymmetric recirculation zone discussed in Section 3 originates from random perturbations which are shedding from the central cylinder and are convected downstream [8][9] . Although the asymmetric time-averaged statistics of the velocity fields are observed, an appropriate choice of the time averaging interval could provide an axisymmetric flow field. Figures 4(a) and 4(b) give the averaged streamwise velocity in an axial plane and in a cross-section (at x/D o = 0.4), which is calculated from 500 instantaneous flow fields (N = 500). The time interval between two successive snapshots is 0.024 s, corresponding to a total sampling time of 12 s. It is seen that the axial velocity field is almost symmetric about the jet central axis in Figs. 4(a) and 4(b), i.e., the averaged flow field contains an axisymmetric recirculation region. Therefore, these 500 snapshots are used in the POD analysis. Figures 5(a) and 5(b) respectively provide the distributions of turbulent kinetic energy across the three-dimensional POD modes and the accumulated energy distributions. In Fig. 5(a), we can see that the first four POD modes all have a turbulent kinetic energy content larger than 1%. In detail, the modes 1 and 2 account for 14.7% and 10.4%, respectively, while the modes 3 and 4 capture 2.1% and 1.3% of the total turbulent kinetic energy. In contrast, the remaining POD modes are weaker as each of them represents less than 1% of the total turbulent kinetic energy. If the Reynolds number is increased, the first four modes would capture less kinetic energy. The reason is that more small-scale turbulent structures appear in the flow field with increasing the Reynolds number, and these fine scale structures also have a contribution to the energy content, which means that the relative contribution of the large scale structures decreases. The dash-dotted line in Fig. 5(b) indicates that the first 82 POD modes reach more than 50% of the total energy. We will illustrate the first four eigenmodes in the following contents, as no noticeable spatial structures can be found in the other ones. The first POD modes in an axial plane and a cross-section (at x/D o = 0.6) are plotted in Figs. 6(a) and 6(b), respectively. In the axial plane, two vortices can be identified between the inner and outer mixing layers. In addition, a strong vortex is located immediately behind the central bluff-body, which leads to a radial movement of the fluid in the recirculation region. This effect is responsible for the shift of the stagnation point as well. The results discussed here are in line with the two-dimensional PIV measurements by Patte-Rouland et al. [8] and by Vanierschot et al. [12,25] . Moreover, the velocity vectors of the first spatial mode display a pair of counter-rotating vortices in the streamwise direction, as shown in Fig. 6(b). The structure of the second POD mode is very similar to the first one, with an angular difference of π/2 in the azimuthal direction, as presented in Figs. 6(c) and 6(d). As a consequence, the two-dimensional POD analyses conducted in previous studies [8,12,25] cannot depict these two modes accurately. The amplitudes of the first two spatial modes as a function of time are plotted in Fig. 7(a). It can be seen that, although the temporal coefficients exhibit small scale fluctuations in time, the curves show a general periodic trend. The normalized cross-correlation coefficient c 12 of a 1 and a 2 is given in Fig. 7(b), where c ij (∆t) = a i (t + ∆t)a j (t) . The main features of the correlation coefficient are the maxima and minima of the curves, showing the periodicity of the motion of the first two modes in the near-field of the annular jet. The positive and negative peaks correspond to a period of around 12 s, demonstrating that the coherent motions are in the lowfrequency regime. From the results above, it can be concluded that the first two POD modes have a similar spatial structure and amplitude and hence form a mode pair. The motion of this pair accounts for 25.1% of the total turbulent kinetic energy. It should be noted that the period obtained from the cross-correlation coefficient may only describe the periodic motion roughly, since the POD analysis only includes sampled data during about one period. More numerical data however would need much more simulation time. To study the dynamical structures in the flow of smaller frequencies, a wavelet analysis of the mode amplitudes a 1 and a 2 is done.
The results are plotted in Figs. 7(c) and 7(d), respectively. The figures show that the temporal coefficients have low-frequency fluctuating components from 0.28 Hz to 8.81 Hz, corresponding to a Strouhal number St = f D h /u 0 between 2.7×10 −3 and 0.083. This dimensionless frequency range is much lower than the one associated with inner shear layer instabilities, which is in the order of 0.15-0.3 [6,16] . The largest magnitude of the velocity is located in the area close to the end of the reversal flow region. As discussed later, this mode represents the stretching and squeezing effects in the radial direction of the reversal flow in its end region. In other words, it depicts the opening and closing of the annular jet. Furthermore, two pairs of counter-rotating vortices result from this mode in the streamwise direction. Similar to the first mode pair, the fourth POD mode has a similar structure as the third one, with an angle difference of π/4 in the azimuthal direction, as presented in Figs. 8(c) and 8(d). Figure 9(a) plots the time evolution of the amplitudes of the modes 3 and 4, and their cross-correlation coefficient c 34 is given in Fig. 9(b). The positive and negative peaks in the curve show that this mode pair has a period of about 6 s, i.e., its period is half the precession period of the first mode pair, which is in line with the result obtained by Vanierschot et al. [25] . This mode pair represents 3.4% of the total turbulent kinetic energy. To reveal the relationship between the POD mode pairs and the asymmetric flow reversal, the velocity field is reconstructed with selected spatial modes as Figures 10(a)-10(d) display the positive and negative iso-surfaces of ψ k (x) for the first four POD modes, respectively. The axial component of those eigenmodes presents a pulsatile motion along the streamwise direction. Also, it verifies the harmonic relation between the spatial modes in the first and second mode pairs. The recirculation regions in the wake flow reconstructed by the first mode pair with maximum and minimum temporal coefficients are given in Figs. 10(e) and 10(f), respectively. The shape of the reversal flow zone in Figs. 10(e) and 10(f) is similar to the one shown in Fig. 2. The velocity fields reconstructed by the first mode pair in the asymmetric and symmetric planes, defined in Fig. 2(b), are illustrated in Figs. 11(a) and 11(b), respectively. It shows two vortices that have different sizes in the asymmetric plane, and two vortices that almost have the same size in the symmetric plane. The flow patterns shown in defined in Fig. 2(b), i.e., it has two symmetric planes. The velocity fields reconstructed by the modes 3 and 4 in these two slices are respectively given in Figs. 12(a) and 12(b). The velocity vectors in both slices show two vortices between the backflow region and the inner shear layer. However, the vortices in Fig. 12(a) are smaller compared with those in Fig. 12(b). The reason is that the recirculation zone is stretched in the radial direction in Fig. 12(a), whereas it is squeezed in the radial direction in Fig. 12 flow fields reconstructed by the modes 3 and 4 display two pairs of counter-rotating vortices in the main flow direction, but they are weaker as the second mode pair captures lower turbulent kinetic energy. The counter-rotating streamwise vortices were detected by the phase averaging method applied to the wake in an annular jet [11] , and they were also observed in a laminar annular jet flow [14] .

Conclusions
In summary, we present the flow dynamics of the wake in a turbulent annular jet at a Reynolds number Re = 8 500. It is found that the time-averaged velocity field is asymmetric although the flow has a symmetric geometry and boundary conditions. The POD analysis of the fluctuating velocity field shows that the influence of the first four eigenmodes on the mean flow cannot be neglected. The first two spatial modes account for 25.1% of the total turbulent kinetic energy, while the modes 3 and 4 only represent 3.4% of the total energy. The first POD mode pair leads to the radial movement of the fluid in the recirculation region. It corresponds to the shift of stagnation point in the radial direction and is responsible for the asymmetric timeaveraged reversal flow region. In contrast, the second mode pair has stretching and squeezing effects on the recirculation region in the radial direction. The flow fields reconstructed by these four spatial modes demonstrate that the recirculation zone precesses around the jet central line at a very low frequency, corresponding to a Strouhal number St less than 0.001. Moreover, the flow fields reconstructed by those two POD mode pairs also present counter-rotating vortices in the streamwise direction downstream of the reversal flow region.