Cantera  2.0
DustyGasTransport.cpp
Go to the documentation of this file.
1 /**
2  * @file DustyGasTransport.cpp
3  * Implementation file for class DustyGasTransport
4  *
5  * @ingroup transportProps
6  */
7 
8 /*
9  * Copyright 2003 California Institute of Technology
10  * See file License.txt for licensing information
11  */
12 
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-20
24 
25 namespace Cantera
26 {
27 
28 //====================================================================================================================
29 DustyGasTransport::DustyGasTransport(thermo_t* thermo) :
30  Transport(thermo),
31  m_mw(0),
32  m_dk(0),
33  m_temp(-1.0),
34  m_multidiff(0,0),
35  m_spwork(0),
36  m_spwork2(0),
37  m_gradP(0.0),
38  m_knudsen_ok(false),
39  m_bulk_ok(false),
40  m_porosity(0.0),
41  m_tortuosity(1.0),
42  m_pore_radius(0.0),
43  m_diam(0.0),
44  m_perm(-1.0),
45  m_gastran(0)
46 {
47 }
48 //====================================================================================================================
50  Transport(),
51  m_mw(0),
52  m_dk(0),
53  m_temp(-1.0),
54  m_multidiff(0,0),
55  m_spwork(0),
56  m_spwork2(0),
57  m_gradP(0.0),
58  m_knudsen_ok(false),
59  m_bulk_ok(false),
60  m_porosity(0.0),
61  m_tortuosity(1.0),
62  m_pore_radius(0.0),
63  m_diam(0.0),
64  m_perm(-1.0),
65  m_gastran(0)
66 {
67  *this = right;
68 }
69 //====================================================================================================================
70 // Assignment operator
71 /*
72  * This is NOT a virtual function.
73  *
74  * @param right Reference to %DustyGasTransport object to be copied
75  * into the current one.
76  */
78 {
79  if (&right == this) {
80  return *this;
81  }
82  Transport::operator=(right);
83 
84  m_mw = right.m_mw;
85  m_d = right.m_d;
86  m_x = right.m_x;
87  m_dk = right.m_dk;
88  m_temp = right.m_temp;
89  m_multidiff = right.m_multidiff;
90  m_spwork = right.m_spwork;
91  m_spwork2 = right.m_spwork2;
92  m_gradP = right.m_gradP;
93  m_knudsen_ok = right.m_knudsen_ok;
94  m_bulk_ok= right.m_bulk_ok;
95  m_porosity = right.m_porosity;
96  m_tortuosity = right.m_tortuosity;
98  m_diam = right.m_diam;
99  m_perm = right.m_perm;
100 
101  // Warning -> gastran may not point to the correct object
102  // after this copy. The routine initialize() must be called
103  if (m_gastran) {
104  delete m_gastran;
105  }
107 
108 
109  return *this;
110 }
111 //====================================================================================================================
113 {
114  delete m_gastran;
115 }
116 //====================================================================================================================
117 // Duplication routine for objects which inherit from %Transport
118 /*
119  * This virtual routine can be used to duplicate %Transport objects
120  * inherited from %Transport even if the application only has
121  * a pointer to %Transport to work with.
122  *
123  * These routines are basically wrappers around the derived copy
124  * constructor.
125  */
127 {
128  DustyGasTransport* tr = new DustyGasTransport(*this);
129  return (dynamic_cast<Transport*>(tr));
130 }
131 //====================================================================================================================
132 // Set the Parameters in the model
133 /*
134  * @param type Type of the parameter to set
135  * 0 - porosity
136  * 1 - tortuosity
137  * 2 - mean pore radius
138  * 3 - mean particle radius
139  * 4 - permeability
140  * @param k Unused int
141  * @param p pointer to double for the input list of parameters
142  *
143  */
144 void DustyGasTransport::setParameters(const int type, const int k, const doublereal* const p)
145 {
146  switch (type) {
147  case 0:
148  setPorosity(p[0]);
149  break;
150  case 1:
151  setTortuosity(p[0]);
152  break;
153  case 2:
154  setMeanPoreRadius(p[0]);
155  break;
156  case 3:
158  break;
159  case 4:
160  setPermeability(p[0]);
161  break;
162  default:
163  throw CanteraError("DustyGasTransport::init", "unknown parameter");
164  }
165 }
166 //====================================================================================================================
167 // Initialization routine called by TransportFactory
168 /*
169  * The DustyGas model is a subordinate model to the gas phase transport model. Here we
170  * set the gas phase models.
171  *
172  * This is a protected routine, so that initialization of the Model must occur within Cantera's setup
173  *
174  * @param phase Pointer to the underlying ThermoPhase model for the gas phase
175  * @param gastr Pointer to the underlying Transport model for transport in the gas phase.
176  */
178 {
179 
180  // constant mixture attributes
181  m_thermo = phase;
182  m_nsp = m_thermo->nSpecies();
183  if (m_gastran != gastr) {
184  if (m_gastran) {
185  delete m_gastran;
186  }
187  m_gastran = gastr;
188  }
189 
190  // make a local copy of the molecular weights
191  m_mw.resize(m_nsp);
192  copy(m_thermo->molecularWeights().begin(), m_thermo->molecularWeights().end(), m_mw.begin());
193 
195  m_d.resize(m_nsp, m_nsp);
196  m_dk.resize(m_nsp, 0.0);
197 
198  m_x.resize(m_nsp, 0.0);
200 
201  // set flags all false
202  m_knudsen_ok = false;
203  m_bulk_ok = false;
204 
205  m_spwork.resize(m_nsp);
206  m_spwork2.resize(m_nsp);
207 }
208 //====================================================================================================================
209 // Private routine to update the dusty gas binary diffusion coefficients
210 /*
211  * The dusty gas binary diffusion coefficients \f$ D^{dg}_{i,j} \f$ are evaluated from the binary
212  * gas-phase diffusion coefficients \f$ D^{bin}_{i,j} \f$ using the following formula
213  *
214  * \f[
215  * D^{dg}_{i,j} = \frac{\phi}{\tau} D^{bin}_{i,j}
216  * \f]
217  *
218  * where \f$ \phi \f$ is the porosity of the media and \f$ \tau \f$ is the tortuosity of the media.
219  *
220  */
222 {
223  if (m_bulk_ok) {
224  return;
225  }
226 
227  // get the gaseous binary diffusion coefficients
229  doublereal por2tort = m_porosity / m_tortuosity;
230  for (size_t n = 0; n < m_nsp; n++) {
231  for (size_t m = 0; m < m_nsp; m++) {
232  m_d(n,m) *= por2tort;
233  }
234  }
235  m_bulk_ok = true;
236 }
237 //====================================================================================================================
238 // Private routine to update the Knudsen diffusion coefficients
239 /*
240  * The Knudsen diffusion coefficients are given by the following form
241  *
242  * \f[
243  * \mathcal{D}^{knud}_k = \frac{2}{3} \frac{r_{pore} \phi}{\tau} \left( \frac{8 R T}{\pi W_k} \right)^{1/2}
244  * \f]
245  *
246  */
248 {
249  if (m_knudsen_ok) {
250  return;
251  }
252  doublereal K_g = m_pore_radius * m_porosity / m_tortuosity;
253  const doublereal TwoThirds = 2.0/3.0;
254  for (size_t k = 0; k < m_nsp; k++) {
255  m_dk[k] = TwoThirds * K_g * sqrt((8.0 * GasConstant * m_temp)/
256  (Pi * m_mw[k]));
257  }
258  m_knudsen_ok = true;
259 }
260 
261 //====================================================================================================================
262 // Private routine to calculate the H matrix
263 /*
264  * The H matrix is the term we have given to the matrix of coefficients in the equation for the molar
265  * fluxes. The matrix must be inverted in order to calculate the molar fluxes.
266  *
267  * The multicomponent diffusion H matrix \f$ H_{k,l} \f$ is given by the following formulas
268  *
269  * \f[
270  * H_{k,l} = - \frac{X_k}{D^e_{k,l}}
271  * \f]
272  * \f[
273  * H_{k,k} = \frac{1}{\mathcal(D)^{e}_{k, knud}} + \sum_{j \ne k}^N{ \frac{X_j}{D^e_{k,j}} }
274  * \f]
275  */
277 {
280  doublereal sum;
281  for (size_t k = 0; k < m_nsp; k++) {
282 
283  // evaluate off-diagonal terms
284  for (size_t l = 0; l < m_nsp; l++) {
285  m_multidiff(k,l) = -m_x[k]/m_d(k,l);
286  }
287 
288  // evaluate diagonal term
289  sum = 0.0;
290  for (size_t j = 0; j < m_nsp; j++) {
291  if (j != k) {
292  sum += m_x[j]/m_d(k,j);
293  }
294  }
295  m_multidiff(k,k) = 1.0/m_dk[k] + sum;
296  }
297 }
298 //====================================================================================================================
299 void DustyGasTransport::getMolarFluxes(const doublereal* const state1,
300  const doublereal* const state2,
301  const doublereal delta,
302  doublereal* const fluxes)
303 {
304 
305  doublereal conc1, conc2;
306 
307  // cbar will be the average concentration between the two points
308  doublereal* const cbar = DATA_PTR(m_spwork);
309  doublereal* const gradc = DATA_PTR(m_spwork2);
310  const doublereal t1 = state1[0];
311  const doublereal t2 = state2[0];
312  const doublereal rho1 = state1[1];
313  const doublereal rho2 = state2[1];
314  const doublereal* const y1 = state1 + 2;
315  const doublereal* const y2 = state2 + 2;
316  doublereal c1sum = 0.0, c2sum = 0.0;
317 
318  for (size_t k = 0; k < m_nsp; k++) {
319  conc1 = rho1 * y1[k] / m_mw[k];
320  conc2 = rho2 * y2[k] / m_mw[k];
321  cbar[k] = 0.5*(conc1 + conc2);
322  gradc[k] = (conc2 - conc1) / delta;
323  c1sum += conc1;
324  c2sum += conc2;
325  }
326 
327  // Calculate the pressures at p1 p2 and pbar
328  doublereal p1 = c1sum * GasConstant * t1;
329  doublereal p2 = c2sum * GasConstant * t2;
330  doublereal pbar = 0.5*(p1 + p2);
331  doublereal gradp = (p2 - p1)/delta;
332  doublereal tbar = 0.5*(t1 + t2);
333 
334  m_thermo->setState_TPX(tbar, pbar, cbar);
335 
337 
338  // Multiply m_multidiff and gradc together and store the result in fluxes[]
339  multiply(m_multidiff, gradc, fluxes);
340 
341  divide_each(cbar, cbar + m_nsp, m_dk.begin());
342 
343  // if no permeability has been specified, use result for
344  // close-packed spheres
345  double b = 0.0;
346  if (m_perm < 0.0) {
347  double p = m_porosity;
348  double d = m_diam;
349  double t = m_tortuosity;
350  b = p*p*p*d*d/(72.0*t*(1.0-p)*(1.0-p));
351  } else {
352  b = m_perm;
353  }
354  b *= gradp / m_gastran->viscosity();
355  scale(cbar, cbar + m_nsp, cbar, b);
356 
357  // Multiply m_multidiff with cbar and add it to fluxes
358  increment(m_multidiff, cbar, fluxes);
359  scale(fluxes, fluxes + m_nsp, fluxes, -1.0);
360 }
361 //====================================================================================================================
362 // Private routine to update the Multicomponent diffusion coefficients that are used in the approximation
363 /*
364  * This routine updates the H matrix and then inverts it.
365  */
367 {
368  // see if temperature has changed
370 
371  // update the mole fractions
373 
374  eval_H_matrix();
375 
376  // invert H
377  int ierr = invert(m_multidiff);
378 
379  if (ierr != 0) {
380  throw CanteraError("DustyGasTransport::updateMultiDiffCoeffs",
381  "invert returned ierr = "+int2str(ierr));
382  }
383 }
384 //====================================================================================================================
385 // Return the Multicomponent diffusion coefficients. Units: [m^2/s].
386 /*
387  * Returns the array of multicomponent diffusion coefficients.
388  *
389  * @param ld The dimension of the inner loop of d (usually equal to m_nsp)
390  * @param d flat vector of diffusion coefficients, fortran ordering.
391  * d[ld*j+i] is the D_ij diffusion coefficient (the diffusion
392  * coefficient for species i due to species j).
393  */
394 void DustyGasTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
395 {
397  for (size_t i = 0; i < m_nsp; i++) {
398  for (size_t j = 0; j < m_nsp; j++) {
399  d[ld*j + i] = m_multidiff(i,j);
400  }
401  }
402 }
403 //====================================================================================================================
404 // Update temperature-dependent quantities within the object
405 /*
406  * The object keeps a value m_temp, which is the temperature at which quantities were last evaluated
407  * at. If the temperature is changed, update Booleans are set false, triggering recomputation.
408  */
410 {
411  if (m_temp == m_thermo->temperature()) {
412  return;
413  }
415  m_knudsen_ok = false;
416  m_bulk_ok = false;
417 }
418 //====================================================================================================================
420 {
422 
423  // add an offset to avoid a pure species condition
424  // (check - this may be unnecessary)
425  for (size_t k = 0; k < m_nsp; k++) {
426  m_x[k] = std::max(MIN_X, m_x[k]);
427  }
428  // diffusion coeffs depend on Pressure
429  m_bulk_ok = false;
430 }
431 //====================================================================================================================
432 // Set the porosity (dimensionless)
433 /*
434  * @param porosity Set the value of the porosity
435  */
436 void DustyGasTransport::setPorosity(doublereal porosity)
437 {
438  m_porosity = porosity;
439  m_knudsen_ok = false;
440  m_bulk_ok = false;
441 }
442 //====================================================================================================================
443 // Set the tortuosity (dimensionless)
444 /*
445  * @param tort Value of the tortuosity
446  */
447 void DustyGasTransport::setTortuosity(doublereal tort)
448 {
449  m_tortuosity = tort;
450  m_knudsen_ok = false;
451  m_bulk_ok = false;
452 }
453 //====================================================================================================================
454 // Set the mean pore radius (m)
455 /*
456  * @param rbar Value of the pore radius ( m)
457  */
459 {
460  m_pore_radius = rbar;
461  m_knudsen_ok = false;
462 }
463 //====================================================================================================================
464 // Set the mean particle diameter
465 /*
466  * @param dbar Set the mean particle diameter (m)
467  */
469 {
470  m_diam = dbar;
471 }
472 //====================================================================================================================
473 // Set the permeability of the media
474 /*
475  * If not set, the value for close-packed spheres will be used by default.
476  *
477  * The value for close-packed spheres is given below, where p is the porosity,
478  * t is the tortuosity, and d is the diameter of the sphere
479  *
480  * \f[
481  * \kappa = \frac{p^3 d^2}{72 t (1 - p)^2}
482  * \f]
483  *
484  * @param B set the permeability of the media (units = m^2)
485  */
487 {
488  m_perm = B;
489 }
490 //====================================================================================================================
491 // Return a reference to the transport manager used to compute the gas
492 // binary diffusion coefficients and the viscosity.
493 /*
494  * @return Returns a reference to the gas transport object
495  */
497 {
498  return *m_gastran;
499 }
500 
501 //====================================================================================================================
502 }