Cantera  2.0
PecosTransport.cpp
Go to the documentation of this file.
1 /**
2  * @file PecosTransport.cpp
3  * Mixture-averaged transport properties.
4  *
5  */
6 
7 #include "cantera/thermo/ThermoPhase.h"
8 #include "cantera/transport/PecosTransport.h"
9 
10 #include "cantera/base/utilities.h"
11 #include "cantera/transport/TransportParams.h"
12 #include "cantera/transport/TransportFactory.h"
13 #include "cantera/base/stringUtils.h"
14 
15 #include "cantera/thermo/IdealGasPhase.h"
16 
17 #include <iostream>
18 using namespace std;
19 
20 namespace Cantera
21 {
22 
23 //////////////////// class PecosTransport methods //////////////
24 
25 PecosTransport::PecosTransport() :
26  m_nsp(0),
27  m_temp(-1.0),
28  m_logt(0.0)
29 {
30 
31 
32 }
33 
34 bool PecosTransport::initGas(GasTransportParams& tr)
35 {
36 
37  // constant substance attributes
38  m_thermo = tr.thermo;
39  m_nsp = m_thermo->nSpecies();
40 
41  // make a local copy of the molecular weights
42  m_mw.resize(m_nsp);
43  copy(m_thermo->molecularWeights().begin(),
44  m_thermo->molecularWeights().end(), m_mw.begin());
45 
46  // copy polynomials and parameters into local storage
47  m_poly = tr.poly;
48  m_visccoeffs = tr.visccoeffs;
49  m_condcoeffs = tr.condcoeffs;
50  m_diffcoeffs = tr.diffcoeffs;
51 
52  m_zrot = tr.zrot;
53  m_crot = tr.crot;
54  m_epsilon = tr.epsilon;
55  m_mode = tr.mode_;
56  m_diam = tr.diam;
57  m_eps = tr.eps;
58  m_alpha = tr.alpha;
59  m_dipoleDiag.resize(m_nsp);
60  for (int i = 0; i < m_nsp; i++) {
61  m_dipoleDiag[i] = tr.dipole(i,i);
62  }
63 
64  m_phi.resize(m_nsp, m_nsp, 0.0);
65  m_wratjk.resize(m_nsp, m_nsp, 0.0);
66  m_wratkj1.resize(m_nsp, m_nsp, 0.0);
67  int j, k;
68  for (j = 0; j < m_nsp; j++)
69  for (k = j; k < m_nsp; k++) {
70  m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]);
71  m_wratjk(k,j) = sqrt(m_wratjk(j,k));
72  m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]);
73  }
74 
75  m_polytempvec.resize(5);
76  m_visc.resize(m_nsp);
77  m_sqvisc.resize(m_nsp);
78  m_cond.resize(m_nsp);
79  m_bdiff.resize(m_nsp, m_nsp);
80 
81  m_molefracs.resize(m_nsp);
82  m_spwork.resize(m_nsp);
83 
84  // set flags all false
85  m_viscmix_ok = false;
86  m_viscwt_ok = false;
87  m_spvisc_ok = false;
88  m_spcond_ok = false;
89  m_condmix_ok = false;
90  m_spcond_ok = false;
91  m_diffmix_ok = false;
92  m_abc_ok = false;
93 
94  // read blottner fit parameters (A,B,C)
95  read_blottner_transport_table();
96 
97  // set specific heats
98  cv_rot.resize(m_nsp);
99  cp_R.resize(m_nsp);
100  cv_int.resize(m_nsp);
101 
102  for (k = 0; k < m_nsp; k++) {
103  cv_rot[k] = tr.crot[k];
104  cp_R[k] = ((IdealGasPhase*)tr.thermo)->cp_R_ref()[k];
105  cv_int[k] = cp_R[k] - 2.5 - cv_rot[k];
106  }
107  return true;
108 }
109 
110 
111 /*********************************************************
112  *
113  * Public methods
114  *
115  *********************************************************/
116 
117 
118 /****************** viscosity ******************************/
119 
120 /**
121  * The viscosity is computed using the Wilke mixture rule.
122  * \f[
123  * \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}.
124  * \f]
125  * Here \f$ \mu_k \f$ is the viscosity of pure species \e k,
126  * and
127  * \f[
128  * \Phi_{k,j} = \frac{\left[1
129  * + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2}
130  * {\sqrt{8}\sqrt{1 + M_k/M_j}}
131  * \f]
132  * @see updateViscosity_T();
133  */
134 doublereal PecosTransport::viscosity()
135 {
136 
137  update_T();
138  update_C();
139 
140  if (m_viscmix_ok) {
141  return m_viscmix;
142  }
143 
144  doublereal vismix = 0.0;
145  int k;
146  // update m_visc and m_phi if necessary
147  if (!m_viscwt_ok) {
148  updateViscosity_T();
149  }
150 
151  multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork));
152 
153  for (k = 0; k < m_nsp; k++) {
154  vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom;
155  }
156  m_viscmix = vismix;
157  return vismix;
158 }
159 
160 /******************* binary diffusion coefficients **************/
161 /*
162  *
163  * Using Ramshaw's self-consistent Effective Binary Diffusion
164  * (1990, J. Non-Equilib. Thermo)
165  * Adding more doxygen would be good here
166  */
167 
168 void PecosTransport::getBinaryDiffCoeffs(const size_t ld, doublereal* const d)
169 {
170  int i,j;
171 
172  update_T();
173 
174  // if necessary, evaluate the binary diffusion coefficents
175  if (!m_bindiff_ok) {
176  updateDiff_T();
177  }
178 
179  doublereal rp = 1.0/pressure_ig();
180  for (i = 0; i < m_nsp; i++)
181  for (j = 0; j < m_nsp; j++) {
182  d[ld*j + i] = rp * m_bdiff(i,j);
183  }
184 }
185 
186 
187 void PecosTransport::getMobilities(doublereal* const mobil)
188 {
189  int k;
190  getMixDiffCoeffs(DATA_PTR(m_spwork));
191  doublereal c1 = ElectronCharge / (Boltzmann * m_temp);
192  for (k = 0; k < m_nsp; k++) {
193  mobil[k] = c1 * m_spwork[k] * m_thermo->charge(k);
194  }
195 }
196 
197 
198 /****************** thermal conductivity **********************/
199 
200 /**
201  * The thermal conductivity is computed using the Wilke mixture rule.
202  * \f[
203  * \k = \sum_s \frac{k_s X_s}{\sum_j \Phi_{s,j} X_j}.
204  * \f]
205  * Here \f$ \k_s \f$ is the conductivity of pure species \e s,
206  * and
207  * \f[
208  * \Phi_{s,j} = \frac{\left[1
209  * + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_s}}\right)}\right]^2}
210  * {\sqrt{8}\sqrt{1 + M_s/M_j}}
211  * \f]
212  * @see updateCond_T();
213  */
214 doublereal PecosTransport::thermalConductivity()
215 {
216  int k;
217  doublereal lambda = 0.0;
218 
219  update_T();
220  update_C();
221 
222  // update m_cond and m_phi if necessary
223  if (!m_spcond_ok) {
224  updateCond_T();
225  }
226  if (!m_condmix_ok) {
227 
228  multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork));
229 
230  for (k = 0; k < m_nsp; k++) {
231  lambda += m_molefracs[k] * m_cond[k]/m_spwork[k]; //denom;
232  }
233 
234  }
235  m_lambda = lambda;
236  return m_lambda;
237 
238 }
239 
240 
241 /****************** thermal diffusion coefficients ************/
242 
243 /**
244  * Thermal diffusion is not considered in this pecos
245  * model. To include thermal diffusion, use transport manager
246  * MultiTransport instead. This methods fills out array dt with
247  * zeros.
248  */
249 void PecosTransport::getThermalDiffCoeffs(doublereal* const dt)
250 {
251  int k;
252  for (k = 0; k < m_nsp; k++) {
253  dt[k] = 0.0;
254  }
255 }
256 
257 /**
258  * @param ndim The number of spatial dimensions (1, 2, or 3).
259  * @param grad_T The temperature gradient (ignored in this model).
260  * @param ldx Leading dimension of the grad_X array.
261  * The diffusive mass flux of species \e k is computed from
262  * \f[
263  * \vec{j}_k = -n M_k D_k \nabla X_k + \frac{\rho_k}{\rho} \sum_r n M_r D_r \nabla X_r
264  * \f]
265  *
266  * This is neglective pressure, forced and thermal diffusion.
267  *
268  */
269 void PecosTransport::getSpeciesFluxes(size_t ndim,
270  const doublereal* const grad_T,
271  size_t ldx, const doublereal* const grad_X,
272  size_t ldf, doublereal* const fluxes)
273 {
274  size_t n = 0;
275  int k;
276 
277  update_T();
278  update_C();
279 
280  getMixDiffCoeffs(DATA_PTR(m_spwork));
281 
282  const vector_fp& mw = m_thermo->molecularWeights();
283  const doublereal* y = m_thermo->massFractions();
284  doublereal rhon = m_thermo->molarDensity();
285 
286  vector_fp sum(ndim,0.0);
287 
288  doublereal correction=0.0;
289  // grab 2nd (summation) term -- still need to multiply by mass fraction (\rho_s / \rho)
290  for (k = 0; k < m_nsp; k++) {
291  correction += rhon * mw[k] * m_spwork[k] * grad_X[n*ldx + k];
292  }
293 
294  for (n = 0; n < ndim; n++) {
295  for (k = 0; k < m_nsp; k++) {
296  fluxes[n*ldf + k] = -rhon * mw[k] * m_spwork[k] * grad_X[n*ldx + k] + y[k]*correction;
297  sum[n] += fluxes[n*ldf + k];
298  }
299  }
300  // add correction flux to enforce sum to zero
301  for (n = 0; n < ndim; n++) {
302  for (k = 0; k < m_nsp; k++) {
303  fluxes[n*ldf + k] -= y[k]*sum[n];
304  }
305  }
306 }
307 
308 /**
309  * Mixture-averaged diffusion coefficients [m^2/s].
310  *
311  * For the single species case or the pure fluid case
312  * the routine returns the self-diffusion coefficient.
313  * This is need to avoid a Nan result in the formula
314  * below.
315  */
316 void PecosTransport::getMixDiffCoeffs(doublereal* const d)
317 {
318 
319  update_T();
320  update_C();
321 
322  // update the binary diffusion coefficients if necessary
323  if (!m_bindiff_ok) {
324  updateDiff_T();
325  }
326 
327  int k, j;
328  doublereal mmw = m_thermo->meanMolecularWeight();
329  doublereal sumxw = 0.0, sum2;
330  doublereal p = pressure_ig();
331  if (m_nsp == 1) {
332  d[0] = m_bdiff(0,0) / p;
333  } else {
334  for (k = 0; k < m_nsp; k++) {
335  sumxw += m_molefracs[k] * m_mw[k];
336  }
337  for (k = 0; k < m_nsp; k++) {
338  sum2 = 0.0;
339  for (j = 0; j < m_nsp; j++) {
340  if (j != k) {
341  sum2 += m_molefracs[j] / m_bdiff(j,k);
342  }
343  }
344  if (sum2 <= 0.0) {
345  d[k] = m_bdiff(k,k) / p;
346  } else {
347  d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2);
348  }
349  }
350  }
351 }
352 
353 void PecosTransport::getMixDiffCoeffsMole(doublereal* const d)
354 {
355  update_T();
356  update_C();
357 
358  // update the binary diffusion coefficients if necessary
359  if (!m_bindiff_ok) {
360  updateDiff_T();
361  }
362 
363  doublereal p = m_thermo->pressure();
364  if (m_nsp == 1) {
365  d[0] = m_bdiff(0,0) / p;
366  } else {
367  for (int k = 0; k < m_nsp; k++) {
368  double sum2 = 0.0;
369  for (int j = 0; j < m_nsp; j++) {
370  if (j != k) {
371  sum2 += m_molefracs[j] / m_bdiff(j,k);
372  }
373  }
374  if (sum2 <= 0.0) {
375  d[k] = m_bdiff(k,k) / p;
376  } else {
377  d[k] = (1 - m_molefracs[k]) / (p * sum2);
378  }
379  }
380  }
381 }
382 
383 void PecosTransport::getMixDiffCoeffsMass(doublereal* const d)
384 {
385  update_T();
386  update_C();
387 
388  // update the binary diffusion coefficients if necessary
389  if (!m_bindiff_ok) {
390  updateDiff_T();
391  }
392 
393  doublereal mmw = m_thermo->meanMolecularWeight();
394  doublereal p = m_thermo->pressure();
395 
396  if (m_nsp == 1) {
397  d[0] = m_bdiff(0,0) / p;
398  } else {
399  for (int k=0; k<m_nsp; k++) {
400  double sum1 = 0.0;
401  double sum2 = 0.0;
402  for (int i=0; i<m_nsp; i++) {
403  if (i==k) {
404  continue;
405  }
406  sum1 += m_molefracs[i] / m_bdiff(k,i);
407  sum2 += m_molefracs[i] * m_mw[i] / m_bdiff(k,i);
408  }
409  sum1 *= p;
410  sum2 *= p * m_molefracs[k] / (mmw - m_mw[k]*m_molefracs[k]);
411  d[k] = 1.0 / (sum1 + sum2);
412  }
413  }
414 }
415 
416 /**
417  * @internal This is called whenever a transport property is
418  * requested from ThermoSubstance if the temperature has changed
419  * since the last call to update_T.
420  */
421 void PecosTransport::update_T()
422 {
423  doublereal t = m_thermo->temperature();
424  if (t == m_temp) {
425  return;
426  }
427  if (t <= 0.0) {
428  throw CanteraError("PecosTransport::update_T",
429  "negative temperature "+fp2str(t));
430  }
431  m_temp = t;
432  m_logt = log(m_temp);
433  m_kbt = Boltzmann * m_temp;
434  m_sqrt_t = sqrt(m_temp);
435  m_t14 = sqrt(m_sqrt_t);
436  m_t32 = m_temp * m_sqrt_t;
437  m_sqrt_kbt = sqrt(Boltzmann*m_temp);
438 
439  // compute powers of log(T)
440  m_polytempvec[0] = 1.0;
441  m_polytempvec[1] = m_logt;
442  m_polytempvec[2] = m_logt*m_logt;
443  m_polytempvec[3] = m_logt*m_logt*m_logt;
444  m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt;
445 
446  // temperature has changed, so polynomial fits will need to be
447  // redone.
448  m_viscmix_ok = false;
449  m_spvisc_ok = false;
450  m_viscwt_ok = false;
451  m_spcond_ok = false;
452  m_diffmix_ok = false;
453  m_bindiff_ok = false;
454  m_abc_ok = false;
455  m_condmix_ok = false;
456 }
457 
458 /**
459  * @internal This is called the first time any transport property
460  * is requested from Mixture after the concentrations
461  * have changed.
462  */
463 void PecosTransport::update_C()
464 {
465  // signal that concentration-dependent quantities will need to
466  // be recomputed before use, and update the local mole
467  // fractions.
468 
469  m_viscmix_ok = false;
470  m_diffmix_ok = false;
471  m_condmix_ok = false;
472 
473  m_thermo->getMoleFractions(DATA_PTR(m_molefracs));
474 
475  // add an offset to avoid a pure species condition
476  int k;
477  for (k = 0; k < m_nsp; k++) {
478  m_molefracs[k] = std::max(Tiny, m_molefracs[k]);
479  }
480 }
481 
482 /*************************************************************************
483  *
484  * methods to update temperature-dependent properties
485  *
486  *************************************************************************/
487 
488 /**
489  *
490  * Update the temperature-dependent parts of the mixture-averaged
491  * thermal conductivity.
492  *
493  * Calculated as,
494  * \f[
495  * k= \mu_s (5/2 * C_{v,s}^{trans} + C_{v,s}^{rot} + C_{v,s}^{vib}
496  * \f]
497  *
498  *
499  */
500 void PecosTransport::updateCond_T()
501 {
502 
503  int k;
504  doublereal fivehalves = 5/2;
505  for (k = 0; k < m_nsp; k++) {
506  // need to add cv_elec in the future
507  m_cond[k] = m_visc[k] * (fivehalves * cv_int[k] + cv_rot[k] + m_thermo->cv_vib(k,m_temp));
508  }
509  m_spcond_ok = true;
510  m_condmix_ok = false;
511 }
512 
513 
514 /**
515  * Update the binary diffusion coefficients. These are evaluated
516  * from the polynomial fits at unit pressure (1 Pa).
517  */
518 void PecosTransport::updateDiff_T()
519 {
520 
521  // evaluate binary diffusion coefficients at unit pressure
522  int i,j;
523  int ic = 0;
524  if (m_mode == CK_Mode) {
525  for (i = 0; i < m_nsp; i++) {
526  for (j = i; j < m_nsp; j++) {
527  m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic]));
528  m_bdiff(j,i) = m_bdiff(i,j);
529  ic++;
530  }
531  }
532  } else {
533  for (i = 0; i < m_nsp; i++) {
534  for (j = i; j < m_nsp; j++) {
535  m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec,
536  m_diffcoeffs[ic]);
537  m_bdiff(j,i) = m_bdiff(i,j);
538  ic++;
539  }
540  }
541  }
542 
543  m_bindiff_ok = true;
544  m_diffmix_ok = false;
545 }
546 
547 
548 /**
549  *
550  * Update the pure-species viscosities. (Pa-s) = (kg/m/sec)
551  *
552  * Using Blottner fit for viscosity. Defines kinematic viscosity
553  * of the form
554  * \f[
555  * \mu_s\left(T\right) = 0.10 \exp\left(A_s\left(\log T\right)^2 + B_s\log T + C_s\right)
556  * \f]
557  * where \f$ A_s \f$, \f$ B_s \f$, and \f$ C_s \f$ are constants.
558  *
559  */
560 void PecosTransport::updateSpeciesViscosities()
561 {
562 
563  // blottner
564  // return 0.10*std::exp(_a*(logT*logT) + _b*logT + _c);
565 
566  int k;
567  // iterate over species, update pure-species viscosity
568  for (k = 0; k < m_nsp; k++) {
569  m_visc[k] = 0.10*std::exp(a[k]*(m_logt*m_logt) + b[k]*m_logt + c[k]);
570  m_sqvisc[k] = sqrt(m_visc[k]);
571  }
572 
573  // time to update mixing
574  m_spvisc_ok = true;
575 }
576 
577 /*
578  * read_blottner_transport_table()
579  * loads up A B and C for blottner fits
580  * hardcoded for air, will need to generalize later
581  */
582 
583 void PecosTransport::read_blottner_transport_table()
584 {
585  // istringstream blot
586  // ("Air 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
587  // "CPAir 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
588  // "N 1.15572000000e-02 6.03167900000e-01 -1.24327495000e+01\n"
589  // "N2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
590  // "CPN2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
591  // "NO 4.36378000000e-02 -3.35511000000e-02 -9.57674300000e+00\n"
592  // "O 2.03144000000e-02 4.29440400000e-01 -1.16031403000e+01\n"
593  // "O2 4.49290000000e-02 -8.26158000000e-02 -9.20194750000e+00\n"
594  // "C -8.3285e-3 0.7703240 -12.7378000\n"
595  // "C2 -8.4311e-3 0.7876060 -13.0268000\n"
596  // "C3 -8.4312e-3 0.7876090 -12.8240000\n"
597  // "C2H -2.4241e-2 1.0946550 -14.5835500\n"
598  // "CN -8.3811e-3 0.7860330 -12.9406000\n"
599  // "CO -0.019527394 1.013295 -13.97873\n"
600  // "CO2 -0.019527387 1.047818 -14.32212\n"
601  // "HCN -2.4241e-2 1.0946550 -14.5835500\n"
602  // "H -8.3912e-3 0.7743270 -13.6653000\n"
603  // "H2 -8.3346e-3 0.7815380 -13.5351000\n"
604  // "e 0.00000000000e+00 0.00000000000e+00 -1.16031403000e+01\n");
605 
606  //
607  // from: AIAA-1997-2474 and Sandia Report SC-RR-70-754
608  //
609  // # Air -- Identical to N2 fit
610  // # N -- Sandia Report SC-RR-70-754
611  // # N2 -- Sandia Report SC-RR-70-754
612  // # CPN2 -- Identical to N2 fit
613  // # NO -- Sandia Report SC-RR-70-754
614  // # O -- Sandia Report SC-RR-70-754
615  // # O2 -- Sandia Report SC-RR-70-754
616  // # C -- AIAA-1997-2474
617  // # C2 -- AIAA-1997-2474
618  // # C3 -- AIAA-1997-2474
619  // # C2H -- wild-ass guess: identical to HCN fit
620  // # CN -- AIAA-1997-2474
621  // # CO -- AIAA-1997-2474
622  // # CO2 -- AIAA-1997-2474
623  // # HCN -- AIAA-1997-2474
624  // # H -- AIAA-1997-2474
625  // # H2 -- AIAA-1997-2474
626  // # e -- Sandia Report SC-RR-70-754
627 
628  istringstream blot
629  ("Air 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
630  "CPAir 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
631  "N 1.15572000000e-02 6.03167900000e-01 -1.24327495000e+01\n"
632  "N2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
633  "CPN2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n"
634  "NO 4.36378000000e-02 -3.35511000000e-02 -9.57674300000e+00\n"
635  "O 2.03144000000e-02 4.29440400000e-01 -1.16031403000e+01\n"
636  "O2 4.49290000000e-02 -8.26158000000e-02 -9.20194750000e+00\n"
637  "C -8.3285e-3 0.7703240 -12.7378000\n"
638  "C2 -8.4311e-3 0.7876060 -13.0268000\n"
639  "C3 -8.4312e-3 0.7876090 -12.8240000\n"
640  "C2H -2.4241e-2 1.0946550 -14.5835500\n"
641  "CN -8.3811e-3 0.7860330 -12.9406000\n"
642  "CO -0.019527394 1.013295 -13.97873\n"
643  "CO2 -0.019527387 1.047818 -14.32212\n"
644  "HCN -2.4241e-2 1.0946550 -14.5835500\n"
645  "H -8.3912e-3 0.7743270 -13.6653000\n"
646  "H2 -8.3346e-3 0.7815380 -13.5351000\n"
647  "e 0.00000000000e+00 0.00000000000e+00 -1.16031403000e+01\n");
648 
649  string line;
650  string name;
651  string ss1,ss2,ss3,ss4,sss;
652  int k;
653  int i = 0;
654 
655  while (std::getline(blot, line)) {
656 
657  istringstream ss(line);
658  std::getline(ss, ss1, ' ');
659  std::getline(ss, ss2, ' ');
660  std::getline(ss, ss3, ' ');
661  std::getline(ss, ss4, ' ');
662  name = ss1;
663 
664  // now put coefficients in correct species
665  for (k = 0; k < m_nsp; k++) {
666  string sss = m_thermo->speciesName(k);
667 
668  // this is the right species index
669  if (sss.compare(ss1) == 0) {
670  a[k] = atof(ss2.c_str());
671  b[k] = atof(ss3.c_str());
672  c[k] = atof(ss4.c_str());
673 
674  // index
675  i++;
676  } else { // default to air
677 
678  a[k] = 0.026;
679  b[k] = 0.3;
680  c[k] = -11.3;
681  }
682 
683  } // done with for loop
684  }
685 
686 
687  // for (k = 0; k < m_nsp; k++)
688  // {
689  // string sss = m_thermo->speciesName(k);
690  // cout << sss << endl;
691  // cout << a[k] << endl;
692  // cout << b[k] << endl;
693  // cout << c[k] << endl;
694  // }
695 
696  // simple sanity check
697  // if(i != m_nsp-1)
698  // {
699  // std::cout << "error\n" << i << std::endl;
700  // }
701 
702 }
703 
704 /**
705  *
706  * Update the temperature-dependent viscosity terms.
707  * Updates the array of pure species viscosities, and the
708  * weighting functions in the viscosity mixture rule.
709  * The flag m_visc_ok is set to true.
710  *
711  */
712 void PecosTransport::updateViscosity_T()
713 {
714  doublereal vratiokj, wratiojk, factor1;
715 
716  if (!m_spvisc_ok) {
717  updateSpeciesViscosities();
718  }
719 
720  // see Eq. (9-5.15) of Reid, Prausnitz, and Poling
721  int j, k;
722  for (j = 0; j < m_nsp; j++) {
723  for (k = j; k < m_nsp; k++) {
724  vratiokj = m_visc[k]/m_visc[j];
725  wratiojk = m_mw[j]/m_mw[k];
726 
727  // Note that m_wratjk(k,j) holds the square root of
728  // m_wratjk(j,k)!
729  factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
730  m_phi(k,j) = factor1*factor1 /
731  (SqrtEight * m_wratkj1(j,k));
732  m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
733  }
734  }
735  m_viscwt_ok = true;
736 }
737 
738 }
739 
Generated by   doxygen 1.8.2