Direct numerical simulation of high-pressure mixing in turbulent jets
Nek Sharan Josette Bellan
physics.flu-dyn physics.comp-ph

Quantum Criticality of Semi-Dirac Fermions in 2+1 Dimensions
    Bruno Uchoa
    Received: date / Accepted: date
============================================================




 
Combustion in automotive and aerospace applications employing diesel,
gas turbine and liquid rocket engines is preceded by injection and
mixing of fuel and oxidizer at high pressures, often exceeding mixture
critical values. Experimental observations indicate that the jets
injected at supercritical pressures exhibit significantly different
dynamics than the jets at subcritical conditions, owing to the lack
of distinct liquid and gas phases in supercritical state. As a result,
the averaged flow quantities such as the potential core length, jet
spatial growth rate and velocity decay profiles differ in the two
conditions, resulting in different mixed-fluid distributions. In this
study, turbulent jet direct numerical simulations (DNS) are performed
to examine the variations in statistics between injection of Nitrogen
(N_2) in Nitrogen (N_2) at subcritical
(perfect-gas) and supercritical conditions. Isothermal round jets
at Reynolds number (Re_D), based on jet diameter (D) and jet
orifice velocity (U_0), of 5000 and Mach number of 0.6
are considered. For mixing analyses, a passive scalar transported
with the flow is examined.






    Injection Mixing Turbulence Supercritical conditions







§ INTRODUCTION


Fuel injection and turbulent mixing at supercritical pressures determines
ignition and combustion in numerous engineering applications. Flow
evolution under such conditions is characterized by strong non-linear
coupling between dynamics, transport coefficients, and thermodynamics.
A model that accounts for these non-linear effects in computation
of the thermodynamic state of the mixture, and the heat and mass fluxes
in a multi-component fluid-flow simulation was proposed by Masi et
al. <cit.>. The goal of the present study is to evaluate
the model in a turbulent free-jet configuration to simulate fuel injection
and mixing in high-pressure(p) combustion chambers of propulsion
systems. 

Validation of numerical simulations at high-p conditions, where
theoretical results are scarce, requires comparisons with experimental
data. However, most experimental studies (e.g., <cit.>)
inject liquid (fuel at subcritical conditions) or high-density fluid
(fuel at supercritical conditions) at high Reynolds numbers (Re_D≥20,000),
and provide qualitative visual information about flow dynamics and
mixing. A numerical simulation of these flows, for direct comparisons
with the experiments, would require several models. For example, models
to accommodate potential two-phase/density-jump regions and to account
for the subgrid-scale fluxes would be required, in addition to the
thermodynamic model. Moreover, appropriate inflow and boundary treatments
that accurately replicate the experimental conditions are also imperative.
Interactions between the models and numerical details make validation
of individual models infeasible in such complex flows at high-Re_D.
Moreover, the lack of quantitative turbulence/mixing statistics measurements
in high-p experiments makes assessment of a model even more challenging. 

To mitigate the model interactions, direct numerical simulations at
Re_D of 5,000 are considered in this study. Jet flow at perfect-gas
conditions, for which theoretical <cit.>
and experimental <cit.>
results exist, is first considered to validate the numerical setup
and to create a database for flow-behavior comparisons against multicomponent
flows at high pressures. To systematically initiate the study, mixing
behavior in the single-species flows is assessed by a virtual passive
scalar transported with the flow, modeling diffusion at unity Schmidt
number (Sc) justifiable under perfect-gas conditions.

Accurate turbulent free-jet flow computation requires a careful choice
of inflow/boundary conditions, domain size, and numerical discretization.
The near-field jet flow evolution is particularly sensitive to the
choice of inflow perturbations, and several studies <cit.>
have examined its influence on turbulence statistics. The jet flow
attains self-similarity after the development region, containing potential
core collapse and transition to turbulence; the axial distance to
attain a self-similar state depends on the inflow perturbations. Moreover,
theoretical results <cit.> suggest that the self-similar
state may also depend on the inflow, requiring the similarity variables
to be appropriately modified for analyses. Although the theoretical
and experimental jet flow studies focus on the flow statistics in
the self-similar region, in practical applications with small combustion
chambers, the jet near-field is equally important because this is
where the phenomena determining the flame occur, and this is where control
must be exercised to influence flow behavior.




§ GOVERNING EQUATIONS


For the single-species flow of interest here, the conservation equations solved in this study are:

    ∂ρ/∂ t+∂/∂ x_j[  ρ
    u_j]  =0,


    ∂/∂ t(  ρ u_i)  +∂/∂ x_j[  ρ u_iu_j+pδ_ij-σ_ij]  =0,


    ∂/∂ t(  ρ e_t)  +∂/∂ x_j[  (  ρ e_t+p)  u_j-u_iσ
    _ij+q_j]  =0,


    ∂/∂ t(  ρξ)  +∂/∂ x_j[  ρξ u_j+J_j]  =0,

where t denotes the time, x is a Cartesian
coordinate, subscripts i and j refer to the spatial coordinates, u_i
is the velocity, p is the pressure, δ_ij is the Kronecker delta, 
e_t=e+u_iu_i/2 is the total energy (i.e., internal
energy, e, plus kinetic energy), ξ∈[0,1] is a virtual passive scalar
transported with the flow, σ_ij is the Newtonian viscous stress tensor

    σ_ij=μ(  2S_ij-23S_kkδ_ij)  ,  S_ij=12(  ∂ u_i∂ x_j
    +∂ u_j∂ x_i)  ,

where μ is the viscosity, S_ij is the strain-rate tensor, and
q_j=-λ∂ T/∂ x_j and J_j=-𝒟∂ξ/∂ x_j 
are the j-direction heat flux and scalar diffusion flux, respectively. λ 
is the thermal conductivity and 𝒟=μ/Sc is the scalar diffusivity,
where Sc denotes the Schmidt number. The injected fluid is assigned a 
scalar value of 1, whereas the chamber fluid a value of 0.

Two jet-flow simulations at conditions summarized in Table <ref>
are performed to examine flow statistics differences between injection
at perfect-gas and supercritical conditions. Only the chamber pressure
p_∞ differs between the two cases.

For the near-atmospheric-p simulation (Case 1), the perfect gas
equation of state is applicable, given by

    p=ρ R_uT/m,

where R_u is the universal gas constant and m is the species
molar mass. The viscosity is modeled as a power law

    μ=μ_R(T/T_R)^n

with n=2/3 and the reference viscosity μ_R=ρ_0U_0D/Re_D,
where ρ_0 and U_0 are the jet-exit fluid density and
velocity, respectively, and the reference temperature T_R=293K.
The thermal conductivity λ=μ C_p/Pr, where Prandtl
number Pr=0.7, the ratio of specific heats γ=1.4,
and the isobaric heat capacity C_p=γ R_u/(γ-1)
is assumed.

For the high-p simulation (Case 2), the governing equations (<ref>)-(<ref>)
are closed using the Peng-Robinson (PR) equation of state (EOS)

    p=R_uT/(v_PR-b_mix)-a_mix/(v_PR^2+2b_mixv_PR-b_mix^2),

where the pressure, p, and temperature, T, are obtained as an
iterative solution from the density, ρ, and internal energy,
e, obtained from the conservation equations <cit.>.
The molar PR volume v_PR=v-v_s, where the
molar volume v=m/ρ. v_s denotes the volume shift
introduced to improve the accuracy of the PR EOS at high pressures
<cit.>. a_mix
and b_mix are obtained from the expressions detailed
in <cit.>.

The physical viscosity, μ_ph, and thermal conductivity,
λ_ph, are calculated using the Lucas method <cit.>
and the Stiel-Thodos method <cit.>,
respectively. The computational viscosity, μ, and thermal conductivity,
λ, are then obtained by scaling μ_ph and
λ_ph with factor ℱ=μ_R/μ_ph,0,
i.e., μ=ℱμ_ph and λ=ℱλ_ph,
to simulate the flow at a specified Reynolds number Re_D. The
inflow physical viscosity, μ_ph,0, is obtained from
the Lucas method using the pressure p_∞ and the average
temperature (T_inj+T_ch)/2, where
the subscripts “inj” and “ch” denote
the injection and chamber conditions, respectively.

To examine the robustness of the above EOS and transport coefficient
models at supercritical conditions, Figure <ref>
compares the density, isobaric heat capacity, and the transport coefficients
μ_ph and λ_ph obtained from the
models against the National Institute of Standards and Technology
(NIST) database <cit.> for N_2 at a
pressure of 50 bar and temperatures ranging from 100 K to 400 K. The
supercritical temperature (T_c) of Nitrogen is 126.2 K. The transport
coefficient models are accurate only at supercritical temperatures
and, thus, the comparison only spans values of T>T_c. As evident,
the models have good agreement with the NIST database, showing their
validity at high-p conditions encountered during Case 2 simulation.
Figure <ref> shows the compressibility factor (Z),
indicating deviation from perfect-gas behavior, of pure Nitrogen for
a temperature range at 50 bar pressure. The compressibility factor
at chamber conditions in Case 2 is 0.994.



















§ NUMERICAL DETAILS


The spatial derivatives are approximated using the sixth-order compact
finite-difference scheme and time integration uses the explicit
fourth-order Runge-Kutta method. The outflow boundary in axial direction
and all lateral boundaries have sponge zones<cit.>
with subsonic non-reflecting outflow Navier-Stokes characteristic
boundary conditions (NSCBC)<cit.> at the boundary
faces. Sponge zones at each outflow boundary have a width of 10%
of the domain length normal to the boundary face. The sponge strength
at each boundary decreases quadratically with distance normal to the
boundary. The performance of one-dimensional NSCBC<cit.>
as well as its three-dimensional extension<cit.>
by inclusion of transverse terms were also evaluated without the sponge
zones; they permit occasional spurious reflections into the domain,
therefore, the use of sponge zones was deemed necessary. To avoid
unphysical accumulation of energy at the highest wavenumber, resulting
from the non-dissipative spatial discretization, the conservative
variables are filtered every five time steps using an eighth-order
filter.

The computational domain extends to 42D_0 in the axial (x-)direction
and 20D_0 in the y- and z-direction including the sponge
zones, as shown in Figure <ref>. 320×288×288
grid points are used in the x× y× z direction, which
is twice the number of grid points in each direction used for DNS
by Boersma<cit.> at similar conditions as Case
1.

The axial grid resolution is chosen to resolve all spatial scales
overwhelmingly responsible for the dissipation. Following the approach
outlined in <cit.>, for a Re_D=5000 jet
simulation, the Kolmogorov length scale η_K=0.0041(x-x_0),
where x_0 denotes the virtual origin. A stretched grid designed
accordingly is used for present simulations. Grid stretching is accounted
for by solving the governing equations in generalized coordinates
<cit.>.

The velocity profile at the jet inflow plane is given by<cit.>

    u(r)=U_0/2(1-tanh[r-r_0/2θ_0]),

where the jet exit radius r_0=D/2 and the momentum thickness
θ_0=0.04r_0 is assumed. The Mach number Ma_0=U_0/c_∞=0.6
is specified, where c_∞ denotes the speed of sound at ambient
conditions. Random perturbations with maximum amplitudes of 0.004U_0
are superimposed on the inflow velocity profile to trigger jet flow
transition to turbulence. No perturbations are added to fields other
than velocity.



















§ JET FLOW RESULTS AND DISCUSSION


Table <ref> summarizes the conditions for the
numerical simulations considered in this study. In both cases, the
injected and chamber fluid temperatures and pressures are the same,
therefore, the jet injects into a chamber that is as dense as the
injected fluid. In other words, Case 1 and 2 represent jets at same
temperature and exit velocity, but in Case 1 a near-atmospheric pressure
jet is injected into similar chamber conditions and in Case 2 a 50
bar pressure jet is injected into similar chamber conditions. Physically,
if the former jet has a Re_D of 5000, the latter jet, based on
the density and viscosity at 50 bar pressure and same exit velocity
U_0 and diameter D as the velocity and length scale, will
have Re_D∼240,000. A DNS at such Re_D is infeasible,
therefore, in this study Re_D is fixed to 5000 for both cases;
the Re_D=5000 condition is enforced using a computational viscosity
μ that is calculated by scaling the physical viscosity by a factor
ℱ, as discussed in Section <ref>.

ℱ≈6.5 (μ_R=ρ_0U_0D/Re_D=1.136×10^-4Pa.s
and μ_ph,0=1.757×10^-5Pa.s)
for Case 1, near atmospheric pressure, and ℱ≈309.4
(μ_R=5.715×10^-3Pa.s and μ_ph,0=1.847×10^-5Pa.s)
for Case 2 at 50 bar pressure. The factor ℱ is higher
in Case 2 because of the higher density ρ_0 at 50 bar that
requires a higher μ_R for a given Reynolds number Re_D.
The physical viscosity μ_ph,0, on the other hand,
remains relatively unchanged with increase in pressure. Simulating
the two jets at same Reynolds number by using a computational viscosity,
as explained above, results in the high-p jet becoming unphysically
more viscous than the atmospheric-p jet. The following comparisons
between the two cases must account for this effect.

Figure <ref> shows the Mach number contours
at tU_0/D≈2500 for both cases. The potential core is comparatively
shorter in the high-p case, likely due to higher computational
viscosity μ, from a higher value of factor ℱ, and
resulting momentum diffusion.

Figure <ref> shows the decay of time-averaged
centerline velocity (U_c) and scalar concentration (ξ_c)
with axial distance. For a self-similar round jet with top-hat exit
velocity profile, the centerline velocity U_c(x) is
given by the empirical relation<cit.>

    U_c(x)/U_0=B/(x-x_0)/D,

where B is a constant and x_0 denotes the virtual origin.
As evident from Figure <ref>(b), downstream
of the potential core collapse, both the time-averaged centerline
velocity and scalar concentration decays as inverse of the axial distance,
where the rate of decay, given by B, is within experimentally observed
range of values <cit.>. The scalar concentration
begins to decay upstream of velocity, consistent with the observation
of Lubbers et al. <cit.>
for a passive scalar diffusing at unity Schmidt number. Despite the
difference in axial location where the velocity or scalar decay begins
between the near-atmospheric-p and high-p cases, the profiles
match asymptotically with axial distance. At the time of reporting,
the scalar statistics for the high-p case did not fully converge;
we expect the dashed red line in Figure <ref>
to asymptotically converge to the solid red line (like the velocity
field shown in blue).

Figure <ref> shows the centerline root-mean-square (r.m.s.)
axial velocity and scalar concentration fluctuation, denoted by u'_rms
and ξ'_rms, respectively, normalized by the time-averaged
centerline values. The centerline r.m.s. values are calculated from

    u'_rms=((u-U_c)^2)^1/2=(u^2-U_c^2)^1/2

and a similar expression for scalar concentration. The overbar, ∙,
denotes a time average at the centerline. The r.m.s. velocity fluctuation
profiles in Figure <ref>(a) compare favorably with the
profile of Crow & Champagne <cit.>
and, similarly, the scalar fluctuation profiles in Figure <ref>(b)
compare well with the distributions of r.m.s. scalar fluctuation in
jets from smooth contraction nozzle shown in Mi et al. <cit.>.
As observed in Figure <ref>, despite
the differences in r.m.s. fluctuations in near-jet regions between
the high-p and near-atmospheric-p cases, the profiles match
asymptotically with axial distance. The fluctuation magnitude signifies
turbulence intensity, which is negligible in the potential core, increases
sharply with collapse of the core, and asymptotes downstream to a
constant value as the flow becomes self-similar. The rise in fluctuation
profiles at shorter axial distance in case of high-pressure jet shows
a shorter potential core. Relatively lower peak fluctuation amplitudes
in high pressure results is likely a manifestation of the higher computational
diffusivity from higher value of μ.




§ CONCLUSIONS


Single-species turbulent jet simulations at different chamber conditions
are performed as part of an effort to understand fuel injection and
fuel-oxidizer mixing at high pressures. Two cases with near-atmospheric
and supercritical chamber pressures, respectively, are analyzed while
keeping inflow, boundary conditions and initial condition identical.
The equation of state and transport coefficient models are chosen
specific to the two conditions. The transport property models at supercritical
conditions are validated against the NIST database values. To avoid
subgrid-scale model errors and its, often difficult to predict, interactions
with thermodynamic calculations in high-Reynolds-number simulations,
DNS at Re_D of 5000 is performed. Profiles of centerline velocity
and scalar concentration mean and r.m.s. fluctuations show favorable
agreement with the experimental data at atmospheric-p conditions.
For comparisons in this study, the near-atmospheric-p and high-p
jets were simulated at the same Reynolds number by using a computational
viscosity adjusted accordingly. Matching the Reynolds number in such
a manner results in the high-p jet becoming more viscous than the
atmospheric-p jet, thus, obscuring mixing statistics' one-to-one
comparisons. An alternative is to match the inflow momentum instead
of the Reynolds number Re_D, keeping the viscosity physical.
This alternative will be a subject of future investigation.




§ ACKNOWLEDGEMENTS

This work was supported by the Army Research Office under the direction of Dr. Ralph Anthenien. The computational resources were provided by the NASA Advanced Supercomputing at Ames Research Center under the T^3 program directed by Dr. Michael Rogers.



