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

Yes!! Composition derivatives working almost entirely.

Browse Source
This commit is contained in:
Ian Bell
2014-08-28 21:52:25 +02:00
parent e7aa3df917
commit d539f5cde8
9 changed files with 400 additions and 184 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/Backends/Helmholtz/FlashRoutines.cpp
View File

@@ -241,7 +241,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
// Set some imput options
SaturationSolvers::mixture_VLE_IO io;
io.sstype = SaturationSolvers::imposed_p;
io.Nstep_max = 5;
io.Nstep_max = 2;
// Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted
long double Tguess = SaturationSolvers::saturation_preconditioner(HEOS, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions);
@@ -258,7 +258,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS)
IO.rhomolar_vap = io.rhomolar_vap;
IO.T = io.T;
IO.p = io.p;
IO.Nstep_max = 25;
IO.Nstep_max = 100;
// Dewpoint
NR.call(HEOS, io.y, io.x, IO);

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

@@ -10,7 +10,7 @@ long double MixtureDerivatives::dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, st
else if(xN_flag == XN_DEPENDENT){
std::vector<long double> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return _HUGE;
if (i==N-1) return 0;
double dar_dxi = HEOS.components[i]->pEOS->baser(HEOS._tau, HEOS._delta) - HEOS.components[N-1]->pEOS->baser(HEOS._tau, HEOS._delta);
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->alphar(HEOS._tau, HEOS._delta);
dar_dxi += (1-2*x[i])*FiNariN;
@@ -34,7 +34,7 @@ long double MixtureDerivatives::d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HE
else if(xN_flag == XN_DEPENDENT){
std::vector<long double> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return _HUGE;
if (i==N-1) return 0;
double d2ar_dxi_dTau = HEOS.components[i]->pEOS->dalphar_dTau(HEOS._tau, HEOS._delta) - HEOS.components[N-1]->pEOS->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);
d2ar_dxi_dTau += (1-2*x[i])*FiNariN;
@@ -59,7 +59,7 @@ long double MixtureDerivatives::d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &
else if(xN_flag == XN_DEPENDENT){
std::vector<long double> &x = HEOS.mole_fractions;
std::size_t N = x.size();
if (i==N-1) return _HUGE;
if (i==N-1) return 0;
double d2ar_dxi_dDelta = HEOS.components[i]->pEOS->dalphar_dDelta(HEOS._tau, HEOS._delta) - HEOS.components[N-1]->pEOS->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);
d2ar_dxi_dDelta += (1-2*x[i])*FiNariN;
@@ -85,8 +85,8 @@ long double MixtureDerivatives::d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS,
std::size_t N = HEOS.mole_fractions.size();
if (i == N-1){ return 0;}
double FiNariN = HEOS.Excess.F[i][N-1]*HEOS.Excess.DepartureFunctionMatrix[i][N-1]->alphar(HEOS._tau, HEOS._delta);
if (i == j) { return FiNariN; }
if (j == N-1){ return _HUGE; }
if (i == j) { return -2*FiNariN; }
if (j == N-1){ return 0; }
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);
@@ -252,18 +252,19 @@ long double MixtureDerivatives::ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend
long double MixtureDerivatives::d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
double line1 = HEOS._delta.pt()*d2alphar_dxi_dDelta(HEOS, j, xN_flag)*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i, xN_flag));
double line2 = -HEOS._delta.pt()*HEOS.dalphar_dDelta()*(1/HEOS._reducing.rhomolar)*(HEOS.Reducing.p->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->drhormolardxi__constxj(HEOS.mole_fractions,j, xN_flag)*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i, xN_flag));
double line3 = HEOS._tau.pt()*d2alphar_dxi_dTau(HEOS, j, xN_flag)*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag);
double line2 = -HEOS._delta.pt()*HEOS.dalphar_dDelta()*(1/HEOS._reducing.rhomolar)*(HEOS.Reducing.p->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j, xN_flag)-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->drhormolardxi__constxj(HEOS.mole_fractions,j, xN_flag)*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i, xN_flag));
double line4 = HEOS._tau.pt()*HEOS.dalphar_dTau()*(1/HEOS._reducing.T)*(HEOS.Reducing.p->d_ndTrdni_dxj__constxi(HEOS.mole_fractions,i,j, xN_flag)-1/HEOS._reducing.T*HEOS.Reducing.p->dTrdxi__constxj(HEOS.mole_fractions, j, xN_flag)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i, xN_flag));
double s = 0;
std::size_t mmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ mmax--; }
for (unsigned int m = 0; m < mmax; m++)
std::size_t kmax = HEOS.mole_fractions.size();
if (xN_flag == XN_DEPENDENT){ kmax--; }
for (unsigned int k = 0; k < kmax; k++)
{
s += HEOS.mole_fractions[m]*d2alphardxidxj(HEOS, j,m, xN_flag);
s += HEOS.mole_fractions[k]*d2alphardxidxj(HEOS, j, k, xN_flag);
}
double line5 = d2alphardxidxj(HEOS, i,j, xN_flag)-dalphar_dxi(HEOS, j, xN_flag)-s;
return line1+line2+line3+line4+line5;
double line5 = d2alphardxidxj(HEOS, i, j, xN_flag)-dalphar_dxi(HEOS, j, xN_flag)-s;
return line1 + line2 + line3 + line4 + line5;
}
long double MixtureDerivatives::nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
@@ -397,6 +398,27 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss3a;
ss3a << "d2alphar_dxi_dDelta, i=" << i;
SECTION(ss3a.str(), "")
{
if (i==1){break;}
double p1 = 101325, dP = 1e-1;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, p1, 300);
double analytic = MixtureDerivatives::d2alphar_dxi_dDelta(rHEOS, i, xN_flag);
rHEOS.update(PT_INPUTS, p1 + dP, 300);
double v1 = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag), delta1 = rHEOS.delta();
rHEOS.update(PT_INPUTS, p1 - dP, 300);
double v2 = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag), delta2 = rHEOS.delta();
double numeric = (v1 - v2)/(delta1 - delta2);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss4;
ss4 << "d_ndalphardni_dTau, i=" << i;
SECTION(ss4.str(), "")
@@ -418,10 +440,33 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss4a;
ss4a << "d2alphar_dxi_dTau, i=" << i;
SECTION(ss4a.str(), "")
{
if (i == 1){ break; }
double p1 = 101325, dT = 1e-2;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(PT_INPUTS, 101325, 300);
double rho1 = rHEOS.rhomolar();
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double analytic = MixtureDerivatives::d2alphar_dxi_dTau(rHEOS, i, xN_flag);
rHEOS.update(DmolarT_INPUTS, rho1, 300 + dT);
double v1 = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag), tau1 = rHEOS.tau();
rHEOS.update(DmolarT_INPUTS, rho1, 300 - dT);
double v2 = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag), tau2 = rHEOS.tau();
double numeric = (v1 - v2)/(tau1 - tau2);
double err = std::abs((numeric-analytic)/analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss5;
ss5 << "dpdxj__constT_V_xi, i=" << i;
SECTION(ss5.str(), "")
{
if (i==1){break;}
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
@@ -429,12 +474,12 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
double rho1 = rHEOS.rhomolar();
double analytic = MixtureDerivatives::dpdxj__constT_V_xi(rHEOS, i, xN_flag);
std::vector<long double> zp = z; /// Copy base composition
zp[i] += dz;
zp[i] += dz; zp[1-i] -= dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v1 = rHEOS.p();
std::vector<long double> zm = z; /// Copy base composition
zm[i] -= dz;
zm[i] -= dz; zm[1-i] += dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v2 = rHEOS.p();
@@ -444,13 +489,13 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-80);
}
CHECK(err < 1e-8);
}
std::ostringstream ss6;
ss6 << "d_dalpharddelta_dxj__constT_V_xi, i=" << i;
SECTION(ss6.str(), "")
{
if (i==1){break;}
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
@@ -460,14 +505,14 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
// Increment mole fraction
std::vector<long double> zp = z; /// Copy base composition
zp[i] += dz;
zp[i] += dz; zp[1-i] -= dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v1 = rHEOS.dalphar_dDelta();
// Decrement mole fraction
std::vector<long double> zm = z; /// Copy base composition
zm[i] -= dz;
zm[i] -= dz; zm[1-i] += dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v2 = rHEOS.dalphar_dDelta();
@@ -477,15 +522,15 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-60);
}
CHECK(err < 1e-8);
}
// These derivatives depend on both the i and j indices
for (std::size_t j = 0; j < z.size(); ++j){
std::ostringstream ss1;
ss1 << "dln_fugacity_coefficient_dxj__constT_p_xi, i=" << i << ", j=" << j;
SECTION(ss1.str(), "")
{
if (j == 1){break;}
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
@@ -495,12 +540,12 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
double rho1 = rHEOS.rhomolar();
double analytic = MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(rHEOS, i, j, xN_flag);
std::vector<long double> zp = z; /// Copy base composition
zp[j] += dz;
zp[j] += dz; zp[1-j] -= dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(PT_INPUTS, p, 300);
double v1 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag);
std::vector<long double> zm = z; /// Copy base composition
zm[j] -= dz;
zm[j] -= dz; zm[1-j] += dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(PT_INPUTS, p, 300);
double v2 = MixtureDerivatives::ln_fugacity_coefficient(rHEOS, i, xN_flag);
@@ -512,10 +557,77 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
CAPTURE(analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss3;
ss2 << "d_ndalphardni_dxj__constT_V_xi, i=" << i << ", j=" << j;
std::ostringstream ss2;
ss2 << "d_ndTrdni_dxj, i=" << i << ", j=" << j;
SECTION(ss2.str(), "")
{
if (j == 1){break;}
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
double rho1 = rHEOS.rhomolar();
double analytic = rHEOS.Reducing.p->d_ndTrdni_dxj__constxi(rHEOS.get_const_mole_fractions(), i, j, xN_flag);
// Increment mole fraction
std::vector<long double> zp = z; /// Copy base composition
zp[j] += dz; zp[1-j] -= dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v1 = rHEOS.Reducing.p->ndTrdni__constnj(rHEOS.get_const_mole_fractions(), i, xN_flag);
// Decrement mole fraction
std::vector<long double> zm = z; /// Copy base composition
zm[j] -= dz; zm[1-j] += dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v2 = rHEOS.Reducing.p->ndTrdni__constnj(rHEOS.get_const_mole_fractions(), i, xN_flag);
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss4;
ss4 << "d_ndrhomolarrdni_dxj, i=" << i << ", j=" << j;
SECTION(ss4.str(), "")
{
if (j == 1){break;}
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
double rho1 = rHEOS.rhomolar();
double analytic = rHEOS.Reducing.p->d_ndrhorbardni_dxj__constxi(rHEOS.get_const_mole_fractions(), i, j, xN_flag);
// Increment mole fraction
std::vector<long double> zp = z; /// Copy base composition
zp[j] += dz; zp[1-j] -= dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v1 = rHEOS.Reducing.p->ndrhorbardni__constnj(rHEOS.get_const_mole_fractions(), i, xN_flag);
// Decrement mole fraction
std::vector<long double> zm = z; /// Copy base composition
zm[j] -= dz; zm[1-j] += dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v2 = rHEOS.Reducing.p->ndrhorbardni__constnj(rHEOS.get_const_mole_fractions(), i, xN_flag);
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-8);
}
std::ostringstream ss3;
ss3 << "d_ndalphardni_dxj__constT_V_xi, i=" << i << ", j=" << j;
SECTION(ss3.str(), "")
{
if (j == 1){break;}
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
@@ -525,14 +637,14 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
// Increment mole fraction
std::vector<long double> zp = z; /// Copy base composition
zp[j] += dz;
zp[j] += dz; zp[1-j] -= dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v1 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i, xN_flag);
// Decrement mole fraction
std::vector<long double> zm = z; /// Copy base composition
zm[j] -= dz;
zm[j] -= dz; zm[1-j] -= dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v2 = MixtureDerivatives::ndalphar_dni__constT_V_nj(rHEOS, i, xN_flag);
@@ -544,7 +656,39 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]")
CAPTURE(analytic);
CHECK(err < 1e-60);
}
std::ostringstream ss3a;
ss3a << "d2alphardxidxj, i=" << i << ", j=" << j;
SECTION(ss3a.str(), "")
{
if (j == 1){break;}
double dz = 1e-6;
rHEOS.specify_phase(iphase_gas);
rHEOS.update(DmolarT_INPUTS, 300, 300);
double rho1 = rHEOS.rhomolar();
double analytic = MixtureDerivatives::d2alphardxidxj(rHEOS,i,j,xN_flag);
// Increment mole fraction
std::vector<long double> zp = z; /// Copy base composition
zp[j] += dz; zp[1-j] -= dz;
rHEOS.set_mole_fractions(zp);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v1 = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag);
// Decrement mole fraction
std::vector<long double> zm = z; /// Copy base composition
zm[j] -= dz; zm[1-j] -= dz;
rHEOS.set_mole_fractions(zm);
rHEOS.update(DmolarT_INPUTS, rho1, 300);
double v2 = MixtureDerivatives::dalphar_dxi(rHEOS, i, xN_flag);
double numeric = (v1 - v2)/(2*dz);
double err = std::abs((numeric-analytic)/analytic);
CAPTURE(numeric);
CAPTURE(analytic);
CHECK(err < 1e-60);
}
}
}

228
src/Backends/Helmholtz/ReducingFunctions.cpp
View File

@@ -146,12 +146,28 @@ ReducingFunction *ReducingFunction::factory(const std::vector<CoolPropFluid*> &c
long double ReducingFunction::d_ndTrdni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
long double s = 0;
for (std::size_t k = 0; k < N; k++)
{
s += x[k]*d2Trdxidxj(x,j,k, xN_flag);
}
return d2Trdxidxj(x,i,j, xN_flag)-dTrdxi__constxj(x,j, xN_flag)-s;
if (xN_flag == XN_INDEPENDENT){
long double s = 0;
for (std::size_t k = 0; k < N; k++)
{
s += x[k]*d2Trdxidxj(x,j,k, xN_flag);
}
return d2Trdxidxj(x,i,j, xN_flag)-dTrdxi__constxj(x,j, xN_flag)-s;
}
else if (xN_flag == XN_DEPENDENT){
if (j == N-1){ return 0;}
long double s = 0;
for (std::size_t k = 0; k < N-1; k++)
{
s += x[k]*d2Trdxidxj(x,j,k, xN_flag);
}
long double val = d2Trdxidxj(x,i,j, xN_flag)-dTrdxi__constxj(x,j, xN_flag)-s;
return val;
}
else{
throw ValueError(format("xN dependency flag invalid"));
}
}
long double ReducingFunction::d_ndrhorbardni_dxj__constxi(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
@@ -217,6 +233,14 @@ long double GERG2008ReducingFunction::dTrdxi__constxj(const std::vector<long dou
{
return dYrdxi__constxj(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
}
long double GERG2008ReducingFunction::d2Trdxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
return d2Yrdxi2__constxj(x, i, beta_T, gamma_T, T_c, Yc_T, xN_flag);
}
long double GERG2008ReducingFunction::d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
return d2Yrdxidxj__constxj(x, i, j, beta_T, gamma_T, T_c, Yc_T, xN_flag);
}
long double GERG2008ReducingFunction::rhormolar(const std::vector<long double> &x)
{
return 1/Yr(x, beta_v, gamma_v, v_c, Yc_v);
@@ -229,6 +253,14 @@ long double GERG2008ReducingFunction::dvrmolardxi__constxj(const std::vector<lon
{
return dYrdxi__constxj(x, i, beta_v, gamma_v, v_c, Yc_v, xN_flag);
}
long double GERG2008ReducingFunction::d2vrmolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
return d2Yrdxi2__constxj(x, i, beta_v, gamma_v, v_c, Yc_v, xN_flag);
}
long double GERG2008ReducingFunction::d2vrmolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
return d2Yrdxidxj__constxj(x, i, j, beta_v, gamma_v, v_c, Yc_v, xN_flag);
}
long double GERG2008ReducingFunction::d2rhormolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
long double rhor = this->rhormolar(x);
@@ -243,105 +275,22 @@ long double GERG2008ReducingFunction::d2rhormolardxidxj(const std::vector<long d
return 2*pow(rhor,(std::size_t)3)*dvrbardxi*dvrbardxj-pow(rhor,(std::size_t)2)*this->d2vrmolardxidxj(x,i,j, xN_flag);
}
long double GERG2008ReducingFunction::d2Trdxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
long double d2Tr_dxi2 = 2*pFluids[i]->pEOS->reduce.T;
for (std::size_t k = 0; k < i; k++)
{
d2Tr_dxi2 += c_Y_ij(k,i,beta_T,gamma_T,T_c)*d2fYkidxi2__constxk(x,k,i,beta_T);
}
for (std::size_t k = i+1; k < N; k++)
{
d2Tr_dxi2 += c_Y_ij(i,k,beta_T,gamma_T,T_c)*d2fYikdxi2__constxk(x,i,k,beta_T);
}
return d2Tr_dxi2;
}
long double GERG2008ReducingFunction::d2Trdxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
if (i == j)
{
return d2Trdxi2__constxj(x, i, xN_flag);
}
else
{
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
return c_Y_ij(i, j, beta_T, gamma_T, T_c)*d2fYijdxidxj(x, i, j, beta_T);
}
}
long double GERG2008ReducingFunction::d2vrmolardxidxj(const std::vector<long double> &x, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag)
{
if (i == j)
{
return d2vrmolardxi2__constxj(x, i, xN_flag);
}
else
{
return c_Y_ij(i, j, beta_v, gamma_v, v_c)*d2fYijdxidxj(x, i, j, beta_v);
}
}
long double GERG2008ReducingFunction::d2vrmolardxi2__constxj(const std::vector<long double> &x, std::size_t i, x_N_dependency_flag xN_flag)
{
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
double d2vrbardxi2 = 2/pFluids[i]->pEOS->reduce.rhomolar;
for (std::size_t k = 0; k < i; k++)
{
d2vrbardxi2 += c_Y_ij(k, i, beta_v, gamma_v, v_c)*d2fYkidxi2__constxk(x, k, i, beta_v);
}
for (std::size_t k = i+1; k < N; k++)
{
d2vrbardxi2 += c_Y_ij(i, k, beta_v, gamma_v, v_c)*d2fYikdxi2__constxk(x, i, k, beta_v);
}
return d2vrbardxi2;
}
long double GERG2008ReducingFunction::Yr(const std::vector<long double> &x, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc)
{
long double Yr = 0;
// if (xN_flag == XN_INDEPENDENT)
// {
for (std::size_t i = 0; i < N; i++)
{
double xi = x[i];
Yr += xi*xi*Yc[i];
// The last term is only used for the pure component, as it is sum_{i=1}^{N-1}sum_{j=1}^{N}
if (i==N-1){ break; }
for (std::size_t i = 0; i < N; i++)
{
double xi = x[i];
Yr += xi*xi*Yc[i];
// The last term is only used for the pure component, as it is sum_{i=1}^{N-1}sum_{j=1}^{N}
if (i==N-1){ break; }
for (std::size_t j = i+1; j < N; j++)
{
Yr += c_Y_ij(i, j, beta, gamma, Y_c_ij)*f_Y_ij(x, i, j, beta);
}
for (std::size_t j = i+1; j < N; j++)
{
Yr += c_Y_ij(i, j, beta, gamma, Y_c_ij)*f_Y_ij(x, i, j, beta);
}
// }
// else if (xN_flag == XN_DEPENDENT){
//
// // A.43 from Gernert, 2014, supplemental information
// for (std::size_t k = 0; k < N-1; ++k)
// {
// double xk = x[k];
// Yr += xk*xk*Yc[k];
// }
// for (std::size_t k = 0; k < N-2; ++k)
// {
// for (std::size_t m = k+1; m < N-1; ++m)
// {
// Yr += c_Y_ij(k, m, beta, gamma, Y_c_ij)*f_Y_ij(x, k, m, beta);
// }
// }
// for (std::size_t k = 0; k < N-1; ++k)
// {
// Yr += c_Y_ij(k, N-1, beta, gamma, Y_c_ij)*f_Y_ij(x, k, N-1, beta);
// }
// double xN = x[N-1];
// Yr += xN*xN*Yc[N-1];
// }
// else{
// throw ValueError(format("xN dependency flag invalid"));
// }
}
return Yr;
}
long double GERG2008ReducingFunction::dYrdxi__constxj(const std::vector<long double> &x, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc, x_N_dependency_flag xN_flag)
@@ -362,6 +311,7 @@ long double GERG2008ReducingFunction::dYrdxi__constxj(const std::vector<long dou
}
else if (xN_flag == XN_DEPENDENT){
// Table S1 from Gernert, 2014, supplemental information
if (i == N-1){return 0.0;}
long double dYr_dxi = 2*x[i]*Yc[i] - 2*x[N-1]*Yc[N-1];
for (std::size_t k = 0; k < i; k++)
{
@@ -384,6 +334,84 @@ long double GERG2008ReducingFunction::dYrdxi__constxj(const std::vector<long dou
throw ValueError(format("xN dependency flag invalid"));
}
}
long double GERG2008ReducingFunction::d2Yrdxi2__constxj(const std::vector<long double> &x, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
long double d2Yr_dxi2 = 2*Yc[i];
for (std::size_t k = 0; k < i; k++)
{
d2Yr_dxi2 += c_Y_ij(k,i,beta,gamma,Y_c_ij)*d2fYkidxi2__constxk(x,k,i,beta);
}
for (std::size_t k = i+1; k < N; k++)
{
d2Yr_dxi2 += c_Y_ij(i,k,beta,gamma,Y_c_ij)*d2fYikdxi2__constxk(x,i,k,beta);
}
return d2Yr_dxi2;
}
else if (xN_flag == XN_DEPENDENT){
// Table S1 from Gernert, 2014, supplemental information
if (i == N-1){return 0.0;}
long double d2Yr_dxi2 = 2*Yc[i] + 2*Yc[N-1];
for (std::size_t k = 0; k < i; k++)
{
d2Yr_dxi2 += c_Y_ji(k, i, beta, gamma, Y_c_ij)*d2fYkidxi2__constxk(x,k,i,beta);
}
for (std::size_t k = i+1; k < N-1; k++)
{
d2Yr_dxi2 += c_Y_ji(i, k, beta, gamma, Y_c_ij)*d2fYkidxi2__constxk(x,i,k,beta);
}
double beta_Y_iN = beta[i][N-1], xN = x[N-1];
d2Yr_dxi2 += 2*c_Y_ji(i, N-1, beta, gamma, Y_c_ij)*(-(x[i]+xN)/(pow(beta_Y_iN,2)*x[i]+xN)+(1-beta_Y_iN*beta_Y_iN)*(xN*xN/pow(beta_Y_iN*beta_Y_iN*x[i]+xN, 2)+((1-beta_Y_iN*beta_Y_iN)*x[i]*xN*xN-beta_Y_iN*beta_Y_iN*x[i]*x[i]*xN)/pow(beta_Y_iN*beta_Y_iN*x[i]+xN, 3)));
for (std::size_t k = 0; k < N-1; ++k)
{
double beta_Y_kN = beta[k][N-1];
d2Yr_dxi2 += 2*c_Y_ji(k, N-1, beta, gamma, Y_c_ij)*x[k]*x[k]*(1-beta_Y_kN*beta_Y_kN)/pow(pow(beta_Y_kN,2)*x[k]+xN, 2)*(xN/(beta_Y_kN*beta_Y_kN*x[k]+xN)-1);
}
return d2Yr_dxi2;
}
else{
throw ValueError(format("xN dependency flag invalid"));
}
}
long double GERG2008ReducingFunction::d2Yrdxidxj__constxj(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc, x_N_dependency_flag xN_flag)
{
if (xN_flag == XN_INDEPENDENT){
if (i == j)
{
return d2Yrdxi2__constxj(x, i, beta, gamma, Y_c_ij, Yc, xN_flag);
}
else
{
// See Table B9 from Kunz Wagner 2012 (GERG 2008)
return c_Y_ij(i, j, beta, gamma, Y_c_ij)*d2fYijdxidxj(x, i, j, beta);
}
}
else if (xN_flag == XN_DEPENDENT){
// Table S1 from Gernert, 2014, supplemental information
if (j == N-1 || i == N-1){ return 0.0; }
if (i == j){ return d2Yrdxi2__constxj(x, i, beta, gamma, Y_c_ij, Yc, xN_flag); }
long double xN = x[N-1];
long double d2Yr_dxidxj = 2*Yc[N-1];
d2Yr_dxidxj += c_Y_ji(i, j, beta, gamma, Y_c_ij)*d2fYijdxidxj(x,i,j,beta);
for (std::size_t k = 0; k < N-1; k++)
{
double beta_Y_kN = beta[k][N-1];
d2Yr_dxidxj += 2*c_Y_ji(k, N-1, beta, gamma, Y_c_ij)*x[k]*x[k]*(1-beta_Y_kN*beta_Y_kN)/pow(beta_Y_kN*beta_Y_kN*x[k]+xN,2)*(xN/(beta_Y_kN*beta_Y_kN*x[k]+xN)-1);
}
double beta_Y_iN = beta[i][N-1], beta_Y_jN = beta[j][N-1];
d2Yr_dxidxj += c_Y_ji(i, N-1, beta, gamma, Y_c_ij)*((1-beta_Y_iN*beta_Y_iN)*(2*x[i]*xN*xN/pow(beta_Y_iN*beta_Y_iN*x[i]+xN, 3)-x[i]*xN/pow(beta_Y_iN*beta_Y_iN*x[i]+xN, 2)) - (x[i]+xN)/(beta_Y_iN*beta_Y_iN*x[i]+xN));
d2Yr_dxidxj += -c_Y_ji(j, N-1, beta, gamma, Y_c_ij)*((1-beta_Y_jN*beta_Y_jN)*(2*x[j]*x[j]*xN*beta_Y_jN*beta_Y_jN/pow(beta_Y_jN*beta_Y_jN*x[j]+xN, 3)-x[j]*xN/pow(beta_Y_jN*beta_Y_jN*x[j]+xN, 2)) - (x[j]+xN)/(beta_Y_jN*beta_Y_jN*x[j]+xN));
return d2Yr_dxidxj;
}
else{
throw ValueError(format("xN dependency flag invalid"));
}
}
long double GERG2008ReducingFunction::dfYkidxi__constxk(const std::vector<long double> &x, std::size_t k, std::size_t i, const STLMatrix &beta)
{

7
src/Backends/Helmholtz/ReducingFunctions.h
View File

@@ -9,7 +9,8 @@ namespace CoolProp{
typedef std::vector<std::vector<long double> > STLMatrix;
enum x_N_dependency_flag{XN_INDEPENDENT, ///< x_N is an independent variable, and not calculated by \f$ x_N = 1-\sum_i x_i\f$
XN_DEPENDENT}; ///< x_N is an dependent variable, calculated by \f$ x_N = 1-\sum_i x_i\f$
XN_DEPENDENT ///< x_N is an dependent variable, calculated by \f$ x_N = 1-\sum_i x_i\f$
};
/// A container for the mixing parameters for CoolProp mixtures
/**
@@ -119,7 +120,7 @@ protected:
STLMatrix gamma_T; ///< \f$ \gamma_{T,ij} \f$ from GERG-2008
std::vector<long double> Yc_T; ///< Vector of critical temperatures for all components
std::vector<long double> Yc_v; ///< Vector of critical molar volumes for all components
std::vector<CoolPropFluid *> pFluids; ///< List of postd::size_ters to fluids
std::vector<CoolPropFluid *> pFluids; ///< List of pointer to fluids
public:
GERG2008ReducingFunction(std::vector<CoolPropFluid *> pFluids, STLMatrix beta_v, STLMatrix gamma_v, STLMatrix beta_T, STLMatrix gamma_T)
@@ -167,6 +168,8 @@ public:
long double Yr(const std::vector<long double> &x, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc);
long double dYrdxi__constxj(const std::vector<long double> &x, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc, x_N_dependency_flag xN_flag);
long double d2Yrdxi2__constxj(const std::vector<long double> &x, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc, x_N_dependency_flag xN_flag);
long double d2Yrdxidxj__constxj(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c_ij, const std::vector<long double> &Yc, x_N_dependency_flag xN_flag);
long double c_Y_ij(std::size_t i, std::size_t j, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c);
long double c_Y_ji(std::size_t j, std::size_t i, const STLMatrix &beta, const STLMatrix &gamma, const STLMatrix &Y_c);
long double f_Y_ij(const std::vector<long double> &x, std::size_t i, std::size_t j, const STLMatrix &beta);

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

@@ -1079,7 +1079,7 @@ void SaturationSolvers::newton_raphson_saturation::check_Jacobian()
// Derivatives with respect to x0
double dx = 1e-5, T = T0;
x = x0; x[0] += dx;// x[1] -= dx;
x = x0; x[0] += dx; x[1] -= dx;
rSatL.set_mole_fractions(x);
rSatL.update_TP_guessrho(T, p, rhomolar_liq); rhomolar_liq = rSatL.rhomolar();
build_arrays(); r1 = r;
@@ -1119,7 +1119,7 @@ void SaturationSolvers::newton_raphson_saturation::call(HelmholtzEOSMixtureBacke
// Build the Jacobian and residual vectors
build_arrays();
check_Jacobian();
//check_Jacobian();
// Solve for the step; v is the step with the contents
// [delta(x_0), delta(x_1), ..., delta(x_{N-2}), delta(spec)]

66
src/Backends/Helmholtz/VLERoutines.h
View File

@@ -233,16 +233,17 @@ namespace SaturationSolvers
rhomolar_liq = _HUGE; rhomolar_vap = _HUGE; T = _HUGE; p = _HUGE;
};
/*! Call the Newton-Raphson VLE Solver
This solver must be passed reasonable guess values for the mole fractions,
densities, etc. You may want to take a few steps of successive substitution
before you start with Newton Raphson.
@param HEOS HelmholtzEOSMixtureBackend instance
@param z Bulk mole fractions [-]
@param z_incipient Initial guesses for the mole fractions of the incipient phase [-]
*/
/** Call the Newton-Raphson VLE Solver
*
* This solver must be passed reasonable guess values for the mole fractions,
* densities, etc. You may want to take a few steps of successive substitution
* before you start with Newton Raphson.
*
* @param HEOS HelmholtzEOSMixtureBackend instance
* @param z Bulk mole fractions [-]
* @param z_incipient Initial guesses for the mole fractions of the incipient phase [-]
* @param IO The input/output data structure
*/
void call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &z_incipient, newton_raphson_saturation_options &IO);
/*! Build the arrays for the Newton-Raphson solve
@@ -287,30 +288,31 @@ namespace SaturationSolvers
rhobar_liq = _HUGE; rhobar_vap = _HUGE; T = _HUGE; p = _HUGE;
};
/*! Call the Newton-Raphson VLE Solver
This solver must be passed reasonable guess values for the mole fractions,
densities, etc. You may want to take a few steps of successive substitution
before you start with Newton Raphson.
@param HEOS Temperature [K]
@param z Pressure [Pa]
@param z Bulk mole fractions [-]
@param K Array of K-factors [-]
*/
/** Call the Newton-Raphson VLE Solver
*
* This solver must be passed reasonable guess values for the mole fractions,
* densities, etc. You may want to take a few steps of successive substitution
* before you start with Newton Raphson.
*
* @param HEOS Temperature [K]
* @param z Bulk mole fractions [-]
* @param K Array of K-factors [-]
* @param IO The input/output data structure
*/
void call(HelmholtzEOSMixtureBackend &HEOS, const std::vector<long double> &z, std::vector<long double> &K, mixture_VLE_IO &IO);
/*! Build the arrays for the Newton-Raphson solve
This method builds the Jacobian matrix, the sensitivity matrix, etc.
@param beta Void fraction [-] (0: bubble, 1: dew)
@param T Temperature [K]
@param p Pressure [Pa]
@param z Bulk mole fractions [-]
@param K Array of K-factors [-]
*/
void build_arrays(HelmholtzEOSMixtureBackend &HEOS, long double beta, long double T, long double rhomolar_liq, const long double rho_vapor, const std::vector<long double> &z, std::vector<long double> & K);
/** \brief Build the arrays for the Newton-Raphson solve
*
* This method builds the Jacobian matrix, the sensitivity matrix, etc.
*
* @param beta Void fraction [-] (0: bubble, 1: dew)
* @param T Temperature [K]
* @param rhomolar_liq Molar density of liquid [mol/m3]
* @param rhomolar_vap Molar density of liquid [mol/m3]
* @param z Bulk mole fractions [-]
* @param K Array of K-factors [-]
*/
void build_arrays(HelmholtzEOSMixtureBackend &HEOS, long double beta, long double T, long double rhomolar_liq, const long double rhomolar_vapor, const std::vector<long double> &z, std::vector<long double> & K);
/** Check the derivatives in the Jacobian using numerical derivatives.
*/