1. Introduction
The utilization of liquid water injection in combustion applications has obtained considerable attention within the scientific community owing to its numerous practical advantages. The interaction between liquid water and the flame is based on three primary effects [
1,
2,
3]:
Cooling effect: This results from the extraction of heat from the hot gases during the evaporation process.
Dilution effect: The presence of water acts as a diluent for the fuel and oxidizer, leading to a reduction in the heat released per unit of mass in the system.
Chemical effect: Water exhibits high efficiency in third-body reactions, influencing the flame structure (this effect is not analyzed in the current study, as it is considered a high-order effect).
The application of water injection in internal combustion engines (automotive and aeronautical engines) has been implemented since the last century (e.g., Pratt and Whitney J57 engine and the 1962 Oldsmobile Jetfire). Historical motivations for such applications include temporary increases in thrust and power output due to the augmented mass processed by the engine. Additional benefits encompass reduced maximum temperatures—in this way, the liquid water acts as thermal protection for structures—and decreased NO
x emissions, as these are generally produced at elevated temperatures [
4]. Previous research has demonstrated the impact of water injection on emissions in diesel engines and gas turbines [
5,
6,
7] and also showed enhanced engine performance. Liquid water injection also holds significance in safety applications, as experimental and numerical studies have illustrated its ability to suppress or mitigate explosions and fires [
8,
9].
Nowadays, a primary focus of combustion research involves reducing CO
2 emissions. In this context, investigations into the combustion of carbon-free fuels, particularly hydrogen, have intensified. Hydrogen combustion offers several advantages, including high availability and inherent renewability [
10,
11], which can contribute to the development of a future green and circular economy [
12,
13]. Additionally, hydrogen has the advantage of a high energy density per unit mass (141 MJ/kg). Despite these advantages, challenges such as low density, high volatility, wide flammability range, and safety concerns arise, which were previously encountered in rocket propulsion, for instance, and are now being considered for distribution and storage. The high diffusivity of hydrogen and tendency for preferential diffusion introduce complexities, leading to intrinsic instabilities, particularly thermodiffusive instability. This instability occurs when the effective Lewis number (
) of the mixture falls below a critical value, as characterized by the dispersion relation proposed by Matalon et al. [
14]. Another criterion for thermodiffusive instability is the Markstein length, representing the effect of the stretch rate on the flame’s displacement speed; a negative value indicates susceptibility to thermodiffusive instability [
15].
The propensity of mixtures with Lewis numbers below a critical threshold to exhibit intrinsic flame instability has been extensively studied [
14,
15,
16,
17,
18,
19,
20,
21,
22,
23] and remains a focal point of community interest. Spherical flames, which are typical configurations in internal combustion engines, enable the analysis of statistically non-zero stretch rates on flame velocities and the consideration of flame propagation characteristics (e.g., Markstein length). The evolution of expanding spherical flames becomes more intricate in the presence of preferential diffusion, as the close correlation between diffusion and topology introduces side effects [
24], as the local focusing (or defocusing) of reactants can generate super-adiabatic temperature regions. Moreover, when spherical flames grow and reach a certain radius, unstable cellular structures begin to form without external forcing. However, the current study excludes consideration of this scenario due to the computational constraints of our setup. The primary emphasis of the current study is on examining the development of turbulent spherical flames in the presence of preferential diffusion and an evaporating liquid phase. The introduction of non-unity Lewis number conditions coupled with a non-zero mean curvature results in a distinct evolution compared to that of a statistically planar flame [
25,
26]. Moreover, the interaction between the flame kernel and turbulence amplifies the growth rate of the average flame radius, as evidenced in a prior investigation [
27]. A growth rate proportional to the radius, rather than a constant rate over time, is observed. While there are instances in the literature where analyses have explored expanding flames interacting with turbulence [
27,
28] and have considered the impact of non-unity Lewis numbers [
29,
30], the authors are not aware of any studies comprehensively examining the simultaneous effects of turbulence, preferential diffusion, and the evaporation of liquid water on the evolution and propagation of spherical flames across various equivalence ratios, Lewis numbers, and water addition conditions.
The subsequent two sections will provide a presentation of the mathematical background and numerical implementation of our study.
Section 4 will delve into the discussion of results and observations, followed by the presentation of conclusions in the final section.
2. Mathematical Background
The method employed in this study to characterize the evolution of flame kernels in the presence of monodispersed water droplets utilizes a hybrid Eulerian–Lagrangian approach with two-way coupling. This approach entails representing the liquid phase as individual droplets treated as point sources and tracing their motion and evolution within the domain (Lagrangian phase). At the same time, the gas phase is characterized through standard conservation equations, with conditions for different fields computed at each grid point (Eulerian phase). Two-way coupling implies that the interaction between the two phases is considered bidirectionally, which means terms related to the other phase are present in the equations for both phases. It is noteworthy that, unlike four-way coupling, which considers particle–particle interaction, two-way coupling neglects this phenomenon. This omission is justifiable due to the sufficiently high average interdistance between the droplets (
, where
is the average interdistance and
is the Kolmogorov length scale), which guarantees the low probability of droplet–droplet collision. Moreover, the Lagrangian description of the particles is possible because the initial diameter of the droplets (
) is significantly smaller than the grid resolution (
) and the Kolmogorov length scale (
). This allows us to avoid resolving the flow at the particle–gas interface [
31,
32,
33]. It has been demonstrated that the present methodology accurately captures the evaporation characteristics in multiphase simulations [
34] when the droplet dimension is below the Kolmogorov length scale. Furthermore, in a study by De Chaisemartin et al. [
35] that analyzed the combustion of dispersed droplets with flames using both a Eulerian–Lagrangian approach and a fully droplet-resolved approach, similar behavior using the two methods was observed. However, the computational cost of the Eulerian–Lagrangian approach was shown to be 10 times lower than that of the fully Eulerian method.
In the Eulerian–Lagrangian formulation, for individual droplets, their positions, velocities, temperatures, and diameters (
,
,
, and
, respectively) are calculated over time using the following expressions:
In the various equations presented, terms associated with the gas field are evident, such as the gas velocity
and the gas temperature
. The droplet temperature equation assumes infinite thermal conductivity for the liquid phase, resulting in a uniform temperature across the particle equal to
. Notably, this equation incorporates the latent heat of vaporization (
), the gaseous specific heat (
), and the thermal Spalding number (
), which is considered identical to the mass transfer Spalding number. The Spalding number is determined as
. The water concentration at the interface
is computed based on equilibrium conditions through the partial pressure derived from the Clausius–Clapeyron equation:
The atmospheric pressure (
) and boiling temperature at that pressure (
) serve as references in this context. The molar masses of the gas phase (
) and liquid phase (
) are denoted as such, with
representing the particle temperature. Equation (
1) incorporates three distinct relaxation times (
,
, and
) associated with the evolution of particle velocity, diameter, and temperature, respectively. The definitions of these relaxation times are as follows:
The relaxation times depend on the properties of both the gaseous and liquid phases. In Equation (
3),
and
represent the mass density and specific heat of the droplets, respectively. The dynamic viscosity of the gaseous phase is denoted as
, while
and
are the Schmidt and Prandtl numbers, respectively, for the gaseous phase. These values are considered identical in Equation (
3). Similarly, the Sherwood and Nusselt numbers are identical and are calculated as
. Finally,
is an aerodynamic resistance coefficient computed as
, where
is the Reynolds number of the droplet (
). The general conservation equation for the gaseous phase is:
The variable
is a general variable and, depending on the equation, assumes a value of (1,
,
e,
,
, or
), while
represents (1,
,
,
,
, or
) as appropriate. Here,
e represents the specific stagnation internal energy. The transport coefficients
in the equation are treated as constants with respect to temperature and composition, and their values are derived from the dynamic viscosity
. The term
is a source/sink term in the gas phase equation: for instance, representing the pressure gradient in the momentum equation. On the other hand, the term
serves as the source/sink term associated with the interaction with droplets and is calculated as follows:
This term is computed individually for each droplet, and its impact is subsequently distributed among the nearest grid points. Mass diffusion of distinct species is addressed through the constant Lewis number approach, with unity Lewis numbers considered for oxygen and water. For hydrogen, the Lewis number (
) is set at
, which corresponds to unburned hydrogen’s Lewis number at a very low equivalence ratio as computed with Cantera [
36]. To isolate the effects of preferential diffusion, additional simulations were conducted with an
artificial hydrogen possessing identical characteristics as real hydrogen but assigned a Lewis number of unity.
While three-dimensional direct numerical simulations (DNSs) with detailed chemistry are currently feasible [
37,
38], such simulations remain computationally demanding, particularly in the context of an extensive parametric study. In this work, an irreversible single-reaction chemical mechanism is therefore considered [
39]. The approach involves adjusting the Zeldovich number (
) and the heat of combustion (
H) based on the local equivalence ratio (
) [
39]. The expressions for the Zeldovich number and heat release are as follows [
39]:
The reaction rate is then calculated via an Arrhenius-type formulation:
Here and in the following,
T indicates the non-dimensional temperature defined as
. In the equation of the chemical mechanism, the heat release parameter, denoted as
, is defined by the expression
. With the selection of
K, the heat release parameter has a value of
. This single-step reaction mechanism has demonstrated its ability to capture the variation in laminar burning velocity with the gaseous equivalence ratio obtained from detailed chemistry simulations [
1]. Several results presented in this article are based on a progress variable denoted as
c. This variable increases monotonically from zero under unburned conditions (
u) and unity in the burned gas under equilibrium conditions (
b). In this context, the progress variable is specifically defined as the normalized oxygen mass fraction:
This definition has already been used in many previous works [
1,
2,
3,
40].
3. Numerical Implementation
Simulations were conducted utilizing a 3D compressible reactive flow solver: SENGA+ [
1,
2,
3,
40]. Spatial discretization employed a 10th-order central finite difference scheme, transitioning to a second-order one-sided finite difference scheme at the domain boundaries. Temporal advancement utilized an explicit third-order low-storage Runge–Kutta method [
41]. In this study, a series of spherical flames with varying equivalence ratios (
and
) were simulated, both with and without water addition (
and
, where
is the ratio of added water to the total mass,
, where
is the added mass of water and
is the sum of the fuel
and air
masses in the unburned part of the simulated volume). Two different types of water droplets, characterized by
and
, were used for water addition. The thickness of the unstretched stoichiometric flame, denoted as
, is defined as:
The subscript
L designates laminar unstretched conditions. The initial droplet temperature is set equal to the unburned gas temperature
for all the cases. Although the initial temperature of the droplets affects the heat needed to reach the water saturation temperature, its effect is much weaker than the cooling effect caused by water evaporation due to large magnitude of the latent heat of evaporation for water. The different parameters for the simulations are presented in
Table 1. The influence of a non-zero mean curvature, which is the case for spherical expanding flames, significantly affects the development and propagation of the flame, particularly under non-unity Lewis number conditions [
25,
26]. The method employed for initializing the flame kernel involves adopting a 1D laminar flame profile with a predetermined equivalence ratio and distributing it in a spherical shape. The kernel is positioned at the center of the domain, and conditions are set as fully burnt up to
. Beyond this point, the remaining 1D flame profile is mapped spherically. However, this procedure introduces discrepancies in thermochemical fields, especially in the presence of preferential diffusion, since the presence of a non-zero average curvature modifies the evolution of the local thermochemical fields due to the effects of diffusion. In the current dataset, the case with
and
is the one that evolves the fastest at the same stoichiometric chemical timescale (
, with
as the thermal diffusivity of the mixture). The stoichiometric
flame is allowed to evolve under laminar conditions until the thermochemical fields are relaxed, and this establishes a target flame radius of around
. In this way, all cases are restarted under laminar conditions. Turbulence is superimposed when the target flame radius is reached, and then the simulations are continued for 2
. The turbulent field, which follows the Batchelor–Townsend energy spectrum, has an integral length scale
and intensity
, which yields a flame in the thin-reaction zone regime. Water addition is initiated at the beginning when the conditions are still laminar throughout the entire domain except for the central part containing the ignited flame. The decision to introduce droplets at the start of the simulations rather than when turbulence is initiated does not significantly impact the subsequent evolution of the flame or the evaporation process. This is a consequence of the fact that the primary mass and heat transfer processes between phases generally occur in the hotter regions of the flame (i.e., the flame and burned gas region), and the droplets are initially positioned outside of these areas. Boundary conditions are established as standard non-reflective inflow–outflow (NSCBC), as there are no pronounced flow directions in this case, in contrast to previous configurations with statistically planar flames [
1,
2,
3]. The domain is cubic with dimensions of
, and the number of grid points is set at
. This ensures a minimum of 10 points to resolve the thermal flame thickness
and at least 2 points for the Kolmogorov length scale
. Similar configurations have been explored in several previous studies [
30,
42].
4. Results
Despite the simulations being conducted for a predetermined duration after the start of turbulence, regardless of the simulation parameters, the duration varies due to the differing time required for each kernel to reach the target radius. Consequently, the concentration level of steam differs slightly between the cases upon the activation of turbulence. Nonetheless, steam concentrations remain low across all cases at the time of turbulence activation because the droplets reside outside the burned gas region, and the differences in steam concentration are minimal for this reason.
The behavior of the progress variable and non-dimensional temperature field are equal in adiabatic and equidiffusive flames (i.e., ). However, our simulations encompass non-adiabatic conditions owing to droplet evaporation. Consequently, we expect a notable difference between the progress variable and the non-dimensional temperature, particularly in proximity to evaporating droplet locations due to their cooling and preferential diffusion effects.
Figure 1 shows iso-surfaces of the progress variable and non-dimensional temperature at
and
. The two iso-surfaces exhibit similar topological characteristics, albeit with differences that are mainly located in regions where droplet evaporation has occurred. The formation of small dimples is attributed to the cooling effect induced by the phase transition of water going from the liquid to the gaseous state. Notably, this cooling effect manifests more prominently in the temperature field, while the progress variable is only indirectly influenced by evaporation and the effects are less discernible. Moreover, the impact of a non-unity Lewis number is less relevant in comparison to non-adiabatic effects in this figure, especially because at stoichiometric conditions, preferential diffusion is less significant compared to leaner conditions, such as
. An essential consequence of preferential diffusion is the occurrence of regions where the temperature surpasses the adiabatic flame temperature corresponding to the specific equivalence ratio. The examination of the non-dimensional temperature field underscores this characteristic trait of lean flames with
.
Figure 2 illustrates the
mid-plane depicting the non-dimensional temperature
fields from the various simulated cases in this study at
after the initiation of turbulence. It is immediately clear that the maximum temperature is much higher in the cases with
and
than in their counterparts without preferential diffusion. In the equidiffusive cases,
is mostly around 1, while there is much more occurrence of
with preferential diffusion, as is also confirmed later on in
Figure 3. Furthermore, it is clear that the case with
and without preferential diffusion exhibits the slowest flame propagation. Moreover, the temperature field in the case with
and
demonstrates significant distortion induced by turbulence. This pronounced wrinkling of the flame surface is attributed to the relatively stronger impact of preferential diffusion under lean conditions, which leads to an increased flame surface area, which is consistent with previous findings [
20,
43,
44]. The importance of preferential diffusion as a function of the equivalence ratio can be evaluated through the effective Lewis number
. This parameter considers the overall behavior of the mixture concerning mass and temperature diffusion. The formulation suggested by Bechtold and Matalon for lean conditions is as follows [
43]:
This definition comprises a weighted mean of the Lewis number of the primary reactants, contingent upon the equivalence ratio via the function
A, which is influenced by the Zeldovich number
.
In
Table 2, the effective Lewis numbers are reported for cases A-F for the two equivalence ratios examined in this study, revealing a discernible decrease in
as
diminishes. Additionally, in
Figure 2, a cooling effect stemming from droplet evaporation is observed, which is dependent on the maximum gas temperature and residence time of the droplets. Consequently, the most pronounced cooling effects are observed in the case with
and
, where the maximum temperature is only slightly less than the case with
and the duration of the interaction of the droplets with the flame and the burned gasses is higher than that of the stoichiometric case with
. Conversely, despite a significantly longer droplet–flame interaction time in the case with
and no preferential diffusion, evaporation remains minimal due to the relatively lower burned gas temperature.
The probability density functions depicted in
Figure 3 substantiate in a quantitative manner the earlier observations based on slices of the temperature field in
Figure 2. Specifically, cases characterized by preferential diffusion under lean conditions manifest regions featuring super-adiabatic temperatures. Analysis of the probability density functions indicates that the cases with
generally yield flames exhibiting higher occurrences of super-adiabatic temperatures than the stoichiometric cases when
. This discrepancy arises due to the much lower effective Lewis number of the mixture, which causes more favorable conditions for this particular preferential diffusion effect. Moreover, the influence of turbulence-induced stretch enhances the consequences of preferential diffusion, which is particularly pronounced in lean scenarios compared to stoichiometric ones due to the favorable effects on topology, as observed by Concetti et al. [
45]. Notably, temperature uniformity is enhanced in instances where preferential diffusion is absent. Probability density functions in the right column attain higher values and lack the extended tails observed in the results in the left column due to the lower thermal non-uniformity in cases with
than in the cases with
. Across nearly all cases, temperature distributions demonstrate only moderate sensitivity to the addition of water, even for the case with the smallest initial diameter. However, the stoichiometric case without preferential diffusion exhibits the most substantial change of the probability density function, as evidenced by the peak shifting to lower temperature values and a reduction in the peak height attributable to thermal non-uniformity arising from water droplet evaporation. Nonetheless, all droplet diameters fail to induce markedly distinct behavior in the temperature distributions, even in instances where evaporation is the most intense among all cases considered in this study. Examination of the
fields offers a qualitative assessment of evaporation levels and their anticipated dilution effects on the flame.
Figure 4 illustrates the fields of gaseous water concentration resulting from evaporation across the various cases. The concentrations remain comparable across different initial droplet dimensions and are contingent solely upon variations in the equivalence ratio and flame propagation speed. The latter factor is particularly significant, as it influences the residence time of particles within the flame’s hot region (i.e., flame and burned gas region), as explained at the beginning of the present section. It can be observed that for the majority of cases, evaporation is insufficient to substantially impact temperature fields in a significant manner. The sole exception is noted in the case with
and the unity Lewis number. Nevertheless, minor cooling effects, which are also discernible in
Figure 3, may still influence flame speed and propagation characteristics.
It is well known that the flame speed is influenced by stretch effects (=strain+curvature effects) via the Markstein relation [
46]. Hence, all propagation-related properties are now plotted as functions of the average radius
. This approach aims to isolate the dependence of the propagation speed on the curvature similarly across all analyzed cases. The average radius is computed from the volume using the following formula:
The normalized total burning rate
S relative to the time of turbulence initiation is defined as follows:
In this way,
S can be simply interpreted as the relative evolution of the total burning rate. In the definition of
S, the stretched laminar flame speed
at the target radius for turbulence initiation and the initial projected flame area
, computed through the initial radius based on volume
, are used for the estimation of the initial total burning rate.
Figure 5 depicts the relative evolution of the total burning rate for cases with and without preferential diffusion, which are presented on the left and right sides, respectively. The presence of preferential diffusion notably enhances the growth of this quantity, particularly under lean conditions. The higher values of
S in lean conditions when
can be elucidated from the effective Lewis number values presented in
Table 2. At the same time, the lean case without preferential diffusion exhibits a notably slower evolution compared to all other cases, which also justifies the limited span of data along the abscissa in the plots. Furthermore, slight differences can be observed between the trends associated with evaporation when comparing cases with and without water addition. However, these deviations remain marginal due to the modest degree of evaporation witnessed in the simulated cases. This is attributed to the limited residence time of droplets within the flame and the limited time spent in the burned gas region. The former directly influences the flame structure, while the latter primarily has indirect effects related to the decrease in thermal expansion. Nonetheless, despite the observed effects being limited, discernible physical trends and phenomena are already apparent. The evolution of the total burning rate is influenced by variations in reactivity, which are reflected in the burning rate per unit area
, where
is the unstretched laminar burning velocity,
is the projected flame area at the instant of the evaluation, the turbulent flame area
, and the turbulent flame speed
. In line with the findings of several previous studies [
23,
44], instances characterized by pronounced preferential diffusion induced by
exhibit significantly greater values of both
and
compared to cases where this phenomenon is absent. Consequently, these characteristics are now observed in order to discern their responses to a non-zero curvature and the introduction of liquid water.
According to Damköhler’s first hypothesis, the value of
, as it is defined in the present work, is expected to equal
under adiabatic and equidiffusive premixed statistically planar flame conditions [
47,
48]. Statistically planar turbulent flames with a Lewis number smaller (larger) than one are characterized by a value of
(
) [
49]. Unity Lewis number flames with a globally negative mean curvature (e.g., Bunsen flames) are characterized by a value of
[
48], whereas flames with a globally positive mean curvature (e.g., flame kernels) are characterized by a value of
[
50]. For curved flames with a non-unity Lewis number, both effects compete with each other [
50]. In addition, as previously noted, the conditions under investigation in this study are non-adiabatic due to water evaporation, and notably, the effects of preferential diffusion are also considered. The values of burning rate per unit area, depicted in
Figure 6, illustrate that when
, the quantity
, reaching values up to
. The cases with
settle to values marginally larger than unity for
and smaller than unity for
. In order to explain these effects, we introduce the Markstein relation:
The displacement speed
is the velocity with which the flame propagates with respect to an initial coincident material surface [
51]. In this context,
denotes the Markstein length, which is typically positive in the absence of thermodiffusive instability, and
K signifies the stretch rate, which describes the variation of flame area over time. The stretch rate is a function of flame curvature, denoted as
, and tangential strain rate, represented by
. In the definitions of curvature and tangential strain rate,
signifies the normal direction of the flame front and is defined as
. The displacement speed, defined in Equation (
13), is a local quantity dependent on the reaction rate, diffusive fluxes, concentration gradients, and evaporation effects, where the latter two are generally negligible compared to the reaction rate and diffusive fluxes [
52].
Upon integration over the computational volume, the contribution of the diffusive flux vanishes [
48], leading to the following equation:
The volume integral of the reaction rate can also be expressed as
. Using the surface averaging introduced by Boger et al. [
53],
, where
denotes the volume average of a general quantity
Q, the left-hand side of Equation (
15) can be rewritten as:
By substituting the density-weighted Markstein relation (Equation (
13)) under the assumption of a constant Markstein length, the following expression is obtained, indicating that the sign of
depends on the product of the sign of the Markstein length and the sign of the surface-averaged stretch rate:
According to the definition of the Markstein length associated with the displacement speed given in [
18],
attains the values
for
, respectively, for
, whereas it remains larger than unity for
, i.e.,
for
, respectively. This explains the qualitative trends observed in
Figure 6. A quantitative assessment is complicated by the fact that Equation (
17) is only valid for small stretch rates; the Markstein length is not a constant (see, e.g., [
51]) and is very sensitive to the definition of the flame speed [
18].
The influence of water on the normalized burning rate per unit area is nearly negligible, as evaporation remains minimal. However, as demonstrated in prior works by the authors [
1,
3], significant cooling effects of water under substantial evaporation rates can considerably impact
, leading to values much lower than unity for both hydrogen and hydrocarbon combustion. Another factor for investigating the influence on flame propagation characteristics is the flame area.
Figure 7 illustrates the temporal evolution of the flame area for cases with preferential diffusion on the left and without preferential diffusion on the right. The trends of the normalized flame area are strongly correlated to those of the normalized total burning rate profiles, highlighting its significant influence on propagation speed.
Additionally, water evaporation marginally impacts the trend of flame area evolution, resulting in diminished area generation due to the cooling effect, which mitigates flame wrinkling induced by turbulence. Another noteworthy characteristic is observed in cases where preferential diffusion is active, where the slope of the curves grows faster with the radius than in the cases without preferential diffusion. In a spherical expanding flame configuration, area increase is influenced by two factors: flame outward expansion, evident in the increase in radius based on volume, and flame surface wrinkling. The intensity of the latter can be estimated by the ratio between the radius estimated from the volume and the radius estimated from the flame area (); when this ratio substantially exceeds unity, it suggests intense surface wrinkling.
Figure 8 illustrates the deviation of this ratio from unity. As anticipated, cases without preferential diffusion generally exhibit slower growth compared to those with preferential diffusion. However, it is noteworthy that the ratio
is higher for lean cases than for stoichiometric cases when
. This is due to the data being plotted over the radial evolution rather than time, with the duration of flame–turbulence interaction being the determining factor for flame wrinkling. Water influences the ratio
in a stronger way in cases with preferential diffusion than in cases with
.
Initial flame kernels are prone to quenching caused by curvature stretch effects and turbulence. Water injection has the potential to further reduce the reaction rate, as found by earlier studies [
1,
45]. Nevertheless, the present results suggest that due to the low evaporation rate of water, water droplets are unlikely to affect the initial flame kernel growth significantly and have a comparably small effect on the conditions for sustained combustion. Nevertheless, for longer interaction times, water droplets are likely to have the potential to reduce flame acceleration or turbulent flame speed for both the preferential diffusion and equidiffusive cases, as has been observed in previous works [
1,
45].