Download this testcase.

Arc-boutement

This example demonstrates the arc-boutement phenomenon in a simple two-dimensional granular configuration. A small chain of circular particles is compressed horizontally against a fixed particle, forming a self-locking structure due to frictional contacts. This effect illustrates how granular materials can support compressive loads through contact force chain.

Keywords

DEM, Rigid Body, Friction

import sys, os, shutil
import numpy as np
from migflow import scontact, time_integration

Output Directory

The output directory is prepared for storing the simulation results.

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)

Particle Problem

We begin by defining the 2D particle problem and configuring the contact solver. The contact geometry is updated dynamically, and the prediction basis is activated to improve convergence of the contact iterations.

p = scontact.ParticleProblem(2)
p.set_fixed_contact_geometry(0)
p.set_predict_basis(1)

Material Properties and Parameters

The particle density, radius, gravity vector and the friction coefficient are defined. A large friction coefficient is assigned between particles to ensure interlocking and prevent sliding.

rho = 1000
r = 0.01
g = np.array([0.0, -9.81])
friction_coefficient = 5
p.set_friction_coefficient(friction_coefficient, "particle", "particle")
p.set_friction_coefficient(0.0, "particle", "wall")

Body Definition

A small chain of five particles is created and initialized as a rigid body. The body is placed slightly above the ground, with each particle in contact with its neighbors. The high inverse inertia value ensures the body behaves as a nearly rigid object.

n = 5  # number of particles in the chain
xp = np.array([0, 3 * r])  # body center position
x_rel = np.arange(n) * 2 * r  # relative particle spacing
x_rel = np.vstack((x_rel, np.zeros_like(x_rel))).T  # particle coordinates
radii = np.full(x_rel.shape[0], r)  # uniform radii
masses = np.full(x_rel.shape[0], np.pi * rho * r**2)  # uniform masses
body = p.init_body(xp, x_rel, masses, radii, material="particle")

Fixed Particle

A fixed particle is added to act as a support against which the body will push. This setup allows the arc-boutement effect to develop between the free body and the fixed particle through frictional contact.

fixed_body = p.add_body(np.array([0.0, r]), 0.0, 0.0)
p.add_particle_to_body(np.array([0.0, 0.0]), r, fixed_body, material="particle")

Boundary Conditions

A horizontal wall is added below the particles to prevent vertical motion.

p.add_boundary_segment(
    np.array([-10 * r, 0.0]), np.array([50 * r, 0.0]), tag="wall", material="wall"
)

Time Integration

The simulation runs for a short time interval with a fixed time step. At each iteration, the gravitational forces are applied to the particles, and the system evolves according to the contact dynamics solver.

t = 0.0
i = 0
dt = 1e-3
tEnd = 250 * dt
outf = 10

while t < tEnd:
    print(f"{i =:4d}, {t:.6g}/{tEnd:.6g}")
    if (i % outf) == 0:
        p.write_mig(outputdir, t)
    forces = np.pi * p.r() ** 2 * rho * g
    time_integration._advance_particles(p, forces, dt, 1, 1e-6 * r)
    t += dt
    i += 1

Plot

python3 -m migflow.plot.migplot output --actors particles