From 8a76c29d2813ba4819cb0cfb30103b574484959f Mon Sep 17 00:00:00 2001 From: Walter Dal'Maz Silva Date: Fri, 4 Jan 2019 09:05:26 +0100 Subject: [PATCH] First release --- CanteraPFR/CanteraPFR.pyx | 217 ++++++++++++++++-- CanteraPFR/doc/source/development.rst | 33 +++ CanteraPFR/doc/source/index.rst | 8 +- CanteraPFR/doc/source/interfaces.rst | 9 + CanteraPFR/doc/source/models.rst | 165 +++++++++++++ CanteraPFR/doc/source/models/adiabatic.rst | 34 --- CanteraPFR/doc/source/models/base.rst | 11 - CanteraPFR/doc/source/models/isothermal.rst | 136 ----------- .../include/CanteraPFR/AdiabaticPFR.hpp | 2 - CanteraPFR/include/CanteraPFR/HeatWallPFR.hpp | 2 - .../include/CanteraPFR/IsothermalPFR.hpp | 2 - CanteraPFR/include/CanteraPFR/SolvePFR.hpp | 5 +- CanteraPFR/src/HeatWallPFR.cpp | 2 +- CanteraPFR/test/test.py | 16 +- README.md | 18 +- 15 files changed, 437 insertions(+), 223 deletions(-) create mode 100644 CanteraPFR/doc/source/development.rst create mode 100644 CanteraPFR/doc/source/interfaces.rst create mode 100644 CanteraPFR/doc/source/models.rst delete mode 100644 CanteraPFR/doc/source/models/adiabatic.rst delete mode 100644 CanteraPFR/doc/source/models/base.rst delete mode 100644 CanteraPFR/doc/source/models/isothermal.rst diff --git a/CanteraPFR/CanteraPFR.pyx b/CanteraPFR/CanteraPFR.pyx index 10248f6..68e9666 100644 --- a/CanteraPFR/CanteraPFR.pyx +++ b/CanteraPFR/CanteraPFR.pyx @@ -1,30 +1,45 @@ # -*- coding: utf-8 -*- + from libcpp.string cimport string from CanteraPFR cimport AdiabaticPFR from CanteraPFR cimport HeatWallPFR from CanteraPFR cimport IsothermalPFR from CanteraPFR cimport SolvePFR +import os import time import ctypes from numpy import arange from numpy import array from pandas import DataFrame -# See https://stackoverflow.com/questions/51044122 -# https://github.com/JohannesBuchner/PyMultiNest/issues/5 +# cdef class Closure: + """ Just-in-time compiler for Python functions. + + Allows a Python function to be called from C/C++ library. + Based on https://stackoverflow.com/questions/51044122. + + TODO + ---- + This class is incoherent, once the interface of `_inner_fun` is fixed + and does not follow what `ftype` requires. This must be wrapped. + + Parameters + ---------- + python_fun : function + Python function to call from C/C++. + ftype : ctypes.CFUNCTYPE + Interface to provide the function. + """ + cdef object python_fun cdef object jit_wrap - def __cinit__(self, python_fun): + def __cinit__(self, python_fun, ftype): self.python_fun = python_fun - ftype = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double) - self.jit_wrap = ftype(self.inner_fun) - - def inner_fun(self, double arg): - return self.python_fun(arg) + self.jit_wrap = ftype(lambda *args: self.python_fun(*args)) cdef func_t get_fun_ptr(self): return (ctypes.addressof(self.jit_wrap))[0] @@ -34,44 +49,81 @@ cdef class PyPFR: """ Python wrapper to plug-flow reactor models. This class is built around several plug-flow reactor models. The selection - of these is made according to the provided model name. + of these is made according to the provided model name: + + - `isothermal`: isothermal PFR (temperature held constant). + - `adiabatic`: adiabatic PFR (only viscous losses). + - `heatwall`: PFR with imposed wall temperature profile. TODO ---- - Provide docstrings to all methods. - Make standard graphical ouput. - Include theory in rst documentation. + 1. Deal with Python strings. + 2. Make standard graphical ouput. + 3. Add interface to user viscosity function. + 4. Make htc a function of gas speed. + 5. Allow vector of absolute tolerances. Parameters ---------- + rtype : str + Key-word for model selection. + mech : str + Path to CTI/XML mechanism file. + phase : str + Phase name in mechanism file. + Di : float + Reactor diameter in meters. + T0 : float + Inlet temperature in kelvin. + p0 : float + Inlet pressure in pascal. + X0 : str + Inlet composition compatible with Cantera in string format. + Q0 : float + Volume flow rate at reference state (assumed to be 273.15 K and + 101325 Pa) given in cubic centimeters per minute. This unit was chosen + because it is the most common used with laboratory scale mass flow + controllers. + htc : float, optional + Wall convective heat transfer coefficient in watts per kelvin per square + meter. If its value is zero, this is equivalent to an adiabatic reactor, + regardless of wall temperature. This parameter is intended to be used + with Tw. Must be supplied if `rtype='heatwall'`. Default is None. + Tw : float or function `float(float)`, optional + Wall temperature profile in kelvin. This can be supplied as a constant + value or as a function of position along reactor axis (start coordinate + at zero). Must be supplied if `rtype='heatwall'`. Default is None. + + Raises + ------ + SystemExit + If reactor type name is not recognized or if provided wall temperature + function is not suitable for use with reactor model `heatwall`. """ + cdef CanteraPFR* obj cdef SolvePFR* sol cdef Closure cl def __cinit__(self, rtype, mech, phase, Di, T0, p0, X0, Q0, htc=None, Tw=None): - # TODO deal with Python strings! # cdef bytearray cpp_mech = bytearray(mech, 'utf8') # cdef string cpp_mech = mech.encode('utf-8') print('*' * 79) if rtype.lower() == 'isothermal': self.obj = new IsothermalPFR(mech, phase, Di, T0, p0, X0, Q0) + elif rtype.lower() == 'adiabatic': self.obj = new AdiabaticPFR(mech, phase, Di, T0, p0, X0, Q0) + elif rtype.lower() == 'heatwall': assert htc is not None, 'Missing heat transfer coefficient' assert Tw is not None, 'Missing wall temperature' - - if isinstance(Tw, float): - self.cl = Closure(lambda x: Tw) - else: - # Some test here! - self.cl = Closure(Tw) - + self.cl = self._get_closure(Tw) self.obj = new HeatWallPFR(mech, phase, Di, T0, p0, X0, Q0, htc, self.cl.get_fun_ptr()) + else: raise SystemExit(f'Unknown reactor type : {rtype}') @@ -83,25 +135,140 @@ cdef class PyPFR: del self.obj del self.sol + def _get_closure(self, Tw): + ftype = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double) + + if isinstance(Tw, float): + return Closure(lambda x: Tw, ftype) + + try: + Tw(0.0) + except Exception as err: + raise SystemExit(f'Provide function is not suitable : {err}') + + return Closure(Tw, ftype) + def set_tolerances(self, rtol, atol): + """ Set solution tolerances. + + Parameters + ---------- + rtol : float + Relative tolerance applied to all variables. + atol : float + Absolute tolerance applied to all variables. + """ self.sol.setTolerances(rtol, atol) def set_max_num_steps(self, num): + """ Set maximum number of integration steps. + + This represents the maximum number of internal steps taken by IDA while + integrating the problem, not the number of outputs made by the user. + + Parameters + ---------- + num : int + Maximum number of integration steps. + """ self.sol.setMaxNumSteps(num) def set_initial_step_size(self, h0): + """ Set length of initial integration step. + + Notice that a value too small here may lead to failure if the maximum + number of steps as provided by :meth:`set_max_num_steps` is not + consistent. Do not use this method if you are not facing failures at + integration start or precision issues. + + Parameters + ---------- + h0 : float + Length of initial integration step. + """ self.sol.setInitialStepSize(h0) - def set_stop_position(self, tstop): - self.sol.setStopPosition(tstop) + def set_stop_position(self, xstop): + """ Set last coordinate to stop integration. + + Eventually IDA solver may trespass the last point to integrate. + If using :meth:`manage_solution`, this is automatically done for you. + + Parameters + ---------- + xstop : float + Coordinate of last point in meters. + """ + self.sol.setStopPosition(xstop) + + def solve(self, xout): + """ Integrate for current time to the provided time. - def solve(self, tout): - return self.sol.solve(tout) + Notice that the solver will fail if `xout` is less than current step + and no test is performed by this class. + + Parameters + ---------- + xout : float + Coordinate of point to return solution. + + Returns + ------- + int + State of solution as defined by IDA constants. + """ + return self.sol.solve(xout) def solution(self, num): + """ Returns solution component by its index. + + Solution array depends on the used model, so no general description can + be provided here. For advanced usage only. + + Parameters + ---------- + num : int + Index of solution component to return. + + Returns + ------- + float + Value of solution component. + """ return self.sol.solution(num) def manage_solution(self, Lr, dx, saveas, outfreq=100): + """ Manage integration of reactorself. + + Parameters + ---------- + Lr : float + Total length of reactor in meters. + dx : float + Step to save results. + saveas : str + Path to results file. + outfreq : int, optional + Count of steps to provide screen output. Default is `100`. + + Raises + ------ + AssertionError + If parameters are unphysical (negative lengths), if step is below + total length, if output file exists, or negative output frequency. + + Returns + ------- + pandas.DataFrame + Table of results for further manipulation. + """ + + assert Lr > 0, 'Reactor length must be positive' + assert dx > 0, 'Saving step must be positive' + assert Lr > dx, 'Reactor length must be above saving step' + assert not os.path.exists(saveas), 'Output file already exists' + assert outfreq > 0, 'Output frequency must be positive' + columns = self.sol.variablesNames() columns = ['x'] + [c.decode('utf-8') for c in columns] sol0 = array([0] + list(self.sol.solutionVector())) @@ -126,3 +293,5 @@ cdef class PyPFR: print(f'Integration took {time.time()-t0:.5f} s') data.to_csv(saveas, index=False) + + return data diff --git a/CanteraPFR/doc/source/development.rst b/CanteraPFR/doc/source/development.rst new file mode 100644 index 0000000..46f38e4 --- /dev/null +++ b/CanteraPFR/doc/source/development.rst @@ -0,0 +1,33 @@ +Development +=========== + +General Practices +----------------- + +Since this project depends on external libraries, it is important to keep a +standard for updates. First, the target programming language for the package is +Python. That said, applications are expected to be conceived in that language. +Nonetheless, the main interface of the library used to develop these models is +provided in C++. This implies that development workflow should follow the path: + +1. Conceive new models in C++: + a. Header files shall be placed under *include/CanteraPFR*. + b. Source files shall be placed under *src*. + c. Test interface for C++ is provided in *test/main*. + d. Edit *Makefile* to include test binary. +2. Provide Cython interface: + a. Add C++ interfaces to *CanteraPFR.pdx*. + b. Implement Python callable interface to new models in *CanteraPFR.pyx*. + c. Add test to new model to *test/test.py*. +3. Build and test project with new model: + a. Run `make` from main directory. + b. Execute *test.py* from its directory. + c. Run `make` from *doc* directory for documentation. + +Functionalities to implement +---------------------------- + +1. Models are not taking cross-section variations into account. +2. Generalize model by incorporating source (surface) terms. +3. Implement test for turbulent limit in momentum loss equation. +4. Provide analytical Jacobian for all models. diff --git a/CanteraPFR/doc/source/index.rst b/CanteraPFR/doc/source/index.rst index cbf9552..1dfe627 100644 --- a/CanteraPFR/doc/source/index.rst +++ b/CanteraPFR/doc/source/index.rst @@ -6,11 +6,11 @@ Welcome to CanteraPFR's documentation! .. toctree:: :maxdepth: 2 - :caption: Contents: + :caption: Package Contents: - models/base - models/adiabatic - models/isothermal + models + interfaces + development references Indices and tables diff --git a/CanteraPFR/doc/source/interfaces.rst b/CanteraPFR/doc/source/interfaces.rst new file mode 100644 index 0000000..73fb639 --- /dev/null +++ b/CanteraPFR/doc/source/interfaces.rst @@ -0,0 +1,9 @@ +Reference API +============= + +.. module:: CanteraPFR + +.. autoclass:: PyPFR + :members: + :undoc-members: + :show-inheritance: diff --git a/CanteraPFR/doc/source/models.rst b/CanteraPFR/doc/source/models.rst new file mode 100644 index 0000000..4f4045e --- /dev/null +++ b/CanteraPFR/doc/source/models.rst @@ -0,0 +1,165 @@ +Plug-Flow Reactor Model +======================= + +A plug-flow reactor (PFR) model represents the limit of a tubular reactor for +which Péclet number in axial direction is too high so that mixing can be +neglected over flow direction (so that convective effects completely shadow +any diffusion of species). That implies segregated flow in which each volume +of fluid behaves as a closed system and travels in a single direction without +any mixing with its neighboring elements of fluid. The composition at each +coordinate of the reactor is homogeneous (a PFR can be seen under some cases +as a series of perfect-stirred reactors) so that any radial effects are +neglected. It is essentially a one-dimensional model and for most purposes is +treated on stationary state, since the non-mixing criterion makes the transient +approach meaningless. This package deals with PFR's that are both constant and +smoothly varying (**to be implemented**) in cross-sectional area (assumed +circular). Each of these cases can be treated according to slightly different +formulations, what includes: + +1. isothermal PFR: is characterized by the absence of energy equation solution. +2. adiabatic PFR: heat released by reactions is not exchanged with environment. +3. imposed wall temperature PFR: simplified convection heat exchanged with environment. + +**Currently no surface exchanges are implemented, that must be a next step. +This will require a way to generalize the interface to place a catalyst or +reacting interface inside the system, so that it can be usable in many cases.** + +Numerical integration +--------------------- + +Solution of a PFR requires at least the coupled integration of mass, momentum, +and species conservation equations. For those cases for which temperature is +solved, energy conservation must be added to this set. These equations are then +coupled by means of ideal gas law. Following, equations solved by this package +are presented. For more details, please check reference [Ref.1]_. + +Throughout the discussion the following symbols convention is adopted: + +- :math:`\rho`: density of gas. +- :math:`\mu`: viscosity of gas. +- :math:`T`: temperature of gas. +- :math:`u`: velocity in axial direction. +- :math:`R`: radius of reactor. +- :math:`A`: cross-section of chamber. +- :math:`V`: volume of chamber. +- :math:`\mathcal{P}_{a}`: perimeter over cross-section :math:`A` of tube. +- :math:`\hat{h}`: global convective heat transfer coefficient. +- :math:`T_{w}`: local wall temperature. +- :math:`C_{p}`: mixture constant pressure heat capacity. +- :math:`h_{k}`: species :math:`k` enthalpy per unit mass. +- :math:`W`: molecular weight, with subscript to represent a species. +- :math:`\bar{W}`: mean molecular weight. +- :math:`Y`: mass fraction, with subscript to represent a species. +- :math:`\omega`: molar rate of production, with subscript to represent a species. +- :math:`CV`: control volume, used as integration limit. +- :math:`CS`: control surface, used as integration limit. + +Mass continuity +--------------- + +Mass continuity equation in the absence of wall (source) terms can be derived +by expanding steady state general compressible continuity equation +:math:`\nabla\cdot{}(\rho{}u)`. The full derivation is provided in [Ref.1]_ and +makes use of Reynolds transport theorem. The expansion of this expression over +:math:`z` axis gives the expression implemented in this package: + +.. math:: + + \rho \frac{\partial u}{\partial z} + + u \frac{\partial \rho}{\partial z} = 0 + +Species conservation +-------------------- + +Species conservation starts by applying Reynolds transport theorem on the overall +system species mass conservation. By using divergence theorem (here specialized +in one dimension) this can be written as: + +.. math:: + + \int_{CV}\frac{\partial \rho{}Y_{k}}{\partial t}dV + + \int_{CS}\rho{}Y_{k}u{}dA = + \int_{CV}\omega_{k}W_{k}dV + +Assuming steady state and no variations in channel: + +.. math:: + + \frac{\partial \rho{}uY_{k}}{\partial z} - \omega_{k}W_{k} = 0 + +The axial derivative can be simplified using the overall mass balance leading to +the expression that is implemented: + +.. math:: + + \rho{}u\frac{\partial Y_{k}}{\partial z} - \omega_{k}W_{k} = 0 + +Momentum conservation +--------------------- + +Momentum conservation is derived for a Newtonian fluid by taking into account +viscous effects (drag over channel walls) and pressure acting at inlet. Viscous +effects are introduced through a model equation as discussed in [Ref.1]_. The +equation is deriving by assuming a smooth pressure gradient and a friction factor +applicable for the laminar limit. Hydraulic diameter proposed in section 16.2 of +[Ref.1]_ is used to simplify the equation to its final form as: + +.. math:: + + \rho{}u\frac{\partial u}{\partial z} + + \frac{\partial p}{\partial z} + + 8\frac{\mu u}{R^2} = 0 + +State equation +-------------- + +Coupling of state variables is made by ideal gas law, here implemented as: + +.. math:: + + p\bar{W} = \rho{}RT + +Energy conservation +------------------- + +Except for the isothermal case, energy equation is included in a general form as +follows. Notice that for adiabatic reactors :math:`\hat{h}=0`, what makes the +convective term disappear from the equation. + +.. math:: + + \rho{}u{}C_{p}\frac{\partial T}{\partial z} + + \sum_{k}h_{k}\omega_{k}W_{k} - + \mathcal{P}_{a}\hat{h}(T_{w}-T) = 0 + +Problem initialization +---------------------- + +Since we intend to solve the problem as a differential algebraic equations +system, the constraining equations for calculation of initial derivatives +is required. This is done by replacing the symbols in previous equations +by their respective initial states. In order to get the same number of +unknowns as equations, implicit total derivative of state equation is added to +the system. These equations form a set of linear equations that can be solved +by standard linear algebra methods. + +.. math:: + + \begin{align*} + \omega_{k,0} W_{k} + & = \rho_{0} u_{0} Y^{\prime}_{k,0} \\ + 0 & = \rho_{0} u^{\prime}_{0} + u_{0} \rho^{\prime}_{0} \\ + -8\frac{\mu_{0} u_{0}}{R^2} + & = \rho_{0} u_{0} u^{\prime}_{0} + p^{\prime}_{0} \\ + 0 & = p_{0}\frac{\bar{W}^{2}_{0}}{W_{k,0}}Y^{\prime}_{k,0} + + RT\rho^{\prime}_{0} + - \bar{W}_{0}p^{\prime}_{0} + + \rho_{0}RT^{\prime}_{0} \\ + -\sum_{k}h_{k,0}\omega_{k,0} W_{k} + +\mathcal{P}_{a}\hat{h}(T_{w,0}-T_{0}) + & = \rho_{0} u_{0} C_{p,0} T^{\prime}_{0} + \end{align*} + +Notice that here equations were expressed in most general form considered in this +package. Temperature derivatives have to be eliminated for isothermal case and +convective heat transfer coefficient set to zero in adiabatic limit. diff --git a/CanteraPFR/doc/source/models/adiabatic.rst b/CanteraPFR/doc/source/models/adiabatic.rst deleted file mode 100644 index b5fdf30..0000000 --- a/CanteraPFR/doc/source/models/adiabatic.rst +++ /dev/null @@ -1,34 +0,0 @@ -AdiabaticPFR Plug-Flow Reactor Model -==================================== - -Model limitations/hypothesis: - - -Functionalities to implement: - - -For more details, please check reference [Ref.1]_. - -Model equations ---------------- - -Throughout the discussion the following symbols convention is adopted: - -- :math:`\rho`: density of gas. -- :math:`\mu`: viscosity of gas. -- :math:`u`: velocity in axial direction. -- :math:`R`: radius of reactor. -- :math:`V`: volume of chamber. -- :math:`A`: cross-section of chamber. -- :math:`W`: molecular weight, subscripted to represent a species. -- :math:`\bar{W}`: mean molecular weight. -- :math:`Y`: mass fraction, subscripted to represent a species. -- :math:`\omega`: molar rate of production, subscripted to represent a species. -- :math:`CV`: control volume, used as integration limit. -- :math:`CS`: control surface, used as integration limit. - -Module documentation --------------------- - -.. autoclass:: CanteraPFR.models._hom_adiabatic.AdiabaticPFR - :members: diff --git a/CanteraPFR/doc/source/models/base.rst b/CanteraPFR/doc/source/models/base.rst deleted file mode 100644 index 79d0ca2..0000000 --- a/CanteraPFR/doc/source/models/base.rst +++ /dev/null @@ -1,11 +0,0 @@ -BasePFR Plug-Flow Reactor Model -=============================== - -Module documentation --------------------- - -.. autoclass:: CanteraPFR.models._base_pfr.BaseSolution - :members: - -.. autoclass:: CanteraPFR.models._base_pfr.BasePFR - :members: diff --git a/CanteraPFR/doc/source/models/isothermal.rst b/CanteraPFR/doc/source/models/isothermal.rst deleted file mode 100644 index e72bd32..0000000 --- a/CanteraPFR/doc/source/models/isothermal.rst +++ /dev/null @@ -1,136 +0,0 @@ -Isothermal Plug-Flow Reactor Model -================================== - -An isothermal plug-flow reactor (IPFR) is characterized by the absence of energy -equation solution in space and a segregated flow in which each volume of fluid -behaves as a closed system and travels in a single direction without any mixing -with its neighboring elements of fluid. In this sense, it is a system presenting -a Péclet number that is too high, so that convective effects completely shadow -any diffusion of species. The model hereafter described assumes, thus, an -one-dimensional flow in a channel of constant or slowly varying section. The -composition at each coordinate of the reactor is homogeneous (a PFR can be seen -under some circumstances as a series of perfect-stirred reactors) so that any -radial effects are neglected. Finally, a plug-flow model is essentially a steady -state model. - -Solution of an IPFR require the coupled integration of mass, momentum, and species -conservation equations. This model implements only gas phase phenomena in such a -way that mass conservation is linked to density and velocity changes only. Species -mole/mass fractions can vary only due to homogeneous processes, which can lead to -changes in average molecular weight and impact directly the momentum and mass -equations. Momentum equation may take viscous effects into account and then be -able to properly represent the pressure loss along the reactor axis. - -Model limitations/hypothesis: - -1. Isothermal system (not suitable for flames). -2. Laminar flow in circular tubular reactor. -3. Constant reactor cross-section. -4. Homogeneous reactions only. -5. Ideal gas equation of state. - -Functionalities to implement: - -1. Model is not taking cross-section variations into account. -2. Generalize model by incorporating source (surface) terms. - -For more details, please check reference [Ref.1]_. - -Model equations ---------------- - -Throughout the discussion the following symbols convention is adopted: - -- :math:`\rho`: density of gas. -- :math:`\mu`: viscosity of gas. -- :math:`u`: velocity in axial direction. -- :math:`R`: radius of reactor. -- :math:`V`: volume of chamber. -- :math:`A`: cross-section of chamber. -- :math:`W`: molecular weight, subscripted to represent a species. -- :math:`\bar{W}`: mean molecular weight. -- :math:`Y`: mass fraction, subscripted to represent a species. -- :math:`\omega`: molar rate of production, subscripted to represent a species. -- :math:`CV`: control volume, used as integration limit. -- :math:`CS`: control surface, used as integration limit. - -Mass continuity equation in the absence of wall (source) terms can be derived -by expanding steady state general compressible continuity equation -:math:`\nabla\cdot{}(\rho{}u)`. The full derivation is provided in [Ref.1]_ and -makes use of Reynolds transport theorem. The expansion of this expression over -:math:`z` axis gives the expression implemented in this model: - -.. math:: - - \rho \frac{\partial u}{\partial z} + - u \frac{\partial \rho}{\partial z} = 0 - -Species conservation starts by applying Reynolds transport theorem on the overall -system species mass conservation. By using divergence theorem (here specialized -in one dimension) this can be written as: - -.. math:: - - \int_{CV}\frac{\partial \rho{}Y_{k}}{\partial t}dV + - \int_{CS}\rho{}Y_{k}u{}dA = - \int_{CV}\omega_{k}W_{k}dV - -Assuming steady state and no variations in channel: - -.. math:: - - \frac{\partial \rho{}uY_{k}}{\partial z} - \omega_{k}W_{k} = 0 - -The axial derivative can be simplified using the overall mass balance leading to -the expression that is implemented: - -.. math:: - - \rho{}u\frac{\partial Y_{k}}{\partial z} - \omega_{k}W_{k} = 0 - -Momentum conservation is derived for a Newtonian fluid by taking into account -viscous effects (drag over channel walls) and pressure acting at inlet. Viscous -effects are introduced through a model equation as discussed in [Ref.1]_. The -equation is deriving by assuming a smooth pressure gradient and a friction factor -applicable for the laminar limit. Hydraulic diameter proposed in section 16.2 of -[Ref.1]_ is used to simplify the equation to its final form as: - -.. math:: - - \rho{}u\frac{\partial u}{\partial z} + - \frac{\partial p}{\partial z} + - 8\frac{\mu u}{R^2} = 0 - -Coupling of state variables is made by ideal gas law, here implemented as: - -.. math:: - - p\bar{W} = \rho{}RT - -Since we intend to solve the problem as a differential algebraic equations -system, the constraining equations for calculation of initial derivatives -need to be provided. This is done by replacing the symbols in previous equations -by their respective initial states. In order to get the same number of -unknowns as equations, implicit total derivative of state equation is added to -the system. These equations form a set of linear equations that can be solved -by standard linear algebra methods. - -.. math:: - - \begin{align*} - 0 & = \rho_{0} u^{\prime}_{0} + u_{0} \rho^{\prime}_{0} \\ - 0 & = \rho_{0} u_{0} Y^{\prime}_{k,0} - \omega_{k,0} W_{k} \\ - 0 & = \rho_{0} u_{0} u^{\prime}_{0} + p^{\prime}_{0} - + 8\frac{\mu_{0} u_{0}}{R^2} \\ - 0 & = RT\rho^{\prime}_{0} + p_{0}\frac{\bar{W}^{2}_{0}}{W_{k,0}} - - \bar{W}_{0}p^{\prime}_{0} - \end{align*} - -Module documentation --------------------- - -.. autoclass:: CanteraPFR.models._hom_isothermal.IsothermalPFRSolution - :members: - -.. autoclass:: CanteraPFR.models._hom_isothermal.IsothermalPFR - :members: diff --git a/CanteraPFR/include/CanteraPFR/AdiabaticPFR.hpp b/CanteraPFR/include/CanteraPFR/AdiabaticPFR.hpp index 54eb49e..4311d04 100644 --- a/CanteraPFR/include/CanteraPFR/AdiabaticPFR.hpp +++ b/CanteraPFR/include/CanteraPFR/AdiabaticPFR.hpp @@ -3,8 +3,6 @@ // // Author : Walter Dal'Maz Silva // Date : December 21st 2018 -// -// TODO provide analytical Jacobian. // *************************************************************************** #ifndef __ADIABATICPFR_HPP__ diff --git a/CanteraPFR/include/CanteraPFR/HeatWallPFR.hpp b/CanteraPFR/include/CanteraPFR/HeatWallPFR.hpp index 01eb576..e25e23c 100644 --- a/CanteraPFR/include/CanteraPFR/HeatWallPFR.hpp +++ b/CanteraPFR/include/CanteraPFR/HeatWallPFR.hpp @@ -3,8 +3,6 @@ // // Author : Walter Dal'Maz Silva // Date : December 31st 2018 -// -// TODO provide analytical Jacobian. // *************************************************************************** #ifndef __HEATWALLPFR_HPP__ diff --git a/CanteraPFR/include/CanteraPFR/IsothermalPFR.hpp b/CanteraPFR/include/CanteraPFR/IsothermalPFR.hpp index b0aa3b6..8b24c24 100644 --- a/CanteraPFR/include/CanteraPFR/IsothermalPFR.hpp +++ b/CanteraPFR/include/CanteraPFR/IsothermalPFR.hpp @@ -3,8 +3,6 @@ // // Author : Walter Dal'Maz Silva // Date : December 21st 2018 -// -// TODO provide analytical Jacobian. // *************************************************************************** #ifndef __ISOTHERMALPFR_HPP__ diff --git a/CanteraPFR/include/CanteraPFR/SolvePFR.hpp b/CanteraPFR/include/CanteraPFR/SolvePFR.hpp index 70ab137..302b650 100644 --- a/CanteraPFR/include/CanteraPFR/SolvePFR.hpp +++ b/CanteraPFR/include/CanteraPFR/SolvePFR.hpp @@ -3,10 +3,6 @@ // // Author : Walter Dal'Maz Silva // Date : December 30th 2018 -// -// TODO -// ==== -// Manage return codes from IDA_Solver.solve in wrapper function. // *************************************************************************** #ifndef __SOLVEPFR_HPP__ @@ -68,6 +64,7 @@ class SolvePFR int solve(doublereal tout) { + // TODO Manage return codes from IDA_Solver.solve. try { return m_solver->solve(tout); diff --git a/CanteraPFR/src/HeatWallPFR.cpp b/CanteraPFR/src/HeatWallPFR.cpp index 38dc4c4..efb10e9 100644 --- a/CanteraPFR/src/HeatWallPFR.cpp +++ b/CanteraPFR/src/HeatWallPFR.cpp @@ -18,7 +18,7 @@ int Cantera::HeatWallPFR::getInitialConditions(const doublereal t0, const doublereal rho0R = rho0 * Cantera::GasConstant; const doublereal rhoUCp = rho0 * m_u0 * m_gas->cp_mass(); doublereal hdot = 0; - doublereal wall = 0 ; //wallHeatExchange(t0, T0); + doublereal wall = wallHeatExchange(t0, T0); m_gas->getMassFractions(y); m_gas->getNetProductionRates(&m_wdot[0]); diff --git a/CanteraPFR/test/test.py b/CanteraPFR/test/test.py index ed3991a..18993be 100644 --- a/CanteraPFR/test/test.py +++ b/CanteraPFR/test/test.py @@ -57,6 +57,19 @@ def case02(): def case03(): + """ Usage of `heatwall` as adiabatic PFR model. """ + + htc = 0 + + rtype = 'heatwall' + prob = PyPFR(rtype, mech, phase, Di, T0, p0, X0, Q0, htc=htc, Tw=T0) + prob.set_tolerances(rtol, atol) + prob.set_max_num_steps(maxsteps) + prob.set_initial_step_size(dx0) + prob.manage_solution(Lr, dx, f'solution_{rtype}_0.csv', outfreq=10) + + +def case04(): """ Usage of PFR model with functional wall temperature. """ htc = 10 @@ -70,7 +83,7 @@ def case03(): prob.manage_solution(Lr, dx, f'solution_{rtype}_2.csv', outfreq=10) -def case04(): +def case05(): """ Usage of PFR model with constant wall temperature. """ htc = 10 @@ -88,3 +101,4 @@ def case04(): case02() case03() case04() +case05() diff --git a/README.md b/README.md index 5c580c4..22ea1bf 100644 --- a/README.md +++ b/README.md @@ -6,6 +6,21 @@ predefined wall profiles variations of PFR. # Installation +The following steps assume you do not dispose of a Cantera (at least 2.4.0 install in +your machine. **Notice that the following will only work if you compiled Cantera +locally with Sundials 2.7.0**. If you already have Cantera, just create a file named +`Makefile.in` in `CanteraPFR` directory (same level as `Makefile`) and write the +following contents to this file: + +``` +EXTERNAL := full path to your Cantera and Sundials-2.7.0 roots (must be the same). +OPTIONS := -g -O3 -std=c++11 +``` + +If Sundials is elsewhere, you will need to edit the `Makefile` itself, and if you +intend to do this I assume you know how to do that. If compiling from Cygwin, also +add `-U__STRICT_ANSI__` to `OPTIONS` line. + ## Linux Install `libopenblas-dev`, `libboost-dev` assuming you are using Ubuntu (or their @@ -52,8 +67,7 @@ line 20 of `get_cantera.sh`. # Documentation -**Outdated/deprecated:** the documentation will continue to be as follows once -Cython modules are written/updated. +**Only after building the project, otherwise Sphinx will not be able to import!** In order to generate the documentation, go to directory `doc/` and run `make`. This will provide a list of types you can compile the documentation. The default