Random PNM with Delaunay Network and Voronoi Geometry using a new material class VoronoiFibers

Getting Started

In this tutorial, we will demonstrate concepts of random pore network modelling outlined in J.T. Gostick's 2013 paper and subsequently used in some of our other papers investigating fuel cell gas diffusion layer compression Tranter et. al 2016 and capillary hysteresis Tranter et. al 2017. Pores are randomly dispersed inside the domain and connections are defined by the nearest neighbour Delaunay triangulation algorithm. The Voronoi diagram is the compliment to the Delaunay trinagulation and is used to replicate the fibres in a porous domain such as a fuel cell GDL. Each pore is surrounded by a cage of fibres which are created from the intersections of the equidistant planes between neighbouring pores. Where the planes intersect, a Voronoi vertex is created.. A throat is defined by a set of vertices which are shared by two neighbouring pores. The throat vertices will be coplanar and the throat normal is the vector from one pore's coordinate to the other. N.B this coordinate is not the pore centroid but can be close if the shape of the cage is near spherical. The vertices are used by the Voronoi geometry model which creates a 3D image of the fibres using a supplied fibre radius. Image analysis is then performed to extract pore and throat sizes by using the convex hull of the pore and throat vertices.

Setting up Network and Geometry

We first import the OpenPNM code, including a utility submodule which has useful functions for working with the Voronoi classes:


In [1]:
import openpnm as op
from openpnm.utils import vertexops as vo
import scipy as sp
import matplotlib.pyplot as plt


C:\Users\Tom\Anaconda3\lib\site-packages\h5py\__init__.py:34: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
  from ._conv import register_converters as _register_converters

Next we instantiate the VoronoiFibers classs with a few parameters: num_points defines the number of pores in the bulk network but some may be removed if the fiber_rad (radius of the fibers in meters) parameter is large compared with the throat diameters as this will lead to occluded throats and possibly isolated pores. The resolution parameter sets the size of the voxels in meters used in image generation from which many of the geometry properties derive. The shape parameter sets the size of the bounding box at which the boundary edges are located and this is positioned with lowest point at the origin. Finally the name parameter is used to subsequently prefix the name of the additional geometry objects that are created in the process.

In this example we set the random seed in the scipy package to make sure that random points are repeatedly generated with the same values for illustration. This is not necessary, otherwise.


In [2]:
scale = 1e-4
wrk = op.core.Workspace()
wrk.loglevel = 50
sp.random.seed(1)
net = op.materials.VoronoiFibers(num_points=100,
                                 fiber_rad=5e-6,
                                 resolution=1e-6,
                                 shape=[scale, scale, scale],
                                 name='test')

We are returned a handle to the network and we can inspect the properties.


In [3]:
print(net)


――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
openpnm.materials.VoronoiFibers: 	test
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
#     Properties                                    Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2     pore.coords                                     778 / 778  
4     throat.conns                                   4077 / 4077 
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
#     Labels                                        Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
1     pore.all                                      778       
2     pore.delaunay                                 226       
3     pore.internal                                 445       
4     pore.surface                                  333       
5     pore.test_del                                 226       
6     pore.test_vor                                 552       
7     pore.voronoi                                  552       
8     throat.all                                    4077      
9     throat.delaunay                               371       
10    throat.interconnect                           2606      
11    throat.internal                               3083      
12    throat.surface                                994       
13    throat.test_del                               371       
14    throat.test_vor                               1100      
15    throat.voronoi                                1100      
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

To get access to the geometry objects we must invoke the project object and retrieve them separately


In [4]:
prj = net.project
print(prj)
del_geom = prj.geometries()['test_del']
vor_geom = prj.geometries()['test_vor']


――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
 Object Name     Object ID                                                        
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
 test            <openpnm.materials.VoronoiFibers object at 0x175bea3b8e0>        
 test_del        <openpnm.materials.VoronoiFibers object at 0x175bea3b518>        
 test_vor        <openpnm.materials.VoronoiFibers object at 0x175a76b9e08>        
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

As you can see we have more than one geometry. This is because the VoronoiFibers class has used the DelaunayVoronoiDual class (covered in another example) which contains pores and throats for the void space in-between fibers assigned to del_geom and also pores and throats for the fibers themselves which can be used for coupled simulations such as heat transfer assigned to vor_geom. We can at this point delete the pores and throats associated with fibers without impacting the traditional pore network as they are separate and connected by the interconnect throats.


In [5]:
from openpnm import topotools as tt
tt.trim(network=net, pores=net.pores('voronoi'))
print(net)


――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
openpnm.materials.VoronoiFibers: 	test
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
#     Properties                                    Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2     pore.coords                                     226 / 226  
4     throat.conns                                    371 / 371  
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
#     Labels                                        Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
1     pore.all                                      226       
2     pore.delaunay                                 226       
3     pore.internal                                 100       
4     pore.surface                                  126       
5     pore.test_del                                 226       
6     pore.test_vor                                 0         
7     pore.voronoi                                  0         
8     throat.all                                    371       
9     throat.delaunay                               371       
10    throat.interconnect                           0         
11    throat.internal                               371       
12    throat.surface                                0         
13    throat.test_del                               371       
14    throat.test_vor                               0         
15    throat.voronoi                                0         
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Note that trimming the voronoi pores has also trimmed the interconnect and voronoi throats as one or both of the pores that the throats connect have been removed, respectively. We can now inspect the material class visually using some functions from vertexops


In [6]:
vo.plot_pore(del_geom, pores=del_geom.pores())


Out[6]:

In [7]:
vo.plot_throat(del_geom, throats=del_geom.throats()[:4])


Out[7]:

It is worth mentioning a little about the boundaries at this point. Behind the scenes all the pores in the network were temporarily reflected about the planes confining the pore space during the Delaunay tessellation. This has the effect of creating throats on the outer confines that are aligned with the boundary planes. The boundary pores are labelled 'surface' and can be further identified by the shared coordinate which will either be zero or the values of the domain shape used to initialize the class. As is custom, with openpnm simulations the fictitious boundary pores are assigned zero volume. The pore-size-distribution can be plotted but will be skewed by including the boundary pores:


In [8]:
plt.close('all')
plt.figure()
del_geom.show_hist(props=['pore.diameter'])
plt.show()


However, the throats connecting the boundary pores to their neigbors have real non-zero dimensions so the throat-size-distribution looks realistic


In [9]:
plt.figure()
del_geom.show_hist(props=['throat.diameter'])
plt.show()


In previous versions of openpnm (1.x) the Voronoi geometry could be generated with the option to use image analysis pore calculations of the pore sizes and other properties, or functions based on using idealized shapes formed from the Voronoi vertices could also be used which is faster but less accurate, especially when dealing with highly anisotropic media. The materials class does not give this option.

Generating the voxel image is a memory intensive process relying on many image analysis routines and it is recommended that a smaller network is tested first on your machine whilst monitoring your system performance to gauge whether larger networks are possible.

The image based approach allows for closer inspection of realistic fibrous structures and a few useful functions are provided for analysis and inspection.


In [10]:
del_geom.plot_porosity_profile()


Out[10]:

In [11]:
del_geom.plot_fiber_slice(plane=[0, 0.5, 0])

Two images are actually stored as private attributes on the delaunay geometry object: _fiber_image and _hull_image. The fiber image is a binary image with 1 representing the fiber and 0 the void space. The hull image represents the convex hulls surrounding each pore and is labelled with the original pore index (which may change after trimming). These images can be exported as tiff stacks and visualized in paraview, as well as being used to populate liquid filled porous networks. An example of such an image is shown below and is taken from Tranter et. al 2016