yade-users team mailing list archive

[Anton Gladky <question253045@xxxxxxxxxxxxxxxxxxxxx>]

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

Anton Gladky posted a new comment:
Hi Dominik,

thanks for the script! I have pushed it to the trunk [1].

[1]
https://github.com/yade/trunk/commit/9030dc1af1a4a85d69a2a9bd7efc6680fcb67a9d

Anton


2014-08-26 20:42 GMT+02:00 Dominik Boemer
<question253045@xxxxxxxxxxxxxxxxxxxxx>:
> Question #253045 on Yade changed:
> https://answers.launchpad.net/yade/+question/253045
>
> Dominik Boemer posted a new comment:
> Hi Anton,
>
> thank you!
>
> Here is what you asked for.  I suppose that resultStatus is set to 0
> before calling the script.  If everything is fine, resultStatus returns
> 0, otherwise 1.  Do not hesitate to change this rule, if I misunderstood
> its meaning.
>
> Regards,
> Dominik
>
>
> #!/usr/bin/env python
> # encoding: utf-8
>
> ################################################################################
> # CONSTITUTIVE LAW TESTING: Law2_ScGeom_ViscElPhys_Basic()
> #
> # Two spheres with velocities (v,0,0) and (-v,0,0) collide.
> # This script checks if:
> #  - the numerical displacement equals the analytical displacement at a certain
> #    time before the spheres separate, for instance in t = 0.05 s
> #    This also implies the consistency of the results in terms of velocity
> #    (and acceleration).
> #  - the normal stiffness and normal damping coefficients have been calculated
> #    correctly in function of the normal coefficient (en) of restitution and the
> #    impact time (tc); this is a consequence of the previous condition.
> #
> # Notice that:
> #  - this script only checks the displacement before the separation because the
> #    analytical solution (given below) supposes that the damping term is still
> #    present, when the spheres separate.  This is however not true in the
> #    numerical model (see source code, prevent appearing of attraction forces
> #    due to a viscous component).
> #  - this script does not check the tangential (or shear) behaviour of
> #    the constitutive law.
> #
>
> from yade import plot,ymport
> import math
>
>
> ################################################################################
> # MATERIAL
>
> rho = 2000 # [kg/m^3] density
> mu  = 0.75 # [-]      friction coefficient
> tc  = 0.2  # [s]      contact time
> en  = 0.3  # [-]      normal restitution coefficient
> et  = 0.2  # [-]      tangential restitution coefficient
>
> frictionAngle = math.atan(mu)
> mat = O.materials.append(ViscElMat(tc=tc,en=en,et=et,
>                                    frictionAngle=frictionAngle,density=rho))
>
>
> ################################################################################
> # GEOMETRY
>
> r   = 0.1  # [m] sphere radius
>
> 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    = 2    # [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)
> ]
>
>
> ################################################################################
> # RUN
>
> O.run(5001,True)
>
>
> ################################################################################
> # COMPARE ANALYTICAL AND NUMERICAL SOLUTIONS
>
> m = 4./3. * math.pi * r**3 * rho # [kg] mass of the sphere
>
>
> # Normal stiffness and damping coefficients according to [Pournin2001]
> meff = m/2
> kn   = 2.0 * meff/tc**2 * (math.pi**2 + math.log(en)**2)
> cn   = -4.0 * meff/tc * math.log(en)
>
>
> # Analytical solution of a linear spring damper system
> omega0      = math.sqrt(kn/m)
> zeta        = cn / (2 * math.sqrt(kn * m))
> omegad      = omega0 * math.sqrt(1 - zeta**2)
> xAnalytical = v/omegad * math.exp(-zeta*omega0*O.time) * math.sin(omegad*O.time)
>
>
> # Comparison (if ok, resultStatus is not incremented)
> tolerance = 0.0001
>
> xNumerical = O.bodies[b2].state.pos[0]
>
> if ((abs(xNumerical-xAnalytical)/xAnalytical)>tolerance):
>   resultStatus += 1
>
> --
> You received this question notification because you are a member of
> yade-users, which is an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to     : yade-users@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users
> More help   : https://help.launchpad.net/ListHelp

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