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

[1D] Correct handling of boundary conditions when energy equation is disabled

Browse source
This commit is contained in:
Ray Speth 2017-03-29 18:46:30 -04:00
parent bfdc2b9e1d
commit 52dbe8c007
3 changed files with 51 additions and 13 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

14
interfaces/cython/cantera/test/test_onedim.py
View file

@ -784,6 +784,20 @@ class TestBurnerFlame(utilities.CanteraTest):
def test_case5(self):
self.solve(phi=1.0, T=400, width=0.2, P=0.01)
def test_fixed_temp(self):
gas = ct.Solution('h2o2.xml')
gas.TPX = 400, 2*ct.one_atm, {'H2':0.7, 'O2':0.5, 'AR':1.5}
sim = ct.BurnerFlame(gas=gas, width=0.05)
sim.burner.mdot = gas.density * 0.15
sim.flame.set_fixed_temp_profile([0, 0.1, 0.9, 1],
[400, 1100, 1100, 500])
sim.energy_enabled = False
sim.solve(loglevel=0, refine_grid=True)
self.assertNear(sim.T[0], 400)
self.assertNear(sim.T[-1], 500)
self.assertNear(max(sim.T), 1100)
class TestImpingingJet(utilities.CanteraTest):
def run_reacting_surface(self, xch4, tsurf, mdot, width):

12
src/oneD/StFlow.cpp
View file

@ -352,7 +352,11 @@ void StFlow::eval(size_t jg, doublereal* xg,
// a result, these residual equations will force the solution
// variables to the values for the boundary object
rsd[index(c_offset_V,0)] = V(x,0);
rsd[index(c_offset_T,0)] = T(x,0);
if (doEnergy(0)) {
rsd[index(c_offset_T,0)] = T(x,0);
} else {
rsd[index(c_offset_T,0)] = T(x,0) - T_fixed(0);
}
rsd[index(c_offset_L,0)] = -rho_u(x,0);
// The default boundary condition for species is zero flux. However,
@ -862,7 +866,11 @@ void AxiStagnFlow::evalRightBoundary(doublereal* x, doublereal* rsd,
// and T, and zero diffusive flux for all species.
rsd[index(0,j)] = rho_u(x,j);
rsd[index(1,j)] = V(x,j);
rsd[index(2,j)] = T(x,j);
if (m_do_energy[j]) {
rsd[index(2,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;

38
src/oneD/boundaries1D.cpp
View file

@ -166,9 +166,11 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg,
// so for finite spreading rate subtract m_V0.
rb[1] -= m_V0;
// The third flow residual is for T, where it is set to T(0). Subtract
// the local temperature to hold the flow T to the inlet T.
rb[2] -= m_temp;
if (m_flow->doEnergy(0)) {
// The third flow residual is for T, where it is set to T(0). Subtract
// the local temperature to hold the flow T to the inlet T.
rb[2] -= m_temp;
}
if (m_flow->fixed_mdot()) {
// The flow domain sets this to -rho*u. Add mdot to specify the mass
@ -193,7 +195,9 @@ void Inlet1D::eval(size_t jg, doublereal* xg, doublereal* rg,
// Array elements corresponding to the flast point in the flow domain
double* rb = rg + loc() - m_flow->nComponents();
rb[1] -= m_V0;
rb[2] -= m_temp; // T
if (m_flow->doEnergy(m_flow->nPoints() - 1)) {
rb[2] -= m_temp; // T
}
rb[0] += m_mdot; // u
for (size_t k = 0; k < m_nsp; k++) {
if (k != m_flow_left->rightExcessSpecies()) {
@ -288,7 +292,9 @@ void Symm1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg,
db[1] = 0;
db[2] = 0;
rb[1] = xb[1] - xb[1 + nc]; // zero dV/dz
rb[2] = xb[2] - xb[2 + nc]; // zero dT/dz
if (m_flow_right->doEnergy(0)) {
rb[2] = xb[2] - xb[2 + nc]; // zero dT/dz
}
}
if (m_flow_left) {
@ -299,7 +305,9 @@ void Symm1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg,
db[1] = 0;
db[2] = 0;
rb[1] = xb[1] - xb[1 - nc]; // zero dV/dz
rb[2] = xb[2] - xb[2 - nc]; // zero dT/dz
if (m_flow_left->doEnergy(m_flow_left->nPoints() - 1)) {
rb[2] = xb[2] - xb[2 - nc]; // zero dT/dz
}
}
}
@ -355,7 +363,9 @@ void Outlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg,
double* xb = x;
double* rb = r;
rb[0] = xb[3];
rb[2] = xb[2] - xb[2 + nc];
if (m_flow_right->doEnergy(0)) {
rb[2] = xb[2] - xb[2 + nc];
}
for (size_t k = c_offset_Y; k < nc; k++) {
rb[k] = xb[k] - xb[k + nc];
}
@ -372,7 +382,9 @@ void Outlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg,
rb[0] = xb[3];
}
rb[2] = xb[2] - xb[2 - nc]; // zero T gradient
if (m_flow_left->doEnergy(m_flow_left->nPoints()-1)) {
rb[2] = xb[2] - xb[2 - nc]; // zero T gradient
}
size_t kSkip = c_offset_Y + m_flow_left->rightExcessSpecies();
for (size_t k = c_offset_Y; k < nc; k++) {
if (k != kSkip) {
@ -459,8 +471,10 @@ void OutletRes1D::eval(size_t jg, doublereal* xg, doublereal* rg,
// zero Lambda
rb[0] = xb[3];
// zero gradient for T
rb[2] = xb[2] - xb[2 + nc];
if (m_flow_right->doEnergy(0)) {
// zero gradient for T
rb[2] = xb[2] - xb[2 + nc];
}
// specified mass fractions
for (size_t k = c_offset_Y; k < nc; k++) {
@ -479,7 +493,9 @@ void OutletRes1D::eval(size_t jg, doublereal* xg, doublereal* rg,
} else {
rb[0] = xb[3]; // zero Lambda
}
rb[2] = xb[2] - m_temp; // zero dT/dz
if (m_flow_left->doEnergy(m_flow_left->nPoints()-1)) {
rb[2] = xb[2] - m_temp; // zero dT/dz
}
size_t kSkip = m_flow_left->rightExcessSpecies();
for (size_t k = c_offset_Y; k < nc; k++) {
if (k != kSkip) {