Skip to content

Solve: chemical equilibrium errors at low temperature when ionized species may be present #1109

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
AlbertoCuadra merged 1 commit into from
Jan 15, 2026
Merged

Solve: chemical equilibrium errors at low temperature when ionized species may be present #1109

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 7 additions & 12 deletions +combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibriumDerivatives.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
function [dNi_T, dN_T, dNi_p, dN_p] = equilibriumDerivatives(J, N, A0, NE, indexGas, indexCondensed, indexElements, H0RT)
function [dNi_T, dN_T, dNi_p, dN_p] = equilibriumDerivatives(J, N, A0, NE, indexGas, indexCondensed, 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
%
Expand All @@ -9,7 +9,6 @@
% 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:
Expand All @@ -21,16 +20,16 @@
% * 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)
% [dNi_T, dN_T, dNi_p, dN_p] = equilibriumDerivatives(J, N, A0, NE, ind, indexGas, indexCondensed, H0RT)

% Equilibrium derivative respect temperature
[dNi_T, dN_T] = equilibrium_dT(J, N, A0, NE, indexGas, indexCondensed, indexElements, H0RT);
[dNi_T, dN_T] = equilibrium_dT(J, N, A0, NE, indexGas, indexCondensed, H0RT);

% Equilibrium derivative respect pressure
[dNi_p, dN_p] = equilibrium_dp(J, N, A0, NE, indexGas, indexCondensed, indexElements);
[dNi_p, dN_p] = equilibrium_dp(J, N, A0, NE, indexGas, indexCondensed);
end

function [dNi_T, dN_T] = equilibrium_dT(J, N, A0, NE, indexGas, indexCondensed, indexElements, H0RT)
function [dNi_T, dN_T] = equilibrium_dT(J, N, A0, NE, indexGas, indexCondensed, 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
%
Expand All @@ -41,7 +40,6 @@
% 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:
Expand All @@ -51,11 +49,10 @@
% * dN_T (float): Thermodynamic derivative of the moles of the mixture respect to temperature
%
% Example:
% [dNi_T, dN_T] = equilibrium_dT(J, N, A0, NE, ind, indexGas, indexCondensed, indexElements, H0RT)
% [dNi_T, dN_T] = equilibrium_dT(J, N, A0, NE, ind, indexGas, indexCondensed, H0RT)

% Definitions
opts.SYM = true; % Options linsolve method: real symmetric
A0 = A0(:, indexElements);

% Initialization
dNi_T = zeros(length(N), 1);
Expand Down Expand Up @@ -83,7 +80,7 @@

end

function [dNi_p, dN_p] = equilibrium_dp(J, N, A0, NE, indexGas, indexCondensed, indexElements)
function [dNi_p, dN_p] = equilibrium_dp(J, N, A0, NE, indexGas, indexCondensed)
% 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
%
Expand All @@ -94,7 +91,6 @@
% 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
%
% Returns:
% Tuple containing
Expand All @@ -107,7 +103,6 @@

% Definitions
opts.SYM = true; % Options linsolve method: real symmetric
A0 = A0(:, indexElements);

% Initialization
dNi_p = zeros(length(N), 1);
Expand Down
48 changes: 41 additions & 7 deletions +combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibriumGibbs.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,10 +64,16 @@
% Find indeces of the species/elements that we have to remove from the stoichiometric matrix A0
% for the sum of elements whose value is <= tolMoles
[A0, indexRemoveSpecies, ind_E, NatomE] = obj.removeElements(NatomE, A0, system.ind_E, obj.tolMoles);

% Check if element E (electron) is present
FLAG_E = ~isempty(ind_E);

% List of indices with nonzero values
[index, indexGas, indexCondensed, indexIons, indexElements, NE, NG, NS] = obj.tempValues(system, NatomE);

% Remove elements with zero atoms from the stoichiometric matrix A0
A0 = A0(:, indexElements);

% Update temp values
if ~isempty(indexRemoveSpecies)
[index, indexCondensed, indexGas, indexIons, NG, NS] = obj.updateTemp(N, indexRemoveSpecies, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE);
Expand Down Expand Up @@ -105,7 +111,7 @@
A0_T = A0';

% Initialize composition vector N
[N, NP] = obj.equilibriumGuess(N, NP, A0_T(indexElements, index0), muRT(index0), NatomE, index0, indexGas_0, indexIons, NG, molesGuess);
[N, NP] = obj.equilibriumGuess(N, NP, A0_T(:, index0), muRT(index0), NatomE, index0, indexGas_0, indexIons, NG, molesGuess);

% Initialization
psi_j = system.propertyVector;
Expand All @@ -118,7 +124,7 @@
x = equilibriumLoopCondensed(x);

% Update matrix J (jacobian) to compute the thermodynamic derivatives
J = update_matrix_J(A0_T(indexElements, :), J22, N, NP, indexGas, indexCondensed, NS - NG, psi_j);
J = update_matrix_J(A0_T, J22, N, NP, indexGas, indexCondensed, NS - NG, psi_j);
J(end, end) = 0;

% Molar enthalpy [J/mol]
Expand All @@ -128,7 +134,7 @@
H0RT = h0 / RT;

% Compute thermodynamic derivates
[dNi_T, dN_T, dNi_p, dN_p] = obj.equilibriumDerivatives(J, N, A0, NE, indexGas, indexCondensed, indexElements, H0RT);
[dNi_T, dN_T, dNi_p, dN_p] = obj.equilibriumDerivatives(J, N, A0, NE, indexGas, indexCondensed, H0RT);

% NESTED FUNCTION
function x = equilibriumLoop
Expand Down Expand Up @@ -157,6 +163,18 @@
% Check singular matrix
if any(isnan(x) | isinf(x))

% Check if index of electron element is defined
if FLAG_E

% Check residual of charge balance
if max( [norm(J(ind_E, :), 1), norm(J(:, ind_E), 1), abs(b(ind_E))] ) < obj.tolE
% Remove element E from matrix
removeElementElectron();
continue
end

end

% Update temp indeces
indexGas = indexGas_0;
indexCondensed = indexCondensed_0;
Expand Down Expand Up @@ -242,8 +260,10 @@
end

% Remove element E from matrix
indexElements(ind_E) = [];
NE = NE - 1;
removeElementElectron();

% Recompute chemical equilibrium without ions
x = equilibriumLoop;
end

function x = equilibriumLoopCondensed(x)
Expand Down Expand Up @@ -281,7 +301,7 @@
end

% Check condensed species
[indexCondensed_add, FLAG_CONDENSED, ~] = obj.equilibriumCheckCondensed(A0(:, indexElements), x(1:NE), W(indexCondensed_check), indexCondensed_check, muRT, NC_max, FLAG_ONE, FLAG_RULE);
[indexCondensed_add, FLAG_CONDENSED, ~] = obj.equilibriumCheckCondensed(A0, x(1:NE), W(indexCondensed_check), indexCondensed_check, muRT, NC_max, FLAG_ONE, FLAG_RULE);

if ~FLAG_CONDENSED
break
Expand Down Expand Up @@ -357,7 +377,7 @@
end

% Check if there were species not considered
[~, FLAG_CONDENSED, dL_dnj] = obj.equilibriumCheckCondensed(A0(:, indexElements), x(1:NE), W(indexCondensed_0), indexCondensed_0, muRT, NC_max, FLAG_ONE, FLAG_RULE);
[~, FLAG_CONDENSED, dL_dnj] = obj.equilibriumCheckCondensed(A0, x(1:NE), W(indexCondensed_0), indexCondensed_0, muRT, NC_max, FLAG_ONE, FLAG_RULE);

% Recompute if there are condensed species that may appear at chemical equilibrium
if FLAG_CONDENSED && any(abs(dL_dnj) > 1e-4)
Expand All @@ -366,6 +386,20 @@

end

function removeElementElectron()
% Remove element E from matrix
A0(:, ind_E) = []; A0_T(ind_E, :) = [];
indexIons = []; indexElements(ind_E) = [];
NatomE(ind_E) = [];
NE = NE - 1;

% Update FLAG_E
FLAG_E = false;

% Update indeces
[~, indexCondensed, indexGas, indexIons, NG, NS] = obj.updateTemp(N, index, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE);
end

end

% SUB-PASS FUNCTIONS
Expand Down
52 changes: 43 additions & 9 deletions +combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibriumHelmholtz.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,15 @@
% for the sum of elements whose value is <= tolMoles
[A0, indexRemoveSpecies, ind_E, NatomE] = obj.removeElements(NatomE, A0, system.ind_E, obj.tolMoles);

% Check if element E (electron) is present
FLAG_E = ~isempty(ind_E);

% List of indices with nonzero values
[index, indexGas, indexCondensed, indexIons, indexElements, NE, NG, NS] = obj.tempValues(system, NatomE);

% Remove elements with zero atoms from the stoichiometric matrix A0
A0 = A0(:, indexElements);

% Update temp values
if ~isempty(indexRemoveSpecies)
[index, indexCondensed, indexGas, indexIons, NG, NS] = obj.updateTemp(N, indexRemoveSpecies, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE);
Expand Down Expand Up @@ -104,7 +110,7 @@
A0_T = A0';

% Initialize composition matrix N [mol, FLAG_CONDENSED]
[N, NP] = obj.equilibriumGuess(N, NP, A0_T(indexElements, index0), muRT(index0), NatomE, index0, indexGas_0, indexIons, NG, molesGuess);
[N, NP] = obj.equilibriumGuess(N, NP, A0_T(:, index0), muRT(index0), NatomE, index0, indexGas_0, indexIons, NG, molesGuess);

% Initialization
psi_j = system.propertyVector;
Expand All @@ -117,9 +123,9 @@
x = equilibriumLoopCondensed(x);

% Update matrix J (jacobian) to compute the thermodynamic derivatives
J = update_matrix_J(A0_T(indexElements, :), N, indexGas, indexCondensed, psi_j);
J = update_matrix_J(A0_T, N, indexGas, indexCondensed, psi_j);
temp_zero = zeros(NS - NG + 1, 1);
J12_2 = [sum(A0_T(indexElements, indexGas) .* N(indexGas)', 2); temp_zero(1:end-1)];
J12_2 = [sum(A0_T(:, indexGas) .* N(indexGas)', 2); temp_zero(1:end-1)];
J = [J, J12_2; J12_2', 0];

% Molar enthalpy [J/mol]
Expand All @@ -129,7 +135,7 @@
H0RT = h0 / RT;

% Compute thermodynamic derivates
[dNi_T, dN_T, dNi_p, dN_p] = obj.equilibriumDerivatives(J, N, A0, NE, indexGas, indexCondensed, indexElements, H0RT);
[dNi_T, dN_T, dNi_p, dN_p] = obj.equilibriumDerivatives(J, N, A0, NE, indexGas, indexCondensed, H0RT);

% NESTED FUNCTION
function x = equilibriumLoop
Expand Down Expand Up @@ -161,6 +167,18 @@
% Check singular matrix
if any(isnan(x) | isinf(x))

% Check if index of electron element is defined
if FLAG_E

% Check residual of charge balance
if max( [norm(J(ind_E, :), 1), norm(J(:, ind_E), 1), abs(b(ind_E))] ) < obj.tolE
% Remove element E from matrix
removeElementElectron();
continue
end

end

% Update temp indeces
indexGas = indexGas_0;
indexCondensed = indexCondensed_0;
Expand Down Expand Up @@ -236,8 +254,10 @@
end

% Remove element E from matrix
indexElements(ind_E) = [];
NE = NE - 1;
removeElementElectron();

% Recompute chemical equilibrium without ions
x = equilibriumLoop;
end

function x = equilibriumLoopCondensed(x)
Expand Down Expand Up @@ -275,7 +295,7 @@
end

% Check condensed species
[indexCondensed_add, FLAG_CONDENSED, ~] = obj.equilibriumCheckCondensed(A0(:, indexElements), x(1:NE), W(indexCondensed_check), indexCondensed_check, muRT, NC_max, FLAG_ONE, FLAG_RULE);
[indexCondensed_add, FLAG_CONDENSED, ~] = obj.equilibriumCheckCondensed(A0, x(1:NE), W(indexCondensed_check), indexCondensed_check, muRT, NC_max, FLAG_ONE, FLAG_RULE);

if ~FLAG_CONDENSED
break
Expand Down Expand Up @@ -341,7 +361,7 @@
end

% Check if there were species not considered
[~, FLAG_CONDENSED, dL_dnj] = obj.equilibriumCheckCondensed(A0(:, indexElements), x(1:NE), W(indexCondensed_0), indexCondensed_0, muRT, NC_max, FLAG_ONE, FLAG_RULE);
[~, FLAG_CONDENSED, dL_dnj] = obj.equilibriumCheckCondensed(A0, x(1:NE), W(indexCondensed_0), indexCondensed_0, muRT, NC_max, FLAG_ONE, FLAG_RULE);

% Recompute if there are condensed species that may appear at chemical equilibrium
if FLAG_CONDENSED && any(abs(dL_dnj) > 1e-4)
Expand All @@ -350,6 +370,20 @@

end

function removeElementElectron()
% Remove element E from matrix
A0(:, ind_E) = []; A0_T(ind_E, :) = [];
indexIons = []; indexElements(ind_E) = [];
NatomE(ind_E) = [];
NE = NE - 1;

% Update FLAG_E
FLAG_E = false;

% Update indeces
[~, indexCondensed, indexGas, indexIons, NG, NS] = obj.updateTemp(N, index, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE);
end

end

% SUB-PASS FUNCTIONS
Expand Down Expand Up @@ -390,7 +424,7 @@
% Compute vector b
bi = N(index)' * A0(index, :);

if any(indexIons)
if ~isempty(indexIons)
bi(ind_E) = NatomE(ind_E);
end

Expand Down