yade-users team mailing list archive

[Alireza Sadeghi <question686024@xxxxxxxxxxxxxxxxxxxxx>]

New question #686024 on Yade:
https://answers.launchpad.net/yade/+question/686024

Hello all,

I want to simulate a sample in triaxial test. The sample has particles with three different sizes. The problem is that when I have only mono-size paricles, I can see the micro strain in paraview, but when I have three different size particles, the paraview does not show me the correct result. I was wondering if you could check my code. 
Thank you very much.

Best Regards,

Alireza

####### code######

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

#==============================================================
#================= 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='aggregate-48'))

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

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= False, young=2e11, 
density=1523.6e2, poisson=0.3, frictionAngle= 0.28, 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','t','t']

mats=['aggregate-48','aggregate-814','aggregate-814']

coke=((1.875e-3,500),(0.9475e-3,1367),(0.50125e-3,10618))

color=((0,0,1),(0,1,0),(1,1,0))

tolerance=[(2e-4),(1e-4),(1e-5)]

for i in range(len(nums)):

     nums[i]=pack.SpherePack()
     nums[i].makeCloud((0,0,0),(2e-1,2e-1,2e-1),rMean=coke[i][0],rRelFuzz=tolerance[i],num=coke[i][1])
     O.bodies.append([utils.sphere(c,r,material=mats[i],color=color[i]) for c,r in nums[i]])


#===============================================================
#=================== 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=1000,file='WallStresses'),
#        VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=1000),
#        qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),
	NewtonIntegrator(damping=0.4,gravity=[0,0,-9.81])
]

yade.qt.Controller(), yade.qt.View()

O.dt=5e-8
#=================================================================
#=================  APPLYING CONFINING PRESSURE   ================
#=================================================================
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.38 and abs(sigmaIso-triax.meanStress)/abs(sigmaIso)<0.01:

    print 'Isotropic strain1:',-triax.strain[0], 'Isotropic strain 2:',-triax.strain[1], 'Isotropic strain 3:',-triax.strain[2]
    break
    print "###      Isotropic state saved      ###"


triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
O.save('RVE-sizeDis-solid-Isoe5-Isopart.yade')

#====================================================================
#=====================   DEVIATORIC LOADING   =======================
#====================================================================

triax.stressMask = 3
triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=-1

O.saveTmp()

#====================================================================
#====================     Micro Strain    ===========================
#====================================================================

TW=TesselationWrapper()
TW.computeVolumes()
TW.volume(20)
TW.setState(0)
step=0
while 1:
    step +=1
    O.run(100,True)
    TW.setState(1)
    TW.defToVtk("strain_"+str(step)+".vtk")
    
O.engines=O.engines+[PyRunner(iterPeriod=20000,command='finishIt()',label='calmEngine')]

def finishIt():
	if abs(-triax.strain[2]) > 0.5:
#		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('RVE-sizeDis-solid-Isoe5-Devpart.yade')
		O.pause()


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

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

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

def history():
  	plot.addData(exx=-triax.strain[0], eyy=-triax.strain[1], ezz=-triax.strain[2], ev=(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], mass=sum(b.state.mass for b in O.bodies if isinstance(b.shape,Sphere)), q=-(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0]), p=triax.meanStress, R=triax.stress(triax.wall_front_id)[2]/sigmaIso, por=triax.porosity, i=O.iter),

plot.saveDataTxt('RVE-sizeDis-solid-Isoe5-Devpart.txt')

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