← Back to team overview

yade-users team mailing list archive

Re: [Question #248996]: Viscoelastic Force Model

 

Question #248996 on Yade changed:
https://answers.launchpad.net/yade/+question/248996

Anton Gladky proposed the following answer:
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

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.