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

[1D] Use most prevalent species to impose sum-of-mass-fractions constraint

Browse source 
Previously, the first species was always used, which could cause problems with
negative mass fractions in cases where the first species has a mass fraction
near zero.
This commit is contained in:
Ray Speth 2016-03-15 20:13:18 -04:00
parent 011a64aaa2
commit 9f3821db82
3 changed files with 55 additions and 19 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -288,6 +288,16 @@ public:
virtual void evalContinuity(size_t j, doublereal* x, doublereal* r,
integer* diag, doublereal rdt) = 0;
//! Index of the species on the left boundary with the largest mass fraction
size_t leftExcessSpecies() const {
return m_kExcessLeft;
}
//! Index of the species on the right boundary with the largest mass fraction
size_t rightExcessSpecies() const {
return m_kExcessRight;
}
protected:
doublereal component(const doublereal* x, size_t i, size_t j) const {
return x[index(i,j)];
@ -487,6 +497,12 @@ protected:
vector_fp m_zfix;
vector_fp m_tfix;
//! Index of species with a large mass fraction at each boundary, for which
//! the mass fraction may be calculated as 1 minus the sum of the other mass
//! fractions
size_t m_kExcessLeft;
size_t m_kExcessRight;
bool m_dovisc;
//! Update the transport properties at grid points in the range from `j0`

21
src/oneD/StFlow.cpp
View file

@ -24,7 +24,9 @@ StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
m_epsilon_right(0.0),
m_do_soret(false),
m_transport_option(-1),
m_do_radiation(false)
m_do_radiation(false),
m_kExcessLeft(0),
m_kExcessRight(0)
{
m_type = cFlowType;
m_points = points;
@ -262,9 +264,14 @@ void StFlow::eval(size_t jg, doublereal* xg,
// ------------ update properties ------------
updateThermo(x, j0, j1);
// update transport properties only if a Jacobian is not being evaluated
if (jg == npos) {
// update transport properties only if a Jacobian is not being evaluated
updateTransport(x, j0, j1);
double* Yleft = x + index(c_offset_Y, jmin);
m_kExcessLeft = distance(Yleft, max_element(Yleft, Yleft + m_nsp));
double* Yright = x + index(c_offset_Y, jmax);
m_kExcessRight = distance(Yright, max_element(Yright, Yright + m_nsp));
}
// update the species diffusive mass fluxes whether or not a
@ -373,7 +380,7 @@ void StFlow::eval(size_t jg, doublereal* xg,
rsd[index(c_offset_Y + k, 0)] =
-(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
}
rsd[index(c_offset_Y, 0)] = 1.0 - sum;
rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum;
} else if (j == m_points - 1) {
evalRightBoundary(x, rsd, diag, rdt);
} else { // interior points
@ -871,8 +878,8 @@ void AxiStagnFlow::evalRightBoundary(doublereal* x, doublereal* rsd,
sum += Y(x,k,j);
rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
}
rsd[index(4,j)] = 1.0 - sum;
diag[index(4,j)] = 0;
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,
@ -924,8 +931,8 @@ void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd,
sum += Y(x,k,j);
rsd[index(k+4,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
}
rsd[index(4,j)] = 1.0 - sum;
diag[index(4,j)] = 0;
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,

37
src/oneD/boundaries1D.cpp
View file

@ -185,8 +185,10 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg,
rb[3] += x[0];
// add the convective term to the species residual equations
for (size_t k = 1; k < m_nsp; k++) {
rb[4+k] += x[0]*m_yin[k];
for (size_t k = 0; k < m_nsp; k++) {
if (k != m_flow_right->leftExcessSpecies()) {
rb[4+k] += x[0]*m_yin[k];
}
}
// if the flow is a freely-propagating flame, mdot is not specified.
@ -204,8 +206,10 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg,
rb[1] -= m_V0;
rb[2] -= x[1]; // T
rb[0] += x[0]; // u
for (size_t k = 1; k < m_nsp; k++) {
rb[4+k] += x[0]*m_yin[k];
for (size_t k = 0; k < m_nsp; k++) {
if (k != m_flow_left->rightExcessSpecies()) {
rb[4+k] += x[0]*m_yin[k];
}
}
}
}
@ -444,9 +448,12 @@ void Outlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg,
}
rb[2] = xb[2] - xb[2 - nc]; // zero T gradient
for (k = 5; k < nc; k++) {
rb[k] = xb[k] - xb[k - nc]; // zero mass fraction gradient
db[k] = 0;
size_t kSkip = 4 + m_flow_left->rightExcessSpecies();
for (k = 4; k < nc; k++) {
if (k != kSkip) {
rb[k] = xb[k] - xb[k - nc]; // zero mass fraction gradient
db[k] = 0;
}
}
}
}
@ -573,9 +580,12 @@ void OutletRes1D::eval(size_t jg, doublereal* xg, doublereal* rg,
rb[0] = xb[3]; // zero Lambda
}
rb[2] = xb[2] - m_temp; // zero dT/dz
for (k = 5; k < nc; k++) {
rb[k] = xb[k] - m_yres[k-4]; // fixed Y
db[k] = 0;
size_t kSkip = m_flow_left->rightExcessSpecies();
for (k = 4; k < nc; k++) {
if (k != kSkip) {
rb[k] = xb[k] - m_yres[k-4]; // fixed Y
db[k] = 0;
}
}
}
}
@ -795,8 +805,11 @@ void ReactingSurf1D::eval(size_t jg, doublereal* xg, doublereal* rg,
rb =r - nc;
xb = x - nc;
rb[2] = xb[2] - x[0]; // specified T
for (size_t nl = 1; nl < m_left_nsp; nl++) {
rb[4+nl] += m_work[nl]*mwleft[nl];
size_t nSkip = m_flow_left->rightExcessSpecies();
for (size_t nl = 0; nl < m_left_nsp; nl++) {
if (nl != nSkip) {
rb[4+nl] += m_work[nl]*mwleft[nl];
}
}
}
}