Immersed Granular Flows

Simulation Overview

A compact, roughly spherical swarm of grains — a "Stokes cloud" — is released from rest in a viscous fluid and left to fall under gravity alone in a fully three-dimensional domain. This is the unresolved CFD-DEM model at its most classical: grains are solved individually with NSCD contacts, the fluid is solved on a coarser mesh with a stabilised P1-P1 finite element scheme (a pressure-Laplacian term keeps the equal-order velocity–pressure pair well-posed), and the two are coupled through a Dallavalle-type drag law corrected by a Wen & Yu porosity function, evaluated with the same semi-implicit scheme used elsewhere on this axis to keep the momentum exchange stable. The fluid mesh is periodically rebuilt and refined tightly around the grains' current position as the cloud moves, so resolution follows the swarm rather than being fixed in space.

As the cloud falls, it drives its own circulating flow: fluid streamlines close around it, elongating the upper part while the lower part flattens. Grains at the rear continually leak out of this closed envelope into a growing trailing tail — the leakage rate set by the cloud's falling speed, its radius, and the grain radius — until the depleted core lets fluid punch straight through, turning the swarm into an open torus that goes on to break into one or more secondary droplets. This configuration releases two clouds with a diagonal offset, reproducing the interacting-cloud case of Constant et al. (2018).

Simulation Overview

The collapse of a granular column is a classic benchmark for immersed granular flows, but here the grains themselves are not disks: each is a rigid line-cluster of 3–4 bonded disks, giving the assembly an orientation that can align, cluster, and jam in ways a circular grain never can. The solid phase is solved with NSCD exactly as elsewhere on this axis; what differs is that the fluid–grain interaction is computed disk by disk and then summed into a single net force and torque on the elongated grain, so that its own rotation naturally emerges from the drag imbalance across its length.

The fluid is solved on a coarser mesh through the same unresolved, porosity-averaged formulation used across the group,

with $\phi$ the porosity and $\mathbf{v}$ the volume-averaged fluid velocity, coupled semi-implicitly to the grains for stability. At each disk, the interaction force combines the pressure gradient with a Dallavalle-type drag corrected for the local porosity (Wen & Yu),

Columns spanning aspect ratios $H_0/L_0$ from 0.125 to 16, dry or immersed, and with either randomly oriented or pre-aligned grains, are released and left to collapse. Across the whole set, the final runout collapses onto a single scaling law driven by the front's kinetic energy at the moment it detaches — grain elongation shifts the prefactor of that law, but not its exponent, meaning shape changes how far a granular avalanche can travel without changing the fundamental physics that governs it.

Simulation Overview

The bay's cross-section is reconstructed as an inclined slope descending into a water-filled basin, with the rock mass represented as a population of large, rigid grains resolved individually by NSCD. The fluid is tracked with the same two-way PFEM–DEM coupling introduced in the Key Concepts above: the mesh is entirely rebuilt from an alpha-shape at every step and advected with the fluid's own velocity, so the free surface can deform without limit as the slide crashes in and the resulting wave races across the bay.

Because the slide is modelled with a comparatively small number of large grains rather than a grain-for-grain reconstruction of the real rockfall, a volume correction factor rescales their density so that the simulated mass matches the historical estimate of the slide directly — keeping the momentum delivered to the water physically consistent with the real event, even though the granular representation itself is necessarily coarse.

Based on wave-type analysis and the front Froude number, the resulting wave is expected to behave as a bore, and the simulation confirms it: the wave does not have time to break before crossing the bay and hitting the opposite shore, where it reflects and produces a secondary sloshing wave of smaller amplitude. Comparing the tracked wave amplitude at a point 885 m from the initial shoreline against recorded data, the Lagrangian formulation captures the main trends of the wave dynamics well. Displaced rocks settle throughout the bay, with the largest deposits forming near the impact zone and near the opposite shore.

Simulation Overview

This case validates the two-way fluid–grain coupling at the fully resolved end of the axis: a rigid 2D cylinder is released from rest in a quiescent viscous fluid and left to fall under gravity, with the mesh refined enough to resolve the boundary layer directly around it (roughly 8 elements across its diameter), rather than reconstructing it through a drag correlation. Results are compared directly against the reference numerical study of Namkoong et al.

With a solid-to-fluid density ratio of only 1.01, the cylinder is close to neutrally buoyant, so its descent is gentle enough for the wake dynamics — rather than gravity alone — to dominate its trajectory. As it accelerates, vortex shedding sets in and begins to modulate the fall.

The streamwise velocity initially overshoots the terminal velocity $U_t$, then is damped down as shedding becomes established; the lateral velocity oscillates with an opposite-sign pattern, its direction set by small mesh-induced asymmetries. Both the wavelength and frequency of these oscillations match the reference closely, confirming that the coupling captures not just the mean settling speed but the full unsteady, wake-driven dynamics of the descent.