Derived strain field in triaxial compression test

Asked by maquankun

I am trying to export the Micro-strain legend on the tutorial after performing a triaxial compression test. https://yade-dem.org/doc/user.html#micro-stress-and-micro-strain. But I can't open it with paraview using the generated strain.vtk file. If you can export the changes in the strain field over time with compression, it will be more suitable for me.Thanks for any help.
from yade import pack,plot,qt,export
nRead=readParamsFromTable(
        num_spheres=1000,
        compFricDegree = 50,
        key='_triax_base_',
        unknownOk=True
)
from yade.params import table
num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.42
compFricDegree = table.compFricDegree
finalFricDegree = 50
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=5e8
mn,mx=Vector3(0,0,0),Vector3(1,2,1) # corners of the initial packing
weiya=400000
poisson=0.1
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(compFricDegree),density=2650,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)
psdSizes,psdCumm =[0.19,0.191,0.195,0.196,0.2],[0.,0.1,0.5,0.6,1.]

sp=pack.SpherePack();
sp.makeCloud(mn,mx,num=num_spheres, psdSizes=psdSizes, psdCumm=psdCumm,distributeMass=True,porosity=0.42)
O.bodies.append([sphere(center, rad, material='spheres') for center, rad in sp])
triax=TriaxialStressController(
        maxMultiplier=1.+2e4/young,
        finalMaxMultiplier=1.+2e3/young,
        thickness = 0,
        stressMask = 7,
        internalCompaction=False,
)
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()],
                [Law2_ScGeom_FrictPhys_CundallStrack()]
        ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax,
        TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
        newton,
]
yade.qt.Controller(), yade.qt.View()
triax.goal1=triax.goal2=triax.goal3=-weiya
while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
  if unb<stabilityThreshold and abs(-weiya-triax.meanStress)/weiya<0.001:
    break
print "### Isotropic state saved ###"
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)
print 'Box Volume: ', triax.boxVolume
print "### Compacted state saved ###"
triax.internalCompaction=False
setContactFriction(radians(finalFricDegree))
triax.stressMask = 5
triax.goal1=-weiya
triax.goal2=rate
triax.goal3=-weiya
newton.damping=0.1
O.saveTmp()
TW=TesselationWrapper()
TW.computeVolumes()
TW.volume(20)
TW.setState(0)
O.run(1000,True)
TW.setState(0)
TW.defToVtk("strain.vtk")

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
Last query:
Last reply:
Revision history for this message
Bruno Chareyre (bruno-chareyre) said :
#1

Hi, the documentation under your link reads:
TW.setState(0)
O.run(1000,True)
TW.setState(1)

Check your script again? :)

Bruno

Revision history for this message
Best Robert Caulk (rcaulk) said :
#2

TW.setState() is used to set states for comparison in order to generate the strain field. If you want to TW to compute a strain field, set the second state correctly (according to the link you provided) using:

TW.setState(1).

>>If you can export the changes in the strain field over time with compression, it will be more suitable for me.Thanks for any help.

You're welcome:

TW=TesselationWrapper()
TW.computeVolumes()
TW.volume(20)
TW.setState(0)
step=0
while 1:
    step +=1
    O.run(1000,True)
    TW.setState(1)
    TW.defToVtk("strain_"+str(step)+".vtk")

Revision history for this message
maquankun (maquankun) said :
#3

Thanks Robert Caulk, that solved my question.