From 3fd77d57ba7819663a1bd36c6600c47c19f65cd2 Mon Sep 17 00:00:00 2001 From: Ian Bell Date: Sat, 30 May 2015 23:49:04 -0600 Subject: [PATCH] Implement all the derivatives needed for critical point determinant --- src/Backends/Helmholtz/MixtureDerivatives.cpp | 49 ++++++++---- src/Backends/Helmholtz/MixtureDerivatives.h | 36 +++++++++ src/Backends/Helmholtz/ReducingFunctions.cpp | 76 ++++++++++--------- 3 files changed, 113 insertions(+), 48 deletions(-) diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp index f7265631..01cd2064 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.cpp +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -328,6 +328,31 @@ CoolPropDbl MixtureDerivatives::dln_fugacity_dxj__constT_rho_xi(HelmholtzEOSMixt return line1 + line2 + line3 + line4; } +CoolPropDbl MixtureDerivatives::ndln_fugacity_i_dnj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ + return Kronecker_delta(i, j)/HEOS.mole_fractions[i]+nd2nalphardnidnj__constT_V(HEOS, i, j, xN_flag); +} +CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_dtau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ + return d_ndalphardni_dTau(HEOS, j, xN_flag) + d_nd_ndalphardni_dnj_dTau__constdelta_x(HEOS, i, j, xN_flag); +} +CoolPropDbl MixtureDerivatives::d_ndln_fugacity_i_dnj_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag){ + return d_ndalphardni_dDelta(HEOS, j, xN_flag) + d_nd_ndalphardni_dnj_dDelta__consttau_x(HEOS, i, j, xN_flag); +} +CoolPropDbl MixtureDerivatives::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){ + return d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j, xN_flag) + d_nd_ndalphardni_dnj_dxk__consttau_delta(HEOS, i, j, k, xN_flag); +} +CoolPropDbl MixtureDerivatives::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){ + double sum = d_ndln_fugacity_i_dnj_dtau__constdelta_x(HEOS, i, j, xN_flag)*ndtaudni__constT_V_nj(HEOS, k, xN_flag) + + d_ndln_fugacity_i_dnj_ddelta__consttau_x(HEOS, i, j, xN_flag)*nddeltadni__constT_V_nj(HEOS, k, xN_flag) + + d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HEOS, i, j, k, xN_flag); + std::size_t mmax = HEOS.mole_fractions.size(); + if (xN_flag == XN_DEPENDENT){ mmax--; } + for (unsigned int m = 0; m < mmax; ++m) + { + sum -= HEOS.mole_fractions[m]*d_ndln_fugacity_i_dnj_ddxk__consttau_delta(HEOS, i, j, m, xN_flag); + } + return sum; +} + CoolPropDbl MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag) { // Gernert 3.115 @@ -525,10 +550,10 @@ CoolPropDbl MixtureDerivatives::d_nd_ndalphardni_dnj_dDelta__consttau_x(Helmholt CoolPropDbl MixtureDerivatives::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) { - double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag)*d_nddeltadni_dxj__constdelta_tau(HEOS, i, k, xN_flag) + double line1 = d_ndalphardni_dDelta(HEOS, i, xN_flag)*d_nddeltadni_dxj__constdelta_tau(HEOS, j, k, xN_flag) + d2_ndalphardni_dxj_dDelta__consttau_xi(HEOS, i, k, xN_flag)*nddeltadni__constT_V_nj(HEOS, j, xN_flag); - double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag)*d_ndtaudni_dxj__constdelta_tau(HEOS, i, k, xN_flag) + double line2 = d_ndalphardni_dTau(HEOS, i, xN_flag)*d_ndtaudni_dxj__constdelta_tau(HEOS, j, k, xN_flag) + d2_ndalphardni_dxj_dTau__constdelta_xi(HEOS, i, k, xN_flag)*ndtaudni__constT_V_nj(HEOS, j, xN_flag); double line3 = d2_ndalphardni_dxj_dxk__constdelta_tau_xi(HEOS, i, j, k, xN_flag) - d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k, xN_flag); @@ -808,7 +833,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") SECTION(ss00.str(),"") { x_N_dependency_flag xN_flag; - for (std::size_t xNxN = 0; xNxN <=1; xNxN++){ + for (std::size_t xNxN = 1; xNxN > 0; xNxN--){ std::ostringstream ss000; std::string xN_string; if (xNxN == 0){ @@ -1404,7 +1429,6 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") // These derivatives depend on i,j, and k indices for (std::size_t k = 0; k < Ncomp; ++k){ - if (xN_flag == XN_DEPENDENT && k == Ncomp-1){ continue; } HelmholtzEOSMixtureBackend & rHEOS_pluszk = (xN_flag == XN_INDEPENDENT) ? *(HEOS_plusz_xNindep[Ncomp][k].get()) : *(HEOS_plusz_xNdep[Ncomp][k].get()); HelmholtzEOSMixtureBackend & rHEOS_minuszk = (xN_flag == XN_INDEPENDENT) ? *(HEOS_minusz_xNindep[Ncomp][k].get()) : *(HEOS_minusz_xNdep[Ncomp][k].get()); @@ -1414,13 +1438,13 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss1; ss1 << "d3Trdxidxjdxk, i=" << i << ", j=" << j << ", k=" << k; SECTION(ss1.str(), "") { - if (xN_flag == XN_INDEPENDENT){ continue; } - if (j == Ncomp-1){ break; } + if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; } double analytic = rHEOS.Reducing->d3Trdxidxjdxk(rHEOS.get_mole_fractions(), i, j, k, xN_flag); double v1 = rHEOS.Reducing->d2Trdxidxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d2Trdxidxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); double numeric = (v1 - v2)/(2*dz); double err = std::abs((numeric-analytic)/analytic); + if (std::abs(numeric) < DBL_EPSILON && std::abs(analytic) < DBL_EPSILON){ err = 0; } CAPTURE(numeric); CAPTURE(analytic); CHECK(err < 1e-7); @@ -1428,8 +1452,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss2; ss2 << "d3rhormolardxidxjdxk, i=" << i << ", j=" << j << ", k=" << k; SECTION(ss2.str(), "") { - if (xN_flag == XN_INDEPENDENT){ continue; } - if (j == Ncomp-1){ break; } + if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; } double analytic = rHEOS.Reducing->d3rhormolardxidxjdxk(rHEOS.get_mole_fractions(), i, j, k, xN_flag); double v1 = rHEOS.Reducing->d2rhormolardxidxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d2rhormolardxidxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); @@ -1442,8 +1465,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss3; ss3 << "d2_ndTrdni_dxj_dxk__constxi, i=" << i << ", j=" << j << ", k=" << k; SECTION(ss3.str(), "") { - if (xN_flag == XN_INDEPENDENT){ continue; } - if (j == Ncomp-1){ break; } + if ((xN_flag == XN_DEPENDENT) && (j == Ncomp-1 || k == Ncomp-1)){ break; } double analytic = rHEOS.Reducing->d2_ndTrdni_dxj_dxk__constxi(rHEOS.get_mole_fractions(), i, j, k, xN_flag); double v1 = rHEOS.Reducing->d_ndTrdni_dxj__constxi(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d_ndTrdni_dxj__constxi(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); @@ -1456,8 +1478,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss4; ss4 << "d2_ndrhorbardni_dxj_dxk__constxi, i=" << i << ", j=" << j << ", k=" << k; SECTION(ss4.str(), "") { - if (xN_flag == XN_INDEPENDENT){ continue; } - if (j == Ncomp-1){ break; } + if ((xN_flag == XN_DEPENDENT) && (j == Ncomp-1 || k == Ncomp-1)){ break; } double analytic = rHEOS.Reducing->d2_ndrhorbardni_dxj_dxk__constxi(rHEOS.get_mole_fractions(), i, j, k, xN_flag); double v1 = rHEOS.Reducing->d_ndrhorbardni_dxj__constxi(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d_ndrhorbardni_dxj__constxi(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); @@ -1470,7 +1491,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss5; ss5 << "d2_PSI_T_dxj_dxk, i=" << i << ", j=" << j << ", k=" << k; SECTION(ss5.str(), "") { - if (j == Ncomp-1){ break; } + if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; } double analytic = rHEOS.Reducing->d2_PSI_T_dxj_dxk(rHEOS.get_mole_fractions(), i, j, k, xN_flag); double v1 = rHEOS.Reducing->d_PSI_T_dxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d_PSI_T_dxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); @@ -1483,7 +1504,7 @@ TEST_CASE("Mixture derivative checks", "[mixtures],[mixture_derivs]") std::ostringstream ss6; ss6 << "d2_PSI_rho_dxj_dxk, i=" << i << ", j=" << j << ", k=" << k; SECTION(ss6.str(), "") { - if (j == Ncomp-1){ break; } + if ((xN_flag == XN_DEPENDENT) && (i == Ncomp-1 || j == Ncomp-1 || k == Ncomp-1)){ break; } double analytic = rHEOS.Reducing->d2_PSI_rho_dxj_dxk(rHEOS.get_mole_fractions(), i, j, k, xN_flag); double v1 = rHEOS.Reducing->d_PSI_rho_dxj(rHEOS_pluszk.get_mole_fractions(), i, j, xN_flag); double v2 = rHEOS.Reducing->d_PSI_rho_dxj(rHEOS_minuszk.get_mole_fractions(), i, j, xN_flag); diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h index c9c6bf8a..9c0f4c6c 100644 --- a/src/Backends/Helmholtz/MixtureDerivatives.h +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -149,6 +149,42 @@ class MixtureDerivatives{ static CoolPropDbl dln_fugacity_i_dtau__constdelta_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag); static CoolPropDbl dln_fugacity_i_ddelta__consttau_x(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, x_N_dependency_flag xN_flag); static CoolPropDbl dln_fugacity_dxj__constT_rho_xi(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, x_N_dependency_flag xN_flag); + + 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 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 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){ + CoolPropDbl RT = HEOS.gas_constant()*HEOS.T(); + return 1/RT*ndln_fugacity_i_dnj__constT_V_xi(HEOS, i, j, xN_flag); + } + static CoolPropDbl n2Aijk(HelmholtzEOSMixtureBackend &HEOS, std::size_t i, std::size_t j, std::size_t k, x_N_dependency_flag xN_flag){ + 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); + } + static CoolPropDbl L1_star(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){ + Eigen::Matrix2d L1; + CoolPropDbl RT = HEOS.gas_constant()*HEOS.T(); + L1(0, 0) = nAij(HEOS, 0, 0, xN_flag); + L1(1, 0) = nAij(HEOS, 1, 0, xN_flag); + L1(0, 1) = L1(1,0); + L1(1, 1) = nAij(HEOS, 1, 1, xN_flag); + return L1.determinant(); + } + static CoolPropDbl M1_star(HelmholtzEOSMixtureBackend &HEOS, x_N_dependency_flag xN_flag){ + Eigen::Matrix2d M1; + CoolPropDbl RT = HEOS.gas_constant()*HEOS.T(); + M1(0, 0) = nAij(HEOS, 0, 0, xN_flag); + M1(0, 1) = nAij(HEOS, 0, 1, xN_flag); + double dL1_dn0 = nAij(HEOS, 0, 0, xN_flag)*n2Aijk(HEOS, 0, 1, 1, xN_flag) + nAij(HEOS, 1, 1, xN_flag)*n2Aijk(HEOS, 0, 0, 0, xN_flag) - 2*nAij(HEOS, 0, 1, xN_flag)*n2Aijk(HEOS, 0, 0, 1, xN_flag); + double dL1_dn1 = nAij(HEOS, 0, 0, xN_flag)*n2Aijk(HEOS, 1, 1, 1, xN_flag) + nAij(HEOS, 1, 1, xN_flag)*n2Aijk(HEOS, 0, 0, 1, xN_flag) - 2*nAij(HEOS, 0, 1, xN_flag)*n2Aijk(HEOS, 0, 1, 1, xN_flag); + M1(1, 0) = dL1_dn0; + M1(1, 1) = dL1_dn1; + return M1.determinant(); + } /** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions * diff --git a/src/Backends/Helmholtz/ReducingFunctions.cpp b/src/Backends/Helmholtz/ReducingFunctions.cpp index 3ad76fd5..c28d0493 100644 --- a/src/Backends/Helmholtz/ReducingFunctions.cpp +++ b/src/Backends/Helmholtz/ReducingFunctions.cpp @@ -337,19 +337,11 @@ CoolPropDbl GERG2008ReducingFunction::d2Yrdxidxj(const std::vector CoolPropDbl d2Yr_dxidxj = 2*Yc[N-1]; d2Yr_dxidxj += c_Y_ij(i, j, beta, gamma, Y_c_ij)*d2fYijdxidxj(x,i,j,beta); - for (std::size_t k = 0; k < N-1; k++) - { - CoolPropDbl beta_Y_kN = beta[k][N-1]; - d2Yr_dxidxj += 2*c_Y_ij(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); - } - { - CoolPropDbl beta_Y_iN = beta[i][N-1]; - d2Yr_dxidxj += c_Y_ij(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)); - } - { - CoolPropDbl beta_Y_jN = beta[j][N-1]; - d2Yr_dxidxj += -c_Y_ij(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)); + for (std::size_t k = 0; k < N-1; k++){ + d2Yr_dxidxj += c_Y_ij(k, N-1, beta, gamma, Y_c_ij)*d2fYkidxi2__constxk(x, k, N-1, beta); } + d2Yr_dxidxj -= c_Y_ij(i, N-1, beta, gamma, Y_c_ij)*d2fYijdxidxj(x, i, N-1, beta); + d2Yr_dxidxj -= c_Y_ij(j, N-1, beta, gamma, Y_c_ij)*d2fYijdxidxj(x, j, N-1, beta); return d2Yr_dxidxj; } else{ @@ -387,28 +379,44 @@ CoolPropDbl GERG2008ReducingFunction::d3Yrdxidxjdxk(const std::vector