Euler–Lagrange numerical simulation of a kerosene droplet mist ignition in air using analytically reduced chemistry

Abstract

High altitude relight is a critical aspect of the aeronautical engine certification and may be addressed with the numerical simulation of two-phase ignition. However, such configurations are stiff and combined with local evaporation may lead to numerical issues. This paper provides several methods to perform two-phase ignition simulations using analytically reduced chemistry in the context of unstructured large Eddy simulation and Euler–Lagrange formalism. Firstly, an exponential formulation combined with a local and dynamic sub-cycling of the stiff chemistry is demonstrated to allow stable integration at the flow time step. Secondly, a particle-bursting method is applied to limit the impact of stiffness induced by the Lagrangian point-source approach in fine meshes. These methods are then applied in the simulation of ignition of a mono-disperse, multi-component kerosene spray in air. The use of the analytically reduced chemistry model enables us to describe in detail the chemical structure of the flame kernel during its formation. Moreover, local increase of fuel concentration is observed as the ignition proceeds which has a large influence on the combustion processes and the flame kernel development.

Keywords

Ignition euler-Lagrange kerosene evaporation analytically reduced chemistry point-source spark high-altitude

Introduction

Ignition at high altitude condition is a critical aspect of aeronautical engine certification. As evidenced by experimental measurements,¹ low pressure ( $P$ = 0.3 bar) and low temperature ( $T$ = 233 K) conditions (labeled LPLT in the following) have a large detrimental influence on ignition. The ignition mechanism in gas turbine can be decomposed in several steps starting from the formation of the sparking to the full ignition of the combustion chamber.² The first step is the generation of the flame kernel by the ignition system and its growth. This step has been largely studied with experimental measurements,^3,4 and direct numerical simulations (DNS).^5–8 To reproduce the behavior of the spark device, an energy deposit is commonly used in numerical studies.⁹ This deposited energy must generate a hot gas region long enough for the fuel evaporation and the chemical combustion run-away to occur.

On one hand, LPLT conditions impact the chemistry of combustion and more specifically the initiation reactions. Therefore, an accurate chemistry description is required to study these phenomena, which is achieved in this work using the analytically reduced chemistry (ARC) formalism.¹⁰ On the other hand, LPLT conditions also influence the two-phase processes such as fuel atomization and droplet evaporation. Hence, the dispersed phase must be accounted for in ignition simulations to fully understand the effect of LPLT.

Because of the strong evaporation-combustion coupling, ignition of two-phase mixtures is a very stiff process. Indeed, fuel evaporation at droplet positions initially induces very localized fuel spots which undergo endothermic pyrolysis and exothermic oxidation, themselves producing localized hot and cold spots. This finally creates very strong fluctuations of evaporation and chemical reaction rates that are difficult to handle numerically. Therefore, new numerical strategies are required for stable computations.

To this purpose, several numerical strategies are proposed in this study. First, chemistry is integrated with a semi-implicit exponential formulation¹¹ associated to a local and dynamic sub-cycling (LDSC) procedure.¹² Then a particle-bursting method (PBM) is used to spatially re-distribute the local fuel evaporation source terms inherent to the point source assumption of the Lagrangian formalism. These methods are first presented in the following sections and then applied to the simulation of a three-dimensional (3D) two-phase ignition.

Chemistry

Chemical kinetics

The kinetic mechanism developed by the CRECK modeling group.¹³ has been used as reference. This detailed mechanism allows an accurate description of carbonated species combustion. The Jet-A1 aeronautical fuel kerosene is emulated using a multi-component surrogate composed of three species: 60% n-dodecane (NC12) that stands for the paraffin behavior, 20% methyl-cyclohexane (MCYC6) representing the cyclic species and 20% xylene (XYL) for the aromatics (in volume).¹⁴

To be affordable in 3D simulations, the number of species and reactions, and the stiffness of the chemistry must be reduced.¹⁵ This is achieved here with the software ARCANE¹⁶ which automatically and successively applies reduction techniques. Firstly, unnecessary species and reactions are discarded using DRGEP method¹⁷ to obtain a skeletal mechanism. The ARC mechanism is finally obtained by applying the Quasi-Steady-State Approximation (QSS).¹⁰ Under this assumption, species with a very short characteristic timescale are considered to have zero net chemical source term, leading to algebraic expressions for their concentration. Therefore, there is no need to solve their conservation equation, which drastically reduces the numerical cost and the chemical stiffness.

The chemical reduction is performed based on the following test cases: (a) (1D) purely gaseous premixed flames at LPLT conditions, targeting the laminar flame speed, the adiabatic temperature and species concentrations of the main products (CO $_{2}$ , H $_{2}$ O and CO); (b) (0D) constant pressure reactors at low pressure to verify the auto-ignition time. The resulting ARC mechanism is composed of 30 transported species, 299 reactions, and 22 QSS species and is labeled S30R299QSS22 in the following. Validation profiles are given on Figure 1 showing the laminar flame speed versus equivalence ratio and auto-ignition time versus initial temperature. A very good agreement is obtained with a detailed CRECK mechanism. Note that this reduced chemistry has been also evaluated in gaseous forced ignition configurations.¹²

Figure 1.

Comparison of analytically reduced chemistry (ARC) scheme (S30R299QSS22) with the reference detailed scheme (CRECK) for laminar flame speed versus equivalence ratio in a 1D premixed flame at $P$ = 0.3 bar and T = 233 K (left) and auto-ignition time versus initial temperature in a 0D constant pressure reactor with stoichiometric mixture at P = 0.3 bar (right).

Chemistry integration

The standard numerical integration of the chemistry in explicit solvers consists of a first-order explicit time integration of the species production rates ${\dot{ω}}_{k}$ and writes:

c_{k}^{n + 1} = c_{k}^{n} + Δ t \times {\dot{ω}}_{k}^{n}

(1)

where

c_{k}

is the concentration of species

k

n

and

n + 1

are successive iterations and

Δ t

is the current time-step of the computation. From equation (1), it is clear that depending on the chemical source term, a too large time-step may lead to a concentration outside

[0, c_{k}^{m a x}]

, where

c_{k}^{m a x}

is the concentration corresponding to pure species

k

in the mixture, or trigger numerical oscillations of species concentrations. Typically, the chemical stiffness encountered with finite rate chemistry approaches requires time-steps of the order of

10^{- 9} - 10^{- 11} s

due to very fast pyrolysis processes and high intermediate radical species reactivity in spark ignition where high temperatures are reached.

To avoid the high numerical cost induced by the explicit integration of stiff ARC mechanisms, the exponential integration method¹¹ is used and summarized below. The species are assumed to be produced or consumed following a first-order differential equation:

\frac{d c_{k}}{d t} = {\dot{ω}}_{k} = A_{k} c_{k} + B_{k}

(2)

where the species production rates can be decomposed into

A_{k} c_{k}

and

B_{k}

, the contributions to destruction and creation, respectively. This assumption is strictly valid only if all reactions are elementary with unity stoichoimetric coefficients for the reactants, which is the case for 99.2% of the reactions in the reduced scheme. The analytical solution is known and can be written in discretized form as:

c_{k}^{n + 1} = c_{k}^{n} + \frac{B_{k}^{n}}{A_{k}^{n}} \exp (A_{k}^{n} Δ t) - \frac{B_{k}^{n}}{A_{k}^{n}}

(3)

A_{k} c_{k} = \sum_{j = 1}^{N_{r e a c}} {\dot{ω}}_{k, j} for all {\dot{ω}}_{k, j} < 0; B_{k} = \sum_{j = 1}^{N_{r e a c}} {\dot{ω}}_{k, j} for all {\dot{ω}}_{k, j} > 0

(4)

where

N_{r e a c}

is the total number of reactions in the ARC mechanism and

{\dot{ω}}_{k, j}

is the chemical rate of species

k

associated with reaction

j

and is computed with Arrhenius laws. The production rates are then re-derived to be explicitly injected in the conservation equations using equation (1).

{\dot{ω}}_{k, e x p o}^{n} = \frac{c_{k}^{n + 1} - c_{k}^{n}}{Δ t}

(5)

In that way, the exponential integration can be seen as a semi-implicit integration with a prediction of the concentration at the iteration

n + 1

. By construction the exponential integration is not mass conservative:

\sum_{k = 1}^{N_{s p e c}} {\dot{ω}}_{k, e x p o} \neq 0

. Thus a correction method¹¹ is added to ensure the mass conservation. One can note that

{\dot{ω}}_{k, e x p o}^{n}

approaches

{\dot{ω}}_{k}^{n}

when the time-step goes to zero in equation (3). Hence conservative issues are reduced using chemical sub-cycling as presented in the following.

In this work, a LDSC procedure is applied to the computation of the chemical source terms.¹² The current time-step of the computation is divided into smaller time-steps ( $Δ t_{c h e m - S C} = Δ t / N_{S u b - C y c}$ ) for the chemistry integration. The number of sub-cycles $N_{S u b - C y c}$ is determined locally and at each iteration based on chemical time scales approximated by the inverse of species production rates. For each species, the maximum allowed time-step corresponds to 1% variation of the species concentrations:

Δ t_{c h e m, k} = \frac{c_{k}}{{\dot{ω}}_{k}} \times 0.01

(6)

In this study, the number of sub-cycles observed in reactive zones is of the order of one hundred and can reach five hundreds for highly reactive unstable regions.

The exponential integration method associated with the LDSC have been evaluated on gaseous ignition configurations and have demonstrated their efficiency.¹² These methods enable to use large flow time-steps with ARC chemistries up to the acoustic Courant-Friedrichs–Lewy and Fourier limits required by the compressible solver, without compromising the accuracy. Therefore, they largely reduce the computational cost of 3D computations.

Two-phase flow

Discrete multi-component evaporation model

Evaporation is a key element of the two-phase ignition as it determines the amount and the composition of the flammable mixture. The use of a chemical scheme composed of several fuel species with their own decomposition pathways enables us to consider the preferential evaporation effects on two-phase ignition. Moreover, the three-components surrogate model for Jet-A1 is consistant with the use of a discrete multi-components (DMC) evaporation model. In this study, the DMC model is built on top of the Spalding model¹⁸ and is briefly recalled here.

The Spalding droplet evaporation rate writes:

\frac{d m_{p}}{d t} = - S h^{*} π d_{p} ρ D \ln (1 + B_{M})

(7)

where

m_{p}

is the mass of the particle and

d_{p}

is the particle diameter. The gaseous density

ρ

and the fuel species diffusivity

D

are estimated through the classical 1/3 rule between the droplet surface and the far field state.

S h^{*}

is a modified Sherwood number proposed by Abramzon and Sirignano¹⁹ to better estimate the mass flux considering the boundary layer around the droplet. Finally, the mass Spalding number

B_{M}

is computed for multi-component mixtures as:

B_{M} = \frac{\sum_{k = 1}^{N_{f u e l}} Y_{k, s} - \sum_{k = 1}^{N_{f u e l}} Y_{k, \infty}}{1 - \sum_{k = 1}^{N_{f u e l}} Y_{k, s}} = \frac{Y_{k, s} - Y_{k, \infty}}{ε_{k} - Y_{k, s}}

(8)

where

Y_{k, s}

and

Y_{k, \infty}

are the species mass fractions of each fuel component respectively at the droplet surface and the far field. Only the components present in the liquid phase are considered in equation (8). Once

B_{M}

is determined, the mass fraction of the fuel component

k

in the evaporation flux

ε_{k}

can be computed and

{\dot{m}}_{p, k}

, the individual component evaporation rates write:

{\dot{m}}_{p, k} = ε_{k} {\dot{m}}_{p}

(9)

The surface mass fractions

Y_{k, s}

required in equation (8) are determined with a liquid equilibrium in equation (10) where

P_{s a t, k}

, the vapor pressure, is obtained by the Clausius-Clapeyron relation.

X_{k, s} P = X_{l i q, k} P_{s a t, k}

(10)

with

P

the pressure at the droplet position and

X_{l i q, k}

the liquid mole fraction of the fuel species k in the droplet. Typically, in a Lagrangian approach, the

Y_{k, \infty}

is taken as the value in the cell containing the droplet.

The DMC evaporation model has been validated in a previous study.²⁰

Point source correction

The standard Lagrangian formalism uses the point source approximation, in which droplets are represented as material points transporting the droplet properties. In this approximation, all source terms produced by droplets are applied to the gas in the cell containing the particle. As a consequence, when using fine meshes the droplet evaporation creates localized peaks in the Eulerian fuel species fields at the droplet position. In the case of ignition, where droplets are placed in a hot environment, the fuel vapor produced in one cell directly pyrolyse before diffusing, resulting in pointwise combustion occurring at the droplet location. This then leads to point-to-point strong gradient variations that are numerically difficult to handle. This problem is well known and a number of regularization methods may be found in the literature.²¹ However, most methods involve surrounding cells of the droplet-containing cell, where the operation of finding neighboring cells can be very costly when computing on an unstructured grid.

In this work, an alternative is proposed that stays efficient on unstructured grids. The PBM is based on the R-parcel technique which consists in computing only numerical particles which represent $R_{p}$ physical particles, instead of computing directly the $R_{p}$ physical particles.²² This R-parcel concept was originally introduced to decrease the computational cost of Euler–Lagrange simulations by taking $R_{p} > 1$ . It is used here with $R_{p} < 1$ to represent one physical particle with a cloud of numerical particles, therefore smoothing out the associated source terms. In the PBM method, the parent physical droplets are burst in $N_{b u r s t}$ child droplets with $R_{p} = 1 / N_{b u r s t}$ . These child droplets are randomly dispersed in a control volume around the initial parent droplet position and form localized droplet clusters. The child numerical particles have the same mass, temperature, and radii than the physical one and thus follow the same evaporation model. However, the mass and heat fluxes transmitted to the gaseous phase are weighted by $R_{p}$ to ensure mass, momentum, and energy conservation. Since the child droplets are dispersed around the physical one, all source terms are automatically spread in the control volume as illustrated in Figure 2. In this study, the control volume is a sphere of radius $R_{c} = 2 Δ x$ where $Δ x$ is the characteristic grid size and $N_{b u r s t} = 100$ to ensure a homogeneous repartition. $R_{c}$ is chosen as the smallest value ensuring the numerical stability and limiting the artificial diffusion of the fuel evaporation source terms.

Figure 2.

Illustration of the particle-bursting method (PBM) resulting droplet cluster with the corresponding evaporation rate field.

3D configuration

Two-phase ignition in a spherical bomb is considered and sketched on Figure 3. The domain radius is $5$ cm, yet, only the center, where combustion processes occur, is refined with a characteristic grid size $Δ x = 100 μ$ m which ensures at least 10 points in the flame front. Thus, no combustion models are required. The resulting unstructured mesh is composed of 4.61 million elements.

Figure 3.

Spherical bomb ignition configuration.

The domain is initially filled at the center with a stoichiometric mixture of droplets in air ( $ϕ_{g, 0} = 0$ and $ϕ_{t o t, 0} = ϕ_{l, 0} = 1$ where $ϕ_{t o t, 0}$ , $ϕ_{g, 0}$ and $ϕ_{l, 0}$ are respectively the total, gaseous and liquid initial equivalence ratios) at rest and at high altitude conditions ( $P = 0.3 b a r, T = 233$ K). The liquid phase is homogeneously distributed with single-size droplets diameter of initially $d_{p} = 10 μ$ m. The resulting liquid volume fraction is $α_{l} = 2.029 \times 10^{- 5}$ . The droplet temperature is initially equal to the gaseous temperature.

Non-reflecting pressure outlet boundary conditions are applied with the Navier-Stokes Characteristic Boundary Conditions formalism²³ to maintain constant pressure in the domain and to evacuate the pressure wave created at ignition.

The ignition is triggered with an energy deposit at the center of the domain. The energy deposit model²⁴ corresponds to a source term $E_{s}$ applied to the energy equation, and writes:

E_{s} (x, y, z, t) = \frac{ϵ_{i}}{(4 π)^{2} σ_{t} σ_{s}^{3}} e^{- \frac{(t - t_{0})^{2}}{2 σ_{t}^{2}}} e^{- \frac{(x - x_{0})^{2} + (y - y_{0})^{2} + (z - z_{0})^{2}}{2 σ_{s}^{2}}}

(11)

with

ϵ_{i} = 200 m J

the total deposited energy chosen to reach a maximum temperature of

3000

K during the deposit. The coordinates (

x_{0}, y_{0}, z_{0}

) correspond to the center of the deposit volume, and

t_{0} = 50 μ

s is the time of maximum energy deposition. Finally,

σ_{s}

and

σ_{t}

are the characteristic size and duration of the deposit respectively, and are derived from the Gaussian distribution so as to reach at a distance

Δ_{s, t} / 2

an amplitude of a hundredth of the amplitude at the deposition center:

Δ_{s, t} = σ_{s, t} \times 4 \sqrt{\ln (10)}

(12)

For this study the deposition duration is

2 Δ_{t} = 100 μ

s, and the radius is

Δ_{s} = 2 cm

The computation is performed with the AVBP code (https://www.cerfacs.fr/avbp7x/) which solves the reactive and compressible Navier-Stokes flow equations with a central finite volume Lax-Wendroff scheme, second order both in time and space.²⁵ The dispersed phase is solved using the Lagrangian formalism. Source terms of mass, momentum and energy from the liquid to gaseous phase are distributed to the closest nodes in the Eulerian grid in a two-way coupling approach.

Results

Overview

The temporal evolutions of maximum temperature and mean heat release rate weighted by the cell volume are displayed on Figure 4. There is a first endothermic ignition phase at $t = 70 μ$ s quickly followed by a large increase of heat release rate at $t = 90 μ$ s which corresponds to the oxidation processes and leads to the flame front formation. Then, the heat release rate stabilizes at a much lower value corresponding to a flame propagation mode. The kernel maximum temperature also stabilizes around $T_{m a x} = 2000$ K after the end of the energy deposit, which indicates an ignition success.

Figure 4.

Temporal evolutions of maximum temperature and mean heat release rate over the whole domain.

In the following sections, the two-phase combustion mechanisms are detailed for the three ignition phases: (a) until $t = 90 μ$ s, early ignition, (b) until $t = 120 μ$ s, flame front formation, and (c) until $t = 600 μ$ s, flame front propagation.

Early ignition phase

The energy deposit heats the gas, which then heats the droplets through conduction. Evaporation starts when the droplet temperature is high enough, and is, therefore, strongest at the deposit center where the temperature is highest, as shown on Figure 5. Due to the small droplet characteristic size, the droplet density is high and the resulting evaporation rate field projected on the Eulerian grid is quite smooth, and only a few empty spots appear. One can note that the PBM approach (similarly, to other point source regularization methods) induces an artificial diffusion of the evaporation source terms which could influence the combustion processes. However, for this case (mono-disperse droplet of size $d_{p} = 10 μ$ m), the resulting mean inter-droplet distance is $I = (V^{1 / 3} / (N_{p t c l}^{1 / 3} - 1) = 3.05 \times 10^{- 4} m$ without PBM. This value is small enough to generate a relatively homogeneous gaseous phase which will first lead to a mixed combustion regime and then to an external group combustion mode.²⁶ Therefore the PBM approach is not expected to drastically influence the combustion regime. The gaseous equivalence ratio (computed from the atomic balance of carbon, hydrogen and oxygen) increases accordingly and reaches stoichiometry at the kernel center where combustion can then start. Note that the equivalence ratio is more sensitive to the liquid phase discrete distribution and does not appear as smooth as evaporation because of diffusion effects towards empty evaporation regions. These regions are highlighted with the iso-surfaces on Figure 5 (right) delimiting the no-evaporation zones in a thickness of $8 Δ x$ centered on the cut-plane. These regions are well correlated with the fluctuations of equivalence ratio.

Figure 5.

Evaporation rate (left) and equivalence ratio (right) cut-plane fields with temperature iso-contours of $T = 1000, 1600$ and $2200$ K, at $t = 80 μ$ s. The iso-surfaces on the right picture correspond to the no-evaporation zones in a thickness of $8 Δ x$ centered on the cut-plane. The visualization window height is 1 cm.

Combustion indeed starts at the kernel center as shown on Figure 6, in a mixed reactor combustion regime, meaning that the thin propagating flame front has not formed yet. This is confirmed by the pool of reactive radical species such as OH found at the center where $T > 2000$ K. Interestingly, this flame kernel is bordered by a spherical endothermic zone around $T = 1600$ K. This endothermic process corresponds to fuel pyrolysis, which occurs everywhere in the reaction zone but is not masked by exothermic combustion only in the lower temperature region. This is well evidenced by the field of nC $_{7}$ H $_{14}$ species which is the first pyrolysis product of the fuel component MCYC6: although produced everywhere in the reaction zone, nC $_{7}$ H $_{14}$ only appears at the border where it is not consumed by oxidation reactions.

Figure 6.

Heat release rate [ $J . m^{- 3} . s^{- 1}$ ] (left), OH mass fraction (center) and nC $_{7}$ H $_{14}$ mass fraction (right) cut-plane fields with temperature iso-contours of $T = 1000, 1600$ K and $2200$ K, at $t = 80 μ$ s. The visualization window height is 1 cm.

The fuel components mass fraction fields displayed on Figure 7 well illustrate preferential evaporation where MCYC6 being the most volatile species evaporates first as soon as $T$ reaches $1000$ K. In the hottest center zone where $T > 1600$ K, this species is immediately pyrolyzed after evaporation and thus disappears. The same mechanism applies to XYL but this less volatile species requires a higher temperature to evaporate and is found only above $\sim 1500$ K. Finally, a negligible quantity of NC12 is observed in the domain (not shown) because the droplet temperature is still too low at that time to enable its evaporation. Indeed, Figure 8 indicates that this species has not started to evaporate at the end of the early ignition phase ( $t = 90 μ$ s). This means that this phase is totally driven by MCYC6 and XYL, which is a major difference with gaseous ignition of the same fuel blend where NC12 would pyrolyse simultaneously with the other components.²⁷

Figure 7.

MCYC6 mass fraction (left) and XYL mass fraction (right) cut-plane fields with temperature iso-contours of $T = 1000, 1600$ K and $2200$ K, at $t = 80 μ$ s. The visualization window height is 1 cm.

Figure 8.

Temporal evolution of the fuel component evaporation rates integrated over the whole domain during the early ignition phase.

Combustion in this early ignition phase has not much progressed even at the kernel center, where the concentrations of the main combustion products CO, CO $_{2}$ and H $_{2}$ O are still low and the oxidizer O $_{2}$ has been barely consumed.

Flame front formation

After the ignition phase at $t = 100 μ$ s, a propagating flame front has formed as shown on Figure 9 with the heat release rate and the radical species OH and CH $_{2}$ O. The flame front is located around the iso-contour $T = 1600$ K which corresponds to the critical temperature for fuel pyrolysis.

Figure 9.

Heat release rate [ $J . m^{- 3} . s^{- 1}$ ] (left), OH mass fraction (center) and CH $_{2}$ O mass fraction (right) cut-plane fields with temperature iso-contours of $T = 1000, 1600$ K and $2200$ K, at $t = 100 μ$ s. The visualization window height is 1 cm.

Similarly to the previous phase, fuel component vapor is visible along the reaction zone in the low temperature side, but this time NC12 has started to evaporate and is present, in much less proportion, about ten time smaller than MCYC6 and XYL (see Figure 10). Hence again at this stage, the flame properties are mainly driven by these two fuel components.

Figure 10.

NC12 (left), MCYC6 (center) and XYL (right) mass fraction cut-plane fields with temperature iso-contours of $T = 1000, 1600$ K and $2200$ K, at $t = 100 μ$ s. The visualization window height is 1 cm.

As shown on Figure 11 (right), the flame front propagates in a mixture around stoichiometry. The two-phase combustion regime is of the kind of weakly evaporation-controlled where the gaseous mixture has reached flammability before reaching the flame.²⁸ Therefore, the flame front propagation is not controlled by the droplet evaporation. The weakly controlled evaporation regime is characterized, for laminar propagating flames, by an evaporation peak located ahead of the flame front in the fresh gases. The situation is different in sparked ignition, where the droplet evaporation is triggered by the energy deposition rather than the flame front. Therefore, the peak evaporation is not located in the fresh gases but remains at the kernel center, until all droplets have evaporated, inducing large values of the equivalence ratio there as shown in Figure11.

Figure 11.

Evaporation rate [ $k g . s^{- 1}$ ] (left) and equivalence ratio (right) cut-plane fields with temperature iso-contours of $T = 1000, 1600$ K and $2200$ K, at $t = 100 μ$ s. The visualization window height is 1 cm.

The evaporated fuel components that evaporate inside the volume enclosed by the flame pyrolyze immediately, leading to the formation of light carbonated species such as CH $_{3}$ and C $_{2}$ H $_{2}$ (see Figure 12). These carbonated species however do not oxidize due to the lack of oxygenated species. Therefore, the temperature is reduced at the kernel center due to both evaporation and pyrolysis. The formation of such rich, not flammable and cooler mixture at the center may be detrimental for the kernel survival and ultimately endanger the ignition process.

Figure 12.

CH $_{3}$ (left), C $_{2}$ H $_{2}$ (center) and O $_{2}$ (right) mass fraction cut-plane fields with temperature iso-contours of $T = 1000, 1600$ and $2200$ K, at $t = 100 μ$ s. The visualization window height is 1 cm.

The number of droplets was initially set to obtain a stoichiometric total equivalence ratio, computed from the mixing with cold air before energy deposit. Since the pressure is constant in the domain, the temperature increase resulting from the energy deposit induces a local decrease of the gas density and an expanding radial flow, which both reduce the oxidizer mass content. Therefore, the two-phase mixture becomes rich and when the droplets start to evaporate, there are fewer oxygen atoms in the gas phase than initially targeted. The resulting total equivalence ratio can be estimated by correcting the initial value with the gas expansion effect as in equation (13), which gives the maximum value which can be reached:

ϕ_{t o t, m a x}^{*} = ϕ_{g, 0} + ϕ_{l, 0} \times \frac{ρ_{0}}{ρ_{m i n}}

(13)

where

ϕ_{t o t, m a x}^{*}

is the rescaled total equivalence ratio before ignition,

ϕ_{g, 0}

and

ϕ_{l, 0}

are respectively the initial gaseous and liquid components of the equivalence ratio and

ρ_{m i n}

and

ρ_{0}

are respectively the minimum and initial gaseous densities. In the present case, the theoretical value

ϕ_{t o t, m a x}^{*} = 12.5

is never reached as it conserves the same minimum density during the entire evaporation. At the end of the evaporation, the value

ϕ_{t o t, m a x} = 8.5

is observed.

Finally, the main combustion products are presented on Figure 13. A large amount of CO is released whereas the CO $_{2}$ production is limited which is characteristic of rich combustion processes. Furthermore, a small depletion of H $_{2}$ O and CO $_{2}$ is observed at the kernel center where the equivalence ratio is the highest. The formation of combustion products and the complete consumption of oxidizer at the kernel center indicate that the combustion is complete at the kernel center. This is also confirmed on Figure 9 where there is no exothermic heat release rate or oxidized radicals at the kernel center.

Figure 13.

H $_{2}$ O (left), CO (center) and CO $_{2}$ (right) mass fraction cut-plane fields with temperature iso-contours of $T = 1000, 1600$ and $2200$ K, at $t = 100 μ$ s. The visualization window height is 1 cm.

Figure 14 shows the temporal evolution of the smallest droplet radius, which is located at the kernel center. At the flame formation phase ( $t = 100 μ$ s) the droplet radius has been reduced by only 15%, which means that evaporation may continue in the next phase. At this time, however, energy deposition stops and the required heat for evaporation will be now taken from the gas. This may lead to either ignition success or failure, depending on the energy balance of the hot gas.

Figure 14.

Temporal evolution of the smallest droplet radius over the whole domain.

Flame front propagation

Jumping to much later time ( $t = 500 μ$ s) a propagating flame is observed, meaning a successful ignition (see Figure 15). Because energy deposition has stopped, the hot gas temperature has decreased slightly below 2000 K but stabilizes there. As may be seen from Figure 14, droplets completed evaporation in the hot gas region at around $t = 175 μ$ s. Thus at $t = 500 μ$ s all endothermic processes have also been completed. However, the resulting very high equivalence ratio, up to 5, does not allow oxidation of the pyrolysis products in this hot zone.

Figure 15.

Heat release rate [ $J . m^{- 3} . s^{- 1}$ ] (left), evaporation rate [ $k g . s^{- 1}$ ] (center) and equivalence ratio (right) cut-plane fields with temperature iso-contours of $T = 1000, 1600$ and $2200$ K, at $t = 500 μ$ s. The visualization window height is 1 cm.

Figure 15 (center) shows that evaporation now only occurs in a spherical zone around the flame front on the cold air side. This means that evaporation is fast enough in the pre-heat zone to feed the propagating flame, which has become therefore purely gaseous. The evaporation zone lies much further from the hot gas than the flame front because evaporation starts at low droplet temperature ( $T_{p} \approx 450 - 500$ K) compared to the gas temperature at the flame front location ( $T = 1600$ K).

Finally, due to the high equivalence ratio resulting from the complete evaporation, the flame evolves in a very rich mixture at the flammability limit ( $ϕ \approx 2 - 3$ ). At this equivalence ratio, the flame is thick and propagates slowly. This reflects on the maximum heat release rate which is one order magnitude lower than it was in the previous phase (see Figure 9). This too rich fresh mixture is therefore detrimental for the kernel development and may lead ultimately to extinction. In practice, however, the turbulent flow may counteract this effect by promoting mixing with fresh air in the pre-heat zone and lowering the too high equivalence ratio.

Conclusion

Thanks to numerical stabilization methods for stiff problems, namely exponential chemistry integration, local and dynamic sub-cycling (LDSC) and particle bursting (PBM), the ignition simulation of a mono-disperse, multi-component kerosene droplet cloud was performed. Accurate chemistry was described with an ARC scheme. In particular, the interest of the PBM was demonstrated, allowing to reduce the evaporation rate source term stiffness associated to the point source approximation.

Results clearly highlighted preferential evaporation effects on ignition in two-phase mixtures, with a dominant role of the most volatile species (methyl-cyclohexane and xylene) in the early phases. Due to the high droplet number density, all two-phase fields appeared quite homogeneous and the flame was weakly perturbed by the dispersed phase. The ARC scheme allowed to describe the endothermic pyrolysis of the fuel components, which was also found to play a role in the ignition scenario. Another important observed behavior was the evolution of the equivalence ratio toward very rich values, due to the decrease of oxygen content associated to the density decrease when the gas heats up, while the fuel content in the droplets does not change. This led eventually to a very rich propagating flame, which may endanger the full ignition process. At the end of the simulated sequence, a two-phase flame is obtained which propagates in a purely gaseous mixture as droplets fully evaporate ahead of the reacting zone.

These results are preliminary and should be extended to evaluate and quantify the influence of the control volume size on the combustion processes and the ignition mechanism. Especially for poly-disperse cases where the artificial diffusion effect may have a large influence. Moreover, the PBM still requires further study to define suitable physical criteria to determine the control volume size and the number of child particles, which also depend on the parent droplet diameter and the mesh discretization respectively. Once validated against experimental measurements, the proposed numerical methodology will allow to consider other and more complex cases, including poly-disperse effects with a realistic spray distribution, turbulence and influence of high-altitude conditions.

Footnotes

Declaration of conflicting interests

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors received the following financial support for this article’s research, authorship, and/or publication: This work was funded by DGAC under the project CEFORA led by Safran Aircraft Engine and Safran Helicopter Engine using GENCI ressources (Grant 2023-A0132B10157).

ORCID iDs

Antoine Pestre

Bénédicte Cuenot

References

Martinos

Zarzalis

Harth

. Analysis of ignition processes at combustors for aero engines at high altitude conditions. In Proc ASME Turbo Expo 2020.

Lefebvre

. Gas turbine combustion. 2nd ed. Taylor & Francis, 1998.

Renou

. An experimental study of freely propagating premixed flames at various lewis numbers. Combust Sci Technol 2001; 162: 347–370.

Cardin

Renou

Cabot

et al. Experimental analysis of laser-induced spark ignition of lean turbulent premixed flames: New insight into ignition transition. Combust Flame 2013; 160: 1414–1427.

Poinsot

. Flame ignition in a premixed turbulent flow. annual research briefs. center for turbulent research, 1991.

Baum

Poinsot

. Effects of mean flow on premixed flame ignition. Combust Sci Technol 1995; 106: 19–39.

Neophytou

Mastorakos

. Simulations of laminar flame propagation in droplet mists. Combust Flame 2009; 156: 1627–1640.

Reddy

Abraham

. Influence of turbulence-kernel intereactions on flame development in lean methane/air mixtures under natural gas-fueled engine conditions. Fuel 2013; 103: 1090–1105.

Collin-Bastiani

Vermorel

Lacour

et al. DNS of spark ignition using Analytically Reduced Chemistry including plasma kinetics. Proc of the Comb Inst 2019; 37: 5057–5064.

10.

Law

. Systematic approach to obtain analytic solutions of quasi steady state species in reduced mechanisms. J Phys Chem A 2006; 110: 13202–13208.

11.

Blanchard

Cuenot

Cazeres

. Chemical modelling for methane oxy-combustion in Liquid Rocket Engines. Acta Astronaut 2022; 190: 98–111.

12.

Pestre

Cuenot

Riber

. Evaluation of numerical methods for explicit chemistry integration and application on DNS of turbulent kerosene ignition at high altitude conditions. In Proc 10th European Combustion Meeting.

13.

Ranzi

Frassoldati

Stagni

et al. Reduced kinetic schemes of complex reaction systems: Fossil and biomass-derived transportation fuels. Int J Chem Kinet 2014; 46: 512–542.

14.

Humer

Seiser

Seshadri

. Experimental investigation of combustion of jet fuels and surrogates in nonpremixed flows. J Propul Power 2011; 27: 847–855.

15.

Law

Yoo

et al. Dynamic stiffness removal for direct numerical simulations. Combust Flame 2009; 156: 1542–1551.

16.

Cazères

Pepiot-Desjardins

Riber

et al. A fully automatic procedure for the analytical reduction of chemical kinetics mechanisms for Computational Fluid Dynamics applications. FUEL 2021; 303: 121247.

17.

Pepiot-Desjardins

Pitsch

. An efficient error-propagation-based reduction method for large chemical kinetic mechanisms. Combust Flame 2008; 154: 67–81.

18.

Spalding

. The combustion of liquid fuels. Symp (International) Combust 1953; 4: 843–864.

19.

Abramzon

Sirignano

. Droplet vaporization model for spray combustion calculations. Int J Heat Mass Transf 1989; 32: 1605–1618.

20.

Shastry

Cazères

Rochette

et al. Numerical study of multicomponent spray flame propagation. Proc of the Comb Inst 2021; 38: 3201–3211.

21.

Poustis

Senoner

Zuzio

et al. Regularization of the Lagrangian point force approximation for deterministic discrete particle simulations. Int J Multiphas Flow 2009; 35: 376–388.

22.

Dukowicz

. A particle-fluid numerical model for liquid sprays. J Comput Phys 1980; 32: 229–253.

23.

Poinsot

Lelef

. Boundary conditions for direct simulations of compressible viscous flows. J Comput Phys 1992; 101: 104–109.

24.

Lacaze

Richardson

Poinsot

. Large eddy simulation of spark ignition in a turbulent methane jet. Combust Flame 2009; 159: 1993–2009.

25.

Lax

Wendroff

. Systems of conservation laws. Commun Pure Appl Math 1960; 13: 217–227.

26.

Nakamura

Akamatsu

Kurose

et al. Combustion mechanism of liquid fuel spray in a gaseous flame. Physics of Fluids 2005; 17: 123301.

27.

Stagni

Esclapez

Govindaraju

et al. The role of preferential evaporation on the ignition of multicomponent fuels in a homogeneous spray/air mixture. Proc of the Comb Inst 2017; 36: 2483–2491.

28.

Rochette

Riber

Cuenot

. Effect of non-zero relative velocity on the flame speed of two-phase laminar flames. Proc of the Comb Inst 2019; 37: 3393–3400.