Calculating formation factor using ohmic conduction algorithm

[karaferner]

Hi all,
I am attempting to use OpenPNM to calculate formation factor/effective transport properties through porous samples based on image data. I'm using the snow2 algorithm to generate a pore network through the solid phase of my image data, and then using the Ohmic Conduction algorithm to get effective electrical conductivity. I have followed some previous related examples, but am getting a formation factor value about an order of magnitude lower than what I am expecting and what is typical for my material. Getting a value of 0.054 here, but 0.295 using meshing and COMSOL.

The voxel size of my image data is 6 nm, but I'm not specifying that anywhere. At first I thought this was an issue, but it shouldn't matter since I'm accounting for the length and inlet area in end calculation.

Is there something else I should be specifying in the phase other than the ["throat.electrical_conductance"]? Not really sure what else could be the issue.

Pasting my code below for any help you could provide.
Thanks in advance!

import numpy as np
import porespy as ps
import openpnm as op
import scipy as sp
import matplotlib.pyplot as plt
from PIL import Image
from skimage import io
import imageio
import math

# Read TIF image stack
im_comp = io.imread('sample_image_data.tif')
comp_solid = im_comp==1

inlets = ps.generators.faces(comp_solid.shape, inlet=1)
outlets = ps.generators.faces(comp_solid.shape, outlet=1)

comp_solid = ps.filters.trim_nonpercolating_paths(comp_solid, inlets=inlets, outlets=outlets)

#Create snow network from image data
snow_comp_solid = ps.networks.snow2(comp_solid)
net_comp_solid,geo_comp_solid = op.io.PoreSpy.import_data(snow_comp_solid.network)

#trim pores and label faces
h = net_comp_solid.check_network_health()
op.topotools.trim(net_comp_solid, pores=h['trim_pores'])
op.topotools.label_faces(net_comp_solid)

#create generic phase
phase_comp_solid = op.phases.GenericPhase(network=net_comp_solid)

#set bulk conductivity = 1
sigma_bulk = 1
phase_comp_solid['throat.electrical_conductance'] = sigma_bulk

#ohmic conduction algorithm
oc_comp_solid = op.algorithms.OhmicConduction(network=net_comp_solid, phase=phase_comp_solid)

#label inlet, outlet
inlets_comp_solid = net_comp_solid.pores("front")
outlets_comp_solid = net_comp_solid.pores("back")

# Dirichlet BCs: inlet= 1 V, outlet = 0 V
vL, vR = 1.0, 0.0
oc_comp_solid.set_value_BC(pores=inlets_comp_solid,values=vL)
oc_comp_solid.set_value_BC(pores=outlets_comp_solid, values=vR)

#run ohmic conduction
oc_comp_solid.run()

r_in = oc_comp_solid.rate(pores=inlets_comp_solid)[0]
r_out = oc_comp_solid.rate(pores=outlets_comp_solid)[0]
print(r_in, r_out)

print(f'Rate: {r_in:.5e}')
print()

# length between inlet and outlet faces
L = comp_solid.shape[1]

# area of inlet (and outlet)
A = np.prod(comp_solid.shape) / L

sigma_eff = (r_in) * L /(A*(vL - vR))

# sigma_eff = formation factor *sigma_bulk
kf = sigma_eff/sigma_bulk

print(f"Formation factor: {kf:.5e}")

2D side view of my porous sample

Ohmic conduction algorithm seems to run properly

final calculations

[jgostick]

This is a lot to digest. My first suggestion is that you specify voxel_size when calling snow2. Some of the properties scale non-linearly with size so maybe this is causing issues.

Another thing that can cause problems with network extractions is the boundary pores. The size of some pores on on the boundary can be quite small, resulting in a constriction there. You can try doing the extraction with boundary_width=0 then finding the boundary pores in the network using coordinates like Ps = pn.coords[:, 0] < pn.coords[:, 0].min()*1.1, which will find all pores within 10% of the start of the x-axis, and Ps = pn.coords[:, 0] > pn.coords[:, 0].max()*0.9 for the other axis.

Lastly, the Porespy.import_data function assumes 'pore.inscribed_diameter' for the 'pore.diameter' (I think?). You can try using the extended and equivalent diameters as well, which are both larger, so might explain the low value you're getting.

Hope that enough to get you started.

[ma-sadeghi]

We had a feature request here to implement a function to rescale snow output based on a new voxel size. @karaferner Maybe you could implement it and make a pull request... of course if interested :)

[karaferner]

Thanks for your responses. I tried your suggestions. Specifying voxel_size in snow2 does scale the pore network appropriately, as far as I can tell when I implement op.topotools.plot_coordinates with size_by=net["pore.diameter"]. However, it does not seem to affect the value of the rate calculated by Ohmic_Conduction algorithm. The request you linked @ma-sadeghi is giving a 404 error, but sounds like others have had similar problems. Not sure that I would be the best person to implement this though.

Similarly, it looks like the Porespy.import_data function actually assumes 'pore.extended_diameter' for the 'pore.diameter', which is the largest of the three. At least in my case, 'pore.extended_diameter' > 'pore.equivalent_diameter' > 'pore.inscribed_diameter'. But even so, when I set the pore network diameters to 'pore.equivalent_diameter' or 'pore.inscribed_diameter' (geo["pore.diameter"] = snow.network["pore.inscribed_diameter"]) then run the algorithm again, the value it gave for the rate was the same.

The only thing that did have any affect on the calculated rate was implementing your suggestion of setting boundary_width=0 when calling the snow2 function, then manually setting the inlet and outlet boundary faces to <10% and >90% of the y-axis, respectively. Doing this only increased the formation factor to 0.074 though, still off by an order of magnitude from what I'm expecting.

[karaferner]

Of course, I'm continuing to work through it myself and don't expect either of you to magically solve my problem. Just replying with what I've found after trying your suggestions. Very much appreciate your time and for providing this platform in the first place.

Thinking that it would be a cool demonstration to show agreement of image-based transport property results between your OpenPNM package and more traditional meshing techniques.

[jgostick]

I really need to turn this integration test into a proper notebook. It reads a berea image and computes both tortuosity (aka formation factor) and permeability with the pore network and shows good agreement with known values.

PS - Your questions are all valid and we do like to be helpful as much as possible. I admit that sometimes I lose my patience with questions that are basically due to a lack of python or numpy knowledge/experience. I am honestly sorry if this makes you uncomfortable asking us for help.

[jgostick]

Also, I don't mean to deflect blame, but are you 100% sure the comsol meshing process worked well? The tetrahedral mesh does not match the voxel surfaces, and often will take "shortcuts". Do you have any independent validation of your numbers? If your image is not too large you can use porespy.simulations.tortuosity_fd or TauFactor, or PyTrax, or Puma, etc. These all work on cubic images directly so remove that uncertainty.