31 double Lprime_m = 0.0;
32 const double c1 = 1./16.04;
34 vector<double> molefracs(nsp);
36 vector<double> cp_0_R(nsp);
39 vector<double> L_i(nsp);
40 vector<double> f_i(nsp);
41 vector<double> h_i(nsp);
42 vector<double> V_k(nsp);
44 m_thermo -> getPartialMolarVolumes(&V_k[0]);
47 for (
size_t i = 0; i <
m_nsp; i++) {
48 double Tc_i = Tcrit_i(i);
49 double Vc_i = Vcrit_i(i);
51 double V_r = V_k[i]/Vc_i;
52 double T_p = std::min(T_r,2.0);
53 double V_p = std::max(0.5,std::min(V_r,2.0));
56 double theta_p = 1.0 + (
m_w_ac[i] - 0.011)*(0.56553
57 - 0.86276*log(T_p) - 0.69852/T_p);
58 double phi_p = (1.0 + (
m_w_ac[i] - 0.011)*(0.38560
59 - 1.1617*log(T_p)))*0.288/Zcrit_i(i);
60 double f_fac = Tc_i*theta_p/190.4;
61 double h_fac = 1000*Vc_i*phi_p/99.2;
63 double mu_0 = 1e-7*(2.90774e6/T_0 - 3.31287e6*pow(T_0,-2./3.)
64 + 1.60810e6*pow(T_0,-1./3.) - 4.33190e5 + 7.06248e4*pow(T_0,1./3.)
65 - 7.11662e3*pow(T_0,2./3.) + 4.32517e2*T_0 - 1.44591e1*pow(T_0,4./3.)
66 + 2.03712e-1*pow(T_0,5./3.));
67 double H = sqrt(f_fac*16.04/
m_mw[i])*pow(h_fac,-2./3.);
68 double mu_i = mu_0*H*
m_mw[i]*c1;
70 L_i_min = min(L_i_min,L_i[i]);
72 double theta_s = 1 + (
m_w_ac[i] - 0.011)*(0.09057 - 0.86276*log(T_p)
73 + (0.31664 - 0.46568/T_p)*(V_p - 0.5));
74 double phi_s = (1 + (
m_w_ac[i] - 0.011)*(0.39490*(V_p - 1.02355)
75 - 0.93281*(V_p - 0.75464)*log(T_p)))*0.288/Zcrit_i(i);
76 f_i[i] = Tc_i*theta_s/190.4;
77 h_i[i] = 1000*Vc_i*phi_s/99.2;
83 for (
size_t i = 0; i <
m_nsp; i++) {
84 for (
size_t j = 0; j <
m_nsp; j++) {
86 double L_ij = 2*L_i[i]*L_i[j]/(L_i[i] + L_i[j] +
Tiny);
87 Lprime_m += molefracs[i]*molefracs[j]*L_ij;
89 double f_ij = sqrt(f_i[i]*f_i[j]);
90 double h_ij = 0.125*pow(pow(h_i[i],1./3.) + pow(h_i[j],1./3.),3.);
92 f_m += molefracs[i]*molefracs[j]*f_ij*h_ij;
93 h_m += molefracs[i]*molefracs[j]*h_ij;
94 mw_m += molefracs[i]*molefracs[j]*sqrt(mw_ij_inv*f_ij)*pow(h_ij,-4./3.);
99 mw_m = pow(mw_m,-2.)*f_m*pow(h_m,-8./3.);
103 double mu_0 = 1e-7*(2.90774e6/T_0 - 3.31287e6*pow(T_0,-2./3.)
104 + 1.60810e6*pow(T_0,-1./3.) - 4.33190e5 + 7.06248e4
105 *pow(T_0,1./3.) - 7.11662e3*pow(T_0,2./3.) + 4.32517e2*T_0
106 - 1.44591e1*pow(T_0,4./3.) + 2.03712e-1*pow(T_0,5./3.));
107 double L_1m = 1944*mu_0;
108 double L_2m = (-2.5276e-4 + 3.3433e-4*pow(1.12 - log(T_0/1.680e2),2))*rho_0;
109 double L_3m = exp(-7.19771 + 85.67822/T_0)*(exp((12.47183
110 - 984.6252*pow(T_0,-1.5))*pow(rho_0,0.1) + (rho_0/0.1617 - 1)
111 *sqrt(rho_0)*(0.3594685 + 69.79841/T_0 - 872.8833*pow(T_0,-2))) - 1.)*1e-3;
112 double H_m = sqrt(f_m*16.04/mw_m)*pow(h_m,-2./3.);
113 double Lstar_m = H_m*(L_1m + L_2m + L_3m);
114 return Lprime_m + Lstar_m;
130 vector<double> PcP(5);
132 vector<double> molefracs(nsp);
143 throw CanteraError(
"HighPressureGasTransport::getBinaryDiffCoeffs",
147 for (
size_t i = 0; i < nsp; i++) {
148 for (
size_t j = 0; j < nsp; j++) {
151 double x_i = std::max(
Tiny, molefracs[i]);
152 double x_j = std::max(
Tiny, molefracs[j]);
155 x_i = x_i/(x_i + x_j);
156 x_j = x_j/(x_i + x_j);
159 double Tr_ij =
m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j));
169 P_corr_ij = setPcorr(Pr_ij, Tr_ij);
180 d[ld*j + i] = P_corr_ij*rp *
m_bdiff(i,j);
203 vector<double> molefracs(nsp);
211 throw CanteraError(
"HighPressureGasTransport::getMultiDiffCoeffs",
214 for (
size_t i = 0; i <
m_nsp; i++) {
215 for (
size_t j = 0; j <
m_nsp; j++) {
218 double x_i = std::max(
Tiny, molefracs[i]);
219 double x_j = std::max(
Tiny, molefracs[j]);
222 double Tr_ij =
m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j));
229 P_corr_ij = setPcorr(Pr_ij, Tr_ij);
252 throw CanteraError(
"HighPressureGasTransport::getMultiDiffCoeffs",
253 "invert returned ierr = {}", ierr);
256 m_lmatrix_soln_ok =
false;
258 double prefactor = 16.0 *
m_temp
261 for (
size_t i = 0; i <
m_nsp; i++) {
262 for (
size_t j = 0; j <
m_nsp; j++) {
263 double c = prefactor/
m_mw[j];
264 d[ld*j + i] = c*molefracs[i]*(m_Lmatrix(i,j) - m_Lmatrix(i,i));
273 double Pc_mix_n = 0.;
274 double Pc_mix_d = 0.;
276 double MW_H =
m_mw[0];
277 double MW_L =
m_mw[0];
283 vector<double> molefracs(nsp);
286 double x_H = molefracs[0];
287 for (
size_t i = 0; i <
m_nsp; i++) {
290 double Tc = Tcrit_i(i);
291 double Tr = tKelvin/Tc;
292 double Zc = Zcrit_i(i);
293 Tc_mix += Tc*molefracs[i];
294 Pc_mix_n += molefracs[i]*Zc;
295 Pc_mix_d += molefracs[i]*Vcrit_i(i);
298 if (
m_mw[i] > MW_H) {
301 }
else if (
m_mw[i] < MW_L) {
308 FP_mix_o += molefracs[i];
309 }
else if (mu_ri < 0.075) {
310 FP_mix_o += molefracs[i]*(1. + 30.55*pow(0.292 - Zc, 1.72));
311 }
else { FP_mix_o += molefracs[i]*(1. + 30.55*pow(0.292 - Zc, 1.72)
312 *fabs(0.96 + 0.1*(Tr - 0.7)));
321 if (spnames[i] ==
"He") {
322 FQ_mix_o += molefracs[i]*FQ_i(1.38,Tr,
m_mw[i]);
323 }
else if (spnames[i] ==
"H2") {
324 FQ_mix_o += molefracs[i]*(FQ_i(0.76,Tr,
m_mw[i]));
325 }
else if (spnames[i] ==
"D2") {
326 FQ_mix_o += molefracs[i]*(FQ_i(0.52,Tr,
m_mw[i]));
328 FQ_mix_o += molefracs[i];
332 double Tr_mix = tKelvin/Tc_mix;
333 double Pc_mix =
GasConstant*Tc_mix*Pc_mix_n/Pc_mix_d;
335 double ratio = MW_H/MW_L;
336 double ksi = pow(
GasConstant*Tc_mix*3.6277*pow(10.0,53.0)/(pow(MW_mix,3)
337 *pow(Pc_mix,4)),1.0/6.0);
339 if (ratio > 9 && x_H > 0.05 && x_H < 0.7) {
340 FQ_mix_o *= 1 - 0.01*pow(ratio,0.87);
344 double Z1m = (0.807*pow(Tr_mix,0.618) - 0.357*exp(-0.449*Tr_mix)
345 + 0.340*exp(-4.058*Tr_mix)+0.018)*FP_mix_o*FQ_mix_o;
350 if (Pr_mix < Pvp_mix/Pc_mix) {
351 double alpha = 3.262 + 14.98*pow(Pr_mix,5.508);
352 double beta = 1.390 + 5.746*Pr_mix;
353 Z2m = 0.600 + 0.760*pow(Pr_mix,alpha) + (0.6990*pow(Pr_mix,beta) -
356 throw CanteraError(
"HighPressureGasTransport::viscosity",
357 "State is outside the limits of the Lucas model, Tr <= 1");
359 }
else if ((Tr_mix > 1.0) && (Tr_mix < 40.0)) {
360 if ((Pr_mix > 0.0) && (Pr_mix <= 100.0)) {
361 double a_fac = 0.001245*exp(5.1726*pow(Tr_mix,-0.3286))/Tr_mix;
362 double b_fac = a_fac*(1.6553*Tr_mix - 1.2723);
363 double c_fac = 0.4489*exp(3.0578*pow(Tr_mix,-37.7332))/Tr_mix;
364 double d_fac = 1.7368*exp(2.2310*pow(Tr_mix,-7.6351))/Tr_mix;
365 double f_fac = 0.9425*exp(-0.1853*pow(Tr_mix,0.4489));
367 Z2m = Z1m*(1 + a_fac*pow(Pr_mix,1.3088)/(b_fac*pow(Pr_mix,f_fac)
368 + pow(1+c_fac*pow(Pr_mix,d_fac),-1)));
370 throw CanteraError(
"HighPressureGasTransport::viscosity",
371 "State is outside the limits of the Lucas model, 1.0 < Tr < 40");
374 throw CanteraError(
"HighPressureGasTransport::viscosity",
375 "State is outside the limits of the Lucas model, Tr > 40");
382 return Z2m*(1 + (FP_mix_o - 1)*pow(Y,-3))*(1 + (FQ_mix_o - 1)
383 *(1/Y - 0.007*pow(log(Y),4)))/(ksi*FP_mix_o*FQ_mix_o);
387double HighPressureGasTransport::Tcrit_i(
size_t i)
398double HighPressureGasTransport::Pcrit_i(
size_t i)
409double HighPressureGasTransport::Vcrit_i(
size_t i)
420double HighPressureGasTransport::Zcrit_i(
size_t i)
431vector<double> HighPressureGasTransport::store(
size_t i,
size_t nsp)
433 vector<double> molefracs(nsp);
435 vector<double> mf_temp(nsp, 0.0);
443double HighPressureGasTransport::FQ_i(
double Q,
double Tr,
double MW)
445 return 1.22*pow(Q,0.15)*(1 + 0.00385*pow(pow(Tr - 12.,2.),1./MW)
446 *fabs(Tr-12)/(Tr-12));
451double HighPressureGasTransport::setPcorr(
double Pr,
double Tr)
453 const static double Pr_lookup[17] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0,
454 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 4.0, 5.0};
455 const static double DP_Rt_lookup[17] = {1.01, 1.01, 1.01, 1.01, 1.01, 1.01,
456 1.01, 1.02, 1.02, 1.02, 1.02, 1.03, 1.03, 1.04, 1.05, 1.06, 1.07};
457 const static double A_ij_lookup[17] = {0.038042, 0.067433, 0.098317,
458 0.137610, 0.175081, 0.216376, 0.314051, 0.385736, 0.514553, 0.599184,
459 0.557725, 0.593007, 0.696001, 0.790770, 0.502100, 0.837452, 0.890390};
460 const static double B_ij_lookup[17] = {1.52267, 2.16794, 2.42910, 2.77605,
461 2.98256, 3.11384, 3.50264, 3.07773, 3.54744, 3.61216, 3.41882, 3.18415,
462 3.37660, 3.27984, 3.39031, 3.23513, 3.13001};
463 const static double C_ij_lookup[17] = {0., 0., 0., 0., 0., 0., 0., 0.141211,
464 0.278407, 0.372683, 0.504894, 0.678469, 0.665702, 0., 0.602907, 0., 0.};
465 const static double E_ij_lookup[17] = {1., 1., 1., 1., 1., 1., 1., 13.45454,
466 14., 10.00900, 8.57519, 10.37483, 11.21674, 1., 6.19043, 1., 1.};
473 frac = (Pr - Pr_lookup[0])/(Pr_lookup[1] - Pr_lookup[0]);
475 for (
int j = 1; j < 17; j++) {
476 if (Pr_lookup[j] > Pr) {
477 frac = (Pr - Pr_lookup[j-1])/(Pr_lookup[j] - Pr_lookup[j-1]);
489 double P_corr_1 = DP_Rt_lookup[Pr_i]*(1.0 - A_ij_lookup[Pr_i]
490 *pow(Tr,-B_ij_lookup[Pr_i]))*(1-C_ij_lookup[Pr_i]
491 *pow(Tr,-E_ij_lookup[Pr_i]));
492 double P_corr_2 = DP_Rt_lookup[Pr_i+1]*(1.0 - A_ij_lookup[Pr_i+1]
493 *pow(Tr,-B_ij_lookup[Pr_i+1]))*(1-C_ij_lookup[Pr_i+1]
494 *pow(Tr,-E_ij_lookup[Pr_i+1]));
495 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.
vector< double > m_mw
Local copy of the species molecular weights.
double m_temp
Current value of the temperature at which the properties in this object are calculated (Kelvin).
virtual void updateDiff_T()
Update the binary diffusion coefficients.
bool m_bindiff_ok
Update boolean for the binary diffusivities at unit pressure.
DenseMatrix m_bdiff
Matrix of binary diffusion coefficients at the reference pressure and the current temperature Size is...
DenseMatrix m_dipole
The effective dipole moment for (i,j) collisions.
vector< double > m_w_ac
Pitzer acentric factor.
void getBinaryDiffCoeffs(const size_t ld, double *const d) override
Returns the matrix of binary diffusion coefficients.
double thermalConductivity() override
Returns the mixture thermal conductivity in W/m/K.
HighPressureGasTransport(ThermoPhase *thermo=0)
default constructor
double viscosity() override
Viscosity of the mixture (kg /m /s)
void getThermalDiffCoeffs(double *const dt) override
Return the thermal diffusion coefficients (kg/m/s)
void getMultiDiffCoeffs(const size_t ld, double *const d) override
Return the Multicomponent diffusion coefficients. Units: [m^2/s].
Class MultiTransport implements multicomponent transport properties for ideal gas mixtures.
bool m_l0000_ok
Boolean indicating viscosity is up to date.
void update_T() override
Update basic temperature-dependent quantities if the temperature has changed.
void eval_L0000(const double *const x)
Evaluate the L0000 matrices.
void updateThermal_T()
Update the temperature-dependent terms needed to compute the thermal conductivity and thermal diffusi...
void update_C() override
Update basic concentration-dependent quantities if the concentrations have changed.
An error indicating that an unimplemented function has been called.
virtual void setMoleFractions(const double *const x)
Set the mole fractions to the specified values.
size_t nSpecies() const
Returns the number of species in the phase.
double temperature() const
Temperature (K).
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
const vector< string > & speciesNames() const
Return a const reference to the vector of species names.
virtual double molarVolume() const
Molar volume (m^3/kmol).
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Base class for a phase with thermodynamic properties.
virtual double critTemperature() const
Critical temperature (K).
virtual void getCp_R_ref(double *cprt) const
Returns the vector of nondimensional constant pressure heat capacities of the reference state at the ...
virtual double critPressure() const
Critical pressure (Pa).
virtual double critVolume() const
Critical volume (m3/kmol).
virtual double satPressure(double t)
Return the saturation pressure given the temperature.
virtual double critCompressibility() const
Critical compressibility (unitless).
ThermoPhase * m_thermo
pointer to the object representing the phase
size_t m_nsp
Number of species.
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...