diff --git a/include/cantera/oneD/Domain1D.h b/include/cantera/oneD/Domain1D.h index 121cf5750..45ad359c1 100644 --- a/include/cantera/oneD/Domain1D.h +++ b/include/cantera/oneD/Domain1D.h @@ -16,6 +16,8 @@ namespace Cantera // domain types const int cFlowType = 50; +const int cFreeFlow = 51; +const int cAxisymmetricStagnationFlow = 52; const int cConnectorType = 100; const int cSurfType = 102; const int cInletType = 104; diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h index deed089bf..b4721ff21 100644 --- a/include/cantera/oneD/StFlow.h +++ b/include/cantera/oneD/StFlow.h @@ -147,9 +147,14 @@ public: virtual void restore(const XML_Node& dom, doublereal* soln, int loglevel); - // overloaded in subclasses virtual std::string flowType() { - return ""; + if (m_type == cFreeFlow) { + return "Free Flame"; + } else if (m_type == cAxisymmetricStagnationFlow) { + return "Axisymmetric Stagnation"; + } else { + throw CanteraError("StFlow::flowType", "Unknown value for 'm_type'"); + } } void solveEnergyEqn(size_t j=npos); @@ -201,7 +206,7 @@ public: } virtual bool fixed_mdot() { - return true; + return (domainType() != cFreeFlow); } void setViscosityFlag(bool dovisc) { m_dovisc = dovisc; @@ -218,13 +223,13 @@ public: integer* mask, doublereal rdt); //! Evaluate all residual components at the right boundary. - virtual void evalRightBoundary(doublereal* x, doublereal* res, - integer* diag, doublereal rdt) = 0; + virtual void evalRightBoundary(double* x, double* res, int* diag, + double rdt); //! Evaluate the residual corresponding to the continuity equation at all //! interior grid points. - virtual void evalContinuity(size_t j, doublereal* x, doublereal* r, - integer* diag, doublereal rdt) = 0; + 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 { @@ -432,6 +437,13 @@ protected: //! to `j1`, based on solution `x`. virtual void updateTransport(doublereal* x, size_t j0, size_t j1); +public: + //! Location of the point where temperature is fixed + double m_zfixed; + + //! Temperature at the point used to fix the flame location + double m_tfixed; + private: vector_fp m_ybar; }; @@ -446,15 +458,7 @@ public: AxiStagnFlow(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1) : StFlow(ph, nsp, points) { m_dovisc = true; - } - - virtual void evalRightBoundary(doublereal* x, doublereal* res, - integer* diag, doublereal rdt); - virtual void evalContinuity(size_t j, doublereal* x, doublereal* r, - integer* diag, doublereal rdt); - - virtual std::string flowType() { - return "Axisymmetric Stagnation"; + m_type = cAxisymmetricStagnationFlow; } }; @@ -465,28 +469,11 @@ public: class FreeFlame : public StFlow { public: - FreeFlame(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1); - virtual void evalRightBoundary(doublereal* x, doublereal* res, - integer* diag, doublereal rdt); - virtual void evalContinuity(size_t j, doublereal* x, doublereal* r, - integer* diag, doublereal rdt); - - virtual std::string flowType() { - return "Free Flame"; + FreeFlame(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1) : + StFlow(ph, nsp, points) { + m_dovisc = false; + m_type = cFreeFlow; } - virtual bool fixed_mdot() { - return false; - } - virtual void _finalize(const doublereal* x); - virtual void restore(const XML_Node& dom, doublereal* soln, int loglevel); - - virtual XML_Node& save(XML_Node& o, const doublereal* const sol); - - //! Location of the point where temperature is fixed - doublereal m_zfixed; - - //! Temperature at the point used to fix the flame location - doublereal m_tfixed; }; } diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 8ae383295..fa4496146 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -26,7 +26,9 @@ StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) : m_do_multicomponent(false), m_do_radiation(false), m_kExcessLeft(0), - m_kExcessRight(0) + m_kExcessRight(0), + m_zfixed(Undef), + m_tfixed(Undef) { m_type = cFlowType; m_points = points; @@ -85,7 +87,6 @@ StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) : gr.push_back(1.0*ng/m_points); } setupGrid(m_points, gr.data()); - setID("stagnation flow"); // Find indices for radiating species m_kRadiating.resize(2, npos); @@ -204,6 +205,29 @@ void StFlow::_finalize(const doublereal* x) if (e) { solveEnergyEqn(); } + + if (domainType() == cFreeFlow) { + // If the domain contains the temperature fixed point, make sure that it + // is correctly set. This may be necessary when the grid has been modified + // externally. + if (m_tfixed != Undef) { + for (size_t j = 0; j < m_points; j++) { + if (z(j) == m_zfixed) { + return; // fixed point is already set correctly + } + } + + for (size_t j = 0; j < m_points - 1; j++) { + // Find where the temperature profile crosses the current + // fixed temperature. + if ((T(x, j) - m_tfixed) * (T(x, j+1) - m_tfixed) <= 0.0) { + m_tfixed = T(x, j+1); + m_zfixed = z(j+1); + return; + } + } + } + } } void StFlow::eval(size_t jg, doublereal* xg, @@ -750,6 +774,11 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) getFloat(ref, "curve"), getFloat(ref, "prune")); refiner().setGridMin(getFloat(ref, "grid_min")); } + + if (domainType() == cFreeFlow) { + getOptionalFloat(dom, "t_fixed", m_tfixed); + getOptionalFloat(dom, "z_fixed", m_zfixed); + } } XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol) @@ -807,6 +836,10 @@ XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol) addFloat(ref, "curve", refiner().maxSlope()); addFloat(ref, "prune", refiner().prune()); addFloat(ref, "grid_min", refiner().gridMin()); + if (m_zfixed != Undef) { + addFloat(flow, "z_fixed", m_zfixed, "m"); + addFloat(flow, "t_fixed", m_tfixed, "K"); + } return flow; } @@ -872,62 +905,7 @@ void StFlow::fixTemperature(size_t j) } } -void AxiStagnFlow::evalRightBoundary(doublereal* x, doublereal* rsd, - integer* diag, doublereal rdt) -{ - size_t j = m_points - 1; - // the boundary object connected to the right of this one may modify or - // replace these equations. The default boundary conditions are zero u, V, - // and T, and zero diffusive flux for all species. - rsd[index(c_offset_U,j)] = rho_u(x,j); - rsd[index(c_offset_V,j)] = V(x,j); - if (m_do_energy[j]) { - rsd[index(c_offset_T,j)] = T(x,j); - } else { - rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); - } - rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); - diag[index(c_offset_L, j)] = 0; - doublereal sum = 0.0; - for (size_t k = 0; k < m_nsp; k++) { - sum += Y(x,k,j); - rsd[index(k+c_offset_Y,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j); - } - rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum; - diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0; -} - -void AxiStagnFlow::evalContinuity(size_t j, doublereal* x, doublereal* rsd, - integer* diag, doublereal rdt) -{ - //---------------------------------------------- - // Continuity equation - // - // Note that this propagates the mass flow rate information to the left - // (j+1 -> j) from the value specified at the right boundary. The - // lambda information propagates in the opposite direction. - // - // d(\rho u)/dz + 2\rho V = 0 - //------------------------------------------------ - rsd[index(c_offset_U,j)] = - -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j] - -(density(j+1)*V(x,j+1) + density(j)*V(x,j)); - - //algebraic constraint - diag[index(c_offset_U, j)] = 0; -} - -FreeFlame::FreeFlame(IdealGasPhase* ph, size_t nsp, size_t points) : - StFlow(ph, nsp, points), - m_zfixed(Undef), - m_tfixed(Undef) -{ - m_dovisc = false; - setID("flame"); -} - -void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd, - integer* diag, doublereal rdt) +void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt) { size_t j = m_points - 1; @@ -935,10 +913,7 @@ void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd, // replace these equations. The default boundary conditions are zero u, V, // and T, and zero diffusive flux for all species. - // zero gradient - rsd[index(c_offset_U,j)] = rho_u(x,j) - rho_u(x,j-1); rsd[index(c_offset_V,j)] = V(x,j); - rsd[index(c_offset_T,j)] = T(x,j) - T(x,j-1); doublereal sum = 0.0; rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1); diag[index(c_offset_L, j)] = 0; @@ -948,76 +923,53 @@ void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd, } rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum; diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0; -} - -void FreeFlame::evalContinuity(size_t j, doublereal* x, doublereal* rsd, - integer* diag, doublereal rdt) -{ - //---------------------------------------------- - // Continuity equation - // - // d(\rho u)/dz + 2\rho V = 0 - //---------------------------------------------- - if (grid(j) > m_zfixed) { - rsd[index(c_offset_U,j)] = - - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1] - - (density(j-1)*V(x,j-1) + density(j)*V(x,j)); - } else if (grid(j) == m_zfixed) { + if (domainType() == cAxisymmetricStagnationFlow) { + rsd[index(c_offset_U,j)] = rho_u(x,j); if (m_do_energy[j]) { - rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); + rsd[index(c_offset_T,j)] = T(x,j); } else { - rsd[index(c_offset_U,j)] = (rho_u(x,j) - - m_rho[0]*0.3); + rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j); } - } else if (grid(j) < m_zfixed) { - rsd[index(c_offset_U,j)] = - - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j] - - (density(j+1)*V(x,j+1) + density(j)*V(x,j)); + } else if (domainType() == cFreeFlow) { + rsd[index(c_offset_U,j)] = rho_u(x,j) - rho_u(x,j-1); + rsd[index(c_offset_T,j)] = T(x,j) - T(x,j-1); } +} + +void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt) +{ //algebraic constraint diag[index(c_offset_U, j)] = 0; -} - -void FreeFlame::_finalize(const doublereal* x) -{ - StFlow::_finalize(x); - // If the domain contains the temperature fixed point, make sure that it - // is correctly set. This may be necessary when the grid has been modified - // externally. - if (m_tfixed != Undef) { - for (size_t j = 0; j < m_points; j++) { - if (z(j) == m_zfixed) { - return; // fixed point is already set correctly - } - } - - for (size_t j = 0; j < m_points - 1; j++) { - // Find where the temperature profile crosses the current - // fixed temperature. - if ((T(x, j) - m_tfixed) * (T(x, j+1) - m_tfixed) <= 0.0) { - m_tfixed = T(x, j+1); - m_zfixed = z(j+1); - return; + //---------------------------------------------- + // Continuity equation + // + // d(\rho u)/dz + 2\rho V = 0 + //---------------------------------------------- + if (domainType() == cAxisymmetricStagnationFlow) { + // Note that this propagates the mass flow rate information to the left + // (j+1 -> j) from the value specified at the right boundary. The + // lambda information propagates in the opposite direction. + rsd[index(c_offset_U,j)] = + -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j] + -(density(j+1)*V(x,j+1) + density(j)*V(x,j)); + } else if (domainType() == cFreeFlow) { + if (grid(j) > m_zfixed) { + rsd[index(c_offset_U,j)] = + - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1] + - (density(j-1)*V(x,j-1) + density(j)*V(x,j)); + } else if (grid(j) == m_zfixed) { + if (m_do_energy[j]) { + rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed); + } else { + rsd[index(c_offset_U,j)] = (rho_u(x,j) + - m_rho[0]*0.3); } + } else if (grid(j) < m_zfixed) { + rsd[index(c_offset_U,j)] = + - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j] + - (density(j+1)*V(x,j+1) + density(j)*V(x,j)); } } } -void FreeFlame::restore(const XML_Node& dom, doublereal* soln, int loglevel) -{ - StFlow::restore(dom, soln, loglevel); - getOptionalFloat(dom, "t_fixed", m_tfixed); - getOptionalFloat(dom, "z_fixed", m_zfixed); -} - -XML_Node& FreeFlame::save(XML_Node& o, const doublereal* const sol) -{ - XML_Node& flow = StFlow::save(o, sol); - if (m_zfixed != Undef) { - addFloat(flow, "z_fixed", m_zfixed, "m"); - addFloat(flow, "t_fixed", m_tfixed, "K"); - } - return flow; -} - } // namespace diff --git a/src/oneD/boundaries1D.cpp b/src/oneD/boundaries1D.cpp index e001d14cb..214c480df 100644 --- a/src/oneD/boundaries1D.cpp +++ b/src/oneD/boundaries1D.cpp @@ -42,7 +42,7 @@ void Bdry1D::_init(size_t n) // check for left and right flow objects if (m_index > 0) { Domain1D& r = container().domain(m_index-1); - if (r.domainType() == cFlowType) { + if (!r.isConnector()) { // flow domain m_flow_left = (StFlow*)&r; m_left_nv = m_flow_left->nComponents(); m_left_points = m_flow_left->nPoints(); @@ -59,7 +59,7 @@ void Bdry1D::_init(size_t n) // if this is not the last domain, see what is connected on the right if (m_index + 1 < container().nDomains()) { Domain1D& r = container().domain(m_index+1); - if (r.domainType() == cFlowType) { + if (!r.isConnector()) { // flow domain m_flow_right = (StFlow*)&r; m_right_nv = m_flow_right->nComponents(); m_right_loc = container().start(m_index+1);