yade-dev team mailing list archive

[Dominik Boemer <dominik.boemer@xxxxxxxxx>]

Hi Anton,

thank you!  In this case, there seems to be a problem in version 1.07.0
(the standard package in Ubuntu 14.04) of Yade.  This bug has thus been
fixed in yadedaily.  If you want to check if everything works correcly,
you could replace "1" by its comment "4/3 * math.pi * r**3 * rho", which
is the real mass of the sphere.  Moreover, the stiffness coefficient has
to be adapted such that the sphere is not going through the wall.  In
this case the ploted curves should overlap during the contact.

I will update my version to yadedaily.

Thank you,
Dominik

-- 
You received this bug notification because you are a member of Yade
developers, which is subscribed to Yade.
https://bugs.launchpad.net/bugs/1360371

Title:
  Law2_ScGeom_ViscElPhys_Basic - Mass Always Equals 1

Status in Yet Another Dynamic Engine:
  New

Bug description:
  Hello everyone,

  I have been testing the constitutive law Law2_ScGeom_ViscElPhys_Basic
  (or rather Law2_Spheres_Viscoelastic_SimpleViscoelastic() in Version
  1.07.0) by comparing the analytical [1] with the numerical solution.
  The test scenario is a sphere that collides with a wall.  I discovered
  that the model (either the constitutive law or the integrator) did not
  take into account the density of the sphere, even if its mass was
  correct; maybe it's just a syntax error.

  To illustrate the problem I created a minimal working example, in
  which the sphere possesses a density of 1e6 kg/m³.  Despite this high
  density, the sphere's displacement is still that of a 1 kg sphere.  In
  this example, I plot the analytical and numerical solutions.  For some
  reason, the inertial force is not influence by a change in density
  (rho).

  [1]
  http://en.wikipedia.org/wiki/Damping#Example:_mass.E2.80.93spring.E2.80.93damper

  See either the following code or the attached file for the MWE.

  Best regards and thank you in advance,
  Dominik

  ============================================================================================

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

  from yade import plot,ymport
  import math

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

  rho = 1e6   # [kg/m^3] density
  kn  = 600    # [N/m]       normal stiffness coefficient
  ks  = 0         # [N/m]       tangential stiffness coefficient
  cn  = 40      # [-]              normal damping coefficient
  cs  = 0         # [-]              tangential damping coefficient
  mu  = 0.75 # [-]              friction coefficient

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

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

  r = 1 # [m] sphere radius

  b1 = O.bodies.append(utils.wall((r,0,0),0,0,material=m1))
  b2 = O.bodies.append(utils.sphere(center=(0,0,0),radius=r,material=m2))

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

  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_Spheres_Viscoelastic_SimpleViscoelastic()]
    ),
    # new law name should be Law2_ScGeom_ViscElPhys_Basic
    # (but not updated in Version 1.07.0)

    NewtonIntegrator(damping=0),

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

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

  m      = 1 #4/3 * math.pi * r**3 * rho (this should be the real mass)
  omega0 = math.sqrt(kn/m)
  zeta   = cn / (2 * math.sqrt(kn * m))
  omegad = omega0 * math.sqrt(1 - zeta**2)

  ################################################################################
  # PLOT NUMERICAL AND ANALYTICAL SOLUTION

  def addPlotData():
    global v
    global omega0
    global zeta
    global omegad
    plot.addData(
      t  = O.time,
      x1 = O.bodies[b2].state.pos[0],
      x2 = v/omegad * math.exp(-zeta*omega0*O.time) * math.sin(omegad*O.time))

  plot.plots={'t':('x1','x2')}

  plot.plot()

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

  O.run(20000)

  ============================================================================================

To manage notifications about this bug go to:
https://bugs.launchpad.net/yade/+bug/1360371/+subscriptions