Setting up a multi-thread laminar flame speed parameter study#

One of the prevailing use case of the freely propagating premixed flame model is to build a flame speed table to be imported by another combustion simulation tools. PyChemkin provides the flexibility to customize the data structure of the flame speed table depending on the simulation goals and the tool. Furthermore, over the years, the chemkin flame speed calculator has derived a set of default solver settings that would greatly improve the convergence performance, especially for those widely adopted hydrocarbon fuel combustion mechanisms. The required input parameters the flame speed calculator are reduced to the composition of the fuel-oxidizer mixture, the initial/inlet pressure and temperature, and the calculation domain.

This tutorial shows the steps of setting up a flame speed parameter study for CH₄-air mixtures at the 5 atmosphere pressure. The predicted flame speed values are compared against the experimental data as a function of the mixture equivalence ratio. The parameter study is performed in the multi-thread mode by using the threading package.

Since the transport processes are critical for flame calculations, the transport data must be included in the mechanism data and preprocessed.

Import PyChemkin packages and start the logger#

import os
from pathlib import Path
import threading
import time

import ansys.chemkin.core as ck  # Chemkin
from ansys.chemkin.core.inlet import Stream  # external gaseous inlet
from ansys.chemkin.core.logger import logger

# Chemkin 1-D premixed freely propagating flame model (steady-state)
from ansys.chemkin.core.premixedflames.premixedflame import (
    FreelyPropagating as FlameSpeed,
)
import matplotlib.pyplot as plt  # plotting
import numpy as np  # number crunching

# check working directory
current_dir = str(Path.cwd())
logger.debug("working directory: " + current_dir)
# set verbose mode
ck.set_verbose(True)
# set interactive mode for plotting the results
# interactive = True: display plot
# interactive = False: save plot as a PNG file
global interactive
interactive = True

Create a flame speed calculator class#

Create a local class that wraps around the actual FlameSpeed class to make the setup of the multi-thread flame speed calculation parameter study more convenient.

class FlameSpeedCalculator:
    """Laminar flame speed calculator with fixed set up parameters."""

    def __init__(self, fresh_mixture: Stream, index: int):
        """Create a Laminar flame speed calculator."""
        """
        Create a Laminar flame speed calculator that instantiates a FlameSpeed object
        with the given fresh (unburnt) mixture condition.

        Parameters
        ----------
            fresh_mixture: Mixture object
                the initial/fresh/unburnt condition
            index: integer
                run index of this flame speed calculator
        """
        # instantiate a flame speed object
        # set up the run and working directory name
        name = "Flame_Speed_" + str(index)
        # instantiate the FlameSpeed object for this run
        self.fs_calculator = FlameSpeed(fresh_mixture, label=name)
        # set the required premixed flame model parameters
        #
        # set the maximum total number of grid points allowed
        # in the calculation (optional)
        # self.fs_calculator.set_max_grid_points(150)
        # define the calculation domain [cm]
        self.fs_calculator.end_position = 1.0
        # set the root directory
        self.root_dir = str(Path.cwd())
        # set the working directory
        self.work_dir = str(Path(self.root_dir) / name)
        # run status
        self.runstatus = -100
        # calculated laminar flame speed [cm/sec]
        self.flame_speed = 0.0

    def run(self):
        """Run the flame speed calculation in a separate working directory."""
        # create or clean up the working directory for this run
        this_work_folder = Path(self.work_dir)
        if this_work_folder.is_dir():
            # directory exists
            for f in this_work_folder.iterdir():
                if f.is_file():
                    # remove any existing file
                    try:
                        f.unlink()
                    except OSError as e:
                        print(f"Error removing {str(f)}: {e}")
        else:
            # create a new directory
            this_work_folder.mkdir()
        # change to the working directory for this run
        os.chdir(self.work_dir)
        # run the flame speed calculation
        self.runstatus = self.fs_calculator.run()
        # extract the laminar flame speed from the solution
        if self.runstatus == 0:
            # postprocess the solutions
            self.fs_calculator.process_solution()
            # get the flame speed value [cm/sec]
            # because the memory is shared, it must be done as soon as
            # the run is finished
            self.flame_speed = self.fs_calculator.get_flame_speed()
        # go back the root directory
        os.chdir(self.root_dir)

    def get_flame_speed(self) -> float:
        """Get the predicted laminar flame speed."""
        """
        Get the predicted laminar flame speed.

        Returns
        -------
            flame_speed: double
                predicted laminar flame speed [cm/sec]
        """
        return self.flame_speed

Set up the flame speed parameter study for multi-threading#

Create a list of FlameSpeedCalculator objects with different initial methane-air equivalence ratios from 0.6 to 1.6. Each object represents one parameter study case and will be assigned to its own thread when the parameter study is executed.

def prepare_multi_thread_runs() -> dict[float, FlameSpeedCalculator]:
    """Set up the flame speed parameter study."""
    """
    Set up the parameter study runs for multi-threading.

    Returns
    -------
        flame_speed_runs: list of FlameSpeedCalculator objects
            flame speed calculation cases
    """
    ##########################################
    # Create an instance of the Chemistry Set
    # ========================================
    # The mechanism loaded is the GRI 3.0 mechanism for methane combustion.
    # The mechanism and its associated data files come with the standard Ansys Chemkin
    # installation under the subdirectory *"/reaction/data"*.
    #

    # set mechanism directory (the default Chemkin mechanism data directory)
    data_dir = Path(ck.ansys_dir) / "reaction" / "data"
    mechanism_dir = data_dir
    # including the full file path is recommended
    chemfile = str(mechanism_dir / "grimech30_chem.inp")
    thermfile = str(mechanism_dir / "grimech30_thermo.dat")
    tranfile = str(mechanism_dir / "grimech30_transport.dat")
    # create a chemistry set based on GRI 3.0
    MyGasMech = ck.Chemistry(
        chem=chemfile, therm=thermfile, tran=tranfile, label="GRI 3.0"
    )

    ##############################
    # Preprocess the Chemistry Set
    # ============================

    # preprocess the mechanism files
    ierror = MyGasMech.preprocess()
    if ierror != 0:
        print("Error: failed to preprocess the mechanism!")
        print(f"       error code = {ierror}")
        exit()

    ########################################################################
    # Set up the CH\ :sub:`4`\ -air mixture for the flame speed calculation
    # ======================================================================
    # Instantiate a stream named ``premixed`` for the inlet gas mixture.
    # This stream  is a mixture with the addition of the
    # inlet flow rate. You can specify the inlet gas properties the same way you
    # set up a ``Mixture``. Here the ``x_by_equivalence_ratio`` method is used.
    # You create the ``fuel`` and the ``air`` mixtures first. Then define the
    # *complete combustion product species* and provide the *additives* composition
    # if applicable. And finally, during the parameter iteration runs, you can
    # simply set different values to ``equivalenceratio`` to create different
    # methane-air mixtures.
    #

    # create the fuel mixture
    fuel = ck.Mixture(MyGasMech)
    # set fuel composition: methane
    fuel.x = [("CH4", 1.0)]
    # setting pressure and temperature condition for the flame speed calculations
    fuel.pressure = 5.0 * ck.P_ATM
    fuel.temperature = 300.0  # inlet temperature

    # create the oxidizer mixture: air
    air = ck.Mixture(MyGasMech)
    air.x = ck.Air.x()
    # setting pressure and temperature is not required in this case
    air.pressure = fuel.pressure
    air.temperature = fuel.temperature

    # create the fuel-air Stream for the premixed flame speed calculation
    premixed = Stream(MyGasMech, label="premixed")
    # products from the complete combustion of the fuel mixture and air
    products = ["CO2", "H2O", "N2"]
    # species mole fractions of added/inert mixture.
    # can also create an additives mixture here
    add_frac = np.zeros(MyGasMech.kk, dtype=np.double)  # no additives: all zeros

    # setting pressure and temperature is not required in this case
    premixed.pressure = fuel.pressure
    premixed.temperature = fuel.temperature

    # set estimated value of the flame mass flux [g/cm2-sec]
    premixed.mass_flowrate = 0.4

    # Set up the flame speed parameter study for multi-threading
    # equivalence ratio for the first case
    phi = 0.6
    # total number of parameter cases
    points = 21
    # equivalence ratio increment
    delta_phi = 0.05
    # set up flame speed calculation runs
    flame_speed_runs: dict[float, FlameSpeedCalculator] = {}
    for i in range(points):
        # create mixture by using the equivalence ratio
        ierror = premixed.x_by_equivalence_ratio(
            MyGasMech, fuel.x, air.x, add_frac, products, equivalenceratio=phi
        )
        # check fuel-oxidizer mixture creation status
        if ierror != 0:
            print(
                "Error: failed to create the methane-air mixture "
                + "for equivalence ratio = "
                + str(phi)
            )
            exit()
        # create a flame speed calculation instance
        flame_speed_runs[phi] = FlameSpeedCalculator(premixed, index=i)
        # update parameter
        phi += delta_phi

    # return the job setup parameters
    return flame_speed_runs

Set up and start the multi-thread runs#

Use the Thread() method to assign the flame speed runs. In this project, each flame speed run/case has its own thread. The target parameter of Thread() method should be the run() method of the FlameSpeedCalculator. Use the start() method to initiate the threads, and the join() method to sync the threads after they are done.

flame_speed_cases = prepare_multi_thread_runs()
# set the start wall time
start_time = time.time()
threads = []
for phi, fsc in flame_speed_cases.items():
    t = threading.Thread(target=fsc.run())
    threads.append(t)
    # start each thread
    t.start()
# wait for all threads to finish
for t in threads:
    t.join()

# compute the total runtime
runtime = time.time() - start_time
print()
print(f"total simulation duration: {runtime} [sec]")
print()

Get the parameter study results#

Use the get_flame_speed() method to get the calculated laminar flame speed values from each run. Then set up the experimental flame speed data for comparison.

points = len(flame_speed_cases.keys())
equival = np.zeros(points, dtype=np.double)
flamespeed = np.zeros_like(equival, dtype=np.double)
for i, case in enumerate(flame_speed_cases.items()):
    # equivalence ratio
    equival[i] = case[0]
    # flame speed calculator case
    fsc = case[1]
    flamespeed[i] = fsc.get_flame_speed()

# experimental data by Kochar
# equivalence ratios
data_equiv = [
    0.7005,
    0.8007,
    0.9009,
    1.001,
    1.1032,
    1.2014,
    1.3014,
]
# methane flame speeds at 5 atm
data_speed = [
    6.906,
    12.0094,
    15.9072,
    19.2376,
    19.6601,
    15.8274,
    10.2925,
]

Plot the premixed flame solution profiles#

Plot the predicted flame speeds against the experimental data.

plt.plot(data_equiv, data_speed, label="data", linestyle="", marker="^", color="blue")
plt.plot(equival, flamespeed, label="GRI 3.0", linestyle="-", color="blue")
plt.legend()
plt.ylabel("Flame Speed [cm/sec]")
plt.xlabel("Equivalence Ratio")

# clean up
ck.done()

# plot results
if interactive:
    plt.show()
else:
    plt.savefig("plot_flame_speed_threading.png", bbox_inches="tight")

Download Jupyter notebook: flame_speed_calculation_threading.ipynb

Download Python source code: flame_speed_calculation_threading.py

Download zipped: flame_speed_calculation_threading.zip

Gallery generated by Sphinx-Gallery