← Back to team overview

yade-users team mailing list archive

[Question #700411]: Aborted (core dumped) occured when I use the Hertz-Minlin contact law with non-zero krot and eta

 

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

python3.8: /builds/yade-dev/trunk/deb/yadedaily/pkg/common/InsertionSortCollider.cpp:778: bool yade::InsertionSortCollider::spatialOverlapPeri(yade::Body::id_t, yade::Body::id_t, yade::Scene*, yade::Vector3i&) const: Assertion `maxima[3 * id1 + axis] - minima[3 * id1 + axis] < .99 * dim' failed.
python3.8: /builds/yade-dev/trunk/deb/yadedaily/pkg/common/InsertionSortCollider.cpp:778: bool yade::InsertionSortCollider::spatialOverlapPeri(yade::Body::id_t, yade::Body::id_t, yade::Scene*, yade::Vector3i&) const: Assertion `maxima[3 * id1 + axis] - minima[3 * id1 + axis] < .99 * dim' failed.
Aborted (core dumped)
I am running a periodic triaxial test using hertz-mindlin contact law. However, the program abort suddenly as soon as I add the rolling stiffness coefficient and rolling plastic limit (as shown above). If these values are set to 0, the simulation can continue.
My code is as follows:
# encoding: utf-8

readParamsFromTable(
num=2500,
)

import glob, os
import numpy as np
from yade.params import table
from yade import pack, plot

isBatch=runningInBatch() #check if run in batch mode
if isBatch:
    print('Running: '+O.tags['description']) #mark the current simulation

num = table.num # number of soil particles
dScaling = 1e3 # density scaling to increase crticial time step
sigmac= -100000 # confining pressure
rate = 0.1 # strain rate
damp = 0.5 # damping coefficient
stabilityRatio = 1.e-3 # threshold for quasi-static condition
stressTolRatio = 1.e-3 # tolerance for stress goal

# corners to define specimen size
m_1=0.004*pow(num/2500,1/3)
m_2=0.008*pow(num/2500,1/3)
mn,mx=Vector3.Zero,Vector3(m_1,m_1,m_2)

# Soil sphere parameters
E=7e9 # micro Young's modulus
v=0.3 # micro Poisson's ratio
kr=0.5 # rolling/bending stiffness coefficient
eta=0.5 # rolling/bending plastic limit
mu = 26 # contact friction angle
rho = 2650*dScaling # soil density

# create materials
spMat = O.materials.append(
   FrictMat(young=E,poisson=v,frictionAngle=radians(mu),density=rho))

# create a cloud of ramdomly positioned spheres
O.periodic = True
O.cell.hSize = Matrix3(mx[0],0,0, 0,mx[1],0, 0,0,mx[2])
sp=pack.SpherePack()
sizes=[.00015,.000225,.0003,.0003625,.00045,.0006,0.00118]
cumm=[0,.215,.3807312826,.601,.8865351132,.9965177017,1]
sp.makeCloud(mn,mx,psdSizes=sizes,psdCumm=cumm,distributeMass=True,seed=1,num=num)

# load spheres to simulation
spIds=sp.toSimulation(material=spMat)

#define engines
O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys(
            krot=kr,
            ktwist=kr,
            eta=eta
        )],
   [Law2_ScGeom_MindlinPhys_Mindlin(includeMoment=True)]
   ),
   PeriTriaxController(label='triax',
      goal=(sigmac,sigmac,sigmac),
      # whether they are strains or stresses
      stressMask=7,
      # type of servo-control
      dynCell=True,
      maxStrainRate=(10*rate,10*rate,10*rate),
      # wait until the unbalanced force goes below this value
      maxUnbalanced=stabilityRatio,
      relStressTol=stressTolRatio,
      doneHook='compactionFinished()'
   ),
   NewtonIntegrator(damping=damp,label='newton'),
   PyRunner(iterPeriod=500,command='addPlotData()',label='data')
]

O.dt=.3*PWaveTimeStep()

def compactionFinished():
    if unbalancedForce()<stabilityRatio:
        print(porosity())
        triax.doneHook='startShearing()'


def startShearing():
    # set the current cell configuration to be the reference one
    O.cell.trsf = Matrix3.Identity
    O.cell.velGrad = Matrix3.Zero
    setRefSe3()
    # set loading type: constant pressure in x,y, 30% compression in z
    triax.stressMask=3
    triax.globUpdate=1
    # allow faster deformation along x,y to better maintain stresses
    triax.maxStrainRate=(10*rate,10*rate,rate)
    triax.goal=(sigmac,sigmac,-.3)
    triax.maxUnbalanced=stabilityRatio
    triax.relStressTol=stressTolRatio
    triax.doneHook='Finished()'

def Finished():
    print(str(num)+'finished')
    O.pause()

def addPlotData():
    if triax.stressMask != 7:
        plot.addData(i=O.iter,t=O.realtime,s11=-triax.stress[0],\
                    s22=-triax.stress[1],s33=-triax.stress[2],e11=triax.strain[0],\
                    e22=triax.strain[1],e33=triax.strain[2],n=porosity(),\
                    V=O.cell.volume,M=triax.mass,Unb=unbalancedForece())
        if isBatch:
            dataName = 'triax2_'+ O.tags['description'] + '.txt'
        else:
            dataName = 'TriaxREV2_'+str(num)+'.txt'
        plot.saveDataTxt(dataName,vars=('s11','s22','s33','e11','e22','e33','n','V','M'))

O.run()
waitIfBatch()


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