Permeability

[HFenice]

If I have the following code, how can I compute the permeability of this network?

import openpnm as op
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as spst
import csv
import pandas
import random
df = pandas.read_csv('.\Tomo1.csv')
print(df)
diametro = df['Diametro'].values
frequenza = df['Frequenza'].values
diametro = diametro.reshape((len(diametro), 1)) / diametro.max()
frequenza = frequenza.reshape((len(frequenza), 1))
plt.plot(diametro,frequenza)
plt.show

pn = op.network.Cubic(shape=[25, 25, 1], spacing=1)

fig = op.topotools.plot_connections(network=pn)
_ = plt.axis('equal')

geo = op.geometry.GenericGeometry(network=pn, pores=pn.Ps, throats=pn.Ts)

random.seed(123) 
geo['pore.diameter'] = np.asarray(random.choices(diametro, frequenza, k = pn.Np)).ravel() 

ig = op.topotools.plot_networkx(pn, scale=1.5)

op.io.VTK.save(network=pn, filename='test_file') # Per visualizzare il porenetwork con paraview

f = op.models.geometry.pore_seed.random 
geo.add_model(propname='pore.seed',
              model=f, num_range=[0.0, 0.99])
f = op.models.misc.from_neighbor_pores
geo.add_model(propname='throat.seed',
              model=f, prop='pore.seed', mode='min')
f = op.models.geometry.throat_size.weibull
geo.add_model(propname='throat.diameter',
              model=f, shape=1.5, scale=.2, loc=0)

geo.show_hist('throat.diameter')

# %% Perform Stokes flow to find Permeability coefficient

# Phase definition

water=op.phases.Water(network=pn)
phys_water=op.physics.Standard(network=pn,phase=water, geometry=geo)

[jgostick]

Here is a link to an integration test that I wrote recently: https://github.com/PMEAL/OpenPNM/blob/remove_setup_method_from_algs/tests/integration/PoreSpyIO_on_Berea.py

It's not a fancy Jupyter notebook, but should walk you through the process. I will convert it to a proper notebook someday. If you feel so inclined, perhaps you could turn it into a step-by-step notebook with explanations etc and contribute it back? :-)

[HFenice]

Yes, of course. Here is my code:

# Importing library
import math
import openpnm as op
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as spst
from scipy.optimize import curve_fit
import csv
import pandas
import random
import openpnm.models as mods
from openpnm.geometry import GenericGeometry
from openpnm.io import PerGeos

# Importing data from Avizo

df = pandas.read_csv('13maggio_prova.csv')

print(df)

diametro = df['Diametro'].values 

frequenza = df['Frequenza'].values 

diametro = diametro.reshape((len(diametro), 1)) / diametro.max() 

frequenza = frequenza.reshape((len(frequenza), 1))

plt.plot(diametro,frequenza)

plt.show


# Network topology

'''
num_points (scalar) – 
The number of points to place in the domain, 
which will become the pore centers after the tessellation is performed.

shape (array_like) – The size of the domain.

'''
pn= op.network.DelaunayVoronoiDual(num_points=len(diametro), shape=[5, 5])


op.topotools.plot_connections(network = pn, throats=pn.throats('boundary'))

op.topotools.plot_coordinates(network=pn, pores=pn.pores('boundary'), 
                              markersize=50, c='r')

# Checking network health
h = pn.check_network_health()
print(h)


# My Geometry object by modifing StickAndBall class

class MyGeo(GenericGeometry):
    def __init__(self, random_state = None, **kwargs):
        super().__init__(**kwargs)
        
        random.seed(random_state)

        self['pore.diameter'] = np.asarray(random.choices(diametro,            
                                                          frequenza, 
                                                          k = pn.Np)).ravel()

        self.add_model(propname='pore.area',
                       model=mods.geometry.pore_cross_sectional_area.sphere,
                       pore_diameter='pore.diameter') #Calculates internal surface area of pore bodies assuming 
                                                      #they are spherical then subtracts the area of the neighboring 
                                                      #throats in a crude way, by simply considering the throat cross-sectional area, 
                                                      #thus not accounting for the actual curvature of the intersection

        self.add_model(propname='pore.volume',
                       model=mods.geometry.pore_volume.sphere, #Calculate pore volume 
                                                               #from diameter assuming a 
                                                               #spherical pore body
                       pore_diameter='pore.diameter')
        
        
        self['throat.diameter'] = np.full((pn.Nt,), 1) #Constant

        self.add_model(propname='throat.endpoints',
                       model=mods.geometry.throat_endpoints.spherical_pores,
                       pore_diameter='pore.diameter',
                       throat_diameter='throat.diameter')

        self.add_model(propname='throat.length',
                       model=mods.geometry.throat_length.piecewise, #Calculates throat length 
                                                                    #from end points and optionally a centroid
                       throat_endpoints='throat.endpoints')

        self.add_model(propname='throat.surface_area',
                       model=mods.geometry.throat_surface_area.cylinder, #Calculate surface area for a cylindrical throat
                       throat_diameter='throat.diameter',
                       throat_length='throat.length')

        self.add_model(propname='throat.volume',
                       model=mods.geometry.throat_volume.cylinder, #Calculate throat volume assuing a cylindrical shape
                       throat_diameter='throat.diameter',
                       throat_length='throat.length')

        self.add_model(propname='throat.area',
                       model=mods.geometry.throat_cross_sectional_area.cylinder, #Calculate throat cross-sectional area for a cylindrical throat
                       throat_diameter='throat.diameter')

        self.add_model(propname='throat.conduit_lengths',
                       model=mods.geometry.throat_length.conduit_lengths,  #Calculates conduit lengths. A conduit is defined as half pore + throat + half pore.
                       throat_endpoints='throat.endpoints',
                       throat_length='throat.length')

geo = MyGeo(network=pn, pores=pn.Ps, throats=pn.Ts, random_state = 123) 
geo.show_hist('pore.diameter')
geo.show_hist('throat.diameter')


# Importing data in Avizo by using PerGeos Class
PerGeos.export_data(network=pn, phases=[], filename='prova')

# Perform Stokes flow to find Permeability coefficient
water = op.phases.Water(network=pn) # Phase Object
phys= op.physics.Standard(network=pn, phase=water, geometry=geo)  # Standard is a  Generic class to generate Physics objects

'''
    Parameters
    ----------
    network : OpenPNM Network object
        The network to which this Physics should be attached

    phase : OpenPNM Phase object
        The Phase object to which this Physics applies

    geometry : OpenPNM Geometry object
        The Geometry object that defines the pores/throats where this Physics
        should be applied.

    name : str, optional
        A unique string name to identify the Physics object, typically same as
        instance name but can be anything.  If left blank, and name will be
        generated that include the class name and a random string.
'''

sf_z = op.algorithms.StokesFlow(network=pn, phase=water) #StokesFlow is a subclass of GenericLinearTransport 
                                                         #to simulate viscous flow

'''
Parameters
       inlets : array_like
            The pores where the inlet pressure boundary conditions were
            applied.  If not given an attempt is made to infer them from the
            algorithm.

        outlets : array_like
            The pores where the outlet pressure boundary conditions were
            applied.  If not given an attempt is made to infer them from the
            algorithm.

        domain_area : scalar, optional
            The area of the inlet (and outlet) boundary faces.  If not given
            then an attempt is made to estimate it, but it is usually
            underestimated.

        domain_length : scalar, optional
            The length of the domain between the inlet and outlet boundary
            faces.  If not given then an attempt is made to estimate it, but it
            is usually underestimated.
'''
'''
Calculates the hydraulic conductance of conduits in network.

A conduit is defined as ( 1/2 pore - full throat - 1/2 pore ).

Parameters
----------
    target : _GenericPhysics
        Physics object with which this model is associated.
    pore_area : str
        Dictionary key of the pore area values.
    throat_area : str
        Dictionary key of the throat area values.
    pore_viscosity : str
        Dictionary key of the pore viscosity values.
    throat_viscosity : str
        Dictionary key of the throat viscosity values.
    conduit_lengths : str
        Dictionary key of the conduit length values.
    conduit_shape_factors : str
        Dictionary key of the conduit hydraulic (flow) shape factor values.
        
Notes
-----
This function requires that all the necessary phase properties already
be calculated.

This function calculates the specified property for the *entire*
network then extracts the values for the appropriate throats at the
end.

This function assumes cylindrical throats with constant cross-section
area. Corrections for different shapes and variable cross-section area
can be imposed by passing the proper conduit_shape_factors argument.
'''
        
Pin = 0
Pout = 101325

sf_z.set_value_BC(pores=pn.pores('top'), values=Pin)
sf_z.set_value_BC(pores=pn.pores('bottom'), values=Pout)
sf_z.run()
A = ((np.pi*(5**2))/4) # Area  in micron
L = 5      # Lenght in micron
mu = water['pore.viscosity'][0] # Water Viscosity in Pa*s
Q = sf_z.rate(pores=pn.pores('top'), mode='group') #Calculates the net rate of material moving 
                                                   #into a given set of pores or throats, micron3/s
'''
 Parameters
        ----------
        pores : array_like
            The pores for which the rate should be calculated
        throats : array_like
            The throats through which the rate should be calculated
        mode : string, optional
            Controls how to return the rate. Options are:
            - *'group'*: (default) Returns the cumulative rate of material
            moving into the given set of pores
            - *'single'* : Calculates the rate for each pore individually  
'''                                          
Kz = (Q*L*mu/(A*(Pin - Pout)))*10e-6 #In m/s

print('The permeability coefficient is:', Kz)

# MIP
hg = op.phases.Mercury(network = pn, name = 'mercury') # Mercury
phys_hg = op.physics.GenericPhysics(network=pn, phase=hg, geometry=geo)  
                                                                       
phys_hg.add_model(propname='throat.entry_pressure',
               model=op.models.physics.capillary_pressure.washburn)    
mip = op.algorithms.Porosimetry(network=pn, phase=hg)

mip.set_inlets(pn.pores('front'))  #Set the locations from which the invader enters the network

mip.run()
mip.plot_intrusion_curve()

But I don't understand how DelaunayVoronoiDual Network works. Is it a random newtork? In particular, for me it's not clear how bottom, top, front, back, left and right pores are assigned.