← Back to team overview

yade-users team mailing list archive

Re: [Question #680144]: plot the deviatoric part of TriaxialStressController simulation

 

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

    Status: Answered => Open

Alireza Sadeghi is still having a problem:
Hello Jerome,

Thank you for your response and help. I will use  https://gitlab.com/yade-dev/trunk. Thank you.
 about my question, I used plot.resetData() but it did not solve my problem. The simulation of triaxial test has two parts. The particles compact with isometric pressure and in the second step we have deviatoric compaction. The problem is, strain in the first step is added to the strain in the second step. but I want to draw strain-Stress diagram just for the second step of triaxial test simulation. my code is in the below of the mail. Thank you very much for your help.

Best Regards

Alireza

#Triaxial compresion for finding a proper size for RVE (stress control
boundary condition)

from yade import utils, plot
from yade import pack, qt
from datetime import datetime


qtr=qt.Renderer()
qtr.bgColor=(1,1,1)

#==============================================================
#================= define the materials =======================
#==============================================================

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion=
1e20, isCohesive= True, young=6.81e8, density=1377.5, poisson=0.3,
frictionAngle= 0.31, fragile=False, label='Coke'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= False, young=4e9, 
density=1523.6, poisson=0.3, frictionAngle= 0.0, fragile=False, label='wall'))


#===============================================================
#=================== define walls ============================
#===============================================================

walls=aabbWalls([(0,0,0),(2e-1,2e-1,2e-1)],thickness=0.0003,oversizeFactor=1.0,material='wall')
wallIds=O.bodies.append(walls)

#===============================================================
#=================== define packing ============================
#===============================================================

nums=['t']

mats=['Coke']

coke=(5e-3,2000)

nums=pack.SpherePack()

nums.makeCloud((0,0,0),(2e-1,2e-1,2e-1),rMean=coke[0],rRelFuzz=1e-4,num=coke[1])

O.bodies.append([utils.sphere(c,r,material=mats[0],color=(0,0,1)) for
c,r in nums])


#===============================================================
#=================== define Engine =============================
#===============================================================

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
	InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
	[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
	[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
	),
#	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
#        VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=100),
#        qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),

	NewtonIntegrator(damping=0.4,gravity=[0,0,0])
]
O.dt=2e-6


#=================================================================
#=================  APPLYING CONFINING PRESSURE   ================
#=================================================================
stabilityThreshold=0.2
sigmaIso=-1e5

triax=TriaxialStressController(
	
	maxMultiplier=1.000,  
	finalMaxMultiplier=1.000,  
	thickness = 0,
	stressMask = 7,
	internalCompaction=False, 
)


O.engines=O.engines+[triax]

triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
triax.wall_back_activated=True


while 1:
  O.run(100, True)
  unb=unbalancedForce()
  meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  print 'mean Stress:',triax.meanStress,'porosity:', triax.porosity,'meanS:',unb
  if triax.porosity < 0.497 and ((sigmaIso-triax.meanStress)/sigmaIso) < 0.001:
    print 'Isotropic strain1:',-triax.strain[0], 'Isotropic strain 2:',-triax.strain[1], 'Isotropic strain 3:',-triax.strain[2]
    break

print "###      Isotropic state saved      ###"

#====================================================================
#=====================   DEVIATORIC LOADING   =======================
#====================================================================
plot.reset()
triax.stressMask = 3
triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=-10
#dataOld = plot.data
#plot.saveDataTxt('isotropic_part_deformableRVE.txt')
plot.resetData()
#plot.plot()

O.saveTmp()


#=================================================================
#======================= data collecting =========================
#=================================================================

O.engines=O.engines+[PyRunner(iterPeriod=16000,command='finishIt()',label='calmEngine')]

def finishIt():
#        if abs(-triax.strain[2]) > 0.2:
#           makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
	   print O.iter
	   print 'time is ', O.time
           print 'porosity:',triax.porosity
	   print 'strain 1:',-triax.strain[0], 'strain 2:',-triax.strain[1], 'strain 3:',-triax.strain[2]
#	   O.save('RVE200-deformable.yade')


def history():
  	plot.addData(exx=-triax.strain[0], eyy=-triax.strain[1], ezz=-triax.strain[2],
  		    ev5=(triax.strain[0]+triax.strain[1]+triax.strain[2]),
		    sxx=-triax.stress(triax.wall_right_id)[0],
		    syy=-triax.stress(triax.wall_top_id)[1],
		    szz=-triax.stress(triax.wall_front_id)[2],
                    q=-(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0]),
                    p=triax.meanStress,
                    R=3*(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0])/triax.meanStress,
                    por=triax.porosity,
		    i=O.iter),


if 1:
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]

else:
  O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')


O.run(50000,True)



plot.plots={('ezz'):('ev5')}
plot.plot()

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.