Performance Assessment of Chemical Kinetics Neural Ordinary Differential Equations in Pairwise Mixing Stirred Reactor

Turbulent combustion lies at the core of numerous technologies for power generation, transportation, propulsion, heating, and other industrial processes. The simulations of turbulent combustion are inherently complex owing to the intricacies of turbulent flows, chemical kinetics, and the interactions between them. The brute-force approach of direct numerical simulation (DNS) involves resolving all spatiotemporal scales in a flow field, which is highly precise but exorbitantly expensive for real-world problems such as internal combustion engines and gas turbine combustors. Alternatively, the averaging or filtering approaches in Reynolds-averaged Navier–Stokes (RANS) and large eddy simulation (LES) offer tractable ways to carry out turbulent combustion simulations. However, averaging or filtering leads to closure problems, particularly for the highly nonlinear chemical reaction source term. In LES, for instance, the computational expense is mitigated by filtering the governing equations so that only large scales of the flow field are resolved while motions at the small scales, generally referred to as subgrid scales, are modeled. However, modeling reactive scalars at the subgrid scales is challenging since the fluctuations at the unresolved small scales, predominantly caused by molecular mixing and reactions, strongly influence the filtered chemical reaction source term and vice versa. An effective approach for this purpose is the filtered density function (FDF) methodology which offers an exact way to formulate the reaction source term in a closed form, eliminating the need for its further modeling [1]. The FDF is essentially the LES equivalent of probability density function (PDF) method [2,3] in RANS; we thus generally refer to both as PDF henceforward.

The joint PDF of reactive scalars carries information about species concentrations and temperature evolution in a highly dimensional state space which makes obtaining the direct numerical solution of its transport equation extremely challenging. Hence, the method of Lagrangian Monte Carlo is often used to solve for the joint PDF of reactive scalars [2–4]. In this approach, the PDF transport equation is cast into an equivalent system of stochastic differential equations (SDEs) which is solved by the Monte Carlo method to determine the composition evolution of Lagrangian particles. The evolution of a large ensemble of such notional particles essentially provides an equivalent solution to the PDF transport equation for the scalars. A common approach that facilitates the solution of SDEs for each particle involves the use of a simple splitting method to decouple the mixing process (due to convection and diffusion) from the chemical reaction effects and to solve the corresponding equations in two separate steps [5–7]. For typical flame calculations, the chemical reaction step requires an enormous number of integrations of the chemical kinetics stiff ordinary differential equations (ODEs) for all particles. Therefore, developing modeling techniques that can effectively reduce stiff chemical kinetics integration cost leads to significant savings in computations, which is the main motivation of the present work.

A prevailing strategy to moderate the computational cost of chemistry calculations is through the reduction of size and stiffness of the detailed chemical mechanism. Most widely employed mechanism reduction methodologies include sensitivity analysis, directed relation graph (DRG) [8], principal component analysis (PCA) [9], computational singular perturbation (CSP) [10], quasi-steady-state approximations (QSSA) [11], partial-equilibrium approximations (PEA) [12], rate-controlled constrained-equilibrium (RCCE) [13–18], and intrinsic low dimensional manifold (ILDM) [19]. Although these approaches assist in alleviating computational demand, the resulting reduced mechanisms often retain some degree of numerical stiffness, necessitating the use of stiff ODE solvers to obtain the solution. Consequently, in some of the previous works, the reduced mechanisms are coupled with storage and retrieval methods such as look-up tables (LUT) [20] and in situ adaptive tabulation (ISAT) [7] to further decrease the cost of chemical kinetics integration. Despite the acceleration in computation time through tabulation methods, their application to PDF simulations is hindered by the necessity to perform direct integration of stiff chemistry at the initial stages of the simulation, which can make the tabulation prohibitively demanding for practical applications. Additionally, as more regions of the composition state space are accessed, the storage requirements and average retrieval times from tables can grow considerably.

In addition to the aforementioned classical modeling techniques, numerous researchers have focused on data-driven approaches based on diverse artificial neural network (ANN) topologies to address the high computational cost of chemistry integration and the increased storage requirement of tabulation methods. The pioneer works of Refs. [21,22] utilized a multilayer perceptron (MLP) topology neural network to tabulate relatively smaller finite-rate chemical kinetics and established its applicability in the simulation of turbulent flames. As the number of species increases, the composition space expands exponentially. To address the increased complexity in such cases, studies [23,24] advanced to clustering algorithms such as self-organizing-map (SOM) [25] and K-means [26] to cluster the composition space into several subdomains in order to manage complex chemical mechanisms. These studies reported good agreement with direct integration calculations, as well as significant savings in computational time and memory requirements. They did, however, identify a significant difficulty in the generation of representative datasets, which made the generalization of ANNs to actual conditions countered in real flames difficult. Later Refs. [27,28] encountered this issue by generating training data through an abstract problem (such as laminar flamelets) spanning the expected composition space and applying the trained ANN to simulate real turbulent flames.

Most prior research studies appear to have primarily focused on using an MLP neural network as a single-step chemistry integrator, with only minor differences in their network architectures. The species concentration and enthalpy at a given time are inputs to the ANN model, and the outputs are the species concentration/chemical source terms at the next time-step. Such ANN models have limited applicability in multistep predictions because they are trained by minimizing single-step prediction errors. The model could accumulate larger errors in longer-time predictions by recursively feeding the outputs from the previous time-step to the ANN as new inputs. This problem was addressed in Ref. [29], which adopted the neural network architecture of the residual network model (ResNet) to minimize multistep prediction error. Unlike MLP ANNs, ResNets are specifically designed for modeling time-series data. Recently, a class of machine learning algorithms known as neural ordinary differential equations (NODE) has emerged as superior alternatives to ResNets [30]. The NODE network is better suited for chemical kinetics integration because they have simpler architecture and better multistep prediction accuracy than simple MLPs and ResNets. Additionally, the NODE network can reduce the stiffness of the system, which can enable the use of an explicit ODE solver resulting in faster computations. The work of Ref. [31] was among the first to demonstrate the use of the NODE network in combustion chemistry. Because of the nonlinearity associated with the dynamics of various species, this study concluded that attempting to learn the source terms for all species concurrently was too unstable. As a result, separate networks were trained for each major species and temperature, with minor species being neglected. Subsequently, a methodical strategy was introduced in Ref. [32] to regularize and construct a single NODE network to learn the dynamics of coupled chemical kinetics. On a similar line, the work of Ref. [33] employed a stiffness-reduced autoencoder and NODE to develop a reduced-order model for methane. The study was specifically focused on training and demonstrating the model performance in a plug flow reactor and a continuous stirred tank reactor.

In the present study, we demonstrate the performance of a novel approach for coupling chemical kinetics NODE network with turbulent mixing to predict reacting systems. The NODE network is solely trained on canonical constant pressure homogeneous reactor data. This is a critical issue to address in order to ensure that the NODE network provides accurate predictions when used for complex turbulent combustion simulations. To our knowledge, no particular study has investigated the performance of the NODE network for chemical kinetics in the presence of the turbulent mixing. Investigations are carried out using a reacting pairwise mixing stirred reactor (PMSR) as the numerical solution procedure in PMSR is analogous to that in PDF methods [5–7,17,34]. The direct integration of detailed kinetic model (DKM) in the chemistry fractional step is replaced by NODE network integration and the computational performance of the NODE network is assessed by comparing the results with the former. Additional tests are performed to determine the sensitivity of the NODE network to numerical parameters such as error tolerances, ODE solvers, and the global time-step size. Although this methodology is directly relevant to simulations using the PDF method, it is also equally applicable to other approaches requiring direct integration of chemical mechanisms while employing splitting schemes such as DNS, unsteady flamelets, conditional moment closure (CMC), eddy dissipation concept (EDC), and partially stirred reactor (PaSR) models [35,36].

The paper is organized as follows: Sec. 2 outlines the basics of PMSR followed by a description of the NODE network-based chemistry integrator framework and NODE network adaptation in PMSR solution methodology. Section 3 discusses the performance and numerical efficiency of the NODE network for H₂-air combustion in PMSR for various configurations and mixing intensities. That includes investigation results pertaining to evaluating the applicability of the NODE network using various ODE solvers, error tolerance, and global time-steps. Finally, Sec. 4 provides the concluding remarks along with an overview of the capabilities of the proposed methodology.

Figure 1 encapsulates the methodology of the framework for the NODE network-based chemistry integrator. The procedure depicted briefly describes (1) the generation of training data by performing simulations of zero-dimensional constant pressure homogeneous reactor, (2) the preprocessing and regularization of training data to facilitate stable training of the NODE network, (3) the NODE network architecture and training methodology, and (4) the workflow to employ optimized NODE network to obtain the solution of Eq. (3). The detailed description is as follows:

1. We select the canonical constant pressure homogeneous reactor to generate the training data as the reaction equation (Eq. (3)) to be solved describes the evolution of species in the same reactor. During training, enthalpy is replaced by temperature as the latter has a nonzero source term and its redundancy offers an accuracy check during testing (by comparing it with the temperature obtained using the mixture composition). The training samples are generated using 30 constant pressure reactor simulations with initial temperatures ranging from 950 K to 1200 K and equivalence ratio in the range of 0.5–1.5 at atmospheric pressure (101,325 Pa) in order to generate representative composition space for real flames. The training of the NODE network depends solely on data and is not sensitive to the choice of ODE solvers used for data generation. Each of the 30 cases has its solution advanced from t = 0 to 10⁻² s with a fixed dt = 10⁻⁶ s, resulting in 10⁴ uniformly spaced points along the temporal composition trajectory for each simulation. As a result, a 3 × 10⁵ sample training dataset is obtained. In this study, we utilize the ODE solver LSODA via scipy toolkit [39] for direct integration of Eq. (3) with DKM description of chemistry. LSODA is a hybrid solver, switching automatically between the nonstiff Adams method and the stiff backward differentiation formula (BDF) solver.

2. The use of a full dataset (3 × 10⁵ data points) for training is unnecessary because most composition states reach equilibrium after ignition. Hence, we strategically sample the subset from the complete solution and regularize it in order to accelerate the NODE network training process. The samples of uneven log-spaced 50 points are chosen of 10⁴ time-steps for each simulation. A half of sampled points are placed before the ignition delay time (τ_ign) and the rest after the ignition. An example of sampling is shown in Fig. 1 for three cases corresponding to constant pressure hydrogen-air reactor at initial temperatures of T₀ = 1000 K, 1100 K, 1200 K, and stoichiometric equivalence ratio. The ignition delay time is defined as the instant when the reactor temperature reaches T₀ + 400 K. The chosen samples seem to capture the system dynamics of ignition curves well and enable training of the NODE network model with manageable computational requirements. The wide range of species concentration magnitudes is another challenge with training datasets for chemical reaction applications; for example, the major species mass fraction values are of the order 10⁻³ − 10⁰, while the minor species mass fractions are of the order of 10⁻³ − 10⁻⁸. Therefore, for network stability, all the composition variables are normalized based on the respective minimum (Φ_min) and maximum (Φ_max) values obtained from the complete training dataset

   $Φ¯=Φ−ΦminΦmax−Φmin$

   (5)
3. Equations (3) and (5) are encapsulated in the neural ODE network and mathematically represented as

   $dΦ¯(t)dt=Net(Φ¯(t),θ,t)$

   (6)

   where the NODE network, $Net(Φ¯(t),θ,t)$⁠, is the representation of the chemical reaction source term of $Φ¯(t)$ by the neural network with network parameters θ. The network consists of three hidden layers constituting 150 neurons each. Each hidden layer is followed by a nonlinear ELU activation function. The explicit Runge–Kutta method of order (4)5 (dopri5) [40] with absolute and relative error tolerances 10⁻⁹ and 10⁻⁷ is chosen, respectively, for integration during the training. Embedding of an explicit solver ensures that the network learns the nonstiff representation of the chemical reaction source term through training. To improve learning performance, the network is trained by sending data in batches. One training batch consists of 50 filtered time instances on the composition space evolution trajectory for 30 initial conditions. As a result, the network learns from all training points at the same time in each epoch. Each epoch begins with a forward pass to calculate the mean absolute error (MAE) of a training batch, followed by backpropagation to obtain loss function derivatives with respect to input parameters. Then, the neural network weights are adjusted to minimize the MAE loss function using the adaptive moment estimation (Adam) optimization technique for gradient descent. The learning rate is automatically varied from 10⁻² to 10⁻⁴ until the MAE is reduced to less than $5%$⁠. We adopt torchdiffeq algorithmic package implemented in PyTorch [30,40] to train the NODE network for chemical kinetics. Further details about NODE network formulation and training details can be found in Ref. [32].
4. Once the network is trained with the desired accuracy, it can be integrated with any ODE solver to obtain its state at any arbitrary time t_N > t₀. For example, with the normalized initial values $Φ¯0$ at time t₀, integrating forward to time t_N using an ODE solver gives

   $Φ¯N=Φ¯0+∫t0tNNet(Φ¯,θ,t)dt$

   (7)

   Subsequently, the normalized predicted values $Φ¯N$ can be converted to absolute composition vector Φ_N using Eq. (5).

At first, we consider a single PMSR test case for the demonstration. The test case has three incoming streams: air (⁠$79%$ N₂, 21% O₂) at 1200 K; H₂ at 500 K; and a pilot stream consisting of an equilibrium mixture corresponding to the inflow of the fuel/air conditions. The mass flowrates of these streams are in the ratio 0.873 : 0.027 : 0.1, respectively. Initially (t = 0), all particles are set to a mixture composition of inlet fuel and oxidizer streams resulting in an average temperature of 1004 K. Pressure is atmospheric throughout. Stochastic simulations with 1000 Monte Carlo particles are conducted over a span of five residence times to reach the statistically stationary state. DKM integration is performed using a hybrid LSODA solver [39] with an absolute error tolerance of η_a = 10⁻⁸ for species mass fraction and 10⁻⁶ for temperature, and relative error tolerance of η_r = 10⁻⁶. Tolerances are chosen based on the solver’s ability to produce converged solutions for a given set of operating conditions. The larger temperature η_a is because temperature values are of 2–3 orders of magnitude larger than those of species mass fractions. The numerical specifications in PMSR calculations are summarized in Table 1. The NODE network is integrated with an explicit ODE solver based on the Runge–Kutta method of order 5(4) (RK45) [39] with optimal error tolerances obtained with the procedure described in the previous study, i.e., η_a = 10⁻⁵ and η_r = 10⁻⁴ [32].

Figure 2(a) compares the evolution of ensemble average temperature and species mass fraction evolution in PMSR obtained with the NODE network and DKM on a log-scale plot. The NODE network predicts the average temperature in the reactor and ignition delay time accurately compared to DKM. The evolution of major species: fuel (H₂), oxidizer (O₂), and product (H₂O), along with the highly reactive OH radical is also well captured by the NODE network. The NODE network agreeably captures the monotonic increase in temperature and product, depletion of oxidizer and fuel, and nonmonotonic behavior of OH along with the fluctuations caused by stochastic particles. The similarities in instantaneous particle temperature predicted by both methods are also evident from Figs. 2(b) and 2(c). There are visible inflowing particles carrying fuel and oxidizer at respective inflow temperatures throughout the simulation. At initial times, the particles are clustered around specified initial temperatures of 1004 K. As these particles mix with higher temperature particles, they ignite and eventually reach a statistically stationary state close to the equilibrium state.

To further demonstrate the applicability of the NODE network in a broader composition space, results are presented for nonpremixed and premixed inlet stream configurations. In a nonpremixed configuration, the fuel and air enter the reactor as separate streams similar to the aforementioned test case. Whereas, in premixed configuration, only one inlet stream carrying the fuel and air mixture is used. For the premixed inlet, the equivalence ratio of the inflowing mixture is varied from lean to rich conditions at T₀ = 1000 K, and the dynamics of temperature evolution are captured as shown in Figs. 3(a) and 3(b). The NODE network shows an excellent comparison in 〈T〉 and root mean square (RMS) of T, $⟨T′2⟩$⁠, for all three equivalence ratios. The global errors in all 〈T〉 and other species mass fractions predominantly stay below $2%$ and $10%$⁠, respectively, as depicted in Fig. 3(c). With a premixed inlet, the majority of particles in the reactor act as an independent constant pressure homogeneous reactor. Hence, 〈T〉 profiles look similar to temperature evolution in the constant pressure homogeneous reactor. Also, the trends in ignition delay (τ_ign) and equilibrium temperatures (T_eq) are consistent with the homogeneous reactor [32]. For example, changing the equivalence ratio of a stoichiometric mixture to a lean or rich mixture results in longer chemical time scales characterized by τ_ign. The stoichiometric mixture has the highest T_eq and fastest ignition delay and the lean mixture has the lowest T_eq. Analyzing the statistics, it is observed that the NODE network exhibits smoother profiles compared to DKM. We speculate that this is due to the removal of small timescales during the stiffness reduction of the underlying chemical kinetics by the NODE network.

Figure 4(a) shows the 〈T〉 evolution in PMSR for lean to rich mixture fractions of inlet streams in nonpremixed configuration. For H₂-air the stoichiometric mixture fraction is Z_st ≈ 0.03. The inlet temperatures also vary according to mixture fraction values, resulting in further expansion of the composition state space occupied by particles. The larger mixture fraction contributes to increased fuel in the reactor, resulting in higher equilibrium temperatures. However, unlike the premixed case, the leaner case ignites faster indicating the importance of temperature in the ignition process due to the nonlinear dependence of reaction rates on temperature. The aggregate effect is faster ignition with increasing average inflow temperature. These trends are clearly evident in Figs. 4(a) and 4(b). Also, the mean and RMS of T compare well with the DKM. As depicted in Fig. 4(c), the largest errors are observed for H radical for fuel-rich conditions. However, the global errors in 〈T〉 and other species mass fractions predominantly stay below $4%$ and $10%$⁠, respectively.

Overall, the performance of the NODE network in Figs. 3 and 4 demonstrates the generalizability of the NODE network for a wide composition space and varying chemical timescales.

In realistic turbulent reacting flows, there is a strong interaction between chemical kinetics and mixing at various scales. Considering that the NODE network is solely trained on constant pressure homogeneous reactor with no mixing, it is important to establish the accuracy and performance of the NODE network when interacting with varying mixing timescales. In PMSR, micro-scale mixing is controlled by mixing timescale. On the other hand, mixing processes occurring at larger scales are characterized by the residence and pairing timescales. In this section, we present the results of the NODE network for varying residence, pairing, and mixing timescales, respectively. Apart from the specific timescales, the remaining numerical specifications of DKM and NODE network simulations presented in this section are the same as those summarized in Table 1.

The effect of residence timescale in PMSR is depicted in Fig. 5. Timescale τ_r is associated with large-scale mixing in the reactor. As observed from Figs. 5(a) and 5(b), NODE network predictions of 〈T〉 and $⟨T′2⟩$ are in good agreement with DKM for all the residence times. Higher residence times correspond to a lower influx of fresh particles and thus, fewer hot particles to ignite the reactor. As a result, ignition delay times are increased with increasing residence times causing a longer time required to reach an equilibrium temperature. It is also clear that the average temperature at a stationary state decreases with decreasing τ_r. This can be attributed to the more frequent replacement of ignited particles with colder inflowing particles which in turn reduces the average temperature of the system. In extreme cases of very small τ_r values, inflow/outflow processes dominate and more particles get removed from the reactor before ignition, causing the average temperature of the reactor to fluctuate around the temperatures of incoming particles (≈1083 K). Figure 5(c) shows the global errors for all species mass fractions and temperatures corresponding to three residence times. In the majority of the cases, global error in averaged species mass fractions remains below $11%$ including highly reactive radicals such as O, H₂O₂, and HO₂. However, the errors in minor species OH and H, two short-lived highly reactive radicals, rise to $25%$ for high τ_r. One possible explanation is that these radicals are known to be good QSSA species candidates which have shorter timescales and remain low in concentration [10]. In the process of removing stiffness, the NODE network loses some information about the fast species. Despite these local errors, the NODE network quite satisfactorily predicts the global ignition characteristics such as equilibrium temperature and ignition delay times.

The pairing timescale controls macro-mixing in PMSR by exchanging the particle partners at random. The smaller pairing times entail the frequent exchange of partners, hence more effective mixing of hot and fresh particles. This results in higher T_eq. This tendency can be observed from Figs. 6(a) and 6(b) with exact similarity. Also, the pairing time has almost no effect on the initial rate of increase in temperature or ignition delays. The global error in scalars remains below $10%$ except in the minor species radicals OH and H similar to the effects of τ_r. However, the 〈T〉 and $⟨T′2⟩$ predicted by the NODE network match well with those from DKM for all three pairing timescales. We also observe a peculiar behavior of sudden drop in T_eq for τ_p = 5 × 10⁻⁴, 10⁻³ s after ignition. This can be attributed to the local extinction of particles because of an inefficient supply of reactants caused by slow macro-mixing. Nevertheless, the NODE network is able to accurately capture this phenomenon as well.

At last, the interaction between micro-mixing and chemistry is studied by varying the mixing timescale. Figures 7(a) and 7(b) show the comparison between the NODE network and DKM predictions of 〈T〉 and $⟨T′2⟩$ for three τ_m values. For small mixing timescales such as τ_m = 10⁻⁴ s, mixing is fast and T_eq reaches the corresponding chemical equilibrium state value, T_eq = 2550 K. However, at larger mixing timescales the equilibrium temperature decreases due to slower reaction progress caused by imperfect micro-mixing. It is also observed that as mixing only tends to redistribute energy among particles at the initial times, τ_m has almost no influence on the initial rate of increase in temperature or ignition delays. Similar trends are noticed for even smaller τ_m values (not shown). The NODE network shows good agreement between the instantaneous mean and RMS temperature values for all the mixing timescales. As shown in Fig. 7(c), the global error in scalars predominantly remains below $10%$ except for OH radical.

Similar results are obtained when experiments are conducted using various combinations of residence, pairing, and mixing timescales between 10⁻⁴ and 10⁻² s. This study indicates that the NODE network accurately captures the combustion characteristics such as ignition, extinction, and statistics of thermochemical scalars. Overall, the NODE network can satisfactorily capture the interaction of chemistry and mixing over a wide range of timescales. The higher than average errors in H and OH radicals could be a result of system stiffness reduction by the NODE network leading to loss of information for very small timescales. However, further investigations are needed for conclusive elucidation.

The computational performance of the NODE network or DKM depends on numerous factors involved in obtaining the solution to Eq. (1) including splitting stepsize, ODE solver type, and error tolerances. The relative (η_r) and absolute (η_a) error tolerances have an impact on how well the adaptive ODE solvers with stepsize control such as LSODA and RK45 perform. The former controls the error in the solution relative to dependent variable values. The latter provides a threshold below which the error is disregarded. This threshold determines the accuracy when the solution approaches zero. In this section, we conduct two studies: first to investigate the performance of the NODE network for varying global time-step or the splitting time-step size, and second to explore the performance scalability of the NODE network with the increased number of chemistry integrations by increasing the number of particles.

Figure 8 shows that irrespective of error tolerance values and ODE solver type, the NODE network is consistently faster than DKM. It requires a fewer number of function evaluations than DKM, indicating a reduction in the stiffness of the system. For DKM simulations with a highly accurate BDF solver the relative error $εr,τign$ tends to increase with the increase in dt, while $εr,Teq$ first increases and then plateaus at $1%$⁠. However, for hybrid LSODA solver, $εr,τign$ increases with dt and $εr,Teq$ has a minimum at $dt=12τm$⁠. The NODE network exhibits similar error behavior for all the solvers and tolerances which establishes the applicability of splitting schemes for the NODE network. This indicates the robustness of the NODE network for various ODE solvers and error tolerances, and their ability to work with larger dt without sacrificing accuracy. The speedup factors of 4 to 6 are recorded with $dt=12τm$ for NODE network cases compared to the corresponding DKM case with BDF solver. Furthermore, even when comparing the DKM case with LSODA, the observed speedup factors range from 2 to 3.

To further assess the scalability of the NODE network with the number of ODE integrations, the number of particles in PMSR is increased from 100 to 10,000 and performance is recorded in Fig. 9 for the same set of ODE solvers and error tolerances as described in Fig. 8. For each simulation, the simulations are progressed in time for 200 time-steps, requiring the total number of ODE integrations of 2 × 10⁴ to 2 × 10⁶. Irrespective of error tolerance values and ODE solvers, the NODE network is consistently faster than the DKM owing to the reduced stiffness of the system. The speedup factors for the NODE network with respect to DKM simulations with BDF solvers remain approximately constant as the number of ODE integrations is increased. However, these factors tend to increase with reduced error tolerances. For example, for the RK45 solvers, the speedup factor of 5 is observed for (η_a, η_r) = (10⁻⁵, 10⁻⁴), and it is increased to a factor of 8 for (η_a, η_r) = (10⁻³, 10⁻³). For the lower order solver RK23, this speedup reaches a factor of 11. The global error in 〈T〉 is the largest for N = 100 which can be partly attributed to the statistical variations of the mean. The mean error reduces to less than $2%$ for N = 10,000, indicating that statistical convergence is achieved with an increase in the number of particles. This statistical convergence can be clearly observed for DKM, as well as the NODE network, in Fig. 9(c).

In summary, this investigation concludes that the NODE network is highly robust, efficient, and accurate in capturing the complex dynamics of chemistry and mixing. NODE network offers great flexibility in the selection of ODE solvers and error tolerances to yield an order of magnitude speedup compared to highly accurate DKM solutions using implicit ODE solvers. This indicates the great potential of NODE in accelerating computations in turbulent combustion applications.

In this study, we illustrate the methodology of coupling chemical kinetics, formulated with NODE, with mixing to simulate mixing/reacting systems in the context of PDF methods. A proof-of-concept is demonstrated in hydrogen-air combustion in the PMSR which corresponds to zero-dimensional PDF calculations. The hydrogen-air chemistry is described by a finite-rate mechanism involving 9 chemical species and 21 reaction steps. Results show that the NODE simulations accurately capture the combustion characteristics of PMSR for a wide range of composition spaces and mixing timescales. It is observed that in extreme conditions such as larger residence and mixing times the NODE introduces relatively larger errors in H and OH radicals which are speculated to be the effect of stiffness reduction. Nevertheless, these errors do not have a significant effect on global combustion characteristics including ignition delay and equilibrium temperatures. Subsequent investigations to evaluate the numerical efficiency of the NODE network affirm the robustness of NODE to larger splitting time-steps, ODE solver type, and error tolerances. Despite the chosen solver and tolerances, the NODE network is shown to be faster than respective DKM simulations. In the best case, with a selection of optimal solvers and tolerances, the NODE network achieves 11 times speedup for comparable accuracy to DKM simulations under specified operating conditions in this study. This analysis demonstrates the accuracy and computational efficiency of the NODE network in describing the dynamics of reacting systems in the presence of mixing. This warrants further development and application of the proposed methodology for large-scale turbulent combustion simulations.

This study is supported in part by the Office of the Vice President for Research at the University of Connecticut through the Research Excellence Program. The authors acknowledge the computing resources provided by the high-performance computing facilities at the University of Connecticut.

This article does not include research in which human participants were involved. Informed consent not applicable. This article does not include any research in which animal participants were involved.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.