Skip to content

Commit bd6a46b

Browse files
AlbertoCuadraAlbertoCuadra
authored
Merge pull request #1068 from CombustionToolbox/develop
Solve: error in `DetonationCJ` method within `DetonationSolver` class when guess provides solutions below tolerance
2 parents 35a05de + 2d28159 commit bd6a46b

7 files changed

Lines changed: 45 additions & 58 deletions

File tree

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

Lines changed: 0 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -6,9 +6,6 @@
66
% * COMPLETE
77
% * HC/O2/N2 EXTENDED
88
% * SOOT FORMATION EXTENDED
9-
% * NASA ALL
10-
% * NASA ALL CONDENSED
11-
% * NASA ALL IONS
129
% * AIR, DISSOCIATED AIR
1310
% * AIR IONS, AIR_IONS
1411
% * IDEAL_AIR, AIR_IDEAL
@@ -146,33 +143,6 @@
146143
'C6H5O_phenoxy', 'C6H5_phenyl', 'C7H7_benzyl', ...
147144
'C7H8', 'C8H8_styrene', 'C10H8_naphthale'};
148145

149-
case 'NASA ALL'
150-
listSpecies = find_species_LS(fieldnames(database.species), ...
151-
{'C', 'N', 'O', 'minus', 'plus', 'Ar', 'H'}, 'any', ...
152-
{'I', 'S', 'L', 'T', 'P', 'F', 'ab', 'W', 'AL','He' ...
153-
'Z', 'X', 'R', 'Os', 'Cr', 'Br', 'G', 'K','Li' ...
154-
'U', 'Co', 'Cu', 'B', 'V', 'Ni', 'Na', 'Mg','Hg' ...
155-
'Mo', 'Ag', 'Nb', 'Cb', 'CL', 'D', 'T', 'M','minus' ...
156-
'Ca', 'Cs', 'Ne', 'Cd', 'Mn', 'cr', 'plus'}, 'all');
157-
158-
case {'NASA ALL CONDENSED', 'NASA ALL_CONDENSED', 'NASA_ALL_CONDENSED'}
159-
listSpecies = find_species_LS(fieldnames(database.species), ...
160-
{'C', 'N', 'O', 'Ar', 'H'}, 'any', ...
161-
{'I', 'S', 'T', 'P', 'F', 'ab', 'W', 'AL','He' ...
162-
'Z', 'X', 'R', 'Os', 'Cr', 'Br', 'G', 'K','Li' ...
163-
'U', 'Co', 'Cu', 'B', 'V', 'Ni', 'Na', 'Mg','Hg' ...
164-
'Mo', 'Ag', 'Nb', 'Cb', 'CL', 'D', 'T', 'M', 'minus', ...
165-
'Ca', 'Cs', 'Ne', 'Cd', 'Mn', 'plus'}, 'all');
166-
167-
case {'NASA ALL IONS', 'NASA ALL_IONS', 'NASA_ALL_IONS'}
168-
listSpecies = find_species_LS(fieldnames(database.species), ...
169-
{'C', 'N', 'O', 'minus', 'plus', 'Ar', 'H'}, 'any', ...
170-
{'I', 'S', 'L', 'T', 'P', 'F', 'ab', 'W', 'AL','He' ...
171-
'Z', 'X', 'R', 'Os', 'Cr', 'Br', 'G', 'K','Li' ...
172-
'U', 'Co', 'Cu', 'B', 'V', 'Ni', 'Na', 'Mg','Hg' ...
173-
'Mo', 'Ag', 'Nb', 'Cb', 'CL', 'D', 'T', 'M', ...
174-
'Ca', 'Cs', 'Ne', 'Cd', 'Mn', 'cr', 'cyclo'}, 'all');
175-
176146
case {'AIR', 'DISSOCIATED AIR'}
177147
listSpecies = {'CO2', 'CO', 'O2', 'N2', 'Ar', 'O', 'O3', ...
178148
'N', 'NO', 'NO2', 'NO3', 'N2O', 'N2O3', ...

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

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,10 @@
2424
% Example:
2525
% P = getPressure(obj, 300, 0.024)
2626

27+
% Definitions
2728
R0 = combustiontoolbox.common.Constants.R0; % Universal gas constant [J/(K mol)]
29+
30+
% Compute pressure [Pa]
2831
pressure = (R0 .* temperature) ./ molarVolume;
2932
end
3033

@@ -41,7 +44,10 @@
4144
% Example:
4245
% V = getVolume(obj, 300, 1e5)
4346

47+
% Definitions
4448
R0 = combustiontoolbox.common.Constants.R0; % Universal gas constant [J/(K mol)]
49+
50+
% Compute molar volume [m3/mol]
4551
molarVolume = R0 * temperature ./ pressure;
4652
end
4753

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

Lines changed: 6 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,7 @@
3939
equilibrateT(obj, mix1, mix2, T, molesGuess);
4040

4141
% Check convergence in case the problemType is TP (defined Temperature and Pressure)
42-
print_convergence(mix2.errorMoles, obj.tolGibbs, mix2.errorMolesIons, obj.tolMultiplierIons, obj.problemType)
42+
checkConvergence(mix2.errorMoles, obj.tolGibbs, mix2.errorMolesIons, obj.tolMultiplierIons, obj.problemType)
4343

4444
% Save error from root finding algorithm
4545
mix2.errorProblem = STOP;
@@ -92,31 +92,19 @@
9292
[x, STOP, molesGuess] = obj.rootMethod(obj, mix1, mix2, attributeName, x0, molesGuess);
9393
end
9494

95-
function print_convergence(STOP, TOL, STOP_ions, TOL_ions, ProblemType)
96-
% Print tolerance error if the convergence criteria was not satisfied
95+
function checkConvergence(STOP, TOL, STOP_ions, TOL_ions, problemType)
96+
% Check tolerance error if the convergence criteria was not satisfied
9797

98-
if ~strcmpi(ProblemType, 'TP')
98+
if ~strcmpi(problemType, 'TP')
9999
return
100100
end
101101

102102
if STOP > TOL
103-
fprintf('***********************************************************\n')
104-
fprintf('Convergence error number of moles: %1.2e\n', STOP);
103+
warning('Convergence error number of moles: %1.2e\n', STOP);
105104
end
106105

107106
if STOP_ions > TOL_ions
108-
fprintf('***********************************************************\n')
109-
fprintf('Convergence error in charge balance: %1.2e\n', STOP_ions);
107+
warning('Convergence error in charge balance: %1.2e\n', STOP_ions);
110108
end
111109

112-
end
113-
114-
function mix = set_volume_SV(obj, mix)
115-
% If the problem type is SV, the product's volume is based on the given v_P/v_R ratio
116-
if ~strcmpi(obj.problemType, 'SV')
117-
return
118-
end
119-
120-
%mix.v = mix.v * obj.PD.vP_vR.value;
121-
fprintf('\nto be clarified\n')
122110
end

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

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -94,8 +94,11 @@
9494
% Molar Gibbs energy [J/mol]
9595
g0([indexGas_0, indexCondensed_0]) = set_g0(system.listSpecies([indexGas_0, indexCondensed_0]), T, system.species);
9696

97+
% Dimensionless Gibbs energy
98+
g0RT = g0/RT;
99+
97100
% Dimensionless chemical potential
98-
muRT = g0/RT;
101+
muRT = g0RT;
99102

100103
% Construction of part of matrix J
101104
J22 = zeros(NS - NG + 1);
@@ -140,7 +143,7 @@
140143
while STOP > obj.tolGibbs && it < itMax
141144
it = it + 1;
142145
% Chemical potentials
143-
muRT(indexGas) = g0(indexGas) / RT + log(N(indexGas) / NP) + log(p);
146+
muRT(indexGas) = g0RT(indexGas) + log(N(indexGas) / NP) + log(p);
144147

145148
% Construction of matrix J
146149
J = update_matrix_J(A0_T, J22, N, NP, indexGas, indexCondensed, NS - NG, psi_j);
@@ -419,5 +422,5 @@
419422

420423
function Delta_ln_nj = update_Delta_ln_nj(A0, pi_i, Delta_NP, muRT, indexGas)
421424
% Compute correction moles of gases
422-
Delta_ln_nj = sum(A0(indexGas, :)' .* pi_i, 1)' + Delta_NP - muRT(indexGas);
425+
Delta_ln_nj = A0(indexGas, :) * pi_i + Delta_NP - muRT(indexGas);
423426
end

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

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -94,8 +94,11 @@
9494
% Molar Gibbs energy [J/mol]
9595
g0([indexGas_0, indexCondensed_0]) = set_g0(system.listSpecies([indexGas_0, indexCondensed_0]), T, system.species);
9696

97+
% Dimensionless Gibbs energy
98+
g0RT = g0/RT;
99+
97100
% Dimensionless chemical potential
98-
muRT = g0/RT;
101+
muRT = g0RT;
99102

100103
% Construction of part of matrix J
101104
A0_T = A0';
@@ -141,7 +144,7 @@
141144
while STOP > obj.tolGibbs && it < itMax
142145
it = it + 1;
143146
% Chemical potential
144-
muRT(indexGas) = g0(indexGas) / RT + log(N(indexGas) * RT / v * 1e-5);
147+
muRT(indexGas) = g0RT(indexGas) + log(N(indexGas) * RT / v * 1e-5);
145148

146149
% Compute total number of moles
147150
NP = sum(N(indexGas));
@@ -393,5 +396,5 @@
393396

394397
function Delta_ln_nj = update_Delta_ln_nj(A0, pi_i, muRT, indexGas)
395398
% Compute correction moles of gases
396-
Delta_ln_nj = sum(A0(indexGas, :)' .* pi_i, 1)' - muRT(indexGas);
399+
Delta_ln_nj = A0(indexGas, :) * pi_i - muRT(indexGas);
397400
end

‎+combustiontoolbox/+shockdetonation/@DetonationSolver/detonationCJ.m‎

Lines changed: 20 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -45,35 +45,52 @@
4545

4646
% Save state
4747
[mix1, mix2] = save_state(mix1, mix2, STOP);
48-
% mix2.T_guess = T_guess;
49-
% mix2.p_guess = p2_guess;
5048

5149
% NESTED-FUNCTIONS
5250
function [T2, p2, STOP, it, T_guess, p2_guess] = solve_cj_detonation(FLAG_FAST)
5351
% Miscellaneous
5452
it = 0; itMax = obj.itMax;
53+
5554
% Initial estimates of p2/p1 and T2/T1
5655
[p2, T2, p2p1, T2T1, STOP] = get_guess(obj, mix1, mix2);
5756
T_guess = T2; p2_guess = p2;
57+
58+
% Check STOP
59+
if STOP < obj.tol0
60+
% Set pressure and temperature of mix2
61+
mix2.p = p2; mix2.T = T2;
62+
63+
% Calculate state given T & p
64+
mix2 = state(obj.equilibriumSolver, mix1, mix2, []);
65+
end
66+
5867
% Check FLAG
5968
if ~FLAG_FAST, guess_moles = []; end
69+
6070
% Loop
6171
while STOP > obj.tol0 && it < itMax
6272
% Update iteration
6373
it = it + 1;
74+
6475
% Construction of the Jacobian matrix and vector b
6576
[J, b, guess_moles] = update_system(obj.equilibriumSolver, mix1, mix2, p2, T2, R0, guess_moles, FLAG_FAST);
77+
6678
% Solve of the linear system J*x = b
6779
x = linsolve(J, b);
80+
6881
% Calculate correction factor
6982
lambda = relax_factor(x);
83+
7084
% Apply correction
7185
[log_p2p1, log_T2T1] = apply_correction(x, p2p1, T2T1, lambda);
86+
7287
% Apply antilog
7388
[p2, T2] = apply_antilog(mix1, log_p2p1, log_T2T1); % [bar] and [K]
89+
7490
% Update ratios
7591
p2p1 = p2 / mix1.p;
7692
T2T1 = T2 / mix1.T;
93+
7794
% Compute STOP criteria
7895
STOP = compute_STOP(x);
7996
end
@@ -131,7 +148,7 @@ function print_guess(T2T1, p2p1, T1, p1, Q)
131148
% Set pressure and temperature of mix2
132149
mix2.p = p2; mix2.T = T2;
133150

134-
% Calculate frozen state given T & p
151+
% Calculate state given T & p
135152
[mix2, r2, dVdT_p, dVdp_T] = state(equilibriumSolver, mix1, mix2, guess_moles);
136153

137154
r2r1 = r2 / r1; % [-]

‎examples/Example_DET_POLAR.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323
DB = NasaDatabase();
2424

2525
% Define chemical system
26-
system = ChemicalSystem(DB, 'hydrogen');
26+
system = ChemicalSystem(DB);
2727

2828
% Initialize mixture
2929
mix = Mixture(system);

0 commit comments

Comments
 (0)