Ordinary Percolation#

OpenPNM contains several percolation algorithms which are central to the multiphase models employed by pore networks. The essential idea is to identify pathways for fluid flow through the network using the entry capillary pressure as a threshold for passage between connected pores. The capillary pressure can either be associated to the pores themselves known as site percolation or the connecting throats known as bond percolation or a mixture of both. OpenPNM provides several models for calculating the entry pressure for a given pore or throat and it generally depends on the size of the pore or throat and the wettability to a particular phase characterised by the contact angle. If a pathway through the network connects pores into clusters that contain both an inlet and an outlet then it is deemed to be percolating.

In this example we will demonstrate Ordinary Percolation which is the fastest and simplest algorithm to run. The number of steps involved in the algorithm is equal to the number of points that are specified in the run method. This can either be an integer, in which case the minimum and maximum capillary entry pressures in the network are used as limits and the integer value is used to create that number of intervals between the limits, or an array of specified pressured can be supplied.

The algorithm progresses incrementally from low pressure to high. At each step, clusters of connected pores are found with entry pressures below the current threshold and those that are not already invaded and connected to an inlet are set to be invaded at this pressure. Therefore the process is quasistatic and represents the steady state saturation that would be achieved if the inlet pressure were to be held at that threshold.

First do our imports

[1]:
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
np.random.seed(10)
from ipywidgets import interact, IntSlider
%matplotlib inline
mpl.rcParams["image.interpolation"] = "None"
ws = op.Workspace()
ws.settings["loglevel"] = 40

Create a 2D Cubic network with standard PSD and define the phase as Water and use Standard physics which implements the washburn capillary pressure relation for throat entry pressure.

[2]:
N = 100
net = op.network.Cubic(shape=[N, N, 1], spacing=2.5e-5)
geom = op.geometry.SpheresAndCylinders(network=net, pores=net.Ps, throats=net.Ts)
water = op.phases.Water(network=net)
phys = op.physics.Standard(network=net, phase=water, geometry=geom)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [2], in <cell line: 4>()
      2 net = op.network.Cubic(shape=[N, N, 1], spacing=2.5e-5)
      3 geom = op.geometry.SpheresAndCylinders(network=net, pores=net.Ps, throats=net.Ts)
----> 4 water = op.phases.Water(network=net)
      5 phys = op.physics.Standard(network=net, phase=water, geometry=geom)

AttributeError: module 'openpnm' has no attribute 'phases'

We can check the model by looking at the model dict on the phys object

[3]:
phys.models['throat.entry_pressure']
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [3], in <cell line: 1>()
----> 1 phys.models['throat.entry_pressure']

NameError: name 'phys' is not defined

Now set up and run the algorithm choosing the left and right sides of the network for inlets and outlets respectively. Because we did not set up the network with boundary pores with zero volume a little warning is given because the starting saturation for the algorithm is not zero. However, this is fine and because the network is quite large the starting saturation is actually quite close to zero.

[4]:
alg = op.algorithms.OrdinaryPercolation(network=net, phase=water)
alg.settings._update({'pore_volume': 'pore.volume',
                      'throat_volume': 'throat.volume'})
alg.set_inlets(pores=net.pores('left'))
alg.set_outlets(pores=net.pores('right'))
alg.run(points=1000)
alg.plot_intrusion_curve()
plt.show()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 alg = op.algorithms.OrdinaryPercolation(network=net, phase=water)
      2 alg.settings._update({'pore_volume': 'pore.volume',
      3                       'throat_volume': 'throat.volume'})
      4 alg.set_inlets(pores=net.pores('left'))

NameError: name 'water' is not defined

The algorithm completes very quickly and the invading phase saturation can be plotted versus the applied boundary pressure.

[5]:
data = alg.get_intrusion_data()
mask = np.logical_and(np.asarray(data.Snwp) > 0.0 , np.asarray(data.Snwp) < 1.0)
mask = np.argwhere(mask).flatten()
pressures = np.asarray(data.Pcap)[mask]
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [5], in <cell line: 1>()
----> 1 data = alg.get_intrusion_data()
      2 mask = np.logical_and(np.asarray(data.Snwp) > 0.0 , np.asarray(data.Snwp) < 1.0)
      3 mask = np.argwhere(mask).flatten()

NameError: name 'alg' is not defined

As the network is 2D and cubic we can easily plot the invading phase configuration at the different invasion steps

[6]:
def plot_saturation(step):
    arg = mask[step]
    Pc = np.ceil(data.Pcap[arg])
    sat = np.around(data.Snwp[arg], 3)
    is_perc = alg.is_percolating(Pc)
    pmask = alg['pore.invasion_pressure'] <= Pc
    im = pmask.reshape([N, N])
    fig, ax = plt.subplots(figsize=[5, 5])
    ax.imshow(im, cmap='Blues');
    title = f"Capillary pressure: {Pc:.0f}, saturation: {sat:.2f}, percolating: {is_perc}"
    ax.set_title(title)
    plt.show()
[7]:
perc_thresh = alg.get_percolation_threshold()
thresh_step = np.argwhere(np.asarray(pressures) == perc_thresh)
slider = IntSlider(min=0, max=len(mask)-1, step=1, value=thresh_step)
interact(plot_saturation, step=slider);
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [7], in <cell line: 1>()
----> 1 perc_thresh = alg.get_percolation_threshold()
      2 thresh_step = np.argwhere(np.asarray(pressures) == perc_thresh)
      3 slider = IntSlider(min=0, max=len(mask)-1, step=1, value=thresh_step)

NameError: name 'alg' is not defined