Numerical investigation of plasma-controlled turbulent jets for mixing enhancement

Plasma-controlled turbulent jets are investigated by means of Implicit Large–Eddy Simulations at a Reynolds number equal to 460,000 (based on the diameter of the jet and the centreline velocity at the nozzle exit). Eight Dielectric Barrier Discharge (DBD) plasma actuators located just before the nozzle exit are used as an active control device with the aim to enhance the mixing of the jet. Four control configurations are presented in this numerical study as well as a reference case with no control and a tripping case where a random forcing is used to destabilize the nozzle boundary layer. Visualisations of the different cases and time-averaged statistics for the different controlled cases are showing strong modifications of the vortex structures downstream of the nozzle exit, with a substantial reduction of the potential core, an increase of the jet radial expansion and an improvement of the mixing properties of the flow.


Introduction
Turbulent jets are used in a variety of industrial applications, such as jet engines and combustion chambers. Greenhouse gases, toxic pollutants, heat ejection and sound radiation emitted from such devices are often detrimental to the environment. Enhancing specific properties of a jet is therefore vital. For instance, improving its mixing property would result in higher thrust for a jet engine and more energy extraction due to a more complete combustion in the combustion chamber. A comprehensive review of research activities for the control of turbulent jets in the last 50 years can be found in Zaman et al. (2011).
Strategies in order to enhance the performance of turbulent jets fall in two control categories, passive and active. Passive control usually involves geometric modifications of the nozzle by using notched nozzles (Pannu and Johannesen, 1976), tabs (Bradbury and Khadem, 1975;Ahuja and Brown, 1989;Samimy et al., 1993) and chevrons (Brausch et al., 2002;Callender et al., 2005;. Even though no external energy is added to the flow, these devices are always present which could result in performance penalties. For example, in the case of an airplane during cruise, a jet engine with chevron nozzles can experience thrust penalty (Zaman et al., 2011;Calkins and Butler, 2004). At first, control studies were directed into varying the geometry of the nozzle, i.e. elliptic (Husain and Hussain, 1983;Hussain and Husain, 1989;Ho and Gutmark, 1987), or rectangular nozzles (Gutmark et al., 1989;Quinn, 1991). Experiments were performed in Ho and Gutmark (1987) using a small-aspect-ratio elliptic nozzle to demonstrate that the entrainment ratio can be greatly enhanced by comparison to a circular or a planar jet. The flow was seen to experience axis switching: the major and the minor axes of the elliptic cross-section of the jet were alternating along the downstream direction resulting in a larger engulfment of the surrounding fluid into the jet. Hot-wire anemometry experiments were performed in Quinn (1991) with rectangular nozzles of aspect ratio 2 and 10. The authors observed an increase in the near field mixing due to higher shear-layer values of the turbulent kinetic energy and the Reynolds stress, as well as a shorter potential core length for a nozzle with aspect ratio 10.
Boeing, General Electric, and NASA have developed serrated edges called chevrons for the back of the nacelle and the engine exhaust nozzle and they found that chevrons can reduce jet noise up to 4 dB, however associated with a reduction of thrust and a loss of 0.25% on fuel consumption (Saiyed et al., 2000;Nesbitt et al., 2002;Bridges and Brown, 2004). Experiments with a coaxial flow test rig were carried out in Callender et al. (2005Callender et al. ( , 2008 to study different chevron nozzles over a wide range of operating conditions. The numbers of lobes and levels of penetration were varied in order to evaluate the impact of these geometric parameters on the noise level. All configurations achieved a reduction of 3-6 dB for the overall Sound Pressure Level (SPL). Calculations of perceived noise level directivity also showed a 4-6 dB reduction at aft angles. It was also observed that the chevron penetration was the primary factor to control the tradeoff between low-frequency reduction and high-frequency SPL increase and to influence the size and intensity of the noise region near the nozzle lip. Hybrid Reynolds-averaged Navier-Stokes (RANS) -Large-Eddy Simulations (LES) of chevron jet flows were performed in Xia et al. (2009) in order to generate noise predictions. It was found that for a Mach number equal to 0.9, the numerical data compare favorably with measurements for the flow field, with encouraging agreement of the predicted far field sound directivity and spectra obtained using the Ffowcs Williams and Hawkings (FWH) surface integral method. The main conclusion was that the chevron penetration angle is a critical parameter to achieve noise mitigation.
To avoid the disadvantages of passive modifications, most of the research related to turbulent jets is nowadays focusing on active control solutions for which energy is only added to the flow when needed. In a similar way to passive control solutions, they are designed to manipulate the topology of the flow, either by provoking instabilities or by directly targeting the destruction or creation of large-scale structures. Various strategies have been studied for active control, such as synthetic or piezoelectric jet actuators (Butler and Calkins, 2003;Low et al., 2010;Önder and Meyers, 2014), secondary jets (Lardeau et al., 2002;Maury et al., 2009;Gautier et al., 2014) and plasma actuators (Samimy et al., 2004;2007a;2007b;Kim et al., 2009;Gaitonde and Samimy, 2010;. Piezoelectric actuators for turbulent jet control were assessed in Butler and Calkins (2003) to investigate if they could alter the turbulent energy distribution. These actuators produced small-scale disturbances in the shear layer just at the nozzle exit. Particle Image Velocimetry (PIV) snapshots showed an increase in vorticity in the near field while an array of microphones showed a shift of the peak values for low emission angles. This shift was attributed to the fact that the large coherent structures were annihilated, in conjunction with a good anisotropy for the small scales. A different approach was taken in Low et al. (2010) where pressure readings downstream of the nozzle were used to conduct open and closed loop control tests using synthetic zero net-mass-flux (ZNMF) actuators radially placed at the nozzle exit. The authors managed to show the ability of their control technique to modify the near field region flow features but with very little impact of the far field noise spectra. The near field of a ZNMF actuated round jet using Direct Numerical Simulations (DNS) was studied in Önder and Meyers (2014) for a low Reynolds number of 2,000 (based on the diameter of the jet and the centreline velocity at the nozzle exit). Strong deformations of the near-field jet region were observed which were very similar to those observed for non-circular jets. These changes were attributed to the self-deformation of the jet's primary vortex rings due to distortions in their azimuthal curvature and by the production of side jets by the development and subsequent detachment of secondary streamwise vortex pairs. Secondary jets with a mass flow rate of 10% of the main jet at an angle of 45°were used in Lardeau et al. (2002) with DNS of a turbulent jet at low Reynolds numbers equal to few thousands. A reduction of the potential core length associated with a significant jet expansion were observed and, surprisingly, the injections of fluid by the secondary jets did not result in a big enstrophy increase for the main jet. The mixing of the flow was shown to increase, especially when the secondary jets were pulsating. Steady and unsteady fluidic actuators, in the form of secondary control jets injecting from the nozzle lip were investigated experimentally in Maury et al. (2009Maury et al. ( , 2011. The actuation comprised an azimuthal distribution of 16 nozzle-lip mounted microjets, injecting fluid at a penetration angle of 60°. Different geometrical configurations were explored by varying the distance between the microjets and in particular microjets that converge in pair. The authors concluded that steady and unsteady forcing affect differently low order statistical moments of the velocity field and that the response of the flow to unsteady forcing appears to comprise a non-linear component, at the main forcing frequency, and two secondary components that compare well with predictions of linear stability theory. Up to 2.4 dB of global sound reduction was reported in these experiments. The same number of microjets in a converging configuration, called fluidevrons, was used by Gautier et al. (2014) to perform DNS for a jet at a low Reynolds number equal to 10,000. It was shown that the secondary jets can destroy the large annular structures generated at the exit of the nozzle, accelerating the generation of smaller scales in the near field, resulting in an increased turbulent kinetic energy. A distinct flow pattern was observed with some ejections from the main jet into the surrounding fluid, and horseshoe vortices generated by the microjets in the near field. As a result, the potential core was observed to increase in length.
Another form of active control for turbulent jet is based on plasma actuators, which are small devices that use high electric potential to accelerate portion of the flow field. Three main types have already been used for turbulent jets, Localised Arc-Filament Plasma Actuators (LAFPA), Plasma Synthetic Jet (PSJ) actuators and more conventional Dielectric Barrier Discharge (DBD) plasma actuators, as seen in Fig. 1. A LAFPA consists of two pin electrodes side-by-side on the wall and when an electric potential difference is applied between them an arc-filament plasma is created. A high temperature, high pressure perturbation is formed, which acts in a similar way to a tab (Samimy et al., 2004). This perturbation can destabilise the boundary layer inside the nozzle, resulting in vorticity generation and triggering of instabilities (Utkin et al., 2006;Samimy et al., 2007b). A PSJ actuator is a zero-net-massflux device mainly composed of two electrodes embedded in a cavity in connection with the external medium, with the help of a small dedicated orifice. By applying a voltage difference, an electrical arc is created between the two electrodes, leading to an increase in the internal energy. Since the air is confined, the temperature and pressure increase very quickly inside the cavity, producing a pulsed air jet (Caruana et al., 2013;Laurendeau et al., 2015). DBD actuators are based on a high electric potential difference between two electrodes separated by a dielectric material (Moreau, 2007;Corke et al., 2010;Thomas et al., 2009). The first electrode is positioned above the wall of the nozzle and exposed to the fluid flow, while the other is embedded inside the wall. When the electric potential is applied, the air flow is ionised above the exposed electrode and then accelerated along the embedded electrode, while it is also drawn nearer the wall, creating a wall jet effect (Moreau, 2007;Corke et al., 2010;Thomas et al., 2009).
LAFPAs were studied extensively by Samimy's group in the US as an active control solution for turbulent jets at high-Reynolds numbers for high-subsonic and supersonic regimes. Only LAFPAs can provide excitation signals of high amplitude and high frequency for high-speed and high-Reynolds-number flow control. The control strategy is based on the excitation of various instabilities and azimuthal modes of the jet.  V. Ioannou, S. Laizet International Journal of Heat and Fluid Flow 70 (2018) 193-205 Modes can be excited by using a step signal with period based on the sinusoid where f is the frequency of the pulse, t is the time, m is the azimuthal mode and θ i the location of the plasma actuator device. Mode zero = m 0 is called the axisymmetric case, flapping or first mixed modes corresponds to m ± 1 and second mixed modes correspond to m ± 2. Several experiments in this research group have explored different Mach numbers, excitation modes, and frequencies by capturing qualitative visualizations, as well as quantitative PIV, pressure, and noise level measurements (Samimy et al., 2004;2007a;Kim et al., 2009;Samimy et al., 2012). It was found that the first flapping mode at a frequency close to the preferred jet column mode ( = St 0.3) has the most effect on reducing the jet core length and increasing the mixing of the jet at the end of the potential core. Overall, the results revealed that the jet flow field and acoustic far field can be dramatically altered, providing a powerful control tool for high-speed and high-Reynoldsnumber (of the order of one million) jets. These trends were also confirmed in Benard et al. (2009) for moderate Reynolds numbers (of the order of 60,000) in a study dedicated to the control of turbulent jets using DBD plasma actuators. These experimental results were then complemented with numerical results by the same research group (Gaitonde and Samimy, 2010;Speth and Gaitonde, 2013) to provide a more detailed understanding of the physics of jet flows. Comparisons with experimental results indicate that the computational model (LES on unstructured meshes with the LAFPAs modelled as a surface heating condition) used in these studies is able to reproduce the main features induced by the actuators. Overall, the numerical results indicated a complex coherent structure generation, evolution, and disintegration process when the LAFPAs are used in order to control high-speed jets. Depending on the mode excitation m, various flow patterns were observed, like distorted elliptic rings alternating on either side of the jet axis for = ± m 1, elliptic vortex rings for = ± m 2, or axisymmetric roll-up events, with vortex rings and braids (Ho and Huerre, 1984) for = m 0. Wind tunnel experiments of turbulent jets with a diffuser nozzle and two DBD actuators were performed by Benard et al. (2008) at air speeds up to 40 m/s. It was observed that DBD actuators could successfully increase jet spreading and turbulent kinetic energy, associated with a decrease of the jet core length for a quasi-steady actuation. Unsteady actuations were also tried at Strouhal numbers ranging from 0.25 to 0.32, with an enhancement of the turbulent kinetic energy. More recently, the response of a isothermal jet to a PSJ actuator at a Mach number equal to 0.6 was investigated in Chedevergne et al. (2015) through PIV measurements and numerical simulations. The authors used a single actuator located at the nozzle lip to produce a synthetic microjet, leading to the generation of a large-scale coherent structure growing into the jet mixing layer. Satisfactory similarities were obtained between their experiments and their simulations using an unsteady RANS approach.
In the present work, turbulence-resolving simulations (ILES) of a turbulent jet at a Reynolds number of 460,000 (based on the diameter of the jet and the centreline velocity at the nozzle exit) are carried out with an array of eight DBD plasma actuators. Statistics and visualisations of the flow field are presented and discussed and the effect of the plasma actuators on the jet flow are investigated from a mixing perspective. A passive scalar study is carried out to determine which configuration of plasma actuators is more promising for the main jet mixing enhancement. For this first study, the number of plasma actuators is fixed as well as their location inside the nozzle.
The paper is organised as follows: first, we describe the numerical methods and flow configuration with details about the tools used to reach a high Reynolds number with high accuracy at a reasonable computational cost. Then, 2D and 3D instantaneous and time-averaged visualisations are presented to show how the DBD plasma actuators affect the flow field. Finally, probability density functions (PDFs) of the passive scalar are presented in order to study mixing enhancement for each controlled case. The paper ends with a conclusion and prospects for future work.

Flow configuration
We consider a turbulent jet at the exit of a rounded nozzle of internal diameter D. Part of the nozzle is included inside the computational domain L x × L y × L z as shown in Fig. 2

. The Reynolds number is
where U c is the velocity of the jet on the centreline at the exit of the nozzle. The coordinate system R is orthonormal with coordinate x in the streamwise direction and coordinates (y, z) in the transverse plane such that = = y z 0 on the centreline. The origin of R is located just at the exit of the nozzle on the centreline of the main jet, at a distance 3D from the upstream side of the computational domain. For the inflow boundary condition at = − x D 3 , the velocity profile inside the nozzle is based on a mean profile 2 which follows the profile in Bres et al. (2015).
A small co-flow velocity equal to 4%U c is imposed around the nozzle in order to avoid dealing with stagnant flow near the outlet, a potential critical situation for any outflow boundary conditions in incompressible flows. The no-slip boundary condition for the outside wall of the nozzle is ensured via the imposition of a laminar boundary profile at the inlet. At the outlet of the computational domain for = x L , x a simplified convection equation is imposed with a convective velocity equal to the local streamwise velocity at the centerline. For the lateral sides ( = ± y L /2 y and = ± z L /2 z ), modified Dirichlet boundary conditions are imposed, following a procedure described in Hasan et al. (2005). The idea is to allow the flow to enter the domain in order to mimic an entrainment mechanism. Note that for the computational domain used in the present study, the aspiration of fluid is marginal (less than 0.01U c ), with no significant impact by comparison to more classic freeslip boundary conditions. The boundary conditions for the passive scalar are similar to the velocity ones except for the inlet boundary conditions with a passive scalar set to zero in the co-flow. The simulations are performed with a timestep of 0.00125D/U c and the averaging of the flow is performed over 200D/U c , after the initial transients have died out.

Numerical methods
The Implicit Large Eddy Simulations (ILES) presented in this V. Ioannou, S. Laizet International Journal of Heat and Fluid Flow 70 (2018) 193-205 numerical study are performed using the high order flow solver Incompact3d (www.incompact3d.com). The following incompressible Navier-Stokes equations are solved, with an extra forcing term to account for the nozzle, the plasma actuators and the random tripping of the boundary layer inside the nozzle: p is the pressure field (for a fluid with a constant density ρ), u(x, t) is the velocity field and ν the kinematic viscosity. Note that convective terms are written in the skew-symmetric form as it allows the reduction of aliasing errors while remaining energy conserving (Kravchenko and Moin, 1997). In addition to the main equations, a passive scalar (ϕ) equation is also solved: where Sc is the Schmidt number (equal to 1 in this numerical study).
x y z mesh nodes, split into 1,024 computational cores. This resolution has been carefully chosen to allow the correct reproduction of the smallscales features of the jet and of the plasma actuators.
The open source flow solver Incompact3d is based on a Cartesian mesh and uses finite-difference 6th-order compact schemes for spatial discretisation and a 3rd-order Adams-Bashforth scheme for time-advancement. The main originality of Incompact3d is that the Poisson equation for the incompressibility of the velocity field is fully solved in spectral space via the use of relevant 3D Fast Fourier transforms (FFTs). With the help of the concept of modified wavenumber (Lele, 1992), the divergence free condition is ensured up to machine accuracy. The pressure mesh is staggered from the velocity one by half a mesh to avoid spurious pressure oscillations observed in a fully collocated approach (Laizet and Lamballais, 2009). The simplicity of the mesh allows an easy implementation of a 2D domain decomposition based on pencils (Laizet and Li, 2011). The computational domain is split into a number of sub-domains (pencils) which are each assigned to an MPI-process. The derivatives and interpolations in the x-direction (y-direction, z-direction) are performed in X-pencils (Y-pencils, Z-pencils), respectively. The 3D FFTs required by the Poisson solver are also broken down as series of 1D FFTs computed in one direction at a time. Global transpositions to switch from one pencil to another are performed with the MPI command MPI_ALLTOALL(V). Incompact3d can scale well with up to hundreds of thousands MPI-processes for simulations with several billion mesh nodes (Laizet and Li, 2011). Further details and validations of this flow solver can be found in Laizet and Lamballais (2009).
An Immersed Boundary Method (IBM) is needed in order to take into account the nozzle inside the computational domain based on a Cartesian mesh. Using the extra forcing term f in the Navier-Stokes equations, a zero velocity can be easily imposed inside the nozzle. For the present study and with the prescribed inlet boundary condition, the loss of continuity on the velocity field is negligible and has a minor impact on the solution. One of the cases in this numerical study includes an artificial random tripping inside the nozzle to trigger turbulence inside the boundary layer. It is based on the method developed by Schlatter and Örlü (2012) for planar boundary layers and adapted here for a round nozzle. The volumetric forcing, with a maximum effect located at = − x D 1.4 , occupies the circumference of the nozzle in the form of sinusoidal waves of period π 10 radians, and is added in the wallnormal momentum equation via the term f. It was suggested in Schlatter and Örlü (2012) that this tripping mechanism can be compared to a physical region with wall roughness, similar to what is typically used in experiments. The transition region generated by this tripping method can be significantly shorter than those generated by other mechanisms, such as Tollmien-Schlichting waves. Note that this forcing was successfully used in Incompact3d for the study of wallshear stress fluctuations in a turbulent boundary layer (Diaz-Daniel et al., 2017).

Implicit LES
Recently, a new method was implemented in Incompact3d in order to perform ILES. It is based on a strategy that introduces a targeted numerical dissipation at the small scales through the discretisation of the second derivatives of the viscous terms (Lamballais et al., 2011;Dairay et al., 2017). It was shown in these studies that it is possible to design a 4th-order finite-difference scheme in order to mimic a subgrid-scale model for LES based on the concept of Spectral Vanishing Viscosity (SVV, see for instance Tadmor, 1989;Karamanos and Karniadakis, 2000), at no extra computational cost.
In Incompact3d, the computation of second derivatives is achieved thanks to the following scheme In the framework of a Fourier analysis, it is well known that a modified square wavenumber ″ k can be related to this scheme with For a conventional 6th-order scheme, the coefficients of the discretisation are = α 2/11, = a 12/11, = b 3/11, = c 0. Using two conditions on the modified square wavenumber ), the scheme 5 can produce the following set of coefficients As explained in Lamballais et al. (2011), the extra-dissipation introduced by this discrete viscous operator can be interpreted as a spectral viscosity expressed as Using this expression, it is quite straighforward to adjust the coefficient in order to mimic the following SVV kernel . This value has been obtained following the procedure described in Dairay et al. (2017). It used a Pao-like spectral closure based on physical arguments to scale the numerical viscosity a priori.

Plasma actuators
The plasma actuator model used in the present study is a phenomenological model for Dielectric-Barrier-Discharge (DBD) plasma V. Ioannou, S. Laizet International Journal of Heat and Fluid Flow 70 (2018) 193-205 actuators which was initially developed by Suzen et al. (2005) and then implemented in Incompact3d for a skin-friction drag reduction study in a channel flow (Mahfoze and Laizet, 2017). The dimensions (in mm) of the DBD plasma actuator which forms the basis of the present model are shown in Fig. 3 (left). This plasma actuator was used in the experiments of Benard et al. (2015) and was used to validate the plasma actuator model implemented in Incompact3d (Brauner et al., 2016). The model is based on solving two equations to get the body forcing = f ρ E, c where E is the electric field and ρ c the charge density. There is one equation for the electric potential and another one for the charge density of the ionized working fluid. These equations are based on the assumption that the time variation of the magnetic field can be neglected, and therefore that the electric field is conservative (Suzen et al., 2005). The equations can be expressed as where ϕ is the electric potential, ϵ r the relative permittivity of the working fluid and λ d the Debye length. The force f is then added to the Navier-Stokes equations after making it non-dimensional. A full description of how to generate f with the values for the dimension of the plasma actuator, ϵ r and λ d as well as the correct set of boundary conditions, can be found in Suzen et al. (2005) and in Mahfoze and Laizet (2017). For the present study, the body forcing is extended azimuthally in order to be included inside the round nozzle at a distance 1D from the nozzle exit. Because this is the first attempt to control a turbulent jet, the location of the actuators with respect to the nozzle exit has been chosen arbitrarily and the influence of the location of the plasma actuators in the nozzle will be investigated in a future study. The forcing term f is non-dimensionalised using ρU 2 D 2 . The values used in the present study are based on an experiment in air ( = ρ 1.225 kg/m 3 ) with a nozzle diameter D equal to 0.20 m for a jet flowing at 33.796 m/s (corresponding to a Reynolds number equal to 460,000).
An array of eight plasma actuators is employed in the present study for the control action as shown in Fig. 3 (right). The full radial extent of the forcing corresponds to a height of 0.25D (discretised with 16 mesh nodes) and a streamwise extent of 0.7D (discretised with 30 mesh nodes) above the nozzle wall, while the 50% intensity contour is located at 0.045D. As a reference height, the sharp boundary layer inside the nozzle has a height of 0.07D. This display was inspired by Samimy's research group with their study of controlled supersonic jets with LAFPA devices (Samimy et al., 2007a;2007b;Kim et al., 2009). Preliminary investigations were carried out in order to determine the intensity of the plasma actuator forcing. The objective is to obtain a significant alteration of the flow field when the plasma actuators are on while having the same flow rate at the nozzle exit as in non-controlled cases. For the present study the increase of flow rate due to the plasma actuators was always kept under 0.75%. 4 cases of turbulent jets with different actuator control configurations are considered and compared with two uncontrolled simulations: one with no forcing corresponding to the baseline case (BS) and the other one with a random forcing to trigger instabilities inside the nozzle boundary layer (TR). The  V. Ioannou, S. Laizet International Journal of Heat and Fluid Flow 70 (2018) 193-205 controlled cases include three simulations with a continuous plasma actuator forcing multiplied by a factor 4, 8 and 12 (simulations I4, I8 and I12, respectively) as well as a simulation with a pulsating signal with a Strouhal number equal to 0.32 (corresponding to the preferred jet column mode, see Hussain and Zaman, 1981;Ho and Huerre, 1984 for more information) with a multiplying factor of 8 and a 50% duty cycle (simulation M0). In this numerical study, the increase in intensity of the forcing only changes its amplitude, not its spatial distribution.

Instantaneous visualisations
In the following, 2D and 3D visualisations are shown in order to get a first glimpse at the effect of the plasma actuator control on the flow field downstream of the nozzle exit. Instantaneous vorticity snapshots of the full flow field are presented in Fig. 4 where significant differences between the controlled and non-controlled cases can be observed. Snapshots from inside the nozzle are also presented in Fig. 5 where the structures generated by the plasma actuators can be seen when they are generated. For the baseline case BS, it can be seen that large ring-like structures are generated after an initial region of almost no shear layer growth with very little turbulence, region located just downstream of the nozzle exit. The large ring-like structures are then breaking down, giving rise to small scale turbulence. For the non-controlled tripping case TR, the initial region of almost no shear layer growth is reduced, but the apparition and the breaking down of the large ring-like structures are still present. The random tripping added in the nozzle boundary layer seems to decrease the strength of the rings due to a more intense small-scale turbulence activity at the nozzle exit, with in particular thin vortex tubes around the ring structure called braids (Ho and Huerre, 1984). The signature of the eight DBD plasma actuators can clearly be seen for the controlled cases when the activation is always on (I4 to I12 cases) with very intense elongated streamwise structures located at the exit of the nozzle up to the point where the shear layer starts to grow. It is interesting to notice that they are grouped by pairs, suggesting that they might be generated between two plasma actuators. Based on Fig. 4, it seems that the turbulence activity is more intense close to the nozzle for the control cases, with only very little vortical structures in the second part of the computational domain (at least for the same isosurface). As expected a pulsating motion can be observed for the M0 simulation with an alternance of turbulent puffs observable up to the end of the computational domain (see black vertical lines in Fig. 4). The frequency of these puffs corresponds to the jet preferred frequency (Hussain and Zaman, 1981;Ho and Huerre, 1984).  V. Ioannou, S. Laizet International Journal of Heat and Fluid Flow 70 (2018) 193-205 beneficial for mixing enhancement. Another notable difference with the non-controlled cases is the turbulence breakdown downstream of the nozzle is happening faster. The intense elongated vortices (braids) observed for the controlled cases can be seen from the inside of the nozzle in Fig. 5 and it can be assumed that they are responsible for the early destruction of the large ring-like structures. The white/grey colour in this figure represents regions where the streamwise velocity is high. It corresponds to the wall jets generated by the plasma actuators. Pairs of intense elongated vortical structures are generated between the plasma actuators by the ejection of streamwise counter-rotating vorticity. They transfer momentum away from the core of the main jet due to their relatively low speed. The higher the intensity of the forcing, the stronger the vorticity ejections are, as evidenced by thicker pairs of elongated tubes and more intense small-scale fluctuations at the nozzle exit, destroying almost instantaneously the large-scale ring structures. The shape of the potential core is discussed in the following section.

Time-averaged visualisations
Time-averaged data are used to investigate the effect of the plasma control on the main characteristic of the potential core and jet spreading. The results are presented in Figs. 6 and 7 where the statistics are averaged in time but also in the azimuthal direction. For the baseline case the length of the potential core is 7.2D with a jet diameter equal to 1.28D at the end of the computational domain. This potential core length is rather high compared to the literature on turbulent jets. This is due to the fact that the boundary layer inside the jet is laminar. The spreading angle of the jet is approximately 10°, after the transition period with no shear layer growth at the exit of the nozzle. The main effect of the tripping (case TR) is a significant decrease in length for the potential core by ≈ 25% when compared to the baseline case and an increase of the diameter of the jet by ≈ 10% at the end of the computational domain. The disturbances introduced in the nozzle by the random tripping are exciting some instabilities inside the nozzle boundary layer causing the turbulence breakdown to occur earlier for the jet. When comparing with the baseline case, increasing the intensity of the continuous forcing for the controlled simulations (I4, I8 and I12 cases) decreases the length of the potential core from 7.2D to 6.3D while at the same time increases the radial size of the jet with a diameter evolving from 1.31D to 1.56D at the end of the computational domain. The length of the potential core is not decreased when the plasma actuator cases are compared with the TR case. It seems that the random perturbations introduced in the nozzle boundary layer in the TR case are also triggering instabilities related to the jet column mode (as well as other modes). However, the potential core is thinner for the pulsating control (M0 case) even if the potential core length is the same as the one obtained in the TR case. These conclusions are consistent with the Fig. 4 where it can be seen that the turbulence levels are quite important in the near field for the controlled cases. As a result, the growth of the diameter is more important and the potential core is shorter. The time-averaged diameter of the jet for the simulation with a pulsating control (M0 case) is similar to the one obtained for case I12 however with a reduced length for the potential core, (5.3D same as case TR). It can be hinted that the large diameter observed for the pulsating case is a consequence of the generation of large puffs of turbulence at the exit of the nozzle as reported in Fig. 4. A reduced length for the potential core and a larger radial size for the jet are indicators of mixing enhancement. It suggests that the simulation with a pulsating control could be a good candidate for mixing enhancement.
Another important feature of the M0 case is the time-averaged radial shape of the jet with a strong initial spreading angle associated Fig. 9. 2D maps of the turbulent kinetic energy (TKE) at locations 0.5D, 1.0D, 2.0D, 4.0D and 6.0D after the nozzle exit.
V. Ioannou, S. Laizet International Journal of Heat and Fluid Flow 70 (2018) 193-205 with a rapid decrease of the diameter of the potential core. This can be attributed to an early transition to turbulence and a very strong intense small-scale activity materialised with the thin elongated structures around the large rings which seems to be able to entrain more surrounding fluid into the jet. It suggests that the flow dynamics for the pulsating case and for the tripping case are different, even if the extent of the potential core is the same. Basically, the jet is becoming fairly large very quickly for the M0 case and it remains large until the end of the computational domain whereas a continuous increase of the jet's radial size is observed for the tripping case. It might suggest that the turbulence puffs are a good mechanism for fluid entrainment and to generate high level of small-scale turbulence activities, even far downstream of the nozzle exit.

First and second order moments
The cross-sectional contours of the time-averaged streamwise velocity at locations 0.5D, 1.0D, 2.0D, 4.0D and 6.0D from the nozzle exit for all the cases are shown in Fig. 8. The imprint of the plasma actuators on the main jet can be directly observed on the mean streamwise velocity with clear distortions visible already 0.5D from the nozzle exit. The instabilities triggered in the nozzle boundary layer for the TR case do not significantly alter the azimuthal symmetry of the flow. The increased jet expansion relative to the baseline case can be confirmed as already suggested in Fig. 5. Very close to the nozzle exit, a pattern with the imprint of the thin elongated vortical structures at the edge of the jet can be seen for the controlled cases, corresponding to the number of plasma actuators. Actually the green imprints can be linked to the regions between two plasma actuators corresponding to the slowly moving elongated streamwise vortical structures observed in Fig. 5. Further downstream, a relaxation trend can be observed very quickly toward more conventional azimuthally quasi-homogeneous jets, except for the pulsating case for which the imprint of the forcing can be felt further downstream with very long thin streamwise vortices at distances 1D to 4D downstream of the nozzle exit located at the edge of the jet. It is interesting to note that for the continuous forcing (I4, I8 and I12 cases), the pattern of the jet at a distance 6D from the nozzle is strongly related to the intensity of the forcing with a circular, diamond and cross pattern observed for the I4, I8 and I12 cases, respectively. For the pulsating case, the pattern is a combination of a large cross and a large diamond.
The previous observations based on the mean velocity field gave us an idea on how the main jet is affected by the plasma actuators. However, it did not provide information about the production of unstable mechanisms resulting from the interaction between the plasma actuators and the main jet. In order to quantify those mechanisms, 2D maps of the turbulent kinetic energy are presented in Fig. 9. It can be seen that the plasma actuators give rise to very high levels of turbulent kinetic energy in the near-nozzle region, especially for the pulsating control case for which the imprint of the elongated structures at the edge of the jet can be observed up to 2D downstream of the nozzle exit. Denoted by a large and intense red region, the interactions between the plasma actuators and the main jet for the pulsating case must be very intense at the edge of the jet through the production of highly energetic instabilities. Overall, the peak intensity region for the turbulent kinetic energy is shifted toward the nozzle for the controlled cases by comparison to the baseline case but the turbulent kinetic energy levels seems to be higher for the baseline case at distances 1D and 2D from the nozzle exit. This is mainly due to an earlier breakdown of the jet into small scales for the controlled cases (as seen in Fig. 4) which means that the turbulent kinetic energy peaks earlier. Once again, the shape of the turbulent kinetic energy in Fig. 9 is strongly dependent of the intensity of the forcing. V. Ioannou, S. Laizet International Journal of Heat and Fluid Flow 70 (2018) 193-205 7. Mixing properties The effects of the control on the mixing properties of the main jet are considered by solving a passive scalar equation. The choice of a relevant criterion for evaluating the mixing efficiency is a rather delicate task. Following Lardeau et al. (2002), we consider the mixing of a passive scalar whose high concentration values of 1 correspond to the main jet inside the nozzle and zero elsewhere. The target for a wellmixed jet is to have an homogeneous scalar field with a mean value as low as possible between 0 and 1. The greater the scalar field expansion in the surrounding fluid, the higher the mixing improvement. In addition, the fluctuations of the passive scalar should be as low as possible for a good mixing. Measurement of mixing efficiency is also quantified by considering the statistical distribution of the passive scalar values on their interval of validity [0, 1].
Similar to the 2D maps for the time-averaged streamwise velocity and turbulence kinetic energy, Figs. 10 and 11 show 2D maps of the mean and fluctuating components of the scalar field at various locations downstream of the nozzle exit. The imprint of the control can be felt for the full computational domain with a visible signature of the braids at the edge of the jet and very distinctive patterns and different shapes for the mean component of the scalar field, depending on the plasma control configuration. The random disturbances inside the nozzle for the TR case and the plasma control actuators generate high values of fluctuations in the shear layer of the jet. As the jet develops, the fluctuations eventually reduce in intensity in the shear layer but expand radially inwards and outwards, phenomenon characterized by the green areas in Fig. 11. It can also be seen that for the baseline case and for the I4 case, a red/yellow ring in the center of the jet can be observed up to a distance 6D from the nozzle exit, suggesting that the scalar fluctuations are still important and not radially homogeneous. It is not an indicator of good mixing. As far as the mean component of the scalar field is concerned, a red region in the center of the main jet can be seen for all cases in Fig. 10 except for the pulsating case. Furthermore, the extent of the green area for that particular case is much larger, suggesting that the scalar field is more homogeneous, a good sign of mixing enhancement. The strong fluctuations in the near field (0.5D and 1D from the nozzle exit) allow a fast expansion of the scalar field outwards of the jet with plenty of time to homogenize.
In order to better quantify the effect of the plasma control on the mixing properties of the jet, probability density functions (PDFs) of the passive scalar obtained for each case at various streamwise sections are presented in Figs. 12 and 13 at the centreline and lipline of the jet. As a reminder, 0 and 1, correspond to ambient and jet fluid respectively. On the centreline at 4.5D from the nozzle exit, the M0 case has a high probability around two values (0.5 and 1). This is more likely due to the 50% duty cycle of the pulsating forcing. Only the I12 case also exhibits two values at this location but the first one (0.85) is much higher than the first one of the M0 case. For all the other, the high probability for the scalar field is for 1 only, similar to the value at the nozzle exit, suggesting that the mixing activity is very low. At a distance 7D downstream of the nozzle exit, all cases are showing a double peak for the PDFs. The M0 case has a very high probability at a value close to 0.5, followed by the TR case with a high probability for a value close to 0.65. All the other cases have their first peak for a value close to 0.8 with highest probabilities when the intensity of the forcing for the plasma actuator is important. For all the cases, the highest probability (more than 0.8) is for the second peak for a value of 1, except for the M0 case for which the probability for the value 1 is only 0.4.
The PDF data generated at the lipline are also very interesting to discuss in line with the mixing properties of the flow. The most important results are related to the size of the Gaussian-like shape for the Fig. 11. 2D maps of the fluctuating component of the scalar field at locations 0.5D, 1.0D, 2.0D, 4.0D and 6.0D after the nozzle exit.
V. Ioannou, S. Laizet International Journal of Heat and Fluid Flow 70 (2018) 193-205 PDFs and its peak value. In our set-up, good mixing can be characterized by a low peak value combined with a narrow shape for the probability function and the same value for the highest probability at the centreline and at the lipline. As hinted by the previous results, the sharpest Gaussian-like shape and the lowest peak are obtained for the M0 case. Furthermore, the M0 case is the only case for which the highest probability value is the same on the centreline and at the lipline at 4.5D and 7D from the nozzle exit, suggesting a very good homogeneity for the scalar field. We can therefore conclude that the pulsating control case is showing a great potential for mixing enhancement by comparison to the other cases.

Conclusions
The effect of four different control solutions based on eight Dielectric Barrier Discharge (DBD) plasma actuators located just before the nozzle exit of a turbulent jet were investigated with the aim to enhance the mixing of the jet. One of the controlled cases is based on a pulsating motion for the forcing at a frequency corresponding to the jet preferred frequency. Data were compared with two non-controlled cases, one with no perturbation inside the nozzle and one with a random tripping to trigger instabilities in the nozzle boundary layer.
The first important result is that the plasma actuators are able to strongly modify the flow field downstream of the nozzle with more or less the same flow rate as the non-controlled cases. The effect of the plasma actuators can easily be seen with ejections of pairs of elongated streamwise vortical structures generated between two plasma actuators. The breakdown to turbulence is therefore happening closer to the nozzle exit by comparison to the baseline case due to the promotion of strong thin elongated structures around the large ring generated at the nozzle of the jet. As a consequence, a reduced length for the potential core and an increase of the size of the jet were observed for the control cases by comparison to the baseline case. The reduction in length of the potential core is not present when comparisons are made with the TR case. However, the shape of the potential core is affected by the plasma actuators with a thinner potential core. It seems to suggest that the shape of the potential core can influence the mixing properties of the jet. It should also be noted that each controlled case is producing a different pattern for the radial shape of the jet which seems to suggest that the intensity of the forcing is an important parameter for control purposes.
The passive scalar study revealed that the pulsating controlled case can enhance mixing by comparison to the non-controlled cases, with a more homogeneous scalar field and low levels of scalar fluctuations. These conclusions were confirmed with PDFs of the scalar field downstream of the lipline with a sharp Gaussian-like shape and the same low  V. Ioannou, S. Laizet International Journal of Heat and Fluid Flow 70 (2018) 193-205 peak value by comparison to the other cases on the centreline and at the lipline. Future studies will investigate the influence of the number of plasma actuators inside the nozzle as well as their location with respect to the nozzle exit. Since the pulsating motion was quite successful, various duty cycles will be investigated at various Reynolds numbers. Different excitation modes will be tested, like the first and second helical modes, following the work of Samimy et al. (2004Samimy et al. ( , 2007a, Kim et al. (2009), Samimy et al. (2012 for supersonic turbulent jets, as well as investigating further the jet column mode in order to find the most effective way to improve the mixing properties of the jet. Finally, the potential of forcing azimuthal modes will be explored in great detail, as they can have an important contribution for triggering instabilities (Glauser et al., 1987;Raman et al., 1994) with the potential to affect the mixing and acoustic properties of the jet.