confining pressure is not stable

Asked by ytang on 2020-11-04

Hi All,

I'm trying to use the Hertz constitutive model to simulate the triaxial tests.
#######################################################
I use the code from the example. [1]
Here is the code:
####################
from yade import pack

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
 num_spheres=20000,# number of spheres
 compFricDegree = 30, # 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.368 #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 = 30 # contact friction during the deviatoric loading
rate=-0.02 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=5e6 # contact stiffness
mn,mx=Vector3(0,0,0),Vector3(1,1,1) # corners of the initial packing

## create materials for spheres and plates
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'))

## 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.3333,num_spheres,False, 0.95,seed=1) #"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.+9e4/young, # spheres growing factor (fast growth)
 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
 thickness = 0,
 stressMask = 7,
 internalCompaction=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

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()]
  [Ip2_FrictMat_FrictMat_MindlinPhys( betan =0.2 , betas = 0.2, krot = 0.2 , eta = 0.5 )],
  [Law2_ScGeom_MindlinPhys_Mindlin()]
 ),
 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()

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

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(-100000-triax.meanStress)/100000<0.001:
    break

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

###################################################
### 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()
 ## 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('compactedState'+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=-100000
triax.goal3=-100000

##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()

#####################################################
### Example of how to record and plot data ###
#####################################################

from yade import plot

### a function saving variables
def history():
 plot.addData(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],
   i=O.iter)
 plot.saveDataTxt('calibration-hz.txt')

if 1:
  ## include a periodic engine calling that function in the simulation loop
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
  ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
else:
  O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')

O.run(100000,True)

#plot.plots={'e22':('s11','s22','s33',None,'ev')}

## display on the screen (doesn't work on VMware image it seems)
#plot.plot()
#######################################################################
#######################################################################
I used these two models: [Law2_ScGeom_FrictPhys_CundallStrack()] and [Law2_ScGeom_MindlinPhys_Mindlin()]
these two models run without errors.
###############################
For both models, I set the confining pressure to 100kPa.
As for the Cundall model, the confining pressure can keep constant. However, for the hertz model, the confining pressures in the x and z direction are not the same. you can find the results here [2].
#################################
All the code is the same except for the constitutive model. I want to know why the confining pressure is not constant.

best,
Yong

References: [1] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/triax-tutorial/script-session1.py
[2] https://www.dropbox.com/sh/twdl0kfx6bfo8mt/AABamjTzasooR7rum0KXc0b4a?dl=0

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-11-04
Last reply:
2020-11-04

This question was reopened

Jan Stránský (honzik) said : #1

Hello,

always try to make your code a MWE [1], M=minimum code w.r.t. the problem. It is:
- easier for others to find problems
- easier for YOU to find the problem

E.g. here, where the problem is confining (isotropic?) pressure, why is there the deviatoric loading?
One reason why pressure is not constant is that the script switched to the deviatoric loading phase. Please check it.

cheers
Jan

[1] https://www.yade-dem.org/wiki/Howtoask

ytang (ytang116) said : #2

Hi Jan,

I know that during the isotropic compression stage the stresses in the three directions are the same. If the stresses are not the same, there will be no deviatoric stage.

As we all know, in a triaxial test, the confining pressure should be a constant value even if in the deviatoric loading stage. (if the deviatoric loading is in the Y direction, the confining pressure in the X and Z directions should be the same value.)

I use two models. all the code is the same except for the constitutive model.
The confining pressure in the Cundall model is the same in the X and Z directions, which is 100kPa, while the confining pressure in the Hertz model is not the same (one is 110kPa, another is 90kPa).

I'm wondering why the confining pressure is not the same during the deviatoric stage?

Jan Stránský (honzik) said : #3

Sorry, I have misinterpreted the question.

A few tips:
- check triax strain rate and compare to triax's max strain rate
- try with different "rate" values
- what are strains compared to the two models?

cheers
Jan

Can you help with this problem?

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

To post a message you must log in.