yade-users team mailing list archive

[chanaka Udaya <question692918@xxxxxxxxxxxxxxxxxxxxx>]

Question #692918 on Yade changed:
https://answers.launchpad.net/yade/+question/692918

Description changed to:
Dear Dr. Robert Caulk,

I'm trying to model three-point bending test in following paper.

"Modeling acoustic emissions inheterogeneous rocks duringtensile
fracture with the DiscreteElement Method"

https://opengeomechanics.centre-
mersenne.org/article/OGEO_2020__2__A2_0.pdf

After applying the model parameters and ran a simulation to obtained the
results as shown in figure 7a of the paper.

The outcome is different from the reported one in the paper. I
understand that maximum load ratio is obtained by the current load
divided by the maximum load during the test. If I compare the vertical
deflections in my simulation and the outcome from the paper, they are
different.( at the peak load paper has deflection value of 0.2mm while
my simulation shows 19 mm)

please see the comparison here,
https://www.dropbox.com/sh/dzqwt14z3jljpi5/AADLoCXg-
2iOMCBYXLVsPyiua?dl=0

Could you please help me to resolve this?

I have followed following steps to model the problem

1. random dense pack algorithm is used to create the particle packing
with radius range mentioned in the paper.

2. apply cpm material to the packing, used interaction detection factor
value of 1.5

3. attached three facetcylinder with diameter of 14mm ( top middle one
for applying load, remaining two at the corner of bottom to simulate the
roller support)

please see the particle assembly and supports here,
https://www.dropbox.com/s/5l09ujh20tz1z9b/assembly.PNG?dl=0

4. Then, constant velocity was imposed on the top middle facetcylinder
to apply loading.

The Minimum working example as follows,

# -*- coding: utf-8 -*-
# Code for three point bending test with crack mouth opening displacement and deflection
from yade import pack

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

from yade import pack, plot,qt,ymport


# Damping:
damp=0.7

#Material parameters
young = 30e9
poisson = .3
frictionAngle = 0.33
sigmaT = 40e6
epsCrackOnset = 3e-4
relDuctility = 30
density = 2600
 # 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 crack opening
B= height of the crack opening
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
'''

L=240e-3 #length of the sample
W=40e-3 #width of the sample
H=80e-3 #height of the sample
R=7e-3 #radius of the rollers/piston(radius of facet clyinder) all mm s (10^-3) scale
C=(12e-3+R)
A=3e-3
B=50e-3

psdmean=1.5e-3
rmin=1.5*psdmean
prepared=0
v=.01*100 # 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))
	SpherePack=pack.randomDensePack(pred,spheresInCell=500,radius=psdmean,rRelFuzz=0.25,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


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

#rollers for support
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)
#interaction detection factor
enlargeFactor=1.5
#colour the particles
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),
	PyRunner(iterPeriod=1,command='apply()',label='apply'),
	PyRunner(iterPeriod=1,command='addPlotData()',label='graph'),
	newton,
]

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

#run 1 step and set interaction detection factor to 1
O.step()
bo1s.aabbEnlargeFactor=ig2ss.interactionDetectionFactor=1.

cylinderIds=[]
for i in range(0,len(piston)):
	cylinderIds.append(piston[i])
	#print piston[i]
#print cylinderIds

'''
#for i in range(0,len(left_roller)):
	#O.bodies[left_roller[i]].state.blockedDOFs='y'
#for i in range(0,len(right_roller)):
	#O.bodies[right_roller[i]].state.blockedDOFs='y'

#for i in range(0,len(piston)):
	#O.bodies[piston[i]].state.blockedDOFs='XYZxyz'
'''
#function to apply load
def apply():
	for i in range(0,len(piston)):
		O.bodies[piston[i]].state.vel[1]=-v


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


plot.plots={'delta_z':('f_2',)}
plot.plot()


from yade import qt
qt.Controller()
qt.View()

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