2D Geometry

[masoodmoghaddam]

After finishing hierarchical porous network, its time to create a suitable geometry. We have two choice: StickAndBall GenericGeometry We had some problem in this section some of them were resulted from network. Now network is in better situation and some problem of geometry generation should be removed. Are there still remaining problems?

[masoodmoghaddam]

I think we have not any consideration about the throats so far. We consider sphere pore as a extruded (in to the paper) cylinder. But what about the throat between to pores? In 3D throat was a pipe with diameter equals to half of minimum diameter of two pores. In 2D I think, its a channel with cross section area like this: height=half of cylinder height (= 0.5 m in this case) width= half of minimum diameter of two cylinder We had not any consideration about height (and half of minimum in width)

[masoodmoghaddam]

I found this in line with our discussion: https://github.com/PMEAL/OpenPNM/issues/740 Please check it

[masoodmoghaddam]

this is another one: https://github.com/PMEAL/OpenPNM/issues/337

[masoodmoghaddam]

In both of these there is a dissociation about 2 network. Jeff said some thing about dual network. Check it out please, is there any example?

[samanjalilian]

I completed the code and 2 Gaussian distribution are defined for pore diameters( merged and not-merged ) but unfortunately the results are not similar to fig 5 of the article .

import openpnm as op import random import numpy as np from math import sqrt import scipy as sp import matplotlib.pyplot as plt from openpnm.utils import PrintableDict, logging, Workspace ws = Workspace() logger = logging.getLogger(name)

------------------------------------------------ Functions --------------------------------------

def Merge_pores(network, pores, Nei,Cordinate, labels=['merged']):

Assert that pores is list of lists

try:
    len(pores[0])
except (TypeError, IndexError):
    pores = [pores]        
N = len(pores)
NBs, XYZs = [], []
for Ps in pores:
    NBs.append(Nei)
    XYZs.append(Cordinate)
op.topotools.extend(network, pore_coords=XYZs, labels=labels)
Pnew = network.Ps[-N::]
pores_set = [set(items) for items in pores]
NBs_set = [set(items) for items in NBs]
ps1, ps2 = [], []    
from itertools import combinations
for i, j in combinations(range(N), 2):
    if not NBs_set[i].isdisjoint(pores_set[j]):
        ps1.append([network.Ps[-N+i]])
        ps2.append([network.Ps[-N+j]])
# Add (possible) connections between the new pores
op.topotools.connect_pores(network, pores1=ps1, pores2=ps2, labels=labels)
# Add connections between the new pores and the rest of the network
op.topotools.connect_pores(network, pores2=sp.split(Pnew, N), pores1=NBs, labels=labels)
# Trim merged pores from the network
op.topotools.trim(network=network, pores=sp.concatenate(pores))
Ps = pn.find_nearby_pores(pores=network.Ps[-len(pores)::], r=PtoP, flatten=True)
op.topotools.trim(network=network, pores=Ps)

def MacroProsity(prosity): Nnp=int((1-prosity)Nni) radiuos=PSRaverageporeaize /1.5#true value for radius is /2 Ps = pn.find_nearby_pores(pores=random.choice(pn.pores(labels=['surface'], mode='not')), r=PtoP2, flatten=True) op.topotools.merge_pores(network=pn, pores=Ps, labels=['merged']) while len(pn.pores(labels=['merged'], mode='not'))>Nnp: Point=random.choice(pn.pores(labels=['surface','merged'], mode='not')) pores = pn.find_nearby_pores(pores=Point, r=radiuos, flatten=True) pores=[i for i in pores if i not in pn.pores(labels=['surface','merged'])] if len(pores)>0: neighbor=pn.find_neighbor_pores(pores=pores, mode='union', flatten=True, include_input=False) cord=pn['pore.coords'][neighbor].mean(axis=0) Merge_pores(network=pn, pores=pores, Nei=neighbor,Cordinate=cord, labels=['merged']) D=(1- len(pn.pores(labels=['merged'], mode='not'))/Nni)/prosity50 print('[','='int(D), '>','.'(50-int(D)),']')

Deleting Pores which placed very close (less than initial pore to pore distance of the grid) to merged porse

op.topotools.trim(network=pn, pores=pn.find_nearby_pores(pores=pn.pores(labels=['merged']), r=PtoP, flatten=True))
# Deleting pores with just one throat
op.topotools.trim(network=pn, pores=[sp.where(pn.num_neighbors(pn.pores()) == 1)])

------------------------------------------------ contorol parameters --------------------------------------

dp=40e-6 #particle diameter PtoP=10e-8 #pore to pore espacing My idea =dp/Numberofpores Numberofpores=int(dp/PtoP) #number of pores on radius phi=0.01 PSR=20 # pore size ratio averageporeaize=1.6e-8 #average pore size R=8.314 radius=PSR*averageporeaize/1.5

------------------------------------------------ circular pore pn generation --------------------------

pn = op.network.Cubic(shape=[Numberofpores, Numberofpores, 1], spacing=PtoP,connectivity=6) middlepoint=pn.Np/2+Numberofpores/2 Ps = pn.find_nearby_pores(pores=(middlepoint), r=dp/2, flatten=True) pn['pore.dummy_1']=[False] * pn.Np pn['pore.dummy_1'][Ps]=True pn['pore.dummy_1'][int(middlepoint)]=True Pss=pn.pores('pore.dummy_1',mode='not')

------------------------------------------------ finding boundary pores --------------------------------------

Pb = pn.find_nearby_pores(pores=(middlepoint), r=(dp/2-PtoP), flatten=True) pn['pore.surface']=[True] * pn.Np pn['pore.surface'][Pb]=False pn['pore.surface'][int(middlepoint)]=False

------------------------------------------------ triming forcircular pore pn generation ---------------------

op.topotools.trim(network=pn, pores=Pss) Nni=pn.Np

------------------------------------------------ Hierarchical pn generation --------------------------------

MacroProsity(prosity=phi)

------------------------------------------------- Geometry(generic geometry) --------------------------------

geom = op.geometry.GenericGeometry(network=pn,pores=pn.Ps, throats=pn.Ts) nanoDistrb=np.random.normal(averageporeaize, averageporeaize/8, len(pn.pores('merged',mode='not'))) macroDistrb=np.random.normal(nanoDistrb.mean()PSR,(radius2-PSRaverageporeaize)/3 , len(pn.pores('merged')))
geom['pore.diameter']=nanoDistrb.mean() geom['pore.diameter'][pn.pores('merged',mode='not')] = nanoDistrb geom['pore.diameter'][pn.pores('merged')]=macroDistrb geom.add_model(propname='throat.diameter', model=op.models.geometry.throat_size.from_neighbor_pores,regen_mode='normal', mode='min') geom['pore.area']=geom['pore.diameter'] #cross section area of a pore geom['pore.volume']=3.14((geom['pore.diameter'])*2)/4 geom.add_model(propname='throat.endpoints', model=op.models.geometry.throat_endpoints.circular_pores, pore_diameter='pore.diameter',throat_diameter='throat.diameter') geom.add_model(propname='throat.conduit_lengths', model=op.models.geometry.throat_length.conduit_lengths) geom['throat.area']=geom['throat.diameter'] geom['throat.length']=geom['throat.conduit_lengths.throat'] geom['throat.volume']=geom['throat.conduit_lengths.throat']geom['throat.diameter'] geom['throat.surface_area']=2geom['throat.length'] geom['pore.surface_area']=3.14geom['pore.diameter'] #surface area of a pore

--------------------------------------------------- Phase --------------------------------------------------

gas = op.phases.GenericPhase(network=pn) gas['pore.diffusivity']=1.76e-5 gas['pore.molecular_weight']=0.01604 gas['pore.temperature'] = 298

--------------------------------------------------- Knudsen diffusion -------------------------------------------

Dab=gas['pore.diffusivity'] Dkn=(pn['pore.diameter']/3)(sqrt((8R300)/(3.14gas['pore.molecular_weight'][1]))) Dtot=1/((1/Dkn)+(1/Dab)) gas['pore.diffusivity']=Dtot throats=pn.map_throats(throats=pn.Ts, origin=pn) gas['throat.diffusivity']=gas.interpolate_data(propname='pore.diffusivity')[throats]

--------------------------------------------------- Physics--------------------------------------------------------

phys_geom = op.physics.GenericPhysics(network=pn, phase=gas , geometry=geom) mod = op.models.physics.diffusive_conductance.ordinary_diffusion phys_geom.add_model(propname='throat.diffusive_conductance', model=mod)

--------------------------------------------------- Damkhler -------------------------------------------------------

Get diffusive conductance values

diff_cond =gas['throat.diffusive_conductance']

Get incidence matrix with diffusive_conductance values as weights

inc_mat = pn.create_incidence_matrix(weights=diff_cond)

Sum the incidence matrix across y-axis

sum_of_diff_cond = inc_mat.sum(axis=1)

Get the number of throats for each pore

num_throats_per_pore = pn.create_incidence_matrix().sum(axis=1)

Calculate average diffusive_conductance

diff_cond_avg = sum_of_diff_cond / num_throats_per_pore

Finally, calculate Damkohler number for each pore

rxn_rate = 1 * pn['pore.surface_area'] # assuming 1st order kinetics Da = rxn_rate /np.array(diff_cond_avg).flatten()

--------------------------------------------------- Source term -------------------------------------------------

j=-5 averageDamkhler=2*10*(j)
Kjeff=averageDamkhler/np.average(Da) gas['pore.concentration']=0 phys_geom['pore.item1'] =-Kjeffgeom['pore.surface_area'] phys_geom['pore.item2'] = 0 phys_geom.add_model(model=op.models.physics.generic_source_term.linear, propname='pore.Reaction', A1='pore.item1', A2='pore.item2', X='pore.concentration')

--------------------------------------------------- Algorithm -------------------------------------------------

alg = op.algorithms.FickianDiffusion(network=pn,phase=gas) alg.setup(phase=gas, quantity='pore.concentration', conductance='throat.diffusive_conductance') alg.set_value_BC(pores=pn.pores('surface'),values=1) alg.set_source(propname='pore.Reaction',pores=pn.pores(labels=['surface','merged'], mode='not')) alg.run() phys_geom.regenerate_models(propnames='pore.Reaction') gas.update(alg.results()) X = gas['pore.concentration'] rate1=phys_geom['pore.Reaction.rate'][geom.pores(labels=['surface','merged'], mode='not')].sum() rate2=alg.rate(pores=pn.pores(labels=['surface','merged'],mode='not'), mode='group') rate3=np.round(np.sum(Kjeffgeom['pore.surface_area'][geom.pores(labels=['surface','merged'], mode='not')] X[pn.pores(labels=['surface','merged'], mode='not')]), 20) print ('phi=' , phi, 'Average pore Damköhler number= 2e' +str(j),'\n', ' rate1=', rate1, '\n',
' rate2=', rate2,'\n', ' rate3=', rate3)

[samanjalilian]

Rate calculated from the code above for phi=0.01, dp=40, PSR=20 and averageDa=2e-5 is 9.53652988907e-07 which is about e-6 times bigger than the rate reported in fig 5

[masoodmoghaddam]

Hi Dr@abbasabbassiaut