yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #28349
[Question #703247]: Run several unloading tests together
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.
