In [ ]:

    

import matplotlib.pyplot as plt
from collections import defaultdict
from chempy import ReactionSystem, Species, Reaction
from chempy.units import to_unitless, SI_base_registry, rescale, default_constants as const, default_units as u
from chempy.kinetics.ode import (
    _create_odesys as create_odesys,
)
%matplotlib inline


    

In [ ]:

    

ads_token = '(ads)'
phases = {'(aq)': 0, ads_token: 1}
ads_comp = -1
vacancy = Species('vacancy(ads)', composition={ads_comp: 1},
                  phase_idx=phases[ads_token], latex_name='vacancy(ads)')

def substance_factory(name):
    if name == vacancy.name:
        return vacancy
    s = Species.from_formula(name, phases=phases)
    if name.endswith(ads_token):
        s.composition[ads_comp] = 1
    return s

rsys = ReactionSystem.from_string("""
vacancy(ads) + H2O2 -> H2O2(ads); 'k_ads_H2O2'
H2O2(ads) + vacancy(ads) -> 2 OH(ads); 'k_split_H2O2ads'
OH(ads) + H2O2 -> H2O + HO2 + vacancy(ads); 'k_abst_OHads_H2O2'
HO2 + HO2 -> O2 + H2O2; 'k_disprop_HO2'
""", substance_factory=substance_factory)
rsys


    

In [ ]:

    

ureg = SI_base_registry.copy()
ureg['length'] = u.dm
odesys, odesys_extra = create_odesys(rsys, unit_registry=ureg)


    

In [ ]:

    

powder_mass = 42*u.g
specific_surf_area = 35*u.m**2/u.g  # e.g. from BET
sample_volume = 250*u.cm3
surf_volume_ratio = powder_mass*specific_surf_area/sample_volume

adsorption_cross_section = 0.162*u.nm**2  # N2 at 77K on active carbon
conc_vacancy = rescale(surf_volume_ratio/(adsorption_cross_section * const.Avogadro_constant), u.molar)
conc_vacancy


    

In [ ]:

    

c0 = defaultdict(lambda: 0*u.M, H2O2=.2*u.M, **{'vacancy(ads)': conc_vacancy})

tend = 3600*u.s
result, result_extra = odesys_extra['unit_aware_solve'](
    (1e-7*u.s, tend), c0,
    dict(
        k_ads_H2O2=.1/u.M/u.s,
        k_split_H2O2ads=.2/u.M/u.s,
        k_abst_OHads_H2O2=.3/u.M/u.s,
        k_disprop_HO2=.4/u.M/u.s
    ), integrator='cvode'
)


    

In [ ]:

    

fig, axes = plt.subplots(1, 2, figsize=(14, 4))
result.plot(ax=axes[0])
result.plot(ax=axes[1], xscale='log', yscale='log')
for ax in axes:
    #ax.set_ylabel('Concentration / M')
    #ax.set_xlabel('Time / min')
    ax.set_ylim(to_unitless([1e-6*u.molar, c0['H2O2']], result_extra['unit_conc']))
    ax.set_xlim(to_unitless([1*u.s, tend], result_extra['unit_time']))