New question #688439 on Yade: https://answers.launchpad.net/yade/+question/688439
Hello, I want to model a sample that is in bellow paper: Simulation of triaxial response of granular materials by modified DEM WANG XiaoLiang & LI JiaChun doi: 10.1007/s11433-014-5605-z December 2014 Vol. 57 No. 12: 2297–2308 I'm trying to model effects of rolling stiffness on macroscopic mechanical parameters like peak friction angle. I used Ip2_CohFrictMat_CohFrictMat_CohFrictPhys and Law2_ScGeom6D_CohFrictPhys_CohesionMoment and i wrote bellow code.The basic parameters like final friction degree, confining pressure,poisson, psd, young modulus,strain rate and porosity are the same as parameters in that paper. I didn't see any effect and changes on shear strength and peak friction angle after setting rolling stiffness(alphaKr). Actullay i couldn't simulate graphs that are in the paper with different alfaKr (for example for alphaKr =0.01, 1). 1) Is there any problem in my code that i can't simulate graphs in that paper? 2) Should i expect that different alphaKr make a difference in peak friction angle of sample? Thanks in advance for your help. ## encoding: utf-8 import matplotlib; matplotlib.rc('axes',grid=True) import matplotlib.pyplot as plt import pylab from yade import pack, qt import numpy as np from numpy import * import math import sys from yade import export from yade import timing from yade import plot import time from math import * ## Define Parameters num_spheres=10000 compFricDegree=20 finalFricDegree=30 confiningS=-100000 # [Pa] RhoW=1000 graindensity=2600 poissonRatio=0.4 youngModulus=270e6 # [Pa] AxialstrainMatrix1=[] DeviatoricStressMatrix1=[] VolumetricStrainMatrix1=[] PrincipalStressRatioMatrix1=[] qPMatrix1=[] AxialstrainMatrix2=[] DeviatoricStressMatrix2=[] VolumetricStrainMatrix2=[] PrincipalStressRatioMatrix2=[] qPMatrix2=[] alphaKrMatrix=[0.01,1] for i in range(0,len(alphaKrMatrix),1): psdSizes=[0.107,0.12,0.132,0.145,0.158,0.171,0.183,0.196,0.209,0.221,0.234] #(mm) psdCumm=[0.001,0.029,0.069,0.121,0.187,0.268,0.367,0.49,0.63,0.8,1.0] #cumulative 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.0045,0.0075,0.0045) #initial box size sp.makeCloud(minCorner=mn,maxCorner=mx,num=num_spheres,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,seed=True) sp.psd(bins=200,mass=True) O.materials.append(CohFrictMat(isCohesive=False,alphaKr=alphaKrMatrix[i],alphaKtw=0,momentRotationLaw=True,young=youngModulus,poisson=poissonRatio,frictionAngle=radians(compFricDegree),density=graindensity,label='spheres')) O.materials.append(FrictMat(young=youngModulus,poisson=poissonRatio,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=False, goal1=confiningS, goal2=confiningS, goal3=confiningS, label="triax" ) 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(traceEnergy=True),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=False,label='cohesiveLaw')] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), triax, NewtonIntegrator(damping=0.3,label="newton"), ] O.dt=utils.PWaveTimeStep() O.dynDt=False while 1: O.run(1000,True) unb=unbalancedForce() e=(triax.porosity)/(1-triax.porosity) print ('unb:',np.round(unb,2),'e: ',np.round(e,3),'meanstress:', np.round(-triax.meanStress,3),' SigmaY: ',np.round(-triax.stress(3)[1],3),' SigmaX: ',np.round(-triax.stress(1)[0],3),' SigmaZ: ',np.round(-triax.stress(5)[2],3)) if unb<0.1 and abs(confiningS-triax.meanStress)/abs(confiningS)<0.01 and e<0.634: break ########################## ## Start test ### ########################## # triaxial conditions triax.internalCompaction=False setContactFriction(radians(finalFricDegree)) triax.stressMask=5 triax.goal1=triax.goal3=confiningS triax.wall_bottom_activated=False cohesiveLaw.always_use_moment_law=True triax.goal2=-5 if i==0: while 1: O.run(1000,True) AxialstrainMatrix1.append((-triax.strain[1])*100) qPMatrix1.append((-triax.stress(triax.wall_top_id)[1]-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)/(-triax.meanStress)) 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.3: break if i==1: while 1: O.run(1000,True) AxialstrainMatrix2.append((-triax.strain[1])*100) qPMatrix2.append((-triax.stress(triax.wall_top_id)[1]-(-triax.stress(triax.wall_right_id)[0]-triax.stress(triax.wall_front_id)[2])/2.0)/(-triax.meanStress)) 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.3: break O.reset() plt.figure() plt.subplot(224) plt.xlabel('Axial Strain (%)') plt.ylabel('q/P') plt.grid(True) plt.plot(AxialstrainMatrix1,qPMatrix1,linestyle='-',color='b',linewidth=1) plt.plot(AxialstrainMatrix2,qPMatrix2,linestyle='-',color='r',linewidth=1) plt.legend(('alfa=0.01','alfa=1'),loc='lower right', shadow=True) plt.show() -- 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 : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp