Predicting Permeability of Berea#

Berea Sandstone Simulation Using PoreSpy and OpenPNM#

The example explains effective permeabilty calculations using PoreSpy and OpenPNM software. The simulation is performed on X-ray tomography image of BereaSandstone. The calculated effective permeablity value can compared with value report in Dong et al.

Start by importing the necessary packages#

[1]:
import os
import imageio
import scipy as sp
import numpy as np
import openpnm as op
import porespy as ps
import matplotlib.pyplot as plt
np.set_printoptions(precision=4)
np.random.seed(10)
%matplotlib inline

Load BreaSandstone Image file#

Give path to image file and load the image. Please note image should be binarized or in boolean format before performing next steps.

[2]:
path = '../../_fixtures/ICL-Sandstone(Berea)/'
file_format = '.tif'
file_name = 'Berea'
file = file_name + file_format
fetch_file = os.path.join(path, file)
im = imageio.mimread(fetch_file)
im = ~np.array(im, dtype=bool)[:250, :250, :250]  # Make image a bit smaller
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Input In [2], in <cell line: 6>()
      4 file = file_name + file_format
      5 fetch_file = os.path.join(path, file)
----> 6 im = imageio.mimread(fetch_file)
      7 im = ~np.array(im, dtype=bool)[:250, :250, :250]

File /usr/share/miniconda3/envs/test/lib/python3.8/site-packages/imageio/v2.py:292, in mimread(uri, format, memtest, **kwargs)
    290 imopen_args = decypher_format_arg(format)
    291 imopen_args["legacy_mode"] = True
--> 292 with imopen(uri, "rI", **imopen_args) as file:
    293     for image in file.iter(**kwargs):
    294         images.append(image)

File /usr/share/miniconda3/envs/test/lib/python3.8/site-packages/imageio/core/imopen.py:110, in imopen(uri, io_mode, plugin, format_hint, legacy_mode, **kwargs)
    108     request.format_hint = format_hint
    109 else:
--> 110     request = Request(uri, io_mode, format_hint=format_hint)
    112 source = "<bytes>" if isinstance(uri, bytes) else uri
    114 # fast-path based on plugin
    115 # (except in legacy mode)

File /usr/share/miniconda3/envs/test/lib/python3.8/site-packages/imageio/core/request.py:248, in Request.__init__(self, uri, mode, format_hint, **kwargs)
    245     raise ValueError(f"Invalid Request.Mode: {mode}")
    247 # Parse what was given
--> 248 self._parse_uri(uri)
    250 # Set extension
    251 if self._filename is not None:

File /usr/share/miniconda3/envs/test/lib/python3.8/site-packages/imageio/core/request.py:388, in Request._parse_uri(self, uri)
    385 if is_read_request:
    386     # Reading: check that the file exists (but is allowed a dir)
    387     if not os.path.exists(fn):
--> 388         raise FileNotFoundError("No such file: '%s'" % fn)
    389 else:
    390     # Writing: check that the directory to write to does exist
    391     dn = os.path.dirname(fn)

FileNotFoundError: No such file: '/home/runner/work/OpenPNM/OpenPNM/docs/_fixtures/ICL-Sandstone(Berea)/Berea.tif'

Confirm image and check image porosity#

Be patient, this might take ~30 seconds (depending on your CPU)

[3]:
# NBVAL_IGNORE_OUTPUT
fig, ax = plt.subplots(1, 3, figsize=(12,5))
ax[0].imshow(im[:, :, 100]);
ax[1].imshow(ps.visualization.show_3D(im));
ax[2].imshow(ps.visualization.sem(im));
ax[0].set_title("Slice No. 100 View");
ax[1].set_title("3D Sketch");
ax[2].set_title("SEM View");
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [3], in <cell line: 3>()
      1 # NBVAL_IGNORE_OUTPUT
      2 fig, ax = plt.subplots(1, 3, figsize=(12,5))
----> 3 ax[0].imshow(im[:, :, 100]);
      4 ax[1].imshow(ps.visualization.show_3D(im));
      5 ax[2].imshow(ps.visualization.sem(im));

NameError: name 'im' is not defined
../../_images/examples_contrib_predicting_effective_permeability_of_berea_8_1.png
[4]:
print(ps.metrics.porosity(im))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 print(ps.metrics.porosity(im))

NameError: name 'im' is not defined

Extract pore network using SNOW algorithm in PoreSpy#

The SNOW algorithm (an accronym for Sub-Network from an Oversegmented Watershed) was presented by Gostick. The algorithm was used to extract pore network from BereaSandstone image.

[5]:
# NBVAL_IGNORE_OUTPUT
resolution = 5.345e-6
snow = ps.networks.snow2(im, voxel_size=resolution)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [5], in <cell line: 3>()
      1 # NBVAL_IGNORE_OUTPUT
      2 resolution = 5.345e-6
----> 3 snow = ps.networks.snow2(im, voxel_size=resolution)

NameError: name 'im' is not defined

Import network in OpenPNM#

The output from the SNOW algorithm above is a plain python dictionary containing all the extracted pore-scale data, but it is NOT yet an OpenPNM network. We need to create an empty network in OpenPNM, then populate it with the data from SNOW:

[6]:
settings = {'pore_shape': 'pyramid',
            'throat_shape': 'cuboid',
            'pore_diameter': 'equivalent_diameter',
            'throat_diameter': 'inscribed_diameter'}
pn, geo = op.io.PoreSpy.import_data(snow.network, settings=settings)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [6], in <cell line: 5>()
      1 settings = {'pore_shape': 'pyramid',
      2             'throat_shape': 'cuboid',
      3             'pore_diameter': 'equivalent_diameter',
      4             'throat_diameter': 'inscribed_diameter'}
----> 5 pn, geo = op.io.PoreSpy.import_data(snow.network, settings=settings)

AttributeError: module 'openpnm.io' has no attribute 'PoreSpy'

Now we can print the network to see how the transferred worked.

Note to developers: We need to ignore the output of the following cell since the number of pores differs depending on whether the code is run on a windows or linux machine.

[7]:
# NBVAL_IGNORE_OUTPUT
print(pn)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [7], in <cell line: 2>()
      1 # NBVAL_IGNORE_OUTPUT
----> 2 print(pn)

NameError: name 'pn' is not defined

Check network health#

Remove isolated pores or cluster of pores from the network by checking it network health. Make sure ALL keys in network health functions have no value.

[8]:
h = pn.check_network_health()
op.topotools.trim(network=pn, pores=h['trim_pores'])
h = pn.check_network_health()
print(h)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 h = pn.check_network_health()
      2 op.topotools.trim(network=pn, pores=h['trim_pores'])
      3 h = pn.check_network_health()

NameError: name 'pn' is not defined

Assign phase#

In this example air is considered as fluid passing through porous channels.

[9]:
air = op.phases.Air(network=pn)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [9], in <cell line: 1>()
----> 1 air = op.phases.Air(network=pn)

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

Assign physics#

[10]:
phys_air = op.physics.Basic(network=pn, phase=air, geometry=geo)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [10], in <cell line: 1>()
----> 1 phys_air = op.physics.Basic(network=pn, phase=air, geometry=geo)

NameError: name 'pn' is not defined

Assign Algorithm and boundary conditions#

Select stokes flow algorithm for simulation and assign dirichlet boundary conditions in top and bottom faces of the network.

[11]:
perm = op.algorithms.StokesFlow(network=pn, phase=air)
perm.set_value_BC(pores=pn.pores('zmax'), values=0)
perm.set_value_BC(pores=pn.pores('zmin'), values=101325)
perm.run()
air.update(perm.results())
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [11], in <cell line: 1>()
----> 1 perm = op.algorithms.StokesFlow(network=pn, phase=air)
      2 perm.set_value_BC(pores=pn.pores('zmax'), values=0)
      3 perm.set_value_BC(pores=pn.pores('zmin'), values=101325)

NameError: name 'pn' is not defined

Calculate effective permeability#

Caclulate effective permeablity using hagen poiseuille equation. Use cross section area and flow length manually from image dimension.

[12]:
resolution = 5.345e-6
Q = perm.rate(pores=pn.pores('zmin'), mode='group')[0]
A = (im.shape[0] * im.shape[1]) * resolution**2
L = im.shape[2] * resolution
mu = air['pore.viscosity'].max()
delta_P = 101325 - 0
K = Q * L * mu / (A * delta_P)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [12], in <cell line: 2>()
      1 resolution = 5.345e-6
----> 2 Q = perm.rate(pores=pn.pores('zmin'), mode='group')[0]
      3 A = (im.shape[0] * im.shape[1]) * resolution**2
      4 L = im.shape[2] * resolution

NameError: name 'perm' is not defined

Note to developers: We need to ignore the output of the following cell since the results are slightly different on different platforms (windows vs linux)

[13]:
# NBVAL_IGNORE_OUTPUT
print(f'The value of K is: {K/0.98e-12*1000:.2f} mD')
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [13], in <cell line: 2>()
      1 # NBVAL_IGNORE_OUTPUT
----> 2 print(f'The value of K is: {K/0.98e-12*1000:.2f} mD')

NameError: name 'K' is not defined