1 1
import numpy as np
2 1
import openpnm as op
3 1
from porespy.filters import trim_nonpercolating_paths
4 1
import collections
5

6

7 1
def tortuosity(im, axis, return_im=False, **kwargs):
8
    r"""
9
    Calculates tortuosity of given image in specified direction
10

11
    Parameters
12
    ----------
13
    im : ND-image
14
        The binary image to analyze with ``True`` indicating phase of interest
15
    axis : int
16
        The axis along which to apply boundary conditions
17
    return_im : boolean
18
        If ``True`` then the resulting tuple contains a copy of the input
19
        image with the concentration profile.
20

21
    Returns
22
    -------
23
    results :  tuple
24
        A named-tuple containing:
25
        * ``tortuosity`` calculated using the ``effective_porosity`` as
26
        * ``effective_porosity`` of the image after applying
27
        ``trim_nonpercolating_paths``.  This removes disconnected
28
        voxels which cause singular matrices.
29
        :math:`(D_{AB}/D_{eff}) \cdot \varepsilon`.
30
        * ``original_porosity`` of the image as given
31
        * ``formation_factor`` found as :math:`D_{AB}/D_{eff}`.
32
        * ``image`` containing the concentration values from the simulation.
33
        This is only returned if ``return_im`` is ``True``.
34

35
    """
36 1
    if axis > (im.ndim - 1):
37 0
        raise Exception("Axis argument is too high")
38
    # Obtain original porosity
39 1
    porosity_orig = im.sum()/im.size
40
    # removing floating pores
41 1
    im = trim_nonpercolating_paths(im, inlet_axis=axis, outlet_axis=axis)
42
    # porosity is changed because of trimmimg floating pores
43 1
    porosity_true = im.sum()/im.size
44 1
    if porosity_true < porosity_orig:
45 0
        print('Caution, True porosity is:', porosity_true,
46
              'and volume fraction filled:',
47
              abs(porosity_orig-porosity_true)*100, '%')
48
    # cubic network generation
49 1
    net = op.network.CubicTemplate(template=im, spacing=1)
50
    # adding phase
51 1
    water = op.phases.Water(network=net)
52 1
    water['throat.diffusive_conductance'] = 1  # dummy value
53
    # running Fickian Diffusion
54 1
    fd = op.algorithms.FickianDiffusion(network=net, phase=water)
55
    # choosing axis of concentration gradient
56 1
    inlets = net['pore.coords'][:, axis] <= 1
57 1
    outlets = net['pore.coords'][:, axis] >= im.shape[axis]-1
58
    # boundary conditions on concentration
59 1
    C_in = 1.0
60 1
    C_out = 0.0
61 1
    fd.set_value_BC(pores=inlets, values=C_in)
62 1
    fd.set_value_BC(pores=outlets, values=C_out)
63
    # Use specified solver if given
64 1
    if 'solver_family' in kwargs.keys():
65 1
        fd.settings.update(kwargs)
66
    else:  # Use pyamg otherwise, if presnet
67 1
        try:
68 1
            import pyamg
69 0
            fd.settings['solver_family'] = 'pyamg'
70 1
        except ModuleNotFoundError:  # Use scipy cg as last resort
71 1
            fd.settings['solver_family'] = 'scipy'
72 1
            fd.settings['solver_type'] = 'cg'
73 1
    op.utils.tic()
74 1
    fd.run()
75 1
    op.utils.toc()
76
    # calculating molar flow rate, effective diffusivity and tortuosity
77 1
    rate_out = fd.rate(pores=outlets)[0]
78 1
    rate_in = fd.rate(pores=inlets)[0]
79 1
    if not np.allclose(-rate_out, rate_in):
80 0
        raise Exception('Something went wrong, inlet and outlet rate do not match')
81 1
    delta_C = C_in - C_out
82 1
    L = im.shape[axis]
83 1
    A = np.prod(im.shape)/L
84 1
    N_A = A/L*delta_C
85 1
    Deff = rate_in/N_A
86 1
    tau = porosity_true/(Deff)
87 1
    result = collections.namedtuple('tortuosity_result', ['tortuosity',
88
                                                          'effective_porosity',
89
                                                          'original_porosity',
90
                                                          'formation_factor',
91
                                                          'image'])
92 1
    result.tortuosity = tau
93 1
    result.formation_factor = 1/Deff
94 1
    result.original_porosity = porosity_orig
95 1
    result.effective_porosity = porosity_true
96 1
    if return_im:
97 0
        conc = np.zeros([im.size, ], dtype=float)
98 0
        conc[net['pore.template_indices']] = fd['pore.concentration']
99 0
        conc = np.reshape(conc, newshape=im.shape)
100 0
        result.image = conc
101
    else:
102 1
        result.image = None
103 1
    return result

Read our documentation on viewing source code .

Loading