Cantera  2.1.2
DustyGasTransport.cpp
Go to the documentation of this file.
1 /**
2  * @file DustyGasTransport.cpp
3  * Implementation file for class DustyGasTransport
4  */
5 
6 /*
7  * Copyright 2003 California Institute of Technology
8  * See file License.txt for licensing information
9  */
10 
14 
15 using namespace std;
16 
17 namespace Cantera
18 {
19 DustyGasTransport::DustyGasTransport(thermo_t* thermo) :
20  Transport(thermo),
21  m_mw(0),
22  m_dk(0),
23  m_temp(-1.0),
24  m_multidiff(0,0),
25  m_spwork(0),
26  m_spwork2(0),
27  m_gradP(0.0),
28  m_knudsen_ok(false),
29  m_bulk_ok(false),
30  m_porosity(0.0),
31  m_tortuosity(1.0),
32  m_pore_radius(0.0),
33  m_diam(0.0),
34  m_perm(-1.0),
35  m_gastran(0)
36 {
37 }
38 
40  Transport(),
41  m_mw(0),
42  m_dk(0),
43  m_temp(-1.0),
44  m_multidiff(0,0),
45  m_spwork(0),
46  m_spwork2(0),
47  m_gradP(0.0),
48  m_knudsen_ok(false),
49  m_bulk_ok(false),
50  m_porosity(0.0),
51  m_tortuosity(1.0),
52  m_pore_radius(0.0),
53  m_diam(0.0),
54  m_perm(-1.0),
55  m_gastran(0)
56 {
57  *this = right;
58 }
59 
61 {
62  if (&right == this) {
63  return *this;
64  }
65  Transport::operator=(right);
66 
67  m_mw = right.m_mw;
68  m_d = right.m_d;
69  m_x = right.m_x;
70  m_dk = right.m_dk;
71  m_temp = right.m_temp;
72  m_multidiff = right.m_multidiff;
73  m_spwork = right.m_spwork;
74  m_spwork2 = right.m_spwork2;
75  m_gradP = right.m_gradP;
76  m_knudsen_ok = right.m_knudsen_ok;
77  m_bulk_ok= right.m_bulk_ok;
78  m_porosity = right.m_porosity;
79  m_tortuosity = right.m_tortuosity;
81  m_diam = right.m_diam;
82  m_perm = right.m_perm;
83 
84  // Warning -> gastran may not point to the correct object
85  // after this copy. The routine initialize() must be called
86  delete m_gastran;
88 
89 
90  return *this;
91 }
92 
93 DustyGasTransport::~DustyGasTransport()
94 {
95  delete m_gastran;
96 }
97 
99 {
100  DustyGasTransport* tr = new DustyGasTransport(*this);
101  return dynamic_cast<Transport*>(tr);
102 }
103 
105 {
106 
107  Transport::setThermo(thermo);
108  m_gastran->setThermo(thermo);
109 }
110 
111 void DustyGasTransport::setParameters(const int type, const int k, const doublereal* const p)
112 {
113  warn_deprecated("DustyGasTransport::setParameters");
114  switch (type) {
115  case 0:
116  setPorosity(p[0]);
117  break;
118  case 1:
119  setTortuosity(p[0]);
120  break;
121  case 2:
122  setMeanPoreRadius(p[0]);
123  break;
124  case 3:
126  break;
127  case 4:
128  setPermeability(p[0]);
129  break;
130  default:
131  throw CanteraError("DustyGasTransport::init", "unknown parameter");
132  }
133 }
134 
136 {
137 
138  // constant mixture attributes
139  m_thermo = phase;
140  m_nsp = m_thermo->nSpecies();
141  if (m_gastran != gastr) {
142  delete m_gastran;
143  m_gastran = gastr;
144  }
145 
146  // make a local copy of the molecular weights
147  m_mw.resize(m_nsp);
148  copy(m_thermo->molecularWeights().begin(), m_thermo->molecularWeights().end(), m_mw.begin());
149 
151  m_d.resize(m_nsp, m_nsp);
152  m_dk.resize(m_nsp, 0.0);
153 
154  m_x.resize(m_nsp, 0.0);
156 
157  // set flags all false
158  m_knudsen_ok = false;
159  m_bulk_ok = false;
160 
161  m_spwork.resize(m_nsp);
162  m_spwork2.resize(m_nsp);
163 }
164 
166 {
167  if (m_bulk_ok) {
168  return;
169  }
170 
171  // get the gaseous binary diffusion coefficients
173  doublereal por2tort = m_porosity / m_tortuosity;
174  for (size_t n = 0; n < m_nsp; n++) {
175  for (size_t m = 0; m < m_nsp; m++) {
176  m_d(n,m) *= por2tort;
177  }
178  }
179  m_bulk_ok = true;
180 }
181 
183 {
184  if (m_knudsen_ok) {
185  return;
186  }
187  doublereal K_g = m_pore_radius * m_porosity / m_tortuosity;
188  const doublereal TwoThirds = 2.0/3.0;
189  for (size_t k = 0; k < m_nsp; k++) {
190  m_dk[k] = TwoThirds * K_g * sqrt((8.0 * GasConstant * m_temp)/
191  (Pi * m_mw[k]));
192  }
193  m_knudsen_ok = true;
194 }
195 
197 {
200  doublereal sum;
201  for (size_t k = 0; k < m_nsp; k++) {
202 
203  // evaluate off-diagonal terms
204  for (size_t l = 0; l < m_nsp; l++) {
205  m_multidiff(k,l) = -m_x[k]/m_d(k,l);
206  }
207 
208  // evaluate diagonal term
209  sum = 0.0;
210  for (size_t j = 0; j < m_nsp; j++) {
211  if (j != k) {
212  sum += m_x[j]/m_d(k,j);
213  }
214  }
215  m_multidiff(k,k) = 1.0/m_dk[k] + sum;
216  }
217 }
218 
219 void DustyGasTransport::getMolarFluxes(const doublereal* const state1,
220  const doublereal* const state2,
221  const doublereal delta,
222  doublereal* const fluxes)
223 {
224  doublereal conc1, conc2;
225 
226  // cbar will be the average concentration between the two points
227  doublereal* const cbar = DATA_PTR(m_spwork);
228  doublereal* const gradc = DATA_PTR(m_spwork2);
229  const doublereal t1 = state1[0];
230  const doublereal t2 = state2[0];
231  const doublereal rho1 = state1[1];
232  const doublereal rho2 = state2[1];
233  const doublereal* const y1 = state1 + 2;
234  const doublereal* const y2 = state2 + 2;
235  doublereal c1sum = 0.0, c2sum = 0.0;
236 
237  for (size_t k = 0; k < m_nsp; k++) {
238  conc1 = rho1 * y1[k] / m_mw[k];
239  conc2 = rho2 * y2[k] / m_mw[k];
240  cbar[k] = 0.5*(conc1 + conc2);
241  gradc[k] = (conc2 - conc1) / delta;
242  c1sum += conc1;
243  c2sum += conc2;
244  }
245 
246  // Calculate the pressures at p1 p2 and pbar
247  doublereal p1 = c1sum * GasConstant * t1;
248  doublereal p2 = c2sum * GasConstant * t2;
249  doublereal pbar = 0.5*(p1 + p2);
250  doublereal gradp = (p2 - p1)/delta;
251  doublereal tbar = 0.5*(t1 + t2);
252 
253  m_thermo->setState_TPX(tbar, pbar, cbar);
254 
256 
257  // Multiply m_multidiff and gradc together and store the result in fluxes[]
258  multiply(m_multidiff, gradc, fluxes);
259 
260  divide_each(cbar, cbar + m_nsp, m_dk.begin());
261 
262  // if no permeability has been specified, use result for
263  // close-packed spheres
264  double b = 0.0;
265  if (m_perm < 0.0) {
266  double p = m_porosity;
267  double d = m_diam;
268  double t = m_tortuosity;
269  b = p*p*p*d*d/(72.0*t*(1.0-p)*(1.0-p));
270  } else {
271  b = m_perm;
272  }
273  b *= gradp / m_gastran->viscosity();
274  scale(cbar, cbar + m_nsp, cbar, b);
275 
276  // Multiply m_multidiff with cbar and add it to fluxes
277  increment(m_multidiff, cbar, fluxes);
278  scale(fluxes, fluxes + m_nsp, fluxes, -1.0);
279 }
280 
282 {
283  // see if temperature has changed
285 
286  // update the mole fractions
288 
289  eval_H_matrix();
290 
291  // invert H
292  int ierr = invert(m_multidiff);
293 
294  if (ierr != 0) {
295  throw CanteraError("DustyGasTransport::updateMultiDiffCoeffs",
296  "invert returned ierr = "+int2str(ierr));
297  }
298 }
299 
300 void DustyGasTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
301 {
303  for (size_t i = 0; i < m_nsp; i++) {
304  for (size_t j = 0; j < m_nsp; j++) {
305  d[ld*j + i] = m_multidiff(i,j);
306  }
307  }
308 }
309 
311 {
312  if (m_temp == m_thermo->temperature()) {
313  return;
314  }
316  m_knudsen_ok = false;
317  m_bulk_ok = false;
318 }
319 
321 {
323 
324  // add an offset to avoid a pure species condition
325  // (check - this may be unnecessary)
326  for (size_t k = 0; k < m_nsp; k++) {
327  m_x[k] = std::max(Tiny, m_x[k]);
328  }
329  // diffusion coeffs depend on Pressure
330  m_bulk_ok = false;
331 }
332 
333 void DustyGasTransport::setPorosity(doublereal porosity)
334 {
335  m_porosity = porosity;
336  m_knudsen_ok = false;
337  m_bulk_ok = false;
338 }
339 
340 void DustyGasTransport::setTortuosity(doublereal tort)
341 {
342  m_tortuosity = tort;
343  m_knudsen_ok = false;
344  m_bulk_ok = false;
345 }
346 
348 {
349  m_pore_radius = rbar;
350  m_knudsen_ok = false;
351 }
352 
354 {
355  m_diam = dbar;
356 }
357 
359 {
360  m_perm = B;
361 }
362 
364 {
365  return *m_gastran;
366 }
367 
368 }
doublereal m_pore_radius
Pore radius (meter)
virtual void getMultiDiffCoeffs(const size_t ld, doublereal *const d)
Return the Multicomponent diffusion coefficients. Units: [m^2/s].
Headers for the DustyGasTransport object, which models transport properties in porous media using the...
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
int invert(DenseMatrix &A, size_t nn)
invert A. A is overwritten with A^-1.
void divide_each(OutputIter x_begin, OutputIter x_end, InputIter y_begin)
Templated divide of each element of x by the corresponding element of y.
Definition: utilities.h:369
virtual void getBinaryDiffCoeffs(const size_t ld, doublereal *const d)
Returns the matrix of binary diffusion coefficients [m^2/s].
void updateTransport_T()
Update temperature-dependent quantities within the object.
virtual void setThermo(thermo_t &thermo)
Specifies the ThermoPhase object.
Transport & gasTransport()
Return a reference to the transport manager used to compute the gas binary diffusion coefficients and...
virtual void setState_TPX(doublereal t, doublereal p, const doublereal *x)
Set the temperature (K), pressure (Pa), and mole fractions.
doublereal m_diam
Particle diameter.
vector_fp m_x
mole fractions
virtual Transport * duplMyselfAsTransport() const
Duplication routine for objects which inherit from Transport.
thermo_t * m_thermo
pointer to the object representing the phase
Class DustyGasTransport implements the Dusty Gas model for transport in porous media.
DenseMatrix m_d
binary diffusion coefficients
DustyGasTransport(thermo_t *thermo=0)
default constructor
Base class for transport property managers.
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:76
DustyGasTransport & operator=(const DustyGasTransport &right)
Assignment operator.
void increment(const DenseMatrix &A, const double *b, double *prod)
Multiply A*b and add it to the result in prod. Uses BLAS routine DGEMV.
vector_fp m_dk
Knudsen diffusion coefficients.
const doublereal Pi
Pi.
Definition: ct_defs.h:51
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:363
void updateMultiDiffCoeffs()
Update the Multicomponent diffusion coefficients that are used in the approximation.
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:519
void eval_H_matrix()
Calculate the H matrix.
doublereal m_perm
Permeability of the media.
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:101
doublereal m_gradP
Pressure Gradient.
void updateTransport_C()
Update concentration-dependent quantities within the object.
void setPorosity(doublereal porosity)
Set the porosity (dimensionless)
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
Definition: DenseMatrix.cpp:71
bool m_knudsen_ok
Update-to-date variable for Knudsen diffusion coefficients.
void multiply(const DenseMatrix &A, const double *const b, double *const prod)
Multiply A*b and return the result in prod. Uses BLAS routine DGEMV.
virtual void setParameters(const int type, const int k, const doublereal *const p)
Set the Parameters in the model.
vector_fp m_spwork2
work space of size m_nsp;
void initialize(ThermoPhase *phase, Transport *gastr)
Initialization routine called by TransportFactory.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
virtual doublereal viscosity()
void setTortuosity(doublereal tort)
Set the tortuosity (dimensionless)
virtual void getMolarFluxes(const doublereal *const state1, const doublereal *const state2, const doublereal delta, doublereal *const fluxes)
Get the molar fluxes [kmol/m^2/s], given the thermodynamic state at two nearby points.
void setMeanPoreRadius(doublereal rbar)
Set the mean pore radius (m)
Transport * m_gastran
Pointer to the transport object for the gas phase.
doublereal m_porosity
Porosity.
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:252
const vector_fp & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:505
doublereal temperature() const
Temperature (K).
Definition: Phase.h:528
doublereal m_tortuosity
Tortuosity.
void setPermeability(doublereal B)
Set the permeability of the media.
bool m_bulk_ok
Update-to-date variable for Binary diffusion coefficients.
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:154
const doublereal Tiny
Small number to compare differences of mole fractions against.
Definition: ct_defs.h:155
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:66
Contains declarations for string manipulation functions within Cantera.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
size_t m_nsp
Number of species.
doublereal m_temp
temperature
virtual void setThermo(thermo_t &thermo)
Specifies the ThermoPhase object.
Header file for class ThermoPhase, the base class for phases with thermodynamic properties, and the text for the Module thermoprops (see Thermodynamic Properties and class ThermoPhase).
void updateBinaryDiffCoeffs()
Private routine to update the dusty gas binary diffusion coefficients.
DenseMatrix m_multidiff
Multicomponent diffusion coefficients.
vector_fp m_mw
Local copy of the species molecular weights.
vector_fp m_spwork
work space of size m_nsp;
void setMeanParticleDiameter(doublereal dbar)
Set the mean particle diameter.
void updateKnudsenDiffCoeffs()
Update the Knudsen diffusion coefficients.