Triaxial compression test displacement field

Asked by maquankun on 2018-08-29

I tried to make the displacement field of the second stage of triaxial compression. I have this answer in the following question, but I have not implemented this method. Here is my code, please tell me where the error is.
https://answers.launchpad.net/yade/+question/478212.
from yade import pack,plot,qt,export
nRead=readParamsFromTable(
        num_spheres=1000,
        compFricDegree = 30,
        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 = 30
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=100000
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.]
                    #[0.01,0.03,0.11,0.13,0.2],[0.,0.1,0.5,0.6,1.]
                    # [0.05,0.07,0.14,0.155,0.2],[0.,0.1,0.5,0.6,1.]
                    #[0.1,0.11,0.16,0.17,0.2],[0.,0.1,0.5,0.6,1.]
                    #[0.15,0.16,0.184,0.19,0.2],[0.,0.1,0.5,0.6,1.]
                    #[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()
vtk = export.VTKExporter('/home/ma/dspl/jan-')
vtk.exportSpheres(useRef=True,what=[('dspl','b.state.displ()')])

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Bruno Chareyre
Solved:
2018-08-31
Last query:
2018-08-31
Last reply:
2018-08-29
Jan Stránský (honzik) said : #1

> please tell me where the error is.

please tell us what the error is.

Cheers
Jan

maquankun (maquankun) said : #2

Thank you for your attention!
When I start running the program, I can only generate a vtk file "-sphere-0000.vtk" at the beginning of compression, and I can't continue to generate files as the program runs.

Jan Stránský (honzik) said : #3

Thanks for clarification :-)

as you have it, the vtk.exportSpheres is called just once, hence you get only one file..
To get vtk files regularly during the simulation, you have to use PyRunner [1] in O.engines and call an export function from it, something like:

####
def doExport():
   vtk.exportSpheres(...)

O.engines = [
   ...
   PyRunner(iterPeriod=1000,command="doExport()"), # you can use realPeriod or virtPeriod, see [1]
]
####

or you can do the run "in cycle":
####
for i in range(20):
   O.run(1000,True)
   vtk.exportSpheres(...)
####

cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.PyRunner

Hi,
Your script executes "vtk.exportSpheres(useRef=True,what=[('dspl','b.state.displ()')])" once, and as a consequence you have one vtk snapshot saved. There is no "error" here.

If you want multiple snapshots you can:

for k in range(N): #or any other control loop
    O.run(M,1)
    vtk.exportSpheres(useRef=True,what=[('dspl','b.state.displ()')])

or insert a PyRunner in O.engines and run as usual:
O.engines=O.engines+[PyRunner(interval=M,command="vtk.exportSpheres(useRef=True,what=[('dspl','b.state.displ()')])")]

Bruno

maquankun (maquankun) said : #5

Sorry, I replied so late, this solved my problem.

maquankun (maquankun) said : #6

Thanks Bruno Chareyre, that solved my question.