diff --git a/src/Backends/Helmholtz/ExcessHEFunction.cpp b/src/Backends/Helmholtz/ExcessHEFunction.cpp index e442f298..164e4c91 100644 --- a/src/Backends/Helmholtz/ExcessHEFunction.cpp +++ b/src/Backends/Helmholtz/ExcessHEFunction.cpp @@ -1,3 +1,5 @@ +#include + #include "ExcessHEFunction.h" #include "mixture_excess_term_JSON.h" @@ -44,12 +46,12 @@ MixtureExcessHELibrary::MixtureExcessHELibrary() for (unsigned int i = 0; i < CAS1V.size(); ++i) { - // Get the vector of CAS numbers + // Get the vector of CAS numbers std::vector CAS; CAS.resize(2); CAS[0] = CAS1V[i]; CAS[1] = CAS2V[i]; - + // Sort the CAS number vector std::sort(CAS.begin(), CAS.end()); @@ -78,7 +80,7 @@ MixtureExcessHELibrary::MixtureExcessHELibrary() { throw ValueError(); } - + // If not in map, add new entry to map with dictionary if (excess_map.find(CAS) == excess_map.end()) { @@ -93,10 +95,10 @@ MixtureExcessHELibrary::MixtureExcessHELibrary() { // Append dictionary to listing excess_map[CAS].push_back(d); - } + } } } - } + } } void ExcessTerm::construct(const std::vector &components) @@ -105,7 +107,7 @@ void ExcessTerm::construct(const std::vector &components) N = components.size(); F.resize(N, std::vector(N, 0)); - DepartureFunctionMatrix.resize(N, std::vector(N, NULL)); + DepartureFunctionMatrix.resize(N); for (unsigned int i = 0; i < N; ++i) { @@ -162,7 +164,7 @@ double ExcessTerm::alphar(double tau, double delta, const std::vectoralphar(tau,delta); } } @@ -174,7 +176,7 @@ double ExcessTerm::dalphar_dTau(double tau, double delta, const std::vectordalphar_dTau(tau,delta); } } @@ -186,7 +188,7 @@ double ExcessTerm::dalphar_dDelta(double tau, double delta, const std::vectordalphar_dDelta(tau,delta); } } @@ -198,7 +200,7 @@ double ExcessTerm::d2alphar_dDelta2(double tau, double delta, const std::vector< for (unsigned int i = 0; i < N-1; i++) { for (unsigned int j = i + 1; j < N; j++) - { + { summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta2(tau,delta); } } @@ -210,7 +212,7 @@ double ExcessTerm::d2alphar_dTau2(double tau, double delta, const std::vectord2alphar_dTau2(tau,delta); } } @@ -222,7 +224,7 @@ double ExcessTerm::d2alphar_dDelta_dTau(double tau, double delta, const std::vec for (unsigned int i = 0; i < N-1; i++) { for (unsigned int j = i + 1; j < N; j++) - { + { summer += x[i]*x[j]*F[i][j]*DepartureFunctionMatrix[i][j]->d2alphar_dDelta_dTau(tau,delta); } } @@ -280,7 +282,7 @@ GERG2008DepartureFunction::GERG2008DepartureFunction(const std::vector & const std::vector &eta,const std::vector &epsilon,const std::vector &beta, const std::vector &gamma, int Npower) { - + /// Break up into power and gaussian terms { std::vector _n(n.begin(), n.begin()+Npower); diff --git a/src/Backends/Helmholtz/ExcessHEFunction.h b/src/Backends/Helmholtz/ExcessHEFunction.h index df1ceec7..99e68f93 100644 --- a/src/Backends/Helmholtz/ExcessHEFunction.h +++ b/src/Backends/Helmholtz/ExcessHEFunction.h @@ -1,6 +1,7 @@ #ifndef EXCESSHE_FUNCTIONS_H #define EXCESSHE_FUNCTIONS_H +#include #include #include "CoolPropFluid.h" @@ -35,7 +36,7 @@ public: d.add_number("F", F[i]); } d.add_number("Npower", cpjson::get_double(val,"Npower")); - + // Terms for the power d.add_double_vector("n", cpjson::get_double_array(val["n"])); d.add_double_vector("d", cpjson::get_double_array(val["d"])); @@ -59,7 +60,7 @@ public: }; }; -/*! +/*! The abstract base class for departure functions for the excess part of the Helmholtz energy */ class DepartureFunction @@ -67,7 +68,7 @@ class DepartureFunction public: DepartureFunction(){}; virtual ~DepartureFunction(){}; - + /// The excess Helmholtz energy of the binary pair /// Pure-virtual function (must be implemented in derived class virtual double alphar(double tau, double delta) = 0; @@ -86,7 +87,7 @@ public: unsigned int N; std::vector > DepartureFunctionMatrix; std::vector > F; - + ExcessTerm(){}; void construct(const std::vector &components); ~ExcessTerm(); @@ -124,4 +125,4 @@ public: }; } /* namespace CoolProp */ -#endif \ No newline at end of file +#endif diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp index 007ce30c..89fc7dde 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.cpp @@ -5,6 +5,8 @@ * Author: jowr */ +#include + #if defined(_MSC_VER) #define _CRTDBG_MAP_ALLOC #define _CRT_SECURE_NO_WARNINGS @@ -67,7 +69,7 @@ void HelmholtzEOSMixtureBackend::set_components(std::vector comp imposed_phase_index = -1; - // Top-level class can hold copies of the base saturation classes, + // Top-level class can hold copies of the base saturation classes, // saturation classes cannot hold copies of the saturation classes if (generate_SatL_and_SatV) { @@ -76,10 +78,6 @@ void HelmholtzEOSMixtureBackend::set_components(std::vector comp SatV.reset(new HelmholtzEOSMixtureBackend(components, false)); SatV->specify_phase(iphase_gas); } - else - { - SatL = NULL; SatV = NULL; - } } void HelmholtzEOSMixtureBackend::set_mole_fractions(const std::vector &mole_fractions) { @@ -154,7 +152,7 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity_dilute(void) { throw NotImplementedError(format("dilute viscosity not implemented for mixtures")); } - + } long double HelmholtzEOSMixtureBackend::calc_viscosity_background() { @@ -203,8 +201,9 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity(void) { // Get reference fluid name std::string fluid_name = component.transport.viscosity_ecs.reference_fluid; + std::vector names(1, fluid_name); // Get a managed pointer to the reference fluid for ECS - std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(std::vector(1, fluid_name))); + std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(names)); // Get the viscosity using ECS return TransportRoutines::viscosity_ECS(*this, *(ref_fluid.get())); } @@ -266,8 +265,9 @@ long double HelmholtzEOSMixtureBackend::calc_conductivity(void) { // Get reference fluid name std::string fluid_name = component.transport.conductivity_ecs.reference_fluid; + std::vector name(1, fluid_name); // Get a managed pointer to the reference fluid for ECS - std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(std::vector(1, fluid_name))); + std::tr1::shared_ptr ref_fluid(new HelmholtzEOSMixtureBackend(name,false)); // Get the viscosity using ECS return TransportRoutines::conductivity_ECS(*this, *(ref_fluid.get())); } @@ -304,7 +304,7 @@ long double HelmholtzEOSMixtureBackend::calc_conductivity(void) } long double lambda_residual = calc_conductivity_background(); - + // Critical part long double lambda_critical = _HUGE; switch(component.transport.conductivity_critical.type) @@ -412,7 +412,7 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double & case PSmass_INPUTS: ///< Pressure in Pa, Entropy in J/kg/K case PUmass_INPUTS: ///< Pressure in Pa, Internal energy in J/kg case HmassSmass_INPUTS: ///< Enthalpy in J/kg, Entropy in J/kg/K - case SmassUmass_INPUTS: ///< Entropy in J/kg/K, Internal energy in J/kg + case SmassUmass_INPUTS: ///< Entropy in J/kg/K, Internal energy in J/kg case DmassHmass_INPUTS: ///< Mass density in kg/m^3, Enthalpy in J/kg case DmassSmass_INPUTS: ///< Mass density in kg/m^3, Entropy in J/kg/K case DmassUmass_INPUTS: ///< Mass density in kg/m^3, Internal energy in J/kg @@ -447,15 +447,15 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(long &input_pair, double & } void HelmholtzEOSMixtureBackend::update(long input_pair, double value1, double value2 ) { - clear(); + clear(); - if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) { - throw ValueError("Mole fractions must be set"); + if (is_pure_or_pseudopure == false && mole_fractions.size() == 0) { + throw ValueError("Mole fractions must be set"); } mass_to_molar_inputs(input_pair, value1, value2); - // Set the mole-fraction weighted gas constant for the mixture + // Set the mole-fraction weighted gas constant for the mixture // (or the pure/pseudo-pure fluid) if it hasn't been set yet gas_constant(); @@ -598,7 +598,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot { long double p_vap = 0.98*static_cast(_pVanc); long double p_liq = 1.02*static_cast(_pLanc); - + if (value < p_vap){ this->_phase = iphase_gas; _Q = -1000; return; } @@ -636,7 +636,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot switch (other) { case iSmolar: - { + { if (value > SatV->calc_smolar()){ this->_phase = iphase_gas; return; } @@ -671,7 +671,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot } } } - + // Determine Q based on the input provided if (!is_pure_or_pseudopure){throw ValueError("possibly two-phase inputs not supported for pseudo-pure for now");} @@ -683,7 +683,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot SaturationSolvers::saturation_T_pure(&HEOS, _T, options); long double Q; - + if (other == iP) { if (value > 100*DBL_EPSILON + HEOS.SatL->p()){ @@ -721,7 +721,7 @@ void HelmholtzEOSMixtureBackend::p_phase_determination_pure_or_pseudopure(int ot this->_phase = iphase_twophase; } _Q = Q; - // Load the outputs + // Load the outputs _p = _Q*HEOS.SatV->p() + (1-_Q)*HEOS.SatL->p(); _rhomolar = 1/(_Q/HEOS.SatV->rhomolar() + (1-_Q)/HEOS.SatL->rhomolar()); return; @@ -747,7 +747,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot _pVanc = components[0]->ancillaries.pV.evaluate(_T); long double p_vap = 0.98*static_cast(_pVanc); long double p_liq = 1.02*static_cast(_pLanc); - + if (value < p_vap){ this->_phase = iphase_gas; _Q = -1000; return; } @@ -785,7 +785,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot switch (other) { case iSmolar: - { + { if (value > SatV->calc_smolar()){ this->_phase = iphase_gas; return; } @@ -820,7 +820,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot } } } - + // Determine Q based on the input provided if (!is_pure_or_pseudopure){throw ValueError("possibly two-phase inputs not supported for pseudo-pure for now");} @@ -832,7 +832,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot SaturationSolvers::saturation_T_pure(&HEOS, _T, options); long double Q; - + if (other == iP) { if (value > 100*DBL_EPSILON + HEOS.SatL->p()){ @@ -870,7 +870,7 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot this->_phase = iphase_twophase; } _Q = Q; - // Load the outputs + // Load the outputs _p = _Q*HEOS.SatV->p() + (1-_Q)*HEOS.SatL->p(); _rhomolar = 1/(_Q/HEOS.SatV->rhomolar() + (1-_Q)/HEOS.SatL->rhomolar()); return; @@ -962,10 +962,10 @@ void HelmholtzEOSMixtureBackend::T_phase_determination_pure_or_pseudopure(int ot // this->_phase = iphase_liquid; // } // else if (Q > 1+100*DBL_EPSILON){ -// this->_phase = iphase_gas; +// this->_phase = iphase_gas; // } // else{ -// this->_phase = iphase_twophase; +// this->_phase = iphase_twophase; // } // _Q = Q; // // Load the outputs @@ -1070,12 +1070,12 @@ void get_dtau_ddelta(HelmholtzEOSMixtureBackend *HEOS, long double T, long doubl // dp/ddelta|tau ddelta = rhor*R*T*(1+2*delta*dalphar_dDelta+pow(delta, 2)*d2alphar_dDelta2); // dp/dtau|delta - dtau = dT_dtau*rho*R*(1+delta*dalphar_dDelta-tau*delta*d2alphar_dDelta_dTau); + dtau = dT_dtau*rho*R*(1+delta*dalphar_dDelta-tau*delta*d2alphar_dDelta_dTau); break; } case iHmolar: // dh/dtau|delta - dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()) + (1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); + dtau = dT_dtau*R*(-pow(tau,2)*(HEOS->d2alpha0_dTau2()+HEOS->d2alphar_dTau2()) + (1+delta*HEOS->dalphar_dDelta()-tau*delta*HEOS->d2alphar_dDelta_dTau())); // dh/ddelta|tau ddelta = rhor*T*R/rho*(tau*delta*HEOS->d2alphar_dDelta_dTau()+delta*HEOS->dalphar_dDelta()+pow(delta,2)*HEOS->d2alphar_dDelta2()); break; @@ -1102,7 +1102,7 @@ void get_dtau_ddelta(HelmholtzEOSMixtureBackend *HEOS, long double T, long doubl long double HelmholtzEOSMixtureBackend::calc_first_partial_deriv_nocache(long double T, long double rhomolar, int Of, int Wrt, int Constant) { long double dOf_dtau, dOf_ddelta, dWrt_dtau, dWrt_ddelta, dConstant_dtau, dConstant_ddelta; - + get_dtau_ddelta(this, T, rhomolar, Of, dOf_dtau, dOf_ddelta); get_dtau_ddelta(this, T, rhomolar, Wrt, dWrt_dtau, dWrt_ddelta); get_dtau_ddelta(this, T, rhomolar, Constant, dConstant_dtau, dConstant_ddelta); @@ -1119,7 +1119,7 @@ long double HelmholtzEOSMixtureBackend::calc_pressure_nocache(long double T, lon SimpleState reducing = calc_reducing_state_nocache(mole_fractions); long double delta = rhomolar/reducing.rhomolar; long double tau = reducing.T/T; - + // Calculate derivative if needed int nTau = 0, nDelta = 1; long double dalphar_dDelta = calc_alphar_deriv_nocache(nTau, nDelta, mole_fractions, tau, delta); @@ -1143,10 +1143,10 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do long double T, value, r, eos, rhomolar; HelmholtzEOSMixtureBackend *HEOS; - solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double value, int other){ + solver_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double value, int other){ this->HEOS = HEOS; this->T = T; this->value = value; this->other = other; }; - double call(double rhomolar){ + double call(double rhomolar){ this->rhomolar = rhomolar; switch(other) { @@ -1159,7 +1159,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do default: throw ValueError(format("Input not supported")); } - + r = eos-value; return r; }; @@ -1176,7 +1176,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do switch(other) { - + case iSmolar: { ymelt = calc_smolar_nocache(_T, rhomelt); @@ -1204,7 +1204,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do default: throw ValueError(); } - + if (is_in_closed_range(ymelt, yc, y)) { long double rhomolar = Brent(resid, rhomelt, rhoc, LDBL_EPSILON, 1e-12, 100, errstring); @@ -1216,7 +1216,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do return rhomolar; } else - { + { throw ValueError(); } } @@ -1270,7 +1270,7 @@ long double HelmholtzEOSMixtureBackend::solver_for_rho_given_T_oneof_HSU(long do long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double p, long double rhomolar_guess) { int phase; - + // Define the residual to be driven to zero class solver_TP_resid : public FuncWrapper1D { @@ -1278,11 +1278,11 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double long double T, p, r, peos, rhomolar, rhor, tau, R_u, delta, dalphar_dDelta; HelmholtzEOSMixtureBackend *HEOS; - solver_TP_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double p){ - this->HEOS = HEOS; this->T = T; this->p = p; this->rhor = HEOS->get_reducing().rhomolar; + solver_TP_resid(HelmholtzEOSMixtureBackend *HEOS, long double T, long double p){ + this->HEOS = HEOS; this->T = T; this->p = p; this->rhor = HEOS->get_reducing().rhomolar; this->tau = HEOS->get_reducing().T/T; this->R_u = HEOS->gas_constant(); }; - double call(double rhomolar){ + double call(double rhomolar){ this->rhomolar = rhomolar; delta = rhomolar/rhor; dalphar_dDelta = HEOS->calc_alphar_deriv_nocache(0, 1, HEOS->get_mole_fractions(), tau, delta); @@ -1306,7 +1306,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp(long double T, long double if (rhomolar_guess < 0) // Not provided { rhomolar_guess = solver_rho_Tp_SRK(T, p, phase); - + if (phase == iphase_gas && rhomolar_guess < 0)// If the guess is bad, probably high temperature, use ideal gas { rhomolar_guess = p/(gas_constant()*T); @@ -1363,7 +1363,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp_SRK(long double T, long do long double m_j = 0.480+1.574*accentric_j-0.176*pow(accentric_j, 2); long double a_j = 0.42747*pow(R_u*Tcj,2)/pcj*pow(1+m_j*(1-sqrt(T/Tcj)),2); - + if (i == j){ k_ij = 0; } @@ -1390,12 +1390,12 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp_SRK(long double T, long do long double rhomolar0 = p/(Z0*R_u*T); long double rhomolar1 = p/(Z1*R_u*T); long double rhomolar2 = p/(Z2*R_u*T); - + // Check if only one solution is positive, return the solution if that is the case if (rhomolar0 > 0 && rhomolar1 <= 0 && rhomolar2 <= 0){ return rhomolar0; } if (rhomolar0 <= 0 && rhomolar1 > 0 && rhomolar2 <= 0){ return rhomolar1; } - if (rhomolar0 <= 0 && rhomolar1 <= 0 && rhomolar2 > 0){ return rhomolar2; } - + if (rhomolar0 <= 0 && rhomolar1 <= 0 && rhomolar2 > 0){ return rhomolar2; } + switch(phase) { case iphase_liquid: @@ -1410,7 +1410,7 @@ long double HelmholtzEOSMixtureBackend::solver_rho_Tp_SRK(long double T, long do } long double HelmholtzEOSMixtureBackend::calc_pressure(void) -{ +{ // Calculate the reducing parameters _delta = _rhomolar/_reducing.rhomolar; _tau = _reducing.T/_T; @@ -1606,7 +1606,7 @@ SimpleState HelmholtzEOSMixtureBackend::calc_reducing_state_nocache(const std::v SimpleState reducing; if (is_pure_or_pseudopure){ reducing = components[0]->pEOS->reduce; - + } else{ reducing.T = Reducing.p->Tr(mole_fractions); @@ -1653,13 +1653,13 @@ long double HelmholtzEOSMixtureBackend::calc_alphar_deriv_nocache(const int nTau else if (nTau == 3 && nDelta == 0){ return components[0]->pEOS->d3alphar_dTau3(tau, delta); } - else + else { throw ValueError(); } } else{ - + std::size_t N = mole_fractions.size(); long double summer = 0; if (nTau == 0 && nDelta == 0){ @@ -1702,13 +1702,13 @@ long double HelmholtzEOSMixtureBackend::calc_alphar_deriv_nocache(const int nTau for (unsigned int i = 0; i < N; ++i){ summer += mole_fractions[i]*components[i]->pEOS->d3alphar_dTau3(tau, delta); } return summer + pExcess.d3alphar_dTau3(tau, delta); }*/ - else + else { throw ValueError(); } } } -long double HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau, const int nDelta, const std::vector &mole_fractions, +long double HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau, const int nDelta, const std::vector &mole_fractions, const long double &tau, const long double &delta, const long double &Tr, const long double &rhor) { if (is_pure_or_pseudopure) @@ -1743,7 +1743,7 @@ long double HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau else if (nTau == 3 && nDelta == 0){ return components[0]->pEOS->d3alpha0_dTau3(tau, delta); } - else + else { throw ValueError(); } @@ -1753,31 +1753,31 @@ long double HelmholtzEOSMixtureBackend::calc_alpha0_deriv_nocache(const int nTau std::size_t N = mole_fractions.size(); long double summer = 0; long double tau_i, delta_i, rho_ci, T_ci; - for (unsigned int i = 0; i < N; ++i){ - rho_ci = components[i]->pEOS->reduce.rhomolar; + for (unsigned int i = 0; i < N; ++i){ + rho_ci = components[i]->pEOS->reduce.rhomolar; T_ci = components[i]->pEOS->reduce.T; tau_i = T_ci*tau/Tr; delta_i = delta*rhor/rho_ci; - if (nTau == 0 && nDelta == 0){ - summer += mole_fractions[i]*(components[i]->pEOS->base0(tau_i, delta_i)+log(mole_fractions[i])); + if (nTau == 0 && nDelta == 0){ + summer += mole_fractions[i]*(components[i]->pEOS->base0(tau_i, delta_i)+log(mole_fractions[i])); } else if (nTau == 0 && nDelta == 1){ - summer += mole_fractions[i]*rhor/rho_ci*components[i]->pEOS->dalpha0_dDelta(tau_i, delta_i); + summer += mole_fractions[i]*rhor/rho_ci*components[i]->pEOS->dalpha0_dDelta(tau_i, delta_i); } else if (nTau == 1 && nDelta == 0){ - summer += mole_fractions[i]*T_ci/Tr*components[i]->pEOS->dalpha0_dTau(tau_i, delta_i); + summer += mole_fractions[i]*T_ci/Tr*components[i]->pEOS->dalpha0_dTau(tau_i, delta_i); } else if (nTau == 0 && nDelta == 2){ - summer += mole_fractions[i]*pow(rhor/rho_ci,2)*components[i]->pEOS->d2alpha0_dDelta2(tau_i, delta_i); + summer += mole_fractions[i]*pow(rhor/rho_ci,2)*components[i]->pEOS->d2alpha0_dDelta2(tau_i, delta_i); } else if (nTau == 1 && nDelta == 1){ - summer += mole_fractions[i]*rhor/rho_ci*T_ci/Tr*components[i]->pEOS->d2alpha0_dDelta_dTau(tau_i, delta_i); + summer += mole_fractions[i]*rhor/rho_ci*T_ci/Tr*components[i]->pEOS->d2alpha0_dDelta_dTau(tau_i, delta_i); } else if (nTau == 2 && nDelta == 0){ - summer += mole_fractions[i]*pow(T_ci/Tr,2)*components[i]->pEOS->d2alpha0_dTau2(tau_i, delta_i); + summer += mole_fractions[i]*pow(T_ci/Tr,2)*components[i]->pEOS->d2alpha0_dTau2(tau_i, delta_i); } - else + else { throw ValueError(); } @@ -1940,7 +1940,7 @@ long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dxj__constT_V_xi( { // 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_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) @@ -2025,7 +2025,7 @@ long double HelmholtzEOSMixtureBackend::mixderiv_d_ndalphardni_dxj__constdelta_t 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); diff --git a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h index eba921a1..5a47040e 100644 --- a/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h +++ b/src/Backends/Helmholtz/HelmholtzEOSMixtureBackend.h @@ -28,7 +28,7 @@ protected: SimpleState _crit; int imposed_phase_index; public: - HelmholtzEOSMixtureBackend(){SatL = NULL; SatV = NULL; imposed_phase_index = -1;}; + HelmholtzEOSMixtureBackend(){imposed_phase_index = -1;}; HelmholtzEOSMixtureBackend(std::vector components, bool generate_SatL_and_SatV = true); HelmholtzEOSMixtureBackend(std::vector &component_names, bool generate_SatL_and_SatV = true); virtual ~HelmholtzEOSMixtureBackend(){}; diff --git a/src/CoolProp.cpp b/src/CoolProp.cpp index ba09aedc..989656dd 100644 --- a/src/CoolProp.cpp +++ b/src/CoolProp.cpp @@ -16,6 +16,8 @@ #endif #endif +#include + #include #include #include @@ -27,6 +29,7 @@ #include "Backends/Helmholtz/Fluids/FluidLibrary.h" #include "Backends/Helmholtz/HelmholtzEOSBackend.h" + namespace CoolProp { diff --git a/src/HumidAirProp.cpp b/src/HumidAirProp.cpp index 16242419..c2f9fad7 100644 --- a/src/HumidAirProp.cpp +++ b/src/HumidAirProp.cpp @@ -2,6 +2,8 @@ #define _CRT_SECURE_NO_WARNINGS #endif +#include + #include "HumidAirProp.h" #include "AbstractState.h" #include "Solvers.h" @@ -9,6 +11,8 @@ #include "Ice.h" #include "CoolProp.h" + + #include #include "math.h" #include "time.h" diff --git a/src/SpeedTest.cpp b/src/SpeedTest.cpp index c37499e1..b6a192c3 100644 --- a/src/SpeedTest.cpp +++ b/src/SpeedTest.cpp @@ -1,3 +1,4 @@ +#include #include "SpeedTest.h" #include "AbstractState.h" #include "DataStructures.h" @@ -10,7 +11,7 @@ void compare_REFPROP_and_CoolProp(std::string fluid, int inputs, double val1, do { time_t t1,t2; double dx = 1/((double)N); - + std::tr1::shared_ptr State(AbstractState::factory("HEOS", fluid)); t1 = clock(); for (std::size_t ii = 0; ii < N; ++ii) @@ -18,7 +19,7 @@ void compare_REFPROP_and_CoolProp(std::string fluid, int inputs, double val1, do State->update(inputs, val1 + ii*d1, val2 + ii*d2); } t2 = clock(); - + double elap = ((double)(t2-t1))/CLOCKS_PER_SEC/((double)N)*1e6; printf("Elapsed time for CoolProp is %g us/call\n",elap); @@ -33,4 +34,4 @@ void compare_REFPROP_and_CoolProp(std::string fluid, int inputs, double val1, do printf("Elapsed time for REFPROP is %g us/call\n",elap); } -} /* namespace CoolProp */ \ No newline at end of file +} /* namespace CoolProp */ diff --git a/src/main.cxx b/src/main.cxx index 389fd7d8..046efb98 100644 --- a/src/main.cxx +++ b/src/main.cxx @@ -17,12 +17,12 @@ using namespace CoolProp; #endif #include "SpeedTest.h" -#include +//#include void generate_melting_curve_data(const char* file_name, const char *fluid_name, double Tmin, double Tmax) { - + FILE *fp; fp = fopen(file_name,"w"); AbstractState *State = AbstractState::factory(std::string("REFPROP"),std::string(fluid_name)); @@ -43,7 +43,7 @@ void generate_melting_curve_data(const char* file_name, const char *fluid_name, } catch(std::exception &e) { - + std::cout << fluid_name << " " << e.what() << std::endl; break; } @@ -51,8 +51,13 @@ void generate_melting_curve_data(const char* file_name, const char *fluid_name, fclose(fp); delete State; } +struct element + { + double d,t,ld; + int l; + }; int main() -{ +{ if (0) { generate_melting_curve_data("Ethylene-I.mlt","ethylene",103.989,110.369); @@ -61,7 +66,7 @@ int main() generate_melting_curve_data("Propylene-II.mlt","propylen",109.6,575); generate_melting_curve_data("ParaHydrogen-I.mlt","parahyd",13.8033,22); generate_melting_curve_data("ParaHydrogen-II.mlt","parahyd",22,2000); - + generate_melting_curve_data("n-Propane.mlt","propane",85.53,2000); generate_melting_curve_data("n-Butane.mlt","butane",134.9,2000); generate_melting_curve_data("n-Pentane.mlt","pentane",143.5,2000); @@ -103,12 +108,12 @@ int main() if (0) { std::vector ss = strsplit(get_global_param_string("FluidsList"),','); - + for (std::vector::iterator it = ss.begin(); it != ss.end(); ++it) { AbstractState *S = AbstractState::factory("HEOS", (*it)); S->update(QT_INPUTS, 0, S->Ttriple()); - std::cout << format("%s %17.15g\n", S->name(), S->p()); + std::cout << format("%s %17.15g\n", S->name().c_str(), S->p()); } } if (1) @@ -137,7 +142,7 @@ int main() double p2 = AS->umolar(); double d2 = AS->rhomolar(); double T2 = AS->delta(); - + double dpdT_constrho2 = (p2-p1)/(T2-T1); AS->update(PT_INPUTS, 101000, 300); @@ -179,13 +184,9 @@ int main() } if (0) { - struct element - { - double d,t,ld; - int l; - }; - double n[] = {0.0125335479355233, 7.8957634722828, -8.7803203303561, 0.31802509345418, -0.26145533859358, -0.0078199751687981, 0.0088089493102134, -0.66856572307965, 0.20433810950965, -6.621260503968699e-005, -0.19232721156002, -0.25709043003438, 0.16074868486251, -0.040092828925807, 3.9343422603254e-007, -7.5941377088144e-006, 0.00056250979351888, -1.5608652257135e-005, 1.1537996422951e-009, 3.6582165144204e-007, -1.3251180074668e-012, -6.2639586912454e-010, -0.10793600908932, 0.017611491008752, 0.22132295167546, -0.40247669763528, 0.58083399985759, 0.0049969146990806, -0.031358700712549, -0.74315929710341, 0.4780732991548, 0.020527940895948, -0.13636435110343, 0.014180634400617, 0.008332650488071301, -0.029052336009585, 0.038615085574206, -0.020393486513704, -0.0016554050063734, 0.0019955571979541, 0.00015870308324157, -1.638856834253e-005, 0.043613615723811, 0.034994005463765, -0.076788197844621, 0.022446277332006, -6.2689710414685e-005, -5.5711118565645e-010, -0.19905718354408, 0.31777497330738, -0.11841182425981 }; - double d[] = {1, 1, 1, 2, 2, 3, 4, 1, 1, 1, 2, 2, 3, 4, 4, 5, 7, 9, 10, 11, 13, 15, 1, 2, 2, 2, 3, 4, 4, 4, 5, 6, 6, 7, 9, 9, 9, 9, 9, 10, 10, 12, 3, 4, 4, 5, 14, 3, 6, 6, 6 }; + + double n[] = {0.0125335479355233, 7.8957634722828, -8.7803203303561, 0.31802509345418, -0.26145533859358, -0.0078199751687981, 0.0088089493102134, -0.66856572307965, 0.20433810950965, -6.621260503968699e-005, -0.19232721156002, -0.25709043003438, 0.16074868486251, -0.040092828925807, 3.9343422603254e-007, -7.5941377088144e-006, 0.00056250979351888, -1.5608652257135e-005, 1.1537996422951e-009, 3.6582165144204e-007, -1.3251180074668e-012, -6.2639586912454e-010, -0.10793600908932, 0.017611491008752, 0.22132295167546, -0.40247669763528, 0.58083399985759, 0.0049969146990806, -0.031358700712549, -0.74315929710341, 0.4780732991548, 0.020527940895948, -0.13636435110343, 0.014180634400617, 0.008332650488071301, -0.029052336009585, 0.038615085574206, -0.020393486513704, -0.0016554050063734, 0.0019955571979541, 0.00015870308324157, -1.638856834253e-005, 0.043613615723811, 0.034994005463765, -0.076788197844621, 0.022446277332006, -6.2689710414685e-005, -5.5711118565645e-010, -0.19905718354408, 0.31777497330738, -0.11841182425981 }; + double d[] = {1, 1, 1, 2, 2, 3, 4, 1, 1, 1, 2, 2, 3, 4, 4, 5, 7, 9, 10, 11, 13, 15, 1, 2, 2, 2, 3, 4, 4, 4, 5, 6, 6, 7, 9, 9, 9, 9, 9, 10, 10, 12, 3, 4, 4, 5, 14, 3, 6, 6, 6 }; double t[] = {-0.5, 0.875, 1, 0.5, 0.75, 0.375, 1, 4, 6, 12, 1, 5, 4, 2, 13, 9, 3, 4, 11, 4, 13, 1, 7, 1, 9, 10, 10, 3, 7, 10, 10, 6, 10, 10, 1, 2, 3, 4, 8, 6, 9, 8, 16, 22, 23, 23, 10, 50, 44, 46, 50 }; double l[] = {0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 6, 6, 6, 6 }; double summer = 0; @@ -216,14 +217,14 @@ int main() summer += (di-lid*pow_delta_li)*exp(ti*log_tau+(di-1)*log_delta-pow_delta_li); } else{ - summer += di*exp(ti*log_tau+(di-1)*log_delta); + summer += di*exp(ti*log_tau+(di-1)*log_delta); } } } double t2 = clock(); double elap = (t2-t1)/CLOCKS_PER_SEC/((double)N)*1e6; printf("%g %g\n",elap, summer); - + } if (0) { @@ -240,7 +241,7 @@ int main() { double T = 300; - + AbstractState *MixRP = AbstractState::factory(std::string("REFPROP"), std::string("propane")); { long N = 100000; @@ -324,9 +325,9 @@ int main() AbstractState *Mix = AbstractState::factory(std::string("CORE"),std::string("Ethane,n-Propane")); Mix->set_mole_fractions(z); - + for (double T = 210; ;T += 0.1) - { + { Mix->update(QT_INPUTS, Q, T); std::cout << format(" %g %g\n",Mix->p(),Mix->T()); } @@ -335,7 +336,7 @@ int main() if(0) { time_t t1,t2; - + std::size_t N = 1000000; AbstractState *State = AbstractState::factory(std::string("CORE"), std::string("Water")); double p = State->p(); @@ -363,7 +364,7 @@ int main() } - + if (0) { @@ -414,5 +415,5 @@ int main() //double sigma = State->surface_tension(); delete State; - } -} \ No newline at end of file + } +}