diff --git a/data/inputs/gri30_ion.cti b/data/inputs/gri30_ion.cti new file mode 100644 index 000000000..e07cd39ca --- /dev/null +++ b/data/inputs/gri30_ion.cti @@ -0,0 +1,88 @@ +units(length='cm', time='s', quantity='mol', act_energy='cal/mol') + +ideal_gas(name='gas', + elements=' O H C N Ar E', + species=['''gri30: H2 H O O2 OH H2O HO2 H2O2 C CH + CH2 CH2(S) CH3 CH4 CO CO2 HCO CH2O CH2OH CH3O + CH3OH C2H C2H2 C2H3 C2H4 C2H5 C2H6 HCCO CH2CO HCCOH + N NH NH2 NH3 NNH NO NO2 N2O HNO CN + HCN H2CN HCNN HCNO HOCN HNCO NCO N2 AR C3H7 + C3H8 CH2CHO CH3CHO''', + 'HCO+ H3O+ E'], + reactions=['gri30: all', 'all'], + transport='Mix', + options=['skip_undeclared_species', 'skip_undeclared_third_bodies'], + initial_state=state(temperature=300.0, pressure=OneAtm)) + +#------------------------------------------------------------------------------- +# Species data +#------------------------------------------------------------------------------- + +species(name = 'HCO+', + atoms = ' H:1 C:1 O:1 E:-1 ', + thermo = ( + NASA( [ 300.00, 1000.00], [ 2.473973600E+00, 8.671559000E-03, + -1.003150000E-05, 6.717052700E-09, -1.787267400E-12, + 9.914660800E+04, 8.175711870E+00] ), + NASA( [ 1000.00, 5000.00], [ 3.741188000E+00, 3.344151700E-03, + -1.239712100E-06, 2.118938800E-10, -1.370415000E-14, + 9.888407800E+04, 2.078613570E+00] ) + ), + transport=gas_transport(geom='linear', + diam=3.59, + well_depth=498.0, + polar=2.5, + rot_relax=0.0), + note = 'J12/70') + +species(name = 'H3O+', + atoms = ' H:3 O:1 E:-1 ', + thermo = ( + NASA( [ 298.15, 1000.00], [ 3.792952700E+00, -9.108540000E-04, + 1.163635490E-05, -1.213648870E-08, 4.261596630E-12, + 7.075124010E+04, 1.471568560E+00] ), + NASA( [ 1000.00, 6000.00], [ 2.496477160E+00, 5.728449200E-03, + -1.839532810E-06, 2.735774390E-10, -1.540939850E-14, + 7.097291130E+04, 7.458507790E+00] ) + ), + transport=gas_transport(geom='nonlinear', + diam=2.605, + well_depth=572.4, + dipole=1.844, + polar=1.5, + rot_relax=2.1), + note = 'TPIS89') + +species(name = 'E', + atoms = ' E:1 ', + thermo = ( + NASA( [ 200.00, 1000.00], [ 2.500000000E+00, 0.000000000E+00, + 0.000000000E+00, 0.000000000E+00, 0.000000000E+00, + -7.453750000E+02, -1.172469020E+01] ), + NASA( [ 1000.00, 6000.00], [ 2.500000000E+00, 0.000000000E+00, + 0.000000000E+00, 0.000000000E+00, 0.000000000E+00, + -7.453750000E+02, -1.172469020E+01] ) + ), + transport=gas_transport(geom='atom', + diam=2.05, + well_depth=145.0, + polar=0.667, + rot_relax=0.0), + note = 'gas L10/92') + +#------------------------------------------------------------------------------- +# Reaction data +#------------------------------------------------------------------------------- + +reaction('CH + O => HCO+ + E', [2.51E+11, 0.0, 1700]) + +reaction('HCO+ + H2O => H3O+ + CO', [1.51E+15, 0.0, 0.0]) + +reaction('H3O+ + E => H2O + H', [2.29E+18, -0.5, 0.0]) + +reaction('H3O+ + E => OH + H + H', [7.95E+21, -1.4, 0.0]) + +reaction('H3O+ + E => H2 + OH', [1.25E+19, -0.5, 0.0]) + +reaction('H3O+ + E => O + H2 + H', [6.0E+17, -0.3, 0.0]) + diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h new file mode 100644 index 000000000..023fe38cb --- /dev/null +++ b/include/cantera/oneD/IonFlow.h @@ -0,0 +1,130 @@ +//! @file IonFlow.h + +// This file is part of Cantera. See License.txt in the top-level directory or +// at http://www.cantera.org/license.txt for license and copyright information. + +#include "Domain1D.h" +#include "cantera/base/Array.h" +#include "cantera/thermo/IdealGasPhase.h" +#include "cantera/kinetics/Kinetics.h" +#include "cantera/oneD/StFlow.h" +#include "cantera/oneD/Sim1D.h" +#include "cantera/IdealGasMix.h" + +namespace Cantera +{ +/** + * A class for ion flow. + * @ingroup onedim + */ +class IonFlow : public FreeFlame +{ +public: + IonFlow(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1); + + //! Turn electric field effect on/off + virtual void enableElectric(bool withElectric); + bool withElectric() const { + return m_do_electric; + } + + virtual void setSolvingPhase(const size_t phase); + + std::vector chargeList() const { + return m_kCharge; + } + + virtual void eval(size_t jg, doublereal* xg, + doublereal* rg, integer* diagg, doublereal rdt); + + virtual void resize(size_t components, size_t points); + + virtual void _finalize(const doublereal* x); + + void solveSpeciesEqn(size_t k=npos); + void fixSpeciesMassFrac(size_t k=npos); + void solvePoissonEqn(size_t j=npos); + void fixElectricPotential(size_t j=npos); + void solveVelocity(size_t j=npos); + void fixVelocity(size_t j=npos); + +protected: + virtual void updateTransport(doublereal* x, size_t j0, size_t j1); + virtual void updateDiffFluxes(const doublereal* x, size_t j0, size_t j1); + virtual void evalPoisson(size_t j, doublereal* x, doublereal* r, integer* diag, doublereal rdt); + virtual void phaseOneDiffFluxes(const doublereal* x, size_t j0, size_t j1); + virtual void phaseTwoDiffFluxes(const doublereal* x, size_t j0, size_t j1); + virtual void phaseThreeDiffFluxes(const doublereal* x, size_t j0, size_t j1); + + bool m_do_electric; + std::vector m_do_velocity; + std::vector m_do_poisson; + + // !electrical properties + vector_int m_speciesCharge; + + // !index of species with charges + std::vector m_kCharge; + + // !index of neutral species + std::vector m_kNeutral; + + // mobility + vector_fp m_mobi; + + // mass fraction of ion by equlibrium + Array2D m_yCharge; + + // IonFlow solving phase + int m_solnPhase; + + // !index of electron + size_t m_kElectron; + + // fixed mass fraction value + vector_fp m_fixedMassFrac; + + // fixed electric potential value + vector_fp m_fixedElecPoten; + + // fixed velocity value + vector_fp m_fixedVelocity; + + //! The fixed electric potential value at point j + doublereal phi_fixed(size_t j) const { + return m_fixedElecPoten[j]; + } + + //! The fixed mass fraction value at point j. + doublereal Y_fixed(size_t k, size_t j) const { + return m_fixedMassFrac[m_points*k+j]; + } + + //! The fixed velocity value at point j + doublereal u_fixed(size_t j) const { + return m_fixedVelocity[j]; + } + + // electric potential + doublereal phi(const doublereal* x, size_t j) const { + return x[index(c_offset_P, j)]; + } + + //electric field + doublereal E(const doublereal* x, size_t j) const { + return -(phi(x,j+1)-phi(x,j))/(z(j+1)-z(j)); + } + + doublereal dEdz(const doublereal* x, size_t j) const { + return 2*(E(x,j)-E(x,j-1))/(z(j+1)-z(j-1)); + } + + // number density + doublereal ND(const doublereal* x, size_t k, size_t j) const { + return Avogadro * m_rho[j] * Y(x,k,j) / m_wt[k]; + } +}; + +} + + diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h index 84cb56e91..0b0d76b7f 100644 --- a/include/cantera/oneD/StFlow.h +++ b/include/cantera/oneD/StFlow.h @@ -23,7 +23,8 @@ const size_t c_offset_U = 0; // axial velocity const size_t c_offset_V = 1; // strain rate const size_t c_offset_T = 2; // temperature const size_t c_offset_L = 3; // (1/r)dP/dr -const size_t c_offset_Y = 4; // mass fractions +const size_t c_offset_P = 4; // electric poisson's equation +const size_t c_offset_Y = 5; // mass fractions class Transport; @@ -179,7 +180,7 @@ public: } //! Change the grid size. Called after grid refinement. - void resize(size_t components, size_t points); + virtual void resize(size_t components, size_t points); virtual void setFixedPoint(int j0, doublereal t0) {} @@ -343,7 +344,7 @@ protected: } //! Update the diffusive mass fluxes. - void updateDiffFluxes(const doublereal* x, size_t j0, size_t j1); + virtual void updateDiffFluxes(const doublereal* x, size_t j0, size_t j1); //--------------------------------------------------------- // member data @@ -414,7 +415,7 @@ protected: //! Update the transport properties at grid points in the range from `j0` //! to `j1`, based on solution `x`. - void updateTransport(doublereal* x, size_t j0, size_t j1); + virtual void updateTransport(doublereal* x, size_t j0, size_t j1); private: vector_fp m_ybar; diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp new file mode 100644 index 000000000..226d71c97 --- /dev/null +++ b/src/oneD/IonFlow.cpp @@ -0,0 +1,461 @@ +//! @file IonFlow.cpp + +// This file is part of Cantera. See License.txt in the top-level directory or +// at http://www.cantera.org/license.txt for license and copyright information. + +#include "cantera/oneD/IonFlow.h" +#include "cantera/oneD/StFlow.h" +#include "cantera/base/ctml.h" +#include "cantera/transport/TransportBase.h" +#include "cantera/numerics/funcs.h" +#include "cantera/oneD/Domain1D.h" + + +using namespace std; + +namespace Cantera +{ + +IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) : + FreeFlame(ph, nsp, points), + m_do_electric(false), + m_solnPhase(1) +{ + // make a local copy of species charge + for (size_t k = 0; k < m_nsp; k++) { + m_speciesCharge.push_back(m_thermo->charge(k)); + } + + // Find indices for charge of species + for (size_t k = 0; k < m_nsp; k++){ + if (m_speciesCharge[k] != 0){ + m_kCharge.push_back(k); + } else { + m_kNeutral.push_back(k); + } + } + + // Find the index of electron + if (m_thermo->speciesIndex("E") < m_nsp ) { + m_kElectron = m_thermo->speciesIndex("E"); + } + + // mass fraction bounds (strict bound for ions) + for (size_t k : m_kCharge) { + setBounds(c_offset_Y+k, -1.0e-20, 1e-5); + } + setBounds(c_offset_P, -1.0e20, 1.0e20); + m_refiner->setActive(c_offset_P, false); + + m_mobi.resize(m_nsp*m_points); + m_do_poisson.resize(m_points,false); + m_do_velocity.resize(m_points,true); +} + +void IonFlow::resize(size_t components, size_t points){ + StFlow::resize(components, points); + m_mobi.resize(m_nsp*m_points); + m_do_species.resize(m_nsp,true); + m_do_poisson.resize(m_points,false); + m_do_velocity.resize(m_points,true); + m_fixedMassFrac.resize(m_points*m_nsp); + m_fixedElecPoten.resize(m_points,0.0); + m_fixedVelocity.resize(m_points); +} + +void IonFlow::updateTransport(doublereal* x, size_t j0, size_t j1) +{ + StFlow::updateTransport(x,j0,j1); + for (size_t j = j0; j < j1; j++) { + setGasAtMidpoint(x,j); + m_trans->getMobilities(&m_mobi[j*m_nsp]); + m_mobi[m_kElectron+m_nsp*j] = 0.4; + m_diff[m_kElectron+m_nsp*j] = 0.4*(Boltzmann * T(x,j)) / ElectronCharge; + } +} +void IonFlow::updateDiffFluxes(const doublereal* x, size_t j0, size_t j1) +{ + if (m_solnPhase == 1) { + phaseOneDiffFluxes(x,j0,j1); + } else if (m_solnPhase == 2) { + phaseTwoDiffFluxes(x,j0,j1); + } else { + phaseThreeDiffFluxes(x,j0,j1); + } +} + +void IonFlow::phaseOneDiffFluxes(const doublereal* x, size_t j0, size_t j1) +{ + for (size_t j = j0; j < j1; j++) { + double wtm = m_wtm[j]; + double rho = density(j); + double dz = z(j+1) - z(j); + double sum = 0.0; + for (size_t k : m_kNeutral) { + m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm); + m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz; + sum -= m_flux(k,j); + } + + // correction flux to insure that \sum_k Y_k V_k = 0. + for (size_t k : m_kNeutral) { + m_flux(k,j) += sum*Y(x,k,j); + } + + // flux for ions + for (size_t k : m_kCharge) { + m_flux(k,j) = 0; + } + } +} + +void IonFlow::phaseTwoDiffFluxes(const doublereal* x, size_t j0, size_t j1) +{ + for (size_t j = j0; j < j1; j++) { + double wtm = m_wtm[j]; + double rho = density(j); + double dz = z(j+1) - z(j); + // mixture-average diffusion + double sum_flux = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm); + m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz; + sum_flux -= m_flux(k,j); + } + + // ambipolar diffusion + double sum_chargeFlux = 0.0; + double sum = 0.0; + for (size_t k : m_kCharge) { + double Xav = 0.5 * (X(x,k,j+1) + X(x,k,j)); + int q_k = m_speciesCharge[k]; + sum_chargeFlux += m_speciesCharge[k] / m_wt[k] * m_flux(k,j); + sum += m_mobi[k+m_nsp*j] * Xav * q_k * q_k; + } + double drift; + double sum_drift = 0.0; + for (size_t k : m_kCharge) { + double Xav = 0.5 * (X(x,k,j+1) + X(x,k,j)); + int q_k = m_speciesCharge[k]; + drift = q_k * q_k * m_mobi[k+m_nsp*j] * Xav / sum; + drift *= -sum_chargeFlux * m_wt[k] / q_k; + m_flux(k,j) += drift; + sum_drift -= drift; + } + + // correction flux + for (size_t k = 0; k < m_nsp; k++) { + m_flux(k,j) += Y(x,k,j) * sum_flux; + } + } +} + +void IonFlow::phaseThreeDiffFluxes(const doublereal* x, size_t j0, size_t j1) +{ + for (size_t j = j0; j < j1; j++) { + double wtm = m_wtm[j]; + double rho = density(j); + double dz = z(j+1) - z(j); + // mixture-average diffusion + double sum = 0.0; + for (size_t k = 0; k < m_nsp; k++) { + m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm); + m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz; + sum -= m_flux(k,j); + } + + // correction flux + for (size_t k = 0; k < m_nsp; k++) { + m_flux(k,j) += Y(x,k,j) * sum; + } + + // ambipolar diffusion + double drift; + double E_ambi = E(x,j); + sum = 0.0; + for (size_t k : m_kCharge) { + double Yav = 0.5 * (Y(x,k,j) + Y(x,k,j+1)); + drift = rho * Yav * E_ambi; + drift *= m_speciesCharge[k] * m_mobi[k+m_nsp*j]; + m_flux(k,j) += drift; + sum -= drift; + } + + // correction drift + for (size_t k : m_kCharge) { + m_flux(k,j) += Y(x,k,j) * sum; + } + } +} + +void IonFlow::enableElectric(bool withElectric) +{ + m_do_electric = withElectric; +} + +void IonFlow::setSolvingPhase(const size_t phase) +{ + m_solnPhase = phase; +} + +void IonFlow::eval(size_t jg, doublereal* xg, + doublereal* rg, integer* diagg, doublereal rdt) +{ + StFlow::eval(jg, xg, rg, diagg, rdt); + // start of local part of global arrays + doublereal* x = xg + loc(); + doublereal* rsd = rg + loc(); + integer* diag = diagg + loc(); + + size_t jmin, jmax; + if (jg == npos) { // evaluate all points + jmin = 0; + jmax = m_points - 1; + } else { // evaluate points for Jacobian + size_t jpt = (jg == 0) ? 0 : jg - firstPoint(); + jmin = std::max(jpt, 1) - 1; + jmax = std::min(jpt+1,m_points-1); + } + // the boundary points are not applied + for (size_t j = jmin; j <= jmax; j++) { + if (j == 0) { + rsd[index(c_offset_P, j)] = -phi(x,j); + diag[index(c_offset_P, j)] = 0; + for ( size_t k : m_kCharge) { + rsd[index(c_offset_Y + k, j)] = Y(x,k,j); + diag[index(c_offset_Y + k, j)] = 0; + } + } else if (j == m_points - 1) { + rsd[index(c_offset_P, j)] = -phi(x,j); + diag[index(c_offset_P, j)] = 0; + for ( size_t k : m_kCharge) { + rsd[index(c_offset_Y + k, j)] = Y(x,k,j); + diag[index(c_offset_Y + k, j)] = 0; + } + } else { + if (!m_do_velocity[j]) { + rsd[index(c_offset_U, j)] = u(x,j) - u_fixed(j); + diag[index(c_offset_U, j)] = 0; + } + for (size_t k = 0; k < m_nsp; k++) { + if (!m_do_species[k]) { + rsd[index(c_offset_Y + k, j)] = Y(x,k,j) - Y_fixed(k,j); + rsd[index(c_offset_Y + k, j)] -= rdt*(Y(x,k,j) - Y_prev(k,j)); + diag[index(c_offset_Y + k, j)] = 1; + } + } + } + } + + // convinent method due to interference + for (size_t j = jmin; j <= jmax; j++) { + if (j == 0) { + rsd[index(c_offset_P, j)] = -phi(x,j); + diag[index(c_offset_P, j)] = 0; + } else if (j == m_points - 1) { + rsd[index(c_offset_P, j)] = -phi(x,j); + diag[index(c_offset_P, j)] = 0; + } else { + if (m_do_poisson[j]) { + evalPoisson(j,x,rsd,diag,rdt); + } else { + rsd[index(c_offset_P, j)] = phi(x,j) - phi_fixed(j); + diag[index(c_offset_P, j)] = 0; + } + } + } +} + +void IonFlow::evalPoisson(size_t j, doublereal* x, doublereal* rsd, integer* diag, doublereal rdt) +{ + //----------------------------------------------- + // Poisson's equation + // + // dE/dz = e/eps_0 * sum(q_k*n_k) + // + // E = -dV/dz + //----------------------------------------------- + doublereal chargeDensity = 0.0; + for (size_t k : m_kCharge) { + chargeDensity += m_speciesCharge[k] * ElectronCharge * ND(x,k,j); + } + rsd[index(c_offset_P, j)] = dEdz(x,j) - chargeDensity / epsilon_0; + diag[index(c_offset_P, j)] = 0; +} + +void IonFlow::solveSpeciesEqn(size_t k) +{ + bool changed = false; + if (k == npos) { + for (size_t i = 0; i < m_nsp; i++) { + if (!m_do_energy[i]) { + changed = true; + } + m_do_species[i] = true; + } + } else { + if (!m_do_species[k]) { + changed = true; + } + m_do_species[k] = true; + } + if (changed) { + needJacUpdate(); + } +} + +void IonFlow::fixSpeciesMassFrac(size_t k) +{ + bool changed = false; + if (k == npos) { + for (size_t i = 0; i < m_nsp; i++) { + if (m_do_species[i]) { + changed = true; + } + m_do_species[i] = false; + } + } else { + if (m_do_species[k]) { + changed = true; + } + m_do_species[k] = false; + } + if (changed) { + needJacUpdate(); + } +} + +void IonFlow::solvePoissonEqn(size_t j) +{ + bool changed = false; + if (j == npos) { + for (size_t i = 0; i < m_points; i++) { + if (!m_do_poisson[i]) { + changed = true; + } + m_do_poisson[i] = true; + } + } else { + if (!m_do_poisson[j]) { + changed = true; + } + m_do_poisson[j] = true; + } + m_refiner->setActive(0, true); + m_refiner->setActive(1, true); + m_refiner->setActive(2, true); + if (changed) { + needJacUpdate(); + } +} + +void IonFlow::fixElectricPotential(size_t j) +{ + bool changed = false; + if (j == npos) { + for (size_t i = 0; i < m_points; i++) { + if (m_do_poisson[i]) { + changed = true; + } + m_do_poisson[i] = false; + } + } else { + if (m_do_poisson[j]) { + changed = true; + } + m_do_poisson[j] = false; + } + m_refiner->setActive(0, false); + m_refiner->setActive(1, false); + m_refiner->setActive(2, false); + if (changed) { + needJacUpdate(); + } +} + +void IonFlow::solveVelocity(size_t j) +{ + bool changed = false; + if (j == npos) { + for (size_t i = 0; i < m_points; i++) { + if (!m_do_velocity[i]) { + changed = true; + } + m_do_velocity[i] = true; + } + } else { + if (!m_do_velocity[j]) { + changed = true; + } + m_do_velocity[j] = true; + } + m_refiner->setActive(0, true); + m_refiner->setActive(1, true); + m_refiner->setActive(2, true); + if (changed) { + needJacUpdate(); + } +} + +void IonFlow::fixVelocity(size_t j) +{ + bool changed = false; + if (j == npos) { + for (size_t i = 0; i < m_points; i++) { + if (m_do_velocity[i]) { + changed = true; + } + m_do_velocity[i] = false; + } + } else { + if (m_do_velocity[j]) { + changed = true; + } + m_do_velocity[j] = false; + } + m_refiner->setActive(0, false); + m_refiner->setActive(1, false); + m_refiner->setActive(2, false); + if (changed) { + needJacUpdate(); + } +} + +void IonFlow::_finalize(const doublereal* x) +{ + FreeFlame::_finalize(x); + + for (size_t k = 0; k < m_nsp; k++) { + bool y = m_do_species[k]; + if (!y) { + for (size_t j = 0; j < m_points; j++) { + m_fixedMassFrac[m_points*k+j] = Y(x,k,j); + } + } + } + + // This method is still not tested + // not sure why you want to return to original state + // if not doing on point zero + bool p = m_do_poisson[0]; + for (size_t j = 0; j < m_points; j++) { + if (!p) { + m_fixedElecPoten[j] = phi(x, j); + } + } + if (p) { + solvePoissonEqn(); + } + + bool v = m_do_velocity[0]; + for (size_t j = 0; j < m_points; j++) { + if (!v) { + m_fixedVelocity[j] = u(x,j); + } + } + if (v) { + solveVelocity(); + } +} + +} \ No newline at end of file diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 7f5c8429f..6b3792191 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -132,7 +132,8 @@ void StFlow::setupGrid(size_t n, const doublereal* z) } } -void StFlow::resetBadValues(double* xg) { +void StFlow::resetBadValues(double* xg) +{ double* x = xg + loc(); for (size_t j = 0; j < m_points; j++) { double* Y = x + m_nv*j + c_offset_Y; @@ -368,10 +369,17 @@ void StFlow::eval(size_t jg, doublereal* xg, -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0)); } rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum; + + // set residual of poisson's equ to zero + rsd[index(c_offset_P, 0)] = x[index(c_offset_P, j)]; } else if (j == m_points - 1) { evalRightBoundary(x, rsd, diag, rdt); + // set residual of poisson's equ to zero + rsd[index(c_offset_P, j)] = x[index(c_offset_P, j)]; } else { // interior points evalContinuity(j, x, rsd, diag, rdt); + // set residual of poisson's equ to zero + rsd[index(c_offset_P, j)] = x[index(c_offset_P, j)]; //------------------------------------------------ // Radial momentum equation @@ -546,6 +554,8 @@ string StFlow::componentName(size_t n) const return "T"; case 3: return "lambda"; + case 4: + return "ePotential"; default: if (n >= c_offset_Y && n < (c_offset_Y + m_nsp)) { return m_thermo->speciesName(n - c_offset_Y); @@ -565,6 +575,8 @@ size_t StFlow::componentIndex(const std::string& name) const return 2; } else if (name=="lambda") { return 3; + } else if (name == "ePotential") { + return 4; } else { for (size_t n=c_offset_Y; n