From cb64aef48f059722415c2bd8a9ea932596713aa6 Mon Sep 17 00:00:00 2001 From: TristanHehnen Date: Thu, 14 Aug 2025 18:39:59 +0200 Subject: [PATCH 1/3] Fixed typos. --- docs/conf.py | 2 +- docs/tutorials/index.rst | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/conf.py b/docs/conf.py index 33763ad..40dc1a3 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -3,7 +3,7 @@ fsp_path = os.path.join("..", "src") sys.path.insert(0, os.path.abspath(fsp_path)) # So firescipy is importable -project = "firescipy" +project = "FireSciPy" author = "Tristan Hehnen, Lukas Arnold" release = "0.1.0" diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst index d799f01..cc24535 100644 --- a/docs/tutorials/index.rst +++ b/docs/tutorials/index.rst @@ -1,7 +1,7 @@ -FireScyPy Tutorials +FireSciPy Tutorials =================== -This section contains hands-on tutorials to help you get started with FireScyPy. +This section contains hands-on tutorials to help you get started with FireSciPy. Pyrolysis Reaction Kinetics --------------------------- From e6b3905a4a228b33c8ee0b53b4e6dfd3ab78db8b Mon Sep 17 00:00:00 2001 From: TristanHehnen Date: Sat, 16 Aug 2025 13:49:57 +0200 Subject: [PATCH 2/3] Add .gitattributes: nbstripout default. --- .gitattributes | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 .gitattributes diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..21e42e9 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,2 @@ +# Strip outputs from Jupyter notebooks +*.ipynb filter=nbstripout From fb6793aeea50dca7b090bfc4af1754fe791186bf Mon Sep 17 00:00:00 2001 From: TristanHehnen Date: Sat, 16 Aug 2025 14:12:51 +0200 Subject: [PATCH 3/3] Added Jupyter notebook with hands-on examples. --- .../notebooks/FireSciPy_KAS_Demo.ipynb | 1511 +++++++++++++++++ 1 file changed, 1511 insertions(+) create mode 100644 docs/tutorials/pyrolysis/notebooks/FireSciPy_KAS_Demo.ipynb diff --git a/docs/tutorials/pyrolysis/notebooks/FireSciPy_KAS_Demo.ipynb b/docs/tutorials/pyrolysis/notebooks/FireSciPy_KAS_Demo.ipynb new file mode 100644 index 0000000..264c8c4 --- /dev/null +++ b/docs/tutorials/pyrolysis/notebooks/FireSciPy_KAS_Demo.ipynb @@ -0,0 +1,1511 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Demo of the Pyrolysis Modeling and Kinetics Computation Capabilities in FireSciPy\n", + "\n", + "This Jupyter notebook demonstrates the modeling and computation capabilites of pyrolysis kinetics in FireSciPy. It contains:\n", + "\n", + "1. A basic example modeling a single-step Arrhenius reaction under a linear temperature program.\n", + "2. A more advanced example where synthetic thermogravimetric data is generated, and activation energy ($E$) is estimated using the Kissinger–Akahira–Sunose (KAS) method.\n", + "3. Estimation of the apparent activation energy ($E_a$) from experimental TGA data using KAS.\n", + "\n", + "These examples show how FireSciPy can be used both for forward modeling and inverse parameter extraction.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import re\n", + "import sys\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import FireSciPy as fsp\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from importlib import reload # Python 3.4+\n", + "from scipy.ndimage import uniform_filter1d\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Order of plot colors\n", + "plt_colors = [\"tab:blue\", \"tab:orange\", \"tab:green\", \n", + " \"tab:red\", \"tab:purple\", \"tab:brown\", \n", + " \"tab:pink\", \"tab:gray\", \"tab:olive\", \n", + " \"tab:cyan\"]\n", + "\n", + "# Global settings for plotting\n", + "plt.rcParams.update({\n", + " 'axes.axisbelow': True, # Keep grid behind plots\n", + " 'figure.autolayout': True, # Equivalent to calling tight_layout()\n", + " 'axes.facecolor': 'white', # Prevents transparent background\n", + " 'grid.alpha': 0.6, # Makes gridlines more readable\n", + " 'font.size': 12 # Set global font size\n", + "})\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Single-Step Pyrolysis Kinetics Modelling\n", + "\n", + "In this section, the basic pyrolysis modelling capabilities of FireSciPy are introduced.\n", + "\n", + "\n", + "### Theoretical Background\n", + "----------------------\n", + "\n", + "Pyrolysis kinetics describe the rates of thermally induced decomposition processes and are commonly studied using thermal analysis techniques. Common methods are thermogravimetric analysis (TGA), Differential Scanning Calorimetry (DSC) or Microscale Combustion Calorimetry (MCC), to name a few. The process rates can be described in terms of the temperature $T$, the extent of conversion $\\alpha$ and the pressure $P$, see formula (1) below:\n", + "\n", + "\n", + "$$\n", + "\\frac{d \\alpha}{dt} = k(T) ~f(\\alpha) ~h(P)\n", + "$$\n", + "\n", + "Where $\\frac{d \\alpha}{dt}$ is the conversion rate, $k(T)$ is the reaction rate constant, $f(\\alpha)$ is the reaction model and $h(P)$ is the pressure dependence. The pressure dependence plays a role when gases can react with the decomposing sample material. During TGA experiments this can be minimized by maintaining an inert atmosphere and ensuring a sufficiently high purge gas flow to swiftly remove evolved gases which might otherwise react with the remaining sample material. Under these conditions the pressure term can than be dropped.\n", + "\n", + "Typically, the reaction rate constant $k(T)$ is expressed as an Arrhenius equation, see formula (2):\n", + "\n", + "\n", + "$$\n", + "k(T) = A ~exp \\left( -\\frac{E}{R ~T}\\right)\n", + "$$\n", + "\n", + "Where $A$ is the pre-exponential factor, $E$ the activation energy and $R$ the gas constant.\n", + "\n", + "Substituting equation (2) into (1) and neglecting $h(P)$ leads to equation (3):\n", + "\n", + "\n", + "$$\n", + "\\frac{d \\alpha}{dt} = A ~exp \\left( -\\frac{E}{R ~T}\\right) ~f(\\alpha)\n", + "$$\n", + "\n", + "Together, $A$, $E$ and $f(\\alpha)$ form the kinetic triplet.\n", + "\n", + "\n", + "More details are available in the recommendations provided by the International Confederation for Thermal Analysis and Calorimetry (ICTAC) Kinetics Committee:\n", + "\n", + "- [ICTAC Kinetics Committee recommendations for performing kinetic computations on thermal analysis data, 2011 (https://doi.org/10.1016/j.tca.2011.03.034)](https://doi.org/10.1016/j.tca.2011.03.034)\n", + "- [ICTAC Kinetics Committee recommendations for collecting experimental thermal analysis data for kinetic computations, 2014 (https://doi.org/10.1016/j.tca.2014.05.036)](https://doi.org/10.1016/j.tca.2014.05.036)\n", + "- [ICTAC Kinetics Committee recommendations for analysis of multi-step kinetics, 2020 (https://doi.org/10.1016/j.tca.2020.178597)](https://doi.org/10.1016/j.tca.2020.178597)\n", + "- [ICTAC Kinetics Committee recommendations for analysis of thermal decomposition kinetics, 2023 (https://doi.org/10.1016/j.tca.2022.179384)](https://doi.org/10.1016/j.tca.2022.179384)\n", + "\n", + "\n", + "### Conversion Modelling in FireSciPy\n", + "---------------------------------\n", + "\n", + "In this introduction, the conversion of a single-step pyrolysis reaction is modelled. As an example, a linear heating rate of $\\beta = 5~K/min$ is chosen for the temperature program. The parameters are taken from [Vyazovkin, Advanced Isoconversional Method, 1997 (https://doi.org/10.1007/BF01983708)](https://doi.org/10.1007/BF01983708).\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Reaction constants\n", + "A = 10**10 / 60 # 1/s\n", + "E = 125.4 * 1000 # J/mol\n", + "alpha0 = 1e-12 # Avoid exactly zero to prevent numerical issues (e.g., division by zero)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "At first, the linear temperature program needs to be created. Let's assume the sample temperature at the start of the experiment is at 300 K and it will end at a temperature of 750 K. The temperature is to change linearly over time according to the heating rate of $\\beta = 5~K/min$. Heating rates can be provided in either $K/min$ or $K/s$, just pass the unit as a string, and FireSciPy will internally handle the unit conversion. The recording frequency of the TGA device is assumed to be such that a spacing of $\\Delta T = 0.5~K$ is achieved. This spacing is a model of the frequency with which a TGA or similar device records the data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Define temperature program (model recording frequency ΔT during TGA experiment)\n", + "n_points = 2*450 + 1 # ΔT = 0.5 K\n", + "\n", + "# Temperatures in Kelvin\n", + "start_temp = 300\n", + "end_temp = 750\n", + "\n", + "\n", + "# Define model temperature program\n", + "beta = 5 # K/min\n", + "\n", + "\n", + "# Create the temperature program\n", + "temp_program = fsp.pyrolysis.modeling.create_linear_temp_program(\n", + " start_temp=start_temp, \n", + " end_temp=end_temp, \n", + " beta=beta, \n", + " beta_unit=\"K/min\", \n", + " steps=n_points)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "Once the temperature program is set up, the conversion can be computed. The kinetics solver needs to be provided with a couple of parameters: the temperature program (`t_array`, `T_array`), the initial conversion level (`alpha0`), the Arrhenius parameters (`A`, `E`) and the reaction model.\n", + "\n", + "FireSciPy comes with a small library of different reaction models like `'nth_order'` or `'Avrami_Erofeev'`. The values for their parameters can be provided as a dictionary.\n", + "In this example, we assume a first-order reaction model $f(\\alpha) = (1 - \\alpha)^n$ with $n~=~1$, which is commonly used in thermal decomposition modeling due to its simplicity.\n", + "Thus, the nth-order model is chosen (`'nth_order'`) with the order set to unity `{'n': 1.0}`. The result are data series of how the conversion (`alpha_sol`) changes over time (`t_sol`).\n", + "\n", + "See also: `fsp.pyrolysis.modeling.create_linear_temp_program`, `fsp.pyrolysis.modeling.solve_kinetics`\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Get time-temperature data series\n", + "time_model = temp_program[\"Time\"]\n", + "temp_model = temp_program[\"Temperature\"]\n", + "\n", + "# Compute conversion\n", + "t_sol, alpha_sol = fsp.pyrolysis.modeling.solve_kinetics(\n", + " t_array=time_model, \n", + " T_array=temp_model, \n", + " A=A, \n", + " E=E,\n", + " alpha0=alpha0,\n", + " R=fsp.constants.GAS_CONSTANT,\n", + " reaction_model='nth_order',\n", + " model_params={'n': 1.0})\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's plot the conversion against the sample temperature to see the result. The code below will produce a plot showing the conversion as a function of temperature. Users are encouraged to run the example locally to see the result." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Plot conversion\n", + "plt.plot(temp_model, alpha_sol,\n", + " label=f\"{beta} K/min\")\n", + "\n", + "\n", + "# Plot meta data\n", + "plt.title(\"Conversion of a Single-Step Pyrolyis Reaction\")\n", + "plt.xlabel(\"Sample Temperature / K\")\n", + "plt.ylabel(\"Conversion ($\\\\alpha$) / -\")\n", + "\n", + "plt.tight_layout()\n", + "plt.legend()\n", + "plt.grid()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For more details, see the [FireSciPy documentation](https://firedynamics.github.io/FireSciPy/)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## KAS Activation Energy Analysis from Simulated TGA\n", + "\n", + "In this example, the conversion data is generated by leaveraging the modelling capabilities of FireSciPy. An example with experimental data is available too.\n", + "\n", + "### Theoretical Background\n", + "----------------------\n", + "\n", + "The conversion rate $\\frac{d \\alpha}{dt}$ of a single-step pyrolysis reaction can be described as presented in equation (4):\n", + "\n", + "\n", + "$$\n", + "\\frac{d \\alpha}{dt} = k(T) ~f(\\alpha)\n", + "$$\n", + "\n", + "and depends on the rate constant $k(T)$ and the reaction model $f(\\alpha)$. Commonly, the rate constant is expressed as an Arrhenius equation (5), where $A$ is the pre-exponential factor, $E$ the activation energy and $R$ the gas constant.\n", + "\n", + "\n", + "$$\n", + "\\frac{d \\alpha}{dt} = A ~exp \\left( -\\frac{E}{R ~T}\\right) ~f(\\alpha)\n", + "$$\n", + "\n", + "Together, $A$, $E$ and $f(\\alpha)$ form the kinetic triplet.\n", + "\n", + "Kinetics computations aim to deduce the elements of the kinetic triplet from experimental data. Here, specifically the determination of the activation energy $E$ is determined. It can be determined using an isoconversional method. These method are based on the observation that for a given level of conversion the reaction rate depends solely on the temperature. Taking the logarithmic derivative of equation (5) leads to equation (6). \n", + "\n", + "\n", + "$$\n", + "\\left[ \\frac{\\partial \\ln \\!\\left(\\dfrac{d\\alpha}{dt}\\right)}{\\partial T^{-1}} \\right]_{\\alpha}\n", + "=\n", + "- \\frac{E_\\alpha}{R}\n", + "$$\n", + "\n", + "At a constant conversion reaction model $f(\\alpha)$ is constant and has no effect. For this reason isoconversional methods are sometimes called \"model-free\". However, this does not mean that $f(\\alpha)$ can be ignored with respect to the kinetics triplet. \n", + "\n", + "For a constant heating rate program, the equation has no analytical solution and the temperature integral needs to be approximated. This is accomplished by recording data from experiments, for example thermogravimetric analysis (TGA). Data recorded at multiple different heating rates is used to determine the respective temperatures for a given level of conversion. The Kissinger–Akahira–Sunose (KAS) is a common isoconversional method, shown in equation (7), with $B=2$ and $C=1$. It has a higher accuracy compared to other well known methods like the Ozawa–Flynn–Wall (OFW) method. \n", + "\n", + "\n", + "$$\n", + "\\ln\\!\\left( \\frac{\\beta_i}{T_{\\alpha,i}^{B}} \\right) = \\text{Const} - C \\frac{E_{\\alpha}}{R T_{\\alpha}}\n", + "$$\n", + "\n", + "The KAS method is further improved by Starink, setting $B=1.92$ and $C=1.0008$. The KAS method is available in FireSciPy and the Starink improvement is the default setup.\n", + "\n", + "\n", + "More details are available in the recommendations provided by the International Confederation for Thermal Analysis and Calorimetry (ICTAC) Kinetics Committee:\n", + "\n", + "- [ICTAC Kinetics Committee recommendations for performing kinetic computations on thermal analysis data, 2011 (https://doi.org/10.1016/j.tca.2011.03.034)](https://doi.org/10.1016/j.tca.2011.03.034)\n", + "- [ICTAC Kinetics Committee recommendations for collecting experimental thermal analysis data for kinetic computations, 2014 (https://doi.org/10.1016/j.tca.2014.05.036)](https://doi.org/10.1016/j.tca.2014.05.036)\n", + "- [ICTAC Kinetics Committee recommendations for analysis of multi-step kinetics, 2020 (https://doi.org/10.1016/j.tca.2020.178597)](https://doi.org/10.1016/j.tca.2020.178597)\n", + "- [ICTAC Kinetics Committee recommendations for analysis of thermal decomposition kinetics, 2023 (https://doi.org/10.1016/j.tca.2022.179384)](https://doi.org/10.1016/j.tca.2022.179384)\n", + "\n", + "\n", + "### Synthetic TGA Data and KAS Estimation\n", + "-------------------------------------------------------------------------\n", + "\n", + "FireSciPy provides functionalities to conduct reaction kinetics computations. In this example, the activation energy $E$ is determined, using the Kissinger-Akaira-Sunose (KAS) isoconversional method. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Reaction constants\n", + "A = 10**10 / 60 # 1/s\n", + "E = 125.4 * 1000 # J/mol\n", + "alpha0 = 1e-12 # Avoid exactly zero to prevent numerical issues (e.g., division by zero)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "FireSciPy introduces an internal data structure for the pyrolysis computations. The goal of this data structure is to keep all the data organised. Furthermore, the different methods of FireSciPy are aware of the structure and can thus easily call the data they need and store their own results. \n", + "\n", + "In this particular example, the reading and processing of experimental data is not discussed. Data series will be generated from modelling, therefore some steps are skipped. The generated data will be placed manually in the data structure in the locations where the processed experimental data would end up. A dedicated example with experimental data is available too.\n", + "\n", + "Below, the skeleton of this data structure is initialised. It also allows to track some meta-data, like which material was investigated, what device was used or who performed the experiments. Furthermore, a dictionary with information on the recorded quantity (signal) is provided. This allows to address mass data from TGA, i.e. `{\"name\": \"Mass\", \"unit\": \"mg\"}`, or heat flow data from DSC, i.e. `{\"name\": \"HeatFlow\", \"unit\": \"W/g\"}`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Initialise data structure for the kinetics assessment\n", + "data_structure = fsp.pyrolysis.kinetics.initialize_investigation_skeleton(\n", + " material=f\"Unobtainium\", \n", + " investigator=\"John Doe, Miskatonic University\", \n", + " instrument=\"ACME TGA 9000\", \n", + " date=\"Stardate: 42.69\", \n", + " notes=\"It has a colour out of space\",\n", + " signal={\"name\": \"Mass\", \"unit\": \"mg\"})\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Basic settings for the temperature programs are defined below. For this example multiple heating rates are chosen and stored in a Python dictionary. The dictionary allows to easily keep track of temperature programs. \n", + "\n", + "Keep in mind, as per the ICTAC recommendations, the range of heating rates should cover at least a factor of 10 from the lowest to the highest. Furthermore, at the very least 3 heating rates should be used for isoconversional computations. Better would be 5 or more. Furthermore, basic isoconversional methods like KAS are built on the assumption that the heating rate is perfectly linear. This is captured here and the heating rates are thus labelled \"nominal\". \n", + "\n", + "The temperature programs all start at $300~K$ and end at $750~K$. Means are provided to adjust the sampling frequency (temperature spacing $\\Delta T$). A few different frequencies are predefined and can easily be commented in to play around with.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Nominal heating rates\n", + "heating_rates = {\n", + " \"1Kmin\": 1,\n", + " \"3Kmin\": 3,\n", + " \"10Kmin\": 10,\n", + " \"15Kmin\": 15,\n", + " \"30Kmin\": 30,\n", + " \"45Kmin\": 45,\n", + " \"60Kmin\": 60\n", + "}\n", + "\n", + "# Data for temperature program (Kelvin)\n", + "T_start = 300\n", + "T_end = 750\n", + "\n", + "# Set resolution for the temperature spacing (ΔT)\n", + "# resolution_factor = 0.5 # ΔT = 1.0 K\n", + "resolution_factor = 1 # ΔT = 0.5 K\n", + "# resolution_factor = 2 # ΔT = 0.25 K\n", + "# resolution_factor = 5 # ΔT = 0.1 K\n", + "\n", + "# Number of points to be used in the arrays to be generated below\n", + "n_points = (int(900 * resolution_factor)) + 1\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the following code block, the individual temperature programs are created, based on the nominal heating rates. This is accomplished with the function `fsp.pyrolysis.modeling.create_linear_temp_program`. \n", + "\n", + "Looping over the above dictionary, a temperature program is created for each nominal heating rate. The time and temperature data series of each program are then extracted. They are used in the pyrolysis solver `fsp.pyrolysis.modeling.solve_kinetics`, together with the n-th order reaction model. The computed conversion ($\\alpha$) is added to the temperature program DataFrame. The column headers of the DataFrame are renamed to match the expected values. Nominal values of the heating rates are store in the data structure, as well as the DataFrame. This information would typically come from the experiments and would then be added differently, see other example.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Modelling conversion and store in data structure\n", + "for hr_label in heating_rates:\n", + " # Compute model heating rates\n", + " beta = heating_rates[hr_label]\n", + " hr_model = fsp.pyrolysis.modeling.create_linear_temp_program(\n", + " start_temp=T_start, \n", + " end_temp=T_end, \n", + " beta=beta, \n", + " beta_unit=\"K/min\", \n", + " steps=n_points)\n", + " \n", + " # Get temperature program\n", + " time_model = fsp.utils.series_to_numpy(hr_model[\"Time\"])\n", + " temp_model = fsp.utils.series_to_numpy(hr_model[\"Temperature\"])\n", + " # Get overview over temperature resolution\n", + " ΔT = temp_model[1] - temp_model[0]\n", + " print(f\"Temperature resolution ({hr_label}): ΔT = {ΔT} K\")\n", + " # Compute conversion for decelerating reaction (n-th order)\n", + " t_sol, alpha_sol = fsp.pyrolysis.modeling.solve_kinetics(\n", + " t_array=time_model, \n", + " T_array=temp_model, \n", + " A=A, \n", + " E=E,\n", + " alpha0=alpha0,\n", + " R=fsp.constants.GAS_CONSTANT,\n", + " reaction_model='nth_order',\n", + " model_params={'n': 1.0}\n", + " )\n", + " \n", + " # Convert to DataFrame\n", + " hr_model = pd.DataFrame(hr_model)\n", + " # Add new column with conversion data\n", + " hr_model[\"alpha\"] = pd.Series(alpha_sol, index=hr_model.index)\n", + " # Rename column headers to match expected values\n", + " hr_model.rename(columns={\n", + " 'Time': 'Time', \n", + " 'Temperature': 'Temperature_Avg', \n", + " 'alpha': 'Alpha'}, \n", + " inplace=True)\n", + " \n", + " # Ensure nested structure to store data, because model \n", + " # data is added manually here in this example and no data\n", + " # processing from experimental TGA curves takes place\n", + " keys = [\"experiments\", \"TGA\", \"constant_heating_rate\", hr_label]\n", + " heating_rate = fsp.pyrolysis.kinetics.ensure_nested_dict(data_structure, keys)\n", + " \n", + " # Add nominal heating rates\n", + " nominal_beta = {\"value\": beta, \"unit\": \"K/min\"}\n", + " heating_rate[\"set_value\"] = nominal_beta\n", + "\n", + " # Add the modeled conversion data to the data structure\n", + " heating_rate[\"conversion\"] = hr_model\n", + " \n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, the desired conversion fractions need to be specified. They are necessary to ensure that the linear fit of the KAS method compares appropriate temperatures for a given level of conversion. During the experiments the data is recorded with a frequency that can be adjusted at the device. The changes that individual data points match up across many heating rates and for all desired conversion levels is very slim. Thus, interpolation of the input data is necessary. This interpolation is conducted with `fsp.pyrolysis.kinetics.compute_conversion_fractions`. As parameters, the data structure, the desired points for the analysis and the setup need to be provided. There is also an option to compute the fractions for selected temperature programs or all that are available.\n", + "\n", + "Commonly, the conversion fractions range between (0.05 < $\\alpha$ < 0.95) or (0.1 < $\\alpha$ < 0.9). The reason being, that at the ends changes in are small over long times and the experimental noise is large. This leads to spurious results. However, the ends should not be rejected flat out. The user needs to assess how significant fluctuations are to not neglect, for example, initial reactions. In this example, a wider range is chosen: (0.01 < $\\alpha$ < 0.99). Since the input data comes from modelling, the noise is low for the provided settings. It will be highlighted below that noise increases towards the end.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Define conversion fractions where to evaluate the activation energy\n", + "# Note: commonly, they range between (0.05 < α < 0.95) or (0.1 < α < 0.9)\n", + "# conversion_levels = np.linspace(0.05, 0.95, 37) # Δα = 2.5\n", + "conversion_levels = np.linspace(0.01, 0.99, 99) # Δα = 1.0\n", + "\n", + "# Compute conversion fractions across all \"experiments\"\n", + "fsp.pyrolysis.kinetics.compute_conversion_levels(\n", + " data_structure, \n", + " desired_levels=conversion_levels, \n", + " setup=\"constant_heating_rate\", \n", + " condition=\"all\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The isoconversional activation energy can now be computed, using the KAS implementation `fsp.pyrolysis.kinetics.compute_Ea_KAS`. It simply needs to be provided with the data structure, all the necessary information should now be in the correct place.\n", + "\n", + "By default, the KAS method implemented in FireSciPy uses the Starink improvement with the parameters $B=1.92$ and $C=1.0008$. For the classical KAS method, they can be changed to $B=2$ and $C=1$. See also: \n", + "- [Starink, The determination of activation energy from linear heating rate experiments: a comparison of the accuracy of isoconversion methods, 2003 (https://doi.org/10.1016/S0040-6031(03)00144-8)](https://doi.org/10.1016/S0040-6031(03)00144-8)\n", + "- [ICTAC Kinetics Committee recommendations for performing kinetic computations on thermal analysis data, 2011 (https://doi.org/10.1016/j.tca.2011.03.034)](https://doi.org/10.1016/j.tca.2011.03.034)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Compute the activation energy using the KAS method\n", + "fsp.pyrolysis.kinetics.compute_Ea_KAS(data_structure, B=1.92, C=1.0008)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, the results can be plotted: the development of the activation energy over the conversion. In gray, the first and last $5\\%$ of the conversion are highlighted. Even using model data as input, it is observable that the noise increases towards the ends. Feel free to play with different sampling rates to see how they affect the result.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "\n", + "# Get the Ea and convert to kJ/mol.\n", + "Ea_results_KAS = data_structure[\"experiments\"][\"TGA\"][\"Ea_results_KAS\"]\n", + "Ea = Ea_results_KAS[\"Ea\"]/1000\n", + "Ea_avg = np.average(Ea)\n", + "\n", + "# Plot the Ea against conversion.\n", + "conv = Ea_results_KAS[\"Conversion\"]\n", + "plt.scatter(conv,\n", + " Ea,\n", + " marker=\".\", s=42,\n", + " facecolors=\"none\",\n", + " edgecolors=\"tab:blue\",\n", + " label=f\"KAS, ΔT = {ΔT} K\"\n", + " )\n", + "\n", + "\n", + "# Plot target value, i.e. model input\n", + "plt.plot([0,1], [125.4, 125.4], \n", + " color=\"black\", linestyle=\"--\",\n", + " label=\"Target, E$_a$=125.40 kJ/mol\"\n", + " )\n", + "\n", + "\n", + "# Shaded areas to indicate first/last 5 % to be typically cut off\n", + "x_min = -0.025\n", + "x_max = 1.025\n", + "ax.axvspan(x_min, 0.05, color='gray', alpha=0.3)\n", + "ax.axvspan(0.95, x_max, color='gray', alpha=0.3)\n", + "\n", + "\n", + "# Plot meta data.\n", + "plt.title(f\"Activation Energy (KAS), E$_a$={Ea_avg:.2f} kJ/mol (avg.)\")\n", + "plt.xlabel(\"Conversion ($\\\\alpha$) / -\")\n", + "plt.ylabel(\"Activation Energy (E$_a$) / kJ/mol\")\n", + "\n", + "plt.xlim(left=x_min, right=x_max)\n", + "# plt.ylim(bottom=125.175, top=125.525)\n", + "plt.ylim(bottom=122.5, top=127.5)\n", + "\n", + "plt.tight_layout()\n", + "plt.legend(loc=\"upper center\")\n", + "plt.grid()\n", + "\n", + "\n", + "# # Save image.\n", + "# plot_label = f\"Ea_Estimate_KAS.png\"\n", + "# plot_path = os.path.join(plot_label)\n", + "# plt.savefig(plot_path, dpi=320, bbox_inches='tight', facecolor='w')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Tracking of Intermediate Steps and Basic Statistics\n", + "\n", + "FireSciPy keeps track of the intermediate steps during the kinetic analysis. Intermediate results are stored in the data structure and can be accessed to check the computation or share condensed results with others.\n", + "\n", + "As an example, the conversion fractions for the different heating rates can easily be plottet as demonstrated below.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Pre-select data sets for convenience\n", + "heating_rate_data = data_structure[\"experiments\"][\"TGA\"][\"constant_heating_rate\"]\n", + "\n", + "# Go over all available heating rate data\n", + "for hr_label in heating_rate_data:\n", + " # Get conversion fractions data series\n", + " conv_frac_data = heating_rate_data[hr_label][\"conversion_fractions\"]\n", + " temp = conv_frac_data[\"Temperature_Avg\"]\n", + " alpha = conv_frac_data[\"Alpha\"]\n", + " # Plot conversion against sample temperature\n", + " plt.plot(temp, alpha, \n", + " linestyle='none', marker=\".\", ms=3,\n", + " label=f\"{hr_label[:-4]} K/min\")\n", + "\n", + "\n", + "# Plot meta data\n", + "plt.title(\"Sampling Rate of Conversion Fractions\")\n", + "plt.xlabel(\"Sample Temperature / K\")\n", + "plt.ylabel(\"Conversion ($\\\\alpha$) / -\")\n", + "\n", + "plt.xlim(left=480, right=730)\n", + "\n", + "plt.tight_layout()\n", + "plt.legend()\n", + "plt.grid()\n", + "\n", + "\n", + "# # Save image.\n", + "# plot_label = f\"ConversionFractions_KAS.png\"\n", + "# plot_path = os.path.join(plot_label)\n", + "# plt.savefig(plot_path, dpi=320, bbox_inches='tight', facecolor='w')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "FireSciPy keeps also track of the individual fits conducted for the conversion faction, by the KAS method. This includes the parameters of the linear fit for each conversion level as well as the respecpective coordinates involved in the fit. These intermediate results are stored together with the primary results of the activation energy, as shown below.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the results of the activation energy\n", + "Ea_KAS = data_structure[\"experiments\"]['TGA'][\"Ea_results_KAS\"]\n", + "\n", + "# Check results\n", + "Ea_KAS.head()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This information can be nicely used to visualise the KAS process. Below, data points and their respective fit are shown for a selection of conversion levels.\n", + "\n", + "Note: Using the Starink improvement, the linear fit is taking place in the space of $ln\\left(\\frac{\\beta}{T^{B}}\\right)$ against $\\frac{1}{T}$, with $B=1.92$. In case a different value for $B$ is chosen (see above) the space is adjusted accordingly. The user has to adjust the axis label accordingly.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Pre-select data sets for convenience\n", + "Ea_KAS = data_structure[\"experiments\"]['TGA'][\"Ea_results_KAS\"]\n", + "\n", + "# Define settings for the plots of the fits\n", + "marker_size = 42\n", + "fit_line = \":\"\n", + "fit_alpha = 0.6\n", + "fit_color = \"black\"\n", + "\n", + "# Choose conversion levels for the plot\n", + "conversion_levels = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]\n", + "\n", + "# Initialise data collection\n", + "conv_temps = np.zeros((len(conversion_levels), len(heating_rates)))\n", + "conv_levels = np.zeros((len(conversion_levels), len(heating_rates)))\n", + "\n", + "for conv_idx, level in enumerate(conversion_levels):\n", + " level_idx = np.abs(Ea_KAS[\"Conversion\"] - level).argmin()\n", + " # Get the first three data points to match previous plot\n", + " conv_temps[conv_idx,:] = np.asarray(Ea_KAS.loc[:,\"x1\":\"x7\"].iloc[level_idx])\n", + " conv_levels[conv_idx,:] = np.asarray(Ea_KAS.loc[:,\"y1\":\"y7\"].iloc[level_idx])\n", + " # Indicate fits\n", + " m_fit = Ea_KAS[\"m_fit\"].iloc[level_idx]\n", + " b_fit = Ea_KAS[\"b_fit\"].iloc[level_idx]\n", + " x_fit = [conv_temps[conv_idx,:][0], conv_temps[conv_idx,:][-1]]\n", + " y_fit = [fsp.utils.linear_model(x_fit[0], m_fit, b_fit),\n", + " fsp.utils.linear_model(x_fit[-1], m_fit, b_fit)]\n", + " \n", + " if conv_idx == 0:\n", + " plot_label = \"Fit\"\n", + " else:\n", + " plot_label = \"_none_\"\n", + " \n", + " plt.plot(x_fit, y_fit,\n", + " linestyle=fit_line,\n", + " alpha=fit_alpha,\n", + " color=fit_color,\n", + " label=plot_label)\n", + "\n", + "\n", + "hr_labels = list(heating_rates)\n", + "# Plot data points by heating rate\n", + "for idx in range(len(conv_temps.T)):\n", + " # Get colour for data series, i.e. heating rate\n", + " plot_colour = plt_colors[idx]\n", + " \n", + " hr_label = hr_labels[idx]\n", + " # Plot data points\n", + " plt.scatter(\n", + " conv_temps.T[idx],\n", + " conv_levels.T[idx],\n", + " marker='o', s=marker_size,\n", + " facecolors='none',\n", + " edgecolors=plot_colour,\n", + " label=f\"{hr_label[:-4]} K/min\")\n", + "\n", + "\n", + "# Plot meta data.\n", + "plt.title(\"Assess Linear Fits of KAS Method\")\n", + "plt.xlabel(\"1/T\")\n", + "plt.ylabel(\"ln($\\\\beta$/T$^{1.92}$)\")\n", + "\n", + "plt.xlim(left=0.00141, right=0.00197)\n", + "plt.ylim(bottom=-12.6, top=-7.9)\n", + "\n", + "plt.tight_layout()\n", + "plt.legend()\n", + "plt.grid()\n", + "\n", + "# # Save image.\n", + "# plot_label = f\"Ea_Estimate_KAS_Fit.png\"\n", + "# plot_path = os.path.join(plot_label)\n", + "# plt.savefig(plot_path, dpi=320, bbox_inches='tight', facecolor='w')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "FireSciPy provides limited means of statistics to assess the quality of the linear regression fits during the KAS computation. This information is stored together with the primary and intermediate results of the activation energy computation.\n", + "\n", + "Specifically, the coefficient of determination ($R^2$) and the root mean square error (RMSE) are available. They can be plotted together as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Pre-select data sets for convenience\n", + "Ea_KAS = data_structure[\"experiments\"]['TGA'][\"Ea_results_KAS\"]\n", + "\n", + "# Reduce number auf data points by using every n-th\n", + "nth = 1\n", + "\n", + "markersize = 42\n", + "\n", + "fig, ax = plt.subplots()\n", + "\n", + "# Plot R² statistics\n", + "x_data = np.asarray(Ea_KAS[\"Conversion\"])[::nth]\n", + "y_data = np.asarray(Ea_KAS[\"R_squared\"])[::nth]\n", + "plt.scatter(x_data,\n", + " y_data,\n", + " marker='.', s=markersize,\n", + " facecolors='none',\n", + " edgecolors='tab:blue',\n", + " label=f\"R² avrg.: {np.average(y_data):.6f}\")\n", + "\n", + "\n", + "# Plot RMSE statistics\n", + "y_data = np.asarray(Ea_KAS[\"RMSE\"])[::nth]\n", + "plt.scatter(x_data,\n", + " y_data,\n", + " marker='.', s=markersize,\n", + " facecolors='none',\n", + " edgecolors='tab:orange',\n", + " label=f\"RMSE avrg.: {np.average(y_data):.6f}\")\n", + "\n", + "\n", + "# Shaded areas to indicate first/last 5 % to be typically cut off\n", + "x_min = -0.025\n", + "x_max = 1.025\n", + "ax.axvspan(x_min, 0.05, color='gray', alpha=0.3, label=\"5%\")\n", + "ax.axvspan(0.95, x_max, color='gray', alpha=0.3)\n", + "\n", + "\n", + "# Plot meta data.\n", + "plt.title(\"Statistics of the KAS Fits\")\n", + "plt.xlabel(\"Conversion, $\\\\alpha$ / -\")\n", + "plt.ylabel('Arbitrary Units / -')\n", + "\n", + "plt.xlim(left=-0.025, right=1.025)\n", + "plt.ylim(bottom=-0.05, top=1.05)\n", + "\n", + "plt.tight_layout()\n", + "plt.legend()\n", + "plt.grid()\n", + "\n", + "# # Save image.\n", + "# plot_label = \"Ea_Estimate_KAS_Statistics.png\"\n", + "# plot_path = os.path.join(plot_label)\n", + "# plt.savefig(plot_path, dpi=320, bbox_inches='tight', facecolor='w')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For this example the results are a bit boring. Using experimental data the statistics are more useful. \n", + "\n", + "For more details, see the [FireSciPy documentation](https://firedynamics.github.io/FireSciPy/).\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Activation Energy Analysis (KAS) from TGA Experiments\n", + "\n", + "In this example, the conversion data is derived from **experimental** thermogravimetric analysis (TGA) data.\n", + "\n", + "### Theoretical Background\n", + "--------------------------\n", + "\n", + "The apparent activation energy $E_a$ is estimated using the well-established isoconversional method known as Kissinger–Akahira–Sunose (KAS). For background on the theory and mathematical formulation, see the corresponding section in the [synthetic conversion data example](mock_kinetics_computation).\n", + "\n", + "We make the distinction here, between the true activation energy $E$ and the apparent activation energy $E_a$. \n", + "\n", + "In micro-scale experiments such as TGA, it is generally not possible to deduce the fundamental reaction mechanisms or individual reaction steps involved in the decomposition. What appears as a single peak in the mass loss rate may in fact result from multiple overlapping reactions, with one step acting as rate-limiting and dominating the overall shape.\n", + "\n", + "As such, the estimated activation energy from these experiments — $E_a$ — is considered apparent, reflecting the net behavior rather than a single elementary step.\n", + "\n", + "The true activation energy is only definitively known in modeling scenarios, where the decomposition process is explicitly prescribed and controlled by the user (e.g., synthetic data or known reaction schemes).\n", + "\n", + "See the ICTAC recommendations below for a detailed discussion of this distinction.\n", + "\n", + "To apply the KAS method meaningfully, experimental data must be recorded under **multiple temperature programs**, either with **linear heating rates** or under **isothermal conditions**. The ICTAC Kinetics Committee recommends using at least three, preferably five, distinct datasets. If linear heating is used, the lowest and highest heating rates should differ by at least a factor of 10.\n", + "\n", + "To ensure high-quality experimental data, a few key aspects must be considered:\n", + "\n", + "- **Thermophysical effects of the sample material**: These influence the heat transfer within the sample and can distort the measured reaction rate. Their impact can be minimized by reducing the sample mass and repeating experiments with decreasing amounts until the results (e.g., mass loss rate curves) can be superimposed. This is especially important in micro-scale TGA experiments.\n", + "\n", + "- **Secondary reactions involving evolved gases**: In some materials, volatile products may react with the remaining solid. This effect can be mitigated by increasing the purge gas flow rate to carry away evolved gases more effectively. Again, comparing mass loss rate curves for different flow rates can help confirm when the flow is sufficient (i.e., when the curves overlap).\n", + "\n", + "In general, **apparent kinetic parameters should not be estimated from a single temperature program**. Multiple kinetic triplets can describe the same conversion curve equally well, making the solution **ambiguous**.\n", + "\n", + "See also:\n", + "- [Single heating rate methods are a faulty approach to pyrolysis kinetics (https://doi.org/10.1007/s13399-022-03735-z)](https://doi.org/10.1007/s13399-022-03735-z)\n", + "- [Computational aspects of kinetic analysis.: Part B: The ICTAC Kinetics Project — the decomposition kinetics of calcium carbonate revisited, or some tips on survival in the kinetic minefield (https://doi.org/10.1016/S0040-6031(00)00444-5)](https://doi.org/10.1016/S0040-6031(00)00444-5)\n", + "\n", + "More detailed guidance is available from the ICTAC Kinetics Committee:\n", + "\n", + "- [ICTAC Kinetics Committee recommendations for performing kinetic computations on thermal analysis data, 2011 (https://doi.org/10.1016/j.tca.2011.03.034)](https://doi.org/10.1016/j.tca.2011.03.034)\n", + "- [ICTAC Kinetics Committee recommendations for collecting experimental thermal analysis data for kinetic computations, 2014 (https://doi.org/10.1016/j.tca.2014.05.036)](https://doi.org/10.1016/j.tca.2014.05.036)\n", + "- [ICTAC Kinetics Committee recommendations for analysis of multi-step kinetics, 2020 (https://doi.org/10.1016/j.tca.2020.178597)](https://doi.org/10.1016/j.tca.2020.178597)\n", + "- [ICTAC Kinetics Committee recommendations for analysis of thermal decomposition kinetics, 2023 (https://doi.org/10.1016/j.tca.2022.179384)](https://doi.org/10.1016/j.tca.2022.179384)\n", + "\n", + "### Experimental TGA Data and KAS Estimation\n", + "--------------------------------------------\n", + "\n", + "In this example, the apparent activation energy $E_a$ is estimated using the KAS method based on experimental TGA data. The data was recorded for a **sample mass of 1 mg**, using multiple linear heating rates.\n", + "\n", + "Other datasets from the same study — including those with different sample masses and complementary microscale calorimeter data — are available at:\n", + "\n", + "- [Thermogravimetric And Microscale Calorimeter Data on Cast PMMA (https://doi.org/10.24355/dbbs.084-202504170956-0)](https://doi.org/10.24355/dbbs.084-202504170956-0)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example we analyse TGA data recorded at different linear heating rates. For each heating rate, three repetitions were conducted. The data is provided in CSV format, one for each repetition and experimental condition. These are plain text files with a `.txt` extension, but they follow a typical tab-delimited CSV structure.\n", + "\n", + "As in the synthetic data example, a Python dictionary is used to store all data related to this investigation. This structure is initialised using the function `initialize_investigation_skeleton`. \n", + "\n", + "Each CSV file is read into a Pandas DataFrame. If the files have meaningful names they can be directly used to label the data. For example, a file named `tga_dynamic_n2_dyn10_powder_1mg_r1.txt` corresponds to a mass of 1mg for a powdered sample, subjected to a 10 K/min heating rate and repetition 1.\n", + "\n", + "The temperatures in the files are reported in degrees Celsius and must be converted to Kelvin before further processing, as FireSciPy's kinetic computations assume Kelvin. \n", + "\n", + "Once converted, each DataFrame is added to the data structure using `add_constant_heating_rate_tga`. This function stores the input data under the appropriate labels for **heating rate** and **repetition**.\n", + "\n", + "Below is the code to read the files, adjust the temperatures, and store the results:\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Path to the experimental data csv files.\n", + "fsp_data_path = os.path.join(\"C:\\\\\", \"path\", \"to\", \"the\", \"CSV_files\")\n", + "exp_root = os.path.join(fsp_data_path, \"docs\", \"tutorials\", \"pyrolysis\", \"data\")\n", + "\n", + "\n", + "# Initialise data structure for the kinetics assessment\n", + "PMMA_data_1mg = fsp.pyrolysis.kinetics.initialize_investigation_skeleton(\n", + " material=f\"PMMA\", \n", + " investigator=\"John Doe, Miskatonic University\", \n", + " instrument=\"TGA/DSC 3+, Mettler Toledo\", \n", + " date=\"Stardate: 42.69\", \n", + " notes=\"Constant heating rates, sample mass 1mg\",\n", + " signal={\"name\": \"Mass\", \"unit\": \"mg\"})\n", + "\n", + "\n", + "for file_name in os.listdir(exp_root):\n", + " if \"powder_1mg_\" in file_name:\n", + " # print(file_name)\n", + " \n", + " # Parse metadata from file name: heating rate and repetition\n", + " name_parts = file_name.split(\"_\")\n", + " hr_value = int(name_parts[3][3:]) # e.g., extracts 10 from 'dyn10'\n", + " hr_label = f\"{hr_value}_Kmin\"\n", + " print(hr_label)\n", + " rep_label = f\"Rep_{name_parts[-1][1]}\" # e.g., 'Rep_1' from '..._r1.txt'\n", + " \n", + " # Read CSV file as Pandas DataFrame\n", + " exp_path = os.path.join(exp_root, file_name)\n", + " exp_df = pd.read_csv(exp_path, header=0, skiprows=[1], \n", + " delimiter=('\\t'), encoding=\"cp858\")\n", + " \n", + " # Adjust temperature to Kelvin\n", + " exp_df[\"ts\"] = exp_df[\"ts\"] + 273.15\n", + " exp_df[\"tr\"] = exp_df[\"tr\"] + 273.15\n", + " \n", + " # Add DataFrame to database\n", + " fsp.pyrolysis.kinetics.add_constant_heating_rate_tga(\n", + " database=PMMA_data_1mg, \n", + " condition=hr_label, \n", + " repetition=rep_label, \n", + " raw_data=exp_df, \n", + " data_type=\"integral\", \n", + " set_value=[hr_value, \"K/min\"])\n", + " \n", + " \n", + "# Adjust column mapping for later functions \n", + "column_mapping = {\n", + " 'time': 't',\n", + " 'temp': 'ts',\n", + " 'signal': 'weight'}\n", + "for hr_label in PMMA_data_1mg[\"experiments\"][\"TGA\"][\"constant_heating_rate\"]:\n", + " # Compute averages and standard deviations per heating rate\n", + " fsp.pyrolysis.kinetics.combine_repetitions(\n", + " database=PMMA_data_1mg, \n", + " condition=hr_label, \n", + " temp_program=\"constant_heating_rate\",\n", + " column_mapping=column_mapping)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, the conversion curves are computed using the function :func:`compute_conversion`. This function operates on the averaged data produced by :func:`combine_repetitions` in the previous step.\n", + "\n", + "After this step, the computed conversion data will be available in the data structure and can be visualized or used in the isoconversional analysis.\n", + "\n", + "To compute conversions for all available heating rates, simply run:\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "fsp.pyrolysis.kinetics.compute_conversion(\n", + " database=PMMA_data_1mg, \n", + " condition=\"all\", \n", + " setup=\"constant_heating_rate\") \n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "Now, the desired conversion levels must be specified to indicate where the apparent activation energy $E_a$ should be evaluated. \n", + "\n", + "Commonly, these levels lie in the range $0.10 < \\alpha < 0.90$ or $0.05 < \\alpha < 0.95$, depending on the quality of the experimental data. In the early and late stages of the reaction, the signal changes more slowly, and even small fluctuations or noise can lead to significant artifacts in the computed activation energy. This is typically visible in the $E_a$ versus conversion plots (see below).\n", + "\n", + "However, these regions should not be discarded blindly. It's important to visually inspect the data to decide whether and to what extent the tails of the conversion curve can be included in the analysis.\n", + "\n", + "The desired conversion levels can be provided as a NumPy array. The function `compute_conversion_levels` performs a linear interpolation of the available data to these levels in preparation for estimating the apparent activation energy.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Define conversion fractions where to evaluate the activation energy\n", + "# Note: commonly, they range between (0.05 < α < 0.95) or (0.1 < α < 0.9)\n", + "# conversion_levels = np.linspace(0.05, 0.95, 37) # Δα = 2.5\n", + "conversion_levels = np.linspace(0.01, 0.99, 99) # Δα = 1.0\n", + "\n", + "\n", + "fsp.pyrolysis.kinetics.compute_conversion_levels(\n", + " database=PMMA_data_1mg, \n", + " desired_levels=conversion_levels,\n", + " setup=\"constant_heating_rate\", \n", + " condition=\"all\")\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "The function `compute_Ea_KAS` performs the final step of the isoconversional method. Internally, it first sorts the available data by heating rate (from lowest to highest) for each conversion level. Then, for each level, it carries out a linear regression according to the KAS method (with the user-defined parameters $B$ and $C$).\n", + "\n", + "The estimated activation energy values $E_a$ are stored in the database, alongside the intermediate fit data and statistics.\n", + "\n", + "This step concludes the KAS-based estimation of the apparent activation energy.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Compute the activation energy using the KAS method\n", + "fsp.pyrolysis.kinetics.compute_Ea_KAS(\n", + " database=PMMA_data_1mg, \n", + " B=1.92, \n", + " C=1.0008)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "Below the result is plotted: the apparent activation energy against the conversion. In gray, the 5 % range at the ends is indicated.\n", + "\n", + "The $E_a$ is computed in J/mol and converted here to kJ/mol, as it is a commen way to use it. \n", + "\n", + "\n", + "\n", + "Below, the apparent activation energy $E_a$ is plotted against the conversion $\\alpha$.\n", + "\n", + "The KAS method estimates $E_a$ in units of J/mol, but it is commonly reported in kJ/mol — so the values are converted accordingly before plotting.\n", + "\n", + "The shaded gray regions indicate the first and last 5% of the conversion range. These are typically excluded from analysis due to higher sensitivity to noise and lower data reliability. In the first 5%, artifacts are visible due to increased noise.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots()\n", + "\n", + "\n", + "# Get the Ea and convert to kJ/mol.\n", + "Ea_results_KAS = PMMA_data_1mg[\"experiments\"][\"TGA\"][\"Ea_results_KAS\"]\n", + "Ea = Ea_results_KAS[\"Ea\"]/1000\n", + "Ea_avg = np.average(Ea)\n", + "\n", + "# Plot the Ea against conversion.\n", + "conv = Ea_results_KAS[\"Conversion\"]\n", + "plt.scatter(conv,\n", + " Ea,\n", + " marker=\".\", s=42,\n", + " facecolors=\"none\",\n", + " edgecolors=\"tab:blue\",\n", + " label=f\"KAS, ΔT = {ΔT} K\")\n", + "\n", + "\n", + "# Shaded areas to indicate first/last 5 % (typically excluded)\n", + "x_min = -0.025\n", + "x_max = 1.025\n", + "ax.axvspan(x_min, 0.05, color='gray', alpha=0.3, label=\"5%\")\n", + "ax.axvspan(0.95, x_max, color='gray', alpha=0.3)\n", + "\n", + "\n", + "# Plot meta data.\n", + "plt.title(f\"Activation Energy (KAS), E$_a$={Ea_avg:.2f} kJ/mol (avg.)\")\n", + "plt.xlabel(\"Conversion ($\\\\alpha$) / -\")\n", + "plt.ylabel(\"Activation Energy (E$_a$) / kJ/mol\")\n", + "\n", + "plt.xlim(left=x_min, right=x_max)\n", + "plt.ylim(bottom=68, top=257)\n", + "\n", + "plt.tight_layout()\n", + "plt.legend(loc=\"lower center\")\n", + "plt.grid()\n", + "\n", + "\n", + "# # Save image.\n", + "# plot_label = f\"Ea_Estimate_KAS.png\"\n", + "# plot_path = os.path.join(plot_label)\n", + "# plt.savefig(plot_path, dpi=320, bbox_inches='tight', facecolor='w')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "FireSciPy also evaluates the quality of each linear fit used in the KAS computation. Specifically, the **root mean square error (RMSE)** and the **coefficient of determination (R²)** are calculated.\n", + "\n", + "- An **RMSE** of 0 indicates a perfect fit.\n", + "- An **R²** of 1 means that the fit perfectly explains the variance in the data.\n", + "\n", + "In the plot below, fluctuations are more pronounced at low levels of conversion. This further supports the common practice of excluding the edges of the conversion range (typically the first and last 5%) in isoconversional analyses.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Pre-select data sets for convenience\n", + "Ea_KAS = PMMA_data_1mg[\"experiments\"]['TGA'][\"Ea_results_KAS\"]\n", + "\n", + "# Optional: Reduce number of data points by using every n-th\n", + "nth = 1\n", + "\n", + "markersize = 42\n", + "\n", + "fig, ax = plt.subplots()\n", + "\n", + "# Plot R² statistics\n", + "x_data = np.asarray(Ea_KAS[\"Conversion\"])[::nth]\n", + "y_data = np.asarray(Ea_KAS[\"R_squared\"])[::nth]\n", + "plt.scatter(x_data,\n", + " y_data,\n", + " marker='.', s=markersize,\n", + " facecolors='none',\n", + " edgecolors='tab:blue',\n", + " label=f\"R² avrg.: {np.average(y_data):.6f}\")\n", + "\n", + "\n", + "# Plot RMSE statistics\n", + "y_data = np.asarray(Ea_KAS[\"RMSE\"])[::nth]\n", + "plt.scatter(x_data,\n", + " y_data,\n", + " marker='.', s=markersize,\n", + " facecolors='none',\n", + " edgecolors='tab:orange',\n", + " label=f\"RMSE avrg.: {np.average(y_data):.6f}\")\n", + "\n", + "\n", + "# Shaded areas to indicate first/last 5 % (typically excluded)\n", + "x_min = -0.025\n", + "x_max = 1.025\n", + "ax.axvspan(x_min, 0.05, color='gray', alpha=0.3, label=\"5%\")\n", + "ax.axvspan(0.95, x_max, color='gray', alpha=0.3)\n", + "\n", + "\n", + "# Plot meta data.\n", + "plt.title(\"Statistics of the KAS Fits\")\n", + "plt.xlabel(\"Conversion, $\\\\alpha$ / -\")\n", + "plt.ylabel('Arbitrary Units / -')\n", + "\n", + "plt.xlim(left=-0.025, right=1.025)\n", + "plt.ylim(bottom=-0.05, top=1.05)\n", + "\n", + "plt.tight_layout()\n", + "plt.legend()\n", + "plt.grid()\n", + "\n", + "# # Save image.\n", + "# plot_label = \"Ea_Estimate_KAS_Statistics.png\"\n", + "# plot_path = os.path.join(plot_label)\n", + "# plt.savefig(plot_path, dpi=320, bbox_inches='tight', facecolor='w')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "Data produced during intermediate steps can also be accessed. For example, the averaged TG curves from the combined data sets can be used to plot normalised mass loss rates (MLR). This is shown below for heating rates of 30 K/min and 60 K/min.\n", + "\n", + "The curves are aligned such that the final mass approaches zero, which is common when no significant residue remains. This zeroing helps compare curves consistently even when initial/final absolute values vary slightly. The MLR is then normalised by the initial mass to allow comparison between different heating rates.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "hr_labels = [\"30_Kmin\", \"60_Kmin\"]\n", + "\n", + "for hr_label in hr_labels:\n", + " # Access combined data\n", + " PMMA_exp = PMMA_data_1mg[\"experiments\"][\"TGA\"][\"constant_heating_rate\"][hr_label][\"combined\"]\n", + "\n", + " # Compute normalised mass loss rate\n", + " time = PMMA_exp[\"Time\"]\n", + " temp_avg = PMMA_exp[\"Temperature_Avg\"]\n", + " mass_avg = PMMA_exp[\"Mass_Avg\"]\n", + "\n", + " mlr = -np.gradient(mass_avg, time, edge_order=2)\n", + " mlr_smooth = uniform_filter1d(mlr, size=9)\n", + " \n", + " # Plot the mass loss rate\n", + " plt.plot(temp_avg, \n", + " mlr_smooth,\n", + " label=f\"{hr_label.split('_')[0]} K/min (1mg)\")\n", + "\n", + "\n", + "# Plot meta data.\n", + "plt.xlabel(\"Sample Temperature / K\")\n", + "plt.ylabel(\"Normalised Mass Loss Rate / 1/s\")\n", + "\n", + "plt.xlim(left=380,right=770)\n", + "plt.ylim(bottom=-0.0005,top=0.0165)\n", + "\n", + "plt.legend()\n", + "plt.grid()\n", + "\n", + "\n", + "# Save image.\n", + "plot_label = f\"NormalisedMLR.png\"\n", + "plot_path = os.path.join(plot_label)\n", + "plt.savefig(plot_path, dpi=320, bbox_inches='tight', facecolor='w')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For more details, see the [FireSciPy documentation](https://firedynamics.github.io/FireSciPy/).\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sensitivity of $E_a$ Estimation to the Number and Range of Heating Rates\n", + "-----------------------------------------------------------------------------\n", + "\n", + "When estimating the apparent activation energy $E_a$, one common question is: \n", + "**\"Why are so many heating rates necessary? Can’t we just use three?\"**\n", + "\n", + "The answer to this question boils down to the fact that the isoconversional methods are fundamentally built on linear regression. For each chosen conversion level, the KAS method fits a straight line to $ln⁡(\\beta/T^B)$ vs. $1/T$. The slope of this line is used to determine the apparent activation energy $E_a$.\n", + "\n", + "According to the [ICTAC Kinetics Committee recommendations (https://doi.org/10.1016/j.tca.2014.05.036)](https://doi.org/10.1016/j.tca.2014.05.036), at the very least three, better are **five or more**, temperature programs should be used. These should span as wide a range as possible. Consider linear heating rates. ICTAC suggests to use a spread of a **factor of 10 or more** between the lowest and highest rate (e.g. 5 K/min to 50 K/min) to ensure robust estimation of $E_a$. This ensures:\n", + "- A wide spread in $1/T$ values, which increases statistical leverage in the fit.\n", + "- Enough data points to reduce sensitivity to noise and improve the robustness of the slope estimate.\n", + "\n", + "If the heating rates are too close together or too few in number, the $1/T$ values cluster. This reduces the fit’s ability to capture the underlying trend, making the slope — and thus $E_a$ — more sensitive to experimental noise and measurement uncertainty.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Initialise data structure for the kinetics assessment\n", + "PMMA_data_1mg_low = fsp.pyrolysis.kinetics.initialize_investigation_skeleton(\n", + " material=f\"PMMA\", \n", + " investigator=\"John Doe, Miskatonic University\", \n", + " instrument=\"TGA/DSC 3+, Mettler Toledo\", \n", + " date=\"Stardate: 42.69\", \n", + " notes=\"Constant heating rates, sample mass 1mg\",\n", + " signal={\"name\": \"Mass\", \"unit\": \"mg\"})\n", + "\n", + "\n", + "# heating rates: 5, 10 K/min\n", + "file_names = [\n", + " \"tga_dynamic_n2_dyn5_powder_1mg_r1.txt\",\n", + " \"tga_dynamic_n2_dyn5_powder_1mg_r2.txt\",\n", + " \"tga_dynamic_n2_dyn5_powder_1mg_r3.txt\",\n", + " \"tga_dynamic_n2_dyn10_powder_1mg_r1.txt\",\n", + " \"tga_dynamic_n2_dyn10_powder_1mg_r2.txt\",\n", + " \"tga_dynamic_n2_dyn10_powder_1mg_r3.txt\"\n", + "]\n", + "\n", + "\n", + "for file_name in file_names:\n", + " if \"powder_1mg_\" in file_name:\n", + " # print(file_name)\n", + " \n", + " # Parse metadata from file name: heating rate and repetition\n", + " name_parts = file_name.split(\"_\")\n", + " hr_value = int(name_parts[3][3:]) # e.g., extracts 10 from 'dyn10'\n", + " hr_label = f\"{hr_value}_Kmin\"\n", + " print(hr_label)\n", + " rep_label = f\"Rep_{name_parts[-1][1]}\" # e.g., 'Rep_1' from '..._r1.txt'\n", + " \n", + " # Read CSV file as Pandas DataFrame\n", + " exp_path = os.path.join(exp_root, file_name)\n", + " exp_df = pd.read_csv(exp_path, header=0, skiprows=[1], \n", + " delimiter=('\\t'), encoding=\"cp858\")\n", + " \n", + " # Adjust temperature to Kelvin\n", + " exp_df[\"ts\"] = exp_df[\"ts\"] + 273.15\n", + " exp_df[\"tr\"] = exp_df[\"tr\"] + 273.15\n", + " \n", + " # Add DataFrame to database\n", + " fsp.pyrolysis.kinetics.add_constant_heating_rate_tga(\n", + " database=PMMA_data_1mg_low, \n", + " condition=hr_label, \n", + " repetition=rep_label, \n", + " raw_data=exp_df, \n", + " data_type=\"integral\", \n", + " set_value=[hr_value, \"K/min\"])\n", + " \n", + " \n", + "# Adjust column mapping for later functions \n", + "column_mapping = {\n", + " 'time': 't',\n", + " 'temp': 'ts',\n", + " 'signal': 'weight'}\n", + "for hr_label in PMMA_data_1mg_low[\"experiments\"][\"TGA\"][\"constant_heating_rate\"]:\n", + " # Compute averages and standard deviations per heating rate\n", + " fsp.pyrolysis.kinetics.combine_repetitions(\n", + " database=PMMA_data_1mg_low, \n", + " condition=hr_label, \n", + " temp_program=\"constant_heating_rate\",\n", + " column_mapping=column_mapping)\n", + " \n", + "\n", + "fsp.pyrolysis.kinetics.compute_conversion(\n", + " database=PMMA_data_1mg_low, \n", + " condition=\"all\", \n", + " setup=\"constant_heating_rate\") \n", + " \n", + "\n", + "# Define conversion fractions where to evaluate the activation energy\n", + "# Note: commonly, they range between (0.05 < α < 0.95) or (0.1 < α < 0.9)\n", + "# conversion_levels = np.linspace(0.05, 0.95, 37) # Δα = 2.5\n", + "conversion_levels = np.linspace(0.01, 0.99, 99) # Δα = 1.0\n", + "\n", + "fsp.pyrolysis.kinetics.compute_conversion_levels(\n", + " database=PMMA_data_1mg_low, \n", + " desired_levels=conversion_levels,\n", + " setup=\"constant_heating_rate\", \n", + " condition=\"all\")\n", + "\n", + "\n", + "# Compute the activation energy using the KAS method\n", + "fsp.pyrolysis.kinetics.compute_Ea_KAS(\n", + " database=PMMA_data_1mg_low, \n", + " B=1.92, \n", + " C=1.0008)\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Pre-select data sets for convenience\n", + "Ea_KAS = PMMA_data_1mg[\"experiments\"]['TGA'][\"Ea_results_KAS\"]\n", + "Ea_KAS_example = PMMA_data_1mg_low[\"experiments\"]['TGA'][\"Ea_results_KAS\"]\n", + "hr_labels = [\"5_Kmin\", \"10_Kmin\", \"20_Kmin\", \"30_Kmin\", \"60_Kmin\"]\n", + "\n", + "# Define settings for the plots of the fits\n", + "marker_size = 42\n", + "fit_line = \":\"\n", + "fit_alpha = 0.8\n", + "fit_color = \"black\"\n", + "\n", + "# Choose conversion levels for the plot\n", + "conversion_levels = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]\n", + "\n", + "# Initialise data collection\n", + "conv_temps = np.zeros((len(conversion_levels), len(hr_labels)))\n", + "conv_levels = np.zeros((len(conversion_levels), len(hr_labels)))\n", + "\n", + "for conv_idx, level in enumerate(conversion_levels):\n", + " level_idx = np.abs(Ea_KAS[\"Conversion\"] - level).argmin()\n", + " # Get the first three data points to match previous plot\n", + " conv_temps[conv_idx,:] = np.asarray(Ea_KAS.loc[:,\"x1\":\"x5\"].iloc[level_idx])\n", + " conv_levels[conv_idx,:] = np.asarray(Ea_KAS.loc[:,\"y1\":\"y5\"].iloc[level_idx])\n", + " # Indicate fits\n", + " m_fit = Ea_KAS[\"m_fit\"].iloc[level_idx]\n", + " b_fit = Ea_KAS[\"b_fit\"].iloc[level_idx]\n", + " x_fit = [conv_temps[conv_idx,:][0], conv_temps[conv_idx,:][-1]]\n", + " y_fit = [fsp.utils.linear_model(x_fit[0], m_fit, b_fit),\n", + " fsp.utils.linear_model(x_fit[-1], m_fit, b_fit)]\n", + " \n", + " if conv_idx == 0:\n", + " plot_label = \"Fit (full)\"\n", + " else:\n", + " plot_label = \"_none_\"\n", + " plt.plot(x_fit, y_fit,\n", + " linestyle=fit_line,\n", + " alpha=fit_alpha,\n", + " color=fit_color,\n", + " label=plot_label)\n", + " \n", + " # Indicate fits, extreme example\n", + " m_fit = Ea_KAS_example[\"m_fit\"].iloc[level_idx]\n", + " b_fit = Ea_KAS_example[\"b_fit\"].iloc[level_idx]\n", + " \n", + " y_fit = [fsp.utils.linear_model(x_fit[0], m_fit, b_fit),\n", + " fsp.utils.linear_model(x_fit[-1], m_fit, b_fit)]\n", + " \n", + " if conv_idx == 0:\n", + " plot_label = \"Fit (5, 10)\"\n", + " else:\n", + " plot_label = \"_none_\"\n", + " plt.plot(x_fit, y_fit,\n", + " linestyle=fit_line,\n", + " alpha=fit_alpha,\n", + " color=\"tab:red\",\n", + " label=plot_label)\n", + "\n", + "# Plot data points by heating rate\n", + "for idx in range(len(conv_temps.T)):\n", + " # Get colour for data series, i.e. heating rate\n", + " plot_colour = plt_colors[idx]\n", + " \n", + " # Plot data points\n", + " plt.scatter(\n", + " conv_temps.T[idx],\n", + " conv_levels.T[idx],\n", + " marker='o', s=marker_size,\n", + " facecolors='none',\n", + " edgecolors=plot_colour,\n", + " label=f\"{hr_labels[idx].split('_')[0]} K/min\")\n", + "\n", + "\n", + "# Plot meta data.\n", + "plt.title(\"Wide Heating Rate Ranges Improve Slope Stability in KAS\")\n", + "plt.xlabel(\"1/T\")\n", + "plt.ylabel(\"ln($\\\\beta$/T$^{1.92}$)\")\n", + "\n", + "plt.xlim(left=0.00141, right=0.00186)\n", + "plt.ylim(bottom=-11.1, top=-7.9)\n", + "\n", + "plt.tight_layout()\n", + "plt.legend()\n", + "plt.grid()\n", + "\n", + "# # Save image.\n", + "# plot_label = f\"Ea_Estimate_KAS_Fit_PMMA1mg_example.png\"\n", + "# plot_path = os.path.join(plot_label)\n", + "# plt.savefig(plot_path, dpi=320, bbox_inches='tight', facecolor='w')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "The figure above illustrates this effect for a few conversion levels. \n", + "- Black dotted lines: fits using the full dataset (five heating rates from 5 to 60 K/min).\n", + "- Red dotted lines: fits using only the two lowest heating rates (5 and 10 K/min).\n", + "\n", + "With the narrower range, the points lie closer together along the $1/T$ axis, and the slope differs noticeably from the full-data case. This illustrates why **both range and redundancy** are critical when designing thermal analysis experiments for kinetic modeling.\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.9" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}