yade-users team mailing list archive

[Elham Hosseinkhani <question688439@xxxxxxxxxxxxxxxxxxxxx>]

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.