# A virtual chemistry framework compatible with thermodynamic and kinetic standards

Matthieu Préteseille, Tan-Phong Luu, Nasser Darabiha, Benoît Fiorina

► To cite this version:

Matthieu Préteseille, Tan-Phong Luu, Nasser Darabiha, Benoît Fiorina. A virtual chemistry framework compatible with thermodynamic and kinetic standards. Proceedings of the Combustion Institute, 2025, 41, pp.105894. 10.1016/j.proci.2025.105894. hal-05316906

HAL Id: hal-05316906

https://hal.science/hal-05316906v1

Submitted on 15 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A virtual chemistry framework compatible with thermodynamic and kinetic standards

<figure><img src="image_3.png" /></figure>

Matthieu Prétesteille<sup>1, ∗</sup>, Tan-Phong Luu<sup>1</sup>, Nasser Darabiha<sup>1</sup>, Benoît Fiorina<sup>1</sup>

Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, 91190, Gif-sur-Yvette, France

## ARTICLE INFO

**Keywords:**

Virtual chemistry

Kinetic mechanism

Optimization

Chemistry reduction

Evolutionary strategy

## ABSTRACT

Accurate combustion modeling is essential for decarbonization, yet simulating the combustion of fuels like sustainable aviation fuels (SAF) remains challenging due to their complex chemistry. Detailed kinetic mechanisms, often involving thousands of reactions, make industrial turbulent combustion simulations computationally prohibitive, necessitating kinetic mechanism reduction techniques. Among these, virtual chemistry offers a promising solution by creating highly-reduced schemes with virtual species and reactions, optimized to reproduce specific flame characteristics. However, its original formalism relies on tabulated stoichiometric coefficients and kinetic rates as functions of the local mixture fraction, requiring source code modifications in traditional reactive flow solvers and limiting its adoption in the combustion community.

This study revisits the virtual chemistry formalism to eliminate parameter tabulation, aiming to develop universal 4-reaction virtual mechanisms for any fuel, directly compatible with traditional software like CHEMKIN and CANTERA. The methodology, divided into two steps, optimizes thermodynamic and chemical properties. Beyond retrieving the mixture's enthalpy and heat capacity, this new thermodynamic optimization also targets the mixture's Gibbs energy, enabling the recovery of equilibrium states across various conditions. Consequently, these novel virtual chemical mechanisms introduce equilibrium reactions, compensating for the lost degrees of freedom due to the elimination of parameter tabulation and bridging the gap between virtual and classical chemistry. Reaction rates are then optimized to match temperature and heat release rate profiles from detailed chemistry.

The formalism's versatility is demonstrated through the optimization of 4-reaction hydrogen and SAF virtual schemes using reference equilibrium and perfectly-stirred reactor computations on CANTERA. By adhering to standard NASA polynomial approximations for thermodynamic properties and an Arrhenius-like closure for kinetic rate constants, integrating virtual chemistry into any reactive flow solver is now straightforward.

## 1. Introduction

In the pursuit of decarbonization, designing efficient net-zero combustion chambers necessitates the use of hydrogen or more complex fuels, such as sustainable aviation fuels. However, accurately modeling their combustion regimes at a reasonable computational cost is only feasible through reduced chemistry [1].

Various routes have been developed to simplify combustion chemistry, namely reduced [2], tabulated [3,4] and optimized global chemistry [5,6]. Among them, the virtual chemistry method [7] is efficient to capture, at low-CPU cost, flame properties that are sensible to complex chemistry phenomena, such as auto-ignition [8] or pollutant production like NOₓ [9] or soot [10]. The principle of the method consists in introducing virtual species and reactions whose thermodynamic

and chemical properties are optimized by optimization algorithms to retrieve properties of reference flames gathered in a learning database.

Like many other global chemistry methods [5], the current mathematical formulation of virtual chemistry requires the tabulation of rate constants as a function of local data such as the mixture fraction [7] to simulate stratified flames, or enthalpy [11] to calculate non-adiabatic combustion regimes. This formulation is highly efficient in terms of accuracy and computational cost. However, its form deviates from the thermodynamic and kinetic packages commonly used in CFD codes such as CHEMKIN [12] or CANTERA [13]. As a result, integrating this kinetic mechanism necessitates modifications to CFD solver source codes, adding complexity to its implementation.

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: matthieu.pretesteille@centralesupelec.fr (M. Prétesteille).</footnote>
The aim of this work is to redefine a new virtual chemistry formalism compatible with standard thermodynamic and kinetic formats whose implementation in CFD codes is straightforward. It requires a fundamental rethinking of the mathematical formalism underlying the virtual chemistry concept, specifically through the introduction of atoms and entropy.

Thermodynamic optimization is presented in Section 2. Thermodynamic properties of virtual species are adjusted to accurately reproduce detailed chemical schemes thermochemical states at equilibrium. A particular emphasis is placed on ensuring that the method satisfies the second law of thermodynamics, leading to standardized thermodynamic properties consistent with NASA's formulations. Results are discussed in Section 3.

Section 4 presents virtual scheme's architecture, whose structure and reaction rates are optimized. Standard Arrhenius laws are retained to model the kinetic rate constants of the set of virtual elementary reactions. The scheme is then trained to recover a set of Perfectly-Stirred Reactors (PSR) reference data, enabling its accurate capture of combustion dynamics. The methodology is then validated on Hydrogen-Air and nDodecane-Air reactive configurations in Section 5. Concluding remarks are finally given in Section 6.

## 2. Thermodynamics optimization

### 2.1. Reduced mixture thermochemical properties

A reduced ensemble of real and virtual species are considered to model fresh and burnt gases of a reactive mixture. Real species are usually the fuel F, the oxidizer Ox, the dilutant D, and, optionally, a selection of products whose concentration should be predicted, such as H₂O or CO₂. Noting $N_p^r$ the number of real combustion products $P^r$, the total number of real species $N^r$ is given by $N^r = 3 + N_p^r$. A virtual species belongs either to the $N_p^v$ virtual products $P^v$ or $N_I^v$ virtual intermediate species $I^v$. The total number of species in the considered gas mixture is therefore $N = N^r + N_p^v + N_I^v$. By definition, intermediates are absent in both fresh and burnt gases. Considering $N_e$ elements in the species, any real or virtual species k is defined by its element composition $E_i$, $i \in [1..N_e]$. Then, $N_I^k$ represents the number of lth element in the kth species.

For any mixture, the thermodynamical state of a species k is characterized by its heat capacity at constant pressure $C_{p,k}^o$, enthalpy $H_k^o$, and entropy $S_k^o$ at standard pressure, modeled by the polynomial forms of temperature T [14]:

$$ \frac{C_{p,k}^o(T)}{R} = \sum_{i=1}^{5} a_i^k T^{i-1} \qquad (1) $$

$$ \frac{H_k^o(T)}{RT} = \sum_{i=1}^{5} \frac{a_i^k}{i} T^{i-1} + \frac{a_6^k}{T} \qquad (2) $$

$$ \frac{S_k^o(T)}{R} = a_1^k \ln(T) + \sum_{i=2}^{5} \frac{a_i^k}{i-1} T^{i-1} + a_7^k. \qquad (3) $$

In these equations, R is the universal gas constant, and $(a_1^k, a_2^k, a_3^k, a_4^k, a_5^k, a_6^k, a_7^k)$ are the so-called NASA7 coefficients of the kth species, which can be real or virtual.

The mean mixture's heat capacity, $\overline{C}_p = \sum_{k=1}^N C_{p,k}^o X_k$, and enthalpy, $\overline{H} = \sum_{k=1}^N H_k^o X_k$, are expressed from Eq. (1), (2), (3) as:

$$ \overline{C}_p = R \sum_{i=1}^{5} \overline{a}_i T^{i-1} \qquad (4) $$

$$ \overline{H} = R \left( \sum_{i=1}^{5} \frac{\overline{a}_i}{i} T^i + \overline{a}_6 \right) \qquad (5) $$

with $\overline{a}_i = \sum_{k=1}^N a_i^k X_k$ symbolizing the ith mixture's NASA7 coefficient.

<figure><img src="image_4.png" /><figcaption>Fig. 1. Optimization protocol for thermodynamics.</figcaption></figure>

Let $\mathcal{R}$ and $\mathcal{V}$ denote the ensemble of real and virtual species that compose the reduced mixture. Decomposing $\overline{a}_i$ in their respective contributions yields:

$$ \overline{a}_i = \sum_{k \in \mathcal{R}} a_i^k X_k + \sum_{k \in \mathcal{V}} a_i^k X_k. \qquad (6) $$

It is important to note that the coefficients $a_i^k$ for real species ($k \in \mathcal{R}$) are known from reference data, while those for virtual species ($k \in \mathcal{V}$) will be optimized.

The mean entropy $\overline{S}$ is also expressed in terms of mixture's NASA7 coefficients as follows:

$$ \overline{S} = R \left( \overline{a}_1 \ln(T) + \sum_{i=2}^{5} \frac{\overline{a}_i}{i-1} T^{i-1} + \overline{a}_7 \right) \qquad (7) $$

with $\overline{a}_7 = \overline{a}_7 - \sum_{k=1}^N [\ln(X_k) + \ln(\frac{P}{P_o})] X_k$.

Mixture Gibbs free energy $\overline{G} = \overline{H} - T\overline{S}$ is derived from Eq. (4), (5), (7) accordingly:

$$ \overline{G} = R \left( \underbrace{\overline{a}_1 T [1 - \ln(T)]}_{\theta_1(T)} + \sum_{i=2}^{5} \underbrace{\overline{a}_i T^i [\frac{1}{i} - \frac{T}{i-1}]}_{\theta_2(T) \text{ to } \theta_5(T)} + \underbrace{\overline{a}_6}_{\theta_6(T)} - \underbrace{\overline{a}_7 T}_{\theta_7(T)} \right). \qquad (8) $$

Finally, the mean molecular weight $\overline{M}$ is defined by $\overline{M} = \sum_{k=1}^N M_k X_k$.

### 2.2. Optimization of virtual products $P^v$

The optimization process of the virtual species properties is described in Fig. 1.

#### Reference detailed thermochemistry

Reference thermodynamic properties of the targeted burnt gases are determined by computing the equilibrium states of a reference detailed chemistry composed of $N_d$ species for a given range of $N_c$ target conditions (equivalence ratio, enthalpy, pressure...). These reference mean thermodynamic states are represented by the vector $\overline{\Phi}^d = (\overline{C}_p^d, \overline{H}^d, \overline{S}^d, \overline{M}^d)$.

#### Virtual species properties

Since intermediate species $I_i$ are absent in the burnt gases, their definition will be addressed later.

The reduced mixture contains $N_p$ combustion products, divided into $N_p^r$ real products (selected by the user) and $N_p^v$ virtual products:

$$ N_p = N_p^r + N_p^v. \qquad (9) $$

Given a chosen set of real products, the number of virtual products $N_p^v$ is determined so as to satisfy the following empirically defined constraints, ensuring the feasibility of the procedure:

* $N_p^v \ge 1$, to allow the optimization to proceed.
• The total number of products satisfies $N_p = N_e$, where $N_e$ denotes the number of elements. This constraint facilitates the determination of the stoichiometric coefficients of the elementary reactions, as discussed later in the kinetics section.

As shown in the results, the condition $N_p = N_e$ does not reduce the optimization quality for the configurations investigated in the present study. In order to capture broader thermodynamic equilibrium situations, the restriction can be lifted, so that $N_p$ may exceed $N_e$.

Atomic compositions and thermodynamic properties of the introduced virtual products $P^v$ are then optimized to reproduce, together with the real species, reference equilibrium states.

The optimization vector $x_{thermo}$ reads accordingly:

$$x_{\text{thermo}} = \left( \underbrace{(a_1^k, a_2^k, a_3^k, a_4^k, a_5^k, a_6^k, a_7^k)}_{\text{Thermodynamic coefficients}}, \underbrace{(N_C^k, N_O^k, N_H^k, \dots)}_{\text{Element composition}}, \dots \forall k \in [1..N_p^v] \right). \quad (10)$$

In contrast to conventional chemistry, where atomic compositions are restricted to integer values (∈ $\mathbb{N}$), each virtual product has here a composition expressed in positive real numbers (∈ $\mathbb{R}_+$). Therefore, classical continuous optimization algorithms such as evolutionary algorithms can be employed. To start the algorithm, a randomly-chosen population is used to initialize $x_{thermo}$.

### Gibbs energy minimization

For this set of parameters, equilibrium states of the reduced mixture are computed via Gibbs energy minimization, yielding species equilibrium molar fractions $X^{eq}$ as dictated by the first and second principles of thermodynamics.

### Mixture's thermochemical equilibrium states

The thermodynamic state of the reduced mixture composed of $N^r$ real and $N_p^v$ virtual species is represented by $\overline{\Phi}^v = (\overline{C}_p^v, \overline{H}_p^v, \overline{S}_p^v, \overline{M}_p^v)$. The thermochemical properties of the $N_p^v$ virtual species are optimized so that $\overline{\Phi}^v$ approximate the reference thermochemical equilibrium states $\overline{\Phi}^d$.

### Fitness function

The single-objective fitness function $f_{thermo}$ of the evolution algorithm is defined as:

$$f_{\text{thermo}} = f_{\text{thermo}}^{\overline{M}} + f_{\text{thermo}}^{\overline{a}} \quad (11)$$

$f_{thermo}$ aims to measure the differences between $\overline{\Phi}^d$ and $\overline{\Phi}^v$ at equilibrium.

The first fitness term, $f_{thermo}^{\overline{M}}$, recovers the reference equilibrium mean molecular weight [7]:

$$f_{\text{thermo}}^{\overline{M}} = \sum_{c=1}^{N_c} \left| \frac{\overline{M}^v(\psi_c) - \overline{M}^d(\psi_c)}{\overline{M}^d(\psi_c)} \right|. \quad (12)$$

In this equation, $N_c$ is the total number of targeted conditions and $\psi_c$ represents the cth target condition.

The second contribution in Eq. (11) aims to capture the temperature dependency of mixture's heat capacity, enthalpy and entropy. The principle is to find the best fit of $(\overline{a}_i)_{i \in [1..6]}$ and $\tilde{a}_7$ intervening in Eq. (4), (5), (7) reproducing the reference polynomial coefficients. Therefore, $f_{thermo}^{\overline{a}}$ reads:

$$f_{\text{thermo}}^{\overline{a}} = \sum_{c=1}^{N_c} \left( \sum_{i=1}^{6} \Theta_i^*(\psi_c) \left| \frac{\overline{a}_i^v - \overline{a}_i^d}{\overline{a}_i^d} \right| (\psi_c) + \Theta_7^*(\psi_c) \left| \frac{\overline{a}_7^v - \overline{a}_7^d}{\overline{a}_7^d} \right| (\psi_c) \right) \quad (13)$$

with weight functions $\Theta_i^*(\psi_c)$ quantifying the relative importance of $(\overline{a}_i)_{i \in [1..6]}$ and $\tilde{a}_7$ on describing the equilibrium thermochemical states

represented by $\overline{G}$ at the computed equilibrium temperature $T_{eq}$. Using notations from Eq. (8), $\forall i \in [1..7]$, $\Theta_i^*$ are expressed following:

$$\forall c \in [1..N_c], \Theta_i^*(\psi_c) = \left| \frac{\Theta_i(T_{eq})}{\sum_{j=1}^{7} \Theta_j(T_{eq})} \right|. \quad (14)$$

## 2.3. Inclusion of real molecules in the products $P^v$

As mentioned above, a distinctive feature of this new formalism is the possibility to include $N_p^r$ real combustion products $P^r$, together with $N_p^v$ virtual products $P^v$. An additional fitness term $f_{thermo}^R$ is then introduced in Eq. (11):

$$f_{\text{thermo}} = f_{\text{thermo}}^{\overline{M}} + f_{\text{thermo}}^{\overline{a}} + f_{\text{thermo}}^{R} \quad (15)$$

with

$$f_{\text{thermo}}^{R} = \sum_{c=1}^{N_c} \sum_{m \in P^r} \left| \frac{X_m^{eq^v}(\psi_c) - X_m^{eq^d}(\psi_c)}{X_m^{eq^d}(\psi_c)} \right| \quad (16)$$

where $X_m^{eq}$ is the mole fraction of species $m$ at chemical equilibrium. This additional term ensures a correct prediction of real products' equilibrium molar fractions. Note that, by including real species in products, the number of $N_p^v$ virtual products required to reach a desired accuracy is reduced. For instance, for hydrocarbon combustion, if one considers H₂O and CO₂ as real products ($N_p^r = 2$) in the mixture, the addition of a single virtual product $P_1$ ($N_p^v = 1$) is sufficient to complete the set of reduced scheme species.

## 2.4. Thermodynamical consistency

Additionally, the fitness function $f_{thermo}$ may be constrained to enforce the thermodynamic consistency of the considered reduced mixture. Consistency is assessed based on two fundamental principles: (i) species heat capacities must remain positive for all temperatures ($\forall T$, $C_{p,k}^\circ(T) > 0$); (ii) species entropy must increase monotonically with temperature ($\forall T$, $\frac{\partial S_k^\circ}{\partial T}(T) > 0$). These conditions must be verified for each virtual product $P_k$ over the combustion temperature range. Any solution violating these principles is discarded, with a large penalty assigned to $f_{thermo}$.

# 3. Thermodynamic optimization results

The optimization of $x_{thermo}$ is carried out using CMA-ES [15], an evolutionary strategy known for its remarkable convergence to global optima. Optimizations are carried out with a population of 20 individuals and an initial standard deviation $\sigma_0 = 0.2$. Equilibrium states are computed using CANTERA's equilibrate function at atmospheric pressure. Three sets of applications are here presented, for both Hydrogen and Dodecane combustion.

## 3.1. H₂-air, with two virtual products

This optimization procedure is first applied to H₂-Air combustion. The reference thermochemistry data align with those used in the state-of-the-art mechanism developed by Zhu et al., which is specifically tailored for H₂ combustion and encompasses 43 species [16]. Elements contained in the products are H and O, as nitrogen chemistry is not considered. According to the $N_p = N_e$ condition, the virtual scheme is made of two virtual products.

Equilibrium temperature is plotted as a function of equivalence ratio in Fig. 2. Learning points are represented by red dots. Validation is achieved by testing the resulting virtual scheme on conditions where the optimization was not applied, represented by blue dots. A perfect agreement with reference data is obtained, yielding an average relative error $\bar{\epsilon}_{rel} \approx 0.06\%$ across both learning and validation points.
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Adiabatic equilibrium temperature $T_{eq}$ for H$_3$-Air combustion over an equivalence ratio range of $\phi \in [0.4, 4]$; Comparison of the results obtained with optimized virtual species with that of detailed mechanism. Fresh gas injection at $T_0 = 300$ K and atmospheric pressure.</figcaption></figure>

### 3.2. $nC_{12}H_{26}$-air, with three virtual products

The method is then validated on a more complex fuel: $nC_{12}H_{26}$. The reference thermochemistry data is now extracted from the 413-species C0-C16 CRECK model [17], designed for heavy hydrocarbon combustion. The optimized reduced scheme, significantly lighter, consists only of $nC_{12}H_{26}$, O$_2$, N$_2$, and three virtual products. The study considers an equivalence ratio range of $\phi \in [0.7, 1.4]$ and a fresh gas temperature $T_0 = 400$ K. Fig. 3 demonstrates very good agreement in the targeted equilibrium thermodynamic mixture properties; heat capacity $C_p$, enthalpy $H$, entropy $S$, and Gibbs energy $G$. For each quantity, the relative errors $\bar{e}_{rel}$ remain largely below 1%. Consistently, the virtual scheme yields excellent predictions of the adiabatic equilibrium temperature $T_{eq}$, as shown in Fig. 4, with an estimated relative error $\bar{e}_{rel} \approx 0.3\%$.

### 3.3. $nC_{12}H_{26}$-air, with two real and one virtual products

The ability of the virtual approach to incorporate real combustion products for $nC_{12}H_{26}$-Air combustion is examined. Here, two real products H$_2$O and CO$_2$, and one virtual product P'$_1$ are considered in the virtual scheme. Therefore, a second optimization for $nC_{12}H_{26}$-Air combustion was performed. H$_2$O and CO$_2$ mass fractions at equilibrium obtained by the virtual scheme are compared to the reference scheme results for different equivalence ratios in Fig. 5. A good agreement between the results is obtained, showing the capacity of the method in the prediction of real species. Although not shown here, the model perfectly retains its ability to accurately reproduce mixture equilibrium states.

## 4. Virtual chemistry optimization

### 4.1. Virtual reaction scheme architecture

**Initiation & termination reactions**

The ensemble of real and virtual combustion products is noted $P = P^r + P^v$. $P$ is therefore composed of $N_p = N_p^r + N_p^v$ species.

A simplified case with a single intermediate species I ($N_p^v = 1$) is here considered. This structure has proven sufficient for retrieving flame properties of interest for various fuels [7]. However, capturing more complex chemical phenomena (such as ignition delay time) requires additional intermediates/reactions, as demonstrated by Pacheco & al. [18].

Inspired by previous virtual chemistry formalisms [7,11], an initiation reaction consumes the fuel F and oxidizer Ox to partly form the optimized intermediate species I and the combustion products in an irreversible reaction:

$$F + S \text{Ox} \longrightarrow \alpha''_{1,1} I + \sum_{P_k \in P} \alpha''_{1,P_k} P_k \quad (R_1)$$

where $S$ is the molar stoichiometric coefficient of the global combustion reaction. Since $N_p = N_e$ has been assumed, the linear system of $N_p$ unknown $\alpha''_{1,P_k}$ is then fully determined from the $N_e$ elemental balance equations. I and $\alpha''_{1,1}$ are optimized.

Finally, I is completely consumed through the following termination reaction:

$$I + \sum_{P_k \in P} \alpha'_{2,P_k} P_k \longrightarrow \sum_{P_k \in P} \alpha''_{2,P_k} P_k \quad (R_2)$$

where stoichiometric coefficients $\alpha'_{2,P_k}$ and $\alpha''_{2,P_k}$ are again determined by solving a linear system to conserve elements. Each species $P_k$ is either a reactant ($\alpha''_{2,P_k} = 0$) or a product ($\alpha'_{2,P_k} = 0$) in the reaction. The above two chemical reactions assume combustion is a complete process. Equilibrium reactions are therefore required to provide the chemical pathways needed to reach the previously optimized equilibrium thermochemical states.

**Equilibrium reactions**

Equilibrium states are characterized by their species mass fractions $Y_k^{eq}$. They are composed of fuel F, oxidizer Ox and dilutant D together with the $N_p$ products. Since D is non-reactive, its equilibrium mass fraction $Y_D^{eq}$ is already known from fresh gases. The remaining $N_p + 2$ equilibrium mass fractions are constrained by $N_e$ elemental conservation equations. Using $N_p = N_e$ again, the number of independent chemical reactions at equilibrium $N_r^{eq}$ is:

$$N_r^{eq} = N_p + 2 - N_e = 2. \quad (17)$$

The proposed equilibrium reactions to complete the virtual mechanism are therefore:

$$\begin{align*} & \text{F} + \sum_{P_k \in P} \alpha'_{3,P_k} P_k & \rightleftharpoons & \sum_{P_k \in P} \alpha''_{3,P_k} P_k \quad (R_3) \\ & \text{Ox} + \sum_{P_k \in P} \alpha'_{4,P_k} P_k & \rightleftharpoons & \sum_{P_k \in P} \alpha''_{4,P_k} P_k \quad (R_4) \end{align*}$$

Again, in each reaction $r \in \{(R_3), (R_4)\}$, stoichiometric coefficients ($\alpha'_{r,P_k}, \alpha''_{r,P_k}$) are computed and imposed to conserve elements. A distinctive feature of the mechanism addresses pyrolysis: for fuels such as heavy hydrocarbons, the equilibrium reaction involving F is made irreversible to allow complete fuel consumption under rich conditions.

The corresponding optimization vector $\mathbf{x}_{\text{Scheme}}$ defining the virtual reaction model reads:

$$\mathbf{x}_{\text{Scheme}} = \left( \underbrace{(a_2^1, a_3^1, a_4^1, a_5^1, a_6^1, a_7^1)}_{\text{Thermodynamic coefficients of I}}, \underbrace{(N_C^1, N_O^1, N_H^1, \dots)}, \underbrace{\alpha''_{1,11}}_{\text{Stoichi. coeff.}} \right). \quad (18)$$

### 4.2. Reaction rates

For a given reaction $r$, the rate of progress $q^r$ is written using forward and backward rate constants $k_f^r$ and $k_b^r$ following:

$$q^r = k_f^r \prod_{k=1}^{N_{\text{react}}} [X_k]^{n_{f,k}} - k_b^r \prod_{k=1}^{N_{\text{prod}}} [X_k]^{n_{b,k}} \quad (19)$$

where $[X_k]$ denotes the concentration of species $k$. $N_{\text{react}}$ and $N_{\text{prod}}$ respectively represent the number of reactants and products in reaction $r$. $n_{f,k}$ and $n_{b,k}$ symbolize the forward and backward orders of reactions of species $k$ in reaction $r$. The reaction orders are either set as the associated stoichiometric coefficient, or optimized, depending on whether the reaction is reversible, or total [13].

The forward rate constant $k_f^r$ uses the classical Arrhenius formula:

$$\forall r \in [1..4], \quad k_f^r(T) = A^r T^{b'} e^{-E_a^r / RT} \quad (20)$$

$A^r$, $b^r$ and $E_a^r$ being respectively the pre-exponential constant, temperature exponent and activation energy of reaction ($R_r$). The backward
<figure><img src="image_6.png" /><figcaption>Fig. 3. Mixture's heat capacity $C_p$, enthalpy $H$, entropy $S$ & Gibbs energy $G$ at chemical equilibrium for $nC_{12}H_{26}$-Air combustion over an equivalence ratio range of $\phi \in [0.7, 1.4]$; Comparison of the results obtained with optimized virtual species with that of detailed mechanism. Associated relative errors across both learning and validation points $\bar{\epsilon}_{rel}$ are provided. Fresh gas injection at $T_0 = 400$ K and atmospheric pressure.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Adiabatic equilibrium temperature $T_{eq}$ for $nC_{12}H_{26}$-Air combustion over an equivalence ratio range of $\phi \in [0.7, 1.4]$; Comparison of the results obtained with optimized virtual species with that of detailed mechanism. Fresh gas injection at $T_0 = 400$ K and atmospheric pressure.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Equilibrium mass fractions of reduced mixture combustion products $H_2O$, $CO_2$ & $P'_1$ for $nC_{12}H_{26}$ combustion over an equivalence ratio range of $\phi \in [0.7, 1.4]$. Predictions of burnt gases mass fractions of $H_2O$, $CO_2$ are compared to the detailed mixture. Fresh gas injection at $T_0 = 400$ K and atmospheric pressure.</figcaption></figure>

rate constant $k'_b$ of any equilibrium reaction is determined by its forward rate constant $k'_f$ and its equilibrium constant $K'_{eq}$ through:

$$k'_{b}(T) = \frac{k'_{f}(T)}{K'_{eq}(T)} = \frac{k'_{f}(T)}{e^{-\Delta G'_{r}(T)/RT}} \quad (21)$$

where $\Delta G'_{r}(T)$ is Gibbs free energy of reaction $r$ at standard conditions. The chemical rates optimization vector $\mathbf{x}_{Arrhenius}$ thus reads:

$$\mathbf{x}_{\text{Arrhenius}} = \left( (A^r, b^r, E_a^r) \dots \forall r \in [1..4], \underbrace{(n_F, n_{Ox}, n_1 \dots)}_{\text{Orders of reaction}} \right). \quad (22)$$

so that the final optimization vector $\mathbf{x}_{chem}$ is:

$$\mathbf{x}_{chem} = \left( \mathbf{x}_{Scheme}, \mathbf{x}_{Arrhenius} \right). \quad (23)$$

### 4.3. Fitness function

The fitness function $f_{chem}$ is formulated to capture both temperature $T$ and heat release rate $\dot{\omega}_T$. It is expressed as:

$$f_{chem} = \gamma_1 f_{chem}^T + \gamma_2 f_{chem}^{\omega_T} \quad (24)$$

with

$$f_{chem}^T = \sum_{c=1}^{N_c} \left| \frac{T^v(\psi_c) - T^d(\psi_c)}{T^d(\psi_c)} \right| \quad (25)$$

$$f_{chem}^{\omega_T} = \sum_{c=1}^{N_c} \left| \frac{\dot{\omega}_T^v(\psi_c) - \dot{\omega}_T^d(\psi_c)}{\dot{\omega}_T^d(\psi_c)} \right|$$

The ratio $\gamma_2/\gamma_1$ is computed to represent the ratio between reactive to thermal thicknesses, ensuring comparable contributions of these two metrics to $f_{chem}$. The constraints defined in Section 2.4 are applied to ensure the thermodynamic consistency of the optimized intermediate I.

## 5. Chemistry optimization results

A reference kinetics database is obtained by computing an ensemble of 0-D Perfectly-Stirred Reactors (PSR). $H_2$-Air reactors are
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Comparison of reference and optimized PSR temperature *T* and heat release rate $\dot{\omega}_T$ profiles at stoichiometry. Top: H₂-Air combustion. Bottom: nC₁₂H₂₆-Air combustion.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Comparison of reference and optimized temperature dependency on equivalence ratio for three fixed residence times $\tau$. Top: H₂-Air combustion. Bottom: nC₁₂H₂₆-Air combustion.</figcaption></figure>

computed using the Zhu & al. detailed mechanism (43 species, 312 reactions) [16] for equivalence ratio ranging in [0.4, 4], while the learning database for nC₁₂H₂₆-Air is computed with the CRECK mechanism (413 species, 14,922 reactions) [17] for $\phi \in [0.7, 1.4]$. Optimizations are performed using CMA-ES with the same hyperparameters as in the thermodynamics optimization step.

Fig. 6 presents the temperature and heat release rate as a function of residence time at stoichiometry. For H₂-Air combustion, the virtual scheme (6 species, 4 reactions) shows a near perfect agreement with detailed chemistry. A comparison is also made with the Boivin skeletal scheme [19] (9 species, 12 reactions), which provides similar results: The relative errors in temperature *T* compared to detailed chemistry at $\phi = 1$ are $\bar{\varepsilon}_{rel} = 2.4\%$ for the virtual scheme, and $\bar{\varepsilon}_{rel} = 3.5\%$ for the Boivin scheme.

For nC₁₂H₂₆-Air combustion, the virtual scheme (7 species, 4 reactions) is evaluated against the detailed mechanism CRECK [17] and the popular BFER global scheme [5]. The 4-reaction virtual scheme largely outperforms BFER in capturing temperature and heat release rate: At stoichiometry, the virtual mechanism reproduces reference heat release rates with $\bar{\varepsilon}_{rel} = 3.6\%$, a dramatic improvement over the 114% average relative error of BFER, as seen in Fig. 6. Notably, BFER employs C₁₀H₂₀ as a surrogate for kerosene, leading to slight differences in fresh gas compositions. Fig. 7 extends the analysis to equivalence ratio variations at three residence times $\tau$ covering different chemical timescales. For both fuels, similar conclusions emerge. A remarkable feature of virtual schemes is their ability to naturally converge towards reference thermochemical equilibrium states for large residence times.
This demonstrates the robustness of reversible reactions in reaching
optimized equilibrium states.

## 6. Conclusions

A new framework has been developed for virtual chemistry, intro-
ducing a compact scheme architecture describing combustion dynamics
towards chemical equilibria without requiring the tabulation of kinetic
parameters. Additionally, the potential of the method to predict the
production of specified molecules, such as H₂O or CO₂, to partially
describe burnt gases has been demonstrated. This innovative formalism
has been applied to both Hydrogen-Air and Dodecane-Air combustion,
resulting in 4-reaction kinetic schemes reproducing the evolution of
temperature and heat release rate of complex detailed mechanisms with
an excellent agreement. By adhering to standard NASA polynomial
approximations for thermodynamic properties and an Arrhenius-like
closure for kinetic rate constants, integrating virtual chemistry into any
reactive flow solver is now straightforward. A forthcoming article will
extend this framework to flame prediction, with a particular emphasis
on chemical stiffness and computational cost estimates to illustrate the
method's efficiency potential. Future developments will also address
pollutant formation within this formalism.

## Novelty and significance statement

This study introduces a new generalized formalism for virtual chem-
istry, creating universal 4-reaction kinetic mechanisms applicable to
various fuel. By optimizing virtual products as linear combinations of
real molecules to match reference thermodynamic properties of burnt
gases, this method accurately reproduces detailed chemistry equilib-
rium states across various conditions. The inclusion of real species like
H₂O and CO₂ extends virtual chemistry's capabilities, enabling predic-
tion of global combustion products. Linking these species, a general
4-reaction kinetic mechanism is then proposed, featuring total reactions
with optimized intermediate radicals and equilibrium reactions to reach
previously optimized states. Optimized using 0-D reactors on Cantera,
the resulting mechanisms for H₂ and nC₁₂H₂₆ demonstrate excellent
agreement with complex chemistry and outperform traditional global
schemes, particularly for complex fuels like nC₁₂H₂₆. By maintaining
a structure identical to detailed mechanisms, these mechanisms can
therefore be seamlessly integrated into any reactive flow solver.

## CRediT authorship contribution statement

**Matthieu Préteseille:** Writing – original draft, Visualization, Validation, Software, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Tan-Phong Luu:** Investigation, Formal analysis, Conceptualization. **Nasser Darabiha:** Writing – original draft, Validation, Supervision, Software, Methodology, Investigation, Conceptualization. **Benoît Fiorina:** Writing – original draft, Validation, Supervision, Resources, Project administration, Methodology, Investigation, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing finan-
cial interests or personal relationships that could have appeared to
influence the work reported in this paper.

## Acknowledgments

The authors acknowledge support from the MONTHY project funded
by the French National Research Agency (ANR) under France 2030
program and reference ANR 22-PEHY-0014. A CC-BY public copyright
license has been applied by the authors to the present document and
will be applied to all subsequent versions up to the Author Accepted
Manuscript arising from this submission, in accordance with the grant's
open access conditions. Authors also acknowledge support from Ariane-
Group. This research was conducted using HPC resources from IDRIS,
TGCC, and CINES under the allocation A0132B10253 by GENCI (Grand
Équipement National de Calcul Intensif). Lastly, the authors acknowl-
edge the 8th Extreme CFD Workshop for providing a collaborative and
stimulating environment that contributed to the development of this
work.

## Appendix A. Supplementary data

The virtual mechanisms for both Hydrogen and Dodecane are at-
tached in supplementary materials, using the CANTERA-compatible
YAML format.

Supplementary material related to this article can be found online
at https://doi.org/10.1016/j.proci.2025.105894.

## References

[1] G. Sorrentino, G.B. Ariemma, F. Ferraro, B. Fiorina, Including detailed chemistry features in the modeling of emerging low-temperature reactive flows: A review on the application to diluted and MILD combustion systems, Appl. Energy Combust. Sci. 20 (2024) 100291.

[2] N. Sekularac, T. Lesaffre, D. Laera, L. Cicquel, Large eddy simulations of n-heptane and n-dodecane binary blends in swirling multi-component spray flames, Proc. Combust. Inst. 40 (1) (2024) 105201.

[3] J. van Oijen, A. Donini, R. Bastiaans, J. ten Thije Boonkkamp, L. de Goey, State-of-the-art in premixed combustion modeling using flamelet generated manifolds, Prog. Energy Combust. Sci. 57 (2016) 30–74.

[4] U. Maas, S. Pope, Simplifying chemical kinetics: Intrinsic low-dimensional manifolds in composition space, Combust. Flame 88 (3) (1992) 239–264.

[5] B. Franzelli, E. Riber, M. Sanjosé, T. Poinsot, A two-step chemical scheme for kerosene-air premixed flames, Combust. Flame 157 (7) (2010) 1364–1373.

[6] A. Millán-Merino, S. Taileb, P. Boivin, A new method for systematic 1-step chemistry reduction applied to hydrocarbon combustion, Proc. Combust. Inst. 39 (1) (2023) 745–753.

[7] M. Cailler, N. Darabiha, B. Fiorina, Development of a virtual optimized chemistry method. Application to hydrocarbon/air combustion, Combust. Flame 211 (2020) 281–302.

[8] A.F. Pacheco, A.A. Oliveira, B. Fiorina, Development of a virtual kinetics chemistry for prediction of ignition delay time, Combust. Flame 246 (2022) 112444.

[9] G. Maio, M. Cailler, A. Cuoci, B. Fiorina, A virtual chemical mechanism for prediction of NO emissions from flames, Combust. Theory Model. 24 (5) (2020) 872–902.

[10] H. Maldonado Colmán, A. Cuoci, N. Darabiha, B. Fiorina, A virtual chemistry model for soot prediction in flames including radiative heat transfer, Combust. Flame 238 (2022) 111879.

[11] G. Maio, M. Cailler, R. Mercier, B. Fiorina, Virtual chemistry for temperature and CO prediction in LES of non-adiabatic turbulent flames, Proc. Combust. Inst. 37 (2) (2019) 2591–2599.

[12] R.J. Kee, F.M. Rupley, J.A. Miller, Chemkin-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics, Tech. Rep., Sandia National Lab. (SNL-CA), Livermore, CA (United States), 1989.

[13] D.G. Goodwin, H.K. Moffat, I. Schoegl, R.L. Speth, B.W. Weber, Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes, 2023, Version 3.0.0.

[14] B.J. McBride, S. Gordon, M.A. Reno, Coefficients for Calculating Thermodynamic and Transport Properties of Individual Species, NASA Technical Memorandum 4513, National Aeronautics and Space Administration, OH (United States), 1993.

[15] N. Hansen, A. Ostermeier, Completely derandomized self-adaptation in evolution strategies, Evol. Comput. 9 (2) (2001) 159–195.
[16] Y. Zhu, H.J. Curran, S. Girhe, Y. Murakami, H. Pitsch, K. Senecal, L. Yang, C.-W. Zhou, The combustion chemistry of ammonia and ammonia/hydrogen mixtures: A comprehensive chemical kinetic modeling study, Combust. Flame 260 (2024) 113239.

[17] E. Ranzi, A. Frassoldati, A. Stagni, M. Pelucchi, A. Cuoci, T. Faravelli, Reduced kinetic schemes of complex reaction systems: Fossil and biomass-derived transportation fuels, Int. J. Chem. Kinet. 46 (9) (2014) 512–542.

[18] A.F. Pacheco, A.A. Oliveira, B. Fiorina, Development of a virtual kinetics chemistry for prediction of ignition delay time, Combust. Flame 246 (2022) 112444.

[19] P. Boivin, Reduced-Kinetic Mechanisms for Hydrogen and Syngas Combustion Including Autoignition (Ph.D. thesis), Universidad Carlos III de Madrid, 2011.