Changing friction angle in triaxial test

Asked by Seti

Hi guys!
I am trying to do parametric study for triaxial test. I have played around the friction angle and alphaKr. However changing of these two parameters have no effect on my results. I am not sure where is my mistake . Would you please advice me.

Thanks heaps,

Here is my code:
from yade import pack
nRead=utils.readParamsFromTable(
        num_spheres=1000,# number of spheres len(O.bodies) to verify: 10006 = 10000 particles + 6 walls is correct
        compFricDegree = 30, # contact friction during the confining phase (1)
        unknownOk=True,
        isoForce=200000, # stress for the isotropic compression phase (1)
        conStress=200000 # confinement stress, for the deviatoric loading session (2)
)
from yade.params import table

num_spheres=table.num_spheres # number of spheres called from table
targetPorosity = 0.8 #the porosity we want for the packing (3 specimens: (Ei,n) = (1,0.382), (2,0.387), (3,0.409) )
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=0.001 # loading rate (strain rate)
damp=0.3 # damping coefficient
stabilityThreshold=0.001 # initial value: 0.001
key='_kozooooooooooooooocki,,200kpa,200kpa,30,30,30e9_0.12,0.12,0.12,alphaKr=0.15,BREAK0.7,.001' # simulation's name here
young=30e9 # contact stiffness k_n/Ds
mn,mx=Vector3(-0.12,-0.12,-0.12),Vector3(0.12,0.12,0.12) # corners of the initial packing
thick = 0.01 # thickness of the plates

## create materials for spheres and plates
O.materials.append(CohFrictMat(alphaKr=0.15,young=young,poisson=0.3,frictionAngle=radians(compFricDegree),momentRotationLaw=True,etaRoll=1,density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

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

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
psdSizes=[0.002,0.003,0.004,0.005,0.006,0.007,0.008] # (sizes or radii of the grains vary from 2mm to 9.5mm)
#psdCumm=[1,9,25,50,69,90,95,100] # the correspondent amount (percentage) of each diameter, uncomment for yade-stable
psdCumm=[0.14,0.28,0.34,0.50,0.65,0.85,1.00] # for the code not using percentage, e.g. yade-daily
sp.makeCloud(mn,mx,0.00575,0.6,num_spheres,False,0.95,seed=1)
sp.toSimulation(material='spheres')
triax=TriaxialStressController(
        maxMultiplier=1.001, # spheres growing factor (fast growth), validated only when internalCompaction = True
        finalMaxMultiplier=1.00001, # spheres growing factor (slow growth), validated only when internalCompaction = True
        thickness = thick,
stressMask = 7,
        #the value of confining stress for the intitial (growth) phase
        goal1=table.isoForce,
        goal2=table.isoForce,
        goal3=table.isoForce,
        max_vel=0.5, # validated only when internalCompaction = False m/s
        internalCompaction=True,)
newton=NewtonIntegrator(damping=damp)

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_FrictPhys (),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],
                [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_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')]
        ),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
        triax,
        TriaxialStateRecorder(iterPeriod=50,file='WallStresses'+key),
        newton
]
### APPLYING CONFINING PRESSURE ###
#######################################

while 1:
  O.run(1000, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
  unb=unbalancedForce()
  #average stress
  #note: triax.stress(k) returns a stress vector, so we need to keep only the normal component
  meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  print 'unbalanced force:',unb,' mean stress: ',meanS
  print 'void ratio=',triax.porosity/(1-triax.porosity), 'porosity=', triax.porosity
  print 'mean stress of engine', triax.meanStress
  print '0.12,0.12,0.12,alphaKr=0.15,BREAK0.7,.001'
  if triax.porosity<0.7:
    break
##################################################
### REACHING A SPECIFIED POROSITY PRECISELY ###
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()
        O.run(500,1)
##############################
### 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))
setContactFriction(radians(30))

#set stress control on x and z, we will impose strain rate on y (5)
triax.stressMask = 5
#now goal2 is the target strain rate
triax.goal2=-rate
#we assign constant stress to the other directions
triax.goal1=table.conStress
triax.goal3=table.conStress

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

##Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()

while triax.strain[1] < 0.25:
  O.run(50); O.wait()

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#1

Hi,
You pasted the tutorial script or almost.
Did you check if changing the parameters of the original script give different results (I know it does, but...)?

To investigate what happens with alphaKr you should check individual interactions and see if they really have a torque computed.

Bruno

Can you help with this problem?

Provide an answer of your own, or ask Seti for more information if necessary.

To post a message you must log in.