github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-04-23 03:00:17 -04:00
Code Issues Packages Projects Releases Wiki Activity

Cubic backend can now calculate critical points finally

Browse Source
This commit is contained in:
Ian Bell
2016-01-31 11:01:47 -07:00
parent ea1485bb6f
commit e4c435281a
5 changed files with 190 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files

103
src/Backends/Cubics/CubicBackend.cpp
View File

@@ -9,4 +9,107 @@ void CoolProp::AbstractCubicBackend::setup(){
}
// Now set the reducing function for the mixture
Reducing.reset(new ConstantReducingFunction(cubic->T_r, cubic->rho_r));
}
void CoolProp::AbstractCubicBackend::get_linear_reducing_parameters(double &rhomolar_r, double &T_r){
// In the case of models where the reducing temperature is not a function of composition (SRK, PR, etc.),
// we need to use an appropriate value for T_r and v_r, here we use a linear weighting
T_r = 0;
double v_r = 0;
const std::vector<double> Tc = cubic->get_Tc(), pc = cubic->get_pc();
std::vector<double> rhoc(2,10139.128); rhoc[1] = 6856.886685;
for (std::size_t i = 0; i < mole_fractions.size(); ++i){
T_r += mole_fractions[i]*Tc[i];
v_r += mole_fractions[i]/rhoc[i];
}
rhomolar_r = 1/v_r;
}
void CoolProp::AbstractCubicBackend::get_critical_point_search_radii(double &R_delta, double &R_tau)
{
// Get the starting values from the base class
CoolProp::HelmholtzEOSMixtureBackend::get_critical_point_search_radii(R_delta, R_tau);
// Now we scale them to get the appropriate
double Tr_GERGlike, rhor_GERGlike;
get_linear_reducing_parameters(rhor_GERGlike, Tr_GERGlike);
R_delta *= rhor_GERGlike/rhomolar_reducing();
R_tau *= T_reducing()/Tr_GERGlike;
}
void CoolProp::AbstractCubicBackend::get_critical_point_starting_values(double &delta0, double &tau0){
// Get the starting values from the base class
CoolProp::HelmholtzEOSMixtureBackend::get_critical_point_starting_values(delta0, tau0);
// The starting tau and delta values can be thought of as being given with the
// GERG reducing values, or tau0 = Tr_GERG/T = 0.xxx and delta0 = rho/rhor_GERG = 0.xxx
// Then we need to multiply by
//
// Tr_cubic/Tr_GERGlike & rhor_GERGlike/rhor_cubic
//
// to get shifted values:
//
// delta0 = rho/rhor_GERG*(rhor_GERGlike/rhor_cubic)
// tau0 = Tr_GERG/T*(Tr_cubic/Tr_GERGlike)
//
double Tr_GERGlike, rhor_GERGlike;
get_linear_reducing_parameters(rhor_GERGlike, Tr_GERGlike);
delta0 *= rhor_GERGlike/rhomolar_reducing();
tau0 *= T_reducing()/Tr_GERGlike;
}
CoolPropDbl CoolProp::AbstractCubicBackend::calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl> & mole_fractions, const CoolPropDbl &tau, const CoolPropDbl &delta){
bool cache_values = true;
HelmholtzDerivatives derivs = residual_helmholtz->all(*this, get_mole_fractions_ref(), cache_values);
switch (nTau){
case 0:
{
switch (nDelta){
case 0: return derivs.alphar;
case 1: return derivs.dalphar_ddelta;
case 2: return derivs.d2alphar_ddelta2;
case 3: return derivs.d3alphar_ddelta3;
case 4: return derivs.d4alphar_ddelta4;
default: throw ValueError(format("nDelta (%d) is invalid",nDelta));
}
break;
}
case 1:
{
switch (nDelta){
case 0: return derivs.dalphar_dtau;
case 1: return derivs.d2alphar_ddelta_dtau;
case 2: return derivs.d3alphar_ddelta2_dtau;
case 3: return derivs.d4alphar_ddelta3_dtau;
default: throw ValueError(format("nDelta (%d) is invalid",nDelta));
}
break;
}
case 2:
{
switch (nDelta){
case 0: return derivs.d2alphar_dtau2;
case 1: return derivs.d3alphar_ddelta_dtau2;
case 2: return derivs.d4alphar_ddelta2_dtau2;
default: throw ValueError(format("nDelta (%d) is invalid",nDelta));
}
}
case 3:
{
switch (nDelta){
case 0: return derivs.d3alphar_dtau3;
case 1: return derivs.d4alphar_ddelta_dtau3;
default: throw ValueError(format("nDelta (%d) is invalid",nDelta));
}
}
case 4:
{
switch (nDelta){
case 0: return derivs.d4alphar_dtau4;
default: throw ValueError(format("nDelta (%d) is invalid",nDelta));
}
}
default: throw ValueError(format("nTau (%d) is invalid",nTau));
}
}

76
src/Backends/Cubics/CubicBackend.h
View File

@@ -42,10 +42,10 @@ public:
bool using_mass_fractions(void){return false;};
bool using_volu_fractions(void){return false;};
void set_mole_fractions(const std::vector<CoolPropDbl> &mole_fractions){throw NotImplementedError("Mole composition has not been implemented.");};
void set_mass_fractions(const std::vector<CoolPropDbl> &mass_fractions){}; // Not implemented, but don't throw any errors
void set_mole_fractions(const std::vector<CoolPropDbl> &mole_fractions){this->mole_fractions = mole_fractions;};
void set_mass_fractions(const std::vector<CoolPropDbl> &mass_fractions){throw NotImplementedError("Mass composition has not been implemented.");};
void set_volu_fractions(const std::vector<CoolPropDbl> &volu_fractions){throw NotImplementedError("Volume composition has not been implemented.");};
const std::vector<CoolPropDbl> & get_mole_fractions(void){throw NotImplementedError("get_mole_fractions composition has not been implemented.");};
const std::vector<CoolPropDbl> & get_mole_fractions(void){ return this->mole_fractions; };
/// Calculate the gas constant in J/mol/K
CoolPropDbl calc_gas_constant(void){
@@ -59,6 +59,18 @@ public:
reducing.rhomolar = cubic->rho_r;
return reducing;
};
/// \brief Get linear mole fraction weighting of the critical molar volumes and temperatures
/// these are used in te
void get_linear_reducing_parameters(double &rhomolar, double &T);
/// Get the the starting values for the critical point evaluation routines
void get_critical_point_starting_values(double &delta0, double &tau0);
/// Get the search radius in delta and tau for the tracer, scaled appropriately for cubic
void get_critical_point_search_radii(double &R_delta, double &R_tau);
CoolPropDbl calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl> & mole_fractions, const CoolPropDbl &tau, const CoolPropDbl &delta);
};
class SRKBackend : public AbstractCubicBackend {
@@ -125,58 +137,80 @@ public:
a.alphar = cubic->alphar(tau, delta, z, 0, 0);
a.dalphar_dtau = cubic->alphar(tau, delta, z, 1, 0);
a.dalphar_ddelta = cubic->alphar(tau, delta, z, 0, 1);
a.d2alphar_dtau2 = cubic->alphar(tau, delta, z, 2, 0);
a.d2alphar_ddelta_dtau = cubic->alphar(tau, delta, z, 1, 1);
a.d2alphar_ddelta2 = cubic->alphar(tau, delta, z, 0, 2);
a.d3alphar_dtau3 = cubic->alphar(tau, delta, z, 3, 0);
a.d3alphar_ddelta_dtau2 = cubic->alphar(tau, delta, z, 2, 1);
a.d3alphar_ddelta2_dtau = cubic->alphar(tau, delta, z, 1, 2);
a.d3alphar_ddelta3 = cubic->alphar(tau, delta, z, 0, 3);
return a;
}
virtual CoolPropDbl dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 0, 0, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 1, 0, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 0, 1, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d3alphar_dxi_dTau2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 2, 0, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d3alphar_dxi_dDelta_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 1, 1, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d3alphar_dxi_dDelta2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 0, 2, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), z, 0, 0, i, j, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d3alphar_dxi_dxj_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), z, 1, 0, i, j, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d3alphar_dxi_dxj_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), z, 0, 1, i, j, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d4alphar_dxi_dTau3(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 3, 0, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d4alphar_dxi_dDelta2_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 1, 2, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d4alphar_dxi_dDelta_dTau2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 2, 1, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d4alphar_dxi_dDelta3(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d_alphar_dxi(HEOS.tau(), HEOS.delta(), z, 0, 3, i, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d4alphar_dxi_dxj_dTau2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), z, 2, 0, i, j, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d4alphar_dxi_dxj_dDelta_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), z, 1, 1, i, j, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d4alphar_dxi_dxj_dDelta2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
return 0;
const std::vector<double> z = std::vector<double>(HEOS.get_mole_fractions().begin(), HEOS.get_mole_fractions().end());
return ACB->get_cubic()->d2_alphar_dxidxj(HEOS.tau(), HEOS.delta(), z, 0, 2, i, j, xN_flag==XN_INDEPENDENT);
}
virtual CoolPropDbl d3alphardxidxjdxk(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){
return 0;
@@ -184,4 +218,4 @@ public:
};
} /* namespace CoolProp */
#endif /* IF97BACKEND_H_ */
#endif /* CUBICBACKEND_H_ */

5
src/Backends/Cubics/GeneralizedCubic.h
View File

@@ -99,6 +99,11 @@ public:
double psi_plus(double delta, const std::vector<double> &x, std::size_t idelta);
double d_psi_plus_dxi(double delta, const std::vector<double> &x, std::size_t idelta, std::size_t i, bool xN_independent);
double d2_psi_plus_dxidxj(double delta, const std::vector<double> &x, std::size_t idelta, std::size_t i, std::size_t j, bool xN_independent);
/// Accessor to return vector of critical tempertures
const std::vector<double> & get_Tc(){return Tc;}
/// Accessor to return vector of critical pressures
const std::vector<double> & get_pc(){return pc;}
};
class PengRobinson : public AbstractCubic

24
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp
View File

@@ -3222,7 +3222,7 @@ public:
this->get_tau_delta(theta, tau, delta, tau_new, delta_new);
// Stop if bounds are exceeded
if (p_MPa > 500 || delta_new > 5){
if (p_MPa > 500){
break;
}
@@ -3256,16 +3256,28 @@ std::vector<CoolProp::CriticalState> HelmholtzEOSMixtureBackend::calc_all_critic
specify_phase(iphase_gas);
std::string errstr;
OneDimObjective resid_L0(*this, 0.5);
double tau0 = 0.66;
double delta0 = _HUGE, tau0 = _HUGE;
get_critical_point_starting_values(delta0, tau0);
OneDimObjective resid_L0(*this, delta0);
resid_L0.call(tau0);
if (resid_L0.deriv(tau0) > 0){
tau0 += 0.1;
tau0 *= 1.1;
}
double tau_L0 = Halley(resid_L0, tau0, 1e-10, 100, errstr);
L0CurveTracer tracer(*this, tau_L0, 0.5);
double T = T_reducing()/tau_L0;
double T0 = T_reducing()/tau0;
double rho = delta0*rhomolar_reducing();
L0CurveTracer tracer(*this, tau_L0, delta0);
double R_delta = 0, R_tau = 0;
get_critical_point_search_radii(R_delta, R_tau);
tracer.R_delta_tracer = R_delta;
tracer.R_tau_tracer = R_tau;
tracer.trace();
// Reset phase to previous value
_phase = old_phase;
return tracer.critical_points;

10
src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h
View File

@@ -99,6 +99,14 @@ public:
CriticalState calc_critical_point(double rho0, double T0);
std::vector<CoolProp::CriticalState> calc_all_critical_points();
virtual void get_critical_point_starting_values(double &delta0, double &tau0){
delta0 = 0.5; // The value of delta where we start searching for crossing with Lstar=0 contour
tau0 = 0.66; // The value of tau where we start searching at delta=delta0
}
/// Get the search radius in delta and tau for the tracer
virtual void get_critical_point_search_radii(double &R_delta, double &R_tau){
R_delta = 0.1; R_tau = 0.1;
}
/// Calculate tangent plane distance for given trial composition w
double calc_tangent_plane_distance(const double T, const double p, const std::vector<double> &w, const double rhomolar_guess);
@@ -326,7 +334,7 @@ public:
std::vector<std::string> calc_fluid_names(void);
void calc_all_alphar_deriv_cache(const std::vector<CoolPropDbl> &mole_fractions, const CoolPropDbl &tau, const CoolPropDbl &delta);
CoolPropDbl calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl> & mole_fractions, const CoolPropDbl &tau, const CoolPropDbl &delta);
virtual CoolPropDbl calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<CoolPropDbl> & mole_fractions, const CoolPropDbl &tau, const CoolPropDbl &delta);
/**
\brief Take derivatives of the ideal-gas part of the Helmholtz energy, don't use any cached values, or store any cached values