Wierd volumetric strain during triaxial tests under different confining pressure

Hello everyone,

I did my triaxial compression tests under different confining pressures. Sigma3 in the horizontal directions is constant all through the test and the cubic sample is compressed in the vertical direction with a constant strain rate. Sigma3 is equal to 100kPa or 400kPa in my tests. Both samples contract first and dilate later as the compression goes on.

As well all know, larger volume contraction and lower volume dilation afterwards will occur through the test under higher sigma3. But my results are just reverse to this. In the case when sigma3 is 100kPa, the volume contraction is larger and the volume dilation is maller afterwards. I really don't know why this happens.

My sample is trimmed from a big gravity deposited sample. My friend obtained his sample by expanding spheres inside a cubic box and found the same weird deformation. We use the same contact law. We have adjusted the elastic modulus and the poisson ratio, but it doesn't work.

I did triaxial tests with PFC with linear contact model under different confining pressures and everything is fine.

Does anyone else meet this kind of problem?

Below is my code:

from yade import ymport
from yade import utils

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

compFricDegree = 26.565 # initial contact friction during the confining phase
finalFricDegree = 26.565 # contact friction during the deviatoric loading
rate=-0.01 # loading rate (strain rate)
damp=0.7 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=5e8 # 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(FrictViscoMat(betan=0.2, young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictViscoMat(betan=0.2, young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))

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

print "Load clumps"
deposited_bodies = ymport.textClumps('cubic_sample_clump_miu0.5_dep_damp0.1_35608_hor_layer_plane.txt', shift=Vector3(0,0,0), material='spheres')
O.bodies.updateClumpProperties()#get more accurate clump masses/volumes/inertia

############################
### 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=False, # If true the confining pressure is generated by growing particles
)

contact = Law2_ScGeom_FrictViscoPhys_CundallStrackVisco(traceEnergy=True)
newton=NewtonIntegrator(damping=damp)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
  [Ip2_FrictViscoMat_FrictViscoMat_FrictViscoPhys()],
  [contact]
 ),
 ## 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,
 newton,
 # save data for Paraview
 VTKRecorder(fileName='post/3d-vtk-',recorders=['all'],iterPeriod = 500000)
]

O.dt=.5*utils.PWaveTimeStep()

#######################################
### APPLYING CONFINING PRESSURE ###
#######################################

## we define the lateral stresses during the test, here the same 100kPa as for the initial confinement.
p_const = 100000

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

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

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

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

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

##set stress control on x and y, we will impose strain rate on z
triax.stressMask = 3
##now goal2 is the target strain rate
triax.goal3=rate

O.run(4000000,True)

O.save('final.yade.gz')