yade-users team mailing list archive

[Chien-Cheng Hung <question691407@xxxxxxxxxxxxxxxxxxxxx>]

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

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

Thanks for your email.

I've tested your microstrain vtk file and it works with my Paraview. So
I think the crash is not related to the Paraview version (btw my
Paraview version is 5.0.1 (64-bit)). Maybe it's related to my code.
Could you help me have a look? I'd appreciate that.

### My simplified code (direct shear simulation of a granular layer)
from yade import pack,plot,export
import numpy as np
import math

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

O.periodic=True

# dimensions of sample
RADIUS=0.25
a=15
b=1
c=1
length=a*(2*RADIUS)
height=length/b
width=length/c
thickness=RADIUS

### Guassian distribution
psdSizes=[.456,.5,.544]
psdCumm=[0,0.5,1]

# friction angles 
wallFRIC=0
boundaryFRIC=0.5
spFRIC=0.5

# boundary conditions
PI=1.e5
SN=5.e6 # normal stress
RATE_NS1=1 # velocity of top plate during compaction
RATE_NS2=1 # velocity of top plate during shear
RATE_shear=1 # shear velocity
roll_stiff=0
roll_fric=0

# simulation control
DAMPSHEAR=0.

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

O.materials.append(FrictMat(density=3000,young=1e8,poisson=0.5,frictionAngle=wallFRIC,label='boxMat'))
O.materials.append(FrictMat(density=3000,young=1e8,poisson=0.5,frictionAngle=boundaryFRIC,label='boundaryMat'))
O.materials.append(FrictMat(density=3000,young=1e8,poisson=0.5,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+3*RADIUS,width),(length,2*height-3*RADIUS,2*width), psdSizes =psdSizes, psdCumm =psdCumm, periodic=True, seed =1)
sp2.makeCloud((0,height+RADIUS,width),(length,height+RADIUS-1e-10,2*width), rMean=RADIUS, periodic=True, seed =1)
sp3.makeCloud((0,2*height-RADIUS,width),(length,2*height-RADIUS-1e-10,2*width), rMean=RADIUS, 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(krot=roll_stiff,eta=roll_fric)],
  [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=(1,1,1),doneHook='triaxDone()',label='triax')
 ,NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')
 ]

def triaxDone():
 global phase
 volRatio=(O.cell.hSize[0,0]*O.cell.hSize[1,1]*O.cell.hSize[2,2])/((O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1])*O.cell.hSize[0,0]*O.cell.hSize[2,2])
 h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
 vol=h*O.cell.hSize[0,0]*O.cell.hSize[2,2]
 contactStress=getStress(vol)
 vol_s=Rmean=Rmax=nbSph=0
 Rmin=1e6
 x_ref=O.bodies[0].state.pos[0]
 for o in O.bodies:
  if isinstance(o.shape,Sphere):
   nbSph+=1
   Rmean+=o.shape.radius
   if o.shape.radius>Rmax: Rmax=o.shape.radius
   if o.shape.radius<Rmin: Rmin=o.shape.radius
   vol_s += 4.*pi/3.*(o.shape.radius)**3
 Rmean=Rmean/nbSph
 n = 1-vol_s/vol
 print ('DONE! iter=',O.iter,'| sample generated: nb spheres',nbSph,', Rmean=',Rmean,', Rratio=',Rmax/Rmin,', porosity=',n,'wall_refPos=',x_ref)
 tt=TriaxialCompressionEngine()
 #tt.setContactProperties(shearFRIC)
 triax.dead=True
 O.pause()


#### If top particle layer contact with wall then give the particle velocity, and vice versa.
def identify_topLayer():
    top = []
    for i in O.interactions.withBody(0):
        if i.id1 != 0:
            top.append(i.id1)
        if i.id2 != 0:
            top.append(i.id2)
    return top

O.run(10000000,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()
  print ('SN=',SN,'| current SN = ',currentSN,' | unbF=',unbF )
  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.005 ):
   stage+=1
   fnPlaten=O.forces.f(0)[1]
   ## the following computes the stiffness of the plate (used for stress control of the top plate)
   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(10000000,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)

### Start shearing
newton.damping=DAMPSHEAR
fixVel.dead = False

refTime = O.time

def fixVelocity(RATE_shear):
    O.bodies[0].state.vel[0] = RATE_shear
    topLayer_id  = identify_topLayer()
    for i in topLayer_id:
        O.bodies[i].state.vel[0] = RATE_shear
    slip = (O.time-refTime)*RATE_shear
    h=O.bodies[0].state.pos[1]-O.bodies[1].state.pos[1]
    ss = slip/h
    if ss > 1:
        O.pause()

O.run(10000,1)

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.