yade-users team mailing list archive

[Dominik Boemer <question253045@xxxxxxxxxxxxxxxxxxxxx>]

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

Dominik Boemer posted a new comment:
Hello Anton,

sorry for contacting you again.  But I think that we missed something.
I have been reading [Pournin2001] and the source code again.  Meanwhile
I have also been testing the code with the script you can find below.
What I found, is that the equivalent mass

const Real massR = 2*mass1*mass2/(mass1+mass2);  // WITH A FACTOR 2!

delivers correct results.  Why?  Pournin's formulas are correct (except
the last expression in equation (22), as I mentionned earlier).  This
means that there has to be another bug in the code.  And this bug, is
the following:

[This is currently in the source code.]
const Real massR = mass1*mass2/(mass1+mass2);
kn1 = kn2 = 1/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR;
cn1 = cn2 = -2.0 /Tc * log(En)*massR;
ks1 = ks2 = 2.0/7.0 /Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(Et),2) )*massR;
cs1 = cs2 = -4.0/7.0 /Tc * log(Et)*massR;

should be replaced by:

[This should be in the source code.]
const Real massR = mass1*mass2/(mass1+mass2);
kn1 = kn2 = 2.0/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR;
cn1 = cn2 = -4.0 /Tc * log(En)*massR;
ks1 = ks2 = 4.0/7.0 /Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(Et),2) )*massR;
cs1 = cs2 = -8.0/7.0 /Tc * log(Et)*massR;

Why did I multiply kn1, kn2, cn1, cn2, ks1, ks2, cs1 and cs2 by 2
instead of multiplying massR by 2 (as it was before)?  Because massR is
given by the expression in [2001Pournin, just below equation (19)].
kn1, kn2, cn1, ... is however not given explicitly by equation
(22)[2001Pournin].  This equation gives the resultant stiffness and
damping coefficients kn, ks, cn and cs (and not the individual
components kn1, kn2, ...).

For instance,

kn = 1.0/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR;

according to [Pournin2001], and as kn1 = kn2,

1/kn = 1/kn1 + 1/kn2 <==> kn = kn1/2 <==> kn1 = 2*kn, which explains why
I multiplied the expressions of kn1, kn2, ... by 2.


To put it simple, we used the expressions of kn, ks, cn, ... in [Pournin2001] for kn1, kn2, ks1, ... which is wrong, as kn1 = 2*kn, kn2 = 2*kn, ks1 = 2*ks (when kn1 = kn2, ks1 = ks2, ...).

As mentionned above, here is the script, which helped me to validate the
constitutive model in the normal direction.  It shows that the contact
time is tc = 0.1 s and that the analytical velocity is -0.3 m/s when tc
= 0.1 s (this is the consequence of the normal restitution coefficient
en = 0.3).

Best regards,
Dominik


################################################################################
# CONSTITUTIVE LAW TESTING

from yade import plot,ymport
import math


################################################################################
# MATERIAL

rho = 1e-1 # [kg/m^3] density
r   = 0.5  # [m]      sphere radius
ks  = 0    # [N/m]    tangential stiffness coefficient
cs  = 0    # [-]      tangential damping coefficient
mu  = 0.75 # [-]      friction coefficient
tc  = 0.1  # [s]      contact time
en  = 0.3  # [-]      normal restitution coefficient

m    = 4./3. * math.pi * r**3 * rho

# m1*m2/(m1+m2)
meff = m/2

# [N/m] normal stiffness coefficient
kn = 2.0 * meff/tc**2 * (math.pi**2 + math.log(en)**2)
print kn

# [-]   normal damping coefficient
cn = -4.0 * meff/tc * math.log(en)
print cn

frictionAngle = math.atan(mu)
mat = O.materials.append(ViscElMat(kn=kn,ks=ks,cn=cn,cs=cs,
                                   frictionAngle=frictionAngle,density=rho))


################################################################################
# GEOMETRY

b1 = O.bodies.append(utils.sphere(center=(2*r,0,0),radius=r,material=mat))
b2 = O.bodies.append(utils.sphere(center=(0,0,0),radius=r,material=mat))


################################################################################
# SIMULATION
O.dt = 1e-5 # [s]   fixed time step
v    = 1    # [m/s] velocity along x-axis

O.bodies[b1].state.vel[0] = -v
O.bodies[b2].state.vel[0] =  v

O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(),
     Ig2_Facet_Sphere_ScGeom(),
     Ig2_Wall_Sphere_ScGeom()],
    [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
    [Law2_ScGeom_ViscElPhys_Basic()]
  ),

  NewtonIntegrator(damping=0),

  PyRunner(command='addPlotData()',iterPeriod=1)
]


################################################################################
# INTERMEDIATE PARAMETERS OF THE ANALYTICAL MODEL

omega0 = math.sqrt(kn/m)
zeta   = cn / (2 * math.sqrt(kn * m))
omegad = omega0 * math.sqrt(1 - zeta**2)


################################################################################
# PLOT NUMERICAL AND ANALYTICAL SOLUTIONS

def addPlotData():
  global v
  global omega0
  global zeta
  global omegad
  plot.addData(
    t  = O.time,
    xNumerical  = O.bodies[b2].state.pos[0],
    xAnalytical = v/omegad * math.exp(-zeta*omega0*O.time) 
                   * math.sin(omegad*O.time),
    vNumerical  = O.bodies[b2].state.vel[0],
    vAnalytical = v/omegad*(-zeta*omega0*math.exp(-zeta*omega0*O.time)
                   * math.sin(omegad*O.time) 
                   + omegad*math.exp(-zeta*omega0*O.time)
                   * math.cos(omegad*O.time)))

plot.plots={'t':('xNumerical','xAnalytical',None,'vNumerical','vAnalytical')}

plot.plot()


################################################################################
# RUN

O.run(20000)

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