Update: generalize EquationState, EquationsStateIdealGas, Mixture, and Species classes to be compatible with non-ideal EoS

+combustiontoolbox/+core/@EquationState/EquationState.m

@@ -4,17 +4,142 @@
     % 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, 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, molarVolume, varargin)
+        pressure = getPressure(obj, temperature, 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, molarFractions, chemicalSystem, varargin)
+
+        % Compute thermodynamic departure functions
+        [heatCapacityPressureDeparture, enthalpyDeparture, entropyDeparture] = getDepartureFunctions(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)
+
+        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

+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,66 @@
             % 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 [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.
+            %
+            % 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

+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,13 @@ function computeComposition(obj)
             % Compute volume [m3]
             if obj.FLAG_VOLUME
                 obj.v = obj.vSpecific * obj.mi;
-                % 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;
+                molarVolume = obj.v / N_gas;
+                % 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
-                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 +1504,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 +1529,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 +1554,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 +1580,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;