Research Axis

NSCD

The NSCD axis develops MigFlow's grain-scale solver: the Non-Smooth Contact Dynamics (NSCD) method used to resolve collisions between rigid particles. NSCD is a "hard-sphere" method — it treats contact as an instantaneous, non-overlapping event and solves directly for the impulse that enforces non-penetration, rather than integrating a contact force over many small time steps. This buys a larger, more robust time step and strict volume conservation of the granular packing, at the cost of solving a genuinely non-smooth (non-differentiable) problem at every contact, and — when many grains touch at once, as in a dense packing — a coupled multi-contact problem across the whole assembly. This axis underpins every other one at MigFlow: it is the solid-phase engine driving immersed granular flows, non-isothermal grain beds, and any coupling with the mesh-based fluid solvers.

Key Concepts & Methodology

Contact Kinematics & the Delassus Operator

Each potential contact defines its own local frame, built from the contact normal $\mathbf{n}$ and two tangent directions, in which the relative velocity between the two bodies and the contact impulse are expressed. Projected into this frame, the effect of an impulse on the relative velocity at the contact point is linear, governed by the Delassus operator $W$ — essentially the combined inverse-mass response of both bodies, seen at the contact point:

$$ \delta\mathbf{u}^+ = \delta\mathbf{u}^- + W\mathbf{r}, \qquad W = W_A + W_B, \qquad W_X = H_X^{\mathsf T} M_X^{-1} H_X, $$

where $\delta\mathbf{u} = \mathbf{u}_A - \mathbf{u}_B$ is the relative velocity at the contact point, $\mathbf{r}$ the contact impulse, $M_X$ the mass matrix of body $X$, and $H_X$ the transformation from global degrees of freedom to the local contact frame. This single relation is what turns "find the forces that prevent overlap" into a solvable local problem at every contact.

Signorini's Condition — Non-Penetration

NSCD is a "hard-sphere" method: grains are never allowed to overlap, so the gap $g(t)$ between two bodies and the normal contact impulse $r_n$ must satisfy a complementarity condition rather than a spring-like force law. Either the bodies are apart and exert nothing on each other, or they are exactly touching and the impulse can be positive — never both partially overlapping and pushing, and never attractive:

$$ g(t) \geq 0, \qquad r_n(t) \geq 0, \qquad r_n(t)\, g(t) = 0. $$

Combined with an inelastic impact rule — the normal relative velocity vanishes as soon as contact occurs, $g(t)=0 \Rightarrow \delta u_n^+ = 0$ — this is what gives the granular packing its strict incompressibility.

Coulomb's Friction Law

Tangentially, the contact impulse is confined to a friction cone set by the coefficient $\mu$: while the contact sticks, the tangential impulse can take any value inside the cone; as soon as the two surfaces slip against each other, the tangential impulse sits exactly on the cone's boundary, opposing the slip direction:

$$ \|\mathbf{r}_t\| \leq \mu r_n, \qquad \|\delta\mathbf{u}_t^+\| \neq 0 \;\Rightarrow\; \mathbf{r}_t = -\mu r_n \frac{\delta\mathbf{u}_t^+}{\|\delta\mathbf{u}_t^+\|}. $$

In two dimensions, whether a contact sticks or slips — and in which direction — can be read off directly by comparing the "stick impulse" (the impulse that would cancel the relative velocity outright) against the friction cone; in three dimensions, the same decision generally requires an iterative search.

Resolving Many Contacts at Once

A single contact is solvable in closed form, but a dense granular packing has many contacts active within the same time step, and solving one changes the velocities that all its neighbours see. MigFlow resolves this coupled problem with a Non-Linear Gauss–Seidel sweep: contacts are visited one at a time, each re-solved against the bodies' current velocities, until every contact's relative-velocity update falls below a tolerance $\xi$,

$$ \big\| \delta\mathbf{u}_k^{\alpha} - \delta\mathbf{u}_{k-1}^{\alpha} \big\|\, \Delta t < \xi. $$

Rather than sweeping the full contact list every iteration, MigFlow keeps a queue: only contacts that have not yet converged — and their immediate neighbours, most likely to be disturbed next — are re-queued, so computational effort concentrates where the packing hasn't settled yet.

Featured Simulation

Cylindrical Column of Elongated Grains

Supplementary documentation
Simulated by

Simulation Overview

A cylindrical 3D column is packed with elongated grains — each one a rigid rod made of 8 spherical sub-grains bonded together — and left to settle and interact under NSCD's usual rules — contacts are resolved as instantaneous, non-overlapping impulses, exactly as described in the axis's Key Concepts above, only now applied to elongated bodies instead of spheres. The elongated shape lets grains interlock and align with one another in ways single spheres cannot, which changes how the packing jams and how forces are transmitted through it.

This 3D exploration follows on from the 2D column-collapse study of Coppin et al. (2023), which used shorter line-clusters of 3–4 disks.

Featured Simulation

Grain Contact Network

Supplementary documentation
Simulated by

Simulation Overview

This simulation studies the drag force exerted on a fixed circular obstacle immersed in a granular flow saturated by a viscous fluid — a 2D DEM–FEM coupling test in which grains are continuously drained at the bottom of the domain and re-injected at the top, streaming past the disk under an imposed inflow. As they do, the contact network between grains is drawn directly: the impulses NSCD computes at every contact — the forces that carry stress through the packing from grain to grain — are made visible as they evolve, rather than staying hidden inside the solver. The total drag force on the obstacle is recorded throughout, with a 2D–3D scaling factor applied to the particle volumes to correct for the reduced dimensionality of the 2D geometry.

Key Parameters

Geometry
Particle Bed Size $0.3 \times 0.3\ \mathrm{m}$
Domain Height $h = 1.4\ \mathrm{m}$
Obstacle Radius $R = 0.025\ \mathrm{m}$
2D–3D Correction Ratio $0.7$
Physical Parameters
Imposed Stream Velocity $-0.05\ \mathrm{m/s}$
Particle Radius $r = 3\ \mathrm{mm}$
Particle / Fluid Density $2500 / 1000\ \mathrm{kg/m^3}$
Fluid Kinematic Viscosity $10^{-6}\ \mathrm{m^2/s}$
Friction Coefficient $\mu = 0.2$
Gravity $g = 9.81\ \mathrm{m/s^2}$
Numerical Parameters
Time Step $dt = 5\times10^{-3}\ \mathrm{s}$
Final Time $15\ \mathrm{s}$
Featured Simulation

Mixing of Cross-Shaped Grain Clusters

Supplementary documentation
Simulated by

Simulation Overview

Here, grains are bonded together into rigid, cross-shaped clusters in 3D — rather than loose individual grains — and coloured to track where they end up over time. Mixing them in a rotating mixer this way makes it possible to directly observe, cluster by cluster, whether the mixer actually blends the different coloured compounds together or leaves them segregated.

Featured Simulation

Density Sorting by Water Jigging

Simulated by

Simulation Overview

Water jigging sorts a binary mixture of grains by density: a pulsating water flow repeatedly resuspends and re-settles the pile, and over many cycles the denser grains work their way to the bottom while the lighter ones rise to the top. This case reproduces a laboratory-scale jig with the same unresolved VANS–NSCD coupling used across this axis — grains solved individually with NSCD, the interstitial fluid on a coarser mesh — with the cyclic injection–drainage of the jig imposed as a pulsating inflow velocity at the bottom of the tank,

$$ \frac{\partial \phi}{\partial t} + \nabla\cdot\mathbf{u} = 0, \qquad \frac{\partial \rho\mathbf{u}}{\partial t} + \nabla\cdot\!\left(\frac{\rho\mathbf{u}\mathbf{u}}{\phi}\right) = \nabla\cdot\big(2\mu\phi\,\mathbf{d}\big) - \phi\nabla p + \mathbf{f} + \phi\rho\mathbf{g}, $$

with $\phi$ here the fluid volume fraction, and the bottom boundary driven by $u_b = \frac{h\pi}{T}\sin\!\big(\frac{2\pi}{T}t\big)$ — a sine matched to the water elevation $h$ and period $T$ measured on the physical jig. Each grain feels this fluid back through a pressure-gradient term and an empirical, porosity-corrected drag force, while contacts between grains are resolved by NSCD with a perfectly inelastic law (no restitution). Density sorting emerges directly from these ingredients: because the drag and buoyancy felt by a grain depend on its density, denser and lighter grains respond differently to the same pulsating flow, and repeated cycles of fluidisation and settling separate them over time.

Key Parameters

Materials
Ceramic $\rho = 3670\ \mathrm{kg/m^3}$
Glass $\rho = 2500\ \mathrm{kg/m^3}$
Plastic (POM) $\rho = 1410\ \mathrm{kg/m^3}$
Density Ratio (glass/plastic) $1.77$
Density Ratio (ceramic/glass) $1.47$
Grains & Domain
Grain Diameter $d_s \in [9, 11]\ \mathrm{mm}$
Number of Grains $\approx 22\,000$
Tank Size $33\times33\times35\ \mathrm{cm}$
Jigging Cycle
Water Elevation $h = 7.0 \pm 0.1\ \mathrm{cm}$
Period $T = 0.776 \pm 0.032\ \mathrm{s}$
Particle Reynolds Number $Re_s \approx 3113$
Numerical
Wall Penalisation Coefficient $C = 600$ no-slip enforcement at lateral walls
Friction Coefficients Explored $\mu \in [0.1, 0.9]$ grain–wall friction found most influential
Real-world Phenomenon

Jigging: Sorting by Gravity Concentration

Jigging is a mineral-processing technique that separates grains by density using nothing but gravity and a pulsating fluid: a granular bed is repeatedly lifted and allowed to settle by cycles of water injection and drainage, and over many cycles the denser material works its way to the bottom of the bed while the lighter material rises to the top.

It is used industrially to concentrate ores and, increasingly, to sort construction and demolition waste into reusable aggregate fractions. The laboratory reference for this simulation is a benchtop water jig (Allmineral Alljig), where binary mixtures of ceramic, glass, and plastic beads were sorted and the resulting layers sampled to track the separation over successive cycles.

Featured Simulation

Granular Flow in a Rotating Mill

Simulated by

Simulation Overview

This case simulates a rotating mill with granular particles inside: a bed of grains sits inside a circular drum, and the drum wall is driven at a prescribed angular velocity. Particles interact with the mill boundaries and among themselves through friction, and NSCD is used here in its purely dry, unresolved-fluid form — no VANS coupling — so the whole behaviour comes from grain–grain and grain–wall contacts under gravity and rotation.

The mill wall itself is not physically rotated: instead, the wall segments are kinematically assigned the velocity field of a rigid rotation about the drum centre at every step, and reset to their reference position before each contact solve. The angular velocity is set from a target Froude number $Fr = \omega^2 R/g$, the standard parameter used to classify the transverse motion regime of a rotating drum (slipping, slumping, rolling, cascading, cataracting, or centrifuging) — see Mellmann (2001) in the Supplementary documentation above.

Key Parameters

Mill & Grains
Mill Radius $R = 0.5\ \mathrm{m}$
Grain Radius $r = 10\ \mathrm{mm}$
Grain Density $\rho = 2500\ \mathrm{kg/m^3}$
Initial Fill Region $1.0\times0.6\ \mathrm{m}$
Rotation & Friction
Froude Number $Fr = 2.5$
Angular Velocity $\omega \approx 3.50\ \mathrm{rad/s}$
Grain–Grain Friction $\mu = 0.8$
Grain–Wall Friction $\mu = 0.8$
Numerical
Time Step $\Delta t = 10^{-3}\ \mathrm{s}$
Final Time $10\ \mathrm{s}$
Featured Simulation

Heat Conduction in a Discharging Silo

Supplementary documentation
Simulated by

Simulation Overview

This example studies heat conduction through a granular material flowing out of a silo, with thermal interaction between particles. A dense bed of polydisperse grains is packed inside a hopper-shaped silo — a rectangular body with a narrow converging outlet — and released to discharge under gravity, with NSCD resolving every grain–grain and grain–wall contact as the packing empties out.

Heat is carried purely through this same contact network: whenever two grains touch, they exchange heat at a rate set by a Hertzian-type contact conductance built from each grain's thermal conductivity, Young's modulus and Poisson's ratio, so a wider or stiffer contact conducts more heat, exactly as it transmits more force. The silo walls are held at a fixed, elevated temperature and act as a heat source without ever storing heat themselves, continuously injecting warmth into the pile at every wall–grain contact; grains start uniformly cold and heat up only as they touch a warm neighbour or the wall. Particles that leave through the outlet are removed from the simulation together with their temperature, so the discharging flow carries the emerging temperature field out of the domain along with the grains themselves.

Key Parameters

Silo Geometry
Silo Height / Width $h = 10\ \mathrm{m}$, $w = 20\ \mathrm{m}$
Outlet Height / Width $h_{\mathrm{out}} = 4\ \mathrm{m}$, $l_{\mathrm{out}} = 2\ \mathrm{m}$
Grain Radius $r \in [0.7, 1]\times0.1\ \mathrm{m}$ polydisperse, deposited by fast-deposit
Number of Grains $5000$
Physical & Thermal Contact Model
Density / Heat Capacity $\rho = 1\ \mathrm{kg/m^3}$, $c_p = 1\ \mathrm{J/(kg\cdot K)}$ normalized, not a calibrated material
Gravity $g = 9.81\ \mathrm{m/s^2}$
Contact Thermal Conductivity $k_s = 1\ \mathrm{W/(m\cdot K)}$ normalized
Young's Modulus / Poisson's Ratio $E = 1\ \mathrm{Pa}$, $\nu = 0.3$ normalized
Thermal Boundary & Initial Conditions
Wall Temperature $T_{\mathrm{wall}} = 1\ \mathrm{K}$ fixed heat source, no self-storage; normalized
Initial Grain Temperature $T_0 = 0\ \mathrm{K}$
Numerical
Time Step $\Delta t = 2.5\times10^{-3}\ \mathrm{s}$
Final Time $5\ \mathrm{s}$
Heat Sub-Iterations $5$ per NSCD step