The source code for this example is test_ucs.py
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
it.command("python-reset-state false")
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 domain extent -1 1 -1 1 -1 1
model 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
model 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
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.
| Was this helpful? ... | PFC 6.0 © 2019, Itasca | Updated: Nov 19, 2021 |