Question #699273 on Yade changed: https://answers.launchpad.net/yade/+question/699273
Description changed to: Hello all, I have a simulation consisting of two spheres (S1, S2) and a clump made of two spheres (C1, C2). The two spheres, S1 and S2, are given the same initial speed but in the opposite directions along the horizontal. The clump is aligned along the vertical and the two spheres, S1 and S2, are made to hit the spheres in the clump at the same time. Could this be a bug in my code or the NewtonIntegrator? I expect the angular momentum to be conserved in the interactions, but it is not. I have a similar problem with clumps of pfacet [1]. It was reported in a bug report related to pfacet [2]. I am posting again because this appears to be an error related to clumps in general. I have put the code below. I am using Yadedaily. Kind regards, Rohit K. John [1] https://answers.launchpad.net/yade/+question/697465 [2] https://gitlab.com/yade-dev/trunk/-/issues/210 # script # ---------------------------------------------------------------------------- from yade.gridpfacet import * from yade import * from yade import geom, utils from yade import plot import sys, os # ---------------------------------------------------------------------------- input parameter # ----------------------------------------------------- clump sphere_young = 50e9 sphere_poisson = 0.3 sphere_friction = radians(30) sphere_offset = 0.3 sphere_radius = 5e-2 sphere_init_vel = 1 sphere_mass = 0.5 sphere_density = sphere_mass / (4/3*pi*sphere_radius**3) # ---------------------------------------------------------------------------- Materials sphere_mat = 'sphere_mat' O.materials.append( FrictMat( young = sphere_young, poisson = sphere_poisson, density = sphere_density, label = sphere_mat, frictionAngle = sphere_friction, ) ) # ---------------------------------------------------------------------------- Engines O.engines = [ ForceResetter(), InsertionSortCollider([ Bo1_Sphere_Aabb(), ]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()], ) ] # ---------------------------------------------------------------------------- objects # ----------------------------------------------------- clump clump_sph_1 = sphere([0, sphere_offset, 0], sphere_radius, material = sphere_mat) clump_sph_2 = sphere([0,-sphere_offset, 0], sphere_radius, material = sphere_mat) clump_sph_1_ID = O.bodies.append(clump_sph_1) clump_sph_2_ID = O.bodies.append(clump_sph_2) clump_clump_ID = O.bodies.clump([clump_sph_1_ID, clump_sph_2_ID]) clump_mmoi = O.bodies[clump_clump_ID].state.inertia # ----------------------------------------------------- spheres sp1 = sphere([ 4*sphere_radius, sphere_offset, 0], sphere_radius, material = sphere_mat) sp2 = sphere([-4*sphere_radius,-sphere_offset, 0], sphere_radius, material = sphere_mat) sp1_ID = O.bodies.append(sp1) sp2_ID = O.bodies.append(sp2) O.bodies[sp1_ID].state.vel = [-sphere_init_vel,0,0] O.bodies[sp2_ID].state.vel = [ sphere_init_vel,0,0] # ---------------------------------------------------------------------------- Additional engines O.engines += [ NewtonIntegrator( gravity = [0,0,0], damping = 0 ), PyRunner( command = 'graph()', iterPeriod = 1000 ), ] # ----------------------------------------------------- Calculating state after collision m = sphere_mass v0 = sphere_init_vel l = sphere_offset iz = clump_mmoi[2] calculated_final_sphere_lin_vel = 0.5*(-4.0*iz*l*m*v0/(iz + 2.0*l**2*m) + 2.0*l*m*v0)/(l*m) calculated_final_clump_ang_vel = 4.0*l*m*v0/(iz + 2.0*l**2*m) calculated_final_clump_ang_mom = clump_mmoi[2]*calculated_final_clump_ang_vel calculated_final_sphere_ang_mom = 2*sphere_offset*sphere_mass*calculated_final_sphere_lin_vel print("clump: ", calculated_final_clump_ang_mom) print("sphere: ", calculated_final_sphere_ang_mom) print("total: ", calculated_final_clump_ang_mom + calculated_final_sphere_ang_mom) # ----------------------------------------------------- plotting plot.plots = {'t':('clump_Lz', 'calculated_clump_final_Lz'), 't1':('Lz', 'calculated_Lz')} def graph(): L_sp1 = get_AngMom(sp1_ID) L_sp2 = get_AngMom(sp2_ID) angMom_clump = O.bodies[clump_clump_ID].state.angMom L_tot = angMom_clump + L_sp1 + L_sp2 plot.addData( t = O.time, t1 = O.time, Lx = L_tot[0], Ly = L_tot[1], Lz = L_tot[2], clump_Lz = angMom_clump[2], calculated = calculated_final_clump_ang_mom + calculated_final_sphere_ang_mom, calculated_clump_final_Lz = calculated_final_clump_ang_mom ) def get_AngMom(id): pos_sp1 = O.bodies[id].state.pos vel_sp1 = O.bodies[id].state.vel angMom_sp1 = O.bodies[id].state.angMom mass_sp1 = O.bodies[id].state.mass return angMom_sp1 + mass_sp1 * pos_sp1.cross(vel_sp1) # ----------------------------------------------------------------------------- simulation controls O.dt = 2e-7#utils.PWaveTimeStep() O.saveTmp() plot.plot() -- 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
