Shock & Detonation Toolbox - Matlab Thermo Fitting Program

Shock & Detonation Toolbox Home

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.

Matlab Program

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).