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