In the standard flamelet approach, the functions T ( ξ , t), Y k ( ξ , t ) define the flame structure if we assume that the structure of the diffusion flame depends only on the mixture fraction ξ and on time t . This is substantiated by making a variable change in the species equations from ( x 1 , x 2 , x 3 , t) to (z, y 2 , y 3 , t) where y 2 and y 3 are spatial variables in planes parallel to iso- ξ surfaces (Poinsot and Veynante, 2005). T = T ξ , t a n d Y K = Y k ξ , t (1)

Neglecting the terms associated with gradients along the flame front ( y 2 and y 3 ) and keeping the terms perpendicular to the flame (along ξ ) as significant we could approximate the flame structure locally as one-dimensional, depending only on time and the coordinate normal to the flame front ( ξ ) as presented in Equation 1. For this assumption to hold the flame must be thin compared to other flow and wrinkling scales in a multi-dimensional flow (Poinsot and Veynante, 2005). This implies that each part of the flame front can be approximated as a small laminar flame, the “flamelet”. Following the above the species mass fraction equation as well as the temperature are written as follows in Equations 2, 3: ρ ∂ Y k ∂ t = ω ̇ k + 1 2 ρ χ ∂ 2 Y k ∂ ξ 2 (2) ρ ∂ T k ∂ t = ω ̇ T + 1 2 ρ χ ∂ 2 T ∂ ξ 2 (3)

The only term depending on spatial variables ( x i ), in the flamelet Equations 2, 3 is the scalar dissipation rate χ which controls mixing as it dictates the gradients of ξ (Poinsot and Veynante, 2005). Once χ is specified, Equations 2, 3 can be solved in the ξ space to provide the flame structure. The usual way to determine its value/distribution is by performing a direct numerical simulation of the problem under consideration or an experiment. If this is not the case, a try error process is followed in order to allocate what is called the extinction scalar dissipation rate over which the flame cannot exist.

The flamelet equations can also be understood as a subset of the Conditional Moment Closure approach, reviewed in Bilger (1993), where the conditioning is performed based on the mixture fraction (Peters, 2001). Rather than focusing on flame surface statistics and the laminar reactive-diffusive structure tied to the flame surface as is the case in the flamelet model, CMC employs conditional moments at a particular location x and time t within the flow field, using the conditional probability density function (Peters, 2001). The flame structure is characterized by the conditional averages of the reactive scalars and the conditional mean Q a is defined as the ensemble average of Y a ( x , t ) under the condition that the associated value ξ ( x , t ) is at the value η (Philip, 2019) as seen in Equation 4: Q α η ; x , t ≡ 〈 Y α x , t ∣ ξ x , t = η 〉 ≡ 〈 Y α ∣ η 〉 (4)

Through this averaging process, the instantaneous reactive scalar value can be recovered into the conditional mean and a deviation from the conditional mean as presented in Equation 5 (Bilger, 1993): Y α x , t = Q α ξ x , t ; x , t + Y ″ α x , t (5)

The core concept of the CMC method involves solving transport equations for the conditional averages of reactive scalars in order to determine the flame structure (Philip, 2019). Furthermore, an exact transport equation for these conditional moments can be derived. That is by adopting a set of assumptions, including high Reynolds number, the primary closure hypothesis, molecular-level Fickian mass diffusion, a Lewis number of unity, and negligible conditional density fluctuations, the equation can be simplified to the form of Equation 6 (Klimenko and Bilger, 1999): ∂ Q α ∂ t + 〈 u ∣ η 〉 ⋅ ∇ Q α = 〈 N ξ ∣ η 〉 ⁢ ∂ 2 Q α ∂ η 2 + 〈 ω ̇ α ∣ η 〉 − ∇ ⋅ 〈 ρ ∣ η 〉 〈 u ″ Y ″ α ∣ η 〉 p η 〈 ρ ∣ η 〉 p η (6)

The unconditional mean of a variable can be determined from its conditional means as follows: ⟨ Y α x , t ⟩ = ∫ 0 1 Q α η ; x , t p η ; x , t   d η (7)

The DCMC equations presented in this sub-section are adapted from the work of Philip (2019). The rationale for using double conditioning is that if most fluctuations around the single-conditional mean can be linked to the second conditional variable, then the fluctuations around the double-conditional mean, determined by two mixture fractions, should be small (Philip, 2019). As singly-conditioned CMC models are limited to cases that are either premixed or predominantly non-premixed (Philip Sitte and Mastorakos, 2019) the double conditioning allows for a wider description of phenomena present in dual fuel systems. The DCMC method extends the parameter space and incorporates the cross-scalar dissipation rate, while remaining a straightforward generalization of the CMC approach. Additionally, it is expected to describe well all combustion modes relevant to dual-fuel engines (Harikrishnan et al., 2024). The doubly conditional space is as follows in Equation 8: D = η 1 , η 2 ∈ R 2 : 0 ≤ η 1 ≤ 1 , 0 ≤ η 2 ≤ 1 (8) The conditional mass fractions are given in Equation 9 and the conditional enthalpy in Equation 10. The conditional temperature equation also exists but only one of the temperature and enthalpy need to be solved choosing the latter: ∂ Q α ∂ t + ∇ ⋅ Q α 〈 u ∣ η 1 , η 2 〉 = Q α ∇ ⋅ 〈 u ∣ η 1 , η 2 〉 − 1 ρ ̄ p ̃ ∇ ⋅ ρ ̄ p ̃ 〈 u Y α ∣ η 1 , η 2 〉 − Q α 〈 u ∣ η 1 , η 2 〉 + 1 ρ ̄ p ̃ ∇ ⋅ ρ ̄ p ̃ 〈 D ∇ Y α ∣ η 1 , η 2 〉 + 〈 N ξ 1 ∣ η 1 , η 2 〉 ∂ 2 Q α ∂ η 1 2 + 2 〈 N ξ 1 ξ 2 ∣ η 1 , η 2 〉 ∂ 2 Q α ∂ η 1 ∂ η 2 + 〈 N ξ 2 ∣ η 1 , η 2 〉 ∂ 2 Q α ∂ η 2 2 + 〈 ω ̇ α ∣ η 1 , η 2 〉 − 〈 ω ̇ c * ∣ η 1 , η 2 〉 ∂ Q α ∂ η 2 (9) ∂ Q h ∂ t + ∇ ⋅ Q h ⁢ 〈 u ∣ η 1 , η 2 〉 = Q h ⁢ ∇ ⋅ 〈 u ∣ η 1 , η 2 〉 − 1 ρ ̄ p ̃ ⁢ ∇ ⋅ ρ ̄ ⁢ p ̃ ⁢ 〈 u h ∣ η 1 , η 2 〉 − Q h ⁢ 〈 u ∣ η 1 , η 2 〉 + 1 ρ ̄ p ̃ ⁢ ∇ ⋅ ρ ̄ ⁢ p ̃ ⁢ 〈 a ∇ h ∣ η 1 , η 2 〉 + 〈 N ξ 1 ∣ η 1 , η 2 〉 ⁢ ∂ 2 Q h ∂ η 1 2 + 2 〈 N ξ 1 ξ 2 ∣ η 1 , η 2 〉 ⁢ ∂ 2 Q h ∂ η 1 ∂ η 2 + 〈 N ξ 2 ∣ η 1 , η 2 〉 ⁢ ∂ 2 Q h ∂ η 2 2 − 〈 ω ̇ c * ∣ η 1 , η 2 〉 ⁢ ∂ Q h ∂ η 2 + 1 ρ ⁢ ∂ p ∂ t | η 1 , η 2 (10) where, according to Harikrishnan et al. (2024), the following are true: ⟨ N ξ 1 ∣ η 1 , η 2 ⟩ , ⟨ N ξ 2 ∣ η 1 , η 2 ⟩ , and ⟨ N ξ 1 ξ 2 ∣ η 1 , η 2 ⟩ are the conditional scalar dissipation rates (SDRs) of ξ 1 , ξ 2 , and the cross-scalar dissipation rate, respectively. The final term in the Q h equation accounts for the pressure ( p ) rise in the constant volume chamber. The term ⟨ ω ̇ α ∣ η 1 , η 2 ⟩ represents the conditional chemical source term while the ⟨ u ∣ η 1 , η 2 ⟩ , p ̃ , and ρ ̄ are the conditional velocity, Favre-filtered probability density function, and filtered density, respectively. Moving on, η 1 and η 2 are the respective sample space variables for the two mixture fractions. Finally, ⟨ ω ̇ c * ∣ η 1 , η 2 ⟩ is the conditional apparent reaction rate term. The apparent reaction rate, encapsulating both the direct chemical source term and effects of the effects of the scalar dissipation rates, is presented in Equation 11: ⟨ ω ̇ c * ∣ η 1 , η 2 ⟩ = 1 ∂ Q ψ ∂ η 2 ⟨ ω ̇ ψ ∣ η 1 , η 2 ⟩ + ⟨ N ξ 1 ∣ η 1 , η 2 ⟩ ∂ 2 Q ψ ∂ η 1 2 + 2 ⟨ N ξ 1 ξ 2 ∣ η 1 , η 2 ⟩ ∂ 2 Q ψ ∂ η 1 ∂ η 2 + ⟨ N ξ 2 ∣ η 1 , η 2 ⟩ ∂ 2 Q ψ ∂ η 2 2 (11)

The conditional expectation of the mixture fraction SDR, ⟨ N ξ ∣ η 1 , η 2 ⟩ , presented in Equation 12 is modeled using Amplitude Mapping Closure (AMC) (O’Brien and Jiang, 1991) of Equation 13. The shape of this bell curve plays a significant role in the diffusion in mixture fraction space due to the zero value it assumes at the boundaries not allowing diffusion to those points from internal values of the mixture fraction. ⟨ N ξ ∣ η 1 , η 2 ⟩ = N 0 G ξ η 1 (12) where, G ξ η = exp − 2 erf − 1 2 η − 1 2 (13) N 0 = N ̃ ξ ∫ 0 1 ∫ 0 1 G ξ η 1 p ̃ η 1 , η 2   d η 1   d η 2 (14) N ξ 1 ≡ D ∇ ξ 1 ⋅ ∇ ξ 1 , N ξ 2 ≡ D ∇ ξ 2 ⋅ ∇ ξ 2 , N ξ 1 ξ 2 ≡ D ∇ ξ 1 ⋅ ∇ ξ 2 (15) where N 0 = ⟨ N ξ ∣ 0.5 , η 2 ⟩ and ⟨ N ξ ∣ η 1 , 0.5 ⟩ respectively using Equation 14. The modeling of the cross-scalar dissipation rate is still a matter of research and constitutes a central input to the DCMC method. The modelling of the latter term of Equation 15 is proposed in Equation 16 by Kronenburg (2004) is utilized. ⟨ N ξ 1 ξ 2 ∣ η 1 , η 2 ⟩ = N ̃ ξ 1 ξ 2 N ̃ ξ 1 N ̃ ξ 2 ⟨ N ξ 1 ∣ η 1 , η 2 ⟩ ⟨ N ξ 2 ∣ η 1 , η 2 ⟩ (16) The ratio N ̃ ξ 1 ξ 2 N ̃ ξ 1 N ̃ ξ 2 will be denoted by N rat in this paper and its magnitude will be varied to examine the consequences to the solution. For partially premixed combustion, the cross-scalar dissipation rate is an important quantity and assumes a positive value for premixed combustion and a negative for non-premixed combustion (Chakraborty et al., 2007). There is very little information in the literature for this quantity whose modelling is a topic for further research.

To initially explore flame structures in dual-fuel systems, as in this paper, a simplification process is followed reducing the DCMC equations to zero dimension and essentially going from three spatial and two mixture fraction dimensions to only two mixture fraction dimensions. The underlying assumptions of the 0D-CMC equations are.

• High Reynolds number that allows for the neglect of some molecular transport terms.

Under these assumptions, the 0D-DCMC equations assume the form of Equations 17, 18. ∂ Q α ∂ t = 〈 N ξ 1 ∣ η 1 , η 2 〉 ∂ 2 Q α ∂ η 1 2 + 2 〈 N ξ 1 ξ 2 ∣ η 1 , η 2 〉 ∂ 2 Q α ∂ η 1 ∂ η 2 + 〈 N ξ 2 ∣ η 1 , η 2 〉 ∂ 2 Q α ∂ η 2 2 + 〈 ω ̇ α ∣ η 1 , η 2 〉 (17) ∂ Q h ∂ t = 〈 N ξ 1 ∣ η 1 , η 2 〉 ∂ 2 Q h ∂ η 1 2 + 2 〈 N ξ 1 ξ 2 ∣ η 1 , η 2 〉 ∂ 2 Q h ∂ η 1 ∂ η 2 + 〈 N ξ 2 ∣ η 1 , η 2 〉 ∂ 2 Q h ∂ η 2 2 + 〈 Q ̇ ∣ η 1 , η 2 〉 (18)

In this paper, 0D-DCMC simulations will be performed to understand the basic structure of dual-fuel autoignition systems, as a first step in a full implementation in a turbulent reacting flow multi-dimensional simulation.

In this section, the initial conditions of the solution of the equations are discussed. Two mixture fractions are used, one for each fuel, ξ 1 for decane, and ξ 2 for ammonia. The initial conditions are chosen as mixing lines. These lines correspond to states where the fuel and oxidizer mix without reacting. Such mixing lines are crucial in capturing ignition and quenching (Poinsot and Veynante, 2005).

With the presence of two fuels, it is not clear what to do in the top right corner. The points ( η 1 , η 2 ) = ( 0,1 ) and ( η 1 , η 2 ) = ( 1,0 ) correspond to pure ammonia and pure diesel respectively whereas the point ( η 1 , η 2 ) = ( 0,0 ) to the pure oxidizer (air). The point ( η 1 , η 2 ) = ( 1,1 ) is a matter of debate as there are multiple options and it is unclear whether this results in different results. As we have two fuels the logical boundary condition is to consider that at that point 50 % ammonia and 50 % diesel (by mass) are present corresponding to a mixture. The resulting distributions are seen in Figures 1, 2. One could also set this point as 100 % ammonia or diesel altering the mathematical formulation of the problem and thus the results. However, the quantitative analysis yields the same conclusions. In conclusion, the current choice is a matter of modelling and is not expected to affect the results if the interpretation is done correctly considering the initial conditions.

As the chemistry is detailed and the c p of each species is different, the initial distributions of species mass fractions and the enthalpy are linear functions of ξ 1 and ξ 2 . However, the temperature is not. In the current analysis, the mass fractions and the enthalpy are used and the temperature distribution is computed from the enthalpy. The same typical approach as in the single mixture fraction case is followed for constructing the mixing lines in each of the four boundaries interpolating between the extreme values at the two corners corresponding to each boundary. The values at the inner nodes are computed from a multivariate interpolation.

The 0D-DCMC equation is a parabolic PDE in two dimensions (the two mixture fractions) and time. In a full 3D simulation in the future, we would add the 3 spatial dimensions, as in Philip, (2019), but this is not done here. The mixture fraction spaces are discretised with a 71 x 71 nonuniform grid covering the whole range of 0–1. The choice of its size is dictated by the trade-off between accuracy and computational effort. The grid is denser around the stoichiometric mixture fraction of each fuel as the flame and most reactions are expected to lie there. Due to the steep gradients and the rapid construction and destruction of species around stoichiometry, a very dense grid is required around there for accuracy. As a result, using a non-uniform grid allows for significant computational improvements. The diffusion in mixture fraction space is solved first with the VODPK solver (Alan, 1983), using the non-stiff Adam’s method which solves the set of non-linear ODEs. Central-differencing schemes for non-uniform grids, which are second-order accurate, are used for the numerical discretisation of the ∂ 2 Q ∂ ξ 2 terms.

Chemical source terms are addressed last using the SpeedCHEM solver (Perini et al., 2012). This solver, referred to as LIBSC, utilizes a sparse Jacobian formulation and precomputed, temperature-dependent properties to greatly enhance computational speed. It also supports the use of an analytical Jacobian, which approximates the Jacobian matrix of the ODE system. For most chemical mechanisms, this approach offers significant performance enhancements. Timesteps equal or smaller than 10 − 6 s are used to minimise operator splitting errors.

This section briefly discusses the selection of the chemical mechanism used for the simulations. Typically, the physical properties of diesel, such as density and evaporation rate, are modeled using decane or dodecane, while heptane is used to mimic diesel’s chemistry. Most current chemical mechanisms for ammonia/hydrocarbon fuels are derived by merging databases for each fuel and as a result in general further research is required in this area.

A recently published chemical scheme, the Aalto75 model (Tolga Kurumus et al., 2024), includes 75 species and 451 reactions and is reported to be well-validated across various experimental datasets, including ignition delay times, laminar flame speeds, and species profiles. However, it was not used in this study as it was released after the project had already begun and is left for future work. Another chemical mechanism for ammonia/n-heptane mixtures was published by (Thorsen et al., 2023), offering good predictions of ignition delay times at pressures up to 100 bar. However, this mechanism was not selected due to its complexity, as it has too many species and reactions, making it very difficult to handle computationally for the current project.

The chemical mechanism ultimately chosen for this study is the one provided by Lin Tay et al. (2017) for decane-ammonia combustion, consisting of 80 species and 374 reactions. The ignition delay times for diesel and ammonia were found to be very similar with the parent mechanisms, indicating that this mechanism is sufficiently robust for use in diesel engine simulations and should contain the key aspects of the problem.

In this section, homogeneous reactor calculations with the use of Cantera are performed. The motivation is twofold: to validate that the chemical mechanism produces reasonable results for homogeneous mixtures and to locate the “most reactive” mixture fraction ξ m r (Mastorakos, 2009) in the case of fuel-air mixing layers. For such cases, the initial condition of temperature is connected with the initial mass fractions through Equations 19–22: Y f u = ξ × Y f u 0 (19) Y O 2 = 1 − ξ × Y O 2 0 (20) Y N 2 = 1 − ξ × Y N 2 0 (21) T = T o x − ξ × T o x − T f u (22)

The temperature of the oxidizer is left as a parameter and the fuel temperature is set constant to 300 K. The calculations are based on a series of batch reactors for the different temperatures and pressures with the mixture fraction being varied with the initial condition as above. For each temperature or pressure, the most reactive mixture fraction is the one producing the smallest ignition delay.

A typical batch reactor simulation for the ignition of homogeneous ammonia/air mixtures is presented in Figure 3 for an oxidiser temperature of 1200 K and a pressure of 2 atm and it is found that ignition takes place at approximately 1.1 ms. We notice that the rapid temperature increase is correlated with the rise of the species OH and H 2 and that N 2 O is rapidly produced and destroyed making it a pre-ignition species, as seen in the right axis of the mass fractions. N 2 O, HNO, NO and NH 2 are quickly produced at the autoignition instant and then quickly destroyed. On the contrary, species like NNH and OH, although again rapidly produced, stay relatively constant after the ignition. Finally, observing the early chemistry we see that NO 2 is produced during the species accumulation period and flattens during the autoignition.

For the case of ammomia/air mixing layers, i.e., where the initial condition of species and temperature are the inert mixing between pure ammonia and pure air, Figures 4, 5 show the variation of ignition delay time with mixture fraction for different pressures and temperatures. These curves serve to locate the most reactive mixture fraction which is the minimum of those U-bell shapes. For low values of mixture fraction, despite the high temperatures, the mixture is too lean while for high mixture fractions, the temperature is low and the mixture too rich, and hence the minimum autoignition time is somewhere between and this characterises the most reactive mixture fraction. The most reactive mixture fraction ξ m r for the given pressure of 20 bar and an oxidizer temperature of 1200 K is 0.11. For reference, the stoichimetric mixture fraction is ξ s t = 0.142 . It is also evident that as the pressure or temperature increase, the most reactive mixture fraction determination becomes less sharp, which suggests that autoignition in non-premixed systems may be occurring almost simultaneously in wide regions in space (i.e., at a range of mixture fractions). Nevertheless, the estimate of ξ m r and the associated autoignition time are important reference values for inhomogeneous systems and can be compared when straining (e.g., through high values of the scalar dissipation) becomes important.

The same procedure is followed for the decane-air mixture. In Figures 6, 7 ignition delay in terms of the mixture fraction ξ m r is presented for different pressures and temperatures. Again, similar trends are observed for lower values of ξ as expected from the stoichiometric ξ s t = 0.0626 . Compared to ammonia, decane ignites quickly for high temperatures and pressures. For example, for 20 bar and T o x = 1200 K, the quickest ignition time for the ammonia-air system is about 35 ms, while for the decane-ar system the quickest ignition delay is about 10 ms. Overall, the higher reactivity of decane is demonstrated in the chemical scheme used. Therefore, the concept that the high reactivity fuel (decane) will autoignite before the ammonia is captured well by the current mechanism.

In this section, the 0D-DCMC exploration takes place to understand the effects of the scalar dissipation rates on the ignition of the decane-ammonia mixture. An overall description is given first, followed by a more detailed presentation of the structure of the autoignition spot and the evolution of the reaction fronts. The sensitivity to scalar dissipation rates is also shown.