github/CoolProp
1
1
Fork 0
mirror of https://github.com/CoolProp/CoolProp.git synced 2026-02-09 21:35:28 -05:00
Code Issues Packages Projects Releases Wiki Activity
Files
daf40b4117d9df623eb17c0cd528f8012acbe9c9
CoolProp/src/Backends/Helmholtz/TransportRoutines.cpp
Ian Bell 68d943c92a Updated reference for Heptane 
Signed-off-by: Ian Bell <ian.h.bell@gmail.com>
2014-06-10 15:24:00 +02:00

965 lines
37 KiB
C++
Raw Blame History

#include "TransportRoutines.h"
#include "CoolPropFluid.h"
namespace CoolProp{
long double TransportRoutines::viscosity_dilute_kinetic_theory(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.is_pure_or_pseudopure)
{
long double Tstar = HEOS.T()/HEOS.components[0]->transport.epsilon_over_k;
long double sigma_nm = HEOS.components[0]->transport.sigma_eta*1e9; // 1e9 to convert from m to nm
long double molar_mass_kgkmol = HEOS.molar_mass()*1000; // 1000 to convert from kg/mol to kg/kmol
// The nondimensional empirical collision integral from Neufeld
// Neufeld, P. D.; Janzen, A. R.; Aziz, R. A. Empirical Equations to Calculate 16 of the Transport Collision Integrals (l,s)*
// for the Lennard-Jones (12-6) Potential. J. Chem. Phys. 1972, 57, 1100-1102
long double OMEGA22 = 1.16145*pow(Tstar, static_cast<long double>(-0.14874))+0.52487*exp(-0.77320*Tstar)+2.16178*exp(-2.43787*Tstar);
// The dilute gas component -
return 26.692e-9*sqrt(molar_mass_kgkmol*HEOS.T())/(pow(sigma_nm, 2)*OMEGA22); // Pa-s
}
else{
throw NotImplementedError("TransportRoutines::viscosity_dilute_kinetic_theory is only for pure and pseudo-pure");
}
}
long double TransportRoutines::viscosity_dilute_collision_integral(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.is_pure_or_pseudopure)
{
// Retrieve values from the state class
CoolProp::ViscosityDiluteGasCollisionIntegralData &data = HEOS.components[0]->transport.viscosity_dilute.collision_integral;
const std::vector<long double> &a = data.a, &t = data.t;
const long double C = data.C, molar_mass = data.molar_mass;
long double S;
// Unit conversions and variable definitions
const long double Tstar = HEOS.T()/HEOS.components[0]->transport.epsilon_over_k;
const long double sigma_nm = HEOS.components[0]->transport.sigma_eta*1e9; // 1e9 to convert from m to nm
const long double molar_mass_kgkmol = molar_mass*1000; // 1000 to convert from kg/mol to kg/kmol
/// Both the collision integral \f$\mathfrak{S}^*\f$ and effective cross section \f$\Omega^{(2,2)}\f$ have the same form,
/// in general we don't care which is used. The are related through \f$\Omega^{(2,2)} = (5/4)\mathfrak{S}^*\f$
/// see Vesovic(JPCRD, 1990) for CO\f$_2\f$ for further information
long double summer = 0, lnTstar = log(Tstar);
for (std::size_t i = 0; i < a.size(); ++i)
{
summer += a[i]*pow(lnTstar,t[i]);
}
S = exp(summer);
// The dilute gas component
return C*sqrt(molar_mass_kgkmol*HEOS.T())/(pow(sigma_nm, 2)*S); // Pa-s
}
else{
throw NotImplementedError("TransportRoutines::viscosity_dilute_collision_integral is only for pure and pseudo-pure");
}
}
long double TransportRoutines::viscosity_dilute_powers_of_T(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.is_pure_or_pseudopure)
{
// Retrieve values from the state class
CoolProp::ViscosityDiluteGasPowersOfT &data = HEOS.components[0]->transport.viscosity_dilute.powers_of_T;
const std::vector<long double> &a = data.a, &t = data.t;
long double summer = 0, T = HEOS.T();
for (std::size_t i = 0; i < a.size(); ++i)
{
summer += a[i]*pow(T, t[i]);
}
return summer;
}
else{
throw NotImplementedError("TransportRoutines::viscosity_dilute_powers_of_T is only for pure and pseudo-pure");
}
}
long double TransportRoutines::viscosity_dilute_collision_integral_powers_of_T(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.is_pure_or_pseudopure)
{
// Retrieve values from the state class
CoolProp::ViscosityDiluteCollisionIntegralPowersOfTstarData &data = HEOS.components[0]->transport.viscosity_dilute.collision_integral_powers_of_Tstar;
const std::vector<long double> &a = data.a, &t = data.t;
long double summer = 0, Tstar = HEOS.T()/data.T_reducing;
for (std::size_t i = 0; i < a.size(); ++i)
{
summer += a[i]*pow(Tstar, t[i]);
}
return data.C*sqrt(HEOS.T())/summer;
}
else{
throw NotImplementedError("TransportRoutines::viscosity_dilute_collision_integral_powers_of_T is only for pure and pseudo-pure");
}
}
long double TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.is_pure_or_pseudopure)
{
CoolProp::ViscosityModifiedBatschinskiHildebrandData &HO = HEOS.components[0]->transport.viscosity_higher_order.modified_Batschinski_Hildebrand;
long double delta = HEOS.rhomolar()/HO.rhomolar_reduce, tau = HO.T_reduce/HEOS.T();
// The first term that is formed of powers of tau (Tc/T) and delta (rho/rhoc)
long double S = 0;
for (unsigned int i = 0; i < HO.a.size(); ++i){
S += HO.a[i]*pow(delta, HO.d1[i])*pow(tau, HO.t1[i])*exp(HO.gamma[i]*pow(delta, HO.l[i]));
}
// For the terms that multiplies the bracketed term with delta and delta0
long double F = 0;
for (unsigned int i = 0; i < HO.f.size(); ++i){
F += HO.f[i]*pow(delta, HO.d2[i])*pow(tau, HO.t2[i]);
}
// for delta_0
long double summer_numer = 0;
for (unsigned int i = 0; i < HO.g.size(); ++i){
summer_numer += HO.g[i]*pow(tau, HO.h[i]);
}
long double summer_denom = 0;
for (unsigned int i = 0; i < HO.p.size(); ++i){
summer_denom += HO.p[i]*pow(tau, HO.q[i]);
}
long double delta0 = summer_numer/summer_denom;
// The higher-order-term component
return S + F*(1/(delta0-delta)-1/delta0); // Pa-s
}
else{
throw NotImplementedError("TransportRoutines::viscosity_higher_order_modified_Batschinski_Hildebrand is only for pure and pseudo-pure");
}
}
long double TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.is_pure_or_pseudopure)
{
// Retrieve values from the state class
CoolProp::ViscosityRainWaterFriendData &data = HEOS.components[0]->transport.viscosity_initial.rainwater_friend;
const std::vector<long double> &b = data.b, &t = data.t;
long double B_eta, B_eta_star;
long double Tstar = HEOS.T()/HEOS.components[0]->transport.epsilon_over_k; // [no units]
long double sigma = HEOS.components[0]->transport.sigma_eta; // [m]
long double summer = 0;
for (unsigned int i = 0; i < b.size(); ++i){
summer += b[i]*pow(Tstar, t[i]);
}
B_eta_star = summer; // [no units]
B_eta = 6.02214129e23*pow(sigma, 3)*B_eta_star; // [m^3/mol]
return B_eta; // [m^3/mol]
}
else{
throw NotImplementedError("TransportRoutines::viscosity_initial_density_dependence_Rainwater_Friend 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));
double sum;
int i,j;
// Dilute-gas component
*mubar_0=100.0*sqrt(Tbar)/(1.67752+2.20462/Tbar+0.6366564/powInt(Tbar,2)-0.241605/powInt(Tbar,3));
//Fill in zeros in H
for (i=0;i<=5;i++)
{
for (j=0;j<=6;j++)
{
H[i][j]=0;
}
}
//Set non-zero parameters of H
H[0][0]=5.20094e-1;
H[1][0]=8.50895e-2;
H[2][0]=-1.08374;
H[3][0]=-2.89555e-1;
H[0][1]=2.22531e-1;
H[1][1]=9.99115e-1;
H[2][1]=1.88797;
H[3][1]=1.26613;
H[5][1]=1.20573e-1;
H[0][2]=-2.81378e-1;
H[1][2]=-9.06851e-1;
H[2][2]=-7.72479e-1;
H[3][2]=-4.89837e-1;
H[4][2]=-2.57040e-1;
H[0][3]=1.61913e-1;
H[1][3]=2.57399e-1;
H[0][4]=-3.25372e-2;
H[3][4]=6.98452e-2;
H[4][5]=8.72102e-3;
H[3][6]=-4.35673e-3;
H[5][6]=-5.93264e-4;
// Finite density component
sum=0;
for (i=0;i<=5;i++)
{
for (j=0;j<=6;j++)
{
sum+=powInt(1/Tbar-1,i)*(H[i][j]*powInt(rhobar-1,j));
}
}
*mubar_1=exp(rhobar*sum);
}
long double TransportRoutines::viscosity_water_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
double x_mu=0.068,qc=1/1.9,qd=1/1.1,nu=0.630,gamma=1.239,zeta_0=0.13,LAMBDA_0=0.06,Tbar_R=1.5, pstar, Tstar, rhostar;
double delta,tau,mubar_0,mubar_1,mubar_2,drhodp,drhodp_R,DeltaChibar,zeta,w,L,Y,psi_D,Tbar,rhobar;
double drhobar_dpbar,drhobar_dpbar_R,R_Water;
pstar = 22.064e6; // [Pa]
Tstar = 647.096; // [K]
rhostar = 322; // [kg/m^3]
Tbar = HEOS.T()/Tstar;
rhobar = HEOS.rhomass()/rhostar;
R_Water = HEOS.gas_constant()/HEOS.molar_mass(); // [J/kg/K]
// Dilute and finite gas portions
visc_Helper(Tbar, rhobar, &mubar_0, &mubar_1);
///************************ Critical Enhancement ************************
delta=rhobar;
// "Normal" calculation
tau=1/Tbar;
drhodp=1/(R_Water*HEOS.T()*(1+2*delta*HEOS.dalphar_dDelta()+delta*delta*HEOS.d2alphar_dDelta2()));
drhobar_dpbar = pstar/rhostar*drhodp;
// "Reducing" calculation
tau=1/Tbar_R;
drhodp_R=1/(R_Water*Tbar_R*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau, delta)));
drhobar_dpbar_R = pstar/rhostar*drhodp_R;
DeltaChibar=rhobar*(drhobar_dpbar-drhobar_dpbar_R*Tbar_R/Tbar);
if (DeltaChibar<0)
DeltaChibar=0;
zeta=zeta_0*pow(DeltaChibar/LAMBDA_0,nu/gamma);
if (zeta<0.3817016416){
Y=1.0/5.0*qc*zeta*powInt(qd*zeta,5)*(1-qc*zeta+powInt(qc*zeta,2)-765.0/504.0*powInt(qd*zeta,2));
}
else
{
psi_D=acos(pow(1+powInt(qd*zeta,2),-1.0/2.0));
w=sqrt(fabs((qc*zeta-1)/(qc*zeta+1)))*tan(psi_D/2.0);
if (qc*zeta>1){
L=log((1+w)/(1-w));
}
else{
L=2*atan(fabs(w));
}
Y=1.0/12.0*sin(3*psi_D)-1/(4*qc*zeta)*sin(2*psi_D)+1.0/powInt(qc*zeta,2)*(1-5.0/4.0*powInt(qc*zeta,2))*sin(psi_D)-1.0/powInt(qc*zeta,3)*((1-3.0/2.0*powInt(qc*zeta,2))*psi_D-pow(fabs(powInt(qc*zeta,2)-1),3.0/2.0)*L);
}
mubar_2=exp(x_mu*Y);
return (mubar_0*mubar_1*mubar_2)/1e6;
}
long double TransportRoutines::viscosity_hydrogen_higher_order_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
long double Tr = HEOS.T()/33.145, rhor = HEOS.keyed_output(CoolProp::iDmass)/90.5;
long double c[] = {0, 6.43449673e-6, 4.56334068e-2, 2.32797868e-1, 9.58326120e-1, 1.27941189e-1, 3.63576595e-1};
return c[1]*pow(rhor,2)*exp(c[2]*Tr+c[3]/Tr+c[4]*pow(rhor,2)/(c[5]+Tr)+c[6]*pow(rhor,6));
}
long double TransportRoutines::viscosity_hexane_higher_order_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
long double Tr = HEOS.T()/507.82, rhor = HEOS.keyed_output(CoolProp::iDmass)/233.182;
// Output is in Pa-s
double c[] = {2.53402335/1e6, -9.724061002/1e6, 0.469437316, 158.5571631, 72.42916856/1e6, 10.60751253, 8.628373915, -6.61346441, -2.212724566};
return pow(rhor,static_cast<long double>(2.0/3.0))*sqrt(Tr)*(c[0]/Tr+c[1]/(c[2]+Tr+c[3]*rhor*rhor)+c[4]*(1+rhor)/(c[5]+c[6]*Tr+c[7]*rhor+rhor*rhor+c[8]*rhor*Tr));
}
long double TransportRoutines::viscosity_heptane_higher_order_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
/// From Michailidou-JPCRD-2014-Heptane
long double Tr = HEOS.T()/540.13, rhor = HEOS.rhomass()/232;
// Output is in Pa-s
double c[] = {0, 22.15000/1e6, -15.00870/1e6, 3.71791/1e6, 77.72818/1e6, 9.73449, 9.51900, -6.34076, -2.51909};
return pow(rhor,2.0L/3.0L)*sqrt(Tr)*(c[1]*rhor+c[2]*pow(rhor,2)+c[3]*pow(rhor,3)+c[4]*rhor/(c[5]+c[6]*Tr+c[7]*rhor+rhor*rhor+c[8]*rhor*Tr));
}
long double TransportRoutines::viscosity_higher_order_friction_theory(HelmholtzEOSMixtureBackend &HEOS)
{
if (HEOS.is_pure_or_pseudopure)
{
CoolProp::ViscosityFrictionTheoryData &F = HEOS.components[0]->transport.viscosity_higher_order.friction_theory;
long double tau = F.T_reduce/HEOS.T(), kii, krrr, kaaa, krr, kdrdr;
double psi1 = exp(tau)-F.c1;
double psi2 = exp(pow(tau,2))-F.c2;
double ki = (F.Ai[0] + F.Ai[1]*psi1 + F.Ai[2]*psi2)*tau;
double ka = (F.Aa[0] + F.Aa[1]*psi1 + F.Aa[2]*psi2)*pow(tau, F.Na);
double kr = (F.Ar[0] + F.Ar[1]*psi1 + F.Ar[2]*psi2)*pow(tau, F.Nr);
double kaa = (F.Aaa[0] + F.Aaa[1]*psi1 + F.Aaa[2]*psi2)*pow(tau, F.Naa);
if (F.Arr.empty()){
krr = 0;
kdrdr = (F.Adrdr[0] + F.Adrdr[1]*psi1 + F.Adrdr[2]*psi2)*pow(tau, F.Nrr);
}
else{
krr = (F.Arr[0] + F.Arr[1]*psi1 + F.Arr[2]*psi2)*pow(tau, F.Nrr);
kdrdr = 0;
}
if (!F.Aii.empty() && !F.Arrr.empty() && !F.Aaaa.empty()){
kii = (F.Aii[0] + F.Aii[1]*psi1 + F.Aii[2]*psi2)*pow(tau, F.Nii);
krrr = (F.Arrr[0] + F.Arrr[1]*psi1 + F.Arrr[2]*psi2)*pow(tau, F.Nrrr);
kaaa = (F.Aaaa[0] + F.Aaaa[1]*psi1 + F.Aaaa[2]*psi2)*pow(tau, F.Naaa);
}
else{
kii = 0;
krrr = 0;
kaaa = 0;
}
double p = HEOS.p()/1e5; // [bar]; 1e5 for conversion from Pa -> bar
double pr = HEOS.T()*HEOS.first_partial_deriv(CoolProp::iP, CoolProp::iT, CoolProp::iDmolar)/1e5; // [bar/K]; 1e5 for conversion from Pa -> bar
double pa = p - pr; //[bar]
double pid = HEOS.rhomolar() * HEOS.gas_constant() * HEOS.T() / 1e5; // [bar]; 1e5 for conversion from Pa -> bar
double deltapr = pr - pid;
double eta_f = ka*pa + kr*deltapr + ki*pid + kaa*pa*pa + kdrdr*deltapr*deltapr + krr*pr*pr + kii*pid*pid + krrr*pr*pr*pr + kaaa*pa*pa*pa;
return eta_f; //[Pa-s]
}
else{
throw NotImplementedError("TransportRoutines::viscosity_higher_order_friction_theory is only for pure and pseudo-pure");
}
}
long double TransportRoutines::viscosity_helium_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
double eta_0,eta_0_slash, eta_E_slash, B,C,D,ln_eta,x;
//
// Arp, V.D., McCarty, R.D., and Friend, D.G.,
// "Thermophysical Properties of Helium-4 from 0.8 to 1500 K with Pressures to 2000 MPa",
// NIST Technical Note 1334 (revised), 1998.
//
// Using Arp NIST report
// Report is not clear on viscosity, referring to REFPROP source code for clarity
// Correlation wants density in g/cm^3; kg/m^3 --> g/cm^3, divide by 1000
long double rho = HEOS.keyed_output(CoolProp::iDmass)/1000.0, T = HEOS.T();
if (T <= 300){
x = log(T);
}
else{
x = log(300.0);
}
// Evaluate the terms B,C,D
B = -47.5295259/x+87.6799309-42.0741589*x+8.33128289*x*x-0.589252385*x*x*x;
C = 547.309267/x-904.870586+431.404928*x-81.4504854*x*x+5.37008433*x*x*x;
D = -1684.39324/x+3331.08630-1632.19172*x+308.804413*x*x-20.2936367*x*x*x;
eta_0_slash = -0.135311743/x+1.00347841+1.20654649*x-0.149564551*x*x+0.012520841*x*x*x;
eta_E_slash = rho*B+rho*rho*C+rho*rho*rho*D;
if (T<=100)
{
ln_eta = eta_0_slash + eta_E_slash;
// Correlation yields viscosity in micro g/(cm-s); to get Pa-s, divide by 10 to get micro Pa-s, then another 1e6 to get Pa-s
return exp(ln_eta)/10.0/1e6;
}
else
{
ln_eta = eta_0_slash + eta_E_slash;
eta_0 = 196*pow(T,static_cast<long double>(0.71938))*exp(12.451/T-295.67/T/T-4.1249);
// Correlation yields viscosity in micro g/(cm-s); to get Pa-s, divide by 10 to get micro Pa-s, then another 1e6 to get Pa-s
return (exp(ln_eta)+eta_0-exp(eta_0_slash))/10.0/1e6;
}
}
long double TransportRoutines::viscosity_R23_hardcoded(HelmholtzEOSMixtureBackend &HEOS)
{
double C1 = 1.3163, //
C2 = 0.1832,
DeltaGstar = 771.23,
rhoL = 32.174,
rhocbar = 7.5114,
Tc = 299.2793,
DELTAeta_max = 3.967,
Ru = 8.31451,
molar_mass = 70.014;
double a[] = {0.4425728, -0.5138403, 0.1547566, -0.02821844, 0.001578286};
double e_k = 243.91, sigma = 0.4278;
double Tstar = HEOS.T()/e_k;
double logTstar = log(Tstar);
double Omega = exp(a[0]+a[1]*logTstar+a[2]*pow(logTstar,2)+a[3]*pow(logTstar,3)+a[4]*pow(logTstar,4));
double eta_DG = 1.25*0.021357*sqrt(molar_mass*HEOS.T())/(sigma*sigma*Omega); // uPa-s
double rhobar = HEOS.rhomolar()/1000; // [mol/L]
double eta_L = C2*(rhoL*rhoL)/(rhoL-rhobar)*sqrt(HEOS.T())*exp(rhobar/(rhoL-rhobar)*DeltaGstar/(Ru*HEOS.T()));
double chi = rhobar - rhocbar;
double tau = HEOS.T() - Tc;
double DELTAeta_c = 4*DELTAeta_max/((exp(chi)+exp(-chi))*(exp(tau)+exp(-tau)));
return (pow((rhoL-rhobar)/rhoL,C1)*eta_DG+pow(rhobar/rhoL,C1)*eta_L+DELTAeta_c)/1e6;
}
long double TransportRoutines::viscosity_dilute_ethane(HelmholtzEOSMixtureBackend &HEOS)
{
double C[] = {0, -3.0328138281, 16.918880086, -37.189364917, 41.288861858, -24.615921140, 8.9488430959, -1.8739245042, 0.20966101390, -9.6570437074e-3};
double OMEGA_2_2 = 0, e_k = 245, Tstar;
Tstar = HEOS.T()/e_k;
for (int i = 1; i<= 9; i++)
{
OMEGA_2_2 += C[i]*pow(Tstar,(i-1)/3.0-1);
}
return 12.0085*sqrt(Tstar)*OMEGA_2_2/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};
double s[] = {0,0,1,0,1,1.5,0,2,0,1,0,1};
double g[] = {0, 0.47177003, -0.23950311, 0.39808301, -0.27343335, 0.35192260, -0.21101308, -0.00478579, 0.07378129, -0.030435255, -0.30435286, 0.001215675};
double sum1 = 0, sum2 = 0, tau = 305.33/HEOS.T(), delta = HEOS.rhomolar()/6870;
for (int i = 1; i<= 9; ++i){
sum1 += g[i]*pow(delta,r[i])*pow(tau,s[i]);
}
for (int i = 10; i<= 11; ++i){
sum2 += g[i]*pow(delta,r[i])*pow(tau,s[i]);
}
return 15.977*sum1/(1+sum2)/1e6;
}
long double TransportRoutines::conductivity_dilute_ratio_polynomials(HelmholtzEOSMixtureBackend &HEOS){
if (HEOS.is_pure_or_pseudopure)
{
// Retrieve values from the state class
CoolProp::ConductivityDiluteRatioPolynomialsData &data = HEOS.components[0]->transport.conductivity_dilute.ratio_polynomials;
long double summer1 = 0, summer2 = 0, Tr = HEOS.T()/data.T_reducing;
for (std::size_t i = 0; i < data.A.size(); ++i)
{
summer1 += data.A[i]*pow(Tr, data.n[i]);
}
for (std::size_t i = 0; i < data.B.size(); ++i)
{
summer2 += data.B[i]*pow(Tr, data.m[i]);
}
return summer1/summer2;
}
else{
throw NotImplementedError("TransportRoutines::conductivity_dilute_ratio_polynomials is only for pure and pseudo-pure");
}
};
long double TransportRoutines::conductivity_residual_polynomial(HelmholtzEOSMixtureBackend &HEOS){
if (HEOS.is_pure_or_pseudopure)
{
// Retrieve values from the state class
CoolProp::ConductivityResidualPolynomialData &data = HEOS.components[0]->transport.conductivity_residual.polynomials;
long double summer = 0, tau = data.T_reducing/HEOS.T(), delta = HEOS.keyed_output(CoolProp::iDmass)/data.rhomass_reducing;
for (std::size_t i = 0; i < data.B.size(); ++i)
{
summer += data.B[i]*pow(tau, data.t[i])*pow(delta, data.d[i]);
}
return summer;
}
else{
throw NotImplementedError("TransportRoutines::conductivity_residual_polynomial is only for pure and pseudo-pure");
}
};
long double TransportRoutines::conductivity_residual_polynomial_and_exponential(HelmholtzEOSMixtureBackend &HEOS){
if (HEOS.is_pure_or_pseudopure)
{
// Retrieve values from the state class
CoolProp::ConductivityResidualPolynomialAndExponentialData &data = HEOS.components[0]->transport.conductivity_residual.polynomial_and_exponential;
long double summer = 0, tau = HEOS.tau(), delta = HEOS.delta();
for (std::size_t i = 0; i < data.A.size(); ++i)
{
summer += data.A[i]*pow(tau, data.t[i])*pow(delta, data.d[i])*exp(-data.gamma[i]*pow(delta,data.l[i]));
}
return summer;
}
else{
throw NotImplementedError("TransportRoutines::conductivity_residual_polynomial_and_exponential is only for pure and pseudo-pure");
}
};
long double TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers(HelmholtzEOSMixtureBackend &HEOS){
if (HEOS.is_pure_or_pseudopure)
{
// Olchowy and Sengers cross-over term
// Retrieve values from the state class
CoolProp::ConductivityCriticalSimplifiedOlchowySengersData &data = HEOS.components[0]->transport.conductivity_critical.Olchowy_Sengers;
double k = data.k,
R0 = data.R0,
nu = data.nu,
gamma = data.gamma,
GAMMA = data.GAMMA,
zeta0 = data.zeta0,
qD = data.qD,
Tc = HEOS.get_reducing().T, // [K]
rhoc = HEOS.get_reducing().rhomolar, // [mol/m^3]
Pcrit = HEOS.get_reducing().p, // [Pa]
Tref, // [K]
cp,cv,delta,num,zeta,mu,pi=M_PI,tau,
OMEGA_tilde,OMEGA_tilde0;
if (ValidNumber(data.T_ref))
Tref = data.T_ref;
else
Tref = 1.5*Tc;
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;
double tau_ref = Tc/Tref;
double dp_drho_ref = HEOS.gas_constant()*Tref*(1+2*delta*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,tau_ref,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,tau_ref,delta));
double Xref = Pcrit/pow(rhoc, 2)*HEOS.rhomolar()/dp_drho_ref*Tref/HEOS.T();
num = X - Xref;
// No critical enhancement if numerator is negative, zero, or just a tiny bit positive due to roundoff
// See also Lemmon, IJT, 2004, page 27
if (num < DBL_EPSILON*10)
return 0.0;
else
zeta = zeta0*pow(num/GAMMA, nu/gamma); //[m]
cp = HEOS.cpmolar(); //[J/mol/K]
cv = HEOS.cvmolar(); //[J/mol/K]
mu = HEOS.viscosity(); //[Pa-s]
OMEGA_tilde = 2.0/pi*((cp-cv)/cp*atan(zeta*qD)+cv/cp*(zeta*qD)); //[-]
OMEGA_tilde0 = 2.0/pi*(1.0-exp(-1.0/(1.0/(qD*zeta)+1.0/3.0*(zeta*qD)*(zeta*qD)/delta/delta))); //[-]
double lambda = HEOS.rhomolar()*cp*R0*k*HEOS.T()/(6*pi*mu*zeta)*(OMEGA_tilde - OMEGA_tilde0); //[W/m/K]
return lambda; //[W/m/K]
}
else{
throw NotImplementedError("TransportRoutines::conductivity_critical_simplified_Olchowy_Sengers is only for pure and pseudo-pure");
}
};
long double TransportRoutines::conductivity_critical_hardcoded_R123(HelmholtzEOSMixtureBackend &HEOS){
double a13 = 0.486742e-2, a14 = -100, a15 = -7.08535;
return a13*exp(a14*pow(HEOS.tau()-1,4)+a15*pow(HEOS.delta()-1,2));
};
long double TransportRoutines::conductivity_critical_hardcoded_CO2_ScalabrinJPCRD2006(HelmholtzEOSMixtureBackend &HEOS){
long double nc = 0.775547504e-3*4.81384, Tr = HEOS.T()/304.1282, alpha, rhor = HEOS.keyed_output(iDmass)/467.6;
static long double a[] = {0.0, 3.0, 6.70697, 0.94604, 0.30, 0.30, 0.39751, 0.33791, 0.77963, 0.79857, 0.90, 0.02, 0.20};
// Equation 6 from Scalabrin
alpha = 1-a[10]*acosh(1+a[11]*pow(pow(1-Tr,2),a[12]));
// Equation 5 from Scalabrin
long double numer = rhor*exp(-pow(rhor,a[1])/a[1]-pow(a[2]*(Tr-1),2)-pow(a[3]*(rhor-1),2));
long double braced = (1-1/Tr)+a[4]*pow(pow(rhor-1,2),0.5/a[5]);
long double denom = pow(pow(pow(braced, 2), a[6]) + pow(pow(a[7]*(rhor-alpha), 2), a[8]),a[9]);
return nc*numer/denom;
}
long double TransportRoutines::conductivity_dilute_hardcoded_CO2(HelmholtzEOSMixtureBackend &HEOS){
double e_k = 251.196, Tstar;
double b[] = {0.4226159, 0.6280115, -0.5387661, 0.6735941, 0, 0, -0.4362677, 0.2255388};
double c[] = {0, 2.387869e-2, 4.350794, -10.33404, 7.981590, -1.940558};
//Vesovic Eq. 31 [no units]
double summer = 0;
for (int i=1; i<=5; i++)
summer += c[i]*pow(HEOS.T()/100.0, 2-i);
double cint_k = 1.0 + exp(-183.5/HEOS.T())*summer;
//Vesovic Eq. 12 [no units]
double r = sqrt(2.0/5.0*cint_k);
// According to REFPROP, 1+r^2 = cp-2.5R. This is unclear to me but seems to suggest that cint/k is the difference
// between the ideal gas specific heat and a monatomic specific heat of 5/2*R. Using the form of cint/k from Vesovic
// does not yield exactly the correct values
Tstar = HEOS.T()/e_k;
//Vesovic Eq. 30 [no units]
summer = 0;
for (int i=0; i<=7; i++)
summer += b[i]/pow(Tstar, i);
double Gstar_lambda = summer;
//Vesovic Eq. 29 [W/m/K]
double lambda_0 = 0.475598e-3*sqrt(HEOS.T())*(1+r*r)/Gstar_lambda;
return lambda_0;
}
long double TransportRoutines::conductivity_dilute_hardcoded_ethane(HelmholtzEOSMixtureBackend &HEOS){
double e_k = 245.0;
double tau = 305.33/HEOS.T(), Tstar = HEOS.T()/e_k;
double fint = 1.7104147-0.6936482/Tstar;
double lambda_0 = 0.276505e-3*(HEOS.calc_viscosity_dilute()*1e6)*(3.75-fint*(tau*tau*HEOS.d2alpha0_dTau2()+1.5)); //[W/m/K]
return lambda_0;
}
long double TransportRoutines::conductivity_dilute_eta0_and_poly(HelmholtzEOSMixtureBackend &HEOS){
if (HEOS.is_pure_or_pseudopure)
{
CoolProp::ConductivityDiluteEta0AndPolyData &E = HEOS.components[0]->transport.conductivity_dilute.eta0_and_poly;
double eta0_uPas = HEOS.calc_viscosity_dilute()*1e6; // [uPa-s]
double summer = E.A[0]*eta0_uPas;
for (std::size_t i=1; i < E.A.size(); ++i)
summer += E.A[i]*pow(static_cast<long double>(HEOS.tau()), E.t[i]);
return summer;
}
else{
throw NotImplementedError("TransportRoutines::conductivity_dilute_eta0_and_poly is only for pure and pseudo-pure");
}
}
long double TransportRoutines::conductivity_hardcoded_water(HelmholtzEOSMixtureBackend &HEOS){
double L[5][6] = {{1.60397357,-0.646013523,0.111443906,0.102997357,-0.0504123634,0.00609859258},
{2.33771842,-2.78843778,1.53616167,-0.463045512,0.0832827019,-0.00719201245},
{2.19650529,-4.54580785,3.55777244,-1.40944978,0.275418278,-0.0205938816},
{-1.21051378,1.60812989,-0.621178141,0.0716373224,0,0},
{-2.7203370,4.57586331,-3.18369245,1.1168348,-0.19268305,0.012913842}};
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;
int i,j;
Tbar = HEOS.T()/Tstar;
rhobar = HEOS.keyed_output(CoolProp::iDmass)/rhostar;
// Dilute gas contribution
lambdabar_0 = sqrt(Tbar)/(2.443221e-3+1.323095e-2/Tbar+6.770357e-3/pow(Tbar,2)-3.454586e-3/pow(Tbar,3)+4.096266e-4/pow(Tbar,4));
sum=0;
for (i=0;i<=4;i++){
for (j=0;j<=5;j++){
sum+=L[i][j]*powInt(1.0/Tbar-1.0,i)*powInt(rhobar-1,j);
}
}
// Finite density contribution
lambdabar_1=exp(rhobar*sum);
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;
double drhodp_Trbar = 1/(R*Tr_bar*Tstar*(1+2*rhobar*HEOS.calc_alphar_deriv_nocache(0,1,HEOS.mole_fractions,1/Tr_bar,delta)+delta*delta*HEOS.calc_alphar_deriv_nocache(0,2,HEOS.mole_fractions,1/Tr_bar,delta)));
double drhobar_dpbar_Trbar = pstar/rhostar*drhodp_Trbar;
double cp = HEOS.cpmass(); // [J/kg/K]
double cv = HEOS.cvmass(); // [J/kg/K]
double cpbar = cp/R; //[-]
double mubar = HEOS.viscosity()/mustar;
double DELTAchibar_T = rhobar*(drhobar_dpbar-drhobar_dpbar_Trbar*Tr_bar/Tbar);
if (DELTAchibar_T<0)
xi = 0;
else
xi = xi_0*pow(DELTAchibar_T/Lambda_0,nu/gamma);
double y = qd_bar*xi;
double Z;
double kappa = cp/cv;
if (y < 1.2e-7)
Z = 0;
else
Z = 2/(pi*y)*(((1-1/kappa)*atan(y)+y/kappa)-(1-exp(-1/(1/y+y*y/3/rhobar/rhobar))));
lambdabar_2 = GAMMA*rhobar*cpbar*Tbar/mubar*Z;
return (lambdabar_0*lambdabar_1+lambdabar_2)*lambdastar;
}
long double TransportRoutines::conductivity_hardcoded_R23(HelmholtzEOSMixtureBackend &HEOS){
double B1 = -2.5370, // [mW/m/K]
B2 = 0.05366, // [mW/m/K^2]
C1 = 0.94215, // [-]
C2 = 0.14914, // [mW/m/K^2]
DeltaGstar = 2508.58, //[J/mol]
rhoL = 68.345, // [mol/dm^3] = [mol/L]
rhocbar = 7.5114, // [mol/dm^3]
DELTAlambda_max = 25, //[mW/m/K]
Ru = 8.31451, // [J/mol/K]
Tc = 299.2793, //[K]
T = HEOS.T(); //[K]
double lambda_DG = B1 + B2*T;
double rhobar = HEOS.rhomolar()/1000; // [mol/L]
double lambda_L = C2*(rhoL*rhoL)/(rhoL-rhobar)*sqrt(T)*exp(rhobar/(rhoL-rhobar)*DeltaGstar/(Ru*T));
double chi = rhobar - rhocbar;
double tau = T - Tc;
double DELTAlambda_c = 4*DELTAlambda_max/((exp(chi)+exp(-chi))*(exp(tau)+exp(-tau)));
return (pow((rhoL-rhobar)/rhoL,C1)*lambda_DG+pow(rhobar/rhoL,C1)*lambda_L+DELTAlambda_c)/1e3;
}
long double TransportRoutines::conductivity_critical_hardcoded_ammonia(HelmholtzEOSMixtureBackend &HEOS){
/*
From "Thermal Conductivity of Ammonia in a Large
Temperature and Pressure Range Including the Critical Region"
by R. Tufeu, D.Y. Ivanov, Y. Garrabos, B. Le Neindre,
Bereicht der Bunsengesellschaft Phys. Chem. 88 (1984) 422-427
*/
double T = HEOS.T(), Tc = 405.4, rhoc = 235, rho;
double LAMBDA=1.2, nu=0.63, gamma =1.24, DELTA=0.50,t,zeta_0_plus=1.34e-10,a_zeta=1,GAMMA_0_plus=0.423e-8;
double pi=3.141592654,a_chi,k_B=1.3806504e-23,X_T,DELTA_lambda,dPdT,eta_B,DELTA_lambda_id,DELTA_lambda_i;
rho = HEOS.keyed_output(CoolProp::iDmass);
t = fabs((T-Tc)/Tc);
a_chi = a_zeta/0.7;
eta_B = (2.60+1.6*t)*1e-5;
dPdT = (2.18-0.12/exp(17.8*t))*1e5; // [Pa-K]
X_T = 0.61*rhoc+16.5*log(t);
// Along the critical isochore (only a function of temperature) (Eq. 9)
DELTA_lambda_i = LAMBDA*(k_B*T*T)/(6*pi*eta_B*(zeta_0_plus*pow(t,-nu)*(1+a_zeta*pow(t,DELTA))))*dPdT*dPdT*GAMMA_0_plus*pow(t,-gamma)*(1+a_chi*pow(t,DELTA));
DELTA_lambda_id = DELTA_lambda_i*exp(-36*t*t);
if (rho < 0.6*rhoc)
{
DELTA_lambda = DELTA_lambda_id*(X_T*X_T)/(X_T*X_T+powInt(0.6*rhoc-0.96*rhoc,2))*powInt(rho,2)/powInt(0.6*rhoc,2);
}
else
{
DELTA_lambda = DELTA_lambda_id*(X_T*X_T)/(X_T*X_T+powInt(rho-0.96*rhoc,2));
}
return DELTA_lambda;
}
long double TransportRoutines::conductivity_hardcoded_helium(HelmholtzEOSMixtureBackend &HEOS){
/*
What an incredibly annoying formulation! Implied coefficients?? Not cool.
*/
double rhoc = 68.0, lambda_e, lambda_c, T = HEOS.T(), rho = HEOS.rhomass();
double summer = 3.739232544/T-2.620316969e1/T/T+5.982252246e1/T/T/T-4.926397634e1/T/T/T/T;
double lambda_0 = 2.7870034e-3*pow(T, 7.034007057e-1)*exp(summer);
double c[]={ 1.862970530e-4,
-7.275964435e-7,
-1.427549651e-4,
3.290833592e-5,
-5.213335363e-8,
4.492659933e-8,
-5.924416513e-9,
7.087321137e-6,
-6.013335678e-6,
8.067145814e-7,
3.995125013e-7};
// Equation 17
lambda_e = (c[0]+c[1]*T+c[2]*pow(T,1/3.0)+c[3]*pow(T,2.0/3.0))*rho
+(c[4]+c[5]*pow(T,1.0/3.0)+c[6]*pow(T,2.0/3.0))*rho*rho*rho
+(c[7]+c[8]*pow(T,1.0/3.0)+c[9]*pow(T,2.0/3.0)+c[10]/T)*rho*rho*log(rho/rhoc);
// Critical component
lambda_c = 0.0;
if (3.5 < T && T < 12)
{
double x0 = 0.392, E1 = 2.8461, E2 = 0.27156, beta = 0.3554, gamma = 1.1743, delta = 4.304, rhoc_crit = 69.158,
Tc = 5.18992, pc = 2.2746e5;
double DeltaT = fabs(1-T/Tc), DeltaRho = fabs(1-rho/rhoc_crit);
double eta = HEOS.viscosity(); // [Pa-s]
double K_T = HEOS.isothermal_compressibility(), K_Tprime, K_Tbar;
double dpdT = HEOS.first_partial_deriv(CoolProp::iP, CoolProp::iT, CoolProp::iDmolar);
double W = pow(DeltaT/0.2,2) + pow(DeltaRho/0.25,2);
if (W > 1)
{
K_Tbar = K_T;
}
else
{
double x = pow(DeltaT/DeltaRho,1/beta);
double h = E1*(1 + x/x0)*pow(1 + E2*pow(1 + x/x0, 2/beta), (gamma-1)/(2*beta));
/**
dh/dx derived using sympy:
E1,x,x0,E2,beta,gamma = symbols('E1,x,x0,E2,beta,gamma')
h = E1*(1 + x/x0)*pow(1 + E2*pow(1 + x/x0, 2/beta), (gamma-1)/(2*beta))
ccode(simplify(diff(h,x)))
*/
double dhdx = E1*(E2*pow((x + x0)/x0, 2/beta)*(gamma - 1)*pow(E2*pow((x + x0)/x0, 2/beta) + 1, (1.0/2.0)*(gamma - 1)/beta) + pow(beta, 2)*pow(E2*pow((x + x0)/x0, 2/beta) + 1, (1.0/2.0)*(2*beta + gamma - 1)/beta))/(pow(beta, 2)*x0*(E2*pow((x + x0)/x0, 2/beta) + 1));
// Right-hand-side of Equation 9
double RHS = pow(DeltaRho,delta-1)*(delta*h-x/beta*dhdx);
K_Tprime = 1/(RHS*pow(rho/rhoc_crit,2)*pc);
K_Tbar = W*K_T + (1-W)*K_Tprime;
}
// 3.4685233d-17 and 3.726229668d0 are "magical" coefficients that are present in the REFPROP source to yield the right values. Not clear why these values are needed.
// Also, the form of the critical term in REFPROP does not agree with Hands paper. EL and MH from NIST are not sure where these coefficients come from.
lambda_c = 3.4685233e-17*3.726229668*sqrt(K_Tbar)*pow(T,2)/rho/eta*pow(dpdT,2)*exp(-18.66*pow(DeltaT,2)-4.25*pow(DeltaRho,4));
}
return lambda_0+lambda_e+lambda_c;
}
long double TransportRoutines::viscosity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference)
{
// Collect some parameters
long double M = HEOS.molar_mass(),
M0 = HEOS_Reference.molar_mass(),
Tc = HEOS.T_critical(),
Tc0 = HEOS_Reference.T_critical(),
rhocmolar = HEOS.rhomolar_critical(),
rhocmolar0 = HEOS_Reference.rhomolar_critical();
// Get a reference to the ECS data
CoolProp::ViscosityECSVariables &ECS = HEOS.components[0]->transport.viscosity_ecs;
// The correction polynomial psi_eta
double psi = 0;
for (std::size_t i=0; i < ECS.psi_a.size(); i++){
psi += ECS.psi_a[i]*pow(HEOS.rhomolar()/ECS.psi_rhomolar_reducing, ECS.psi_t[i]);
}
// The dilute gas portion for the fluid of interest [Pa-s]
long double eta_dilute = viscosity_dilute_kinetic_theory(HEOS);
/// \todo To be solved for...
// TODO: To be solved for...
long double theta = 1;
long double phi = 1;
// The equivalent substance reducing ratios
long double f = Tc/Tc0*theta;
long double h = rhocmolar0/rhocmolar*phi; // Must be the ratio of MOLAR densities!!
// To be solved for
long double T0 = HEOS.T()/f;
long double rhomolar0 = HEOS.rhomolar()*h;
// Update the reference fluid with the conformal state
HEOS_Reference.update(DmolarT_INPUTS, rhomolar0*psi, T0);
// The reference fluid's contribution to the viscosity [Pa-s]
long double eta_resid = HEOS_Reference.calc_viscosity_background();
// The F factor
long double F_eta = sqrt(f)*pow(h, static_cast<long double>(-2.0/3.0))*sqrt(M/M0);
// The total viscosity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
long double eta = eta_dilute + eta_resid*F_eta;
return eta;
}
long double TransportRoutines::conductivity_ECS(HelmholtzEOSMixtureBackend &HEOS, HelmholtzEOSMixtureBackend &HEOS_Reference)
{
// Collect some parameters
long double M = HEOS.molar_mass(),
M_kmol = M*1000,
M0 = HEOS_Reference.molar_mass(),
Tc = HEOS.T_critical(),
Tc0 = HEOS_Reference.T_critical(),
rhocmolar = HEOS.rhomolar_critical(),
rhocmolar0 = HEOS_Reference.rhomolar_critical(),
R_u = HEOS.gas_constant(),
R = HEOS.gas_constant()/HEOS.molar_mass(), //[J/kg/K]
R_kJkgK = R_u/M_kmol;
// Get a reference to the ECS data
CoolProp::ConductivityECSVariables &ECS = HEOS.components[0]->transport.conductivity_ecs;
// The correction polynomial psi_eta in rho/rho_red
double psi = 0;
for (std::size_t i=0; i < ECS.psi_a.size(); ++i){
psi += ECS.psi_a[i]*pow(HEOS.rhomolar()/ECS.psi_rhomolar_reducing, ECS.psi_t[i]);
}
// The correction polynomial f_int in T/T_red
double fint = 0;
for (std::size_t i=0; i < ECS.f_int_a.size(); ++i){
fint += ECS.f_int_a[i]*pow(HEOS.T()/ECS.f_int_T_reducing, ECS.f_int_t[i]);
}
// The dilute gas density for the fluid of interest [uPa-s]
long double eta_dilute = viscosity_dilute_kinetic_theory(HEOS)*1e6;
// The mass specific ideal gas constant-pressure specific heat [J/kg/K]
long double cp0 = HEOS.calc_cpmolar_idealgas()/HEOS.molar_mass();
// The internal contribution to the thermal conductivity [W/m/K]
long double lambda_int = fint*eta_dilute*(cp0 - 2.5*R)/1e3;
// The dilute gas contribution to the thermal conductivity [W/m/K]
long double lambda_dilute = 15.0e-3/4.0*R_kJkgK*eta_dilute;
/// \todo To be solved for...
// TODO: To be solved for...
long double theta = 1;
long double phi = 1;
// The equivalent substance reducing ratios
long double f = Tc/Tc0*theta;
long double h = rhocmolar0/rhocmolar*phi; // Must be the ratio of MOLAR densities!!
// To be solved for
long double T0 = HEOS.T()/f;
long double rhomolar0 = HEOS.rhomolar()*h;
// Update the reference fluid with the conformal state
HEOS_Reference.update(DmolarT_INPUTS, rhomolar0*psi, T0);
// The reference fluid's contribution to the conductivity [W/m/K]
long double lambda_resid = HEOS_Reference.calc_conductivity_background();
// The F factor
long double F_lambda = sqrt(f)*pow(h, static_cast<long double>(-2.0/3.0))*sqrt(M0/M);
// The critical contribution from the fluid of interest [W/m/K]
long double lambda_critical = conductivity_critical_simplified_Olchowy_Sengers(HEOS);
// The total thermal conductivity considering the contributions of the fluid of interest and the reference fluid [Pa-s]
long double lambda = lambda_int + lambda_dilute + lambda_resid*F_lambda + lambda_critical;
return lambda;
}
}; /* namespace CoolProp */
Reference in New Issue View Git Blame Copy Permalink