calibration of CPM model

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