ThermoRefit

author: abdurrahman.imren @ WA, USA Nov 2021

ThermoRefit corrects discontinuities found at the midpoint temperature between the low- and high-temperature-range property polynomial coefficients given in CHEMKIN-II NASA polynomial format¹. If discontinuity is found, a constrained linear least squares fitting procedure is applied to enforce C1 continuity. ThermoRefit has been verified with the chemical mechanisms listed in Imren (2022)². (View only access 👍).

Two input files are required: "chem.inp" (we read only species information) and "therm.dat". Here, thermodynamic data (see given example "therm.dat") are specified in two temperature ranges with seven coefficients for each. Standard-state molar heat capacities at constant pressure, molar enthalpy, and molar entropy are calculated using expressions

• C_p/R = a₁ + a₂T + a₃T² + a₄T³ + a₅T⁴
• H/RT = a₁ + (a₂/2)T + (a₃/3)T² + (a₄/4)T³ + (a₅/5)T⁴ + a₆/T
• S/R = a₁ ln(T) + a₂T + (a₃/2)T² + (a₄/3)T³ + (a₅/4)T⁴ + a₇ .

After the first line which begins with "THERMO ALL", temperature ranges for two sets of coefficients are specified: lowest T, midpoint T, and highest T (Tmin-common, Tmid-common, Tmax-common). These temperature values are common for all species and overridden for each specie at columns 46-76 of the specie-description line where atomic symbols and formula are defined (Tmin-specie, Tmid-specie, Tmax-specie). If any specie-temperature value is not given, corresponding common-temperature value is used.

ThermoRefit can check discontinuities at both midpoint temperature values:

1. If exists at specie mid temperature (Tmid-specie)
2. If exists at common mid temperature (Tmid-common) when Tmid-specie is different than Tmid-common.

Second test is specifically for software where only common Tmid is used and specie specific temperature inputs are ignored. Here, we check if extrapolated curves (from Tmid-specie to Tmid-common) would create artificial discontinuity. An illustration is given below.

If "refit" is decided ("atol=1e-4; rtol=1e-4" @ line 133 in "ThermoRefit.py"), Tmid-specie value is replaced by the common Tmid value. To use this feature: set "TmidRefit=True" @ line 71-72 in "ThermoRefit.py". Note that proper mechanisms (not having discontinuities at Tmid-specie) still might require "refit" procedure because of implementation differences mentioned above. This is demonstrated using an example chemical mechanism ("chem.inp" and "therm.dat" files). Running ThermoRefit with "TmidRefit=False" option for this mechanism indicates no discontinuities at Tmid-specie. However, if "True" option is activated, discontinuities become obvious (see plots in "plots" folder generated by ThermoRefit). Note that "TmidRefit=False" option is default and might not need to be changed for many cases.

Linear solver is "dgglse" from SciPy LAPACK interface. Implementation is low level. Applying equality constraints and specific index treatments of the variables at midpoint temperature values made high level coding be difficult. Contributions are welcome.

Equality Constraints:

1. At left boundary: C_p(Tmin)
2. At right boundary: C_p(Tmax)
3. Interboundary: C0 for C_p(Tmid)
4. Interboundary: C1 for C_p(Tmid)
5. Interboundary: C2 for C_p(Tmid)

For H(Tmid) and S(Tmid), only C0 continuity is applied (i.e., a₆ and a₇ coefficients are calculated, respectively.)

ThermoRefit Run Options:

ThermoRefit was written in Python 3.9 using NumPy 1.21.1, SciPy 1.7.0, and Matplotlib 3.4.2 packages. It was tested on both Windows 11 and Ubuntu 20.04.05 OSs.

• "chem.inp" and "therm.dat" are required files. For OpenFOAM output, "capitalize_letters.py" script can be used.
• As explained above, two types of discontinuity checks can be made by setting "TmidRefit=True or False" @ line 71-72 in "ThermoRefit.py" file.
• Default tolerance values are atol=1e-4 and rtol=1e-4, which can be set @ line 133 in "ThermoRefit.py". For faster runs, user can specify atol=1e-1 if solution accuracy is considered to be sufficient.
• Output files are "newtherm.dat" and "therm_duplicates.dat". Only species listed in "chem.inp" are written in "newtherm.dat". Plots of refitted species are collected in "plots" folder.

In ThermoRefit, we aim to catch user-induced errors in "therm.dat" as much as possible. Once diagnosed, we want to stop and see where we made mistakes. Tested mechanisms are listed in Imren (2022).

Features of ThermoRefit:

• while reading "chem.inp", duplicates are not allowed (we are using set() data structure)
• can read small letters and UTF-8 characters
• we read "chem.inp" in order to pull out only listed species from "therm.dat"
• if specie in "chem.inp" can't be found in "therm.dat", we raise "ERROR"
• while reading "therm.dat", duplicates are updated with the first entry
• if there exists any missing specie temperature data (Tmin, Tmax, Tmid), we fill it with given common temperature values
• if specie phase is not given or its position changed, we do not fill it with "G" - default gas phase. We raise "ERROR" to correct its position and manually check given data to ensure that neighbor temperature inputs are OK
• while reading coefficients, we strictly apply two temperature-ranges. If extra coefficient line is found, "WARNING" is raised and that line is ignored. If there exists a missing line or coefficient or can't read any coefficient, we raise "ERROR"
• if specie Tbreak >= Tmax, corresponding data is flagged to "no refit"
• OpenFOAM CK format might not accept exclamation marks (or any user entry) after column 80 in "newtherm.dat". Feel free to remove !RF comments added after refitting.

Footnotes

1. Kee, R.J., Rupley, F.M., , Miller, J.A.: CHEMKIN-II: A FORTRAN chemical kinetics package for the analysis of gas-phase chemical kinetics. Technical Report SAND89-8009, Sandia National Laboratories (1990). ↩

2. Imren, A.: A Detailed Error Quantification Analysis of Extrapolation-Based Stiff ODE Solvers for Combustion CFD, Flow, Turbulence and Combustion (2022). ↩