New question #707088 on Yade: https://answers.launchpad.net/yade/+question/707088
Hi, I'm trying to create a simulation that would have large spherical ball drop and hit the clump plate repeatedly 100 times. and calculate the average z and mass every time. Also, use export.text to get position of the packing every 10th time. However, the calculation would calculated every time the large ball drop instead of calculating in after it hit the plate. ---------------------------------------------------------------------------------This is my code------------------------------------------------------------------------------------------------ import random import math from yade import geom, pack, utils, plot, ymport, export import numpy as np # Define cylinder parameters diameter = 0.102 height = 0.18 center = (0, 0, height/2) # create cylindrical body with radius 0.102 m and height 0.064 m cylinder = geom.facetCylinder(center=center, radius=diameter/2, height=height, segmentsNumber=80, wallMask=6) # add cylinder to simulation O.bodies.append(cylinder) # add sphere packing O.bodies.append(ymport.textExt('packing8.txt',format='x_y_z_r')) # materials Properties gravel = CohFrictMat(young = 1e7, poisson = 0.25, density = 2700, label = 'gravel') asphalt_binder = CohFrictMat(young = 1e7, poisson = 0.25, density = 1060, frictionAngle = radians(40), normalCohesion = 5e4, shearCohesion = 5e4, label = 'asphalt_binder') weight = CohFrictMat(young = 1e7, poisson = 0.25, density = 11450,frictionAngle = radians(0), label = 'weight') # add properties O.materials.append(gravel) O.materials.append(asphalt_binder) O.materials.append(weight) # give color and properties to shpere for body in O.bodies: if not isinstance(body.shape, Sphere): continue if body.shape.radius == 0.01575/2 : body.shape.color = (0,0,1) #blue body.material = gravel if body.shape.radius == 0.011/2: body.shape.color = (1,0,0) #red body.material = gravel if body.shape.radius == 0.007125/2: body.shape.color = (0,1,0) #green body.material = gravel if body.shape.radius == 0.003555/2: body.shape.color = (1,1,0) #yellow body.material = gravel if body.shape.radius == 0.00160/2 : body.shape.color = (1,0,1) #magenta body.material = gravel if body.shape.radius == 0.0008/2 : body.shape.color = (0,0,0) #black body.material = asphalt_binder # add clump plate clump_bodies = ymport.textExt('clump8.txt',format='x_y_z_r') # plate properties clump_plate = CohFrictMat(density = 7500, label = 'clump_plate') # add properties O.materials.append(clump_plate) # define layer total_clump_bodies = len(clump_bodies) bodies_per_layer = total_clump_bodies / 8 for i, clump_body in enumerate(clump_bodies): layer_number = i // bodies_per_layer # Calculate the layer number for the clump body if layer_number % 2 == 0: color = (1, 0, 0) # red for even-numbered layers else: color = (1, 1, 1) # white for odd-numbered layers clump_body.shape.color = color clump_body.material = clump_plate O.bodies.appendClumped(clump_bodies) # z-coordinate for clump clump_z = np.mean([clump_body.state.pos[2] for clump_body in clump_bodies]) # clump (center of mass) z-coordinate # create large ball O.bodies.append(sphere((0, 0, clump_z + 0.25), 0.1/2)) # give color and properties to shpere for body in O.bodies: if not isinstance(body.shape, Sphere): continue if body.shape.radius == 0.10/2 : body.shape.color = (0.5,0.5,0.5) #grey body.material = weight # define original condition x = 0 y = 0 window = 0.01575/2 def calculate_zmax(x, y, window): zmax = float('-inf') # Initialize zmax to negative infinity # Define the square region x_min = x - window x_max = x + window y_min = y - window y_max = y + window # Iterate over all bodies in the simulation for body in O.bodies: if isinstance(body.shape, Sphere) and (body.material == gravel or body.material == asphalt_binder or body.material == weight): sphere_x, sphere_y, sphere_z = body.state.pos # Get the position of the sphere sphere_radius = body.shape.radius # Check if the sphere is within the square region if x_min <= sphere_x <= x_max and y_min <= sphere_y <= y_max: z = sphere_z + sphere_radius # Calculate the z-coordinate of the top of the sphere # Update zmax if the current z-coordinate is higher if z > zmax: zmax = z return zmax def hundredblows(): # Perform 100 blows of the large ball for blow_number in range(100): # Perform the blow O.run() # Adjust the number of iterations as needed # Calculate clump_z clump_z = np.mean([clump_body.state.pos[2] for clump_body in clump_bodies]) # clump (center of mass) z-coordinate # Remove spheres above clump_z to_remove = [] # List to store the bodies to be removed for body in O.bodies: if isinstance(body.shape, Sphere) and (body.material == gravel or body.material == asphalt_binder): if body.state.pos[2] > clump_z: to_remove.append(body) for body in to_remove: O.bodies.erase(body.id) # Calculate and print the average z-coordinate random_points = [(0, 0), (0.025, 0), (0, 0.025), (-0.025, 0), (0, -0.025)] zmax_values = [] # List to store the zmax values for point in random_points: x, y = point zmax = calculate_zmax(x, y, window) zmax_values.append(zmax) average_zmax = sum(zmax_values) / len(zmax_values) print("Average z-coordinate: {}".format(average_zmax)) # Calculate total mass excluding clump total_mass = 0.0 clump_sphere_ids = set() # Create a set of all sphere IDs in the clumps for clump in clump_bodies: if isinstance(clump.shape, Sphere): clump_sphere_ids.add(clump.id) for body in O.bodies: if isinstance(body.shape, Sphere) and body.id not in clump_sphere_ids and body.material != weight: volume = (4/3) * math.pi * (body.shape.radius**3) mass = body.material.density * volume total_mass += mass print("Total mass excluding Clump:", total_mass) # Export positions to a text file if (blow_number + 1) % 10 == 0: export.text("dense_{}.txt".format(blow_number + 1)) 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), # call the checkUnbalanced function (defined below) every 600 seconds PyRunner(command='checkUnbalanced()', realPeriod=30) ] O.dt = PWaveTimeStep() def checkUnbalanced(): if unbalancedForce() < .1: O.pause() # Call hundredblows() function hundredblows() # Run the simulation O.run() -- 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