—
Tube furnace – N₂ temperature profiles vs. Mei et al. (2019).
Reproduces Fig. 2 of Mei et al. (2019) Combust. Flame 200, 14–25 (doi: 10.1016/j.combustflame.2018.11.010): axial N₂ gas temperature profiles for 0, 7 and 13 L/min in a tube furnace at T_wall = 1673 K. Digitised reference points are overlaid for comparison.
This module also implements the ctwrap simulation interface used
by parametric sweep scripts.
ctwrap interface#
DEFAULTS_YAMLExample default parameters in YAML format (see also
defaults.parametersinmodels/tube_furnace/tube_furnace.yaml). Parsed bydefaults().run(parameters, **kwargs) -> dictBuild a
LagrangianPFRReactor, advance through the tube, and return a one-row result dict. Passsolver=dict to control integrator settings (e.g.use_preconditioner).
Output columns#
diameter_mmtube inner diameter [mm]T_outlet_Cgas temperature at tube exit [degC]X_H2H2 mole fraction at outletX_C2H2C2H2 mole fraction at outletn_stepsCVODE steps takenelapsed_swall-clock time [s]
import io
import sys
import time
from pathlib import Path
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from ruamel.yaml import YAML
# Make sure bloc is importable when run from examples/
_HERE = Path(__file__).parent.absolute()
RESULTS_DIR = _HERE / "Results"
sys.path.insert(0, str(_HERE.parent))
from bloc.reactors import LagrangianPFRNetwork, LagrangianPFRReactor # noqa: E402
from bloc.utils import get_mechanism_path, get_validation_data # noqa: E402
# ---------------------------------------------------------------------------
# Example default parameters (YAML-like; parsed by defaults() for ctwrap)
# ---------------------------------------------------------------------------
DEFAULTS_YAML = """
parameters:
mechanism: gri30.yaml # must include transport (e.g. Mix); CRECK kinetics has none
diameter: 0.100 # [m]
total_length: 1.0 # [m]
T_wall_K: 1873.15 # [K] 1600 degC
T_ambient_K: 300.0 # [K]
entry_leg: 0.15 # [m] cold leg at inlet
exit_leg: 0.15 # [m] cold leg at outlet
entry_zone: 0.20 # [m] ramp up
plateau_zone: 0.70 # [m]
kappa_grey: 0.5 # [1/m]
flow_slm: 0.0 # [SLM]
composition: "N2:1"
T_inlet_K: 298.15 # [K]
P_inlet_Pa: 101325.0 # [Pa]
"""
def defaults():
"""Return default parameters for ctwrap (parsed from :data:`DEFAULTS_YAML`)."""
yaml_loader = YAML(typ="safe")
return yaml_loader.load(io.StringIO(DEFAULTS_YAML))
# ---------------------------------------------------------------------------
# Helper: trapezoidal wall temperature profile
# ---------------------------------------------------------------------------
def _wall_T_fn(params):
"""Return a callable ``T_wall(x) -> float [K]`` for the given parameters."""
T_wall_K = params["T_wall_K"]
T_ambient_K = params["T_ambient_K"]
entry_zone = params["entry_zone"]
plateau_zone = params["plateau_zone"]
total_length = params["total_length"]
entry_leg = params.get("entry_leg", 0.0)
exit_leg = params.get("exit_leg", 0.0)
ramp_start = entry_leg
plateau_start = entry_leg + entry_zone
plateau_end = plateau_start + plateau_zone
ramp_end = total_length - exit_leg
def fn(x):
x = min(max(x, 0.0), total_length)
if x < entry_leg:
return T_ambient_K
if x >= total_length - exit_leg:
return T_ambient_K
if x < plateau_start:
frac = (x - ramp_start) / entry_zone if entry_zone > 0 else 1.0
return T_ambient_K + (T_wall_K - T_ambient_K) * frac
if x <= plateau_end:
return T_wall_K
exit_zone_len = ramp_end - plateau_end
frac = (x - plateau_end) / exit_zone_len if exit_zone_len > 0 else 1.0
return T_wall_K + (T_ambient_K - T_wall_K) * frac
return fn
# ---------------------------------------------------------------------------
# Shared setup
# ---------------------------------------------------------------------------
def _setup(parameters, wall_T_fn=None, solver=None):
"""Unpack *parameters*, build reactor and network.
Parameters
----------
parameters : dict
Physics parameters (see :func:`defaults` and :data:`DEFAULTS_YAML`).
wall_T_fn : callable or None, optional
Override the trapezoidal wall temperature profile. When ``None``
(default), :func:`_wall_T_fn` is built from *parameters*. Pass a
custom ``T_wall(x) -> float [K]`` to use measured wall data directly.
solver : dict or None, optional
LagrangianPFRNetwork solver settings (from the ``solver:`` section of
``tube_furnace.yaml``). Recognised keys: ``use_preconditioner``
(bool), ``rtol`` (float), ``atol`` (float). Defaults apply when
``None`` or when individual keys are absent.
Returns
-------
network : LagrangianPFRNetwork
Ready to advance (not yet run).
"""
mechanism_path = get_mechanism_path(parameters["mechanism"])
diameter = float(parameters["diameter"]) # [m]
total_length = float(parameters["total_length"]) # [m]
kappa_grey = float(parameters["kappa_grey"]) # [1/m]
flow_slm = float(parameters["flow_slm"]) # [SLM]
composition = parameters["composition"] # e.g. "N2:1"
T_inlet_K = float(parameters["T_inlet_K"]) # [K]
P_inlet_Pa = float(parameters["P_inlet_Pa"]) # [Pa]
gas_inlet = ct.Solution(mechanism_path)
gas_inlet.TPX = 273.15, ct.one_atm, composition
mass_flow_rate = flow_slm * 1e-3 / 60.0 * gas_inlet.density_mass
gas_inlet.TPX = T_inlet_K, P_inlet_Pa, composition
A_cs = np.pi / 4.0 * diameter**2
gas = ct.Solution(mechanism_path)
gas.TPX = T_inlet_K, P_inlet_Pa, composition
reactor = LagrangianPFRReactor(
gas,
clone=False,
diameter=diameter,
kappa_grey=kappa_grey,
mass_flow_rate=mass_flow_rate,
)
reactor.volume = A_cs * total_length
if wall_T_fn is None:
wall_T_fn = _wall_T_fn(parameters)
_s = solver or {}
network = LagrangianPFRNetwork(
reactor,
total_length,
mass_flow_rate,
wall_T_fn,
use_preconditioner=bool(_s.get("use_preconditioner", True)),
rtol=float(_s.get("rtol", 1e-4)),
atol=float(_s.get("atol", 1e-12)),
)
return network
# ---------------------------------------------------------------------------
# ctwrap run function
# ---------------------------------------------------------------------------
def run(parameters, wall_T_fn=None, solver=None, **kwargs):
"""Run one tube furnace Lagrangian simulation.
Called by ``ctwrap.SimulationHandler`` for each point in the parametric
sweep. The ``parameters`` dict is pre-populated from ``defaults`` and
then overridden with the current strategy value.
Parameters
----------
parameters : dict
Physics parameters (see :func:`defaults` and :data:`DEFAULTS_YAML`).
wall_T_fn : callable or None, optional
Override the trapezoidal wall temperature profile (see :func:`_setup`).
solver : dict or None, optional
LagrangianPFRNetwork solver settings. Passed through to
:func:`_setup`. See :data:`tube_furnace.yaml` ``solver:`` section.
Returns
-------
dict
``{"sim": ct.ReactorNet, "sim_extra": dict, "res_dic": dict}``
* ``sim`` — solved ReactorNet for topology-based extraction.
* ``sim_extra`` — Lagrangian scalars not recoverable from *sim* post-advance:
``E_conv_J``, ``E_rad_J``, ``t_res_s``, ``mass_flow_rate_kg_s``,
``n_steps``, ``elapsed_s``, ``diameter_mm``.
* ``res_dic`` — legacy flat dict (backward-compat; deprecated in Phase 3).
"""
network = _setup(parameters, wall_T_fn=wall_T_fn, solver=solver)
t0 = time.perf_counter()
network.advance(record_profiles=True)
elapsed = time.perf_counter() - t0
t_res = float(network.net.time)
e_conv = float(network.E_conv)
e_rad = float(network.E_rad)
mdot = float(network.mass_flow_rate)
sim_extra = {
"E_conv_J": e_conv,
"E_rad_J": e_rad,
"t_res_s": t_res,
"mass_flow_rate_kg_s": mdot,
"n_steps": int(network.n_steps),
"elapsed_s": float(elapsed),
"diameter_mm": float(network.reactor.diameter * 1e3),
"Fo_D_min": float(network.Fo_D_min),
"h_conv_avg_W_m2_K": float(network.h_conv_avg),
# LagrangianPFRNetwork object — exposes draw/plot methods for the
# Calculation Note Detail sheet via SIM_FIGURE_GENERATORS.
"network": network,
}
return {
"sim": network.net,
"sim_extra": sim_extra,
"res_dic": _build_result(network, elapsed),
}
def get_wall_profile(parameters, n_points: int = 200) -> tuple[np.ndarray, np.ndarray]:
"""Return (x_cm, T_wall_K) for the trapezoidal wall profile (no flow)."""
total_length = float(parameters["total_length"])
wall_T_fn = _wall_T_fn(parameters)
x = np.linspace(0, total_length, n_points)
T_wall_K = np.array([wall_T_fn(xi) for xi in x])
return x * 1e2, T_wall_K
def run_with_profiles(
parameters,
wall_T_fn=None,
solver=None,
) -> tuple[pd.Series, np.ndarray, np.ndarray, np.ndarray, np.ndarray, float]:
"""Run simulation and return result plus (x_cm, T_wall_C, T_gas_C, Y_Cs, Y_C_feed).
Parameters
----------
parameters : dict
Physics parameters (see :func:`defaults` and :data:`DEFAULTS_YAML`).
wall_T_fn : callable or None, optional
Override the trapezoidal wall temperature profile (see :func:`_setup`).
Returns
-------
result : pandas.Series
Result metrics for this case.
x_cm : np.ndarray
Axial positions [cm].
T_wall_C : np.ndarray
Wall temperature profile [°C].
T_gas_C : np.ndarray
Gas temperature profile [°C].
Y_Cs : np.ndarray
Solid carbon mass fraction profile along the tube.
Y_C_feed : float
Inlet (feed) total carbon elemental mass fraction. Use to compute
carbon yield as Y_Cs / Y_C_feed when Y_C_feed > 0.
"""
network = _setup(parameters, wall_T_fn=wall_T_fn, solver=solver)
Y_C_feed = network.reactor.phase.elemental_mass_fraction("C")
t0 = time.perf_counter()
network.advance(verbose=True, record_profiles=True)
elapsed = time.perf_counter() - t0
x_cm = network.x_arr * 1e2
T_wall_C = network.T_wall_arr - 273.15
T_gas_C = network.T_gas_arr - 273.15
Y_Cs = network.Y_Cs_arr
return (
pd.Series(_build_result(network, elapsed)),
x_cm,
T_wall_C,
T_gas_C,
Y_Cs,
Y_C_feed,
)
def _build_result(network, elapsed):
"""Build result dict from a completed *network* run."""
gas_out = network.reactor.phase
def _X(species):
return (
gas_out.X[gas_out.species_index(species)]
if species in gas_out.species_names
else 0.0
)
t_res_actual = network.net.time
P_conv_avg = network.E_conv / t_res_actual if t_res_actual > 0 else 0.0
P_rad_avg = network.E_rad / t_res_actual if t_res_actual > 0 else 0.0
# Specific Energy Output: wall-to-gas energy transfer per unit mass of gas.
# SEO = (E_conv + E_rad) / (mdot * t_res) [MJ/kg]
m_total = network.mass_flow_rate * t_res_actual
SEO_MJ_kg = (
(network.E_conv + network.E_rad) / m_total / 1e6
if m_total > 0
else float("nan")
)
return {
"diameter_mm": network.reactor.diameter * 1e3,
"T_outlet_C": gas_out.T - 273.15,
"X_CH4": _X("CH4"),
"X_H2": _X("H2"),
"X_C2H2": _X("C2H2"),
"X_C2H4": _X("C2H4"),
"X_C6H6": _X("C6H6"),
"n_steps": network.n_steps,
"elapsed_s": elapsed,
"E_conv_J": network.E_conv,
"E_rad_J": network.E_rad,
"P_conv_avg_W": P_conv_avg,
"P_rad_avg_W": P_rad_avg,
"SEO_MJ_kg": SEO_MJ_kg,
"Fo_D_min": network.Fo_D_min,
"h_conv_avg_W_m2_K": network.h_conv_avg,
}
# ---------------------------------------------------------------------------
Demo / validation cells — guarded so ctwrap can import this file as a pure module without executing the Mei et al. validation run.
if __name__ == "__main__":
# Loading Mei et al. (2019) validation data
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Geometry parameters and digitised reference curves are stored in a single
# YAML file under ``tests/validation_data/``.
# :func:`~bloc.utils.get_validation_data` resolves the path from the repo root
# so it works identically from tests and from example scripts.
_mei_data = get_validation_data("mei2019_fig2.yaml")
_MEI_PARAMS = _mei_data["parameters"]
_REF = _mei_data["reference"]
# %%
# Building the wall temperature function
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# The measured 0 L/min static profile is the best available description of the
# actual furnace wall temperature. It is used directly as the boundary
# condition for both flow simulations so that the model wall matches the
# reference exactly.
_x_wall_m = np.array(_REF["static_wall"]["x_cm"]) * 1e-2
_T_wall_K_ref = np.array(_REF["static_wall"]["T_K"])
x_wall_cm = _x_wall_m * 1e2
T_wall_K = _T_wall_K_ref
def _mei_wall_T_fn(x):
"""Linearly interpolate the measured static wall temperature profile."""
return float(
np.interp(
x,
_x_wall_m,
_T_wall_K_ref,
left=_T_wall_K_ref[0],
right=_T_wall_K_ref[-1],
)
)
# %%
# Running the simulations
# ~~~~~~~~~~~~~~~~~~~~~~~
# Both flow cases use the measured wall profile as boundary condition.
_, x_7_cm, _, T_7_gas_C, Y_Cs_7, Y_C_feed = run_with_profiles(
{**_MEI_PARAMS, "flow_slm": 7.0}, wall_T_fn=_mei_wall_T_fn
)
_, x_13_cm, _, T_13_gas_C, Y_Cs_13, _ = run_with_profiles(
{**_MEI_PARAMS, "flow_slm": 13.0}, wall_T_fn=_mei_wall_T_fn
)
# Carbon yield profile = solid C mass fraction / initial feed carbon mass fraction
carbon_yield_7 = np.divide(
Y_Cs_7, Y_C_feed, out=np.zeros_like(Y_Cs_7), where=Y_C_feed > 0
)
carbon_yield_13 = np.divide(
Y_Cs_13, Y_C_feed, out=np.zeros_like(Y_Cs_13), where=Y_C_feed > 0
)
# %%
# Comparing simulated profiles with Mei et al. (2019) reference
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# The model (forced convection only) reproduces the key trends:
# heating ramp, plateau, and thermal inertia at higher flow rates.
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(x_wall_cm, T_wall_K, color="black", lw=2, label="Model – static (wall)")
ax.plot(
x_7_cm, T_7_gas_C + 273.15, color="red", lw=2, ls="--", label="Model – 7 L/min"
)
ax.plot(
x_13_cm,
T_13_gas_C + 273.15,
color="blue",
lw=2,
ls="-.",
label="Model – 13 L/min",
)
ax.plot(
_REF["static_wall"]["x_cm"],
_REF["static_wall"]["T_K"],
"ks",
ms=5,
label="Mei 2019 – static",
)
ax.plot(
_REF["flow_7lpm"]["x_cm"],
_REF["flow_7lpm"]["T_K"],
"ro",
ms=5,
mfc="none",
label="Mei 2019 – 7 L/min",
)
ax.plot(
_REF["flow_13lpm"]["x_cm"],
_REF["flow_13lpm"]["T_K"],
"b^",
ms=5,
mfc="none",
label="Mei 2019 – 13 L/min",
)
ax.set_xlabel("Axial position [cm]")
ax.set_ylabel("Temperature [K]")
ax.set_title("N₂ temperature profiles – Mei et al. (2019) Fig. 2")
ax.legend(loc="lower center", ncol=2, fontsize=7.5)
ax.set_xlim(0, 140)
ax.set_ylim(200, 1900)
ax.grid(True, alpha=0.3)
plt.tight_layout()
# Save figure to Results folder next to this script
RESULTS_DIR.mkdir(exist_ok=True)
fig_path = RESULTS_DIR / "tube_furnace_mei2019_fig2.png"
fig.savefig(fig_path, dpi=150, bbox_inches="tight")
print(f"Figure saved to {fig_path}")
# Second figure: carbon yield profile along the tube (when feed contains carbon)
if Y_C_feed > 0:
fig_y, ax_y = plt.subplots(figsize=(8, 5))
ax_y.plot(
x_7_cm,
carbon_yield_7 * 100.0,
color="red",
lw=2,
ls="--",
label="7 L/min",
)
ax_y.plot(
x_13_cm,
carbon_yield_13 * 100.0,
color="blue",
lw=2,
ls="-.",
label="13 L/min",
)
ax_y.set_xlabel("Axial position [cm]")
ax_y.set_ylabel("Carbon yield [%]")
ax_y.set_title("Carbon yield along the tube")
ax_y.set_xlim(0, 140)
ax_y.set_ylim(0, 100)
ax_y.grid(True, alpha=0.3)
ax_y.legend(loc="lower right")
plt.tight_layout()
yield_path = RESULTS_DIR / "tube_furnace_carbon_yield.png"
fig_y.savefig(yield_path, dpi=150, bbox_inches="tight")
print(f"Carbon-yield figure saved to {yield_path}")