yade-dev team mailing list archive

[Dominik Boemer <dominik.boemer@xxxxxxxxx>]

Public bug reported:

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)

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

** Affects: yade
     Importance: Undecided
         Status: New

** Attachment added: "sphere_wall_analytic.py"
   https://bugs.launchpad.net/bugs/1360371/+attachment/4184658/+files/sphere_wall_analytic.py

** Description changed:

  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 solution.  For some reason,
- the inertial force is not influence by a change in density (rho).
+ 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))
+                                   frictionAngle=frictionAngle))
  m2 = O.materials.append(ViscElMat(kn=kn,ks=ks,cn=cn,cs=cs,
-                                   frictionAngle=frictionAngle,density=rho))
- 
+                                   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)
+   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),
+   NewtonIntegrator(damping=0),
  
-   PyRunner(command='addPlotData()',iterPeriod=1)
+   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))
- 
+   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)
  
  ============================================================================================

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