yade-users team mailing list archive

[Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>]

Hi Chiara,

Even with linear law, it was giving strange results, see modified script 
attached. You need to at least damp the cell deformation.
For HM, I don't see any obvious explanation of the problem. Can you get 
stable systems with similar properties, timestep, and with rigid boundaries?
I would double check the values of relative velocities and see if they 
have something wrong.

Good luck!

Bruno


On 27/07/10 11:21, chiara modenese wrote:
> 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
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to     : yade-users@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users
> More help   : https://help.launchpad.net/ListHelp
>    


--
_______________
Bruno Chareyre
Associate Professor
ENSE³ - Grenoble INP
Lab. 3SR
BP 53 - 38041, Grenoble cedex 9 - France
Tél : +33 4 56 52 86 21
Fax : +33 4 76 82 70 43
________________

#!/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.e5 # [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(),
	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.5),
	NewtonIntegrator(damping=0.3,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.save('temp');
#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()