# Thermochemical equilibrium module#

In this section, you will find the documentation of the kernel of the code, which is used to obtain the chemical equilibrium composition for a desired thermochemical transformation, e.g., constant enthalpy and pressure. It also includes routines to compute chemical equilibrium assuming a complete combustion and the calculation of the thermodynamic derivates. The code stems from the minimization of the free energy of the system by using Lagrange multipliers combined with a Newton-Raphson method, upon condition that initial gas properties are defined by two functions of states, e.g., temperature and pressure.

Note

The kernel of the code is based on Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

## Thermodynamic derivatives#

All thermodynamic first derivatives can be obtained with any set of three independent first derivatives . Combustion Toolbox computes all thermodynamic first derivatives from $$(\partial \text{ln } V/\partial \text{ln } T)_p$$, $$(\partial \text{ln } V/\partial \text{ln } p)_T$$, and $$(\partial h/\partial T)_p = c_p$$. Considering the ideal equation of state

\begin{equation} pV = nRT \end{equation}

and applying logarithms to both sides

\begin{equation} \text{ln } V = \text{ln } n + \text{ln } R + \text{ln } T - \text{ln } p \end{equation}

\begin{equation} \left(\dfrac{\partial \text{ln } V }{\partial \text{ln } T}\right)_p = 1 + \left(\dfrac{\partial \text{ln } n }{\partial \text{ln } T}\right)_p, \end{equation}
\begin{equation} \left(\dfrac{\partial \text{ln } V }{\partial \text{ln } p}\right)_T = -1 + \left(\dfrac{\partial \text{ln } n }{\partial \text{ln } p}\right)_T. \end{equation}

To compute $$c_p$$ we have to distinguish between the frozen contribution and the reaction contribution

\begin{equation} c_p = c_{p,f} + c_{p,r} \end{equation}

given by the following relations

\begin{equation} c_{p,f} = \sum\limits_{j=1}^{\text{NS}} n_j c_{p,f}^\circ, \end{equation} \begin{equation} c_{p,r} = \dfrac{1}{T}\left[\sum\limits_{j=1}^{\text{NS}} [1 + \delta_j(n_j - 1)] h_j^\circ\left(\dfrac{\partial \eta_j}{\partial \text{ln } T} \right)\right], \end{equation}

with $$\eta_j = \text{ln } n_j$$ and $$\delta_j = 1$$ for $$j=1,\dots,NG$$ (non-condensed species), and $$\eta_j = n_j$$ and $$\delta_j = 0$$ for $$j = NG + 1, \dots, NS$$ (condensed species).

### Derivatives with respect to temperature#

\begin{aligned} \delta_j\left(\dfrac{\partial \eta_j }{\partial \text{ln } T}\right)_p - \sum\limits_{i = 1}^{\text{NE}} a_{ij} \left(\dfrac{\partial \pi_i }{\partial \text{ln } T}\right)_p - \delta_j\left(\dfrac{\partial \text{ln } n}{\partial \text{ln } T}\right)_p &= \dfrac{h_j^\circ}{RT}, \quad &\text{for } j=1, \dots, \text{NS}\\ \sum\limits_{j = 1}^{\text{NS}} a_{ij} [1 + \delta_j(n_j - 1)] \left(\dfrac{\partial \eta_j }{\partial \text{ln } T}\right)_p &= 0, \quad &\text{for } i=1, \dots, \text{NE}\\ \sum\limits_{j = 1}^{\text{NG}} n_j \left(\dfrac{\partial \eta_j }{\partial \text{ln } T}\right)_p - n \left(\dfrac{\partial \text{ln } n}{\partial \text{ln } T}\right)_p&= 0, \end{aligned}

### Derivatives with respect to pressure#

\begin{aligned} \delta_j\left(\dfrac{\partial \eta_j }{\partial \text{ln } p}\right)_T - \sum\limits_{i = 1}^{\text{NE}} a_{ij} \left(\dfrac{\partial \pi_i }{\partial \text{ln } p}\right)_T - \delta_j\left(\dfrac{\partial \text{ln } n}{\partial \text{ln } p}\right)_T &= -\delta, \quad &\text{for } j=1, \dots, \text{NS}\\ \sum\limits_{j = 1}^{\text{NS}} a_{ij} [1 + \delta_j(n_j - 1)] \left(\dfrac{\partial \eta_j }{\partial \text{ln } p}\right)_T &= 0, \quad &\text{for } i=1, \dots, \text{NE}\\ \sum\limits_{j = 1}^{\text{NG}} n_j \left(\dfrac{\partial \eta_j }{\partial \text{ln } p}\right)_T - n \left(\dfrac{\partial \text{ln } n}{\partial \text{ln } p}\right)_T&= 0. \end{aligned}

Routines
complete_combustion(self, mix, phi)#

Solve chemical equilibrium for CHNO mixtures assuming a complete combustion using the equilibrium constants method

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• mix (struct) – Properties of the initial mixture

• phi (float) – Equivalence ratio [-]

Returns:

Tuple containing

• moles (float): Equilibrium composition [moles] at defined temperature

• species (str): Species considered in the complemte combustion model

Example

[moles, species] = complete_combustion(self, self.PS.strR{1}, 0.5)

equilibrate(self, mix1, pP, varargin)#

Obtain properties at equilibrium for the set thermochemical transformation

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• mix1 (struct) – Properties of the initial mixture

• pP (float) – Pressure [bar]

Optional Args:

mix2 (struct): Properties of the final mixture (previous calculation)

Returns:

mix2 (struct) – Properties of the final mixture

Example

mix2 = equilibrate(self, self.PS.strR{1}, 1.01325)

equilibrate_T(self, mix1, pP, TP, varargin)#

Obtain equilibrium properties and composition for the given temperature [K] and pressure [bar]

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• mix1 (struct) – Properties of the initial mixture

• pP (float) – Pressure [bar]

• TP (float) – Temperature [K]

Optional Args:

guess_moles (float): Mixture composition [mol] of a previous computation

Returns:

mix2 (struct) – Properties of the final mixture

Example

mix2 = equilibrate(self, self.PS.strR{1}, 1.01325, 3000)

equilibrate_T_tchem(self, mix1, pP, TP)#

Obtain equilibrium properties and composition for the given temperature [K] and pressure [bar] assuming a calorically perfect gas

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• mix1 (struct) – Properties of the initial mixture

• pP (float) – Pressure [bar]

• TP (float) – Temperature [K]

Returns:

mix2 (struct) – Properties of the final mixture assuming a calorically perfect gas

Example

mix2 = equilibrate_T_tchem(self, self.PS.strR{1}, 1.01325, 3000)

equilibrium_dT(J, N0, A0, NE, ind_nswt, ind_swt, ind_elem, H0RT)#

Obtain thermodynamic derivative of the moles of the species and of the moles of the mixture respect to temperature from a given composition [moles] at equilibrium

Parameters:
• J (float) – Matrix J to solve the linear system J*x = b

• N0 (float) – Equilibrium composition [moles]

• A0 (float) – Stoichiometric matrix

• NE (float) – Temporal total number of elements

• ind (float) – Temporal index of species in the final mixture

• ind_nswt (float) – Temporal index of gaseous species in the final mixture

• ind_swt (float) – Temporal index of condensed species in the final mixture

• ind_elem (float) – Temporal index of elements in the final mixture

• H0RT (float) – Dimensionless enthalpy

Returns:

Tuple containing

• dNi_T (float): Thermodynamic derivative of the moles of the species respect to temperature

• dN_T (float): Thermodynamic derivative of the moles of the mixture respect to temperature

Example

[dNi_T, dN_T] = equilibrium_dT(J, N0, A0, NE, ind, ind_nswt, ind_swt, ind_elem, H0RT)

equilibrium_dT_large(self, moles, T, A0, NG, NS, NE, ind, ind_nswt, ind_swt, ind_E)#

Obtain thermodynamic derivative of the moles of the species and of the moles of the mixture respect to temperature from a given composition [moles] at equilibrium

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• moles (float) – Equilibrium composition [moles]

• T (float) – Temperature [K]

• A0 (float) – Stoichiometric matrix

• NG (float) – Temporal total number of gaseous species

• NS (float) – Temporal total number of species

• NE (float) – Temporal total number of elements

• ind (float) – Temporal index of species in the final mixture

• ind_nswt (float) – Temporal index of gaseous species in the final mixture

• ind_swt (float) – Temporal index of condensed species in the final mixture

• ind_E (float) – Temporal index of elements in the final mixture

Returns:

Tuple containing

• dNi_T (float): Thermodynamic derivative of the moles of the species respect to temperature

• dN_T (float): Thermodynamic derivative of the moles of the mixture respect to temperature

Example

[dNi_T, dN_T] = equilibrium_dT_large(self, moles, T, A0, NG, NS, NE, ind, ind_nswt, ind_swt, ind_E)

equilibrium_dp(J, N0, A0, NE, ind_nswt, ind_swt, ind_elem)#

Obtain thermodynamic derivative of the moles of the species and of the moles of the mixture respect to pressure from a given composition [moles] at equilibrium

Parameters:
• J (float) – Matrix J to solve the linear system J*x = b

• N0 (float) – Equilibrium composition [moles]

• A0 (float) – Stoichiometric matrix

• NE (float) – Temporal total number of elements

• ind (float) – Temporal index of species in the final mixture

• ind_nswt (float) – Temporal index of gaseous species in the final mixture

• ind_swt (float) – Temporal index of condensed species in the final mixture

• ind_elem (float) – Temporal index of elements in the final mixture

Returns:

Tuple containing

• dNi_p (float): Thermodynamic derivative of the moles of the species respect to pressure

• dN_p (float): Thermodynamic derivative of the moles of the mixture respect to pressure

Example

[dNi_p, dN_p] = equilibrium_dp(J, N0, A0, NE, ind, ind_nswt, ind_swt, ind_elem)

equilibrium_dp_large(self, N0, A0, NG, NS, NE, ind, ind_nswt, ind_swt, ind_E)#

Obtain thermodynamic derivative of the moles of the species and of the moles of the mixture respect to pressure from a given composition [moles] at equilibrium

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• N0 (float) – Equilibrium composition [moles]

• A0 (float) – Stoichiometric matrix

• NG (float) – Temporal total number of gaseous species

• NS (float) – Temporal total number of species

• NE (float) – Temporal total number of elements

• ind (float) – Temporal index of species in the final mixture

• ind_nswt (float) – Temporal index of gaseous species in the final mixture

• ind_swt (float) – Temporal index of condensed species in the final mixture

• ind_E (float) – Temporal index of elements in the final mixture

Returns:

Tuple containing

• dNi_p (float): Thermodynamic derivative of the moles of the species respect to pressure

• dN_p (float): Thermodynamic derivative of the moles of the mixture respect to pressure

Example

[dNi_p, dN_p] = equilibrium_dp_large(self, N0, A0, NG, NS, NE, ind, ind_nswt, ind_swt, ind_E)

equilibrium_gibbs(self, pP, TP, mix1, guess_moles)#

Obtain equilibrium composition [moles] for the given temperature [K] and pressure [bar]. The code stems from the minimization of the free energy of the system by using Lagrange multipliers combined with a Newton-Raphson method, upon condition that initial gas properties are defined by temperature and pressure.

The algorithm implemented take advantage of the sparseness of the upper left submatrix obtaining a matrix J of size NE + NS - NG + 1.

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• pP (float) – Pressure [bar]

• TP (float) – Temperature [K]

• mix1 (struct) – Properties of the initial mixture

• guess_moles (float) – Mixture composition [mol] of a previous computation

Returns:

Tuple containing

• N0 (float): Composition matrix [n_i, FLAG_CONDENSED_i] for the given temperature [K] and pressure [bar] at equilibrium

• dNi_T (float): Thermodynamic derivative of the moles of the species respect to temperature

• dN_T (float): Thermodynamic derivative of the moles of the mixture respect to temperature

• dNi_p (float): Thermodynamic derivative of the moles of the species respect to pressure

• dN_p (float): Thermodynamic derivative of the moles of the mixture respect to pressure

• ind (float): List of chemical species indices

• STOP (float): Relative error in moles of species [-]

Examples

• N0 = equilibrium_gibbs(self, 1.01325, 3000, self.PS.strR{1}, [])

• [N0, dNi_T, dN_T, dNi_p, dN_p, ind, STOP, STOP_ions] = equilibrium_gibbs(self, 1.01325, 3000, self.PS.strR{1}, [])

equilibrium_gibbs_eos(self, pP, TP, mix1, guess_moles)#

Obtain equilibrium composition [moles] for the given temperature [K] and pressure [bar]. The code stems from the minimization of the free energy of the system by using Lagrange multipliers combined with a Newton-Raphson method, upon condition that initial gas properties are defined by temperature and pressure.

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• pP (float) – Pressure [bar]

• TP (float) – Temperature [K]

• mix1 (struct) – Properties of the initial mixture

• guess_moles (float) – Mixture composition [mol] of a previous computation

Returns:

Tuple containing

• N0 (float): Equilibrium composition [moles] for the given temperature [K] and pressure [bar]

• dNi_T (float): Thermodynamic derivative of the moles of the species respect to temperature

• dN_T (float): Thermodynamic derivative of the moles of the mixture respect to temperature

• dNi_p (float): Thermodynamic derivative of the moles of the species respect to pressure

• dN_p (float): Thermodynamic derivative of the moles of the mixture respect to pressure

• STOP (float): Relative error in moles of species [-]

• STOP_ions (float): Relative error in moles of ionized species [-]

Example

[N0, dNi_T, dN_T, dNi_p, dN_p, STOP, STOP_ions] = equilibrium_gibbs_eos(self, pP, TP, mix1, guess_moles)

equilibrium_gibbs_large(self, pP, TP, mix1, guess_moles)#

Obtain equilibrium composition [moles] for the given temperature [K] and pressure [bar]. The code stems from the minimization of the free energy of the system by using Lagrange multipliers combined with a Newton-Raphson method, upon condition that initial gas properties are defined by temperature and pressure.

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• pP (float) – Pressure [bar]

• TP (float) – Temperature [K]

• mix1 (struct) – Properties of the initial mixture

• guess_moles (float) – Mixture composition [mol] of a previous computation

Returns:

Tuple containing

• N0 (float): Equilibrium composition [moles] for the given temperature [K] and pressure [bar]

• dNi_T (float): Thermodynamic derivative of the moles of the species respect to temperature

• dN_T (float): Thermodynamic derivative of the moles of the mixture respect to temperature

• dNi_p (float): Thermodynamic derivative of the moles of the species respect to pressure

• dN_p (float): Thermodynamic derivative of the moles of the mixture respect to pressure

• ind (float): List of chemical species indices

• STOP (float): Relative error in moles of species [-]

• STOP_ions (float): Relative error in moles of ionized species [-]

Examples

• N0 = equilibrium_gibbs_large(self, 1.01325, 3000, self.PS.strR{1}, [])

• [N0, dNi_T, dN_T, dNi_p, dN_p, ind, STOP, STOP_ions] = equilibrium_gibbs_large(self, 1.01325, 3000, self.PS.strR{1}, [])

equilibrium_helmholtz(self, vP, TP, mix1, guess_moles)#

Obtain equilibrium composition [moles] for the given temperature [K] and volume [m3]. The code stems from the minimization of the free energy of the system by using Lagrange multipliers combined with a Newton-Raphson method, upon condition that initial gas properties are defined by temperature and volume.

The algorithm implemented take advantage of the sparseness of the upper left submatrix obtaining a matrix J of size NE + NS - NG.

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• vP (float) – Volume [m3]

• TP (float) – Temperature [K]

• mix1 (struct) – Properties of the initial mixture

• guess_moles (float) – Mixture composition [mol] of a previous computation

Returns:

Tuple containing

• N0 (float): Equilibrium composition [moles] for the given temperature [K] and pressure [bar]

• dNi_T (float): Thermodynamic derivative of the moles of the species respect to temperature

• dN_T (float): Thermodynamic derivative of the moles of the mixture respect to temperature

• dNi_p (float): Thermodynamic derivative of the moles of the species respect to pressure

• dN_p (float): Thermodynamic derivative of the moles of the mixture respect to pressure

• ind (float): List of chemical species indices

• STOP (float): Relative error in moles of species [-]

• STOP_ions (float): Relative error in moles of ionized species [-]

Examples

• N0 = equilibrium_helmholtz(self, 0.0716, 3000, self.PS.strR{1}, [])

• [N0, dNi_T, dN_T, dNi_p, dN_p, ind, STOP, STOP_ions] = equilibrium_helmholtz(self, 0.0716, 3000, self.PS.strR{1}, [])

equilibrium_helmholtz_eos(self, vP, TP, mix1, guess_moles)#

Obtain equilibrium composition [moles] for the given temperature [K] and volume [m3]. The code stems from the minimization of the free energy of the system by using Lagrange multipliers combined with a Newton-Raphson method, upon condition that initial gas properties are defined by temperature and volume.

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• vP (float) – Volume [m3]

• TP (float) – Temperature [K]

• mix1 (struct) – Properties of the initial mixture

• guess_moles (float) – Mixture composition [mol] of a previous computation

Returns:

Tuple containing

• N0 (float): Equilibrium composition [moles] for the given temperature [K] and pressure [bar]

• dNi_T (float): Thermodynamic derivative of the moles of the species respect to temperature

• dN_T (float): Thermodynamic derivative of the moles of the mixture respect to temperature

• dNi_p (float): Thermodynamic derivative of the moles of the species respect to pressure

• dN_p (float): Thermodynamic derivative of the moles of the mixture respect to pressure

• STOP (float): Relative error in moles of species [-]

• STOP_ions (float): Relative error in moles of ionized species [-]

Example

[N0, dNi_T, dN_T, dNi_p, dN_p, STOP, STOP_ions] = equilibrium_helmholtz_eos(self, vP, TP, mix1, guess_moles)

equilibrium_helmholtz_large(self, vP, TP, mix1, guess_moles)#

Obtain equilibrium composition [moles] for the given temperature [K] and volume [m3]. The code stems from the minimization of the free energy of the system by using Lagrange multipliers combined with a Newton-Raphson method, upon condition that initial gas properties are defined by temperature and volume.

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:
• self (struct) – Data of the mixture, conditions, and databases

• vP (float) – Volume [m3]

• TP (float) – Temperature [K]

• mix1 (struct) – Properties of the initial mixture

• guess_moles (float) – Mixture composition [mol] of a previous computation

Returns:

Tuple containing

• N0 (float): Equilibrium composition [moles] for the given temperature [K] and pressure [bar]

• dNi_T (float): Thermodynamic derivative of the moles of the species respect to temperature

• dN_T (float): Thermodynamic derivative of the moles of the mixture respect to temperature

• dNi_p (float): Thermodynamic derivative of the moles of the species respect to pressure

• dN_p (float): Thermodynamic derivative of the moles of the mixture respect to pressure

• ind (float): List of chemical species indices

• STOP (float): Relative error in moles of species [-]

• STOP_ions (float): Relative error in moles of ionized species [-]

Examples

• N0 = equilibrium_helmholtz_large(self, 0.0716, 3000, self.PS.strR{1}, [])

• [N0, dNi_T, dN_T, dNi_p, dN_p, ind, STOP, STOP_ions] = equilibrium_helmholtz_large(self, 0.0716, 3000, self.PS.strR{1}, [])

1. McBride, Bonnie J. Computer program for calculation of complex chemical equilibrium compositions and applications. Vol. 2. NASA Lewis Research Center, 1996.