From 52dbe8c007921769287127701017599e37460598 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 29 Mar 2017 18:46:30 -0400 Subject: [PATCH] [1D] Correct handling of boundary conditions when energy equation is disabled --- interfaces/cython/cantera/test/test_onedim.py | 14 +++++++ src/oneD/StFlow.cpp | 12 +++++- src/oneD/boundaries1D.cpp | 38 +++++++++++++------ 3 files changed, 51 insertions(+), 13 deletions(-) diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py index 37c7e9dc8..945d8f62f 100644 --- a/interfaces/cython/cantera/test/test_onedim.py +++ b/interfaces/cython/cantera/test/test_onedim.py @@ -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): diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 673fcae14..7f5c8429f 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -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; diff --git a/src/oneD/boundaries1D.cpp b/src/oneD/boundaries1D.cpp index c8ef339ad..677dc1f4f 100644 --- a/src/oneD/boundaries1D.cpp +++ b/src/oneD/boundaries1D.cpp @@ -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) {