Transient Fickian Diffusion#

The package OpenPNM allows for the simulation of many transport phenomena in porous media such as Stokes flow, Fickian diffusion, advection-diffusion, transport of charged species, etc. Transient and steady-state simulations are both supported. An example of a transient Fickian diffusion simulation through a Cubic pore network is shown here.

First, OpenPNM is imported.

[1]:
import numpy as np
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
np.random.seed(10)
%matplotlib inline
np.set_printoptions(precision=5)

Define new workspace and project#

[2]:
ws = op.Workspace()
ws.settings["loglevel"] = 40
proj = ws.new_project()

Generate a pore network#

An arbitrary Cubic 3D pore network is generated consisting of a layer of \(29\times13\) pores with a constant pore to pore centers spacing of \({10}^{-4}{m}\).

[3]:
shape = [13, 29, 1]
net = op.network.Cubic(shape=shape, spacing=1e-4, project=proj)

Create a geometry#

Here, a geometry, corresponding to the created network, is created. The geometry contains information about the size of pores and throats in the network such as length and diameter, etc. OpenPNM has many prebuilt geometries that represent the microstructure of different materials such as Toray090 carbon papers, sand stone, electrospun fibers, etc. In this example, a simple geometry known as SpheresAndCylinders that assigns random diameter values to pores throats, with certain constraints, is used.

[4]:
geo = op.geometry.SpheresAndCylinders(network=net, pores=net.Ps, throats=net.Ts)

Add a phase#

Then, a phase (water in this example) is added to the simulation and assigned to the network. The phase contains the physical properties of the fluid considered in the simulation such as the viscosity, etc. Many predefined phases as available on OpenPNM.

[5]:
phase = op.phases.Water(network=net)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [5], in <cell line: 1>()
----> 1 phase = op.phases.Water(network=net)

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

Add a physics#

Next, a physics object is defined. The physics object stores information about the different physical models used in the simulation and is assigned to specific network, geometry and phase objects. This ensures that the different physical models will only have access to information about the network, geometry and phase objects to which they are assigned. In fact, models (such as Stokes flow or Fickian diffusion) require information about the network (such as the connectivity between pores), the geometry (such as the pores and throats diameters), and the phase (such as the diffusivity coefficient).

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

NameError: name 'phase' is not defined

The diffusivity coefficient of the considered chemical species in water is also defined.

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

NameError: name 'phase' is not defined

Defining a new model#

The physical model, consisting of Fickian diffusion, is defined and attached to the physics object previously defined.

[8]:
mod = op.models.physics.diffusive_conductance.ordinary_diffusion
phys.add_model(propname='throat.diffusive_conductance', model=mod, regen_mode='normal')
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [8], in <cell line: 2>()
      1 mod = op.models.physics.diffusive_conductance.ordinary_diffusion
----> 2 phys.add_model(propname='throat.diffusive_conductance', model=mod, regen_mode='normal')

NameError: name 'phys' is not defined

Define a transient Fickian diffusion algorithm#

Here, an algorithm for the simulation of transient Fickian diffusion is defined. It is assigned to the network and phase of interest to be able to retrieve all the information needed to build systems of linear equations.

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

NameError: name 'phase' is not defined

Add boundary conditions#

Next, Dirichlet boundary conditions are added over the back and front boundaries of the network.

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

NameError: name 'fd' is not defined

Define initial conditions#

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.

Setup the transient algorithm settings#

The settings of the transient algorithm are updated here. When calling alg.run, you can pass the following arguments: - x0: initial conditions - tspan: integration time span - saveat: the interval at which the solution is to be stored

[11]:
x0 = 0.2
tspan = (0, 100)
saveat = 5

Run the algorithm#

The algorithm is run here.

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

NameError: name 'fd' is not defined

Post process and export the results#

Once the simulation is successfully performed. The solution at every time steps is stored within the algorithm object. The algorithm’s stored information is printed here.

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

NameError: name 'fd' is not defined

Note that the solutions at every exported time step contain the @ character followed by the time value. Here the solution is exported after each \(5s\) in addition to the final time step which is not a multiple of \(5\) in this example.

To print the solution at \(t=10s\)

[15]:
soln(10)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [15], in <cell line: 1>()
----> 1 soln(10)

NameError: name 'soln' is not defined

The solution is here stored in the phase before export.

[16]:
phase.update(fd.results())
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [16], in <cell line: 1>()
----> 1 phase.update(fd.results())

NameError: name 'phase' is not defined

Visialization using Matplotlib#

One can perform post processing and visualization using the exported files on an external software such as Paraview. Additionally, the Pyhton library Matplotlib can be used as shown here to plot the concentration color map at steady-state.

[17]:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

c = fd.x.reshape(shape)
fig, ax = plt.subplots(figsize=(6, 6))
im = ax.imshow(c[:,:,0])
ax.set_title('concentration (mol/m$^3$)')
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 [17], in <cell line: 4>()
      1 import matplotlib.pyplot as plt
      2 from mpl_toolkits.axes_grid1 import make_axes_locatable
----> 4 c = fd.x.reshape(shape)
      5 fig, ax = plt.subplots(figsize=(6, 6))
      6 im = ax.imshow(c[:,:,0])

NameError: name 'fd' is not defined