From 6fd0b7faa683e73b071e23275128da6f400c0925 Mon Sep 17 00:00:00 2001 From: AlbertoCuadra Date: Thu, 26 Mar 2026 14:36:29 +0100 Subject: [PATCH 1/6] Update: include Peng-Robinson data in the database --- .../Peng-Robinson/addPengRobinsonProperties.m | 42 +++++++++++++++++++ .../eos/Peng-Robinson/assign_DB_PR_eos.m | 36 ---------------- 2 files changed, 42 insertions(+), 36 deletions(-) create mode 100644 databases/eos/Peng-Robinson/addPengRobinsonProperties.m delete mode 100644 databases/eos/Peng-Robinson/assign_DB_PR_eos.m diff --git a/databases/eos/Peng-Robinson/addPengRobinsonProperties.m b/databases/eos/Peng-Robinson/addPengRobinsonProperties.m new file mode 100644 index 000000000..f5e0197cb --- /dev/null +++ b/databases/eos/Peng-Robinson/addPengRobinsonProperties.m @@ -0,0 +1,42 @@ +function speciesArray = addPengRobinsonProperties(speciesArray) + % Assign parameters for pure species for Peng-Robinson EoS. + % + % References: + % [1] THERMODYNAMIC AND THERMOPHYSICAL PROPERTIES OF HUMID AIR BY USING CUBIC PENG-ROBINSON EOS + % [2] PROPERTIES OF GASES, Isidoro Martinez + + % 1. Initialize all existing species in the struct with NaN. + % This prevents [] (empty) values which would crash vectorized logical checks. + fnames = fieldnames(speciesArray); + for i = 1:length(fnames) + speciesArray.(fnames{i}).Tcritical = NaN; + speciesArray.(fnames{i}).Pcritical = NaN; + speciesArray.(fnames{i}).acentricFactor = NaN; + end + + % 2. Assign known values + speciesArray = addSpecies(speciesArray, 'Ar', 150.87, 48.98, 0.000); + speciesArray = addSpecies(speciesArray, 'CH4', 190.56, 45.99, 0.011); + speciesArray = addSpecies(speciesArray, 'C2H4', 282.35, 50.42, 0.087); + speciesArray = addSpecies(speciesArray, 'C2H6', 305.32, 48.72, 0.099); + speciesArray = addSpecies(speciesArray, 'C3H8', 369.83, 42.48, 0.152); + speciesArray = addSpecies(speciesArray, 'CO', 132.86, 34.94, 0.045); + speciesArray = addSpecies(speciesArray, 'CO2', 304.13, 73.75, 0.225); + speciesArray = addSpecies(speciesArray, 'O2', 154.581, 50.429, 0.0222); + speciesArray = addSpecies(speciesArray, 'N2', 126.19, 33.96, 0.037); + speciesArray = addSpecies(speciesArray, 'NO', 180.15, 64.80, 0.588); + speciesArray = addSpecies(speciesArray, 'NO2', 431.20, 101.32, 0.851); + speciesArray = addSpecies(speciesArray, 'H2', 33.19, 13.13, -0.219); + speciesArray = addSpecies(speciesArray, 'H2O', 647.14, 220.64, 0.344); + speciesArray = addSpecies(speciesArray, 'H2ObLb', 647.14, 220.64, 0.344); +end + +% SUB-PASS FUNCTIONS +function speciesArray = addSpecies(speciesArray, species, Tcritical, Pcritical, acentricFactor) + % Only assign if the species name actually exists as a field in the struct + assert(isfield(speciesArray, species), 'Species "%s" not found in the input struct. Check for typos or missing entries.', species); + + speciesArray.(species).Tcritical = Tcritical; % [K] + speciesArray.(species).Pcritical = Pcritical; % [bar] + speciesArray.(species).acentricFactor = acentricFactor; % [-] +end \ No newline at end of file diff --git a/databases/eos/Peng-Robinson/assign_DB_PR_eos.m b/databases/eos/Peng-Robinson/assign_DB_PR_eos.m deleted file mode 100644 index 2028d28d2..000000000 --- a/databases/eos/Peng-Robinson/assign_DB_PR_eos.m +++ /dev/null @@ -1,36 +0,0 @@ -function DB = assign_DB_PR_eos(DB) - % Assing parameters for pure species considering Peng-Robinson Equation - % of State (EoS) - % - % - % References: - % [1] THERMODYNAMIC AND THERMOPHYSICAL PROPERTIES OF HUMID AIR BY USING CUBIC PENG-ROBINSON EOS - % [2] PROPERTIES OF GASES, Isidoro Martinez - - DB = add_species(DB, 'Ar', 151, 48.6 * 1e5, 0); % Ref. [2] - - DB = add_species(DB, 'CH4', 190.6, 46 * 1e5, 0.008); - DB = add_species(DB, 'C2H4', 283, 51.2 * 1e5, 0.085); % Ref. [2] - DB = add_species(DB, 'C2H6', 305.4, 48.84 * 1e5, 0.098); - DB = add_species(DB, 'C3H8', 369.8, 42.46 * 1e5, 0.152); - - DB = add_species(DB, 'CO', 132.86, 34.94 * 1e5, 0.0497); % Ref. [1] - DB = add_species(DB, 'CO2', 304.2, 73.76 * 1e5, 0.225); - - DB = add_species(DB, 'O2', 154.6, 50.4 * 1e5, 0.025); % Ref. [1] - DB = add_species(DB, 'N2', 126.2, 33.9 * 1e5, 0.039); % Ref. [1] - DB = add_species(DB, 'NO', 180, 64.8 * 1e5, 0.588); % Ref. [1] - DB = add_species(DB, 'NO2', 261.9, 101 * 1e5, 0.834); % Ref. [1] - - DB = add_species(DB, 'H2', 33, 13.2 * 1e5, -0.22); % Ref. [2] - DB = add_species(DB, 'H2O', 647.13, 22.34 * 1e5, 0.344); - DB = add_species(DB, 'H2ObLb', 647.13, 22.34 * 1e5, 0.344); -end - -% SUB-PASS FUNCTIONS -function DB = add_species(DB, species, T_c, p_c, w) - % Coefficients - DB.(species).critical_temperature = T_c; % [K] - DB.(species).critical_pressure = p_c; % [Pa] - DB.(species).acentric_factor = w; % [-] -end From 19fe258a9a5b994d3ee816ae91083a9c0013bea4 Mon Sep 17 00:00:00 2001 From: AlbertoCuadra Date: Thu, 26 Mar 2026 14:37:22 +0100 Subject: [PATCH 2/6] Add: include parameters required in cubic EoS --- +combustiontoolbox/+core/@Species/Species.m | 46 +++++++++++---------- 1 file changed, 25 insertions(+), 21 deletions(-) diff --git a/+combustiontoolbox/+core/@Species/Species.m b/+combustiontoolbox/+core/@Species/Species.m index 5695dffaf..2d210e1da 100644 --- a/+combustiontoolbox/+core/@Species/Species.m +++ b/+combustiontoolbox/+core/@Species/Species.m @@ -4,27 +4,31 @@ % See also: :mat:func:`Database`, :mat:func:`NasaDatabase` properties - name % Name chemical species - fullname % Fullname chemical species - refCode % Reference date code - comments % Additional comments from database - formula % Chemical formula - W % Molecular weight [kg/mol] - hf % Enthalpy of formation at Tref from its reference species in their standard state [J/mol] - hftoh0 % Enthalpy of formation at Tref relative to molar enthalpy at 0 K for standard state [J/mol] - ef % Internal energy of formation [J/mol] - phase % Phase - T % Temperature [K] - Tref % Temperature reference [K] - Trange % Temperature range intervals - Tintervals % Number of temperature intervals - Texponents % Exponents polynomials - a % Coefficients a polynomials - b % Coefficients b polynomials - cpcurve % Gridded interpolant object with specific heat at constant pressure of the individual species - h0curve % Gridded interpolant object with enthalpy of the individual species - s0curve % Gridded interpolant object with entropy of the individual species - g0curve % Gridded interpolant object with Gibbs free energy of the individual species + name % Name chemical species + fullname % Fullname chemical species + refCode % Reference date code + comments % Additional comments from database + formula % Chemical formula + W % Molecular weight [kg/mol] + hf % Enthalpy of formation at Tref from its reference species in their standard state [J/mol] + hftoh0 % Enthalpy of formation at Tref relative to molar enthalpy at 0 K for standard state [J/mol] + ef % Internal energy of formation [J/mol] + phase % Phase + T % Temperature [K] + Tref % Temperature reference [K] + Trange % Temperature range intervals + Tintervals % Number of temperature intervals + Texponents % Exponents polynomials + a % Coefficients a polynomials + b % Coefficients b polynomials + Tcritical % Critical temperature [K] + Pcritical % Critical pressure [bar] + Vcritical % Critical molar volume [m3/mol] + acentricFactor % Acentric factor [-] + cpcurve % Gridded interpolant object with specific heat at constant pressure of the individual species + h0curve % Gridded interpolant object with enthalpy of the individual species + s0curve % Gridded interpolant object with entropy of the individual species + g0curve % Gridded interpolant object with Gibbs free energy of the individual species end properties (Access = ?combustiontoolbox.databases.Database) From 95272c4602fe800445d80744d24a4c30f78fcb4a Mon Sep 17 00:00:00 2001 From: AlbertoCuadra Date: Thu, 26 Mar 2026 14:38:47 +0100 Subject: [PATCH 3/6] Update: rewrite functions to be compatible with cubic EoS --- .../+core/@EquationState/EquationState.m | 43 +++++++++++-- .../EquationStateIdealGas.m | 62 ++++++++++++++----- +combustiontoolbox/+core/@Mixture/Mixture.m | 47 +++++++++----- .../@ShockSolver/shockIncident.m | 2 +- 4 files changed, 117 insertions(+), 37 deletions(-) diff --git a/+combustiontoolbox/+core/@EquationState/EquationState.m b/+combustiontoolbox/+core/@EquationState/EquationState.m index 4b6b34043..60afbbdbb 100644 --- a/+combustiontoolbox/+core/@EquationState/EquationState.m +++ b/+combustiontoolbox/+core/@EquationState/EquationState.m @@ -4,17 +4,50 @@ % mixture using a specified equation of state. % % Subclasses must implement the following abstract methods: - % * getPressure(temperature, molarVolume, varargin) - % * getVolume(temperature, pressure, varargin) + % * getPressure(obj,temperature, molarVolume, varargin) + % * getVolume(obj, temperature, pressure, varargin) + % * getDepartureFunctions(obj, temperature, pressure, varargin) + % * getVolumeDerivatives(obj, temperature, pressure, molarVolume, varargin) % % See also: :mat:func:`EquationStateIdealGas`, :mat:func:`Mixture` - + methods (Abstract) % Compute pressure [Pa] given the temperature and molar volume. - pressure = getPressure(obj, temperature, molarVolume, varargin) + pressure = getPressure(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) % Compute molar volume [m3/mol] given the temperature and pressure. - molarVolume = getVolume(obj, temperature, pressure, varargin) + molarVolume = getVolume(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + + % Compute thermodynamic departure functions + [heatCapacityPressureDeparture, enthalpyDeparture, entropyDeparture] = getDepartureFunctions(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + + % Compute dimensionless volume derivatives: (dlnV/dlnT)_p and (dlnV/dlnP)_T + [dlnVdT_p, dlnVdp_T] = getVolumeDerivatives(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + end + + methods (Access = public, Static) + + function Z = getCompressibilityFactor(temperature, pressure, molarVolume) + % Compute compressibility factor + % + % Args: + % temperature (float): Temperature of the mixture [K] + % pressure (float): Pressure of the mixture [Pa] + % molarVolume (float): Molar volume of the mixture [m3/mol] + % + % Returns: + % Z (float): Compressibility factor [-] + % + % Example: + % Z = EquationState.getCompressibilityFactor(300, 1e5, 0.024) + + % Definitions + R0 = combustiontoolbox.common.Constants.R0; % Universal gas constant [J/(K mol)] + + % Compute compressibility factor + Z = (pressure * molarVolume) / (R0 * temperature); + end + end end \ No newline at end of file diff --git a/+combustiontoolbox/+core/@EquationStateIdealGas/EquationStateIdealGas.m b/+combustiontoolbox/+core/@EquationStateIdealGas/EquationStateIdealGas.m index 096c577e5..4d304159e 100644 --- a/+combustiontoolbox/+core/@EquationStateIdealGas/EquationStateIdealGas.m +++ b/+combustiontoolbox/+core/@EquationStateIdealGas/EquationStateIdealGas.m @@ -1,17 +1,22 @@ classdef EquationStateIdealGas < combustiontoolbox.core.EquationState - % The :mat:func:`EquationStateIdealGas` class implements the ideal gas law. + % The :mat:func:`EquationStateIdealGas` class implements the ideal gas equation of state. % - % The :mat:func:`EquationStateIdealGas` object can be used to compute - % the pressure and molar volume of a mixture assuming ideal gas behavior. + % The :mat:func:`EquationStateIdealGas` object can be used to compute the pressure, + % molar volume, dimensionless volume derivatives, and thermodynamic departure functions + % of a mixture assuming ideal gas behavior. % % Example: % eos = EquationStateIdealGas(); % % See also: :mat:func:`EquationState` and :mat:func:`Mixture` - methods (Access = public, Static) + properties (Constant, Access = private) + R0 = combustiontoolbox.common.Constants.R0; % Universal gas constant [J/(K mol)] + end + + methods (Access = public) - function pressure = getPressure(temperature, molarVolume, varargin) + function pressure = getPressure(obj, temperature, molarVolume, ~, ~, varargin) % Compute pressure [Pa] using the ideal gas law. % % Args: @@ -24,14 +29,11 @@ % Example: % P = getPressure(obj, 300, 0.024) - % Definitions - R0 = combustiontoolbox.common.Constants.R0; % Universal gas constant [J/(K mol)] - % Compute pressure [Pa] - pressure = (R0 .* temperature) ./ molarVolume; + pressure = (obj.R0 .* temperature) ./ molarVolume; end - function molarVolume = getVolume(temperature, pressure, varargin) + function molarVolume = getVolume(obj, temperature, pressure, ~, ~, varargin) % Compute molar volume [m3/mol] using the ideal gas law. % % Args: @@ -44,11 +46,43 @@ % Example: % V = getVolume(obj, 300, 1e5) - % Definitions - R0 = combustiontoolbox.common.Constants.R0; % Universal gas constant [J/(K mol)] - % Compute molar volume [m3/mol] - molarVolume = R0 * temperature ./ pressure; + molarVolume = obj.R0 * temperature ./ pressure; + end + + function [dVdT_p, dVdp_T] = getVolumeDerivatives(~, ~, ~, ~, ~, ~, varargin) + % Compute dimensionless volume derivatives for the mixture assuming frozen chemistry. + % + % Returns: + % Tuple containing + % + % * dVdT_p (float): Logarithmic derivative of volume with respect to temperature at constant pressure [-] + % * dVdp_T (float): Logarithmic derivative of volume with respect to pressure at constant temperature [-] + % + % Example: + % [dVdT_p, dVdp_T] = getVolumeDerivatives(obj) + + % For an ideal gas (V = RT/p), the dimensionless derivatives are exactly 1 and -1. + dVdT_p = 1; + dVdp_T = -1; + end + + function [heatCapacityPressureDeparture, enthalpyDeparture, entropyDeparture] = getDepartureFunctions(~, ~, ~, ~, ~, ~, varargin) + % Compute thermodynamic departure functions for an ideal gas. For an ideal gas, all departure functions are zero. + % + % Returns: + % Tuple containing + % + % * heatCapacityPressureDeparture (float): Heat capacity at constant pressure departure [J/(mol-K)] + % * enthalpyDeparture (float): Enthalpy departure [J/mol] + % * entropyDeparture (float): Entropy departure [J/(mol-K)] + % + % Example: + % [dcp, dh, ds] = getDepartureFunctions(obj) + + heatCapacityPressureDeparture = 0; + enthalpyDeparture = 0; + entropyDeparture = 0; end end diff --git a/+combustiontoolbox/+core/@Mixture/Mixture.m b/+combustiontoolbox/+core/@Mixture/Mixture.m index ad46f2b54..8946750fc 100644 --- a/+combustiontoolbox/+core/@Mixture/Mixture.m +++ b/+combustiontoolbox/+core/@Mixture/Mixture.m @@ -1007,7 +1007,7 @@ % Compute molar volume [m3/mol] from specific volume [m3/kg] vMolar = vSpecific2vMolar(obj, obj.vSpecific, obj.quantity, obj.quantity(obj.indexGas)); % Compute pressure in Pascals [Pa] using the equationState - pressure = obj.equationState.getPressure(obj.T, vMolar, obj.chemicalSystem.listSpecies, obj.quantity / sum(obj.quantity)); + pressure = obj.equationState.getPressure(obj.T, vMolar, obj.quantity / sum(obj.quantity), obj.chemicalSystem); % Convert pressure to [bar] obj.p = pressure * combustiontoolbox.common.Units.Pa2bar; end @@ -1467,7 +1467,7 @@ function computeComposition(obj) % Definitions temperature = obj.T; - pressure = obj.p; + pressure = obj.p; pressure_Pa = pressure * combustiontoolbox.common.Units.bar2Pa; R0 = combustiontoolbox.common.Constants.R0; % Universal gas constant [J/(K mol)] system = obj.chemicalSystem; propertiesMatrix = system.propertiesMatrix; % Properties matrix @@ -1489,10 +1489,12 @@ function computeComposition(obj) % Compute volume [m3] if obj.FLAG_VOLUME obj.v = obj.vSpecific * obj.mi; + molarVolume = obj.v / N_gas; % Update pressure [bar] if specific volume is given - obj.p = obj.equationState.getPressure(temperature, obj.v / N_gas, obj.chemicalSystem.listSpecies, obj.Xi) * combustiontoolbox.common.Units.Pa2bar; + obj.p = obj.equationState.getPressure(temperature, obj.v / N_gas, obj.Xi, obj.chemicalSystem) * combustiontoolbox.common.Units.Pa2bar; else - obj.v = obj.equationState.getVolume(temperature, pressure * combustiontoolbox.common.Units.bar2Pa, obj.chemicalSystem.listSpecies, obj.Xi) * N_gas; + molarVolume = obj.equationState.getVolume(temperature, pressure_Pa, obj.Xi, obj.chemicalSystem); + obj.v = molarVolume * N_gas; end % Compute specific volume [m3/kg] @@ -1501,8 +1503,16 @@ function computeComposition(obj) % Compute density [kg/m3] obj.rho = 1 / obj.vSpecific; + % Compute thermodynamic departure (real gas effects) [psiDeparture = psiNonIdealGas - psiDepartureIdealGas] + [cp_dep_molar, h_dep_molar, s_dep_molar] = obj.equationState.getDepartureFunctions(... + temperature, pressure_Pa, molarVolume, obj.Xi, obj.chemicalSystem); + + obj.cp = obj.cp + cp_dep_molar * N_gas; + obj.h = obj.h + h_dep_molar * N_gas; + obj.s0 = obj.s0 + s_dep_molar * N_gas; + % Compute internal energy [J] - obj.e = obj.h - N_gas * R0 * temperature; + obj.e = obj.h - pressure_Pa * obj.v; % Compute thermal internal energy [J] obj.DeT = obj.e - obj.ef; @@ -1518,17 +1528,21 @@ function computeComposition(obj) % Compute Gibbs energy [J] obj.g = obj.h - obj.T * obj.s; - + + % Retrieve frozen dimensionless derivatives: (dlnV/dlnT)_p, (dlnV/dlnP)_T + [dVdT_p_frozen, dVdp_T_frozen] = obj.equationState.getVolumeDerivatives(... + temperature, pressure_Pa, molarVolume, obj.Xi, obj.chemicalSystem); + % Compute frozen component of the specific heats [J/K] obj.cp_f = obj.cp; - obj.cv_f = obj.cp_f - R0 * N_gas; + obj.cv_f = obj.cp_f + (pressure_Pa * obj.v / temperature) * (dVdT_p_frozen^2) / dVdp_T_frozen; % Compute thermodynamic derivatives, cp, cv, gamma, and speed % of sound considering chemical reaction if obj.FLAG_REACTION % Set thermodynamic derivatives - obj.dVdT_p = obj.dN_T + 1; % [-] - obj.dVdp_T = obj.dN_p - 1; % [-] + obj.dVdT_p = obj.dN_T + dVdT_p_frozen; % [-] + obj.dVdp_T = obj.dN_p + dVdp_T_frozen; % [-] if ~any(isnan(obj.dNi_T)) && ~any(isinf(obj.dNi_T)) % Definitions @@ -1539,14 +1553,14 @@ function computeComposition(obj) obj.cp = obj.cp_f + sum(h0_j / temperature .* (1 + delta .* (Ni - 1)) .* obj.dNi_T, 'omitnan'); % Compute specific heat at constant volume [J/K] - obj.cv = obj.cp + (N_gas * R0 * obj.dVdT_p^2) / obj.dVdp_T; + obj.cv = obj.cp + (pressure_Pa * obj.v / temperature) * (obj.dVdT_p^2) / obj.dVdp_T; % Compute Adibatic index [-] obj.gamma = obj.cp / obj.cv; obj.gamma_s = -obj.gamma / obj.dVdp_T; % Compute sound velocity [m/s] - obj.sound = sqrt(obj.gamma_s * R0 * obj.T / obj.W); + obj.sound = sqrt(obj.gamma_s * pressure_Pa * obj.vSpecific); % Compute Mach number if ~isempty(obj.u) @@ -1565,19 +1579,18 @@ function computeComposition(obj) % of sound considering frozen chemistry % Set thermodynamic derivatives [-] - obj.dVdT_p = 1; - obj.dVdp_T = -1; + obj.dVdT_p = dVdT_p_frozen; + obj.dVdp_T = dVdp_T_frozen; % Compute specific heat at constant volume [J/K] - obj.cv = obj.cp - R0 * N_gas; + obj.cv = obj.cv_f; % Compute adibatic index [-] obj.gamma = obj.cp / obj.cv; - obj.gamma_s = obj.gamma; + obj.gamma_s = -obj.gamma / obj.dVdp_T; % Compute sound velocity [m/s] - obj.sound = sqrt(obj.gamma * R0 * obj.T / obj.W); - + obj.sound = sqrt(obj.gamma_s * pressure_Pa * obj.vSpecific); % Compute Mach number if ~isempty(obj.u) obj.mach = obj.u / obj.sound; diff --git a/+combustiontoolbox/+shockdetonation/@ShockSolver/shockIncident.m b/+combustiontoolbox/+shockdetonation/@ShockSolver/shockIncident.m index 512372440..87f5ee066 100644 --- a/+combustiontoolbox/+shockdetonation/@ShockSolver/shockIncident.m +++ b/+combustiontoolbox/+shockdetonation/@ShockSolver/shockIncident.m @@ -176,7 +176,7 @@ h2 = mix2.h / mix2.mi; % [J/kg] cp2 = mix2.cp / mix2.mi; % [J/(K-kg)] - alpha = (W1 * u1^2) / (R0 * T1); + alpha = (r1 * u1^2) / (p1 * 1e5); J1 = -r1r2 * alpha * dVdp_T - p2p1; J2 = -r1r2 * alpha * dVdT_p; b1 = p2p1 - 1 + alpha * (r1 / r2 - 1); From dd95ada79295041e067f40908232f3539e263e49 Mon Sep 17 00:00:00 2001 From: AlbertoCuadra Date: Thu, 26 Mar 2026 19:10:16 +0100 Subject: [PATCH 4/6] Update: generalize `EquationState` and `EquationStateIdealGas` --- .../+core/@EquationState/EquationState.m | 108 ++++++++++++++++-- .../EquationStateIdealGas.m | 23 ++++ 2 files changed, 123 insertions(+), 8 deletions(-) diff --git a/+combustiontoolbox/+core/@EquationState/EquationState.m b/+combustiontoolbox/+core/@EquationState/EquationState.m index 60afbbdbb..f4db9461c 100644 --- a/+combustiontoolbox/+core/@EquationState/EquationState.m +++ b/+combustiontoolbox/+core/@EquationState/EquationState.m @@ -4,25 +4,117 @@ % mixture using a specified equation of state. % % Subclasses must implement the following abstract methods: - % * getPressure(obj,temperature, molarVolume, varargin) - % * getVolume(obj, temperature, pressure, varargin) - % * getDepartureFunctions(obj, temperature, pressure, varargin) - % * getVolumeDerivatives(obj, temperature, pressure, molarVolume, varargin) + % * getPressure(obj,temperature, molarVolume, molarFractions, chemicalSystem, varargin) + % * getVolume(obj, temperature, pressure, molarFractions, chemicalSystem, varargin) + % * getDepartureFunctions(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + % * getPressureDerivativesDimensional(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) % % See also: :mat:func:`EquationStateIdealGas`, :mat:func:`Mixture` methods (Abstract) % Compute pressure [Pa] given the temperature and molar volume. - pressure = getPressure(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + pressure = getPressure(obj, temperature, molarVolume, molarFractions, chemicalSystem, varargin) % Compute molar volume [m3/mol] given the temperature and pressure. - molarVolume = getVolume(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + molarVolume = getVolume(obj, temperature, pressure, molarFractions, chemicalSystem, varargin) % Compute thermodynamic departure functions [heatCapacityPressureDeparture, enthalpyDeparture, entropyDeparture] = getDepartureFunctions(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) - % Compute dimensionless volume derivatives: (dlnV/dlnT)_p and (dlnV/dlnP)_T - [dlnVdT_p, dlnVdp_T] = getVolumeDerivatives(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + % Compute dimensional pressure derivatives: (dP/dV)_T and (dP/dT)_V + [dPdV_T, dPdT_V] = getPressureDerivativesDimensional(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin); + end + + methods (Access = public) + + function [dPdV_T, dPdT_V] = getPressureDerivatives(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + % Compute dimensionless (logarithmic) partial pressure derivatives for the mixture + % + % Args: + % obj (EquationState): Equation of state object + % temperature (float): Temperature of the mixture [K] + % pressure (float): Pressure of the mixture [Pa] + % molarVolume (float): Molar volume of the mixture [m3/mol] + % molarFractions (float): Molar fractions of the species in the mixture + % chemicalSystem (ChemicalSystem): Chemical system object containing species data + % + % Returns: + % Tuple containing + % + % * dPdV_T (float): Logarithmic derivative of pressure with respect to volume at constant temperature [-] + % * dPdT_V (float): Logarithmic derivative of pressure with respect to temperature at constant volume [-] + + % Get dimensional pressure derivatives from the specific EoS implementation + [dPdV_T, dPdT_V] = obj.getPressureDerivativesDimensional(temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin{:}); + + % Convert to dimensionless (logarithmic) pressure derivatives + % (dlnP/dlnV)_T = (V/P) * (dP/dV)_T + % (dlnP/dlnT)_V = (T/P) * (dP/dT)_V + dPdV_T = (molarVolume / pressure) * dPdV_T; + dPdT_V = (temperature / pressure) * dPdT_V; + end + + function [dVdT_p, dVdp_T] = getVolumeDerivatives(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + % Compute dimensionless (logarithmic) volume derivatives for the mixture assuming frozen chemistry + % + % Args: + % obj (EquationStatePengRobinson): Equation of state object + % temperature (float): Temperature of the mixture [K] + % pressure (float): Pressure of the mixture [Pa] + % molarVolume (float): Molar volume of the mixture [m3/mol] + % molarFractions (float): Molar fractions of the species in the mixture + % chemicalSystem (ChemicalSystem): Chemical system object containing species data + % + % Returns: + % Tuple containing + % + % * dVdT_p (float): Logarithmic derivative of volume with respect to temperature at constant pressure [-] + % * dVdp_T (float): Logarithmic derivative of volume with respect to pressure at constant temperature [-] + % + % Example: + % [dVdT_p, dVdp_T] = getVolumeDerivatives(obj, 300, 1e5, 0.024, [0.5, 0.5], chemicalSystem) + + % Compute dimensional pressure derivatives + [dPdV_T, dPdT_V] = obj.getPressureDerivativesDimensional(temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin{:}); + + % Convert to volume derivatives using the chain rule: + % + % (dlnV/dlnT)_p = -(T/V) * (dP/dT)_V / (dP/dV)_T + % (dlnV/dlnP)_T = (P/V) * (1 / (dP/dV)_T) + dVdT_p = (temperature / molarVolume) * (-dPdT_V / dPdV_T); + dVdp_T = (pressure / molarVolume) * (1 / dPdV_T); + end + + function [dVdT_p, dVdp_T] = getVolumeDerivativesDimensional(obj, temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin) + % Compute dimensional volume derivatives for the mixture assuming frozen chemistry + % + % Args: + % obj (EquationState): Equation of state object + % temperature (float): Temperature of the mixture [K] + % pressure (float): Pressure of the mixture [Pa] + % molarVolume (float): Molar volume of the mixture [m3/mol] + % molarFractions (float): Molar fractions of the species in the mixture + % chemicalSystem (ChemicalSystem): Chemical system object containing species data + % + % Returns: + % Tuple containing + % + % * dVdT_p (float): Partial derivative of volume with respect to temperature at constant pressure [m3/(mol-K)] + % * dVdp_T (float): Partial derivative of volume with respect to pressure at constant temperature [m3/(mol-Pa)] + % + % Example: + % [dVdT_p, dVdp_T] = getVolumeDerivativesDimensional(obj, 300, 1e5, 0.024, [0.5, 0.5], chemicalSystem) + + % Get dimensional pressure derivatives from the specific EoS implementation + [dPdV_T, dPdT_V] = obj.getPressureDerivativesDimensional(temperature, pressure, molarVolume, molarFractions, chemicalSystem, varargin{:}); + + % Compute dimensional volume derivatives + % (dV/dT)_p = -(dP/dT)_V / (dP/dV)_T + % (dV/dp)_T = 1 / (dP/dV)_T + dVdT_p = -dPdT_V / dPdV_T; + dVdp_T = 1 / dPdV_T; + end + end methods (Access = public, Static) diff --git a/+combustiontoolbox/+core/@EquationStateIdealGas/EquationStateIdealGas.m b/+combustiontoolbox/+core/@EquationStateIdealGas/EquationStateIdealGas.m index 4d304159e..a6bc8f15a 100644 --- a/+combustiontoolbox/+core/@EquationStateIdealGas/EquationStateIdealGas.m +++ b/+combustiontoolbox/+core/@EquationStateIdealGas/EquationStateIdealGas.m @@ -50,6 +50,29 @@ molarVolume = obj.R0 * temperature ./ pressure; end + function [dPdV_T, dPdT_V] = getPressureDerivativesDimensional(obj, temperature, ~, molarVolume, ~, ~, varargin) + % Compute dimensional partial pressure derivatives for an ideal gas + % + % Args: + % obj (EquationStateIdealGas): Equation of state object + % temperature (float): Temperature of the mixture [K] + % pressure (float): Pressure of the mixture [Pa] + % molarVolume (float): Molar volume of the mixture [m3/mol] + % + % Returns: + % Tuple containing + % + % * dPdV_T (float): Partial derivative of pressure with respect to volume at constant temperature [Pa/(m3/mol)] + % * dPdT_V (float): Partial derivative of pressure with respect to temperature at constant volume [Pa/K] + % + % Example: + % [dPdV_T, dPdT_V] = getPressureDerivativesDimensional(obj, 300, 1e5, 0.024) + + % Compute dimensional pressure derivatives for Ideal Gas: P = R0*T/V + dPdV_T = -(obj.R0 * temperature) / (molarVolume^2); + dPdT_V = obj.R0 / molarVolume; + end + function [dVdT_p, dVdp_T] = getVolumeDerivatives(~, ~, ~, ~, ~, ~, varargin) % Compute dimensionless volume derivatives for the mixture assuming frozen chemistry. % From fb33a367aedd054eaaefe291cc6c19816227cef8 Mon Sep 17 00:00:00 2001 From: AlbertoCuadra Date: Thu, 26 Mar 2026 19:21:40 +0100 Subject: [PATCH 5/6] Solve: typo --- +combustiontoolbox/+utils/+display/interpreterLabel.m | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/+combustiontoolbox/+utils/+display/interpreterLabel.m b/+combustiontoolbox/+utils/+display/interpreterLabel.m index 0e9dab20f..b87de04d9 100644 --- a/+combustiontoolbox/+utils/+display/interpreterLabel.m +++ b/+combustiontoolbox/+utils/+display/interpreterLabel.m @@ -374,11 +374,11 @@ property_unit = ''; case {'dvdtp', 'dvdt_p'} property_name = ''; - property_latex = '(\rm{d}v/\rm{d}T)_p'; + property_latex = '(\rm{d\,ln}v/\rm{d\,ln}T)_p'; property_unit = ''; case {'dvdpt', 'dvdp_t'} property_name = ''; - property_latex = '(\rm{d}v/\rm{d}p)_T'; + property_latex = '(\rm{d\,ln}v/\rm{d\,ln}p)_T'; property_unit = ''; case 'theta' property_name = 'Deflection angle'; From 6c0643548298cf8397a1bb6c5821f56d4206d305 Mon Sep 17 00:00:00 2001 From: AlbertoCuadra Date: Thu, 26 Mar 2026 19:54:26 +0100 Subject: [PATCH 6/6] Solve: typo --- +combustiontoolbox/+core/@Mixture/Mixture.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/+combustiontoolbox/+core/@Mixture/Mixture.m b/+combustiontoolbox/+core/@Mixture/Mixture.m index 8946750fc..50f462ddf 100644 --- a/+combustiontoolbox/+core/@Mixture/Mixture.m +++ b/+combustiontoolbox/+core/@Mixture/Mixture.m @@ -1490,8 +1490,9 @@ function computeComposition(obj) if obj.FLAG_VOLUME obj.v = obj.vSpecific * obj.mi; molarVolume = obj.v / N_gas; - % Update pressure [bar] if specific volume is given - obj.p = obj.equationState.getPressure(temperature, obj.v / N_gas, obj.Xi, obj.chemicalSystem) * combustiontoolbox.common.Units.Pa2bar; + % Update pressure if specific volume is given + obj.p = obj.equationState.getPressure(temperature, molarVolume, obj.Xi, obj.chemicalSystem) * combustiontoolbox.common.Units.Pa2bar; + pressure_Pa = obj.p * combustiontoolbox.common.Units.bar2Pa; else molarVolume = obj.equationState.getVolume(temperature, pressure_Pa, obj.Xi, obj.chemicalSystem); obj.v = molarVolume * N_gas;