OpenPNM_Code

"""# Installing OpenPNM and PoreSpy"""

!pip install porespy #PoreSpy installation already contains OpenPNM

"""# Importing the required libraries"""

import os
import imageio
import scipy as sp
import numpy as np
import openpnm as op
import porespy as ps
import pandas as pd
import pickle as pk
import matplotlib.pyplot as plt
import math
from porespy.tools import randomize_colors
from skimage.morphology import binary_dilation
#ps.visualization.set_mpl_style()
op.visualization.set_mpl_style()
np.set_printoptions(precision=4)

"""#Importing the image into the environment"""

path ='/content/pores.tiff'
file_name = 'pores'
file_format = '.tif'
file = file_name + file_format
image = os.path.join(path, file)
fetch_file = os.path.join(path, file)
im = imageio.mimread('/content/pores.tiff', memtest=False)
im = np.array(im, dtype=bool)
print('The input Micro-CT image size is: {}'.format(im.shape))
fig, ax = plt.subplots(1, 3, figsize=(12,5))
ax[0].imshow(im[:, :, 300],cmap='gray');
ax[0].set_title("Slice No. 300 View");
ax[1].imshow(im[:, :, 200],cmap='gray');
ax[1].set_title("Slice No. 200 View");
ax[2].imshow(im[:, :, 100],cmap='gray');
ax[2].set_title("Slice No. 100 View");

"""# Using the Snow algorithm"""

resolution = 12e-6 #meters
net=ps.networks.snow2(im, voxel_size=resolution, accuracy='standard', sigma=0.4, r_max=4)

"""# Creation of the total pore network (connected and isolated)"""

proj_1 = op.Project()
pn_tudo = op.network.Network(name=file_name, project=proj_1)
pn_tudo.update(net.network) # Fills 'pn' with data from 'net'
print(pn_tudo)

"""# Creation of the geometry of the total pore and throat network (connected and isolated)"""

mods = op.models.collections.geometry.spheres_and_cylinders
pn_tudo.add_model_collection(mods)
pn_tudo.regenerate_models()
print(pn_tudo)

geo_tudo = pn_tudo
print(geo_tudo)

"""# Exporting geometry"""

op.io.project_to_vtk(proj_1, filename='network_tudo', fill_nans=None, fill_infs=None)

"""# Export to csv of the properties of interest from the network pores"""

pore_volume = 1e18*geo_tudo['pore.volume']
pore_volume = pd.DataFrame(pore_volume, columns = ['pore.volume'])

pore_extended_diameter = 1e06*geo_tudo['pore.extended_diameter']
pore_extended_diameter = pd.DataFrame(pore_extended_diameter, columns = ['pore.extended_diameter'])

pore_surface_area = 1e12*geo_tudo['pore.surface_area']
pore_surface_area = pd.DataFrame(pore_surface_area, columns = ['pore.surface_area'])

pore_equivalent_diameter = 1e06*geo_tudo['pore.equivalent_diameter']
pore_equivalent_diameter = pd.DataFrame(pore_equivalent_diameter, columns = ['pore.equivalent_diameter'])

pore_inscribed_diameter = 1e06*geo_tudo['pore.inscribed_diameter']
pore_inscribed_diameter = pd.DataFrame(pore_inscribed_diameter, columns = ['pore.inscribed_diameter'])

pore_coordination_number = geo_tudo['pore.coordination_number']
pore_coordination_number = pd.DataFrame(pore_coordination_number , columns = ['pore.coordination_number'])

pore_data_sheet_tudo = pd.concat([pore_inscribed_diameter, pore_equivalent_diameter, pore_volume, pore_surface_area, pore_coordination_number, pore_extended_diameter], axis=1, join="inner")
pore_data_sheet_tudo = pore_data_sheet_tudo.round(2)
names = ['pore.inscribed_diameter','pore.equivalent_diameter','pore.volume','pore.surface_area', 'pore_coordination_number', 'pore.extended_diameter']
np.savetxt("pore_properties.csv", pore_data_sheet_tudo, header=','.join(names), fmt='% 1.2f', delimiter=',') 
pore_data_sheet_tudo

"""# Export to csv of the properties of interest from the network throats"""

throat_volume = 1e18*geo_tudo['throat.volume']
throat_volume = pd.DataFrame(throat_volume, columns = ['throat.volume'])

throat_total_length = 1e06*geo_tudo['throat.total_length']
throat_total_length = pd.DataFrame(throat_total_length, columns = ['throat.total_length'])

throat_cross_sectional_area = 1e12*geo_tudo['throat.cross_sectional_area']
throat_cross_sectional_area = pd.DataFrame(throat_cross_sectional_area, columns = ['throat.cross_sectional_area'])

throat_equivalent_diameter = 1e06*geo_tudo['throat.equivalent_diameter']
throat_equivalent_diameter = pd.DataFrame(throat_equivalent_diameter, columns = ['throat.equivalent_diameter'])

throat_inscribed_diameter = 1e06*geo_tudo['throat.inscribed_diameter']
throat_inscribed_diameter = pd.DataFrame(throat_inscribed_diameter, columns = ['throat.inscribed_diameter'])

throat_data_sheet_tudo = pd.concat([throat_inscribed_diameter, throat_equivalent_diameter, throat_cross_sectional_area, throat_total_length, throat_volume], axis=1, join="inner")
throat_data_sheet_tudo = throat_data_sheet_tudo.round(2)
names = ['throat.inscribed_diameter','throat.equivalent_diameter','throat.cross_sectional_area','throat.total_length', 'throat.volume']
np.savetxt("throat_properties.csv", throat_data_sheet_tudo, header=','.join(names), fmt='% 1.2f', delimiter=',') 
throat_data_sheet_tudo

"""# Creating graphs of properties of interest"""

import seaborn as sns
#ploting my network details
fig = plt.figure(figsize=(20, 10))

plt.subplot(2,3,1)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e06*geo_tudo['pore.inscribed_diameter'], stat='percent', bins=20, color='Blue', cbar=True, cbar_kws=dict(shrink=.75), legend=True, label='Pore Inscribed Diameter')
plt.ylim(0,35)
plt.xlim(0,850)
plt.xlabel('Pore inscribed diameter (µm)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,2)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e06*geo_tudo['pore.equivalent_diameter'], stat='percent', bins=20, color='Blue', legend=True, cbar=True, cbar_kws=dict(shrink=.75), label = 'Pore Equivalent Diameter')
plt.ylim(0,35)
plt.xlim(0,850)
plt.xlabel('Pore equivalent diameter (µm)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,3)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e18*geo_tudo['pore.volume'], stat='percent', bins=20, color='Blue', legend=True, cbar=True, cbar_kws=dict(shrink=.75), label = 'Pore Volume')
#plt.ylim(0,100)
#plt.xlim(0,1e8)
plt.xlabel('Pore volume (µm³)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,4)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e06*geo_tudo['throat.inscribed_diameter'], stat='percent', bins=20, color='Blue', cbar=True, cbar_kws=dict(shrink=.75), legend=True, label='Throat Inscribed Diameter')
plt.ylim(0,35)
plt.xlim(0,850)
plt.xlabel('Throat inscribed diameter (µm)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,5)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e06*geo_tudo['throat.equivalent_diameter'], stat='percent', bins=20, color='Blue', legend=True, cbar=True, cbar_kws=dict(shrink=.75), label = 'Throat Equivalent Diameter')
plt.ylim(0,35)
plt.xlim(0,850)
plt.xlabel('Throat equivalent diameter (µm)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,6)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e18*geo_tudo['throat.volume'], stat='percent', bins=20, color='Blue', legend=True, cbar=True, cbar_kws=dict(shrink=.75), label = 'Throat Volume')
plt.ylim(0,100)
#plt.xlim(0,2e6)
plt.xlabel('Throat Volume (µm³)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.savefig("detailed plot.jpg")

"""# Creation of the second pore network"""

proj = op.Project()
pn = op.network.Network(name=file_name, project=proj)
pn.update(net.network) # Fills 'pn' with data from 'net'
print(pn)

mods = op.models.collections.geometry.spheres_and_cylinders
pn.add_model_collection(mods)
pn.regenerate_models()
print(pn)

"""# Exclusion of isolated pores"""

h = op.utils.check_network_health(pn)
op.topotools.trim(network=pn, pores=h['disconnected_pores'])
print(h)

"""# Creating the geometry of the connected pore networks"""

geo = pn
print(geo)

"""# Export of geometry with connected pores"""

op.io.project_to_vtk(proj, filename='connected_network', fill_nans=None, fill_infs=None)

"""# Export to csv of the properties of interest from the pores of the connected network"""

pore_volume = 1e18*geo['pore.volume']
pore_volume = pd.DataFrame(pore_volume, columns = ['pore.volume'])

pore_extended_diameter = 1e06*geo['pore.extended_diameter']
pore_extended_diameter = pd.DataFrame(pore_extended_diameter, columns = ['pore.extended_diameter'])

pore_surface_area = 1e12*geo['pore.surface_area']
pore_surface_area = pd.DataFrame(pore_surface_area, columns = ['pore.surface_area'])

pore_equivalent_diameter = 1e06*geo['pore.equivalent_diameter']
pore_equivalent_diameter = pd.DataFrame(pore_equivalent_diameter, columns = ['pore.equivalent_diameter'])

pore_inscribed_diameter = 1e06*geo['pore.inscribed_diameter']
pore_inscribed_diameter = pd.DataFrame(pore_inscribed_diameter, columns = ['pore.inscribed_diameter'])

pore_coordination_number = geo['pore.coordination_number']
pore_coordination_number = pd.DataFrame(pore_coordination_number , columns = ['pore.coordination_number'])

pore_data_sheet = pd.concat([pore_inscribed_diameter, pore_equivalent_diameter, pore_volume, pore_surface_area, pore_coordination_number, pore_extended_diameter], axis=1, join="inner")
pore_data_sheet = pore_data_sheet.round(2)
names = ['pore.inscribed_diameter','pore.equivalent_diameter','pore.volume','pore.surface_area','pore.coordination_number', 'pore.extended_diameter']
np.savetxt("connected_pores_properties.csv", pore_data_sheet, header=','.join(names), fmt='% 1.2f', delimiter=',') 
pore_data_sheet

"""# Export to csv of the properties of interest from the throats of the connected network"""

throat_volume = 1e18*geo['throat.volume']
throat_volume = pd.DataFrame(throat_volume, columns = ['throat.volume'])

throat_total_length = 1e06*geo['throat.total_length']
throat_total_length = pd.DataFrame(throat_total_length, columns = ['throat.total_length'])

throat_cross_sectional_area = 1e12*geo['throat.cross_sectional_area']
throat_cross_sectional_area = pd.DataFrame(throat_cross_sectional_area, columns = ['throat.cross_sectional_area'])

throat_equivalent_diameter = 1e06*geo['throat.equivalent_diameter']
throat_equivalent_diameter = pd.DataFrame(throat_equivalent_diameter, columns = ['throat.equivalent_diameter'])

throat_inscribed_diameter = 1e06*geo['throat.inscribed_diameter']
throat_inscribed_diameter = pd.DataFrame(throat_inscribed_diameter, columns = ['throat.inscribed_diameter'])

throat_data_sheet = pd.concat([throat_inscribed_diameter, throat_equivalent_diameter, throat_cross_sectional_area, throat_total_length, throat_volume], axis=1, join="inner")
throat_data_sheet = throat_data_sheet.round(2)
names = ['throat.inscribed_diameter','throat.equivalent_diameter','throat.cross_sectional_area','throat.total_length', 'throat.volume']
np.savetxt("connected_throat_properties.csv", throat_data_sheet, header=','.join(names), fmt='% 1.2f', delimiter=',') 
throat_data_sheet

"""# Creating graphs of the connected network properties of interest"""

import seaborn as sns
#ploting my network details
fig = plt.figure(figsize=(20, 10))

plt.subplot(2,3,1)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e06*geo['pore.inscribed_diameter'], stat='percent', bins=20, color='Blue', cbar=True, cbar_kws=dict(shrink=.75), legend=True, label='Pore Inscribed Diameter')
plt.ylim(0,35)
plt.xlim(0,900)
plt.xlabel('Pore inscribed diameter (µm)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,2)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e06*geo['pore.equivalent_diameter'], stat='percent', bins=20, color='Blue', legend=True, cbar=True, cbar_kws=dict(shrink=.75), label = 'Pore Equivalent Diameter')
plt.ylim(0,35)
plt.xlim(0,900)
plt.xlabel('Pore equivalent diameter (µm)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,3)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e18*geo['pore.volume'], stat='percent', bins=20, color='Blue', legend=True, cbar=True, cbar_kws=dict(shrink=.75), label = 'Pore Volume')
plt.ylim(0,100)
#plt.xlim(0,2e8)
plt.xlabel('Pore volume (µm³)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,4)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e06*geo['throat.inscribed_diameter'], stat='percent', bins=20, color='Blue', cbar=True, cbar_kws=dict(shrink=.75), legend=True, label='Throat Inscribed Diameter')
plt.ylim(0,35)
plt.xlim(0,900)
plt.xlabel('Throat inscribed diameter (µm)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,5)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e06*geo['throat.equivalent_diameter'], stat='percent', bins=20, color='Blue', legend=True, cbar=True, cbar_kws=dict(shrink=.75), label = 'Throat Equivalent Diameter')
plt.ylim(0,35)
plt.xlim(0,900)
plt.xlabel('Throat equivalent diameter (µm)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.subplot(2,3,6)
newBins = [0,56,60,65,70,80,85,90,95,100]
sns.histplot(data = 1e18*geo['throat.volume'], stat='percent', bins=20, color='Blue', legend=True, cbar=True, cbar_kws=dict(shrink=.75), label = 'Throat Volume')
plt.ylim(0,100)
plt.xlim(0,1e7)
plt.xlabel('Throat Volume (µm³)')
plt.ylabel('Relative frequency (%)')
plt.legend()

plt.savefig("detailed plot.jpg")

"""# Use of the Stokes algorithm to calculate absolute permeability"""

phase_1 = op.phase.Water(network=pn)
#phase_1.add_model_collection(op.models.collections.phase.air)
phase_1.add_model_collection(op.models.collections.physics.standard)
#mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
#phase_1.add_model(propname='throat.hydraulic_conductance', model=mod)
phase_1.regenerate_models()
inlet = pn.pores('ymin')
outlet = pn.pores('ymax')
flow = op.algorithms.StokesFlow(network=pn, phase=phase_1)
flow.set_value_BC(pores=inlet, values=1)
flow.set_value_BC(pores=outlet, values=0)
flow.run()
phase_1.update(flow.soln)

pi=math.pi
Q = flow.rate(pores=inlet, mode='group')[0]
A = (pi*(im.shape[1]*im.shape[1])/4) * resolution**2
L = im.shape[2] * resolution
mu = phase_1['pore.viscosity'].max()
delta_P = (1 - 0)
K = abs((Q * L * mu) / (A * delta_P))
print(f'The value of K is: {K/0.98e-12:.2f} D')
print(f'The value of K is: {K/0.98e-12*1000:.2f} mD')

"""# Creation of the oil phase"""

oil = op.phase.Phase(network=pn)
oil['pore.viscosity'] = 0.547e-3 #Pa.s
oil['pore.contact_angle'] = 120
oil['pore.surface_tension'] = 0.0483 #N/m
mod = op.models.physics.capillary_pressure.washburn
oil.add_model(propname='throat.entry_pressure', model=mod)
import openpnm.models.physics as pmods
oil.add_model(propname='throat.hydraulic_conductance',
              model=pmods.hydraulic_conductance.hagen_poiseuille)
oil.regenerate_models()

"""# Creation of the water phase"""

water = op.phase.Water(network=pn)
water.add_model_collection(op.models.collections.phase.water)
water.add_model_collection(op.models.collections.physics.basic)
mod = op.models.physics.capillary_pressure.washburn
water.add_model(propname='throat.entry_pressure', model=mod)
water.add_model(propname='throat.hydraulic_conductance',
              model=pmods.hydraulic_conductance.hagen_poiseuille)
water.regenerate_models()


oil['pore.temperature'][0], water['pore.temperature'][0], oil['pore.contact_angle'][0], water['pore.contact_angle'][0]

ip = op.algorithms.InvasionPercolation(network=pn, phase=oil)
Finlets_init = pn.pores('ymin')
Finlets=([Finlets_init[x] for x in range(0, len(Finlets_init), 2)])
ip.set_inlet_BC(pores=Finlets)
ip.run()

def sat_occ_update(network, nwp, wp, ip, i):
    r"""
        Calculates the saturation of each phase using the invasion
        sequence from either invasion percolation.
        Parameters
        ----------
        network: network
        nwp : phase
            non-wetting phase
        wp : phase
            wetting phase
        ip : IP
            invasion percolation (ran before calling this function)
        i: int
            The invasion_sequence limit for masking pores/throats that
            have already been invaded within this limit range. The
            saturation is found by adding the volume of pores and thorats
            that meet this sequence limit divided by the bulk volume.
    """
    pore_mask = ip['pore.invasion_sequence'] < i
    throat_mask = ip['throat.invasion_sequence'] < i
    sat_p = np.sum(network['pore.volume'][pore_mask])
    sat_t = np.sum(network['throat.volume'][throat_mask])
    sat1 = sat_p + sat_t
    bulk = network['pore.volume'].sum() + network['throat.volume'].sum()
    sat = sat1/bulk
    nwp['pore.occupancy'] = pore_mask
    nwp['throat.occupancy'] = throat_mask
    wp['throat.occupancy'] = 1-throat_mask
    wp['pore.occupancy'] = 1-pore_mask
    return sat

def Rate_calc(network, phase, inlet, outlet, conductance):
    phase.regenerate_models()
    St_p = op.algorithms.StokesFlow(network=network, phase=phase)
    St_p.settings._update({'conductance' : conductance})
    St_p.set_value_BC(pores=inlet, values=1)
    St_p.set_value_BC(pores=outlet, values=0)
    St_p.run()
    val = np.abs(St_p.rate(pores=inlet, mode='group'))
    return val

"""# Definition of entry and exit flow"""

flow_in = pn.pores('ymin')
flow_out = pn.pores('ymax')

model_mp_cond = op.models.physics.multiphase.conduit_conductance
oil.add_model(model=model_mp_cond, propname='throat.conduit_hydraulic_conductance',
              throat_conductance='throat.hydraulic_conductance', mode='medium', regen_mode='deferred')
water.add_model(model=model_mp_cond, propname='throat.conduit_hydraulic_conductance',
              throat_conductance='throat.hydraulic_conductance', mode='medium', regen_mode='deferred')

"""# Calculation of relative permeabilities"""

Snwp_num=10
flow_in = pn.pores('ymin')
flow_out = pn.pores('ymax')
max_seq = np.max([np.max(ip['pore.invasion_sequence']),
          np.max(ip['throat.invasion_sequence'])])
start = max_seq//Snwp_num
stop = max_seq
step = max_seq//Snwp_num
Snwparr = []
relperm_nwp = []
relperm_wp = []

for i in range(start, stop, step):
    oil.regenerate_models();
    water.regenerate_models();
    sat = sat_occ_update(network=pn, nwp=oil, wp=water, ip=ip, i=i)
    Snwparr.append(sat)
    Rate_abs_nwp = Rate_calc(pn, oil, flow_in, flow_out, conductance = 'throat.hydraulic_conductance')
    Rate_abs_wp = Rate_calc(pn, water, flow_in, flow_out, conductance = 'throat.hydraulic_conductance')
    Rate_enwp = Rate_calc(pn, oil, flow_in, flow_out, conductance = 'throat.conduit_hydraulic_conductance')
    Rate_ewp = Rate_calc(pn, water, flow_in, flow_out, conductance = 'throat.conduit_hydraulic_conductance')
    relperm_nwp.append(Rate_enwp/Rate_abs_nwp)
    relperm_wp.append(Rate_ewp/Rate_abs_wp)

data = {'Saturation':np.hstack(Snwparr),'Kr_oil':np.hstack(relperm_nwp),'Kr_water':np.hstack((relperm_wp))}
DF=pd.DataFrame(data).round(6)
DF

plt.figure(figsize=[5,5])
plt.plot(DF['Saturation'], DF['Kr_oil'], '*-', label='Kr_oil', color='red')
plt.plot(DF['Saturation'], DF['Kr_water'], 'o-', label='Kr_water', color='blue')
plt.xlabel('So')
plt.ylabel('Kr')
#plt.xticks(np.arange(0.4, 1.1, 0.1))
plt.title('Relative Permeability Curves')
plt.legend()