no softening behavior of cpm in tension test

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