yade-dev team mailing list archive

[Dominik Boemer <dominik.boemer@xxxxxxxxx>]

Hi everybody,

I checked the script with yadedaily and everything is fine.  In
conclusion, the bug exists in Yade 1.07.0 but not in yadedaily.  Here is
the script, if anyone wants to test the constitutive law:

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

from yade import plot,ymport
import math

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

rho = 1e-1 # [kg/m^3] density
kn  = 1e3  # [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_ScGeom_ViscElPhys_Basic()]
  ),

  NewtonIntegrator(damping=0),

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


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

kn = kn/2
cn = cn/2

m      = 4./3. * math.pi * r**3 * rho
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),
    v1 = O.bodies[b2].state.vel[0])


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

plot.plot()


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

O.run(20000)

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