yade-users team mailing list archive

[Luc Sibille <question688439@xxxxxxxxxxxxxxxxxxxxx>]

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

Luc Sibille proposed the following answer:
Did you check if there is actually a rolling moment not nil at the 
contacts between spheres?
Luc

Le 09/02/2020 à 10:23, Elham Hosseinkhani a écrit :
> 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()
>

-- 
Luc Sibille
Université Grenoble Alpes / IUT1 de Grenoble
Laboratoire 3SR: Sols, Solides, Structures, Risques

Tel lab.: +33 (0)4 76 82 63 48
Tel IUT: +33 (0)4 76 82 53 36

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