New question #703247 on Yade:
https://answers.launchpad.net/yade/+question/703247
I want to do a series of biaxial unloading tests.
Before the unloading, I do a biaxial loading test, and save the model in many
".bz2" files when the vertical strain reaches different values. The obtained
files are: save0.00010yade.bz2, save0.00020yade.bz2, save0.00030yade.bz2...
Now, I want to unload the sample from these values of vertical strain, until
the vertical stress reduces to the confining pressure. Because there are too
many .bz2 files, I hope the code could call these files one by one, and after
one unloading test is completed, another test would be done automatically.
I first write a code that can do an unloading test once at a time, and it
works. Then, I put the code in a “for” loop. However, it seems that the
“O.engines” does not start when the unloading test is to be done, because no
information about the running time is shown in the terminal window and no file
is outputted (I introduce two PyRunners into O.engines).
Could you help me solve this problem?
##########################################
# unicode: UTF-8
# For 2D biaxial simulation
# 19/09/2022
#########################################
### Defining parameters and variables ###
#########################################
#Material constants
FrictionAngle = 35.000000000000000
Damp = 0.5
Porosity = 0.0
#Confining variables
ConfPress = -1.0e5
#Loading control
LoadRate = 0.01
#time calculation
startT = O.time
endT = O.time
timeSpent = endT - startT
########################################
#import necessary packages
from yade import pack,plot,os,timing
#record of e22
recorde =
['0.00010','0.00020','0.00031','0.00041','0.00051','0.00061','0.00071','0.00082','0.00092','0.00102',\
'0.00112','0.00122','0.00132','0.00143','0.00153','0.00163','0.00173','0.00183','0.00194','0.00204',\
'0.00214','0.00224','0.00234','0.00245','0.00255','0.00265','0.00275','0.00286','0.00296','0.00306',\
'0.00316','0.00326','0.00337','0.00347','0.00357','0.00367','0.00377','0.00388','0.00398','0.00408',\
'0.00418','0.00428','0.00439','0.00449','0.00459','0.00469','0.00479','0.00490','0.00500']
for iii in recorde:
#folders for post-processing
try:
os.mkdir('contact-unload'+str(iii))
except:
pass
#restart
O.load('save'+str(iii)+'yade.bz2') # obtained from "work
calculation_dt=0.00001"
for k in O.bodies:
k.state.vel=(0,0,0) # set the initial speed to be 0
k.state.angVel=(0,0,0)
triax1=O.engines[4]
#O.run(100,1)
def Output():
e22=-triax1.strain[1]
#particle info
g =
open('./contact-unload'+str(iii)+'/pInfo'+'{:.5f}'.format(e22)+'.txt', 'w')
print('particle information at Iter %d' % O.iter, file=g)
g.write('left right bottom top\n')
print(O.bodies[0].state.pos[0],O.bodies[1].state.pos[0],O.bodies[2].state.pos[1],O.bodies[3].state.pos[1],file=g)
g.write('ID x y radius mass fx fy disx disy rotz velx vely omega\n')
for b in O.bodies:
if isinstance(b.shape,Sphere):
print(b.id, b.state.pos[0], b.state.pos[1], b.shape.radius,
b.state.mass, O.forces.f(b.id)[0], O.forces.f(b.id)[1], b.state.displ()[0],
b.state.displ()[1], b.state.rot()[2], b.state.vel[0], b.state.vel[1],
b.state.angVel[2], file=g)
g.close()
Output()
#################################
#####Defining triaxil engines####
#################################
triax1=TriaxialStressController(
width0=0.9994224497343608,
height0=1.4991128144248458,
depth0=0.0400000000000000,
wall_back_activated = False, #for 2d simulation, z-axis is along the
direction of front and back
wall_front_activated = False,
thickness = 0.001,
#maxMultiplier = 1.002,
internalCompaction = False, # If true the confining pressure is
generated by growing particles
max_vel = 0.1,
stressMask = 5,
computeStressStrainInterval = 10,
goal1 = ConfPress,
goal2 = LoadRate,
#goal3 = ConfPress,
)
ini_e22 = -triax1.strain[1]
ini_e22i = -triax1.strain[1]
kkk = 1
# print('ini_e22',ini_e22) # The initial strain is always 0, before the
engine runs
O.usesTimeStepper=True
# O.dt=2e-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_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax1,
newton,
PyRunner(iterPeriod=100,command='endCheck()',label='check'),
PyRunner(command='addPlotData()',iterPeriod=10,label='record'),
# VTKRecorder(iterPeriod=10000,recorders=['all'],fileName='./post/p1-')
# TriaxialStateRecorder(iterPeriod=2000,file='WallStresses')
]
print('ini_e22',ini_e22)
# Simulation stop conditions defination
def endCheck():
#unb=unbalancedForce()
e22=-triax1.strain[1]
if e22<0.0 or triax1.stress(triax1.wall_top_id)[1]>ConfPress:
O.pause()
endT = O.time
timeSpent = endT - startT
print('iii',iii)
print('Time using:',timeSpent)
plot.saveDataTxt('unloadinglog'+str(iii)+'.txt')
# collect history of data
def addPlotData():
global ini_e22i,ini_e22
global kkk
e22=-triax1.strain[1]
if kkk:
ini_e22i = -triax1.strain[1]
ini_e22 = -triax1.strain[1]
kkk = 0
#print('e22',e22,'ini_e22',ini_e22,'ini_e22i',ini_e22i)
if abs(ini_e22i) - abs(e22) < 0.00025:
test = abs(ini_e22) - abs(e22) > 0.000001
elif abs(ini_e22i) - abs(e22) < 0.08:
test = abs(ini_e22) - abs(e22) > 0.0001
if test:
unb = unbalancedForce()
mStress =
-(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2.0
area =
(O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])
Porosity = 1.0-sum(pi*b.shape.radius**2 for b in O.bodies if
isinstance(b.shape,Sphere))/area
plot.addData(e11=-triax1.strain[0], e22=-triax1.strain[1],
e33=-triax1.strain[2],
ev=-triax1.strain[0]-triax1.strain[1]-triax1.strain[2],
s11=-triax1.stress(triax1.wall_right_id)[0],
s22=-triax1.stress(triax1.wall_top_id)[1],
s33=-triax1.stress(triax1.wall_front_id)[2],
ub=unbalancedForce(),
dstress=(-triax1.stress(triax1.wall_top_id)[1])-(-triax1.stress(triax1.wall_right_id)[0]),
disx=triax1.width,
disy=triax1.height,
disz=triax1.depth,
i=O.iter,
porosity=Porosity,
cnumber=avgNumInteractions(),
total=O.energy.total(),
**O.energy
)
print('axial strain',e22,'unbalanced force:',unb,'mean stress:
',mStress,'s11',-triax1.stress(triax1.wall_right_id)[0],'s22',-triax1.stress(triax1.wall_top_id)[1],'coordination
number: ',avgNumInteractions(),'porosity: ',Porosity)
Output()
ini_e22 = e22
if iii==0:
O.run(); #O.wait()
--
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 : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
