diff --git a/data/inputs/gri30_ion.cti b/data/inputs/gri30_ion.cti index e07cd39ca..933010ee7 100644 --- a/data/inputs/gri30_ion.cti +++ b/data/inputs/gri30_ion.cti @@ -2,21 +2,163 @@ 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 + species=['H2 O2 H2O CH4 CO CO2 N2', + '''gri30: H O OH HO2 H2O2 C CH + CH2 CH2(S) CH3 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 AR C3H7 C3H8 CH2CHO CH3CHO''', 'HCO+ H3O+ E'], reactions=['gri30: all', 'all'], - transport='Mix', + transport='Ion', options=['skip_undeclared_species', 'skip_undeclared_third_bodies'], initial_state=state(temperature=300.0, pressure=OneAtm)) #------------------------------------------------------------------------------- # Species data #------------------------------------------------------------------------------- +# The values of polarizability of H2, O2, H2O, CH4, CO, CO2, and N2 are from +# the supplementary material of Han, Jie, et al. "Numerical modelling of ion +# transport in flames." Combustion Theory and Modelling 19.6 (2015): 744-772. +# DOI: 10.1080/13647830.2015.1090018 + +species(name = "H2", + atoms = " H:2 ", + thermo = ( + NASA( [ 200.00, 1000.00], [ 2.344331120E+00, 7.980520750E-03, + -1.947815100E-05, 2.015720940E-08, -7.376117610E-12, + -9.179351730E+02, 6.830102380E-01] ), + NASA( [ 1000.00, 3500.00], [ 3.337279200E+00, -4.940247310E-05, + 4.994567780E-07, -1.795663940E-10, 2.002553760E-14, + -9.501589220E+02, -3.205023310E+00] ) + ), + transport = gas_transport( + geom = "linear", + diam = 2.92, + well_depth = 38.00, + polar = 0.455, + rot_relax = 280.00), + note = '''The value of polarizability is from the supplementary + material of Han, Jie, et al. "Numerical modelling of ion + transport in flames." Combustion Theory and Modelling + 19.6 (2015): 744-772. DOI: 10.1080/13647830.2015.1090018''' + ) + +species(name = "O2", + atoms = " O:2 ", + thermo = ( + NASA( [ 200.00, 1000.00], [ 3.782456360E+00, -2.996734160E-03, + 9.847302010E-06, -9.681295090E-09, 3.243728370E-12, + -1.063943560E+03, 3.657675730E+00] ), + NASA( [ 1000.00, 3500.00], [ 3.282537840E+00, 1.483087540E-03, + -7.579666690E-07, 2.094705550E-10, -2.167177940E-14, + -1.088457720E+03, 5.453231290E+00] ) + ), + transport = gas_transport( + geom = "linear", + diam = 3.46, + well_depth = 107.40, + polar = 1.131, + rot_relax = 3.80), + note = "TPIS89" + ) + +species(name = "H2O", + atoms = " H:2 O:1 ", + thermo = ( + NASA( [ 200.00, 1000.00], [ 4.198640560E+00, -2.036434100E-03, + 6.520402110E-06, -5.487970620E-09, 1.771978170E-12, + -3.029372670E+04, -8.490322080E-01] ), + NASA( [ 1000.00, 3500.00], [ 3.033992490E+00, 2.176918040E-03, + -1.640725180E-07, -9.704198700E-11, 1.682009920E-14, + -3.000429710E+04, 4.966770100E+00] ) + ), + transport = gas_transport( + geom = "nonlinear", + diam = 2.60, + well_depth = 572.40, + dipole = 1.84, + polar = 1.053, + rot_relax = 4.00), + note = "L 8/89" + ) + +species(name = "CH4", + atoms = " C:1 H:4 ", + thermo = ( + NASA( [ 200.00, 1000.00], [ 5.149876130E+00, -1.367097880E-02, + 4.918005990E-05, -4.847430260E-08, 1.666939560E-11, + -1.024664760E+04, -4.641303760E+00] ), + NASA( [ 1000.00, 3500.00], [ 7.485149500E-02, 1.339094670E-02, + -5.732858090E-06, 1.222925350E-09, -1.018152300E-13, + -9.468344590E+03, 1.843731800E+01] ) + ), + transport = gas_transport( + geom = "nonlinear", + diam = 3.75, + well_depth = 141.40, + polar = 2.60, + rot_relax = 13.00), + note = "L 8/88" + ) + +species(name = "CO", + atoms = " C:1 O:1 ", + thermo = ( + NASA( [ 200.00, 1000.00], [ 3.579533470E+00, -6.103536800E-04, + 1.016814330E-06, 9.070058840E-10, -9.044244990E-13, + -1.434408600E+04, 3.508409280E+00] ), + NASA( [ 1000.00, 3500.00], [ 2.715185610E+00, 2.062527430E-03, + -9.988257710E-07, 2.300530080E-10, -2.036477160E-14, + -1.415187240E+04, 7.818687720E+00] ) + ), + transport = gas_transport( + geom = "linear", + diam = 3.65, + well_depth = 98.10, + polar = 1.95, + rot_relax = 1.80), + note = "TPIS79" + ) + +species(name = "CO2", + atoms = " C:1 O:2 ", + thermo = ( + NASA( [ 200.00, 1000.00], [ 2.356773520E+00, 8.984596770E-03, + -7.123562690E-06, 2.459190220E-09, -1.436995480E-13, + -4.837196970E+04, 9.901052220E+00] ), + NASA( [ 1000.00, 3500.00], [ 3.857460290E+00, 4.414370260E-03, + -2.214814040E-06, 5.234901880E-10, -4.720841640E-14, + -4.875916600E+04, 2.271638060E+00] ) + ), + transport = gas_transport( + geom = "linear", + diam = 3.76, + well_depth = 244.00, + polar = 2.65, + rot_relax = 2.10), + note = "L 7/88" + ) + +species(name = "N2", + atoms = " N:2 ", + thermo = ( + NASA( [ 300.00, 1000.00], [ 3.298677000E+00, 1.408240400E-03, + -3.963222000E-06, 5.641515000E-09, -2.444854000E-12, + -1.020899900E+03, 3.950372000E+00] ), + NASA( [ 1000.00, 5000.00], [ 2.926640000E+00, 1.487976800E-03, + -5.684760000E-07, 1.009703800E-10, -6.753351000E-15, + -9.227977000E+02, 5.980528000E+00] ) + ), + transport = gas_transport( + geom = "linear", + diam = 3.62, + well_depth = 97.53, + polar = 1.76, + rot_relax = 4.00), + note = "121286" + ) species(name = 'HCO+', atoms = ' H:1 C:1 O:1 E:-1 ', @@ -31,9 +173,11 @@ species(name = 'HCO+', transport=gas_transport(geom='linear', diam=3.59, well_depth=498.0, - polar=2.5, - rot_relax=0.0), - note = 'J12/70') + polar=1.356), + note = '''The polarizability is from Han, Jie, et al. + "Numerical modelling of ion transport in flames." + ,and the rest of the parameters are from its neutral + counterpart HCO''') species(name = 'H3O+', atoms = ' H:3 O:1 E:-1 ', @@ -46,12 +190,12 @@ species(name = 'H3O+', 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') + diam=3.15, + well_depth=106.2, + dipole=1.417, + polar=0.897), + note = '''The transport parameters are from Han, Jie, et al. + "Numerical modelling of ion transport in flames."''') species(name = 'E', atoms = ' E:1 ', @@ -66,9 +210,8 @@ species(name = 'E', transport=gas_transport(geom='atom', diam=2.05, well_depth=145.0, - polar=0.667, - rot_relax=0.0), - note = 'gas L10/92') + polar=0.667), + note = 'The transport parameters are not used in IonGasTransport') #------------------------------------------------------------------------------- # Reaction data diff --git a/include/cantera/transport/GasTransport.h b/include/cantera/transport/GasTransport.h index 224205263..997b0276a 100644 --- a/include/cantera/transport/GasTransport.h +++ b/include/cantera/transport/GasTransport.h @@ -144,12 +144,13 @@ protected: //! @name Initialization //! @{ - //! Prepare to build a new kinetic-theory-based transport manager for + //! Setup parameters for a new kinetic-theory-based transport manager for //! low-density gases - /*! - * Uses polynomial fits to Monchick & Mason collision integrals. - */ - void setupMM(); + virtual void setupCollisionParameters(); + + //! Setup range for polynomial fits to collision integrals of + //! Monchick & Mason + void setupCollisionIntegral(); //! Read the transport database /*! @@ -180,29 +181,35 @@ protected: */ void fitCollisionIntegrals(MMCollisionInt& integrals); - //! Generate polynomial fits to the viscosity, conductivity, and - //! the binary diffusion coefficients + //! Generate polynomial fits to the viscosity and conductivity /*! * If CK_mode, then the fits are of the form * \f[ * \log(\eta(i)) = \sum_{n = 0}^3 a_n(i) (\log T)^n * \f] - * and - * \f[ - * \log(D(i,j)) = \sum_{n = 0}^3 a_n(i,j) (\log T)^n - * \f] * Otherwise the fits are of the form * \f[ * \eta(i)/sqrt(k_BT) = \sum_{n = 0}^4 a_n(i) (\log T)^n * \f] - * and + * + * @param integrals interpolator for the collision integrals + */ + virtual void fitProperties(MMCollisionInt& integrals); + + //! Generate polynomial fits to the binary diffusion coefficients + /*! + * If CK_mode, then the fits are of the form + * \f[ + * \log(D(i,j)) = \sum_{n = 0}^3 a_n(i,j) (\log T)^n + * \f] + * Otherwise the fits are of the form * \f[ * D(i,j)/sqrt(k_BT)) = \sum_{n = 0}^4 a_n(i,j) (\log T)^n * \f] * * @param integrals interpolator for the collision integrals */ - void fitProperties(MMCollisionInt& integrals); + virtual void fitDiffCoeffs(MMCollisionInt& integrals); //! Second-order correction to the binary diffusion coefficients /*! diff --git a/include/cantera/transport/IonGasTransport.h b/include/cantera/transport/IonGasTransport.h new file mode 100644 index 000000000..779d97f64 --- /dev/null +++ b/include/cantera/transport/IonGasTransport.h @@ -0,0 +1,95 @@ +/** + * @file IonGasTransport.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. + +#ifndef CT_ION_GAS_TRANSPORT_H +#define CT_ION_GAS_TRANSPORT_H + +#include "MixTransport.h" + +namespace Cantera +{ +//! Class IonGasTransport implements Stockmayer-(n,6,4) model for transport of ions. +/*! + * As implemented here, only binary transport between netrals and ions is considered + * for calculating mixture-average diffusion coefficients and mobilities. When + * polarizability is not provide for an ion, LJ model is used instead of n64 model. + * Only neutral species are considered for thermal conductivity and viscousity. + * + * References for Stockmayer-(n,6,4) model: + * + * 1. Selle, Stefan, and Uwe Riedel. "Transport properties of ionized species." + * Annals of the New York Academy of Sciences 891.1 (1999): 72-80. + * 2. Selle, Stefan, and Uwe Riedel. "Transport coefficients of reacting air at +* high temperatures." 38th Aerospace Sciences Meeting and Exhibit. 1999. + * 3. Han, Jie, et al. "Numerical modelling of ion transport in flames." + * Combustion Theory and Modelling 19.6 (2015): 744-772. + * DOI: 10.1080/13647830.2015.1090018 + * 4. Chiflikian, R. V. "The analog of Blanc’s law for drift velocities + * of electrons in gas mixtures in weakly ionized plasma." + * Physics of Plasmas 2.10 (1995): 3902-3909. + * 5. Viehland, L. A., et al. "Tables of transport collision integrals for + * (n, 6, 4) ion-neutral potentials." Atomic Data and Nuclear Data Tables + * 16.6 (1975): 495-514. + * @ingroup tranprops + */ +class IonGasTransport : public MixTransport +{ +public: + IonGasTransport(); + + virtual std::string transportType() const { + return "Ion"; + } + + virtual void init(thermo_t* thermo, int mode, int log_level); + + //! Viscosity of the mixture (kg/m/s). + //! Only Neutral species contribute to Viscosity. + virtual double viscosity(); + + //! Returns the mixture thermal conductivity (W/m/K). + //! Only Neutral species contribute to therrmal conductivity. + virtual double thermalConductivity(); + +protected: + //! setup parameters for n64 model + void setupN64(); + + //! Generate polynomial fits to the binary diffusion coefficients. + //! Use Stockmayer-(n,6,4) model for collision between charged and neutral species. + virtual void fitDiffCoeffs(MMCollisionInt& integrals); + + /*! + * Collision integral of omega11 of n64 collision model. + * The collision integral was fitted by Han et al. using the table + * by Viehlan et al. + * Note: Han release the range to 1000, but Selle suggested that + * a high temperature model is needed for T* > 10. + */ + double omega11_n64(const double tstar, const double gamma); + + virtual void getMixDiffCoeffs(doublereal* const d); + + //! electrical properties + vector_int m_speciesCharge; + + //! index of ions (exclude electron.) + std::vector m_kIon; + + //! index of neutral species + std::vector m_kNeutral; + + //! index of electron + size_t m_kElectron; + + //! parameter of omega11 of n64 + DenseMatrix m_gamma; +}; + +} + +#endif diff --git a/include/cantera/transport/MixTransport.h b/include/cantera/transport/MixTransport.h index 0c8e0449f..79fed77e5 100644 --- a/include/cantera/transport/MixTransport.h +++ b/include/cantera/transport/MixTransport.h @@ -151,7 +151,7 @@ public: virtual void init(thermo_t* thermo, int mode=0, int log_level=0); -private: +protected: //! Calculate the pressure from the ideal gas law doublereal pressure_ig() const { return (m_thermo->molarDensity() * GasConstant * diff --git a/src/transport/GasTransport.cpp b/src/transport/GasTransport.cpp index b9e7aae61..277b4dbc9 100644 --- a/src/transport/GasTransport.cpp +++ b/src/transport/GasTransport.cpp @@ -277,8 +277,10 @@ void GasTransport::init(thermo_t* thermo, int mode, int log_level) m_nsp = m_thermo->nSpecies(); m_mode = mode; m_log_level = log_level; + // set up Monchick and Mason collision integrals - setupMM(); + setupCollisionParameters(); + setupCollisionIntegral(); m_molefracs.resize(m_nsp); m_spwork.resize(m_nsp); @@ -299,15 +301,9 @@ void GasTransport::init(thermo_t* thermo, int mode, int log_level) m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]); } } - - // set flags all false - m_visc_ok = false; - m_viscwt_ok = false; - m_spvisc_ok = false; - m_bindiff_ok = false; } -void GasTransport::setupMM() +void GasTransport::setupCollisionParameters() { m_epsilon.resize(m_nsp, m_nsp, 0.0); m_delta.resize(m_nsp, m_nsp, 0.0); @@ -332,7 +328,6 @@ void GasTransport::setupMM() m_poly[i].resize(m_nsp); } - double tstar_min = 1.e8, tstar_max = 0.0; double f_eps, f_sigma; for (size_t i = 0; i < m_nsp; i++) { @@ -346,11 +341,6 @@ void GasTransport::setupMM() // the effective well depth for (i,j) collisions m_epsilon(i,j) = sqrt(m_eps[i]*m_eps[j]); - // The polynomial fits of collision integrals vs. T* - // will be done for the T* from tstar_min to tstar_max - tstar_min = std::min(tstar_min, Boltzmann * m_thermo->minTemp()/m_epsilon(i,j)); - tstar_max = std::max(tstar_max, Boltzmann * m_thermo->maxTemp()/m_epsilon(i,j)); - // the effective dipole moment for (i,j) collisions m_dipole(i,j) = sqrt(m_dipole(i,i)*m_dipole(j,j)); @@ -370,7 +360,19 @@ void GasTransport::setupMM() m_delta(j,i) = m_delta(i,j); } } +} +void GasTransport::setupCollisionIntegral() +{ + double tstar_min = 1.e8, tstar_max = 0.0; + for (size_t i = 0; i < m_nsp; i++) { + for (size_t j = i; j < m_nsp; j++) { + // The polynomial fits of collision integrals vs. T* + // will be done for the T* from tstar_min to tstar_max + tstar_min = std::min(tstar_min, Boltzmann * m_thermo->minTemp()/m_epsilon(i,j)); + tstar_max = std::max(tstar_max, Boltzmann * m_thermo->maxTemp()/m_epsilon(i,j)); + } + } // Chemkin fits the entire T* range in the Monchick and Mason tables, // so modify tstar_min and tstar_max if in Chemkin compatibility mode if (m_mode == CK_Mode) { @@ -662,7 +664,29 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) } } - mxerr = 0.0, mxrelerr = 0.0; + fitDiffCoeffs(integrals); +} + +void GasTransport::fitDiffCoeffs(MMCollisionInt& integrals) +{ + // number of points to use in generating fit data + const size_t np = 50; + int degree = (m_mode == CK_Mode ? 3 : 4); + double dt = (m_thermo->maxTemp() - m_thermo->minTemp())/(np-1); + vector_fp tlog(np); + vector_fp w(np), w2(np); + + // generate array of log(t) values + for (size_t n = 0; n < np; n++) { + double t = m_thermo->minTemp() + dt*n; + tlog[n] = log(t); + } + + // vector of polynomial coefficients + vector_fp c(degree + 1), c2(degree + 1); + double err, relerr, + mxerr = 0.0, mxrelerr = 0.0; + vector_fp diff(np + 1); m_diffcoeffs.clear(); for (size_t k = 0; k < m_nsp; k++) { diff --git a/src/transport/IonGasTransport.cpp b/src/transport/IonGasTransport.cpp new file mode 100644 index 000000000..5b2db60b8 --- /dev/null +++ b/src/transport/IonGasTransport.cpp @@ -0,0 +1,336 @@ +//! @file IonGasTransport.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/transport/IonGasTransport.h" +#include "cantera/numerics/polyfit.h" +#include "cantera/base/stringUtils.h" +#include "MMCollisionInt.h" + +namespace Cantera +{ +IonGasTransport::IonGasTransport() : + m_kElectron(npos) +{ +} + +void IonGasTransport::init(thermo_t* thermo, int mode, int log_level) +{ + m_thermo = thermo; + m_nsp = m_thermo->nSpecies(); + m_mode = mode; + if (m_mode == CK_Mode) { + throw CanteraError("IonGasTransport::init(thermo, mode, log_level)", + "mode = CK_Mode, which is an outdated lower-order fit."); + } + m_log_level = log_level; + // 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 the index of electron + if (m_thermo->speciesIndex("E") != npos ) { + m_kElectron = m_thermo->speciesIndex("E"); + } + + // Find indices for charge of species + for (size_t k = 0; k < m_nsp; k++) { + if (m_speciesCharge[k] != 0){ + if (k != m_kElectron) { + m_kIon.push_back(k); + } + } else { + m_kNeutral.push_back(k); + } + } + // set up Monchick and Mason parameters + setupCollisionParameters(); + // set up n64 parameters + setupN64(); + // setup collision integrals + setupCollisionIntegral(); + m_molefracs.resize(m_nsp); + m_spwork.resize(m_nsp); + m_visc.resize(m_nsp); + m_sqvisc.resize(m_nsp); + m_phi.resize(m_nsp, m_nsp, 0.0); + m_bdiff.resize(m_nsp, m_nsp); + m_cond.resize(m_nsp); + + // make a local copy of the molecular weights + m_mw = m_thermo->molecularWeights(); + + m_wratjk.resize(m_nsp, m_nsp, 0.0); + m_wratkj1.resize(m_nsp, m_nsp, 0.0); + for (size_t j = 0; j < m_nsp; j++) { + for (size_t k = j; k < m_nsp; k++) { + m_wratjk(j,k) = sqrt(m_mw[j]/m_mw[k]); + m_wratjk(k,j) = sqrt(m_wratjk(j,k)); + m_wratkj1(j,k) = sqrt(1.0 + m_mw[k]/m_mw[j]); + } + } +} + +double IonGasTransport::viscosity() +{ + update_T(); + update_C(); + + if (m_visc_ok) { + return m_viscmix; + } + + double vismix = 0.0; + // update m_visc and m_phi if necessary + if (!m_viscwt_ok) { + updateViscosity_T(); + } + + multiply(m_phi, m_molefracs.data(), m_spwork.data()); + + for (size_t k : m_kNeutral) { + vismix += m_molefracs[k] * m_visc[k]/m_spwork[k]; //denom; + } + m_viscmix = vismix; + return vismix; +} + +double IonGasTransport::thermalConductivity() +{ + update_T(); + update_C(); + if (!m_spcond_ok) { + updateCond_T(); + } + if (!m_condmix_ok) { + doublereal sum1 = 0.0, sum2 = 0.0; + for (size_t k : m_kNeutral) { + sum1 += m_molefracs[k] * m_cond[k]; + sum2 += m_molefracs[k] / m_cond[k]; + } + m_lambda = 0.5*(sum1 + 1.0/sum2); + m_condmix_ok = true; + } + return m_lambda; +} + +void IonGasTransport::fitDiffCoeffs(MMCollisionInt& integrals) +{ + GasTransport::fitDiffCoeffs(integrals); + + // number of points to use in generating fit data + const size_t np = 50; + int degree = 4; + double dt = (m_thermo->maxTemp() - m_thermo->minTemp())/(np-1); + vector_fp tlog(np); + vector_fp w(np); + + // generate array of log(t) values + for (size_t n = 0; n < np; n++) { + double t = m_thermo->minTemp() + dt*n; + tlog[n] = log(t); + } + + // vector of polynomial coefficients + vector_fp c(degree + 1); + double err = 0.0, relerr = 0.0, + mxerr = 0.0, mxrelerr = 0.0; + + vector_fp diff(np + 1); + // The array order still not ideal + for (size_t k = 0; k < m_nsp; k++) { + for (size_t j = k; j < m_nsp; j++) { + if (m_alpha[k] == 0.0 || m_alpha[j] == 0.0 || + k == m_kElectron || j == m_kElectron) { + continue; + } + if (m_speciesCharge[k] == 0) { + if (m_speciesCharge[j] == 0) { + continue; + } + } else { + if (m_speciesCharge[j] != 0) { + continue; + } + } + for (size_t n = 0; n < np; n++) { + double t = m_thermo->minTemp() + dt*n; + double eps = m_epsilon(j,k); + double tstar = Boltzmann * t/eps; + double sigma = m_diam(j,k); + double om11 = omega11_n64(tstar, m_gamma(j,k)); + double diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,j)) + * pow(Boltzmann * t, 1.5) / (Pi * sigma * sigma * om11); + + diff[n] = diffcoeff/pow(t, 1.5); + w[n] = 1.0/(diff[n]*diff[n]); + } + polyfit(np, degree, tlog.data(), diff.data(), w.data(), c.data()); + + for (size_t n = 0; n < np; n++) { + double val, fit; + double t = exp(tlog[n]); + double pre = pow(t, 1.5); + val = pre * diff[n]; + fit = pre * poly4(tlog[n], c.data()); + err = fit - val; + relerr = err/val; + mxerr = std::max(mxerr, fabs(err)); + mxrelerr = std::max(mxrelerr, fabs(relerr)); + } + size_t sum = k * (k + 1) / 2; + m_diffcoeffs[k*m_nsp+j-sum] = c; + if (m_log_level >= 2) { + writelog(m_thermo->speciesName(k) + "__" + + m_thermo->speciesName(j) + ": [" + vec2str(c) + "]\n"); + } + } + } + + if (m_log_level) { + writelogf("Maximum binary diffusion coefficient absolute error:" + " %12.6g\n", mxerr); + writelogf("Maximum binary diffusion coefficient relative error:" + "%12.6g", mxrelerr); + } +} + +void IonGasTransport::setupN64() +{ + m_gamma.resize(m_nsp, m_nsp, 0.0); + for (size_t i : m_kIon) { + for (size_t j : m_kNeutral) { + if (m_alpha[j] != 0.0 && m_alpha[i] != 0.0) { + double r_alpha = m_alpha[i] / m_alpha[j]; + // save a copy of polarizability in Angstrom + double alphaA_i = m_alpha[i] * 1e30; + double alphaA_j = m_alpha[j] * 1e30; + // The ratio of dispersion to induction forces + double xi = alphaA_i / (m_speciesCharge[i] * m_speciesCharge[i] * + (1.0 + pow(2 * r_alpha, 2./3.)) * sqrt(alphaA_j)); + + // the collision diameter + double K1 = 1.767; + double kappa = 0.095; + m_diam(i,j) = K1 * (pow(m_alpha[i], 1./3.) + pow(m_alpha[j], 1./3.)) / + pow(alphaA_i * alphaA_j * (1.0 + 1.0 / xi), kappa); + + // The original K2 is 0.72, but Han et al. suggested that K2 = 1.44 + // for better fit. + double K2 = 1.44; + double epsilon = K2 * ElectronCharge * ElectronCharge * + m_speciesCharge[i] * m_speciesCharge[i] * + m_alpha[j] * (1.0 + xi) / + (8 * Pi * epsilon_0 * pow(m_diam(i,j),4)); + if (epsilon != 0.0) { + m_epsilon(i,j) = epsilon; + } + + // Calculate dipersion coefficient and quadrupole polarizability + // from curve fitting if not available. + // Neutrals + if (m_disp[j] == 0.0) { + m_disp[j] = exp(1.8846*log(alphaA_j)-0.4737)* 1e-50; + } + if (m_quad_polar[j] == 0.0) { + m_quad_polar[j] = 2.0 * m_disp[j]; + } + // Ions + if (m_disp[i] == 0.0) { + if (m_speciesCharge[i] > 0) { + m_disp[i] = exp(1.8853*log(alphaA_i)+0.2682)* 1e-50; + } else { + m_disp[i] = exp(3.2246*log(alphaA_i)-3.2397)* 1e-50; + } + } + + // The binary dispersion coefficient is determined by the combination rule + // Reference: + // Tang, K. T. "Dynamic polarizabilities and van der Waals coefficients." + // Physical Review 177.1 (1969): 108. + double C6 = 2.0 * m_disp[i] * m_disp[j] / + (1.0/r_alpha * m_disp[i] + r_alpha * m_disp[j]); + + m_gamma(i,j) = (2.0 / pow(m_speciesCharge[i],2) * C6 + m_quad_polar[j]) / + (m_alpha[j] * m_diam(i,j) * m_diam(i,j));//Dimensionless + + // properties are symmetric + m_diam(j,i) = m_diam(i,j); + m_epsilon(j,i) = m_epsilon(i,j); + m_gamma(j,i) = m_gamma(i,j); + } + } + } +} + +double IonGasTransport::omega11_n64(const double tstar, const double gamma) +{ + double logtstar = log(tstar); + double om11 = 0.0; + if (tstar < 0.01) { + throw CanteraError("IonGasTransport::omega11_n64(tstar, gamma)", + "tstar = {} is smaller than 0.01", tstar); + } else if (tstar <= 0.04) { + // for interval 0.01 to 0.04, SSE = 0.006; R^2 = 1; RMSE = 0.020 + om11 = 2.97 - 12.0 * gamma + - 0.887 * logtstar + + 3.86 * gamma * gamma + - 6.45 * gamma * logtstar + - 0.275 * logtstar * logtstar + + 1.20 * gamma * gamma * logtstar + - 1.24 * gamma * logtstar * logtstar + - 0.164 * pow(logtstar,3); + } else if (tstar <= 1000) { + // for interval 0.04 to 1000, SSE = 0.282; R^2 = 1; RMSE = 0.033 + om11 = 1.22 - 0.0343 * gamma + + (-0.769 + 0.232 * gamma) * logtstar + + (0.306 - 0.165 * gamma) * logtstar * logtstar + + (-0.0465 + 0.0388 * gamma) * pow(logtstar,3) + + (0.000614 - 0.00285 * gamma) * pow(logtstar,4) + + 0.000238 * pow(logtstar,5); + } else { + throw CanteraError("IonGasTransport::omega11_n64(tstar, gamma)", + "tstar = {} is larger than 1000", tstar); + } + return om11; +} + +void IonGasTransport::getMixDiffCoeffs(double* const d) +{ + update_T(); + update_C(); + + // update the binary diffusion coefficients if necessary + if (!m_bindiff_ok) { + updateDiff_T(); + } + + double mmw = m_thermo->meanMolecularWeight(); + double p = m_thermo->pressure(); + if (m_nsp == 1) { + d[0] = m_bdiff(0,0) / p; + } else { + for (size_t k = 0; k < m_nsp; k++) { + if (k == m_kElectron) { + d[k] = 0.4 * m_kbt / ElectronCharge; + } else { + double sum2 = 0.0; + for (size_t j : m_kNeutral) { + if (j != k) { + sum2 += m_molefracs[j] / m_bdiff(j,k); + } + } + if (sum2 <= 0.0) { + d[k] = m_bdiff(k,k) / p; + } else { + d[k] = (mmw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2); + } + } + } + } +} + +} + diff --git a/src/transport/MixTransport.cpp b/src/transport/MixTransport.cpp index 9eae08ee4..066f4a555 100644 --- a/src/transport/MixTransport.cpp +++ b/src/transport/MixTransport.cpp @@ -24,10 +24,6 @@ void MixTransport::init(ThermoPhase* thermo, int mode, int log_level) { GasTransport::init(thermo, mode, log_level); m_cond.resize(m_nsp); - - // set flags all false - m_spcond_ok = false; - m_condmix_ok = false; } void MixTransport::getMobilities(doublereal* const mobil) diff --git a/src/transport/TransportFactory.cpp b/src/transport/TransportFactory.cpp index a81dd15cf..dc193f5d7 100644 --- a/src/transport/TransportFactory.cpp +++ b/src/transport/TransportFactory.cpp @@ -7,6 +7,7 @@ #include "cantera/transport/MultiTransport.h" #include "cantera/transport/MixTransport.h" #include "cantera/transport/UnityLewisTransport.h" +#include "cantera/transport/IonGasTransport.h" #include "cantera/transport/SolidTransport.h" #include "cantera/transport/DustyGasTransport.h" #include "cantera/transport/SimpleTransport.h" @@ -50,6 +51,7 @@ TransportFactory::TransportFactory() reg("UnityLewis", []() { return new UnityLewisTransport(); }); reg("Mix", []() { return new MixTransport(); }); reg("Multi", []() { return new MultiTransport(); }); + reg("Ion", []() { return new IonGasTransport(); }); m_synonyms["CK_Mix"] = "Mix"; m_synonyms["CK_Multi"] = "Multi"; reg("HighP", []() { return new HighPressureGasTransport(); });