Skip to content

Commit 789e1ae

Browse files
AlbertoCuadraAlbertoCuadra
authored
Merge pull request #1115 from CombustionToolbox/nonideal
Add: include `EquationStatePengRobinson` class and validation with Cantera under thermally perfect gas assumption
2 parents 6182d50 + 1cd3108 commit 789e1ae

5 files changed

Lines changed: 486 additions & 8 deletions

File tree

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

Lines changed: 343 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,343 @@
1+
classdef EquationStatePengRobinson < combustiontoolbox.core.EquationState
2+
% The :mat:class:`EquationStatePengRobinson` class implements the
3+
% Peng-Robinson equation of state for real gases.
4+
%
5+
% Example:
6+
% eos = EquationStatePengRobinson();
7+
%
8+
% See also: :mat:class:`EquationState`, :mat:class:`Mixture`
9+
10+
properties (Access = public)
11+
tol0 = 1e-8; % Tolerance for root finding
12+
end
13+
14+
properties (Access = private)
15+
cachedListSpecies % Cell array of strings to validate the cache
16+
temperatureCritical % Critical temperatures of all species [K]
17+
pressureCritical % Critical pressures of all species [Pa]
18+
acentricFactor % Acentric factors of all species [-]
19+
FLAG_VALID % Flag array identifying species with valid PR data
20+
end
21+
22+
properties (Constant, Access = private)
23+
R0 = combustiontoolbox.common.Constants.R0; % Universal gas constant [J/(K mol)]
24+
end
25+
26+
methods (Access = public)
27+
28+
function pressure = getPressure(obj, temperature, molarVolume, molarFractions, chemicalSystem, varargin)
29+
% Compute pressure [Pa] using the Peng-Robinson equation of state, namely:
30+
%
31+
% .. math::
32+
%
33+
% `P = \\frac{RT}{V - b} - \\frac{a}{V^2 + 2bV - b^2},`
34+
%
35+
% where :math:`a` and :math:`b` are mixture parameters computed using
36+
% van der Waals one-fluid mixing rules.
37+
%
38+
%
39+
% Args:
40+
% obj (EquationStatePengRobinson): Equation of state object
41+
% temperature (float): Temperature of the mixture [K]
42+
% molarVolume (float): Molar volume of the mixture [m3/mol]
43+
% molarFractions (float): Molar fractions of the species in the mixture
44+
% chemicalSystem (ChemicalSystem): Chemical system object containing species data
45+
%
46+
% Returns:
47+
% pressure (float): Pressure of the mixture [Pa]
48+
%
49+
% Example:
50+
% P = getPressure(obj, 300, 0.024, [0.5, 0.5], chemicalSystem)
51+
52+
% Compute mixture parameters
53+
[a_mix, b_mix, ~, ~] = obj.getMixtureParameters(temperature, molarFractions, chemicalSystem);
54+
55+
% Compute pressure [Pa]
56+
pressure = (obj.R0 * temperature) / (molarVolume - b_mix) - ...
57+
a_mix / (molarVolume^2 + 2 * b_mix * molarVolume - b_mix^2);
58+
end
59+
60+
function molarVolume = getVolume(obj, temperature, pressure, molarFractions, chemicalSystem, varargin)
61+
% Compute gas-phase molar volume [m3/mol] by solving the Peng-Robinson
62+
% cubic equation of state for the given temperature and pressure
63+
%
64+
% Args:
65+
% obj (EquationStatePengRobinson): Equation of state object
66+
% temperature (float): Temperature of the mixture [K]
67+
% pressure (float): Pressure of the mixture [Pa]
68+
% molarFractions (float): Molar fractions of the species in the mixture
69+
% chemicalSystem (ChemicalSystem): Chemical system object containing species data
70+
%
71+
% Returns:
72+
% molarVolume (float): Molar volume of the mixture [m3/mol]
73+
%
74+
% Example:
75+
% V = getVolume(obj, 300, 1e5, [0.5, 0.5], chemicalSystem)
76+
77+
% Compute mixture parameters
78+
[a_mix, b_mix, ~, ~] = obj.getMixtureParameters(temperature, molarFractions, chemicalSystem);
79+
80+
% Dimensionless coefficients A and B
81+
A = (a_mix * pressure) / (obj.R0^2 * temperature^2);
82+
B = (b_mix * pressure) / (obj.R0 * temperature);
83+
84+
% Cubic coefficients for Z^3 + c2*Z^2 + c1*Z + c0 = 0
85+
coeffs = [1.0, -(1.0 - B), (A - 2*B - 3*B^2), -(A*B - B^2 - B^3)];
86+
87+
% Solve for Z and pick the largest real root (gas phase)
88+
Z_roots = roots(coeffs);
89+
Z_real = real(Z_roots(abs(imag(Z_roots)) < obj.tol0));
90+
91+
if isempty(Z_real)
92+
error('EquationStatePengRobinson:getVolume', 'No real roots found for Z.');
93+
end
94+
95+
Z_gas = max(Z_real);
96+
97+
% Compute molar volume [m3/mol]
98+
molarVolume = (Z_gas * obj.R0 * temperature) / pressure;
99+
end
100+
101+
function [dPdV_T, dPdT_V] = getPressureDerivativesDimensional(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin)
102+
% Compute dimensional partial pressure derivatives for the mixture assuming frozen chemistry
103+
%
104+
% Args:
105+
% obj (EquationStatePengRobinson): Equation of state object
106+
% temperature (float): Temperature of the mixture [K]
107+
% pressure (float): Pressure of the mixture [Pa]
108+
% molarVolume (float): Molar volume of the mixture [m3/mol]
109+
% molarFractions (float): Molar fractions of the species in the mixture
110+
% chemicalSystem (ChemicalSystem): Chemical system object containing species data
111+
%
112+
% Returns:
113+
% Tuple containing
114+
%
115+
% * dPdV_T (float): Partial derivative of pressure with respect to volume at constant temperature [Pa/(m3/mol)]
116+
% * dPdT_V (float): Partial derivative of pressure with respect to temperature at constant volume [Pa/K]
117+
%
118+
% Example:
119+
% [dPdV_T, dPdT_V] = getPressureDerivativesDimensional(obj, 300, 1e5, 0.024, [0.5, 0.5], chemicalSystem)
120+
121+
% Compute mixture parameters
122+
[a_mix, b_mix, dadT_mix, ~] = obj.getMixtureParameters(temperature, molarFractions, chemicalSystem);
123+
124+
% Compute dimensional pressure derivatives
125+
dPdV_T = -(obj.R0 * temperature) / (molarVolume - b_mix)^2 + (2 * a_mix * (molarVolume + b_mix)) / (molarVolume^2 + 2*b_mix*molarVolume - b_mix^2)^2;
126+
dPdT_V = obj.R0 / (molarVolume - b_mix) - dadT_mix / (molarVolume^2 + 2*b_mix*molarVolume - b_mix^2);
127+
end
128+
129+
function [heatCapacityPressureDeparture, enthalpyDeparture, entropyDeparture] = getDepartureFunctions(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin)
130+
% Compute thermodynamic departure functions for the mixture using the Peng-Robinson equation of state
131+
%
132+
% Args:
133+
% obj (EquationStatePengRobinson): Equation of state object
134+
% temperature (float): Temperature of the mixture [K]
135+
% pressure (float): Pressure of the mixture [Pa]
136+
% molarVolume (float): Molar volume of the mixture [m3/mol]
137+
% molarFractions (float): Molar fractions of the species in the mixture
138+
% chemicalSystem (ChemicalSystem): Chemical system object containing species data
139+
%
140+
% Returns:
141+
% Tuple containing
142+
%
143+
% * heatCapacityPressureDeparture (float): Heat capacity at constant pressure departure [J/(mol-K)]
144+
% * enthalpyDeparture (float): Enthalpy departure [J/mol]
145+
% * entropyDeparture (float): Entropy departure [J/(mol-K)]
146+
%
147+
% Example:
148+
% [dcp, dh, ds] = getDepartureFunctions(obj, 300, 1e5, 0.024, [0.5, 0.5], chemicalSystem)
149+
150+
% Compute mixture parameters
151+
[a_mix, b_mix, dadT_mix, d2adT2_mix] = obj.getMixtureParameters(temperature, molarFractions, chemicalSystem);
152+
153+
% If mixture behaves ideally (b_mix is zero), return zeros for all departure functions
154+
if b_mix < 1e-15
155+
heatCapacityPressureDeparture = 0;
156+
enthalpyDeparture = 0;
157+
entropyDeparture = 0;
158+
return
159+
end
160+
161+
Z = obj.getCompressibilityFactor(temperature, pressure, molarVolume);
162+
B = (b_mix * pressure) / (obj.R0 * temperature);
163+
164+
% Common terms for departure functions
165+
arg = (Z + (1 + sqrt(2)) * B) / (Z + (1 - sqrt(2)) * B);
166+
logTerm = log(max(arg, 1e-12));
167+
denom = 2 * sqrt(2) * b_mix;
168+
169+
% Enthalpy departure [J/mol]
170+
enthalpyDeparture = obj.R0 * temperature * (Z - 1) + ((temperature * dadT_mix - a_mix) / denom) * logTerm;
171+
172+
% Entropy departure [J/(mol-K)]
173+
entropyDeparture = obj.R0 * log(max(Z - B, 1e-12)) + (dadT_mix / denom) * logTerm;
174+
175+
% Heat capacity at constant volume departure [J/(mol-K)]
176+
heatCapacityVolumeDeparture = (temperature * d2adT2_mix / denom) * logTerm;
177+
178+
% Compute pressure derivatives
179+
[dPdV_T, dPdT_V] = obj.getPressureDerivativesDimensional(temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin{:});
180+
181+
% Heat capacity at constant pressure departure [J/(mol-K)]
182+
heatCapacityPressureDeparture = heatCapacityVolumeDeparture + (-temperature * (dPdT_V^2) / dPdV_T) - obj.R0;
183+
end
184+
185+
end
186+
187+
methods (Access = private)
188+
189+
function initializeCache(obj, chemicalSystem)
190+
% Cache the database values once to avoid field lookups during iterative solver loops
191+
%
192+
% Args:
193+
% obj (EquationStatePengRobinson): Equation of state object
194+
% chemicalSystem (ChemicalSystem): Chemical system object containing species data
195+
%
196+
% Example:
197+
% obj.initializeCache(chemicalSystem)
198+
199+
% Definitions
200+
listSpecies = chemicalSystem.listSpecies;
201+
numSpecies = chemicalSystem.numSpecies;
202+
203+
% Preallocate arrays
204+
Tc = zeros(1, numSpecies);
205+
Pc = zeros(1, numSpecies);
206+
omega = zeros(1, numSpecies);
207+
208+
% Extract Species objects from the chemical system
209+
species = chemicalSystem.species;
210+
for i = 1:numSpecies
211+
name = listSpecies{i};
212+
Tc(i) = species.(name).Tcritical;
213+
Pc(i) = species.(name).Pcritical;
214+
omega(i) = species.(name).acentricFactor;
215+
end
216+
217+
obj.temperatureCritical = Tc; % [K]
218+
obj.pressureCritical = Pc * 1e5; % [Pa]
219+
obj.acentricFactor = omega; % [-]
220+
221+
% Identify species with valid PR data (not NaN and > 0)
222+
obj.FLAG_VALID = ~isnan(Tc) & (Tc > 0) & ~isnan(Pc) & (Pc > 0);
223+
224+
% Cache the list of species to validate future calls
225+
obj.cachedListSpecies = listSpecies;
226+
end
227+
228+
function [a_mix, b_mix, dadT_mix, d2adT2_mix] = getMixtureParameters(obj, temperature, molarFractions, chemicalSystem)
229+
% Compute mixture parameters using van der Waals one-fluid mixing rules
230+
%
231+
% Args:
232+
% obj (EquationStatePengRobinson): Equation of state object
233+
% temperature (float): Temperature of the mixture [K]
234+
% molarFractions (float): Molar fractions of the species in the mixture
235+
% chemicalSystem (ChemicalSystem): Chemical system object containing species data
236+
%
237+
% Returns:
238+
% Tuple containing
239+
%
240+
% * a_mix (float): Mixture attraction parameter [J-m3/mol^2]
241+
% * b_mix (float): Mixture co-volume parameter [m3/mol]
242+
% * dadT_mix (float): First temperature derivative of a_mix [J-m3/(mol^2 K)]
243+
% * d2adT2_mix (float): Second temperature derivative of a_mix [J-m3/(mol^2 K^2)]
244+
%
245+
% Example:
246+
% [a_mix, b_mix, dadT_mix, d2adT2_mix] = getMixtureParameters(obj, 300, [0.5, 0.5], chemicalSystem)
247+
248+
% Rebuild cache if empty or if species list changed/reordered
249+
if isempty(obj.cachedListSpecies) || ~isequal(obj.cachedListSpecies, chemicalSystem.listSpecies)
250+
obj.initializeCache(chemicalSystem);
251+
end
252+
253+
% Find species that are both active (>0) AND have valid PR data
254+
FLAG_ACTIVE = (molarFractions(:) > 0) & obj.FLAG_VALID(:);
255+
256+
% If no real species are present, mixture is purely ideal
257+
if ~any(FLAG_ACTIVE)
258+
a_mix = 0; b_mix = 0; dadT_mix = 0; d2adT2_mix = 0;
259+
return;
260+
end
261+
262+
% Extract data for active species
263+
X_active = molarFractions(FLAG_ACTIVE);
264+
Tc = obj.temperatureCritical(FLAG_ACTIVE);
265+
Pc = obj.pressureCritical(FLAG_ACTIVE);
266+
omega = obj.acentricFactor(FLAG_ACTIVE);
267+
268+
% Pure species parameters
269+
kappa = 0.37464 + 1.54226 * omega - 0.26992 * omega.^2;
270+
Tr = temperature ./ Tc;
271+
sqrtTr = sqrt(Tr);
272+
alpha = (1 + kappa .* (1 - sqrtTr)).^2;
273+
274+
a_i = 0.45724 * (obj.R0^2 * Tc.^2 ./ Pc) .* alpha;
275+
b_i = 0.07780 * (obj.R0 * Tc ./ Pc);
276+
277+
% Temperature derivatives
278+
dalpha_dT = -kappa ./ (sqrtTr .* Tc) .* (1 + kappa .* (1 - sqrtTr));
279+
d2alpha_dT2 = kappa .* (kappa + 1) ./ (2 * Tc.^2 .* Tr.^(3/2));
280+
281+
a0 = 0.45724 * (obj.R0^2 * Tc.^2 ./ Pc);
282+
da_dT_i = a0 .* dalpha_dT;
283+
d2a_dT2_i = a0 .* d2alpha_dT2;
284+
285+
% Apply mixing rules
286+
b_mix = dot(X_active, b_i);
287+
288+
sqrt_a_i = sqrt(max(a_i, 1e-20));
289+
S1 = dot(X_active, sqrt_a_i);
290+
a_mix = S1^2;
291+
292+
S2 = dot(X_active, da_dT_i ./ (2 * sqrt_a_i));
293+
dadT_mix = 2 * S1 * S2;
294+
295+
S3 = dot(X_active, (d2a_dT2_i ./ (2 * sqrt_a_i)) - (da_dT_i.^2 ./ (4 * sqrt_a_i.^3)));
296+
d2adT2_mix = 2 * S2^2 + 2 * S1 * S3;
297+
end
298+
299+
function [temperatureCritical_mix, pressureCritical_mix, acentricFactor_mix] = getPseudoCriticalProperties(obj, molarFractions, chemicalSystem)
300+
% Computes pseudo-critical properties for multi-component mixtures
301+
%
302+
% Args:
303+
% obj (EquationStatePengRobinson): Equation of state object
304+
% molarFractions (float): Molar fractions of the species in the mixture
305+
% chemicalSystem (ChemicalSystem): Chemical system object containing species data
306+
%
307+
% Returns:
308+
% Tuple containing
309+
%
310+
% * temperatureCritical_mix (float): Pseudo-critical temperature of the mixture [K]
311+
% * pressureCritical_mix (float): Pseudo-critical pressure of the mixture [Pa]
312+
% * acentricFactor_mix (float): Pseudo-critical acentric factor of the mixture [-]
313+
%
314+
% Example:
315+
% [temperatureCritical_mix, pressureCritical_mix, acentricFactor_mix] = getPseudoCriticalProperties(obj, [0.5, 0.5], chemicalSystem)
316+
317+
if isempty(obj.cachedListSpecies) || ~isequal(obj.cachedListSpecies, chemicalSystem.listSpecies)
318+
obj.initializeCache(chemicalSystem);
319+
end
320+
321+
% Definitions
322+
mask = (molarFractions(:) > 0) & obj.FLAG_VALID(:);
323+
Xi = molarFractions(mask);
324+
Tc_i = obj.temperatureCritical(mask);
325+
Pc_i = obj.pressureCritical(mask);
326+
omega_i = obj.acentricFactor(mask);
327+
328+
% a and b at critical point (alpha = 1)
329+
a_i_tc = 0.45724 * (obj.R0^2 * Tc_i.^2 ./ Pc_i);
330+
b_i = 0.07780 * (obj.R0 * Tc_i ./ Pc_i);
331+
332+
% Mixture parameters at critical condition
333+
a_mix_tc = ( dot(Xi, sqrt(a_i_tc)) )^2;
334+
b_mix = dot(Xi, b_i);
335+
336+
% Back-calculate pseudo-critical T and P
337+
temperatureCritical_mix = (a_mix_tc * 0.07780) / (b_mix * 0.45724 * obj.R0);
338+
pressureCritical_mix = (0.07780 * obj.R0 * temperatureCritical_mix) / b_mix;
339+
acentricFactor_mix = dot(Xi, omega_i);
340+
end
341+
342+
end
343+
end

0 commit comments

Comments
 (0)