New question #695734 on Yade:
https://answers.launchpad.net/yade/+question/695734
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()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
