From 3af1787662a2a0b8767624334918f239b8511b5c Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Thu, 25 Sep 2003 20:37:58 +0000 Subject: [PATCH] added DustyGasTransport --- Cantera/src/transport/DustyGasTransport.cpp | 214 ++++++++++++++++++++ Cantera/src/transport/DustyGasTransport.h | 165 +++++++++++++++ Cantera/src/transport/Makefile.in | 2 +- Cantera/src/transport/TransportBase.h | 1 + Cantera/src/transport/TransportFactory.cpp | 12 +- 5 files changed, 391 insertions(+), 3 deletions(-) create mode 100644 Cantera/src/transport/DustyGasTransport.cpp create mode 100644 Cantera/src/transport/DustyGasTransport.h diff --git a/Cantera/src/transport/DustyGasTransport.cpp b/Cantera/src/transport/DustyGasTransport.cpp new file mode 100644 index 000000000..34bb33808 --- /dev/null +++ b/Cantera/src/transport/DustyGasTransport.cpp @@ -0,0 +1,214 @@ +/** + * + * @file DustyGasTransport.cpp + * Implementation file for class DustyGasTransport + * + * @ingroup transportProps + * + * $Author$ + * $Date$ + * $Revision$ + * + * Copyright 2003 California Institute of Technology + * See file License.txt for licensing information + * + */ + + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "DustyGasTransport.h" +#include "ctlapack.h" +#include "../../ext/math/gmres.h" + +#include "DenseMatrix.h" +#include "polyfit.h" +#include "utilities.h" +#include "TransportParams.h" +#include "IdealGasPhase.h" + +#include "TransportFactory.h" + +#include + +/** + * Mole fractions below MIN_X will be set to MIN_X when computing + * transport properties. + */ + +#define MIN_X 1.e-20 + + +namespace Cantera { + + + //////////////////// class DustyGasTransport methods ////////////// + + + DustyGasTransport::DustyGasTransport(thermo_t* thermo) + : Transport(thermo), + m_temp(-1.0), + m_porosity(0.0), + m_tortuosity(1.0), + m_pore_radius(0.0), + m_diam(0.0), + m_perm(-1.0), + m_gastran(0) + {} + + + void DustyGasTransport::initialize(ThermoPhase* phase, Transport* gastr) { + + // constant mixture attributes + m_thermo = phase; + m_nsp = m_thermo->nSpecies(); + m_tmin = m_thermo->minTemp(); + m_tmax = m_thermo->maxTemp(); + + m_gastran = gastr; + + // make a local copy of the molecular weights + m_mw.resize(m_nsp); + copy(m_thermo->molecularWeights().begin(), + m_thermo->molecularWeights().end(), m_mw.begin()); + + m_multidiff.resize(m_nsp, m_nsp); + m_d.resize(m_nsp, m_nsp); + m_dk.resize(m_nsp, 0.0); + m_x.resize(m_nsp); + + // set flags all false + m_knudsen_ok = false; + m_bulk_ok = false; + + // some work space + m_spwork.resize(m_nsp); + } + + + /******************* binary diffusion coefficients **************/ + + + void DustyGasTransport::updateBinaryDiffCoeffs() { + if (m_bulk_ok) return; + int n,m; + // get the gaseous binary diffusion coefficients + m_gastran->getBinaryDiffCoeffs(m_nsp, m_d.begin()); + doublereal por2tort = m_porosity / m_tortuosity; + for (n = 0; n < m_nsp; n++) + for (m = 0; m < m_nsp; m++) + m_d(n,m) *= por2tort; + m_bulk_ok = true; + } + + void DustyGasTransport::updateKnudsenDiffCoeffs() { + if (m_knudsen_ok) return; + doublereal K_g = m_pore_radius * m_porosity / m_tortuosity; + const doublereal FourThirds = 4.0/3.0; + for (int k = 0; k < m_nsp; k++) { + m_dk[k] = FourThirds * K_g * sqrt((8.0 * GasConstant * m_temp)/ + (Pi * m_mw[k])); + } + m_knudsen_ok = true; + } + + + void DustyGasTransport::eval_H_matrix() { + updateBinaryDiffCoeffs(); + updateKnudsenDiffCoeffs(); + int k,l,j; + doublereal sum; + for (k = 0; k < m_nsp; k++) { + + // evaluate off-diagonal terms + for (l = 0; l < m_nsp; l++) m_multidiff(k,l) = -m_x[k]/m_d(k,l); + + // evaluate diagonal term + sum = 0.0; + for (j = 0; j < m_nsp; j++) sum += m_x[j]/m_d(k,j); + m_multidiff(k,k) = 1.0/m_dk[k] + sum; + } + } + + void DustyGasTransport::getMolarFluxes(const double* grad_conc, + double grad_P, double* fluxes) { + updateMultiDiffCoeffs(); + copy(grad_conc, grad_conc + m_nsp, m_spwork.begin()); + multiply(m_multidiff, m_spwork.begin(), fluxes); + m_thermo->getConcentrations(m_spwork.begin()); + divide_each(m_spwork.begin(), m_spwork.end(), m_dk.begin()); + + // if no permeability has been specified, use result for + // close-packed spheres + double b = 0.0; + if (m_perm < 0.0) { + double p = m_porosity; + double d = m_diam; + double t = m_tortuosity; + b = p*p*p*d*d/(72.0*t*(1.0-p)*(1.0-p)); + } + else { + b = m_perm; + } + b *= grad_P / m_gastran->viscosity(); + scale(m_spwork.begin(), m_spwork.end(), m_spwork.begin(), b); + increment(m_multidiff, m_spwork.begin(), fluxes); + scale(fluxes, fluxes + m_nsp, fluxes, -1.0); + } + + void DustyGasTransport::updateMultiDiffCoeffs() { + + // see if temperature has changed + updateTransport_T(); + + // update the mole fractions + updateTransport_C(); + + eval_H_matrix(); + + // invert H + int ierr = invert(m_multidiff, m_nsp); + if (ierr != 0) { + throw CanteraError("DustyGasTransport::updateMultiDiffCoeffs", + "invert returned ierr = "+int2str(ierr)); + } + } + + void DustyGasTransport::getMultiDiffCoeffs(int ld, doublereal* d) { + int i,j; + updateMultiDiffCoeffs(); + for (i = 0; i < m_nsp; i++) { + for (j = 0; j < m_nsp; j++) { + d[ld*j + i] = m_multidiff(i,j); + } + } + } + + + /** + * Update temperature-dependent quantities. + */ + void DustyGasTransport::updateTransport_T() + { + if (m_temp == m_thermo->temperature()) return; + m_temp = m_thermo->temperature(); + m_knudsen_ok = false; + m_bulk_ok = false; + } + + void DustyGasTransport::updateTransport_C() + { + m_thermo->getMoleFractions(m_x.begin()); + + // add an offset to avoid a pure species condition + // (check - this may be unnecessary) + int k; + for (k = 0; k < m_nsp; k++) { + m_x[k] = fmaxx(MIN_X, m_x[k]); + } + } +} diff --git a/Cantera/src/transport/DustyGasTransport.h b/Cantera/src/transport/DustyGasTransport.h new file mode 100644 index 000000000..0062a227c --- /dev/null +++ b/Cantera/src/transport/DustyGasTransport.h @@ -0,0 +1,165 @@ +/** + * + * @file DustyGasTransport.h + * Interface for class DustyGasTransport + * + */ + +// Copyright 2003 California Institute of Technology + + +#ifndef CT_DUSTYGASTRAN_H +#define CT_DUSTYGASTRAN_H + + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + + +// STL includes +#include +#include +#include +#include +#include + +using namespace std; + +// Cantera includes +#include "TransportBase.h" +#include "../DenseMatrix.h" + + +namespace Cantera { + + + class TransportParams; + + class DustyGasTransport : public Transport { + + public: + + + virtual ~DustyGasTransport() {} + + // overloaded base class methods + virtual int model() { return cDustyGasTransport; } + + virtual void setParameters(int type, int k, doublereal* p) { + switch(type) { + case 0: + setPorosity(p[0]); break; + case 1: + setTortuosity(p[0]); break; + case 2: + setMeanPoreRadius(p[0]); break; + case 3: + setMeanParticleDiameter(p[0]); break; + case 4: + setPermeability(p[0]); break; + default: + throw CanteraError("DustyGasTransport::init", + "unknown parameter"); + } + } + + virtual void getBinaryDiffCoeffs(int ld, doublereal* d); + virtual void getMultiDiffCoeffs(int ld, doublereal* d); + + + // new methods + + void getMolarFluxes(const double* grad_conc, + double grad_P, double* fluxes); + + void setPorosity(doublereal porosity) { + m_porosity = porosity; + m_knudsen_ok = false; + m_bulk_ok = false; + } + + void setTortuosity(doublereal tort) { + m_tortuosity = tort; + m_knudsen_ok = false; + m_bulk_ok = false; + } + + void setMeanPoreRadius(doublereal rbar) { + m_pore_radius = rbar; + m_knudsen_ok = false; + } + + void setMeanParticleDiameter(doublereal dbar) { + m_diam = dbar; + } + + void setPermeability(doublereal B) { + m_perm = B; + } + + /** + * @internal + */ + virtual bool init(TransportParams& tr); + + void updateTransport_T(); + void updateTransport_C(); + + friend class TransportFactory; + + protected: + + void updateBinaryDiffCoeffs(); + void updateMultiDiffCoeffs(); + void updateKnudsenDiffCoeffs(); + void eval_H_matrix(); + /// default constructor + DustyGasTransport(thermo_t* thermo=0); + void initialize(ThermoPhase* phase, Transport* gastr); + private: + + + // mixture attributes + int m_nsp; + doublereal m_tmin, m_tmax; + vector_fp m_mw; + + // property values + DenseMatrix m_d; + vector_fp m_visc; + vector_fp m_x; + vector_fp m_dk; + doublereal m_temp; + + // H matrix quantities + DenseMatrix m_multidiff; + + // work space + vector_fp m_spwork; + + bool m_knudsen_ok; + bool m_bulk_ok; + + doublereal m_porosity; + doublereal m_tortuosity; + doublereal m_pore_radius; + doublereal m_diam; + doublereal m_perm; + + Transport* m_gastran; + + doublereal pressure_ig() { + return m_thermo->molarDensity() * GasConstant * m_thermo->temperature(); + } + }; +} +#endif + + + + + + diff --git a/Cantera/src/transport/Makefile.in b/Cantera/src/transport/Makefile.in index 04bdeeda4..e9119edd6 100644 --- a/Cantera/src/transport/Makefile.in +++ b/Cantera/src/transport/Makefile.in @@ -17,7 +17,7 @@ CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) # stirred reactors OBJS = TransportFactory.o MultiTransport.o MixTransport.o MMCollisionInt.o \ - SolidTransport.o + SolidTransport.o DustyGasTransport.o CXX_INCLUDES = -I.. LIB = @buildlib@/libtransport.a diff --git a/Cantera/src/transport/TransportBase.h b/Cantera/src/transport/TransportBase.h index f898b3916..a8a8998ec 100755 --- a/Cantera/src/transport/TransportBase.h +++ b/Cantera/src/transport/TransportBase.h @@ -47,6 +47,7 @@ namespace Cantera { const int cMixtureAveraged = 210; const int CK_MixtureAveraged = 211; const int cSolidTransport = 300; + const int cDustyGasTransport = 400; class XML_Writer; diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp index 94825a891..00cb121a7 100755 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -18,7 +18,7 @@ #include "MultiTransport.h" #include "MixTransport.h" #include "SolidTransport.h" - +#include "DustyGasTransport.h" #include "TransportFactory.h" #include "polyfit.h" @@ -196,6 +196,7 @@ namespace Cantera { m_models["Mix"] = cMixtureAveraged; m_models["Multi"] = cMulticomponent; m_models["Solid"] = cSolidTransport; + m_models["DustyGas"] = cDustyGasTransport; m_models["None"] = 0; } @@ -216,7 +217,8 @@ namespace Cantera { if (transportModel == "") return new Transport; vector_fp state; - Transport* tr = 0; + Transport *tr = 0, *gastr = 0; + DustyGasTransport* dtr = 0; phase->saveState(state); switch(m_models[transportModel]) { @@ -242,6 +244,12 @@ namespace Cantera { tr = new SolidTransport; tr->setThermo(*phase); break; + case cDustyGasTransport: + tr = new DustyGasTransport; + gastr = new MixTransport; + dtr = (DustyGasTransport*)tr; + dtr->initialize(phase, gastr); + break; default: throw CanteraError("newTransport","unknown transport model"); }