Yade

Dynamic Compaction Can't Repeat

Asked by Huan

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()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Huan
Solved:
Last query:
Last reply:
Revision history for this message
Karol Brzezinski (kbrzezinski) said : 		#1

Hi,

Honestly, I don't understand the question. Could you please rephrase it?

Secondly, I thought that running the code would be helpful, but I the packing is stored in the external file 'packing8.txt'.

Please try to simplify your example (MWE[1])

Cheers,
Karol

[1] https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Huan (huan-liu) said : 		#2

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()