Added melting lines for many fluids - nearly all that are in REFPROP

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-06-09 15:48:52 +02:00
parent 2b32a18e13
commit 64946e6017
28 changed files with 809 additions and 112 deletions

View File

@@ -311,7 +311,7 @@ protected:
EOS.reduce.rhomolar = cpjson::get_double(reducing_state,"rhomolar");
EOS.reduce.p = cpjson::get_double(reducing_state,"p");
/// todo: define limits of EOS better
/// \todo: define limits of EOS better
EOS.limits.Tmin = cpjson::get_double(satminL_state, "T");
EOS.Ttriple = EOS.limits.Tmin;
@@ -804,6 +804,71 @@ protected:
fluid.transport.epsilon_over_k = Tc/1.3593; // [K]
}
void parse_melting_line(rapidjson::Value &melting_line, CoolPropFluid & fluid)
{
fluid.ancillaries.melting_line.T_m = cpjson::get_double(melting_line, "T_m");
fluid.ancillaries.melting_line.BibTeX = cpjson::get_string(melting_line, "BibTeX");
if (melting_line.HasMember("type"))
{
std::string type = cpjson::get_string(melting_line, "type");
if (!type.compare("Simon"))
{
rapidjson::Value &parts = melting_line["parts"];
for (rapidjson::Value::ValueIterator itr = parts.Begin(); itr != parts.End(); ++itr)
{
MeltingLinePiecewiseSimonSegment data;
data.a = cpjson::get_double((*itr),"a");
data.c = cpjson::get_double((*itr),"c");
data.T_min = cpjson::get_double((*itr),"T_min");
data.T_max = cpjson::get_double((*itr),"T_max");
data.T_0 = cpjson::get_double((*itr),"T_0");
data.p_0 = cpjson::get_double((*itr),"p_0");
fluid.ancillaries.melting_line.simon.parts.push_back(data);
}
fluid.ancillaries.melting_line.type = MeltingLineVariables::MELTING_LINE_SIMON_TYPE;
}
else if (!type.compare("polynomial_in_Tr"))
{
rapidjson::Value &parts = melting_line["parts"];
for (rapidjson::Value::ValueIterator itr = parts.Begin(); itr != parts.End(); ++itr)
{
MeltingLinePiecewisePolynomialInTrSegment data;
data.a = cpjson::get_long_double_array((*itr),"a");
data.t = cpjson::get_long_double_array((*itr),"t");
data.T_min = cpjson::get_double((*itr),"T_min");
data.T_max = cpjson::get_double((*itr),"T_max");
data.T_0 = cpjson::get_double((*itr),"T_0");
data.p_0 = cpjson::get_double((*itr),"p_0");
fluid.ancillaries.melting_line.polynomial_in_Tr.parts.push_back(data);
}
fluid.ancillaries.melting_line.type = MeltingLineVariables::MELTING_LINE_POLYNOMIAL_IN_TR_TYPE;
}
else if (!type.compare("polynomial_in_Theta"))
{
rapidjson::Value &parts = melting_line["parts"];
for (rapidjson::Value::ValueIterator itr = parts.Begin(); itr != parts.End(); ++itr)
{
MeltingLinePiecewisePolynomialInThetaSegment data;
data.a = cpjson::get_long_double_array((*itr),"a");
data.t = cpjson::get_long_double_array((*itr),"t");
data.T_min = cpjson::get_double((*itr),"T_min");
data.T_max = cpjson::get_double((*itr),"T_max");
data.T_0 = cpjson::get_double((*itr),"T_0");
data.p_0 = cpjson::get_double((*itr),"p_0");
fluid.ancillaries.melting_line.polynomial_in_Theta.parts.push_back(data);
}
fluid.ancillaries.melting_line.type = MeltingLineVariables::MELTING_LINE_POLYNOMIAL_IN_THETA_TYPE;
}
else{
throw ValueError(format("melting line type [%s] is not understood for fluid %s", type.c_str(), fluid.name.c_str()));
}
}
else{
throw ValueError(format("melting line does not have \"type\" for fluid %s", fluid.name.c_str()));
}
};
/// Parse the critical state for the given EOS
void parse_states(rapidjson::Value &states, CoolPropFluid & fluid)
{
@@ -936,6 +1001,16 @@ public:
parse_surface_tension(fluid_json["ANCILLARIES"]["surface_tension"], fluid);
}
// Melting line
if (!(fluid_json["ANCILLARIES"].HasMember("melting_line"))){
if (get_debug_level() > 0){
std::cout << format("Melting line curves are missing for fluid [%s]\n", fluid.name.c_str()) ;
}
}
else{
parse_melting_line(fluid_json["ANCILLARIES"]["melting_line"], fluid);
}
// Parse the environmental parameters
if (!(fluid_json.HasMember("ENVIRONMENTAL"))){
if (get_debug_level() > 0){

View File

@@ -115,6 +115,17 @@ long double HelmholtzEOSMixtureBackend::calc_molar_mass(void)
}
return summer;
}
long double HelmholtzEOSMixtureBackend::calc_melt_p_T(long double T)
{
if (is_pure_or_pseudopure)
{
return components[0]->ancillaries.melting_line.evaluate(iP, iT, T);
}
else
{
throw NotImplementedError(format("calc_melt_p_T not implemented for mixtures"));
}
}
long double HelmholtzEOSMixtureBackend::calc_surface_tension(void)
{
if (is_pure_or_pseudopure)
@@ -1591,8 +1602,8 @@ long double HelmholtzEOSMixtureBackend::calc_speed_sound(void)
long double d2ar_dDelta2 = d2alphar_dDelta2();
long double d2ar_dDelta_dTau = d2alphar_dDelta_dTau();
long double d2ar_dTau2 = d2alphar_dTau2();
long double R_u = static_cast<long double>(_gas_constant);
long double mm = static_cast<long double>(_molar_mass);
long double R_u = gas_constant();
long double mm = molar_mass();
// Get speed of sound
_speed_sound = sqrt(R_u*_T/mm*(1+2*_delta.pt()*dar_dDelta+pow(_delta.pt(),2)*d2ar_dDelta2 - pow(1+_delta.pt()*dar_dDelta-_delta.pt()*_tau.pt()*d2ar_dDelta_dTau,2)/(pow(_tau.pt(),2)*(d2ar_dTau2 + d2a0_dTau2))));

View File

@@ -17,10 +17,10 @@ class FlashRoutines;
class HelmholtzEOSMixtureBackend : public AbstractState {
protected:
std::vector<CoolPropFluid*> components; ///< The components that are in use
bool is_pure_or_pseudopure; ///< A flag for whether the substance is a pure or pseudo-pure fluid (true) or a mixture (false)
std::vector<long double> mole_fractions; ///< The mole fractions of the components
std::vector<long double> mole_fractions_liq, ///< The mole fractions of the saturated liquid
std::vector<long double> mole_fractions_liq, ///< The mole fractions of the saturated liquid
mole_fractions_vap, ///< The mole fractions of the saturated vapor
K, ///< The K factors for the components
lnK; ///< The natural logarithms of the K factors of the components
@@ -37,7 +37,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; //
// Helmholtz EOS backend uses mole fractions
bool using_mole_fractions(){return true;}
@@ -46,7 +46,7 @@ public:
std::vector<long double> &get_K(){return K;};
std::vector<long double> &get_lnK(){return lnK;};
shared_ptr<HelmholtzEOSMixtureBackend> SatL, SatV; ///<
shared_ptr<HelmholtzEOSMixtureBackend> SatL, SatV; ///<
void update(long input_pair, double value1, double value2);
@@ -69,19 +69,19 @@ public:
void set_excess_term();
/// Set the mole fractions
/**
/**
@param mole_fractions The vector of mole fractions of the components
*/
void set_mole_fractions(const std::vector<long double> &mole_fractions);
const std::vector<long double> &get_mole_fractions(){return mole_fractions;};
/// Set the mass fractions
/**
/**
@param mass_fractions The vector of mass fractions of the components
*/
void set_mass_fractions(const std::vector<long double> &mass_fractions){throw std::exception();};
long double calc_molar_mass(void);
long double calc_gas_constant(void);
@@ -118,6 +118,7 @@ public:
long double calc_d2alpha0_dDelta_dTau(void);
long double calc_d2alpha0_dTau2(void);
long double calc_melt_p_T(long double T);
long double calc_surface_tension(void);
long double calc_viscosity(void);
long double calc_viscosity_dilute(void);
@@ -137,7 +138,7 @@ public:
std::string calc_name(void);
long double calc_alphar_deriv_nocache(const int nTau, const int nDelta, const std::vector<long double> & mole_fractions, const long double &tau, const long double &delta);
/**
\brief Take derivatives of the ideal-gas part of the Helmholtz energy, don't use any cached values, or store any cached values
@@ -163,7 +164,7 @@ public:
\sa Table B5, GERG 2008 from Kunz Wagner, JCED, 2012
*/
long double calc_alpha0_deriv_nocache(const int nTau, const int nDelta, const std::vector<long double> & mole_fractions, const long double &tau, const long double &delta, const long double &Tr, const long double &rhor);
void calc_reducing_state(void);
SimpleState calc_reducing_state_nocache(const std::vector<long double> & mole_fractions);
@@ -183,11 +184,11 @@ public:
void mass_to_molar_inputs(long &input_pair, double &value1, double &value2);
// ***************************************************************
// ***************************************************************
// ***************************************************************
// ************* PHASE DETERMINATION ROUTINES ******************
// ***************************************************************
// ***************************************************************
// ***************************************************************
void T_phase_determination_pure_or_pseudopure(int other, long double value);
void p_phase_determination_pure_or_pseudopure(int other, long double value);
void DmolarP_phase_determination();
@@ -197,8 +198,8 @@ public:
// ***************************************************************
// ******************* SOLVER ROUTINES *************************
// ***************************************************************
// ***************************************************************
// ***************************************************************
long double solver_rho_Tp(long double T, long double p, long double rho_guess = -1);
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);
@@ -310,7 +311,7 @@ public:
\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
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]
@@ -320,7 +321,7 @@ public:
/// Gernert Equation 3.115
/// Catch test provided
long double mixderiv_dln_fugacity_coefficient_dxj__constT_p_xi(int i, int j);
/// Gernert Equation 3.130
/// Catch test provided
long double mixderiv_dpdxj__constT_V_xi(int j);
@@ -390,7 +391,7 @@ public:
\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}\\
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}
@@ -404,7 +405,7 @@ public:
\f]
*/
long double mixderiv_nddeltadni__constT_V_nj(int i);
/// GERG 2004 Monograph equation 7.49
/** The derivative term
\f[