Polyhedral particle convergence issue

Summary

Hello,

I am trying to run a simulation for rock slope stability analysis. The setup is a jointed wedge (32 polyhedral particles) resting on a 26.5° inclined slope.

To initialize in-situ stresses (or achieve initial equilibrium), I first set the friction angle to 89° and the damping to 0.8. This follows a procedure similar to what is done in 3DEC: using high friction and damping to suppress excess oscillations so that the unbalanced force should decay toward zero. The full code is pasted below.

However, in my case the unbalanced force does not drop below ~1.0, and the particles slowly creep downslope even with these very high friction and damping values. I tested with yade-batch across a range of time steps and normal stiffness values (see params.txt), but the problem persists.

Geometry setup:

• The slope (fixed) is loaded from slope.json.
• The sliding wedge (32 particles) is loaded from sliding32.json.
• At the start of the simulation, particles are only in touching contact and should not be applying forces on each other.
• To confirm this, I run the simulation for a few hundred steps with gravity set to zero. If any forces appear, the script terminates with an error message.

Main question:

Could this creeping behavior be simply due to the time step not being small enough? With my current parameters (dt=1e-5 to 1e-7), I expected the model to reach equilibrium. Is there something fundamental I am missing in how YADE handles equilibrium for tightly packed polyhedra?

Appreciate your help. Thanks in advance.

Full YADE code:

readParamsFromTable(
        kn = 1e8,
        dt = 1e-5,
        noTableOk=True     # use default values if not run in batch
)
from yade.params import table
kn = table.kn
dt = table.dt
print("Parameters: kn=", kn, " dt=", dt)
from yade import plot, polyhedra_utils, export
import numpy as np
import shutil
import json

# ========= Material Properties =========
m = PolyhedraMat(label = 'rockMat')
m.density = 2600  #kg/m^3
ks = kn/4.0
m.young = 2*kn  #Pa
m.poisson = 2*ks / m.young
m.frictionAngle = radians(89)  # Extremely high friction angle to arrest displacement during the first stage of the simulation

mat = O.materials.append(m)

# ========= Load JSON files with geometry of particles =========
with open("slope.json", 'r') as file:
        boundFile = json.load(file)

for i in range(len(boundFile)):
        vertices0 = boundFile[i]["vertices"]
        bound = polyhedra_utils.polyhedra(O.materials['rockMat'], v = vertices0, fixed=True)
        O.bodies.append(bound)

with open("sliding32.json", 'r') as file:
        slideFile = json.load(file)

for i in range(len(slideFile)):
        vertices0 = slideFile[i]["vertices"]
        slide = polyhedra_utils.polyhedra(O.materials['rockMat'], v = vertices0)
        O.bodies.append(slide)

# ========= Create output directories for vtk visualization for each simulation =========        
dir = f'outputs_kn_{int(kn):.0e}_dt_{(dt):.0e}/'
if os.path.exists(dir):
    shutil.rmtree(dir)
os.makedirs(dir)
vtkExporter = export.VTKExporter(dir)
vtkExporter.exportPolyhedra(what=dict(particleVelocity='b.state.vel',n='b.id', force="O.forces.f(b.id)"))

fileName = dir

# ========= Helper function to calculate average displacement of the sliding mass =========        
def calcAvgPos():
        i = 0
        pos = np.zeros(3)
        for b in O.bodies:
                if b.dynamic:
                        pos += b.state.pos
                        i = i+1
        return pos/float(i)

# ========= VTK output function =========        
def myVtk():
        vtkExporter.exportPolyhedra(what=dict(particleVelocity='b.state.vel',n='b.id', force="O.forces.f(b.id)"))

g = 9.81
O.dt = dt
vtkOut = math.floor(0.1/O.dt)
damp = 0.8 # note the high damping to damp out excess ocsillations in the system
newton = NewtonIntegrator(damping=0.8, gravity=(0.0, 0.0, -g))
O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Polyhedra_Aabb()]),
        InteractionLoop(
                [Ig2_Polyhedra_Polyhedra_PolyhedraGeom()],
                [Ip2_PolyhedraMat_PolyhedraMat_PolyhedraPhys(label = 'myPhys')],
                [Law2_PolyhedraGeom_PolyhedraPhys_Volumetric()]
        ),
        newton,
        PyRunner(iterPeriod=vtkOut, command = 'myVtk()')
]

# ========= Run the model with 0 gravity to detect contacts at the beginning of the simulation =========
newton.gravity = (0.0,0.0,0.0)
for step in range(200):
    O.step()        
    for i in O.interactions:
        b1, b2 = i.id1, i.id2
        f1 = O.forces.f(b1)
        f2 = O.forces.f(b2)
        print(f"Contact {b1}-{b2}: Force on {b1}={f1}, Force on {b2}={f2}")
        if i.isReal:
            print("Particles are in contact when gravity is not applied. Terminate.")
            quit()
        
# ========= Calculate initial average position of the sliding mass =========
pos0 = calcAvgPos()
outFile = os.path.join(dir, "unbal.txt") # store time and unbalanced force to check for convergence
newton.gravity = (0.0,0.0,-g)

# ========= Run the model with high damoing and friction angle so the unbalanced force decays =========
with open(outFile, 'w') as file:
    while True:
        O.run(100, True)
        ub = unbalancedForce()
        print(ub, flush=True)
        file.write(f"{O.time},{ub}\n")
        pos = calcAvgPos() # calculate current average position of the sliding mass
        if np.linalg.norm(pos - pos0) > 0.5: # termination criterion
            print("Equilibrium NOT achieved")
            quit()
        if ub < 1e-5:
            print("Equilibrium achieved")
            O.save(os.path.join(dir, 'sliding_eq.xml'))
            break

utils.waitIfBatch()

Code:params.txt

sliding32.json

slope.json