The analysis of the near-field dynamics of an aircraft wake and its interaction with an engine jet exhaust is of primary interest in applications covering a wide spectrum of aerospace technology. Examples range from the characterization of the structure of persistent and hazardous trailing vortices during take-off and landing phases  , to the investigation of the impact of pollutant emissions on the atmospheric environment  .
The phenomenology of the aircraft wake is usually described in terms of the downstream distance behind the aircraft wing: the wake evolution can be divided into four characteristic regions as proposed by Jacquin et al.  , namely:
− the near-field wake ( , being the mean aerodynamic wing chord) is characterized by the roll-up of the vorticity sheet emanating from the wing trailing edge and the presence of several concentrated vortices (wing tip, flap, nacelle, fuselage, horizontal tailplane, ∙∙∙);
− the extended near-field wake ( ) where the roll-up completes and the primary co-rotating vortices interact and ultimately merge;
− the mid- to far-field wake ( ) is generally composed of a pair of symmetric counter-rotating vortices that descend in the atmosphere due to their mutual velocity induction. In this region the vortices develop linear instabilities such as the Crow instability  , although other unsteady processes may take place such as vortex wandering, a poorly understood phenomenon  , or Crouch-type instabilities  that occur when the wake is composed of more than one pair of vortices;
− the dispersion regime ( ) is characterized by the final decay of the vortex system―after the occurrence of the instabilities in the far-field―which is controlled by the local ambient conditions (i.e. atmospheric turbulence and thermal stratification).
The present study is an attempt to understand the dynamics of the wake flow in the near- to the extended near-field, its intrinsic instabilities and the effects of jet engine exhaust on the development of the wake vortex structure.
In the initial stages of its formation, the aircraft wake is a complex vortex system composed of multiple interacting counter- and co-rotating vortices. The presence of two co-rotating vortices is a common feature in the wake of an aircraft in high-lift configuration. A pair of distinct two-dimensional co-rotating vortices experiences merging   , under the condition that the ratio of the core size and the distance between the vortices a/b exceeds a certain critical threshold value  . Although, initially the ratio a/b of the primary co-rotating vortices is often smaller than this value, the vortex core size is increased by the roll-up of the vorticity sheet. Furthermore, the velocity induced by the remnants of the vorticity sheet and the multiple vortices present in the wake leads to the decrease of the distance between the vortices. This increases the ratio a/b resulting in the merger of the vortices.
Depending on their separation distance and the respective core sizes, two co-rotating vortices can be unstable due to the strain field that each vortex exerts on the other. This results in a short-wavelength instability that is characterized by a three-dimensional sinusoidal deformation of the core and exists both for counter-rotating  ) and co-rotating vortices    ). The first objective of this study is to determine if a realistic aircraft wake permits the development of short-wavelength instabilities and how the instability mechanism affects the final vortex structure.
The qualitative features of the interaction between an engine jet and a single wake vortex were illustrated by Miake-Lye et al.  who identified two distinct phases. During the first few seconds after the emission the jets rapidly mix with ambient air (jet regime) while the vorticity shed from the wings rolls up into a pair of trailing vortices. Later on, the dynamics are dominated by the entrainment of the jets into the vortex flow (interaction regime). This pioneering analysis was followed by the analytical and experimental studies by Jacquin and Garnier  , Gerz and Ehret  and Brunet et al.  . A Direct Numerical Simulation (DNS) of jet/vortex dynamics and mixing was performed by Ferreira Gago et al.  who observed the development of counter-rotating structures of azimuthal vorticity in a region around the primary vortex where the engine exhaust gases concentrate. Paoli et al.   performed large-eddy simulations (LES) of jet/vortex interaction for different engine/wake vortex configurations. They showed that the dynamics are strongly dependent on their relative position and velocity ratio: if the jet and the vortex are close together, the jet blows inside the vortex, while, if they are initially well separated, the jet is entrained by the vortex and wraps around its core. The emergence of azimuthal vorticity structures induced by the jet axial velocity made the flow-field highly tridimensional, especially in the blowing case. In both cases, the initial vortex induced velocity field was described by an analytical relation for an axisymmetric vortex.
This suggests the possibility to extend the analysis to different configurations with complex distributions of the initial vorticity. Experimental investigations of the interaction between a jet and the vorticity sheet shed by the wing were performed by Wang and Zaman  who observed that the jet experienced stretching and vertical compression due to the dynamics of the neighboring vortices which ultimately increased the decay of axial velocity. An interesting tentative to simulate the jet/vortex sheet interaction in a realistic aircraft wake was recently done by Fares et al.  with a sophisticated Reynolds averaged Navier-Stokes solver, keeping, however, the intrinsic limitations of steady flow assumption.
Due to the restriction of direct numerical simulations to low Reynolds number flows, large-eddy simulations are a natural candidate to represent these inherently unsteady phenomena at a high Reynolds number, and has then been used in the present study.
The paper is organized as follows: the governing equations and the modeling are first described; the results of the 2D stable and unstable near-field dynamics as well as its interaction with an engine jet exhaust follow; conclusions and the main outcomes are finally given.
2. Governing Equation, Modeling and Numerical Tool
In the LES approach the Navier-Stokes equations are filtered spatially, such that any variable may be decomposed into a resolved (or large scale) part and a non-resolved (or subgrid-scale) part , with . This procedure may be obtained by a convolution integral of the variable with a filter function depending on a filter width . Practically, the filter width is simply given by the computational mesh cell size . For compressible flows, Favre-filtered variables, defined as , with are used. The dimensionless Favre-filtered equations are:
where the subgrid-scale (SGS) stress tensor , and the SGS heat flux are to be modeled, and where the following classical approximations have been made:
− The Favre-filtered shear stress tensor is identified with the filtered shear stress tensor
− The Favre-filtered heat flux is identified with the filtered heat flux
− The filtered kinetic energy term in the energy equation is approximated by , where is the kinetic energy.
The Favre-filtered passive scalar equation is:
The SGS momentum, , the SGS heat flux, , and the SGS scalar flux, , are modeled through subgrid-scale eddy-viscosity concept:
where is the SGS dynamic viscosity, is the large scale strain rate tensor and is the turbulent Schmidt number; while is the turbulent Prandtl number, defining the modified temperature , where is the specific heat at constant volume.
The SGS viscosity model is based on the Structure Function model  initially developed in spectral space (effective viscosity model) and then transposed into the physical space. The expression of the Structure Function is
where is the cutoff length, and where denotes spatial averaging, here over the sphere of radius . As the information brought by the model is local in space, it leads to a poor estimation of the kinetic energy at the cutoff, which may be improved by a suitable filtering  , in order to remove the influence of the large scales on the SGS viscosity. The procedure defined by Ducros et al.  is to apply (possibly n times) a discrete Laplacian high-pass filter to the velocity field before calculating the Structure Function. The optimum value of n found by Ducros et al.  for their simulations is . This value has also been used here. Finally, the Filtered Structure Function model reads:
where the superscript (n) indicates that the filter has been applied n times. The value of α used here is . The Structure Function model formulation of Métais and Lesieur  in the spectral space insures that the SGS viscosity vanishes when there is no energy at the cutoff wavelength. This property is particularly important for the simulation of transitional flows as discussed in recent high Reynolds LES of jet/vortex interactions  and of the elliptic stability of a vortex pair  . The reference quantities used for nondimensionalization are reference length ; density ; velocity , the reference speed of sound; pressure ; temperature ; dynamic viscosity ; and specific heat . Expect when explicitly mentioned, hereafter all variables are shown normalized by these reference quantities. In particular, the non-dimensional velocity
where M is the local Mach number. Hence, for a flow with small temperature variations ( ), the non-dimensional velocity corresponds effectively to the local Mach number, .
Numerical Tool NTMIX3D and Initialization Procedure
The numerical code   used in this study is a parallel, three-dimensional, finite differences Navier-Stokes solver. The space discretization is performed by a sixth-order compact scheme  for both convective and viscous terms. Time integration is performed by means of a three-stage third-order Runge-Kutta method.
Temporal simulations were carried out to analyze the evolution of the near-field wake dynamics and mixing. This is based on the assumption of a locally parallel flow, which means that the gradients of the mean flow in the axial direction are neglected over the short distance corresponding to the axial dimension of the simulation domain. Instabilities developing in the simulated flow are automatically of convective nature, and we may not capture absolute instabilities   .
Periodic boundary conditions are used along the vortex/jet axis z and in the upper-lower boundaries (y-direction), while symmetry boundary conditions were used in the spanwise x-direction. The boundary in Figure 1 is an effective symmetry axis (only half of the aircraft wake is simulated); the use of symmetry on the opposite side is a numerical artifact necessary to keep any incoming flow from entering the computational domain. This implies the presence of an infinite number of “image vortices”, however, the influence of
Figure 1. The initial condition with the isocontours of the x-component of the velocity containing the measurement window.
those vortices was significantly reduced by placing the boundary sufficiently far from the region where the vortex dynamics takes place (as verified a posteriori in the simulations).
The initial condition for the wing-generated wake used in the computations is obtained by 5-hole probe data from a windtunnel measurement behind a half-model of a typical transport aircraft. The data comes from a measurement plane located at , where is the downstream distance normalized by the wing span of the model B. In the temporal simulations, the downstream distance is simply related to the physical time of the wake by where is the free-stream speed of the model. As the size of the measurement plane was too small to simulate the wake vortex dynamics in the near-field to extended near-field, interpolation and extrapolation procedures were required. This procedure provided the initial condition on a computational grid containing the measurement window that was sufficiently large to avoid the influence of the domain boundaries. An interpolation routine was used to interpolate the measured velocity field on the mesh with the desired spatial resolution. Furthermore, an extrapolation routine was used to reconstruct the aerodynamic flow field outside of the measurement window. This was done by using analytical functions that describe the velocity field of the vortices outside of the measurement window, while minimizing the velocity gradients to avoid the generation of spurious vorticity at its borders. The resulting fields are denoted by and , respectively. As an example, Figure 1 shows the distribution of .
The Reynolds number of the measurements based on the circulation contained by one half of the wake has been retained for the numerical simulations, i.e. .
3. Results and Discussion
3.1. Two-Dimensional Results
A two-dimensional computation was first carried out to study the near- to extended near-field dynamics. The size of the computational domain size used for this simulation is and the number of grid points is -note that only half of the aircraft wake was simulated as described above.
For the sake of validation, the numerical results at (see Figure 2) were compared to data in a measurement plane at the same downstream distance (not shown here). It was verified that the positions of the main vortices agree with those measured, however the finer structures visible in the numerical results were not observed in the experimental data, because of the weak resolution of the measurement images.
At , the two closely spaced vortices that are horizontally aligned at experience a strong interaction and finally coalesce. This occurs at and results in a vortex wake (Figure 3) that is composed of two primary co-rotating vortices (denoted by I and II) and a smaller secondary vortex (III) that orbits around the main “dipole”.
The dynamics of the primary co-rotating vortices is characterized by the merging process, which is generally a fast, two-dimensional and stable phenomenon as described in detail by Cerretelli et al.  .
Figure 2. Isocontours of vorticity magnitude for the computed flow field at a downstream distance of .
Figure 3. Isocontours of vorticity magnitude for the flow field obtained by a two-dimensional computation at a downstream distance of .
Meunier et al.  identified a critical ratio above which two-dimensional merging occurs, for a symmetric dipole of Gaussian vortices. The primary co-rotating vortex pair observed in Figure 3 is characterized by vortex core ratios and and the circulation ratio between the two vortices is found to be . It is expected that a dipole of Gaussian vortices with such characteristics will not rapidly merge at a flight Reynolds number and the vortices will rotate around each other.
Under certain circumstances, depending on the Reynolds number and the arrangement of the two co-rotating vortices ( ), a system composed of a co-rotating vortex system may be subject to the development of the three-dimensional elliptical instability (see a recent review by Kerswell  ), as a result the merging process becomes unstable. This phenomenon is analyzed in the next Section.
3.2. Three-Dimensional Results―No Jet
The remnants of the vorticity sheet and the concentrated vortices orbiting around the primary dipole of the wake system bring the co-rotating vortices rapidly together (through velocity induction) and ultimately cause the merging of the vortex system. This merging process is often fast in a real aircraft wake and usually stable in the absence of strong perturbations. One might therefore think to trigger the elliptical instability by injecting a significant amount of energy in the corresponding unstable mode using some forcing techniques. On the other hand, one of the motivations of this study was to understand if the elliptical instability can emerge “naturally” in an aircraft wake. Thus, instead of forcing selectively the elliptical instability mode by injecting strong perturbations, a weak random white noise is imposed to simulate the intrinsic dynamics of the flow that is susceptible to instabilities. The unstable wavelength corresponding to the elliptical instability being unknown a priori, the axial length of the computational domain was chosen to be significantly larger than the core size a, namely one wing span (note that , see e.g. Kerswell  ). Taking an axial domain length that is larger than 5 instability wavelengths limits confinement effects caused by the computational domain  , which allows the unstable mode to naturally amplify.
The simulation were performed in a computational domain given by and a number of gridpoints in the three directions giving a total of approximately gridpoints. This large number of gridpoints is necessary to capture the short-wavelength/elliptical instability with sufficient resolution to avoid numerical diffusion and damping effects on the development of the instability.
In order to trigger the flow instability, a weak random white noise is superimposed to the base flow as follows:
where is the amplitude of the perturbations and rand a random generator with .
The global dynamics is initially governed by the same phenomena discussed in the previous Section, namely the roll-up of the vorticity sheet, leading to the generation of several concentrated vortices. Subsequently, a reduction of the number of vortices through the coalescence of co-rotating vortices is observed. However, some of the smaller vortical structures are deformed three-dimensionally, become turbulent and fall apart. At a downstream distance of approximately the vortex system is composed of two primary co-rotating vortices and a secondary vortex orbiting around it.
This is shown in Figure 4 which reports the contours of the vorticity magnitude averaged in the axial (z-) direction. The primary co-rotating vortex pair is characterized by the two core ratios (normalized by their separation distance) and , respectively, and by a circulation ratio .
At a distance larger than the secondary vortex (III) shows a sinusoidal short-wavelength deformation of the vortex core indicating the elliptical instability (see Figure 5(a)). The elliptical instability is explained by the fact that an elliptical flow is unstable to three-dimensional disturbances  . The presence of the co-rotating dipole exerts a significant external strain on the vortex III causing the streamlines to deform in an elliptical manner, which
Figure 4. Isocontours of vorticity magnitude for the mean flow field obtained by a three-dimensional computation at a downstream distance of .
Figure 5. Isosurfaces and isocontours of vorticity magnitude as a function of downstream distance . Top panel: ; bottom panel: .
makes it unstable. The occurrence of the elliptic instability is suspected because of the short-wavelength disturbance, which is of the order of the core size. Furthermore, a spectral decomposition of the kinetic energy in the axial direction shows the exponential amplification of the energy contained by the Fourier mode which corresponds to the wavelength of the elliptical instability (see Figure 6). The instability saturates in vortex III resulting in a local transition to turbulence leaving behind a completely noncoherent vortex structure (see also Figure 5(a)).
At a downstream distance larger than the dynamics of the system is mainly governed by the behaviour of the two co-rotating vortices I and II. They rotate around each other and, due to the relatively small core ratios ( ), the elliptical instability has time to develop. A clear sinusoidal deformation of the vortex two vortex cores is observed in Figure 5(b), which is somewhat asymmetric resulting from the unequal core to spacing ratios of the individual vortices. The strong deformation of the vortex cores due to the instability leads to the exchange of vorticity between the two vortices and results in the unstable merger.
The occurrence of the elliptical instability is confirmed by a spectral analysis. The spectral decomposition in the axial direction shows the amplification of the energy contained by the most unstable mode corresponding to the elliptical wavelength, which is depicted for vortex I and II in Figure 7. After a transient regime one can clearly observe an exponential growth which indicates the linear regime of the elliptic instability. In this phase, the extracted growth rate (normalized by the turnover time of the vortex system) is approximately the same in the two vortices, i.e. , for the modes corresponding to the most amplified wavelengths in each vortex (respectively and ). An explicit formula was obtained by Le Dizès et al.  through a theoretical analysis that predicts the temporal growth
Figure 6. Evolution of the most amplified Fourier mode ( ) of the perturbation energy contained in vortex III.
rate of the elliptical instability occurring in each vortex as a function of the circulation and , the vortex core radii and , and the Reynolds number . Evaluating the parameters and of the resulting vortex system at a distance and inserting them in the Formula (6.1) of Le Dizès et al.  provides the nondimensional growth rate . This theoretical prediction agrees fairly well with the results obtained in the simulations of the aircraft wake as shown in Figure 8. The slight difference can be explained by the fact that the theory assumes the two vortices are parallel and Gaussian  , while the structure of the vortex core in the simulations resembles the multiple core scale structure as suggested by Jacquin et al.  . This is further confirmed by the analysis of Fabre et al.  who stated that the properties of the cooperative short-wavelength instabilities of vortices representative of aircraft wakes may be differ from those of Gaussian vortices.
Figure 7. Evolution of the most amplified Fourier mode ( ) of the perturbation energy contained in vortices I and II.
Figure 8. Comparison between simulation and theory  for vortex I and II in the configuration at . The simulated growth rates are shown by circles, the first three instability branches are displayed by solid lines (the vertical dash-dotted lines indicate the wavelengths ( ) that were accessible in the simulation).
As a final remark, the development of the elliptic instability during the merging process modifies the final vortex structure. In particular, as shown in Figure 9 and Figure 10, it causes the formation of a larger vortex core and leads to a decrease of the tangential velocity.
3.3. Three-Dimensional Results―Jet
In the previous Section, the formation and evolution of an aircraft wake was analyzed by focusing on the development of the intrinsic elliptic instability of the vortex system. In this Section, a different situation is considered where the vortex sheet interacts with an engine jet. The aim is first to find out whether and how the jet modifies the roll-up process of the vorticity sheet and, secondly, to analyze the mixing of an exhaust passive scalar in the wake.
Figure 9. Isocontours of u at a downstream distance of . Top panel: 2D simulations of the vortex sheet; bottom panel: 3D simulations of the vortex sheet, averaged in the axial direction (the line where the cut through the core was made is indicated).
Figure 10. Velocity profiles through the center of the merged vortex at (data extracted at locations indicated by the line in Figure 9): solid line, 2D simulation; dashed line, 3D simulation.
The computational domain is and consists of gridpoints which are equally spaced in each of the three directions. The first two dimensions are the same as the previous Section, while the axial length corresponds to the most unstable jet velocity profile in the simulations of Michalke and Hermann  (see also Ferreira Gago et al.  , Paoli et al.  and Paoli et al.  ). The jet axial velocity and the scalar field (representing any exhaust passive species) are initialized with a tanh
profile, , where is the jet radius, is the radial distance from the center and is the
displacement thickness of the jet (see Paoli et al.  for details). With these profiles, one has
where subscripts a and e indicate, respectively, the free-stream conditions and the conditions at the center of the exhaust jet. In the present study, there is no co-flow, and , while the exhaust velocity and scalar are, respectively, and .
The velocity field in Equation (13) is then added to Equation (12) to trigger the jet instability and its transition to turbulence, while, at the same time, it interacts with the vorticity sheet. Two main features characterizes this interaction as shown in Figure 11 which reports the isosurface of the vorticity
Figure 11. Isosurface of obtained by a three-dimensional computation of the jet/sheet system, at a downstream distance of .
magnitude at a late stage of the roll-up process, . First, the jet is entrained around the multiple vortex system arising from the roll-up. Secondly, the jet continuously exchanges its axial momentum with the sheet by vortex stretching, causing the formation of complex three-dimensional vortical structures. Note that the figure displays the vorticity level ( being the instantaneous maximum vorticity) which identifies the core of a axisymmetric Lamb-Oseen vortex. Indeed, the core of the merged vortex is clearly visible as well as the fine-scale eddies, remnants of the initial sheet vorticity and the jet turbulence, which are still wrapping around the core. This was already observed in previous simulations of jet/vortex interactions   , although in those cases, only a few azimuthal vortical structures (vortex rings) were formed around the core and only came from jet-generated turbulence. In the present case, the wake dynamics is more complex because the jet interacts with a vorticity sheet rather than a single well-formed vortex.
This strong momentum exchange and the turbulent diffusion induced by the jet affects the tangential velocity field of the final vortex. For example, the isocontour lines reported in Figure 12 shows the decrease of the (z-average) u-component of the velocity and the increase of the vortex core. The latter is due to the conservation of circulation (not shown): as the peak azimuthal velocity decreases the core size must increase to keep constant. The above picture is confirmed by the u-profiles in Figure 13, taken along the line marked in Figure 12, which shows that, at the u-peak is about half the value it had in the 2D case.
To conclude the analysis of the jet/sheet interaction, the isocontour lines of the passive scalar Y are reported in Figure 14 at . The figure shows
Figure 12. Isocontours of u at a downstream distance of . Top panel: 2D simulations of the vortex sheet; bottom panel: 3D simulations of the jet/sheet system, averaged in the axial direction (the i-line through the core is also indicated).
Figure 13. Velocity profiles through the center of the merged vortex at (data extracted at locations indicated by the line in Figure 12): solid line, 2D sheet simulation; dashed line, 3D jet/sheet simulation.
Figure 14. Isocontours of scalar field at a downstream distance of : top, 2D simulation of the vortex sheet; down, 3D simulations of the jet/sheet system, averaged in the axial direction (the initial location of the jet is also represented by thick lines).
that the combined effects of turbulent diffusion of the jet and its entrainment within the vortex sheet, enhances the dispersion of exhaust gases in the wake. This results in a larger plume area and causes a stronger dilution of the scalar field, which reaches a value of approximately 0.005, at a final stage of the merging process, compared to 0.476 in the the 2D laminar case.
The formation and near-field dynamics of an aircraft wake was studied by means of Large Eddy Simulations, using available experimental data of a wing-generated vorticity sheet. The roll-up process and the stability properties were first analyzed by imposing a random white noise and the occurrence of the short-wavelength instability in the principal co-rotating vortex pair was evidenced. This caused the unstable merging of the vortices and modified the core of the final vortex.
Furthermore, the roll-up and the formation of the wake were simulated by adding a model jet flow. This affected the dynamics of the primary co-rotating vortices during the roll-up, mainly by vortex stretching due to the momentum exchange with the vorticity sheet. The consequences for the merging were the decrease of the azimuthal velocity and the enlargement of the vortex core. Finally, the jet/sheet interaction affected the exhaust passive scalar dilution, leading to more efficient mixing and a larger plume area at the end of the merging process.
The authors wish to thank both CINES (Centre Informatique National de l’Enseignement Supérieur) and IDRIS (Institut du Développement et des Ressources en Informatique Scientifique) for providing the CPU hours. The authors also thank Airbus-Deutschland for providing the wind tunnel data.