TUTORIAL_composition-dependent surface adsorption.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Composition-dependent surface adsorption\n",
    "\n",
    "One may have a system whereby the adsorption of a gas-phase species is dependent on the composition of the particle or film surface. MultilayerPy can facilitate this in the model building and optimisation process. \n",
    "\n",
    "Let's build the same simple reactive system used in other tutorials: `oleic acid + ozone --> products`\n",
    "\n",
    "Now, let's say that ozone does not adsorb to the product molecules as readily as it does to pure oleic acid. This is accounted for in a separate `alpha_s0` parameter for ozone on a product surface (the surface accommodation coefficient on a free surface). Composition-dependent surface adsorption is defined when defining specific model components.\n",
    "\n",
    "Surface-dependent surface adsorption is defined below:\n",
    "\n",
    "$$\n",
    "\\alpha_{s,0,X_{i}} = \\sum_{j} \\alpha_{s,0,X_{i},Y_{j}}\\theta_{ss,Y_{j}}\n",
    "$$\n",
    "\n",
    "Where $\\alpha_{s,0,X_{i}}$ is the overall surface accommodation coefficient for gas $X_{i}$, $\\alpha_{s,0,X_{i},Y_{j}}$ is the surface accommodation coefficent of $X_{i}$ on a surface of pure $Y_{j}$ and $\\theta_{ss,Y_{j}}$ is the fractional surface area of $Y_{j}$.\n",
    "\n",
    "The initial code is the same as in other tutorial examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import necessary packages\n",
    "import numpy as np\n",
    "import multilayerpy\n",
    "import multilayerpy.build as build \n",
    "import multilayerpy.simulate as simulate\n",
    "import multilayerpy.optimize as optimize\n",
    "\n",
    "# check the version\n",
    "print(multilayerpy.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import classes needed to build the model\n",
    "from multilayerpy.build import ModelType, ReactionScheme, DiffusionRegime, ModelComponent\n",
    "\n",
    "# define the model type (spherical geometry)\n",
    "mod_type = ModelType('km-sub','spherical')\n",
    "\n",
    "# build the reaction tuple list, in this case only 1 tuple in the list (for 1 reaction)\n",
    "# component 1 (oleic acid) reacts with component 2 (ozone)\n",
    "reaction_tuple_list = [(1,2)]\n",
    "\n",
    "# build the product tuple list, only component 3 (products) is a product\n",
    "# a tuple with a single value inside is defined (value,)\n",
    "product_tuple_list = [(3,)]\n",
    "\n",
    "# now construct the reaction scheme\n",
    "# we can give it a name and define the nuber of components as below\n",
    "reaction_scheme = ReactionScheme(mod_type,name='Oleic acid ozonolysis',\n",
    "                                                   reactants=reaction_tuple_list,\n",
    "                                                products=product_tuple_list)\n",
    "\n",
    "# let's print out a representation of the reaction scheme\n",
    "reaction_scheme.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adding in model components\n",
    "# oleic acid\n",
    "OA = ModelComponent(1,reaction_scheme,name='Oleic acid')\n",
    "\n",
    "#================================\n",
    "\n",
    "# ozone, declare that it is a gas\n",
    "# ***HERE IS WHERE COMPOSITION-DEPENDENT ADSORPTION IS DEFINED***\n",
    "O3 = ModelComponent(2,reaction_scheme,gas=True,name='Ozone',comp_dependent_adsorption=True) \n",
    "\n",
    "#================================\n",
    "\n",
    "# products\n",
    "prod = ModelComponent(3,reaction_scheme, name='Reaction products')\n",
    "\n",
    "# collect into a dictionary\n",
    "model_components_dict = {'1':OA,\n",
    "                        '2':O3,\n",
    "                        '3':prod}\n",
    "\n",
    "# making the diffusion dictionary\n",
    "diff_dict = {'1' : None,\n",
    "             '2': None,\n",
    "             '3':None}  \n",
    "\n",
    "# make diffusion regime\n",
    "diff_regime = DiffusionRegime(mod_type,model_components_dict,diff_dict=diff_dict)\n",
    "\n",
    "# call it to build diffusion code ready for the builder\n",
    "diff_regime()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now construct the model and identify the required parameters\n",
    "\n",
    "# import ModelBuilder class\n",
    "from multilayerpy.build import ModelBuilder\n",
    "\n",
    "# create the model object\n",
    "model = ModelBuilder(reaction_scheme,model_components_dict,diff_regime)\n",
    "\n",
    "# build the model. Will save a file, don't include the date in the model filename\n",
    "model.build(date=False)\n",
    "\n",
    "# print out the parameters required for the model to run\n",
    "print(model.req_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Additional required input parameters\n",
    "\n",
    "Notice that `alpha_s_0_2_1` and `alpha_s_0_2_3` are included in the list of required parameters. `alpha_s_0_2_1` stands for \"surface adsorption coefficient of `2` on pure `1` and has a value between 0 and 1\". These parameters would not be required if composition-dependent adsorption was not imposed for component `2`.\n",
    "\n",
    "These parameters need to be supplied when running the model. For the purpose of this tutorial, let's say `alpha_s_0_2_3` is  small (`4.2e-5`) compared with `alpha_s_0_2_1` (`4.2e-4`). In other words, ozone will adsorb less to the surface as more product is produced. The diffusivity of all components has also been decreased so that we are in a more surface process-limited regime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build and run the simulation\n",
    "\n",
    "# import the Simulate class\n",
    "from multilayerpy.simulate import Simulate\n",
    "\n",
    "# import the Parameter class\n",
    "from multilayerpy.build import Parameter\n",
    "\n",
    "# make the parameter dictionary\n",
    "param_dict = {'delta_3':Parameter(1e-7),  # cm\n",
    "              'alpha_s_0_2_1':Parameter(4.2e-4),  \n",
    "              'alpha_s_0_2_3':Parameter(4.2e-5),\n",
    "              'delta_2':Parameter(0.4e-7),  # cm\n",
    "              'Db_2':Parameter(1e-12),  # cm2 s-1\n",
    "              'delta_1':Parameter(0.8e-7),  # cm\n",
    "              'Db_1':Parameter(1e-12),  # cm2 s-1\n",
    "              'Db_3':Parameter(1e-12),  # cm2 s-1\n",
    "              'k_1_2':Parameter(1.7e-15),  # cm3 s-1\n",
    "              'H_2':Parameter(4.8e-4),  # mol cm-3 atm-1\n",
    "              'Xgs_2': Parameter(7.0e13),  # cm-3\n",
    "              'Td_2': Parameter(1e-2),  # s\n",
    "              'w_2':Parameter(3.6e4),  # cm s-1\n",
    "              'T':Parameter(298.0),  # K\n",
    "              'k_1_2_surf':Parameter(6.0e-12),  # cm2 s-1\n",
    "             }\n",
    "\n",
    "\n",
    "# make the simulate object with the model and parameter dictionary\n",
    "sim = Simulate(model,param_dict)\n",
    "\n",
    "# define required parameters\n",
    "n_layers = 10\n",
    "rp = 0.2e-4 # radius in cm\n",
    "time_span = [0,40] # in s\n",
    "n_time = 999 # number of timepoints to save to output\n",
    "\n",
    "#spherical V and A\n",
    "# use simulate.make_layers function\n",
    "V, A, layer_thick = simulate.make_layers(mod_type,n_layers,rp)\n",
    "\n",
    "# initial conc. of everything\n",
    "\n",
    "bulk_conc_dict = {'1':1.21e21,'2':0,'3':0} # key=model component number, value=bulk conc, units: cm-3\n",
    "surf_conc_dict = {'1':9.68e13,'2':0,'3':0} # key=model component number, value=surf conc, units: cm-2\n",
    "\n",
    "y0 = simulate.initial_concentrations(mod_type,bulk_conc_dict,surf_conc_dict,n_layers) \n",
    "    \n",
    "# now run the model\n",
    "output = sim.run(n_layers,rp,time_span,n_time,V,A,layer_thick,Y0=y0)\n",
    "\n",
    "%matplotlib inline\n",
    "# plot the model\n",
    "fig = sim.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the bulk concentrations\n",
    "sim.plot_bulk_concs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison with and without composition-dependent surface adsorption\n",
    "\n",
    "The next block of code will compare with and without composition-dependent surface adsorption. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make a an empty list for simulate objects, append to later\n",
    "sim_objs = []\n",
    "\n",
    "# non-composition-dependent adsorption ozone\n",
    "O3 = ModelComponent(2,reaction_scheme,gas=True,name='Ozone')   \n",
    "\n",
    "# composition-dependent adsorption ozone\n",
    "O3_comp_depend_ads = ModelComponent(2,reaction_scheme,gas=True,name='Ozone_comp_dependent',comp_dependent_adsorption=True)\n",
    "\n",
    "# 2 parameter dicts for the two models...\n",
    "\n",
    "# without composition-dependent adsorption\n",
    "param_dict_without = {'delta_3':Parameter(1e-7),  # cm\n",
    "              'alpha_s_0_2':Parameter(4.2e-4),\n",
    "              'delta_2':Parameter(0.4e-7),  # cm\n",
    "              'Db_2':Parameter(1e-12),  # cm2 s-1\n",
    "              'delta_1':Parameter(0.8e-7),  # cm\n",
    "              'Db_1':Parameter(1e-12),  # cm2 s-1\n",
    "              'Db_3':Parameter(1e-12),  # cm2 s-1\n",
    "              'k_1_2':Parameter(1.7e-15),  # cm3 s-1\n",
    "              'H_2':Parameter(4.8e-4),  # mol cm-3 atm-1\n",
    "              'Xgs_2': Parameter(7.0e13),  # cm-3\n",
    "              'Td_2': Parameter(1e-2),  # s\n",
    "              'w_2':Parameter(3.6e4),  # cm s-1\n",
    "              'T':Parameter(298.0),  # K\n",
    "              'k_1_2_surf':Parameter(6.0e-12),  # cm2 s-1\n",
    "             }\n",
    "\n",
    "# with composition-dependent adsorption\n",
    "param_dict_with = {'delta_3':Parameter(1e-7),  # cm\n",
    "              'alpha_s_0_2_1':Parameter(4.2e-4),\n",
    "              'alpha_s_0_2_3':Parameter(4.2e-6),\n",
    "              'delta_2':Parameter(0.4e-7),  # cm\n",
    "              'Db_2':Parameter(1e-12),  # cm2 s-1\n",
    "              'delta_1':Parameter(0.8e-7),  # cm\n",
    "              'Db_1':Parameter(1e-12),  # cm2 s-1\n",
    "              'Db_3':Parameter(1e-12),  # cm2 s-1\n",
    "              'k_1_2':Parameter(1.7e-15),  # cm3 s-1\n",
    "              'H_2':Parameter(4.8e-4),  # mol cm-3 atm-1\n",
    "              'Xgs_2': Parameter(7.0e13),  # cm-3\n",
    "              'Td_2': Parameter(1e-2),  # s\n",
    "              'w_2':Parameter(3.6e4),  # cm s-1\n",
    "              'T':Parameter(298.0),  # K\n",
    "              'k_1_2_surf':Parameter(6.0e-12),  # cm2 s-1\n",
    "             }\n",
    "\n",
    "O3_model_component_objects = [O3,O3_comp_depend_ads]\n",
    "param_dict_list = [param_dict_without,param_dict_with]\n",
    "mod_name_ext = ['without','with']\n",
    "\n",
    "# for each \"version\" of ozone (without/with composition-dependent adsorption),run the model\n",
    "for index, component in enumerate(O3_model_component_objects): \n",
    "    \n",
    "    O3_test = component\n",
    "    model_comps_dict = {'1':OA,\n",
    "                        '2':O3_test,\n",
    "                        '3':prod}\n",
    "    \n",
    "    # diffusion regime and reaction scheme are the same and defined previously in this notebook, so I'll skip those steps\n",
    "    # (reusing building blocks)\n",
    "    \n",
    "     # create the model object for this model iteration\n",
    "    model = ModelBuilder(reaction_scheme,model_comps_dict,diff_regime)\n",
    "    model.build(date=False,name_extention=mod_name_ext[index])\n",
    "    \n",
    "    #print(model.req_params) # uncomment to print required parameters\n",
    "    #print(param_dict[\"alpha_s_0_2\"]) # uncomment to check value of alpha_s_0_2\n",
    "    \n",
    "    # make simulate object and run it, use same params as before\n",
    "    # The extra alpha_S0_2_1 and alpha_S0_2_3 params will be ignored by the non-composition-dependent adsorption simulation\n",
    "    \n",
    "    # make the simulate object with the model and parameter dictionary\n",
    "    sim = Simulate(model,param_dict_list[index])\n",
    "\n",
    "    # define required parameters\n",
    "    n_layers = 10\n",
    "    rp = 0.2e-4 # radius in cm\n",
    "    time_span = [0,40] # in s\n",
    "    n_time = 999 # number of timepoints to save to output\n",
    "    \n",
    "\n",
    "\n",
    "    #spherical V and A\n",
    "    # use simulate.make_layers function\n",
    "    V, A, layer_thick = simulate.make_layers(mod_type,n_layers,rp)\n",
    "\n",
    "    # initial conc. of everything\n",
    "\n",
    "    bulk_conc_dict = {'1':1.21e21,'2':0,'3':0} # key=model component number, value=bulk conc\n",
    "    surf_conc_dict = {'1':9.68e13,'2':0,'3':0} # key=model component number, value=surf conc\n",
    "\n",
    "    y0 = simulate.initial_concentrations(mod_type,bulk_conc_dict,surf_conc_dict,n_layers) \n",
    "\n",
    "    # now run the model\n",
    "    output = sim.run(n_layers,rp,time_span,n_time,V,A,layer_thick,Y0=y0)\n",
    "\n",
    "    # append the simulate object to the list\n",
    "    sim_objs.append(sim)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting without composition-dependent adsoprtion\n",
    "sim = sim_objs[0]\n",
    "fig = sim.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting with composition-dependent adsorption\n",
    "sim = sim_objs[1]\n",
    "fig = sim.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary\n",
    "\n",
    "It is clear that with composition-dependent diffusion the decay of oleic acid at the surface of the particle slows down as more product appears at the surface. This appears to also slow down the total rate of oleic acid loss, suggesting surface adsorption in this system affect the decay of the condensed phase reactant somewhat. \n",
    "\n",
    "Surface-dependent surface adsorption can be extended to more complex systems. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}