initial import

Cantera/src/zeroD/FlowReactor.cpp

@@ -0,0 +1,146 @@
+ /**
+ *  @file ReactorZND.cpp
+ *
+ *  A zero-dimensional reactor
+ */
+ 
+// Copyright 2001  California Institute of Technology
+
+
+#ifdef WIN32
+#pragma warning(disable:4786)
+#pragma warning(disable:4503)
+#endif
+
+#include "FlowReactor.h"
+
+using namespace Cantera;
+
+namespace CanteraZeroD {
+
+
+    FlowReactor::FlowReactor() : Reactor(), m_fctr(1.0e10), 
+                                 m_speed0(0.0) {}
+
+    // overloaded method of FuncEval. Called by the integrator to
+    // get the initial conditions.
+    void FlowReactor::getInitialConditions(double t0, size_t leny, double* y) 
+    {
+        m_init = true;
+        if (m_mix == 0) {
+            writelog("Error: reactor is empty.\n");
+            return;
+        }  
+        m_time = t0;
+        m_mix->restoreState(m_state);
+
+        m_mix->getMassFractions(y+2);
+            
+        y[0] = 0.0; // distance
+        
+        // set the second component to the initial speed
+        y[1] = m_speed0;
+    }
+
+    /*
+     *  Must be called before calling method 'advance'
+     */
+    void FlowReactor::initialize(doublereal t0) {
+        m_mix->restoreState(m_state);
+        m_nv = m_nsp + 2;
+        m_init = true;
+    }
+    
+    void FlowReactor::updateState(doublereal* y) {
+
+        // Set the mass fractions and  density of the mixture.
+        m_dist         = y[0];
+        m_speed        = y[1];
+        doublereal* mss = y + 2;
+        //        doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0);
+        m_mix->setMassFractions(mss);
+        
+        doublereal rho = m_rho0 * m_speed0/m_speed;
+
+        // assumes frictionless
+        doublereal pmom = m_P0 - rho*m_speed*m_speed;
+
+        doublereal hmom;
+        // assumes adiabatic
+        if (m_energy) {
+            hmom = m_h0 - 0.5*m_speed*m_speed;
+            m_thermo->setState_HP(hmom, pmom);
+        }
+        else {
+            m_thermo->setState_TP(m_T, pmom);
+        }
+        m_mix->saveState(m_state);
+    }
+
+
+    /*
+     * Called by the integrator to evaluate ydot given y at time 'time'.
+     */
+    void FlowReactor::evalEqs(doublereal time, doublereal* y, 
+        doublereal* ydot, doublereal* params) 
+    {
+        m_time = time;
+        m_mix->restoreState(m_state);
+
+        double mult;
+        Kinetics* kin;
+        int n, npar;
+
+        // process sensitivity parameters
+        if (params) {
+            npar = nSensParams();
+            for (n = 0; n < npar; n++) {
+                mult = m_kin->multiplier(m_pnum[n]);
+                m_kin->setMultiplier(m_pnum[n], mult*params[n]);
+            }
+        }
+
+        // distance equation
+        ydot[0] = m_speed;
+
+        // speed equation. Set m_fctr to a large value, so that rho*u is
+        // held fixed
+        ydot[1] = m_fctr*(m_speed0 - m_mix->density()*m_speed/m_rho0);
+
+        /* species equations */
+        const doublereal* mw = m_mix->molecularWeights().begin();
+
+        if (m_chem) {
+            m_kin->getNetProductionRates(ydot+2);   // "omega dot"
+        }
+        else {
+            fill(ydot + 2, ydot + 2 + m_nsp, 0.0);
+        }
+        doublereal rrho = 1.0/m_mix->density();
+        for (n = 0; n < m_nsp; n++) {
+            ydot[n+2] *= mw[n]*rrho;
+        }
+
+        // reset sensitivity parameters
+        if (params) {
+            npar = nSensParams();
+            for (n = 0; n < npar; n++) {
+                mult = m_kin->multiplier(m_pnum[n]);
+                m_kin->setMultiplier(m_pnum[n], mult/params[n]);
+            }
+        }
+
+        
+    }
+
+
+    int FlowReactor::componentIndex(string nm) const {
+        if (nm == "X") return 0;
+        if (nm == "U") return 1;
+        // check for a gas species name
+        int k = m_mix->speciesIndex(nm);
+        if (k >= 0) return k + 2;
+        else return -1;
+    }
+
+}

Cantera/src/zeroD/FlowReactor.h

@@ -0,0 +1,84 @@
+/**
+ *  @file FlowReactor.h
+ *
+ * $Author$
+ * $Revision$
+ * $Date$
+ */
+
+// Copyright 2001  California Institute of Technology
+
+#ifndef CT_FLOWREACTOR_H
+#define CT_FLOWREACTOR_H
+
+#ifdef WIN32
+#pragma warning(disable:4786)
+#pragma warning(disable:4503)
+#endif
+
+#include "Reactor.h"
+
+namespace CanteraZeroD {
+
+    /**
+     */
+
+    class FlowReactor : public Reactor {
+
+    public:
+
+        /**
+         * Default constructor.
+         */
+        FlowReactor();
+
+        /**
+         * Destructor. 
+         */
+        virtual ~FlowReactor(){}
+        
+        virtual int type() const { return FlowReactorType; }
+
+        //-----------------------------------------------------
+
+        virtual void getInitialConditions(doublereal t0, size_t leny, 
+            doublereal* y);
+
+
+        //-----------------------------------------------------
+        virtual int neq() { return m_nv; }
+
+        virtual void initialize(doublereal t0 = 0.0);
+	virtual void evalEqs(doublereal t, doublereal* y, 
+            doublereal* ydot, doublereal* params);
+        virtual void updateState(doublereal* y);
+
+        void setMassFlowRate(doublereal mdot) {
+            m_rho0 = m_thermo->density();
+            m_speed = mdot/m_rho0;
+            m_speed0 = m_speed;
+            m_T = m_thermo->temperature();
+            m_P0 = m_thermo->pressure() + m_rho0*m_speed*m_speed;
+            m_h0 = m_thermo->enthalpy_mass() + 0.5*m_speed*m_speed;
+        }
+
+        void setTimeConstant(doublereal tau) {
+            m_fctr = 1.0/tau;
+        }
+
+        double speed() const { return m_speed; }
+        double distance() const { return m_dist; }
+        virtual int componentIndex(string nm) const;
+
+    protected:
+        
+        doublereal m_speed, m_dist, m_T;
+        doublereal m_fctr;
+        doublereal m_rho0, m_speed0, m_P0, m_h0;
+
+    private:
+    };
+}
+
+#endif
+