no softening behavior of cpm in tension test

Asked by Zhoufeng Shi on 2021-02-25

hello everyone,
I am trying to calibrate the concrete by using cpm. But I find that no matter how I adjust the parameters (eg, epsCrackOnset, relductility), the uniaxial tension test in model always have a linear elastic state and then followed by a sharp drop in stress-strain diagram, the tangent of the stress decrease is almost 90 degree. But in experiment data, there should be a softening part after the peak stress. The model seems so too 'brittle'. Is it because of the interaction enlarge factor? I took 1.5 as recommended.
I try to reduce the factor to 1.1, 1.05,etc and the model start to have softening behavior but the elastic modulus has a huge decrease from 180 to 50 MPa. and I have to increase the young's to a very large value in order to match the elastic modulus, but then I found excessive long hardening curve before the peak, which should be there. Any advice or thoughts is welcomed! (the code attach below may take 2-3hr to finish at current strainrate). Thanks!!!
#########################code begin###########################

###############get dense pack of particles########################################
from yade import pack,plot
# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
 num_spheres=10000,# number of spheres
 compFricDegree = 30, # contact friction during the confining phase
 key='_triax_base_', # put you simulation's name here
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.55 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=-0.02 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=5e6 # contact stiffness
mn,mx=Vector3(0,-.0075,0),Vector3(0.7,.0075,.15) # corners of the initial packing (0,-.0075,0),(0.7,.0075,.15)

O.switchScene();O.resetThisScene()
## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=aabbWalls([(0,-.1,0),(.05,.1,.3)],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(.05,0,.3),-1,0.3,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
spp=O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
for s in spp:
  O.bodies[s].state.blockDOFs='yXZ'
triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
 ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
 ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
 ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
 ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
 stressMask = 7,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [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)
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton
]
#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal3=-500
triax.goal2=0
while 1:
  O.run(1000, True)
  ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  print ('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<stabilityThreshold and abs(-1000/3-triax.meanStress)/1000/3<0.001:
    break
sp=SpherePack(); sp.fromSimulation()

print ("### Isotropic state saved ###")
O.switchScene()
##define parameters
readParamsFromTable(strainRateTension=.1,sphereRadius=.0005,young=1.5e9,sigmaT=4e6,frictionAngle=atan(0.6),relDuctility= 1.1,omegaThreshold=0.99,epsCrackOnset=.025,damping=.6,sampleNum=1,damLaw = 1,relFuzz=0,factor=1.5,strength_law=1,std=0.4,low=0.5,high=1.5)
from yade.params.table import *
import numpy as np
from yade import pack,plot
isoPrestress=0
poisson= 0.1
pred=pack.inAlignedBox((-1,-1,-1),(1,1,1))-pack.inAlignedBox((-.005,-.0075,-.075),(.005,.0075,-.033))
idConcrete=O.materials.append(CpmMat(young=young,poisson=poisson,frictionAngle=frictionAngle,density=3120,sigmaT=sigmaT,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress,relDuctility = relDuctility,damLaw = damLaw))
cement=sp.toSimulation(material=idConcrete)
################################################################### start the simulation################################
###obtain the information of the model
volume = 0
NumOfParticle = 0
for i in cement:
  volume += pi * (O.bodies[i].shape.radius) ** 2
  NumOfParticle += 1
print(f'density = {volume/(.05*.3)}')
print(f'number of particle = {NumOfParticle}')
bb=uniaxialTestFeatures()
negIds,posIds,axis,crossSectionArea=bb['negIds'],bb['posIds'],bb['axis'],bb['area']
O.engines=[
  ForceResetter(),
  InsertionSortCollider([
    Bo1_Sphere_Aabb(aabbEnlargeFactor=factor,label='bo1s'),],verletDist=.05*sphereRadius), ## expand the collision detection to creat initial interactions
      ##
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=factor,label='ig2ss'), ## create collision geometry
      Ig2_Facet_Sphere_ScGeom()],
    [Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=1)], ## creat collision phys
    [Law2_ScGeom_CpmPhys_Cpm(omegaThreshold=omegaThreshold)],
   ),
  GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3),
  NewtonIntegrator(damping=damping,label='damper'),
  CpmStateUpdater(realPeriod=.5), ## update CpmState of bodies based on interactions and can change the bodies's color depending on the damage
  UniaxialStrainer(strainRate=strainRateTension,axis=axis,asymmetry=0,posIds=posIds,negIds=negIds,crossSectionArea=crossSectionArea,blockDisplacements=False,blockRotations=False,setSpeeds=True,label='strainer'),## apply strain on the both end of the body along the axis and the strain speed is set at the beginning directly
  #PyRunner(iterPeriod=100,command='crackCount()'),
  PyRunner(iterPeriod=100,command='addPlotData()'),
  PyRunner(realPeriod=5,command='stopIfDamaged()',label='damageChecker'), ## set the stop criteria for the simulation
  #qt.SnapshotEngine(fileBase='3d',iterPeriod=1000,label='snapshot'),
  #PyRunner(command='finish()',iterPeriod=20000
]
#sigma=strainer.avgStress+isoPrestress
O.step() ## go one step to creat interactions
bo1s.aabbEnlargeFactor=1 ## reset the interaction radius
ig2ss.interactionDetectionFactor=1

### reinforcing the boundary
dim=utils.aabbExtrema()
if strainRateTension>0:
  layerSize=1/6
  height=dim[1][axis]-dim[0][axis]
  for b in O.bodies:
    if isinstance(b.shape,Sphere):
      if (b.state.pos[axis]<(dim[0][axis]+layerSize*height)) or (b.state.pos[axis]> (dim[1][axis]-layerSize*height)):
        b.shape.color=(1,1,1)
  for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
      if O.bodies[i.id1].shape.color==(1,1,1) or O.bodies[i.id2].shape.color==(1,1,1):
        i.phys.neverDamage=True

def stopIfDamaged():
  if O.iter<2 or 'sigma' not in plot.data : return
  epsilon = plot.data['eps']
  sigma=plot.data['sigma']
  extremum=max(sigma)
  ##Gf calculator
  global Gf
  if abs(sigma[-1]/extremum)>0.2:
    Gf += 0.5*(sigma[-1]+sigma[-2])*(epsilon[-1]-epsilon[-2])
  ##elastic modulus calculator
  up=round(sigma.index(extremum)*0.5)
  down=round(sigma.index(extremum)*0.125)
  elasticModulus = (sigma[up]-sigma[down])/(epsilon[up]-epsilon[down])
  if abs(sigma[-1]/extremum)<0.2: ## stop the test if the sigma is too small or strain is too large
    O.pause()

def addPlotData():
  plot.addData(i=O.iter,sigma=strainer.avgStress,eps=strainer.strain)
plot.plots={'eps':('sigma')}
plot.plot()
O.saveTmp()
O.run()
waitIfBatch()

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2021-03-06
Last reply:
2021-03-06
Jan Stránský (honzik) said : #1

Hello,

> no matter how I adjust the parameters (eg, epsCrackOnset, relductility)

please provide the specific values

> The model seems so too 'brittle'. Is it because of the interaction enlarge factor?

no, it is because of parameters / combination of parameters. interaction enlarge factor is just one of them, not so important itself (IMO).

> I try to reduce the factor to 1.1, 1.05,etc and the model start to have softening behavior but the elastic modulus has a huge decrease from 180 to 50 MPa. and I have to increase the young's to a very large value in order to match the elastic modulus ...

'young' itself has nothing like "very large value". If the resulting macroscopic elastic modulus is OK, than the 'young' has just the right value.
Of course, changing young influences almost anything else and triggers need of recalibration of other parameters.

> but then I found excessive long hardening curve before the peak

hardening might be because "diagonal" interactions are already in nonlinear mode (plastic sliding), but damage is not yet present.
If you are interested, try different decreased values of sigmaT.

> readParamsFromTable(...,relDuctility= 1.1,...,epsCrackOnset=.025,...)

try relDuctility 2, 5, 10, 100, 1000, ...
relDuctility=1.1 should result in very brittle material :-)
0.025 seems relatively large (but depends on your needs)

Concerning calibration, have a look at [1], section 3.3

cheers
Jan

[1] https://www.researchgate.net/publication/45194493_Cohesive_Particle_Model_using_the_Discrete_Element_Method_on_the_Yade_Platform

Zhoufeng Shi (zhoufeng-s) said : #2

Hi, Jane
Thanks for your prompt response!

>please provide the specific values

Actually I've been stuck in this for nearly a month. I've tried quite a lot combinations like epsCrackOnset = 0.02, relductility= 1.05, 1.5, 2, 10, 50, 100, 300, 500. but from my simulation results, it seems that by changing relductility only affect the hardening process and after peak it always goes straight down. (Do you know how to post picture here? So that I may share some of my simulation results)

>hardening might be because "diagonal" interactions are already in nonlinear mode (plastic sliding), but damage is not yet present.
If you are interested, try different decreased values of sigmaT.

I found sigmaT didn't have much influence in the tension test so I just remain it untouched. I will try different values tonight to see how it affects the results.

>Concerning calibration, have a look at [1], section 3.3

Thanks for the paper, honestly I was following this paper to do the calibration process. The calibration of the elastic part went smoothly but when it came to the characteristic length KnGf/ft2 which is the step to determine the relductility, the 'too brittle' problem prevent me from moving on.

PS: my stress-strain diagram is very similar to the [2] fig6.
Regards,
Zhoufeng

[2].Scholtès, L., & Donzé, F. V. (2013). A DEM model for soft and hard rocks: role of grain interlocking on strength. Journal of the Mechanics and Physics of Solids, 61(2), 352-369.

Zhoufeng Shi (zhoufeng-s) said : #3

hi, Jan
sorry for type your name wrong last post:(
I tried changing sigmaT value and it indeed affect the hardening a lot. However, the behavior after peak stress is still showing brittle as long as I keep the interaction enlarge factor as 1.5. Also, I notice that there is nearly no crack until it reaches peak stress and all suddenly appears at that moment, could this phenomenon relate to my problem?
Cheers,
Zhoufeng

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

> the behavior after peak stress is still showing brittle as long as I keep the interaction enlarge factor as 1.5.

interaction enlarge factor should not play so important role, maybe opposite - that value close to 1 could be problematic.
What value of relDuctility do you use for these results?

> I notice that there is nearly no crack until it reaches peak stress and all suddenly appears at that moment, could this phenomenon relate to my problem?

yes :-)
the described behavior fits together with the described stress-strain diagram

cheers
Jan

Zhoufeng Shi (zhoufeng-s) said : #5

Hi Jan,
Yes, you are right, the interaction factor didn't make much difference on the overall behavior of the model. Previously my model was packed by randomDensePack function and it seemed a bit loose.
I try to manipulate the bond strength like the code I attached below so that some of the bonds may broke at early stage. And the cracks do appear earlier but still the shape of stress-strain curve looks the same except the peak stress & strain changes, and the relDuctivity seems only affect the hardening process.
for i in O.interactions:
   if strength_law == 0:
     alpha=np.random.uniform(low,high)
   elif strength_law == 1:
     alpha=np.random.normal(1,std)
   if alpha<0.01:
     alpha = 0.01
   i.phys.epsCrackOnset *= alpha
   i.phys.undamagedCohesion *= alpha

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

> the relDuctivity seems only affect the hardening process.

in the model design, relDuctility should affect the softening behavior :-)
If it only affects the hardening, then there is some problem of normal/shear parameters interplay, and IMO only detailed analysis (many interactions in many stages) can reveal the reason / solution...

cheers
Jan

Zhoufeng Shi (zhoufeng-s) said : #7

Hi, Jan
Thanks for the suggestions. I've been try to use the same way to establish other models such as 3-point bending, compression scenario and they all look good to the real experiment data except this uniaxial tension test. I guess I'll start with two particles and gradually increase the number to see at which point the softening behavior disappear.
regards,
Zhoufeng

Zhoufeng Shi (zhoufeng-s) said : #8

Hi, Jane
problem is solved. it seems that 0.9 is too high a value for omegaThreshold especially when the relductivity is large (eg.over 10). so when the value is increased to 0.99999 or something like that, there will be softening curve for the stress-strain diagram. But I am also wondering if there's way to control the omega threshold so that I can define different omegaThreshold for each interaction, for example, a uniform distribution between (0.5,1)?
Cheers,
Zhoufeng

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

> it seems that 0.9 is too high a value for omegaThreshold

sorry, haven't noticed this. For sure, 0.9 is too little (not too high) value. Use 1.0 (default) for these kind of tests, then decrease the value if needed. Because of the exponential character of the softening law and significant non-linearity between damage value [0-1] and relative strength value [0-1], the values like 0.999999 could be OK. 0.9 is really too little.

cheers
Jan

Zhoufeng Shi (zhoufeng-s) said : #10

thanks for the response.
but it seems that crack (bond broken) only happened after peak stress which is abnormal because both in experiment and PFC simulation, the crack should start earlier than that and slightly after the peak capacity, a surge in number of cracks is expected. But in my case, the crack number is 0 until after the peak and then gradually increase to a platform.

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

> crack (bond broken) only happened
> crack should start earlier
> the crack number is 0
> ...

depends very much on the definition of "crack".
For omegaThreshold=1, there is no real crack, because no interaction is deleted, but some of them has already no (or negligible) strength. You can try "one interaction" tension test for some values of omegaThreshold to see the effect. For 0.9, the interaction is (I think) deleted shortly after the peak, resulting in very brittle behavior.
On the other hand, there should be a zone with very localized strain, which can be represented as a crack (even no interation is really cracked).

What you can try:
- play with model parameters
- play with parameter randomization
- use another material model
- ...

> in ... PFC simulation

Do you know the material model used? How is it different from Yade CPM? Maybe it is a better choice for this scenario.

There is also a possibility that CPM is simply not suitable for your scenario.

cheers
Jan

Zhoufeng Shi (zhoufeng-s) said : #12

hi,Jan
Thanks for your suggestions. The model in PFC is parallel bond model and I think the difference is that in parallel bond model its moment is contributing to the bond broken but it doesnt have any plastic softening behavior in contact-level. I'll try randomize some parameters but I don't know how to randomize omegaThreshold because I cant find the attribute for that in dir(Interactions). Could you give me some hints for that?
Regards,
Zhoufeng

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

omegaThreshold is currently a Law2 attribute, i.e. one global value..

Can you help with this problem?

Provide an answer of your own, or ask Zhoufeng Shi for more information if necessary.

To post a message you must log in.