Volume Averaged Navier-Stokes

Volume-Averaged Navier–Stokes Equations

The fluid is not solved around each grain individually. Instead, following Anderson's framework, the pointwise Navier–Stokes equations are passed through a filtering (averaging) operation over a kernel large enough to span several grains, producing a fluid volume fraction $\epsilon$ (with $\phi = 1-\epsilon$ the solid fraction) and a handful of closure terms standing in for the unresolved grain-scale physics. Once those closures are specified, the resulting Volume-Averaged Navier–Stokes (VANS) system reads:

where $\mathbf v$ is the averaged grain velocity and $\boldsymbol\tau$ the viscous stress. Two of the unclosed terms are simplified directly: the residual (Reynolds-like) stress from velocity fluctuations is neglected at the low particle Reynolds numbers considered, and the filtered viscous stress is kept as the plain Newtonian $\epsilon\boldsymbol\tau \approx \eta(\nabla\mathbf u+\nabla\mathbf u^{\mathsf T})$ — any extra dissipation from dense packing is left to the discrete grains themselves, rather than folded into an effective fluid viscosity, so it isn't counted twice. The remaining term, the grain–fluid interaction force $\mathbf f_d$, is where the real closure modelling happens.

This is a Dallavalle-type drag law, corrected for neighbouring grains through the Di Felice voidage exponent and valid from the dilute to the dense limit.

Stabilized Finite Element Discretization

Velocity and pressure are both discretised with linear (P1–P1) continuous Galerkin shape functions on the same mesh — simple and cheap, but equal-order interpolation does not satisfy the inf–sup (LBB) condition on its own, which would otherwise produce spurious pressure oscillations. Two consistent stabilisation terms are added to the weak form to make the discretisation well-posed: a Pressure-Stabilising Petrov–Galerkin (PSPG) term restores a solvable pressure problem, a Streamline-Upwind Petrov–Galerkin (SUPG) term stabilises the advection-dominated regime.

Representing the Grains on the Fluid Mesh

Three quantities carry the whole grain–fluid interaction: the solid fraction $\phi$ (with $\epsilon=1-\phi$ the fluid fraction), the solid flux $\phi\mathbf{v}$, and the interaction force $\phi\mathbf{f}$. Since the fluid is already discretised with the Finite Element Method, an $L^2$ projection is the natural way to bring each of them from the grain scale onto the mesh — testing the grain-scale quantity against a nodal test function $\tau$ and integrating it over each grain's own volume $V_p$:

Left as it is, this projection can overshoot or undershoot the true solid fraction; MigFlow lumps the porosity's mass matrix so that $\phi$ — and with it $\epsilon$ — stays bounded between 0 and 1 at every node, rather than only in an averaged sense.

How the grain-volume integral $\int_{V_p}$ is actually evaluated is where accuracy used to break down. The original Particle-Centroid Method (PCM) evaluates it at a single point — the grain's centre — which is cheap but breaks down whenever a grain spans several elements: historically this forced a rule of thumb, keeping the grain diameter at least three times smaller than the mesh size, and still required an unphysical artificial pressure-diffusion term just to keep the pressure field stable. MigFlow now integrates instead: each grain is decomposed into the sub-regions defined by its intersection with every mesh element it overlaps, and the projection is evaluated over these exact sub-elements rather than a single point. Because each sub-region's shape is known exactly, any linear field — the shape functions $\tau$ included — can be integrated over it exactly, with no quadrature error. This removes the artificial diffusion entirely, spreads the interaction force smoothly as a grain crosses element boundaries instead of jumping between them, and keeps the void fraction bounded and physically consistent — valid precisely because NSCD already guarantees grains never overlap.

Semi-Implicit Coupling

The drag force $\mathbf f_d$ enters the momentum equation as a source term that can dwarf every other one, so coupling it explicitly to the grain motion would force impractically small time steps to stay stable. Instead, the grain velocity is predicted implicitly, within the same linear system as the fluid unknowns:

which, once substituted back into the drag force, effectively divides the drag coefficient by $(1+\gamma^n\Delta t/m)$: the stiffer the local coupling, the more it self-limits, so the drag stays bounded and always dissipative no matter how small the local void fraction gets. A floor $\epsilon_{\min} = 10^{-8}$ is enforced for the rare cells fully covered by a grain, keeping the pressure problem well-posed even in that limit — the same mechanism behind volume-penalisation methods.

A Two-Dimensional Approximation

Most simulations are run in two dimensions, with spherical grains modelled as infinite cylinders — discs in the simulation plane. This is a genuine approximation, not a free simplification: a random packing of discs settles at a packing fraction $\phi \approx 0.82$, well above the $\phi \approx 0.64$ of real spheres, since discs cannot nest on top of one another out of plane. To compensate, the fluid solver uses a smaller effective grain diameter when computing the local volume fraction $\epsilon$, while the grains' own dynamics still uses their true diameter and are restricted to three degrees of freedom rather than six. The result captures the essential physics well enough to be reused, with calibration, in later chapters — but it is not intended to be quantitatively three-dimensional.