Fickian Diffusion with Reaction

OpenPNM is capable of simulating chemical reactions in pores by adding source and sink terms. This example shows how to add source and sink terms to a steady state fickian diffusion simulation.

import warnings
import scipy as sp
import numpy as np
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
import matplotlib.pyplot as plt
np.set_printoptions(precision=5)
np.random.seed(10)
%matplotlib inline

Create Network and Other Objects

Start by creating the network, geometry, phase and physics objects as usual:

pn = op.network.Cubic(shape=[40, 40], spacing=1e-4)
geo = op.geometry.SpheresAndCylinders(network=pn, pores=pn.Ps, throats=pn.Ts)
gas = op.phases.Air(network=pn)
phys = op.physics.Standard(network=pn, phase=gas, geometry=geo)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [2], in <cell line: 3>()
      1 pn = op.network.Cubic(shape=[40, 40], spacing=1e-4)
      2 geo = op.geometry.SpheresAndCylinders(network=pn, pores=pn.Ps, throats=pn.Ts)
----> 3 gas = op.phases.Air(network=pn)
      4 phys = op.physics.Standard(network=pn, phase=gas, geometry=geo)

AttributeError: module 'openpnm' has no attribute 'phases'

Now add the source and sink models to the physics object. In this case we’ll think of the as chemical reactions. We’ll add one source term and one sink term, meaning one negative reaction rate and one positive reaction rate

gas['pore.concentration'] = 0
phys['pore.sinkA'] = -1e-10
phys['pore.sinkb'] = 1
phys.add_model(propname='pore.sink', model=op.models.physics.generic_source_term.power_law,
               A1='pore.sinkA', A2='pore.sinkb', X='pore.concentration')
phys['pore.srcA'] = 1e-11
phys['pore.srcb'] = 1
phys.add_model(propname='pore.source', model=op.models.physics.generic_source_term.power_law,
               A1='pore.srcA', A2='pore.srcb', X='pore.concentration')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [3], in <cell line: 1>()
----> 1 gas['pore.concentration'] = 0
      2 phys['pore.sinkA'] = -1e-10
      3 phys['pore.sinkb'] = 1

NameError: name 'gas' is not defined

Setup Fickian Diffusion with Sink Terms

Now we setup a FickianDiffusion algorithm, with concentration boundary conditions on two side, and apply the sink term to 3 pores:

rx = op.algorithms.FickianDiffusion(network=pn, phase=gas)
rx.set_source(propname='pore.sink', pores=[420, 820, 1220])
rx.set_value_BC(values=1, pores=pn.pores('front'))
rx.set_value_BC(values=1, pores=pn.pores('back'))
rx.run()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 rx = op.algorithms.FickianDiffusion(network=pn, phase=gas)
      2 rx.set_source(propname='pore.sink', pores=[420, 820, 1220])
      3 rx.set_value_BC(values=1, pores=pn.pores('front'))

NameError: name 'gas' is not defined

Plot Distributions

We can use the plot_connections and plot_coordinates to get a quick view of the pore network, with colors and sizes scaled appropriately.

pn['pore.concentration'] = rx['pore.concentration']
fig, ax = plt.subplots(figsize=[15, 15])
op.topotools.plot_connections(network=pn, color='k', alpha=.2, ax=ax)
op.topotools.plot_coordinates(network=pn, color_by=pn['pore.concentration'],
                              size_by=pn['pore.diameter'], cmap='plasma',
                              ax=ax, markersize=100)
_ = plt.axis('off')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [5], in <cell line: 1>()
----> 1 pn['pore.concentration'] = rx['pore.concentration']
      2 fig, ax = plt.subplots(figsize=[15, 15])
      3 op.topotools.plot_connections(network=pn, color='k', alpha=.2, ax=ax)

NameError: name 'rx' is not defined

Plot Distributions as Heatmaps

Because the network is a 2D cubic, it is convenient to visualize it as an image, so we reshape the ‘pore.concentration’ array that is produced by the FickianDiffusion algorithm upon running, and turn it into a colormap representing concentration in each pore.

pn['pore.concentration'] = rx['pore.concentration']
fig, ax = plt.subplots(figsize=[15, 15])
op.topotools.plot_connections(network=pn, color='k', alpha=0, ax=ax)
op.topotools.plot_coordinates(network=pn, color_by=pn['pore.concentration'], cmap='plasma',
                              ax=ax, markersize=420, marker='s')
_ = plt.axis('off')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [6], in <cell line: 1>()
----> 1 pn['pore.concentration'] = rx['pore.concentration']
      2 fig, ax = plt.subplots(figsize=[15, 15])
      3 op.topotools.plot_connections(network=pn, color='k', alpha=0, ax=ax)

NameError: name 'rx' is not defined

Setup Fickian Diffusion with Source Terms

Similarly, for the source term:

rx = op.algorithms.FickianDiffusion(network=pn, phase=gas)
rx.set_source(propname='pore.source', pores=[420, 820, 1220])
rx.set_value_BC(values=1, pores=pn.pores('front'))
rx.set_value_BC(values=1, pores=pn.pores('back'))
rx.run()

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [7], in <cell line: 1>()
----> 1 rx = op.algorithms.FickianDiffusion(network=pn, phase=gas)
      2 rx.set_source(propname='pore.source', pores=[420, 820, 1220])
      3 rx.set_value_BC(values=1, pores=pn.pores('front'))

NameError: name 'gas' is not defined

And plotting the result as a color map:

pn['pore.concentration'] = rx['pore.concentration'] - 1
fig, ax = plt.subplots(figsize=[15, 15])
op.topotools.plot_connections(network=pn, color='k', alpha=0, ax=ax)
op.topotools.plot_coordinates(network=pn, color_by=pn['pore.concentration'], cmap='plasma',
                              ax=ax, markersize=420, marker='s')
_ = plt.axis('off')

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 pn['pore.concentration'] = rx['pore.concentration'] - 1
      2 fig, ax = plt.subplots(figsize=[15, 15])
      3 op.topotools.plot_connections(network=pn, color='k', alpha=0, ax=ax)

NameError: name 'rx' is not defined