Issue of Dt>Dp in extracted network (using PoreSpy's snow2) for computing throat hydraulic conductance

[VictordHaas]

Dear OpenPNM team,

Currently I am rewriting an existing network extraction and network property extraction notebook (e.g. for obtaining porosity, coordination number, permeability) from OpenPNM V2.2.0 to OpenPNM V3.0 (dev). The network is obtained from X-ray u-CT (.tif extension) and details on the network extraction are given at the end of this question. Although the network extraction notebook was easily rewritten to V3.0, I am having some problems with extracting the permeability in all cartesian axes in V3.0. This problem is due to the fact that (some of) the throat diameters (Dt) are larger than the pore diameters (Dp) in the extracted network. At the end I pose two solutions, and I would like to have your feedback on them for which would be the best to implement, and if possible, a solution from your side.

In both versions of the network property extraction notebook, permeability is obtained from the results of the stokes flow algorithm (via Darcy’s law). In OpenPNM V2.2.0 the following snippet displays the initialization of the phase, geometry and algorithm, and running of the latter. Here “net” is the extracted network, saved according to V2.2.0 standards.

The op.physics.Standard module (line 222) makes use of the .hydraulic_conductance.hagen_poiseuille method, which in turn uses the throat_length.conduit_lengths method (line 219), for computing the hydraulic conductance used in the Stokes algorithm. Here no issues for a Dt > Dp is presented, whereby the stokes algorithm can be ran and the permeability can be obtained.

For OpenPNM V3.0, the code is displayed in the following snippet. Here “net” is the extracted network saved according to V3.0 standards.

The physics collection (line 229) makes use of the .hydraulic_conductance.generic_hydraulic method for computing the hydraulic conductance. One could also run the .hydraulic_conductance.hagen_poiseuille method like done in V2.2.0, but this model also runs the generic_hydraulic method to obtain its results. The generic_hydraulic method requires the shape factors, which are computed in op.models.geometry.hydraulic_size_factors.spheres_and_cylinders, which in turn uses the throat.conduit_length method. However, in V3.0 the conduit_length model contains an exception for the case that Dt>Dp inside an if statement:

Subsequently, I commented out the if statement (line 56-57) in op.models.geometry.conduit_lengths.spheres_and_cylinders:

and obtained Nans, due to the negative square root, in L1, L2 and Lt (D1, D2 and Dt contained no Nans), and therefore also in the conduit lengths and in the shape factors.

Possible solutions:
The first solution is by assigning throats the smallest diameter of the pores to which it connects. This will, however, significantly deform the extracted network. To limit this, a mask could be used to select only those throats whose diameters is larger than its neighboring pores (which occurs at 1304/ 38734 throats in the extracted network). I have followed this strategy but it has not worked for me yet. Below the network extraction details, the progress is reported.

The second solution would be to write a conduit lengths function similar to that in OpenPNM V2.2.0, which does not suffer from the fact that Dt>Dp.

Lastly, could I use a functionality in PoreSpy snow2 which circumvents assigning Dt>Dp?

Thank you in advance for your response.

Victor

Network extraction details:
The network is loaded in using skimage.io.imread, inverted using skimage.util.invert and finally converted to a binary image in the boolean from (true’s as void-phase and false’s as solid-phase). The actual network extraction is then done using PoreSpy’s snow2 method. Following the disconnected and isolated pores are trimmed from the network, and it is saved as a .pnm and .vtk file. The following figure shows the .vtk visualization in ParaView (including an overlay image.tif).

Further, assigning either equivalent or inscribed diameter as the pore and throat diameters gave rise to the Dt>Dp.

Edit
I followed (as of now) the first solution strategy: using a mask to select only those throats whose diameters are larger than their neighboring pores, and changing their throat diameters to be equal to that of the smallest connecting pore, whereby I can now compute the throat hydraulic size factors :).

However, I still obtained an error due to the throat hydraulic conductance containing a few Nans (in total 4), and further inspection also indicated that the throat hydraulic conductance contains negative numbers. Both arise from the size factors S1 and S2, computed in .hydraulic_size_factors.spheres_and_cylinders.

Here L1 and L2 contain negative numbers (which seems unphysical) and the Nans are introduced via the arctanh added to b.

[jgostick]

Firstly, I am very sorry for the delayed response here. My github notifications seem to still not working propertly, and/or my email inbox fills up so fast that it just slipped by.

To answer your question, we are familiar with your problem and have implemented other solutions. Specifically, you can use the intersecting cones (or pyramid) models which do care if the throat is bigger. We also have a hybrid version of these two models which uses cones and cylinders except in cases where the cones touch then it switches to the intersection models.

The ultimate problem occurs upstream in the network extraction which I'm guess is on an image with high porosity and some anisotropy rigth? in these cases the act of drawing inscribed spheres inside flattened pore regions results in pores that are as small as the flatter dimension. Coupled with the fact that the throats can be large, this results in the issue you are seeing.