yade-users team mailing list archive

[Alireza Sadeghi <question686273@xxxxxxxxxxxxxxxxxxxxx>]

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

Hello All,

I want to simulate a triaxial test of sample of granular material with real shape. I am using clumps instead of spheres. I want to plot micro strain of sampel in paraview. Is it possible to get micro strain when the particles are clumps?
Because my code does not work.

Thank you very much for your help in advance.

Best Regards,

Alireza


########my code#########
O.reset()

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

#==================================COKE AGGREGATE CLUMP TEMPLATES================================================================
	
radz1=[0.355155e-3,0.505113e-3,0.397713e-3,0.465286e-3,0.484395e-3,0.394534e-3,0.493151e-3,0.487328e-3,0.613619e-3,0.413455e-3]
poz1= [[1.13418e-3,-0.703895e-3,-1.20338e-3],[-0.390408e-3,0.476061e-3,-0.150612e-3],[-0.556545e-3,0.451341e-3,1.1495e-3],[-0.633942e-3,0.498253e-3,0.348231e-3],[0.0256934e-3,0.388855e-3,-0.733445e-3],[-0.218563e-3,0.504478e-3,1.54117e-3],[-0.319601e-3,0.104778e-3,0.742895e-3],[0.650678e-3,-0.76675e-3,-0.289908e-3],[0.0113115e-3,-0.207684e-3,0.00255944e-3],[0.594902e-3,-0.301473e-3,-0.878654e-3]]		 
template1= []
template1.append(clumpTemplate(relRadii=radz1,relPositions=poz1))

radz2=[0.330164e-3,0.504115e-3,0.399587e-3,0.614205e-3,0.466444e-3,0.495302e-3,0.394324e-3,0.486898e-3,0.444037e-3,0.489396e-3]
poz2= [[1.16386e-3,-0.70706e-3,-1.25222e-3],[-0.39596e-3,0.482385e-3,-0.144097e-3],[-0.554928e-3,0.448475e-3,1.14685e-3],[-0.00361766e-3,-0.198211e-3, 0.00106559e-3],[-0.633362e-3,0.490424e-3,0.357391e-3],[0.0434148e-3,0.367924e-3,-0.736319e-3],[-0.218749e-3,0.504703e-3,1.54165e-3],[-0.311777e-3, 0.101954e-3,0.750235e-3],[0.621565e-3,-0.779387e-3,-0.179029e-3],[0.711114e-3,-0.503202e-3,-0.795602e-3]]		 
template2= []
template2.append(clumpTemplate(relRadii=radz2,relPositions=poz2))

radz3=[0.959101e-3,0.774782e-3,0.924682e-3,0.882986e-3,0.572207e-3,0.875338e-3,0.499054e-3,0.582586e-3,1.03184e-3,0.561428e-3]
poz3=[[-0.742455e-3,-0.539002e-3,0.030705e-3],[0.800775e-3,1.00778e-3,0.366095e-3],[-0.217508e-3,-0.423924e-3,-0.671114e-3],[0.65016e-3,0.741405e-3,0.305558e-3],[-0.770717e-3,0.423006e-3,0.33259e-3],[0.20475e-3,-0.707669e-3,-0.332745e-3],[-0.0356616e-3,1.34122e-3,0.703809e-3],[0.64889e-3,0.340976e-3,-0.456228e-3],[0.258817e-3,0.13357e-3,0.178886e-3],[0.35466e-3,1.36934e-3,0.583507e-3]]
template3= []
template3.append(clumpTemplate(relRadii=radz3,relPositions=poz3))		 
		 
radz4=[1.18842e-3,1.0381e-3,0.664711e-3,0.531753e-3,0.853808e-3,0.789023e-3,0.717061e-3,1.0081e-3,0.967001e-3,0.535189e-3]
poz4=[[-0.0354242e-3,-0.245642e-3,-0.0514299e-3],[0.444228e-3,-0.342353e-3,-0.327639e-3],[-0.815227e-3,1.41647e-3,-0.246128e-3],[1.38493e-3,-0.405625e-3,-0.69207e-3],[0.608272e-3,0.0171296e-3,0.863712e-3],[-0.191982e-3,-1.05752e-3,-0.341411e-3],[-1.08258e-3,-0.336001e-3,-0.482846e-3],[-0.396508e-3,0.47029e-3,0.0820364e-3],[0.486211e-3,0.345623e-3,0.492724e-3],[-0.0502768e-3,0.977767e-3,0.467401e-3]]
template4= []
template4.append(clumpTemplate(relRadii=radz4,relPositions=poz4))	


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

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


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'))

#========= creating walls ======================

#fc=yade.geom.facetCylinder((0,0,0.05),12.7e-3, 0.1, orientation=Quaternion((0, 0, 1), 0), segmentsNumber=100, wallMask=7, angleRange=None, closeGap=False, radiusTopInner=-1, radiusBottomInner=-1,material='wall')
#O.bodies.append(fc)


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


############################
###   DEFINING ENGINES   ###
############################


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'),
	NewtonIntegrator(damping=0.4,gravity=[0,0,-10])
	,PyRunner(command='calm()',iterPeriod=10,label='calmEngine')
]

O.dt=1e-7


#====================================================Generating the aggregates===================================================

coke=((1.875e-3,100),(0.9475e-3,367),(0.50125e-3,500),(0.50125e-3,500))

nums=['t','t','t','t']

temps=[template1,template2,template3,template4]


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



for i in range(len(nums)):
    nums[i]=pack.SpherePack()
    nums[i].makeCloud((-0.045,-0.045,-0.045),(0.045,0.045,0.045),rMean=coke[i][0],rRelFuzz=0.0,num=coke[i][1])
    O.bodies.append([utils.sphere(c,r,material=mats[i]) for c,r in nums[i]])
    O.bodies.replaceByClumps(temps[i],[1.0],discretization=5)


O.step()
print 'number of interactions  is ',len([i for i in O.interactions])
print ''

print ''
print ''


############################
###   Display settings   ###
############################

for x in range(len(O.bodies)):
	if (O.bodies[x]):
                if isinstance(O.bodies[x].shape,Sphere):
			if O.bodies[x].isStandalone:
				O.bodies[x].shape.color=(0.2,0.3,0.6)
			else:
				O.bodies[x].shape.color=(0.7,0.63,0.0)
		else:
			O.bodies[x].shape.color=(0.2,0.3,0.1)


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


O.engines=O.engines+[PyRunner(command='check_it()',iterPeriod=500,label='checkEngine')]



def check_it():
	print 'number of interactions is ',len([i for i in O.interactions])
	pp=max(i.geom.penetrationDepth for i in O.interactions)
	print 'max pen is ',pp
	print 'unbalanced force is:   ',unbalancedForce()
	print ''




#O.run(5000,True)

#===============================================
#=============== Compaction ====================
#===============================================


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


O.engines=O.engines+[triax]

triax.goal1=-1.0e5
triax.goal2=-1.0e5
triax.goal3=-1.0e5
triax.wall_back_activated=True
sigmaIso=-1.0e5
O.dt=2e-6

calmEngine.dead=True
checkEngine.dead=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.385 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 = -0.1



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


TW=TesselationWrapper()
TW.computeVolumes()
TW.volume(20)
TW.setState(0)
step=0
while 1:
    step +=1
    O.run(1000,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.