Detailed Description

Peng-Robinson equation of state.

This class implements the Peng-Robinson equation of state (PR-EOS), a widely used cubic equation of state for describing the behaviour of real gases, particularly hydrocarbons. It accounts for non-ideal behaviour of fluids over a range of temperatures and pressures.

The equation is given in terms of the molar density \(\rho\) as:

\[ P = \frac{R T \rho}{1 - b \rho} - \frac{a \rho^2}{1 + 2b \rho - b^2 \rho^2} \]

where \(P\) is the pressure, \(T\) the temperature, \(\rho\) the molar density, \(R\) the universal gas constant, and \(a\), \(b\) are substance-specific parameters.

The parameters \(a\) and \(b\) are computed from the critical temperature \(T_c\) in Kelvin, critical pressure \(p_c\) in Pascal, and (dimensionless) acentric factor \(\omega\) as follows:

\[ a = 0.457235 \frac{R^2 T_c^2}{p_c} \]

\[ b = 0.077796 \frac{R T_c}{p_c} \]

The Peng-Robinson equation is applicable for a wide range of substances (gases and liquids), particularly hydrocarbons, at conditions ranging from subcritical to supercritical. The EOS is not suitable for solid phases or very low-temperature applications where real gases behave ideally.

All input parameters (temperature, pressure, density) are assumed to be in SI units:

Temperature \(T\) in Kelvin [K]
Pressure \(P\) in Pascal [Pa]
Mass density \(\rho\) in \([kg/m^3]\)
The resulting properties will also follow SI units.

Original source: D.-Y. Peng and D.B. Robinson, "A New Two-Constant Equation of State," Industrial & Engineering Chemistry Fundamentals, vol. 15, pp. 59-64, 1976.

Definition at line 58 of file PengRobinson.h.

#include <PengRobinson.h>

Inheritance diagram for MaterialPropertyLib::PengRobinson:

[legend]

Collaboration diagram for MaterialPropertyLib::PengRobinson:

[legend]

Public Member Functions
	PengRobinson (const double Tc, const double pc, const double omega)

PropertyDataType	value (MaterialPropertyLib::VariableArray const &variable_array, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override

PropertyDataType	dValue (MaterialPropertyLib::VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const override

Public Member Functions inherited from MaterialPropertyLib::Property
virtual	~Property ()

virtual PropertyDataType	initialValue (ParameterLib::SpatialPosition const &pos, double const t) const

virtual PropertyDataType	value () const

virtual PropertyDataType	value (VariableArray const &variable_array, VariableArray const &variable_array_prev, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

virtual PropertyDataType	dValue (VariableArray const &variable_array, VariableArray const &variable_array_prev, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

virtual PropertyDataType	d2Value (VariableArray const &variable_array, Variable const variable1, Variable const variable2, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const
	Default implementation: 2nd derivative of any constant property is zero.

virtual void	setProperties (std::vector< std::unique_ptr< Phase > > const &phases)
	Default implementation:

void	setScale (std::variant< Medium , Phase , Component * > scale)

template<typename T >
T	initialValue (ParameterLib::SpatialPosition const &pos, double const t) const

template<typename T >
T	value () const

template<typename T >
T	value (VariableArray const &variable_array, VariableArray const &variable_array_prev, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

template<typename T >
T	value (VariableArray const &variable_array, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

template<typename T >
T	dValue (VariableArray const &variable_array, VariableArray const &variable_array_prev, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

template<typename T >
T	dValue (VariableArray const &variable_array, Variable const variable, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

template<typename T >
T	d2Value (VariableArray const &variable_array, Variable const &variable1, Variable const &variable2, ParameterLib::SpatialPosition const &pos, double const t, double const dt) const

Private Attributes
const double	Tc_
	Critical temperature.

const double	pc_
	Critical pressure.

const double	omega_
	Acentric factor.

double	a_
	Parameter a (cohesion pressure)

double	b_
	Parameter b (covolume)

Additional Inherited Members
Protected Attributes inherited from MaterialPropertyLib::Property
std::string	name_

PropertyDataType	value_
	The single value of a property.

PropertyDataType	dvalue_

std::variant< Medium , Phase , Component * >	scale_

Constructor & Destructor Documentation

◆ PengRobinson()

MaterialPropertyLib::PengRobinson::PengRobinson	(	const double	Tc,
		const double	pc,
		const double	omega )

inlineexplicit

Definition at line 61 of file PengRobinson.h.

        : Tc_(Tc), pc_(pc), omega_(omega)
    {
        const double gas_constant =
            MaterialLib::PhysicalConstant::IdealGasConstant;
        a_ = 0.457235 * gas_constant * gas_constant * Tc_ * Tc_ / pc_;
        b_ = 0.077796 * gas_constant * Tc_ / pc_;
    }

References a_, b_, MaterialLib::PhysicalConstant::IdealGasConstant, pc_, and Tc_.

Member Function Documentation

◆ dValue()

PropertyDataType MaterialPropertyLib::PengRobinson::dValue	(	MaterialPropertyLib::VariableArray const &	variable_array,
		Variable const	variable,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt ) const

overridevirtual

This virtual method will compute the property derivative value based on the variables that are passed as arguments with the default implementation using empty variables array for the previous time step.

The default implementation of this method only returns the property value derivative without altering it.

Reimplemented from MaterialPropertyLib::Property.

Definition at line 57 of file PengRobinson.cpp.

{
    if (variable == Variable::temperature)
    {
        const double temperature = variable_array.temperature;
        const double epsilon = 1.e-6;
 
        MaterialPropertyLib::VariableArray perturbed = variable_array;
        // Increase temperature by +epsilon
        perturbed.temperature = temperature + epsilon;
        const double rho_plus = std::get<double>(value(perturbed, pos, t, dt));
 
        // Decrease temperature by -epsilon
        perturbed.temperature = temperature - epsilon;
        const double rho_minus = std::get<double>(value(perturbed, pos, t, dt));
 
        // Calculate the central difference quotient
        return (rho_plus - rho_minus) / (2 * epsilon);
    }
 
    if (variable == Variable::gas_phase_pressure)
    {
        const double pressure = variable_array.gas_phase_pressure;
        const double epsilon = 1.e-3;
 
        MaterialPropertyLib::VariableArray perturbed = variable_array;
        // Increase pressure by +epsilon
        perturbed.gas_phase_pressure = pressure + epsilon;
        const double rho_plus = std::get<double>(value(perturbed, pos, t, dt));
 
        // Decrease pressure by -epsilon
        perturbed.gas_phase_pressure = pressure - epsilon;
        const double rho_minus = std::get<double>(value(perturbed, pos, t, dt));
 
        // Calculate the central difference quotient
        return (rho_plus - rho_minus) / (2 * epsilon);
    }
 
    OGS_FATAL(
        "PengRobinson::dValue is implemented for derivatives with respect to "
        "gas phase pressure or temperature only.");
 
    return 0.;
}

References MaterialPropertyLib::gas_phase_pressure, MaterialPropertyLib::VariableArray::gas_phase_pressure, OGS_FATAL, MaterialPropertyLib::temperature, MaterialPropertyLib::VariableArray::temperature, and MaterialPropertyLib::Property::value().

◆ value()

PropertyDataType MaterialPropertyLib::PengRobinson::value	(	MaterialPropertyLib::VariableArray const &	variable_array,
		ParameterLib::SpatialPosition const &	pos,
		double const	t,
		double const	dt ) const

overridevirtual

This virtual method will compute the property value based on the variables that are passed as arguments with the default implementation using empty variables array for the previous time step.

Reimplemented from MaterialPropertyLib::Property.

Definition at line 22 of file PengRobinson.cpp.

{
    const double gas_constant = MaterialLib::PhysicalConstant::IdealGasConstant;
    const double pressure = variable_array.gas_phase_pressure;
    const double temperature = variable_array.temperature;
    const double molar_mass = variable_array.molar_mass;
 
    const double Tr = temperature / Tc_;  // reduced Temperature
 
    const double kappa = 0.37464 + 1.5422 * omega_ - 0.26992 * omega_ * omega_;
 
    const double alpha = boost::math::pow<2>(1 + kappa * (1 - std::sqrt(Tr)));
 
    // EOS in the form: 0 = rho^3 + z1*rho^2 + z2*rho + z3
    const double denominator =
        b_ *
        (pressure * b_ * b_ + b_ * gas_constant * temperature - a_ * alpha);
 
    const double z1 =
        (molar_mass * a_ * alpha - 3 * molar_mass * b_ * b_ * pressure -
         2 * molar_mass * gas_constant * temperature * b_) /
        denominator;
    const auto z2 = (molar_mass * molar_mass *
                     (b_ * pressure - gas_constant * temperature)) /
                    denominator;
    const auto z3 =
        (molar_mass * molar_mass * molar_mass * pressure) / denominator;
 
    MathLib::CubicSolver cubic_solver_(1., z1, z2, z3);
    return cubic_solver_.smallestPositiveRealRoot();
}

References a_, MaterialPropertyLib::alpha, b_, MaterialPropertyLib::VariableArray::gas_phase_pressure, MaterialLib::PhysicalConstant::IdealGasConstant, MaterialPropertyLib::molar_mass, MaterialPropertyLib::VariableArray::molar_mass, omega_, MathLib::CubicSolver::smallestPositiveRealRoot(), Tc_, MaterialPropertyLib::temperature, and MaterialPropertyLib::VariableArray::temperature.

Member Data Documentation

◆ a_

double MaterialPropertyLib::PengRobinson::a_

private

Parameter a (cohesion pressure)

Definition at line 87 of file PengRobinson.h.

Referenced by PengRobinson(), and value().

◆ b_

double MaterialPropertyLib::PengRobinson::b_

private

Parameter b (covolume)

Definition at line 89 of file PengRobinson.h.

Referenced by PengRobinson(), and value().

◆ omega_

const double MaterialPropertyLib::PengRobinson::omega_

private

Acentric factor.

Definition at line 85 of file PengRobinson.h.

Referenced by value().

◆ pc_

const double MaterialPropertyLib::PengRobinson::pc_

private

Critical pressure.

Definition at line 83 of file PengRobinson.h.

Referenced by PengRobinson().

◆ Tc_

const double MaterialPropertyLib::PengRobinson::Tc_

private

Critical temperature.

Definition at line 81 of file PengRobinson.h.

Referenced by PengRobinson(), and value().

The documentation for this class was generated from the following files:

MaterialLib/MPL/Properties/PengRobinson.h
MaterialLib/MPL/Properties/PengRobinson.cpp

Detailed Description

Public Member Functions

Private Attributes

Additional Inherited Members

Constructor & Destructor Documentation

◆ PengRobinson()

Member Function Documentation

◆ dValue()

◆ value()

Member Data Documentation

◆ a_

◆ b_

◆ omega_

◆ pc_

◆ Tc_