How to get Micro-strain field from a 2D simulation in YADE

Hello Leonard,

I think that the main problem is 2D formulation. I would try to prepare thicker samples (artificially replicate your 2D "layers").

Cheers,
Karol

Hi Karol,

Thank you for your reply.

Could you please give more hints? Sorry that I didn't get the answer.

Thanks
Leonard

Hi Leonard,

The algorithm in Yade builds a 3D tesselation. It is hard for the algorithm to build it if all the particles lie on the same plane. Maybe you could add some 'thickness' to your simulation by creating three identical sets of particles on three parallel planes.

Of course, it is not efficient to do the simulation this way. But you can store the information about your particles ("one layer") in separate files during simulation [1]. After this, you could prepare another script just for loading the particles into the simulation at certain positions. However, in the second script, you can load the previously stored 'layer' three times (two times with shifted position, using the shift option [2]). That way, you will have some "thickness" of the sample.

Cheers,
Karol

[1] https://yade-dem.org/doc/yade.export.html?highlight=export#yade.export.text
[2] https://yade-dem.org/doc/yade.ymport.html?highlight=ymport#module-yade.ymport

Hi Leonard,

Try this (please note that I am using python 3, so the call pf print function is modified):
##########
################### 2D MWE #############
from __future__ import division
from yade import pack, plot, export, ymport

num_spheres=1000

rate=-0.01
damp=0.6
stabilityThreshold=0.001
young=5e6
confinement=6.7e3

mn,mx=Vector3(0,0,0.02),Vector3(1,2,0.02)

O.materials.append(FrictMat(young=young, poisson=0.5, frictionAngle=radians(30), density=2600, label='spheres'))
O.materials.append(FrictMat(young=young, poisson=0.5, frictionAngle=0, density=0, label='walls'))
walls=aabbWalls([Vector3(0,0,0),Vector3(1,2,0.04)],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.75,seed=1)

O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

Gl1_Sphere.quality=3

for b in O.bodies:
        if isinstance(b.shape,Sphere):
                b.state.blockedDOFs='zXY'
                b.shape.color=[1,1,1]

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

newton=NewtonIntegrator(damping=damp)

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Box_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom(), Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
        GlobalStiffnessTimeStepper(active=1, timeStepUpdateInterval=100, timestepSafetyCoefficient=0.8),
        triax,
        newton
]

triax.goal1=triax.goal2=triax.goal3=-confinement

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print('unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1],'left:',-triax.stress(triax.wall_left_id)[0],'front:',-triax.stress(triax.wall_front_id)[2])
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
          break

print("### state 1 completed ###")

######## export files and run simulation
export.text('state_0.txt')## export files state 0
O.run(500,True)
export.text('state_1.txt')## export files state 1

######################
layerThickness = 0.015*2 # this is more or less diameter of your particle
O.bodies.clear()
for i in range(3):
    bb = ymport.text('state_0.txt', shift = Vector3(0,0,i*layerThickness))
    O.bodies.append(bb)
TW=TesselationWrapper()
TW.triangulate()
TW.computeVolumes()
TW.setState(0)

O.bodies.clear()
for i in range(3):
    bb = ymport.text('state_1.txt', shift = Vector3(0,0,i*layerThickness))
    O.bodies.append(bb)
TW.setState(1)
TW.defToVtk("strain_pseudo2D.vtk")

#######
Cheers,
Karol

Thanks Karol, this is a great idea!

Cheers,
Leonard