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

Merge branch 'master' of https://github.com/coolprop/coolprop

Browse Source
This commit is contained in:
Ian Bell
2016-04-04 20:06:21 -06:00
parent 130c6ed456 396aed1b3d
commit 567cbd5c3d
6 changed files with 247 additions and 28 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
dev/buildbot/master/master.cfg
View File

@@ -77,13 +77,15 @@ slave_commons = dict(notify_on_missing=["ian.h.bell@gmail.com", "jowr@ipu.dk"],
c['slaves'] = [BuildSlave("linux-slave", pass_dict["linux-slave"], **slave_commons),
BuildSlave("OSX-slave", pass_dict["OSX-slave"], **slave_commons),
BuildSlave("windows-slave", pass_dict["windows-slave"], **slave_commons),
BuildSlave("windows-DTU-slave", pass_dict["windows-DTU-slave"], **slave_commons),
BuildSlave("OSX-IPU-slave", pass_dict["OSX-IPU-slave"], **slave_commons),
BuildSlave("OSX-IPU-slave", pass_dict["OSX-IPU-slave"], **slave_commons),
BuildSlave("Linux64-IPU-slave", pass_dict["Linux64-IPU-slave"], **slave_commons),
BuildSlave("Linux32-IPU-slave", pass_dict["Linux32-IPU-slave"], **slave_commons),
# BuildSlave("LinuxWeb-IPU-slave", pass_dict["LinuxWeb-IPU-slave"], **slave_commons),
# BuildSlave("docker-slave", pass_dict["docker-slave"], **slave_commons),
]
# The Windows builder is a dedicated machine with 4 cores.
slave_commons["max_builds"] = 1 # Should be 3, but there are problems with the upload step
c['slaves'].extend([
BuildSlave("windows-DTU-slave", pass_dict["windows-DTU-slave"], **slave_commons)
])
# 'slavePortnum' defines the TCP port to listen on for connections from slaves.
# This must match the value configured into the buildslaves (with their

23
src/Backends/Helmholtz/FlashRoutines.cpp
View File

@@ -11,28 +11,7 @@
namespace CoolProp{
template<class T> T g_RachfordRice(const std::vector<T> &z, const std::vector<T> &lnK, T beta)
{
// g function from Rachford-Rice
T summer = 0;
for (std::size_t i = 0; i < z.size(); i++)
{
T Ki = exp(lnK[i]);
summer += z[i]*(Ki-1)/(1-beta+beta*Ki);
}
return summer;
}
template<class T> T dgdbeta_RachfordRice(const std::vector<T> &z, const std::vector<T> &lnK, T beta)
{
// derivative of g function from Rachford-Rice with respect to beta
T summer = 0;
for (std::size_t i = 0; i < z.size(); i++)
{
T Ki = exp(lnK[i]);
summer += -z[i]*pow((Ki-1)/(1-beta+beta*Ki),2);
}
return summer;
}
void FlashRoutines::PT_flash_mixtures(HelmholtzEOSMixtureBackend &HEOS)
{

25
src/Backends/Helmholtz/FlashRoutines.h
View File

@@ -27,6 +27,31 @@ and not pollute the HelmholtzEOSMixtureBackend namespace
*/
class FlashRoutines{
public:
template<class T> T
static g_RachfordRice(const std::vector<T> &z, const std::vector<T> &lnK, T beta)
{
// g function from Rachford-Rice
T summer = 0;
for (std::size_t i = 0; i < z.size(); i++)
{
T Ki = exp(lnK[i]);
summer += z[i]*(Ki-1)/(1-beta+beta*Ki);
}
return summer;
}
template<class T> T
static dgdbeta_RachfordRice(const std::vector<T> &z, const std::vector<T> &lnK, T beta)
{
// derivative of g function from Rachford-Rice with respect to beta
T summer = 0;
for (std::size_t i = 0; i < z.size(); i++)
{
T Ki = exp(lnK[i]);
summer += -z[i]*pow((Ki-1)/(1-beta+beta*Ki),2);
}
return summer;
}
/// Flash for given pressure and (molar) quality
/// @param HEOS The HelmholtzEOSMixtureBackend to be used

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

@@ -3144,10 +3144,17 @@ CoolProp::CriticalState HelmholtzEOSMixtureBackend::calc_critical_point(double r
double min_eigenvalue = es.eigenvalues().real().minCoeff();
CriticalState critical;
critical.stable = (min_eigenvalue > 0); // TODO: stability analysis
critical.T = _critical.T;
critical.p = _critical.p;
critical.rhomolar = _critical.rhomolar;
if (_critical.p < 0){
critical.stable = false;
}
else{
// Otherwise we try to check stability with TPD-based analysis
StabilityRoutines::StabilityEvaluationClass stability_tester(*this);
critical.stable = stability_tester.is_stable();
}
return critical;
}

167
src/Backends/Helmholtz/VLERoutines.cpp
View File

@@ -4,6 +4,7 @@
#include "MatrixMath.h"
#include "MixtureDerivatives.h"
#include "Configuration.h"
#include "FlashRoutines.h"
namespace CoolProp {
@@ -1699,5 +1700,171 @@ void SaturationSolvers::newton_raphson_twophase::build_arrays()
error_rms = r.norm(); // Square-root (The R in RMS)
}
class RachfordRiceResidual: public FuncWrapper1DWithDeriv{
private:
const std::vector<double> &z, &lnK;
public:
RachfordRiceResidual(const std::vector<double> &z, const std::vector<double> &lnK) : z(z), lnK(lnK) {};
double call(double beta){ return FlashRoutines::g_RachfordRice(z, lnK, beta); }
double deriv(double beta){ return FlashRoutines::dgdbeta_RachfordRice(z, lnK, beta); }
};
void StabilityRoutines::StabilityEvaluationClass::trial_compositions(){
x.resize(z.size()); y.resize(z.size()); lnK.resize(z.size()); K.resize(z.size());
double g0 = 0, g1 = 0, beta = -1;
for (int i = 0; i < static_cast<int>(z.size()); ++i){
// Calculate the K-factor
lnK[i] = SaturationSolvers::Wilson_lnK_factor(HEOS, HEOS.T(), HEOS.p(), i);
K[i] = exp(lnK[i]);
g0 += z[i]*(K[i]-1); // The summation for beta = 0
g1 += z[i]*(1-1/K[i]); // The summation for beta = 1
}
// Now see what to do about the g(0) and g(1) values
// -----
//
if (g0 < 0){
beta = 0; // Assumed to be at bubble-point temperature
}
else if (g1 > 0){
beta = 1; // Assumed to be at the dew-point temperature
}
else{
// Need to iterate to find beta that makes g of Rachford-Rice zero
RachfordRiceResidual resid(z, lnK);
beta = Brent(resid, 0, 1, DBL_EPSILON, 1e-10, 100);
}
// Get the compositions from given value for beta, K, z
SaturationSolvers::x_and_y_from_K(beta, K, z, x, y);
normalize_vector(x);
normalize_vector(y);
}
void StabilityRoutines::StabilityEvaluationClass::successive_substitution(int num_steps){
// ----
// Do a few steps of successive substitution
// ----
HEOS.SatL->set_mole_fractions(x);
HEOS.SatV->set_mole_fractions(y);
HEOS.SatL->calc_reducing_state();
HEOS.SatV->calc_reducing_state();
// First use cubic as a guess for the density of liquid and vapor phases
rhomolar_liq = HEOS.SatL->solver_rho_Tp_SRK(HEOS.T(), HEOS.p(), iphase_liquid); // [mol/m^3]
rhomolar_vap = HEOS.SatV->solver_rho_Tp_SRK(HEOS.T(), HEOS.p(), iphase_gas); // [mol/m^3]
// Use Peneloux volume translation to shift liquid volume
// As in Horstmann :: doi:10.1016/j.fluid.2004.11.002
double summer_c = 0, v_SRK = 1/rhomolar_liq;
for (std::size_t i = 0; i < z.size(); ++i){
// Get the parameters for the cubic EOS
CoolPropDbl Tc = HEOS.get_fluid_constant(i, iT_critical),
pc = HEOS.get_fluid_constant(i, iP_critical),
rhomolarc = HEOS.get_fluid_constant(i, irhomolar_critical);
CoolPropDbl R = 8.3144598;
summer_c += z[i]*(0.40768*R*Tc/pc*(0.29441 - pc/(rhomolarc*R*Tc)));
}
rhomolar_liq = 1/(v_SRK - summer_c);
for (int step_count = 0; step_count < num_steps; ++step_count){
// Set the composition
HEOS.SatL->set_mole_fractions(x); HEOS.SatV->set_mole_fractions(y);
HEOS.SatL->calc_reducing_state(); HEOS.SatV->calc_reducing_state();
// Re-calculate the density
HEOS.SatL->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_liq);
rhomolar_liq = HEOS.SatL->rhomolar();
HEOS.SatV->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_vap);
rhomolar_vap = HEOS.SatV->rhomolar();
// Calculate the new K-factors from the fugacity coefficients
double g0 = 0, g1 = 0;
for (std::size_t i = 0; i < z.size(); ++i){
lnK[i] = log(HEOS.SatL->fugacity_coefficient(i)/HEOS.SatV->fugacity_coefficient(i));
K[i] = exp(lnK[i]);
g0 += z[i]*(K[i]-1); // The summation for beta = 0
g1 += z[i]*(1-1/K[i]); // The summation for beta = 1
}
// TODO: this logic for g0 < 0 and g1 > 0 seems incorrect; beta then becomes by definition in 0 <= beta <= 1
if (g0 < 0){
beta = 0; // Assumed to be at bubble-point temperature
}
else if (g1 > 0){
beta = 1; // Assumed to be at the dew-point temperature
}
else{
// Need to iterate to find beta that makes g of Rachford-Rice zero
RachfordRiceResidual resid(z, lnK);
beta = Brent(resid, 0, 1, DBL_EPSILON, 1e-10, 100);
}
// Get the compositions from given values for beta, K, z
SaturationSolvers::x_and_y_from_K(beta, K, z, x, y);
normalize_vector(x);
normalize_vector(y);
}
}
void StabilityRoutines::StabilityEvaluationClass::check_stability(){
std::vector<double> tpdL, tpdH;
_stable = true;
if (0 <= beta && beta <= 1){
// Calculate the tpd and the Gibbs energy difference (Gernert, 2014, Eqs. 20-22)
this->tpd_liq = HEOS.tangent_plane_distance(HEOS.T(), HEOS.p(), x, rhomolar_liq);
this->tpd_vap = HEOS.tangent_plane_distance(HEOS.T(), HEOS.p(), y, rhomolar_vap);
this->DELTAG_nRT = (1-beta)*tpd_liq + beta*(tpd_vap);
// If any of these cases are met, conclusively unstable, stop!
if (this->tpd_liq < 0 || this->tpd_vap < 0 || this->DELTAG_nRT < 0){
_stable = false; return;
}
}
// Ok, we aren't sure about stability, need to keep going with the full tpd analysis
// Generate light and heavy test compositions (Gernert, 2014, Eq. 23)
std::vector<double> xL(z.size()), xH(z.size());
for (std::size_t i = 0; i < z.size(); ++i){
xL[i] = z[i]*K[i];
xH[i] = z[i]/K[i];
}
normalize_vector(xL);
normalize_vector(xH);
// For each composition, use successive substitution to try to evaluate stability
for (int step_count = 0; step_count < 100; ++step_count){
// Set the composition
HEOS.SatL->set_mole_fractions(xH); HEOS.SatV->set_mole_fractions(xL);
HEOS.SatL->calc_reducing_state(); HEOS.SatV->calc_reducing_state();
// Re-calculate the density
HEOS.SatL->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_liq);
rhomolar_liq = HEOS.SatL->rhomolar();
HEOS.SatV->update_TP_guessrho(HEOS.T(), HEOS.p(), rhomolar_vap);
rhomolar_vap = HEOS.SatV->rhomolar();
// Calculate and store TPD values
double tpd_L = HEOS.tangent_plane_distance(HEOS.T(), HEOS.p(), xL, rhomolar_vap);
double tpd_H = HEOS.tangent_plane_distance(HEOS.T(), HEOS.p(), xH, rhomolar_liq);
tpdL.push_back(tpd_L); tpdH.push_back(tpd_H);
if (tpd_L < 0 || tpd_H < 0){_stable = false; return;}
// Calculate the new composition from the fugacity coefficients
for (std::size_t i = 0; i < z.size(); ++i){
// TODO: this seems weird; what fugacity coefficient is meant to be used in numerator of Gernert, 2014, Eq. 24 ?
xL[i] = z[i]*HEOS.fugacity_coefficient(i)/HEOS.SatV->fugacity_coefficient(i);
xH[i] = z[i]*HEOS.fugacity_coefficient(i)/HEOS.SatL->fugacity_coefficient(i);
}
normalize_vector(xL);
normalize_vector(xH);
}
}
} /* namespace CoolProp*/

41
src/Backends/Helmholtz/VLERoutines.h
View File

@@ -433,6 +433,45 @@ namespace SaturationSolvers
};
};
} /* namespace CoolProp*/
namespace StabilityRoutines{
/** \brief Evaluate phase stability
* Based on the work of Gernert et al., J. Chem. Thermodyn., 2014 http://dx.doi.org/10.1016/j.fluid.2014.05.012
*/
class StabilityEvaluationClass{
protected:
HelmholtzEOSMixtureBackend &HEOS;
std::vector<double> lnK, K, x, y;
const std::vector<double> &z;
double rhomolar_liq, rhomolar_vap, beta, tpd_liq, tpd_vap, DELTAG_nRT;
private:
bool _stable;
public:
StabilityEvaluationClass(HelmholtzEOSMixtureBackend &HEOS)
: HEOS(HEOS), z(HEOS.get_mole_fractions_doubleref()), rhomolar_liq(-1), rhomolar_vap(-1), beta(-1), tpd_liq(10000), tpd_vap(100000), DELTAG_nRT(10000), _stable(false) {};
/** \brief Calculate trial compositions
*/
void trial_compositions();
/** \brief Successive substitution
*/
void successive_substitution(int num_steps);
/** \brief Check stability
* 1. Check stability by looking at tpd', tpd'' and \f$ \Delta G/(nRT)\f$
* 2. Do a full TPD analysis
*/
void check_stability();
/** \brief Return best estimate for the stability of the point
*/
bool is_stable(){
trial_compositions();
successive_substitution(3);
check_stability();
return _stable;
}
};
}; /* namespace StabilityRoutines*/
}; /* namespace CoolProp*/
#endif