diff --git a/include/cantera/oneD/Sim1D.h b/include/cantera/oneD/Sim1D.h index 450a5d00f13..684b596dd6a 100644 --- a/include/cantera/oneD/Sim1D.h +++ b/include/cantera/oneD/Sim1D.h @@ -196,11 +196,11 @@ class Sim1D : public OneDim double fixedTemperatureLocation(); //! ------ One and Two-Point flame control methods - //! Set the left control point location. This is used for one or two + //! Set the left control point location. This is used for two //! point flame control. void setLeftControlPoint(double temperature); - //! Set the right control point location. This is used for one or two + //! Set the right control point location. This is used for two //! point flame control. void setRightControlPoint(double temperature); //! ------------------- diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h index 2c82e7b5768..e39d1bd2380 100644 --- a/include/cantera/oneD/StFlow.h +++ b/include/cantera/oneD/StFlow.h @@ -242,14 +242,21 @@ class StFlow : public Domain1D //! the value of the solution at these points is fixed. The values of the control //! points are dictated and thus serve as a boundary condition that affects the //! solution of the governing equations in the 1D domain. The imposition of fixed - //! points in the domain means that the original set of governing equations' boundary + //! points in the domain means that the original set of governing equations' boundary //! conditions would over-specify the problem. Thus, the boundary conditions are changed //! to reflect the fact that the control points are serving as internal boundary conditions. //! + //! In this method, the imposition of the two internal boundary conditions requires that two other + //! boundary conditions be changed. The first is the boundary condition for the continuity equation + //! at the left boundary, which is changed to be a value that is derived from the solution at the + //! left boundary. The second is the continuity boundary condition at the right boundary, which is + //! also determined from the flow solution by using the oxidizer axial velocity equation variable to + //! compute the mass flux at the right boundary. + //! //! This method is based on the work of M. Nishioka, C.K. Law, and T. Takeno (1996) titled - //! "A Flame-Controlling Continuation Method for Generating S-Curve Responses with + //! "A Flame-Controlling Continuation Method for Generating S-Curve Responses with //! Detailed Chemistry" - + //! The current left control point temperature double leftControlPointTemperature() const { if (m_twoPointControl && (m_zLeft != Undef)) { @@ -782,7 +789,7 @@ class StFlow : public Domain1D //! ------------- - + private: vector m_ybar; }; diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp index fc08c012a81..fd346b488c0 100644 --- a/src/oneD/Boundary1D.cpp +++ b/src/oneD/Boundary1D.cpp @@ -200,6 +200,11 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg, m_mdot = m_flow->density(0) * xb[c_offset_U]; } else if (m_flow->isStrained()) { //axisymmetric flow if (m_flow->twoPointControlEnabled()) { + // When using two-point control, the mass flow rate at the left inlet is not + // specified. Instead, the mass flow rate is dictated by the velocity at the + // left inlet, which comes from the U variable. The default boundary condition specified + // in the StFlow.cpp file already handles this case. We only need to update the stored + // value of m_mdot so that other equations that use the quantity are consistent. m_mdot = m_flow->density(0)*xb[c_offset_U]; } else { // The flow domain sets this to -rho*u. Add mdot to specify the mass @@ -237,12 +242,13 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg, if (m_flow->twoPointControlEnabled()) {// For point control adjustments // At the right boundary, the mdot is dictated by the velocity at the - // right boundary, which comes from the Uo variable. + // right boundary, which comes from the Uo variable. The variable Uo is + // the left-moving velocity and has a negative value, so the mass flow has + // to be negated to give a positive value when using Uo. m_mdot = -(m_flow->density(last_index) * xb[c_offset_Uo]); - rb[c_offset_U] += m_mdot; + rb[c_offset_U] += m_mdot; } 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++) { diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 2920b6b5535..3b31c29d24f 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -412,7 +412,7 @@ void StFlow::evalContinuity(double* x, double* rsd, int* diag, rsd[index(c_offset_U,jmin)] = -(rho_u(x,jmin + 1) - rho_u(x,jmin))/m_dz[jmin] -(density(jmin + 1)*V(x,jmin + 1) + density(jmin)*V(x,jmin)); - + diag[index(c_offset_U,jmin)] = 0; // Algebraic constraint } @@ -436,7 +436,7 @@ void StFlow::evalContinuity(double* x, double* rsd, int* diag, // 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)); - + diag[index(c_offset_U, j)] = 0; // Algebraic constraint } } else if (m_isFree) { // "free-flow" @@ -505,7 +505,7 @@ void StFlow::evalLambda(double* x, double* rsd, int* diag, } return; } - + if (jmin == 0) { // left boundary if (m_twoPointControl) { rsd[index(c_offset_L, jmin)] = lambda(x,jmin+1) - lambda(x,jmin); @@ -522,7 +522,7 @@ void StFlow::evalLambda(double* x, double* rsd, int* diag, // j0 and j1 are constrained to only interior points size_t j0 = std::max(jmin, 1); size_t j1 = std::min(jmax, m_points - 2); - double epsilon = 1e-5; // Precision threshold for being 'equal' to a coordinate + double epsilon = 1e-8; // Precision threshold for being 'equal' to a coordinate for (size_t j = j0; j <= j1; j++) { // interior points if (m_twoPointControl) { if (std::abs(grid(j) - m_zLeft) < epsilon ) { @@ -590,15 +590,13 @@ void StFlow::evalUo(double* x, double* rsd, int* diag, if (jmin == 0) { // left boundary // Because the Uo equation is used for two-point control, the boundary // for Uo is located in the domain interior at the right control point, - // thus at the boundary, the + // thus at the boundary, the rsd[index(c_offset_Uo,jmin)] = Uo(x,jmin+1) - Uo(x,jmin); } if (jmax == m_points - 1) { // right boundary if(m_twoPointControl) { rsd[index(c_offset_Uo, jmax)] = Uo(x,jmax) - Uo(x,jmax-1); - } else { - rsd[index(c_offset_Uo, jmax)] = Uo(x,jmax); } diag[index(c_offset_Uo, jmax)] = 0; } @@ -606,7 +604,7 @@ void StFlow::evalUo(double* x, double* rsd, int* diag, // j0 and j1 are constrained to only interior points size_t j0 = std::max(jmin, 1); size_t j1 = std::min(jmax, m_points - 2); - double epsilon = 1e-5; // Precision threshold for being 'equal' to a coordinate + double epsilon = 1e-8; // Precision threshold for being 'equal' to a coordinate for (size_t j = j0; j <= j1; j++) { // interior points if (m_twoPointControl) { if (std::abs(grid(j) - m_zRight) < epsilon) {