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

Implemented saturation P 1D solver for low pressure for some fluids like n-Propane

Browse Source 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-13 11:24:27 +02:00
parent fb23000f98
commit a84edcdb57
7 changed files with 125 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/AbstractState.cpp
View File

@@ -115,6 +115,7 @@ bool AbstractState::clear() {
this->_speed_sound.clear();
this->_hmolar.clear();
this->_smolar.clear();
this->_gibbsmolar.clear();
this->_logp.clear();
this->_logrhomolar.clear();

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

@@ -155,22 +155,29 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
options.use_logdelta = false;
double increment = 0.2;
for (double omega = 1.0; omega > 0; omega -= increment){
try{
options.omega = omega;
// Actually call the solver
SaturationSolvers::saturation_PHSU_pure(&HEOS, HEOS._p, options);
// If you get here, there was no error, all is well
break;
}
catch(std::exception &e){
if (omega < 1.1*increment){
throw;
}
}
}
try{
for (double omega = 1.0; omega > 0; omega -= increment){
try{
options.omega = omega;
// Actually call the solver
SaturationSolvers::saturation_PHSU_pure(&HEOS, HEOS._p, options);
// If you get here, there was no error, all is well
break;
}
catch(std::exception &e){
if (omega < 1.1*increment){
throw;
}
}
}
}
catch(std::exception &){
// We may need to polish the solution at low pressure
SaturationSolvers::saturation_P_pure_1D_T(&HEOS, HEOS._p, options);
}
// Load the outputs
HEOS._phase = iphase_twophase;

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

@@ -121,13 +121,15 @@ void HelmholtzEOSMixtureBackend::update_states(void)
{
CoolPropFluid &component = *(components[0]);
EquationOfState &EOS = component.EOSVector[0];
long double rho, T;
// Clear the state class
clear();
// Calculate the new enthalpy and entropy values
update(DmolarT_INPUTS, EOS.hs_anchor.rhomolar, EOS.hs_anchor.T);
EOS.hs_anchor.hmolar = hmolar();
EOS.hs_anchor.smolar = smolar();
// Clear again just to be sure
clear();
}
@@ -1764,6 +1766,9 @@ long double HelmholtzEOSMixtureBackend::calc_umolar(void)
return static_cast<long double>(_umolar);
}
else{
throw ValueError(format("phase is invalid"));
}
}
long double HelmholtzEOSMixtureBackend::calc_cvmolar(void)
{
@@ -1833,7 +1838,34 @@ long double HelmholtzEOSMixtureBackend::calc_speed_sound(void)
return static_cast<double>(_speed_sound);
}
long double HelmholtzEOSMixtureBackend::calc_gibbsmolar(void)
{
if (isTwoPhase())
{
_gibbsmolar = _Q*SatV->gibbsmolar() + (1 - _Q)*SatL->gibbsmolar();
return static_cast<long double>(_gibbsmolar);
}
else if (isHomogeneousPhase())
{
// Calculate the reducing parameters
_delta = _rhomolar/_reducing.rhomolar;
_tau = _reducing.T/_T;
// Calculate derivatives if needed, or just use cached values
long double ar = alphar();
long double a0 = alpha0();
long double dar_dDelta = dalphar_dDelta();
long double R_u = gas_constant();
// Get molar entropy
_gibbsmolar = R_u*_T*(1 + a0 + ar +_delta.pt()*dar_dDelta);
return static_cast<long double>(_gibbsmolar);
}
else{
throw ValueError(format("phase is invalid"));
}
}
long double HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(int i)
{
return exp(mixderiv_ln_fugacity_coefficient(i));

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

@@ -104,6 +104,7 @@ public:
long double calc_pressure(void);
long double calc_cvmolar(void);
long double calc_cpmolar(void);
long double calc_gibbsmolar(void);
long double calc_cpmolar_idealgas(void);
long double calc_pressure_nocache(long double T, long double rhomolar);
long double calc_smolar(void);

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

@@ -5,6 +5,52 @@
namespace CoolProp {
void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS, long double p, saturation_PHSU_pure_options &options){
// Define the residual to be driven to zero
class solver_resid : public FuncWrapper1D
{
public:
HelmholtzEOSMixtureBackend *HEOS;
long double r, p, rhomolar_liq, rhomolar_vap, value, T, gL, gV;
int other;
solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double p, long double rhomolar_liq_guess, long double rhomolar_vap_guess)
: HEOS(HEOS), p(p), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){};
double call(double T){
this->T = T;
// Recalculate the densities using the current guess values
rhomolar_liq = HEOS->SatL->solver_rho_Tp(T, p, rhomolar_liq);
rhomolar_vap = HEOS->SatV->solver_rho_Tp(T, p, rhomolar_vap);
// Set the densities in the saturation classes
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 r;
};
};
solver_resid resid(HEOS, p, options.rhoL, options.rhoV);
if (!ValidNumber(options.T)){throw ValueError("options.T is not valid in saturation_P_pure_1D_T");};
if (!ValidNumber(options.rhoL)){throw ValueError("options.rhoL is not valid in saturation_P_pure_1D_T");};
if (!ValidNumber(options.rhoV)){throw ValueError("options.rhoV is not valid in saturation_P_pure_1D_T");};
std::string errstr;
long double Tmax = std::min(options.T + 1, static_cast<long double>(HEOS->T_critical()-1e-6));
long double Tmin = std::max(options.T - 1, static_cast<long double>(HEOS->Ttriple()+1e-6));
BoundedSecant(resid, options.T, Tmin, Tmax, 0.5, 1e-11, 100, errstr);
int r =3;
}
void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, long double specified_value, saturation_PHSU_pure_options &options)
{
/*
@@ -199,9 +245,12 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l
{
throw SolutionError(format("saturation_PHSU_pure solver T < 0"));
}
if (iter > 200){
if (iter > 25){
// Set values back into the options structure for use in next solver
options.rhoL = rhoL; options.rhoV = rhoV; options.T = T;
// Error out
std::string info = get_parameter_information(specified_parameter, "short");
throw SolutionError(format("saturation_PHSU_pure solver did not converge after 200 iterations for %s=%Lg current error is %Lg", info.c_str(), specified_value, error));
throw SolutionError(format("saturation_PHSU_pure solver did not converge after 25 iterations for %s=%Lg current error is %Lg", info.c_str(), specified_value, error));
}
}
while (error > 1e-10);

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

@@ -59,7 +59,7 @@ namespace SaturationSolvers
bool use_guesses, ///< True to start off at the values specified by rhoL, rhoV, T
use_logdelta; ///< True to use partials with respect to log(delta) rather than delta
int specified_variable;
long double omega, rhoL, rhoV, pL, pV;
long double omega, rhoL, rhoV, pL, pV, T;
saturation_PHSU_pure_options(){ specified_variable = IMPOSED_INVALID_INPUT; use_guesses = true; omega = 1.0; }
};
/**
@@ -67,6 +67,16 @@ namespace SaturationSolvers
*/
void saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, long double specified_value, saturation_PHSU_pure_options &options);
/* \brief This is a backup saturation_p solver for the case where the Newton solver cannot approach closely enough the solution
*
* This is especially a problem at low pressures where the truncation error occurs, especially in the saturated vapor side
*
* @param HEOS The Helmholtz EOS backend instance to be used
* @param p Imposed pressure in kPa
* @param options Options to be passed to the function (at least T, rhoL and rhoV must be provided)
*/
void saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS, long double p, saturation_PHSU_pure_options &options);
void successive_substitution(HelmholtzEOSMixtureBackend &HEOS,
const long double beta,
long double T,