—
Hot spots example: PFR with exponential temperature decay and nucleation analysis.
from pathlib import Path
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import yaml
from bloc import carbon_utils_nobili, ctutils
from bloc.reactors import solve_fixed_time_temperature_profile_pfr, switch_mechanism
from bloc.utils import get_mechanism_path
SCRIPT_DIR = Path(__file__).resolve().parent
RESULTS_BASE = SCRIPT_DIR / "results"
SKIP_SIMULATION = False
def _get_results_dir(inp):
"""Return results subfolder for this mechanism: results/{mech}/."""
mech_string = inp["mech"].replace(".yaml", "")
return RESULTS_BASE / mech_string
def _get_time_ms_str(inp):
"""Return time_exponential in ms as string for filenames, e.g. '10' for 0.01 s."""
time_ms = int(round(inp["time_exponential"] * 1000))
return f"{time_ms}ms"
def _load_inputs(yaml_path=None):
"""Load inputs from hot_spots_inputs.yaml and resolve mechanism path."""
if yaml_path is None:
yaml_path = SCRIPT_DIR / "hot_spots_inputs.yaml"
with open(yaml_path, "r") as f:
inp = yaml.safe_load(f)
return inp
def run(inp):
"""Run the hot spots simulation."""
# Unpack Temperature profile
temp_profile_exponential = inp["temp_profile_exponential"]
initial_temperature = inp["initial_temperature"] # K, intialtemperature
final_temperature = inp["final_temperature"] # K, final temperature
time_exponential = inp["time_exponential"] # s
dt = inp["dt"]
n_time_exponential = inp[
"n_time_exponential"
] # number of exponential times for the temperature profile
# Exponential temperature profile
if temp_profile_exponential:
print(
f"Using exponential temperature profile with time_exponential = {time_exponential:.2e} s"
)
def temperature_profile_exponential(
t, initial_temperature, final_temperature, time_exponential
):
return (-final_temperature + initial_temperature) * np.exp(
-t / time_exponential
) + final_temperature
temperature_profile_time = np.linspace(
0, n_time_exponential * time_exponential, 100
)
temperature_profile = temperature_profile_exponential(
temperature_profile_time,
initial_temperature,
final_temperature,
time_exponential,
)
temperature_profile = np.vstack((temperature_profile_time, temperature_profile))
# Equilibrate at initial temperature with Fincke
print(
f"Equilibrating at initial temperature {initial_temperature} K with Caltech mechanism"
)
gas = ct.Solution(get_mechanism_path("Caltech.yaml"))
gas.TPX = initial_temperature, inp["P_bar"], inp["inlet_composition"]
gas.equilibrate("TP")
# Report major 5 species with mole fraction
top5_idx = np.argsort(gas.X)[-5:][::-1]
print(
f"Major 5 species at initial temperature {initial_temperature} K (mole fraction):"
)
for idx in top5_idx:
print(f" {gas.species_names[idx]}: {gas.X[idx]:.6e}")
print(
f"Equilibration completed at temperature {gas.T} K, now switching to {inp['mech']} mechanism..."
)
gas = switch_mechanism(gas, inp["mech"], verbose=True)
# Run PFR with fixed time-temperature profile
print("----------------------PFR---------------------------------")
print(f"Using time step dt = {dt:.2e} s")
results = solve_fixed_time_temperature_profile_pfr(
inlet_composition=gas.X,
P_bar=inp["P_bar"],
mechanism=inp["mech"],
temperature_profile_time=temperature_profile,
dt=dt,
name="hot_spots",
verbose=2,
)
# Save results as CSV under results/{mech}/states_{time_ms}ms.csv
states = results["states"]
out_dir = _get_results_dir(inp)
out_dir.mkdir(parents=True, exist_ok=True)
time_ms = _get_time_ms_str(inp)
states_csv = out_dir / f"states_{time_ms}.csv"
states.save(states_csv, overwrite=True)
states_df = pd.read_csv(states_csv)
states_df["residence_time"] = states.t
states_df.to_csv(states_csv)
if inp["mech"] == "CRECK_Nobili2024.yaml":
states_df = pd.read_csv(states_csv)
mass_carbon_rates_by_class = results["mass_carbon_rates_by_class"]
reaction_rates_by_class = results["reaction_rates_by_class"]
for cls, rates in mass_carbon_rates_by_class.items():
states_df[f"mass_carbon_rate_{cls}"] = rates
for cls, rates in reaction_rates_by_class.items():
states_df[f"reaction_rate_{cls}"] = rates
states_df.to_csv(states_csv)
return results
def load_results(inp):
"""Load results from results/{mech}/states_{time_ms}ms.csv."""
out_dir = _get_results_dir(inp)
time_ms = _get_time_ms_str(inp)
states_csv = out_dir / f"states_{time_ms}.csv"
if not states_csv.exists():
raise FileNotFoundError(f"States file {states_csv} not found")
gas = ctutils.get_cached_solution(get_mechanism_path(inp["mech"]))
states = ct.SolutionArray(gas)
states.read_csv(states_csv)
results = {"states": states, "gas": gas}
return results
def _load_precursors_species(inp, results):
"""Load precursors species and names depending on the mechanism."""
gas = results["gas"]
# Caltech mechanism: A2, A2R5, A3, A3R5, A4
if inp["mech"] == "Caltech.yaml":
precursors_species = ["A2", "A2R5", "A3", "A3R5", "A4", "A4R5"]
elif inp["mech"] == "CRECK_Nobili2024.yaml":
species_groups = carbon_utils_nobili.generate_species_groups_nobili(gas)
precursors_species = species_groups["liquid_particles"]
else:
precursors_species = []
return precursors_species
def plot_results(inp, results):
# Unpack results
gas = results["gas"]
states = results["states"]
precursor_species = _load_precursors_species(inp, results)
out_dir = _get_results_dir(inp)
out_dir.mkdir(parents=True, exist_ok=True)
time_ms = _get_time_ms_str(inp)
# Find states.mass_carbon_rate_Nucleation peak time and temperature
if inp["mech"] == "CRECK_Nobili2024.yaml":
nucleation_peak_time = states.t[np.argmax(states.mass_carbon_rate_Nucleation)]
# Figure settings: ax1, ax2, ax3 share time (ms); ax4 uses temperature (K) separately
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(10, 8), sharex=False)
ax2.sharex(ax1)
ax3.sharex(ax1)
fig.suptitle(
f"Hot Spot temperature relaxation\n"
f"Exponential decay (tau={time_ms}) simulated with {inp['mech'].replace('.yaml', '')}"
)
# PLOT 1: Plot temperature profile
ax1.plot(states.t * 1000, states.T, label="Temperature", c="tab:red")
ax1.set_xlim(0, np.max(states.t) * 1000)
ax1.set_ylabel("Temperature (K)")
ax1.set_title("Temperature profile")
if inp["mech"] == "CRECK_Nobili2024.yaml":
ax1.axvline(
nucleation_peak_time * 1000,
color="tab:blue",
linestyle="--",
label="Nucleation peak",
)
ax1.legend(fontsize="small", loc="best", ncol=2)
# PLOT 2: Plot 5 major species by final mass fraction
species_names = gas.species_names
y_start = states.Y[0]
y_end = states.Y[-1]
top5_indices_start = np.argsort(y_start)[-5:][::-1]
top5_indices_end = np.argsort(y_end)[-5:][::-1]
for idx in top5_indices_start:
ax2.plot(
states.t * 1000,
states(species_names[idx]).Y,
label=species_names[idx],
linestyle="--",
)
for idx in top5_indices_end:
ax2.plot(
states.t * 1000, states(species_names[idx]).Y, label=species_names[idx]
)
mass_fraction_precursors = np.zeros(len(states.t))
for species in precursor_species:
if species in species_names:
mass_fraction_precursors += states(species).Y[:, 0]
ax2.plot(
states.t * 1000,
mass_fraction_precursors,
label="Liquid particles"
if inp["mech"] == "CRECK_Nobili2024.yaml"
else "Precursors",
)
ax2.set_ylabel("Mass fraction")
ax2.set_yscale("log")
ax2.set_ylim(1e-6, 1)
ax2.set_xlim(0, np.max(states.t) * 1000)
ax2.set_title("Species concentrations")
ax2.legend(fontsize="small", loc="best", ncol=2)
# PLOT 3: Plot mass carbon rates by class
if inp["mech"] == "CRECK_Nobili2024.yaml":
ax3.plot(
states.t * 1000, states.mass_carbon_rate_Nucleation, label="Nucleation"
)
ax3.plot(
states.t * 1000,
states.mass_carbon_rate_SurfaceGrowthPAH,
label="SurfaceGrowthPAH",
)
ax3.plot(
states.t * 1000,
states.mass_carbon_rate_SurfaceGrowthRadicals,
label="SurfaceGrowthRadicals",
)
ax3.plot(states.t * 1000, states.mass_carbon_rate_HACA, label="HACA")
ax3.plot(
states.t * 1000, states.mass_carbon_rate_Coalescence, label="Coalescence"
)
ax3.axvline(
nucleation_peak_time * 1000,
color="black",
linestyle="--",
label="Nucleation peak",
)
ax3.set_xlabel("Time (ms)")
ax3.set_ylabel("Mass carbon rate (kg/m^3/s)")
ax3.set_yscale("log")
ax3.set_ylim(1e-6, 1)
ax3.set_xlim(0, np.max(states.t) * 1000)
ax3.set_title("Mass carbon rates by class")
ax3.legend(fontsize="small", loc="best", ncol=2)
# PLOT 4: Histogram of nucleated mass vs temperature (mass per time step binned by T)
if inp["mech"] == "CRECK_Nobili2024.yaml":
dt = np.diff(np.concatenate(([0], states.t)))
mass_per_step = np.array(states.mass_carbon_rate_Nucleation) * dt
T_history = np.asarray(states.T)
T_min, T_max = T_history.min(), T_history.max()
n_bins = max(30, int((T_max - T_min) / 50)) # at least 30 bins, or ~50 K width
bin_edges = np.linspace(T_min, T_max, n_bins + 1)
mass_hist, _ = np.histogram(T_history, bins=bin_edges, weights=mass_per_step)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
bin_width = np.diff(bin_edges)
ax4.bar(
bin_centers,
mass_hist,
width=bin_width * 0.9,
align="center",
color="tab:blue",
edgecolor="gray",
)
ax4.set_xlabel("Temperature (K)")
ax4.set_ylabel("Nucleated mass (kg/m³)")
ax4.set_title("Nucleation mass distribution vs temperature")
ax4.set_xlim(T_min, T_max)
fig.tight_layout()
fig.savefig(out_dir / f"results_{time_ms}.png", bbox_inches="tight")
if __name__ == "__main__":
print("---------------- Running hot spots simulation ----------------")
inp = _load_inputs()
if SKIP_SIMULATION:
states_path = _get_results_dir(inp) / f"states_{_get_time_ms_str(inp)}.csv"
print(f"Skipping simulation and loading results from {states_path}")
results = load_results(inp)
else:
results = run(inp)
results = load_results(inp)
plot_results(inp, results)