Cantera  3.1.0a1
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 https://cantera.org/license.txt for license and copyright information.
8 
12 #include "cantera/base/utilities.h"
13 
14 namespace Cantera
15 {
16 
18 {
19  Transport::init(phase);
20  // constant mixture attributes
21  m_thermo = phase;
22  m_nsp = m_thermo->nSpecies();
23  if (m_gastran.get() != gastr) {
24  m_gastran.reset(gastr);
25  }
26 
27  // make a local copy of the molecular weights
29 
32  m_dk.resize(m_nsp, 0.0);
33 
34  m_x.resize(m_nsp, 0.0);
35  m_thermo->getMoleFractions(m_x.data());
36 
37  // set flags all false
38  m_knudsen_ok = false;
39  m_bulk_ok = false;
40 
41  m_spwork.resize(m_nsp);
42  m_spwork2.resize(m_nsp);
43 }
44 
46 {
47  if (m_bulk_ok) {
48  return;
49  }
50 
51  // get the gaseous binary diffusion coefficients
52  m_gastran->getBinaryDiffCoeffs(m_nsp, m_d.ptrColumn(0));
53  double por2tort = m_porosity / m_tortuosity;
54  for (size_t n = 0; n < m_nsp; n++) {
55  for (size_t m = 0; m < m_nsp; m++) {
56  m_d(n,m) *= por2tort;
57  }
58  }
59  m_bulk_ok = true;
60 }
61 
63 {
64  if (m_knudsen_ok) {
65  return;
66  }
67  double K_g = m_pore_radius * m_porosity / m_tortuosity;
68  for (size_t k = 0; k < m_nsp; k++) {
69  m_dk[k] = 2.0/3.0 * K_g * sqrt((8.0 * GasConstant * m_temp)/
70  (Pi * m_mw[k]));
71  }
72  m_knudsen_ok = true;
73 }
74 
76 {
79  for (size_t k = 0; k < m_nsp; k++) {
80  // evaluate off-diagonal terms
81  for (size_t j = 0; j < m_nsp; j++) {
82  m_multidiff(k,j) = -m_x[k]/m_d(k,j);
83  }
84 
85  // evaluate diagonal term
86  double sum = 0.0;
87  for (size_t j = 0; j < m_nsp; j++) {
88  if (j != k) {
89  sum += m_x[j]/m_d(k,j);
90  }
91  }
92  m_multidiff(k,k) = 1.0/m_dk[k] + sum;
93  }
94 }
95 
96 void DustyGasTransport::getMolarFluxes(const double* const state1,
97  const double* const state2,
98  const double delta,
99  double* const fluxes)
100 {
101  // cbar will be the average concentration between the two points
102  double* const cbar = m_spwork.data();
103  double* const gradc = m_spwork2.data();
104  const double t1 = state1[0];
105  const double t2 = state2[0];
106  const double rho1 = state1[1];
107  const double rho2 = state2[1];
108  const double* const y1 = state1 + 2;
109  const double* const y2 = state2 + 2;
110  double c1sum = 0.0, c2sum = 0.0;
111 
112  for (size_t k = 0; k < m_nsp; k++) {
113  double conc1 = rho1 * y1[k] / m_mw[k];
114  double conc2 = rho2 * y2[k] / m_mw[k];
115  cbar[k] = 0.5*(conc1 + conc2);
116  gradc[k] = (conc2 - conc1) / delta;
117  c1sum += conc1;
118  c2sum += conc2;
119  }
120 
121  // Calculate the pressures at p1 p2 and pbar
122  double p1 = c1sum * GasConstant * t1;
123  double p2 = c2sum * GasConstant * t2;
124  double pbar = 0.5*(p1 + p2);
125  double gradp = (p2 - p1)/delta;
126  double tbar = 0.5*(t1 + t2);
127  m_thermo->setState_TPX(tbar, pbar, cbar);
129 
130  // Multiply m_multidiff and gradc together and store the result in fluxes[]
131  multiply(m_multidiff, gradc, fluxes);
132  for (size_t k = 0; k < m_nsp; k++) {
133  cbar[k] /= m_dk[k];
134  }
135 
136  // if no permeability has been specified, use result for
137  // close-packed spheres
138  double b = 0.0;
139  if (m_perm < 0.0) {
140  double p = m_porosity;
141  double d = m_diam;
142  double t = m_tortuosity;
143  b = p*p*p*d*d/(72.0*t*(1.0-p)*(1.0-p));
144  } else {
145  b = m_perm;
146  }
147  b *= gradp / m_gastran->viscosity();
148  scale(cbar, cbar + m_nsp, cbar, b);
149 
150  // Multiply m_multidiff with cbar and add it to fluxes
151  increment(m_multidiff, cbar, fluxes);
152  scale(fluxes, fluxes + m_nsp, fluxes, -1.0);
153 }
154 
156 {
157  // see if temperature has changed
159 
160  // update the mole fractions
162  eval_H_matrix();
163 
164  // invert H
165  int ierr = invert(m_multidiff);
166  if (ierr != 0) {
167  throw CanteraError("DustyGasTransport::updateMultiDiffCoeffs",
168  "invert returned ierr = {}", ierr);
169  }
170 }
171 
172 void DustyGasTransport::getMultiDiffCoeffs(const size_t ld, double* const d)
173 {
175  for (size_t i = 0; i < m_nsp; i++) {
176  for (size_t j = 0; j < m_nsp; j++) {
177  d[ld*j + i] = m_multidiff(i,j);
178  }
179  }
180 }
181 
183 {
184  if (m_temp == m_thermo->temperature()) {
185  return;
186  }
188  m_knudsen_ok = false;
189  m_bulk_ok = false;
190 }
191 
193 {
194  m_thermo->getMoleFractions(m_x.data());
195 
196  // add an offset to avoid a pure species condition
197  // (check - this may be unnecessary)
198  for (size_t k = 0; k < m_nsp; k++) {
199  m_x[k] = std::max(Tiny, m_x[k]);
200  }
201  // diffusion coeffs depend on Pressure
202  m_bulk_ok = false;
203 }
204 
205 void DustyGasTransport::setPorosity(double porosity)
206 {
207  m_porosity = porosity;
208  m_knudsen_ok = false;
209  m_bulk_ok = false;
210 }
211 
213 {
214  m_tortuosity = tort;
215  m_knudsen_ok = false;
216  m_bulk_ok = false;
217 }
218 
220 {
221  m_pore_radius = rbar;
222  m_knudsen_ok = false;
223 }
224 
226 {
227  m_diam = dbar;
228 }
229 
231 {
232  m_perm = B;
233 }
234 
236 {
237  return *m_gastran;
238 }
239 
240 }
Headers for the DustyGasTransport object, which models transport properties in porous media using the...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
double * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:203
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:66
void resize(size_t n, size_t m, double v=0.0) override
Resize the matrix.
Definition: DenseMatrix.cpp:60
vector< double > m_mw
Local copy of the species molecular weights.
double m_diam
Particle diameter.
void getMolarFluxes(const double *const state1, const double *const state2, const double delta, double *const fluxes) override
Get the molar fluxes [kmol/m^2/s], given the thermodynamic state at two nearby points.
DenseMatrix m_multidiff
Multicomponent diffusion coefficients.
bool m_bulk_ok
Update-to-date variable for Binary diffusion coefficients.
vector< double > m_x
mole fractions
vector< double > m_dk
Knudsen diffusion coefficients.
void updateTransport_T()
Update temperature-dependent quantities within the object.
bool m_knudsen_ok
Update-to-date variable for Knudsen diffusion coefficients.
void initialize(ThermoPhase *phase, Transport *gastr)
Initialization routine called by TransportFactory.
vector< double > m_spwork
work space of size m_nsp;
void updateBinaryDiffCoeffs()
Private routine to update the dusty gas binary diffusion coefficients.
double m_perm
Permeability of the media.
double m_tortuosity
Tortuosity.
void eval_H_matrix()
Calculate the H matrix.
void updateTransport_C()
Update concentration-dependent quantities within the object.
void updateMultiDiffCoeffs()
Update the Multicomponent diffusion coefficients that are used in the approximation.
void updateKnudsenDiffCoeffs()
Update the Knudsen diffusion coefficients.
void setMeanParticleDiameter(double dbar)
Set the mean particle diameter.
void setTortuosity(double tort)
Set the tortuosity (dimensionless)
void setPorosity(double porosity)
Set the porosity (dimensionless)
vector< double > m_spwork2
work space of size m_nsp;
Transport & gasTransport()
Return a reference to the transport manager used to compute the gas binary diffusion coefficients and...
unique_ptr< Transport > m_gastran
Pointer to the transport object for the gas phase.
void setPermeability(double B)
Set the permeability of the media.
double m_pore_radius
Pore radius (meter)
void getMultiDiffCoeffs(const size_t ld, double *const d) override
Return the Multicomponent diffusion coefficients. Units: [m^2/s].
void setMeanPoreRadius(double rbar)
Set the mean pore radius (m)
DenseMatrix m_d
binary diffusion coefficients
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:231
double temperature() const
Temperature (K).
Definition: Phase.h:562
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:434
const vector< double > & molecularWeights() const
Return a const reference to the internal vector of molecular weights.
Definition: Phase.cpp:395
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:390
virtual void setState_TPX(double t, double p, const double *x)
Set the temperature (K), pressure (Pa), and mole fractions.
Definition: ThermoPhase.cpp:85
Base class for transport property managers.
Definition: Transport.h:72
ThermoPhase * m_thermo
pointer to the object representing the phase
Definition: Transport.h:420
virtual void init(ThermoPhase *thermo, int mode=0, int log_level=0)
Initialize a transport manager.
Definition: Transport.h:407
size_t m_nsp
Number of species.
Definition: Transport.h:423
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
Definition: utilities.h:104
const double GasConstant
Universal Gas Constant [J/kmol/K].
Definition: ct_defs.h:120
const double Pi
Pi.
Definition: ct_defs.h:68
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:564
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.
const double Tiny
Small number to compare differences of mole fractions against.
Definition: ct_defs.h:173
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.
int invert(DenseMatrix &A, size_t nn)
invert A. A is overwritten with A^-1.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...