Sphere goes through the facet wall in Harmonic vibration

Asked by gaoxuesong

Hi, in my case there are lots of spheres in one box and a piston in the bottom of the box. Then the piston does harmonic vibration. I find when the vibration amplitude goes to high, the spheres go through the piston, which is constructed by facet. I want to know the reason and fix it. In addition, the density is amplified by a factor of 1e9 to increase the calculation speed. And also I want to know how to define the word of "rigid" in yade. In my case, the box is rigid and the spheres are relatively soft. Does it make sense if i set the young's modulus of the box material as 100 times of that of the sphere material?

Following is my code,

# PhysicalParameters
muS = 0.57735
muF = 0.17632
FricAngleS = math.atan(muS)
FricAngleF = math.atan(muF)
densitys_ac = 8150*1e9
densityw_ac = 2700*1e9
yongmodu = 195e9
poisn = 0.3

# material model
matSph = CohFrictMat(density=densitys_ac, young=yongmodu, poisson=poisn, frictionAngle=FricAngleS, momentRotationLaw=True)
SMat = O.materials.append(matSph)
matFacet = CohFrictMat(density=densityw_ac, young=yongmodu, poisson=poisn, frictionAngle=FricAngleF, momentRotationLaw=True)
FMat = O.materials.append(matFacet)

## box and piston
O.bodies.append(geom.facetBox((0.5e-3,0.5e-3,0.5e-3),(0.5e-3,0.5e-3,0.5e-3), wallMask=31, material=FMat))
Cylinder1IDs=O.bodies.append(geom.facetBox((0.5e-3,0.5e-3,0.15e-3),(0.5e-3,0.5e-3,0.2e-3), wallMask=63, color=(0,1,0), wire=False, material=FMat))

# spreading process motions first harmonic vibratiion then stop
def change_motion():
    if O.time > 0.4:
        harmEngineP1.dead = False
        harmEngineP1.A = (0.0,0.0,0.35e-3*0.25)
    elif O.time > 0.5:
        harmEngineP1.A = (0, 0, 0)
        harmEngineP1.f = (0, 0, 0)

#define engines:
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
       # handle sphere+sphere and facet+sphere collisions
         [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
         [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
          [Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
    ),
    PyRunner(iterPeriod=100,command='change_motion()'),
    HarmonicMotionEngine(dead=True,label='harmEngineP1', A=(0.0,0.0,1.5e-5), f=(0.0,0.0,60.0), ids=Cylinder1IDs),
    NewtonIntegrator(damping=0.75, exactAsphericalRot=True, gravity=(0,0,-9.81)),
# VTKRecorder(iterPeriod=2000, recorders=['spheres','colors'], fileName='./vtu/t3-SS316L_60um.vtu'),
]
## generation particles
sp = pack.SpherePack()
sp.makeCloud((0,0,0.35e-3), (1.0e-3,1.0e-3,1.2e-3), rMean=40e-6, rRelFuzz=40e-6)
sp.toSimulation(material=SMat)
## time step
O.dt = 0.85*utils.PWaveTimeStep()
O.usesTimeStepper = False

O.run()

Thanks,
Xuesong

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
gaoxuesong
Solved:
Last query:
Last reply:
Revision history for this message
Jan Stránský (honzik) said :
#1

Hello,

> Then the piston does harmonic vibration. I find when the vibration amplitude goes to high, the spheres go through the piston, which is constructed by facet. I want to know the reason and fix it.

There might be several issues:
a) the piston period is to short w.r.t O.dt that the piston simply "overjumps" the particles, not interacting with them at all
b) the stiffness/mass/timestep combination is such that the repulsive force acts so little time, has so liitle magnitude etc. that the piston is not able to "push away" the particles.
c) Linear contact law also does not help b) as for really high penetration, the force is "too low"

> In addition, the density is amplified by a factor of 1e9 to increase the calculation speed.

by increasing density, you can get higher critical time step. BUT you also change the physics of the system (see point b above).
It always depends if and which values are acceptable. 1e9 * density sounds really too much..
If you want to get "realistic" results of a dynamic system, **in my opinion** this is not the way to go..

Somewhat less drastic scenario is to increase mass of some fraction (e.g. 5% of mass) of the smallest particles.
It increases time step (determined by the smallest particles), but preserves dynamic properties of the "most" of the system.

> And also I want to know how to define the word of "rigid" in yade. In my case, the box is rigid and the spheres are relatively soft.
> Does it make sense if i set the young's modulus of the box material as 100 times of that of the sphere material?

the stiffness of interaction is computed as a harmonic average of stiffness of individual bodies, something like
E=2*E1*E2/(E1+E2)
The limit for E2->infinity is E=2*E1. For E2=100*E1, E=1.98, so it could be "close enough" to rigid.

cheers
Jan

Revision history for this message
gaoxuesong (260582472-9) said :
#2

Hi Jan,

Thanks for your answers. For the density magnification, in my case, the diameter of the sphere is around 80 um. If i use the normal density, it costs about 4 hours to simulate virtual time of 55 ms. The total time i need is 3s, so it is very time-consuming. In my method, the last step is mass shrinking, where the mass of all the particles is gradually decreased to the 1/1000 of the used value. Does it make sense or do you have good advice?

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

> it is very time-consuming

yes, sometimes / usually DEM is like that :-)

> the diameter of the sphere is around 80 um
> sp.makeCloud((0,0,0.35e-3), (1.0e-3,1.0e-3,1.2e-3), rMean=40e-6, rRelFuzz=40e-6)

the problem is that with rMean=rRelFuzz, the radius may be arbitrarily small (you can print the minimal value), making the PWaveTimeStep very small, too. Thereofore I suggested to change mass only by the smallest particles (but with manually modified mass, you cannot use PWaveTimeStep but a "hand made" estimation)

> the last step is mass shrinking, where the mass of all the particles is gradually decreased to the 1/1000 of the used value. Does it make sense or do you have good advice?

sorry, I did not get this point at all..

cheers
Jan

Revision history for this message
gaoxuesong (260582472-9) said :
#4

>the problem is that with rMean=rRelFuzz, the radius may be arbitrarily small (you can print the minimal value), making the PWaveTimeStep very small, too.

Actually i run an example case where all the particles have the same diameter of 100 um and the time step is set as 1e7. As soon as the simulation starts, the particles explode and flow out.

> the last step is mass shrinking, where the mass of all the particles is gradually decreased to the 1/1000 of the used value. Does it make sense or do you have good advice?

Sorry for confusing words used. The following is an example,

### a function to reduce the mass gradually
redt = 0.5 ### how long for one mass reduction
retime = 1.0 ### the time when mass reduction happens
mrecof = 0.1 ## the mass reduction cofficient
def mass_relaxation():
    global redt
    global retime
    if O.time > retime:
        retime = retime + redt
        for eb in O.bodies:
            if isinstance(eb.shape, Sphere):
                eb.state.mass = eb.state.mass*mrecof
        print('change mass', O.time)

Then the function of mass_relaxation() is invoked by pyRunner Engine periodically.

Thanks.

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

Hi
>the density is amplified by a factor of 1e9 to increase the calculation speed

It will not work. With such mass the particles will be nearly immobile on short time scale, and that's probably why the piston goes through.
Or you would have to scale down frequency by the sqrt of this multiplier, hence exactly the same number of time iterations in total.
Bruno

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

> time step is set as 1e7. As soon as the simulation starts, the particles explode and flow out.

it seems like the simulation is unstable = has too large time step. What is the value of PWaveTimeStep?

Concerning mass reduction, what does "last step" in "the last step is mass shrinking" mean? At the end of simulation? Now I understand what you mean, but don't see the point / advantage fro its usage..

As said, any material parameters modifications "to make the simulation faster" changes the physics. If it is suitable is very individual and there is no general answer if it is good or bad.

cheers
Jan

Revision history for this message
gaoxuesong (260582472-9) said :
#7

> the time step
    The PWaveTimeStep is 1.02e-08.
    By the way, the default option of Yade is in parallel and uses all of the available cores? In the tutorial, it says, "By default, each job uses all available cores for itself, which causes jobs to be effectively run in parallel. Number of cores per job can be globally changed via the --job-threads option"
   If i just start a job by "yade pyname.py", is it executed in parallel and uses the available cores? If so, how many cores does it use when i start a new job without designating the cores to use? To grab some cores from the first job?

> the meaning of the mass reduction
    For my case, the packing density is an important factor to affect the following process. Larger mass makes denser packing. So this is the meaning of the step. Yes. I agree with you the mass modification changes the physical process. But if we actually care about the random feature of particle packing not the exact the particle distribution.

Thanks,
Xuesong

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

> The PWaveTimeStep is 1.02e-08.
> time step is set as 1e7

1e7 or 1e-7? in any case, it is > PWaveTimeStep and most likely makes the simulation unstable with the consequences you see

> By the way, the default option of Yade is in parallel and uses all of the available cores? In the tutorial, it says, "By default, each job uses all available cores for itself, which causes jobs to be effectively run in parallel. Number of cores per job can be globally changed via the --job-threads option"
> If i just start a job by "yade pyname.py", is it executed in parallel and uses the available cores? If so, how many cores does it use when i start a new job without designating the cores to use? To grab some cores from the first job?

This is unrelated to the original problem, next time please open a new question for a new question [1]
Also next time please refer "the tutorial" etc. I see it in User's manual, being a different part of documentation than tutorial.

"all available cores" you are referring holds for batch execution (yade-batch).
yade pyname.py runs by default using 1 core.

cheers
Jan

[1] https://yade-dem.org/wiki/Howtoask

Revision history for this message
gaoxuesong (260582472-9) said :
#9

Thanks, my question has been solved.