WIP: 1d point control

include/cantera/oneD/Sim1D.h

@@ -252,6 +252,14 @@ class Sim1D : public OneDim
     //! Return location of the point where temperature is fixed
     double fixedTemperatureLocation();
 
+    //! Set the fuel internal boundary location. This is used for one/two-point
+    //! flame control for stagnation flows.
+    void setFuelInternalBoundary(double temperature);
+
+    //! Set the oxidizer internal boundary location. This is used for one/two-point
+    //! flame control for stagnation flows.
+    void setOxidizerInternalBoundary(double temperature);
+
     /**
      * Set grid refinement criteria. If dom >= 0, then the settings
      * apply only to the specified domain.  If dom < 0, the settings

include/cantera/oneD/StFlow.h

@@ -27,7 +27,8 @@ enum offset
     , c_offset_T //! temperature
     , c_offset_L //! (1/r)dP/dr
     , c_offset_E //! electric poisson's equation
-    , c_offset_Y //! mass fractions
+    , c_offset_Uo //! oxidizer axial velocity
+    , c_offset_Y //! mass fractions, additions after this can cause issues, add them before
 };
 
 class Transport;
@@ -155,6 +156,35 @@ class StFlow : public Domain1D
         return m_fixedtemp[j];
     }
 
+    //! The current fuel internal boundary temperature
+    double fuelInternalBoundaryTemperature() const {
+        if ((m_onePointControl || m_twoPointControl) && (m_zFuel != Undef)) {
+            return m_tFuel;
+        }
+    }
+
+    //! Set the temperature in the fuel side internal boundary
+    void setFuelInternalBoundaryTemperature(double tFuel) {
+        if ((m_onePointControl || m_twoPointControl) && (m_zFuel != Undef)) {
+            m_tFuel = tFuel;
+        }
+    }
+
+    //! The current oxidizer side internal boundary temperature
+    double oxidInternalBoundaryTemperature() const {
+        if (m_twoPointControl && (m_zOxid != Undef)) {
+            return m_tOxid;
+        }
+    }
+
+    //! Set the temperature in the oxidizer side internal boundary
+    void setOxidInternalBoundaryTemperature(doublereal tOxid) {
+        if (m_twoPointControl && (m_zOxid != Undef)) {
+            m_tOxid = tOxid;
+        }
+    }
+
+
     //! @}
 
     virtual std::string componentName(size_t n) const;
@@ -298,12 +328,20 @@ class StFlow : public Domain1D
 
     //! Evaluate all residual components at the right boundary.
     virtual void evalRightBoundary(double* x, double* res, int* diag,
-                                   double rdt);
+                                   double rdt, size_t j);
+
+    //! Evaluate all residual components at the left boundary.
+    virtual void evalLeftBoundary(double* x, double* res, int* diag,
+                                   double rdt, size_t j);
+
+    //! Evaluate all residual components at the left boundary.
+    virtual void evalInterior(double* x, double* res, int* diag,
+                                   double rdt, size_t j);
 
     //! Evaluate the residual corresponding to the continuity equation at all
     //! interior grid points.
-    virtual void evalContinuity(size_t j, double* x, double* r,
-                                int* diag, double rdt);
+    //virtual void evalContinuity(size_t j, double* x, double* r,
+    //                            int* diag, double rdt);
 
     //! Index of the species on the left boundary with the largest mass fraction
     size_t leftExcessSpecies() const {
@@ -315,6 +353,36 @@ class StFlow : public Domain1D
         return m_kExcessRight;
     }
 
+    //! Set the status of one-point flame control
+    void enableOnePointControl(bool onePointControl) {
+        if (m_twoPointControl && onePointControl) {
+            throw CanteraError("StFlow::enableOnePointControl",
+                               "Two different point control techniques can not be set simultaneously");
+        } else {
+            m_onePointControl = onePointControl;
+        }
+    }
+
+    //! Set the status of two-point flame control
+    void enableTwoPointControl(bool twoPointControl) {
+        if (m_onePointControl && twoPointControl) {
+            throw CanteraError("StFlow::enableTwoPointControl",
+                               "Two different point control techniques can not be set simultaneously");
+        } else {
+            m_twoPointControl = twoPointControl;
+        }
+    }
+
+    //! Get the status of one point flame control
+    bool onePointControlEnabled() {
+        return m_onePointControl;
+    }
+
+    //! Get the status of two point flame control
+    bool twoPointControlEnabled() {
+        return m_twoPointControl;
+    }
+
 protected:
     virtual AnyMap getMeta() const;
     virtual void setMeta(const AnyMap& state);
@@ -338,6 +406,13 @@ class StFlow : public Domain1D
     //! after updateProperties is called.
     virtual void evalResidual(double* x, double* rsd, int* diag,
                               double rdt, size_t jmin, size_t jmax);
+    void computeRadiation(double* x, size_t jmin, size_t jmax);
+    void evalEnergyEquation(double* x, double* rsd, int* diag, double rdt, size_t j);
+    void evalSpeciesEquation(double* x, double* rsd, int* diag, double rdt, size_t j);
+    void evalContinuityEquation(double* x, double* rsd, int* diag, double rdt, size_t j);
+    void evalRadialMomentumEquation(double* x, double* rsd, int* diag, double rdt, size_t j);
+    void evalLEquation(double* x, double* rsd, int* diag, double rdt, size_t j);
+    void evalUoEquation(double* x, double* rsd, int* diag, double rdt, size_t j);
 
     /**
      * Update the thermodynamic properties from point j0 to point j1
@@ -385,6 +460,10 @@ class StFlow : public Domain1D
         return x[index(c_offset_L, j)];
     }
 
+    double Uo(const double* x, size_t j) const {
+        return x[index(c_offset_Uo, j)];
+    }
+
     doublereal Y(const doublereal* x, size_t k, size_t j) const {
         return x[index(c_offset_Y + k, j)];
     }
@@ -527,13 +606,30 @@ class StFlow : public Domain1D
     //! to `j1`, based on solution `x`.
     virtual void updateTransport(doublereal* x, size_t j0, size_t j1);
 
+    //! Flag for one-point  & two-point flame control
+    bool m_onePointControl = false;
+    bool m_twoPointControl = false;
+
 public:
     //! Location of the point where temperature is fixed
     double m_zfixed = Undef;
 
     //! Temperature at the point used to fix the flame location
     double m_tfixed = -1.0;
 
+
+    //! The location of the internal boundary at fuel side
+    double m_zFuel = Undef;
+
+    //! The fixed temperature at the fuel side internal boundary
+    double m_tFuel = Undef;
+
+    //! The location of the internal boundary at oxidizer side
+    double m_zOxid = Undef;
+
+    //! The fixed temperature at the oxidizer side internal boundary
+    double m_tOxid = Undef;
+
 private:
     vector_fp m_ybar;
 };

interfaces/cython/cantera/_onedim.pxd

@@ -74,6 +74,10 @@ cdef extern from "cantera/oneD/StFlow.h":
         void setPressure(double)
         void enableRadiation(cbool)
         cbool radiationEnabled()
+        void enableOnePointControl(cbool)
+        cbool onePointControlEnabled()
+        void enableTwoPointControl(cbool)
+        cbool twoPointControlEnabled()
         double radiativeHeatLoss(size_t)
         double pressure()
         void setFixedTempProfile(vector[double]&, vector[double]&)
@@ -87,6 +91,10 @@ cdef extern from "cantera/oneD/StFlow.h":
         double rightEmissivity()
         void solveEnergyEqn()
         void fixTemperature()
+        double fuelInternalBoundaryTemperature()
+        double oxidInternalBoundaryTemperature()
+        void setFuelInternalBoundaryTemperature(double)
+        void setOxidInternalBoundaryTemperature(double)
         cbool doEnergy(size_t)
         void enableSoret(cbool) except +translate_exception
         cbool withSoret()
@@ -141,6 +149,8 @@ cdef extern from "cantera/oneD/Sim1D.h":
         void setFixedTemperature(double) except +translate_exception
         double fixedTemperature()
         double fixedTemperatureLocation()
+        void setFuelInternalBoundary(double) except +translate_exception
+        void setOxidizerInternalBoundary(double) except +translate_exception
         void setInterrupt(CxxFunc1*) except +translate_exception
         void setTimeStepCallback(CxxFunc1*)
         void setSteadyCallback(CxxFunc1*)

interfaces/cython/cantera/_onedim.pyx

@@ -515,6 +515,20 @@ cdef class _FlowBase(Domain1D):
         def __set__(self, P):
             self.flow.setPressure(P)
 
+    property tFuel:
+        """ Fuel side internal boundary temperature [K] """
+        def __get__(self):
+            return self.flow.fuelInternalBoundaryTemperature()
+        def __set__(self, T):
+            self.flow.setFuelInternalBoundaryTemperature(T)
+
+    property tOxid:
+        """ Oxidizer side internal boundary temperature [K] """
+        def __get__(self):
+            return self.flow.oxidInternalBoundaryTemperature()
+        def __set__(self, T):
+            self.flow.setOxidInternalBoundaryTemperature(T)
+
     property transport_model:
         """
         Get/set the transport model used for calculating transport properties.
@@ -582,6 +596,20 @@ cdef class _FlowBase(Domain1D):
             else:
                 self.flow.fixTemperature()
 
+    property onePointControl_enabled:
+        """ Determines whether or not to enable one point flame control"""
+        def __get__(self):
+            return self.flow.onePointControlEnabled()
+        def __set__(self, enable):
+            self.flow.enableOnePointControl(<cbool>enable)
+
+    property twoPointControl_enabled:
+        """ Determines whether or not to enable two point flame control"""
+        def __get__(self):
+            return self.flow.twoPointControlEnabled()
+        def __set__(self, enable):
+            self.flow.enableTwoPointControl(<cbool>enable)
+
     def set_fixed_temp_profile(self, pos, T):
         """Set the fixed temperature profile. This profile is used
         whenever the energy equation is disabled.
@@ -1581,6 +1609,18 @@ cdef class Sim1D:
         def __set__(self, T):
             self.sim.setFixedTemperature(T)
 
+    def set_fuel_side_boundary(self, T):
+        """
+        Set the fuel side internal boundary with approximate temperature.
+        """
+        self.sim.setFuelInternalBoundary(T)
+
+    def set_oxid_side_boundary(self, T):
+        """
+        Set the xoidizer side internal boundary with approximate temperature.
+        """
+        self.sim.setOxidizerInternalBoundary(T)
+
     property fixed_temperature_location:
         """
         Return the location of the point where temperature is fixed for a freely

interfaces/cython/cantera/onedim.py

@@ -274,6 +274,28 @@ def boundary_emissivities(self, epsilon):
             raise ValueError("Boundary emissivities must both be set at the same time.")
         self.flame.boundary_emissivities = epsilon[0], epsilon[1]
 
+    @property
+    def onePointControl_enabled(self):
+        """
+        Get/Set whether or not to active one point flame control.
+        """
+        return self.flame.onePointControl_enabled
+
+    @onePointControl_enabled.setter
+    def onePointControl_enabled(self, enable):
+        self.flame.onePointControl_enabled = enable
+
+    @property
+    def twoPointControl_enabled(self):
+        """
+        Get/Set whether or not to active two point flame control.
+        """
+        return self.flame.twoPointControl_enabled
+
+    @twoPointControl_enabled.setter
+    def twoPointControl_enabled(self, enable):
+        self.flame.twoPointControl_enabled = enable
+
     @property
     def grid(self):
         """ Array of grid point positions along the flame. """
@@ -288,6 +310,24 @@ def P(self):
     def P(self, P):
         self.flame.P = P
 
+    @property
+    def tFuel(self):
+        """ Get/Set the fuel side internal boundary temperature [K] """
+        return self.flame.tFuel
+
+    @tFuel.setter
+    def tFuel(self, T):
+        self.flame.tFuel = T
+
+    @property
+    def tOxid(self):
+        """ Get/Set the oxidizer side internal boundary temperature [K] """
+        return self.flame.tOxid
+
+    @tOxid.setter
+    def tOxid(self, T):
+        self.flame.tOxid = T
+
     @property
     def T(self):
         """ Array containing the temperature [K] at each grid point. """
@@ -325,6 +365,12 @@ def E(self):
             raise AttributeError(
                 "Electric field is only defined for transport model 'ionized_gas'.")
         return self.profile(self.flame, 'eField')
+    @property
+    def Uo(self):
+        """
+        Array containing the oxidizer velocity [m/s] at each point.
+        """
+        return self.profile(self.flame, 'Uo')
 
     def elemental_mass_fraction(self, m):
         r"""

src/oneD/Boundary1D.cpp

@@ -192,11 +192,17 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
             // The flow domain sets this to -rho*u. Add mdot to specify the mass
             // flow rate.
             rb[c_offset_L] += m_mdot;
+        } else if (m_flow->onePointControlEnabled()) {
+            m_mdot = m_flow->density(0)*xb[c_offset_U];
+            rb[c_offset_Uo] += xb[c_offset_U];
+        } else if (m_flow->twoPointControlEnabled()) {
+            m_mdot = m_flow->density(0)*xb[c_offset_U];
         } else {
             // if the flow is a freely-propagating flame, mdot is not specified.
             // Set mdot equal to rho*u, and also set lambda to zero.
             m_mdot = m_flow->density(0) * xb[c_offset_U];
             rb[c_offset_L] = xb[c_offset_L];
+            rb[c_offset_Uo] = xb[c_offset_Uo];
         }
 
         // add the convective term to the species residual equations
@@ -210,11 +216,22 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
         // right inlet
         // Array elements corresponding to the last point in the flow domain
         double* rb = rg + loc() - m_flow->nComponents();
+        double* xb = xg + loc() - m_flow->nComponents();
+        size_t last_index = m_flow->nPoints() - 1;
         rb[c_offset_V] -= m_V0;
-        if (m_flow->doEnergy(m_flow->nPoints() - 1)) {
+        if (m_flow->doEnergy(last_index)) {
             rb[c_offset_T] -= m_temp; // T
         }
-        rb[c_offset_U] += m_mdot; // u
+        // Point-control
+        if (m_flow->onePointControlEnabled() || m_flow->twoPointControlEnabled()) {
+            m_mdot = -(m_flow->density(last_index) * xb[c_offset_Uo]);
+            rb[c_offset_U] += m_mdot; // u
+            rb[c_offset_Uo] += 0;
+        } else {
+            rb[c_offset_U] += m_mdot;
+            rb[c_offset_Uo] += m_mdot/m_flow->density(last_index);
+        }
+
         for (size_t k = 0; k < m_nsp; k++) {
             if (k != m_flow_left->rightExcessSpecies()) {
                 rb[c_offset_Y+k] += m_mdot * m_yin[k];