Question #707088 on Yade changed: https://answers.launchpad.net/yade/+question/707088
Status: Needs information => Solved Huan confirmed that the question is solved: Thank Karol, I resolve my issue. --------------------------------Resolve of 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('packing_loose.txt',format='x_y_z_r')) # materials Properties gravel = CohFrictMat(young = 1e6, poisson = 0.25, density = 2700, label = 'gravel') asphalt_binder = CohFrictMat(young = 3e5, poisson = 0.25, density = 1060, frictionAngle = radians(40), normalCohesion = 1e4, shearCohesion = 1e4, label = 'asphalt_binder') weight = CohFrictMat(young = 1.9e8, poisson = 0.28, density = 8645,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('clump_loose.txt',format='x_y_z_r') # plate properties clump_plate = CohFrictMat(young = 1.9e8, poisson = 0.28, 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 layer_7_bodies = [] 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 # Add clump body to layer 7 list if layer_number == 7: layer_7_bodies.append(clump_body) # Calculate the mean z-coordinate of layer 7 center_top_clump_surface = np.mean([clump_body.state.pos[2] for clump_body in layer_7_bodies]) O.bodies.appendClumped(clump_bodies) # clump center of mass clump_z = np.mean([clump_body.state.pos[2] for clump_body in clump_bodies]) # clump (center of mass) z-coordinate 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 ] # set dt to stabilize the force in packing O.dynDt = False # function for stabalize def impaction(): # stabalize the packing O.dt = PWaveTimeStep() print("dt=",O.dt) O.run(2000,True) print("2000") # create large ball idball = O.bodies.append(sphere((0, 0, center_top_clump_surface + 0.05721), 0.1/2)) O.bodies[idball].material = weight O.bodies[idball].shape.color = (0.5, 0.5, 0.5) s = 0.445 g = 9.81 v = math.sqrt(2*g*s) O.bodies[idball].state.vel = (0,0,-v) # stabilize when ball impact packing O.run(5000,True) print("7000") O.run(5000,True) print("12000") O.run(5000,True) print("17000") # remove ball after impact O.bodies.erase(idball) print("weight removed") # list to store the bodies to be removed to_remove = [] # remove asphalt and aggregate that is above clump 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) # check unbalanced force print("checking_unbalanced") i_unbalanced = unbalancedForce() print("i_unbalanced =",i_unbalanced) tocontinue = True while (tocontinue): O.run(1000, True) n_unbalanced = unbalancedForce() print("n_unbalanced=",n_unbalanced) if (n_unbalanced < 0.01) and ((abs(i_unbalanced-n_unbalanced)/i_unbalanced) < 0.01): tocontinue = False else: i_unbalanced = n_unbalanced # 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): 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 # Randomly select 5 points random_points = [(0, 0), (0.025, 0), (0, 0.025), (-0.025, 0), (0, -0.025)] zmax_values = [] # List to store the zmax values average_zmax = 0.0 def calculate_average_zmax(random_points, zmax_values): global average_zmax # Declare average_zmax as a global variable 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 mass of individual spheres excluding clumped spheres total_mass = 0.0 def calculate_mass(): global total_mass # Declare total_mass as a global variable mass = 0.0 # Iterate over the bodies and calculate mass for individual spheres for body in O.bodies: if isinstance(body.shape, Sphere) and (body.material == gravel or body.material == asphalt_binder): volume = (4/3) * math.pi * (body.shape.radius**3) packing_mass = body.material.density * volume mass += packing_mass total_mass = mass print("Total mass of packing: {}".format(total_mass)) return total_mass def impaction_loop(): with open("output_d04_co001.txt", "w") as file: for looping in range(1): impaction() calculate_zmax(x, y, window) calculate_average_zmax(random_points, zmax_values) calculate_mass() file.write("Loop {}: Average height= {}\n".format(looping + 1, average_zmax)) file.write("Loop {}: Total mass= {}\n".format(looping + 1, total_mass)) file.write("\n") # Export positions of the packing every tenth iteration if looping == 0 : export.text('packing_positions_{}.txt'.format(looping + 1)) print("Output saved to 'output_d04_co001.txt' file.") impaction_loop() -- 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