PFV model is not working without any issue

Asked by Sam

Hiya All,

I wrote a code that is supposed to simulate a Triaxial_Undrained test of a saturated sphere packing within a cylindrical specimen. The The code run without any issue but it's not working properly. Can anyone help me to make this code work?
 The code is copied below:

# import corrosponding yade modulus

from yade import pack

from yade import geom

from yade import plot

from yade import ymport

from yade import qt

from yade.gridpfacet import *

import gts, os.path, locale

locale.setlocale(locale.LC_ALL, 'en_US.UTF-8')

import sys,os

from yade import plot,pack,export

sys.path.append(os.pardir)

########################################################################

# A. DEFINING VARIABLES, MATERIALS' PROPERTIES AND PACKING #

########################################################################

# A.a). define variables

key = 'Triaxial_Undrained' # file name to be saved

young=550e6 # normal contact stiffness

compFricDegree = 1.8 # initial contact friction during the confining phase

finalFricDegree = 43 # contact friction during the deviatoric loading

poisson = 0.3 # shear-to-normal stiffness ratio

isoStress = 110000 # confining stress

conStress = 100000 # confining stress for deviatoric loading stage

width = 1.4e-1 # sample width

height = 2.8e-1 # target sample height(after consolidation)

height_0 = 3.2e-1 # initial sample height

num_spheres=500 # number of spheres

R_p = 0.0084 # mean particle radius

rCoff = 10 # thickness of top and bot sphere cap (based on rParticle)

rParticle = 0.02e-1 # membrane grid seed size

alpha = 8

rate = 0.1 # loading rate (strain rate)

damp = 0.3 # damping coefficient

targetPorosity = 0.43 # target porosity

thresholdvalue = 0.05 # threshold unbalance force

final_rate = 0.1 # strain rate for deviator loading

thresholdstrain = 0.06 # threshold axial strain for terminate

enlargefactor = 1.00

# A.b). create materials for sand spheres and plates

Sand = O.materials.append(CohFrictMat(

young=young,poisson=poisson,frictionAngle=radians(compFricDegree),

alphaKr=0.25,alphaKtw=0,etaRoll=0.005,etaTwist=0,

normalCohesion=5e6,shearCohesion=5e6,

momentRotationLaw=True,density=2650,label='spheres'

))

# A.c). create membrane materials

207

GridMat = O.materials.append(CohFrictMat(

young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),

alphaKr=0,alphaKtw=0,etaRoll=0,etaTwist=0,

normalCohesion=1e9,shearCohesion=1e9,

momentRotationLaw=True,label='gridNodeMat'

))

pFacetMat = O.materials.append(FrictMat(

young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),label='pFacetMat'

))

# A.d). create TOP & BOT plate materials

frictMat = O.materials.append(FrictMat(

young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),label='frictMat'

))

###################################################################

# B. DEFINING GLOBAL ENGINES #

###################################################################

#**********************************************************************#

O.engines=[

ForceResetter(),

InsertionSortCollider([

Bo1_Sphere_Aabb(),

Bo1_PFacet_Aabb(),

Bo1_Facet_Aabb(),

Bo1_GridConnection_Aabb()

]),

InteractionLoop(

[

Ig2_Sphere_Sphere_ScGeom6D(),

Ig2_GridNode_GridNode_GridNodeGeom6D(),

Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),

Ig2_Sphere_PFacet_ScGridCoGeom(),

Ig2_Facet_Sphere_ScGeom6D()

],

[

Ip2_FrictMat_FrictMat_FrictPhys(),

Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")

],

[

Law2_ScGeom_FrictPhys_CundallStrack(),

Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm

=True,always_use_moment_law=False,label='cohesiveLaw'),

Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),

Law2_ScGridCoGeom_FrictPhys_CundallStrack(),

],label="iloop"

),

GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8),

NewtonIntegrator(gravity=(0,0,0),damping=0.1,label='newton')

]

#**********************************************************************#

O.engines[2].lawDispatcher.functors[1].always_use_moment_law=False

O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False

#######################################################################

# C. GENERATING PACKING #

#######################################################################

# C.a). generate random dense sphere pack

sp=pack.SpherePack()

pred = pack.inCylinder((0,0,0),(0,0,height_0),.5*width)

sp.makeCloud((0,0,0),(.3,.3,2),rMean=0.0083,rRelFuzz=0.1)

#sp = pack.randomDensePack( pred,spheresInCell=num_spheres,radius=R_p,rRelFuzz=0.3, returnSpherePack=True,memoDbg=True),#memoizeDb='/tmp/loosePackings11.sqlite')

sp.toSimulation(color=(0,1,1),material=Sand)

triax=TriaxialStressController(

 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)

 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)

 thickness = 0,

 stressMask = 7,

 max_vel = 0.005,

 internalCompaction=True, # If true the confining pressure is generated by growing particles

)

newton=NewtonIntegrator(damping=0.2)

O.engines=[

 ForceResetter(),

 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),

 InteractionLoop(

  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],

  [Ip2_FrictMat_FrictMat_FrictPhys()],

  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"

 ),

 FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section

 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),

 triax,

 newton

]

triax.goal1=triax.goal2=triax.goal3=-10000

while 1:

  O.run(1000, True)

  unb=unbalancedForce()

  if unb<0.001 and abs(-10000-triax.meanStress)/10000<0.001:

    break

setContactFriction(radians(finalFricDegree))

sand=sp.toSimulation(color=(0,1,1),material=Sand)

# C.b). define different sections of sphere pack

bot = [O.bodies[s] for s in sand if O.bodies[s].state.pos[2]<rParticle*rCoff]

top = [O.bodies[s] for s in sand if O.bodies[s].state.pos[2]>height_0-rParticle*rCoff]

tot = [O.bodies[s] for s in sand if O.bodies[s].state.pos[2]<=height_0]

# C.c). detect the position of particles in top & bot layer

top_limit = 0

top_id = 0

for s in top:

 if s.state.pos[2]>=top_limit:

  top_limit = s.state.pos[2]

  top_id = s.id

  bot_limit = height_0

  bot_id = 0

for s in bot:

 if s.state.pos[2]<=bot_limit:

  bot_limit = s.state.pos[2]

  bot_id = s.id

# C.d). create facet wall around particle packing

facets = []

nw = 45

nh = 15

rCyl2 = .5*width / cos(pi/float(nw))

for r in xrange(nw):

 for h in xrange(nh):

  v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),

  rCyl2*sin(2*pi*(r+0)/float(nw)), height_0*(h+0)/float(nh) )

  v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),

  rCyl2*sin(2*pi*(r+1)/float(nw)), height_0*(h+0)/float(nh) )

  v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),

  rCyl2*sin(2*pi*(r+1)/float(nw)), height_0*(h+1)/float(nh) )

  v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),

  rCyl2*sin(2*pi*(r+0)/float(nw)), height_0*(h+1)/float(nh) )

  f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat)

  f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat)

  facets.extend((f1,f2))

wall = O.bodies.append(facets)

# C.e). define different sections of facet wall

for b in wall:

 O.bodies[b].state.blockedDOFs = 'xyzXYZ'

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

# C.f). create top & bot facet plate

facets3 = []

nw=45

rCyl2 = (.6*width+2*rParticle) / cos(pi/float(nw))

for r in xrange(nw):

 if r%2==0:

  v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),

  rCyl2*sin(2*pi*(r+0)/float(nw)), height_0 )

  v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),

  rCyl2*sin(2*pi*(r+1)/float(nw)), height_0 )

  v3 = Vector3( rCyl2*cos(2*pi*(r+2)/float(nw)),

  rCyl2*sin(2*pi*(r+2)/float(nw)), height_0 )

  v4 = Vector3( 0, 0, height_0 )

  f1 = facet((v1,v2,v4),color=(0,0,0),material=frictMat)

  f2 = facet((v2,v3,v4),color=(0,0,0),material=frictMat)

  facets3.extend((f1,f2))

  topcap = O.bodies.append(facets3)

  facets3 = []

for r in xrange(nw):

 if r%2==0:

  v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),

  rCyl2*sin(2*pi*(r+0)/float(nw)), 0 )

  v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),

  rCyl2*sin(2*pi*(r+1)/float(nw)), 0 )

  v3 = Vector3( rCyl2*cos(2*pi*(r+2)/float(nw)),

  rCyl2*sin(2*pi*(r+2)/float(nw)), 0 )

  v4 = Vector3( 0, 0, 0 )

  f1 = facet((v1,v2,v4),color=(0,0,0),material=frictMat)

  f2 = facet((v2,v3,v4),color=(0,0,0),material=frictMat)

  facets3.extend((f1,f2))

  botcap = O.bodies.append(facets3)

# C.g). define top & bot wall id

for s in topcap:

 top_id = s

 bot_id = 0

for s in botcap:

 bot_id = s

# D.h). calculate porosity

V_sand = 0

num_sand = 0

for b in sand:

 r = O.bodies[b].shape.radius

 V_sand += 1.3333333*3.1416*r*r*r

 num_sand +=1

 porosity = 1-V_sand/(.25*width*width*3.1416*height_0)

 print 'v_sand= ',V_sand,' number of sand: ',num_sand,'porosity is: ',porosity

 O.pause()

triax.stressMask=2

triax.goal1=triax.goal3=0

triax.internalCompaction=False

triax.wall_bottom_activated=False

#load

triax.goal2=-11000; O.run(2000,1)

#unload

triax.goal2=-10000; O.run(2000,1)

#load

triax.goal2=-11000; O.run(2000,1)

e22=triax.strain[1]

#unload

triax.goal2=-10000; O.run(2000,1)

e22=e22-triax.strain[1]

modulus = 1000./abs(e22)

#######################################################################

# D. DEFINING ADD-ON FUNCTIONS #

#######################################################################

# D.a). a function for saving variables

def plotAddData():

 f1 = sum(O.forces.f(b)[2] for b in topcap)

 f2 = sum(O.forces.f(b)[2] for b in botcap)

 f11 = sum(O.forces.f(b.id)[2] for b in top)

 f22 = sum(O.forces.f(b.id)[2] for b in bot)

 fa = abs(.5*(f2-f1))

 fa1 = abs(.5*(f22-f11))

 e = (top[0].state.displ()[2] - bot[0].state.displ()[2])/height_0-e_ini

 f4 = 0

 r_cum = 0

 count = 0

 Area = 0

for f in membrane_grid:

 f.shape.color = Vector3(0,0,0)

 x = f.state.pos[0]

 y = f.state.pos[1]

 dist = math.sqrt(x*x+y*y)

 n = Vector3(x/dist,y/dist,0)

 z = f.state.pos[2]

 a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e)

if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]:

 count += 1

 r_local = dist

 r_cum += r_local

 r_avg = r_cum/count-rParticle

 Area = r_avg*2*3.1416*(O.bodies[top_id].state.pos[2] -

 O.bodies[bot_id].state.pos[2])

 area_avg = Area/count

 s = fa/(3.1416*r_avg*r_avg)

 s1 = fa1/(3.1416*r_avg*r_avg)

for b in membrane_grid:

 x = b.state.pos[0]

 y = b.state.pos[1]

 dist = math.sqrt(x*x+y*y)

 n = Vector3(x/dist,y/dist,0)

 a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e)

 z = b.state.pos[2]

 if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]:

  f_local = O.forces.f(b.id)

 length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2])

 cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length

 p_normal = (length*cos_theta/a)

 f4 += (p_normal)

 p = abs(f4/count/1000)

 h = O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2]

 VV = h*r_avg*r_avg*3.1416

 dV = VV-V_ini

 ev = -((O.bodies[top_id].state.pos[2] -

 O.bodies[bot_id].state.pos[2])*r_avg*r_avg*3.1416-V_ini)/V_ini

 er = (r_avg-R_avg)/R_avg

 if (abs(e*100)>thresholdstrain*100):

  O.pause()

plot.addData(

 i = O.iter,

 q = (abs(s)-p*1000)/1000,

 q1 = (abs(s1)-conStress)/1000,

 p = p,

 u = conStress/1000-p,

 e = -e*100,

 ev = ev*100,

 )

for b in tot:

 b.shape.color=scalarOnColorScale(b.state.rot().norm(),0,pi/2.)

#return (dV,e)

# D.b). a function for adding force (servo-controlled of lateral wall)

def addforce():

 h_sample = O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2]

#print 'height is ',h_sample

 r_cum = 0

 count = 0

 f4 = 0

for b in membrane_grid:

 x = b.state.pos[0]

 y = b.state.pos[1]

 dist = math.sqrt(x*x+y*y)

 n = Vector3(x/dist,y/dist,0)

 a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e)

 z = b.state.pos[2]

 if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]:

  f_local = O.forces.f(b.id)

 length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2])

 cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length

 p_normal = (length*cos_theta/a)

 f4 += (p_normal)

 count += 1

 p = abs(f4/count/1000)

for f in membrane_grid:

 f.shape.color = Vector3(0,0,0)

 x = f.state.pos[0]

 y = f.state.pos[1]

 dist = math.sqrt(x*x+y*y)

 n = Vector3(x/dist,y/dist,0)

 z = f.state.pos[2]

 a = 2*3.1416*dist/(360/alpha)*6*rParticle*(1+e)

 if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]:

  r_local = dist

  r_cum += r_local

  r_avg = r_cum/count-rParticle

  Volume = r_avg*r_avg*3.1416*h_sample

  delV = Volume - V_ini

for b in topcap:

 O.bodies[b].state.blockedDOFs = 'xyzXYZ'

 O.bodies[b].state.vel = (0,0,-vel_a)

for b in botcap:

 O.bodies[b].state.blockedDOFs = 'xyzXYZ'

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

for f in membrane_grid:

 f.shape.color = Vector3(0,0,0)

 x = f.state.pos[0]

 y = f.state.pos[1]

 dist = math.sqrt(x*x+y*y)

 n = Vector3(x/dist,y/dist,0)

 z = f.state.pos[2]

 a = 2*3.1416*dist/(360/alpha)*6*rParticle#*(1+e)

 if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]:

  f.state.blocked = 'z'

  f.state.vel = (dist/h_sample)*vel_a*n

  if delV>0:

   f.state.vel = -3.0*(dist/h_sample)*vel_a*n

  else:

   f.state.vel = 3.0*(dist/h_sample)*vel_a*n

else:

   f.state.vel = 0*n

# D.c). a function for recording data

def checkrecord():

 plot.saveDataTxt('results_'+key)

# D.d). a function used for consolidation

def confining():

 e_ini = (top[0].state.displ()[2] - bot[0].state.displ()[2])/height_0

 f1 = sum(O.forces.f(b)[2] for b in topcap)

 f2 = sum(O.forces.f(b)[2] for b in botcap)

 f4 = 0

 r_cum = 0

 count = 0

 a = 2*3.1416*(.5*width+rParticle)/(360/alpha)*6*rParticle

for f in membrane_grid:

 f.shape.color = Vector3(0,0,0)

 x = f.state.pos[0]

 y = f.state.pos[1]

 dist = math.sqrt(x*x+y*y)

 n = Vector3(x/dist,y/dist,0)

 a = 2*3.1416*dist/(360/alpha)*6*rParticle

 z = f.state.pos[2]

 f.state.vel = -vel_ini_r*n

for b in topcap:

 O.bodies[b].state.blockedDOFs = 'xyzXYZ'

 O.bodies[b].state.vel = (0,0,-vel_ini_a)

for b in botcap:

 O.bodies[b].state.blockedDOFs = 'xyzXYZ'

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

for b in membrane_grid:

 x = b.state.pos[0]

 y = b.state.pos[1]

 dist = math.sqrt(x*x+y*y)

 n = Vector3(x/dist,y/dist,0)

 a = 2*3.1416*dist/(360/alpha)*6*rParticle

 z = b.state.pos[2]

 if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]:

  f_local = O.forces.f(b.id)

  length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2])

  cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length

  p_normal = (length*cos_theta/a)

  f4 += (p_normal)

  r_cum += dist

  count += 1

  r_avg = r_cum/count-rParticle

  fa = abs(.5*(f2-f1))

  p_r = f4/count/1000

  p_a = fa/(3.1416*0.25*width*width)/1000

  e_ini2 = ((O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2])- height_0)/height_0

  V_ini = (O.bodies[top_id].state.pos[2] -

  O.bodies[bot_id].state.pos[2])*r_avg*r_avg*3.1416

  R_avg = r_avg

  H_ini = O.bodies[top_id].state.pos[2] - O.bodies[bot_id].state.pos[2]

  porosity = 1-V_sand/((R_avg)*(R_avg)*3.1416*H_ini)

#return (p_r,p_a,e_ini,V_ini,R_avg,H_ini,porosity)

flow.dead=0

flow.defTolerance=0.3

flow.meshUpdateInterval=200

flow.useSolver=3

flow.permeabilityFactor=1

flow.viscosity=10

flow.bndCondIsPressure=[0,0,1,1,0,0]

flow.bndCondValue=[0,0,1,0,0,0]

flow.boundaryUseMaxMin=[0,0,0,0,0,0]

O.dt=0.1e-3

O.dynDt=False

O.run(1,1)

Qin = flow.getBoundaryFlux(2)

Qout = flow.getBoundaryFlux(3)

permeability = abs(Qin)/1.e-4 #size is one, we compute K=V/H

print "Qin=",Qin," Qout=",Qout," permeability=",permeability

flow.bndCondIsPressure=[0,0,0,1,0,0]

flow.bndCondValue=[0,0,0,0,0,0]

newton.damping=0

Cv=permeability*modulus/1e4

zeroTime=O.time

zeroe22 = - triax.strain[1]

dryFraction=0.05 #the top layer is affected by drainage on a certain depth, we account for it here

drye22 = 1000/modulus*dryFraction

wetHeight=1*(1-dryFraction)

# D.e). a function for stablization

def stable():

  for f in membrane_grid:

   x = f.state.pos[0]

   y = f.state.pos[1]

   dist = math.sqrt(x*x+y*y)

   n = Vector3(x/dist,y/dist,0)

   a = 2*3.1416*dist/(360/alpha)*6*rParticle

   z = f.state.pos[2]

   f.state.blockedDOFs = 'xyzXYZ'

   f.state.vel = 0*n

for b in topcap:

 O.bodies[b].state.blockedDOFs = 'xyzXYZ'

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

for b in botcap:

 O.bodies[b].state.blockedDOFs = 'xyzXYZ'

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

def compress():

 for b in wall:

  O.bodies[b].state.blockedDOFs = 'xyzXYZ'

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

#######################################################################

# E. APPLYING CONFINING STRESS TO FLEXIBLE MEMBRANE #

#######################################################################

# E.a). adding corrosponding python function

O.engines=O.engines+[

PyRunner(command='compress()',iterPeriod=1),

PyRunner(command='plotAddData',iterPeriod=1)

]

# E.b). compress until target porosity

vel_ini_a = rate*height_0

vel_ini_r = rate*height_0

for b in topcap:

 O.bodies[b].state.blockedDOFs = 'xyzXYZ'

 O.bodies[b].state.vel = (0,0,-vel_ini_a)

for b in botcap:

 O.bodies[b].state.blockedDOFs = 'xyzXYZ'

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

while 1:

 O.run(100,True)

 h = (O.bodies[top_id].state.pos[2]-O.bodies[bot_id].state.pos[2])

 V = h*0.25*width*width*3.1416

 porosity = 1-V_sand/V

 print 'porosity: ',porosity, ' height: ', h

 if (porosity <= targetPorosity):

  print 'compression stage finished!'

  break

h_ini = h # record height after compression

# E.c). ADD CONFINING STRESS

# E.c.1). remove facet wall

for b in wall:

 O.bodies.erase(b)

# E.c.2). create membrane around particle packing by reading mesh file

shiftfactor = O.bodies[bot_id].state.pos[2]-((height-h_ini)/2)

nodesIds,cylIds,pfIds = gtsPFacet( 'Mesh_cylinder.gts',shift=(0,0,shiftfactor),scale=1,radius=rParticle,wire=False,fixed=False, materialNodes='gridNodeMat',material='pFacetMat',color=Vector3(0.5,1,0.5)

)

# E.c.3). define different sections of membrane

membrane_grid = [O.bodies[s] for s in nodesIds ]

membrane_pfacet = [O.bodies[s] for s in pfIds]

# E.c.4). run one interaction

for f in membrane_grid:

 f.shape.color = Vector3(0,0,0)

 x = f.state.pos[0]

 y = f.state.pos[1]

 dist = math.sqrt(x*x+y*y)

 n = Vector3(x/dist,y/dist,0)

 z = f.state.pos[2]

if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]:

    f.state.vel = -vel_ini_r*n

    O.engines[2].physDispatcher.functors[1].setCohesionNow=True

while 1:

  O.run(1,True)

  break

# E.c.5). redefine engine

#**********************************************************************#

O.engines=[

ForceResetter(),

InsertionSortCollider([

Bo1_Sphere_Aabb(),

Bo1_PFacet_Aabb(),

Bo1_Facet_Aabb(),

Bo1_GridConnection_Aabb()

]),

InteractionLoop(

[

Ig2_Sphere_Sphere_ScGeom6D(),

Ig2_GridNode_GridNode_GridNodeGeom6D(),

Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),

Ig2_Sphere_PFacet_ScGridCoGeom(),

Ig2_Facet_Sphere_ScGeom6D()

],

[

Ip2_FrictMat_FrictMat_FrictPhys(),

Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")

],

[

Law2_ScGeom_FrictPhys_CundallStrack(),

Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm

=True,always_use_moment_law=False,label='cohesiveLaw'),

Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),

Law2_ScGridCoGeom_FrictPhys_CundallStrack(),

],label="iloop"

),

GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8),

NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton'),

PyRunner(command='confining()',iterPeriod=1),

PyRunner(command='plotAddData',iterPeriod=1)

]

#**********************************************************************#

# set final friction angle, enable cohesion

setContactFriction(radians(finalFricDegree))

O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True

O.engines[2].physDispatcher.functors[1].setCohesionNow=True

O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False

# E.c.6). confining

# some initial parameters

p_a = 0

p_r = 0

e_ini = 0

V_ini = 0

R_avg = 0

H_ini = 0

porosity = 0

# velocity

vel_ini_a = 0.05*rate*height_0

vel_ini_r = 0.05*rate*height_0

# loops (fast-slow) for reaching target confining stress

while 1:

 O.run(10,True)

 (p_r,p_a,e_ini,V_ini,R_avg,H_ini,porosity)=confining()

 p_r = abs(p_r)

 pressure = max(p_r,p_a)

 dif = p_r-p_a

 if (p_a > isoStress/1000):

  vel_ini_a = -abs(vel_ini_a)

 if (p_r > isoStress/1000):

  vel_ini_r = -abs(vel_ini_r)

else:

  vel_ini_r = abs(vel_ini_r)

   #elif (p_a <= isoStress/1000):

if (p_r > isoStress/1000):

   vel_ini_a = abs(vel_ini_a)

   vel_ini_r = -abs(vel_ini_r)

else:

  if (pressure<0.9*isoStress/1000):

   if dif > 5:

    vel_ini_a = 1.05*abs(vel_ini_a)

   #elif dif < -5:

    vel_ini_r = 1.05*abs(vel_ini_r)

if (pressure>=0.85*isoStress/1000 and pressure<=isoStress/1000):

 if dif > 1:

  if dif > 5:

   vel_ini_a = 1.5*abs(vel_ini_a)

 else:

  vel_ini_a = 1.01*abs(vel_ini_a)

elif dif < -1:

   if dif < -5:

    vel_ini_r = 1.5*abs(vel_ini_r)

else:

 vel_ini_r = 1.01*abs(vel_ini_r)

 mean = (p_r+p_a)/2

 unb=unbalancedForce()

 print 'p= ',p_r,' q= ',p_a,' porosity= ',porosity,' unbalanced force: ',unb

if abs(isoStress/1000-mean)/(isoStress/1000)<0.15 and abs(dif) <15:

  print 'initial strain: ',e_ini

  print 'initial volume: ',V_ini

  print 'Confining stage I is finished!'

    #break

while 1:

 for f in membrane_grid:

  f.state.blockedDOFs = 'xyzXYZ'

O.run(1,True)

(p_r,p_a,e_ini,V_ini,R_avg,H_ini,porosity)=confining()

p_r = abs(p_r)

pressure = max(p_r,p_a)

dif = p_r-p_a

if (p_a > isoStress/1000):

 vel_ini_a = -abs(vel_ini_a)

if (p_r > isoStress/1000):

 vel_ini_r = -abs(vel_ini_r)

else:

 vel_ini_r = abs(vel_ini_r)

#elif (p_a <= isoStress/1000):

if (p_r > isoStress/1000):

  vel_ini_a = abs(vel_ini_a)

  vel_ini_r = -abs(vel_ini_r)

else:

 if (pressure<0.9*isoStress/1000):

  if dif > 5:

   vel_ini_a = 1.05*abs(vel_ini_a)

#elif dif < -5:

 vel_ini_r = 1.05*abs(vel_ini_r)

if (pressure>=0.9*isoStress/1000 and pressure<=isoStress/1000):

 if dif > 1:

  if dif > 5:

   vel_ini_a = 1.1*abs(vel_ini_a)

else:

 vel_ini_a = 1.01*abs(vel_ini_a)

#elif dif < -1:

 if dif < -5:

  vel_ini_r = 1.1*abs(vel_ini_r)

 else:

  vel_ini_r = 1.01*abs(vel_ini_r)

 mean = (p_r+p_a)/2

 unb=unbalancedForce()

print 'p= ',p_r,' q= ',p_a,' porosity= ',porosity,' unbalanced force: ',unb

if abs(isoStress/1000-mean)/(isoStress/1000)<0.05 and abs(dif) <10:

 print 'initial strain: ',e_ini

 print 'initial volume: ',V_ini

 print 'Confining stage II is finished!'

 #break

print 'V_ini= ',V_ini

# E.c.7). stablize

#**********************************************************************#

O.engines=[

ForceResetter(),

InsertionSortCollider([

Bo1_Sphere_Aabb(),

Bo1_PFacet_Aabb(),

Bo1_Facet_Aabb(),

Bo1_GridConnection_Aabb()

]),

InteractionLoop(

[

Ig2_Sphere_Sphere_ScGeom6D(),

Ig2_GridNode_GridNode_GridNodeGeom6D(),

Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),

Ig2_Sphere_PFacet_ScGridCoGeom(),

Ig2_Facet_Sphere_ScGeom6D()

],

[

Ip2_FrictMat_FrictMat_FrictPhys(),

Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")

],

[

Law2_ScGeom_FrictPhys_CundallStrack(),

Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm

=True,always_use_moment_law=False,label='cohesiveLaw'),

Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),

Law2_ScGridCoGeom_FrictPhys_CundallStrack(),

],label="iloop"

),

GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8),

NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton'),

PyRunner(command='stable()',iterPeriod=1),

PyRunner(command='plotAddData',iterPeriod=1)

]

#**********************************************************************#

O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True

O.engines[2].physDispatcher.functors[1].setCohesionNow=True

O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False

vel_a = abs(vel_ini_a)

while 1:

 O.run(100,True)

 unb=unbalancedForce()

 print 'unbalanced force: ',unb

 e_ini = (top[0].state.displ()[2] - bot[0].state.displ()[2])/height_0

 f1 = sum(O.forces.f(b)[2] for b in topcap)

 f2 = sum(O.forces.f(b)[2] for b in botcap)

 f4 = 0

 r_cum = 0

 count = 0

for b in membrane_grid:

 x = b.state.pos[0]

 y = b.state.pos[1]

 dist = math.sqrt(x*x+y*y)

 n = Vector3(x/dist,y/dist,0)

 a = 2*3.1416*dist/(360/alpha)*6*rParticle#*(1+e)

 z = b.state.pos[2]

if z<=O.bodies[top_id].state.pos[2] and z>=O.bodies[bot_id].state.pos[2]:

    f_local = O.forces.f(b.id)

    length = math.sqrt(f_local[0]*f_local[0]+f_local[1]*f_local[1]+f_local[2]*f_local[2])

    cos_theta = (n[0]*f_local[0]+n[1]*f_local[1]+n[2]*f_local[2])/length

p_normal = (length*cos_theta/a)

f4 += (p_normal)

r_cum += dist

count += 1

r_avg = r_cum/count-rParticle

fa = abs(.5*(f2-f1))

p_r = f4/count/1000

p_a = fa/(3.1416*r_avg*r_avg)/1000

print 'pr=', p_r, ' pa=',p_a

if unb <= thresholdvalue:
#break

#######################################################################

# F. APPLYING DEVIATOR LOADING #

#######################################################################

# F.a). redifine engines

    O.engines=[
    ForceResetter(),
    InsertionSortCollider([
                   Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargefactor),
                   Bo1_PFacet_Aabb(),
                   Bo1_Facet_Aabb(),
                   Bo1_GridConnection_Aabb()
 ]),

InteractionLoop(

[

Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=enlargefactor),

Ig2_GridNode_GridNode_GridNodeGeom6D(),

Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),

Ig2_Sphere_PFacet_ScGridCoGeom(),

Ig2_Facet_Sphere_ScGeom6D()

],

[

Ip2_FrictMat_FrictMat_FrictPhys(),

Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")

],

[

Law2_ScGeom_FrictPhys_CundallStrack(),

Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm

=True,always_use_moment_law=False,label='cohesiveLaw'),

Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),

Law2_ScGridCoGeom_FrictPhys_CundallStrack(),

],label="iloop"

),

GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8),

NewtonIntegrator(gravity=(0,0,0),damping=0.3,label='newton')

]

O.engines[2].lawDispatcher.functors[1].always_use_moment_law=True

O.engines[2].physDispatcher.functors[1].setCohesionNow=True

O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False

O.engines=O.engines+[

PyRunner(command='addforce()',iterPeriod=1,label='force'),

PyRunner(command='plotAddData()',iterPeriod=1,label='recorder'),

PyRunner(command='checkrecord()',realPeriod=10,label='checker')

]

# F.b). define the velocity of membrane walls to maintain the volume constant condition

vel_a = final_rate*abs(vel_ini_a)

vel_r = vel_a*.5*width/height

Vel_r = vel_r

conStress = p_r*1000

#######################################################################

# G. UTILITIES #

#######################################################################

# G.a). time step (recommanded by YADE)

O.dt=0.3*PWaveTimeStep()

t = O.dt

# G.b). funtion for plot

plot.plots={'e':('q','p'),'e':('u','ev')}

plot.plot()

O.saveTmp()

O.timingEnabled=1

def consolidation(Tv): #see your soil mechanics handbook...

 U=1

 for k in range(50):

  M=pi/2*(2*k+1)

  U=U-2/M**2*exp(-M**2*Tv)

 return U

triax.goal2=-11000

from yade import plot

def history():

   plot.addData(e22=-triax.strain[1]-zeroe22,e22_theory=drye22+(1-dryFraction)*consolidation((O.time-zeroTime)*Cv/wetHeight**2)*1000./modulus,t=O.time,p=flow.getPorePressure((0.5,0.1,0.5)),s22=-triax.stress(3)[1]-10000)

from yade import plot

plot.plots={'t':('e22','e22_theory',None,'s22','p')}

plot.plot()

O.saveTmp()

O.timingEnabled=1

from yade import timing

print "starting oedometer simulation"

O.run(200,1)

timing.stats()

####################################################################

# G. RUN #

####################################################################

print '========================='

print "start triaxial simulation"

print '========================='

O.run()

Cheers
Sam

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said :
#1

Hello Sam,

>The The code run without any issue but it's not working properly.

I'm sure we can help, but you need to help us help you first. Please review [1] and then supply the supplementary information.

Have a nice day!

Robert

[1]https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Sam (sambahmani) said :
#2

Hi Robert,

Thanks for the reply. When I run the code, I get this one:

Welcome to Yade 2018.02b-290bf6a54e~xenial
TCP python prompt on localhost:9000, auth cookie `yaksus'
XMLRPC info provider on http://localhost:21000
Running script MINPFV.py

You can find the MWE pf the code as follows:

# import corrosponding yade modulus

from yade import pack

from yade import geom

from yade import plot

from yade import ymport

from yade import qt

from yade.gridpfacet import *

import gts, os.path, locale

locale.setlocale(locale.LC_ALL, 'en_US.UTF-8')

import sys,os

from yade import plot,pack,export

sys.path.append(os.pardir)

########################################################################

# A. DEFINING VARIABLES, MATERIALS' PROPERTIES AND PACKING #

########################################################################

# A.a). define variables

key = 'Triaxial_Undrained' # file name to be saved

young=550e6 # normal contact stiffness

compFricDegree = 1.8 # initial contact friction during the confining phase

finalFricDegree = 43 # contact friction during the deviatoric loading

poisson = 0.3 # shear-to-normal stiffness ratio

isoStress = 110000 # confining stress

conStress = 100000 # confining stress for deviatoric loading stage

width = 1.4e-1 # sample width

height = 2.8e-1 # target sample height(after consolidation)

height_0 = 3.2e-1 # initial sample height

num_spheres=500 # number of spheres

R_p = 0.0084 # mean particle radius

rCoff = 10 # thickness of top and bot sphere cap (based on rParticle)

rParticle = 0.02e-1 # membrane grid seed size

alpha = 8

rate = 0.1 # loading rate (strain rate)

damp = 0.3 # damping coefficient

targetPorosity = 0.43 # target porosity

thresholdvalue = 0.05 # threshold unbalance force

final_rate = 0.1 # strain rate for deviator loading

thresholdstrain = 0.06 # threshold axial strain for terminate

enlargefactor = 1.00

# A.b). create materials for sand spheres and plates

Sand = O.materials.append(CohFrictMat(

young=young,poisson=poisson,frictionAngle=radians(compFricDegree),

alphaKr=0.25,alphaKtw=0,etaRoll=0.005,etaTwist=0,

normalCohesion=5e6,shearCohesion=5e6,

momentRotationLaw=True,density=2650,label='spheres'

))

# A.c). create membrane materials

207

GridMat = O.materials.append(CohFrictMat(

young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),

alphaKr=0,alphaKtw=0,etaRoll=0,etaTwist=0,

normalCohesion=1e9,shearCohesion=1e9,

momentRotationLaw=True,label='gridNodeMat'

))

pFacetMat = O.materials.append(FrictMat(

young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),label='pFacetMat'

))

# A.d). create TOP & BOT plate materials

frictMat = O.materials.append(FrictMat(

young=100e6,poisson=0.3,density=2650,frictionAngle=radians(0),label='frictMat'

))

###################################################################

# B. DEFINING GLOBAL ENGINES #

###################################################################

#**********************************************************************#

O.engines=[

ForceResetter(),

InsertionSortCollider([

Bo1_Sphere_Aabb(),

Bo1_PFacet_Aabb(),

Bo1_Facet_Aabb(),

Bo1_GridConnection_Aabb()

]),

InteractionLoop(

[

Ig2_Sphere_Sphere_ScGeom6D(),

Ig2_GridNode_GridNode_GridNodeGeom6D(),

Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),

Ig2_Sphere_PFacet_ScGridCoGeom(),

Ig2_Facet_Sphere_ScGeom6D()

],

[

Ip2_FrictMat_FrictMat_FrictPhys(),

Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")

],

[

Law2_ScGeom_FrictPhys_CundallStrack(),

Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm

=True,always_use_moment_law=False,label='cohesiveLaw'),

Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),

Law2_ScGridCoGeom_FrictPhys_CundallStrack(),

],label="iloop"

),

GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=25,timestepSafetyCoefficient=0.8),

NewtonIntegrator(gravity=(0,0,0),damping=0.1,label='newton')

]

#**********************************************************************#

O.engines[2].lawDispatcher.functors[1].always_use_moment_law=False

O.engines[2].physDispatcher.functors[1].setCohesionOnNewContacts=False

#######################################################################

# C. GENERATING PACKING #

#######################################################################

# C.a). generate random dense sphere pack

sp=pack.SpherePack()

pred = pack.inCylinder((0,0,0),(0,0,height_0),.5*width)

sp.makeCloud((0,0,0),(.3,.3,2),rMean=0.0083,rRelFuzz=0.1)

#sp = pack.randomDensePack( pred,spheresInCell=num_spheres,radius=R_p,rRelFuzz=0.3, returnSpherePack=True,memoDbg=True),#memoizeDb='/tmp/loosePackings11.sqlite')

sp.toSimulation(color=(0,1,1),material=Sand)

triax=TriaxialStressController(

 maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)

 finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)

 thickness = 0,

 stressMask = 7,

 max_vel = 0.005,

 internalCompaction=True, # If true the confining pressure is generated by growing particles

)

newton=NewtonIntegrator(damping=0.2)

O.engines=[

 ForceResetter(),

 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),

 InteractionLoop(

  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],

  [Ip2_FrictMat_FrictMat_FrictPhys()],

  [Law2_ScGeom_FrictPhys_CundallStrack()],label="iloop"

 ),

 FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section

 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),

 triax,

 newton

]

triax.goal1=triax.goal2=triax.goal3=-10000

while 1:

  O.run(1000, True)

  unb=unbalancedForce()

  if unb<0.001 and abs(-10000-triax.meanStress)/10000<0.001:

    break

Revision history for this message
Robert Caulk (rcaulk) said :
#3

Hello Sam,

That is not an error. Everything appears to be working properly.

Cheers,

Robert

Revision history for this message
Sam (sambahmani) said :
#4

Robert,

you are completely right and That's my question. I run the model with 50 particles 3 days ago and nothing happened. that's little bit weird, isn't it?

Cheers
Sam

Revision history for this message
Robert Caulk (rcaulk) said :
#5

Hello Sam,

We are happy to help you. But you have not provided enough information for us to do so. It seems you missed [1] so I will summarize and rewrite [1] here for you:

What are you expecting to happen?
How do you know "nothing happened?"
Which example script are you basing your script on?
What have you changed in comparison between your script and the example script?

Cheers,

Robert

[1]https://www.yade-dem.org/wiki/Howtoask

Revision history for this message
Launchpad Janitor (janitor) said :
#6

This question was expired because it remained in the 'Needs information' state without activity for the last 15 days.