Throat length vs Direct vs Total throat lenght

[JorgeESantos]

Greetings.

I've been making some Stokes simulations on Networks generated with Snow over 3D Binary images and, in the process of tunning my Network to understand how different parameters affect the obtained permeability, I noticed something that I cannot clarify by myself.

Assuming that I am treating my Network pores and throat as "Spheres and Cylinders", the throat.length value to be used in the hydraulic size factors calculation can be achieved by:

Net.add_model(propname='throat.length', model=pnm.models.geometry.throat_length.spheres_and_cylinders)

This value is the same obtained by the method:

Net.add_model(propname='conduits_lenght', model=pnm.models.geometry.conduit_lengths.spheres_and_cylinders)

My question is regarding how is this value calculated, since the available links to the respective scripts do not provide any helpful information. Moreover, and more importantly, why are the obtained values for throat.lenghts or conduit.lenghts, achieved by these models, so different from the ones labeled as "throat.direct_lenght" and "throat.total_lenght"??

https://openpnm.org/_modules/openpnm/models/geometry/conduit_lengths/_funcs.html#spheres_and_cylinders

Best regards, Jorge E. Santos.

[EvgenyUgolkov]

@JorgeESantos Hello. Could you share the example of your code with me? I try to create my own, but all the time i am confused with the old examples of syntax and functionality from the previous 2.8 release on the website, so i have a weak understanding of how we should follow the process and how the functions look like. The working example will help me a lot

[JorgeESantos]

Hello there. Sure thing. Follow my running script below.

Consider the M0 variable a 3D binary numpy array representing the foreground of the porous media which I intend to segment and analyze. I apply the SNOW algorithm to extract the Network, assign the geometric parameters I need to calculate the hydraulic size factors, define a phase (water) with a pore-scale physics model (Hagen-Poisseulle) and I run the flow algorithm.

`from functions import TPMS_engine
import numpy as np

M, voxel_l = TPMS_engine([2,2,2], elem=40, min_dimen=6.5e-3, iT=0.365, meshdata=False)
M0 = np.logical_not(M)*1

from porespy.networks import snow2
import openpnm as pnm
#Getting the border tolerance
snow_net = snow2(M0, boundary_width=0, accuracy='high', voxel_size=voxel_l)
Net = pnm.io.network_from_porespy(snow_net['network'])
Net['pore.diameter'] = Net['pore.inscribed_diameter']
pnm.topotools.find_surface_pores(Net)
border_tol = Net['pore.diameter']/2*Net['pore.surface']
border_tol = np.amax(border_tol[border_tol != 0.0])
xN = border_tol/voxel_l

#region Network for flow
pad_x = 3
snow_net = snow2(M0, boundary_width=[[pad_x, pad_x], 0, 0], accuracy='high', voxel_size=voxel_l)
Net = pnm.io.network_from_porespy(snow_net['network'])

#Network manipulation/tunning
Net['pore.diameter'] = Net['pore.extended_diameter']
#Net['pore.diameter'] = Net['pore.equivalent_diameter']
Net['throat.diameter'] = Net['throat.inscribed_diameter']
Net.add_model(propname='throat.length',
model=pnm.models.geometry.throat_length.spheres_and_cylinders)
Net['pore.volume'] = Net['pore.region_volume']
Net.add_model(propname='throat.lens_volume',
model=pnm.models.geometry.throat_volume.lens)
Net.add_model(propname='throat.total_volume',
model=pnm.models.geometry.throat_volume.cylinder)
Net.add_model(propname='throat.volume',
model=pnm.models.misc.difference,
props=['throat.total_volume', 'throat.lens_volume'])
Net.add_model(propname='pore.cross_sectional_area',
model=pnm.models.geometry.pore_cross_sectional_area.sphere)
Net.add_model(propname='throat.cross_sectional_area',
model=pnm.models.geometry.throat_cross_sectional_area.cylinder)
Net.add_model(propname='throat.hydraulic_size_factors',
model=pnm.models.geometry.hydraulic_size_factors.spheres_and_cylinders)

#region Boundary conditions and labels
max_x = np.amax(Net['pore.coords'][:, 0])
min_x = np.amin(Net['pore.coords'][:, 0])
pores_left = Net['pore.coords'][:, 0] <= min_x + border_tol
pores_right = Net['pore.coords'][:, 0] >= max_x - border_tol
N_pores_left = len(np.argwhere(pores_left == True))
N_pores_right = len(np.argwhere(pores_right == True))
Net.set_label(pores=pores_left, label='left')
Net.set_label(pores=pores_right, label='right')

scaff_Bleft = Net['pore.coords'][:, 0] >= np.amin([i for i in Net['pore.coords'][:, 0]]) + pad_xvoxel_l
scaff_Bright = Net['pore.coords'][:, 0] <= np.amax([i for i in Net['pore.coords'][:, 0]]) - pad_xvoxel_l
scaff_coords_bool = scaff_Bright*scaff_Bleft #mascara booleana das coordenadas do scaffold
Net.set_label(pores=scaff_coords_bool, label='scaffold')
scaff_x_coords = np.where(scaff_coords_bool == True, Net['pore.coords'][:, 0], 0.0)
scaff_x_max = np.amax(scaff_x_coords)
scaff_x_min = np.amin(scaff_x_coords[scaff_x_coords != 0.0])

scaff_pores_right = Net['pore.coords'][:, 0] >= scaff_x_max - border_tol
scaff_pores_left = Net['pore.coords'][:, 0] <= scaff_x_min + border_tol
scaff_pores_left = scaff_pores_left * scaff_coords_bool
scaff_pores_right = scaff_pores_right * scaff_coords_bool
N_scaff_left = len(np.argwhere(scaff_pores_left == True))
N_scaff_right = len(np.argwhere(scaff_pores_right == True))
Net.set_label(pores=scaff_pores_left, label='scaff_left')
Net.set_label(pores=scaff_pores_right, label='scaff_right')
#endregion

#Network health check
Net.add_model(propname='pore.cluster_number',
model=pnm.models.network.cluster_number)
Net.add_model(propname='pore.cluster_size',
model=pnm.models.network.cluster_size)
health = pnm.utils.check_network_health(Net)
print('Network processed')
#endregion

#region Phase and physics definition
water = pnm.phase.Water(network=Net, name='water')
water.add_model_collection(pnm.models.collections.physics.basic)
water['pore.temperature'] = 295.15
water['throat.temperature'] = 295.15 #24ºC = 297.15K
#density
water.add_model(propname='pore.density',
model=pnm.models.phase.density.water_correlation)
#viscosity
water.add_model(propname='pore.viscosity',
model=pnm.models.phase.viscosity.water_correlation)
#surface tension
water.add_model(propname='pore.surface_tension',
model=pnm.models.phase.surface_tension.water_correlation)
water.add_model(propname='throat.diffusive_size_factors',
model=pnm.models.geometry.diffusive_size_factors.spheres_and_cylinders)
water.add_model(propname='throat.hydraulic_conductance',
model=pnm.models.physics.hydraulic_conductance.hagen_poiseuille)
#water.add_model(propname='throat.hydraulic_conductance',
#model=pnm.models.physics.hydraulic_conductance.generic_hydraulic)
water.add_model(propname='throat.capillary_pressure',
model=pnm.models.physics.capillary_pressure.washburn)
print('Phase and physics set')
#endregion

#region Setting Flow Algorithm
Q_in = 1/60*1e-6 #m^3/s
flow = pnm.algorithms.StokesFlow(network=Net, phase=water)
flow.set_rate_BC(pores=Net.pores('left'), rates=Q_in)
flow.set_value_BC(pores=Net.pores('right'), values=0)
#print(flow)
flow.run()
#Permeability calculation
#P_scaffold = Net['pore.scaffold']*flow['pore.pressure']
P_scaffold_left = Net['pore.scaff_left']*flow['pore.pressure']
#P_scaffold_right = Net['pore.scaff_right']*flow['pore.pressure']
dP2 = np.mean(P_scaffold_left[P_scaffold_left != 0.0]) - np.mean(P_scaffold_right[P_scaffold_right != 0.0])

#Q_inlet = flow. Rate(pores=Net.pores('scaff_left'), mode='group')[0]
Q_inlet = flow. Rate(pores=Net.pores('left'), mode='group')[0]
L = pnm.topotools.get_domain_length(Net, inlets=Net.pores('scaff_left'), outlets=Net.pores('scaff_right'))
A = pnm.topotools.get_domain_area(Net, inlets=Net.pores('scaff_left'), outlets=Net.pores('scaff_right'))
Lg = min(M0.shape)*voxel_l
As = Lg**2

u = water['pore. Viscosity'][0]
#K1 = (Q_in * u * L) / (A * dP2)
K2 = (Q_inlet * u * L) / (A * dP2) #full pnm
K4 = (Q_inlet * u * Lg) / (As * dP2)
#experimental comparison
#K = (Q_in * u * L) / (A * P_left_mean)
#K = (Q_inlet * u * L) / (A * P_left_face)

print('Flow simulation performed')`

Hope it helps you out. Regards, Jorge E. Santos.

[EvgenyUgolkov]

@JorgeESantos Many thanks!