# Three-point bending unnotched beam

Asked by chanaka Udaya on 2020-09-15

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)

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

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

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

#rollers for support

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'),
newton,
]

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
qtr.bgColor=(1,1,1)

#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'
'''
def apply():
for i in range(0,len(piston)):
O.bodies[piston[i]].state.vel=-v

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.plots={'delta_z':('f_2',)}
plot.plot()

qt.Controller()
qt.View()

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2020-09-15
2020-09-17
 Robert Caulk (rcaulk) said on 2020-09-15: #1

Hello,

I am unfamiliar with CPM and its parameters, such as relductility, crackonset, etc. The model used in the paper you reference is JCFpm . I do not expect CPM and JCFpm to exhibit the same behavior at all. If you wish to reproduce the curves in the paper you cite, you will need to use JCFpm instead.

We have a script associated with the publication you reference and you can find it here for now  (later, it will be in the master branch).

Cheers,

Robert

 chanaka Udaya (chanaka-udaya) said on 2020-09-17: #2

Dear Dr. Robert Caluk,

Many thanks for sharing the code.

I was able to run the simulation using the provided script with JCFPM model.(result-https://www.dropbox.com/s/7tx32am1b33pbwb/JCFPM-fail-matched.PNG?dl=0 )

Later, I have used CPM model with same particle packing file (240x80mmBeam_1.5mmrad.sphres)provided by you. (result- https://www.dropbox.com/s/cbfcut6up5dvrcp/cpm-same%20packing.PNG?dl=0 )

It worked !

Then, I tried to use randomdensepackfunction to generate the packing instead of your packing file.

This is the line I changed
-------------------------------------------------
numSpheres=10000
mn,mx=Vector3(-0.04, -.140,-0.020),Vector3(0.04, .140, 0.020)
pred = pack.inAlignedBox(mn, mx)
sp = pack.randomDensePack(pred, radius=0.0015, spheresInCell=numSpheres, rRelFuzz=0.25, material='JCFmat', returnSpherePack=True)
sp.toSimulation()
----------------------------------------------

Then I tried to run the program again with above changes and after some time program stopped with the error "Core dumped-Aborted")
what could be the reason.

The working CPM model script for attached image as follows,

-------------------------------------------------------------------------------------------------------------------------------

from yade import ymport, utils, pack, export
from pylab import *
import os
import math
import numpy as np
import numpy.linalg as la
import time
timeStr = time.strftime('%m-%d-%Y')

conf = 0,
weibullCutOffMax=10,
weibullCutOffMin=0.1,
xSectionShape = 4,
xSectionScale = 0,

)

mn,mx=Vector3(-0.04, -.140,-0.020),Vector3(0.04, .140, 0.020)

young = 30e9
poisson = 0.3
targetPorosity = 0.55
density = 5000
relDuctility=30
finalFricDegree = 19
sigmaT = 10e6
epsCrackOnset=(10e6/30e9)
iterper=1000
cohesion=40e6
dispVel = -0.02 # m/s
highFric = 45

identifier='shape-'+str(xSectionShape)+'-'+timeStr
outputDir='out_'+identifier
if not os.path.exists(outputDir):
os.mkdir(outputDir)
if not os.path.exists('txt'):
os.mkdir('txt')

output = './'+outputDir+'/'+identifier

sigmaT=cohesion,epsCrackOnset=epsCrackOnset,neverDamage=False,label='JCFmat',relDuctility=0,))
#### preprocessing to get dimensions

dim=utils.aabbExtrema()
xinf=dim
xsup=dim
X=xsup-xinf
yinf=dim
ysup=dim
Y=ysup-yinf
zinf=dim
zsup=dim
Z=zsup-zinf

r=X/15.

O.reset()

# Rigid blocks
block1 = O.bodies.append(geom.facetCylinder(center=(xinf-0.971*r,yinf+2*r,0),radius=r,height=Z,orientation=Quaternion((0, 0, 0), 90),segmentsNumber=20,wire=False,material=wallMat))

block2 = O.bodies.append(geom.facetCylinder(center=(xinf-0.915*r,ysup-2*r,0),radius=r,height=Z,orientation=Quaternion((0, 0, 0), 90),segmentsNumber=20,wire=False,material=wallMat))

p0=O.bodies[piston].state.pos

sigmaT=cohesion,epsCrackOnset=epsCrackOnset,neverDamage=False,label='JCFmat',relDuctility=0,))
# Specimen

AEfile = 'AEcounts_'+identifier+'.txt'

f = open('txt/'+AEfile, 'w')
f.write('time iteration AEcount P deflection pistonDisp\n')
f.close

# remove rigid block DOFs
for i in block1:
O.bodies[i].state.blockedDOFs='xyzXYZ'

for i in block2:
O.bodies[i].state.blockedDOFs='xyzXYZ'

# set engine list
O.engines=[
ForceResetter(),
InteractionLoop(
[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=10, label='jcf')],
[Law2_ScGeom_CpmPhys_Cpm(label='interactionLaw'),Law2_ScGeom_FrictPhys_CundallStrack()]
),

GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.4, defaultDt=0.1*utils.PWaveTimeStep()),
TranslationEngine(ids=piston,translationAxis=(1,0,0), velocity=dispVel, label='pistonEngine'),
TranslationEngine(ids=block1,translationAxis=(1,0,0), velocity=0, label='block1Engine'),
TranslationEngine(ids=block2,translationAxis=(1,0,0), velocity=0, label='block2Engine'),
NewtonIntegrator(damping=0.2)
]

# options for AE model
#interactionLaw.clusterMoments=True
#interactionLaw.neverErase=True
#interactionLaw.useIncrementalForm=True
#interactionLaw.useStrainEnergy = True

O.dt = 0.004 * utils.PWaveTimeStep()

# collect data for active plotting and post processing

# engage blocks to ensure force balance and remove dynamics of system

def ensureBlockEngagement():
Pblock1, Pblock2 = checkBlockEngagement()
P = getPistonForce()
maxVel = 0.01
if Pblock1 != P/2.:
multiplier = (P/2. - Pblock1)/(P/2.)
block1Engine.velocity = multiplier*maxVel
if Pblock2 != P/2.:
multiplier = (P/2. - Pblock2)/(P/2.)
block2Engine.velocity = multiplier*maxVel

def getPistonForce():
P=0
for i in piston:
P += la.norm(O.forces.f(i))
return P

def getPistonDisp():
pistonDisp = O.bodies[piston].state.displ().norm()
block1Disp = O.bodies[block1].state.displ().norm()
block2Disp = O.bodies[block2].state.displ().norm()
return pistonDisp + (block1Disp + block2Disp)/2.

def checkBlockEngagement():
Pblock1 = Pblock2 = 0
for i in block1:
Pblock1+= la.norm(O.forces.f(i))
for i in block2:
Pblock2+= la.norm(O.forces.f(i))
return Pblock1, Pblock2

def stressStrainHist():
i=O.iter,
time = O.time,
disp = -O.time*dispVel,
P = getPistonForce(),
pistonVel = pistonEngine.velocity,
PistonDisp = getPistonDisp(),
Pblock1 = checkBlockEngagement(),
Pblock2 = checkBlockEngagement())
plot.saveDataTxt('txt/data'+identifier+'.txt',vars= ('i', 'disp','P', 'pistonVel','PistonDisp'))
#momentFile = 'moments_'+identifier+'.txt'

#AEcount=0
#if os.path.isfile(momentFile):
#AEcount = sum(1 for line in open(momentFile))

#f = open('txt/'+AEfile, 'a')
#P, pistonDisp = plot.data['P'],plot.data['PistonDisp']
#f.write('%g %g %g %g %g\n' % (O.time, O.iter, AEcount, P[-1], pistonDisp[-1]))
#f.close

def stopifDamaged():
P=plot.data['P']
if O.iter > 10000:
if max(P)>2000 and P[-1] < 0.5*max(P):
print('failure reached')
sigma_t = 3*max(P)*(2*(ysup-2*r))/(2*Z*X**2)
print("Tensile strength=",sigma_t/1e6)
O.pause()

plot.plots={'PistonDisp':('P'), 'i':('Pblock1','Pblock2')}
plot.plot(subPlots=True)

O.step()
### initializes the interaction detection factor
SSgeom.interactionDetectionFactor=-1.
Saabb.aabbEnlargeFactor=-1.

O.run()

waitIfBatch()
----------------------------------------------------------------------------------------------------------

I got a bit confusion about following lines,

# Rigid blocks
block1 = O.bodies.append(geom.facetCylinder(center=(xinf-0.971*r,yinf+2*r,0),radius=r,height=Z,orientation=Quaternion((0, 0, 0), 90),segmentsNumber=20,wire=False,material=wallMat))

block2 = O.bodies.append(geom.facetCylinder(center=(xinf-0.915*r,ysup-2*r,0),radius=r,height=Z,orientation=Quaternion((0, 0, 0), 90),segmentsNumber=20,wire=False,material=wallMat))

from there, I can see

xinf-0.971*r
xinf-0.915*r
xsup+0.943*r

should this be changed to

xinf-1*r
xinf-1*r
xsup+1*r

this or, I would like to know why there are different coefficients(0.971, 0.915,0.943)?

May I know how the packing is generated also?

 chanaka Udaya (chanaka-udaya) said on 2020-09-17: #3

This is the error messeage I got when I used the random dense pack algorithm

https://www.dropbox.com/s/47vh8b0fft2wh0s/error.PNG?dl=0

 Robert Caulk (rcaulk) said on 2020-09-17: #4

Hello,

Please stop linking to dropbox. I'll point you, once again after our personal communications, to our posting guidelines .

Launchpad will not let files be attached to questions and answers. As a workaround many users tend to provide links to external repositories. However removing the need for external files is part of creating a good MWE, hence external links on Yade launchpad are declared poor practice and should be avoided.

External links mean the potential responder will need to click multiple times just to read the question, while he/she is possibly reading it offline (via emails) and on an arbitrary device (e.g. cell phone). This potential responder might then simply move to next question (some of them even systematically do so).
The persistence of external links over time is not guaranteed, it will make the archive incomplete.
External links are not properly crawled into by search engines

Showing a script should be possible in launchpad, since the scripts are supposed to be short (see previous section on MWE).
If the input of one script is the output of another script, then both scripts should be merged into one. That is: avoid sending multiple scripts just like you avoid sending input files to a script. If needed, the method to save and reload in one single script is illustrated here.
Error logs (be it compilation or runtime error) can be long, sometimes, yet the relevant lines are not many. If the log is very long please identify the error which comes first and paste only that one (together with user's input command).
Input files with particle coordinates or input mesh are not needed, in most cases. They can be replaced by writing coordinates directly in the script or using makeCloud or similar functions. Most problems are not due to a specific number of particles or the specific shape of an imported geometry. Instead, they can be reproduced with just a few facets and a few spheres (or clumps or whatever is in a specific situation).
Although removing external links may need a bit more work before sending the question, this extra time spent by the person asking means less time lost by the person answering. This is just fair. Moreover, by doing so you may well understand the problem by yourself in many cases."

That being said, the example now shows how to generate the packing .

>this or, I would like to know why there are different coefficients(0.971, 0.915,0.943)?

These just accelerate the simulation by placing the blocks as close as possible to the rough surface without touching. You can place them further away and wait for them to approach the beam if you like. These values are unique to each packing.

Cheers,

Robert

 chanaka Udaya (chanaka-udaya) said on 2020-09-17: #5

Dear Dr. Rober Caulk,

Sorry for the inconvinence caused. Next I will not place external links.

I will give a go with this code and come back to you.

Thanks again for the prompt response.

Chanaka