Want to connect all the particles in an aggregate using cohesive bond

Hi Varun,
I see the main problem in how to correctly create cohesive interactions
between particles with CohFrictMat. Unfortunately, I am not familiar with
this model. You have to wait for somebody experienced with the model :-) or
you can search CohFrictMat in previous questions [1]..
cheers
Jan

[1] https://answers.launchpad.net/yade

2016-05-04 23:42 GMT+02:00 VG <email address hidden>:

> New question #293295 on Yade:
> https://answers.launchpad.net/yade/+question/293295
>
> Hello everyone,
>
> Here is a description of the problem I am trying to solve:
>
> -- I have a granular sample with spherical particles of diameters D1, D2,
> D3 etc. All of such big particles (lets call them agglomerates) are
> composed of smaller particles of diameter 'd', since d is the smallest
> particle size of interest.
>
> -- Initially, within each big spherical particle, all the smaller
> particles should be densely packed and connected with cohesive bonds
> (CohFrictMat). As the loading is applied, the initially defined cohesive
> bonds within each big sphere will break progressively leading to smaller
> fragments.
>
> Now, the problem I am facing: Initially, within a single big sphere
> itself, sometimes there are particles which are shown not to have any
> cohesive interaction, leaving the small individual particle free to fly.
> How can I alleviate this problem ?
>
> This is a follow-up question to
> https://answers.launchpad.net/yade/+question/292846
>
> Here is the script:
>
>
> from yade import pack
>
>
> #############################################################################
> # Set up run
>
> #############################################################################
> run_name="PeriodicCohesive_MWE"
> data_root_dir="."
>
>
>
>
> #############################################################################
> # Materials
>
> #############################################################################
> plate_material=FrictMat(
> young=200e9
> ,poisson=0.3
> ,density=8000.
> ,frictionAngle=radians(30)
> ,label='plate_mat')
>
> O.materials.append(plate_material)
>
>
> sample_material=CohFrictMat(
> young=4e9
> ,poisson=0.25
> ,density=1400
> ,frictionAngle=radians(30)
> ,normalCohesion=1e8*1.2
> ,shearCohesion=.4e8*1.2
> ,momentRotationLaw=True
> ,label='sample_mat')
> O.materials.append(sample_material)
>
>
>
>
> #############################################################################
> # Component dimensions and operating condition
>
> #############################################################################
> # Granular material dimension
> sample_diameter=2e-4
> #sample_diameter=126e-6
> sample_radius=sample_diameter/2.0
> # Sub-particle dimension
> particle_diameter=74e-6
> particle_radius=particle_diameter/2.
>
>
> #############################################################################
> # grinding plate dimension
>
> #############################################################################
>
> rotvel=2./3.*pi*(1.5+0.5)*.254
>
>
> #############################################################################
> # Periodic Geometry
>
> #############################################################################
>
>
> # Set up periodic boundary conditions
> O.periodic=True
> xExt=4*sample_diameter
> yExt=3.*sample_diameter*2 #to block the periodicity in y direction
> zExt=xExt
> xLim=xExt
> yLim=yExt/4
> zLim=zExt
> O.cell.hSize=Matrix3(
> xExt, 0, 0,
> 0, yExt, 0,
> 0, 0, zExt)
>
>
> length=xExt
> height=yExt
> width=zExt
>
> # Top and bottom plate thickness
> thickness=0.1*height
>
>
>
>
>
> bottomBoxes = []
> for ix in (0,1,2):
> for iz in (0,1,2):
> bottomBoxes.append(box( # create 3x3 boxes with 1/3 cell
> size
> center=(xExt/6.*(1+2*ix),yLim -
> thickness/2.0,zExt/6.*(1+2*iz))
> ,extents=(xExt/6.,thickness/2.0,zExt/6.)
> ,wire=False
> ,material='plate_mat'
> ))
>
> bottom_id,bottom_ids = O.bodies.appendClumped(bottomBoxes) # bottom_id is
> the clump id,
>
>
> O.bodies[bottom_id].state.blockedDOFs='xyzXYZ'
>
>
>
>
> #############################################################################
> # Particle Packing
>
> #############################################################################
>
> min_corner= (0,yLim,0)
> max_corner= (xLim, yExt-yLim, zLim)
>
> sp=pack.SpherePack()
> sp.makeCloud( min_corner,max_corner, rMean=sample_radius, periodic=False,
> seed=1)
>
> print "Generated ",len(sp)," aggregates"
>
> ###########################################
> # Sample
> ###########################################
> for s in sp:
> sphere=pack.inSphere((s[0][0],s[0][1],s[0][2]),s[1])
> sp1=pack.randomDensePack(
> sphere
> ,spheresInCell=2000
> ,radius=particle_radius
> ,memoizeDb='/tmp/triaxPackCache.sqlite'
> ,returnSpherePack=True
> )
>
> sp1.toSimulation(material='sample_mat',color=(0.9,0.8,0.6))
> print 'Generated ',len(sp1),' particles'
>
> Gl1_Sphere(stripes=True)
>
> #########
> # Top
> #########
> topBoxes = []
> for ix in (0,1,2):
> for iz in (0,1,2):
> topBoxes.append(box( # create 3x3 boxes with 1/3 cell size
> center=(xExt/6.*(1+2*ix),yExt - yLim +
> thickness/2.0,zExt/6.*(1+2*iz))
> ,extents=(xExt/6.,thickness/2.0,zExt/6.)
> ,wire=False
> ,material='plate_mat'
> ))
>
> top_id,top_ids = O.bodies.appendClumped(topBoxes) # top_id is the clump id,
>
>
> O.bodies[top_id].state.blockedDOFs='xzXYZ'
>
>
>
> plate_downforce=-0.036
>
> O.forces.addF(top_id,(0,plate_downforce,0),permanent=True)
>
> O.bodies[top_id].state.vel[0]= rotvel
>
>
>
> #############################################################################
> # Run the simulation
>
> #############################################################################
> r_int = 1.0 # interaction radius
>
> O.dt=0.5*PWaveTimeStep()
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([
> Bo1_Sphere_Aabb(aabbEnlargeFactor=r_int, label='bo1s')
> ,Bo1_Box_Aabb()
> ], allowBiggerThanPeriod=True
> ),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=r_int,
> label='ig2ss')
> ,Ig2_Box_Sphere_ScGeom6D()
> ],
> [Ip2_FrictMat_FrictMat_FrictPhys()
>
> ,Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp")
> ],
> [Law2_ScGeom_FrictPhys_CundallStrack()
> ,Law2_ScGeom6D_CohFrictPhys_CohesionMoment()
> ]
> ), # End InteractionLoop
> NewtonIntegrator(damping=0.8,gravity=(0.,0.,0.)),
>
> ]
>
> O.step()
>
> def addBodyToAggreg(body,aggreg): # auxiliary function, add body
> [yade.Body instance] and all its neighbors into aggreg (python set instance)
> if body.id in aggreg: # do nothing if b is already in aggreg ...
> return
> aggreg.add(body.id) # ... otherwise add it to aggreg
> intrs = body.intrs()
> for i in intrs: # and add also all its neighbors ...
> if not isinstance(i.phys,CohFrictPhys): # ... but only
> that connected with CohFrictPhys interactions
> continue
> if i.phys.cohesionBroken: # ... but only not yet broken
> continue
> # if not i.isReal:
> # continue
> i2 = i.id1 if i.id2==body.id else i.id2 # choose the
> other body of interaction
> b2 = O.bodies[i2]
> addBodyToAggreg(b2,aggreg) # and add it and all its
> neighbors to aggreg
>
> def aggregs(): # actual function to detect standalone aggregates
> ids = set(b.id for b in O.bodies if isinstance(b.shape,Sphere)) #
> first make a set of all spheres ids
> ret = []
> while len(ids)>0: # while there are still some particles not
> assigned to any aggregate ...
> i = ids.pop() # ... choose one random ...
> b = O.bodies[i]
> a = set() # ... create new aggregate (set of sphere ids)
> addBodyToAggreg(b,a) # ... and add the sphere together
> with all its neigbors to aggregate
> for bid in a: # delete all used ids from ids
> ids.discard(bid)
> ret.append(a)
> return ret
>
> aggs = aggregs()
> for a in aggs:
> print a
>
>
> I would really appreciate your help.
>
> Thanks
> Varun
>
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Thanks Jan!
 Thanks to Bruno, I actually found part of the solution here: https://answers.launchpad.net/yade/+question/256364
 However, there is something additional described below which I am trying to find an answer for.

There are two things which I am doing now in addition to the above script:
 1. Using a value of 1.1 for aabbEnlargeFactor and interactionDetectionFactor for collision and interaction respectively.
 2. Using very high normal and shear cohesion values in the first iteration: ~1e20 and then resetting these values to the original values.

This results in the correct number of cohesive interactions in each aggregate. However, with the application of the load, these cohesive bonds don't seem to break anymore.

I am using the commands in the following sequence:

-- In the initial definition of material properties for sample material, set very high values of cohesion ~1e20
-- Run one iteration using O.step()
-- Reset the value of cohesion parameters by:
       sample_material.normalCohesion = 1e8*1.2
       sample_material.shearCohesion = 0.4*1e8*1.2
-- O.run()

I also checked the value of sample_material.normalCohesion from command line interface, and it shows the correct value.

Is there a problem in the sequence of commands ? I think I am pretty close to finding a solution but not there yet.
May be when I am setting very high values of normal and shear cohesion initially, some parameters are calculated internally and those need to be reset as well after the first iteration ?

Hi varun,

I did not go through your script. But from your description in #3, this may be due to the value of cohesion.
The value of shearCohesion/normalCohesion within the existing interactions will not change if you just change the value. You need another setCohesionNow to reset it.
Since the very high cohesion value result in the correct number of cohesive interactions, you may want to find a proper value that make the initial interaction number correct and these bonds could be broken in the following process.

Cheers,
Huanran

Hello Huanran,

Thanks for your response, that makes sense.

So it seems like I can not reset the value of shearCohesion/normalCohesion for existing interactions created with higher initial values of these parameters.
If I do another setCohesionNow with lower (original) values of cohesion parameters, then it would not lead to the correct number of cohesive interactions that I want.

Also, I would prefer to keep the cohesion parameters at their original value after the initial interactions, since otherwise a different value will alter the breakage behavior.

Thanks
Varun

Hello Varun,

If you do another setCohesionNow, the cohesion for the existing interactions will change after another timestep. You can write a minimal script to verify it or check that in the source code.

If you prefer to keep to keep the cohesion parameters at their original value (this make sense), you'd better choose a proper value, not too high (bonds wouldn't break) or too low (seems no bonds at all).

Cheers,
Huanran

Thanks Huanran, I will look into it further on the lines of what you suggested.

Thanks Bruno, that solved the issue I was facing.

Also, can you please recommend me a good reference for the CohFrictMat and cohesive model ?

Thanks Bruno Chareyre, that solved my question.

Regarding the above solution, it should be

i.phys.unp = i.geom.penetrationDepth

instead of

i.phys.unp = 0

#to make it initially force free