Adding blocks not in contact affect the results of previously created blocks

Hello YADE developers,

I'm relatively new to YADE, though I have experience with DEM in general. I'm running a simple test case using polyhedral particles to simulate block behavior on an inclined plane, and I'm encountering unexpected results when introducing a second block.

The setup is straightforward:

• A slope inclined at 20°.
• A single cubical polyhedral block (edge length 1 m) is placed on the slope.
• Initialization is done with high damping (0.8) and a high friction angle (89°), similar to common 3DEC practices.
• The simulation runs until the unbalanced force drops below a threshold.
• After initialization, damping is set to 0.0 and the friction angle is set to 21°, so no motion should occur (as 21° > slope angle).

This works perfectly with one block centered at (1.5, 0, 0.5). However, when I add a second identical block, centered at (3.5, 0, 0.5)—so there is no contact between the two—the behavior changes. Both blocks start sliding significantly, even though they should remain stationary.

I’ve confirmed that the second block is added before any call to O.run(), and all settings are symmetric. This unexpected motion seems to indicate some unintended side effect or internal interaction being triggered.

I'm including the full script below for reproducibility. Any insights would be greatly appreciated.

Thanks in advance for your help!

from yade import plot, polyhedra_utils, export
import numpy as np
import shutil

# Simulation - sliding cubical (polyhedral) block on a 20 degrees inclined slope
        
alpha = radians(20.0)
g = 9.81
mod = 1e11

m = PolyhedraMat(label = 'rockMat')
m.density = 2600  #kg/m^3
kn = 1e9
ks = kn/4.0
m.young = 2*kn  #Pa
m.poisson = 2*ks / m.young
m.frictionAngle = radians(89)  #rad


mat = O.materials.append(m)

vertices0 = [[15, 2, 0], [15, -2, 0], [0, -2, 0], [0, 2, 0], [0, 2, -1], [0, -2, -1], [15, -2, -1], [15, 2, -1]]
bound = polyhedra_utils.polyhedra(O.materials['rockMat'], v = vertices0, fixed=True)
O.bodies.append(bound)

vertices = [[-0.5, 0.5, 0.5], [-0.5, -0.5, 0.5], [-0.5, -0.5, -0.5], [-0.5, 0.5, -0.5], [0.5, 0.5, -0.5], [0.5, -0.5, -0.5], [0.5, -0.5, 0.5], [0.5, 0.5, 0.5]]
slide = polyhedra_utils.polyhedra(O.materials['rockMat'], v = vertices)
slide.state.pos = (1.5, 0., 0.5)
O.bodies.append(slide)

slide2 = polyhedra_utils.polyhedra(O.materials['rockMat'], v = vertices)
slide2.state.pos = (3.5, 0., 0.5)
O.bodies.append(slide2)

for b in O.bodies:
        if b.dynamic:
                print(b.state.pos)
        
dir = 'outputs2/'
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

def myVtk():
        vtkExporter.exportPolyhedra(what=dict(particleVelocity='b.state.vel',n='b.id', force="O.forces.f(b.id)"))
        
def saveToFile(fileName):
        i = 0
        for b in O.bodies:
                if b.dynamic:
                        with open(fileName+'/disp' +str(i) + '.txt', 'a') as file:
                                file.write(str(O.time) + ',' + str(b.state.pos[0]) + ',' + str(b.state.pos[1]) + ',' + str(b.state.pos[2]) + '\n')
                        i = i+1

def contactData():
        for i in O.interactions: print(i.geom.contactPoint)

newton = NewtonIntegrator(damping=0.80, gravity=(g*np.sin(alpha), 0.0, -g*np.cos(alpha)))
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=1000, command = 'contactData()')
]

O.dt = 1e-4

# Initialization step - eliminate excess spring oscillations by running the model with high friction angle and damping
while 1:
        O.run(100,True)
        print(unbalancedForce())
        if unbalancedForce() < 1e-10:
                print("Equilibrium achieved!")
                break

O.resetTime()

# Create output file to track displacement of the block
i = 0
for b in O.bodies:
        if b.dynamic:
                with open(fileName+'/disp' +str(i) + '.txt', 'a') as file:
                        file.write('t,dispX,dispY,dispZ' + '\n')
                        file.write('0' + ',' + str(b.state.pos[0]) + ',' + str(b.state.pos[1]) + ',' + str(b.state.pos[2]) + '\n')
                i = i+1

# Reset damping and friction
newton.damping = 0.0
O.engines += [
        PyRunner(iterPeriod=1000, command = 'myVtk()'),
        PyRunner(iterPeriod=1000, command = 'saveToFile(fileName)')   
]

newPhi = radians(20.5)
tanPhi = tan(newPhi)
O.materials['rockMat'].frictionAngle = newPhi

for i in O.interactions:
        i.phys.tangensOfFrictionAngle = tanPhi

# Run simulation for 10 seconds real time. Since the inclination is 20 degrees and the friction angle is set to 21, no displacment should occur
O.saveTmp()
totalTime = 10
steps = totalTime / O.dt
O.run(int(steps), 1)
utils.waitIfBatch()