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

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

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.

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.

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*}

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

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.

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.

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.

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

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.

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.

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.

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