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 <[email protected]>: > 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 : [email protected] > 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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
