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

First implementation of critical point solver (much slower than on windows?)

Browse Source 
``` C++
    std::vector<std::string> names; names.push_back("Methane"); names.push_back("H2S");
    shared_ptr<CoolProp::HelmholtzEOSMixtureBackend> mix(new CoolProp::HelmholtzEOSMixtureBackend(names));
    mix->specify_phase(iphase_gas);
    std::vector<CoolPropDbl> z(2,0);
    for (double x0 = 0.01; x0 < 0.5; x0 += 0.01)
    {
        z[0] = x0; z[1] = 1-x0;
        mix->set_mole_fractions(z);
        std::vector<CoolProp::SimpleState> critpts = mix->find_all_critical_points();
        int rr = 4;
    }
```
This commit is contained in:
Ian Bell
2015-08-03 22:04:44 -06:00
parent a0ebb5f96c
commit fed41d4518
2 changed files with 131 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -2945,5 +2945,135 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0)
_critical.rhomolar = x[1]*rhomolar_reducing();
}
class OneDimObjective : public FuncWrapper1DWithDeriv
{
public:
CoolProp::HelmholtzEOSMixtureBackend &HEOS;
double delta, R_tau, R_delta;
OneDimObjective(HelmholtzEOSMixtureBackend &HEOS, double delta0) : HEOS(HEOS), delta(delta0) {};
double call(double tau){
double rhomolar = HEOS.rhomolar_reducing()*delta, T = HEOS.T_reducing()/tau;
HEOS.update_DmolarT_direct(rhomolar, T);
double L1 = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT).determinant();
//std::cout << L1 << std::endl;
return L1;
}
double deriv(double tau){
std::size_t N = HEOS.get_mole_fractions().size();
Eigen::MatrixXd adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT), N),
dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau);
return (adjL*dLdTau).trace();
};
};
class L0CurveTracer : public FuncWrapper1DWithDeriv
{
public:
CoolProp::HelmholtzEOSMixtureBackend &HEOS;
double delta, tau, M1_last, R_tau, R_delta;
std::vector<CoolProp::SimpleState> critical_points;
int N_critical_points;
L0CurveTracer(HelmholtzEOSMixtureBackend &HEOS, double tau0, double delta0) : HEOS(HEOS), delta(delta0), tau(tau0), N_critical_points(0)
{
R_delta = 0.1; R_tau = 0.1;
};
/***
\brief Calculate the value of L1
@param theta The angle
@param tau The old value of tau
@param delta The old value of delta
@param tau_new The new value of tau
@param delta_new The new value of delta
*/
void get_tau_delta(const double theta, const double tau, const double delta, double &tau_new, double &delta_new){
tau_new = tau + R_tau*cos(theta);
delta_new = delta + R_delta*sin(theta);
};
/***
\brief Calculate the value of L1
@param theta The angle
*/
double call(double theta){
double tau_new, delta_new;
this->get_tau_delta(theta, tau, delta, tau_new, delta_new);
double rhomolar = HEOS.rhomolar_reducing()*delta_new, T = HEOS.T_reducing()/tau_new;
HEOS.update_DmolarT_direct(rhomolar, T);
double L1 = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT).determinant();
//std::cout << "call: " << theta << " " << L1 << std::endl;
return L1;
};
/***
\brief Calculate the first partial derivative of L1 with respect to the angle
@param theta The angle
*/
double deriv(double theta){
std::size_t N = HEOS.get_mole_fractions().size();
Eigen::MatrixXd adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT), N),
dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau),
dLdDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta);
double dL1_dtau = (adjL*dLdTau).trace(), dL1_ddelta = (adjL*dLdDelta).trace();
return -R_tau*sin(theta)*dL1_dtau + R_delta*cos(theta)*dL1_ddelta;
};
void trace(){
double theta;
static std::string errstr;
for (int i = 0; i < 100; ++i){
if (i == 0){
// In the first iteration, search all angles in the positive delta direction using a
// bounded solver
theta = Brent(this, 0, M_PI, DBL_EPSILON, 1e-10, 100, errstr);
}
else{
// In subsequent iterations, you already have an excellent guess for the direction to
// be searching in, use Newton's method to refine the solution
theta = Newton(this, theta, 1e-10, 100, errstr);
}
// Calculate the second criticality condition
double M1 = MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT).determinant();
double p_MPa = HEOS.p()/1e6;
// Calculate the new tau and delta at the new point
double tau_new, delta_new;
this->get_tau_delta(theta, tau, delta, tau_new, delta_new);
// Stop if bounds are exceeded
if (p_MPa > 500 || delta_new > 5){
break;
}
// If the sign of M1 and the previous value of M1 have different signs, it means that
// you have bracketed a critical point, run the full critical point solver to
// find the critical point and store it
if (i > 0 && M1*M1_last < 0){
double rhomolar = HEOS.rhomolar_reducing()*delta_new, T = HEOS.T_reducing()/tau_new;
HEOS.calc_critical_point(rhomolar, T);
CoolProp::SimpleState crit = HEOS.get_state("critical");
critical_points.push_back(crit);
N_critical_points++;
std::cout << HEOS.get_mole_fractions()[0] << " " << crit.rhomolar << " " << crit.T << " " << p_MPa << std::endl;
}
// Update the storage values
this->tau = tau_new;
this->delta = delta_new;
this->M1_last = M1;
}
};
};
std::vector<CoolProp::SimpleState> HelmholtzEOSMixtureBackend::find_all_critical_points()
{
std::string errstr;
OneDimObjective resid_L0(*this, 0.5);
double tau_L0 = Newton(resid_L0, 1.3, 1e-10, 100, errstr);
L0CurveTracer tracer(*this, tau_L0, 0.5);
tracer.trace();
return tracer.critical_points;
}
} /* namespace CoolProp */

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

@@ -83,6 +83,7 @@ public:
CoolPropDbl calc_first_two_phase_deriv_splined(parameters Of, parameters Wrt, parameters Constant, CoolPropDbl x_end);
void calc_critical_point(double rho0, double T0);
std::vector<CoolProp::SimpleState> find_all_critical_points();
/// 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();