Question #706332 on Yade changed: https://answers.launchpad.net/yade/+question/706332
Description changed to: import random import math from yade import geom, pack, utils import pandas as pd # Define cylinder parameters center = (0, 0, 0) cyl_radius = 0.102 cyl_height = 0.064 # create cylindrical body with radius 0.102 m and height 0.064 m cylinder = geom.facetCylinder(center=center, radius=cyl_radius, height=cyl_height, segmentsNumber=80, wallMask=6) # define material properties mat = utils.PolyhedraMat() mat.density = 2600 # kg/m^3 mat.young = 1E6 # Pa mat.poisson = 20000 / 1E6 mat.frictionAngle = 0.6 # rad # assign material to each body in the cylinder for body in cylinder: body.bodyMat = mat # add cylinder to simulation O.bodies.append(cylinder) # Define box dimensions box_dimensions = (0.1, 0.1, 0.1) # Create box at origin box = geom.facetBox(center=(0,0,0), extents=box_dimensions, wallMask=31) # Define material properties for box mat = O.materials.append(FrictMat(density=2600, young=1E6, poisson=0.3, frictionAngle=0.6)) # Assign material to each facet of box for body in box: body.bodyMat = mat # Add box to simulation O.bodies.append(box) # define sphere parameters and number of spheres radius = 0.01575 num = 6 radius1 = 0.01175 num1 = 7 radius2 = 0.011 num2 = 8 radius3 = 0.01025 num3 = 8 radius4 = 0.0083125 num4 = 31 radius5 = 0.007125 num5 = 31 radius6 = 0.0059375 num6 = 32 radius7 = 0.0041525 num7 = 25 radius8 = 0.003555 num8 = 26 radius9 = 0.0029575 num9 = 26 radius10 = 0.002065 num10 = 25 radius11 = 0.00177 num11 = 26 radius12 = 0.001475 num12 = 26 radius13 = 0.001035 num13 = 155 radius14 = 0.00089 num14 = 155 radius15 = 0.000745 num15 = 156 radius16 = 0.00045 num16 = 9735 ## create empty sphere packing sp = pack.SpherePack() sp1 = pack.SpherePack() sp2 = pack.SpherePack() sp3 = pack.SpherePack() sp4 = pack.SpherePack() sp5 = pack.SpherePack() sp6 = pack.SpherePack() sp7 = pack.SpherePack() sp8 = pack.SpherePack() sp9 = pack.SpherePack() sp10 = pack.SpherePack() sp11 = pack.SpherePack() sp12 = pack.SpherePack() sp13 = pack.SpherePack() sp14 = pack.SpherePack() sp15 = pack.SpherePack() sp16 = pack.SpherePack() # generate randomly sphere sp.makeCloud((-0.055,-0.055,0.0105), (0.055,0.055,0.0105), rMean=radius, rRelFuzz=0, num=num) sp1.makeCloud((-0.055,-0.055,0.0380), (0.055,0.055,0.0380), rMean=radius1, rRelFuzz=0, num=num1) sp2.makeCloud((-0.055,-0.055,0.0608), (0.055,0.055,0.0608), rMean=radius2, rRelFuzz=0, num=num2) sp3.makeCloud((-0.055,-0.055,0.0820), (0.055,0.055,0.0820), rMean=radius3, rRelFuzz=0, num=num3) sp4.makeCloud((-0.0605,-0.0605,0.1006), (0.0605,0.0605,0.1006), rMean=radius4, rRelFuzz=0, num=num4) sp5.makeCloud((-0.0615,-0.0615,0.1160), (0.0615,0.0615,0.1160), rMean=radius5, rRelFuzz=0, num=num5) sp6.makeCloud((-0.0605,-0.0605,0.1291), (0.0605,0.0605,0.1291), rMean=radius6, rRelFuzz=0, num=num6) sp7.makeCloud((-0.0605,-0.0605,0.1392), (0.0605,0.0605,0.1392), rMean=radius7, rRelFuzz=0, num=num7) sp8.makeCloud((-0.0605,-0.0605,0.1469), (0.0605,0.0605,0.1469), rMean=radius8, rRelFuzz=0, num=num8) sp9.makeCloud((-0.0605,-0.0605,0.1534), (0.0605,0.0605,0.1534), rMean=radius9, rRelFuzz=0, num=num9) sp10.makeCloud((-0.0605,-0.0605,0.1584), (0.0605,0.0605,0.1584), rMean=radius10, rRelFuzz=0, num=num10) sp11.makeCloud((-0.0605,-0.0605,0.1622), (0.0605,0.0605,0.1622), rMean=radius11, rRelFuzz=0, num=num11) sp12.makeCloud((-0.0605,-0.0605,0.1655), (0.0605,0.0605,0.1655), rMean=radius12, rRelFuzz=0, num=num12) sp13.makeCloud((-0.0605,-0.0605,0.1680), (0.0605,0.0605,0.1680), rMean=radius13, rRelFuzz=0, num=num13) sp14.makeCloud((-0.0605,-0.0605,0.1699), (0.0605,0.0605,0.1699), rMean=radius14, rRelFuzz=0, num=num14) sp15.makeCloud((-0.0605,-0.0605,0.1715), (0.0605,0.0605,0.1715), rMean=radius15, rRelFuzz=0, num=num15) sp16.makeCloud((-0.0605,-0.0605,0.1727), (0.0605,0.0605,0.1727), rMean=radius16, rRelFuzz=0, num=num16) # add the sphere pack to the simulation sp.toSimulation() sp1.toSimulation() sp2.toSimulation() sp3.toSimulation() sp4.toSimulation() sp5.toSimulation() sp6.toSimulation() sp7.toSimulation() sp8.toSimulation() sp9.toSimulation() sp10.toSimulation() sp11.toSimulation() sp12.toSimulation() sp13.toSimulation() sp14.toSimulation() sp15.toSimulation() sp16.toSimulation() # add interaction of gravity of sphere and cylinder O.engines = [ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]), InteractionLoop( # handle sphere+sphere and facet+sphere collisions [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()] ), NewtonIntegrator(gravity=(0,0,-9.81), damping=0.4), ] O.dt = 0.1 * PWaveTimeStep() # run the simulation for 1000 steps O.run(1000, wait=True) for body in O.bodies: if not isinstance(body.shape, Sphere): continue if body.shape.radius == radius: #SP body.shape.color = (0,0,1) #blue if body.shape.radius == radius1: #SP1 body.shape.color = (1,0,0) #red if body.shape.radius == radius2: #SP2 body.shape.color = (0,1,0) #green if body.shape.radius == radius3: #SP3 body.shape.color = (1,1,1) #white if body.shape.radius == radius4: #SP4 body.shape.color = (1,1,0) #yellow if body.shape.radius == radius5: #SP5 body.shape.color = (1,0.5,0) #orange if body.shape.radius == radius6: #SP6 body.shape.color = (0.5,0,1) #purple if body.shape.radius == radius7: #SP7 body.shape.color = (0.5,0.5,0.5) #gray if body.shape.radius == radius8: #SP8 body.shape.color = (1,0.5,0.5) #pink if body.shape.radius == radius9: #SP9 body.shape.color = (0,0.85,0.8) #turquoise if body.shape.radius == radius10: #SP10 body.shape.color = (0.6,0.4,0.2) #brown if body.shape.radius == radius11: #SP11 body.shape.color = (0.5,0,0) #maroon if body.shape.radius == radius12: #SP12 body.shape.color = (0.5,0.5,0) #olive if body.shape.radius == radius13: #SP13 body.shape.color = (1,0,1) #magenta if body.shape.radius == radius14: #SP14 body.shape.color = (0,1,1) #cyan if body.shape.radius == radius15: #SP15 body.shape.color = (1,0.84,0) #gold if body.shape.radius == radius16: #SP16 body.shape.color = (0,0,0) #black # calculate the volume of the packing volume_packing = 0 num_spheres = 0 for b in O.bodies: if isinstance(b.shape, yade.wrapper.Sphere): volume_packing += 4/3 * math.pi * b.shape.radius**3 num_spheres += 1 # calculate the volume of the cylinder volume_cylinder = math.pi * cyl_radius**2 * cyl_height # calculate the porosity and porosity percentage porosity = (volume_cylinder - volume_packing) / volume_cylinder porosity_percent = porosity * 100 print("Number of spheres:", "{:d}".format(num_spheres)) print("V Packing:", "{:e}".format(volume_packing)) print("V Cylinder:", "{:e}".format(volume_cylinder)) print("Porosity:", "{:.2f}".format(porosity)) print("Porosity:", "{:.2f}%".format(porosity_percent)) # print the number of spheres in each layer print("Number of spheres in each layer:") print(f"Layer 1: {len(sp)}") print(f"Layer 2: {len(sp1)}") print(f"Layer 3: {len(sp2)}") print(f"Layer 4: {len(sp3)}") print(f"Layer 5: {len(sp4)}") print(f"Layer 6: {len(sp5)}") print(f"Layer 7: {len(sp6)}") print(f"Layer 8: {len(sp7)}") print(f"Layer 9: {len(sp8)}") print(f"Layer 10: {len(sp9)}") print(f"Layer 11: {len(sp10)}") print(f"Layer 12: {len(sp11)}") print(f"Layer 13: {len(sp12)}") print(f"Layer 14: {len(sp13)}") print(f"Layer 15: {len(sp14)}") print(f"Layer 16: {len(sp15)}") print(f"Layer 17: {len(sp16)}") # create a dictionary to store simulation data data = {'Number of spheres': [num_spheres], 'V Packing': [volume_packing], 'V Cylinder': [volume_cylinder], 'Porosity': [porosity], 'Porosity %': [porosity_percent]} # create a DataFrame from the data df = pd.DataFrame(data) # append the data to an existing file or create a new file if it doesn't exist with open('simulation_data.csv', 'a') as file: df.to_csv(file, header=file.tell() == 0, index=False) 1. Is it possible to only allows the particles to fall into the container without bouncing out after contact whether the contact might be against other spheres or the container? In the simulation, I tried to generate a box and try to contain all the sphere into the cylinder, but it didn't work as I expected it to -- 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