Plate for three point bending load application

Asked by chanaka Udaya on 2020-08-12

Dear all,

I'm trying to simulate three-point bending test in Yade.
I have prepared a sample with random dense pack.

In order to apply the load as a piston, I need to have a plate with 3mm in width and 50mm length parallel to xz plane in the middle of the beam.

and another two plates as rollers with the same dimensions at the bottom left and right locations?

How can I add a plate with a specific dimension?

Previously I have used facetcylinder to apply load as well as to use as rollers.

Now I need to use a single plate instead of this facet cylinder. could you please tell me how to do that?

this is my existing code,

# -*- coding: utf-8 -*-
# Code for three point bending test with crack mouth opening displacement and deflection
#Developed by S.M.C.U. Senanayake, Department of Civil Engineering, Monash University, Australia

from yade import pack

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

from yade import pack, plot,qt,ymport
import os
import csv

# Damping:
damp=0.7
stabilityThreshold=0.001 # we test unbalancedForce against this value in different loops (see below)
# material parameters
#young = 6.75e8
#poisson = .2
#frictionAngle = 0.6
#sigmaT = 1.5e6
#epsCrackOnset = 16.4e-3
#relDuctility = 1.05
#density = 2600

young = 3e9
poisson = .3
frictionAngle = 0.5
sigmaT = .3e6
epsCrackOnset = .9e-3
relDuctility = 30
density = 2600
# prestress
# axial strain rate
rate=0.01 # loading rate (strain rate)
 # stress path
'''
# 3 point bending test sample preparation parameters
L= length of the sample
W= width of the sample
H= height of the sample
C=clearance from two ends(canteliver length
A= width of the notch
B= height of the notch
R=radius of the rollers/piston(radius of facet clyinder) all mm s (10^-3) scale
Length along x axis, height along y axis and width along z axis
'''
#mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

L=500e-3
W=50e-3
H=100e-3
C=50e-3
A=3e-3
B=50e-3
R=7e-3
psdmean=3e-3
rmin=1.5*psdmean
prepared=0
v=.1 # applied velocity on piston

concMat=O.materials.append(CpmMat(young=young,poisson=poisson,density=density,damLaw=1,frictionAngle=frictionAngle,
sigmaT=sigmaT,epsCrackOnset=epsCrackOnset,neverDamage=False,label='concrete',relDuctility=relDuctility,))

frictMat = O.materials.append(CpmMat(young=100*young,poisson=poisson,frictionAngle=0,sigmaT=0,epsCrackOnset=epsCrackOnset,relDuctility=1.0001,neverDamage=True,density=3000,label='walls'))

if prepared<0.5:
 pred=pack.inAlignedBox((0,0,0),(L,H,W))-pack.inAlignedBox((.5*(L-A),0,0),(.5*(L+A),B,W))
 #SpherePack=pack.randomDensePack(pred,spheresInCell=11000,radius=psdmean,rRelFuzz=0.2,color=(0,0,1),material=concMat,returnSpherePack=True)
 SpherePack=pack.randomDensePack(pred,spheresInCell=500,radius=psdmean,rRelFuzz=0.2,color=(0,0,1),material=concMat,returnSpherePack=True)
 sand=SpherePack.toSimulation() # assign all particles to sand object
else:
 O.bodies.append(ymport.textExt('tpb_sample.txt',format='x_y_z_r',material=concMat))
 print ('I have previously prepared packing!!')
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.shape.color = (0.2,0.6,0.6)# give a colour
# In order to extract the ids of each particle in the left and right hand corner of the notch and make a list out of them
leftids=[]
for i in O.bodies:
 if i.state.pos[0]>(0.5*(L-A)-2*rmin) and i.state.pos[0]<(0.5*(L-A)):
  if i.state.pos[1]>0 and i.state.pos[1]<(2*rmin):
   i.shape.color=(1,0,0)
   leftids.append(i.id)
   #left=[O.bodies[i] in sand]

rightids=[]
c=0
for i in O.bodies:
 if i.state.pos[0]>(0.5*(L+A)) and i.state.pos[0]<(0.5*(L+A)+2*rmin):
  if i.state.pos[1]>0 and i.state.pos[1]<(2*rmin):
   i.shape.color=(1,0,0)
   c=c+1
   print i
   print i.id
   rightids.append(i.id)
   #right=[O.bodies[i] in sand]
print c
print rightids
#print right
# assign to a lh object all particles in the left hand side interested region of the notch
lh = [O.bodies[s] for s in sand if (O.bodies[s].state.pos[0]>(0.5*(L-A)-2*rmin) and O.bodies[s].state.pos[0]<(0.5*(L-A)) and O.bodies[s].state.pos[1]>0 and O.bodies[s].state.pos[1]<(2*rmin))]
rh=[O.bodies[s] for s in sand if (O.bodies[s].state.pos[0]>(0.5*(L+A)) and O.bodies[s].state.pos[0]<(0.5*(L+A)+2*rmin) and O.bodies[s].state.pos[1]>0 and O.bodies[s].state.pos[1]<(2*rmin))]
for i in rh:
 print i.id
# assign global variable to total numbe rof particles in each block
global left_part
left_part=len(leftids)
global right_part
right_part=len(rightids)

piston=O.bodies.append(geom.facetCylinder(center=(.5*L,H+R,.5*W),radius=R,height=W,orientation=Quaternion((0, 0, 1), pi/2),segmentsNumber=20,wire=False,material=frictMat)) #piston radius is 15mm

left_roller=O.bodies.append(geom.facetCylinder(center=(C,-R,.5*W),radius=R,height=W,orientation=Quaternion((0, 0, 1), pi/2),segmentsNumber=20,wire=False,material=frictMat))

right_roller=O.bodies.append(geom.facetCylinder(center=(L-C,-R,.5*W),radius=R,height=W,orientation=Quaternion((0, 0, 1), pi/2),segmentsNumber=20,wire=False,material=frictMat))
newton=NewtonIntegrator(damping=damp)
enlargeFactor=1.5
for b in O.bodies:
 if isinstance(b.shape,Sphere):
  b.shape.color = (0.2,0.6,0.6)
print len(O.bodies)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop(
  [
   Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ig2ss'),
   Ig2_Facet_Sphere_ScGeom()
  ],
  [
   Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=O.iter+5),
   Ip2_FrictMat_CpmMat_FrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys(),
  ],
  [
   Law2_ScGeom_CpmPhys_Cpm(),
   Law2_ScGeom_FrictPhys_CundallStrack(),
  ],
 ),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 #TranslationEngine(ids=[piston.id],translationAxis=[0,0,1],velocity=1),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3),
 #VTKRecorder(fileName='/home/ud/vtkdata/3pb_test', recorders=['all'], iterPeriod=100),
 PyRunner(iterPeriod=1,command='apply()',label='apply'),
 PyRunner(iterPeriod=1,command='addPlotData()',label='graph'),
 #PyRunner(iterPeriod=50,command='pr()',label='ph'),
 newton,
 #CpmStateUpdater(iterPeriod=100),
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
qtr=yade.qt.Renderer()
qtr.bgColor=(1,1,1)
yade.qt.Controller(), yade.qt.View()
O.step()
bo1s.aabbEnlargeFactor=ig2ss.interactionDetectionFactor=1.
cylinderIds=[]
for i in range(0,len(piston)):
 cylinderIds.append(piston[i])
 #print piston[i]
#print cylinderIds

def apply():
 for i in range(0,len(piston)):
  O.bodies[piston[i]].state.vel[1]=-v
# a function to find crack mouth opening displacement
def cmod():
 sumlh=0
 sumrh=0
 for i in lh:
  dp_left=i.state.pos[0]-i.state.refPos[0]# left displacment
  sumlh+=dp_left # addig all the displacement in x direction components

 for j in rh:
  dp_right=j.state.pos[0]-j.state.refPos[0]
  sumrh+=dp_right

 d1=sumlh/left_part # average displacment in left hand side of the notch
 d2=sumrh/right_part
 cmodi=abs(d1)+abs(d2) # cmod by adding right ahand side and left hand side
 return (cmodi)

def addPlotData():
 f_2=utils.sumForces(ids=cylinderIds,direction =(0,1,0))# measure total force induced on the cylinder
 delta_z=v*O.time
 crackopen=cmod()
 #print f_2
 plot.addData(t = O.time,f_2 = f_2,delta_z=delta_z,crackopen=crackopen)

#O.engines=O.engines+[PyRunner(command='pr()',iterPeriod=50)]
def pr():
 presentcohesive_count = 0
 for i in O.interactions:
  if hasattr(i.phys, 'isCohesive'):
   if i.phys.isCohesive == 0:
    presentcohesive_count+=1
    a=i.id1
    b=i.id2
    O.bodies[a].shape.color = (1,0,0)
    O.bodies[b].shape.color = (1,0,0)
    #print (i.id1, "-", i.id2)
 noncohesive_count=presentcohesive_count

#O.engines=O.engines+[TranslationEngine(ids=[piston.id],translationAxis=[0,0,1],velocity=1)]
#plot.plots={'delta_z':('f_2',), 'crackopen':('f_2 ')}
plot.plots={'crackopen':('f_2',)}
plot.plot()

#print (len(O.bodies))
from yade import qt
qt.Controller()
qt.View()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2020-08-15
Last query:
2020-08-15
Last reply:
2020-08-13
Jan Stránský (honzik) said : #1

Hello,
you can use two facets to form a rectangle.
Or a box instead of facets.
Is this what you are looking for?

> Now I need to use a single plate instead of this facet cylinder

what is the reason for this need? I don't expect much difference.. Using randomDensePack, the surface is unrealistically rough anyway..

cheers
Jan

chanaka Udaya (chanaka-udaya) said : #2

Dear Dr. Jan,

yes I can use a box.many thanks!

>>Using randomDensePack, the surface is unrealistically rough anyway..

That's the problem. initially I wanted to prepare sample with radius expansion method. and I have used 6 walls to create the beam with radius expansion method. But, to have a rectangular notch with 3mm width in the middle was difficult. notch was irregular shape.

the procedure as follows,
After creation of the beam using radius expansion method. I have exported particle positions and radius in to a separate file. Next, imported into the new file with particle location and radius and appended in to the new simulation using yade.ymport function.Afterwards created a bounding box with 3mm with and 50mm height and 50m length. the particles which have their centrod inside the bounding box , deleted to obtain the notch.

Is there a way to create a rectangular notch in the middle with nearly smoother boundary?

Thanks
Chanaka

Best Jan Stránský (honzik) said : #3

> Is there a way to create a rectangular notch in the middle with nearly smoother boundary?

not really I am afraid.. You can modify the surface e.g. by changing size of the boundary particles to form a smoother surface, but then the notch is different from the rest of the specimen..

It will always be a tradeoff between physical accuracy, resources etc.
E.g. you can use half-sized particles to make the surface "twice smoother" (in the limit of infinitesimally small particles, the surface is smooth), but the computation time increase...

cheers
Jan

chanaka Udaya (chanaka-udaya) said : #4

Thnaks. I have used radius expansion method with finer mesh.

chanaka Udaya (chanaka-udaya) said : #5

Thanks Jan Stránský, that solved my question.