micro strain in clumps

Asked by Alireza Sadeghi on 2019-11-26

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

Question information

Language:
English Edit question
Status:
Open
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-11-27
Last reply:
2019-11-27
Robert Caulk (rcaulk) said : #1

>Because my code does not work.

Hello please share with us what happens when you execute your code. Error? If so what is the error? Would you please be kind enough to review and adhere to our forum guidelines? You can find them here [1]

Thank you very much,

Robert

[1]https://www.yade-dem.org/wiki/Howtoask

Alireza Sadeghi (asadeghime) said : #2

Dear Robert,

Thank you very much for your answer. It does not show any errors in the terminal. Just I have something similar this:

.
.
.
10398 : Vh==NULL!! id=10398 Point=-0.0135973 0.0130332 0.00710087 rad=0
10409 : Vh==NULL!! id=10409 Point=0.00605096 -0.00551408 -0.00619671 rad=0
Triangulated Grains : 8687
sym_grad_u_total_g (wrong averaged strain):
0.451431 -0.0578464 -0.166245
-0.0578464 0.562565 -0.0421092
-0.166245 -0.0421092 0.16872

Total volume = 0.168629, grad_u =
0.451431 -0.0917767 -0.254586
-0.0239162 0.562565 -0.0416519
-0.0779029 -0.0425665 0.16872

sym_grad_u (true average strain):
0.451431 -0.0578464 -0.166245
-0.0578464 0.562565 -0.0421092
-0.166245 -0.0421092 0.16872

Macro strain :
0 0 0
0 0 0
0 0 0

I can not open the vtk files in paraview. Thanks a lot.

Best Regards,

Alireza

Can you help with this problem?

Provide an answer of your own, or ask Alireza Sadeghi for more information if necessary.

To post a message you must log in.