yade-users team mailing list archive

[Jan Stránský <question274242@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Open => Answered

Jan Stránský proposed the following answer:
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 <question274242@xxxxxxxxxxxxxxxxxxxxx>:

> 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     : yade-users@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-users
> More help   : https://help.launchpad.net/ListHelp
>

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.