research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[1D] Move functions from FreeFlame and AxiStagnFlow into StFlow

Browse source 
This makes it possible to implement alternative constitutive relations
(e.g. ionized or non-ideal gases) as a derived class from StFlow and have them
support all of the standard flame configurations (freely propagating, burner
stabilized, counterflow).
This commit is contained in:
Ray Speth 2018-08-01 15:25:02 -04:00
parent 045f3d37bf
commit c1067aa6e9
4 changed files with 102 additions and 161 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
include/cantera/oneD/Domain1D.h
View file

@ -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;

61
include/cantera/oneD/StFlow.h
View file

@ -147,9 +147,14 @@ public:
virtual void restore(const XML_Node& dom, doublereal* soln,
int loglevel);
// overloaded in subclasses
virtual std::string flowType() {
return "<none>";
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;
};
}

196
src/oneD/StFlow.cpp
View file

@ -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

4
src/oneD/boundaries1D.cpp
View file

@ -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);