← 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

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.


Follow ups