plot the deviatoric part of TriaxialStressController simulation

Asked by Alireza Sadeghi on 2019-04-09

Hello All,

I am using TriaxialStressController to simulate a triaxial test. I want to plot stress in the z-direction vs strain in the z-direction. my code is similar to [1], except without growing of particles. The strain which is plotted in the code has isotropic compaction part too. I was wondering if I could remove an isotropic compaction part and reset the strain and plot just deviatoric part. Thank you very much for your help in advance.

Best Regards

Alireza

[1] http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/examples/triax-tutorial/script-session1.py

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jérôme Duriez
Solved:
2019-04-11
Last query:
2019-04-11
Last reply:
2019-04-11
Jérôme Duriez (jduriez) said : #1

Hi,

plot.resetData() [*] should do the trick.
Obviously you can also export your data in a text file [**] and analyse/plot them as you wish in your favorite data plotting software.

[*] https://yade-dev.gitlab.io/trunk/yade.plot.html?highlight=resetdata#yade.plot.resetData

[**] https://yade-dev.gitlab.io/trunk/yade.plot.html?highlight=resetdata#yade.plot.saveDataTxt

Jérôme Duriez (jduriez) said : #2

PS: the way you browse source code (through http://bazaar.launchpad.net/~yade-pkg/*) is completely outdated, unfortunately. One now has to go to https://gitlab.com/yade-dev/trunk.

Note .py examples scripts are also installed on your computer when you install yade packages (and they also are when you compile from sources, obviously...).

Ask Linux terminal to "locate script-session1.py"

Alireza Sadeghi (asadeghime) said : #3

Hello Jerome,

Thank you for your response and help. I will use https://gitlab.com/yade-dev/trunk. Thank you.
 about my question, I used plot.resetData() but it did not solve my problem. The simulation of triaxial test has two parts. The particles compact with isometric pressure and in the second step we have deviatoric compaction. The problem is, strain in the first step is added to the strain in the second step. but I want to draw strain-Stress diagram just for the second step of triaxial test simulation. my code is in the below of the mail. Thank you very much for your help.

Best Regards

Alireza

#Triaxial compresion for finding a proper size for RVE (stress control boundary condition)

from yade import utils, plot
from yade import pack, qt
from datetime import datetime

qtr=qt.Renderer()
qtr.bgColor=(1,1,1)

#==============================================================
#================= define the materials =======================
#==============================================================

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= True, young=6.81e8, density=1377.5, poisson=0.3, frictionAngle= 0.31, fragile=False, label='Coke'))

O.materials.append(CohFrictMat(normalCohesion= 1e20, shearCohesion= 1e20, isCohesive= False, young=4e9,
density=1523.6, poisson=0.3, frictionAngle= 0.0, fragile=False, label='wall'))

#===============================================================
#=================== define walls ============================
#===============================================================

walls=aabbWalls([(0,0,0),(2e-1,2e-1,2e-1)],thickness=0.0003,oversizeFactor=1.0,material='wall')
wallIds=O.bodies.append(walls)

#===============================================================
#=================== define packing ============================
#===============================================================

nums=['t']

mats=['Coke']

coke=(5e-3,2000)

nums=pack.SpherePack()

nums.makeCloud((0,0,0),(2e-1,2e-1,2e-1),rMean=coke[0],rRelFuzz=1e-4,num=coke[1])

O.bodies.append([utils.sphere(c,r,material=mats[0],color=(0,0,1)) for c,r in nums])

#===============================================================
#=================== define Engine =============================
#===============================================================

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop([Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom6D()],
 [Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
 [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment()]
 ),
# TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
# VTKRecorder(fileName='3d-vtk-',recorders=['all'],iterPeriod=100),
# qt.SnapshotEngine(fileBase='3d-',iterPeriod=200,label='snapshot'),

 NewtonIntegrator(damping=0.4,gravity=[0,0,0])
]
O.dt=2e-6

#=================================================================
#================= APPLYING CONFINING PRESSURE ================
#=================================================================
stabilityThreshold=0.2
sigmaIso=-1e5

triax=TriaxialStressController(

 maxMultiplier=1.000,
 finalMaxMultiplier=1.000,
 thickness = 0,
 stressMask = 7,
 internalCompaction=False,
)

O.engines=O.engines+[triax]

triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=sigmaIso
triax.wall_back_activated=True

while 1:
  O.run(100, True)
  unb=unbalancedForce()
  meanS=(triax.stress(triax.wall_right_id)[0]+triax.stress(triax.wall_top_id)[1]+triax.stress(triax.wall_front_id)[2])/3
  print 'mean Stress:',triax.meanStress,'porosity:', triax.porosity,'meanS:',unb
  if triax.porosity < 0.497 and ((sigmaIso-triax.meanStress)/sigmaIso) < 0.001:
    print 'Isotropic strain1:',-triax.strain[0], 'Isotropic strain 2:',-triax.strain[1], 'Isotropic strain 3:',-triax.strain[2]
    break

print "### Isotropic state saved ###"

#====================================================================
#===================== DEVIATORIC LOADING =======================
#====================================================================
plot.reset()
triax.stressMask = 3
triax.goal1=sigmaIso
triax.goal2=sigmaIso
triax.goal3=-10
#dataOld = plot.data
#plot.saveDataTxt('isotropic_part_deformableRVE.txt')
plot.resetData()
#plot.plot()

O.saveTmp()

#=================================================================
#======================= data collecting =========================
#=================================================================

O.engines=O.engines+[PyRunner(iterPeriod=16000,command='finishIt()',label='calmEngine')]

def finishIt():
# if abs(-triax.strain[2]) > 0.2:
# makeVideo(snapshot.snapshots,'3d.mpeg',fps=10,bps=10000)
    print O.iter
    print 'time is ', O.time
           print 'porosity:',triax.porosity
    print 'strain 1:',-triax.strain[0], 'strain 2:',-triax.strain[1], 'strain 3:',-triax.strain[2]
# O.save('RVE200-deformable.yade')

def history():
   plot.addData(exx=-triax.strain[0], eyy=-triax.strain[1], ezz=-triax.strain[2],
        ev5=(triax.strain[0]+triax.strain[1]+triax.strain[2]),
      sxx=-triax.stress(triax.wall_right_id)[0],
      syy=-triax.stress(triax.wall_top_id)[1],
      szz=-triax.stress(triax.wall_front_id)[2],
                    q=-(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0]),
                    p=triax.meanStress,
                    R=3*(triax.stress(triax.wall_front_id)[2]-triax.stress(triax.wall_right_id)[0])/triax.meanStress,
                    por=triax.porosity,
      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(50000,True)

plot.plots={('ezz'):('ev5')}
plot.plot()

Best Jérôme Duriez (jduriez) said : #4

If your question is about resetting the reference length values for strain estimation (so that strain values start from 0 in the 2nd phase), the solution is to update width0 / height0 / depth0 attributes of TriaxialStressController engine, as in [1] (it's true this is not done in triax-tutorial/script-session1.py).

[1] https://gitlab.com/yade-dev/trunk/blob/master/examples/FluidCouplingPFV/drainage-2PFV-Yuan_and_Chareyre_2017.py#L81 and two following lines

Alireza Sadeghi (asadeghime) said : #5

Thank you dear Jerome, This solved my problem.

Best Regards

Alireza

Alireza Sadeghi (asadeghime) said : #6

Thanks Jérôme Duriez, that solved my question.