yade-users team mailing list archive

[jessica sjah <question216096@xxxxxxxxxxxxxxxxxxxxx>]

New question #216096 on Yade:
https://answers.launchpad.net/yade/+question/216096

Dear yade user,

Could anyone please help me to see my code whether is it correct or not to make a biaxial undrained test in 2D (without periodic condition)?

I'm using "ThreeDTriaxialEngine", so I can use "strainRate1=strainRate2" to avoid the volume change (Please correct me if my idea is wrong).

Is there anyone have an experience in doing the triaxial/biaxial undrained test? What is the best way to model it?

Thank you for your help and your sharing.

Here is my code:

# -*- coding: utf-8 -*-
from yade import pack,log, plot
from utils import *

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

num_spheres		= 500			# number of spheres
compFricDegree 		= 30 			# contact friction during the confining phase
finalFricDegree 	= 30 			# contact friction during the deviatoric loading
rate			= 0.01 			# loading rate (strain rate)
damp			= 0.1 			# damping coefficient
stabilityThreshold	= 0.001 		# unbalancedForce
key			= '_biax_undrained_' 	# simulation's name 
young			= 5e6			# contact stiffness
mn,mx=Vector3(0,0,0),Vector3(1,1,1) 		# corners of the initial packing
thick 			= 0.01 			# thickness of the plates

## 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=utils.aabbWalls([mn,mx],thickness=thick,material='walls')
wallIds=O.bodies.append(walls)

## generate a random loose particles packing
sp=pack.SpherePack()
sp.makeCloud(Vector3(0,0,0.5),Vector3(1,1,0.5) ,-1,0.3333,num_spheres,False, 0.95)

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

# Blocked certain degress of freedom to make 2D-Model in plane-XY
for k in O.bodies:
 if isinstance(k.shape, Sphere): k.state.blockedDOFs='zXY'

# initial timestep
O.dt=.5*utils.PWaveTimeStep() 
O.usesTimeStepper=True

############################
###   DEFINING ENGINES   ###
############################

triax=ThreeDTriaxialEngine(
	## ThreeDTriaxialEngine will be used to control stress and strain. It controls particles size and plates positions.
	maxMultiplier=1.005, # spheres growing factor (fast growth)
	finalMaxMultiplier=1.002, # spheres growing factor (slow growth)
	thickness = thick,
	stressControl_1 = False, #switch stress/strain control
	stressControl_2 = False,
	## The stress used for (isotropic) internal compaction
	sigma_iso = 10000,
	## Independant stress values for anisotropic loadings
	sigma1=10000,
	sigma2=10000,
	internalCompaction=True, # If true the confining pressure is generated by growing particles
	Key=key, # passed to the engine so that the output file will have the correct name
)

newton=NewtonIntegrator(damping=damp)
unb=unbalancedForce()

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()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
	triax,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),
	newton,
        PyRunner(command='addPlotData()',iterPeriod=100)
]

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

def addPlotData():
   plot.addData(unbalanced=unb,i=O.iter,
      s_right=triax.stress(triax.wall_right_id)[0],s_top=triax.stress(triax.wall_top_id)[1],s_front=triax.stress(triax.wall_front_id)[2],
      s_left=triax.stress(triax.wall_left_id)[0],s_bot=triax.stress(triax.wall_bottom_id)[1],s_back=triax.stress(triax.wall_back_id)[2],
      exx=triax.strain[0],eyy=triax.strain[1],ezz=triax.strain[2],axial_strain=triax.strain[1],
      dev_stress=triax.stress(triax.wall_top_id)[1]-triax.stress(triax.wall_right_id)[0], 
      conf_stress=triax.stress(triax.wall_right_id)[0],
      vol_strain= triax.strain[0]+triax.strain[1]
   )

# define what to plot
plot.plots={'i':('unbalanced'),'i ':('s_right','s_top','s_front','s_left','s_bot','s_back'),' i':('exx','eyy','ezz'),
   ' axial_strain':('dev_stress', 'conf_stress'),' axial_strain ':('vol_strain'),
}
# show the plot
plot.plot()

#Display spheres with 2 colors for seeing rotations better
Gl1_Sphere.stripes=1
yade.qt.Controller(), yade.qt.View()

#######################################
###   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
  if unb<stabilityThreshold and abs(meanS-triax.sigma_iso)/triax.sigma_iso<0.001:
    break

O.save('confinedState'+key+'.yade.gz')
print "###      Isotropic 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)
triax.setContactProperties(finalFricDegree)

#set independant stress control on each axis
triax.stressControl_1=triax.stressControl_2=True

# We turn all these flags true, else boundaries will be fixed
triax.wall_bottom_activated=True
triax.wall_top_activated=True
triax.wall_left_activated=True
triax.wall_right_activated=True
triax.wall_back_activated=True
triax.wall_front_activated=True

#If we want a triaxial loading at imposed strain rate, let's assign srain rate instead of stress
##...  we are tired of typing "True" and "False", we use implicit conversion from integer to boolean
triax.stressControl_2=0
##... fix a maximum strain rate to go progressivly to the desired stress state in direction 2
triax.strainRate1=0.01
triax.strainRate2=0.01

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.