1. Introduction
A variable-density jet is a flow induced by the continuous release of a fluid of density 
$\rho _s$ from a circular source, of radius 
$r_s$, at a steady uniform velocity 
$w_s$ within a still environment of density 
$\rho _0$. For sufficiently high Reynolds number 
$Re_s={w_s r_s}/{\nu _s}$ (
$\nu _s$ is the kinematic viscosity of the fluid discharged) the flow is fully turbulent, and the growth of inertial instabilities gives rise to complex and multiscale dynamics (see figure 1). The process leads to the entrainment of ambient fluid within the jet that ultimately results in the mixing and dilution of the discharged fluid. These releases are referred to as jets (instead of plumes or fountains) because their dynamics is dominated by momentum (rather than buoyancy), which implies that their Richardson number at the source 
$Ri_s = {g r_s (\rho _s - \rho _0) / \rho _0 }/{w_s^2} \rightarrow 0$. We can therefore assume that gravitational effects are negligible in the flow dynamics.
Instantaneous density ratio 
$\rho / \rho_0$ of a turbulent jet with 
$\rho _s / \rho _0 = 0.2$ (large-eddy simulation): white corresponds to 
$\rho / \rho _0 =1$ while black to 
$\rho / \rho _0 = 0.2$. Superimposed are the averaged value of the streamwise velocity (in red), the top-hat velocity 
$w_m$ and the radius 
$r_m$ at varying axial distance 
$z$ from the source.
Examples of variable-density turbulent jets (and plumes) can be found in several areas of engineering and the physical sciences, e.g. ocean outfalls, exhausts of engines, volcanic eruptions, releases of chemicals from pipes and reservoirs. A robust approach in the modelling of these flows is provided by the plume equations (see e.g. Priestley & Ball Reference Priestley and Ball1955; Morton, Taylor & Turner Reference Morton, Taylor and Turner1956), which provide an integral description of the flow, in which spatially averaged variables depend only on the distance from the source 
$z$ (figure 1). In this framework, the integration of the equation of conservation of (time-averaged) mass over the radial coordinate 
$r$ (see figure 1) leads to
\begin{equation} \frac{\mathrm{d} G}{\mathrm{d} z} = 2 \rho_0 \alpha w_m r_m, \end{equation}
where 
$G \equiv 2\int _{0}^{\infty } \overline {\rho w}r\,\mathrm {d} r$ is the time-averaged mass flux (
$w$ is the jet velocity along the streamwise direction 
$z$, 
$\rho$ is the jet density and overbar denotes Reynolds averages), 
$\alpha$ is the entrainment coefficient and 
$r_m$ and 
$w_m$ are characteristic scales of the jet width and velocity, respectively (their definition is given in § 2).
Our focus is the parameter 
$\alpha$, controlling the dilution of the jet with the ambient fluid. A reference study providing experimental measurements of the entrainment process within variable-density releases is that of Ricou & Spalding (Reference Ricou and Spalding1961), hereafter referred to as R&S. Their results suggest a dependence of 
$\alpha$ on a local jet density scale 
$\rho _m (z)$ (also defined in § 2), and a scaling of the form
\begin{equation} {\alpha = \alpha_0 \left( \frac{\rho_m}{\rho_0} \right)^{\beta}} ,\end{equation}
with 
$\beta =1/2$ and 
$\alpha _0$ the entrainment coefficient for a pure (iso-density) jet, which they set equal to 0.08, which is slightly larger than that estimated in other studies (Papanicolaou & List Reference Papanicolaou and List1988; Wang & Law Reference Wang and Law2002; Ezzamel, Salizzoni & Hunt Reference Ezzamel, Salizzoni and Hunt2015; van Reeuwijk et al. Reference van Reeuwijk, Salizzoni, Hunt and Craske2016). According to (1.2), the rate of dilution of the jet is therefore increased for an increasing density ratio 
${{\rho _m}/{\rho _0}}$. This result was subsequently interpreted physically by Morton (Reference Morton1965) in terms of the influence of the Reynolds stress (i.e. of the turbulent kinetic energy (TKE) production) on the entrainment process. Almost forty years after the seminal work by R&S, Rooney & Linden (Reference Rooney and Linden1996) and Woods (Reference Woods1997) adopted an entrainment parametrisation on the basis of (1.2) with 
$\beta =1/2$, showing that the dynamics of a so-called non-Boussinesq plume (i.e. characterised by large density variations between plume and ambient fluid) could be ultimately modelled in analogy with the classical plume theory. This led to the adoption of self-similar solutions for non-Boussinesq plumes (van den Bremer & Hunt Reference van den Bremer and Hunt2010), that could be subsequently applied to releases in unconfined and confined environments (Vauquelin Reference Vauquelin2015; Jiang et al. Reference Jiang, Creyssels, Hunt and Salizzoni2019; Mehaddi et al. Reference Mehaddi, Boulet, Koutaiba, Vauquelin and Candelier2021).
R&S's variable-density jets have been studied using experiments (e.g. Panchapakesan & Lumley Reference Panchapakesan and Lumley1993; Gharbi, Amielh & Anselmet Reference Gharbi, Amielh and Anselmet1995; Amielh et al. Reference Amielh, Djeridane, Anselmet and Fulachier1996; Djeridane et al. Reference Djeridane, Amielh, Anselmet and Fulachier1996; Pietri, Amielh & Anselmet Reference Pietri, Amielh and Anselmet2000; O'Hern et al. Reference O'Hern, Weckman, Gerhart, Tieszen and Schefer2005; Ai, Law & Yu Reference Ai, Law and Yu2006; Gerashchenko & Prestridge Reference Gerashchenko and Prestridge2015; Charonko & Prestridge Reference Charonko and Prestridge2017; Viggiano et al. Reference Viggiano, Dib, Ali, Mastin, Cal and Solovitz2018) and numerical simulations (e.g. Gharbi et al. Reference Gharbi, Ruffin, Anselmet and Schiestel1996; García-Villalba, Fröhlich & Rodi Reference García-Villalba, Fröhlich and Rodi2006; Wang et al. Reference Wang, Frohlich, Michelassi and Rodi2008). However, the determination of the dependence of the entrainment coefficient on the density ratio was beyond the scope of those studies. None, therefore, could confirm (or invalidate) the scaling provided by (1.2).
The issue was tackled recently by Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023), who estimated velocity and density statistics within a low-density jet (
$\rho _m / \rho _0 = 0.4$) and an iso-density jet by means of laboratory experiments and large-eddy simulations. Their results highlight the role of a reduced density in the jet dynamics, notably by enhancing the TKE production and modifying the shape of the mean streamwise velocity. Despite this, their results clearly indicate that a reduced density has no major effect in modifying the entrainment coefficient 
$\alpha$, whose evolution for increasing distance from the source does not differ significantly from that observed in an iso-density jet (therefore contradicting the R&S scaling). The aim of this study is to build upon the recent results by Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023) and extend their analysis to a broader set of source density conditions (
$\rho _s$), carrying out an exhaustive investigation into the dependence of flow dynamics and the entrainment process in variable-density jets.
We begin by revisiting the scaling (1.2) proposed by R&S, showing analytically that their experimental results are not conclusive in determining the value of the exponent 
$\beta$. We then present the details of the large-eddy simulations (§ 4) of jets with varying density ratios (spanning over two orders of magnitude) and the entrainment decomposition (§ 3), which we subsequently use to interpret the numerical results. Next, we analyse the evolution of integral quantities characterising the entrainment dynamics, integrating in our discussion the evidence provided by recent experimental data on light (Salizzoni et al. Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023) and heavy (Charonko & Prestridge Reference Charonko and Prestridge2017) jets, and ultimately investigate the dependence of 
$\alpha$ on the density ratio (§ 5).
2. Revisiting the experimental results of R&S
By employing a meticulous experimental approach, R&S successfully estimated the entrained ambient air flux within a variable-density jet. This estimation was achieved using a porous screen positioned at the jet border, therefore avoiding challenging (especially at that time) measurements of local fluid velocities and densities within the jet. Their conclusion was that the mass fluxes 
$G$ of the different releases, normalised by their source value 
$G_s$, collapsed when plotted against vertical coordinate rescaled on 
${( {\rho _s}/{\rho _0} )^{-1/2}}$ (see figure 2c,d). Coupled with the assumption of a constant momentum flux, corresponding to negligible gravitational effects on the flow dynamics, their results led to the conclusion that 
$\beta = 1/2$ in (1.2).
Results of analytical solutions for 
$z/r_s \leq 30$ (a,b) and for 
$z/r_s \leq 900$ (c,d). Dotted lines: R&S model, 
$\beta =1/2$, (2.4) for 
$\rho _s/\rho _0 = 0.069$ (in red) and for 
$\rho _s/\rho _0 = 1.5$ (in black). Solid green line: (2.4) with 
$\rho _s = \rho _0$. Dashed lines: solution of the differential equation (2.3) with 
$\alpha = \alpha _0 = 0.08$, the value originally used by R&S (see text), and 
$\beta = 0$. Experimental results of Ricou & Spalding (Reference Ricou and Spalding1961) are also shown in (c,d) (symbols).
In this section we demonstrate that the experiments of R&S do not provide sufficient information to determine the exponent 
$\beta$ in (1.2) because, to leading order, the value of 
$\beta$ does not affect the scaling of 
$G$ in the far field.
To show this, we solve (1.1) analytically, assuming that 
$\alpha$ scales as (1.2) without specifying 
$\beta$. We then compare solutions for two values of 
$\beta$ in the near field (
$z/r_s \leq 30$; figure 2a,b), where variable-density effects are large, and in the far field (
$z/r_s$ up to 900). The latter range was investigated by R&S (figure 2c,d), which allows for a direct comparison between the analytical solution and their data.
In what follows, the characteristic velocity, radius and density of the jet are defined as 
$w_m = M/G$, 
$r_m = (Q G/M)^{1/2}$ and 
$\rho _m = G/Q$, respectively, as a function of the mass 
$G$, volume 
$Q$ and momentum 
$M$ fluxes (Salizzoni et al. Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023). Adopting Favre averages (denoted by a tilde and defined as 
$\tilde {\chi }=\overline {\chi \rho }/ \bar {\rho }$) the fluxes are expressed as 
$G = 2\int _{0}^{\infty } \bar {\rho } \tilde {w}r\,\mathrm {d} r = \rho _m w_m r_m^2$, 
$Q = 2\int _{0}^{\infty }\tilde {w}r\,\mathrm {d} r = w_m r_m^2$ and 
$M = 2\int _{0}^{\infty }\bar {\rho }\tilde {w}^{2}r\,\mathrm {d} r = \rho _m (w_m r_m)^2$.
2.1. Analytical solutions for the mass flux 
$G$
Since we are assuming negligible buoyancy effects, which means that the momentum flux is conserved, i.e. 
$M \equiv \rho _m w_m^2 r_m^2 = \rho _s w_s^2 r_s^2$ (with 
$r_s = r_m(z=0)$, 
$w_s = w_m(z=0)$ and 
$\rho _s = \rho _m(z=0)$), the mass balance equation (1.1) can be rewritten as
\begin{equation} \frac{\mathrm{d} (G/G_s)}{\mathrm{d} (z/r_s)} =2 \alpha \left ( \frac{\rho_m}{\rho_0}\right )^{{-}1/2} \left (\frac{\rho_s}{\rho_0} \right )^{{-}1/2}. \end{equation}
By further assuming a low-Mach-number flow at high Reynolds number (i.e. negligible viscous effects) it can be shown (Rooney & Linden Reference Rooney and Linden1996) that the flow has a divergence-free velocity field, which, in turn, implies that the density deficit flux is conserved, such that 
$(\rho _0 -\rho _m) w_m r_m^2 = (\rho _0 -\rho _s) w_s r_s^2$, or, expressed in terms of 
$G$,
\begin{equation} \frac{\rho_0}{\rho_m} = 1+ \varDelta_s \frac{G_s}{G}, \end{equation}
where 
$\varDelta _s = (\rho _0-\rho _s)/\rho _s$ and 
$G_s= \rho _s w_s r_s^2$. Combining (2.1) and (2.2), and assuming that 
$\alpha$ exhibits a power-law dependence of the form (1.2), leads to
\begin{equation} \frac{\mathrm{d} (G/G_s) }{\mathrm{d} (z/r_s)} = 2 \alpha_0 \left (\frac{\rho_s}{\rho_0} \right )^{{-}1/2} \left [ 1 + \varDelta_s \frac{G_s}{G} \right]^{1/2-\beta}. \end{equation} For 
$\beta = 1/2$, integrating (2.3) gives
\begin{equation} \frac{G}{G_s} = 1+ 2 \alpha_0 \left (\frac{\rho_s}{\rho_0} \right )^{{-}1/2} \frac{z}{r_s},\quad \mbox{for } \beta = 1/2, \end{equation}
which is the reference relation used by R&S to rescale their experimental results and demonstrate collapse onto a single non-dimensional curve. Indeed, (2.4) exhibits a dependence on the source conditions (
$r_s$ and 
$\rho _s$) only, and not on the local jet radius 
$r_m$, density 
$\rho _m$ or velocity 
$w_m$, which R&S were not able to measure using their experimental apparatus.
When considering the general case of 
$\beta \neq 1/2$, 
$G/G_s$ is not linear in 
$z/r_s$ and therefore does not scale with 
$({\rho _s} / {\rho _0} )^{-1/2}$. However, we show hereafter that, irrespective of the 
$\beta$ value, 
$G/G_s$ behaves as 
$({\rho _s} / {\rho _0} )^{-1/2} {z}/{r_s}$ at large distances from the source, i.e. 
$z \gg r_s$, corresponding to the limit 
$G_s / G \rightarrow 0$. By defining 
$\epsilon = G_s / G$, the differential equation for 
$\epsilon$ is
\begin{equation} \frac{\mathrm{d} \epsilon}{\mathrm{d} (z/r_s)} + 2 \alpha_0 \left (\frac{\rho_s}{\rho_0} \right )^{{-}1/2} \epsilon^2 \left [ 1 + \varDelta_s \epsilon \right]^{-\gamma} = 0, \end{equation}
with 
$\gamma = \beta - 1/2$. Separating 
$\epsilon$ and 
$z/r_s$ and integrating leads to
\begin{equation} \int_{\epsilon}^1 \left [ 1 + \varDelta_s x \right]^{\gamma} \frac{\mathrm{d}\kern0.07em x}{x^2} = 2 \alpha_0 \left (\frac{\rho_s}{\rho_0} \right )^{{-}1/2} \frac{z}{r_s}. \end{equation}Formally, the integrand of the left-hand side of (2.6) can be expanded as a power series which can be integrated term by term:
\begin{align} \int_{\epsilon}^1
\left[ 1 + \varDelta_s x \right]^{\gamma} \frac{\mathrm{d}\kern0.07em
x}{x^2} &= \int_{\epsilon}^1 \left[ 1 + \sum_{n=1}^\infty
\binom{\gamma}{n} \left ( \varDelta_s x \right )^n \right] \frac{\mathrm{d}\kern0.07em x}{x^2}
\end{align}
\begin{align} &= \frac{1}{\epsilon} -1 - \gamma \varDelta_s \ln \epsilon + \sum_{n=2}^\infty \binom{\gamma}{n} \frac{\left ( \varDelta_s \right )^n}{ n-1} \left ( 1 - \epsilon^{n-1} \right ). \end{align} At large distances from the source, 
$\epsilon = G_s / G \rightarrow 0$ the leading-order contribution from (2.8) is 
$1 / \epsilon = G/G_s$ which from (2.6) implies that
\begin{equation} \frac{G}{G_s} \sim 2 \alpha_0 \left (\frac{\rho_s}{\rho_0} \right )^{{-}1/2} \frac{z}{r_s},\quad \mbox{for any } \beta \text{ with } G_s / G \rightarrow 0.\end{equation} This limit can be easily confounded with the exact solution for 
$\beta =1/2$ (equation (2.4)). In other words, when performing measurements in the far field where 
$G / G_s \gg 1$, we expect 
${G}/ {G_s}$ to scale with 
$({\rho _s} / {\rho _0} )^{-1/2}$ regardless of the value of the exponent 
$\beta$.
2.2. Near- and far-field behaviour and comparison with R&S data
In order to analyse the dependence of the mass flux 
$G$ on the exponent 
$\beta$, we focus on the near-field region, i.e. 
$z/ r_s \leq 30$ (for which R&S did not obtain data), and compare the R&S model (
$\beta =1/2$, (2.4)) and solution of the differential equation (2.3) with 
$\beta =0$. The results are plotted in figure 2(a) for the three density ratios investigated by R&S: 
$\rho _s / \rho _0 = 1.5$ (carbon dioxide in air), 
$\rho _s / \rho _0 = 1$ (air in air) and 
$\rho _s / \rho _0 = 0.069$ (hydrogen in air).
The results, plotted in figure 2(a), show that, for density ratios slightly different from unity (i.e. carbon dioxide release), the value of 
$\beta$ has very little effect on the mass flux 
$G(z)/G_s$. Significant differences can instead be seen for the case of hydrogen, for which the two solutions (
$\beta =0$ and 
$\beta =1/2$) lead to differences of almost 50 % after a few source radii from the source. Consequently, the different solutions do not collapse on a single curve, when re-scaling the data as in R&S (figure 2b).
Turning to the far field, and including R&S data in our analysis, we obtain a different picture (see figure 2c,d). As previously pointed out, the solutions obtained for the two different 
$\beta$ are almost indistinguishable when plotted over large distances from the source, i.e. up to 
$z / r_s = 900$. Therefore, very good agreement is observed between the experimental data and both analytical solutions, one obtained using the R&S model (
$\beta = 1/2$, dotted lines) and the other using an entrainment coefficient that is independent of the density ratio (
$\beta =0$, dashed lines). Both solutions collapse, as a first approximation, on a single curve. Indeed, the discrepancies between the two solutions in figure 2(d) are lower than the uncertainties observed in the experimental data.
We conclude that the solution of the differential equation (2.3) is sensitive to the value of 
$\beta$ in the near field but not in the far field. Therefore, the collapse of the experimental data on a single curve in the far field according to (2.9), as observed by R&S, does not imply that the entrainment coefficient scales as (1.2), with 
$\beta =1/2$.
3. The entrainment relation
Based on early work of Priestley & Ball (Reference Priestley and Ball1955), recent studies (Kaminski, Tait & Carazzo Reference Kaminski, Tait and Carazzo2005; Ezzamel et al. Reference Ezzamel, Salizzoni and Hunt2015; van Reeuwijk & Craske Reference van Reeuwijk and Craske2015; Craske, Salizzoni & van Reeuwijk Reference Craske, Salizzoni and van Reeuwijk2017) on iso-density jets and Boussinesq plumes have shown that, by combining the mass, momentum and mean kinetic energy balance equations, a relation for the entrainment coefficient, expressed as a function of local flow statistics, can be obtained. As discussed by Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023), this framework can be conveniently extended to the more general case of large density difference between an ambient fluid and the release, by using Favre averages.
To leading order, the entrainment relation for a momentum-dominated release (in which the role of gravity is negligible) can be written as (Salizzoni et al. Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023)
\begin{equation} \alpha_E = \underbrace{- \frac{\rho_m}{\rho_0} \frac{{\delta}_m}{2{\gamma}_m}}_{\alpha_{prod}} + \underbrace{\frac{\rho_m}{\rho_0} r_m \frac{\text{d} }{\text{d} z} \left( \ln {{\gamma}_m^{1/2}} \right) }_{\alpha_{shape}} ,\end{equation}where
\begin{equation} {\gamma}_m = \frac{2}{\rho_m w_m^3 r_m^2} \int_0^\infty \bar{\rho}\tilde{w}^3 r \,\text{d} r, \quad {\delta}_m = \frac{4}{\rho_m w_m^3 r_m} \int_0^\infty \bar{\rho}\widetilde{w''u''} \frac{\partial\bar{w}}{\partial r} r\, \text{d} r, \end{equation}
and double primes denote fluctuations from the Favre averages. The parameters 
$\gamma _{m}$ and 
$\delta _{m}$ correspond to the dimensionless mean flux and production of kinetic energy, respectively, and are referred to as ‘profile coefficients’ because their value depends on the radial dependence of the primitive flow variables.
As evidenced by (3.1), the entrainment coefficient is determined by two contributions (both weighted by the ratio between the top-hat jet density 
${\rho _m}$, varying with the streamwise coordinate, and the ambient density 
${\rho _0}$): (i) 
$\alpha _{prod}$, directly linked to the ratio of production of TKE to the mean kinetic energy flux, and (ii) 
$\alpha _{shape}$, which depends on the change of the form of the radial profiles of mean streamwise velocity. A priori, there is no reason not to assume a possible dependence of 
$\alpha _{prod}$ and 
$\alpha _{shape}$ on 
$\rho _m/\rho _0$, which is effectively taken as unity in the far field and therefore typically ignored.
Note that the estimate of the entrainment coefficient as given by (3.1) is independent of that provided by (1.2), denoted here as 
$\alpha$, and linked to the rate of change of the mass flux 
$G(z)$ for increasing distance from the source. However, to fulfil basic mass, momentum and kinetic energy balances, the coefficients 
$\alpha$ and 
$\alpha _E$ have to coincide at any distance from the source.
As also revealed by (3.1), the evolution of the entrainment coefficient can be captured by analysing three main variables of the flow: the mean streamwise velocity 
$\tilde {w}$, the Reynolds stress 
$\widetilde {w''u''}$ and the density 
$\bar {\rho }$. The determination of these flow statistics by means of laboratory experiments faces practical and metrological difficulties (as discussed in Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023)), which can be conveniently circumvented by performing numerical simulations.
4. Numerical simulations
Numerical simulations of variable-density jets in a quiescent environment were performed using the numerical computational code CALIF
$^3$S (developed at the Institut de Radioprotection et de Sûreté Nucléaire-IRSN), dedicated to three-dimensional simulations of turbulent and slightly compressible flows (low-Mach-number approach).
The code uses large-eddy simulation to solve the Favre-filtered Navier–Stokes equations (mass and momentum balance), along with species transport equations, in Cartesian coordinates. The subgrid-scale Reynolds stresses are modelled by the Vreman subgrid-scale model (Vreman Reference Vreman2004), and a gradient diffusion hypothesis (SGDH) is used to close the problem with a turbulent Schmidt number 
$Sc_t$ set to 
$0.7$. As molecular viscosity is disregarded in this formalism, the Reynolds number is undefined. The code has been previously used to study flows characterised by large density differences such as light jets (Salizzoni et al. Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023) and heavy and light fountains (Vaux et al. Reference Vaux, Mehaddi, Vauquelin and Candelier2019).
The domain is a cube of size 
$40 r_s$, with 
$r_s$ the jet radius at the source. In the central part of the domain (a rectangular cuboid with square section of side length 
$10 r_s$ centred on the origin) a refined grid size of 
$r_s/14$ is used. The grid is then progressively stretched for increasing distance from the jet axis, as in Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023). In the vertical direction the grid size is equal to 
$r_s/8$. For the time discretisation, a Courant–Friedrichs–Lewy (CFL) number close to unity has been imposed for each calculation even though time step sizes for which CFL numbers greater than one are allowed with the use of implicit schemes. Each simulation lasted 
$1000 T$ where 
$T= r_s / w_s$. Results for the first 
$400 T$ were discarded and flow statistics were computed over an interval of 
$600 T$. At the source, the inlet flow was perturbed with the method presented by Jarrin et al. (Reference Jarrin, Benhamadouche, Laurence and Prosser2006) in order to trigger the transition to turbulence. All numerical simulations have been performed with the same computational settings (grid and duration), spatially uniform source conditions and inlet perturbation, except for the density ratio 
$\rho _s/\rho _0$.
The source is set flush with the bottom solid boundary of the computational domain from which the jet emerges at the centre with a top-hat inlet mean velocity profile. Since we consider jets in an infinite (open) environment, the computational domain must be bounded by artificial boundary conditions which perturb as little as possible the flow in the interior of the domain. The boundary conditions used in our simulations are based on the control of kinetic energy, allowing one to distinguish between the flow that leaves the domain and the flow that enters it (Bruneau & Fabrie Reference Bruneau and Fabrie1994, Reference Bruneau and Fabrie1996; Bruneau Reference Bruneau2000). We have performed eight numerical simulations of turbulent jets (three light releases, four heavy releases and an iso-density jet) with a density ratio 
$\rho _s / \rho _0$ spanning over almost two orders of magnitude, i.e. 
$0.2 \leq \rho _s / \rho _0 \leq 10$. Our experience is that spatially uniform source conditions naturally promote the decay of the potential core, although the difference we observed in supplementary simulations of light jets using spatially non-uniform source conditions did not significantly affect the results.
In order to avoid any effects of variable density associated with the gravitational field (i.e. buoyancy effects), we considered a situation without gravity, by imposing a null gravitational acceleration in the momentum budget. Despite this, we will still refer to as ‘heavy’ the high-density releases (
$\rho _s > \rho _0$) and the low-density releases (
$\rho _s < \rho _0$) as ‘light’.
A detailed verification of the reliability of the code in reproducing the first- and second-order flow statistics was performed in Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023), who compared the large-eddy simulation with experimental results for an air–helium jet with initial density ratio 
$\rho _s / \rho _0 = 0.4$ and for an iso-density jet. Their study included an analysis of the effect on the flow statistics of slightly different boundary conditions as well as different subgrid models. It was shown that the results were not sensitive to the subgrid model and that boundary conditions could impact local flow statistics near the source, but do not significantly affect the integral flow statistics.
5. Results
5.1. Local and integral flow statistics
An overview of the results of the eight simulations is provided in figure 3. Here we show the vertical evolution of integral flow variables as a function of the distance from the source. Figure 3(a) presents the variation of the density ratio 
$\rho _m / \rho _0$ due to the progressive dilution with the ambient fluid: its decay for heavy releases and its increase for the light ones. This dilution process is accompanied by an enhancement of the jet radius 
$r_m$ (figure 3b), a decay of the jet velocity 
$w_m$ (figure 3c) and an increase of the volume 
$Q$ and mass 
$G$ fluxes (figure 3d,e).
Density ratio (
$a$), jet width (
$b$), jet velocity (
$c$), volume (
$d$) and mass (
$e$) fluxes as a function of the streamwise 
$z$ coordinate. Eight jets have been simulated, from the lightest to the heaviest: 
${\rho _s} / {\rho _0}=$ 0.2 (
$\bigtriangleup$), 0.4 (
$*$), 0.66 (
$\times$), 1 (
$\circ$), 1.5 (
$\star$), 2.5 (
$\diamond$), 5 (
$\triangledown$), 10 (
$\square$).
Figure 3(e) demonstrates that the heavier the release, the lower the growth rate of the non-dimensional mass flux 
$G/G_s$. Conversely, the non-dimensional volume flux 
$Q/ Q_s$ is higher for heavier releases (figure 3d). Indeed, if the entrained fluid has a density which is much smaller than the jet density, it contributes little to the enhancement of the mass flux 
$G$. These heavy releases dilute more slowly than their lighter counterparts, as is evidenced by a reduced growth rate for the jet radius 
$r_m$ and a reduced decay of the jet velocity 
$w_m$.
As discussed in § 3, besides the evolution of the jet density 
$\rho _{m}$ and the mass flux 
$G$, the key elements underlying the turbulent entrainment process (3.1) are the evolutions of the radial profiles of the mean vertical velocity 
$\tilde {w}(r,z)$ and of the Reynolds stress 
$\widetilde {w''u''}(r,z)$. We therefore focus on these two flow variables, whose radial profiles are presented in figure 4, for three of the simulated jets: a light jet (
$\rho _s / \rho _0 = 0.2$), an iso-density jet (
$\rho _s / \rho _0 = 1$) and a heavy jet (
$\rho _s / \rho _0 = 5$). The light jet clearly exhibits a more rapid modification of the radial profiles of both 
$\tilde {w}$ and 
$\widetilde {w''u''}$, compared with the other two (figure4a,d). The uniform conditions imposed at the source are almost immediately lost, even in the core of the jet, for 
$z / r_s > 5$. Therefore, the extent of the potential core tends to be smaller as the density ratio is reduced. Analysing the case of the heavy jet, we can observe a flat velocity profile in the core of the jet that persists up to 
$z / r_s = 15$ (see figure 4a–c). Corresponding profiles of the Reynolds stress (see figure 4d–f) exhibit zero value in the core of the jet that persists much longer within heavy releases, and one peak that increases in magnitude for the light release.
Radial profiles of mean vertical velocity 
$\tilde {w}$ (a–c) and Reynolds stress 
$\widetilde {u''w''}$ (d–f) for a light jet with (
$\rho _s/\rho _0 = 0.2$), an iso-density jet and a heavy jet (
$\rho _s/\rho _0 = 5$). Symbols as in figure 3. From left to right: 
$z = 5r_s$, 
$15r_s$ and 
$30 r_s$.
As a consequence of the behaviour shown in figure 4, the profile coefficients 
$\gamma _m$ and 
$\delta _m$ exhibit significant variation (see figure 5). It is worth noting that, since the simulations are characterised by the same boundary conditions at the source, all variations observed in the profile coefficients 
$\gamma _m$ and 
$\delta _m$ are solely due to the varying density ratios. In the case of light releases, the dimensionless mean kinetic energy flux 
$\gamma _m$ increases sharply in the near field, exceeding the asymptotic value 
$\gamma _m =4/3$, which corresponds to a Gaussian radial velocity profile. Its peak is just below 
$z / r_s = 10$, after which it reduces in value. The maximum value reached by 
$\gamma _m$ depends on the density ratio at the source: the smaller the density ratio, the higher the value of the peak. Conversely, for heavy jets, moving away from the source, the value of 
$\gamma _m$ increases smoothly and tends progressively to the asymptotic value for an iso-density jet, which is far from being attained at the top of our computational domains for most of the cases analysed. The larger the density ratio at the source, the smoother the increase of 
$\gamma _m$ with distance from the source.
Vertical evolution of the (a) 
${\gamma }_m$ and (b) 
${\delta }_m$ profile coefficients. Symbols as in figure 3.
The behaviour of the non-dimensional TKE production term 
$- \delta _m$ is qualitatively similar to that of the dimensionless mean kinetic energy flux term 
$\gamma _m$. For light releases, we observe a rapid increase up to a peak and a subsequent decrease. Heavy releases exhibit instead a progressive enhancement towards the asymptotic value. Physically this means that light jets are characterised by a larger TKE production in the near field, when compared with iso-density jets, and that in heavy jets the TKE production is instead inhibited.
5.2. Entrainment
According to (3.1), the variations of 
$\gamma _m$ and 
$\delta _m$, observed for the jets with varying density ratios at the source, are expected to have a direct impact on the entrainment rate. Depending on the initial source conditions, these affect the entrainment coefficient 
$\alpha _E$ in different ways by inducing variations on both the production and the shape terms (
$\alpha _E = \alpha _{prod}+\alpha _{shape}$). Figure 6 shows the evolution with distance from the source of the entrainment coefficients 
$\alpha$, as defined by (1.1), and 
$\alpha _E$, 
$\alpha _{prod}$ and 
$\alpha _{shape}$, as defined by (3.1) for six (light and heavy) releases.
Vertical evolution of the entrainment coefficients 
$\alpha$ (solid black lines), 
$\alpha _E$ (filled symbols), 
$\alpha _{prod}$ (open symbols) and 
$\alpha _{shape}$ (dashed lines), for light jets (a,b), iso-density jet (c) and three heavy jets (d–f). Also shown: experimental data of Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023) represented by black diamonds in (b,c) and Charonko & Prestridge (Reference Charonko and Prestridge2017) (black pentagrams) in (c,e).
For all cases the values of the two coefficients 
$\alpha$ and 
$\alpha _E$ (represented by solid black lines and filled symbols, respectively, in figure 6a–d), which provide independent means of evaluating entrainment, are in good agreement. This means that, as previously observed for Boussinesq plumes (van Reeuwijk et al. Reference van Reeuwijk, Salizzoni, Hunt and Craske2016), neglecting higher-order terms is a reasonable approximation for the estimate of 
$\alpha _E$ (see Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023) for further details). The only differences between the two can be observed very close to the source for the two heaviest releases (
${{\rho _s}/{\rho _0}= 5 ; 10}$). Figure 6(b,c) shows good agreement with the experimental estimate of 
$\alpha$ by Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023), for the case 
$\rho _s / \rho _0 = 0.4$ and for the iso-density jet. Further insight on the numerical estimates of the 
$\alpha _{prod}$ and 
$\alpha _{shape}$ coefficients (which are not plotted here to limit the number of symbols on the graph) and the reliability of large-eddy simulation data in reproducing the flow statistics are provided in Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023).
In the case of light releases (
${{\rho _s}/{\rho _0}= 0.2 ; 0.4}$) most of the variations of both coefficients 
$\gamma _m$ and 
$\delta _m$ occur within 10 source radii (see figure 5), attaining peak values that exceed significantly those for an iso-density jet. Their influence on the production and shape terms of 
$\alpha _E$ is, however, weighted by the local value of the density ratio 
$\rho _m / \rho _0$ (see (3.1)). This implies that, even if 
$\gamma _m$ and 
$-\delta _m$ increase with 
$z$ more rapidly for a light jet than for an iso-density jet (see figure 5a–c), the evolution of 
$\alpha _{prod}$ is very similar. At the same time, 
$\alpha _{shape}$ varies very little and is close to 0 for the light jet further than 10
$r_s$, indicating no changes to the shape of the profiles. As a result, and as already evidenced by Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023), the evolution of the entrainment coefficient 
$\alpha _E=\alpha _{prod}+\alpha _{shape}$ for the light jet is not markedly different from that of an iso-density jet. In both cases of a light and iso-density jet, 
$\alpha _{prod}$ rapidly becomes the most significant contribution to the global entrainment.
The heavy jet cases, i.e. 
${{\rho _s}/{\rho _0}= 1.5 ; 5 ; 10}$, behave differently (see figure 6d–f). The contribution of the shape term 
$\alpha _{shape}$ decreases much more slowly and is not negligible within almost the whole computational domain. In the very near field (
$z \leq 5 r_s$), due to the high value of 
$\alpha _{shape}$, the entrainment coefficient is clearly enhanced compared with that of an iso-density jet. When 
$z\geq 5 r_s$, 
$\alpha$ and 
$\alpha _E$ follow the same trend as 
$\alpha _{shape}$, i.e. a slight decrease, then they reach a value close to 
$\alpha _0 \approx 0.07$, which is in close agreement with the estimate provided by direct numerical simulations of iso-density jets (van Reeuwijk et al. Reference van Reeuwijk, Salizzoni, Hunt and Craske2016). The production term for heavy jets is instead reduced in the very near field (
$z \leq 4 r_s$) due to the reduced TKE production.
It is clear from the analysis above that the evolution of the profile coefficients and entrainment contributions depend strongly on the density ratio and its distance to the source. In order to obtain a clearer picture of the effect of the density ratio, we remove the near-source region 
$z/r_s<15$ from the analysis, as this region (for the case of the heavy releases) is strongly affected by the breakdown of the potential core, which implies strong contributions to 
$\alpha _{shape}$. As can be seen in figure 6, this discards contributions of 
$\alpha _{shape}$ in all except the heaviest cases. Even with the removal of data for 
$z/r_s<10$, the range of density ratios remains considerable and spans 
$0.5<\rho _m/\rho _0 < 4$.
When plotting 
$\gamma _m$ and 
$- \delta _m$ as a function of the density ratio 
$\rho _m / \rho _0$, all simulations draw clear one-to-one dependencies of the two coefficients on the density ratio (figure 7a,b), which are well fitted by a power law with exponents 
$-0.17$ and 
$-1.4$ for 
$\gamma _m$ and 
$\delta _m$, respectively. Also plotted in figure 7 are the experimental results by Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023) for a light jet (
$\rho _s/ \rho _0=0.4$) and an iso-density jet, and estimates inferred from the experiments by Charonko & Prestridge (Reference Charonko and Prestridge2017) for a heavy jet. Even though distributed over a narrower range of values of 
$\rho _m/ \rho _0$, the experimental results correspond closely to the numerical data, confirming that the descending trends of 
$\gamma _m$ and 
$- \delta _m$ with 
$\rho _m/ \rho _0$ observed here are a robust feature of these flows.
Log–log plot of the profile coefficients 
${\gamma }_m$ (a) and 
$-{\delta }_m$ (b) as a function of the density ratio. Symbols as in figure 3. Black pentagrams: Charonko & Prestridge (Reference Charonko and Prestridge2017). Black diamonds: Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023). Dashed lines represent power-law behaviours with slopes 
$-0.17$ for 
${\gamma }_m$ (a) and 
$-1.4$ for 
${\delta }_m$ (b).
Finally, we show in figure 8 the dependence of 
$\alpha$ on 
$\rho _m / \rho _0$ for all simulations (over the range 
$10 < z/r_s<30$) and including the experimental estimates by Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023) and Charonko & Prestridge (Reference Charonko and Prestridge2017). The black dashed line shows the prediction of 
$\alpha \approx \alpha _{prod}$ based on the power-law scaling observed in figure 7(a,b), indicating a decreasing trend with a much smaller exponent than those for 
$\gamma _m$ and 
$\delta _m$. Indeed, 
$\alpha _{prod} = (\rho _m / \rho _0) \delta _m / (2 \gamma _m)$, and 
$\delta _m / \gamma _m \sim (\rho _m/\rho _0)^{-1.23}$ (see figure 7), implying that 
$\alpha _{prod} \sim (\rho _m/\rho _0)^{-0.23}$, i.e. only very weakly dependent on 
$\rho _m / \rho _0$. This slight decreasing trend is lost when taking into account the contribution of 
$\alpha _{shape}$. Although the data range is too limited to make definite statements, no clear dependence of 
$\alpha$ on the density ratio is discernible.
Entrainment coefficient as a function of the density ratio. Symbols as in figures 3 and 7. Shown as black diamonds are the experimental data of Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023) and as black stars the experimental data of Charonko & Prestridge (Reference Charonko and Prestridge2017). The black dotted line shows 
$\alpha = 0.07$, the dashed black line represents the power-law behaviour of 
$\alpha _{prod}$, while the solid red line represents the R&S scaling for 
$\alpha$ (i.e. (1.2), with 
$\alpha _0=0.068$).
Figure 8 also shows (red solid line) the relation provided by the R&S scaling (1.2), which is evidently not in accordance with the numerical and experimental data. To get insight on the mismatch between the data and the R&S scaling, we can investigate the latter's validity by manipulating the entrainment relation (3.1). According to (3.1), a sufficient condition for the R&S result 
$\beta = 1/2$ to hold is for the TKE production to scale as 
$\alpha _{prod} \sim ({\rho _m} / {\rho _0})^{1/2}$ and the contribution of the shape term 
$\alpha _{shape}$ to be negligible.
Although the shape term 
$\alpha _{shape}$ vanishes for the light jets, this is not the case for the heavy jets, which are instead characterised by a persistent potential core, whose extent increases for increasing density ratios. In the case of heavy releases, it is therefore evident that density differences affect the flow in several distinct respects that are likely to be difficult/impossible to decouple in the near field.
If self-similarity had to be recovered, therefore, this would be only for low-density releases, given their tendency to squeeze the potential core in few source radii. Similarity solutions should, however, be consistent with the evolution of 
$\alpha _{prod}$, whose scaling, determined here based on empirical evidence, goes as 
$({\rho _m} / {\rho _0})^{-0.23}$, i.e. different from that initially suggested by R&S. Strictly speaking, this power-law behaviour of 
$\alpha _{prod}$ with density ratio may not be valid outside the range of density ratios studied.
6. Conclusions
We studied the effect of a variable density 
$\rho _m$ on the turbulent entrainment of ambient fluid, of density 
$\rho _0$, within axisymmetric jets, focusing on the evolution of the entrainment coefficient 
$\alpha$ and its dependence on the density ratio 
$\rho _m / \rho _0$.
Early laboratory experiments conducted by R&S quantified entrainment in the far field of variable-density jets and suggested a scaling described by 
$\alpha \propto (\rho _m / \rho _0)^{1/2}$. However, recent experimental and numerical analyses that focused on jets’ near field (Salizzoni et al. Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023) have challenged this scaling. In order to address this apparent contradiction, we revisited the original experimental data by R&S in § 2. We demonstrated that, since the R&S experiments focused on the far-field behaviour of the jet, their data do not provide sufficient evidence for a definitive conclusion about the dependence of 
$\alpha$ on 
$\rho _m / \rho _0$.
To investigate the scaling of 
$\alpha$ on 
$\rho _m / \rho _0$, we performed large-eddy simulations, verified against experimental data for two reference cases (Salizzoni et al. Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023), of jets having a density at the source varying over two orders of magnitude. Following Salizzoni et al. (Reference Salizzoni, Vaux, Creyssels, Amielh, Pietri and Anselmet2023), the interpretation of the results relies on the entrainment decomposition by van Reeuwijk & Craske (Reference van Reeuwijk and Craske2015), whose formulation was conveniently extended to deal with large-density releases. This framework allows the entrainment coefficient to be expressed as a function of two integral quantities, 
$\gamma _m$ and 
$\delta _m$ expressing the effect of a varying shape of the mean radial velocity profiles and that of the TKE production, respectively. Numerical results show that an increasing density ratio 
$(\rho _m / \rho _0 )$ has a major effect in reducing both 
$\gamma _m$ and 
$\delta _m$. This tendency is in very good agreement with that depicted by recent experimental data, obtained over a narrower rage of 
$\rho _m / \rho _0$. These modifications reflect the tendency of heavier releases to delay, compared with lighter ones, the TKE production and therefore the modification of the radial profiles of the mean velocity. Despite this, their combined effect in the overall entrainment process tends to cancel, so that the entrainment coefficient turns out to be relatively insensitive to variations of the density ratio 
$(\rho _m / \rho _0 )$.
Our observations are, therefore, not consistent with the scaling 
$\alpha \sim (\rho _m / \rho _0 )^{1/2}$ originally proposed by R&S. This has significant implications for the mathematical modelling of variable-density jets, since the R&S entrainment relation is indeed the only one consistent with full self-similarity of the flow (for similarity solutions, see e.g. van den Bremer & Hunt (Reference van den Bremer and Hunt2010)). Therefore, our findings suggest that variable-density jets are not expected to be self-similar in those (typically near-field) parts of the flow where the effects of large density differences are significant.
However, the simulations also clearly show that variable-density effects on entrainment are compounded by entrance effects, particularly the existence of a potential core. We showed that the length of the potential core has a strong dependence on the density ratio, which, in turn, has a direct effect on the profile coefficient 
$\gamma _m$ which affects 
$\alpha$. Entrance effects cannot be expected to have any universality since they are intimately related to the manner in which the fluid is injected into the domain and will therefore be configuration-specific. For this reason, it will be challenging to make further progress regarding the scaling of 
$\alpha$ without resorting to idealised simulations. One way to circumvent the potential-core issue is to perform variable-density simulations of a planar temporal jet (e.g. Cimarelli et al. Reference Cimarelli, Mollicone, van Reeuwijk and De Angelis2021). For this flow, the jet evolves in time rather than in space, and it will be possible to study the dependence of 
$\alpha$ on the density ratio once the initial transients have passed.
In a more general way, there is a clear need for further experimental and numerical work on variable-density releases covering a large range of source conditions. This is required to widen our knowledge of the dynamics of these flows, in the case of both momentum-dominated releases (variable-density jets) and buoyancy-dominated releases, usually referred to in the literature as non-Boussinesq plumes (Rooney & Linden Reference Rooney and Linden1996; Woods Reference Woods1997; Lanzini et al. Reference Lanzini, Marro, Creyssels and Salizzoni2024).
Funding
M.v.R. acknowledges financial support of the Engineering and Physical Sciences Research Council (EPSRC) through the project ‘Multi-scale dynamics at the turbulent/non-turbulent interface of jets and plumes’ – grant EP/R042640/1.
Declaration of interests
The authors report no conflict of interest.
