yade-users team mailing list archive

[Huan <question707088@xxxxxxxxxxxxxxxxxxxxx>]

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.