yade-users team mailing list archive

[Fu zuoguang <question226352@xxxxxxxxxxxxxxxxxxxxx>]

Question #226352 on Yade changed:
https://answers.launchpad.net/yade/+question/226352

    Status: Answered => Open

Fu zuoguang is still having a problem:
Dear Jan Stránský:
     Thanks for helping me last time and I have some new problems to ask your help, which are about getting output of the deviator stress in biaxial compression simulation. I wanna get the reaults of deviator stress using these script :

### fundamental details of application ###
# unicode: UTF-8 
Filename='2D-simulation'
from yade import pack,os
################################# preprocessing for simulation ##########################################  
### prescribing variables and functions for simulation controller ###
# material defination
spheremat = O.materials.append(ViscElMat(kn=4e5,ks=4e5,cn=.0,cs=.0,density=1500,frictionAngle=5.8))
wallmat = O.materials.append(ViscElMat(kn=4e5,ks=4e5,cn=.0,cs=.0,density=1500,frictionAngle=5.8))
# walls defination
mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0.1)
wallIds=O.bodies.append(utils.aabbWalls([mn,mx],thickness=.0001,material=wallmat))
# ThreeDTriaxialEngine defination for initial-state determination(the first calculation step)
triax01=ThreeDTriaxialEngine(
	wall_bottom_id=wallIds[2],wall_top_id=wallIds[3],
	wall_left_id=wallIds[0],wall_right_id=wallIds[1],
	wall_back_id=wallIds[4],wall_front_id=wallIds[5],
	wall_front_activated = False,wall_back_activated = False,
	internalCompaction=False, 
	stressControl_1 = True, stressControl_2 = True,stressControl_3 = True,
	computeStressStrainInterval =10,
	sigma_iso = 1.25e4,
	sigma1 = 1.25e4,
	sigma2 = 1.25e4,
	sigma3 = 1.25e4,
	strainRate1 = 0.03,strainRate2 = 0.03,
)
# ThreeDTriaxialEngine defination for triaxial compression(the second calculation step)
triax02=ThreeDTriaxialEngine(
	wall_bottom_id=wallIds[2],wall_top_id=wallIds[3],
	wall_left_id=wallIds[0],wall_right_id=wallIds[1],
	wall_back_id=wallIds[4],wall_front_id=wallIds[5],
	wall_front_activated = False,wall_back_activated = False,
	internalCompaction=False, 
	stressControl_1 = True, stressControl_2 = False,stressControl_3 = True,
	computeStressStrainInterval =10,
	sigma_iso = 1.25e4,
	sigma1 = 1.25e4,
#	sigma2 = 3,
#	sigma3 = 1.25e4,
	strainRate1 = -0.072,strainRate2 = 0.214,
)
# Simulation stop controller defination 
'''#def checkUnbalanced():
    unb=unbalancedForce()
    meanS=(triax01.stress(triax01.wall_right_id)[0]+triax01.stress(triax01.wall_top_id)[1])/2
    q=unb
    r=abs(meanS-triax01.sigma_iso)/triax01.sigma_iso
    if q<0.01 and r<1e-4:
       O.pause()
       O.save('first-step.xml'.format(Filename))'''
################################# control flow for simulation ##########################################  
# particles generation
O.periodic=1
O.cell.setBox(0.07,0.07,0.07)
sp=pack.SpherePack()
sp.makeCloud((0,0,.05),(0.07,0.07,.05),rMean=0.0008,rRelFuzz=0.3,num=840,periodic=True)
sp.toSimulation(material=spheremat)
# determining colors for particles in different aeras of the cell
for b in O.bodies:
    if isinstance(b.shape,Sphere):
         pos = b.state.pos
         if pos[0] <0.035 and pos[1] < 0.035: b.shape.color = (1,0,0)       # area 1
         elif pos[0] >= 0.035 and pos[1] <0.035: b.shape.color = (0,1,0)    # area 2
         elif pos[0] >= 0.035 and pos[1] >= 0.035: b.shape.color = (0,0,1)  # area 3
         elif pos[0] < 0.035 and pos[1] >= 0.035: b.shape.color = (0,10,1)  # area 4
         else: b.shape.pos = (1,1,0)         
O.periodic=0
# blockedDOFs
for b in O.bodies:
	if isinstance(b.shape,Sphere):
		 b.state.blockedDOFs='zXY'
# Simulation assembly for the first step
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
		[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
		[Law2_ScGeom_ViscElPhys_Basic()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax01,
	NewtonIntegrator(damping=.1),
#	PyRunner(command='checkUnbalanced()',iterPeriod=200)
]
# first step of simulation startting with a correct inheriting for the next step
O.dt = 2e-4
O.run(21300,True);
for b in O.bodies:
      b.state.vel=(0,0,0)
O.save('first-step.xml')
O.wait()
f = file("/home/fzg/fu/result-first.dat",'w')
f.write("varibles='X-refcordinate','Y-refcordinate','Z-refcordinate'\n\n")
f.write('%-16s %-16s %-16s\n'%('X-refcordinate','Y-refcordinate','Z-refcordinate'))
for b in O.bodies:
    if isinstance(b.shape,Sphere):
       refpos= b.state.refPos
       poscu = b.state.pos
       displ = b.state.displ()      
   #    pos0 = poscu[0]-0.035
    #   pos1 = poscu[1]-0.035
     #  pos2 = poscu[2]-0.035
       rad = b.shape.radius
       strainrate = b.state.vel
 #      vel0 = strainrate[0]*pos0
  #     vel1 = strainrate[1]*pos1
       f.write('%-16g %-16g %-16g\n'%(poscu[0],poscu[1],poscu[2]))
f.write('################################ ends ##########################################')
f.close()
f=file("/home/fzg/fu/macro.dat",'a')
f.write('%-16s %-16s %-16s %-16s %-16s %-16s %-16s\n'%('iters','time','area','axial-strain','volumetric-strain','porosity','DeviatorS')),
f.close()
# loading inheriting
O.load('first-step.xml')
# Simulation assembly for the second step
def macroresults():
	wall_left = O.bodies[0].state.pos[0]
	wall_right = O.bodies[1].state.pos[0]
	wall_bottom = O.bodies[2].state.pos[1]
	wall_top = O.bodies[3].state.pos[1]
	wall_back = O.bodies[4].state.pos[2]
	wall_front = O.bodies[5].state.pos[2]
	x_dimension = wall_right - wall_left
	y_dimension = wall_top - wall_bottom
	#z_dimension = wall_front - wall_back
	area = x_dimension * y_dimension
        volumetricstrain=(area-0.0021919)/0.0021919
	#area = (x_dimension-thick) * (y_dimension-thick) ## 'thick' is the thickness of the walls
	#vol = x_dimension * y_dimension * z_dimension
	porosity = 1-sum(pi*b.shape.radius**2 for b in O.bodies if isinstance(b.shape,Sphere))/area
        t=(O.iter-21300)*1e-4
        axialstrain =-(y_dimension-0.0468177)/0.0468177
        DeviatorS= (triax02.stress(triax02.wall_top_id)[1]-triax02.stress(triax02.wall_right_id)[0])
        f=file("/home/fzg/fu/macro.dat",'a')
        f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g\n'%(O.iter,t,area,axialstrain,volumetricstrain,porosity,DeviatorS)),
        f.close()
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
		[Ip2_ViscElMat_ViscElMat_ViscElPhys()],
		[Law2_ScGeom_ViscElPhys_Basic()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax02,
	NewtonIntegrator(damping=.1),
        PyRunner(iterPeriod=500,command='macroresults()')
]
# second step of simulation startting 
O.dt = 0.5e-4
macroresults()
O.run(16000,True);
# whole task over
O.save('final-step.xml'.format(Filename));
O.wait()
macroresults()
##################################################
f = file("/home/fzg/fu/result.dat",'w')
f.write('##  This is the result data of 2D simulation\n\n')
f.write('##  There are 12 types of varibles in this data as follows:\n\n')
f.write("varibles='X-refcordinate','Y-refcordinate','Z-refcordinate','X-cordinate','Y-cordinate','Z-cordinate','Radius','StrainRate-xx','StrainRate-yy','Velocity-xx','Velocity-yy','Ids'\n\n")
f.write('%-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s\n'%('X-refcordinate','Y-refcordinate','Z-refcordinate','X-cordinate','Y-cordinate','Z-cordinate','Radius','Velocity-xx','Velocity-yy','resi-Velocity-xx','resi-Velocity-yy','Ids'))
for b in O.bodies:
    if isinstance(b.shape,Sphere):
       refpos= b.state.refPos
       poscu = b.state.pos
       displ = b.state.displ()      
   #    pos0 = poscu[0]-0.035
    #   pos1 = poscu[1]-0.035
     #  pos2 = poscu[2]-0.035
       rad = b.shape.radius
       strainrate = b.state.vel
       resivel0 = strainrate[0]*(1.035-poscu[0])
       resivel1 = strainrate[1]*(1.035-poscu[1])
       f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g\n'%(refpos[0],refpos[1],refpos[2],poscu[0],poscu[1],poscu[2],rad,strainrate[0],strainrate[1],resivel0,resivel0,b.id))
f.write('################################ ends ##########################################')
f.close()
def rename():
    global Filename
    os.rename("/home/fzg/fu/result.dat","/home/fzg/fu/{0}.plt".format(Filename))
rename()

But the output text bring me a surprise that almost all the numbers of
deviator stress are 0, I can not determine where is wrong and seek your
help!

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.