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

Dear all,
I'm having a problem regarding setting stress mask in triaxial example. I've 
read other posts related to my problem but it still exists. I need STRAIN 
controlled loading in the 2nd phase of trixial test. I want to apply a similar 
strian rate in all 3 directions in order to achieve a mean stress of 4 MPa. 
I've set the stress mask value to 0 and 7 respectively but still the walls 
don't move. How should I modify the code to get all of the walls moving with 
the same rate. Thanks in advance for your kind answers. My script is as follows:

from yade import pack,plot
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab
############################################
###   DEFINING VARIABLES AND MATERIALS   ###
############################################

psdSizes=[.107*0.86,.132*0.86,.216*0.86,.399*0.86,.885*0.86,.97*0.86,1.113*0.86,1.359*0.86,1.633*0.86,1.993*0.86,2.433*0.86]
psdCumm=[0.*100,.0643*100,.1969*100,.1969*100,.2*100,.2988*100,.4004*100,.5515*100,.6965*100,.8476*100,1.*100]
pylab.plot(psdSizes,psdCumm)
pylab.xscale('log')
pylab.gcf().savefig('psd.pdf')
compFricDegree = 0, # contact friction during the confining phase
finalFricDegree = 26.5 # contact friction during the deviatoric loading
rate=0.1 # loading rate (strain rate)
damp=0.5 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in 
different loops (see below)
key='_Cu=3_' # put you simulation's name here
young=52e9 # contact stiffness

## create materials for spheres and plates
O.materials.append(FrictMat(young=52e9,poisson=.35,frictionAngle=radians(0),density=2000,label='spheres'))
O.materials.append(FrictMat(young=52e9,poisson=.2,frictionAngle=0,density=0,label='walls'))

## create walls around the packing
walls=utils.aabbWalls([(0,0,0),(13.244,13.244,13.244)],oversizeFactor=1)
wallIds=O.bodies.append(walls)

## use a SpherePack object to generate a random loose particles packing
sp=pack.SpherePack();
sp.particleSD2(radii=psdSizes,passing=psdCumm,numSph=15000,cloudPorosity=0.585,seed=1)
O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in 
sp])

############################
###   DEFINING ENGINES   ###
############################

triax=TriaxialStressController(
        maxMultiplier=1.0000025, # spheres growing factor (fast growth)
        finalMaxMultiplier=1.00000025, # spheres growing factor (slow growth)
        thickness = 0,
        ## switch stress/strain control using a bitmask. What is a bitmask, 
huh?!
        ## Say x=1 if stess is controlled on x, else x=0. Same for for y and z, 
which are 1 or 0.
        ## Then an integer uniquely defining the combination of all these tests 
is: mask = x*1 + y*2 + z*4
        ## to put it differently, the mask is the integer whose binary 
representation is xyz, i.e.
        ## "100" (1) means "x", "110" (3) means "x and y", "111" (7) means "x 
and y and z", etc.
        stressMask = 7,
        internalCompaction=True, # If true the confining pressure is generated 
by growing particles
)
newton=NewtonIntegrator(damping=damp)

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_FrictMat_FrictMat_MindlinPhys()],
                [Law2_ScGeom_MindlinPhys_MindlinDeresiewitz()]
        ),
        
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax,
        TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+key),
        newton
]

#######################################
###   APPLYING CONFINING PRESSURE   ###
#######################################

#the value of (isotropic) confining stress defines the target stress to be 
applied in all three directions
triax.goal1=triax.goal2=triax.goal3=100000

while 1:
  O.run(2500, True)
  #the global unbalanced force on dynamic bodies, thus excluding boundaries, 
which are not at equilibrium
  unb=unbalancedForce()
  fi = open('%s.xls'%"radius of each body",'w')
  fi.write('%s\n'%' ')
  fi.write('%s\n'%' ')
  fi.write('%s\n'%' ')
  fi.write('%s\n'%' ')  
  fi.write('%s\n'%' ')
  fi.write('%s\n'%' ')
  fi.write('%s\n'%' ')
  fi.write('%s\n'%' ')
  fi.write('%s\n'%' ')
  for m in O.bodies:
        if isinstance(m.shape,Sphere): 
                        radii=m.shape.radius
                        fi.write('%s\n'%str(radii))
  print 'unb-force: ',unb,'   mean stress: ',triax.meanStress,'  Iteration: 
',O.iter
  if unb<stabilityThreshold and abs(100000-triax.meanStress)/100000<0.01:
    break

O.save('confinedState'+key+'.yade.gz')
print "###      Isotropic state saved      ###",'     Calculation Time=   
',O.realtime

#############################
##   DEVIATORIC LOADING   ###
#############################


#We move to deviatoric loading, let us turn internal compaction off to keep 
particles sizes constant
triax.internalCompaction=False

# Change contact friction (remember that decreasing it would generate 
instantaneous instabilities)
setContactFriction(radians(finalFricDegree))

#set stress control on x and z, we will impose strain rate on y
triax.stressMask = 7
#now goal2 is the target strain rate
triax.goal2=-rate
# we define the lateral stresses during the test, here the same 10kPa as for 
the initial confinement.
triax.goal1=-rate
triax.goal3=-rate
thickness=1
#Save temporary state in live memory. This state will be reloaded from the 
interface with the "reload" button.
O.saveTmp()


-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.

_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to     : [email protected]
Unsubscribe : https://launchpad.net/~yade-users
More help   : https://help.launchpad.net/ListHelp

Reply via email to