Comparative study of turbulent round jet flows through direct numerical simulation at medium–high Reynolds numbers

Turbulent jets are essential in jet propulsion systems, where understanding and optimizing jet dynamics can improve engine performance and fuel efficiency. In addition, turbulent round jets have unique properties that make them useful for validating turbulence models and simulation techniques. In particular, these flows exhibit statistical self-similarity in the far-field, where radial profiles of different statistical quantities can be scaled to converge to a single similarity function. Round jet flows describe fluid flows ejected from a circular inlet entering an ambient fluid of the same fluid leading to turbulence and mixing (see Fig. 1).

The concept of universal self-similarity in turbulent round jets, widely discussed in the literature, suggests that this state should be reached asymptotically and independent of initial conditions (Townsend, 1976). However, theoretical analyses indicate that self-similar profiles may depend on initial conditions (George, 1989). Experimental studies have shown variability in self-similar profiles between different experiments (Wygnanski and Fiedler, 1969; Panchapakesan and Lumley, 1993; Hussein , 1994), often attributed to measurement inaccuracies or environmental disturbances.

Experimental studies have also investigated scalar properties in turbulent jets. Dowling and Dimotakis (1990) studied the turbulent concentration field at Re = 5000, 16 000, and 40 000 and showed nearly independent behavior with respect to the Reynolds number near the centerline for the scaled probability density function (PDF) of the jet fluid concentration. Antonia and Mi (1993) measured the average temperature dissipation using parallel cold wires at $Re=19000$ and a Péclet number of Pe = 83 and found the components of the mean temperature dissipation to be nearly equal, with slight deviations from isotropy.

Xu and Antonia (2002) compared jets exiting from a smooth contraction nozzle and a pipe at $Re=86000$⁠, demonstrating that a contraction jet reaches self-preservation earlier than a fully developed pipe jet. This is in contrast to Ferdman (2000), who found that a pipe jet reaches self-preservation earlier. Both studies found that a fully developed pipe profile results in a lower far-field decay rate compared to a top-hat velocity distribution.

Darisse et al. (2015) collected data on a slightly heated round jet at $Re=140 000$⁠, measuring pure and mixed moments up to the third order, allowing for the determination of the TKE and passive scalar transport balances. Different assumptions for the dissipation and pressure diffusion terms led to some discrepancies in the results.

The advent of supercomputers in the 1980s made direct numerical simulations (DNS) of fluid flows possible, providing high-fidelity representations of turbulent flows without relying on turbulence models. One of the earliest DNS of a turbulent round jet was performed by Boersma (1998), who studied the effect of inflow conditions on self-similarity up to $z/D=45$ at Re = 2400. However, they were unable to resolve the far-field due to limited computational resources, so the fluctuations have not yet fully reached self-similarity. As an extension of this work, Lubbers (2001) studied the self-similarity of a passive scalar concentration at Re = 2000 and a Schmidt number Sc = 1 in a box length $z/D=40$⁠. The statistics were extracted over $80D/Ub$ time units, where Ub refers to the bulk velocity at the inlet. The results show that the mean concentration in the far-field is self-similar. However, the root mean square of the concentration fluctuations is not self-similar.

Taub (2013) conducted a DNS at Re = 2000 at the same box length where he extracted statistics up to the third order, including the TKE terms and the terms of the Reynolds stress transport equations terms, directly. The data are compared extensively with previous studies. They also conclude that there are inconsistencies in the dissipation profiles due to different approximations in previous experimental studies, which may also be due to the Reynolds number dependence of dissipation (Bogey and Bailly, 2006).

In Babu and Mahesh (2005), a DNS of a turbulent jet at Re = 2400 and Sc = 1 is performed. The data are averaged over $1400D/Ub$ time units. The instantaneous radial profiles of the velocity and passive scalar exhibit similarities and alternate between “top-hat” and “triangle” profiles, both spatially and temporally. These profiles are influenced by entrainment from the free-stream into the jet, resulting in a mean Gaussian profile as a function of r. Diffusion-dominated regions are observed, occurring closer to the jet center and as “brush-like” regions near the jet edge. The width of these brush-like regions decreases with increasing Re, suggesting a transition in mixing behavior.

Babu and Mahesh (2005) performed a DNS of a turbulent jet at Re = 2400 and Sc = 1, observing instantaneous radial profiles of velocity and passive scalar that alternate between “top-hat” and “triangle” profiles, influenced by entrainment from the free-stream. Gilliland (2012) focused on scalar intermittency in a turbulent round jet through DNS at Re = 2400 and LES at $Re=68000$⁠, emphasizing the importance of external intermittency in scalar mixing.

Recently, we have conducted a DNS of a spatially evolving turbulent round jet flow at Re = 3500 (Nguyen and Oberlack, 2024a) with an additional passive scalar at Pr = 0.71 (Nguyen and Oberlack, 2024c). A fully developed turbulent pipe flow profile has been utilized as an inlet by running two simulations in parallel. High-quality statistics, comprising mean quantities, Reynolds stresses, heat fluxes, mixed moments, Reynolds stress budgets, TKE budgets, and PDFs of the axial velocity and passive scalar, were taken on-the-fly over $75 000D/Ub$ time units, achieving an excellent convergence of statistics. All statistics show self-similarity across the range of $z/D=25−65$⁠. The PDFs on the centerline show almost complete Gaussian behavior with strong tails appearing away from the centerline. In addition, it was found that a fully turbulent velocity profile as an inlet can have similar effects on dissipation as a larger Reynolds number. In a companion paper (Nguyen and Oberlack, 2024b), the data are used to validate a Lie symmetry analysis of a turbulent round jet flow, which revealed that the full momentum integral is crucial for the generation of turbulent scaling laws that allow for a variation of the decay behavior through the inlet condition, just as in George (1989). The symmetry-based turbulent scaling laws for pure hydrodynamics and the addition of passive scalars have been generalized in (Nguyen and Oberlack, 2024c).

The present study provides statistics of a spatially evolving turbulent round jet DNS at Re = 7000 under identical conditions of the DNS in Nguyen and Oberlack (2024a; 2024c) to ensure comparability. The statistics have been averaged over $30 000D/Ub$ time units. These results are compared with previous studies, particularly with Nguyen and Oberlack (2024a; 2024c), to investigate the effects of varying Reynolds numbers on turbulent jet flow characteristics and scalar properties.

In the study of fluid dynamics, the behavior of fluid flow is fundamentally described by the Navier–Stokes equations (NSEs). These equations, derived from the conservation of mass and momentum, provide a mathematical framework for understanding the motion of fluid substances in a wide range of scenarios. In the following, incompressibility is implied, which assumes that density variations are negligible, i.e., density is constant.

The NSEs (2) and (3) and the passive scalar transport equation (4) are solved directly using Nek5000, a computational fluid dynamics (CFD) code which has been developed by Fischer (2008). Nek5000 uses a spectral element method (Patera, 1984), that capitalizes on the accuracy of spectral methods, resulting in a high order of accuracy and a fast convergence of solution to the exact solution. The solver is designed to make efficient use of parallel computing architectures, using Message Passing Interface (MPI) for parallelization and allowing scalability to a large number of processors. In the DNS of the turbulent round jet, a polynomial order of N = 7 was used for the spatial discretization. This choice of polynomial order optimizes computational efficiency with accuracy in resolving turbulent structures over a range of length scales present in the flow. In addition, the implicit second-order backward differentiation formula (BDF2) is employed as a time-stepping scheme. The simulation used 96 000 cores and ran for approximately 50 × 106 core hours.

To efficiently derive statistical moments, a statistics toolbox (Massaro , 2024) has been used. This toolbox computes 44 variables, which can be used to calculate mean velocities, Reynolds stress, turbulent kinetic energy (TKE) budgets, and Reynolds stress budgets. Initially designed to compute statistical velocity moments up to the third order on-the-fly, it has been expanded to handle passive scalar and mixed moments for the present DNS. To obtain accurate statistics, the averages have been taken over $30 000D/Ub$ time units. The output for a single set of statistical data is approximately 415GB. Taking advantage of the spectral code, the statistics have been extracted at arbitrary points. In addition, azimuthal direction averaging was performed during post-processing, receiving statistics in the r and z directions.

The computational domain has been described in Nguyen and Oberlack (2024a) for the DNS of a turbulent round jet flow at Re = 3500. The overall mesh geometry has not been changed for the present DNS, and only the resolution has been adjusted for the larger Reynolds number. The computational domain is split into two subdomains. A fully developed turbulent pipe flow is utilized as an inlet condition for the turbulent round jet flow, which is run in parallel on one of the subdomains. The main computational domain, where the fully turbulent pipe flow velocity profile is interpolated on, incorporates the spatially developing turbulent round jet flow. To ensure feasibility and the ability to identify differences between the present simulation and the one in Nguyen and Oberlack (2024a), the Reynolds number Re = 7000 has been selected, which is twice as high as in the aforementioned work. The Prandtl number Pr = 0.71 for air has been chosen for the passive scalar field, which yields the Péclet number Pe = 4970. Table I summarizes the parameters of the present DNS and that of Nguyen and Oberlack (2024a).

The computational domain of the pipe flow extends 5D axially with a diameter D. A cross-section of this domain has 756 cells, while in the axial direction there are 72 cells. Near the wall, the usual refinement has been applied. With N = 7, the computational domain of the pipe flow contains around 74 × 106 grid points. A quarter of the cross-section is showcased in Fig. 2.

At the wall, the boundary conditions (BCs) are no-slip and impermeable wall, and at $z/D=−5$ and $z/D=0$⁠, a periodic BC has been employed. To generate turbulence, a small perturbation, being a sine function at 10% of the maximum value of the laminar pipe flow profile, has been superposed to a laminar pipe flow profile. All quantities have been non-dimensionalized by setting the pressure gradient and viscosity such that the bulk velocity is 1 and the Reynolds number is constant. This simulation is run until the velocity profile is statistically stationary, reaching the state depicted in Fig. 3. It can be observed that the large scales are dominant in the center of the flow.

The computational domain of the round jet flow is a truncated cone with a diameter of 4D at $z/D=0$ and 64D at $z/D=75$⁠. The diameter increases linearly in z-direction to capture the spreading of the spatially evolving round jet flow. Additionally, it is ensured that the lateral boundaries do not influence the jet flow.

The cross-section of this mesh encases the pipe flow mesh in the center of an annulus mesh (see Fig. 4). The mesh in Fig. 4 is extruded and spreads linearly until $z/D=75$⁠. The linear spreading factors in that the length scales increase further away from the orifice. Therefore, the mesh can be coarser in the far-field compared to the near-field. In z-direction, the mesh is stretched geometrically with a factor of 1.003 which leads to 275 cells in the axial direction. Since the ambient region does not need a fine mesh, the annulus mesh also stretches geometrically with a factor of 1.06 which results in 25 cells radially. In sum, the jet flow mesh comprises 900 000 cells, which, under consideration of N = 7, yields 1.24 × 109 degrees of freedom.

The boundary conditions (BCs) at the orifice interpolate the velocity field at $z/D=0$ from the pipe mesh onto the jet flow mesh at $z/D=0$⁠. On the annulus, the BC is a no-slip and impermeable wall with the option of inducing a co-flow. All the other boundaries employ the open boundary as described in Dong (2014), which addresses instability issues of turbulent round jet flow simulations, among others. At the orifice, the passive scalar is set to Θ = 1 as a Dirichlet BC. The annulus at $z/D=0$ is also set as a no-slip and impermeable wall, while the boundaries use a similar open boundary for passive scalars (Liu , 2020). The BCs of the velocity and passive scalar field are summarized in Table II.

The time step has been set with a target Courant number of 2.5. The initial time step satisfies this condition with $Δt=2.5·10−3$⁠. To ensure the stability of the simulation, the time step is being adjusted dynamically such that the Courant number is constant.

In Nguyen and Oberlack (2024a), the DNS at Re = 3500 has been extensively validated against data from previous studies, while in Nguyen and Oberlack (2024c) the study of the additional passive scalar is discussed. Therefore, the focus here is the comparison of the present DNS to Nguyen and Oberlack (2024a; 2024c).

The inverse of the mean axial centerline velocity can be viewed in Fig. 5. The data are compared to previous experimental and numerical studies. Further, an inset figure shows a magnified view of the near-field region. It can be observed that the potential core, which is defined as the region up to where the mean axial centerline velocity reaches $95%$ (Jambunathan , 1992), is larger for the present DNS compared to the one in Nguyen and Oberlack (2024a). However, they are both smaller than the other studies with the length being $z/D=3–4$ compared to $z/D=6–9.5$ which amounts to a decrease of 30%–70%. This can be attributed to the fully turbulent pipe flow inlet, which generates initial turbulence intensities to be viewed later. The different sizes of the potential core may be explained by the turbulent pipe flow profile. At higher Reynolds numbers, the steep near-wall gradients of turbulent pipe flows increase. Therefore, stronger shear layers with sharper gradients result near the orifice of the jet. As a consequence, the “coalescence” of the shear layers to form the jet profile is delayed in the z-direction, resulting in longer potential cores. The top-hat profiles of previous studies at lower Reynolds numbers have sharper gradients near the orifice compared to the present Re = 3500 and Re = 7000 cases, leading to even larger potential cores. In the far-field, it can be observed that the inverse velocity dips slightly. This is due to the boundary effects. Therefore, this region is ignored when adjusting the parameters.

The parameters of Eq. (6) have been determined by minimizing the sum of the squared residuals relative to the DNS data on the interval from $z/D=15–65$⁠. The far-field between $z/D=65–75$ has been ignored so that boundary effects are left out. The present DNS decays at $Bu=5.25$ compared to $Bu=5.15$ in Nguyen and Oberlack (2024a) and $z0=2$ compared to $z0=0$⁠, respectively. The shift of the virtual origin influences the self-similarity for all quantities in the following. Due to the larger virtual origin, self-similarity sets in further away from the orifice in the present DNS than the DNS at the smaller Reynolds number. However, even with the higher Reynolds number, where a lengthening of the potential core is observed, it remains significantly shorter than in previous studies, indicating a fundamental change in flow behavior due to the initial turbulent conditions. This earlier onset can result in more predictable flow, which is beneficial for engineering applications where consistent flow behavior is desired and Reynolds numbers are typically very high.

Similarly, the inverse of the mean passive scalar is depicted in Fig. 6, where again an inset figure shows a magnified view of the near-field region. Also here, the potential core is larger for the present DNS compared to the previous DNS at Re = 3500 in Nguyen and Oberlack (2024c) and smaller than other studies. The size difference compared to other studies is not as prominent meaning that the scalar field is less sensitive to the initial turbulent conditions than the velocity field. Still, the potential core is affected by the Reynolds number independent of thermal diffusivity, which has remained constant at both Reynolds numbers. The parameters of Eq. (6) have also been determined by minimizing the sum of the squared residuals showing that $BΘ=4.9$ compared to $BΘ=5.0$ from Nguyen and Oberlack (2024c). The virtual origin at $z0=2.7$ is close to the virtual origin of the axial velocity determined earlier. Also here, the virtual origin moves to a larger value due to the increased Reynolds number. Figures 5 and 6 show highly converged statistics of the present DNS compared to previous studies.

The invariant $Uz¯̃(η=0)=BuDU¯0=U¯z(r=0,z)(z–z0)$ is depicted in Fig. 9. It shows that self-similarity is observed well into the far-field for larger z. It is compared to previous studies showcasing a highly converged behavior and is slightly larger than that of the DNS at Re = 3500. This implies a faster decay rate at a larger Reynolds number. Additionally, self-similarity is reached at a later stage than in the DNS at Re = 3500.

In this section, the components of the turbulence intensity on the centerline are compared with various studies.

Figures 10–12 show the radial and axial components of the turbulence intensity and the normalized root mean square (rms) of the passive scalar fluctuation, respectively. The regions in these figures have to be distinguished in the following. The region in the far-field near the boundary shall be ignored due to numerical effects. In the far-field region, where self-similarity occurs, it can be observed that the turbulence intensity is larger for higher Reynolds numbers. This is to be expected and is well known. It is noted that self-similarity occurs further downstream compared to first order moments which has also been observed in Wygnanski and Fiedler (1969).

In the near-field region, similar to Nguyen and Oberlack (2024a), in Figs. 10 and 11, the components have nonzero turbulence intensity. This is due to the fully developed turbulent pipe velocity profile, which introduces turbulence. This leads to the smaller transition region, i.e., a smaller potential core, which has been observed in Figs. 5 and 6 (see also Bogey and Bailly, 2009). In contrast, the other studies depicted in the figure have used top-hat profiles instead. Taub (2013) has introduced small nonphysical perturbations onto that profile, while Bogey and Bailly (2006) utilized divergence free vortex rings. Both did not initiate fully turbulent behavior compared to the present inlet condition. Compared to Nguyen and Oberlack (2024a), where the initial maxima are at $z/D=8$ for both turbulence intensity components, this is presently not observed, leading to a smoother transition to self-similarity. At the higher Reynolds numbers, the dissipation of the TKE is larger (see later in Fig. 17). This dissipation may lead to a stronger reduction in the turbulence intensity peaks as localized regions of high turbulence are more quickly dissipated and distributed throughout the flow. In addition, since the energy is more evenly distributed over a wider range of scales, this may also imply an alteration in noise generation and mixing in practical applications. Further, the potential core length is larger for the Re = 7000 case (see Fig. 5) which suggests a longer duration needed for turbulence to develop downstream on the centerline since the thinner shear layers at higher Reynolds numbers that develop on the edges of the jet connect more smoothly. This might also contribute to the more gradual development of the turbulence intensities. It can be observed that the components of the turbulence intensity are marginally larger than at Re = 3500 in the self-similar region when the jet is fully developed.

In Fig. 12, where the rms of the passive scalar fluctuation $θ2¯$ is shown, the initial maximum observed in the DNS data at Re = 3500 (Nguyen and Oberlack, 2024c) at $z/D=6$ is even more dominant than in the turbulence intensities and is smoothed out for the present DNS data. The reasons have been mentioned above. It follows the experimental results of Birch (1978) at $Re=16000$ quite well before reaching a constant value at $RΘΘ/Θc=0.22$⁠, while the data of Birch (1978) show an increasing trend. In the far-field, the behavior of the present DNS reaches a higher value compared to the DNS at Re = 3500.

The passive scalar fluctuation is presented in Fig. 14 at different distances z/D from the orifice. The profiles collapse later than $z/D=30$ while the value on the centerline is at $θ2¯/Θc2=0.05$ which is slightly larger than in Nguyen and Oberlack (2024c) at $θ2¯/Θc2=0.048$⁠. However, both simulations exhibit the same value on the off-axis peak. Darisse et al. (2015) reported $θ2¯/Θc2=0.04$ on the centerline and $θ2¯/Θc2=0.05$ at the off-axis peak. The difference on the off-axis peak might be due to the different inlet conditions, as explained for the Reynolds stress.

The turbulent heat fluxes in Fig. 15 reveal a maximum value of $urθ¯max/(Uz,cΘc)=2.2$ and $uzθ¯max/(Uz,cΘc)=3.2$⁠. Further, the value at η = 0 is at $uzθ¯/(Uz,cΘc)=0.027$⁠. Both maximum values are also observed in the DNS at Re = 3500 (Nguyen and Oberlack, 2024c). Only the centerline value of the axial heat flux is slightly larger at $uzθ¯/(Uz,cΘc)=0.027$ compared to $uzθ¯/(Uz,cΘc)=0.025$ at Re = 3500. Darisse et al. (2015) have conducted their experiment at $Re=16000$ and received $uzθ¯/(Uz,cΘc)=0.024$ on the centerline, which suggests that a turbulent velocity inlet might increase the axial turbulent heat flux on the centerline. Additionally, a larger Reynolds number increases the maximum of the axial heat flux.

The dissipation of the present DNS is compared to the DNS at Re = 3500 and various previous studies in Fig. 17. It can be observed that the qualitative behavior of the profiles is similar. Further, the dissipation of the present DNS is larger compared to the DNS at Re = 3500, implying growth of dissipation with increasing Reynolds numbers.

Further, the third-order mixed statistics $RiΘΘ$ have been extracted from the present DNS. The radial profiles are showcased in Fig. 18. The maximum of these moments at $urθ2¯max/(Uz,cΘc2)=2.7·10–3$ and $uzθ2¯max/(Uz,cΘc2)=2·10–3$ are found in the present DNS and also in the comparative DNS at Re = 3500 (Nguyen and Oberlack, 2024c). The minimum are also equivalent in both at $urθ2¯min/(Uz,cΘc2)=–1.6·10–3$ and $uzθ2¯min/(Uz,cΘc2)=–0.9·10–3$⁠, although $uzθ2¯$ is not fully converged near the centerline. At the centerline, $uzθ2¯/(Uz,cΘc2)=–0.6·10–3$⁠, which is slightly larger in magnitude than $uzθ2¯/(Uz,cΘc2)=–0.5·10–3$ at Re = 3500 (Nguyen and Oberlack, 2024c). Similar to the axial heat flux mentioned above, a larger Reynolds number increases the maximum of the moment correlated with the axial velocity $RzΘΘ$⁠.

PDFs have been generated for the axial velocity Uz and the passive scalar Θ by collecting the quantity at various points every $0.025D/Ub$ over a span of $30 000D/Ub$⁠. The data are then sorted into bins, where data from specific r and z values are grouped into the same PDF, contributing to the smoothness of the resulting PDF. After binning, the probability density for each bin is calculated by counting the number of data points within each bin. To ensure the PDF reflects true probabilities, these densities are normalized so they sum to 1.

On the centerline, both the PDF f of axial velocity Uz and the passive scalar Θ exhibit self-similar properties. In Fig. 19, the PDF of both, normalized by the centerline value, can be viewed. Both collapse well and follow a centered Gaussian distribution, indicated by a black line. A Gaussian distribution is only dependent on the mean, which, due to the normalization, is 1, and the standard deviation. The standard deviation of both quantities is given by $uz2¯$ and $θ2¯$⁠, respectively. In Figs. 11 and 12, we find that in the far-field, if the area close to the outflow boundary is ignored, the higher Reynolds number case has a slightly larger standard deviation. This differentiates the PDF of Uz and Θ on the centerline for different Reynolds numbers. This is expected since it is well known that at higher Reynolds numbers, the turbulence intensities typically increase (Kumar and Sharma, 2024).

In Fig. 20, the PDFs f of the normalized Uz and Θ are provided for three distances from the orifice for varying η. Just as in Nguyen and Oberlack (2024a) and Nguyen and Oberlack (2024c), these deviate from a Gaussian distribution with increasing η and approach a delta distribution due to the non-turbulent ambient region. The profiles collapse for both quantities. Their skewness S and kurtosis K are exhibited in Fig. 21 for the axial velocity and Fig. 22 for the passive scalar. It is noted that a Gaussian distribution is described by S = 0 and K = 3 and is denoted by a dashed line in these figures. The S and K profiles collapse well for the present DNS and the one at Re = 3500. Just as in Nguyen and Oberlack (2024a), the Uz based skewness in Fig. 21 increases with K = 3 until $η=0.09$ before deviating from a Gaussian distribution. For Θ in Fig. 22, at η = 0, the kurtosis is $K≈3$ while it is slightly super-Gaussian, i.e., K > 3, for Re = 3500 (Nguyen and Oberlack, 2024c). The PDF of the passive scalar is slightly negatively skewed on the centerline compared to the axial velocity PDF. This may be due to the reason that velocities can be any real number, while passive scalars can only be positive.

In this work, a large-scale DNS of a spatially evolving turbulent round jet at Re = 7000 has been conducted to retrieve velocity and passive scalar statistics. The present DNS is directly compared to a DNS at Re = 3500 (Nguyen and Oberlack, 2024a; 2024c), whose conditions only differ by their Reynolds number and where both utilize a turbulent pipe flow profile as an inlet condition. Therefore, the numerical box also extends 75D axially and up to 65D radially. This allows the investigation of how turbulence characteristics evolve with increasing Reynolds number, where other effects can be excluded. A vast collection of statistical quantities has been extracted from the present DNS by averaging over $30 000D/Ub$ time units ranging from first to tenth order moments. From these, the mean, Reynolds stress, heat flux, Reynolds stress budgets, TKE budgets, third-order mixed moments, and PDFs of the axial velocity and passive scalar are discussed and compared to previous and the comparative study. This comprehensive data set, which is openly available, provides a reliable basis for the research community to conduct further analysis and validation of turbulence models and theories.

The results show that the turbulent round jet achieves self-similarity for several statistical quantities, although at a later stage compared to the DNS at Re = 3500. The present study places particular emphasis on the impact of using a fully developed turbulent pipe flow profile as the inlet condition instead of top-hat profiles or laminar pipe flow profiles. Here, it is demonstrated how the turbulent pipe flow profile affects the flow dependent on the Reynolds number. The virtual origin shifts downstream when increasing the Reynolds number. Additionally, the magnitude of the Reynolds normal stresses is increased at higher Reynolds numbers, while the Reynolds shear stress is unaffected. Using a fully turbulent pipe velocity profile as the inlet reduces the size of the potential core for both Reynolds numbers, resulting in an earlier onset of decay with a larger potential core as the Reynolds number increases. In the near-field, a larger Reynolds number smooths the initial maximum visible in the rms of the passive scalar fluctuation and turbulence intensities at Re = 3500. On the centerline, the second and third-order mixed moments are large compared to previous studies, which is attributed to the turbulent velocity inlet. In addition, the axial second and third moments are increased at a larger Reynolds number when comparing the DNS at Re = 3500 and Re = 7000. The TKE budgets show larger production and dissipation for the present DNS compared to the DNS at Re = 3500. A Gaussian distribution describes the PDF of the axial velocity and the passive scalar at the centerline well. At higher Reynolds numbers, the standard deviation increases, which is expected since it is known that turbulence intensities, i.e., the standard deviations, typically increase at higher Reynolds numbers. At the centerline, the passive scalar PDF is slightly negatively skewed compared to the axial velocity PDF.

The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project under the Project No. pn73fu by providing computing time on the GCS Supercomputer SUPERMUC-NG at Leibniz Supercomputing Centre (www.lrz.de). Further, C.T.N. acknowledges the funding from the Studienstiftung des deutschen Volkes (Academic Scholarship Foundation) and M.O. gratefully acknowledges partial funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—SPP 2410 Hyperbolic Balance Laws in Fluid Mechanics: Complexity, Scales, Randomness (CoScaRa), within the Project “Approximation Methods for Statistical Conservation Laws of Hyperbolically Dominated Flow” under Project No. 526024901.

The authors have no conflicts to disclose.

Cat Tuong Nguyen: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). Martin Oberlack: Methodology (equal); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are openly available in the TUdatalib Repository of TU Darmstadt at https://doi.org/10.48328/tudatalib-1572.

Further, the Reynolds stress budgets have been directly calculated from the DNS and are displayed in Fig. 23 at $z/D=35$⁠. They are received in Cartesian coordinates and have been transformed to their cylindrical form and with subsequent spatially averaging in $φ$-direction. Although the data are available at various distances from the orifice, they are only shown here at $z/D=35$ due to readability.