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 and McBride [1994].

Thermodynamic derivatives#

All thermodynamic first derivatives can be obtained with any set of three independent first derivatives [McBride, 1996]. 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}

is readily seen that

\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}, [])