github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-01-22 20:38:01 -05:00
Code Issues Packages Projects Releases Wiki Activity

Excess function now using cached values in Departure function matrix; closes #689

Browse Source
This commit is contained in:
Ian Bell
2015-06-05 00:34:57 -06:00
parent 9340a6e846
commit aa55db2d42
3 changed files with 134 additions and 144 deletions
Download Patch File Download Diff File Expand all files Collapse all files

153
src/Backends/Helmholtz/ExcessHEFunction.h
View File

@@ -20,19 +20,26 @@ class DepartureFunction
public:
DepartureFunction(){};
virtual ~DepartureFunction(){};
ResidualHelmholtzGeneralizedExponential phi;
HelmholtzDerivatives derivs;
/// The excess Helmholtz energy of the binary pair
/// Pure-virtual function (must be implemented in derived class
virtual double alphar(double tau, double delta) = 0;
virtual double dalphar_dDelta(double tau, double delta) = 0;
virtual double d2alphar_dDelta2(double tau, double delta) = 0;
virtual double d2alphar_dDelta_dTau(double tau, double delta) = 0;
virtual double dalphar_dTau(double tau, double delta) = 0;
virtual double d2alphar_dTau2(double tau, double delta) = 0;
virtual double d3alphar_dTau3(double tau, double delta) = 0;
virtual double d3alphar_dDelta_dTau2(double tau, double delta) = 0;
virtual double d3alphar_dDelta2_dTau(double tau, double delta) = 0;
virtual double d3alphar_dDelta3(double tau, double delta) = 0;
void update(double tau, double delta){
derivs.reset(0.0);
phi.all(tau, delta, derivs);
};
double alphar(){ return derivs.alphar;};
double dalphar_dDelta(){ return derivs.dalphar_ddelta;};
double dalphar_dTau(){ return derivs.dalphar_dtau;};
double d2alphar_dDelta2(){return derivs.d2alphar_ddelta2;};
double d2alphar_dDelta_dTau(){return derivs.d2alphar_ddelta_dtau;};
double d2alphar_dTau2(){return derivs.d2alphar_dtau2;};
double d3alphar_dTau3(){ return derivs.d3alphar_dtau3; };
double d3alphar_dDelta_dTau2(){ return derivs.d3alphar_ddelta_dtau2; };
double d3alphar_dDelta2_dTau(){ return derivs.d3alphar_ddelta2_dtau; };
double d3alphar_dDelta3(){ return derivs.d3alphar_ddelta3; };
};
/** \brief The departure function used by the GERG-2008 formulation
@@ -44,10 +51,7 @@ public:
* It is symmetric so \f$\alphar^r_{ij} = \alphar^r_{ji}\f$
*/
class GERG2008DepartureFunction : public DepartureFunction
{
protected:
bool using_gaussian;
ResidualHelmholtzGeneralizedExponential phi;
{
public:
GERG2008DepartureFunction(){};
GERG2008DepartureFunction(const std::vector<double> &n,const std::vector<double> &d,const std::vector<double> &t,
@@ -64,11 +68,9 @@ public:
}
if (n.size() == Npower)
{
using_gaussian = false;
}
else
{
using_gaussian = true;
std::vector<CoolPropDbl> _n(n.begin()+Npower, n.end());
std::vector<CoolPropDbl> _d(d.begin()+Npower, d.end());
std::vector<CoolPropDbl> _t(t.begin()+Npower, t.end());
@@ -80,17 +82,6 @@ public:
}
};
~GERG2008DepartureFunction(){};
double alphar(double tau, double delta){return phi.base(tau, delta);};
double dalphar_dDelta(double tau, double delta){return phi.dDelta(tau, delta);};
double d2alphar_dDelta_dTau(double tau, double delta){return phi.dDelta_dTau(tau, delta);};
double dalphar_dTau(double tau, double delta){return phi.dTau(tau, delta);};
double d2alphar_dDelta2(double tau, double delta){return phi.dDelta2(tau, delta);};
double d2alphar_dTau2(double tau, double delta){return phi.dTau2(tau, delta);};
double d3alphar_dTau3(double tau, double delta){ return phi.dTau3(tau, delta); };
double d3alphar_dDelta_dTau2(double tau, double delta){ return phi.dDelta_dTau2(tau, delta); };
double d3alphar_dDelta2_dTau(double tau, double delta){ return phi.dDelta2_dTau(tau, delta); };
double d3alphar_dDelta3(double tau, double delta){ return phi.dDelta3(tau, delta); };
};
/** \brief A polynomial/exponential departure function
@@ -103,8 +94,6 @@ public:
*/
class ExponentialDepartureFunction : public DepartureFunction
{
protected:
ResidualHelmholtzGeneralizedExponential phi;
public:
ExponentialDepartureFunction(){};
ExponentialDepartureFunction(const std::vector<double> &n, const std::vector<double> &d,
@@ -117,19 +106,6 @@ public:
phi.add_Power(_n, _d, _t, _l);
};
~ExponentialDepartureFunction(){};
double alphar(double tau, double delta){return phi.base(tau, delta);};
double dalphar_dDelta(double tau, double delta){return phi.dDelta(tau, delta);};
double d2alphar_dDelta_dTau(double tau, double delta){return phi.dDelta_dTau(tau, delta);};
double dalphar_dTau(double tau, double delta){return phi.dTau(tau, delta);};
double d2alphar_dDelta2(double tau, double delta){return phi.dDelta2(tau, delta);};
double d2alphar_dTau2(double tau, double delta){return phi.dTau2(tau, delta);};
double d3alphar_dTau3(double tau, double delta){ return phi.dTau3(tau, delta); };
double d3alphar_dDelta_dTau2(double tau, double delta){ return phi.dDelta_dTau2(tau, delta); };
double d3alphar_dDelta2_dTau(double tau, double delta){ return phi.dDelta2_dTau(tau, delta); };
double d3alphar_dDelta3(double tau, double delta){ return phi.dDelta3(tau, delta); };
};
typedef shared_ptr<DepartureFunction> DepartureFunctionPointer;
@@ -152,233 +128,244 @@ public:
DepartureFunctionMatrix[i].resize(N);
}
};
/// Update the internal cached derivatives in each departure function
void update(double tau, double delta){
for (std::size_t i = 0; i < N; i++){
for (std::size_t j = i + 1; j < N; j++){
DepartureFunctionMatrix[i][j]->update(tau, delta);
}
for (std::size_t j = 0; j < i; j++){
DepartureFunctionMatrix[i][j]->update(tau, delta);
}
}
}
double alphar(double tau, double delta, const std::vector<CoolPropDbl> &x)
double alphar(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->alphar(tau,delta);
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->alphar();
}
}
return summer;
}
double dalphar_dDelta(double tau, double delta, const std::vector<CoolPropDbl> &x)
double dalphar_dDelta(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dDelta(tau,delta);
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dDelta();
}
}
return summer;
}
double d2alphar_dDelta2(double tau, double delta, const std::vector<CoolPropDbl> &x)
double d2alphar_dDelta2(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta2(tau,delta);
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta2();
}
}
return summer;
};
double d2alphar_dDelta_dTau(double tau, double delta, const std::vector<CoolPropDbl> &x)
double d2alphar_dDelta_dTau(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta_dTau(tau,delta);
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta_dTau();
}
}
return summer;
}
double dalphar_dTau(double tau, double delta, const std::vector<CoolPropDbl> &x)
double dalphar_dTau(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dTau(tau,delta);
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dTau();
}
}
return summer;
};
double d2alphar_dTau2(double tau, double delta, const std::vector<CoolPropDbl> &x)
double d2alphar_dTau2(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N-1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dTau2(tau,delta);
summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dTau2();
}
}
return summer;
};
double d3alphar_dTau3(double tau, double delta, const std::vector<CoolPropDbl> &x)
double d3alphar_dTau3(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N - 1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->d3alphar_dTau3(tau, delta);
summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->d3alphar_dTau3();
}
}
return summer;
};
double d3alphar_dDelta_dTau2(double tau, double delta, const std::vector<CoolPropDbl> &x)
double d3alphar_dDelta_dTau2(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N - 1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->d3alphar_dDelta_dTau2(tau, delta);
summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->d3alphar_dDelta_dTau2();
}
}
return summer;
};
double d3alphar_dDelta2_dTau(double tau, double delta, const std::vector<CoolPropDbl> &x)
double d3alphar_dDelta2_dTau(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N - 1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->d3alphar_dDelta2_dTau(tau, delta);
summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->d3alphar_dDelta2_dTau();
}
}
return summer;
};
double d3alphar_dDelta3(double tau, double delta, const std::vector<CoolPropDbl> &x)
double d3alphar_dDelta3(const std::vector<CoolPropDbl> &x)
{
double summer = 0;
for (std::size_t i = 0; i < N - 1; i++)
{
for (std::size_t j = i + 1; j < N; j++)
{
summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->d3alphar_dDelta3(tau, delta);
summer += x[i] * x[j] * F[i][j] * DepartureFunctionMatrix[i][j]->d3alphar_dDelta3();
}
}
return summer;
};
double dalphar_dxi(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i)
double dalphar_dxi(const std::vector<CoolPropDbl> &x, std::size_t i)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->alphar(tau,delta);
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->alphar();
}
}
return summer;
};
double d2alphardxidxj(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j)
double d2alphardxidxj(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j)
{
if (i != j)
{
return F[i][j]*DepartureFunctionMatrix[i][j]->alphar(tau,delta);
return F[i][j]*DepartureFunctionMatrix[i][j]->alphar();
}
else
{
return 0;
}
};
double d3alphar_dxi_dxj_dDelta(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j)
double d3alphar_dxi_dxj_dDelta(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j)
{
if (i != j)
{
return F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dDelta(tau, delta);
return F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dDelta();
}
else
{
return 0;
}
};
double d3alphar_dxi_dxj_dTau(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j)
double d3alphar_dxi_dxj_dTau(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j)
{
if (i != j)
{
return F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dTau(tau, delta);
return F[i][j]*DepartureFunctionMatrix[i][j]->dalphar_dTau();
}
else
{
return 0;
}
};
double d3alphardxidxjdxk(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k)
double d3alphardxidxjdxk(const std::vector<CoolPropDbl> &x, std::size_t i, std::size_t j, std::size_t k)
{
return 0;
};
double d2alphar_dxi_dTau(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i)
double d2alphar_dxi_dTau(const std::vector<CoolPropDbl> &x, std::size_t i)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dTau(tau,delta);
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dTau();
}
}
return summer;
};
double d2alphar_dxi_dDelta(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i)
double d2alphar_dxi_dDelta(const std::vector<CoolPropDbl> &x, std::size_t i)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dDelta(tau,delta);
summer += x[k]*F[i][k]*DepartureFunctionMatrix[i][k]->dalphar_dDelta();
}
}
return summer;
};
double d3alphar_dxi_dDelta2(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i)
double d3alphar_dxi_dDelta2(const std::vector<CoolPropDbl> &x, std::size_t i)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta2(tau, delta);
summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta2();
}
}
return summer;
};
double d3alphar_dxi_dTau2(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i)
double d3alphar_dxi_dTau2(const std::vector<CoolPropDbl> &x, std::size_t i)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dTau2(tau, delta);
summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dTau2();
}
}
return summer;
};
double d3alphar_dxi_dDelta_dTau(double tau, double delta, const std::vector<CoolPropDbl> &x, std::size_t i)
double d3alphar_dxi_dDelta_dTau(const std::vector<CoolPropDbl> &x, std::size_t i)
{
double summer = 0;
for (std::size_t k = 0; k < N; k++)
{
if (i != k)
{
summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau(tau, delta);
summer += x[k] * F[i][k] * DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau();
}
}
return summer;

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

@@ -943,7 +943,7 @@ void HelmholtzEOSMixtureBackend::pre_update(CoolProp::input_pairs &input_pair, C
void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double value1, double value2 )
{
if (get_debug_level() > 10){std::cout << format("%s (%d): update called with (%d: (%s), %g, %g)",__FILE__,__LINE__, input_pair, get_input_pair_short_desc(input_pair).c_str(), value1, value2) << std::endl;}
CoolPropDbl ld_value1 = value1, ld_value2 = value2;
pre_update(input_pair, ld_value1, ld_value2);
value1 = ld_value1; value2 = ld_value2;
@@ -2413,17 +2413,18 @@ void HelmholtzEOSMixtureBackend::calc_all_alphar_deriv_cache(const std::vector<C
summer_dDelta_dTau2 += xi*derivs.d3alphar_ddelta_dtau2;
summer_dTau3 += xi*derivs.d3alphar_dtau3;
}
_alphar = summer_base + Excess.alphar(tau, delta, mole_fractions);
_dalphar_dDelta = summer_dDelta + Excess.dalphar_dDelta(tau, delta, mole_fractions);
_dalphar_dTau = summer_dTau + Excess.dalphar_dTau(tau, delta, mole_fractions);
_d2alphar_dDelta2 = summer_dDelta2 + Excess.d2alphar_dDelta2(tau, delta, mole_fractions);
_d2alphar_dDelta_dTau = summer_dDelta_dTau + Excess.d2alphar_dDelta_dTau(tau, delta, mole_fractions);
_d2alphar_dTau2 = summer_dTau2 + Excess.d2alphar_dTau2(tau, delta, mole_fractions);
Excess.update(tau, delta);
_alphar = summer_base + Excess.alphar(mole_fractions);
_dalphar_dDelta = summer_dDelta + Excess.dalphar_dDelta(mole_fractions);
_dalphar_dTau = summer_dTau + Excess.dalphar_dTau(mole_fractions);
_d2alphar_dDelta2 = summer_dDelta2 + Excess.d2alphar_dDelta2(mole_fractions);
_d2alphar_dDelta_dTau = summer_dDelta_dTau + Excess.d2alphar_dDelta_dTau(mole_fractions);
_d2alphar_dTau2 = summer_dTau2 + Excess.d2alphar_dTau2(mole_fractions);
_d3alphar_dDelta3 = summer_dDelta3 + Excess.d3alphar_dDelta3(tau, delta, mole_fractions);
_d3alphar_dDelta2_dTau = summer_dDelta2_dTau + Excess.d3alphar_dDelta2_dTau(tau, delta, mole_fractions);
_d3alphar_dDelta_dTau2 = summer_dDelta_dTau2 + Excess.d3alphar_dDelta_dTau2(tau, delta, mole_fractions);
_d3alphar_dTau3 = summer_dTau3 + Excess.d3alphar_dTau3(tau, delta, mole_fractions);
_d3alphar_dDelta3 = summer_dDelta3 + Excess.d3alphar_dDelta3(mole_fractions);
_d3alphar_dDelta2_dTau = summer_dDelta2_dTau + Excess.d3alphar_dDelta2_dTau(mole_fractions);
_d3alphar_dDelta_dTau2 = summer_dDelta_dTau2 + Excess.d3alphar_dDelta_dTau2(mole_fractions);
_d3alphar_dTau3 = summer_dTau3 + Excess.d3alphar_dTau3(mole_fractions);
//_d4alphar_dDelta4 = derivs.d4alphar_ddelta4;
//_d4alphar_dDelta3_dTau = derivs.d4alphar_ddelta3_dtau;
//_d4alphar_dDelta2_dTau2 = derivs.d4alphar_ddelta2_dtau2;
@@ -2477,27 +2478,27 @@ CoolPropDbl HelmholtzEOSMixtureBackend::calc_alphar_deriv_nocache(const int nTau
CoolPropDbl summer = 0;
if (nTau == 0 && nDelta == 0){
for (unsigned int i = 0; i < N; ++i){ summer += mole_fractions[i]*components[i].EOS().baser(tau, delta); }
return summer + Excess.alphar(tau, delta, mole_fractions);
return summer + Excess.alphar(mole_fractions);
}
else if (nTau == 0 && nDelta == 1){
for (unsigned int i = 0; i < N; ++i){ summer += mole_fractions[i]*components[i].EOS().dalphar_dDelta(tau, delta); }
return summer + Excess.dalphar_dDelta(tau, delta, mole_fractions);
return summer + Excess.dalphar_dDelta(mole_fractions);
}
else if (nTau == 1 && nDelta == 0){
for (unsigned int i = 0; i < N; ++i){ summer += mole_fractions[i]*components[i].EOS().dalphar_dTau(tau, delta); }
return summer + Excess.dalphar_dTau(tau, delta, mole_fractions);
return summer + Excess.dalphar_dTau(mole_fractions);
}
else if (nTau == 0 && nDelta == 2){
for (unsigned int i = 0; i < N; ++i){ summer += mole_fractions[i]*components[i].EOS().d2alphar_dDelta2(tau, delta); }
return summer + Excess.d2alphar_dDelta2(tau, delta, mole_fractions);
return summer + Excess.d2alphar_dDelta2(mole_fractions);
}
else if (nTau == 1 && nDelta == 1){
for (unsigned int i = 0; i < N; ++i){ summer += mole_fractions[i]*components[i].EOS().d2alphar_dDelta_dTau(tau, delta); }
return summer + Excess.d2alphar_dDelta_dTau(tau, delta, mole_fractions);
return summer + Excess.d2alphar_dDelta_dTau(mole_fractions);
}
else if (nTau == 2 && nDelta == 0){
for (unsigned int i = 0; i < N; ++i){ summer += mole_fractions[i]*components[i].EOS().d2alphar_dTau2(tau, delta); }
return summer + Excess.d2alphar_dTau2(tau, delta, mole_fractions);
return summer + Excess.d2alphar_dTau2(mole_fractions);
}
/*else if (nTau == 0 && nDelta == 3){
for (unsigned int i = 0; i < N; ++i){ summer += mole_fractions[i]*components[i].EOS().d3alphar_dDelta3(tau, delta); }
@@ -3000,6 +3001,7 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0)
M1 = MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT).determinant();
std::vector<double> o(2);
o[0] = L1; o[1] = M1;
double p = HEOS.p();
return o;
};
std::vector<std::vector<double> > Jacobian(const std::vector<double> &x)
@@ -3016,13 +3018,13 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0)
J0(0,0) = (adjL*dLdTau).trace();
J0(0,1) = (adjL*dLdDelta).trace();
//std::cout << J0 << std::endl;
std::cout << J0 << std::endl;
std::vector<double> r0 = call(x);
// Build the Jacobian by column
for (std::size_t i = 0; i < N; ++i)
{
xp = x;
epsilon = 0.0001;
epsilon = 0.00001;
xp[i] += epsilon;
r = call(xp);
@@ -3036,6 +3038,7 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0)
};
};
Resid resid(*this);
CoolPropDbl Tr = T_reducing();
std::vector<double> x, tau_delta(2); tau_delta[0] = T_reducing()/T0; tau_delta[1] = rho0/rhomolar_reducing();
std::vector<double> td2 = tau_delta;
std::string errstr;

84
src/Backends/Helmholtz/MixtureDerivatives.cpp
View File

@@ -5,19 +5,19 @@ namespace CoolProp{
CoolPropDbl MixtureDerivatives::dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return HEOS.components[i].EOS().baser(HEOS._tau, HEOS._delta) + HEOS.Excess.dalphar_dxi(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
return HEOS.components[i].EOS().baser(HEOS.tau(), HEOS.delta()) + HEOS.Excess.dalphar_dxi(HEOS.mole_fractions, i);
}
else if(xN_flag == XN_DEPENDENT){
std::vector<CoolPropDbl> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i == N-1) return 0;
double dar_dxi = HEOS.components[i].EOS().baser(HEOS._tau, HEOS._delta) - HEOS.components[N-1].EOS().baser(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->alphar(HEOS._tau, HEOS._delta);
double dar_dxi = HEOS.components[i].EOS().baser(HEOS.tau(), HEOS.delta()) - HEOS.components[N-1].EOS().baser(HEOS.tau(), HEOS.delta());
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->alphar();
dar_dxi += (1-2*x[i])*FiNariN;
for (std::size_t k = 0; k < N-1; ++k){
if (i == k) continue;
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->alphar(HEOS._tau, HEOS._delta);
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->alphar(HEOS._tau, HEOS._delta);
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->alphar();
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->alphar();
dar_dxi += x[k]*(Fikarik - FiNariN - FkNarkN);
}
return dar_dxi;
@@ -29,19 +29,19 @@ CoolPropDbl MixtureDerivatives::dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, st
CoolPropDbl MixtureDerivatives::d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return HEOS.components[i].EOS().dalphar_dTau(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dTau(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
return HEOS.components[i].EOS().dalphar_dTau(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dTau(HEOS.mole_fractions, i);
}
else if(xN_flag == XN_DEPENDENT){
std::vector<CoolPropDbl> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return 0;
double d2ar_dxi_dTau = HEOS.components[i].EOS().dalphar_dTau(HEOS._tau, HEOS._delta) - HEOS.components[N-1].EOS().dalphar_dTau(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dTau(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dTau();
d2ar_dxi_dTau += (1-2*x[i])*FiNariN;
for (std::size_t k = 0; k < N-1; ++k){
if (i==k) continue;
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->dalphar_dTau(HEOS._tau, HEOS._delta);
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->dalphar_dTau(HEOS._tau, HEOS._delta);
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->dalphar_dTau();
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->dalphar_dTau();
d2ar_dxi_dTau += x[k]*(Fikarik - FiNariN - FkNarkN);
}
return d2ar_dxi_dTau;
@@ -54,19 +54,19 @@ CoolPropDbl MixtureDerivatives::d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HE
CoolPropDbl MixtureDerivatives::d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return HEOS.components[i].EOS().dalphar_dDelta(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dDelta(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
return HEOS.components[i].EOS().dalphar_dDelta(HEOS.tau(), HEOS.delta()) + HEOS.Excess.d2alphar_dxi_dDelta(HEOS.mole_fractions, i);
}
else if(xN_flag == XN_DEPENDENT){
std::vector<CoolPropDbl> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return 0;
double d2ar_dxi_dDelta = HEOS.components[i].EOS().dalphar_dDelta(HEOS._tau, HEOS._delta) - HEOS.components[N-1].EOS().dalphar_dDelta(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dDelta(HEOS._tau, HEOS._delta);
double d2ar_dxi_dDelta = HEOS.components[i].EOS().dalphar_dDelta(HEOS.tau(), HEOS.delta()) - HEOS.components[N-1].EOS().dalphar_dDelta(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dDelta();
d2ar_dxi_dDelta += (1-2*x[i])*FiNariN;
for (std::size_t k = 0; k < N-1; ++k){
if (i==k) continue;
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->dalphar_dDelta(HEOS._tau, HEOS._delta);
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->dalphar_dDelta(HEOS._tau, HEOS._delta);
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->dalphar_dDelta();
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->dalphar_dDelta();
d2ar_dxi_dDelta += x[k]*(Fikarik - FiNariN - FkNarkN);
}
return d2ar_dxi_dDelta;
@@ -78,19 +78,19 @@ CoolPropDbl MixtureDerivatives::d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &
CoolPropDbl MixtureDerivatives::d3alphar_dxi_dDelta2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return HEOS.components[i].EOS().d2alphar_dDelta2(HEOS._tau, HEOS._delta) + HEOS.Excess.d3alphar_dxi_dDelta2(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
return HEOS.components[i].EOS().d2alphar_dDelta2(HEOS.tau(), HEOS.delta()) + HEOS.Excess.d3alphar_dxi_dDelta2(HEOS.mole_fractions, i);
}
else if (xN_flag == XN_DEPENDENT){
std::vector<CoolPropDbl> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return 0;
double d3ar_dxi_dDelta2 = HEOS.components[i].EOS().d2alphar_dDelta2(HEOS._tau, HEOS._delta) - HEOS.components[N-1].EOS().d2alphar_dDelta2(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->d2alphar_dDelta2(HEOS._tau, HEOS._delta);
double d3ar_dxi_dDelta2 = HEOS.components[i].EOS().d2alphar_dDelta2(HEOS.tau(), HEOS.delta()) - HEOS.components[N-1].EOS().d2alphar_dDelta2(HEOS.tau(), HEOS.delta());
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->d2alphar_dDelta2();
d3ar_dxi_dDelta2 += (1-2*x[i])*FiNariN;
for (std::size_t k = 0; k < N-1; ++k){
if (i==k) continue;
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->d2alphar_dDelta2(HEOS._tau, HEOS._delta);
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->d2alphar_dDelta2(HEOS._tau, HEOS._delta);
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->d2alphar_dDelta2();
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->d2alphar_dDelta2();
d3ar_dxi_dDelta2 += x[k]*(Fikarik - FiNariN - FkNarkN);
}
return d3ar_dxi_dDelta2;
@@ -102,19 +102,19 @@ CoolPropDbl MixtureDerivatives::d3alphar_dxi_dDelta2(HelmholtzEOSMixtureBackend
CoolPropDbl MixtureDerivatives::d3alphar_dxi_dTau2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return HEOS.components[i].EOS().d2alphar_dTau2(HEOS._tau, HEOS._delta) + HEOS.Excess.d3alphar_dxi_dTau2(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
return HEOS.components[i].EOS().d2alphar_dTau2(HEOS.tau(), HEOS.delta()) + HEOS.Excess.d3alphar_dxi_dTau2(HEOS.mole_fractions, i);
}
else if (xN_flag == XN_DEPENDENT){
std::vector<CoolPropDbl> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return 0;
double d3ar_dxi_dTau2 = HEOS.components[i].EOS().d2alphar_dTau2(HEOS._tau, HEOS._delta) - HEOS.components[N-1].EOS().d2alphar_dTau2(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->d2alphar_dTau2(HEOS._tau, HEOS._delta);
double d3ar_dxi_dTau2 = HEOS.components[i].EOS().d2alphar_dTau2(HEOS.tau(), HEOS.delta()) - HEOS.components[N-1].EOS().d2alphar_dTau2(HEOS.tau(), HEOS.delta());
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->d2alphar_dTau2();
d3ar_dxi_dTau2 += (1-2*x[i])*FiNariN;
for (std::size_t k = 0; k < N-1; ++k){
if (i==k) continue;
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->d2alphar_dTau2(HEOS._tau, HEOS._delta);
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->d2alphar_dTau2(HEOS._tau, HEOS._delta);
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->d2alphar_dTau2();
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->d2alphar_dTau2();
d3ar_dxi_dTau2 += x[k]*(Fikarik - FiNariN - FkNarkN);
}
return d3ar_dxi_dTau2;
@@ -126,19 +126,19 @@ CoolPropDbl MixtureDerivatives::d3alphar_dxi_dTau2(HelmholtzEOSMixtureBackend &H
CoolPropDbl MixtureDerivatives::d3alphar_dxi_dDelta_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return HEOS.components[i].EOS().d2alphar_dDelta_dTau(HEOS._tau, HEOS._delta) + HEOS.Excess.d3alphar_dxi_dDelta_dTau(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i);
return HEOS.components[i].EOS().d2alphar_dDelta_dTau(HEOS.tau(), HEOS.delta()) + HEOS.Excess.d3alphar_dxi_dDelta_dTau(HEOS.mole_fractions, i);
}
else if (xN_flag == XN_DEPENDENT){
std::vector<CoolPropDbl> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return 0;
double d3ar_dxi_dDelta_dTau = HEOS.components[i].EOS().d2alphar_dDelta_dTau(HEOS._tau, HEOS._delta) - HEOS.components[N-1].EOS().d2alphar_dDelta_dTau(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->d2alphar_dDelta_dTau(HEOS._tau, HEOS._delta);
double d3ar_dxi_dDelta_dTau = HEOS.components[i].EOS().d2alphar_dDelta_dTau(HEOS.tau(), HEOS.delta()) - HEOS.components[N-1].EOS().d2alphar_dDelta_dTau(HEOS.tau(), HEOS.delta());
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->d2alphar_dDelta_dTau();
d3ar_dxi_dDelta_dTau += (1-2*x[i])*FiNariN;
for (std::size_t k = 0; k < N-1; ++k){
if (i==k) continue;
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau(HEOS._tau, HEOS._delta);
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->d2alphar_dDelta_dTau(HEOS._tau, HEOS._delta);
double Fikarik = HEOS.Excess.F[i][k]*HEOS.Excess.DepartureFunctionMatrix[i][k]->d2alphar_dDelta_dTau();
double FkNarkN = HEOS.Excess.F[k][N-1]*HEOS.Excess.DepartureFunctionMatrix[k][N-1]->d2alphar_dDelta_dTau();
d3ar_dxi_dDelta_dTau += x[k]*(Fikarik - FiNariN - FkNarkN);
}
return d3ar_dxi_dDelta_dTau;
@@ -151,16 +151,16 @@ CoolPropDbl MixtureDerivatives::d3alphar_dxi_dDelta_dTau(HelmholtzEOSMixtureBack
CoolPropDbl MixtureDerivatives::d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return 0 + HEOS.Excess.d2alphardxidxj(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i, j);
return 0 + HEOS.Excess.d2alphardxidxj(HEOS.mole_fractions, i, j);
}
else if(xN_flag == XN_DEPENDENT){
std::size_t N = HEOS.mole_fractions.size();
if (i == N-1 || j == N-1){ return 0;}
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->alphar(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->alphar();
if (i == j) { return -2*FiNariN; }
double Fijarij = HEOS.Excess.F[i][j]*HEOS.Excess.DepartureFunctionMatrix[i][j]->alphar(HEOS._tau, HEOS._delta);
double FjNarjN = HEOS.Excess.F[j][N-1]*HEOS.Excess.DepartureFunctionMatrix[j][N-1]->alphar(HEOS._tau, HEOS._delta);
double Fijarij = HEOS.Excess.F[i][j]*HEOS.Excess.DepartureFunctionMatrix[i][j]->alphar();
double FjNarjN = HEOS.Excess.F[j][N-1]*HEOS.Excess.DepartureFunctionMatrix[j][N-1]->alphar();
return Fijarij - FiNariN - FjNarjN;
}
else{
@@ -170,15 +170,15 @@ CoolPropDbl MixtureDerivatives::d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS,
CoolPropDbl MixtureDerivatives::d3alphar_dxi_dxj_dDelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return 0 + HEOS.Excess.d3alphar_dxi_dxj_dDelta(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i, j);
return 0 + HEOS.Excess.d3alphar_dxi_dxj_dDelta(HEOS.mole_fractions, i, j);
}
else if (xN_flag == XN_DEPENDENT){
std::size_t N = HEOS.mole_fractions.size();
if (i == N-1 || j == N-1){ return 0; }
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dDelta(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dDelta();
if (i == j) { return -2*FiNariN; }
double Fijarij = HEOS.Excess.F[i][j]*HEOS.Excess.DepartureFunctionMatrix[i][j]->dalphar_dDelta(HEOS._tau, HEOS._delta);
double FjNarjN = HEOS.Excess.F[j][N-1]*HEOS.Excess.DepartureFunctionMatrix[j][N-1]->dalphar_dDelta(HEOS._tau, HEOS._delta);
double Fijarij = HEOS.Excess.F[i][j]*HEOS.Excess.DepartureFunctionMatrix[i][j]->dalphar_dDelta();
double FjNarjN = HEOS.Excess.F[j][N-1]*HEOS.Excess.DepartureFunctionMatrix[j][N-1]->dalphar_dDelta();
return Fijarij - FiNariN - FjNarjN;
}
else{
@@ -188,15 +188,15 @@ CoolPropDbl MixtureDerivatives::d3alphar_dxi_dxj_dDelta(HelmholtzEOSMixtureBacke
CoolPropDbl MixtureDerivatives::d3alphar_dxi_dxj_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return 0 + HEOS.Excess.d3alphar_dxi_dxj_dTau(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i, j);
return 0 + HEOS.Excess.d3alphar_dxi_dxj_dTau(HEOS.mole_fractions, i, j);
}
else if (xN_flag == XN_DEPENDENT){
std::size_t N = HEOS.mole_fractions.size();
if (i == N-1 || j == N-1){ return 0; }
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dTau(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->dalphar_dTau();
if (i == j) { return -2*FiNariN; }
double Fijarij = HEOS.Excess.F[i][j]*HEOS.Excess.DepartureFunctionMatrix[i][j]->dalphar_dTau(HEOS._tau, HEOS._delta);
double FjNarjN = HEOS.Excess.F[j][N-1]*HEOS.Excess.DepartureFunctionMatrix[j][N-1]->dalphar_dTau(HEOS._tau, HEOS._delta);
double Fijarij = HEOS.Excess.F[i][j]*HEOS.Excess.DepartureFunctionMatrix[i][j]->dalphar_dTau();
double FjNarjN = HEOS.Excess.F[j][N-1]*HEOS.Excess.DepartureFunctionMatrix[j][N-1]->dalphar_dTau();
return Fijarij - FiNariN - FjNarjN;
}
else{
@@ -206,7 +206,7 @@ CoolPropDbl MixtureDerivatives::d3alphar_dxi_dxj_dTau(HelmholtzEOSMixtureBackend
CoolPropDbl MixtureDerivatives::d3alphardxidxjdxk(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
return 0 + HEOS.Excess.d3alphardxidxjdxk(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i, j, k);
return 0 + HEOS.Excess.d3alphardxidxjdxk(HEOS.mole_fractions, i, j, k);
}
else if (xN_flag == XN_DEPENDENT){
return 0;