plastic dissipation

Hello,

I use "Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True,label='plastic')"
> to calculate the plastic dissipation but i have a message of error

> "FATAL /home/hamdi/YADE/sources/trunk/pkg/common/InteractionLoop.cpp:133
> action: None of given Law2 functors can handle interaction #2+931, types
> geom:ScGeom=1 and phys:JCFpmPhys=9 (LawDispatcher::getFunctor2D returned
> empty functor)"

As the error says, you use JCFpmMat and JCFpmPhys,
but Law2_ScGeom_FrictPhys_CundallStrack expects FrictPhys. Since you have
no other Law2 defined, "
 None of given Law2 functors can handle interaction" with JCFpmPhys

Two solutions
1) Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM (I am not sure how it
is with energy tracking)
2) Use FrictMat, Ip2_FrictMat_FrictMat_FrictPhys,
Law2_ScGeom_FrictPhys_CundallStrack

cheers
Jan

2015-12-03 13:47 GMT+01:00 Yor1 <email address hidden>:

> New question #276604 on Yade:
> https://answers.launchpad.net/yade/+question/276604
>
> Hello,
>
> I want to calculate the plastic dissipation in triaxial test.
> I use
> "Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True,label='plastic')" to
> calculate the plastic dissipation but i have a message of error
>
> "FATAL /home/hamdi/YADE/sources/trunk/pkg/common/InteractionLoop.cpp:133
> action: None of given Law2 functors can handle interaction #2+931, types
> geom:ScGeom=1 and phys:JCFpmPhys=9 (LawDispatcher::getFunctor2D returned
> empty functor)"
>
> 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'
>
> #### Simulation Control
> DAMP=0.4 # numerical damping
> saveData=100 # data record interval
> iterMax=60000 # maximum number of iteration (to be adjusted)
> saveVTK=10000 # Vtk files record interval
>
> #### Boundary Conditions
> confinement=-1e6
> #uniaxial_stress=-1e6
> delta_stress=-1e6
> stress_max=-200e6
> 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))
>
> ### set a color to the spheres
> for o in O.bodies:
> if isinstance(o.shape,Sphere):
> o.shape.color=(0.7,0.5,0.3)
>
>
> #---------------- 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')],
>
> [Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy=True,label='plastic')]
> ),
> TriaxialStressController(internalCompaction=False,label='triax'),
>
> GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.4,
> defaultDt=0.1*utils.PWaveTimeStep()),
> NewtonIntegrator(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()
>
> #---------------- SIMULATION STARTS HERE
>
> #### manage interaction detection factor during the first timestep and
> then set default interaction range ((cf. A DEM model for soft and hard
> rock, Scholtes & Donze, JMPS 2013))
> O.step();
> ### initializes the interaction detection factor
> SSgeom.interactionDetectionFactor=-1.
> Saabb.aabbEnlargeFactor=-1.
>
> #### coordination number verification and reinforcement of boundary
> particles
> 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.01
>
> 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.save(OUT+'_isotropicState.yade.gz')
>
> #### APPLYING DEVIATORIC LOADING
>
> #### !!! if you want to block boundary particles to simulate extreme
> friction on platens -> not sure that it works!
> #for o in O.bodies:
> #if isinstance(o.shape,Sphere):
> #if (o.state.pos[1])<(yinf+2*Rmean) or (o.state.pos[1])>(ysup-2*Rmean):
> #o.state.blockedDOFs+='xz'
> #o.shape.color=(1,1,1)
>
> #### Do you want friction on loading platens? (Rk: not much effect as the
> interparticle friction is set to the minimum value-> 10 here for the
> particles)
> O.bodies[wallIds[2]].mat.frictionAngle=radians(30)
> O.bodies[wallIds[3]].mat.frictionAngle=radians(30)
>
> triax.stressMask=7
> triax.goal1=confinement
> triax.goal2=confinement
> triax.goal3=confinement
> triax.max_vel=1
>
> for i in range(0,int(-1e-6*stress_max)-1):
> if ( abs(triax.goal2) < abs(stress_max) ):
> O.run(200,True)
> triax.goal2+=delta_stress
>
> triax.stressMask=5
> triax.goal1=confinement
> triax.goal2=strainRate
> triax.goal3=confinement
> triax.max_vel=1
>
>
> #triax.goal2=triax.goal2+strainRate
>
>
> 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
>

Hello Jan,

I use the formula of plastic dissipation used in sources.
In fact, i integrate the formula in the recorder like this :

tensCks=shearCks=cks=cks0=0
e10=e20=e30=0

def recorder():
    E=0
    Ub=0
    Uc=0
    Wp=0
    Wpla=0

    global tensCks, shearCks, e10,e20,e30
    tensCks=0
    shearCks=0
    for o in O.bodies:
 tensCks+=o.state.tensBreak
 shearCks+=o.state.shearBreak

    for i in O.interactions:
       #Ft=i.phys.shearForce
       #du=i.geom.shearInc
       #kt=i.phys.ks
       #Fn=i.phys.normalForce
       #Ftt=Ft-kt*du
       #maxFs=(Fn.norm()*tan(math.pi*FRICT/180)

      if((i.phys.shearForce-i.phys.ks*i.geom.shearInc).norm() > i.phys.normalForce.norm()*tan(math.pi*FRICT/180) ):
                Ftt=i.phys.shearForce-i.phys.ks*i.geom.shearInc
   ratio = i.phys.normalForce.norm()*tan(math.pi*FRICT/180) / Ftt.norm()
      trialForce=Ftt
      Ftt *= ratio
      e=((1/i.phys.ks)*(trialForce-Ftt)).dot(Ftt)
          Wpla=Wpla+e
      if not i.isReal : continue
      if isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Sphere):
        E+=0.5*(i.phys.normalForce.squaredNorm()/i.phys.kn + i.phys.shearForce.squaredNorm()/i.phys.ks)

    yade.plot.addData( t=O.time
   ,i=O.iter
   ,e1=triax.strain[0]-e10
   ,e2=triax.strain[1]-e20
   ,e3=triax.strain[2]-e30
   ,s1=triax.stress(triax.wall_right_id)[0]
   ,s2=triax.stress(triax.wall_top_id)[1]
   ,s3=triax.stress(triax.wall_front_id)[2]
   ,tc=0.5*tensCks,sc=0.5*shearCks,unbF=utils.unbalancedForce()
   ,E=E #Energie elastique
   ,Wp=triax.externalWork #Travail des contraintes appliquees aux limites de echantillon
   ,Wk=utils.kineticEnergy() #Energie cinetique
   ,Wd=2*utils.kineticEnergy()*DAMP*0.1*utils.PWaveTimeStep() #Energie dissipee par amortissement
       ,Wpla=Wpla #energie plastique
    )
    plot.saveDataTxt(OUT)

Hello,

Finally, do you compute plastic dissipation using traceEnergy=True of Law2_ScGeom_FrictPhys_CundallStrack ? (1)
Or do you compute it yourself in recorders ? (2)

If (1), it is not possible with "JCFpm" simulations (as Jan said)
If (2), get rid once for all of Law2_ScGeom_FrictPhys_CundallStrack which is useless, and problematic in this case (as Jan said)
;-)

Jerome

Hello Jerome,

I compute the plastic dissipation by myself in the recorders.
I think that it is possible to calculate the plastic dissipation in "JCFpm"simulations if we insert the formula
of plastic dissipation in "JointedCohesiveFrictionalPM.cpp/hpp" (using the Mohr-Coulomb criterion) and the results will be obtained in the file "cracks_____.txt"

what do you think ?

Jabrane

Yes, it is possible to modify the code in JointedCohesiveFrictionalPM.cpp/hpp to output some energy data in cracks_____.txt, and it would probably be useful.

BUT in case of JCFpm simulations and broken interactions, dissipated energy during plastic sliding is not an adequate variable.... The file cracks___.txt is specifically about broken interactions that lose their cohesive feature. At that failure point, there is indeed some change in energy, but in a brittle sense : the stored elastic energy suddenly decreases (for shear failure, in case there is an overlap) or vanishes (for tensile failure or shear failure without overlap)..

In case of JCFpm, sliding plastic dissipation makes sense for non cohesive interactions that reached the frictional threshold only. Such interactions are completely not related with the file cracks___.txt

Jerome

PS : we are diverging from the initial question, here.. ;-)