CFD–DEM: Coupled Fluid & Particle Dynamics

• CFD–DEM is an Eulerian–Lagrangian framework that models fluids as continua and particles via discrete elements for high-fidelity multiphase simulations.
• It employs various formulations—from unresolved and coarse-grained to particle-resolved methods—to accurately capture hydrodynamic interactions, drag forces, and porosity effects.
• CFD–DEM is applied to sediment transport, fluidized beds, and granular flows, advancing reactor design, geophysical analysis, and multiphysics integration.

Computational Fluid Dynamics–Discrete Element Method (CFD–DEM) is an Eulerian–Lagrangian simulation framework in which the carrier fluid is modeled as a continuum by CFD while individual particles are advanced by DEM. In the formulations represented across recent arXiv literature, CFD–DEM spans unresolved volume-averaged models, resolved immersed-boundary or fictitious-domain models, semi-resolved coarse-grained variants, and mesh-less implementations, with applications ranging from sediment transport and fluidized beds to saturated granular failure, hydraulic conveying, reacting particles, heat transfer, and erosive wear (Chhushyabaga et al., 1 Aug 2025).

1. Scope and methodological variants

CFD–DEM is not a single algorithm but a family of couplings that share the same division of roles: the fluid is solved in an Eulerian frame, while particles are treated as Lagrangian bodies obeying Newton’s laws. In unresolved formulations, the fluid equations are volume-averaged over cells larger than the particles, so particle effects enter through porosity, solid volume fraction, and interphase momentum exchange terms. In resolved formulations, the flow around each particle is spatially resolved on the CFD grid, and hydrodynamic forces are obtained from local pressure and viscous fields rather than from drag closures (Fonceca et al., 2021).

This distinction is fundamental. Unresolved CFD–DEM is the dominant choice for large particle counts and dense engineering or geophysical systems, because it replaces particle-scale hydrodynamics by closure relations. Resolved CFD–DEM, including immersed-boundary and fictitious-domain variants, is used when near-particle flow detail is indispensable, for example in confined settling or very-narrow hydraulic conveying, where each particle spans multiple fluid cells (Filho et al., 2021). Between these extremes, semi-resolved approaches attempt to reduce grid dependence and recover finer local volumetric effects without paying the full cost of particle-resolved flow; examples include point-cloud-based two-step mapping and Voronoi-based two-grid coarse-graining (Liu et al., 11 Jun 2025, Che et al., 2022). Mesh-less CFD–DEM, represented by SPH–DEM, replaces the finite-volume grid by SPH particles while retaining the same unresolved locally averaged fluid model (Markauskas et al., 2016).

A recurrent misconception is that CFD–DEM is defined only by drag-coupled unresolved finite-volume solvers. The recent literature shows a broader landscape: OpenFOAM-, MFIX-, CFDEM-, SediFoam-, and finite-element-based implementations coexist, and the choice among them depends on the required fidelity in hydrodynamics, contact mechanics, coarse-graining, and multiphysics closure (Sun et al., 2015, Geitani et al., 2022).

2. Governing equations and phase descriptions

A representative unresolved formulation writes fluid mass conservation in terms of the fluid volume fraction $\varepsilon_f$ and the volume-averaged fluid velocity $\mathbf{v}_f$ as

$$\frac{\partial (\varepsilon_f \rho_f)}{\partial t} + \nabla \cdot (\varepsilon_f \rho_f \mathbf{v}_f) = 0,$$

with fluid momentum in volume-averaged form,

$$\frac{D}{Dt}(\varepsilon_f\rho_f\mathbf{v}_f) = \nabla \cdot \mathbf{S}_f + \varepsilon_f\rho_f\mathbf{f} - \sum_{m=1}^{M} \mathbf{I}_{fm},$$

where $\mathbf{S}_f$ is the fluid stress tensor and $\mathbf{I}_{fm}$ is the interphase momentum exchange term (Chhushyabaga et al., 1 Aug 2025). Closely related locally averaged forms are used in sediment transport and unresolved DEM–FVM/DEM–SPH studies, with porosity entering both continuity and momentum and with $\mathbf{F}^{fp}$ or $\mathbf{S}_{pf}$ acting as the fluid–particle source term (Sun et al., 2015, Markauskas et al., 2016).

Particles satisfy translational and rotational balance laws. A representative DEM formulation is

$$m^{(i)} \frac{d\mathbf{V}^{(i)}(t)}{dt} = \mathbf{F}^{(i)}_T(t), \qquad I^{(i)} \frac{d\boldsymbol{\omega}^{(i)}(t)}{dt} = \mathbf{T}^{(i)}(t),$$

with total force composed of gravity, fluid drag, and contact forces (Chhushyabaga et al., 1 Aug 2025). Contact models in the cited literature are predominantly soft-sphere. Depending on the application, they are Hertzian with Coulomb friction, or linear spring–dashpot with Coulomb friction, often with damping calibrated by a restitution coefficient (Chhushyabaga et al., 1 Aug 2025, Sun et al., 2015).

Hydrodynamic loading depends on the chosen coupling level. In unresolved schemes, drag is often written in relaxation-time or correlation form, for example

$$\mathbf{F}_d^{(i)} = m_p \frac{\mathbf{V}_p - \mathbf{u}_f}{\tau_p},$$

with $\mathbf{v}_f$0 determined by fluid properties, particle diameter, and local porosity through a drag correlation such as Gidaspow-, Wen–Yu-, Syamlal–O’Brien-, or Di Felice-type closures (Chhushyabaga et al., 1 Aug 2025, Che et al., 2022). In resolved CFDEM coupling, the fluid–particle interaction force is reconstructed from local pressure gradients and viscous terms over the particle occupation region,

$\mathbf{v}_f$1

so buoyancy is embedded in the pressure field rather than added as a separate DEM term (Fonceca et al., 2021).

3. Coupling, porosity mapping, and time integration

The technical core of CFD–DEM is the bidirectional transfer between particles and fluid cells. In a tightly coupled Eulerian–Lagrangian scheme, DEM deposits particle volumes and interaction forces onto the CFD mesh, producing $\mathbf{v}_f$2, solid volume fraction, and momentum source terms; CFD then interpolates fluid velocity and pressure back to particle positions to evaluate drag and related forces (Chhushyabaga et al., 1 Aug 2025). In sediment transport solvers such as SediFoam, this mapping is performed by diffusion-based coarse-graining, which yields smooth fields of solid volume fraction, solid-phase velocity, and fluid–particle interaction source on the CFD mesh and permits CFD cells smaller than particle diameters while maintaining mesh-independent averaged fields (Sun et al., 2015).

Because porosity enters both fluid equations and drag closures, the mapping problem is not a secondary implementation detail. Conventional direct particle-to-grid schemes can become strongly grid dependent, unstable as particles move across fluid grids, and unable to capture pore-fluid-pressure effects in very dense granular systems. Two recent alternatives target precisely this issue. One introduces a point-based coarse-graining layer using a multi-layer Fibonacci point cloud: particles are first converted into smooth point-cloud continuum fields independent of the CFD mesh, and only then projected to the fluid grid (Liu et al., 11 Jun 2025). Another defines local porosity for each particle through a radical Voronoi tessellation and transfers these values to the CFD mesh using a fine regular point cloud; this conserves porosity data, is reasonably grid-independent, and yields a relatively smooth porosity field, with improved prediction of fluid–particle interaction forces in polydisperse systems (Che et al., 2022).

Temporal coupling is equally consequential. DEM contact dynamics impose time scales far smaller than fluid time scales, so DEM is commonly subcycled within each CFD step. A local adaptive particle-advance strategy based on orthogonal recursive bisection resolves collisional time scales only for particles near potential collision partners and was reported to be $\mathbf{v}_f$3 faster than traditional explicit methods for problems that involve both dense and dilute regions while maintaining the same level of accuracy (Sitaraman et al., 2018). This suggests that the dominant computational bottleneck in many CFD–DEM simulations is not only fluid solving, but the global imposition of the smallest collision time step on the entire particle population.

4. Particle representation, constitutive enrichment, and multiphysics extensions

The simplest CFD–DEM representation uses rigid spheres, but several lines of work extend particle kinematics or internal physics without abandoning the Eulerian–Lagrangian structure. For natural sediment, one approach represents disk-, blade-, rod-shaped and equant grains with a small number of non-overlapping spherical components rigidly bonded into composite particles. Fluid forces are computed on each component sphere and summed, which restores physically consistent fluid torques while retaining sphere-based contact detection and drag/lift closures (Sun et al., 2016). This is a particle-representation change rather than a change to the fluid equations, but it materially alters settling, incipient motion, and transport behavior.

At the opposite end, bonded-particle representations convert a continuum solid object into a deformable assembly of bonded DEM spheres. In erosive-wear simulations, a bond network reproduces material properties at the microscale, while fluid acts on the bonded particles through CFD–DEM coupling. Surface deformation, cracking, and material removal then emerge from bond damage and fracture, allowing the geometry to evolve during erosion rather than being updated by an external wear law (Nguyen et al., 15 Feb 2025). This effectively embeds a discrete structural model inside CFD–DEM and creates a three-way fluid–structure–particle interaction.

Multiphysics extensions are now substantial. For reacting particles, intra-particle transport and reaction can dominate cost, so reduced-order models based on proper orthogonal decomposition and Galerkin projection have been developed for DEM/CFD calcination; replacing finite-volume particle models by reduced models reduced overall reactor simulation time by about $\mathbf{v}_f$4 in the treated example (Reineking et al., 2023). For thermal gas–solid systems, an Eulerian solid-phase heat-transfer model has been proposed in which particle motion remains DEM-based, but solid conduction is solved as a field on the Eulerian grid. This makes solid-phase conduction independent of DEM contact stiffness and strongly compatible with coarse-grained DEM, addressing a longstanding inconsistency of contact-based thermal DEM models (Imatani et al., 16 Jun 2025).

5. Applications, rheology, and validation practice

The application range of CFD–DEM now spans both canonical benchmarks and strongly nonlinear process- or hazard-driven flows. In sediment transport, CFD–DEM has been used to simulate flat bed in motion, small dune, vortex dune, and suspended transport regimes in a single framework, with validation against experimental and numerical benchmark data for transport rate, dune geometry, migration speed, concentration profiles, and Reynolds stresses (Sun et al., 2015). Shape-resolving bonded-sphere CFD–DEM further captures terminal velocity, threshold of incipient motion, and transport rate of natural sediments, showing that non-spherical grains rotate and translate differently from equivalent spheres (Sun et al., 2016).

In granular geomechanics, fully coupled CFD–DEM has been used to study collapse and runout of dense and loose granular columns under dry and fully saturated conditions. These simulations show pore-pressure-controlled differences between dilation-stabilized dense assemblies and compaction-driven loose assemblies, and further relate local friction to inertial and viscous numbers in rate-strengthening regimes (Chhushyabaga et al., 1 Aug 2025). A related pressurization-driven study of confined granular shear layers under imposed shear stress shows that instability is governed by the coupled evolution of effective stress, drainage, dilation or compaction, hydraulic connectivity, and granular fabric; viscous-number scaling organizes part of the low-$\mathbf{v}_f$5 creeping response, but not onto a unique local rheology (Chhushyabaga et al., 4 Jun 2026). This demonstrates one of CFD–DEM’s distinctive strengths: it can connect contact-network evolution, porosity redistribution, and pore-pressure fields within a single simulation.

Industrial and transport applications provide additional validation modes. Resolved CFD–DEM has been applied to hydraulic conveying through a very-narrow $\mathbf{v}_f$6 elbow with $\mathbf{v}_f$7, where it reveals particle settling, persistent crystal-like lattices in the elbow region at smaller velocities, and the corresponding contact networks (Filho et al., 2021). Fluidized-bed CFD–DEM remains a canonical benchmark for reduced-order modeling; a non-intrusive PODI ROM trained on CFD–DEM snapshots was shown to reconstruct Eulerian fluid volume fraction and particle observables for a parametrized fluidized bed, with local PODI outperforming global PODI in a Stokes-number study (Hajisharifi et al., 2023).

Benchmarking practice in this literature is notably heterogeneous but systematic. Common tests include single-particle sedimentation, hindered settling, Ergun pressure-drop tests, fluidized beds, spouted beds, current-induced dune formation, immersed granular column collapse, and pore-pressure evolution against analytical poromechanics solutions (Liu et al., 11 Jun 2025, Chhushyabaga et al., 1 Aug 2025). This diversity reflects the fact that CFD–DEM accuracy depends simultaneously on hydrodynamic closure, contact mechanics, coarse-graining, and temporal integration; no single benchmark is sufficient.

6. Limitations, controversies, and current directions

Several limitations recur across the field. First, unresolved CFD–DEM remains sensitive to how porosity and interaction terms are mapped. Direct grid-based mapping can produce strong grid dependence and temporal oscillations as particles cross cell boundaries, and in very dense granular systems it can fail to capture subtle volumetric changes that control pore pressure (Liu et al., 11 Jun 2025). Second, unresolved formulations necessarily average pore-scale hydrodynamics; in confined shear layers or dense suspensions, lubrication, local pore geometry, and anisotropic permeability are not explicitly resolved, even when pressure diffusion and drainage dominate the mechanics (Chhushyabaga et al., 4 Jun 2026).

Resolved methods address part of this deficit, but at high cost and with their own restrictions. The open literature emphasizes that to capture fluid-induced rotation of non-spherical grains, the CFD mesh must resolve velocity gradients at the particle scale; otherwise torques are under-predicted and particles slide rather than roll (Sun et al., 2016). Even in resolved CFDEM settling, confinement corrections can be qualitatively correct while quantitatively underestimating creeping-flow drag near walls (Fonceca et al., 2021). Mesh-less SPH–DEM avoids meshing but exhibits pressure fluctuations, whereas DEM–FVM gives smoother pressure fields but can suffer from mesh-induced porosity jumps (Markauskas et al., 2016). These are not contradictions so much as manifestations of different discretization errors.

A further controversy concerns what CFD–DEM should be expected to predict constitutively. The granular-flow studies reviewed here show that local friction can correlate with inertial or viscous numbers, but in creeping, heterogeneous, drainage-sensitive regimes no unique local rheology emerges (Chhushyabaga et al., 4 Jun 2026). Likewise, thermal and reactive extensions show that naively coupling new physics at the contact level can introduce spurious dependence on DEM numerical parameters such as spring stiffness; moving those closures to Eulerian or reduced-order descriptions is an active corrective trend (Imatani et al., 16 Jun 2025, Reineking et al., 2023).

Current development directions therefore target both fidelity and scalability. These include monolithic finite-element CFD–DEM with dynamically load-balanced parallelization and support for high order schemes (Geitani et al., 2022), local adaptive DEM time stepping (Sitaraman et al., 2018), non-intrusive reduced-order models for parametrized CFD–DEM (Hajisharifi et al., 2023), and grid-robust two-step coarse-graining strategies (Liu et al., 11 Jun 2025). Taken together, these efforts indicate that contemporary CFD–DEM is evolving from a fixed unresolved drag-coupled methodology into a broader computational framework in which representation, mapping, multiphysics closure, and solver architecture are all active research variables.