[Yade-users] [Question #688421]: Radius changes in packing
New question #688421 on Yade:
https://answers.launchpad.net/yade/+question/688421
Hello everyone,
I want to define a packing of spheres with the following properties:
r_min=0.0001(m)
r_avr=0.0002(m)
r_max=0.0003(m)
number of particles = 2 (at least, It can be more than this)
cube size = 0.02 * 0.02 * 0.02 (m3)
I am using the makeCloud to do so, the problem is when I print the minimum and
maximum value of radii, I realized they have been changed because of
"internalCompaction=True".
My goal is to not decrease the sample dimensions, not increase the particles'
radii, just increase the number of particles to get the sample and target
porosity. Which obviously my code works based on increasing the radius now.
Based on the document of YADE, there are 2 options for packing which are moving
the walls or increasing the size of particles that we can implement those by
maxWallVelocity and finalMaxMultiplier functions. My desire is to not change
these 2 and instead, increase the number of particles. Is that possible to do
so?
My code:
# TRIAXIAL PROBLEM, Y IS THE VERTICAL AXIS, X IS THE RIGHT AXIS, Z IS
THE FRONT AXIS #
print ('** START **')
import numpy as np
import time
import datetime, os
start_time=time.time()
from datetime import datetime
import math
from yade import qt, export, utils
from yade import pack
# DEFINING VARIABLES #
print (' DEFINING VARIABLES ')
nRead=readParamsFromTable(
num_spheres=2,
compFricDegree = 29,
key='_triax_',
unknownOk=True
)
from yade.params import table
num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.4
compFricDegree = table.compFricDegree
finalFricDegree = 29
IP=100 # iteration period to record data and stuff
micro_record_iterPeriod=IP
ORN=3000 # O.Run Number of iterations
micro_record_enable_normal_branch=True
micro_record_float_type = np.float32
damp=0.2
thick=0
stabilityThreshold=0.01
PCPC=0.0001 # Precision of Confining Pressure Convergence
r_min=0.1*1e-3 # m
d_min=2*r_min # m
r_max=0.3*1e-3 # m
d_max=2*r_max # m
r_avr=(r_min+r_max)/2 # m
d_avr=2*r_avr # m
r_fuz=(r_max/r_avr)-1 # m
Kn=10e8*(d_avr) ### FIXME
Kt=10e8*(d_avr) ### FIXME
young=Kn/r_avr # 2 (E r1 E r2 / E r1 + E r2) >>> E = Kn/r_avr
#young=5e6 # contact stiffness
poisson=Kn/Kt # Kt/Kn
#poisson=0.3 # Kt/Kn
Ls=0.02 # mlength of specimen ### FIXME
L_REV=7*(d_avr) # m
if Ls < L_REV:
sys.exit("*** ERROR! The specimen's dimension is too samll! ***")
elif Ls==L_REV:
print ("*** This is the minimum specimen's dimension you can take! ***")
else:
print ("*** The specimen's dimension is good enough! ***")
mn,mx=Vector3(0,0,0),Vector3(Ls,Ls,Ls)
Vt=-1*1e-3 # m/s # negative sign describes the compression direction
strainRate=Vt/Ls
#strainRate=-0.01 # %/sec ### FIXME
target_strain=0.2 ### FIXME
print ("The target strain has been set to:", target_strain)
sigmaIso=-5e5 # Pa ### FIXME
particleDensity=2000#kg/m3
#particleDensity=2600
# DEFINING MATERIALS #
print (' DEFINING MATERIALS ')
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=particleDensity,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls'))
# DEFINING PACKING #
print (' DEFINING PACKING ')
walls=aabbWalls([mn,mx],thickness=thick,material='walls')
for w in walls:w.shape.radius=0
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
clumps=False
sp.makeCloud(mn,mx,r_avr,r_fuz,num_spheres,False, 0.95,seed=1)
#sp.makeCloud(mn,mx,-1,0.,num_spheres,False, 0.95,seed=1) #"seed" make the
"random" generation always the same
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
from yade import export
os.mkdir('3axresults')
os.chdir('/home/ehsan/Desktop/3axresults')
export.text('InitialPackingData')
# DEFINING TRIAXIAL TEST #
print (' DEFINING TRIAXIAL TEST ')
triax=TriaxialStressController(
maxMultiplier=1.+2e4/young,
finalMaxMultiplier=1.+2e3/young,
thickness = thick,
stressMask = 7,
internalCompaction=True,
)
# DEFINING FUNCTIONS #
print (' DEFINING FUNCTIONS ')
from yade import plot
def history():
plot.addData(
e11 = -triax.strain[0],
e22 = -triax.strain[1],
e33 = -triax.strain[2],
ev = -triax.strain[0]-triax.strain[1]-triax.strain[2],
s11 = -triax.stress(triax.wall_right_id)[0],
s22 = -triax.stress(triax.wall_top_id)[1],
s33 = -triax.stress(triax.wall_front_id)[2],
i = O.iter,
t = O.time, # virtual (yade)
Re: [Yade-users] [Question #688400]: Energy of the body is very high, defying gravity and normal stress - cpm material
Question #688400 on Yade changed:
https://answers.launchpad.net/yade/+question/688400
Akm gave more information on the question:
I will post the full script here. I will attach the drive links for the
particles because I exported those particles using a gts surface. I did
the sample preparation in a separate file using radial expansion. So the
particle positions are here-kindly download those text files.
''https://drive.google.com/drive/folders/1M2rKs-
vnR3uHC_ybv2vliYLhvO2g0RTm?usp=sharing''
The complete script is as follows:
from yade import pack, qt, export, ymport, plot
#Default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
intfactor = 1.2, #Max value of 1.2 is preferred as it can be reverted back to
1 after a step.
Shear_Vel=0.1,
S_length=0.060,
S_area=0.060*0.030,
Sig0=300e3, #300kN/m2
kn=600e3, #300kPa/mm
plotperiod=500,
topname='/home/akm/Documents/Yade/install/bin/Full-10k-Top.txt',
botname='/home/akm/Documents/Yade/install/bin/Full-10k-Bot.txt',
spname='/home/akm/Documents/Yade/install/bin/snapshots-trial/plot',
vtkname='/home/akm/Documents/Yade/install/bin/vtk-trial/s',
)
from yade.params.table import *
F_init=Sig0*S_area #Initial Normal force
#MATERIAL SPECIFICATION
mat1=CpmMat(
young = 32e9,
density=2400,
poisson = 0.30,
frictionAngle = radians(35),
epsCrackOnset = 5e-1,
relDuctility = 1.1,
sigmaT = 1.2e6,
)
mat2=CpmMat(
young = 5e9,
density=1940,
poisson = 0.30,
frictionAngle = radians(35),
epsCrackOnset = 5e-2,
relDuctility = 1.0001,
sigmaT = 1.2e1,
)
O.materials.append(mat1)
O.materials.append(mat2)
#UPPER PARTICLE IMPORT:
upbox_particles=ymport.text(topname,material=mat2, color=(1,0,0))
O.bodies.append(upbox_particles)
#LOWER PARTICLE IMPORT:
lowbox_particles=ymport.text(botname,material=mat1, color=(0,1,0))
O.bodies.append(lowbox_particles)
###
##GEOMETRY CREATION :
O.materials.append(FrictMat(frictionAngle=0,density=0,label='walls'))
#Dimension of Packing
d1,d2=Vector3(.01,.01,.01),Vector3(.07,.0783,.04) # corners of the initial
packing
xmin=d1[0]
xmax=d2[0]
ymin=d1[1]
ymax=d2[1]
zmin=d1[2]
zmax=d2[2]
ht_lower=0.03
ht_upper=0.03
asp_ht=0.0083
#Shear plate for shearing in the +X direction
P_Shear=utils.box((xmin,ymin+(ht_lower/2)-(ht_lower/20),zmin+(zmax-zmin)/2 ),
(0,(ht_lower/2),(zmax-zmin)/2 ),fixed=True, wire=True, color=(0,0,1))
O.bodies.append(P_Shear)
P_Shear.state.blockedDOFs = "xyzXYZ"
#Surrounding setup geometry
U_box_1=geom.facetBox((xmin+(xmax-xmin)/2,
ymin+ht_lower+asp_ht+(asp_ht+ht_upper)/2,zmin+(zmax-zmin)/2 ),
((xmax-xmin)/2,(asp_ht+ht_upper)/2,(zmax-zmin)/2 ), wallMask=2) # +X Boundary
constraints
U_box_2=geom.facetBox((xmin+(xmax-xmin)/2,
ymin+ht_lower+(asp_ht*2+ht_upper)/2,zmin+(zmax-zmin)/2 ),
((xmax-xmin)/2,(asp_ht*2+ht_upper)/2,(zmax-zmin)/2 ), wallMask=1) # -X Boundary
constraints
L_box=geom.facetBox((xmin+(xmax-xmin)/2,
ymin+(asp_ht+ht_lower)/2,zmin+(zmax-zmin)/2 ),
((xmax-xmin)/2,(asp_ht+ht_lower)/2,(zmax-zmin)/2 ), wallMask=6)
Cover_box=geom.facetBox(((xmin+xmax)/2,(ymin+ymax)/2,(zmin+zmax)/2
),((xmin+xmax+0.06)/2,(ymin+ymax+0.02)/2,(zmin+zmax+0.02)/2 ) ,wallMask=63)
O.bodies.append(U_box_1)
O.bodies.append(U_box_2)
O.bodies.append(L_box)
O.bodies.append(Cover_box)
idlist=[]
idlist.append(P_Shear.id)
for i in L_box:
idlist.append(i.id) #ids for the translation motion of the box
for i in O.bodies:
if isinstance(i.shape,Sphere):
if i.mat==mat1:
i.state.blockedDOFs = "yzXYZ" #blocking the mat1 movement to ensure that
there is no moment at the end point
#Adding a plate at the top:
global plate
plate_pos=max([b.state.pos[1]+(b.shape.radius) for b in O.bodies if
isinstance(b.shape,Sphere)])
plate=utils.box((xmin+(xmax-xmin)/2,plate_pos,zmin+(zmax-zmin)/2),
((xmax-xmin)/2 - xmin/60, 0, (zmax-zmin)/2),fixed=True,
wire=False,color=(0,0,1))
O.bodies.append(plate)
print('Plate added')
plate.state.blockedDOFs = "xyzXYZ"
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intfactor,label='bo1s'),Bo1_Facet_Aabb(),Bo1_Wall_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intfactor,label='ig2s'),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys(),
Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=1)],
[Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGeom_CpmPhys_Cpm()]
),
NewtonIntegrator(gravity=(0,-9.81,0),damping=0.5),
GlobalStiffnessTimeStepper(active=True,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8),
PyRunner(iterPeriod=1,command='addPlotData()'),
]
O.step()
bo1s.aabbEnlargeFactor=1
ig2s.interactionDetectionFactor=1
def addPlotData():
Sf = -O.forces.f(P_Shear.id)[0] #Negative sign since the force applied is in
-X direction
Hor_dspl = 1000*P_Shear.state.displ()[0] #Hor disp in mm
Shear_stress = Sf/S_area
Vert_dspl = 1000*plate.state.displ
Re: [Yade-users] [Question #688400]: Energy of the body is very high, defying gravity and normal stress - cpm material
Question #688400 on Yade changed:
https://answers.launchpad.net/yade/+question/688400
Akm gave more information on the question:
I will post the full script here. I will attach the drive links for the
particles because I exported those particles using a gts surface. I did
the sample preparation in a separate file using radial expansion. So the
particle positions are here-kindly download those text files.
''https://drive.google.com/drive/folders/1M2rKs-
vnR3uHC_ybv2vliYLhvO2g0RTm?usp=sharing''
The complete script is as follows:
from yade import pack, qt, export, ymport, plot
#Default parameters or from table
readParamsFromTable(noTableOk=True, # unknownOk=True,
intfactor = 1.2, #Max value of 1.2 is preferred as it can be reverted
back to 1 after a step.
Shear_Vel=0.1,
S_length=0.060,
S_area=0.060*0.030,
Sig0=300e3, #300kN/m2
kn=600e3, #300kPa/mm
F_init=Sig0*S_area,
#Initial Normal force
plotperiod=500,
topname='/home/akm/Documents/Yade/install/bin/Full-10k-Top.txt',
botname='/home/akm/Documents/Yade/install/bin/Full-10k-Bot.txt',
spname='/home/akm/Documents/Yade/install/bin/snapshots-trial/plot',
vtkname='/home/akm/Documents/Yade/install/bin/vtk-trial/s',
)
from yade.params.table import *
#MATERIAL SPECIFICATION
mat1=CpmMat(
young = 32e9,
density=2400,
poisson = 0.30,
frictionAngle = radians(35),
epsCrackOnset = 5e-1,
relDuctility = 1.1,
sigmaT = 1.2e6,
)
mat2=CpmMat(
young = 5e9,
density=1940,
poisson = 0.30,
frictionAngle = radians(35),
epsCrackOnset = 5e-2,
relDuctility = 1.0001,
sigmaT = 1.2e1,
)
O.materials.append(mat1)
O.materials.append(mat2)
#UPPER PARTICLE IMPORT:
upbox_particles=ymport.text(topname,material=mat2, color=(1,0,0))
O.bodies.append(upbox_particles)
#LOWER PARTICLE IMPORT:
lowbox_particles=ymport.text(botname,material=mat1, color=(0,1,0))
O.bodies.append(lowbox_particles)
###
##GEOMETRY CREATION :
O.materials.append(FrictMat(frictionAngle=0,density=0,label='walls'))
#Dimension of Packing
d1,d2=Vector3(.01,.01,.01),Vector3(.07,.0783,.04) # corners of the initial
packing
xmin=d1[0]
xmax=d2[0]
ymin=d1[1]
ymax=d2[1]
zmin=d1[2]
zmax=d2[2]
ht_lower=0.03
ht_upper=0.03
asp_ht=0.0083
#Shear plate for shearing in the +X direction
P_Shear=utils.box((xmin,ymin+(ht_lower/2)-(ht_lower/20),zmin+(zmax-zmin)/2 ),
(0,(ht_lower/2),(zmax-zmin)/2 ),fixed=True, wire=True, color=(0,0,1))
O.bodies.append(P_Shear)
P_Shear.state.blockedDOFs = "xyzXYZ"
#Surrounding setup geometry
U_box_1=geom.facetBox((xmin+(xmax-xmin)/2,
ymin+ht_lower+asp_ht+(asp_ht+ht_upper)/2,zmin+(zmax-zmin)/2 ),
((xmax-xmin)/2,(asp_ht+ht_upper)/2,(zmax-zmin)/2 ), wallMask=2) # +X
Boundary constraints
U_box_2=geom.facetBox((xmin+(xmax-xmin)/2,
ymin+ht_lower+(asp_ht*2+ht_upper)/2,zmin+(zmax-zmin)/2 ),
((xmax-xmin)/2,(asp_ht*2+ht_upper)/2,(zmax-zmin)/2 ), wallMask=1) # -X
Boundary constraints
L_box=geom.facetBox((xmin+(xmax-xmin)/2,
ymin+(asp_ht+ht_lower)/2,zmin+(zmax-zmin)/2 ),
((xmax-xmin)/2,(asp_ht+ht_lower)/2,(zmax-zmin)/2 ), wallMask=6)
Cover_box=geom.facetBox(((xmin+xmax)/2,(ymin+ymax)/2,(zmin+zmax)/2
),((xmin+xmax+0.06)/2,(ymin+ymax+0.02)/2,(zmin+zmax+0.02)/2 ) ,wallMask=63)
O.bodies.append(U_box_1)
O.bodies.append(U_box_2)
O.bodies.append(L_box)
O.bodies.append(Cover_box)
idlist=[]
idlist.append(P_Shear.id)
for i in L_box:
idlist.append(i.id) #ids for the translation motion
of the box
for i in O.bodies:
if isinstance(i.shape,Sphere):
if i.mat==mat1:
i.state.blockedDOFs = "yzXYZ" #blocking the mat1
movement to ensure that there is no moment at the end point
#Adding a plate at the top:
global plate
plate_pos=max([b.state.pos[1]+(b.shape.radius) for b in O.bodies if
isinstance(b.shape,Sphere)])
plate=utils.box((xmin+(xmax-xmin)/2,plate_pos,zmin+(zmax-zmin)/2),
((xmax-xmin)/2 - xmin/60, 0, (zmax-zmin)/2),fixed=True,
wire=False,color=(0,0,1))
O.bodies.append(plate)
print('Plate added')
plate.state.blockedDOFs = "xyzXYZ"
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intfactor,label='bo1s'),Bo1_Facet_Aabb(),Bo1_Wall_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intfactor,label='ig2s'),Ig2_Facet_Sph
Re: [Yade-users] [Question #688388]: On the interaction range of the JCFPm
Question #688388 on Yade changed: https://answers.launchpad.net/yade/+question/688388 Status: Open => Answered Robert Caulk proposed the following answer: Hello, Interaction range is not handled by JCFpm, it is handled by the collider [1] and the interaction detection [2]. Here is an example usage [3]. Cheers, Robert [1]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Bo1_Sphere_Aabb.aabbEnlargeFactor [2]https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Ig2_Sphere_Sphere_ScGeom.interactionDetectionFactor [3]https://gitlab.com/yade-dev/trunk/blob/master/examples/concrete/brazilian.py#L62 -- 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
Re: [Yade-users] [Question #687871]: Outputs of the Triaxial code by Bruno Chareyre
Question #687871 on Yade changed: https://answers.launchpad.net/yade/+question/687871 Status: Open => Expired Launchpad Janitor expired the question: This question was expired because it remained in the 'Open' state without activity for the last 15 days. -- 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
