Triaxial test comparison between periodic boundary and rigid boundary

Hello,

> I carried out ...

a MWE [1] would be great, most likely necessary, to provide serious answer.
Or at least more detailed information on "triaxial compression tests" (loading scenario, load values, what is prescribed (stress, strain), .... )

> the two cases use ... same strain rate (I think)

you should know, not think

> the two cases use ... the same of sample size

definitely not.
For the rigid walls, the sample dimension is a true dimension.
For the size of periodic cell, it just an arbitrary sample out of **infinite** medium.
If you "repeat" the periodic cell (like "2x2x2" or "3x3x3" sample), you should get the same results as for the one cell.

> but is it normal to have such a big difference?

it might be normal, or not :-)

> Or did I miss some factors which can also lead to the difference?

friction of walls?

Cheers
Jan

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

Hello Leonard,

When you say "same initial void ratio", does "initial" mean just before
the triaxial compression?
Similarly did you check the coordination number just before the triaxial
compression is the same for the two cases?

Best,
Luc

Le 01/09/2021 à 18:35, Leonard a écrit :
> New question #698594 on Yade:
> https://answers.launchpad.net/yade/+question/698594
>
> Hi,
> I carried out triaxial compression tests at two boundary conditions (one case with rigid wall and one case with periodic boundary), the two cases use the same material, the same particle size distribution and paticle number, same contact law and the same initial void ratio, at same confining pressure, same strain rate (I think), and almost the same of sample size. But I got really different mechanical response, especially for the volumetric behaviour, the periodic one shows contractive behaviour, while the rigid case shows dilative behaviour. I thoungh I have controlled all the conditions (except boundary condition) the same.
> It is not surprise that there are some difference, but is it normal to have such a big difference? Or did I miss some factors which can also lead to the difference?
>
> Thanks!
> Leonard
>

--
Luc Sibille
Université Grenoble Alpes / IUT1 de Grenoble
Laboratoire 3SR: Sols, Solides, Structures, Risques

Tel lab.: +33 (0)4 76 82 63 48
Tel IUT: +33 (0)4 76 82 53 36

Hello Jan and Luc,
Thanks very much for your reply.

Here I made two MWE for rigid wall case and periodic boundary case respectively. It may take around 8 min for each of the simulation, or you can also find the results of stress-strain and volumetric strain evolution at [1].

>same strain rate (I think)
Because in rigid wall which uses TriaxialStressController, I can specify the loading rate by triax.goal. But as far as I understand, it seems that for PeriTriaxController, it can only specify the loading rate by maxStrainRate. So I made triax.strain=1 (rigid case) and triax.maxStrainRate = (1., 1., 1.) in periodic case. Do you know how to specify the strain rate in PeriTriaxController?

>the two cases use ... the same of sample size
Now I know they are not the same. Just explain here why I though they were same before. If you run the periodic MWE below, it prints out the cell size just before the deviatoric loading, this size is almost the same as the size in rigid case.

>initial porostiy
Sorry for not clearly describe it. Yes, the initial porosity I mentioned is the porosity at the begining of the deviatoric loading. In the following two MWE, they are both 0.358 getting by function of porosity(). But note that, in rigid boundary case, I also print the porosity at the same state by using triax.porosity, which retures 0.33.

>coordination number
Thanks Luc, that is a good point to mention. In the two MWE, I print the coordination number (CN) just before the deviatoric loading. And the CN of rigid boundary is 7.28, while the CN in periodic is 6.3. Which means there are more contacts in rigid wall case, indicating a denser state in general. Then my question is why they have same inputtings (porosity, particle number...) but they have different CN? And another question is how can I make the same sample under the two boundary conditions? It would be conflict if I control them as the same from porosity and from CN.

################################### rigidWall case ##########################
from __future__ import division
from yade import pack, plot
import numpy as np
num_spheres=7000# number of spheres
targetPorosity = 0.33
compFricDegree = 30
finalFricDegree = 35
damp=0.6
stabilityThreshold=0.001
confinement=100e3
mn,mx=Vector3(0,0,0),Vector3(0.07,0.14,0.07) #determine the size of the sample
young=5e6
## create materials for spheres and plates
MatWall=O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0,density=0,label='walls'))
MatSand = O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(30),\
           density=2650.0,label='sand'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='sand') for center,rad in sp])

Gl1_Sphere.quality=3
triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 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()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax,
 TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
 newton,
 PyRunner(iterPeriod=100,command='history()',label='recorder'),
 PyRunner(command='stopIfDamaged()',iterPeriod=200,label="endSimulation"),
]

recorder.dead=True
endSimulation.dead=True
Gl1_Sphere.stripes=0
triax.goal1=triax.goal2=triax.goal3=-confinement

def getCN():
 # get coordination number
 CN=0
 for i in O.bodies:
  if isinstance(i.shape,Sphere):
   CN +=len(i.intrs())
 Z = CN/ 7000.0
 return Z

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbF:',unb,' meanStress: ',-triax.meanStress,'top:',-triax.stress(triax.wall_top_id)[1]
  if unb<stabilityThreshold and abs(-confinement-triax.meanStress)/confinement<0.0001:
   break

print "### state 1 Isotropic completed ###"

import sys
while triax.porosity - targetPorosity>0.0001:
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)

print "### state 2 Reach target porosity completed ###"
print "particle radii is: ", O.bodies[2000].shape.radius ## here get the particle radius which will be used as input in makecloud in periodic boundary case.
print 'triax.porosity at 100 kPa',triax.porosity,'porosity()',porosity()
print ("coordination number at 100 kPa is",getCN())
topStress_initial=-triax.stress(triax.wall_top_id)[1]
bottomStress_initial=-triax.stress(triax.wall_bottom_id)[1]

# start deviatoric loading
triax.internalCompaction=False
triax.wall_bottom_activated=False
recorder.dead=False
endSimulation.dead=False
triax.stressMask=5
triax.goal2=-1
triax.goal1=-confinement
triax.goal3=-confinement
setContactFriction(radians(finalFricDegree))

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],
   q=-triax.stress(triax.wall_top_id)[1]-100000,)

def stopIfDamaged():
 if -triax.strain[1]>0.4:
  O.pause()
  plot.saveDataTxt('results_rigidWall.txt')

plot.plots={'e22':('q',),'e22 ':('ev')}
plot.plot()
###################################### end of the 1st MWE #############

#################### periodic case #######################3
from __future__ import print_function
sigmaIso = -100e3
from yade import pack, qt, plot

O.periodic = True
O.cell.setBox(.09,.18,.09)
young=5e6
compFricDegree=23 # using this value can make the porosity at 100 kpa be 0.358, which is the same as that of rigid wall
finalFricDegree = 35
MatSand = O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(30),density=2650.0,label='sand'))
sp = pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.1,0.2,0.1)
sp.makeCloud(mn,mx,0.0025,0,7000,False, 0.95,seed=1)
sp.toSimulation(material="sand")
# here set contact friction angle to reach a low porosity at target confining pressure.
setContactFriction(radians(compFricDegree))
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
        PeriTriaxController(
                label='triax',
                goal=(sigmaIso, sigmaIso, sigmaIso),
                stressMask=7,
                dynCell=True,
                maxStrainRate=(10, 10, 10),
                maxUnbalanced=.1,
                relStressTol=1e-3,
                doneHook='compactionFinished()'
        ),
        NewtonIntegrator(damping=.6),
        PyRunner(command='addPlotData()', iterPeriod=100,label="recorder"),
]
O.dt = .5 * PWaveTimeStep()
recorder.dead=True

def addPlotData():
 plot.addData(
         unbalanced=unbalancedForce(),
         i=O.iter,
         sxx=-triax.stress[0],
         syy=-triax.stress[1],q=-triax.stress[1]-100000,
         szz=-triax.stress[2],
         exx=-triax.strain[0],ev=triax.strain[0]+triax.strain[1]+triax.strain[2],
         eyy=-triax.strain[1],
         ezz=-triax.strain[2],

 )

O.trackEnergy = True

plot.plots = {
        'eyy': ('q',),
        'eyy ': ('ev'),}
plot.plot()

def compactionFinished():
 print ("porosity at 100 kPa is",porosity())
 print ("coordination number at 100 kPa is",getCN())
 print ("O.cell.size is", O.cell.size)
 setContactFriction(radians(finalFricDegree))
 O.pause()
 O.cell.trsf = Matrix3.Identity
 triax.goal = (sigmaIso,-0.4, sigmaIso)
 triax.stressMask = 5
 triax.maxStrainRate = (1., 1., 1.)
 triax.doneHook = 'triaxFinished()'
 triax.maxUnbalanced = 10
 O.step() # if I didnot use O.step(), the initial eyy will not start from 0.
 recorder.dead=False
def triaxFinished():
 print('Finished')
 O.pause()
 plot.saveDataTxt('results_periodic.txt')
def getCN():
 CN=0
 for i in O.bodies:
  if isinstance(i.shape,Sphere):
   CN +=len(i.intrs())
 Z = CN/ 7000.0
 return Z
##################################################

Thanks
Leonard
[1]https://we.tl/t-xlkvIPB9zx
ps: both of the MWE are modified from YADE trunk.

This question was expired because it remained in the 'Open' state without activity for the last 15 days.