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

Preliminary work on two-phase derivatives and splines

Browse Source 
Added two-phase derivative, works!
Spline for drho/dh|p works
Spline for drho/dh|p NOT WORKING

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2015-01-01 22:58:05 -05:00
parent f28e2be4fa
commit 9e5c3d6035
5 changed files with 272 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -2545,7 +2545,26 @@ long double HelmholtzEOSMixtureBackend::calc_d3alpha0_dTau3(void)
const int nTau = 3, nDelta = 0;
return calc_alpha0_deriv_nocache(nTau, nDelta, mole_fractions, _tau, _delta, _reducing.T, _reducing.rhomolar);
}
long double HelmholtzEOSMixtureBackend::calc_first_saturation_deriv(parameters Of1, parameters Wrt1, HelmholtzEOSMixtureBackend &SatL, HelmholtzEOSMixtureBackend &SatV)
{
// Derivative of temperature w.r.t. pressure ALONG the saturation curve
long double dTdP_sat = T()*(1/SatV.rhomolar()-1/SatL.rhomolar())/(SatV.hmolar()-SatL.hmolar());
// "Trivial" inputs
if (Of1 == iT && Wrt1 == iP){ return dTdP_sat;}
else if (Of1 == iP && Wrt1 == iT){ return 1/dTdP_sat;}
// Derivative taken with respect to T
else if (Wrt1 == iT){
return first_partial_deriv(Of1, iT, iP) + first_partial_deriv(Of1, iP, iT)/dTdP_sat;
}
// Derivative taken with respect to p
else if (Wrt1 == iP){
return first_partial_deriv(Of1, iP, iT) + first_partial_deriv(Of1, iT, iP)*dTdP_sat;
}
else{
throw ValueError(format("Not possible to take first saturation derivative with respect to %s", get_parameter_information(Wrt1,"short").c_str()));
}
}
long double HelmholtzEOSMixtureBackend::calc_first_saturation_deriv(parameters Of1, parameters Wrt1)
{
// Derivative of temperature w.r.t. pressure ALONG the saturation curve
@@ -2596,5 +2615,143 @@ long double HelmholtzEOSMixtureBackend::calc_second_saturation_deriv(parameters
}
}
long double HelmholtzEOSMixtureBackend::calc_second_two_phase_deriv(parameters Of, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2)
{
if (Of == iDmolar && Wrt1 == iHmolar && Constant1 == iP && Wrt2 == iP && Constant2 == iHmolar){
parameters h_key = iHmolar, rho_key = iDmolar, p_key = iP;
// taking the derivative of (drho/dv)*(dv/dh|p) with respect to p with h constant
long double dv_dh_constp = calc_first_two_phase_deriv(rho_key,h_key,p_key)/(-POW2(rhomolar()));
long double drhomolar_dp__consth = first_two_phase_deriv(rho_key, p_key, h_key);
// Calculate the derivative of dvdh|p with respect to p at constant h
long double dhL_dp_sat = SatL->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
long double dhV_dp_sat = SatV->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
long double drhoL_dp_sat = SatL->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
long double drhoV_dp_sat = SatV->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
long double numerator = 1/SatV->keyed_output(rho_key) - 1/SatL->keyed_output(rho_key);
long double denominator = SatV->keyed_output(h_key) - SatL->keyed_output(h_key);
long double dnumerator = -1/POW2(SatV->keyed_output(rho_key))*drhoV_dp_sat - 1/POW2(SatL->keyed_output(rho_key))*drhoL_dp_sat;
long double ddenominator = dhV_dp_sat - dhL_dp_sat;
long double d_dvdh_dp__consth = (denominator*dnumerator - numerator*ddenominator)/POW2(denominator);
return -POW2(rhomolar())*d_dvdh_dp__consth + dv_dh_constp*(-2*rhomolar())*drhomolar_dp__consth;
}
else{
throw ValueError();
}
}
long double HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant)
{
if (Of == iDmolar && Wrt == iHmolar && Constant == iP){
return -POW2(rhomolar())*(1/SatV->rhomolar() - 1/SatL->rhomolar())/(SatV->hmolar() - SatL->hmolar());
}
else if (Of == iDmass && Wrt == iHmass && Constant == iP){
return -POW2(rhomass())*(1/SatV->rhomass() - 1/SatL->rhomass())/(SatV->hmass() - SatL->hmass());
}
else if (Of == iDmolar && Wrt == iP && Constant == iHmolar){
// v = 1/rho; dvdrho = -rho^2; dvdrho = -1/rho^2
long double dvdrhoL = -1/POW2(SatL->rhomolar());
long double dvdrhoV = -1/POW2(SatV->rhomolar());
long double dvL_dp = dvdrhoL*SatL->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
long double dvV_dp = dvdrhoV*SatV->calc_first_saturation_deriv(iDmolar, iP, *SatL, *SatV);
long double dhL_dp = SatL->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
long double dhV_dp = SatV->calc_first_saturation_deriv(iHmolar, iP, *SatL, *SatV);
long double dxdp_h = (Q()*dhV_dp + (1 - Q())*dhL_dp)/(SatL->hmolar() - SatV->hmolar());
long double dvdp_h = dvL_dp + dxdp_h*(1/SatV->rhomolar() - 1/SatL->rhomolar()) + Q()*(dvV_dp - dvL_dp);
return -POW2(rhomolar())*dvdp_h;
}
else if (Of == iDmass && Wrt == iP && Constant == iHmass){
// v = 1/rho; dvdrho = -rho^2; dvdrho = -1/rho^2
long double dvdrhoL = -1/POW2(SatL->rhomass());
long double dvdrhoV = -1/POW2(SatV->rhomass());
long double dvL_dp = dvdrhoL*SatL->calc_first_saturation_deriv(iDmass, iP, *SatL, *SatV);
long double dvV_dp = dvdrhoV*SatV->calc_first_saturation_deriv(iDmass, iP, *SatL, *SatV);
long double dhL_dp = SatL->calc_first_saturation_deriv(iHmass, iP, *SatL, *SatV);
long double dhV_dp = SatV->calc_first_saturation_deriv(iHmass, iP, *SatL, *SatV);
long double dxdp_h = (Q()*dhV_dp + (1 - Q())*dhL_dp)/(SatL->hmass() - SatV->hmass());
long double dvdp_h = dvL_dp + dxdp_h*(1/SatV->rhomass() - 1/SatL->rhomass()) + Q()*(dvV_dp - dvL_dp);
return -POW2(rhomass())*dvdp_h;
}
else{
throw ValueError("These inputs are not supported to calc_first_two_phase_deriv");
}
}
long double HelmholtzEOSMixtureBackend::calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, long double x_end)
{
shared_ptr<HelmholtzEOSMixtureBackend> Liq(new HelmholtzEOSMixtureBackend(this->get_components())),
End(new HelmholtzEOSMixtureBackend(this->get_components()));
Liq->specify_phase(iphase_liquid);
Liq->update(DmolarT_INPUTS, SatL->rhomolar(), SatL->T());
End->update(QT_INPUTS, x_end, SatL->T());
bool mass_based_inputs = ((Of == iDmass) && ((Wrt == iHmass && Constant == iP) || (Wrt == iP && Constant == iHmass) || (Wrt == iDmass && Constant == iDmass)));
bool mole_based_inputs = ((Of == iDmolar) && ((Wrt == iHmolar && Constant == iP) || (Wrt == iP && Constant == iHmolar) || (Wrt == iDmolar && Constant == iDmolar)));
if (mass_based_inputs || mole_based_inputs)
{
parameters p_key, h_key, rho_key;
if (mass_based_inputs){
rho_key = iDmass; h_key = iHmass; p_key = iP;
}
else{
rho_key = iDmolar; h_key = iHmolar; p_key = iP;
}
long double Delta = Q()*(SatV->keyed_output(h_key) - SatL->keyed_output(h_key));
long double Delta_end = End->keyed_output(h_key) - SatL->keyed_output(h_key);
// At the end of the zone to which spline is applied
long double drho_dh_end = End->calc_first_two_phase_deriv(rho_key, h_key, p_key);
long double rho_end = End->keyed_output(rho_key);
// Faking single-phase
long double rho_liq = Liq->keyed_output(rho_key);
long double drho_dh_liq = Liq->first_partial_deriv(rho_key, h_key, p_key);
// Spline coordinates
long double a = 1/POW3(Delta_end) * (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq + drho_dh_end));
long double b = 3/POW2(Delta_end) * (-rho_liq + rho_end) - 1/Delta_end * (drho_dh_end + 2 * drho_dh_liq);
long double c = drho_dh_liq;
long double d = rho_liq;
// Either the spline value or drho/dh|p can be directly evaluated now
long double rho_spline = a*POW3(Delta) + b*POW2(Delta) + c*Delta + d;
long double d_rho_spline_dh__constp = 3*a*POW2(Delta) + 2*b*Delta + c;
if ((Wrt == iDmass || Wrt == iDmolar) && (Constant == iDmass || Constant == iDmolar)){
return rho_spline;
}
if ((Wrt == iHmass || Wrt == iHmolar) && Constant == iP){
return d_rho_spline_dh__constp;
}
// Its drho/dp|h
// ... calculate some more things
// Faking single-phase
long double drhodhdp_liq = Liq->second_partial_deriv(rho_key, h_key, p_key, p_key, h_key);
long double drho_dp_liq = Liq->first_partial_deriv(rho_key, p_key, h_key);
// Derivatives *along* the saturation curve using the special internal method
long double dhL_dp_sat = SatL->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
long double dhV_dp_sat = SatV->calc_first_saturation_deriv(h_key, p_key, *SatL, *SatV);
long double drhoL_dp_sat = SatL->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
// Derivatives at the end point
long double drho_dp_end = End->calc_first_two_phase_deriv(rho_key, p_key, h_key);
long double drhodhdp_end = End->calc_second_two_phase_deriv(rho_key, h_key, p_key, p_key, h_key);
long double dDelta_end_dp = x_end * (dhL_dp_sat - dhV_dp_sat);
long double dDelta_dp = -dhL_dp_sat;
// First pressure derivative at constant h of the coefficients
long double da_dp = (-6/POW4(Delta_end)) * dDelta_end_dp * (rho_liq - rho_end) + (2/POW3(Delta_end)) * (drho_dp_liq - drho_dp_end) + (1/POW2(Delta_end)) * (drhodhdp_liq + drhodhdp_end) - (2/POW3(Delta_end)) * (drho_dh_liq + drho_dh_end)*dDelta_end_dp;
long double db_dp = (-6/POW3(Delta_end)) * dDelta_end_dp * (-rho_liq + rho_end) + (3/POW2(Delta_end)) * (-drho_dp_liq + drho_dp_end) + (1/POW2(Delta_end)) * dDelta_end_dp * (drho_dh_end + 2 * drho_dh_liq) - (1/Delta_end) * (drhodhdp_end + 2*drhodhdp_liq);
long double dc_dp = drhodhdp_liq;
long double dd_dp = drho_dp_liq;
long double d_rho_spline_dp__consth = (3*a*POW2(Delta) + 2*b*Delta + c)*dDelta_dp + POW3(Delta)*da_dp + POW2(Delta)*db_dp + Delta*dc_dp + dd_dp;
return d_rho_spline_dp__consth;
}
}
} /* namespace CoolProp */

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

@@ -69,7 +69,11 @@ public:
long double calc_ODP();
long double calc_first_saturation_deriv(parameters Of1, parameters Wrt1);
long double calc_first_saturation_deriv(parameters Of1, parameters Wrt1, HelmholtzEOSMixtureBackend &SatL, HelmholtzEOSMixtureBackend &SatV);
long double calc_second_saturation_deriv(parameters Of1, parameters Wrt1, parameters Of2, parameters Wrt2);
long double calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant);
long double calc_second_two_phase_deriv(parameters Of, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2);
long double calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, long double x_end);
/// Calculate the phase once the state is fully calculated but you aren't sure if it is liquid or gas or ...
void recalculate_singlephase_phase();

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

@@ -923,7 +923,7 @@ 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;
long double rhoL = _HUGE, rhoV = _HUGE, error = 999, DeltavL, DeltavV, pL, pV, p, rhoL0, rhoV0;
long double rhoL = _HUGE, rhoV = _HUGE, error = 999, DeltavL, DeltavV, pL, pV, p;
int iter=0, small_step_count = 0;
try
@@ -947,8 +947,6 @@ void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend &HE
rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T);
rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T);
p = HEOS.get_components()[0]->ancillaries.pV.evaluate(T);
rhoL0 = rhoL;
rhoV0 = rhoV;
CoolProp::SimpleState &crit = HEOS.get_components()[0]->crit;
CoolProp::SimpleState &tripleL = HEOS.get_components()[0]->triple_liquid;
@@ -1021,7 +1019,6 @@ void SaturationSolvers::saturation_T_pure_Maxwell(HelmholtzEOSMixtureBackend &HE
do{
pL = SatL->p(); pV = SatV->p();
long double vL = 1/SatL->rhomolar(), vV = 1/SatV->rhomolar();
long double gL = SatL->gibbsmolar(), gV = SatV->gibbsmolar();
// Get alpha, the pressure derivative with volume at constant T
// Given by (dp/drho|T)*drhodv
long double alphaL = SatL->first_partial_deriv(iP, iDmolar, iT)*(-POW2(SatL->rhomolar()));