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

[Transport] Make better use of local variables

Browse source
This commit is contained in:
Ray Speth 2016-06-27 17:11:00 -04:00
parent fd4cbb8718
commit a8c6fe09d9
12 changed files with 196 additions and 342 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

11
src/transport/DustyGasTransport.cpp
View file

@ -136,9 +136,8 @@ void DustyGasTransport::updateKnudsenDiffCoeffs()
return;
}
doublereal K_g = m_pore_radius * m_porosity / m_tortuosity;
const doublereal TwoThirds = 2.0/3.0;
for (size_t k = 0; k < m_nsp; k++) {
m_dk[k] = TwoThirds * K_g * sqrt((8.0 * GasConstant * m_temp)/
m_dk[k] = 2.0/3.0 * K_g * sqrt((8.0 * GasConstant * m_temp)/
(Pi * m_mw[k]));
}
m_knudsen_ok = true;
@ -148,7 +147,6 @@ void DustyGasTransport::eval_H_matrix()
{
updateBinaryDiffCoeffs();
updateKnudsenDiffCoeffs();
doublereal sum;
for (size_t k = 0; k < m_nsp; k++) {
// evaluate off-diagonal terms
for (size_t j = 0; j < m_nsp; j++) {
@ -156,7 +154,7 @@ void DustyGasTransport::eval_H_matrix()
}
// evaluate diagonal term
sum = 0.0;
double sum = 0.0;
for (size_t j = 0; j < m_nsp; j++) {
if (j != k) {
sum += m_x[j]/m_d(k,j);
@ -171,7 +169,6 @@ void DustyGasTransport::getMolarFluxes(const doublereal* const state1,
const doublereal delta,
doublereal* const fluxes)
{
doublereal conc1, conc2;
// cbar will be the average concentration between the two points
doublereal* const cbar = m_spwork.data();
doublereal* const gradc = m_spwork2.data();
@ -184,8 +181,8 @@ void DustyGasTransport::getMolarFluxes(const doublereal* const state1,
doublereal c1sum = 0.0, c2sum = 0.0;
for (size_t k = 0; k < m_nsp; k++) {
conc1 = rho1 * y1[k] / m_mw[k];
conc2 = rho2 * y2[k] / m_mw[k];
double conc1 = rho1 * y1[k] / m_mw[k];
double conc2 = rho2 * y2[k] / m_mw[k];
cbar[k] = 0.5*(conc1 + conc2);
gradc[k] = (conc2 - conc1) / delta;
c1sum += conc1;

67
src/transport/GasTransport.cpp
View file

@ -154,7 +154,6 @@ doublereal GasTransport::viscosity()
void GasTransport::updateViscosity_T()
{
doublereal vratiokj, wratiojk, factor1;
if (!m_spvisc_ok) {
updateSpeciesViscosities();
}
@ -162,11 +161,11 @@ void GasTransport::updateViscosity_T()
// see Eq. (9-5.15) of Reid, Prausnitz, and Poling
for (size_t j = 0; j < m_nsp; j++) {
for (size_t k = j; k < m_nsp; k++) {
vratiokj = m_visc[k]/m_visc[j];
wratiojk = m_mw[j]/m_mw[k];
double vratiokj = m_visc[k]/m_visc[j];
double wratiojk = m_mw[j]/m_mw[k];
// Note that m_wratjk(k,j) holds the square root of m_wratjk(j,k)!
factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
double factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j);
m_phi(k,j) = factor1*factor1 / (sqrt(8.0) * m_wratkj1(j,k));
m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
}
@ -507,8 +506,6 @@ void GasTransport::makePolarCorrections(size_t i, size_t j,
void GasTransport::fitCollisionIntegrals(MMCollisionInt& integrals)
{
double dstar;
// Chemkin fits to sixth order polynomials
int degree = (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE);
if (m_log_level) {
@ -523,11 +520,7 @@ void GasTransport::fitCollisionIntegrals(MMCollisionInt& integrals)
for (size_t i = 0; i < m_nsp; i++) {
for (size_t j = i; j < m_nsp; j++) {
// Chemkin fits only delta* = 0
if (m_mode != CK_Mode) {
dstar = m_delta(i,j);
} else {
dstar = 0.0;
}
double dstar = (m_mode != CK_Mode) ? m_delta(i,j) : 0.0;
// if a fit has already been generated for delta* = m_delta(i,j),
// then use it. Otherwise, make a new fit, and add m_delta(i,j) to
@ -577,7 +570,7 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
if (m_log_level && m_log_level < 2) {
writelog("*** polynomial coefficients not printed (log_level < 2) ***\n");
}
double sqrt_T, visc, err, relerr,
double visc, err, relerr,
mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0;
if (m_log_level) {
@ -590,9 +583,6 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
}
}
double cp_R, cond, w_RT, f_int, A_factor, B_factor, c1, cv_rot, cv_int,
f_rot, f_trans, om11, diffcoeff;
const vector_fp& mw = m_thermo->molecularWeights();
for (size_t k = 0; k < m_nsp; k++) {
for (size_t n = 0; n < np; n++) {
@ -600,32 +590,31 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
m_thermo->setTemperature(t);
vector_fp cp_R_all(m_thermo->nSpecies());
m_thermo->getCp_R_ref(&cp_R_all[0]);
cp_R = cp_R_all[k];
double cp_R = cp_R_all[k];
double tstar = Boltzmann * t/ m_eps[k];
sqrt_T = sqrt(t);
double sqrt_T = sqrt(t);
double om22 = integrals.omega22(tstar, m_delta(k,k));
om11 = integrals.omega11(tstar, m_delta(k,k));
double om11 = integrals.omega11(tstar, m_delta(k,k));
// self-diffusion coefficient, without polar corrections
diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,k)) *
pow((Boltzmann * t), 1.5)/
(Pi * m_sigma[k] * m_sigma[k] * om11);
double diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,k)) *
pow((Boltzmann * t), 1.5)/
(Pi * m_sigma[k] * m_sigma[k] * om11);
// viscosity
visc = 5.0/16.0 * sqrt(Pi * mw[k] * Boltzmann * t / Avogadro) /
(om22 * Pi * m_sigma[k]*m_sigma[k]);
// thermal conductivity
w_RT = mw[k]/(GasConstant * t);
f_int = w_RT * diffcoeff/visc;
cv_rot = m_crot[k];
A_factor = 2.5 - f_int;
B_factor = m_zrot[k] + 2.0/Pi * (5.0/3.0 * cv_rot + f_int);
c1 = 2.0/Pi * A_factor/B_factor;
cv_int = cp_R - 2.5 - cv_rot;
f_rot = f_int * (1.0 + c1);
f_trans = 2.5 * (1.0 - c1 * cv_rot/1.5);
cond = (visc/mw[k])*GasConstant*(f_trans * 1.5
double f_int = mw[k]/(GasConstant * t) * diffcoeff/visc;
double cv_rot = m_crot[k];
double A_factor = 2.5 - f_int;
double B_factor = m_zrot[k] + 2.0/Pi * (5.0/3.0 * cv_rot + f_int);
double c1 = 2.0/Pi * A_factor/B_factor;
double cv_int = cp_R - 2.5 - cv_rot;
double f_rot = f_int * (1.0 + c1);
double f_trans = 2.5 * (1.0 - c1 * cv_rot/1.5);
double cond = (visc/mw[k])*GasConstant*(f_trans * 1.5
+ f_rot * cv_rot + f_int * cv_int);
if (m_mode == CK_Mode) {
@ -660,7 +649,7 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
val = exp(spvisc[n]);
fit = exp(poly3(tlog[n], c.data()));
} else {
sqrt_T = exp(0.5*tlog[n]);
double sqrt_T = exp(0.5*tlog[n]);
val = sqrt_T * pow(spvisc[n],2);
fit = sqrt_T * pow(poly4(tlog[n], c.data()),2);
}
@ -677,7 +666,7 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
val = exp(spcond[n]);
fit = exp(poly3(tlog[n], c2.data()));
} else {
sqrt_T = exp(0.5*tlog[n]);
double sqrt_T = exp(0.5*tlog[n]);
val = sqrt_T * spcond[n];
fit = sqrt_T * poly4(tlog[n], c2.data());
}
@ -723,18 +712,16 @@ void GasTransport::fitProperties(MMCollisionInt& integrals)
mxerr = 0.0, mxrelerr = 0.0;
vector_fp diff(np + 1);
double eps, sigma;
for (size_t k = 0; k < m_nsp; k++) {
for (size_t j = k; j < m_nsp; j++) {
for (size_t n = 0; n < np; n++) {
double t = m_thermo->minTemp() + dt*n;
eps = m_epsilon(j,k);
double eps = m_epsilon(j,k);
double tstar = Boltzmann * t/eps;
sigma = m_diam(j,k);
om11 = integrals.omega11(tstar, m_delta(j,k));
diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,j)) *
pow(Boltzmann * t, 1.5) /
(Pi * sigma * sigma * om11);
double sigma = m_diam(j,k);
double om11 = integrals.omega11(tstar, m_delta(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);
// 2nd order correction
// NOTE: THIS CORRECTION IS NOT APPLIED

57
src/transport/HighPressureGasTransport.cpp
View file

@ -132,7 +132,6 @@ void HighPressureGasTransport::getThermalDiffCoeffs(doublereal* const dt)
void HighPressureGasTransport::getBinaryDiffCoeffs(const size_t ld, doublereal* const d)
{
doublereal P_corr_ij, Tr_ij, Pr_ij;
vector_fp PcP(5);
size_t nsp = m_thermo->nSpecies();
vector_fp molefracs(nsp);
@ -161,9 +160,10 @@ void HighPressureGasTransport::getBinaryDiffCoeffs(const size_t ld, doublereal*
x_j = x_j/(x_i + x_j);
//Calculate Tr and Pr based on mole-fraction-weighted crit constants:
Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j));
Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j));
double Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j));
double Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j));
double P_corr_ij;
if (Pr_ij < 0.1) {
// If pressure is low enough, no correction is needed:
P_corr_ij = 1;
@ -204,7 +204,6 @@ void HighPressureGasTransport::getMultiDiffCoeffs(const size_t ld, doublereal* c
// Correct the binary diffusion coefficients for high-pressure effects; this
// is basically the same routine used in 'getBinaryDiffCoeffs,' above:
doublereal P_corr_ij, Tr_ij, Pr_ij;
size_t nsp = m_thermo->nSpecies();
vector_fp molefracs(nsp);
m_thermo->getMoleFractions(&molefracs[0]);
@ -225,9 +224,10 @@ void HighPressureGasTransport::getMultiDiffCoeffs(const size_t ld, doublereal* c
doublereal x_j = std::max(Tiny, molefracs[j]);
x_i = x_i/(x_i+x_j);
x_j = x_j/(x_i+x_j);
Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j));
Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j));
double Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j));
double Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j));
double P_corr_ij;
if (Pr_ij < 0.1) {
P_corr_ij = 1;
}else {
@ -260,14 +260,12 @@ void HighPressureGasTransport::getMultiDiffCoeffs(const size_t ld, doublereal* c
m_l0000_ok = false; // matrix is overwritten by inverse
m_lmatrix_soln_ok = false;
doublereal pres = m_thermo->pressure();
doublereal prefactor = 16.0 * m_temp
*m_thermo->meanMolecularWeight()/(25.0*pres);
doublereal c;
*m_thermo->meanMolecularWeight()/(25.0*m_thermo->pressure());
for (size_t i = 0; i < m_nsp; i++) {
for (size_t j = 0; j < m_nsp; j++) {
c = prefactor/m_mw[j];
double c = prefactor/m_mw[j];
d[ld*j + i] = c*molefracs[i]*(m_Lmatrix(i,j) - m_Lmatrix(i,i));
}
}
@ -280,7 +278,6 @@ doublereal HighPressureGasTransport::viscosity()
double Pc_mix_n = 0.;
double Pc_mix_d = 0.;
double MW_mix = m_thermo->meanMolecularWeight();
doublereal x_H, Tc, Zc, Tr, Afac, Z1m, Z2m;
double MW_H = m_mw[0];
double MW_L = m_mw[0];
doublereal FP_mix_o = 0;
@ -291,13 +288,13 @@ doublereal HighPressureGasTransport::viscosity()
vector_fp molefracs(nsp);
m_thermo->getMoleFractions(&molefracs[0]);
x_H = molefracs[0];
double x_H = molefracs[0];
for (size_t i = 0; i < m_nsp; i++) {
// Calculate pure-species critical constants and add their contribution
// to the mole-fraction-weighted mixture averages:
Tc = Tcrit_i(i);
Tr = tKelvin/Tc;
Zc = Zcrit_i(i);
double Tc = Tcrit_i(i);
double Tr = tKelvin/Tc;
double Zc = Zcrit_i(i);
Tc_mix += Tc*molefracs[i];
Pc_mix_n += molefracs[i]*Zc; //numerator
Pc_mix_d += molefracs[i]*Vcrit_i(i); //denominator
@ -345,17 +342,15 @@ doublereal HighPressureGasTransport::viscosity()
*pow(Pc_mix,4)),1.0/6.0);
if (ratio > 9 && x_H > 0.05 && x_H < 0.7) {
Afac = 1 - 0.01*pow(ratio,0.87);
} else {
Afac = 1;
FQ_mix_o *= 1 - 0.01*pow(ratio,0.87);
}
FQ_mix_o *= Afac;
// Calculate Z1m
Z1m = (0.807*pow(Tr_mix,0.618) - 0.357*exp(-0.449*Tr_mix)
+ 0.340*exp(-4.058*Tr_mix)+0.018)*FP_mix_o*FQ_mix_o;
double Z1m = (0.807*pow(Tr_mix,0.618) - 0.357*exp(-0.449*Tr_mix)
+ 0.340*exp(-4.058*Tr_mix)+0.018)*FP_mix_o*FQ_mix_o;
// Calculate Z2m:
double Z2m;
if (Tr_mix <= 1.0) {
if (Pr_mix < Pvp_mix/Pc_mix) {
doublereal alpha = 3.262 + 14.98*pow(Pr_mix,5.508);
@ -396,9 +391,8 @@ doublereal HighPressureGasTransport::viscosity()
// Pure species critical properties - Tc, Pc, Vc, Zc:
doublereal HighPressureGasTransport::Tcrit_i(size_t i)
{
size_t nsp = m_thermo->nSpecies();
// Store current molefracs and set temp molefrac of species i to 1.0:
vector_fp molefracs = store(i,nsp);
vector_fp molefracs = store(i, m_thermo->nSpecies());
double tc = m_thermo->critTemperature();
// Restore actual molefracs:
@ -408,9 +402,8 @@ doublereal HighPressureGasTransport::Tcrit_i(size_t i)
doublereal HighPressureGasTransport::Pcrit_i(size_t i)
{
size_t nsp = m_thermo->nSpecies();
// Store current molefracs and set temp molefrac of species i to 1.0:
vector_fp molefracs = store(i,nsp);
vector_fp molefracs = store(i, m_thermo->nSpecies());
double pc = m_thermo->critPressure();
// Restore actual molefracs:
@ -420,9 +413,8 @@ doublereal HighPressureGasTransport::Pcrit_i(size_t i)
doublereal HighPressureGasTransport::Vcrit_i(size_t i)
{
size_t nsp = m_thermo->nSpecies();
// Store current molefracs and set temp molefrac of species i to 1.0:
vector_fp molefracs = store(i,nsp);
vector_fp molefracs = store(i, m_thermo->nSpecies());
double vc = m_thermo->critVolume();
// Restore actual molefracs:
@ -432,9 +424,8 @@ doublereal HighPressureGasTransport::Vcrit_i(size_t i)
doublereal HighPressureGasTransport::Zcrit_i(size_t i)
{
size_t nsp = m_thermo->nSpecies();
// Store current molefracs and set temp molefrac of species i to 1.0:
vector_fp molefracs = store(i,nsp);
vector_fp molefracs = store(i, m_thermo->nSpecies());
double zc = m_thermo->critCompressibility();
// Restore actual molefracs:
@ -446,12 +437,8 @@ vector_fp HighPressureGasTransport::store(size_t i, size_t nsp)
{
vector_fp molefracs(nsp);
m_thermo->getMoleFractions(&molefracs[0]);
vector_fp mf_temp(nsp);
for (size_t j = 0; j < nsp; j++) {
if (j == i) {
mf_temp[j] = 1;
} else {mf_temp[j] = 0;}
}
vector_fp mf_temp(nsp, 0.0);
mf_temp[i] = 1;
m_thermo->setMoleFractions(&mf_temp[0]);
return molefracs;
}

4
src/transport/LTPspecies.cpp
View file

@ -16,7 +16,6 @@ namespace Cantera
class LTPError : public CanteraError
{
public:
//! Constructor is a wrapper around CanteraError
/*!
* @param msg Informative message
@ -44,8 +43,7 @@ static void getArrhenius(const XML_Node& node,
// parse the children for the A, b, and E components.
A = getFloat(node, "A", "toSI");
b = getFloat(node, "b");
E = getFloat(node, "E", "actEnergy");
E /= GasConstant;
E = getFloat(node, "E", "actEnergy") / GasConstant;
}
LTPspecies::LTPspecies(const XML_Node* const propNode, const std::string name,

76
src/transport/LiquidTranInteraction.cpp
View file

@ -35,20 +35,16 @@ LiquidTranInteraction::LiquidTranInteraction(TransportPropertyType tp_ind) :
LiquidTranInteraction::~LiquidTranInteraction()
{
size_t kmax = m_Aij.size();
for (size_t k = 0; k < kmax; k++) {
for (size_t k = 0; k < m_Aij.size(); k++) {
delete m_Aij[k];
}
kmax = m_Bij.size();
for (size_t k = 0; k < kmax; k++) {
for (size_t k = 0; k < m_Bij.size(); k++) {
delete m_Bij[k];
}
kmax = m_Hij.size();
for (size_t k = 0; k < kmax; k++) {
for (size_t k = 0; k < m_Hij.size(); k++) {
delete m_Hij[k];
}
kmax = m_Sij.size();
for (size_t k = 0; k < kmax; k++) {
for (size_t k = 0; k < m_Sij.size(); k++) {
delete m_Sij[k];
}
}
@ -60,19 +56,16 @@ void LiquidTranInteraction::init(const XML_Node& compModelNode,
size_t nsp = thermo->nSpecies();
m_Dij.resize(nsp, nsp, 0.0);
m_Eij.resize(nsp, nsp, 0.0);
std::string speciesA;
std::string speciesB;
size_t num = compModelNode.nChildren();
for (size_t iChild = 0; iChild < num; iChild++) {
for (size_t iChild = 0; iChild < compModelNode.nChildren(); iChild++) {
XML_Node& xmlChild = compModelNode.child(iChild);
std::string nodeName = lowercase(xmlChild.name());
if (nodeName != "interaction") {
throw CanteraError("TransportFactory::getLiquidInteractionsTransportData",
"expected <interaction> element and got <" + nodeName + ">");
}
speciesA = xmlChild.attrib("speciesA");
speciesB = xmlChild.attrib("speciesB");
string speciesA = xmlChild.attrib("speciesA");
string speciesB = xmlChild.attrib("speciesB");
size_t iSpecies = m_thermo->speciesIndex(speciesA);
if (iSpecies == npos) {
throw CanteraError("TransportFactory::getLiquidInteractionsTransportData",
@ -417,8 +410,7 @@ doublereal LTI_Log_MoleFracs::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
}
}
}
value = exp(value);
return value;
return exp(value);
}
void LTI_Pairwise_Interaction::setParameters(LiquidTransportParams& trParam)
@ -436,22 +428,12 @@ void LTI_Pairwise_Interaction::setParameters(LiquidTransportParams& trParam)
doublereal LTI_Pairwise_Interaction::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
{
size_t nsp = m_thermo->nSpecies();
vector_fp molefracs(nsp);
m_thermo->getMoleFractions(&molefracs[0]);
doublereal value = 0;
throw LTPmodelError("Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense.");
return value;
}
doublereal LTI_Pairwise_Interaction::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
{
size_t nsp = m_thermo->nSpecies();
vector_fp molefracs(nsp);
m_thermo->getMoleFractions(&molefracs[0]);
doublereal value = 0;
throw LTPmodelError("Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense.");
return value;
}
void LTI_Pairwise_Interaction::getMatrixTransProp(DenseMatrix& mat, doublereal* speciesValues)
@ -478,18 +460,17 @@ void LTI_Pairwise_Interaction::getMatrixTransProp(DenseMatrix& mat, doublereal*
void LTI_StefanMaxwell_PPN::setParameters(LiquidTransportParams& trParam)
{
size_t nsp = m_thermo->nSpecies();
size_t nsp2 = nsp*nsp;
m_ionCondMix = 0;
m_ionCondMixModel = trParam.ionConductivity;
m_ionCondSpecies.resize(nsp,0);
m_mobRatMix.resize(nsp,nsp,0.0);
m_mobRatMixModel.resize(nsp2);
m_mobRatSpecies.resize(nsp2);
m_mobRatMixModel.resize(nsp*nsp);
m_mobRatSpecies.resize(nsp*nsp);
m_selfDiffMix.resize(nsp,0.0);
m_selfDiffMixModel.resize(nsp);
m_selfDiffSpecies.resize(nsp);
for (size_t k = 0; k < nsp2; k++) {
for (size_t k = 0; k < nsp*nsp; k++) {
m_mobRatMixModel[k] = trParam.mobilityRatio[k];
m_mobRatSpecies[k].resize(nsp,0);
}
@ -501,7 +482,7 @@ void LTI_StefanMaxwell_PPN::setParameters(LiquidTransportParams& trParam)
for (size_t k = 0; k < nsp; k++) {
LiquidTransportData& ltd = trParam.LTData[k];
m_ionCondSpecies[k] = ltd.ionConductivity;
for (size_t j = 0; j < nsp2; j++) {
for (size_t j = 0; j < nsp*nsp; j++) {
m_mobRatSpecies[j][k] = ltd.mobilityRatio[j];
}
for (size_t j = 0; j < nsp; j++) {
@ -512,22 +493,12 @@ void LTI_StefanMaxwell_PPN::setParameters(LiquidTransportParams& trParam)
doublereal LTI_StefanMaxwell_PPN::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
{
size_t nsp = m_thermo->nSpecies();
vector_fp molefracs(nsp);
m_thermo->getMoleFractions(&molefracs[0]);
doublereal value = 0;
throw LTPmodelError("Calling LTI_StefanMaxwell_PPN::getMixTransProp does not make sense.");
return value;
}
doublereal LTI_StefanMaxwell_PPN::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
{
size_t nsp = m_thermo->nSpecies();
vector_fp molefracs(nsp);
m_thermo->getMoleFractions(&molefracs[0]);
doublereal value = 0;
throw LTPmodelError("Calling LTI_StefanMaxwell_PPN::getMixTransProp does not make sense.");
return value;
}
void LTI_StefanMaxwell_PPN::getMatrixTransProp(DenseMatrix& mat, doublereal* speciesValues)
@ -588,16 +559,13 @@ void LTI_StefanMaxwell_PPN::getMatrixTransProp(DenseMatrix& mat, doublereal* spe
double vM = viS[anion[0]];
double zP = charges[cation[0]];
double zM = charges[anion[0]];
doublereal xA, xB, eps;
doublereal inv_vP_vM_MutualDiff;
vector_fp dlnActCoeffdlnN_diag;
dlnActCoeffdlnN_diag.resize(neut_molefracs.size(),0.0);
vector_fp dlnActCoeffdlnN_diag(neut_molefracs.size(),0.0);
marg_thermo->getdlnActCoeffdlnN_diag(&dlnActCoeffdlnN_diag[0]);
xA = neut_molefracs[neutMolIndex[cation[0]]];
xB = neut_molefracs[neutMolIndex[cation[1]]];
eps = (1-m_mobRatMix(cation[1],cation[0]))/(xA+xB*m_mobRatMix(cation[1],cation[0]));
inv_vP_vM_MutualDiff = (xA*(1-xB+dlnActCoeffdlnN_diag[neutMolIndex[cation[1]]])/m_selfDiffMix[cation[1]]+xB*(1-xA+dlnActCoeffdlnN_diag[neutMolIndex[cation[0]]])/m_selfDiffMix[cation[0]]);
double xA = neut_molefracs[neutMolIndex[cation[0]]];
double xB = neut_molefracs[neutMolIndex[cation[1]]];
double eps = (1-m_mobRatMix(cation[1],cation[0]))/(xA+xB*m_mobRatMix(cation[1],cation[0]));
double inv_vP_vM_MutualDiff = (xA*(1-xB+dlnActCoeffdlnN_diag[neutMolIndex[cation[1]]])/m_selfDiffMix[cation[1]]+xB*(1-xA+dlnActCoeffdlnN_diag[neutMolIndex[cation[0]]])/m_selfDiffMix[cation[0]]);
mat.resize(nsp, nsp, 0.0);
mat(cation[0],cation[1]) = mat(cation[1],cation[0]) = (1+vM/vP)*(1+eps*xB)*(1-eps*xA)*inv_vP_vM_MutualDiff-zP*zP*Faraday*Faraday/GasConstant/temp/m_ionCondMix/vol;
@ -607,22 +575,12 @@ void LTI_StefanMaxwell_PPN::getMatrixTransProp(DenseMatrix& mat, doublereal* spe
doublereal LTI_StokesEinstein::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight)
{
size_t nsp = m_thermo->nSpecies();
vector_fp molefracs(nsp);
m_thermo->getMoleFractions(&molefracs[0]);
doublereal value = 0;
throw LTPmodelError("Calling LTI_StokesEinstein::getMixTransProp does not make sense.");
return value;
}
doublereal LTI_StokesEinstein::getMixTransProp(std::vector<LTPspecies*> LTPptrs)
{
size_t nsp = m_thermo->nSpecies();
vector_fp molefracs(nsp);
m_thermo->getMoleFractions(&molefracs[0]);
doublereal value = 0;
throw LTPmodelError("Calling LTI_StokesEinstein::getMixTransProp does not make sense.");
return value;
}
void LTI_StokesEinstein::setParameters(LiquidTransportParams& trParam)

21
src/transport/LiquidTransport.cpp
View file

@ -563,14 +563,8 @@ void LiquidTransport::set_Grad_X(const doublereal* grad_X)
doublereal LiquidTransport::getElectricConduct()
{
vector_fp gradT(m_nDim,0.0);
vector_fp gradX(m_nDim * m_nsp);
vector_fp gradV(m_nDim);
for (size_t i = 0; i < m_nDim; i++) {
for (size_t k = 0; k < m_nsp; k++) {
gradX[ i*m_nDim + k] = 0.0;
}
gradV[i] = 1.0;
}
vector_fp gradX(m_nDim * m_nsp, 0.0);
vector_fp gradV(m_nDim, 1.0);
set_Grad_T(&gradT[0]);
set_Grad_X(&gradX[0]);
@ -887,9 +881,8 @@ void LiquidTransport::updateHydrodynamicRadius_T()
void LiquidTransport::update_Grad_lnAC()
{
doublereal grad_T;
for (size_t k = 0; k < m_nDim; k++) {
grad_T = m_Grad_T[k];
double grad_T = m_Grad_T[k];
size_t start = m_nsp*k;
m_thermo->getdlnActCoeffds(grad_T, &m_Grad_X[start], &m_Grad_lnAC[start]);
for (size_t i = 0; i < m_nsp; i++) {
@ -900,12 +893,10 @@ void LiquidTransport::update_Grad_lnAC()
}
}
}
return;
}
void LiquidTransport::stefan_maxwell_solve()
{
doublereal tmp;
m_B.resize(m_nsp, m_nDim, 0.0);
m_A.resize(m_nsp, m_nsp, 0.0);
@ -1002,7 +993,7 @@ void LiquidTransport::stefan_maxwell_solve()
m_A(i,i) = 0.0;
for (size_t j = 0; j < m_nsp; j++) {
if (j != i) {
tmp = m_molefracs_tran[j] * m_bdiff(i,j);
double tmp = m_molefracs_tran[j] * m_bdiff(i,j);
m_A(i,i) -= tmp;
m_A(i,j) = tmp;
}
@ -1044,7 +1035,7 @@ void LiquidTransport::stefan_maxwell_solve()
m_A(i,i) = 0.0;
for (size_t j = 0; j < m_nsp; j++) {
if (j != i) {
tmp = m_molefracs_tran[j] * m_bdiff(i,j);
double tmp = m_molefracs_tran[j] * m_bdiff(i,j);
m_A(i,i) -= tmp;
m_A(i,j) = tmp;
}
@ -1084,7 +1075,7 @@ void LiquidTransport::stefan_maxwell_solve()
m_A(i,i) = 0.0;
for (size_t j = 0; j < m_nsp; j++) {
if (j != i) {
tmp = m_molefracs_tran[j] * m_bdiff(i,j);
double tmp = m_molefracs_tran[j] * m_bdiff(i,j);
m_A(i,i) -= tmp;
m_A(i,j) = tmp;
}

8
src/transport/LiquidTransportData.cpp
View file

@ -80,14 +80,10 @@ LiquidTransportData::~LiquidTransportData()
delete ionConductivity;
for (size_t k = 0; k < mobilityRatio.size(); k++) {
if (mobilityRatio[k]) {
delete mobilityRatio[k];
}
delete mobilityRatio[k];
}
for (size_t k = 0; k < selfDiffusion.size(); k++) {
if (selfDiffusion[k]) {
delete selfDiffusion[k];
}
delete selfDiffusion[k];
}
delete thermalCond;

59
src/transport/MMCollisionInt.cpp
View file

@ -20,13 +20,12 @@ double MMCollisionInt::delta[8] = {0.0, 0.25, 0.50, 0.75, 1.0,
doublereal quadInterp(doublereal x0, doublereal* x, doublereal* y)
{
doublereal dx21, dx32, dx31, dy32, dy21, a;
dx21 = x[1] - x[0];
dx32 = x[2] - x[1];
dx31 = dx21 + dx32;
dy32 = y[2] - y[1];
dy21 = y[1] - y[0];
a = (dx21*dy32 - dy21*dx32)/(dx21*dx31*dx32);
double dx21 = x[1] - x[0];
double dx32 = x[2] - x[1];
double dx31 = dx21 + dx32;
double dy32 = y[2] - y[1];
double dy21 = y[1] - y[0];
double a = (dx21*dy32 - dy21*dx32)/(dx21*dx31*dx32);
return a*(x0 - x[0])*(x0 - x[1]) + (dy21/dx21)*(x0 - x[1]) + y[1];
}
@ -237,7 +236,7 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level)
writelogf("T*_max = %g\n", tstar[m_nmax + 1]);
}
m_logTemp.resize(37);
doublereal rmserr, e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0;
doublereal e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0;
if (m_loglevel > 0) {
writelog("Collision integral fits at each tabulated T* vs. delta*.\n"
@ -253,7 +252,7 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level)
m_logTemp[i] = log(tstar[i+1]);
vector_fp c(DeltaDegree+1);
rmserr = fitDelta(0, i, DeltaDegree, c.data());
double rmserr = fitDelta(0, i, DeltaDegree, c.data());
if (log_level > 3) {
writelogf("\ndelta* fit at T* = %.6g\n", tstar[i+1]);
writelog("omega22 = [" + vec2str(c) + "]\n");
@ -319,12 +318,13 @@ doublereal MMCollisionInt::fitDelta(int table, int ntstar, int degree, doublerea
doublereal MMCollisionInt::omega22(double ts, double deltastar)
{
int i;
for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
for (i = 0; i < 37; i++) {
if (ts < tstar22[i]) {
break;
}
int i1, i2;
i1 = std::max(i - 1, 0);
i2 = i1+3;
}
int i1 = std::max(i - 1, 0);
int i2 = i1+3;
if (i2 > 36) {
i2 = 36;
i1 = i2 - 3;
@ -346,9 +346,8 @@ doublereal MMCollisionInt::astar(double ts, double deltastar)
for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
break;
}
int i1, i2;
i1 = std::max(i - 1, 0);
i2 = i1+3;
int i1 = std::max(i - 1, 0);
int i2 = i1+3;
if (i2 > 36) {
i2 = 36;
i1 = i2 - 3;
@ -370,9 +369,8 @@ doublereal MMCollisionInt::bstar(double ts, double deltastar)
for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
break;
}
int i1, i2;
i1 = std::max(i - 1, 0);
i2 = i1+3;
int i1 = std::max(i - 1, 0);
int i2 = i1+3;
if (i2 > 36) {
i2 = 36;
i1 = i2 - 3;
@ -394,9 +392,8 @@ doublereal MMCollisionInt::cstar(double ts, double deltastar)
for (i = 0; i < 37; i++) if (ts < tstar22[i]) {
break;
}
int i1, i2;
i1 = std::max(i - 1,0);
i2 = i1+3;
int i1 = std::max(i - 1,0);
int i2 = i1+3;
if (i2 > 36) {
i2 = 36;
i1 = i2 - 3;
@ -415,12 +412,11 @@ doublereal MMCollisionInt::cstar(double ts, double deltastar)
void MMCollisionInt::fit_omega22(int degree, doublereal deltastar,
doublereal* o22)
{
int i, n = m_nmax - m_nmin + 1;
int n = m_nmax - m_nmin + 1;
vector_fp values(n);
doublereal rmserr;
vector_fp w(n);
doublereal* logT = &m_logTemp[m_nmin];
for (i = 0; i < n; i++) {
for (int i = 0; i < n; i++) {
if (deltastar == 0.0) {
values[i] = omega22_table[8*(i + m_nmin)];
} else {
@ -428,7 +424,7 @@ void MMCollisionInt::fit_omega22(int degree, doublereal deltastar,
}
}
w[0]= -1.0;
rmserr = polyfit(n, degree, logT, values.data(), w.data(), o22);
double rmserr = polyfit(n, degree, logT, values.data(), w.data(), o22);
if (m_loglevel > 0 && rmserr > 0.01) {
writelogf("Warning: RMS error = %12.6g in omega_22 fit"
"with delta* = %12.6g\n", rmserr, deltastar);
@ -438,12 +434,11 @@ void MMCollisionInt::fit_omega22(int degree, doublereal deltastar,
void MMCollisionInt::fit(int degree, doublereal deltastar,
doublereal* a, doublereal* b, doublereal* c)
{
int i, n = m_nmax - m_nmin + 1;
int n = m_nmax - m_nmin + 1;
vector_fp values(n);
doublereal rmserr;
vector_fp w(n);
doublereal* logT = &m_logTemp[m_nmin];
for (i = 0; i < n; i++) {
for (int i = 0; i < n; i++) {
if (deltastar == 0.0) {
values[i] = astar_table[8*(i + m_nmin + 1)];
} else {
@ -451,9 +446,9 @@ void MMCollisionInt::fit(int degree, doublereal deltastar,
}
}
w[0]= -1.0;
rmserr = polyfit(n, degree, logT, values.data(), w.data(), a);
double rmserr = polyfit(n, degree, logT, values.data(), w.data(), a);
for (i = 0; i < n; i++) {
for (int i = 0; i < n; i++) {
if (deltastar == 0.0) {
values[i] = bstar_table[8*(i + m_nmin + 1)];
} else {
@ -463,7 +458,7 @@ void MMCollisionInt::fit(int degree, doublereal deltastar,
w[0]= -1.0;
rmserr = polyfit(n, degree, logT, values.data(), w.data(), b);
for (i = 0; i < n; i++) {
for (int i = 0; i < n; i++) {
if (deltastar == 0.0) {
values[i] = cstar_table[8*(i + m_nmin + 1)];
} else {

119
src/transport/MultiTransport.cpp
View file

@ -219,23 +219,20 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_
}
// copy grad_X to fluxes
const doublereal* gx;
for (size_t n = 0; n < ndim; n++) {
gx = grad_X + ldx*n;
const double* gx = grad_X + ldx*n;
copy(gx, gx + m_nsp, fluxes + ldf*n);
fluxes[jmax + n*ldf] = 0.0;
}
// solve the equations
solve(m_aa, fluxes, ndim, ldf);
size_t offset;
doublereal pp = pressure_ig();
// multiply diffusion velocities by rho * V to create mass fluxes, and
// restore the gradx elements that were modified
for (size_t n = 0; n < ndim; n++) {
offset = n*ldf;
size_t offset = n*ldf;
for (size_t i = 0; i < m_nsp; i++) {
fluxes[i + offset] *= rho * y[i] / pp;
}
@ -244,7 +241,7 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_
// thermal diffusion
if (addThermalDiffusion) {
for (size_t n = 0; n < ndim; n++) {
offset = n*ldf;
size_t offset = n*ldf;
doublereal grad_logt = grad_T[n]/m_temp;
for (size_t i = 0; i < m_nsp; i++) {
fluxes[i + offset] -= m_spwork[i]*grad_logt;
@ -259,7 +256,7 @@ void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* s
double* x1 = m_spwork1.data();
double* x2 = m_spwork2.data();
double* x3 = m_spwork3.data();
size_t n, nsp = m_thermo->nSpecies();
size_t nsp = m_thermo->nSpecies();
m_thermo->restoreState(nsp+2, state1);
double p1 = m_thermo->pressure();
double t1 = state1[0];
@ -273,7 +270,7 @@ void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* s
double p = 0.5*(p1 + p2);
double t = 0.5*(state1[0] + state2[0]);
for (n = 0; n < nsp; n++) {
for (size_t n = 0; n < nsp; n++) {
x3[n] = 0.5*(x1[n] + x2[n]);
}
m_thermo->setState_TPX(t, p, x3);
@ -380,10 +377,9 @@ void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
doublereal prefactor = 16.0 * m_temp
* m_thermo->meanMolecularWeight()/(25.0 * p);
doublereal c;
for (size_t i = 0; i < m_nsp; i++) {
for (size_t j = 0; j < m_nsp; j++) {
c = prefactor/m_mw[j];
double c = prefactor/m_mw[j];
d[ld*j + i] = c*m_molefracs[i]*
(m_Lmatrix(i,j) - m_Lmatrix(i,i));
}
@ -430,12 +426,10 @@ void MultiTransport::updateThermal_T()
updateDiff_T();
// evaluate polynomial fits for A*, B*, C*
doublereal z;
int ipoly;
for (size_t i = 0; i < m_nsp; i++) {
for (size_t j = i; j < m_nsp; j++) {
z = m_logt - m_log_eps_k(i,j);
ipoly = m_poly[i][j];
double z = m_logt - m_log_eps_k(i,j);
int ipoly = m_poly[i][j];
if (m_mode == CK_Mode) {
m_om22(i,j) = poly6(z, m_omega22_poly[ipoly].data());
m_astar(i,j) = poly6(z, m_astar_poly[ipoly].data());
@ -456,18 +450,15 @@ void MultiTransport::updateThermal_T()
m_abc_ok = true;
// evaluate the temperature-dependent rotational relaxation rate
doublereal tr, sqtr;
for (size_t k = 0; k < m_nsp; k++) {
tr = m_eps[k]/ m_kbt;
sqtr = m_sqrt_eps_k[k] / m_sqrt_t;
double tr = m_eps[k]/ m_kbt;
double sqtr = m_sqrt_eps_k[k] / m_sqrt_t;
m_rotrelax[k] = std::max(1.0,m_zrot[k]) * m_frot_298[k]/Frot(tr, sqtr);
}
doublereal d;
doublereal c = 1.2*GasConstant*m_temp;
for (size_t k = 0; k < m_nsp; k++) {
d = c * m_visc[k] * m_astar(k,k)/m_mw[k];
m_bdiff(k,k) = d;
m_bdiff(k,k) = c * m_visc[k] * m_astar(k,k)/m_mw[k];
}
// Calculate the internal heat capacities by subtracting off the translational contributions
@ -520,11 +511,10 @@ void MultiTransport::eval_L0000(const doublereal* const x)
void MultiTransport::eval_L0010(const doublereal* const x)
{
doublereal prefactor = 1.6*m_temp;
doublereal sum, wj, xj;
for (size_t j = 0; j < m_nsp; j++) {
xj = x[j];
wj = m_mw[j];
sum = 0.0;
double xj = x[j];
double wj = m_mw[j];
double sum = 0.0;
for (size_t i = 0; i < m_nsp; i++) {
m_Lmatrix(i,j + m_nsp) = - prefactor * x[i] * xj * m_mw[i] *
(1.2 * m_cstar(j,i) - 1.0) /
@ -551,26 +541,22 @@ void MultiTransport::eval_L1010(const doublereal* x)
{
const doublereal fiveover3pi = 5.0/(3.0*Pi);
doublereal prefactor = (16.0*m_temp)/25.0;
doublereal constant1, wjsq, constant2, constant3, constant4,
fourmj, threemjsq, sum, sumwij;;
doublereal term1, term2;
for (size_t j = 0; j < m_nsp; j++) {
// get constant terms that depend on just species "j"
constant1 = prefactor*x[j];
wjsq = m_mw[j]*m_mw[j];
constant2 = 13.75*wjsq;
constant3 = m_crot[j]/m_rotrelax[j];
constant4 = 7.5*wjsq;
fourmj = 4.0*m_mw[j];
threemjsq = 3.0*m_mw[j]*m_mw[j];
sum = 0.0;
double constant1 = prefactor*x[j];
double wjsq = m_mw[j]*m_mw[j];
double constant2 = 13.75*wjsq;
double constant3 = m_crot[j]/m_rotrelax[j];
double constant4 = 7.5*wjsq;
double fourmj = 4.0*m_mw[j];
double threemjsq = 3.0*m_mw[j]*m_mw[j];
double sum = 0.0;
for (size_t i = 0; i < m_nsp; i++) {
sumwij = m_mw[i] + m_mw[j];
term1 = m_bdiff(i,j) * sumwij*sumwij;
term2 = fourmj*m_astar(i,j)*(1.0 + fiveover3pi*
(constant3 +
(m_crot[i]/m_rotrelax[i]))); // see Eq. (12.125)
double sumwij = m_mw[i] + m_mw[j];
double term1 = m_bdiff(i,j) * sumwij*sumwij;
double term2 = fourmj*m_astar(i,j)*(1.0 + fiveover3pi*
(constant3 + (m_crot[i]/m_rotrelax[i]))); // see Eq. (12.125)
m_Lmatrix(i+m_nsp,j+m_nsp) = constant1*x[i]*m_mw[i] /(m_mw[j]*term1) *
(constant2 - threemjsq*m_bstar(i,j)
@ -588,25 +574,21 @@ void MultiTransport::eval_L1010(const doublereal* x)
void MultiTransport::eval_L1001(const doublereal* x)
{
doublereal prefactor = 32.00*m_temp/(5.00*Pi);
doublereal constant, sum;
size_t n2 = 2*m_nsp;
int npoly = 0;
for (size_t j = 0; j < m_nsp; j++) {
// collect terms that depend only on "j"
if (hasInternalModes(j)) {
constant = prefactor*m_mw[j]*x[j]*m_crot[j]/(m_cinternal[j]*m_rotrelax[j]);
sum = 0.0;
double constant = prefactor*m_mw[j]*x[j]*m_crot[j]/(m_cinternal[j]*m_rotrelax[j]);
double sum = 0.0;
for (size_t i = 0; i < m_nsp; i++) {
// see Eq. (12.127)
m_Lmatrix(i+m_nsp,j+n2) = constant * m_astar(j,i) * x[i] /
m_Lmatrix(i+m_nsp,j+2*m_nsp) = constant * m_astar(j,i) * x[i] /
((m_mw[j] + m_mw[i]) * m_bdiff(j,i));
sum += m_Lmatrix(i+m_nsp,j+n2);
sum += m_Lmatrix(i+m_nsp,j+2*m_nsp);
}
npoly++;
m_Lmatrix(j+m_nsp,j+n2) += sum;
m_Lmatrix(j+m_nsp,j+2*m_nsp) += sum;
} else {
for (size_t i = 0; i < m_nsp; i++) {
m_Lmatrix(i+m_nsp,j+n2) = 0.0;
m_Lmatrix(i+m_nsp,j+2*m_nsp) = 0.0;
}
}
}
@ -614,64 +596,57 @@ void MultiTransport::eval_L1001(const doublereal* x)
void MultiTransport::eval_L0001()
{
size_t n2 = 2*m_nsp;
for (size_t j = 0; j < m_nsp; j++) {
for (size_t i = 0; i < m_nsp; i++) {
m_Lmatrix(i,j+n2) = 0.0;
m_Lmatrix(i,j+2*m_nsp) = 0.0;
}
}
}
void MultiTransport::eval_L0100()
{
size_t n2 = 2*m_nsp;
for (size_t j = 0; j < m_nsp; j++) {
for (size_t i = 0; i < m_nsp; i++) {
m_Lmatrix(i+n2,j) = 0.0; // see Eq. (12.123)
m_Lmatrix(i+2*m_nsp,j) = 0.0; // see Eq. (12.123)
}
}
}
void MultiTransport::eval_L0110()
{
size_t n2 = 2*m_nsp;
for (size_t j = 0; j < m_nsp; j++) {
for (size_t i = 0; i < m_nsp; i++) {
m_Lmatrix(i+n2,j+m_nsp) = m_Lmatrix(j+m_nsp,i+n2); // see Eq. (12.123)
m_Lmatrix(i+2*m_nsp,j+m_nsp) = m_Lmatrix(j+m_nsp,i+2*m_nsp); // see Eq. (12.123)
}
}
}
void MultiTransport::eval_L0101(const doublereal* x)
{
const doublereal fivepi = 5.00*Pi;
const doublereal eightoverpi = 8.0 / Pi;
doublereal prefactor = 4.00*m_temp;
size_t n2 = 2*m_nsp;
doublereal constant1, constant2, diff_int, sum;
for (size_t i = 0; i < m_nsp; i++) {
if (hasInternalModes(i)) {
// collect terms that depend only on "i"
constant1 = prefactor*x[i]/m_cinternal[i];
constant2 = 12.00*m_mw[i]*m_crot[i] /
(fivepi*m_cinternal[i]*m_rotrelax[i]);
sum = 0.0;
double constant1 = 4*m_temp*x[i]/m_cinternal[i];
double constant2 = 12*m_mw[i]*m_crot[i] /
(5*Pi*m_cinternal[i]*m_rotrelax[i]);
double sum = 0.0;
for (size_t k = 0; k < m_nsp; k++) {
// see Eq. (12.131)
diff_int = m_bdiff(i,k);
m_Lmatrix(k+n2,i+n2) = 0.0;
double diff_int = m_bdiff(i,k);
m_Lmatrix(k+2*m_nsp,i+2*m_nsp) = 0.0;
sum += x[k]/diff_int;
if (k != i) sum += x[k]*m_astar(i,k)*constant2 /
(m_mw[k]*diff_int);
if (k != i) {
sum += x[k]*m_astar(i,k)*constant2 / (m_mw[k]*diff_int);
}
}
// see Eq. (12.130)
m_Lmatrix(i+n2,i+n2) =
- eightoverpi*m_mw[i]*x[i]*x[i]*m_crot[i] /
m_Lmatrix(i+2*m_nsp,i+2*m_nsp) =
- 8/Pi*m_mw[i]*x[i]*x[i]*m_crot[i] /
(m_cinternal[i]*m_cinternal[i]*GasConstant*m_visc[i]*m_rotrelax[i])
- constant1*sum;
} else {
for (size_t k = 0; k < m_nsp; k++) {
m_Lmatrix(i+n2,i+n2) = 1.0;
m_Lmatrix(i+2*m_nsp,i+2*m_nsp) = 1.0;
}
}
}

9
src/transport/SimpleTransport.cpp
View file

@ -176,9 +176,7 @@ bool SimpleTransport::initLiquid(LiquidTransportParams& tr)
m_viscSpecies.resize(m_nsp);
m_coeffVisc_Ns.clear();
m_coeffVisc_Ns.resize(m_nsp);
std::string spName = m_thermo->speciesName(0);
for (size_t k = 0; k < m_nsp; k++) {
spName = m_thermo->speciesName(k);
LiquidTransportData& ltd = tr.LTData[k];
m_coeffVisc_Ns[k] = ltd.viscosity;
ltd.viscosity = 0;
@ -189,7 +187,6 @@ bool SimpleTransport::initLiquid(LiquidTransportParams& tr)
m_coeffLambda_Ns.clear();
m_coeffLambda_Ns.resize(m_nsp);
for (size_t k = 0; k < m_nsp; k++) {
spName = m_thermo->speciesName(k);
LiquidTransportData& ltd = tr.LTData[k];
m_coeffLambda_Ns[k] = ltd.thermalCond;
ltd.thermalCond = 0;
@ -201,7 +198,7 @@ bool SimpleTransport::initLiquid(LiquidTransportParams& tr)
m_coeffDiff_Ns.clear();
m_coeffDiff_Ns.resize(m_nsp);
for (size_t k = 0; k < m_nsp; k++) {
spName = m_thermo->speciesName(k);
string spName = m_thermo->speciesName(k);
LiquidTransportData& ltd = tr.LTData[k];
m_coeffDiff_Ns[k] = ltd.speciesDiffusivity;
ltd.speciesDiffusivity = 0;
@ -280,7 +277,6 @@ void SimpleTransport::getSpeciesViscosities(doublereal* const visc)
void SimpleTransport::getBinaryDiffCoeffs(size_t ld, doublereal* d)
{
double bdiff;
update_T();
// if necessary, evaluate the species diffusion coefficients
@ -291,8 +287,7 @@ void SimpleTransport::getBinaryDiffCoeffs(size_t ld, doublereal* d)
for (size_t i = 0; i < m_nsp; i++) {
for (size_t j = 0; j < m_nsp; j++) {
bdiff = 0.5 * (m_diffSpecies[i] + m_diffSpecies[j]);
d[i*m_nsp+j] = bdiff;
d[i*m_nsp+j] = 0.5 * (m_diffSpecies[i] + m_diffSpecies[j]);
}
}
}

3
src/transport/SolidTransport.cpp
View file

@ -160,8 +160,7 @@ void SolidTransport::getMobilities(doublereal* const mobil)
void SolidTransport::getMixDiffCoeffs(doublereal* const d)
{
size_t nsp = m_thermo->nSpecies();
for (size_t k = 0; k < nsp; k++) {
for (size_t k = 0; k < m_thermo->nSpecies(); k++) {
d[k] = 0.0;
}
}

104
src/transport/TransportFactory.cpp
View file

@ -100,26 +100,19 @@ std::string TransportFactory::modelName(int model)
LTPspecies* TransportFactory::newLTP(const XML_Node& trNode, const std::string& name,
TransportPropertyType tp_ind, thermo_t* thermo)
{
LTPspecies* ltps = 0;
std::string model = lowercase(trNode["model"]);
switch (m_LTRmodelMap[model]) {
case LTP_TD_CONSTANT:
ltps = new LTPspecies_Const(trNode, name, tp_ind, thermo);
break;
return new LTPspecies_Const(trNode, name, tp_ind, thermo);
case LTP_TD_ARRHENIUS:
ltps = new LTPspecies_Arrhenius(trNode, name, tp_ind, thermo);
break;
return new LTPspecies_Arrhenius(trNode, name, tp_ind, thermo);
case LTP_TD_POLY:
ltps = new LTPspecies_Poly(trNode, name, tp_ind, thermo);
break;
return new LTPspecies_Poly(trNode, name, tp_ind, thermo);
case LTP_TD_EXPT:
ltps = new LTPspecies_ExpT(trNode, name, tp_ind, thermo);
break;
return new LTPspecies_ExpT(trNode, name, tp_ind, thermo);
default:
throw CanteraError("TransportFactory::newLTP","unknown transport model: " + model);
ltps = new LTPspecies(&trNode, name, tp_ind, thermo);
}
return ltps;
}
LiquidTranInteraction* TransportFactory::newLTI(const XML_Node& trNode,
@ -127,54 +120,52 @@ LiquidTranInteraction* TransportFactory::newLTI(const XML_Node& trNode,
LiquidTransportParams& trParam)
{
LiquidTranInteraction* lti = 0;
thermo_t* thermo = trParam.thermo;
std::string model = trNode["model"];
switch (m_LTImodelMap[model]) {
switch (m_LTImodelMap[trNode["model"]]) {
case LTI_MODEL_SOLVENT:
lti = new LTI_Solvent(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
break;
case LTI_MODEL_MOLEFRACS:
lti = new LTI_MoleFracs(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
break;
case LTI_MODEL_MASSFRACS:
lti = new LTI_MassFracs(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
break;
case LTI_MODEL_LOG_MOLEFRACS:
lti = new LTI_Log_MoleFracs(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
break;
case LTI_MODEL_PAIRWISE_INTERACTION:
lti = new LTI_Pairwise_Interaction(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
lti->setParameters(trParam);
break;
case LTI_MODEL_STEFANMAXWELL_PPN:
lti = new LTI_StefanMaxwell_PPN(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
lti->setParameters(trParam);
break;
case LTI_MODEL_STOKES_EINSTEIN:
lti = new LTI_StokesEinstein(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
lti->setParameters(trParam);
break;
case LTI_MODEL_MOLEFRACS_EXPT:
lti = new LTI_MoleFracs_ExpT(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
break;
case LTI_MODEL_NOTSET:
case LTI_MODEL_NONE:
case LTI_MODEL_MULTIPLE:
lti = new LiquidTranInteraction(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
break;
default:
// @TODO make sure we can throw an error here with existing datasets and tests before changing code
lti = new LiquidTranInteraction(tp_ind);
lti->init(trNode, thermo);
lti->init(trNode, trParam.thermo);
}
return lti;
}
@ -269,9 +260,7 @@ void TransportFactory::setupLiquidTransport(thermo_t* thermo, int log_level,
trParam.log_level = log_level;
// Get the molecular weights and load them into trParam
trParam.mw.resize(nsp);
copy(trParam.thermo->molecularWeights().begin(),
trParam.thermo->molecularWeights().end(), trParam.mw.begin());
trParam.mw = trParam.thermo->molecularWeights();
// Resize all other vectors in trParam
trParam.LTData.resize(nsp);
@ -282,7 +271,7 @@ void TransportFactory::setupLiquidTransport(thermo_t* thermo, int log_level,
trParam.diff_Dij.resize(nsp,nsp);
trParam.radius_Aij.resize(nsp,nsp);
XML_Node root, log;
XML_Node log;
// Note that getLiquidSpeciesTransportData just populates the pure species transport data.
getLiquidSpeciesTransportData(species_database, log, trParam.thermo->speciesNames(), trParam);
@ -302,21 +291,17 @@ void TransportFactory::setupSolidTransport(thermo_t* thermo, int log_level,
// constant mixture attributes
trParam.thermo = thermo;
trParam.nsp_ = trParam.thermo->nSpecies();
size_t nsp = trParam.nsp_;
trParam.tmin = thermo->minTemp();
trParam.tmax = thermo->maxTemp();
trParam.log_level = log_level;
// Get the molecular weights and load them into trParam
trParam.mw.resize(nsp);
copy(trParam.thermo->molecularWeights().begin(),
trParam.thermo->molecularWeights().end(), trParam.mw.begin());
XML_Node root, log;
trParam.mw = trParam.thermo->molecularWeights();
// getSolidTransportData() populates the phase transport models like
// electronic conductivity thermal conductivity, interstitial diffusion
if (phase_database->hasChild("transport")) {
XML_Node log;
XML_Node& transportNode = phase_database->child("transport");
getSolidTransportData(transportNode, log, thermo->name(), trParam);
}
@ -347,8 +332,6 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML
const std::vector<std::string> &names,
LiquidTransportParams& trParam)
{
std::string name;
// Create a map of species names versus liquid transport data parameters
std::map<std::string, LiquidTransportData> datatable;
@ -363,8 +346,7 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML
// listed in 'names'.
for (size_t i = 0; i < nsp; i++) {
const XML_Node& sp = *xspecies[i];
name = sp["name"];
vector_fp vCoeff;
string name = sp["name"];
// Species with no 'transport' child are skipped. However, if that
// species is in the list, it will throw an exception below.
@ -377,7 +359,6 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML
data.speciesName = name;
data.mobilityRatio.resize(nsp*nsp,0);
data.selfDiffusion.resize(nsp,0);
ThermoPhase* temp_thermo = trParam.thermo;
size_t num = trNode.nChildren();
for (size_t iChild = 0; iChild < num; iChild++) {
XML_Node& xmlChild = trNode.child(iChild);
@ -385,10 +366,10 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML
switch (m_tranPropMap[nodeName]) {
case TP_VISCOSITY:
data.viscosity = newLTP(xmlChild, name, m_tranPropMap[nodeName], temp_thermo);
data.viscosity = newLTP(xmlChild, name, m_tranPropMap[nodeName], trParam.thermo);
break;
case TP_IONCONDUCTIVITY:
data.ionConductivity = newLTP(xmlChild, name, m_tranPropMap[nodeName], temp_thermo);
data.ionConductivity = newLTP(xmlChild, name, m_tranPropMap[nodeName], trParam.thermo);
break;
case TP_MOBILITYRATIO: {
for (size_t iSpec = 0; iSpec< nBinInt; iSpec++) {
@ -397,8 +378,8 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML
size_t loc = specName.find(":");
std::string firstSpec = specName.substr(0,loc);
std::string secondSpec = specName.substr(loc+1);
size_t index = temp_thermo->speciesIndex(firstSpec)+nsp*temp_thermo->speciesIndex(secondSpec);
data.mobilityRatio[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], temp_thermo);
size_t index = trParam.thermo->speciesIndex(firstSpec)+nsp*trParam.thermo->speciesIndex(secondSpec);
data.mobilityRatio[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], trParam.thermo);
};
};
break;
@ -406,8 +387,8 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML
for (size_t iSpec = 0; iSpec< nsp; iSpec++) {
XML_Node& propSpecNode = xmlChild.child(iSpec);
std::string specName = propSpecNode.name();
size_t index = temp_thermo->speciesIndex(specName);
data.selfDiffusion[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], temp_thermo);
size_t index = trParam.thermo->speciesIndex(specName);
data.selfDiffusion[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], trParam.thermo);
};
};
break;
@ -415,25 +396,25 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML
data.thermalCond = newLTP(xmlChild,
name,
m_tranPropMap[nodeName],
temp_thermo);
trParam.thermo);
break;
case TP_DIFFUSIVITY:
data.speciesDiffusivity = newLTP(xmlChild,
name,
m_tranPropMap[nodeName],
temp_thermo);
trParam.thermo);
break;
case TP_HYDRORADIUS:
data.hydroRadius = newLTP(xmlChild,
name,
m_tranPropMap[nodeName],
temp_thermo);
trParam.thermo);
break;
case TP_ELECTCOND:
data.electCond = newLTP(xmlChild,
name,
m_tranPropMap[nodeName],
temp_thermo);
trParam.thermo);
break;
default:
throw CanteraError("getLiquidSpeciesTransportData","unknown transport property: " + nodeName);
@ -451,12 +432,11 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vector<const XML
if (it == datatable.end()) {
throw TransportDBError(0,"No transport data found for species " + names[i]);
}
LiquidTransportData& trdat = it->second;
// Now, transfer these objects into LTData in the correct phase index
// order by calling the default copy constructor for
// LiquidTransportData.
trParam.LTData.push_back(trdat);
trParam.LTData.push_back(it->second);
}
}
@ -475,14 +455,12 @@ void TransportFactory::getLiquidInteractionsTransportData(const XML_Node& transp
try {
size_t nsp = trParam.nsp_;
size_t nBinInt = nsp*(nsp-1)/2;
size_t num = transportNode.nChildren();
for (size_t iChild = 0; iChild < num; iChild++) {
for (size_t iChild = 0; iChild < transportNode.nChildren(); iChild++) {
//tranTypeNode is a type of transport property like viscosity
XML_Node& tranTypeNode = transportNode.child(iChild);
std::string nodeName = tranTypeNode.name();
trParam.mobilityRatio.resize(nsp*nsp,0);
trParam.selfDiffusion.resize(nsp,0);
ThermoPhase* temp_thermo = trParam.thermo;
if (tranTypeNode.name() == "compositionDependence") {
std::string modelName = tranTypeNode.attrib("model");
@ -514,7 +492,7 @@ void TransportFactory::getLiquidInteractionsTransportData(const XML_Node& transp
size_t loc = specName.find(":");
string firstSpec = specName.substr(0,loc);
string secondSpec = specName.substr(loc+1);
size_t index = temp_thermo->speciesIndex(firstSpec)+nsp*temp_thermo->speciesIndex(secondSpec);
size_t index = trParam.thermo->speciesIndex(firstSpec)+nsp*trParam.thermo->speciesIndex(secondSpec);
trParam.mobilityRatio[index] = newLTI(propSpecNode,
m_tranPropMap[nodeName],
trParam);
@ -525,7 +503,7 @@ void TransportFactory::getLiquidInteractionsTransportData(const XML_Node& transp
for (size_t iSpec = 0; iSpec< nsp; iSpec++) {
XML_Node& propSpecNode = compDepNode.child(iSpec);
string specName = propSpecNode.name();
size_t index = temp_thermo->speciesIndex(specName);
size_t index = trParam.thermo->speciesIndex(specName);
trParam.selfDiffusion[index] = newLTI(propSpecNode,
m_tranPropMap[nodeName],
trParam);
@ -591,39 +569,37 @@ void TransportFactory::getSolidTransportData(const XML_Node& transportNode,
const std::string phaseName,
SolidTransportData& trParam)
{
size_t num = transportNode.nChildren();
for (size_t iChild = 0; iChild < num; iChild++) {
for (size_t iChild = 0; iChild < transportNode.nChildren(); iChild++) {
//tranTypeNode is a type of transport property like viscosity
XML_Node& tranTypeNode = transportNode.child(iChild);
std::string nodeName = tranTypeNode.name();
ThermoPhase* temp_thermo = trParam.thermo;
//tranTypeNode contains the interaction model
switch (m_tranPropMap[nodeName]) {
case TP_IONCONDUCTIVITY:
trParam.ionConductivity = newLTP(tranTypeNode, phaseName,
m_tranPropMap[nodeName],
temp_thermo);
trParam.thermo);
break;
case TP_THERMALCOND:
trParam.thermalConductivity = newLTP(tranTypeNode, phaseName,
m_tranPropMap[nodeName],
temp_thermo);
trParam.thermo);
break;
case TP_DEFECTDIFF:
trParam.defectDiffusivity = newLTP(tranTypeNode, phaseName,
m_tranPropMap[nodeName],
temp_thermo);
trParam.thermo);
break;
case TP_DEFECTCONC:
trParam.defectActivity = newLTP(tranTypeNode, phaseName,
m_tranPropMap[nodeName],
temp_thermo);
trParam.thermo);
break;
case TP_ELECTCOND:
trParam.electConductivity = newLTP(tranTypeNode, phaseName,
m_tranPropMap[nodeName],
temp_thermo);
trParam.thermo);
break;
default:
throw CanteraError("getSolidTransportData","unknown transport property: " + nodeName);