Repository: PMEAL/OpenPNM
Pull request #1199 (2019-01-16 18:41:24)

Statoil update


Patch:
diff --git a/openpnm/core/Base.py b/openpnm/core/Base.py
index 29615dbf4..b8ffdd940 100644
--- a/openpnm/core/Base.py
+++ b/openpnm/core/Base.py
@@ -258,7 +258,8 @@ class Base(dict):
             vals = {}
             keys = self.keys()
             vals.update({k: self.get(k) for k in keys if k.startswith(key + '.')})
-        elif any([k.startswith(key + '.') for k in self.keys(mode='all', deep=True)]):
+        elif any([k.startswith(key + '.') for k in self.keys(mode='all',
+                                                             deep=True)]):
             # Create a subdict of values in subdomains by interleaving
             vals = {}
             keys = self.keys(mode='all', deep=True)
@@ -1134,7 +1135,7 @@ class Base(dict):
             units = [a.units.__str__() for a in arrs if hasattr(a, 'units')]
             if len(units) > 0:
                 if len(set(units)) == 1:
-                        temp_arr *= sp.array([1])*getattr(unyt, units[0])
+                    temp_arr *= sp.array([1])*getattr(unyt, units[0])
                 else:
                     raise Exception('Units on the interleaved array are not equal')
         return temp_arr
@@ -1167,28 +1168,22 @@ class Base(dict):
         array([1.5, 2.5])
 
         """
-        mro = self._mro()
-        if 'GenericNetwork' in mro:
-            net = self
-            Ts = net.throats()
-            Ps = net.pores()
+        boss = self.project.find_full_domain(self)
+        net = self.project.network
+        if boss is self:
+            Ts = boss.throats()
+            Ps = boss.pores()
             label = 'all'
-        elif ('GenericPhase' in mro) or ('GenericAlgorithm' in mro):
-            net = self.project.network
-            Ts = net.throats()
-            Ps = net.pores()
-            label = 'all'
-        elif ('GenericGeometry' in mro) or ('GenericPhysics' in mro):
-            net = self.project.network
-            Ts = net.throats(self.name)
-            Ps = net.pores(self.name)
+        else:
+            Ts = boss.throats(self.name)
+            Ps = boss.pores(self.name)
             label = self.name
         if propname.startswith('throat'):
             # Upcast data to full network size
-            temp = sp.ones((net.Nt,))*sp.nan
+            temp = sp.ones((boss.Nt,))*sp.nan
             temp[Ts] = self[propname]
             data = temp
-            temp = sp.ones((net.Np,))*sp.nan
+            temp = sp.ones((boss.Np,))*sp.nan
             for pore in Ps:
                 neighborTs = net.find_neighbor_throats(pore)
                 neighborTs = net.filter_by_label(throats=neighborTs,
diff --git a/openpnm/core/Subdomain.py b/openpnm/core/Subdomain.py
index 1f5399d96..94d0fb454 100644
--- a/openpnm/core/Subdomain.py
+++ b/openpnm/core/Subdomain.py
@@ -25,8 +25,8 @@ class Subdomain(Base):
     +---------------------+---------------------------------------------------+
 
     The Project object has two methods, ``check_geometry_health`` and
-    ``check_physics_health`` that look to make sure all locations are assigned
-    to one and only one Geometry and/or Physics.
+    ``check_physics_health`` that look to make sure all locations are
+    assigned to one and only one Geometry and/or Physics.
 
     """
 
diff --git a/openpnm/io/Statoil.py b/openpnm/io/Statoil.py
index 5d8b6f0f7..13ee8c48c 100644
--- a/openpnm/io/Statoil.py
+++ b/openpnm/io/Statoil.py
@@ -82,8 +82,14 @@ class Statoil(GenericIO):
                          'throat.length', 'throat.volume',
                          'throat.clay_volume']
         # Add link2 props to net
-        net['throat.length'] = sp.array(link2['throat.length'])
+        cl_t = sp.array(link2['throat.length'])
+        net['throat.length'] = cl_t
+        net['throat.conduit_lengths.throat'] = cl_t
         net['throat.volume'] = sp.array(link2['throat.volume'])
+        cl_p1 = sp.array(link2['throat.pore1_length'])
+        net['throat.conduit_lengths.pore1'] = cl_p1
+        cl_p2 = sp.array(link2['throat.pore2_length'])
+        net['throat.conduit_lengths.pore2'] = cl_p2
         net['throat.clay_volume'] = sp.array(link2['throat.clay_volume'])
         # ---------------------------------------------------------------------
         # Parse the node1 file
@@ -119,6 +125,10 @@ class Statoil(GenericIO):
         net['pore.radius'] = sp.array(node2['pore.radius'])
         net['pore.shape_factor'] = sp.array(node2['pore.shape_factor'])
         net['pore.clay_volume'] = sp.array(node2['pore.clay_volume'])
+        net['throat.area'] = ((net['throat.radius']**2) /
+                              (4.0*net['throat.shape_factor']))
+        net['pore.area'] = ((net['pore.radius']**2) /
+                            (4.0*net['pore.shape_factor']))
 
         if network is None:
             network = GenericNetwork()
diff --git a/openpnm/models/physics/hydraulic_conductance.py b/openpnm/models/physics/hydraulic_conductance.py
index a92d0d787..8e2b965c2 100644
--- a/openpnm/models/physics/hydraulic_conductance.py
+++ b/openpnm/models/physics/hydraulic_conductance.py
@@ -7,6 +7,7 @@ r"""
 """
 
 import scipy as _sp
+import numpy as np
 
 
 def hagen_poiseuille(target,
@@ -164,6 +165,118 @@ def hagen_poiseuille_power_law(
                                pore_pressure=pore_pressure)
 
 
+def valvatne_blunt(target,
+                   pore_viscosity='pore.viscosity',
+                   throat_viscosity='throat.viscosity',
+                   pore_shape_factor='pore.shape_factor',
+                   throat_shape_factor='throat.shape_factor',
+                   pore_area='pore.area',
+                   throat_area='throat.area',
+                   conduit_lengths='throat.conduit_lengths'):
+    r"""
+    Calculate the single phase hydraulic conductance of conduits in network,
+    where a conduit is ( 1/2 pore - full throat - 1/2 pore ) according to [1].
+    Function has been adapted for use with the Statoil imported networks and
+    makes use of the shape factor in these networks to apply Hagen-Poiseuille
+    flow for conduits of different shape classes: Triangular, Square and
+    Circular [2].
+
+    Parameters
+    ----------
+    target : OpenPNM Object
+        The object which this model is associated with. This controls the
+        length of the calculated array, and also provides access to other
+        necessary properties.
+
+    pore_viscosity : string
+        Dictionary key of the pore viscosity values
+
+    throat_viscosity : string
+        Dictionary key of the throat viscosity values
+
+    pore_shape_factor : string
+        Dictionary key of the pore geometric shape factor values
+
+    throat_shape_factor : string
+        Dictionary key of the throat geometric shape factor values
+
+    pore_area : string
+        Dictionary key of the pore area values. The pore area is calculated
+        using following formula:
+            pore_area = (pore_radius ** 2) / (4 * pore_shape_factor)
+        Where theoratical value of pore_shape_factor in circular tube is
+        calculated using following formula:
+            pore_shape_factor = pore_area / perimeter **2 = 1/4π
+
+    throat_area : string
+        Dictionary key of the throat area values. The throat area is calculated
+        using following formula:
+            throat_area = (throat_radius ** 2) / (4 * throat_shape_factor)
+        Where theoratical value of throat_shape_factor in circular tube is
+        calculated using following formula:
+            throat_shape_factor = throat_area / perimeter **2 = 1/4π
+
+    conduit_lengths : string
+        Dictionary key of the throat conduit lengths
+
+    Returns
+    -------
+    g : ND-array
+        Array containing hydraulic conductance values for conduits in the
+        geometry attached to the given physics object.
+
+    References
+    ----------
+    [1] Valvatne, Per H., and Martin J. Blunt. "Predictive pore‐scale modeling
+    of two‐phase flow in mixed wet media." Water Resources Research 40,
+    no. 7 (2004).
+    [2] Patzek, T. W., and D. B. Silin (2001), Shape factor and hydraulic
+    conductance in noncircular capillaries I. One-phase creeping flow,
+    J. Colloid Interface Sci., 236, 295–304.
+    """
+    network = target.project.network
+    mu_p = target[pore_viscosity]
+    try:
+        mu_t = target[throat_viscosity]
+    except KeyError:
+        mu_t = target.interpolate_data(pore_viscosity)
+    # Throat Portion
+    Gt = network[throat_shape_factor]
+    tri = Gt <= np.sqrt(3)/36.0
+    circ = Gt >= 0.07
+    square = ~(tri | circ)
+    ntri = np.sum(tri)
+    nsquare = np.sum(square)
+    ncirc = np.sum(circ)
+    kt = np.ones_like(Gt)
+    kt[tri] = 3.0/5.0
+    kt[square] = 0.5623
+    kt[circ] = 0.5
+    # Pore Portions
+    Gp = network[pore_shape_factor]
+    tri = Gp <= np.sqrt(3)/36.0
+    circ = Gp >= 0.07
+    square = ~(tri | circ)
+    ntri += np.sum(tri)
+    nsquare += np.sum(square)
+    ncirc += np.sum(circ)
+    kp = np.ones_like(Gp)
+    kp[tri] = 3.0/5.0
+    kp[square] = 0.5623
+    kp[circ] = 0.5
+    gp = kp*(network[pore_area]**2)*Gp/mu_p
+    gt = kt*(network[throat_area]**2)*Gt/mu_t
+    conns = network['throat.conns']
+    l1 = network[conduit_lengths + '.pore1']
+    lt = network[conduit_lengths + '.throat']
+    l2 = network[conduit_lengths + '.pore2']
+    # Resistors in Series
+    value = (l1/gp[conns[:, 0]] +
+             lt/gt +
+             l2/gp[conns[:, 1]])
+    return 1/value
+
+
 def classic_hagen_poiseuille(target,
                              pore_diameter='pore.diameter',
                              pore_viscosity='pore.viscosity',
@@ -219,7 +332,7 @@ def classic_hagen_poiseuille(target,
     # Get shape factor
     try:
         sf = network[shape_factor]
-    except:
+    except KeyError:
         sf = _sp.ones(network.num_throats())
     sf[_sp.isnan(sf)] = 1.0
     gt = (1/sf)*_sp.pi*(tdia)**4/(128*tlen*mut)
@@ -266,7 +379,7 @@ def generic_conductance(target, transport_type, pore_area, throat_area,
 
     Returns
     -------
-    g : ndarray
+    g : ND-array
         Array containing conductance values for conduits in the geometry
         attached to the given physics object.
 


Test patch:
diff --git a/tests/unit/models/physics/HydraulicConductanceTest.py b/tests/unit/models/physics/HydraulicConductanceTest.py
index d4128a97b..aa3f23403 100644
--- a/tests/unit/models/physics/HydraulicConductanceTest.py
+++ b/tests/unit/models/physics/HydraulicConductanceTest.py
@@ -1,6 +1,7 @@
 import openpnm as op
 from numpy.testing import assert_approx_equal
 import scipy as _sp
+import numpy as np
 
 
 class HydraulicConductanceTest:
@@ -37,14 +38,24 @@ class HydraulicConductanceTest:
         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)
+                                              phase=self.air,
+                                              geometry=self.geo)
         mod = op.models.physics.hydraulic_conductance.classic_hagen_poiseuille
         self.phys.add_model(propname='throat.conductance',
-                             model=mod)
+                            model=mod)
         assert _sp.allclose(a=self.phys['throat.conductance'][0],
                             b=1330.68207684)
 
+    def test_valvatne_blunte(self):
+        mod = op.models.physics.hydraulic_conductance.valvatne_blunt
+        sf = np.sqrt(3)/36.0
+        self.geo['pore.shape_factor'] = np.ones(self.geo.Np)*sf
+        self.geo['throat.shape_factor'] = np.ones(self.geo.Nt)*sf
+        self.phys.add_model(propname='throat.valvatne_conductance', model=mod)
+        actual = self.phys['throat.valvatne_conductance'].mean()
+        assert_approx_equal(actual, desired=1030.9826)
+
+
 if __name__ == '__main__':
 
     t = HydraulicConductanceTest()