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

Worked on the mixture diffusion coefficients. Not done yet

Browse source
This commit is contained in:
Harry Moffat 2009-02-19 16:19:50 +00:00
parent f05b4f51e9
commit 1ecc55e64e
2 changed files with 123 additions and 52 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

116
Cantera/src/transport/LiquidTransport.cpp
View file

@ -23,7 +23,7 @@ using namespace std;
* Mole fractions below MIN_X will be set to MIN_X when computing
* transport properties.
*/
#define MIN_X 1.e-20
#define MIN_X 1.e-14
namespace Cantera {
@ -443,6 +443,35 @@ namespace Cantera {
}
}
void LiquidTransport::getSpeciesDiffusiveMassFluxes(doublereal* const fluxes) {
int n, k;
update_temp();
update_conc();
getMixDiffCoeffs(DATA_PTR(m_spwork));
const array_fp& mw = m_thermo->molecularWeights();
const doublereal* const y = m_thermo->massFractions();
const doublereal rhon = m_thermo->molarDensity();
// Unroll wrt ndim
vector_fp sum(m_nDim,0.0);
for (n = 0; n < m_nDim; n++) {
for (k = 0; k < m_nsp; k++) {
fluxes[n*m_nsp + k] = -rhon * mw[k] * m_spwork[k] * m_Grad_X[n*m_nsp + k];
sum[n] += fluxes[n*m_nsp + k];
}
}
// add correction flux to enforce sum to zero
for (n = 0; n < m_nDim; n++) {
for (k = 0; k < m_nsp; k++) {
fluxes[n*m_nsp + k] -= y[k]*sum[n];
}
}
}
/**
* Mixture-averaged diffusion coefficients [m^2/s].
*
@ -463,24 +492,28 @@ namespace Cantera {
int k, j;
doublereal mmw = m_thermo->meanMolecularWeight();
doublereal sumxw = 0.0, sum2;
doublereal p = m_press;
doublereal sumxw_tran = 0.0;
doublereal sum2;
if (m_nsp == 1) {
d[0] = m_bdiff(0,0) / p;
d[0] = m_bdiff(0,0);
} else {
for (k = 0; k < m_nsp; k++) sumxw += m_molefracs[k] * m_mw[k];
for (k = 0; k < m_nsp; k++) {
sumxw_tran += m_molefracs_tran[k] * m_mw[k];
}
for (k = 0; k < m_nsp; k++) {
sum2 = 0.0;
for (j = 0; j < m_nsp; j++) {
if (j != k) {
sum2 += m_molefracs[j] / m_bdiff(j,k);
sum2 += m_molefracs_tran[j] / m_bdiff(j,k);
}
}
if (sum2 <= 0.0) {
d[k] = m_bdiff(k,k) / p;
} else {
d[k] = (sumxw - m_molefracs[k] * m_mw[k])/(p * mmw * sum2);
}
// Because we use m_molefracs_tran, sum2 must be positive definate
// if (sum2 <= 0.0) {
// d[k] = m_bdiff(k,k);
// } else {
d[k] = (sumxw_tran - m_molefracs_tran[k] * m_mw[k])/(mmw * sum2);
// }
}
}
}
@ -569,12 +602,15 @@ namespace Cantera {
if (iStateNew != m_iStateMF) {
qReturn = false;
m_thermo->getMoleFractions(DATA_PTR(m_molefracs));
m_thermo->getMoleFractions(DATA_PTR(m_concentrations));
m_thermo->getConcentrations(DATA_PTR(m_concentrations));
double ctot = 0.0;
for (int k = 0; k < m_nsp; k++) {
m_molefracs[k] = fmaxx(MIN_X, m_molefracs[k]);
m_molefracs[k] = fmaxx(0.0, m_molefracs[k]);
m_molefracs_tran[k] = fmaxx(MIN_X, m_molefracs[k]);
ctot += m_concentrations[k];
}
dens_ = m_thermo->density();
meanMolecularWeight_ = m_thermo->meanMolecularWeight();
double ctotmin = 0.0;
for (int k = 0; k < m_nsp; k++) {
m_concentrations[k]= fmaxx(ctotmin, m_concentrations[k]);
@ -597,6 +633,12 @@ namespace Cantera {
// We formulate the directional derivative
/*
* We only calculate the change in ac due to composition.
* The pressure and the temperature are taken care of in
* other parts of the expression.
*
*/
void LiquidTransport::update_Grad_lnAC() {
int k;
@ -632,25 +674,12 @@ namespace Cantera {
}
for (k = 0; k < m_nsp; k++) {
m_Grad_lnAC[m_nsp * a + k] = sum * (lnActCoeffMolarDelta_[k] - log(actCoeffMolar_[k])) / mag;
m_Grad_lnAC[m_nsp * a + k] =
sum * (lnActCoeffMolarDelta_[k] - log(actCoeffMolar_[k])) / mag;
}
}
m_thermo->setMoleFractions(DATA_PTR(m_molefracs));
double Tbase = m_thermo->temperature();
double T_new = Tbase - 1.0E-6;
m_thermo->setTemperature(T_new);
m_thermo->getActivityCoefficients(DATA_PTR(lnActCoeffMolarDelta_));
double *dlnActCoeffdT = &Xdelta_[0];
for (k = 0; k < m_nsp; k++) {
dlnActCoeffdT[k] == (lnActCoeffMolarDelta_[k] - log(actCoeffMolar_[k]))/(-1.0E-6);
}
for (int a = 0; a < m_nDim; a++) {
for (k = 0; k < m_nsp; k++) {
m_Grad_lnAC[m_nsp * a + k] += dlnActCoeffdT[k] * m_Grad_T[a];
}
}
m_thermo->setTemperature(Tbase);
}
/*************************************************************************
@ -780,11 +809,7 @@ namespace Cantera {
//! grab a local copy of the molecular weights
const vector_fp& M = m_thermo->molecularWeights();
//! get the mean molecular weight of the mixture
//double M_mix = m_thermo->meanMolecularWeight();
/*
* Update the concentrations in the mixture.
*/
@ -792,15 +817,24 @@ namespace Cantera {
double T = m_thermo->temperature();
m_thermo->getEntropy_R(DATA_PTR(entropy_R_specSS_));
m_thermo->getStandardVolumes(DATA_PTR(volume_specSS_));
m_thermo->getStandardVolumes(DATA_PTR(volume_specPM_));
m_thermo->getActivityCoefficients(DATA_PTR(actCoeffMolar_));
/*
* Calculate the electrochemical potential gradient. This is the
* driving force for relative diffusional transport.
*
* Here we calculate c_i * grad (mu_i), p. 297 Newman
* Here we calculate
*
* c_i * (grad (mu_i) + S_i grad T - M_i / dens * grad P
*
* This is Eqn. 13-1 p. 318 Newman. The original equation is from
* Hershfeld, Curtis, and Bird.
*
* S_i is the partial molar entropy of species i. This term will cancel
* out a lot of the grad T terms in grad (mu_i), therefore simplifying
* the expression.
*
* Ok I think there may be many ways to do this. One way is to do it via basis
* functions, at the nodes, as a function of the variables in the problem.
@ -808,15 +842,12 @@ namespace Cantera {
* For calculation of molality based thermo systems, we current get
* the molar based values. This may change.
*
*
*/
for (i = 0; i < m_nsp; i++) {
for (a = 0; a < VIM; a++) {
m_ck_Grad_mu[a*m_nsp + i] =
m_chargeSpecies[i] * m_concentrations[i] * Faraday * m_Grad_V[a]
+ m_concentrations[i] * GasConstant * entropy_R_specSS_[i] * m_Grad_T[a]
+ m_concentrations[i] * volume_specSS_[i] * m_Grad_P[a]
+ m_concentrations[i] * GasConstant * m_Grad_T[a] * log(actCoeffMolar_[i] * m_molefracs[i])
+ m_concentrations[i] * (volume_specPM_[i] - M[i]/dens_) * m_Grad_P[a]
+ m_concentrations[i] * GasConstant * T * m_Grad_lnAC[a*m_nsp+i] / actCoeffMolar_[i]
+ concTot_ * GasConstant * T * m_Grad_X[a*m_nsp+i];
}
@ -829,7 +860,8 @@ namespace Cantera {
double lnmnaught = log(mnaught);
for (i = 1; i < m_nsp; i++) {
for (a = 0; a < VIM; a++) {
m_ck_Grad_mu[a*m_nsp + i] -= m_concentrations[i] * GasConstant * m_Grad_T[a] * lnmnaught;
m_ck_Grad_mu[a*m_nsp + i] -=
m_concentrations[i] * GasConstant * m_Grad_T[a] * lnmnaught;
}
}
}

59
Cantera/src/transport/LiquidTransport.h
View file

@ -33,6 +33,10 @@ namespace Cantera {
const int LVISC_WILKES = 1;
const int LVISC_MIXTUREAVG = 2;
const int LDIFF_MIXDIFF_UNCORRECTED = 0;
const int LDIFF_MIXDIFF_FLUXCORRECTED = 1;
const int LDIFF_MULTICOMP_STEFANMAXWELL = 2;
class TransportParams;
@ -302,6 +306,20 @@ namespace Cantera {
int ldx, const doublereal* grad_X,
int ldf, doublereal* fluxes);
//! Return the species diffusive mass fluxes
/*!
*
*
*
* @param ndim The number of spatial dimensions (1, 2, or 3).
* @param grad_T The temperature gradient (ignored in this model).
* @param ldx Leading dimension of the grad_X array.
* The diffusive mass flux of species \e k is computed from
*
*
*/
virtual void getSpeciesDiffusiveMassFluxes(doublereal* const fluxes);
/**
* @param ndim The number of spatial dimensions (1, 2, or 3).
* @param grad_T The temperature gradient (ignored in this model).
@ -445,9 +463,13 @@ namespace Cantera {
//! Array of Binary Diffusivities
/*!
* This has a size equal to nsp x nsp
* It is a symmetric matrix.
* D_ii is undefined.
* Depends on the temperature. We have set the pressure dependence
* to zero for this liquid phase constituitve model
*
* This has a size equal to nsp x nsp
* It is a symmetric matrix.
* D_ii is the self diffusion coefficient. D_ii is not
* needed except for when there is one species in the mixture.
*
* units m2/sec
*/
@ -458,8 +480,8 @@ namespace Cantera {
* Viscosity of the species
* Length = number of species
*
* Depends on the temperature and perhaps pressure, but
* not the species concentrations
* Depends on the temperature. We have set the pressure dependence
* to zero for this liquid phase constituitve model
*
* controlling update boolean -> m_visc_temp_ok
*/
@ -496,11 +518,23 @@ namespace Cantera {
//! Local copy of the mole fractions of the species in the phase
/*!
* The mole fractions here are assumed to be bounded by 0.0, and 1.0
* and they are assumed to add up to one.
* Update info?
* length = m_nsp
*/
vector_fp m_molefracs;
//! Mole fraction vector
/*!
* The mole fractions here are assumed to be bounded by MIN_X and 1.0
* and they may not be assumed to add up to one.
*
* Update info?
* length = m_nsp
*/
vector_fp m_molefracs_tran;
vector_fp Xdelta_;
//! Local copy of the concentrations of the species in the phase
@ -512,6 +546,8 @@ namespace Cantera {
//! Local copy of the total concentration
doublereal concTot_;
doublereal meanMolecularWeight_;
doublereal dens_;
//! Local copy of the charge of each species
/*!
@ -519,9 +555,8 @@ namespace Cantera {
*/
vector_fp m_chargeSpecies;
vector_fp entropy_R_specSS_;
vector_fp volume_specSS_;
vector_fp volume_specPM_;
vector_fp actCoeffMolar_;
@ -601,9 +636,13 @@ namespace Cantera {
doublereal m_press;
//! Solution of the flux system
Array2D m_flux;
/*!
* This is the mass flux of species k
* in units of kg m-3 s-1.
*/
Array2D m_flux;
//! saved value of the mixture thermal conductivity
//! Saved value of the mixture thermal conductivity
doublereal m_lambda;
//! Saved value of the mixture viscosity