yade-users team mailing list archive

Thread
Date
Re: [Question #482413]: TriaxialStressController lose control of stress

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: Wei Cao <question482413@xxxxxxxxxxxxxxxxxxxxx>
Date: Sat, 25 Feb 2017 22:28:26 -0000
Reply-to: question482413@xxxxxxxxxxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx
Question #482413 on Yade changed:
https://answers.launchpad.net/yade/+question/482413

    Status: Answered => Open

Wei Cao is still having a problem:
Hi, thanks again. I modified my code to make sure that the sample is at equilibrium before applying the deviatic stress, by adding:
{
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
}


Could you helpt me to look at the code? I only want to know why the strain control give me so strange results.  It really bothers me for a long time. Thank you！

Here is my code:

##############################################################################################################################
#Authors: Luc Sibille  luc.sibille@xxxxxxxxxxxxxxx
#         Franck Lomine  franck.lomine@xxxxxxxxxxxxxx
#
#   Script to simulate buoyancy in a granular assembly with the DEM-LBM coupling.
#   This script was written for a practical work in the OZ ALERT school in Grenoble 2011.
#   Script updated in 2014
##############################################################################################################################


from yade import pack, timing, utils, qt, plot
############################################
###   DEFINING VARIABLES AND MATERIALS   ###
############################################

# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
	num_spheres=10000,# 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.5 #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.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
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))


## control to open the file
Not_Open_flag = 1
if Not_Open_flag:
   stress_info = open("/home/caowei/Documents/myYade/trunk/examples/stress_info.txt",'a')  
   contact_info = open("/home/caowei/Documents/myYade/trunk/examples/contactOri.txt",'a') 
   particle_info = open("/home/caowei/Documents/myYade/trunk/examples/particle.txt",'a')
   Not_Open_flag = 0

sigmaIso = -1e5

rMean = 0.035  #mean radius (m) of solid particles

Gravity=-0.1 #gravity (m.s-2) that will be applied to solid particles

#definition of the simulation domain (position of walls), the domain is a 3D domain but particles will be generated on a unique plane at z=0, be carefull that wall normal to z directions do not touch particles.
lowerCornerW = (0.001,0.001,0.001)
upperCornerW = (0.701,0.701,0.701)

#definitions of the domain for particles: positions in z direction is set to 0 to force a 2D granular assembly 
#lowerCornerS = (0.00001,0.00001,0)
#upperCornerS = (0.00701,0.00501,0)

lowerCornerS = (0.001,0.001,0.001)
upperCornerS = (0.701,0.701,0.701)

nbSpheres = 5000 #this not the real number of particles but number of
times where we try to insert a new particle, the real number of
particles can be much lower!

def addPlotData():
   plot.addData(unbalanced=unbalancedForce(),i=O.iter,
      sxx=triax.stress(triax.wall_right_id)[0],syy=triax.stress(triax.wall_top_id)[1],szz=triax.stress(triax.wall_front_id)[2],
      exx=triax.strain[0],eyy=triax.strain[1],ezz=triax.strain[2],
      # save all available energy data
      Etot=O.energy.total(),**O.energy
   )

   stress_info.write("%6.2f " %(triax.stress(triax.wall_right_id)[0]))
   stress_info.write("%6.2f " %(triax.stress(triax.wall_top_id)[1]))
   stress_info.write("%6.2f " %(triax.stress(triax.wall_front_id)[2]))
   #print("stress:")
   #print(triax.stress(triax.wall_right_id)[0])
   #print(triax.stress(triax.wall_top_id)[1])
   #print(triax.stress(triax.wall_front_id)[2])
   #print("strain:")
   #print(triax.strain[0])
   #print(triax.strain[1])
   #print(triax.strain[2])
   for i in range(3):
      stress_info.write("%6.2f " %(triax.strain[i]))
   stress_info.write("\n")


   if O.iter % 20000 == 0:

      for i in O.interactions: 
         contact_info.write("%d " %(i.id1))	
         contact_info.write("%d " %(i.id2))	
         contact_info.write("%.4f " %(i.phys.normalForce[0]))
         contact_info.write("%.4f " %(i.phys.normalForce[1]))
         contact_info.write("%.4f " %(i.phys.normalForce[2]))
         contact_info.write("%.4f " %(i.phys.shearForce[0]))
         contact_info.write("%.4f " %(i.phys.shearForce[1]))
         contact_info.write("%.4f " %(i.phys.shearForce[2]))
         contact_info.write("\n")

      for i in range(len(O.bodies)): 
         if i > 5:
            particle_info.write("%d " %(i))
            if O.bodies[i].isClump:
               particle_info.write("\n")
            else:
               particle_info.write("%f " %(O.bodies[i].state.pos[0]))	
               particle_info.write("%f " %(O.bodies[i].state.pos[1]))	
               particle_info.write("%f " %(O.bodies[i].state.pos[2]))	
               particle_info.write("%f " %(O.bodies[i].shape.radius))	
               particle_info.write("%d " %(O.bodies[i].clumpId))	
               particle_info.write("\n")


# enable energy tracking in the code
O.trackEnergy=True

# define what to plot
plot.plots={'i':('unbalanced',),'i ':('sxx','syy','szz'),' i':('exx','eyy','ezz'),
   # energy plot
   ' i ':(O.energy.keys,None,'Etot'),
}
# show the plot
plot.plot()

def compactionFinished():
   # set the current cell configuration to be the reference one
   # O.cell.trsf=Matrix3.Identity
   triax.depth0 = triax.depth
   triax.height0 = triax.height
   triax.width0 = triax.width
   triax.strain[0] = triax.strain[1] = triax.strain[2] = 0


   # change control type: keep constant confinement in x,y, 20% compression in z
   triax.stressMask = 0
   triax.goal1 =  0.001
   triax.goal2 =  0.001
   triax.goal3 = -0.002

   count = 0
   countN = 2
   while count <= countN:

      for i in range(100):
         O.step()
      
      if triax.stress(triax.wall_front_id)[2] < -1.1e5:
         if triax.goal3 < 0:
            count = count + 1
         triax.goal1 = -0.001
         triax.goal2 = -0.001
         triax.goal3 =  0.002

      if triax.stress(triax.wall_front_id)[2] > -0.9e5:
         if triax.goal3 > 0:
            count = count + 1
         triax.goal1 =  0.001
         triax.goal2 =  0.001
         triax.goal3 = -0.002
     
   stress_info.close()
   contact_info.close()
   particle_info.close()
   print 'Finished'
   O.pause()
   # allow faster deformation along x,y to better maintain stresses
   # triax.maxStrainRate=(0.05,0.05,.1)
   # next time, call triaxFinished instead of compactionFinished
   # triax.doneHook = 'triaxFinished()'
   # do not wait for stabilization before calling triaxFinished
   # triax.maxUnbalanced =1 0
   # O.run(10000, True)

#def triaxFinished():
   #stress_info.close()
   #contact_info.close() 
   #particle_info.close()
   #print 'Finished'
   #O.pause()

newton=NewtonIntegrator(damping=damp)

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
   # call this function when goal is reached and the packing is stable
)
## Build of the engine vector
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()]
	),

	triax,	
	newton,
	## 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),
	## NewtonIntegrator(damping=0.2,gravity=(Gravity,0.,0.),label="ENewton"),
	PyRunner(command='addPlotData()',iterPeriod=100),
]


###############################################################################################################
# Creation of 6 walls at the limit of the simulation domain

# defintiion of the material for walls
O.materials.append(FrictMat(young=50e6,poisson=.5,frictionAngle=0.0,density=3000,label='walls'))

#lowerCornerW = (0.001,0.001,0.001)
#upperCornerW = (0.701,0.701,0.701)

mn,mx=Vector3(-0.1,-0.1,-0.1),Vector3(0.8,0.8,0.8) # corners of the
initial packing

## create materials for spheres and plates

O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=0.0,density=0,label='wallsss'))

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



###############################################################################################################
# Creation of the assembly of particles

# defintiion of the material for particles
O.materials.append(FrictMat(young=50e6,poisson=.5,frictionAngle=35*3.1415926535/180,density=3000,label='spheres'))
## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
#sp.makeCloud(lowerCornerS,upperCornerS,-1,0.33,nbSpheres,False, 0.5,seed=1) #"seed" make the "random" generation always the same
#sp.makeCloud(lowerCornerS,upperCornerS,-1,0.33,nbSpheres,False, 0.5,seed=1)

#sp.makeCloud(lowerCornerS,upperCornerS,rMean=.03,seed=1)
c1=pack.SpherePack([((0,0,0),.03),((.03,0,0),.03)])
sp.makeClumpCloud(lowerCornerS,upperCornerS,[c1],periodic=True,num=5000)
sp.toSimulation(material='spheres')


O.dt=.5*PWaveTimeStep()
# for the samples of this size (diametre is around 1.8cm), the O.dt is 8.8e-5 s, which is about a hundred times of 
# the small sample, and this leads to much faster calculation speed.

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller()

#######################################
###   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 = sigmaIso


print"___________________"
print"PHASE 1"
print "Settlement of particle under gravity only, fluid flow d"

triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 1
O.run(2000, True)

triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 2
O.run(2000, True)
triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 4

O.run(2000, True)
triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 7
O.run(2000, True)
triax.goal1 = triax.goal2 = triax.goal3 = sigmaIso / 10 * 10
O.run(2000, True)


stabilityThreshold = 0.001
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
print(O.iter)
#compactionFinished()


triax.depth0 = triax.depth
triax.height0 = triax.height
triax.width0 = triax.width
# triax.strain[0] = triax.strain[1] = triax.strain[2] = 0


# change control type: keep constant confinement in x,y, 20% compression in z
triax.stressMask = 0
triax.goal1 =  0.0001
triax.goal2 =  0.0001
triax.goal3 = -0.0002

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.