From fa9b9374cfe8c11349866a475c4ebf535f47368a Mon Sep 17 00:00:00 2001 From: BangShiuh Date: Mon, 30 Apr 2018 20:05:15 -0400 Subject: [PATCH] [1D] Add polyfit for electron transport profile --- include/cantera/oneD/IonFlow.h | 29 +++++------- interfaces/cython/cantera/test/test_onedim.py | 2 +- src/oneD/IonFlow.cpp | 47 ++++++++++--------- 3 files changed, 40 insertions(+), 38 deletions(-) diff --git a/include/cantera/oneD/IonFlow.h b/include/cantera/oneD/IonFlow.h index 39c457f02..3fa83404b 100644 --- a/include/cantera/oneD/IonFlow.h +++ b/include/cantera/oneD/IonFlow.h @@ -59,8 +59,9 @@ public: * If in the future the class GasTranport is improved, this method may * be discard. This method specifies this profile. */ - void setElectronTransport(vector_fp& zfixed, vector_fp& diff_e_fixed, - vector_fp& mobi_e_fixed); + void setElectronTransport(vector_fp& tfix, + vector_fp& diff_e, + vector_fp& mobi_e); protected: /*! @@ -81,9 +82,6 @@ protected: //! flag for importing transport of electron bool m_import_electron_transport; - //! flag for overwrite transport of electron or not - bool m_overwrite_eTransport; - //! electrical properties vector_int m_speciesCharge; @@ -93,9 +91,9 @@ protected: //! index of neutral species std::vector m_kNeutral; - //! fixed transport profile of electron - vector_fp m_elecMobility; - vector_fp m_elecDiffCoeff; + //! coefficients of polynomial fitting of fixed electron transport profile + vector_fp m_mobi_e_fix; + vector_fp m_diff_e_fix; //! mobility vector_fp m_mobility; @@ -113,11 +111,6 @@ protected: //! fixed electric potential value vector_fp m_fixedElecPoten; - //! fixed electron transport values - vector_fp m_ztfix; - vector_fp m_diff_e_fix; - vector_fp m_mobi_e_fix; - //! The fixed electric potential value at point j double phi_fixed(size_t j) const { return m_fixedElecPoten[j]; @@ -142,9 +135,13 @@ protected: return Avogadro * m_rho[j] * Y(x,k,j) / m_wt[k]; } - //! total number density - double ND_t(size_t j) const { - return Avogadro * m_rho[j] / m_wtm[j]; + //! total charge density + double rho_e(double* x, size_t j) const { + double chargeDensity = 0.0; + for (size_t k : m_kCharge) { + chargeDensity += m_speciesCharge[k] * ElectronCharge * ND(x,k,j); + } + return chargeDensity; } }; diff --git a/interfaces/cython/cantera/test/test_onedim.py b/interfaces/cython/cantera/test/test_onedim.py index d28625601..12919494c 100644 --- a/interfaces/cython/cantera/test/test_onedim.py +++ b/interfaces/cython/cantera/test/test_onedim.py @@ -950,4 +950,4 @@ class TestIonFlame(utilities.CanteraTest): self.sim.solve(loglevel=0, stage=2, enable_energy=True) # Regression test - self.assertNear(max(self.sim.E), 113.5274, 1e-3) + self.assertNear(max(self.sim.E), 114.4623, 1e-3) diff --git a/src/oneD/IonFlow.cpp b/src/oneD/IonFlow.cpp index f24f69cb7..267a2e9db 100644 --- a/src/oneD/IonFlow.cpp +++ b/src/oneD/IonFlow.cpp @@ -8,6 +8,7 @@ #include "cantera/base/ctml.h" #include "cantera/transport/TransportBase.h" #include "cantera/numerics/funcs.h" +#include "cantera/numerics/polyfit.h" using namespace std; @@ -17,7 +18,6 @@ namespace Cantera IonFlow::IonFlow(IdealGasPhase* ph, size_t nsp, size_t points) : FreeFlame(ph, nsp, points), m_import_electron_transport(false), - m_overwrite_eTransport(true), m_stage(1), m_inletVoltage(0.0), m_outletVoltage(0.0), @@ -56,8 +56,6 @@ void IonFlow::resize(size_t components, size_t points){ m_do_species.resize(m_nsp,true); m_do_poisson.resize(m_points,false); m_fixedElecPoten.resize(m_points,0.0); - m_elecMobility.resize(m_points); - m_elecDiffCoeff.resize(m_points); } void IonFlow::updateTransport(double* x, size_t j0, size_t j1) @@ -66,14 +64,11 @@ void IonFlow::updateTransport(double* x, size_t j0, size_t j1) for (size_t j = j0; j < j1; j++) { setGasAtMidpoint(x,j); m_trans->getMobilities(&m_mobility[j*m_nsp]); - if (m_overwrite_eTransport && (m_kElectron != npos)) { - if (m_import_electron_transport) { - m_mobility[m_kElectron+m_nsp*j] = m_elecMobility[j]; - m_diff[m_kElectron+m_nsp*j] = m_elecDiffCoeff[j]; - } else { - m_mobility[m_kElectron+m_nsp*j] = 0.4; - m_diff[m_kElectron+m_nsp*j] = 0.4*(Boltzmann * T(x,j)) / ElectronCharge; - } + if (m_import_electron_transport) { + size_t k = m_kElectron; + double tlog = log(m_thermo->temperature()); + m_mobility[k+m_nsp*j] = poly5(tlog, m_mobi_e_fix.data()); + m_diff[k+m_nsp*j] = poly5(tlog, m_diff_e_fix.data()); } } } @@ -184,6 +179,12 @@ void IonFlow::evalResidual(double* x, double* rsd, int* diag, for (size_t j = jmin; j <= jmax; j++) { if (j == 0) { + // enforcing the flux for charged species is difficult + // since charged species are also affected by electric + // force, so Neumann boundary condition is used. + for (size_t k : m_kCharge) { + rsd[index(c_offset_Y + k, 0)] = Y(x,k,0) - Y(x,k,1); + } rsd[index(c_offset_P, j)] = m_inletVoltage - phi(x,j); diag[index(c_offset_P, j)] = 0; } else if (j == m_points - 1) { @@ -197,11 +198,7 @@ void IonFlow::evalResidual(double* x, double* rsd, int* diag, // // E = -dV/dz //----------------------------------------------- - double 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; + rsd[index(c_offset_P, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0; diag[index(c_offset_P, j)] = 0; } } @@ -257,13 +254,21 @@ void IonFlow::fixElectricPotential(size_t j) } } -void IonFlow::setElectronTransport(vector_fp& zfixed, vector_fp& diff_e_fixed, - vector_fp& mobi_e_fixed) +void IonFlow::setElectronTransport(vector_fp& tfix, vector_fp& diff_e, + vector_fp& mobi_e) { - m_ztfix = zfixed; - m_diff_e_fix = diff_e_fixed; - m_mobi_e_fix = mobi_e_fixed; m_import_electron_transport = true; + size_t degree = 5; + size_t n = tfix.size(); + vector_fp tlog; + for (size_t i = 0; i < n; i++) { + tlog.push_back(log(tfix[i])); + } + vector_fp w(n, -1.0); + m_diff_e_fix.resize(degree + 1); + m_mobi_e_fix.resize(degree + 1); + polyfit(n, degree, tlog.data(), diff_e.data(), w.data(), m_diff_e_fix.data()); + polyfit(n, degree, tlog.data(), mobi_e.data(), w.data(), m_mobi_e_fix.data()); } void IonFlow::_finalize(const double* x)