Yade

how to define friction at the boundaries in triaxial test

Asked by stephen

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=15,dynamic=1))
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"]

thanks in advance.

best regards
stephen

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Christian Jakob (jakob-ifgt) said : 		#1

Hi stephen,

Welcome to yade community.

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:

O.bodies.append(sphere([0,0,0],radius=R, material = SphereMat, fixed=True))
O.bodies.append(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=15, material = BoundaryMat, dynamic=1))

Hope it helps,

Christian

Revision history for this message
stephen (ifeanyi2295) said : 		#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.
thanks in advanced.

#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.write(str(0)+"\n") #adding zero deviatoric stress
devi.close()

strain=open("strain.txt","w")
#strain.write(str(0)+"\n") #adding zero strain
strain.close()

sigma11=open("sigma11.txt","w") # Axial stress file
#strain.write(str(0)+"\n") #adding zero strain
sigma11.close()

shear11=open("shear11.txt","w")# shear stress file. shear stress = ((sigma1-sigma2)/2)
#strain.write(str(0)+"\n") #adding zero strain
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")
#strain.write(str(0)+"\n") #adding zero strain
hop.close()

from yade import pack, plot

young=1e6 #young modulus of the spheres [Pa] (1e6 kPa=1e9 Pa)
pr=0.25 #poisson's ratios of the spheres
fa=28.6 #friction angle of the spheres [rad] 0.2 radian=11.46 degree, 0.5 radian= 28.6 degree, 0.4 rad=22.9 degreee, 0.8 rad=45.8 degree
fb=28.6 #friction at the boundaries in degree
number=1500

##create materials for spheres and plates

idSand=O.materials.append(FrictMat(young=young,poisson=pr,frictionAngle=radians(fa),label='idSand')) #defining sphere properties
BoundaryMat = FrictMat(young=young,poisson=pr,frictionAngle=fb)

O.materials.append(BoundaryMat)

O.bodies.append(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=15, material = BoundaryMat, dynamic=1))
O.bodies.append(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=15, dynamic=1))

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

  if utils.unbalancedForce()<0.1: O.bodies.append(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=15,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)

   checker.command='load()'

print "loading is starting!"

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

def load():

  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:

     print "loading is stopped"
     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()'
     print "### Axial Loading is Starting ###"

#### 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='addPlotData()',iterPeriod=200)]

  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")
  uzunluk=leng.readlines()

  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")
  #numbers=length.readlines()

  #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.pause() #this stops yade!

O.run()

Revision history for this message
Christian Jakob (jakob-ifgt) said : 		#3

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

##### 1. problem:

idSand=O.materials.append(FrictMat(young=young,poisson=pr,frictionAngle=radians(fa),label='idSand'))
#frictionAngle=fa !!!

##### 2. problem:

O.bodies.append(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=15, material = BoundaryMat, dynamic=1))
#this is ok

O.bodies.append(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=15, dynamic=1))
#missing material = BoundaryMat !!!

O.bodies.append(utils.geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=16,dynamic=0))
#missing material = BoundaryMat !!!

O.bodies.append(utils.geom.facetBox((0.035,0.035,0.105),(0.035,0.035,0.035),wallMask=15))
#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])
#instead of sp.toSimulation()

regards,

christian

Revision history for this message
stephen (ifeanyi2295) said : 		#4

hi christian,
thanks for your observation,
i have made the necessary corrections and get the best output that correspond to exertly what i was expecting.
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.write(str(0)+"\n") #adding zero deviatoric stress
devi.close()

strain=open("strain.txt","w")
#strain.write(str(0)+"\n") #adding zero strain
strain.close()

sigma11=open("sigma11.txt","w") # Axial stress file
#strain.write(str(0)+"\n") #adding zero strain
sigma11.close()

shear11=open("shear11.txt","w")# shear stress file. shear stress = ((sigma1-sigma2)/2)
#strain.write(str(0)+"\n") #adding zero strain
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")
#strain.write(str(0)+"\n") #adding zero strain
hop.close()

from yade import pack, plot

young_modulus=1e6 #young modulus of the spheres [Pa] (1e6 kPa=1e9 Pa)
poison_ratio=0.25 #poisson's ratios of the spheres
fa=0.5 #friction angle of the spheres [rad] 0.2 radian=11.46 degree, 0.5 radian= 28.6 degree, 0.4 rad=22.9 degreee, 0.8 rad=45.8 degree
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(geom.facetBox((0.035,0.035,0.035),(0.035,0.035,0.035),wallMask=15, material = BoundaryMat, dynamic=1))
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)

   checker.command='load()'

print "loading is starting!"

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

def load():

  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:

     print "loading is stopped"
     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()'
     print "### Axial Loading is Starting ###"

#### 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='addPlotData()',iterPeriod=200)]

  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")
  uzunluk=leng.readlines()

  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")
  #numbers=length.readlines()

  #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.pause() #this stops yade!

O.run()

Revision history for this message
Jérôme Duriez (jduriez) said : 		#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 ??

Revision history for this message
Christian Jakob (jakob-ifgt) said : 		#6

Hi,

> 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?.

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

Can you help with this problem?

Provide an answer of your own, or ask stephen for more information if necessary.

To post a message you must log in.