New question #697465 on Yade: https://answers.launchpad.net/yade/+question/697465
Hello, In the following simulation, two sphere are fired in opposite direction at a pfacet box clump. Ideally after collision, the box starts rotating and the spheres move in the opposite direction (opposite to their corresponding initial velocity) such that the angular momentum of the whole system is conserved. I want to change the inertia and the mass of the pfacet box clump, so I did that using the lines O.bodies[target_clump_ID].state.inertia = Vector3([10,10,10]) O.bodies[target_clump_ID].state.mass = 1 When I do this, the angular momentum of the system is not conserved. Am I not seeing something here, have I made a mistake? Moreover, in the comment #20 of [1], it was found that for pfacet clump to have correct dynamics, we should set the mass of the cylinder connections to 0 (as the mass is supposed to be concentrated at the nodes). But when I do that I am no longer able to access the angular momentum of the clump using O.bodies[target_clump_ID].state.angMom as it returns Vector3(0,0,0). How do I fix this? The only work around I found was to set the mass close to zero (1e-10) instead of 0. But is there a proper fix for this? Kind regards, Rohit K. John https://answers.launchpad.net/yade/+question/695558 # GTS file # -------------------------------------------------------------------------------------------------------------------------------------------------------- 14 36 24 GtsSurface GtsFace GtsEdge GtsVertex 0.5 0.5 0.5 0.5 0.5 -0.5 0.5 -0.5 0.5 0.5 -0.5 -0.5 -0.5 0.5 0.5 -0.5 0.5 -0.5 -0.5 -0.5 0.5 -0.5 -0.5 -0.5 0.5 0.0 0.0 0.0 -0.5 0.0 0.0 0.0 0.5 -0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 -0.5 6 8 2 6 1 2 8 7 3 4 5 6 3 7 1 3 8 4 7 5 5 1 4 2 9 3 2 9 4 9 1 9 10 7 4 10 8 10 3 10 11 3 5 11 7 11 1 11 12 6 7 12 8 12 5 12 13 2 5 13 6 13 1 13 14 8 2 14 6 14 4 14 7 21 23 7 17 20 6 25 28 9 33 36 8 13 16 3 29 32 12 14 15 5 15 13 3 16 14 9 18 19 4 19 17 5 20 18 11 22 24 10 23 22 8 24 21 4 26 27 1 27 25 10 28 26 6 30 31 2 31 29 11 32 30 2 34 35 1 35 33 12 36 34 # script # -------------------------------------------------------------------------------------------------------------------------------------------------------- from yade.gridpfacet import * from yade import * from yade import geom, utils from yade import plot import sys, os sys.path.append(".") # ---------------------------------------------------------------------------- input parameter # ----------------------------------------------------- target target_young = 50e9 target_density = 1000 target_poisson = 0.3 target_friction = radians(30) p_radius = 5e-2 # ---------------------------------------------------------------------------- Materials target_int_mat = 'pfacet_int_mat' target_ext_mat = 'pfacet_ext_mat' O.materials.append( FrictMat( young = target_young, poisson = target_poisson, density = target_density, label = target_ext_mat, frictionAngle = target_friction, ) ) O.materials.append( CohFrictMat( young = target_young, poisson = target_poisson, density = target_density, label = target_int_mat, frictionAngle = target_friction, normalCohesion = 3e100, shearCohesion = 3e100, momentRotationLaw = True, )) # ---------------------------------------------------------------------------- Engines O.engines = [ ForceResetter(), InsertionSortCollider([ Bo1_GridConnection_Aabb(), Bo1_PFacet_Aabb(), Bo1_Sphere_Aabb(), ]), InteractionLoop( [ Ig2_PFacet_PFacet_ScGeom(), Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), Ig2_GridNode_GridNode_GridNodeGeom6D(), Ig2_GridConnection_PFacet_ScGeom(), Ig2_Sphere_PFacet_ScGridCoGeom(), ], [ Ip2_FrictMat_FrictMat_FrictPhys(), Ip2_CohFrictMat_CohFrictMat_CohFrictPhys( setCohesionNow = True, setCohesionOnNewContacts = False ), ], [ Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGridCoGeom_FrictPhys_CundallStrack(), Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), ], ) ] # ---------------------------------------------------------------------------- objects # ----------------------------------------------------- target ( pnode, pcyl, pfacet ) = gtsPFacet( 'cube.gts', radius = p_radius, shift = (0,0,0), scale = 1, wire = False, fixed = False, color = [0.1,0.5,0.1], materialNodes = 'pfacet_int_mat', material = 'pfacet_ext_mat', ) target_ids = pnode + pcyl + pfacet # for i in pcyl: # O.bodies[i].state.mass = 0 # See https://answers.launchpad.net/yade/+question/695558 comment #20 for why I am doing this # If this is uncommented, then we are not able to access the angular momentum of the clump using # O.bodies[target_clump_ID].state.angMom target_clump_ID = O.bodies.clump(target_ids) O.bodies[target_clump_ID].state.inertia = Vector3([10,10,10]) O.bodies[target_clump_ID].state.mass = 1 print(target_clump_ID) # ----------------------------------------------------- spheres sp1 = sphere([ 0.61, 0.3, 0], 5e-2, material = target_ext_mat) sp2 = sphere([-0.61, -0.3, 0], 5e-2, material = target_ext_mat) sp1_ID = O.bodies.append(sp1) sp2_ID = O.bodies.append(sp2) O.bodies[sp1_ID].state.vel = [-1,0,0] O.bodies[sp2_ID].state.vel = [ 1,0,0] # ---------------------------------------------------------------------------- Additional engines ids = target_clump_ID O.engines += [ NewtonIntegrator(gravity = [0,0,0], damping = 0), PyRunner(command = 'graph()', iterPeriod = 1000) ] # ----------------------------------------------------- plotting plot.plots = {'t':('wx', 'wy', 'wz'), 't1':('Lz')} def graph(): w = O.bodies[ids].state.angVel pos_sp1 = O.bodies[sp1_ID].state.pos vel_sp1 = O.bodies[sp1_ID].state.vel angMom_sp1 = O.bodies[sp1_ID].state.angMom mass_sp1 = O.bodies[sp1_ID].state.mass L_sp1 = angMom_sp1 + mass_sp1 * pos_sp1.cross(vel_sp1) pos_sp2 = O.bodies[sp2_ID].state.pos vel_sp2 = O.bodies[sp2_ID].state.vel angMom_sp2 = O.bodies[sp2_ID].state.angMom mass_sp2 = O.bodies[sp2_ID].state.mass L_sp2 = angMom_sp2 + mass_sp2 * pos_sp2.cross(vel_sp2) angVel_clump = O.bodies[target_clump_ID].state.angVel mmoi_clump = O.bodies[target_clump_ID].state.inertia mmoi_mat_clump = Matrix3([ mmoi_clump[0], 0, 0, 0, mmoi_clump[1], 0, 0, 0, mmoi_clump[2] ]) L_clump = mmoi_mat_clump * angVel_clump L_tot = L_clump + L_sp1 + L_sp2 angMom_clump = O.bodies[target_clump_ID].state.angMom # L_tot = angMom_clump + L_sp1 + L_sp2 print(angMom_clump) plot.addData( t = O.time, t1 = O.time, wx = w[0], wy = w[1], wz = w[2], Lx = L_tot[0], Ly = L_tot[1], Lz = L_tot[2], ) plot.plot() # ---------------------------------------------------------------------------- GUI controller from yade import qt qt.Controller() qtv = qt.View() qtr = qt.Renderer() qtr.light2 = True qtr.lightPos = Vector3(1200,1500,500) qtr.bgColor = [1,1,1] qtv.ortho = True # orthographic view in viewport Gl1_Sphere.stripes = True O.dt = 2e-7#utils.PWaveTimeStep() 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 : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp