# Transient Fickian Diffusion with Reaction#

OpenPNM supports adding reaction terms to both steady state and transient simulations. OpenPNM already includes many different source term models that can be added to simulate a reaction. In this example, we show how to add a powerlaw source term model to a transient fickian diffusion simulation.

Start by importing openpnm

[1]:

import openpnm as op
%config InlineBackend.figure_formats = ['svg']
import matplotlib.pyplot as plt
import numpy as np
ws = op.Workspace()
proj = ws.new_project()
np.random.seed(10)


## Define network, geometry, and phase objects#

[2]:

shape = [25, 25, 1]
net = op.network.Cubic(shape=shape, spacing=1e-2)
geo = op.geometry.SpheresAndCylinders(network=net, pores=net.Ps, throats=net.Ts)
phase = op.phases.Water(network=net)

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

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


You can adjust the diffusion coefficient by calling the pore.diffusivity dictionary key on the phase object and then setting it’s value.

[3]:

phase['pore.diffusivity'] = 2e-07

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [3], in <cell line: 1>()
----> 1 phase['pore.diffusivity'] = 2e-07

NameError: name 'phase' is not defined


## Define physics object#

Here, we will use Standard physics which already includes many standard models. You could also use Generic physics and add neccessary models like throat.diffusive_conductance after using phys.add_model() method.

[4]:

phys = op.physics.Standard(network=net, phase=phase, geometry=geo)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 phys = op.physics.Standard(network=net, phase=phase, geometry=geo)

NameError: name 'phase' is not defined


## Add reaction model#

Add reaction model to physics object. For this example we use the powerlaw reaction model which is of the form:

\begin{equation*} \ A_1X^{A_2} \end{equation*}

[5]:

phase['pore.concentration'] = 0
phys['pore.rxnA'] = -1e-10
phys['pore.rxnb'] = 1
phys.add_model(propname='pore.reaction', model=op.models.physics.generic_source_term.power_law,
A1='pore.rxnA', A2='pore.rxnb', X='pore.concentration')

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

NameError: name 'phase' is not defined


## Define transient fickian diffusion object#

[6]:

tfd = op.algorithms.TransientFickianDiffusion(network=net, phase=phase)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [6], in <cell line: 1>()
----> 1 tfd = op.algorithms.TransientFickianDiffusion(network=net, phase=phase)

NameError: name 'phase' is not defined


## Setup the transient algorithm settings#

To do this you can call the tfd.setup() method and use it to set settings such as transient solver scheme, final time, time step, tolerance, etc. The cranknicolson scheme used here is the most accurate but slowest. Other time schemes are implicit which is faster but less accurate and steady which gives the steady state solution.

[7]:

x0 = 0
tspan = (0, 2000)
saveat = 100


## Set value boundary conditions#

In this example we set the concentraton of the front pores to 1 and the concentration of the back pores to 0.

[8]:

tfd.set_value_BC(pores=net.pores('front'), values=1)
tfd.set_value_BC(pores=net.pores('back'), values=0)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 tfd.set_value_BC(pores=net.pores('front'), values=1)
2 tfd.set_value_BC(pores=net.pores('back'), values=0)

NameError: name 'tfd' is not defined


## Set initial condition#

Initial conditions must be specified when alg.run is called as alg.run(x0=x0), where x0 could either be a scalar (in which case it’ll be broadcasted to all pores), or an array.

## Set source term#

In this example, we apply the source term to pores 212 and 412.

[9]:

tfd.set_source(propname='pore.reaction', pores=[212, 412])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [9], in <cell line: 1>()
----> 1 tfd.set_source(propname='pore.reaction', pores=[212, 412])

NameError: name 'tfd' is not defined


## Run the simulation#

[10]:

tfd.run(x0=x0, tspan=tspan, saveat=saveat)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [10], in <cell line: 1>()
----> 1 tfd.run(x0=x0, tspan=tspan, saveat=saveat)

NameError: name 'tfd' is not defined


## Visualize results#

Ater simulation runs we can visualize results using a colour plot from matplotlib.pyplot. Here we visualize the results at the final time t_final by using the pore.concentration key on the algorithm object.

[11]:

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

c = tfd['pore.concentration'].reshape(shape)
fig, ax = plt.subplots(figsize=(5, 5))
im = ax.imshow(c[:,:,0])
ax.set_title(f'concentration (mol/m$^3$) @ t = {tspan[1]}')
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="4%", pad=0.1)
plt.colorbar(im, cax=cax);

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [11], in <cell line: 4>()
1 import matplotlib.pyplot as plt
2 from mpl_toolkits.axes_grid1 import make_axes_locatable
----> 4 c = tfd['pore.concentration'].reshape(shape)
5 fig, ax = plt.subplots(figsize=(5, 5))
6 im = ax.imshow(c[:,:,0])

NameError: name 'tfd' is not defined


If we print the TransientFickianDiffusion object we can see a list of the object’s properties. Notice how the concentration is recorded here for each output concentration, t_output.

[12]:

print(tfd)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [12], in <cell line: 1>()
----> 1 print(tfd)

NameError: name 'tfd' is not defined


We can visualize intermediate concentration profiles using a colour plot but use pore.concentration@100 (or similar) as the dictionary key on the algorithm object.

[13]:

c = tfd.soln(100).reshape(shape)
fig, ax = plt.subplots(figsize=(5, 5))
im = ax.imshow(c[:,:,0])
ax.set_title(f'concentration (mol/m$^3$) @ t = 100')
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="4%", pad=0.1)
plt.colorbar(im, cax=cax);

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [13], in <cell line: 1>()
----> 1 c = tfd.soln(100).reshape(shape)
2 fig, ax = plt.subplots(figsize=(5, 5))
3 im = ax.imshow(c[:,:,0])

NameError: name 'tfd' is not defined