yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #25718
[Question #697461]: Triaxial Dense no Peak
New question #697461 on Yade:
https://answers.launchpad.net/yade/+question/697461
Hello everyone,
I am a PhD Student at INSA de lyon and I am making a rectangular triaxial simulation of a real experiment that exists here.
I have managed to create a simulation following the many examples in this forum and the triax-tutorial/session1 file. But I am having problems simulating a dense soil behavior (hardening - peak - softening). I tried the following methods to create a dense specimen:
1 - Insertion through cloud with 0 friction angle and then apply radius expansion,
2 - Insertion through cloud with a high friction angle (45) and then apply radius expansion. Next gradual reduction of the friction angle (from 45 to 0.1) with piston activation to reduce the porosity
3 - Insertion through cloud with a high friction angle (45) and then apply radius expansion. Next gradual reduction of the friction angle (from 45 to 0.1) with radius expansion to reduce the porosity
In all three cases I cannot get a peak stress, it seems like the soil is not loose but also not dense enough. In the 2 and 3 cases I manage to reduce the porosity from 0.45 to 0.4 .
In this Google drive folder you will be able to find the script and also a Image of the behavior of my simulations.
https://drive.google.com/drive/folders/1Fdw4PusTu-Li1z45kz1JrEfQeYGyga69?usp=sharing
PS: in the image the porosity is not updated after the start of the real triaxial test thus leading to a constant porosity.
PS 2 : I am a little too obssessed with automating things, so I hope my script does not confuse you.
Thank you very much for your time,
GONCALVES DE O C Joao
joao.goncalves-de-oliveira-chueire@xxxxxxxxxxxx
SCRIPT ::
# encoding: utf-8
#
# Made by GONCALVES DE OLIVEIRA CHUEIRE Joao Augusto as part of PhD at INSA de Lyon
#
# The initial packing composed of a cloud of specified granulometry (radii and percentages)
# The triaxial consists of 3 stages:
#
# 1. Calibration : after grains insertion their size is increased untill sigStar value
# is reached to create the sample asked.
# 2. Isotropic consolidation : piston advancing until sigCons is reached in all directions;
# 3. Triaxial test : Constant-strain compression along the z-axis, while maintaining constant
# stress (sigCons) laterally. Run untill a certain vertical strain level is reached;
#
# Controlling of strain and stresses is performed via TriaxialStressController,
# of which parameters determine type of control and also stability condition .
# The changes between each phase is made with the checkState function and the
# variable stt_dict.
#
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Variables #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Geometry
w=0.10 # specimen targeted widith in 0 direction
l=0.35 # specimen targeted length in 1 direction
h=0.45 # specimen targeted height in 2 direction
# Material Parameters
#grain
gDsty = 2600 # grain density
gYoung = 1e9 # grain young modulus
gPRat = 0.3 # grain poisson ratio
gFrAng = 30 # grain friciton angle
#wall
wDsty = 0 # wall density
wYoung = 50e9 # wall young modulus
wPRat = 0.3 # wall poisson ratio
wFrAng = 0.0 # wall friciton angle
# Simulation Parameters - [convention : positive is traction]
#insertion parameters
insT ='NEW' # insertion type : 'NEW' or 'LOAD'
tPo = 0.40 # targeted porosity
rfz = 0.33 # radius delta factor
nPart = 3000 # nb of particles
#triaxial parameters
sigStar =-90e3 # calibration stress
sigCons =-100e3 # consolidation stress
sigComp =-100e3 # compressison stress
thk = 0.001 # piston thickness
lRate =-0.05 # compression strain rate
dplSpd = 0.1 # piston maximal displacement speed
eMax = 0.15 # strain value to stop execution
#general
Damp = 0.2 # damp coef
# Save path - Matlab & VTK
svPath='generalvariables/'
vtkPath='vtk/'
#acceptable errors to pass to next phase
uFR = 0.03 # unbalanced force on Radius Increase phase
uFS = 0.01 # unbalanced force on Calibration phase
uFC = 0.001 # unbalanced force on Consolidation phase
stS = 5000 # stress on Calibration phase
stC = 500 # stress on Consolidation phase
mnI = 20e3 # minimal number of timesteps to nextphase
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Start files #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Auto generated variables
startT = O.time # time calculation
stt_dict={'val': 0} # state variable used on stateCheck function
# Imports
from datetime import date
import errno
from yade import pack, qt, plot,os,timing
#Create folder
try:
os.mkdir(svPath)
except OSError as exc:
if exc.errno != errno.EEXIST:
raise
pass
try:
os.mkdir(vtkPath)
except OSError as exc:
if exc.errno != errno.EEXIST:
raise
pass
#Create StressForce files
fnm=svPath+'servoforce.txt'
f=open(fnm,'w')
f.write('Yade triaxial simulation - Piston files \n')
f.write('Executed at '+date.today().strftime("%d/%m/%Y")+' by GONCALVES Joao\n')
f.write('Step\tSx\tSy\tSz\tUz\tUy1\tUy2\tUx1\tUx2 \n' )
f.close()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Materials and geometries #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# first choose a high gFrAng0
if gFrAng>45:
gFrAng0=float(gFrAng)
else:
gFrAng0=45.
#Create materials for spheres and plates
SphereMat = O.materials.append(FrictMat(young = gYoung, poisson = gPRat,
frictionAngle = radians(gFrAng0), density = gDsty))
WallMat = O.materials.append(FrictMat(young = wYoung, poisson = wPRat,
frictionAngle = radians(wFrAng), density = wDsty))
if insT=='NEW':
# Create box - box are created before spheres, thus the triax_ID is not redefined
mn,mx = Vector3(-w/2,0.,0.),Vector3(w/2,l,h)
wallIds=O.bodies.append(aabbWalls([mn,mx],thickness=thk,oversizeFactor=2,material=WallMat))
#Changing following creating type
sp=pack.SpherePack()
sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=-1,rRelFuzz=rfz,
num=nPart,periodic=False)
sp.toSimulation(material = SphereMat)
elif insT=='LOAD':
#Try to load file
loadNm=insT.upper()+'ins.yade.bz2'
try:
sp.load(loadNm)
except:
print('#**#** File not found in the root of the execution script **#**#')
exit()
else:
print('#**#** New/Load Preference not set **#**#')
exit()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Engines #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Triaxial Controller
triax1= TriaxialStressController(
# specify target values and whether they are strains or stresses
goal1=sigStar,
goal2=sigStar,
goal3=sigStar,
## 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 zyx, i.e.
## "001" (1) means "x", "011" (3) means "x and y", "111" (7) means "x and y and z", etc.
stressMask=7,
thickness = thk,
max_vel=dplSpd, # velocity
wall_back_activated = False, #turn off movment on the wall direction 2(z) movment +
internalCompaction = True, #If true the confining pressure is generated by growing particles
maxMultiplier=1.+1.5e5/gYoung, # spheres growing factor (fast growth)
finalMaxMultiplier=1.+4e3/gYoung,
computeStressStrainInterval = 10,
)
# Contact friction
newton=NewtonIntegrator(damping=Damp)
O.dt = 0.25*PWaveTimeStep()
#Stuff?
O.usesTimeStepper=True
O.timingEnabled=True
#O.trackEnergy=True # enable energy tracking in the code
# Start engine
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()],
),
#GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.4),
triax1,
newton,
PyRunner(iterPeriod=1000,command='addPlotData()'),
PyRunner(iterPeriod=1000,command='cmdUpdate()'),
PyRunner(iterPeriod=2000,command='checkState()'),
PyRunner(iterPeriod=2000,command='saveFiles()'),
#VTKRecorder(iterPeriod=5000,fileName=vtkPath+'vt',recorders=['spheres','bstresses','coordNumber','mass','velocity'])
#TriaxialStateRecorder(iterPeriod=2000,file=svPath+'WallStresses')
]
# Add plots
plot.plots={'i':('ub',),'i ':('s1','s2','s3'),' i':('e1','e2','e3'),' i ':('por'),
# energy plot' i ':(O.energy.keys,None,'totE'),
}
# Plot
plot.plot()
if insT=='LOAD':
# If load jump directly to the consolidation
prepareCalibration()
prepareConsolidation()
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Main Functions #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def checkState():
global stt_dict,uFR,uFS,uFC,stS,stC,mnI,emax,gFrAng0,Porosity
# Function that controls the triaxial engine behavior in three phases
if stt_dict['val']==0:
# Diam change : increase the diameter of the particles to create a sample
# where are particles are in contact. Using a high gFrAngle0
unb=unbalancedForce()
strXY =(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2
strZ=triax1.stress(triax1.wall_front_id)[2]
if unb<uFR and abs(sigStar-strXY)<stS and abs(sigStar-strZ)<stS and O.iter>mnI:
prepareCalibration()
if stt_dict['val']==1:
# Calibration : gradually lower the value of gFrAng0 until the target porosity is reached
unb=unbalancedForce()
if unb<uFS:
if Porosity<=tPo:
prepareConsolidation()
else:
if gFrAng0==0.20:
prepareConsolidation()
print('Targeted porosity could not be reached, continuing with p=' + str(Porosity))
return
elif gFrAng0>5.:
gFrAng0=5.
elif gFrAng0==5.:
gFrAng0=1.
elif gFrAng0==1.:
gFrAng0=0.20
O.materials[0].frictionAngle=radians(gFrAng0)
print('Ang='+str(gFrAng0)+' Porosity :'+ str(Porosity))
setContactFriction(radians(gFrAng0))
stt_dict.update({'consoSt':O.iter})
elif stt_dict['val']==2:
# Consolidation : insotropic control untill all pistons are applying the same
# stress to the grains.
unb=unbalancedForce()
strXY =(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2
strZ=triax1.stress(triax1.wall_front_id)[2]
if unb<uFC and abs(sigCons-strXY)<stC and abs(sigCons-strZ)<stC and (O.iter-stt_dict['consoSt'])>mnI:
prepareCompression()
else:
# Compression : piston direction 2 control changed to strian. Run untill the
# maximal strain defined above is applyed.
e2=-triax1.strain[2]
if e2>eMax:
O.save(svPath+insT.upper()+'compEnd.yade.bz2')
stt_dict.update({'compEnd':O.iter})
print('#**#** End Compression at '+ str(stt_dict['compEnd'])+' **#**#')
# Prepare to end simulation
#Save plot data
plot.saveDataTxt(svPath+'PlotData.txt.bz2')
#Save last files
saveFiles()
#create file for Matlab
endFile()
O.pause()
def prepareCalibration():
global stt_dict
# Count nb of spheres
nbSph=len(O.bodies)-6
# Save variables for later usage
stt_dict.update({'val': 1,'consoSt':O.iter,'nbS':nbSph})
# add gravity
newton.gravity = (0.0,0.0,-9.81)
# Prind in cmd line
print('#**#** End Diametre change at '+ str(O.iter)+' **#**#')
def prepareConsolidation():
global stt_dict
# Prepare for consolidation
# Save state
O.save(svPath+insT.upper()+'ins.yade.bz2')
# Save variables for later usage
stt_dict.update({'val': 2,'consoSt':O.iter,'Pz':O.bodies[5].state.pos[2],
'Py1':O.bodies[3].state.pos[1],'Py2': O.bodies[2].state.pos[1],
'Px1':O.bodies[1].state.pos[0],'Px2':O.bodies[0].state.pos[0]})
# Change piston states
triax1.goal1 = sigCons
triax1.goal2 = sigCons
triax1.goal3 = sigCons
triax1.depth0 = triax1.depth
triax1.height0 = triax1.height
triax1.width0 = triax1.width
triax1.internalCompaction = False
#update material values
O.materials[0].frictionAngle=radians(gFrAng)
setContactFriction(radians(gFrAng))
# restart plot data
#if insT!='LOAD':
#plot.saveDataTxt(svPath+'InsertionPlotData.txt.bz2')
#plot.resetData()
# Prind in cmd line
print('#**#** End Calibration at '+ str(stt_dict['consoSt'])+' **#**#')
print('w: '+str(triax1.depth)+'\t l: '+str(triax1.height)+'\t h: '+ str(triax1.width))
def prepareCompression():
global stt_dict
# Save variables
O.save(svPath+insT.upper()+'consoEnd.yade.bz2')
stt_dict.update({'val': 3,'consoEnd':O.iter})
# Change piston states
triax1.goal1 = sigComp
triax1.goal2 = sigComp
triax1.goal3 = lRate
triax1.stressMask=3 # integer for binary 011
# Prind in cmd line
print('#**#** End Consolidation at '+ str(stt_dict['consoEnd'])+' **#**#')
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Execution Data #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def cmdUpdate():
sxy=(-triax1.stress(triax1.wall_right_id)[0]-triax1.stress(triax1.wall_top_id)[1])/2
print('I: %d\t ubf: %.4f\tsxy: %.0f\tsz: %.0f\tez: %.3f'%(O.iter,unbalancedForce(),
sxy,-triax1.stress(triax1.wall_front_id)[2],-triax1.strain[2]))
# Plots preparation
def addPlotData():
global Porosity
#mStress = -(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1]+
# triax1.stress(triax1.wall_front_id)[2])/3.0
#Vol = ((O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*
# (O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])*
# (O.bodies[5].state.pos[2]-O.bodies[4].state.pos[2]))
# Volume spheres
Vsph = sum(4.0/3.0*pi*b.shape.radius**3 for b in O.bodies if isinstance(b.shape,Sphere))
Porosity = 1-Vsph/(h*w*l)
plot.addData(e1=-triax1.strain[0], e2=-triax1.strain[1], e3=-triax1.strain[2],
ev=-triax1.strain[0]-triax1.strain[1]-triax1.strain[2],
s1=-triax1.stress(triax1.wall_right_id)[0],
s2=-triax1.stress(triax1.wall_top_id)[1],
s3=-triax1.stress(triax1.wall_front_id)[2],
ub=unbalancedForce(),
i=O.iter,
por=Porosity
#totE=O.energy.total(),**O.energy
#dstress=(-triax1.stress(triax1.wall_top_id)[1])-(-triax1.stress(triax1.wall_right_id)[0]),
#disx=triax1.width,
#disy=triax1.height,
#disz=triax1.depth,
#porosity=Porosity,
#ke=utils.kineticEnergy(),
)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
############################# Save Functions #############################
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
def saveGrains():
# Save grains details into a txt file for matlab to read
global stt_dict
# per grain stress calculation following https://yade-dem.org/doc/user.html#micro-stress-and-micro-strain
TW=TesselationWrapper()
TW.setState()
TW.computeVolumes()
stress=bodyStressTensors()
# start writing
fnm=svPath+'grains'+str(O.iter)+'.txt'
f=open(fnm,'w')
f.write('Yade triaxial simulation - Grains files for step'+str(O.iter)+'\n')
f.write('Executed at '+date.today().strftime("%d/%m/%Y")+' by GONCALVES Joao\n')
f.write('Number of grains : '+str(stt_dict['nbS'])+'\n')
f.write('ID\tX\tY\tZ\tRadius\tUx\tUy\tUz\tSx\tSy\tSz\tSxy\tSxz\tSyz \n' )
for k in O.bodies:
if isinstance(k.shape, Sphere):
stp = stress[k.id]*4.*pi/3.*k.shape.radius**3/TW.volume(k.id)
[ID,r,pos,vel]=[k.id,k.shape.radius,k.state.pos,k.state.vel]
f.write(str(ID)+'\t'+str(pos[0])+'\t'+str(pos[1])+'\t'+str(pos[2])+
'\t'+str(r)+'\t'+str(vel[0])+'\t'+str(vel[1])+'\t'+str(vel[2])+
'\t'+str(stp[0,0])+'\t'+str(stp[1,1])+'\t'+str(stp[2,2])+
'\t'+str((stp[0,1]+stp[1,0])/2)+'\t'+str((stp[0,2]+stp[2,0])/2)+
'\t'+str((stp[1,2]+stp[2,1])/2)+'\n'
)
f.close()
def saveContacts():
# Save contacts details into a txt file for matlab to read
global stt_dict
fnm=svPath+'contForce'+str(O.iter)+'.txt'
f=open(fnm,'w')
f.write('Yade triaxial simulation - Contact files step'+str(O.iter)+'\n')
f.write('Executed at '+date.today().strftime("%d/%m/%Y")+' by GONCALVES Joao\n')
f.write('Number of grains : '+str(stt_dict['nbS'])+'\n')
f.write('ID1\tID2\tFx\tFy\tFz \n' )
for k in O.interactions:
[ID1,ID2,force]=[k.id1,k.id2,k.phys.normalForce+k.phys.shearForce]
f.write(str(ID1)+'\t'+str(ID2)+'\t'+str(force[0])+'\t'+str(force[1])+'\t'+str(force[2])+'\n')
f.close()
def savePiston():
# Save pistons details into a txt file for matlab to read
fnm=svPath+'servoforce.txt'
f=open(fnm,'a')
f.write(str(O.iter)+'\t'+str(triax1.stress(triax1.wall_right_id)[0])+'\t'+str(triax1.stress(triax1.wall_top_id)[1])+'\t'+
str(triax1.stress(triax1.wall_front_id)[2])+'\t'+str(O.bodies[5].state.pos[2])+'\t'+
str(O.bodies[3].state.pos[1])+'\t'+str(O.bodies[2].state.pos[1])+'\t'+
str(O.bodies[1].state.pos[0])+'\t'+str(O.bodies[0].state.pos[0])+'\n'
)
f.close()
def saveFiles():
# Command to lauch all three save files in one pyrunner.
# Only launch after the calibration part
if stt_dict['val']>1:
saveGrains()
saveContacts()
savePiston()
def endFile():
# Print file with important information for matlab
global stt_dict
print('Time taken:',(O.time - startT))
#YADEtoMatlab file
f=open('YADEtoMatlab.txt','w')
f.write('The following values give the Matlab app information about the simulaiton\n')
f.write('Executed at '+date.today().strftime("%d/%m/%Y")+' by GONCALVES Joao\n')
f.write('Number of grains : '+str(stt_dict['nbS'])+'\n')
f.write(' \n')
f.write('ExpWidth = '+str(triax1.width0)+'\n' )
f.write('ExpLength = '+str(triax1.height0)+'\n' )
f.write('ExpHeight = '+str(triax1.depth0)+'\n' )
f.write('Timestep = '+str(O.dt)+'\n' )
f.write('StartCons = '+str(stt_dict['consoSt'])+'\n' )
f.write('EndCons = '+str(stt_dict['consoEnd'])+'\n' )
f.write('EndComp = '+str(stt_dict['compEnd'])+'\n' )
f.write('vtkFolder = '+svPath[:-1]+'\n' )
f.write('genFolder = '+svPath[:-1]+'\n' )
f.write('PositionPz = '+str(stt_dict['Pz'])+'\n' )
f.write('PositionPy1 = '+str(stt_dict['Py1'])+'\n' )
f.write('PositionPy2 = '+str(stt_dict['Py2'])+'\n' )
f.write('PositionPx1 = '+str(stt_dict['Px1'])+'\n' )
f.write('PositionPx2 = '+str(stt_dict['Px2'])+'\n' )
f.write(' \n')
f.close()
############################# Auto Run #############################
#O.run(); O.wait()
#For info :
# %F.Pt => F stand for flag (nb of characters), P stand for precision (nb of val after .) and t for the
#type of data (d signed integer decimal, floating point decimal, e floating point exponential format)
#f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16d %-16g\n'%(pos[0],pos[1],pos[2],rad,displ[0],displ[1],displ[2],b.id,b.material.density))
--
You received this question notification because your team yade-users is
an answer contact for Yade.
