yade-users team mailing list archive

[Chien-Cheng Hung <question695851@xxxxxxxxxxxxxxxxxxxxx>]

Question #695851 on Yade changed:
https://answers.launchpad.net/yade/+question/695851

    Status: Needs information => Open

Chien-Cheng Hung gave more information on the question:
Hi Janek,

I've tried to simplify my script as possible as I can. The MWE is below.
Thank you for your time.

Cheers,
Chien-Cheng

###
from yade import pack,plot,export
import numpy as np
import math

sp1=pack.SpherePack()
sp2=pack.SpherePack()
sp3=pack.SpherePack()

O.periodic=True

RADIUS1=0.25    
RADIUS2=0.125
length=5*(2*RADIUS1)
height=5*(2*RADIUS1)
width=5*(2*RADIUS1)
thickness=RADIUS1

# boundary conditions
PI=1.e5
SN=5.e6
RATE_NS1=0.1
RATE_NS2=100
RATE_shear=1

# friction angles 
wallFRIC=0
boundaryFRIC=0.5 # during compaction (controls porosity)
spFRIC=0.5

# simulation control
DAMPSHEAR=0.

##### Microproperties

O.cell.hSize=Matrix3(length,0,0,0,3*height,0,0,0,width)

O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=wallFRIC,label='boxMat'))
O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=boundaryFRIC,label='boundaryMat'))
O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=spFRIC,label='sphereMat'))

upBox = utils.box(center=(length/2,2*height+thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2,height-thickness,1.5*width),orientation=Quaternion(1,0,0,0),extents=(2*length,thickness/2,width),fixed=1,wire=False,color=(1,0,0),material='boxMat')

O.bodies.append([upBox,lowBox])

sp1.makeCloud((0,height+1*RADIUS1,width),(length,2*height-1*RADIUS1,2*width), rMean=RADIUS2, periodic=True, seed =1)
sp2.makeCloud((0,height+0.1*thickness,width),(length,height+0.1*thickness-1e-10,2*width), rMean=RADIUS2, periodic=True, seed =1)
sp3.makeCloud((0,2*height-0.1*RADIUS1,width),(length,2*height-0.1*RADIUS1-1e-10,2*width), rMean=RADIUS2, periodic=True, seed =1)

sphere_id = O.bodies.append([utils.sphere(s[0],s[1],color=(0,0,1),material='sphereMat') for s in sp1])
bottomLayer_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,0,1),material='boundaryMat') for s in sp2])
topLayer_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,0,1),material='boundaryMat') for s in sp3])

effCellVol=(O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2]
volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/effCellVol

O.engines=[
 ForceResetter()
 ,InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb()],verletDist=-0.1,allowBiggerThanPeriod=True)
 ,InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D()],
  [Ip2_FrictMat_FrictMat_MindlinPhys()],
  [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True)]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=-1)
 ,PyRunner(command='fixVelocity(RATE_shear)',iterPeriod=1,label='fixVel',dead=True)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(10,10,10),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ,VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=1000,label='recVTK',dead=True)
 ]

def triaxDone():
 triax.dead=True
 O.pause()

O.run(10000000000,1)

stage=0
stiff=fnPlaten=currentSN=0.
def servo():
 global stage,stiff,fnPlaten,currentSN
 if stage==0:
  currentSN=O.forces.f(0)[1]/(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  unbF=unbalancedForce()
  boundaryVel=copysign(min(RATE_NS1,abs(0.5*(currentSN-SN))),currentSN-SN)
  O.bodies[0].state.vel[1]=boundaryVel
  if ( (abs(currentSN-SN)/SN)<0.001 and unbF<0.001 ):
   stage+=1
   fnPlaten=O.forces.f(0)[1]
   for i in O.interactions.withBody(O.bodies[0].id):
    stiff+=i.phys.kn
   O.pause()
 if stage==1:
  fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  boundaryVel=copysign(min(RATE_NS2,abs(0.333*(O.forces.f(0)[1]-fnDesired)/stiff/O.dt)),O.forces.f(0)[1]-fnDesired)
  O.bodies[0].state.vel[1]=boundaryVel
  
O.engines = O.engines[:5]+[PyRunner(command='servo()',iterPeriod=1,label='servo')]+O.engines[5:]

O.run(10000000000,1)

#### fixed bottom bottomLayer_id
for i in bottomLayer_id:
    O.bodies[i].state.blockedDOFs='xyzXYZ'
    O.bodies[i].state.dynamics = False
    O.bodies[i].state.vel = Vector3(0,0,0)

recVTK.dead=False
newton.damping=DAMPSHEAR
fixVel.dead = False

initial_particle_pos = O.bodies[topLayer_id[0]].state.pos[0]

def fixVelocity(RATE_shear):
    O.bodies[0].state.vel[0] = RATE_shear
    for i in topLayer_id:
        O.bodies[i].state.vel[0] = RATE_shear
    slip = O.bodies[topLayer_id[0]].state.pos[0] - initial_particle_pos
    h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
    ss = slip/h
    if ss > 0.5:
        O.pause()
        
O.run(1000000000,1)
###

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.