[Yade-users] [Question #289705]: Minimal cohesive simulation

[Hotei]

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

Hello. Can I modify the following script (modified from 
https://answers.dogfood.paddev.net/yade/+question/254316), in order to the 
upper sphere drags the lower one, by means of a cohesive contact? Thanks a lot 
in advance

from yade import *
import sys, time, random, os, gts, math
from yade import ymport
from yade import pack
from yade import qt
from yade import bodiesHandling

O.materials.append(CohFrictMat(
young=3e8,
poisson=0.3,
frictionAngle=radians(30),
density=2600,
isCohesive=True,
alphaKr=0.031,
alphaKtw=0.031,
momentRotationLaw=False,
etaRoll=5.0,
label='granular_material',
normalCohesion = 10.**9.,
fragile = False))

s = utils.sphere(center=[0,0,0],radius=1)
p = utils.sphere(center=[0,0,2],radius=1)

vel = 0.0001

O.bodies.append(s)
O.bodies.append(p)

O.bodies[1].state.vel=(0,0,vel)

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D()],
[Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
[Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]),
NewtonIntegrator(damping=0.2)
]

O.dt = utils.PWaveTimeStep()
O.step()

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

Re: [Yade-users] [Question #289662]: ZeroDivisionError: float division by zero

[Jan Stránský]

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

Status: Open => Answered

Jan Stránský proposed the following answer:
Hi Jabrane,

please provide also the file X1Y2Z1_20k.spheres such that we can test it
ourselves.

Also the error says exactly on which line the error occurred, could you
please send the line alone in next message (such that we do not need to go
through the whole file)? Or just verify that the error is really at the
energy computation line and not related to some other divisions in the
scripts.

Thanks
Jan


2016-03-31 11:43 GMT+02:00 Yor1 :

> 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=200 # maximum number of iteration (to be adjusted)
> saveVTK=1 # 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
> 

[Yade-users] [Question #289663]: Total elastic shear energy in Hertz-Mindlin model

[Jiangtao Wei]

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

Dear all,
My question is how to calculate the total shear elastic energy when using 
Hertz-Mindlin model. In the code, MindlinPhys_Mindlin.shearEnergy only gives 
the incremental shear elastic energy. 

Thanks you.


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

[Yade-users] [Question #289662]: ZeroDivisionError: float division by zero

[Yor1]

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=200 # maximum number of iteration (to be adjusted)
saveVTK=1 # 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