Download this testcase.

Isotropic Compression of a Dry Granular Material

This example simulates the isotropic compression of a granular assembly confined between four moving rigid walls. The walls apply an increasing force over time, compressing the particles uniformly.

Keywords

DEM, Compression, Granular Material

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

np.random.seed(0)

Output Directory

Create a clean output directory for simulation results.

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

Particle Problem

A 2D DEM problem is created.

p = scontact.ParticleProblem(2)

Material Properties

In this section, we define the physical properties of the granular material.

  • rhop : particle density (kg/m³)
  • r : reference particle radius

To introduce polydispersity and obtain a more realistic packing, particle radii are sampled uniformly between 0.5·r and 1.0·r.

np_max limits the maximum number of particles that will be generated.

rhop = 2600
r = 1e-4
np_max = int(1e6)
radius = np.random.uniform(0.5, 1.0, np_max) * r

Geometry Definition

The granular domain is rectangular and defined by:
  • width = w
  • height = h

It is enclosed by four rigid walls (top, bottom, left, right), which will later apply compressive loading to the particle assembly.

Domain dimensions scale with the particle radius r. Wall sizes are set larger than the domain to ensure full coverage of the granular domain.

h = 40 * 2 * r
w = 40 * 2 * r

vertical_wall_size = 1.5 * h
horizontal_wall_size = 1.5 * h

Particles Initial Conditions

Particle centers are initialized on a regular 2D grid with spacing 2r, ensuring no initial overlap between grains.

The grid spans the rectangular granular domain of width w and height h. For each grid position, a particle is created using its assigned radius and computed mass (mass = ρ * π * r²).

x = np.arange(r, w - r, 2 * r)
y = np.arange(r, h - r, 2 * r)
xp = np.array(np.meshgrid(x, y)).T.reshape(-1, 2)
for xi, ri in zip(xp, radius):
    p.add_particle(xi, ri, np.pi * ri**2 * rhop, "Sand")

Walls Initial Conditions

Four rigid walls (top, bottom, left, right) are created around the domain. Each wall is modeled as a rigid body with a sufficiently large mass to ensure that the particle–wall interactions behave correctly during compression.

Wall mass is chosen proportional to the total particle mass:
mass_wall ≈ Σ (ρ * π * r_i²)

Walls are represented by line segments attached to these rigid bodies. Rotational inertia is disabled (invertI = 0) to prevent wall rotation. A small offset eps shifts the walls slightly outside the particle domain.

wall_mass = np.sum(p.r() ** 2 * np.pi * rhop)

eps = 1 * r

x0 = np.array([-vertical_wall_size / 2, 0])
x1 = np.array([vertical_wall_size / 2, 0])

top_body = p.add_body(np.array([w / 2, h + eps]), invertI=0, invertM=1 / wall_mass)
p.add_segment_to_body(x0, x1, top_body, "Top", "wall")

bottom_body = p.add_body(np.array([w / 2, -eps]), invertI=0, invertM=1 / wall_mass)
p.add_segment_to_body(x0, x1, bottom_body, "Bottom", "wall")

y0 = np.array([0, horizontal_wall_size / 2])
y1 = np.array([0, -horizontal_wall_size / 2])

left_body = p.add_body(np.array([-eps, h / 2]), invertI=0, invertM=1 / wall_mass)
p.add_segment_to_body(y0, y1, left_body, "Left", "wall")

right_body = p.add_body(np.array([w + eps, h / 2]), invertI=0, invertM=1 / wall_mass)
p.add_segment_to_body(y0, y1, right_body, "Right", "wall")

Simulation Parameters

This section defines the mechanical interaction properties and initializes the force fields acting on particles and walls.

Friction coefficients:
  • Particle-Particle (Sand-Sand)
  • Particle-Wall (Sand-wall)
External forces:
  • fext : external forces applied to each particle (none initially)
  • fbody : external forces applied to rigid bodies (walls) at t = 0
p.set_friction_coefficient(0.4, "Sand", "Sand")
p.set_friction_coefficient(0.4, "Sand", "wall")

fext = np.zeros((p.n_particles(), 2))
fbody = np.zeros((p.n_bodies(), 2))

Time Integration of the Simulation

The simulation is advanced in time using a fixed time step. Numerical parameters are defined below, and the initial state of the system is written to the output directory.

Parameters:
  • outf : number of iterations between output files
  • dt : time step for numerical integration
  • tEnd : final simulation time
  • t : initial simulation time
  • ii : iteration counter

Computational Loop: At each time step, the following operations are performed:

  1. A ramp function gradually increases the compressive forces applied to the four rigid walls at the start of the simulation to prevent sudden shocks. Once the simulation reaches the midpoint, the force reaches its target value and remains constant for the remainder of the run.

  2. Forces are applied symmetrically to the top, bottom, left, and right walls.

  3. Particle positions and velocities are advanced using the DEM solver via time_integration._advance_particles

    • Updates particle positions from velocities
    • Computes contact forces between particles and walls
    • Resolves collisions using sub-iterations for stability
    • Incorporates external and wall-applied forces
    • Sets a maximum allowed particle overlap (contact tolerance) to ensure numerical stability

  4. Mean wall velocity is subtracted to the walls to prevent numerical drift.

  5. Simulation state is periodically written to the output directory for visualization and post-processing.

outf = 10
dt = 5e-5
tEnd = 0.1
t = 0
i = 0


def ramp(t, tStop, tStart=0, value=10):
    if t < tStart:
        return 0
    elif t > tStop:
        return value
    else:
        return (t - tStart) / (tStop - tStart) * value


while t < tEnd:
    print(f"{i:4d}, {t:.6g}/{tEnd:.6g}, {dt:.6g}")
    if i % outf == 0:
        p.write_mig(outputdir, t)

    body_force = wall_mass * ramp(t, tEnd / 2, 0, 20)
    fbody[top_body, 1] = -body_force
    fbody[bottom_body, 1] = +body_force
    fbody[left_body, 0] = +body_force
    fbody[right_body, 0] = -body_force

    time_integration._advance_particles(
        p, fext, dt, min_nsub=1, contact_tol=1e-4 * 2 * r, fbody=fbody
    )

    # Remove numerical drift: set mean wall velocity to zero
    mean_velocity = np.mean(p.body_velocity(), axis=0)
    p.body_velocity()[:] -= mean_velocity

    t += dt
    i += 1

Plot

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