Cantera  2.0
LiquidTranInteraction.cpp
1 /**
2  * @file LiquidTransportParams.cpp
3  * Source code for liquid mixture transport property evaluations.
4  */
5 
6 #include "cantera/transport/LiquidTransportParams.h"
7 #include <iostream>
8 #include "cantera/thermo/IonsFromNeutralVPSSTP.h"
9 #include "cantera/thermo/MargulesVPSSTP.h"
10 #include "cantera/base/stringUtils.h"
11 
12 #include <stdlib.h>
13 using namespace std;
14 using namespace ctml;
15 
16 namespace Cantera
17 {
18 
19 /**
20  * Exception thrown if an error is encountered while reading the
21  * transport database.
22  */
23 class LTPError : public CanteraError
24 {
25 public:
26  LTPError(std::string msg)
27  : CanteraError("LTPspecies",
28  "error parsing transport data: "
29  + msg + "\n") {}
30 };
31 
32 /**
33  * Exception thrown if an error is encountered while reading the
34  * transport database.
35  */
36 class LTPmodelError : public CanteraError
37 {
38 public:
39  LTPmodelError(std::string msg)
40  : CanteraError("LTPspecies",
41  "error parsing transport data: "
42  + msg + "\n") {}
43 };
44 
45 
46 // Constructor
47 /*
48  * @param tp_ind Index indicating transport property type (i.e. viscosity)
49  */
50 LiquidTranInteraction::LiquidTranInteraction(TransportPropertyType tp_ind) :
51  m_model(LTI_MODEL_NOTSET),
52  m_property(tp_ind)
53 {
54 }
55 
56 LiquidTranInteraction::~LiquidTranInteraction()
57 {
58  size_t kmax = m_Aij.size();
59  for (size_t k = 0; k < kmax; k++) {
60  if (m_Aij[k]) {
61  delete m_Aij[k];
62  }
63  }
64  kmax = m_Bij.size();
65  for (size_t k = 0; k < kmax; k++) {
66  if (m_Bij[k]) {
67  delete m_Bij[k];
68  }
69  }
70  kmax = m_Hij.size();
71  for (size_t k = 0; k < kmax; k++) {
72  if (m_Hij[k]) {
73  delete m_Hij[k];
74  }
75  }
76  kmax = m_Sij.size();
77  for (size_t k = 0; k < kmax; k++) {
78  if (m_Sij[k]) {
79  delete m_Sij[k];
80  }
81  }
82 }
83 
84 //====================================================================================================================
85 
86 void LiquidTranInteraction::init(const XML_Node& compModelNode,
87  thermo_t* thermo)
88 {
89 
90  m_thermo = thermo;
91 
92  size_t nsp = thermo->nSpecies();
93  m_Dij.resize(nsp, nsp, 0.0);
94  m_Eij.resize(nsp, nsp, 0.0);
95  /*
96  m_Aij.resize(nsp);
97  m_Bij.resize(nsp);
98  m_Hij.resize(nsp);
99  m_Sij.resize(nsp);
100  for (int k = 0; k < nsp; k++ ){
101  (*m_Aij[k]).resize(nsp, nsp, 0.0);
102  (*m_Bij[k]).resize(nsp, nsp, 0.0);
103  (*m_Hij[k]).resize(nsp, nsp, 0.0);
104  (*m_Sij[k]).resize(nsp, nsp, 0.0);
105  }
106  */
107 
108  std::string speciesA;
109  std::string speciesB;
110 
111  size_t num = compModelNode.nChildren();
112  for (size_t iChild = 0; iChild < num; iChild++) {
113  XML_Node& xmlChild = compModelNode.child(iChild);
114  std::string nodeName = lowercase(xmlChild.name());
115  if (nodeName != "interaction") {
116  throw CanteraError("TransportFactory::getLiquidInteractionsTransportData",
117  "expected <interaction> element and got <" + nodeName + ">");
118  }
119  speciesA = xmlChild.attrib("speciesA");
120  speciesB = xmlChild.attrib("speciesB");
121  size_t iSpecies = m_thermo->speciesIndex(speciesA);
122  if (iSpecies == npos) {
123  throw CanteraError("TransportFactory::getLiquidInteractionsTransportData",
124  "Unknown species " + speciesA);
125  }
126  size_t jSpecies = m_thermo->speciesIndex(speciesB);
127  if (jSpecies == npos) {
128  throw CanteraError("TransportFactory::getLiquidInteractionsTransportData",
129  "Unknown species " + speciesB);
130  }
131  /* if (xmlChild.hasChild("Aij" ) ) {
132  m_Aij(iSpecies,jSpecies) = getFloat(xmlChild, "Aij", "toSI" );
133  m_Aij(jSpecies,iSpecies) = m_Aij(iSpecies,jSpecies) ;
134  }*/
135 
136  if (xmlChild.hasChild("Eij")) {
137  m_Eij(iSpecies,jSpecies) = getFloat(xmlChild, "Eij", "actEnergy");
138  m_Eij(iSpecies,jSpecies) /= GasConstant;
139  m_Eij(jSpecies,iSpecies) = m_Eij(iSpecies,jSpecies) ;
140  }
141 
142  if (xmlChild.hasChild("Aij")) {
143  vector_fp poly;
144  // poly0 = getFloat(poly, xmlChild, "Aij", "toSI" );
145  getFloatArray(xmlChild, poly, true, "toSI", "Aij");
146  // if (!poly.size() ) poly.push_back(poly0);
147  while (m_Aij.size()<poly.size()) {
148  DenseMatrix* aTemp = new DenseMatrix();
149  aTemp->resize(nsp, nsp, 0.0);
150  m_Aij.push_back(aTemp);
151  }
152  for (int i = 0; i < (int)poly.size(); i++) {
153  (*m_Aij[i])(iSpecies,jSpecies) = poly[i];
154  //(*m_Aij[i])(jSpecies,iSpecies) = (*m_Aij[i])(iSpecies,jSpecies) ;
155  }
156  }
157 
158  if (xmlChild.hasChild("Bij")) {
159  vector_fp poly;
160  getFloatArray(xmlChild, poly, true, "toSI", "Bij");
161  //if (!poly.size() ) poly.push_back(poly0);
162  while (m_Bij.size() < poly.size()) {
163  DenseMatrix* bTemp = new DenseMatrix();
164  bTemp->resize(nsp, nsp, 0.0);
165  m_Bij.push_back(bTemp);
166  }
167  for (size_t i=0; i<poly.size(); i++) {
168  (*m_Bij[i])(iSpecies,jSpecies) = poly[i];
169  //(*m_Bij[i])(jSpecies,iSpecies) = (*m_Bij[i])(iSpecies,jSpecies) ;
170  }
171  }
172 
173  if (xmlChild.hasChild("Hij")) {
174  vector_fp poly;
175  // poly0 = getFloat(poly, xmlChild, "Hij", "actEnergy" );
176  getFloatArray(xmlChild, poly, true, "actEnergy", "Hij");
177  // if (!poly.size() ) poly.push_back(poly0);
178  while (m_Hij.size()<poly.size()) {
179  DenseMatrix* hTemp = new DenseMatrix();
180  hTemp->resize(nsp, nsp, 0.0);
181  m_Hij.push_back(hTemp);
182  }
183  for (size_t i=0; i<poly.size(); i++) {
184  (*m_Hij[i])(iSpecies,jSpecies) = poly[i];
185  (*m_Hij[i])(iSpecies,jSpecies) /= GasConstant;
186  //(*m_Hij[i])(jSpecies,iSpecies) = (*m_Hij[i])(iSpecies,jSpecies) ;
187  }
188  }
189 
190  if (xmlChild.hasChild("Sij")) {
191  vector_fp poly;
192  // poly0 = getFloat(poly, xmlChild, "Sij", "actEnergy" );
193  getFloatArray(xmlChild, poly, true, "actEnergy", "Sij");
194  // if (!poly.size() ) poly.push_back(poly0);
195  while (m_Sij.size()<poly.size()) {
196  DenseMatrix* sTemp = new DenseMatrix();
197  sTemp->resize(nsp, nsp, 0.0);
198  m_Sij.push_back(sTemp);
199  }
200  for (size_t i=0; i<poly.size(); i++) {
201  (*m_Sij[i])(iSpecies,jSpecies) = poly[i];
202  (*m_Sij[i])(iSpecies,jSpecies) /= GasConstant;
203  //(*m_Sij[i])(jSpecies,iSpecies) = (*m_Sij[i])(iSpecies,jSpecies) ;
204  }
205  }
206 
207  /*0 if (xmlChild.hasChild("Sij" ) ) {
208  m_Sij(iSpecies,jSpecies) = getFloat(xmlChild, "Sij", "toSI" );
209  m_Sij(iSpecies,jSpecies) /= GasConstant;
210  //m_Sij(jSpecies,iSpecies) = m_Sij(iSpecies,jSpecies) ;
211  }*/
212 
213  if (xmlChild.hasChild("Dij")) {
214  m_Dij(iSpecies,jSpecies) = getFloat(xmlChild, "Dij", "toSI");
215  m_Dij(jSpecies,iSpecies) = m_Dij(iSpecies,jSpecies) ;
216  }
217  }
218 }
219 
220 // Copy constructor
221 LiquidTranInteraction::LiquidTranInteraction(const LiquidTranInteraction& right)
222 {
223  *this = right; //use assignment operator to do other work
224 }
225 
226 // Assignment operator
227 LiquidTranInteraction& LiquidTranInteraction::operator=(const LiquidTranInteraction& right)
228 {
229  if (&right != this) {
230  m_model = right.m_model;
231  m_property = right.m_property;
232  m_thermo = right.m_thermo;
233  //m_trParam = right.m_trParam;
234  m_Aij = right.m_Aij;
235  m_Bij = right.m_Bij;
236  m_Eij = right.m_Eij;
237  m_Hij = right.m_Hij;
238  m_Sij = right.m_Sij;
239  m_Dij = right.m_Dij;
240  }
241  return *this;
242 }
243 
244 
245 //====================================================================================================================
246 LTI_Solvent::LTI_Solvent(TransportPropertyType tp_ind) :
247  LiquidTranInteraction(tp_ind)
248 {
249  m_model = LTI_MODEL_SOLVENT;
250 }
251 //====================================================================================================================
252 
253 doublereal LTI_Solvent::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
254 {
255 
256  size_t nsp = m_thermo->nSpecies();
257  doublereal temp = m_thermo->temperature();
258  vector_fp molefracs(nsp);
259  m_thermo->getMoleFractions(&molefracs[0]);
260 
261  doublereal value = 0.0;
262 
263  //if weightings are specified, use those
264  if (speciesWeight) {
265  for (size_t k = 0; k < nsp; k++) {
266  //molefracs[k] = molefracs[k];
267  // should be: molefracs[k] = molefracs[k]*speciesWeight[k]; for consistency, but weight(solvent)=1?
268  }
269  } else {
270  throw CanteraError("LTI_Solvent::getMixTransProp","You should be specifying the speciesWeight");
271  /* //This does not follow directly a solvent model
272  //although if the solvent mole fraction is dominant
273  //and the other species values are given or zero,
274  //it should work.
275  for (int k = 0; k < nsp; k++) {
276  value += speciesValues[k] * molefracs[k];
277  }*/
278  }
279 
280  for (size_t i = 0; i < nsp; i++) {
281  //presume that the weighting is set to 1.0 for solvent and 0.0 for everything else.
282  value += speciesValues[i] * speciesWeight[i];
283  if (i == 0) {
284  AssertTrace(speciesWeight[i] == 1.0);
285  } else {
286  AssertTrace(speciesWeight[i] == 0.0);
287  }
288  for (size_t j = 0; j < nsp; j++) {
289  for (size_t k = 0; k < m_Aij.size(); k++) {
290  value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i], (int) k);
291  }
292  for (size_t k = 0; k < m_Bij.size(); k++) {
293  value += molefracs[i]*molefracs[j]*(*m_Bij[k])(i,j)*temp*pow(molefracs[i], (int) k);
294  }
295  }
296  }
297 
298  return value;
299 }
300 //====================================================================================================================
301 doublereal LTI_Solvent::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
302 {
303 
304  size_t nsp = m_thermo->nSpecies();
305  doublereal temp = m_thermo->temperature();
306  vector_fp molefracs(nsp);
307  m_thermo->getMoleFractions(&molefracs[0]);
308 
309  doublereal value = 0.0;
310 
311  for (size_t k = 0; k < nsp; k++) {
312  //molefracs[k] = molefracs[k];
313  // should be: molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight(); for consistency, but weight(solvent)=1?
314  }
315 
316  for (size_t i = 0; i < nsp; i++) {
317  //presume that the weighting is set to 1.0 for solvent and 0.0 for everything else.
318  value += LTPptrs[i]->getSpeciesTransProp() * LTPptrs[i]->getMixWeight();
319  for (size_t j = 0; j < nsp; j++) {
320  for (size_t k = 0; k < m_Aij.size(); k++) {
321  value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i], (int) k);
322  }
323  for (size_t k = 0; k < m_Bij.size(); k++) {
324  value += molefracs[i]*molefracs[j]*(*m_Bij[k])(i,j)*temp*pow(molefracs[i], (int) k);
325  }
326  }
327  }
328 
329  return value;
330 }
331 //====================================================================================================================
332 void LTI_Solvent::getMatrixTransProp(DenseMatrix& mat, doublereal* speciesValues)
333 {
334  mat = (*m_Aij[0]);
335 }
336 //====================================================================================================================
337 
338 
339 doublereal LTI_MoleFracs::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
340 {
341 
342  size_t nsp = m_thermo->nSpecies();
343  doublereal temp = m_thermo->temperature();
344  vector_fp molefracs(nsp);
345  m_thermo->getMoleFractions(&molefracs[0]);
346 
347  doublereal value = 0;
348 
349  //if weightings are specified, use those
350  if (speciesWeight) {
351  for (size_t k = 0; k < nsp; k++) {
352  molefracs[k] = molefracs[k]*speciesWeight[k];
353  }
354  } else {
355  throw CanteraError("LTI_MoleFracs::getMixTransProp","You should be specifying the speciesWeight");
356  }
357 
358  for (size_t i = 0; i < nsp; i++) {
359  value += speciesValues[i] * molefracs[i];
360  for (size_t j = 0; j < nsp; j++) {
361  for (size_t k = 0; k < m_Aij.size(); k++) {
362  value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i], (int) k);
363  }
364  for (size_t k = 0; k < m_Bij.size(); k++) {
365  value += molefracs[i]*molefracs[j]*(*m_Bij[k])(i,j)*temp*pow(molefracs[i], (int) k);
366  }
367  }
368  }
369 
370  return value;
371 }
372 
373 
374 doublereal LTI_MoleFracs::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
375 {
376 
377  size_t nsp = m_thermo->nSpecies();
378  doublereal temp = m_thermo->temperature();
379  vector_fp molefracs(nsp);
380  m_thermo->getMoleFractions(&molefracs[0]);
381 
382  doublereal value = 0;
383 
384  for (size_t k = 0; k < nsp; k++) {
385  molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight();
386  }
387 
388  for (size_t i = 0; i < nsp; i++) {
389  value += LTPptrs[i]->getSpeciesTransProp() * molefracs[i];
390  for (size_t j = 0; j < nsp; j++) {
391  for (size_t k = 0; k < m_Aij.size(); k++) {
392  value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i], (int) k);
393  }
394  for (size_t k = 0; k < m_Bij.size(); k++) {
395  value += molefracs[i]*molefracs[j]*(*m_Bij[k])(i,j)*temp*pow(molefracs[i], (int) k);
396  }
397  }
398  }
399  return value;
400 }
401 
402 
403 doublereal LTI_MassFracs::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
404 {
405 
406  size_t nsp = m_thermo->nSpecies();
407  doublereal temp = m_thermo->temperature();
408  vector_fp massfracs(nsp);
409  m_thermo->getMassFractions(&massfracs[0]);
410 
411  doublereal value = 0;
412 
413  //if weightings are specified, use those
414  if (speciesWeight) {
415  for (size_t k = 0; k < nsp; k++) {
416  massfracs[k] = massfracs[k]*speciesWeight[k];
417  }
418  } else {
419  throw CanteraError("LTI_MassFracs::getMixTransProp","You should be specifying the speciesWeight");
420  }
421 
422  for (size_t i = 0; i < nsp; i++) {
423  value += speciesValues[i] * massfracs[i];
424  for (size_t j = 0; j < nsp; j++) {
425  for (size_t k = 0; k < m_Aij.size(); k++) {
426  value += massfracs[i]*massfracs[j]*(*m_Aij[k])(i,j)*pow(massfracs[i], (int) k);
427  }
428  for (size_t k = 0; k < m_Bij.size(); k++) {
429  value += massfracs[i]*massfracs[j]*(*m_Bij[k])(i,j)*temp*pow(massfracs[i], (int) k);
430  }
431  }
432  }
433 
434  return value;
435 }
436 
437 
438 doublereal LTI_MassFracs::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
439 {
440 
441  size_t nsp = m_thermo->nSpecies();
442  doublereal temp = m_thermo->temperature();
443  vector_fp massfracs(nsp);
444  m_thermo->getMassFractions(&massfracs[0]);
445 
446  doublereal value = 0;
447 
448  for (size_t k = 0; k < nsp; k++) {
449  massfracs[k] = massfracs[k]*LTPptrs[k]->getMixWeight();
450  }
451 
452  for (size_t i = 0; i < nsp; i++) {
453  value += LTPptrs[i]->getSpeciesTransProp() * massfracs[i];
454  for (size_t j = 0; j < nsp; j++) {
455  for (size_t k = 0; k < m_Aij.size(); k++) {
456  value += massfracs[i]*massfracs[j]*(*m_Aij[k])(i,j)*pow(massfracs[i], (int) k);
457  }
458  for (size_t k = 0; k < m_Bij.size(); k++) {
459  value += massfracs[i]*massfracs[j]*(*m_Bij[k])(i,j)*temp*pow(massfracs[i], (int) k);
460  }
461  }
462  }
463 
464  return value;
465 }
466 
467 
468 
469 
470 doublereal LTI_Log_MoleFracs::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
471 {
472 
473  size_t nsp = m_thermo->nSpecies();
474  doublereal temp = m_thermo->temperature();
475  vector_fp molefracs(nsp);
476  m_thermo->getMoleFractions(&molefracs[0]);
477 
478 
479 
480  doublereal value = 0;
481 
482  //if weightings are specified, use those
483  if (speciesWeight) {
484  for (size_t k = 0; k < nsp; k++) {
485  molefracs[k] = molefracs[k]*speciesWeight[k];
486  }
487  } else {
488  throw CanteraError("LTI_Log_MoleFracs::getMixTransProp","You probably should have a speciesWeight when you call getMixTransProp to convert ion mole fractions to molecular mole fractions");
489  }
490 
491  for (size_t i = 0; i < nsp; i++) {
492  value += log(speciesValues[i]) * molefracs[i];
493  for (size_t j = 0; j < nsp; j++) {
494  for (size_t k = 0; k < m_Hij.size(); k++) {
495  value += molefracs[i]*molefracs[j]*(*m_Hij[k])(i,j)/temp*pow(molefracs[i], (int) k);
496  //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl;
497  }
498  for (size_t k = 0; k < m_Sij.size(); k++) {
499  value -= molefracs[i]*molefracs[j]*(*m_Sij[k])(i,j)*pow(molefracs[i], (int) k);
500  //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl;
501  }
502  }
503  }
504 
505  value = exp(value);
506  return value;
507 }
508 
509 
510 doublereal LTI_Log_MoleFracs::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
511 {
512 
513  size_t nsp = m_thermo->nSpecies();
514  doublereal temp = m_thermo->temperature();
515  vector_fp molefracs(nsp);
516  m_thermo->getMoleFractions(&molefracs[0]);
517 
518 
519  doublereal value = 0;
520 
521  //if weightings are specified, use those
522 
523  for (size_t k = 0; k < nsp; k++) {
524  molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight();
525  }
526 
527  for (size_t i = 0; i < nsp; i++) {
528  value += log(LTPptrs[i]->getSpeciesTransProp()) * molefracs[i];
529  for (size_t j = 0; j < nsp; j++) {
530  for (size_t k = 0; k < m_Hij.size(); k++) {
531  value += molefracs[i]*molefracs[j]*(*m_Hij[k])(i,j)/temp*pow(molefracs[i], (int) k);
532  //cout << "1 = " << molefracs[i]+molefracs[j] << endl;
533  //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl;
534  }
535  for (size_t k = 0; k < m_Sij.size(); k++) {
536  value -= molefracs[i]*molefracs[j]*(*m_Sij[k])(i,j)*pow(molefracs[i], (int) k);
537  //cout << "1 = " << molefracs[i]+molefracs[j] << endl;
538  //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl;
539  }
540  }
541  }
542 
543  value = exp(value);
544  // cout << ", viscSpeciesA = " << LTPptrs[0]->getSpeciesTransProp() << endl;
545  //cout << ", viscSpeciesB = " << LTPptrs[1]->getSpeciesTransProp() << endl;
546  //cout << "value = " << value << " FINAL" << endl;
547  return value;
548 }
549 
550 
551 
552 
553 
554 void LTI_Pairwise_Interaction::setParameters(LiquidTransportParams& trParam)
555 {
556  size_t nsp = m_thermo->nSpecies();
557  m_diagonals.resize(nsp, 0);
558 
559  for (size_t k = 0; k < nsp; k++) {
560  Cantera::LiquidTransportData& ltd = trParam.LTData[k];
561  if (ltd.speciesDiffusivity) {
562  m_diagonals[k] = ltd.speciesDiffusivity;
563  }
564  }
565 }
566 
567 doublereal LTI_Pairwise_Interaction::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
568 {
569 
570  size_t nsp = m_thermo->nSpecies();
571  vector_fp molefracs(nsp);
572  m_thermo->getMoleFractions(&molefracs[0]);
573 
574  doublereal value = 0;
575 
576  throw LTPmodelError("Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense.");
577 
578  return value;
579 }
580 
581 
582 doublereal LTI_Pairwise_Interaction::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
583 {
584 
585  size_t nsp = m_thermo->nSpecies();
586  vector_fp molefracs(nsp);
587  m_thermo->getMoleFractions(&molefracs[0]);
588 
589  doublereal value = 0;
590 
591  throw LTPmodelError("Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense.");
592 
593  return value;
594 }
595 
596 void LTI_Pairwise_Interaction::getMatrixTransProp(DenseMatrix& mat, doublereal* speciesValues)
597 {
598 
599  size_t nsp = m_thermo->nSpecies();
600  doublereal temp = m_thermo->temperature();
601  vector_fp molefracs(nsp);
602  m_thermo->getMoleFractions(&molefracs[0]);
603 
604  mat.resize(nsp, nsp, 0.0);
605  for (size_t i = 0; i < nsp; i++)
606  for (size_t j = 0; j < i; j++) {
607  mat(i,j) = mat(j,i) = exp(m_Eij(i,j) / temp) / m_Dij(i,j);
608  }
609 
610  for (size_t i = 0; i < nsp; i++)
611  if (mat(i,i) == 0.0 && m_diagonals[i]) {
612  mat(i,i) = 1.0 / m_diagonals[i]->getSpeciesTransProp() ;
613  }
614 }
615 
616 
617 void LTI_StefanMaxwell_PPN::setParameters(LiquidTransportParams& trParam)
618 {
619  size_t nsp = m_thermo->nSpecies();
620  size_t nsp2 = nsp*nsp;
621  //vector<std
622 
623  m_ionCondMix = 0;
624  m_ionCondMixModel = trParam.ionConductivity;
625  //trParam.ionConductivity = 0;
626  m_ionCondSpecies.resize(nsp,0);
627  m_mobRatMix.resize(nsp,nsp,0.0);
628  m_mobRatMixModel.resize(nsp2);
629  m_mobRatSpecies.resize(nsp2);
630  m_selfDiffMix.resize(nsp,0.0);
631  m_selfDiffMixModel.resize(nsp);
632  m_selfDiffSpecies.resize(nsp);
633 
634  for (size_t k = 0; k < nsp2; k++) {
635  m_mobRatMixModel[k] = trParam.mobilityRatio[k];
636  //trParam.mobilityRatio[k] = 0;
637  m_mobRatSpecies[k].resize(nsp,0);
638  }
639  for (size_t k = 0; k < nsp; k++) {
640  m_selfDiffMixModel[k] = trParam.selfDiffusion[k];
641  //trParam.selfDiffusion[k] = 0;
642  m_selfDiffSpecies[k].resize(nsp,0);
643  }
644 
645  for (size_t k = 0; k < nsp; k++) {
646  Cantera::LiquidTransportData& ltd = trParam.LTData[k];
647  m_ionCondSpecies[k] = ltd.ionConductivity;
648  //ltd.ionConductivity = 0;
649  for (size_t j = 0; j < nsp2; j++) {
650  m_mobRatSpecies[j][k] = ltd.mobilityRatio[j];
651  //ltd.mobilityRatio[j] = 0;
652  }
653  for (size_t j = 0; j < nsp; j++) {
654  m_selfDiffSpecies[j][k] = ltd.selfDiffusion[j];
655  //ltd.selfDiffusion[j] = 0;
656  }
657  }
658 }
659 
660 doublereal LTI_StefanMaxwell_PPN::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
661 {
662 
663  size_t nsp = m_thermo->nSpecies();
664  vector_fp molefracs(nsp);
665  m_thermo->getMoleFractions(&molefracs[0]);
666 
667  doublereal value = 0;
668 
669  throw LTPmodelError("Calling LTI_StefanMaxwell_PPN::getMixTransProp does not make sense.");
670 
671  return value;
672 }
673 
674 
675 doublereal LTI_StefanMaxwell_PPN::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
676 {
677 
678  size_t nsp = m_thermo->nSpecies();
679  vector_fp molefracs(nsp);
680  m_thermo->getMoleFractions(&molefracs[0]);
681 
682  doublereal value = 0;
683 
684  throw LTPmodelError("Calling LTI_StefanMaxwell_PPN::getMixTransProp does not make sense.");
685 
686  return value;
687 }
688 
689 void LTI_StefanMaxwell_PPN::getMatrixTransProp(DenseMatrix& mat, doublereal* speciesValues)
690 {
691  //CAL
692 
693  IonsFromNeutralVPSSTP* ions_thermo = dynamic_cast<IonsFromNeutralVPSSTP*>(m_thermo);
694  size_t nsp = m_thermo->nSpecies();
695  if (nsp != 3) {
696  throw CanteraError("LTI_StefanMaxwell_PPN::getMatrixTransProp","Function may only be called with a 3-ion system");
697  }
698  doublereal temp = m_thermo->temperature();
699  vector_fp molefracs(nsp);
700  m_thermo->getMoleFractions(&molefracs[0]);
701  vector_fp neut_molefracs;
702  ions_thermo->getNeutralMolecMoleFractions(neut_molefracs);
703  vector<size_t> cation;
704  vector<size_t> anion;
705  ions_thermo->getCationList(cation);
706  ions_thermo->getAnionList(anion);
707 
708  // Reaction Coeffs and Charges
709  std::vector<double> viS(6);
710  std::vector<double> charges(3);
711  std::vector<size_t> neutMolIndex(3);
712  ions_thermo->getDissociationCoeffs(viS,charges,neutMolIndex);
713 
714  if (anion.size() != 1) {
715  throw CanteraError("LTI_StefanMaxwell_PPN::getMatrixTransProp","Must have one anion only for StefanMaxwell_PPN");
716  }
717  if (cation.size() != 2) {
718  throw CanteraError("LTI_StefanMaxwell_PPN::getMatrixTransProp","Must have two cations of equal charge for StefanMaxwell_PPN");
719  }
720  if (charges[cation[0]] != charges[cation[1]]) {
721  throw CanteraError("LTI_StefanMaxwell_PPN::getMatrixTransProp","Cations must be of equal charge for StefanMaxwell_PPN");
722  }
723 
724  m_ionCondMix = m_ionCondMixModel->getMixTransProp(m_ionCondSpecies);
725 
726  MargulesVPSSTP* marg_thermo = dynamic_cast<MargulesVPSSTP*>(ions_thermo->neutralMoleculePhase_);
727  doublereal vol = m_thermo->molarVolume();
728 
729  size_t k = 0;
730  for (size_t j = 0; j < nsp; j++) {
731  for (size_t i = 0; i < nsp; i++) {
732  if (m_mobRatMixModel[k]) {
733  m_mobRatMix(i,j) = m_mobRatMixModel[k]->getMixTransProp(m_mobRatSpecies[k]);
734  if (m_mobRatMix(i,j) > 0.0) {
735  m_mobRatMix(j,i) = 1.0/m_mobRatMix(i,j);
736  }
737  }
738  k++;
739  }
740  }
741 
742 
743  for (k = 0; k < nsp; k++) {
744  m_selfDiffMix[k] = m_selfDiffMixModel[k]->getMixTransProp(m_selfDiffSpecies[k]);
745  }
746 
747  //! @todo Suspicious implicit conversion from double to int.
748  int vP = max(viS[cation[0]],viS[cation[1]]);
749  int vM = viS[anion[0]];
750  int zP = charges[cation[0]];
751  int zM = charges[anion[0]];
752  doublereal xA, xB, eps;
753  doublereal inv_vP_vM_MutualDiff;
754  vector_fp dlnActCoeffdlnN_diag;
755  dlnActCoeffdlnN_diag.resize(neut_molefracs.size(),0.0);
756  marg_thermo->getdlnActCoeffdlnN_diag(&dlnActCoeffdlnN_diag[0]);
757 
758  xA = neut_molefracs[neutMolIndex[cation[0]]];
759  xB = neut_molefracs[neutMolIndex[cation[1]]];
760  eps = (1-m_mobRatMix(cation[1],cation[0]))/(xA+xB*m_mobRatMix(cation[1],cation[0]));
761  inv_vP_vM_MutualDiff = (xA*(1-xB+dlnActCoeffdlnN_diag[neutMolIndex[cation[1]]])/m_selfDiffMix[cation[1]]+xB*(1-xA+dlnActCoeffdlnN_diag[neutMolIndex[cation[0]]])/m_selfDiffMix[cation[0]]);
762 
763  mat.resize(nsp, nsp, 0.0);
764  mat(cation[0],cation[1]) = mat(cation[1],cation[0]) = (1+vM/vP)*(1+eps*xB)*(1-eps*xA)*inv_vP_vM_MutualDiff-zP*zP*Faraday*Faraday/GasConstant/temp/m_ionCondMix/vol;
765  mat(cation[0],anion[0]) = mat(anion[0],cation[0]) = (1+vP/vM)*(-eps*xB*(1-eps*xA)*inv_vP_vM_MutualDiff)-zP*zM*Faraday*Faraday/GasConstant/temp/m_ionCondMix/vol;
766  mat(cation[1],anion[0]) = mat(anion[0],cation[1]) = (1+vP/vM)*(eps*xA*(1+eps*xB)*inv_vP_vM_MutualDiff)-zP*zM*Faraday*Faraday/GasConstant/temp/m_ionCondMix/vol;
767 
768 }
769 
770 
771 doublereal LTI_StokesEinstein::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
772 {
773 
774  size_t nsp = m_thermo->nSpecies();
775  vector_fp molefracs(nsp);
776  m_thermo->getMoleFractions(&molefracs[0]);
777 
778  doublereal value = 0;
779 
780  throw LTPmodelError("Calling LTI_StokesEinstein::getMixTransProp does not make sense.");
781 
782  return value;
783 }
784 
785 
786 doublereal LTI_StokesEinstein::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
787 {
788 
789  size_t nsp = m_thermo->nSpecies();
790  vector_fp molefracs(nsp);
791  m_thermo->getMoleFractions(&molefracs[0]);
792 
793  doublereal value = 0;
794 
795  throw LTPmodelError("Calling LTI_StokesEinstein::getMixTransProp does not make sense.");
796 
797  return value;
798 }
799 
800 
801 
802 void LTI_StokesEinstein::setParameters(LiquidTransportParams& trParam)
803 {
804  size_t nsp = m_thermo->nSpecies();
805  m_viscosity.resize(nsp, 0);
806  m_hydroRadius.resize(nsp, 0);
807  for (size_t k = 0; k < nsp; k++) {
808  Cantera::LiquidTransportData& ltd = trParam.LTData[k];
809  m_viscosity[k] = ltd.viscosity;
810  m_hydroRadius[k] = ltd.hydroRadius;
811  }
812 }
813 
814 void LTI_StokesEinstein::getMatrixTransProp(DenseMatrix& mat, doublereal* speciesValues)
815 {
816  size_t nsp = m_thermo->nSpecies();
817  doublereal temp = m_thermo->temperature();
818 
819  vector_fp viscSpec(nsp);
820  vector_fp radiusSpec(nsp);
821 
822  for (size_t k = 0; k < nsp; k++) {
823  viscSpec[k] = m_viscosity[k]->getSpeciesTransProp() ;
824  radiusSpec[k] = m_hydroRadius[k]->getSpeciesTransProp() ;
825  }
826 
827  mat.resize(nsp,nsp, 0.0);
828  for (size_t i = 0; i < nsp; i++)
829  for (size_t j = 0; j < nsp; j++) {
830  mat(i,j) = (6.0 * Pi * radiusSpec[i] * viscSpec[j]) / GasConstant / temp;
831  }
832 }
833 
834 doublereal LTI_MoleFracs_ExpT::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
835 {
836 
837  size_t nsp = m_thermo->nSpecies();
838  doublereal temp = m_thermo->temperature();
839  vector_fp molefracs(nsp);
840  m_thermo->getMoleFractions(&molefracs[0]);
841 
842  doublereal value = 0;
843 
844  //if weightings are specified, use those
845  if (speciesWeight) {
846  for (size_t k = 0; k < nsp; k++) {
847  molefracs[k] = molefracs[k]*speciesWeight[k];
848  }
849  } else {
850  throw CanteraError("LTI_MoleFracs_ExpT::getMixTransProp","You should be specifying the speciesWeight");
851  }
852 
853  for (size_t i = 0; i < nsp; i++) {
854  value += speciesValues[i] * molefracs[i];
855  for (size_t j = 0; j < nsp; j++) {
856  for (size_t k = 0; k < m_Aij.size(); k++) {
857  value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i], (int) k)*exp((*m_Bij[k])(i,j)*temp);
858  }
859  }
860  }
861 
862  return value;
863 }
864 
865 
866 doublereal LTI_MoleFracs_ExpT::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
867 {
868 
869  size_t nsp = m_thermo->nSpecies();
870  doublereal temp = m_thermo->temperature();
871  vector_fp molefracs(nsp);
872  m_thermo->getMoleFractions(&molefracs[0]);
873 
874  doublereal value = 0;
875 
876  for (size_t k = 0; k < nsp; k++) {
877  molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight();
878  }
879 
880  for (size_t i = 0; i < nsp; i++) {
881  value += LTPptrs[i]->getSpeciesTransProp() * molefracs[i];
882  for (size_t j = 0; j < nsp; j++) {
883  for (size_t k = 0; k < m_Aij.size(); k++) {
884  value += molefracs[i]*molefracs[j]*(*m_Aij[k])(i,j)*pow(molefracs[i], (int) k)*exp((*m_Bij[k])(i,j)*temp);
885  }
886  }
887  }
888  return value;
889 }
890 
891 
892 
893 } //namespace Cantera
Generated by   doxygen 1.8.2