mirror of
https://github.com/CoolProp/CoolProp.git
synced 2026-04-23 03:00:17 -04:00
Merge branch 'master' of https://github.com/coolprop/coolprop
This commit is contained in:
Ian Bell
@@ -290,6 +290,9 @@ protected:
|
||||
|
||||
virtual long double calc_first_saturation_deriv(parameters Of1, parameters Wrt1){throw NotImplementedError("calc_first_saturation_deriv is not implemented for this backend");};
|
||||
virtual long double calc_second_saturation_deriv(parameters Of1, parameters Wrt1, parameters Of2, parameters Wrt2){throw NotImplementedError("calc_second_saturation_deriv is not implemented for this backend");};
|
||||
virtual long double calc_first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant){throw NotImplementedError("calc_first_two_phase_deriv is not implemented for this backend");};
|
||||
virtual long double calc_second_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant, parameters Wrt2, parameters Constant2){throw NotImplementedError("calc_second_two_phase_deriv is not implemented for this backend");};
|
||||
virtual long double calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, long double x_end){throw NotImplementedError("calc_first_two_phase_deriv_splined is not implemented for this backend");};
|
||||
|
||||
virtual long double calc_saturated_liquid_keyed_output(parameters key){throw NotImplementedError("calc_saturated_liquid_keyed_output is not implemented for this backend");};
|
||||
virtual long double calc_saturated_vapor_keyed_output(parameters key){throw NotImplementedError("calc_saturated_vapor_keyed_output is not implemented for this backend");};
|
||||
@@ -571,6 +574,38 @@ public:
|
||||
* */
|
||||
long double second_saturation_deriv(parameters Of1, parameters Wrt1, parameters Of2, parameters Wrt2){return calc_second_saturation_deriv(Of1,Wrt1,Of2,Wrt2);};
|
||||
|
||||
/**
|
||||
* @brief Calculate the first "two-phase" derivative as described by Thorade and Sadaat, EAS, 2013
|
||||
*
|
||||
* Implementing the algorithms and ideas of:
|
||||
* Matthis Thorade, Ali Saadat, "Partial derivatives of thermodynamic state properties for dynamic simulation",
|
||||
* Environmental Earth Sciences, December 2013, Volume 70, Issue 8, pp 3497-3503
|
||||
*
|
||||
* Spline evaluation is as described in:
|
||||
* S Quoilin, I Bell, A Desideri, P Dewallef, V Lemort,
|
||||
* "Methods to increase the robustness of finite-volume flow models in thermodynamic systems",
|
||||
* Energies 7 (3), 1621-1640
|
||||
*
|
||||
* \note Not all derivatives are supported!
|
||||
*
|
||||
* @param Of The parameter to be derived
|
||||
* @param Wrt The parameter that the derivative is taken with respect to
|
||||
* @param Constant The parameter that is held constant
|
||||
* @param type_flag A flag describing how the derivative should be calculated, either normal or using splines
|
||||
* @return
|
||||
*/
|
||||
double first_two_phase_deriv(parameters Of, parameters Wrt, parameters Constant){
|
||||
return calc_first_two_phase_deriv(Of, Wrt, Constant);
|
||||
};
|
||||
|
||||
double second_two_phase_deriv(parameters Of, parameters Wrt1, parameters Constant1, parameters Wrt2, parameters Constant2){
|
||||
return calc_second_two_phase_deriv(Of, Wrt1, Constant1, Wrt2, Constant2);
|
||||
};
|
||||
|
||||
double first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, double x_end){
|
||||
return calc_first_two_phase_deriv_splined(Of, Wrt, Constant, x_end);
|
||||
};
|
||||
|
||||
// ----------------------------------------
|
||||
// Phase envelope for mixtures
|
||||
// ----------------------------------------
|
||||
|
||||
@@ -2560,7 +2560,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
|
||||
@@ -2611,5 +2630,183 @@ 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) || (Wrt2 == iHmolar && Constant2 == iP && Wrt1 == iP && Constant1 == 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 if (Of == iDmass && ((Wrt1 == iHmass && Constant1 == iP && Wrt2 == iP && Constant2 == iHmass) || (Wrt2 == iHmass && Constant2 == iP && Wrt1 == iP && Constant1 == iHmass))){
|
||||
parameters h_key = iHmass, rho_key = iDmass, p_key = iP;
|
||||
long double rho = keyed_output(rho_key);
|
||||
// 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(rho));
|
||||
long double drho_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(rho)*d_dvdh_dp__consth + dv_dh_constp*(-2*rho)*drho_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__constp = Liq->first_partial_deriv(rho_key, h_key, p_key);
|
||||
|
||||
// Spline coordinates a, b, c, d
|
||||
long double Abracket = (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end));
|
||||
long double a = 1/POW3(Delta_end) * Abracket;
|
||||
long double b = 3/POW2(Delta_end) * (-rho_liq + rho_end) - 1/Delta_end * (drho_dh_end + 2 * drho_dh_liq__constp);
|
||||
long double c = drho_dh_liq__constp;
|
||||
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;
|
||||
}
|
||||
|
||||
// It's drho/dp|h
|
||||
// ... calculate some more things
|
||||
|
||||
// At the saturated state
|
||||
long double rhoL = SatL->keyed_output(rho_key);
|
||||
long double rhoV = SatV->keyed_output(rho_key);
|
||||
long double hL = SatL->keyed_output(h_key);
|
||||
long double hV = SatV->keyed_output(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);
|
||||
long double drhoV_dp_sat = SatV->calc_first_saturation_deriv(rho_key, p_key, *SatL, *SatV);
|
||||
|
||||
long double drho_dp_end = POW2(End->keyed_output(rho_key))*(x_end/POW2(rhoV)*drhoV_dp_sat + (1-x_end)/POW2(rhoL)*drhoL_dp_sat);
|
||||
|
||||
// Faking single-phase
|
||||
long double drho_dp__consth_liq = Liq->first_partial_deriv(rho_key, p_key, h_key);
|
||||
long double d2rhodhdp_liq = Liq->second_partial_deriv(rho_key, h_key, p_key, p_key, h_key); // ?
|
||||
|
||||
// Derivatives at the end point
|
||||
long double drho_dp__consth_end = End->calc_first_two_phase_deriv(rho_key, p_key, h_key);
|
||||
long double d2rhodhdp_end = End->calc_second_two_phase_deriv(rho_key, h_key, p_key, p_key, h_key);
|
||||
|
||||
// Reminder:
|
||||
// Delta = Q()*(hV-hL) = h-hL
|
||||
// Delta_end = x_end*(hV-hL);
|
||||
long double d_Delta_dp__consth = -dhL_dp_sat;
|
||||
long double d_Delta_end_dp__consth = x_end*(dhV_dp_sat - dhL_dp_sat);
|
||||
|
||||
// First pressure derivative at constant h of the coefficients a,b,c,d
|
||||
// long double Abracket = (2*rho_liq - 2*rho_end + Delta_end * (drho_dh_liq__constp + drho_dh_end));
|
||||
long double d_Abracket_dp_consth = (2*drhoL_dp_sat - 2*drho_dp__consth_end + Delta_end*(d2rhodhdp_liq + d2rhodhdp_end) + d_Delta_dp__consth*(drho_dh_liq__constp + drho_dh_end));
|
||||
long double da_dp = 1/POW3(Delta_end)*d_Abracket_dp_consth + Abracket*(-3/POW4(Delta_end)*d_Delta_end_dp__consth);
|
||||
long double db_dp = - 6/POW3(Delta_end)*d_Delta_end_dp__consth*(rho_end - rho_liq)
|
||||
+ (3/POW2(Delta_end))*(drho_dp__consth_end - drhoL_dp_sat)
|
||||
+ (1/POW2(Delta_end)*d_Delta_end_dp__consth) * (drho_dh_end + 2*drho_dh_liq__constp)
|
||||
- (1/Delta_end) * (d2rhodhdp_end + 2*d2rhodhdp_liq);
|
||||
long double dc_dp = d2rhodhdp_liq;
|
||||
long double dd_dp = drhoL_dp_sat;
|
||||
|
||||
long double d_rho_spline_dp__consth = (3*a*POW2(Delta) + 2*b*Delta + c)*d_Delta_dp__consth + POW3(Delta)*da_dp + POW2(Delta)*db_dp + Delta*dc_dp + dd_dp;
|
||||
|
||||
return d_rho_spline_dp__consth;
|
||||
}
|
||||
else{
|
||||
throw ValueError("inputs to calc_first_two_phase_deriv_splined are currently invalid");
|
||||
}
|
||||
}
|
||||
|
||||
} /* namespace CoolProp */
|
||||
|
||||
|
||||
@@ -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();
|
||||
|
||||
@@ -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()));
|
||||
|
||||
@@ -3,6 +3,7 @@
|
||||
#include "AbstractState.h"
|
||||
#include "DataStructures.h"
|
||||
#include "../Backends/Helmholtz/HelmholtzEOSMixtureBackend.h"
|
||||
#include "../Backends/Helmholtz/HelmholtzEOSBackend.h"
|
||||
// ############################################
|
||||
// TESTS
|
||||
// ############################################
|
||||
@@ -1511,6 +1512,122 @@ TEST_CASE("Check the second saturation derivatives", "[second_saturation_partial
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("Check the first two-phase derivative", "[first_two_phase_deriv]")
|
||||
{
|
||||
const int number_of_pairs = 4;
|
||||
struct pair {parameters p1, p2, p3;};
|
||||
pair pairs[number_of_pairs] = {{iDmass, iP, iHmass}, {iDmolar, iP, iHmolar},
|
||||
{iDmolar, iHmolar, iP}, {iDmass, iHmass, iP}};
|
||||
shared_ptr<CoolProp::HelmholtzEOSBackend> AS(new CoolProp::HelmholtzEOSBackend("n-Propane"));
|
||||
for (std::size_t i = 0; i < number_of_pairs; ++i)
|
||||
{
|
||||
// See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU
|
||||
std::ostringstream ss1;
|
||||
ss1 << "for (" << get_parameter_information(pairs[i].p1,"short") << ", " << get_parameter_information(pairs[i].p2,"short") << ", " << get_parameter_information(pairs[i].p3,"short") << ")";
|
||||
SECTION(ss1.str(),"")
|
||||
{
|
||||
AS->update(QT_INPUTS, 0.3, 300);
|
||||
long double numerical;
|
||||
long double analytical = AS->first_two_phase_deriv(pairs[i].p1, pairs[i].p2, pairs[i].p3);
|
||||
CAPTURE(analytical);
|
||||
|
||||
long double out1, out2;
|
||||
long double v2base, v3base;
|
||||
v2base = AS->keyed_output(pairs[i].p2);
|
||||
v3base = AS->keyed_output(pairs[i].p3);
|
||||
long double v2plus = v2base*1.001;
|
||||
long double v2minus = v2base*0.999;
|
||||
CoolProp::input_pairs input_pair1 = generate_update_pair(pairs[i].p2, v2plus, pairs[i].p3, v3base, out1, out2);
|
||||
AS->update(input_pair1, out1, out2);
|
||||
long double v1 = AS->keyed_output(pairs[i].p1);
|
||||
CoolProp::input_pairs input_pair2 = generate_update_pair(pairs[i].p2, v2minus, pairs[i].p3, v3base, out1, out2);
|
||||
AS->update(input_pair2, out1, out2);
|
||||
long double v2 = AS->keyed_output(pairs[i].p1);
|
||||
|
||||
numerical = (v1 - v2)/(v2plus - v2minus);
|
||||
CAPTURE(numerical);
|
||||
CHECK(std::abs(numerical/analytical-1) < 1e-4);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("Check the second two-phase derivative", "[second_two_phase_deriv]")
|
||||
{
|
||||
SECTION("d2rhodhdp",""){
|
||||
shared_ptr<CoolProp::HelmholtzEOSBackend> AS(new CoolProp::HelmholtzEOSBackend("n-Propane"));
|
||||
AS->update(QT_INPUTS, 0.3, 300);
|
||||
long double analytical = AS->second_two_phase_deriv(iDmolar, iHmolar, iP, iP, iHmolar);
|
||||
CAPTURE(analytical);
|
||||
long double pplus = AS->p()*1.001, pminus = AS->p()*0.999, h = AS->hmolar();
|
||||
AS->update(HmolarP_INPUTS, h, pplus);
|
||||
long double v1 = AS->first_two_phase_deriv(iDmolar, iHmolar, iP);
|
||||
AS->update(HmolarP_INPUTS, h, pminus);
|
||||
long double v2 = AS->first_two_phase_deriv(iDmolar, iHmolar, iP);
|
||||
long double numerical = (v1 - v2)/(pplus - pminus);
|
||||
CAPTURE(numerical);
|
||||
CHECK(std::abs(numerical/analytical-1) < 1e-6);
|
||||
}
|
||||
SECTION("d2rhodhdp using mass",""){
|
||||
shared_ptr<CoolProp::HelmholtzEOSBackend> AS(new CoolProp::HelmholtzEOSBackend("n-Propane"));
|
||||
AS->update(QT_INPUTS, 0.3, 300);
|
||||
long double analytical = AS->second_two_phase_deriv(iDmass, iHmass, iP, iP, iHmass);
|
||||
CAPTURE(analytical);
|
||||
long double pplus = AS->p()*1.001, pminus = AS->p()*0.999, h = AS->hmass();
|
||||
AS->update(HmassP_INPUTS, h, pplus);
|
||||
long double v1 = AS->first_two_phase_deriv(iDmass, iHmass, iP);
|
||||
AS->update(HmassP_INPUTS, h, pminus);
|
||||
long double v2 = AS->first_two_phase_deriv(iDmass, iHmass, iP);
|
||||
long double numerical = (v1 - v2)/(pplus - pminus);
|
||||
CAPTURE(numerical);
|
||||
CHECK(std::abs(numerical/analytical-1) < 1e-6);
|
||||
}
|
||||
}
|
||||
|
||||
TEST_CASE("Check the first two-phase derivative using splines", "[first_two_phase_deriv_splined]")
|
||||
{
|
||||
const int number_of_pairs = 4;
|
||||
struct pair {parameters p1, p2, p3;};
|
||||
pair pairs[number_of_pairs] = {
|
||||
{iDmass, iP, iHmass},
|
||||
{iDmolar, iP, iHmolar},
|
||||
{iDmolar, iHmolar, iP},
|
||||
{iDmass, iHmass, iP}
|
||||
};
|
||||
shared_ptr<CoolProp::HelmholtzEOSBackend> AS(new CoolProp::HelmholtzEOSBackend("n-Propane"));
|
||||
for (std::size_t i = 0; i < number_of_pairs; ++i)
|
||||
{
|
||||
// See https://groups.google.com/forum/?fromgroups#!topic/catch-forum/mRBKqtTrITU
|
||||
std::ostringstream ss1;
|
||||
ss1 << "for (" << get_parameter_information(pairs[i].p1,"short") << ", " << get_parameter_information(pairs[i].p2,"short") << ", " << get_parameter_information(pairs[i].p3,"short") << ")";
|
||||
SECTION(ss1.str(),"")
|
||||
{
|
||||
AS->update(QT_INPUTS, 0.2, 300);
|
||||
long double numerical;
|
||||
long double analytical = AS->first_two_phase_deriv_splined(pairs[i].p1, pairs[i].p2, pairs[i].p3, 0.3);
|
||||
CAPTURE(analytical);
|
||||
|
||||
long double out1, out2;
|
||||
long double v2base, v3base;
|
||||
v2base = AS->keyed_output(pairs[i].p2);
|
||||
v3base = AS->keyed_output(pairs[i].p3);
|
||||
long double v2plus = v2base*1.001;
|
||||
long double v2minus = v2base*0.999;
|
||||
|
||||
CoolProp::input_pairs input_pair1 = generate_update_pair(pairs[i].p2, v2plus, pairs[i].p3, v3base, out1, out2);
|
||||
AS->update(input_pair1, out1, out2);
|
||||
long double v1 = AS->first_two_phase_deriv_splined(pairs[i].p1, pairs[i].p1, pairs[i].p1, 0.3);
|
||||
|
||||
CoolProp::input_pairs input_pair2 = generate_update_pair(pairs[i].p2, v2minus, pairs[i].p3, v3base, out1, out2);
|
||||
AS->update(input_pair2, out1, out2);
|
||||
long double v2 = AS->first_two_phase_deriv_splined(pairs[i].p1, pairs[i].p1, pairs[i].p1, 0.3);
|
||||
|
||||
numerical = (v1 - v2)/(v2plus - v2minus);
|
||||
CAPTURE(numerical);
|
||||
CHECK(std::abs(numerical/analytical-1) < 1e-10);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
TEST_CASE("Test that HS solver works for a few fluids", "[HS_solver]")
|
||||
{
|
||||
|
||||
@@ -1,16 +1,10 @@
|
||||
include lualib.mk
|
||||
|
||||
COOLPROP_CFLAGS ?= $(shell $(PKGCONFIG) --cflags CoolProp)
|
||||
COOLPROP_LDFLAGS ?= $(shell $(PKGCONFIG) --libs-only-L CoolProp)
|
||||
COOLPROP_LDLIBS ?= $(or $(shell $(PKGCONFIG) --libs-only-l CoolProp), -lCoolProp)
|
||||
COOLPROP_INCDIR ?= $(shell $(PKGCONFIG) --variable=includedir CoolProp)
|
||||
COOLPROP_LIBDIR ?= $(shell $(PKGCONFIG) --variable=libdir CoolProp)
|
||||
COOLPROP_HEADER ?= $(or $(COOLPROP_INCDIR)
|
||||
|
||||
CFLAGS ?= -g -O2 -Wall -Wextra -Wswitch-enum -Wwrite-strings -Wshadow
|
||||
XCFLAGS += -std=c99 -pedantic-errors -fPIC
|
||||
XCFLAGS += $(LUA_CFLAGS) $(COOLPROP_CFLAGS)
|
||||
XLDFLAGS += $(GUMBO_LDFLAGS) $(COOLPROP_LDLIBS)
|
||||
|
||||
all: coolprop/capi.so
|
||||
|
||||
|
||||
@@ -47,11 +47,19 @@ cdef class AbstractState:
|
||||
cpdef double molar_mass(self) except *
|
||||
cpdef double keyed_output(self, constants_header.parameters) except *
|
||||
|
||||
## ----------------------------------------
|
||||
## Derivatives
|
||||
## ----------------------------------------
|
||||
|
||||
cpdef long double first_partial_deriv(self, constants_header.parameters, constants_header.parameters, constants_header.parameters) except *
|
||||
cpdef long double second_partial_deriv(self, constants_header.parameters, constants_header.parameters, constants_header.parameters, constants_header.parameters, constants_header.parameters) except *
|
||||
cpdef long double first_saturation_deriv(self, constants_header.parameters, constants_header.parameters) except *
|
||||
cpdef long double second_saturation_deriv(self, constants_header.parameters, constants_header.parameters, constants_header.parameters, constants_header.parameters) except *
|
||||
|
||||
cpdef double first_two_phase_deriv(self, constants_header.parameters Of, constants_header.parameters Wrt, constants_header.parameters Constant) except *
|
||||
cpdef double second_two_phase_deriv(self, constants_header.parameters Of, constants_header.parameters Wrt1, constants_header.parameters Constant1, constants_header.parameters Wrt2, constants_header.parameters Constant2) except *
|
||||
cpdef double first_two_phase_deriv_splined(self ,constants_header.parameters Of, constants_header.parameters Wrt, constants_header.parameters Constant, double x_end) except *
|
||||
|
||||
cpdef double melting_line(self, int, int, double) except *
|
||||
cpdef bool has_melting_line(self) except *
|
||||
|
||||
|
||||
@@ -126,6 +126,15 @@ cdef class AbstractState:
|
||||
cpdef long double second_saturation_deriv(self, constants_header.parameters OF1 , constants_header.parameters WRT1, constants_header.parameters OF2, constants_header.parameters WRT2) except *:
|
||||
""" Get the second derivative along the saturation curve - wrapper of c++ function :cpapi:`CoolProp::AbstractState::second_saturation_deriv` """
|
||||
return self.thisptr.second_saturation_deriv(OF1, WRT1, OF2, WRT2)
|
||||
cpdef double first_two_phase_deriv(self, constants_header.parameters Of, constants_header.parameters Wrt, constants_header.parameters Constant) except *:
|
||||
""" Get the first two-phase derivative - wrapper of C++ function :cpapi:`CoolProp::AbstractState::first_two_phase_deriv` """
|
||||
return self.thisptr.first_two_phase_deriv(Of, Wrt, Constant)
|
||||
cpdef double second_two_phase_deriv(self, constants_header.parameters Of1, constants_header.parameters Wrt1, constants_header.parameters Constant1, constants_header.parameters Wrt2, constants_header.parameters Constant2) except *:
|
||||
""" Get the second two-phase derivative - wrapper of C++ function :cpapi:`CoolProp::AbstractState::second_two_phase_deriv` """
|
||||
return self.thisptr.second_two_phase_deriv(Of1, Wrt1, Constant1, Wrt2, Constant2)
|
||||
cpdef double first_two_phase_deriv_splined(self, constants_header.parameters Of, constants_header.parameters Wrt, constants_header.parameters Constant, double x_end) except *:
|
||||
""" Get the first two-phase derivative using splines - wrapper of C++ function :cpapi:`CoolProp::AbstractState::first_two_phase_deriv_splined` """
|
||||
return self.thisptr.first_two_phase_deriv_splined(Of, Wrt, Constant, x_end)
|
||||
|
||||
## ----------------------------------------
|
||||
## Melting Line
|
||||
|
||||
@@ -64,6 +64,10 @@ cdef extern from "AbstractState.h" namespace "CoolProp":
|
||||
long double first_saturation_deriv(constants_header.parameters, constants_header.parameters) except+ValueError
|
||||
long double second_saturation_deriv(constants_header.parameters, constants_header.parameters, constants_header.parameters, constants_header.parameters) except+ValueError
|
||||
|
||||
double first_two_phase_deriv(constants_header.parameters Of, constants_header.parameters Wrt, constants_header.parameters Constant) except+ValueError
|
||||
double second_two_phase_deriv(constants_header.parameters Of, constants_header.parameters Wrt1, constants_header.parameters Constant1, constants_header.parameters Wrt2, constants_header.parameters Constant2) except+ValueError
|
||||
double first_two_phase_deriv_splined(constants_header.parameters Of, constants_header.parameters Wrt, constants_header.parameters Constant, double x_end) except+ValueError
|
||||
|
||||
void set_mole_fractions(vector[double]) except+ValueError
|
||||
void set_mass_fractions(vector[double]) except+ValueError
|
||||
void set_volu_fractions(vector[double]) except+ValueError
|
||||
|
||||
Reference in New Issue
Block a user