added DustyGasTransport

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 <iostream>
+
+/** 
+ * 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]);
+        }
+    }
+}

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 <vector>
+#include <string>
+#include <map>
+#include <numeric>
+#include <algorithm>
+
+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
+
+
+
+
+
+

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

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;
 

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");
         }