Not sure about utils.PWaveTimeStep()

Asked by Christian Sommerfeld

Hi,

i'm not sure if i have to use the utils.PWaveTimeStep() in my case. I want to simulate a facetbox which is rotating around the z-axes and enclosed by a pack of spheres . Therefore i want to record the force which acts on the facetbox.
The ducumentation says that the utils.PWaveTimeStep() makes some simplifications, like
• all particles are spherical and have the same radius R;
• the sphere’s material has the same E and ρ

All of my spheres have the same radius, but i use a facet as well. And i have seen that the utils.PWaveTimeStep() in some examples with facets and spheres. So i'm still wondering if i can use the utils.PWaveTimeStep() in my case? Or should i use GlobalStiffnessTimeStepper because i use a facet?

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

Hello Christian,

> i'm not sure if i have to use the utils.PWaveTimeStep() in my case. I want
> to simulate a facetbox which is rotating around the z-axes and enclosed by
> a pack of spheres . Therefore i want to record the force which acts on the
> facetbox.
>

if you use non-dynamic facets, then it is ok to use utils.PWaveTimeStep.
The critical time step is computed for minimal possible (reasonable)
interaction length (i.e. sphere radius), so also interaction sphere-facet
(where the interaction length should equals radius) should behave stable
for such time step.

The requirement of spherical bodies is due to internal implementation..

> The ducumentation says that the utils.PWaveTimeStep() makes some
> simplifications, like
> * all particles are spherical and have the same radius R;
> * the sphere's material has the same E and ρ
>
>
where did you find these information? I was not able to find it :-)

> All of my spheres have the same radius, but i use a facet as well. And i
> have seen that the utils.PWaveTimeStep() in some examples with facets and
> spheres. So i'm still wondering if i can use the utils.PWaveTimeStep() in
> my case? Or should i use GlobalStiffnessTimeStepper because i use a facet?
>

then it depends on the speed of your facets. If for example for the wave
critical time step the facets would move through 3 spheres, then you should
reduce it. But GlobalStiffnessStepper would not help in this case (as from
the stiffness point of view the pwave time step will be ok)

I hope I formulated the answer in understandable way, if not, don't
hesitate to ask again :-)

Jan

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

Hi Christian,

As long as the facets are not dynamic, you don't need to consider facets
in timestep definition.

PWaveTimeStep can't be directly compared to GlobalStiffnessTimeStepper
as you.
You can consider PWaveTimeStep as an order of magnitude of the expected
critical timestep (it at least scales correctly with sizes, modulus,
density). Since it is only an order of magnitude, it's usually
multiplied by a safety factor after trial and error for each type of
simulation.

GlobalStiffnessTimeStep computes the exact critical timestep for current
configuration, but it gives nothing as long as there is no contact, so
that it can't be used to define dt in the early steps of most simulations.

So there is not to choose between them, most scripts use both: the first
for early steps, the second to set dt to its maximum acceptable value in
the rest of the simulation (hence save time). Both work equally well
with facets.

Last note: nothing can replace your own tests. It is not uncommon to
read a "no it doesn't work" answer in this list when the thing actually
works (Yade code is so huge that nobody knows the entire code now). So,
when you wonder if something works, the best is to try it before asking.

Bruno

On 24/07/12 12:45, Christian Sommerfeld wrote:
> New question #204014 on Yade:
> https://answers.launchpad.net/yade/+question/204014
>
> Hi,
>
> i'm not sure if i have to use the utils.PWaveTimeStep() in my case. I want to simulate a facetbox which is rotating around the z-axes and enclosed by a pack of spheres . Therefore i want to record the force which acts on the facetbox.
> The ducumentation says that the utils.PWaveTimeStep() makes some simplifications, like
> • all particles are spherical and have the same radius R;
> • the sphere’s material has the same E and ρ
>
> All of my spheres have the same radius, but i use a facet as well. And i have seen that the utils.PWaveTimeStep() in some examples with facets and spheres. So i'm still wondering if i can use the utils.PWaveTimeStep() in my case? Or should i use GlobalStiffnessTimeStepper because i use a facet?
>
>

--
_______________
Bruno Chareyre
Associate Professor
ENSE³ - Grenoble INP
Lab. 3SR
BP 53
38041 Grenoble cedex 9
Tél : +33 4 56 52 86 21
Fax : +33 4 76 82 70 43
________________

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

Hi,

thank you guys!

@Jan: I found it in the DEM-Formulation at the point "Estimation of \Dtcr by wave propagation speed"

My programm is still working but i'm not sure if i get the correct results.

So now i understand that i have to use both (PWaveTimeStep and PWaveTimeStep) in my case. Because my facet is dynamic.
I can give a short summary of my simulation to make it more understandable:

I created a ball and a box which are located in a cylinder and the box and the box rotates around the z axis enclosed by the spheres. That's what i said before.

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

cylinder=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))

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='mat_spheres')
            O.bodies.append(sp1)

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]

   ),
  NewtonIntegrator(damping=0.4, gravity=(0,0,-9.81)),
 RotationEngine(rotateAroundZero=True,zeroPoint=(0,0,0),rotationAxis=(0,0,1),angularVelocity=100*(2*pi/60),ids=box,label='rotor')
]

O.dt=0.2*utils.PWaveTimeStep()

Should i add the GlobalStiffnessTimeStepper to my code?

Christian

Revision history for this message
Jan Stránský (honzik) said :
#4

Hi Christian,

I was only checking the utils module documentation, thanks for the
location.

your facets are not actually dynamic (as far as I see from the script),
because you apply RotationEngine on box and thus prescribe all their
motion, while cylinder is fixed from definition. So from this point of view
PVaweTimeStep is ok.

Anyway, it is always a good idea to use GlobalTimeStepper at least once in
the beginning of your simulation (but, as Brune mentioned, after a few
steps), because PVaweTimeStep computes time step from "young" of spheres,
but each material law handles it in different way (this is also mentioned
in DEM formulation in Estimade dT_cr by vawe propagation), so it might be a
bit dangerous - but you will definitely know it if it was your case :-)

to conclude, use both PVaweTimeStepper and GlobalTimeStepper if you are not
sure about your time step :-) one other option is to set dT manually
O.dt = sqrt( (density * 4/3.*pi*r*r) / young ) # for CpmMat_CpmPhys, in
your case (FrictMat_FrictPhys) it would be different :-)

Jan

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

Hi,

thank you very much! So i will try it first with the PWaveTimeStep and add PWaveTimeStep later to compare the results.

But now i'm interested in the option to set dt manually, where can i find some advices to do it in this way? Because i didn't find so many formulas for FrictMat_FrictPhys in the documentation to calculate different things for it!

Christian

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

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

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

Hi,

thank you very much! So i will try it first with the PWaveTimeStep and add PWaveTimeStep later to compare the results.

But now i'm interested in the option to set dt manually, where can i find some advices to do it in this way? Because i didn't find so many formulas for FrictMat_FrictPhys in the documentation to calculate different things for it!

Christian

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

A general rule in explicit DEM:
- If dt is too large, the scheme is unstable. Simulation explodes in a
few steps, it is impossible to not notice the problem.
- If the scheme is stable, then it is at the same time accurate (i.e. in
most cases you will not notice differences with smaller dt).

Exception: high velocity impacts, maybe your case if the rotational
speed of the facets is high. In this case, it can be stable BUT innacurate.
Depends on your rotational velocity, none of the timesteppers will help
to decide in this situation (well, they still give an upper bound for dt).

Bruno

Revision history for this message
Jan Stránský (honzik) said :
#9

>
> thank you very much! So i will try it first with the PWaveTimeStep and
> add PWaveTimeStep later to compare the results.
>
> But now i'm interested in the option to set dt manually, where can i
> find some advices to do it in this way? Because i didn't find so many
> formulas for FrictMat_FrictPhys in the documentation to calculate
> different things for it!
>

For the case of FrictMat_FrictPhys it is mentioned in [1]. The critical
time step for discrete systems (mass spring systems), which Yade actually
simulate, is given by

dT = sqrt( M / k_N )
where M=density*4/3.*pi*r*r*r is mass of the particle and k_N [N/m]
normal stiffness of contacts / interactions / bonds. Specifically for
Frict_phys it is (according to [1] )
k_N = 2*young*r

good luck and if you have more questions, just ask :-)

Jan

[1]
https://yade-dem.org/doc/formulation.html#estimation-of-by-wave-propagation-speed

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

Thank you guys! I will try it with your advices. My velocity isn't so high. So i think there isn't a problem with a high velocity impact.

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

Hi,

i tried to find all of the calculations in the yade code! And i found a lot of formulas for different calculations. But i don't know exactly where i can find the class Law2_ScGeom_FrictPhys_CundallStrack because there is for example a class Law2_Dem3Dof_CSPhys_CundallStrack in the CundallStrack.hpp/cpp.

Furthermore i don't know where i can find a calculation for the mass from the density value. But the computer have to know how to calculate it.

Christian

Revision history for this message
Jan Stránský (honzik) said :
#12

Hello Christian,

i tried to find all of the calculations in the yade code! And i found a
> lot of formulas for different calculations. But i don't know exactly
> where i can find the class Law2_ScGeom_FrictPhys_CundallStrack because
> there is for example a class Law2_Dem3Dof_CSPhys_CundallStrack in the
> CundallStrack.hpp/cpp.
>
>
pkg/dem/ElasticContactLaw.xpp

hint: you can use linux "grep" command (assuming you are in yade top
directory)
grep -R "class Law2_ScGeom_FrictPhys_CundallStrack" .

> Furthermore i don't know where i can find a calculation for the mass
> from the density value. But the computer have to know how to calculate
> it.
>
>
the mass can be assigned independently on density, but if it is computed
from density, it is done in py/utils.py file, in _commonBodySetup function
(this function is used by utils.sphere or utils.chainedCylinder functions
for example).

hope it will help you
Jan

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

Hello Jan,

sorry for my late response! Thanks a lot for your help!