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.

../../../../../../../_images/p3d-examples-ribbonblender-geometry.png

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).

../../../../../../../_images/BlenderBalls.png

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.

../../../../../../../_images/BlenderBallsResults.png

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.

../../../../../../../_images/BlenderClumps.png

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.

../../../../../../../_images/BlenderClumpsResults.png

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.

../../../../../../../_images/ClusterModel.png

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.

../../../../../../../_images/SizeDistribution.png

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.

../../../../../../../_images/BlenderClusters.png

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 '../blender_geometry/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 prop dp_nratio 0.8 dp_mode 1

wall import from-geom 'rotor1' clean id 1
wall import from-geom 'rotor2' clean id 2
wall import from-geom 'rotor3' clean id 3
wall import from-geom 'rotor4' clean id 4
wall import from-geom 'rotor5' clean id 5
wall import from-geom 'rotor6' clean id 6
wall import from-geom 'rotor7' clean id 7
wall import from-geom 'rotor8' clean id 8

wall import from-geom 'blade1' clean id 9
wall import from-geom 'blade2' clean id 10
wall import from-geom 'blade3' clean id 11
wall import from-geom 'blade4' clean id 12
wall import from-geom 'blade5' clean id 13
wall import from-geom 'blade6' clean id 14
wall import from-geom 'blade7' clean id 15
wall import from-geom 'blade8' clean id 16
wall group 'facetassembly' facet
wall group 'wallassembly'

wall import from-geom '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 def 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 def 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 def 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 prop 'kn' [knB_] 'ks' [ksB_] 'fric' [fric_]
wall prop 'kn' [knW_] 'ks' [ksW_] 'fric' [fric_*0.5]
ball attr density [dens_] damp 0.7
model cycle 2000 calm 500
model mechanical timestep scale
model solve ratio-average 1e-3
model mechanical timestep auto

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 def 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 '../blender_geometry/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-geom 'rotor1' clean id 1
wall import from-geom 'rotor2' clean id 2
wall import from-geom 'rotor3' clean id 3
wall import from-geom 'rotor4' clean id 4
wall import from-geom 'rotor5' clean id 5
wall import from-geom 'rotor6' clean id 6
wall import from-geom 'rotor7' clean id 7
wall import from-geom 'rotor8' clean id 8

wall import from-geom 'blade1' clean id 9
wall import from-geom 'blade2' clean id 10
wall import from-geom 'blade3' clean id 11
wall import from-geom 'blade4' clean id 12
wall import from-geom 'blade5' clean id 13
wall import from-geom 'blade6' clean id 14
wall import from-geom 'blade7' clean id 15
wall import from-geom 'blade8' clean id 16
wall group 'facetassembly' facet 
wall group 'wallassembly'

wall import from-geom '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 def 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 inh off ...  
                     ks 10e4 inh off fric 0.18
contact cmat default type Pebble-Facet ...
                     model Linear property kn 10e5 inh off ...
                     ks 10e5 inh 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'       ...
                        az 90 90          ...
                        el 0 0            ...
                        tilt 0 0          ...
                 box [lbx] [ubx] [lby] [uby] [lbz] [lbz+0.3*dz]
 
fish def 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 def 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 attr density [dens_] damp 0.7
model cycle 2000 calm 500
model mechanical timestep scale
model solve ratio-average 1e-1
model mechanical timestep auto
clump attr damp 0.0

wall attr rotation-center [rp_] range group 'wallassembly'
wall attr 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 '../blender_geometry/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-geom 'rotor1' clean id 1
wall import from-geom 'rotor2' clean id 2
wall import from-geom 'rotor3' clean id 3
wall import from-geom 'rotor4' clean id 4
wall import from-geom 'rotor5' clean id 5
wall import from-geom 'rotor6' clean id 6
wall import from-geom 'rotor7' clean id 7
wall import from-geom 'rotor8' clean id 8

wall import from-geom 'blade1' clean id 9
wall import from-geom 'blade2' clean id 10
wall import from-geom 'blade3' clean id 11
wall import from-geom 'blade4' clean id 12
wall import from-geom 'blade5' clean id 13
wall import from-geom 'blade6' clean id 14
wall import from-geom 'blade7' clean id 15
wall import from-geom 'blade8' clean id 16
wall group 'facetassembly' facet
wall group 'wallassembly'

wall import from-geom '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 inh off ...  
                     ks 10e4 inh off fric 0.18
contact cmat default type pebble-facet ...
                     model linear property kn 10e5 inh off ...
                     ks 10e5 inh 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'     ...
                     az 90 90            ...
                     el 0 0              ...
                     tilt 0 0            ...  
                 box [lbx] [ubx] [lby] [uby] [lbz] [lbz+0.3*dz]

fish def 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 def 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 attr 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 inh off ks 10e4 inh off fric 0.18 inh off
contact cmat 	default type ball-facet   model Linear ... 
        property kn 10e5 inh off ks 10e5 inh off fric 0.09 inh off
contact cmat	add 1 model linearPBond ... 
        property fric 0.0 kn 0.0 inh off ks 0.0 inh off dp_nratio 0.2 ...
        range fish inclump 

model clean
contact cmat apply

[applybonds]

fish def 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 attr density 2000 damp 0.7

fish callback add restorefriction event bond_break

model mechanical timestep auto
model cycle 2000 calm 500
ball attr damp 0.0

wall attribute rotation-center [rp_] range group 'wallassembly'
wall attribute spin [spin_] range group 'wallassembly'

fish def initializeTime
  mp_age_c = mech.time.total
  global tinc = 0.01
  global tnext = tinc
end
[initializeTime]

fish def 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] rad [prad] pos [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 def 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 def 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 def 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

Endnotes

[1]

To view this project in PFC, use the program menu.

Help ▼ Examples…

  ⮡   PFC
    ⮡   Examples
      ⮡   RibbonBlender/Balls/Clumps/Clusters
        ⮡   Blender   _balls.p3prj/_clumps.p3prj/_clusters.p3prj