yade-users team mailing list archive

[Chien-Cheng Hung <question689919@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Chien-Cheng Hung is still having a problem:
Dear Jan,

I simplified my script and hope that you could have a look.

#######
from __future__ import print_function
import builtins
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.025
a=15
b=1
c=1
length=a*(2*RADIUS)
height=length/b
width=length/c
thickness=RADIUS
psdSizes=[.0456,.05,.0544]
psdCumm=[0,0.5,1]

# friction angles
wallFRIC=0
spFRIC=26.6

# boundary conditions
PI=1.e5
SN=5.e6
RATE_shear=1

# simulation control
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=radians(wallFRIC),label='boxMat'))
O.materials.append(FrictMat(density=2500,young=5.5e10,poisson=0.25,frictionAngle=radians(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=True,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=True,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='sphereMat') for s in sp2])
topLayer_id = O.bodies.append([utils.sphere(s[0],s[1],color=(1,0,1),material='sphereMat') 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

class SomeCalculator:
    def __init__(self):
        self.cumulativeData = dict((a,0) for a in range(2,10))
        self.currentData = dict(self.cumulativeData)
    def collect(self):
        """Collect data from simulation"""
        # acumulate cumulative data
        for key,valve in self.currentData.items():
           self.cumulativeData[key] += valve
        # compute new current data
        self.currentData = dict((a,self.getDataFromOneParticle(a)) for a in range(2,10))
    def getDataFromOneParticle(self,a):
        b = O.bodies[a]
        totalHeat = 0
        for i in b.intrs():
            O.interactions.erase(i.id1,0)
            O.interactions.erase(i.id1,1)
            O.interactions.erase(0,i.id2)
            O.interactions.erase(1,i.id2)
            if not i.isReal: continue
            penetrationDepth = i.geom.penetrationDepth                     
            radius1 = i.geom.refR1                                         
            radius2 = i.geom.refR2                                         
            effectiveRadius = (radius1*radius2)/(radius1+radius2)          
            contactRadius = math.sqrt(penetrationDepth*effectiveRadius)    
            contactArea = np.pi*(contactRadius**2)                         
            relativeVelocity = i.geom.shearInc.norm()/O.dt                 
            generatedHeat = (i.phys.shearForce.norm()/contactArea)/(i.phys.normalForce.norm()/contactArea) * i.phys.normalForce.norm() * relativeVelocity   
            totalHeat = totalHeat + generatedHeat
        return totalHeat
    def pprint(self):
        print()
        print("current step")
        print("============")
        for key,valve in self.currentData.items():
            print(key,valve)
        print("cumulative")
        print("==========")
        for key,valve in self.cumulativeData.items():
            print(key,valve)
    def getCurrent(self,b):
        return self.currentData[b.id]
    def getCumulative(self,b):
        return self.cumulativeData[b.id]

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()]
 )
 ,GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=1,timestepSafetyCoefficient=0.8,defaultDt=utils.PWaveTimeStep())
 ,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)
 ,PyRunner(iterPeriod=1,command="someCalculator.collect()",label='collect',dead=True)
 ,PyRunner(iterPeriod=1,command="someCalculator.pprint()",label='pprint',dead=True)
 ,PyRunner(iterPeriod=1,command="vtkExport()",label='vtkExp',dead=True)
 ]

someCalculator = SomeCalculator()
 
def triaxDone():
 triax.dead=True
 O.pause()
 
vtk = export.VTKExporter("heatData")
def vtkExport():
    vtk.exportSpheres(what=dict(current="someCalculator.getCurrent(b)",cumulative="someCalculator.getCumulative(b)"))
    
builtins.someCalculator = someCalculator

O.run(1000000,1)
O.step()

#### Applying normal stress
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.01,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]
   print ('Normal stress =',currentSN,' | unbF=',unbF)
   for i in O.interactions.withBody(O.bodies[0].id):
    stiff+=i.phys.kn
   print ('DONE! iter=',O.iter)
   O.pause()
 if stage==1:
  fnDesired=SN*(O.cell.hSize[0,0]*O.cell.hSize[2,2])
  boundaryVel=copysign(min(0.1,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[:4]+[PyRunner(command='servo()',iterPeriod=10,label='servo')]+O.engines[4:]

O.run(1000000,1)

vtkExp.dead=False
collect.dead=False
pprint.dead=False
newton.damping=DAMPSHEAR
fixVel.dead = False

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
        
O.run(100,1)
###

Thank you for taking your time!

Cheers,
Chien-Cheng

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