New question #706821 on Yade:
https://answers.launchpad.net/yade/+question/706821
Hi! I want to model the particle breakage of clumps during the one-dimensional
compression in a cylindrical boundary. My steps are as followed:
Step 1: I import 800 GTS files into Yade and generated one-on-one 800 clumps.
That is to say that each GTS corresponds to a unique clump.
Step 2: I create a cylindrical boundary to cover those clumps and provide two
plates (bottom one fixed, top one applied with force). Note that the force
applied on the top plate is a multi-load, which varies with the strain of the
sample.
Step 3: I assign JCFpm material for each clump and material of steel for the
boundary and plates.
Step 4: record what we want to know and implement simulation.
Here, my questions are as follows:
1) Is my whole procedure appropriate to simulate the one-dimensional
compression while considering particle breakage?
2) For Step 2, I would like to know how to assign multi-loads is correct. It
seems currently the top plate is not quasi-static during loading. Do I need to
clear the speed for every number of calculation steps?
3) For the whole simulation, how do we know how many clumps are broken and can
we record the isolated clumps at the final?
The MWM is presented below (I don't know how to provide the gts files without
an external link, so maybe randomly generating some clumps can replace them.):
from __future__ import print_function
from yade import pack, ymport
import gts, os.path, locale, plot, random
import numpy as np
locale.setlocale(
locale.LC_ALL, 'en_US.UTF-8'
) #gts is locale-dependend. If, for example, german locale is used,
gts.read()-function does not import floats normally
############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
# The following lines are used parameter definitions
readParamsFromTable(
# directory
dir_gts = '/home/kyle/Yade/install/bin/gts', # directory of storing gts
files
dir_vtk = '/home/kyle/Yade/install/bin/vtk', # directory of storing vtk
files
dir_txt = '/home/kyle/Yade/install/bin/txt', # directory of storing txt
files
fileType = '.gts', # type of file
# type of generating spheres ['hexa','dense']
genType = 'hexa',
# material parameters
young_w = 210e9, # Young's modulus of plates and walls
young_g = 180e6, # Young's modulus of grains [1]
den_w = 7874, # density of plates and walls
den_g = 1162, # density of grains
poi_w = 0.3, # possion's ratio of plates and walls
poi_g = 0.25, # possion's ratio of grains
friAngle_w = 15, # friction angle of plates and walls
friAngle_g = 19, # friction angle of grains
# Parameters of cylindrical walls
x_cyl = 0.0106, # x-coordinate of center of cylindrical walls
y_cyl = 0.01005, # y-coordinate of center of cylindrical walls
z_cyl = 0, # z-coordinate of center of cylindrical walls
r_cyl = 0.0097, # radius of the cylinder
h_cyl = 0.04151, # height of the cylinder
n_cyl = 100, # number of boxes forming cylinder
t_cyl = 0.0001, # thickness of the boxes
#Parameters of lower and upper plates (square box)
w_plate = 0.015, # width of plates
t_plate = 0.0001, # thickness of plates
# applied stress
stress = [12500, 25000, 50000, 100000, 200000, 300000, 400000, 800000,
1600000, 3200000], # target applied stress
# target displacement
displacement = [3.6e-5, 6.8e-5, 0.00014, 0.00028, 0.00047, 0.000642,
0.000796, 0.001342, 0.002441, 0.003989] # target displacement
)
from yade.params.table import *
# create materials for spheres and walls
wallMat = O.materials.append(FrictMat(young = young_w, poisson = poi_w,
frictionAngle = friAngle_w, density = den_w, label = 'walls'))
#sphereMat = O.materials.append(FrictMat(young = young_g, poisson = poi_g,
frictionAngle = friAngle_g, density = den_g, label = 'spheres'))
def sphereMat():
return JCFpmMat(
type=1,
young = young_g,
frictionAngle = radians(friAngle_g),
density = den_g,
poisson = poi_g,
tensileStrength = 1e6,
cohesion=1e6,
jointNormalStiffness=1e7,
jointShearStiffness=1e7,
jointCohesion=1e6,
jointFrictionAngle=radians(20),
jointDilationAngle=radians(0)
) ## Rq: density needs to be adapted as porosity of real rock is
different to granular assembly due to difference in porosity
(utils.sumForces(baseBodies,(0,1,0))/(Z*X) should be equal to Gamma*g*h with
h=Y, g=9.82 and Gamma=2700 kg/m3
###############################################
### IMPORTING GRAINS AND CREATING CLUMPS ###
###############################################
##### a function that searching all .gts files
def SearchFiles(directory, ftype):
fileList = []
for root, subDirs, files in os.walk(directory):
for fileName in files:
if fileName.endswith(ftype):
fileList.append(os.path.join(root, fileName))
return fileList
#####
# import gts and obatin file list
fileList = SearchFiles(dir_gts,fileType)
##### a function that determining the radius of spheres
def createSphere(pred, ndivmin, spacing, genType):
#dim = pred.dim()
#minDim = min(dim[0], dim[1], dim[2])
aabb = pred.aabb()
dim0 = aabb[1][0] - aabb[0][0]
bodyList = []
ndiv = ndivmin
if genType == 'hexa':
while len(bodyList)==0:
radius = dim0 / ndiv # get some characteristic
dimension, use it for radius
bodyList = O.bodies.append(pack.regularHexa(pred,
radius = radius, gap = -radius/4, color =
[random.random(),random.random(),random.random()], material = sphereMat))
ndiv = ndiv + spacing
else:
while len(bodyList)==0:
radius = dim0 / ndiv # get some characteristic
dimension, use it for radius
sp =
pack.randomDensePack(pred,radius=radius,returnSpherePack=True, color =
[random.random(),random.random(),random.random()], material = sphereMat)
bodyList = sp.toSimulation()
ndiv = ndiv + spacing
return bodyList
#####
# import gts and create clump for each .gts file
sphereList = []
for gtsName in fileList:
surf = gts.read(open(gtsName))
if surf.is_closed():
print(gtsName)
pred = pack.inGtsSurface(surf)
bodyList = createSphere(pred, 20, 2, genType)
idClump = O.bodies.clump(bodyList)
del surf
del pred
del bodyList
O.bodies.updateClumpProperties()
print('Number of bodies:', len(O.bodies))
###################################################################
### CREATING CYLINDRICAL WALLS, BOTTOM PLATE, AND TOP PLATE ###
###################################################################
##### a function that creating cylindrical walls
def cylSurf(center,radius,height,nSegments=12,thick=0,**kw):
"""creates cylinder made of boxes. Axis is parallel with z+
center ... center of bottom disc
radius ... radius of cylinder
height ... height of cylinder
nSegments ... number of boxes forming cylinder
thick ... thickness of the boxes
**kw ... keywords passed to box function (e.g. material, color ...)
"""
# just converts (a,b,c) to Vector3(a,b,c) to enable "center + ..."
center = Vector3(center)
# nSegments angles to define individual boxes, from 0 to 2*pi
angles = [i*2*pi/nSegments for i in range(nSegments)]
# centers of boxes. Vector 3(radius,0,0) rotated by angle 'a'
# z coordinate is .5*height
# center + shift = center of corresponding box
pts = [center + Vector3(radius*cos(a),radius*sin(a),.5*height) for a in
angles] # centers of plates
ret = []
for i,pt in enumerate(pts): # for each box center create corresponding box
# "length" of the box along circumference, cylinder circumference /
nSegments
l = pi*radius/nSegments
# extents, half dimensions along x,y,z axis
es = (.5*thick,l,.5*height)
# by default box is axis-aligned, we need to rotate it along z axis
(0,0,1) by angle i*2*pi/nSegments
ori = Quaternion((0,0,1),i*2*pi/nSegments)
# have a look at utils.box documentation. pt=center of box, es=extents
(half dimension), ori=orientation
ret.append(box(pt,es,ori,**kw))
return ret
#####
# create cylindrical walls and fix them
surf = cylSurf([x_cyl, y_cyl, z_cyl], r_cyl, h_cyl, nSegments = n_cyl, thick =
t_cyl, material = wallMat, wire=True)
O.bodies.append(surf)
for b in surf: b.state.blockedDOFs = 'xyzXYZ'
for b in surf: b.state.vel = (0, 0, 0)
# create the bottom plate and fix it
bottom_plate = O.bodies.append(box([x_cyl, y_cyl, z_cyl], (w_plate, w_plate,
t_plate), material = wallMat))
O.bodies[bottom_plate].state.vel = (0, 0, 0)
O.bodies[bottom_plate].state.blockedDOFs = 'xyzXYZ'
# create the top plate
top_plate=O.bodies.append(box([x_cyl, y_cyl, h_cyl], (w_plate, w_plate,
t_plate), material = wallMat))
# apply initial force
force = np.dot(stress, w_plate * w_plate) # corresponding applied force
global deltaF
deltaF = force[0]
O.forces.setPermF(top_plate,(0, 0, -deltaF))
# target displacement
global numIter
numIter = 0
global maxIter
maxIter = 9
global dispT
dispT = displacement[numIter]
global disp
disp = 0
global dispE
dispE = displacement[-1]
# calculate initial void ratio
vol0 = h_cyl * pi * r_cyl * r_cyl
volG0 = sum([b.state.mass for b in O.bodies if b.isClump])/den_g
volV0 = vol0 - volG0
e0 = volV0/volG0
print(volG0, e0)
print('numIter', 'dispT', 'disp', 'disp/dispT', 'velTopPlate')
############################
### DEFINING ENGINES ###
############################
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1),
Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1),
Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys(),
Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
[Law2_ScGeom_FrictPhys_CundallStrack(),
Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True)]),
GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.8),
PyRunner(iterPeriod = 10,command = "checkDisplacement(displacement,
force, h_cyl)"),
VTKRecorder(iterPeriod = 500, recorders = ['jcfpm', 'cracks',
'spheres', 'boxes', 'stress', 'id', 'clumpId', 'coordNumber', 'colors'],
fileName = dir_vtk + '/'),
NewtonIntegrator(gravity = [0, 0, -9.81], damping = 0.4)
]
O.dt= 0.1 * PWaveTimeStep()
# check whether the target displacement has already met. if so, update applied
stress
def checkDisplacement(displacement, force, h_cyl):
global numIter
global maxIter
global disp
global dispT
global dispE
global deltaF
# at the very start, applied stress is the first of value of array
'stress'
t = O.time
H = O.bodies[top_plate].state.pos[2]
disp = h_cyl - H
velTopPlate = O.bodies[top_plate].state.vel[2]
if disp > dispT:
if numIter < maxIter:
numIter = numIter + 1
dispT = displacement[numIter]
#deltaF = force[numIter] - force[numIter - 1]
O.forces.setPermF(top_plate,(0, 0, -force[numIter]))
else:
print(numIter, dispT, disp, disp/dispT, velTopPlate,
O.forces.permF(top_plate)[2], O.forces.f(top_plate)[2])
f = open(dir_txt + '/data.txt', 'a')
f.write('{} {} {} {} {} {} {}\n'.format(O.iter,
numIter, dispT, disp, velTopPlate, O.forces.permF(top_plate)[2], H))
f.close()
plot.saveDataTxt(O.tags['d.id'] + '.txt')
O.pause()
print(numIter, dispT, disp, disp/dispT, velTopPlate,
O.forces.permF(top_plate)[2], O.forces.f(top_plate)[2])
plot.addData(t=t,disp=disp)
f = open(dir_txt + '/data.txt', 'a')
f.write('{} {} {} {} {} {} {}\n'.format(O.iter, numIter, dispT, disp,
velTopPlate, O.forces.permF(top_plate)[2], H))
f.close()
##########################
### RUN SIMULATION ###
##########################
plot.plots = {'t':('disp')}
plot.plot()
O.run()
--
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
