Compression Test

Hello Jabrane,

with prescribed strain rate, you only somehow prescribe boundary conditions
(not motion of all particles). Yade runs dynamic simulation, so (since
motion of all particles is not prescribed) there will always be stress
oscilations...

In what order is the oscilation? One solution is to prescribe lower strain
rate. In the limit case of infinitesimally small rate it should converge to
the static simulation with no oscilation..

cheers
Jan

PS: if you precribe motion of all particles, the stress oscilation would
vanish. But in that case there is not much to compute..

2015-11-16 11:07 GMT+01:00 Yor1 <email address hidden>:

> New question #274242 on Yade:
> https://answers.launchpad.net/yade/+question/274242
>
> Hello,
>
> I try to calculate the balance energy in a compression test.
> The problem that i have oscillations in the stress despite i fixed a
> strain rate.
> The confinement is equal to 1 MPa.
> This is my code.
>
> from yade import ymport, utils , plot
> import math
>
> #---------------- SIMULATIONS DEFINED HERE (assembly, material, boundary
> conditions)
>
> #### packing (previously constructed)
> PACKING='X1Y2Z1_20k'
> OUT=PACKING+'_compressionTest_5MPa_r002_frictWall'
>
> #### Simulation Control
> DAMP=0.4 # numerical damping
> saveData=100 # data record interval
> iterMax=2000000 # maximum number of iteration (to be adjusted)
> saveVTK=1000 # Vtk files record interval
>
> #### Boundary Conditions
> confinement=-1e6
> strainRate=-0.01
>
> #### Material microproperties -> Lac du Bonnet granite (cf. A DEM model
> for soft and hard rock, Scholtes & Donze, JMPS 2013)
> intR=1.4464 # allows near neighbour interaction and coordination number
> K=13 (determined with coordinationNumber.py -> to be adjusted for every
> packing)
> DENS=4000 # could be adapted to match material density:
> dens_DEM=dens_rock*(V_rock/V_particles)=dens_rock*1/(1-poro_DEM) -> poro?
> YOUNG=65e9
> FRICT=10
> ALPHA=0.4
> TENS=8e6
> COH=160e6
>
> #### material definition
> def sphereMat(): return JCFpmMat(type=1,density=DENS,young=YOUNG,poisson =
> ALPHA,frictionAngle=radians(FRICT),tensileStrength=TENS,cohesion=COH)
> def wallMat(): return
> JCFpmMat(type=0,density=DENS,young=YOUNG,frictionAngle=radians(0))
>
> #### preprocessing to get dimensions
>
> O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
>
> dim=utils.aabbExtrema()
> xinf=dim[0][0]
> xsup=dim[1][0]
> X=xsup-xinf
> yinf=dim[0][1]
> ysup=dim[1][1]
> Y=ysup-yinf
> zinf=dim[0][2]
> zsup=dim[1][2]
> Z=zsup-zinf
>
> R=0
> Rmax=0
> numSpheres=0.
> for o in O.bodies:
> if isinstance(o.shape,Sphere):
> numSpheres+=1
> R+=o.shape.radius
> if o.shape.radius>Rmax:
> Rmax=o.shape.radius
> Rmean=R/numSpheres
>
> O.reset() # all previous lines were for getting dimensions of the packing
> to create walls at the right positions (below) because walls have to be
> genrated after spheres for FlowEngine
>
> #### now we construct the surrounding walls with right dimensions
> ### walls
>
> mn,mx=Vector3(xinf+0.1*Rmean,yinf+0.1*Rmean,zinf+0.1*Rmean),Vector3(xsup-0.1*Rmean,ysup-0.1*Rmean,zsup-0.1*Rmean)
>
> walls=utils.aabbWalls(oversizeFactor=1.5,extrema=(mn,mx),thickness=min(X,Y,Z)/100.,material=wallMat)
> wallIds=O.bodies.append(walls)
>
> ### packing
>
> beam=O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))
>
> g=Vector3(0,0,-9.8) ###gravitational acceleration
>
> #---------------- ENGINES DEFINED HERE
>
> #### simulation is defined here (DEM loop, interaction law, servo control,
> recording, etc...)
> O.engines=[
> ForceResetter(),
>
> InsertionSortCollider([Bo1_Box_Aabb(),Bo1_Sphere_Aabb(aabbEnlargeFactor=intR,label='Saabb')]),
> InteractionLoop(
>
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=intR,label='SSgeom'),Ig2_Box_Sphere_ScGeom()],
>
> [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
>
> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT,label='interactionLaw')]
> ),
> TriaxialStressController(internalCompaction=False,label='triax'),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4,
> defaultDt=0.1*utils.PWaveTimeStep()),
> NewtonIntegrator(gravity=g,damping=DAMP,label="newton"),
>
> PyRunner(iterPeriod=int(saveData),initRun=True,command='recorder()',label='data'),
>
> VTKRecorder(iterPeriod=int(saveVTK),initRun=True,fileName=OUT+'-',recorders=['spheres','jcfpm','cracks'],Key=OUT,label='vtk')
> ]
>
> plot.plots={'i':('s1','s2','s3')}
> plot.plot()
>
>
> O.step();
>
> SSgeom.interactionDetectionFactor=-1.
> Saabb.aabbEnlargeFactor=-1.
>
> numSSlinks=0
> numCohesivelinks=0
> numFrictionalLinks=0
> for i in O.interactions:
> if not i.isReal : continue
> if isinstance(O.bodies[i.id1].shape,Sphere) and
> isinstance(O.bodies[i.id2].shape,Sphere):
> numSSlinks+=1
> if i.phys.isCohesive :
> numCohesivelinks+=1
> else :
> numFrictionalLinks+=1
>
> print "nbSpheres=", numSpheres," | coordination number =",
> 2.0*numCohesivelinks/numSpheres
>
>
> #### APPLYING ISOTROPIC LOADING
> triax.stressMask=7
> triax.goal1=confinement
> triax.goal2=confinement
> triax.goal3=confinement
> triax.max_vel=0.001
>
> while 1:
> if confinement==0:
> O.run(1000,True) # to stabilize the system
> break
> O.run(100,True)
> unb=unbalancedForce()
> #note: triax.stress(k) returns a stress vector, so we need to keep only
> the normal component
>
> meanS=abs(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<0.005 and abs(meanS-abs(confinement))/abs(confinement)<0.001:
> O.run(1000,True) # to stabilize the system
> e10=triax.strain[0]
> e20=triax.strain[1]
> e30=triax.strain[2]
> break
>
> O.bodies[wallIds[2]].mat.frictionAngle=radians(30)
> O.bodies[wallIds[3]].mat.frictionAngle=radians(30)
>
> triax.stressMask=7
> triax.goal1=confinement
> triax.goal2=strainRate
> triax.goal3=confinement
> triax.max_vel=1
>
> O.run(iterMax)
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi Jan,

In fact when i delete the gravity from the code
NewtonIntegrator(gravity=g,damping=DAMP,label="newton") ------> NewtonIntegrator(damping=DAMP,label="newton")

it runs well and i don't have the oscillation.

The question is how can i integrate the gravity in the code (without having oscillations) to calculate potential energy.

Jabrane.

Hi Jabrane,

well, this is a good situation in my point of view :-) comming to our
previous discussions: Is gravity really important in your compression test?
If not (my opinion), just don't introduce it. The work done by gravity (for
normal materials) is (or at least should be) negligible to other kinds of
energies and works.

Concerning energies. Even if you run explicit dynamic simulation, I would
say experimentally the situation is almost static. Therefore the
contribution of kinetic energy should also not be very high. Usually also
such test is simulated without gravity (as it has almost zero effect).

To get feeling about different energies, you can create a mental
experiment, just one mass point (with mass m) on a vertical spring
(stiffness k) oscilating with certain amplitude "a" around equilibrium
position. After a short simple math you ends with maximum difference in
gravity potantial energy 2*m*g*a, maximum difference in potential elastic
energy of the spring k*a^2 and maximum kinetic energy (1/2)*k*a^2. For the
same amplitude a and for increasing stiffness, you wil get less and less
contribution of gravity. Normaly for triaxial test of usual materials,
"stiffness is much higher than displacements" and gravity basicle has no
effect.

Similar experiment you can do for static loading applying force F=k*u, u
being spring elongation. Then the change of gravity potential energy is
m*g*u and change of spring elastic potential energy is (1/2)*k*u^2. Again,
for higher stiffness the contribution of gravity derceases..

cheers
Jan

2015-11-16 13:46 GMT+01:00 Yor1 <email address hidden>:

> Question #274242 on Yade changed:
> https://answers.launchpad.net/yade/+question/274242
>
> Yor1 posted a new comment:
> Hi Jan,
>
> In fact when i delete the gravity from the code
> NewtonIntegrator(gravity=g,damping=DAMP,label="newton") ------>
> NewtonIntegrator(damping=DAMP,label="newton")
>
> it runs well and i don't have the oscillation.
>
> The question is how can i integrate the gravity in the code (without
> having oscillations) to calculate potential energy.
>
> Jabrane.
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Assigning gravity for a compression test is a little bit strange. If it plays a role there will be a significant stress gradient and thus the results can't be interpreted properly: the 1MPa confinment is irrelevant since there is not a unique state of stress.

I want to simulate a pillar in underground mines.
Also, i want to study the energetic behavior to predict the rockburst.
In YADE, i develop a script to calculate the release of energy. At this current stage i want to validate my calcul of energy
compared to UDEC with a simple case.

Jabrane.