New question #690122 on Yade: https://answers.launchpad.net/yade/+question/690122
Good morning, I'm new to YADE and I'm using it for a part of my master thesis. I want to model a granular damper. In the case study the phases are the following 1) The spheres are inserted in a parallelepiped box (approx. 2cmx1cmx0.5cm) 2) The vibration of the all system makes the box shear cyclically. (code at the end) IMG to schematize the system ---> https://we.tl/t-0e9ZcKzmhp The idea is to build before all the objects (box as facets and sphere pack), then to compress the pack of spheres by imposing a motion on the cover and finally to excite the system with HarmonicRotationEngine. I use Hertz-Mindlin contact law (code below). I'm facing some "physical" problems: 1) I want to get the total energy given to the system, so I save the Torque aruond z and the rotation to integrate over an half cycle. Is there another way to get this total energy directly from YADE? Am I taking the wrong data or interpreting them bad? 2) From Law2_ScGeom_MindlinPhys_Mindlin I take the terms of dissipation but what I noticed is the loss rate (dissipated energy per cycle) is higher than the stored energy per cycle estimated before. I think there must be some problem with how I compute the stored energy or within the simulation. 3) I want to relate the confinement pressure (F_closing in code) with the dissipation in this system, so I expect that at low pressure friction is not effective and at very high pressure is not effective as well, with a sort of optimum in the middle. What I see is that I can increase the pressure on the pack more and more (up to complete compenetration of the sphere layers) and the dissipation by friction only increases to not physical values? What is the reason of this? Am I missing something in controlling? 4) Are there other ways to achieve the my result? If yes, are there some examples to referr to? Thank you very much for the helping and for the care you all have in this community. I really hope that you can help me with this problems that make me stuck. I hope as well that I did respect all the rules for making questions on this forum. Thank you again and regards. Stefano Simone ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- The code is the following: from yade import plot,utils from yade import pack import math import random import numpy #parameters f = 250 period = 1/f wrad = 2*pi*f ang_ampl = 0.04 e_n = 0.63 e_s = 0.63 frictionAngle = 0.197 num_per = 1 g = 9.81 F_closing= 15000 #materials #O.materials.append(FrictMat(young = 210e9, poisson = 0.35, frictionAngle = 0.7, density = 7800, label = 'spmat')) O.materials.append(FrictMat(young = 70e9, poisson = 0.22, frictionAngle = 0.197, density = 2500 , label = 'spmat')) #sphere of glass O.materials.append(FrictMat(young = 210e9, poisson = 0.35, frictionAngle = 0, density = 0, label = 'wmat')) #bodies # sphere packing creation pred = pack.inAlignedBox((-0.010,-0.005,0), (.010,0.005,0.005)) O.bodies.append(pack.regularHexa(pred, radius=0.00075, gap=0.00001, material = 'spmat')) # random distribution of radii Dev = 0.0001 media = 0.00075 val_max = media + Dev val_min = media - Dev radii = numpy.linspace(val_min, val_max, num=50, endpoint=True, retstep=False) for b in O.bodies: if isinstance(b.shape,Sphere): b.shape.radius = random.choice(radii) # total volume and number of the spheres V_spheres = 0 number = 0 for b in O.bodies: if isinstance(b.shape, Sphere): V_spheres = V_spheres + 4/3*pi*b.shape.radius**3 number = number + 1 # Container Down = O.bodies.append(geom.facetBox((0,0,0),(0.02,0.005,0))) Right = O.bodies.append(geom.facetBox((0.01,-0.005,0.0025), (0,0.005*2,0.0025))) Left = O.bodies.append(geom.facetBox((-0.01,-0.005,0.0025), (0,0.005*2,0.0025))) Bottom = O.bodies.append(geom.facetBox((0,-0.0050,0.005), (0.01,0,0.005))) Upper = O.bodies.append(geom.facetBox((0,0.005,0.005), (0.014,0,0.005))) Cover = O.bodies.append(geom.facetBox((0,0,0.005),(0.02,0.005,0))) #engines O.engines = [ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_MindlinPhys(frictAngle = 0.05, en=e_n , es=e_s, label='Ip2')], [Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=False, calcEnergy=True, label='Mindlin')] ), NewtonIntegrator(gravity=(0, 0, -g), damping=0.1), TranslationEngine(dead = True, translationAxis=[0, 0, 1], velocity=-0.002, ids=Cover, label='pres'), HarmonicRotationEngine(dead=True, f =f, A = ang_ampl*pi, fi = 0*pi/2 , ids=Right, rotateAroundZero=True, rotationAxis=(0, 0, 1), zeroPoint=(0.01, -0.005, 0), label = 'rotR'), HarmonicRotationEngine(dead=True, f =f, A = ang_ampl*pi, fi = 0*pi/2 , ids=Left, rotateAroundZero=True, rotationAxis=(0, 0, 1), zeroPoint=(-0.01, -0.0050, 0), label = 'rotL'), #HarmonicMotionEngine(dead=True, f =[f,0,0], fi = [-pi/2*0, 0, 0],A = [sin(ang_ampl*pi),0,0] , ids=Upper, label = 'rotU'), PyRunner(dead = False, command = 'startPressure()', iterPeriod = 100, label = 'phase'), PyRunner(dead = False, command = 'plotData()', iterPeriod = 1000, label = 'data') ] O.dt = PWaveTimeStep() O.trackEnergy = True def startPressure(): if O.iter > 250000 and unbalancedForce() < 0.05: pres.dead = False pres.velocity = -0.05 phase.command = 'stopPressure()' phase.iterPeriod = 100 def stopPressure(): top_sphere = max([b.state.pos[2] for b in O.bodies if isinstance(b.shape, Sphere)]) if (O.bodies[Cover[0]].state.pos[2] - top_sphere)< 0.0005: pres.velocity = -0.05 if abs((O.forces.f(O.bodies[Cover[0]].id)[2] + O.forces.f(O.bodies[Cover[1]].id)[2])/(1))*1 >= F_closing: #if O.bodies[Cover[0]].state.pos[2] <=0.5: Ip2.frictAngle = frictionAngle O.materials[0].frictionAngle = frictionAngle O.materials[1].frictionAngle = frictionAngle Ip2.en = 0.63 Ip2.es = 0.63 print(((O.forces.f(O.bodies[Cover[0]].id)[2] + O.forces.f(O.bodies[Cover[1]].id)[2])/(1))*1) pres.dead = True O.bodies[Cover[0]].state.blockedDOFs = 'xyzXYZ' O.bodies[Cover[0]].state.vel = (0, 0, 0) O.bodies[Cover[1]].state.blockedDOFs = 'xyzXYZ' O.bodies[Cover[1]].state.vel = (0,0,0) phase.command = 'Rotation_control()' def Rotation_control(): if unbalancedForce() < 0.02 and O.time > 0.4 : rotR.dead = False rotR.zeroPoint[0] = O.bodies[Right[0]].state.pos[0] rotL.dead = False #rotU.dead = False data.dead = False def plotData(): plot.addData( t = O.time, Edf = Mindlin.frictionDissipation, Eds = Mindlin.shearDampDissip, Edn = Mindlin.normDampDissip, Etot = Mindlin.normElastEnergy() + Mindlin.shearEnergy, Ed = Mindlin.frictionDissipation + Mindlin.shearDampDissip + Mindlin.normDampDissip, v_av = numpy.average([b.state.vel.norm() for b in O.bodies if isinstance(b.shape, Sphere)]), TorqueR = O.forces.t(O.bodies[Right[0]].id)[2] + O.forces.t(O.bodies[Right[1]].id)[2], TorqueL = O.forces.t(O.bodies[Left[0]].id)[2] + O.forces.t(O.bodies[Left[1]].id)[2], ForceUP = -O.forces.f(O.bodies[Cover[0]].id)[2] - O.forces.f(O.bodies[Cover[1]].id)[2], v_UP = O.bodies[Cover[0]].state.vel[2], Angle = O.bodies[Right[0]].state.ori[2], AngVel = O.bodies[Right[0]].state.angVel[2], Pressure = ((O.forces.f(O.bodies[Cover[0]].id)[2] + O.forces.f(O.bodies[Cover[1]].id)[2])/(0.0002))*1e-6, Unbalanced = unbalancedForce(), **O.energy ) # if O.iter > 1000000: # plot.saveDataTxt('Side.dat', vars=('t', 'Edf', 'Eds', 'Edn', 'TorqueL','TorqueR','ForceUP', 'v_UP','AngVel','Angle','Unbalanced','Pressure')) plot.plots = {'t':('Ed','kinetic','Etot'), 't ':'TorqueR', 't ':'Angle', 't ': 'TorqueL', 't ':'v_UP'} plot.plot() O.saveTmp() -- You received this question notification because your team yade-users 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
