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

Basic algorithm is working for P,Y for vapor - more checking needed, but basically works

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-07 01:50:09 +02:00
parent 252ddb1dd3
commit a5bf52b52f
6 changed files with 193 additions and 86 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -116,14 +116,14 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.is_pure_or_pseudopure)
{
if (HEOS.components[0]->pEOS->pseudo_pure){
// It is a psedo-pure mixture
throw NotImplementedError("PQ_flash not implemented for pseudo-pure fluids");
throw NotImplementedError("PQ_flash not implemented for pseudo-pure fluids yet");
}
else{
// ------------------
// It is a pure fluid
// ------------------
// Set some imput options
SaturationSolvers::saturation_PHSU_pure_options options;
@@ -371,7 +371,8 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
throw NotImplementedError("PHSU_D_flash not ready for mixtures");
}
}
void FlashRoutines::HSU_P_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, int other, long double T0, long double rhomolar0)
void FlashRoutines::HSU_P_flash_singlephase_Newton(HelmholtzEOSMixtureBackend &HEOS, int other, long double T0, long double rhomolar0)
{
double A[2][2], B[2][2];
long double y;
@@ -463,8 +464,61 @@ void FlashRoutines::HSU_P_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, in
HEOS.update(DmolarT_INPUTS, rhoc*delta, Tc/tau);
}
void FlashRoutines::HSU_P_flash_singlephase_Brent(HelmholtzEOSMixtureBackend &HEOS, int other, long double value, long double Tmin, long double Tmax)
{
if (!ValidNumber(HEOS._p)){throw ValueError("value for p in HSU_P_flash_singlephase_Brent is invalid");};
if (!ValidNumber(value)){throw ValueError("value for other in HSU_P_flash_singlephase_Brent is invalid");};
class solver_resid : public FuncWrapper1D
{
public:
HelmholtzEOSMixtureBackend *HEOS;
long double r, eos, p, value, T, rhomolar;
int other;
int iter;
long double r0, r1, T1, T0;
solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double p, long double value, int other) : HEOS(HEOS), p(p), value(value), other(other){iter = 0;};
double call(double T){
this->T = T;
// Run the solver with T,P as inputs;
HEOS->update(PT_INPUTS, p, T);
// Get the value of the desired variable
eos = HEOS->keyed_output(other);
// Difference between the two is to be driven to zero
r = eos - value;
// Store values for later use if there are errors
if (iter == 0){
r0 = r; T0 = T;
}
else if (iter == 1){
r1 = r; T1 = T;
}
else{
r0 = r1; T0 = T1;
r1 = r; T1 = T;
}
iter++;
std::cout << format("%g %g\n",T,r);
return r;
};
};
solver_resid resid(&HEOS, HEOS._p, value, other);
std::string errstr;
Brent(resid, Tmin, Tmax, DBL_EPSILON, 1e-12, 100, errstr);
}
// P given and one of H, S, or U
void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
{
bool saturation_called = false;
long double value;
if (HEOS.imposed_phase_index > -1)
{
// Use the phase defined by the imposed phase
@@ -474,31 +528,57 @@ void FlashRoutines::HSU_P_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
{
if (HEOS.is_pure_or_pseudopure)
{
// Find the phase, while updating all internal variables possible
switch (other)
{
case iSmolar:
HEOS.p_phase_determination_pure_or_pseudopure(iSmolar, HEOS._smolar); break;
value = HEOS.smolar(); HEOS.p_phase_determination_pure_or_pseudopure(iSmolar, value, saturation_called); break;
case iHmolar:
HEOS.p_phase_determination_pure_or_pseudopure(iHmolar, HEOS._hmolar); break;
value = HEOS.hmolar(); HEOS.p_phase_determination_pure_or_pseudopure(iHmolar, value, saturation_called); break;
case iUmolar:
HEOS.p_phase_determination_pure_or_pseudopure(iUmolar, HEOS._umolar); break;
value = HEOS.umolar(); HEOS.p_phase_determination_pure_or_pseudopure(iUmolar, value, saturation_called); break;
default:
throw ValueError(format("Input for other [%s] is invalid", get_parameter_information(other, "long").c_str()));
}
if (HEOS.isHomogeneousPhase())
{
// Now we use the single-phase solver to find T,rho given P,Y using a
// bounded 1D solver by adjusting T and using given value of p
long double Tmin, Tmax;
switch(HEOS._phase)
{
case iphase_gas:
{
Tmax = HEOS.Tmax();
if (saturation_called){ Tmin = HEOS.SatV->T() + 1e-2;}else{Tmin = HEOS._TVanc.pt() + 0.5;}
break;
}
case iphase_liquid:
{
if (saturation_called){ Tmax = HEOS.SatL->T() + 1e-2;}else{Tmax = HEOS._TLanc.pt() - 0.5;}
Tmin = HEOS.Tmin() + 1; // or melting curve data
break;
}
case iphase_supercritical:
{
Tmax = HEOS.Tmax();
Tmin = HEOS.Tmin();
break;
}
default:
{ throw ValueError(format("Not a valid homogeneous state")); }
}
HSU_P_flash_singlephase_Brent(HEOS, other, value, Tmin, Tmax);
HEOS._Q = -1;
}
else{throw ValueError("Can't be non homogeneous phase at this point");}
}
else
{
throw NotImplementedError("HSU_P_flash does not support mixtures (yet)");
}
}
if (HEOS.isHomogeneousPhase() && !ValidNumber(HEOS._T))
{
long double T0 = 300, rho0 = 100;
HSU_P_flash_singlephase(HEOS, other, T0, rho0);
HEOS._Q = -1;
}
}
}
void FlashRoutines::DHSU_T_flash(HelmholtzEOSMixtureBackend &HEOS, int other)
{

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

@@ -54,7 +54,15 @@ public:
/// @param other The index for the other input from CoolProp::parameters; allowed values are iHmolar, iSmolar, iUmolar
/// @param T0 The initial guess value for the temperature [K]
/// @param rho0 The initial guess value for the density [mol/m^3]
static void HSU_P_flash_singlephase(HelmholtzEOSMixtureBackend &HEOS, int other, long double T0, long double rhomolar0);
static void HSU_P_flash_singlephase_Newton(HelmholtzEOSMixtureBackend &HEOS, int other, long double T0, long double rhomolar0);
/// The single-phase flash routine for the pairs (P,H), (P,S), and (P,U). Similar analysis is needed
/// @param HEOS The HelmholtzEOSMixtureBackend to be used
/// @param other The index for the other input from CoolProp::parameters; allowed values are iHmolar, iSmolar, iUmolar
/// @param value The value of the other input
/// @param Tmin The lower temperature limit [K]
/// @param Tmax The higher temperature limit [K]
static void HSU_P_flash_singlephase_Brent(HelmholtzEOSMixtureBackend &HEOS, int other, long double value, long double Tmin, long double Tmax);
/// A generic flash routine for the pairs (D,P), (D,H), (D,S), and (D,U). Similar analysis is needed
/// @param HEOS The HelmholtzEOSMixtureBackend to be used

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

@@ -590,11 +590,12 @@ long double HelmholtzEOSMixtureBackend::calc_dCvirial_dT()
long double dtau_dT =-get_reducing().T/pow(_T,2);
return 1/pow(get_reducing().rhomolar,2)*calc_alphar_deriv_nocache(1,2,mole_fractions,_tau,1e-12)*dtau_dT;
}
void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int other, long double value)
void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int other, long double value, bool &saturation_called)
{
/*
Determine the phase given p and one other state variable
*/
saturation_called = false;
// Reference declaration to save indexing
CoolPropFluid &component = *(components[0]);
@@ -754,47 +755,47 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
}
break;
}
default:
{
// If it is not density, update the saturation states, needed to calculate other inputs
SatV->update(DmolarT_INPUTS, rho_vap, _T);
SatL->update(DmolarT_INPUTS, rho_liq, _T);
//default:
//{
// // If it is not density, update the saturation states, needed to calculate other inputs
// SatV->update(DmolarT_INPUTS, _rhoVanc, _TVanc);
// SatL->update(DmolarT_INPUTS, _rhoLanc, _TLanc);
switch (other)
{
case iSmolar:
{
if (value > SatV->calc_smolar()){
this->_phase = iphase_gas; return;
}
if (value < SatL->calc_smolar()){
this->_phase = iphase_liquid; return;
}
break;
}
case iHmolar:
{
if (value > SatV->calc_hmolar()){
this->_phase = iphase_gas; return;
}
else if (value < SatL->calc_hmolar()){
this->_phase = iphase_liquid; return;
}
}
case iUmolar:
{
if (value > SatV->calc_umolar()){
this->_phase = iphase_gas; return;
}
else if (value < SatL->calc_umolar()){
this->_phase = iphase_liquid; return;
}
break;
}
default:
throw ValueError(format("invalid input for other to T_phase_determination_pure_or_pseudopure"));
}
}
// switch (other)
// {
// case iSmolar:
// {
// if (value > SatV->calc_smolar()){
// this->_phase = iphase_gas; return;
// }
// if (value < SatL->calc_smolar()){
// this->_phase = iphase_liquid; return;
// }
// break;
// }
// case iHmolar:
// {
// if (value > SatV->calc_hmolar()){
// this->_phase = iphase_gas; return;
// }
// else if (value < SatL->calc_hmolar()){
// this->_phase = iphase_liquid; return;
// }
// }
// case iUmolar:
// {
// if (value > SatV->calc_umolar()){
// this->_phase = iphase_gas; return;
// }
// else if (value < SatL->calc_umolar()){
// this->_phase = iphase_liquid; return;
// }
// break;
// }
// default:
// throw ValueError(format("invalid input for other to T_phase_determination_pure_or_pseudopure"));
// }
//}
}
}
@@ -803,39 +804,36 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
// Actually have to use saturation information sadly
// For the given pressure, find the saturation state
// Run the saturation routines to determine the saturation densities and pressures
FlashRoutines::PQ_flash(*this);
HelmholtzEOSMixtureBackend HEOS(components);
HEOS._p = this->_p;
HEOS._Q = 0; // ?? What is the best to do here? Doesn't matter for our purposes since pure fluid
FlashRoutines::PQ_flash(HEOS);
// We called the saturation routines, so HEOS.SatL and HEOS.SatV are now updated
// with the saturated liquid and vapor values, which can therefore be used in
// the other solvers
saturation_called = true;
long double Q;
// if (other == iT)
// {
// if (value > 100*DBL_EPSILON + HEOS.SatL->p()){
// this->_phase = iphase_liquid; _Q = -1000; return;
// }
// else if (value < HEOS.SatV->p()-100*DBL_EPSILON){
// this->_phase = iphase_gas; _Q = 1000; return;
// }
// else{
// throw ValueError(format("subcrit T, funny p"));
// }
// }
switch (other)
{
case iDmolar:
Q = (1/value-1/SatL->rhomolar())/(1/SatV->rhomolar()-1/SatL->rhomolar()); break;
Q = (1/value-1/HEOS.SatL->rhomolar())/(1/HEOS.SatV->rhomolar()-1/HEOS.SatL->rhomolar()); break;
case iSmolar:
Q = (value - SatL->smolar())/(SatV->smolar() - SatL->smolar()); break;
Q = (value - HEOS.SatL->smolar())/(HEOS.SatV->smolar() - HEOS.SatL->smolar()); break;
case iHmolar:
Q = (value - SatL->hmolar())/(SatV->hmolar() - SatL->hmolar()); break;
Q = (value - HEOS.SatL->hmolar())/(HEOS.SatV->hmolar() - HEOS.SatL->hmolar()); break;
case iUmolar:
Q = (value - SatL->umolar())/(SatV->umolar() - SatL->umolar()); break;
Q = (value - HEOS.SatL->umolar())/(HEOS.SatV->umolar() - HEOS.SatL->umolar()); break;
default:
throw ValueError(format("bad input for other"));
}
// Update the states
this->SatL->update(DmolarT_INPUTS, HEOS.SatL->rhomolar(), HEOS.SatL->T());
this->SatV->update(DmolarT_INPUTS, HEOS.SatV->rhomolar(), HEOS.SatV->T());
if (Q < -100*DBL_EPSILON){
this->_phase = iphase_liquid; _Q = -1000; return;
this->_phase = iphase_liquid; _Q = -1000; return;
}
else if (Q > 1+100*DBL_EPSILON){
this->_phase = iphase_gas; _Q = 1000; return;
@@ -843,10 +841,11 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
else{
this->_phase = iphase_twophase;
}
_Q = Q;
// Load the outputs
_p = _Q*SatV->p() + (1-_Q)*SatL->p();
_rhomolar = 1/(_Q/SatV->rhomolar() + (1-_Q)/SatL->rhomolar());
_p = _Q*HEOS.SatV->p() + (1-_Q)*HEOS.SatL->p();
_rhomolar = 1/(_Q/HEOS.SatV->rhomolar() + (1-_Q)/HEOS.SatL->rhomolar());
return;
}
else if (_p < components[0]->pEOS->ptriple)
@@ -856,6 +855,8 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot
}
void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int other, long double value)
{
if (!ValidNumber(value)){
throw ValueError(format("value to T_phase_determination_pure_or_pseudopure is invalid"));};
// T is known, another input P, T, H, S, U is given (all molar)
if (_T < _crit.T)
{
@@ -958,10 +959,10 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot
if (other == iP)
{
if (value > 100*DBL_EPSILON + HEOS.SatL->p()){
if (value > HEOS.SatL->p()*(100*DBL_EPSILON + 1)){
this->_phase = iphase_liquid; _Q = -1000; return;
}
else if (value < HEOS.SatV->p()-100*DBL_EPSILON){
else if (value < HEOS.SatV->p()*(1 - 100*DBL_EPSILON)){
this->_phase = iphase_gas; _Q = 1000; return;
}
else{

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

@@ -36,6 +36,8 @@ public:
friend class FlashRoutines; // Allows the static methods in the FlashRoutines class to have access to all the protected members and methods of this class
friend class TransportRoutines; // Allows the static methods in the TransportRoutines class to have access to all the protected members and methods of this class
//friend class MixtureDerivatives; //
// Helmholtz EOS backend uses mole fractions
@@ -196,7 +198,7 @@ public:
// ***************************************************************
// ***************************************************************
void T_phase_determination_pure_or_pseudopure(int other, long double value);
void p_phase_determination_pure_or_pseudopure(int other, long double value);
void p_phase_determination_pure_or_pseudopure(int other, long double value, bool &saturation_called);
void DmolarP_phase_determination();

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

@@ -24,7 +24,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS->SatL,
SatV = HEOS->SatV;
long double T, rhoL,rhoV;
long double T, rhoL, rhoV, pL, pV;
long double deltaL=0, deltaV=0, tau=0, error;
int iter=0, specified_parameter;
@@ -79,8 +79,8 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
SatL->update(DmolarT_INPUTS, rhoL, T);
SatV->update(DmolarT_INPUTS, rhoV, T);
long double pL = SatL->p();
long double pV = SatV->p();
pL = SatL->p();
pV = SatV->p();
// These derivatives are needed for both cases
long double alpharL = SatL->alphar();
@@ -205,6 +205,8 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
}
}
while (error > 1e-10);
HEOS->SatL->update(DmolarT_INPUTS, rhoL, T);
HEOS->SatV->update(DmolarT_INPUTS, rhoV, T);
}
void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long double rhomolar, saturation_D_pure_options &options)
{

16
src/main.cxx
View File

@@ -18,7 +18,7 @@ using namespace CoolProp;
#endif
#include "SpeedTest.h"
#include "HumidAirProp.h"
#include "CoolPropLib.h"
//#include "CoolPropLib.h"
#include "crossplatform_shared_ptr.h"
@@ -59,7 +59,21 @@ struct element
int main()
{
#if 1
double T0 = 295;
double dT = 0.010011;
double pV = CoolProp::PropsSI("P","T",T0,"Q",1,"Water");
double hV = CoolProp::PropsSI("Hmolar","T",T0 + dT,"P",pV,"Water");
std::cout << get_global_param_string("errstring");
double TV = CoolProp::PropsSI("T","Hmolar",hV,"P",pV,"Water");
std::cout << get_global_param_string("errstring");
double ahV = CoolProp::PropsSI("P", "Q", 0, "T", 373.124, "Water");
std::cout << get_global_param_string("parameter_list") << std::endl;
#endif
#if 0
double pm = CoolProp::PropsSI("pmax","T",-100,"Q",0.0000000000000000e+00,"REFPROP::Water");
double Tm = PropsSI("Tmax","T",-100,"Q",0.0000000000000000e+00,"REFPROP::Water");
double Tt = PropsSI("P","T",-100,"Q",0.0000000000000000e+00,"AceticAcid");