New question #270336 on Yade:
https://answers.launchpad.net/yade/+question/270336
Hello everyone,
I want to calculate the potential energy of the beam and the work performed on
the beam in 3-points bending test.
To calculate the energy potential in each particle i use this formula:
'Ub(particle) = m(mass of particle)*g*d(displacement of particle)
and the work performed in each particle on the boundary of the model, i use
this formula:
'W(particle)=F(force exerced on the particle)*d(displacement of particle)
My question is how can i determine the force exerced on each particle and the
displacement of each particle?
Thanks so much
Jabrane
from yade import ymport, utils, pack, export
from yade import plot
from pylab import *
import math
#---------------- SIMULATIONS DEFINED HERE (assembly, material, boundary
conditions)
#### packing (previously constructed)
PACKING='X240Y80Z30_400'
OUT=PACKING+'_flexion_3_points'
#### Simulation Control
DAMP=0.4 # numerical damping
saveData=1000 # data record interval
iterMax=10000 # maximum number of iteration (to be adjusted)
saveVTK=1000 # Vtk files record interval
#### Material microproperties -> Lac du Bonnet granite (cf. A DEM model for
soft and hard rock, Scholtes & Donze, JMPS 2013)
intR=1.6684# allows near neighbour interaction and coordination number K=13
(determined with coordinationNumber.py -> to be adjusted for every packing)
DENS=4000 # could be adapted to match material density:
dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> poro?
YOUNG=65e9
FRICT=10
ALPHA=0.4
TENS=8e6
COH=160e6
#### material definition
sphereMat = JCFpmMat(type=1,density=DENS,young=YOUNG,poisson =
ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
wallMat = JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))
for mat in (sphereMat,wallMat):
O.materials.append(mat) # then wallMat will be used if material is not
specified
#### preprocessing to get dimensions
O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0)))
dim=utils.aabbExtrema()
xinf=dim[0][0]
xsup=dim[1][0]
X=xsup-xinf
yinf=dim[0][1]
ysup=dim[1][1]
Y=ysup-yinf
zinf=dim[0][2]
zsup=dim[1][2]
Z=zsup-zinf
r=X/15.
#### preprocessing to get dimensions
R=0
Rmax=0
numSpheres=0.
for o in O.bodies:
if isinstance(o.shape,Sphere):
numSpheres+=1
R+=o.shape.radius
if o.shape.radius>Rmax:
Rmax=o.shape.radius
Rmean=R/numSpheres
O.reset() # all previous lines were for getting dimensions of the packing to
create walls at the right positions (below) because walls have to be genrated
after spheres for FlowEngine
O.bodies.append(geom.facetCylinder(center=(xinf+X/12.,yinf-0.9*r,Z/2.),radius=r,height=Z,orientation=Quaternion((1,
0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # bas gauche
O.bodies.append(geom.facetCylinder(center=(xinf+11*X/12.,yinf-0.9*r,Z/2.),radius=r,height=Z,orientation=Quaternion((1,
0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # bas droite
piston=O.bodies.append(geom.facetCylinder(center=(0.5*X,Y+0.9*r,Z/2.),radius=r,height=Z,dynamic=False,orientation=Quaternion((1,
0, 0), 0),segmentsNumber=20,wire=False,material=wallMat)) # haut
p0=O.bodies[piston[0]].state.pos[1]
### packing
beam=O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
# Definition of the facets for joint's geometry
alpha = 90 #angle of the crack with horizontal (pi/4.)
c = 0.375*Y # length of the crack
# 4 points pour l'entaille
ptA = gts.Vertex( X/2., c, 8./7.*Z)
ptB = gts.Vertex( X/2., 0, 8./7.*Z)
ptApr = gts.Vertex(X/2., c, -1./7.*Z)
ptBpr = gts.Vertex(X/2., 0, -1./7.*Z)
e1 = gts.Edge(ptA,ptB)
e2 = gts.Edge(ptA,ptApr)
e3 = gts.Edge(ptApr,ptB)
f1 = gts.Face(e1,e2,e3)
e4 = gts.Edge(ptB,ptBpr)
e5 = gts.Edge(ptBpr,ptApr)
f2 = gts.Face(e4,e5,e3)
s1 = gts.Surface()
s1.add(f1)
s1.add(f2)
facet = gtsSurface2Facets(s1,wire = False,material=wallMat)
O.bodies.append(facet)
## Determinate the particle in the boundary
limite=[]
limite=range(0, len(beam))
for i in range(0, len(beam)):
if X/3. <= O.bodies[beam[i]].state.pos[0] <= 2*X/3. or Y/3. <=
O.bodies[beam[i]].state.pos[1] <= 2*Y/3. or Z/3. <=
O.bodies[beam[i]].state.pos[2] <= 2*Z/3. :
limite[i] = 0
else:
limite[i] = beam[i]
for i in range(0, len(limite)):
try:
limite.remove(0)
except ValueError:
pass
## Determinate the initial position of the particle on the boundary
poslimx=[]
poslimy=[]
poslimz=[]
poslimx=range(0,len(limite))
poslimy=range(0,len(limite))
poslimz=range(0,len(limite))
for i in range(0,len(limite)):
poslimx[i]=O.bodies[limite[i]].state.pos[0]
for i in range(0,len(limite)):
poslimy[i]=O.bodies[limite[i]].state.pos[1]
for i in range(0,len(limite)):
poslimz[i]=O.bodies[limite[i]].state.pos[2]
### set a color to the spheres
for o in O.bodies:
if isinstance(o.shape,Sphere):
o.shape.color=(0.7,0.5,0.3)
#---------------- ENGINES DEFINED HERE
#### simulation is defined here (DEM loop, interaction law, servo control,
recording, etc...)
##### simulation piston's motion
for i in range(0,len(piston)):
O.bodies[piston[i]].state.vel[1]=-1
##### simulation beam's grains motion
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb'),Bo1_Facet_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Facet_Sphere_ScGeom()],
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4,
defaultDt=0.1*utils.PWaveTimeStep()),
NewtonIntegrator(damping=DAMP,label="newton"),
PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','bstresses','jcfpm','cracks'],Key=OUT,label='vtk')
]
## Determinate the initial position of all particles of the model
posx=[]
posy=[]
posz=[]
posx=range(0,len(beam))
posy=range(0,len(beam))
posz=range(0,len(beam))
for i in range(0,len(beam)):
posx[i]=O.bodies[beam[i]].state.pos[0]
for i in range(0,len(beam)):
posy[i]=O.bodies[beam[i]].state.pos[1]
for i in range(0,len(beam)):
posz[i]=O.bodies[beam[i]].state.pos[2]
## Calculate the work on the boundary of model
W=[]
W=range(0,len(limite))
## W est le travail effectue sur les billes du limite de modele
#for i in range(0, len(W)):
#
W[i]=utils.sumFacetNormalForces(ids=limite[i],axis=0)*abs(O.bodies[limite[i]].state.pos[0]-poslimx[i])+utils.sumFacetNormalForces(ids=limite[i],axis=1)*abs(O.bodies[limite[i]].state.pos[1]-poslimy[i])+utils.sumFacetNormalForces(ids=limite[i],axis=2)*abs(O.bodies[limite[i]].state.pos[2]-poslimz[i])
## gravitational acceleration
g=[0,0,-9.8]
## calculate the potential energy of model
Ub=[]
Ub=range(0,len(beam))
## Ub is the potential energy of the model
for i in range(0, len(Ub)):
Ub[i]=O.bodies[beam[i]].state.mass*abs(g[2])*abs(O.bodies[beam[i]].state.pos[2]-posz[i])
# if you want to plot during simulation
plot.plots={'i':('f')}
#plot.plot()
#---------------- SIMULATION STARTS HERE
#### manage interaction detection factor during the first timestep and then set
default interaction range ((cf. A DEM model for soft and hard rock, Scholtes &
Donze, JMPS 2013))
O.step();
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.
O.run(iterMax)
--
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
