yade-users team mailing list archive

[Xin Wu <question226198@xxxxxxxxxxxxxxxxxxxxx>]

New question #226198 on Yade:
https://answers.launchpad.net/yade/+question/226198

Dear Sir,
I have now some Problems with Clump. I will try Box-Experiment with Clump and write the Code. Butthe Prigramme stop wth warn Clump::updateProperties: NaNs in eigen-decomposition of inertia matrix?!. Can anyone help me. Thank you.

from yade import utils
from numpy import linspace
from numpy import arange
import gts
import itertools
from yade import pack

rMean=.0096
rRelFuzz=.0016
maxLoad=4500

frictionAngleSt=radians(35)
frictionAngleBo=radians(23.5)

a=0.05
id_Mat=O.materials.append(FrictMat(density=2650,young=175e6,poisson=0.3,frictionAngle=frictionAngleBo))
Mat=O.materials[id_Mat]
steel=O.materials.append(FrictMat(young=210e9, poisson=.25,frictionAngle=frictionAngleSt,density=8000))

from yade import pack, plot,qt

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(0.3,0.3,0.3250),rMean=.0096,rRelFuzz=.0016,periodic=False)
O.bodies.append([sphere(a,r,material=Mat) for a,r in sp])

print len(sp),' particles generated.'

def makeClumpTemplate():
    relRadList1 = [.007,.0036]
    relPosList1 = [[0,0,0],[0.009,0,0]]
    relRadList2 = [.0036,0.007,.0036]
    relPosList2 = [[0.0060,0,0],[0,0,0],[-0.0060,0,0]]
    templates= []
    templates.append(clumpTemplate(relRadii=relRadList1,relPositions=relPosList1))
    templates.append(clumpTemplate(relRadii=relRadList2,relPositions=relPosList2))
    O.bodies.replaceByClumps(templates,[.2,.3])


fIDSI=O.bodies.append(utils.geom.facetBox((.15,.15,.135),(.15,.15,.045),wallMask=15,material=steel))
fIDSII=O.bodies.append(utils.geom.facetBox((.15,.15,.045),(.15,.15,.045),wallMask=31,material=steel))

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   NewtonIntegrator(damping=0.7,gravity=[0,0,-9.810]),
   qt.SnapshotEngine(iterPeriod=60000,fileBase='/home/wuxin/Desktop/Clump/78Fr07-',label='snapshooter'),
   PyRunner(command='checkUnbalanced()',realPeriod=1,label='checker'),
   TranslationEngine(translationAxis=[1,0,0],velocity=0,ids=fIDSI,label='Transl'),
]
O.dt=.25*utils.PWaveTimeStep()

def AngVel():
    for b in O.bodies:
        if isinstance(b.shape,Sphere):
           b.state.angVel[1]=0          

def checkUnbalanced():
    if O.iter<240000: return 
    fIDSIII=O.bodies.append(utils.wall(max([b.state.pos[2]+b.shape.radius for b in O.bodies if isinstance(b.shape,Sphere)]),axis=2,sense=-1,material=spheremat))
    global plate        
    plate=O.bodies[fIDSIII]  
    plate.state.vel=(0,0,-.025)
    O.engines=O.engines+[PyRunner(command='addPlotData()',iterPeriod=1000)]
    checker.command='unloadPlate()'
    checker.command='makeClumpTemplate()' 

def unloadPlate():
   if abs(O.forces.f(plate.id)[2])>=maxLoad:
      O.engines=O.engines+[PyRunner(command='ServoContorl()',iterPeriod=1)]
      O.engines=O.engines+[PyRunner(command='AngVel()',iterPeriod=1)]
      checker.command='stopUnloading()'

def stopUnloading():
    if abs(O.forces.f(plate.id)[2])>(maxLoad-1):
       if abs(O.forces.f(plate.id)[2])<(maxLoad+1):
          Transl.velocity=0.0010

def ServoContorl():
    Fzs=O.forces.f(plate.id)[2]
    if Fzs==0 :
       Fzs=0.000000000000001
    global G
    global KN
    KN=0.000000000000000001
    for i in plate.intrs() :
        KN=KN+i.phys.kn
    G=a/(KN*(O.dt))
    plate.state.vel[2]=(G*(Fzs-maxLoad))
 
def addPlotData():
    Fz=O.forces.f(plate.id)[2]
    F = 0
    global i
    for i in fIDSI:
        F += O.forces.f(i)[0]
        PX1=(O.bodies[i].state.pos[0]-O.bodies[i].state.refPos[0])
        SF=((0.30-PX1)*0.30)
    plot.addData(t=O.time,Fx=(-F/SF),PX=PX1,Fz=Fz,w=plate.state.pos[2]-plate.state.refPos[2],i=O.iter)
    if (O.bodies[i].state.pos[0]-O.bodies[i].state.refPos[0])>0.050:
       plot.saveDataTxt('175e6p03RW235.txt',vars=('i','PX','Fx','w','Fz'))
       O.pause()

plot.plots={'i':('w','Fz',),'PX':('Fx',)}
plot.plot()

qt.Controller()
qt.View()
r=qt.Renderer()

O.run()

utils.waitIfBatch()

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.