calibration of CPM model

Asked by Othman Sh on 2019-11-13

Hello,

I'm modeling the elastic stress-strain behavior of porous concrete cylinder under compression. I have experimental results of a porous concrete with a porosity = 0.25 and the slope of the elastic stress-strain curve (elastic modulus) is about 2e6 psi and the strength (max stress) is 1600 psi (the curve is linear elastic). I'm trying to calibrate my DEM model (single size spheres, utils.porosity()=0.25, sphere contact model is CPM) to get similar E (elastic modulus of the whole packing) and strength as experimental results but I couldn't (I only want to match the elastic portion of the stress-strain curve). If I get to the 1600 strength level, the E would be way too high (~1e9 psi). When I play with the CPM parameters I could get the E close to 2e6 but the strength will be way too low (less than 200 psi). The CPM parameters that I'm playing with are:

elastic modulus of the contact
Poisson's ratio
frictionAngle
epsCrackOnset
relDuctility
sigmaT

Note that every run I use ymport and import exactly same sphere packing. I tried using packing with spheres with different sizes but the results were not much different. Could you please provide some recommendations to have my DEM model reach similar strength and E of the experimental results? Is there a cohesive contact model in Yade that can fit my case better that CPM?

My code is shown below. it consists of two codes, one for packing the spheres and the other one is for compression test simulation.

Thanks
Othman

----------------------------------------
Compression test code:

from yade import plot,pack, export, ymport
from pprint import pprint

#specimen geometry
radiuscyl=.05

################ Material model parameters ################

intRadius=1.5
pervconc=O.materials.append(CpmMat(

 young = 3.6e6,
 poisson = 0.1,
 frictionAngle = atan(1.2),
 epsCrackOnset = 3e-2,
 relDuctility = 60,
 sigmaT = 1e6,
))

################ spheres ################

spheres=O.bodies.append(ymport.text('/home/RMC1_B25.txt')) #please change the file location in order to import the packing from the packing code provided below
yade.qt.View()
cylheight=max(utils.aabbDim())*39.37 #dimensions in inches for plotting
print ('cylheight',cylheight)
################ creating disks ################

#Wall
zMin,zMax=[pt[2] for pt in aabbExtrema()]
wallIDs = O.bodies.append([wall((0,0,z),2) for z in (zMin,zMax)])
walls = wallMin,wallMax = [O.bodies[i] for i in wallIDs]
#Applying force by moving the walls in constant displacement (m/second)
wallMin.state.vel = (0,0,0.000001)
wallMax.state.vel = (0,0,-0.005)

################ engines ################
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='is2aabb'),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intRadius,label='ss2sc'),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_CpmMat_CpmMat_CpmPhys()],
  [Law2_ScGeom_CpmPhys_Cpm(epsSoft=1e-3)],
 ),
 NewtonIntegrator(damping=.2),
 CpmStateUpdater(realPeriod=.5),
 PyRunner(iterPeriod=100,command='addPlotData()',initRun=True),
 PyRunner(iterPeriod=100,command='stopIfDamaged()'),
]

################ stop condition ################
def stopIfDamaged():
 if O.iter < 1000: # do nothing at the beginning
  return
 fMax = max(plot.data["f"])
 f = plot.data["f"][-1]
 if f/fMax < .01:
  print "Damaged, stopping."
  print "f'c = ",max(plot.data["stress"])
  print "max force = ", max (plot.data ["f"])
  O.pause()
  plot.saveDataTxt('R1-3-5.txt')

################ plot stuff ################
# see how to fix this: plot.saveDataTxt('bbb.txt.bz2')
def addPlotData():
 # forces of walls. f1 is "down", f2 is "up" (f1 needs to be negated for evlauation)
 f1,f2 = [O.forces.f(i)[2] for i in wallIDs]
 f1 *= -1
 # average force (in lbf)
 f = .5*(f1+f2)*.224809
 # displacement (2 times each wall) (inches)
 wall = O.bodies[wallIDs[0]]
 displacement = 39.37*2*wall.state.displ()[2]
 # stress (psi)
 stress = f/(1550*pi*radiuscyl**2)
 # store values
 yade.plot.addData(
  t = O.time,
  i = O.iter,
  displacement = displacement,
  f1 = f1,
  f2 = f2,
  f = f,
  strain=displacement/cylheight,
  stress = stress,

 )

plot.plots={('strain'):('stress',)}

O.dt = 0.
O.step(); # to create initial contacts
# now reset the interaction radius and go ahead
ss2sc.interactionDetectionFactor=1
is2aabb.aabbEnlargeFactor=1

O.dt = .8*PWaveTimeStep()

plot.plot()

O.run()

----------------------------------------
sphere packing code:

from yade import pack, export, ymport

targetp = .25 ##specify the targeted porosity

##specimen geometry
radiuscyl=.05
heightcyl=0.203*4.32
SphereRadius = .005205
# material parameters
O.materials.append(FrictMat(young = 5e10, poisson = 0.15,frictionAngle = atan(.2), density=1920))

############################ spheres #############################
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(.3,.3,1.5),rMean=SphereRadius,rRelFuzz=0,periodic=True,seed=1) ##
#### cylinder extraction
pred=pack.inCylinder((.2,.2,0),(.2,.2,heightcyl),radiuscyl)
spFilter=filterSpherePack(pred,sp,Material=Material, returnSpherePack=True)
spFilter.toSimulation()
############################ facets #############################

facets=geom.facetCylinder((.2,.2,heightcyl/2),radiuscyl,heightcyl,segmentsNumber=150,wallMask=4)

cylinder=O.bodies.append(facets)
yade.qt.View()

##creating disks

d1=geom.facetCylinder((.2,.2,heightcyl),radiuscyl-.0001,0,segmentsNumber=150,wallMask=1)
d2=geom.facetCylinder((.2,.2,0),radiuscyl-.0001,0,segmentsNumber=150,wallMask=1)

disk1IDs= O.bodies.append(d1)
disk2IDs= O.bodies.append(d2)

for i in disk1IDs:
 body= O.bodies[i]
 body.state.vel = (0,0,-20)

for n in disk2IDs:
 body= O.bodies[n]
 body.state.vel = (0,0,0)

############################ compaction #############################
O.dt=.5*utils.PWaveTimeStep()
enlargeFactor=1.5
O.engines=[
 ForceResetter(),
 InsertionSortCollider([
  Bo1_Sphere_Aabb(),
  Bo1_Facet_Aabb()
 ]),
 InteractionLoop(
  [
   Ig2_Sphere_Sphere_ScGeom(),
   Ig2_Facet_Sphere_ScGeom(),
  ],
  [

   Ip2_FrictMat_FrictMat_FrictPhys(),
   Ip2_FrictMat_FrictMat_FrictPhys(),
  ],
  [

   Law2_ScGeom_FrictPhys_CundallStrack(),
  ],
 ),

 NewtonIntegrator(damping=.3),

 PyRunner(iterPeriod=50,command='stop()'),

]
O.run()

def stop():
   if utils.porosity()<targetp:
 O.pause()
 print 'Finished'
 for i in disk1IDs: O.bodies.erase(i)
 for i in disk2IDs: O.bodies.erase(i)
 for i in cylinder: O.bodies.erase(i)
 print ('unbalanced forces = ', utils.unbalancedForce())
 print ('cylinder dimensions = ', utils.aabbDim()) # this shows the dimensions of the cylinder (Diameter in x-axis, Diameter in y-axis, height in z-axis)
 print ('porosity = ', utils.porosity())
 print ('No. of spheres = ', len (O.bodies))
 O.wait()
 export.textExt('RMC1_B25.txt','x_y_z_r')

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Othman Sh
Solved:
2019-11-14
Last query:
2019-11-14
Last reply:
2019-11-13
Jan Stránský (honzik) said : #1

Hello,

elastic behavior is controlled by young and poisson parameters. These values should be calibrated first and left unchanged for the inelastic part of the calibration.

The inelastic behavior is controlled by the rest of parameters you mentioned.
If you do not find good strength for no parameter combination, maybe the material model is not suitable for your case.

There might be several reasons, e.g. the strength is correct, but after failure, the pores are collapsed and new particles comes to contact "shadowing" the strength..

Could you post the final compaction file (I have tried your code, but after relatively long time without finish I have quit it)?

cheers
Jan

Othman Sh (othman-sh) said : #2

Hi Jan,

Thanks for your reply. Please find the packing file in the link below.

https://drive.google.com/open?id=1yriqxJrPVWz6Vx58p31zVKYLRwEeZY-M

Thanks,
Othman

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

Hi,

thanks for the data.
For such case, fitting the strength should not be a problem IMO..

> O.dt = .8*PWaveTimeStep()

this is very sharp dt, consider using a lower value (as your particles overlap a lot, having the effect that their "size" from interaction point of view is actually smaller).
For a computed strength, try to re-simulate it with half and/or quarter dt value to verify if it is correct. If the value would differ significantly, then you should use smaller dt.

just a note, epsSoft should be negative if you want it to have effect [1]

cheers
Jan

[1] https://gitlab.com/yade-dev/trunk/blob/master/pkg/dem/ConcretePM.cpp#L400

Othman Sh (othman-sh) said : #4

Thanks Jan,

I try what you recommended and see what happen.

Regards,
Othman