yade-users team mailing list archive

[Akm <question688400@xxxxxxxxxxxxxxxxxxxxx>]

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_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()[1]	
	Force_Stif=kn*Vert_dspl*S_area  					#Force(N)=Stiffness(N/m2/mm)*Displacement(mm)*Area(m2)
	Force_Reqd=Force_Stif+F_init						#Final force on plate = Force due to CNS + Intial force on the plate(SigmaO)
	Forceonplate=O.forces.f(plate.id)[1]
	if Forceonplate<Force_Reqd:
		plate.state.vel=(0,-2,0)
	elif Forceonplate>Force_Reqd:
		plate.state.vel=(0,1,0)
	Norm_stress_applied=(Forceonplate)/S_area
	Norm_stress_calculated=Force_Reqd/S_area

	yade.plot.addData(
		Hor_Disp = Hor_dspl,
		Shear_stress = Shear_stress,
		Ver_Disp = Vert_dspl,
		Norm_stress_applied = Norm_stress_applied,
		Norm_stress_calculated=Norm_stress_calculated
	)

# note the space in 'i ' so that it does not overwrite the 'i' entry
plot.plots={'Hor_Disp':('Shear_stress'),'Hor_Disp  ':('Norm_stress_applied','Norm_stress_calculated'),}
plot.plot()

#EXTRA ENGINES:
O.engines=O.engines+[TranslationEngine(ids=idlist,translationAxis=[1,0,0], velocity=Shear_Vel)]
O.engines=O.engines+[SnapshotEngine(fileBase=spname,iterPeriod=plotperiod)]
#O.engines=O.engines+[VTKRecorder(fileName=vtkname,recorders=['cpm','all'],iterPeriod=500)]

yade.qt.View()
O.saveTmp()

#####Thanks in advance guys#####

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.