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.


In [1]:
import warnings
import scipy as sp
import numpy as np
import openpnm as op
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.


In [2]:
from pathlib import Path
path = Path('../fixtures/ICL-Sandstone(Berea)/')
project = op.io.Statoil.load(path=path, prefix='Berea')
pn = project.network
pn.name = 'berea'

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:


In [3]:
print(pn)


――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
openpnm.network.GenericNetwork : berea
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
#     Properties                                    Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
1     pore.area                                      6298 / 6298 
2     pore.clay_volume                               6298 / 6298 
3     pore.coords                                    6298 / 6298 
4     pore.radius                                    6298 / 6298 
5     pore.shape_factor                              6298 / 6298 
6     pore.volume                                    6298 / 6298 
7     throat.area                                   12098 / 12098
8     throat.clay_volume                            12098 / 12098
9     throat.conduit_lengths.pore1                  12098 / 12098
10    throat.conduit_lengths.pore2                  12098 / 12098
11    throat.conduit_lengths.throat                 12098 / 12098
12    throat.conns                                  12098 / 12098
13    throat.length                                 12098 / 12098
14    throat.radius                                 12098 / 12098
15    throat.shape_factor                           12098 / 12098
16    throat.total_length                           12098 / 12098
17    throat.volume                                 12098 / 12098
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
#     Labels                                        Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
1     pore.all                                      6298      
2     pore.inlets                                   201       
3     pore.outlets                                  246       
4     throat.all                                    12098     
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――

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


In [4]:
project.export_data(filename='imported_statoil')

The resulting network is shown below:

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.


In [5]:
print('Number of pores before trimming: ', pn.Np)
h = pn.check_network_health()
op.topotools.trim(network=pn, pores=h['trim_pores'])
print('Number of pores after trimming: ', pn.Np)


Number of pores before trimming:  6298
Number of pores after trimming:  6004

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:


In [6]:
water = op.phases.Water(network=pn)

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.


In [7]:
water.add_model(propname='throat.hydraulic_conductance',
                model=op.models.physics.hydraulic_conductance.valvatne_blunt)

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):


In [8]:
print(water['throat.hydraulic_conductance'])


[1.1825e-13 1.6357e-13 4.3332e-14 ... 1.6230e-11 7.9766e-12 6.7575e-11]

Run StokesFlow Algorithm

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


In [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()

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

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} $$


In [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)


[1.1566e-12]