diff --git a/src/openpnm/core/_base2.py b/src/openpnm/core/_base2.py index f98adc99db..ff52a8340c 100644 --- a/src/openpnm/core/_base2.py +++ b/src/openpnm/core/_base2.py @@ -526,7 +526,10 @@ def get_conduit_data(self, propname): Parameters ---------- propname : str - The dictionary key of the property to fetch. + The dictionary key of the property to fetch. Can be 'pore.diameter' + or just 'diameter'. To fetch different propnames for pores and throats + supply them both in a list, like + `['pore.diameter', 'throat.inscribed_diameter']`. Returns ------- @@ -535,8 +538,13 @@ def get_conduit_data(self, propname): for pore1, throat, and pore2 respectively. """ - poreprop = 'pore.' + propname.split('.', 1)[-1] - throatprop = 'throat.' + propname.split('.', 1)[-1] + if isinstance(propname, str): + poreprop = 'pore.' + propname.split('.', 1)[-1] + throatprop = 'throat.' + propname.split('.', 1)[-1] + elif isinstance(propname, list): + poreprop = propname[0] if propname[0].startswith('pore') else propname[1] + throatprop = \ + propname[0] if propname[0].startswith('throat') else propname[1] conns = self.network.conns try: T = self[throatprop] diff --git a/src/openpnm/models/geometry/conduit_lengths/_funcs.py b/src/openpnm/models/geometry/conduit_lengths/_funcs.py index d9db7423c3..4e90f5037e 100644 --- a/src/openpnm/models/geometry/conduit_lengths/_funcs.py +++ b/src/openpnm/models/geometry/conduit_lengths/_funcs.py @@ -1,7 +1,7 @@ import numpy as np - from openpnm.models.geometry import _geodocs + __all__ = [ "spheres_and_cylinders", "circles_and_rectangles", @@ -29,8 +29,6 @@ def spheres_and_cylinders( Calculates conduit lengths in the network assuming pores are spheres and throats are cylinders. - A conduit is defined as ( 1/2 pore - full throat - 1/2 pore ). - Parameters ---------- %(network)s @@ -44,9 +42,16 @@ def spheres_and_cylinders( of the pore-throat-pore conduits. The array is formatted as ``[pore1, throat, pore2]``. + Notes + ----- + The cylindrical throat is assumed to overlap the spherical pore bodies creating + a hemispherical cap at the area of intersection. The length of the cap is + assigned to the throat. + """ L_ctc = _get_L_ctc(network) - D1, Dt, D2 = network.get_conduit_data(pore_diameter.split(".", 1)[-1]).T + D1, Dt, D2 = \ + network.get_conduit_data(propname=[pore_diameter, throat_diameter]).T # Handle the case where Dt > Dp if (Dt > D1).any() or (Dt > D2).any(): @@ -81,7 +86,9 @@ def spheres_and_cylinders( @_geodocs def circles_and_rectangles( - network, pore_diameter="pore.diameter", throat_diameter="throat.diameter" + network, + pore_diameter="pore.diameter", + throat_diameter="throat.diameter", ): r""" Calculates conduit lengths in the network assuming pores are circles @@ -107,13 +114,14 @@ def circles_and_rectangles( return spheres_and_cylinders( network=network, pore_diameter=pore_diameter, - throat_diameter=throat_diameter, ) @_geodocs def cones_and_cylinders( - network, pore_diameter="pore.diameter", throat_diameter="throat.diameter" + network, + pore_diameter="pore.diameter", + throat_diameter="throat.diameter", ): r""" Calculates conduit lengths in the network assuming pores are cones @@ -132,7 +140,8 @@ def cones_and_cylinders( """ L_ctc = _get_L_ctc(network) - D1, Dt, D2 = network.get_conduit_data(pore_diameter.split(".", 1)[-1]).T + D1, Dt, D2 = \ + network.get_conduit_data(propname=[pore_diameter, throat_diameter]).T L1 = D1 / 2 L2 = D2 / 2 @@ -147,7 +156,11 @@ def cones_and_cylinders( @_geodocs -def intersecting_cones(network, pore_coords="pore.coords", throat_coords="throat.coords"): +def intersecting_cones( + network, + pore_coords="pore.coords", + throat_coords="throat.coords", +): r""" Calculates conduit lengths in the network assuming pores are cones that intersect. Therefore, the throat is the cross sectional plane where @@ -176,7 +189,10 @@ def intersecting_cones(network, pore_coords="pore.coords", throat_coords="throat @_geodocs def hybrid_cones_and_cylinders( - network, pore_diameter="pore.diameter", throat_coords="throat.coords" + network, + pore_diameter="pore.diameter", + throat_diameter="throat.diameter", + throat_coords="throat.coords", ): r""" Calculates conduit lengths in the network assuming pores are cones @@ -188,6 +204,7 @@ def hybrid_cones_and_cylinders( ---------- %(network)s %(Dp)s + %(Dt)s %(Tcoords)s Returns @@ -195,7 +212,8 @@ def hybrid_cones_and_cylinders( """ L_ctc = _get_L_ctc(network) - D1, Dt, D2 = network.get_conduit_data(pore_diameter.split(".", 1)[-1]).T + D1, Dt, D2 = \ + network.get_conduit_data(propname=[pore_diameter, throat_diameter]).T L1 = D1 / 2 L2 = D2 / 2 @@ -216,7 +234,9 @@ def hybrid_cones_and_cylinders( @_geodocs def trapezoids_and_rectangles( - network, pore_diameter="pore.diameter", throat_diameter="throat.diameter" + network, + pore_diameter="pore.diameter", + throat_diameter="throat.diameter", ): r""" Calculates conduit lengths in the network assuming pores are @@ -246,7 +266,9 @@ def trapezoids_and_rectangles( @_geodocs def intersecting_trapezoids( - network, pore_coords="pore.coords", throat_coords="throat.coords" + network, + pore_coords="pore.coords", + throat_coords="throat.coords", ): r""" Calculates conduit lengths in the network assuming pores are @@ -276,7 +298,9 @@ def intersecting_trapezoids( @_geodocs def hybrid_trapezoids_and_rectangles( - network, pore_diameter="pore.diameter", throat_coords="throat.coords" + network, + pore_diameter="pore.diameter", + throat_coords="throat.coords", ): r""" Calculates conduit lengths in the network assuming pores are @@ -306,7 +330,9 @@ def hybrid_trapezoids_and_rectangles( @_geodocs def pyramids_and_cuboids( - network, pore_diameter="pore.diameter", throat_diameter="throat.diameter" + network, + pore_diameter="pore.diameter", + throat_diameter="throat.diameter", ): r""" Calculates conduit lengths in the network assuming pores are truncated @@ -331,7 +357,9 @@ def pyramids_and_cuboids( @_geodocs def intersecting_pyramids( - network, pore_coords="pore.coords", throat_coords="throat.coords" + network, + pore_coords="pore.coords", + throat_coords="throat.coords", ): r""" Calculates conduit lengths in the network assuming pores are @@ -356,7 +384,9 @@ def intersecting_pyramids( @_geodocs def hybrid_pyramids_and_cuboids( - network, pore_diameter="pore.diameter", throat_coords="throat.coords" + network, + pore_diameter="pore.diameter", + throat_coords="throat.coords", ): r""" Calculates conduit lengths in the network assuming pores are truncated @@ -382,7 +412,9 @@ def hybrid_pyramids_and_cuboids( @_geodocs def cubes_and_cuboids( - network, pore_diameter="pore.diameter", throat_diameter="throat.diameter" + network, + pore_diameter="pore.diameter", + throat_diameter="throat.diameter", ): r""" Calculates conduit lengths in the network assuming pores are cubes @@ -399,7 +431,8 @@ def cubes_and_cuboids( """ L_ctc = _get_L_ctc(network) - D1, Dt, D2 = network.get_conduit_data(pore_diameter.split(".", 1)[-1]).T + D1, Dt, D2 = \ + network.get_conduit_data(propname=[pore_diameter, throat_diameter]).T L1 = D1 / 2 L2 = D2 / 2 @@ -416,7 +449,9 @@ def cubes_and_cuboids( @_geodocs def squares_and_rectangles( - network, pore_diameter="pore.diameter", throat_diameter="throat.diameter" + network, + pore_diameter="pore.diameter", + throat_diameter="throat.diameter", ): r""" Calculates conduit lengths in the network assuming pores are squares diff --git a/tests/unit/models/physics/DiffusiveConductanceTest.py b/tests/unit/models/physics/DiffusiveConductanceTest.py index 5a114bc0db..660c857daa 100644 --- a/tests/unit/models/physics/DiffusiveConductanceTest.py +++ b/tests/unit/models/physics/DiffusiveConductanceTest.py @@ -120,6 +120,142 @@ def test_taylor_aris_diffusion(self): desired = np.array([0.121193, 0.126131, 0.118578]) assert_allclose(actual, desired, rtol=1e-5) + def test_spheres_and_cylinders(self): + pn = op.network.Cubic(shape=[2, 1, 1]) + Dp1, Dp2 = 0.5, 0.25 + Dt = 0.1 + pn['pore.diameter'] = [Dp1, Dp2] + pn['throat.diameter'] = Dt + + f = op.models.geometry.diffusive_size_factors.spheres_and_cylinders + pn.add_model(propname='throat.diffusive_size_factors', model=f) + Scalc = pn['throat.diffusive_size_factors'].flatten() + # Pore 1 + Lp1 = (Dp1/2)*np.cos(np.arcsin(Dt/Dp1)) + Sp1 = 1/(np.pi*(Dp1/2))*np.arctanh(Lp1/(Dp1/2)) + assert_allclose(1/Sp1, Scalc[0], rtol=1e-7) + # Pore 2 + Lp2 = (Dp2/2)*np.cos(np.arcsin(Dt/Dp2)) + Sp2 = 1/(np.pi*(Dp2/2))*np.arctanh(Lp2/(Dp2/2)) + assert_allclose(1/Sp2, Scalc[2], rtol=1e-7) + # Throat + Lt = 1 - Lp1 - Lp2 + St = Lt/(np.pi*(Dt/2)**2) + assert_allclose(1/St, Scalc[1], rtol=1e-6) + # Finally check conductance + f = op.models.physics.diffusive_conductance.generic_diffusive + phase = op.phase.Phase(network=pn) + phase['pore.diffusivity'] = 1.0 + phase['throat.diffusivity'] = 1.0 + phase.add_model(propname='throat.diffusive_conductance', model=f) + Gcalc = phase['throat.diffusive_conductance'][0] + G = (Sp1 + St + Sp2)**-1 + assert_allclose(G, Gcalc, rtol=1e-7) + + def test_spheres_and_cylinders_oversized(self): + pn = op.network.Cubic(shape=[2, 1, 1]) + Dp1, Dp2 = 1.75, 1.25 + Dt = 0.1 + pn['pore.diameter'] = [Dp1, Dp2] + pn['throat.diameter'] = Dt + + f = op.models.geometry.diffusive_size_factors.spheres_and_cylinders + pn.add_model(propname='throat.diffusive_size_factors', model=f) + Scalc = pn['throat.diffusive_size_factors'].flatten() + # Pore 1 + Lp1 = (Dp1/2)*np.cos(np.arcsin(Dt/Dp1)) + Sp1 = 1/(np.pi*(Dp1/2))*np.arctanh(Lp1/(Dp1/2)) + assert_allclose(1/Sp1, Scalc[0], rtol=1e-7) + # Pore 2 + Lp2 = (Dp2/2)*np.cos(np.arcsin(Dt/Dp2)) + Sp2 = 1/(np.pi*(Dp2/2))*np.arctanh(Lp2/(Dp2/2)) + assert_allclose(1/Sp2, Scalc[2], rtol=1e-7) + # Throat + Lt = 1 - Lp1 - Lp2 + St = Lt/(np.pi*(Dt/2)**2) + assert_allclose(1/St, Scalc[1], rtol=1e-6) + # Finally check conductance + f = op.models.physics.diffusive_conductance.generic_diffusive + phase = op.phase.Phase(network=pn) + phase['pore.diffusivity'] = 1.0 + phase['throat.diffusivity'] = 1.0 + phase.add_model(propname='throat.diffusive_conductance', model=f) + Gcalc = phase['throat.diffusive_conductance'][0] + G = (Sp1 + St + Sp2)**-1 + assert_allclose(G, Gcalc, rtol=1e-7) + + def test_cones_and_cylinders(self): + pn = op.network.Cubic(shape=[2, 1, 1]) + Dp1, Dp2 = 0.5, 0.25 + Dt = 0.1 + pn['pore.diameter'] = [Dp1, Dp2] + pn['throat.diameter'] = Dt + + f = op.models.geometry.diffusive_size_factors.cones_and_cylinders + pn.add_model(propname='throat.diffusive_size_factors', model=f) + Scalc = pn['throat.diffusive_size_factors'].flatten() + # Pore 1 + Sp1 = (Dp1/2)/(np.pi*(Dp1/2*Dt/2)) + assert_allclose(1/Sp1, Scalc[0], rtol=1e-7) + # Pore 2 + Sp2 = (Dp2/2)/(np.pi*(Dp2/2*Dt/2)) + assert_allclose(1/Sp2, Scalc[2], rtol=1e-7) + # Throat + Lt = 1 - Dp1/2 - Dp2/2 + St = Lt/(np.pi*(Dt/2)**2) + assert_allclose(1/St, Scalc[1], rtol=1e-6) + # Finally check conductance + f = op.models.physics.diffusive_conductance.generic_diffusive + phase = op.phase.Phase(network=pn) + phase['pore.diffusivity'] = 1.0 + phase['throat.diffusivity'] = 1.0 + phase.add_model(propname='throat.diffusive_conductance', model=f) + Gcalc = phase['throat.diffusive_conductance'][0] + G = (Sp1 + St + Sp2)**-1 + assert_allclose(G, Gcalc, rtol=1e-7) + + def test_intersecting_cones(self): + pn = op.network.Cubic(shape=[2, 1, 1]) + Dp1, Dp2 = 0.5, 0.25 + Dt = 0.1 + pn['pore.diameter'] = [Dp1, Dp2] + pn['throat.diameter'] = Dt + pn['throat.coords'] = [[1.0, 0.5, 0.5]] + + f = op.models.geometry.diffusive_size_factors.intersecting_cones + pn.add_model(propname='throat.diffusive_size_factors', model=f) + Scalc = pn['throat.diffusive_size_factors'].flatten() + # Pore 1 + Lp1 = 0.5 + Sp1 = (Lp1)/(np.pi*(Dp1/2*Dt/2)) + assert_allclose(1/Sp1, Scalc[0], rtol=1e-7) + # Pore 2 + Lp2 = 0.5 + Sp2 = (Lp2)/(np.pi*(Dp2/2*Dt/2)) + assert_allclose(1/Sp2, Scalc[2], rtol=1e-7) + # Throat + St = 0 + # Finally check conductance + f = op.models.physics.diffusive_conductance.generic_diffusive + phase = op.phase.Phase(network=pn) + phase['pore.diffusivity'] = 1.0 + phase['throat.diffusivity'] = 1.0 + phase.add_model(propname='throat.diffusive_conductance', model=f) + Gcalc = phase['throat.diffusive_conductance'][0] + G = (Sp1 + St + Sp2)**-1 + assert_allclose(G, Gcalc, rtol=1e-7) + + + + + + + + + + + + if __name__ == '__main__':