Implemented viscosity for cyclohexane from Tariq, 2014.

Added unit tests, reformulated initial density term

Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
This commit is contained in:
Ian Bell
2014-08-21 00:50:57 +02:00
parent 3c3ffa66c4
commit b4ff7a8cc5
8 changed files with 305 additions and 28 deletions

View File

@@ -378,8 +378,8 @@ protected:
if (!target.compare("Ethane")){
fluid.transport.viscosity_dilute.type = CoolProp::ViscosityDiluteVariables::VISCOSITY_DILUTE_ETHANE; return;
}
else if (!target.compare("Ethane")){
fluid.transport.viscosity_dilute.type = CoolProp::ViscosityDiluteVariables::VISCOSITY_DILUTE_ETHANE; return;
else if (!target.compare("Cyclohexane")){
fluid.transport.viscosity_dilute.type = CoolProp::ViscosityDiluteVariables::VISCOSITY_DILUTE_CYCLOHEXANE; return;
}
else{
throw ValueError(format("hardcoded dilute viscosity [%s] is not understood for fluid %s",target.c_str(),fluid.name.c_str()));
@@ -430,16 +430,33 @@ protected:
};
/// Parse the transport properties
void parse_initial_density_viscosity(rapidjson::Value &dilute, CoolPropFluid & fluid)
void parse_initial_density_viscosity(rapidjson::Value &initial_density, CoolPropFluid & fluid)
{
std::string type = cpjson::get_string(dilute, "type");
std::string type = cpjson::get_string(initial_density, "type");
if (!type.compare("Rainwater-Friend")){
// Get a reference to the entry in the fluid instance
CoolProp::ViscosityRainWaterFriendData &RF = fluid.transport.viscosity_initial.rainwater_friend;
// Load up the values
RF.b = cpjson::get_long_double_array(dilute["b"]);
RF.t = cpjson::get_long_double_array(dilute["t"]);
RF.b = cpjson::get_long_double_array(initial_density["b"]);
RF.t = cpjson::get_long_double_array(initial_density["t"]);
// Set the type flag
fluid.transport.viscosity_initial.type = CoolProp::ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_RAINWATER_FRIEND;
}
else if (!type.compare("empirical")){
// Get a reference to the entry in the fluid instance
CoolProp::ViscosityInitialDensityEmpiricalData &EM = fluid.transport.viscosity_initial.empirical;
// Load up the values
EM.n = cpjson::get_long_double_array(initial_density["n"]);
EM.d = cpjson::get_long_double_array(initial_density["d"]);
EM.t = cpjson::get_long_double_array(initial_density["t"]);
EM.T_reducing = cpjson::get_double(initial_density,"T_reducing");
EM.rhomolar_reducing = cpjson::get_double(initial_density,"rhomolar_reducing");
// Set the type flag
fluid.transport.viscosity_initial.type = CoolProp::ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_EMPIRICAL;
}
else{
throw ValueError(format("type [%s] is not understood for fluid %s",type.c_str(),fluid.name.c_str()));

View File

@@ -202,6 +202,8 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity_dilute(void)
eta_dilute = TransportRoutines::viscosity_dilute_collision_integral_powers_of_T(*this); break;
case ViscosityDiluteVariables::VISCOSITY_DILUTE_ETHANE:
eta_dilute = TransportRoutines::viscosity_dilute_ethane(*this); break;
case ViscosityDiluteVariables::VISCOSITY_DILUTE_CYCLOHEXANE:
eta_dilute = TransportRoutines::viscosity_dilute_cyclohexane(*this); break;
default:
throw ValueError(format("dilute viscosity type [%d] is invalid for fluid %s", components[0]->transport.viscosity_dilute.type, name().c_str()));
}
@@ -220,13 +222,29 @@ long double HelmholtzEOSMixtureBackend::calc_viscosity_background()
}
long double HelmholtzEOSMixtureBackend::calc_viscosity_background(long double eta_dilute)
{
// Residual part
long double B_eta_initial = TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(*this);
long double rho = rhomolar();
long double initial_part = eta_dilute*B_eta_initial*rho;
long double initial_part = 0.0;
switch(components[0]->transport.viscosity_initial.type){
case ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_RAINWATER_FRIEND:
{
long double B_eta_initial = TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(*this);
long double rho = rhomolar();
initial_part = eta_dilute*B_eta_initial*rho;
break;
}
case ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_EMPIRICAL:
{
initial_part = TransportRoutines::viscosity_initial_density_dependence_empirical(*this);
break;
}
case ViscosityInitialDensityVariables::VISCOSITY_INITIAL_DENSITY_NOT_SET:
{
break;
}
}
// Higher order terms
long double delta_eta_h;
long double delta_eta_h = 0.0;
switch(components[0]->transport.viscosity_higher_order.type)
{
case ViscosityHigherOrderVariables::VISCOSITY_HIGHER_ORDER_BATSCHINKI_HILDEBRAND:
@@ -592,7 +610,9 @@ void HelmholtzEOSMixtureBackend::mass_to_molar_inputs(CoolProp::input_pairs &inp
case DmassHmass_INPUTS: input_pair = DmolarHmolar_INPUTS; value1 /= mm; value2 *= mm; break;
case DmassSmass_INPUTS: input_pair = DmolarSmolar_INPUTS; value1 /= mm; value2 *= mm; break;
case DmassUmass_INPUTS: input_pair = DmolarUmolar_INPUTS; value1 /= mm; value2 *= mm; break;
default: break;
}
break;
}
default:
return;
@@ -625,6 +645,7 @@ void HelmholtzEOSMixtureBackend::update(CoolProp::input_pairs input_pair, double
long double ld_value1 = value1, ld_value2 = value2;
pre_update(input_pair, ld_value1, ld_value2);
value1 = ld_value1; value2 = ld_value2;
switch(input_pair)
{

View File

@@ -162,6 +162,29 @@ long double TransportRoutines::viscosity_initial_density_dependence_Rainwater_Fr
}
}
long double TransportRoutines::viscosity_initial_density_dependence_empirical(HelmholtzEOSMixtureBackend &HEOS)
{
// Inspired by the form from Tariq, JPCRD, 2
if (HEOS.is_pure_or_pseudopure)
{
// Retrieve values from the state class
CoolProp::ViscosityInitialDensityEmpiricalData &data = HEOS.components[0]->transport.viscosity_initial.empirical;
const std::vector<long double> &n = data.n, &d = data.d, &t = data.t;
long double tau = data.T_reducing/HEOS.T(); // [no units]
long double delta = HEOS.rhomolar()/data.rhomolar_reducing; // [no units]
long double summer = 0;
for (unsigned int i = 0; i < n.size(); ++i){
summer += n[i]*pow(delta, d[i])*pow(tau, t[i]);
}
return summer; // [Pa-s]
}
else{
throw NotImplementedError("TransportRoutines::viscosity_initial_density_dependence_empirical is only for pure and pseudo-pure");
}
}
static void visc_Helper(double Tbar, double rhobar, double *mubar_0, double *mubar_1)
{
std::vector<std::vector<long double> > H(6,std::vector<long double>(7,0));
@@ -436,6 +459,13 @@ long double TransportRoutines::viscosity_dilute_ethane(HelmholtzEOSMixtureBacken
return 12.0085*sqrt(Tstar)*OMEGA_2_2/1e6; //[Pa-s]
}
long double TransportRoutines::viscosity_dilute_cyclohexane(HelmholtzEOSMixtureBackend &HEOS)
{
// From Tariq, JPCRD, 2014
long double T = HEOS.T();
long double S_eta = exp(-1.5093 + 364.87/T - 39537/pow(T, 2)); //[nm^2]
return 0.19592*sqrt(T)/S_eta/1e6; //[Pa-s]
}
long double TransportRoutines::viscosity_ethane_higher_order_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
double r[] = {0,1,1,2,2,2,3,3,4,4,1,1};
@@ -531,7 +561,7 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(
rhoc = HEOS.get_reducing_state().rhomolar, // [mol/m^3]
Pcrit = HEOS.get_reducing_state().p, // [Pa]
Tref, // [K]
cp,cv,delta,num,zeta,mu,pi=M_PI,tau,
cp,cv,delta,num,zeta,mu,pi=M_PI,
OMEGA_tilde,OMEGA_tilde0;
if (ValidNumber(data.T_ref))
@@ -541,7 +571,6 @@ long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(
delta = HEOS.delta();
tau = HEOS.tau();
double dp_drho = HEOS.gas_constant()*HEOS.T()*(1+2*delta*HEOS.dalphar_dDelta()+delta*delta*HEOS.d2alphar_dDelta2());
double X = Pcrit/pow(rhoc,2)*HEOS.rhomolar()/dp_drho;
@@ -660,7 +689,7 @@ long double TransportRoutines::conductivity_hardcoded_water(HelmholtzEOSMixtureB
double lambdabar_0,lambdabar_1,lambdabar_2,rhobar,Tbar,sum;
double Tstar=647.096,rhostar=322,pstar=22064000,lambdastar=1e-3,mustar=1e-6;
double tau,xi;
double xi;
int i,j;
Tbar = HEOS.T()/Tstar;
@@ -680,7 +709,6 @@ long double TransportRoutines::conductivity_hardcoded_water(HelmholtzEOSMixtureB
double nu=0.630,GAMMA =177.8514,gamma=1.239,xi_0=0.13,Lambda_0=0.06,Tr_bar=1.5,
qd_bar=1/0.4,pi=3.141592654, delta = HEOS.delta(), R=461.51805;//J/kg/K
tau=1/Tbar;
double drhodp = 1/(R*HEOS.T()*(1+2*rhobar*HEOS.dalphar_dDelta()+rhobar*rhobar*HEOS.d2alphar_dDelta2()));
double drhobar_dpbar = pstar/rhostar*drhodp;

View File

@@ -67,12 +67,23 @@ public:
IMPORTANT: This function returns \f$B_{\eta}\f$, not \f$\eta_{RF}\f$
*/
static long double viscosity_initial_density_dependence_Rainwater_Friend(HelmholtzEOSMixtureBackend &HEOS);
/**
* \brief An empirical form for the initial density dependence
*
* Given by the polynomial-like form
* \f[
* \eta^1 = \sum_i n_i\delta^{d_i}\tau^{t_i}
* \f]
* where the output is in Pa-s
*/
static long double viscosity_initial_density_dependence_empirical(HelmholtzEOSMixtureBackend &HEOS);
/**
\brief The modified Batschinski-Hildebrand contribution to the viscosity
\f[
\Delta\eta = \displaystyle\sum_{i}a_{i}\delta^{d1_i}\tau^{t1_j}+\left(\displaystyle\sum_{i}f_i\delta^{d2_i}\tau^{t2_i}\right)\left(\frac{1}{\delta_0(\tau)-\delta}-\frac{1}{\delta_0(\tau)}\right)
\Delta\eta = \displaystyle\sum_{i}a_{i}\delta^{d1_i}\tau^{t1_j}\exp(\gamma_i\delta^{l_i})+\left(\displaystyle\sum_{i}f_i\delta^{d2_i}\tau^{t2_i}\right)\left(\frac{1}{\delta_0(\tau)-\delta}-\frac{1}{\delta_0(\tau)}\right)
\f]
where \f$\tau = T_c/T\f$ and \f$\delta = \rho/\rho_c\f$
\f[
@@ -83,6 +94,7 @@ public:
static long double viscosity_higher_order_modified_Batschinski_Hildebrand(HelmholtzEOSMixtureBackend &HEOS);
static long double viscosity_dilute_ethane(HelmholtzEOSMixtureBackend &HEOS);
static long double viscosity_dilute_cyclohexane(HelmholtzEOSMixtureBackend &HEOS);
static long double viscosity_water_hardcoded(HelmholtzEOSMixtureBackend &HEOS);
static long double viscosity_helium_hardcoded(HelmholtzEOSMixtureBackend &HEOS);

View File

@@ -206,6 +206,21 @@ vel("Ethane", "T", 100, "Dmolar", 21330, "V", 878.6e-6, 1e-2),
vel("Ethane", "T", 430, "Dmolar", 12780, "V", 58.70e-6, 1e-2),
vel("Ethane", "T", 500, "Dmolar", 11210, "V", 48.34e-6, 1e-2),
// From Xiang, JPCRD, 2006
vel("Methanol", "T", 600, "Dmass", 800.23, "V", 0.1888e-3, 1e-4),
vel("Methanol", "T", 600, "Dmass", 833.20, "V", 0.2092e-3, 1e-4),
vel("Methanol", "T", 600, "Dmass", 861.37, "V", 0.2279e-3, 1e-4),
vel("Methanol", "T", 600, "Dmass", 908.33, "V", 0.2634e-3, 1e-4),
vel("Methanol", "T", 620, "Dmass", 788.58, "V", 0.1779e-3, 1e-4),
vel("Methanol", "T", 620, "Dmass", 822.14, "V", 0.1972e-3, 1e-4),
vel("Methanol", "T", 620, "Dmass", 850.77, "V", 0.2148e-3, 1e-4),
vel("Methanol", "T", 620, "Dmass", 898.48, "V", 0.2477e-3, 1e-4),
vel("Methanol", "T", 630, "Dmass", 782.76, "V", 0.1729e-3, 1e-4),
vel("Methanol", "T", 630, "Dmass", 811.06, "V", 0.1917e-3, 1e-4),
vel("Methanol", "T", 630, "Dmass", 840.11, "V", 0.2088e-3, 1e-4),
vel("Methanol", "T", 630, "Dmass", 888.50, "V", 0.2405e-3, 1e-4),
// From REFPROP 9.1 since no data provided
vel("n-Butane", "T", 150, "Q", 0, "V", 0.0013697657668, 1e-4),
vel("n-Butane", "T", 400, "Q", 1, "V", 1.2027464524762453e-005, 1e-4),
@@ -214,6 +229,17 @@ vel("IsoButane", "T", 400, "Q", 1, "V", 1.4761041187617117e-005, 2e-4),
vel("R134a", "T", 175, "Q", 0, "V", 0.0017558494524138289, 1e-4),
vel("R134a", "T", 360, "Q", 1, "V", 1.7140264998576107e-005, 1e-4),
// From Tariq, JPCRD, 2014
vel("Cyclohexane", "T", 300, "Dmolar", 1e-10, "V", 7.058e-6, 1e-4),
vel("Cyclohexane", "T", 300, "Dmolar", 0.0430e3, "V", 6.977e-6, 1e-4),
vel("Cyclohexane", "T", 300, "Dmolar", 9.1756e3, "V", 863.66e-6, 1e-4),
vel("Cyclohexane", "T", 300, "Dmolar", 9.9508e3, "V", 2850.18e-6, 1e-4),
vel("Cyclohexane", "T", 500, "Dmolar", 1e-10, "V", 11.189e-6, 1e-4),
vel("Cyclohexane", "T", 500, "Dmolar", 6.0213e3, "V", 94.842e-6, 1e-4),
vel("Cyclohexane", "T", 500, "Dmolar", 8.5915e3, "V", 380.04e-6, 1e-4),
vel("Cyclohexane", "T", 700, "Dmolar", 1e-10, "V", 15.093e-6, 1e-4),
vel("Cyclohexane", "T", 700, "Dmolar", 7.4765e3, "V", 176.749e-6, 1e-4),
};
class TransportValidationFixture
@@ -255,7 +281,7 @@ TEST_CASE_METHOD(TransportValidationFixture, "Compare viscosities against publis
CAPTURE(el.in2);
CAPTURE(el.v2);
CHECK_NOTHROW(set_pair(el.in1, el.v1, el.in2, el.v2));
get_value(CoolProp::iviscosity);
CHECK_NOTHROW(get_value(CoolProp::iviscosity));
CAPTURE(el.expected);
CAPTURE(actual);
CHECK(fabs(actual/el.expected-1) < el.tol);