pc_curve

This single method can convert either (a) images containing invasion size or (b) invasion pressure results, such as those produced by porosimetry and drainage, respectively.

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

/opt/hostedtoolcache/Python/3.8.16/x64/lib/python3.8/site-packages/openpnm/algorithms/_invasion_percolation.py:358: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def _find_trapped_pores(inv_seq, indices, indptr, outlets):  # pragma: no cover

Generate an image and apply the porosimetry function:

np.random.seed(0)
im = ps.generators.blobs(shape=[200, 200], porosity=0.6)
drainage = ps.filters.porosimetry(im)

sizes

The porosimetry function returns an image indicating the radius of the capillary meniscus that at which is was invaded. In this case the Washburn equation is used to convert sizes to capillary pressures, surface tension (sigma) and contact (theta) must be provided, as well as voxel size (voxel_size) if the sizes array is in units of voxels.

data = ps.metrics.pc_curve(im=im, sizes=drainage, voxel_size=1e-5)

The function results a Results object, whose attributes can be inspected by printing:

print(data)

――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
Results of pc_curve generated at Tue Jun  6 13:56:25 2023
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
pc                        [407.8608703277795, 815.721740655559, 924.1694802519894, 1047.0350190039558, 1186.2351597260008, 1343.9415383726086, 1522.6145033338332, 1725.041498880861, 1954.3805515746424, 2214.2095380610917, 2508.582002870727, 2842.090397026047, 3219.9377240305935, 3648.018710975678, 4133.011771100392, 4682.483192495964, 5305.005178383415, 6809.339411869592, 7714.620502905222, 11218.720341000993, 14399.999999999998]
snwp                      [0.0, 0.04955537590945837, 0.07877930476960388, 0.0881568310428456, 0.09280517380759903, 0.21883589329021827, 0.2791430881164107, 0.3557801131770412, 0.4190784155214228, 0.6601455133387227, 0.8425222312045271, 0.8719482619240096, 0.9545675020210186, 0.9654405820533549, 0.9690379951495554, 0.9706548100242522, 0.9867016976556184, 0.9870654810024252, 0.9881568310428456, 0.9896928051738076, 0.9897736459175425]
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

Matplotlib can be used to plot the curve:

plt.plot(data.pc, data.snwp, 'b-o')
plt.xlabel('Capillary Pressure [Pa]')
plt.ylabel('Non-wetting Phase Saturation');

Or, other options:

fig, ax = plt.subplots(1, 2, figsize=[12, 6])
ax[0].step(data.pc, data.snwp, 'b-o', where='post')
ax[0].set_xlabel('Capillary Pressure [Pa]')
ax[0].set_ylabel('Non-wetting Phase Saturation')
ax[1].semilogx(data.pc, data.snwp, 'r-o')
ax[1].set_xlabel('Capillary Pressure [Pa]')
ax[1].set_ylabel('Non-wetting Phase Saturation');

pc

In this case no additional information is needed since the sizes have been converted to capillary pressure already, either by a previous function, or by hand.

sigma = 0.072
theta = 180
voxel_size = 1e-5
pc = -2*sigma*np.cos(np.deg2rad(theta))/(drainage*voxel_size)

data = ps.metrics.pc_curve(im=im, pc=pc)
plt.plot(data.pc, data.snwp, 'b-o')
plt.xlabel('Capillary Pressure [Pa]')
plt.ylabel('Non-wetting Phase Saturation');

Note that some voxels in pc had a value of np.inf since they were uninvaded. The pc_curve function replaces these with 2*pressure.max() to improve the plotting.