Merge pull request #2853 from PMEAL/more_robust_conduit_lengths

More robust conduit length models #bug

openpnm/models/geometry/conduit_lengths/_funcs.py

@@ -51,15 +51,31 @@ def spheres_and_cylinders(
     # Handle the case where Dt > Dp
     if (Dt > D1).any() or (Dt > D2).any():
         _raise_incompatible_data()
+
+    # If spheres do not overlap:
     L1 = np.sqrt(D1**2 - Dt**2) / 2
     L2 = np.sqrt(D2**2 - Dt**2) / 2
+    Lt = L_ctc - (L1 + L2)
+
+    if np.any(Lt < 0):  # Find pores that touch/overlap
+        # Find diameter of overlap between spheres
+        d = L_ctc
+        R1 = D1/2
+        R2 = D2/2
+        # Check distance to the intersection
+        L1_int = (d**2 - R2**2 + R1**2) / (2*d)
+        if np.any(L1_int < 0) or np.any(L1_int > L_ctc):
+            raise Exception('The pores overlap too much')
+        D_int = 2/(2*d) * np.sqrt(4*d**2 * R1**2 - (d**2 - R2**2 + R1**2)**2)
+        mask = D_int > Dt
+        if np.any(mask):
+            d = d[mask]
+            R1 = R1[mask]
+            R2 = R2[mask]
+            L1[mask] = (d**2 - R2**2 + R1**2) / (2*d)
+            L2[mask] = L_ctc[mask] - L1[mask]
+            Lt[mask] = 1e-15
 
-    # Handle throats w/ overlapping pores
-    _L1 = (4 * L_ctc**2 + D1**2 - D2**2) / (8 * L_ctc)
-    mask = L_ctc - 0.5 * (D1 + D2) < 0
-    L1[mask] = _L1[mask]
-    L2[mask] = (L_ctc - L1)[mask]
-    Lt = np.maximum(L_ctc - (L1 + L2), 1e-15)
     return np.vstack((L1, Lt, L2)).T
 
 

tests/integration/Flow_in_straight_conduit.py

@@ -0,0 +1,55 @@
+# This test checks to ensure that pores which overlap give the same transport
+# rates as pores which do not but are the same size
+if __name__ == "__main__":
+
+    import openpnm as op
+    import numpy as np
+
+    mu = 1.0
+    Pin = 1.0
+    Pout = 0.0
+    diameters = [0.9, 1.0, 1.1, 1.5]
+    Q = np.pi*(np.array(diameters)/2)**4/(8*mu*4) * (Pin - Pout)
+
+    for i, D in enumerate(diameters):
+        pn5 = op.network.Cubic(shape=[5, 1, 1], spacing=1)
+        pn5['pore.diameter'] = D
+        pn5['throat.diameter'] = D
+        pn5.regenerate_models()
+        pn5.add_model(
+            propname='throat.hydraulic_size_factors',
+            model=op.models.geometry.hydraulic_size_factors.spheres_and_cylinders)
+
+        pn3 = op.network.Cubic(shape=[3, 1, 1], spacing=2)
+        pn3['pore.diameter'] = D
+        pn3['throat.diameter'] = D
+        pn3.add_model(
+            propname='throat.hydraulic_size_factors',
+            model=op.models.geometry.hydraulic_size_factors.spheres_and_cylinders)
+        pn3.regenerate_models()
+
+        ph1 = op.phase.Phase(network=pn3)
+        ph1['pore.viscosity'] = mu
+        ph1.add_model(
+            propname='throat.hydraulic_conductance',
+            model=op.models.physics.hydraulic_conductance.generic_hydraulic)
+
+        ph2 = op.phase.Phase(network=pn5)
+        ph2['pore.viscosity'] = mu
+        ph2.add_model(
+            propname='throat.hydraulic_conductance',
+            model=op.models.physics.hydraulic_conductance.generic_hydraulic)
+
+        sf1 = op.algorithms.StokesFlow(network=pn3, phase=ph1)
+        sf1.set_value_BC(pores=pn3.pores('xmin'), values=Pin)
+        sf1.set_value_BC(pores=pn3.pores('xmax'), values=Pout)
+        sf1.run()
+        r1 = sf1.rate(pores=pn3.pores('xmin'))
+
+        sf2 = op.algorithms.StokesFlow(network=pn5, phase=ph1)
+        sf2.set_value_BC(pores=pn5.pores('xmin'), values=Pin)
+        sf2.set_value_BC(pores=pn5.pores('xmax'), values=Pout)
+        sf2.run()
+        r2 = sf2.rate(pores=pn5.pores('xmin'))
+        assert np.allclose(r1, r2, atol=0, rtol=1e-10)
+        assert np.allclose(r1, Q[i], atol=0, rtol=1e-10)

tests/unit/models/geometry/ConduitDiffusiveSizeFactorsTest.py

@@ -21,10 +21,9 @@ def test_spheres_and_cylinders(self):
                     'throat': SF[:, 1],
                     'pore2': SF[:, 2]}
         S_desired = {
-            'pore1': array([0.93875728, 0.97459088]),
-            'throat': array([1.25663706e+14, 4.05813214e-01]),
-            'pore2': array([0.82884551, 0.94886745])
-        }
+            'pore1': array([1.06932839, 0.97459088]),
+            'throat': array([4.02746504, 0.40581321]),
+            'pore2': array([0.97459088, 0.94886745])}
         for k, v in S_actual.items():
             assert_allclose(v, S_desired[k])
 
@@ -34,10 +33,9 @@ def test_circles_and_rectangles(self):
                     'throat': SF[:, 1],
                     'pore2': SF[:, 2]}
         S_desired = {
-            'pore1': array([1.53390798, 1.80140852]),
-            'throat': array([4.00000000e+14, 1.29174358e+00]),
-            'pore2': array([1.65097536, 2.07781253])
-        }
+            'pore1': array([1.62474893, 1.80140852]),
+            'throat': array([12.81981939,  1.29174358]),
+            'pore2': array([1.80140852, 2.07781253])}
         for k, v in S_actual.items():
             assert_allclose(v, S_desired[k])
 
@@ -112,10 +110,9 @@ def test_ncylinders_in_series(self):
                     'throat': SF[:, 1],
                     'pore2': SF[:, 2]}
         S_desired = {
-            'pore1': array([0.49260457, 0.77443401]),
-            'throat': array([1.25663706e+14, 4.05813214e-01]),
-            'pore2': array([0.563723, 0.84600371])
-        }
+            'pore1': array([0.69688454, 0.77443401]),
+            'throat': array([4.02746504, 0.40581321]),
+            'pore2': array([0.77443401, 0.84600371])}
         for k, v in S_actual.items():
             assert_allclose(v, S_desired[k])
 

tests/unit/models/geometry/ConduitHydraulicSizeFactorsTest.py

@@ -21,10 +21,9 @@ def test_spheres_and_cylinders(self):
                     'throat': SF[:, 1],
                     'pore2': SF[:, 2]}
         S_desired = {
-            'pore1': array([0.01068901, 0.01195694]),
-            'throat': array([6.28318531e+11, 2.02906607e-03]),
-            'pore2': array([0.00771764, 0.00917032])
-        }
+            'pore1': array([0.01655401, 0.01195694]),
+            'throat': array([0.02013733, 0.00202907]),
+            'pore2': array([0.01195694, 0.00917032])}
         for k, v in S_actual.items():
             assert_allclose(v, S_desired[k], rtol=1e-5)
 
@@ -34,10 +33,9 @@ def test_circles_and_rectangles(self):
                     'throat': SF[:, 1],
                     'pore2': SF[:, 2]}
         S_desired = {
-            'pore1': array([0.06563123, 0.06697876]),
-            'throat': array([5.33333333e+12, 1.72232477e-02]),
-            'pore2': array([0.05072058, 0.05686537])
-        }
+            'pore1': array([0.08485281, 0.06697876]),
+            'throat': array([0.17093093, 0.01722325]),
+            'pore2': array([0.06697876, 0.05686537])}
         for k, v in S_actual.items():
             assert_allclose(v, S_desired[k], rtol=1e-5)
 
@@ -112,10 +110,9 @@ def test_ncylinders_in_series(self):
                     'throat': SF[:, 1],
                     'pore2': SF[:, 2]}
         S_desired = {
-            'pore1': array([0.0020471, 0.00612218]),
-            'throat': array([6.28318531e+11, 2.02906607e-03]),
-            'pore2': array([0.00261812, 0.00661558])
-        }
+            'pore1': array([0.00506726, 0.00612218]),
+            'throat': array([0.02013733, 0.00202907]),
+            'pore2': array([0.00612218, 0.00661558])}
         for k, v in S_actual.items():
             assert_allclose(v, S_desired[k], rtol=1e-5)
 

tests/unit/models/geometry/ConduitLengthsTest.py

@@ -24,12 +24,23 @@ def test_spheres_and_cylinders(self):
         with pytest.raises(Exception):
             L_actual = mods.spheres_and_cylinders(self.net)
         self.net["pore.diameter"][1] = 0.9
-        # Overlapping pores
+        # Slightly Overlapping pores
         self.net["pore.diameter"][0] = 1.2
         L_actual = mods.spheres_and_cylinders(self.net)
+        L_desired = np.array([[0.56568542, 0.03120169, 0.40311289],
+                              [0.40311289, 0.30965898, 0.28722813]])
+        assert_allclose(L_actual, L_desired)
+        # Moderately overlapping pores
+        self.net["throat.diameter"][0] = 0.2
+        L_actual = mods.spheres_and_cylinders(self.net)
         L_desired = np.array([[5.78750000e-01, 1.00000000e-15, 4.21250000e-01],
                               [4.03112887e-01, 3.09658980e-01, 2.87228132e-01]])
         assert_allclose(L_actual, L_desired)
+        self.net["throat.diameter"][0] = 0.4
+        # Strongly overlapping pores
+        self.net["pore.diameter"][0] = 2.5
+        with pytest.raises(Exception):
+            mods.spheres_and_cylinders(self.net)
         self.net["pore.diameter"][0] = 0.5
 
     def test_circles_and_rectangles(self):
@@ -40,14 +51,25 @@ def test_circles_and_rectangles(self):
         # Incompatible data with model assumptions
         self.net["pore.diameter"][1] = 0.3
         with pytest.raises(Exception):
-            L_actual = mods.circles_and_squares(self.net)
+            L_actual = mods.spheres_and_cylinders(self.net)
         self.net["pore.diameter"][1] = 0.9
-        # Overlapping pores
+        # Slightly Overlapping pores
         self.net["pore.diameter"][0] = 1.2
-        L_actual = mods.circles_and_rectangles(self.net)
+        L_actual = mods.spheres_and_cylinders(self.net)
+        L_desired = np.array([[0.56568542, 0.03120169, 0.40311289],
+                              [0.40311289, 0.30965898, 0.28722813]])
+        assert_allclose(L_actual, L_desired)
+        # Moderately overlapping pores
+        self.net["throat.diameter"][0] = 0.2
+        L_actual = mods.spheres_and_cylinders(self.net)
         L_desired = np.array([[5.78750000e-01, 1.00000000e-15, 4.21250000e-01],
                               [4.03112887e-01, 3.09658980e-01, 2.87228132e-01]])
         assert_allclose(L_actual, L_desired)
+        self.net["throat.diameter"][0] = 0.4
+        # Strongly overlapping pores
+        self.net["pore.diameter"][0] = 2.5
+        with pytest.raises(Exception):
+            mods.spheres_and_cylinders(self.net)
         self.net["pore.diameter"][0] = 0.5
 
     def test_cones_and_cylinders(self):