In [1]:
def update_occupancy(res, phase, p_map, t_map):
r'''
Custom function to take the results of the percolation algorithm and update phase occupancy
'''
phase['pore.occupancy'] = 1 - res['pore.occupancy'][p_map]
phase['throat.occupancy'] = 1 - res['throat.occupancy'][t_map]
In [2]:
def _calc_eff_prop(self):
r"""
Returns the main parameters for calculating the effective property
in a linear transport equation. It also checks for the proper
boundary conditions, inlets and outlets.
"""
import numpy as np
if self.settings['quantity'] not in self.keys():
raise Exception('The algorithm has not been run yet. Cannot ' +
'calculate effective property.')
# Determine boundary conditions by analyzing algorithm object
inlets, outlets = self._get_inlets_and_outlets()
Ps = np.isfinite(self['pore.bc_value'])
BCs = np.unique(self['pore.bc_value'][Ps])
Dx = np.abs(np.diff(BCs))
net = self.project.network
# Fetch area and length of domain
[amax, bmax, cmax] = np.max(net['pore.coords'], axis=0)
[amin, bmin, cmin] = np.min(net['pore.coords'], axis=0)
lx = amax-amin
ly = bmax-bmin
lz = cmax-cmin
A = lx*ly
L = lz
flow = self.rate(pores=inlets)
D = np.sum(flow)*L/A/Dx
return D
In [3]:
def bulk_diffusion_wu(physics,
network,
phase,
geometry,
propname,
diffusivity = 'pore.diffusivity',
molar_density = 'pore.molar_density',
throat_diameter = 'throat.diameter',
throat_length = 'throat.length',
pore_diameter = 'pore.diameter',
**params):
r"""
Calculate the diffusive conductance of throats in network (instead of a
conduit) based on the areas
Parameters
----------
network : OpenPNM Network Object
phase : OpenPNM Phase Object
The phase of interest
Notes
-----
This function requires that all the necessary phase properties already be
calculated.
"""
#ct = phase.get_data(prop='molar_density',throats='all',mode='interpolate')
#Interpolate pore values to throats
DABt = phase.interpolate_data(propname='pore.diffusivity')
#Find g for full throat
tdia = network[throat_diameter]
tlen = network[throat_length]
gt = (sp.pi*DABt*tdia**2)/(tlen*4)
g = gt[geometry.throats()]
phase[propname]=g
In [4]:
def simulation(n=8, npts=51):
r'''
Run a percolation and diffusion simulation
Parameters
----------
n : int
number of pores along each direction of the network used for diffusion.
Percolation network has length 2n along the z-axis
npts : int
number of percolation points
'''
import openpnm as op
import numpy as np
print('-' * 80)
print('Running Sim with Net Size =', n, 'num points', npts)
print('-' * 80)
Lc = 25e-6
pn1 = op.network.Cubic(shape=[n,n,2*n], spacing=Lc)
#code to run if boundaries was set to false
Ps = pn1.pores()
Ts = pn1.throats()
geo1 = op.geometry.GenericGeometry(network=pn1, pores=Ps,throats=Ts)
low = .5e-6
high = 9.5e-6
geo1['pore.diameter'] = 24e-6
radii = low + np.random.rand(pn1.num_throats())*(high - low)
geo1['throat.diameter'] = radii*2
geo1.add_model(propname='pore.volume', model=gm.pore_volume.sphere)
geo1.add_model(propname='pore.area', model=gm.pore_area.sphere)
geo1.add_model(propname='throat.length', model=gm.throat_length.classic)
geo1.add_model(propname='throat.volume', model=gm.throat_volume.cylinder)
geo1.add_model(propname='throat.area', model=gm.throat_area.cylinder)
water = op.phases.Water(network = pn1)
Ps = geo1.pores()
Ts = geo1.throats()
phys_water = op.physics.Standard(network=pn1, phase=water, geometry=geo1)
IP_1 = op.algorithms.InvasionPercolation(network=pn1)
IP_1.setup(phase=water)
inlets = pn1.pores('bottom')
IP_1.set_inlets(pores=inlets)
IP_1.run()
z_vals = np.unique(pn1['pore.coords'][:, 2])[int(n/2):int(3*n/2)]
pore_map = np.in1d(pn1['pore.coords'][:, 2], z_vals)
throat_map = pn1.find_neighbor_throats(pn1.pores()[pore_map], mode='intersection')
pn2 = op.network.Cubic(shape=[n,n,n], spacing=Lc)
prj2 = pn2.project
Ps = pn2.pores()
Ts = pn2.throats()
geo2 = op.geometry.GenericGeometry(network=pn2, pores=Ps, throats=Ts)
for prop in geo1.props():
if 'pore' in prop:
geo2[prop] = geo1[prop][pore_map]
else:
geo2[prop] = geo1[prop][throat_map]
air = op.phases.Air(network = pn2)
phys_air = op.physics.Standard(network=pn2, phase=air, geometry=geo2)
del phys_air['throat.diffusive_conductance']
bulk_diffusion_wu(physics = phys_air, network=pn2, phase=air, geometry=geo2, propname='throat.diffusive_conductance')
phys_water.regenerate_models()
phys_air.regenerate_models()
x_values = []
y_values = []
res = IP_1.results(0.0)
update_occupancy(res, air, pore_map, throat_map)
phys_air.add_model(model=op.models.physics.multiphase.conduit_conductance,
propname='throat.conduit_diffusive_conductance',
throat_conductance='throat.diffusive_conductance',
mode='strict')
bulk_diffusivity = air['pore.diffusivity']
bottom_boundary = pn2['pore.bottom']
top_boundary = pn2['pore.top']
for x in range(npts):
IPsat = float(x)/float(npts)
res = IP_1.results(IPsat)
update_occupancy(res, air, pore_map, throat_map)
phys_air.regenerate_models(propnames=['throat.conduit_diffusive_conductance'])
mask = air['throat.conduit_diffusive_conductance'] == air['throat.diffusive_conductance']
am = pn2.create_adjacency_matrix(weights=mask, fmt='coo')
if tt.ispercolating(am,
bottom_boundary,
top_boundary,
'bond'):
print('Step', x, 'is percolating')
Fickian_alg = op.algorithms.FickianDiffusion(network=pn2, phase=air)
Fickian_alg.set_value_BC(values=0.6, pores=top_boundary)
Fickian_alg.set_value_BC(values=0.2, pores=bottom_boundary)
Fickian_alg.setup(conductance='throat.conduit_diffusive_conductance')
Fickian_alg.run()
effective_diffusivity = _calc_eff_prop(Fickian_alg)/np.mean(air['pore.molar_density'])
prj2.purge_object(Fickian_alg)
else:
print('Step', x, 'is NOT percolating')
effective_diffusivity = 1e-12
#calculation of saturation
p_vol = pn2['pore.volume']
saturation = 1 - np.sum(air['pore.occupancy']*p_vol)/np.sum(p_vol)
Deff = (effective_diffusivity/bulk_diffusivity)[0]
print('Step', x, 'Saturation', "%.3f" % saturation, 'Diffusivity', "%.3f" % Deff)
x_values.append(saturation)
y_values.append(Deff)
return x_values, y_values