design_reactor_benchmark

Skill: Design Reactor / Chemical Equilibrium Benchmark

Purpose

Create a structured test matrix for evaluating Gibbs energy minimization solvers across reaction systems, conditions, and difficulty levels.

When to Use

• Starting a paper on chemical equilibrium algorithms (Gibbs reactor)
• Benchmarking Jacobian formulations or solver improvements
• Comparing Gibbs minimization against reference solutions (JANAF/NASA CEA)
• Evaluating convergence for reactive systems with trace species

Benchmark Design Procedure

Step 1: Select Reaction Systems

Choose test systems that span different thermochemical challenges:

System	Feed Components	Key Products	Nc	Challenge
Claus (direct)	H2S, O2, N2	H2O, S8, SO2	5+	Sulfur precipitation, trace species
Claus (two-stage)	H2S, O2 → SO2; then H2S + SO2	S, H2O	5+	Multi-reactor, intermediate species
Methane combustion	CH4, O2, N2	CO2, H2O, CO, NO	7+	High temperature, many products
Steam methane reforming	CH4, H2O	CO, H2, CO2	5	Endothermic, equilibrium-limited
Water-gas shift	CO, H2O	CO2, H2	4	Temperature-sensitive equilibrium
Ammonia synthesis	N2, H2	NH3	3	High pressure, sparse products
CO2 hydrogenation	CO2, H2	CH3OH, H2O, CO	5	Catalyst-dependent selectivity
Iron sulfide corrosion	Fe, H2S	FeS, H2	4	Solid product formation
Sour gas sweetening	H2S, CO2, CH4, MEA	Various	8+	Acid gas + amine chemistry

Step 2: Define Condition Sweeps

For each system, define the parameter space:

import numpy as np

def generate_reactor_conditions(system):
    """Generate test conditions for a reaction system."""
    cases = []

    # Temperature sweep (most important for equilibrium)
    T_values = np.linspace(system["T_min_K"], system["T_max_K"], system["n_T"])

    # Pressure sweep
    P_values = np.logspace(
        np.log10(system["P_min_bara"]),
        np.log10(system["P_max_bara"]),
        system["n_P"]
    )

    # Feed composition perturbations
    for T in T_values:
        for P in P_values:
            # Stoichiometric feed
            cases.append({"T_K": float(T), "P_bara": float(P),
                          "feed": system["stoichiometric_feed"],
                          "label": "stoichiometric"})
            # Excess reactant A
            cases.append({"T_K": float(T), "P_bara": float(P),
                          "feed": system["excess_A_feed"],
                          "label": "excess_A"})
            # Excess reactant B
            cases.append({"T_K": float(T), "P_bara": float(P),
                          "feed": system["excess_B_feed"],
                          "label": "excess_B"})
    return cases