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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
