Repository: PMEAL/OpenPNM
Pull request #1234 (2019-04-09 20:02:41)

classic pore scale model addition 
This PR add classic pore scale model in openpnm v2.1 to reproduce statoil simulations.


Patch:
diff --git a/openpnm/models/physics/diffusive_conductance.py b/openpnm/models/physics/diffusive_conductance.py
index 1f4cd07b9..84b92a2ca 100644
--- a/openpnm/models/physics/diffusive_conductance.py
+++ b/openpnm/models/physics/diffusive_conductance.py
@@ -276,3 +276,72 @@ def generic_conductance(target, transport_type, pore_area, throat_area,
                         ' "diffusion" or "taylor_aris_diffusion"')
     # Apply shape factors and calculate the final conductance
     return (1/gt/SFt + 1/g1/SF1 + 1/g2/SF2)**(-1)
+
+
+def classic_ordinary_diffusion(target,
+                               pore_molar_density='pore.molar_density',
+                               pore_diffusivity='pore.diffusivity',
+                               pore_area='pore.area',
+                               pore_diameter='pore.diameter',
+                               throat_area='throat.area',
+                               throat_length='throat.length',
+                               throat_diameter='throat.diameter',
+                               shape_factor='throat.shape_factor',
+                               **kwargs):
+    r"""
+    Calculate the diffusive conductance of conduits in network, where a
+    conduit is ( 1/2 pore - full throat - 1/2 pore ) based on the areas
+    Parameters
+    ----------
+    network : OpenPNM Network Object
+    phase : OpenPNM Phase Object
+        The phase of interest
+    Notes
+    -----
+    (1) This function requires that all the necessary phase properties already
+    be calculated.
+    (2) This function calculates the specified property for the *entire*
+    network then extracts the values for the appropriate throats at the end.
+    """
+    network = target.project.network
+    throats = network.map_throats(throats=target.Ts, origin=target)
+    phase = target.project.find_phase(target)
+    # Get Nt-by-2 list of pores connected to each throat
+    Ps = network['throat.conns']
+    # Get properties in every pore in the network
+    parea = network[pore_area]
+    pdia = network[pore_diameter]
+    # Get the properties of every throat
+    tdia = network[throat_diameter]
+    tarea = _sp.pi * (tdia / 2) ** 2
+    tlen = network[throat_length]
+    # Interpolate pore phase property values to throats
+    DABt = phase.interpolate_data(propname=pore_diffusivity)[throats]
+    ct = phase.interpolate_data(propname=pore_molar_density)[throats]
+    # Get pore lengths
+    plen1 = (0.5 * pdia[Ps[:, 0]])
+    plen2 = (0.5 * pdia[Ps[:, 1]])
+    # Remove any non-positive lengths
+    plen1[plen1 <= 0] = 1e-12
+    plen2[plen2 <= 0] = 1e-12
+    # Find g for half of pore 1
+    gp1 = ct * DABt * parea[Ps[:, 0]] / plen1
+    gp1[_sp.isnan(gp1)] = _sp.inf
+    gp1[~(gp1 > 0)] = _sp.inf  # Set 0 conductance pores (boundaries) to inf
+    # Find g for half of pore 2
+    gp2 = ct * DABt * parea[Ps[:, 1]] / plen2
+    gp2[_sp.isnan(gp2)] = _sp.inf
+    gp2[~(gp2 > 0)] = _sp.inf  # Set 0 conductance pores (boundaries) to inf
+    # Find g for full throat, remove any non-positive lengths
+    tlen[tlen <= 0] = 1e-12
+    # Get shape factor
+    try:
+        sf = network[shape_factor]
+    except:
+        sf = _sp.ones(network.num_throats())
+    sf[_sp.isnan(sf)] = 1.0
+    gt = (1 / sf) * ct * DABt * tarea / tlen
+    # Set 0 conductance pores (boundaries) to inf
+    gt[~(gt > 0)] = _sp.inf
+    value = (1 / gt + 1 / gp1 + 1 / gp2) ** (-1)
+    return value
diff --git a/openpnm/models/physics/hydraulic_conductance.py b/openpnm/models/physics/hydraulic_conductance.py
index e409b0271..a92d0d787 100644
--- a/openpnm/models/physics/hydraulic_conductance.py
+++ b/openpnm/models/physics/hydraulic_conductance.py
@@ -164,6 +164,70 @@ def hagen_poiseuille_power_law(
                                pore_pressure=pore_pressure)
 
 
+def classic_hagen_poiseuille(target,
+                             pore_diameter='pore.diameter',
+                             pore_viscosity='pore.viscosity',
+                             throat_length='throat.length',
+                             throat_diameter='throat.diameter',
+                             shape_factor='throat.shape_factor',
+                             **kwargs):
+    r"""
+    Calculates the hydraulic conductivity of throat assuming cylindrical
+    geometry using the Hagen-Poiseuille model
+
+    Parameters
+    ----------
+    network : OpenPNM Network Object
+
+    phase : OpenPNM Phase Object
+
+    Notes
+    -----
+    This function calculates the specified property for the *entire* network
+    then extracts the values for the appropriate throats at the end.
+
+    """
+    network = target.project.network
+    throats = network.map_throats(throats=target.Ts, origin=target)
+    # Get Nt-by-2 list of pores connected to each throat
+    Ps = network['throat.conns']
+    # Get properties in every pore in the network
+    phase = target.project.find_phase(target)
+    mup = phase[pore_viscosity]
+    mut = phase.interpolate_data(propname=pore_viscosity)[throats]
+    pdia = network[pore_diameter]
+    # Get pore lengths
+    plen1 = (0.5*pdia[Ps[:, 0]])
+    plen2 = (0.5*pdia[Ps[:, 1]])
+    # Remove any non-positive lengths
+    plen1[plen1 <= 0] = 1e-12
+    plen2[plen2 <= 0] = 1e-12
+    # Find g for half of pore 1
+    gp1 = _sp.pi*(pdia[Ps[:, 0]])**4/(128*plen1*mut)
+    gp1[_sp.isnan(gp1)] = _sp.inf
+    gp1[~(gp1 > 0)] = _sp.inf  # Set 0 conductance pores (boundaries) to inf
+
+    # Find g for half of pore 2
+    gp2 = _sp.pi*(pdia[Ps[:, 1]])**4/(128*plen2*mut)
+    gp2[_sp.isnan(gp2)] = _sp.inf
+    gp2[~(gp2 > 0)] = _sp.inf  # Set 0 conductance pores (boundaries) to inf
+    # Find g for full throat
+    tdia = network[throat_diameter]
+    tlen = network[throat_length]
+    # Remove any non-positive lengths
+    tlen[tlen <= 0] = 1e-12
+    # Get shape factor
+    try:
+        sf = network[shape_factor]
+    except:
+        sf = _sp.ones(network.num_throats())
+    sf[_sp.isnan(sf)] = 1.0
+    gt = (1/sf)*_sp.pi*(tdia)**4/(128*tlen*mut)
+    gt[~(gt > 0)] = _sp.inf  # Set 0 conductance pores (boundaries) to inf
+    value = (1/gt + 1/gp1 + 1/gp2)**(-1)
+    return value
+
+
 def generic_conductance(target, transport_type, pore_area, throat_area,
                         pore_diffusivity, throat_diffusivity,
                         conduit_lengths, conduit_shape_factors, **kwargs):
diff --git a/openpnm/physics/Classic.py b/openpnm/physics/Classic.py
new file mode 100644
index 000000000..7a979fe73
--- /dev/null
+++ b/openpnm/physics/Classic.py
@@ -0,0 +1,54 @@
+from openpnm.utils import Workspace, logging
+from openpnm.physics import GenericPhysics
+from openpnm.models import physics as mods
+logger = logging.getLogger(__name__)
+ws = Workspace()
+
+
+class Classic(GenericPhysics):
+    r"""
+    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.
+
+    """
+
+    def __init__(self, project=None, network=None, phase=None,
+                 geometry=None, settings={}, **kwargs):
+
+        # Define some default settings
+        self.settings.update({'prefix': 'phys'})
+        # Overwrite with user supplied settings, if any
+        self.settings.update(settings)
+
+        # Deal with network or project arguments
+        if network is not None:
+            if project is not None:
+                assert network is project.network
+            else:
+                project = network.project
+
+        super().__init__(project=project, phase=phase, geometry=geometry,
+                         **kwargs)
+
+        self.add_model(propname='throat.hydraulic_conductance',
+                       model=mods.hydraulic_conductance.classic_hagen_poiseuille)
+        self.add_model(propname='throat.diffusive_conductance',
+                       model=mods.diffusive_conductance.classic_ordinary_diffusion)
+        self.add_model(propname='throat.entry_pressure',
+                       model=mods.capillary_pressure.washburn)
diff --git a/openpnm/physics/__init__.py b/openpnm/physics/__init__.py
index 52b600de5..b02ef8871 100644
--- a/openpnm/physics/__init__.py
+++ b/openpnm/physics/__init__.py
@@ -25,3 +25,4 @@ classes that include common pore-scale physical models
 
 from .GenericPhysics import GenericPhysics
 from .Standard import Standard
+from .Classic import Classic


Test patch:
diff --git a/tests/unit/models/physics/DiffusiveConductanceTest.py b/tests/unit/models/physics/DiffusiveConductanceTest.py
index 2ad42b693..b6835d009 100644
--- a/tests/unit/models/physics/DiffusiveConductanceTest.py
+++ b/tests/unit/models/physics/DiffusiveConductanceTest.py
@@ -47,6 +47,18 @@ class DiffusiveConductanceTest:
         actual = sp.around(actual, decimals=5)
         assert_allclose(actual, desired)
 
+    def test_classic_ordinary_diffusion(self):
+        self.geo['pore.diameter'] = 1.0
+        self.geo['throat.diameter'] = 1.0
+        self.geo['throat.length'] = 1e-9
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                              phase=self.air,
+                                              geometry=self.geo)
+        mod = op.models.physics.diffusive_conductance.classic_ordinary_diffusion
+        self.phys.add_model(propname='throat.conductance', model=mod)
+        assert sp.allclose(a=self.phys['throat.conductance'][0],
+                           b=0.00084552)
 
 if __name__ == '__main__':
 
diff --git a/tests/unit/models/physics/HydraulicConductanceTest.py b/tests/unit/models/physics/HydraulicConductanceTest.py
index acc1cdd7e..d4128a97b 100644
--- a/tests/unit/models/physics/HydraulicConductanceTest.py
+++ b/tests/unit/models/physics/HydraulicConductanceTest.py
@@ -1,5 +1,6 @@
 import openpnm as op
 from numpy.testing import assert_approx_equal
+import scipy as _sp
 
 
 class HydraulicConductanceTest:
@@ -30,6 +31,19 @@ class HydraulicConductanceTest:
         actual = self.phys['throat.hydraulic_conductance'].mean()
         assert_approx_equal(actual, desired=1421.0262776)
 
+    def test_classic_hagen_poiseuille(self):
+        self.geo['pore.diameter'] = 1.0
+        self.geo['throat.diameter'] = 1.0
+        self.geo['throat.length'] = 1.0e-9
+        self.air = op.phases.Air(network=self.net)
+        self.phys = op.physics.GenericPhysics(network=self.net,
+                                                   phase=self.air,
+                                                   geometry=self.geo)
+        mod = op.models.physics.hydraulic_conductance.classic_hagen_poiseuille
+        self.phys.add_model(propname='throat.conductance',
+                             model=mod)
+        assert _sp.allclose(a=self.phys['throat.conductance'][0],
+                            b=1330.68207684)
 
 if __name__ == '__main__':