Which physical parameters i need and how to add them

Asked by Christian Sommerfeld

Hi,

I want to create a simulation with the basic DEM model from Cundall and Strack. But i don find any calculation for it.
Therefore i have the following question:
Is it correct to input young's modulus, poission's ratio and the friction angle for the law of cundall and strack and the density for the calculation of the gravity force?
I have done it in this way:
O.materials.append(FrictMat(young=5e4,poisson=0.5,frictionAngle=.6,density=2400))

with the following engine:
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
      [Ip2_2xFrictMat_CSPhys()],
      [Law2_Dem3Dof_CSPhys_CundallStrack()]
   ),
   GravityEngine(gravity=(0,0,-9.81)),
 RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(0,0,1),angularVelocity=50*(2*pi/60),ids=quader,label='rotor'),
   NewtonIntegrator(damping=0.4),
 PyRunner(iterPeriod=20,command="myAddPlotData()"),
 VTKRecorder(iterPeriod=100,fileName='/tmp/QuaderAufnahme-',recorders=['spheres','facets'])
]

Furthermore i am still wondering why the simulation only apply the last convention of my material, when i write it in this way:

FacetMat=O.materials.append(FrictMat(young=5e4,poisson=0.5,frictionAngle=.6,density=2400))
SphereMat=O.materials.append(FrictMat(young=5e10,poisson=0.5,frictionAngle=.6,density=10000))

oriBody = Quaternion(Vector3(0,0,1),(math.pi))
quader=O.bodies.append(geom.facetBox((.5,0,.5),(.3,.3,.3),oriBody, wallMask=63, color=(0,0,1), dynamic=True, material=FacetMat))

sphereRadius = 0.06
nbSpheres = (10,10,30)
for i in xrange(nbSpheres[0]):
    for j in xrange(nbSpheres[1]):
        for k in xrange(nbSpheres[2]):
            x = (i*2 - nbSpheres[0])*sphereRadius*1.1
            y = (j*2 - nbSpheres[1])*sphereRadius*1.1
            z = (k*2 - nbSpheres[2])*sphereRadius*1.1
            sp1=utils.sphere([x,y,z+10],sphereRadius,color= (1,0,0),material=SphereMat)
            O.bodies.append(sp1)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Christian Sommerfeld
Solved:
Last query:
Last reply:

This question was reopened

Revision history for this message
Chiara Modenese (chiara-modenese) said :
#1

Hi Christian,

Could you please (briefly) explain here what material you aim to simulate and what kind of test you want to run with it? It would help us to address your problem.

Chiara

On 20 Jul 2012, at 11:51, Christian Sommerfeld wrote:

> New question #203709 on Yade:
> https://answers.launchpad.net/yade/+question/203709
>
> Hi,
>
> I want to create a simulation with the basic DEM model from Cundall and Strack. But i don find any calculation for it.
> Therefore i have the following question:
> Is it correct to input young's modulus, poission's ratio and the friction angle for the law of cundall and strack and the density for the calculation of the gravity force?
> I have done it in this way:
> O.materials.append(FrictMat(young=5e4,poisson=0.5,frictionAngle=.6,density=2400))
>
> with the following engine:
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
> InteractionLoop(
> [Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
> [Ip2_2xFrictMat_CSPhys()],
> [Law2_Dem3Dof_CSPhys_CundallStrack()]
> ),
> GravityEngine(gravity=(0,0,-9.81)),
> RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(0,0,1),angularVelocity=50*(2*pi/60),ids=quader,label='rotor'),
> NewtonIntegrator(damping=0.4),
> PyRunner(iterPeriod=20,command="myAddPlotData()"),
> VTKRecorder(iterPeriod=100,fileName='/tmp/QuaderAufnahme-',recorders=['spheres','facets'])
> ]
>
> Furthermore i am still wondering why the simulation only apply the last convention of my material, when i write it in this way:
>
> FacetMat=O.materials.append(FrictMat(young=5e4,poisson=0.5,frictionAngle=.6,density=2400))
> SphereMat=O.materials.append(FrictMat(young=5e10,poisson=0.5,frictionAngle=.6,density=10000))
>
> oriBody = Quaternion(Vector3(0,0,1),(math.pi))
> quader=O.bodies.append(geom.facetBox((.5,0,.5),(.3,.3,.3),oriBody, wallMask=63, color=(0,0,1), dynamic=True, material=FacetMat))
>
> sphereRadius = 0.06
> nbSpheres = (10,10,30)
> for i in xrange(nbSpheres[0]):
> for j in xrange(nbSpheres[1]):
> for k in xrange(nbSpheres[2]):
> x = (i*2 - nbSpheres[0])*sphereRadius*1.1
> y = (j*2 - nbSpheres[1])*sphereRadius*1.1
> z = (k*2 - nbSpheres[2])*sphereRadius*1.1
> sp1=utils.sphere([x,y,z+10],sphereRadius,color= (1,0,0),material=SphereMat)
> O.bodies.append(sp1)
>
>
>
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

Revision history for this message
Jérôme Duriez (jduriez) said :
#2

> Is it correct to input young's modulus, poission's ratio and the friction angle for the law of cundall and strack and the density for the calculation of the gravity force?

Yes, not only correct, but also required. The law you plan to use : Law2_Dem3Dof_CSPhys_CundallStrack, requires interaction data of "CSPhys" type (see the name). Data that are computed "FrictMat" type (see the previous "engine" Ip2_2xFrictMat_CSPhys() ).

Hence you need to define all data for FrictMat, otherwise default values will be used.

Note that, for simulations with "Cundall-Strack" law, there is also "Law2_ScGeom_FrictPhys_CundallStrack", for which you should probably find several examples e.g. in the scripts.

Revision history for this message
Jérôme Duriez (jduriez) said :
#3

Read "Data that are computed FROM MATERIAL DATA OF "FrictMat" type" 3rd line of my answer....

Revision history for this message
Chiara Modenese (chiara-modenese) said :
#4

Sorry, I get now the type of question. Christian, follow Jerome's suggestions, you should be fine then. Chiara

On 20 Jul 2012, at 13:11, jduriez wrote:

> Question #203709 on Yade changed:
> https://answers.launchpad.net/yade/+question/203709
>
> jduriez proposed the following answer:
>> Is it correct to input young's modulus, poission's ratio and the
> friction angle for the law of cundall and strack and the density for the
> calculation of the gravity force?
>
> Yes, not only correct, but also required. The law you plan to use :
> Law2_Dem3Dof_CSPhys_CundallStrack, requires interaction data of "CSPhys"
> type (see the name). Data that are computed "FrictMat" type (see the
> previous "engine" Ip2_2xFrictMat_CSPhys() ).
>
> Hence you need to define all data for FrictMat, otherwise default values
> will be used.
>
> Note that, for simulations with "Cundall-Strack" law, there is also
> "Law2_ScGeom_FrictPhys_CundallStrack", for which you should probably
> find several examples e.g. in the scripts.
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

Revision history for this message
Christian Sommerfeld (sommerfeld-m) said :
#5

Hi Chiara,

for the spheres we want to use a ceramic material and for the quader/box we want to use steel as it's material. First of all we want to simulate a elastic behaviour, therefore i chose FricMat as the best Material to handle this problem. Furthermore we want to use create a simulation with the law of cundall and strack to compare it afterwards with other law, which are maybe more suitable.

It is not such a difficult test. The spheres and the box are inside a cylinder, whereby there are so many spheres that they enclose the box. And the box executes a rotational movement around the z-axes and is also fixed. Now we want to measure the agent force on the box.

Our programm is still working, but i am not sure if i take the correct physical parameters for this law, because i nly found the calculation for kN and kS in the PhD of Smilauer but there is no reference to a law.

And i can not add different materials to the spheres and the box.

My programm is the following:

# import yade modules that we will use below
from yade import pack, plot

# PhysicalParameters

FacetMat=O.materials.append(FrictMat(young=5e4,poisson=0.5,frictionAngle=.6,density=2400))

oriBody = Quaternion(Vector3(0,0,1),(math.pi))
quader=O.bodies.append(geom.facetBox((.5,0,.5),(.3,.3,.3),oriBody, wallMask=63, color=(0,0,1), dynamic=True, material=FacetMat))

# show all bodies of the facet
for l in O.bodies:
 g=l.id
 print g

# create Cylinder from facets
cyl=O.bodies.append(utils.geom.facetCylinder((0,0,1),1.5,2,orientation=Quaternion((1, 0, 0),0),segmentsNumber=20,wallMask=6, color=(0,1,0), fixed=True, material=FacetMat))

SphereMat=O.materials.append(FrictMat(young=5e7,poisson=0.5,frictionAngle=.6,density=1000))

# create block of spheres
sphereRadius = 0.06
nbSpheres = (10,10,30)
#nbSpheres = (50,50,50)
for i in xrange(nbSpheres[0]):
    for j in xrange(nbSpheres[1]):
        for k in xrange(nbSpheres[2]):
            x = (i*2 - nbSpheres[0])*sphereRadius*1.1
            y = (j*2 - nbSpheres[1])*sphereRadius*1.1
            z = (k*2 - nbSpheres[2])*sphereRadius*1.1
            sp1=utils.sphere([x,y,z+10],sphereRadius,color= (1,0,0),material=SphereMat)
            O.bodies.append(sp1)

sphereRadius = 0.06
nbSpheres = (10,10,30)
#nbSpheres = (50,50,50)
for i in xrange(nbSpheres[0]):
    for j in xrange(nbSpheres[1]):
        for k in xrange(nbSpheres[2]):
            x = (i*2 - nbSpheres[0])*sphereRadius*1.1
            y = (j*2 - nbSpheres[1])*sphereRadius*1.1
            z = (k*2 - nbSpheres[2])*sphereRadius*1.1
            sp2=utils.sphere([x,y,z+4],sphereRadius,color=(1,0,0),material=SphereMat)
            O.bodies.append(sp2)

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      # handle sphere+sphere and facet+sphere collisions
      [Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
      #[Ip2_FrictMat_FrictMat_FrictPhys()],
      [Ip2_2xFrictMat_CSPhys()],
      [Law2_Dem3Dof_CSPhys_CundallStrack()]
      #[Law2_L3Geom_FrictPhys_ElPerfPl()]
   ),
   GravityEngine(gravity=(0,0,-9.81)),
 RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(0,0,1),angularVelocity=50*(2*pi/60),ids=quader,label='rotor'),
   NewtonIntegrator(damping=0.4),
 PyRunner(iterPeriod=20,command="myAddPlotData()"),
 VTKRecorder(iterPeriod=100,fileName='/tmp/QuaderAufnahme-',recorders=['spheres','facets'])
]

O.dt=1*utils.PWaveTimeStep()

def myAddPlotData():
 force=plot.addData(t=O.time,F=O.forces.f(0).norm()+O.forces.f(1).norm()+O.forces.f(2).norm()+O.forces.f(3).norm()+O.forces.f(4).norm()+O.forces.f(5).norm()+O.forces.f(6).norm()+O.forces.f(7).norm()+O.forces.f(8).norm()+O.forces.f(9).norm()+O.forces.f(10).norm()+O.forces.f(11).norm())

plot.plots={'t':('F')}
plot.plot()

Revision history for this message
Chiara Modenese (chiara-modenese) said :
#6

Hi Christian,

Thanks for details.

You can attribute whatever material you want to whatever body. For instance:

# create two materials
O.materials.append(..., label='mat_spheres');
O.materials.append(..., label='mat_box');

# assign materials to different bodies
O.bodies.append(..., material='mat_spheres')
O.bodies.append(..., material='mat_box')

Please see: https://yade-dem.org/doc/yade.utils.html?highlight=utils.sphere#yade.utils.sphere
You do something similar in your code (but using the material instance instead of a label like I do), it should work too.

As for the choice of your input parameters, that is never easy to do. To start with, you could take some initial values from the literature for a similar material. Then you can carry out a parametric study to assess the influence of stiffness on your test, for instance.

Good luck with your test,
Chiara

--
Chiara Modenese, BSc MSc
DPhil(PhD) Candidate in Engineering Science
Department of Engineering Science, University of Oxford
Parks Road, Oxford, OX1 3PJ, UK
email: <email address hidden>

On 20 Jul 2012, at 13:35, Christian Sommerfeld wrote:

> Question #203709 on Yade changed:
> https://answers.launchpad.net/yade/+question/203709
>
> Christian Sommerfeld posted a new comment:
> Hi Chiara,
>
> for the spheres we want to use a ceramic material and for the quader/box
> we want to use steel as it's material. First of all we want to simulate
> a elastic behaviour, therefore i chose FricMat as the best Material to
> handle this problem. Furthermore we want to use create a simulation
> with the law of cundall and strack to compare it afterwards with other
> law, which are maybe more suitable.
>
> It is not such a difficult test. The spheres and the box are inside a
> cylinder, whereby there are so many spheres that they enclose the box.
> And the box executes a rotational movement around the z-axes and is also
> fixed. Now we want to measure the agent force on the box.
>
> Our programm is still working, but i am not sure if i take the correct
> physical parameters for this law, because i nly found the calculation
> for kN and kS in the PhD of Smilauer but there is no reference to a law.
>
> And i can not add different materials to the spheres and the box.
>
> My programm is the following:
>
>
> # import yade modules that we will use below
> from yade import pack, plot
>
> # PhysicalParameters
>
>
> FacetMat=O.materials.append(FrictMat(young=5e4,poisson=0.5,frictionAngle=.6,density=2400))
>
>
>
> oriBody = Quaternion(Vector3(0,0,1),(math.pi))
> quader=O.bodies.append(geom.facetBox((.5,0,.5),(.3,.3,.3),oriBody, wallMask=63, color=(0,0,1), dynamic=True, material=FacetMat))
>
> # show all bodies of the facet
> for l in O.bodies:
> g=l.id
> print g
>
> # create Cylinder from facets
> cyl=O.bodies.append(utils.geom.facetCylinder((0,0,1),1.5,2,orientation=Quaternion((1, 0, 0),0),segmentsNumber=20,wallMask=6, color=(0,1,0), fixed=True, material=FacetMat))
>
>
> SphereMat=O.materials.append(FrictMat(young=5e7,poisson=0.5,frictionAngle=.6,density=1000))
>
> # create block of spheres
> sphereRadius = 0.06
> nbSpheres = (10,10,30)
> #nbSpheres = (50,50,50)
> for i in xrange(nbSpheres[0]):
> for j in xrange(nbSpheres[1]):
> for k in xrange(nbSpheres[2]):
> x = (i*2 - nbSpheres[0])*sphereRadius*1.1
> y = (j*2 - nbSpheres[1])*sphereRadius*1.1
> z = (k*2 - nbSpheres[2])*sphereRadius*1.1
> sp1=utils.sphere([x,y,z+10],sphereRadius,color= (1,0,0),material=SphereMat)
> O.bodies.append(sp1)
>
> sphereRadius = 0.06
> nbSpheres = (10,10,30)
> #nbSpheres = (50,50,50)
> for i in xrange(nbSpheres[0]):
> for j in xrange(nbSpheres[1]):
> for k in xrange(nbSpheres[2]):
> x = (i*2 - nbSpheres[0])*sphereRadius*1.1
> y = (j*2 - nbSpheres[1])*sphereRadius*1.1
> z = (k*2 - nbSpheres[2])*sphereRadius*1.1
> sp2=utils.sphere([x,y,z+4],sphereRadius,color=(1,0,0),material=SphereMat)
> O.bodies.append(sp2)
>
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
> InteractionLoop(
> # handle sphere+sphere and facet+sphere collisions
> [Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
> #[Ip2_FrictMat_FrictMat_FrictPhys()],
> [Ip2_2xFrictMat_CSPhys()],
> [Law2_Dem3Dof_CSPhys_CundallStrack()]
> #[Law2_L3Geom_FrictPhys_ElPerfPl()]
> ),
> GravityEngine(gravity=(0,0,-9.81)),
> RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(0,0,1),angularVelocity=50*(2*pi/60),ids=quader,label='rotor'),
> NewtonIntegrator(damping=0.4),
> PyRunner(iterPeriod=20,command="myAddPlotData()"),
> VTKRecorder(iterPeriod=100,fileName='/tmp/QuaderAufnahme-',recorders=['spheres','facets'])
> ]
>
> O.dt=1*utils.PWaveTimeStep()
>
> def myAddPlotData():
> force=plot.addData(t=O.time,F=O.forces.f(0).norm()+O.forces.f(1).norm()+O.forces.f(2).norm()+O.forces.f(3).norm()+O.forces.f(4).norm()+O.forces.f(5).norm()+O.forces.f(6).norm()+O.forces.f(7).norm()+O.forces.f(8).norm()+O.forces.f(9).norm()+O.forces.f(10).norm()+O.forces.f(11).norm())
>
> plot.plots={'t':('F')}
> plot.plot()
>
> --
> 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 address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp

Revision history for this message
Christian Sommerfeld (sommerfeld-m) said :
#7

Hi,

first of all thank you for your response!

But i have a last question to Jerome's response. I found the other basic law for DEM in the user's manual as well.

InteractionLoop(
   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
   [Ip2_FrictMat_FrictMat_FrictPhys()],
   [Law2_ScGeom_FrictPhys_CundallStrack()]
)

and with regard to the law which i chose

InteractionLoop(
      [Ig2_Sphere_Sphere_Dem3DofGeom(),Ig2_Facet_Sphere_Dem3DofGeom()],
      [Ip2_2xFrictMat_CSPhys()],
      [Law2_Dem3Dof_CSPhys_CundallStrack()]
   )

i am a little bit confused!
Because i understand that the ScGeom is only useable for a Sphere and a non-spherical body and not if i want to determinate the force between to spheres? Furthermore i don't know what is the different between CSPhys and FrictPhys, because i don't find any description about CSPhys ?

Revision history for this message
Jérôme Duriez (jduriez) said :
#8

Hello,

No fear to have about ScGeom, it is useable for both (sphere-sphere) and (sphere-non sphere) ! But it's true I would not have understood this only whith the doc https://www.yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ScGeom. Do not hesitate to improve it !

And according to the documentation (https://www.yade-dem.org/doc/yade.wrapper.html#iphys), CSPhys indeed does not exist... Is this script, where you found it, working ? Open a bug (concerning the script or the doc) anyway !

Revision history for this message
Christian Sommerfeld (sommerfeld-m) said :
#9

Hi Jerome,

thank you for your advice! I will try it with the ScGeom and the suggested engine for the law of cundall and strack.

My script is still working with the CSPhys and i have seen this kind of physics in the documentation a few weeks ago.

Thanks a lot for all your advices!

Cheers,
Christian

Revision history for this message
Christian Sommerfeld (sommerfeld-m) said :
#10

Oh sorry i forgot to ask you what do you mean with "Hence you need to define all data for FrictMat"? Does it mean that i have to define G_over_E and the coeff_dech for example?

Sorry for that stupid question!

Revision history for this message
Jérôme Duriez (jduriez) said :
#11

No, I meant "all variables belonging to FrictMat (unless you precisely know what you're doing)". So density, young, poisson, and friction angle and that's all.

 These variables "G_over_E", or "coeff_dech" for example belong to other material types (derived from FrictMat).

Revision history for this message
Christian Sommerfeld (sommerfeld-m) said :
#12

Oh yeah i only use this kinds of parameters. I was only a liitle bit confused about your answer, sorry for that!

Now me programm is still working fine, but the ubuntu terminal told me that kind of error:

ERROR /build/buildd/yade-daily-1+3049+27~oneiric1/pkg/common/GravityEngines.cpp:19 action: GravityEngine is deprecated, consider using Newton::gravity instead (unless gravitational energy has to be tracked - not implemented with the newton attribute)

when i use

GravityEngine(gravity=[0,0,-9.81])

is that a problem, because they use the gravital engine in the same way in the examples?

Revision history for this message
Luc Scholtès (luc) said :
#13

Dear Christian,

As mentionned by yade (in your terminal), GravityEngine is now deprecated and you should use NewtonIntegrator to apply gravity in your simulation like this:

NewtonIntegrator(damping=yourValueOfDamping, gravity=(0,yourValueOfGravity,0))

The examples using GravityEngine are deprecated, but everything should work with GravityEngine anyway.

Cheers

  Luc

Revision history for this message
Christian Sommerfeld (sommerfeld-m) said :
#14

Hi Guys,

thank you very much for your advices!

Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#15

@Jerome
>I would not have understood this only whith the doc

Thanks for seeing that. The documentation was correct in the source but it was compiled incorrectly and diplayed with missing sentences. I'll fix that. The original text reads (ScGeom.hpp:53):

"Class representing geometry of a contact point between two bodies. It is more general than sphere-sphere contact even though it is primarily focused on spheres interactions (reason for the 'Sc' namming); it is also used for representing contacts of a Sphere with non-spherical bodies (Facet, Plane, Box, ChainedCylinder)"