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