# How to make some stress plates consist of 3 circles (that have radius gradually (0.5,1,2)?

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
import math
from pylab import rand #for sand color

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
#[b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]
O.bodies.state.blockedDOFs='xyzXYZ'
O.bodies.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 #the last particles is the boulder
O.materials.append(FrictMat(young=1e8,poisson=0.2,density=2700,frictionAngle=25))
O.bodies.state.mass=2000

if abs(O.forces.f(boulder.id)):
O.run()

b=O.bodies
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>=-2) and (m.state.pos<=2):
if (m.state.pos>=-1) and (m.state.pos<=1) and (m.state.pos<-.5):
ForceP=abs(O.forces.f(m.id))
StP=ForceP/area1
ForceP1=abs(O.forces.f(m.id))
StP1=ForceP1/area2
ForceP2=abs(O.forces.f(m.id))
StP2=ForceP2/area3

if O.iter>1000:

globals().update(locals())

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

plot.plot(subPlots=True)

O.saveTmp()

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2021-03-15: #1

Hello,

I am sorry, but I could not understand at all, what you want to achieve..
Maybe you could provide a link to a picture or scheme?

> 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.

what is "circle of stress plate"?
what is "contact force"?

> 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.

what is "the facet"?

A few notes to the code:

I got "IndexError: Body id out of range."

Try to avoid this structure, put the PyRunner directly to O.engines=[...] definition

> PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),

using iterPeriod or virtPeriod is preferable over realPeriod

> global boulder
> boulder=O.bodies
> O.bodies.state.mass=2000

"global boulder" is ugly, not necessary and misleading here.
You can make a global 'boulder' variable
boulder = O.bodies
after sp.toSimulation(), without "ugly" global keyword.
Or just use O.bodies, is not much longer and maybe even more clear what it is and you are already using it this way (inside checkUnbalanced) anyway..

> if abs(O.forces.f(boulder.id)):
> O.run()

seems like it is not needed, O.run should not be there (the simulation is running already), so you can set checker.command='stress_rad1()' directly in checkUnbalanced

cheers
Jan

 Revision history for this message Defri (daredefri) said on 2021-03-15: #2

1. The following link are the images about what I want to achieve https://cutt.ly/GzVIYuA (I am trying to get stress plate due to calculate the stress increment from boulder)

2. what is "circle of stress plate"?
what is "contact force"?

I attach the link that is the images to describe what I need to achieve

3. what is "the facet"?
I attach the link that is the images to describe what I need to achieve

4. I got "IndexError: Body id out of range."

Yes, it is because the boulder which is O.bodies has not appeared yet, when the boulder has appeared, then IndexError becomes stopped.
I have tried to change it, and I hope you could help me to check, do you think it is normal now to run it?

5. -PyRunner
-realperiod
-global boulder
-O.run

I have solved all by following your directions

THANK YOU

this is my existing code now:

##Sphere Cylinder pack
import math
from pylab import rand #for sand color

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
b=O.bodies
b.state.blockedDOFs='xyzXYZ'
b.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()',iterPeriod=100,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
O.materials.append(FrictMat(young=1e8,poisson=0.2,density=2700,frictionAngle=25))
O.bodies.state.mass=2000

abs(O.forces.f(boulder.id))

b1=O.bodies
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>=-2) and (m.state.pos<=2):
if (m.state.pos>=-1) and (m.state.pos<=1) and (m.state.pos<-.5):
ForceP=abs(O.forces.f(m.id))
StP=ForceP/area1
ForceP1=abs(O.forces.f(m.id))
StP1=ForceP1/area2
ForceP2=abs(O.forces.f(m.id))
StP2=ForceP2/area3

if O.iter>1000:

globals().update(locals())

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

plot.plot(subPlots=True)

O.saveTmp()

I hope you could understand what I really want to achieve
Cheers!

 Revision history for this message Jan Stránský (honzik) said on 2021-03-15: #3

Thanks for more information. I am too tired now to dive into the problem, just briefly:
You want to measure time evolution of stress at certain points at the bottom plane?
Jan

 Revision history for this message Defri (daredefri) said on 2021-03-15: #4

thank you for your response Mr. Stránský
I will wait whenever you're ready

The important thing is the stress (that derive from Force and Area of the certain point) and the time evolution as well (like we can put in the graph for x axis is time, y is stress) so the graph will show each of stress.

Cheers!

 Revision history for this message  Jan Stránský (honzik) said on 2021-03-16: #5

Ok, now it is clear what you want - measure stress at certain point and certain time :-)
Then, if you measure it in circles or whatever shape, how often etc. is easy extension.

First of all, in the Discrete element method, everything is discrete, i.e. you have just forces, no stress.
You can somehow convert forces to stress. Here, you should know what you want and how you want to calculate it (!), what approach is suitable, which not....

Some brainstorming of possible solutions:
- select "bottom spheres" and estimate the plate stress from the stresses of these spheres
- discretize the bottom plate with more fine mesh and then you can measure the stress from the forces acting on the smaller "elements"
- leave one big plate, but estimate the stress from "nearby" interactions.

A very scientific approach would be taking all of them and compare them :-)

The last option seems to me as the for now easiest approach and also a nice approach in general.
something like:
###
##Make a floor
bottom = wall((0,0,-1.5),axis=2)
O.bodies.append(bottom)

def getBottomStress(point):
interactions = bottom.intrs() # get all bottom interactions
interactions = [i for i in interactions if (i.geom.contactPoint-point).norm() < limitDistance] # filtering out far away points
cps = [i.geom.contactPoint for i in interactions] # contact points
fs = [i.phys.normalForce + i.phys.shearForce for i in interactions] # forces
stress = somehowConvertForcesAndContactPointsToStress(fs,cps) # here you need to know what you want
return stress

# some running
stress34 = getBottomStress((3,4,0))
###

cheers
Jan

 Revision history for this message Defri (daredefri) said on 2021-03-18: #6

Thank you for your help Mr. Stránský

I am using suggestion for your easiest and general approach (last part).
Apparently I am still having many failure note in terminal.
It might be from my side, don't know how to implement my code after I use your suggestion syntax.

Thank you anyway.

 Revision history for this message Jan Stránský (honzik) said on 2021-03-18: #7

> Apparently I am still having many failure note in terminal.

please provide actual code and the errors.
I wrote the code "by heart" (hence the "something like" note), there might be some problems and the need of modifications. It was meant to show the core idea..

cheers
Jan

 Revision history for this message Defri (daredefri) said on 2021-03-18: #8

The actual code:

##Sphere Cylinder pack
import math
from pylab import rand #for sand color

print('the ids of the facets:',facets)

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
b=O.bodies
b.state.blockedDOFs='xyzXYZ'
b.state.vel=(0,0,0)
O.materials.append(FrictMat(young=1e7,poisson=.3,density=2700,frictionAngle=18))

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

##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()',iterPeriod=100,label='checker'),
PyRunner(command='getbottomstress()')
#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
O.materials.append(FrictMat(young=1e8,poisson=0.2,density=2700,frictionAngle=25))
O.bodies.state.mass=2000
checker.command='getbottomstress()'

def getbottomstress():
b1=O.bodies
interactions = bottom.intrs()
interactions = [i for i in interactions if (i.geom.area1).norm()]
cp = [i.geom.area1 for i in interactions]
fs = [i.phys.normalForce + i.phys.shearForce for i in interactions]
StP = (fs,cp)/area1

if O.iter>1000:

globals().update(locals())

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

plot.plot(subPlots=True)

O.saveTmp()
##############################

What I am trying to do based on my code, is that I am trying to build a circle stress plate, then I expect to get the stress by dividing the Force and the area of circle stress plate.
I am trying to find the stress by assigning the boulder (O.bodies) to the spheres, and expecting to get the stress, and generate the value by making graph.

Turns out I got the failure note:
AttributeError Traceback (most recent call last)

95 interactions = bottom.intrs()
---> 96 interactions = [i for i in interactions if (i.geom.area1).norm()]
97 cp = [i.geom.area1 for i in interactions]
98 fs = [i.phys.normalForce + i.phys.shearForce for i in interactions]

95 interactions = bottom.intrs()
---> 96 interactions = [i for i in interactions if (i.geom.area1).norm()]
97 cp = [i.geom.area1 for i in interactions]
98 fs = [i.phys.normalForce + i.phys.shearForce for i in interactions]

AttributeError: 'ScGeom' object has no attribute 'area1'

Thank you
Cheers!

 Revision history for this message Jan Stránský (honzik) said on 2021-03-18: #9

The key point here is that you should know what and how you want to compute.

Randomly modifying the code (sorry, but the last code seems like that) is rarely the way to go.

So, let's take the "last option".
You know all the interactions (i.e. essentially all the forces and their positions) acting on the bottom plate.
How to convert this discrete informations to continuum-like stress?
In general, this is not unique transformation and you can have many possible ways.
One crude option is to define a region, select forces acting in the region, sum them and divide by the area of the region.
But it will be very sensitive to the cases "at the boundary" - a force just inside or just outside, a very small difference in action point could result in a jump change in the resulting stress.

So a better option might (or might not :-) be to do some "smoothing" - weighted average with weights being distances from the region boundary.

We can help you how to code it, maybe also give some suggestions, but you have do the decision what/how to compute.

> interactions = [i for i in interactions if (i.geom.area1).norm()]

i.geom is ScGeom . It does not have 'area' nor 'area1'
What should the condition "if (i.geom.area1).norm()" do?

> cp = [i.geom.area1 for i in interactions]
> fs = [i.phys.normalForce + i.phys.shearForce for i in interactions]
> StP = (fs,cp)/area1

cp and fs are lists of vectors. cp is list of contact points, fs is list of forces.
What should "(fs,cp)" do?
Now, "(fs,cp)" would be tuple of lists of vectors, something like ([Vector3(...), Vector3(...), ...], [Vector3(...), Vector3(...), ...])
Dividing this by area1 would result in error

> circle stress plate

Do you want an average stress inside the entire circle?
Or average stress at the boundary of the circle?
Or ... ?

> I am willing to add more circle stress plate

once you have one circle, more circles is easy task

cheers
Jan

 Revision history for this message Defri (daredefri) said on 2021-03-18: #10

> What should the condition "if (i.geom.area1).norm()" do?
I thought to derive the contact point we introduce the area with that syntax, in order to get the stress

> What should "(fs,cp)" do?
Now, "(fs,cp)" would be tuple of lists of vectors, something like ([Vector3(...), Vector3(...), ...], [Vector3(...), Vector3(...), ...])
Dividing this by area1 would result in error

I thought it will have such a "coordinate" like x y, represented by fs and cp value then I will get the Force value.

>Do you want an average stress inside the entire circle?
Or average stress at the boundary of the circle?
Or ... ?

the entire circle, and then after adding another circle, I want to derive also average stress inside the next circle (which has a slight bigger radius)

Thank you

 Revision history for this message Jan Stránský (honzik) said on 2021-03-18: #11

> AttributeError: 'ScGeom' object has no attribute 'area1'

the error message should be clear.
Always try to get rid of errors in the very first place. The rest of the function (getBottomStress) has no meaning if you have error at the very beginning of it.
As discussed, in DEM (D = discrete), everything is discrete, so interaction.geom.area (area being not so much discrete) is likely not to be available (which actually is the case).

> if (i.geom.area1).norm()

even i.geom.area1 exists (if assuming area, why area1 and not area?), what would be the meaning of "area.norm()"? And what "if area.norm()", which would always evaluate to True?

> with that syntax

apart from syntax, code has some context and semantics. First, try to understand the code before modifying it.

> I thought ...

you can always put some prints to see what the variables are (in case of no errors, of course).

Have a try. But as discussed, this is a "crude" approach, I **personally** would add some "smoothing".
###
# center ... center of the stress circle
interactions = bottom.intrs()
interactions = [i for i in interactions if (i.geom.contactPoint-center).norm() < radius]
fns = [i.phys.normalForce for i in interactions]
fnTot = sum(fn.norm() for fn in fns)
stress = fnTot / area
return stress
# TODO shear stress?
###

cheers
Jan

 Revision history for this message Bruno Chareyre (bruno-chareyre) said on 2021-03-18: #12

On this basis you could generate stress for each (x,y) on a grid (possibly with square domains instead of circles, with square size = grid size).

Then the average on any subdomain (e.g. a circle) would be the average on the contained grid points. That's just one of many ways.

I would say you need to separate problems to help people helping you:
1/ suppose you have a list of forces and (x,y,z) coordinates of the corresponding contact points. What do you want to do of this (mathematics)?
2/ how to obtain such a list of forces and (x,y,z) coordinates of the corresponding contact points with yade (and syntax problems)

Bruno

 Revision history for this message Defri (daredefri) said on 2021-04-05: #13

Thanks Jan Stránský, that solved my question.