Absolute Permeability Calculation Issue (help me please)

[je-ha-jung]

Hello. I am currently facing difficulties in calculating absolute permeability. I have successfully implemented a pore network model using binarized images, but the next step—calculating absolute permeability—has been challenging. Others seem to obtain reasonable results, but I am getting abnormal values. I couldn’t find an example in OpenPNM or PoreSpy to resolve this issue, and it has been quite frustrating. Could you help me identify any steps I might be missing?

Below is my code. I set the top as the inlet and the bottom as the outlet to calculate permeability, adding a sealing process to prevent lateral flow. However, the sides aren’t completely sealed, and a small amount of fluid leaks out. When comparing the inflow and outflow volumes, they should ideally have identical magnitudes with opposite signs, but that’s not the case here, indicating possible leakage. This leakage might be causing the abnormal permeability values I am repeatedly observing. I’m not sure where to start troubleshooting, so I’m seeking assistance. How can I achieve accurate permeability values?

this is my code

import os
import glob
import imageio.v2 as imageio
import porespy as ps
import numpy as np
import openpnm as op
import matplotlib.pyplot as plt
from skimage.segmentation import watershed
from scipy import ndimage as spim
from porespy.tools import randomize_colors
from porespy.filters import find_peaks, trim_saddle_points, trim_nearby_peaks
from PIL import Image
from vtk.util import numpy_support
import vtk

folder_path = r"C:\Users\user\Desktop\lab\포어 네트워크 모델\bentheimer_binary_image_cut_300_tif_invert"

tif_files = glob.glob(os.path.join(folder_path, "*.tif"))
tif_files.sort()
sample_img = Image.open(tif_files[0])
width, height = sample_img.size

stack = np.zeros((len(tif_files), height, width), dtype=np.uint8)

for i, tif_file in enumerate(tif_files):
img = Image.open(tif_file)
stack[i, :, :] = np.array(img)

im = stack == 255

resolution = 0.0581 snow_output = ps.networks.snow2(im, voxel_size=resolution) pn = op.io.network_from_porespy(snow_output.network) print(pn) geo = op.models.collections.geometry.spheres_and_cylinders pn.add_model_collection(geo) pn.regenerate_models()

phase = op.phase.Phase(network=pn)
phase['pore.viscosity'] = 0.001
phase['pore.density'] = 1000
phase.add_model_collection(op.models.collections.physics.basic)
phase.regenerate_models()

percentile_top = 90 inlet = np.where(coords[:, 1] >= np.percentile(coords[:, 1], percentile_top))[0] inlet = inlet[inlet < pn.Np] pn.set_label('top', pores=inlet)

percentile_bottom = 10
outlet = np.where(coords[:, 1] <= np.percentile(coords[:, 1], percentile_bottom))[0]
outlet = outlet[outlet < pn.Np]
pn.set_label('bottom', pores=outlet)

print("Step: Checking network connectivity...")
pore_mask = np.ones(pn.Np, dtype=bool)
clustering = op.topotools.find_clusters(network=pn, mask=pore_mask)

pore_clusters = clustering.pore_labels
largest_cluster_label = np.argmax(np.bincount(pore_clusters))
isolated_pores = np.where(pore_clusters != largest_cluster_label)[0]

if len(isolated_pores) > 0:
print(f"고립된 포어가 {len(isolated_pores)}개 발견되었습니다. 제거합니다.")
op.topotools.trim(network=pn, pores=isolated_pores)
else:
print("고립된 포어가 없습니다.")

flow = op.algorithms.StokesFlow(network=pn, phase=phase) flow.clear_BCs()

flow.set_value_BC(pores=inlet, values=10, mode='overwrite') flow.set_value_BC(pores=outlet, values=0, mode='overwrite')

percentile_side = 5
left_pores = np.where(coords[:, 0] <= np.percentile(coords[:, 0], percentile_side))[0]
right_pores = np.where(coords[:, 0] >= np.percentile(coords[:, 0], 100 - percentile_side))[0]

front_pores = np.where(coords[:, 2] <= np.percentile(coords[:, 2], percentile_side))[0]
back_pores = np.where(coords[:, 2] >= np.percentile(coords[:, 2], 100 - percentile_side))[0]

side_pores = np.concatenate([left_pores, right_pores, front_pores, back_pores])
side_pores = side_pores[side_pores < pn.Np]

for pore in side_pores:
flow.set_rate_BC(pores=[pore], rates=0, mode='overwrite')

try:
flow.run()
phase.update(flow.soln)

inlet_flow_rate = flow.rate(pores=pn.pores('top'), mode='group')[0]
outlet_flow_rate = flow.rate(pores=pn.pores('bottom'), mode='group')[0]
print(f"Flow rate at the inlet (top): {inlet_flow_rate:.4e} m^3/s")
print(f"Flow rate at the outlet (bottom): {outlet_flow_rate:.4e} m^3/s")

if abs(inlet_flow_rate + outlet_flow_rate) > 1e-6:
    print("Warning: Flow imbalance detected. Inlet and outlet flow rates do not match!")

Q = flow.rate(pores=pn.pores('top'), mode='group')[0] A = (300 * 300) * resolution**2 L = 300 * resolution Delta_P = 10 mu = 0.001

K = Q * L * mu / (A * Delta_P)

print(f'The value of K is: {K/0.98e-12*1000:.2f} mD')

except Exception as e:
print(f"Error running the flow algorithm: {e}")

this is my result image