Finding Effective Diffusivity and Tortuosity of a Network#
In this example we show how to find the effective diffusivity and tortuosity of a cubic network. A similar procedure can be used on more complicated extracted networks from Porespy. The effective diffusivity is defined as the diffusion coefficient multiplied by the porosity and divided by the tortuosity.
The tortuosity is an attempt to correct for the windy path of matter moving through a porous substance. The classical definition is the actual length divided by the direct length. In reality, there are many actual pathways that matter could move through. This means that the actual length should more precisely be thought of as the average path length.
import openpnm as op
import matplotlib.pyplot as plt
import numpy as np
op.visualization.set_mpl_style()
np.random.seed(10)
%matplotlib inline
np.set_printoptions(precision=5)
Network#
For this example we generate a cubic network.
shape = [10, 10, 1]
spacing = 1e-5
net = op.network.Cubic(shape=shape, spacing=spacing)
Geometry#
Next, we need to add geometry models to the generated network. OpenPNM has collections of geometry models ready to add to the network object. Choosing the right geometry models is important to represent the microstructure of different materials such as Toray090 carbon papers, sand stone, electrospun fibers, etc. For now, we will use a simple collection of geometry models for a spheres_and_cylinders geometry. These geometry models assume pores are spheres and throats are cylinders.
geo = op.models.collections.geometry.spheres_and_cylinders
net.add_model_collection(geo, domain='all')
net.regenerate_models()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[3], line 3
1 geo = op.models.collections.geometry.spheres_and_cylinders
2 net.add_model_collection(geo, domain='all')
----> 3 net.regenerate_models()
File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/openpnm/core/_models.py:476, in ModelsMixin2.regenerate_models(self, propnames, exclude)
474 for item in propnames:
475 try:
--> 476 self.run_model(item)
477 except KeyError as e:
478 msg = (f"{item} was not run since the following property"
479 f" is missing: {e}")
File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/openpnm/core/_models.py:507, in ModelsMixin2.run_model(self, propname, domain)
505 if item.startswith(propname+"@"):
506 _, _, domain = item.partition("@")
--> 507 self.run_model(propname=propname, domain=domain)
508 else: # domain was given explicitly
509 domain = domain.split('.', 1)[-1]
File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/openpnm/core/_models.py:521, in ModelsMixin2.run_model(self, propname, domain)
519 if 'domain' not in inspect.getfullargspec(mod_dict['model']).args:
520 _ = kwargs.pop('domain', None)
--> 521 vals = mod_dict['model'](self, **kwargs)
522 if isinstance(vals, dict): # Handle models that return a dict
523 for k, v in vals.items():
File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/openpnm/models/misc/_neighbor_lookups.py:123, in from_neighbor_pores(target, prop, mode, ignore_nans)
121 network = target.network
122 throats = target.Ts
--> 123 P12 = network.find_connected_pores(throats)
124 pvalues = target[prop][P12]
125 if ignore_nans:
File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/openpnm/network/_network.py:443, in Network.find_connected_pores(self, throats, flatten, mode)
391 r"""
392 Return a list of pores connected to the given list of throats
393
(...)
440
441 """
442 Ts = self._parse_indices(throats)
--> 443 pores = topotools.find_connected_sites(bonds=Ts, network=self,
444 flatten=flatten, logic=mode)
445 return pores
File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/openpnm/topotools/_graphtools.py:37, in find_connected_sites(bonds, **kwargs)
36 def find_connected_sites(bonds, **kwargs):
---> 37 return queries.find_connected_nodes(inds=bonds, **kwargs)
File /opt/hostedtoolcache/Python/3.10.18/x64/lib/python3.10/site-packages/openpnm/_skgraph/queries/_funcs.py:167, in find_connected_nodes(network, inds, flatten, logic)
165 temp = temp.astype(float)
166 temp[inds] = np.nan
--> 167 temp = np.reshape(a=temp, newshape=[len(edges), 2], order='F')
168 neighbors = temp
169 else:
TypeError: reshape() got some positional-only arguments passed as keyword arguments: 'a'
A collection of geometry models should now be added to the network object. Make sure to regenerate_models after adding any new models to the network. To view the models that have been added to the network object print the network as follows. Notice how throat.diameter and pore.volume as well as other shape factor models have been added to the network.
print(net)
══════════════════════════════════════════════════════════════════════════════
net : <openpnm.network.Cubic at 0x1bb3553c630>
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# Properties Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
1 pore.coordination_number 100 / 100
2 pore.coords 100 / 100
3 pore.diameter 100 / 100
4 pore.max_size 100 / 100
5 pore.seed 100 / 100
6 pore.volume 100 / 100
7 throat.conns 180 / 180
8 throat.cross_sectional_area 180 / 180
9 throat.diameter 180 / 180
10 throat.diffusive_size_factors 180 / 180
11 throat.hydraulic_size_factors 180 / 180
12 throat.length 180 / 180
13 throat.lens_volume 180 / 180
14 throat.max_size 180 / 180
15 throat.spacing 180 / 180
16 throat.total_volume 180 / 180
17 throat.volume 180 / 180
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# Labels Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
1 pore.all 100
2 pore.back 10
3 pore.front 10
4 pore.left 10
5 pore.right 10
6 pore.surface 36
7 throat.all 180
8 throat.surface 36
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Phase#
Now we add a phase object to our simulation. A phase object contains thermophysical information about the working fluid in the simulation. For this simulation, we will use air as our working fluid. OpenPNM has other built-in phase objects to choose from.
air = op.phase.Air(network=net)
Physics#
Next we need to add physics models to the phase object. Here we add the basic collection of physics models. This includes a diffusive_conductance model which we will need for our Fickian diffusion simulation. It is easy to remove any physics models you don’t need such as throat.entry_pressure.
phys = op.models.collections.physics.basic
del phys['throat.entry_pressure']
air.add_model_collection(phys)
air.regenerate_models()
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:
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.
inlet = net.pores('left')
outlet = net.pores('right')
C_in, C_out = [10, 5]
fd.set_value_BC(pores=inlet, values=C_in)
fd.set_value_BC(pores=outlet, values=C_out)
Now, it’s time to run the algorithm. This is done by calling the run method attached to the algorithm object.
fd.run();
Visualize the Results#
Now that we know the quantity for which FickianDiffusion was solved, let’s take a look at the results. We use the plot_coordinates and plot_connections available under visualization. Throat data is interpolated from pore data.
pc = fd['pore.concentration']
tc = fd.interpolate_data(propname='throat.concentration')
d = net['pore.diameter']
fig, ax = plt.subplots(figsize=[5, 5])
op.visualization.plot_coordinates(network=net, color_by=pc, size_by=d, markersize=400, ax=ax)
op.visualization.plot_connections(network=net, color_by=tc, linewidth=3, ax=ax)
_ = plt.axis('off')
Calculate the Effective Diffusivity#
We can determine the effective diffusivity of the network by solving Fick’s law as shown below.
To determine the molar flow rate we can use the rate method attached to the algorithm object. The molar flow rate going through the boundary pores at the inlet is calculated. This gives the molar flow rate in units of moles per second.
rate_inlet = fd.rate(pores=inlet)[0]
print(f'Molar flow rate: {rate_inlet:.5e} mol/s')
Molar flow rate: 3.88846e-11 mol/s
Now let’s calculate the effective diffusivity
A = (shape[1] * shape[2])*(spacing**2)
L = shape[0]*spacing
D_eff = rate_inlet * L / (A * (C_in - C_out))
print("{0:.6E}".format(D_eff))
7.776925E-07
Calculate the Tortuosity#
Tortuosity is most easily determined as a fitting factor to transport data. We can use the effective diffusivity measured from the Fickian diffusion simulation and the known diffusivity of air to solve for the tortuosity.
But first the porosity of the network must be calculated. Here is how we can calculate the porosity. For more information on how to measure the porosity please refer to the porosity notebook.
V_p = net['pore.volume'].sum()
V_t = net['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.429901E-02
D_AB = air['pore.diffusivity'][0]
tau = e * D_AB / D_eff
print('The tortuosity is:', "{0:.6E}".format(tau))
The tortuosity is: 2.270571E+00