Adiabatic flame temperature - methane and oxygen¶

To calculate the adiabatic flame temperature of the CH4 + O2 reaction, the total enthalpy balance of this reaction is set to 0 Joule.

The necessary thermodynamic data is stored in the file methane_combustion.dat. We are solely interested in the gaseous phase constituents.

Loading datafile and preparing the calculation¶

A ThermochemicalSystem is typically instantiated by loading a .dat or .cst file. It has functions that allow to get information about system components, phases etc.

However, contrary to the Fortran/C/other language libraries, in ChemApp for Python, there is little need to find e.g. indexes of phases or constituents, since they can be directly adressed by names.

To set the units is a good precaution. The inputs to the settings are 'guarded' by using enumerations that only allow certain limits of inputs. Therefore, we need to import the respective enums for pressure, volume, temperature, energy and amounts.

In [ ]:

from chemapp.friendly import Units

# Make sure the units are set as expected
from chemapp.core import (
    PressureUnit,
    VolumeUnit,
    TemperatureUnit,
    EnergyUnit,
    AmountUnit,
)

Units.set(
    T=TemperatureUnit.C,
    P=PressureUnit.atm,
    V=VolumeUnit.dm3,
    A=AmountUnit.mol,
    E=EnergyUnit.J,
)

Setup of the calculation¶

In order to setup, we use the EqulibriumCalculator from chemapp.friendly. The EquilibriumCalculator is

In [ ]:

from chemapp.friendly import EquilibriumCalculation, StreamCalculation
from chemapp.friendly import ConditionVariable
import chemapp.basic as ChemAppBasic

The friendly layer allows to e.g. set the incoming amount of a phase constituent, using the routine set_IA_pc(Phase, Constituent, Amount) in an intuitive way. Please note that we do not need to acquire any indexes to address phases or constituents, even though we could.

In [ ]:

EquilibriumCalculation.set_IA_pc("Gas", "CH4", 0.3)
EquilibriumCalculation.set_IA_pc("Gas", "O2", 0.7)
# EquilibriumCalculation.set_eq_P(1.0)

Out[ ]:

Intuitiveness of function names can get very complicated in case of more complicated calculations. Therefore, it is at all times possible to 'fall back' on the basic layer, which gives access to the full ChemApp API. Mixing the friendly and basic layer is working seamlessly.

To set up the appropriate enthalpy condition, we use the basic layer API function of tqsetc, which is also used behind the curtains of the friendly layer. Here, indexes have to be used for components, phases and constituents.

The ConditionVariable enumeration provides access to all settable condition types in ChemApp.

In [ ]:

# condition for the adiabatic point H(reaction) == 0
# -1, -1 for phase and constituents == System wide
ChemAppBasic.tqsetc(ConditionVariable.H, -1, -1, 0.0)

Out[ ]:

These conditions are all that are needed to perform a target calculation to find the temperature for which (if possible!) all conditions are satisfied. Again, we could use both the EquilibriumCalculation or just as well the tqce routine exposed by the basic layer.

To give a reasonable starting value for the calculation potentially increases calculation performance and convergence significantly in more complicated calculations.

In [ ]:

# perform the calculation with a 'suggestion' of 1000C as the solution
StreamCalculation.calculate_eq_T(1000.0)

flame_temperature = StreamCalculation.get_eq_T()
print(f"The resulting adiabatic flame temperature is {flame_temperature:.2f}°C")

The resulting adiabatic flame temperature is 2822.71°C

Determining the maximum temperature¶

In the example above, we set the system composition in amounts of CH4 and O2, and then solved the thermodynamic equilibrium with a search for the temperature at which the system enthalpy is equal to zero.

Now, we want to look at these temperatures over a wider range of compositions. We will use the power of some 3rd party modules to make quick visualisations possible.

We use both numpy and matplotlib to access easy math functions and operations as well as easy and powerful plotting options.

First, we generate a range of points at which we want to calculate the temperature. Used in this way, it is similar to the input field for A in the Equilib Module, but only scratches the surface of the power of numpy.

Note that we are not using 0 and 1, but just values very close. This is because neither pure CH4 nor pure O2 has an adiabatic flame temperature.

In [ ]:

A_list = np.linspace(
    0.1e-6, # lower end
    1.0,  # upper end
    num=21,  # number of points
)
# to save the results, an empty list
temperatures = []

A simple iteration over the aforementioned range and appending the resulting temperature to a list, using the same set of conditions setting routines as well as calculation routine. Notably, only the conditions that need to be updated are updated in the loop. The condition on enthalpy is maintained throughout all the calculations.

In [ ]:

for A in A_list:
    EquilibriumCalculation.set_IA_pc("Gas", "CH4", A)
    EquilibriumCalculation.set_IA_pc("Gas", "O2", 1.0 - A)
    StreamCalculation.calculate_eq_T(1000.0)

    temperatures.append(StreamCalculation.get_eq_T())

This gives us an unimpressive list of numbers.

In [ ]:

print(f"Calculated temperatures: \n{temperatures}")

Calculated temperatures: 
[25.002991221208163, 1262.8414195712353, 2113.406486948618, 2479.6343157889733, 2655.544818355322, 2759.459824247294, 2822.7052935764755, 2852.514795785148, 2842.4718952655603, 2770.3181885350537, 2602.461212392842, 2320.2580227766384, 1906.057660319348, 1394.3597286843883, 1153.5457621869682, 1046.4701289174955, 977.1447052685186, 924.0315793709973, 878.0360065866811, 834.480279835854, 790.1096260213686]

Which we can easily turn into a plot which is supplied by the matplotlib package, again only scratching the surface of the abilities and possibilities to create beautiful and efficient plots.

Determining the composition of the maximum temperature¶

Leveraging the many modules available with Python, we will employ the scipy.optimize module, namely the minimize_scalar function, which is a blackbox optimizer that can search for minimal points in arbitrary linear parameter spaces.

To be able to optimize the temperature as function of A, which expresses the ratio between CH4 and O2, we simply write a function that takes A as a parameter and returns the negative value of the temperature, in order to minimize the outcome of the function.

In [ ]:

def negative_adiabatic_temperature(A):
    EquilibriumCalculation.set_IA_pc("Gas", "CH4", A)
    EquilibriumCalculation.set_IA_pc("Gas", "O2", 1.0 - A)
    StreamCalculation.calculate_eq_T(1000.0)

    return -StreamCalculation.get_eq_T()

After having the function created, it is a simple one-liner to call the minimization of that function. There are a few options for the method parameter of that routine, here we use bounded so that we can set the upper and lower limits for A to ~ 0 - 1.

In [ ]:

result = minimize_scalar(
    negative_adiabatic_temperature,
    bounds=[1.0e-4, 1.0 - 1e-4],
    method="bounded",
)

# The result object contains every information we want
t_max = -1.0 * result.fun
A_max = result.x
print(
    f"The maximum temperature of {t_max:.2f}°C "
    f"is reached for {A_max:.3f}mol CH4 + {1-A_max:.3f}mol O2."
)

# plotting it quickly
plt.ylabel("Temperature [°C]")
plt.xlabel("A·CH4 + (1-A)·O2")
plt.plot(A_list, temperatures)
plt.axvline(x=A_max, ymin=0, ymax=1, color="orange")
plt.show()

The maximum temperature of 2854.35°C is reached for 0.365mol CH4 + 0.635mol O2.