Fickian Diffusion and Tortuosity#

In this example, we will learn how to perform Fickian diffusion on a Cubic network. The algorithm works fine with every other network type, but for now we want to keep it simple. Refere to Tutorials > Network for more details on different network types.

[1]:
import numpy as np
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(10)
ws = op.Workspace()
ws.settings["loglevel"] = 40
np.set_printoptions(precision=5)

Generating network#

First, we need to generate a Cubic network. For now, we stick to a 2d network, but you might as well try it in 3d!

[2]:
shape = [1, 10, 10]
spacing = 1e-5
net = op.network.Cubic(shape=shape, spacing=spacing)

Adding geometry#

Next, we need to add a geometry to the generated network. A geometry contains information about size of the pores/throats in a network. OpenPNM has tons of prebuilt geometries that represent the microstructure of different materials such as Toray090 carbon papers, sand stone, electrospun fibers, etc. For now, we’ll stick to a simple geometry called SpheresAndCylinders that assigns random values to pore/throat diameters.

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

Adding phase#

Next, we need to add a phase to our simulation. A phase object(s) contain(s) thermophysical information about the working fluid(s) in the simulation. OpenPNM has tons of prebuilt phases as well! For this simulation, we use air as our working fluid.

[4]:
air = op.phase.Air(network=net)

Adding physics#

Finally, we need to add a physics. A physics object contains information about the working fluid in the simulation that depend on the geometry of the network. A good example is diffusive conductance, which not only depends on the thermophysical properties of the working fluid, but also depends on the geometry of pores/throats. OpenPNM includes a pre-defined physics class called Standard which as the name suggests contains all the standard pore-scale models to get you going:

[5]:
phys_air = op.physics.Standard(network=net, phase=air, geometry=geom)

Performing Fickian diffusion#

Now that everything’s set up, it’s time to perform our Fickian diffusion simulation. For this purpose, we need to add the FickianDiffusion algorithm to our simulation. Here’s how we do it:

[6]:
fd = op.algorithms.FickianDiffusion(network=net, phase=air)

Note that network and phase are required parameters for pretty much every algorithm we add, since we need to specify on which network and for which phase we want to run the algorithm.

Adding boundary conditions#

Next, we need to add some boundary conditions to the simulation. By default, OpenPNM assumes zero flux for the boundary pores.

[7]:
inlet  = net.pores('front')
outlet = net.pores('back')
C_in = 1.0
C_out = 0.0
fd.set_value_BC(pores=inlet, values=C_in)
fd.set_value_BC(pores=outlet, values=C_out)

set_value_BC applies the so-called “Dirichlet” boundary condition to the specified pores. Note that unless you want to apply a single value to all of the specified pores (like we just did), you must pass a list (or ndarray) as the values parameter.

Running the algorithm#

Now, it’s time to run the algorithm. This is done by calling the run method attached to the algorithm object.

[8]:
fd.run()
[8]:
{'pore.concentration': SteadyStateSolution([1.     , 1.     , 1.     , 1.     , 1.     , 1.     ,
                      1.     , 1.     , 1.     , 1.     , 0.91247, 0.90179,
                      0.90063, 0.91127, 0.91092, 0.89035, 0.87522, 0.87787,
                      0.85195, 0.84409, 0.8086 , 0.80099, 0.77881, 0.82206,
                      0.83878, 0.81892, 0.8155 , 0.77651, 0.80579, 0.79607,
                      0.71218, 0.69453, 0.70386, 0.70642, 0.70437, 0.68794,
                      0.687  , 0.66227, 0.69233, 0.73643, 0.63349, 0.63437,
                      0.59717, 0.60263, 0.57961, 0.56621, 0.57162, 0.59177,
                      0.55168, 0.57888, 0.50755, 0.4831 , 0.44376, 0.47502,
                      0.51439, 0.51734, 0.50933, 0.49063, 0.45289, 0.4257 ,
                      0.39291, 0.35487, 0.35859, 0.33091, 0.37506, 0.43326,
                      0.3751 , 0.34388, 0.32536, 0.30339, 0.24265, 0.25011,
                      0.25119, 0.25372, 0.22365, 0.24942, 0.21214, 0.21747,
                      0.20779, 0.20799, 0.08587, 0.12816, 0.14487, 0.1239 ,
                      0.11533, 0.11314, 0.10743, 0.11603, 0.1247 , 0.1043 ,
                      0.     , 0.     , 0.     , 0.     , 0.     , 0.     ,
                      0.     , 0.     , 0.     , 0.     ])}

Post processing#

When an algorithm is successfully run, the results are attached to the same object. To access the results, you need to know the quantity for which the algorithm was solving. For instance, FickianDiffusion solves for the quantity pore.concentration, which is somewhat intuitive. However, if you ever forget it, or wanted to manually check the quantity, you can take a look at the algorithm settings:

[9]:
print(fd.settings)
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Settings                            Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
name                                fick_01
prefix                              fick
uuid                                18c98e7a-f442-4a4a-af0e-8e34279a17a2
cache                               True
conductance                         throat.diffusive_conductance
f_rtol                              1e-06
newton_maxiter                      5000
phase                               phase_01
quantity                            pore.concentration
relaxation_quantity                 1.0
sources                             []
variable_props                      TypedSet()
x_rtol                              1e-06
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Visualizing#

Now that we know the quantity for which FickianDiffusion was solved, let’s take a look at the results:

[10]:
c = fd['pore.concentration']
r = fd.rate(throats=net.Ts, mode='single')
d = net['pore.diameter']
fig, ax = plt.subplots(figsize=[4, 4])
op.topotools.plot_coordinates(network=net, color_by=c, size_by=d, markersize=400, ax=ax)
op.topotools.plot_connections(network=net, color_by=r, linewidth=3, ax=ax)
_ = plt.axis('off')
../../../_images/examples_simulations_steady_state_fickian_diffusion_and_tortuosity_22_0.svg

Calculating flux#

You might as well be interested in calculating the mass flux from a boundary! This is easily done in OpenPNM via calling the rate method attached to the algorithm. Let’s see how it works:

[11]:
rate_inlet = fd.rate(pores=inlet)[0]
print(f'Mass flow rate from inlet: {rate_inlet:.5e} mol/s')
Mass flow rate from inlet: 7.67614e-12 mol/s

We can determine the effective diffusivity of the network by solving Fick’s law:

\[D_{eff} = \frac{N_A L}{ A \Delta C}\]
[12]:
A = (shape[0] * shape[1])*(spacing**2)
L = shape[2]*spacing
D_eff = rate_inlet * L / (A * (C_in - C_out))
print("{0:.6E}".format(D_eff))
7.676140E-07

And the formation factor can be found since the diffusion coefficient of open air is known:

\[F = \frac{D_{AB}}{D_{eff}}\]
[13]:
D_AB = air['pore.diffusivity'][0]
F = D_AB / D_eff
print('The formation factor is: ', "{0:.6E}".format(F))
The formation factor is:  2.693473E+01

The tortuosity is defined as follows:

\[\frac{D_{eff}}{D_{AB}} = \frac{\varepsilon}{\tau} \rightarrow \tau = \varepsilon \frac{ D_{AB}}{D_{eff}}\]

Note that finding the tortuosity requires knowing the porosity, which is annoyingly difficult to calculate accurately, so here we will just gloss over the details.

[14]:
V_p = geom['pore.volume'].sum()
V_t = geom['throat.volume'].sum()
V_bulk = np.prod(shape)*(spacing**3)
e = (V_p + V_t) / V_bulk
print('The porosity is: ', "{0:.6E}".format(e))
The porosity is:  8.546959E-02
[15]:
tau = e * D_AB / D_eff
print('The tortuosity is:', "{0:.6E}".format(tau))
The tortuosity is: 2.302101E+00