yade-users team mailing list archive

[Sina Jafari <question247465@xxxxxxxxxxxxxxxxxxxxx>]

Question #247465 on Yade changed:
https://answers.launchpad.net/yade/+question/247465

Sina Jafari gave more information on the question:
here is my full script:

# -*- coding: utf-8 -*-
#*************************************************************************
#  Copyright (C) 2010 by Bruno Chareyre                                  *
#  bruno.chareyre_at_grenoble-inp.fr                                     *
#                                                                        *
#  This program is free software; it is licensed under the terms of the  *
#  GNU General Public License v2 or later. See file LICENSE for details. *
#*************************************************************************

## This script details the simulation of a triaxial test on sphere packings using Yade
## See the associated pdf file for detailed exercises
## the algorithms presented here have been used in published papers, namely:
## * Chareyre et al. 2002 (http://www.geosyntheticssociety.org/Resources/Archive/GI/src/V9I2/GI-V9-N2-Paper1.pdf)
## * Chareyre and Villard 2005 (https://yade-dem.org/w/images/1/1b/Chareyre&Villard2005_licensed.pdf)
## * Scholtès et al. 2009 (http://dx.doi.org/10.1016/j.ijengsci.2008.07.002)
## * Tong et al.2012 (http://dx.doi.org/10.2516/ogst/2012032)
##
## Most of the ideas were actually developped during my PhD.
## If you want to know more on micro-macro relations evaluated by triaxial simulations
## AND if you can read some french, it is here: http://tel.archives-ouvertes.fr/docs/00/48/68/07/PDF/Thesis.pdf

from yade import pack,plot,qt
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab
############################################
###   DEFINING VARIABLES AND MATERIALS   ###
############################################
key="K"
num_spheres=8000
psdSizes,psdCumm=[0.26*0.866,0.388*0.866,0.536*0.866,0.706*0.866,0.976*0.866,1.333*0.866,1.757*0.866,2.458*0.866,2.771*0.866,3.124*0.866,3.608*0.866,4.019*0.866,4.424*0.866,4.811*0.866,5.556*0.866,6.35*0.866],[0,1.865,3.657,6.12,9.25,14.4,18.88,29.18,35.45,42.4,52.2,61.857,71.716,80.898,90.522,100]
#targetPorosity = 0.387 #the porosity we want for the packing
compFricDegree = 26.5 # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 26.5 # contact friction during the deviatoric loading
targetPorosity=0.3 ### INJA BARAYE MAX ACN va MIN e
rate=0.0002 # loading rate (strain rate)
damp=0.25 # damping coefficient!!!!!!!!!!
stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
young=520e6 # contact stiffness!!!CHANGED!!!
mn,mx=Vector3(0,0,0),Vector3(44.25,44.25,44.25) # corners of the initial packing


## create materials for spheres and plates
O.materials.append(FrictMat(young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2000,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))

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

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack()
sp.particleSD2(radii=psdSizes,passing=psdCumm,numSph=8000,cloudPorosity=0.57)
O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
#walls=aabbWalls(material='walls',oversizeFactor=1)
#wallIds=O.bodies.append(walls)
#or alternatively (higher level function doing exactly the same):
#sp.toSimulation(material='spheres')

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

triax=TriaxialStressController(
	## ThreeDTriaxialEngine 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=True, # If true the confining pressure is generated by growing particles
)

newton=NewtonIntegrator(damping=damp)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_MindlinPhys()],
		[Law2_ScGeom_MindlinPhys_Mindlin()]
	),
	## 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=10,timestepSafetyCoefficient=0.7),
	triax,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses.dat'),
	qt.SnapshotEngine(fileBase="snap",iterPeriod=0,label='snapshoter'),	
	newton,
	

]

########################################
####   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=100000
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.01:
    break

O.save('confinedState'+'.yade.gz')
print "###      Isotropic state saved      ###"
print 'ACN=',utils.avgNumInteractions(),'Porosity=',utils.voxelPorosityTriaxial(triax),'Calculation Time(Sec)=',O.realtime

##############################
###   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)
setContactFriction(radians(finalFricDegree))

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

#set stress control on x and z, we will impose strain rate on y
triax.stressMask = 0
#now goal2 is the target strain rate
triax.goal2=-rate
# we define three lateral stresses during the test, here the same 10kPa as for the initial confinement.
triax.goal1=-rate
triax.goal3=-rate

#we can change damping here. What is the effect in your opinion?
#newton.damping=0.1

#Save temporary state in live memory. This state will be reloaded from the interface with the "reload" button.
O.saveTmp()
#####################################################
###    Example of how to record and plot data     ###
#####################################################

from yade import plot
import pylab
import numpy as np
import os

#########################################################################
#########################################################################
#########################################################################

class StressChecker():  
 dStress=nextStress=100000	
 def Sintrhisto(self):
		stress=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3		
		axis=2 
		ax1,ax2=(axis+0)%3,(axis+1)%3 
		angles,forces=[],[]
		for z in O.interactions:
			if not z.isReal: continue
			if z.id1<6 or z.id2<6: continue
			norm=z.geom.normal
			if norm[ax1]==0:
				angle=0
				force=z.phys.shearForce.norm()
			else:		
				angle=atan(norm[ax2]/norm[ax1])
				force=z.phys.shearForce.norm()
			angles.append(angle+pi/2)
			forces.append(force*10e-6)
		pylab.figure()
		values,bins=numpy.histogram(angles,weights=forces,bins=30)
		subp=pylab.subplot(111,polar=True)	
		pylab.bar(left=bins[:-1],height=values,width=np.pi/(1.05*30),alpha=.7,label=['xy'])
		pylab.xlabel('Shear Force Histogram-XY Plane')
		pylab.plot()
		pylab.savefig("interaction histogram-shear")
		FName = "interaction histogram-shear.png"
		NFName = "interaction histogram-shear-%.2f kPa.png"%float(stress/1000)
		os.rename(FName,NFName)

 def output(self):
		stress=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
		if abs(stress) > self.nextStress:
			self.nextStress += self.dStress
			self.Sintrhisto()

checker=StressChecker()
# include a periodic engine calling that function in the simulation loop
O.engines=O.engines[0:6]+[PyRunner(iterPeriod=20,command='checker.output()')]+O.engines[6:8]
#plot.plot()
#O.run(100,True)
#O.run(1000000)
####  PLAY THE SIMULATION HERE WITH "PLAY" BUTTON OR WITH THE COMMAND O.run(N)  #####

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