Yade

Free fall simulation did not have good results

Asked by Xiaohu Tao

Hello, I am a beginner of Yade. I am trying to couple the dem with LBM (Using Bouzidi boundary condition, some work like Franck Lominé's model). But when I test the free fall simulation, the Yade did not show a good result to me. I'm a little confused.

The following is the code I typed in yade:
---------------------------------------------------------------------------------------------------------------------
    s=utils.sphere(center=[0,0,0],radius=1)

    O.bodies.append(s)

    O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(
    [Ig2_Sphere_Sphere_L3Geom_Inc(),Ig2_Wall_Sphere_L3Geom_Inc()],
    [Ip2_FrictMat_FrictMat_FrictPhys()],
    [Law2_L3Geom_FrictPhys_ElPerfPl()]
    ),

    NewtonIntegrator(damping=0,gravity=(0,0.,-9.81),label="ENewton"),
    ]

    O.dt=1
    O.bodies[0].state.vel #Result: Vector3(0,0,0)
    O.bodies[0].state.pos #Result: Vector3(0,0,0)

    O.step() #1st iter
    O.bodies[0].state.vel #Result: Vector3(0,0,-9.81)
    O.bodies[0].state.pos #Result: Vector3(0,0,-9.81) this doesn't make sense, I think it should be (0,0,-4.905)

    O.step() #2nd iter
    O.bodies[0].state.vel #Result: Vector3(0,0,-19.62)
    O.bodies[0].state.pos #Result: Vector3(0,0,-29.43) this doesn't make sense, I think it should be (0,0,-19.62)

    O.step() #3rd iter
    O.bodies[0].state.vel #Result: Vector3(0,0,-29.43)
    O.bodies[0].state.pos #Result: Vector3(0,0,-58.86) this doesn't make sense, I think it should be (0,0,-44.145)
---------------------------------------------------------------------------------------------

Thanks in advance, I appreciate your help.

Xiaohu

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
Last query:
Last reply:
Revision history for this message
Jérôme Duriez (jduriez) said : 		#1

Hello,

The good thing is that these results seem completely consistent with the doc at https://yade-dem.org/doc/formulation.html#motion-integration.

Motion integration in DEM in general and YADE in particular involves a finite difference scheme, which is necessarily a discrete (in time) approximation of true phenomena being continuous in time.

Thus, like with any other finite difference scheme, the precision of our (classical ?) scheme with respect to theoretical equations directly depends on the size of the time increment (compared with some time scale for the study).

Here, your time scale is pretty much the same as your time increment, hence you can not expect to be satisfied with your results.

On the other hand, you would get a 1% precision if you would do the comparison at t = 100 = 100 * O.dt (or at t = 1, but using O.dt = 0.01)

Note that during "true" YADE simulations (with dynamic bodies), all this discussion would be transparent since the timestep would be drastically reduced for stability reasons.

PS : thank you for your minimal code, note that utils. prefix is unnecessary: you can just type sphere(..) thanks to the way the utils YADE module is imported at startup.
Note also that the L3Geom description of contact is said to be experimental in the doc [*]. I think most users use ScGeom and related classes, and do not think there remains anyone around that could maintain/help with L3Geom.

[*] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.L3Geom

Revision history for this message
Xiaohu Tao (jingje) said : 		#2

Hi, Jérôme

Thank you so much for the detailed explanation.

I read the background of yade. Can I understand in this way? "U" is the position, U(0) is the position at Iter 0, U(1) is the position at Iter 1, and U(2) is the position at Iter 2. "V" is the mean velocity and "a" is the acceleration. in my code a=gravity. △t=1, U(0)=(0,0,0)

So

V(1)=V(0)+a*△t while "V(0)=0", V(1)=(0,0,-9.81)*1=(0,0,-9.81), this is just the mean velocity, NOT the instantaneous velocity at time 0+△t

U(1)=U(0)+V(1)*△t=(0,0,0)+(0,0,-9.81)*1=(0,0,-9.81)

V(2)=V(1)+a*△t =(0,0,-9.81)+(0,0,-9.81)*1=(0,0,-19.62)

U(2)=U(1)+V(2)*△t=(0,0,-9.81)+(0,0,-19.62)=(0,0,-29.43)

------------------------------------------------------------------------------------------------------------------------

In my understanding：(1) the value of V(0)=0 is the key of the error, so the error is formed in the first iter. (2) the velocity showed by yade is the mean velocity in one iter.

Xiaohu

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

I think you correctly understood https://yade-dem.org/doc/formulation.html#motion-integration in terms of notation, except we use funny symbols (+/- with a circle around, instead of just 0,+,-) for velocity because you have to consider the velocity as being a value *at the middle* of a time step.
See the explanations in "Let us recall..." of https://yade-dem.org/doc/formulation.html#motion-integration

You can call it "mean velocity" (during one timestep) if you want, this point of view is exact for a linear evolution of the velocity such as here.

Then your V(1) or V(2) are rather V(1/2) and V(3/2) (and V(0) is V(-1/2)..) Indeed, YADE does not give you the instantaneous velocity value at time △t

To conclude with my feelings with respect to your (1)/(2)

(1) I do not see any other error outside the one necessarily introduced by the finite difference approximation

(2) Yes if you want. More precisely it is the value one half time step earlier (this "circled minus superscript" value).

Revision history for this message
Xiaohu Tao (jingje) said : 		#4

Hi, Jérôme, thank you so much and I got it. Your patience and explaining is much appreciated.--xiaohu

Revision history for this message
Xiaohu Tao (jingje) said : 		#5

Thanks Jérôme Duriez, that solved my question.