From 76526a6aa3b7d4455bf1dfeabe5cf13d2af758ce Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Fri, 4 Jul 2003 06:37:42 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/clib/src/ct.cpp | 4 + Cantera/clib/src/ct.h | 1 + .../matlab/cantera/@Transport/setParameters.m | 5 + .../@Transport/setThermalConductivity.m | 9 + Cantera/matlab/cantera/examples/flame2.m | 4 +- .../cantera/private/transportmethods.cpp | 20 + Cantera/python/Cantera/OneDim.py | 1 - Cantera/python/Cantera/SolidTransport.py | 34 + Cantera/python/Cantera/Transport.py | 8 +- Cantera/python/Cantera/solid.py | 46 ++ Cantera/python/Cantera/solids.py | 2 +- Cantera/python/src/cttransport_methods.cpp | 15 + Cantera/python/src/methods.h | 1 + Cantera/src/ImplicitSurfChem.cpp | 15 +- Cantera/src/ImplicitSurfChem.h | 9 +- Cantera/src/InterfaceKinetics.cpp | 63 +- Cantera/src/InterfaceKinetics.h | 13 +- Cantera/src/Kinetics.h | 12 + Cantera/src/Makefile.in | 4 +- Cantera/src/MetalPhase.h | 84 +++ Cantera/src/RateCoeffMgr.h | 4 + Cantera/src/RxnRates.h | 4 + Cantera/src/SolidCompound.cpp | 50 ++ Cantera/src/SolidCompound.h | 187 ++++++ Cantera/src/SurfPhase.h | 2 +- Cantera/src/ThermoFactory.cpp | 35 +- Cantera/src/importCTML.cpp | 53 +- Cantera/src/mix_defs.h | 2 + Cantera/src/oneD/Domain1D.h | 8 +- Cantera/src/oneD/MultiNewton.cpp | 15 +- Cantera/src/oneD/Sim1D.cpp | 6 +- Cantera/src/oneD/Solid1D.cpp | 634 ++++++++++++++++++ Cantera/src/oneD/Solid1D.h | 402 +++++++++++ Cantera/src/oneD/refine.cpp | 10 +- Cantera/src/transport/SolidTransport.cpp | 84 +++ Cantera/src/transport/SolidTransport.h | 81 +++ Cantera/src/units.h | 14 +- data/inputs/elements.xml | 193 +++--- examples/cxx/flame1.cpp | 60 ++ include/Interface.h | 53 ++ include/onedim.h | 2 +- 41 files changed, 2086 insertions(+), 163 deletions(-) create mode 100644 Cantera/matlab/cantera/@Transport/setParameters.m create mode 100644 Cantera/matlab/cantera/@Transport/setThermalConductivity.m create mode 100644 Cantera/python/Cantera/SolidTransport.py create mode 100644 Cantera/python/Cantera/solid.py create mode 100644 Cantera/src/MetalPhase.h create mode 100644 Cantera/src/SolidCompound.cpp create mode 100644 Cantera/src/SolidCompound.h create mode 100644 Cantera/src/oneD/Solid1D.cpp create mode 100644 Cantera/src/oneD/Solid1D.h create mode 100644 Cantera/src/transport/SolidTransport.cpp create mode 100644 Cantera/src/transport/SolidTransport.h create mode 100644 examples/cxx/flame1.cpp create mode 100644 include/Interface.h diff --git a/Cantera/clib/src/ct.cpp b/Cantera/clib/src/ct.cpp index 28fbd5041..7b4f38ca4 100755 --- a/Cantera/clib/src/ct.cpp +++ b/Cantera/clib/src/ct.cpp @@ -728,6 +728,10 @@ extern "C" { catch (CanteraError) { return -1; } } + int DLL_EXPORT trans_setParameters(int n, int type, int k, double* d) { + try { trans(n)->setParameters(type, k, d); return 0;} + catch (CanteraError) { return -1; } + } //-------------------- Functions --------------------------- diff --git a/Cantera/clib/src/ct.h b/Cantera/clib/src/ct.h index 7b6b5ce8b..aa76f3a92 100755 --- a/Cantera/clib/src/ct.h +++ b/Cantera/clib/src/ct.h @@ -108,6 +108,7 @@ extern "C" { int DLL_IMPORT trans_getMixDiffCoeffs(int n, int ld, double* d); int DLL_IMPORT trans_getBinDiffCoeffs(int n, int ld, double* d); int DLL_IMPORT trans_getMultiDiffCoeffs(int n, int ld, double* d); + int DLL_IMPORT trans_setParameters(int n, int type, int k, double* d); int DLL_IMPORT import_phase(int nth, int nxml, char* id); int DLL_IMPORT import_kinetics(int nxml, char* id, diff --git a/Cantera/matlab/cantera/@Transport/setParameters.m b/Cantera/matlab/cantera/@Transport/setParameters.m new file mode 100644 index 000000000..8e254aae7 --- /dev/null +++ b/Cantera/matlab/cantera/@Transport/setParameters.m @@ -0,0 +1,5 @@ +function setParameters(tr, type, k, p) +% SETPARAMETERS - set parameters. +% +v = trans_get(tr.id, 31, type, k, p); + diff --git a/Cantera/matlab/cantera/@Transport/setThermalConductivity.m b/Cantera/matlab/cantera/@Transport/setThermalConductivity.m new file mode 100644 index 000000000..ce698b79a --- /dev/null +++ b/Cantera/matlab/cantera/@Transport/setThermalConductivity.m @@ -0,0 +1,9 @@ +function setThermalConductivity(tr, lam) +% SETTHERMALCONDUCTIVITY - Set the thermal conductivity. +% +% This method can only be used with transport models that +% support directly setting the value of the thermal +% conductivity. +% +setParameters(tr, 1, 0, lam); + diff --git a/Cantera/matlab/cantera/examples/flame2.m b/Cantera/matlab/cantera/examples/flame2.m index 998a2086d..dd3cb0b79 100644 --- a/Cantera/matlab/cantera/examples/flame2.m +++ b/Cantera/matlab/cantera/examples/flame2.m @@ -10,8 +10,8 @@ t0 = cputime; % record the starting time % parameter values p = oneatm; % pressure tin = 300.0; % inlet temperature -mdot_o = 0.24; % air, kg/m^2/s -mdot_f = 0.08; % fuel, kg/m^2/s +mdot_o = 0.72; % air, kg/m^2/s +mdot_f = 0.24; % fuel, kg/m^2/s rxnmech = 'gri30.xml'; % reaction mechanism file transport = 'Mix'; % transport model diff --git a/Cantera/matlab/cantera/private/transportmethods.cpp b/Cantera/matlab/cantera/private/transportmethods.cpp index 0bc92e891..942429378 100644 --- a/Cantera/matlab/cantera/private/transportmethods.cpp +++ b/Cantera/matlab/cantera/private/transportmethods.cpp @@ -67,6 +67,7 @@ void reportError(); mexErrMsgTxt("unknown Transport method"); } } + else if (job < 30) { nsp = getInt(prhs[3]); plhs[0] = mxCreateNumericMatrix(nsp,nsp,mxDOUBLE_CLASS,mxREAL); @@ -80,6 +81,25 @@ void reportError(); mexErrMsgTxt("unknown Transport method"); } } + + // set parameters + else if (job < 40) { + double* params; + int typ, k; + switch (job) { + case 31: + typ = getInt(prhs[3]); + k = getInt(prhs[4]); + params = mxGetPr(prhs[5]); + iok = trans_setParameters(n, typ, k, params); + break; + default: + mexErrMsgTxt("unknown Transport method"); + } + plhs[0] = mxCreateNumericMatrix(1,1,mxDOUBLE_CLASS,mxREAL); + h = mxGetPr(plhs[0]); + *h = double(iok); + } else { mexErrMsgTxt("unknown Transport method"); } diff --git a/Cantera/python/Cantera/OneDim.py b/Cantera/python/Cantera/OneDim.py index f836d3d73..973d1ce56 100755 --- a/Cantera/python/Cantera/OneDim.py +++ b/Cantera/python/Cantera/OneDim.py @@ -404,7 +404,6 @@ class OneDim: """ self.setNewtonOptions(max_jac_age = self._opt['ts_jac_age']) - print 'max jac age = ',self._opt['ts_jac_age'] if loglevel > 0: print_heading('Begin time integration.\n\n') diff --git a/Cantera/python/Cantera/SolidTransport.py b/Cantera/python/Cantera/SolidTransport.py new file mode 100644 index 000000000..91f18ee37 --- /dev/null +++ b/Cantera/python/Cantera/SolidTransport.py @@ -0,0 +1,34 @@ +""" +Transport properties for solids. + +This class implements a simple model for the diffusion coefficients and +the thermal conductivity of a solid. The diffusion coefficients have +modified Arrhenius form, and the thermal conductivity is constant. +All parameters are user-specified, not computed from a physical model. + +Examples: + + >>> tr = SolidTransport(solid_phase) + >>> tr.setThermalConductivity(0.5) # W/m/K + >>> tr.setDiffCoeff(species = "OxygenIon", A = 2.0, n = 0.0, E = 700.0) + +Note that the diffusion coefficient is computed from D = A * T^n * +exp(-E/t) in m^2/s. Diffusion coefficients for unspecified species are +set to zero. + +""" + +from Cantera.Transport import Transport + +class SolidTransport(Transport): + def __init__(self, phase = None): + Transport.__init__(self, model = "Solid", phase = phase) + + def setThermalConductivity(self, lamb): + self.setParameters(1, 0, [lamb, 0.0]) + + def setDiffCoeff(self, species = "", A = 0.0, n = 0.0, E = 0.0): + k = self._phase.speciesIndex(species) + self.setParameters(0, k, [A, n, E]) + + diff --git a/Cantera/python/Cantera/Transport.py b/Cantera/python/Cantera/Transport.py index 6caabf970..71883f971 100755 --- a/Cantera/python/Cantera/Transport.py +++ b/Cantera/python/Cantera/Transport.py @@ -1,4 +1,5 @@ import _cantera +from Numeric import asarray class Transport: @@ -80,5 +81,10 @@ class Transport: def multiDiffCoeffs(self): return _cantera.tran_multiDiffCoeffs(self.__tr_id, - self.trnsp) + self.trnsp) + + def setParameters(self, type, k, params): + return _cantera.tran_setParameters(self.__tr_id, + type, k, asarray(params)) + diff --git a/Cantera/python/Cantera/solid.py b/Cantera/python/Cantera/solid.py new file mode 100644 index 000000000..da360e03c --- /dev/null +++ b/Cantera/python/Cantera/solid.py @@ -0,0 +1,46 @@ +""" +""" + +import string +import os + +from constants import * +from ThermoPhase import ThermoPhase +from Kinetics import Kinetics +from SolidTransport import SolidTransport +import XML +import _cantera + +class Solid(ThermoPhase, Kinetics, SolidTransport): + """ + """ + + def __init__(self, src="", root=None): + + self.ckin = 0 + self._owner = 0 + self.verbose = 1 + + # get the 'phase' element + s = XML.find_XML(src=src, root=root, name="phase") + + # get the equation of state model + ThermoPhase.__init__(self, xml_phase=s) + + # get the kinetics model + Kinetics.__init__(self, xml_phase=s, phases=[self]) + + SolidTransport.__init__(self, phase=self) + + #self.setState_TP(300.0, OneAtm) + + + def __repr__(self): + return _cantera.phase_report(self._phase_id, self.verbose) + + + def __del__(self): + SolidTransport.__del__(self) + Kinetics.__del__(self) + ThermoPhase.__del__(self) + diff --git a/Cantera/python/Cantera/solids.py b/Cantera/python/Cantera/solids.py index 9123e57b5..6be3da64d 100755 --- a/Cantera/python/Cantera/solids.py +++ b/Cantera/python/Cantera/solids.py @@ -2,7 +2,7 @@ from solution import Solution -def Solid(import_file="", id="", +def Solid(src="", kmodel=1, transport=None): return Solution(import_file=import_file, thermo_db="", diff --git a/Cantera/python/src/cttransport_methods.cpp b/Cantera/python/src/cttransport_methods.cpp index 00953e9c0..bcffda9aa 100644 --- a/Cantera/python/src/cttransport_methods.cpp +++ b/Cantera/python/src/cttransport_methods.cpp @@ -30,6 +30,21 @@ py_transport_delete(PyObject *self, PyObject *args) } +static PyObject* +py_setParameters(PyObject *self, PyObject *args) { + int n, k, typ; + PyObject* parray; + if (!PyArg_ParseTuple(args, "iiiO:py_setParameters", + &n, &typ, &k, &parray)) return NULL; + + PyArrayObject* a = (PyArrayObject*)parray; + double* xd = (double*)a->data; + int ok = trans_setParameters(n, typ, k, xd); + if (ok < 0) return reportError(ok); + return Py_BuildValue("i",ok); +} + + static PyObject* py_viscosity(PyObject *self, PyObject *args) { int n; diff --git a/Cantera/python/src/methods.h b/Cantera/python/src/methods.h index 2597488a7..b87e0546c 100644 --- a/Cantera/python/src/methods.h +++ b/Cantera/python/src/methods.h @@ -75,6 +75,7 @@ static PyMethodDef ct_methods[] = { {"tran_binaryDiffCoeffs", py_binaryDiffCoeffs, METH_VARARGS}, {"tran_mixDiffCoeffs", py_mixDiffCoeffs, METH_VARARGS}, {"tran_multiDiffCoeffs", py_multiDiffCoeffs, METH_VARARGS}, + {"tran_setParameters", py_setParameters, METH_VARARGS}, {"get_Cantera_Error", ct_get_cantera_error, METH_VARARGS}, {"ct_print", ct_print, METH_VARARGS}, diff --git a/Cantera/src/ImplicitSurfChem.cpp b/Cantera/src/ImplicitSurfChem.cpp index 567d32253..6e91017e1 100755 --- a/Cantera/src/ImplicitSurfChem.cpp +++ b/Cantera/src/ImplicitSurfChem.cpp @@ -21,12 +21,12 @@ namespace Cantera { - ImplicitSurfChem::ImplicitSurfChem(SurfKinetics& kin) + ImplicitSurfChem::ImplicitSurfChem(InterfaceKinetics& kin) : FuncEval(), m_kin(&kin), m_integ(0), - m_atol(1.e-15), m_rtol(1.e-7), m_maxstep(0.0) + m_atol(1.e-14), m_rtol(1.e-7), m_maxstep(0.0) { m_integ = new CVodeInt; - m_surf = &kin.sphase(); + m_surf = (SurfPhase*)&kin.thermo(kin.nPhases()-1); // use backward differencing, with a full Jacobian computed // numerically, and use a Newton linear iterator @@ -35,7 +35,7 @@ namespace Cantera { m_integ->setProblemType(DENSE + NOJAC); m_integ->setIterator(Newton_Iter); m_nsp = m_surf->nSpecies(); - m_work.resize(m_kin->nTotal()); + m_work.resize(m_kin->nTotalSpecies()); } // overloaded method of FuncEval. Called by the integrator to @@ -69,11 +69,14 @@ namespace Cantera { doublereal rs0 = 1.0/m_surf->siteDensity(); m_kin->getNetProductionRates(m_work.begin()); int k; - int kbulk = m_kin->nTotal() - m_nsp; + int kbulk = m_kin->nTotalSpecies() - m_nsp; doublereal sum = 0.0; - for (k = 0; k < m_nsp; k++) { + for (k = 1; k < m_nsp; k++) { ydot[k] = m_work[kbulk + k] * rs0 * m_surf->size(k); + sum -= ydot[k]; } + //if (sum < 0.0) sum = 0.0; + ydot[0] = sum; } } diff --git a/Cantera/src/ImplicitSurfChem.h b/Cantera/src/ImplicitSurfChem.h index d80f976b6..72be5bdf0 100755 --- a/Cantera/src/ImplicitSurfChem.h +++ b/Cantera/src/ImplicitSurfChem.h @@ -21,7 +21,8 @@ #include "FuncEval.h" #include "CVode.h" -#include "surfKinetics.h" +#include "InterfaceKinetics.h" +#include "SurfPhase.h" namespace Cantera { @@ -37,7 +38,7 @@ namespace Cantera { /** * Constructor. */ - ImplicitSurfChem(SurfKinetics& kin); + ImplicitSurfChem(InterfaceKinetics& kin); /** @@ -90,8 +91,8 @@ namespace Cantera { */ void updateState(doublereal* y); - SurfacePhase* m_surf; - SurfKinetics* m_kin; + SurfPhase* m_surf; + InterfaceKinetics* m_kin; int m_nsp; Integrator* m_integ; // pointer to integrator doublereal m_atol, m_rtol; // tolerances diff --git a/Cantera/src/InterfaceKinetics.cpp b/Cantera/src/InterfaceKinetics.cpp index 95c0f16ab..48aa69b6d 100644 --- a/Cantera/src/InterfaceKinetics.cpp +++ b/Cantera/src/InterfaceKinetics.cpp @@ -111,7 +111,7 @@ namespace Cantera { void SurfPhase:: setElectricPotential(doublereal V) { for (int k = 0; k < m_kk; k++) { - m_pe[k] = charge(k)*Faraday; + m_pe[k] = charge(k)*Faraday*V; } _updateThermo(true); } @@ -125,7 +125,7 @@ namespace Cantera { } void SurfPhase:: - getCoverages(doublereal* theta) { + getCoverages(doublereal* theta) const { getConcentrations(theta); for (int k = 0; k < m_kk; k++) { theta[k] *= size(k)/m_n0; @@ -147,7 +147,7 @@ namespace Cantera { m_s0[k] *= GasConstant; m_cp0[k] *= GasConstant; deltaE = m_pe[k]; - m_h0[k] += deltaE; + //m_h0[k] += deltaE; m_mu0[k] = m_h0[k] - tnow*m_s0[k]; } m_tlast = tnow; @@ -166,6 +166,7 @@ namespace Cantera { InterfaceKinetics(thermo_t* thermo) : Kinetics(thermo), m_kk(0), + m_redo_rates(false), m_nirrev(0), m_nrev(0), m_finalized(false) @@ -177,12 +178,14 @@ namespace Cantera { void InterfaceKinetics:: _update_rates_T() { doublereal T = thermo().temperature(); - if (T != m_kdata->m_temp) { + if (T != m_kdata->m_temp || m_redo_rates) { doublereal logT = log(T); m_rates.update(T, logT, m_kdata->m_rfn.begin()); + correctElectronTransferRates(m_kdata->m_rfn.begin()); m_kdata->m_temp = T; updateKc(); m_kdata->m_ROP_ok = false; + m_redo_rates = false; } }; @@ -212,10 +215,14 @@ namespace Cantera { doublereal rrt = 1.0/rt; int np = nPhases(); for (n = 0; n < np; n++) { + // cout << n << "start = " << m_start[n] << endl; thermo(n).getStandardChemPotentials(m_mu0.begin() + m_start[n]); nsp = thermo(n).nSpecies(); for (k = 0; k < nsp; k++) { + //cout << ik << "mu0 = " << m_mu0[ik] << endl; m_mu0[ik] -= rt*thermo(n).logStandardConc(k); + m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k); + //cout << ik << "mu0 = " << m_mu0[ik] << endl; ik++; } } @@ -228,7 +235,9 @@ namespace Cantera { for (i = 0; i < m_nrev; i++) { irxn = m_revindex[i]; + //cout << "rev " << irxn << " " << m_rkc[irxn] << endl; m_rkc[irxn] = exp(m_rkc[irxn]*rrt); + //cout << "rev " << irxn << " " << m_rkc[irxn] << endl; } for(i = 0; i != m_nirrev; ++i) { @@ -251,7 +260,14 @@ namespace Cantera { thermo(n).getStandardChemPotentials(m_mu0.begin() + m_start[n]); nsp = thermo(n).nSpecies(); for (k = 0; k < nsp; k++) { + //cout << thermo(n).id() << " " << thermo(n).speciesName(k) + // << " " << m_mu0[ik] << endl; m_mu0[ik] -= rt*thermo(n).logStandardConc(k); + m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k); + //if (thermo(n).charge(k) != 0.0) { + // cout << thermo(n).id() << " " << thermo(n).speciesName(k) + // << " " << m_phi[n] << " " << thermo(n).charge(k) << endl; + //} ik++; } } @@ -268,6 +284,41 @@ namespace Cantera { } + /** + * Get the equilibrium constants of all reactions, whether + * reversible or not. + */ + void InterfaceKinetics::correctElectronTransferRates(doublereal* kf) { + int i; + + int n, nsp, k, ik=0; + doublereal rt = GasConstant*thermo(0).temperature(); + doublereal rrt = 1.0/rt; + int np = nPhases(); + for (n = 0; n < np; n++) { + nsp = thermo(n).nSpecies(); + for (k = 0; k < nsp; k++) { + m_pot[ik] = Faraday*thermo(n).charge(k)*m_phi[n]; + ik++; + } + } + fill(m_rwork.begin(), m_rwork.begin() + m_ii, 0.0); + m_reactantStoich.decrementReactions(m_pot.begin(), m_rwork.begin()); + m_revProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin()); + m_irrevProductStoich.incrementReactions(m_pot.begin(), m_rwork.begin()); + doublereal eamod, ea; + for (i = 0; i < m_ii; i++) { + //loc = m_index[i].second; + //if (loc >= 0) { + // const Arrhenius& r = m_rates.rateCoeff(m_index[i].second); + // ea = GasConstant*r.activationEnergy_R(); + eamod = 0.5*m_rwork[i]; + if (m_index[i].second >= 0) kf[i] *= exp(-eamod*rrt); + } + } + + + void InterfaceKinetics::updateROP() { _update_rates_T(); @@ -312,7 +363,6 @@ namespace Cantera { m_kdata->m_ROP_ok = true; } - void InterfaceKinetics:: addReaction(const ReactionData& r) { @@ -436,9 +486,12 @@ namespace Cantera { m_prxn.resize(m_kk); m_conc.resize(m_kk); m_mu0.resize(m_kk); + m_pot.resize(m_kk, 0.0); + m_phi.resize(np,0.0); } void InterfaceKinetics::finalize() { + m_rwork.resize(nReactions()); m_finalized = true; } diff --git a/Cantera/src/InterfaceKinetics.h b/Cantera/src/InterfaceKinetics.h index b6bccb617..20d8377d9 100644 --- a/Cantera/src/InterfaceKinetics.h +++ b/Cantera/src/InterfaceKinetics.h @@ -66,6 +66,11 @@ namespace Cantera { virtual int ID() { return cInterfaceKinetics; } + void setElectricPotential(int n, doublereal V) { + m_phi[n] = V; + m_redo_rates = true; + } + virtual doublereal reactantStoichCoeff(int k, int i) const { return m_rrxn[k][i]; } @@ -153,7 +158,7 @@ namespace Cantera { < m_revindex.end()) return true; else return false; } - + void correctElectronTransferRates(doublereal* kf); void _update_rates_T(); void _update_rates_C(); @@ -162,7 +167,8 @@ namespace Cantera { int m_kk; Rate1 m_rates; - + bool m_redo_rates; + mutable map > m_index; vector m_irrev; @@ -191,6 +197,9 @@ namespace Cantera { vector_fp m_conc; vector_fp m_mu0; + vector_fp m_phi; + vector_fp m_pot; + vector_fp m_rwork; private: diff --git a/Cantera/src/Kinetics.h b/Cantera/src/Kinetics.h index 14650ebe2..92d6e9093 100755 --- a/Cantera/src/Kinetics.h +++ b/Cantera/src/Kinetics.h @@ -251,6 +251,18 @@ namespace Cantera { return -2; } + thermo_t& speciesPhase(string nm) { + int np = m_thermo.size(); + int k; + string id; + for (int n = 0; n < np; n++) { + k = thermo(n).speciesIndex(nm); + if (k >= 0) return thermo(n); + } + throw CanteraError("speciesPhase", "unknown species "+nm); + } + + /** * Prepare to add reactions. */ diff --git a/Cantera/src/Makefile.in b/Cantera/src/Makefile.in index 9610e030c..d0cb9c177 100755 --- a/Cantera/src/Makefile.in +++ b/Cantera/src/Makefile.in @@ -26,13 +26,13 @@ BASE = Elements.o Constituents.o stringUtils.o misc.o importCTML.o plots.o xml.o Phase.o DenseMatrix.o ctml.o funcs.o ctvector.o phasereport.o # thermodynamic properties -THERMO = $(BASE) ThermoPhase.o IdealGasPhase.o ConstDensityThermo.o SpeciesThermoFactory.o ThermoFactory.o +THERMO = $(BASE) ThermoPhase.o IdealGasPhase.o ConstDensityThermo.o SolidCompound.o SpeciesThermoFactory.o ThermoFactory.o # homogeneous kinetics KINETICS = GRI_30_Kinetics.o KineticsFactory.o GasKinetics.o FalloffFactory.o GasKineticsWriter.o $(THERMO) # heterogeneous kinetics -HETEROKIN = InterfaceKinetics.o $(THERMO) +HETEROKIN = InterfaceKinetics.o ImplicitSurfChem.o $(THERMO) # support for importing from Chemkin-compatible reaction mechanisms CK = $(KINETICS) diff --git a/Cantera/src/MetalPhase.h b/Cantera/src/MetalPhase.h new file mode 100644 index 000000000..f0a37372e --- /dev/null +++ b/Cantera/src/MetalPhase.h @@ -0,0 +1,84 @@ +/** + * + * @file MetalPhase.h + * + */ + +/* $Author$ + * $Date$ + * $Revision$ + * + * Copyright 2003 California Institute of Technology + * + */ + + +#ifndef CT_METALPHASE_H +#define CT_METALPHASE_H + +//#include "ct_defs.h" +#include "mix_defs.h" +#include "ThermoPhase.h" +#include "SpeciesThermo.h" + +namespace Cantera { + + /** + * @ingroup thermoprops + * + * Class MetalPhase represents electrons in a metal. + * + */ + class MetalPhase : public ThermoPhase { + + public: + + MetalPhase() {} + + virtual ~MetalPhase() {} + + /** + * Equation of state flag. + */ + virtual int eosType() const { return cMetal; } + + virtual doublereal enthalpy_mole() const { return 0.0; } + virtual doublereal intEnergy_mole() const { return 0.0; } + virtual doublereal entropy_mole() const { return 0.0; } + virtual doublereal gibbs_mole() const { return 0.0; } + virtual doublereal cp_mole() const { return 0.0; } + virtual doublereal cv_mole() const { return 0.0; } + + virtual void getChemPotentials(doublereal* mu) const { + mu[0] = 0.0; + } + + virtual void getStandardChemPotentials(doublereal* mu0) const { + mu0[0] = 0.0; + } + + virtual void getActivityConcentrations(doublereal* c) const { + c[0] = 1.0; + } + + virtual doublereal standardConcentration(int k=0) const { + return 1.0; + } + + virtual doublereal logStandardConc(int k=0) const { + return 0.0; + } + + protected: + + private: + + }; +} + +#endif + + + + + diff --git a/Cantera/src/RateCoeffMgr.h b/Cantera/src/RateCoeffMgr.h index 333de48da..3ca127509 100755 --- a/Cantera/src/RateCoeffMgr.h +++ b/Cantera/src/RateCoeffMgr.h @@ -78,6 +78,10 @@ namespace Cantera { return -1; } + const R& rateCoeff(int loc) const { + return m_rates[loc]; + } + void update_C(const doublereal* c) { TYPENAME_KEYWORD vector::iterator b = m_rates.begin(); TYPENAME_KEYWORD vector::iterator e = m_rates.end(); diff --git a/Cantera/src/RxnRates.h b/Cantera/src/RxnRates.h index 2d837d6b5..8ece6bdaf 100755 --- a/Cantera/src/RxnRates.h +++ b/Cantera/src/RxnRates.h @@ -42,6 +42,10 @@ namespace Cantera { s << ");" << endl; } + doublereal activationEnergy_R() const { + return m_E; + } + protected: doublereal m_logA, m_b, m_E; }; diff --git a/Cantera/src/SolidCompound.cpp b/Cantera/src/SolidCompound.cpp new file mode 100644 index 000000000..099a44794 --- /dev/null +++ b/Cantera/src/SolidCompound.cpp @@ -0,0 +1,50 @@ +/** + * + * @file IdealGasPhase.cpp + * + */ + +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "ct_defs.h" +#include "mix_defs.h" +#include "SolidCompound.h" +#include "SpeciesThermo.h" + +namespace Cantera { + + void SolidCompound::initThermo() { + m_kk = nSpecies(); + if (m_kk > 1) { + throw CanteraError("initThermo", + "solid compounds may only contain one species."); + } + doublereal tmin = m_spthermo->minTemp(); + doublereal tmax = m_spthermo->maxTemp(); + if (tmin > 0.0) m_tmin = tmin; + if (tmax > 0.0) m_tmax = tmax; + m_p0 = refPressure(); + + int leng = m_kk; + m_h0_RT.resize(leng); + m_cp0_R.resize(leng); + m_s0_R.resize(leng); + } + + + void SolidCompound::_updateThermo() const { + doublereal tnow = temperature(); + if (m_tlast != tnow) { + m_spthermo->update(tnow, m_cp0_R.begin(), m_h0_RT.begin(), + m_s0_R.begin()); + m_tlast = tnow; + } + } +} + + + + diff --git a/Cantera/src/SolidCompound.h b/Cantera/src/SolidCompound.h new file mode 100644 index 000000000..148060e6a --- /dev/null +++ b/Cantera/src/SolidCompound.h @@ -0,0 +1,187 @@ +/** + * + * @file SolidPhase.h + * + */ + +/* $Author$ + * $Date$ + * $Revision$ + * + * Copyright 2001 California Institute of Technology + * + */ + + +#ifndef CT_SOLIDPHASE_H +#define CT_SOLIDPHASE_H + +//#include "ct_defs.h" +#include "mix_defs.h" +#include "ThermoPhase.h" +#include "SpeciesThermo.h" + +namespace Cantera { + + /** + * @ingroup thermoprops + * + * Class IdealGasPhase represents low-density gases that obey the + * ideal gas equation of state. It derives from class ThermoPhase, + * and overloads the virtual methods defined there with ones that + * use expressions appropriate for ideal gas mixtures. + * + */ + class SolidCompound : public ThermoPhase { + + public: + + SolidCompound(): m_tlast(-1.0), m_tmin(0.0), m_tmax(0.0), + m_press(OneAtm), m_p0(OneAtm) {} + + virtual ~SolidCompound() {} + + /** + * Equation of state flag. Returns the value cIdealGas, defined + * in mix_defs.h. + */ + virtual int eosType() const { return cSolidCompound; } + + + /** + * @name Molar Thermodynamic Properties + * @{ + */ + + /** + * Molar enthalpy. Units: J/kmol. + */ + virtual doublereal enthalpy_mole() const { + double hh = intEnergy_mole() + m_press / molarDensity(); + return hh; + } + + /** + * Molar internal energy. J/kmol. + */ + virtual doublereal intEnergy_mole() const { + _updateThermo(); + // cout << "intEnergy: " << m_h0_RT[0] << " " << m_p0/molarDensity() + // << endl; + return GasConstant * temperature() * m_h0_RT[0] + - m_p0 / molarDensity(); + } + + /** + * Molar entropy. Units: J/kmol/K. + */ + virtual doublereal entropy_mole() const { + _updateThermo(); + //cout << "s/r = " << m_s0_R[0] << endl; + return GasConstant * m_s0_R[0]; + } + + + virtual doublereal gibbs_mole() const { + return enthalpy_mole() - temperature() * entropy_mole(); + } + + + /** + * Molar heat capacity at constant pressure. Units: J/kmol/K. + */ + virtual doublereal cp_mole() const { + _updateThermo(); + return GasConstant * m_cp0_R[0]; + } + + /** + * Molar heat capacity at constant volume. Units: J/kmol/K. + */ + virtual doublereal cv_mole() const { + return cp_mole(); + } + + //@} + + + /** + * @name Mechanical Equation of State + * @{ + */ + + /** + * Pressure. Units: Pa. + */ + virtual doublereal pressure() const { + return m_press; + } + + /** + * Set the pressure at constant temperature. Units: Pa. + */ + virtual void setPressure(doublereal p) { + m_press = p; + } + + //@} + + + virtual void getChemPotentials(doublereal* mu) const { + mu[0] = gibbs_mole(); + } + + virtual void getStandardChemPotentials(doublereal* mu0) const { + mu0[0] = gibbs_mole(); + //cout << m_h0_RT[0] << " " << m_s0_R[0] << endl; + //cout << "std chem pot = " << mu0[0] << endl; + } + + /** + * This method returns the array of generalized + * concentrations. For a solid compound, there is only one + * species, and the generalized concentration is 1.0. + */ + virtual void getActivityConcentrations(doublereal* c) const { + c[0] = 1.0; + } + + /** + * The standard concentration. This is defined as the concentration + * by which the generalized concentration is normalized to produce + * the activity. + */ + virtual doublereal standardConcentration(int k=0) const { + return 1.0; + } + + virtual doublereal logStandardConc(int k=0) const { + return 0.0; + } + + virtual void initThermo(); + + + protected: + + int m_kk; + doublereal m_tmin, m_tmax, m_p0, m_press; + + mutable doublereal m_tlast; + mutable array_fp m_h0_RT; + mutable array_fp m_cp0_R; + mutable array_fp m_s0_R; + + private: + + void _updateThermo() const; + }; + +} + +#endif + + + + + diff --git a/Cantera/src/SurfPhase.h b/Cantera/src/SurfPhase.h index fd29a7681..d76bf033e 100644 --- a/Cantera/src/SurfPhase.h +++ b/Cantera/src/SurfPhase.h @@ -52,7 +52,7 @@ namespace Cantera { void setSiteDensity(doublereal n0); void setElectricPotential(doublereal V); void setCoverages(const doublereal* theta); - void getCoverages(doublereal* theta); + void getCoverages(doublereal* theta) const; protected: diff --git a/Cantera/src/ThermoFactory.cpp b/Cantera/src/ThermoFactory.cpp index d13ce1bcc..0f6101e9d 100644 --- a/Cantera/src/ThermoFactory.cpp +++ b/Cantera/src/ThermoFactory.cpp @@ -21,15 +21,19 @@ #include "IdealGasPhase.h" #include "ConstDensityThermo.h" #include "SurfPhase.h" +#include "MetalPhase.h" +#include "SolidCompound.h" #include "importCTML.h" namespace Cantera { ThermoFactory* ThermoFactory::__factory = 0; - static int ntypes = 3; - static string _types[] = {"IdealGas", "Incompressible", "Surface"}; - static int _itypes[] = {cIdealGas, cIncompressible, cSurf}; + static int ntypes = 5; + static string _types[] = {"IdealGas", "Incompressible", + "Surface", "Metal", "SolidCompound"}; + static int _itypes[] = {cIdealGas, cIncompressible, + cSurf, cMetal, cSolidCompound}; ThermoPhase* ThermoFactory::newThermoPhase(string model) { @@ -55,6 +59,14 @@ namespace Cantera { th = new SurfPhase; break; + case cMetal: + th = new MetalPhase; + break; + + case cSolidCompound: + th = new SolidCompound; + break; + default: throw CanteraError("newThermo", "newThermo: unknown equation of state: "+model); @@ -99,22 +111,23 @@ namespace Cantera { switch (ieos) { case cIdealGas: - // Ideal gas th = new IdealGasPhase; break; case cIncompressible: - // getFloats(eos,d); - //dens = d["density"]; th = new ConstDensityThermo; - //th->setParameters(1, &dens); break; case cSurf: - //getFloats(eos,d); - //dens = d["site_density"]; th = new SurfPhase; - //th->setParameters(1, &dens); + break; + + case cMetal: + th = new MetalPhase; + break; + + case cSolidCompound: + th = new SolidCompound; break; default: @@ -122,8 +135,10 @@ namespace Cantera { "newThermo: unknown equation of state: "+eostype); } th->setSpeciesThermo(spthermo); + // import the phase specification importPhase(node, th); + return th; } diff --git a/Cantera/src/importCTML.cpp b/Cantera/src/importCTML.cpp index a4bb8d21f..00cd4743c 100755 --- a/Cantera/src/importCTML.cpp +++ b/Cantera/src/importCTML.cpp @@ -33,6 +33,7 @@ using namespace std; #include "ReactionData.h" #include "global.h" #include "stringUtils.h" +#include "GasKineticsWriter.h" #include "xml.h" #include "ctml.h" @@ -40,6 +41,8 @@ using namespace ctml; #include +GasKineticsWriter* writer = 0; + namespace Cantera { typedef vector nodeset_t; @@ -191,15 +194,15 @@ namespace Cantera { * Define a map and get all of the floats in the * current XML species block */ - map fd; - getFloats(s, fd); + //map fd; + //getFloats(s, fd); + doublereal chrg = 0.0; + if (s.hasChild("charge")) chrg = getFloat(s, "charge"); + doublereal sz = 1.0; + if (s.hasChild("size")) sz = getFloat(s, "size"); - /* - * Set a default for the size parameter to one - */ - if (fd["size"] == 0.0) fd["size"] = 1.0; p.addUniqueSpecies(s["name"], ecomp.begin(), - fd["charge"], fd["size"]); + chrg, sz); // get thermo XML_Node& thermo = s.child("thermo"); @@ -391,7 +394,8 @@ namespace Cantera { vector_fp coeff(3); getArrhenius(c, order, coeff[0], coeff[1], coeff[2]); if (order == 0) order = nr; - if (order == nr || rdata.reactionType == THREE_BODY_RXN) + if (order == nr || rdata.reactionType == THREE_BODY_RXN + || rdata.reactionType == ELEMENTARY_RXN) chigh = coeff; else if (order == nr + 1) clow = coeff; else { @@ -401,6 +405,15 @@ namespace Cantera { "wrong Arrhenius coeff order"); } } +// else if (nm == "Stick") { +// vector_fp coeff(3); +// string spname = c["species"]; +// ThermoPhase& th = kin.speciesPhase(spname); +// int isp = th.speciesIndex(spname); +// double mw = th.molecularWeights()[isp]; +// cbar = sqrt((8.0*GasConstant)/(Pi*mw)); + // + // } else if (nm == "falloff") { getFalloff(c, rdata); } @@ -501,11 +514,21 @@ namespace Cantera { "wrong equation of state type"); } } + else if (eos["model"] == "SolidCompound") { + if (th->eosType() == cSolidCompound) { + doublereal rho = getFloat(eos, "density", "-"); + th->setDensity(rho); + } + else { + throw CanteraError("importCTML", + "wrong equation of state type"); + } + } else if (eos["model"] == "Surface") { if (th->eosType() == cSurf) { - map d; + //map d; //getFloats(eos, d); - doublereal n = fpValue(eos("site_density")); + doublereal n = getFloat(eos, "site_density", "-"); th->setParameters(1, &n); } else { @@ -685,12 +708,16 @@ namespace Cantera { getRateCoefficient(r.child("rateCoeff"), kin, rdata); kin.addReaction(rdata); + //if (writer) writer->addReaction(rdata); return true; } bool installReactionArrays(XML_Node& p, Kinetics& kin, string default_phase) { + + writer = new GasKineticsWriter; + vector rarrays; int itot = 0; p.getChildren("reactionArray",rarrays); @@ -738,6 +765,12 @@ namespace Cantera { } } kin.finalize(); + ofstream fwrite("mech.cpp"); + //writer->writeGetNetProductionRates(cout, kin.nTotalSpecies(), + // kin.nReactions()); + fwrite.close(); + delete writer; + writer = 0; return true; } diff --git a/Cantera/src/mix_defs.h b/Cantera/src/mix_defs.h index 335b49424..4613b0d1a 100755 --- a/Cantera/src/mix_defs.h +++ b/Cantera/src/mix_defs.h @@ -13,6 +13,8 @@ namespace Cantera { const int cIdealGas = 1; const int cIncompressible = 2; const int cSurf = 3; + const int cMetal = 4; + const int cSolidCompound = 5; // kinetic manager types const int cGasKinetics = 2; diff --git a/Cantera/src/oneD/Domain1D.h b/Cantera/src/oneD/Domain1D.h index 6707b77e1..6ecdde204 100644 --- a/Cantera/src/oneD/Domain1D.h +++ b/Cantera/src/oneD/Domain1D.h @@ -107,10 +107,10 @@ namespace Cantera { m_nv = nv; m_max.resize(m_nv, 0.0); m_min.resize(m_nv, 0.0); - m_rtol_ss.resize(m_nv, 0.0); - m_atol_ss.resize(m_nv, 0.0); - m_rtol_ts.resize(m_nv, 0.0); - m_atol_ts.resize(m_nv, 0.0); + m_rtol_ss.resize(m_nv, 1.0e-8); + m_atol_ss.resize(m_nv, 1.0e-15); + m_rtol_ts.resize(m_nv, 1.0e-8); + m_atol_ts.resize(m_nv, 1.0e-15); m_points = np; m_z.resize(np, 0.0); m_slast.resize(m_nv * m_points, 0.0); diff --git a/Cantera/src/oneD/MultiNewton.cpp b/Cantera/src/oneD/MultiNewton.cpp index 637410393..c02ae41eb 100644 --- a/Cantera/src/oneD/MultiNewton.cpp +++ b/Cantera/src/oneD/MultiNewton.cpp @@ -97,14 +97,17 @@ namespace Cantera { const doublereal* step, OneDim& r) const { doublereal f, sum = 0.0;//, fmx = 0.0; int n; - int nd = r.nDomains(); + int nd = r.nDomains(); for (n = 0; n < nd; n++) { f = norm_square(x + r.start(n), step + r.start(n), r.domain(n)); sum += f; - // if (f > fmx) fmx = f; + // cout << "n = " << n << " f = " << f << endl; } + //cout << "sum = " << sum << endl; + //cout << "r.size() = " << r.size() << endl; sum /= r.size(); + //cout << "sum = " << sum << " " << sqrt(sum) << endl; return sqrt(sum); } @@ -121,7 +124,11 @@ namespace Cantera { for (n = 0; n < sz; n++) { step[n] = -step[n]; } + jac.solve(sz, step, step); + +#undef DEBUG_STEP #ifdef DEBUG_STEP + bool ok = false; Domain1D* d; if (!ok) { for (n = 0; n < sz; n++) { @@ -132,10 +139,10 @@ namespace Cantera { r.pointDomain(n)->componentName(n - d->loc() - nvd*pt) << " " << x[n] << " " << step[n] << endl; } - if (!ok) throw "not ok"; + //if (!ok) throw "not ok"; } #endif - jac.solve(sz, step, step); + } diff --git a/Cantera/src/oneD/Sim1D.cpp b/Cantera/src/oneD/Sim1D.cpp index a36fcf0ae..da35dc12f 100644 --- a/Cantera/src/oneD/Sim1D.cpp +++ b/Cantera/src/oneD/Sim1D.cpp @@ -336,9 +336,9 @@ namespace Cantera { } } } - //else { - // throw CanteraError("refine","keepPoint is false at m = "+int2str(m)); - //} + else { + throw CanteraError("refine","keepPoint is false at m = "+int2str(m)); + } } dsize.push_back(znew.size() - nstart); } diff --git a/Cantera/src/oneD/Solid1D.cpp b/Cantera/src/oneD/Solid1D.cpp new file mode 100644 index 000000000..139995ceb --- /dev/null +++ b/Cantera/src/oneD/Solid1D.cpp @@ -0,0 +1,634 @@ +/** + * @file Solid.cpp + */ + +/* + * $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2003 California Institute of Technology + + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include +#include + +#include "Solid1D.h" +#include "../ArrayViewer.h" +#include "ctml.h" +#include "MultiJac.h" + +using namespace ctml; + +namespace Cantera { + + + int Solid1D::c_T_loc = 0; + int Solid1D::c_phi_loc = 1; + int Solid1D::c_Y_loc = 2; + + + Solid1D::Solid1D(ThermoPhase* ph, int nsp, int points) : + Domain1D(nsp+1, points), + m_nsp(nsp), + m_thermo(0), + m_kin(0), + m_trans(0), + m_jac(0), + m_ok(false), + { + m_type = cSolidType; + + m_points = points; + m_thermo = ph; + + if (ph == 0) return; // used to create a dummy object + + int nsp2 = m_thermo->nSpecies(); + if (nsp2 != m_nsp) { + m_nsp = nsp2; + Domain1D::resize(m_nsp+1, points); + } + + + // make a local copy of the species molecular weight vector + m_wt = m_thermo->molecularWeights(); + + m_nv = m_nsp + 1; + + // turn off the energy equation at all points + m_do_energy.resize(m_points,false); + m_do_gauss.resize(m_points,false); + m_do_species.resize(m_nsp,false); + + m_diff.resize(m_nsp*m_points); + m_flux.resize(m_nsp,m_points); + m_wdot.resize(m_nsp,m_points, 0.0); + m_ybar.resize(m_nsp); + + + //-------------- default solution bounds -------------------- + + vector_fp vmin(m_nv), vmax(m_nv); + + // temperature bounds + vmin[c_T_loc] = 200.0; + vmax[c_T_loc]= 1.e9; + + // mass fraction bounds + int k; + for (k = 0; k < m_nsp; k++) { + vmin[c_Y_loc + k] = -1.0e-5; + vmax[c_Y_loc + k] = 1.0e5; + } + setBounds(vmin.size(), vmin.begin(), vmax.size(), vmax.begin()); + + + //-------------------- default error tolerances ---------------- + vector_fp rtol(m_nv, 1.0e-8); + vector_fp atol(m_nv, 1.0e-15); + setTolerances(rtol.size(), rtol.begin(), atol.size(), atol.begin(),false); + setTolerances(rtol.size(), rtol.begin(), atol.size(), atol.begin(),true); + + //-------------------- grid refinement ------------------------- + m_refiner->setActive(c_T_loc, false); + + vector_fp gr; + for (int ng = 0; ng < m_points; ng++) gr.push_back(1.0*ng/m_points); + setupGrid(m_points, gr.begin()); + setID("solid"); + } + + + /** + * Change the grid size. Called after grid refinement. + */ + void Solid1D::resize(int points) { + Domain1D::resize(m_nv, points); + + m_rho.resize(m_points, 0.0); + m_wtm.resize(m_points, 0.0); + m_cp.resize(m_points, 0.0); + m_tcon.resize(m_points, 0.0); + m_cdens.resize(m_points, 0.0); + m_diff.resize(m_nsp*m_points); + + m_flux.resize(m_nsp,m_points); + m_wdot.resize(m_nsp,m_points, 0.0); + + m_do_energy.resize(m_points,false); + m_do_gauss.resize(m_points,false); + m_do_species.resize(m_nsp,false); + + m_fixedtemp.resize(m_points); + m_fixedphi.resize(m_points); + + m_dz.resize(m_points-1); + m_z.resize(m_points); + } + + + + void Solid1D::setupGrid(int n, const doublereal* z) { + resize(n); + int j; + m_z[0] = z[0]; + for (j = 1; j < m_points; j++) { + m_z[j] = z[j]; + m_dz[j-1] = m_z[j] - m_z[j-1]; + } + } + + + /** + * Install a transport manager. + */ + void Solid1D::setTransport(Transport& trans) { + m_trans = &trans; + + if (m_trans->model() != cSolidTransport) { + throw CanteraError("setTransport","unknown transport model."); + } + + + /** + * Set the solid object state to be consistent with the solution at + * point j. + */ + void Solid1D::setThermoState(const doublereal* x,int j) { + m_thermo->setTemperature(T(x,j)); + m_thermo->setElectricPotential(phi(x,j)); + const doublereal* yy = x + m_nv*j + 1; + m_thermo->setMassFractions_NoNorm(yy); + } + + + /** + * Set the state to be consistent with the solution at the + * midpoint between j and j + 1. + */ + void Solid1D::setStateAtMidpoint(const doublereal* x,int j) { + m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1))); + m_thermo->setElectricPotential(0.5*(phi(x,j)+phi(x,j+1))); + const doublereal* yyj = x + m_nv*j + 1; + const doublereal* yyjp = x + m_nv*(j+1) + 1; + for (int k = 0; k < m_nsp; k++) + m_ybar[k] = 0.5*(yyj[k] + yyjp[k]); + m_thermo->setMassFractions_NoNorm(m_ybar.begin()); + } + + + + void Solid1D::eval(int jg, doublereal* xg, + doublereal* rg, integer* diagg, doublereal rdt) { + + // if evaluating a Jacobian, and the global point is outside + // the domain of influence for this domain, then skip + // evaluating the residual + if (jg >=0 && (jg < firstPoint() - 1 || jg > lastPoint() + 1)) return; + + // if evaluating a Jacobian, compute the steady-state residual + if (jg >= 0) rdt = 0.0; + + // start of local part of global arrays + doublereal* x = xg + loc(); + doublereal* rsd = rg + loc(); + integer* diag = diagg + loc(); + + int jmin, jmax, jpt; + jpt = jg - firstPoint(); + + if (jg < 0) { // evaluate all points + jmin = 0; + jmax = m_points - 1; + } + else { // evaluate points for Jacobian + jmin = max(jpt-1, 0); + jmax = min(jpt+1,m_points-1); + } + + // properties are computed for grid points from j0 to j1 + int j0 = max(jmin-1,0); + int j1 = min(jmax+1,m_points-1); + + + int j, k; + + + //----------------------------------------------------- + // update properties + //----------------------------------------------------- + + // thermodynamic properties only if a Jacobian is + // not being evaluated + if (jpt < 0) + updateThermo(x, j0, j1); + + // update transport properties only if a Jacobian is + // not being evaluated + if (jpt < 0) + updateTransport(x, j0, j1); + + // update the species diffusive mass fluxes whether or not a + // Jacobian is being evaluated + updateDiffFluxes(x, j0, j1); + + for (j = j0; j <= j1; j++) { + setThermoState(j); + m_cdens[j] = m_thermo->chargeDensity(); + } + + //---------------------------------------------------- + // evaluate the residual equations at all required + // grid points + //---------------------------------------------------- + + for (j = jmin; j <= jmax; j++) { + + + //---------------------------------------------- + // left boundary + //---------------------------------------------- + + if (j == 0) { + rsd[index(c_T_loc,0)] = T(x,0); + rsd[index(c_phi_loc,0)] = phi(x,0); + + // The default boundary condition for species is zero + // flux. However, the boundary object may modify + // this. + for (k = 0; k < m_nsp; k++) { + rsd[index(c_Y_loc + k, 0)] = - m_flux(k,0); + } + } + + + //---------------------------------------------- + // + // right boundary + // + //---------------------------------------------- + + else if (j == m_points - 1) { + rsd[index(c_T_loc,j)] = T(x,j); + rsd[index(c_phi_loc,j)] = phi(x,j); + doublereal sum = 0.0; + for (k = 0; k < m_nsp; k++) { + sum += Y(x,k,j); + rsd[index(k+c_Y_loc,j)] = m_flux(k,j-1); + } + rsd[index(c_Y_loc,j)] = 1.0 - sum; + diag[index(c_Y_loc,j)] = 0; + } + + + //------------------------------------------ + // interior points + //------------------------------------------ + + else { + + //------------------------------------------------- + // Species equations + // + // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k + // + //------------------------------------------------- + getWdot(x,j); + + doublereal convec, diffus; + for (k = 0; k < m_nsp; k++) { + diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) + /(z(j+1) - z(j-1)); + rsd[index(c_Y_loc + k, j)] + = (m_wt[k]*(wdot(k,j) ) + - diffus)/m_rho[j] + - rdt*(Y(x,k,j) - Y_prev(k,j)); + diag[index(c_Y_loc + k, j)] = 1; + } + + //----------------------------------------------- + // energy equation + //----------------------------------------------- + + if (m_do_energy[j]) { + + rsd[index(c_T_loc, j)] = - divHeatFlux(x,j); + rsd[index(c_T_loc, j)] /= (m_rho[j]*m_cp[j]); + + rsd[index(c_T_loc, j)] -= rdt*(T(x,j) - T_prev(j)); + diag[index(c_T_loc, j)] = 1; + } + + //---------------------------------------------- + // Gauss's equation + //---------------------------------------------- + rsd[index(c_phi_loc, j)] = m_cdens[j] - divDisplCurr(x,j); + + } + + // residual equations if the energy or species equations + // are disabled + + if (!m_do_energy[j]) { + rsd[index(c_T_loc, j)] = T(x,j) - T_fixed(j); + diag[index(c_T_loc, j)] = 0; + } + if (!m_do_gauss[j]) { + rsd[index(c_phi_loc, j)] = phi(x,j) - phi_fixed(j); + diag[index(c_phi_loc, j)] = 0; + } + } + } + + + + /** + * Update the transport properties at grid points in the range + * from j0 to j1, based on solution x. + */ + void Surf1D::updateTransport(doublereal* x,int j0, int j1) { + int j; + for (j = j0; j < j1; j++) { + setStateAtMidpoint(x,j); + m_trans->getMixDiffCoeffs(m_diff.begin() + j*m_nsp); + m_tcon[j] = m_trans->thermalConductivity(); + } + } + + + /** + * Print the solution. + */ + void Solid1D::showSolution(const doublereal* x) { + int nn = m_nv/5; + int i, j, n; + char* buf = new char[100]; + + // The mean molecular weight is needed to convert + updateThermo(x, 0, m_points-1); + + for (i = 0; i < nn; i++) { + drawline(); + sprintf(buf, "\n z "); + writelog(buf); + for (n = 0; n < 5; n++) { + sprintf(buf, " %10s ",componentName(i*5 + n).c_str()); + writelog(buf); + } + drawline(); + for (j = 0; j < m_points; j++) { + sprintf(buf, "\n %10.4g ",m_z[j]); + writelog(buf); + for (n = 0; n < 5; n++) { + sprintf(buf, " %10.4g ",component(x, i*5+n,j)); + writelog(buf); + } + } + writelog("\n"); + } + int nrem = m_nv - 5*nn; + drawline(); + sprintf(buf, "\n z "); + writelog(buf); + for (n = 0; n < nrem; n++) { + sprintf(buf, " %10s ", componentName(nn*5 + n).c_str()); + writelog(buf); + } + drawline(); + for (j = 0; j < m_points; j++) { + sprintf(buf, "\n %10.4g ",m_z[j]); + writelog(buf); + for (n = 0; n < nrem; n++) { + sprintf(buf, " %10.4g ",component(x, nn*5+n,j)); + writelog(buf); + } + } + writelog("\n"); + } + + + /** + * Update the diffusive mass fluxes. + */ + void Solid1D::updateDiffFluxes(const doublereal* x, int j0, int j1) { + int j, k, m; + doublereal sum, wtm, rho, dz, gradlogT, s; + doublereal dphidz, a1; + for (j = j0; j < j1; j++) { + sum = 0.0; + wtm = m_wtm[j]; + rho = density(j); + dz = z(j+1) - z(j); + dphidz = (phi(x,j+1) - phi(x,j))/dz; + a1 = rho*Faraday*dphidz/(GasConstant * T(x,j)); + for (k = 0; k < m_nsp; k++) { + m_flux(k,j) = m_wt[k]*(rho*m_diff[k+m_nsp*j]/wtm); + m_flux(k,j) *= (X(x,k,j) - X(x,k,j+1))/dz; + m_flux(k,j) += a1*0.5*(Y(x,k,j) + + Y(x,k,j+1))*m_diff[k+m_nsp*j]*m_charge[k]; + sum -= m_flux(k,j); + } + for (k = 0; k < m_nsp; k++) m_flux(k,j) += Y(x,k,j)*sum; + } + break; + } + + + void Solid1D::outputTEC(ostream &s, const doublereal* x, + string title, int zone) { + int j,k; + s << "TITLE = \"" + title + "\"" << endl; + s << "VARIABLES = \"Z (m)\"" << endl; + s << "\"T (K)\"" << endl; + s << "\"phi (V)\"" << endl; + + for (k = 0; k < m_nsp; k++) { + s << "\"" << m_thermo->speciesName(k) << "\"" << endl; + } + s << "ZONE T=\"c" << zone << "\"" << endl; + s << " I=" << m_points << ",J=1,K=1,F=POINT" << endl; + s << "DT=(SINGLE SINGLE"; + for (k = 0; k < m_nsp; k++) s << " SINGLE"; + s << " )" << endl; + for (j = 0; j < m_points; j++) { + s << z(j) << " "; + for (k = 0; k < m_nv; k++) { + s << component(x, k, j) << " "; + } + s << endl; + } + } + + + string Solid1D::componentName(int n) const { + switch(n) { + case c_T_loc: return "T"; + case c_phi_loc: return "phi"; + default: + if (n >= (int) 1 && n < (int) (c_Y_loc + m_nsp)) { + return m_thermo->speciesName(n - 1); + } + else + return ""; + } + } + + + void Solid1D::restore(XML_Node& dom, doublereal* soln) { + + vector ignored; + int nsp = m_thermo->nSpecies(); + vector_int did_species(nsp, 0); + + vector str; + dom.getChildren("string",str); + int nstr = str.size(); + for (int istr = 0; istr < nstr; istr++) { + XML_Node& nd = *str[istr]; + writelog(nd["title"]+": "+nd.value()+"\n"); + } + + map params; + getFloats(dom, params); + + vector d; + dom.child("grid_data").getChildren("floatArray",d); + int nd = d.size(); + + vector_fp x; + int n, np, j, ks, k; + string nm; + bool readgrid = false, wrote_header = false; + for (n = 0; n < nd; n++) { + XML_Node& fa = *d[n]; + nm = fa["title"]; + if (nm == "z") { + getFloatArray(fa,x,false); + np = x.size(); + writelog("Grid contains "+int2str(np)+ + " points.\n"); + readgrid = true; + + // note that setupGrid also resizes the domain. + setupGrid(np, x.begin()); + } + } + if (!readgrid) { + throw CanteraError("Solid1D::restore", + "domain contains no grid points."); + } + + writelog("Importing datasets:\n"); + for (n = 0; n < nd; n++) { + XML_Node& fa = *d[n]; + nm = fa["title"]; + getFloatArray(fa,x,false); + if (nm == "z") { + ; // already read grid + } + else if (nm == "T") { + writelog("temperature "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) + soln[index(c_T_loc,j)] = x[j]; + + // For fixed-temperature simulations, use the imported temperature profile by default. + // If this is not desired, call setFixedTempProfile *after* restoring the solution. + vector_fp zz(np); + for (int jj = 0; jj < np; jj++) zz[jj] = (grid(jj) - zmin())/(zmax() - zmin()); + setFixedTempProfile(zz, x); + } + else goto error; + } + else if (nm == "phi") { + writelog("potential "); + if ((int) x.size() == np) { + for (j = 0; j < np; j++) + soln[index(c_phi_loc,j)] = x[j]; + } + else goto error; + } + else if (m_thermo->speciesIndex(nm) >= 0) { + writelog(nm+" "); + if ((int) x.size() == np) { + k = m_thermo->speciesIndex(nm); + did_species[k] = 1; + for (j = 0; j < np; j++) + soln[index(k+c_Y_loc,j)] = x[j]; + } + } + else + ignored.push_back(nm); + } + + if (ignored.size() != 0) { + writelog("\n\n"); + writelog("Ignoring datasets:\n"); + int nn = ignored.size(); + for (int n = 0; n < nn; n++) { + writelog(ignored[n]+" "); + } + } + + for (ks = 0; ks < nsp; ks++) { + if (did_species[ks] == 0) { + if (!wrote_header) { + writelog("Missing data for species:\n"); + wrote_header = true; + } + writelog(m_thermo->speciesName(ks)+" "); + } + } + + return; + error: + throw CanteraError("Solid1D::restore","Data size error"); + } + + + + void Solid1D::save(XML_Node& o, doublereal* sol) { + int k; + + ArrayViewer soln(m_nv, m_points, sol + loc()); + + XML_Node& flow = (XML_Node&)o.addChild("domain"); + flow.addAttribute("type",flowType()); + flow.addAttribute("id",m_id); + flow.addAttribute("points",m_points); + flow.addAttribute("components",m_nv); + + if (m_desc != "") addString(flow,"description",m_desc); + XML_Node& gv = flow.addChild("grid_data"); + addFloatArray(gv,"z",m_z.size(),m_z.begin(), + "m","length"); + vector_fp x(soln.nColumns()); + + soln.getRow(c_T_loc,x.begin()); + addFloatArray(gv,"T",x.size(),x.begin(),"K","temperature",0.0); + + soln.getRow(c_phi_loc,x.begin()); + addFloatArray(gv,"phi",x.size(),x.begin(),"V","potential",0.0); + + for (k = 0; k < m_nsp; k++) { + soln.getRow(c_Y_loc+k,x.begin()); + addFloatArray(gv,m_thermo->speciesName(k), + x.size(),x.begin(),"","massFraction",0.0,1.0); + } + } + + + void Solid1D::setJac(MultiJac* jac) { + m_jac = jac; + } + + +} diff --git a/Cantera/src/oneD/Solid1D.h b/Cantera/src/oneD/Solid1D.h new file mode 100644 index 000000000..b8dcbbac3 --- /dev/null +++ b/Cantera/src/oneD/Solid1D.h @@ -0,0 +1,402 @@ +§/** + * @file Solid1D.h + * + */ + +/* + * $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2001 California Institute of Technology + +#ifndef CT_SOLID1D_H +#define CT_SOLID1D_H + +#include "../transport/TransportBase.h" +#include "Domain1D.h" +#include "../Array.h" +#include "../sort.h" +#include "../ThermoPhase.h" +#include "../Kinetics.h" +#include "../funcs.h" + + +namespace Cantera { + + class MultiJac; + + + //----------------------------------------------------------- + // Class Solid1D + //----------------------------------------------------------- + + + /** + * A class for one-dimensional reacting solids with current + * transport. This class implements the one-dimensional + * similarity solution for a chemically-reacting, axisymmetric, + * stagnation-point flow. + */ + class Solid1D : public Domain1D { + + public: + + + //------------------------------------------ + // constants + //------------------------------------------ + + /** + * Offsets of solution components in the solution array. + */ + const unsigned int c_phi_loc; // electric potential + const unsigned int c_T_loc; // temperature + const unsigned int c_C_loc; // concentrations + + + //-------------------------------- + // construction and destruction + //-------------------------------- + + // Constructor. + Solid1D(ThermoPhase* ph = 0, int nsp = 1, int points = 1); + + /// Destructor. + virtual ~Solid1D(){} + + + /** + * @name Problem Specification + */ + //@{ + + virtual void setupGrid(int n, const doublereal* z); + + thermo_t& phase() { return *m_thermo; } + kinetics_t& kinetics() { return *m_kin; } + + /** + * Set the thermo manager. + */ + void setThermo(thermo_t& th) { + m_thermo = &th; + } + + /// set the kinetics manager + void setKinetics(kinetics_t& kin) { m_kin = &kin; } + + /// set the transport manager + void setTransport(Transport& trans); + + virtual void setState(int point, const doublereal* state) { + setTemperature(point, state[c_T_loc]); + setElectricPotential(point, state[c_phi_loc); + int k; + for (k = 0; k < m_nsp; k++) { + setConcentration(point, k, state[c_C_loc+k]); + } + } + + + virtual void _getInitialSoln(doublereal* x) { + int k, j; + for (j = 0; j < m_points; j++) { + x[index(c_T_loc,j)] = T_fixed(j); + x[index(c_phi_loc,j)] = phi_fixed(j); + for (k = 0; k < m_nsp; k++) { + x[index(c_C_loc+k,j)] = C_fixed(k,j); + } + } + } + + virtual void _finalize(const doublereal* x) { + int k, j; + doublereal zz, tt; + int nz = m_zfix.size(); + bool e = m_do_energy[0]; + for (j = 0; j < m_points; j++) { + if (e || nz == 0) + setTemperature(j, T(x, j)); + else { + zz = (z(j) - z(0))/(z(m_points - 1) - z(0)); + tt = linearInterp(zz, m_zfix, m_tfix); + setTemperature(j, tt); + } + setElectricPotential(j, phi(x,j)); + for (k = 0; k < m_nsp; k++) { + setConcentration(j, k, C(x, k, j)); + } + } + if (e) solveEnergyEqn(); + } + + + void setFixedTempProfile(vector_fp& zfixed, vector_fp& tfixed) { + m_zfix = zfixed; + m_tfix = tfixed; + } + + /** + * Set the temperature fixed point at grid point j, and + * disable the energy equation so that the solution will be + * held to this value. + */ + void setTemperature(int j, doublereal t) { + m_fixedtemp[j] = t; + m_do_energy[j] = false; + } + + /** + * Set the electric potential fixed point at grid point j, and + * disable Gauss's equation so that the solution will be + * held to this value. + */ + void setElectricPotential(int j, doublereal phi) { + m_fixedphi[j] = phi; + m_do_gauss[j] = false; + } + + /** + * Set the mass fraction fixed point for species k at grid + * point j, and disable the species equation so that the + * solution will be held to this value. + */ + void setConcentration(int j, int k, doublereal c) { + m_fixedc(k,j) = c; + m_do_species[k] = true; // false; + } + + /** + * The fixed temperature value at point j. + */ + doublereal T_fixed(int j) const {return m_fixedtemp[j];} + + /** + * The fixed potential value at point j. + */ + doublereal phi_fixed(int j) const {return m_fixedphi[j];} + + /** + * The fixed mass fraction value of species k at point j. + */ + doublereal C_fixed(int k, int j) const {return m_fixedc(k,j);} + + virtual string componentName(int n) const; + + void setDielectricConstant(doublereal e) { m_eps = e; } + doublereal dielectricConstant() { return e; } + + /** + * Write a Tecplot zone corresponding to the current solution. + * May be called multiple times to generate animation. + */ + void outputTEC(ostream &s, const doublereal* x, + string title, int zone); + + virtual void showSolution(const doublereal* x); + + virtual void save(XML_Node& o, doublereal* sol); + + virtual void restore(XML_Node& dom, doublereal* soln); + + // overloaded in subclasses + virtual string solidType() { return ""; } + + void solveEnergyEqn(int j=-1) { + if (j < 0) + for (int i = 0; i < m_points; i++) + m_do_energy[i] = true; + else + m_do_energy[j] = true; + m_refiner->setActive(c_T_loc, true); + needJacUpdate(); + } + + void fixTemperature(int j=-1) { + if (j < 0) + for (int i = 0; i < m_points; i++) { + m_do_energy[i] = false; + } + else m_do_energy[j] = false; + m_refiner->setActive(c_T_loc, false); + needJacUpdate(); + } + + void solveGaussEqn(int j=-1) { + if (j < 0) + for (int i = 0; i < m_points; i++) + m_do_gauss[i] = true; + else + m_do_gauss[j] = true; + m_refiner->setActive(c_phi_loc, true); + needJacUpdate(); + } + + void fixElectricPotential(int j=-1) { + if (j < 0) + for (int i = 0; i < m_points; i++) { + m_do_gauss[i] = false; + } + else m_do_gauss[j] = false; + m_refiner->setActive(c_phi_loc, false); + needJacUpdate(); + } + + bool doSpecies(int k) { return m_do_species[k]; } + bool doEnergy(int j) { return m_do_energy[j]; } + bool doGauss(int j) { return m_do_gauss[j]; } + + void solveSpecies(int k=-1) { + if (k == -1) { + for (int i = 0; i < m_nsp; i++) + m_do_species[i] = true; + } + else m_do_species[k] = true; + needJacUpdate(); + } + + void fixSpecies(int k=-1) { + if (k == -1) { + for (int i = 0; i < m_nsp; i++) + m_do_species[i] = false; + } + else m_do_species[k] = false; + needJacUpdate(); + } + + void resize(int points); + + void setJac(MultiJac* jac); + void setThermoState(const doublereal* x,int j); + void setStateAtMidpoint(const doublereal* x,int j); + + + protected: + + doublereal component(const doublereal* x, int i, int j) const { + doublereal xx = x[index(i,j)]; + return xx; + } + + doublereal wdot(int k, int j) const {return m_wdot(k,j);} + + /// write the net production rates at point j into array m_wdot + void getWdot(doublereal* x,int j) { + setThermoState(x,j); + m_kin->getNetProductionRates(&m_wdot(0,j)); + } + + /** + * update the thermodynamic properties from point + * j0 to point j1 (inclusive), based on solution x. + */ + void updateThermo(const doublereal* x, int j0, int j1) { + int j; + for (j = j0; j <= j1; j++) { + setThermoState(x,j); + m_cp[j] = m_thermo->cp_mass(); + } + } + + + //-------------------------------- + // solution components + //-------------------------------- + + + doublereal T(const doublereal* x,int j) const { + return x[index(c_T_loc, j)]; + } + + doublereal& T(doublereal* x,int j) {return x[index(c_T_loc, j)];} + + doublereal T_prev(int j) const {return prevSoln(c_T_loc, j);} + + doublereal C(const doublereal* x,int k, int j) const { + return x[index(c_C_loc + k, j)]; + } + + doublereal& C(doublereal* x,int k, int j) { + return x[index(c_C_loc + k, j)]; + } + + doublereal C_prev(int k, int j) const { + return prevSoln(c_C_loc + k, j); + } + + doublereal flux(int k, int j) const { + return m_flux(k, j); + } + + doublereal phi(doublereal* x, j) { + return x[index(c_phi_loc, j)]; + } + + doublereal divHeatFlux(const doublereal* x, int j) const { + doublereal c1 = m_tcon[j-1]*(T(x,j) - T(x,j-1)); + doublereal c2 = m_tcon[j]*(T(x,j+1) - T(x,j)); + return -2.0*(c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/(z(j+1) - z(j-1)); + } + + doublereal divDisplCurr(const doublereal* x, int j) const { + doublereal c1 = (phi(x,j) - phi(x,j-1)); + doublereal c2 = (phi(x,j+1) - phi(x,j)); + return -2.0*m_eps*epsilon_0* + (c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/(z(j+1) - z(j-1)); + } + + void updateDiffFluxes(const doublereal* x, int j0, int j1); + + //--------------------------------------------------------- + // + // member data + // + //--------------------------------------------------------- + + doublereal m_eps // relative dielectric constant + + // grid parameters + vector_fp m_dz; + + // mixture thermo properties + vector_fp m_cdens; + + // transport properties + vector_fp m_tcon; + vector_fp m_diff; + Array2D m_flux; + + // production rates + Array2D m_wdot; + + int m_nsp; + + thermo_t* m_thermo; + kinetics_t* m_kin; + Transport* m_trans; + + MultiJac* m_jac; + + bool m_ok; + + // flags + vector m_do_energy; + vector m_do_species; + vector m_do_gauss; + + // fixed T and Y values + Array2D m_fixedy; + Array2D m_fixedphi; + vector_fp m_fixedtemp; + vector_fp m_zfix; + vector_fp m_tfix; + + private: + + vector_fp m_cbar; + }; +} + +#endif diff --git a/Cantera/src/oneD/refine.cpp b/Cantera/src/oneD/refine.cpp index 5f400a864..ec5ea517e 100644 --- a/Cantera/src/oneD/refine.cpp +++ b/Cantera/src/oneD/refine.cpp @@ -58,10 +58,10 @@ namespace Cantera { m_nv = m_domain->nComponents(); // check consistency - if (n != m_domain->nPoints()) return -1; + if (n != m_domain->nPoints()) throw CanteraError("analyze","inconsistent"); - if (n >= m_npmax) return 0; + if (n >= m_npmax) throw CanteraError("analyze","max points"); /** * find locations where cell size ratio is too large. @@ -88,7 +88,7 @@ namespace Cantera { // } for (int i = 0; i < m_nv; i++) { - //cout << i << " " << m_nv << " " << m_active[i] << endl; + cout << i << " " << m_nv << " " << m_active[i] << endl; if (m_active[i]) { name = m_domain->componentName(i); @@ -129,7 +129,7 @@ namespace Cantera { m_c[name] = 1; if (int(m_loc.size()) + n > m_npmax) goto done; } - if (r >= 0.0) { + if (r >= -1.0) { m_keep[j] = 1; m_keep[j+1] = 1; } @@ -155,7 +155,7 @@ namespace Cantera { m_loc[j+1] = 1; if (int(m_loc.size()) + n > m_npmax) goto done; } - if (r >= 0.0) { + if (r >= -1.0) { m_keep[j+1] = 1; } //cout << "at point " << j << " slope r = " diff --git a/Cantera/src/transport/SolidTransport.cpp b/Cantera/src/transport/SolidTransport.cpp new file mode 100644 index 000000000..dd3c499a3 --- /dev/null +++ b/Cantera/src/transport/SolidTransport.cpp @@ -0,0 +1,84 @@ +/** + * + * @file SolidTransport.cpp + */ + +/* $Author$ + * $Revision$ + * $Date$ + */ + +// copyright 2003 California Institute of Technology + + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +#include "SolidTransport.h" + +#include "utilities.h" +#include + + +namespace Cantera { + + //////////////////// class SolidTransport methods ////////////// + + SolidTransport::SolidTransport() {} + + void SolidTransport::setParameters(int n, int k, double* p) { + switch (n) { + case 0: + m_sp.push_back(k); + m_Adiff.push_back(p[0]); + m_Ndiff.push_back(p[1]); + m_Ediff.push_back(p[2]); + m_nmobile = m_sp.size(); + break; + case 1: + m_lam = p[0]; + break; + default: + ; + } + } + + + /********************************************************* + * + * Public methods + * + *********************************************************/ + + + void SolidTransport::getMobilities(doublereal* mobil) { + int k; + getMixDiffCoeffs(mobil); + doublereal t = m_thermo->temperature(); + int nsp = m_thermo->nSpecies(); + doublereal c1 = ElectronCharge / (Boltzmann * t); + for (k = 0; k < nsp; k++) { + mobil[k] *= c1 * m_thermo->charge(k); + } + } + + + doublereal SolidTransport::thermalConductivity() { + return m_lam; + } + + + void SolidTransport::getMixDiffCoeffs(doublereal* d) { + doublereal temp = m_thermo->temperature(); + int nsp = m_thermo->nSpecies(); + int k; + for (k = 0; k < nsp; k++) d[k] = 0.0; + for (k = 0; k < m_nmobile; k++) { + d[m_sp[k]] = + m_Adiff[k] * pow(temp, m_Ndiff[k]) * exp(-m_Ediff[k]/temp); + } + } +} diff --git a/Cantera/src/transport/SolidTransport.h b/Cantera/src/transport/SolidTransport.h new file mode 100644 index 000000000..b03ff69ed --- /dev/null +++ b/Cantera/src/transport/SolidTransport.h @@ -0,0 +1,81 @@ +/** + * + * @file SolidTransport.h + * Header file defining class SolidTransport + */ + +/* $Author$ + * $Revision$ + * $Date$ + */ + +// Copyright 2003 California Institute of Technology + + +#ifndef CT_SOLIDTRAN_H +#define CT_SOLIDTRAN_H + + +// turn off warnings under Windows +#ifdef WIN32 +#pragma warning(disable:4786) +#pragma warning(disable:4503) +#endif + +// STL includes +#include +#include +#include +#include +#include + +using namespace std; + +// Cantera includes +#include "TransportBase.h" +#include "../DenseMatrix.h" + +namespace Cantera { + + + /** + * Class SolidTransport implements transport + * properties for solids. + */ + class SolidTransport : public Transport { + + public: + + virtual ~SolidTransport() {} + + virtual int model() { return cSolidTransport; } + + virtual doublereal thermalConductivity(); + virtual void getMixDiffCoeffs(doublereal* d); + virtual void getMobilities(doublereal* mobil); + virtual void setParameters(int n, int k, doublereal* p); + + friend class TransportFactory; + + protected: + + /// default constructor + SolidTransport(); + + private: + + int m_nmobile; // number of mobile species + vector_fp m_Adiff; + vector_fp m_Ndiff; + vector_fp m_Ediff; + vector_int m_sp; + doublereal m_lam; + }; +} +#endif + + + + + + diff --git a/Cantera/src/units.h b/Cantera/src/units.h index 1e70f1e9d..f599a444f 100644 --- a/Cantera/src/units.h +++ b/Cantera/src/units.h @@ -22,7 +22,7 @@ namespace Cantera { if (units == "") return 1.0; doublereal f = 1.0, fctr; int tsize; - string u = units, tok; + string u = units, tok, tsub; int k; char action = '-'; while (1 > 0) { @@ -33,18 +33,22 @@ namespace Cantera { tok = u; tsize = tok.size(); if (tok[tsize - 1] == '2') { - fctr = m_u[tok.substr(0,tsize-2)]; + tsub = tok.substr(0,tsize-1); + fctr = m_u[tsub]; fctr *= fctr; } else if (tok[tsize - 1] == '3') { - fctr = m_u[tok.substr(0,tsize-2)]; + tsub = tok.substr(0,tsize-1); + fctr = m_u[tsub]; fctr *= fctr*fctr; } - else + else { + tsub = tok; fctr = m_u[tok]; + } if (fctr == 0) - throw CanteraError("toSI","unknown unit: "+tok); + throw CanteraError("toSI","unknown unit: "+tsub); if (action == '-') f *= fctr; else if (action == '/') f /= fctr; if (k < 0) break; diff --git a/data/inputs/elements.xml b/data/inputs/elements.xml index d0139e2b5..e092e826c 100644 --- a/data/inputs/elements.xml +++ b/data/inputs/elements.xml @@ -1,100 +1,101 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/cxx/flame1.cpp b/examples/cxx/flame1.cpp new file mode 100644 index 000000000..e5163aefd --- /dev/null +++ b/examples/cxx/flame1.cpp @@ -0,0 +1,60 @@ + +#include "Cantera.h" +#include "IdealGasMix.h" +#include "transport.h" + +main() { + + // create the gas object + + IdealGasMix gas("gri30.xml"); + doublereal temp = 500.0; + doublereal pres = 2.0*OneAtm; + gas.setState_TPX(temp, pres, "CH4:1.0, O2:2.0, N2:7.52"); + + // create a transport manager that implements + // mixture-averaged transport properties + + Transport* tr = newTransportMgr("Mix", &gas); + + + //============= build each domain ======================== + + + //-------- step 1: create the stagnation flow ------------- + + StFlow flow(&gas); + + // create an initial grid + doublereal z[] = {0.0, 0.05, 0.1, 0.15, 0.2}; + flow.setupGrid(5, z); + + // specify the objects to use to compute kinetic rates and + // transport properties + flow.setKinetics(&gas); + flow.setTransport(&tr); + + flow.setPressure(0.05*OneAtm); + + + + //------- step 2: create the inlet ----------------------- + + Inlet1D inlet; + inlet.setMoleFractions("CH4:1, O2:2, N2:7.52"); + inlet.setMdot(0.1); + + + //------- step 3: create the surface --------------------- + + Surf1D surf; + + //=================== create the container and insert the domains ===== + + vector domains; + domains.push_back(inlet); + domains.push_back(flow); + domains.push_back(surf); + + OneDim flamesim(domains); +} diff --git a/include/Interface.h b/include/Interface.h new file mode 100644 index 000000000..a79dd647b --- /dev/null +++ b/include/Interface.h @@ -0,0 +1,53 @@ +#ifndef CXX_INTERFACE +#define CXX_INTERFACE + +#include + +#include "kernel/SurfPhase.h" +#include "kernel/InterfaceKinetics.h" +#include "kernel/importCTML.h" + +namespace Cantera { + + class Interface : + public SurfPhase, public InterfaceKinetics + { + public: + Interface(string infile, string id, vector phases) + : m_ok(false), m_r(0) { + string path = findInputFile(infile); + ifstream fin(path.c_str()); + if (!fin) { + throw CanteraError("Interface","could not open " + +path+" for reading."); + } + + m_r = new XML_Node("-"); + m_r->build(fin); + + XML_Node* x = find_XML("", m_r, id, "", ""); + if (!x) + throw CanteraError("Interface","error in find_XML"); + + importPhase(*x, this); + phases.push_back(this); + importKinetics(*x, phases, this); + m_ok = true; + } + + + virtual ~Interface() {} + + bool operator!() { return !m_ok;} + bool ready() { return m_ok; } + + protected: + bool m_ok; + XML_Node* m_r; + + private: + }; +} + + +#endif diff --git a/include/onedim.h b/include/onedim.h index b12f8f06b..9dc433542 100755 --- a/include/onedim.h +++ b/include/onedim.h @@ -3,7 +3,7 @@ #include "kernel/oneD/Sim1D.h" #include "kernel/oneD/OneDim.h" -#include "kernel/oneD/Resid1D.h" +#include "kernel/oneD/Domain1D.h" #include "kernel/oneD/Inlet1D.h" #include "kernel/oneD/MultiNewton.h" #include "kernel/oneD/MultiJac.h"