Computational fluid dynamics (CFD) modeling of reacting flows, such as those encountered in gas turbine and internal combustion engines, are computationally demanding owing to the complex interactions among multiple physico-chemical phenomena and the need to resolve a wide range of spatiotemporal scales governing the evolution of a large number of reactive scalars (chemical species). In particular, modeling of detailed chemical kinetics presents a major bottleneck, which is governed by a stiff system of coupled ordinary differential equations (ODEs) and characterized by high condition number of the corresponding chemical Jacobian matrix (Shampine, 1993). In addition, as the complexity of fuel chemistry increases, the dimensionality and stiffness of the ODE system also become more prohibitive (Lu and Law, 2009). To address these computational challenges, kinetic model reduction is typically performed (Lu and Law, 2005; Jones and Rigopoulos, 2005; Valorani et al., 2006; Maas and Pope, 1992), but it often leads to less reliable description of chemical kinetics.

In this context, data-driven machine learning (ML) techniques (Dana, 1987; Wold et al., 1987) have been extensively explored in an effort to emulate chemical kinetics and accelerate detailed finite-rate chemistry computations. Some of these approaches reduce the dimensionality of the reaction system via linear projection onto an appropriate basis typically identified through Principal Component Analysis (PCA). Artificial neural network (ANNs) are, then, employed for regression of the PC source terms and transport coefficients (Owoyele et al., 2017; Kumar et al., 2023). On the other hand, ANN-based approaches have also been used for predicting chemical source terms directly from the thermochemical state (Christo et al., 1996; Blasco et al., 1998; Sen et al., 2010; Ranade et al., 2019; Wan et al., 2020). More recently, novel deep learning architectures have been explored to capture the temporal evolution of chemical kinetics and combustion simulations, such as Deep Operator Network (DeepONet) (Kumar and Echekki, 2024) and Fourier Neural Operator (FNO) (Zhang et al., 2024).

Despite the potential benefits of employing ML methods, When they are trained in an a-priori or offline setting and then coupled with a numerical solver, predicted solutions may diverge or become unstable. Due to the non-linearity of the combustion process, even minor predictive errors can escalate into significant discrepancies in the temporal evolution of thermochemical state due to error accumulation. To handle this issue, an alternative and more robust data-driven technique for chemical source term computations based on neural ordinary differential equations (NeuralODEs, NODEs) (Chen et al., 2018) was developed by the authors, known as ChemNODE (Owoyele and Pal, 2021). It combines the chemical source term predictions with ODE integration in an a-posteriori learning paradigm, where the source terms predicted by the neural network are passed to the ODE solver for time integration, and the neural network weights are optimized to minimize the loss computed between the predicted and ground truth thermochemical states (comprised of species mass fractions and thermodynamic variables). A key advantage of this approach is that NODEs learn continuous-time dynamics which can be integrated using existing ODE solvers. This ensures that the predicted thermochemical state, even after a long-time horizon, remains adherent to the ground truth solution trajectory. For relatively large chemical kinetic mechanisms, the coupling of a non-linear autoencoder (AE) to perform dimensionality reduction and a NODE to evolve the dynamics in the lower-dimensional latent space has shown promise (Vijayarangan et al., 2024). A similar approach has also been pursued in the field of astro-chemistry (Sun Tang and Turk, 2022; Maes et al., 2024).

Another major limitation of traditional black-box ML techniques applied to the modeling of chemical kinetics is that they do not inherently incorporate conservation laws, which can adversely impact simulation accuracy and hinder the integration of ML surrogate models with multidimensional CFD solvers. With the advent of scientific machine learning, the practice of combining domain science-specific constraints to embed physics into neural network training has emerged lately (Raissi et al., 2019). This is achieved through regularization of the loss function, in which physical laws are embedded in the learning process as soft constraints (Kumar Tadbhagya et al., 2023; Almeldein and Van Dam, 2023; Kumar et al., 2025; Weng et al., 2025; Kercher and Votsmeier, 2025). This has been shown to improve model accuracy and generalization, especially in scenarios where traditional methods struggle with noisy data and high-dimensional problems governed by parameterized differential equations. Another way to enforce conservation laws is by enforcing a hard constraint in the neural network architecture through constraint layers (Sturm and Wexler, 2022; Mohan et al., 2023).

In light of the above discussion, the overarching goal of the present work is to demonstrate an end-to-end physics-constrained AE-NODE framework called Phy-ChemNODE for data-driven modeling and temporal emulation of stiff chemical kinetics for hydrocarbon fuels. Unlike previous work (Vijayarangan et al., 2024; Sun Tang and Turk, 2022; Maes et al., 2024), the present study demonstrates that this integrated approach improves physical consistency of the resulting data-driven model with mass conservation laws and accelerates training convergence, while simultaneously reducing temporal stiffness of the chemical kinetic system. The remainder of the paper is organized as follows: the Phy-ChemNODE framework is first outlined along with details of the physics-constrained formulation and training methodology in Sections 2, 3, respectively. Subsequently, results from a-posteriori proof-of-concept studies are discussed in Section 4. Finally, the Conclusions section summarizes the major findings and directions for future work.

In combustion CFD simulations, it is a common numerical approach to decouple finite-rate chemistry from transport using operator splitting. Chemistry is solved (independently from advective and diffusive transport) within each computational grid cell considered as a homogeneous reactor, which is equivalent to solving a system of stiff ODEs. The temporal evolution of N s reactive scalars (chemical species) can be defined by: d Y k d t = ω k ̇ ρ , k = 1,2,3 , … , N s (1) where Y k is the mass fraction of k t h species ( N s being the total number of species), ω k ̇ is the corresponding chemical source term computed using law of mass action, and ρ refers to mixture density. The temporal evolution of temperature (T) is governed by a similar ODE defined as follows: d T d t = − ∑ k = 1 N s h k ω k ̇ ρ c p , k = 1,2,3 , … , N s (2) where h k and c p refer to the enthalpy of k t h species and mixture-averaged constant pressure specific heat, respectively. To calculate the source terms, one needs to account for several elementary reactions involving production and consumption of multiple species. As the chemical mechanism becomes larger, the number of chemical species and elementary reactions also increase (Lu and Law, 2009). This leads to prohibitive computational costs since all chemical time scales must be fully resolved. In the NODE-based data-driven framework (Owoyele and Pal, 2021), the expensive physics-based computation of chemical source terms is replaced by a neural network, which can be described as Equation 3 below: d Φ d t = f Φ , t ; Θ (3) where Φ = [ T , Y 1 , Y 2 , … , Y N s ] is the vector of thermochemical state (temperature and species mass fractions), and f ( Φ , t ; Θ ) is a feedforward neural network parameterized by weights Θ . For larger chemical mechanisms, Φ increases in dimensionality and stiffness. To address the high dimensionality, a non-linear AE is coupled with the NODE for dimensionality reduction, so that the NODE learns the temporal evolution of the dynamical system in a reduced-order latent space obtained from the non-linear projection of the AE. A schematic of the coupled AE-NODE data-driven modeling framework is shown in Figure 1. The model training process is posed as an optimization problem of determining the optimal parameters of the encoder ( φ ) , NODE ( h ( z ) ) , and the decoder ( ψ ) networks in an end-to-end manner, that minimizes the loss function defined as: L Phy−ChemNODE = λ rec L rec + L data + λ z L z + ∑ j = 1 N e l λ e l − j L e l − j (4) where the reconstruction loss L rec = L ( Φ , Φ ̃ ) measures the loss between ground truth ( Φ ) and corresponding encoder-decoder mapping ( Φ ̃ = ψ ( φ ( Φ ) ) , the data loss L data = L ( Φ , Φ ̂ ) measures the loss between ground truth and encoder + NODE + decoder prediction ( Φ ^ ) , and the latent loss L z = L ( z ̄ , z ) measures the loss between encoder mapping of ground truth ( z ̄ = φ ( Φ ) ) and encoder + NODE prediction. Each loss term is chosen to be in mean absolute error (MAE) form. It is noted that the loss terms L rec and L z ensure that the encoder and decoder mappings are bijective or unique. The loss function also contains elemental mass conservation constraints (Kumar Tadbhagya et al., 2023; Kumar et al., 2025), defined in Equation 5 below: L e l − j = 1 N ∑ i = 1 N ∑ k = 1 N s N j k W j | Y k , i − Y ̂ k , i | W k (5) where L e l − j refers to the loss associated with mass conservation of element j (in the chemical system with a total of N e l elements). Y ̂ k , i and Y k , i correspond to the AE + NODE predicted and ground truth mass fractions of k t h species, respectively. W j is the atomic mass of element j , N j k is the number of atoms of element j in k t h species, W k is the molecular weight of k t h species, and N is the number of training data points. Lastly, the weights λ rec , λ z , and λ e l − j in Equation 4 balance the contributions from the different loss terms.

For proof-of-concept demonstration of Phy-ChemNODE, an autoigniting methane-oxygen ( CH 4 - O 2 ) zero-dimensional (0D) homogeneous reactor is considered, at a constant pressure of 20 atm that corresponds to practical high-pressure isobaric operating conditions for stationary gas turbine engines. The detailed chemical mechanism (Gregory et al., 2016) consists of 32 species and 266 chemical reactions. The ground truth data for model training is generated using Cantera (Goodwin et al., 2009), which solves the coupled ODE system (Equations 1, 2) with detailed chemistry. The thermodynamic and composition space chosen for data generation comprises 9 equispaced initial temperatures in the range T i = [1600 K, 2000 K] and 11 equivalence ratios within ϕ = [ 1.0 , 1.5 ] resulting in a total of 99 initial conditions. Each of these initial conditions is integrated to chemical equilibrium. Ground truth data is generated by selecting time instants such that the change in any of the thermochemical scalars (species mass fractions or temperature) between two successive time instants is greater than 1% of their corresponding overall ranges of variation, and then downsampling 200 points from the selected time instants. A 70 % , 20 % , 10 % random split (based on initial conditions) is used to obtain the training, validation, and test datasets, respectively. The AE-NODE model is initialized with the same initial conditions (during training) as the physics-based simulations. The input to the encoder is the vector Φ containing the temperature and species mass fractions, which is scaled using the maximum and minimum of the training data, respectively, and the output is a vector in the latent space ( z ) . The decoder has the same dense architecture as the encoder, with an input size equal to the latent dimension ( d i m ( z ) = 4 ) (chosen based on a sensitivity study discussed in Section 4) and the output size equal to the physical space vector ( d i m ( Φ ) = 33 ) . Both the encoder and decoder have 5 hidden layers with 64 neurons each and Exponential Linear Unit (ELU) activation function. The NODE has the same input and output dimensions as the latent space ( d i m ( z ) = 4 ) , and a 4 hidden-layer dense network with 64 neurons in each hidden layer and ELU activation function, to model the chemical source terms in the lower-dimensional latent space. The output layers for the encoder, NODE, and the decoder are considered to be linear. The above discussed AE-NODE architecture was finalized based on multiple ablation studies, targeting a good balance between dimensionality reduction and overall model accuracy. Although more exhaustive hyperparameter optimization can be performed to further optimize the number of latent variables, it is out of scope of the current work.

The forward pass through the NODE requires time integration, for which a 4th order explicit Runge-Kutta (RK) solver is used. Once the time integration is completed, the temporal trajectories are mapped back to physical space and the loss is computed using Equation 4. The gradients for updating the neural network parameters are calculated using backward adjoint automatic differentiation and ADAM optimizer with exponential learning rate decay (every 200 epochs) is used. The model is trained for 10000 epochs. To ensure that all the loss terms are of similar magnitude, λ rec = 5.0 is used. Further, λ z = 0.05 and λ e l − H = λ e l − C = λ e l − O = 0.5 are chosen based on a hyperparameter sweep (discussed in Section 4). The training framework was implemented in Julia programming language using Flux.jl (Innes et al., 2018) library and the model was trained on 2 AMD EPYC 7713 64-core processors for a walltime of 96 h.

To determine the optimal weighting of terms in the loss function and the latent space dimensionality, hyperparameter studies were carried out. Figure 2a compares the loss terms computed on the validation set (post training) corresponding to different λ z values for λ rec = 5 , λ e l − j = 0.5 , and fixed size of the latent space (dim ( z ) = 4 ), and Figure 2b shows a similar comparison for varying size of the latent space (dim ( z ) ) with λ rec = 5 , λ z = 0.05 , and λ e l − j = 0.5 . Moreover, to assess the impact of adding elemental mass conservation constraints to the training loss function, Figure 2c compares the decay of data loss (during training) between the cases trained with ( λ ele−j = 0.5 ) and without elemental mass constraints ( λ ele−j = 0 ) on the validation set. It can be clearly seen that incorporating the soft constraints in the loss function results in lower loss for the same number of epochs, thereby enabling more efficient model training and faster training convergence.

The trained Phy-ChemNODE model is then deployed for predicting the temporal evolution of thermochemical scalars in a-posteriori autoregressive tests. Figure 3 plots the temporal evolution of temperature and a subset of species mass fractions ( CH 4 , CO, CO 2 , OH, and O 2 ) for an initial condition (corresponding to T 0 = 1600 K) in the training set ( ϕ = 1.0 ) and test set ( ϕ = 1.1 ) , where ground truth data is indicated by solid lines and the predicted Phy-ChemNODE solutions are shown in markers. Overall, great agreement can be observed between the predictions and ground truth data. Figure 4 shows the temporal evolution of a few intermediate species for another set of initial conditions in the training set ( T 0 = 1650 K, ϕ = 1.0 ) and the test set ( T 0 = 1700 K , ϕ = 1.05 ), again demonstrating high accuracy. As further quantification of the accuracy of Phy-ChemNODE framework, Figure 5 shows the test set MAEs for a subset of the thermochemical scalars including both major and minor species, scaled by their corresponding data ranges, indicating that Phy-ChemNODE performs extremely well in capturing the temporal dynamics.

Based on inference on an Intel i7-1165G7 workstation with 16 cores, Phy-ChemNODE yields speedups of 6x and 860x over the full chemical mechanism in terms of overall simulation walltime, when deployed with implicit (backward differentiation formula (BDF)) and explicit (RK45) solvers, respectively. Figure 6 plots the latent space temporal dynamics corresponding to the initial conditions of T 0 = 1650 K, ϕ = 1.0 . It is evident that the evolution of latent variables is much smoother than that of thermochemical scalars in the physical space earlier shown in Figure 4a. This indicates that the coupled AE-NODE approach significantly reduces the temporal stiffness of the chemical kinetic system. Similar stiffness reduction in the latent space was also observed by Nair et al. (2025) in case of AE-NODE models of advection-dominated dynamical systems. Lastly, Figure 7 shows the predicted temporal evolution of C, H, and O mass fractions from autoregressive tests corresponding to certain initial conditions from the training and test sets. Evidently, the model trained with elemental mass constraints in the loss function (Phy-ChemNODE/PCNODE) conserves the elemental mass fractions during deployment much better than the MAE-trained case without constraints.

In future studies, multiple avenues for further extension of the Phy-ChemNODE framework will be pursued. These include: (a) efficient scaling of the training workflow to wider ranges of thermodynamic conditions (including multiple pressures), larger kinetic mechanisms, and multicomponent fuels; (b) integration of Phy-ChemNODE with CFD solvers and demonstration of accelerated multidimensional reacting flow simulations on modern high-performance computing (HPC) platforms; (c) demonstration for constant-volume combustion; (d) incorporation of uncertainty quantification; and (e) exploration of training methodologies to enhance model out-of-distribution generalizability.

In this work, an end-to-end physics-constrained AE-NODE framework (PhyChemNODE) was introduced for accelerated temporal emulation of stiff chemical kinetics targeting complex hydrocarbon fuel combustion. The deep learning approach employed a non-linear AE for dimensionality reduction of the thermochemical state and utilized a NODE to learn temporal dynamics of the kinetic system in the low-dimensional latent space obtained from the AE. In addition, elemental mass conservation constraints were included in the loss function during training of the data-driven model to ensure that total mass and mass of each elemental species are conserved. Proof-of-concept studies were performed for homogeneous autoignition of a methane-oxygen mixture over a range of composition and thermodynamic conditions at high pressure. The results showed that the physics-based constraints not only improve physical consistency of the resulting data-driven model, but also enhance training efficiency. In addition, Phy-ChemNODE achieved 6–860 × speedup relative to the detailed chemical mechanism depending on the type of ODE solver (implicit or explicit) used for numerical integration during autoregressive inference tests. Temporal evolution of the latent variables was visualized and it was found that the coupled AE-NODE approach leads to reduced temporal stiffness of the chemical kinetic system. In future work, Phy-ChemNODE will be further scaled to larger kinetic mechanisms and wider ranges of thermodynamic conditions (including pressure), and demonstrated for multidimensional combustion CFD simulations. In addition, uncertainty quantification techniques will be incorporated and training methodologies will be explored to enhance out-of-distribution generalizability.