Question #688400 “Energy of the body is very high, defying gravi...” : Questions : Yade

Revision history for this message

Akm (arunkumarmurali) said on 2020-02-01:

#2

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()[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#####

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)