yade-users team mailing list archive

[Yor1 <question276604@xxxxxxxxxxxxxxxxxxxxx>]

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.