yade-users team mailing list archive

[Chien-Cheng Hung <question696150@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Chien-Cheng Hung is still having a problem:
Hi Bruno,

Thanks for your quick reply.

I've minimized my script as much as possible. I also tested it and the
running time should be within 30 seconds.

Below is the MWE.

###
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=3*(2*RADIUS1)
height=3*(2*RADIUS1)
width=3*(2*RADIUS1)
PI=1.e5
SN=5.e6
RATE_shear=1
wallFRIC=0
boundaryFRIC=0.5 # during compaction (controls porosity)
spFRIC=0.5
DAMPSHEAR=0.

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+RADIUS1,1.5*width),orientation=Quaternion(1,0,0,0),extents=(5*length,RADIUS1/2,5*width),fixed=1,wire=False,color=(1,0,0),material='boxMat')
lowBox = utils.box(center=(length/2,height-RADIUS1,1.5*width),orientation=Quaternion(1,0,0,0),extents=(5*length,RADIUS1/2,5*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*RADIUS1,width),(length,height+0.1*RADIUS1-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)
 ,PeriTriaxController(dynCell=True,mass=10,maxUnbalanced=1e-3,relStressTol=1e-4,stressMask=7,goal=(-PI/volRatio,-PI/volRatio,-PI/volRatio),globUpdate=1,maxStrainRate=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ,PyRunner(command='fixVelocity(RATE_shear)',iterPeriod=1,label='fixVel',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(0.1,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(100,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)

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.1:
        O.pause()
        
O.run(1000000000,1)
 
### Micro-strain analysis
TW=TesselationWrapper()
TW.triangulate()        
TW.computeVolumes()     
TW.volume(10)           
TW.setState(0)          
O.run(100,True)         
TW.setState(1)          
TW.defToVtk("strain.vtk")

###

Cheers,
Chien-Cheng

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