Author Archives: yanggao

Numerical approximation of derivatives

Suppose we are numerically approximating the second-order derivative \frac{d^2y}{dx^2}

A common way is to use finite difference scheme which uses Taylor expansion.
1. Perform Taylor expansion for y_{i+1}, y_{i}, y_{i-1} .etc
2. Assemble the Taylor series and cancel the high-order terms to construct the desired-order scheme, i.e. finding a set of coefficients which makes the summation of high-order terms to be zero.

For example, the second-order central difference scheme:
\frac{d^2y}{dx^2} \approx \frac{y_{i+1}-y_{i-1}}{\Delta x^2}

An alternative way to approximate the derivative.
1. Express y as y=a+bx+cx^2
We can see c=\frac{d^2y}{dx^2}, so the objective is to find the coefficient c.

2. Use the current data to fit this second-order polynomial (a parabola, actually a curve fitting problem)
y_1 = a+bx_1+cx_1^2
y_2 = a+bx_2+cx_2^2
y_3 = a+bx_3+cx_3^2
...
y_n = a+bx_n+cx_n^2

3. Solve this fitting problem
Since we have three unknown coefficients a, b, and c, we only need three data points. However, we can choose more data points to fit this curve, which is actually least square fitting. Now, let us write the above equations in a general vector form: Ax=B
where

    \[ A= \begin{bmatrix} 1 & x_{1} & x_{1}^2 & \\ 1 & x_{2} & x_{2}^2 & \\ 1 & x_{3} & x_{3}^2 & \\ \hdotsfor{3} \\ 1 & x_{n} & x_{n}^2 & \\ \end{bmatrix} , \quad x= \begin{bmatrix} a\\ b\\ c\\ \end{bmatrix} , \quad \textrm{and} \quad B= \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \dots \\ y_n\\ \end{bmatrix} . \]

This is an over-determined system if n>3, which we can approximate it with least square fitting.
See https://en.wikipedia.org/wiki/Overdetermined_system

4. Solution
The approximated/fitted solution is: x=(A^TA)^{-1}A^TB, where the second-order derivative \frac{d^2y}{dx^2}=c=x[3]

A Python function to extract reaction information from Cantera

GitHub repo link: https://github.com/yanggaome/Python-functions-based-on-Cantera

Author: Yang Gao

Email: yanggao.me@gmail.com

Extract reaction information from Cantera functions

including: reaction rate parameters, third-body species indices and efficiencies, fall-off parameters

Reaction types supported

Support simple reactions, third-body reactions, fall-off reactions, and PLOG reactions.

Do not support SRI, CHEBY, and chemically activated reactions so far.

How to call the function

reacget_reaction_info(g, i)

g is the initialized chemistry

e.g. g = ct.Solution(‘h2.cti’)

i the reaction index – 1 (Python is zero based)

funtion returns a structure reac

Reaction rate information in reac

reac.RA, reac.RB, reac.RE : forward temperature ABE factors

RE is converted and divided by RU, so

kf = AT**Bexp(-E/T)

RA is kmole based, and converted to mole based in the test_get_reaction_info.py based on reaction order

Make sure you are using the correct unit

Third body information

reac.ITHB: number of enhanced third-body species, type: int

reac.NKTB: indices of enhanced third-body species, type: list

reac.AIK: efficiencies of enhanced third-body species, type: list

Fall off reaction information

reac.Fall: first 3-> low pressure limit ABE, 4th-7th: alpha, T3, T1, T2 (T2 is only for 7 parameters Troe)

Reaction type flag (default value is false)

reac.isReversible = False

reac.isThirdbody = False

reac.isFalloff = False

reac.isChemical = False

reac.isPLOG = False

reac.isSimple = False

reac.isLindemann = False

reac.isTroe = False

reac.isTroe6 = False

reac.isTroe7 = False

Instructions for using reduced kinetic models in CHEMKIN

Updates for ANSYS Chemkin 17.0:
As Chemkin (Reaction Design) has been a part of ANSYS, a new version, ANSYS Chemkin 17.0, has been released.
To compile the user-programmed subroutines for ANSYS Chemkin 17.0, Visual Studio 2012 together with Intel Parallel Studio XE 2015 are needed.
The versions of Visual Studio and Intel compiler are very important. Otherwise, you will have some errors in compiling your user-defined subroutines.
======
This is an instruction for using reduced combustion kinetic models developed by Prof. Tianfeng Lu.
To download the reduced models, you can go to:
http://www.engr.uconn.edu/~tlu/mechs/mechs.htm
You can download the reaction model package from the above website, including chem.inp, therm.dat, tran.dat, and a reaction rate subroutine ckwyp.f (FORTRAN code).
======
1, Configure the compiler environment

The following configuration based on:
Windows 64-bit system
CHEMKIN is installed in C:\Program Files\Reaction\chemkin15112_win64
======
Compiler and VS: Intel Parallel Studio XE 2015 with Visual Studio 2012
Intel Fortran Compiler is included in Intel Parallel Studio XE 2015

======
Set up the environment variable FLEXLM_ARCH for the compiler:
Open the Intel Fortran compiler command window, type: SET FLEXLM_ARCH=win64
Then check it, type: ECHO %FLEXLM_ARCH%
It should return: win64

Before you make any changes, I strongly suggest to back up three folders in the installed folder of CHEMKIN
C:\Program Files\Reaction\chemkin15112_win64\

They are: bin, lib and user_routines.

2, Revise the input file: chem.inp

Add a keyword “USRPROD” in the line of REACTIONS as follows
REACTIONS MOLES CAL/MOLE USRPROD

This is to tell CHEMKIN to use user-provide reaction rate subroutine.

3, Revise the cklib_user_routines.f file

Backup cklib_user_routines.f
cklib_user_routines.f file is located in the installed folder of CHEMKIN
C:\Program Files\Reaction\chemkin15112_win64\user_routines.

Revise the subroutine CKUPROD in cklib_user_routines.f, which is the user supplied subroutine to calculate species production rate using T (temperature) and C (species concentration).
======
SUBROUTINE CKUPROD (LOUT, NKK, KTYP, T, C, ICKWRK, RCKWRK, RPROD, IFLAG)
IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N)
INCLUDE ‘user_routines_interface.inc’
DIMENSION ICKWRK(*), RCKWRK(*), C(*), T(*)
C In the original codes, there is no Y, add Y(NKK) here
DIMENSION RPROD(NKK), KTYP(NKK), Y(NKK)
PARAMETER (LIPAR=1000, LRPAR=1000, KMAX=100, IMAX=500)
DIMENSION IPAR(LIPAR), RPAR(LRPAR), WT(KMAX)
C
C Obtain pressure from T and C
CALL CKPC (T, C, ICKWRK, RCKWRK, P)
C Convert C to Y
CALL CKCTY(C, ICKWRK, RCKWRK, Y)

C Compute RPROD, moles/(cm**3*sec), using P, T, Y
CALL CKWYP_USER(P,T(1),Y,ICKWRK, RCKWRK,RPROD)

C
RETURN
END
======
Copy the code in the ckwyp.f file for the specific reduced model to the end of cklib_user_routines.f.

Because CKWYP and CKWYR are default subroutine names in CHEMKIN, to avoid compiling errors, you should change the name “CKWYP” to “CKWYP_USER”. If “CKWYR” exists, you should also change it to “CKWYR_USER”.

SUBROUTINE CKWYP_USER(P,T,Y,ICKWRK, RCKWRK,RPROD)
SUBROUTINE CKWYR_USER(RHO,T,Y,ICKWRK, RCKWRK,RPROD)

Note that, ICKWRK, RCKWRK are not used in the CKWYP_USER. Just keep consistent, leave it as it is.

4, Recompile and link the new routines

Open the Intel Fortran compiler command window, go to the location of user_routines folder:
C:\Program Files\Reaction\chemkin15112_win64\user_routines

Type:
nmake –i –f user_routines_pc.mak

Then go to drivers_cpp folder by:
Type: cd ..\drivers_cpp
Next step will compile the specific program you want to use.
Then take aurora, which is for 0-D auto-ignition calculation, for example
Before you type this command, you can go to ..\bin folder to rename or delete the auroradll.dll, auroradll.exp and auroradll.lib. But make sure you have backup.
Then type:
nmake –i –f drivers_cpp_pc.mak ..\bin\auroradll.dll

This command will generate a new aurora program (0-D auto-ignition) for the reduced model in ..\bin, including three files (auroradll.dll, auroradll.exp, auroradll.lib).

If you want to use other applications, just replace auroradll.dllwith other xxx.dll in this command.

After the above steps, you can use the reduced model as a regular one in CHEMKIN. But for each different reduced model, you should use different cklib_user_routines.f file, recompile and link the user subroutine.

For more detailed information, refer to API.PDF document in CHEMKIN, Chapter 2: User Supplemental Programming.

Modify memory limits in ANSYS CHEMKIN17.0

ECHO %CKJAVAMEMORY%
SET CKJAVAMEMORY=-Xms4g -Xmx8g

The syntax is case-sensitive.

s: starting memory
x: maximum memory

4g, 8g: memory specification in GB unit, can also be specified in MB unit, e.g. 3584m.

For more information, it is recommended to refer to Getting_Started.pdf document, section 1.4.