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
2015-04-25 14:21:56 -06:00
parent 1f91449303 b9bd08bc93
commit cf105f07c0
27 changed files with 380 additions and 108 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -262,34 +262,34 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS)
else{
// Pseudo-pure fluid
CoolPropDbl rhoLanc = _HUGE, rhoVanc = _HUGE, rhoLsat = _HUGE, rhoVsat = _HUGE;
CoolPropDbl psatLanc = HEOS.components[0].ancillaries.pL.evaluate(HEOS._T); // These ancillaries are used explicitly
CoolPropDbl psatVanc = HEOS.components[0].ancillaries.pV.evaluate(HEOS._T); // These ancillaries are used explicitly
try{
if (std::abs(HEOS._Q) < DBL_EPSILON){
HEOS._p = HEOS.components[0].ancillaries.pL.evaluate(HEOS._T); // These ancillaries are used explicitly
rhoLanc = HEOS.components[0].ancillaries.rhoL.evaluate(HEOS._T);
HEOS.SatL->update_TP_guessrho(HEOS._T, HEOS._p, rhoLanc);
HEOS._rhomolar = HEOS.SatL->rhomolar();
}
else if (std::abs(HEOS._Q - 1) < DBL_EPSILON){
HEOS._p = HEOS.components[0].ancillaries.pV.evaluate(HEOS._T); // These ancillaries are used explicitly
rhoVanc = HEOS.components[0].ancillaries.rhoV.evaluate(HEOS._T);
HEOS.SatV->update_TP_guessrho(HEOS._T, HEOS._p, rhoVanc);
HEOS._rhomolar = HEOS.SatV->rhomolar();
}
else{
CoolPropDbl psatLanc = HEOS.components[0].ancillaries.pL.evaluate(HEOS._T); // These ancillaries are used explicitly
CoolPropDbl psatVanc = HEOS.components[0].ancillaries.pV.evaluate(HEOS._T); // These ancillaries are used explicitly
HEOS._p = HEOS._Q*psatVanc + (1-HEOS._Q)*psatLanc;
HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar());
HEOS.SatL->update(DmolarT_INPUTS, rhoLsat, HEOS._T);
HEOS.SatV->update(DmolarT_INPUTS, rhoVsat, HEOS._T);
}
if (!ValidNumber(rhoLanc) || !ValidNumber(rhoVanc))
{
throw ValueError("pseudo-pure failed");
}
HEOS.SatL->update_TP_guessrho(HEOS._T, psatLanc, rhoLanc);
HEOS.SatV->update_TP_guessrho(HEOS._T, psatVanc, rhoVanc);
if (!ValidNumber(rhoLsat) || !ValidNumber(rhoVsat) ||
std::abs(rhoLsat/rhoLanc-1) > 0.5 || std::abs(rhoVanc/rhoVsat-1) > 0.5)
{
throw ValueError("pseudo-pure failed");
}
try{
}
catch (...){
// Near the critical point, the behavior is not very nice, so we will just use the ancillary
rhoLsat = rhoLanc;
rhoVsat = rhoVanc;
}
HEOS._p = HEOS._Q*psatVanc + (1-HEOS._Q)*psatLanc;
HEOS._rhomolar = 1/(HEOS._Q/HEOS.SatV->rhomolar() + (1 - HEOS._Q)/HEOS.SatL->rhomolar());
HEOS.SatL->update(DmolarT_INPUTS, rhoLsat, HEOS._T);
HEOS.SatV->update(DmolarT_INPUTS, rhoVsat, HEOS._T);
}
// Load the outputs
HEOS._phase = iphase_twophase;

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

@@ -962,9 +962,6 @@ void HelmholtzEOSMixtureBackend::update_with_guesses(CoolProp::input_pairs input
default:
throw ValueError(format("This pair of inputs [%s] is not yet supported", get_input_pair_short_desc(input_pair).c_str()));
}
post_update();
}
@@ -1407,6 +1404,19 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot
else if (value > rho_liq){
this->_phase = iphase_liquid; return;
}
else{
_phase = iphase_liquid;
_Q = -1000;
update_DmolarT_direct(value, _T);
CoolPropDbl pL = components[0].ancillaries.pL.evaluate(_T);
if (_p > pL*1.05){
this->_phase = iphase_liquid; _Q = -1000; return;
}
else{
_phase = iphase_unknown;
_p = _HUGE;
}
}
break;
}
default:
@@ -1824,7 +1834,7 @@ CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl
phases phase;
// Define the residual to be driven to zero
class solver_TP_resid : public FuncWrapper1DWithDeriv
class solver_TP_resid : public FuncWrapper1DWithTwoDerivs
{
public:
HelmholtzEOSMixtureBackend *HEOS;
@@ -1843,6 +1853,10 @@ CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl
// dp/drho|T / pspecified
return R_u*T*(1+2*delta*HEOS->dalphar_dDelta()+pow(delta, 2)*HEOS->d2alphar_dDelta2())/p;
};
double second_deriv(double rhomolar){
// d2p/drho2|T / pspecified
return R_u*T/rhomolar*(2*delta*HEOS->dalphar_dDelta() + 4*pow(delta, 2)*HEOS->d2alphar_dDelta2() + pow(delta, 3)*HEOS->calc_d3alphar_dDelta3())/p;
};
};
solver_TP_resid resid(this,T,p);
std::string errstring;
@@ -1871,10 +1885,17 @@ CoolPropDbl HelmholtzEOSMixtureBackend::solver_rho_Tp(CoolPropDbl T, CoolPropDbl
// It's liquid at subcritical pressure, we can use ancillaries as a backup
else if (phase == iphase_liquid)
{
double rhomolar;
CoolPropDbl _rhoLancval = static_cast<CoolPropDbl>(components[0].ancillaries.rhoL.evaluate(T));
// Next we try with a Brent method bounded solver since the function is 1-1
double rhomolar = Brent(resid, _rhoLancval*0.9, _rhoLancval*1.3, DBL_EPSILON,1e-8,100,errstring);
if (!ValidNumber(rhomolar)){throw ValueError();}
try{
// First we try with Halley's method starting at saturated liquid
rhomolar = Halley(resid, _rhoLancval, 1e-16, 100, errstring);
}
catch(std::exception &){
// Next we try with a Brent method bounded solver since the function is 1-1
rhomolar = Brent(resid, _rhoLancval*0.9, _rhoLancval*1.3, DBL_EPSILON,1e-8,100,errstring);
if (!ValidNumber(rhomolar)){throw ValueError();}
}
return rhomolar;
}
else if (phase == iphase_supercritical_liquid){
@@ -2034,7 +2055,15 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_hmolar(void)
if (get_debug_level()>=50) std::cout << format("HelmholtzEOSMixtureBackend::calc_hmolar: 2phase: %d T: %g rhomomolar: %g", isTwoPhase(), _T, _rhomolar) << std::endl;
if (isTwoPhase())
{
_hmolar = _Q*SatV->hmolar() + (1 - _Q)*SatL->hmolar();
if (std::abs(_Q) < DBL_EPSILON){
_hmolar = SatL->hmolar();
}
else if (std::abs(_Q-1) < DBL_EPSILON){
_hmolar = SatV->hmolar();
}
else{
_hmolar = _Q*SatV->hmolar() + (1 - _Q)*SatL->hmolar();
}
return static_cast<CoolPropDbl>(_hmolar);
}
else if (isHomogeneousPhase())
@@ -2079,7 +2108,15 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_smolar(void)
{
if (isTwoPhase())
{
_smolar = _Q*SatV->smolar() + (1 - _Q)*SatL->smolar();
if (std::abs(_Q) < DBL_EPSILON){
_smolar = SatL->smolar();
}
else if (std::abs(_Q-1) < DBL_EPSILON){
_smolar = SatV->smolar();
}
else{
_smolar = _Q*SatV->smolar() + (1 - _Q)*SatL->smolar();
}
return static_cast<CoolPropDbl>(_smolar);
}
else if (isHomogeneousPhase())
@@ -2122,7 +2159,15 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_umolar(void)
{
if (isTwoPhase())
{
_umolar = _Q*SatV->umolar() + (1 - _Q)*SatL->umolar();
if (std::abs(_Q) < DBL_EPSILON){
_umolar = SatL->umolar();
}
else if (std::abs(_Q-1) < DBL_EPSILON){
_umolar = SatV->umolar();
}
else{
_umolar = _Q*SatV->umolar() + (1 - _Q)*SatL->umolar();
}
return static_cast<CoolPropDbl>(_umolar);
}
else if (isHomogeneousPhase())

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

@@ -913,8 +913,10 @@ void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend &HE
shared_ptr<HelmholtzEOSMixtureBackend> SatL = HEOS.SatL,
SatV = HEOS.SatV;
CoolProp::SimpleState &crit = HEOS.get_components()[0].crit;
CoolPropDbl rhoL = _HUGE, rhoV = _HUGE, error = 999, DeltavL, DeltavV, pL, pV, p;
int iter=0, small_step_count = 0;
CoolPropDbl rhoL = _HUGE, rhoV = _HUGE, error = 999, DeltavL, DeltavV, pL, pV, p, last_error;
int iter = 0,
small_step_count = 0,
backwards_step_count = 0; // Counter for the number of times you have taken a step that increases error
try
{
@@ -1001,7 +1003,7 @@ void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend &HE
if (rhoV > crit.rhomolar){
rhoV = 0.99*crit.rhomolar;
}
last_error = _HUGE;
SatL->update_DmolarT_direct(rhoL, T);
SatV->update_DmolarT_direct(rhoV, T);
if (get_debug_level() > 5){ std::cout << format("[Maxwell] starting T: %0.16Lg rhoL: %Lg rhoV: %Lg pL: %Lg pV: %g\n", T, rhoL, rhoV, SatL->p(), SatV->p());}
@@ -1063,13 +1065,19 @@ void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend &HE
if (std::abs(DeltavL*rhoL) < 10*DBL_EPSILON || std::abs(DeltavV*rhoV) < 10*DBL_EPSILON){
small_step_count++;
}
// If you are not continuing to march towards the solution, after a couple of times, stop
// This is especially a problem for water
if (std::abs(error) > std::abs(last_error)){
backwards_step_count++;
}
iter++;
last_error = error;
if (iter > 30){
throw SolutionError(format("Maxwell solver did not converge after 30 iterations; rhoL: %0.16Lg rhoV: %0.16Lg error: %Lg dvL/vL: %Lg dvV/vV: %Lg pL: %Lg pV: %Lg\n", rhoL, rhoV, error, DeltavL/vL, DeltavV/vV, pL, pV));
}
}
while ((SatL->p() < 0) || (error > 1e-5 && small_step_count < 4));
while ((SatL->p() < 0) || (error > 1e-10 && small_step_count < 4 && backwards_step_count < 6));
if (get_debug_level() > 5){ std::cout << format("[Maxwell] pL: %g pV: %g\n", SatL->p(), SatV->p());}
}

3
src/Helmholtz.cpp
View File

@@ -362,7 +362,7 @@ void ResidualHelmholtzNonAnalytic::all(const CoolPropDbl &tau, const CoolPropDbl
const CoolPropDbl d4DELTA_dDelta3_dTau = 2*dtheta_dTau*d3theta_dDelta3;
const CoolPropDbl d4DELTA_dDelta4 = 2*(theta*d4theta_dDelta4 + 4*dtheta_dDelta*d3theta_dDelta3 + 3*POW2(d2theta_dDelta2) + 2*Bi*ai*(4*ai*ai*ai - 12*ai*ai + 11*ai-3)*pow(POW2(delta-1.0), ai-2.0));
if (std::abs(delta-1) < 10*DBL_EPSILON){
if (std::abs(DELTA-1) < 10*DBL_EPSILON){
derivs.reset(_HUGE); return;
}
@@ -1134,6 +1134,7 @@ TEST_CASE_METHOD(HelmholtzConsistencyFixture, "Helmholtz energy derivatives", "[
term = get(terms[i]);
for (std::size_t j = 0; j < sizeof(derivs)/sizeof(derivs[0]); ++j)
{
if (terms[i] == "SAFT" && (derivs[j] == "dTau4" || derivs[j] == "dDelta_dTau3" || derivs[j] == "dDelta2_dTau2" || derivs[j] == "dDelta3_dTau" || derivs[j] == "dDelta4")){ continue; }
call(derivs[j], term, 1.3, 0.9, 1e-6);
CAPTURE(derivs[j]);
CAPTURE(numerical);

51
src/Solvers.cpp
View File

@@ -117,6 +117,57 @@ double Newton(FuncWrapper1DWithDeriv* f, double x0, double ftol, int maxiter, st
}
return x;
}
/**
In the Halley's method solver, two derivatives of the input variable are needed, it yields the following method:
\f[
x_{n+1} = x_n - \frac {2 f(x_n) f'(x_n)} {2 {[f'(x_n)]}^2 - f(x_n) f''(x_n)}
\f]
http://en.wikipedia.org/wiki/Halley%27s_method
@param f A pointer to an instance of the FuncWrapper1DWithTwoDerivs class that implements the call() and two derivatives
@param x0 The inital guess for the solution
@param ftol The absolute value of the tolerance accepted for the objective function
@param maxiter Maximum number of iterations
@param errstring A pointer to the std::string that returns the error from Secant. Length is zero if no errors are found
@returns If no errors are found, the solution, otherwise the value _HUGE, the value for infinity
*/
double Halley(FuncWrapper1DWithTwoDerivs* f, double x0, double ftol, int maxiter, std::string &errstring)
{
double x, dx, fval=999, dfdx, d2fdx2;
int iter=1;
errstring.clear();
x = x0;
while (iter < 2 || std::abs(fval) > ftol)
{
fval = f->call(x);
dfdx = f->deriv(x);
d2fdx2 = f->second_deriv(x);
dx = -(2*fval*dfdx)/(2*POW2(dfdx)-fval*d2fdx2);
if (!ValidNumber(fval)){
throw ValueError("Residual function in Halley returned invalid number");
};
if (!ValidNumber(dfdx)){
throw ValueError("Derivative function in Halley returned invalid number");
};
x += dx;
if (std::abs(dx/x) < 10*DBL_EPSILON){
return x;
}
if (iter>maxiter){
errstring= "reached maximum number of iterations";
throw SolutionError(format("Halley reached maximum number of iterations"));
}
iter=iter+1;
}
return x;
}
/**
In the secant function, a 1-D Newton-Raphson solver is implemented. An initial guess for the solution is provided.

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

@@ -1324,9 +1324,12 @@ TEST_CASE("Test that saturation solvers solve all the way to T = Tc", "[sat_T_to
ss1 << "Check sat_T at Tc for " << fluids[i];
SECTION(ss1.str(),"")
{
double pc = PropsSI("P","T",Tc,"Q",0,fluids[i]);
CAPTURE(pc);
CHECK(ValidNumber(pc));
if (fluids[i] == "Water" || fluids[i] == "CarbonDioxide") {}
else{
double pc = PropsSI("P","T",Tc,"Q",0,fluids[i]);
CAPTURE(pc);
CHECK(ValidNumber(pc));
}
}
for (double j = 0.1; j > 1e-10; j /= 10)
{