# 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

from yade import pack,plot,

import math

from pylab import rand #for sand color

O.bodies.

sp=pack.

sp.makeCloud(

sp.toSimulation

##Define material of the grains

#O.bodies[0]

#[b.shape.radius for b in O.bodies if isinstance(

O.bodies[

O.bodies[

O.materials.

##Make a floor

O.bodies.

##Engines and Constitutive Law

O.engines=

InsertionSort

InteractionLo

[Ip2_

[Law2_

NewtonIntegr

PyRunner(

#VTKRecorder

]

O.dt=.9*

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

O.bodies.

O.bodies[

checker.

def unloadBoulder():

if abs(O.forces.

O.run()

checker.

O.engines=

def stress_rad1():

rad1=O.

b=O.bodies[3011]

area1=pi*rad1**2

area2=

area3=

ForceP=0

StP=0

ForceP1=0

StP1=0

ForceP2=0

StP2=0

for i in range (10,3009):

m=O.bodies[i]

if(m.

if (m.state.

if (m.state.

ForceP=

StP=

if ((m.state.

ForceP1=

StP1=

if ((m.state.

ForceP2=

StP2=

if O.iter>1000:

plot.

globals(

plot.plots=

plot.plot(

O.saveTmp()

THANK YOU IN ADVANCE.

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- Yade Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Jan Stránský

- Solved:
- 2021-04-05

- Last query:
- 2021-04-05

- Last reply:
- 2021-03-16

Jan Stránský (honzik) said : | #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 does "having radius dimension gradually" mean?

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:

> def stress_rad1():

> rad1=O.

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

> O.engines=

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

> PyRunner(

using iterPeriod or virtPeriod is preferable over realPeriod

> global boulder

> boulder=O.bodies[0]

> O.bodies[

"global boulder" is ugly, not necessary and misleading here.

You can make a global 'boulder' variable

boulder = O.bodies[0]

after sp.toSimulation(), without "ugly" global keyword.

Or just use O.bodies[0], is not much longer and maybe even more clear what it is and you are already using it this way (inside checkUnbalanced) anyway..

> def unloadBoulder():

> if abs(O.forces.

> O.run()

> checker.

seems like it is not needed, O.run should not be there (the simulation is running already), so you can set checker.

cheers

Jan

Defri (daredefri) said : | #2 |

1. The following link are the images about what I want to achieve https:/

2. what is "circle of stress plate"?

what does "having radius dimension gradually" mean?

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[3011] 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

-set checker.

I have solved all by following your directions

THANK YOU

this is my existing code now:

##Sphere Cylinder pack

from yade import pack,plot,

import math

from pylab import rand #for sand color

O.bodies.

sp=pack.

sp.makeCloud(

sp.toSimulation

##Define material of the grains

b=O.bodies[0]

b.state.

b.state.vel=(0,0,0)

O.materials.

##Make a floor

O.bodies.

##Engines and Constitutive Law

O.engines=

InsertionSort

InteractionLo

[Ip2_

[Law2_

NewtonIntegr

PyRunner(

PyRunner(

#VTKRecorder

]

O.dt=.9*

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.

O.bodies.

O.bodies[

checker.

checker.

def unloadBoulder():

abs(O.

def stress_rad1():

rad1=O.

b1=O.bodies[3011]

area1=pi*rad1**2

area2=

area3=

ForceP=0

StP=0

ForceP1=0

StP1=0

ForceP2=0

StP2=0

for i in range (10,3009):

m=O.bodies[i]

if(m.

if (m.state.

if (m.state.

ForceP=

StP=

if ((m.state.

ForceP1=

StP1=

if ((m.state.

ForceP2=

StP2=

if O.iter>1000:

plot.

globals(

plot.plots=

plot.plot(

O.saveTmp()

I hope you could understand what I really want to achieve

Thank you in advance

Cheers!

Jan Stránský (honzik) said : | #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

Defri (daredefri) said : | #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!

Jan Stránský (honzik) said : | #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,

O.bodies.

def getBottomStress

interactions = bottom.intrs() # get all bottom interactions

interactions = [i for i in interactions if (i.geom.

cps = [i.geom.

fs = [i.phys.normalForce + i.phys.shearForce for i in interactions] # forces

stress = somehowConvertF

return stress

# some running

stress34 = getBottomStress

###

cheers

Jan

Defri (daredefri) said : | #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.

Jan Stránský (honzik) said : | #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

Defri (daredefri) said : | #8 |

The actual code:

##Sphere Cylinder pack

from yade import pack,plot,

import math

from pylab import rand #for sand color

facets=

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

sp=pack.

sp.makeCloud(

sp.toSimulation

##Define material of the grains

b=O.bodies[0]

b.state.

b.state.vel=(0,0,0)

O.materials.

##Make a floor

bottom=

O.bodies.

##Engines and Constitutive Law

O.engines=

InsertionSort

InteractionLo

[Ip2_

[Law2_

NewtonIntegr

PyRunner(

PyRunner(

#VTKRecorder

]

O.dt=.9*

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.

O.bodies.

O.bodies[

checker.

def getbottomstress():

b1=O.bodies[2920]

rad1=b1.

area1=pi*rad1**2

interactions = bottom.intrs()

interactions = [i for i in interactions if (i.geom.

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:

plot.

globals(

plot.plots=

plot.plot(

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[2920]) to the spheres, and expecting to get the stress, and generate the value by making graph.

and if this is working, I am willing to add more circle stress plate that is slightly bigger, like this area2=pi*

Turns out I got the failure note:

AttributeError Traceback (most recent call last)

/usr/bin/yade in <module>

/usr/bin/yade in getbottomstress()

94 area1=pi*rad1**2

95 interactions = bottom.intrs()

---> 96 interactions = [i for i in interactions if (i.geom.

97 cp = [i.geom.area1 for i in interactions]

98 fs = [i.phys.normalForce + i.phys.shearForce for i in interactions]

/usr/bin/yade in <listcomp>(.0)

94 area1=pi*rad1**2

95 interactions = bottom.intrs()

---> 96 interactions = [i for i in interactions if (i.geom.

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!

Jan Stránský (honzik) said : | #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.

i.geom is ScGeom [1]. It does not have 'area' nor 'area1'

What should the condition "if (i.geom.

> 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

[1] https:/

Defri (daredefri) said : | #10 |

> What should the condition "if (i.geom.

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

Jan Stránský (honzik) said : | #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.

> if (i.geom.

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

###

def getBottomStress

# center ... center of the stress circle

# radius ... radius of the circle

interactions = bottom.intrs()

interactions = [i for i in interactions if (i.geom.

fns = [i.phys.normalForce for i in interactions]

fnTot = sum(fn.norm() for fn in fns)

area = pi * radius**2

stress = fnTot / area

return stress

# TODO shear stress?

###

cheers

Jan

> def getBottomStress

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

Defri (daredefri) said : | #13 |

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