yade-users team mailing list archive

[weijy <question404135@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi All,

I am a beginner and have been trying to remodel below paper:

[1] Scholt E S L and Donz E F E D E. 2013. A DEM model for soft and hard rocks: Role of grain interlocking on strength [J]. J. Mech. Phys. Solids.,61(2): 352-369.

In chapter 2.2 interaction range is introduced. I've also read the relavent contents in Yade documentation: 

***********Yade document is posted for reference, you may just ignore and read my question*******************
Some material models (such as the concrete model) rely on initial interaction network which is denser than geometrical
contact of spheres: spheres radii as enlarged by a dimensionless factor called interaction radius (or interaction ratio)
to create this initial network. This is done typically in this way (see examples/concrete/uniax.py for an example):
1. Approximate collision detection is adjusted so that approximate contacts are detected also between particles within
the interaction radius. This consists in setting value of Bo1_Sphere_Aabb.aabbEnlargeFactor to the interaction
radius value.
2. The geometry functor (Ig2) would normally say that there is no contact if given 2 spheres that are not in contact.
Therefore, the same value as for Bo1_Sphere_Aabb.aabbEnlargeFactor must be given to it. (Either Ig2_Sphere_Sphere_Dem3DofGeom.distFactor or Ig2_Sphere_Sphere_ScGeom.interactionDetectionFactor,
depending on the functor that is in use.
Note that only Sphere + Sphere interactions are supported; there is no parameter analogous to distFactor in
Ig2_Facet_Sphere_Dem3DofGeom. This is on purpose, since the interaction radius is meaningful in bulk material
represented by sphere packing, whereas facets usually represent boundary conditions which should be exempt
from this dense interaction network.
3. Run one single step of the simulation so that the initial network is created.
4. Reset interaction radius in both Bo1 and Ig2 functors to their default value again.
5. Continue the simulation; interactions that are already established will not be deleted (the Law2 functor in
usepermitting).
In code, such scenario might look similar to this one (labeling is explained in Labeling things):

intRadius=1.5
O.engines=[
ForceResetter(),
InsertionSortCollider([
# enlarge here
Bo1_Sphere_Aabb(aabbEnlargeFactor=intRadius,label='bo1s'),
Bo1_Facet_Aabb(),
]),
InteractionLoop(
[
# enlarge here
Ig2_Sphere_Sphere_Dem3DofGeom(distFactor=intRadius,label='ig2ss'),
Ig2_Facet_Sphere_Dem3DofGeom(),
],
[Ip2_CpmMat_CpmMat_CpmPhys()],
[Law2_Dem3DofGeom_CpmPhys_Cpm(epsSoft=0)], # deactivated
),
NewtonIntegrator(damping=damping,label='damper'),
]
# run one single step
O.step()
# reset interaction radius to the default value
# see documentation of those attributes for the meaning of negative values
bo1s.aabbEnlargeFactor=-1
ig2ss.distFactor=-1
# now continue simulation
O.run()
*********Relavent Yade documentation ends here************

To remodel the triaxial compression process in [1], I tried to run a simpler triaxial test using CohFrictMat at first. There are 500 particles, and parameters goes as follows:

young = 5e6, normalcohesion = shearcohesion = 4500, confining pressure = 10kPa

Particles are packed randomly in a 0.15*0.15*0.3 box. Then confining pressure is reached by internal compaction. After a stable status is reached, count the number of all interactions. Then the aabbenlargefunctor and interactiondetection factor is altered to 1.25(or 1.5). Next step, count the total number of interactions again and the two factors are set to 1. 
This method didn't work correctly:

1. Whatever the two factor is, 1.25 or 1.5, total interactions before and after alteration show little difference. 

2. Besides, after the alteration, the utils.avgNumInteractions() is printed instantly and always turn out to be less than 7. According to [1], this number should be around 10 when intRaius is 1.25, and around 14 when intRadius is 1.5.

3. Stress and strain data is recorded, stress-strain curve show no obvious difference under different intRadius.

The code is like this:

*********************code starts here******************************
from yade import pack

num_spheres=500
targetPorosity = 0.44
compFricDegree = 0 
finalFricDegree = 20 
rate=-0.02 # loading rate (strain rate)
stabilityThreshold=0.01
young=5e6 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(0.15,0.3,0.15)
intRadius=1.25 # interaction range coefficient
confiningPressure=1e4  
normalAdhesion = 4500
shearAdhesion = 4500

x=0.00
num_broken = 0
num_intr = 0
num_coh = 0
num_coh_ini = 0
num_intr_unknown = 0
recordstep = 2000
coordn = 0.00
k=0

O.materials.append(CohFrictMat(young=young, poisson=0.333,frictionAngle=radians(compFricDegree),normalCohesion = normalAdhesion, shearCohesion = shearAdhesion, density=2600, label='spheres'))
O.materials.append(FrictMat(young=young, poisson=0.333,frictionAngle=0,density=0,label='walls'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
	maxMultiplier=1.+4e-3,
	finalMaxMultiplier=1.+1e-4,
	stressMask = 7,
	internalCompaction=True,
)

newton=NewtonIntegrator(damping=0.2)
triax.goal1=triax.goal2=triax.goal3=-confiningPressure

O.engines=[
    ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(label = 'bo1s'),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom6D(label = 'ig2ss'),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label = 'ip2cc')],
		[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
    triax,
    newton,
    PyRunner(command = 'isotropic_compaction()', iterPeriod = 1000, label = 'checker'),
    VTKRecorder (
        fileName = '/home/weijy/triax-test/triax/Paraview/1-', 
        recorders = ['spheres','boxes','intr'], 
        iterPeriod = recordstep, label = 'export')
]

def isotropic_compaction():
    global intRadius
    global num_intr
    unb=unbalancedForce()
    print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
    if unb<stabilityThreshold and abs(-confiningPressure-triax.meanStress)/confiningPressure<0.001:
        num_intr = 0
        bo1s.aabbEnlargeFactor = intRadius     ###THIS IS WHERE I ALTER THE FACTOR, MEANWHILE COHESION IS SET#####
        ig2ss.interactionDetectionFactor = intRadius
        ip2cc.setCohesionNow = True
        for i in O.interactions:
            num_intr = num_intr + 1
        print '*************************************************'
        print 'total interactions before reset intR is ', num_intr
        num_intr = 0
        checker.iterPeriod = 1                    ####ITERPERIOD IS SET TO 1 SO FACTORS CAN BE ALTERED AGAIN AT NEXT STEP######
        print 'current step is ',O.iter
        print "###      Isotropic Compaction Accomplished      ###"
        checker.command = 'set_interaction()'
        
def set_interaction():
    global coordn
    global num_intr
    global num_coh_ini
    global num_intr_unknown
    print '************************************************'
    print 'current step is ',O.iter
#    print "###      Setting Interactions      ###"
    bo1s.aabbEnlargeFactor = -1                                         ###########NOW THE FACTORS ARE SET TO 1###########
    ig2ss.interactionDetectionFactor = -1
    coordn = utils.avgNumInteractions()
    print coordn
    num_intr=0
    num_coh_ini=0
    num_intr_unknown=0
    for i in O.interactions:
        num_intr = num_intr + 1
        try:
            if i.phys.shearAdhesion > 0:
                num_coh_ini = num_coh_ini + 1
#                print i,' .phys.shearAdhesion = ', i.phys.shearAdhesion
#                print i,' .phys.normalAdhesion = ', i.phys.normalAdhesion
#            if i.phys.shearAdhesion == 0:
#                num_coh
        except:
            num_intr_unknown = num_intr_unknown + 1
    
    print '\n num_coh_ini = ', num_coh_ini
    print '\n num_intr_unknown =  ', num_intr_unknown
    print '\n total interactions after reset intR is ', num_intr
    
    num_intr = 0
    
    print "###    Interactions Established      ###"
    checker.command = 'start_compression()'


def start_compression():
    print '****************************************************8'
    print 'current step is ',O.iter
    global rate
    global finalFricDegree
    global confiningPressure
    triax.internalCompaction=False
    setContactFriction(radians(finalFricDegree))
    triax.stressMask = 5
    triax.goal2= rate
    triax.goal1= -confiningPressure
    triax.goal3= -confiningPressure
    newton.damping = 0.4

    O.engines = O.engines + [TriaxialStateRecorder(file = '/home/weijy/triax-test/triax/TriaxialState-test', iterPeriod = recordstep)]
    O.engines = O.engines + [PyRunner(command='check_contact()',iterPeriod = recordstep)]

    a = open('/home/weijy/triax-test/triax/result', "a")
    a.write('*********************************************************************\n')
    a.write('number_spheres = '+str(num_spheres)+'   contact friction = '+str(finalFricDegree)+'   load rate = '+str(rate)+'   intradius = '+str(intRadius)+'\n')
    a.write('young = '+str(young)+'   confiningpressure = '+str(confiningPressure)+'   c = '+str(normalAdhesion)+'   t = '+str(shearAdhesion)+'   N = '+str(coordn)+'\n')
    a.write('*********************************************************************\n\n')
    a.write('i e22 p q ev porosity num_broken num_intr\n')
    print "###    Start Compression      ###"
    a.close()
    checker.iterPeriod = 2000
    checker.command = 'check_stop()'
    O.pause()
    
def check_contact():
    global a
    global num_coh
    global num_broken
    global num_coh_ini
    global x
    global num_intr
    for i in O.interactions:       
        try:
            x = i.phys.normalAdhesion
            num_intr = num_intr + 1
#            print '\n normaladhesion = ', i.phys.normalAdhesion
#            print ' normalforce = ', i.phys.normalForce
#            print ' shearadhesion = ', i.phys.shearAdhesion
#            print ' shearforce = ', i.phys.shearForce
#            print ' kn= ',i.phys.kn
#            print ' ks= ',i.phys.ks
#            print ' tangensOfFrictionAngle = ' , i.phys.tangensOfFrictionAngle
            if x > 0:
                num_coh = num_coh + 1
        except:
            x = 1
            
    num_broken = num_coh_ini - num_coh
    print '\n\n iter =  ', O.iter    
    print ' num_coh_ini = ', num_coh_ini
    print ' num_broken = ' ,num_broken    
    print ' totalinteracions between particles= ' , num_intr

    e11=-triax.strain[0]
    e22=-triax.strain[1]
    e33=-triax.strain[2]
    ev=triax.strain[0]+triax.strain[1]+triax.strain[2]
    s11=-triax.stress(triax.wall_right_id)[0]
    s22=-triax.stress(triax.wall_top_id)[1]
    s33=-triax.stress(triax.wall_front_id)[2]
    p=(s11+s22+s33)/3
    q=s22-(s11+s33)/2
    poro=triax.porosity
    num_broken = float(num_broken) / num_coh_ini
    a = open('/home/weijy/triax-test/triax/result', "a")
    a.write(str(O.iter)+' '+str(e22)+' '+str(p)+' '+str(q)+' '+str(ev)+' '+str(poro)+' '+str(num_broken)+' '+str(num_intr)+'\n')
    
    num_intr = 0
    num_coh = 0
            
    
def check_stop():
    if -triax.strain[1] > 0.3:
        print 'End'
        O.pause()

********************code ends****************************





















-- 
You received this question notification because your team yade-users is
an answer contact for Yade.