Capillary pressure inconsistency between pnextract and openPNM

[RayLee-uoa]

Hi prof,

I am currently working on the capillary pressure prediction of a Bentheimer sandstone sample and have encountered a problem. I tried to use pnextract from Blunt's team and openPNM.
Figure 1 shows the Pc curves generated by the two methods:

We can find the Pc curves obtained from the two methods are quite different!
By comparing to the lab measurement on a similar Bentheimer sandstone sample (figure 2), it can be found that the pnextract's result is more reasonable. The lab measurement is taken from Aghabozorgi, S., Rostami, B. An Investigation of Polymer Adsorption in Porous Media Using Pore Network Modelling. Transp Porous Med 115, 169–187 (2016). https://doi.org/10.1007/s11242-016-0760-5

So I am pretty sure there is something wrong with the way I was processing the openPNM.
below is the code lines of openPNM:

import numpy as np
import openpnm as op
import porespy as ps
import os
import pandas as pd
import matplotlib.pyplot as plt
op.visualization.set_mpl_style()
op.Workspace().settings.loglevel=40
from tabulate import tabulate
from pathlib import Path
np.set_printoptions(precision=6)
np.random.seed(10)

resolution = 18.0e-6
name = 'BS'
nx = 75
ny = 75
nz = 75

print('reading geometry')
dat = np.loadtxt('geo.dat')
im = np.reshape(dat,(nx,ny,nz),order='F')
im = im==0

phi = ps.metrics.porosity(im)

print('start generating pore-network')
net = ps.networks.snow2(phases = im, voxel_size=resolution)
pn = op.io.network_from_porespy(net.network)
pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders)
pn.regenerate_models()
print(pn)

h = op.utils.check_network_health(pn)
op.topotools.trim(network=pn, pores=h['disconnected_pores'])
h = op.utils.check_network_health(pn)
print(h)

hg = op.phase.Mercury(network=pn)
f = op.models.physics.capillary_pressure.washburn
hg.add_model(propname='throat.entry_pressure',
             model=f, 
             surface_tension='throat.surface_tension',
             contact_angle='throat.contact_angle',
             diameter='throat.diameter',)
hg.regenerate_models()

mip = op.algorithms.Drainage(network=pn, phase=hg)
mip.set_inlet_BC(pores=pn.pores(['xmin', 'xmax']))
mip.run(pressures=np.logspace(4, 6))
print(mip)
mip.set_outlet_BC(pores=pn.pores('xmax'), mode='overwrite')
mip.apply_trapping()
data = mip.pc_curve()

pressure = np.array(data.pc)
saturation = np.array(data.snwp)
result_pc = np.column_stack((saturation,pressure))
hed_pc = ['Snw', 'Pc']
hed_pc = np.array(hed_pc)
result_pc = np.vstack((hed_pc, result_pc))
df_pc = pd.DataFrame(result_pc)
df_pc.to_excel('pc_curve.xlsx',index=False,header=False)

Though the size of the test model is 75x75x75 (which is small), bigger block sizes that I tested before faced the same problem. Could you please give me some guidance of how to solve this problem?

Best wishes
Ray

[jgostick]

This is probably due to the way boundary pores are being extracted. snow2 puts pores on the edges in the form of thin regions, 3 voxels thick I think. The size of the 'inscribed sphere' in these pores is only going to be 3 voxels. Depending on how the throat diameter is computed it may be affected. Try doing an extraction with no boundary pores added, then find the 'internal' pores for boundaries using logic like Ps = pn.coords[:, 0] < 1.1*pn.coords[:, 0].min(). Or use the 'extended' diameter instead of the inscribed one. Basically, just make sure the surface throats are correctly sized.