Cantera  3.1.0a1
IonGasTransport.cpp
Go to the documentation of this file.
1 //! @file IonGasTransport.cpp
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5 
10 #include "cantera/base/utilities.h"
11 #include "cantera/base/global.h"
12 #include "MMCollisionInt.h"
13 
14 namespace Cantera
15 {
16 
17 void IonGasTransport::init(ThermoPhase* thermo, int mode, int log_level)
18 {
19  m_thermo = thermo;
20  m_nsp = m_thermo->nSpecies();
21  m_mode = mode;
22  if (m_mode == CK_Mode) {
23  throw CanteraError("IonGasTransport::init",
24  "mode = CK_Mode, which is an outdated lower-order fit.");
25  }
26  m_log_level = log_level;
27  // make a local copy of species charge
28  for (size_t k = 0; k < m_nsp; k++) {
29  m_speciesCharge.push_back(m_thermo->charge(k));
30  }
31 
32  // Find the index of electron
33  size_t nElectronSpecies = 0;
34  for (size_t k = 0; k < m_nsp; k++) {
35  if (m_thermo->molecularWeight(k) == getElementWeight("E") &&
36  m_thermo->charge(k) == -1)
37  {
38  m_kElectron = k;
39  nElectronSpecies++;
40  }
41  }
42  if (nElectronSpecies > 1) {
43  throw CanteraError("IonGasTransport::init",
44  "Found multiple electron species.");
45  }
46 
47  // Find indices for charge of species
48  for (size_t k = 0; k < m_nsp; k++) {
49  if (m_speciesCharge[k] != 0){
50  if (k != m_kElectron) {
51  m_kIon.push_back(k);
52  }
53  } else {
54  m_kNeutral.push_back(k);
55  }
56  }
57  // set up O2/O2- collision integral [A^2]
58  // Data taken from Prager (2005)
59  const vector<double> temp{300.0, 400.0, 500.0, 600.0, 800.0, 1000.0,
60  1200.0, 1500.0, 2000.0, 2500.0, 3000.0, 4000.0};
61  const vector<double> om11_O2{120.0, 107.0, 98.1, 92.1, 83.0, 77.0,
62  72.6, 67.9, 62.7, 59.3, 56.7, 53.8};
63  vector<double> w(temp.size(),-1);
64  int degree = 5;
65  m_om11_O2.resize(degree + 1);
66  polyfit(temp.size(), degree, temp.data(), om11_O2.data(),
67  w.data(), m_om11_O2.data());
68  // set up Monchick and Mason parameters
70  // set up n64 parameters
71  setupN64();
72  // setup collision integrals
74  m_molefracs.resize(m_nsp);
75  m_spwork.resize(m_nsp);
76  m_visc.resize(m_nsp);
77  m_sqvisc.resize(m_nsp);
78  m_phi.resize(m_nsp, m_nsp, 0.0);
80  m_cond.resize(m_nsp);
81 
82  // make a local copy of the molecular weights
84 
85  m_wratjk.resize(m_nsp, m_nsp, 0.0);
86  m_wratkj1.resize(m_nsp, m_nsp, 0.0);
87  for (size_t j = 0; j < m_nsp; j++) {
88  for (size_t k = j; k < m_nsp; k++) {
89  m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
90  m_wratjk(k,j) = sqrt(m_wratjk(j,k));
91  m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
92  }
93  }
94 }
95 
97 {
98  update_T();
99  update_C();
100 
101  if (m_visc_ok) {
102  return m_viscmix;
103  }
104 
105  double vismix = 0.0;
106  // update m_visc and m_phi if necessary
107  if (!m_viscwt_ok) {
109  }
110 
111  multiply(m_phi, m_molefracs.data(), m_spwork.data());
112 
113  for (size_t k : m_kNeutral) {
114  vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
115  }
116  m_viscmix = vismix;
117  return vismix;
118 }
119 
121 {
122  update_T();
123  update_C();
124  if (!m_spcond_ok) {
125  updateCond_T();
126  }
127  if (!m_condmix_ok) {
128  double sum1 = 0.0, sum2 = 0.0;
129  for (size_t k : m_kNeutral) {
130  sum1 += m_molefracs[k] * m_cond[k];
131  sum2 += m_molefracs[k] / m_cond[k];
132  }
133  m_lambda = 0.5*(sum1 + 1.0/sum2);
134  m_condmix_ok = true;
135  }
136  return m_lambda;
137 }
138 
140 {
141  vector<double> mobi(m_nsp);
142  getMobilities(&mobi[0]);
143  double p = m_thermo->pressure();
144  double sum = 0.0;
145  for (size_t k : m_kIon) {
146  double ND_k = m_molefracs[k] * p / m_kbt;
147  sum += ND_k * std::abs(m_speciesCharge[k]) * ElectronCharge * mobi[k];
148  }
149  if (m_kElectron != npos) {
150  sum += m_molefracs[m_kElectron] * p / m_kbt *
151  ElectronCharge * mobi[m_kElectron];
152  }
153  return sum;
154 }
155 
157 {
158  GasTransport::fitDiffCoeffs(integrals);
159 
160  // number of points to use in generating fit data
161  const size_t np = 50;
162  int degree = 4;
163  double dt = (m_thermo->maxTemp() - m_thermo->minTemp())/(np-1);
164  vector<double> tlog(np);
165  vector<double> w(np);
166 
167  // generate array of log(t) values
168  for (size_t n = 0; n < np; n++) {
169  double t = m_thermo->minTemp() + dt*n;
170  tlog[n] = log(t);
171  }
172 
173  // vector of polynomial coefficients
174  vector<double> c(degree + 1);
175  double err = 0.0, relerr = 0.0,
176  mxerr = 0.0, mxrelerr = 0.0;
177 
178  vector<double> diff(np + 1);
179  // The array order still not ideal
180  for (size_t k = 0; k < m_nsp; k++) {
181  for (size_t j = k; j < m_nsp; j++) {
182  if (m_alpha[k] == 0.0 || m_alpha[j] == 0.0 ||
183  k == m_kElectron || j == m_kElectron) {
184  continue;
185  }
186  if (m_speciesCharge[k] == 0) {
187  if (m_speciesCharge[j] == 0) {
188  continue;
189  }
190  } else {
191  if (m_speciesCharge[j] != 0) {
192  continue;
193  }
194  }
195  for (size_t n = 0; n < np; n++) {
196  double t = m_thermo->minTemp() + dt*n;
197  double eps = m_epsilon(j,k);
198  double tstar = Boltzmann * t/eps;
199  double sigma = m_diam(j,k);
200  double om11 = omega11_n64(tstar, m_gamma(j,k))
201  * Pi * sigma * sigma;
202  // Stockmayer-(n,6,4) model is not suitable for collision
203  // between O2/O2- due to resonant charge transfer.
204  // Therefore, the experimental collision data is used instead.
205  if ((k == m_thermo->speciesIndex("O2-") ||
206  j == m_thermo->speciesIndex("O2-")) &&
207  (k == m_thermo->speciesIndex("O2") ||
208  j == m_thermo->speciesIndex("O2"))) {
209  om11 = poly5(t, m_om11_O2.data()) / 1e20;
210  }
211  double diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,j))
212  * pow(Boltzmann * t, 1.5) / om11;
213 
214  diff[n] = diffcoeff/pow(t, 1.5);
215  w[n] = 1.0/(diff[n]*diff[n]);
216  }
217  polyfit(np, degree, tlog.data(), diff.data(), w.data(), c.data());
218 
219  for (size_t n = 0; n < np; n++) {
220  double val, fit;
221  double t = exp(tlog[n]);
222  double pre = pow(t, 1.5);
223  val = pre * diff[n];
224  fit = pre * poly4(tlog[n], c.data());
225  err = fit - val;
226  relerr = err/val;
227  mxerr = std::max(mxerr, fabs(err));
228  mxrelerr = std::max(mxrelerr, fabs(relerr));
229  }
230  size_t sum = k * (k + 1) / 2;
231  m_diffcoeffs[k*m_nsp+j-sum] = c;
232  if (m_log_level >= 2) {
233  writelog(m_thermo->speciesName(k) + "__" +
234  m_thermo->speciesName(j) + ": [" + vec2str(c) + "]\n");
235  }
236  }
237  }
238 
239  if (m_log_level) {
240  writelogf("Maximum binary diffusion coefficient absolute error:"
241  " %12.6g\n", mxerr);
242  writelogf("Maximum binary diffusion coefficient relative error:"
243  "%12.6g", mxrelerr);
244  }
245 }
246 
248 {
249  m_gamma.resize(m_nsp, m_nsp, 0.0);
250  for (size_t i : m_kIon) {
251  for (size_t j : m_kNeutral) {
252  if (m_alpha[j] != 0.0 && m_alpha[i] != 0.0) {
253  double r_alpha = m_alpha[i] / m_alpha[j];
254  // save a copy of polarizability in Angstrom
255  double alphaA_i = m_alpha[i] * 1e30;
256  double alphaA_j = m_alpha[j] * 1e30;
257  // The ratio of dispersion to induction forces
258  double xi = alphaA_i / (m_speciesCharge[i] * m_speciesCharge[i] *
259  (1.0 + pow(2 * r_alpha, 2./3.)) * sqrt(alphaA_j));
260 
261  // the collision diameter
262  double K1 = 1.767;
263  double kappa = 0.095;
264  m_diam(i,j) = K1 * (pow(m_alpha[i], 1./3.) + pow(m_alpha[j], 1./3.)) /
265  pow(alphaA_i * alphaA_j * (1.0 + 1.0 / xi), kappa);
266 
267  // The original K2 is 0.72, but Han et al. suggested that K2 = 1.44
268  // for better fit.
269  double K2 = 1.44;
270  double epsilon = K2 * ElectronCharge * ElectronCharge *
272  m_alpha[j] * (1.0 + xi) /
273  (8 * Pi * epsilon_0 * pow(m_diam(i,j),4));
274  if (epsilon != 0.0) {
275  m_epsilon(i,j) = epsilon;
276  }
277 
278  // Calculate dispersion coefficient and quadrupole polarizability
279  // from curve fitting if not available.
280  // Neutrals
281  if (m_disp[j] == 0.0) {
282  m_disp[j] = exp(1.8846*log(alphaA_j)-0.4737)* 1e-50;
283  }
284  if (m_quad_polar[j] == 0.0) {
285  m_quad_polar[j] = 2.0 * m_disp[j];
286  }
287  // Ions
288  if (m_disp[i] == 0.0) {
289  if (m_speciesCharge[i] > 0) {
290  m_disp[i] = exp(1.8853*log(alphaA_i)+0.2682)* 1e-50;
291  } else {
292  m_disp[i] = exp(3.2246*log(alphaA_i)-3.2397)* 1e-50;
293  }
294  }
295 
296  // The binary dispersion coefficient is determined by the combination rule
297  // Reference:
298  // Tang, K. T. "Dynamic polarizabilities and van der Waals coefficients."
299  // Physical Review 177.1 (1969): 108.
300  double C6 = 2.0 * m_disp[i] * m_disp[j] /
301  (1.0/r_alpha * m_disp[i] + r_alpha * m_disp[j]);
302 
303  m_gamma(i,j) = (2.0 / pow(m_speciesCharge[i],2) * C6 + m_quad_polar[j]) /
304  (m_alpha[j] * m_diam(i,j) * m_diam(i,j));//Dimensionless
305 
306  // properties are symmetric
307  m_diam(j,i) = m_diam(i,j);
308  m_epsilon(j,i) = m_epsilon(i,j);
309  m_gamma(j,i) = m_gamma(i,j);
310  }
311  }
312  }
313 }
314 
315 double IonGasTransport::omega11_n64(const double tstar, const double gamma)
316 {
317  double logtstar = log(tstar);
318  double om11 = 0.0;
319  if (tstar < 0.01) {
320  throw CanteraError("IonGasTransport::omega11_n64",
321  "tstar = {} is smaller than 0.01", tstar);
322  } else if (tstar <= 0.04) {
323  // for interval 0.01 to 0.04, SSE = 0.006; R^2 = 1; RMSE = 0.020
324  om11 = 2.97 - 12.0 * gamma
325  - 0.887 * logtstar
326  + 3.86 * gamma * gamma
327  - 6.45 * gamma * logtstar
328  - 0.275 * logtstar * logtstar
329  + 1.20 * gamma * gamma * logtstar
330  - 1.24 * gamma * logtstar * logtstar
331  - 0.164 * pow(logtstar,3);
332  } else if (tstar <= 1000) {
333  // for interval 0.04 to 1000, SSE = 0.282; R^2 = 1; RMSE = 0.033
334  om11 = 1.22 - 0.0343 * gamma
335  + (-0.769 + 0.232 * gamma) * logtstar
336  + (0.306 - 0.165 * gamma) * logtstar * logtstar
337  + (-0.0465 + 0.0388 * gamma) * pow(logtstar,3)
338  + (0.000614 - 0.00285 * gamma) * pow(logtstar,4)
339  + 0.000238 * pow(logtstar,5);
340  } else {
341  throw CanteraError("IonGasTransport::omega11_n64",
342  "tstar = {} is larger than 1000", tstar);
343  }
344  return om11;
345 }
346 
348 {
349  update_T();
350  update_C();
351 
352  // update the binary diffusion coefficients if necessary
353  if (!m_bindiff_ok) {
354  updateDiff_T();
355  }
356 
357  double mmw = m_thermo->meanMolecularWeight();
358  double p = m_thermo->pressure();
359  if (m_nsp == 1) {
360  d[0] = m_bdiff(0,0) / p;
361  } else {
362  for (size_t k = 0; k < m_nsp; k++) {
363  if (k == m_kElectron) {
364  d[k] = 0.4 * m_kbt / ElectronCharge;
365  } else {
366  double sum2 = 0.0;
367  for (size_t j : m_kNeutral) {
368  if (j != k) {
369  sum2 += m_molefracs[j] / m_bdiff(j,k);
370  }
371  }
372  if (sum2 <= 0.0) {
373  d[k] = m_bdiff(k,k) / p;
374  } else {
375  d[k] = (mmw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2);
376  }
377  }
378  }
379  }
380 }
381 
382 void IonGasTransport::getMobilities(double* const mobi)
383 {
384  update_T();
385  update_C();
386 
387  // update the binary diffusion coefficients if necessary
388  if (!m_bindiff_ok) {
389  updateDiff_T();
390  }
391  double p = m_thermo->pressure();
392  for (size_t k = 0; k < m_nsp; k++) {
393  if (k == m_kElectron) {
394  mobi[k] = 0.4;
395  } else {
396  mobi[k] = 0.0;
397  }
398  }
399  for (size_t k : m_kIon) {
400  double sum = 0.0;
401  for (size_t j : m_kNeutral) {
402  double bmobi = m_bdiff(k,j) * ElectronCharge / m_kbt;
403  sum += m_molefracs[j] / bmobi;
404  }
405  mobi[k] = 1.0 / sum / p;
406  }
407 }
408 
409 }
Monchick and Mason collision integrals.
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
Definition: DenseMatrix.cpp:60
virtual void setupCollisionParameters()
Setup parameters for a new kinetic-theory-based transport manager for low-density gases.
vector< double > m_mw
Local copy of the species molecular weights.
Definition: GasTransport.h:334
vector< double > m_molefracs
Vector of species mole fractions.
Definition: GasTransport.h:297
vector< double > m_quad_polar
Quadrupole polarizability.
Definition: GasTransport.h:542
bool m_visc_ok
Update boolean for mixture rule for the mixture viscosity.
Definition: GasTransport.h:303
DenseMatrix m_wratkj1
Holds square roots of molecular weight ratios.
Definition: GasTransport.h:349
virtual void fitDiffCoeffs(MMCollisionInt &integrals)
Generate polynomial fits to the binary diffusion coefficients.
vector< double > m_disp
Dispersion coefficient.
Definition: GasTransport.h:539
virtual void updateDiff_T()
Update the binary diffusion coefficients.
vector< double > m_sqvisc
vector of square root of species viscosities sqrt(kg /m /s).
Definition: GasTransport.h:354
DenseMatrix m_wratjk
Holds square roots of molecular weight ratios.
Definition: GasTransport.h:343
bool m_bindiff_ok
Update boolean for the binary diffusivities at unit pressure.
Definition: GasTransport.h:312
DenseMatrix m_epsilon
The effective well depth for (i,j) collisions.
Definition: GasTransport.h:509
DenseMatrix m_diam
hard-sphere diameter for (i,j) collision
Definition: GasTransport.h:500
vector< double > m_spwork
work space length = m_kk
Definition: GasTransport.h:322
int m_mode
Type of the polynomial fits to temperature.
Definition: GasTransport.h:316
virtual void updateViscosity_T()
Update the temperature-dependent viscosity terms.
double m_viscmix
Internal storage for the viscosity of the mixture (kg /m /s)
Definition: GasTransport.h:300
vector< double > m_visc
vector of species viscosities (kg /m /s).
Definition: GasTransport.h:326
DenseMatrix m_bdiff
Matrix of binary diffusion coefficients at the reference pressure and the current temperature Size is...
Definition: GasTransport.h:392
vector< vector< double > > m_diffcoeffs
Polynomial fits to the binary diffusivity of each species.
Definition: GasTransport.h:388
double m_kbt
Current value of Boltzmann constant times the temperature (Joules)
Definition: GasTransport.h:364
DenseMatrix m_reducedMass
This is the reduced mass of the interaction between species i and j.
Definition: GasTransport.h:491
int m_log_level
Level of verbose printing during initialization.
Definition: GasTransport.h:545
bool m_viscwt_ok
Update boolean for the weighting factors for the mixture viscosity.
Definition: GasTransport.h:306
vector< double > m_alpha
Polarizability of each species in the phase.
Definition: GasTransport.h:468
DenseMatrix m_phi
m_phi is a Viscosity Weighting Function. size = m_nsp * n_nsp
Definition: GasTransport.h:319
void setupCollisionIntegral()
Setup range for polynomial fits to collision integrals of Monchick & Mason .
void fitDiffCoeffs(MMCollisionInt &integrals) override
Generate polynomial fits to the binary diffusion coefficients.
void getMobilities(double *const mobi) override
The mobilities for ions in gas.
void init(ThermoPhase *thermo, int mode, int log_level) override
Initialize a transport manager.
size_t m_kElectron
index of electron
double thermalConductivity() override
Returns the mixture thermal conductivity (W/m/K).
void getMixDiffCoeffs(double *const d) override
The mixture transport for ionized gas.
vector< double > m_om11_O2
polynomial of the collision integral for O2/O2-
double omega11_n64(const double tstar, const double gamma)
Collision integral of omega11 of n64 collision model.
DenseMatrix m_gamma
parameter of omega11 of n64
void setupN64()
setup parameters for n64 model
double viscosity() override
Viscosity of the mixture (kg/m/s).
vector< size_t > m_kIon
index of ions (exclude electron.)
double electricalConductivity() override
The electrical conductivity (Siemens/m).
vector< size_t > m_kNeutral
index of neutral species
vector< double > m_speciesCharge
electrical properties
Calculation of Collision integrals.
void update_T() override
Update the internal parameters whenever the temperature has changed.
bool m_spcond_ok
Update boolean for the species thermal conductivities.
Definition: MixTransport.h:175
double m_lambda
Internal storage for the calculated mixture thermal conductivity.
Definition: MixTransport.h:172
void updateCond_T()
Update the temperature dependent parts of the species thermal conductivities.
void update_C() override
Update the internal parameters whenever the concentrations have changed.
vector< double > m_cond
vector of species thermal conductivities (W/m /K)
Definition: MixTransport.h:166
bool m_condmix_ok
Update boolean for the mixture rule for the mixture thermal conductivity.
Definition: MixTransport.h:178
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:231
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:655
string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:142
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:395
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
Definition: Phase.cpp:129
double molecularWeight(size_t k) const
Molecular weight of species k.
Definition: Phase.cpp:383
virtual double pressure() const
Return the thermodynamic pressure (Pa).
Definition: Phase.h:580
double charge(size_t k) const
Dimensionless electrical charge of a single molecule of species k The charge is normalized by the the...
Definition: Phase.h:538
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:390
virtual double minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid.
Definition: ThermoPhase.h:451
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
Definition: ThermoPhase.h:504
ThermoPhase * m_thermo
pointer to the object representing the phase
Definition: Transport.h:420
size_t m_nsp
Number of species.
Definition: Transport.h:423
ThermoPhase & thermo()
Phase object.
Definition: Transport.h:103
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
string vec2str(const vector< double > &v, const string &fmt, const string &sep)
Convert a vector to a string (separated by commas)
Definition: stringUtils.cpp:33
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:195
void writelog(const string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:175
R poly5(D x, R *c)
Templated evaluation of a polynomial of order 5.
Definition: utilities.h:141
R poly4(D x, R *c)
Evaluates a polynomial of order 4.
Definition: utilities.h:153
double polyfit(size_t n, size_t deg, const double *xp, const double *yp, const double *wp, double *pp)
Fits a polynomial function to a set of data points.
Definition: polyfit.cpp:14
const double Boltzmann
Boltzmann constant [J/K].
Definition: ct_defs.h:84
const double epsilon_0
Permittivity of free space [F/m].
Definition: ct_defs.h:137
const double ElectronCharge
Elementary charge [C].
Definition: ct_defs.h:90
const double Pi
Pi.
Definition: ct_defs.h:68
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:180
double getElementWeight(const string &ename)
Get the atomic weight of an element.
Definition: Elements.cpp:251
void multiply(const DenseMatrix &A, const double *const b, double *const prod)
Multiply A*b and return the result in prod. Uses BLAS routine DGEMV.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...