Invasion pressure < entry pressure for throats in output from algorithms.Porosimetry

[1minus1]

Apologies if this is due to a lack of understanding of the algorithm, but in analyzing the output of a porosimetry simulation on a network extracted from a porespy cylinders image, I am finding that throats are being opened at pressures less than their entry_pressure as defined with the phase physics and associated model. The plot generated by the code below shows the throat.invasion_pressure and throat.entry_pressure as a function of throat.invasion_sequence. I would think that for all throats opened at a pressure step, the entry_pressure would have to be <= invasion_pressure. Blue points on the generated plot that sit above the red points represent points invaded at pressures less than their entry pressure.

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

im = ps.generators.cylinders([200,200,200],5,250,0,90,None)
print(ps.metrics.porosity(im))

snow_output = ps.networks.snow(im,
                    voxel_size=1,
                    boundary_faces=['left','right','front','back','top','bottom'],
                    marching_cubes_area=False)
pore_locations = snow_output['pore.coords']
pore_diameters = snow_output['pore.diameter']
throat_diameters = snow_output['throat.diameter']

ws = op.Workspace()
prj = ws.new_project()
pn = op.network.GenericNetwork(project=prj)
pn.update(snow_output)
geo = op.geometry.GenericGeometry(project=prj,pores=pn.Ps,throats=pn.Ts)

hg = op.phases.Mercury(network = pn)
phys_hg = op.physics.GenericPhysics(network=pn, phase=hg, geometry=geo)
phys_hg.add_model(propname='throat.entry_pressure',
                model=op.models.physics.capillary_pressure.washburn,
                contact_angle='pore.contact_angle',
                surface_tension='pore.surface_tension')

mip = op.algorithms.Porosimetry(project=prj)

throat_entry_pressure = phys_hg["throat.entry_pressure"]
pressure_points = np.logspace(np.log10(0.5*np.min(throat_entry_pressure)), np.log10(2*np.max(throat_entry_pressure)),50)

mip.setup(phase = hg)
mip.set_inlets(pn.pores(['top']),overwrite=True)
mip.run(points = pressure_points)
these_results=mip.results()

throat_invasion_pressure = these_results['throat.invasion_pressure']
throat_invasion_sequence = mip['throat.invasion_sequence']

plt.scatter(throat_invasion_sequence,throat_invasion_pressure,s=10,label='Invasion Pressure',c='r')
plt.xlabel('Throat invasion sequence')
plt.ylabel('Pressure')
plt.scatter(throat_invasion_sequence,throat_entry_pressure,s=5,label='Entry Pressure',c='b')
plt.legend()
plt.show()

[1minus1]

I added some visualization to show how the connectivity of opened pores and throats doesn't make sense for an intrusion process. Like in the image below, where there are open pores not connected to any open throats, and vice versa. In this example, percolation is coming from the top.

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

def plot_network(network, pore_mask=None, throat_mask=None,title=None):
    pore_locations = network['pore.coords']
    pore_diameters = network['pore.diameter']
    throat_diameters = network['throat.diameter']
    throat_conns = network['throat.conns']

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')

    minx = np.min([loc[0] for loc in pore_locations])
    miny = np.min([loc[1] for loc in pore_locations])
    minz = np.min([loc[2] for loc in pore_locations])

    maxx = np.max([loc[0] for loc in pore_locations])
    maxy = np.max([loc[1] for loc in pore_locations])
    maxz = np.max([loc[2] for loc in pore_locations])

    rangex = maxx-minx
    rangey = maxy-miny
    rangez = maxz-minz

    ax.set_xlim3d(minx-.1*rangex,maxx+.1*rangex)
    ax.set_ylim3d(miny-.1*rangey,maxy+.1*rangey)
    ax.set_zlim3d(minz-.1*rangez,maxz+.1*rangez)

    for thispoint,thisdiam,masked in zip(pore_locations,pore_diameters,pore_mask):
        if masked:
            ax.scatter(thispoint[0],thispoint[1],thispoint[2],s=2*thisdiam,c='b')
            ax.view_init(elev=1., azim=45)

    for conn,diam,masked in zip(throat_conns,throat_diameters,throat_mask):
        if masked:
            x = [pore_locations[conn[0]][0],pore_locations[conn[1]][0]]
            y = [pore_locations[conn[0]][1],pore_locations[conn[1]][1]]
            z = [pore_locations[conn[0]][2],pore_locations[conn[1]][2]]
            ax.plot3D(x,y,z,c='r',lw=2)
    if title is not None:
        plt.title(title)
    plt.show()

im = ps.generators.cylinders([100,100,100],5,100,0,90,None)
print(ps.metrics.porosity(im))

snow_output = ps.networks.snow(im,
                    voxel_size=1,
                    boundary_faces=['left','right','front','back','top','bottom'],
                    marching_cubes_area=False)
pore_locations = snow_output['pore.coords']
pore_diameters = snow_output['pore.diameter']
throat_diameters = snow_output['throat.diameter']

ws = op.Workspace()
prj = ws.new_project()
pn = op.network.GenericNetwork(project=prj)
pn.update(snow_output)
geo = op.geometry.GenericGeometry(project=prj,pores=pn.Ps,throats=pn.Ts)

hg = op.phases.Mercury(network = pn)
phys_hg = op.physics.GenericPhysics(network=pn, phase=hg, geometry=geo)
phys_hg.add_model(propname='throat.entry_pressure',
                model=op.models.physics.capillary_pressure.washburn,
                contact_angle='pore.contact_angle',
                surface_tension='pore.surface_tension')

mip = op.algorithms.Porosimetry(project=prj)

throat_entry_pressure = phys_hg["throat.entry_pressure"]
pressure_points = np.logspace(np.log10(0.5*np.min(throat_entry_pressure)), np.log10(2*np.max(throat_entry_pressure)),50)

mip.setup(phase = hg)
mip.set_inlets(pn.pores(['top']),overwrite=True)
mip.run(points = pressure_points)
these_results=mip.results()

throat_invasion_pressure = these_results['throat.invasion_pressure']
throat_invasion_sequence = mip['throat.invasion_sequence']
pore_invasion_sequence = mip['pore.invasion_sequence']

plt.scatter(throat_invasion_sequence,throat_invasion_pressure,s=10,label='Invasion Pressure',c='r')
plt.xlabel('Throat invasion sequence')
plt.ylabel('Pressure')
plt.scatter(throat_invasion_sequence,throat_entry_pressure,s=5,label='Entry Pressure',c='b')
plt.legend()
plt.show()

charttitle = lambda poreseq, throatseq: "poreseq = {}, throatseq = {}".format(poreseq,throatseq)

throatseq=None
for poreseq in np.unique(pore_invasion_sequence):
    if throatseq is None:
        plot_network(pn,pore_mask=pore_invasion_sequence<=poreseq,throat_mask=throat_invasion_sequence<0, title=charttitle(poreseq,throatseq))
    else:
        plot_network(pn,pore_mask=pore_invasion_sequence<=poreseq,throat_mask=throat_invasion_sequence<=throatseq, title=charttitle(poreseq,throatseq))

    if throatseq is None:
        throatseq = 0
    else:
        throatseq += 1
    plot_network(pn,pore_mask=pore_invasion_sequence<=poreseq,throat_mask=throat_invasion_sequence<=throatseq, title=charttitle(poreseq,throatseq))

[jgostick]

Thanks for this extremely detailed report! I have to admit I don't quite understand your graph. I would've plotted is as invasion pressure vs entry pressure and looked for values below the 45 degree line. you seem to have a lot of throat data points so not sure what you're only plotting 25 invasion steps. Anyway, I'm trying to investigate on my own.

I am playing with it and have not seen any issues. Here is what I've tried:

1. Using a basic cubic network with a stick and ball geometry, I found no throats that were invaded at a pressure below their calculated entry pressure
2. Since you're doing this on an extracted network, I fudged the throat sizes so that some are bigger than their neighboring pores, which is a problem encountered in extracted networks sometimes.
3. When I do network extraction on a 3D blobs image, I AM finding some throats that are larger than their neighboring pores, but NOT finding any throats that are invaded too soon in the invasion simulation.

[jgostick]

Ahhh, I wasn't paying enough attention...you used OrdinaryPercolation, not InvasionPercolation. When I run this, I do get the error. Will keep digging.

[jgostick]

This is fixed in PR #1509.