Cantera  2.0
TransportFactory.cpp
Go to the documentation of this file.
1 /**
2  * @file TransportFactory.cpp
3  *
4  * Implementation file for class TransportFactory.
5  */
6 
7 #include "cantera/thermo/ThermoPhase.h"
8 
9 // known transport models
10 #include "cantera/transport/MultiTransport.h"
11 #include "cantera/transport/MixTransport.h"
12 #include "cantera/transport/SolidTransport.h"
13 #include "cantera/transport/DustyGasTransport.h"
14 #include "cantera/transport/SimpleTransport.h"
15 #include "cantera/transport/LiquidTransport.h"
16 #include "cantera/transport/AqueousTransport.h"
17 #include "cantera/transport/TransportFactory.h"
18 
19 #include "cantera/numerics/polyfit.h"
20 #include "MMCollisionInt.h"
21 #include "cantera/base/xml.h"
22 #include "cantera/base/XML_Writer.h"
23 #include "cantera/transport/TransportParams.h"
24 #include "cantera/transport/LiquidTransportParams.h"
25 #include "cantera/transport/LiquidTranInteraction.h"
26 #include "cantera/base/global.h"
27 #include "cantera/thermo/IdealGasPhase.h"
28 #include "cantera/base/ctml.h"
29 #include "cantera/base/stringUtils.h"
30 
31 #include <cstdio>
32 #include <cstring>
33 #include <fstream>
34 
35 using namespace std;
36 
37 
38 //! polynomial degree used for fitting collision integrals
39 //! except in CK mode, where the degree is 6.
40 #define COLL_INT_POLY_DEGREE 8
41 
42 namespace Cantera
43 {
44 /////////////////////////// constants //////////////////////////
45 //@ \cond
46 const doublereal ThreeSixteenths = 3.0/16.0;
47 const doublereal TwoOverPi = 2.0/Pi;
48 const doublereal FiveThirds = 5.0/3.0;
49 //@ \endcond
50 
51 
52 //====================================================================================================================
53 TransportFactory* TransportFactory::s_factory = 0;
54 
55 // declaration of static storage for the mutex
56 mutex_t TransportFactory::transport_mutex;
57 
58 ////////////////////////// exceptions /////////////////////////
59 
60 //====================================================================================================================
61 //! Exception thrown if an error is encountered while reading the transport database
62 class TransportDBError : public CanteraError
63 {
64 public:
65  //! Default constructor
66  /*!
67  * @param linenum inputs the line number
68  * @param msg String message to be sent to the user
69  */
70  TransportDBError(int linenum, std::string msg) :
71  CanteraError("getTransportData", "error reading transport data: " + msg + "\n") {
72  }
73 };
74 //====================================================================================================================
75 
76 //////////////////// class TransportFactory methods //////////////
77 
78 
79 void TransportFactory::getBinDiffCorrection(doublereal t,
80  const GasTransportParams& tr, MMCollisionInt& integrals,
81  size_t k, size_t j, doublereal xk, doublereal xj,
82  doublereal& fkj, doublereal& fjk)
83 {
84  doublereal w1, w2, wsum, sig1, sig2, sig12, sigratio, sigratio2,
85  sigratio3, tstar1, tstar2, tstar12,
86  om22_1, om22_2, om11_12, astar_12, bstar_12, cstar_12,
87  cnst, wmwp, sqw12, p1, p2, p12, q1, q2, q12;
88 
89  w1 = tr.mw[k];
90  w2 = tr.mw[j];
91  wsum = w1 + w2;
92  wmwp = (w1 - w2)/wsum;
93  sqw12 = sqrt(w1*w2);
94 
95  sig1 = tr.sigma[k];
96  sig2 = tr.sigma[j];
97  sig12 = 0.5*(tr.sigma[k] + tr.sigma[j]);
98  sigratio = sig1*sig1/(sig2*sig2);
99  sigratio2 = sig1*sig1/(sig12*sig12);
100  sigratio3 = sig2*sig2/(sig12*sig12);
101 
102  tstar1 = Boltzmann * t / tr.eps[k];
103  tstar2 = Boltzmann * t / tr.eps[j];
104  tstar12 = Boltzmann * t / sqrt(tr.eps[k] * tr.eps[j]);
105 
106  om22_1 = integrals.omega22(tstar1, tr.delta(k,k));
107  om22_2 = integrals.omega22(tstar2, tr.delta(j,j));
108  om11_12 = integrals.omega11(tstar12, tr.delta(k,j));
109  astar_12 = integrals.astar(tstar12, tr.delta(k,j));
110  bstar_12 = integrals.bstar(tstar12, tr.delta(k,j));
111  cstar_12 = integrals.cstar(tstar12, tr.delta(k,j));
112 
113  cnst = sigratio * sqrt(2.0*w2/wsum) * 2.0 *
114  w1*w1/(wsum * w2);
115  p1 = cnst * om22_1 / om11_12;
116 
117  cnst = (1.0/sigratio) * sqrt(2.0*w1/wsum) * 2.0*w2*w2/(wsum*w1);
118  p2 = cnst * om22_2 / om11_12;
119 
120  p12 = 15.0 * wmwp*wmwp + 8.0*w1*w2*astar_12/(wsum*wsum);
121 
122  cnst = (2.0/(w2*wsum))*sqrt(2.0*w2/wsum)*sigratio2;
123  q1 = cnst*((2.5 - 1.2*bstar_12)*w1*w1 + 3.0*w2*w2
124  + 1.6*w1*w2*astar_12);
125 
126  cnst = (2.0/(w1*wsum))*sqrt(2.0*w1/wsum)*sigratio3;
127  q2 = cnst*((2.5 - 1.2*bstar_12)*w2*w2 + 3.0*w1*w1
128  + 1.6*w1*w2*astar_12);
129 
130  q12 = wmwp*wmwp*15.0*(2.5 - 1.2*bstar_12)
131  + 4.0*w1*w2*astar_12*(11.0 - 2.4*bstar_12)/(wsum*wsum)
132  + 1.6*wsum*om22_1*om22_2/(om11_12*om11_12*sqw12)
133  * sigratio2 * sigratio3;
134 
135  cnst = 6.0*cstar_12 - 5.0;
136  fkj = 1.0 + 0.1*cnst*cnst *
137  (p1*xk*xk + p2*xj*xj + p12*xk*xj)/
138  (q1*xk*xk + q2*xj*xj + q12*xk*xj);
139  fjk = 1.0 + 0.1*cnst*cnst *
140  (p2*xk*xk + p1*xj*xj + p12*xk*xj)/
141  (q2*xk*xk + q1*xj*xj + q12*xk*xj);
142 }
143 //=============================================================================================================================
144 // Corrections for polar-nonpolar binary diffusion coefficients
145 /*
146  * Calculate corrections to the well depth parameter and the
147  * diameter for use in computing the binary diffusion coefficient
148  * of polar-nonpolar pairs. For more information about this
149  * correction, see Dixon-Lewis, Proc. Royal Society (1968).
150  *
151  * @param i Species one - this is a bimolecular correction routine
152  * @param j species two - this is a bimolecular correction routine
153  * @param tr Database of species properties read in from the input xml file.
154  * @param f_eps Multiplicative correction factor to be applied to epsilon(i,j)
155  * @param f_sigma Multiplicative correction factor to be applied to diam(i,j)
156  */
157 void TransportFactory::makePolarCorrections(size_t i, size_t j,
158  const GasTransportParams& tr, doublereal& f_eps, doublereal& f_sigma)
159 {
160 
161  // no correction if both are nonpolar, or both are polar
162  if (tr.polar[i] == tr.polar[j]) {
163  f_eps = 1.0;
164  f_sigma = 1.0;
165  return;
166  }
167 
168  // corrections to the effective diameter and well depth
169  // if one is polar and one is non-polar
170 
171  size_t kp = (tr.polar[i] ? i : j); // the polar one
172  size_t knp = (i == kp ? j : i); // the nonpolar one
173 
174  doublereal d3np, d3p, alpha_star, mu_p_star, xi;
175  d3np = pow(tr.sigma[knp],3);
176  d3p = pow(tr.sigma[kp],3);
177  alpha_star = tr.alpha[knp]/d3np;
178  mu_p_star = tr.dipole(kp,kp)/sqrt(d3p * tr.eps[kp]);
179  xi = 1.0 + 0.25 * alpha_star * mu_p_star * mu_p_star *
180  sqrt(tr.eps[kp]/tr.eps[knp]);
181  f_sigma = pow(xi, -1.0/6.0);
182  f_eps = xi*xi;
183 }
184 //=============================================================================================================================
185 /*
186  TransportFactory(): default constructor
187 
188  The default constructor for this class sets up
189  m_models[], a mapping between the string name
190  for a transport model and the integer name.
191 */
192 TransportFactory::TransportFactory() :
193  m_verbose(false)
194 {
195  m_models["Mix"] = cMixtureAveraged;
196  m_models["Multi"] = cMulticomponent;
197  m_models["Solid"] = cSolidTransport;
198  m_models["DustyGas"] = cDustyGasTransport;
199  m_models["CK_Multi"] = CK_Multicomponent;
200  m_models["CK_Mix"] = CK_MixtureAveraged;
201  m_models["Liquid"] = cLiquidTransport;
202  m_models["Aqueous"] = cAqueousTransport;
203  m_models["Simple"] = cSimpleTransport;
204  m_models["User"] = cUserTransport;
205  m_models["None"] = None;
206  //m_models["Radiative"] = cRadiative;
207 
208  m_tranPropMap["viscosity"] = TP_VISCOSITY;
209  m_tranPropMap["ionConductivity"] = TP_IONCONDUCTIVITY;
210  m_tranPropMap["mobilityRatio"] = TP_MOBILITYRATIO;
211  m_tranPropMap["selfDiffusion"] = TP_SELFDIFFUSION;
212  m_tranPropMap["thermalConductivity"] = TP_THERMALCOND;
213  m_tranPropMap["speciesDiffusivity"] = TP_DIFFUSIVITY;
214  m_tranPropMap["hydrodynamicRadius"] = TP_HYDRORADIUS;
215  m_tranPropMap["electricalConductivity"] = TP_ELECTCOND;
216 
217  m_LTRmodelMap[""] = LTP_TD_CONSTANT;
218  m_LTRmodelMap["constant"] = LTP_TD_CONSTANT;
219  m_LTRmodelMap["arrhenius"] = LTP_TD_ARRHENIUS;
220  m_LTRmodelMap["coeffs"] = LTP_TD_POLY;
221  m_LTRmodelMap["exptemp"] = LTP_TD_EXPT;
222 
223  m_LTImodelMap[""] = LTI_MODEL_NOTSET;
224  m_LTImodelMap["none"] = LTI_MODEL_NONE;
225  m_LTImodelMap["solvent"] = LTI_MODEL_SOLVENT;
226  m_LTImodelMap["moleFractions"] = LTI_MODEL_MOLEFRACS;
227  m_LTImodelMap["massFractions"] = LTI_MODEL_MASSFRACS;
228  m_LTImodelMap["logMoleFractions"] = LTI_MODEL_LOG_MOLEFRACS;
229  m_LTImodelMap["pairwiseInteraction"] = LTI_MODEL_PAIRWISE_INTERACTION;
230  m_LTImodelMap["stefanMaxwell_PPN"] = LTI_MODEL_STEFANMAXWELL_PPN;
231  m_LTImodelMap["moleFractionsExpT"] = LTI_MODEL_MOLEFRACS_EXPT;
232 }
233 
234 
235 // This static function deletes the statically allocated instance.
236 void TransportFactory::deleteFactory()
237 {
238  ScopedLock transportLock(transport_mutex);
239  if (s_factory) {
240  delete s_factory;
241  s_factory = 0;
242  }
243 }
244 
245 /*
246  make one of several transport models, and return a base class
247  pointer to it. This method operates at the level of a
248  single transport property as a function of temperature
249  and possibly composition.
250 */
251 LTPspecies* TransportFactory::newLTP(const XML_Node& trNode, std::string& name,
252  TransportPropertyType tp_ind, thermo_t* thermo)
253 {
254  LTPspecies* ltps = 0;
255  std::string model = lowercase(trNode["model"]);
256  switch (m_LTRmodelMap[model]) {
257  case LTP_TD_CONSTANT:
258  ltps = new LTPspecies_Const(trNode, name, tp_ind, thermo);
259  break;
260  case LTP_TD_ARRHENIUS:
261  ltps = new LTPspecies_Arrhenius(trNode, name, tp_ind, thermo);
262  break;
263  case LTP_TD_POLY:
264  ltps = new LTPspecies_Poly(trNode, name, tp_ind, thermo);
265  break;
266  case LTP_TD_EXPT:
267  ltps = new LTPspecies_ExpT(trNode, name, tp_ind, thermo);
268  break;
269  default:
270  throw CanteraError("newLTP","unknown transport model: " + model);
271  ltps = new LTPspecies(&trNode, name, tp_ind, thermo);
272  }
273  return ltps;
274 }
275 
276 /*
277  make one of several transport models, and return a base class
278  pointer to it. This method operates at the level of a
279  single mixture transport property. Individual species
280  transport properties are addressed by the LTPspecies
281  returned by newLTP
282 */
283 LiquidTranInteraction* TransportFactory::newLTI(const XML_Node& trNode,
284  TransportPropertyType tp_ind,
285  LiquidTransportParams& trParam)
286 {
287  LiquidTranInteraction* lti = 0;
288 
289  thermo_t* thermo = trParam.thermo;
290 
291  std::string model = trNode["model"];
292  switch (m_LTImodelMap[model]) {
293  case LTI_MODEL_SOLVENT:
294  lti = new LTI_Solvent(tp_ind);
295  lti->init(trNode, thermo);
296  break;
297  case LTI_MODEL_MOLEFRACS:
298  lti = new LTI_MoleFracs(tp_ind);
299  lti->init(trNode, thermo);
300  break;
301  case LTI_MODEL_MASSFRACS:
302  lti = new LTI_MassFracs(tp_ind);
303  lti->init(trNode, thermo);
304  break;
305  case LTI_MODEL_LOG_MOLEFRACS:
306  lti = new LTI_Log_MoleFracs(tp_ind);
307  lti->init(trNode, thermo);
308  break;
309  case LTI_MODEL_PAIRWISE_INTERACTION:
310  lti = new LTI_Pairwise_Interaction(tp_ind);
311  lti->init(trNode, thermo);
312  lti->setParameters(trParam);
313  break;
314  case LTI_MODEL_STEFANMAXWELL_PPN:
315  lti = new LTI_StefanMaxwell_PPN(tp_ind);
316  lti->init(trNode, thermo);
317  lti->setParameters(trParam);
318  break;
319  case LTI_MODEL_STOKES_EINSTEIN:
320  lti = new LTI_StokesEinstein(tp_ind);
321  lti->init(trNode, thermo);
322  lti->setParameters(trParam);
323  break;
324  case LTI_MODEL_MOLEFRACS_EXPT:
325  lti = new LTI_MoleFracs_ExpT(tp_ind);
326  lti->init(trNode, thermo);
327  break;
328  default:
329  // throw CanteraError("newLTI","unknown transport model: " + model );
330  lti = new LiquidTranInteraction(tp_ind);
331  lti->init(trNode, thermo);
332  }
333  return lti;
334 }
335 
336 /*
337  make one of several transport models, and return a base class
338  pointer to it.
339 */
340 Transport* TransportFactory::newTransport(std::string transportModel,
341  thermo_t* phase, int log_level)
342 {
343 
344  if (transportModel == "") {
345  return new Transport;
346  }
347 
348 
349  vector_fp state;
350  Transport* tr = 0, *gastr = 0;
351  DustyGasTransport* dtr = 0;
352  phase->saveState(state);
353 
354  switch (m_models[transportModel]) {
355  case None:
356  tr = new Transport;
357  break;
358  case cMulticomponent:
359  tr = new MultiTransport;
360  initTransport(tr, phase, 0, log_level);
361  break;
362  case CK_Multicomponent:
363  tr = new MultiTransport;
364  initTransport(tr, phase, CK_Mode, log_level);
365  break;
366  case cMixtureAveraged:
367  tr = new MixTransport;
368  initTransport(tr, phase, 0, log_level);
369  break;
370  case CK_MixtureAveraged:
371  tr = new MixTransport;
372  initTransport(tr, phase, CK_Mode, log_level);
373  break;
374  case cSolidTransport:
375  tr = new SolidTransport;
376  tr->setThermo(*phase);
377  break;
378  case cDustyGasTransport:
379  tr = new DustyGasTransport;
380  gastr = new MultiTransport;
381  initTransport(gastr, phase, 0, log_level);
382  dtr = (DustyGasTransport*)tr;
383  dtr->initialize(phase, gastr);
384  break;
385  case cSimpleTransport:
386  tr = new SimpleTransport();
387  initLiquidTransport(tr, phase, log_level);
388  tr->setThermo(*phase);
389  break;
390  case cLiquidTransport:
391  tr = new LiquidTransport;
392  initLiquidTransport(tr, phase, log_level);
393  tr->setThermo(*phase);
394  break;
395  case cAqueousTransport:
396  tr = new AqueousTransport;
397  initLiquidTransport(tr, phase, log_level);
398  tr->setThermo(*phase);
399  break;
400  default:
401  throw CanteraError("newTransport","unknown transport model: " + transportModel);
402  }
403  phase->restoreState(state);
404  return tr;
405 }
406 
407 /*
408  make one of several transport models, and return a base class
409  pointer to it.
410 */
411 Transport* TransportFactory::newTransport(thermo_t* phase, int log_level)
412 {
413  XML_Node& phaseNode=phase->xml();
414  /*
415  * Find the Thermo XML node
416  */
417  if (!phaseNode.hasChild("transport")) {
418  throw CanteraError("TransportFactory::newTransport",
419  "no transport XML node");
420  }
421  XML_Node& transportNode = phaseNode.child("transport");
422  std::string transportModel = transportNode.attrib("model");
423  if (transportModel == "") {
424  throw CanteraError("TransportFactory::newTransport",
425  "transport XML node doesn't have a model string");
426  }
427  return newTransport(transportModel, phase,log_level);
428 }
429 
430 //====================================================================================================================
431 // Prepare to build a new kinetic-theory-based transport manager for low-density gases
432 /*
433  * This class fills up the GastransportParams structure for the current phase
434  *
435  * Uses polynomial fits to Monchick & Mason collision integrals. store then in tr
436  *
437  * @param flog Reference to the ostream for writing log info
438  * @param transport_database Reference to a vector of pointers containing the
439  * transport database for each species
440  * @param thermo Pointer to the %ThermoPhase object
441  * @param mode Mode -> Either it's CK_Mode, chemkin compatibility mode, or it is not
442  * We usually run with chemkin compatibility mode turned off.
443  * @param log_level log level
444  * @param tr GasTransportParams structure to be filled up with information
445  */
446 void TransportFactory::setupMM(std::ostream& flog, const std::vector<const XML_Node*> &transport_database,
447  thermo_t* thermo, int mode, int log_level, GasTransportParams& tr)
448 {
449 
450  // constant mixture attributes
451  tr.thermo = thermo;
452  tr.nsp_ = tr.thermo->nSpecies();
453  size_t nsp = tr.nsp_;
454 
455  tr.tmin = thermo->minTemp();
456  tr.tmax = thermo->maxTemp();
457  tr.mw.resize(nsp);
458  tr.log_level = log_level;
459 
460  copy(tr.thermo->molecularWeights().begin(), tr.thermo->molecularWeights().end(), tr.mw.begin());
461 
462  tr.mode_ = mode;
463  tr.epsilon.resize(nsp, nsp, 0.0);
464  tr.delta.resize(nsp, nsp, 0.0);
465  tr.reducedMass.resize(nsp, nsp, 0.0);
466  tr.dipole.resize(nsp, nsp, 0.0);
467  tr.diam.resize(nsp, nsp, 0.0);
468  tr.crot.resize(nsp);
469  tr.zrot.resize(nsp);
470  tr.polar.resize(nsp, false);
471  tr.alpha.resize(nsp, 0.0);
472  tr.poly.resize(nsp);
473  tr.sigma.resize(nsp);
474  tr.eps.resize(nsp);
475 
476  XML_Node root, log;
477  getTransportData(transport_database, log, tr.thermo->speciesNames(), tr);
478 
479  for (size_t i = 0; i < nsp; i++) {
480  tr.poly[i].resize(nsp);
481  }
482 
483  doublereal ts1, ts2, tstar_min = 1.e8, tstar_max = 0.0;
484  doublereal f_eps, f_sigma;
485 
486  DenseMatrix& diam = tr.diam;
487  DenseMatrix& epsilon = tr.epsilon;
488 
489  for (size_t i = 0; i < nsp; i++) {
490  for (size_t j = i; j < nsp; j++) {
491  // the reduced mass
492  tr.reducedMass(i,j) = tr.mw[i] * tr.mw[j] / (Avogadro * (tr.mw[i] + tr.mw[j]));
493 
494  // hard-sphere diameter for (i,j) collisions
495  diam(i,j) = 0.5*(tr.sigma[i] + tr.sigma[j]);
496 
497  // the effective well depth for (i,j) collisions
498  epsilon(i,j) = sqrt(tr.eps[i]*tr.eps[j]);
499 
500  // The polynomial fits of collision integrals vs. T*
501  // will be done for the T* from tstar_min to tstar_max
502  ts1 = Boltzmann * tr.tmin/epsilon(i,j);
503  ts2 = Boltzmann * tr.tmax/epsilon(i,j);
504  if (ts1 < tstar_min) {
505  tstar_min = ts1;
506  }
507  if (ts2 > tstar_max) {
508  tstar_max = ts2;
509  }
510 
511  // the effective dipole moment for (i,j) collisions
512  tr.dipole(i,j) = sqrt(tr.dipole(i,i)*tr.dipole(j,j));
513 
514  // reduced dipole moment delta* (nondimensional)
515  doublereal d = diam(i,j);
516  tr.delta(i,j) = 0.5 * tr.dipole(i,j)*tr.dipole(i,j)
517  / (epsilon(i,j) * d * d * d);
518 
519  makePolarCorrections(i, j, tr, f_eps, f_sigma);
520  tr.diam(i,j) *= f_sigma;
521  epsilon(i,j) *= f_eps;
522 
523  // properties are symmetric
524  tr.reducedMass(j,i) = tr.reducedMass(i,j);
525  diam(j,i) = diam(i,j);
526  epsilon(j,i) = epsilon(i,j);
527  tr.dipole(j,i) = tr.dipole(i,j);
528  tr.delta(j,i) = tr.delta(i,j);
529  }
530  }
531 
532  // Chemkin fits the entire T* range in the Monchick and Mason tables,
533  // so modify tstar_min and tstar_max if in Chemkin compatibility mode
534 
535  if (mode == CK_Mode) {
536  tstar_min = 0.101;
537  tstar_max = 99.9;
538  }
539 
540 
541  // initialize the collision integral calculator for the desired
542  // T* range
543 #ifdef DEBUG_MODE
544  if (m_verbose) {
545  tr.xml->XML_open(flog, "collision_integrals");
546  }
547 #endif
548  MMCollisionInt integrals;
549  integrals.init(tr.xml, tstar_min, tstar_max, log_level);
550  fitCollisionIntegrals(flog, tr, integrals);
551 #ifdef DEBUG_MODE
552  if (m_verbose) {
553  tr.xml->XML_close(flog, "collision_integrals");
554  }
555 #endif
556  // make polynomial fits
557 #ifdef DEBUG_MODE
558  if (m_verbose) {
559  tr.xml->XML_open(flog, "property fits");
560  }
561 #endif
562  fitProperties(tr, integrals, flog);
563 #ifdef DEBUG_MODE
564  if (m_verbose) {
565  tr.xml->XML_close(flog, "property fits");
566  }
567 #endif
568 }
569 
570 //====================================================================================================================
571 // Prepare to build a new transport manager for liquids assuming that
572 // viscosity transport data is provided in Arrhenius form.
573 /*
574  * @param flog Reference to the ostream for writing log info
575  * @param thermo Pointer to the %ThermoPhase object
576  * @param log_level log level
577  * @param trParam LiquidTransportParams structure to be filled up with information
578  */
579 void TransportFactory::setupLiquidTransport(std::ostream& flog, thermo_t* thermo, int log_level,
580  LiquidTransportParams& trParam)
581 {
582 
583  const std::vector<const XML_Node*> & species_database = thermo->speciesData();
584  const XML_Node* phase_database = &thermo->xml();
585 
586  // constant mixture attributes
587  trParam.thermo = thermo;
588  trParam.nsp_ = trParam.thermo->nSpecies();
589  size_t nsp = trParam.nsp_;
590 
591  trParam.tmin = thermo->minTemp();
592  trParam.tmax = thermo->maxTemp();
593  trParam.log_level = log_level;
594 
595  // Get the molecular weights and load them into trParam
596  trParam.mw.resize(nsp);
597  copy(trParam.thermo->molecularWeights().begin(),
598  trParam.thermo->molecularWeights().end(), trParam.mw.begin());
599 
600  // Resize all other vectors in trParam
601  trParam.LTData.resize(nsp);
602 
603  // Need to identify a method to obtain interaction matrices.
604  // This will fill LiquidTransportParams members visc_Eij, visc_Sij
605  // trParam.visc_Eij.resize(nsp,nsp);
606  // trParam.visc_Sij.resize(nsp,nsp);
607  trParam.thermalCond_Aij.resize(nsp,nsp);
608  trParam.diff_Dij.resize(nsp,nsp);
609  trParam.radius_Aij.resize(nsp,nsp);
610 
611  XML_Node root, log;
612  // Note that getLiquidSpeciesTransportData just populates the pure species transport data.
613  getLiquidSpeciesTransportData(species_database, log, trParam.thermo->speciesNames(), trParam);
614 
615  // getLiquidInteractionsTransportData() populates the
616  // species-species interaction models parameters
617  // like visc_Eij
618  if (phase_database->hasChild("transport")) {
619  XML_Node& transportNode = phase_database->child("transport");
620  getLiquidInteractionsTransportData(transportNode, log, trParam.thermo->speciesNames(), trParam);
621  }
622 }
623 //====================================================================================================================
624 // Initialize an existing transport manager
625 /*
626  * This routine sets up an existing gas-phase transport manager.
627  * It calculates the collision integrals and calls the initGas() function to
628  * populate the species-dependent data structure.
629  *
630  * @param tr Pointer to the Transport manager
631  * @param thermo Pointer to the ThermoPhase object
632  * @param mode Chemkin compatible mode or not. This alters the specification of the
633  * collision integrals. defaults to no.
634  * @param log_level Defaults to zero, no logging
635  *
636  * In DEBUG_MODE, this routine will create the file transport_log.xml
637  * and write informative information to it.
638  */
639 void TransportFactory::initTransport(Transport* tran,
640  thermo_t* thermo, int mode, int log_level)
641 {
642  ScopedLock transportLock(transport_mutex);
643  const std::vector<const XML_Node*> & transport_database = thermo->speciesData();
644 
645  GasTransportParams trParam;
646 #ifdef DEBUG_MODE
647  if (log_level == 0) {
648  m_verbose = 0;
649  }
650  ofstream flog("transport_log.xml");
651  trParam.xml = new XML_Writer(flog);
652  if (m_verbose) {
653  trParam.xml->XML_open(flog, "transport");
654  }
655 #else
656  // create the object, but don't associate it with a file
657  std::ostream& flog(std::cout);
658 #endif
659  // set up Monchick and Mason collision integrals
660  setupMM(flog, transport_database, thermo, mode, log_level, trParam);
661  // do model-specific initialization
662  tran->initGas(trParam);
663 #ifdef DEBUG_MODE
664  if (m_verbose) {
665  trParam.xml->XML_close(flog, "transport");
666  }
667  // finished with log file
668  flog.close();
669 #endif
670  return;
671 }
672 //====================================================================================================================
673 
674 /* Similar to initTransport except uses LiquidTransportParams
675  class and calls setupLiquidTransport().
676 */
677 void TransportFactory::initLiquidTransport(Transport* tran,
678  thermo_t* thermo,
679  int log_level)
680 {
681 
682  LiquidTransportParams trParam;
683 #ifdef DEBUG_MODE
684  ofstream flog("transport_log.xml");
685  trParam.xml = new XML_Writer(flog);
686  if (m_verbose) {
687  trParam.xml->XML_open(flog, "transport");
688  }
689 #else
690  // create the object, but don't associate it with a file
691  std::ostream& flog(std::cout);
692 #endif
693  setupLiquidTransport(flog, thermo, log_level, trParam);
694  // do model-specific initialization
695  tran->initLiquid(trParam);
696 #ifdef DEBUG_MODE
697  if (m_verbose) {
698  trParam.xml->XML_close(flog, "transport");
699  }
700  // finished with log file
701  flog.close();
702 #endif
703  return;
704 
705 }
706 
707 void TransportFactory::fitCollisionIntegrals(ostream& logfile,
708  GasTransportParams& tr,
709  MMCollisionInt& integrals)
710 {
711  vector_fp::iterator dptr;
712  doublereal dstar;
713  size_t nsp = tr.nsp_;
714  int mode = tr.mode_;
715  size_t i, j;
716 
717  // Chemkin fits to sixth order polynomials
718  int degree = (mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE);
719 #ifdef DEBUG_MODE
720  if (m_verbose) {
721  tr.xml->XML_open(logfile, "tstar_fits");
722  tr.xml->XML_comment(logfile, "fits to A*, B*, and C* vs. log(T*).\n"
723  "These are done only for the required dstar(j,k) values.");
724  if (tr.log_level < 3) {
725  tr.xml->XML_comment(logfile, "*** polynomial coefficients not printed (log_level < 3) ***");
726  }
727  }
728 #endif
729  for (i = 0; i < nsp; i++) {
730  for (j = i; j < nsp; j++) {
731  // Chemkin fits only delta* = 0
732  if (mode != CK_Mode) {
733  dstar = tr.delta(i,j);
734  } else {
735  dstar = 0.0;
736  }
737 
738  // if a fit has already been generated for
739  // delta* = tr.delta(i,j), then use it. Otherwise,
740  // make a new fit, and add tr.delta(i,j) to the list
741  // of delta* values for which fits have been done.
742 
743  // 'find' returns a pointer to end() if not found
744  dptr = find(tr.fitlist.begin(), tr.fitlist.end(), dstar);
745  if (dptr == tr.fitlist.end()) {
746  vector_fp ca(degree+1), cb(degree+1), cc(degree+1);
747  vector_fp co22(degree+1);
748  integrals.fit(logfile, degree, dstar,
749  DATA_PTR(ca), DATA_PTR(cb), DATA_PTR(cc));
750  integrals.fit_omega22(logfile, degree, dstar,
751  DATA_PTR(co22));
752  tr.omega22_poly.push_back(co22);
753  tr.astar_poly.push_back(ca);
754  tr.bstar_poly.push_back(cb);
755  tr.cstar_poly.push_back(cc);
756  tr.poly[i][j] = static_cast<int>(tr.astar_poly.size()) - 1;
757  tr.fitlist.push_back(dstar);
758  }
759 
760  // delta* found in fitlist, so just point to this
761  // polynomial
762  else {
763  tr.poly[i][j] = static_cast<int>((dptr - tr.fitlist.begin()));
764  }
765  tr.poly[j][i] = tr.poly[i][j];
766  }
767  }
768 #ifdef DEBUG_MODE
769  if (m_verbose) {
770  tr.xml->XML_close(logfile, "tstar_fits");
771  }
772 #endif
773 }
774 //====================================================================================================================
775 
776 
777 
778 /*********************************************************
779  *
780  * Read Transport Database
781  *
782  *********************************************************/
783 
784 /*
785  Read transport property data from a file for a list of species.
786  Given the name of a file containing transport property
787  parameters and a list of species names, this method returns an
788  instance of TransportParams containing the transport data for
789  these species read from the file.
790 */
791 void TransportFactory::getTransportData(const std::vector<const XML_Node*> &xspecies,
792  XML_Node& log, const std::vector<std::string> &names, GasTransportParams& tr)
793 {
794  std::string name;
795  int geom;
796  std::map<std::string, GasTransportData> datatable;
797  doublereal welldepth, diam, dipole, polar, rot;
798 
799  size_t nsp = xspecies.size();
800 
801  // read all entries in database into 'datatable' and check for
802  // errors. Note that this procedure validates all entries, not
803  // only those for the species listed in 'names'.
804 
805  std::string val, type;
806  map<std::string, int> gindx;
807  gindx["atom"] = 100;
808  gindx["linear"] = 101;
809  gindx["nonlinear"] = 102;
810  int linenum = 0;
811  for (size_t i = 0; i < nsp; i++) {
812  const XML_Node& sp = *xspecies[i];
813  name = sp["name"];
814  // std::cout << "Processing node for " << name << std::endl;
815 
816  // put in a try block so that species with no 'transport'
817  // child are skipped, instead of throwing an exception.
818  try {
819  XML_Node& tr = sp.child("transport");
820  ctml::getString(tr, "geometry", val, type);
821  geom = gindx[val] - 100;
822  map<std::string, doublereal> fv;
823 
824  welldepth = ctml::getFloat(tr, "LJ_welldepth");
825  diam = ctml::getFloat(tr, "LJ_diameter");
826  dipole = ctml::getFloat(tr, "dipoleMoment");
827  polar = ctml::getFloat(tr, "polarizability");
828  rot = ctml::getFloat(tr, "rotRelax");
829 
830  GasTransportData data;
831  data.speciesName = name;
832  data.geometry = geom;
833  if (welldepth >= 0.0) {
834  data.wellDepth = welldepth;
835  } else throw TransportDBError(linenum,
836  "negative well depth");
837 
838  if (diam > 0.0) {
839  data.diameter = diam;
840  } else throw TransportDBError(linenum,
841  "negative or zero diameter");
842 
843  if (dipole >= 0.0) {
844  data.dipoleMoment = dipole;
845  } else throw TransportDBError(linenum,
846  "negative dipole moment");
847 
848  if (polar >= 0.0) {
849  data.polarizability = polar;
850  } else throw TransportDBError(linenum,
851  "negative polarizability");
852 
853  if (rot >= 0.0) {
854  data.rotRelaxNumber = rot;
855  } else throw TransportDBError(linenum,
856  "negative rotation relaxation number");
857 
858  datatable[name] = data;
859  } catch (CanteraError& err) {
860  err.save();
861  }
862  }
863 
864  for (size_t i = 0; i < tr.nsp_; i++) {
865 
866  GasTransportData& trdat = datatable[names[i]];
867 
868  // 'datatable' returns a default TransportData object if
869  // the species name is not one in the transport database.
870  // This can be detected by examining 'geometry'.
871  if (trdat.geometry < 0) {
872  throw TransportDBError(0,"no transport data found for species "
873  + names[i]);
874  }
875 
876  // parameters are converted to SI units before storing
877 
878  // rotational heat capacity / R
879  switch (trdat.geometry) {
880  case 0:
881  tr.crot[i] = 0.0; // monatomic
882  break;
883  case 1:
884  tr.crot[i] = 1.0; // linear
885  break;
886  default:
887  tr.crot[i] = 1.5; // nonlinear
888  }
889 
890 
891  tr.dipole(i,i) = 1.e-25 * SqrtTen * trdat.dipoleMoment;
892 
893  if (trdat.dipoleMoment > 0.0) {
894  tr.polar[i] = true;
895  } else {
896  tr.polar[i] = false;
897  }
898 
899  // A^3 -> m^3
900  tr.alpha[i] = 1.e-30 * trdat.polarizability;
901 
902  tr.sigma[i] = 1.e-10 * trdat.diameter;
903 
904  tr.eps[i] = Boltzmann * trdat.wellDepth;
905  tr.zrot[i] = std::max(1.0, trdat.rotRelaxNumber);
906 
907  }
908 }
909 
910 /*
911  Read transport property data from a file for a list of species.
912  Given the name of a file containing transport property
913  parameters and a list of species names, this method returns an
914  instance of TransportParams containing the transport data for
915  these species read from the file.
916 */
917 void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML_Node*> &xspecies,
918  XML_Node& log,
919  const std::vector<std::string> &names,
920  LiquidTransportParams& trParam)
921 {
922  std::string name;
923  /*
924  Create a map of species names versus liquid transport data parameters
925  */
926  std::map<std::string, LiquidTransportData> datatable;
927  std::map<std::string, LiquidTransportData>::iterator it;
928 
929  // Store the number of species in the phase
930  size_t nsp = trParam.nsp_;
931 
932  // Store the number of off-diagonal symmetric interactions between species in the phase
933  size_t nBinInt = nsp*(nsp-1)/2;
934 
935  // read all entries in database into 'datatable' and check for
936  // errors. Note that this procedure validates all entries, not
937  // only those for the species listed in 'names'.
938  for (size_t i = 0; i < nsp; i++) {
939  const XML_Node& sp = *xspecies[i];
940  name = sp["name"];
941  vector_fp vCoeff;
942 
943  // Species with no 'transport' child are skipped. However, if that species is in the list,
944  // it will throw an exception below.
945  try {
946  if (sp.hasChild("transport")) {
947  XML_Node& trNode = sp.child("transport");
948 
949  // Fill datatable with LiquidTransportData objects for error checking
950  // and then insertion into LiquidTransportData objects below.
951  LiquidTransportData data;
952  data.speciesName = name;
953  data.mobilityRatio.resize(nsp*nsp,0);
954  data.selfDiffusion.resize(nsp,0);
955  ThermoPhase* temp_thermo = trParam.thermo;
956 
957  size_t num = trNode.nChildren();
958  for (size_t iChild = 0; iChild < num; iChild++) {
959  XML_Node& xmlChild = trNode.child(iChild);
960  std::string nodeName = xmlChild.name();
961 
962  switch (m_tranPropMap[nodeName]) {
963  case TP_VISCOSITY:
964  data.viscosity = newLTP(xmlChild, name, m_tranPropMap[nodeName], temp_thermo);
965  break;
966  case TP_IONCONDUCTIVITY:
967  data.ionConductivity = newLTP(xmlChild, name, m_tranPropMap[nodeName], temp_thermo);
968  break;
969  case TP_MOBILITYRATIO: {
970  for (size_t iSpec = 0; iSpec< nBinInt; iSpec++) {
971  XML_Node& propSpecNode = xmlChild.child(iSpec);
972  std::string specName = propSpecNode.name();
973  size_t loc = specName.find(":");
974  std::string firstSpec = specName.substr(0,loc);
975  std::string secondSpec = specName.substr(loc+1);
976  size_t index = temp_thermo->speciesIndex(firstSpec.c_str())+nsp*temp_thermo->speciesIndex(secondSpec.c_str());
977  data.mobilityRatio[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], temp_thermo);
978  };
979  };
980  break;
981  case TP_SELFDIFFUSION: {
982  for (size_t iSpec = 0; iSpec< nsp; iSpec++) {
983  XML_Node& propSpecNode = xmlChild.child(iSpec);
984  std::string specName = propSpecNode.name();
985  size_t index = temp_thermo->speciesIndex(specName.c_str());
986  data.selfDiffusion[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], temp_thermo);
987  };
988  };
989  break;
990  case TP_THERMALCOND:
991  data.thermalCond = newLTP(xmlChild,
992  name,
993  m_tranPropMap[nodeName],
994  temp_thermo);
995  break;
996  case TP_DIFFUSIVITY:
997  data.speciesDiffusivity = newLTP(xmlChild,
998  name,
999  m_tranPropMap[nodeName],
1000  temp_thermo);
1001  break;
1002  case TP_HYDRORADIUS:
1003  data.hydroRadius = newLTP(xmlChild,
1004  name,
1005  m_tranPropMap[nodeName],
1006  temp_thermo);
1007  break;
1008  case TP_ELECTCOND:
1009  data.electCond = newLTP(xmlChild,
1010  name,
1011  m_tranPropMap[nodeName],
1012  temp_thermo);
1013 
1014  break;
1015  default:
1016  throw CanteraError("getLiquidSpeciesTransportData","unknown transport property: " + nodeName);
1017  }
1018 
1019  }
1020  datatable.insert(pair<std::string, LiquidTransportData>(name,data));
1021  }
1022  } catch (CanteraError& err) {
1023  err.save();
1024  throw err;
1025  }
1026  }
1027 
1028  trParam.LTData.clear();
1029  for (size_t i = 0; i < trParam.nsp_; i++) {
1030  /*
1031  Check to see that we have a LiquidTransportData object for all of the
1032  species in the phase. If not, throw an error.
1033  */
1034  it = datatable.find(names[i]);
1035  if (it == datatable.end()) {
1036  throw TransportDBError(0,"No transport data found for species " + names[i]);
1037  }
1038  LiquidTransportData& trdat = it->second;
1039 
1040  /*
1041  Now, transfer these objects into LTData in the correct phase index order by
1042  calling the default copy constructor for LiquidTransportData.
1043  */
1044  trParam.LTData.push_back(trdat);
1045  }
1046 }
1047 
1048 
1049 /*
1050  Read transport property data from a file for interactions
1051  between species in a liquid.
1052  Given the name of a file containing transport property
1053  parameters and a list of species names, this method returns an
1054  instance of TransportParams containing the transport data for
1055  these species read from the file.
1056 */
1057 void TransportFactory::getLiquidInteractionsTransportData(const XML_Node& transportNode,
1058  XML_Node& log,
1059  const std::vector<std::string> &names,
1060  LiquidTransportParams& trParam)
1061 {
1062  try {
1063 
1064  size_t nsp = trParam.nsp_;
1065  size_t nBinInt = nsp*(nsp-1)/2;
1066 
1067  size_t num = transportNode.nChildren();
1068  for (size_t iChild = 0; iChild < num; iChild++) {
1069  //tranTypeNode is a type of transport property like viscosity
1070  XML_Node& tranTypeNode = transportNode.child(iChild);
1071  std::string nodeName = tranTypeNode.name();
1072 
1073  trParam.mobilityRatio.resize(nsp*nsp,0);
1074  trParam.selfDiffusion.resize(nsp,0);
1075  ThermoPhase* temp_thermo = trParam.thermo;
1076 
1077  if (tranTypeNode.hasChild("compositionDependence")) {
1078  //compDepNode contains the interaction model
1079  XML_Node& compDepNode = tranTypeNode.child("compositionDependence");
1080  switch (m_tranPropMap[nodeName]) {
1081  break;
1082  case TP_VISCOSITY:
1083  trParam.viscosity = newLTI(compDepNode, m_tranPropMap[nodeName], trParam);
1084  break;
1085  case TP_IONCONDUCTIVITY:
1086  trParam.ionConductivity = newLTI(compDepNode,
1087  m_tranPropMap[nodeName],
1088  trParam);
1089  break;
1090  case TP_MOBILITYRATIO: {
1091  for (size_t iSpec = 0; iSpec< nBinInt; iSpec++) {
1092  XML_Node& propSpecNode = compDepNode.child(iSpec);
1093  string specName = propSpecNode.name();
1094  size_t loc = specName.find(":");
1095  string firstSpec = specName.substr(0,loc);
1096  string secondSpec = specName.substr(loc+1);
1097  size_t index = temp_thermo->speciesIndex(firstSpec.c_str())+nsp*temp_thermo->speciesIndex(secondSpec.c_str());
1098  trParam.mobilityRatio[index] = newLTI(propSpecNode,
1099  m_tranPropMap[nodeName],
1100  trParam);
1101  };
1102  };
1103  break;
1104  case TP_SELFDIFFUSION: {
1105  for (size_t iSpec = 0; iSpec< nsp; iSpec++) {
1106  XML_Node& propSpecNode = compDepNode.child(iSpec);
1107  string specName = propSpecNode.name();
1108  size_t index = temp_thermo->speciesIndex(specName.c_str());
1109  trParam.selfDiffusion[index] = newLTI(propSpecNode,
1110  m_tranPropMap[nodeName],
1111  trParam);
1112  };
1113  };
1114  break;
1115  case TP_THERMALCOND:
1116  trParam.thermalCond = newLTI(compDepNode,
1117  m_tranPropMap[nodeName],
1118  trParam);
1119  break;
1120  case TP_DIFFUSIVITY:
1121  trParam.speciesDiffusivity = newLTI(compDepNode,
1122  m_tranPropMap[nodeName],
1123  trParam);
1124  break;
1125  case TP_HYDRORADIUS:
1126  trParam.hydroRadius = newLTI(compDepNode,
1127  m_tranPropMap[nodeName],
1128  trParam);
1129  break;
1130  case TP_ELECTCOND:
1131  trParam.electCond = newLTI(compDepNode,
1132  m_tranPropMap[nodeName],
1133  trParam);
1134  break;
1135  default:
1136  throw CanteraError("getLiquidInteractionsTransportData","unknown transport property: " + nodeName);
1137 
1138  }
1139  }
1140  /* Allow a switch between mass-averaged, mole-averaged
1141  * and solvent specified reference velocities.
1142  * XML code within the transportProperty node
1143  * (i.e. within <viscosity>) should read as follows
1144  * <velocityBasis basis="mass"> <!-- mass averaged -->
1145  * <velocityBasis basis="mole"> <!-- mole averaged -->
1146  * <velocityBasis basis="H2O"> <!-- H2O solvent -->
1147  */
1148  if (tranTypeNode.hasChild("velocityBasis")) {
1149  std::string velocityBasis =
1150  tranTypeNode.child("velocityBasis").attrib("basis");
1151  if (velocityBasis == "mass") {
1152  trParam.velocityBasis_ = VB_MASSAVG;
1153  } else if (velocityBasis == "mole") {
1154  trParam.velocityBasis_ = VB_MOLEAVG;
1155  } else if (trParam.thermo->speciesIndex(velocityBasis) > 0) {
1156  trParam.velocityBasis_ = static_cast<int>(trParam.thermo->speciesIndex(velocityBasis));
1157  } else {
1158  int linenum = __LINE__;
1159  throw TransportDBError(linenum, "Unknown attribute \"" + velocityBasis + "\" for <velocityBasis> node. ");
1160  }
1161  }
1162  }
1163  } catch (CanteraError& err) {
1164  std::cout << err.what() << std::endl;
1165  }
1166  return;
1167 }
1168 
1169 /*********************************************************
1170  *
1171  * Polynomial fitting
1172  *
1173  *********************************************************/
1174 
1175 void TransportFactory::fitProperties(GasTransportParams& tr,
1176  MMCollisionInt& integrals, std::ostream& logfile)
1177 {
1178  doublereal tstar;
1179  int ndeg = 0;
1180 #ifdef DEBUG_MODE
1181  char s[100];
1182 #endif
1183  // number of points to use in generating fit data
1184  const size_t np = 50;
1185 
1186  int mode = tr.mode_;
1187  int degree = (mode == CK_Mode ? 3 : 4);
1188 
1189  doublereal t, om22;
1190  doublereal dt = (tr.tmax - tr.tmin)/(np-1);
1191  vector_fp tlog(np), spvisc(np), spcond(np);
1192  doublereal val, fit;
1193 
1194  vector_fp w(np), w2(np);
1195 
1196  // generate array of log(t) values
1197  for (size_t n = 0; n < np; n++) {
1198  t = tr.tmin + dt*n;
1199  tlog[n] = log(t);
1200  }
1201 
1202  // vector of polynomial coefficients
1203  vector_fp c(degree + 1), c2(degree + 1);
1204 
1205 
1206  // fit the pure-species viscosity and thermal conductivity for
1207  // each species
1208 #ifdef DEBUG_MODE
1209  if (tr.log_level < 2 && m_verbose) {
1210  tr.xml->XML_comment(logfile,
1211  "*** polynomial coefficients not printed (log_level < 2) ***");
1212  }
1213 #endif
1214  doublereal sqrt_T, visc, err, relerr,
1215  mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0;
1216 
1217 #ifdef DEBUG_MODE
1218  if (m_verbose) {
1219  tr.xml->XML_open(logfile, "viscosity");
1220  tr.xml->XML_comment(logfile,"Polynomial fits for viscosity");
1221  if (mode == CK_Mode) {
1222  tr.xml->XML_comment(logfile,"log(viscosity) fit to cubic "
1223  "polynomial in log(T)");
1224  } else {
1225  sprintf(s, "viscosity/sqrt(T) fit to "
1226  "polynomial of degree %d in log(T)",degree);
1227  tr.xml->XML_comment(logfile,s);
1228  }
1229  }
1230 #endif
1231 
1232  doublereal cp_R, cond, w_RT, f_int, A_factor, B_factor,
1233  c1, cv_rot, cv_int, f_rot, f_trans, om11;
1234  doublereal diffcoeff;
1235 
1236  for (size_t k = 0; k < tr.nsp_; k++) {
1237  for (size_t n = 0; n < np; n++) {
1238  t = tr.tmin + dt*n;
1239 
1240  tr.thermo->setTemperature(t);
1241  cp_R = ((IdealGasPhase*)tr.thermo)->cp_R_ref()[k];
1242 
1243  tstar = Boltzmann * t/ tr.eps[k];
1244  sqrt_T = sqrt(t);
1245  om22 = integrals.omega22(tstar, tr.delta(k,k));
1246  om11 = integrals.omega11(tstar, tr.delta(k,k));
1247 
1248  // self-diffusion coefficient, without polar
1249  // corrections
1250  diffcoeff = ThreeSixteenths *
1251  sqrt(2.0 * Pi/tr.reducedMass(k,k)) *
1252  pow((Boltzmann * t), 1.5)/
1253  (Pi * tr.sigma[k] * tr.sigma[k] * om11);
1254 
1255  // viscosity
1256  visc = FiveSixteenths
1257  * sqrt(Pi * tr.mw[k] * Boltzmann * t / Avogadro) /
1258  (om22 * Pi * tr.sigma[k]*tr.sigma[k]);
1259 
1260  // thermal conductivity
1261  w_RT = tr.mw[k]/(GasConstant * t);
1262  f_int = w_RT * diffcoeff/visc;
1263  cv_rot = tr.crot[k];
1264 
1265  A_factor = 2.5 - f_int;
1266  B_factor = tr.zrot[k] + TwoOverPi
1267  *(FiveThirds * cv_rot + f_int);
1268  c1 = TwoOverPi * A_factor/B_factor;
1269  cv_int = cp_R - 2.5 - cv_rot;
1270 
1271  f_rot = f_int * (1.0 + c1);
1272  f_trans = 2.5 * (1.0 - c1 * cv_rot/1.5);
1273 
1274  cond = (visc/tr.mw[k])*GasConstant*(f_trans * 1.5
1275  + f_rot * cv_rot + f_int * cv_int);
1276 
1277  if (mode == CK_Mode) {
1278  spvisc[n] = log(visc);
1279  spcond[n] = log(cond);
1280  w[n] = -1.0;
1281  w2[n] = -1.0;
1282  } else {
1283  // the viscosity should be proportional
1284  // approximately to sqrt(T); therefore,
1285  // visc/sqrt(T) should have only a weak
1286  // temperature dependence. And since the mixture
1287  // rule requires the square root of the
1288  // pure-species viscosity, fit the square root of
1289  // (visc/sqrt(T)) to avoid having to compute
1290  // square roots in the mixture rule.
1291  spvisc[n] = sqrt(visc/sqrt_T);
1292 
1293  // the pure-species conductivity scales
1294  // approximately with sqrt(T). Unlike the
1295  // viscosity, there is no reason here to fit the
1296  // square root, since a different mixture rule is
1297  // used.
1298  spcond[n] = cond/sqrt_T;
1299  w[n] = 1.0/(spvisc[n]*spvisc[n]);
1300  w2[n] = 1.0/(spcond[n]*spcond[n]);
1301  }
1302  }
1303  polyfit(np, DATA_PTR(tlog), DATA_PTR(spvisc),
1304  DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c));
1305  polyfit(np, DATA_PTR(tlog), DATA_PTR(spcond),
1306  DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c2));
1307 
1308  // evaluate max fit errors for viscosity
1309  for (size_t n = 0; n < np; n++) {
1310  if (mode == CK_Mode) {
1311  val = exp(spvisc[n]);
1312  fit = exp(poly3(tlog[n], DATA_PTR(c)));
1313  } else {
1314  sqrt_T = exp(0.5*tlog[n]);
1315  val = sqrt_T * pow(spvisc[n],2);
1316  fit = sqrt_T * pow(poly4(tlog[n], DATA_PTR(c)),2);
1317  }
1318  err = fit - val;
1319  relerr = err/val;
1320  if (fabs(err) > mxerr) {
1321  mxerr = fabs(err);
1322  }
1323  if (fabs(relerr) > mxrelerr) {
1324  mxrelerr = fabs(relerr);
1325  }
1326  }
1327 
1328  // evaluate max fit errors for conductivity
1329  for (size_t n = 0; n < np; n++) {
1330  if (mode == CK_Mode) {
1331  val = exp(spcond[n]);
1332  fit = exp(poly3(tlog[n], DATA_PTR(c2)));
1333  } else {
1334  sqrt_T = exp(0.5*tlog[n]);
1335  val = sqrt_T * spcond[n];
1336  fit = sqrt_T * poly4(tlog[n], DATA_PTR(c2));
1337  }
1338  err = fit - val;
1339  relerr = err/val;
1340  if (fabs(err) > mxerr_cond) {
1341  mxerr_cond = fabs(err);
1342  }
1343  if (fabs(relerr) > mxrelerr_cond) {
1344  mxrelerr_cond = fabs(relerr);
1345  }
1346  }
1347  tr.visccoeffs.push_back(c);
1348  tr.condcoeffs.push_back(c2);
1349 
1350 #ifdef DEBUG_MODE
1351  if (tr.log_level >= 2 && m_verbose) {
1352  tr.xml->XML_writeVector(logfile, " ", tr.thermo->speciesName(k),
1353  c.size(), DATA_PTR(c));
1354  }
1355 #endif
1356  }
1357 #ifdef DEBUG_MODE
1358  if (m_verbose) {
1359  sprintf(s, "Maximum viscosity absolute error: %12.6g", mxerr);
1360  tr.xml->XML_comment(logfile,s);
1361  sprintf(s, "Maximum viscosity relative error: %12.6g", mxrelerr);
1362  tr.xml->XML_comment(logfile,s);
1363  tr.xml->XML_close(logfile, "viscosity");
1364 
1365 
1366  tr.xml->XML_open(logfile, "conductivity");
1367  tr.xml->XML_comment(logfile,"Polynomial fits for conductivity");
1368  if (mode == CK_Mode)
1369  tr.xml->XML_comment(logfile,"log(conductivity) fit to cubic "
1370  "polynomial in log(T)");
1371  else {
1372  sprintf(s, "conductivity/sqrt(T) fit to "
1373  "polynomial of degree %d in log(T)",degree);
1374  tr.xml->XML_comment(logfile,s);
1375  }
1376  if (tr.log_level >= 2)
1377  for (size_t k = 0; k < tr.nsp_; k++) {
1378  tr.xml->XML_writeVector(logfile, " ", tr.thermo->speciesName(k),
1379  degree+1, DATA_PTR(tr.condcoeffs[k]));
1380  }
1381  sprintf(s, "Maximum conductivity absolute error: %12.6g", mxerr_cond);
1382  tr.xml->XML_comment(logfile,s);
1383  sprintf(s, "Maximum conductivity relative error: %12.6g", mxrelerr_cond);
1384  tr.xml->XML_comment(logfile,s);
1385  tr.xml->XML_close(logfile, "conductivity");
1386 
1387  // fit the binary diffusion coefficients for each species pair
1388 
1389  tr.xml->XML_open(logfile, "binary_diffusion_coefficients");
1390  tr.xml->XML_comment(logfile, "binary diffusion coefficients");
1391  if (mode == CK_Mode)
1392  tr.xml->XML_comment(logfile,"log(D) fit to cubic "
1393  "polynomial in log(T)");
1394  else {
1395  sprintf(s, "D/T**(3/2) fit to "
1396  "polynomial of degree %d in log(T)",degree);
1397  tr.xml->XML_comment(logfile,s);
1398  }
1399  }
1400 #endif
1401 
1402  mxerr = 0.0, mxrelerr = 0.0;
1403  vector_fp diff(np + 1);
1404  doublereal eps, sigma;
1405  for (size_t k = 0; k < tr.nsp_; k++) {
1406  for (size_t j = k; j < tr.nsp_; j++) {
1407  for (size_t n = 0; n < np; n++) {
1408 
1409  t = tr.tmin + dt*n;
1410 
1411  eps = tr.epsilon(j,k);
1412  tstar = Boltzmann * t/eps;
1413  sigma = tr.diam(j,k);
1414  om11 = integrals.omega11(tstar, tr.delta(j,k));
1415 
1416  diffcoeff = ThreeSixteenths *
1417  sqrt(2.0 * Pi/tr.reducedMass(k,j)) *
1418  pow((Boltzmann * t), 1.5)/
1419  (Pi * sigma * sigma * om11);
1420 
1421 
1422  // 2nd order correction
1423  // NOTE: THIS CORRECTION IS NOT APPLIED
1424  doublereal fkj, fjk;
1425  getBinDiffCorrection(t, tr, integrals, k, j, 1.0, 1.0, fkj, fjk);
1426  //diffcoeff *= fkj;
1427 
1428 
1429  if (mode == CK_Mode) {
1430  diff[n] = log(diffcoeff);
1431  w[n] = -1.0;
1432  } else {
1433  diff[n] = diffcoeff/pow(t, 1.5);
1434  w[n] = 1.0/(diff[n]*diff[n]);
1435  }
1436  }
1437  polyfit(np, DATA_PTR(tlog), DATA_PTR(diff),
1438  DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c));
1439 
1440  doublereal pre;
1441  for (size_t n = 0; n < np; n++) {
1442  if (mode == CK_Mode) {
1443  val = exp(diff[n]);
1444  fit = exp(poly3(tlog[n], DATA_PTR(c)));
1445  } else {
1446  t = exp(tlog[n]);
1447  pre = pow(t, 1.5);
1448  val = pre * diff[n];
1449  fit = pre * poly4(tlog[n], DATA_PTR(c));
1450  }
1451  err = fit - val;
1452  relerr = err/val;
1453  if (fabs(err) > mxerr) {
1454  mxerr = fabs(err);
1455  }
1456  if (fabs(relerr) > mxrelerr) {
1457  mxrelerr = fabs(relerr);
1458  }
1459  }
1460  tr.diffcoeffs.push_back(c);
1461 #ifdef DEBUG_MODE
1462  if (tr.log_level >= 2 && m_verbose) {
1463  tr.xml->XML_writeVector(logfile, " ", tr.thermo->speciesName(k)
1464  + "__"+tr.thermo->speciesName(j), c.size(), DATA_PTR(c));
1465  }
1466 #endif
1467  }
1468  }
1469 #ifdef DEBUG_MODE
1470  if (m_verbose) {
1471  sprintf(s,"Maximum binary diffusion coefficient absolute error:"
1472  " %12.6g", mxerr);
1473  tr.xml->XML_comment(logfile,s);
1474  sprintf(s, "Maximum binary diffusion coefficient relative error:"
1475  "%12.6g", mxrelerr);
1476  tr.xml->XML_comment(logfile,s);
1477  tr.xml->XML_close(logfile, "binary_diffusion_coefficients");
1478  }
1479 #endif
1480 }
1481 //====================================================================================================================
1482 // Create a new transport manager instance.
1483 /*
1484  * @param transportModel String identifying the transport model to be instantiated, defaults to the empty string
1485  * @param thermo ThermoPhase object associated with the phase, defaults to null pointer
1486  * @param loglevel int containing the Loglevel, defaults to zero
1487  * @param f ptr to the TransportFactory object if it's been malloced.
1488  *
1489  * @ingroup transportProps
1490  */
1491 Transport* newTransportMgr(std::string transportModel, thermo_t* thermo, int loglevel, TransportFactory* f)
1492 {
1493  if (f == 0) {
1494  f = TransportFactory::factory();
1495  }
1496  Transport* ptr = f->newTransport(transportModel, thermo, loglevel);
1497  /*
1498  * Note: We delete the static s_factory instance here, instead of in
1499  * appdelete() in misc.cpp, to avoid linking problems involving
1500  * the need for multiple cantera and transport library statements
1501  * for applications that don't have transport in them.
1502  */
1503  return ptr;
1504 }
1505 //====================================================================================================================
1506 // Create a new transport manager instance.
1507 /*
1508  * @param thermo ThermoPhase object associated with the phase, defaults to null pointer
1509  * @param loglevel int containing the Loglevel, defaults to zero
1510  * @param f ptr to the TransportFactory object if it's been malloced.
1511  *
1512  * @ingroup transportProps
1513  */
1514 Transport* newDefaultTransportMgr(thermo_t* thermo, int loglevel, TransportFactory* f)
1515 {
1516  if (f == 0) {
1517  f = TransportFactory::factory();
1518  }
1519  Transport* ptr = f->newTransport(thermo, loglevel);
1520  /*
1521  * Note: We delete the static s_factory instance here, instead of in
1522  * appdelete() in misc.cpp, to avoid linking problems involving
1523  * the need for multiple cantera and transport library statements
1524  * for applications that don't have transport in them.
1525  */
1526  return ptr;
1527 }
1528 //====================================================================================================================
1529 }
1530 
Generated by   doxygen 1.8.2