Cantera  3.1.0a1
Loading...
Searching...
No Matches
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, int log_level)
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 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++) {
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
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 *
152 }
153 return sum;
154}
155
157{
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
315double 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
382void 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.
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.
int m_log_level
Level of verbose printing during initialization.
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.
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.
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: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.
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: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)
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 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: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...