UCS Test With Imported Geometry

In this example, a compression test is performed on an assembly of balls for which the initial packing geometry is imported from a file.

The procedure is as follows:

Read the packing geometry from a file and recreate the balls in PFC
Initialize the mechanical properties. In this example, the linear parallel bond contact model is used.
Apply a constant downward velocity to particles on the top of the specimen and fix the z coordinate of particles on the lower boundary of the specimen.
Record the stress and strain during cycling and plot the results.

Note

This example shows a simple method for conducting a UCS test on a bonded particle specimen. See http://www.itascacg.com/material-modeling-support for more information on material testing with PFC.

First, we import the required modules and give some commands to setup the model domain. We also set the default slots of the Contact Model Assignment Table. In this example the linear parallel bond is used for all contacts.

import itasca as it

from itasca import ballarray as ba

from itasca import ball as balls

from itasca import contact as contacts

import numpy as np

it.command("""

model new

model large-strain on

model domain extent -1 1 -1 1 -1 1

contact cmat default model linearpbond property pb_ten 12.0e6 pb_coh 10.0e6 pb_fa 20.0 fric 0.5  ...

                   method deformability    emod 8.5e9 kratio 2.5                         ...

                   method pb_deformability emod 8.5e9 kratio 2.5                         ...

                   proximity 4e-4

""")

Particles are created at the locations specified in the file: balls.txt. This is a comma separated value (.csv) file. We can read this as a numpy array with the np.loadtxt function. This file contains one line per ball, with the following format:

the first column stores the radius of the ball,

the three following columns store the x,y and z coordinates,

the last column tells us if this particle has a boundary condition applied.

box_len = 38.1e-3

data = np.loadtxt("balls.txt",delimiter=",")

for datum in data:

    r, p, fix = datum[0],datum[1:4], datum[4]

    b = balls.create(r, p)

    if fix == 1.0:

        if b.pos_z() > 0.0:

            b.set_extra(1, 1)

        else:

            b.set_extra(1, 2)

    else:

        b.set_extra(1,0)

next we set some ball attributes, clean the model to create the initial network of contacts, and bond the active contacts.

it.command("""

ball attribute density 2000.0 damp 0.6

model clean

contact method bond gap 4e-4

contact property lin_mode 1

contact cmat default proximity 0.0

model clean all

""")

print ("This model has {} contacts".format(it.contact.count()))

output:

This model has 3914 contacts

Check the properties of a contact to make sure everything makes sense. Below, a loop is made over the properties of the first contact and the property name and value are printed.

for name, value in it.contact.find(it.BallBallContact, 1).props().items():

    print ("{:15}: {}".format(name, value))

output:

kn             : 25721879.177202348
ks             : 10288751.67088094
fric           : 0.5
lin_force      : vec3(( 0.000000e+00, 0.000000e+00, 0.000000e+00 ))
lin_slip       : False
lin_mode       : 1
rgap           : 0.0
emod           : 8499999999.999999
kratio         : 2.5
dp_nratio      : 0.0
dp_sratio      : 0.0
dp_mode        : 0
dp_force       : vec3(( 0.000000e+00, 0.000000e+00, 0.000000e+00 ))
pb_state       : 3
pb_rmul        : 1.0
pb_kn          : 1917737629125.7866
pb_ks          : 767095051650.3147
pb_mcf         : 1.0
pb_ten         : 12000000.0
pb_shear       : 10000000.0
pb_coh         : 10000000.0
pb_fa          : 20.0
pb_sigma       : 0.0
pb_tau         : 0.0
pb_force       : vec3(( 0.000000e+00, 0.000000e+00, 0.000000e+00 ))
pb_moment      : vec3(( 0.000000e+00, 0.000000e+00, 0.000000e+00 ))
pb_radius      : 0.002066245
pb_emod        : 8500000000.0
pb_kratio      : 2.5
user_area      : 0.0

Now we apply the boundary conditions to particles on the upper and lower surface. The z-component of the velocity of the particles on the upper and lower boundaries are fixed, and the top particles are assigned a non-zero value.

vel = 1e-2

top_count = 0

bottom_count = 0

for b in balls.list():

    if b.extra(1) == 1:

        b.set_vel_z(-vel)

        b.set_fix(3, True)

        top_count += 1

    elif b.extra(1) == 2:

        b.set_vel_z(0)

        b.set_fix(3,True)

        bottom_count += 1

print ("{} top boundary particles and {} bottom boundary particles".format(top_count, bottom_count))

output:

112 top boundary particles and 110 bottom boundary particles

Create an array masks for the top and bottom boundaries. This is used later to calculate the average z force on the boundary particles.

top_mask = np.array([b.extra(1) == 1 for b in balls.list()])

bottom_mask  = np.array([b.extra(1) == 2 for b in balls.list()])

boundary_mask = np.logical_or(top_mask, bottom_mask)

top_mask = np.array([b.extra(1) == 1 for b in balls.list()])

bottom_mask  = np.array([b.extra(1) == 2 for b in balls.list()])

boundary_mask = np.logical_or(top_mask, bottom_mask)

r"""Create some (empty) lists to store the time and unbalanced forces

during the UCS test. """

strain = []

stress = []

top_mask = np.array([b.extra(1) == 1 for b in balls.list()])

bottom_mask  = np.array([b.extra(1) == 2 for b in balls.list()])

boundary_mask = np.logical_or(top_mask, bottom_mask)

r"""Create some (empty) lists to store the time and unbalanced forces

during the UCS test. """

strain = []

stress = []

r"""We use the matplotlib.pylab module to create a plot of the stress versus strain response of the model

during the UCS test. """

import matplotlib.pylab as plt

area = box_len**2

plt.plot(strain, stress)

plt.title("Stress vs Strain")

plt.xlabel("Strain [%]")

plt.ylabel("Stress [MPa]")

plt.grid(True)

plt.draw()

We define a function to record the values of the strain and the stress on the boundaries, and update the plot window every 100 cycles. This function to be called during the PFC cycle sequence.

def store_force(*args):

    """This function is called during |pfc| cycling to record the stress

    strain curve.

    """

    if it.cycle() % 100: return

    strain.append(vel*it.mech_age()*100.0/box_len)

    stress.append(abs(ba.force_unbal()[boundary_mask][:,2]).sum()/2.0/area/1e6)

    plt.gca().clear()

    plt.xlabel("Strain [%]")

    plt.ylabel("Stress [MPa]")

    plt.grid(True)

    plt.plot(strain, stress)

    plt.draw()

it.set_callback("store_force", 43.0)

cycle PFC

it.command("model cycle 40000")

Finally, we compute the failure strength and the Young’s modulus, and print them on the plot.

speak = np.amax(stress)

ipeak = np.argmax(stress)

epeak = strain[ipeak]/100.0

i50 = int(0.5*ipeak)

s50 = stress[i50]*1e6

e50 = strain[i50]/100.0

emod50 = s50 / e50 / 1e9

txtE = r'$E_{{50}}={:.1f} GPa$'.format(emod50)

txtS = r'$UCS = {:.1f} MPa$'.format(speak)

plt.text(0.25*100*e50,0.5*speak, txtE,fontsize=12)

plt.text(0.75*100*epeak,1.01*speak, txtS,fontsize=12)

plt.draw()

plt.savefig('p3d-test-ucs.png')

plt.close('all')

pfcpythonmodule/test3d/examples/test_ucs/p3d-test-ucs.*

Example File

The source code for the examples here is test_ucs.py.

⇄

⇐ array_interface.py | test_ucs.py ⇒

Was this helpful? ...

Updated: Sep 26, 2024