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

[Transport] Add class IonGasTransport

Browse source
This commit is contained in:
BangShiuh 2018-04-26 10:49:09 -04:00 • committed by Ray Speth
parent b7e32e4604
commit 3a0f46eb56
8 changed files with 654 additions and 51 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

175
data/inputs/gri30_ion.cti
View file

@ -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
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 N2 AR C3H7
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

33
include/cantera/transport/GasTransport.h
View file

@ -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
/*!

95
include/cantera/transport/IonGasTransport.h Normal file
View file

@ -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<size_t> m_kIon;
//! index of neutral species
std::vector<size_t> m_kNeutral;
//! index of electron
size_t m_kElectron;
//! parameter of omega11 of n64
DenseMatrix m_gamma;
};
}
#endif

2
include/cantera/transport/MixTransport.h
View file

@ -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 *

52
src/transport/GasTransport.cpp
View file

@ -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)
}
}
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++) {

336
src/transport/IonGasTransport.cpp Normal file
View file

@ -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);
}
}
}
}
}
}

4
src/transport/MixTransport.cpp
View file

@ -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)

2
src/transport/TransportFactory.cpp
View file

@ -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(); });