Vortex Core Dynamics in a Swirling Jet Near Vortex Breakdown

The dynamics of the vortex core are investigated in a swirling jet at swirl numbers in the range of the critical swirl number for vortex breakdown. Vortex breakdown, a bifurcation in the structure of a swirling jet, results in the establishment of a stagnation point and recirculation region along the centerline of the jet. In this study, we explore the dynamics of the swirling jet near vortex breakdown. Investigation of time-averaged velocity fields and profiles leads to the identification of three flow regimes: pre-breakdown, near-breakdown, and post-breakdown. Velocity fields in these regimes are further analyzed using dynamic mode decomposition, Rankine-vortex fitting, and proper orthogonal decomposition to characterize jet dynamics with a particular focus on the development of the recirculation region characteristic of vortex breakdown. A precessing vortex core is also identified in the post-breakdown regime and its behavior is discussed.


Nomenclature
= an infinite-dimensional linear operator = initial amplitude of the kth dunamic mode D = nozzle diameter DMD = dynamic mode decomposition = the kth unit vector pdf = probability distribution function PIV = particle image velocimetry POD = proper orthogonal decomposition R = vortex core radius = residual vector r = radial coordinate = companion matrix with eigenvalues and eigenvectors that approximate those of S = geometric swirl number t = time u = velocity VB = vortex breakdown = matrix of velocity snapshots p through q x = axial coordinate

II. Experimental Configuration
The experimental facility used in this study is illustrated in Figure 1. The facility is separated into three main components: a two-microphone injector nozzle, a variable-angle, radial-entry swirler, and a settling chamber fitted with speakers that provide longitudinal forcing. Air enters the 15.25 cm diameter settling chamber at the base of the experiment, where it passes through two perforated plates to reduce the impact of incoming turbulence. Flow proceeds upwards into the swirler chamber, where it is turned inwards towards the radial-entry swirler. The swirler is made up of eight NACA 0025 airfoils, each with a height and chord length of 2.54 cm. To allow for variations in swirl, a stepper motor at the base of the experiment controls the angle of these airfoils. In this configuration, airfoil blade angle can be varied between 65° and -65° with a resolution of 2.5°. A centerbody in the swirler directs flow upwards into the injector nozzle and tapers to a point 8.9 cm upstream of the nozzle exit. In the nozzle itself, flow passes two pressure transducers (PCB Model 113B28), mounted 6.92 cm and 1.84 cm upstream of American Institute of Aeronautics and Astronautics the exit. These positions correspond to 2.54 cm and 7.62 cm downstream of the swirler, respectively. Air emerges from the 2.54 cm diameter nozzle as an unconfined jet. Throughout this study, the volumetric flow rate is maintained at 30 SCFM with a maximum deviation of 0.5 SCFM. This results in a bulk flow velocity of 28 m/s and a Reynolds number of 35,000 based on the nozzle diameter. Flow rate is monitored using a thermal mass flow meter (Thermal Instruments model 600-9). Signals from the flow meter and pressure transducers are sampled at 20 kHz for 3 seconds using a National Instruments CompactRIO system.
Stereo particle image velocimetry (PIV) is used to measure three components of the velocity field in the r-x plane, as indicated by the axes shown in Figure 1. In this 'r-x configuration,' radial and axial velocities are in-plane, while azimuthal velocity is out-of-plane. PIV images are captured by two Photron SA5 CMOS high-speed cameras in double frame mode. The cameras are arranged in a forward-forward configuration and are equipped with Scheimpflug adapters. Images are taken at 5 kHz with a variable interframe time of 17 -23 µs depending on swirl number. A Hawk/Darwin Duo Nd-YAG, 532 nm, 60 W laser is used to illuminate the motion of aluminum oxide tracer particles. These particles have a nominal diameter of 1 -2 µm and can track flow perturbatons up to 4000 Hz.
Velocity fields are calculated in DaVis 8.3.1. Multi-pass cross-correlation is conducted without pre-processing or masking. A 32x32 pixel interrogation window with 50% overlap is used in the first pass and a 16x16 pixel window with 50% overlap is used in the two subsequent passes. Two methods are used to reject spurious vectors during postprocessing. First, a vector is removed and replaced if its magnitude is greater than 3 times the RMS of the surrounding vectors. Second, DaVis employs universal outlier detection to identify and replace outlying vectors.
Six levels of swirl are examined in this study, corresponding to blade angles of 0°, 15°, 30°, 35°, 40° and 45°. Swirl is quantified based on geometric swirl number, S, as defined by Lefebvre; 43 the relationship between angle and swirl number is listed in Table 1.

A. Dynamic Mode Decomposition
The identification of coherent structures was conducted using dynamic mode decomposition (DMD), as introduced to the study of flow dynamics by Schmid. 44 DMD allows for the decomposition of a velocity field time series into spatial basis functions, , each oscillating at a unique complex frequency, . At its core, DMD relies on the assumption that there exists an operator, , that linearly maps each velocity snapshot, , to each subsequent snapshot, +1 , such that For the purposes of DMD, sequential velocity snapshots, separated by a time step Δ , are represented as columns in a data matrix, 1 , where American Institute of Aeronautics and Astronautics While the eigenvalues and corresponding eigenvectors of can be computed at this stage, the application of DMD to experimental data requires additional care due to the presence of noise and turbulence. Thus, a pre-processing step is necessary in order to eliminate low-energy flow behavior that could obscure any coherent structures that may be extracted via DMD. To accomplish this, proper orthogonal decomposition (POD) 45 , a decomposition method commonly used to extract energetic coherent structures from flow fields, is used as a filter prior to the application of DMD. This filtering technique was derived from the higher-order DMD method of Le Clainche and Vega. 46 Prior to the approximation of , POD is applied to the data matrix 1 via the singular value decomposition (SVD), such that 1 = * . The matrices and contain the spatial and temporal content, respectively, of the original data matrix. The diagonal matrix contains the singular values, , of the data matrix. The square of these singular values is a measure of the fluctuating energy of the corresponding POD mode. To eliminate modes that contain noise and turbulence, the SVD is truncated, retaining modes according to the following relation: where is the fraction of total turbulent kinetic energy that is contained in the removed modes. The selection of this truncation parameter is nontrivial, and is discussed in further detail in the following section. Using the results of the truncated SVD, a "reduced" snapshot matrix is found as ̂1 =̂̂ * , and it is this matrix that is used to find the DMD modes. In order to approximate the companion matrix in Equation 4, new data matrices ̂1 −1 and ̂2 are constructed. SVD is applied to the former, yielding ̂1 −1 = * . Thus, ̃, a dimension-reduced approximation of , is be found as ̃= * ̃2 −1 . The eigenvalues and eigenvectors of are approximated by the eigenvalues and eigenvectors of ̃, which are found by solving the eigenvalue problem below: To construct the DMD modes, , the eigenvectors found above are projected back into the space of the original data matrix via Equation 7 and the corresponding eigenvalues are transformed into complex frequencies via Equation 8: The real parts of these complex frequencies are the temporal growth rates of each mode, and the imaginary parts are the corresponding frequencies. The initial amplitudes of the DMD modes, , is calculated: Where is the matrix containing all of the DMD modes, , and the superscript '+' denotes the Moore-Penrose pseudoinverse. Thus, the original time series data can be approximately reconstructed as a linear combination of the resulting DMD modes:

B. Truncation Parameter Selection
The truncation parameter, , as introduced in Eq. 5, must be carefully selected in order to maintain the integrity of the final set of modes. Truncating too much could remove relevant physical motions, but truncating too little results in poor performance of the DMD algorithm and nonphysical results. According to Le Clainche and Vega, 46 optimal DMD performance is achieved when the resulting number of modes, , is equal to the dimension of the subspace spanned by the modes, . The values and are referred to as the spectral complexity and spatial complexity, respectively. Without POD filtering, attempts to apply DMD to the raw velocity snapshot data obtained from PIV resulted in a complexity mismatch where < ; this situation also occurs for small values of the truncation parameter American Institute of Aeronautics and Astronautics ( < 1%). Above a threshold of approximately 1%, which varied slightly for each swirl number considered, the spectral and spatial complexity of the resulting DMD modes was no longer mismatched. Thus, the lower limit for was set by the value at which the spectral and spatial complexities are approximately matched, i.e., = .
To determine a set threshold for each swirl number considered, a balance was struck between practical analytical considerations and energy retention. Generally, for the same value of , jets with higher swirl number require more modes to capture the same fraction of fluctuating energy. Thus, a single threshold could not be selected for all cases. Instead, a value for the truncation parameter was selected for each of the three flow regimes identified in Section IV. This ensures that the two cases within each regime maintain approximately the same number of modes and that motions that may be present in one case are not filtered out completely in another case in the same flow regime. The values selected for are listed in Table 2.

A. Time-Averaged Flow Structure
The six swirl numbers examined in this study were selected in order to understand the behavior of vortex breakdown as it forms and eventually dominates the structure of the swirling jet. Generally, the non-swirling jet (S = 0) is used as a baseline with which to compare the behavior of the swirling jets. As swirl is added to the jet, vortex breakdown begins to manifest as a region of diminished axial velocity in the center of the jet. As can be seen in the time-averaged axial flow fields presented in Figure 2, this region grows in size and travels upstream as the level of swirl is increased. At intermediate levels of swirl, i.e., S = 0.38 and 0.46, a region of stagnation intermittently forms in the flow field. In Figure 2, this appears as a time-averaged velocity deficit along the jet centerline that increases with swirl number. At the highest swirl number considered in this study, S = 0.66, time-averaged stagnation is visible and vortex breakdown fully and consistently dominates the flow field along the centerline. The axial velocity deficit that forms with increasing swirl number is better visualized in the time-averaged axial velocity profiles in Figure 3. At S = 0, the axial velocity profile is relatively flat, peaking at the jet centerline. As the level of swirl increases, the profile becomes bimodal, peaking on either side of the centerline. The centerline velicity deficit continues to increase as swirl increases. Notably, this deficit increases significantly between S = 0.46 and 0.56, suggesting the emergence of a strong stagnation region consistent with vortex breakdown. At S = 0.66, the average axial velocity along the centerline is near zero, indicating the onset of time-averaged stagnation and the consistent presence of vortex breakdown in the jet. American Institute of Aeronautics and Astronautics

B. Pre-Breakdown Analysis
The temporal growth rates provided by DMD allow for the classification of modes as stable or unstable. To characterize the results of DMD, mode spectra such as those found in Figure 5 are generated. These spectra allow for the identification of particularly unstable modes that may lend further insight into the flow processes taking place. Notably, the growth rates of all modes considered are negative. This is likely due to averaging effects, as the timescales of the modal dynamics are much shorter than the PIV time series from which the modes are extracted. While no quantitative results can be drawn from the resulting growth rates, modes can still be identified as more or less stable based on their relative position in the spectrum. Modes that seem to be local outliers, especially those with larger growth rates, are of particular interest, as it is these modes that dominate the flow field. One such local outlier is the mean flow field, which appears in every mode specrtrum with a frequency of 0 Hz and a growth rate near zero. Unlike American Institute of Aeronautics and Astronautics with POD, the mean cannot be subtracted from each velocity snapshot prior to applying DMD. While the same mode shapes and frequencies will be identified if the mean is removed, all corresponding growth rates will be zero and stability information will be lost.
Another criterion used to evaluate the DMD modes is a measure of coherence calculated by projecting each dynamic mode onto the POD basis found during the POD filtering process. This coherence value is shown as the coloring of the points plotted in the mode spectra. Modes with high coherence (in red) contain large-scale, energetic structures while those with low coherence (in blue) contain small-scale structures that are less energetic. As Schmid 44 notes, however, modes dynamic modes with a small coherence of the projection onto the POD basis can still represent dynamically important motions.
The mode spectra for the swirl numbers analyzed in the pre-breakdown regime are shown in Figure 5. The notable modes identified for S = 0 largely contain axisymmetric behavior similar to that shown in Figure 6a, as is characteristic of non-swirling jets. 7 Vortex rings formed by the axial Kelvin-Helmholtz instability are found at a wide range of frequencies, notably at 734 Hz, as shown in Figure 6b. The pairing of these vortices as they convect away from the nozzle results in the larger-scale motions found in modes with lower frequency (not shown). This behavior is also seen at S = 0.18, but is complicated by the presence of swirl. As swirl is added to the jet, these vortex rings tilt, particularly farther downstream. This can be seen in Figure 7.

C. Near-Breakdown Modal Analysis
As swirl is increased into the near-breakdown regime, the jet approaches the onset of vortex breakdown. This is characterized by an increase in intermittent behavior. To understand the general nature of the fluctuations present in the flow field, outlier events in axial velocity are considered in Figure 8 for S = 0.38. Figure 8   These fluctuations are further characterized in this regime using DMD. In particular, the emergence of the central velocity-deficit region characteristic of vortex breakdown is evidenced by the increased modal activity along the jet centerline across a range of frequencies at S = 0.38. Figure 11 illustrates the emergence of the central shear layer formed due to vortex breakdown. Upstream of x/D = 3, the 560 Hz mode features axisymmetric motion similar to that seen in Figure 6b in the pre-breakdown regime. Downstream of x/D = 3, activity shifts to the centerline of the jet, suggesting the presence of a central shear layer in which vortical structures can form. While many modes depict asymmetric activity along the centerline similar to that in Figure 11, no modes contain clear axisymmetric or antisymmetric motion. This suggests that vortex breakdown has not fully set in and that a recirculation zone has not stabilized along the centerline of the flow. At S = 0.46, motions along the centerline are clearer and better defined, as in Figure 10, which depicts axisymmetric motion along the centerline taking place at a frequency of 265 Hz. To better illustrate the distinction between modal activity in the outer shear layer and activity in the inner shear layer, the timeaveraged locations of maximum/minimum vorticity are overlayed on the mode in Figure 10. These locations correspond to the shear layers, with the leftmost and rightmost lines indicating the outer shear layer and the central lines indicating the inner shear layer. Both axisymmetric motions, as in Figure 10, and antisymmetric motions (not shown) appear along the centerline, suggest the onset of a more stable and consistent recirculation zone as compared to S = 0.46.

D. Near-Breakdown Vortex Core Dynamics
Independently of the DMD decomposition, we extracted information from the azimuthal and radial velocities related to the dynamics of the vortex core. The starting point is the observation of a significant mean radial velocity along the jet's axis, as seen in Figure 12. We interpreted these velocities as the result of the vortex core being outside the plane of the PIV measurements. This hypothesis leads to the azimuthal component (swirl) appearing along the ydirection in the PIV plane. We use this hypothesis to identify a vortex core location. A vortex with its axis parallel to the plane of measurements will have two components at any given point: along the r axis ("radial") and along the z axis ("azimuthal"). For an ideal steady line vortex, the velocity would be tangential around the vortex axis, with magnitude varying with distance from the axis. Using a Rankine vortex model, i.e., a rigid-body vortex core surrounded by inviscid induced motion, two velocity vectors identify the vortex axis location at the intersection of the lines normal to the velocity. Because of measurement uncertainties and ambient turbulence, we calculated a best fit for axis location, using two coordinates, r and z, and rotation speed based on vectors at small |r|. In a first attempt, the fit was calculated based on the seven points straddling the r = 0 line; an improved algorithm, keeping the best fit from a set of sliding four-point-windows near r = 0, yielded the results described below. Vortex core reconstruction was calculated at each x/D < 3 for successive frames. The r-x-z-t location of the vortex core was calculated, as well as the vortex strength using a rigid body model. Two representative views, at a particular instant, are shown in Figure 13 and Figure 14. The top view in Figure 13, with the axial location color-coded from blue to red, shows that the vortex core wobbles in a region centered at about z/D = .05, r/D = 0, i.e., just in front of the PIV axis. In the 3D projection view shown in Figure 14, the vortex strength is color-coded from green (high roatation rate) to blue (low rotation rate); good consistency of the vortex strength along the core seems to validate the algorithm, the model and the convergence.
Large excursions from the vortex center are observed mostly near the nozzle and at x/D > 2; these appear physically realistic, but the quality of the algorithm is likely to be less accurate at large distances from the centerline.  Lack of convergence at isolated x/D levels, resulting in very large excursions from the frame of view, was skipped in the plots. It is unclear if the vortex wobble is due to buffeting by the turbulence, by vortex instability, by incomplete convergence, or a combination of such effects. Convergence of the optimization procedure at each x/D is excellent for S = 0.38 and 0.46, and deteriorates significantly at larger swirl numbers. Poor convergence was diagnosed as large values of the penalty function and outlier excursions from the geometric center of the jet, and the corresponding level (x/D) was deleted from the vortex geometry. Fewer than three such points were noted for S = 0.38 and S = 0.46. Probability distributions (pdfs) of the vortex core location were calculated for x/D < 2.5 over all 5000 data frames. Two representative pdfs are shown in Figure 15 for S = 0.38; results for S = 0.46 are not very different. We see that the core location is slighlty off the PIV plane, and drifts by a few millimeters at successive axial locations. The width of the distribution increases as x/D increases. This observation was made quantitative by calculating the first and second moments of the pdf, with the results shown in Figure 16. The first moments determine the mean r and z coordinates of the vortex core at each level x/D, whereas the second moment quantifies the spread of the distribution around the central point. Except for the immediate vicinity of the nozzle, the trends are very smooth. Furthermore, the trends for S = 0.46 and S = 0.56 are very similar, which gives confidence in the statistical reliability of the approach. To validate the physical relevance of this vortex wobble, a point was selected on the vortex near x/D = 0.4 and its z-coordinate was calculated as a function of time using the Rankine vortex model. This signal was then crosscorrelated with the radial velocity in the entire field via a process termed space-time cross-correlation. Figure 17 shows a sequence of snapshots of the space-time cross-correlation. The reference signal was taken as the zc coordinate (out of the PIV plane) of the vortex core location close to the dump plane. The radial component of velocity is crosscorrelated with zc at every (r,x) point. The color scale is common to all lags. The sequence shows the propagation of fluctuations in two ways. Most obvious is the region of saturated color along the centerline near x/D = 0, indicating that the radial velocity is very strongly correlated with out-of-plane core location for short negative and positive lags, establishing the consistency of our analysis. We see that the vortex core is slanted relative to the axis r = 0. We also see weakly-correlated but slowly and consistently propagating red and grey regions along |r/D| ≈ 1 and x/D ≈ 3 with

E. Post-Breakdown Modal Analysis
As swirl is increased into the post-breakdown regime, the velocity fluctuation distributions illustrate a significant structural change. These distributions are shown in Figure 18 and use the same magnitude thresholds as those for S = 0.38 shown in Figure 8. The emergence of a stable central shear layer, as visualized in Figure 4, is evident in the appearance of strong axial fluctuations along the centerline in Figure 18a. However, the central shear layer has little impact on the spatial distribution of low-freuqency fluctuations, seen in Figure 18b. These fluctuations remain concentrated in the outer shear layer, though this shear layer does exhibit additional vortex spreading as compared to Figure 8b. Most notably, the jet contains significantly more high-frequency activity, shown in Figure 18b. This activity is concentrated farther upstream than at lower swirl, extending all the way to the nozzle. In addition, radial velocity fluctuations are more present than at lower swirl, suggesting a general increase in flow field fluctuation. Similar trends were found for S = 0.66. At S = 0.56, both antisymmetric and axisymmetric motions appear in both shear layers. Notably, a highly damped, antisymmetric, short length-scale motion appears at a frequency of 884 Hz. This mode manifests as a local outlier in the frequency spectrum in Figure 19a and is shown in Figure 20. A similarly-structured mode appears at a frequency just below it, at 842 Hz (not shown). This motion is characteristic of a precessing vortex core (PVC), as has been shown in other works. 38,47 The highly damped nature of these modes suggests that the PVC, which is a global instability, has not fully developed in the flow. At S = 0.66, a clearer and more coherent PVC is present in the flow at a frequency of 858 Hz. This mode is a clear local outlier in the mode spectrum in Figure 19b and is shown in Figure  21. Unlike the highly damped motion seen at S = 0.56, the corresponding PVC mode at S = 0.66 is significantly more unstable than the modes around it. Just above the dominant PVC frequency, two modes with decreasing growth rate but similar mode shape are observed. These motions are likely due to the weak nature of the PVC at this swirl number. While strong, stable PVCs manifest as limit-cycle oscillators with a single frequency, it is possible that weaker PVCs are more susceptible to fluctuating phenomenon such as turbulence. These fluctuations affect parameters such as local swirl number, which could impact the frequency of the PVC and lead to a frequency that varies in time and thus appears as multiple modes.   While DMD suggests that a weak PVC exists in the flow at S = 0.66, conclusions about the energy content of the PVC cannot be drawn directly from DMD. Instead, POD is used to quantify the fluctuating energy contained in the PVC motion to verify its relative strength in the flow field. As POD modes are energy-ordered, low-numbered modes are dominant in the flow field. POD modes 1 through 6 each contain between 1.3% and 2.2% of the total fluctuating energy and capture broadband, low-frequency activity related to vortex roll-up in the shear layers. Modes 7 through 10, which each contain around 1% of the total turbulent kinetic energy, are illustrated in Figure 22 and show signs of the PVC identified using DMD in both their frequency spectra and mode shapes, which can be compared to Figure  20a and Figure 21a. The spectra of modes 7 through 10 contain wide peaks around 850 Hz, which is in general agreement with the frequency of 858 Hz identified by DMD. In addition, the mode shapes roughly correspond to those identified using DMD, suggesting that both methods have identified the same motion. The relative energy content of the PVC modes produced by POD confirms the assertion that the PVC present at S = 0.66 is relatively weak. Though a PVC is present at this swirl number, the flow field continues to be dominated by the vortical structures contained in POD modes 1 through 6 produced due to axial and azimuthal Kelvin-Helmholtz instability and the inertial instability caused by the rotating core.  Previous work 38 has shown that at high swirl numbers, i.e., S > 1.05, stronger PVCs exist than those presented in this study. Strong PVCs are characterized by high amplitude, narrow band peaks in POD frequency spectra. This contrasts the broadband response of developing PVCs, as shown in Figure 16. From an energy standpoint, strong PVCs are always contained in the highest energy modes of the flow field, typically just modes one and two, which have significantly higher energy than any following modes. In previous studies using the same experimental configuration, 38,48 the in-nozzle pressure transucers have been used to identify the presence of a PVC. At high levels of swirl, PVCs appear as narrow-band peaks in the measured pressure spectra, as indicated by the arrow in Figure 23. For S = 0.66, shown in red in Figure 23, no such peak can be identified despite the indication of DMD and POD that a PVC is indeed present. The small peak present in both spectra at around 1050 Hz is a system frequency present at all swirl numbers, and is not indicative of a vortical motion in the flow field.

IV. Conclusions
This study has examined the dynamics of a turbulent swirling jet at a range of swirl numbers from S = 0 to S = 0.66. Three flow regimes were identified based on the presence of VB in the jet: pre-breakdown, near-breakdown, and post-breakdown. The dynamics of the jet in each regime were examined using DMD, which revealed the development of motions along the jet centerline as vortex breakdown set in in the near-breakdown regime. The dynamics identified by DMD in this turbulent study closely followed those identified in laminar studies of vortex breakdown: axisymmetric motions preceded tilting vortices, which were ultimately followed by the appearance of both axisymmetric and antisymmetric motions in a recirculation zone along the jet centerline. The location of these motions is apparent in overall analysis of strong velocity fluctuations. High-frequency fluctuations were shown to be present closer to the nozzle, while low-frequency fluctuations appear farther downstream. An apparent wobbling of the jet was also identified in the near-breakdown regime and was characterized by fitting a Rankine vortex model to the velocity fields. This wobble is not frequently discussed in swirling jet studies, and its effects could warrant further study as its presence coincides with the onset of vortex breakdown. PVC-like motions were found at a swirl number of S = 0.56, and a weak PVC was identified by DMD in the post-breakdown regime at a swirl number of S = 0.66. Its energy content was analyzed using POD.
Notably, the PVC identified by DMD was found at a significantly lower swirl number than PVCs previously identified in the same flow field. PVC-like motions were identified at a swirl number as low as S = 0.56, far below the point at which PVCs were observed using nozzle pressure measurements in previous studies. This suggests that innozzle pressure measurements are not sufficient means with which to identify weak PVCs, and that further analysis of the velocity fields is needed to confirm or deny the presence of a PVC. In particular, the DMD algorithm developed for this study is a convenient tool for this analysis due to its ability to identify single-frequency motions like PVCs. Future work will employ an improved DMD algorithm at higher swirl numbers to study the dynamics of stronger PVCs.