Loading [MathJax]/extensions/tex2jax.js
Cantera  3.2.0a1
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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
11#include "cantera/base/global.h"
12#include "MMCollisionInt.h"
13
14namespace Cantera
15{
16
17void IonGasTransport::init(ThermoPhase* thermo, int mode)
18{
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 // make a local copy of species charge
27 for (size_t k = 0; k < m_nsp; k++) {
28 m_speciesCharge.push_back(m_thermo->charge(k));
29 }
30
31 // Find the index of electron
32 size_t nElectronSpecies = 0;
33 for (size_t k = 0; k < m_nsp; k++) {
35 m_thermo->charge(k) == -1)
36 {
37 m_kElectron = k;
38 nElectronSpecies++;
39 }
40 }
41 if (nElectronSpecies > 1) {
42 throw CanteraError("IonGasTransport::init",
43 "Found multiple electron species.");
44 }
45
46 // Find indices for charge of species
47 for (size_t k = 0; k < m_nsp; k++) {
48 if (m_speciesCharge[k] != 0){
49 if (k != m_kElectron) {
50 m_kIon.push_back(k);
51 }
52 } else {
53 m_kNeutral.push_back(k);
54 }
55 }
56 // set up O2/O2- collision integral [A^2]
57 // Data taken from Prager (2005)
58 const vector<double> temp{300.0, 400.0, 500.0, 600.0, 800.0, 1000.0,
59 1200.0, 1500.0, 2000.0, 2500.0, 3000.0, 4000.0};
60 const vector<double> om11_O2{120.0, 107.0, 98.1, 92.1, 83.0, 77.0,
61 72.6, 67.9, 62.7, 59.3, 56.7, 53.8};
62 vector<double> w(temp.size(),-1);
63 int degree = 5;
64 m_om11_O2.resize(degree + 1);
65 polyfit(temp.size(), degree, temp.data(), om11_O2.data(),
66 w.data(), m_om11_O2.data());
67 // set up Monchick and Mason parameters
69 // set up n64 parameters
70 setupN64();
71 // setup collision integrals
73 m_molefracs.resize(m_nsp);
74 m_spwork.resize(m_nsp);
75 m_visc.resize(m_nsp);
76 m_sqvisc.resize(m_nsp);
77 m_phi.resize(m_nsp, m_nsp, 0.0);
79 m_cond.resize(m_nsp);
80
81 // make a local copy of the molecular weights
83
86 for (size_t j = 0; j < m_nsp; j++) {
87 for (size_t k = j; k < m_nsp; k++) {
88 m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
89 m_wratjk(k,j) = sqrt(m_wratjk(j,k));
90 m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
91 }
92 }
93}
94
96{
97 update_T();
98 update_C();
99
100 if (m_visc_ok) {
101 return m_viscmix;
102 }
103
104 double vismix = 0.0;
105 // update m_visc and m_phi if necessary
106 if (!m_viscwt_ok) {
108 }
109
110 multiply(m_phi, m_molefracs.data(), m_spwork.data());
111
112 for (size_t k : m_kNeutral) {
113 vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
114 }
115 m_viscmix = vismix;
116 return vismix;
117}
118
120{
121 update_T();
122 update_C();
123 if (!m_spcond_ok) {
124 updateCond_T();
125 }
126 if (!m_condmix_ok) {
127 double sum1 = 0.0, sum2 = 0.0;
128 for (size_t k : m_kNeutral) {
129 sum1 += m_molefracs[k] * m_cond[k];
130 sum2 += m_molefracs[k] / m_cond[k];
131 }
132 m_lambda = 0.5*(sum1 + 1.0/sum2);
133 m_condmix_ok = true;
134 }
135 return m_lambda;
136}
137
139{
140 vector<double> mobi(m_nsp);
141 getMobilities(&mobi[0]);
142 double p = m_thermo->pressure();
143 double sum = 0.0;
144 for (size_t k : m_kIon) {
145 double ND_k = m_molefracs[k] * p / m_kbt;
146 sum += ND_k * std::abs(m_speciesCharge[k]) * ElectronCharge * mobi[k];
147 }
148 if (m_kElectron != npos) {
149 sum += m_molefracs[m_kElectron] * p / m_kbt *
151 }
152 return sum;
153}
154
156{
158
159 // number of points to use in generating fit data
160 const size_t np = 50;
161 int degree = 4;
162 double dt = (m_thermo->maxTemp() - m_thermo->minTemp())/(np-1);
163 vector<double> tlog(np);
164 vector<double> w(np);
165
166 // generate array of log(t) values
167 for (size_t n = 0; n < np; n++) {
168 double t = m_thermo->minTemp() + dt*n;
169 tlog[n] = log(t);
170 }
171
172 // vector of polynomial coefficients
173 vector<double> c(degree + 1);
174 double err = 0.0, relerr = 0.0,
175 mxerr = 0.0, mxrelerr = 0.0;
176
177 vector<double> diff(np + 1);
178 // The array order still not ideal
179 for (size_t k = 0; k < m_nsp; k++) {
180 for (size_t j = k; j < m_nsp; j++) {
181 if (m_alpha[k] == 0.0 || m_alpha[j] == 0.0 ||
182 k == m_kElectron || j == m_kElectron) {
183 continue;
184 }
185 if (m_speciesCharge[k] == 0) {
186 if (m_speciesCharge[j] == 0) {
187 continue;
188 }
189 } else {
190 if (m_speciesCharge[j] != 0) {
191 continue;
192 }
193 }
194 for (size_t n = 0; n < np; n++) {
195 double t = m_thermo->minTemp() + dt*n;
196 double eps = m_epsilon(j,k);
197 double tstar = Boltzmann * t/eps;
198 double sigma = m_diam(j,k);
199 double om11 = omega11_n64(tstar, m_gamma(j,k))
200 * Pi * sigma * sigma;
201 // Stockmayer-(n,6,4) model is not suitable for collision
202 // between O2/O2- due to resonant charge transfer.
203 // Therefore, the experimental collision data is used instead.
204 if ((k == m_thermo->speciesIndex("O2-") ||
205 j == m_thermo->speciesIndex("O2-")) &&
206 (k == m_thermo->speciesIndex("O2") ||
207 j == m_thermo->speciesIndex("O2"))) {
208 om11 = poly5(t, m_om11_O2.data()) / 1e20;
209 }
210 double diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,j))
211 * pow(Boltzmann * t, 1.5) / om11;
212
213 diff[n] = diffcoeff/pow(t, 1.5);
214 w[n] = 1.0/(diff[n]*diff[n]);
215 }
216 polyfit(np, degree, tlog.data(), diff.data(), w.data(), c.data());
217
218 for (size_t n = 0; n < np; n++) {
219 double val, fit;
220 double t = exp(tlog[n]);
221 double pre = pow(t, 1.5);
222 val = pre * diff[n];
223 fit = pre * poly4(tlog[n], c.data());
224 err = fit - val;
225 relerr = err/val;
226 mxerr = std::max(mxerr, fabs(err));
227 mxrelerr = std::max(mxrelerr, fabs(relerr));
228 }
229 size_t sum = k * (k + 1) / 2;
230 m_diffcoeffs[k*m_nsp+j-sum] = c;
231 }
232 }
233 m_fittingErrors["diff-coeff-max-abs-error"] =
234 std::max(m_fittingErrors.getDouble("diff-coeff-max-abs-error", 0.0),
235 mxerr);
236 m_fittingErrors["diff-coeff-max-rel-error"] =
237 std::max(m_fittingErrors.getDouble("diff-coeff-max-rel-error", 0.0),
238 mxrelerr);
239}
240
242{
243 m_gamma.resize(m_nsp, m_nsp, 0.0);
244 for (size_t i : m_kIon) {
245 for (size_t j : m_kNeutral) {
246 if (m_alpha[j] != 0.0 && m_alpha[i] != 0.0) {
247 double r_alpha = m_alpha[i] / m_alpha[j];
248 // save a copy of polarizability in Angstrom
249 double alphaA_i = m_alpha[i] * 1e30;
250 double alphaA_j = m_alpha[j] * 1e30;
251 // The ratio of dispersion to induction forces
252 double xi = alphaA_i / (m_speciesCharge[i] * m_speciesCharge[i] *
253 (1.0 + pow(2 * r_alpha, 2./3.)) * sqrt(alphaA_j));
254
255 // the collision diameter
256 double K1 = 1.767;
257 double kappa = 0.095;
258 m_diam(i,j) = K1 * (pow(m_alpha[i], 1./3.) + pow(m_alpha[j], 1./3.)) /
259 pow(alphaA_i * alphaA_j * (1.0 + 1.0 / xi), kappa);
260
261 // The original K2 is 0.72, but Han et al. suggested that K2 = 1.44
262 // for better fit.
263 double K2 = 1.44;
264 double epsilon = K2 * ElectronCharge * ElectronCharge *
266 m_alpha[j] * (1.0 + xi) /
267 (8 * Pi * epsilon_0 * pow(m_diam(i,j),4));
268 if (epsilon != 0.0) {
269 m_epsilon(i,j) = epsilon;
270 }
271
272 // Calculate dispersion coefficient and quadrupole polarizability
273 // from curve fitting if not available.
274 // Neutrals
275 if (m_disp[j] == 0.0) {
276 m_disp[j] = exp(1.8846*log(alphaA_j)-0.4737)* 1e-50;
277 }
278 if (m_quad_polar[j] == 0.0) {
279 m_quad_polar[j] = 2.0 * m_disp[j];
280 }
281 // Ions
282 if (m_disp[i] == 0.0) {
283 if (m_speciesCharge[i] > 0) {
284 m_disp[i] = exp(1.8853*log(alphaA_i)+0.2682)* 1e-50;
285 } else {
286 m_disp[i] = exp(3.2246*log(alphaA_i)-3.2397)* 1e-50;
287 }
288 }
289
290 // The binary dispersion coefficient is determined by the combination rule
291 // Reference:
292 // Tang, K. T. "Dynamic polarizabilities and van der Waals coefficients."
293 // Physical Review 177.1 (1969): 108.
294 double C6 = 2.0 * m_disp[i] * m_disp[j] /
295 (1.0/r_alpha * m_disp[i] + r_alpha * m_disp[j]);
296
297 m_gamma(i,j) = (2.0 / pow(m_speciesCharge[i],2) * C6 + m_quad_polar[j]) /
298 (m_alpha[j] * m_diam(i,j) * m_diam(i,j));//Dimensionless
299
300 // properties are symmetric
301 m_diam(j,i) = m_diam(i,j);
302 m_epsilon(j,i) = m_epsilon(i,j);
303 m_gamma(j,i) = m_gamma(i,j);
304 }
305 }
306 }
307}
308
309double IonGasTransport::omega11_n64(const double tstar, const double gamma)
310{
311 double logtstar = log(tstar);
312 double om11 = 0.0;
313 if (tstar < 0.01) {
314 throw CanteraError("IonGasTransport::omega11_n64",
315 "tstar = {} is smaller than 0.01", tstar);
316 } else if (tstar <= 0.04) {
317 // for interval 0.01 to 0.04, SSE = 0.006; R^2 = 1; RMSE = 0.020
318 om11 = 2.97 - 12.0 * gamma
319 - 0.887 * logtstar
320 + 3.86 * gamma * gamma
321 - 6.45 * gamma * logtstar
322 - 0.275 * logtstar * logtstar
323 + 1.20 * gamma * gamma * logtstar
324 - 1.24 * gamma * logtstar * logtstar
325 - 0.164 * pow(logtstar,3);
326 } else if (tstar <= 1000) {
327 // for interval 0.04 to 1000, SSE = 0.282; R^2 = 1; RMSE = 0.033
328 om11 = 1.22 - 0.0343 * gamma
329 + (-0.769 + 0.232 * gamma) * logtstar
330 + (0.306 - 0.165 * gamma) * logtstar * logtstar
331 + (-0.0465 + 0.0388 * gamma) * pow(logtstar,3)
332 + (0.000614 - 0.00285 * gamma) * pow(logtstar,4)
333 + 0.000238 * pow(logtstar,5);
334 } else {
335 throw CanteraError("IonGasTransport::omega11_n64",
336 "tstar = {} is larger than 1000", tstar);
337 }
338 return om11;
339}
340
342{
343 update_T();
344 update_C();
345
346 // update the binary diffusion coefficients if necessary
347 if (!m_bindiff_ok) {
348 updateDiff_T();
349 }
350
351 double mmw = m_thermo->meanMolecularWeight();
352 double p = m_thermo->pressure();
353 if (m_nsp == 1) {
354 d[0] = m_bdiff(0,0) / p;
355 } else {
356 for (size_t k = 0; k < m_nsp; k++) {
357 if (k == m_kElectron) {
358 d[k] = 0.4 * m_kbt / ElectronCharge;
359 } else {
360 double sum2 = 0.0;
361 for (size_t j : m_kNeutral) {
362 if (j != k) {
363 sum2 += m_molefracs[j] / m_bdiff(j,k);
364 }
365 }
366 if (sum2 <= 0.0) {
367 d[k] = m_bdiff(k,k) / p;
368 } else {
369 d[k] = (mmw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2);
370 }
371 }
372 }
373 }
374}
375
376void IonGasTransport::getMobilities(double* const mobi)
377{
378 update_T();
379 update_C();
380
381 // update the binary diffusion coefficients if necessary
382 if (!m_bindiff_ok) {
383 updateDiff_T();
384 }
385 double p = m_thermo->pressure();
386 for (size_t k = 0; k < m_nsp; k++) {
387 if (k == m_kElectron) {
388 mobi[k] = 0.4;
389 } else {
390 mobi[k] = 0.0;
391 }
392 }
393 for (size_t k : m_kIon) {
394 double sum = 0.0;
395 for (size_t j : m_kNeutral) {
396 double bmobi = m_bdiff(k,j) * ElectronCharge / m_kbt;
397 sum += m_molefracs[j] / bmobi;
398 }
399 mobi[k] = 1.0 / sum / p;
400 }
401}
402
403}
Monchick and Mason collision integrals.
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
double getDouble(const string &key, double default_) const
If key exists, return it as a double, otherwise return default_.
Definition AnyMap.cpp:1580
Base class for exceptions thrown by Cantera classes.
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
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.
vector< double > m_molefracs
Vector of species mole fractions.
vector< double > m_quad_polar
Quadrupole polarizability.
bool m_visc_ok
Update boolean for mixture rule for the mixture viscosity.
DenseMatrix m_wratkj1
Holds square roots of molecular weight ratios.
virtual void fitDiffCoeffs(MMCollisionInt &integrals)
Generate polynomial fits to the binary diffusion coefficients.
vector< double > m_disp
Dispersion coefficient.
virtual void updateDiff_T()
Update the binary diffusion coefficients.
vector< double > m_sqvisc
vector of square root of species viscosities sqrt(kg /m /s).
DenseMatrix m_wratjk
Holds square roots of molecular weight ratios.
bool m_bindiff_ok
Update boolean for the binary diffusivities at unit pressure.
DenseMatrix m_epsilon
The effective well depth for (i,j) collisions.
DenseMatrix m_diam
hard-sphere diameter for (i,j) collision
vector< double > m_spwork
work space length = m_kk
int m_mode
Type of the polynomial fits to temperature.
virtual void updateViscosity_T()
Update the temperature-dependent viscosity terms.
double m_viscmix
Internal storage for the viscosity of the mixture (kg /m /s)
vector< double > m_visc
vector of species viscosities (kg /m /s).
DenseMatrix m_bdiff
Matrix of binary diffusion coefficients at the reference pressure and the current temperature Size is...
vector< vector< double > > m_diffcoeffs
Polynomial fits to the binary diffusivity of each species.
double m_kbt
Current value of Boltzmann constant times the temperature (Joules)
DenseMatrix m_reducedMass
This is the reduced mass of the interaction between species i and j.
bool m_viscwt_ok
Update boolean for the weighting factors for the mixture viscosity.
vector< double > m_alpha
Polarizability of each species in the phase.
DenseMatrix m_phi
m_phi is a Viscosity Weighting Function. size = m_nsp * n_nsp
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.
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).
void init(ThermoPhase *thermo, int mode) override
Initialize a transport manager.
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.
double m_lambda
Internal storage for the calculated mixture thermal conductivity.
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)
bool m_condmix_ok
Update boolean for the mixture rule for the mixture thermal conductivity.
size_t nSpecies() const
Returns the number of species in the phase.
Definition Phase.h:232
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
Definition Phase.h:656
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:581
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:539
Base class for a phase with thermodynamic properties.
virtual double minTemp(size_t k=npos) const
Minimum temperature for which the thermodynamic data for the species or phase are valid.
virtual double maxTemp(size_t k=npos) const
Maximum temperature for which the thermodynamic data for the species are valid.
ThermoPhase * m_thermo
pointer to the object representing the phase
Definition Transport.h:434
size_t m_nsp
Number of species.
Definition Transport.h:437
ThermoPhase & thermo()
Phase object.
Definition Transport.h:103
AnyMap m_fittingErrors
Maximum errors associated with fitting pure species transport properties.
Definition Transport.h:440
This file contains definitions for utility functions and text for modules, inputfiles and logging,...
R poly4(D x, R *c)
Evaluates a polynomial of order 4.
Definition utilities.h:153
R poly5(D x, R *c)
Templated evaluation of a polynomial of order 5.
Definition utilities.h:141
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:595
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...