[Yade-users] [Question #289705]: Minimal cohesive simulation
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
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
> isinstance(O.bodies[i.id2].shape,Box)):
> E+=0.5*(i.phys.
[Yade-users] [Question #289663]: Total elastic shear energy in Hertz-Mindlin model
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
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
triax.goal3=confine
Re: [Yade-users] [Question #288670]: What's XSMP error?
Question #288670 on Yade changed: https://answers.launchpad.net/yade/+question/288670 Status: Open => Expired Launchpad Janitor expired the question: This question was expired because it remained in the 'Open' state without activity for the last 15 days. -- 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
