1. Introduction
Hydrogen (H
$_2$
) has been considered a promising solution to improve the performance and reduce emissions of a range of energy and propulsion systems, such as off-road internal combustion engines (ICEs) and industrial gas turbines (Akal, Öztuna & Büyükakın Reference Akal, Öztuna and Büyükakın2020; Hwang, Maharjan & Cho Reference Hwang, Maharjan and Cho2023; Manigandan et al. Reference Manigandan, Praveenkumar, Ryu, Verma and Pugazhendhi2023). In combustion systems such as H
$_2$
ICEs with port fuel injection (PFI), hydrogen is injected into and mixed with air through a jet in crossflow (JICF) configuration. The mixing process is key to achieving the desired level of fuel–air premixing that can strongly impact engine performance and emissions characteristics. Hence it is important to better understand the mixing process of H
$_2$
JICF systems, which would enable improved injector and engine designs towards next-generation hydrogen engines.
Computational fluid dynamics (CFD) provides a unique tool for characterising fuel–air mixing processes in JICFs. Previous numerical studies on JICFs have demonstrated that both direct numerical simulations (DNS), which resolve all turbulent scales, and large-eddy simulations (LES), which resolve large turbulent scales with sufficiently small grid size, are able to predict the key physical structures in JICF and scalar fields accurately (Muppidi & Mahesh Reference Muppidi and Mahesh2007, Reference Muppidi and Mahesh2008; Bodart et al. Reference Bodart, Coletti, Bermejo-Moreno and Eaton2013; Ryan et al. Reference Ryan, Bodart, Folkersma, Elkins and Eaton2017; Zhang & Yang Reference Zhang and Yang2017; Uyanwaththa et al. Reference Uyanwaththa, Malalasekera, Hargrave and Dubal2019). For example, Muppidi & Mahesh (2007, Reference Muppidi and Mahesh2008) performed DNS studies for a transverse JICF, and they showed that DNS predictions agreed well with experimental measurement of the mean velocity field, turbulent intensities and scalar fields. Bodart et al. (Reference Bodart, Coletti, Bermejo-Moreno and Eaton2013) conducted an LES study for an inclined JICF, and showed good agreement between their LES results and experimental data on velocity and scalar fields. However, from a practical standpoint, both LES and particularly DNS are still computationally too expensive to be applied directly to complex industrial problems. Therefore, Reynolds-averaged Navier–Stokes (RANS) based approaches, due to their low computational cost, remain the predominant model used in the iterative engine design process.
In the RANS framework, an additional model for the turbulent scalar fluxes is required to close the scalar transport equations, which will significantly influence the accuracy of prediction for the mixing process. Among various existing turbulent scalar flux models, the simplest and most commonly used is the gradient diffusion hypothesis (GDH), which models the turbulent scalar fluxes as
\begin{align} \overline {u_i^{\prime }C^{\prime }} = - D_t \frac {\partial \overline {C}}{\partial x_i}, \end{align}
where
$u_i$
is the
$i$
th component of velocity, the over-bar represents the time-averaging operation, and
$C$
is an arbitrary scalar of interest, such as species mass fraction or temperature. The turbulent diffusivity
$D_t$
is usually calculated through the Reynolds analogy with a fixed turbulent Schmidt number
$ \textit{Sc}_t$
(or Prandtl number
$ \textit{Pr}_t$
if
$C$
is temperature) as
\begin{align} D_t = \nu _t/Sc_t , \end{align}
where
$\nu _t$
is the eddy viscosity. Clearly, the accuracy of the mixing model relies on two terms,
$\nu _t$
and
$ \textit{Sc}_t$
. Typically,
$ \textit{Sc}_t$
is set in the range
$0.7 \leqslant Sc_t \leqslant 1.0$
, while in many cases,
$ \textit{Sc}_t$
is considered a tunable model constant that is determined case by case. However, many studies (Alvarez, Jones & Seoud Reference Alvarez, Jones and Seoud1993; He et al. Reference He, Guo, Hsu, Brankovic, Syed and Liu1999; Galeazzo et al. Reference Galeazzo, Donnert, Habisreuther, Zarzalis, Valdes and Krebs2011; Ivanova, Noll & Aigner Reference Ivanova, Noll and Aigner2011) have found that this typical range of
$ \textit{Sc}_t$
tends to under-predict the mixing process in JICF problems, and lower values of
$ \textit{Sc}_t$
in the range 0.2–0.3 were recommended by trial and error to achieve better agreement with experimental data. In addition, Ivanova et al. (Reference Ivanova, Noll and Aigner2011) found that the RANS model also under-predicted the turbulent kinetic energy in JICF, indicating a possible under-prediction of turbulent viscosity
$\nu _t$
. In this context, the extremely low
$ \textit{Sc}_t$
recommended (i.e.
$ \textit{Sc}_t=0.2{-}0.3$
) might be just an artefact associated with the under-prediction in
$\nu _t$
. In particular, Ivanova, Noll & Aigner (Reference Ivanova, Noll and Aigner2013) found that the optimal
$ \textit{Sc}_t$
value depends on the RANS turbulence model used. When employing the shear-stress transport model, which under-predicted the turbulence level, the optimal
$ \textit{Sc}_t$
was 0.2. In contrast, when a curvature correction modification was applied to properly reproduce the turbulence, an optimal value of 0.5 was found. More interestingly, this optimal
$ \textit{Sc}_t$
of 0.5 matches the
$ \textit{Sc}_t$
derived from their LES results. A similar optimal
$ \textit{Sc}_t$
of 0.6 was also reported in Ling et al. (Reference Ling, Ryan, Bodart and Eaton2016) by extracting
$ \textit{Sc}_t$
from LES results conducted by Bodart et al. (Reference Bodart, Coletti, Bermejo-Moreno and Eaton2013). However, a smaller optimal
$ \textit{Sc}_t$
of 0.35 was observed in the LES study conducted by Ryan et al. (Reference Ryan, Bodart, Folkersma, Elkins and Eaton2017) for a skewed JICF. It is therefore evident that when using the GDH, accurate descriptions of both the turbulence model and the Schmidt number are critical.
Another key simplification in the GDH is that it assumes an isotropic turbulent diffusivity, which was also believed to cause errors in mixing prediction for inhomogeneous flows (Lakehal Reference Lakehal2002; Liu, Zhu & Bai Reference Liu, Zhu and Bai2008). To allow for anisotropy of
$D_t$
, Daly & Harlow (Reference Daly and Harlow1970) proposed the generalised gradient diffusion hypothesis (GGDH) by including Reynolds stress components in the definition of
$D_t$
. Abe & Suga (Reference Abe and Suga2001) further proposed the higher-order generalised gradient diffusion hypothesis (HOGGDH) to better predict the streamwise component in fully developed turbulent channel flows by using quadratic products of the Reynolds stresses. Ling et al. (Reference Ling, Ryan, Bodart and Eaton2016) evaluated the performance of the GDH, GGDH and HOGGDH in the RANS calculation, with the velocity field provided by LES for an inclined JICF. They found that switching from the GDH to higher-order models would improve scalar concentration predictions to a certain extent, but the effect was significantly smaller compared to the tuning of the model parameters. Similar observations were made in the LES study by Ryan et al. (Reference Ryan, Bodart, Folkersma, Elkins and Eaton2017), which highlighted the importance of model constant optimisation.
In all the above-mentioned JICF studies, the density ratio between the jet and crossflow was set to unity with constant molecular transport properties, and the mixing process was simply treated as a passive scalar transport problem. However, for practical systems with H
$_2$
JICF, the density ratio is much smaller than 1 as H
$_2$
is much lighter than air. In addition, all the physical properties, such as density, viscosity and diffusivity, are functions of the local H
$_2$
concentration. Besides, the high diffusivity of H
$_2$
can lead to a very small molecular Schmidt number in the H
$_2$
-lean regions, which may further impact the turbulent mixing process. Nevertheless, to the best of our knowledge, currently there is no numerical study particularly targeting the mixing process within the H
$_2$
JICF configuration that is directly relevant to practical engine applications, despite the vital importance and the unique characteristics of H
$_2$
.
Hence the objectives of this paper are threefold. The first objective is to generate a high-resolution, fully-resolved DNS dataset for an H
$_2$
JICF configuration that is representative of key aspects of the intake system in a realistic H
$_2$
-fuelled off-road PFI engine. The second is to use the DNS data to assess the turbulence and mixing models in RANS and LES approaches that are widely used in industry for engine CFD. The third is to systematically identify fundamental mechanisms where current RANS mixing models fail, and to provide quantitative evidence for future model development.
2. Methodology
2.1. Physical problem: an inclined H
$_2$
JICF
An inclined H
$_2$
JICF configuration is considered in the present study. As shown in figures 1 and 2, gaseous H
$_2$
is injected into air crossflow with rectangular cross-section from the side wall through a round jet inclined at 60
$^{\circ }$
with respect to the wall. The jet diameter is
$D = 19$
mm, and the dimensions of the crossflow channel are 12
$D$
(streamwise,
$x$
direction)
$\times$
2.6
$D$
(transverse,
$y$
direction)
$\times$
3.6
$D$
(spanwise,
$z$
direction). The centre of the H
$_2$
injection hole is located 2
$D$
away from the channel inlet, and 1
$D$
away from the bottom wall. The temperatures of the H
$_2$
jet and air crossflow inlets are both 325 K. The pressure is 4.3 bar throughout the domain. The bulk velocities of the crossflow and jet are
$U_{\textit{bulk}}=1.66$
m s−1 and
$U_{\textit{jet}}=7.06$
m s−1, respectively. Consequently, Reynolds numbers for the jet and crossflow inlets (based on the hydraulic diameter) are 4600 and 22 100, respectively. More details about the layout of the H
$_2$
jet tube can be found in figure 2. Note that the origin of coordinates is placed at a location that is below the centre of the jet exit, as indicated in figure 2. This configuration is designed to maintain the same relative geometry and flow characteristics as in the real intake port of a heavy-duty H
$_2$
ICE, while the jet and crossflow inlet velocities are reduced to allow for fully resolved DNS.
Schematic for the flow configuration, computation domain and boundary conditions. The distribution of the H
$_2$
mass fraction in a specific horizontal plane cutting through the centre of the H
$_2$
jet hole is also presented for reference.
Figure 1. Long description
A diagram of the flow configuration, computation domain, and boundary conditions for hydrogen mixing in a jet in crossflow system. The diagram shows a rectangular domain with air inflow on the left side and outflow on the right side. Hydrogen inflow is introduced through a circular jet located on the top surface of the domain. The domain dimensions are labeled as 12D in length, 2.6D in width, and 3.6D in height. The walls of the domain are labeled as nonslip walls. The diagram also includes a color scale representing the mass fraction of hydrogen (Y_H2) ranging from 0 to 1. The distribution of the hydrogen mass fraction in a specific horizontal plane cutting through the center of the hydrogen jet hole is presented for reference.
Three views of the H
$_2$
jet tube and air channel.
2.2. The DNS configuration
The three-dimensional, multi-species, low-Mach-number Navier–Stokes equations that are solved in the DNS are given as follows:
\begin{align} \frac {{\rm D} \rho ^*}{{\rm D} t^*} + \rho ^*\, \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\textit {u}}^*=0, \\[-28pt] \nonumber \end{align}
\begin{align} \rho ^* \left ( \frac {\partial \boldsymbol{\textit {u}}^*}{\partial t^*} + \boldsymbol{\textit {u}}^* \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\textit {u}}^* \right )=-\boldsymbol{\nabla }p^{\prime *} + \frac {1}{Re_0}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }^*, \\[-28pt] \nonumber \end{align}
\begin{align} \rho ^* c_p^* \left ( \frac {\partial T^*}{\partial t^*} + \boldsymbol{\textit {u}}^* \boldsymbol{\cdot }\boldsymbol{\nabla }T^* \right )=\frac {\gamma _0 -1}{\gamma _0} \frac {{\rm D}p^*}{{\rm D}t^*}-\frac {1}{Re_0\, Pr_0} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q}^*- \left ( \sum _{n_s} \rho ^* c_{p,k}^*Y_k \boldsymbol{V}_k^* \right ) \boldsymbol{\cdot }\boldsymbol{\nabla }T^*, \\[-28pt] \nonumber \end{align}
\begin{align} \rho ^* \left ( \frac {\partial Y_k}{\partial t^*} + \boldsymbol{\textit {u}}^* \boldsymbol{\cdot }\boldsymbol{\nabla }Y_k \right )= -\frac {1}{Re_0\, Sc_0} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J}_k^*. \\[0pt] \nonumber \end{align}
Here, the superscript
$^*$
indicates that a variable is non-dimensionalised by its reference value. Specifically,
$\rho$
is the mixture density,
$\boldsymbol{u}$
is the fluid velocity,
$p$
and
$p'$
are thermodynamic and hydrodynamic pressures, respectively,
$T$
is temperature,
$Y_k$
is mass fraction of species
$k$
,
$c_p$
and
$c_{p,k}$
are the constant pressure specific heat of the mixture and the species
$k$
, respectively,
$n_s$
is the total number of species,
$Re_0$
,
$ \textit{Pr}_0$
and
$ \textit{Sc}_0$
are the Reynolds, Prandtl and Schmidt numbers, respectively, and
$\gamma _0$
is the specific heat ratio, defined on the reference quantities. The viscous stress tensor
$\boldsymbol{\tau }^*$
, the heat flux
$\boldsymbol{q}^*$
, and the diffusion flux
$\boldsymbol{J}_k^*$
are expressed as
\begin{align} \boldsymbol{\tau }^* = \mu ^* \left ( \boldsymbol{\nabla }\boldsymbol{u}^* + \boldsymbol{\nabla }\boldsymbol{u}^{*{\rm T}}-2/3\,\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^* \mathrm{\boldsymbol{I}} \right ), \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{q}^* = -\lambda ^*\, \boldsymbol{\nabla }T^*, \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{J}_k^* = \rho ^*Y_k \boldsymbol{V}_k^*, \\[-6pt] \nonumber \end{align}
where
$\mu$
is the dynamic viscosity,
${\boldsymbol{I}}$
is the identity matrix, and
$\lambda$
is the thermal conductivity. The diffusion velocity
$\boldsymbol{V}_k^*$
is given by
\begin{align} Y_k \boldsymbol{V}_k^* = Y_k \big ( \boldsymbol{V}_k^{*0} + \boldsymbol{V}_c^{*} \big ), \end{align}
where the mixture-averaged diffusion velocity
$\boldsymbol{V}_k^{*0}$
can be expressed as
\begin{align} Y_k \boldsymbol{V}_k^{*0} = -D_k^* \left ( \boldsymbol{\nabla }Y_k + Y_k\, \boldsymbol{\nabla }\ln \bar M \right ), \end{align}
and
$\boldsymbol{V}_c^*=-\sum _{n_s} Y_k \boldsymbol{V}_k^{*0}$
is the correction velocity to ensure
$\sum _{n_s} Y_k \boldsymbol{V}_k^* = 0$
. Here,
$D_k$
is the mixture-averaged mass diffusivity of the species
$k$
, and
$\bar M$
is the mean molecular weight of the mixture.
The thermodynamic properties of species, e.g.
$c_{p,k}$
, are computed using the NASA polynomials (Gordon & McBride Reference Gordon and McBride1976). The transport properties of species – i.e. the dynamic viscosity
$\mu _k$
, the thermal conductivity
$\lambda _k$
, and the binary diffusion coefficients
$D_{i,j}$
– are calculated from the Chapman–Enskog solution of kinetic gas theory (Poling, Prausnitz & O’Connell Reference Poling, Prausnitz and O’Connell2001; Kee, Coltrin & Glarborg Reference Kee, Coltrin and Glarborg2005).
The corresponding properties of gas mixture are then given by
\begin{align} \bar M = \sum _k X_k M_k, \\[-28pt] \nonumber \end{align}
\begin{align} c_p=\sum _k Y_k c_{p,k}, \\[-28pt] \nonumber \end{align}
\begin{align} \mu =\sum _i \frac {\mu _i X_i}{\sum _k \varPhi _{i,k} X_k}, \quad \varPhi _{i,k}=\left (8\left (1+\frac {M_i}{M_k}\right )\right )^{-{1}/{2}} \left (1+\left ( \frac {M_j}{M_k}\right )^{{1}/{4}}\left ( \frac {\mu _k}{\mu _j}\right )^{{1}/{4}} \right )^2, \\[-28pt] \nonumber \end{align}
\begin{align} \lambda = \frac {1}{2}\left (\sum _i X_i \lambda _i + \frac {1}{\sum _i {X_i}/{\lambda _i}} \right ), \\[-28pt] \nonumber \end{align}
\begin{align} D_{k}=\frac {1-X_k}{\sum _{i \neq k}X_i/D_{i,k}}, \\[-12pt] \nonumber \end{align}
where
$X_k$
and
$M_k$
are the mole fraction and molecular weight of species
$k$
, respectively.
To close this system, the equation of state for an ideal gas is employed to obtain the density directly:
\begin{align} \rho = \frac {p \bar M}{\mathcal{R} T}, \end{align}
where
$\mathcal{R}$
is the universal gas constant. With the low Mach number assumption (Tomboulides, Lee & Orszag Reference Tomboulides, Lee and Orszag1997), the hydrodynamic pressure
$p'=p'(\boldsymbol{x},t)$
is decoupled from the thermodynamic pressure
$p$
, which is set to be constant in this study.
The DNS of the H
$_2$
JICF are conducted using the spectral element CFD code Nek5000 (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2019), developed at Argonne National Laboratory. Nek5000 is high-order by design, and highly scalable on modern computing machines. It has been applied to a wide range of fluid dynamic problems (Schmitt et al. Reference Schmitt, Frouzakis, Wright, Tomboulides and Boulouchos2016; Giannakopoulos et al. Reference Giannakopoulos, Frouzakis, Boulouchos, Fischer and Tomboulides2017), including DNS of various wall-bounded flows (El Khoury et al. Reference Khoury, George, Schlatter, Noorani, Fischer, Brethouwer and Johansson2013; Hosseini et al. Reference Hosseini, Vinuesa, Schlatter, Hanifi and Henningson2016), wall-resolved LES of film cooling flows (Nunno et al. Reference Nunno, Wu, Ameen, Pal, Kundu, Abouhussein, Peet, Joly and Cocks2022), and DNS of turbulent reacting flows (Xu et al. Reference Xu, Ameen, Pal and Som2022; Wang et al. Reference Wang, Xu and Scarcelli2025a
,
Reference Wang, Xu and Scarcellib
).
The computational domain is discretised by a non-uniform mesh consisting of 268 000 spectral elements. To ensure that all turbulent length scales are well resolved, seventh-order polynomials are used in each element (namely, eight grid points along each spatial direction). Local refinement is applied to all walls so that the non-dimensional wall-normal distance satisfies
$y^+\lt 1$
at all grid points adjacent to the wall. Consequently, there are approximately 93 million grid points in total. A grid convergence study is provided in the supplementary material. Time integration is based on the semi-implicit second-order backward differentiation method with an operator-integrator factor scheme (Patel et al. Reference Patel, Fischer, Min and Tomboulides2017), allowing variable time-stepping with CFL number up to 2. Constant-velocity inflow that is normal to the surface is applied to the H
$_2$
jet and air channel inlet boundaries with prescribed values, while a no-slip boundary condition with fixed temperature 325 K is applied to all walls. Note that the turbulent inflow boundary condition has also been considered, and the results are consistent with those obtained using the constant-velocity inflow. Therefore, the latter is adopted to ensure identical inflow conditions across DNS, LES and RANS. Detailed comparisons of turbulent inflow effects are provided in the supplementary material. The time-averaging and Favre-averaging statistics are collected for 10 flow-through times (FTTs) after the flow has reached quasi-statistically steady state. The FTT is defined as the time required for the injected H
$_2$
to travel from the inlet to the outlet of the channel. For each FTT, it takes approximately 36 hours when using 4096 CPUs (AMD EPYC 7713 64-Core 2.0 GHz processor).
2.3. The LES and RANS configurations
The same problem is modelled in LES and RANS using CONVERGE CFD code version 4.0 (Richards, Senecal & Pomraning Reference Richards, Senecal and Pomraning2024), which is a finite-volume solver with the capability of simulating both incompressible and compressible flows in complex geometries. In this work, the compressible Navier–Stokes conservation equations are solved. A second-order central difference scheme is applied for the spatial discretisation of momentum, mass and energy transport equations, while a first-order upwind scheme is used for turbulence model related equations. For all governing equations, a first-order implicit Euler time integration is applied. For RANS simulation, the renormalisation group
$k{-}\epsilon$
model (Yakhot & Orszag Reference Yakhot and Orszag1986) is used, while for LES, the dynamic structure model (Pomraning & Rutland Reference Pomraning and Rutland2002) is chosen. Both turbulence models have been widely used in engine CFD simulations (Yosri et al. Reference Yosri, Ho, Meulemans, Talei, Gordon, Brear, Cosby and Lacey2021, Reference Yosri, Kar, Talei, Gordon, Brear and Khosravi2023; Wu, Torelli & Pei Reference Wu, Torelli and Pei2023; Wang et al. Reference Wang, Scarcelli, Bestel, Demir and Srna2024). Higher-order schemes, including second-order spatial discretisation for turbulence quantities (MUSCL scheme with MINMOD limiter for LES; second-order central scheme for RANS) and second-order time integration (Crank–Nicolson scheme) have also been tested, but showed only marginal influence on both LES and RANS results. Detailed comparisons are provided in the supplementary material.
Boundary conditions similar to those in DNS have been used for both LES and RANS. In particular, the same inlet velocity as in DNS but with a low turbulence intensity is applied to the two inlet boundaries. Regarding the wall treatments, the Werner and Wengle wall model (Werner & Wengle Reference Werner and Wengle1993) is employed for LES, while the scalable wall function (Richards et al. Reference Richards, Senecal and Pomraning2024) is used for RANS.
CONVERGE relies on a base Cartesian grid and a cut-cell approach to generate a computational grid at runtime, and features adaptive mesh refinement (AMR) and fixed embedding (FE) for local mesh refinement. In both LES and RANS, base grid size 4 mm is used, with 1 mm cells at walls. The AMR is applied to the entire domain to resolve local flow and mixing fields, with refinement criteria based on the local sub-grid velocity 0.02 m s−1 and sub-grid H
$_2$
mass fraction 0.003. Here, the sub-grid scale of a field is defined as the second spatial derivative of this field. For RANS, the refinement levels of FE and AMR are both 3 (denoted as ‘Level = 3’), corresponding to minimum cell size 0.5 mm, which is the typical mesh resolution used for engine simulations. For LES, a finer mesh is used, with the levels of FE and AMR both being 4 (denoted by ‘Level = 4’), leading to minimum cell size 0.25 mm. The resulting
$y^+$
is approximately 10 and 5 for RANS and LES, respectively.
For LES, the total cell count is 8.9 million, and each simulation takes 131 hours for each FTT when using 432 CPUs. For RANS, the total cell count is 1.2 millions and each simulations takes 16 hours for each FTT when using 216 CPUs. The turbulent statistics are collected for 5 FTTs after the flow has reached the statistically steady state.
2.4. Grid convergence in LES and RANS
To assess the mesh quality for LES, an index of quality (
$ \textit{IQ}_{\textit{LES}}$
) proposed by Celik, Cehreli & Yavuz (Reference Celik, Cehreli and Yavuz2005) is computed for two mesh resolutions (Level = 3 and Level = 4) as
\begin{align} \textit{IQ}_{\textit{LES}}=\frac {1}{1+\alpha _\nu ({\nu _{t,\textit{eff}}}/{\nu })^n}, \end{align}
where
$\alpha _\nu$
is constant 0.05, and
$n = 0.53$
. Here,
$\nu _{t,\textit{eff}}$
is the effective viscosity, which is the sum of the molecular viscosity (
$\nu$
) and the eddy viscosity (
$\nu _t$
) (Lien et al. Reference Lien, Li, Pati, Sadiki and Hasse2024). Figure 3 shows the instantaneous distributions of
$ \textit{IQ}_{\textit{LES}}$
at the horizontal plane
$z/D=1$
for these two meshes. As suggested by Celik et al. (Reference Celik, Cehreli and Yavuz2005),
$ \textit{IQ}_{\textit{LES}}$
values of 0.75–0.85 can be considered adequate for most engineering applications. In this sense, both meshes are sufficient for the current configuration since
$ \textit{IQ}_{\textit{LES}}$
is higher than 0.85 in most regions for both meshes. An alternative quality index based on the fraction of resolved turbulent kinetic energy also confirms that more than 95 % of the turbulent kinetic energy is resolved in most regions for both meshes, as shown in the supplementary material. Nevertheless, the Level = 4 mesh indeed shows higher quality than the Level = 3 mesh, especially at the leeward side of the H
$_2$
jet. This results in substantial differences in the distribution of time-averaged H
$_2$
mass fraction at the outlet vertical plane between the two meshes, as shown in figure 4. It is therefore suggested that the Level = 3 mesh is not sufficient to achieve an accurate prediction on the mixing process in LES, thus the Level = 4 mesh is used.
Instantaneous distributions of
$\textit{IQ}_{\textit{LES}}$
for two LES mesh resolutions at the horizontal plane
$z/D=1$
.
Unlike LES, the RANS prediction of the time-averaged H
$_2$
mass fraction at the outlet vertical plane is relatively insensitive to the mesh resolution, as shown in figure 5, suggesting grid convergence at Level = 3. Considering that Level = 3 mesh is much cheaper for the purpose of practical engine simulations, the Level = 3 mesh is used for RANS.
Time-averaged H
$_2$
mass fraction distributions at the outlet vertical plane
$x/D = 10$
for two LES mesh resolutions.
Time-averaged H
$_2$
mass fraction distributions at the outlet vertical plane
$x/D = 10$
for two RANS mesh resolutions.
3. Results and discussion
3.1. Comparison between DNS, LES and RANS
In this subsection, simulation results from DNS, LES and RANS are compared in terms of both flow and H
$_2$
mass fractions fields. The DNS serve as the benchmark for assessing the performance of LES and RANS models in predicting the flow dynamics and H
$_2$
mixing process.
3.1.1. Flow field predictions
Figure 6 shows instantaneous and time-averaged distributions of the velocity magnitude at the horizontal plane
$z/D = 1$
, predicted by DNS, LES and RANS. From the instantaneous field, DNS show strong turbulence generated by the shear layer between the H
$_2$
jet and air crossflow, resulting in a large amount of fine flow structure downstream of the jet. The LES, as expected, capture the large structure in the flow field. In contrast, RANS only predicts the main flow pattern near the jet exit without any fine flow structure being resolved. For the time-averaged flow field, the results from the three numerical approaches are similar in the region close to the jet exit. This is reasonable since all small structures have been filtered out by the averaging operation.
(a) Instantaneous and (b) time-averaged velocity magnitude distributions at the horizontal plane
$z/D = 1$
, predicted by DNS, LES and RANS.
Figure 6. Long description
Panel A: Instantaneous fields. The image shows three horizontal planes depicting instantaneous velocity magnitude distributions. The top plane is labeled DNS, the middle plane is labeled LES, and the bottom plane is labeled RANS. Each plane displays a color gradient representing velocity magnitude in meters per second, ranging from 0 to 14. The DNS plane shows the most detailed and turbulent flow patterns, followed by the LES plane with slightly smoother patterns, and the RANS plane with the smoothest and least detailed patterns. Panel B: Time-averaged fields. The image shows three horizontal planes depicting time-averaged velocity magnitude distributions. The top plane is labeled DNS, the middle plane is labeled LES, and the bottom plane is labeled RANS. Each plane displays a color gradient representing velocity magnitude in meters per second, ranging from 0 to 14. The DNS plane shows a more detailed and structured flow pattern compared to the LES and RANS planes. The LES plane exhibits a smoother flow pattern, while the RANS plane shows the smoothest and most simplified flow pattern.
Figures 7, 8 and 9 further show the time-averaged velocity components (
$\overline {u}$
,
$\overline {v}$
and
$\overline {w}$
) at different vertical planes from jet exit (
$x/D = 0$
) to the crossflow outlet (
$x/D = 10$
), predicted by DNS, LES and RANS. It is found that LES has a much better agreement with DNS than RANS. Specifically, in DNS and LES, the initial high
$\overline {u}$
jet quickly evolves into a ring-shaped cloud, then slowly transitions into a round-shaped diffused cloud. In the RANS results, however,
$\overline {u}$
decays along the streamwise direction at a much slower rate, resulting in a band-shaped high
$\overline {u}$
region near the bottom wall at downstream locations (
$x/D \geqslant 2$
). Besides, considerable differences in the distributions of
$\overline {v}$
and
$\overline {w}$
are also observed near the bottom left corner of the planes
$x/D = 4$
and 6. However, it is worth mentioning that despite the differences observed above, the velocity components predicted by RANS are generally in a good agreement with those from DNS, at least qualitatively, especially compared with the predictions for Reynolds stress components, which are shown in figures 10, 11 and 12.
Time-averaged velocity component
$\overline {u}/U_{\textit{jet}}$
at different vertical planes with
$x/D$
values 0, 2, 4, 6, 8 and 10, predicted by DNS, LES and RANS.
Time-averaged velocity component
$\overline {v}/U_{\textit{jet}}$
at different vertical planes with
$x/D$
values 0, 2, 4, 6, 8 and 10, predicted by DNS, LES and RANS.
Figure 8. Long description
Panel-wise Structure: Panel A: description, Panel B: description, Panel C: description.
Time-averaged velocity component
$\overline {w}/U_{\textit{jet}}$
at different vertical planes with
$x/D$
values 0, 2, 4, 6, 8 and 10, predicted by DNS, LES and RANS.
Figures 10, 11 and 12 show the time-averaged Reynolds stress components (
$\overline {u'u'}$
,
$\overline {v'v'}$
and
$\overline {u'v'}$
) at different vertical planes. Similar to previous results, LES still show comparable predictions to DNS regarding Reynolds stress components. However, dramatic differences are seen between the results from RANS and DNS. The same phenomena are also observed for the distributions of
$\overline {w'w'}$
,
$\overline {v'w'}$
and
$\overline {w'u'}$
, which are not plotted here for brevity. The interested reader may refer to the supplementary material for detailed results.
Time-averaged Reynolds stress component
$\overline {u'u'}/U_{\textit{jet}}^2$
at different vertical planes with
$x/D$
values 0, 2, 4, 6, 8 and 10, predicted by DNS, LES and RANS.
Figure 10. Long description
Panel-wise Structure: Panel A: description Panel B: description Panel C: description Each panel description must include: -Graph or visual type -Axis labels and units -Legends or markers -Data patterns or trends -Visible annotations.
Time-averaged Reynolds stress component
$\overline {v'v'}/U_{\textit{jet}}^2$
at different vertical planes with
$x/D$
values 0, 2, 4, 6, 8 and 10, predicted by DNS, LES and RANS.
Figure 11. Long description
Panel-wise Structure: Panel A: description, Panel B: description, Panel C: description.
Time-averaged Reynolds stress component
$\overline {u'v'}/U_{\textit{jet}}^2$
at different vertical planes with
$x/D$
values 0, 2, 4, 6, 8 and 10, predicted by DNS, LES and RANS.
To summarise, the flow field results indicate that LES perform well in predicting both the mean flow velocity and the Reynolds stress. In contrast, RANS does not predict the Reynolds stress well, while its prediction of the mean flow velocity is somewhat satisfactory.
3.1.2. The H
$_2$
distribution predictions
The predictions of the H
$_2$
mixing process by DNS, LES and RANS are further examined. Figure 13 shows the instantaneous and time-averaged distributions of the H
$_2$
mass fraction in the horizontal plane
$z/D = 1$
. Similar to figure 6, DNS show fine mixing structures in the instantaneous H
$_2$
mass fraction field, which is captured by LES to some extent. For the time-averaged field, it can be clearly seen that the LES prediction agrees better with DNS data, compared to RANS, especially at the periphery of the jet and downstream regions. In particular, RANS predicts a clear separation between the primary H
$_2$
jet and the secondary branch near the bottom wall, resulting in a low-H
$_2$
region on the leeward side of the jet. Such a low-H
$_2$
region is not observed in DNS or LES, suggesting that the mixing between these two H
$_2$
branches is under-predicted by the current RANS model.
(a) Instantaneous and (b) time-averaged H
$_2$
mass fraction distributions at the horizontal plane
$z/D = 1$
, predicted by DNS, LES and RANS.
Figure 13. Long description
Panel A: Instantaneous fields. The image shows three visualizations of instantaneous hydrogen mass fraction distributions. The top visualization, labeled DNS, displays a highly detailed and turbulent distribution. The middle visualization, labeled LES, shows a less detailed but still turbulent distribution. The bottom visualization, labeled RANS, presents a smoother and more simplified distribution. Panel B: Time-averaged fields. The image shows three visualizations of time-averaged hydrogen mass fraction distributions. The top visualization, labeled DNS, displays a smooth and gradual distribution. The middle visualization, labeled LES, shows a similar smooth distribution with slightly less detail. The bottom visualization, labeled RANS, presents the smoothest and most simplified distribution. The color scale at the bottom indicates the hydrogen mass fraction, ranging from 0 to 1.
Time-averaged H
$_2$
mass fraction distributions at different vertical planes with
$x/D$
values 0, 2, 4, 6, 8 and 10, predicted by DNS, LES and RANS.
To comprehensively evaluate the performance of LES and RANS models in predicting H
$_2$
mixing in JICF, figure 14 compares the time-averaged H
$_2$
mass fraction distributions at different vertical planes spanning
$x/D = 0{-}10$
. At
$x/D = 2$
near the jet exit, a kidney-shaped H
$_2$
distribution is observed in DNS, due to the presence of the counter-rotating vortex pair. However, since the H
$_2$
jet centre is closer to the bottom wall than the top wall, and there is an overall downward velocity due to the vertical section of the injection tube, the mixing process in the lower part is more suppressed compared to the upper part. Hence the lower half of the kidney-shaped profile is larger and features higher H
$_2$
concentration compared with the upper half. As the mixture moves downstream, H
$_2$
continues to mix with surrounding air. Consequently, the H
$_2$
cloud grows while the peak value of H
$_2$
mass fraction keeps decreasing. The LES results again show excellent agreement with DNS. A similar trend is observed in RANS as well; however, compared with DNS and LES, RANS predicts a higher concentration region at the bottom wall, and shows less mixing with the air in the top part of the channel, consistent with the velocity field.
In particular, at the outlet plane (
$x/D = 10$
, the rightmost column in figure 14), the kidney-shaped profile has almost completely vanished in DNS and LES results, with H
$_2$
occupying the majority of the lower half corner of the channel, resulting in a fair homogeneous H
$_2$
in that region. On the contrary, RANS under-predicts H
$_2$
mixing, as indicated by the smaller H
$_2$
cloud and higher peak value of H
$_2$
mass fraction. This is further quantified in figure 15, where the time-averaged H
$_2$
mass fraction profiles along the transverse (
$y$
) direction at different heights of the outlet vertical plane are plotted. Compared with DNS and LES, at
$z = 0H$
height, RANS predicts a much higher H
$_2$
mass fraction in the region near the side wall (from which the jet is injected into the crossflow), but much lower H
$_2$
mass fraction in the region far away from jet. At higher locations, such over-prediction from RANS shifted to the centre region (
$0.5\lt y/D\lt2$
). In contrast, the agreement between DNS and LES is excellent, with only marginal deviation near the wall.
Figure 16 further quantifies the maximum (MAX) and root mean square (RMS) of time-averaged H
$_2$
mass fraction at different vertical planes predicted by DNS, LES and RANS. It is seen that LES predictions are in very good agreement with DNS for both MAX and RMS along the streamwise direction. Although a similar trend along the streamwise direction is captured, RANS significantly over-predicts the MAX and RMS of H
$_2$
mass fraction at each downstream location, indicating that H
$_2$
is less mixed with air crossflow in RANS compared with LES and DNS.
Time-averaged H
$_2$
mass fraction profiles along the transverse (
$y$
) direction at different heights (i.e.
$z/H = 0, 1/8, 1/4, 3/8$
, where
$H= 3.6D$
is the dimension of the channel in the spanwise (
$z$
) direction) in the outlet vertical plane (
$x/D = 10$
).
Figure 15. Long description
The image contains four line graphs labeled (a), (b), (c), and (d), each depicting H2 mass fraction profiles along the transverse direction (y/D) at different heights (z) in the outlet vertical plane. Panel A: The line graph shows H2 mass fraction at z = 0 H. The x-axis represents y/D ranging from 0 to 2.5, and the y-axis represents H2 mass fraction ranging from 0 to 0.4. Three lines are plotted: a solid black line for DNS, a dashed red line for LES, and a dotted blue line for RANS. The DNS and LES lines show a similar trend with a peak around y/D = 1.5, while the RANS line shows a higher peak around y/D = 1. Panel B: The line graph shows H2 mass fraction at z = 2/8 H. The x-axis represents y/D ranging from 0 to 2.5, and the y-axis represents H2 mass fraction ranging from 0 to 0.15. The DNS and LES lines show a similar trend with a peak around y/D = 1.5, while the RANS line shows a higher peak around y/D = 1. Panel C: The line graph shows H2 mass fraction at z = 1/8 H. The x-axis represents y/D ranging from 0 to 2.5, and the y-axis represents H2 mass fraction ranging from 0 to 0.25. The DNS and LES lines show a similar trend with a peak around y/D = 1.5, while the RANS line shows a higher peak around y/D = 1. Panel D: The line graph shows H2 mass fraction at z = 3/8 H. The x-axis represents y/D ranging from 0 to 2.5, and the y-axis represents H2 mass fraction ranging from 0 to 0.12. The DNS and LES lines show a similar trend with a peak around y/D = 1.5, while the RANS line shows a higher peak around y/D = 1.
The (a) maximum and (b) RMS of the time-averaged H
$_2$
mass fraction at different vertical planes.
The above results clearly demonstrate that the current RANS model lacks the capability to satisfactorily predict the mixing process between the H
$_2$
jet and air crossflow. Considering the significant computational cost benefit of RANS compared to LES and DNS, it is of great practical relevance to gain deeper understanding of the root causes for the deficiencies in the RANS mixing model.
3.2. Assessment of turbulent species flux models in RANS
To capture turbulent mixing in RANS, the turbulent species flux model is required. In this subsection, the commonly used turbulent species flux model in RANS, i.e. the GDH model, will be thoroughly examined based on DNS data. Specifically, the ground-truth values of all turbulent transport properties that are involved in turbulent species flux model will be directly extracted based on the DNS data, and compared with those used in RANS. The goal is to identify the origins for the large deficiency of the GDH model observed in the previous subsection. Then the performance of the GDH model in characterising the anisotropy of turbulent diffusivity is evaluated to provide insights for further model improvement.
3.2.1. Turbulent transport properties
In RANS simulations for multi-component flows, the turbulent species flux term
$\widetilde {Y_k^{\prime \prime }u_j^{\prime \prime }}$
is usually modelled based on the GDH as
\begin{align} \widetilde {Y_k^{\prime \prime }u_i^{\prime \prime }} = - D_t \frac {\partial \widetilde {Y_k}}{\partial x_i}, \end{align}
where
$Y_k$
is the mass fraction of species
$k$
,
$u_i$
is the flow velocity in the
$i$
direction, and
$D_t$
is the turbulent diffusivity. Note that unlike (1.1), where time averaging is employed, Favre averaging is adopted here since there is substantial density variation in the current H
$_2$
JICF system. Therefore,
$\widetilde {f}$
denotes Favre averaging of quantity
$f$
, and
$f^{\prime \prime }=f-\widehat{f}$
. Typically, the turbulent diffusivity
$D_t$
is assumed to be correlated with the turbulent kinematic viscosity
$\nu _t$
through a turbulent Schmidt number
$ \textit{Sc}_t$
(see (1.2)). A default value
$ \textit{Sc}_t = 0.78$
is typically used in engine CFD applications.
The turbulent kinematic viscosity
$\nu _t$
is introduced for the closure of the Reynolds stress term in the momentum equation as follows:
\begin{align} -\widetilde {u_i^{\prime \prime }u_j^{\prime \prime }}=2\nu _t S_{ij}-\frac {2}{3}\delta _{ij}\left(k+\nu _t \frac {\partial \widetilde {u_i}}{\partial x_i}\right), \end{align}
where
$ S_{ij}=({1}/{2})( {\partial \widetilde {u_i}}/{\partial x_j}+{\partial \widetilde {u_j}}/{\partial x_i})$
is strain-rate tensor, and
$ k=({1/}{2})\,\widetilde {u_i^{\prime \prime }u_i^{\prime \prime }}$
is the turbulent kinetic energy. In RANS,
$\nu _t=C_\mu k^2/\varepsilon$
is adopted, where
$C_\mu$
is a model constant.
In (3.1) and (3.2), all terms except the turbulent transport properties (
$D_t$
and
$\nu _t$
) can be directly obtained from the DNS data by the time-averaging and Favre-averaging operations. As such, the representative values of
$D_t$
and
$\nu _t$
can be inferred using a ‘least squares sense’ optimisation method based on (3.1) and (3.2), as follows:
\begin{align} D_t = \frac {-\left(\dfrac {\partial \widetilde {Y_k}}{\partial x_j}\right)\widetilde {Y_k^{\prime \prime }u_j^{\prime \prime }}}{\dfrac {\partial \widetilde {Y_k}}{\partial x_j}\dfrac {\partial \widetilde {Y_k}}{\partial x_j}}, \\[-28pt] \nonumber \end{align}
\begin{align} v_t = \frac {\left(2S_{ij}-\dfrac {2}{3}\delta _{ij}\dfrac {\partial \widetilde {u_i}}{\partial x_i}\right)\left(-\widetilde {u_i^{\prime \prime }u_j^{\prime \prime }}+\dfrac {2}{3}\delta _{ij}k\right)}{\left(2S_{ij}-\dfrac {2}{3}\delta _{ij}\dfrac {\partial \widetilde {u_i}}{\partial x_i}\right)\left(2S_{ij}-\dfrac {2}{3}\delta _{ij}\dfrac {\partial \widetilde {u_i}}{\partial x_i}\right)}. \\[0pt] \nonumber \end{align}
Then the turbulent Schmidt number
$ \textit{Sc}_t$
can obtained as
\begin{align} Sc_t=\nu _t / D_t. \end{align}
Figure 17 compares the turbulent diffusivity
$D_t$
between DNS and RANS in vertical planes. Note that in the region with a small gradient of H
$_2$
mass fraction,
$D_t$
is ill-defined, thus regions with
$|\boldsymbol{\nabla }\widetilde {Y_k}|\lt 0.05$
are blanked. The DNS results show a high diffusivity near the jet core throughout the domain, and relatively a low diffusivity in the interface between H
$_2$
and surrounding air. For RANS, although it qualitatively predicts the distribution of
$D_t$
, the magnitude of
$D_t$
is significantly under-predicted throughout the domain. This clearly explains the under-prediction of the H
$_2$
mixing process by RANS as shown in previous sections.
Normalised turbulent diffusivity
$D_t/(U_{\textit{jet}}D)$
from DNS and RANS at different vertical planes. Regions with
$|\boldsymbol{\nabla }\widetilde {Y_{H_2}}|\lt 0.05$
are blanked.
Figure 17. Long description
Panel A: A heat map titled DNS shows normalized turbulent diffusivity at different vertical planes. The x-axis is labeled x/D, the y-axis is labeled y/D, and the z-axis is labeled z/D. The color scale ranges from 0 to 0.02, with darker colors indicating lower values and lighter colors indicating higher values. The heat map shows a series of cross-sectional slices at different x/D positions, revealing the distribution and intensity of turbulent diffusivity. Panel B: A heat map titled RANS shows normalized turbulent diffusivity at different vertical planes. The x-axis is labeled x/D, the y-axis is labeled y/D, and the z-axis is labeled z/D. The color scale ranges from 0 to 0.02, with darker colors indicating lower values and lighter colors indicating higher values. The heat map shows a series of cross-sectional slices at different x/D positions, revealing the distribution and intensity of turbulent diffusivity.
The value of
$D_t$
is determined by both
$\nu _t$
and
$ \textit{Sc}_t$
. In order to reveal the key contributing factors for the difference in
$D_t$
between DNS and RANS, figures 18 and 19 show the comparisons between DNS and RANS on
$\nu _t$
and
$ \textit{Sc}_t$
, respectively. Note that for
$ \textit{Sc}_t$
, only the results from DNS are plotted since
$ \textit{Sc}_t$
is a constant in RANS (
$ \textit{Sc}_{t,\textit{RANS}} = 0.78$
). Figure 18 shows a phenomenon similar to that observed in figure 17. That is, the magnitude of
$\nu _t$
is significantly under-predicted by RANS, although it predicts the shape of the overall distribution, i.e. high
$\nu _t$
near the jet core and low
$\nu _t$
along the interface. Such under-prediction is consistently observed when a second-order spatial scheme for turbulent quantities or a finer mesh is used, as shown in the supplementary material. Figure 19 plots the ratio of
$ \textit{Sc}_t$
between DNS and RANS (i.e.
$ \textit{Sc}_{t,\textit{DNS}}/Sc_{t,\textit{RANS}}$
). It is seen that
$ \textit{Sc}_{t,\textit{DNS}}/Sc_{t,\textit{RANS}}$
varies considerably throughout the domain, with most regions having
$ \textit{Sc}_{t,\textit{DNS}}/Sc_{t,\textit{RANS}}\lt1$
, and some small regions near the jet exit with
$ \textit{Sc}_{t,\textit{DNS}}/Sc_{t,\textit{RANS}}\gt2$
. Hence the under-prediction of
$D_t$
in RANS is largely attributed to the inaccuracies in the predictions of both
$\nu _t$
and
$ \textit{Sc}_t$
.
The above results indicate that the under-prediction of RANS on the H
$_2$
mixing process are attributed to both the under-prediction of turbulent viscosity
$\nu _t$
and the high value of
$ \textit{Sc}_t$
used in RANS. The former aspect has been rarely reported in the literature. However, the present study suggests that this should not be neglected in the context of model development for the purpose of improving models for H
$_2$
mixing. More detailed investigation on how to improve
$\nu _t$
prediction is out of the scope of this work, and merits further investigation in the future.
Normalised turbulent viscosity
$\nu _t/(U_{\textit{jet}}D)$
from DNS and RANS at different vertical planes. Regions with
$|\boldsymbol{\nabla }\widetilde {Y_{H_2}}|\lt 0.05$
are blanked.
Turbulent Schmidt number
$ \textit{Sc}_k$
from DNS at different vertical planes. Regions with
$|\boldsymbol{\nabla }\widetilde {Y_{H_2}}|\lt 0.05$
are blanked.
Figure 19. Long description
A heat map represents the distribution of turbulent Schmidt number from DNS at different vertical planes. The heat map is structured in a 3D grid layout with axes labeled as z over D, x over D, and y over D. The color scale ranges from 0 to 2, with blue indicating lower values and red indicating higher values. Regions with blanked values are not displayed. The heat map shows varying intensities of turbulent Schmidt number across different vertical planes, with some regions exhibiting higher values than others. The overall trend indicates a concentration of higher values in specific areas, while other regions show lower values.
To quantitatively illustrate the differences in
$D_t$
,
$\nu _t$
and
$ \textit{Sc}_t$
over the entire domain, figure 20 further shows the probability density functions (PDFs) of normalised
$D_t$
, normalised
$v_t$
, and
$ \textit{Sc}_k$
, which are collected from the three-dimensional domain. Note that the regions with
$|\boldsymbol{\nabla }\widetilde {Y_k}|\lt 0.05$
are excluded to rule out non-physical values. It is seen that the PDFs of
$D_t$
and
$\nu _t$
derived from DNS span a wider range towards higher values than those from RANS. In the meantime, the peak of the
$ \textit{Sc}_t$
PDF obtained from DNS is considerably smaller than the default value used in the current RANS model (0.44 versus 0.78). These values are consistent with observations in figures 17, 18 and 19. It is also worth mentioning that the most probable value of
$ \textit{Sc}_t$
observed in the present study is within the range reported in the literature (Ivanova et al. Reference Ivanova, Noll and Aigner2013; Ling et al. Reference Ling, Ryan, Bodart and Eaton2016; Ryan et al. Reference Ryan, Bodart, Folkersma, Elkins and Eaton2017), where the optimal value of
$ \textit{Sc}_t$
was derived based on LES data rather than DNS.
Comparison of the PDFs of normalised
$D_t$
, normalised
$v_t$
, and
$ \textit{Sc}_k$
between DNS and RANS. For
$ \textit{Sc}_k$
, only DNS results are plotted. The data are collected from the entire domain with
$|\boldsymbol{\nabla }\widetilde {Y_k}|\gt 0.05$
.
Turbulent species fluxes directly obtained from DNS and predicted by the GDH model with two different
$ \textit{Sc}_t$
values at the horizontal plane
$z/D = 1$
.
Figure 21. Long description
Panel A: A heat map showing flux in the x direction. The x-axis represents the horizontal distance, and the y-axis represents the vertical distance. The color scale ranges from -0.10 to 0.10, with blue indicating lower values and red indicating higher values. The DNS model shows a distinct pattern with higher values concentrated at the bottom left. The GDH default model shows a similar pattern but with less intensity. The GDH corrected model shows a more diffused pattern with some areas of higher intensity. Panel B: A heat map showing flux in the y direction. The x-axis represents the horizontal distance, and the y-axis represents the vertical distance. The color scale ranges from -0.10 to 0.10, with blue indicating lower values and red indicating higher values. The DNS model shows a distinct pattern with higher values concentrated at the bottom left. The GDH default model shows a similar pattern but with less intensity. The GDH corrected model shows a more diffused pattern with some areas of higher intensity. Panel C: A heat map showing flux in the z direction. The x-axis represents the horizontal distance, and the y-axis represents the vertical distance. The color scale ranges from -0.10 to 0.10, with blue indicating lower values and red indicating higher values. The DNS model shows a distinct pattern with higher values concentrated at the bottom left. The GDH default model shows a similar pattern but with less intensity. The GDH corrected model shows a more diffused pattern with some areas of higher intensity.
In figure 21, the species fluxes predicted from the GDH model with two different values of
$ \textit{Sc}_t$
– one is the default value 0.78, and the other is 0.44, which is the most portable value derived from DNS data (see figure 20
c) – are directly compared against DNS results. It is seen that the GDH model with default
$ \textit{Sc}_t$
has poor performance on predicting the species fluxes. Specifically, it under-predicts the positive species fluxes in most regions, and negative species fluxes in the
$y$
and
$z$
directions in the leeward of the jet, while it also over-predicts the negative species fluxes in the
$x$
direction in the windward of the jet. By using a corrected
$ \textit{Sc}_t$
, the performance of the GDH model is improved – however, only to some limited extent, as there are still some regions showing large discrepancies compared with DNS results. For example, the over-prediction for the negative species flux in the
$x$
direction is amplified by the corrected
$ \textit{Sc}_t$
. These results clearly demonstrate the importance of accounting for anisotropy for the turbulent mixing model in the context of the H
$_2$
jet mixing process, which will be the focus of the next subsubsection.
3.2.2. Anisotropy of
$ \textit{Sc}_t$
When adopting the GDH formulation to model turbulent species flux (i.e. (3.1) and (1.2)), it is implicitly assumed that the turbulent diffusivity
$D_t$
is an isotropic property. To evaluate this assumption, the anisotropic components of
$D_t$
in vector form, i.e.
$ \textit{Sc}_{t,x}$
,
$ \textit{Sc}_{t,y}$
and
$ \textit{Sc}_{t,z}$
, are separately derived from DNS data as
\begin{align} Sc_{t,x}=-\nu _t \frac {\partial \widetilde {Y_k}}{\partial x}/\widetilde {Y_k^{\prime \prime }u_x^{\prime \prime }}, \quad Sc_{t,y}=-\nu _t \frac {\partial \widetilde {Y_k}}{\partial y}/\widetilde {Y_k^{\prime \prime }u_y^{\prime \prime }}, \quad Sc_{t,z}=-\nu _t \frac {\partial \widetilde {Y_k}}{\partial z}/\widetilde {Y_k^{\prime \prime }u_z^{\prime \prime }}. \end{align}
Figure 22 plots the distributions of these three anisotropic components of
$ \textit{Sc}_t$
at different vertical planes. There are some regions with negative
$ \textit{Sc}_{t,i}$
, indicating the existence of counter-gradient diffusion, but they are limited to small regions and thus do not represent a significant contribution to the overall H
$_2$
fluxes. The more important observation is that
$ \textit{Sc}_{t,x}$
is much smaller than the other two components. In most regions,
$ \textit{Sc}_{t,x}$
is a small value, close to zero. By contrast,
$ \textit{Sc}_{t,y}$
varies dramatically, especially on the left-hand side of the jet (
$y/D\lt 1$
) where regions with
$ \textit{Sc}_{t,y}\gt 1$
are observed, while
$ \textit{Sc}_{t,z}$
shows a relatively homogeneous distribution at approximately 0.5 in regions with
$x/D\gt 4$
. These results clearly show the deficiency of isotropic
$D_t$
in the GDH.
Anisotropic components of
$ \textit{Sc}_t$
obtained from DNS at different vertical planes. Regions are blanked for
$|\boldsymbol{\nabla }\widetilde {Y_{H_2}}|\lt 0.05$
.
To quantitatively measure the anisotropy of turbulent diffusivity
$D_t$
, the misalignment angle
$\theta$
between the species flux vectors directly obtained from DNS (
$\boldsymbol{f}_{\textit{dns}}$
) and predicted by the GDH model (
$\boldsymbol{f}_{model}$
, e.g. the right-hand side of (3.1)) is calculated as
\begin{align} \theta = \arccos\left (\frac {|\boldsymbol{f}_{\textit{dns}}\boldsymbol{\cdot }\boldsymbol{f}_{model}|}{|\boldsymbol{f}_{\textit{dns}}|\,|\boldsymbol{f}_{model}|}\right). \end{align}
The analysis based on this quantity has the benefit of being independent of the choices of the model constants used since they only act as scaling factors and will not influence
$\theta$
. A smaller
$\theta$
would suggest better alignment between these two flux vectors, thus less anisotropy in
$D_t$
.
Figure 23(a) shows the distribution of
$\theta$
of the GDH at different vertical planes. It is seen that the GDH predicts large values of
$\theta$
(near 60
$^{\circ }$
) in most regions, especially near the jet exit, where
$\theta$
can be as high as 90
$^{\circ }$
. Figure 24 plots the PDF of
$\theta$
collected from the three-dimensional domain. It is seen that
$\theta$
of the GDH spans a wide range, from 0
$^{\circ }$
to 90
$^{\circ }$
, and the PDF peaks at approximately 37
$^{\circ }$
, which is considerably larger than 0. These results clearly demonstrate that the anisotropy in
$ \textit{Sc}_t$
is not negligible.
The distributions of misalignment angle
$\theta$
of different models at different vertical planes. Regions are blanked for
$|\boldsymbol{\nabla }\widetilde {Y_k}|\lt0.05$
.
Figure 23. Long description
Panel A: A 3D plot labeled GDH shows the distribution of misalignment angles at different vertical planes. The x-axis is labeled x/D, the y-axis is labeled y/D, and the z-axis is labeled z/D. The color scale ranges from blue to red, representing angles from 0 to 90 degrees. The plot shows distinct regions of varying angles, with higher angles indicated by red and lower angles by blue. Panel B: A 3D plot labeled GGDH displays the distribution of misalignment angles for a different model. The axes and color scale are the same as in Panel A. The distribution pattern is similar but shows slight variations in angle regions compared to Panel A. Panel C: A 3D plot labeled HOGGDH illustrates the distribution of misalignment angles for another model. The axes and color scale remain consistent with the previous panels. This plot also shows distinct regions of varying angles, with noticeable differences in the distribution pattern compared to Panels A and B.
The PDF of misalignment angle
$\theta$
collected from the three-dimensional domain with
$|\boldsymbol{\nabla }\widetilde {Y_k}|\gt 0.05$
.
As mentioned in the Introduction, the GGDH model proposed by Daly & Harlow (Reference Daly and Harlow1970) allows for anisotropy between components of the turbulent diffusivity
$D_t$
. This was achieved by including Reynolds stress components in the definition of
$D_t$
as
\begin{align} \widetilde {Y_k^{\prime \prime }u_i^{\prime \prime }} = - C_{\textit{GGDH}}\, \tau\, \widetilde {u_i^{\prime \prime }u_j^{\prime \prime }} \frac {\partial \widetilde {Y_k}}{\partial x_j}, \end{align}
where
$C_{\textit{GGDH}}$
is a model constant, and
$\tau$
is the turbulent time scale. Typically,
$\tau$
is taken to be
$k/\varepsilon$
, where
$k$
and
$\varepsilon$
are the turbulent kinetic energy and dissipation rate, respectively. In the present study, the definition
$\tau =\nu _t/k$
is adopted instead for the convenience of DNS data post-processing. Substituting this definition into (3.8) yields
\begin{align} \widetilde {Y_k^{\prime \prime }u_i^{\prime \prime }} = - C_{\textit{GGDH}}\, \frac {\nu _t}{k}\widetilde {u_i^{\prime \prime }u_j^{\prime \prime }} \frac {\partial \widetilde {Y_k}}{\partial x_j}. \end{align}
As such, the GGDH model can be regarded as replacing the constant scalar
$1/Sc_t$
used in the GDH by a matrix
$C_{\textit{GGDH}}\,\widetilde {u_i^{\prime \prime }u_j^{\prime \prime }}/{k}$
.
The HOGGDH model proposed by Abe & Suga (Reference Abe and Suga2001) can further improve the accuracy in predicting the streamwise component of the scalar fluxes. This was achieved by using quadratic products of the Reynolds stresses as
\begin{align} \widetilde {Y_k^{\prime \prime }u_i^{\prime \prime }} = - C_{\textit{HOGGDH}}\, \frac {\tau }{k} \widetilde {u_i^{\prime \prime }u_l^{\prime \prime }}\, \widetilde {u_l^{\prime \prime }u_j^{\prime \prime }} \frac {\partial \widetilde {Y_k}}{\partial x_j}, \end{align}
where
$C_{\textit{HOGGDH}}$
is a model constant. By adopting
$\tau =\nu _t/k$
, (3.10) can be written as
\begin{align} \widetilde {Y_k^{\prime \prime }u_i^{\prime \prime }} = - C_{\textit{HOGGDH}}\, \nu _t \frac {\widetilde {u_i^{\prime \prime }u_l^{\prime \prime }}}{k} \frac {\widetilde {u_l^{\prime \prime }u_j^{\prime \prime }}}{k} \frac {\partial \widetilde {Y_k}}{\partial x_j}. \end{align}
As such, the HOGGDH model can be regarded as replacing the constant scalar
$1/Sc_t$
in the GDH by a matrix
$C_{\textit{HOGGDH}}\,\widetilde {u_i^{\prime \prime }u_l^{\prime \prime }}\,\widetilde {u_l^{\prime \prime }u_j^{\prime \prime }}/{k^2}$
.
The same analysis of the misalignment angle
$\theta$
performed for the GDH is also applied to the GGDH and HOGGDH by employing (3.9) and (3.11) in (3.7), respectively. The results are presented in figures 23 and 24 for comparison with the GDH. It is observed that the GGDH significantly reduces the misalignment angle. Specifically, near the jet exit (
$x/D = 0,2$
),
$\theta$
decreases from approximately 90
$^{\circ }$
to approximately 60
$^{\circ }$
, although some local regions of high
$\theta$
are still present. In the downstream regions, the improvement is more pronounced, with a large portion of the plane exhibiting
$\theta$
of 15
$^{\circ }$
or less. The HOGGDH provides additional improvement over the GGDH, with
$\theta$
reduced across all planes. In the upstream regions,
$\theta$
of around 10
$^{\circ }$
is observed. Interestingly, the HOGGDH has worsened the alignment in a few local regions where
$\theta$
is always high in all three models.
As shown in figure 24, the PDF of
$\theta$
shifts towards smaller values for the GGDH and HOGGDH compared with the GDH, and the most probable value of
$\theta$
decreases to 18
$^{\circ }$
and 8
$^{\circ }$
, respectively. These results demonstrate the superiority of the GGDH and HOGGDH over the GDH in predicting the H
$_2$
JICF mixing process. In future work, a posteriori evaluations of these higher-order models would be of great interest.
4. Conclusions
In this work, a numerical study that employs DNS, LES and RANS approaches was conducted for an H
$_2$
jet in a crossflow based on a geometry representative of key aspects of PFI in an H
$_2$
-fuelled heavy-duty ICE. Based on the DNS data, the performance of LES and RANS on predicting the turbulent flow fields and mixing process in this specific configuration was comprehensively evaluated. It was found that LES has very good performance on predicting not only the mean velocity components but also the Reynolds stress components. In contrast, RANS significantly under-predicts all components of Reynolds stress, although it shows relatively good performance on predicting the mean flow field. Regarding the H
$_2$
mixing prediction, LES results also show excellent agreement with DNS results, while RANS still significantly under-predicts the mixing process.
The underlying reasons for the poor performance of RANS were systematically investigated by extracting turbulent transport properties used in the RANS approach from DNS data. Results showed that the turbulent diffusivity predicted by RANS is much smaller than that derived from DNS, which explains the under-prediction of the RANS mixing process. The under-prediction of turbulent diffusivity is found to be attributed to both the over-prediction of turbulent Schmidt number and the under-prediction on turbulent viscosity. By further analysing the anisotropic components of the turbulent Schmidt number and the misalignment angle between turbulent species fluxes directly obtained from DNS and those predicted by the mixing model used in RANS, the commonly used assumption of isotropic turbulent diffusivity was demonstrated to be invalid for the present configuration. In addition, higher-order species flux models that account for anisotropy were shown to exhibit superior performance over the default isotropic model.
These findings on anisotropic turbulent diffusivity and the relative contributions of turbulent viscosity versus turbulent Schmidt number errors have broader implications for turbulent mixing modelling in variable-density flows beyond the specific application studied here. The present study and the accompanying DNS dataset provide a foundation for developing improved models of turbulent viscosity and turbulent diffusivity in the RANS framework.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2025.11102.
Acknowledgements
The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (‘Argonne’). Argonne, a US Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The US Government retains for itself, and others acting on its behalf, a paid-up non-exclusive, irrevocable worldwide licence in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).
Funding
This work is performed under a Cooperative Research and Development Agreement (CRADA) between Argonne, Caterpillar Inc., and Convergent Science, Inc. The authors would like to thank G. Singh and N. Hansford, respectively Program Manager and Technology Manager at the Vehicle Technologies Office, Energy Efficiency and Renewable Energy, Off-road, Rail, Maritime, and Aviation (ORMA) Program, for supporting this work, and Caterpillar and Convergent Science for co-funding this work. The authors would also like to acknowledge the computing resources on the Bebop and Improv clusters provided by the Laboratory Computing Resource Center (LCRC) at Argonne National Laboratory.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data will be made available upon reasonable request.
2
2
2
IQLES
z/D=1
2
x/D=10
2
x/D=10
z/D=1
u¯/Ujet
x/D
v¯/Ujet
x/D
w¯/Ujet
x/D
u′u′¯/Ujet2
x/D
v′v′¯/Ujet2
x/D
u′v′¯/Ujet2
x/D
2
z/D=1
2
x/D
2
y
z/H=0,1/8,1/4,3/8
H=3.6D
z
x/D=10
2
Dt/(UjetD)
|∇YH2~|<0.05
νt/(UjetD)
|∇YH2~|<0.05
Sck
|∇YH2~|<0.05
Dt
vt
Sck
Sck
|∇Yk~|>0.05
Sct
z/D=1
Sct
|∇YH2~|<0.05
θ
|∇Yk~|<0.05
θ
|∇Yk~|>0.05
