Eplastic.py

##### part 1 : generate a soil sample by gravity deposit

from yade import pack,plot
from yade import utils
import random


finalFricDegree = 30
O.trackEnergy=True
young=100e6
rmean=0.1

rmin=0.05
rmax=0.15

mn=Vector3(0,0,-2*rmean)
mx=Vector3(7*2*rmean,5*2*rmean,2*rmean)

number=30
fallingHeight=5

O.materials.append(CohFrictMat(young=100e6,poisson=0.3,frictionAngle=radians(0),density=2600,normalCohesion=0,shearCohesion=0,alphaKr=0.5,alphaKtw=0.01,momentRotationLaw=False,isCohesive=False,etaRoll=-1,label='spheres'))

O.materials.append(FrictMat(young=100e6,poisson=0.3,frictionAngle=0,density=0,label='walls'))


walls=utils.aabbWalls([mn,mx],thickness=0.5*rmean, oversizeFactor=1.0, material='walls')
wallIds=O.bodies.append(walls)


ssp=pack.SpherePack()

for i in range(1,number):
    r =rmin + (rmax-rmin)*(-0.15*log(1-random.random()))**1.33333
    if r>rmax:
     r =rmin + (rmax-rmin)*(-0.15*log(1-random.random()))**1.33333
    mbegin=(0,0,-1*r)
    mend=(mx[0],mx[1],1*r)
    ssp.makeCloud(minCorner=mbegin,maxCorner=mend,rMean=r,rRelFuzz=0,num=1)

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


newton=NewtonIntegrator(damping=0.9,kinSplit=True,gravity=(0,0,0))

O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
                [Law2_ScGeom_FrictPhys_CundallStrack(label='law1'),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
                              useIncrementalForm=False,
                              always_use_moment_law=False,
                              label='law2')]
    ),
    GravityEngine(gravity=(0,-9.81,0)),
    newton

]

for b in O.bodies:
    if isinstance(b.shape,Sphere): b.state.blockedDOFs='zXY' # blocked movement along Z


O.dt=0.5*utils.PWaveTimeStep()

O.trackEnergy=True


while 1:
    O.run(1000,True)
    unb=unbalancedForce()
    if unb < 0.001:
        break

## reset boundaries
Mz=max([b.state.pos[1]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)])
h=Mz/4

mn=(0,0,-2*rmean)
mxafter=Vector3(mx[0],Mz+2*rmean,2*rmean)

O.bodies.erase(0)
O.bodies.erase(1)
O.bodies.erase(2)
O.bodies.erase(3)
O.bodies.erase(4)
O.bodies.erase(5)

O.bodies.append(utils.aabbWalls([mn,mxafter],thickness=0.1*rmean, oversizeFactor=1.0, material='walls'))

O.bodies.erase(3)
O.bodies.erase(4)
O.bodies.erase(5)

## add friction and moment law
O.materials[0].frictionAngle=radians(finalFricDegree)
O.materials[0].momentRotationLaw=True
O.engines[2].lawDispatcher.functors[1].always_use_moment_law = True

O.engines[4].damping=0.000

O.run(1000, True)
O.saveTmp()


#### part 3: add the boulder and conduct impact modelling, track energy dissipations

boulderradius=1.5*rmean

sp2=pack.SpherePack()

mbegin=(0.5*mxafter[0]*1.0-boulderradius,mxafter[1],-1*boulderradius)
mend=(0.5*mxafter[0]*1.0+boulderradius,mxafter[1]+2
*boulderradius,1*boulderradius)
sp2.makeCloud(minCorner=mbegin,maxCorner=mend,rMean=boulderradius,rRelFuzz=0,num=1)
O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp2])

O.bodies[3].state.blockedDOFs='zXY'

O.bodies[3].state.vel=Vector3(0,-5,0)

O.dt=0.1*utils.PWaveTimeStep()

for i in range (1,2000):
    O.run(20,True)
    if -O.forces.f(3)[1]< O.bodies[3].state.mass*9.81:
      tt=O.time
      break


O.run(1000,True)
tt=O.time

O.saveTmp()

tt=O.time
f=open('epTrack.txt','w')
O.step()
e=0
aa=O.energy['plastDissip']

for k in range (1,100):

    for i in O.interactions:
        Ft=i.phys.shearForce
        du=i.geom.shearInc
        kt=i.phys.ks
        kn=i.phys.kn
        Fn=i.phys.normalForce
            Ftt=Ft-kt*du
        maxFs=Fn.norm()*i.phys.tangensOfFrictionAngle
        if(Ftt.norm() > maxFs):
                ratio = maxFs / Ftt.norm()
                trialForce=Ftt
                Ftt *= ratio
                e=e+((1/kt)*(trialForce-Ftt)).dot(Ftt)

    O.step()
    E=O.energy['plastDissip']-aa
    f.write('%f %f %f \n' %(O.time-tt,e,E))


f.close()