Self-similarity of fluid residence time statistics in a turbulent round jet

Abstract

Fluid residence time is a key concept in the understanding and design of chemically reacting flows. In order to investigate how turbulent mixing affects the residence time distribution within a flow, this study examines statistics of fluid residence time from a direct numerical simulation (DNS) of a statistically stationary turbulent round jet with a jet Reynolds number of 7290. The residence time distribution in the flow is characterised by solving transport equations for the residence time of the jet fluid and for the jet fluid mass fraction. The product of the jet fluid residence time and the jet fluid mass fraction, referred to as the mass-weighted stream age, gives a quantity that has stationary statistics in the turbulent jet. Based on the observation that the statistics of the mass fraction and velocity are self-similar downstream of an initial development region, the transport equation for the jet fluid residence time is used to derive a model describing a self-similar profile for the mean of the mass-weighted stream age. The self-similar profile predicted is dependent on, but different from, the self-similar profiles for the mass fraction and the axial velocity. The DNS data confirm that the first four moments and the shape of the one-point probability density function of mass-weighted stream age are indeed self-similar, and that the model derived for the mean mass-weighted stream-age profile provides a useful approximation. Using the self-similar form of the moments and probability density functions presented it is therefore possible to estimate the local residence time distribution in a wide range of practical situations in which fluid is introduced by a high-Reynolds-number jet of fluid.

1 Introduction

In many industrial, environmental and biological flows, the extent to which chemical, physical or biological processes progress within a given system depends upon the residence time of fluid within the system. The residence time distribution for a system provides information that may, in some situations, be related to the distribution of the progress of processes within the system (Langmuir 1908; Danckwerts 1953; Spalding 1958). The residence time distribution is used for modelling processes in the food industry (Fayolle, Belhamri & Flick 2013), mineral processing (Mulenga & Chimwani 2013), wastewater management (Moullec et al. 2008) and combustion (Costa, Silva & Azevedo 2003). Understanding of the relationship between the dynamics of a flow and the residence time distributions throughout a flow therefore might offer valuable insights into the effects of flow features on the progress of processes taking place within the fluid.

In special cases involving steady flow and first-order kinetic processes there is a unique relationship between the residence time distribution and the extent of progress of the kinetic processes (Langmuir 1908). In flows involving more complex kinetic processes, generally there is not a direct relationship between the residence time distribution and the extent of progress of kinetic processes, but the residence time distribution may still provide information regarding the likely extent of these processes. Combustion is an important example of a complex kinetic process in which the extent of the combustion reactions and of pollutant formation in some cases can be related to fluid residence time. As a first example, waste incinerators rely upon adequate residence time of combustion gases within a high-temperature region of the combustor in order to ensure burn-out of pollutants: in particular, a minimum one-second residence time above $800\,^{\circ }\text{C}$ for small-scale incinerators is specified in the best practice guide by the World Health Organization (Batterman 2004). As a second example, residence time is the principal indicator of where and when an ignitable mixture will autoignite, and Gomet, Robin & Mura (2012) used residence time statistics to predict the locus of autoignition in a supersonic non-premixed flow by comparing the residence time with the ignition delay time of the fuel–air mixture. A further example where residence time is used to characterise the extent of physical processes in combustion is given by Nambully et al. (2014), who used the mean residence time to explain the accumulation of preferential diffusion effects in a bluff-body stabilised turbulent flame. In contrast to the preceding examples that use residence time to characterise the extent of kinetic processes, recent works in the context of turbulent reacting flow modelling (Bilger, Kim & Martin 2004; Grout 2007; Enjalbert, Domingo & Vervisch 2012) have used information concerning the fluid residence time as a modelling input in order to improve modelling of the rate of progress of kinetic processes. These examples serve to emphasise that the fluid residence time distribution has a wide domain of utility, not restricted to simple situations involving steady-flow reactors and first-order kinetics.

Several definitions of fluid residence time appear in the literature and a clarification of the differences is provided here. Danckwerts’ usage of residence time (Danckwerts 1953) leads to its definition as the cumulative time that a conserved element of material has spent within a specified domain. A conserved element of material can be considered as a distinct physical object whose substance is not changed by the processes taking place in the flow, such as an atom in the absence of nuclear reaction, or a molecule that does not undergo chemical change. The instantaneous value of the residence time of an element of material is described in the following as the elementary age and denoted by $a^{\ast }(t)$ , with the corresponding mass denoted by $m^{\ast }$ . A number of experimental and numerical approaches have been developed previously in order to characterise and model residence time distributions, as described in the popular text book by Levenspiel (1999) and in a recent review by Nauman (2008), and these methods are applied routinely in, among others, the process industry. The probability density function of the elementary age at a point within a turbulent flow is not readily amenable to experimental measurement due to the difficulty of tracking the motion of individual elements of material such as molecules or atoms. In principle, direct numerical simulations describing the motion of individual elements of material, such as molecular dynamic simulations, provide a means to evaluate distributions of elementary age within a turbulent flow. However, molecular dynamic simulations of practical turbulent flow remain computationally infeasible. Instead, elementary-age distributions can be modelled by Lagrangian tracking of material elements (Mackley & Saraiva 1999), provided that the effects of Brownian and turbulent motions are accounted for accurately, at least in a statistical sense. Rather than simulating the molecular dynamics directly, random walk procedures (Mackley & Saraiva 1999; Nauman 2008) and Langevin equations (Langevin 1908; Haworth & Pope 1986) have been used to take account of molecular and turbulent contributions to the transport of the elementary age.

In contrast to the elementary age, the fluid age ( $a_{M}$ ) is a continuum property of a fluid that is defined by a transport equation in the form proposed by Sandberg (1981). Ghirelli & Leckner (2004) generalised Sandberg’s equation for variable property flow. For flow of a fluid consisting of one or more chemical species with equal diffusivities, the fluid-age transport equation is given by

(1.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}a_{M}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}a_{M}=\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}{\mathcal{D}}\unicode[STIX]{x1D735}a_{M})+1,\end{eqnarray}$$

where $\boldsymbol{u}$ , $\unicode[STIX]{x1D70C}$ and ${\mathcal{D}}$ are the velocity vector, the density and the molecular diffusivity of the fluid. The first term in (1.1) represents unsteadiness of the local fluid age, the second term represents advection by the convective velocity $\boldsymbol{u}$ , the third term represents the effects of diffusive flux of material between regions with different values of fluid age, and the final source term represents the process of age increasing due to the passage of time. Sandberg & Sjöberg (1983) and Balo & Cloirec (2000) used this equation in conjunction with turbulence closure models to simulate the residence time in turbulent flows. Note that the fluid age may be interpreted as an instantaneous mass-weighted average of the elementary age of atoms in an infinitesimal volume within the flow (Ghirelli & Leckner 2004).

In a turbulent flow, the one-point one-time probability density function of the elementary age has two contributions: the first is due to the distribution of elementary age present within one realisation of the turbulent flow; and the second is due to turbulent fluctuations in the flow. In contrast, the one-point one-time probability density function of fluid age shows variations only due to turbulent fluctuations in the flow. The fluid age is therefore a useful metric for isolating the effect of turbulent flow features on the residence time distribution. Accordingly, the concept of fluid age is employed for the purposes of modelling or analysis in several turbulent combustion studies (Bilger et al. 2004; Grout 2007; Enjalbert et al. 2012; Nambully et al. 2014).

In order to characterise the residence time of fluid from different sources or streams, Ghirelli & Leckner (2004) proceeded to derive a further transport equation that defines stream age $a_{i}$ , which can be interpreted as the mass-weighted elementary age of material originating from the $i$ th stream of fluid. The stream age is undefined when the mass fraction $\unicode[STIX]{x1D709}_{i}$ , i.e. the fraction of mass originating from stream $i$ , is equal to zero. Avoiding singularities where $\unicode[STIX]{x1D709}_{i}=0$ , Ghirelli & Leckner (2004) defined the stream age via a transport equation for the mass-weighted stream age $\unicode[STIX]{x1D6F7}_{i}=\unicode[STIX]{x1D709}_{i}a_{i}$ ,

(1.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F7}_{i}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D6F7}_{i})=\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}{\mathcal{D}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}_{i})+\unicode[STIX]{x1D709}_{i}.\end{eqnarray}$$

The mass-weighted stream age and fluid age are related by $\sum _{i=1}^{N}\unicode[STIX]{x1D6F7}_{i}=a_{M}$ where $N$ is the total number of streams considered. Note that (1.1) and (1.2) are valid given the assumption that all species have the same diffusivity ${\mathcal{D}}$ , although a more general set of equations without this restriction is presented by Ghirelli & Leckner (2004). Appendix A provides an alternative to the Ghirelli & Leckner (2004) derivation of (1.1) and (1.2) that makes explicit the link between the dispersion of elementary age by Brownian motion and the diffusion of the fluid age or stream age.

Provided that all of the conserved elements of material being tracked have the same mass, it follows that there is a linear relationship between the mass-weighted and mole-weighted stream age. In chemically reacting flows it is pertinent to track the residence time of either atoms or molecules. Since atoms and molecules have fixed masses related to their elemental composition, the mass-weighted and mole-weighted stream-age distributions convey equivalent information. A mass-based presentation is retained in the following for consistency with Ghirelli & Leckner (2004).

The mass fraction $\unicode[STIX]{x1D709}_{i}$ is governed by the following transport equation:

(1.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{i}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D709}_{i})=\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}{\mathcal{D}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}_{i}),\end{eqnarray}$$

where the sum of all the mass fractions is unity, i.e. $\sum _{i=1}^{N}\unicode[STIX]{x1D709}_{i}=1$ .

The stream-age equation is useful because a broad set of chemical processes require the presence of more than one reactant, so that reaction progress may be related to the time since individual streams were introduced, rather than the fluid age of the mixture (Mouangue et al. 2014). For example, Gomet et al. (2012) used an ensemble-averaged form of the stream-age equation in order to model the different residence times of fuel, air and fuel–air mixtures in an autoigniting flow, and subsequently obtained accurate predictions for the ignition locations.

Despite the long-standing use of residence times for analysis of systems in process engineering, there is a lack of information concerning the characteristics of residence time distributions at locations within turbulent flows. In order to provide fundamental information concerning the statistical distribution of fluid residence time within turbulent flows, this study analyses stream-age statistics obtained from new direct numerical simulation (DNS) data for a turbulent jet flow. Given the lack of previous studies of this nature, we investigate the turbulent round jet because it is a canonical example of a turbulent free shear flow and it exhibits features found in many flows in engineering and nature.

The statistically stationary turbulent round jet has been studied extensively. A key feature of such jets is that the velocity (Wygnanski & Fiedler 1969; Hussein, Capp & George 1994) and mass fraction (Mi, Nathan & Nobes 2001) statistics become self-similar downstream from an initial development region when normalised by centreline values. Various measurements confirm the self-similarity of the velocity statistics up to the fourth-order moments of the velocity vector (Panchapakesan & Lumley 1993). Furthermore, the approximate solution for this self-similarity can be derived theoretically for the mean velocity by using a mixing length hypothesis (Tollmien 1926) or an eddy viscosity model across the flow (Schlichting & Gersten 2000). As the jet develops further downstream, the momentum flux remains constant, but the mass flow rate increases. The increase of the mass flow rate implies that the jet entrains surrounding fluid. In order to measure the mass entrainment rate, Ricou & Spalding (1961) devised a porous-walled cylindrical chamber that controlled the radial inflow of fluid. They determined the mass entrainment rate from the flow rate through the porous wall that was required in order to maintain a uniform pressure along the axis of the jet. Last, recent studies on the statistically stationary turbulent jet focused on the effect of the jet boundary condition. Xu & Antonia (2002) tested various nozzle configurations and obtained different coefficients for half-radius profile and entrainment rates, confirming the dependence of the self-similar jet profile on the inlet boundary conditions. Although the half-radius profiles and the entrainment rates depend on the inflow conditions, the self-similar profiles of velocities are found to be independent of the inlet conditions. Measurements with pipe inlets, converging nozzles and fractal nozzles all show comparable self-similar profiles after the initial development region.

The primary objective of this study is to investigate the implications of the self-similarities of velocity and mass fraction statistics for the fluid-age distribution in turbulent round jets. To this end, the following sections introduce the numerical simulations and analysis used in order to obtain residence time statistics for a turbulent round jet; demonstrate self-similarity of the mass-weighted stream-age statistics; and derive a model that provides an estimate for the self-similar profiles of mean mass-weighted stream age as a function of the profiles of the mean velocities and mass fraction.

2 Numerical simulation

The simulation configuration consists of a round jet of turbulent fluid issuing from a flat plate into a quiescent environment. The injected fluid is a perfect gas with the same temperature and density as the ambient fluid. The jet Reynolds number $Re=(U_{0}D)/\unicode[STIX]{x1D708}$ is 7290 and the Mach number is 0.304, based on the bulk velocity $U_{0}$ , jet diameter $D$ and kinematic viscosity $\unicode[STIX]{x1D708}$ . The Taylor Reynolds number on the centreline at $x/D=15$ is calculated to be $Re_{\unicode[STIX]{x1D706}}=35$ . Two streams are considered in the simulation: the jet fluid with a subscript $j$ and the ambient fluid. Transport equations for the mass fraction of fluid originating from the jet ( $\unicode[STIX]{x1D709}_{j}$ ) and for the fluid age ( $a_{M}$ ) are solved with the jet fluid assigned zero age as it exits from the injector.

The mass fraction of jet fluid in the ambient fluid is initialised to zero and arbitrarily we initialise the age of the ambient fluid to zero at the start of the simulation. This implies that the stream age for the ambient fluid is uniform and equal to the time $t$ , and that the mass-weighted stream age of the jet fluid ( $\unicode[STIX]{x1D6F7}_{j}$ ) is given by

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{j}=a_{M}-(1-\unicode[STIX]{x1D709}_{j})t.\end{eqnarray}$$

The mass fraction ( $\unicode[STIX]{x1D709}_{j}$ ) and the stream age ( $a_{j}$ ) of jet fluid at the inlet are set to unity and zero, respectively, across the inlet. The domain in which the source term in the age equations is applied is the entire simulation domain downstream of the inlet. The jet inlet mean velocity is prescribed as a smoothed top-hat profile. The mean value of axial velocity at the inlet is uniform out to a radius of $r=0.475D$ and then smoothly drops to zero at $r=0.5D$ following a half-cosine function. Away from the jet inlet ( $r\geqslant 0.5D$ ), a no-slip wall boundary condition is imposed at $x=0$ , where $x$ is the streamwise coordinate. Pseudo-turbulent homogeneous isotropic velocity fluctuations are superimposed at the inlet using the digital filter method (Touber & Sandham 2009) with a low turbulence intensity of 1.7 %. The low turbulence intensity is chosen to be similar to the nozzle inlet measurements of Mi et al. (2001), where the measured inlet turbulent intensity was  ${\sim}1\,\%$ .

All the other boundaries are non-reflecting outlets (Poinsot & Lele 1992) with a small buffer region (Sandberg & Sandham 2006) at the downstream outlet boundary where artificial viscosity is used to damp the characteristic waves in order to minimise spurious acoustic reflection from the boundaries. All scalar diffusivities ( ${\mathcal{D}}$ ) are set by assuming Lewis numbers equal to unity and a Prandtl number equal to 0.72. Viscosity changes with temperature following Sutherland’s law (Sutherland 1893).

The flow is simulated with the compressible DNS code HiPSTAR, developed at the University of Southampton and recently used for a jet noise study (Sandberg & Tester 2016). A fourth-order central finite difference scheme with energy-conserving boundary schemes of the same order (Nordström & Carpenter 1998) is used in the longitudinal and the radial directions, while a pseudo-spectral method is used in the circumferential direction. The singularity of the axis resulting from the cylindrical coordinate system is treated using parity conditions (Sandberg 2011). A fourth-order low-memory Runge–Kutta scheme (Kennedy, Carpenter & Lewis 2000) is used for time advancement. In addition, skew-symmetric splitting of the nonlinear terms is used to enhance the stability (Kennedy & Gruber 2008).

For the computational mesh, a structured and stretched grid is used, adapted from a previous round jet study (Sandberg, Sandham & Suponitsky 2012) where $Re=6300$ . The grid spacing, $\unicode[STIX]{x1D6E5}$ , from the previous study was refined using Reynolds number scaling as $\unicode[STIX]{x1D6E5}/D\propto Re^{-3/4}$ (Antonia, Satyaprakash & Hussain 1980), since the turbulent Reynolds numbers in the jet scale with the jet Reynolds number $Re$ to a good approximation. In the radial direction, the grid is most refined near the edge of the jet inlet ( $r=D/2$ ) where the velocity and scalar gradients are the greatest, and 145 points are assigned radially within the jet diameter. In the axial direction, the grid is most refined near the inlet and gradually stretched moving downstream. In the circumferential direction, 64 Fourier modes are used, corresponding to 130 physical collocation points. The grid consists of $3020\times 834\times 130$ structured nodes, spanning axially from $x=0{-}60D$ and radially from $r=0{-}30D$ .

Figure 1. The ratios between the grid spacing and the calculated Kolmogorov length scale along (a) the axial direction and (b) the radial direction.

In order to resolve all the turbulent flow structures, the grid spacing should have the same order of magnitude as the Kolmogorov length scale ( $L_{K}$ ). In particular, Yeung & Pope (1989) recommended use of $\unicode[STIX]{x1D6E5}/L_{K}\approx 2.1$ . Figure 1(a,b) shows the ratios between the grid spacing and the calculated Kolmogorov length scale. $\unicode[STIX]{x0394}x$ and $\unicode[STIX]{x0394}r$ vary as the grid is stretched, and the Kolmogorov length scale also changes spatially as noted by Antonia et al. (1980). The values of $\unicode[STIX]{x1D6E5}/L_{K}$ along the axial direction at $r/D=0$ and $r/D=0.5$ and across the radial direction at $x/D=0.5$ and $x/D=15$ shown in figure 1(a,b) remain below 2.5, confirming that the grid resolution is adequate.

In order to accelerate the development of the statistically stationary jet flow field, the flow is simulated for 540 jet times ( $\unicode[STIX]{x1D70F}=D/U_{0}$ ) using a computational mesh with half the resolution of the final grid. By 540 jet times the first- and second-order statistics in the first 30 diameters of the domain show that the simulation has reached a statistically stationary state. Then, the coarse-grid solution is interpolated onto the final mesh, and the simulation continued over an additional $80\unicode[STIX]{x1D70F}$ , confirming that a statistically stationarity state is re-established. Subsequently, statistics are sampled over 380 jet times. The converged turbulent jet simulation also displays self-similarity downstream of 40 jet diameters as discussed below.

Figure 2. Radial profiles of (a) the mean axial velocity, (b) the axial velocity fluctuation, (c) the axial–radial Reynolds stress and (d) the energy spectrum of axial velocity.

Figure 2(a–c) shows the first and second moments of the axial and radial velocity components from the simulation versus the scaled radius ( $\unicode[STIX]{x1D702}=r/(x-x_{0})$ ), where $x_{0}$ is a virtual origin. The experiment of Panchapakesan & Lumley (1993) was chosen for the comparison because they used a similar boundary condition to that in the present DNS, involving a near-laminar jet ( $u_{rms}\approx 0.5\,\%$ ) with a top-hat velocity profile from a convergent nozzle issuing through a flat plate. All the mean quantities of the current simulation show adequate agreement with the laboratory measurements. Panchapakesan & Lumley (1993) reported that the higher moments reach a self-similar state for $x/D\geqslant 70$ . We noted that the statistics achieved self-similarity already after $x/D\geqslant 15$ for the first moment and $x/D\geqslant 25$ for the second moment, and note that this might be attributable to the lower Reynolds number of the DNS, since the jet Reynolds number is 11 000 in the measurements reported by Panchapakesan & Lumley (1993).

Figure 2(d) shows the energy spectrum of the axial velocity component, $E_{11}$ . The energy spectrum is obtained from a single-point measurement and the temporal variation of velocity is converted to a spatial dependency employing Taylor’s hypothesis, i.e. the turbulent fluctuations are assumed to move downstream with the averaged velocity. As $E_{11}$ is the energy spectrum of a single component, the peak value occurs near $k=0$ (Pope 2000). Over the wavenumber range $kD=2{-}20$ , $k^{-5/3}$ scaling is obtained, indicating that the flow is locally isotropic. The range of the spectrum also confirms that the Kolmogorov length scales are resolved in the simulation.

Figure 3. The radial variation of (a) the mean mass fraction and (b) the root mean square fluctuation of mass fraction.

Figure 3 shows the mean and the root mean square fluctuations of the jet fluid mass fraction. The contracting nozzle experiments of Mi et al. (2001) are shown for comparison as their inlet velocity conditions match approximately with those in the current DNS. The mean of the current DNS shows good agreement with the measurements, while the level of fluctuations predicted is slightly lower than the level in the laboratory observations. The difference in the level of fluctuations may be attributable to the differences in inlet conditions, including a 50 K temperature difference between the jet fluid and 288 K ambient fluid in the experiment by Mi et al. (2001).

The statistical distributions of the fluid properties can be further characterised by their higher moments. The third-order (skewness) and the fourth-order (kurtosis) moments of jet fluid mass fraction are defined as follows:

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \mathit{skewness}=\frac{\overline{(\unicode[STIX]{x1D709}_{j}-\overline{\unicode[STIX]{x1D709}_{j}})^{3}}}{(\unicode[STIX]{x1D709}_{j}^{rms})^{3}}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \mathit{kurtosis}=\frac{\overline{(\unicode[STIX]{x1D709}_{j}-\overline{\unicode[STIX]{x1D709}_{j}})^{4}}}{(\unicode[STIX]{x1D709}_{j}^{rms})^{4}}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D709}_{j}^{rms}$ is the root mean square (r.m.s.) value. The skewness is a measure of symmetry for a given distribution – zero for a symmetric distribution, negative indicating a longer tail on the left of the distribution and positive indicating a longer tail on the right. The kurtosis is a measure of tailedness or peakedness – large if a distribution is spread out (or more influenced by tails) and small if a distribution has a more pronounced peak (or it is less influenced by tails). For reference, the skewness and the kurtosis of a Gaussian distribution are zero and three, respectively. Figure 4 shows the skewness and the kurtosis of the jet fluid mass fraction over the scaled radius at different axial locations. In addition, the measurements from Mi et al. (2001) have been added for comparison, showing close agreement between the DNS and the measurements. Note that measurements from Mi et al. (2001) are extracted between $x/D=25{-}65$ , while the DNS data are extracted between $x/D=10{-}30$ , and both the experimental and the numerical data lack statistical convergence in the outer portion of the jet ( $\unicode[STIX]{x1D702}\geqslant 0.2$ ).

Figure 4. The skewness (a) and the kurtosis (b) of the mass fraction over the scaled radius at different axial locations.

The normalised profiles of statistical quantities evaluated from the DNS are expected to be self-similar downstream of an initial development region. A degree of variation in the normalised profiles is evident between different axial locations, and this may be attributed partly to incomplete statistical convergence and partly to ongoing development of the flow in the region shown, $x/D=15{-}40$ . The statistics reported are evaluated over 380 jet times; however, the integral time scale of the turbulent flow increases approximately quadratically in the downstream direction, and higher-order statistical moments are strongly influenced by relatively infrequent intermittency events in the radially outer portion of the jet. Therefore the statistical convergence of the results varies through the flow, reducing with axial distance and for higher moments. Statistical moments of second order and higher are reported downstream of either $x/D=20$ or 25, which is similar to the upstream location $x/D=25$ at which Mi et al. (2001) report self-similarity of the first four moments of the mass fraction. Figure 4 also indicates that the statistics obtained from the present DNS exhibit a similar level of scatter compared to the previous experimental study by Mi et al. (2001), and this is considered adequate for the purposes of characterising statistics of the residence time distribution.

3 Results

In this section we first present an a priori analysis concerning the self-similarity of the mass-weighted stream-age profiles. Second, we demonstrate that the mass-weighted stream-age distribution is statistically stationary and, downstream of an initial development region, self-similar with respect to the scaled radius of the jet. Last, we will use the self-similar properties of the velocity and mass fraction fields to develop a model for the mean mass-weighted stream age, which we test using the DNS data.

3.1 A priori analysis of the mass-weighted stream-age profiles

The mass-weighted stream age ( $\unicode[STIX]{x1D6F7}_{j}$ ) is a passive scalar that depends on flow variables $\boldsymbol{u}$ and $\unicode[STIX]{x1D709}_{j}$ whose moments are known to be self-similar with scaled radius $\unicode[STIX]{x1D702}$ . However, it does not immediately follow that the mean mass-weighted stream age $\overline{\unicode[STIX]{x1D6F7}_{j}}$ is self-similar with $\unicode[STIX]{x1D702}$ for two reasons. First, we need to confirm that the statistics of $\unicode[STIX]{x1D6F7}_{j}$ are statistically stationary. It is shown in § 3.2 that the statistics of $\unicode[STIX]{x1D6F7}_{j}$ are indeed stationary, whereas those of the fluid age are not. Second, there are unclosed turbulent flux terms in the transport equation for $\overline{\unicode[STIX]{x1D6F7}_{j}}$ and more needs to be known about these before presuming that profiles of $\unicode[STIX]{x1D6F7}_{j}$ will be self-similar. The second point is illustrated below by demonstrating that certain forms for the turbulent flux terms would give profiles of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ that are not a function of the scaled radius, and therefore not self-similar.

The self-similar flow variables of $\bar{u}$ , $\bar{v}$ and $\overline{\unicode[STIX]{x1D709}_{j}}$ in turbulent round jets can be expressed as follows:

(3.1a-c ) $$\begin{eqnarray}\bar{u}=\frac{A}{x}f(\unicode[STIX]{x1D702}),\quad \bar{v}=\frac{A}{x}g(\unicode[STIX]{x1D702})\quad \text{and}\quad \overline{\unicode[STIX]{x1D709}_{j}}=\frac{B}{x}h(\unicode[STIX]{x1D702}),\end{eqnarray}$$

where $A$ and $B$ are constants known as decay ratios and $\unicode[STIX]{x1D702}=r/x$ . Note that the virtual origin, $x_{0}$ , is omitted in this section for brevity.

Next, assuming that $\unicode[STIX]{x1D6F7}_{j}$ is statistically stationary, the averaged equation for mass-weighted stream age can be written from (1.2). Substituting (3.1) into the average of (1.2) leads to

(3.2) $$\begin{eqnarray}f(\unicode[STIX]{x1D702})\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D6F7}_{j}}}{\unicode[STIX]{x2202}x}+g(\unicode[STIX]{x1D702})\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D6F7}_{j}}}{\unicode[STIX]{x2202}r}=-\frac{x}{A}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}(\overline{u^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }})-\frac{x}{Ar}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }})+\frac{B}{A}h(\unicode[STIX]{x1D702}),\end{eqnarray}$$

where the molecular diffusion terms are neglected for a simpler presentation. Equation (3.2) contains two unclosed turbulent flux terms, so the generic form for $\overline{\unicode[STIX]{x1D6F7}_{j}}$ cannot be determined. In order to illustrate how these unclosed terms influence the profile of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ , we construct hypothetical profiles for $\overline{u^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$ and $\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$ as follows:

(3.3a,b ) $$\begin{eqnarray}\overline{u^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}=F_{1}+F_{2}+F_{3}\quad \text{and}\quad \overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}=G_{1},\end{eqnarray}$$

where $F_{1}$ , $F_{2}$ , $F_{3}$ and $G_{1}$ are defined as the solutions of

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle F_{1}(x,r)=\int _{0}^{x/r}\frac{B}{\unicode[STIX]{x1D702}^{\prime }}h\left(\frac{1}{\unicode[STIX]{x1D702}^{\prime }}\right)\text{d}\unicode[STIX]{x1D702}^{\prime }, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}F_{2}(x,r)}{\unicode[STIX]{x2202}x}=A\unicode[STIX]{x1D6FC}x^{\unicode[STIX]{x1D6FC}-2}f\left(\frac{r}{x}\right)\unicode[STIX]{x1D711}\left(\frac{r}{x^{\unicode[STIX]{x1D6FD}}}\right), & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}F_{3}(x,r)}{\unicode[STIX]{x2202}x}=-A\unicode[STIX]{x1D6FD}x^{\unicode[STIX]{x1D6FC}-2}f\left(\frac{r}{x}\right)\frac{r}{x^{\unicode[STIX]{x1D6FD}}}\unicode[STIX]{x1D711}^{\prime }\left(\frac{r}{x^{\unicode[STIX]{x1D6FD}}}\right), & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{r}\frac{\unicode[STIX]{x2202}(rG_{1}(x,r))}{\unicode[STIX]{x2202}r}=Ax^{\unicode[STIX]{x1D6FC}-2}\frac{x}{r}g\left(\frac{r}{x}\right)\frac{r}{x^{\unicode[STIX]{x1D6FD}}}\unicode[STIX]{x1D711}^{\prime }\left(\frac{r}{x^{\unicode[STIX]{x1D6FD}}}\right), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D702}^{\prime }$ is a dummy variable, $\unicode[STIX]{x1D711}^{\prime }$ represents the derivative of $\unicode[STIX]{x1D711}$ , $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ are constants and $\unicode[STIX]{x1D711}$ is a function. The functions $F_{1}$ , $F_{2}$ , $F_{3}$ and $G_{1}$ are constructed in such a way that the right-hand side of (3.2) has a general form $\sum x^{(\unicode[STIX]{x1D6FC}-1)}C_{1}(r/x)C_{2}(r/x^{\unicode[STIX]{x1D6FD}})$ , where $C_{X}$ denotes a generic function. Hence, given the assumed turbulent closure terms, $\overline{\unicode[STIX]{x1D6F7}_{j}}$ becomes

(3.8) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D6F7}_{j}}=x^{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D711}\left(\frac{r}{x^{\unicode[STIX]{x1D6FD}}}\right).\end{eqnarray}$$

As $\unicode[STIX]{x1D6FD}$ is an arbitrary constant, $\overline{\unicode[STIX]{x1D6F7}_{j}}$ can have a dependence other than on $\unicode[STIX]{x1D702}$ . This example illustrates that further information about the form of the turbulent flux terms in (3.2) is required before presuming whether or not $\overline{\unicode[STIX]{x1D6F7}_{j}}$ is self-similar with  $\unicode[STIX]{x1D702}$ .

3.2 Evolution of the mass-weighted stream age

Fluid age increases with residence time, causing the stream age of the ambient fluid around the jet to increase continuously. However, as the jet fluid issues from the inlet with age equal to zero and as the turbulent jet reaches a statistically stationary state, the influx of new jet fluid balances the ageing of the jet fluid inside the turbulent jet. Figure 5(a) shows the mass-weighted stream age of the jet fluid $\unicode[STIX]{x1D6F7}_{j}$ at two positions in the flow over a period of $350\unicode[STIX]{x1D70F}$ , confirming that the evolution of $\unicode[STIX]{x1D6F7}_{j}$ achieves a statistically stationary state. By comparison, figure 5(b) shows that the fluid age at the same locations increases over time due to the contribution of the uniformly increasing stream age of the ambient fluid.

Figure 5. (a) Temporal profile of mass-weighted stream age at selected locations marked in figure 6(a,b); (b) temporal profile of fluid age at the same locations.

Having established that the mass-weighted stream age of the jet fluid is statistically stationary, samples of $\unicode[STIX]{x1D6F7}_{j}$ from around the statistically homogeneous circumferential direction are accumulated over time in order to obtain ensemble statistics as a function of the axial and radial positions. Figure 6 shows an instantaneous and an ensemble-averaged $\unicode[STIX]{x1D6F7}_{j}$ field on a cross-section through the jet centreline, as well as the radial profiles of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ at selected axial locations. Comparison of the instantaneous and ensemble-averaged fields of $\unicode[STIX]{x1D6F7}_{j}$ indicates that, similar to the velocity and mass fraction fields, the width of the $\unicode[STIX]{x1D6F7}_{j}$ profile appears to increase approximately linearly in the downstream direction. The ensemble average $\overline{\unicode[STIX]{x1D6F7}_{j}}$ increases monotonically in the downstream direction along the centreline. Since the mean mass fraction $\overline{\unicode[STIX]{x1D709}_{j}}$ decreases downstream, the increase of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ indicates that the fractional increase of stream age $\overline{a_{j}}$ exceeds the fractional decrease of mass fraction along the centreline. Conversely, the ensemble average $\overline{\unicode[STIX]{x1D709}_{j}}$ decreases monotonically in the radial direction, revealing that the radial decrease of mean mass fraction outweighs the radial increase in stream age.

Figure 6. (a) Instantaneous mass-weighted stream age, $\unicode[STIX]{x1D6F7}_{j}$ , (b) ensemble-averaged $\overline{\unicode[STIX]{x1D6F7}_{j}}$ fields on a section through the jet centreline and (c) the radial variation of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ at different axial locations.

Downstream of an initial development region ( $x/D\geqslant 7$ ), $\overline{\unicode[STIX]{x1D6F7}_{j}}$ increases approximately linearly along the centreline of the jet, as shown in figure 7(a). The linear increase of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ along the centreline can be explained by considering a Reynolds-averaged form of (1.2) for statistically stationary flow:

(3.9) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D6F7}_{j}}}{\unicode[STIX]{x2202}x}=\frac{\overline{\unicode[STIX]{x1D709}_{j}}}{\bar{u}}+R,\end{eqnarray}$$

where the $R$ term includes all radial transport terms and molecular and turbulent transport in the axial direction. The detailed contributions of budget terms for $\overline{\unicode[STIX]{x1D6F7}_{j}}$ are discussed in § 3.3. However, it is instructive to consider that since $\overline{\unicode[STIX]{x1D709}_{j,c}}$ and $\overline{u_{c}}$ both decrease linearly along the centreline of the jet, as shown in figure 7(a), the first term on the right-hand side of (3.9) is approximately constant in the developed region of the jet. Figure 7(b) shows that the first term on the right-hand side of (3.9) provides the dominant contribution to the centreline gradient of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ , and therefore explains why $\overline{\unicode[STIX]{x1D6F7}_{j}}$ increases approximately linearly along the jet centreline.

Figure 7. Axial variation of: (a) the centreline average of mass-weighted stream age $\overline{\unicode[STIX]{x1D6F7}_{j,c}}$ (right axis) and the reciprocals of the centreline average mixture fraction $1/\overline{\unicode[STIX]{x1D709}_{j,c}}$ and axial velocity $1/\overline{u_{c}}$ (left axis); (b) the axial gradient of $\overline{\unicode[STIX]{x1D6F7}_{j,c}}$ and the ratio of $\overline{\unicode[STIX]{x1D709}_{j,c}}$ and $\overline{u_{c}}$ along the centreline.

3.3 Self-similarity of mass-weighted stream age

In this section, the radial profiles of mass-weighted stream-age statistics are presented. All values are normalised by relevant centreline mean values, $\overline{\unicode[STIX]{x1D709}_{j,c}}$ , $\overline{u_{c}}$ or $\overline{\unicode[STIX]{x1D6F7}_{j,c}}$ , and the radial position is scaled by the axial distance downstream from the virtual origin, $\unicode[STIX]{x1D702}=r/(x-x_{0})$ . The virtual origin, $x_{0}$ , was obtained from figure 7(a) by extrapolating $\overline{u_{c}}$ upstream from the linear region. All data presented in this section are sampled in the range $15\leqslant x/D\leqslant 40$ and $15\leqslant x/D\leqslant 30$ for the first and second moments, respectively.

Figure 8 shows the radial profiles of the mean and r.m.s. mass-weighted stream age, along with the turbulent fluxes $\overline{u^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$ and $\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$ that appear in the axisymmetric, steady-state, constant-density and Reynolds-averaged transport equation for $\overline{\unicode[STIX]{x1D6F7}_{j}}$ :

(3.10) $$\begin{eqnarray}\overline{u}\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D6F7}_{j}}}{\unicode[STIX]{x2202}x}+\overline{v}\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D6F7}_{j}}}{\unicode[STIX]{x2202}r}={\mathcal{D}}\left[\frac{\unicode[STIX]{x2202}^{2}\overline{\unicode[STIX]{x1D6F7}_{j}}}{\unicode[STIX]{x2202}x^{2}}+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D6F7}_{j}}}{\unicode[STIX]{x2202}r}\right)\right]-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\overline{u^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}\right)-\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}\right)+\overline{\unicode[STIX]{x1D709}_{j}}.\end{eqnarray}$$

All of these radial profiles appear self-similar for the range of downstream positions shown. The radial profile of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ is compared with the radial profiles of the mean axial velocity and mean jet mass fraction in figure 8(a). The mean axial velocity and mean mass fraction display approximately Gaussian profiles, with the jet fluid mass fraction profile slightly wider than the axial velocity profile, consistent with previous observations for passive scalars with order-unity Schmidt numbers (e.g. Dowling & Dimotakis 1990). In contrast, the radial profile of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ has a flatter peak than the Gaussian profile, but a similar overall width to the mass fraction profile. The shape of the profile of the r.m.s. fluctuations of $\unicode[STIX]{x1D6F7}_{j}$ in figure 8(b) also differs from the shape of the $u^{rms}$ and $\unicode[STIX]{x1D709}_{j}^{rms}$ profiles shown in figures 2(b) and 3(b), respectively. In each case the location of the peak fluctuation level is close to the location of the maximum radial gradient of the corresponding mean quantity, which is close to the point where the variance production term due to the mean radial gradient is greatest. Due to the broad-peaked radial profile of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ , the location of the peak r.m.s. value $\unicode[STIX]{x1D6F7}_{j}^{rms}$ is radially outward from the location of the peak values of $u^{rms}$ and $\unicode[STIX]{x1D709}_{j}^{rms}$ . Panels 8(c) and 8(d) show the axial and radial turbulent fluxes $\overline{u^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$ and $\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$ versus the scaled radius. These profiles are qualitatively similar to the $\overline{u^{\prime }v^{\prime }}$ profile shown in figure 2(c), but the radial location of the peak turbulent fluxes lies at approximately $\unicode[STIX]{x1D702}=0.12$ , roughly mid-way between the location of the peak r.m.s. value of the velocity fluctuations at $\unicode[STIX]{x1D702}=0.06$ and the peak r.m.s. value of the $\unicode[STIX]{x1D6F7}_{j}$ fluctuations at $\unicode[STIX]{x1D702}=0.18$ .

Figure 8. Radial variation of (a) the normalised mean mass-weighted stream age $\overline{\unicode[STIX]{x1D6F7}_{j}}$ over the range $x/D=15{-}40$ ; (b) the normalised $\unicode[STIX]{x1D6F7}_{j}$ r.m.s. over the range $x/D=15{-}30$ ; (c) the normalised axial turbulent flux $\overline{u^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$ over the range $x/D=15{-}30$ ; and (d) the normalised radial turbulent flux $\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$ over the range $x/D=15{-}30$ .

The budget of (3.10) is shown in figure 9 versus the scaled radius. Since the moments appearing in (3.10) are self-similar with respect to the scaled radius, the data for each value of the scaled radius are normalised by the centreline value $\overline{u_{c}}\overline{\unicode[STIX]{x1D6F7}_{j,c}}/(x-x_{0})$ and averaged over axial positions in the range $x/D=15{-}40$ . The relative magnitudes of the terms in (3.10) show that axial convection and the source term are the dominant contributions, consistent with the centreline analysis shown in figure 7(b). The next most significant terms are the radial transport terms due to the mean convection $\overline{v}\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D6F7}_{j}}/\unicode[STIX]{x2202}r$ and turbulent flux $-(1/r)\unicode[STIX]{x2202}(r\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }})/\unicode[STIX]{x2202}r$ . Diffusive transport of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ and axial transport by turbulent velocity fluctuations play a relatively minor role. The relative importance of the terms in the $\overline{\unicode[STIX]{x1D6F7}_{j}}$ budget is similar to the relative magnitudes of the corresponding terms in the mean axial velocity or mean mass fraction equations. The mean convection terms and the source term $\overline{\unicode[STIX]{x1D709}_{j}}$ are closed at first order, and the relative magnitude of diffusive transport of mean properties is expected to decrease in high-Reynolds-number turbulent flow and may be neglected; however, additional modelling is required in order to close the turbulent flux terms. Noting that the axial turbulent flux term is small compared to the radial turbulent flux, a model for the radial turbulent flux is developed in § 3.4.

Figure 9. The budget of (3.10) versus the scaled radius. Budget terms are normalised by $\overline{u_{c}}\overline{\unicode[STIX]{x1D6F7}_{j,c}}/(x-x_{0})$ and averaged over normalised axial distances $x/D=15{-}40$ .

3.4 Closure model for $\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$

In this section, we propose closure models for the radial turbulent flux of the mass-weighted stream age, $\overline{v^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }}$ . The axial and radial turbulent flux terms in (3.10) for the Reynolds average of $\unicode[STIX]{x1D6F7}_{j}$ are unclosed and both require modelling. However, figure 9 illustrates that the axial turbulent flux term, $\overline{u^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }}$ , has a minor influence on the overall transport budget, suggesting that it may be neglected, and here we focus only on developing a closure for the radial turbulent flux term $\overline{v^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }}$ in terms of first-order moments. Since the velocity, mass fraction and mass-weighted stream-age statistics are all self-similar, it is convenient to develop the closure model in terms of self-similar variables that are normalised by centreline values and that depend only on the scaled radius  $\unicode[STIX]{x1D702}$ . Due to the very similar shapes of the normalised mean mass fraction and normalised mean axial velocity profiles in the turbulent jet shown in figure 8, a number of different model forms can be used to approximate $\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}$ . Initially two model forms are considered:

(3.11) $$\begin{eqnarray}\text{Model 1:}~\frac{\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}}{\overline{u_{c}}\overline{\unicode[STIX]{x1D6F7}_{j,c}}}=-\unicode[STIX]{x1D6FD}_{1}\frac{\unicode[STIX]{x2202}(\overline{\unicode[STIX]{x1D6F7}_{j}}/\overline{\unicode[STIX]{x1D6F7}_{j,c}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\end{eqnarray}$$

and

(3.12) $$\begin{eqnarray}\text{Model 2:}~\frac{\overline{v^{\prime }\unicode[STIX]{x1D6F7}_{j}^{\prime }}}{\overline{u_{c}}\overline{\unicode[STIX]{x1D6F7}_{j,c}}}=-\unicode[STIX]{x1D6FD}_{2}\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x2202}(\bar{u}/\overline{u_{c}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{i}$ are modelling constants. Model 1 is proposed based on the eddy viscosity model for $\overline{u^{\prime }v^{\prime }}$ (Boussinesq 1877; Prandtl 1925). Model 2 was obtained by trying various combinations of flow variables. However, Model 2 has the undesirable property, from a modelling perspective, that it contains neither $v$ nor $\unicode[STIX]{x1D6F7}_{j}$ . Hence a third model was developed from Model 2 by using the continuity equation to obtain an expression in terms of $\bar{v}$ and $\overline{\unicode[STIX]{x1D6F7}_{j}}$ .

The Reynolds-averaged axisymmetric steady-state constant-density continuity equation can be rewritten in terms of the self-similarity variable as

(3.13) $$\begin{eqnarray}\frac{\bar{u}}{\overline{u_{c}}}+\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x2202}(\bar{u}/\overline{u_{c}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}-\frac{1}{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\left\{\unicode[STIX]{x1D702}\frac{\bar{v}}{\overline{u_{c}}}\right\}=0.\end{eqnarray}$$

Due to the similarity between the axial velocity and the $\overline{\unicode[STIX]{x1D6F7}_{j}}$ profiles, approximating $\bar{u}/\overline{u_{c}}$ with $\overline{\unicode[STIX]{x1D6F7}_{j}}/\overline{\unicode[STIX]{x1D6F7}_{j,c}}$ leads to

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x2202}(\bar{u}/\overline{u_{c}})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\approx \frac{1}{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\left\{\unicode[STIX]{x1D702}\frac{\bar{v}}{\overline{u_{c}}}\right\}-\frac{\overline{\unicode[STIX]{x1D6F7}_{j}}}{\overline{\unicode[STIX]{x1D6F7}_{j,c}}}.\end{eqnarray}$$

This approximation constitutes Model 3 as follows:

(3.15) $$\begin{eqnarray}\text{Model 3:}~\frac{\overline{v^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }}}{\overline{u_{c}}\overline{\unicode[STIX]{x1D6F7}_{j,c}}}=-\unicode[STIX]{x1D6FD}_{3}\left(\frac{1}{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\left(\unicode[STIX]{x1D702}\frac{\bar{v}}{\overline{u_{c}}}\right)-\frac{\overline{\unicode[STIX]{x1D6F7}_{j}}}{\overline{\unicode[STIX]{x1D6F7}_{j,c}}}\right).\end{eqnarray}$$

Figure 10 shows the comparison between $\overline{v^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }}$ and the three models proposed in (3.11), (3.12) and (3.15). The present data indicate which forms for the turbulent flux model are plausible; however, data for a wider range of flow conditions should be used to establish the general validity of the modelling. Model 1 differs substantially from the observed profile of $\overline{v^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }}$ , but Models 2 and 3 give close agreement with the actual profile. Models 2 and 3 are very close to one another, supporting the validity of the approximation made in (3.14). From this comparison, it is clear that either Model 2 or Model 3 provides a good approximation. As Model 2 contains neither $v$ nor $\unicode[STIX]{x1D6F7}_{j}$ , Model 3 would be preferred.

Figure 10. Comparison of the models proposed in (3.11), (3.12) and (3.15) with the DNS data of normalised radial flux of mass-weighted stream age, $\overline{v^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }}/\overline{u_{c}}\overline{\unicode[STIX]{x1D6F7}_{j,c}}$ (shaded area), plotted against the scaled radius,  $\unicode[STIX]{x1D702}$ . Fitting parameters are $\unicode[STIX]{x1D6FD}_{1}=0.0015$ , $\unicode[STIX]{x1D6FD}_{2}=0.016$ and $\unicode[STIX]{x1D6FD}_{3}=0.01$ .

3.5 Theoretical analysis of the self-similar profile of $\overline{\unicode[STIX]{x1D6F7}_{j}}$

In this section, the self-similar properties of the turbulent round jet and the closure developed in § 3.4 are used to develop a closed-form solution for the self-similar profile of the mean mass-weighted stream age, $\overline{\unicode[STIX]{x1D6F7}_{j}}$ . In agreement with previous studies (Hussein et al. 1994; Mi et al. 2001), figure 7 shows that, in the self-similar region of the flow, the centreline mean axial velocity $\overline{u_{c}}$ and mean mass fraction $\overline{\unicode[STIX]{x1D709}_{j,c}}$ vary in inverse proportion to the distance $x^{\prime }$ downstream from the virtual origin at $x_{0}$ . Consequently the jet spreads linearly (Pope 2000) and exhibits a constant entrainment rate (Ricou & Spalding 1961). The axial and radial mean velocities and the mean mass fraction may then be expressed as

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{u}=\frac{A}{x^{\prime }}f(\unicode[STIX]{x1D702}), & \displaystyle\end{eqnarray}$$

(3.17) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{v}=\frac{A}{x^{\prime }}g(\unicode[STIX]{x1D702}) & \displaystyle\end{eqnarray}$$

and

(3.18) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D709}_{j}}=\frac{B}{x^{\prime }}h(\unicode[STIX]{x1D702}),\end{eqnarray}$$

where $f(\unicode[STIX]{x1D702})$ , $g(\unicode[STIX]{x1D702})$ and $h(\unicode[STIX]{x1D702})$ are the self-similar shape functions for $\overline{u}$ , $\overline{v}$ and $\overline{\unicode[STIX]{x1D709}_{j}}$ respectively.

The mean mass-weighted stream age is assumed to have an $n$ -polynomial axial dependence with a coefficient $\unicode[STIX]{x1D6FC}$ , and its self-similar shape function is denoted by $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ :

(3.19) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D6F7}_{j}}=\unicode[STIX]{x1D6FC}(x^{\prime })^{n}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702}).\end{eqnarray}$$

According to figure 8, the radial turbulent flux of $\unicode[STIX]{x1D6F7}_{j}$ is also self-similar when normalised by the centreline values. Combining (3.17) and (3.19), $\overline{v^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }}$ is written as

(3.20) $$\begin{eqnarray}\overline{v^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }}=\overline{u_{c}}\overline{\unicode[STIX]{x1D6F7}_{j,c}}\unicode[STIX]{x1D703}(\unicode[STIX]{x1D702})=A\unicode[STIX]{x1D6FC}(x^{\prime })^{n-1}\unicode[STIX]{x1D703}(\unicode[STIX]{x1D702}),\end{eqnarray}$$

where $\unicode[STIX]{x1D703}(\unicode[STIX]{x1D702})$ is the shape function for $\overline{v^{\prime }\unicode[STIX]{x1D6F7}^{\prime }}$ .

Neglecting the axial turbulent flux and molecular transport, equation (3.10) for the Reynolds-averaged transport of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ can be rewritten as

(3.21) $$\begin{eqnarray}\overline{u}\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D6F7}_{j}}}{\unicode[STIX]{x2202}x}+\overline{v}\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D6F7}_{j}}}{\unicode[STIX]{x2202}r}=-\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r\overline{v^{\prime }{\unicode[STIX]{x1D6F7}_{j}}^{\prime }})+\overline{\unicode[STIX]{x1D709}_{j}}.\end{eqnarray}$$

Substituting (3.11), (3.12) and (3.15) into (3.21) and transforming the equation with respect to $\unicode[STIX]{x1D702}$ leads to

(3.22) $$\begin{eqnarray}x^{n-1}\left[n\,f(\unicode[STIX]{x1D702})\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})-\unicode[STIX]{x1D702}\,f(\unicode[STIX]{x1D702})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+g(\unicode[STIX]{x1D702})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\right]=-x^{n-1}\frac{1}{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}[\unicode[STIX]{x1D702}\unicode[STIX]{x1D703}(\unicode[STIX]{x1D702})]+\frac{B}{A\unicode[STIX]{x1D6FC}}h(\unicode[STIX]{x1D702}).\end{eqnarray}$$

Given that the flow is self-similar, the solution of (3.22) is required to be a function of $\unicode[STIX]{x1D702}$ only. Since the coefficients $A$ , $B$ and $\unicode[STIX]{x1D6FC}$ have no $x$ -dependence, $n$ must equal unity. This linear growth of $\overline{\unicode[STIX]{x1D6F7}}$ in the axial direction is consistent with the result shown in figure 7(a).

Next, this equation is closed by incorporating Model 3 from § 3.4, which is rewritten in the self-similar form as

(3.23) $$\begin{eqnarray}\text{Model 1:}~\unicode[STIX]{x1D703}(\unicode[STIX]{x1D702})=-\unicode[STIX]{x1D6FD}_{1}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}.\end{eqnarray}$$

Substituting (3.23) into (3.22) leads to

(3.24) $$\begin{eqnarray}f(\unicode[STIX]{x1D702})\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})-\unicode[STIX]{x1D702}f(\unicode[STIX]{x1D702})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+g(\unicode[STIX]{x1D702})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=\frac{\unicode[STIX]{x1D6FD}_{1}}{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\left(\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\right)+\frac{B}{A\unicode[STIX]{x1D6FC}}h(\unicode[STIX]{x1D702}),\end{eqnarray}$$

which is rewritten in the self-similar form as

(3.25) $$\begin{eqnarray}\text{Model 3:}~\unicode[STIX]{x1D703}(\unicode[STIX]{x1D702})=-\unicode[STIX]{x1D6FD}_{3}\left(\frac{1}{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D702}g(\unicode[STIX]{x1D702}))}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}-\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})\right).\end{eqnarray}$$

Substituting (3.25) into (3.22) leads to

(3.26) $$\begin{eqnarray}f(\unicode[STIX]{x1D702})\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})-\unicode[STIX]{x1D702}\,f(\unicode[STIX]{x1D702})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+g(\unicode[STIX]{x1D702})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}=\frac{\unicode[STIX]{x1D6FD}_{3}}{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}(\unicode[STIX]{x1D702}g(\unicode[STIX]{x1D702}))-\unicode[STIX]{x1D702}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})\right]+\frac{B}{A\unicode[STIX]{x1D6FC}}h(\unicode[STIX]{x1D702}),\end{eqnarray}$$

which is a closed-form equation for $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ , the shape function for the mean mass-weighted stream age $\overline{\unicode[STIX]{x1D6F7}_{j}}$ . This equation is an ordinary differential equation of $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ , so once $f(\unicode[STIX]{x1D702})$ , $g(\unicode[STIX]{x1D702})$ , $h(\unicode[STIX]{x1D702})$ and a single point of $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ are known, the entire solution of $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ can be calculated.

Figure 11. Comparison of the mass-weighted stream-age shape function profiles over the scaled radius between the current DNS data and solutions from (a) (3.24) and (b) (3.26). The parameters used are $\unicode[STIX]{x1D6FD}_{1}=0.0038$ , $\unicode[STIX]{x1D6FD}_{3}=0.01$ and $B/(A\unicode[STIX]{x1D6FC})=1.33$ .

Figure 12. (a) The skewness and (b) the kurtosis of mass-weighted stream age over the scaled radius.

Figure 11 shows the comparison of $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ from the DNS data and $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ evaluated from (3.24) and (3.26). In addition, the axial velocity shape function $f(\unicode[STIX]{x1D702})$ is shown for reference. The evaluation of the equations uses DNS data for $f(\unicode[STIX]{x1D702})$ , $g(\unicode[STIX]{x1D702})$ , $h(\unicode[STIX]{x1D702})$ and a single point of $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ , and then (3.24) and (3.26) are solved numerically. Note that due to the difficulties associated with $\unicode[STIX]{x1D702}=0$ (i.e. $1/\unicode[STIX]{x1D702}$ on the right-hand side of the equations), the $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ profiles are evaluated starting away from the centreline at $\unicode[STIX]{x1D702}=0.001$ and $\unicode[STIX]{x1D702}=0.04$ for (3.24) and (3.26), respectively. For (3.26), the slope of the calculated $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ at $\unicode[STIX]{x1D702}=0.04$ does not necessarily match the DNS data as the solutions are obtained by imposing a single value of $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ at a single point  $\unicode[STIX]{x1D702}$ . In addition, the solutions are very sensitive to small changes in the parameters $\unicode[STIX]{x1D6FD}_{1}$ , $\unicode[STIX]{x1D6FD}_{3}$ and $B/A\unicode[STIX]{x1D6FC}$ . Nevertheless, the comparison shows good agreement between the DNS and $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D702})$ calculated using Models 1 and 3.

3.6 Statistical distributions of mass-weighted stream age

Figure 12 shows the radial profiles of the skewness and the kurtosis of $\unicode[STIX]{x1D6F7}_{j}$ at various locations in the self-similar region of the jet. Overall, the skewness and kurtosis profiles show similar trends to those for the mass fraction shown in figure 4. Near $\unicode[STIX]{x1D702}=0$ , the weak negative skewness indicates that the distribution has a slightly longer tail on the left – making large negative fluctuations of $\unicode[STIX]{x1D6F7}_{j}$ more likely than large positive fluctuations. Away from $\unicode[STIX]{x1D702}=0$ the skewness increases parabolically with radius and switches to positive values, indicating that the distribution of $\unicode[STIX]{x1D6F7}_{j}$ has a relatively high probability for $\unicode[STIX]{x1D6F7}_{j}$ values several standard deviations greater than the mean in the radially outer portion of the flow.

The kurtosis takes values of less than 5 at the jet centreline, which, combined with the low magnitude of the skewness at the jet centreline, indicates that the $\unicode[STIX]{x1D6F7}_{j}$ distribution is approximately Gaussian near the centreline. The kurtosis values remain similar until $\unicode[STIX]{x1D702}=0.2$ , beyond which they increase rapidly, corresponding to a switch to bi-modal probability density function shapes.

Figure 13(a,b) shows the spectra of $\unicode[STIX]{x1D6F7}_{j}$ over the wavenumber. In figure 13(a), the Batchelor spectrum is calculated in the same way as calculating $E_{11}$ . The scaling of $k^{-5/3}$ is also observed over the range $k=3{-}30$ , which is similar to the velocity spectrum as shown in figure 2. Figure 13(b) shows the power spectrum by direct Fourier transform of the time series of $\unicode[STIX]{x1D6F7}_{j}$ at $x/D=15$ and $r/D=1.1$ . In the range $\unicode[STIX]{x1D714}=3{-}10$ , $\unicode[STIX]{x1D714}^{-4/3}$ scaling is observed. The spectrum shows that there are no particular narrow band fluctuations in  $\unicode[STIX]{x1D6F7}_{j}$ .

Figure 13. (a) Energy spectrum of mass-weighted stream age and (b) power spectrum by direct Fourier transform of the mass-weighted stream age at $x/D=15$ and $r/D=1.1$ .

Figure 14 shows the probability density function (pdf) of mass-weighted stream age on the centreline at axial positions in the range $x/D=5{-}35$ . The centreline $\unicode[STIX]{x1D6F7}_{j}$ pdf remains mono-modal but its peak value and variance evolve in the axial direction. Standardising the $\unicode[STIX]{x1D6F7}_{j}$ pdf by subtracting the mean and dividing by the standard deviation gives pdf shapes in the scaled sample space ( $(\unicode[STIX]{x1D6F7}_{j}-\overline{\unicode[STIX]{x1D6F7}_{j}})/\unicode[STIX]{x1D6F7}_{j}^{rms}$ ) that are self-similar between axial locations, for each value of scaled radius. This collapse is consistent with the self-similarities of the $\unicode[STIX]{x1D6F7}_{j}$ moments demonstrated above.

Figure 14. The probability density function of mass-weighted stream age on the centreline at axial positions in the range $x/D=5{-}35$ .

Figure 15(a–d) plots the self-similar profiles of the standardised $\unicode[STIX]{x1D6F7}_{j}$ pdf at $\unicode[STIX]{x1D702}=0$ , 0.08, 0.16 and 0.20 respectively, showing the transition from an approximately Gaussian distribution with a slight negative skewness at the centreline, through an approximately symmetrical but bi-modal distribution at $\unicode[STIX]{x1D702}=0.16$ , to a strongly positively skewed distribution at $\unicode[STIX]{x1D702}=0.20$ . This transition of the $\unicode[STIX]{x1D6F7}_{j}$ pdf shape with radius is qualitatively similar to the behaviour of the jet fluid mass fraction pdf (Gampert et al. 2013), in which the transition of the pdf shape is attributed to the effect of external intermittency increasing with radius, with the peak at low values of $\unicode[STIX]{x1D709}_{j}$ or $\unicode[STIX]{x1D6F7}_{j}$ corresponding to instances of nominally non-turbulent ambient fluid engulfed by the jet, and the peak at higher values of $\unicode[STIX]{x1D709}_{j}$ or $\unicode[STIX]{x1D6F7}_{j}$ corresponding to turbulent mixture. The self-similarity of the $\unicode[STIX]{x1D6F7}_{j}$ pdf demonstrated here substantially simplifies the task of modelling the distribution of $\unicode[STIX]{x1D6F7}_{j}$ in turbulent flow.

Figure 15. The probability density function of mass-weighted stream age on the scaled sample space at (a)  $\unicode[STIX]{x1D702}=0$ , (b)  $\unicode[STIX]{x1D702}=0.08$ , (c)  $\unicode[STIX]{x1D702}=0.16$ and (d)  $\unicode[STIX]{x1D702}=0.20$ .

4 Conclusion

The statistical distribution of the residence time within a turbulent flow is investigated using direct numerical simulation of a turbulent round jet. The jet issues from a smooth pipe with a top-hat velocity profile and a jet Reynolds number of 7290. The instantaneous local residence time and mass fraction of fluid originating from the jet are simulated using transport equations for fluid age and for mass fraction. The product of the mass fraction and the fluid age – described as the mass-weighted stream age $\unicode[STIX]{x1D6F7}_{j}$ – is found to be statistically stationary, and analysis of this quantity leads to the following key conclusions:

1. (i) The radial variations of the $\unicode[STIX]{x1D6F7}_{j}$ statistics are self-similar downstream of an initial development region when normalised by centreline values – the first four moments of $\unicode[STIX]{x1D6F7}_{j}$ exhibit self-similar profiles, and the one-point probability density functions of $\unicode[STIX]{x1D6F7}_{j}$ from different axial positions show collapse when standardised and grouped by the scaled radius.

2. (ii) The spatial variations of the mean $\unicode[STIX]{x1D6F7}_{j}$ is influenced most by the counteracting effects of convection by the mean axial velocity and increase of the fluid age due to the passage of time. Radial transport by the mean radial velocity and by the radial turbulent flux provide the next most significant contributions to the transport of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ , whereas the molecular transport and the axial turbulent flux are negligible.

3. (iii) The radial turbulent flux of $\unicode[STIX]{x1D6F7}_{j}$ is self-similar and modelling is developed in terms of the self-similar Reynolds-averaged radial velocity and $\unicode[STIX]{x1D6F7}_{j}$ profiles, providing a closure for the Reynolds-averaged transport equation for  $\unicode[STIX]{x1D6F7}_{j}$ .

4. (iv) Consistency of the Reynolds-averaged equation for $\unicode[STIX]{x1D6F7}_{j}$ with the self-similar properties of the flow requires that the mean centreline value $\overline{\unicode[STIX]{x1D6F7}_{j}}$ increases linearly along the jet axis, and this behaviour is confirmed by the observations from the DNS.

5. (v) The modelled Reynolds-averaged transport equation for $\unicode[STIX]{x1D6F7}_{j}$ written in terms of the radial similarity coordinate $\unicode[STIX]{x1D702}$ provides a closed equation for the self-similar shape function of $\overline{\unicode[STIX]{x1D6F7}_{j}}$ , and the predicted shape function shows good agreement with the observations from the DNS.

The residence time distribution provides information that can be related to the extent of progress of kinetic processes in a flow. The self-similar properties of the residence time distribution that are demonstrated by this investigation simplify the task of modelling the residence time distribution, thereby facilitating the development of modelling based on residence time for kinetically limited processes, such as for autoignition in turbulent fuel jets.