yade-users team mailing list archive

[Ricardo Lorenzoni <question696302@xxxxxxxxxxxxxxxxxxxxx>]

Question #696302 on Yade changed:
https://answers.launchpad.net/yade/+question/696302

    Status: Needs information => Open

Ricardo Lorenzoni gave more information on the question:
Hi. Thanks for your answer.

figure 1 - https://pasteboard.co/JV5qiMz.png
In this first image, you can see the experimental results, where the curve shows that as more weight is added to the sample, the effect on the compacting of the grain mass will decrease.

figure 2 - https://pasteboard.co/JV5mPSf.png
In this second image, the results of the computer simulations are plotted, in which, we can see that the effect of increasing weight on the compaction of the particles is linear.

I have already carried out a series of tests with several contact and
physical models, in none of them I was able to replicate the curve
behavior, as shown in Figure 1.

Here is an example of the base code I am using for the simulations.
the variable "altura" indicates the height of the grain mass on the sample analyzed and, is set to 10 m. 

sample code with cundallStrack contact model:

import sys, math, csv, os
from minieigen import *
from yade import pack, plot, utils, ymport, export

matWall=CohFrictMat(density=2700,young=7e10,poisson=0.334,label="parede")
idparede=O.materials.append(matWall)

r_mean=0.002963
altura=10


lx=0
hx=0.1
ly=0
hy=0.1
lz=0
hz=0.15

p1 = (lx, ly, lz)
p2 = (lx, hy, lz)
p3 = (hx, ly, lz)
p4 = (hx, hy, lz)

p5 = (lx, ly, hz)
p6 = (lx, hy, hz)
p7 = (hx, ly, hz)
p8 = (hx, hy, hz)

p9 = (lx, ly, hz+0.05)
p10= (lx, hy, hz+0.05)
p11= (hx, ly, hz+0.05)
p12= (hx, hy, hz+0.05)


# DEFINICAO DAS FACES

#base
base1=utils.facet([p1,p2,p3],mask=1,wire=True,material=idparede)
base2=utils.facet([p2,p3,p4],mask=1,wire=True,material=idparede)
O.bodies.append(base1)
O.bodies.append(base2)

#Wall1
parede11=utils.facet([p1,p2,p5],mask=1,wire=True,material=idparede)
parede12=utils.facet([p2,p5,p6],mask=1,wire=True,material=idparede)
O.bodies.append(parede11)
O.bodies.append(parede12)

#Wall2
parede21=utils.facet([p1,p3,p5],mask=1,wire=True,material=idparede)
parede22=utils.facet([p3,p5,p7],mask=1,wire=True,material=idparede)
O.bodies.append(parede21)
O.bodies.append(parede22)

#Wall3
parede31=utils.facet([p3,p4,p7],mask=1,wire=True,material=idparede)
parede32=utils.facet([p4,p7,p8],mask=1,wire=True,material=idparede)
O.bodies.append(parede31)
O.bodies.append(parede32)

#Wall4
parede41=utils.facet([p2,p4,p6],mask=1,wire=True,material=idparede)
parede42=utils.facet([p4,p6,p8],mask=1,wire=True,material=idparede)
O.bodies.append(parede41)
O.bodies.append(parede42)

#Wall5
parede51=utils.facet([p5,p6,p9],mask=1,wire=True,material=idparede)
parede52=utils.facet([p6,p9,p10],mask=1,wire=True,material=idparede)
O.bodies.append(parede51)
O.bodies.append(parede52)

#Wall6
parede61=utils.facet([p6,p8,p12],mask=1,wire=True,material=idparede)
parede62=utils.facet([p12,p10,p6],mask=1,wire=True,material=idparede)
O.bodies.append(parede61)
O.bodies.append(parede62)

#Wall7
parede71=utils.facet([p7,p8,p11],mask=1,wire=True,material=idparede)
parede72=utils.facet([p8,p11,p12],mask=1,wire=True,material=idparede)
O.bodies.append(parede71)
O.bodies.append(parede72)

#Wall8
parede81=utils.facet([p5,p7,p9],mask=1,wire=True,material=idparede)
parede82=utils.facet([p7,p9,p11],mask=1,wire=True,material=idparede)
O.bodies.append(parede81)
O.bodies.append(parede82)

particles=CohFrictMat(density=1163,young=2.6e6,poisson=0.25,alphaKr=0.05,frictionAngle=radians(22.6),label="particle")
idparticle=O.materials.append(particles)

O.bodies.append(sphere((0.01,0.01,0.01),radius=r_mean,material=idparticle))

# DEFINICAO DOS MOTORES
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()],verletDist=.05*.29e-3),
   InteractionLoop(
      # handle sphere+sphere and facet+sphere collisions
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()] #modelo linear de contato
   ),
   NewtonIntegrator(gravity=(0,0,-9.81),damping=0.04),
   BoxFactory(maxParticles=-1, extents=(hx/2,hy/2,0.01), center=(hx/2,hy/2,hz+0.03), rMin=r_mean, rMax=r_mean, vMin=100, vMax=100, vAngle=0, massFlowRate=0.01*16.6, normal=(0,0,0), label='factory', materialId=idparticle),
   PyRunner(virtPeriod=0.02, command='verificaAltura()'),
]

global pesomassa
pesomassa=0
def calculaPeso():
	global pesomassa
	for b in O.bodies:
		if isinstance(b.shape, Sphere):
			pesomassa+=b.state.mass	
	pesometro=(altura*((100*pesomassa)/15))
	adicionaPeso(pesometro)
	
global executa
executa=True

def verificaAltura():
	global executa
	if executa:	
		count=0
		for b in O.bodies:
			if isinstance(b.shape, Sphere):
				if (b.state.vel[2]<0.05 and b.state.pos[2]>0.15 and b.state.pos[2]<0.16):
					count+=1
		if (count>16.6*4):
			executa=False
			O.engines[4].dead=True
			ajusta()
		

def ajusta():
	bodiesToBeDeleted=[]
	for b in O.bodies:
		if b.state.pos[2]>hz or b.state.pos[1]<ly or b.state.pos[1]>hy or b.state.pos[0]>hx or b.state.pos[0]<lx: # artificial condition for particles deletion,
	    		bodiesToBeDeleted.append(b)
	for b in bodiesToBeDeleted:
  		O.bodies.erase(b.id)
	calculaPeso()

def calculaporosidade():
#	array=[]
	if os.path.exists('analise_porosidade_variando_young_por_altura_'+str(altura)+'.csv'):
		conf='a+'
	else:
		conf='w'
	with open('analise_porosidade_variando_young_por_altura_'+str(altura)+'.csv', conf,newline='') as csvfile:
		writer=csv.writer(csvfile)
		if conf=='w':
			writer.writerow(["Porosidade","Tempo","eixo x/y","altura"])
			writer.writerow([str(voxelPorosity(800,((hx/5)*2,(hy/5)*2,0.01),((hx/5)*4,(hy/5)*4,0.04))),str(O.time),str(hx),str(altura)])
		else:
			writer.writerow([str(voxelPorosity(800,((hx/5)*2,(hy/5)*2,0.01),((hx/5)*4,(hy/5)*4,0.04))),str(O.time),str(hx),str(altura)])

global amax
amax=0
def calc_altura():
	global amax
	for b in O.bodies:
		if isinstance(b.shape, Sphere):
			if b.state.pos[2]>amax:
				amax=b.state.pos[2]

def adicionaPeso(peso):
	calc_altura()
	peso1=[lx+0.0001,ly+0.0001,amax+0.003]
	peso2=[lx+0.0001,hy-0.0001,amax+0.003]
	peso3=[hx-0.0001,ly+0.0001,amax+0.003]
	peso4=[hx-0.0001,hy-0.0001,amax+0.003]
	facepeso1=yade.utils.facet([peso1,peso2,peso3],dynamic=True,fixed=False )
	facepeso2=yade.utils.facet([peso2,peso3,peso4],dynamic=True,fixed=False)
	facepeso1.state.mass = peso/2 # (!)
	facepeso1.state.inertia = (0,0,0) # (!)	
	facepeso2.state.mass = peso/2 # (!)
	facepeso2.state.inertia = (0,0,0) # (!)
	O.bodies.appendClumped([facepeso1,facepeso2])
	O.bodies[len(O.bodies)-1].state.blockedDOFs="xyXYZ"
	O.bodies[len(O.bodies)-2].state.blockedDOFs="xyXYZ"
	O.bodies[len(O.bodies)-3].state.blockedDOFs="xyXYZ"
	momento=O.time+0.01
	O.engines=O.engines+[PyRunner(command='O.pause()',virtPeriod=momento)]
	O.engines=O.engines+[PyRunner(virtPeriod=momento, command='calculaporosidade()')]


O.dt=0.9*utils.PWaveTimeStep()
O.saveTmp()
O.run()
O.wait()


for each physical and contact model, I changed the properties and materials of the particles and walls (when I thought it was necessary) and also, the contact and physical models, as described below.

elPerfPl Model:

wall -> FrictMat(density=2700,young=7e10,poisson=0.334,label="parede")
particles -> FrictMat(density=1163,young=2.6e6,poisson=0.25,frictionAngle=radians(22.6),label="particle")
InteractionLoop(
      [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom(),Ig2_Wall_Sphere_L3Geom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_L3Geom_FrictPhys_ElPerfPl()]
)

HertzMindlin Model:

wall -> ViscElMat(density=2700,young=7e10,poisson=0.334,label="parede")
particles -> ViscElMat(density=1163,young=2.6e6,poisson=0.25,frictionAngle=radians(22.6),label="particle")
InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_MindlinPhys()],
      [Law2_ScGeom_MindlinPhys_Mindlin()] #modelo linear de contato
)

ViscElPhys_Basic Model:

wall -> ViscElMat(density=2700,young=7e10,poisson=0.334,label="parede")
particles -> ViscElMat(density=1163,young=1.0e7,poisson=0.25,frictionAngle=radians(22.6),label="particle")
InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
      [Law2_ScGeom_ViscElPhys_Basic()] #modelo linear de contato
)

About using non-spherical particles like clumps, this option is in my
plans but, for preliminarly tests, spherical particles have less
computational costs to simulate, and, we hope that spherical particles
could be used for the study with a good similarity to the real particles
behavior.

Thanks for your help. :)

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.