Shock & Detonation Toolbox - Matlab Thermo Fitting Program
We have implemented a Matlab demonstration program to fit tabular thermodynamic data with a constrainted least squares fit. With this program, output for direct use in a Cantera mechanism is created. Additionally, output in the CHEMKIN format is created. Thus, a user may creat new mechanisms or append this output to existing mechanisms, replacing or creating thermodynamic coefficients for a species.
Numerical Solution Methods for Shock and Detonation Jump Conditions
Documentation discussing algorithm used for this demo, See specifically the "Thermodynamic Property Polynomial Representation" appendix.
The Matlab program is run from the script demo interp thermo.m, which then calls the function interp
thermo, located in interp thermo.m. Input for the program is the LINGRAF thermodynamic data for
the example molecule 2-Butenal and user-defined values input at the Matlab prompt. The thermodynamic
input data is given a Matlab matrix in the script twobutenal.m. The output from the program fitting program
includes the coefficients, optimization diagnostics, goodness of fit measures, plots showing
the fits, and an output file, CH3CHCHCHO output , in the CHEMKIN format that can be used in Cantera
with the ck2cti utlility to create the .cti format. Also, an Output file in the .cti format is created and
can be manually pasted into a .cti file.
By independently obtaining thermodynamic
data from molecular modeling and statistical mechanics software such as LINGRAF, one may easily adapt
this Matlab program for constructing fits for any desired species.
The user-defined input for the fitting program includes the temperature ranges, the standard state enthalpy
of formation and standard state entropy, the species name, and the species molecular composition.
The molecular composition is specified by typing at the Matlab prompt at run time. For this example, there are 3 types of atoms, four carbon, six hydrogen, and one oxygen atom so the input is "3 'Enter' C 'Enter' 4 'Enter' H 'Enter' 6 'Enter' O 'Enter' 1 'Enter".
Example fits, for two-butenal appear below, with Tmin = 200 K, Tmid = 2500
K, and Tmax = 5000 K. The piecewise function was fitted using the constrained least-squares method in Matlab.
All three piecewise functions for the non-dimensional properties,
cp(T)/R, h(T)/RT, and s0(T)/R (pressure independent entropy), were simultaneously optimized and constrained for the best overall fit.
The constraints include continuity for h(T)/RT and s0(T)/R and continuity for both cp(T)/R and its
derivatives. Also, the boundary point constraints for cp(T)/R and h(T)/RT are included.
Figure: Comparison of Properties for 2-Butenal (CH3CHCHCHO) calculated from the statistical mechanics
representation (points) to the piecewise polynomial fit (solid lines).