Part 1: Import Networks from Statoil Files#

This example explains how to use the OpenPNM.Utilies.IO.Statoil class to import a network produced by the Maximal Ball network extraction code developed by Martin Blunt’s group at Imperial College London. The code is available from him upon request, but they offer a small library of pre-extracted networks on their website.

[1]:
import warnings
import scipy as sp
import numpy as np
import openpnm as op
%config InlineBackend.figure_formats = ['svg']
np.set_printoptions(precision=4)
np.random.seed(10)
%matplotlib inline

The following assumes that the folder containing the ‘dat’ files is in a directory called ‘fixtures’ in the same directory as this script. You can also enter a full path to the files.

[2]:
from pathlib import Path
path = Path('../_fixtures/ICL-Sandstone(Berea)/')
project = op.io.Statoil.import_data(path=path, prefix='Berea')
pn = project.network
pn.name = 'berea'
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [2], in <cell line: 3>()
      1 from pathlib import Path
      2 path = Path('../_fixtures/ICL-Sandstone(Berea)/')
----> 3 project = op.io.Statoil.import_data(path=path, prefix='Berea')
      4 pn = project.network
      5 pn.name = 'berea'

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

This import class extracts all the information contained in the ‘Statoil’ files, such as sizes, locations and connectivity. Note that the io classes return a project object, and the network itself can be accessed using the network attribute. The following printout display which information was contained in the file:

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

NameError: name 'pn' is not defined

At this point, the network can be visualized in Paraview. A suitable ‘.vtp’ file can be created with:

[4]:
op.io.VTK.export_data(network=pn, filename='imported_statoil')
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [4], in <cell line: 1>()
----> 1 op.io.VTK.export_data(network=pn, filename='imported_statoil')

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

The resulting network is shown below:

7bef11ee84d54d3f8ad3ea557c6e7daf

Clean up network topology#

Although it’s not clear in the network image, there are a number of isolated and disconnected pores in the network. These are either naturally part of the sandstone, or artifacts of the Maximal Ball algorithm. In any event, these must be removed before proceeding since they cause problems for the matrix solver. The easiest way to find these is to use the check_network_health function on the network object. This will return a dictionary with several key attributes, including a list of which pores are isolated. These can then be trimmed using the trim function in the topotools module.

[5]:
print('Number of pores before trimming: ', pn.Np)
h = pn.check_network_health()
op.topotools.trim(network=pn, pores=np.hstack([h['isolated_pores'],h['disconnected_pores']]))
print('Number of pores after trimming: ', pn.Np)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [5], in <cell line: 1>()
----> 1 print('Number of pores before trimming: ', pn.Np)
      2 h = pn.check_network_health()
      3 op.topotools.trim(network=pn, pores=np.hstack([h['isolated_pores'],h['disconnected_pores']]))

NameError: name 'pn' is not defined

Dealing with Inlet and Outlet Pores#

When importing Statoil networks, OpenPNM must perform some ‘optimizations’ to make the network compatible. The main problem is that the original network contains a large number of throats connecting actual internal pores to fictitious ‘reservoir’ pores. OpenPNM strips away all these throats since ‘headless throats’ break the graph theory representation. OpenPNM then labels the real internal pores as either ‘inlet’ or ‘outlet’ if they were connected to one of these fictitious reservoirs.

It is fairly simple to add a new pores to each end of the domain and stitch tehm to the internal pores labelled ‘inlet’ and ‘outlet’, but this introduces a subsequent complication that the new pores don’t have any geometry properties. For this example, we will not add boundary pores, but just the pores on the inlet and outlet faces.

Part 2: Calculating Permeability of the Network#

Setup Geometry, Phase, and Physics Objects#

In OpenPNM 2+ it is optional to define Geometry and Physics objects (These are really only necessary for simulations with diverse geometrical properties in different regions, resulting in different physical processes in each region, such as multiscale networks for instance). It is still necessary to define Phase objects:

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

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

Apply Pore-Scale Models#

We must add the hagen-poiseuille model for calculating the conductance. In OpenPNM 2+ it is possible to add Physics models to Phase objects, which is often simpler than than applying the same model to multiple Physics.

[7]:
water.add_model(propname='throat.hydraulic_conductance',
                model=op.models.physics.hydraulic_conductance.valvatne_blunt)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [7], in <cell line: 1>()
----> 1 water.add_model(propname='throat.hydraulic_conductance',
      2                 model=op.models.physics.hydraulic_conductance.valvatne_blunt)

NameError: name 'water' is not defined

Recall that boundary pores and throats had no geometrical properties associated with them, so the hydraulic conductances of boundary throats will be undefined (filled with NaNs):

[8]:
print(water['throat.hydraulic_conductance'])
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [8], in <cell line: 1>()
----> 1 print(water['throat.hydraulic_conductance'])

NameError: name 'water' is not defined

Run StokesFlow Algorithm#

Finally, we can create a StokesFlow object to run some fluid flow simulations:

[9]:
flow = op.algorithms.StokesFlow(network=pn, phase=water)
flow.set_value_BC(pores=pn.pores('inlets'), values=200000)
flow.set_value_BC(pores=pn.pores('outlets'), values=100000)
flow.run()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [9], in <cell line: 1>()
----> 1 flow = op.algorithms.StokesFlow(network=pn, phase=water)
      2 flow.set_value_BC(pores=pn.pores('inlets'), values=200000)
      3 flow.set_value_BC(pores=pn.pores('outlets'), values=100000)

NameError: name 'pn' is not defined

The resulting pressure field can be visualized in Paraview, giving the following:

0c40bc5a0a73435e8523ac6ee5efcc8a

Determination of Permeability Coefficient#

The way to calculate K is the determine each of the values in Darcy’s law manually and solve for K, such that

\[K = \frac{Q\mu L} {\Delta P A}\]
[10]:
# Get the average value of the fluid viscosity
mu = np.mean(water['pore.viscosity'])
# Specify a pressure difference (in Pa)
delta_P = 100000
# Using the rate method of the StokesFlow algorithm
Q = np.absolute(flow.rate(pores=pn.pores('inlets')))
# Because we know the inlets and outlets are at x=0 and x=X
Lx = np.amax(pn['pore.coords'][:, 0]) - np.amin(pn['pore.coords'][:, 0])
A = Lx*Lx  # Since the network is cubic Lx = Ly = Lz
K = Q*mu*Lx/(delta_P*A)
print(K)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [10], in <cell line: 2>()
      1 # Get the average value of the fluid viscosity
----> 2 mu = np.mean(water['pore.viscosity'])
      3 # Specify a pressure difference (in Pa)
      4 delta_P = 100000

NameError: name 'water' is not defined