Repository: PMEAL/OpenPNM
Pull request #1314 (2019-07-25 03:42:41)

Lens volume models
This PR adds two simple models for calculating the volume of lens and pendular rings due to overlap or not of cylindrical throats and spherical pores.  The lens model is meant to be subtracted from the throat volume in the overlapping case (when throat.endpoints are used) and the pendular_rings model is meant to be added in the case when sphere and throats just barely contact each other.  


Patch:
diff --git a/openpnm/models/geometry/throat_volume.py b/openpnm/models/geometry/throat_volume.py
index 5a4953c8e..fa6f22c41 100644
--- a/openpnm/models/geometry/throat_volume.py
+++ b/openpnm/models/geometry/throat_volume.py
@@ -130,3 +130,98 @@ def rectangle(target, throat_length='throat.length',
     intersection of the throat with a spherical pore body.
     """
     return target[throat_length] * target[throat_diameter]
+
+
+def lens(target, throat_diameter='throat.diameter',
+         pore_diameter='pore.diameter'):
+    r"""
+    Calculates the volume residing the hemispherical caps formed by the
+    intersection between cylindrical throats and spherical pores.
+
+    This volume should be subtracted from throat volumes if the throat lengths
+    were found using throat end points.
+
+    Parameters
+    ----------
+    throat_diameter : string
+        The dictionary keys containing the array with the throat diameter
+        values.
+    pore_diameter : string
+        The dictionary keys containing the array with the pore diameter
+        values.
+
+    Returns
+    -------
+    volume : ND-array
+        The volume that should be subtracted from each throat volume to prevent
+        double counting the volume of overlapping area.
+
+    Notes
+    -----
+    This model does NOT consider the possibility that multiple throats might
+    overlap in the same location which could happen if throats are large and
+    connectivity is random.
+
+    See Also
+    --------
+    pendular_ring
+    """
+    network = target.network
+    conns = network['throat.conns']
+    Rp = target[pore_diameter]
+    Rt = target[throat_diameter]
+    a = _sp.atleast_2d(Rt).T
+    q = _sp.arcsin(a/Rp[conns])
+    b = Rp[conns]*_sp.cos(q)
+    h = Rp[conns] - b
+    V = 1/6*_sp.pi*h*(3*a**2 + h**2)
+    return _sp.sum(V, axis=1)
+
+
+def pendular_ring(target, throat_diameter='throat.diameter',
+                  pore_diameter='pore.diameter'):
+    r"""
+    Calculates the volume of the pendular rings residing between the end of
+    a cylindrical throat and spherical pores that are in contact but not
+    overlapping.
+
+    This volume should be added to the throat volume if the throat length was
+    found as the center-to-center distance less the pore radii.
+
+    Parameters
+    ----------
+    throat_diameter : string
+        The dictionary keys containing the array with the throat diameter
+        values.
+    pore_diameter : string
+        The dictionary keys containing the array with the pore diameter
+        values.
+
+    Returns
+    -------
+    volume : ND-array
+        The volume that should be added to each throat volume to account for
+        under-represented void volume at the pore-throat junctions.
+
+    Notes
+    -----
+    This model does NOT consider the possibility that multiple throats might
+    overlap in the same location which could happen if throats are large and
+    connectivity is random.
+
+    See Also
+    --------
+    lens
+    """
+    network = target.network
+    conns = network['throat.conns']
+    Rp = target[pore_diameter]
+    Rt = target[throat_diameter]
+    a = _sp.atleast_2d(Rt).T
+    q = _sp.arcsin(a/Rp[conns])
+    b = Rp[conns]*_sp.cos(q)
+    h = Rp[conns] - b
+    Vlens = 1/6*_sp.pi*h*(3*a**2 + h**2)
+    Vcyl = _sp.pi*(a)**2*h
+    V = Vcyl - Vlens
+    return _sp.sum(V, axis=1)


Test patch:
diff --git a/tests/unit/models/geometry/ThroatVolumeTest.py b/tests/unit/models/geometry/ThroatVolumeTest.py
index 51f6263e1..07aab9b2e 100644
--- a/tests/unit/models/geometry/ThroatVolumeTest.py
+++ b/tests/unit/models/geometry/ThroatVolumeTest.py
@@ -18,6 +18,22 @@ class ThroatVolumeTest:
         self.geo['throat.length'] = 1.0
         self.geo['throat.area'] = 0.03
 
+    def test_lens_and_pendular_ring(self):
+        net = op.network.Cubic(shape=[2, 1, 1])
+        net['pore.diameter'] = 0.5
+        net['throat.diameter'] = 0.25
+        mod = op.models.geometry.throat_volume.lens
+        net.add_model(propname='throat.lens_volume',
+                      model=mod)
+        Vlens = net['throat.lens_volume']
+        assert Vlens == 2*(0.006733852203712552)
+        mod = op.models.geometry.throat_volume.pendular_ring
+        net.add_model(propname='throat.ring_volume',
+                      model=mod)
+        Vcyl = 2*(0.01315292522620208)
+        Vring = net['throat.ring_volume']
+        assert (Vcyl - Vring) == 2*(0.006733852203712552)
+
     def test_cylinder(self):
         self.geo.add_model(propname='throat.volume',
                            model=mods.cylinder,