yade-users team mailing list archive

[chiara modenese <c.modenese@xxxxxxxxx>]

Hi guys (Bruno ?),

I ask you a bit of help/advice. I am testing the PBC with HM law. I attach a
simple script. With the linear elastic law everything works fine. With HM,
when the wanted isotropic stress level is about to be reached, I see
particles going crazy and kinetic energy terribly increasing. Pleas just run
the script to reproduce it (look at qt.View and at the plots as well). On
the other hand, if I block particle rotations it works fine (just set
block_particles_rot=True in the py script and run it again). Why? Any
advice? I tried to play with time step, strainRate and so on but I could get
a stable solution.

Thanks for help!

Chiara

#!/usr/local/bin/yade-trunk
# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-

"""
@Note:
It seems there is a problem with particle rotations. If rotations are fixed and damping (local or global) is present, then everything seems fine.
Otherwise something is wrong. In the contact model? In PBC?
"""

from matplotlib import rc
rc('font',**{'family':'serif','serif':['Times']})
rc('text', usetex=True)

# Set PBC
O.periodic=True

from math import *
# Script to perform a typical triaxial test
#-----------------------------------------------------
# simulation parameters
#-----------------------------------------------------
# get the size of the cube with a higher initial porosity, then compact the sample
avg_radius		= 5e-3 # [m]
number_of_spheres	= 500
V_ball 			= 4./3.*pi*pow(avg_radius,3)  # [m^3]
poro_init		= 0.90
delta 			= pow(number_of_spheres*V_ball/((1-poro_init)),1./3.)
#-----------------------------------------------------
O.cell.refSize		= (delta,delta,delta)
young_spheres		= 1.e4 # [Pa]
local_poisson		= 0.25
frictionAngle		= atan(0.5) # friction angle [rad]
std_dev_radius		= 0.2
betan_damp		= 0.2 # contact normal damping coeff (HM only)
betas_damp		= 0.2 # contact shear damping coeff (HM only)
block_particles_rot	= False
global_damp		= 0.0 # Cundall's non-viscous damping (use only if no contact damping is present)
#-----------------------------------------------------
# material
#-----------------------------------------------------
O.materials.append(FrictMat(    young=young_spheres,
				poisson=local_poisson,
				frictionAngle=frictionAngle, # [rad]
				density=1.e3,
				label='granular'))
#-----------------------------------------------------
# geometry
#-----------------------------------------------------
from yade import pack
from yade import utils
# add the cloud of spheres
sphere_cloud=pack.SpherePack()
sphere_cloud.makeCloud( minCorner=(0,0,0),
			maxCorner=O.cell.refSize,
			rMean=avg_radius,
			rRelFuzz=std_dev_radius,
			num=number_of_spheres,
			periodic=True,
			porosity=-1)
O.bodies.append([utils.sphere(	center,
				rad,
				material='granular')
				for center,rad in sphere_cloud])
#-----------------------------------------------------
# list of engines
#-----------------------------------------------------
O.engines=[
	ForceResetter(),
	BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InsertionSortCollider(nBins=5,sweepLength=avg_radius*0.05),
	InteractionDispatchers(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()], 

		#----------------------- linear spring law 
		#[Ip2_FrictMat_FrictMat_FrictPhys()],
		#[Law2_ScGeom_FrictPhys_Basic(label='contactLaw',traceEnergy=False)]

		#----------------------- non linear HM law
		[Ip2_FrictMat_FrictMat_MindlinPhys()],
		[Law2_ScGeom_MindlinPhys_Mindlin(betan=betan_damp,betas=betas_damp,useDamping=True,label='contactLaw')]
	),
	GlobalStiffnessTimeStepper(timeStepUpdateInterval=10),
	PeriTriaxController(dynCell=True,mass=0.3,maxUnbalanced=0.01,relStressTol=0.02,goal=[-1e4,-1e4,-1e4],stressMask=7,globUpdate=5,maxStrainRate=(.1,.1,.1),doneHook='triaxDone()',reversedForces=True,label='triax',growDamping=0.),
	NewtonIntegrator(damping=global_damp,homotheticCellResize=2),
	PeriodicPythonRunner(iterPeriod=100,command='saveData()')
]
#-----------------------------------------------------
# kinematic constrain
#-----------------------------------------------------
def block_rot():
	for b in O.bodies:
		b.state.blockedDOFs=['rx','ry','rz'] # block particles rotations (to increase the peak strenght)
if (block_particles_rot==True): block_rot(); # eventually block rotations
O.dt=.1*utils.PWaveTimeStep() # initial timestep, then GSTS will be used

from yade import qt
qt.Controller(); qt.View()
#-----------------------------------------------------
# plots and other functions
#-----------------------------------------------------
from math import *
from yade import plot

plot.plots={'Iter':(('sigma1','r'),('sigma2','m'),('sigma3','b')),'Time':('Unbalanced'),'I':(('strain1','r'),('strain2','m'),('strain3','b')),'step':('Ek')}

def saveData():
	plot.addData(sigma1=triax.stress[0],sigma2=triax.stress[1],sigma3=triax.stress[2], # stresses
		     strain1=triax.strain[0],strain2=triax.strain[1],strain3=triax.strain[2], # strains
		     Iter=O.iter,Time=O.time,I=O.iter,step=O.iter,
		     Unbalanced=triax.currUnbalanced, # unbalanced force
		     Ek=utils.kineticEnergy()) # kinetic energy

plot.live=True; plot.plot()
O.run();

phase=0
def triaxDone():
	global phase
	if phase==0: # isotropic phase
		print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
		print 'Now maintain constant lateral stresses and increase epsilon z'
		triax.goal=(-1e4,-1e4,-1.5)
		triax.stressMask=3 # strain controll in the vertical direction
		phase+=1
	elif phase==1:
		print 'Here we are: stress',triax.stress,'strain',triax.strain,'stiffness',triax.stiff
		print 'Done, pausing now.'
		O.pause()