SetCohesionNow & setCohesionOnNewContacts /triaxial test

[Seti]

Hello all,

can I turn on setCohesionNow, setCohesionOnNewContacts in triaxial test as well? if yes, why I set these parameters the porosity & friction does not change and I face with below error? while without these consideration code works.

 "Friction: 33.25 porosity: 1.0python: malloc.c:3720: _int_malloc: Assertion `(unsigned long) (size) >= (unsigned long) (nb)' failed.
Aborted (core dumped)"

the triaxial script is not different from the original one, I just copy it here if you need to know about the inputs.

Thanks,
Seti

from yade import pack,plot

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
 num_spheres=1000,# number of spheres
 compFricDegree =35, # contact friction during the confining phase
 key='_triax_base_', # put you simulation's name here
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.42 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 35# contact friction during the deviatoric loading
rate=-0.005 # loading rate (strain rate)
damp=0.3 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=100e6# contact stiffness
mn,mx=Vector3(0,0,0),Vector3(0.09,0.18,0.09) # corners of the initial packing

## create materials for spheres and plates
O.materials.append(CohFrictMat(alphaKr=0.5,young=young,poisson=0.09,frictionAngle=radians(33.5),normalCohesion=7.5e3,shearCohesion=2.25e3,momentRotationLaw=True,etaRoll=0.001,density=2600,isCohesive=True,label='spheres'))
O.materials.append(CohFrictMat(young=young,poisson=0,frictionAngle=radians(0),density=0,label='walls'))

## create walls around the packing
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()

clumps=False #turn this true for the same example with clumps
if clumps:
 ## approximate mean rad of the futur dense packing for latter use
 volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
 mean_rad = pow(0.09*volume/num_spheres,0.3333)
 ## define a unique clump type (we could have many, see clumpCloud documentation)
 c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
 ## generate positions and input them in the simulation
 sp.makeClumpCloud(mn,mx,[c1],periodic=False)
 sp.toSimulation(material='spheres')
 O.bodies.updateClumpProperties()#get more accurate clump masses/volumes/inertia
else:
 sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
 #sp.makeCloud(mn,mx,0.066,num_spheres) #"seed" make the "random" generation always the same
 O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
 #or alternatively (higher level function doing exactly the same):
 #sp.toSimulation(material='spheres')

############################
### DEFINING ENGINES ###
############################

triax=TriaxialStressController(
 ## TriaxialStressController will be used to control stress and strain. It controls particles size and plates positions.
 ## this control of boundary conditions was used for instance in http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 ## switch stress/strain control using a bitmask. What is a bitmask, huh?!
 ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, which are 1 or 0.
 ## Then an integer uniquely defining the combination of all these tests is: mask = x*1 + y*2 + z*4
 ## to put it differently, the mask is the integer whose binary representation is xyz, i.e.
 ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x and y and z", etc.
 stressMask = 7,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)
########################################
#Modified engine
##################################
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys (),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow = True, setCohesionOnNewContacts = True,label="cohesiveIp")],
                [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom_CohFrictPhys_CohesionMoment(
   useIncrementalForm=True, #useIncrementalForm is turned on as we want plasticity on the contact moments
   always_use_moment_law=False, #if we want "rolling" friction even if the contact is not cohesive (or cohesion is broken), we will have to turn this true somewhere
   label='cohesiveLaw')]
        ),
        ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax,
        TriaxialStateRecorder(iterPeriod=100,file='150,damp0.8,rate 0.005,NEW50,alphaKr=0.5,young=100e6,poisson=0.09,frictionAngle=radians(50),normalCohesion=7.5e10,shearCohesion=2.25e10,etaRoll=0.025,density=2600,wall35,'+key),
        newton
]
##########################################################
#O.engines=[
 #ForceResetter(),
 #InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 #InteractionLoop(
  #[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  #[Ip2_FrictMat_FrictMat_FrictPhys()],
  #[Law2_ScGeom_FrictPhys_CundallStrack()]
 #),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 #GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 #triax,
 #TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 #newton
#]
#############################
#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

## UNCOMMENT THE FOLLOWING SECTIONS ONE BY ONE
## DEPENDING ON YOUR EDITOR, IT COULD BE DONE
## BY SELECTING THE CODE BLOCKS BETWEEN THE SUBTITLES
## AND PRESSING CTRL+SHIFT+D
#if nRead==0: yade.qt.Controller(), yade.qt.View()
print 'Number of elements: ', len(O.bodies)
print 'Box Volume: ', triax.boxVolume
#######################################
### APPLYING CONFINING PRESSURE ###
#######################################

#the value of (isotropic) confining stress defines the target stress to be applied in all three directions
triax.goal1=triax.goal2=triax.goal3=-150000

#while 1:
  #O.run(1000, True)
  ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  #unb=unbalancedForce()
  #print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  #if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
    #break

#O.save('confinedState'+key+'.yade.gz')
#print "### Isotropic state saved ###"

###################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
###################################################

### We will reach a prescribed value of porosity with the REFD algorithm
### (see http://dx.doi.org/10.2516/ogst/2012032 and
### http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)

import sys #this is only for the flush() below
while triax.porosity>targetPorosity:
 ## we decrease friction value and apply it to all the bodies and contacts
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 ## while we run steps, triax will tend to grow particles as the packing
 ## keeps shrinking as a consequence of decreasing friction. Consequently
 ## porosity will decrease
 O.run(500,1)

O.save('compactedStateBEL20,young=63.9e8'+key+'.yade.gz')
print "### Compacted state saved ###"

##############################
### DEVIATORIC LOADING ###
##############################

##We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False

## Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

##set stress control on x and z, we will impose strain rate on y
triax.stressMask = 5
##now goal2 is the target strain rate
triax.goal2=rate
## we define the lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-150000
triax.goal3=-150000

##we can change damping here. What is the effect in your opinion?
newton.damping=0.1

[Bruno Chareyre (bruno-chareyre)]

Hi, I realize that there are plenty differences between your script and the example script, depsite your statement that it is similar to the original.
Please be more accurate on which change caused the crash, and provide us with a way to pinpoint it (e.g. triggering the crash or not by changing just one line).
Bruno

[Seti (seti)]

Hi Bruno, I just tried to change material to cohesion one consequently added cohesive IP and law, would you please advise me if it is not correct?

[Bruno Chareyre (bruno-chareyre)]

Hi, what you describe in #2 sounds ok, in principle.
Since there are many other changes compared to the original script I can't tell if this is really the origin of your problem.
You may try to start again from the original script and re-apply each change to see when the problem appears.
B

[Seti (seti)]

Hi Bruno,
The changes are as per below: ( As mentioned before just material property and engine)

Can you please advise?

original code:

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
....
.
.
.

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
 newton
]

###############

Current code:

O.materials.append(CohFrictMat(alphaKr=0.5,young=young,poisson=0.25,frictionAngle=radians(36.28),normalCohesion=7.5e3,shearCohesion=2.25e3,momentRotationLaw=True,etaRoll=0.001,density=2600,isCohesive=True,label='spheres'))
O.materials.append(CohFrictMat(young=young,poisson=0.25,frictionAngle=radians(36.28),density=0,label='walls'))
......
.
.
.

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],
  #[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp")],
                [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom_CohFrictPhys_CohesionMoment(
   useIncrementalForm=True, #useIncrementalForm is turned on as we want plasticity on the contact moments
   always_use_moment_law=True, #if we want "rolling" friction even if the contact is not cohesive (or cohesion is broken), we will have to turn this true somewhere
   label='cohesiveLaw')]
  #[Law2_ScGeom_CohFrictPhys_CohesionMoment(
   #useIncrementalForm=True, #useIncrementalForm is turned on as we want plasticity on the contact moments
   #always_use_moment_law=True, #if we want "rolling" friction even if the contact is not cohesive (or cohesion is broken), we will have to turn this true somewhere
   #label='cohesiveLaw')]
        ),
        ## We will use the global stiffness of each body to determine an optimal timestep (see https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax,
        TriaxialStateRecorder(iterPeriod=100,file='150,damp0.3,rate 0.005,NEW50,alphaKr=0.5,young=150e6,poisson=0.25,frictionAngle=radians(36.28),,etaRoll=0.025,density=2600,wall36.28,'+key),
        newton
]

[Seti (seti)]

When I am turning on:
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,label="cohesiveIp")],

instead of

Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],

I will face with error!

Regards,
Seti

[Bruno Chareyre (bruno-chareyre)]

Interesting. There is still more than one difference between the two lines. Could you find which change shows the problem?
B

[Janek Kozicki (cosurgi)]

moved from https://bugs.launchpad.net/yade/+bug/1663899 to https://gitlab.com/yade-dev/trunk/issues/33