Ribbon Blender
Problem Statement
Note
The project file for this example may be viewed/run in PFC.[1] The data files used are shown at the end of this example.
Ribbon blenders are frequently used in manufacturing and industrial applications. The mixer has a central shaft with mixing blades angled in different ways. These look like ribbons of metal wrapped around the shaft. The geometry is able to move parts of the mixture in different directions at the same time, ensuring that all ingredients are blended.
The discrete element method is well-suited for this kind of application. A ribbon blender is modeled using PFC3D, and simulations are performed using three kinds of particles: balls, clumps, and clusters of balls (the latter allows simulation of particle breaking during mixing).
Balls Simulation
The geometry import command is used to import the geometry of the container, the shaft and the mixing blades. This command allows the direct use of CAD files (.dxf and .stl formats are supported) to include complex-shaped elements in the model that will interact with the particles. The file “importBlenderGeometry.dat” shows the use of the command.
Faceted walls are generated for each imported geometry, using the wall import from-geometry command. The blender structure issued from these operations is shown in Figure 1.
Figure 1: Modeling a ribbon blender.
An assembly of balls is then generated (see the ball distribute command) within the boundaries defined by the wall positions. Due to the complex geometry of the blender structure, particles are primarily generated within a box that exceeds the model boundaries in some zone. Then the position of each ball is revisited and the particle itself deleted whenever it is found to occupy a zone outside the model limits. This is done with the functions inBlender and inBody, which control whether each particle stays outside the model boundaries or whether it intersects the blender blades. The FISH function wall.inside and the range FISH commands are used to accomplish these operations. The former checks whether a particle falls within a wall or not; the latter allows automated detection of the particles to be deleted, according
to the criteria specified in the functions inBlender and inBody.
Finally, particles are allowed to settle inside the blender and the simulation is ready to start. A rotational velocity is defined for the walls and four sets of balls are grouped according to their initial positions in order to evaluate the mixing efficiency and effects during the simulation (see Figure 2).
Figure 2: Modeling a ribbon blender to mix spherical particles.
All of the commands are listed in the file “BlenderBalls.dat”.
The mixing efficiency can be evaluated from a qualitative point of view. After a number of revolutions of the blades, we observe the result of mixing the four initial groups. A screenshot of the simulation after a few revolutions of the blades is shown in Figure 3.
Figure 3: Side view of the blender after a few steps.
Clumps Simulation
Similar simulations can be done with non-spherical particles, using the clump logic. In particular, in this example, the geometry import command will be used to create a clump template. Thus, the “ellipsoid.stl” file is imported and used to define the clump template geometry, as it is shown in “BlenderClumps.dat”.
The geometry of the blender is imported using the same procedure that was used for the simulation with the spherical particles (see previous section). A clump template is created, and an assembly of clumps is generated within the model boundaries. Clumps are made of pebbles that are rigidly connected, giving the final shape to the particle. Once the cloud of clumps has been created, the inBlender and inBody functions control whether any pebble stays outside of the model boundaries or whether they intersect the blender blades. In such a case, the entire clump is deleted (the range FISH command is used again). Finally, clumps settle inside the blender and the simulation is ready to start.
Figure 4: Modeling a ribbon blender to mix elliptic-shaped particles.
A screenshot of the simulation after a few revolutions of the blades is shown in Figure 5.
Figure 5: Side view of the blender after few steps.
Clusters Simulation
In order to simulate granulation, the clumps are replaced by clusters of balls. The pebbles are bonded together with parallel bonds. Such parallel bonds are characterized by user-defined tensile strength and cohesion, allowing for breakage during the simulation.
Figure 6: Cluster of balls.
Clumps are generated and allowed to settle, as shown in the previous section. Subsequently, the clumps are replaced by clusters of balls. Each operation is listed in “ClusterMakeCluster.dat”. The function replace goes through all of the existing clumps in the model and calls the function make_cluster, which contains the operations needed to substitute each pebble for a ball. Lists of balls comprising each cluster are created. These lists are then ordered to ensure that two consecutive balls in the list are neighbors in the cluster (see the order_clusters function). This allows one to easily install parallel bonds between neighboring balls, avoiding any accidental bonding of balls that overlap even if they are not neighbors, as shown in Figure 6.
The function apply_bonds is invoked to install the parallel bonds whenever two balls are found to be part of the same cluster and are neighbors. Normal and shear stiffness of the bond, together with tensile and shear strength force, are defined. A linear model with friction is set at contacts between clusters. If a bond breaks during the simulation, the contact model between formerly bonded balls is switched to the linear type. A linearly increasing stiffness (see the lin_inc keyword of the linear model, in the “Linear Model Properties” section) is defined to avoid the sudden repulsive forces that could occur after breakage. A callback event bond_change is called whenever a bond breaks (see the “Linear Parallel Bond Model” section); the function restore_friction is called at the same moment, and the contact model is updated.
All of the operations for this new simulation are listed in “BlenderClusters.dat”.
Periodically during the simulation, the size_distribution function is called to compute the size distribution curve. This allows tracking the evolution of the system and breaking clusters into smaller particles. The size_distribution function is listed in “ClusterSizeDistribution.dat”.
It is useful to look at the evolution of the grain size distribution. The evolution of the grain size distribution curve is shown in Figure 7. Note that the diameters plotted are the equivalent radius of the sphere that has the same volume as the total volume of the cluster.
Figure 7: Cluster simulation. Evolution of the grain size distribution curve.
A screenshot of the simulation after a few revolutions of the blades is shown in Figure 8. Particles are colored by their velocity. The strongest collisions and the breaking of particles seem to occur predominantly in the central zone of the blender.
Figure 8: Cluster simulation, side view of the blender. Particles velocity.
It would be interesting to build a movie of such a simulation. A series of images can be generated easily using the ball results or clump results command. See the “Hopper Discharge” tutorial problem to see an example of its use.
Data Files
BlenderBalls.dat
; fname: BlenderBalls.dat
model new
model large-strain on
wall resolution full
program call 'ImportBlenderGeometry'
model domain extent [lbx-0.05*dx] [ubx+0.05*dx] [lby-0.05*dy] ...
[uby+0.05*dy] [lbz-0.05*dz] [ubz+dz] condition destroy
contact cmat default model linear property dp_nratio 0.8 dp_mode 1
wall import from-geometry 'rotor1' clean id 1
wall import from-geometry 'rotor2' clean id 2
wall import from-geometry 'rotor3' clean id 3
wall import from-geometry 'rotor4' clean id 4
wall import from-geometry 'rotor5' clean id 5
wall import from-geometry 'rotor6' clean id 6
wall import from-geometry 'rotor7' clean id 7
wall import from-geometry 'rotor8' clean id 8
wall import from-geometry 'blade1' clean id 9
wall import from-geometry 'blade2' clean id 10
wall import from-geometry 'blade3' clean id 11
wall import from-geometry 'blade4' clean id 12
wall import from-geometry 'blade5' clean id 13
wall import from-geometry 'blade6' clean id 14
wall import from-geometry 'blade7' clean id 15
wall import from-geometry 'blade8' clean id 16
wall group 'facetassembly' facet
wall group 'wallassembly'
wall import from-geometry 'bin' clean id 17
wall activeside top range position-z [ubz-0.05*ubz] [ubz+0.05*ubz]
wall activeside top range name 'facetassembly'
model save 'blender'
fish define setup(rpm)
global rp_ = vector(0,0,0) ;rotation point
global spin_ = vector(0,-1,0)*2*math.pi*rpm/60.0 ;axis of rotation
global knB_ = 1e4
global ksB_ = 1e4
global knW_ = 1e6
global ksW_ = 1e6
global fric_ = 0.18 ;friction
global dens_ = 2000 ;density
end
[setup(60)]
;generate some balls in the box
fish define inBlender(pos,b)
inBlender = 0
if type.pointer.id(b) = ball.typeid
wp = wall.find(17);
inBlender = wall.inside(wp,ball.pos(b))
endif
end
fish define inBody(pos,b)
inBody = 0
if type.pointer.id(b) = ball.typeid
local tmp = 0
loop foreach local wp wall.list
if wall.id(wp) # 17
tmp = wall.inside(wp,ball.pos(b))
if tmp = 1
inBody = 1
exit
endif
endif
endloop
endif
end
ball distribute porosity 0.45 ...
resolution 1.5e-2 ...
bin 1 ...
radius 1.0 1.5 ...
volume-fraction 1.0 ...
box [lbx] [ubx] [lby] [uby] [lbz] [lbz+0.3*dz] ...
range fish inBlender ...
fish inBody not
model gravity 0,0,-9.81
ball property 'kn' [knB_] 'ks' [ksB_] 'fric' [fric_]
wall property 'kn' [knW_] 'ks' [ksW_] 'fric' [fric_*0.5]
ball attribute density [dens_] damp 0.7
model cycle 2000 calm 500
model mechanical timestep scale
model solve ratio-average 1e-3
model mechanical timestep automatic
wall attribute rotation-center [rp_] range group 'wallassembly'
wall attribute spin [spin_] range group 'wallassembly'
;set up groups for mixing
ball group 'FrontLeft' range position-x [lbx-0.05*dx] 0.0 ...
position-y [lby-0.05*dy] 0.5
ball group 'FrontRight' range position-x 0.0 [ubx+0.05*dx] ...
position-y [lby-0.05*dy] 0.5
ball group 'BackLeft' range position-x [lbx-0.05*dx] 0.0 ...
position-y 0.5 [uby+0.05*dy]
ball group 'BackRight' range position-x 0.0 [ubx+0.05*dx] ...
position-y 0.5 [uby+0.05*dy]
model save 'ball_model'
model solve time 1.0
model save 'end_blenderball'
;=========================================================================
; eof: BlenderBalls.dat
importBlenderGeometry.dat
; fname: ImportBlenderGeometry.dat
geometry import 'rotor1.stl'
geometry import 'rotor2.stl'
geometry import 'rotor3.stl'
geometry import 'rotor4.stl'
geometry import 'rotor5.stl'
geometry import 'rotor6.stl'
geometry import 'rotor7.stl'
geometry import 'rotor8.stl'
geometry import 'blade1.stl'
geometry import 'blade2.stl'
geometry import 'blade3.stl'
geometry import 'blade4.stl'
geometry import 'blade5.stl'
geometry import 'blade6.stl'
geometry import 'blade7.stl'
geometry import 'blade8.stl'
geometry import 'bin.stl'
;calculate the dimensions of the box needed
fish define cdim
lbx = 1e5
lby = 1e5
lbz = 1e5
ubx = -1e5
uby = -1e5
ubz = -1e5
local gs = geom.set.find(17)
loop foreach local gn geom.node.list(gs)
local pos = geom.node.pos(gn)
if lbx > comp.x(pos)
lbx = comp.x(pos)
endif
if lby > comp.y(pos)
lby = comp.y(pos)
endif
if lbz > comp.z(pos)
lbz = comp.z(pos)
endif
if ubx < comp.x(pos)
ubx = comp.x(pos)
endif
if uby < comp.y(pos)
uby = comp.y(pos)
endif
if ubz < comp.z(pos)
ubz = comp.z(pos)
endif
endloop
dx = ubx - lbx
dy = uby - lby
dz = ubz - lbz
end
[cdim]
;=========================================================================
; eof: ImportBlenderGeometry.dat
BlenderClumps.dat
; fname: BlenderClumps.dat
model new
model large-strain on
wall resolution full
program call 'ImportBlenderGeometry'
model domain extent [lbx-0.05*dx] [ubx+0.05*dx] [lby-0.05*dy] ...
[uby+0.05*dy] [lbz-0.05*dz] [ubz+dz] condition destroy
wall import from-geometry 'rotor1' clean id 1
wall import from-geometry 'rotor2' clean id 2
wall import from-geometry 'rotor3' clean id 3
wall import from-geometry 'rotor4' clean id 4
wall import from-geometry 'rotor5' clean id 5
wall import from-geometry 'rotor6' clean id 6
wall import from-geometry 'rotor7' clean id 7
wall import from-geometry 'rotor8' clean id 8
wall import from-geometry 'blade1' clean id 9
wall import from-geometry 'blade2' clean id 10
wall import from-geometry 'blade3' clean id 11
wall import from-geometry 'blade4' clean id 12
wall import from-geometry 'blade5' clean id 13
wall import from-geometry 'blade6' clean id 14
wall import from-geometry 'blade7' clean id 15
wall import from-geometry 'blade8' clean id 16
wall group 'facetassembly' facet
wall group 'wallassembly'
wall import from-geometry 'bin' clean id 17
wall activeside top range position-z [ubz-0.05*ubz] [ubz+0.05*ubz]
wall activeside top range name 'facetassembly'
model save 'blender'
fish define setup(rpm)
rp_ = vector(0,0,0)
spin_ = vector(0,-1,0)*2*math.pi*rpm/60.0
knB_ = 1e4
ksB_ = 1e4
knW_ = 1e8
ksW_ = 1e8
fric_ = 0.18
dens_ = 2000
end
[setup(60)]
geometry import 'ellipsoid.stl'
clump template create name 'pill' geometry 'ellipsoid' ...
calculate-surface bubblepack ratio 0.5 distance 150
contact cmat default type Pebble-Pebble ...
model Linear property kn 10e4 inheritance off ...
ks 10e4 inheritance off fric 0.18
contact cmat default type Pebble-Facet ...
model Linear property kn 10e5 inheritance off ...
ks 10e5 inheritance off fric 0.09
clump distribute porosity 0.45 ...
resolution 5e-2 ...
number-bins 1 ...
diameter ...
bin 1 size 1.0 1.20 ...
volume-fraction 1 ...
template 'pill' ...
azimuth 90 90 ...
elevation 0 0 ...
tilt 0 0 ...
box [lbx] [ubx] [lby] [uby] [lbz] [lbz+0.3*dz]
fish define inBlender(pos,c)
inBlender = 1
if type.pointer.id(c) = clump.typeid
loop foreach local pb clump.pebblelist(c)
local tmp = 0
wp = wall.find(17);
tmp = wall.inside(wp,clump.pebble.pos(pb))
if tmp = 0
inBlender = 0
exit
endif
endloop
endif
end
clump delete range fish inBlender not
fish define inBody(pos,c)
inBody = 0
if type.pointer.id(c) = clump.typeid
loop foreach local pb clump.pebblelist(c)
local tmp = 0
loop foreach local wp wall.list
if wall.id(wp) # 17
tmp = wall.inside(wp,clump.pebble.pos(pb))
if tmp = 1
inBody = 1
exit
endif
endif
endloop
endloop
endif
end
clump delete range fish inBody
model gravity 0,0,-9.81
clump attribute density [dens_] damp 0.7
model cycle 2000 calm 500
model mechanical timestep scale
model solve ratio-average 1e-1
model mechanical timestep automatic
clump attribute damp 0.0
wall attribute rotation-center [rp_] range group 'wallassembly'
wall attribute spin [spin_] range group 'wallassembly'
;set up groups for mixing
clump group 'FrontLeft' range ...
position-x [lbx-0.05*dx] 0.0 position-y [lby-0.05*dy] 0.5
clump group 'FrontRight' range ...
position-x 0.0 [ubx+0.05*dx] position-y [lby-0.05*dy] 0.5
clump group 'BackLeft' range ...
position-x [lbx-0.05*dx] 0.0 position-y 0.5 [uby+0.05*dy]
clump group 'BackRight' range ...
position-x 0.0 [ubx+0.05*dx] position-y 0.5 [uby+0.05*dy]
model save 'clump_model'
model solve time 0.5
model save 'end_blenderclump'
;=========================================================================
; eof: BlenderClumps.dat
BlenderClusters.dat
; fname: BlenderClusters.dat
model new
model large-strain on
wall resolution full
program call 'ClusterMakeCluster'
program call 'ClusterSizeDistribution'
program call 'ImportBlenderGeometry'
geometry import 'ellipsoid.stl'
model domain extent [lbx-0.05*dx] [ubx+0.05*dx] [lby-0.05*dy] ...
[uby+0.05*dy] [lbz-0.05*dz] [ubz+dz] condition destroy
wall import from-geometry 'rotor1' clean id 1
wall import from-geometry 'rotor2' clean id 2
wall import from-geometry 'rotor3' clean id 3
wall import from-geometry 'rotor4' clean id 4
wall import from-geometry 'rotor5' clean id 5
wall import from-geometry 'rotor6' clean id 6
wall import from-geometry 'rotor7' clean id 7
wall import from-geometry 'rotor8' clean id 8
wall import from-geometry 'blade1' clean id 9
wall import from-geometry 'blade2' clean id 10
wall import from-geometry 'blade3' clean id 11
wall import from-geometry 'blade4' clean id 12
wall import from-geometry 'blade5' clean id 13
wall import from-geometry 'blade6' clean id 14
wall import from-geometry 'blade7' clean id 15
wall import from-geometry 'blade8' clean id 16
wall group 'facetassembly' facet
wall group 'wallassembly'
wall import from-geometry 'bin' clean id 17
wall activeside top range position-z [ubz-0.05*ubz] [ubz+0.05*ubz]
wall activeside top range name 'facetassembly'
fish define setup(rpm)
rp_ = vector(0,0,0)
spin_ = vector(0,-1,0)*2*math.pi*rpm/60.0
knB_ = 1e4
ksB_ = 1e4
knW_ = 1e8
ksW_ = 1e8
fric_ = 0.18
dens_ = 2000
end
[setup(60)]
clump template create name 'pill' ...
geometry 'ellipsoid' ...
calculate-surface bubblepack ...
ratio 0.5 ...
distance 150
contact cmat default type pebble-pebble ...
model linear property kn 10e4 inheritance off ...
ks 10e4 inheritance off fric 0.18
contact cmat default type pebble-facet ...
model linear property kn 10e5 inheritance off ...
ks 10e5 inheritance off fric 0.09
clump distribute porosity 0.45 ...
number-bins 1 ...
diameter ...
resolution 1.0 ...
bin 1 size 0.05 0.075 ...
volume-fraction 1 ...
template 'pill' ...
azimuth 90 90 ...
elevation 0 0 ...
tilt 0 0 ...
box [lbx] [ubx] [lby] [uby] [lbz] [lbz+0.3*dz]
fish define inBlender(pos,c)
inBlender = 1
if type.pointer.id(c) = clump.typeid
loop foreach local pb clump.pebblelist(c)
local tmp = 0
wp = wall.find(17)
tmp = wall.inside(wp,clump.pebble.pos(pb))
if tmp = 0
inBlender = 0
exit
endif
endloop
endif
end
clump delete range fish inBlender not
fish define inBody(pos,c)
inBody = 0
if type.pointer.id(c) = clump.typeid
loop foreach local pb clump.pebblelist(c)
local tmp = 0
loop foreach local wp wall.list
if wall.id(wp) # 17
tmp = wall.inside(wp,clump.pebble.pos(pb))
if tmp = 1
inBody = 1
exit
endif
endif
endloop
endloop
endif
end
clump delete range fish inBody
model gravity 0,0,-9.81
clump attribute density [dens_] damp 0.7
model cycle 2000 calm 500
[ini_clusters]
[replace]
[order_clusters]
[list_clusters]
fish define inclump(extra, ct)
inclump = 0
if type.pointer.id(ct) = contact.typeid("ball-ball")
if ball.extra(contact.end1(ct),1) == ball.extra(contact.end2(ct),1) then
inclump = 1
exit
endif
endif
end
contact cmat default type ball-ball model linear ...
property kn 10e4 inheritance off ks 10e4 inheritance off fric 0.18 inheritance off
contact cmat default type ball-facet model Linear ...
property kn 10e5 inheritance off ks 10e5 inheritance off fric 0.09 inheritance off
contact cmat add 1 model linearPBond ...
property fric 0.0 kn 0.0 inheritance off ks 0.0 inheritance off dp_nratio 0.2 ...
range fish inclump
model clean
contact cmat apply
[applybonds]
fish define restorefriction(arr)
ct=arr(1)
;local oa = out(string.build("extra1 %1 - ...
; extra2 %2",ball.extra(c_end1(ct),1), ...
; ball.extra(c_end2(ct),1)))
if ball.extra(contact.end1(ct),1) == ball.extra(contact.end2(ct),1)
io.out(string.build("Contact broken! %1",contact.prop(ct,'pb_state')))
contact.model(ct) = 'Linear'
contact.prop(ct,'kn') = 10e4
contact.prop(ct,'ks') = 10e4
contact.prop(ct,'lin_mode') = 1
idc = idc + 1
blist = array.create(1)
blist(1) = contact.end1(ct)
size = array.size(blist,1)
ind = 0
stop = 0
loop while stop = 0
newbonds = 0
newvisited = 0
loop local i (1,size+ind)
;local oppz = out(string.build("blist(i) %1",blist(i)))
if ball.extra(blist(i),2) = 0 ;isNOTvisited
newvisited = newvisited + 1
loop foreach local con ball.contactmap(blist(i))
if contact.model(con) # "Linear"
if contact.prop(con,'pb_state')=3
if ball.extra(contact.end1(con),2) = 0 ;not visited
if ball.extra(contact.end2(con),2) = 0 ; " "
newbonds = newbonds + 1
;local ff = out(string.build...
; ("ball ids = %1 %2",ball.id...
; (c_end1(con)),ball.id(c_end2(con))))
endif
endif
endif
endif
endloop
ball.extra(blist(i),2) = 1 ;isvisited
ball.extra(blist(i),1) = idc
endif
endloop
;local ops = out(string.build("newbonds %1",newbonds))
if newbonds#0
oldblist = blist
blist = array.create(size+ind+newbonds)
loop local k (1,size+ind)
blist(k) = oldblist(k)
endloop
loop local ii (size+ind+1-newvisited,size+ind)
clist = ball.contactmap(blist(ii))
loop foreach local jj clist
if contact.model(jj) # "linear"
if contact.prop(jj,'pb_state')=3
if ball.extra(contact.end1(jj),2) = 0 ;not visited
ind = ind + 1
blist(size+ind) = contact.end1(jj)
else if ball.extra(contact.end2(jj),2) = 0 ;not visited
ind = ind + 1
blist(size+ind) = contact.end2(jj)
endif
endif
endif
endloop
endloop
else
stop = 1
endif
endloop
endif
loop foreach local bp ball.list
ball.extra(bp,2) = 0 ;not visited
endloop
end
ball attribute density 2000 damp 0.7
fish callback add restorefriction event bond_break
model mechanical timestep automatic
model cycle 2000 calm 500
ball attribute damp 0.0
wall attribute rotation-center [rp_] range group 'wallassembly'
wall attribute spin [spin_] range group 'wallassembly'
fish define initializeTime
mp_age_c = mech.time.total
global tinc = 0.01
global tnext = tinc
end
[initializeTime]
fish define buildPSD
whilestepping
global temps = mech.time.total-mp_age_c
if temps > tnext then
sizedistribution(20)
tnext = tnext + tinc
endif
end
[sizedistribution(20)]
model save 'cluster_model'
model solve time 0.1
model save 'end_blendercluster'
;=========================================================================
; eof: BlenderClusters.dat
ClusterMakeCluster.dat
; fname: ClusterMakeCluster.dat
fish define ini_clusters
global cluster_head = null
global bnumber=0
clt = clump.template.find('pill')
loop foreach p clump.template.pebblelist(clt)
bnumber = bnumber + 1
endloop
io.out(string.build("number of pebbles = %1", bnumber))
end
fish define make_cluster(clump)
; replace a clump with a set of balls
newclus = memory.create(bnumber+3)
memory(newclus) = cluster_head
cluster_head = newclus
global idc = clump.id(clump)
local num=1
memory(newclus+num) = idc
num=num+1
local bfirst
loop foreach local p clump.pebblelist(clump)
global prad = clump.pebble.radius(p)
global ppos = clump.pebble.pos(p)
global idp = clump.pebble.id(p)
command
ball create id [idp] radius [prad] position [ppos]
endcommand
num = num + 1
local bp = ball.find(idp)
if num>3 then
if comp.x(ball.pos(bp)) < comp.x(ball.pos(bfirst)) then
bfirst = bp
endif
else
bfirst = bp
endif
ball.extra(bp,1) = idc
memory(newclus+num) = bp
endloop
memory(newclus+2) = bfirst
clump.delete(clump)
end
fish define replace
loop foreach clump clump.list
make_cluster(clump)
endloop
end
fish define order_clusters
cluster = cluster_head
loop while cluster # null
newcluster = array.create(9)
newcluster(1) = memory(cluster+2)
loop j(2,9)
referenceball = newcluster(j-1)
dist=1000
loop i(1,bnumber)
thisball = memory(cluster+2+i)
newdist = comp.x(ball.pos(thisball))-comp.x(ball.pos(referenceball))
io.out(string.build("newdist = %1",newdist))
if newdist > 0
if newdist < dist
nearestball = thisball
dist = newdist
endif
endif
endloop
newcluster(j) = nearestball
endloop
loop i(1,bnumber)
memory(cluster+i+2) = newcluster(i)
endloop
cluster = memory(cluster)
endloop
end
fish define list_clusters
cluster = cluster_head
loop while cluster # null
local idc = memory(cluster+1)
local bleft = memory(cluster+2)
io.out(string.build("--cluster with ID %1 : ",idc))
io.out(string.build("--left most ball id : %1",ball.id(bleft)))
loop i(1,bnumber)
bp = memory(cluster+i+2)
io.out(string.build(" - ball id %1",ball.id(bp)))
endloop
cluster = memory(cluster)
endloop
end
fish define applybonds
cluster = cluster_head
loop while cluster # null
loop i(1,bnumber-1)
bp = memory(cluster+i+2)
bp_next=memory(cluster+i+3)
local clist = ball.contactmap(bp)
section
loop foreach local con clist
ok=0
if contact.end1(con) = bp_next
ok=1
else if contact.end2(con) = bp_next
ok=1
endif
if ok=1
local arg = array.create(1,2)
arg(1,1) = 'gap'
arg(1,2) = 0.0
contact.method(con,"bond",arg)
contact.prop(con,"pb_kn")=1e12
contact.prop(con,"pb_ks")=1e12
contact.prop(con,"pb_ten")=1e8
contact.prop(con,"pb_coh")=1e6
contact.prop(con,"pb_rmul")=0.8
exit section
endif
endloop
endsection
endloop
cluster = memory(cluster)
endloop
end
;=========================================================================
; eof: ClusterMakeCluster.dat
ClusterSizeDistribution.dat
; fname: ClusterSizeDistribution.dat
fish define clumpvolumes
;local oo45 = out(string.build("nclumps = %1",nclumps))
global volumes = array.create(idc)
loop local i (1,array.size(volumes,1))
volumes(i)=0.0
endloop
loop foreach bp ball.list
volumes(ball.extra(bp,1)) = ...
volumes(ball.extra(bp,1)) + 1.3333*math.pi*ball.radius(bp)^3
endloop
end
fish define sizedistribution(intervals)
command
[clumpvolumes]
endcommand
io.out(string.build("size volumes = %1", array.size(volumes,1)))
nentries = array.size(volumes,1)
radius = array.create(nentries)
minD = 1000
maxD = -1
local nonzeroentries = 0
loop local i (1,array.size(volumes,1))
radius(i) = ((3*volumes(i))/(4*math.pi))^(0.33333)
if radius(i)#0
nonzeroentries = nonzeroentries + 1
minD = math.min(radius(i)*2,minD)
maxD = math.max(radius(i)*2,maxD)
endif
;local oox = out(string.build("current value of radius = %1", radius[i]))
endloop
io.out(string.build("minD = %1 , maxD = %2", minD, maxD))
increment = (maxD-minD)/intervals
passing = array.create(intervals+1)
aveD = array.create(intervals+1)
t = string.build("sizedist%1",nentries)
dist = table.get(t)
loop local j (1,array.size(aveD,1))
aveD(j) = minD + increment*(j-1)
passing(j) = 0.0
endloop
table(t,aveD(1))=passing(1)
loop local p (2,array.size(passing,1))
pass = 0
loop local r (1,array.size(radius,1))
d = radius(r)*2
if d#0
if d<=aveD(p)
;local oo52 = ...
;out(string.build("d = %1 , aveD = %2, pass = %3", d, aveD[p+1],pass))
pass = pass+1
endif
endif
endloop
;local oo23 = ...
;out(string.build("nonzeroentries = %1 , pass = %2", nonzeroentries, pass))
passing(p) = 100*(float(pass)/nonzeroentries)
table(t,aveD(p))=passing(p)
endloop
loop local tt (1,array.size(aveD,1))
io.out(string.build("passing[%1] = %2", tt, passing(tt)))
endloop
end
;==============================================================================
; eof: SizeDistribution.dat
Endnote
| Was this helpful? ... | Itasca Software © 2024, Itasca | Updated: Nov 20, 2024 |
