# how to define friction at the boundaries in triaxial test

hi ,
please how can one define the friction at the top and bottom of the specimen in a Triaxial test simulation of a granular material.

the center ,extension and the wall mask of my model are as below.

O.bodies.append(utils.geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=16,dynamic=0))#bottom of the sample, which is not allowed to move. because of that dynamic=0 [O means "False"]

best regards
stephen

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2013-09-05
2013-09-05
 Christian Jakob (jakob-ifgt) said on 2013-09-04: #1

Hi stephen,

In yade you have to define materials (e.g. one material for particles and one for boundaries):

### properties:
density = 2650
shear_modulus = 2e5
poisson_ratio = 0.15
young_modulus = 2*shear_modulus*(1+poisson_ratio)

### define materials:
SphereMat = FrictMat(young=young_modulus,poisson=poisson_ratio,density=density,frictionAngle=0.5)
BoundaryMat = FrictMat(young=young_modulus,poisson=poisson_ratio,frictionAngle=1)

O.materials.append(SphereMat)
O.materials.append(BoundaryMat)

Then you can use the materials in the generation process:

Hope it helps,

Christian

 stephen (ifeanyi2295) said on 2013-09-04: #2

hi Jakob,
thanks for your prompt response and great contribution.
however, i made some changes and try to run my script again, but the result i was geting is not what i was expected. here is the detailed copy of the script for your suggestions and corrections.

#file=open("data1.txt","w")
#file.close()

#stra=open("strain.txt","w")
#stra.write(str(0)+"\n") #because we need to add the zero strain. #look at later, may be we should also add zero force???
#stra.close()

devi=open("deviatoric.txt","w") #deviatoric stresses file
devi.close()

strain=open("strain.txt","w")
strain.close()

sigma11=open("sigma11.txt","w") # Axial stress file
sigma11.close()

shear11=open("shear11.txt","w")# shear stress file. shear stress = ((sigma1-sigma2)/2)
shear11.close()

epsellion=("epsellion.txt","w")
#epsellion.close()

file=open("sigma1.txt","w")
file1=open("sigma2.txt","w")
file.close()
file1.close()

hop=open("leng.txt","w")
hop.close()

young=1e6 #young modulus of the spheres [Pa] (1e6 kPa=1e9 Pa)
pr=0.25 #poisson's ratios of the spheres
fb=28.6 #friction at the boundaries in degree
number=1500

##create materials for spheres and plates

BoundaryMat = FrictMat(young=young,poisson=pr,frictionAngle=fb)

O.materials.append(BoundaryMat)

O.bodies.append(utils.geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=16,dynamic=0))#bottom of the sample, which is not allowed to move. because of that dynamic=0 [O means "False"]

O.bodies.append(utils.geom.facetBox((0.035,0.035,0.105),(0.035,0.035,0.035),wallMask=15))#!!!! in order to control the particles!!! after the end of the stabilization, we will remove it

psdSizes,psdCumm=[0.00525,0.00603,0.0105,0.01206,0.02099],[0,0.1,0.5,0.6,1.0] #d50=10.5mm, cu=2

mincorner,maxcorner=Vector3(0.001,0.001,0.001),Vector3(0.07,0.07,0.14) #min and max coordinates of the wall

sp=pack.SpherePack() #burda pack modulunun Spherepack modulunu cagiriyo, bi altinda da o cagirdinin,piramit gibi lan iste!
sp.makeCloud(mincorner,maxcorner,psdSizes=psdSizes,num=number,psdCumm=psdCumm,distributeMass=1)
#sp.makeCloud(mincorner,maxcorner,rMean=0.005,rRelFuzz=0.3)
Gl1_Sphere.stripes=True #showing the spheres in two colors

sp.toSimulation()#color=(0,0,1))

unb=utils.unbalancedForce()

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
InteractionLoop(

[Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_L3Geom_FrictPhys_ElPerfPl()]
),
GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(damping=0.7), #after the first phase we are going to change it to 0.5
PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),

]

O.dt=0.3*utils.PWaveTimeStep()

#### D E P O S I T I O N ####

print '### Deposion is Started ###'

def checkUnbalanced():

print 'iteration number:',O.iter,'unbalanced force:', unbalancedForce()
if O.iter<25000: return
# the rest will be run only if unbalanced is < .1 (stabilized packing)
if utils.unbalancedForce()>0.1: return

checker.command='erasingwalls()'
print "Depostion is Completed and the Top Part of the Membran is Added!"

def erasingwalls():
print "### Additional Walls are Erased! ###"
O.bodies.erase(10)
O.bodies.erase(11)
O.bodies.erase(12)
O.bodies.erase(13)
O.bodies.erase(14)
O.bodies.erase(15)
O.bodies.erase(16)
O.bodies.erase(17)

aim=100000 #PASCAL!!

vel=0.03 #applyied velocity
#0.003 calisii

#### C O N F I N I N G P R E S S U R E ####

O.materials[idSand].young=1e9
O.engines=O.engines+[PyRunner(command='writingDatas()',iterPeriod=100000)]
print ('CONFINING PRESSURE')
print 'stress at the top :',(abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))

print 'stress at the front:',(abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the back:',(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the right:',(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the left:',(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

if (abs((O.forces.f(O.bodies[number+18].id)[2])+(O.forces.f(O.bodies[number+19].id)[2])))/((abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))))<aim:

O.bodies[(number+18)].state.vel=(0,0,(-vel)) #top
O.bodies[(number+19)].state.vel=(0,0,(-vel))

if (abs((O.forces.f(O.bodies[number+18].id)[2])+(O.forces.f(O.bodies[number+19].id)[2])))/((abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))))>aim:
O.bodies[(number+18)].state.vel=(0,0,vel/2) #top
O.bodies[(number+19)].state.vel=(0,0,vel/2)

if (abs((O.forces.f(O.bodies[2].id)[0])+(O.forces.f(O.bodies[3].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))<aim:
O.bodies[2].state.vel=(-vel,0,0) #front
O.bodies[3].state.vel=(-vel,0,0)
if (abs((O.forces.f(O.bodies[2].id)[0])+(O.forces.f(O.bodies[3].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim:
O.bodies[2].state.vel=(vel/2,0,0) #front
O.bodies[3].state.vel=(vel/2,0,0)

if (abs((O.forces.f(O.bodies[0].id)[0])+(O.forces.f(O.bodies[1].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))<aim:
O.bodies[0].state.vel=(vel,0,0) #behind
O.bodies[1].state.vel=(vel,0,0)
if (abs((O.forces.f(O.bodies[0].id)[0])+(O.forces.f(O.bodies[1].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim:
O.bodies[0].state.vel=(-vel/2,0,0) #behind
O.bodies[1].state.vel=(-vel/2,0,0)

if (abs((O.forces.f(O.bodies[6].id)[1])+(O.forces.f(O.bodies[7].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[6].state.vel=(0,-vel,0) #right
O.bodies[7].state.vel=(0,-vel,0)
if (abs((O.forces.f(O.bodies[6].id)[1])+(O.forces.f(O.bodies[7].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[6].state.vel=(0,(vel/2),0) #right
O.bodies[7].state.vel=(0,(vel/2),0)

if (abs((O.forces.f(O.bodies[4].id)[1])+(O.forces.f(O.bodies[5].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[4].state.vel=(0,vel,0) #left
O.bodies[5].state.vel=(0,vel,0)
if (abs((O.forces.f(O.bodies[4].id)[1])+(O.forces.f(O.bodies[5].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[4].state.vel=(0,-vel/2,0) #left
O.bodies[5].state.vel=(0,-vel/2,0)

if (abs((O.forces.f(O.bodies[number+18].id)[2])+(O.forces.f(O.bodies[number+19].id)[2])))/((abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))))>aim and (abs((O.forces.f(O.bodies[2].id)[0])+(O.forces.f(O.bodies[3].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim and (abs((O.forces.f(O.bodies[0].id)[0])+(O.forces.f(O.bodies[1].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim and (abs((O.forces.f(O.bodies[6].id)[1])+(O.forces.f(O.bodies[7].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim and (abs((O.forces.f(O.bodies[4].id)[1])+(O.forces.f(O.bodies[5].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:

uzun=(abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))
print uzun
hip=open("leng.txt","a")
hip.write(str(uzun)+"\n")
hip.close()
checker.command='axial()'

#### A X I A L L O A D I N G ####

def axial():
O.materials[idSand].young=1e9

O.engines=O.engines+[PyRunner(command='writingLength()',iterPeriod=150000)]

print 'stress at the top :',(abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))),'iteration number:',O.iter

print 'stress at the front:',(abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the back:',(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the right:',(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the left:',(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

if (abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))<4*aim:
O.bodies[number+18].state.vel=(0,0,-2*(vel)) #top
O.bodies[number+19].state.vel=(0,0,-2*(vel))

if (abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[2].state.vel=((-vel/4),0,0) #front
O.bodies[3].state.vel=((-vel/4),0,0)

if (abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[2].state.vel=((vel/4),0,0) #front
O.bodies[3].state.vel=((vel/4),0,0)
O.bodies[number+18].state.vel=(0,0,-(vel)) #top
O.bodies[number+19].state.vel=(0,0,-(vel))

if (abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[0].state.vel=((vel/4),0,0) #behind
O.bodies[1].state.vel=((vel/4),0,0)

if (abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[0].state.vel=((-vel/4),0,0) #behind
O.bodies[1].state.vel=((-vel/4),0,0)
O.bodies[number+18].state.vel=(0,0,-(vel)) #top
O.bodies[number+19].state.vel=(0,0,-(vel))

if (abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[6].state.vel=(0,(-vel/4),0) #right
O.bodies[7].state.vel=(0,(-vel/4),0)

if (abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[6].state.vel=(0,(vel/4),0) #right
O.bodies[7].state.vel=(0,(vel/4),0)
O.bodies[number+18].state.vel=(0,0,-(vel)) #top
O.bodies[number+19].state.vel=(0,0,-(vel))

if (abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[4].state.vel=(0,(vel/4),0) #left
O.bodies[5].state.vel=(0,(vel/4),0)

if (abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[4].state.vel=(0,(-vel/4),0) #left
O.bodies[5].state.vel=(0,(-vel/4),0)
O.bodies[number+18].state.vel=(0,0,-(vel)) #top
O.bodies[number+19].state.vel=(0,0,-(vel))

if (abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))>4*aim:

print "### Triaxial Test is Completed ###"
print "### ENDE GUT ALLES GUT ###"
print "### Thank you JAMAL ###"

checker.command='strain()'

def writingDatas():

x=(abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))

y=(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

file=open("sigma1.txt","a")
file1=open("sigma2.txt","a")

file.write(str(x)+"\n")
file1.write(str(y)+"\n")
#+"\n")
file.close()
file1.close()

for b in O.bodies: b.shape.color=utils.scalarOnColorScale(b.state.rot().norm(),0,pi/2.)
# rot() gives rotation vector between reference and current position

def writingLength():

leng=open("leng.txt","r")

first=(float(uzunluk[0]))
leng.close()

if (abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))<4*aim:

s=((0.07)-abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))/(0.07)

E=((0.07)-abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1])))/(0.07)

shr=((((abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))-(((abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))/4)))/1000)/2
sig=((((abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))+(((abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))/4)))/1000)

y=(((abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))-(((abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))/4)))/1000

devi=open("deviatoric.txt","a")
devi.write(str(y)+"\n")
devi.close()

strain=open("strain.txt","a")
strain.write(str(s)+"\n")
strain.close()

epsellion=open("epsellion.txt","a")
epsellion.write(str(E)+"\n")
epsellion.close()

sigma11=open("sigma11.txt","a")
sigma11.write(str(sig)+"\n")
sigma11.close()

shear11=open("shear11.txt","a")
shear11.write(str(shr)+"\n")
shear11.close()

print ("strain:",s," ","epsellion:",E," ","deviatoric:",y)

def strain():

# length=open("length.txt","r")

#for i in range(0,(len(numbers)-1)):
#"stra" defined as the file of strains calculated for each iteration.the calculation starts with the end of axial loading.
# stra=open("strain.txt","a")
#strain=(abs(float(numbers[i]))-(float(numbers[i+1])))/(float(numbers[i]))
#stra.write(str(strain)+"\n")
#stra.close()

#print ("### calculation of the strain completed ###")
print ("_____________________________________________________________________")

O.run()

 Christian Jakob (jakob-ifgt) said on 2013-09-04: #3

in your script i see some inconsistencies, that can lead to unexpected results:

##### 1. problem:

#frictionAngle=fa !!!

##### 2. problem:

#this is ok

#missing material = BoundaryMat !!!

#missing material = BoundaryMat !!!

#missing material = BoundaryMat

##### 3.problem:

sp=pack.SpherePack()
sp.makeCloud(mincorner,maxcorner,psdSizes=psdSizes,num=number,psdCumm=psdCumm,distributeMass=1)

sp.toSimulation()

you did not set material for spheres. but i do not know how to do this with sp.toSimulation()
i prefer to use:

O.bodies.append([sphere(c,r,material=idSand) for c,r in sp])

regards,

christian

 stephen (ifeanyi2295) said on 2013-09-05: #4

hi christian,
how ever, i have some additional development during the simulation. an error named " FATAL /build/buildd/yade-0.70.0/core/ThreadRunner.cpp:31 run: Exception occured:
Body #37 has velocity==NaN! "
occured and this automatically terminates the simulation.

do you know anything that could be the cause of this problem and how to fix it?.
i have checked the similar post here but the problem still persist.

about setting the material for sphere, i used this command to do that

O.materials[idSand].young=1e6.

thanks for your great assistance....am greatful

best regards
stephen

you may also wish to go through the recent script shown below

#file=open("data1.txt","w")
#file.close()

#stra=open("strain.txt","w")
#stra.write(str(0)+"\n") #because we need to add the zero strain. #look at later, may be we should also add zero force???
#stra.close()

devi=open("deviatoric.txt","w") #deviatoric stresses file
devi.close()

strain=open("strain.txt","w")
strain.close()

sigma11=open("sigma11.txt","w") # Axial stress file
sigma11.close()

shear11=open("shear11.txt","w")# shear stress file. shear stress = ((sigma1-sigma2)/2)
shear11.close()

epsellion=("epsellion.txt","w")
#epsellion.close()

file=open("sigma1.txt","w")
file1=open("sigma2.txt","w")
file.close()
file1.close()

hop=open("leng.txt","w")
hop.close()

young_modulus=1e6 #young modulus of the spheres [Pa] (1e6 kPa=1e9 Pa)
poison_ratio=0.25 #poisson's ratios of the spheres
fb=0.5 #friction at the boundaries in
number=1400

##create materials for spheres and plates

idSand=O.materials.append(FrictMat(young=young_modulus,poisson=poison_ratio,frictionAngle=fa,label='idSand')) #defining sphere properties

#SphereMat = FrictMat(young=young_modulus,poisson=poisson_ratio,density=density,frictionAngle=0.5)

BoundaryMat = FrictMat(young=young_modulus,poisson=poison_ratio,frictionAngle=fb)

O.materials.append(BoundaryMat)

O.bodies.append(utils.geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=16, material = BoundaryMat,dynamic=0))#bottom of the sample, which is not allowed to move. because of that dynamic=0 [O means "False"]

O.bodies.append(utils.geom.facetBox((0.035,0.035,0.105),(0.035,0.035,0.035),wallMask=15, material = BoundaryMat))#!!!! in order to control the particles!!! after the end of the stabilization, we will remove it

psdSizes,psdCumm=[0.00525,0.00603,0.0105,0.01206,0.02099],[0,0.1,0.5,0.6,1.0] #d50=10.5mm, cu=2

mincorner,maxcorner=Vector3(0.001,0.001,0.001),Vector3(0.07,0.07,0.14) #min and max coordinates of the wall

sp=pack.SpherePack() #burda pack modulunun Spherepack modulunu cagiriyo, bi altinda da o cagirdinin,piramit gibi lan iste!
sp.makeCloud(mincorner,maxcorner,psdSizes=psdSizes,num=number,psdCumm=psdCumm,distributeMass=1)
#sp.makeCloud(mincorner,maxcorner,rMean=0.005,rRelFuzz=0.3)
Gl1_Sphere.stripes=True #showing the spheres in two colors

sp.toSimulation()#color=(0,0,1))

unb=utils.unbalancedForce()

O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
InteractionLoop(

[Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_L3Geom_FrictPhys_ElPerfPl()]
),
GravityEngine(gravity=(0,0,-9.81)),
NewtonIntegrator(damping=0.8), #after the first phase we are going to change it to 0.5
PyRunner(command='checkUnbalanced()',realPeriod=2,label='checker'),

]

O.dt=0.3*utils.PWaveTimeStep()

#### D E P O S I T I O N ####

print '### Deposion is Started ###'

def checkUnbalanced():

print 'iteration number:',O.iter,'unbalanced force:', unbalancedForce()
if O.iter<150000: return
# the rest will be run only if unbalanced is < .1 (stabilized packing)
if utils.unbalancedForce()>0.1: return

if utils.unbalancedForce()<0.1: O.bodies.append(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=32,material = BoundaryMat, dynamic=1))

checker.command='erasingwalls()'
print "Depostion is Completed and the Top Part of the Membran is Added!"

def erasingwalls():
print "### Additional Walls are Erased! ###"
O.bodies.erase(10)
O.bodies.erase(11)
O.bodies.erase(12)
O.bodies.erase(13)
O.bodies.erase(14)
O.bodies.erase(15)
O.bodies.erase(16)
O.bodies.erase(17)

aim=100000 #PASCAL!!

vel=0.003 #applyied velocity
#0.003 calisii

#### C O N F I N I N G P R E S S U R E ####

O.materials[idSand].young=1e6
O.engines=O.engines+[PyRunner(command='writingDatas()',iterPeriod=100000)]
print ('CONFINING PRESSURE')
print 'stress at the top :',(abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))

print 'stress at the front:',(abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the back:',(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the right:',(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the left:',(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

if (abs((O.forces.f(O.bodies[number+18].id)[2])+(O.forces.f(O.bodies[number+19].id)[2])))/((abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))))<aim:

O.bodies[(number+18)].state.vel=(0,0,(-vel)) #top
O.bodies[(number+19)].state.vel=(0,0,(-vel))

if (abs((O.forces.f(O.bodies[number+18].id)[2])+(O.forces.f(O.bodies[number+19].id)[2])))/((abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))))>aim:
O.bodies[(number+18)].state.vel=(0,0,vel/2) #top
O.bodies[(number+19)].state.vel=(0,0,vel/2)

if (abs((O.forces.f(O.bodies[2].id)[0])+(O.forces.f(O.bodies[3].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))<aim:
O.bodies[2].state.vel=(-vel,0,0) #front
O.bodies[3].state.vel=(-vel,0,0)
if (abs((O.forces.f(O.bodies[2].id)[0])+(O.forces.f(O.bodies[3].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim:
O.bodies[2].state.vel=(vel/2,0,0) #front
O.bodies[3].state.vel=(vel/2,0,0)

if (abs((O.forces.f(O.bodies[0].id)[0])+(O.forces.f(O.bodies[1].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))<aim:
O.bodies[0].state.vel=(vel,0,0) #behind
O.bodies[1].state.vel=(vel,0,0)
if (abs((O.forces.f(O.bodies[0].id)[0])+(O.forces.f(O.bodies[1].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim:
O.bodies[0].state.vel=(-vel/2,0,0) #behind
O.bodies[1].state.vel=(-vel/2,0,0)

if (abs((O.forces.f(O.bodies[6].id)[1])+(O.forces.f(O.bodies[7].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[6].state.vel=(0,-vel,0) #right
O.bodies[7].state.vel=(0,-vel,0)
if (abs((O.forces.f(O.bodies[6].id)[1])+(O.forces.f(O.bodies[7].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[6].state.vel=(0,(vel/2),0) #right
O.bodies[7].state.vel=(0,(vel/2),0)

if (abs((O.forces.f(O.bodies[4].id)[1])+(O.forces.f(O.bodies[5].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[4].state.vel=(0,vel,0) #left
O.bodies[5].state.vel=(0,vel,0)
if (abs((O.forces.f(O.bodies[4].id)[1])+(O.forces.f(O.bodies[5].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[4].state.vel=(0,-vel/2,0) #left
O.bodies[5].state.vel=(0,-vel/2,0)

if (abs((O.forces.f(O.bodies[number+18].id)[2])+(O.forces.f(O.bodies[number+19].id)[2])))/((abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))))>aim and (abs((O.forces.f(O.bodies[2].id)[0])+(O.forces.f(O.bodies[3].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim and (abs((O.forces.f(O.bodies[0].id)[0])+(O.forces.f(O.bodies[1].id)[0])))/((abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))>aim and (abs((O.forces.f(O.bodies[6].id)[1])+(O.forces.f(O.bodies[7].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim and (abs((O.forces.f(O.bodies[4].id)[1])+(O.forces.f(O.bodies[5].id)[1])))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:

uzun=(abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))
print uzun
hip=open("leng.txt","a")
hip.write(str(uzun)+"\n")
hip.close()
checker.command='axial()'

#### A X I A L L O A D I N G ####

def axial():
O.materials[idSand].young=1e6

O.engines=O.engines+[PyRunner(command='writingLength()',iterPeriod=150000)]

print 'stress at the top :',(abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1]))),'iteration number:',O.iter

print 'stress at the front:',(abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the back:',(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the right:',(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

print 'stress at the left:',(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

if (abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))<4*aim:
O.bodies[number+18].state.vel=(0,0,-2*(vel)) #top
O.bodies[number+19].state.vel=(0,0,-2*(vel))

if (abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[2].state.vel=((-vel/4),0,0) #front
O.bodies[3].state.vel=((-vel/4),0,0)

if (abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[2].state.vel=((vel/4),0,0) #front
O.bodies[3].state.vel=((vel/4),0,0)
O.bodies[number+18].state.vel=(0,0,-(vel)) #top
O.bodies[number+19].state.vel=(0,0,-(vel))

if (abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[0].state.vel=((vel/4),0,0) #behind
O.bodies[1].state.vel=((vel/4),0,0)

if (abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[0].state.vel=((-vel/4),0,0) #behind
O.bodies[1].state.vel=((-vel/4),0,0)
O.bodies[number+18].state.vel=(0,0,-(vel)) #top
O.bodies[number+19].state.vel=(0,0,-(vel))

if (abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[6].state.vel=(0,(-vel/4),0) #right
O.bodies[7].state.vel=(0,(-vel/4),0)

if (abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[6].state.vel=(0,(vel/4),0) #right
O.bodies[7].state.vel=(0,(vel/4),0)
O.bodies[number+18].state.vel=(0,0,-(vel)) #top
O.bodies[number+19].state.vel=(0,0,-(vel))

if (abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))<aim:
O.bodies[4].state.vel=(0,(vel/4),0) #left
O.bodies[5].state.vel=(0,(vel/4),0)

if (abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))>aim:
O.bodies[4].state.vel=(0,(-vel/4),0) #left
O.bodies[5].state.vel=(0,(-vel/4),0)
O.bodies[number+18].state.vel=(0,0,-(vel)) #top
O.bodies[number+19].state.vel=(0,0,-(vel))

if (abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))>4*aim:

print "### Triaxial Test is Completed ###"
print "### ENDE GUT ALLES GUT ###"
print "### Thank you JAMAL ###"

checker.command='strain()'

def writingDatas():

x=(abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))

y=(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))

file=open("sigma1.txt","a")
file1=open("sigma2.txt","a")

file.write(str(x)+"\n")
file1.write(str(y)+"\n")
#+"\n")
file.close()
file1.close()

for b in O.bodies: b.shape.color=utils.scalarOnColorScale(b.state.rot().norm(),0,pi/2.)
# rot() gives rotation vector between reference and current position

def writingLength():

leng=open("leng.txt","r")

first=(float(uzunluk[0]))
leng.close()

if (abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))<4*aim:

s=((0.07)-abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))/(0.07)

E=((0.07)-abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1])))/(0.07)

shr=((((abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))-(((abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))/4)))/1000)/2
sig=((((abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))+(((abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))/4)))/1000)

y=(((abs(O.forces.f(O.bodies[number+18].id)[2])+abs(O.forces.f(O.bodies[number+19].id)[2]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[4].state.pos[1])-(O.bodies[7].state.pos[1])))-(((abs(O.forces.f(O.bodies[2].id)[0])+abs(O.forces.f(O.bodies[3].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[0].id)[0])+abs(O.forces.f(O.bodies[1].id)[0]))/(abs((O.bodies[5].state.pos[1])-(O.bodies[6].state.pos[1]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[6].id)[1])+abs(O.forces.f(O.bodies[7].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2])))+(abs(O.forces.f(O.bodies[4].id)[1])+abs(O.forces.f(O.bodies[5].id)[1]))/(abs((O.bodies[2].state.pos[0])-(O.bodies[1].state.pos[0]))*abs((O.bodies[8].state.pos[2])-(O.bodies[number+19].state.pos[2]))))/4)))/1000

devi=open("deviatoric.txt","a")
devi.write(str(y)+"\n")
devi.close()

strain=open("strain.txt","a")
strain.write(str(s)+"\n")
strain.close()

epsellion=open("epsellion.txt","a")
epsellion.write(str(E)+"\n")
epsellion.close()

sigma11=open("sigma11.txt","a")
sigma11.write(str(sig)+"\n")
sigma11.close()

shear11=open("shear11.txt","a")
shear11.write(str(shr)+"\n")
shear11.close()

print ("strain:",s," ","epsellion:",E," ","deviatoric:",y)

def strain():

# length=open("length.txt","r")

#for i in range(0,(len(numbers)-1)):
#"stra" defined as the file of strains calculated for each iteration.the calculation starts with the end of axial loading.
# stra=open("strain.txt","a")
#strain=(abs(float(numbers[i]))-(float(numbers[i+1])))/(float(numbers[i]))
#stra.write(str(strain)+"\n")
#stra.close()

#print ("### calculation of the strain completed ###")
print ("_____________________________________________________________________")

O.run()

 Jérôme Duriez (jduriez) said on 2013-09-05: #5

Hello Stephen,

- It seems to me that maybe, these "material" are not perfectly clear to you (otherwise, please excuse my answer). You can define in your Yade-simulation (this "O" of O. ...) a list of different "materials" defined by their type, and the values of the corresponding parameters

But this list is useless until you do not affect one material of the list to one "body" (=element) of your simulation. So when you type O.materials[idSand].young=1e6, you only change one parameter value of one material of the list.

Are you sure it is the material that you affected to your spheres ? Note that you still use sp.toSimulation() (cf Christian's remark)

- Second, I fear problems with this "idSand" that is both a label of a material, and the return value of O.bodies.append() function, see your line :
idSand=O.materials.append(..........,label='idSand'))

I'm not yet enough python/yade - aware to be sure there is a problem, but I'm not sure that your line
O.materials[idSand].young=1e6.
does what you want ??

 Christian Jakob (jakob-ifgt) said on 2013-09-05: #6

Hi,

> how ever, i have some additional development during the simulation. an error named
> Body #37 has velocity==NaN! "
> occured and this automatically terminates the simulation.
>
> do you know anything that could be the cause of this problem and how to fix it?.

1.) follow statements/remarks of jduriez and do something like:

idSand=O.materials.append(..........,label='sandLabel'))

Then you can do:

sandLabel.young=1e6

or

O.materials[idSand].young=1e6

which is doing the same ...

But changing the property of a material will not change material properties of existing bodies.
Do you want to change material properties during the run?

2.) If your simulation is dynamical, time step could be too high. Decrease it either directly with a fixed O.dt= ...? or decrease timestepSafetyCoefficient in GlobalStiffnessTimeStepper engine.

Additionally I want to say, that GravityEngine is deprecated. You can set it directly in the NewtonIntegrator

NewtonIntegrator(gravity=(0,0,-9.81),damping=0.7,label='integrator')

and remove this line:

GravityEngine(gravity=(0,0,-9.81)),

Next time please send a short example script, that shows the problem. It is quite time consuming investigating a 1000-line script including 950 irrelevant lines.

Regards,

Christian