New question #693129 on Yade:
https://answers.launchpad.net/yade/+question/693129
Dear all,
I am trying to model a drained triaxial test in dry sand. I used Mindlin
contact law for contact between spheres. Also, I want to consider bending
moment and activate includeMoment=True in Law2_ScGeom_MindlinPhys_Mindlin().
But i couldn't make a pack or continue shear loading after making pack and
consolidation(I see Error nan in terminal). I need to model rolling resistance
to get experimental shear strength. Therefore, i gave some value to eta for
considering plastic bending moment.
Would you please look at my script or run it and inform me about mistakes if
there are?
Thank you for your help,
Elham
# encoding: utf-8
import matplotlib; matplotlib.rc('axes',grid=True)
import matplotlib.pyplot as plt
from yade import pack
import pylab
from numpy import *
import numpy as np
from math import *
utils.readParamsFromTable(seed=1,num_spheres=1000,compFricDegree =3.0)
from yade.params import table
seed=table.seed
num_spheres=table.num_spheres
compFricDegree = table.compFricDegree
confiningS=-50000
rate=1
Young=200 #MPa
Poisson=0.3
## creat a packing with a specific particle side distribution (PSD)
psdSizes,psdCumm=[0.075,0.107,0.153,0.192,0.211,0.232,0.252,0.279,0.303,0.462],[0.035,0.065,0.246,0.5,0.6,0.7,0.8,0.9,0.965,1]
psdSizesArray=np.array(psdSizes)
psdSizesMeter=psdSizesArray*0.001 #Convert the size of particles to meter
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.002,0.004,0.002)
sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=seed)
O.materials.append(FrictMat(young=20e7,poisson=.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=20e7,poisson=.3,frictionAngle=0,density=0,label='frictionless'))
walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)
O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in
sp])
triax=TriaxialStressController(
internalCompaction=True,
goal1=confiningS,
goal2=confiningS,
goal3=confiningS,
max_vel=10,
label="triax"
)
newton=NewtonIntegrator(damping=0.4)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys(en=0.9,es=0.9,eta=0.8,krot=1)],
[Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True,label='mindlin')]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
triax,
newton
]
while 1:
O.run(1000,True)
unb=unbalancedForce()
print
('Porosity:',triax.porosity,'unb:',unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1))
if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
break
print ("Particle Distribution is Stable")
print ('Porosity:',triax.porosity)
#############################
## REACH NEW EQU. STATE ###
#############################
finalFricDegree = 26.5 # contact friction during the 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))
while 1:
O.run(1000,True)
unb=unbalancedForce()
print
('Porosity:',triax.porosity,'unb:',unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1))
if unb<0.001 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
break
triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width
print ("Time for the deviator loads")
print ('Porosity:',triax.porosity)
##########################
## Start test ###
##########################
AxialstrainMatrix1=[]
DeviatoricStressMatrix1=[]
VolumetricStrainMatrix1=[]
PrincipalStressRatioMatrix1=[]
ESecanti=[]
triax.stressMask=5
triax.goal2=-rate
triax.goal1=confiningS
triax.goal3=confiningS
key='E_'+str(Young)+'_noo_'+str(Poisson)
file=open('Results'+key+'.txt',"w")
while 1:
O.run(1000,True)
AxialstrainMatrix1.append((-triax.strain[1]))
DeviatoricStressMatrix1.append((-triax.stress(triax.wall_top_id)[1]-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)*1e-3)
PrincipalStressRatioMatrix1.append(abs((-triax.stress(3)[1])/((-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)))
VolumetricStrainMatrix1.append((triax.strain[0]+triax.strain[1]+triax.strain[2]))
ESecanti.append((-triax.stress(triax.wall_top_id)[1]-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)*1e-3/(-triax.strain[1]))
print ('eyy=',triax.strain[1],'Sxx=',
triax.stress(0)[0],'Syy=',triax.stress(3)[1], 'Szz=',triax.stress(4)[2],'e=',
(triax.porosity)/(1-triax.porosity))
if triax.strain[1] <-0.03:
break
file.write(str(AxialstrainMatrix1)+" "+str(DeviatoricStressMatrix1)+"
"+str(PrincipalStressRatioMatrix1)+" "+str(VolumetricStrainMatrix1)+" "+"\n")
file.close()
--
You received this question notification because your team yade-users 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
