← Back to team overview

yade-users team mailing list archive

[Question #698135]: PFV model is not working without any issue

 

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

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


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.