Welcome to ask.woodem.org. You may post when you login through your GitHub account.

How to insert spheres in two times

0 votes
16 views
asked Jul 27 by alanbds (140 points)
from woo.core import *
from woo.dem import *
from minieigen import *
import woo.models, woo

woo.master.usesApi=10102

S=woo.master.scene=Scene(fields=[DemField(gravity=(0,0,-10))])

# contact model plus material for particles
model=woo.models.ContactModelSelector(name='linear',damping=.4,mats=[FrictMat(young=1e7,ktDivKn=.2,tanPhi=.4)])
# same material for walls as for particles, but without friction
S.lab.wallMat=model.mats[0].deepcopy(tanPhi=0.)
# intial box size
box=AlignedBox3((0,0,0),(2,2,3))
subballast=AlignedBox3((0,0,0),(2,2,1.5))
ballast=AlignedBox3((0,0,1.5),(2,2,3))
#boxGeneration=AlignedBox3((0,0,1.5),(3,3,3))
# use something else than spheres here, if you like
gen=PsdSphereGenerator(psdPts=[(.1,0),(.1,.2),(.2,.9),(.3,1)])
ctrlStep=100

# wall box without the top (will be added later)
S.dem.par.add(Wall.makeBox(box=box,which=(1,1,1,1,1,0),mat=S.lab.wallMat))

S.engines=DemField.minimalEngines(model=model)+[
    # generate initial particles, until massRate cannot be achieved, then turn yourself off (dead)
    BoxInlet(stepPeriod=ctrlStep,box=subballast,generator=gen,materials=[model.mats[0]],massRate=1000,maxMass=0,maxAttempts=1000,atMaxAttempts='dead',label='inlet'),
	BoxInlet(stepPeriod=ctrlStep,box=ballast,generator=gen,materials=[model.mats[0]],massRate=1000,maxMass=0,maxAttempts=1000,atMaxAttempts='dead',label='inlet'),
    # check what's happening periodically
    PyRunner(ctrlStep,'checkProgress(S)')
]

# stage: 0=gravity deposition, 1=oedometric loading, 2=oedometric unloading, 3=box loading
S.lab.stage=0
# avoid warning when we change the value
S.lab._setWritable('stage')
# things that we want to look at during the simualtion
S.plot.plots={'dz':('Fz',)}

S.saveTmp()

# callback function from PyRunner
def checkProgress(S):
    # if the top wall is there, save its position and acting force
    if S.lab.stage>0: S.plot.addData(dz=S.lab.z0-S.lab.topWall.pos[2],Fz=S.lab.topWall.f[2],i=S.step,t=S.time)
    # check progress
    # if inlet is off and unbalanced force is already low, go ahead
    if S.lab.stage==0 and S.lab.inlet.dead and woo.utils.unbalancedForce(S)<0.05:
        S.lab.stage+=1
        # find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
        z0=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
        # place Wall just above the top
        S.lab.topWall=Wall.make(z0,axis=2,sense=-1,mat=S.lab.wallMat)
        # save the initial z as reference
        S.lab.z0=z0
        # assign constant velocity to the wall
        S.lab.topWall.vel=(0,0,-.05)
        # add to the simulation
        S.dem.par.add(S.lab.topWall)
    # when the wall goes down to z=1, make it move upwards again
    elif S.lab.stage==1 and S.lab.topWall.pos[2]<1:
        S.lab.stage+=1
        # reverse the loading sense
        S.lab.topWall.vel*=-1 
    # move until the initial position is reached again     
    elif S.lab.stage==2 and S.lab.topWall.pos[2]>S.lab.z0:
        # everything done, nothing will happen with stage==3
        S.lab.stage+=1
        # stop the wall
        S.lab.topWall.vel=(0,0,0) 
        # remove the wall
        S.dem.par.remove(S.lab.topWall.id)
        # find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
        z0=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])        
        # create the box
        zBox=z0+0.25
        S.lab.Box=woo.triangulated.box(dim=(.5,.5,.5),center=(1,1,zBox),which=(1,1,1,1,1,1),fixed=True,mat=S.lab.wallMat,wire=False)
        # add to the simulation
        S.lab.boxNode=Node(pos=(1,1,zBox),dem=ClumpData(blocked='xyzXYZ'))
        S.dem.par.addClumped(S.lab.Box,centralNode=S.lab.boxNode)
        # assign velocity to the box
        S.lab.boxNode.dem.impose=HarmonicOscillation(freq=100.,amp=.02,dir=(0,0,1))        


So... If you run this code, you'll a box with sphere generated and an axial compression after the forces are balanced.

Well, I need to generate the spheres in two times. But the way I did this code I don't know how to make it. I need to insert the spheres from 0 - 1.5 (subballast) and them make a compression on it, after that, insert the spheres from 1.5 - 3.0 (ballast) and them make a compression on it, after that I insert the other box upon the spheres, like it is writen there.

If you guys have any tips for me, please, feel free to help haha

FYI i'm a begginer on programing

commented Aug 17 by eudoxos (44,970 points)

Hi Alan, I was away for a longer time, sorry for the delay. I need some clarification. You need to have 2-layer compact packing, with gravity compaction? Predefined thicknesses of the layer, or it is enough for you to have the ratio between them something given? Are they different material or do they differ in the PSD only? Thanks for clarifications, I will then look for the best way to do it. Cheers! Vaclav

commented Aug 31 by alanbds (140 points)
  1. I need to compact 2-layers with a Wall.make, to represent the compression of the layer;
  2. I have to 2 thicknesses pre-defined, but for now, I just need a code that I can play with thickness;
  3. They are the same material, but with different gradations.

I've done this so far:

from woo.core import *
from woo.dem import *
from minieigen import *
import woo.models, woo

woo.master.usesApi=10102

S=woo.master.scene=Scene(fields=[DemField(gravity=(0,0,-10))])

# contact model plus material for particles
model=woo.models.ContactModelSelector(name='linear',damping=.4,mats=[FrictMat(young=1e7,ktDivKn=.2,tanPhi=.4)])
# same material for walls as for particles, but without friction
S.lab.wallMat=model.mats[0].deepcopy(tanPhi=0.)
# intial box size
box=AlignedBox3((0,0,0),(2,2,4))
lastro=AlignedBox3((0,0,1.5),(2,2,3.5))
slastro=AlignedBox3((0,0,0),(2,2,2))
# use something else than spheres here, if you like
gen=PsdSphereGenerator(psdPts=[(.1,0),(.1,.2),(.2,.9),(.3,1)])
ctrlStep=100

# wall box without the top (will be added later)
S.dem.par.add(Wall.makeBox(box=box,which=(1,1,1,1,1,0),mat=S.lab.wallMat))

S.engines=DemField.minimalEngines(model=model)+[
    # generate initial particles, until massRate cannot be achieved, then turn yourself off (dead)
    BoxInlet(stepPeriod=ctrlStep,box=slastro,generator=gen,materials=[model.mats[0]],massRate=1000,maxMass=0,maxAttempts=1000,atMaxAttempts='dead',label='inletL'),
    BoxInlet(stepPeriod=ctrlStep,box=lastro,generator=gen,materials=[model.mats[0]],massRate=1000,maxMass=0,maxAttempts=1000,atMaxAttempts='dead',label='inletSL',dead=True),
    # check what's happening periodically
    PyRunner(ctrlStep,'checkProgress(S)')
]

# stage: 0=gravity deposition, 1=oedometric loading, 2=oedometric unloading, 3=box loading
S.lab.stage=0
# avoid warning when we change the value
S.lab._setWritable('stage')
# things that we want to look at during the simualtion
S.plot.plots={'dz0':('Fz0',)}
#S.plot.plots={'dz1':('Fz1',)}

S.saveTmp()

# callback function from PyRunner
def checkProgress(S):
    # if the top wall is there, save its position and acting force
    if S.lab.stage>0:
        S.plot.addData(dz0=S.lab.z0-S.lab.topWall.pos[2],Fz0=S.lab.topWall.f[2],i=S.step,t=S.time)
        #S.plot.addData(dz1=S.lab.z1-S.lab.topWall1.pos[2],Fz1=S.lab.topWall1.f[2],i=S.step,t=S.time)
        # nao  consegui achar o erro entre Z0 e Z1 so consegui plotar z0
        # check progress
        # if inlet is off and unbalanced force is already low, go ahead
    if S.lab.stage==0 and S.lab.inletL.dead and woo.utils.unbalancedForce(S)<0.05:
        S.lab.stage+=1
        # find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
        z0=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
        # place Wall just above the top
        S.lab.topWall=Wall.make(z0,axis=2,sense=-1,mat=S.lab.wallMat)
        print('z0' , z0)
        # save the initial z as reference
        S.lab.z0=z0
        # assign constant velocity to the wall
        S.lab.topWall.vel=(0,0,-.05)
        # add to the simulation
        S.dem.par.add(S.lab.topWall)
        # when the wall goes down to z=1, make it move upwards again
    elif S.lab.stage==1 and S.lab.topWall.pos[2]<1:
        S.lab.stage+=1
        # reverse the loading sense
        S.lab.topWall.vel*=-1
        # move until the initial position is reached again
    elif S.lab.stage==2 and S.lab.topWall.pos[2]>S.lab.z0:
        S.lab.stage+=1
        #S.plot.addData(dz1=S.lab.z1-S.lab.topWall1.pos[2],Fz1=S.lab.topWall1.f[2],i=S.step,t=S.time) VERIFICAR VERIFICAR VERIFICAR
        S.lab.topWall.vel=(0,0,0)
        S.dem.par.remove(S.lab.topWall.id)
        S.lab.inletSL.dead=False
    elif S.lab.stage==3 and S.lab.inletSL.dead and woo.utils.unbalancedForce(S)<0.05:
        S.lab.stage+=1
        # find the topmost particle (only spheroids (i.e. Spheres, Ellipsoids, Capsules) have equivRadius)
        z1=max([p.pos[2]+p.shape.equivRadius for p in S.dem.par if p.shape.equivRadius>0])
        print('z1' , z1)
        # place Wall just above the top
        S.lab.topWall1=Wall.make(z1,axis=2,sense=-1,mat=S.lab.wallMat)
        # save the initial z as reference
        S.lab.z1=z1
        # assign constant velocity to the wall
        S.lab.topWall1.vel=(0,0,-.05)
        # add to the simulation
        S.dem.par.add(S.lab.topWall1)
        # when the wall goes down to z=1, make it move upwards again
    elif S.lab.stage==4 and S.lab.topWall1.pos[2]<2:
        S.lab.stage+=1
        # reverse the loading sense
        S.lab.topWall1.vel*=-1
        # move until the initial position is reached again
    elif S.lab.stage==5 and S.lab.topWall1.pos[2]>S.lab.z1:
        # everything done, nothing will happen with stage==3
        S.lab.stage+=1
        S.lab.topWall.vel=(0,0,0)
        S.dem.par.remove(S.lab.topWall1.id)
        # and stop the simulation as well

I need some help to plot the graphic z1: Fz1, from line 39. Didn't get the informations to plot that on the documentation.

commented Aug 31 by eudoxos (44,970 points)

Wow, very nice :) addPlotData always adds a full row of data (filled with NaN for unspecified values), so it should be called only once at the step where you add plot data. Try something like this at the beginning of checkProgress:

if S.lab.stage in (1,2): S.plot.addData(dz=S.lab.z0-S.lab.topWall.pos[2],Fz0=S.lab.topWall.f[2],i=S.step,t=S.time)
elif S.lab.stage in (4,5): S.plot.addData(dz=S.lab.z1-S.lab.topWall1.pos[2],Fz1=S.lab.topWall1.f[2],i=S.step,t=S.time)

For testing, you can set young=5e5, S.dtSafety=.9, goes faster.
wall force plot

As another optimization, this deactivates the wall as soon as it has zero vertical force:

elif S.lab.stage==2 and S.lab.topWall.f[2]==0:
# ...
elif S.lab.stage==5 and S.lab.topWall1.f[2]==0:

1 Answer

0 votes
answered Sep 13 by eudoxos (44,970 points)

For the answer, see the last comment https://ask.woodem.org/index.php/887/how-to-insert-spheres-in-two-times?show=892#c892. You will also have to adjust S.plot.plots=... line to something like {'dz':('Fz0','Fz1')}, of course.

...