Introduction

In this project the combustion of a hydrogen-oxygen mixture is studied. If a high precision in the time resolution of this reaction is wanted, regarding all intermediate products, it is not sufficient just to regard the simplified process 2H2 + O2 -> 2H2O

Instead a lot of intermediate states like H, O, OH, HO2, H2O2 and their equilibrium reactions need to be considered.

From this a set of first order differential equations can be derived. These can be integrated using suitable algorithms.

This work provides a tool set that is capable of reading in and interpreting the equilibrium reaction coefficients in a somewhat human-readable format (see FORMAT TODO), to derive the rate equations, to solve the differential equation in an adjusted way and to write the time evolution of the involved initial, intermediate and final substances.

I also used this project as an opportunity to update my compiler and to dive in to the world of C++11 [7]. Therefore I am using some features that are only implemented by newer C++ compilers [12]. The make-file requires the GCC compiler in version 4.7 or higher. (Tested with version 4.8.)

This text was written using the facilities of the documentation systen Doxygen. Thus it consists of two parts:

Chapter 1 contains the lab report for this project
Chapter 2 & 3 contain the documentation of the written code as produced by Doxygen.

Chemical formalism

In this program the modified Arrhenius equation is used for calculation of the reaction constants:

k = A T^n exp -E/RT.

This differs from the original Arrhenius equation [4],

k = A exp -E/RT,

in the additional temperature difference that can occur in some processes.

The reaction rates are then calculated in the usual manner from the reaction constants:

r = k [R_1]^{n_1} [R_2]^{n_2} …,

where [R_i] denotes the concentration of the reagent i and n_i its multiplicity in the reaction formula.

The data for some of the reactions were designed to use in the Troe formalism [11] [6]. However, for reasons of simplicity the shorter formula

k = k_0 k_∞ / (k_0 + k_∞)

was used in this project. The classes Reaction and ReactionParser have included nonetheless the functionality to handle the full Troe formalism. Only the maintaining class CombustionProcess needs to be modified to activate the usage of this formalism.

Data basis

As a data basis the data for the involved reactions should be used.

The list of reactions stated in the task sheet for this project turned out to be incomplete in terms of giving all required coefficients. Thus it was looked for the original source of the data. The trace was:

The picture in the task sheet was copied from a picture in a talk [1].
This talk took the picture from a table in the dissertation [3].
[3] simply copied the table from [2].
[2] took the data from [8], but omitted some important data.

In most of the papers it was not clearly displayed where the data was taken from. This is a story of copying data without copying all the data or explicitely stating, where the data comes from.

So finally the data from table I in [8] (which is fortunately also mentioned in the task sheet) was used.

The change of (Gibbs-)energy was not denoted in [8] for the reactions (though I later found out that it could be calculated from the data in table II.) Thus the data from the older paper [9] was used, where the energy change is state explicitely per-reaction. The energy change is denoted as Delta H in [9], but since the values were mostly the same as in [8] there were taken to be Delta G .

Results

See PDF document.

The results were obtained using the following calls of the programm:

./h2combustion -T 1900 -d 1e-14 -p 1e-2 -t 1e-7 -o T1900ap2t7 -a
./h2combustion -T 2900 -d 1e-14 -p 1e-2 -t 1e-9 -o T2900ap2t9 -a
./h2combustion -T 900 -d 1e-14 -p 1e-2 -t 1e-7 -o T900ap2t7 -a

./h2combustion -T 2900 -d 1e-14 -p 1e-2 -t 1e-9 -o T2900d14t9
./h2combustion -T 900 -d 1e-14 -p 1e-2 -t 1e-7 -o T900d14t7

./h2combustion -T 1900 -d 1e-14 -p 1e-2 -t 7.65e-10 -o T1900adiabatic -a --adiabatic

The values 10^-14 and 10^-2 for the timestep and the accuracy respectively have been chosen to be the highest values that leave the integration to be stable for the chosen temperatures.

Usage of the program

Compile

A make-file is available. Enter make to compile the project. Enter make clean to remove all generated binary and temporary files. The make-file assumes G++ version 4.7 or higher [5].

Run

The application, once successfully compiled, can be called by simply typing ./h2combustion The program is the run with some sane parameters, reads its input from the files reactions.dat and initialconditions.dat and writes its output to data.dat.

Commandline options

-r <filename>: Specify the file of reactions.
-i <filename>: Specify the file of initial conditions.
-o <filename>: Specify the output file.
-o stdout: Print output on screen (Unix only.)
-T <value>: Specify the initial temperature.
-t <value>: Specify the integration time.
-d <value>: Specify the time step of integration.
-a: Select adaptive time step integration.
-p <value>: Specify the accuracy in adaptive time step integration.
--adiabatic: Allowing the temperature to change in the integration (non-isothermic).

Remarks about the design

The design, as can be seen, is clearly influenced by the style of object-oriented-programming (OOP). It makes heavy use of techniques like encapsulation, inheritance and polymorphism.

However, I have read the Stroustroup [10] meanwhile and thought about some things I would change if I were to design this project from the beginning again.

Obviously the function diffeq() in the class Integrable and derived classes is a frequently called function. It is declared virtual, so the integrator does not need to now what special kind of system it is dealing with. However, virtual function calls can be hardly be optimised by an compiler and always have some overhead due to the v-table mechanism. It would have been better (in the sense of more efficient) to implement the integrator as a generic class with a template parameter of the system, which could be represented by a functor for example (and thus allowing for function pointers as well, for the differential equation.) This would be much better for a compiler to optimise, and would probably also lead to slightly more readable code.