diff --git a/src/Backends/Helmholtz/FlashRoutines.cpp b/src/Backends/Helmholtz/FlashRoutines.cpp index e0c23337..7885f85d 100644 --- a/src/Backends/Helmholtz/FlashRoutines.cpp +++ b/src/Backends/Helmholtz/FlashRoutines.cpp @@ -75,7 +75,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS) options.omega = omega; // Actually call the solver - SaturationSolvers::saturation_T_pure(&HEOS, HEOS._T, options); + SaturationSolvers::saturation_T_pure(HEOS, HEOS._T, options); // If you get here, there was no error, all is well break; @@ -91,7 +91,7 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS) catch(std::exception &){ try{ // We may need to polish the solution at low pressure - SaturationSolvers::saturation_T_pure_1D_P(&HEOS, T, options); + SaturationSolvers::saturation_T_pure_1D_P(HEOS, T, options); } catch(std::exception &){ throw; @@ -144,10 +144,10 @@ void FlashRoutines::QT_flash(HelmholtzEOSMixtureBackend &HEOS) options.Nstep_max = 20; // Get an extremely rough guess by interpolation of ln(p) v. T curve where the limits are mole-fraction-weighted - long double pguess = SaturationSolvers::saturation_preconditioner(&HEOS, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions); + long double pguess = SaturationSolvers::saturation_preconditioner(HEOS, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions); // Use Wilson iteration to obtain updated guess for pressure - pguess = SaturationSolvers::saturation_Wilson(&HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess); + pguess = SaturationSolvers::saturation_Wilson(HEOS, HEOS._Q, HEOS._T, SaturationSolvers::imposed_T, HEOS.mole_fractions, pguess); // Actually call the successive substitution solver SaturationSolvers::successive_substitution(HEOS, HEOS._Q, HEOS._T, pguess, HEOS.mole_fractions, HEOS.K, options); @@ -210,7 +210,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS) options.omega = omega; // Actually call the solver - SaturationSolvers::saturation_PHSU_pure(&HEOS, HEOS._p, options); + SaturationSolvers::saturation_PHSU_pure(HEOS, HEOS._p, options); // If you get here, there was no error, all is well break; @@ -225,7 +225,7 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS) } catch(std::exception &){ // We may need to polish the solution at low pressure - SaturationSolvers::saturation_P_pure_1D_T(&HEOS, HEOS._p, options); + SaturationSolvers::saturation_P_pure_1D_T(HEOS, HEOS._p, options); } // Load the outputs @@ -250,10 +250,10 @@ void FlashRoutines::PQ_flash(HelmholtzEOSMixtureBackend &HEOS) io.Nstep_max = 20; // 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); + long double Tguess = SaturationSolvers::saturation_preconditioner(HEOS, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions); // Use Wilson iteration to obtain updated guess for temperature - Tguess = SaturationSolvers::saturation_Wilson(&HEOS, HEOS._Q, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions, Tguess); + Tguess = SaturationSolvers::saturation_Wilson(HEOS, HEOS._Q, HEOS._p, SaturationSolvers::imposed_p, HEOS.mole_fractions, Tguess); // Actually call the successive substitution solver SaturationSolvers::successive_substitution(HEOS, HEOS._Q, Tguess, HEOS._p, HEOS.mole_fractions, HEOS.K, io); @@ -350,14 +350,14 @@ void FlashRoutines::PHSU_D_flash(HelmholtzEOSMixtureBackend &HEOS, int other) if (HEOS._rhomolar > HEOS._crit.rhomolar) { options.imposed_rho = SaturationSolvers::saturation_D_pure_options::IMPOSED_RHOL; - SaturationSolvers::saturation_D_pure(&HEOS, HEOS._rhomolar, options); + SaturationSolvers::saturation_D_pure(HEOS, HEOS._rhomolar, options); // SatL and SatV have the saturation values Sat = HEOS.SatL; } else { options.imposed_rho = SaturationSolvers::saturation_D_pure_options::IMPOSED_RHOV; - SaturationSolvers::saturation_D_pure(&HEOS, HEOS._rhomolar, options); + SaturationSolvers::saturation_D_pure(HEOS, HEOS._rhomolar, options); // SatL and SatV have the saturation values Sat = HEOS.SatV; } diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 49da191e..aa22829a 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -26,6 +26,7 @@ #include "VLERoutines.h" #include "FlashRoutines.h" #include "TransportRoutines.h" +#include "MixtureDerivatives.h" static int deriv_counter = 0; @@ -1132,7 +1133,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot // Run the saturation routines to determine the saturation densities and pressures HelmholtzEOSMixtureBackend HEOS(components); SaturationSolvers::saturation_T_pure_options options; - SaturationSolvers::saturation_T_pure(&HEOS, _T, options); + SaturationSolvers::saturation_T_pure(HEOS, _T, options); long double Q; @@ -1991,7 +1992,7 @@ long double HelmholtzEOSMixtureBackend::calc_gibbsmolar(void) } long double HelmholtzEOSMixtureBackend::calc_fugacity_coefficient(int i) { - return exp(mixderiv_ln_fugacity_coefficient(i)); + return exp(MixtureDerivatives::ln_fugacity_coefficient(*this, i)); } SimpleState HelmholtzEOSMixtureBackend::calc_reducing_state_nocache(const std::vector & mole_fractions) @@ -2291,228 +2292,4 @@ long double HelmholtzEOSMixtureBackend::calc_d3alpha0_dTau3(void) } -long double HelmholtzEOSMixtureBackend::mixderiv_dalphar_dxi(int i) -{ - return components[i]->pEOS->baser(_tau, _delta) + Excess.dalphar_dxi(_tau, _delta, mole_fractions, i); -} -long double HelmholtzEOSMixtureBackend::mixderiv_d2alphar_dxi_dTau(int i) -{ - return components[i]->pEOS->dalphar_dTau(_tau, _delta) + Excess.d2alphar_dxi_dTau(_tau, _delta, mole_fractions, i); -} -long double HelmholtzEOSMixtureBackend::mixderiv_d2alphar_dxi_dDelta(int i) -{ - return components[i]->pEOS->dalphar_dDelta(_tau, _delta) + Excess.d2alphar_dxi_dDelta(_tau, _delta, mole_fractions, i); -} -long double HelmholtzEOSMixtureBackend::mixderiv_d2alphardxidxj(int i, int j) -{ - return 0 + Excess.d2alphardxidxj(_tau, _delta, mole_fractions, i, j); -} - -long double HelmholtzEOSMixtureBackend::mixderiv_ln_fugacity_coefficient(int i) -{ - return alphar() + mixderiv_ndalphar_dni__constT_V_nj(i)-log(1+_delta.pt()*dalphar_dDelta()); -} -long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_dT__constrho_n(int i) -{ - double dtau_dT = -_tau.pt()/_T; //[1/K] - return (dalphar_dTau() + mixderiv_d_ndalphardni_dTau(i)-1/(1+_delta.pt()*dalphar_dDelta())*(_delta.pt()*d2alphar_dDelta_dTau()))*dtau_dT; -} -long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_drho__constT_n(int i) -{ - double ddelta_drho = 1/_reducing.rhomolar; //[m^3/mol] - return (dalphar_dDelta() + mixderiv_d_ndalphardni_dDelta(i)-1/(1+_delta.pt()*dalphar_dDelta())*(_delta.pt()*d2alphar_dDelta2()+dalphar_dDelta()))*ddelta_drho; -} -long double HelmholtzEOSMixtureBackend::mixderiv_dnalphar_dni__constT_V_nj(int i) -{ - // GERG Equation 7.42 - return alphar() + mixderiv_ndalphar_dni__constT_V_nj(i); -} -long double HelmholtzEOSMixtureBackend::mixderiv_d2nalphar_dni_dT(int i) -{ - return -_tau.pt()/_T*(dalphar_dTau() + mixderiv_d_ndalphardni_dTau(i)); -} -long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_dT__constp_n(int i) -{ - double T = _reducing.T/_tau.pt(); - long double R_u = static_cast(_gas_constant); - return mixderiv_d2nalphar_dni_dT(i) + 1/T-mixderiv_partial_molar_volume(i)/(R_u*T)*mixderiv_dpdT__constV_n(); -} -long double HelmholtzEOSMixtureBackend::mixderiv_partial_molar_volume(int i) -{ - return -mixderiv_ndpdni__constT_V_nj(i)/mixderiv_ndpdV__constT_n(); -} - -long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_dp__constT_n(int i) -{ - // GERG equation 7.30 - long double R_u = static_cast(_gas_constant); - double partial_molar_volume = mixderiv_partial_molar_volume(i); // [m^3/mol] - double term1 = partial_molar_volume/(R_u*_T); // m^3/mol/(N*m)*mol = m^2/N = 1/Pa - double term2 = 1.0/p(); - return term1 - term2; -} - -long double HelmholtzEOSMixtureBackend::mixderiv_dln_fugacity_coefficient_dxj__constT_p_xi(int i, int j) -{ - // Gernert 3.115 - long double R_u = static_cast(_gas_constant); - // partial molar volume is -dpdn/dpdV, so need to flip the sign here - return mixderiv_d2nalphar_dni_dxj__constT_V(i,j) - mixderiv_partial_molar_volume(i)/(R_u*_T)*mixderiv_dpdxj__constT_V_xi(j); -} -long double HelmholtzEOSMixtureBackend::mixderiv_dpdxj__constT_V_xi(int j) -{ - // Gernert 3.130 - long double R_u = static_cast(_gas_constant); - return _rhomolar*R_u*_T*(mixderiv_ddelta_dxj__constT_V_xi(j)*dalphar_dDelta()+_delta.pt()*mixderiv_d_dalpharddelta_dxj__constT_V_xi(j)); -} - -long double HelmholtzEOSMixtureBackend::mixderiv_d_dalpharddelta_dxj__constT_V_xi(int j) -{ - // Gernert Equation 3.134 (Catch test provided) - return d2alphar_dDelta2()*mixderiv_ddelta_dxj__constT_V_xi(j) - + d2alphar_dDelta_dTau()*mixderiv_dtau_dxj__constT_V_xi(j) - + mixderiv_d2alphar_dxi_dDelta(j); -} - -long double HelmholtzEOSMixtureBackend::mixderiv_dalphar_dxj__constT_V_xi(int j) -{ - //Gernert 3.119 (Catch test provided) - return dalphar_dDelta()*mixderiv_ddelta_dxj__constT_V_xi(j)+dalphar_dTau()*mixderiv_dtau_dxj__constT_V_xi(j)+mixderiv_dalphar_dxi(j); -} -long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dxj__constT_V_xi(int i, int j) -{ - // Gernert 3.118 - return mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(i,j) - + mixderiv_ddelta_dxj__constT_V_xi(j)*mixderiv_d_ndalphardni_dDelta(i) - + mixderiv_dtau_dxj__constT_V_xi(j)*mixderiv_d_ndalphardni_dTau(i); -} -long double HelmholtzEOSMixtureBackend::mixderiv_ddelta_dxj__constT_V_xi(int j) -{ - // Gernert 3.121 (Catch test provided) - return -_delta.pt()/_reducing.rhomolar*Reducing.p->drhormolardxi__constxj(mole_fractions,j); -} -long double HelmholtzEOSMixtureBackend::mixderiv_dtau_dxj__constT_V_xi(int j) -{ - // Gernert 3.122 (Catch test provided) - return 1/_T*Reducing.p->dTrdxi__constxj(mole_fractions,j); -} - -long double HelmholtzEOSMixtureBackend::mixderiv_dpdT__constV_n() -{ - long double R_u = static_cast(_gas_constant); - return _rhomolar*R_u*(1+_delta.pt()*dalphar_dDelta()-_delta.pt()*_tau.pt()*d2alphar_dDelta_dTau()); -} -long double HelmholtzEOSMixtureBackend::mixderiv_dpdrho__constT_n() -{ - long double R_u = static_cast(_gas_constant); - return R_u*_T*(1+2*_delta.pt()*dalphar_dDelta()+pow(_delta.pt(),2)*d2alphar_dDelta2()); -} -long double HelmholtzEOSMixtureBackend::mixderiv_ndpdV__constT_n() -{ - long double R_u = static_cast(_gas_constant); - return -pow(_rhomolar,2)*R_u*_T*(1+2*_delta.pt()*dalphar_dDelta()+pow(_delta.pt(),2)*d2alphar_dDelta2()); -} -long double HelmholtzEOSMixtureBackend::mixderiv_ndpdni__constT_V_nj(int i) -{ - // Eqn 7.64 and 7.63 - long double R_u = static_cast(_gas_constant); - double ndrhorbar_dni__constnj = Reducing.p->ndrhorbardni__constnj(mole_fractions,i); - double ndTr_dni__constnj = Reducing.p->ndTrdni__constnj(mole_fractions,i); - double summer = 0; - for (unsigned int k = 0; k < mole_fractions.size(); ++k) - { - summer += mole_fractions[k]*mixderiv_d2alphar_dxi_dDelta(k); - } - double nd2alphar_dni_dDelta = _delta.pt()*d2alphar_dDelta2()*(1-1/_reducing.rhomolar*ndrhorbar_dni__constnj)+_tau.pt()*d2alphar_dDelta_dTau()/_reducing.T*ndTr_dni__constnj+mixderiv_d2alphar_dxi_dDelta(i)-summer; - return _rhomolar*R_u*_T*(1+_delta.pt()*dalphar_dDelta()*(2-1/_reducing.rhomolar*ndrhorbar_dni__constnj)+_delta.pt()*nd2alphar_dni_dDelta); -} - -long double HelmholtzEOSMixtureBackend::mixderiv_ndalphar_dni__constT_V_nj(int i) -{ - double term1 = _delta.pt()*dalphar_dDelta()*(1-1/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions,i)); - double term2 = _tau.pt()*dalphar_dTau()*(1/_reducing.T)*Reducing.p->ndTrdni__constnj(mole_fractions,i); - - double s = 0; - for (unsigned int k = 0; k < mole_fractions.size(); k++) - { - s += mole_fractions[k]*mixderiv_dalphar_dxi(k); - } - double term3 = mixderiv_dalphar_dxi(i); - return term1 + term2 + term3 - s; -} -long double HelmholtzEOSMixtureBackend::mixderiv_ndln_fugacity_coefficient_dnj__constT_p(int i, int j) -{ - long double R_u = static_cast(_gas_constant); - return mixderiv_nd2nalphardnidnj__constT_V(j, i) + 1 - mixderiv_partial_molar_volume(j)/(R_u*_T)*mixderiv_ndpdni__constT_V_nj(i); -} -long double HelmholtzEOSMixtureBackend::mixderiv_nddeltadni__constT_V_nj(int i) -{ - return _delta.pt()-_delta.pt()/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions, i); -} -long double HelmholtzEOSMixtureBackend::mixderiv_ndtaudni__constT_V_nj(int i) -{ - return _tau.pt()/_reducing.T*Reducing.p->ndTrdni__constnj(mole_fractions, i); -} -long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(int i, int j) -{ - double line1 = _delta.pt()*mixderiv_d2alphar_dxi_dDelta(j)*(1-1/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions, i)); - double line2 = -_delta.pt()*dalphar_dDelta()*(1/_reducing.rhomolar)*(Reducing.p->d_ndrhorbardni_dxj__constxi(mole_fractions, i, j)-1/_reducing.rhomolar*Reducing.p->drhormolardxi__constxj(mole_fractions,j)*Reducing.p->ndrhorbardni__constnj(mole_fractions,i)); - double line3 = _tau.pt()*mixderiv_d2alphar_dxi_dTau(j)*(1/_reducing.T)*Reducing.p->ndTrdni__constnj(mole_fractions, i); - double line4 = _tau.pt()*dalphar_dTau()*(1/_reducing.T)*(Reducing.p->d_ndTrdni_dxj__constxi(mole_fractions,i,j)-1/_reducing.T*Reducing.p->dTrdxi__constxj(mole_fractions,j)*Reducing.p->ndTrdni__constnj(mole_fractions, i)); - double s = 0; - for (unsigned int m = 0; m < mole_fractions.size(); m++) - { - s += mole_fractions[m]*mixderiv_d2alphardxidxj(j,m); - } - double line5 = mixderiv_d2alphardxidxj(i,j)-mixderiv_dalphar_dxi(j)-s; - return line1+line2+line3+line4+line5; -} -long double HelmholtzEOSMixtureBackend::mixderiv_nd2nalphardnidnj__constT_V(int i, int j) -{ - double line0 = mixderiv_ndalphar_dni__constT_V_nj(j); // First term from 7.46 - double line1 = mixderiv_d_ndalphardni_dDelta(i)*mixderiv_nddeltadni__constT_V_nj(j); - double line2 = mixderiv_d_ndalphardni_dTau(i)*mixderiv_ndtaudni__constT_V_nj(j); - double summer = 0; - for (unsigned int k = 0; k < mole_fractions.size(); k++) - { - summer += mole_fractions[k]*mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(i, k); - } - double line3 = mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(i, j)-summer; - return line0 + line1 + line2 + line3; -} -long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dDelta(int i) -{ - // The first line - double term1 = (_delta.pt()*d2alphar_dDelta2()+dalphar_dDelta())*(1-1/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions, i)); - - // The second line - double term2 = _tau.pt()*d2alphar_dDelta_dTau()*(1/_reducing.T)*Reducing.p->ndTrdni__constnj(mole_fractions, i); - - // The third line - double term3 = mixderiv_d2alphar_dxi_dDelta(i); - for (unsigned int k = 0; k < mole_fractions.size(); k++) - { - term3 -= mole_fractions[k]*mixderiv_d2alphar_dxi_dDelta(k); - } - return term1 + term2 + term3; -} - -long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dTau(int i) -{ - // The first line - double term1 = _delta.pt()*d2alphar_dDelta_dTau()*(1-1/_reducing.rhomolar*Reducing.p->ndrhorbardni__constnj(mole_fractions, i)); - - // The second line - double term2 = (_tau.pt()*d2alphar_dTau2()+dalphar_dTau())*(1/_reducing.T)*Reducing.p->ndTrdni__constnj(mole_fractions, i); - - // The third line - double term3 = mixderiv_d2alphar_dxi_dTau(i); - for (unsigned int k = 0; k < mole_fractions.size(); k++) - { - term3 -= mole_fractions[k]*mixderiv_d2alphar_dxi_dTau(k); - } - return term1 + term2 + term3; -} - - } /* namespace CoolProp */ diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index f1133ac0..e404cf92 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -41,9 +41,7 @@ public: friend class FlashRoutines; // Allows the static methods in the FlashRoutines class to have access to all the protected members and methods of this class friend class TransportRoutines; // Allows the static methods in the TransportRoutines class to have access to all the protected members and methods of this class - - - //friend class MixtureDerivatives; // + friend class MixtureDerivatives; // Allows the static methods in the MixtureDerivatives class to have access to all the protected members and methods of this class // Helmholtz EOS backend uses mole fractions bool using_mole_fractions(){return true;} @@ -241,246 +239,6 @@ public: long double solver_rho_Tp_SRK(long double T, long double p, int phase); long double solver_for_rho_given_T_oneof_HSU(long double T, long double value, int other); - - - - // *************************************************************** - // *************************************************************** - // ***************** MIXTURE DERIVATIVES *********************** - // *************************************************************** - // *************************************************************** - - - long double mixderiv_dalphar_dxi(int i); - long double mixderiv_d2alphar_dxi_dTau(int i); - long double mixderiv_d2alphar_dxi_dDelta(int i); - long double mixderiv_d2alphardxidxj(int i, int j); - - /*! \brief GERG 2004 Monograph equation 7.61 - * - * The derivative term - \f[ - \left(\frac{\partial p}{\partial T} \right)_{V,\bar n} = \rho R(1+\delta \alpha_{\delta}^r-\delta \tau \alpha^r_{\delta\tau}) - \f] - - */ - long double mixderiv_dpdT__constV_n(); - - long double mixderiv_dpdrho__constT_n(); - - /** \brief GERG 2004 Monograph equation 7.62 - * - * The derivative term - \f[ - n\left(\frac{\partial p}{\partial V} \right)_{T,\bar n} = -\rho^2 RT(1+2\delta \alpha_{\delta}^r+\delta^2\alpha^r_{\delta\delta}) - \f] - */ - long double mixderiv_ndpdV__constT_n(); - - /** \brief GERG 2004 Monograph equation 7.63 - * - * The derivative term - \f[ - n\left(\frac{\partial p}{\partial n_i} \right)_{T,V,n_j} = \rho RT\left[1+\delta\alpha_{\delta}^r\left[2- \frac{1}{\rho_r}\cdot n\left( \frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] +\delta\cdot n\left(\frac{\partial\alpha_{\delta}^r}{\partial n_i}\right)_{T,V,n_j}\right] - \f] - - */ - long double mixderiv_ndpdni__constT_V_nj(int i); - - /** GERG 2004 monograph Eqn. 7.32 - - * The partial molar volume - \f[ - \hat v_i = \left( \frac{\partial V}{\partial n_i}\right)_{T,p,n_j} = \frac{-\left(\dfrac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\dfrac{\partial p}{\partial V}\right)_{T,\bar n}} - \f] - */ - long double mixderiv_partial_molar_volume(int i); - - /*! - Natural logarithm of the fugacity coefficient - */ - long double mixderiv_ln_fugacity_coefficient(int i); - - /*! - Derivative of the natural logarithm of the fugacity coefficient with respect to T - */ - long double mixderiv_dln_fugacity_coefficient_dT__constrho_n(int i); - - /*! - Derivative of the natural logarithm of the fugacity coefficient with respect to T - */ - long double mixderiv_dln_fugacity_coefficient_drho__constT_n(int i); - - /// GERG 2004 Monograph Eqn. 7.29 - /** The derivative term - \f[ - \left(\frac{\partial \ln \phi_i}{\partial T} \right)_{p,\bar n} = \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} + \frac{1}{T}-\frac{\hat v}{RT}\left(\frac{\partial p}{\partial T}\right)_{V,\bar n} - \f] - */ - long double mixderiv_dln_fugacity_coefficient_dT__constp_n(int i); - - /// Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions - /*! The derivative term - \f[ - n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} - \f] - which is equal to - \f{eqnarray*}{ - n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} &=& \delta \phi^r_{\delta}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\ - && +\tau \phi^r_{\tau}\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\ - && +\phi^r_{x_i}-\sum_{k=1}^{N}x_k\phi^r_{x_k} - \f} - */ - long double mixderiv_ndalphar_dni__constT_V_nj(int i); - - /// GERG Equation 7.42 - long double mixderiv_dnalphar_dni__constT_V_nj(int i); - - /// GERG 2004 Monograph Eqn. 7.30 - /** - The derivative term - \f[ - \left(\frac{\partial \ln \phi_i}{\partial p} \right)_{T,\bar n} = \frac{\hat v_i}{RT}-\frac{1}{p} - \f] - */ - long double mixderiv_dln_fugacity_coefficient_dp__constT_n(int i); - - ///GERG 2004 Monograph Equation 7.31 - /** The derivative term - \f[ - n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1+\frac{n}{RT}\frac{\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}} - \f] - which is also equal to - \f[ - n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1-\frac{\hat v_i}{RT}\left[n\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\right] - \f] - */ - long double mixderiv_ndln_fugacity_coefficient_dnj__constT_p(int i, int j); - - /// Gernert Equation 3.115 - /** The derivative term - \f[ - \left(\frac{\partial \ln \phi_i}{\partial x_j}\right)_{T,p,x_{k\neq j}} = \left(\frac{\partial^2n\alpha^r}{\partial x_j \partial n_i} \right)_{T,V}+\frac{1}{RT}\frac{\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_{k\neq i}}\left(\frac{\partial p}{\partial x_j}\right)_{T,V,x_{k\neq j}}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}} - \f] - */ - long double mixderiv_dln_fugacity_coefficient_dxj__constT_p_xi(int i, int j); - - /// Gernert Equation 3.130 - /** The derivative term - \f[ - \left(\frac{\partial p}{\partial x_j} \right)_{T,V,x_{k\neq j}} = \rho RT\left(-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_{k\neq j}} \delta\alpha_{\delta}^r + \delta\left(\frac{\partial}{\partial x_j}\left(\left( \frac{\partial \alpha^r}{\partial \delta}\right)_{\tau,\bar x}\right)\right)_{T,V,x_{k\neq j}}\right) - \f] - */ - long double mixderiv_dpdxj__constT_V_xi(int j); - - /// Gernert Equation 3.117 - /** The derivative term - \f[ - \left(\frac{\partial^2n\alpha^r}{\partial x_i\partial n_j} \right)_{T,V} = \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_{j\neq i}}\right)\right)_{T,V,x_{k\neq j}} +\left(\frac{\partial \alpha^r}{\partial x_j}\right)_{T,V,x_{k\neq j}} - \f] - */ - long double mixderiv_d2nalphar_dxi_dxj__constT_V(int i, int j){ return mixderiv_d_ndalphardni_dxj__constT_V_xi(i, j) + mixderiv_dalphar_dxj__constT_V_xi(j);}; - - /// Gernert Equation 3.119 - /// Catch test provided - long double mixderiv_dalphar_dxj__constT_V_xi(int j); - - /// Gernert Equation 3.118 - /// Catch test provided - long double mixderiv_d_ndalphardni_dxj__constT_V_xi(int i, int j); - - /// Gernert Equation 3.134 - /// Catch test provided - long double mixderiv_d_dalpharddelta_dxj__constT_V_xi(int j); - - /// Gernert Equation 3.121 - /// Catch test provided - long double mixderiv_ddelta_dxj__constT_V_xi(int j); - - /// Gernert Equation 3.122 - /// Catch test provided - long double mixderiv_dtau_dxj__constT_V_xi(int j); - - /// GERG 2004 Monograph, equations 7.44 and 7.51 - /** The derivative term - \f[ - \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = \left( \frac{\partial}{\partial T}\left(\frac{\partial n \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right)_{V,\bar n} - \f] - \f[ - \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = -\frac{\tau}{T}\left[\alpha_{\tau}^r +\left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\bar x}\right] - \f] - */ - long double mixderiv_d2nalphar_dni_dT(int i); - - /// GERG 2004 Monograph Equation 7.51 and Table B4, Kunz, JCED, 2012 - /** The derivative term - \f{eqnarray*}{ - \frac{\partial }{\partial \tau} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right) &=& \delta \phi^r_{\delta\tau}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\ - && +(\tau \phi^r_{\tau\tau}+\phi^r_{\tau})\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\ - && +\phi^r_{x_i\tau}-\sum_{k=1}^{N}x_k\phi^r_{x_k\tau} - \f} - */ - long double mixderiv_d_ndalphardni_dTau(int i); - - /** \brief GERG 2004 Monograph Equation 7.50 and Table B4, Kunz, JCED, 2012 - * - The derivative term - \f{eqnarray*}{ - \left(\frac{\partial }{\partial \delta} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right)\right)_{\tau,\bar x} &=& (\alpha_{\delta}^r+\delta\alpha_{\delta\delta}^r)\left[1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j} \right] \\ - &+&\tau\alpha^r_{\delta\tau}\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\ - &+&\phi^r_{\delta x_i}-\sum_{k=1}^{N}x_k\phi^r_{\delta x_k} - \f} - */ - long double mixderiv_d_ndalphardni_dDelta(int i); - - /** \brief GERG 2004 Monograph equation 7.41 - * The derivative term - \f[ - n\left(\frac{\partial^2n\alpha^r}{\partial n_i \partial n_j} \right)_{T,V} = n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i} - \f] - and - GERG 2004 Monograph equation 7.46: - \f[ - n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i} = n\left( \frac{\partial \alpha^r}{\partial n_j}\right)_{T,V,n_i} + n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i} - \f] - GERG 2004 Monograph equation 7.47: - \f{eqnarray*}{ - n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i} &=& \left( \frac{\partial}{\partial \delta}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\delta}{\partial n_j}\right)_{T,V,n_i}\\ - &+& \left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\tau}{\partial n_j}\right)_{T,V,n_i}\\ - &+& \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}-\sum_{k=1}^{N}x_k \left( \frac{\partial}{\partial x_k}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}\\ - \f} - */ - long double mixderiv_nd2nalphardnidnj__constT_V(int i, int j); - - /** \brief GERG 2004 Monograph equation 7.48 - The derivative term - \f[ - n\left(\frac{\partial \delta}{\partial n_i} \right)_{T,V,n_j} = \delta - \frac{\delta}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j} - \f] - */ - long double mixderiv_nddeltadni__constT_V_nj(int i); - - /** \brief GERG 2004 Monograph equation 7.49 - * - The derivative term - \f[ - n\left(\frac{\partial \tau}{\partial n_i} \right)_{T,V,n_j} = \frac{\tau}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j} - \f] - */ - long double mixderiv_ndtaudni__constT_V_nj(int i); - - /** \brief GERG 2004 Monograph equation 7.52 - * - The derivative term - \f{eqnarray*}{ - \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i} &=& \delta\alpha_{\delta x_j}^{r}\left[ 1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\ - &-& \delta\alpha_{\delta}^{r}\frac{1}{\rho_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\ - &+& \tau\alpha_{\tau x_j}^r\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\ - &+& \tau\alpha_{\tau}^{r}\frac{1}{T_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{T_r}\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right] \\ - &+& \alpha_{x_ix_j}^r-\alpha_{x_j}^r-\sum_{m=1}^Nx_m\alpha_{x_jx_m}^r - \f} - */ - long double mixderiv_d_ndalphardni_dxj__constdelta_tau_xi(int i, int j); }; } /* namespace CoolProp */ diff --git a/src/Backends/Helmholtz/MixtureDerivatives.cpp b/src/Backends/Helmholtz/MixtureDerivatives.cpp new file mode 100644 index 00000000..2a71bef2 --- /dev/null +++ b/src/Backends/Helmholtz/MixtureDerivatives.cpp @@ -0,0 +1,229 @@ +#include "MixtureDerivatives.h" + +namespace CoolProp{ + +long double MixtureDerivatives::dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + return HEOS.components[i]->pEOS->baser(HEOS._tau, HEOS._delta) + HEOS.Excess.dalphar_dxi(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i); +} +long double MixtureDerivatives::d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + return HEOS.components[i]->pEOS->dalphar_dTau(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dTau(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i); +} +long double MixtureDerivatives::d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + return HEOS.components[i]->pEOS->dalphar_dDelta(HEOS._tau, HEOS._delta) + HEOS.Excess.d2alphar_dxi_dDelta(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i); +} +long double MixtureDerivatives::d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, int i, int j) +{ + return 0 + HEOS.Excess.d2alphardxidxj(HEOS._tau, HEOS._delta, HEOS.mole_fractions, i, j); +} + +long double MixtureDerivatives::ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i)-log(1+HEOS._delta.pt()*HEOS.dalphar_dDelta()); +} +long double MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + double dtau_dT = -HEOS._tau.pt()/HEOS._T; //[1/K] + return (HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i)-1/(1+HEOS._delta.pt()*HEOS.dalphar_dDelta())*(HEOS._delta.pt()*HEOS.d2alphar_dDelta_dTau()))*dtau_dT; +} +long double MixtureDerivatives::dln_fugacity_coefficient_drho__constT_n(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + double ddelta_drho = 1/HEOS._reducing.rhomolar; //[m^3/mol] + return (HEOS.dalphar_dDelta() + d_ndalphardni_dDelta(HEOS, i)-1/(1+HEOS._delta.pt()*HEOS.dalphar_dDelta())*(HEOS._delta.pt()*HEOS.d2alphar_dDelta2()+HEOS.dalphar_dDelta()))*ddelta_drho; +} +long double MixtureDerivatives::dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + // GERG Equation 7.42 + return HEOS.alphar() + ndalphar_dni__constT_V_nj(HEOS, i); +} +long double MixtureDerivatives::d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + return -HEOS._tau.pt()/HEOS._T*(HEOS.dalphar_dTau() + d_ndalphardni_dTau(HEOS, i)); +} +long double MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + double T = HEOS._reducing.T/HEOS._tau.pt(); + long double R_u = HEOS.gas_constant(); + return d2nalphar_dni_dT(HEOS, i) + 1/T-partial_molar_volume(HEOS, i)/(R_u*T)*dpdT__constV_n(HEOS); +} +long double MixtureDerivatives::partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + return -ndpdni__constT_V_nj(HEOS, i)/ndpdV__constT_n(HEOS); +} + +long double MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + // GERG equation 7.30 + long double R_u = HEOS.gas_constant(); + double partial_molar_volumeval = partial_molar_volume(HEOS, i); // [m^3/mol] + double term1 = partial_molar_volumeval/(R_u*HEOS._T); // m^3/mol/(N*m)*mol = m^2/N = 1/Pa + double term2 = 1.0/HEOS.p(); + return term1 - term2; +} + +long double MixtureDerivatives::dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j) +{ + // Gernert 3.115 + long double R_u = HEOS.gas_constant(); + // partial molar volume is -dpdn/dpdV, so need to flip the sign here + return d2nalphar_dxi_dnj__constT_V(HEOS, i, j) - partial_molar_volume(HEOS, i)/(R_u*HEOS._T)*dpdxj__constT_V_xi(HEOS, j); +} +long double MixtureDerivatives::dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j) +{ + // Gernert 3.130 + long double R_u = HEOS.gas_constant(); + return HEOS._rhomolar*R_u*HEOS._T*(ddelta_dxj__constT_V_xi(HEOS, j)*HEOS.dalphar_dDelta()+HEOS._delta.pt()*d_dalpharddelta_dxj__constT_V_xi(HEOS, j)); +} + +long double MixtureDerivatives::d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j) +{ + // Gernert Equation 3.134 (Catch test provided) + return HEOS.d2alphar_dDelta2()*ddelta_dxj__constT_V_xi(HEOS, j) + + HEOS.d2alphar_dDelta_dTau()*dtau_dxj__constT_V_xi(HEOS, j) + + d2alphar_dxi_dDelta(HEOS, j); +} + +long double MixtureDerivatives::dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j) +{ + //Gernert 3.119 (Catch test provided) + return HEOS.dalphar_dDelta()*ddelta_dxj__constT_V_xi(HEOS, j)+HEOS.dalphar_dTau()*dtau_dxj__constT_V_xi(HEOS, j)+dalphar_dxi(HEOS, j); +} +long double MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j) +{ + // Gernert 3.118 + return d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i,j) + + ddelta_dxj__constT_V_xi(HEOS, j)*d_ndalphardni_dDelta(HEOS, i) + + dtau_dxj__constT_V_xi(HEOS, j)*d_ndalphardni_dTau(HEOS, i); +} +long double MixtureDerivatives::ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j) +{ + // Gernert 3.121 (Catch test provided) + return -HEOS._delta.pt()/HEOS._reducing.rhomolar*HEOS.Reducing.p->drhormolardxi__constxj(HEOS.mole_fractions,j); +} +long double MixtureDerivatives::dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j) +{ + // Gernert 3.122 (Catch test provided) + return 1/HEOS._T*HEOS.Reducing.p->dTrdxi__constxj(HEOS.mole_fractions,j); +} + +long double MixtureDerivatives::dpdT__constV_n(HelmholtzEOSMixtureBackend &HEOS) +{ + long double R_u = HEOS.gas_constant(); + return HEOS._rhomolar*R_u*(1+HEOS._delta.pt()*HEOS.dalphar_dDelta()-HEOS._delta.pt()*HEOS._tau.pt()*HEOS.d2alphar_dDelta_dTau()); +} +long double MixtureDerivatives::dpdrho__constT_n(HelmholtzEOSMixtureBackend &HEOS) +{ + long double R_u = HEOS.gas_constant(); + return R_u*HEOS._T*(1+2*HEOS._delta.pt()*HEOS.dalphar_dDelta()+pow(HEOS._delta.pt(),2)*HEOS.d2alphar_dDelta2()); +} +long double MixtureDerivatives::ndpdV__constT_n(HelmholtzEOSMixtureBackend &HEOS) +{ + long double R_u = HEOS.gas_constant(); + return -pow(HEOS._rhomolar,2)*R_u*HEOS._T*(1+2*HEOS._delta.pt()*HEOS.dalphar_dDelta()+pow(HEOS._delta.pt(),2)*HEOS.d2alphar_dDelta2()); +} +long double MixtureDerivatives::ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + // Eqn 7.64 and 7.63 + long double R_u = HEOS.gas_constant(); + double ndrhorbar_dni__constnj = HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i); + double ndTr_dni__constnj = HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions,i); + double summer = 0; + for (unsigned int k = 0; k < HEOS.mole_fractions.size(); ++k) + { + summer += HEOS.mole_fractions[k]*d2alphar_dxi_dDelta(HEOS, k); + } + double nd2alphar_dni_dDelta = HEOS._delta.pt()*HEOS.d2alphar_dDelta2()*(1-1/HEOS._reducing.rhomolar*ndrhorbar_dni__constnj)+HEOS._tau.pt()*HEOS.d2alphar_dDelta_dTau()/HEOS._reducing.T*ndTr_dni__constnj+d2alphar_dxi_dDelta(HEOS, i)-summer; + return HEOS._rhomolar*R_u*HEOS._T*(1+HEOS._delta.pt()*HEOS.dalphar_dDelta()*(2-1/HEOS._reducing.rhomolar*ndrhorbar_dni__constnj)+HEOS._delta.pt()*nd2alphar_dni_dDelta); +} + +long double MixtureDerivatives::ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + double term1 = HEOS._delta.pt()*HEOS.dalphar_dDelta()*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i)); + double term2 = HEOS._tau.pt()*HEOS.dalphar_dTau()*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions,i); + + double s = 0; + for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++) + { + s += HEOS.mole_fractions[k]*dalphar_dxi(HEOS, k); + } + double term3 = dalphar_dxi(HEOS, i); + return term1 + term2 + term3 - s; +} +long double MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, int i, int j) +{ + long double R_u = HEOS.gas_constant(); + return nd2nalphardnidnj__constT_V(HEOS, j, i) + 1 - partial_molar_volume(HEOS, j)/(R_u*HEOS._T)*ndpdni__constT_V_nj(HEOS, i); +} +long double MixtureDerivatives::nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + return HEOS._delta.pt()-HEOS._delta.pt()/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i); +} +long double MixtureDerivatives::ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + return HEOS._tau.pt()/HEOS._reducing.T*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i); +} +long double MixtureDerivatives::d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j) +{ + double line1 = HEOS._delta.pt()*d2alphar_dxi_dDelta(HEOS, j)*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i)); + double line2 = -HEOS._delta.pt()*HEOS.dalphar_dDelta()*(1/HEOS._reducing.rhomolar)*(HEOS.Reducing.p->d_ndrhorbardni_dxj__constxi(HEOS.mole_fractions, i, j)-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->drhormolardxi__constxj(HEOS.mole_fractions,j)*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions,i)); + double line3 = HEOS._tau.pt()*d2alphar_dxi_dTau(HEOS, j)*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i); + double line4 = HEOS._tau.pt()*HEOS.dalphar_dTau()*(1/HEOS._reducing.T)*(HEOS.Reducing.p->d_ndTrdni_dxj__constxi(HEOS.mole_fractions,i,j)-1/HEOS._reducing.T*HEOS.Reducing.p->dTrdxi__constxj(HEOS.mole_fractions, j)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i)); + double s = 0; + for (unsigned int m = 0; m < HEOS.mole_fractions.size(); m++) + { + s += HEOS.mole_fractions[m]*d2alphardxidxj(HEOS, j,m); + } + double line5 = d2alphardxidxj(HEOS, i,j)-dalphar_dxi(HEOS, j)-s; + return line1+line2+line3+line4+line5; +} +long double MixtureDerivatives::nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, int i, int j) +{ + double line0 = ndalphar_dni__constT_V_nj(HEOS, j); // First term from 7.46 + double line1 = d_ndalphardni_dDelta(HEOS, i)*nddeltadni__constT_V_nj(HEOS, j); + double line2 = d_ndalphardni_dTau(HEOS, i)*ndtaudni__constT_V_nj(HEOS, j); + double summer = 0; + for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++) + { + summer += HEOS.mole_fractions[k]*d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, k); + } + double line3 = d_ndalphardni_dxj__constdelta_tau_xi(HEOS, i, j)-summer; + return line0 + line1 + line2 + line3; +} +long double MixtureDerivatives::d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + // The first line + double term1 = (HEOS._delta.pt()*HEOS.d2alphar_dDelta2()+HEOS.dalphar_dDelta())*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i)); + + // The second line + double term2 = HEOS._tau.pt()*HEOS.d2alphar_dDelta_dTau()*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i); + + // The third line + double term3 = d2alphar_dxi_dDelta(HEOS, i); + for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++) + { + term3 -= HEOS.mole_fractions[k]*d2alphar_dxi_dDelta(HEOS, k); + } + return term1 + term2 + term3; +} + +long double MixtureDerivatives::d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, int i) +{ + // The first line + double term1 = HEOS._delta.pt()*HEOS.d2alphar_dDelta_dTau()*(1-1/HEOS._reducing.rhomolar*HEOS.Reducing.p->ndrhorbardni__constnj(HEOS.mole_fractions, i)); + + // The second line + double term2 = (HEOS._tau.pt()*HEOS.d2alphar_dTau2()+HEOS.dalphar_dTau())*(1/HEOS._reducing.T)*HEOS.Reducing.p->ndTrdni__constnj(HEOS.mole_fractions, i); + + // The third line + double term3 = d2alphar_dxi_dTau(HEOS, i); + for (unsigned int k = 0; k < HEOS.mole_fractions.size(); k++) + { + term3 -= HEOS.mole_fractions[k]*d2alphar_dxi_dTau(HEOS, k); + } + return term1 + term2 + term3; +} + + +} /* namespace CoolProp */ diff --git a/src/Backends/Helmholtz/MixtureDerivatives.h b/src/Backends/Helmholtz/MixtureDerivatives.h new file mode 100644 index 00000000..db7861df --- /dev/null +++ b/src/Backends/Helmholtz/MixtureDerivatives.h @@ -0,0 +1,296 @@ +/** + * This file contains derivatives needed in the mixture model. The derivatives are quite nasty, and there + * are a lot of them, so they are put in this file for cleanness. The MixtureDerivatives class is a friend class + * of the HelmholtzEOSMixtureBackend, so it can access all the private members of the HelmholtzEOSMixtureBackend + * class +*/ + +// *************************************************************** +// *************************************************************** +// ***************** MIXTURE DERIVATIVES *********************** +// *************************************************************** +// *************************************************************** + +#ifndef MIXTURE_DERIVATIVES_H +#define MIXTURE_DERIVATIVES_H + +#include "HelmholtzEOSMixtureBackend.h" + +namespace CoolProp{ + +/** +This class is a friend class of HelmholtzEOSMixtureBackend, therefore the +static methods contained in it have access to the private and +protected variables in the HelmholtzEOSMixtureBackend instance. + +In this way the derivative routines can be kept in their own separate file +and not pollute the HelmholtzEOSMixtureBackend namespace +*/ +class MixtureDerivatives{ + public: + + static long double dalphar_dxi(HelmholtzEOSMixtureBackend &HEOS, int i); + static long double d2alphar_dxi_dTau(HelmholtzEOSMixtureBackend &HEOS, int i); + static long double d2alphar_dxi_dDelta(HelmholtzEOSMixtureBackend &HEOS, int i); + static long double d2alphardxidxj(HelmholtzEOSMixtureBackend &HEOS, int i, int j); + + /** \brief GERG 2004 Monograph equation 7.61 + * + * The derivative term + * \f[ + * \left(\frac{\partial p}{\partial T} \right)_{V,\bar n} = \rho R(1+\delta \alpha_{\delta}^r-\delta \tau \alpha^r_{\delta\tau}) + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double dpdT__constV_n(HelmholtzEOSMixtureBackend &HEOS); + + static long double dpdrho__constT_n(HelmholtzEOSMixtureBackend &HEOS); + + /** \brief GERG 2004 Monograph equation 7.62 + * + * The derivative term + * \f[ + * n\left(\frac{\partial p}{\partial V} \right)_{T,\bar n} = -\rho^2 RT(1+2\delta \alpha_{\delta}^r+\delta^2\alpha^r_{\delta\delta}) + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double ndpdV__constT_n(HelmholtzEOSMixtureBackend &HEOS); + + /** \brief GERG 2004 Monograph equation 7.63 + * + * The derivative term + * \f[ + * n\left(\frac{\partial p}{\partial n_i} \right)_{T,V,n_j} = \rho RT\left[1+\delta\alpha_{\delta}^r\left[2- \frac{1}{\rho_r}\cdot n\left( \frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] +\delta\cdot n\left(\frac{\partial\alpha_{\delta}^r}{\partial n_i}\right)_{T,V,n_j}\right] + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double ndpdni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief GERG 2004 monograph Eqn. 7.32 + * + * The partial molar volume + * \f[ + * \hat v_i = \left( \frac{\partial V}{\partial n_i}\right)_{T,p,n_j} = \frac{-\left(\dfrac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\dfrac{\partial p}{\partial V}\right)_{T,\bar n}} + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double partial_molar_volume(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief Natural logarithm of the fugacity coefficient + * + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double ln_fugacity_coefficient(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief Derivative of the natural logarithm of the fugacity coefficient with respect to T + * + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double dln_fugacity_coefficient_dT__constrho_n(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief Derivative of the natural logarithm of the fugacity coefficient with respect to T + * + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double dln_fugacity_coefficient_drho__constT_n(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief GERG 2004 Monograph Eqn. 7.29 + * + * The derivative term + * + * \f[ + * \left(\frac{\partial \ln \phi_i}{\partial T} \right)_{p,\bar n} = \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} + \frac{1}{T}-\frac{\hat v}{RT}\left(\frac{\partial p}{\partial T}\right)_{V,\bar n} + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double dln_fugacity_coefficient_dT__constp_n(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief Table B4, Kunz, JCED, 2012 for the original term and the subsequent substitutions + * + * The derivative term + * \f[ + * n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} + * \f] + * which is equal to + * \f{eqnarray*}{ + * n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} &=& \delta \phi^r_{\delta}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\ + * && +\tau \phi^r_{\tau}\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\ + * && +\phi^r_{x_i}-\sum_{k=1}^{N}x_k\phi^r_{x_k} + * \f} + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double ndalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i); + + /// GERG Equation 7.42 + static long double dnalphar_dni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief GERG 2004 Monograph Eqn. 7.30 + * + * The derivative term + * \f[ + * \left(\frac{\partial \ln \phi_i}{\partial p} \right)_{T,\bar n} = \frac{\hat v_i}{RT}-\frac{1}{p} + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double dln_fugacity_coefficient_dp__constT_n(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief GERG 2004 Monograph Equation 7.31 + * + * The derivative term + * \f[ + * n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1+\frac{n}{RT}\frac{\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_j}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}} + * \f] + * which is also equal to + * \f[ + * n\left(\frac{\partial \ln \phi_i}{\partial n_j}\right)_{T,p} = n\left(\frac{\partial^2n\alpha^r}{\partial n_j \partial n_i} \right)_{T,V}+1-\frac{\hat v_i}{RT}\left[n\left(\frac{\partial p}{\partial n_j}\right)_{T,V,n_i}\right] + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double ndln_fugacity_coefficient_dnj__constT_p(HelmholtzEOSMixtureBackend &HEOS, int i, int j); + + /** \brief Gernert Equation 3.115 + * + * The derivative term + * \f[ + * \left(\frac{\partial \ln \phi_i}{\partial x_j}\right)_{T,p,x_{k\neq j}} = \left(\frac{\partial^2n\alpha^r}{\partial x_j \partial n_i} \right)_{T,V}+\frac{1}{RT}\frac{\left(\frac{\partial p}{\partial n_i}\right)_{T,V,n_{k\neq i}}\left(\frac{\partial p}{\partial x_j}\right)_{T,V,x_{k\neq j}}}{\left(\frac{\partial p}{\partial V}\right)_{T,\bar n}} + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double dln_fugacity_coefficient_dxj__constT_p_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j); + + /** \brief Gernert Equation 3.130 + * + * The derivative term + * \f[ + * \left(\frac{\partial p}{\partial x_j} \right)_{T,V,x_{k\neq j}} = \rho RT\left(-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_{k\neq j}} \delta\alpha_{\delta}^r + \delta\left(\frac{\partial}{\partial x_j}\left(\left( \frac{\partial \alpha^r}{\partial \delta}\right)_{\tau,\bar x}\right)\right)_{T,V,x_{k\neq j}}\right) + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double dpdxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j); + + /** \brief Gernert Equation 3.117 + * + * The derivative term + * \f[ + * \left(\frac{\partial^2n\alpha^r}{\partial x_i\partial n_j} \right)_{T,V} = \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_{j\neq i}}\right)\right)_{T,V,x_{k\neq j}} +\left(\frac{\partial \alpha^r}{\partial x_j}\right)_{T,V,x_{k\neq j}} + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double d2nalphar_dxi_dnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, int i, int j){ return MixtureDerivatives::d_ndalphardni_dxj__constT_V_xi(HEOS, i, j) + MixtureDerivatives::dalphar_dxj__constT_V_xi(HEOS, j);}; + + /// Gernert Equation 3.119 + /// Catch test provided + static long double dalphar_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j); + + /// Gernert Equation 3.118 + /// Catch test provided + static long double d_ndalphardni_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j); + + /// Gernert Equation 3.134 + /// Catch test provided + static long double d_dalpharddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j); + + /// Gernert Equation 3.121 + /// Catch test provided + static long double ddelta_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j); + + /// Gernert Equation 3.122 + /// Catch test provided + static long double dtau_dxj__constT_V_xi(HelmholtzEOSMixtureBackend &HEOS, int j); + + /** \brief GERG 2004 Monograph, equations 7.44 and 7.51 + * + * The derivative term + * \f[ + * \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = \left( \frac{\partial}{\partial T}\left(\frac{\partial n \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right)_{V,\bar n} + * \f] + * \f[ + * \left(\frac{\partial^2n\alpha^r}{\partial T\partial n_i} \right)_{V,n_j} = -\frac{\tau}{T}\left[\alpha_{\tau}^r +\left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\bar x}\right] + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double d2nalphar_dni_dT(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief GERG 2004 Monograph Equation 7.51 and Table B4, Kunz, JCED, 2012 + * + * The derivative term + * \f{eqnarray*}{ + * \frac{\partial }{\partial \tau} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right) &=& \delta \phi^r_{\delta\tau}\left[ 1-\frac{1}{\rho_r}\left[\left(\frac{\partial \rho_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial \rho_r}{\partial x_k}\right)_{x_j} \right]\right]\\ + * && +(\tau \phi^r_{\tau\tau}+\phi^r_{\tau})\frac{1}{T_r}\left[\left(\frac{\partial T_r}{\partial x_i}\right)_{x_j} - \sum_{k=1}^N x_k\left(\frac{\partial T_r}{\partial x_k}\right)_{x_j} \right]\\ + * && +\phi^r_{x_i\tau}-\sum_{k=1}^{N}x_k\phi^r_{x_k\tau} + * \f} + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double d_ndalphardni_dTau(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief GERG 2004 Monograph Equation 7.50 and Table B4, Kunz, JCED, 2012 + * + * The derivative term + * \f{eqnarray*}{ + * \left(\frac{\partial }{\partial \delta} \left( n\left(\frac{\partial \phi^r}{\partial n_i} \right)_{T,V,n_j} \right)\right)_{\tau,\bar x} &=& (\alpha_{\delta}^r+\delta\alpha_{\delta\delta}^r)\left[1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j} \right] \\ + * &+&\tau\alpha^r_{\delta\tau}\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\ + * &+&\phi^r_{\delta x_i}-\sum_{k=1}^{N}x_k\phi^r_{\delta x_k} + * \f} + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double d_ndalphardni_dDelta(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief GERG 2004 Monograph equation 7.41 + * + * The derivative term + * \f[ + * n\left(\frac{\partial^2n\alpha^r}{\partial n_i \partial n_j} \right)_{T,V} = n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i} + * \f] + * and + * GERG 2004 Monograph equation 7.46: + * \f[ + * n\left( \frac{\partial}{\partial n_j}\left(\frac{\partial n\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)_{T,V,n_i} = n\left( \frac{\partial \alpha^r}{\partial n_j}\right)_{T,V,n_i} + n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i} + * \f] + * GERG 2004 Monograph equation 7.47: + * \f{eqnarray*}{ + * n\left( \frac{\partial}{\partial n_j}\left(n\left(\frac{\partial \alpha^r}{\partial n_i}\right)_{T,V,n_j} \right) \right)_{T,V,n_i} &=& \left( \frac{\partial}{\partial \delta}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\delta}{\partial n_j}\right)_{T,V,n_i}\\ + * &+& \left( \frac{\partial}{\partial \tau}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\tau,\bar x}\cdot n\left(\frac{\partial\tau}{\partial n_j}\right)_{T,V,n_i}\\ + * &+& \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}-\sum_{k=1}^{N}x_k \left( \frac{\partial}{\partial x_k}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i}\\ + * \f} + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double nd2nalphardnidnj__constT_V(HelmholtzEOSMixtureBackend &HEOS, int i, int j); + + /** \brief GERG 2004 Monograph equation 7.48 + * + * The derivative term + * \f[ + * n\left(\frac{\partial \delta}{\partial n_i} \right)_{T,V,n_j} = \delta - \frac{\delta}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i} \right)_{n_j} + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double nddeltadni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief GERG 2004 Monograph equation 7.49 + * + * The derivative term + * \f[ + * n\left(\frac{\partial \tau}{\partial n_i} \right)_{T,V,n_j} = \frac{\tau}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i} \right)_{n_j} + * \f] + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double ndtaudni__constT_V_nj(HelmholtzEOSMixtureBackend &HEOS, int i); + + /** \brief GERG 2004 Monograph equation 7.52 + * + * The derivative term + * \f{eqnarray*}{ + * \left( \frac{\partial}{\partial x_j}\left(n\left(\frac{\partial\alpha^r}{\partial n_i}\right)_{T,V,n_j}\right)\right)_{\delta,\tau,x_i} &=& \delta\alpha_{\delta x_j}^{r}\left[ 1-\frac{1}{\rho_r}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\ + * &-& \delta\alpha_{\delta}^{r}\frac{1}{\rho_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{\rho_r}\left(\frac{\partial \rho_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial \rho_r}{\partial n_i}\right)_{n_j}\right] \\ + * &+& \tau\alpha_{\tau x_j}^r\frac{1}{T_r}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\\ + * &+& \tau\alpha_{\tau}^{r}\frac{1}{T_r}\left[ \left(\frac{\partial}{\partial x_j}\left(n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right)\right)_{x_i}-\frac{1}{T_r}\left(\frac{\partial T_r}{\partial x_j}\right)_{x_i}\cdot n\left(\frac{\partial T_r}{\partial n_i}\right)_{n_j}\right] \\ + * &+& \alpha_{x_ix_j}^r-\alpha_{x_j}^r-\sum_{m=1}^Nx_m\alpha_{x_jx_m}^r + * \f} + * @param HEOS The HelmholtzEOSMixtureBackend to be used + */ + static long double d_ndalphardni_dxj__constdelta_tau_xi(HelmholtzEOSMixtureBackend &HEOS, int i, int j); + +}; /* class MixtureDerivatives */ + +} /* namepsace CoolProp*/ +#endif \ No newline at end of file diff --git a/src/Backends/Helmholtz/VLERoutines.cpp b/src/Backends/Helmholtz/VLERoutines.cpp index 2ec1b03a..ad80b6f2 100644 --- a/src/Backends/Helmholtz/VLERoutines.cpp +++ b/src/Backends/Helmholtz/VLERoutines.cpp @@ -2,10 +2,11 @@ #include "HelmholtzEOSMixtureBackend.h" #include "VLERoutines.h" #include "MatrixMath.h" +#include "MixtureDerivatives.h" namespace CoolProp { -void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend *HEOS, parameters ykey, long double y){ +void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend &HEOS, parameters ykey, long double y){ class inner_resid : public FuncWrapper1D{ public: @@ -35,9 +36,9 @@ void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend *HEOS, pa long double r, T, rhomolar_liq, rhomolar_vap, value, p, gL, gV, rhomolar_crit; int other; - outer_resid(HelmholtzEOSMixtureBackend *HEOS, CoolProp::parameters ykey, long double y) - : HEOS(HEOS), ykey(ykey), y(y){ - rhomolar_crit = HEOS->rhomolar_critical(); + outer_resid(HelmholtzEOSMixtureBackend &HEOS, CoolProp::parameters ykey, long double y) + : HEOS(&HEOS), ykey(ykey), y(y){ + rhomolar_crit = HEOS.rhomolar_critical(); }; double call(double rhomolar_vap){ this->y = y; @@ -67,13 +68,13 @@ void SaturationSolvers::saturation_critical(HelmholtzEOSMixtureBackend *HEOS, pa }; outer_resid resid(HEOS, iT, y); - double rhomolar_crit = HEOS->rhomolar_critical(); + double rhomolar_crit = HEOS.rhomolar_critical(); std::string errstr; Brent(&resid, rhomolar_crit*(1-1e-8), rhomolar_crit*0.5, DBL_EPSILON, 1e-9, 20, errstr); } -void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options) +void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_options &options) { // Define the residual to be driven to zero @@ -85,8 +86,8 @@ void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double r, T, rhomolar_liq, rhomolar_vap, value, p, gL, gV; int other; - solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double rhomolar_liq_guess, long double rhomolar_vap_guess) - : HEOS(HEOS), T(T), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){}; + solver_resid(HelmholtzEOSMixtureBackend &HEOS, long double T, long double rhomolar_liq_guess, long double rhomolar_vap_guess) + : HEOS(&HEOS), T(T), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){}; double call(double p){ this->p = p; // Recalculate the densities using the current guess values @@ -118,13 +119,13 @@ void SaturationSolvers::saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, Secant(resid, options.p, options.p*1.1, 1e-10, 100, errstr); } catch(std::exception &){ - long double pmax = std::min(options.p*1.03, static_cast(HEOS->p_critical()+1e-6)); - long double pmin = std::max(options.p*0.97, static_cast(HEOS->p_triple()-1e-6)); + long double pmax = std::min(options.p*1.03, static_cast(HEOS.p_critical()+1e-6)); + long double pmin = std::max(options.p*0.97, static_cast(HEOS.p_triple()-1e-6)); Brent(resid, pmin, pmax, LDBL_EPSILON, 1e-8, 100, errstr); } } -void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS, long double p, saturation_PHSU_pure_options &options){ +void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend &HEOS, long double p, saturation_PHSU_pure_options &options){ // Define the residual to be driven to zero class solver_resid : public FuncWrapper1D @@ -135,8 +136,8 @@ void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS, long double r, p, rhomolar_liq, rhomolar_vap, value, T, gL, gV; int other; - solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double p, long double rhomolar_liq_guess, long double rhomolar_vap_guess) - : HEOS(HEOS), p(p), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){}; + solver_resid(HelmholtzEOSMixtureBackend &HEOS, long double p, long double rhomolar_liq_guess, long double rhomolar_vap_guess) + : HEOS(&HEOS), p(p), rhomolar_liq(rhomolar_liq_guess), rhomolar_vap(rhomolar_vap_guess){}; double call(double T){ this->T = T; // Recalculate the densities using the current guess values @@ -160,12 +161,12 @@ void SaturationSolvers::saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS, if (!ValidNumber(options.rhoV)){throw ValueError("options.rhoV is not valid in saturation_P_pure_1D_T");}; std::string errstr; - long double Tmax = std::min(options.T + 2, static_cast(HEOS->T_critical()-1e-6)); - long double Tmin = std::max(options.T - 2, static_cast(HEOS->Ttriple()+1e-6)); + long double Tmax = std::min(options.T + 2, static_cast(HEOS.T_critical()-1e-6)); + long double Tmin = std::max(options.T - 2, static_cast(HEOS.Ttriple()+1e-6)); Brent(resid, Tmin, Tmax, LDBL_EPSILON, 1e-11, 100, errstr); } -void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, long double specified_value, saturation_PHSU_pure_options &options) +void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, long double specified_value, saturation_PHSU_pure_options &options) { /* This function is inspired by the method of Akasaka: @@ -179,10 +180,10 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l std::vector negativer(3,_HUGE), v; std::vector > J(3, std::vector(3,_HUGE)); - HEOS->calc_reducing_state(); - const SimpleState & reduce = HEOS->get_reducing_state(); - shared_ptr SatL = HEOS->SatL, - SatV = HEOS->SatV; + HEOS.calc_reducing_state(); + const SimpleState & reduce = HEOS.get_reducing_state(); + shared_ptr SatL = HEOS.SatL, + SatV = HEOS.SatV; long double T, rhoL, rhoV, pL, pV; long double deltaL=0, deltaV=0, tau=0, error; @@ -196,7 +197,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l // Invert liquid density ancillary to get temperature // TODO: fit inverse ancillaries too try{ - T = HEOS->get_components()[0]->ancillaries.pL.invert(specified_value); + T = HEOS.get_components()[0]->ancillaries.pL.invert(specified_value); } catch(std::exception &e) { @@ -207,11 +208,11 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l { throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid",options.specified_variable)); } - if (T > HEOS->T_critical()-1){ T -= 1; } + if (T > HEOS.T_critical()-1){ T -= 1; } // Evaluate densities from the ancillary equations - rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T); - rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T); + rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T); + rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T); // Apply a single step of Newton's method to improve guess value for liquid // based on the error between the gas pressure (which is usually very close already) @@ -372,7 +373,7 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, l } while (error > 1e-9); } -void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long double rhomolar, saturation_D_pure_options &options) +void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend &HEOS, long double rhomolar, saturation_D_pure_options &options) { /* This function is inspired by the method of Akasaka: @@ -386,10 +387,10 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long std::vector r(2,_HUGE), v; std::vector > J(2, std::vector(2,_HUGE)); - HEOS->calc_reducing_state(); - const SimpleState & reduce = HEOS->get_reducing_state(); - shared_ptr SatL = HEOS->SatL, - SatV = HEOS->SatV; + HEOS.calc_reducing_state(); + const SimpleState & reduce = HEOS.get_reducing_state(); + shared_ptr SatL = HEOS.SatL, + SatV = HEOS.SatV; long double T, rhoL,rhoV; long double deltaL=0, deltaV=0, tau=0, error, p_error; @@ -402,16 +403,16 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long { // Invert liquid density ancillary to get temperature // TODO: fit inverse ancillaries too - T = HEOS->get_components()[0]->ancillaries.rhoL.invert(rhomolar); - rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T); + T = HEOS.get_components()[0]->ancillaries.rhoL.invert(rhomolar); + rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T); rhoL = rhomolar; } else if (options.imposed_rho == saturation_D_pure_options::IMPOSED_RHOV) { // Invert vapor density ancillary to get temperature // TODO: fit inverse ancillaries too - T = HEOS->get_components()[0]->ancillaries.rhoV.invert(rhomolar); - rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T); + T = HEOS.get_components()[0]->ancillaries.rhoV.invert(rhomolar); + rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T); rhoV = rhomolar; } else @@ -532,7 +533,7 @@ void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long throw SolutionError(format("saturation_D_pure solver abs error on p [%Lg] > limit [%Lg]", p_error, p_error_limit)); } } -void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options) +void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_options &options) { // Set some imput options SaturationSolvers::saturation_T_pure_Akasaka_options _options; @@ -552,7 +553,7 @@ void SaturationSolvers::saturation_T_pure(HelmholtzEOSMixtureBackend *HEOS, long SaturationSolvers::saturation_T_pure_1D_P(HEOS, T, options); } } -void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_Akasaka_options &options) +void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_Akasaka_options &options) { // Start with the method of Akasaka @@ -566,11 +567,11 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE Ancillary equations are used to get a sensible starting point */ - HEOS->calc_reducing_state(); - const SimpleState & reduce = HEOS->get_reducing_state(); - long double R_u = HEOS->calc_gas_constant(); - shared_ptr SatL = HEOS->SatL, - SatV = HEOS->SatV; + HEOS.calc_reducing_state(); + const SimpleState & reduce = HEOS.get_reducing_state(); + long double R_u = HEOS.calc_gas_constant(); + shared_ptr SatL = HEOS.SatL, + SatV = HEOS.SatV; long double rhoL,rhoV,JL,JV,KL,KV,dJL,dJV,dKL,dKV; long double DELTA, deltaL=0, deltaV=0, error, PL, PV, stepL, stepV; @@ -589,13 +590,13 @@ void SaturationSolvers::saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HE // Use the density ancillary function as the starting point for the solver // If very close to the critical temp, evaluate the ancillaries for a slightly lower temperature - if (T > 0.99*HEOS->get_reducing_state().T){ - rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T-0.1); - rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T-0.1); + if (T > 0.99*HEOS.get_reducing_state().T){ + rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T-0.1); + rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T-0.1); } else{ - rhoL = HEOS->get_components()[0]->ancillaries.rhoL.evaluate(T); - rhoV = HEOS->get_components()[0]->ancillaries.rhoV.evaluate(T); + rhoL = HEOS.get_components()[0]->ancillaries.rhoL.evaluate(T); + rhoV = HEOS.get_components()[0]->ancillaries.rhoV.evaluate(T); // Apply a single step of Newton's method to improve guess value for liquid // based on the error between the gas pressure (which is usually very close already) @@ -738,16 +739,16 @@ void SaturationSolvers::successive_substitution(HelmholtzEOSMixtureBackend &HEOS for (std::size_t i = 0; i < N; ++i) { - ln_phi_liq[i] = HEOS.SatL->mixderiv_ln_fugacity_coefficient(i); - ln_phi_vap[i] = HEOS.SatV->mixderiv_ln_fugacity_coefficient(i); + ln_phi_liq[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatL.get()), i); + ln_phi_vap[i] = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatV.get()), i); if (options.sstype == imposed_p){ - deriv_liq = HEOS.SatL->mixderiv_dln_fugacity_coefficient_dT__constp_n(i); - deriv_vap = HEOS.SatV->mixderiv_dln_fugacity_coefficient_dT__constp_n(i); + deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatL.get()), i); + deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dT__constp_n(*(HEOS.SatV.get()), i); } else if (options.sstype == imposed_T){ - deriv_liq = HEOS.SatL->mixderiv_dln_fugacity_coefficient_dp__constT_n(i); - deriv_vap = HEOS.SatV->mixderiv_dln_fugacity_coefficient_dp__constT_n(i); + deriv_liq = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatL.get()), i); + deriv_vap = MixtureDerivatives::dln_fugacity_coefficient_dp__constT_n(*(HEOS.SatV.get()), i); } else {throw ValueError();} @@ -811,15 +812,15 @@ void SaturationSolvers::newton_raphson_VLE_GV::resize(unsigned int N) // Last entry is 1 neg_dFdS[N+1] = 1.0; } -void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtureBackend *HEOS, const std::vector &z, std::vector &K, mixture_VLE_IO &IO) +void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtureBackend &HEOS, const std::vector &z, std::vector &K, mixture_VLE_IO &IO) { // Reset all the variables and resize pre_call(); std::size_t N = K.size(); resize(N); - shared_ptr SatL(new HelmholtzEOSMixtureBackend(HEOS->get_components())), - SatV(new HelmholtzEOSMixtureBackend(HEOS->get_components())); + shared_ptr SatL(new HelmholtzEOSMixtureBackend(HEOS.get_components())), + SatV(new HelmholtzEOSMixtureBackend(HEOS.get_components())); SatL->specify_phase(iphase_liquid); SatV->specify_phase(iphase_gas); @@ -878,7 +879,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::check_Jacobian(HelmholtzEOSMixtur std::cout << vec_to_string(get_col(J,N+1),"%12.11f") << std::endl; } -void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend *HEOS, const std::vector &z, std::vector &K, mixture_VLE_IO &IO) +void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend &HEOS, const std::vector &z, std::vector &K, mixture_VLE_IO &IO) { int iter = 0; @@ -886,8 +887,8 @@ void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend * pre_call(); resize(K.size()); - shared_ptr SatL(new HelmholtzEOSMixtureBackend(HEOS->get_components())), - SatV(new HelmholtzEOSMixtureBackend(HEOS->get_components())); + shared_ptr SatL(new HelmholtzEOSMixtureBackend(HEOS.get_components())), + SatV(new HelmholtzEOSMixtureBackend(HEOS.get_components())); SatL->specify_phase(iphase_liquid); // So it will always just use single-phase solution SatV->specify_phase(iphase_gas); // So it will always just use single-phase solution @@ -931,7 +932,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::call(HelmholtzEOSMixtureBackend * IO.y = &y; // Mole fractions in vapor } -void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureBackend *HEOS, long double beta, long double T, long double rhomolar_liq, const long double rhomolar_vap, const std::vector &z, std::vector &K) +void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureBackend &HEOS, long double beta, long double T, long double rhomolar_liq, const long double rhomolar_vap, const std::vector &z, std::vector &K) { // Step 0: // -------- @@ -941,6 +942,8 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB // Set the mole fractions in the classes SatL->set_mole_fractions(x); SatV->set_mole_fractions(y); + + HelmholtzEOSMixtureBackend &rSatV = *SatV, &rSatL = *SatL; // Update the liquid and vapor classes SatL->update(DmolarT_INPUTS, rhomolar_liq, T); @@ -958,22 +961,22 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB // For the residuals F_i for (unsigned int i = 0; i < N; ++i) { - long double ln_phi_liq = SatL->mixderiv_ln_fugacity_coefficient(i); - long double phi_iT_liq = SatL->mixderiv_dln_fugacity_coefficient_dT__constrho_n(i); - dlnphi_drho_liq[i] = SatL->mixderiv_dln_fugacity_coefficient_drho__constT_n(i); + long double ln_phi_liq = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatL.get()), i); + long double phi_iT_liq = MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(*(HEOS.SatL.get()), i); + dlnphi_drho_liq[i] = MixtureDerivatives::dln_fugacity_coefficient_drho__constT_n(*(HEOS.SatL.get()), i); for (unsigned int j = 0; j < N; ++j) { // I think this is wrong. - phi_ij_liq[j] = SatL->mixderiv_ndln_fugacity_coefficient_dnj__constT_p(i,j) + (SatL->mixderiv_partial_molar_volume(i)/(SatL->gas_constant()*T)-1/p)*SatL->mixderiv_ndpdni__constT_V_nj(i); // 7.126 from GERG monograph + phi_ij_liq[j] = MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(rSatL, i, j) + (MixtureDerivatives::partial_molar_volume(rSatL, i)/(SatL->gas_constant()*T)-1/p)*MixtureDerivatives::ndpdni__constT_V_nj(rSatL, i); // 7.126 from GERG monograph } - long double ln_phi_vap = SatV->mixderiv_ln_fugacity_coefficient(i); - long double phi_iT_vap = SatV->mixderiv_dln_fugacity_coefficient_dT__constrho_n(i); - dlnphi_drho_vap[i] = SatV->mixderiv_dln_fugacity_coefficient_drho__constT_n(i); + long double ln_phi_vap = MixtureDerivatives::ln_fugacity_coefficient(*(HEOS.SatV.get()), i); + long double phi_iT_vap = MixtureDerivatives::dln_fugacity_coefficient_dT__constrho_n(*(HEOS.SatV.get()), i); + dlnphi_drho_vap[i] = MixtureDerivatives::dln_fugacity_coefficient_drho__constT_n(*(HEOS.SatV.get()), i); for (unsigned int j = 0; j < N; ++j) { // I think this is wrong. - phi_ij_vap[j] = SatV->mixderiv_ndln_fugacity_coefficient_dnj__constT_p(i,j) + (SatV->mixderiv_partial_molar_volume(i)/(SatV->gas_constant()*T)-1/p)*SatV->mixderiv_ndpdni__constT_V_nj(i); ; // 7.126 from GERG monograph + phi_ij_vap[j] = MixtureDerivatives::ndln_fugacity_coefficient_dnj__constT_p(rSatV, i,j) + (MixtureDerivatives::partial_molar_volume(rSatV, i)/(SatV->gas_constant()*T)-1/p)*MixtureDerivatives::ndpdni__constT_V_nj(rSatV, i); ; // 7.126 from GERG monograph } r[i] = log(K[i]) + ln_phi_vap - ln_phi_liq; @@ -1009,12 +1012,12 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB r[N+1] = p_liq-p_vap; for (unsigned int j = 0; j < N; ++j) { - J[N+1][j] = HEOS->gas_constant()*T*K[j]*z[j]/pow(1-beta+beta*K[j],(int)2)*((1-beta)*dlnphi_drho_vap[j]+beta*dlnphi_drho_liq[j]); + J[N+1][j] = HEOS.gas_constant()*T*K[j]*z[j]/pow(1-beta+beta*K[j],(int)2)*((1-beta)*dlnphi_drho_vap[j]+beta*dlnphi_drho_liq[j]); } // dF_{N+1}/d(ln(T)) - J[N+1][N] = T*(SatL->mixderiv_dpdT__constV_n() - SatV->mixderiv_dpdT__constV_n()); + J[N+1][N] = T*(MixtureDerivatives::dpdT__constV_n(*(HEOS.SatL.get())) - MixtureDerivatives::dpdT__constV_n(*(HEOS.SatV.get()))); // dF_{N+1}/d(ln(rho')) - J[N+1][N+1] = rhomolar_liq*SatL->mixderiv_dpdrho__constT_n(); + J[N+1][N+1] = rhomolar_liq*MixtureDerivatives::dpdrho__constT_n(*(HEOS.SatL.get())); // Flip all the signs of the entries in the residual vector since we are solving Jv = -r, not Jv=r // Also calculate the rms error of the residual vector at this step @@ -1027,7 +1030,7 @@ void SaturationSolvers::newton_raphson_VLE_GV::build_arrays(HelmholtzEOSMixtureB error_rms = sqrt(error_rms); // Square-root (The R in RMS) } -void PhaseEnvelope::PhaseEnvelope_GV::build(HelmholtzEOSMixtureBackend *HEOS, const std::vector &z, std::vector &K, SaturationSolvers::mixture_VLE_IO &IO) +void PhaseEnvelope::PhaseEnvelope_GV::build(HelmholtzEOSMixtureBackend &HEOS, const std::vector &z, std::vector &K, SaturationSolvers::mixture_VLE_IO &IO) { // Use the residual function based on ln(K_i), ln(T) and ln(rho') as independent variables. rho'' is specified SaturationSolvers::newton_raphson_VLE_GV NRVLE; diff --git a/src/Backends/Helmholtz/VLERoutines.h b/src/Backends/Helmholtz/VLERoutines.h index 05b4d06c..1ebca6dc 100644 --- a/src/Backends/Helmholtz/VLERoutines.h +++ b/src/Backends/Helmholtz/VLERoutines.h @@ -45,14 +45,14 @@ namespace SaturationSolvers @param p Pressure [Pa] @param i Index of component [-] */ - static long double Wilson_lnK_factor(HelmholtzEOSMixtureBackend *HEOS, long double T, long double p, int i){ - EquationOfState *EOS = (HEOS->get_components())[i]->pEOS; + static long double Wilson_lnK_factor(HelmholtzEOSMixtureBackend &HEOS, long double T, long double p, int i){ + EquationOfState *EOS = (HEOS.get_components())[i]->pEOS; return log(EOS->reduce.p/p)+5.373*(1 + EOS->accentric)*(1-EOS->reduce.T/T); }; - void saturation_D_pure(HelmholtzEOSMixtureBackend *HEOS, long double rhomolar, saturation_D_pure_options &options); - void saturation_T_pure(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options); - void saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_Akasaka_options &options); + void saturation_D_pure(HelmholtzEOSMixtureBackend &HEOS, long double rhomolar, saturation_D_pure_options &options); + void saturation_T_pure(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_options &options); + void saturation_T_pure_Akasaka(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_Akasaka_options &options); /** */ @@ -67,7 +67,7 @@ namespace SaturationSolvers /** */ - void saturation_PHSU_pure(HelmholtzEOSMixtureBackend *HEOS, long double specified_value, saturation_PHSU_pure_options &options); + void saturation_PHSU_pure(HelmholtzEOSMixtureBackend &HEOS, long double specified_value, saturation_PHSU_pure_options &options); /* \brief This is a backup saturation_p solver for the case where the Newton solver cannot approach closely enough the solution * @@ -77,7 +77,7 @@ namespace SaturationSolvers * @param p Imposed pressure in kPa * @param options Options to be passed to the function (at least T, rhoL and rhoV must be provided) */ - void saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend *HEOS, long double p, saturation_PHSU_pure_options &options); + void saturation_P_pure_1D_T(HelmholtzEOSMixtureBackend &HEOS, long double p, saturation_PHSU_pure_options &options); /* \brief This is a backup saturation_T solver for the case where the Newton solver cannot approach closely enough the solution * @@ -87,7 +87,7 @@ namespace SaturationSolvers * @param T Imposed temperature in K * @param options Options to be passed to the function (at least p, rhoL and rhoV must be provided) */ - void saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend *HEOS, long double T, saturation_T_pure_options &options); + void saturation_T_pure_1D_P(HelmholtzEOSMixtureBackend &HEOS, long double T, saturation_T_pure_options &options); /* \brief A robust but slow solver in the very-near-critical region * @@ -102,7 +102,7 @@ namespace SaturationSolvers * @param ykey The CoolProp::parameters key to be imposed - one of iT or iP * @param y The value for the imposed variable */ - void saturation_critical(HelmholtzEOSMixtureBackend *HEOS, CoolProp::parameters ykey, long double y); + void saturation_critical(HelmholtzEOSMixtureBackend &HEOS, CoolProp::parameters ykey, long double y); void successive_substitution(HelmholtzEOSMixtureBackend &HEOS, const long double beta, @@ -127,8 +127,8 @@ namespace SaturationSolvers std::vector *K; HelmholtzEOSMixtureBackend *HEOS; - WilsonK_resid(HelmholtzEOSMixtureBackend *HEOS, double beta, double imposed_value, int input_type, const std::vector &z, std::vector &K){ - this->z = &z; this->K = &K; this->HEOS = HEOS; this->beta = beta; this->input_type = input_type; + WilsonK_resid(HelmholtzEOSMixtureBackend &HEOS, double beta, double imposed_value, int input_type, const std::vector &z, std::vector &K){ + this->z = &z; this->K = &K; this->HEOS = &HEOS; this->beta = beta; this->input_type = input_type; if (input_type == imposed_T){ this->T = imposed_value; } @@ -145,19 +145,19 @@ namespace SaturationSolvers T = input_value; // Iterate on temperature, pressure imposed } for (unsigned int i = 0; i< (*z).size(); i++) { - (*K)[i] = exp(Wilson_lnK_factor(HEOS,T,p,i)); + (*K)[i] = exp(Wilson_lnK_factor(*HEOS,T,p,i)); summer += (*z)[i]*((*K)[i]-1)/(1-beta+beta*(*K)[i]); } return summer; }; }; - inline double saturation_preconditioner(HelmholtzEOSMixtureBackend *HEOS, double input_value, int input_type, const std::vector &z) + inline double saturation_preconditioner(HelmholtzEOSMixtureBackend &HEOS, double input_value, int input_type, const std::vector &z) { double ptriple = 0, pcrit = 0, Ttriple = 0, Tcrit = 0; for (unsigned int i = 0; i < z.size(); i++) { - EquationOfState *EOS = (HEOS->get_components())[i]->pEOS; + EquationOfState *EOS = (HEOS.get_components())[i]->pEOS; ptriple += EOS->sat_min_liquid.p*z[i]; pcrit += EOS->reduce.p*z[i]; @@ -175,14 +175,14 @@ namespace SaturationSolvers } else{ throw ValueError();} } - inline double saturation_Wilson(HelmholtzEOSMixtureBackend *HEOS, double beta, double input_value, int input_type, const std::vector &z, double guess) + inline double saturation_Wilson(HelmholtzEOSMixtureBackend &HEOS, double beta, double input_value, int input_type, const std::vector &z, double guess) { double T; std::string errstr; // Find first guess for T using Wilson K-factors - WilsonK_resid Resid(HEOS, beta, input_value, input_type, z, HEOS->get_K()); + WilsonK_resid Resid(HEOS, beta, input_value, input_type, z, HEOS.get_K()); T = Secant(Resid, guess, 0.001, 1e-10, 100, errstr); if (!ValidNumber(T)){throw ValueError("saturation_p_Wilson failed to get good T");} @@ -234,7 +234,7 @@ namespace SaturationSolvers @param z Bulk mole fractions [-] @param K Array of K-factors [-] */ - void call(HelmholtzEOSMixtureBackend *HEOS, const std::vector &z, std::vector &K, mixture_VLE_IO &IO); + void call(HelmholtzEOSMixtureBackend &HEOS, const std::vector &z, std::vector &K, mixture_VLE_IO &IO); /*! Build the arrays for the Newton-Raphson solve @@ -246,11 +246,11 @@ namespace SaturationSolvers @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 &z, std::vector & K); + void build_arrays(HelmholtzEOSMixtureBackend &HEOS, long double beta, long double T, long double rhomolar_liq, const long double rho_vapor, const std::vector &z, std::vector & K); /** Check the derivatives in the Jacobian using numerical derivatives. */ - void check_Jacobian(HelmholtzEOSMixtureBackend *HEOS, const std::vector &z, std::vector &K, mixture_VLE_IO &IO); + void check_Jacobian(HelmholtzEOSMixtureBackend &HEOS, const std::vector &z, std::vector &K, mixture_VLE_IO &IO); }; }; @@ -299,7 +299,7 @@ namespace PhaseEnvelope public: PhaseEnvelopeLog bubble, dew; - void build(HelmholtzEOSMixtureBackend *HEOS, const std::vector &z, std::vector &K, SaturationSolvers::mixture_VLE_IO &IO); + void build(HelmholtzEOSMixtureBackend &HEOS, const std::vector &z, std::vector &K, SaturationSolvers::mixture_VLE_IO &IO); }; };