yade-dev team mailing list archive

["Alexander Eulitz \[Eugen\]" <1360371@xxxxxxxxxxxxxxxxxx>]

Thanks for pointing that out and for presenting the script!

-- 
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