yade-users team mailing list archive

[Anton Gladky <gladky.anton@xxxxxxxxx>]

Hi Andy,

thanks for raising the question. Sometimes it is important to
check the correctness of basic functionality.

See an attached, slightly modified, script and diagram. I put the first
particle (upper) not at the height 10m, but directly on the
border of the second one, adding an initial velocity,  which
is almost the same, like you have before the impact. Also I have
added the plot function, so you can see all parameters during simulation.

As you can see, there are no attractive forces, if there is no penetration.
The normal force is 0.

Regards


Anton


2014-05-20 21:26 GMT+02:00 Andy Torres <question248996@xxxxxxxxxxxxxxxxxxxxx>:
> Question #248996 on Yade changed:
> https://answers.launchpad.net/yade/+question/248996
>
> Andy Torres posted a new comment:
> Hi Anton,
> thanks for the quick anser. In the code,  there is an "if"
> l.35  if (computeForceTorqueViscEl(_geom, _phys, I, force, torque1, torque2) and (I->isActive)) {...}
> l.44  else {scene->interactions->requestErase(I); return;}
> should this lines remove the interaction? (if so, seems like its not working)
> why YADE keep the interaction when the particles are not in contact?
>
>
> i just try this code:(with yadedaily  Yade 1.07.0-116-c00c469~jessie)
>
> #!/usr/bin/env python
> # -*- coding: utf-8 -*-
> from ya#!/usr/bin/env python
> # -*- coding: utf-8 -*-
> from yade import pack, plot
> import numpy as np
>
> Corr='-test'
>
> Arch1='energ'+Corr+'.out'
> deltaT=5e-5
> grav=-9.81
> zi=10.0
> sumado=Vector3(0.0,0.0,.0)
> # "PARTÍCULAS"
> idmatVisEl=O.materials.append(ViscElMat(kn=1e9,cn=300e4,ks=1.4e4,cs=200e4,frictionAngle=0.36,label='matVisEl'))
>
> EsfMo = O.bodies.append([utils.sphere(center=(0,0,zi),material='matVisEl', radius=1.0,fixed=False)])
> EsFix = O.bodies.append([utils.sphere(center=(0,0,0), material='matVisEl',radius=1.0,fixed=True)])
>
> O.engines=[
>    ForceResetter(),
>    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
>    InteractionLoop(
>       [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
>       [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
>       [Law2_ScGeom_ViscElPhys_Basic()]
>    ),
>         NewtonIntegrator(gravity=(0,0,grav),damping=0.0),
>         PyRunner(command='dissip(sumado)',iterPeriod=1,label='potDisip'),#Cuando comienzo a calcular las interacciones
> ]
> def dissip(sumado):
>         g = open (Arch1,'a')
>         v=O.bodies[0].state.vel[2]
>         z=O.bodies[0].state.pos[2]
>         m=O.bodies[0].state.mass
>         Fn=0.0
>         for i in O.interactions:
>                 PD=i.geom.penetrationDepth
>                 YFn=i.phys.normalForce
>                 Fn=YFn[2]
>                 g.write(str(O.iter)+'\t'+str(z)+'\t'+str(v)+'\t'+'\t'+str(Fn)+'\t'+str(PD)+'\n')
>         g.closed
> O.dt=deltaT
> O.saveTmp()de import pack, plot
> import numpy as np
>
> Corr='-test'
>
> Arch1='energ'+Corr+'.out'
> deltaT=5e-5
> grav=-9.81
> zi=10.0
> sumado=Vector3(0.0,0.0,.0)
> # "PARTÍCULAS"
> idmatVisEl=O.materials.append(ViscElMat(kn=1e9,cn=300e4,ks=1.4e4,cs=200e4,frictionAngle=0.36,label='matVisEl'))
>
> EsfMo = O.bodies.append([utils.sphere(center=(0,0,zi),material='matVisEl', radius=1.0,fixed=False)])
> EsFix = O.bodies.append([utils.sphere(center=(0,0,0), material='matVisEl',radius=1.0,fixed=True)])
>
> O.engines=[
>    ForceResetter(),
>    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
>    InteractionLoop(
>       [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
>       [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
>       [Law2_ScGeom_ViscElPhys_Basic()]
>    ),
>         NewtonIntegrator(gravity=(0,0,grav),damping=0.0),
>         PyRunner(command='dissip(sumado)',iterPeriod=1,label='potDisip'),#Cuando comienzo a calcular las interacciones
> ]
> def dissip(sumado):
>         g = open (Arch1,'a')
>         v=O.bodies[0].state.vel[2]
>         z=O.bodies[0].state.pos[2]
>         m=O.bodies[0].state.mass
>         Fn=0.0
>         for i in O.interactions:
>                 PD=i.geom.penetrationDepth
>                 YFn=i.phys.normalForce
>                 Fn=YFn[2]
>                 g.write(str(O.iter)+'\t'+str(z)+'\t'+str(v)+'\t'+'\t'+str(Fn)+'\t'+str(PD)+'\n')
>         g.closed
> O.dt=deltaT
> O.saveTmp()
>
> and using gnuplot, i plot:
>
> Kn=5e8
> cn=1.5e6
> dt=5e-5
> m=4188.79
> g=9.81
> set term wxt 0
> #set terminal eps
> #set output "NormalForce.eps"
> set title "Penetration Depth, normal Velcity and Force Vs time"
> set xrange[25000:45000]
> set grid
> plot \
>  "energ-test.out" every 15 u 1:($5) ps 0.5 pt 1 title "PD",\
>  "energ-test.out" every 15 u 1:($3/100) ps 0.5 title "normVel/100",\
>  "energ-test.out" every 25 u 1:($4/1e8) ps 0.5 pt 1 title "Fn from yade/1e8",\
>  "energ-test.out" every 25 u 1:((-$5*Kn+$3*cn)/1e8) ps 0.5 pt 4 title "(calculated Fn=-PD*Kn+normVel*Cn)/1e8"
>
> Regards,
> Andy
>
> --
> 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     : yade-users@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users
> More help   : https://help.launchpad.net/ListHelp

Attachment: dia.png
Description: PNG image

Attachment: sim-data_Sphere.data.bz2
Description: BZip2 compressed data

Attachment: sim-data_Sphere_ant.gnuplot
Description: application/gnuplot

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from yade import pack, plot
import numpy as np

Corr='-test'

Arch1='energ'+Corr+'.out'
deltaT=5e-5
grav=-9.81
zi=10.0
sumado=Vector3(0.0,0.0,.0)
# "PARTÍCULAS"
idmatVisEl=O.materials.append(ViscElMat(kn=1e9,cn=300e4,ks=1.4e4,cs=200e4,frictionAngle=0.36,label='matVisEl'))

EsfMo = O.bodies.append(utils.sphere(center=(0,0,2.0),material='matVisEl', radius=1.0,fixed=False))
EsFix = O.bodies.append(utils.sphere(center=(0,0,0), material='matVisEl',radius=1.0,fixed=True))

O.bodies[EsfMo].state.vel=Vector3(0,0,-12.528367)

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
      [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
      [Law2_ScGeom_ViscElPhys_Basic()]
   ),
        NewtonIntegrator(gravity=(0,0,grav),damping=0.0),
        PyRunner(command='dissip(sumado)',iterPeriod=1,label='potDisip'),#Cuando comienzo a calcular las interacciones
        PyRunner(command='addPlotData()',iterPeriod=1),
]

from yade import qt

qt.View()


# Function to add data to plot
def addPlotData():
  posZ = O.bodies[EsfMo].state.pos[2]
  velZ = O.bodies[EsfMo].state.vel[2]
  PDcalc = 2.0 - posZ
  PDgeom = 0.0
  
  try:
    PDgeom = O.interactions[0,1].geom.penetrationDepth
  except:
    pass
  
  Fn = 0.0
  try:
    Fn = O.interactions[0,1].phys.normalForce[2]
  except:
    pass
    
  FnCalculated = 5e8*PDgeom - 1.5e6*velZ
  
  plot.addData(it1=O.iter, it2=O.iter, it3=O.iter, PDcalc = PDcalc,  PDgeom = PDgeom, velZ = velZ, Fn = Fn, FnYade =  -Fn, FnCalculated = FnCalculated)
  

plot.plots={'it1':('PDcalc','PDgeom'), 'it2':('velZ'), 'it3':('Fn'), 'PDgeom':('FnCalculated', 'FnYade')}; plot.plot()

def dissip(sumado):
        g = open (Arch1,'a')
        v=O.bodies[0].state.vel[2]
        z=O.bodies[0].state.pos[2]
        m=O.bodies[0].state.mass
        Fn=0.0
        for i in O.interactions:
                PD=i.geom.penetrationDepth
                YFn=i.phys.normalForce
                Fn=YFn[2]
                g.write(str(O.iter)+'\t'+str(z)+'\t'+str(v)+'\t'+'\t'+str(Fn)+'\t'+str(PD)+'\n')
        g.closed
O.dt=deltaT

O.run(200, True)
O.wait
plot.saveGnuplot('sim-data_Sphere')