New question #701369 on Yade:
https://answers.launchpad.net/yade/+question/701369
i am trying to stimulate the small strain level, i set the loadrate,
timestep,and max_vel very small. but i cannot get a constant small strain
stiffness, the curve of stress-strain always has fluctuation. i cannot figure
it out . here is my code
FrictionAngle = 35
Damp = 0.25
#Confining variables
ConfPress = -2.0e5
#Loading control
LoadRate = -0.001
O.load('calm.yade.bz2')
########################################
#import necessary packages
from yade import pack,plot,os,timing
triax=O.engines[4]
#################################
#####Defining triaxil engines####
#################################
triax=TriaxialStressController(
thickness = 0.001,
internalCompaction = False, # If true the confining pressure is generated
by growing particles
max_vel = 1e-20,
stressMask = 5,
goal1 = ConfPress,
goal2 = LoadRate,
goal3 = ConfPress,
)
ini_e22a = -triax.strain[1]
ini_e22b = -triax.strain[1]
#O.usesTimeStepper=True
O.dt=1e-6
newton=NewtonIntegrator(damping=Damp)
setContactFriction(radians(FrictionAngle))
O.trackEnergy=True
O.timingEnabled=True
###engine
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys()],
[Law2_ScGeom_MindlinPhys_Mindlin()]
),
#GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
newton,
PyRunner(iterPeriod=1,command='endCheck()',label='check'),
PyRunner(command='addPlotData()',iterPeriod=10,label='record'),
#VTKRecorder(iterPeriod=100,recorders=['all'],fileName='./post/p1-'),
TriaxialStateRecorder(iterPeriod=100,file='WallStresses3')
]
# Simulation stop conditions defination
def endCheck():
unb=unbalancedForce()
e22=-triax.strain[1]
if abs(e22)>1.5e-5:
O.pause()
O.save('e0.05.yade.bz2')
# collect history of data
def addPlotData():
global ini_e22a
global ini_e22b
e22=-triax.strain[1]
unb = unbalancedForce()
mStress =
-(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3.0
volume =
(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])
Porosity=1-sum(4*pi*b.shape.radius**3/3 for b in O.bodies if
isinstance(b.shape,Sphere))/volume
plot.addData(e11=-triax.strain[0], e22=-triax.strain[1],
e33=-triax.strain[2],
ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
s11=-triax.stress(triax.wall_right_id)[0],
s22=-triax.stress(triax.wall_top_id)[1],
s33=-triax.stress(triax.wall_front_id)[2],
ub=unbalancedForce(),
dstress=(-triax.stress(triax.wall_top_id)[1])-(-triax.stress(triax.wall_right_id)[0]),
disx=triax.width,
disy=triax.height,
disz=triax.depth,
i=O.iter,
porosity=Porosity,
cnumber=avgNumInteractions(),
ke=utils.kineticEnergy(),
totale=O.energy.total()
)
print ('unbalanced force: %f, mean stress: %f, s11: %f, s22: %f, s33: %f,
coordination number: %f, porosity: %f'
%(unb,mStress,-triax.stress(triax.wall_right_id)[0],-triax.stress(triax.wall_top_id)[1],-triax.stress(triax.wall_front_id)[2],avgNumInteractions(),Porosity))
plot.saveDataTxt('loadinglog3.txt.bz2')
##################################################
######Defining output for postprocessing##########
##################################################
plot.plots={'e22':('s22',None,'ev')}
plot.plot()
--
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
