yade-users team mailing list archive

[Gary Pekmezi <question295531@xxxxxxxxxxxxxxxxxxxxx>]

New question #295531 on Yade:
https://answers.launchpad.net/yade/+question/295531

This was an issue I ran into before with the Cundall-Strack law, and I am now also seeing with Hertz-Mindlin law.
If I have:

triax=TriaxialStressController(...)

then the "Wallstresses" file shows all stresses and porosity as "NaN" for Hertz-Mindlin, or I get a core dump for Cundall-Strack. This happens only when I set "includeMoment" to True for Hertz-Mindlin, or "setCohesionNow" to True for Cundall-Strack.  It's fairly simple to get my models to run, all I have to do is to switch to :

triax=TriaxialCompressionEngine(...)

So is the moment law and cohesion not supported with TriaxialStressController, is it a bug, or (probably most likely) am I doing something wrong?

I should add that I tested this with the latest version released last week and saw the same behavior

I am pasting my Hertz-Mindlin input below for completeness.

Thanks,
Gary

============================================================================================

import math, os, sys
from yade import plot
from yade import pack
############################################
###                     DEFINE VARIABLES AND MATERIALS               ###
############################################
readParamsFromTable(krot=1.0,noTableOk=True)
from yade.params.table import *
desc='description'
key=O.tags[desc]
print key, krot

os.mkdir(key)

finalFricDegree=14.5
rate=-25.0 # loading rate (strain rate)
damp=0.2 # damping coefficient
stabilityThreshold=0.02
triaxialThreshold=0.005
young=20e9
poisson=0.33
triaxial_p = -5.0e6 

sp=pack.SpherePack()
sp.load('sp.txt')
print "Saved Confined State Loaded"
box=sp.aabb()
print box
print (sp.psd(bins=10))
# #######################################
# ###           DEFINE THE SCENE      ###
# #######################################
O.materials.append(FrictMat(young=young,
                            poisson=poisson,
                            frictionAngle=radians(finalFricDegree),
                            density=2600,
                            label='spheres'))
O.materials.append(FrictMat(young=young,
                            poisson=poisson,
                            frictionAngle=0,
                            density=0,
                            label='walls'))
walls=aabbWalls(box,thickness=0,oversizeFactor=2.0,material='walls')
wallIds=O.bodies.append(walls)
sp.toSimulation(material='spheres')
############################
###   DEFINING ENGINES   ###
############################
triax=TriaxialStressController(                     #replace with#TriaxialCompressionEngine(#
  maxMultiplier=1.005, # spheres growing factor (fast growth)
  finalMaxMultiplier=1.00005, # spheres growing factor (slow growth)
  thickness = -1,
  stressMask = 7,
  internalCompaction=True,
  autoCompressionActivation=False,
  sigmaIsoCompaction=triaxial_p,
  sigmaLateralConfinement=triaxial_p,
  )
triax.goal1=triax.goal2=triax.goal3=confinement_p
          
newton=NewtonIntegrator(damping=damp)
          
vtkrec=VTKRecorder(recorders=['spheres','boxes','stress'],)
vtkrec.fileName=key+'/'+'vtkHC'
          
O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
               [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
	       [Ip2_FrictMat_FrictMat_MindlinPhys(krot=krot,eta=1.0)],
	       [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True,label="mindlinLaw2")],           ),
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        triax,
        newton,
        TriaxialStateRecorder(iterPeriod=100,file=key+'/'+'WallStresses'),
        vtkrec,                                                                                      ]
#####################################################
###    Record and plot data     ###
#####################################################
def history():
        plot.addData(e1=triax.strain[1]*-100.0,
                    ev=-100.0*triax.volumetricStrain,
                    pDiff=pDiff,
                    i=O.iter)
vtkrec.fileName=key+'/'+'vtkTriaxial'
f=open(key+'/'+'./zout.txt', 'w+')
# #######################################
# ###   APPLYING TRIAXIAL PRESSURE   ###
# #######################################
triax.goal1=triax.goal2=triax.goal3=triaxial_p
while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  triDiff=abs((triax.goal1-triax.meanStress)/(triax.goal1))
  print >>f,'unb:',unb,' triDiff: ',triDiff,' meanS: ',triax.meanStress,' porosity: ',triax.porosity,'pMult',triax.previousMultiplier
  sys.stdout.flush()
  if unb<stabilityThreshold and triDiff<triaxialThreshold:
    break
# #######################################
# ###        DEVIATORIC LOADING       ###
# #######################################
#setContactFriction(radians(finalFricDegree))
triax.internalCompaction=False
triax.width0=triax.width  
triax.height0=triax.height
triax.depth0=triax.depth
triax.stressMask = 5
triax.goal2=rate
triax.goal1=triaxial_p
triax.goal3=triaxial_p

while 1:
  O.run(1000,True)
  s11=-triax.stress(triax.wall_right_id)[0]/1000000.0
  s22=-triax.stress(triax.wall_top_id)[1]/1000000.0
  s33=-triax.stress(triax.wall_front_id)[2]/1000000.0
  pDiff=(0.5*((s11-s22)**2+(s11-s33)**2+(s22-s33)**2))**0.5
  history()
  plot.saveDataTxt(key+'/'+'results')
  sys.stdout.flush()
  if triax.strain[1]<(-0.16):
    break

O.save(key+'/'+'triaxialState.yade.gz')
print "###      Triaxial state saved      ###"

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.