Using chemkin files for *.cti creation

[Matthis Huhn]

Dear community,

I found this amazing page: https://www.detchem.com/mechanisms

There are a lot of diferent mechanism that can be used in cantera. I think even the base "CH4-O2 SURFACE MECHANISM ON PT"*.cti file is an import from this database.

I really want to make a new methanation mechanism work in cantera and I downloaded the chemkin and termo files from there. But I cannot solve the importing part.

I used the https://cantera.org/tutorials/ck2cti-tutorial.html tutorial and tried some stuff but ended up googeling and posting questions on researchgate. Maybe somebody is out here andis able to help me with it.

Thanks a bunch and take care.

Matthis

[Bryan Weber]

Hi,

Can you please be more specific about the error messages you're getting by copy-pasting them into a message?

Best,

Bryan

[Matthis Huhn]

Hi,

certainly:

I use the following command:

ck2cti --thermo=thermdata.dat --input=Methanation_Ni_CHEMKIN.inp --permissive

and I get the following output:

cantera.ck2cti.InputParseError: Unexpected token "H2+Ni+Ni" in reaction expression "H2+Ni+Ni=>H-Ni+H-Ni

The reaction expression is like the first expression in the file, so something is really messed up here.

In the past I used the exports from: https://cearun.grc.nasa.gov/ThermoBuild/ for other calculations and it worked there, so I think I confident, that everything 'under the hood' is set in the right order.

Thanks a lot and take care

matthis

[Matthis Huhn]

Here's the full output:

INFO:root:Error reading reaction entry starting on line 31:

WARNING:root:
ERROR: Unable to parse 'Methanation_Ni_CHEMKIN.inp' near line 31:

Traceback (most recent call last):
  File "C:\Users\ **** \Anaconda3\envs\Modell\lib\site-packages\cantera\ck2cti.py", line 1217, in readKineticsEntry
    locs[j] = int(token), 'coeff'
ValueError: invalid literal for int() with base 10: 'H2+Ni+Ni'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:\Users\ **** \Anaconda3\envs\Modell\lib\site-packages\cantera\ck2cti.py", line 1220, in readKineticsEntry
    locs[j] = float(token), 'coeff'
ValueError: could not convert string to float: 'H2+Ni+Ni'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:\Users\****\Anaconda3\envs\Modell\Scripts\ck2cti-script.py", line 10, in <module>
    sys.exit(script_entry_point())
  File "C:\Users\****\Anaconda3\envs\Modell\lib\site-packages\cantera\ck2cti.py", line 2334, in script_entry_point
    main(sys.argv[1:])
  File "C:\Users\ **** \Anaconda3\envs\Modell\lib\site-packages\cantera\ck2cti.py", line 2307, in main
    permissive=permissive)
  File "C:\Users\ **** \Anaconda3\envs\Modell\lib\site-packages\cantera\ck2cti.py", line 2200, in convertMech
    parser.loadChemkinFile(inputFile)
  File "C:\Users\ **** \Anaconda3\envs\Modell\lib\site-packages\cantera\ck2cti.py", line 1869, in loadChemkinFile
    reaction, revReaction = self.readKineticsEntry(kinetics, surface)
  File "C:\Users\ **** \Anaconda3\envs\Modell\lib\site-packages\cantera\ck2cti.py", line 1222, in readKineticsEntry
    raise InputParseError('Unexpected token "{0}" in reaction expression "{1}".'.format(token, original_reaction))

cantera.ck2cti.InputParseError: Unexpected token "H2+Ni+Ni" in reaction expression "H2+Ni+Ni=>H-Ni+H-Ni

".

Thanks a lot again

[Ray Speth]

Hi Matthis,

That input file doesn’t seem to be complete - it’s just a REACTIONS section, but for a surface reaction mechanism, you also need something to define the species in a SITE section, usually in the same input file, and something to define the gas phase species, usually in a separate file. If you have input files that include these sections, then the usual way to convert this with ck2cti (or better, ck2yaml) would be

ck2yaml --thermo=thermdata.dat --input=gas-definition.inp --surface=surface-definition.inp

Regards,
Ray

​

[kurthi DBFZ]

Hey Ray,

I am having the same problem. Could you provide an example of the used syntax in those files? I can get the 'thermo' file and the 'species' from detchem.com but the surface reaction file with the 'SITE' section.

Thank you in advance and take care

[Ray Speth]

Hi,

The example files we use for testing purposes are surface1.inp (which contains SITE, THERMO ALL, and REACTIONS sections) and surface1-gas.inp (which contains ELEMENTS, SPECIES, THERMO, and REACTIONS sections). This can be converted to the Cantera format by running:

ck2yaml --input=surfac1-gas.inp --surface=surface1.inp --output=surface1.yaml

If you have a separate file containing the thermo data, that can be specified with the --thermo option.

I’m not familiar with how these inputs are currently handled by Chemkin. Does it allow the inputs to be spread across even more files?

Regards,
Ray

​

[kurthi DBFZ]

Hi Ray,

Amazing!  Thank you very much. Your example file works like a charm and I can start from there!

I'll post the solution here, when I'm done!

Cheers and take care and please keep up the great work! You're making a difference!

[kurthi DBFZ]

yep,

I'm done and it was quite a walk because I didn't listen to the Tutorial (https://cantera.org/tutorials/ck2cti-tutorial.html).

• Some of the column specific positions are not important but the last integer must be on column 80
• The phase of the species must be on column  45
• Because I dealt with surface reactions it was important that I had two different input files and the catalytic surface had to begin with "_"... I don't know why but I duplicated it from Ray's files

My way was as follows:

1. create two files: One with the surface reaction, where all the species and elements are declared on top and one with the gas phase with the exact same syntac but of course without the associated species
   
2. Watch out for the columns in the files and keep close to the tutorial
3. When it comes to reaction mechanism avoid an "&" on the end of the line
4. Convert the files to *.cti or *.yaml with:
   

python -m cantera.ck2yaml --input=surfaceNi-gas.inp --surface=surfaceNi.inp --output=surfaceNi.yaml --permissive

Thank you very much it really helped me a lot, Ray!

Kurth

[kurthi DBFZ]

Hey guys,

I was succcesfull with the conversion of the file but I reran some tests and implemented the *.yaml and *.cti file in the pfr.py example file from cantera and now I'm confident that the files are somehow incorrect.

The reaction rates from the file are at all states '0'. I checked all the reaction mechanisms in the files and the dimension, values and number of reaction are right.

Do you guys have any other idea how my file can be wrong?

I'm in contact with the authors of the source but they couldn't help me with this one.

Any help is absolutely welcome

Thank you very much

kurth

[kurthi DBFZ]

Hey,

I think I can specify my problem a bit further: 

How can I use surface reactions in the custom description of the system, where I don't define a cantera reactor and thus don't declare the reacting surface?

My system is a 1-D PFR reactor solved with the solvers from the custom.py. The system is the CO2 methanation on Ni-catalyst. I implemented the yaml file as described above and tested it with the "methane_pox_on_pt" reactor and got reasonable results. But I don't see how I can implement the whole system of gas phase, surface phase and reacting interface within the custom implementation.

When I use:

`self.surf.get_net_production_rates(self.gas)`

I get only the desciption of the adorption of hydrogen and no reaction occurs. I need to use the custom ODE solvers from scipy because of the complexity of my system.

Do you guys have any better solution or any idea how I can advance my model?

have a great day and thanks in advance

kurth

[Ray Speth]

Hi Kurth,

Can you share the YAML file you generated and any of the code you’ve tried writing? I think some of the details there would be helpful in understanding what’s going on.

The function you mentioned, surf.get_net_production_rates, is specifically meant for calculating just the production rates in the specified phase (as written, self.gas), based on reactions occurring on the surface. So that should only include absorption and desorption reactions. The production rates of gas phase species due to gas phase reactions would be returned by self.gas.net_production_rates, and the production rates for surface species due to surface reactions would be returned by self.surf.get_net_production_rates(self.surf).

Regards,
Ray

​

[kurthi DBFZ]

Hi Ray,

 

thanks for your reply and your kind help!

 

attached you'll find the .yaml file I created. I moified the coverage of the catalyst and the coverage but everything else

When I run   ` self.gas.net_productio_rates ` I get the error:

```

***********************************************************************

NotImplementedError thrown by Kinetics::updateROP:

Not implemented.

***********************************************************************

```

The function call is the same as in the example:

 

 

```

class ReactorOde(object):

def __init__(self, gas, surf, u_0, dx, d_react):

 

self.gas = gas

self.surf = surf

 

self.mem_porosity = 0.38

self.x_out = dict( zip( self.gas.species_names, np.zeros(len(self.gas.species_names)) ) ) # Konzentration der Komponenten ausserhalb der Membran

self.d_0 = 2e-9 # pore size of membrane in m

self.L_mem = 35e-6 # Thickness of membrane in m

self.p_out = 1 # Druck ausserhalb der Membran in bar

 

self.epsilon_bed = 0.39 # Porosität des Bettes des Reaktors/void fraction

 

# Parameters of the ODE system and auxiliary data are stored in the

# ReactorOde object.

 

self.P = gas.P

self.T = gas.T

self.u_0 = u_0/self.epsilon_bed # von leerrohrgeschwindigkeit auf effektive Geschwindigkeit

self.dx = dx

self.d_react = d_react

 

# Ein dict Objekt erzeugen mit den Namen der Komponenten und dann als zweiten Teil, das Objekt der properties Klasse als reine Komponente. Dies basiert auf dem

# Thermp package von pypi

# self.properties = dict(zip( self.gas.species_names ,[( Chemical(i) ) for i in self.gas.species_names ] ) )

 

 

def __call__(self, t, y):

"""the ODE function, y' = f(t,y) """

# Kühlung

T_w = 350 + 273.15 # Wandtemperatur in K

U = 2000 # W/(m² * K)

 

# State vector is [T, Y_1, Y_2, ... Y_K]

self.gas.set_unnormalized_mole_fractions(y[1:]/self.gas.P)

 

self.gas.TP = y[0], self.gas.P

rho = self.gas.density_mole # kmol/m³

# wdot = self.kinetik # wdot ist in kmol/(m³ * s) Kinetik ist in mol/(s* kg_cat), deswegen wird der Wirkungsgrad und Dichte des Kats verwendet

 

Vo_Ar_fact = 4/self.d_react # Faktor, um die Oberfläche auf das Volumen zu beziehen in m²/m³

alpha_mem = 0.2 # Theoretischer Enhancementfaktor

alpha_heat = 1 # Idee einer besseren Wärmeübertragung durch die Feststoffschüttung

 

wdot = self.gas.net_production_rates

 

m_mem = 0

# m_mem = self.dgm_cantera * Vo_Ar_fact * alpha_mem # membran permeabilität auch in kmol/(m²*s), Durch Membranfaktor umformen auf kmol/m³-s

## Energiebilanz

 

dTdx = ( -(np.dot(self.gas.partial_molar_enthalpies, (wdot)) + (U * Vo_Ar_fact * ( gas.T - T_w ) * alpha_heat ) )/ # J/(s*m³)

(rho * gas.cp_mole) * self.u_0) # J/(m³*K) ---> # Energiebilanz in K/m

### Massenbilanz

dpdx = ((wdot-m_mem) * const.R * gas.T * 1000)/ (self.u_0) - \

(y[1:]/ gas.T * dTdx ) # 1/m # Es muss an dieser Stelle in 1/m sein, weil in diesem Fall der Ort diskretisiert wird

return np.hstack((dTdx, dpdx))

```

 

I can clearly see that something in the .yaml file is not set up correctly, but just following the tutorials I can't get my head around it. 

Thanks in advance and have a great day!

 

Kurth

[kurthi DBFZ]

Maybe in addition to that: When I use the same input data and take the .cti file, I don't get that error, I just end up with a reaction rate of zero. But when I use the same yaml or cti file in the surface reaction example provided by cantera, I get decent results. So I'm a bit puzzled if the error is in the yaml (or cti) file or in my method to call the reaction rates...

Thank you guys

Kurth 

[Ray Speth]

Hi Kurth,

Your YAML file doesn’t specify any gas phase reactions. Is this correct? In that case, the gas object won’t have kinetics model attached to it, so all of the kinetics-related properties will throw the NotImplementedError exception that you see.

The plug flow reactor with surface chemistry is actually a bit of a thorny problem, and so far we haven’t been able to get a reliable solver built into Cantera. One of the complexities is that it’s no longer just an ODE system, but a system of differential-algebraic equations, since the equations for the surface species don’t contain spatial derivatives of the coverages, and these equations seem to be very prone to problems with stiffness.

You might be interested in looking at the 1D PFR with surface chemistry Jupyter Notebook, which implements the system in Python. The one complexity I know of in this case is that the scikits.odes.dae solver that it uses can be a little difficult to install, depending on your environment. There was some work towards rewriting the existing FlowReactor model to allow surface chemistry, using the Sundials IDA solver. That can be seen in this pull request, but that work has been stalled.

Regards,
Ray

​

[kurthi DBFZ]

Hey Ray,

I think I have to stop working on the idea of implementing my problem in the standard cantera models, since the system is clearly a DAE system right now and as you told me that can't be solved right now using cantera.

and yep all the reaction take place on the surface and no reaction is defined in the gas phase. It's the same in the methane_pox_on_pt.yaml file and I get decent results in the surf_pfr.py file with my yaml file and with the methane_pox_on_pt.yaml file. So I thought this couldn't be the problem, because both system deal with "just surface chemistry".

This is why I was surprised by the "not implemented" error in my custom.py file.

Thanks for the advice with the advanced notebook example. I think I will start with this one and to be frank I'm having problems installing the scikits.odes on my work's windows machine but will try it at home on my fedora system.

Looking at the notebook I couldn't find the .cti file that you're using there to get a grip on how to implement my system to this problem and to run it myself. Can you help me with that one? Couldn't find it in the 2.4 repo on github neither.

Thanks again and best regards

Kurth

[Ray Speth]

Hi Kurth,

The input files are stored in the Git repository that contains all of the Jupyter examples. You can find the input file for this specific one here: https://github.com/Cantera/cantera-jupyter/tree/main/reactors/data.

Regards,

Ray

[Ronny Mönnig]

Solved:

[Ronny Mönnig]

driver Code used:

import numpy as np
import cantera as ct
import matplotlib.pyplot as plt
import csv

#mech = 'Methanation_Ni_b.yaml'
mech = 'surfaceNi.yaml'


# import the model for gas-Si-N interface and adjacent bulk phases
gas_NIO_interface = ct.Interface(mech, 'NI_SURFACE')
print(gas_NIO_interface.adjacent)
gas = gas_NIO_interface.adjacent['gas']

# Set the initial conditions
T0 = 600  # K
p0 = ct.one_atm
gas.TPX = T0, p0, "CO2:0.25, H2:0.75"
gas_NIO_interface.TP = T0, p0
D = 2.5e-3  # diameter [m]
Ac = np.pi * D**2 / 4  # cross section [m²]
u0 = 1.5/100  # initial velocity [m s⁻1]

reactor = ct.FlowReactor(gas)
reactor.area = Ac
reactor.mass_flow_rate = gas.density * u0 * Ac
reactor.energy_enabled = False

rsurf = ct.ReactorSurface(gas_NIO_interface, reactor)
net = ct.ReactorNet([reactor])
soln = ct.SolutionArray(gas, extra=['x', 'speed', 'surf_coverages'])
#kN = gas_si_n_interface.kinetics_species_index('N(D)')
#kSi = gas_si_n_interface.kinetics_species_index('Si(D)')


output_data = []
while net.distance < 10.0:
    #print(net.distance, rsurf.coverages)
    net.step()
    wdot = rsurf.kinetics.net_production_rates
    soln.append(TDY=reactor.thermo.TDY,
                x=net.distance,
                speed=reactor.speed,
                surf_coverages=rsurf.coverages)
                #N_dep=wdot[kN],
                #Si_dep=wdot[kSi])
   
    #output_data.append([net.distance,soln.speed] + list(gas.X)+ list(rsurf.coverages))
    output_data.append([net.distance] + list(gas.X)+ list(rsurf.coverages))

output_filename = 'surfpfr1_out.csv'
with open(output_filename, 'w', newline="") as outfile:
    writer = csv.writer(outfile)
    writer.writerow(['Distance (mm)'] + gas.species_names + gas_NIO_interface.species_names)
    writer.writerows(output_data)

print("Results saved to '{0}'".format(output_filename))