github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-09 05:15:45 -05:00
Code Issues Packages Projects Releases Wiki Activity

phases is now an enum, added (non-working) critical region VLE solver

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-15 17:13:00 +02:00
parent 3e0864e620
commit 61bc24e7b1
9 changed files with 194 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
src/AbstractState.cpp
View File

@@ -12,6 +12,7 @@
#include "Backends/Helmholtz/HelmholtzEOSBackend.h"
#include "Backends/Incompressible/IncompressibleBackend.h"
#include "Backends/Helmholtz/Fluids/FluidLibrary.h"
#include "Backends/Tabular/TabularBackends.h"
namespace CoolProp {
@@ -46,9 +47,11 @@ AbstractState * AbstractState::factory(const std::string &backend, const std::st
{
return new IncompressibleBackend(fluid_string);
}
else if (!backend.compare("BRINE"))
else if (backend.find("TTSE&") == 0)
{
throw ValueError("BRINE backend not yet implemented");
// Will throw if there is a problem with this backend
shared_ptr<AbstractState> AS(factory(backend.substr(5), fluid_string));
return new TTSEBackend(*AS.get());
}
else if (!backend.compare("TREND"))
{

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

@@ -6,7 +6,7 @@ namespace CoolProp{
void FlashRoutines::PT_flash(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.imposed_phase_index == -1) // If no phase index is imposed (see set_components function)
if (HEOS.imposed_phase_index == iphase_not_imposed) // If no phase index is imposed (see set_components function)
{
// Find the phase, while updating all internal variables possible
HEOS.T_phase_determination_pure_or_pseudopure(iP, HEOS._p);
@@ -25,6 +25,7 @@ void FlashRoutines::PT_flash(HelmholtzEOSMixtureBackend &HEOS)
void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
{
long double T = HEOS._T;
if (HEOS.is_pure_or_pseudopure)
{
// The maximum possible saturation temperature
@@ -67,8 +68,13 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
}
}
catch(std::exception &){
// We may need to polish the solution at low pressure
SaturationSolvers::saturation_T_pure_1D_P(&HEOS, HEOS._T, options);
try{
// We may need to polish the solution at low pressure
SaturationSolvers::saturation_T_pure_1D_P(&HEOS, T, options);
}
catch(std::exception &){
SaturationSolvers::saturation_critical(&HEOS, iT, T);
}
}
// Load the outputs
HEOS._phase = iphase_twophase;
@@ -263,7 +269,7 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
std::string errstring;
if (HEOS.imposed_phase_index > -1)
if (HEOS.imposed_phase_index != iphase_not_imposed)
{
// Use the phase defined by the imposed phase
HEOS._phase = HEOS.imposed_phase_index;
@@ -467,7 +473,7 @@ void FlashRoutines::HSU_P_flash_singlephase_Newton(HelmholtzEOSMixtureBackend &H
break;
}
default:
break;
throw ValueError("other variable in HSU_P_flash_singlephase_Newton is invalid");
}
//First index is the row, second index is the column
@@ -563,7 +569,7 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
{
bool saturation_called = false;
long double value;
if (HEOS.imposed_phase_index > -1)
if (HEOS.imposed_phase_index != iphase_not_imposed)
{
// Use the phase defined by the imposed phase
HEOS._phase = HEOS.imposed_phase_index;
@@ -640,7 +646,7 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
}
void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
{
if (HEOS.imposed_phase_index > -1)
if (HEOS.imposed_phase_index != iphase_not_imposed)
{
// Use the phase defined by the imposed phase
HEOS._phase = HEOS.imposed_phase_index;

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

@@ -74,7 +74,7 @@ void HelmholtzEOSMixtureBackend::set_components(std::vector<CoolPropFluid*> comp
set_excess_term();
}
imposed_phase_index = -1;
imposed_phase_index = iphase_not_imposed;
// Top-level class can hold copies of the base saturation classes,
// saturation classes cannot hold copies of the saturation classes
@@ -581,7 +581,6 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double &
case DmassSmass_INPUTS: input_pair = DmolarSmolar_INPUTS; value1 /= mm; value2 *= mm; break;
case DmassUmass_INPUTS: input_pair = DmolarUmolar_INPUTS; value1 /= mm; value2 *= mm; break;
}
}
default:
return;
@@ -589,6 +588,8 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double &
}
void HelmholtzEOSMixtureBackend::update(long input_pair, double value1, double value2 )
{
if (get_debug_level() > 0){std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)",__FILE__,__LINE__, input_pair, get_input_pair_short_desc(input_pair).c_str(), value1, value2) << std::endl;}
clear();
if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) {
@@ -1305,7 +1306,7 @@ long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv_nocache(long do
return (dOf_dtau*dConstant_ddelta-dOf_ddelta*dConstant_dtau)/(dWrt_dtau*dConstant_ddelta-dWrt_ddelta*dConstant_dtau);
}
long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv(int Of, int Wrt, int Constant)
long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant)
{
return calc_first_partial_deriv_nocache(_T, _rhomolar, Of, Wrt, Constant);
}
@@ -1465,7 +1466,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do
}
long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double p, long double rhomolar_guess)
{
int phase;
phases phase;
// Define the residual to be driven to zero
class solver_TP_resid : public FuncWrapper1D
@@ -1495,10 +1496,14 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double
solver_TP_resid resid(this,T,p);
std::string errstring;
if (imposed_phase_index > -1)
// Check if the phase is imposed
if (imposed_phase_index != iphase_not_imposed)
// Use the imposed phase index
phase = imposed_phase_index;
else
// Use the phase index in the class
phase = _phase;
if (rhomolar_guess < 0) // Not provided
{
// Calculate a guess value using SRK equation of state

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

@@ -24,10 +24,10 @@ protected:
lnK; ///< The natural logarithms of the K factors of the components
SimpleState _crit;
int imposed_phase_index;
phases imposed_phase_index;
int N; ///< Number of components
public:
HelmholtzEOSMixtureBackend(){imposed_phase_index = -1; _phase = iphase_unknown;};
HelmholtzEOSMixtureBackend(){imposed_phase_index = iphase_not_imposed; _phase = iphase_unknown;};
HelmholtzEOSMixtureBackend(std::vector<CoolPropFluid*> components, bool generate_SatL_and_SatV = true);
HelmholtzEOSMixtureBackend(std::vector<std::string> &component_names, bool generate_SatL_and_SatV = true);
virtual ~HelmholtzEOSMixtureBackend(){};
@@ -47,7 +47,7 @@ public:
bool has_melting_line(){ return is_pure_or_pseudopure && components[0]->ancillaries.melting_line.enabled();};
long double calc_melting_line(int param, int given, long double value);
int calc_phase(void){return _phase;};
phases calc_phase(void){return _phase;};
const CoolProp::SimpleState &calc_state(const std::string &state);
@@ -73,7 +73,7 @@ public:
\brief Specify the phase - this phase will always be used in calculations
@param phase_index The index from CoolProp::phases
*/
void specify_phase(int phase_index){imposed_phase_index = phase_index; _phase = phase_index;};
void specify_phase(phases phase_index){imposed_phase_index = phase_index; _phase = phase_index;};
void set_reducing_function();
void set_excess_term();
@@ -193,7 +193,7 @@ public:
\left(\frac{\partial A}{\partial B}\right)_C = \frac{\left(\frac{\partial A}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial A}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta}{\left(\frac{\partial B}{\partial \tau}\right)_\delta\left(\frac{\partial C}{\partial \delta}\right)_\tau-\left(\frac{\partial B}{\partial \delta}\right)_\tau\left(\frac{\partial C}{\partial \tau}\right)_\delta}
\f]
*/
long double calc_first_partial_deriv(int Of, int Wrt, int Constant);
long double calc_first_partial_deriv(parameters Of, parameters Wrt, parameters Constant);
/**
This version doesn't use any cached values

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

@@ -4,8 +4,77 @@
#include "MatrixMath.h"
namespace CoolProp {
void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend *HEOS, parameters ykey, long double y){
class inner_resid : public FuncWrapper1D{
public:
HelmholtzEOSMixtureBackend *HEOS;
long double T, desired_p, rhomolar_liq, calc_p, rhomolar_crit;
inner_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double desired_p)
: HEOS(HEOS), T(T), desired_p(desired_p){};
double call(double rhomolar_liq){
this->rhomolar_liq = rhomolar_liq;
HEOS->SatL->update(DmolarT_INPUTS, rhomolar_liq, T);
calc_p = HEOS->SatL->p();
std::cout << format("inner p: %0.16Lg; res: %0.16Lg", calc_p, calc_p - desired_p) << std::endl;
return calc_p - desired_p;
}
};
// Define the outer residual to be driven to zero - this is the equality of
// Gibbs function for both co-existing phases
class outer_resid : public FuncWrapper1D
{
public:
void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options){
HelmholtzEOSMixtureBackend *HEOS;
parameters ykey;
long double y;
long double r, T, rhomolar_liq, rhomolar_vap, value, p, gL, gV, rhomolar_crit;
int other;
outer_resid(HelmholtzEOSMixtureBackend *HEOS, CoolProp::parameters ykey, long double y)
: HEOS(HEOS), ykey(ykey), y(y){
rhomolar_crit = HEOS->rhomolar_critical();
};
double call(double rhomolar_vap){
this->y = y;
// Calculate the other variable (T->p or p->T) for given vapor density
if (ykey == iT){
T = y;
HEOS->SatV->update(DmolarT_INPUTS, rhomolar_vap, y);
this->p = HEOS->SatV->p();
std::cout << format("outer p: %0.16Lg",this->p) << std::endl;
inner_resid inner(HEOS, T, p);
std::string errstr2;
rhomolar_liq = Brent(inner, rhomolar_crit*1.5, rhomolar_crit*(1+1e-8), LDBL_EPSILON, 1e-10, 100, errstr2);
}
HEOS->SatL->update(DmolarT_INPUTS, rhomolar_liq, T);
HEOS->SatV->update(DmolarT_INPUTS, rhomolar_vap, T);
// Calculate the Gibbs functions for liquid and vapor
gL = HEOS->SatL->gibbsmolar();
gV = HEOS->SatV->gibbsmolar();
// Residual is difference in Gibbs function
r = gL - gV;
return this->p;
};
};
outer_resid resid(HEOS, iT, y);
double rhomolar_crit = HEOS->rhomolar_critical();
std::string errstr;
Brent(&resid, rhomolar_crit*(1-1e-8), rhomolar_crit*0.5, DBL_EPSILON, 1e-9, 20, errstr);
}
void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options)
{
// Define the residual to be driven to zero
class solver_resid : public FuncWrapper1D

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

@@ -15,8 +15,8 @@ namespace SaturationSolvers
};
struct saturation_T_pure_options{
bool use_guesses; ///< true to start off at the values specified by rhoL, rhoV
long double omega, rhoL, rhoV, pL, pV, p;
saturation_T_pure_options(){omega = _HUGE; rhoV = _HUGE; rhoL = _HUGE; rhoL = _HUGE; pV = _HUGE, pL = _HUGE;}
long double omega, rhoL, rhoV, pL, pV, p, T;
saturation_T_pure_options(){omega = _HUGE; rhoV = _HUGE; rhoL = _HUGE; rhoL = _HUGE; pV = _HUGE, pL = _HUGE; T = _HUGE;}
};
struct saturation_D_pure_options{
@@ -89,6 +89,21 @@ namespace SaturationSolvers
*/
void saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options);
/* \brief A robust but slow solver in the very-near-critical region
*
* This solver operates in the following fashion:
* 1. Using a bounded interval for rho'':[rhoc, rhoc-??], guess a value for rho''
* 2. For guessed value of rho'' and given value of T, calculate p
* 3. Using a Brent solver on the other co-existing phase (rho'), calculate the (bounded) value of rho' that yields the same pressure
* 4. Use another outer Brent solver on rho'' to enforce the same Gibbs function between liquid and vapor
* 5. Fin.
*
* @param HEOS The Helmholtz EOS backend instance to be used
* @param ykey The CoolProp::parameters key to be imposed - one of iT or iP
* @param y The value for the imposed variable
*/
void saturation_critical(HelmholtzEOSMixtureBackend *HEOS, CoolProp::parameters ykey, long double y);
void successive_substitution(HelmholtzEOSMixtureBackend &HEOS,
const long double beta,
long double T,

41
src/Tests/CoolProp-Tests.cpp
View File

@@ -793,4 +793,45 @@ TEST_CASE("Tests for solvers in P,Y flash using Water", "[flash],[PH],[PS],[PU]"
}
}
//TEST_CASE("Test that states agree with CoolProp", "[states]")
//{
// std::vector<std::string> fluids = strsplit(CoolProp::get_global_param_string("fluids_list"),',');
// for (std::size_t i = 0; i < fluids.size(); ++i)
// {
// shared_ptr<CoolProp::AbstractState> AS(CoolProp::AbstractState::factory("HEOS",fluids[i]));
//
// CoolProp::SimpleState &triple_liquid = AS->get_state("triple_liquid");
// CoolProp::SimpleState &triple_vapor = AS->get_state("triple_vapor");
//
// // See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU
// std::ostringstream ss1;
// ss1 << "Check state for triple_liquid for " << fluids[i];
// SECTION(ss1.str(),"")
// {
// std::string note = "The enthalpy and entropy are hardcoded in the fluid JSON files. They MUST agree with the values calculated by the EOS";
// AS->update(CoolProp::DmolarT_INPUTS, hs_anchor.rhomolar, hs_anchor.T);
// CAPTURE(hs_anchor.hmolar);
// CAPTURE(hs_anchor.smolar);
// double EOS_hmolar = AS->hmolar();
// double EOS_smolar = AS->smolar();
// CHECK( std::abs(EOS_hmolar - hs_anchor.hmolar) < 1e-3);
// CHECK( std::abs(EOS_smolar - hs_anchor.smolar) < 1e-3);
// }
//
// // See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU
// std::ostringstream ss2;
// ss2 << "Check state for triple_vapor for " << fluids[i];
// SECTION(ss2.str(),"")
// {
// std::string note = "The enthalpy and entropy are hardcoded in the fluid JSON files. They MUST agree with the values calculated by the EOS";
// AS->update(CoolProp::DmolarT_INPUTS, hs_anchor.rhomolar, hs_anchor.T);
// CAPTURE(hs_anchor.hmolar);
// CAPTURE(hs_anchor.smolar);
// double EOS_hmolar = AS->hmolar();
// double EOS_smolar = AS->smolar();
// CHECK( std::abs(EOS_hmolar - hs_anchor.hmolar) < 1e-3);
// CHECK( std::abs(EOS_smolar - hs_anchor.smolar) < 1e-3);
// }
// }
//}
#endif