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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
