Cantera  2.0
AqueousTransport.cpp
Go to the documentation of this file.
1 /**
2  * @file AqueousTransport.cpp
3  * Transport properties for aqueous systems
4  */
5 
8 
13 
15 
17 
18 #include <iostream>
19 #include <cstdio>
20 
21 using namespace std;
22 
23 /**
24  * Mole fractions below MIN_X will be set to MIN_X when computing
25  * transport properties.
26  */
27 #define MIN_X 1.e-20
28 
29 namespace Cantera
30 {
31 
32 
33 //====================================================================================================================
34 AqueousTransport::AqueousTransport() :
35  m_tmin(-1.0),
36  m_tmax(100000.),
37  m_iStateMF(-1),
38  m_temp(-1.0),
39  m_logt(0.0),
40  m_sqrt_t(-1.0),
41  m_t14(-1.0),
42  m_t32(-1.0),
43  m_sqrt_kbt(-1.0),
44  m_press(-1.0),
45  m_lambda(-1.0),
46  m_viscmix(-1.0),
47  m_viscmix_ok(false),
48  m_viscwt_ok(false),
49  m_spvisc_ok(false),
50  m_diffmix_ok(false),
51  m_bindiff_ok(false),
52  m_spcond_ok(false),
53  m_condmix_ok(false),
54  m_mode(-1000),
55  m_debug(false),
56  m_nDim(1)
57 {
58 }
59 
60 //====================================================================================================================
61 // Initialize the object
62 /*
63  * This is where we dimension everything.
64  */
66 {
67 
68  // constant substance attributes
69  m_thermo = tr.thermo;
70  m_nsp = m_thermo->nSpecies();
71  m_tmin = m_thermo->minTemp();
72  m_tmax = m_thermo->maxTemp();
73 
74  // make a local copy of the molecular weights
75  m_mw.resize(m_nsp);
76  copy(m_thermo->molecularWeights().begin(),
77  m_thermo->molecularWeights().end(), m_mw.begin());
78 
79  // copy polynomials and parameters into local storage
80  //m_visccoeffs = tr.visccoeffs;
81  //m_condcoeffs = tr.condcoeffs;
82  //m_diffcoeffs = tr.diffcoeffs;
83  cout << "In AqueousTransport::initLiquid we need to replace" << endl
84  << "LiquidTransportParams polynomial coefficients with" << endl
85  << "those in LiquidTransportData as in SimpleTransport." << endl;
86 
87  m_mode = tr.mode_;
88 
89  m_phi.resize(m_nsp, m_nsp, 0.0);
90 
91 
92  m_wratjk.resize(m_nsp, m_nsp, 0.0);
93  m_wratkj1.resize(m_nsp, m_nsp, 0.0);
94  for (size_t j = 0; j < m_nsp; j++)
95  for (size_t k = j; k < m_nsp; k++) {
96  m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
97  m_wratjk(k,j) = sqrt(m_wratjk(j,k));
98  m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
99  }
100 
101  m_polytempvec.resize(5);
102  m_visc.resize(m_nsp);
103  m_sqvisc.resize(m_nsp);
104  m_cond.resize(m_nsp);
105  m_bdiff.resize(m_nsp, m_nsp);
106 
107  m_molefracs.resize(m_nsp);
108  m_spwork.resize(m_nsp);
109 
110  // resize the internal gradient variables
111  m_Grad_X.resize(m_nDim * m_nsp, 0.0);
112  m_Grad_T.resize(m_nDim, 0.0);
113  m_Grad_V.resize(m_nDim, 0.0);
114  m_Grad_mu.resize(m_nDim * m_nsp, 0.0);
115 
116 
117  // set all flags to false
118  m_viscmix_ok = false;
119  m_viscwt_ok = false;
120  m_spvisc_ok = false;
121  m_spcond_ok = false;
122  m_condmix_ok = false;
123  m_spcond_ok = false;
124  m_diffmix_ok = false;
125 
126  return true;
127 }
128 //====================================================================================================================
129 /*
130  * The viscosity is computed using the Wilke mixture rule.
131  * \f[
132  * \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}.
133  * \f]
134  * Here \f$ \mu_k \f$ is the viscosity of pure species \e k,
135  * and
136  * \f[
137  * \Phi_{k,j} = \frac{\left[1
138  * + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2}
139  * {\sqrt{8}\sqrt{1 + M_k/M_j}}
140  * \f]
141  * @see updateViscosity_T();
142  */
144 {
145 
146  update_T();
147  update_C();
148 
149  if (m_viscmix_ok) {
150  return m_viscmix;
151  }
152 
153  // update m_visc[] and m_phi[] if necessary
154  if (!m_viscwt_ok) {
156  }
157 
159 
160  m_viscmix = 0.0;
161  for (size_t k = 0; k < m_nsp; k++) {
162  m_viscmix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
163  }
164  return m_viscmix;
165 }
166 //====================================================================================================================
167 // Returns the pure species viscosities
168 /*
169  *
170  * Controlling update boolean = m_viscwt_ok
171  *
172  * @param visc Vector of species viscosities
173  */
174 void AqueousTransport::getSpeciesViscosities(doublereal* const visc)
175 {
177  copy(m_visc.begin(), m_visc.end(), visc);
178 }
179 //====================================================================================================================
180 void AqueousTransport::getBinaryDiffCoeffs(const size_t ld, doublereal* const d)
181 {
182  update_T();
183 
184  // if necessary, evaluate the binary diffusion coefficients
185  // from the polynomial fits
186  if (!m_bindiff_ok) {
187  updateDiff_T();
188  }
189  doublereal pres = m_thermo->pressure();
190 
191  doublereal rp = 1.0/pres;
192  for (size_t i = 0; i < m_nsp; i++)
193  for (size_t j = 0; j < m_nsp; j++) {
194  d[ld*j + i] = rp * m_bdiff(i,j);
195  }
196 }
197 //====================================================================================================================
198 // Get the electrical Mobilities (m^2/V/s).
199 /*
200  * This function returns the mobilities. In some formulations
201  * this is equal to the normal mobility multiplied by faraday's constant.
202  *
203  * Frequently, but not always, the mobility is calculated from the
204  * diffusion coefficient using the Einstein relation
205  *
206  * \f[
207  * \mu^e_k = \frac{F D_k}{R T}
208  * \f]
209  *
210  * @param mobil_e Returns the mobilities of
211  * the species in array \c mobil_e. The array must be
212  * dimensioned at least as large as the number of species.
213  */
214 void AqueousTransport::getMobilities(doublereal* const mobil)
215 {
217  doublereal c1 = ElectronCharge / (Boltzmann * m_temp);
218  for (size_t k = 0; k < m_nsp; k++) {
219  mobil[k] = c1 * m_spwork[k];
220  }
221 }
222 //====================================================================================================================
223 void AqueousTransport::getFluidMobilities(doublereal* const mobil)
224 {
226  doublereal c1 = 1.0 / (GasConstant * m_temp);
227  for (size_t k = 0; k < m_nsp; k++) {
228  mobil[k] = c1 * m_spwork[k];
229  }
230 }
231 //====================================================================================================================
232 void AqueousTransport::set_Grad_V(const doublereal* const grad_V)
233 {
234  for (size_t a = 0; a < m_nDim; a++) {
235  m_Grad_V[a] = grad_V[a];
236  }
237 }
238 //====================================================================================================================
239 void AqueousTransport::set_Grad_T(const doublereal* const grad_T)
240 {
241  for (size_t a = 0; a < m_nDim; a++) {
242  m_Grad_T[a] = grad_T[a];
243  }
244 }
245 //====================================================================================================================
246 void AqueousTransport::set_Grad_X(const doublereal* const grad_X)
247 {
248  size_t itop = m_nDim * m_nsp;
249  for (size_t i = 0; i < itop; i++) {
250  m_Grad_X[i] = grad_X[i];
251  }
252 }
253 //====================================================================================================================
254 /*
255  * The thermal conductivity is computed from the following mixture rule:
256  * \[
257  * \lambda = 0.5 \left( \sum_k X_k \lambda_k
258  * + \frac{1}{\sum_k X_k/\lambda_k}\right)
259  * \]
260  */
262 {
263  update_T();
264  update_C();
265 
266  if (!m_spcond_ok) {
267  updateCond_T();
268  }
269  if (!m_condmix_ok) {
270  doublereal sum1 = 0.0, sum2 = 0.0;
271  for (size_t k = 0; k < m_nsp; k++) {
272  sum1 += m_molefracs[k] * m_cond[k];
273  sum2 += m_molefracs[k] / m_cond[k];
274  }
275  m_lambda = 0.5*(sum1 + 1.0/sum2);
276  }
277  return m_lambda;
278 }
279 //====================================================================================================================
280 // Return a vector of Thermal diffusion coefficients [kg/m/sec].
281 /*
282  * The thermal diffusion coefficient \f$ D^T_k \f$ is defined
283  * so that the diffusive mass flux of species <I>k<\I> induced by the
284  * local temperature gradient is given by the following formula
285  *
286  * \f[
287  * M_k J_k = -D^T_k \nabla \ln T.
288  * \f]
289  *
290  * The thermal diffusion coefficient can be either positive or negative.
291  *
292  * In this method we set it to zero.
293  *
294  * @param dt On return, dt will contain the species thermal
295  * diffusion coefficients. Dimension dt at least as large as
296  * the number of species. Units are kg/m/s.
297  */
298 void AqueousTransport::getThermalDiffCoeffs(doublereal* const dt)
299 {
300  for (size_t k = 0; k < m_nsp; k++) {
301  dt[k] = 0.0;
302  }
303 }
304 //====================================================================================================================
305 // Get the species diffusive mass fluxes wrt to the specified solution averaged velocity,
306 // given the gradients in mole fraction and temperature
307 /*
308  * Units for the returned fluxes are kg m-2 s-1.
309  *
310  * Usually the specified solution average velocity is the mass averaged velocity.
311  * This is changed in some subclasses, however.
312  *
313  * @param ndim Number of dimensions in the flux expressions
314  * @param grad_T Gradient of the temperature
315  * (length = ndim)
316  * @param ldx Leading dimension of the grad_X array
317  * (usually equal to m_nsp but not always)
318  * @param grad_X Gradients of the mole fraction
319  * Flat vector with the m_nsp in the inner loop.
320  * length = ldx * ndim
321  * @param ldf Leading dimension of the fluxes array
322  * (usually equal to m_nsp but not always)
323  * @param fluxes Output of the diffusive mass fluxes
324  * Flat vector with the m_nsp in the inner loop.
325  * length = ldx * ndim
326  */
327 void AqueousTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_T,
328  size_t ldx, const doublereal* const grad_X,
329  size_t ldf, doublereal* const fluxes)
330 {
331  set_Grad_T(grad_T);
332  set_Grad_X(grad_X);
333  getSpeciesFluxesExt(ldf, fluxes);
334 }
335 //====================================================================================================================
336 // Return the species diffusive mass fluxes wrt to the specified averaged velocity,
337 /*
338  * This method acts similarly to getSpeciesFluxesES() but
339  * requires all gradients to be preset using methods set_Grad_X(), set_Grad_V(), set_Grad_T().
340  * See the documentation of getSpeciesFluxesES() for details.
341  *
342  * units = kg/m2/s
343  *
344  * Internally, gradients in the in mole fraction, temperature
345  * and electrostatic potential contribute to the diffusive flux
346  *
347  * The diffusive mass flux of species \e k is computed from the following formula
348  *
349  * \f[
350  * j_k = - \rho M_k D_k \nabla X_k - Y_k V_c
351  * \f]
352  *
353  * where V_c is the correction velocity
354  *
355  * \f[
356  * V_c = - \sum_j {\rho M_j D_j \nabla X_j}
357  * \f]
358  *
359  * @param ldf Stride of the fluxes array. Must be equal to or greater than the number of species.
360  * @param fluxes Output of the diffusive fluxes. Flat vector with the m_nsp in the inner loop.
361  * length = ldx * ndim
362  */
363 void AqueousTransport::getSpeciesFluxesExt(size_t ldf, doublereal* const fluxes)
364 {
365  update_T();
366  update_C();
367 
369 
370  const vector_fp& mw = m_thermo->molecularWeights();
371  const doublereal* y = m_thermo->massFractions();
372  doublereal rhon = m_thermo->molarDensity();
373  // Unroll wrt ndim
374  vector_fp sum(m_nDim,0.0);
375  for (size_t n = 0; n < m_nDim; n++) {
376  for (size_t k = 0; k < m_nsp; k++) {
377  fluxes[n*ldf + k] = -rhon * mw[k] * m_spwork[k] * m_Grad_X[n*m_nsp + k];
378  sum[n] += fluxes[n*ldf + k];
379  }
380  }
381  // add correction flux to enforce sum to zero
382  for (size_t n = 0; n < m_nDim; n++) {
383  for (size_t k = 0; k < m_nsp; k++) {
384  fluxes[n*ldf + k] -= y[k]*sum[n];
385  }
386  }
387 }
388 //====================================================================================================================
389 /**
390  * Mixture-averaged diffusion coefficients [m^2/s].
391  *
392  * For the single species case or the pure fluid case
393  * the routine returns the self-diffusion coefficient.
394  * This is need to avoid a Nan result in the formula
395  * below.
396  */
397 void AqueousTransport::getMixDiffCoeffs(doublereal* const d)
398 {
399 
400  update_T();
401  update_C();
402 
403  // update the binary diffusion coefficients if necessary
404  if (!m_bindiff_ok) {
405  updateDiff_T();
406  }
407 
408  size_t k, j;
409  doublereal mmw = m_thermo->meanMolecularWeight();
410  doublereal sumxw = 0.0, sum2;
411  doublereal p = m_press;
412  if (m_nsp == 1) {
413  d[0] = m_bdiff(0,0) / p;
414  } else {
415  for (k = 0; k < m_nsp; k++) {
416  sumxw += m_molefracs[k] * m_mw[k];
417  }
418  for (k = 0; k < m_nsp; k++) {
419  sum2 = 0.0;
420  for (j = 0; j < m_nsp; j++) {
421  if (j != k) {
422  sum2 += m_molefracs[j] / m_bdiff(j,k);
423  }
424  }
425  if (sum2 <= 0.0) {
426  d[k] = m_bdiff(k,k) / p;
427  } else {
428  d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2);
429  }
430  }
431  }
432 }
433 
434 //====================================================================================================================
435 // Handles the effects of changes in the Temperature, internally
436 // within the object.
437 /*
438  * This is called whenever a transport property is
439  * requested.
440  * The first task is to check whether the temperature has changed
441  * since the last call to update_T().
442  * If it hasn't then an immediate return is carried out.
443  *
444  * @internal
445  */
447 {
448  doublereal t = m_thermo->temperature();
449  if (t == m_temp) {
450  return;
451  }
452  if (t < 0.0) {
453  throw CanteraError("AqueousTransport::update_T",
454  "negative temperature "+fp2str(t));
455  }
456 
457  // Compute various functions of temperature
458  m_temp = t;
459  m_logt = log(m_temp);
460  m_kbt = Boltzmann * m_temp;
461  m_sqrt_t = sqrt(m_temp);
462  m_t14 = sqrt(m_sqrt_t);
463  m_t32 = m_temp * m_sqrt_t;
464  m_sqrt_kbt = sqrt(Boltzmann*m_temp);
465 
466  // compute powers of log(T)
467  m_polytempvec[0] = 1.0;
468  m_polytempvec[1] = m_logt;
470  m_polytempvec[3] = m_logt*m_logt*m_logt;
471  m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt;
472 
473  // temperature has changed, so polynomial temperature
474  // interpolations will need to be reevaluated.
475  // Set all of these flags to false
476  m_viscmix_ok = false;
477  m_spvisc_ok = false;
478  m_viscwt_ok = false;
479  m_spcond_ok = false;
480  m_diffmix_ok = false;
481  m_bindiff_ok = false;
482  m_condmix_ok = false;
483 
484  // For now, for a concentration redo also
485  m_iStateMF = -1;
486 }
487 //====================================================================================================================
488 /**
489  * @internal This is called the first time any transport property
490  * is requested from Mixture after the concentrations
491  * have changed.
492  */
494 {
495 
496  doublereal pres = m_thermo->pressure();
497  // Check for changes in the mole fraction vector.
498  //int iStateNew = m_thermo->getIStateMF();
499  //if (iStateNew == m_iStateMF) {
500  // if (pres == m_press) {
501  // return;
502  // }
503  // } else {
504  // m_iStateMF = iStateNew;
505  //}
506  m_press = pres;
507 
508  // signal that concentration-dependent quantities will need to
509  // be recomputed before use, and update the local mole
510  // fractions.
511 
512  m_viscmix_ok = false;
513  m_diffmix_ok = false;
514  m_condmix_ok = false;
515 
517 
518  // add an offset to avoid a pure species condition or
519  // negative mole fractions. MIN_X is 1.0E-20, a value
520  // which is below the additive machine precision of mole fractions.
521  for (size_t k = 0; k < m_nsp; k++) {
523  }
524 }
525 //====================================================================================================================
526 /*
527  * Update the temperature-dependent parts of the mixture-averaged
528  * thermal conductivity.
529  */
531 {
532  if (m_mode == CK_Mode) {
533  for (size_t k = 0; k < m_nsp; k++) {
534  m_cond[k] = exp(dot4(m_polytempvec, m_condcoeffs[k]));
535  }
536  } else {
537  for (size_t k = 0; k < m_nsp; k++) {
539  }
540  }
541  m_spcond_ok = true;
542  m_condmix_ok = false;
543 }
544 //====================================================================================================================
545 /*
546  * Update the binary diffusion coefficients. These are evaluated
547  * from the polynomial fits at unit pressure (1 Pa).
548  */
550 {
551 
552  // evaluate binary diffusion coefficients at unit pressure
553  size_t ic = 0;
554  if (m_mode == CK_Mode) {
555  for (size_t i = 0; i < m_nsp; i++) {
556  for (size_t j = i; j < m_nsp; j++) {
557  m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic]));
558  m_bdiff(j,i) = m_bdiff(i,j);
559  ic++;
560  }
561  }
562  } else {
563  for (size_t i = 0; i < m_nsp; i++) {
564  for (size_t j = i; j < m_nsp; j++) {
566  m_diffcoeffs[ic]);
567  m_bdiff(j,i) = m_bdiff(i,j);
568  ic++;
569  }
570  }
571  }
572 
573  m_bindiff_ok = true;
574  m_diffmix_ok = false;
575 }
576 //====================================================================================================================
577 /*
578  * Update the pure-species viscosities.
579  */
581 {
582  if (m_mode == CK_Mode) {
583  for (size_t k = 0; k < m_nsp; k++) {
584  m_visc[k] = exp(dot4(m_polytempvec, m_visccoeffs[k]));
585  m_sqvisc[k] = sqrt(m_visc[k]);
586  }
587  } else {
588  for (size_t k = 0; k < m_nsp; k++) {
589  // the polynomial fit is done for sqrt(visc/sqrt(T))
591  m_visc[k] = (m_sqvisc[k]*m_sqvisc[k]);
592  }
593  }
594  m_spvisc_ok = true;
595 }
596 //====================================================================================================================
597 /*
598  * Update the temperature-dependent viscosity terms.
599  * Updates the array of pure species viscosities, and the
600  * weighting functions in the viscosity mixture rule.
601  * The flag m_visc_ok is set to true.
602  */
604 {
605  doublereal vratiokj, wratiojk, factor1;
606 
607  if (!m_spvisc_ok) {
609  }
610 
611  // see Eq. (9-5.15) of Reid, Prausnitz, and Poling
612  for (size_t j = 0; j < m_nsp; j++) {
613  for (size_t k = j; k < m_nsp; k++) {
614  vratiokj = m_visc[k]/m_visc[j];
615  wratiojk = m_mw[j]/m_mw[k];
616 
617  // Note that m_wratjk(k,j) holds the square root of
618  // m_wratjk(j,k)!
619  factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
620  m_phi(k,j) = factor1*factor1 /
621  (SqrtEight * m_wratkj1(j,k));
622  m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
623  }
624  }
625  m_viscwt_ok = true;
626 }
627 //====================================================================================================================
628 /*
629  * This function returns a Transport data object for a given species.
630  *
631  */
633 {
635  td.speciesName = m_thermo->speciesName(kSpecies);
636 
637 
638  return td;
639 }
640 //====================================================================================================================
641 /*
642  *
643  * Solve for the diffusional velocities in the Stefan-Maxwell equations
644  *
645  */
647 {
648  size_t VIM = 2;
649  m_B.resize(m_nsp, VIM);
650  // grab a local copy of the molecular weights
651  const vector_fp& M = m_thermo->molecularWeights();
652 
653 
654  // get the mean molecular weight of the mixture
655  //double M_mix = m_thermo->meanMolecularWeight();
656 
657 
658  // get the concentration of the mixture
659  //double rho = m_thermo->density();
660  //double c = rho/M_mix;
661 
662 
664 
665  double T = m_thermo->temperature();
666 
667 
668  /* electrochemical potential gradient */
669  for (size_t i = 0; i < m_nsp; i++) {
670  for (size_t a = 0; a < VIM; a++) {
671  m_Grad_mu[a*m_nsp + i] = m_chargeSpecies[i] * Faraday * m_Grad_V[a]
672  + (GasConstant*T/m_molefracs[i]) * m_Grad_X[a*m_nsp+i];
673  }
674  }
675 
676  /*
677  * Just for Note, m_A(i,j) refers to the ith row and jth column.
678  * They are still fortran ordered, so that i varies fastest.
679  */
680  switch (VIM) {
681  case 1: /* 1-D approximation */
682  m_B(0,0) = 0.0;
683  for (size_t j = 0; j < m_nsp; j++) {
684  m_A(0,j) = 1.0;
685  }
686  for (size_t i = 1; i < m_nsp; i++) {
687  m_B(i,0) = m_concentrations[i] * m_Grad_mu[i] / (GasConstant * T);
688  for (size_t j = 0; j < m_nsp; j++) {
689  if (j != i) {
690  m_A(i,j) = m_molefracs[i] / (M[j] * m_DiffCoeff_StefMax(i,j));
691  m_A(i,i) -= m_molefracs[j] / (M[i] * m_DiffCoeff_StefMax(i,j));
692  } else if (j == i) {
693  m_A(i,i) = 0.0;
694  }
695  }
696  }
697 
698  //! invert and solve the system Ax = b. Answer is in m_B
699  solve(m_A, m_B.ptrColumn(0));
700 
701  m_flux = m_B;
702 
703 
704  break;
705  case 2: /* 2-D approximation */
706  m_B(0,0) = 0.0;
707  m_B(0,1) = 0.0;
708  for (size_t j = 0; j < m_nsp; j++) {
709  m_A(0,j) = 1.0;
710  }
711  for (size_t i = 1; i < m_nsp; i++) {
712  m_B(i,0) = m_concentrations[i] * m_Grad_mu[i] / (GasConstant * T);
713  m_B(i,1) = m_concentrations[i] * m_Grad_mu[m_nsp + i] / (GasConstant * T);
714  for (size_t j = 0; j < m_nsp; j++) {
715  if (j != i) {
716  m_A(i,j) = m_molefracs[i] / (M[j] * m_DiffCoeff_StefMax(i,j));
717  m_A(i,i) -= m_molefracs[j] / (M[i] * m_DiffCoeff_StefMax(i,j));
718  } else if (j == i) {
719  m_A(i,i) = 0.0;
720  }
721  }
722  }
723 
724  //! invert and solve the system Ax = b. Answer is in m_B
725  //solve(m_A, m_B);
726 
727  m_flux = m_B;
728 
729 
730  break;
731 
732  case 3: /* 3-D approximation */
733  m_B(0,0) = 0.0;
734  m_B(0,1) = 0.0;
735  m_B(0,2) = 0.0;
736  for (size_t j = 0; j < m_nsp; j++) {
737  m_A(0,j) = 1.0;
738  }
739  for (size_t i = 1; i < m_nsp; i++) {
740  m_B(i,0) = m_concentrations[i] * m_Grad_mu[i] / (GasConstant * T);
741  m_B(i,1) = m_concentrations[i] * m_Grad_mu[m_nsp + i] / (GasConstant * T);
742  m_B(i,2) = m_concentrations[i] * m_Grad_mu[2*m_nsp + i] / (GasConstant * T);
743  for (size_t j = 0; j < m_nsp; j++) {
744  if (j != i) {
745  m_A(i,j) = m_molefracs[i] / (M[j] * m_DiffCoeff_StefMax(i,j));
746  m_A(i,i) -= m_molefracs[j] / (M[i] * m_DiffCoeff_StefMax(i,j));
747  } else if (j == i) {
748  m_A(i,i) = 0.0;
749  }
750  }
751  }
752 
753  //! invert and solve the system Ax = b. Answer is in m_B
754  //solve(m_A, m_B);
755 
756  m_flux = m_B;
757 
758 
759  break;
760  default:
761  printf("unimplemented\n");
762  throw CanteraError("routine", "not done");
763  break;
764  }
765 
766 
767 }
768 //====================================================================================================================
769 }
770 //======================================================================================================================