Yade

effect of rolling stiffness on peak friction angle of sample

Asked by Elham Hosseinkhani

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()

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Elham Hosseinkhani
Solved:
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said : 		#1

>>I want to model a sample that is in bellow paper:
>> Is there any problem in my code that i can't simulate graphs in that paper?

Have you tried contacting the corresponding authors for assistance with this [1]? It is their responsibility to make sure their plots are reproducible. They may be using different contact laws than the ones you've selected.

[1] LI JiaChun, email: <email address hidden>; WANG XiaoLiang, email: wangxiaoLiang52086@126.com

Revision history for this message
Elham Hosseinkhani (ehosseinkhani) said : 		#2

Thanks Robert Caulk for your help,

>> Have you tried contacting the corresponding authors for assistance with this?

Actually, based on relations of contact laws that are in the paper, i found that i need these two contact laws for simulating rolling stiffness in my sample. But i'll send an email for paper authors and ask for more help.

Revision history for this message
Luc Sibille (luc-sibille) said : 		#3

Dear Elham,

Please note that the two papers [1] and [2] (below) are based on the
contact law you want to use.
It has been shown in paper [1] that the peak friction angle does not
depend on alphaKr when alphaKr is sufficiently high. Maybe you are in
this latter case...?

Regards,
Luc

[1] Discrete numerical modeling of loose soil with spherical particles
and interparticle rolling friction
RA Hosn, L Sibille, N Benahmed, B Chareyre
Granular matter 19 (1), 4
http://hal.univ-grenoble-alpes.fr/hal-01951842/document

[2] Quantitative prediction of discrete element models on complex
loading paths
Luc Sibille, Pascal Villard, Félix Darve, Rodaina Aboul Hosn
International Journal for Numerical and Analytical Methods in Geomechanics
43(5): 858-887
https://hal.archives-ouvertes.fr/hal-02059431/document

--
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

Revision history for this message
Elham Hosseinkhani (ehosseinkhani) said : 		#4

Hello Luc Sibille, and thanks for your comment.

I've studied paper [1] before and I've observed that point you said (peak friction angle does not depend on alphaKr when alphaKr is sufficiently high). But my problem is about the amount of peak friction angle before considering rolling stiffness (I mean without considering CohesionMomment) and after that. when I use just FrictPhys_CundallStrack as the contact law, with the Inter-Particle friction angle of 30 degree, I can't get peak friction angle of experimental data that is about 40 degree. Because of this problem, I tried to model a sample with more stiffness until getting more peak friction angle. After using CohesionMoment and setting alpaKr (for example alphaKr=0.5 or 1), I couldn't see any changes.
I should say that I know about using more Inter-Particle friction angle( higher than 30 Degree like 50, 60,...) and getting more peak friction angle in the macroscopic response. I've tried this dependency an observed exact results. But I want to see effects of rolling parameters like stiffness and strength in macroscopic parameters.
Thanks Luc Sibille for your suggestion for more studying. I also will study paper [2].

>>[1] Discrete numerical modeling of loose soil with spherical particles
and interparticle rolling friction
RA Hosn, L Sibille, N Benahmed, B Chareyre

>>[2] Quantitative prediction of discrete element models on complex
loading paths
Luc Sibille, Pascal Villard, Félix Darve, Rodaina Aboul Hosn

Revision history for this message
Luc Sibille (luc-sibille) said : 		#5

> But my problem is about the amount of peak friction angle before considering rolling stiffness (I mean without considering CohesionMomment) and after that. when I use just FrictPhys_CundallStrack as the contact law, with the Inter-Particle friction angle of 30 degree, I can't get peak friction angle of experimental data that is about 40 degree.
This is a classical result: with spherical particles and without rolling
resistance, internal peak friction angle relatively high (typically 40
deg) cannot be reproduced.
> Because of this problem, I tried to model a sample with more stiffness until getting more peak friction angle. After using CohesionMoment and setting alpaKr (for example alphaKr=0.5 or 1), I couldn't see any changes.
Again, stiffnesses, as normal, tangential or rolling stiffnesses, if
their are chosen sufficient high should not affect the failure
properties as the peak friction angle. So it may be normal that peak
friction angle is not affected by such stiffnesses. Then peak friction
angle is more strongly affected by the contact friction angle and
plastic rolling moment coefficient (or rolling friction coefficient)
etaRoll (and not alphaKr).

Besides, be sure that the rolling resistance is activated. It is not
because you use Law2_ScGeom6D_CohFrictPhys_CohesionMoment that rolling
resistance is necessarily activated: rolling resistance can be enable or
disable depending on some arguments in
Law2_ScGeom6D_CohFrictPhys_CohesionMoment and/or CohFrictMat

Luc

--
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

Revision history for this message
Elham Hosseinkhani (ehosseinkhani) said : 		#6

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()

Revision history for this message
Luc Sibille (luc-sibille) said : 		#7

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

Revision history for this message
Elham Hosseinkhani (ehosseinkhani) said : 		#8

Hello,
I found the problem of my code and now it works properly. I could reproduce some results in some papers.
I had to change the way of making pack. I used uniform distribution and activate this option; "always_use_moment_law=True", in contact law during making pack and deviatoric loading. I thank from Luc Sibille for his help.