diff --git a/Cantera/clib/src/ctfunc.cpp b/Cantera/clib/src/ctfunc.cpp index 7f1680540..27f94c5ee 100755 --- a/Cantera/clib/src/ctfunc.cpp +++ b/Cantera/clib/src/ctfunc.cpp @@ -59,6 +59,9 @@ extern "C" { "not enough Arrhenius coefficients"); r = new Arrhenius1(n, params); } + else if (type == PeriodicFuncType) { + r = new PeriodicFunc(*_func(n), params[0]); + } else if (type == SumFuncType) { r = new Func1Sum(*_func(n), *_func(m)); } diff --git a/Cantera/clib/src/ctreactor.cpp b/Cantera/clib/src/ctreactor.cpp index b863dc0cc..f6769e8d4 100755 --- a/Cantera/clib/src/ctreactor.cpp +++ b/Cantera/clib/src/ctreactor.cpp @@ -94,10 +94,10 @@ extern "C" { return 0; } - int DLL_EXPORT reactor_setInitialTime(int i, double t) { - _reactor(i)->setInitialTime(t); - return 0; - } + //int DLL_EXPORT reactor_setInitialTime(int i, double t) { + // _reactor(i)->setInitialTime(t); + // return 0; + //} int DLL_EXPORT reactor_setThermoMgr(int i, int n) { _reactor(i)->setThermoMgr(*_th(n)); @@ -111,17 +111,17 @@ extern "C" { return 0; } - int DLL_EXPORT reactor_advance(int i, double t) { - try { - _reactor(i)->advance(t); - return 0; - } - catch (CanteraError) {return -1;} - } + //int DLL_EXPORT reactor_advance(int i, double t) { + // try { + // _reactor(i)->advance(t); + // return 0; + // } + // catch (CanteraError) {return -1;} + //} - double DLL_EXPORT reactor_step(int i, double t) { - return _reactor(i)->step(t); - } + // double DLL_EXPORT reactor_step(int i, double t) { + // return _reactor(i)->step(t); + //} double DLL_EXPORT reactor_time(int i) { return _reactor(i)->time(); @@ -165,48 +165,6 @@ extern "C" { return 0; } -// int DLL_EXPORT reactor_setArea(int i, double a) { -// reactor_t* r = _reactor(i); -// if (r->type() == ReactorType) ((Reactor*)r)->setArea(a); -// return 0; -// } - -// int DLL_EXPORT reactor_setExtTemp(int i, double t) { -// reactor_t* r = _reactor(i); -// if (r->type() == ReactorType) ((Reactor*)r)->setExtTemp(t); -// return 0; -// } - -// int DLL_EXPORT reactor_setExtRadTemp(int i, double t) { -// reactor_t* r = _reactor(i); -// if (r->type() == ReactorType) ((Reactor*)r)->setExtRadTemp(t); -// return 0; -// } - -// int DLL_EXPORT reactor_setVDotCoeff(int i, double v) { -// reactor_t* r = _reactor(i); -// if (r->type() == ReactorType) ((Reactor*)r)->setVDotCoeff(v); -// return 0; -// } - -// int DLL_EXPORT reactor_setHeatTransferCoeff(int i, double h) { -// reactor_t* r = _reactor(i); -// if (r->type() == ReactorType) ((Reactor*)r)->setHeatTransferCoeff(h); -// return 0; -// } - -// int DLL_EXPORT reactor_setEmissivity(int i, double eps) { -// reactor_t* r = _reactor(i); -// if (r->type() == ReactorType) ((Reactor*)r)->setEmissivity(eps); -// return 0; -// } - -// int DLL_EXPORT reactor_setExtPressure(int i, double p) { -// reactor_t* r = _reactor(i); -// if (r->type() == ReactorType) ((Reactor*)r)->setExtPressure(p); -// return 0; -// } - // reactor networks @@ -248,11 +206,16 @@ extern "C" { _reactornet(i)->advance(t); return 0; } - catch (CanteraError) {return -1;} + catch (...) {return -1;} } double DLL_EXPORT reactornet_step(int i, double t) { - return _reactornet(i)->step(t); + try { + return _reactornet(i)->step(t); + } + catch (...) { + return DERR; + } } @@ -263,8 +226,8 @@ extern "C" { switch (type) { case MFC_Type: r = new MassFlowController(); break; - case PressureReg_Type: - r = new PressureRegulator(); break; + case PressureController_Type: + r = new PressureController(); break; case Valve_Type: r = new Valve(); break; default: @@ -278,14 +241,6 @@ extern "C" { return 0; } - int DLL_EXPORT flowdev_copy(int i) { - return Cabinet::cabinet()->newCopy(i); - } - - int DLL_EXPORT flowdev_assign(int i, int j) { - return Cabinet::cabinet()->assign(i,j); - } - int DLL_EXPORT flowdev_install(int i, int n, int m) { try { bool ok = _flowdev(i)->install(*_reactor(n), *_reactor(m) ); @@ -297,26 +252,19 @@ extern "C" { } } - double DLL_EXPORT flowdev_massFlowRate(int i) { - return _flowdev(i)->massFlowRate(); - } - - double DLL_EXPORT flowdev_setpoint(int i) { - return _flowdev(i)->setpoint(); - } - - int DLL_EXPORT flowdev_setSetpoint(int i, double v) { - _flowdev(i)->setSetpoint(v); + int DLL_EXPORT flowdev_setMaster(int i, int n) { + if (_flowdev(i)->type() == PressureController_Type) { + ((PressureController*)_flowdev(i))->setMaster(_flowdev(n)); + } return 0; } - int DLL_EXPORT flowdev_setGains(int i, int n, double* gains) { - _flowdev(i)->setGains(n, gains); - return 0; + double DLL_EXPORT flowdev_massFlowRate(int i, double time) { + return _flowdev(i)->massFlowRate(time); } - int DLL_EXPORT flowdev_getGains(int i, int n, double* gains) { - _flowdev(i)->getGains(n, gains); + int DLL_EXPORT flowdev_setMassFlowRate(int i, double mdot) { + _flowdev(i)->setMassFlowRate(mdot); return 0; } @@ -330,20 +278,6 @@ extern "C" { return 0; } - int DLL_EXPORT flowdev_reset(int i) { - _flowdev(i)->reset(); - return 0; - } - - int DLL_EXPORT flowdev_update(int i) { - _flowdev(i)->update(); - return 0; - } - - double DLL_EXPORT flowdev_maxError(int i) { - return _flowdev(i)->maxError(); - } - int DLL_EXPORT flowdev_ready(int i) { bool ok = _flowdev(i)->ready(); if (ok) return 1; @@ -427,8 +361,13 @@ extern "C" { return 0; } - int DLL_EXPORT wall_setExpansionRate(int i, int n) { - _wall(i)->setExpansionRate(_func(n)); + int DLL_EXPORT wall_setVelocity(int i, int n) { + _wall(i)->setVelocity(_func(n)); + return 0; + } + + int DLL_EXPORT wall_setEmissivity(int i, double epsilon) { + _wall(i)->setEmissivity(epsilon); return 0; } diff --git a/Cantera/clib/src/ctreactor.h b/Cantera/clib/src/ctreactor.h index 92eb5aa1f..8fa53a5e6 100755 --- a/Cantera/clib/src/ctreactor.h +++ b/Cantera/clib/src/ctreactor.h @@ -10,12 +10,12 @@ extern "C" { int DLL_IMPORT reactor_copy(int i); int DLL_IMPORT reactor_assign(int i, int j); int DLL_IMPORT reactor_setInitialVolume(int i, double v); - int DLL_IMPORT reactor_setInitialTime(int i, double t); + // int DLL_IMPORT reactor_setInitialTime(int i, double t); int DLL_IMPORT reactor_setEnergy(int i, int eflag); int DLL_IMPORT reactor_setThermoMgr(int i, int n); int DLL_IMPORT reactor_setKineticsMgr(int i, int n); - int DLL_IMPORT reactor_advance(int i, double t); - double DLL_IMPORT reactor_step(int i, double t); + //int DLL_IMPORT reactor_advance(int i, double t); + //double DLL_IMPORT reactor_step(int i, double t); double DLL_IMPORT reactor_time(int i); double DLL_IMPORT reactor_mass(int i); double DLL_IMPORT reactor_volume(int i); @@ -26,15 +26,6 @@ extern "C" { double DLL_IMPORT reactor_pressure(int i); double DLL_IMPORT reactor_massFraction(int i, int k); - //int DLL_IMPORT reactor_setArea(int i, double a); - //int DLL_IMPORT reactor_setExtTemp(int i, double t); - //int DLL_IMPORT reactor_setExtRadTemp(int i, double t); - //int DLL_IMPORT reactor_setVDotCoeff(int i, double v); - //int DLL_IMPORT reactor_setHeatTransferCoeff(int i, double h); - //int DLL_IMPORT reactor_setEmissivity(int i, double eps); - //int DLL_IMPORT reactor_setExtPressure(int i, double p); - //int DLL_IMPORT reactor_setEnergy(int i, int eflag); - int DLL_IMPORT reactornet_new(); int DLL_IMPORT reactornet_del(int i); int DLL_IMPORT reactornet_copy(int i); @@ -46,19 +37,12 @@ extern "C" { int DLL_IMPORT flowdev_new(int type); int DLL_IMPORT flowdev_del(int i); - //int DLL_IMPORT flowdev_copy(int i); - //int DLL_IMPORT flowdev_assign(int i, int j); int DLL_IMPORT flowdev_install(int i, int n, int m); - double DLL_IMPORT flowdev_massFlowRate(int i); - double DLL_IMPORT flowdev_setpoint(int i); - int DLL_IMPORT flowdev_setSetpoint(int i, double v); - //int DLL_IMPORT flowdev_setGains(int i, int n, double* gains); - //int DLL_IMPORT flowdev_getGains(int i, int n, double* gains); + int DLL_IMPORT flowdev_setMaster(int i, int n); + double DLL_IMPORT flowdev_massFlowRate(int i, double time); + int DLL_IMPORT flowdev_setMassFlowRate(int i, double mdot); int DLL_IMPORT flowdev_setParameters(int i, int n, double* v); int DLL_IMPORT flowdev_setFunction(int i, int n); - //int DLL_IMPORT flowdev_reset(int i); - //int DLL_IMPORT flowdev_update(int i); - //double DLL_IMPORT flowdev_maxError(int i); int DLL_IMPORT flowdev_ready(int i); int DLL_IMPORT wall_new(int type); @@ -75,7 +59,8 @@ extern "C" { int DLL_IMPORT wall_setHeatTransferCoeff(int i, double u); int DLL_IMPORT wall_setHeatFlux(int i, int n); int DLL_IMPORT wall_setExpansionRateCoeff(int i, double k); - int DLL_IMPORT wall_setExpansionRate(int i, int n); + int DLL_IMPORT wall_setVelocity(int i, int n); + int DLL_IMPORT wall_setEmissivity(int i, double epsilon); int DLL_IMPORT wall_ready(int i); } diff --git a/Cantera/matlab/Makefile.in b/Cantera/matlab/Makefile.in index 8e7e676f6..07afebb94 100644 --- a/Cantera/matlab/Makefile.in +++ b/Cantera/matlab/Makefile.in @@ -18,6 +18,7 @@ SRCS = cantera/private/ctmethods.cpp \ cantera/private/kineticsmethods.cpp \ cantera/private/transportmethods.cpp \ cantera/private/reactormethods.cpp \ + cantera/private/reactornetmethods.cpp \ cantera/private/wallmethods.cpp \ cantera/private/flowdevicemethods.cpp \ cantera/private/onedimmethods.cpp \ @@ -60,6 +61,7 @@ install: @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Solution @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@XML_Node/private @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Reactor/private + @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@ReactorNet/private @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Wall/private @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@FlowDevice/private @INSTALL@ -d @prefix@/matlab/toolbox/cantera/cantera/@Func/private @@ -94,6 +96,10 @@ install: @prefix@/matlab/toolbox/cantera/cantera/@Reactor cd cantera/@Reactor/private; @INSTALL@ *.m \ @prefix@/matlab/toolbox/cantera/cantera/@Reactor/private + cd cantera/@ReactorNet; @INSTALL@ *.m \ + @prefix@/matlab/toolbox/cantera/cantera/@ReactorNet + cd cantera/@ReactorNet/private; @INSTALL@ *.m \ + @prefix@/matlab/toolbox/cantera/cantera/@ReactorNet/private cd cantera/@Wall; @INSTALL@ *.m \ @prefix@/matlab/toolbox/cantera/cantera/@Wall cd cantera/@Wall/private; @INSTALL@ *.m \ diff --git a/Cantera/matlab/cantera/private/ctmethods.cpp b/Cantera/matlab/cantera/private/ctmethods.cpp index f4993d176..5a95da858 100644 --- a/Cantera/matlab/cantera/private/ctmethods.cpp +++ b/Cantera/matlab/cantera/private/ctmethods.cpp @@ -24,6 +24,7 @@ const int PHASE_CLASS = 30; const int KINETICS_CLASS = 40; const int TRANSPORT_CLASS = 50; const int REACTOR_CLASS = 60; +const int REACTORNET_CLASS = 65; const int WALL_CLASS = 70; const int FLOWDEVICE_CLASS = 80; const int ONEDIM_CLASS = 90; @@ -54,6 +55,9 @@ void transportmethods( int nlhs, mxArray *plhs[], int nrhs, void reactormethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); +void reactornetmethods( int nlhs, mxArray *plhs[], int nrhs, + const mxArray *prhs[] ); + void wallmethods( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ); @@ -93,6 +97,8 @@ extern "C" { transportmethods(nlhs, plhs, nrhs, prhs); break; case REACTOR_CLASS: reactormethods(nlhs, plhs, nrhs, prhs); break; + case REACTORNET_CLASS: + reactornetmethods(nlhs, plhs, nrhs, prhs); break; case WALL_CLASS: wallmethods(nlhs, plhs, nrhs, prhs); break; case FLOWDEVICE_CLASS: diff --git a/Cantera/matlab/cantera/private/reactormethods.cpp b/Cantera/matlab/cantera/private/reactormethods.cpp index 0e2a4e968..75788696b 100644 --- a/Cantera/matlab/cantera/private/reactormethods.cpp +++ b/Cantera/matlab/cantera/private/reactormethods.cpp @@ -46,18 +46,18 @@ case 4: iok = reactor_setInitialVolume(i, v); break; - case 5: - iok = reactor_setInitialTime(i, v); - break; + // case 5: + //iok = reactor_setInitialTime(i, v); + //break; case 6: iok = reactor_setThermoMgr(i, int(v)); break; case 7: iok = reactor_setKineticsMgr(i, int(v)); break; - case 8: - iok = reactor_advance(i, v); - break; + //case 8: + //iok = reactor_advance(i, v); + //break; case 9: iok = reactor_setEnergy(i, int(v)); break; @@ -76,9 +76,9 @@ else if (job < 40) { switch (job) { - case 21: - r = reactor_step(i, v); - break; + //case 21: + //r = reactor_step(i, v); + //break; case 22: r = reactor_time(i); break; diff --git a/Cantera/python/Cantera/Func.py b/Cantera/python/Cantera/Func.py index efceaf973..07e8ba657 100644 --- a/Cantera/python/Cantera/Func.py +++ b/Cantera/python/Cantera/Func.py @@ -13,7 +13,11 @@ import types class Func1: - """Base class for functions of one variable.""" + """A class for functors of one variable. + + A Functor is an object that behaves like a function. Class 'Func1' + is the base class from which several functor classes derive. These + classes are designed to be used with the Cantera kernel. """ def __init__(self, typ, n, coeffs=[]): self.n = n @@ -69,44 +73,76 @@ class Func1: return RatioFunction(other, self) def func_id(self): + """Return the integer index used internally to access the kernel-level object.""" return self._func_id class Polynomial(Func1): - """A polynomial. The degree is determined by the number of coefficients - supplied. Examples: - - p1 = Polynomial([1.0, -2.0, 3.0]) # 3t^2 - 2t + 1 - p2 = Polynomial([6.0, 8.0]) # 8t + 6 + """A polynomial. + Instances of class 'Polynomial' evaluate + \f[ + f(t) = \sum_{n = 0}^N a_n t^n. + \f] + The coefficients are supplied as a list, beginning with \f$a_N\f$ and ending with \f$a_0\f$. + >>> p1 = Polynomial([1.0, -2.0, 3.0]) # 3t^2 - 2t + 1 + >>> p2 = Polynomial([6.0, 8.0]) # 8t + 6 """ def __init__(self, coeffs=[]): Func1.__init__(self, 2, len(coeffs)-1, coeffs) class Gaussian(Func1): - """A Gaussian pulse. + """A Gaussian pulse. Instances of class 'Gaussian' evaluate + \f[ + f(t) = A \exp[-(t - t_0) / \tau] + \f] + where + \f[ + \tau = \frac{\mbox{FWHM}}{2.0\sqrt{\log(2.0)}} + \f] + Here FWHM denotes the full width at half maximum. """ - def __init__(self, A = 0.0, t0 = 0.0, FWHM = 0.0): - coeffs = array([A, t0, 0.5*FWHM], 'd') + def __init__(self, A = 0.0, t0 = 0.0, FWHM = 1.0): + """ + + A - Peak value. + + t0 - time at which pulse is centered. + + FWHM - full width at half-maximum. + + """ + coeffs = array([A, t0, FWHM], 'd') Func1.__init__(self, 4, 0, coeffs) + class Fourier(Func1): + """Fourier series. Instances of class 'Fourier' evaluate the Fourier series + \f[ + f(t) = \frac{a_0}{2} + \sum_{n=1}^N [a_n \cos(n\omega t) + b_n \sin(n \omega t)] + \f] + where + \f[ + a_n = \int_{-\pi/\omega}^{\pi/\omega} f(t) \cos(n \omega t) dt + \f] + and + \f[ + b_n = \int_{-\pi/\omega}^{\pi/\omega} f(t) \sin(n \omega t) dt. + \f] + The function \f$ f(t) \f$ must be periodic, with period \f$ T = 2\pi/\omega \f$. + >>> coeffs = [(a0, b0), (a1, b1), (a2, b2)] + >>> f = Fourier(omega, coeffs) + Note that b0 must be specified, but is not + used. The value of b0 is arbitrary. """ - Fourier series. - - f(t) = a[0]/2 + sum_{i=1}^n [a[i]*cos(n*omega*t) + b[i]*sin(n*omega*t)] - - Note that b[0] must be specified for symmetry with 'a', but is not - used. - - Example: - - coeffs = [(a0, b0), (a1, b1), (a2, b2)] - f = Fourier(omega, coeffs) - """ - - def __init__(self, omega, c): - cc = asarray(c,'d') + def __init__(self, omega, coefficients): + """ + omega - fundamental frequency [radians/sec]. + + coefficients - List of (a,b) pairs, beginning with \f$n = 0\f$. + + """ + cc = asarray(coefficients,'d') n, m = cc.shape if m <> 2: raise CanteraError('provide (a, b) for each term') @@ -115,17 +151,21 @@ class Fourier(Func1): class Arrhenius(Func1): - """Sum of modified Arrhenius terms. - + """Sum of modified Arrhenius terms. Instances of class 'Arrhenius' evaluate + \f[ f(T) = \sum_{i=1}^n A_n T^{b_n}\exp(-E_n/T) - + \f] Example: - f = Arrhenius([(a0, b0, e0), (a1, b1, e1)]) + >>> f = Arrhenius([(a0, b0, e0), (a1, b1, e1)]) """ - def __init__(self, c): - cc = asarray(c,'d') + def __init__(self, coefficients): + """ + coefficients - sequence of \f$(A, b, E)\f$ triplets. + + """ + cc = asarray(coefficients,'d') n, m = cc.shape if m <> 3: raise CanteraError('Three Arrhenius parameters (A, b, E) required.') @@ -134,35 +174,87 @@ class Arrhenius(Func1): def Const(value): - """Constant function.""" + """Constant function. + >>> f = Const(4.0) # evaluates f(t) = 4.0. + """ return Polynomial([value]) +class PeriodicFunction(Func1): + def __init__(self, func, T): + Func1.__init__(self, 50, func._func_id(), array([T],'d')) + # functions that combine two functions class SumFunction(Func1): - """f = f1 + f2""" + """Sum of two functions. + Instances of class SumFunction evaluate the sum of two supplied functors. + It is not necessary to explicitly create an instance of SumFunction, since + the addition operator of the base class is overloaded to return a SumFunction + instance. + >>> f1 = Polynomial([2.0, 1.0]) + >>> f2 = Polynomial([3.0, -5.0]) + >>> f3 = f1 + f2 # functor to evaluate (2t + 1) + (3t - 5) + In this example, object 'f3' is a functor of class'SumFunction' that calls f1 and f2 + and returns their sum. + """ + def __init__(self, f1, f2): + """ + f1 - first functor. + + f2 - second functor. + """ self.f1 = f1 self.f2 = f2 self.n = -1 - self._func_id = _cantera.func_newcombo(20, f1.func_id(), f2.func_id()) + self._func_id = _cantera.func_newcombo(20, f1._func_id(), f2._func_id()) class ProdFunction(Func1): - """f = f1 * f2""" + """Product of two functions. + Instances of class ProdFunction evaluate the product of two supplied functors. + It is not necessary to explicitly create an instance of 'ProdFunction', since + the multiplication operator of the base class is overloaded to return a 'ProdFunction' + instance. + >>> f1 = Polynomial([2.0, 1.0]) + >>> f2 = Polynomial([3.0, -5.0]) + >>> f3 = f1 * f2 # functor to evaluate (2t + 1)*(3t - 5) + In this example, object 'f3' is a functor of class'ProdFunction' that calls f1 and f2 + and returns their product. + """ def __init__(self, f1, f2): + """ + f1 - first functor. + + f2 - second functor. + """ self.f1 = f1 self.f2 = f2 self.n = -1 - self._func_id = _cantera.func_newcombo(30, f1.func_id(), f2.func_id()) + self._func_id = _cantera.func_newcombo(30, f1._func_id(), f2._func_id()) class RatioFunction(Func1): - """f = f1 / f2""" + """Ratio of two functions. + Instances of class RatioFunction evaluate the ratio of two supplied functors. + It is not necessary to explicitly create an instance of 'RatioFunction', since + the division operator of the base class is overloaded to return a RatioFunction + instance. + >>> f1 = Polynomial([2.0, 1.0]) + >>> f2 = Polynomial([3.0, -5.0]) + >>> f3 = f1 / f2 # functor to evaluate (2t + 1)/(3t - 5) + In this example, object 'f3' is a functor of class'RatioFunction' that calls f1 and f2 + and returns their ratio. + """ def __init__(self, f1, f2): + """ + f1 - first functor. + + f2 - second functor. + """ self.f1 = f1 self.f2 = f2 self.n = -1 - self._func_id = _cantera.func_newcombo(40, f1.func_id(), f2.func_id()) + self._func_id = _cantera.func_newcombo(40, f1._func_id(), f2._func_id()) diff --git a/Cantera/python/Cantera/Phase.py b/Cantera/python/Cantera/Phase.py index 63e23efd3..180031ddd 100755 --- a/Cantera/python/Cantera/Phase.py +++ b/Cantera/python/Cantera/Phase.py @@ -36,32 +36,29 @@ class Phase: def phase_id(self): return self._phase_id -## def addElement(self, symbol, atomicWeight): -## """Add an element.""" -## _cantera.phase_addelement(self._phase_id,symbol,atomicWeight) - -## def addSpecies(self, name, atoms = [], -## thermoType = 'poly', -## thermoCoeffs = []): -## _cantera.phase_addSpecies(self._phase_id, name, array(atoms,'d'), -## thermoType, array(thermoCoeffs,'d'), -## minTemp, maxTemp, refPressure) - def nElements(self): """Number of elements.""" return _cantera.phase_nelements(self._phase_id) - def atomicWeights(self): + def atomicWeights(self, elements = []): """Array of element molar masses [kg/kmol].""" - return _cantera.phase_getarray(self._phase_id,1) - + atw = _cantera.phase_getarray(self._phase_id,1) + if elements: + ae = [] + m = 0 + for e in elements: + m = self.elementIndex(e) + ae.append(atw[m]) + return Numeric.asarray(ae) + else: + return atw + def nSpecies(self): """Number of species.""" return _cantera.phase_nspecies(self._phase_id) def nAtoms(self, species = -1, element = -1): """Number of atoms of element 'element' in species 'species'. - The element and species may be specified by name or by number.""" try: m = self.elementIndex(element) @@ -97,18 +94,20 @@ class Phase: """Mean molar mass [kg/kmol].""" return _cantera.phase_meanmolwt(self._phase_id) - def molecularWeights(self): + def molarMasses(self, species = []): + """Array of species molar masses [kg/kmol].""" + mm = _cantera.phase_getarray(self._phase_id,22) + return self.selectSpecies(mm, species) + + def molecularWeights(self, species = []): """Array of species molar masses [kg/kmol]. DEPRECATED: use molarMasses""" - return _cantera.phase_getarray(self._phase_id,22) + return self.molarMasses(species) - def molarMasses(self): - """Array of species molar masses [kg/kmol].""" - return _cantera.phase_getarray(self._phase_id,22) - - def moleFractions(self): + def moleFractions(self, species = []): """Species mole fraction array.""" - return _cantera.phase_getarray(self._phase_id,20) + x = _cantera.phase_getarray(self._phase_id,20) + return self.selectSpecies(x, species) def moleFraction(self, species=-1): """Mole fraction of a species, referenced by name or @@ -119,12 +118,15 @@ class Phase: k = self.speciesIndex(species) return _cantera.phase_molefraction(self._phase_id,k) - def massFractions(self): + + def massFractions(self, species = []): """Species mass fraction array.""" - return _cantera.phase_getarray(self._phase_id,21) + y = _cantera.phase_getarray(self._phase_id,21) + return self.selectSpecies(y, species) + def massFraction(self, species=-1): - """Mass fraction of one or more species, referenced by name or + """Mass fraction of one species, referenced by name or index number. >>> ph.massFraction(4) >>> ph.massFraction('CH4') @@ -219,7 +221,7 @@ class Phase: if type(x) == types.StringType: _cantera.phase_setstring(self._phase_id,2,x) else: - _cantera.phase_setarray(self._phase_id,2,norm,x) + _cantera.phase_setarray(self._phase_id,2,norm,Numeric.asarray(x)) def setState_TRX(self, t, rho, x): """Set the temperature, density, and mole fractions.""" @@ -230,7 +232,7 @@ class Phase: def setState_TRY(self, t, rho, y): """Set the temperature, density, and mass fractions.""" self.setTemperature(t) - self.setMassFractions(x) + self.setMassFractions(y) self.setDensity(rho) def setState_TR(self, t, rho): @@ -238,3 +240,14 @@ class Phase: self.setTemperature(t) self.setDensity(rho) + def selectSpecies(self, f, sp): + if sp: + fs = [] + k = 0 + for s in sp: + k = self.speciesIndex(s) + fs.append(f[k]) + return Numeric.asarray(fs) + else: + return f + diff --git a/Cantera/python/Cantera/Reactor.py b/Cantera/python/Cantera/Reactor.py index 66a554758..38a9a8589 100644 --- a/Cantera/python/Cantera/Reactor.py +++ b/Cantera/python/Cantera/Reactor.py @@ -20,14 +20,23 @@ class ReactorBase: 2 = Reservoir). """ self.__reactor_id = _cantera.reactor_new(type) + self._type = type self._inlets = [] self._outlets = [] self._walls = [] + self._reservoirs = [] self._name = name self._verbose = verbose self.insert(contents) self._setInitialVolume(volume) self._setEnergy(energy) + if self._verbose: + print 'Created '+self._name + print ' Volume = ',volume,' m^3' + if energy <> 'on': + print ' Temperature will be held constant' + print ' Initial State:' + print contents def __del__(self): @@ -38,12 +47,14 @@ class ReactorBase: def __str__(self): s = self._name + s += ':\n Volume = '+`self.volume()` if self._contents: - s += ": \n"+`self._contents` + s += "\n"+`self._contents` return s def __repr__(self): s = self._name + s += ':\n Volume = '+`self.volume()` if self._contents: s += ": \n"+`self._contents` return s @@ -206,27 +217,33 @@ class ReactorBase: """ return self._walls - def _addInlet(self, inlet): + def _addInlet(self, inlet, other): """For internal use. Store a reference to 'inlet' so that it will not be deleted before this object.""" self._inlets.append(inlet) + if self._type == 1 and other._type == 2: + self._reservoirs.append(other) - def _addOutlet(self, outlet): + def _addOutlet(self, outlet, other): """For internal use. Store a reference to 'outlet' so that it will not be deleted before this object.""" self._outlets.append(outlet) + if self._type == 1 and other._type == 2: + self._reservoirs.append(other) - def _addWall(self, wall): + def _addWall(self, wall, other): """For internal use. Store a reference to 'wall' so that it will not be deleted before this object.""" self._walls.append(wall) + if self._type == 1 and other._type == 2: + self._reservoirs.append(other) def syncContents(self): """Set the state of the object representing the reactor contents to the current reactor state. >>> r = Reactor(gas) >>> (statements that change the state of object 'gas') - >>> r.syncContents(self) + >>> r.syncContents() After this statement, the state of object 'gas' is synchronized with the reactor state. See 'contents'. @@ -356,6 +373,7 @@ class Reservoir(ReactorBase): + #------------------ FlowDevice --------------------------------- class FlowDevice: @@ -388,21 +406,9 @@ class FlowDevice: """ return _cantera.flowdev_ready(self.__fdev_id) - def massFlowRate(self): + def massFlowRate(self, time = -999.0): """Mass flow rate (kg/s). """ - return _cantera.flowdev_massFlowRate(self.__fdev_id) - - def setSetpoint(self, v): - """ - Deprecated. Set the set point. - """ - _cantera.flowdev_setSetpoint(self.__fdev_id, v) - - def setpoint(self): - """ - Deprecated. The setpoint value. - """ - return _cantera.flowdev_setpoint(self.__fdev_id) + return _cantera.flowdev_massFlowRate(self.__fdev_id, time) def install(self, upstream, downstream): """ @@ -413,8 +419,8 @@ class FlowDevice: if self._verbose: print print self._name+': installing between '+upstream.name()+' and '+downstream.name() - upstream._addOutlet(self) - downstream._addInlet(self) + upstream._addOutlet(self, downstream) + downstream._addInlet(self, upstream) _cantera.flowdev_install(self.__fdev_id, upstream.reactor_id(), downstream.reactor_id()) def _setParameters(self, c): @@ -422,7 +428,12 @@ class FlowDevice: n = len(params) return _cantera.flowdev_setParameters(self.__fdev_id, n, params) + def setFunction(self, f): + _cantera.flowdev_setFunction(self.__fdev_id, f.func_id()) + def flowdev_id(self): + return self.__fdev_id + _mfccount = 0 class MassFlowController(FlowDevice): @@ -487,14 +498,18 @@ class MassFlowController(FlowDevice): """ if self._verbose: print self._name+': setting mdot to '+`mdot`+' kg/s' - self.setSetpoint(mdot) + if type(mdot) == types.InstanceType: + self.setFunction(mdot) + else: + _cantera.flowdev_setMassFlowRate(self.flowdev_id(), mdot) + def set(self, mdot = 0.0): """Set the mass flow rate [kg/s]. >>> mfc.set(mdot = 0.2) """ - self.setSetpoint(mdot) + self._setMassFlowRate(mdot) _valvecount = 0 @@ -535,7 +550,7 @@ class Valve(FlowDevice): """ def __init__(self, upstream=None, downstream=None, - name='', Kv = 0.0, verbose=0): + name='', Kv = 0.0, mdot0 = 0.0, verbose=0): """ upstream - upstream reactor or reservoir. @@ -559,19 +574,89 @@ class Valve(FlowDevice): FlowDevice.__init__(self,3,name,verbose) if upstream and downstream: self.install(upstream, downstream) - self.setValveCoeff(Kv) + self.setValveCoeff(Kv, mdot0) - def setValveCoeff(self, v): + def setValveCoeff(self, Kv, mdot0 = 0.0): """Set or reset the valve coefficient \f$ K_v \f$.""" - vv = zeros(1,'d') - vv[0] = v + vv = zeros(2,'d') + vv[0] = Kv + vv[1] = mdot0 if self._verbose: print - print self._name+': setting valve coefficient to '+`v`+' kg/Pa-s' + print self._name+': setting valve coefficient to '+`Kv`+' kg/Pa-s' self._setParameters(vv) + def _setValveCharacteristic(self, f): + """Set or reset the valve characteristics. + """ + if type(f) == types.InstanceType: + self.setFunction(f) + else: + raise CanteraError("Wrong type for valve characteristic function.") + def set(self, Kv = -1.0, mdot = 0.0, F = None): + if F: + self.setFunction(F) + if Kv > 0.0: + self.setValveCoeff(Kv, mdot0 = mdot) + + + +_pccount = 0 + +class PressureController(FlowDevice): + + def __init__(self, upstream=None, downstream=None, + name='', master = None, Kv = 0.0, verbose=0): + """ + upstream - upstream reactor or reservoir. + + downstream - downstream reactor or reservoir. + + name - name used to identify the valve in output. + If no name is specified, it defaults to 'Valve_n', where n is an + integer assigned in the order the Valve object + was created. + + Kv - the constant in the mass flow rate equation. + + verbose - if set to a positive integer, additional diagnostic + information will be printed. + + """ + global _pccount + if name == '': + name = 'PressureController_'+`_pccount` + _pccount += 1 + FlowDevice.__init__(self,2,name,verbose) + if upstream and downstream: + self.install(upstream, downstream) + self.setPressureCoeff(Kv) + self.setMaster(master) + + + def setPressureCoeff(self, Kv): + """Set or reset the pressure coefficient \f$ K_v \f$.""" + vv = zeros(1,'d') + vv[0] = Kv + if self._verbose: + print + print self._name+': setting pressure coefficient to '+`Kv`+' kg/Pa-s' + self._setParameters(vv) + + def setMaster(self, master): + _cantera.flowdev_setMaster(self.flowdev_id(), + master.flowdev_id()) + + def set(self, Kv = -1.0, master = None): + if master: + self.setMaster(master) + if Kv > 0.0: + self.setPressureCoeff(Kv) + + + #------------- Wall --------------------------- _wallcount = 0 @@ -583,7 +668,7 @@ class Wall: """ def __init__(self, left=None, right=None, name = '', A = 1.0, K = 0.0, U = 0.0, - Q = None, Vdot = None, + Q = None, velocity = None, kinetics = [None, None]): typ = 0 self.__wall_id = _cantera.wall_new(typ) @@ -601,7 +686,7 @@ class Wall: raise CanteraError('both left and right reactors must be specified.') self.setArea(A) self.setExpansionRateCoeff(K) - self.setExpansionRate(Vdot) + self.setVelocity(velocity) self.setHeatTransferCoeff(U) self.setHeatFlux(Q) @@ -641,6 +726,13 @@ class Wall: """ return _cantera.wall_setHeatTransferCoeff(self.__wall_id, u) + def setEmissivity(self, epsilon): + """ + Set the emissivity. + The radiative heat flux through the wall is computed from + \f[ q_r = \epsion \sigma (T_\ell^4 - T_r^4) \f] + """ + def setHeatFlux(self, qfunc=None): """ Specify the time-dependent heat flux function [W/m2]. @@ -655,13 +747,13 @@ class Wall: resulting from a unit pressure drop.""" _cantera.wall_setExpansionRateCoeff(self.__wall_id, k) - def setExpansionRate(self, vfunc=None): + def setVelocity(self, vfunc=None): """ - Specify the volumetric expansion rate function [m^3/s]. + Specify the velocity function [m/s]. """ n = 0 if vfunc: n = vfunc.func_id() - _cantera.wall_setExpansionRate(self.__wall_id, n) + _cantera.wall_setVelocity(self.__wall_id, n) def vdot(self): """Rate of volume change [m^3]. A positive value corresponds @@ -676,8 +768,8 @@ class Wall: return _cantera.wall_Q(self.__wall_id) def install(self, left, right): - left._addWall(self) - right._addWall(self) + left._addWall(self, right) + right._addWall(self, left) _cantera.wall_install(self.__wall_id, left.reactor_id(), right.reactor_id()) diff --git a/Cantera/python/Cantera/ThermoPhase.py b/Cantera/python/Cantera/ThermoPhase.py index d922af7ca..f75526194 100644 --- a/Cantera/python/Cantera/ThermoPhase.py +++ b/Cantera/python/Cantera/ThermoPhase.py @@ -106,48 +106,53 @@ class ThermoPhase(Phase): """ The pressure [Pa].""" return _cantera.thermo_getfp(self._phase_id,7) - def chemPotentials(self): + def chemPotentials(self, species = []): """Species chemical potentials. This method returns an array containing the species chemical potentials [J/kmol]. The expressions used to compute these depend on the model implemented by the underlying kernel thermo manager.""" - return _cantera.thermo_getarray(self._phase_id,20) + mu = _cantera.thermo_getarray(self._phase_id,20) + return self.selectSpecies(mu, species) - def enthalpies_RT(self): + def enthalpies_RT(self, species = []): """Pure species non-dimensional enthalpies. This method returns an array containing the pure-species standard-state enthalpies divided by RT. For gaseous species, these values are ideal gas enthalpies.""" - return _cantera.thermo_getarray(self._phase_id,23) - - def entropies_R(self): + hrt = _cantera.thermo_getarray(self._phase_id,23) + return self.selectSpecies(hrt, species) + + def entropies_R(self, species = []): """Pure species non-dimensional entropies. This method returns an array containing the pure-species standard-state entropies divided by R. For gaseous species, these values are ideal gas entropies.""" - return _cantera.thermo_getarray(self._phase_id,24) + sr = _cantera.thermo_getarray(self._phase_id,24) + return self.selectSpecies(sr, species) - def gibbs_RT(self): + def gibbs_RT(self, species = []): """Pure species non-dimensional Gibbs free energies. This method returns an array containing the pure-species standard-state Gibbs free energies divided by R. For gaseous species, these are ideal gas values.""" - return (_cantera.thermo_getarray(self._phase_id,23) + grt = (_cantera.thermo_getarray(self._phase_id,23) - _cantera.thermo_getarray(self._phase_id,24)) + return self.selectSpecies(grt, species) - def cp_R(self): + def cp_R(self, species = []): """Pure species non-dimensional heat capacities at constant pressure. This method returns an array containing the pure-species standard-state heat capacities divided by R. For gaseous species, these values are ideal gas heat capacities.""" - return _cantera.thermo_getarray(self._phase_id,25) + cpr = _cantera.thermo_getarray(self._phase_id,25) + return self.selectSpecies(cpr, species) def setPressure(self, p): diff --git a/Cantera/python/src/ctreactor_methods.cpp b/Cantera/python/src/ctreactor_methods.cpp index d0b728746..36da789f9 100644 --- a/Cantera/python/src/ctreactor_methods.cpp +++ b/Cantera/python/src/ctreactor_methods.cpp @@ -32,17 +32,17 @@ py_reactor_setInitialVolume(PyObject *self, PyObject *args) return Py_BuildValue("i",0); } -static PyObject* -py_reactor_setInitialTime(PyObject *self, PyObject *args) -{ - int n; - double t; - if (!PyArg_ParseTuple(args, "id:reactor_setInitialTime", &n, &t)) - return NULL; - int iok = reactor_setInitialTime(n, t); - if (iok < 0) return reportError(iok); - return Py_BuildValue("i",0); -} +// static PyObject* +// py_reactor_setInitialTime(PyObject *self, PyObject *args) +// { +// int n; +// double t; +// if (!PyArg_ParseTuple(args, "id:reactor_setInitialTime", &n, &t)) +// return NULL; +// int iok = reactor_setInitialTime(n, t); +// if (iok < 0) return reportError(iok); +// return Py_BuildValue("i",0); +// } static PyObject* py_reactor_setEnergy(PyObject *self, PyObject *args) @@ -79,27 +79,27 @@ py_reactor_setKineticsMgr(PyObject *self, PyObject *args) return Py_BuildValue("i",0); } -static PyObject* -py_reactor_advance(PyObject *self, PyObject *args) -{ - int n; - double t; - if (!PyArg_ParseTuple(args, "id:reactor_advance", &n, &t)) - return NULL; - int iok = reactor_advance(n, t); - if (iok < 0) return reportError(iok); - return Py_BuildValue("i",0); -} +// static PyObject* +// py_reactor_advance(PyObject *self, PyObject *args) +// { +// int n; +// double t; +// if (!PyArg_ParseTuple(args, "id:reactor_advance", &n, &t)) +// return NULL; +// int iok = reactor_advance(n, t); +// if (iok < 0) return reportError(iok); +// return Py_BuildValue("i",0); +// } -static PyObject* -py_reactor_step(PyObject *self, PyObject *args) -{ - int n; - double t; - if (!PyArg_ParseTuple(args, "id:reactor_step", &n, &t)) - return NULL; - return Py_BuildValue("d",reactor_step(n, t)); -} +// static PyObject* +// py_reactor_step(PyObject *self, PyObject *args) +// { +// int n; +// double t; +// if (!PyArg_ParseTuple(args, "id:reactor_step", &n, &t)) +// return NULL; +// return Py_BuildValue("d",reactor_step(n, t)); +// } static PyObject* py_reactor_time(PyObject *self, PyObject *args) @@ -224,34 +224,46 @@ py_flowdev_install(PyObject *self, PyObject *args) return Py_BuildValue("i",0); } +static PyObject* +py_flowdev_setMaster(PyObject *self, PyObject *args) +{ + int n, m; + if (!PyArg_ParseTuple(args, "ii:flowdev_setMaster", &n, &m)) + return NULL; + int iok = flowdev_setMaster(n, m); + if (iok < 0) return reportError(iok); + return Py_BuildValue("i",0); +} + static PyObject* py_flowdev_massFlowRate(PyObject *self, PyObject *args) { int n; - if (!PyArg_ParseTuple(args, "i:flowdev_massFlowRate", &n)) + double t; + if (!PyArg_ParseTuple(args, "id:flowdev_massFlowRate", &n, &t)) return NULL; - double mdot = flowdev_massFlowRate(n); + double mdot = flowdev_massFlowRate(n, t); return Py_BuildValue("d",mdot); } -static PyObject* -py_flowdev_setpoint(PyObject *self, PyObject *args) -{ - int n; - if (!PyArg_ParseTuple(args, "i:flowdev_setpoint", &n)) - return NULL; - double v = flowdev_setpoint(n); - return Py_BuildValue("d",v); -} +// static PyObject* +// py_flowdev_setpoint(PyObject *self, PyObject *args) +// { +// int n; +// if (!PyArg_ParseTuple(args, "i:flowdev_setpoint", &n)) +// return NULL; +// double v = flowdev_setpoint(n); +// return Py_BuildValue("d",v); +// } static PyObject* -py_flowdev_setSetpoint(PyObject *self, PyObject *args) +py_flowdev_setMassFlowRate(PyObject *self, PyObject *args) { int n; - double v; - if (!PyArg_ParseTuple(args, "id:flowdev_setSetpoint", &n, &v)) + double mdot; + if (!PyArg_ParseTuple(args, "id:flowdev_setMassFlowRate", &n, &mdot)) return NULL; - int iok = flowdev_setSetpoint(n, v); + int iok = flowdev_setMassFlowRate(n, mdot); if (iok < 0) return reportError(iok); return Py_BuildValue("i",0); } @@ -404,6 +416,18 @@ py_wall_setHeatTransferCoeff(PyObject *self, PyObject *args) return Py_BuildValue("i",0); } +static PyObject* +py_wall_setEmissivity(PyObject *self, PyObject *args) +{ + int n; + double epsilon; + if (!PyArg_ParseTuple(args, "id:wall_setEmissivity", &n, &epsilon)) + return NULL; + int iok = wall_setEmissivity(n,epsilon); + if (iok < 0) return reportError(iok); + return Py_BuildValue("i",0); +} + static PyObject* py_wall_setExpansionRateCoeff(PyObject *self, PyObject *args) { @@ -417,12 +441,12 @@ py_wall_setExpansionRateCoeff(PyObject *self, PyObject *args) } static PyObject* -py_wall_setExpansionRate(PyObject *self, PyObject *args) +py_wall_setVelocity(PyObject *self, PyObject *args) { int n, m; - if (!PyArg_ParseTuple(args, "ii:wall_setExpansionRate", &n, &m)) + if (!PyArg_ParseTuple(args, "ii:wall_setVelocity", &n, &m)) return NULL; - int iok = wall_setExpansionRate(n,m); + int iok = wall_setVelocity(n,m); if (iok < 0) return reportError(iok); return Py_BuildValue("i",0); } diff --git a/Cantera/python/src/methods.h b/Cantera/python/src/methods.h index 80df15188..b17bc4c53 100644 --- a/Cantera/python/src/methods.h +++ b/Cantera/python/src/methods.h @@ -192,21 +192,22 @@ static PyMethodDef ct_methods[] = { {"bndry_setmdot", py_bndry_setmdot, METH_VARARGS}, {"flowdev_ready", py_flowdev_ready, METH_VARARGS}, - {"reactor_setInitialTime", py_reactor_setInitialTime, METH_VARARGS}, + //{"reactor_setInitialTime", py_reactor_setInitialTime, METH_VARARGS}, {"reactornet_setInitialTime", py_reactornet_setInitialTime, METH_VARARGS}, {"flowdev_new", py_flowdev_new, METH_VARARGS}, {"flowdev_massFlowRate", py_flowdev_massFlowRate, METH_VARARGS}, {"flowdev_del", py_flowdev_del, METH_VARARGS}, - {"flowdev_setpoint", py_flowdev_setpoint, METH_VARARGS}, + // {"flowdev_setpoint", py_flowdev_setpoint, METH_VARARGS}, {"reactor_temperature", py_reactor_temperature, METH_VARARGS}, - {"flowdev_setSetpoint", py_flowdev_setSetpoint, METH_VARARGS}, + {"flowdev_setMassFlowRate", py_flowdev_setMassFlowRate, METH_VARARGS}, {"flowdev_install", py_flowdev_install, METH_VARARGS}, + {"flowdev_setMaster", py_flowdev_setMaster, METH_VARARGS}, {"reactor_setThermoMgr", py_reactor_setThermoMgr, METH_VARARGS}, {"reactor_setEnergy", py_reactor_setEnergy, METH_VARARGS}, {"reactor_volume", py_reactor_volume, METH_VARARGS}, {"reactor_time", py_reactor_time, METH_VARARGS}, - {"reactor_advance", py_reactor_advance, METH_VARARGS}, - {"reactor_step", py_reactor_step, METH_VARARGS}, + // {"reactor_advance", py_reactor_advance, METH_VARARGS}, + //{"reactor_step", py_reactor_step, METH_VARARGS}, {"reactornet_addreactor", py_reactornet_addreactor, METH_VARARGS}, {"reactornet_advance", py_reactornet_advance, METH_VARARGS}, {"reactornet_step", py_reactornet_step, METH_VARARGS}, @@ -235,7 +236,8 @@ static PyMethodDef ct_methods[] = { {"wall_new", py_wall_new, METH_VARARGS}, {"wall_vdot", py_wall_vdot, METH_VARARGS}, {"wall_del", py_wall_del, METH_VARARGS}, - {"wall_setExpansionRate", py_wall_setExpansionRate, METH_VARARGS}, + {"wall_setVelocity", py_wall_setVelocity, METH_VARARGS}, + {"wall_setEmissivity", py_wall_setEmissivity, METH_VARARGS}, {"wall_setExpansionRateCoeff", py_wall_setExpansionRateCoeff, METH_VARARGS}, {"wall_ready", py_wall_ready, METH_VARARGS}, diff --git a/Cantera/src/Func1.h b/Cantera/src/Func1.h index 93b32e496..bf0b7b125 100644 --- a/Cantera/src/Func1.h +++ b/Cantera/src/Func1.h @@ -28,7 +28,8 @@ namespace Cantera { const int DiffFuncType = 25; const int ProdFuncType = 30; const int RatioFuncType = 40; - + const int PeriodicFuncType = 50; + /** * Base class for 'functor' classes that evaluate a function of * one variable. @@ -47,10 +48,10 @@ namespace Cantera { class Gaussian : public Func1 { public: - Gaussian(double A, double t0, double tau) { + Gaussian(double A, double t0, double fwhm) { m_A = A; m_t0 = t0; - m_tau = tau; + m_tau = fwhm/(2.0*sqrt(log(2.0))); } virtual ~Gaussian() {} virtual doublereal eval(doublereal t) { @@ -177,6 +178,29 @@ namespace Cantera { vector_fp m_A, m_b, m_E; }; + /** + * Periodic function. Takes any function and makes it + * periodic with period T. + */ + class PeriodicFunc : public Func1 { + public: + PeriodicFunc(Func1& f, doublereal T) { + m_func = &f; + m_period = T; + } + virtual ~PeriodicFunc() {} + virtual doublereal eval(doublereal t) { + int np = int(t/m_period); + doublereal time = t - np*m_period; + return m_func->eval(time); + } + protected: + Func1* m_func; + doublereal m_period; + + private: + }; + /** * Sum of two functions. diff --git a/Cantera/src/IdealGasPhase.h b/Cantera/src/IdealGasPhase.h index 3df19b1b0..f00df4451 100644 --- a/Cantera/src/IdealGasPhase.h +++ b/Cantera/src/IdealGasPhase.h @@ -228,6 +228,14 @@ namespace Cantera { } + virtual doublereal isothermalCompressibility() { + return -1.0/pressure(); + } + + virtual doublereal thermalExpansionCoeff() { + return 1.0/temperature(); + } + // new methods defined here const array_fp& enthalpy_RT() const { diff --git a/Cantera/src/NasaThermo.h b/Cantera/src/NasaThermo.h index d706c7fbe..b1da5b626 100755 --- a/Cantera/src/NasaThermo.h +++ b/Cantera/src/NasaThermo.h @@ -86,11 +86,7 @@ namespace Cantera { vector_fp chigh(7); copy(c + 8, c + 15, chigh.begin()); - // Make cp exactly continuous at tmid by offsetting the - // high temp fit. - doublereal cplow = poly4(tmid, clow+2); - doublereal cphigh = poly4(tmid, chigh.begin()+2); - chigh[0] += (cplow - cphigh); + checkContinuity(tmid, clow, chigh.begin()); m_high[igrp-1].push_back(NasaPoly1(index, tmid, thigh, pref, chigh.begin())); @@ -105,7 +101,6 @@ namespace Cantera { m_low_map[index] = &m_low[igrp-1].back(); } - /** * update the properties for only one species. */ @@ -193,6 +188,52 @@ namespace Cantera { doublereal m_p0; int m_ngroups; mutable vector_fp m_t; + + private: + + + /// Check the continuity of properties at the midpoint + /// temperature, and adjust the high-T coefficients to + /// make the properties exactly continuous at Tmid. + void checkContinuity(double tmid, const doublereal* clow, + doublereal* chigh) { + + // heat capacity + doublereal cplow = poly4(tmid, clow+2); + doublereal cphigh = poly4(tmid, chigh+2); + doublereal delta = cplow - cphigh; + if (fabs(delta/cplow) > 0.001) { + writelog("WARNING: discontinuity in cp/R detected at Tmid = " + +fp2str(tmid)+"\n"); + writelog(" Adjusting high-temperature coefficient to fix.\n"); + } + chigh[2] += cplow - cphigh; + + // enthalpy + doublereal hrtlow = enthalpy_RT(tmid, clow); + doublereal hrthigh = enthalpy_RT(tmid, chigh); + chigh[0] += tmid*(hrtlow - hrthigh); + + // entropy + doublereal srlow = entropy_R(tmid, clow); + doublereal srhigh = entropy_R(tmid, chigh); + chigh[1] += srlow - srhigh; + } + + /// for internal use by checkContinuity + doublereal enthalpy_RT(double t, const doublereal* c) { + return c[2] + 0.5*c[3]*t + OneThird*c[4]*t*t + + 0.25*c[5]*t*t*t + 0.2*c[6]*t*t*t*t + + c[0]/t; + } + + /// for internal use by checkContinuity + doublereal entropy_R(double t, const doublereal* c) { + return c[2]*log(t) + c[3]*t + 0.5*c[4]*t*t + + OneThird*c[5]*t*t*t + 0.25*c[6]*t*t*t*t + + c[1]; + } + }; } @@ -200,7 +241,10 @@ namespace Cantera { #endif // $Log$ -// Revision 1.2 2003-11-01 04:50:35 dggoodwin +// Revision 1.3 2004-04-22 21:44:36 dggoodwin +// *** empty log message *** +// +// Revision 1.2 2003/11/01 04:50:35 dggoodwin // *** empty log message *** // // Revision 1.1.1.1 2003/04/14 17:57:51 dggoodwin diff --git a/Cantera/src/PureFluidPhase.h b/Cantera/src/PureFluidPhase.h index b18d7bff8..9b17965a3 100644 --- a/Cantera/src/PureFluidPhase.h +++ b/Cantera/src/PureFluidPhase.h @@ -149,6 +149,13 @@ namespace Cantera { mu[0] = gibbs_mole(); } + virtual doublereal isothermalCompressibility() { + return m_sub->isothermalCompressibility(); + } + + virtual doublereal thermalExpansionCoeff() { + return m_sub->thermalExpansionCoeff(); + } tpx::Substance& TPX_Substance() { return *m_sub; } diff --git a/Cantera/src/ThermoPhase.cpp b/Cantera/src/ThermoPhase.cpp index 4e2071343..0292ff50d 100644 --- a/Cantera/src/ThermoPhase.cpp +++ b/Cantera/src/ThermoPhase.cpp @@ -115,13 +115,15 @@ namespace Cantera { for (int n = 0; n < 50; n++) { dt = (u - intEnergy_mass())/cv_mass(); if (dt > 100.0) dt = 100.0; - else if (dt < -100.0) dt = -100.0; + else if (dt < -100.0) dt = -100.0; setTemperature(temperature() + dt); if (fabs(dt) < tol) { return; } } - throw CanteraError("setState_UV","no convergence. dt = " + fp2str(dt)); + throw CanteraError("setState_UV", + "no convergence. dt = " + fp2str(dt)+"\n" + +"u = "+fp2str(u)+" v = "+fp2str(v)+"\n"); } void ThermoPhase::setState_SP(doublereal s, doublereal p, diff --git a/Cantera/src/ThermoPhase.h b/Cantera/src/ThermoPhase.h index 63afe557e..2848b4690 100755 --- a/Cantera/src/ThermoPhase.h +++ b/Cantera/src/ThermoPhase.h @@ -569,7 +569,14 @@ namespace Cantera { */ virtual void setParameters(int n, doublereal* c) {} - + virtual doublereal isothermalCompressibility() { + err("isothermalCompressibility"); return -1.0; + } + + virtual doublereal thermalExpansionCoeff() { + err("thermalExpansionCoeff()"); return -1.0; + } + //--------------------------------------------------------- /// @name Critical state properties. diff --git a/Cantera/src/ctexceptions.h b/Cantera/src/ctexceptions.h index 23d985e28..16217fbee 100755 --- a/Cantera/src/ctexceptions.h +++ b/Cantera/src/ctexceptions.h @@ -24,6 +24,7 @@ namespace Cantera { CanteraError() {} CanteraError(string proc, string msg) { setError(proc, msg); + writelog("Throwing CanteraError. "+proc+" "+msg+"\n"); //m_msg = msg; } virtual ~CanteraError(){} diff --git a/Cantera/src/zeroD/FlowDevice.cpp b/Cantera/src/zeroD/FlowDevice.cpp index 1a6eaafba..b23ffcb6d 100644 --- a/Cantera/src/zeroD/FlowDevice.cpp +++ b/Cantera/src/zeroD/FlowDevice.cpp @@ -34,14 +34,16 @@ namespace Cantera { return true; } - void FlowDevice::setFunction(Func1* f) {} + void FlowDevice::setFunction(Func1* f) { + m_func = f; + } /** * Mass flow rate of outlet species k. Returns zero if this * species is not present in the upstream mixture. */ - doublereal FlowDevice::massFlowRate(int k) { + doublereal FlowDevice::outletSpeciesMassFlowRate(int k) { if (k < 0 || k >= m_nspout) return 0.0; int ki = m_out2in[k]; if (ki < 0) return 0.0; diff --git a/Cantera/src/zeroD/FlowDevice.h b/Cantera/src/zeroD/FlowDevice.h index a0ff54bd5..7baddbdaa 100644 --- a/Cantera/src/zeroD/FlowDevice.h +++ b/Cantera/src/zeroD/FlowDevice.h @@ -24,7 +24,7 @@ namespace Cantera { class Func1; const int MFC_Type = 1; - const int PressureReg_Type = 2; + const int PressureController_Type = 2; const int Valve_Type = 3; /** @@ -43,79 +43,92 @@ namespace Cantera { public: /// Constructor - FlowDevice() : m_mdot(0.0), + FlowDevice() : m_mdot(0.0), m_func(0), m_type(0), m_nspin(0), m_nspout(0), m_in(0), m_out(0) {} /// Destructor (does nothing) virtual ~FlowDevice(){} - /// Copy constructor. - FlowDevice(const FlowDevice& a) : m_in(a.m_in), m_out(a.m_out) {} + // /// Copy constructor. +// FlowDevice(const FlowDevice& a) : m_in(a.m_in), m_out(a.m_out) {} - /// Assignment operator - FlowDevice& operator=(const FlowDevice& a) { - if (this == &a) return *this; - m_in = a.m_in; - m_out = a.m_out; - return *this; - } - - /** - * Mass flow rate (kg/s). May be overloaded in derived - * classes. - */ - virtual doublereal massFlowRate() {return m_mdot;} - - doublereal massFlowRate(int k); - virtual doublereal enthalpy_mass(); +// /// Assignment operator +// FlowDevice& operator=(const FlowDevice& a) { +// if (this == &a) return *this; +// m_in = a.m_in; +// m_out = a.m_out; +// return *this; +// } - /** - * Setpoint. Default = 0.0. - */ - virtual doublereal setpoint() { warn("setpoint"); return 0.0; } - - /* Update the internal state, if necessary. By default this method - * does nothing, but may be overloaded for devices that have a - * state. - */ - virtual void update() {warn("update");} - - /* Reset the device. By default this method does nothing, but - * may be overloaded for devices that have a state that depends on - * past history. - */ - virtual void reset() {warn("reset");} + int type() { return m_type; } /** - * Set the setpoint. May be changed at any time. By default, - * this does nothing. + * Mass flow rate (kg/s). */ - virtual void setSetpoint(doublereal value) {warn("setSetpoint");} - - /** - * Set the controller gains. Returns false if the number of - * gains is too small, or if an illegal value is specified. - */ - virtual bool setGains(int n, const doublereal* gains) { - warn("setGains"); - return true; + doublereal massFlowRate(double time = -999.0) { + if (time != -999.0) updateMassFlowRate(time); + return m_mdot; } - /** - * Get the controller gains. Returns false if the 'gains' - * array is too small. - */ - virtual bool getGains(int n, doublereal* gains) { - warn("getGains"); - return true; - } + // Update the mass flow rate at time 'time'. This must be + // overloaded in subclassess to update m_mdot. + virtual void updateMassFlowRate(doublereal time) {} - /** - * Maximum difference between input and setpoint since - * last call to 'reset'. - */ - virtual doublereal maxError() {warn("maxError"); return 0.0;} + // mass flow rate of outlet species k + doublereal outletSpeciesMassFlowRate(int k); + + // specific enthalpy + doublereal enthalpy_mass(); + +// /** +// * Setpoint. Default = 0.0. +// */ +// virtual doublereal setpoint() { warn("setpoint"); return 0.0; } + + +// /* Update the internal state, if necessary. By default this method +// * does nothing, but may be overloaded for devices that have a +// * state. +// */ +// virtual void update() {warn("update");} + + +// /* Reset the device. By default this method does nothing, but +// * may be overloaded for devices that have a state that depends on +// * past history. +// */ +// virtual void reset() {warn("reset");} + +// /** +// * Set the setpoint. May be changed at any time. By default, +// * this does nothing. +// */ +// virtual void setSetpoint(doublereal value) {warn("setSetpoint");} + +// /** +// * Set the controller gains. Returns false if the number of +// * gains is too small, or if an illegal value is specified. +// */ +// virtual bool setGains(int n, const doublereal* gains) { +// warn("setGains"); +// return true; +// } + +// /** +// * Get the controller gains. Returns false if the 'gains' +// * array is too small. +// */ +// virtual bool getGains(int n, doublereal* gains) { +// warn("getGains"); +// return true; +// } + +// /** +// * Maximum difference between input and setpoint since +// * last call to 'reset'. +// */ +// virtual doublereal maxError() {warn("maxError"); return 0.0;} /** * Install a flow device between two reactors. @@ -125,7 +138,6 @@ namespace Cantera { bool install(ReactorBase& in, ReactorBase& out); virtual bool ready() { return (m_in != 0 && m_out != 0); } - doublereal m_mdot; /// Return a reference to the upstream reactor. ReactorBase& in() const { return *m_in; } @@ -139,11 +151,16 @@ namespace Cantera { copy(coeffs, coeffs + n, m_coeffs.begin()); } - virtual void setFunction(Func1* f); + void setFunction(Func1* f); + void setMassFlowRate(doublereal mdot) {m_mdot = mdot;} + protected: + doublereal m_mdot; + Func1* m_func; vector_fp m_coeffs; + int m_type; private: diff --git a/Cantera/src/zeroD/Reactor.cpp b/Cantera/src/zeroD/Reactor.cpp index c578f4ba4..1d2077c8a 100644 --- a/Cantera/src/zeroD/Reactor.cpp +++ b/Cantera/src/zeroD/Reactor.cpp @@ -26,7 +26,6 @@ namespace Cantera { Reactor::Reactor() : ReactorBase(), FuncEval(), m_kin(0), - m_integ(0), m_temp_atol(1.e-11), m_maxstep(0.0), m_vdot(0.0), @@ -35,13 +34,15 @@ namespace Cantera { m_chem(true), m_energy(true) { +#ifdef INCL_REACTOR_INTEG m_integ = new CVodeInt; // use backward differencing, with a full Jacobian computed // numerically, and use a Newton linear iterator m_integ->setMethod(BDF_Method); m_integ->setProblemType(DENSE + NOJAC); - m_integ->setIterator(Newton_Iter); + m_integ->setIterator(Newton_Iter); +#endif } @@ -79,7 +80,8 @@ namespace Cantera { for (int m = 0; m < m_nwalls; m++) { surf = m_wall[m]->surface(m_lr[m]); if (surf) { - surf->getCoverages(y+loc); + m_wall[m]->getCoverages(m_lr[m], y + loc); + //surf->getCoverages(y+loc); loc += surf->nSpecies(); } } @@ -96,12 +98,13 @@ namespace Cantera { for (int w = 0; w < m_nwalls; w++) if (m_wall[w]->surface(m_lr[w])) m_nv += m_wall[w]->surface(m_lr[w])->nSpecies(); +#ifdef INCL_REACTOR_INTEG m_atol.resize(neq()); fill(m_atol.begin(), m_atol.end(), 1.e-15); m_integ->setTolerances(m_rtol, neq(), m_atol.begin()); m_integ->setMaxStep(m_maxstep); m_integ->initialize(t0, *this); - +#endif m_enthalpy = m_thermo->enthalpy_mass(); m_pressure = m_thermo->pressure(); m_intEnergy = m_thermo->intEnergy_mass(); @@ -140,33 +143,26 @@ namespace Cantera { doublereal* mss = y + 2; doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0); m_mix->setMassFractions(mss); + m_mix->setDensity(mass/m_vol); doublereal temp = temperature(); mix.setTemperature(temp); if (m_energy) { - doublereal u_mass = u/mass; // specific int. energy - doublereal delta; - - do { - delta = -(m_thermo->intEnergy_mass() - - u_mass)/m_thermo->cv_mass(); - temp += delta; - mix.setTemperature(temp); - } - while (fabs(delta) > m_temp_atol); + m_thermo->setState_UV(u/mass,m_vol/mass); + temp = mix.temperature(); //mix.setTemperature(temp); } - mix.setTemperature(temp); - m_state[0] = temp; + //m_state[0] = temp; int loc = m_nsp + 2; SurfPhase* surf; for (int m = 0; m < m_nwalls; m++) { surf = m_wall[m]->surface(m_lr[m]); if (surf) { - surf->setTemperature(temp); - surf->setCoverages(y+loc); + // surf->setTemperature(temp); + //surf->setCoverages(y+loc); + m_wall[m]->setCoverages(m_lr[m], y+loc); loc += surf->nSpecies(); } } @@ -175,7 +171,7 @@ namespace Cantera { m_enthalpy = m_thermo->enthalpy_mass(); m_pressure = m_thermo->pressure(); m_intEnergy = m_thermo->intEnergy_mass(); - + //m_kappa = m_thermo->isothermalCompressibility(); m_mix->saveState(m_state); } @@ -217,6 +213,8 @@ namespace Cantera { rs0 = 1.0/surf->siteDensity(); nk = surf->nSpecies(); sum = 0.0; + surf->setTemperature(m_state[0]); + m_wall[i]->syncCoverages(m_lr[i]); kin->getNetProductionRates(m_work.begin()); ns = kin->surfacePhaseIndex(); surfloc = kin->kineticsSpeciesIndex(0,ns); @@ -281,12 +279,16 @@ namespace Cantera { int n; doublereal mdot_out; for (i = 0; i < m_nOutlets; i++) { - mdot_out = m_outlet[i]->massFlowRate(); + mdot_out = m_outlet[i]->massFlowRate(time); for (n = 0; n < m_nsp; n++) { ydot[2+n] -= mdot_out * mf[n]; } - if (m_energy) + if (m_energy) { + // cout << "before = " << ydot[0] << endl; ydot[0] -= mdot_out * enthalpy; + //cout << mdot_out << " " << enthalpy << endl; + //cout << "after = " << ydot[0] << endl; + } } @@ -294,12 +296,13 @@ namespace Cantera { doublereal mdot_in; for (i = 0; i < m_nInlets; i++) { - mdot_in = m_inlet[i]->massFlowRate(); + mdot_in = m_inlet[i]->massFlowRate(time); for (n = 0; n < m_nsp; n++) { - ydot[2+n] += m_inlet[i]->massFlowRate(n); + ydot[2+n] += m_inlet[i]->outletSpeciesMassFlowRate(n); } - if (m_energy) + if (m_energy) { ydot[0] += mdot_in * m_inlet[i]->enthalpy_mass(); + } } } } diff --git a/Cantera/src/zeroD/Reactor.h b/Cantera/src/zeroD/Reactor.h index 0a34078c8..5c9aba71e 100644 --- a/Cantera/src/zeroD/Reactor.h +++ b/Cantera/src/zeroD/Reactor.h @@ -69,7 +69,11 @@ namespace Cantera { /** * Destructor. Deletes the integrator. */ - virtual ~Reactor(){ delete m_integ; } + virtual ~Reactor(){ +#ifdef INCL_REACTOR_INTEG + delete m_integ; +#endif +} virtual int type() const { return ReactorType; } @@ -84,6 +88,7 @@ namespace Cantera { * @param time Final time (s). */ virtual void advance(doublereal time) { +#ifdef INCL_REACTOR_INTEG if (!m_init) { setMaxStep(time); initialize(); @@ -92,9 +97,14 @@ namespace Cantera { m_time = time; updateState(m_integ->solution()); m_mix->saveState(m_state); +#else + throw CanteraError("Reactor::advance", + "Reactor::advance is deprecated. Use ReactorNet::advance"); +#endif } virtual double step(doublereal time) { +#ifdef INCL_REACTOR_INTEG if (!m_init) { setMaxStep(time); initialize(); @@ -103,6 +113,10 @@ namespace Cantera { updateState(m_integ->solution()); m_mix->saveState(m_state); return m_time; +#else + throw CanteraError("Reactor::step", + "Reactor::step is deprecated. Use ReactorNet::step"); +#endif } /** @@ -128,53 +142,6 @@ namespace Cantera { m_maxstep = maxstep; } - // /** -// * Set the reactor surface area [m$^2$]. Can be changed at any time. -// */ -// void setArea(doublereal area) { -// m_area = area; -// } - -// /** -// * Set the external temperature \f$ T_0 \f$ -// * used for heat loss calculations. -// * The heat loss rate is calculated from -// * \f[ -// * \dot Q_{out} = h A (T - T_0) + \epsilon A (T^4 - T_{0,R}^4). -// * \f] -// * @see setArea, setEmissivity, setExtRadTemp -// */ -// void setExtTemp(doublereal ts) { -// m_ext_temp = ts; -// if (!m_trad_set) m_ext_temp4 = ts*ts*ts*ts; -// } - -// /** -// * Set the external temperature for radiation. By default, this -// * is the same as the temperature set by setExtTemp. But if -// * setExtRadTemp is called, then subsequent of calls to -// * setExtTemp do not modify the value set here. -// */ -// void setExtRadTemp(doublereal tr) { -// m_ext_temp4 = tr*tr*tr*tr; -// } - -// void setHeatTransferCoeff(doublereal h) { -// m_h = h; -// } - -// void setVDotCoeff(doublereal k) { -// m_kv = k; -// } - -// void setEmissivity(doublereal emis) { -// m_emis = emis; -// } - -// void setExtPressure(doublereal p0) { -// m_p0 = p0; -// } - void disableChemistry() { m_chem = false; } void enableChemistry() { m_chem = true; } @@ -208,61 +175,6 @@ namespace Cantera { virtual void initialize(doublereal t0 = 0.0); void evalEqs(doublereal t, doublereal* y, doublereal* ydot); - /** - * @name Methods to specify simulation options. - * These virtual methods may be overloaded in - * derived classes to implement models for heat gain/loss, - * surface chemistry, and compression/expansion. - */ - //@{ - - /** - * Initialize the boundary conditions, if necessary. This - * method does nothing, but may be overloaded in derived classes if - * initialization is needed. - */ - // virtual void initBC() {} - - /** - * Evaluate the reactor boundary conditions. This procedure is - * called during integration to evaluate the rate of volume - * change \f$ dV/dt \f$ [m^3/s], the heat loss rate [W], and - * the species production rates due to surface chemistry. - * - * It may be overloaded in derived classes to implement other - * boundary conditions. If not overloaded, this routine - * implements the following boundary conditions. - * - * The rate of volume change is - * \f[ - * dV/dt = K ( P - P_{ext}) - * \f] - * where K is set in procedure setVDotCoeff. - * - * - * The heat loss rate is calculated from - * \f[ - * \dot Q_{out} = h A (T - T_0) + \epsilon A (T^4 - T_{0,R}^4). - * \f] - * @see setArea, setEmissivity, setExtRadTemp - */ - -// virtual void evalBC(doublereal& vdot, -// doublereal& heatLossRate, doublereal* sdot) { -// doublereal t = m_mix->temperature(); - -// //m_p0 = m_env->pressure(); -// vdot = m_kv * (m_thermo->pressure()/m_p0 - 1.0)*m_vol0; -// heatLossRate = m_area * ( -// m_h * (t - m_ext_temp) -// + m_emis * StefanBoltz * (t*t*t*t - m_ext_temp4) -// ); -// } - - //@} - - //----------------------------------------------------- - /** * Set the mixture to a state consistent with solution * vector y. @@ -273,21 +185,15 @@ namespace Cantera { protected: Kinetics* m_kin; - // ReactorBase* m_env; - // Thermo* m_thermo; Integrator* m_integ; // pointer to integrator doublereal m_temp_atol; // tolerance on T doublereal m_maxstep; // max step size doublereal m_vdot, m_Q; - // doublereal m_emis, m_h, m_area; - //doublereal m_ext_temp, m_ext_temp4; - //doublereal m_kv, m_p0; vector_fp m_atol; doublereal m_rtol; vector_fp m_work; vector_fp m_sdot; // surface production rates - //bool m_trad_set; bool m_chem; bool m_energy; int m_nv; diff --git a/Cantera/src/zeroD/ReactorBase.cpp b/Cantera/src/zeroD/ReactorBase.cpp index c0591a783..774d9f939 100644 --- a/Cantera/src/zeroD/ReactorBase.cpp +++ b/Cantera/src/zeroD/ReactorBase.cpp @@ -32,6 +32,7 @@ namespace Cantera { m_enthalpy(0.0), m_intEnergy(0.0), m_pressure(0.0), + m_kappa(0.0), m_nwalls(0) { m_name = name; @@ -50,9 +51,11 @@ namespace Cantera { m_thermo = &thermo; m_nsp = m_mix->nSpecies(); m_mix->saveState(m_state); + m_rho0 = m_thermo->density(); m_enthalpy = m_thermo->enthalpy_mass(); m_intEnergy = m_thermo->intEnergy_mass(); m_pressure = m_thermo->pressure(); + m_kappa = m_thermo->isothermalCompressibility(); } void ReactorBase::addInlet(FlowDevice& inlet) { diff --git a/Cantera/src/zeroD/ReactorBase.h b/Cantera/src/zeroD/ReactorBase.h index d1afd0daf..66f5af50a 100644 --- a/Cantera/src/zeroD/ReactorBase.h +++ b/Cantera/src/zeroD/ReactorBase.h @@ -134,6 +134,8 @@ namespace Cantera { doublereal mass() const { return m_vol * density(); } const doublereal* massFractions() const { return m_state.begin() + 2; } doublereal massFraction(int k) const { return m_state[k+2]; } + doublereal compressibility() const { return m_kappa; } + //@} int error(string msg) const { @@ -142,7 +144,7 @@ namespace Cantera { } protected: - + int m_nsp; thermo_t* m_mix; thermo_t* m_thermo; @@ -160,7 +162,9 @@ namespace Cantera { vector_int m_lr; int m_nwalls; string m_name; - + double m_rho0; + double m_kappa; + private: void tilt(string method="") const { diff --git a/Cantera/src/zeroD/ReactorNet.cpp b/Cantera/src/zeroD/ReactorNet.cpp index 741c9390f..bc96aea55 100644 --- a/Cantera/src/zeroD/ReactorNet.cpp +++ b/Cantera/src/zeroD/ReactorNet.cpp @@ -4,7 +4,8 @@ namespace Cantera { ReactorNet::ReactorNet() : FuncEval(), m_nr(0), m_nreactors(0), - m_integ(0), m_init(false), m_nv(0), m_rtol(1.0e-6), + m_integ(0), m_init(false), + m_nv(0), m_rtol(1.0e-6), m_verbose(false) { m_integ = new CVodeInt; @@ -23,6 +24,9 @@ namespace Cantera { m_nv = 0; m_reactors.clear(); m_nreactors = 0; + if (m_verbose) { + writelog("Initializing reactor network.\n"); + } for (n = 0; n < m_nr; n++) { if (m_r[n]->type() == ReactorType) { m_r[n]->initialize(t0); @@ -54,7 +58,7 @@ namespace Cantera { void ReactorNet::advance(doublereal time) { if (!m_init) { - m_maxstep = time; + m_maxstep = time - m_time; initialize(); } m_integ->integrate(time); @@ -64,7 +68,7 @@ namespace Cantera { double ReactorNet::step(doublereal time) { if (!m_init) { - m_maxstep = time; + m_maxstep = time - m_time; initialize(); } m_time = m_integ->step(time); @@ -75,10 +79,15 @@ namespace Cantera { void ReactorNet::eval(doublereal t, doublereal* y, doublereal* ydot) { int n; int start = 0; - updateState(y); - for (n = 0; n < m_nreactors; n++) { - m_reactors[n]->evalEqs(t, y + start, ydot + start); - start += m_size[n]; + try { + updateState(y); + for (n = 0; n < m_nreactors; n++) { + m_reactors[n]->evalEqs(t, y + start, ydot + start); + start += m_size[n]; + } + } + catch (CanteraError) { + showErrors(cout); } } diff --git a/Cantera/src/zeroD/Wall.cpp b/Cantera/src/zeroD/Wall.cpp index 5b1c61fe6..b835fafd4 100644 --- a/Cantera/src/zeroD/Wall.cpp +++ b/Cantera/src/zeroD/Wall.cpp @@ -9,7 +9,7 @@ namespace Cantera { Wall::Wall() : m_left(0), m_right(0), - m_area(0.0), m_k(0.0), m_rrth(0.0), + m_area(0.0), m_k(0.0), m_rrth(0.0), m_emiss(0.0), m_vf(0), m_qf(0) { for (int n = 0; n < 2; n++) { m_chem[n] = 0; @@ -41,6 +41,8 @@ namespace Cantera { if (ileft >= 0) { m_surf[0] = (SurfPhase*)&left->thermo(ileft); m_nsp[0] = m_surf[0]->nSpecies(); + m_leftcov.resize(m_nsp[0]); + m_surf[0]->getCoverages(m_leftcov.begin()); } } if (right) { @@ -48,6 +50,8 @@ namespace Cantera { if (iright >= 0) { m_surf[1] = (SurfPhase*)&right->thermo(iright); m_nsp[1] = m_surf[1]->nSpecies(); + m_rightcov.resize(m_nsp[1]); + m_surf[1]->getCoverages(m_rightcov.begin()); } } if (ileft < 0 || iright < 0) { @@ -61,11 +65,14 @@ namespace Cantera { * The volume rate of change is given by * \f[ \dot V = K A (P_{left} - P_{right}) + F(t) \f] * where \f$ F(t) \f$ is a specified function of time. + * + * This method is used by class Reactor to compute the + * rate of volume change of the reactor. */ doublereal Wall::vdot(doublereal t) { double rate1 = m_k * m_area * (m_left->pressure() - m_right->pressure()); - if (m_vf) rate1 += m_vf->eval(t); + if (m_vf) rate1 += m_area * m_vf->eval(t); return rate1; } @@ -78,8 +85,33 @@ namespace Cantera { doublereal Wall::Q(doublereal t) { double q1 = (m_area * m_rrth) * (m_left->temperature() - m_right->temperature()); + if (m_emiss > 0.0) { + double tl = m_left->temperature(); + double tr = m_right->temperature(); + q1 += m_area * StefanBoltz * (tl*tl*tl*tl - tr*tr*tr*tr); + } if (m_qf) q1 += m_area * m_qf->eval(t); return q1; } + void Wall::setCoverages(int leftright, const doublereal* cov) { + if (leftright == 0) + copy(cov, cov + m_nsp[0], m_leftcov.begin()); + else + copy(cov, cov + m_nsp[1], m_rightcov.begin()); + } + + void Wall::getCoverages(int leftright, doublereal* cov) { + if (leftright == 0) + copy(m_leftcov.begin(), m_leftcov.end(), cov); + else + copy(m_rightcov.begin(), m_rightcov.end(), cov); + } + + void Wall::syncCoverages(int leftright) { + if (leftright == 0) + m_surf[0]->setCoverages(m_leftcov.begin()); + else + m_surf[1]->setCoverages(m_rightcov.begin()); + } } diff --git a/Cantera/src/zeroD/Wall.h b/Cantera/src/zeroD/Wall.h index 43d3e9eea..46fecfba9 100644 --- a/Cantera/src/zeroD/Wall.h +++ b/Cantera/src/zeroD/Wall.h @@ -57,8 +57,13 @@ namespace Cantera { /// Set the overall heat transfer coefficient [W/m^2/K]. void setHeatTransferCoeff(doublereal U) { m_rrth = U; } - /** Set the rate of volume change to a specified function.*/ - void setExpansionRate(Func1* f=0) {if (f) m_vf = f;} + void setEmissivity(doublereal epsilon) { m_emiss = epsilon; } + + // /** Set the rate of volume change to a specified function.*/ + // void setExpansionRate(Func1* f=0) {if (f) m_vf = f;} + + /** Set the piston velocity to a specified function. */ + void setVelocity(Func1* f=0) {if (f) m_vf = f;} /** * Set the expansion rate coefficient. @@ -98,6 +103,13 @@ namespace Cantera { return m_chem[leftright]; } + void setCoverages(int leftright, const doublereal* cov); + + void getCoverages(int leftright, doublereal* cov); + + void syncCoverages(int leftright); + + protected: vector_fp m_coeffs; @@ -108,8 +120,10 @@ namespace Cantera { SurfPhase* m_surf[2]; int m_nsp[2]; doublereal m_area, m_k, m_rrth; + doublereal m_emiss; Func1 *m_vf; Func1 *m_qf; + vector_fp m_leftcov, m_rightcov; private: diff --git a/Cantera/src/zeroD/flowControllers.h b/Cantera/src/zeroD/flowControllers.h index 1123ccadd..b9e3f9b96 100644 --- a/Cantera/src/zeroD/flowControllers.h +++ b/Cantera/src/zeroD/flowControllers.h @@ -21,141 +21,107 @@ #include "FlowDevice.h" #include "ReactorBase.h" -#include "PID_Controller.h" +//#include "PID_Controller.h" #include "../Func1.h" namespace Cantera { - /////////////////////////////////////////////////////////////// - - - /** - * A base class for devices that do not use closed-loop control. - * This is defined only for convenience, in order to overload - * virtual methods of FlowDevice that print warnings with ones - * that do nothing. - */ - class NoController : public FlowDevice { - public: - - NoController() {} - virtual ~NoController() {} - - NoController(const NoController& a) - : FlowDevice(a) {} - - NoController& operator=(const NoController& a) { - return *this; - } - - // unneeded methods - virtual void update() {} - virtual void reset() {} - virtual bool setGains(int n, const doublereal* gains) {return true;} - virtual bool getGains(int n, doublereal* gains) {return true;} - virtual doublereal maxError() { return 0.0; } - virtual doublereal setpoint() { return 0.0; } - virtual void setSetpoint(doublereal mdot) { } - virtual bool ready() { - return FlowDevice::ready(); - } - - protected: - private: - }; - - - ////////////////////////////////////////////////////////// - - /** * A class for mass flow controllers. The mass flow rate is constant, * independent of any other parameters. */ - class MassFlowController : public NoController { + class MassFlowController : public FlowDevice { public: - MassFlowController() {} - virtual ~MassFlowController() {} - - MassFlowController(const MassFlowController& a) - : NoController(a) {} - - MassFlowController& operator=(const MassFlowController& a) { - if (this == &a) return *this; - m_mdot = a.m_mdot; - return *this; + MassFlowController() : FlowDevice() { + m_type = MFC_Type; } - virtual doublereal setpoint() { return m_mdot; } - virtual void setSetpoint(doublereal mdot) { m_mdot = mdot; } + virtual ~MassFlowController() {} + virtual bool ready() { return FlowDevice::ready() && m_mdot >= 0.0; } - protected: - private: - }; - - - class UserValve : public NoController { - public: - - UserValve() : m_func(0) {} - virtual ~UserValve() {} - - UserValve(const UserValve& a) : NoController(a) {} - - UserValve& operator=(const UserValve& a) { - if (this == &a) return *this; - m_func = a.m_func; - return *this; - } - - virtual bool ready() { - return FlowDevice::ready() && m_func != 0; - } - - virtual void setFunction(Func1* f) { m_func = f; } - - virtual doublereal massFlowRate() { - return m_func->eval(in().pressure() - out().pressure()); + /// If a function of time has been specified for + /// mdot, then update the stored mass flow rate. + /// Otherwise, mdot is a constant, and does not need + /// updating. + virtual void updateMassFlowRate(doublereal time) { + if (m_func) m_mdot = m_func->eval(time); + if (m_mdot < 0.0) m_mdot = 0.0; } protected: - Func1* m_func; private: }; - /** * A class for mass flow controllers. The mass flow rate is constant, * independent of any other parameters. */ - class Valve : public NoController { + class PressureController : public FlowDevice { public: - Valve() {} - virtual ~Valve() {} - - Valve(const Valve& a) : NoController(a) {} - - Valve& operator=(const Valve& a) { - if (this == &a) return *this; - m_mdot = a.m_mdot; - return *this; + PressureController() : FlowDevice(), m_master(0) { + m_type = PressureController_Type; } + virtual ~PressureController() {} + + virtual bool ready() { + return FlowDevice::ready() && m_master != 0; + } + + void setMaster(FlowDevice* master) { + m_master = master; + } + + virtual void updateMassFlowRate(doublereal time) { + doublereal master_mdot = m_master->massFlowRate(time); + m_mdot = master_mdot + m_coeffs[0]*(in().pressure() - + out().pressure()); + if (m_mdot < 0.0) m_mdot = 0.0; + } + + protected: + FlowDevice* m_master; + + private: + }; + + + /// Valve objects supply a mass flow rate that is a function of the + /// pressure drop across the valve. The default behavior is a linearly + /// proportional to the pressure difference. Note that + /// real valves do not have this behavior, so this class + /// does not model real, physical valves. + class Valve : public FlowDevice { + public: + + Valve() : FlowDevice() { + m_type = Valve_Type; + } + + virtual ~Valve() {} + virtual bool ready() { return FlowDevice::ready() && m_coeffs.size() >= 1; } - virtual doublereal massFlowRate() { - m_mdot = m_coeffs[0]* (in().pressure() - out().pressure()); + /// Compute the currrent mass flow rate, based on + /// the pressure difference. + virtual void updateMassFlowRate(doublereal time) { + double delta_P = in().pressure() - out().pressure(); + if (m_func) { + m_mdot = m_func->eval(delta_P); + } + else { + m_mdot = m_coeffs[0]*delta_P; + } if (m_mdot < 0.0) m_mdot = 0.0; - return m_mdot; } protected: @@ -163,57 +129,6 @@ namespace Cantera { private: }; - - - /** - * A PressureRegulator is a device that controls the pressure - * of the upstream reactor by regulating the mass flow rate. - */ - class PressureRegulator : public FlowDevice { - - public: - - PressureRegulator() {} - virtual ~PressureRegulator() {} - PressureRegulator(const PressureRegulator& p) : m_pid(p.m_pid) {} - PressureRegulator& operator=(const PressureRegulator& p) { - if (this == &p) return *this; - m_pid = p.m_pid; - return *this; - } - - // overloaded virtual methods - virtual void setSetpoint(doublereal p0) { m_pid.setpoint(-p0); } - virtual doublereal setpoint() { return -m_pid.setpoint(); } - virtual bool ready() { - return FlowDevice::ready() && m_pid.ready(); } - virtual void reset() { - m_pid.reset(in().time()-1.e-12, -in().pressure()); - } - virtual void update() { - m_pid.update(in().time(), -in().pressure()); - } - - virtual bool setGains(int n, const doublereal* gains) { - return m_pid.setGains(n, gains); - } - - virtual bool getGains(int n, doublereal* gains) { - return m_pid.getGains(n, gains); - } - - virtual doublereal maxError() { return m_pid.maxError(); } - - virtual doublereal massFlowRate() { - m_mdot = m_pid.output(-in().pressure()); - return m_mdot; - } - - protected: - - private: - PID_Controller m_pid; - }; } #endif diff --git a/configure b/configure index 66a44aca6..72790f46e 100755 --- a/configure +++ b/configure @@ -183,7 +183,7 @@ LAPACK_FTN_STRING_LEN_AT_END='y' CXX=${CXX:=g++} # C++ compiler flags -CXXFLAGS=${CXXFLAGS:="-O0 -g -Wall"} +CXXFLAGS=${CXXFLAGS:="-O0 -Wall"} # the C++ flags required for linking #LCXX_FLAGS= diff --git a/data/inputs/ptcombust.cti b/data/inputs/ptcombust.cti index 4cfd1a643..43154883c 100644 --- a/data/inputs/ptcombust.cti +++ b/data/inputs/ptcombust.cti @@ -1,5 +1,6 @@ # -# see http://reaflow.iwr.uni-heidelberg.de/~Olaf.Deutschmann/ for more about this mechanism +# see http://reaflow.iwr.uni-heidelberg.de/~Olaf.Deutschmann/ for +# more about this mechanism # #---------------------------------------------------------------------! #*********************************************************************** @@ -28,9 +29,17 @@ units(length = "cm", time = "s", quantity = "mol", act_energy = "J/mol") + +# +# Define a gas mixture with species imported from GRI-Mech. +# Reactions will be imported from GRI-Mech 3.0, as long as they +# don't involve species not declared here. Transport properties +# will be computed using a mixture-averaged model. +# ideal_gas(name = "gas", elements = "O H C N Ar", - species = """gri30: H2 H O O2 OH H2O HO2 H2O2 + species = """gri30: H2 H O O2 OH + H2O HO2 H2O2 C CH CH2 CH2(S) CH3 CH4 CO CO2 HCO CH2O CH2OH CH3O CH3OH C2H C2H2 C2H3 C2H4 C2H5 C2H6 HCCO CH2CO HCCOH AR N2""", @@ -55,12 +64,12 @@ ideal_interface(name = "Pt_surf", coverages = 'O(S):0.0, PT(S):0.5, H(S):0.5') ) -#------------------------------------------------------------------------------- +#----------------------------------------------------------------------------- # Species data # # Note that reactions 12-14 are reversible, and therefore require thermo # data -#------------------------------------------------------------------------------- +#----------------------------------------------------------------------------- species(name = "PT(S)", atoms = " Pt:1 ", @@ -202,22 +211,23 @@ species(name = "O(S)", # Reaction 1 surface_reaction("H2 + 2 PT(S) => 2 H(S)", [4.45790E+10, 0.5, 0], - order = "PT(S):1") + order = "PT(S):1") # Reaction 2 surface_reaction( "2 H(S) => H2 + 2 PT(S)", - Arrhenius(3.70000E+21, 0, 67400, coverage = ['H(S)', 0.0, 0.0, -6000.0])) + Arrhenius(3.70000E+21, 0, 67400, + coverage = ['H(S)', 0.0, 0.0, -6000.0])) # Reaction 3 surface_reaction( "H + PT(S) => H(S)", stick(1.00000E+00, 0, 0)) # Reaction 4 surface_reaction( "O2 + 2 PT(S) => 2 O(S)", Arrhenius(1.80000E+21, -0.5, 0), - options = 'duplicate') + options = 'duplicate') # Reaction 5 surface_reaction( "O2 + 2 PT(S) => 2 O(S)", stick(2.30000E-02, 0, 0), - options = 'duplicate') + options = 'duplicate') # Reaction 6 surface_reaction( "2 O(S) => O2 + 2 PT(S)", @@ -266,7 +276,8 @@ surface_reaction( "CH4 + 2 PT(S) => CH3(S) + H(S)", [4.63340E+20, 0.5, 0], order = "PT(S):2.3") # Reaction 20 -surface_reaction( "CH3(S) + PT(S) => CH2(S)s + H(S)", [3.70000E+21, 0, 20000]) +surface_reaction( "CH3(S) + PT(S) => CH2(S)s + H(S)", + [3.70000E+21, 0, 20000]) # Reaction 21 surface_reaction( "CH2(S)s + PT(S) => CH(S) + H(S)", [3.70000E+21, 0, 20000]) diff --git a/ext/tpx/Sub.h b/ext/tpx/Sub.h index c0c0d78d7..cf20e223b 100755 --- a/ext/tpx/Sub.h +++ b/ext/tpx/Sub.h @@ -133,7 +133,7 @@ namespace tpx { return T*(s2 - s1)/(2.0*dt); } - virtual double thermExpCoeff() { + virtual double thermalExpansionCoeff() { double Tsave = T, dt = 1.e-4*T; double p0 = P(); Set(TP, Tsave - dt, p0); @@ -144,6 +144,16 @@ namespace tpx { return (v2 - v1)/((v2 + v1)*dt); } + virtual double isothermalCompressibility() { + double Psave = P(), dp = 1.e-4*Psave; + Set(TP, T, Psave - dp); + double v1 = v(); + Set(TP, T, Psave + dp); + double v2 = v(); + Set(TP, T, Psave); + return -(v2 - v1)/((v2 + v1)*dp); + } + // saturation properties