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

Merge pull request #778 from CoolProp/mixture_derivs

Browse Source 
Mixture derivs
This commit is contained in:
Ian Bell
2015-08-16 17:44:42 -06:00
parent ba7d5a9c97 30c5abc35c
commit 381f613164
5 changed files with 945 additions and 162 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -2932,34 +2932,46 @@ void HelmholtzEOSMixtureBackend::calc_critical_point(double rho0, double T0)
};
std::vector<std::vector<double> > Jacobian(const std::vector<double> &x)
{
double epsilon;
std::size_t N = x.size();
std::vector<double> r, xp;
std::vector<std::vector<double> > J(N, std::vector<double>(N, 0));
Eigen::MatrixXd J0(N, N), adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)),
dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau),
dLdDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta),
adjM = adjugate(MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT)),
dMdTau = dLdTau, dMdDelta = dLdDelta;
Eigen::MatrixXd adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)),
adjM = adjugate(MixtureDerivatives::Mstar(HEOS, XN_INDEPENDENT)),
dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau),
dLdDelta = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iDelta),
dMdTau = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iTau),
dMdDelta = MixtureDerivatives::dMstar_dX(HEOS, XN_INDEPENDENT, iDelta);
J0(0,0) = (adjL*dLdTau).trace();
J0(0,1) = (adjL*dLdDelta).trace();
//std::cout << J0 << std::endl;
std::vector<double> r0 = call(x);
J[0][0] = (adjL*dLdTau).trace();
J[0][1] = (adjL*dLdDelta).trace();
J[1][0] = (adjM*dMdTau).trace();
J[1][1] = (adjM*dMdDelta).trace();
return J;
}
/// Not used, for testing purposes
std::vector<std::vector<double> > numerical_Jacobian(const std::vector<double> &x)
{
std::size_t N = x.size();
std::vector<double> rplus, rminus, xp;
std::vector<std::vector<double> > J(N, std::vector<double>(N, 0));
double epsilon = 0.0001;
// Build the Jacobian by column
for (std::size_t i = 0; i < N; ++i)
{
xp = x;
epsilon = 0.00001;
xp[i] += epsilon;
r = call(xp);
rplus = call(xp);
xp[i] -= 2*epsilon;
rminus = call(xp);
for (std::size_t j = 0; j < N; ++j)
{
J[j][i] = (r[j]-r0[j])/epsilon;
J[j][i] = (rplus[j]-rminus[j])/(2*epsilon);
}
}
std::cout << J[0][0] << " " << J[0][1] << std::endl;
std::cout << J[1][0] << " " << J[1][1] << std::endl;
return J;
};
};
@@ -2987,13 +2999,11 @@ public:
return L1;
}
double deriv(double tau){
std::size_t N = HEOS.get_mole_fractions().size();
Eigen::MatrixXd adjL = adjugate(MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT)),
dLdTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau);
return (adjL*dLdTau).trace();
};
double second_deriv(double tau){
std::size_t N = HEOS.get_mole_fractions().size();
Eigen::MatrixXd Lstar = MixtureDerivatives::Lstar(HEOS, XN_INDEPENDENT),
dLstardTau = MixtureDerivatives::dLstar_dX(HEOS, XN_INDEPENDENT, iTau),
d2LstardTau2 = MixtureDerivatives::d2Lstar_dX2(HEOS, XN_INDEPENDENT, iTau, iTau),
@@ -3050,7 +3060,6 @@ public:
@param theta The angle
*/
double deriv(double theta){
std::size_t N = HEOS.get_mole_fractions().size();
double dL1_dtau = (adjLstar*dLstardTau).trace(), dL1_ddelta = (adjLstar*dLstardDelta).trace();
return -R_tau*sin(theta)*dL1_dtau + R_delta*cos(theta)*dL1_ddelta;
};

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

File diff suppressed because it is too large Load Diff

115
src/Backends/Helmholtz/MixtureDerivatives.h
View File

@@ -71,7 +71,6 @@ class MixtureDerivatives{
static CoolPropDbl d4alphar_dxi_dxj_dDelta2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d4alphar_dxi_dxj_dDelta_dTau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d4alphar_dxi_dxj_dTau2(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d4alphardxidxjdxk(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.63
@@ -162,9 +161,14 @@ class MixtureDerivatives{
static CoolPropDbl ndln_fugacity_i_dnj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d_ndln_fugacity_i_dnj_dtau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_ndln_fugacity_i_dnj_dtau2__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_ndln_fugacity_i_dnj_ddelta2__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d_ndln_fugacity_i_dnj_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl nAij(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){
@@ -175,6 +179,51 @@ class MixtureDerivatives{
CoolPropDbl RT = HEOS.gas_constant()*HEOS.T();
return 1/RT*nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HEOS, i, j, k, xN_flag) - nAij(HEOS, i, j, xN_flag);
}
static CoolPropDbl d_nAij_dX(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag, parameters WRT)
{
if (WRT == iTau){
return 1/(HEOS.gas_constant()*HEOS.T_reducing())*MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(HEOS, i, j, xN_flag)
+ 1/(HEOS.gas_constant()*HEOS.T())*MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag);
}
else if (WRT == iDelta){
return 1/(HEOS.gas_constant()*HEOS.T())*MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag);
}
else{
throw ValueError(format("wrong WRT"));
}
}
static CoolPropDbl d_n2Aijk_dX(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag, parameters WRT){
CoolPropDbl RT = HEOS.gas_constant()*HEOS.T();
if (WRT == iTau){
double summer = 0;
summer += d2_ndln_fugacity_i_dnj_dtau2__constdelta_x(HEOS, i, j, xN_flag)*ndtaudni__constT_V_nj(HEOS, k, xN_flag);
summer += d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag)*d_ndtaudni_dTau(HEOS, k, xN_flag);
summer += d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(HEOS, i, j, xN_flag)*nddeltadni__constT_V_nj(HEOS, k, xN_flag);
summer += d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(HEOS, i, j, k, xN_flag);
std::size_t mmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ mmax--; }
for (std::size_t m = 0; m < mmax; ++m){
summer -= HEOS.mole_fractions[m]*d2_ndln_fugacity_i_dnj_dxk_dTau__constdelta(HEOS, i, j, m, xN_flag);
}
return 1/RT*summer + 1/(HEOS.gas_constant()*HEOS.T_reducing())*nd_ndln_fugacity_i_dnj_dnk__constT_V_xi(HEOS,i,j,k,xN_flag) - d_nAij_dX(HEOS, i, j, xN_flag, WRT);
}
else if (WRT== iDelta){
double summer = 0;
summer += d2_ndln_fugacity_i_dnj_ddelta_dtau__constx(HEOS, i, j, xN_flag)*ndtaudni__constT_V_nj(HEOS, k, xN_flag);
summer += d2_ndln_fugacity_i_dnj_ddelta2__consttau_x(HEOS, i, j, xN_flag)*nddeltadni__constT_V_nj(HEOS, k, xN_flag);
summer += d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag)*d_nddeltadni_dDelta(HEOS, k, xN_flag);
summer += d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau(HEOS, i, j, k, xN_flag);
std::size_t mmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ mmax--; }
for (std::size_t m = 0; m < mmax; ++m){
summer -= HEOS.mole_fractions[m]*d2_ndln_fugacity_i_dnj_dxk_dDelta__consttau(HEOS, i, j, m, xN_flag);
}
return 1/RT*summer - d_nAij_dX(HEOS, i, j, xN_flag, iDelta);
}
else{
return _HUGE;
}
}
static Eigen::MatrixXd Lstar(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){
std::size_t N = HEOS.mole_fractions.size();
Eigen::MatrixXd L;
@@ -198,16 +247,7 @@ class MixtureDerivatives{
Eigen::MatrixXd dLstar_dX(N, N);
for (std::size_t i = 0; i < N; ++i){
for (std::size_t j = i; j < N; ++j){
if (WRT == iTau){
dLstar_dX(i, j) = 1/(HEOS.gas_constant()*HEOS.T_reducing())*MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(HEOS, i, j, xN_flag)
+ 1/(HEOS.gas_constant()*HEOS.T())*MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag);
}
else if (WRT == iDelta){
dLstar_dX(i, j) = 1/(HEOS.gas_constant()*HEOS.T())*MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag);
}
else{
throw ValueError(format("dLstar_dX invalid WRT"));
}
dLstar_dX(i, j) = d_nAij_dX(HEOS, i, j, xN_flag, WRT);
}
}
// Fill in the symmetric elements
@@ -244,7 +284,9 @@ class MixtureDerivatives{
static Eigen::MatrixXd Mstar(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){
std::size_t N = HEOS.mole_fractions.size();
Eigen::MatrixXd L = Lstar(HEOS, xN_flag), M = L;
Eigen::MatrixXd L = Lstar(HEOS, xN_flag),
M = L,
adjL = adjugate(L);
// Last row
for (std::size_t i = 0; i < N; ++i){
@@ -252,18 +294,39 @@ class MixtureDerivatives{
for (std::size_t j = 0; j < N; ++j){
for (std::size_t k = j; k < N; ++k){
n2dLdni(j, k) = n2Aijk(HEOS, j, k, i, xN_flag);
// Fill in the symmetric elements
n2dLdni(k, j) = n2dLdni(j, k);
}
}
// Fill in the symmetric elements
for (std::size_t j = 0; j < N; ++j){
for (std::size_t k = 0; k < j; ++k){
n2dLdni(j, k) = n2dLdni(k, j);
}
}
M(N-1, i) = (adjugate(L)*n2dLdni).trace();
M(N-1, i) = (adjL*n2dLdni).trace();
}
return M;
}
static Eigen::MatrixXd dMstar_dX(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag, parameters WRT){
std::size_t N = HEOS.mole_fractions.size();
Eigen::MatrixXd L = Lstar(HEOS, xN_flag),
dL_dX = dLstar_dX(HEOS, xN_flag, WRT),
dMstar = dL_dX,
adjL = adjugate(L),
d_adjL_dX = adjugate_derivative(L, dL_dX);
// Last row in the d(Mstar)/d(X) requires derivatives of
for (std::size_t i = 0; i < N; ++i){
Eigen::MatrixXd n2dLdni(N, N), d_n2dLdni_dX(N, N);
for (std::size_t j = 0; j < N; ++j){
for (std::size_t k = j; k < N; ++k){
n2dLdni(j, k) = n2Aijk(HEOS, j, k, i, xN_flag);
d_n2dLdni_dX(j, k) = d_n2Aijk_dX(HEOS, j, k, i, xN_flag, WRT);
// Fill in the symmetric elements
n2dLdni(k, j) = n2dLdni(j, k);
d_n2dLdni_dX(k, j) = d_n2dLdni_dX(j, k);
}
}
dMstar(N-1, i) = (n2dLdni*d_adjL_dX + adjL*d_n2dLdni_dX).trace();
}
return dMstar;
}
/** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions
*
@@ -497,6 +560,8 @@ class MixtureDerivatives{
*/
static CoolPropDbl d2_ndalphardni_dxj_dTau__constdelta_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d3_ndalphardni_dxj_dTau2__constdelta_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d3_ndalphardni_dxj_dDelta2__consttau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d3_ndalphardni_dxj_dDelta_dTau__constxi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.41
*
@@ -540,6 +605,8 @@ class MixtureDerivatives{
*/
static CoolPropDbl d_nd_ndalphardni_dnj_dTau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_nd_ndalphardni_dnj_dTau2__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_nd_ndalphardni_dnj_dDelta2__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_nd_ndalphardni_dnj_dDelta_dTau__constx(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/* \brief \f$\delta\f$ derivative of GERG 2004 7.47
*
@@ -550,6 +617,9 @@ class MixtureDerivatives{
*
*/
static CoolPropDbl d_nd_ndalphardni_dnj_dxk__consttau_delta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_nd_ndalphardni_dnj_dxk_dTau__constdelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_nd_ndalphardni_dnj_dxk_dDelta__consttau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.48
*
@@ -567,6 +637,8 @@ class MixtureDerivatives{
static CoolPropDbl d_nddeltadni_dxj__constdelta_tau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_nddeltadni_dxj_dDelta__consttau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.49
*
* The derivative term
@@ -583,6 +655,8 @@ class MixtureDerivatives{
static CoolPropDbl d_ndtaudni_dxj__constdelta_tau(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
static CoolPropDbl d2_ndtaudni_dxj_dTau__constdelta(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag);
/** \brief GERG 2004 Monograph equation 7.52
*
* The derivative term
@@ -605,6 +679,9 @@ class MixtureDerivatives{
*/
static CoolPropDbl d2_ndalphardni_dxj_dxk__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl d3_ndalphardni_dxj_dxk_dTau__constdelta_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
static CoolPropDbl d3_ndalphardni_dxj_dxk_dDelta__consttau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag);
}; /* class MixtureDerivatives */
} /* namepsace CoolProp*/