yade-users team mailing list archive

[Yor1 <question289662@xxxxxxxxxxxxxxxxxxxxx>]

New question #289662 on Yade:
https://answers.launchpad.net/yade/+question/289662

Hello,

I try to simulate triaxial compression test on fractured sample in order to determine the energy balance.
To calculate the elastic energy i use this formula:
Elastic energy=Fn^2/kn + Fs^2/ks with Fn, Fs normal and shear forces and kn, ks normal and shear stiffnesses.

When i run the code i have this error message:
ZeroDivisionError: float division by zero.

Best regards.
Jabrane.

This is the code:

from yade import ymport, utils , plot
import math

PACKING='X1Y2Z1_20k'
OUT=PACKING+'_compressionTest'


DAMP=0.4 # numerical damping
saveData=100 # data record interval
iterMax=2000000 # maximum number of iteration (to be adjusted)
saveVTK=10000 # Vtk files record interval


confinement=-1e6
#uniaxial_stress=-1e6
#delta_stress=-1e6
#stress_max=-200e6
strainRate=-0.02


intR=1.455
DENS=4000
YOUNG=65e9 
FRICT=10
ALPHA=0.4
TENS=8e6 
COH=160e6


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))


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() 

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)

O.bodies.append(ymport.text(PACKING+'.spheres',scale=1.,shift=Vector3(0,0,0),material=sphereMat))

import gts

c=X/4
alpha=math.pi/4

ptA =  gts.Vertex( X/2. + c/2.*cos(alpha), Y/2. + c/2.*sin(alpha), 8./7.*Z)
ptB =  gts.Vertex( X/2. - c/2.*cos(alpha), Y/2. - c/2.*sin(alpha), 8./7.*Z)
ptApr = gts.Vertex(X/2. + c/2.*cos(alpha), Y/2. + c/2.*sin(alpha), -1./7.*Z)
ptBpr = gts.Vertex(X/2. - c/2.*cos(alpha), Y/2. - c/2.*sin(alpha), -1./7.*Z)
e1 = gts.Edge(ptA,ptB)
e2 = gts.Edge(ptA,ptApr)
e3 = gts.Edge(ptApr,ptB)
f1 = gts.Face(e1,e2,e3)

e4 = gts.Edge(ptB,ptBpr)
e5 = gts.Edge(ptBpr,ptApr)
f2 = gts.Face(e4,e5,e3)

s1 = gts.Surface()
s1.add(f1)
s1.add(f2)
facet = gtsSurface2Facets(s1,wire = False,material=wallMat)
O.bodies.append(facet)

execfile('identifBis.py')

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(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')
]


tensCks=shearCks=cks=cks0=0

def recorder():
   
    E=0

    global tensCks, shearCks, e10,e20,e30
    tensCks=0
    shearCks=0
    e10=0
    e20=0
    e30=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)) or (isinstance(O.bodies[i.id1].shape,Box) and isinstance(O.bodies[i.id2].shape,Sphere)) or (isinstance(O.bodies[i.id1].shape,Sphere) and isinstance(O.bodies[i.id2].shape,Box)):
       	        E+=0.5*(i.phys.normalForce.squaredNorm()/i.phys.kn + i.phys.shearForce.squaredNorm()/i.phys.ks)

  
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()
  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

triax.stressMask=5
triax.goal1=confinement
triax.goal2=strainRate
triax.goal3=confinement
triax.max_vel=1

O.run(iterMax)


and this is the code names identifBis.py:

interactionRadius=intR;
O.engines=[

	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),Bo1_Facet_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg'),Ig2_Facet_Sphere_ScGeom()],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True,label='interactionLaw')]
	),
	NewtonIntegrator(damping=1)

]

############################ timestep + opening yade windows
O.dt=0.001*utils.PWaveTimeStep()

# from yade import qt
# v=qt.Controller()
# v=qt.View()

############################ Identification spheres on joint
#### color set for particles on joint
jointcolor1=(0,1,0)
jointcolor2=(1,0,0)
jointcolor3=(0,0,1)
jointcolor4=(1,1,1)
jointcolor5=(0,0,0)

#### first step-> find spheres on facet
O.step();

for i in O.interactions:
    ##if not i.isReal : continue
    ### Rk: facet are only stored in id1 
    if isinstance(O.bodies[i.id1].shape,Facet) and isinstance(O.bodies[i.id2].shape,Sphere): 
	vertices=O.bodies[i.id1].shape.vertices
	normalRef=vertices[0].cross(vertices[1]) # defines the normal to the facet normalRef
	nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
	normalFacetSphere=i.geom.normal # geom.normal is oriented from id1 to id2 -> normalFacetSphere from facet (i.id1) to sphere (i.id2)

	if O.bodies[i.id2].state.onJoint==False : ## particles has not yet been identified as belonging to a joint plane
	    O.bodies[i.id2].state.onJoint=True
	    O.bodies[i.id2].state.joint=1
	    O.bodies[i.id2].shape.color=jointcolor1
	    if nRef.dot(normalFacetSphere)>=0 :
		O.bodies[i.id2].state.jointNormal1=nRef
	    elif nRef.dot(normalFacetSphere)<0 :
		O.bodies[i.id2].state.jointNormal1=-nRef
	elif O.bodies[i.id2].state.onJoint==True :  ## particles has already been identified as belonging to, at least, 1 facet
	    if O.bodies[i.id2].state.joint==1 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) : ## particles has already been identified as belonging to only 1 facet
		O.bodies[i.id2].state.joint=2
		O.bodies[i.id2].shape.color=jointcolor2
		if nRef.dot(normalFacetSphere)>=0 :
		    O.bodies[i.id2].state.jointNormal2=nRef
		elif nRef.dot(normalFacetSphere)<0 :
		    O.bodies[i.id2].state.jointNormal2=-nRef
	    elif O.bodies[i.id2].state.joint==2 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05): ## particles has already been identified as belonging to more than 1 facet
		O.bodies[i.id2].state.joint=3
		O.bodies[i.id2].shape.color=jointcolor3
		if nRef.dot(normalFacetSphere)>=0 :
		    O.bodies[i.id2].state.jointNormal3=nRef
		elif nRef.dot(normalFacetSphere)<0 :
		    O.bodies[i.id2].state.jointNormal3=-nRef
	    elif O.bodies[i.id2].state.joint==3 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal3.cross(nRef)).norm()>0.05):
		O.bodies[i.id2].state.joint=4
		O.bodies[i.id2].shape.color=jointcolor5

#### second step -> find spheres interacting with spheres on facet (could be executed in the same timestep as step 1?)
for j in O.interactions:
    #if not i.isReal : continue
    ## Rk: facet are only stored in id1 
    if isinstance(O.bodies[j.id1].shape,Facet) and isinstance(O.bodies[j.id2].shape,Sphere): 
	vertices=O.bodies[j.id1].shape.vertices
	normalRef=vertices[0].cross(vertices[1]) # defines the normal to the facet normalRef
	nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
	if ((O.bodies[j.id2].state.jointNormal1.cross(nRef)).norm()<0.05) :
	    jointNormalRef=O.bodies[j.id2].state.jointNormal1
	elif ((O.bodies[j.id2].state.jointNormal2.cross(nRef)).norm()<0.05) :
	    jointNormalRef=O.bodies[j.id2].state.jointNormal2
	elif ((O.bodies[j.id2].state.jointNormal3.cross(nRef)).norm()<0.05) :
	    jointNormalRef=O.bodies[j.id2].state.jointNormal3
	else : continue
	facetCenter=O.bodies[j.id1].state.pos
	#### seek for each sphere interacting with the identified sphere j.id2
	for n in O.interactions.withBody(j.id2) :
            if isinstance(O.bodies[n.id1].shape,Sphere) and isinstance(O.bodies[n.id2].shape,Sphere):
                if j.id2==n.id1: # the sphere that was detected on facet (that is, j.id2) is id1 of interaction n
                    sphOnF=n.id1
                    othSph=n.id2
                elif j.id2==n.id2: # here, this sphere that was detected on facet (that is, j.id2) is id2 of interaction n
                    sphOnF=n.id2
                    othSph=n.id1
		facetSphereDir=(O.bodies[othSph].state.pos-facetCenter)
		if O.bodies[othSph].state.onJoint==True :
		    if O.bodies[othSph].state.joint==3 and ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].state.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].state.jointNormal3.cross(jointNormalRef)).norm()>0.05):
			O.bodies[othSph].state.joint=4
			O.bodies[othSph].shape.color=jointcolor5
		    elif O.bodies[othSph].state.joint==2 and ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].state.jointNormal2.cross(jointNormalRef)).norm()>0.05):
			O.bodies[othSph].state.joint=3
			if facetSphereDir.dot(jointNormalRef)>=0:
			    O.bodies[othSph].state.jointNormal3=jointNormalRef
			elif facetSphereDir.dot(jointNormalRef)<0:
			    O.bodies[othSph].state.jointNormal3=-jointNormalRef 
		    elif O.bodies[othSph].state.joint==1 and ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
			O.bodies[othSph].state.joint=2
			if facetSphereDir.dot(jointNormalRef)>=0:
			    O.bodies[othSph].state.jointNormal2=jointNormalRef
			elif facetSphereDir.dot(jointNormalRef)<0:
			    O.bodies[othSph].state.jointNormal2=-jointNormalRef
		elif  O.bodies[othSph].state.onJoint==False :
		    O.bodies[othSph].state.onJoint=True
		    O.bodies[othSph].state.joint=1
		    O.bodies[othSph].shape.color=jointcolor4
		    if facetSphereDir.dot(jointNormalRef)>=0:
			O.bodies[othSph].state.jointNormal1=jointNormalRef
		    elif facetSphereDir.dot(jointNormalRef)<0:
			O.bodies[othSph].state.jointNormal1=-jointNormalRef


O.resetTime()
O.interactions.clear()
print '\n IdentificationSpheresOnJoint executed ! Spheres onJoint (and so on...) detected, facets deleted\n\n'




		

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