yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07875
[Branch ~yade-dev/yade/trunk] Rev 2921: - Add material from the "discrete mechanics" school in Grenoble 2011.
------------------------------------------------------------
revno: 2921
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Wed 2011-09-14 16:14:26 +0200
message:
- Add material from the "discrete mechanics" school in Grenoble 2011.
added:
scripts/triax-tutorial/
scripts/triax-tutorial/README
scripts/triax-tutorial/script-session1.py
scripts/triax-tutorial/script-session2.py
scripts/triax-tutorial/session1-basic-triax.pdf
scripts/triax-tutorial/session2-stress-probing.pdf
scripts/triax-tutorial/triax-table
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== added directory 'scripts/triax-tutorial'
=== added file 'scripts/triax-tutorial/README'
--- scripts/triax-tutorial/README 1970-01-01 00:00:00 +0000
+++ scripts/triax-tutorial/README 2011-09-14 14:14:26 +0000
@@ -0,0 +1,4 @@
+The material in this folder was used for the numerical sessions of Alert's "Olek Zinkiewicz" course on discrete mechanics, Grenoble, 2011. They show basic simulation of triaxial tests (session 1) and advanced boundary conditions for stress-strain probing (session 2).
+Authors: Nejib Hada, Luc Sibille, Bruno Chareyre.
+Yade's version should be bzr2830 or later, although some earlier versions could work as well.
+
\ No newline at end of file
=== added file 'scripts/triax-tutorial/script-session1.py'
--- scripts/triax-tutorial/script-session1.py 1970-01-01 00:00:00 +0000
+++ scripts/triax-tutorial/script-session1.py 2011-09-14 14:14:26 +0000
@@ -0,0 +1,188 @@
+# -*- coding: utf-8 -*-
+from yade import pack,log
+from utils import *
+############################################
+### DEFINING VARIABLES AND MATERIALS ###
+############################################
+
+# The following 5 lines will be used later for batch execution
+nRead=utils.readParamsFromTable(
+ num_spheres=500,# number of spheres
+ compFricDegree = 30, # contact friction during the confining phase
+ unknownOk=True
+)
+from yade.params import table
+
+num_spheres=table.num_spheres# number of spheres
+compFricDegree = table.compFricDegree # contact friction during the confining phase
+finalFricDegree = 30 # contact friction during the deviatoric loading
+rate=0.02 # loading rate (strain rate)
+damp=0.2 # damping coefficient
+stabilityThreshold=0.01 # we test unbalancedForce against this value in different loops (see below)
+key='_triax_base_' # put you simulation's name here
+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)
+
+## use a SpherePack object to generate a random loose particles packing
+sp=pack.SpherePack()
+sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95)
+
+## approximate mean rad of the futur dense packing for latter use (as an exercise: you can compute exact
+volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
+mean_rad = pow(0.09*volume/num_spheres,0.3333)
+
+
+clumps=False
+if clumps:
+ c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
+ sp.makeClumpCloud((0,0,0),(1,1,1),[c1],periodic=False)
+ O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
+ standalone,clumps=sp.getClumps()
+ for clump in clumps:
+ O.bodies.clump(clump)
+ for i in clump[1:]: O.bodies[i].shape.color=O.bodies[clump[0]].shape.color
+ #sp.toSimulation()
+else:
+ O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
+
+O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away
+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.+2e4/young, # spheres growing factor (fast growth)
+ finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
+ thickness = thick,
+ stressControl_1 = False, #switch stress/strain control
+ stressControl_2 = False,
+ stressControl_3 = False,
+ ## The stress used for (isotropic) internal compaction
+ sigma_iso = 10000,
+ ## Independant stress values for anisotropic loadings
+ sigma1=10000,
+ sigma2=10000,
+ sigma3=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)
+
+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
+]
+
+#Display spheres with 2 colors for seeing rotations better
+Gl1_Sphere.stripes=0
+if nRead==0: 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=triax.stressControl_3=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
+#triax.stressControl_2=0 #we are tired of typing "True" and "False", we use implicit conversion from integer to boolean
+#triax.strainRate2=rate
+#triax.strainRate1=100*rate
+#triax.strainRate3=100*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()
+
+#####################################################
+### Exemple of how to record and plot data ###
+#####################################################
+
+#from yade import plot
+
+### a function saving variables
+#def history():
+ #plot.addData(e11=triax.strain[0], e22=triax.strain[1], e33=triax.strain[2],
+ #s11=triax.stress(triax.wall_right_id)[0],
+ #s22=triax.stress(triax.wall_top_id)[1],
+ #s33=triax.stress(triax.wall_front_id)[2],
+ #i=O.iter)
+
+#if 1:
+ ## include a periodic engine calling that function in the simulation loop
+ #O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
+ ##O.engines.insert(4,PyRunner(iterPeriod=20,command='history()',label='recorder'))
+#else:
+ ## With the line above, we are recording some variables twice. We could in fact replace the previous
+ ## TriaxialRecorder
+ ## by our periodic engine. Uncomment the following line:
+ #O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')
+
+#O.run(100,True)
+
+### declare what is to plot. "None" is for separating y and y2 axis
+#plot.plots={'i':('e11','e22','e33',None,'s11','s22','s33')}
+
+##display on the screen (doesn't work on VMware image it seems)
+##plot.plot()
+
+## In that case we can still save the data to a text file at the the end of the simulation, with:
+#plot.saveDataTxt('results'+key)
+##or even generate a script for gnuplot. Open another terminal and type "gnuplot plotScriptKEY.gnuplot:
+#plot.saveGnuplot('plotScript'+key)
=== added file 'scripts/triax-tutorial/script-session2.py'
--- scripts/triax-tutorial/script-session2.py 1970-01-01 00:00:00 +0000
+++ scripts/triax-tutorial/script-session2.py 2011-09-14 14:14:26 +0000
@@ -0,0 +1,228 @@
+# -*- coding: utf-8 -*-
+from yade import pack,log
+from utils import *
+
+num_spheres=500
+## corners of the initial packing
+mn,mx=Vector3(0,0,0),Vector3(1,1,1)
+thick = 0.01
+compFricDegree = 2
+rate=0.2
+damp=0.1
+stabilityThreshold=0.001
+key='_define_a_name_'
+
+## create material #0, which will be used as default
+O.materials.append(FrictMat(young=5e6,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
+O.materials.append(FrictMat(young=5e6,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)
+
+sp=pack.SpherePack()
+sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95)
+
+volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
+mean_rad = pow(0.09*volume/num_spheres,0.3333)
+
+clumps=False
+if clumps:
+ c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
+ sp.makeClumpCloud((0,0,0),(1,1,1),[c1],periodic=False)
+ O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
+ standalone,clumps=sp.getClumps()
+ for clump in clumps:
+ O.bodies.clump(clump)
+ for i in clump[1:]: O.bodies[i].shape.color=O.bodies[clump[0]].shape.color
+ #sp.toSimulation()
+else:
+ O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])
+
+O.dt=.5*utils.PWaveTimeStep() # initial timestep, to not explode right away
+O.usesTimeStepper=True
+
+triax=ThreeDTriaxialEngine(
+ maxMultiplier=1.005,
+ finalMaxMultiplier=1.002,
+ thickness = thick,
+ stressControl_1 = False,
+ stressControl_2 = False,
+ stressControl_3 = False,
+ ## The stress used for (isotropic) internal compaction
+ sigma_iso = 10000,
+ ## Independant stress values for anisotropic loadings
+ sigma1=10000,
+ sigma2=10000,
+ sigma3=10000,
+ internalCompaction=True,
+ Key=key,
+)
+
+newton=NewtonIntegrator(damping=damp)
+
+O.engines=[
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()],nBins=-1,verletDist=-mean_rad*0.06,binCoeff=2),
+ 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
+]
+
+#Display spheres with 2 colors for seeing rotations better
+Gl1_Sphere.stripes=0
+yade.qt.Controller(), yade.qt.View()
+
+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('compressedState'+key+'.xml')
+print "### Isotropic state saved ###"
+
+#let us turn internal compaction off...
+triax.internalCompaction=False
+
+#
+triax.setContactProperties(30)
+
+#... and make stress control independant on each axis
+triax.stressControl_1=triax.stressControl_2=triax.stressControl_3=True
+# We have to 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
+triax.stressControl_2=0 #we are tired of typing "True" and "False", we use implicit conversion from integer to boolean
+triax.strainRate2=0.01
+triax.strainRate1=triax.strainRate3=1000.0
+
+#Else if we want imposed stress increments, etc...
+
+
+##First perform a transverse isotropic compression (or a stress controlled drained triaxial compression) to reach an initial state from where stress probes will be applied
+##... need to active stress control in 3 directions
+#triax.stressControl_1=triax.stressControl_2=triax.stressControl_3=True
+##... choose the value of axial stress where we want to stop the compression
+#triax.sigma2=12000
+##... fix a maximum strain rate to go progressivly to the desired stress state in direction 2
+#triax.strainRate2=0.01
+##... fix a high value of maximum strain rate in radial direction to be sure to keep in any conditions a constant confining pressure
+#triax.strainRate1=triax.strainRate3=1000.0
+
+
+#while 1:
+ #O.run(100, True)
+ ##the global unbalanced force on dynamic bodies, thus excluding boundaries, which are not at equilibrium
+ #unb=unbalancedForce()
+ ##note: triax.stress(k) returns a stress vector, so we need to keep only the normal component
+ #axialS=triax.stress(triax.wall_top_id)[1]
+ #print 'unbalanced force:',unb,' sigma2: ',axialS
+ #if unb<stabilityThreshold and abs(axialS-triax.sigma2)/triax.sigma2<0.001:
+ #break
+
+#O.save('anisotropicState'+key+'.xml')
+
+
+##Perform stress probes from the anisotropic state
+
+#dSnorm = 100.0 #norm of the stress increment
+#nbProbes = 16 #number of stress directions tested
+#rampIte = 20 #nb iterations to increase the stress state until the final desired stress value
+##LUC: je fixe des nombres d'iterations c'est moins elegant que de chercher explicitement un etat d'equilibre mais ca permet de poursuivre le calcul meme si un etat de contrainte n'est pas correctement atteint pour un stress probe et qu'il est difficile de se stabiliser a cet etat de contrainte (i.e. attendre longtemps...)
+#stabIte = 5000 #nb iterations to stabilize sample after reaching the final stress value
+
+
+## an array for saving stress increments and strain responses; arrays are in "numpy" extension
+#import numpy
+#probings=numpy.zeros((3,nbProbes))
+
+#def increment(dsr=0,dsa=1):
+ #for ite in range(rampIte):# progressivaly increase of stress state
+ #O.run(20, True)
+ ##incrementation of stress state
+ #triax.sigma2 = initSa+dsa/rampIte*ite
+ #triax.sigma1 = triax.sigma3 = initSr+dsr/rampIte*ite
+ #print triax.sigma1, triax.sigma2
+
+ ## fix the stress value for stabilization at the final state
+ #triax.sigma2 = initSa+dsa
+ #triax.sigma1 = triax.sigma3 = initSr+dsr
+
+ #while 1:
+ #O.run(100, True)
+ #unb=unbalancedForce()
+ #print 'unbalanced force:',unb,' strain: ',triax.strain
+ #if unb<stabilityThreshold: break
+
+
+
+## loop over all the stress directions
+#for i in range(nbProbes):
+
+ ## computation of the stress direction of the current stress probe
+ #alphaS = 2*pi/nbProbes*(i-1)
+ #print 'stress probe nb:',i,' stress direction (deg): ',degrees(alphaS)
+
+ ## computation of the stress increment in the axial direction
+ #dSa = dSnorm*sin(alphaS)
+
+ ## computation of the stress increment in the radial direction
+ #dSr = dSnorm*cos(alphaS)/sqrt(2.0)
+
+ ##Load the initial anisotropic state before running a new stress probe
+ #O.load('anisotropicState'+key+'.xml')
+ ##We redefine the "triax" label, else it would point to inactive engine from previous simulation that is still in memory
+ #triax=O.engines[4]
+
+ #initSa=triax.sigma2 #save of the initial axial stress
+ #initSr=triax.sigma1 #save of the initial radial stress
+
+ ## define the final stress state to be reached
+ #finalSa = initSa+dSa
+ #finalSr = initSr+dSr
+
+ ##... need to active stress control in 3 directions if not yet done
+ #triax.stressControl_1=triax.stressControl_2=triax.stressControl_3=True
+
+ ## fix a high value of maximum strain rate, the progressive loading will be done by progressively increasing the desired stress state at each iteration
+ #triax.strainRate1=triax.strainRate2=triax.strainRate3=1000.0
+
+ #increment(dSr,dSa)
+ #probings[:,i]=triax.strain
+
+
+
+##open a file for writing probing results
+#a=open('probings'+key,'w')
+#for i in range(nbProbes): a.write(str(probings[0][i])+' '+str(probings[1][i])+' '+str(probings[2][i])+'\n')
+#a.close()
+
+##plot
+#from pylab import *
+#plot(probings[0,:]*sqrt(2),probings[1,:],'bo--')
+#ylabel(r'$\epsilon_{22}$')
+#xlabel(r'$\epsilon_{11} \times \sqrt{2}$')
+#title('response envelop')
+#grid()
+
+#show()
+
=== added file 'scripts/triax-tutorial/session1-basic-triax.pdf'
Binary files scripts/triax-tutorial/session1-basic-triax.pdf 1970-01-01 00:00:00 +0000 and scripts/triax-tutorial/session1-basic-triax.pdf 2011-09-14 14:14:26 +0000 differ
=== added file 'scripts/triax-tutorial/session2-stress-probing.pdf'
Binary files scripts/triax-tutorial/session2-stress-probing.pdf 1970-01-01 00:00:00 +0000 and scripts/triax-tutorial/session2-stress-probing.pdf 2011-09-14 14:14:26 +0000 differ
=== added file 'scripts/triax-tutorial/triax-table'
--- scripts/triax-tutorial/triax-table 1970-01-01 00:00:00 +0000
+++ scripts/triax-tutorial/triax-table 2011-09-14 14:14:26 +0000
@@ -0,0 +1,4 @@
+!OMP_NUM_THREADS description num_spheres compFricDegree
+1 sim1 500 4
+1 sim2 500 8
+1 sim3 500 16
\ No newline at end of file
