yade-users team mailing list archive

[Elham Hosseinkhani <question688439@xxxxxxxxxxxxxxxxxxxxx>]

Question #688439 on Yade changed:
https://answers.launchpad.net/yade/+question/688439

Elham Hosseinkhani posted a new comment:
Hello,
Unfortunately,I couldn't get the results that I expected based on the papers. I changed my script and modify the script in YADE examples. I added CohesionMomentLaw and necessary parameters like alphaKr and etaRoll. I get PrincipalStressRatio  about 2.4 with this parameters in this script. But I expect more PrincipalStressRatio after using CohesionMomentLaw.

Could you please run this script and check PrincipalStressRatio? 
Is there any problem that i couldn't get changes in PrincipalStressRatio and peak friction angle? 

# -*- coding: utf-8 -*-
from yade import pack


nRead=readParamsFromTable(
	num_spheres=10000,
	young = 500e6,
	poisson = 0.4,
	compFricDegree = 30, 
	alphaKr = 1.25,
	etaRoll = 0.55,
	finalFricDegree = 30,
	targetPorosity = 0.4,
	confining =-300000,
	unknownOk=True
)
from yade.params import table
num_spheres=table.num_spheres
young = table.young
poisson = table.poisson
compFricDegree = table.compFricDegree 
alphaKr = table.alphaKr
etaRoll = table.etaRoll
finalFricDegree = table.finalFricDegree 
targetPorosity = table.targetPorosity
confining = table.confining
rate=-0.02 
damp=0.2 

key='_P='+str(confining)
stabilityThreshold=0.01 
mn,mx=Vector3(0,0,0),Vector3(1,2,1) 

O.materials.append(CohFrictMat(young=young,poisson=poisson,density=2600,frictionAngle=radians(compFricDegree),normalCohesion=1e6,shearCohesion=1e6,alphaKr=alphaKr,alphaKtw=1e-10,momentRotationLaw=True,etaRoll=etaRoll,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='walls'))

walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()

clumps=False 
if clumps:

 volume = (mx[0]-mn[0])*(mx[1]-mn[1])*(mx[2]-mn[2])
 mean_rad = pow(0.09*volume/num_spheres,0.3333)
 c1=pack.SpherePack([((-0.2*mean_rad,0,0),0.5*mean_rad),((0.2*mean_rad,0,0),0.5*mean_rad)])
 sp.makeClumpCloud(mn,mx,[c1],periodic=False)
 sp.toSimulation(material='spheres')
 O.bodies.updateClumpProperties()
else:
 sp.makeCloud(mn,mx,-1,0.3333,num_spheres,False, 0.95,seed=1) #"seed" make the "random" generation always the same
 O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
	maxMultiplier=1.+2e5/young, 
	finalMaxMultiplier=1.+2e4/young, 
	thickness = 0,
	stressMask = 7,
	internalCompaction=True, 
)

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_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label="cohesiveIp")],
		[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(
		useIncrementalForm=True,
		always_use_moment_law=False,
		label='cohesiveLaw')]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'),
	newton
]

Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

triax.goal1=triax.goal2=triax.goal3=confining

while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print ('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<stabilityThreshold and abs(confining-triax.meanStress)/abs(confining)<0.01:
    break

import sys
while triax.porosity>targetPorosity:
	compFricDegree = 0.95*compFricDegree
	setContactFriction(radians(compFricDegree))
	print ("\r Friction: ",compFricDegree," porosity:",triax.porosity)
	sys.stdout.flush()
	O.run(500,1)

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

triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))
triax.stressMask = 5
triax.goal2=rate
triax.goal1=confining
triax.goal3=confining
O.engines[2].lawDispatcher.functors[1].always_use_moment_law = True
#cohesiveLaw.always_use_moment_law=True
triax.wall_bottom_activated=False
newton.damping=0.1
O.saveTmp()

from yade import plot

def history():
	plot.addData(e11=-triax.strain[0], e22=-triax.strain[1], e33=-triax.strain[2],
			ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
			s11=-triax.stress(triax.wall_right_id)[0],
			s22=-triax.stress(triax.wall_top_id)[1],
			s33=-triax.stress(triax.wall_front_id)[2],
			PrincipalStressRatio=abs((-triax.stress(3)[1])/((-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2)),
			q=-triax.stress(triax.wall_top_id)[1]-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0,
			i=O.iter)

if 1:
  O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
else:
  O.engines[4]=PyRunner(iterPeriod=20,command='history()',label='recorder')

O.run(1000,True)
plot.plots={'e22':('PrincipalStressRatio',None,'ev')}
plot.plot()

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