Cantera  2.4.0
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 http://www.cantera.org/license.txt for license and copyright information.
5 
9 #include "MMCollisionInt.h"
10 
11 namespace Cantera
12 {
13 IonGasTransport::IonGasTransport() :
14  m_kElectron(npos)
15 {
16 }
17 
18 void IonGasTransport::init(thermo_t* thermo, int mode, int log_level)
19 {
20  m_thermo = thermo;
21  m_nsp = m_thermo->nSpecies();
22  m_mode = mode;
23  if (m_mode == CK_Mode) {
24  throw CanteraError("IonGasTransport::init(thermo, mode, log_level)",
25  "mode = CK_Mode, which is an outdated lower-order fit.");
26  }
27  m_log_level = log_level;
28  // make a local copy of species charge
29  for (size_t k = 0; k < m_nsp; k++) {
30  m_speciesCharge.push_back(m_thermo->charge(k));
31  }
32 
33  // Find the index of electron
34  if (m_thermo->speciesIndex("E") != npos ) {
36  }
37 
38  // Find indices for charge of species
39  for (size_t k = 0; k < m_nsp; k++) {
40  if (m_speciesCharge[k] != 0){
41  if (k != m_kElectron) {
42  m_kIon.push_back(k);
43  }
44  } else {
45  m_kNeutral.push_back(k);
46  }
47  }
48  // set up Monchick and Mason parameters
50  // set up n64 parameters
51  setupN64();
52  // setup collision integrals
54  m_molefracs.resize(m_nsp);
55  m_spwork.resize(m_nsp);
56  m_visc.resize(m_nsp);
57  m_sqvisc.resize(m_nsp);
58  m_phi.resize(m_nsp, m_nsp, 0.0);
60  m_cond.resize(m_nsp);
61 
62  // make a local copy of the molecular weights
64 
65  m_wratjk.resize(m_nsp, m_nsp, 0.0);
66  m_wratkj1.resize(m_nsp, m_nsp, 0.0);
67  for (size_t j = 0; j < m_nsp; j++) {
68  for (size_t k = j; k < m_nsp; k++) {
69  m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
70  m_wratjk(k,j) = sqrt(m_wratjk(j,k));
71  m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
72  }
73  }
74 }
75 
77 {
78  update_T();
79  update_C();
80 
81  if (m_visc_ok) {
82  return m_viscmix;
83  }
84 
85  double vismix = 0.0;
86  // update m_visc and m_phi if necessary
87  if (!m_viscwt_ok) {
89  }
90 
91  multiply(m_phi, m_molefracs.data(), m_spwork.data());
92 
93  for (size_t k : m_kNeutral) {
94  vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
95  }
96  m_viscmix = vismix;
97  return vismix;
98 }
99 
101 {
102  update_T();
103  update_C();
104  if (!m_spcond_ok) {
105  updateCond_T();
106  }
107  if (!m_condmix_ok) {
108  doublereal sum1 = 0.0, sum2 = 0.0;
109  for (size_t k : m_kNeutral) {
110  sum1 += m_molefracs[k] * m_cond[k];
111  sum2 += m_molefracs[k] / m_cond[k];
112  }
113  m_lambda = 0.5*(sum1 + 1.0/sum2);
114  m_condmix_ok = true;
115  }
116  return m_lambda;
117 }
118 
120 {
121  GasTransport::fitDiffCoeffs(integrals);
122 
123  // number of points to use in generating fit data
124  const size_t np = 50;
125  int degree = 4;
126  double dt = (m_thermo->maxTemp() - m_thermo->minTemp())/(np-1);
127  vector_fp tlog(np);
128  vector_fp w(np);
129 
130  // generate array of log(t) values
131  for (size_t n = 0; n < np; n++) {
132  double t = m_thermo->minTemp() + dt*n;
133  tlog[n] = log(t);
134  }
135 
136  // vector of polynomial coefficients
137  vector_fp c(degree + 1);
138  double err = 0.0, relerr = 0.0,
139  mxerr = 0.0, mxrelerr = 0.0;
140 
141  vector_fp diff(np + 1);
142  // The array order still not ideal
143  for (size_t k = 0; k < m_nsp; k++) {
144  for (size_t j = k; j < m_nsp; j++) {
145  if (m_alpha[k] == 0.0 || m_alpha[j] == 0.0 ||
146  k == m_kElectron || j == m_kElectron) {
147  continue;
148  }
149  if (m_speciesCharge[k] == 0) {
150  if (m_speciesCharge[j] == 0) {
151  continue;
152  }
153  } else {
154  if (m_speciesCharge[j] != 0) {
155  continue;
156  }
157  }
158  for (size_t n = 0; n < np; n++) {
159  double t = m_thermo->minTemp() + dt*n;
160  double eps = m_epsilon(j,k);
161  double tstar = Boltzmann * t/eps;
162  double sigma = m_diam(j,k);
163  double om11 = omega11_n64(tstar, m_gamma(j,k));
164  double diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,j))
165  * pow(Boltzmann * t, 1.5) / (Pi * sigma * sigma * om11);
166 
167  diff[n] = diffcoeff/pow(t, 1.5);
168  w[n] = 1.0/(diff[n]*diff[n]);
169  }
170  polyfit(np, degree, tlog.data(), diff.data(), w.data(), c.data());
171 
172  for (size_t n = 0; n < np; n++) {
173  double val, fit;
174  double t = exp(tlog[n]);
175  double pre = pow(t, 1.5);
176  val = pre * diff[n];
177  fit = pre * poly4(tlog[n], c.data());
178  err = fit - val;
179  relerr = err/val;
180  mxerr = std::max(mxerr, fabs(err));
181  mxrelerr = std::max(mxrelerr, fabs(relerr));
182  }
183  size_t sum = k * (k + 1) / 2;
184  m_diffcoeffs[k*m_nsp+j-sum] = c;
185  if (m_log_level >= 2) {
186  writelog(m_thermo->speciesName(k) + "__" +
187  m_thermo->speciesName(j) + ": [" + vec2str(c) + "]\n");
188  }
189  }
190  }
191 
192  if (m_log_level) {
193  writelogf("Maximum binary diffusion coefficient absolute error:"
194  " %12.6g\n", mxerr);
195  writelogf("Maximum binary diffusion coefficient relative error:"
196  "%12.6g", mxrelerr);
197  }
198 }
199 
201 {
202  m_gamma.resize(m_nsp, m_nsp, 0.0);
203  for (size_t i : m_kIon) {
204  for (size_t j : m_kNeutral) {
205  if (m_alpha[j] != 0.0 && m_alpha[i] != 0.0) {
206  double r_alpha = m_alpha[i] / m_alpha[j];
207  // save a copy of polarizability in Angstrom
208  double alphaA_i = m_alpha[i] * 1e30;
209  double alphaA_j = m_alpha[j] * 1e30;
210  // The ratio of dispersion to induction forces
211  double xi = alphaA_i / (m_speciesCharge[i] * m_speciesCharge[i] *
212  (1.0 + pow(2 * r_alpha, 2./3.)) * sqrt(alphaA_j));
213 
214  // the collision diameter
215  double K1 = 1.767;
216  double kappa = 0.095;
217  m_diam(i,j) = K1 * (pow(m_alpha[i], 1./3.) + pow(m_alpha[j], 1./3.)) /
218  pow(alphaA_i * alphaA_j * (1.0 + 1.0 / xi), kappa);
219 
220  // The original K2 is 0.72, but Han et al. suggested that K2 = 1.44
221  // for better fit.
222  double K2 = 1.44;
223  double epsilon = K2 * ElectronCharge * ElectronCharge *
225  m_alpha[j] * (1.0 + xi) /
226  (8 * Pi * epsilon_0 * pow(m_diam(i,j),4));
227  if (epsilon != 0.0) {
228  m_epsilon(i,j) = epsilon;
229  }
230 
231  // Calculate dipersion coefficient and quadrupole polarizability
232  // from curve fitting if not available.
233  // Neutrals
234  if (m_disp[j] == 0.0) {
235  m_disp[j] = exp(1.8846*log(alphaA_j)-0.4737)* 1e-50;
236  }
237  if (m_quad_polar[j] == 0.0) {
238  m_quad_polar[j] = 2.0 * m_disp[j];
239  }
240  // Ions
241  if (m_disp[i] == 0.0) {
242  if (m_speciesCharge[i] > 0) {
243  m_disp[i] = exp(1.8853*log(alphaA_i)+0.2682)* 1e-50;
244  } else {
245  m_disp[i] = exp(3.2246*log(alphaA_i)-3.2397)* 1e-50;
246  }
247  }
248 
249  // The binary dispersion coefficient is determined by the combination rule
250  // Reference:
251  // Tang, K. T. "Dynamic polarizabilities and van der Waals coefficients."
252  // Physical Review 177.1 (1969): 108.
253  double C6 = 2.0 * m_disp[i] * m_disp[j] /
254  (1.0/r_alpha * m_disp[i] + r_alpha * m_disp[j]);
255 
256  m_gamma(i,j) = (2.0 / pow(m_speciesCharge[i],2) * C6 + m_quad_polar[j]) /
257  (m_alpha[j] * m_diam(i,j) * m_diam(i,j));//Dimensionless
258 
259  // properties are symmetric
260  m_diam(j,i) = m_diam(i,j);
261  m_epsilon(j,i) = m_epsilon(i,j);
262  m_gamma(j,i) = m_gamma(i,j);
263  }
264  }
265  }
266 }
267 
268 double IonGasTransport::omega11_n64(const double tstar, const double gamma)
269 {
270  double logtstar = log(tstar);
271  double om11 = 0.0;
272  if (tstar < 0.01) {
273  throw CanteraError("IonGasTransport::omega11_n64(tstar, gamma)",
274  "tstar = {} is smaller than 0.01", tstar);
275  } else if (tstar <= 0.04) {
276  // for interval 0.01 to 0.04, SSE = 0.006; R^2 = 1; RMSE = 0.020
277  om11 = 2.97 - 12.0 * gamma
278  - 0.887 * logtstar
279  + 3.86 * gamma * gamma
280  - 6.45 * gamma * logtstar
281  - 0.275 * logtstar * logtstar
282  + 1.20 * gamma * gamma * logtstar
283  - 1.24 * gamma * logtstar * logtstar
284  - 0.164 * pow(logtstar,3);
285  } else if (tstar <= 1000) {
286  // for interval 0.04 to 1000, SSE = 0.282; R^2 = 1; RMSE = 0.033
287  om11 = 1.22 - 0.0343 * gamma
288  + (-0.769 + 0.232 * gamma) * logtstar
289  + (0.306 - 0.165 * gamma) * logtstar * logtstar
290  + (-0.0465 + 0.0388 * gamma) * pow(logtstar,3)
291  + (0.000614 - 0.00285 * gamma) * pow(logtstar,4)
292  + 0.000238 * pow(logtstar,5);
293  } else {
294  throw CanteraError("IonGasTransport::omega11_n64(tstar, gamma)",
295  "tstar = {} is larger than 1000", tstar);
296  }
297  return om11;
298 }
299 
301 {
302  update_T();
303  update_C();
304 
305  // update the binary diffusion coefficients if necessary
306  if (!m_bindiff_ok) {
307  updateDiff_T();
308  }
309 
310  double mmw = m_thermo->meanMolecularWeight();
311  double p = m_thermo->pressure();
312  if (m_nsp == 1) {
313  d[0] = m_bdiff(0,0) / p;
314  } else {
315  for (size_t k = 0; k < m_nsp; k++) {
316  if (k == m_kElectron) {
317  d[k] = 0.4 * m_kbt / ElectronCharge;
318  } else {
319  double sum2 = 0.0;
320  for (size_t j : m_kNeutral) {
321  if (j != k) {
322  sum2 += m_molefracs[j] / m_bdiff(j,k);
323  }
324  }
325  if (sum2 <= 0.0) {
326  d[k] = m_bdiff(k,k) / p;
327  } else {
328  d[k] = (mmw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2);
329  }
330  }
331  }
332  }
333 }
334 
335 void IonGasTransport::getMobilities(double* const mobi)
336 {
337  double p = m_thermo->pressure();
338  for (size_t k = 0; k < m_nsp; k++) {
339  if (k == m_kElectron) {
340  mobi[k] = 0.4;
341  } else {
342  mobi[k] = 0.0;
343  }
344  }
345  for (size_t k : m_kIon) {
346  double sum = 0.0;
347  for (size_t j : m_kNeutral) {
348  double bmobi = m_bdiff(k,j) * ElectronCharge / m_kbt;
349  sum += m_molefracs[j] / bmobi;
350  }
351  mobi[k] = 1.0 / sum / p;
352  }
353 }
354 
355 }
356 
virtual void getMixDiffCoeffs(double *const d)
The mixture transport for ionized gas.
const vector_fp & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:437
vector_fp m_cond
vector of species thermal conductivities (W/m /K)
Definition: MixTransport.h:174
bool m_visc_ok
Update boolean for mixture rule for the mixture viscosity.
Definition: GasTransport.h:258
std::string vec2str(const vector_fp &v, const std::string &fmt, const std::string &sep)
Convert a vector to a string (separated by commas)
Definition: stringUtils.cpp:34
std::vector< size_t > m_kIon
index of ions (exclude electron.)
size_t speciesIndex(const std::string &name) const
Returns the index of a species named &#39;name&#39; within the Phase object.
Definition: Phase.cpp:175
std::vector< vector_fp > m_diffcoeffs
Polynomial fits to the binary diffusivity of each species.
Definition: GasTransport.h:350
DenseMatrix m_gamma
parameter of omega11 of n64
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
bool m_condmix_ok
Update boolean for the mixture rule for the mixture thermal conductivity.
Definition: MixTransport.h:186
virtual void fitDiffCoeffs(MMCollisionInt &integrals)
Generate polynomial fits to the binary diffusion coefficients.
virtual void update_C()
Update the internal parameters whenever the concentrations have changed.
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:153
bool m_bindiff_ok
Update boolean for the binary diffusivities at unit pressure.
Definition: GasTransport.h:267
thermo_t * m_thermo
pointer to the object representing the phase
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:266
size_t m_kElectron
index of electron
virtual doublereal minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid...
Definition: ThermoPhase.h:131
DenseMatrix m_wratjk
Holds square roots of molecular weight ratios.
Definition: GasTransport.h:298
DenseMatrix m_bdiff
Matrix of binary diffusion coefficients at the reference pressure and the current temperature Size is...
Definition: GasTransport.h:354
const doublereal Pi
Pi.
Definition: ct_defs.h:51
virtual void updateViscosity_T()
Update the temperature-dependent viscosity terms.
Monk and Monchick collision integrals.
bool m_spcond_ok
Update boolean for the species thermal conductivities.
Definition: MixTransport.h:183
R poly4(D x, R *c)
Evaluates a polynomial of order 4.
Definition: utilities.h:467
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:93
int m_mode
Type of the polynomial fits to temperature.
Definition: GasTransport.h:271
vector_int m_speciesCharge
electrical properties
virtual void update_T()
Update the internal parameters whenever the temperature has changed.
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
Definition: DenseMatrix.cpp:69
void updateCond_T()
Update the temperature dependent parts of the species thermal conductivities.
virtual double thermalConductivity()
Returns the mixture thermal conductivity (W/m/K).
virtual void updateDiff_T()
Update the binary diffusion coefficients.
vector_fp m_visc
vector of species viscosities (kg /m /s).
Definition: GasTransport.h:281
vector_fp m_spwork
work space length = m_kk
Definition: GasTransport.h:277
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:191
virtual double viscosity()
Viscosity of the mixture (kg/m/s).
vector_fp m_disp
Dispersion coefficient.
Definition: GasTransport.h:497
void setupCollisionIntegral()
Setup range for polynomial fits to collision integrals of Monchick & Mason.
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.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:65
std::vector< size_t > m_kNeutral
index of neutral species
double omega11_n64(const double tstar, const double gamma)
bool m_viscwt_ok
Update boolean for the weighting factors for the mixture viscosity.
Definition: GasTransport.h:261
DenseMatrix m_phi
m_phi is a Viscosity Weighting Function. size = m_nsp * n_nsp
Definition: GasTransport.h:274
thermo_t & thermo()
virtual doublereal pressure() const
Return the thermodynamic pressure (Pa).
Definition: ThermoPhase.h:245
virtual void init(thermo_t *thermo, int mode, int log_level)
Initialize a transport manager.
DenseMatrix m_reducedMass
This is the reduced mass of the interaction between species i and j.
Definition: GasTransport.h:449
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:174
virtual void fitDiffCoeffs(MMCollisionInt &integrals)
Generate polynomial fits to the binary diffusion coefficients.
doublereal m_lambda
Internal storage for the calculated mixture thermal conductivity.
Definition: MixTransport.h:180
vector_fp m_sqvisc
vector of square root of species viscosities sqrt(kg /m /s).
Definition: GasTransport.h:309
DenseMatrix m_epsilon
The effective well depth for (i,j) collisions.
Definition: GasTransport.h:467
virtual void getMobilities(double *const mobi)
The mobilities for ions in gas.
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:157
doublereal meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition: Phase.h:661
DenseMatrix m_diam
hard-sphere diameter for (i,j) collision
Definition: GasTransport.h:458
virtual void setupCollisionParameters()
Setup parameters for a new kinetic-theory-based transport manager for low-density gases...
Contains declarations for string manipulation functions within Cantera.
size_t m_nsp
Number of species.
DenseMatrix m_wratkj1
Holds square roots of molecular weight ratios.
Definition: GasTransport.h:304
const doublereal epsilon_0
Permittivity of free space in F/m.
Definition: ct_defs.h:106
Calculation of Collision integrals.
doublereal m_viscmix
Internal storage for the viscosity of the mixture (kg /m /s)
Definition: GasTransport.h:255
doublereal m_kbt
Current value of Boltzmann constant times the temperature (Joules)
Definition: GasTransport.h:319
void setupN64()
setup parameters for n64 model
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
virtual doublereal maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
Definition: ThermoPhase.h:184
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
int m_log_level
Level of verbose printing during initialization.
Definition: GasTransport.h:503
doublereal 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:577
const doublereal Boltzmann
Boltzmann&#39;s constant [J/K].
Definition: ct_defs.h:76
vector_fp m_quad_polar
Quadrupole polarizability.
Definition: GasTransport.h:500
vector_fp m_molefracs
Vector of species mole fractions.
Definition: GasTransport.h:252
vector_fp m_alpha
Polarizability of each species in the phase.
Definition: GasTransport.h:426
vector_fp m_mw
Local copy of the species molecular weights.
Definition: GasTransport.h:289