22 double HighPressureGasTransport::thermalConductivity()
26 double Lprime_m = 0.0;
27 const double c1 = 1./16.04;
28 size_t nsp = m_thermo->nSpecies();
29 vector<double> molefracs(nsp);
30 m_thermo->getMoleFractions(&molefracs[0]);
31 vector<double> cp_0_R(nsp);
32 m_thermo->getCp_R_ref(&cp_0_R[0]);
34 vector<double> L_i(nsp);
35 vector<double> f_i(nsp);
36 vector<double> h_i(nsp);
37 vector<double> V_k(nsp);
39 m_thermo -> getPartialMolarVolumes(&V_k[0]);
42 for (
size_t i = 0; i < m_nsp; i++) {
43 double Tc_i = Tcrit_i(i);
44 double Vc_i = Vcrit_i(i);
45 double T_r = m_thermo->temperature()/Tc_i;
46 double V_r = V_k[i]/Vc_i;
47 double T_p = std::min(T_r,2.0);
48 double V_p = std::max(0.5,std::min(V_r,2.0));
51 double theta_p = 1.0 + (m_w_ac[i] - 0.011)*(0.56553
52 - 0.86276*log(T_p) - 0.69852/T_p);
53 double phi_p = (1.0 + (m_w_ac[i] - 0.011)*(0.38560
54 - 1.1617*log(T_p)))*0.288/Zcrit_i(i);
55 double f_fac = Tc_i*theta_p/190.4;
56 double h_fac = 1000*Vc_i*phi_p/99.2;
57 double T_0 = m_temp/f_fac;
58 double mu_0 = 1e-7*(2.90774e6/T_0 - 3.31287e6*pow(T_0,-2./3.)
59 + 1.60810e6*pow(T_0,-1./3.) - 4.33190e5 + 7.06248e4*pow(T_0,1./3.)
60 - 7.11662e3*pow(T_0,2./3.) + 4.32517e2*T_0 - 1.44591e1*pow(T_0,4./3.)
61 + 2.03712e-1*pow(T_0,5./3.));
62 double H = sqrt(f_fac*16.04/m_mw[i])*pow(h_fac,-2./3.);
63 double mu_i = mu_0*H*m_mw[i]*c1;
64 L_i[i] = mu_i*1.32*
GasConstant*(cp_0_R[i] - 2.5)/m_mw[i];
65 L_i_min = min(L_i_min,L_i[i]);
67 double theta_s = 1 + (m_w_ac[i] - 0.011)*(0.09057 - 0.86276*log(T_p)
68 + (0.31664 - 0.46568/T_p)*(V_p - 0.5));
69 double phi_s = (1 + (m_w_ac[i] - 0.011)*(0.39490*(V_p - 1.02355)
70 - 0.93281*(V_p - 0.75464)*log(T_p)))*0.288/Zcrit_i(i);
71 f_i[i] = Tc_i*theta_s/190.4;
72 h_i[i] = 1000*Vc_i*phi_s/99.2;
78 for (
size_t i = 0; i < m_nsp; i++) {
79 for (
size_t j = 0; j < m_nsp; j++) {
81 double L_ij = 2*L_i[i]*L_i[j]/(L_i[i] + L_i[j] +
Tiny);
82 Lprime_m += molefracs[i]*molefracs[j]*L_ij;
84 double f_ij = sqrt(f_i[i]*f_i[j]);
85 double h_ij = 0.125*pow(pow(h_i[i],1./3.) + pow(h_i[j],1./3.),3.);
86 double mw_ij_inv = (m_mw[i] + m_mw[j])/(2*m_mw[i]*m_mw[j]);
87 f_m += molefracs[i]*molefracs[j]*f_ij*h_ij;
88 h_m += molefracs[i]*molefracs[j]*h_ij;
89 mw_m += molefracs[i]*molefracs[j]*sqrt(mw_ij_inv*f_ij)*pow(h_ij,-4./3.);
94 mw_m = pow(mw_m,-2.)*f_m*pow(h_m,-8./3.);
96 double rho_0 = 16.04*h_m/(1000*m_thermo->molarVolume());
97 double T_0 = m_temp/f_m;
98 double mu_0 = 1e-7*(2.90774e6/T_0 - 3.31287e6*pow(T_0,-2./3.)
99 + 1.60810e6*pow(T_0,-1./3.) - 4.33190e5 + 7.06248e4
100 *pow(T_0,1./3.) - 7.11662e3*pow(T_0,2./3.) + 4.32517e2*T_0
101 - 1.44591e1*pow(T_0,4./3.) + 2.03712e-1*pow(T_0,5./3.));
102 double L_1m = 1944*mu_0;
103 double L_2m = (-2.5276e-4 + 3.3433e-4*pow(1.12 - log(T_0/1.680e2),2))*rho_0;
104 double L_3m = exp(-7.19771 + 85.67822/T_0)*(exp((12.47183
105 - 984.6252*pow(T_0,-1.5))*pow(rho_0,0.1) + (rho_0/0.1617 - 1)
106 *sqrt(rho_0)*(0.3594685 + 69.79841/T_0 - 872.8833*pow(T_0,-2))) - 1.)*1e-3;
107 double H_m = sqrt(f_m*16.04/mw_m)*pow(h_m,-2./3.);
108 double Lstar_m = H_m*(L_1m + L_2m + L_3m);
109 return Lprime_m + Lstar_m;
112 void HighPressureGasTransport::getThermalDiffCoeffs(
double*
const dt)
123 void HighPressureGasTransport::getBinaryDiffCoeffs(
const size_t ld,
double*
const d)
125 vector<double> PcP(5);
126 size_t nsp = m_thermo->nSpecies();
127 vector<double> molefracs(nsp);
128 m_thermo->getMoleFractions(&molefracs[0]);
138 throw CanteraError(
"HighPressureGasTransport::getBinaryDiffCoeffs",
141 double rp = 1.0/m_thermo->pressure();
142 for (
size_t i = 0; i < nsp; i++) {
143 for (
size_t j = 0; j < nsp; j++) {
146 double x_i = std::max(
Tiny, molefracs[i]);
147 double x_j = std::max(
Tiny, molefracs[j]);
150 x_i = x_i/(x_i + x_j);
151 x_j = x_j/(x_i + x_j);
154 double Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j));
155 double Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j));
164 P_corr_ij = setPcorr(Pr_ij, Tr_ij);
175 d[ld*j + i] = P_corr_ij*rp * m_bdiff(i,j);
180 void HighPressureGasTransport::getMultiDiffCoeffs(
const size_t ld,
double*
const d)
197 size_t nsp = m_thermo->nSpecies();
198 vector<double> molefracs(nsp);
199 m_thermo->getMoleFractions(&molefracs[0]);
206 throw CanteraError(
"HighPressureGasTransport::getMultiDiffCoeffs",
209 for (
size_t i = 0; i < m_nsp; i++) {
210 for (
size_t j = 0; j < m_nsp; j++) {
213 double x_i = std::max(
Tiny, molefracs[i]);
214 double x_j = std::max(
Tiny, molefracs[j]);
217 double Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j));
218 double Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j));
224 P_corr_ij = setPcorr(Pr_ij, Tr_ij);
230 m_bdiff(i,j) *= P_corr_ij;
233 m_bindiff_ok =
false;
241 eval_L0000(molefracs.data());
245 int ierr =
invert(m_Lmatrix, m_nsp);
247 throw CanteraError(
"HighPressureGasTransport::getMultiDiffCoeffs",
248 "invert returned ierr = {}", ierr);
251 m_lmatrix_soln_ok =
false;
253 double prefactor = 16.0 * m_temp
254 *m_thermo->meanMolecularWeight()/(25.0*m_thermo->pressure());
256 for (
size_t i = 0; i < m_nsp; i++) {
257 for (
size_t j = 0; j < m_nsp; j++) {
258 double c = prefactor/m_mw[j];
259 d[ld*j + i] = c*molefracs[i]*(m_Lmatrix(i,j) - m_Lmatrix(i,i));
264 double HighPressureGasTransport::viscosity()
268 double Pc_mix_n = 0.;
269 double Pc_mix_d = 0.;
270 double MW_mix = m_thermo->meanMolecularWeight();
271 double MW_H = m_mw[0];
272 double MW_L = m_mw[0];
275 double tKelvin = m_thermo->temperature();
276 double Pvp_mix = m_thermo->satPressure(tKelvin);
277 size_t nsp = m_thermo->nSpecies();
278 vector<double> molefracs(nsp);
279 m_thermo->getMoleFractions(&molefracs[0]);
281 double x_H = molefracs[0];
282 for (
size_t i = 0; i < m_nsp; i++) {
285 double Tc = Tcrit_i(i);
286 double Tr = tKelvin/Tc;
287 double Zc = Zcrit_i(i);
288 Tc_mix += Tc*molefracs[i];
289 Pc_mix_n += molefracs[i]*Zc;
290 Pc_mix_d += molefracs[i]*Vcrit_i(i);
293 if (m_mw[i] > MW_H) {
296 }
else if (m_mw[i] < MW_L) {
300 double mu_ri = 52.46*100000*m_dipole(i,i)*m_dipole(i,i)
303 FP_mix_o += molefracs[i];
304 }
else if (mu_ri < 0.075) {
305 FP_mix_o += molefracs[i]*(1. + 30.55*pow(0.292 - Zc, 1.72));
306 }
else { FP_mix_o += molefracs[i]*(1. + 30.55*pow(0.292 - Zc, 1.72)
307 *fabs(0.96 + 0.1*(Tr - 0.7)));
315 vector<string> spnames = m_thermo->speciesNames();
316 if (spnames[i] ==
"He") {
317 FQ_mix_o += molefracs[i]*FQ_i(1.38,Tr,m_mw[i]);
318 }
else if (spnames[i] ==
"H2") {
319 FQ_mix_o += molefracs[i]*(FQ_i(0.76,Tr,m_mw[i]));
320 }
else if (spnames[i] ==
"D2") {
321 FQ_mix_o += molefracs[i]*(FQ_i(0.52,Tr,m_mw[i]));
323 FQ_mix_o += molefracs[i];
327 double Tr_mix = tKelvin/Tc_mix;
328 double Pc_mix =
GasConstant*Tc_mix*Pc_mix_n/Pc_mix_d;
329 double Pr_mix = m_thermo->pressure()/Pc_mix;
330 double ratio = MW_H/MW_L;
331 double ksi = pow(
GasConstant*Tc_mix*3.6277*pow(10.0,53.0)/(pow(MW_mix,3)
332 *pow(Pc_mix,4)),1.0/6.0);
334 if (ratio > 9 && x_H > 0.05 && x_H < 0.7) {
335 FQ_mix_o *= 1 - 0.01*pow(ratio,0.87);
339 double Z1m = (0.807*pow(Tr_mix,0.618) - 0.357*exp(-0.449*Tr_mix)
340 + 0.340*exp(-4.058*Tr_mix)+0.018)*FP_mix_o*FQ_mix_o;
345 if (Pr_mix < Pvp_mix/Pc_mix) {
346 double alpha = 3.262 + 14.98*pow(Pr_mix,5.508);
347 double beta = 1.390 + 5.746*Pr_mix;
348 Z2m = 0.600 + 0.760*pow(Pr_mix,alpha) + (0.6990*pow(Pr_mix,beta) -
351 throw CanteraError(
"HighPressureGasTransport::viscosity",
352 "State is outside the limits of the Lucas model, Tr <= 1");
354 }
else if ((Tr_mix > 1.0) && (Tr_mix < 40.0)) {
355 if ((Pr_mix > 0.0) && (Pr_mix <= 100.0)) {
356 double a_fac = 0.001245*exp(5.1726*pow(Tr_mix,-0.3286))/Tr_mix;
357 double b_fac = a_fac*(1.6553*Tr_mix - 1.2723);
358 double c_fac = 0.4489*exp(3.0578*pow(Tr_mix,-37.7332))/Tr_mix;
359 double d_fac = 1.7368*exp(2.2310*pow(Tr_mix,-7.6351))/Tr_mix;
360 double f_fac = 0.9425*exp(-0.1853*pow(Tr_mix,0.4489));
362 Z2m = Z1m*(1 + a_fac*pow(Pr_mix,1.3088)/(b_fac*pow(Pr_mix,f_fac)
363 + pow(1+c_fac*pow(Pr_mix,d_fac),-1)));
365 throw CanteraError(
"HighPressureGasTransport::viscosity",
366 "State is outside the limits of the Lucas model, 1.0 < Tr < 40");
369 throw CanteraError(
"HighPressureGasTransport::viscosity",
370 "State is outside the limits of the Lucas model, Tr > 40");
377 return Z2m*(1 + (FP_mix_o - 1)*pow(Y,-3))*(1 + (FQ_mix_o - 1)
378 *(1/Y - 0.007*pow(log(Y),4)))/(ksi*FP_mix_o*FQ_mix_o);
382 double HighPressureGasTransport::Tcrit_i(
size_t i)
385 vector<double> molefracs = store(i, m_thermo->nSpecies());
387 double tc = m_thermo->critTemperature();
389 m_thermo->setMoleFractions(&molefracs[0]);
393 double HighPressureGasTransport::Pcrit_i(
size_t i)
396 vector<double> molefracs = store(i, m_thermo->nSpecies());
398 double pc = m_thermo->critPressure();
400 m_thermo->setMoleFractions(&molefracs[0]);
404 double HighPressureGasTransport::Vcrit_i(
size_t i)
407 vector<double> molefracs = store(i, m_thermo->nSpecies());
409 double vc = m_thermo->critVolume();
411 m_thermo->setMoleFractions(&molefracs[0]);
415 double HighPressureGasTransport::Zcrit_i(
size_t i)
418 vector<double> molefracs = store(i, m_thermo->nSpecies());
420 double zc = m_thermo->critCompressibility();
422 m_thermo->setMoleFractions(&molefracs[0]);
426 vector<double> HighPressureGasTransport::store(
size_t i,
size_t nsp)
428 vector<double> molefracs(nsp);
429 m_thermo->getMoleFractions(&molefracs[0]);
430 vector<double> mf_temp(nsp, 0.0);
432 m_thermo->setMoleFractions(&mf_temp[0]);
438 double HighPressureGasTransport::FQ_i(
double Q,
double Tr,
double MW)
440 return 1.22*pow(Q,0.15)*(1 + 0.00385*pow(pow(Tr - 12.,2.),1./MW)
441 *fabs(Tr-12)/(Tr-12));
446 double HighPressureGasTransport::setPcorr(
double Pr,
double Tr)
448 const static double Pr_lookup[17] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0,
449 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 4.0, 5.0};
450 const static double DP_Rt_lookup[17] = {1.01, 1.01, 1.01, 1.01, 1.01, 1.01,
451 1.01, 1.02, 1.02, 1.02, 1.02, 1.03, 1.03, 1.04, 1.05, 1.06, 1.07};
452 const static double A_ij_lookup[17] = {0.038042, 0.067433, 0.098317,
453 0.137610, 0.175081, 0.216376, 0.314051, 0.385736, 0.514553, 0.599184,
454 0.557725, 0.593007, 0.696001, 0.790770, 0.502100, 0.837452, 0.890390};
455 const static double B_ij_lookup[17] = {1.52267, 2.16794, 2.42910, 2.77605,
456 2.98256, 3.11384, 3.50264, 3.07773, 3.54744, 3.61216, 3.41882, 3.18415,
457 3.37660, 3.27984, 3.39031, 3.23513, 3.13001};
458 const static double C_ij_lookup[17] = {0., 0., 0., 0., 0., 0., 0., 0.141211,
459 0.278407, 0.372683, 0.504894, 0.678469, 0.665702, 0., 0.602907, 0., 0.};
460 const static double E_ij_lookup[17] = {1., 1., 1., 1., 1., 1., 1., 13.45454,
461 14., 10.00900, 8.57519, 10.37483, 11.21674, 1., 6.19043, 1., 1.};
468 frac = (Pr - Pr_lookup[0])/(Pr_lookup[1] - Pr_lookup[0]);
470 for (
int j = 1; j < 17; j++) {
471 if (Pr_lookup[j] > Pr) {
472 frac = (Pr - Pr_lookup[j-1])/(Pr_lookup[j] - Pr_lookup[j-1]);
484 double P_corr_1 = DP_Rt_lookup[Pr_i]*(1.0 - A_ij_lookup[Pr_i]
485 *pow(Tr,-B_ij_lookup[Pr_i]))*(1-C_ij_lookup[Pr_i]
486 *pow(Tr,-E_ij_lookup[Pr_i]));
487 double P_corr_2 = DP_Rt_lookup[Pr_i+1]*(1.0 - A_ij_lookup[Pr_i+1]
488 *pow(Tr,-B_ij_lookup[Pr_i+1]))*(1-C_ij_lookup[Pr_i+1]
489 *pow(Tr,-E_ij_lookup[Pr_i+1]));
490 return P_corr_1*(1.0-frac) + P_corr_2*frac;
Interface for class HighPressureGasTransport.
ThermoPhase object for the ideal gas equation of state - workhorse for Cantera (see Thermodynamic Pro...
Interface for class MultiTransport.
Header file defining class TransportFactory (see TransportFactory)
Base class for exceptions thrown by Cantera classes.
An error indicating that an unimplemented function has been called.
const double GasConstant
Universal Gas Constant [J/kmol/K].
Namespace for the Cantera kernel.
const double Tiny
Small number to compare differences of mole fractions against.
int invert(DenseMatrix &A, size_t nn)
invert A. A is overwritten with A^-1.
const double BigNumber
largest number to compare to inf.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...