How to do sensitvity analysis using cantera

[potta sai chaitanya]

I am performing microkinetic modeling (MKM) of CO₂ methanation using Cantera, and I want to calculate the sensitivity coefficients of each elementary reaction with respect to the partial pressure of CH₄.”
=========================================================================

CODE

=========================================================================from logging import exception

import csv

import cantera as ct

import numpy as np

import matplotlib.pyplot as plt

def composition_solver(tc):

    cm = 0.01

    minute = 60.0

    dia = 9 / 1000

    porosity = 0.2

    vol_flow_rate = (250 * 10**(-6)) / 60

    velocity = (4 * vol_flow_rate) / (np.pi * dia**2)

    area = np.pi * (dia**2) / 4

    m = 2.5  # mg

    bulkdensity = 2106453

    length = (m * 0.001) / (bulkdensity * area)

    cat_area_per_vol = 7.23 * 10**6

    NReactors = 1002

    file = 'Methanation_Ni.yaml'

    t = tc + 273.15

    surf = ct.Interface(file, 'Ni-surface')

    surf.TP = t, ct.one_atm

    gas = surf.adjacent['gas']

    gas.TPX = t, ct.one_atm, 'CO2:0.08, H2:0.32, N2:0.60'

    rlen = length / (NReactors - 1)

    rvol = area * rlen

    cat_area = cat_area_per_vol * rvol

    mass_flow_rate = velocity * gas.density * area * porosity

    r = ct.IdealGasReactor(gas, energy='off')

    r.volume = rvol

    upstream = ct.Reservoir(gas)

    downstream = ct.Reservoir(gas)

    rsurf = ct.ReactorSurface(surf, r, A=cat_area)

    m = ct.MassFlowController(upstream, r, mdot=mass_flow_rate)

    v = ct.PressureController(r, downstream, primary=m, K=1e-5)

    sim = ct.ReactorNet([r])

    sim.max_err_test_fails = 20

    sim.rtol = 1.0e-8

    sim.atol = 1.0e-18

    co2_in = gas['CO2'].X[0]

    for n in range(NReactors):

        sim.reinitialize()

        upstream.syncState()

        sim.advance_to_steady_state()

        if n == NReactors - 1:

            dist = n * rlen * 1e3

            co2_out = r.thermo['CO2'].X[0]

            conversion = ((co2_in - co2_out) / co2_in)*100

            mole_fractions = list(r.thermo.X)

            coverages = list(rsurf.kinetics.coverages)

            rxn_rates = list(rsurf.kinetics.net_rates_of_progress)

    return (dist, tc, r.thermo.P / ct.one_atm, conversion,

            mole_fractions, coverages, rxn_rates,

            gas.species_names, surf.species_names)

temperatures =[]

all_results = []

conversions = []

for temp in np.linspace(200,501,57):

    print(f"\nSimulating at {temp}°C...")

    try:

      (dist, temp_out, pressure, conversion, mole_fractions,

       coverages, rxn_rates, gas_names, surf_names) = composition_solver(temp)

      conversions.append(conversion)

      result = ([dist, temp_out, pressure, conversion] +

              mole_fractions + coverages + rxn_rates)

      all_results.append(result)

      temperatures.append(temp)

      print(f"Final CO2 conversion at {temp}°C: {conversion:.3f}")

    except Exception as e:

        print(f"An error occurred at temperature {temp}°C: {e}")

        continue

output_filename = 'methanation_results_all_temps.csv'

with open(output_filename, 'w', newline='') as csvfile:

    writer = csv.writer(csvfile)

    headers = (['Distance (mm)', 'T (C)', 'P (atm)', 'CO2 Conversion'] +

               [f'X_{s}' for s in gas_names] +

               [f'Coverage_{s}' for s in surf_names] +

               [f'Rate_{i}' for i in range(len(rxn_rates))])

    writer.writerow(headers)

    writer.writerows(all_results)

print(f"\nAll results saved to '{output_filename}'")

print("\nSummary of CO2 Conversions:")

for temp, conv in zip(temperatures, conversions):

    print(f"Temperature: {temp}°C, CO2 Conversion: {conv:.3f}")

# Plot conversions vs temperature

plt.figure(figsize=(8, 6))

plt.plot(temperatures, conversions, 'bo-', label='CO2 Conversion')

plt.xlabel('Temperature (°C)')

plt.ylabel('CO2 Conversion')

plt.title('CO2 Conversion vs Temperature')

plt.grid(True)

plt.legend()

plt.savefig('conversion_vs_temperature.png')

plt.show()