Skip to content

Commit 40d3b54

Browse files
AlbertoCuadraAlbertoCuadra
authored
Merge pull request #1123 from CombustionToolbox/develop
Solve: error updating thermo state with chemical equilibrium mixtures and error initial guess with ionized species at very low temperature
2 parents 1aae625 + 221223c commit 40d3b54

5 files changed

Lines changed: 95 additions & 46 deletions

File tree

‎+combustiontoolbox/+core/@ChemicalSystem/ChemicalSystem.m‎

Lines changed: 25 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -401,13 +401,37 @@
401401
obj = updatePropertiesMatrixCompositionFast(obj, moles, index);
402402
end
403403

404+
function obj = setPropertiesMatrixThermo(obj, T, index)
405+
% Update temperature-dependent properties without changing composition
406+
%
407+
% Args:
408+
% obj (ChemicalSystem): ChemicalSystem object
409+
% T (float): Temperature [K]
410+
% index (float): Vector with the indexes of the species to update
411+
%
412+
% Returns:
413+
% obj (ChemicalSystem): ChemicalSystem object with updated thermodynamic properties
414+
415+
if nargin < 3
416+
index = find(obj.propertiesMatrix(:, obj.ind_ni));
417+
end
418+
419+
obj.propertiesMatrix(:, obj.ind_hi:obj.numProperties) = 0;
420+
421+
if isempty(index)
422+
return
423+
end
424+
425+
obj = updatePropertiesMatrixThermo(obj, obj.listSpecies(index), T, index, true);
426+
end
427+
404428
function obj = clean(obj)
405429
% Set temperature-dependent matrix properties to zero
406430
obj.propertiesMatrix(:, 5:end) = 0;
407431
end
408432

409433
function obj = cleanMoles(obj)
410-
% Set temperature-dependent matrix properties to zero
434+
% Set moles vector to zero
411435
obj.propertiesMatrix(:, obj.ind_ni) = 0;
412436
end
413437

‎+combustiontoolbox/+core/@Mixture/Mixture.m‎

Lines changed: 14 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1021,21 +1021,24 @@
10211021
% obj (Mixture): Mixture object with updated properties
10221022

10231023
% Check if initial state is defined (temperature, pressure/volume, and composition)
1024-
if ~sum(obj.quantity) && ~obj.T || (~obj.p && ~obj.vSpecific)
1024+
if ~obj.T || (~obj.p && ~obj.vSpecific)
10251025
return
10261026
end
10271027

1028-
if obj.FLAG_VOLUME
1029-
% Compute molar volume [m3/mol] from specific volume [m3/kg]
1030-
vMolar = vSpecific2vMolar(obj, obj.vSpecific, obj.quantity, obj.quantity(obj.indexGas));
1031-
% Compute pressure in Pascals [Pa] using the equationState
1032-
pressure = obj.equationState.getPressure(obj.T, vMolar, obj.quantity / sum(obj.quantity), obj.chemicalSystem);
1033-
% Convert pressure to [bar]
1034-
obj.p = pressure * combustiontoolbox.common.Units.Pa2bar;
1028+
% Initialize from the recipe only if no state composition exists yet
1029+
system = obj.chemicalSystem;
1030+
currentMoles = system.propertiesMatrix(:, system.ind_ni);
1031+
1032+
if sum(currentMoles) > 0
1033+
system.setPropertiesMatrixThermo(obj.T, find(currentMoles));
1034+
elseif sum(obj.quantity)
1035+
system.setPropertiesMatrixInitialIndex(obj.listSpecies, obj.quantity, obj.T, obj.indexSpecies);
1036+
else
1037+
return
10351038
end
1036-
1037-
% Assign values to the propertiesMatrix
1038-
obj.chemicalSystem.setPropertiesMatrixInitialIndex(obj.listSpecies, obj.quantity, obj.T, obj.indexSpecies);
1039+
1040+
% Refresh composition properties from the state composition
1041+
computeComposition(obj);
10391042

10401043
% Compute thermodynamic properties
10411044
computeThermodynamics(obj);

‎+combustiontoolbox/+equilibrium/@EquilibriumSolver/EquilibriumSolver.m‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@
3535
itMaxIons = 30 % Max number of iterations - charge balance (ions)
3636
itMaxRecursion = 30 % Max number of recursion iterations - e.g., re-check condendensed species
3737
slackGuess = 1e-14 % Initial guess of the slack variables for condensed species
38-
temperatureIons = 0 % Minimum temperature [K] to consider ionized species
38+
temperatureIons = 1500 % Minimum temperature [K] to consider ionized species
3939
tol0 = 1e-3 % Tolerance of the root finding algorithm
4040
itMax = 30 % Max number of iterations - root finding method
4141
rootMethod = @newton % Root finding method [newton (2nd order), steff (2nd order), or nsteff (3rd order)]
@@ -415,7 +415,7 @@ function printTime(obj)
415415
NG = length(indexGas);
416416
NS = length(index);
417417
end
418-
418+
419419
function [index, indexCondensed, indexGas, indexIons, NG, NS, N] = updateTemp(N, index, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE)
420420
% Update temporal values
421421

‎+combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibriumGibbs.m‎

Lines changed: 27 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -73,11 +73,17 @@
7373

7474
% Remove elements with zero atoms from the stoichiometric matrix A0
7575
A0 = A0(:, indexElements);
76+
A0_T = A0';
7677

7778
% Update temp values
7879
if ~isempty(indexRemoveSpecies)
7980
[index, indexCondensed, indexGas, indexIons, NG, NS] = obj.updateTemp(N, indexRemoveSpecies, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE);
8081
end
82+
83+
% Remove ionized species below the configured temperature threshold
84+
if FLAG_E && ~isempty(indexIons) && T < obj.temperatureIons
85+
removeIons();
86+
end
8187

8288
% Remove condensed species with temperature out of bounds
8389
indexCondensed = system.filterSpeciesTemperatureRange(T, indexCondensed, NS - NG, false);
@@ -108,7 +114,6 @@
108114

109115
% Construction of part of matrix J
110116
J22 = zeros(NS - NG + 1);
111-
A0_T = A0';
112117

113118
% Initialize composition vector N
114119
[N, NP] = obj.equilibriumGuess(N, NP, A0_T(:, index0), muRT(index0), NatomE, index0, indexGas_0, indexIons, NG, molesGuess);
@@ -168,8 +173,8 @@
168173

169174
% Check residual of charge balance
170175
if max( [norm(J(ind_E, :), 1), norm(J(:, ind_E), 1), abs(b(ind_E))] ) < obj.tolE
171-
% Remove element E from matrix
172-
removeElementElectron();
176+
% Remove ionized species and element E from matrix
177+
removeIons();
173178
continue
174179
end
175180

@@ -259,8 +264,8 @@
259264
return
260265
end
261266

262-
% Remove element E from matrix
263-
removeElementElectron();
267+
% Remove ionized species and element E from matrix
268+
removeIons();
264269

265270
% Recompute chemical equilibrium without ions
266271
x = equilibriumLoop;
@@ -386,18 +391,25 @@
386391

387392
end
388393

389-
function removeElementElectron()
390-
% Remove element E from matrix
391-
A0(:, ind_E) = []; A0_T(ind_E, :) = [];
392-
indexIons = []; indexElements(ind_E) = [];
393-
NatomE(ind_E) = [];
394-
NE = NE - 1;
394+
function removeIons()
395+
% Remove ionized species and the electron element from the active set
396+
if ~isempty(indexIons)
397+
N(indexIons) = 0;
398+
indexGas(ismember(indexGas, indexIons)) = [];
399+
indexIons = [];
400+
end
395401

396-
% Update FLAG_E
397-
FLAG_E = false;
402+
if FLAG_E
403+
A0(:, ind_E) = []; A0_T(ind_E, :) = [];
404+
indexElements(ind_E) = [];
405+
NatomE(ind_E) = [];
406+
NE = NE - 1;
407+
FLAG_E = false;
408+
end
398409

399-
% Update indeces
400-
[~, indexCondensed, indexGas, indexIons, NG, NS] = obj.updateTemp(N, index, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE);
410+
index = [indexGas, indexCondensed];
411+
NG = length(indexGas);
412+
NS = length(index);
401413
end
402414

403415
end

‎+combustiontoolbox/+equilibrium/@EquilibriumSolver/equilibriumHelmholtz.m‎

Lines changed: 27 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -73,12 +73,18 @@
7373

7474
% Remove elements with zero atoms from the stoichiometric matrix A0
7575
A0 = A0(:, indexElements);
76+
A0_T = A0';
7677

7778
% Update temp values
7879
if ~isempty(indexRemoveSpecies)
7980
[index, indexCondensed, indexGas, indexIons, NG, NS] = obj.updateTemp(N, indexRemoveSpecies, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE);
8081
end
8182

83+
% Remove ionized species below the configured temperature threshold
84+
if FLAG_E && ~isempty(indexIons) && T < obj.temperatureIons
85+
removeIons();
86+
end
87+
8288
% Remove condensed species with temperature out of bounds
8389
indexCondensed = system.filterSpeciesTemperatureRange(T, indexCondensed, NS - NG, false);
8490

@@ -106,9 +112,6 @@
106112
% Dimensionless chemical potential
107113
muRT = g0RT;
108114

109-
% Construction of part of matrix J
110-
A0_T = A0';
111-
112115
% Initialize composition matrix N [mol, FLAG_CONDENSED]
113116
[N, NP] = obj.equilibriumGuess(N, NP, A0_T(:, index0), muRT(index0), NatomE, index0, indexGas_0, indexIons, NG, molesGuess);
114117

@@ -172,8 +175,8 @@
172175

173176
% Check residual of charge balance
174177
if max( [norm(J(ind_E, :), 1), norm(J(:, ind_E), 1), abs(b(ind_E))] ) < obj.tolE
175-
% Remove element E from matrix
176-
removeElementElectron();
178+
% Remove ionized species and element E from matrix
179+
removeIons();
177180
continue
178181
end
179182

@@ -253,8 +256,8 @@
253256
return
254257
end
255258

256-
% Remove element E from matrix
257-
removeElementElectron();
259+
% Remove ionized species and element E from matrix
260+
removeIons();
258261

259262
% Recompute chemical equilibrium without ions
260263
x = equilibriumLoop;
@@ -370,18 +373,25 @@
370373

371374
end
372375

373-
function removeElementElectron()
374-
% Remove element E from matrix
375-
A0(:, ind_E) = []; A0_T(ind_E, :) = [];
376-
indexIons = []; indexElements(ind_E) = [];
377-
NatomE(ind_E) = [];
378-
NE = NE - 1;
376+
function removeIons()
377+
% Remove ionized species and the electron element from the active set
378+
if ~isempty(indexIons)
379+
N(indexIons) = 0;
380+
indexGas(ismember(indexGas, indexIons)) = [];
381+
indexIons = [];
382+
end
379383

380-
% Update FLAG_E
381-
FLAG_E = false;
384+
if FLAG_E
385+
A0(:, ind_E) = []; A0_T(ind_E, :) = [];
386+
indexElements(ind_E) = [];
387+
NatomE(ind_E) = [];
388+
NE = NE - 1;
389+
FLAG_E = false;
390+
end
382391

383-
% Update indeces
384-
[~, indexCondensed, indexGas, indexIons, NG, NS] = obj.updateTemp(N, index, indexCondensed, indexGas, indexIons, NP, NG, NS, SIZE);
392+
index = [indexGas, indexCondensed];
393+
NG = length(indexGas);
394+
NS = length(index);
385395
end
386396

387397
end

0 commit comments

Comments
 (0)