yade-users team mailing list archive

[Grace Mejico <question685401@xxxxxxxxxxxxxxxxxxxxx>]

New question #685401 on Yade:
https://answers.launchpad.net/yade/+question/685401

Good evening, 
I'm from Peru and I'm part of an investigation that tries to represent the pirca performance at seismic action. The pirca is a dry stone retaining very used here.
For making the model, I was trying to use spherepack with clumps to represent the blocks of rock form.

First, If i use a volume like a  regular prism using vector,it works perfectly but what i want is a trapezoidal form (cross section) of the dry stone retaining wall. I tried many times but i couldnt find a way to make it work,thats why im sharing the code i was working on and hopefully i can have a hand with this problem. Is there a way to put the clumps in a volume with a trapezoidal form as cross section?

Second,I also want to recognize the clump that has the greatest displacement thats why i wrote that code . I would like to know if there is another way to do it.

listTopBodies=[]
for z1 in O.bodies:
	if isinstance(z1.shape,Sphere) and z1.state.pos[2]<=H+0.1 and z1.state.pos[2]>H-3.000*r and z1.state.pos[0]>(L/2.000)-2.000*r and z1.state.pos[0]<(L/2.000)+2.000*r and z1.state.pos[1]>B/2.000:
		listTopBodies.append(z1.id)
i=max(listTopBodies)

def addPlotData():
	plot.addData(Ypos=O.bodies[i].state.pos[1],Zpos=O.bodies[i].state.pos[2],iteracion=O.iter,unbal=unbalancedForce())
plot.plots={'Ypos':('Zpos'),'iteracion':('unbal')}

--------------------------------------------------------------------------------------------------------------------------------------------------------------------
 Here i share the complete code for the model:

# MODELO NUMERICO CON CLUMPS *************************
from yade import pack, qt, plot
from numpy import arange
import numpy as np
import random
# Datos de entrada
r=0.07 
H=1 #altura del muro
B=.45 #base mayor
b=.25 #base menor
L=2 #longitud del muro
#EN SERVICIO
EmpujeTop=2400 #Pa
EmpujeBottom=9100 #Pa
damp=0.7
aceloutofplane=0
ii=0
# Definicion del material mat1 para la pirca **********
mat1 = JCFpmMat()
mat1.cohesion = 0 #32.37e3 #Pa
mat1.density = 2233.2 #kg/m3
mat1.frictionAngle = radians(58.74) #rad
mat1.poisson = 0.18
mat1.young = 37.04e9 #Pa

# DEFINICION DE LA JUNTA  ******************************************************
mat1.jointCohesion = 232500 
mat1.jointFrictionAngle = radians(40)
mat1.jointNormalStiffness= 1.75e8
mat1.jointShearStiffness=0.875e8
#mat1.jointDilationAngle = 0
#mat1.jointTensileStrength = 0
O.materials.append(mat1)

# GEOMETRIA ********************************************************************
a=(B-b)*.5
Vol1=pack.inParallelepiped(o=(0,0,0), a=(L,0,0), b=(0,B,0), c=(0,0,H))
Vol2=pack.inParallelepiped(o=(L,B+a,0), a=(0,B+a,0), b=(0,a+b,H), c=(L,a+b,H))
Vol3=pack.inParallelepiped(o=(L,-a,0), a=(0,-a,0), b=(0,a,H), c=(L,a,H))
Vol4=Vol1-Vol2-Vol3-Vol3

# Conjunto de esferas representadas como centros y radios
sp=pack.SpherePack()
biesfera=pack.SpherePack([((0,0,0),.033),((.033,0,0),.033)]) # clump de dos esferas
triesfera=pack.SpherePack([((0,0,0),.04),((.02,.02,0),.04),((0,.02,.02),.04)]) # clump de tres esferas
tetraesfera=pack.SpherePack([((0,0,0),.033),((.033,0,0),.033),((0,.033,0),.033),((0,0,.033),.033)]) # clump decuatro esferas
sp.makeClumpCloud(Vol4,[biesfera,triesfera, tetraesfera]) # Empaquetamiento suelto y aleatoriode clumps
sp.toSimulation() # llama las esferas a la simulacion
#Imprime la cantidad total de cuerpos
Q=len(O.bodies)
print 'numero de cuerpos', Q
# Definicion del piso
piso=utils.wall(position=Vector3(0,0,0), axis=2, sense=0, material=mat1, color=(0.5,0.5,0.5))
O.bodies.append(piso)

# Definicion de las fuerzas en aplicacion
listSphereA=[]
listSphereB=[]
listSphereC=[]
listSphereD=[]
for x1 in O.bodies:
	if isinstance(x1.shape,Sphere) and x1.state.pos[1]<3*r and x1.state.pos[2]<=H and x1.state.pos[2]>3*H/4.000:
		listSphereA.append(x1.id)
for x2 in O.bodies:
	if isinstance(x2.shape,Sphere) and x2.state.pos[1]<3*r and x2.state.pos[2]<=3*H/4.000+0.1 and x2.state.pos[2]>H/2.000+0.1:
		listSphereB.append(x2.id)
for x3 in O.bodies:
	if isinstance(x3.shape,Sphere) and x3.state.pos[1]<3*r and x3.state.pos[2]<=H/2.000+0.1 and x3.state.pos[2]>H/4.000+0.1:
		listSphereC.append(x3.id)
for x4 in O.bodies:
	if isinstance(x4.shape,Sphere) and x4.state.pos[1]<3*r and x4.state.pos[2]<=H/4.000+0.1 and x4.state.pos[2]>0.000+0.1:
		listSphereD.append(x4.id)
ha=1*H/8.000
hb=3*H/8.000
hc=5*H/8.000
hd=7*H/8.000
m=EmpujeBottom-EmpujeTop
Ea=m*ha+EmpujeTop
Eb=m*hb+EmpujeTop
Ec=m*hc+EmpujeTop
Ed=m*hd+EmpujeTop

Fa=Ea*(H/4.000)*L/len(listSphereA)
Fb=Eb*(H/4.000)*L/len(listSphereB)
Fc=Ec*(H/4.000)*L/len(listSphereC)
Fd=Ed*(H/4.000)*L/len(listSphereD)

print 'Fzas en c/ esfera (N)', Fa, Fb, Fc, Fd

def aplicarFuerzaA():
	for i in listSphereA:
		O.forces.setPermF(i,(0,Fa,0))
def aplicarFuerzaB():
	for i in listSphereB:
		O.forces.setPermF(i,(0,Fb,0))
def aplicarFuerzaC():
	for i in listSphereC:
		O.forces.setPermF(i,(0,Fc,0))
def aplicarFuerzaD():
	for i in listSphereD:
		O.forces.setPermF(i,(0,Fd,0))

# CONDICIONES DE BORDE *********************************************************
Base=utils.wall(0.1,axis=2,sense=0)
Borde1=utils.wall(0.22, axis=0, sense=0)
Borde2=utils.wall(L-.22, axis=0, sense=0)
Borde3=utils.wall(-.0, axis=1, sense=0)
O.bodies.append(Base)
O.bodies.append(Borde1)
O.bodies.append(Borde2)
O.bodies.append(Borde3)

listSphereBase=[]

for y1 in O.bodies:
	if isinstance(y1.shape,Sphere) and y1.state.pos[2]<=3*r+0.1:
		listSphereBase.append(y1.id)
#for y2 in O.bodies:
#	if isinstance(y2.shape,Sphere) and y2.state.pos[0]<=2.5*r:
#		listSphereBorde1.append(y2.id)
#for y3 in O.bodies:
#	if isinstance(y3.shape,Sphere) and y3.state.pos[0]>L-2.5*r:
#		listSphereBorde2.append(y3.id)
for i in listSphereBase:
		O.bodies[i].state.blockedDOFs='xyzXYZ'
#for i in listSphereBorde1:
#		O.bodies[i].state.blockedDOFs='xyzXYZ'
#for i in listSphereBorde2:
#		O.bodies[i].state.blockedDOFs='xyzXYZ'

def checkUnbalanced():
	print ('iter %d, unbalanced forces = %f'%(O.iter, utils.unbalancedForce()))
#	if unbalancedForce()<.00001:

listTopBodies=[]
for z1 in O.bodies:
	if isinstance(z1.shape,Sphere) and z1.state.pos[2]<=H+0.1 and z1.state.pos[2]>H-3.000*r and z1.state.pos[0]>(L/2.000)-2.000*r and z1.state.pos[0]<(L/2.000)+2.000*r and z1.state.pos[1]>B/2.000:
		listTopBodies.append(z1.id)
i=max(listTopBodies)

def addPlotData():
	plot.addData(Ypos=O.bodies[i].state.pos[1],Zpos=O.bodies[i].state.pos[2],iteracion=O.iter,unbal=unbalancedForce())
plot.plots={'Ypos':('Zpos'),'iteracion':('unbal')}

# MOTOR DEL ALGORITMO DEM *******************************
O.engines=[
		ForceResetter(),
		InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
		InteractionLoop(
			[Ig2_Sphere_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
			[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
			[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM()]
		),
		NewtonIntegrator(gravity=(0,aceloutofplane,-9.81),damping=damp),
		PyRunner(firstIterRun=10000,command='aplicarFuerzaA()'),
		PyRunner(firstIterRun=10000,command='aplicarFuerzaB()'),
		PyRunner(firstIterRun=10000,command='aplicarFuerzaC()'),
		PyRunner(firstIterRun=10000,command='aplicarFuerzaD()'),
		PyRunner(iterPeriod=10000,command='checkUnbalanced()'),
		PyRunner(iterPeriod=1,command='addPlotData()')
]

# DETALLES FINALES *************************************************************
O.dt=0.5*PWaveTimeStep() # establece el delta de tiempo como una fraccion del tiempo critico
plot.plot(subPlots=True) # llama a los sub plots
plot.live=True # ploteo en tiempo real
qt.Controller() # abre la ventana del controlador
V=qt.View() # abre la ventada de la vista
R=yade.qt.Renderer() # llama a la renderizacion
R.bgColor=(1.,1.,1.) # definel color blan
V.screenSize = (900,500) # tamano de pantalla


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