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:

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

Standard ranges by system:

Step 3: Add Stress Cases

Stress cases specific to Gibbs minimization:

1. Trace species challenge: Initial guess with < 1e-8 mol for a species that should be significant at equilibrium (tests step sizing)
2. Near-complete reaction: Conditions where one reactant is nearly consumed (tests handling of near-zero mole numbers)
3. Inert dilution: 0%, 50%, 90%, 99% N2 dilution (tests scaling with dilution)
4. Adiabatic vs isothermal: Same feed at both modes (tests energy balance coupling)
5. Phase boundary: Conditions near condensation/solidification of products
6. High-temperature dissociation: T > 2000 K where minor species become significant
7. Ill-conditioned Jacobian: Systems with species spanning 10+ orders of magnitude

Step 4: Define Metrics

Every reactor benchmark case must record:

Step 5: Reference Solutions

For validation, compare against:

For each system, obtain at least 5 reference points at different conditions.

Step 6: Generate Config File

Output benchmark_config.json:

{
  "benchmark_id": "gibbs_reactor_2026",
  "created": "2026-03-31",
  "solver_variants": [
    {
      "name": "baseline",
      "description": "GibbsReactor with default settings",
      "settings": {"minIterations": 100, "adaptiveStepSize": false}
    },
    {
      "name": "optimized",
      "description": "GibbsReactor with adaptive step + min_iter=3",
      "settings": {"minIterations": 3, "adaptiveStepSize": true}
    }
  ],
  "reaction_systems": [
    {
      "name": "claus_direct",
      "feed_components": {"H2S": 0.10, "oxygen": 0.05, "nitrogen": 0.85},
      "expected_products": ["H2O", "S8", "SO2"],
      "T_range_K": [400, 1200],
      "P_range_bara": [1, 50],
      "n_T": 10,
      "n_P": 5,
      "perturbations": ["stoichiometric", "excess_H2S", "excess_O2", "diluted"],
      "reference_source": "NASA CEA"
    }
  ],
  "stress_cases": {
    "trace_species": 10,
    "near_complete_reaction": 10,
    "high_dilution": 10,
    "ill_conditioned": 5
  }
}

Python Execution Pattern for Reactor Benchmarks

from tools.neqsim_bootstrap import get_jneqsim
jneqsim = get_jneqsim()

SystemSrkEos = jneqsim.thermo.system.SystemSrkEos
ProcessSystem = jneqsim.process.processmodel.ProcessSystem
Stream = jneqsim.process.equipment.stream.Stream
GibbsReactor = jneqsim.process.equipment.reactor.GibbsReactor

def run_gibbs_reactor_case(case, mode="isothermal"):
    """Run a single Gibbs reactor case and return metrics."""
    import time

    fluid = SystemSrkEos(case["T_K"], case["P_bara"])
    for comp, frac in case["feed"].items():
        fluid.addComponent(comp, frac)
    fluid.setMixingRule("classic")
    fluid.setMultiPhaseCheck(True)

    feed = Stream("feed", fluid)
    feed.setFlowRate(1000.0, "kg/hr")
    feed.run()

    reactor = GibbsReactor("gibbs", feed)
    if mode == "adiabatic":
        reactor.setEnergyMode(
            jneqsim.process.equipment.reactor.GibbsReactor.EnergyMode.ADIABATIC)

    # Configure solver
    reactor.setMinIterations(3)
    reactor.setUseAdaptiveStepSize(True)

    process = ProcessSystem()
    process.add(feed)
    process.add(reactor)

    t0 = time.perf_counter_ns()
    try:
        process.run()
        elapsed_ms = (time.perf_counter_ns() - t0) / 1e6

        out_fluid = reactor.getOutletStream().getFluid()
        return {
            "converged": reactor.hasConverged(),
            "iterations": reactor.getActualIterations(),
            "cpu_time_ms": round(elapsed_ms, 3),
            "outlet_T_K": round(float(out_fluid.getTemperature()), 2),
            "n_phases": int(out_fluid.getNumberOfPhases()),
            "error": None
        }
    except Exception as e:
        elapsed_ms = (time.perf_counter_ns() - t0) / 1e6
        return {
            "converged": False,
            "cpu_time_ms": round(elapsed_ms, 3),
            "error": str(e)
        }

Checklist

• At least 3 reaction systems selected
• Temperature and pressure ranges cover the relevant equilibrium region
• Stress cases include trace species and ill-conditioned scenarios
• Reference solutions from NASA CEA or published data are identified
• Element balance verification is included in metrics
• Both isothermal and adiabatic modes are tested
• Total case count is between 200 and 1000
• Random seed is fixed for composition perturbations