A Eulerian Multi-Fluid Model for High-Speed Evaporating Sprays

: Advancements in internal combustion technology, such as efﬁciency improvements and the usage of new complex fuels, are often coupled with developments of suitable numerical tools for predicting the complex dynamic behavior of sprays. Therefore, this work presents a Eulerian multi-ﬂuid model specialized for the dynamic behavior of dense evaporating liquid fuel sprays. The introduced model was implemented within the open-source OpenFOAM library, which is constantly gaining popularity in both industrial and academic settings. Therefore, it represents an ideal framework for such development. The presented model employs the classes method and advanced interfacial momentum transfer models. The droplet breakup is considered using the enhanced WAVE breakup model, where the mass taken from the parent droplets is distributed among child classes using a triangular distribution. Furthermore, the complex thermal behavior within the moving droplets is considered using a parabolic temperature proﬁle and an effective thermal conductivity approach. This work includes an uncertainty estimation analysis (for both spatial and temporal resolutions) for the developed solver. Furthermore, the solver was validated against two ECN Spray A conditions (evaporating and non-evaporating). Overall, the presented results show the capability of the implemented model to successfully predict the complex dynamic behavior of dense liquid sprays for the selected operating conditions.


Introduction
Although fuel spray modeling is not a new or unknown problem, continuous advancements in efficiency improvements and new fuels require new modeling solutions with improved functionality and accuracy. Therefore, numerical tools for predicting the complex, dynamic behavior of sprays impact the development of improved internal combustion (IC) engines. The utilization of numerical simulations greatly affects the development and improvement of multiple other complex engineering challenges, e.g., furnaces [1], fluidized beds [2], fuel cells [3], and carbon capture and storage [4].
As previously mentioned, the two primary development efforts in IC technology are improving the efficiency of engines and developments regarding new fuels and engines which can run on them. Efficiency improvements are mainly related to the compression ignition regime and increasing the compression ratio [5], which often contributes to an increase in nitrogen oxide emissions. The partially premixed combustion approach [6] is one of the various strategies that should allow high efficiencies and minimize emissions. Even the spark-ignited engines have several methods for reducing emissions and increasing efficiencies, e.g., the corona ignition condition [7], where the electron dissociation reaction replaces the classic hydrocarbon oxidation, and partial fuel stratification combustion strategies [8]. Optimal fuel-air mixing is a prerequisite for emissions reduction, and is directly coupled to the flow conditions occurring in the nozzle. Injection velocity controls the atomization and droplet breakup, directly affecting the liquid core penetration and evaporation rate.
Modern engines are also required to operate on a broad range of different fuels, for example, biofuels [19], solar fuels [20] (where solar energy is stored as synthetic chemical fuels, i.e., methanol or ethanol), and "smart" fuels which have additives for improving the properties of the fuel (e.g., ignition control [21] and reduction of emissions [22]). Furthermore, surrogate fuels [23,24] are developed to mimic the selected properties of targeted real fuels. Increasing demand in the heavy-duty market resulted in the utilization of naphtha and heavy fuel oils as alternative fuels [25,26].
Computational fluid dynamics (CFD) proposes numerous multiphase approaches which are suitable for modeling liquid sprays. The most detailed and computationally most expensive is the direct numerical simulation (DNS) method [27][28][29][30]. DNS does not need any sub-models to reproduce detailed spray behavior, such as atomization, secondary breakup, and turbulence coupling between the liquid and gas phases. DNS is still not viable for standard engineering calculations, but it is constantly gaining popularity due to the increasing HPC resource availability.
The most popular approach is the Lagrangian method, where the group of droplets (parcels) are traced in a Lagrangian manner. In contrast, the continuous phase is represented as a continuum in the Eulerian coordinates. Despite its popularity and multiple benefits, this approach has some disadvantages, e.g., problems with modeling the near-nozzle region where the spray is denser [27,31], pronounced sensitivity to cell resolution [32], and problems with numerical instabilities [33]. Furthermore, the droplet phase and the continuous phase are calculated in a decoupled manner.
The third approach-the Euler-Euler approach-describes both the droplet disperse phase and the continuous phase as interpenetrating continua in the Eulerian coordinates. Therein, all phases are described using properly derived conservation equations [34,35], but some of the small-scale phenomena are neglected due to the averaging procedure. This method has various formulations, but this study focuses on the Eulerian multi-fluid model. This approach couples the population balance equation (PBE) [36] with the averaged momentum and continuity equations, which allows modeling of polydisperse flows. For the discretization of the PBE, a classes method is employed. The PBE is divided into an arbitrary number of droplet classes with predefined diameters. Consequently, all droplet classes have their phase continuity and momentum equations, which introduce higher precision due to the allowed spatial and temporal variance in velocity and interfacial momentum transfer models (highly dependent on the dispersed element size). The mixture continuity equation is used to derive the mixture pressure equation, and all phases, including the single continuous phase, share the same mixture pressure. The higher precision and resolution of the results come at a price. The numbers of equations and calculations that need to be executed are directly coupled to the selected number of droplet classes. Particular problems encountered by the engineering and academic communities require different levels of accuracy and execution speed. Therefore, all three aforementioned approaches are actively used and developed [27,28,[37][38][39][40][41][42][43][44].
This research presents the upgrade of the previously developed and published Eulerian multi-fluid model (specialized for dense spray applications) with evaporation capability. The developed opensource framework can now successfully predict the dynamic behavior of dense evaporating liquid fuel sprays, i.e., the atomization and secondary breakup in the dense near-nozzle region, and the evaporation and mixing phenomena which dominantly occur in the dilute part of the spray. This work is a significant upgrade to the previously developed and published model [45][46][47]. The addition of an energy conservation equation, species transfer, and density variance allow the implementation of an evaporation model. Furthermore, to capture the finite thermal conductivity of droplets, the model also includes a parabolic temperature profile within droplets, and to account for the internal flow within the droplets, an effective thermal conductivity approach is employed. The developed model was tested and validated with the ECN Spray A experimental measurements [48][49][50]. The development was done within the OpenFOAM library. To the authors' knowledge, such a detailed modeling approach has not yet been applied with the Eulerian method using the multi-fluid formulation. The proposed model allows straightforward upgrades in the future, allowing simulation of even more complex fuel behavior (e.g., pronounced multi-component behavior).
The remainder of the paper is organized into four sections. Section 2 gives a detailed description of the mathematical model, i.e., updates and new implementations performed within the previously developed model. Section 3 presents the employed numerical procedure. Section 4 introduces the selected test cases and presents the numerical results, including uncertainty analysis and two validation cases for evaporating and nonevaporating conditions. Section 5 provides a conclusion and comments regarding the implemented model.

Mathematical Model
This work is a significant upgrade to the previously developed and published model, which was employed to predict the dynamic behavior of bubbly flows and non-evaporating dense sprays [45][46][47].
In this work, the model is further upgraded with evaporation handling, which requires implementing a species transfer equation and energy equations for all phases. Due to evaporation, the continuous phase exhibits significant changes in chemical composition, which directly influence the thermophysical properties. Substantial temperature differences also impact the thermophysical properties of the liquid fuel. Therefore, the phase continuity equations and momentum equations were re-implemented into a compressible formulation. Furthermore, the mass transfer due to droplet breakup is distributed among multiple child classes using a triangular distribution. The previously developed turbulence model is updated with an algebraic model, which improves the coupling of the droplet turbulence variables with the continuous phase turbulence.
This section presents the Eulerian multi-fluid model, which is specialized for highspeed evaporating sprays. The following sub-sections emphasize the updates and new developments within the described solver. The details regarding the employed averaging procedures, pressure-velocity coupling, and phase-intensive continuity and momentum equations can be found in [35,[45][46][47]51,52]. The following sub-sections use the finite volume notation of Weller [51].

The Phase-Intensive Momentum Equation
Following the approach presented by [51], the compressible phase-intensive momentum equation can be generalized for a multi-fluid formulation: where U ϕ is the density-weighted ensemble averaged phase velocity (for phase ϕ), p gives the mixture pressure shared by all phases, ρ ϕ is the phase density, α ϕ is the phase fraction, R eff ϕ gives the viscous and turbulent stress, and g denotes the gravitational acceleration. Vectors M ϕ and S Mϕ denote the averaged interfacial momentum transfer term and net momentum source term due to mass transfer caused by droplet breakup and evaporation.
The interfacial momentum transfer term M ϕ is responsible for exchanging momentum between the continuous gas phase and droplet classes. The momentum exchange is taken into account using the turbulent dispersion and drag. The momentum exchange term for the i-th droplet class reads: In Equation (2), the subscript c denotes the continuous phase, and the subscript d, i signifies the i-th droplet phase. The U r,i term denotes the relative velocity, d i is the droplet diameter, k c is turbulence kinetic energy of the continuous phase, and C td,i gives the turbulent dispersion coefficient, which is calculated using the approach described by [53,54]. The droplet drag coefficient C d,i is implemented following the procedure described by [55], which takes into account large deformations of droplets: where the droplet drag is blended between the drag of an ideal sphere C d,sphere,i and a disc, which has a 3.6 times larger drag. The blending factor y i , i.e., the normalized distortion parameter, is evaluated using the Taylor-Analogy (TAB) model [56]. The drag of an ideal sphere C d,sphere,i is evaluated using the correlation given by [57], which takes into account the local phase fraction correction for dense spray regions: where Re d,i gives the Reynolds number for the i-th droplet class. The total interfacial momentum transfer to the continuous phase is given by: where n droplets is the selected number of droplet classes.

Phase Continuity Equation
Following the approach presented by [51], the compressible phase continuity equation is generalized for a multi-fluid formulation: where U r,i,j gives the relative velocity between phases i and j, and U denotes the ensemble averaged mixture velocity: and d i ρ i dt is given by: In Equation (6), the net source term S i describes the mass transfer between the phases. In this work, the mass exchange between phases is taken into account using the breakup and evaporation model. The net mass exchange term for the i-th droplet class is given by: where B B,d,i gives the droplet birth rate, and D B,d,i denotes the droplet death rate due to breakup. The D E,d,i term gives the droplet death rate due to evaporation. The net source term for the continuous phase is given by:

Breakup Model
Using the previously developed functionality [46], the aerodynamic interaction between the gas and high-speed liquid phases is taken into account using the WAVE breakup model [58][59][60][61][62]. The interaction develops and increases disturbances on the droplet's surface, which eventually lead to the breakup of parent droplets into smaller child droplets. The size of child droplets (predicted stable radius) r s is directly proportional to the wavelength Λ: In this work, the proportionality coefficient B 0 is taken to be 0.61. The mass loss per unit volume (death rate) D B,d,i of the parent class i due to the breakup process is given by: where the rate of change of the i-th class radius dr i dt is evaluated as: In Equation (13), τ i denotes the breakup time of the i-th phase: where B 1 is a model constant and Ω i denotes the estimated wavelength. Further details regarding the numerical implementation of the WAVE model for the multi-fluid formulation are given in [46]. In this work, the WAVE model is further upgraded by adding a probability distribution for the mass transfer to the child droplet classes. The implementation was done following the procedure described by [63]. In [63], the authors presented that this approach prevents the Eulerian multi-fluid model from overestimating the generation of small droplets. The distribution function smears the transfer of mass to multiple child droplet classes. This work employs the triangular distribution function, which is straightforward to implement because all parameters are directly available. Figure 1 gives an example of the approach employed for distributing mass transfer calculated with the WAVE breakup model. The mass taken from the parent class is distributed among the smaller droplet classes using a triangular distribution. The parent class's lower boundary gives the upper limit b, and the lower limit a is given by the lower boundary of the smallest class. The distribution's mode value c is given by the target diameter, estimated by the WAVE model d target = 2r s . In Figure 1 the area given by the blue bars represents the weighting factors of individual child classes. To satisfy the mass conservation criterion, the sum of all weighting factors must be unity. In Figure 1 the red bar is given only to represent the location and bounds of the parent droplet class. The probability density function (PDF) of the selected triangular distribution is given by:

Turbulence Model
Due to the relatively low computation cost and reasonable accuracy, the standard kmodel [64] is employed to account for turbulence effects within the simulation. To improve the model's predictive capabilities for spray applications, the model utilizes the round jet correction of Pope [65], where C 1 = 1.6. The droplet phase turbulence is coupled to the continuous phase turbulence values using an algebraic model [66,67]: where k d,i is the turbulent kinetic energy of the droplet phase, and the frequency ω i is given by: and the relaxation time τ i : where ν c denotes the kinematic viscosity of the continuous phase. The characteristic macroscopic length scale L x can be written as: where C µ is the turbulence model constant, and c denotes the rate of dissipation of turbulent kinetic energy (of the continuous phase). Furthermore, the droplet phase eddy viscosity ν t d,i is evaluated using the following expression: where ν t c is the eddy viscosity of the continuous phase.

Species Transfer
The fuel vapor species transfer equation is given by: where Y 1 is the mass fraction of the fuel vapor in the continuous phase, D Y 1 is the binary diffusion coefficient of fuel vapor, Sc t is the turbulent Schmidt number, and S Y 1 gives the source term due to the evaporation of droplet classes.

Energy Equation
Following the approach given by [67], the energy equation for the continuous phase is given by: where h c denotes the static enthalpy of the continuous phase, c p,c is the specific heat capacity, and κ eff c is the effective thermal conductivity. The net enthalpy source term S h,c describes the heat rate supplied to the gas phase, including both the convective heat transfer and evaporation. The droplet phase energy equation is given by: where S h,d,i accounts for the energy transfer due to droplet breakup, convective heat transfer, and evaporation. Equation (23) does not include the energy diffusion term because a parabolic temperature profile model [68] considers the finite thermal conductivity within the droplet classes. Following the guidelines described in [68], the droplet surface temperature is calculated as: whereṙ d,i is the derivative of the droplet radius (with respect to time), L d,i is the latent heat of vaporization, and ζ is given by: ψ = 1 + 0.2ζ, and T d,i is defined as: In Equation (25), Nu is the Nusselt number. Furthermore, to take into account the intensive recirculation occurring inside the moving droplets, which is caused by the surface friction, an effective conductivity approach is adopted. Abramzon and Sirignano [69] suggested a practical approach where the droplet thermal conductivity is replaced by effective conductivity: where the factor χ is calculated as: χ = 1.86 + 0.86 tanh 2.245 log 10 (Pe d,i /30) .
In Equation (27), Pe d,i denotes the Peclet number of the i-th droplet phase. The implemented parabolic temperature profile and effective conductivity model allow a more accurate prediction of the droplet thermal behavior [70]. Better prediction of the droplet temperature improves the evaluation of droplet surface properties that directly influence the phenomena occurring on the droplet's surface, e.g., evaporation.
The temperature-dependent thermophysical properties of the liquid fuel and vapor were evaluated using the expressions available in [71,72].

Evaporation Model
For modeling the droplet evaporation process, a hydrodynamic approach is employed. This approach assumes that the fuel vapor near the droplet surface is saturated all the time. Therefore, the droplet evaporation rate is equal to the rate of the vapor diffusion going from the droplet surface to the surrounding gas [73]. This approach focuses more on modeling the diffusion process than the detachment process of molecules from the droplet surface, which is much more challenging. This work utilizes the Abramzon and Sirignano hydrodynamic model [69]. The authors employed the film theory to consider the convective transport caused by the relative velocity between the droplet and ambient gas. Furthermore, the model is applicable for non-unitary Lewis number (in the gas film) cases.
The proposed model gives the following relations for the instantaneous rate of droplet evaporation:ṁ In this model, the averaged variables, i.e., the ones with overlines, are evaluated at a reference fuel vapor concentration and temperature using the "1/3 rule" [69]. In Equation (29), Sh denotes the dimensionless Sherwood number and B M gives the Spalding mass transfer number: In Equation (30), c pF is the average specific heat capacity of vapor in the film and B T denotes the Spalding heat transfer number: The non-dimensional parameters are defined as [69]: and where the correction factors F M and F T are given by: and The non-dimensional parameters for non-evaporating droplets (Nu 0 and Sh 0 ) are implemented as [69,74]: and where: Assuming ideal mixing and by applying Raoult's law, the mass fuel vapor fraction at the droplet surface is evaluated as [71,73]: where M F and M c give the molar masses of vapor and surrounding gas. However, neglecting real gas behavior (at high temperatures and pressures) introduces deviations because the ideal behavior assumption does not consider molecular interactions (attractive and repulsive forces) [15,73,75]. In Equation (40), the p Fs term denotes the saturation vapor pressure, and for n-dodecane, it can be estimated using the following expression [73]: During the execution of the solution procedure for the implemented evaporation model, F T is estimated using the old (from the previous iteration or time-step) B T value. The new B T value is evaluated using: where Φ is given by: In Equation (43), Le indicates the dimensionless Lewis number. Both Nu * and B T are re-evaluated until the difference |B new T − B old T | is below the desired accuracy [69]. The heat penetrating into the droplet is calculated as: The droplet death rate D E,d,i due to evaporation can be evaluated using the following expression [76]: and similarly, the heat penetrating the droplet class is given by:

Numerical Approach
The described model was implemented in foam-extend, a community-driven fork of OpenFOAM. Hence, the implemented equations were discretized using a collocated cellcentered finite volume method (FVM) [77,78]. The developed solution algorithm employs the PISO loop [79] and the selected solution procedure for each time step is presented in Algorithm 1. In this work, the mass fraction, turbulence variables, phase fraction variables, and energy variables were advected using a linear upwind-biased approximation. The momentum variables used the gamma scheme [80]. Time derivatives employed the first-order implicit Euler scheme. Laplacians, gradients, and cell-to-face interpolations were calculated using a linear interpolation. The solution of the pressure equation employed the selection algebraic multigrid algorithm [81] and the Gauss-Seidel smoother [82]. The solution procedure for the energy equations, phase continuity equations, turbulence model, and species transfer used the bi-conjugate gradient method, which was preconditioned by DILU [83].

Results and Discussion
This section presents the performed uncertainty estimation analysis (both grid and time step sensitivity) of the developed solver. Furthermore, the solver was tested for two ECN Spray A conditions (evaporating and non-evaporating) and compared to available experimental measurements [48,50].
The ECN Spray A test case included a fuel injector with a nominal nozzle outlet diameter of 0.09 mm, and n-dodecane as the injected fuel. In this work, the injection process was modeled using the blob injection model [62]-i.e., large droplets with a diameter similar to the nozzle size were injected into the domain. The velocity of the blobs was estimated from the fuel injection rate. The selected fuel injection curve [84] is given in Figure 2. The employed two-dimensional wedge domain is presented in Figure 3. The inlet is located on the left side and has the same size as the selected nozzle. The remaining surfaces (except the axisymmetric wedge planes) were treated as open boundaries. Furthermore, the inlet side of the geometry is scaled down, compared to the outlet side, to increase the grid density in the radial direction in the near nozzle region. The domain size is selected to be sufficiently large to reduce the influence of boundary conditions on the numerical results.

Uncertainty Analysis
The employed uncertainty analysis followed the procedures for unsteady flows described in [85] and was performed using the ReFRESCO application [86]. To quantify the uncertainty of the implemented solver, multiple simulations were carried out for the same flow conditions by systematically varying the density of the employed computational grid and the selected time step.
The initial grid was defined to have two cells per nozzle radius in the radial direction. With increasing distance from the nozzle location, the refinement was progressively decreased, which resulted in 26 radial cells in total. In the streamwise direction, the first cell (near the nozzle exit) had a width of 0.25 mm, and it was gradually increased to 2 mm over a length of 0.1 m (total length of the domain), resulting in 118 cells. The initial grid is presented in Figure 3. The remaining three coarser grids were obtained by uniformly decreasing the number of cells in the streamwise direction, but the number and distribution of radial cells remained the same. However, references [32,87] suggest that even finer grids are required to resolve the near-nozzle and mixing layer accurately.
The uncertainty analysis was performed for the liquid penetration length. The penetration length was calculated as the streamwise distance from the nozzle outlet, where the cumulative liquid fuel mass reached 98.5% of the total liquid mass located in the computational domain. The uncertainty analysis employed the evaporation conditions for the Spray A test case, which will be described in a more detailed manner in Section 4.2.2. The input data for the analysis are presented in Table 1. Table 1. Input data for the uncertainty estimation analysis. The values give the liquid spray penetration length (in millimeters) at 0.2 ms (after start of injection) for various grids and time steps.  Table 2 presents the output of the analysis following the notation given by [85]. φ 1 gives the input value of the most refined level, φ 0 is the extrapolated value, and U φ denotes the uncertainty estimate. The p and q variables present the achieved accuracy in space and time. Considering the previously described numerical settings in Section 3, the achieved second-order accuracy in space and the first-order accuracy in time were expected. The relatively high uncertainty value suggests that the liquid penetration length (for high threshold values) is sensitive to the employed spatial and temporal resolution. However, similar behavior was observed with various codes and models presented during the ECN Workshop [88].

Validation
The implemented model was tested and validated for evaporating and non-evaporating Spray A conditions. Both cases employed the finest computational grid and smallest time step from the previous section.

Non-Evaporating Conditions
This section demonstrates the predictive capability of the implemented breakup model by comparing the resulting droplet population in the dense part of the spray at non-evaporating conditions. The liquid fuel (n-dodecane) was injected into a vessel at 2.0 MPa and 300 K [50]. The resulting Sauter mean diameter (SMD) curve is compared to the available experimental measurements [50]. The comparison is shown in Figure 4. The given SMD curve was calculated as in [67]:  The SMD data were extracted using an axial sampling line after the SMD profile reached a steady state. The numerical results are in good agreement with the experimentally measured SMD curve, especially in the near-nozzle region where the SMD curve undergoes a rapid decline. However, the numerical results failed to predict the slight increase in the more stable part of the spray, but similar behavior was also reproduced with the Lagrangian solvers [89]. Overall, the implemented breakup model successfully predicted the rapid decline of the SMD curve and the stable droplet size in the farther part of the spray.
The obtained penetration plot for the liquid phase is shown in Figure 5. Following the approach presented in [89], the non-evaporating liquid penetration profile is compared to the experimental measurements given by [49,88]. As the comparison in Figure 5 indicates, the numerical results are in good agreement with the experimental measurements. The visual representation of the liquid spray (at t = 1.4 ms) is given in Figure 6.

Evaporating Conditions
In this test case, the liquid fuel at 363 K was injected into a vessel at 6.0 MPa and 900 K with 0% oxygen (non-reacting condition) [48]. These conditions, coupled with the previously described breakup regime, caused intense evaporation of the liquid phase. Therefore, the experimental measurements predicted the liquid penetration around 10 mm, which meant that practically all injected liquid mass evaporated in the first 10 mm in the axial direction.
The obtained penetration plot for the liquid phase and fuel vapor is shown in Figure 7. Here, the fuel vapor penetration is defined as the maximum streamwise distance from the nozzle outlet, where the mass fraction of the fuel vapor is 0.1%. As seen in Figure 7, the numerical results are in very good agreement with experimental measurements for both penetration curves. The model could successfully predict the fuel vapor penetration during the whole duration of fuel injection. Furthermore, the implemented model effectively predicted stable liquid penetration around 10 mm, but it gave a slight overprediction during the initial stabilization period.
The visual representation of the vapor penetration and liquid spray (at t = 1.4 ms) is given in Figure 8. The comparison of mixture fraction profiles at two different radial positions is given by Figures 9 and 10. Figure 9 gives the comparison at the streamwise location z = 25 mm in the radial direction. The "Gaussian" mixture fraction profile was predicted well, especially the peak value, but there was a slight underprediction for the outer part of the jet.  Figure 10 gives the comparison at the location z = 45 mm. Again, the shape of the "Gaussian" mixture fraction profile is in good agreement with the measurements, but there was a slight underprediction over the entire sampling radius. Overall, the implemented model successfully predicted the dynamic behavior of dense liquid sprays for the selected operating conditions.
The presented results were obtained on a desktop workstation with the AMD EPYC 7302 CPU. The validation case (finest computational grid and smallest time step) took approximately 4.5 h using five CPU cores. Since the solver is still in the development phase, there is a lot of potential for enhancing the performance of the implemented model.

Conclusions
The development of an advanced simulation framework for the dynamic behavior of dense evaporating liquid fuel sprays was presented. This work gave a detailed description of upgrades performed on the previously developed and published model, including updating the solver to a compressible formulation, implementing energy equations for all phases, and adding species transfer and evaporation functionality. Special attention was given to modeling the thermal phenomena occurring inside the moving droplets. The presented Eulerian multi-fluid model was implemented within the open-source foamextend library, a community-driven fork of OpenFOAM.
The validation section showed that the implemented solver could accurately predict the atomization process and secondary breakup in the dense region of the spray. Furthermore, it also correctly predicted the evaporation and mixing phenomena, which are more pronounced in the dilute part of the spray. The validation was performed for evaporating and non-evaporating ECN Spray conditions. The verification analysis proved that the solver behaves consistently, provided our numerical settings and computational grids. Therefore, the developed solver represents a stable foundation for further development. The solver is planned to be upgraded with multi-component functionality in future work, which should enable simulations of even more complex fuels.

Conflicts of Interest:
The authors declare no conflict of interest.