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

Description changed to:
Dear All

I am trying to simulate 3 circles of stress plates that have radius
dimension gradually, which has relation with contact force and the area
of each circle.

In order to apply the force that corresponds to each area of circles, I
need to have 3 stress plates with radius 0.5,1, and 2 in XY plane in the
center of the facet.

How can I add those stress plates and how to make some relations among
of each contact force and the area of circle (stress plate)?

Could you please tell me how to do that?

this is my existing code:

##Sphere Cylinder pack
from yade import pack,plot,utils,export
import math
from pylab import rand #for sand color

O.bodies.append(geom.facetBox((0,0,0),(5,2,1.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((-5,-2,-1.5),(5,2,1.5),rMean=.14,rRelFuzz=.1,num=3000)
sp.toSimulation(color=(0.6+0.15*rand(),0.5+0.15*rand(),0.15+0.15*rand()))

##Define material of the grains
#O.bodies[0]
#[b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]
O.bodies[0].state.blockedDOFs='xyzXYZ'
O.bodies[0].state.vel=(0,0,0)
O.materials.append(FrictMat(young=1e7,poisson=.3,density=2700,frictionAngle=18))

##Make a floor
O.bodies.append(wall((0,0,-1.5),axis=2))

##Engines and Constitutive Law
O.engines=[ForceResetter(),
        
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]), 
#move for sphere, facet, wall
        
InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
 #interaction between them
   [Ip2_FrictMat_FrictMat_FrictPhys()], # ip2 list (just one list!)
   [Law2_ScGeom_FrictPhys_CundallStrack()]), # law2 list
   NewtonIntegrator(damping=.3,gravity=[0,0,-9.81]),
   PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),
   #VTKRecorder(fileName='tes',recorders=['all'],iterPeriod=100),
   ]

O.dt=.9*PWaveTimeStep()

def checkUnbalanced():
        if O.iter<1000:return # at the very start, unbalanced force can be low 
as there is only few contacts, but it does not mean the packing                 
                                                         is stable
        global boulder  #without this line, the BOULDER variable would only 
exist inside this function
        boulder=O.bodies[0]      #the last particles is the boulder
        
O.materials.append(FrictMat(young=1e8,poisson=0.2,density=2700,frictionAngle=25))
        O.bodies.append(utils.sphere(center=(0,0,5),radius=2,color=[0,1,1]))
        O.bodies[0].state.mass=2000
        checker.command='unloadBoulder()'

def unloadBoulder():
        if abs(O.forces.f(boulder.id)[0]):
                O.run() 
        checker.command='stress_rad1()'
                
O.engines=O.engines+[PyRunner(command='stress_rad1()',iterPeriod=50)]

def stress_rad1():
        rad1=O.bodies[3011].shape.radius
        b=O.bodies[3011]
        area1=pi*rad1**2
        area2=pi*((rad1+.8)**2-rad1**2)
        area3=pi*((rad1+1.2)**2-(rad1+.8)**2)
        ForceP=0
        StP=0
        ForceP1=0
        StP1=0
        ForceP2=0
        StP2=0
        for i in range (10,3009):
                m=O.bodies[i]
                if(m.state.pos[1]>=-2) and (m.state.pos[1]<=2):
                        if (m.state.pos[0]>=-1) and (m.state.pos[0]<=1) and 
(m.state.pos[2]<-.5):
                                if 
(m.state.pos[0]**2+m.state.pos[1]**2)**.5<=rad1:     #1st area
                                                ForceP=abs(O.forces.f(m.id)[2])
                                                StP=ForceP/area1
                                if 
((m.state.pos[0]**2+m.state.pos[1]**2)**.5<=rad1+.8) and 
((m.state.pos[0]**2+m.state.pos[1]**2)**.5>rad1):  #2nd area
                                                ForceP1=abs(O.forces.f(m.id)[2])
                                                StP1=ForceP1/area2
                                if 
((m.state.pos[0]**2+m.state.pos[1]**2)**.5<=rad1+1.2) and 
((m.state.pos[0]**2+m.state.pos[1]**2)**.5>rad1+.8):  #3rd area
                                                ForceP2=abs(O.forces.f(m.id)[2])
                                                StP2=ForceP2/area3
                                                
        if O.iter>1000:
                
plot.addData(z=b.state.pos[2],Stress=StP,Stress1=StP1,Stress2=StP2,i=O.iter,t=O.time)
                
globals().update(locals())

plot.plots={'i':('z'),'t':('Stress','Stress1','Stress2')}

plot.plot(subPlots=True)

O.saveTmp()


THANK YOU IN ADVANCE.

-- 
You received this question notification because your team yade-users 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

Reply via email to