How to create a pore network with a given porosity and pore throat distribution?

[Camelot18]

Dear Prof. Gostick,

Is it possible to generate a pore network with a specific porosity and pore throat size distribution in OpenPNM? It seems that we can only adjust the pore throat distribution function and obtain the porosity and the corresponding histogram, can we do the reverse? Or do we have to create such pore networks in Porespy by creating images and using the snow algorithm?

Also, how can I add more pores with known coordinates and diameters in the existing pore network?

Hope to hear from you soon.

Best,
Peter

[jgostick]

Hi @Camelot18
There is not function for this, but it's not too hard to do. You create a 10x10x10 network, play with the psd parameters, check the properites like porosity, pc curve, permeability, etc, then repeat. There is a bit too much going on here to make it a built in function. This example gives a partial tutorial.

As fo your second question, you'll want to use the op.topotools.extend function. It allows you to add pores and throats. Hopefully the docstring is sufficient to get you started.

[Camelot18]

I tried many times to extract the network via snow (used different values for nbins and r_min), but it seems that the throats volume is always greater than pore volume. This is the reverse to the result from openpnm. How can I get more pore volume via snow to create a resonable REV?

[jgostick]

Hi @Camelot18, I just looked at the code for snow2 to refresh my memory, and I can tell you that snow2 does not return throat volume values. These are probably being added by the pore-scale models in openpnm. In the transition to version 3 of openpnm we trimmed out a lot of code, to make the maintenance more managable. The new preferred approach is to not apply any geometry pore scale models. This way you just need to tell openpnm which of the snow extracted arrays to use for each property, like pn['pore.diameter'] = pn['pore.inscribed_diameter']. We made this process a bit more manual so that we could stop making assumptions for users, since some might want do use the equivalent diameter or whatever.

I will add a todo that we should make an OpenPNM example on working with the snow data.

[Camelot18]

Yes, snow2 doesn't return those results, and I used pore-scale models to add them.

net.add_model(propname='throat.length',
             model=op.models.geometry.throat_length.hybrid_cones_and_cylinders)
conns = net['throat.conns']

model=op.models.geometry.pore_volume.sphere
net.add_model(propname='pore.volume',
             model=model,
             regen_mode='normal')
model=op.models.geometry.throat_volume.cylinder
net.add_model(propname='throat.volume',
             model=model,
             regen_mode='normal')

The op.models.geometry.throat_length.spheres_and_cylinders doesn't work and I have to try .hybrid_cones_and_cylinders instead. Is there anything wrong with this method? How can I get a network with more pores rather than throats?

[Camelot18]

Another question about the vtk export, I used to write the vtk files by op.io.VTK.export_data(network, filename) but now it becomes op.io.project_to_vtk(project, filename). How to get a project for this since I only have network after snow2?

[Camelot18]

The vtk issue has been fixed.

[Camelot18]

Dear Professor,

Thank you for the kind reply.
Is there any difference between openpnm and porespy in creating a porous network with a given porosity and distribution? In porespy, the snow algorithm returns parameters like pore.inscribed_diameter and throat.inscribed_diameter, should I just equal them with pore.diameter and throat.diameter just like this example?

As for the op.topotools.extend function, can I add my own label name instead of pores? Otherwise, I'm afraid that they are messed up with the previous pores. Actually, I need to obtain this csv form where the coordinates and radius are functions of those of pores in the network I created.

Looking forward to your reply.

Best,
Peter

[jgostick]

Hi @Camelot18
Firstly, when I answered this question I thought it was on the OpenPNM repo, so sorry for the confusion. There is a lot of interaction between porespy and openpnm. Porespy is all about image analysis, including network extraction using snow, OpenPNM has network generator functions, based on lattices and also Deluanay/Voronoi tessellations. You could try to make an image in porespy that looks the way you want, or you could use OpenPNM to just generate something that is statistically equivalent.

RE: Labels, I think you can send a label for the new pores/throats, but if not then it's also easy to add a label like 'pore.original' prior to calling extend so you'll know the new pores/throats by their lack of the 'pore.original' label.

[jgostick]

I should also add that I am actually fine tuning that "spheres from csv" function as we speak. I will push the code to a new branch on the porespy repo and you can browse it. Just give me a few days.

[Camelot18]

Dear Professor,

Thanks again for your quick reply.
So the openpnm will directly generate a pore network based on the input distribution while porespy creates an image of spheres controlled by porosity and distribution function, and the rest part of the image will be extracted via snow to form the network? The distribution of spheres in the image would determine the pore throat distribution (or they are the same), right?

I tried porespy before in a lower version where I used op.io.PoreSpy.import_data(snow.network) to bring the parameters into openpnm. I obtained all the pore diameters, throat diameters, and throat lengths. But after I updated porespy/openpnm, and used op.io.network_from_porespy(snow.network) to do the same thing but only retuned with parameters like pore.inscribed_diameter. Should I use this example to get those parameters I need for openpnm?

Best,
Peter

[jgostick]

The first part of your above question is correct.

For the second part, we have stopped assuming that users want 'pore.inscribed_diameter' as their 'pore.diameter'. The user must do this themselves. This can be done with:

del pn.models['pore.diameter']
pn['pore.diameter'] = pn['pore.inscribed_diameter']

The first step removes any models that may be present which compute pore diameter, as these would over write the hand written values that you assign on the next step.

[Camelot18]

Thanks for your kind reply.

I want to build a pore network where the pore size distribution is a norm function(loc=50um,scale=10um) and the porosity is 0.4. The image size is 300300300, and the resolution is 44um/voxel. I tried in porespy via ps.generators.polydisperse_spheres and snow, the porosity and distribution, as well as the percentage of pores volume in the total pore space are not correct. Could you please help me check it?

import matplotlib.pyplot as plt
import numpy as np
import openpnm as op
import porespy as ps
import scipy.stats as spst

shape = [300, 300,300]
porosity = 0.4
resolution = 44e-6
dist = spst.norm(loc=50e-6/resolution, scale=10e-6/resolution)
nbins = 5
r_min = 10e-6/resolution

im = ps.generators.polydisperse_spheres(shape, porosity, dist,nbins,r_min)
graphx = ps.visualization.show_3D(im)
plt.imshow(im[:, :, 50]);
im_poro = ps.metrics.porosity(im)
print(f"Porosity from image: {im_poro*100:.1f}%")

fig = plt.subplots(figsize=[4, 4])
plt.imshow(graphx)
plt.axis(True)
plt.title(f'porosity = {im.sum()/im.size:.2f}')

snow = ps.networks.snow2(im, voxel_size=1)
net = op.io.network_from_porespy(snow.network)
del net.models['pore.diameter']
del net.models['throat.diameter']
del net.models['throat.coords']

net['pore.diameter'] = net['pore.inscribed_diameter']
net['throat.diameter'] = net['throat.inscribed_diameter']
net['throat.coords'] = net['throat.global_peak']
net.add_model(propname='throat.length',
             model=op.models.geometry.throat_length.hybrid_cones_and_cylinders)
conns = net['throat.conns']
model=op.models.geometry.pore_volume.sphere
net.add_model(propname='pore.volume',
             model=model,
             regen_mode='normal')
model=op.models.geometry.throat_volume.cylinder
net.add_model(propname='throat.volume',
             model=model,
             regen_mode='normal')
net.regenerate_models()
print(net)

fig, ax = plt.subplots(1,1,figsize=[8, 8])
ax = op.visualization.plot_coordinates(ax=fig,
                              network=net,
                              size_by=net["pore.inscribed_diameter"],
                              color_by=net["pore.inscribed_diameter"],
                              markersize=200)
ax = op.visualization.plot_connections(network=net,ax=fig,
                                   size_by=net["throat.inscribed_diameter"],
                              color_by=net["throat.inscribed_diameter"])

# PSD histogram
plt.figure(figsize=(12, 8))

plt.subplot(221)
plt.hist(net['pore.diameter']*44e-6, bins=15, edgecolor='k')
plt.title('Pore Diameter')

plt.subplot(222)
plt.hist(net['pore.coordination_number'], bins=15, edgecolor='k')
plt.title('Pore Coordination Number')

plt.subplot(223)
plt.hist(net['throat.diameter']*44e-6, bins=15, edgecolor='k')
plt.title('Throat Diameter')

plt.subplot(224)
plt.hist(net['throat.length']*44e-6, bins=15, edgecolor='k')
plt.title('Throat Length')
plt.show()

Vol_void = np.sum(net['pore.volume'])+np.sum(net['throat.volume'])
print(Vol_void)
Vol_bulk = 300**3 # from the image
pnm_poro = Vol_void / Vol_bulk
print(f"Porosity from pnm: {pnm_poro*100:.1f}%")

print( f"pores percentage:{np.sum(net['pore.volume'])/Vol_void*100:.2f}%")
print(f"throats percentage:{np.sum(net['throat.volume'])/Vol_void*100:.2f}%")

However, the pores percentage is about 16%, while throats account for 84%. The pores are supposed to be dominant in the pore space.

Best,
Peter

[Camelot18]

I tried to add diameters for the newly added pores, but it didn't work, can you please check my code below?

op.topotools.extend(network=pn, coords=hydrate_coords , labels=['hydrate'])
pn.add_model(propname='hydrate.diameter',
               model=0.5*pn['pore.diameter'],
               domain = 'hydrate', regen_mode='normal')

Thanks in advance.

[jgostick]

The problem above is that you're not assigning an actual model to the modle argument. See my above comment for the correct way to do what you're trying to do.

[Camelot18]

Thanks, I think I figured it out.

[Camelot18]

Actually, there are still some problems when I add more pores with diameters; the code below shows how I deal with it.

op.topotools.extend(network=net, coords=[(x,y,z)], labels=['hydrate'])
del net.models['pore.diameter']
net['pore.hydrate'] = r_h3*net['pore.diameter']

The size of the returned net['pore.hydrate'] is twice of net['pore.diameter'], and half of it is Nans. How to solve this?