Cantera  2.4.0
DustyGasTransport.cpp
Go to the documentation of this file.
1 /**
2  * @file DustyGasTransport.cpp
3  * Implementation file for class DustyGasTransport
4  */
5 
6 // This file is part of Cantera. See License.txt in the top-level directory or
7 // at http://www.cantera.org/license.txt for license and copyright information.
8 
11 
12 using namespace std;
13 
14 namespace Cantera
15 {
16 DustyGasTransport::DustyGasTransport(thermo_t* thermo) :
17  Transport(thermo),
18  m_temp(-1.0),
19  m_gradP(0.0),
20  m_knudsen_ok(false),
21  m_bulk_ok(false),
22  m_porosity(0.0),
23  m_tortuosity(1.0),
24  m_pore_radius(0.0),
25  m_diam(0.0),
26  m_perm(-1.0)
27 {
28 }
29 
31 {
33  m_gastran->setThermo(thermo);
34 }
35 
37 {
38  // constant mixture attributes
39  m_thermo = phase;
40  m_nsp = m_thermo->nSpecies();
41  if (m_gastran.get() != gastr) {
42  m_gastran.reset(gastr);
43  }
44 
45  // make a local copy of the molecular weights
47 
50  m_dk.resize(m_nsp, 0.0);
51 
52  m_x.resize(m_nsp, 0.0);
53  m_thermo->getMoleFractions(m_x.data());
54 
55  // set flags all false
56  m_knudsen_ok = false;
57  m_bulk_ok = false;
58 
59  m_spwork.resize(m_nsp);
60  m_spwork2.resize(m_nsp);
61 }
62 
64 {
65  if (m_bulk_ok) {
66  return;
67  }
68 
69  // get the gaseous binary diffusion coefficients
70  m_gastran->getBinaryDiffCoeffs(m_nsp, m_d.ptrColumn(0));
71  doublereal por2tort = m_porosity / m_tortuosity;
72  for (size_t n = 0; n < m_nsp; n++) {
73  for (size_t m = 0; m < m_nsp; m++) {
74  m_d(n,m) *= por2tort;
75  }
76  }
77  m_bulk_ok = true;
78 }
79 
81 {
82  if (m_knudsen_ok) {
83  return;
84  }
85  doublereal K_g = m_pore_radius * m_porosity / m_tortuosity;
86  for (size_t k = 0; k < m_nsp; k++) {
87  m_dk[k] = 2.0/3.0 * K_g * sqrt((8.0 * GasConstant * m_temp)/
88  (Pi * m_mw[k]));
89  }
90  m_knudsen_ok = true;
91 }
92 
94 {
97  for (size_t k = 0; k < m_nsp; k++) {
98  // evaluate off-diagonal terms
99  for (size_t j = 0; j < m_nsp; j++) {
100  m_multidiff(k,j) = -m_x[k]/m_d(k,j);
101  }
102 
103  // evaluate diagonal term
104  double sum = 0.0;
105  for (size_t j = 0; j < m_nsp; j++) {
106  if (j != k) {
107  sum += m_x[j]/m_d(k,j);
108  }
109  }
110  m_multidiff(k,k) = 1.0/m_dk[k] + sum;
111  }
112 }
113 
114 void DustyGasTransport::getMolarFluxes(const doublereal* const state1,
115  const doublereal* const state2,
116  const doublereal delta,
117  doublereal* const fluxes)
118 {
119  // cbar will be the average concentration between the two points
120  doublereal* const cbar = m_spwork.data();
121  doublereal* const gradc = m_spwork2.data();
122  const doublereal t1 = state1[0];
123  const doublereal t2 = state2[0];
124  const doublereal rho1 = state1[1];
125  const doublereal rho2 = state2[1];
126  const doublereal* const y1 = state1 + 2;
127  const doublereal* const y2 = state2 + 2;
128  doublereal c1sum = 0.0, c2sum = 0.0;
129 
130  for (size_t k = 0; k < m_nsp; k++) {
131  double conc1 = rho1 * y1[k] / m_mw[k];
132  double conc2 = rho2 * y2[k] / m_mw[k];
133  cbar[k] = 0.5*(conc1 + conc2);
134  gradc[k] = (conc2 - conc1) / delta;
135  c1sum += conc1;
136  c2sum += conc2;
137  }
138 
139  // Calculate the pressures at p1 p2 and pbar
140  doublereal p1 = c1sum * GasConstant * t1;
141  doublereal p2 = c2sum * GasConstant * t2;
142  doublereal pbar = 0.5*(p1 + p2);
143  doublereal gradp = (p2 - p1)/delta;
144  doublereal tbar = 0.5*(t1 + t2);
145  m_thermo->setState_TPX(tbar, pbar, cbar);
147 
148  // Multiply m_multidiff and gradc together and store the result in fluxes[]
149  multiply(m_multidiff, gradc, fluxes);
150  divide_each(cbar, cbar + m_nsp, m_dk.begin());
151 
152  // if no permeability has been specified, use result for
153  // close-packed spheres
154  double b = 0.0;
155  if (m_perm < 0.0) {
156  double p = m_porosity;
157  double d = m_diam;
158  double t = m_tortuosity;
159  b = p*p*p*d*d/(72.0*t*(1.0-p)*(1.0-p));
160  } else {
161  b = m_perm;
162  }
163  b *= gradp / m_gastran->viscosity();
164  scale(cbar, cbar + m_nsp, cbar, b);
165 
166  // Multiply m_multidiff with cbar and add it to fluxes
167  increment(m_multidiff, cbar, fluxes);
168  scale(fluxes, fluxes + m_nsp, fluxes, -1.0);
169 }
170 
172 {
173  // see if temperature has changed
175 
176  // update the mole fractions
178  eval_H_matrix();
179 
180  // invert H
181  int ierr = invert(m_multidiff);
182  if (ierr != 0) {
183  throw CanteraError("DustyGasTransport::updateMultiDiffCoeffs",
184  "invert returned ierr = {}", ierr);
185  }
186 }
187 
188 void DustyGasTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
189 {
191  for (size_t i = 0; i < m_nsp; i++) {
192  for (size_t j = 0; j < m_nsp; j++) {
193  d[ld*j + i] = m_multidiff(i,j);
194  }
195  }
196 }
197 
199 {
200  if (m_temp == m_thermo->temperature()) {
201  return;
202  }
204  m_knudsen_ok = false;
205  m_bulk_ok = false;
206 }
207 
209 {
210  m_thermo->getMoleFractions(m_x.data());
211 
212  // add an offset to avoid a pure species condition
213  // (check - this may be unnecessary)
214  for (size_t k = 0; k < m_nsp; k++) {
215  m_x[k] = std::max(Tiny, m_x[k]);
216  }
217  // diffusion coeffs depend on Pressure
218  m_bulk_ok = false;
219 }
220 
221 void DustyGasTransport::setPorosity(doublereal porosity)
222 {
223  m_porosity = porosity;
224  m_knudsen_ok = false;
225  m_bulk_ok = false;
226 }
227 
228 void DustyGasTransport::setTortuosity(doublereal tort)
229 {
230  m_tortuosity = tort;
231  m_knudsen_ok = false;
232  m_bulk_ok = false;
233 }
234 
236 {
237  m_pore_radius = rbar;
238  m_knudsen_ok = false;
239 }
240 
242 {
243  m_diam = dbar;
244 }
245 
247 {
248  m_perm = B;
249 }
250 
252 {
253  return *m_gastran;
254 }
255 
256 }
doublereal m_pore_radius
Pore radius (meter)
const vector_fp & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:437
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::unique_ptr< Transport > m_gastran
Pointer to the transport object for the gas phase.
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:269
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...
doublereal temperature() const
Temperature (K).
Definition: Phase.h:601
virtual void setState_TPX(doublereal t, doublereal p, const doublereal *x)
Set the temperature (K), pressure (Pa), and mole fractions.
Definition: ThermoPhase.cpp:95
doublereal m_diam
Particle diameter.
vector_fp m_x
mole fractions
thermo_t * m_thermo
pointer to the object representing the phase
DenseMatrix m_d
binary diffusion coefficients
Base class for transport property managers.
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:266
STL namespace.
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:292
void updateMultiDiffCoeffs()
Update the Multicomponent diffusion coefficients that are used in the approximation.
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:93
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:69
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.
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:65
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)
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:466
doublereal m_porosity
Porosity.
thermo_t & thermo()
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:130
const doublereal Tiny
Small number to compare differences of mole fractions against.
Definition: ct_defs.h:143
const doublereal GasConstant
Universal Gas Constant. [J/kmol/K].
Definition: ct_defs.h:64
Contains declarations for string manipulation functions within Cantera.
size_t m_nsp
Number of species.
doublereal m_temp
temperature
virtual void setThermo(thermo_t &thermo)
Specifies the ThermoPhase object.
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:8
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.