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

Merging the changes that had been occurring on the change branch

Browse source 
liquidTransportDevelop into the main trunk. This seemed like a good time
to do it. And, there were some fundamental bugs that were fixed just on
the change branch.
This commit is contained in:
Harry Moffat 2009-10-12 19:26:15 +00:00
commit 2ee219c50b
51 changed files with 2777 additions and 472 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

21
Cantera/src/base/xml.cpp
View file

@ -482,6 +482,27 @@ namespace Cantera {
}
}
void XML_Node::clear() {
int n = static_cast<int>(m_children.size());
for (int i = 0; i < n; i++) {
if (m_children[i]) {
if (m_children[i]->parent() == this) {
delete m_children[i];
m_children[i] = 0;
}
}
}
m_value.clear();
m_childindex.clear();
m_attribs.clear();
m_children.clear();
m_nchildren = 0;
m_iscomment = false;
m_linenum = 0;
}
// Add a child node to the current node containing a comment
/*
* Child node will have the name, "comment".

9
Cantera/src/base/xml.h
View file

@ -380,6 +380,13 @@ namespace Cantera {
*/
std::string attrib(const std::string & attr) const;
//! Clear the current node and everything under it
/*!
* The value, attributes and children are all zeroed. The name and the
* parent information is kept.
*/
void clear();
private:
//! Returns a changeable value of the attributes map for the current node
/*!
@ -644,7 +651,7 @@ namespace Cantera {
//! into the destination XML_Node tree, doing a complete copy
//! as we go.
/*!
* Note this is a const function becuase the current XML_Node and
* Note this is a const function because the current XML_Node and
* its children isn't altered by this operation.
*
* @param node_dest This is the XML node to receive the information

0
Cantera/src/kinetics/Enhanced3BConc.h Executable file → Normal file
View file

0
Cantera/src/kinetics/FalloffFactory.cpp Executable file → Normal file
View file

0
Cantera/src/kinetics/FalloffFactory.h Executable file → Normal file
View file

0
Cantera/src/kinetics/FalloffMgr.h Executable file → Normal file
View file

0
Cantera/src/kinetics/GRI_30_Kinetics.cpp Executable file → Normal file
View file

0
Cantera/src/kinetics/GRI_30_Kinetics.h Executable file → Normal file
View file

0
Cantera/src/kinetics/GasKinetics.cpp Executable file → Normal file
View file

0
Cantera/src/kinetics/GasKinetics.h Executable file → Normal file
View file

0
Cantera/src/kinetics/GasKineticsWriter.cpp Executable file → Normal file
View file

0
Cantera/src/kinetics/GasKineticsWriter.h Executable file → Normal file
View file

0
Cantera/src/kinetics/Group.cpp Executable file → Normal file
View file

0
Cantera/src/kinetics/Group.h Executable file → Normal file
View file

0
Cantera/src/kinetics/ImplicitChem.cpp Executable file → Normal file
View file

0
Cantera/src/kinetics/ImplicitChem.h Executable file → Normal file
View file

0
Cantera/src/kinetics/ImplicitSurfChem.cpp Executable file → Normal file
View file

0
Cantera/src/kinetics/ImplicitSurfChem.h Executable file → Normal file
View file

0
Cantera/src/kinetics/Kinetics.h Executable file → Normal file
View file

0
Cantera/src/kinetics/RateCoeffMgr.h Executable file → Normal file
View file

0
Cantera/src/kinetics/ReactionData.h Executable file → Normal file
View file

0
Cantera/src/kinetics/ReactionPath.cpp Executable file → Normal file
View file

0
Cantera/src/kinetics/ReactionPath.h Executable file → Normal file
View file

0
Cantera/src/kinetics/RxnRates.h Executable file → Normal file
View file

0
Cantera/src/kinetics/StoichManager.h Executable file → Normal file
View file

0
Cantera/src/kinetics/ThirdBodyMgr.h Executable file → Normal file
View file

0
Cantera/src/kinetics/reaction_defs.h Executable file → Normal file
View file

17
Cantera/src/numerics/Integrator.h Executable file → Normal file
View file

@ -64,7 +64,14 @@ namespace Cantera {
public:
virtual ~Integrator() {}
//! Default Constructor
Integrator()
{
}
//! Destructor
virtual ~Integrator() {
}
/** Set or reset the number of equations. */
//virtual void resize(int n)=0;
@ -76,10 +83,12 @@ namespace Cantera {
* @param abstol array of N absolute tolerance values
*/
virtual void setTolerances(doublereal reltol, int n,
doublereal* abstol) { warn("setTolerances"); }
doublereal* abstol) {
warn("setTolerances");
}
/**
* Set error tolerances.
//! Set error tolerances.
/*!
* @param reltol scalar relative tolerance
* @param abstol scalar absolute tolerance
*/

11
Cantera/src/thermo/ThermoFactory.cpp
View file

@ -418,6 +418,17 @@ namespace Cantera {
", is not a phase element.");
}
/*
* In this section of code, we get the reference to the
* phase xml tree within the ThermoPhase object. Then,
* we clear it and fill it with the current information that
* we are about to use to construct the object. We will then
* be able to resurrect the information later by calling xml().
*/
XML_Node &phaseNode_XML = th->xml();
phaseNode_XML.clear();
phase.copy(&phaseNode_XML);
// set the id attribute of the phase to the 'id' attribute
// in the XML tree.
th->setID(phase.id());

62
Cantera/src/transport/AqueousTransport.cpp
View file

@ -19,6 +19,7 @@
#include "utilities.h"
#include "TransportParams.h"
#include "LiquidTransportParams.h"
#include "TransportFactory.h"
#include "ctlapack.h"
@ -72,7 +73,7 @@ namespace Cantera {
/*
* This is where we dimension everything.
*/
bool AqueousTransport::init(TransportParams& tr) {
bool AqueousTransport::initLiquid( LiquidTransportParams& tr ) {
// constant substance attributes
m_thermo = tr.thermo;
@ -86,15 +87,11 @@ namespace Cantera {
m_thermo->molecularWeights().end(), m_mw.begin());
// copy polynomials and parameters into local storage
m_poly = tr.poly;
m_visccoeffs = tr.visccoeffs;
m_condcoeffs = tr.condcoeffs;
m_diffcoeffs = tr.diffcoeffs;
m_mode = tr.mode;
m_diam = tr.diam;
m_eps = tr.eps;
m_alpha = tr.alpha;
m_mode = tr.mode_;
m_phi.resize(m_nsp, m_nsp, 0.0);
@ -184,6 +181,7 @@ namespace Cantera {
/******************* binary diffusion coefficients **************/
//================================================================================================
void AqueousTransport::getBinaryDiffCoeffs(const int ld, doublereal* const d) {
int i,j;
@ -200,32 +198,52 @@ namespace Cantera {
d[ld*j + i] = rp * m_bdiff(i,j);
}
}
//================================================================================================
// Get the electrical Mobilities (m^2/V/s).
/*
* This function returns the mobilities. In some formulations
* this is equal to the normal mobility multiplied by faraday's constant.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^e_k = \frac{F D_k}{R T}
* \f]
*
* @param mobil_e Returns the mobilities of
* the species in array \c mobil_e. The array must be
* dimensioned at least as large as the number of species.
*/
void AqueousTransport::getMobilities(doublereal* const mobil) {
// this needs to be checked out.
int k;
getMixDiffCoeffs(DATA_PTR(m_spwork));
doublereal c1 = ElectronCharge / (Boltzmann * m_temp);
for (k = 0; k < m_nsp; k++) {
mobil[k] = c1 * m_spwork[k] * m_thermo->charge(k);
mobil[k] = c1 * m_spwork[k];
}
}
//================================================================================================
void AqueousTransport::getFluidMobilities(doublereal* const mobil) {
getMixDiffCoeffs(DATA_PTR(m_spwork));
doublereal c1 = 1.0 / (GasConstant * m_temp);
for (int k = 0; k < m_nsp; k++) {
mobil[k] = c1 * m_spwork[k];
}
}
//================================================================================================
void AqueousTransport::set_Grad_V(const doublereal* const grad_V) {
for (int a = 0; a < m_nDim; a++) {
m_Grad_V[a] = grad_V[a];
}
}
//================================================================================================
void AqueousTransport::set_Grad_T(const doublereal* const grad_T) {
for (int a = 0; a < m_nDim; a++) {
m_Grad_T[a] = grad_T[a];
}
}
//================================================================================================
void AqueousTransport::set_Grad_X(const doublereal* const grad_X) {
int itop = m_nDim * m_nsp;
for (int i = 0; i < itop; i++) {
@ -585,17 +603,15 @@ namespace Cantera {
* This function returns a Transport data object for a given species.
*
*/
struct GasTransportData AqueousTransport::
getGasTransportData(int kSpecies)
struct LiquidTransportData AqueousTransport::
getLiquidTransportData(int kSpecies)
{
struct GasTransportData td;
struct LiquidTransportData td;
td.speciesName = m_thermo->speciesName(kSpecies);
td.wellDepth = m_eps[kSpecies] / Boltzmann;
td.diameter = m_diam(kSpecies, kSpecies) * 1.0E10;
td.polarizability = m_alpha[kSpecies] * 1.0E30;
/* NEEDS WORK
td.hydroradius = ???;
*/
return td;
}

51
Cantera/src/transport/AqueousTransport.h
View file

@ -11,13 +11,15 @@
#ifndef CT_AQUEOUSTRAN_H
#define CT_AQYEOUSTRAN_H
#define CT_AQUEOUSTRAN_H
using namespace std;
// Cantera includes
#include "TransportBase.h"
#include "DenseMatrix.h"
#include "TransportParams.h"
#include "LiquidTransportParams.h"
#include <vector>
@ -29,7 +31,7 @@ using namespace std;
namespace Cantera {
class TransportParams;
class LiquidTransportParams;
//! Class AqueousTransport implements mixture-averaged transport
@ -133,7 +135,7 @@ namespace Cantera {
virtual ~AqueousTransport() {}
//! Return the model id for this transport parameterization
virtual int model() { return cAqueousTransport; }
virtual int model() const { return cAqueousTransport; }
//! overloaded base class methods
@ -192,12 +194,43 @@ namespace Cantera {
*/
virtual void getMixDiffCoeffs(doublereal* const d);
//! Get the Mobilities
//! Get the Electrical mobilities (m^2/V/s).
/*!
* @param mobil
* This function returns the electrical mobilities. In some formulations
* this is equal to the normal mobility multiplied by faraday's constant.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^e_k = \frac{F D_k}{R T}
* \f]
*
* @param mobil_e Returns the mobilities of
* the species in array \c mobil_e. The array must be
* dimensioned at least as large as the number of species.
*/
virtual void getMobilities(doublereal* const mobil);
virtual void getMobilities(doublereal* const mobil_e);
//! Get the fluid mobilities (s kmol/kg).
/*!
* This function returns the fluid mobilities. Usually, you have
* to multiply Faraday's constant into the resulting expression
* to general a species flux expression.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^f_k = \frac{D_k}{R T}
* \f]
*
* @param mobil_f Returns the mobilities of
* the species in array \c mobil_f. The array must be
* dimensioned at least as large as the number of species.
*/
virtual void getFluidMobilities(doublereal* const mobil_f);
//! Specify the value of the gradient of the voltage
/*!
@ -275,7 +308,7 @@ namespace Cantera {
* @param tr Transport parameters for all of the species
* in the phase.
*/
virtual bool init(TransportParams& tr);
virtual bool initLiquid( LiquidTransportParams& tr );
friend class TransportFactory;
@ -286,7 +319,7 @@ namespace Cantera {
*
* @param k Species number to obtain the properties about.
*/
struct GasTransportData getGasTransportData(int k);
struct LiquidTransportData getLiquidTransportData(int k);
//! Solve the stefan_maxell equations for the diffusive fluxes.

2
Cantera/src/transport/DustyGasTransport.h
View file

@ -40,7 +40,7 @@ namespace Cantera {
//---------------------------------------------------------
// overloaded base class methods
virtual int model() { return cDustyGasTransport; }
virtual int model() const { return cDustyGasTransport; }
virtual void setParameters(const int type, const int k, const doublereal* const p);

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

@ -39,10 +39,6 @@ namespace Cantera {
m_iStateMF(-1),
m_temp(-1.0),
m_logt(0.0),
m_sqrt_t(-1.0),
m_t14(-1.0),
m_t32(-1.0),
m_sqrt_kbt(-1.0),
m_press(-1.0),
m_lambda(-1.0),
m_viscmix(-1.0),
@ -67,10 +63,6 @@ namespace Cantera {
m_iStateMF(-1),
m_temp(-1.0),
m_logt(0.0),
m_sqrt_t(-1.0),
m_t14(-1.0),
m_t32(-1.0),
m_sqrt_kbt(-1.0),
m_press(-1.0),
m_lambda(-1.0),
m_viscmix(-1.0),
@ -100,36 +92,35 @@ namespace Cantera {
m_tmin = right.m_tmin;
m_tmax = right.m_tmax;
m_mw = right.m_mw;
viscCoeffsVector_ = right.viscCoeffsVector_;
m_condcoeffs = right.m_condcoeffs;
m_visc_A = right.m_visc_A;
m_visc_logA = right.m_visc_logA;
m_visc_n = right.m_visc_n;
m_visc_Tact = right.m_visc_Tact;
m_visc_Eij = right.m_visc_Eij;
m_visc_Sij = right.m_visc_Sij;
m_thermCond_A = right.m_thermCond_A;
m_thermCond_n = right.m_thermCond_n;
m_thermCond_Tact = right.m_thermCond_Tact;
m_hydrodynamic_radius = right.m_hydrodynamic_radius;
m_diffcoeffs = right.m_diffcoeffs;
m_Grad_X = right.m_Grad_X;
m_Grad_T = right.m_Grad_T;
m_Grad_V = right.m_Grad_V;
m_ck_Grad_mu = right.m_ck_Grad_mu;
m_bdiff = right.m_bdiff;
viscSpecies_ = right.viscSpecies_;
m_sqvisc = right.m_sqvisc;
m_cond = right.m_cond;
m_polytempvec = right.m_polytempvec;
m_viscSpecies = right.m_viscSpecies;
m_logViscSpecies = right.m_logViscSpecies;
m_condSpecies = right.m_condSpecies;
m_iStateMF = -1;
m_molefracs = right.m_molefracs;
m_concentrations = right.m_concentrations;
m_chargeSpecies = right.m_chargeSpecies;
m_DiffCoeff_StefMax = right.m_DiffCoeff_StefMax;
viscosityModel_ = right.viscosityModel_;
m_phi = right.m_phi;
m_wratjk = right.m_wratjk;
m_wratkj1 = right.m_wratkj1;
m_B = right.m_B;
m_A = right.m_A;
m_eps = right.m_eps;
m_temp = right.m_temp;
m_logt = right.m_logt;
m_sqrt_t = right.m_sqrt_t;
m_t14 = right.m_t14;
m_t32 = right.m_t32;
m_sqrt_kbt = right.m_sqrt_kbt;
m_press = right.m_press;
m_flux = right.m_flux;
m_lambda = right.m_lambda;
@ -173,32 +164,34 @@ namespace Cantera {
copy(m_thermo->molecularWeights().begin(),
m_thermo->molecularWeights().end(), m_mw.begin());
// copy polynomials and parameters into local storage
viscCoeffsVector_ = tr.viscCoeffsVector_;
m_condcoeffs = tr.condcoeffs;
// copy parameters into local storage
m_visc_A = tr.visc_A ;
m_visc_n = tr.visc_n ;
m_visc_Tact = tr.visc_Tact ;
//The following two are not yet filled in LiquidTransportParams
m_visc_Eij = tr.visc_Eij ;
m_visc_Sij = tr.visc_Sij ;
//save logarithm of pre-exponential for easier computation
m_visc_logA.resize(m_nsp);
for ( int i = 0; i < m_nsp; i++ )
m_visc_logA[i] = log( m_visc_A[i] );
m_thermCond_A = tr.thermCond_A ;
m_thermCond_n = tr.thermCond_n ;
m_thermCond_Tact = tr.thermCond_Tact ;
m_hydrodynamic_radius = tr.hydroRadius ;
//m_diffcoeffs = tr.diffcoeffs;
m_mode = tr.mode;
m_diam = tr.diam;
m_eps = tr.eps;
m_mode = tr.mode_;
m_phi.resize(m_nsp, m_nsp, 0.0);
m_wratjk.resize(m_nsp, m_nsp, 0.0);
m_wratkj1.resize(m_nsp, m_nsp, 0.0);
int j, k;
for (j = 0; j < m_nsp; j++)
for (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]);
}
m_polytempvec.resize(5);
viscSpecies_.resize(m_nsp);
m_sqvisc.resize(m_nsp);
m_cond.resize(m_nsp);
m_viscSpecies.resize(m_nsp);
m_logViscSpecies.resize(m_nsp);
m_condSpecies.resize(m_nsp);
m_bdiff.resize(m_nsp, m_nsp);
m_molefracs.resize(m_nsp);
@ -249,7 +242,7 @@ namespace Cantera {
if (m_visc_mix_ok) return m_viscmix;
// update viscSpecies_[] and m_phi[] if necessary
// update m_viscSpecies[] if necessary
if (!m_visc_temp_ok) {
updateViscosity_temp();
}
@ -258,16 +251,27 @@ namespace Cantera {
updateViscosities_conc();
}
/* We still need to implement interaction parameters */
/* This constant viscosity model has no input */
if (viscosityModel_ == LVISC_CONSTANT) {
return m_viscmix;
} else if (viscosityModel_ == LVISC_MIXTUREAVG) {
m_viscmix = dot_product(viscSpecies_, m_molefracs);
} else if (viscosityModel_ == LVISC_WILKES) {
multiply(m_phi, DATA_PTR(m_molefracs), DATA_PTR(m_spwork));
m_viscmix = 0.0;
for (int k = 0; k < m_nsp; k++) {
m_viscmix += m_molefracs[k] * viscSpecies_[k]/m_spwork[k];
}
err("constant viscosity not implemented for LiquidTransport.");
//return m_viscmix;
} else if (viscosityModel_ == LVISC_AVG_ENERGIES) {
m_viscmix = exp( dot_product(m_logViscSpecies, m_molefracs) );
} else if (viscosityModel_ == LVISC_INTERACTION) {
// log_visc_mix = sum_i (X_i log_visc_i) + sum_i sum_j X_i X_j G_ij
double interaction = dot_product(m_logViscSpecies, m_molefracs);
for ( int i = 0; i < m_nsp; i++ )
for ( int j = 0; j < i; j++ )
interaction += m_molefracs[i] * m_molefracs[j]
* ( m_visc_Sij(i,j) + m_visc_Eij(i,j) / m_temp );
m_viscmix = exp( interaction );
}
return m_viscmix;
@ -278,7 +282,7 @@ namespace Cantera {
if (!m_visc_temp_ok) {
updateViscosity_temp();
}
copy(viscSpecies_.begin(), viscSpecies_.end(), visc);
copy(m_viscSpecies.begin(), m_viscSpecies.end(), visc);
}
@ -301,41 +305,79 @@ namespace Cantera {
d[ld*j + i] = rp * m_bdiff(i,j);
}
}
//================================================================================================
// Get the electrical Mobilities (m^2/V/s).
/*
* This function returns the mobilities. In some formulations
* this is equal to the normal mobility multiplied by faraday's constant.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^e_k = \frac{F D_k}{R T}
* \f]
*
* @param mobil_e Returns the mobilities of
* the species in array \c mobil_e. The array must be
* dimensioned at least as large as the number of species.
*/
void LiquidTransport::getMobilities(doublereal* const mobil) {
// this needs to be checked out.
int k;
getMixDiffCoeffs(DATA_PTR(m_spwork));
doublereal c1 = ElectronCharge / (Boltzmann * m_temp);
for (k = 0; k < m_nsp; k++) {
mobil[k] = c1 * m_spwork[k] * m_thermo->charge(k);
mobil[k] = c1 * m_spwork[k];
}
}
//================================================================================================
//! Get the fluid mobilities (s kmol/kg).
/*!
* This function returns the fluid mobilities. Usually, you have
* to multiply Faraday's constant into the resulting expression
* to general a species flux expression.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^f_k = \frac{D_k}{R T}
* \f]
*
*
* @param mobil_f Returns the mobilities of
* the species in array \c mobil. The array must be
* dimensioned at least as large as the number of species.
*/
void LiquidTransport::getFluidMobilities(doublereal* const mobil_f) {
getMixDiffCoeffs(DATA_PTR(m_spwork));
doublereal c1 = 1.0 / (GasConstant * m_temp);
for (int k = 0; k < m_nsp; k++) {
mobil_f[k] = c1 * m_spwork[k];
}
}
//================================================================================================
void LiquidTransport::set_Grad_V(const doublereal* const grad_V) {
for (int a = 0; a < m_nDim; a++) {
m_Grad_V[a] = grad_V[a];
}
}
//================================================================================================
void LiquidTransport::set_Grad_T(const doublereal* const grad_T) {
for (int a = 0; a < m_nDim; a++) {
m_Grad_T[a] = grad_T[a];
}
}
void LiquidTransport::set_Grad_X(const doublereal* const grad_X) {
int itop = m_nDim * m_nsp;
for (int i = 0; i < itop; i++) {
m_Grad_X[i] = grad_X[i];
}
update_Grad_lnAC();
}
//================================================================================================
void LiquidTransport::set_Grad_X(const doublereal* const grad_X) {
int itop = m_nDim * m_nsp;
for (int i = 0; i < itop; i++) {
m_Grad_X[i] = grad_X[i];
}
update_Grad_lnAC();
}
//================================================================================================
/****************** thermal conductivity **********************/
/*
@ -356,8 +398,8 @@ namespace Cantera {
if (!m_cond_mix_ok) {
doublereal sum1 = 0.0, sum2 = 0.0;
for (int k = 0; k < m_nsp; k++) {
sum1 += m_molefracs[k] * m_cond[k];
sum2 += m_molefracs[k] / m_cond[k];
sum1 += m_molefracs[k] * m_condSpecies[k];
sum2 += m_molefracs[k] / m_condSpecies[k];
}
m_lambda = 0.5*(sum1 + 1.0/sum2);
m_cond_mix_ok = true;
@ -439,35 +481,6 @@ 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].
*
@ -542,18 +555,6 @@ namespace Cantera {
m_temp = t;
m_logt = log(m_temp);
m_kbt = Boltzmann * m_temp;
m_sqrt_t = sqrt(m_temp);
m_t14 = sqrt(m_sqrt_t);
m_t32 = m_temp * m_sqrt_t;
m_sqrt_kbt = sqrt(Boltzmann*m_temp);
// compute powers of log(T)
// -> may move this
m_polytempvec[0] = 1.0;
m_polytempvec[1] = m_logt;
m_polytempvec[2] = m_logt*m_logt;
m_polytempvec[3] = m_logt*m_logt*m_logt;
m_polytempvec[4] = m_logt*m_logt*m_logt*m_logt;
// temperature has changed so temp flags are flipped
m_visc_temp_ok = false;
@ -689,18 +690,20 @@ namespace Cantera {
*/
void LiquidTransport::updateCond_temp() {
int k;
/*
if (m_mode == CK_Mode) {
for (k = 0; k < m_nsp; k++) {
m_cond[k] = exp(dot4(m_polytempvec, m_condcoeffs[k]));
m_condSpecies[k] = exp(m_condcoeffs[k]);
}
} else {
for (k = 0; k < m_nsp; k++) {
m_cond[k] = m_sqrt_t * dot5(m_polytempvec, m_condcoeffs[k]);
m_condSpecies[k] = m_sqrt_t * m_condcoeffs[k];
}
}
m_cond_temp_ok = true;
m_cond_mix_ok = false;
*/
}
@ -711,12 +714,12 @@ namespace Cantera {
void LiquidTransport::updateDiff_temp() {
// evaluate binary diffusion coefficients at unit pressure
int i,j;
int ic = 0;
/*
if (m_mode == CK_Mode) {
for (i = 0; i < m_nsp; i++) {
for (j = i; j < m_nsp; j++) {
m_bdiff(i,j) = exp(dot4(m_polytempvec, m_diffcoeffs[ic]));
m_bdiff(i,j) = exp(m_diffcoeffs[ic]);
m_bdiff(j,i) = m_bdiff(i,j);
ic++;
}
@ -725,8 +728,7 @@ namespace Cantera {
else {
for (i = 0; i < m_nsp; i++) {
for (j = i; j < m_nsp; j++) {
m_bdiff(i,j) = m_temp * m_sqrt_t*dot5(m_polytempvec,
m_diffcoeffs[ic]);
m_bdiff(i,j) = m_temp * m_sqrt_t*m_diffcoeffs[ic];
m_bdiff(j,i) = m_bdiff(i,j);
ic++;
}
@ -735,6 +737,7 @@ namespace Cantera {
m_diff_temp_ok = true;
m_diff_mix_ok = false;
*/
}
@ -754,38 +757,16 @@ namespace Cantera {
*/
void LiquidTransport::updateViscosity_temp() {
int k;
doublereal vratiokj, wratiojk, factor1;
if (m_mode == CK_Mode) {
for (k = 0; k < m_nsp; k++) {
viscSpecies_[k] = exp(dot4(m_polytempvec, viscCoeffsVector_[k]));
m_sqvisc[k] = sqrt(viscSpecies_[k]);
}
for (k = 0; k < m_nsp; k++) {
m_logViscSpecies[k] = m_visc_logA[k] + m_visc_n[k] * m_logt
+ m_visc_Tact[k] / m_temp ;
m_viscSpecies[k] = exp( m_logViscSpecies[k] );
}
else {
for (k = 0; k < m_nsp; k++) {
// the polynomial fit is done for sqrt(visc/sqrt(T))
m_sqvisc[k] = m_t14*dot5(m_polytempvec, viscCoeffsVector_[k]);
viscSpecies_[k] = (m_sqvisc[k]*m_sqvisc[k]);
}
}
// see Eq. (9-5.15) of Reid, Prausnitz, and Poling
int j;
for (j = 0; j < m_nsp; j++) {
for (k = j; k < m_nsp; k++) {
vratiokj = viscSpecies_[k]/viscSpecies_[j];
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);
m_phi(k,j) = factor1*factor1 /
(SqrtEight * m_wratkj1(j,k));
m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk);
}
}
//for (k = 0; k < m_nsp; k++) {
//m_viscSpecies[k] = m_visc_A[k] * exp( m_visc_n[k] * m_logt
// + m_visc_Tact[k] / m_temp );
//}
m_visc_temp_ok = true;
m_visc_mix_ok = false;
}
@ -953,4 +934,20 @@ namespace Cantera {
}
}
}
/**
* Throw an exception if this method is invoked.
* This probably indicates something is not yet implemented.
*/
doublereal LiquidTransport::err(std::string msg) const {
throw CanteraError("Liquid Transport Class",
"\n\n\n**** Method "+ msg +" not implemented in model "
+ int2str(model()) + " ****\n"
"(Did you forget to specify a transport model?)\n\n\n");
return 0.0;
}
}

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

@ -25,13 +25,14 @@ using namespace std;
// Cantera includes
#include "TransportBase.h"
#include "DenseMatrix.h"
#include "TransportParams.h"
#include "LiquidTransportParams.h"
namespace Cantera {
const int LVISC_CONSTANT = 0;
const int LVISC_WILKES = 1;
const int LVISC_MIXTUREAVG = 2;
const int LVISC_INTERACTION = 1;
const int LVISC_AVG_ENERGIES = 2;
const int LDIFF_MIXDIFF_UNCORRECTED = 0;
const int LDIFF_MIXDIFF_FLUXCORRECTED = 1;
@ -39,7 +40,7 @@ namespace Cantera {
class TransportParams;
class LiquidTransportParams;
//! Class LiquidTransport implements mixture-averaged transport
@ -135,12 +136,21 @@ namespace Cantera {
class LiquidTransport : public Transport {
public:
//! default constructor
//! Default constructor.
/*!
* This requires call to initLiquid(LiquidTransportParams& tr)
* after filling LiquidTransportParams to complete instantiation.
* The filling of LiquidTransportParams is currently carried out
* in the TransportFactory class, but might be moved at some point.
*
* @param thermo ThermoPhase object holding species information.
* @param ndim Number of spatial dimensions.
*/
LiquidTransport(thermo_t* thermo = 0, int ndim = 1);
//!Copy Constructor for the %LiquidThermo object.
/*!
* @param right ThermoPhase to be copied
* @param right %LiquidTransport to be copied
*/
LiquidTransport(const LiquidTransport &right);
@ -148,8 +158,8 @@ namespace Cantera {
/*!
* This is NOT a virtual function.
*
* @param right Reference to %ThermoPhase object to be copied into the
* current one.
* @param right Reference to %LiquidTransport object to be copied
* into the current one.
*/
LiquidTransport& operator=(const LiquidTransport& right);
@ -169,8 +179,21 @@ namespace Cantera {
//! virtual destructor
virtual ~LiquidTransport() {}
//! Initialize the transport object
/*!
* Here we change all of the internal dimensions to be sufficient.
* We get the object ready to do property evaluations.
*
* @param tr Transport parameters for all of the species
* in the phase.
*/
virtual bool initLiquid(LiquidTransportParams& tr);
friend class TransportFactory;
//! Return the model id for this transport parameterization
virtual int model() {
virtual int model() const {
return cLiquidTransport;
}
@ -178,17 +201,14 @@ namespace Cantera {
//! Returns the viscosity of the solution
/*!
* The viscosity is computed using the Wilke mixture rule.
* The viscosity is computed using mixture averaging plus
* any information on interaction parameters
* \f[
* \mu = \sum_k \frac{\mu_k X_k}{\sum_j \Phi_{k,j} X_j}.
* \mu = \sum_k {\mu_k X_k} {\sum_j \sum_k {G_{j,k} X_k X_j} }.
* \f]
* Here \f$ \mu_k \f$ is the viscosity of pure species \e k,
* and
* \f[
* \Phi_{k,j} = \frac{\left[1
* + \sqrt{\left(\frac{\mu_k}{\mu_j}\sqrt{\frac{M_j}{M_k}}\right)}\right]^2}
* {\sqrt{8}\sqrt{1 + M_k/M_j}}
* \f]
* and \f$ G_{k,j} \f$ is the interaction parameter.
* @see updateViscosity_T();
*
* Controlling update boolean m_viscmix_ok
@ -197,25 +217,11 @@ namespace Cantera {
//! Returns the pure species viscosities
/*!
*
*
* The pure species viscosities are to be given in an Arrhenius
* form in accordance with activated-jump-process dominated transport.
*/
virtual void getSpeciesViscosities(doublereal* const visc);
virtual void getThermalDiffCoeffs(doublereal* const dt);
//! Return the thermal conductivity of the solution
/*!
* The thermal conductivity is computed from the following mixture rule:
* \f[
* \lambda = 0.5 \left( \sum_k X_k \lambda_k
* + \frac{1}{\sum_k X_k/\lambda_k}\right)
* \f]
*
* Controlling update boolean = m_condmix_ok
*/
virtual doublereal thermalConductivity();
//! Returns the binary diffusion coefficients
/*!
* @param ld
@ -231,11 +237,56 @@ namespace Cantera {
virtual void getMixDiffCoeffs(doublereal* const d);
//! Get the Mobilities
virtual void getThermalDiffCoeffs(doublereal* const dt);
//! Return the thermal conductivity of the solution
/*!
* @param mobil
* The thermal conductivity is computed from the following mixture rule:
* \f[
* \lambda = 0.5 \left( \sum_k X_k \lambda_k
* + \frac{1}{\sum_k X_k/\lambda_k}\right)
* \f]
*
* Controlling update boolean = m_condmix_ok
*/
virtual void getMobilities(doublereal* const mobil);
virtual doublereal thermalConductivity();
//! Get the Electrical mobilities (m^2/V/s).
/*!
* This function returns the mobilities. In some formulations
* this is equal to the normal mobility multiplied by faraday's constant.
*
* The mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^e_k = \frac{F D_k}{R T}
* \f]
*
* @param mobil_e Returns the electrical mobilities of
* the species in array \c mobil_e. The array must be
* dimensioned at least as large as the number of species.
*/
virtual void getMobilities(doublereal* const mobil_e);
//! Get the fluid mobilities (s kmol/kg).
/*!
* This function returns the fluid mobilities. Usually, you have
* to multiply Faraday's constant into the resulting expression
* to general a species flux expression.
*
* The mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^f_k = \frac{D_k}{R T}
* \f]
*
* @param mobil_f Returns the fluid mobilities of
* the species in array \c mobil_f. The array must be
* dimensioned at least as large as the number of species.
*/
virtual void getFluidMobilities(doublereal* const mobil_f);
//! Specify the value of the gradient of the voltage
/*!
@ -306,19 +357,6 @@ 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).
@ -331,20 +369,6 @@ namespace Cantera {
virtual void getSpeciesFluxesExt(int ldf, doublereal* fluxes);
//! Initialize the transport object
/*!
* Here we change all of the internal dimensions to be sufficient.
* We get the object ready to do property evaluations.
*
* @param tr Transport parameters for all of the species
* in the phase.
*/
virtual bool initLiquid(LiquidTransportParams& tr);
friend class TransportFactory;
//! Solve the stefan_maxell equations for the diffusive fluxes.
void stefan_maxwell_solve();
@ -368,19 +392,34 @@ namespace Cantera {
*/
vector_fp m_mw;
//! Polynomial coefficients of the viscosity
/*!
* These express the temperature dependendence of the pures
* species viscosities.
*/
std::vector<vector_fp> viscCoeffsVector_;
//! Pure species viscosities in Arrhenius temperature-dependent form.
vector_fp m_visc_A;
vector_fp m_visc_logA; //logarithm of coefficient
vector_fp m_visc_n;
vector_fp m_visc_Tact;
//! Polynomial coefficients of the conductivities
/*!
* These express the temperature dependendence of the pures
* species conductivities
//! Molecular interaction energies associated with viscosity
/**
* These multiply the viscosity according to
* \f[ exp( \sum_{i} \sum_{j} X_i X_j E_{i,j} / T \f].
*/
vector<vector_fp> m_condcoeffs;
DenseMatrix m_visc_Eij;
//! Molecular interaction entropies associated with viscosity
/**
* These multiply the viscosity according to
* \f[ exp( \sum_{i} \sum{j} X_i X_j S_{i,j} \f].
*/
DenseMatrix m_visc_Sij;
//! Pure species thermal conductivities in Arrhenius temperature-dependent form.
vector_fp m_thermCond_A;
vector_fp m_thermCond_n;
vector_fp m_thermCond_Tact;
//! Species hydrodynamic radius
vector_fp m_hydrodynamic_radius;
//! Polynomial coefficients of the binary diffusion coefficients
/*!
@ -389,7 +428,6 @@ namespace Cantera {
* added.
*/
vector<vector_fp> m_diffcoeffs;
//! Internal value of the gradient of the mole fraction vector
/*!
@ -472,9 +510,9 @@ namespace Cantera {
*/
DenseMatrix m_bdiff;
//! Species viscosities
//! Species viscosities and their logarithm
/*!
* Viscosity of the species
* Viscosity of the species and its logarithm
* Length = number of species
*
* Depends on the temperature. We have set the pressure dependence
@ -482,19 +520,8 @@ namespace Cantera {
*
* controlling update boolean -> m_visc_temp_ok
*/
vector_fp viscSpecies_;
//! Sqrt of the species viscosities
/*!
* The sqrt(visc) is used in the mixing formulas
* Length = m_nsp
*
* Depends on the temperature and perhaps pressure, but
* not the species concentrations
*
* controlling update boolean m_visc_temp_ok
*/
vector_fp m_sqvisc;
vector_fp m_viscSpecies;
vector_fp m_logViscSpecies;
//! Internal value of the species individual thermal conductivities
/*!
@ -505,10 +532,7 @@ namespace Cantera {
*
* controlling update boolean -> m_cond_temp_ok
*/
vector_fp m_cond;
//! Polynomials of the log of the temperature
vector_fp m_polytempvec;
vector_fp m_condSpecies;
//! State of the mole fraction vector.
int m_iStateMF;
@ -592,31 +616,12 @@ namespace Cantera {
*/
int viscosityModel_;
//! viscosity weighting functions
DenseMatrix m_phi;
//! Matrix of the ratios of the species molecular weights
/*!
* m_wratjk(i,j) = (m_mw[j]/m_mw[k])**0.25
*/
DenseMatrix m_wratjk;
//! Matrix of the ratios of the species molecular weights
/*!
* m_wratkj1(i,j) = (1.0 + m_mw[k]/m_mw[j])**0.5
*/
DenseMatrix m_wratkj1;
//! RHS to the stefan-maxwell equation
DenseMatrix m_B;
//! Matrix for the stefan maxwell equation.
DenseMatrix m_A;
//! Internal storage for the species LJ well depth
vector_fp m_eps;
//! Current Temperature -> locally storred
/*!
* This is used to test whether new temperature computations
@ -630,21 +635,6 @@ namespace Cantera {
//! Current value of kT
doublereal m_kbt;
//! Current Temperature **0.5
doublereal m_sqrt_t;
//! Current Temperature **0.25
doublereal m_t14;
//! Current Temperature **1.5
doublereal m_t32;
//! Current temperature function
/*!
* This is equal to sqrt(Boltzmann * T)
*/
doublereal m_sqrt_kbt;
//! Current value of the pressure
doublereal m_press;
@ -745,6 +735,17 @@ namespace Cantera {
* Either 1, 2, or 3
*/
int m_nDim;
private:
//! Throw an exception if this method is invoked.
/*!
* This probably indicates something is not yet implemented.
*
* @pram msg Indicates the member function which is not implemented
*/
doublereal err(std::string msg) const;
};
}
#endif

80
Cantera/src/transport/LiquidTransportData.h Normal file
View file

@ -0,0 +1,80 @@
/**
* @file TransportFactory.h
* Header file defining class TransportFactory
* (see \link Cantera::TransportFactory TransportFactory\endlink)
*/
/*
* $Author: hkmoffa $
* $Date: 2008/12/24 18:19:01 $
* $Revision: 1.14 $
*
*
*
*/
#ifndef CT_LIQUIDTRANSPORTDATA_H
#define CT_LIQUIDTRANSPORTDATA_H
// STL includes
#include <vector>
#include <string>
#include <iostream>
#include <new>
// Cantera includes
#include "ct_defs.h"
#include "TransportBase.h"
#include "FactoryBase.h"
namespace Cantera {
enum LiquidTR_Model {
LTR_MODEL_NOTSET=-1,
LTR_MODEL_CONSTANT,
LTR_MODEL_ARRHENIUS,
LTR_MODEL_COEFF
};
class LiquidTransportData {
public:
LiquidTransportData() :
speciesName("-"),
model_hydroradius(LTR_MODEL_NOTSET),
hydroradius(-1.0),
model_viscosity(LTR_MODEL_NOTSET),
model_thermalCond(LTR_MODEL_NOTSET),
model_speciesDiffusivity(LTR_MODEL_NOTSET)
{
}
std::string speciesName;
//! Model type for the hydroradius
LiquidTR_Model model_hydroradius;
//! Actual value of the hydroradius
doublereal hydroradius;
//! Model type for the hydroradius
LiquidTR_Model model_viscosity;
vector_fp viscCoeffs;
//! Model type for the hydroradius
LiquidTR_Model model_thermalCond;
vector_fp thermalCondCoeffs;
//! Model type for the hydroradius
LiquidTR_Model model_speciesDiffusivity;
vector_fp speciesDiffusivityCoeffs;
};
}
#endif

75
Cantera/src/transport/LiquidTransportParams.h
View file

@ -5,27 +5,62 @@
#include "ct_defs.h"
#include "TransportBase.h"
#include "TransportParams.h"
#include "LiquidTransportData.h"
#include "xml.h"
#include "XML_Writer.h"
namespace Cantera {
/**
*
* Holds transport data. Used by TransportFactory.
*
* Holds transport model parameters relevant to transport in
* liquids for which activated jump processes limit transport
* (giving Arrhenius type transport properties).
* Used by TransportFactory.
*/
class LiquidTransportParams {
class LiquidTransportParams :public TransportParams {
public:
LiquidTransportParams() : thermo(0), xml(0) {}
virtual ~LiquidTransportParams();
int nsp;
LiquidTransportParams() {}
~LiquidTransportParams() {}
// phase_t* mix;
thermo_t* thermo;
vector_fp mw;
//section for liquid transport properties
//Arrhenius parameters for transport coefficients:
//!Arrhenius pre-exponential parameter for viscosity.
vector_fp visc_A;
//!Temperature exponent for viscosity.
vector_fp visc_n;
//!Arrhenius activation temperature for viscosity.
vector_fp visc_Tact;
//!Arrhenius pre-exponential parameter for thermal conductivity.
vector_fp thermCond_A;
//!Temperature exponent for thermal conductivity.
vector_fp thermCond_n;
//!Arrhenius activation temperature for thermal conductivity.
vector_fp thermCond_Tact;
//! Energies of molecular interaction associated with viscosity.
/**
* These multiply the mixture viscosity by
* \f[ \exp( \sum_{i} \sum_{j} X_i X_j ( S_{i,j} + E_{i,j} / T ) ) \f].
*
* The overall formula for the logarithm of the mixture viscosity is
*
* \f[ \ln \eta_{mix} = \sum_i X_i \ln \eta_i
* + \sum_i \sum_j X_i X_j ( S_{i,j} + E_{i,j} / T ) \f].
*/
DenseMatrix visc_Eij;
//! Entropies of molecular interaction associated with viscosity.
DenseMatrix visc_Sij;
//Hydrodynamic radius of transported molecule
vector_fp hydroRadius;
//! Coefficients for the limiting conductivity of ions
//! in solution: A_k
@ -47,26 +82,8 @@ namespace Cantera {
vector_fp B_k_cond;
// polynomial fits
std::vector<vector_fp> viscCoeffsVector_;
std::vector<vector_fp> condcoeffs;
std::vector<vector_fp> diffcoeffs ;
std::vector<Cantera::LiquidTransportData> LTData;
std::vector<bool> polar;
//vector_fp alpha;
vector_fp fitlist;
vector_fp eps;
vector_fp sigma;
DenseMatrix reducedMass;
DenseMatrix diam;
DenseMatrix epsilon;
DenseMatrix dipole;
DenseMatrix delta;
doublereal tmax, tmin;
int mode;
XML_Writer* xml;
int log_level;
};
}

6
Cantera/src/transport/Makefile.in
View file

@ -35,11 +35,13 @@ CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG)
# Base Transport Object Files
TRAN_OBJ = TransportFactory.o MultiTransport.o MixTransport.o MMCollisionInt.o \
SolidTransport.o DustyGasTransport.o TransportBase.o WaterTransport.o
SolidTransport.o DustyGasTransport.o TransportBase.o WaterTransport.o \
SimpleTransport.o
TRAN_H = TransportFactory.h MultiTransport.h MixTransport.h \
MMCollisionInt.h SolidTransport.h DustyGasTransport.h \
TransportBase.h L_matrix.h TransportParams.h WaterTransport.h
TransportBase.h L_matrix.h TransportParams.h WaterTransport.h \
SimpleTransport.h LiquidTransportData.h
ifeq ($(do_electro),1)
do_issp = 1

4
Cantera/src/transport/MixTransport.cpp
View file

@ -50,7 +50,7 @@ namespace Cantera {
}
bool MixTransport::init(TransportParams& tr) {
bool MixTransport::initGas( GasTransportParams& tr ) {
// constant substance attributes
m_thermo = tr.thermo;
@ -72,7 +72,7 @@ namespace Cantera {
m_zrot = tr.zrot;
m_crot = tr.crot;
m_epsilon = tr.epsilon;
m_mode = tr.mode;
m_mode = tr.mode_;
m_diam = tr.diam;
m_eps = tr.eps;
m_alpha = tr.alpha;

6
Cantera/src/transport/MixTransport.h Executable file → Normal file
View file

@ -37,7 +37,7 @@ using namespace std;
namespace Cantera {
class TransportParams;
class GasTransportParams;
/**
* Class MixTransport implements mixture-averaged transport
@ -51,7 +51,7 @@ namespace Cantera {
virtual ~MixTransport() {}
virtual int model() { return cMixtureAveraged; }
virtual int model() const { return cMixtureAveraged; }
//! Viscosity of the mixture
/*!
@ -120,7 +120,7 @@ namespace Cantera {
* @param tr Transport parameters for all of the species
* in the phase.
*/
virtual bool init(TransportParams& tr);
virtual bool initGas( GasTransportParams& tr );
friend class TransportFactory;

4
Cantera/src/transport/MultiTransport.cpp
View file

@ -142,7 +142,7 @@ namespace Cantera {
}
bool MultiTransport::init(TransportParams& tr) {
bool MultiTransport::initGas( GasTransportParams& tr ) {
// constant mixture attributes
//m_phase = tr.mix;
@ -167,7 +167,7 @@ namespace Cantera {
m_zrot = tr.zrot;
m_crot = tr.crot;
m_epsilon = tr.epsilon;
m_mode = tr.mode;
m_mode = tr.mode_;
m_diam = tr.diam;
m_eps = tr.eps;
m_alpha = tr.alpha;

6
Cantera/src/transport/MultiTransport.h Executable file → Normal file
View file

@ -40,7 +40,7 @@ namespace Cantera {
TRANSOLVE_LU
};
class TransportParams;
class GasTransportParams;
/////////////////////////////////////////////////////////////
@ -85,7 +85,7 @@ namespace Cantera {
virtual ~MultiTransport();
// overloaded base class methods
virtual int model() {
virtual int model() const {
if (m_mode == CK_Mode)
return CK_Multicomponent;
else
@ -162,7 +162,7 @@ namespace Cantera {
/**
* @internal
*/
virtual bool init(TransportParams& tr);
virtual bool initGas( GasTransportParams& tr );
/**

758
Cantera/src/transport/SimpleTransport.cpp Normal file
View file

@ -0,0 +1,758 @@
/**
* @file SimpleTransport.cpp
* Simple mostly constant transport properties
*/
/*
* $Revision: 1.10 $
* $Date: 2009/03/24 20:44:30 $
*/
#include "ThermoPhase.h"
#include "SimpleTransport.h"
#include "utilities.h"
#include "LiquidTransportParams.h"
#include "TransportFactory.h"
#include "ctlapack.h"
#include <iostream>
using namespace std;
/**
* Mole fractions below MIN_X will be set to MIN_X when computing
* transport properties.
*/
#define MIN_X 1.e-14
namespace Cantera {
//================================================================================================
SimpleTransport::SimpleTransport(thermo_t* thermo, int ndim) :
Transport(thermo, ndim),
m_nsp(0),
tempDepType_(0),
compositionDepType_(0),
useHydroRadius_(false),
doMigration_(0),
m_tmin(-1.0),
m_tmax(100000.),
m_iStateMF(-1),
concTot_(0.0),
m_temp(-1.0),
m_press(-1.0),
m_lambda(-1.0),
m_viscmix(-1.0),
m_visc_mix_ok(false),
m_visc_temp_ok(false),
m_diff_mix_ok(false),
m_diff_temp_ok(false),
m_cond_temp_ok(false),
m_cond_mix_ok(false),
m_nDim(1)
{
}
//================================================================================================
SimpleTransport::SimpleTransport(const SimpleTransport &right) :
Transport(),
m_nsp(0),
tempDepType_(0),
compositionDepType_(0),
useHydroRadius_(false),
doMigration_(0),
m_tmin(-1.0),
m_tmax(100000.),
m_iStateMF(-1),
m_temp(-1.0),
m_press(-1.0),
m_lambda(-1.0),
m_viscmix(-1.0),
m_visc_mix_ok(false),
m_visc_temp_ok(false),
m_diff_mix_ok(false),
m_diff_temp_ok(false),
m_cond_temp_ok(false),
m_cond_mix_ok(false),
m_nDim(1)
{
/*
* Use the assignment operator to do the brunt
* of the work for the copy construtor.
*/
*this = right;
}
//================================================================================================
SimpleTransport& SimpleTransport::operator=(const SimpleTransport& right) {
if (&right != this) {
return *this;
}
Transport::operator=(right);
m_nsp = right.m_nsp;
tempDepType_ = right.tempDepType_;
compositionDepType_ = right.compositionDepType_;
useHydroRadius_ = right.useHydroRadius_;
doMigration_ = right.doMigration_;
m_tmin = right.m_tmin;
m_tmax = right.m_tmax;
m_mw = right.m_mw;
m_coeffVisc_Ns = right.m_coeffVisc_Ns;
m_coeffLambda_Ns = right.m_coeffLambda_Ns;
m_coeffDiff_Ns = right.m_coeffDiff_Ns;
m_Grad_X = right.m_Grad_X;
m_Grad_T = right.m_Grad_T;
m_Grad_P = right.m_Grad_P;
m_Grad_V = right.m_Grad_V;
m_diffSpecies = right.m_diffSpecies;
m_viscSpecies = right.m_viscSpecies;
m_condSpecies = right.m_condSpecies;
m_iStateMF = -1;
m_molefracs = right.m_molefracs;
m_concentrations = right.m_concentrations;
concTot_ = right.concTot_;
meanMolecularWeight_ = right.meanMolecularWeight_;
dens_ = right.dens_;
m_chargeSpecies = right.m_chargeSpecies;
m_temp = right.m_temp;
m_press = right.m_press;
m_lambda = right.m_lambda;
m_viscmix = right.m_viscmix;
m_spwork = right.m_spwork;
m_visc_mix_ok = false;
m_visc_temp_ok = false;
m_diff_mix_ok = false;
m_diff_temp_ok = false;
m_cond_temp_ok = false;
m_cond_mix_ok = false;
m_nDim = right.m_nDim;
return *this;
}
//================================================================================================
Transport *SimpleTransport::duplMyselfAsTransport() const {
SimpleTransport* tr = new SimpleTransport(*this);
return (dynamic_cast<Transport *>(tr));
}
//================================================================================================
// Initialize the object
/*
* This is where we dimension everything.
*/
bool SimpleTransport::initLiquid(LiquidTransportParams& tr) {
int k;
// constant substance attributes
m_thermo = tr.thermo;
m_nsp = m_thermo->nSpecies();
m_tmin = m_thermo->minTemp();
m_tmax = m_thermo->maxTemp();
// make a local copy of the molecular weights
m_mw.resize(m_nsp);
copy(m_thermo->molecularWeights().begin(),
m_thermo->molecularWeights().end(), m_mw.begin());
/*
* Get the input Viscosities
*/
m_viscSpecies.resize(m_nsp);
m_coeffVisc_Ns.clear();
m_coeffVisc_Ns.resize(m_nsp);
Cantera::LiquidTransportData &ltd0 = tr.LTData[0];
LiquidTR_Model vm0 = ltd0.model_viscosity;
if (vm0 == LTR_MODEL_CONSTANT) {
tempDepType_ = 0;
} else if (vm0 == LTR_MODEL_ARRHENIUS) {
tempDepType_ = 1;
} else if (vm0 == LTR_MODEL_NOTSET) {
throw CanteraError("SimpleTransport::initLiquid",
"Viscosity Model is not set in the input file");
} else {
throw CanteraError("SimpleTransport::initLiquid",
"Viscosity Model is not handled by this object");
}
for (k = 0; k < m_nsp; k++) {
Cantera::LiquidTransportData &ltd = tr.LTData[k];
LiquidTR_Model vm = ltd.model_viscosity;
if (vm != vm0) {
throw CanteraError(" SimpleTransport::initLiquid",
"different viscosity models");
}
vector_fp &kentry = m_coeffVisc_Ns[k];
kentry = ltd.viscCoeffs;
}
/*
* Get the input thermal conductivities
*/
m_condSpecies.resize(m_nsp);
m_coeffLambda_Ns.clear();
m_coeffLambda_Ns.resize(m_nsp);
LiquidTR_Model cm0 = ltd0.model_thermalCond;
if (cm0 != vm0) {
throw CanteraError("SimpleTransport::initLiquid",
"Conductivity model is not the same as the viscosity model");
}
for (k = 0; k < m_nsp; k++) {
Cantera::LiquidTransportData &ltd = tr.LTData[k];
LiquidTR_Model cm = ltd.model_thermalCond;
if (cm != cm0) {
throw CanteraError(" SimpleTransport::initLiquid",
"different thermal conductivity models");
}
vector_fp &kentry = m_coeffLambda_Ns[k];
kentry = ltd.thermalCondCoeffs;
}
/*
* Get the input species diffusivities
*/
useHydroRadius_ = false;
m_diffSpecies.resize(m_nsp);
m_coeffDiff_Ns.clear();
m_coeffDiff_Ns.resize(m_nsp);
LiquidTR_Model dm0 = ltd0.model_speciesDiffusivity;
if (dm0 != vm0) {
if (dm0 == LTR_MODEL_NOTSET) {
LiquidTR_Model rm0 = ltd0.model_hydroradius;
if (rm0 != vm0) {
throw CanteraError("SimpleTransport::initLiquid",
"hydroradius model is not the same as the viscosity model");
} else {
useHydroRadius_ = true;
}
}
for (k = 0; k < m_nsp; k++) {
Cantera::LiquidTransportData &ltd = tr.LTData[k];
LiquidTR_Model dm = ltd.model_speciesDiffusivity;
if (dm == LTR_MODEL_NOTSET) {
LiquidTR_Model rm = ltd.model_hydroradius;
if (rm != vm0) {
throw CanteraError("SimpleTransport::initLiquid",
"hydroradius model is not the same as the viscosity model");
}
if (rm != LTR_MODEL_CONSTANT) {
throw CanteraError("SimpleTransport::initLiquid",
"hydroradius model is not constant");
}
vector_fp &kentry = m_coeffHydroRadius_Ns[k];
kentry.push_back(ltd.hydroradius);
} else {
if (dm != dm0) {
throw CanteraError(" SimpleTransport::initLiquid",
"different thermal conductivity models");
}
vector_fp &kentry = m_coeffDiff_Ns[k];
kentry = ltd.speciesDiffusivityCoeffs;
}
}
}
m_molefracs.resize(m_nsp);
m_concentrations.resize(m_nsp);
m_chargeSpecies.resize(m_nsp);
for (k = 0; k < m_nsp; k++) {
m_chargeSpecies[k] = m_thermo->charge(k);
}
m_spwork.resize(m_nsp);
// resize the internal gradient variables
m_Grad_X.resize(m_nDim * m_nsp, 0.0);
m_Grad_T.resize(m_nDim, 0.0);
m_Grad_P.resize(m_nDim, 0.0);
m_Grad_V.resize(m_nDim, 0.0);
// set all flags to false
m_visc_mix_ok = false;
m_visc_temp_ok = false;
m_cond_temp_ok = false;
m_cond_mix_ok = false;
m_diff_temp_ok = false;
m_diff_mix_ok = false;
return true;
}
//================================================================================================
// Returns the mixture viscosity of the solution
/*
* The viscosity is computed using the general mixture rules
* specified in the variable compositionDepType_.
*
* Solvent-only:
* \f[
* \mu = \mu_0
* \f]
* Mixture-average:
* \f[
* \mu = \sum_k {\mu_k X_k}
* \f]
*
* Here \f$ \mu_k \f$ is the viscosity of pure species \e k.
*
* @see updateViscosity_T();
*/
doublereal SimpleTransport::viscosity() {
update_T();
update_C();
if (m_visc_mix_ok) return m_viscmix;
// update m_viscSpecies[] if necessary
if (!m_visc_temp_ok) {
updateViscosity_T();
}
if (compositionDepType_ == 0) {
m_viscmix = m_viscSpecies[0];
} else if (compositionDepType_ == 1) {
m_viscmix = 0.0;
for (int k = 0; k < m_nsp; k++) {
m_viscmix += m_viscSpecies[k] * m_molefracs[k];
}
}
m_visc_mix_ok = true;
return m_viscmix;
}
//================================================================================================
void SimpleTransport::getSpeciesViscosities(doublereal* visc) {
update_T();
if (!m_visc_temp_ok) {
updateViscosity_T();
}
copy(m_viscSpecies.begin(), m_viscSpecies.end(), visc);
}
//================================================================================================
void SimpleTransport::getBinaryDiffCoeffs(int ld, doublereal* d) {
int i, j;
double bdiff;
update_T();
// if necessary, evaluate the species diffusion coefficents
// from the polynomial fits
if (!m_diff_temp_ok) updateDiff_T();
for (i = 0; i < m_nsp; i++) {
for (j = 0; j < m_nsp; j++) {
bdiff = 0.5 * (m_diffSpecies[i] + m_diffSpecies[j]);
d[i*m_nsp+j] = bdiff;
}
}
}
//================================================================================================
// Get the electrical Mobilities (m^2/V/s).
/*
* This function returns the mobilities. In some formulations
* this is equal to the normal mobility multiplied by faraday's constant.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^e_k = \frac{F D_k}{R T}
* \f]
*
* @param mobil_e Returns the mobilities of
* the species in array \c mobil_e. The array must be
* dimensioned at least as large as the number of species.
*/
void SimpleTransport::getMobilities(doublereal* const mobil) {
int k;
getMixDiffCoeffs(DATA_PTR(m_spwork));
doublereal c1 = ElectronCharge / (Boltzmann * m_temp);
for (k = 0; k < m_nsp; k++) {
mobil[k] = c1 * m_spwork[k];
}
}
//================================================================================================
// Get the fluid mobilities (s kmol/kg).
/*
* This function returns the fluid mobilities. Usually, you have
* to multiply Faraday's constant into the resulting expression
* to general a species flux expression.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^f_k = \frac{D_k}{R T}
* \f]
*
*
* @param mobil_f Returns the mobilities of
* the species in array \c mobil. The array must be
* dimensioned at least as large as the number of species.
*/
void SimpleTransport::getFluidMobilities(doublereal* const mobil_f) {
int k;
getMixDiffCoeffs(DATA_PTR(m_spwork));
doublereal c1 = 1.0 / (GasConstant * m_temp);
for (k = 0; k < m_nsp; k++) {
mobil_f[k] = c1 * m_spwork[k];
}
}
//================================================================================================
void SimpleTransport::set_Grad_V(const doublereal* const grad_V) {
doMigration_ = false;
for (int a = 0; a < m_nDim; a++) {
m_Grad_V[a] = grad_V[a];
if (fabs(grad_V[a]) > 1.0E-13) doMigration_ = true;
}
}
//================================================================================================
void SimpleTransport::set_Grad_T(const doublereal* const grad_T) {
for (int a = 0; a < m_nDim; a++) {
m_Grad_T[a] = grad_T[a];
}
}
//================================================================================================
void SimpleTransport::set_Grad_X(const doublereal* const grad_X) {
int itop = m_nDim * m_nsp;
for (int i = 0; i < itop; i++) {
m_Grad_X[i] = grad_X[i];
}
}
//================================================================================================
// Returns the mixture thermal conductivity of the solution
/*
* The thermal is computed using the general mixture rules
* specified in the variable compositionDepType_.
*
* Solvent-only:
* \f[
* \lambda = \lambda_0
* \f]
* Mixture-average:
* \f[
* \lambda = \sum_k {\lambda_k X_k}
* \f]
*
* Here \f$ \lambda_k \f$ is the thermal conductivity of pure species \e k.
*
* @see updateCond_T();
*/
doublereal SimpleTransport::thermalConductivity() {
update_T();
update_C();
if (!m_cond_temp_ok) {
updateCond_T();
}
if (!m_cond_mix_ok) {
if (compositionDepType_ == 0) {
m_lambda = m_condSpecies[0];
} else if (compositionDepType_ == 1) {
m_lambda = 0.0;
for (int k = 0; k < m_nsp; k++) {
m_lambda += m_condSpecies[k] * m_molefracs[k];
}
}
m_cond_mix_ok = true;
}
return m_lambda;
}
//================================================================================================
/*
* Thermal diffusion is not considered in this mixture-averaged
* model. To include thermal diffusion, use transport manager
* MultiTransport instead. This methods fills out array dt with
* zeros.
*/
void SimpleTransport::getThermalDiffCoeffs(doublereal* const dt) {
for (int k = 0; k < m_nsp; k++) {
dt[k] = 0.0;
}
}
//================================================================================================
/**
* @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
*
* \f[
* \vec{j}_k = -n M_k D_k \nabla X_k.
* \f]
*/
void SimpleTransport::getSpeciesFluxes(int ndim,
const doublereal* grad_T,
int ldx, const doublereal* grad_X,
int ldf, doublereal* fluxes) {
set_Grad_T(grad_T);
set_Grad_X(grad_X);
getSpeciesFluxesExt(ldf, fluxes);
}
//================================================================================================
// Return the species diffusive mass fluxes wrt to
// the mass averaged velocity.
/*
*
* units = kg/m2/s
*
* Internally, gradients in the in mole fraction, temperature
* and electrostatic potential contribute to the diffusive flux
*
*
* The diffusive mass flux of species \e k is computed from the following
* formula
*
* \f[
* j_k = - M_k z_k u^f_k F c_k \nabla \Psi - c M_k D_k \nabla X_k - Y_k V_c
* \f]
*
* where V_c is the correction velocity
*
* \f[
* V_c = - \sum_j {M_k z_k u^f_k F c_k \nabla \Psi + c M_j D_j \nabla X_j}
* \f]
*
* @param ldf stride of the fluxes array. Must be equal to
* or greater than the number of species.
* @param fluxes Vector of calculated fluxes
*/
void SimpleTransport::getSpeciesFluxesExt(int ldf, doublereal* fluxes) {
int n, k;
AssertThrow(ldf >= m_nsp ,"SimpleTransport::getSpeciesFluxesExt: Stride must be greater than m_nsp");
update_T();
update_C();
getMixDiffCoeffs(DATA_PTR(m_spwork));
const array_fp& mw = m_thermo->molecularWeights();
const doublereal* y = m_thermo->massFractions();
doublereal conc = m_thermo->molarDensity();
// Unroll wrt ndim
vector_fp sum(m_nDim, 0.0);
if (doMigration_) {
for (n = 0; n < m_nDim; n++) {
for (k = 0; k < m_nsp; k++) {
fluxes[n*ldf + k] = -conc * mw[k] * m_spwork[k] * m_Grad_X[n*m_nsp + k];
sum[n] += fluxes[n*ldf + k];
}
}
} else {
double FRT = ElectronCharge / (Boltzmann * m_temp);
for (n = 0; n < m_nDim; n++) {
for (k = 0; k < m_nsp; k++) {
fluxes[n*ldf + k] = -conc * mw[k] * m_spwork[k] *
( m_Grad_X[n*m_nsp + k] + FRT * m_molefracs[k] * m_chargeSpecies[k] * m_Grad_V[n*m_nsp + k]);
sum[n] += fluxes[n*ldf + 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*ldf + k] -= y[k]*sum[n];
}
}
}
//================================================================================================
// Mixture-averaged diffusion coefficients [m^2/s].
/*
* Returns the simple diffusion coefficients input into the model. Nothing fancy here.
*/
void SimpleTransport::getMixDiffCoeffs(doublereal* const d) {
update_T();
update_C();
// update the binary diffusion coefficients if necessary
if (!m_diff_temp_ok) {
updateDiff_T();
}
for (int k = 0; k < m_nsp; k++) {
d[k] = m_diffSpecies[k];
}
}
//================================================================================================
// Handles the effects of changes in the mixture concentration
/*
* This is called for every interface call to check whether
* the concentrations have changed. Concentrations change
* whenever the pressure or the mole fraction has changed.
* If it has changed, the recalculations should be done.
*
* Note this should be a lightweight function since it's
* part of all of the interfaces.
*
* @internal
*/
bool SimpleTransport::update_C() {
// If the pressure has changed then the concentrations
// have changed.
doublereal pres = m_thermo->pressure();
bool qReturn = true;
if (pres != m_press) {
qReturn = false;
m_press = pres;
}
int iStateNew = m_thermo->stateMFNumber();
if (iStateNew != m_iStateMF) {
qReturn = false;
m_thermo->getMoleFractions(DATA_PTR(m_molefracs));
m_thermo->getConcentrations(DATA_PTR(m_concentrations));
concTot_ = 0.0;
for (int k = 0; k < m_nsp; k++) {
m_molefracs[k] = fmaxx(0.0, m_molefracs[k]);
concTot_ += m_concentrations[k];
}
dens_ = m_thermo->density();
meanMolecularWeight_ = m_thermo->meanMolecularWeight();
}
if (qReturn) {
return false;
}
// Mixture stuff needs to be evaluated
m_visc_mix_ok = false;
m_diff_mix_ok = false;
m_cond_mix_ok = false;
return true;
}
//================================================================================================
/**
* Update the temperature-dependent parts of the mixture-averaged
* thermal conductivity.
*/
void SimpleTransport::updateCond_T() {
int k;
if (tempDepType_ == 0) {
for (k = 0; k < m_nsp; k++) {
Coeff_T_ &coeff = m_coeffLambda_Ns[k];
m_condSpecies[k] = coeff[0];
}
} else if (tempDepType_ == 1) {
for (k = 0; k < m_nsp; k++) {
Coeff_T_ &coeff = m_coeffLambda_Ns[k];
m_condSpecies[k] = coeff[0] * pow(m_temp,coeff[1]) * exp(-coeff[2]/m_temp);
}
}
m_cond_temp_ok = true;
m_cond_mix_ok = false;
}
//================================================================================================
/**
* Update the species diffusion coefficients.
*/
void SimpleTransport::updateDiff_T() {
int k;
if (useHydroRadius_) {
if (tempDepType_ == 0) {
for (k = 0; k < m_nsp; k++) {
Coeff_T_ &coeff = m_coeffDiff_Ns[k];
m_diffSpecies[k] = coeff[0];
}
} else if (tempDepType_ == 1) {
for (k = 0; k < m_nsp; k++) {
Coeff_T_ &coeff = m_coeffDiff_Ns[k];
m_diffSpecies[k] = coeff[0] * pow(m_temp,coeff[1]) * exp(-coeff[2]/m_temp);
}
}
} else {
double visc = viscosity();
double RT = GasConstant * m_temp;
for (k = 0; k < m_nsp; k++) {
Coeff_T_ &coeff = m_coeffHydroRadius_Ns[k];
double rad = coeff[0];
m_diffSpecies[k] = RT / (6.0 * Pi * visc * rad);
}
}
m_diff_temp_ok = true;
m_diff_mix_ok = false;
}
//================================================================================================
/**
* Update the pure-species viscosities.
*/
void SimpleTransport::updateViscosities_C() {
}
//================================================================================================
/**
* Update the temperature-dependent viscosity terms.
* Updates the array of pure species viscosities, and the
* weighting functions in the viscosity mixture rule.
* The flag m_visc_ok is set to true.
*/
void SimpleTransport::updateViscosity_T() {
int k;
if (tempDepType_ == 0) {
for (k = 0; k < m_nsp; k++) {
Coeff_T_ &coeff = m_coeffVisc_Ns[k];
m_viscSpecies[k] = coeff[0];
}
} else if (tempDepType_ == 1) {
for (k = 0; k < m_nsp; k++) {
Coeff_T_ &coeff = m_coeffVisc_Ns[k];
m_viscSpecies[k] = coeff[0] * pow(m_temp,coeff[1]) * exp(-coeff[2]/m_temp);
}
}
m_visc_temp_ok = true;
m_visc_mix_ok = false;
}
//=================================================================================================
bool SimpleTransport::update_T()
{
doublereal t = m_thermo->temperature();
if (t == m_temp) return false;
if (t < 0.0) {
throw CanteraError("SimpleTransport::update_T",
"negative temperature "+fp2str(t));
}
// Compute various functions of temperature
m_temp = t;
// temperature has changed, so polynomial temperature
// interpolations will need to be reevaluated.
// Set all of these flags to false
m_visc_mix_ok = false;
m_visc_temp_ok = false;
m_cond_temp_ok = true;
m_cond_mix_ok = false;
m_diff_mix_ok = false;
m_diff_temp_ok = false;
return true;
}
//================================================================================================
/**
* Throw an exception if this method is invoked.
* This probably indicates something is not yet implemented.
*/
doublereal SimpleTransport::err(std::string msg) const {
throw CanteraError("SimpleTransport Class",
"\n\n\n**** Method "+ msg +" not implemented in model "
+ int2str(model()) + " ****\n"
"(Did you forget to specify a transport model?)\n\n\n");
return 0.0;
}
//================================================================================================
}
//================================================================================================

718
Cantera/src/transport/SimpleTransport.h Normal file
View file

@ -0,0 +1,718 @@
/**
*
* @file SimpleTransport.h
* Header file defining class SimpleTransport
*/
/*
* $Revision: 1.9 $
* $Date: 2009/03/27 18:24:39 $
*/
#ifndef CT_SIMPLETRAN_H
#define CT_SIMPLETRAN_H
// STL includes
#include <vector>
#include <string>
#include <map>
#include <numeric>
#include <algorithm>
using namespace std;
// Cantera includes
#include "TransportBase.h"
#include "DenseMatrix.h"
#include "TransportParams.h"
#include "LiquidTransportParams.h"
namespace Cantera {
class LiquidTransportParams;
//! Class LiquidTransport implements mixture-averaged transport
//! properties for liquid phases.
/*!
* The model is based on that
* described by Newman, Electrochemical Systems
*
* The velocity of species i may be described by the
* following equation p. 297 (12.1)
*
* \f[
* c_i \nabla \mu_i = R T \sum_j \frac{c_i c_j}{c_T D_{ij}}
* (\mathbf{v}_j - \mathbf{v}_i)
* \f]
*
* This as written is degenerate by 1 dof.
*
* To fix this we must add in the definition of the mass averaged
* velocity of the solution. We will call the simple bold-faced
* \f$\mathbf{v} \f$
* symbol the mass-averaged velocity. Then, the relation
* between \f$\mathbf{v}\f$ and the individual species velocities is
* \f$\mathbf{v}_i\f$
*
* \f[
* \rho_i \mathbf{v}_i = \rho_i \mathbf{v} + \mathbf{j}_i
* \f]
* where \f$\mathbf{j}_i\f$ are the diffusional fluxes of species i
* with respect to the mass averaged velocity and
*
* \f[
* \sum_i \mathbf{j}_i = 0
* \f]
*
* and
*
* \f[
* \sum_i \rho_i \mathbf{v}_i = \rho \mathbf{v}
* \f]
*
* Using these definitions, we can write
*
* \f[
* \mathbf{v}_i = \mathbf{v} + \frac{\mathbf{j}_i}{\rho_i}
* \f]
*
*
* \f[
* c_i \nabla \mu_i = R T \sum_j \frac{c_i c_j}{c_T D_{ij}}
* (\frac{\mathbf{j}_j}{\rho_j} - \frac{\mathbf{j}_i}{\rho_i})
* = R T \sum_j \frac{1}{D_{ij}}
* (\frac{x_i \mathbf{j}_j}{M_j} - \frac{x_j \mathbf{j}_i}{M_i})
* \f]
*
* The equations that we actually solve are
*
* \f[
* c_i \nabla \mu_i =
* = R T \sum_j \frac{1}{D_{ij}}
* (\frac{x_i \mathbf{j}_j}{M_j} - \frac{x_j \mathbf{j}_i}{M_i})
* \f]
* and we replace the 0th equation with the following:
*
* \f[
* \sum_i \mathbf{j}_i = 0
* \f]
*
* When there are charged species, we replace the rhs with the
* gradient of the electrochemical potential to obtain the
* modified equation
*
* \f[
* c_i \nabla \mu_i + c_i F z_i \nabla \Phi
* = R T \sum_j \frac{1}{D_{ij}}
* (\frac{x_i \mathbf{j}_j}{M_j} - \frac{x_j \mathbf{j}_i}{M_i})
* \f]
*
* With this formulation we may solve for the diffusion velocities,
* without having to worry about what the mass averaged velocity
* is.
*
* <H2> Viscosity Calculation </H2>
*
* The viscosity calculation may be broken down into two parts.
* In the first part, the viscosity of the pure species are calculated
* In the second part, a mixing rule is applied, based on the
* Wilkes correlation, to yield the mixture viscosity.
*
*
*
*/
class SimpleTransport : public Transport {
public:
typedef vector_fp Coeff_T_;
//! Default constructor.
/*!
* This requires call to initLiquid(LiquidTransportParams& tr)
* after filling LiquidTransportParams to complete instantiation.
* The filling of LiquidTransportParams is currently carried out
* in the TransportFactory class, but might be moved at some point.
*
* @param thermo ThermoPhase object holding species information.
* @param ndim Number of spatial dimensions.
*/
SimpleTransport(thermo_t* thermo = 0, int ndim = 1);
//!Copy Constructor for the %LiquidThermo object.
/*!
* @param right %LiquidTransport to be copied
*/
SimpleTransport(const SimpleTransport &right);
//! Assignment operator
/*!
* This is NOT a virtual function.
*
* @param right Reference to %LiquidTransport object to be copied
* into the current one.
*/
SimpleTransport& operator=(const SimpleTransport& right);
//! Duplication routine for objects which inherit from
//! %Transport
/*!
* This virtual routine can be used to duplicate %Transport objects
* inherited from %Transport even if the application only has
* a pointer to %Transport to work with.
*
* These routines are basically wrappers around the derived copy
* constructor.
*/
virtual Transport *duplMyselfAsTransport() const;
//! virtual destructor
virtual ~SimpleTransport() {}
//! Initialize the transport object
/*!
* Here we change all of the internal dimensions to be sufficient.
* We get the object ready to do property evaluations.
*
* @param tr Transport parameters for all of the species
* in the phase.
*/
virtual bool initLiquid(LiquidTransportParams& tr);
friend class TransportFactory;
//! Return the model id for this transport parameterization
virtual int model() const {
return cSimpleTransport;
}
//! overloaded base class methods
//! Returns the mixture viscosity of the solution
/*!
* The viscosity is computed using the general mixture rules
* specified in the variable compositionDepType_.
*
* Solvent-only:
* \f[
* \mu = \mu_0
* \f]
* Mixture-average:
* \f[
* \mu = \sum_k {\mu_k X_k}
* \f]
*
* Here \f$ \mu_k \f$ is the viscosity of pure species \e k.
*
* @see updateViscosity_T();
*/
virtual doublereal viscosity();
//! Returns the pure species viscosities
/*!
* The pure species viscosities are to be given in an Arrhenius
* form in accordance with activated-jump-process dominated transport.
*/
virtual void getSpeciesViscosities(doublereal* const visc);
//! Returns the binary diffusion coefficients
/*!
* @param ld
* @param d
*/
virtual void getBinaryDiffCoeffs(const int ld, doublereal* const d);
//! Get the Mixture diffusion coefficients
/*!
* @param d vector of mixture diffusion coefficients
* units = m2 s-1. length = number of species
*/
virtual void getMixDiffCoeffs(doublereal* const d);
//! Return the thermal diffusion coefficients
/*!
* These are all zero for this simple implementaion
*
* @param dt thermal diffusion coefficients
*/
virtual void getThermalDiffCoeffs(doublereal* const dt);
//! Returns the mixture thermal conductivity of the solution
/*!
* The thermal is computed using the general mixture rules
* specified in the variable compositionDepType_.
*
* Controlling update boolean = m_condmix_ok
*
* Units are in W/m/K or equivalently kg m / s3 / K
*
* Solvent-only:
* \f[
* \lambda = \lambda_0
* \f]
* Mixture-average:
* \f[
* \lambda = \sum_k {\lambda_k X_k}
* \f]
*
* Here \f$ \lambda_k \f$ is the thermal conductivity of pure species \e k.
*
* @see updateCond_T();
*/
virtual doublereal thermalConductivity();
//! Get the electrical Mobilities (m^2/V/s).
/*!
* This function returns the mobilities. In some formulations
* this is equal to the normal mobility multiplied by faraday's constant.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^e_k = \frac{F D_k}{R T}
* \f]
*
* @param mobil_e Returns the mobilities of
* the species in array \c mobil_e. The array must be
* dimensioned at least as large as the number of species.
*/
virtual void getMobilities(doublereal* const mobil_e);
//! Get the fluid mobilities (s kmol/kg).
/*!
* This function returns the fluid mobilities. Usually, you have
* to multiply Faraday's constant into the resulting expression
* to general a species flux expression.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^f_k = \frac{D_k}{R T}
* \f]
*
*
* @param mobil_f Returns the mobilities of
* the species in array \c mobil. The array must be
* dimensioned at least as large as the number of species.
*/
virtual void getFluidMobilities(doublereal* const mobil_f);
//! Specify the valpdaue of the gradient of the voltage
/*!
*
* @param grad_V Gradient of the voltage (length num dimensions);
*/
virtual void set_Grad_V(const doublereal* const grad_V);
//! Specify the value of the gradient of the temperature
/*!
*
* @param grad_V Gradient of the temperature (length num dimensions);
*/
virtual void set_Grad_T(const doublereal* const grad_T);
//! Specify the value of the gradient of the MoleFractions
/*!
*
* @param grad_X Gradient of the mole fractions(length nsp * num dimensions);
*/
virtual void set_Grad_X(const doublereal* const grad_X);
/**
* @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 getSpeciesFluxes(int ndim,
const doublereal* grad_T,
int ldx, const doublereal* grad_X,
int ldf, doublereal* fluxes);
//! Return the species diffusive mass fluxes wrt to
//! the mass averaged velocity,
/*!
*
* units = kg/m2/s
*
* Internally, gradients in the in mole fraction, temperature
* and electrostatic potential contribute to the diffusive flux
*
*
* The diffusive mass flux of species \e k is computed from the following
* formula
*
* \f[
* j_k = - \rho M_k D_k \nabla X_k - Y_k V_c
* \f]
*
* where V_c is the correction velocity
*
* \f[
* V_c = - \sum_j {\rho M_j D_j \nabla X_j}
* \f]
*
* @param ldf stride of the fluxes array. Must be equal to
* or greater than the number of species.
* @param fluxes Vector of calculated fluxes
*/
virtual void getSpeciesFluxesExt(int ldf, doublereal* fluxes);
protected:
//! Handles the effects of changes in the Temperature, internally
//! within the object.
/*!
* This is called whenever a transport property is requested.
* The first task is to check whether the temperature has changed
* since the last call to update_T().
* If it hasn't then an immediate return is carried out.
*
* @internal
*
* @return Returns true if the temperature has changed, and false otherwise
*/
virtual bool update_T();
//! Handles the effects of changes in the mixture concentration
/*!
* This is called for every interface call to check whether
* the concentrations have changed. Concentrations change
* whenever the pressure or the mole fraction has changed.
* If it has changed, the recalculations should be done.
*
* Note this should be a lightweight function since it's
* part of all of the interfaces.
*
* @internal
*/
virtual bool update_C();
//! Update the temperature-dependent viscosity terms.
//! Updates the array of pure species viscosities, and the
//! weighting functions in the viscosity mixture rule.
/*!
* The flag m_visc_temp_ok is set to true.
*/
void updateViscosity_T();
//! Update the temperature-dependent parts of the mixture-averaged
//! thermal conductivity.
void updateCond_T();
//! Update the concentration parts of the viscosities
/*!
* Internal routine is run whenever the update_boolean
* is false. This routine will calculate
* internal values for the species viscosities.
*
* @internal
*/
void updateViscosities_C();
//! Update the binary diffusion coefficients wrt T.
/*!
* These are evaluated
* from the polynomial fits at unit pressure (1 Pa).
*/
void updateDiff_T();
private:
//! Number of species in the mixture
int m_nsp;
//! Temperature dependence type
/*!
* The following coefficients are allowed to have simple
* temperature dependencies:
* mixture viscosity
* mixture thermal conductivity
* diffusitivy
*
* Types of temperature dependencies:
* 0 - Independent of temperature (only one implemented so far)
* 1 - extended arrhenius form
* 2 - power law form
*/
int tempDepType_;
//! Composition dependence of the transport properties
/*!
* The following coefficients are allowed to have simple
* composition dependencies
* mixture viscosity
* mixture thermal conductivity
*
*
* Types of composition dependencies
* 0 - Solvent values (i.e., species 0) contributes only
* 1 - linear combination of mole fractions;
*/
int compositionDepType_;
bool useHydroRadius_;
//! Boolean indicating whether electro-migration term should be
//! added
/*!
*
*/
bool doMigration_;
//! Minimum temperature applicable to the transport property eval
doublereal m_tmin;
//! Maximum temperature applicable to the transport property evaluator
doublereal m_tmax;
//! Local Copy of the molecular weights of the species
/*!
* Length is Equal to the number of species in the mechanism.
*/
vector_fp m_mw;
//! Pure species viscosities in Arrhenius temperature-dependent form.
std::vector<Coeff_T_> m_coeffVisc_Ns;
//! Pure species thermal conductivities in Arrhenius temperature-dependent form.
/*!
*
*/
std::vector<Coeff_T_> m_coeffLambda_Ns;
//! Pure species viscosities in Arrhenius temperature-dependent form.
std::vector<Coeff_T_> m_coeffDiff_Ns;
std::vector<Coeff_T_> m_coeffHydroRadius_Ns;
//! Internal value of the gradient of the mole fraction vector
/*!
* Note, this is the only gradient value that can and perhaps
* should reflect the true state of the mole fractions in the
* application solution vector. In other words no cropping or
* massaging of the values to make sure they are above zero
* should occur. - developing ....
*
* m_nsp is the number of species in the fluid
* k is the species index
* n is the dimensional index (x, y, or z). It has a length
* equal to m_nDim
*
* m_Grad_X[n*m_nsp + k]
*/
vector_fp m_Grad_X;
//! Internal value of the gradient of the Temperature vector
/*!
* Generally, if a transport property needs this
* in its evaluation it will look to this place
* to get it.
*
* No internal property is precalculated based on gradients.
* Gradients are assumed to be freshly updated before
* every property call.
*/
vector_fp m_Grad_T;
//! Internal value of the gradient of the Pressure vector
/*!
* Generally, if a transport property needs this
* in its evaluation it will look to this place
* to get it.
*
* No internal property is precalculated based on gradients.
* Gradients are assumed to be freshly updated before
* every property call.
*/
vector_fp m_Grad_P;
//! Internal value of the gradient of the Electric Voltage
/*!
* Generally, if a transport property needs this
* in its evaluation it will look to this place
* to get it.
*
* No internal property is precalculated based on gradients.
* Gradients are assumed to be freshly updated before
* every property call.
*/
vector_fp m_Grad_V;
// property values
//! Vector of Species Diffusivities
/*!
* Depends on the temperature. We have set the pressure dependence
* to zero for this liquid phase constituitve model
*
* units m2/s
*/
vector_fp m_diffSpecies;
//! Species viscosities
/*!
* Viscosity of the species
* Length = number of species
*
* Depends on the temperature. We have set the pressure dependence
* to zero for this model
*
* controlling update boolean -> m_visc_temp_ok
*/
vector_fp m_viscSpecies;
//! Internal value of the species individual thermal conductivities
/*!
* Then a mixture rule is applied to get the solution conductivities
*
* Depends on the temperature and perhaps pressure, but
* not the species concentrations
*
* controlling update boolean -> m_cond_temp_ok
*/
vector_fp m_condSpecies;
//! State of the mole fraction vector.
int m_iStateMF;
//! 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 exactly. This mole
* fraction vector comes from the ThermoPhase object. Derivative
* quantities from this are referred to as bounded.
*
* Update info?
* length = m_nsp
*/
vector_fp m_molefracs;
//! Local copy of the concentrations of the species in the phase
/*!
* The concentrations are consistent with the m_molefracs
* vector which is bounded and sums to one.
*
* Update info?
* length = m_nsp
*/
vector_fp m_concentrations;
//! Local copy of the total concentration.
/*!
* This is consistent with the m_concentrations[] and
* m_molefracs[] vector.
*/
doublereal concTot_;
//! Mean molecular weight
doublereal meanMolecularWeight_;
//! Density
doublereal dens_;
//! Local copy of the charge of each species
/*!
* Contains the charge of each species (length m_nsp)
*/
vector_fp m_chargeSpecies;
//! Current Temperature -> locally storred
/*!
* This is used to test whether new temperature computations
* should be performed.
*/
doublereal m_temp;
//! Current value of the pressure
doublereal m_press;
//! Saved value of the mixture thermal conductivity
doublereal m_lambda;
//! Saved value of the mixture viscosity
doublereal m_viscmix;
//! work space
/*!
* Length is equal to m_nsp
*/
vector_fp m_spwork;
private:
//! Boolean indicating that the top-level mixture viscosity is current
/*!
* This is turned false for every change in T, P, or C.
*/
bool m_visc_mix_ok;
//! Boolean indicating that weight factors wrt viscosity is current
bool m_visc_temp_ok;
//! Boolean indicating that mixture diffusion coeffs are current
bool m_diff_mix_ok;
//! Boolean indicating that binary diffusion coeffs are current
bool m_diff_temp_ok;
//! Flag to indicate that the pure species conductivities
//! are current wrt the temperature
bool m_cond_temp_ok;
//! Boolean indicating that mixture conductivity is current
bool m_cond_mix_ok;
//! Number of dimensions
/*!
* Either 1, 2, or 3
*/
int m_nDim;
private:
//! Throw an exception if this method is invoked.
/*!
* This probably indicates something is not yet implemented.
*
* @pram msg Indicates the member function which is not implemented
*/
doublereal err(std::string msg) const;
};
}
#endif

2
Cantera/src/transport/SolidTransport.h
View file

@ -46,7 +46,7 @@ namespace Cantera {
public:
virtual ~SolidTransport() {}
virtual int model() { return cSolidTransport; }
virtual int model() const { return cSolidTransport; }
virtual doublereal thermalConductivity();
virtual void getMixDiffCoeffs(doublereal* const d);

36
Cantera/src/transport/TransportBase.cpp
View file

@ -9,10 +9,12 @@
#include "ThermoPhase.h"
#include "LiquidTransport.h"
#include "ctexceptions.h"
#include "utilities.h"
#include "LiquidTransportParams.h"
#include "TransportFactory.h"
#include "stringUtils.h"
#include "ctlapack.h"
@ -31,8 +33,6 @@ namespace Cantera {
//////////////////// class LiquidTransport methods //////////////
Transport::Transport(thermo_t* thermo, int ndim) :
m_thermo(thermo),
m_ready(false),
@ -109,5 +109,37 @@ namespace Cantera {
{
err("setParameters");
}
void Transport::setThermo(thermo_t& thermo) {
if (!ready()) {
m_thermo = &thermo;
m_nmin = m_thermo->nSpecies();
}
else
throw CanteraError("Transport::setThermo",
"the phase object cannot be changed after "
"the transport manager has been constructed.");
}
doublereal Transport::err(std::string msg) const {
throw CanteraError("Transport Base Class",
"\n\n\n**** Method "+ msg +" not implemented in model "
+ int2str(model()) + " ****\n"
"(Did you forget to specify a transport model?)\n\n\n");
return 0.0;
}
void Transport::finalize() {
if (!ready())
m_ready = true;
else
throw CanteraError("Transport::finalize",
"finalize has already been called.");
}
}

130
Cantera/src/transport/TransportBase.h Executable file → Normal file
View file

@ -28,6 +28,8 @@
namespace Cantera {
class TransportParams;
class GasTransportParams;
class LiquidTransportParams;
const int CK_Mode = 10;
@ -43,6 +45,7 @@ namespace Cantera {
const int cFtnTransport = 600;
const int cLiquidTransport = 700;
const int cAqueousTransport = 750;
const int cSimpleTransport = 770;
const int cRadiativeTransport = 800;
const int cWaterTransport = 721;
@ -107,7 +110,7 @@ namespace Cantera {
* virtual method returns an integer flag that identifies the
* transport model implemented. The base class returns 0.
*/
virtual int model() {return 0;}
virtual int model() const {return 0;}
/**
* Phase object. Every transport manager is designed to compute
@ -157,6 +160,14 @@ namespace Cantera {
virtual doublereal viscosity()
{ return err("viscosity"); }
//! Returns the pure species viscosities
/*!
* The units are Pa-s and the length is the number of species
*
* @param visc Vector of viscosities
*/
virtual void getSpeciesViscosities(doublereal* const visc)
{ err("getSpeciesViscosities"); }
/**
* The bulk viscosity in Pa-s. The bulk viscosity is only
@ -167,9 +178,11 @@ namespace Cantera {
virtual doublereal bulkViscosity()
{ return err("bulkViscosity"); }
/**
* The thermal conductivity in W/m/K.
//! Returns the mixture thermal conductivity in W/m/K.
/*!
* Units are in W / m K or equivalently kg m / s3 K
*
* @return returns thermal conductivity in W/m/K.
*/
virtual doublereal thermalConductivity()
{ return err("thermalConductivity"); }
@ -180,14 +193,48 @@ namespace Cantera {
virtual doublereal electricalConductivity()
{ return err("electricalConductivity"); }
/**
* Electrical mobilities (m^2/V/s). Returns the mobilities of
* the species in array \c mobil. The array must be
* dimensioned at least as large as the number of species.
//! Get the Electrical mobilities (m^2/V/s).
/*!
* This function returns the mobilities. In some formulations
* this is equal to the normal mobility multiplied by faraday's constant.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^e_k = \frac{F D_k}{R T}
* \f]
*
*
* @param mobil_e Returns the mobilities of
* the species in array \c mobil_e. The array must be
* dimensioned at least as large as the number of species.
*/
virtual void getMobilities(doublereal* const mobil)
virtual void getMobilities(doublereal* const mobil_e)
{ err("getMobilities"); }
//! Get the fluid mobilities (s kmol/kg).
/*!
* This function returns the fluid mobilities. Usually, you have
* to multiply Faraday's constant into the resulting expression
* to general a species flux expression.
*
* Frequently, but not always, the mobility is calculated from the
* diffusion coefficient using the Einstein relation
*
* \f[
* \mu^f_k = \frac{D_k}{R T}
* \f]
*
*
* @param mobil_f Returns the mobilities of
* the species in array \c mobil. The array must be
* dimensioned at least as large as the number of species.
*/
virtual void getFluidMobilities(doublereal* const mobil_f)
{ err("getFluidMobilities"); }
//@}
@ -221,6 +268,39 @@ namespace Cantera {
err("getSpeciesFluxes");
}
//! Get the species diffusive mass fluxes wrt to
//! the mass averaged velocity,
//! given the gradients in mole fraction, temperature
//! and electrostatic potential.
/*!
* Units for the returned fluxes are kg m-2 s-1.
*
* @param ndim Number of dimensions in the flux expressions
* @param grad_T Gradient of the temperature
* (length = ndim)
* @param ldx Leading dimension of the grad_X array
* (usually equal to m_nsp but not always)
* @param grad_X Gradients of the mole fraction
* Flat vector with the m_nsp in the inner loop.
* length = ldx * ndim
* @param ldf Leading dimension of the fluxes array
* (usually equal to m_nsp but not always)
* @param grad_Phi Gradients of the electrostatic potential
* (length = ndim)
* @param fluxes Output of the diffusive mass fluxes
* Flat vector with the m_nsp in the inner loop.
* length = ldx * ndim
*/
virtual void getSpeciesFluxesES(int ndim,
const doublereal* grad_T,
int ldx,
const doublereal* grad_X,
int ldf,
const doublereal* grad_Phi,
doublereal* fluxes) {
getSpeciesFluxes( ndim, grad_T, ldx, grad_X, ldf, fluxes );
}
/**
* Get the molar fluxes [kmol/m^2/s], given the thermodynamic
* state at two nearby points.
@ -238,9 +318,9 @@ namespace Cantera {
* Get the mass fluxes [kg/m^2/s], given the thermodynamic
* state at two nearby points.
* @param state1 Array of temperature, density, and mass
* fractions for state 1.
* fractions for state 1.
* @param state2 Array of temperature, density, and mass
* fractions for state 2.
* fractions for state 2.
* @param delta Distance from state 1 to state 2 (m).
*/
virtual void getMassFluxes(const doublereal* state1,
@ -263,9 +343,13 @@ namespace Cantera {
{ err("getThermalDiffCoeffs"); }
/**
* Binary diffusion coefficients [m^2/s].
*/
//! Returns the matrix of binary diffusion coefficients [m^2/s].
/*!
* @param ld Inner stride for writing the two dimension diffusion
* coefficients into a one dimensional vector
* @param d Diffusion coefficient matrix (must be at least m_k * m_k
* in length.
*/
virtual void getBinaryDiffCoeffs(const int ld, doublereal* const d)
{ err("getBinaryDiffCoeffs"); }
@ -313,8 +397,22 @@ namespace Cantera {
/**
* Called by TransportFactory to set parameters.
*/
virtual bool init(TransportParams& tr)
{ err("init"); return false; }
//virtual bool init(TransportParams& tr)
//{ err("init"); return false; }
/**
* Called by TransportFactory to set parameters.
*/
virtual bool initGas( GasTransportParams& tr )
{ err("initGas"); return false; }
/**
* Called by TransportFactory to set parameters.
*/
virtual bool initLiquid( LiquidTransportParams& tr )
{ err("initLiquid"); return false; }
/**

544
Cantera/src/transport/TransportFactory.cpp Executable file → Normal file
View file

@ -20,6 +20,7 @@
#include "MixTransport.h"
#include "SolidTransport.h"
#include "DustyGasTransport.h"
#include "SimpleTransport.h"
#ifdef WITH_IDEAL_SOLUTIONS
#include "LiquidTransport.h"
@ -95,31 +96,21 @@ namespace Cantera {
#endif
};
//////////////////// class Transport methods /////////////////////
void Transport::setThermo(thermo_t& thermo) {
if (!ready()) {
m_thermo = &thermo;
m_nmin = m_thermo->nSpecies();
}
else
throw CanteraError("Transport::setThermo",
"the phase object cannot be changed after "
"the transport manager has been constructed.");
}
void Transport::finalize() {
if (!ready())
m_ready = true;
else
throw CanteraError("Transport::finalize",
"finalize has already been called.");
}
doublereal Transport::err(string msg) const {
throw NotImplemented(msg);
//return 0.0;
}
/**
* getArrhenius() parses the xml element called Arrhenius.
* The Arrhenius expression is
* \f[ k = A T^(b) exp (-E_a / RT). \f]
*/
static void getArrhenius(const XML_Node& node,
doublereal& A, doublereal& b, doublereal& E) {
/* parse the children for the A, b, and E conponents.
*/
A = getFloat(node, "A", "toSI");
b = getFloat(node, "b");
E = getFloat(node, "E", "actEnergy");
E /= GasConstant;
}
//////////////////// class TransportFactory methods //////////////
@ -148,7 +139,7 @@ namespace Cantera {
* @note This method is not used currently.
*/
void TransportFactory::getBinDiffCorrection(doublereal t,
const TransportParams& tr, int k, int j, doublereal xk, doublereal xj,
const GasTransportParams& tr, int k, int j, doublereal xk, doublereal xj,
doublereal& fkj, doublereal& fjk) {
doublereal w1, w2, wsum, sig1, sig2, sig12, sigratio, sigratio2,
@ -220,7 +211,7 @@ namespace Cantera {
* correction, see Dixon-Lewis, Proc. Royal Society (1968).
*/
void TransportFactory::makePolarCorrections(int i, int j,
const TransportParams& tr, doublereal& f_eps, doublereal& f_sigma) {
const GasTransportParams& tr, doublereal& f_eps, doublereal& f_sigma) {
// no correction if both are nonpolar, or both are polar
if (tr.polar[i] == tr.polar[j]) {
@ -262,6 +253,9 @@ namespace Cantera {
m_models["DustyGas"] = cDustyGasTransport;
m_models["CK_Multi"] = CK_Multicomponent;
m_models["CK_Mix"] = CK_MixtureAveraged;
m_models["Liquid"] = cLiquidTransport;
m_models["Aqueous"] = cAqueousTransport;
m_models["Simple"] = cSimpleTransport;
m_models["User"] = cUserTransport;
m_models["None"] = None;
//m_models["Radiative"] = cRadiative;
@ -301,7 +295,7 @@ namespace Cantera {
* make one of several transport models, and return a base class
* pointer to it.
*/
Transport* TransportFactory::newTransport(string transportModel,
Transport* TransportFactory::newTransport(std::string transportModel,
thermo_t* phase, int log_level) {
if (transportModel == "") return new Transport;
@ -342,6 +336,11 @@ namespace Cantera {
dtr = (DustyGasTransport*)tr;
dtr->initialize(phase, gastr);
break;
case cSimpleTransport:
tr = new SimpleTransport();
initLiquidTransport(tr, phase, log_level);
tr->setThermo(*phase);
break;
#ifdef WITH_IDEAL_SOLUTIONS
case cLiquidTransport:
tr = new LiquidTransport;
@ -357,27 +356,48 @@ namespace Cantera {
break;
#endif
default:
throw CanteraError("newTransport","unknown transport model");
throw CanteraError("newTransport","unknown transport model: " + transportModel);
}
phase->restoreState(state);
return tr;
}
/**
* make one of several transport models, and return a base class
* pointer to it.
*/
Transport* TransportFactory::newTransport(thermo_t* phase, int log_level) {
XML_Node &phaseNode=phase->xml();
/*
* Find the Thermo XML node
*/
if (!phaseNode.hasChild("transport")) {
throw CanteraError("TransportFactory::newTransport",
"no transport XML node");
}
XML_Node& transportNode = phaseNode.child("transport");
string transportModel = transportNode.attrib("model");
if (transportModel == "") {
throw CanteraError("TransportFactory::newTransport",
"transport XML node doesn't have a model string");
}
return newTransport(transportModel, phase,log_level);
}
/**
/**
* Prepare to build a new kinetic-theory-based transport manager
* for low-density gases. Uses polynomial fits to Monchick & Mason
* collision integrals.
*/
void TransportFactory::setupMM(std::ostream &flog,
const std::vector<const XML_Node*> &transport_database,
thermo_t* thermo, int mode, int log_level, TransportParams& tr) {
thermo_t* thermo, int mode, int log_level, GasTransportParams& tr) {
// constant mixture attributes
tr.thermo = thermo;
tr.nsp = tr.thermo->nSpecies();
int nsp = tr.nsp;
tr.nsp_ = tr.thermo->nSpecies();
int nsp = tr.nsp_;
tr.tmin = thermo->minTemp();
tr.tmax = thermo->maxTemp();
@ -387,7 +407,7 @@ namespace Cantera {
copy(tr.thermo->molecularWeights().begin(),
tr.thermo->molecularWeights().end(), tr.mw.begin());
tr.mode = mode;
tr.mode_ = mode;
tr.epsilon.resize(nsp, nsp, 0.0);
tr.delta.resize(nsp, nsp, 0.0);
tr.reducedMass.resize(nsp, nsp, 0.0);
@ -456,7 +476,6 @@ namespace Cantera {
}
}
// Chemkin fits the entire T* range in the Monchick and Mason tables,
// so modify tstar_min and tstar_max if in Chemkin compatibility mode
@ -496,30 +515,155 @@ namespace Cantera {
}
/**
* Prepare to build a new transport manager for liquids assuming that
* viscosity transport data is provided in Arhennius form.
*/
void TransportFactory::setupLiquidTransport(std::ostream &flog,
const std::vector<const XML_Node*> &transport_database,
thermo_t* thermo, int log_level, LiquidTransportParams& trParam) {
// constant mixture attributes
trParam.thermo = thermo;
trParam.nsp_ = trParam.thermo->nSpecies();
int nsp = trParam.nsp_;
trParam.tmin = thermo->minTemp();
trParam.tmax = thermo->maxTemp();
trParam.mw.resize(nsp);
trParam.log_level = log_level;
copy(trParam.thermo->molecularWeights().begin(),
trParam.thermo->molecularWeights().end(), trParam.mw.begin());
//trParam.epsilon.resize(nsp, nsp, 0.0);
//trParam.delta.resize(nsp, nsp, 0.0);
//trParam.reducedMass.resize(nsp, nsp, 0.0);
//trParam.dipole.resize(nsp, nsp, 0.0);
//trParam.diam.resize(nsp, nsp, 0.0);
//trParam.polar.resize(nsp, false);
//trParam.poly.resize(nsp);
//trParam.sigma.resize(nsp);
//trParam.eps.resize(nsp);
XML_Node root, log;
getLiquidTransportData(transport_database, log,
trParam.thermo->speciesNames(), trParam);
//int i, j;
//for (i = 0; i < nsp; i++) trParam.poly[i].resize(nsp);
//doublereal ts1, ts2, tstar_min = 1.e8, tstar_max = 0.0;
//doublereal f_eps, f_sigma;
//DenseMatrix& diam = trParam.diam;
//DenseMatrix& epsilon = trParam.epsilon;
//for (i = 0; i < nsp; i++)
// {
// for (j = i; j < nsp; j++)
// {
// // the reduced mass
// trParam.reducedMass(i,j) =
// trParam.mw[i] * trParam.mw[j] / (Avogadro * (trParam.mw[i] + trParam.mw[j]));
//
// // hard-sphere diameter for (i,j) collisions
// diam(i,j) = 0.5*(trParam.sigma[i] + trParam.sigma[j]);
//
// // the effective well depth for (i,j) collisions
// epsilon(i,j) = sqrt(trParam.eps[i]*trParam.eps[j]);
//
// // The polynomial fits of collision integrals vs. T*
// // will be done for the T* from tstar_min to tstar_max
// ts1 = Boltzmann * trParam.tmin/epsilon(i,j);
// ts2 = Boltzmann * trParam.tmax/epsilon(i,j);
// if (ts1 < tstar_min) tstar_min = ts1;
// if (ts2 > tstar_max) tstar_max = ts2;
//
// // the effective dipole moment for (i,j) collisions
// trParam.dipole(i,j) = sqrt(trParam.dipole(i,i)*trParam.dipole(j,j));
//
// // reduced dipole moment delta* (nondimensional)
// doublereal d = diam(i,j);
// trParam.delta(i,j) = 0.5 * trParam.dipole(i,j)*trParam.dipole(i,j)
// / (epsilon(i,j) * d * d * d);
//
// makePolarCorrections(i, j, trParam, f_eps, f_sigma);
// trParam.diam(i,j) *= f_sigma;
// epsilon(i,j) *= f_eps;
//
// // properties are symmetric
// trParam.reducedMass(j,i) = trParam.reducedMass(i,j);
// diam(j,i) = diam(i,j);
// epsilon(j,i) = epsilon(i,j);
// trParam.dipole(j,i) = trParam.dipole(i,j);
// trParam.delta(j,i) = trParam.delta(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 (mode == CK_Mode) {
// tstar_min = 0.101;
// tstar_max = 99.9;
//}
// initialize the collision integral calculator for the desired
// T* range
//#ifdef DEBUG_MODE
// if (m_verbose) {
// trParam.xml->XML_open(flog, "collision_integrals");
// }
//#endif
// m_integrals = new MMCollisionInt;
// m_integrals->init(trParam.xml, tstar_min, tstar_max, log_level);
// fitCollisionIntegrals(flog, trParam);
//#ifdef DEBUG_MODE
// if (m_verbose) {
// trParam.xml->XML_close(flog, "collision_integrals");
// }
//#endif
// // make polynomial fits
//#ifdef DEBUG_MODE
// if (m_verbose) {
// trParam.xml->XML_open(flog, "property fits");
// }
//#endif
// fitProperties(trParam, flog);
//#ifdef DEBUG_MODE
// if (m_verbose) {
// trParam.xml->XML_close(flog, "property fits");
// }
//#endif
}
void TransportFactory::initTransport(Transport* tran,
thermo_t* thermo, int mode, int log_level) {
const std::vector<const XML_Node*> & transport_database = thermo->speciesData();
TransportParams tr;
GasTransportParams trParam;
#ifdef DEBUG_MODE
ofstream flog("transport_log.xml");
tr.xml = new XML_Writer(flog);
trParam.xml = new XML_Writer(flog);
if (m_verbose) {
tr.xml->XML_open(flog, "transport");
trParam.xml->XML_open(flog, "transport");
}
#else
// create the object, but don't associate it with a file
std::ostream &flog(std::cout);
#endif
// set up Monchick and Mason collision integrals
setupMM(flog, transport_database, thermo, mode, log_level, tr);
setupMM(flog, transport_database, thermo, mode, log_level, trParam);
// do model-specific initialization
tran->init(tr);
tran->initGas(trParam);
#ifdef DEBUG_MODE
if (m_verbose) {
tr.xml->XML_close(flog, "transport");
trParam.xml->XML_close(flog, "transport");
}
// finished with log file
flog.close();
@ -528,11 +672,37 @@ namespace Cantera {
}
void
TransportFactory::initLiquidTransport(Transport* tran,
thermo_t* thermo,
int log_level) {
/** Similar to initTransport except uses LiquidTransportParams
* class and calls setupLiquidTransport().
*/
void TransportFactory::initLiquidTransport(Transport* tran,
thermo_t* thermo,
int log_level) {
const std::vector<const XML_Node*> & transport_database = thermo->speciesData();
LiquidTransportParams trParam;
#ifdef DEBUG_MODE
ofstream flog("transport_log.xml");
trParam.xml = new XML_Writer(flog);
if (m_verbose) {
trParam.xml->XML_open(flog, "transport");
}
#else
// create the object, but don't associate it with a file
std::ostream &flog(std::cout);
#endif
setupLiquidTransport(flog, transport_database, thermo, log_level, trParam);
// do model-specific initialization
tran->initLiquid(trParam);
#ifdef DEBUG_MODE
if (m_verbose) {
trParam.xml->XML_close(flog, "transport");
}
// finished with log file
flog.close();
#endif
return;
}
@ -546,12 +716,12 @@ namespace Cantera {
void TransportFactory::fitCollisionIntegrals(ostream& logfile,
TransportParams& tr) {
GasTransportParams& tr) {
vector_fp::iterator dptr;
doublereal dstar;
int nsp = tr.nsp;
int mode = tr.mode;
int nsp = tr.nsp_;
int mode = tr.mode_;
int i, j;
// Chemkin fits to sixth order polynomials
@ -630,7 +800,7 @@ namespace Cantera {
* these species read from the file.
*/
void TransportFactory::getTransportData(const std::vector<const XML_Node*> &xspecies,
XML_Node& log, const std::vector<std::string> &names, TransportParams& tr)
XML_Node& log, const std::vector<std::string> &names, GasTransportParams& tr)
{
string name;
int geom;
@ -653,6 +823,7 @@ namespace Cantera {
for (i = 0; i < nsp; i++) {
const XML_Node& sp = *xspecies[i];
name = sp["name"];
// std::cout << "Processing node for " << name << std::endl;
// put in a try block so that species with no 'transport'
// child are skipped, instead of throwing an exception.
@ -698,7 +869,7 @@ namespace Cantera {
}
}
for (i = 0; i < tr.nsp; i++) {
for (i = 0; i < tr.nsp_; i++) {
GasTransportData& trdat = datatable[names[i]];
@ -743,6 +914,267 @@ namespace Cantera {
}
}
/**
* Read transport property data from a file for a list of species.
* Given the name of a file containing transport property
* parameters and a list of species names, this method returns an
* instance of TransportParams containing the transport data for
* these species read from the file.
*/
void TransportFactory::getLiquidTransportData( const std::vector<const XML_Node*> &xspecies,
XML_Node& log,
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;
doublereal A_visc, n_visc, Tact_visc, hydrodynamic_radius;
doublereal A_thcond, n_thcond, Tact_thcond;
doublereal A_spdiff, n_spdiff, Tact_spdiff;
int nsp = static_cast<int>(xspecies.size());
std::cout << "Size of xspecies " << nsp << std::endl;
// read all entries in database into 'datatable' and check for
// errors. Note that this procedure validates all entries, not
// only those for the species listed in 'names'.
int linenum = 0;
int i;
for (i = 0; i < nsp; i++) {
const XML_Node& sp = *xspecies[i];
name = sp["name"];
vector_fp vCoeff;
// std::cout << "Processing node for " << name << std::endl;
// put in a try block so that species with no 'transport'
// child are skipped, instead of throwing an exception.
try {
if (sp.hasChild("transport")) {
XML_Node& trNode = sp.child("transport");
// Fill datatable with LiquidTransportData objects for error checking
// and then insertion into LiquidTransportData objects below.
LiquidTransportData data;
data.speciesName = name;
/*
* hydrodynamic radius
*
* format:
* <hydrodynamic_radius model="Constant"> 3.0 </hydrodynamic_radius>
* <hydrodynamic_radius> 3.0 </hydrodynamic_radius>
*/
if (trNode.hasChild("hydrodynamic_radius")) {
XML_Node& hnode = trNode.child("hydrodynamic_radius");
std::string model = lowercase(hnode["model"]);
if (model == "" || model == "constant") {
hydrodynamic_radius = hnode.fp_value();
if (hydrodynamic_radius > 0.0) data.hydroradius = hydrodynamic_radius;
else throw TransportDBError(linenum,
"negative or zero hydrodynamic radius");
data.model_hydroradius = LTR_MODEL_CONSTANT;
} else {
throw CanteraError(" TransportFactory::getLiquidTransportData",
"Unknown model for hydrodynamic_radius:" + model);
}
}
/*
* viscosity
*
* format:
* <viscosity model="Constant"> 3.0 </viscosity>
* <viscosity> 3.0 </viscosity>
* <viscosity model="Arrhenius">
* <A units="Pa S"> 1.0 </A>
* <b> 2.0 </b>
* <E units="kcal/gmol"> 3.0 </E>
* </viscosity>
*
* <viscosity model="Coeff">
* <float_array> 0.0. 1.0, 2.0, 3.0, 4.0 </float_array>
* </viscosity>
*
*/
if (trNode.hasChild("viscosity")) {
XML_Node& vnode = trNode.child("viscosity");
std::string model = lowercase(vnode["model"]);
if (model == "" || model == "constant") {
A_visc = vnode.fp_value();
if (A_visc > 0.0) (data.viscCoeffs).push_back(A_visc);
else throw TransportDBError(linenum,
"negative or zero viscosity");
data.model_viscosity = LTR_MODEL_CONSTANT;
} else if (model == "arrhenius") {
getArrhenius(vnode, A_visc, n_visc, Tact_visc);
if (A_visc <= 0.0) {
throw TransportDBError(linenum, "negative or zero viscosity");
}
(data.viscCoeffs).push_back(A_visc);
(data.viscCoeffs).push_back(n_visc);
(data.viscCoeffs).push_back(Tact_visc);
data.model_viscosity = LTR_MODEL_ARRHENIUS;
} else if (model == "coeff") {
getFloatArray(vnode, vCoeff, true);
data.viscCoeffs = vCoeff;
vCoeff.clear();
data.model_viscosity = LTR_MODEL_COEFF;
} else {
throw CanteraError(" TransportFactory::getLiquidTransportData",
"Unknown model for viscosity:" + vnode["model"]);
}
}
/*
* thermal_conductivity
*
* format:
* <thermal_conductivity model="Constant"> 3.0 </thermal_conductivity>
* <thermal_conductivity> 3.0 </thermal_conductivity>
* <thermal_conductivity model="Arrhenius">
* <A units="Pa S"> 1.0 </A>
* <b> 2.0 </b>
* <E units="kcal/gmol"> 3.0 </E>
* </thermal_conductivity>
*
* <thermal_conductivity model="Coeff">
* <float_array> 0.0. 1.0, 2.0, 3.0, 4.0 </float_array>
* </thermal_conductivity>
*
*/
if (trNode.hasChild("thermal_conductivity")) {
XML_Node& tnode = trNode.child("thermal_conductivity");
std::string model = lowercase(tnode["model"]);
if (model == "" || model == "constant") {
A_thcond = tnode.fp_value();
if (A_thcond > 0.0) (data.thermalCondCoeffs).push_back(A_thcond);
else throw TransportDBError(linenum,
"negative or zero thermal_conductivity");
data.model_thermalCond = LTR_MODEL_CONSTANT;
} else if (model == "arrhenius") {
getArrhenius(tnode, A_thcond, n_thcond, Tact_thcond);
if (A_thcond <= 0.0) {
throw TransportDBError(linenum, "negative or zero thermal_conductivity");
}
(data.thermalCondCoeffs).push_back(A_thcond);
(data.thermalCondCoeffs).push_back(n_thcond);
(data.thermalCondCoeffs).push_back(Tact_thcond);
data.model_thermalCond = LTR_MODEL_ARRHENIUS;
} else if (model == "coeff") {
getFloatArray(tnode, vCoeff, true);
data.thermalCondCoeffs = vCoeff;
vCoeff.clear();
data.model_thermalCond = LTR_MODEL_COEFF;
} else {
throw CanteraError(" TransportFactory::getLiquidTransportData",
"Unknown model for thermal_conductivity:" + tnode["model"]);
}
}
/*
* speciesDiffusivity
*
* format:
* <speciesDiffusivity model="Constant"> 3.0 </speciesDiffusivity>
* <speciesDiffusivity> 3.0 </speciesDiffusivity>
* <speciesDiffusivity model="Arrhenius">
* <A units="Pa S"> 1.0 </A>
* <b> 2.0 </b>
* <E units="kcal/gmol"> 3.0 </E>
* </speciesDiffusivity>
*
* <speciesDiffusivity model="Coeff">
* <float_array> 0.0. 1.0, 2.0, 3.0, 4.0 </float_array>
* </speciesDiffusivity>
*
*/
if (trNode.hasChild("speciesDiffusivity")) {
XML_Node& dnode = trNode.child("speciesDiffusivity");
std::string model = lowercase(dnode["model"]);
if (model == "" || model == "constant") {
A_spdiff = dnode.fp_value();
if (A_spdiff > 0.0) (data.speciesDiffusivityCoeffs).push_back(A_spdiff);
else throw TransportDBError(linenum,
"negative or zero speciesDiffusivity");
data.model_speciesDiffusivity = LTR_MODEL_CONSTANT;
} else if (model == "arrhenius") {
getArrhenius(dnode, A_spdiff, n_spdiff, Tact_spdiff);
if (A_spdiff <= 0.0) {
throw TransportDBError(linenum, "negative or zero speciesDiffusivity");
}
(data.speciesDiffusivityCoeffs).push_back(A_spdiff);
(data.speciesDiffusivityCoeffs).push_back(n_spdiff);
(data.speciesDiffusivityCoeffs).push_back(Tact_spdiff);
data.model_speciesDiffusivity = LTR_MODEL_ARRHENIUS;
} else if (model == "coeff") {
getFloatArray(dnode, vCoeff, true);
data.speciesDiffusivityCoeffs = vCoeff;
data.model_speciesDiffusivity = LTR_MODEL_COEFF;
} else {
throw CanteraError(" TransportFactory::getLiquidTransportData",
"Unknown model for speciesDiffusivity:" + dnode["model"]);
}
}
datatable[name] = data;
}
}
catch(CanteraError) {
;
}
}
trParam.LTData.clear();
for (i = 0; i < trParam.nsp_; i++) {
LiquidTransportData& trdat = datatable[names[i]];
// 'datatable' returns a default TransportData object if
// the species name is not one in the transport database.
// This can be detected by examining 'geometry'.
if (trdat.viscCoeffs[0] < 0) {
throw TransportDBError(0,"no transport data found for species "
+ names[i]);
}
// parameters should be converted to SI units before storing
if (trdat.viscCoeffs.size() > 0) {
trParam.visc_A[i] = trdat.viscCoeffs[0] ;
}
if (trdat.viscCoeffs.size() > 2) {
trParam.visc_n[i] = trdat.viscCoeffs[1] ;
trParam.visc_Tact[i] = trdat.viscCoeffs[2] ;
}
if (trdat.thermalCondCoeffs.size() > 0) {
trParam.thermCond_A[i] = trdat.thermalCondCoeffs[0] ;
}
if (trdat.thermalCondCoeffs.size() > 2) {
trParam.thermCond_n[i] = trdat.thermalCondCoeffs[1] ;
trParam.thermCond_Tact[i] = trdat.thermalCondCoeffs[2] ;
}
// Angstroms -> meters
trParam.hydroRadius[i] = 1.e-10 * trdat.hydroradius;
/*
* this is a much more general way to handle the transfer
* -> calling the default copy constructor for LiquidTransportData
*/
trParam.LTData.push_back(trdat);
}
// Need to identify a method to obtain interaction matrices.
// This will fill LiquidTransportParams members visc_Eij, visc_Sij
trParam.visc_Eij.resize(trParam.nsp_,trParam.nsp_);
cout << "No support for species viscosity interactions in TransportFactory.cpp" << endl;
}
/*********************************************************
*
@ -772,7 +1204,7 @@ namespace Cantera {
* D(i,j)/sqrt(k_BT)) = \sum_{n = 0}^4 a_n(i,j) (\log T)^n
* \f]
*/
void TransportFactory::fitProperties(TransportParams& tr,
void TransportFactory::fitProperties(GasTransportParams& tr,
ostream& logfile) {
doublereal tstar;
int k, j, n, ndeg = 0;
@ -782,7 +1214,7 @@ namespace Cantera {
// number of points to use in generating fit data
const int np = 50;
int mode = tr.mode;
int mode = tr.mode_;
int degree = (mode == CK_Mode ? 3 : 4);
doublereal t, om22;
@ -834,7 +1266,7 @@ namespace Cantera {
c1, cv_rot, cv_int, f_rot, f_trans, om11;
doublereal diffcoeff;
for (k = 0; k < tr.nsp; k++)
for (k = 0; k < tr.nsp_; k++)
{
for (n = 0; n < np; n++) {
t = tr.tmin + dt*n;
@ -971,7 +1403,7 @@ namespace Cantera {
tr.xml->XML_comment(logfile,s);
}
if (tr.log_level >= 2)
for (k = 0; k < tr.nsp; k++) {
for (k = 0; k < tr.nsp_; k++) {
tr.xml->XML_writeVector(logfile, " ", tr.thermo->speciesName(k),
degree+1, DATA_PTR(tr.condcoeffs[k]));
}
@ -999,9 +1431,9 @@ namespace Cantera {
mxerr = 0.0, mxrelerr = 0.0;
vector_fp diff(np + 1);
doublereal eps, sigma;
for (k = 0; k < tr.nsp; k++)
for (k = 0; k < tr.nsp_; k++)
{
for (j = k; j < tr.nsp; j++) {
for (j = k; j < tr.nsp_; j++) {
ipoly = tr.poly[k][j];
for (n = 0; n < np; n++) {

89
Cantera/src/transport/TransportFactory.h Executable file → Normal file
View file

@ -34,6 +34,7 @@
#include "ct_defs.h"
#include "TransportBase.h"
#include "FactoryBase.h"
#include "LiquidTransportData.h"
#if defined(THREAD_SAFE_CANTERA)
#include <boost/thread/mutex.hpp>
@ -61,9 +62,11 @@ namespace Cantera {
doublereal rotRelaxNumber;
};
// forward references
class MMCollisionInt;
class TransportParams;
class GasTransportParams;
class LiquidTransportParams;
class XML_Node;
@ -121,23 +124,35 @@ namespace Cantera {
* single instance.
*/
virtual ~TransportFactory();
/// Build a new transport manager
//! Build a new transport manager using a transport manager
//! that may not be the same as in the phase description
/*!
* @param model String name for the transport manager
* @param thermo ThermoPhase object
* @param log_level log level
*/
virtual Transport*
newTransport(std::string model="", thermo_t* thermo=0, int log_level=0);
newTransport(std::string model, thermo_t* thermo, int log_level=0);
//! Build a new transport manager using the default transport manager
//! in the phase description
/*!
* @param thermo ThermoPhase object
* @param log_level log level
*/
virtual Transport*
newTransport(thermo_t* thermo, int log_level=0);
/// Initialize an existing transport manager
virtual void initTransport(Transport* tr,
thermo_t* thermo=0, int mode=0, int log_level=0);
thermo_t* thermo, int mode=0, int log_level=0);
/// Initialize an existing transport manager for liquid phase
virtual void initLiquidTransport(Transport* tr,
thermo_t* thermo=0,
thermo_t* thermo,
int log_level=0);
private:
//! Static instance of the factor -> This is the only instance of this
@ -147,8 +162,6 @@ namespace Cantera {
static boost::mutex transport_mutex ;
#endif
//! The constructor is private; use static method factory() to
//! get a pointer to a factory instance
/*!
@ -161,34 +174,54 @@ namespace Cantera {
void getTransportData(const std::vector<const XML_Node*> &db,
XML_Node& log, const std::vector<std::string>& names,
TransportParams& tr);
GasTransportParams& tr);
//! Read transport property data from a file for a list of species.
/*!
*
* Given the name of a file containing transport property
* parameters and a list of species names, this method returns an
* instance of TransportParams containing the transport data for
* these species read from the file.
*
*/
void getLiquidTransportData(const std::vector<const XML_Node*> &db,
XML_Node& log, const std::vector<std::string>& names,
LiquidTransportParams& tr);
/** Generate polynomial fits to viscosity, conductivity, and
* binary diffusion coefficients */
void fitProperties(TransportParams& tr, std::ostream & logfile);
void fitProperties(GasTransportParams& tr, std::ostream & logfile);
/// Generate polynomial fits to collision integrals
void fitCollisionIntegrals(std::ostream & logfile,
TransportParams& tr);
GasTransportParams& tr);
void setupMM(std::ostream &flog, const std::vector<const XML_Node*> &transport_database,
thermo_t* thermo, int mode, int log_level,
TransportParams& tr);
GasTransportParams& tr);
void setupLiquidTransport(std::ostream &flog, const std::vector<const XML_Node*> &transport_database,
thermo_t* thermo, int log_level,
LiquidTransportParams& tr);
/// Second-order correction to the binary diffusion coefficients
void getBinDiffCorrection(doublereal t,
const TransportParams& tr, int k, int j,
const GasTransportParams& tr, int k, int j,
doublereal xk, doublereal xj,
doublereal& fkj, doublereal& fjk);
/// Corrections for polar-nonpolar binary diffusion coefficients
void makePolarCorrections(int i, int j,
const TransportParams& tr, doublereal& f_eps,
const GasTransportParams& tr, doublereal& f_eps,
doublereal& f_sigma);
//! Boolean indicating whether to turn on verbose printing
bool m_verbose;
@ -205,8 +238,8 @@ namespace Cantera {
* Create a new transport manager instance.
* @ingroup transportProps
*/
inline Transport* newTransportMgr(std::string transportModel="",
thermo_t* thermo=0, int loglevel=0,
inline Transport* newTransportMgr(std::string transportModel = "",
thermo_t* thermo = 0, int loglevel=0,
TransportFactory* f=0) {
if (f == 0) {
f = TransportFactory::factory();
@ -218,7 +251,25 @@ namespace Cantera {
* the need for multiple cantera and transport library statements
* for applications that don't have transport in them.
*/
//TransportFactory::deleteFactory();
return ptr;
}
/**
* Create a new transport manager instance.
* @ingroup transportProps
*/
inline Transport* newDefaultTransportMgr(thermo_t* thermo, int loglevel=0,
TransportFactory* f=0) {
if (f == 0) {
f = TransportFactory::factory();
}
Transport* ptr = f->newTransport(thermo, loglevel);
/*
* Note: We delete the static s_factory instance here, instead of in
* appdelete() in misc.cpp, to avoid linking problems involving
* the need for multiple cantera and transport library statements
* for applications that don't have transport in them.
*/
return ptr;
}

53
Cantera/src/transport/TransportParams.h
View file

@ -3,7 +3,6 @@
#include <vector>
#include "ct_defs.h"
#include "TransportBase.h"
#include "xml.h"
@ -11,10 +10,9 @@
namespace Cantera {
/**
*
* Holds transport data. Used by TransportFactory.
*
/**
* Base class to hold transport model parameters.
* Used by TransportFactory.
*/
class TransportParams {
@ -22,18 +20,43 @@ namespace Cantera {
TransportParams() : thermo(0), xml(0) {}
virtual ~TransportParams();
int nsp;
int nsp_;
// phase_t* mix;
thermo_t* thermo;
vector_fp mw;
// polynomial fits
std::vector<vector_fp> visccoeffs;
std::vector<vector_fp> condcoeffs;
std::vector<vector_fp> diffcoeffs;
vector_fp polytempvec;
//temperature-fit viscosity
std::vector<vector_fp> visccoeffs;
//temperature-fit heat conduction
std::vector<vector_fp> condcoeffs;
//temperature-fit diffusivity
std::vector<vector_fp> diffcoeffs;
vector_fp polytempvec;
//minimum and maximum temperatures for parameter fits
doublereal tmax, tmin;
int mode_;
XML_Writer* xml;
int log_level;
};
/**
* Holds transport model parameters relevant to transport in ideal
* gases with a kinetic theory of gases derived transport model.
* Used by TransportFactory.
*/
class GasTransportParams : public TransportParams {
public:
GasTransportParams() {}
~GasTransportParams() {}
std::vector<std::vector<int> > poly;
std::vector<vector_fp > omega22_poly;
std::vector<vector_fp > astar_poly;
@ -53,11 +76,9 @@ namespace Cantera {
DenseMatrix epsilon;
DenseMatrix dipole;
DenseMatrix delta;
doublereal tmax, tmin;
int mode;
XML_Writer* xml;
int log_level;
};
}
#endif
#endif //CT_TRANSPORTPARAMS_H

2
Cantera/src/transport/WaterTransport.h
View file

@ -81,7 +81,7 @@ namespace Cantera {
virtual ~WaterTransport();
//! Return the model id for this transport parameterization
virtual int model() {
virtual int model() const {
return cWaterTransport;
}

3
configure vendored
View file

File diff suppressed because one or more lines are too long