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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
