The curse of modeling detailed chemistry

In this article I will tell the story behind how I managed to go from no knowledge about detailed chemical kinetics to sort of mastering it in a few months. By the end of the article I will illustrate a practical problem solved with a Python package I wrote after recycling some code from my PhD thesis.

Example of acetylene pyrolysis simulation results in a PFR.

If you are studying Chemical Engineering or a related field you have probably heard about ideal reactors and how they can be used for solving real world problems. The mathematics behind ideal reactors is quite interesting because several cases allow for analytical solutions of the differential equations behind the scenes et voilà, all your life’s problems are solved! Although this scenario is as ideal as the reactors are during an exam, practical cases go far beyond the analytical solutions.

Before proceeding with the actual content, let me introduce you to my background because maybe that will make my approach to reactor engineering clearer. First of all, I am not a Chemical Engineer and I hadn’t studied reactor models since my high school years. One morning I receive a confirmation that I had been accepted to a PhD in Materials Science at Université de Lorraine, France. Great! I had applied to that position knowingly that about 50–60% of the proposed research content matched my background only — I am a Materials Engineer. Now I had a lot of new stuff to learn regarding the about 40–50% of Chemical Engineering involved in my thesis proposal description! It was 2013, I quit my job as a Product Engineer at a major Oil & Gas company and move to France.

First steps

The project I got involved required me to study half the time about metallurgy, the other half being dedicated to the gas phase processes linked to the treatment of the alloys I studied. Today I would not recommend anyone to accept a PhD position as generalist like that unless your target is to work in classical industry later on. If you target is really innovative fields and startups, do one thing or the other, but not a generalist PhD! Well, so I finish that with no regrets, it was too late to say I would not pursue it according to my morals. My first task was to investigate the pyrolysis of acetylene under conditions relevant to the process of steel carburizing, i.e. pressures below 5000 Pa and temperatures in the order of 1200 K (I will only use SI units for physical formality), in a tubular reactor. A typical task for someone in the field of pyrolysis or combustion, not to me — at least at that time. I needed at least some reading about reactor engineering.

After some online research I found the great textbook by H. Scott Fogler and immediately bought it. It was some sort of discovery of a whole new field to me. In a few weeks I finished reading most chapters and working out examples. Given the flow controllers I had a the lab it was possible to compute the range of gas speeds I should expect inside the reactor, and thus use my recently acquired knowledge about Péclet number to determine that a Plug Flow Reactor moel would represent quite well my experimental system given the diffusivity of the main species from the problem I found on the literature. I felt ready to proceed — and I was wrong!

If you have ever got in touch with combustion or pyrolysis of gas phase products, probably you have already understood my mistake. Such problems involve rather complex kinetic mechanisms, with dozens to hundreds of chemical species and maybe thousands of reactions. Everything I had learned with the book I had just finished represented only the extreme surface of a much deeper ocean! For modeling detailed chemical kinetics systems much more is required. That moment I talked to our IT guys in the lab and asked what tools for this sort of problem were available — and we didn’t have a Chemkin lisence, the long term default tool for this sort of problems. OK, I sort of like doing difficult stuff and told myself to start developing a similar tool. By that time I had no knowledge of the huge scientific open source community and my programming skills were limited to some Matlab.

I learned C++ and after a few days of coding I managed to interpret the 240-species mechanism text file describing acetylene pyrolysis that was provided by one of my references and computing reaction rates that matched reported values. That seemed fairly well, but not well enough. Once you have forward reaction rates, there is still the need to interpret standardized polynomials (generally the so-called NASA polynomials) to compute backward processes from equilibrium hypothesis, and perform energy balance. Task accomplished! It was time to integrate the system! But wait, my background in numerical Mathematics was very limited — and I didn’t know that. Without any further reading, I started coding a simple semi-implicit Euler’s method solver, and obviously failed! What you have to know before integrating large chemical kinetics systems is that they are generally stiff. Stiff ordinary differential equations are difficult to integrate for reasons you can check in the link. I had no time to develop and validate a stiff equation solver, time to move on.

Who would save me?

I started a deep online research of tools that were capable of treating my problem. For some weird reason, the first open source tool I found was the old — and sort of abandoned — Sandia Laboratory’s TChem. The library works quite well and you can learn a lot with its documentation — really recommend the reading. The main issue with using this package is that you will probably be one of the few users around the globe, and generally you will have to modify the source code in case of any bugs or if you need extra functionalities because there is no one you can talk to. Furthermore, you need to have good knowledge of library compilation in order to use the tool. By the other side, there is a thing TChem put me in touch with that pays for all the pain I had first time using the library: that day I got to know Sundials. Sundials is possibly the most useful library I have ever used. This suite of solvers can be used in any field and has great performance. With Sundials I finally used my home-made code and verified I had the same results as TChem, but my code still had too many bugs to pursue the calculations with it in a safe manner.

That moment I knew I was close to getting ready to attack my real world problem. The initial hurry to integrate the equations was gone. The path reading the book, coding a lot, and using third-party libraries thought me how hard it was to have this sort of problem solved — at least the first time you do it by yourself. It is curious the fact that that same day, while searching for other stuff online, I discovered about Cantera. That moment I almost cried — seriously. Cantera was everything I needed and have a great community and support behind it. And Cantera also uses Sundials for integrating the equations as a default solver!

Life with Cantera

A whole new era during my PhD studies started with Cantera. I no longer needed to develop everything by myself, so I could spend more time simulating scenarios that could be interesting for planning my experimental tests. Furthermore, Cantera has an excellent Python interface, so all the work of coding, compiling and executing was replaced by scripting and executing. What a time saver! It has been about 4.5 years that I am using this amazing tool, and a lot has improved during this time. Since I am no longer in academic research, sadly it is difficult for me finding some time to contribute with official Cantera releases. Nonetheless, I am often promoting the community and this week I finally managed to organize some code used during my thesis to provide some plug flow reactor (PFR) models with help of Cantera — currently Cantera does not provide any PFR model as a built-in feature, but it is quite easy to build your own if you have some programming background.

The package CanteraPFR is available on GitHub and you can simulate any gas phase process with one of its PFR models in just a few lines of code. Some examples are available here. Just to illustrate the simplicity of using this package to integrate a PFR, the following minimum script presents the integration of mechanism.cti (an hypothetical gas phase mechanism) in an isothermal scenario. Inlet composition is provided according to Cantera’s format string and the only non-SI unit is the flow rate. This value must be provided in sccm because this is the most common laboratory scale flow rate unit (so that most users do not need to apply temperature corrections). Solution of a simple problem can be done just be creating a PyPFR object and then calling its method manage_solution, which will provide you a CSV file containing the atmosphere composition every dx meters as well as the velocity, density and pressure along the reactor axis.

# -*- coding: utf-8 -*-
from CanteraPFR import PyPFR
mechf = 'mechanism.cti' # Mechanism file
phase = 'gas' # Name of phase in mechanism
# Composition string (as accepted by Cantera)
X0 = 'N2:0.75, C2H2:0.25'
T0 = 1200.0 # Inlet temperature [K]
P0 = 5000.0 # Inlet pressure [Pa]
Q0 = 500.0 # Flow rate [sccm]
Di = 0.025 # Reactor diameter [m]
Lr = 1.00 # Reactor length [m]
dx = 0.01 # Steps to write results [m]
prob = PyPFR('isothermal', mechf, phase, Di, T0, p0, X0, Q0) prob.manage_solution(Lr, dx, 'solution.csv')


Modeling of detailed chemical kinetics may take huge efforts of beginners if they are not supplied with the correct tools. Current state-of-art open source libraries, specially Cantera can provide a much steeper learning curve and reliable results. While you learn Cantera you have the bonus of improving your programing skills, what is quite interesting. Hope you liked the post and if you find CanteraPFR implementation useful, please let me know — or if you wish to request any modifications and new features. Have a great learning!

They told me I was a Materials Scientist. Well, nowadays a feel like this is too narrow. I am a Scientist, that’s all. Let’s attack the numerical problems!

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store