# Pore size

Hi,

I am working on a similar problem as explained in the question . I am trying implement the solution suggested by Robert in #7 by using print all the vertices .
I have compacted spherical particles into a cylinder shape, and want estimates/measure pore size in the compact. I was running the suggestion by print all vertices..
When I am using this I am getting following error message:

"ArgumentError: Python argument types in
FlowEngineT.printVertices()
did not match C++ signature:
"

My question how should implement this code to extract the vertices and type of the void/pore space? What should I use as an input?

Best regards,
Mithu

MWE:

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 15 16:29:55 2021

@author: njb18198
"""

#!/usr/bin/env python
#encoding: ascii

# Testing of the Deformation Enginge with Luding Contact Law
# Modified Oedometric Test
# The reference paper [Haustein2017]
from __future__ import print_function
from yade import utils, plot, timing
import pandas as pd
import numpy as np
from PIL import Image
from scipy.interpolate import interp1d
from csv import writer
import os
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import csv
from matplotlib.pyplot import figure
from pylab import *
from scipy.optimize import curve_fit

o = Omega()
save=save
# Physical parameters
fr = 0.41
rho = 1561
Diameter = 7.9e-5
D=Diameter
r1 = Diameter/2
#r2 = Diameter/2
k1 = 10000
kp = 140000
kc = k1 * 0.1
ks = k1 * 0.1
DeltaPMax = Diameter/3.0
Chi1 = 0.34

if save==0:
g=-9.81
elif save==1:
g=0

o.dt = 1.0e-8
particleMass = (4.0/3.0)*math.pi*r1*r1*r1*rho
r_tab=r_tab*1e-3 #real size
h_tab=h_tab*1e-3
m_tab=m_tab*1e-3
v_tab=math.pi*(r_tab**2)*h_tab
#m_tot=4e-06
m_tot=(v_1mm/v_tab)*m_tab
no_p=m_tot/particleMass
#no_p=8000
Vi1 = math.sqrt(k1/particleMass)*DeltaPMax*Chi1
PhiF1=0.999
#PhiF1 = DeltaPMax*(kp-k1)*(r1+r2)/(kp*2*r1*r2)

Cyl_height=0.005

Comp_press_up= comp_press
Comp_force_up=Comp_press_up*cross_area

Comp_press_lp= comp_press
Comp_force_lp=Comp_press_lp*cross_area

compression_data_save=[]
sc_por_15=2
#sc_por_2=2
#sc_por_1=1
data_to_save=[comp_press/1e6,round(no_p),rho]
compression_data_save.append(data_to_save)

mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-7.5*Diameter,-7.5*Diameter,-30*Diameter),(7.5*Diameter,7.5*Diameter,5.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=round(no_p))
sp.toSimulation(material=mat1)

vtkRecorder = VTKRecorder(fileName='vtkRecorder_'+str(tab_porosity),recorders=['all'])
tab_porosity=tab_porosity
tab_height=tab_height

##Single particle swelling model
def model(r,t,Q_max,rho_t,rho_w,r_0,Diff):
Q=((rho_w*(r**3))/(rho_t*(r_0**3)))-(rho_w/rho_t)+1;
drdt =((Diff*rho_t)/(r*rho_w))*((Q_max-Q)/Q);
return drdt
P=[1.45,rho,1000,396.39e-12]

o.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
Bo1_Wall_Aabb(),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
Ig2_Facet_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom()],
[Ip2_LudingMat_LudingMat_LudingPhys()],
[Law2_ScGeom_LudingPhys_Basic()]
),
NewtonIntegrator(damping=0.1, gravity=[0, 0, g]),
PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
#DeformControl(label="DefControl")
]

def checkForce():
if O.iter < 3000000:
return
if unbalancedForce() > 1:
return
global upper_punch
for i in upper_punch:
body= O.bodies[i]
body.state.vel = (0,0,-0.02)
global lower_punch
for n in lower_punch:
body= O.bodies[n]
body.state.vel = (0,0,0.02)
O.engines = O.engines + [PyRunner(command='storeData()', iterPeriod=25000)]+ [PyRunner(command='saveData()', iterPeriod=100000)]

force_up=0
for i in upper_punch:
body= O.bodies[i]
force_up=force_up+abs(O.forces.f(body.id))
force_lp=0
for n in lower_punch:
body = O.bodies[n]
force_lp = force_lp + abs(O.forces.f(body.id))
if ((force_up > Comp_force_up) and (force_lp > Comp_force_lp)):
for i in upper_punch:
body= O.bodies[i]
body.state.vel = (0,0,0.04)
for n in lower_punch:
body= O.bodies[n]
body.state.vel = (0,0,-0.04)

#force_up=0
#for i in upper_punch:
#body= O.bodies[i]
#force_up=force_up+(O.forces.f(body.id))
force_lp=0
for n in lower_punch:
body = O.bodies[n]
force_lp = force_lp + abs(O.forces.f(body.id))
if force_lp==0:
for i in lower_punch:
body= O.bodies[i]
body.state.vel = (0,0,0)
#if ((force_up==0) and (force_lp==0)):
for i in upper_punch:
body=O.bodies[i]
pos_up=body.state.pos
for i in lower_punch:
body=O.bodies[i]
pos_lp=body.state.pos
if pos_up> pos_lp+utils.aabbDim()+0.0002:
for j in upper_punch: O.bodies.erase(j)
fCheck.command = 'Savecheck()'

def Savecheck():
if save==1:
utils.saveVars('lower_punch',lower_punch=lower_punch)
save_filename='PH101_'+str(tab_porosity)+'_'+str(tab_height)+'mm.xml'
o.save(save_filename)
o.pause()
if save==0:
o.engines = o.engines+[PyRunner(command='ParticleSwelling()', iterPeriod=100000)]

if save==1:
O.run()
# when running with yade-batch, the script must not finish until the simulation is done fully
# this command will wait for that (has no influence in the non-batch mode)
waitIfBatch()
g=-9.81

if save==0:
for b in O.bodies:
if isinstance(b.shape, Sphere):
#b.state.blockedDOFs = 'xyXY'

oldm=b.state.mass
oldI=b.state.inertia

m=oldm*3./4./r
b.state.mass=m

b.state.inertia = 15./16./r*oldI #inertia with respect to x and y axes are not used and the computation here is wrong
b.state.inertia = 15./16./r*oldI #inertia with respect to x and y axes are not used and the computation here is wrong
b.state.inertia = 15./16./r*oldI #only inertia with respect to z axis is usefull
o.dt=1e-6))
for j in lower_punch: O.bodies.erase(j)
for b in O.bodies:
if isinstance(b.shape, Sphere):
continue
else:
O.bodies.erase(b.id)
g=0
FlowEngine.printVertices()

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
 Revision history for this message Robert Caulk (rcaulk) said on 2022-03-21: #1

Hello,

You need to first run a step with FlowEngine included in the O.engines list, then you can use any of the methods associated with that engine with:

flow.printVertices()

(assuming you have labeled it 'flow')

Please note that what you have posted requires us to install PIL, scipy, pandas etc all just to test a very basic functionality of Yade. MWE means you narrow it down to the issue - 98% of what you posted has nothing to do with the issue at hand.

-Robert

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2022-03-21 #2

Hi Robert,

Thank you for your help, I implemented your suggestion and defined as PoreSize in the MWE. However, I get an empty vertices dataset.

I want track the pore size change over time during swelling (radius increase) of the particle. Do I need include a flow in there over time, or can still use it by define dead=1 in the engine? How do change the name of the vertices.txt file so it differ for each iteration?

Best regards,
Mithu

Sorry for the code, here is the new MWE:
from __future__ import print_function
from yade import utils, plot, timing

o = Omega()
o.dt = 1.0e-8

mat1 = O.materials.append(LudingMat(frictionAngle=fr, density=rho, k1=k1, kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0))

# Spheres for compression and walls
sp=pack.SpherePack()
sp.makeCloud((-7.5*Diameter,-7.5*Diameter,-30*Diameter),(7.5*Diameter,7.5*Diameter,5.0*Diameter), rMean=Diameter/2.0,rRelFuzz=0.18,num=round(no_p))
sp.toSimulation(material=mat1)

o.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
Bo1_Wall_Aabb(),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
Ig2_Facet_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom()],
[Ip2_LudingMat_LudingMat_LudingPhys()],
[Law2_ScGeom_LudingPhys_Basic()]
),
NewtonIntegrator(damping=0.1, gravity=[0, 0, g]),
PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
#DeformControl(label="DefControl")
]

def checkForce():
if O.iter < 3000000:
return
if unbalancedForce() > 1:
return

if ((force_up > Comp_force_up) and (force_lp > Comp_force_lp)):
for i in upper_punch:
body= O.bodies[i]
body.state.vel = (0,0,0.04)
for n in lower_punch:
body= O.bodies[n]
body.state.vel = (0,0,-0.04)

if pos_up> pos_lp+utils.aabbDim()+0.0002:
for j in upper_punch: O.bodies.erase(j)
fCheck.command = 'Savecheck()'

def Savecheck():
if save==0:
o.engines = o.engines+[PyRunner(command='PoreSize()', iterPeriod=100000)]

def PoreSize():
flow.printVertices()

if save==0:
for b in O.bodies:
if isinstance(b.shape, Sphere):
#b.state.blockedDOFs = 'xyXY'

oldm=b.state.mass
oldI=b.state.inertia

m=oldm*3./4./r
b.state.mass=m

b.state.inertia = 15./16./r*oldI #inertia with respect to x and y axes are not used and the computation here is wrong
b.state.inertia = 15./16./r*oldI #inertia with respect to x and y axes are not used and the computation here is wrong
b.state.inertia = 15./16./r*oldI #only inertia with respect to z axis is usefull
o.dt=1e-6

 Revision history for this message Robert Caulk (rcaulk) said on 2022-03-22: #3

Hello,

>>can still use it by define dead=1.However, I get an empty vertices dataset.

If you want FlowEngine to triangulate your packing, then it should not be dead. If FlowEngine is dead, then it is not doing anything, i.e. not computing a triangulation, hence your empty dataset. If you do not want to compute fluid flow but you still want a triangulation, then you should look into using TesselationWrapper .

>>pore size change over time during swelling (radius increase) of the particle.

printVertices() is simply giving you the vertices of the triangulation, it is not giving you any pore sizes. You would need to post process the vertices yourself to get the pore sizes.

I see you are using "o.engines" where you should be using "O.engines" (and plenty other places where you use o.* instead of O.*).

Cheers,

Robert

 Revision history for this message Jan Stránský (honzik) said on 2022-03-22: #4

@ Robert
> I see you are using "o.engines" where you should be using "O.engines"

I noticed this, too. In the code, this structure is used:
o = Omega()

not very "mainstream", but works..

Cheers
Jan

 Revision history for this message Robert Caulk (rcaulk) said on 2022-03-22: #5

Ah voilà ;)

 Revision history for this message Mithushan Soundaranathan (mithushan93) said on 2022-03-22: #6

Hi Robert,

Sorry for asking, how can obtain the pore size using tesselationwrapper. Do I use volume function in the this engine than substract the sphere volume. Than how does I account for overlapping sphere? Could you please point me in the right direction.

Best regards,
Mithu

 Revision history for this message Robert Caulk (rcaulk) said on 2022-03-23: #7

Hello,

TesselationWrapper can "saveState" which appears to print the vertex positions . Here is some example usage for microstrain .

Cheers,

Robert