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.

EquilibriumSolver#

class EquilibriumSolver(varargin)#

Bases: handle

The EquilibriumSolver() class is used to compute the composition at the equilibrium of multi-component gas mixtures that undergo canonical thermochemical transformations from an initial state (reactants), defined by its initial composition, temperature, and pressure, to a final state (products), defined by a set of chemical species (in gaseous—included ions—or pure condensed phase).

The EquilibriumSolver() object can be initialized as follows:

solver = EquilibriumSolver('problemType', problemType, ...)

Here problemType represents the acronym of the problem to be solved (see below). Additional optional parameters can be provided to customize the solver’s behavior.

Problem types:

  • TP: Equilibrium composition at defined temperature and pressure

  • TV: Equilibrium composition at defined temperature and specific volume

  • HP: Equilibrium composition at defined enthalpy and pressure

  • EV: Equilibrium composition at defined internal energy and specific volume

  • SP: Equilibrium composition at defined entropy and pressure

  • SV: Equilibrium composition at defined entropy and specific volume

See also: Mixture(), solve(), solveArray(), report()

EquilibriumSolver(varargin)#

Constructor

FLAG_CACHE = 'true'#

Flag to clear cache after calculations

FLAG_EOS = 'false'#

Flag to use non-ideal Equation of States (EoS)

FLAG_EXTRAPOLATE = 'true'#

Flag indicating linear extrapolation of the polynomials fits

FLAG_FAST = 'true'#

Flag indicating use guess composition of the previous computation

FLAG_FROZEN = 'false'#

Flag to consider a calorically imperfect gas with frozen chemistry

FLAG_REPORT = 'false'#

Flag to postprocess all the results with predefined plots

FLAG_RESULTS = 'true'#

Flag to print results

FLAG_TCHEM_FROZEN = 'false'#

Flag to consider a thermochemically frozen gas (calorically perfect gas)

FLAG_TIME = 'true'#

Flag to print elapsed time

getGpoint(mix1, mix2, attributeName, x0, molesGuess)#

Get fixed point of a function based on the chemical transformation

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix1 (Mixture) – Properties of the initial mixture

  • mix2 (Mixture) – Properties of the final mixture

  • attributeName (str) – Fieldname in Problem Description (PD)

  • x0 (float) – Guess temperature [K]

Returns:

Tuple containing

  • gpoint (float): Fixed point of the function [J] (HP, EV) or [J/K] (SP, SV)

  • gpointRelative (float): Fixed relative point of the function [J] (HP, EV) or [J/K] (SP, SV)

getPartialDerivative(mix)#

Get value of the partial derivative for the set problem type [J/K] (HP, EV) or [J/K^2] (SP, SV)

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix (Mixture) – Mixture object

Returns:

value (float) – Value of the partial derivative for the set problem type [J/K] (HP, EV) or [J/K^2] (SP, SV)

static getPoint(x_vector, f_vector)#

Get point using the regula falsi method

Parameters:

  • x_vector (float) – Guess temperature [K]

  • f_vector (struct) – evaluated functions [J] (HP, EV) or [J/K] (SP, SV)

Returns:

point (float) – Point of the function [K]

static getPointAitken(x0, g_vector)#

Get fixed point of a function based on the chemical transformation using the Aitken acceleration method

Parameters:

  • x0 (float) – Initial guess for temperature [K]

  • g_vector (struct) – Fixed points of the function [J] (HP, EV) or [J/K] (SP, SV)

Returns:

point (float) – Point of the function [K]

getRatioNewton(mix1, mix2, attributeName, x, molesGuess)#

Get the residual of f, its derivative with temperature, and the relative value of the residual

itMax = '30'#

Max number of iterations - root finding method

itMaxGibbs = '70'#

Max number of iterations - Gibbs/Helmholtz minimization method

itMaxIons = '30'#

Max number of iterations - charge balance (ions)

newton(mix1, mix2, attributeName, x0, molesGuess)#

Find the temperature [K] (root) for the set chemical transformation at equilibrium using the second-order Newton-Raphson method

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix1 (Mixture) – Properties of the initial mixture

  • mix2 (Mixture) – Properties of the final mixture

  • pP (float) – Pressure [bar]

  • attributeName (char) – Attribute name of the problem type

  • x0 (float) – Initial guess for temperature [K]

  • molesGuess (float) – Initial guess for moles in the final mixture

Returns:

Tuple containing

  • x (float): Temperature at equilibrium [K]

  • STOP (float): Relative error [-]

  • molesGuess (float): Updated guess for moles in the final mixture

nsteff(mix1, mix2, attributeName, x0, molesGuess)#

Find the temperature [K] (root) for the set chemical transformation at equilibrium using the third-order Newton-Steffensen method

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix1 (Mixture) – Properties of the initial mixture

  • mix2 (Mixture) – Properties of the final mixture

  • attributeName (char) – Attribute name of the problem type

  • x0 (float) – Initial guess for temperature [K]

  • molesGuess (float) – Initial guess for moles in the final mixture

Returns:

Tuple containing

  • x (float): Temperature at equilibrium [K]

  • STOP (float): Relative error [-]

  • molesGuess (float): Updated guess for moles in the final mixture

plot(mixArray, varargin)#

Plot results

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mixArray (Mixture) – Array of Mixture objects

Optional Args:

  • mixArray_i (Mixture): Array of Mixture objects

Returns:

ax2 (Axes) – Axes object with properties of the mixtures

Examples

  • plot(EquilibriumSolver(), mixArray);

  • plot(EquilibriumSolver(), mixArray, mixArray2);

plotConfig = None#

PlotConfig object

printError(it, T, STOP)#

Print error of the method if the number of iterations is greater than maximum iterations allowed

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • it (float) – Number of iterations executed in the method

  • T (float) – Temperature [K]

  • STOP (float) – Relative error [-]

printTime()#

Print execution time

Parameters:

obj (EquilibriumSolver) – EquilibriumSolver object

problemType = None#

Problem type [TP, TV, HP, EV, SP, SV]

regulaGuess(mix1, mix2, attributeName)#

Find a estimate of the temperature for the set chemical equilibrium transformation using the regula falsi method

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix1 (Mixture) – Properties of the initial mixture

  • mix2 (Mixture) – Properties of the final mixture

  • attributeName (char) – Attribute name of the problem type

Returns:

x0 (float) – Initial guess for temperature [K]

static relaxFactorCondensed(NP, N, psi_j, Delta_nj, indexCondensed, NG, NS, SIZE, tau, RT)#

Compute and apply relaxation factor for condensed species

static relaxFactorGas(NP, nj_gas, Delta_ln_nj, Delta_ln_NP)#

Compute relaxation factor for gas species

static removeElements(NatomE, A0, ind_E, tol)#

Find zero sum elements and remove them from stoichiometrix matrix

Parameters:

  • NatomE (float) – Number of atoms of each element

  • A0 (float) – Stoichiometrix matrix

  • ind_E (float) – Index of electron element

  • tol (float) – Tolerance

Returns:

** A0 (float)* – Stoichiometrix matrix updated * indexRemoveSpecies (float): Index of species to be removed * ind_E (float): Index of electron element updated * NatomE (float): Number of atoms of each element updated

report(mixArray, varargin)#

Postprocess all the results with predefined plots

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mixArray (Mixture) – Array of Mixture objects

Optional args:

  • mixArray_i (Mixture): Array of Mixture objects

Returns:

ax (Axes) – Axes object with properties of the mixtures

Examples

  • ax = report(EquilibriumSolver(), mixArray);

  • ax = report(EquilibriumSolver(), mixArray, mixArray2);

rootMethod = '@newton'#

Root finding method [newton (2nd order), steff (2nd order), or nsteff (3rd order)]

root_T0 = '3000'#

Temperature guess [K] if it’s outside previous range - root finding method

root_T0_l = '1000'#

First temperature guess [K] left branch - root finding method

root_T0_r = '3000'#

First temperature guess [K] right branch - root finding method

slackGuess = '1e-14'#

Initial guess of the slack variables for condensed species

solve(mix, varargin)#

Obtain chemical equilibrium composition and thermodynamic properties

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix (Mixture) – Mixture object containing initial state

Optional Args:

  • mixGuess (Mixture): Mixture object from previous calculation

Returns:

mix (Mixture) – Mixture object updated with equilibrium composition and properties

Examples

  • mix = solve(EquilibriumSolver(), mix);

  • mix = solve(EquilibriumSolver(), mix, mixGuess);

solveArray(mixArray, varargin)#

Obtain chemical equilibrium composition and thermodynamic properties for an array of mixture values

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mixArray (Mixture) – Array of Mixture objects containing initial states

Returns:

mixArray (Mixture) – Array of Mixture objects updated with equilibrium compositions and properties

Example

  • mixArray = solveArray(EquilibriumSolver(), mixArray);

steff(mix1, mix2, attributeName, x0, molesGuess)#

Find the temperature [K] (root) for the set chemical transformation at equilibrium using the Steffenson-Aitken method

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix1 (Mixture) – Properties of the initial mixture

  • mix2 (Mixture) – Properties of the final mixture

  • attributeName (char) – Attribute name of the problem type

  • x0 (float) – Initial guess for temperature [K]

  • molesGuess (float) – Initial guess for moles in the final mixture

Returns:

Tuple containing

  • x (float): Temperature at equilibrium [K]

  • STOP (float): Relative error [-]

  • molesGuess (float): Updated guess for moles in the final mixture

static tempValues(system, NatomE)#

List of indices with nonzero values and lengths

temperatureIons = '0'#

Minimum temperature [K] to consider ionized species

time = None#

Elapsed time [s]

tol0 = '1e-3'#

Tolerance of the root finding algorithm

tolE = '1e-6'#

Tolerance of the mass balance

tolGibbs = '1e-6'#

Tolerance of the Gibbs/Helmholtz minimization method

tolMoles = '1e-14'#

Tolerance of the composition of the mixture

tolMolesGuess = '1e-6'#

Tolerance of the molar composition of the mixture (guess)

tolMultiplierIons = '1e-4'#

Tolerance of the dimensionless Lagrangian multiplier - ions

tolTau = '1e-25'#

Tolerance of the slack variables for condensed species

static updateTemp(N, index, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE)#

Update temporal values

equilibrate(obj, mix2, varargin)#

Obtain properties at equilibrium for the given thermochemical transformation

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix2 (Mixture) – Mixture considering a thermochemical frozen gas

Optional Args:

  • mix2Guess (Mixture): Mixture object of a previous calculation

Returns:

mix2 (Mixture) – Mixture at chemical equilibrium for the given thermochemical transformation

Examples

  • mix2 = equilibrate(EquilibriumSolver(), mix2)

  • mix2 = equilibrate(EquilibriumSolver(), mix2, mix2Guess)

equilibrateT(obj, mix1, mix2, T, varargin)#

Obtain equilibrium properties and composition for the given temperature [K] and pressure [bar] / specific volume [m3/kg]

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix1 (Mixture) – Properties of the initial mixture

  • mix2 (Mixture) – Properties of the final mixture

  • T (float) – Temperature [K]

Optional Args:

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

Returns:

mix2 (Mixture) – Properties of the final mixture

Examples

mix2 = equilibrateT(EquilibriumSolver(), mix1, mix2, 3000) mix2 = equilibrateT(EquilibriumSolver(), mix1, mix2, 3000, molesGuess)

equilibriumCheckCondensed(A0, pi_i, W, indexCondensed, muRT, NC_max, FLAG_ONE, FLAG_RULE)#

Check condensed species that may appear at chemical equilibrium

Parameters:

  • A0 (float) – Stoichiometric matrix

  • pi_i (float) – Dimensionless Lagrange multiplier vector

  • W (float) – Molecular mass [kg/mol] vector

  • indexCondensed (float) – Index condensed species to be considered

  • muRT (float) – List of chemical species indices in gaseous phase

  • NC_max (float) – Maximum number of condensed species (Gibbs phase rule)

  • FLAG_ONE (bool) – Flag indicating to include condensed species in the system one by one

  • FLAG_RULE (bool) – Flag indicating to include condensed species in the system up to the maximum number of condensed species that satisfy the Gibbs phase rule

Returns:

Tuple containing

  • indexCondensed (float): Index condensed species that may appear at chemical equilibrium

  • FLAG_CONDENSED (bool): Flag indicating additional condensed species that may appear at chemical equilibrium

  • dL_dnj (float): Vapor pressure test vector of the species that may appear at chemical equilibrium

equilibriumCheckIons(obj, N, A0, ind_E, indexGas, indexIons)#

Check convergence of ionized species in the mixture

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • N (float) – Mixture composition [mol]

  • A0 (float) – Stoichiometric matrix

  • ind_E (float) – Index of electron element

  • indexGas (float) – List of chemical species indices in gaseous phase

  • indexIons (float) – List of ionized chemical species indices

Returns:

Tuple containing

  • N (float): Mixture composition [mol]

  • STOP (float): Relative error in charge balance [-]

  • FLAG_ION (bool): Flag indicating if ionized species are present in the mixture

equilibriumDerivatives(J, N, A0, NE, indexGas, indexCondensed, indexElements, H0RT)#

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

Parameters:

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

  • N (float) – Mixture composition [mol]

  • A0 (float) – Stoichiometric matrix

  • NE (float) – Temporal total number of elements

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

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

  • indexElements (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

  • 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_T, dN_T, dNi_p, dN_p] = equilibriumDerivatives(J, N, A0, NE, ind, indexGas, indexCondensed, indexElements, H0RT)

equilibriumGibbs(obj, system, p, T, mix, molesGuess)#

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 function is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311 and in Leal, A. M., Kulik, D. A., Kosakowski, G., & Saar, M. O. (2016). Computational methods for reactive transport modeling: An extended law of mass-action, xLMA, method for multiphase equilibrium calculations. Advances in Water Resources, 96, 405-422.

Parameters:

  • obj (EquilibriumSolver) – Equilibrium solver object

  • system (ChemicalSystem) – Chemical system object

  • p (float) – Pressure [bar]

  • T (float) – Temperature [K]

  • mix1 (struct) – Properties of the initial mixture

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

Returns:

Tuple containing

  • N (float): Mixture composition [mol]

  • 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

  • index (float): List of chemical species indices

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

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

  • h0 (float): Molar enthalpy [J/mol]

Examples

  • N = EquilibriumSolver().equilibriumGibbs(system, 1.01325, 3000, mix, [])

  • [N, dNi_T, dN_T, dNi_p, dN_p, index, STOP, STOP_ions, h0] = EquilibriumSolver().equilibriumGibbs(system, 1.01325, 3000, mix, [])

equilibriumGuess(N, NP, A0, muRT, b0, index, indexGas, indexIons, NG, molesGuess)#

Initialize molar composition from a previous calculation or using the simplex algorithm

Parameters:

  • N (float) – Mixture composition [mol]

  • NP (float) – Total moles of gaseous species [mol]

  • A0 (float) – Stoichiometric matrix

  • muRT (float) – Dimensionless chemical potentials

  • b0 (float) – Moles of each element

  • index (float) – Index of chemical species

  • indexGas (float) – Index of gaseous species

  • indexIons (float) – Index of ionized species

  • NG (float) – Number of gaseous species

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

Returns:

Tuple containing

  • N (float): Mixture composition [mol]

  • NP (float): Total moles of gaseous species [mol]

Example

[N, NP] = equilibriumGuess(N, NP, A0, muRT, b0, index, indexGas, indexIons, NG, molesGuess)

equilibriumHelmholtz(obj, system, v, T, mix, molesGuess)#

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 function is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311 and in Leal, A. M., Kulik, D. A., Kosakowski, G., & Saar, M. O. (2016). Computational methods for reactive transport modeling: An extended law of mass-action, xLMA, method for multiphase equilibrium calculations. Advances in Water Resources, 96, 405-422.

Parameters:

  • obj (EquilibriumSolver) – Equilibrium solver object

  • system (ChemicalSystem) – Chemical system object

  • v (float) – Volume [m3]

  • T (float) – Temperature [K]

  • mix1 (Mixture) – Properties of the initial mixture

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

Returns:

Tuple containing

  • N (float): Mixture composition [moles]

  • 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

  • index (float): List of chemical species indices

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

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

  • h0 (float): Molar enthalpy [J/mol]

Examples

  • N = EquilibriumSolver().equilibriumHelmholtz(system, 1, 3000, mix, [])

  • [N, dNi_T, dN_T, dNi_p, dN_p, index, STOP, STOP_ions] = EquilibriumSolver().equilibriumHelmholtz(system, 1, 3000, mix, [])