Cantera  2.0
LiquidTransport.cpp
Go to the documentation of this file.
1 /**
2  * @file LiquidTransport.cpp
3  * Mixture-averaged transport properties for ideal gas mixtures.
4  */
5 #include "cantera/thermo/ThermoPhase.h"
6 #include "cantera/transport/LiquidTransport.h"
7 
8 #include "cantera/base/utilities.h"
9 #include "cantera/transport/LiquidTransportParams.h"
10 #include "cantera/transport/TransportFactory.h"
11 #include "cantera/numerics/ctlapack.h"
12 #include "cantera/base/stringUtils.h"
13 
14 #include <iostream>
15 #include <cstdio>
16 
17 using namespace std;
18 
19 /**
20  * Mole fractions below MIN_X will be set to MIN_X when computing
21  * transport properties.
22  */
23 #define MIN_X 1.e-14
24 
25 
26 namespace Cantera
27 {
28 
29 //////////////////// class LiquidTransport methods //////////////
30 
31 
32 LiquidTransport::LiquidTransport(thermo_t* thermo, int ndim) :
33  Transport(thermo, ndim),
34  m_nsp2(0),
35  m_tmin(-1.0),
36  m_tmax(100000.),
37  m_viscMixModel(0),
38  m_ionCondMixModel(0),
39  m_lambdaMixModel(0),
40  m_diffMixModel(0),
41  m_radiusMixModel(0),
42  m_iStateMF(-1),
43  concTot_(0.0),
44  concTot_tran_(0.0),
45  dens_(0.0),
46  m_temp(-1.0),
47  m_press(-1.0),
48  m_lambda(-1.0),
49  m_viscmix(-1.0),
50  m_ionCondmix(-1.0),
51  m_mobRatMix(0),
52  m_selfDiffMix(0),
53  m_visc_mix_ok(false),
54  m_visc_temp_ok(false),
55  m_visc_conc_ok(false),
56  m_ionCond_mix_ok(false),
57  m_ionCond_temp_ok(false),
58  m_ionCond_conc_ok(false),
59  m_mobRat_mix_ok(false),
60  m_mobRat_temp_ok(false),
61  m_mobRat_conc_ok(false),
62  m_selfDiff_mix_ok(false),
63  m_selfDiff_temp_ok(false),
64  m_selfDiff_conc_ok(false),
65  m_radi_mix_ok(false),
66  m_radi_temp_ok(false),
67  m_radi_conc_ok(false),
68  m_diff_mix_ok(false),
69  m_diff_temp_ok(false),
70  m_lambda_temp_ok(false),
71  m_lambda_mix_ok(false),
72  m_mode(-1000),
73  m_debug(false),
74  m_nDim(1)
75 {
76 }
77 
78 
79 LiquidTransport::LiquidTransport(const LiquidTransport& right) :
80  Transport(right.m_thermo, right.m_nDim),
81  m_nsp2(0),
82  m_tmin(-1.0),
83  m_tmax(100000.),
84  m_viscMixModel(0),
85  m_ionCondMixModel(0),
86  m_lambdaMixModel(0),
87  m_diffMixModel(0),
88  m_radiusMixModel(0),
89  m_iStateMF(-1),
90  concTot_(0.0),
91  concTot_tran_(0.0),
92  dens_(0.0),
93  m_temp(-1.0),
94  m_press(-1.0),
95  m_lambda(-1.0),
96  m_viscmix(-1.0),
97  m_ionCondmix(-1.0),
98  m_mobRatMix(0),
99  m_selfDiffMix(0),
100  m_visc_mix_ok(false),
101  m_visc_temp_ok(false),
102  m_visc_conc_ok(false),
103  m_ionCond_mix_ok(false),
104  m_ionCond_temp_ok(false),
105  m_ionCond_conc_ok(false),
106  m_mobRat_mix_ok(false),
107  m_mobRat_temp_ok(false),
108  m_mobRat_conc_ok(false),
109  m_selfDiff_mix_ok(false),
110  m_selfDiff_temp_ok(false),
111  m_selfDiff_conc_ok(false),
112  m_radi_mix_ok(false),
113  m_radi_temp_ok(false),
114  m_radi_conc_ok(false),
115  m_diff_mix_ok(false),
116  m_diff_temp_ok(false),
117  m_lambda_temp_ok(false),
118  m_lambda_mix_ok(false),
119  m_mode(-1000),
120  m_debug(false),
121  m_nDim(1)
122 {
123  /*
124  * Use the assignment operator to do the brunt
125  * of the work for the copy constructor.
126  */
127  *this = right;
128 }
129 
130 LiquidTransport& LiquidTransport::operator=(const LiquidTransport& right)
131 {
132  if (&right == this) {
133  return *this;
134  }
135  Transport::operator=(right);
136  m_nsp2 = right.m_nsp2;
137  m_tmin = right.m_tmin;
138  m_tmax = right.m_tmax;
139  m_mw = right.m_mw;
140  m_viscTempDep_Ns = right.m_viscTempDep_Ns;
141  m_ionCondTempDep_Ns = right.m_ionCondTempDep_Ns;
142  m_mobRatTempDep_Ns = right.m_mobRatTempDep_Ns;
143  m_selfDiffTempDep_Ns = right.m_selfDiffTempDep_Ns;
144  m_lambdaTempDep_Ns = right.m_lambdaTempDep_Ns;
145  m_diffTempDep_Ns = right.m_diffTempDep_Ns;
146  m_radiusTempDep_Ns = right.m_radiusTempDep_Ns;
147  m_hydrodynamic_radius = right.m_hydrodynamic_radius;
148  m_Grad_X = right.m_Grad_X;
149  m_Grad_T = right.m_Grad_T;
150  m_Grad_V = right.m_Grad_V;
151  m_Grad_mu = right.m_Grad_mu;
152  m_bdiff = right.m_bdiff;
153  m_viscSpecies = right.m_viscSpecies;
154  m_ionCondSpecies = right.m_ionCondSpecies;
155  m_mobRatSpecies = right.m_mobRatSpecies;
156  m_selfDiffSpecies = right.m_selfDiffSpecies;
157  m_hydrodynamic_radius = right.m_hydrodynamic_radius;
158  m_lambdaSpecies = right.m_lambdaSpecies;
159  m_viscMixModel = right.m_viscMixModel;
160  m_ionCondMixModel = right.m_ionCondMixModel;
161  m_mobRatMixModel = right.m_mobRatMixModel;
162  m_selfDiffMixModel = right.m_selfDiffMixModel;
163  m_lambdaMixModel = right.m_lambdaMixModel;
164  m_diffMixModel = right.m_diffMixModel;
165  m_iStateMF = -1;
166  m_massfracs = right.m_massfracs;
167  m_massfracs_tran = right.m_massfracs_tran;
168  m_molefracs = right.m_molefracs;
169  m_molefracs_tran = right.m_molefracs_tran;
170  m_concentrations = right.m_concentrations;
171  m_actCoeff = right.m_actCoeff;
172  m_Grad_lnAC = right.m_Grad_lnAC;
173  m_chargeSpecies = right.m_chargeSpecies;
174  m_B = right.m_B;
175  m_A = right.m_A;
176  m_temp = right.m_temp;
177  m_press = right.m_press;
178  m_flux = right.m_flux;
179  m_Vdiff = right.m_Vdiff;
180  m_lambda = right.m_lambda;
181  m_viscmix = right.m_viscmix;
182  m_ionCondmix = right.m_ionCondmix;
183  m_mobRatMix = right.m_mobRatMix;
184  m_selfDiffMix = right.m_selfDiffMix;
185  m_spwork = right.m_spwork;
186  m_visc_mix_ok = false;
187  m_visc_temp_ok = false;
188  m_visc_conc_ok = false;
189  m_ionCond_mix_ok = false;
190  m_ionCond_temp_ok = false;
191  m_ionCond_conc_ok = false;
192  m_mobRat_mix_ok = false;
193  m_mobRat_temp_ok = false;
194  m_mobRat_conc_ok = false;
195  m_selfDiff_mix_ok = false;
196  m_selfDiff_temp_ok = false;
197  m_selfDiff_conc_ok = false;
198  m_radi_mix_ok = false;
199  m_radi_temp_ok = false;
200  m_radi_conc_ok = false;
201  m_diff_mix_ok = false;
202  m_diff_temp_ok = false;
203  m_lambda_temp_ok = false;
204  m_lambda_mix_ok = false;
205  m_mode = right.m_mode;
206  m_debug = right.m_debug;
207  m_nDim = right.m_nDim;
208 
209  return *this;
210 }
211 
212 
213 Transport* LiquidTransport::duplMyselfAsTransport() const
214 {
215  LiquidTransport* tr = new LiquidTransport(*this);
216  return (dynamic_cast<Transport*>(tr));
217 }
218 
219 LiquidTransport::~LiquidTransport()
220 {
221 
222  //These are constructed in TransportFactory::newLTP
223  for (size_t k = 0; k < m_nsp; k++) {
224  if (m_viscTempDep_Ns[k]) {
225  delete m_viscTempDep_Ns[k];
226  }
227  if (m_ionCondTempDep_Ns[k]) {
228  delete m_ionCondTempDep_Ns[k];
229  }
230  for (size_t l = 0; l < m_nsp; l++) {
231  if (m_selfDiffTempDep_Ns[l][k]) {
232  delete m_selfDiffTempDep_Ns[l][k];
233  }
234  }
235  for (size_t l=0; l < m_nsp2; l++) {
236  if (m_mobRatTempDep_Ns[l][k]) {
237  delete m_mobRatTempDep_Ns[l][k];
238  }
239  }
240  if (m_lambdaTempDep_Ns[k]) {
241  delete m_lambdaTempDep_Ns[k];
242  }
243  if (m_radiusTempDep_Ns[k]) {
244  delete m_radiusTempDep_Ns[k];
245  }
246  if (m_diffTempDep_Ns[k]) {
247  delete m_diffTempDep_Ns[k];
248  }
249  //These are constructed in TransportFactory::newLTI
250  if (m_selfDiffMixModel[k]) {
251  delete m_selfDiffMixModel[k];
252  }
253  }
254 
255  for (size_t k = 0; k < m_nsp2; k++) {
256  if (m_mobRatMixModel[k]) {
257  delete m_mobRatMixModel[k];
258  }
259  }
260 
261  if (m_viscMixModel) {
262  delete m_viscMixModel;
263  }
264  if (m_ionCondMixModel) {
265  delete m_ionCondMixModel;
266  }
267  if (m_lambdaMixModel) {
268  delete m_lambdaMixModel;
269  }
270  if (m_diffMixModel) {
271  delete m_diffMixModel;
272  }
273  //if ( m_radiusMixModel ) delete m_radiusMixModel;
274 
275 }
276 
277 // Initialize the transport object
278 /*
279  * Here we change all of the internal dimensions to be sufficient.
280  * We get the object ready to do property evaluations.
281  * A lot of the input required to do property evaluations is
282  * contained in the LiquidTransportParams class that is
283  * filled in TransportFactory.
284  *
285  * @param tr Transport parameters for all of the species
286  * in the phase.
287  */
288 bool LiquidTransport::initLiquid(LiquidTransportParams& tr)
289 {
290 
291  // constant substance attributes
292  m_thermo = tr.thermo;
293  tr.thermo = 0;
294  m_velocityBasis = tr.velocityBasis_;
295  m_nsp = m_thermo->nSpecies();
296  m_nsp2 = m_nsp*m_nsp;
297  m_tmin = m_thermo->minTemp();
298  m_tmax = m_thermo->maxTemp();
299 
300  // make a local copy of the molecular weights
301  m_mw.resize(m_nsp, 0.0);
302  copy(m_thermo->molecularWeights().begin(),
303  m_thermo->molecularWeights().end(), m_mw.begin());
304 
305  /*
306  * Get the input Viscosities, and stuff
307  */
308  m_viscSpecies.resize(m_nsp, 0.0);
309  m_viscTempDep_Ns.resize(m_nsp, 0);
310  m_ionCondSpecies.resize(m_nsp, 0.0);
311  m_ionCondTempDep_Ns.resize(m_nsp, 0);
312  m_mobRatTempDep_Ns.resize(m_nsp2);
313  m_mobRatMixModel.resize(m_nsp2);
314  m_mobRatSpecies.resize(m_nsp2, m_nsp, 0.0);
315  m_mobRatMix.resize(m_nsp2,0.0);
316  m_selfDiffTempDep_Ns.resize(m_nsp);
317  m_selfDiffMixModel.resize(m_nsp);
318  m_selfDiffSpecies.resize(m_nsp, m_nsp, 0.0);
319  m_selfDiffMix.resize(m_nsp,0.0);
320  for (size_t k=0; k < m_nsp; k++) {
321  m_selfDiffTempDep_Ns[k].resize(m_nsp, 0);
322  }
323  for (size_t k=0; k < m_nsp2; k++) {
324  m_mobRatTempDep_Ns[k].resize(m_nsp, 0);
325  }
326  m_lambdaSpecies.resize(m_nsp, 0.0);
327  m_lambdaTempDep_Ns.resize(m_nsp, 0);
328  m_hydrodynamic_radius.resize(m_nsp, 0.0);
329  m_radiusTempDep_Ns.resize(m_nsp, 0);
330 
331  //first populate mixing rules and indices
332  for (size_t k = 0; k < m_nsp; k++) {
333  m_selfDiffMixModel[k] = tr.selfDiffusion[k];
334  tr.selfDiffusion[k] = 0;
335  }
336  for (size_t k = 0; k < m_nsp2; k++) {
337  m_mobRatMixModel[k] = tr.mobilityRatio[k];
338  tr.mobilityRatio[k] = 0;
339  }
340 
341  //for each species, assign viscosity model and coefficients
342  for (size_t k = 0; k < m_nsp; k++) {
343  Cantera::LiquidTransportData& ltd = tr.LTData[k];
344  m_viscTempDep_Ns[k] = ltd.viscosity;
345  ltd.viscosity = 0;
346  m_ionCondTempDep_Ns[k] = ltd.ionConductivity;
347  ltd.ionConductivity = 0;
348  for (size_t j = 0; j < m_nsp2; j++) {
349  m_mobRatTempDep_Ns[j][k] = ltd.mobilityRatio[j];
350  ltd.mobilityRatio[j] = 0;
351  }
352  for (size_t j = 0; j < m_nsp; j++) {
353  m_selfDiffTempDep_Ns[j][k] = ltd.selfDiffusion[j];
354  ltd.selfDiffusion[j] = 0;
355  }
356  m_lambdaTempDep_Ns[k] = ltd.thermalCond;
357  ltd.thermalCond = 0;
358  m_radiusTempDep_Ns[k] = ltd.hydroRadius;
359  ltd.hydroRadius = 0;
360  }
361 
362 
363  /*
364  * Get the input Species Diffusivities
365  * Note that species diffusivities are not what is needed.
366  * Rather the Stefan Boltzmann interaction parameters are
367  * needed for the current model. This section may, therefore,
368  * be extraneous.
369  */
370  m_diffTempDep_Ns.resize(m_nsp, 0);
371  //for each species, assign viscosity model and coefficients
372  for (size_t k = 0; k < m_nsp; k++) {
373  Cantera::LiquidTransportData& ltd = tr.LTData[k];
374  if (ltd.speciesDiffusivity != 0) {
375  cout << "Warning: diffusion coefficient data for "
376  << m_thermo->speciesName(k)
377  << endl
378  << "in the input file is not used for LiquidTransport model."
379  << endl
380  << "LiquidTransport model uses Stefan-Maxwell interaction "
381  << endl
382  << "parameters defined in the <transport> input block."
383  << endl;
384  }
385  }
386 
387  /*
388  * Here we get interaction parameters from LiquidTransportParams
389  * that were filled in TransportFactory::getLiquidInteractionsTransportData
390  * Interaction models are provided here for viscosity, thermal conductivity,
391  * species diffusivity and hydrodynamics radius (perhaps not needed in the
392  * present class).
393  */
394 
395 
396  m_viscMixModel = tr.viscosity;
397  tr.viscosity = 0;
398 
399  m_ionCondMixModel = tr.ionConductivity;
400  tr.ionConductivity = 0;
401  //m_mobRatMixModel = tr.mobilityRatio;
402 
403  m_lambdaMixModel = tr.thermalCond;
404  tr.thermalCond = 0;
405 
406  m_diffMixModel = tr.speciesDiffusivity;
407  tr.speciesDiffusivity = 0;
408 
409 
410  m_bdiff.resize(m_nsp,m_nsp, 0.0);
411 
412  //Don't really need to update this here.
413  //It is updated in updateDiff_T()
414  m_diffMixModel->getMatrixTransProp(m_bdiff);
415 
416 
417  m_mode = tr.mode_;
418 
419  m_massfracs.resize(m_nsp, 0.0);
420  m_massfracs_tran.resize(m_nsp, 0.0);
421  m_molefracs.resize(m_nsp, 0.0);
422  m_molefracs_tran.resize(m_nsp, 0.0);
423  m_concentrations.resize(m_nsp, 0.0);
424  m_actCoeff.resize(m_nsp, 0.0);
425  m_chargeSpecies.resize(m_nsp, 0.0);
426  for (size_t i = 0; i < m_nsp; i++) {
427  m_chargeSpecies[i] = m_thermo->charge(i);
428  }
429  m_volume_spec.resize(m_nsp, 0.0);
430  m_Grad_lnAC.resize(m_nsp, 0.0);
431  m_spwork.resize(m_nsp, 0.0);
432 
433  // resize the internal gradient variables
434  m_Grad_X.resize(m_nDim * m_nsp, 0.0);
435  m_Grad_T.resize(m_nDim, 0.0);
436  m_Grad_V.resize(m_nDim, 0.0);
437  m_Grad_mu.resize(m_nDim * m_nsp, 0.0);
438 
439  m_flux.resize(m_nsp, m_nDim, 0.0);
440  m_Vdiff.resize(m_nsp, m_nDim, 0.0);
441 
442 
443  // set all flags to false
444  m_visc_mix_ok = false;
445  m_visc_temp_ok = false;
446  m_visc_conc_ok = false;
447  m_ionCond_mix_ok = false;
448  m_ionCond_temp_ok = false;
449  m_ionCond_conc_ok = false;
450  m_mobRat_mix_ok = false;
451  m_mobRat_temp_ok = false;
452  m_mobRat_conc_ok = false;
453  m_selfDiff_mix_ok = false;
454  m_selfDiff_temp_ok = false;
455  m_selfDiff_conc_ok = false;
456  m_radi_temp_ok = false;
457  m_radi_conc_ok = false;
458  m_lambda_temp_ok = false;
459  m_lambda_mix_ok = false;
460  m_diff_temp_ok = false;
461  m_diff_mix_ok = false;
462 
463  return true;
464 }
465 
466 
467 
468 /****************** viscosity ******************************/
469 
470 // Returns the viscosity of the solution
471 /*
472  * The viscosity calculation is handled by subclasses of
473  * LiquidTranInteraction as specified in the input file.
474  * These in turn employ subclasses of LTPspecies to
475  * determine the individual species viscosities.
476  */
477 doublereal LiquidTransport::viscosity()
478 {
479 
480  update_T();
481  update_C();
482 
483  if (m_visc_mix_ok) {
484  return m_viscmix;
485  }
486 
487  ////// LiquidTranInteraction method
488  m_viscmix = m_viscMixModel->getMixTransProp(m_viscTempDep_Ns);
489 
490  return m_viscmix;
491 }
492 
493 // Returns the pure species viscosities for all species
494 /*
495  * The pure species viscosities are evaluated using the
496  * appropriate subclasses of LTPspecies as specified in the
497  * input file.
498  *
499  * @param visc array of length "number of species"
500  * to hold returned viscosities.
501  */
502 void LiquidTransport::getSpeciesViscosities(doublereal* const visc)
503 {
504  update_T();
505  if (!m_visc_temp_ok) {
506  updateViscosity_T();
507  }
508  copy(m_viscSpecies.begin(), m_viscSpecies.end(), visc);
509 }
510 
511 /****************** ionConductivity ******************************/
512 
513 // Returns the ionic conductivity of the solution
514 /*
515  * The ionConductivity calculation is handled by subclasses of
516  * LiquidTranInteraction as specified in the input file.
517  * These in turn employ subclasses of LTPspecies to
518  * determine the individual species ionic conductivities.
519  */
520 doublereal LiquidTransport:: ionConductivity()
521 {
522 
523  update_T();
524  update_C();
525 
526  if (m_ionCond_mix_ok) {
527  return m_ionCondmix;
528  }
529 
530  ////// LiquidTranInteraction method
531  m_ionCondmix = m_ionCondMixModel->getMixTransProp(m_ionCondTempDep_Ns);
532 
533  return m_ionCondmix;
534 
535  /*
536  // update m_ionCondSpecies[] if necessary
537  if (!m_ionCond_temp_ok) {
538  updateIonConductivity_T();
539  }
540 
541  if (!m_ionCond_conc_ok) {
542  updateIonConductivity_C();
543  }
544  */
545 }
546 
547 // Returns the pure species ionic conductivities for all species
548 /*
549  * The pure species ionic conductivities are evaluated using the
550  * appropriate subclasses of LTPspecies as specified in the
551  * input file.
552  *
553  * @param ionCond array of length "number of species"
554  * to hold returned ionic conductivities.
555  */
556 void LiquidTransport::getSpeciesIonConductivity(doublereal* ionCond)
557 {
558  update_T();
559  if (!m_ionCond_temp_ok) {
560  updateIonConductivity_T();
561  }
562  copy(m_ionCondSpecies.begin(), m_ionCondSpecies.end(), ionCond);
563 }
564 
565 /****************** mobilityRatio ******************************/
566 
567 // Returns the mobility ratios of the solution
568 /*
569  * The mobility ratio calculation is handled by subclasses of
570  * LiquidTranInteraction as specified in the input file.
571  * These in turn employ subclasses of LTPspecies to
572  * determine the individual species mobility ratios.
573  */
574 void LiquidTransport:: mobilityRatio(doublereal* mobRat)
575 {
576 
577  update_T();
578  update_C();
579 
580  // LiquidTranInteraction method
581  if (!m_mobRat_mix_ok) {
582  for (size_t k = 0; k < m_nsp2; k++) {
583  if (m_mobRatMixModel[k]) {
584  m_mobRatMix[k] = m_mobRatMixModel[k]->getMixTransProp(m_mobRatTempDep_Ns[k]);
585  if (m_mobRatMix[k] > 0.0) {
586  m_mobRatMix[k / m_nsp + m_nsp * (k % m_nsp)] = 1.0 / m_mobRatMix[k]; // Also must be off diagonal: k%(1+n)!=0, but then m_mobRatMixModel[k] shouldn't be initialized anyway
587  }
588  }
589  }
590  }
591  for (size_t k = 0; k < m_nsp2; k++) {
592  mobRat[k] = m_mobRatMix[k];
593  }
594 }
595 
596 // Returns the pure species mobility ratios for all species
597 /*
598  * The pure species mobility ratios are evaluated using the
599  * appropriate subclasses of LTPspecies as specified in the
600  * input file.
601  *
602  * @param mobRat array of length "number of species"
603  * to hold returned mobility ratio.
604  */
605 void LiquidTransport::getSpeciesMobilityRatio(doublereal** mobRat)
606 {
607  update_T();
608  if (!m_mobRat_temp_ok) {
609  updateMobilityRatio_T();
610  }
611  for (size_t k = 0; k < m_nsp2; k++) {
612  for (size_t j = 0; j < m_nsp; j++) {
613  mobRat[k][j] = m_mobRatSpecies(k,j);
614  }
615  }
616 }
617 //====================================================================================================================
618 // Returns the self diffusion coefficients of the species in the phase
619 /*
620  * The self diffusion coefficient is the diffusion coefficient of a tracer species
621  * at the current temperature and composition of the species. Therefore,
622  * the dilute limit of transport is assumed for the tracer species.
623  * The effective formula may be calculated from the stefan-maxwell formulation by
624  * adding another row for the tracer species, assigning all D's to be equal
625  * to the respective species D's, and then taking the limit as the
626  * tracer species mole fraction goes to zero. The corresponding flux equation
627  * for the tracer species k in units of kmol m-2 s-1 is.
628  *
629  * \f[
630  * J_k = - D^{sd}_k \frac{C_k}{R T} \nabla \mu_k
631  * \f]
632  *
633  * The derivative is taken at constant T and P.
634  *
635  * The self diffusion calculation is handled by subclasses of
636  * LiquidTranInteraction as specified in the input file.
637  * These in turn employ subclasses of LTPspecies to
638  * determine the individual species self diffusion coeffs.
639  *
640  * @param selfDiff Vector of self-diffusion coefficients
641  * Length = number of species in phase
642  * units = m**2 s-1
643  */
644 void LiquidTransport::selfDiffusion(doublereal* const selfDiff)
645 {
646  update_T();
647  update_C();
648  if (!m_selfDiff_mix_ok) {
649  for (size_t k = 0; k < m_nsp; k++) {
650  m_selfDiffMix[k] = m_selfDiffMixModel[k]->getMixTransProp(m_selfDiffTempDep_Ns[k]);
651  }
652  }
653  for (size_t k = 0; k < m_nsp; k++) {
654  selfDiff[k] = m_selfDiffMix[k];
655  }
656 }
657 //====================================================================================================================
658 // Returns the pure species self diffusion for all species
659 /*
660  * The pure species self diffusion coeffs are evaluated using the
661  * appropriate subclasses of LTPspecies as specified in the
662  * input file.
663  *
664  * @param selfDiff array of size "number of species"^2
665  * to hold returned self diffusion.
666  */
667 void LiquidTransport::getSpeciesSelfDiffusion(doublereal** selfDiff)
668 {
669  update_T();
670  if (!m_selfDiff_temp_ok) {
671  updateSelfDiffusion_T();
672  }
673  for (size_t k=0; k<m_nsp; k++) {
674  for (size_t j=0; j < m_nsp; j++) {
675  selfDiff[k][j] = m_selfDiffSpecies(k,j);
676  }
677  }
678 }
679 
680 //===============================================================
681 // Returns the hydrodynamic radius for all species
682 /*
683  * The species hydrodynamic radii are evaluated using the
684  * appropriate subclasses of LTPspecies as specified in the
685  * input file.
686  *
687  * @param radius array of length "number of species"
688  * to hold returned radii.
689  */
690 void LiquidTransport::getSpeciesHydrodynamicRadius(doublereal* const radius)
691 {
692  update_T();
693  if (!m_radi_temp_ok) {
694  updateHydrodynamicRadius_T();
695  }
696  copy(m_hydrodynamic_radius.begin(), m_hydrodynamic_radius.end(), radius);
697 
698 }
699 
700 //================================================================
701 
702 // Return the thermal conductivity of the solution
703 /*
704  * The thermal conductivity calculation is handled by subclasses of
705  * LiquidTranInteraction as specified in the input file.
706  * These in turn employ subclasses of LTPspecies to
707  * determine the individual species thermal condictivities.
708  */
709 doublereal LiquidTransport::thermalConductivity()
710 {
711 
712  update_T();
713  update_C();
714 
715  if (!m_lambda_mix_ok) {
716  m_lambda = m_lambdaMixModel->getMixTransProp(m_lambdaTempDep_Ns);
717  m_cond_mix_ok = true;
718  }
719 
720  return m_lambda;
721 }
722 
723 
724 /****************** thermal diffusion coefficients ************/
725 
726 // Return the thermal diffusion coefficients
727 /*
728  * These are all zero for this simple implementaion
729  *
730  * @param dt thermal diffusion coefficients
731  */
732 void LiquidTransport::getThermalDiffCoeffs(doublereal* const dt)
733 {
734  for (size_t k = 0; k < m_nsp; k++) {
735  dt[k] = 0.0;
736  }
737 }
738 
739 /******************* binary diffusion coefficients **************/
740 
741 
742 // Returns the binary diffusion coefficients
743 /*
744  * The binary diffusion coefficients are specified in the input
745  * file through the LiquidTransportInteractions class. These
746  * are the binary interaction coefficients employed in the
747  * Stefan-Maxwell equation.
748  *
749  * @param ld number of species in system
750  * @param d vector of binary diffusion coefficients
751  * units = m2 s-1. length = ld*ld = (number of species)^2
752  */
753 void LiquidTransport::getBinaryDiffCoeffs(size_t ld, doublereal* d)
754 {
755  if (ld != m_nsp)
756  throw CanteraError("LiquidTransport::getBinaryDiffCoeffs",
757  "First argument does not correspond to number of species in model.\nDiff Coeff matrix may be misdimensioned");
758  update_T();
759 
760  // if necessary, evaluate the binary diffusion coefficients
761  // from the polynomial fits
762  if (!m_diff_temp_ok) {
763  updateDiff_T();
764  }
765 
766  for (size_t i = 0; i < m_nsp; i++) {
767  for (size_t j = 0; j < m_nsp; j++) {
768  //if (!( ( m_bdiff(i,j) > 0.0 ) | ( m_bdiff(i,j) < 0.0 ))){
769  // throw CanteraError("LiquidTransport::getBinaryDiffCoeffs ",
770  // "m_bdiff has zero entry in non-diagonal.");}
771  d[ld*j + i] = 1.0 / m_bdiff(i,j);
772 
773  }
774  }
775 }
776 //================================================================================================
777 // Get the Electrical mobilities (m^2/V/s).
778 /*
779  * The electrical mobilities are not well defined
780  * in the context of LiquidTransport because the Stefan Maxwell
781  * equation is solved. Here the electrical mobilities
782  * are calculated from the mixture-averaged
783  * diffusion coefficients through a call to getMixDiffCoeffs()
784  * using the Einstein relation
785  *
786  * \f[
787  * \mu^e_k = \frac{F D_k}{R T}
788  * \f]
789  *
790  * Note that this call to getMixDiffCoeffs() requires
791  * a solve of the Stefan Maxwell equation making this
792  * determination of the mixture averaged diffusion coefficients
793  * a {\em slow} method for obtaining diffusion coefficients.
794  *
795  * Also note that the Stefan Maxwell solve will be based upon
796  * the thermodynamic state (including gradients) most recently
797  * set. Gradients can be set specifically using set_Grad_V,
798  * set_Grad_X and set_Grad_T or through calls to
799  * getSpeciesFluxes, getSpeciesFluxesES, getSpeciesVdiff,
800  * getSpeciesVdiffES, etc.
801  *
802  * @param mobil_e Returns the electrical mobilities of
803  * the species in array \c mobil_e. The array must be
804  * dimensioned at least as large as the number of species.
805  */
806 void LiquidTransport::getMobilities(doublereal* const mobil)
807 {
808  getMixDiffCoeffs(DATA_PTR(m_spwork));
809  doublereal c1 = ElectronCharge / (Boltzmann * m_temp);
810  for (size_t k = 0; k < m_nsp; k++) {
811  mobil[k] = c1 * m_spwork[k];
812  }
813 }
814 //================================================================================================
815 // Get the fluid mobilities (s kmol/kg).
816 /*
817  * The fluid mobilities are not well defined
818  * in the context of LiquidTransport because the Stefan Maxwell
819  * equation is solved. Here the fluid mobilities
820  * are calculated from the mixture-averaged
821  * diffusion coefficients through a call to getMixDiffCoeffs()
822  * using the Einstein relation
823  *
824  * \f[
825  * \mu^f_k = \frac{D_k}{R T}
826  * \f]
827  *
828  * Note that this call to getMixDiffCoeffs() requires
829  * a solve of the Stefan Maxwell equation making this
830  * determination of the mixture averaged diffusion coefficients
831  * a {\em slow} method for obtaining diffusion coefficients.
832  *
833  * Also note that the Stefan Maxwell solve will be based upon
834  * the thermodynamic state (including gradients) most recently
835  * set. Gradients can be set specifically using set_Grad_V,
836  * set_Grad_X and set_Grad_T or through calls to
837  * getSpeciesFluxes, getSpeciesFluxesES, getSpeciesVdiff,
838  * getSpeciesVdiffES, etc.
839  *
840  * @param mobil_f Returns the fluid mobilities of
841  * the species in array \c mobil_f. The array must be
842  * dimensioned at least as large as the number of species.
843  */
844 void LiquidTransport::getFluidMobilities(doublereal* const mobil_f)
845 {
846  getMixDiffCoeffs(DATA_PTR(m_spwork));
847  doublereal c1 = 1.0 / (GasConstant * m_temp);
848  for (size_t k = 0; k < m_nsp; k++) {
849  mobil_f[k] = c1 * m_spwork[k];
850  }
851 }
852 //==============================================================
853 // Specify the value of the gradient of the temperature
854 /*
855  * @param grad_T Gradient of the temperature (length num dimensions);
856  */
857 void LiquidTransport::set_Grad_T(const doublereal* const grad_T)
858 {
859  for (size_t a = 0; a < m_nDim; a++) {
860  m_Grad_T[a] = grad_T[a];
861  }
862 }
863 //==============================================================
864 // Specify the value of the gradient of the voltage
865 /*
866  *
867  * @param grad_V Gradient of the voltage (length num dimensions);
868  */
869 void LiquidTransport::set_Grad_V(const doublereal* const grad_V)
870 {
871  for (size_t a = 0; a < m_nDim; a++) {
872  m_Grad_V[a] = grad_V[a];
873  }
874 }
875 //==============================================================
876 // Specify the value of the gradient of the MoleFractions
877 /*
878  *
879  * @param grad_X Gradient of the mole fractions(length nsp * num dimensions);
880  */
881 void LiquidTransport::set_Grad_X(const doublereal* const grad_X)
882 {
883  size_t itop = m_nDim * m_nsp;
884  for (size_t i = 0; i < itop; i++) {
885  m_Grad_X[i] = grad_X[i];
886  }
887 }
888 //==============================================================
889 
890 // Compute the mixture electrical conductivity from
891 // the Stefan-Maxwell equation.
892 /*
893  * To compute the mixture electrical conductance, the Stefan
894  * Maxwell equation is solved for zero species gradients and
895  * for unit potential gradient, \f$ \nabla V \f$.
896  * The species fluxes are converted to current by summing over
897  * the charge-weighted fluxes according to
898  * \f[
899  * \vec{i} = \sum_{i} z_i F \rho \vec{V_i} / W_i
900  * \f]
901  * where \f$ z_i \f$ is the charge on species i,
902  * \f$ F \f$ is Faradays constant, \f$ \rho \f$ is the density,
903  * \f$ W_i \f$ is the molecular mass of species i.
904  * The conductance, \f$ \kappa \f$ is obtained from
905  * \f[
906  * \kappa = \vec{i} / \nabla V.
907  * \f]
908  */
909 doublereal LiquidTransport::getElectricConduct()
910 {
911  doublereal gradT = 0.0;
912  vector_fp gradX(m_nDim * m_nsp);
913  vector_fp gradV(m_nDim);
914  for (size_t i = 0; i < m_nDim; i++) {
915  for (size_t k = 0; k < m_nsp; k++) {
916  gradX[ i*m_nDim + k] = 0.0;
917  }
918  gradV[i] = 1.0;
919  }
920 
921  set_Grad_T(&gradT);
922  set_Grad_X(&gradX[0]);
923  set_Grad_V(&gradV[0]);
924 
925  vector_fp fluxes(m_nsp * m_nDim);
926  doublereal current;
927 
928  getSpeciesFluxesExt(m_nDim, &fluxes[0]);
929 
930  //sum over species charges, fluxes, Faraday to get current
931  // Since we want the scalar conductivity, we need only consider one-dim
932  for (size_t i = 0; i < 1; i++) {
933  current = 0.0;
934  for (size_t k = 0; k < m_nsp; k++) {
935  current += m_chargeSpecies[k] * Faraday * fluxes[k] / m_mw[k];
936  }
937  //divide by unit potential gradient
938  current /= - gradV[i];
939  }
940  return current;
941 }
942 
943 // Compute the electric current density in A/m^2
944 /*
945  * The electric current is computed first by computing the
946  * species diffusive fluxes using the Stefan Maxwell solution
947  * and then the current, \f$ \vec{i} \f$ by summing over
948  * the charge-weighted fluxes according to
949  * \f[
950  * \vec{i} = \sum_{i} z_i F \rho \vec{V_i} / W_i
951  * \f]
952  * where \f$ z_i \f$ is the charge on species i,
953  * \f$ F \f$ is Faradays constant, \f$ \rho \f$ is the density,
954  * \f$ W_i \f$ is the molecular mass of species i.
955  *
956  * @param ndim The number of spatial dimensions (1, 2, or 3).
957  * @param grad_T The temperature gradient (ignored in this model).
958  * @param ldx Leading dimension of the grad_X array.
959  * @param grad_T The temperature gradient (ignored in this model).
960  * @param ldf Leading dimension of the grad_V and current vectors.
961  * @param grad_V The electrostatic potential gradient.
962  * @param current The electric current in A/m^2.
963  */
964 void LiquidTransport::getElectricCurrent(int ndim,
965  const doublereal* grad_T,
966  int ldx,
967  const doublereal* grad_X,
968  int ldf,
969  const doublereal* grad_V,
970  doublereal* current)
971 {
972 
973  set_Grad_T(grad_T);
974  set_Grad_X(grad_X);
975  set_Grad_V(grad_V);
976 
977  vector_fp fluxes(m_nsp * m_nDim);
978 
979  getSpeciesFluxesExt(ldf, &fluxes[0]);
980 
981  //sum over species charges, fluxes, Faraday to get current
982  for (size_t i = 0; i < m_nDim; i++) {
983  current[i] = 0.0;
984  for (size_t k = 0; k < m_nsp; k++) {
985  current[i] += m_chargeSpecies[k] * Faraday * fluxes[k] / m_mw[k];
986  }
987  //divide by unit potential gradient
988  }
989 }
990 
991 // Get the species diffusive velocities wrt to
992 // the averaged velocity,
993 // given the gradients in mole fraction and temperature
994 /*
995  * The average velocity can be computed on a mole-weighted
996  * or mass-weighted basis, or the diffusion velocities may
997  * be specified as relative to a specific species (i.e. a
998  * solvent) all according to the velocityBasis input parameter.
999  *
1000  * Units for the returned fluxes are kg m-2 s-1.
1001  *
1002  * @param ndim Number of dimensions in the flux expressions
1003  * @param grad_T Gradient of the temperature
1004  * (length = ndim)
1005  * @param ldx Leading dimension of the grad_X array
1006  * (usually equal to m_nsp but not always)
1007  * @param grad_X Gradients of the mole fraction
1008  * Flat vector with the m_nsp in the inner loop.
1009  * length = ldx * ndim
1010  * @param ldf Leading dimension of the fluxes array
1011  * (usually equal to m_nsp but not always)
1012  * @param Vdiff Output of the diffusive velocities.
1013  * Flat vector with the m_nsp in the inner loop.
1014  * length = ldx * ndim
1015  */
1016 void LiquidTransport::getSpeciesVdiff(size_t ndim,
1017  const doublereal* grad_T,
1018  int ldx, const doublereal* grad_X,
1019  int ldf, doublereal* Vdiff)
1020 {
1021  set_Grad_T(grad_T);
1022  set_Grad_X(grad_X);
1023  getSpeciesVdiffExt(ldf, Vdiff);
1024 }
1025 
1026 /*
1027  * @param ndim The number of spatial dimensions (1, 2, or 3).
1028  * @param grad_T The temperature gradient (ignored in this model).
1029  * @param ldx Leading dimension of the grad_X array.
1030  * The diffusive mass flux of species \e k is computed from
1031  *
1032  * \f[
1033  * \vec{j}_k = -n M_k D_k \nabla X_k.
1034  * \f]
1035  */
1036 void LiquidTransport::getSpeciesVdiffES(size_t ndim,
1037  const doublereal* grad_T,
1038  int ldx,
1039  const doublereal* grad_X,
1040  int ldf,
1041  const doublereal* grad_V,
1042  doublereal* Vdiff)
1043 {
1044  set_Grad_T(grad_T);
1045  set_Grad_X(grad_X);
1046  set_Grad_V(grad_V);
1047  getSpeciesVdiffExt(ldf, Vdiff);
1048 }
1049 
1050 // Return the species diffusive mass fluxes wrt to
1051 // the averaged velocity in [kmol/m^2/s].
1052 /*
1053  *
1054  * The diffusive mass flux of species \e k is computed
1055  * using the Stefan-Maxwell equation
1056  * \f[
1057  * X_i \nabla \mu_i
1058  * = RT \sum_i \frac{X_i X_j}{D_{ij}}
1059  * ( \vec{V}_j - \vec{V}_i )
1060  * \f]
1061  * to determine the diffusion velocity and
1062  * \f[
1063  * \vec{N}_i = C_T X_i \vec{V}_i
1064  * \f]
1065  * to determine the diffusion flux. Here \f$ C_T \f$ is the
1066  * total concentration of the mixture [kmol/m^3], \f$ D_{ij} \f$
1067  * are the Stefa-Maxwell interaction parameters in [m^2/s],
1068  * \f$ \vec{V}_{i} \f$ is the diffusion velocity of species \e i,
1069  * \f$ \mu_i \f$ is the electrochemical potential of species \e i.
1070  *
1071  * Note that for this method, there is no argument for the
1072  * gradient of the electric potential (voltage). Electric
1073  * potential gradients can be set with set_Grad_V() or
1074  * method getSpeciesFluxesES() can be called.x
1075  *
1076  * The diffusion velocity is relative to an average velocity
1077  * that can be computed on a mole-weighted
1078  * or mass-weighted basis, or the diffusion velocities may
1079  * be specified as relative to a specific species (i.e. a
1080  * solvent) all according to the \verbatim <velocityBasis>
1081  * \endverbatim input parameter.
1082 
1083  * @param ndim The number of spatial dimensions (1, 2, or 3).
1084  * @param grad_T The temperature gradient (ignored in this model).
1085  * (length = ndim)
1086  * @param ldx Leading dimension of the grad_X array.
1087  * (usually equal to m_nsp but not always)
1088  * @param grad_X Gradients of the mole fraction
1089  * Flat vector with the m_nsp in the inner loop.
1090  * length = ldx * ndim
1091  * @param ldf Leading dimension of the fluxes array
1092  * (usually equal to m_nsp but not always)
1093  * @param grad_Phi Gradients of the electrostatic potential
1094  * length = ndim
1095  * @param fluxes Output of the diffusive mass fluxes
1096  * Flat vector with the m_nsp in the inner loop.
1097  * length = ldx * ndim
1098  */
1099 void LiquidTransport::getSpeciesFluxes(size_t ndim,
1100  const doublereal* const grad_T,
1101  size_t ldx, const doublereal* const grad_X,
1102  size_t ldf, doublereal* const fluxes)
1103 {
1104  set_Grad_T(grad_T);
1105  set_Grad_X(grad_X);
1106  getSpeciesFluxesExt(ldf, fluxes);
1107 }
1108 
1109 // Return the species diffusive mass fluxes wrt to
1110 // the averaged velocity in [kmol/m^2/s].
1111 /*
1112  *
1113  * The diffusive mass flux of species \e k is computed
1114  * using the Stefan-Maxwell equation
1115  * \f[
1116  * X_i \nabla \mu_i
1117  * = RT \sum_i \frac{X_i X_j}{D_{ij}}
1118  * ( \vec{V}_j - \vec{V}_i )
1119  * \f]
1120  * to determine the diffusion velocity and
1121  * \f[
1122  * \vec{N}_i = C_T X_i \vec{V}_i
1123  * \f]
1124  * to determine the diffusion flux. Here \f$ C_T \f$ is the
1125  * total concentration of the mixture [kmol/m^3], \f$ D_{ij} \f$
1126  * are the Stefa-Maxwell interaction parameters in [m^2/s],
1127  * \f$ \vec{V}_{i} \f$ is the diffusion velocity of species \e i,
1128  * \f$ \mu_i \f$ is the electrochemical potential of species \e i.
1129  *
1130  * The diffusion velocity is relative to an average velocity
1131  * that can be computed on a mole-weighted
1132  * or mass-weighted basis, or the diffusion velocities may
1133  * be specified as relative to a specific species (i.e. a
1134  * solvent) all according to the \verbatim <velocityBasis>
1135  * \endverbatim input parameter.
1136 
1137  * @param ndim The number of spatial dimensions (1, 2, or 3).
1138  * @param grad_T The temperature gradient (ignored in this model).
1139  * (length = ndim)
1140  * @param ldx Leading dimension of the grad_X array.
1141  * (usually equal to m_nsp but not always)
1142  * @param grad_X Gradients of the mole fraction
1143  * Flat vector with the m_nsp in the inner loop.
1144  * length = ldx * ndim
1145  * @param ldf Leading dimension of the fluxes array
1146  * (usually equal to m_nsp but not always)
1147  * @param grad_Phi Gradients of the electrostatic potential
1148  * length = ndim
1149  * @param fluxes Output of the diffusive mass fluxes
1150  * Flat vector with the m_nsp in the inner loop.
1151  * length = ldx * ndim
1152  */
1153 void LiquidTransport::getSpeciesFluxesES(size_t ndim,
1154  const doublereal* grad_T,
1155  int ldx,
1156  const doublereal* grad_X,
1157  int ldf,
1158  const doublereal* grad_V,
1159  doublereal* fluxes)
1160 {
1161  set_Grad_T(grad_T);
1162  set_Grad_X(grad_X);
1163  set_Grad_V(grad_V);
1164  getSpeciesFluxesExt(ldf, fluxes);
1165 }
1166 
1167 // Return the species diffusive velocities relative to
1168 // the averaged velocity.
1169 /*
1170  * This method acts similarly to getSpeciesVdiffES() but
1171  * requires all gradients to be preset using methods
1172  * set_Grad_X(), set_Grad_V(), set_Grad_T().
1173  * See the documentation of getSpeciesVdiffES() for details.
1174  *
1175  * @param ldf Leading dimension of the Vdiff array.
1176  * @param Vdiff Output of the diffusive velocities.
1177  * Flat vector with the m_nsp in the inner loop.
1178  * length = ldx * ndim
1179  */
1180 void LiquidTransport::getSpeciesVdiffExt(size_t ldf, doublereal* Vdiff)
1181 {
1182  stefan_maxwell_solve();
1183 
1184  for (size_t n = 0; n < m_nDim; n++) {
1185  for (size_t k = 0; k < m_nsp; k++) {
1186  Vdiff[n*ldf + k] = m_Vdiff(k,n);
1187  }
1188  }
1189 }
1190 
1191 // Return the species diffusive fluxes relative to
1192 // the averaged velocity.
1193 /*
1194  * This method acts similarly to getSpeciesFluxesES() but
1195  * requires all gradients to be preset using methods
1196  * set_Grad_X(), set_Grad_V(), set_Grad_T().
1197  * See the documentation of getSpeciesFluxesES() for details.
1198  *
1199  * units = kg/m2/s
1200  *
1201  * @param ldf Leading dimension of the Vdiff array.
1202  * @param fluxes Output of the diffusive fluxes.
1203  * Flat vector with the m_nsp in the inner loop.
1204  * length = ldx * ndim
1205  */
1206 void LiquidTransport::getSpeciesFluxesExt(size_t ldf, doublereal* fluxes)
1207 {
1208  stefan_maxwell_solve();
1209 
1210  for (size_t n = 0; n < m_nDim; n++) {
1211  for (size_t k = 0; k < m_nsp; k++) {
1212  fluxes[n*ldf + k] = m_flux(k,n);
1213  }
1214  }
1215 }
1216 
1217 // Get the Mixture diffusion coefficients [m^2/s]
1218 /*
1219  * The mixture diffusion coefficients are not well defined
1220  * in the context of LiquidTransport because the Stefan Maxwell
1221  * equation is solved. Here the mixture diffusion coefficients
1222  * are defined according to Ficks law:
1223  * \f[
1224  * X_i \vec{V_i} = -D_i \nabla X_i.
1225  * \f]
1226  * Solving Ficks Law for \f$ D_i \f$ gives a mixture diffusion
1227  * coefficient
1228  * \f[
1229  * D_i = - X_i \vec{V_i} / ( \nabla X_i ).
1230  * \f]
1231  * If \f$ \nabla X_i = 0 \f$ this is undefined and the
1232  * nonsensical value -1 is returned.
1233  *
1234  * Note that this evaluation of \f$ \vec{V_i} \f$ requires
1235  * a solve of the Stefan Maxwell equation making this
1236  * determination of the mixture averaged diffusion coefficients
1237  * a {\em slow} method for obtaining diffusion coefficients.
1238  *
1239  * Also note that the Stefan Maxwell solve will be based upon
1240  * the thermodynamic state (including gradients) most recently
1241  * set. Gradients can be set specifically using set_Grad_V,
1242  * set_Grad_X and set_Grad_T or through calls to
1243  * getSpeciesFluxes, getSpeciesFluxesES, getSpeciesVdiff,
1244  * getSpeciesVdiffES, etc.
1245  *
1246  * @param d vector of mixture diffusion coefficients
1247  * units = m2 s-1. length = number of species
1248  */
1249 void LiquidTransport::getMixDiffCoeffs(doublereal* const d)
1250 {
1251 
1252  stefan_maxwell_solve();
1253 
1254  for (size_t n = 0; n < m_nDim; n++) {
1255  for (size_t k = 0; k < m_nsp; k++) {
1256  if (m_Grad_X[n*m_nsp + k] != 0.0) {
1257  d[n*m_nsp + k] = - m_Vdiff(k,n) * m_molefracs[k]
1258  / m_Grad_X[n*m_nsp + k];
1259  } else {
1260  //avoid divide by zero with nonsensical response
1261  d[n*m_nsp + k] = - 1.0;
1262  }
1263  }
1264  }
1265 }
1266 
1267 
1268 // Handles the effects of changes in the Temperature, internally
1269 // within the object.
1270 /*
1271  * This is called whenever a transport property is
1272  * requested.
1273  * The first task is to check whether the temperature has changed
1274  * since the last call to update_T().
1275  * If it hasn't then an immediate return is carried out.
1276  *
1277  * @internal
1278  */
1279 bool LiquidTransport::update_T()
1280 {
1281  // First make a decision about whether we need to recalculate
1282  doublereal t = m_thermo->temperature();
1283  if (t == m_temp) {
1284  return false;
1285  }
1286 
1287  // Next do a reality check on temperature value
1288  if (t < 0.0) {
1289  throw CanteraError("LiquidTransport::update_T()",
1290  "negative temperature "+fp2str(t));
1291  }
1292 
1293  // Compute various direct functions of temperature
1294  m_temp = t;
1295 
1296  // temperature has changed so temp flags are flipped
1297  m_visc_temp_ok = false;
1298  m_ionCond_temp_ok = false;
1299  m_mobRat_temp_ok = false;
1300  m_selfDiff_temp_ok = false;
1301  m_radi_temp_ok = false;
1302  m_diff_temp_ok = false;
1303  m_lambda_temp_ok = false;
1304 
1305  // temperature has changed, so polynomial temperature
1306  // interpolations will need to be reevaluated.
1307  // This means that many concentration
1308  m_visc_conc_ok = false;
1309  m_ionCond_conc_ok = false;
1310  m_mobRat_conc_ok = false;
1311  m_selfDiff_conc_ok = false;
1312 
1313  // Mixture stuff needs to be evaluated
1314  m_visc_mix_ok = false;
1315  m_ionCond_mix_ok = false;
1316  m_mobRat_mix_ok = false;
1317  m_selfDiff_mix_ok = false;
1318  m_diff_mix_ok = false;
1319  m_lambda_mix_ok = false; //(don't need it because a lower lvl flag is set
1320  return true;
1321 }
1322 
1323 
1324 // Handles the effects of changes in the mixture concentration
1325 /*
1326  * This is called for every interface call to check whether
1327  * the concentrations have changed. Concentrations change
1328  * whenever the pressure or the mole fraction has changed.
1329  * If it has changed, the recalculations should be done.
1330  *
1331  * Note this should be a lightweight function since it's
1332  * part of all of the interfaces.
1333  *
1334  * @internal
1335  */
1336 bool LiquidTransport::update_C()
1337 {
1338  // If the pressure has changed then the concentrations
1339  // have changed.
1340  doublereal pres = m_thermo->pressure();
1341  bool qReturn = true;
1342  if (pres != m_press) {
1343  qReturn = false;
1344  m_press = pres;
1345  }
1346  int iStateNew = m_thermo->stateMFNumber();
1347  if (iStateNew != m_iStateMF) {
1348  qReturn = false;
1349  m_thermo->getMassFractions(DATA_PTR(m_massfracs));
1350  m_thermo->getMoleFractions(DATA_PTR(m_molefracs));
1351  m_thermo->getConcentrations(DATA_PTR(m_concentrations));
1352  concTot_ = 0.0;
1353  concTot_tran_ = 0.0;
1354  for (size_t k = 0; k < m_nsp; k++) {
1355  m_molefracs[k] = std::max(0.0, m_molefracs[k]);
1356  m_molefracs_tran[k] = std::max(MIN_X, m_molefracs[k]);
1357  m_massfracs_tran[k] = std::max(MIN_X, m_massfracs[k]);
1358  concTot_tran_ += m_molefracs_tran[k];
1359  concTot_ += m_concentrations[k];
1360  }
1361  dens_ = m_thermo->density();
1362  meanMolecularWeight_ = m_thermo->meanMolecularWeight();
1363  concTot_tran_ *= concTot_;
1364  }
1365  if (qReturn) {
1366  return false;
1367  }
1368 
1369  // signal that concentration-dependent quantities will need to
1370  // be recomputed before use, and update the local mole
1371  // fractions.
1372  m_visc_conc_ok = false;
1373  m_ionCond_conc_ok = false;
1374  m_mobRat_conc_ok = false;
1375  m_selfDiff_conc_ok = false;
1376 
1377  // Mixture stuff needs to be evaluated
1378  m_visc_mix_ok = false;
1379  m_ionCond_mix_ok = false;
1380  m_mobRat_mix_ok = false;
1381  m_selfDiff_mix_ok = false;
1382  m_diff_mix_ok = false;
1383  m_lambda_mix_ok = false;
1384 
1385  return true;
1386 }
1387 
1388 /*************************************************************************
1389  *
1390  * methods to update species temperature-dependent properties
1391  *
1392  *************************************************************************/
1393 
1394 /**
1395  * Update the temperature-dependent parts of the species
1396  * thermal conductivity internally using calls to the
1397  * appropriate LTPspecies subclass.
1398  */
1399 void LiquidTransport::updateCond_T()
1400 {
1401  for (size_t k = 0; k < m_nsp; k++) {
1402  m_lambdaSpecies[k] = m_lambdaTempDep_Ns[k]->getSpeciesTransProp() ;
1403  }
1404  m_lambda_temp_ok = true;
1405  m_lambda_mix_ok = false;
1406 }
1407 
1408 
1409 // Update the binary Stefan-Maxwell diffusion coefficients
1410 // wrt T using calls to the appropriate LTPspecies subclass
1411 void LiquidTransport::updateDiff_T()
1412 {
1413  m_diffMixModel->getMatrixTransProp(m_bdiff);
1414  m_diff_temp_ok = true;
1415  m_diff_mix_ok = false;
1416 }
1417 
1418 
1419 // Update the pure-species viscosities functional dependence on concentration.
1420 void LiquidTransport::updateViscosities_C()
1421 {
1422  m_visc_conc_ok = true;
1423 }
1424 
1425 
1426 /*
1427  * Updates the array of pure species viscosities internally
1428  * using calls to the appropriate LTPspecies subclass.
1429  * The flag m_visc_ok is set to true.
1430  *
1431  * Note that for viscosity, a positive activation energy
1432  * corresponds to the typical case of a positive argument
1433  * to the exponential so that the Arrhenius expression is
1434  *
1435  * \f[
1436  * \mu = A T^n \exp( + E / R T )
1437  * \f]
1438  */
1439 void LiquidTransport::updateViscosity_T()
1440 {
1441  for (size_t k = 0; k < m_nsp; k++) {
1442  m_viscSpecies[k] = m_viscTempDep_Ns[k]->getSpeciesTransProp() ;
1443  }
1444  m_visc_temp_ok = true;
1445  m_visc_mix_ok = false;
1446 }
1447 
1448 
1449 // Update the pure-species ionic conductivities functional dependence on concentration.
1450 void LiquidTransport::updateIonConductivity_C()
1451 {
1452  m_ionCond_conc_ok = true;
1453 }
1454 
1455 
1456 /*
1457  * Updates the array of pure species ionic conductivities internally
1458  * using calls to the appropriate LTPspecies subclass.
1459  * The flag m_ionCond_ok is set to true.
1460  */
1461 void LiquidTransport::updateIonConductivity_T()
1462 {
1463  for (size_t k = 0; k < m_nsp; k++) {
1464  m_ionCondSpecies[k] = m_ionCondTempDep_Ns[k]->getSpeciesTransProp() ;
1465  }
1466  m_ionCond_temp_ok = true;
1467  m_ionCond_mix_ok = false;
1468 }
1469 
1470 
1471 // Update the pure-species mobility ratios functional dependence on concentration.
1472 void LiquidTransport::updateMobilityRatio_C()
1473 {
1474  m_mobRat_conc_ok = true;
1475 }
1476 
1477 /*
1478  * Updates the array of pure species mobility ratios internally
1479  * using calls to the appropriate LTPspecies subclass.
1480  * The flag m_mobRat_ok is set to true.
1481  */
1482 void LiquidTransport::updateMobilityRatio_T()
1483 {
1484  for (size_t k = 0; k < m_nsp2; k++) {
1485  for (size_t j = 0; j < m_nsp; j++) {
1486  m_mobRatSpecies(k,j) = m_mobRatTempDep_Ns[k][j]->getSpeciesTransProp();
1487  }
1488  }
1489  m_mobRat_temp_ok = true;
1490  m_mobRat_mix_ok = false;
1491 }
1492 
1493 
1494 // Update the pure-species self diffusion functional dependence on concentration.
1495 void LiquidTransport::updateSelfDiffusion_C()
1496 {
1497  m_selfDiff_conc_ok = true;
1498 }
1499 
1500 
1501 /*
1502  * Updates the array of pure species self diffusion internally
1503  * using calls to the appropriate LTPspecies subclass.
1504  * The flag m_selfDiff_ok is set to true.
1505  */
1506 void LiquidTransport::updateSelfDiffusion_T()
1507 {
1508  for (size_t k = 0; k < m_nsp2; k++) {
1509  for (size_t j = 0; j < m_nsp; j++) {
1510  m_selfDiffSpecies(k,j) = m_selfDiffTempDep_Ns[k][j]->getSpeciesTransProp() ;
1511  }
1512  }
1513  m_selfDiff_temp_ok = true;
1514  m_selfDiff_mix_ok = false;
1515 }
1516 //=============================================================================================================
1517 void LiquidTransport::updateHydrodynamicRadius_C()
1518 {
1519  m_radi_conc_ok = true;
1520 }
1521 //=============================================================================================================
1522 // Update the temperature-dependent hydrodynamic radius terms
1523 // for each species internally using calls to the
1524 // appropriate LTPspecies subclass
1525 void LiquidTransport::updateHydrodynamicRadius_T()
1526 {
1527  for (size_t k = 0; k < m_nsp; k++) {
1528  m_hydrodynamic_radius[k] = m_radiusTempDep_Ns[k]->getSpeciesTransProp() ;
1529  }
1530  m_radi_temp_ok = true;
1531  m_radi_mix_ok = false;
1532 }
1533 
1534 void LiquidTransport::update_Grad_lnAC()
1535 {
1536  doublereal grad_T;
1537  vector_fp grad_lnAC(m_nsp), grad_X(m_nsp);
1538  // IonsFromNeutralVPSSTP * tempIons = dynamic_cast<IonsFromNeutralVPSSTP *> m_thermo;
1539  //MargulesVPSSTP * tempMarg = dynamic_cast<MargulesVPSSTP *> (tempIons->neutralMoleculePhase_);
1540 
1541 
1542  //m_thermo->getdlnActCoeffdlnX( DATA_PTR(grad_lnAC) );
1543  for (size_t k = 0; k < m_nDim; k++) {
1544  grad_T = m_Grad_T[k];
1545  grad_X.assign(m_Grad_X.begin()+m_nsp*k,m_Grad_X.begin()+m_nsp*(k+1));
1546  m_thermo->getdlnActCoeffds(grad_T, DATA_PTR(grad_X), DATA_PTR(grad_lnAC));
1547  for (size_t i = 0; i < m_nsp; i++)
1548  if (m_molefracs[i] < 1.e-15) {
1549  grad_lnAC[i] = 0;
1550  } else {
1551  grad_lnAC[i] += grad_X[i]/m_molefracs[i];
1552  }
1553  copy(grad_lnAC.begin(),grad_lnAC.end(),m_Grad_lnAC.begin()+m_nsp*k);
1554  // std::cout << k << " m_Grad_lnAC = " << m_Grad_lnAC[k] << std::endl;
1555  }
1556 
1557  return;
1558 }
1559 //====================================================================================================================
1560 /*
1561  *
1562  * Solve for the diffusional velocities in the Stefan-Maxwell equations
1563  *
1564  */
1565 // Solve the stefan_maxell equations for the diffusive fluxes.
1566 /*
1567  * The diffusive mass flux of species \e k is computed
1568  * using the Stefan-Maxwell equation
1569  * \f[
1570  * X_i \nabla \mu_i
1571  * = RT \sum_i \frac{X_i X_j}{D_{ij}}
1572  * ( \vec{V}_j - \vec{V}_i )
1573  * \f]
1574  * to determine the diffusion velocity and
1575  * \f[
1576  * \vec{N}_i = C_T X_i \vec{V}_i
1577  * \f]
1578  * to determine the diffusion flux. Here \f$ C_T \f$ is the
1579  * total concentration of the mixture [kmol/m^3], \f$ D_{ij} \f$
1580  * are the Stefa-Maxwell interaction parameters in [m^2/s],
1581  * \f$ \vec{V}_{i} \f$ is the diffusion velocity of species \e i,
1582  * \f$ \mu_i \f$ is the electrochemical potential of species \e i.
1583  *
1584  * The diffusion velocity is relative to an average velocity
1585  * that can be computed on a mole-weighted
1586  * or mass-weighted basis, or the diffusion velocities may
1587  * be specified as relative to a specific species (i.e. a
1588  * solvent) all according to the \verbatim <velocityBasis>
1589  * \endverbatim input parameter.
1590  *
1591  * One of the Stefan Maxwell equations is replaced by the appropriate
1592  * definition of the mass-averaged velocity, the mole-averaged velocity
1593  * or the specification that velocities are relative to that
1594  * of one species.
1595  */
1596 void LiquidTransport::stefan_maxwell_solve()
1597 {
1598  doublereal tmp;
1599  m_B.resize(m_nsp, m_nDim, 0.0);
1600  m_A.resize(m_nsp, m_nsp, 0.0);
1601 
1602  //! grab a local copy of the molecular weights
1603  const vector_fp& M = m_thermo->molecularWeights();
1604  //! grad a local copy of the ion molar volume (inverse total ion concentration)
1605  const doublereal vol = m_thermo->molarVolume();
1606 
1607  /*
1608  * Update the temperature, concentrations and diffusion coefficients in the mixture.
1609  */
1610  update_T();
1611  update_C();
1612  if (!m_diff_temp_ok) {
1613  updateDiff_T();
1614  }
1615 
1616  double T = m_thermo->temperature();
1617 
1618  update_Grad_lnAC() ;
1619 
1620  //m_thermo->getStandardVolumes(DATA_PTR(m_volume_spec));
1621  m_thermo->getActivityCoefficients(DATA_PTR(m_actCoeff));
1622 
1623  /*
1624  * Calculate the electrochemical potential gradient. This is the
1625  * driving force for relative diffusional transport.
1626  *
1627  * Here we calculate
1628  *
1629  * X_i * (grad (mu_i) + S_i grad T - M_i / dens * grad P
1630  *
1631  * This is Eqn. 13-1 p. 318 Newman. The original equation is from
1632  * Hershfeld, Curtis, and Bird.
1633  *
1634  * S_i is the partial molar entropy of species i. This term will cancel
1635  * out a lot of the grad T terms in grad (mu_i), therefore simplifying
1636  * the expression.
1637  *
1638  * Ok I think there may be many ways to do this. One way is to do it via basis
1639  * functions, at the nodes, as a function of the variables in the problem.
1640  *
1641  * For calculation of molality based thermo systems, we current get
1642  * the molar based values. This may change.
1643  *
1644  * Note, we have broken the symmetry of the matrix here, due to
1645  * consideratins involving species concentrations going to zero.
1646  *
1647  */
1648  for (size_t i = 0; i < m_nsp; i++) {
1649  for (size_t a = 0; a < m_nDim; a++) {
1650  m_Grad_mu[a*m_nsp + i] =
1651  m_chargeSpecies[i] * Faraday * m_Grad_V[a]
1652  //+ (m_volume_spec[i] - M[i]/dens_) * m_Grad_P[a]
1653  + GasConstant * T * m_Grad_lnAC[a*m_nsp+i];
1654  }
1655  }
1656 
1657  if (m_thermo->activityConvention() == cAC_CONVENTION_MOLALITY) {
1658  int iSolvent = 0;
1659  double mwSolvent = m_thermo->molecularWeight(iSolvent);
1660  double mnaught = mwSolvent/ 1000.;
1661  double lnmnaught = log(mnaught);
1662  for (size_t i = 1; i < m_nsp; i++) {
1663  for (size_t a = 0; a < m_nDim; a++) {
1664  m_Grad_mu[a*m_nsp + i] -=
1665  m_molefracs[i] * GasConstant * m_Grad_T[a] * lnmnaught;
1666  }
1667  }
1668  }
1669 
1670  /*
1671  * Just for Note, m_A(i,j) refers to the ith row and jth column.
1672  * They are still fortran ordered, so that i varies fastest.
1673  */
1674 
1675  double condSum1;
1676  switch (m_nDim) {
1677  case 1: /* 1-D approximation */
1678 
1679  m_B(0,0) = 0.0;
1680  //equation for the reference velocity
1681  for (size_t j = 0; j < m_nsp; j++) {
1682  if (m_velocityBasis == VB_MOLEAVG) {
1683  m_A(0,j) = m_molefracs_tran[j];
1684  } else if (m_velocityBasis == VB_MASSAVG) {
1685  m_A(0,j) = m_massfracs_tran[j];
1686  } else if ((m_velocityBasis >= 0)
1687  && (m_velocityBasis < static_cast<int>(m_nsp)))
1688  // use species number m_velocityBasis as reference velocity
1689  if (m_velocityBasis == static_cast<int>(j)) {
1690  m_A(0,j) = 1.0;
1691  } else {
1692  m_A(0,j) = 0.0;
1693  }
1694  else
1695  throw CanteraError("LiquidTransport::stefan_maxwell_solve",
1696  "Unknown reference velocity provided.");
1697  }
1698  for (size_t i = 1; i < m_nsp; i++) {
1699  m_B(i,0) = m_Grad_mu[i] / (GasConstant * T);
1700  m_A(i,i) = 0.0;
1701  for (size_t j = 0; j < m_nsp; j++) {
1702  if (j != i) {
1703  //if ( !( m_bdiff(i,j) > 0.0 ) )
1704  //throw CanteraError("LiquidTransport::stefan_maxwell_solve",
1705  // "m_bdiff has zero entry in non-diagonal.");
1706  tmp = m_molefracs_tran[j] * m_bdiff(i,j);
1707  m_A(i,i) -= tmp;
1708  m_A(i,j) = + tmp;
1709  }
1710  }
1711  }
1712 
1713  //! invert and solve the system Ax = b. Answer is in m_B
1714  solve(m_A, m_B);
1715 
1716  /*
1717  condSum2 = m_chargeSpecies[1]*m_chargeSpecies[1]*m_molefracs_tran[1]*m_bdiff(2,3) +
1718  m_chargeSpecies[2]*m_chargeSpecies[2]*m_molefracs_tran[2]*m_bdiff(1,3) +
1719  m_chargeSpecies[3]*m_chargeSpecies[3]*m_molefracs_tran[3]*m_bdiff(1,2);
1720  condSum1 = m_molefracs_tran[1]*m_bdiff(1,2)*m_bdiff(1,3) +
1721  m_molefracs_tran[2]*m_bdiff(2,3)*m_bdiff(1,2) +
1722  m_molefracs_tran[3]*m_bdiff(1,3)*m_bdiff(2,3);
1723  condSum2 = condSum2/condSum1*Faraday*Faraday/GasConstant/T/vol;
1724  */
1725 
1726  condSum1 = 0;
1727  for (size_t i = 0; i < m_nsp; i++) {
1728  condSum1 -= Faraday*m_chargeSpecies[i]*m_B(i,0)*m_molefracs_tran[i]/vol;
1729  }
1730 
1731  /*
1732  Check Mobility Ratio of Cations
1733  cout << "mobility ratio = " << m_chargeSpecies[1]*(m_B(1,0)-m_B(2,0))/m_chargeSpecies[0]/(m_B(0,0)-m_B(2,0)) << endl;
1734  */
1735 
1736  // cout << condSum1 << " = " << condSum2 << endl;
1737 
1738 
1739  break;
1740  case 2: /* 2-D approximation */
1741  m_B(0,0) = 0.0;
1742  m_B(0,1) = 0.0;
1743  //equation for the reference velocity
1744  for (size_t j = 0; j < m_nsp; j++) {
1745  if (m_velocityBasis == VB_MOLEAVG) {
1746  m_A(0,j) = m_molefracs_tran[j];
1747  } else if (m_velocityBasis == VB_MASSAVG) {
1748  m_A(0,j) = m_massfracs_tran[j];
1749  } else if ((m_velocityBasis >= 0)
1750  && (m_velocityBasis < static_cast<int>(m_nsp)))
1751  // use species number m_velocityBasis as reference velocity
1752  if (m_velocityBasis == static_cast<int>(j)) {
1753  m_A(0,j) = 1.0;
1754  } else {
1755  m_A(0,j) = 0.0;
1756  }
1757  else
1758  throw CanteraError("LiquidTransport::stefan_maxwell_solve",
1759  "Unknown reference velocity provided.");
1760  }
1761  for (size_t i = 1; i < m_nsp; i++) {
1762  m_B(i,0) = m_Grad_mu[i] / (GasConstant * T);
1763  m_B(i,1) = m_Grad_mu[m_nsp + i] / (GasConstant * T);
1764  m_A(i,i) = 0.0;
1765  for (size_t j = 0; j < m_nsp; j++) {
1766  if (j != i) {
1767  //if ( !( m_bdiff(i,j) > 0.0 ) )
1768  //throw CanteraError("LiquidTransport::stefan_maxwell_solve",
1769  // "m_bdiff has zero entry in non-diagonal.");
1770  tmp = m_molefracs_tran[j] * m_bdiff(i,j);
1771  m_A(i,i) -= tmp;
1772  m_A(i,j) = + tmp;
1773  }
1774  }
1775  }
1776 
1777  //! invert and solve the system Ax = b. Answer is in m_B
1778  solve(m_A, m_B);
1779 
1780 
1781  break;
1782 
1783  case 3: /* 3-D approximation */
1784  m_B(0,0) = 0.0;
1785  m_B(0,1) = 0.0;
1786  m_B(0,2) = 0.0;
1787  //equation for the reference velocity
1788  for (size_t j = 0; j < m_nsp; j++) {
1789  if (m_velocityBasis == VB_MOLEAVG) {
1790  m_A(0,j) = m_molefracs_tran[j];
1791  } else if (m_velocityBasis == VB_MASSAVG) {
1792  m_A(0,j) = m_massfracs_tran[j];
1793  } else if ((m_velocityBasis >= 0)
1794  && (m_velocityBasis < static_cast<int>(m_nsp)))
1795  // use species number m_velocityBasis as reference velocity
1796  if (m_velocityBasis == static_cast<int>(j)) {
1797  m_A(0,j) = 1.0;
1798  } else {
1799  m_A(0,j) = 0.0;
1800  }
1801  else
1802  throw CanteraError("LiquidTransport::stefan_maxwell_solve",
1803  "Unknown reference velocity provided.");
1804  }
1805  for (size_t i = 1; i < m_nsp; i++) {
1806  m_B(i,0) = m_Grad_mu[i] / (GasConstant * T);
1807  m_B(i,1) = m_Grad_mu[m_nsp + i] / (GasConstant * T);
1808  m_B(i,2) = m_Grad_mu[2*m_nsp + i] / (GasConstant * T);
1809  m_A(i,i) = 0.0;
1810  for (size_t j = 0; j < m_nsp; j++) {
1811  if (j != i) {
1812  //if ( !( m_bdiff(i,j) > 0.0 ) )
1813  //throw CanteraError("LiquidTransport::stefan_maxwell_solve",
1814  // "m_bdiff has zero entry in non-diagonal.");
1815  tmp = m_molefracs_tran[j] * m_bdiff(i,j);
1816  m_A(i,i) -= tmp;
1817  m_A(i,j) = + tmp;
1818  }
1819  }
1820  }
1821 
1822  //! invert and solve the system Ax = b. Answer is in m_B
1823  solve(m_A, m_B);
1824 
1825  break;
1826  default:
1827  printf("unimplemented\n");
1828  throw CanteraError("routine", "not done");
1829  break;
1830  }
1831 
1832  for (size_t a = 0; a < m_nDim; a++) {
1833  for (size_t j = 0; j < m_nsp; j++) {
1834  m_Vdiff(j,a) = m_B(j,a);
1835  m_flux(j,a) = concTot_ * M[j] * m_molefracs_tran[j] * m_B(j,a);
1836  }
1837  }
1838 }
1839 //====================================================================================================================
1840 // Throw an exception indicating something is not yet implemented.
1841 /*
1842  * @param msg String with an informative message
1843  */
1844 doublereal LiquidTransport::err(std::string msg) const
1845 {
1846  throw CanteraError("LiquidTransport::err()",
1847  "\n\n\n**** Method "+ msg +" not implemented in model "
1848  + int2str(model()) + " ****\n"
1849  "(Did you forget to specify a transport model?)\n\n\n");
1850  return 0.0;
1851 }
1852 //====================================================================================================================
1853 }
1854 //======================================================================================================================
Generated by   doxygen 1.8.2