1. Introduction
Particle suspensions are abundant in industrial and natural systems, from cosmetic and food products to human blood. Dense particle suspensions can exhibit a wide range of complex dynamical behaviours, owing to particle interactions mediated by the suspending fluid. From a continuum modelling perspective, their effective rheology depends not only on the local shear rate, but also, among other things, on the local particle volume fraction and particle pressure – a macroscopic representation of isotropic interparticle forces (Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). One common feature of particle suspensions is that a critical stress must be met to overcome friction and allow particles to rearrange; that is, suspensions have a frictional ‘yield stress’. A simple representation of this behaviour within a continuum framework is to impose that, below the yield stress, no shear can occur and the solid fraction is equal to a critical ‘jamming fraction’. Above the yield stress, the suspension can flow and the volume fraction drops below the jamming fraction. Numerous experimental measurements have captured the divergence of the effective suspension viscosity as the particle fraction approaches the jamming fraction from below (Zarraga, Hill & Leighton Reference Zarraga, Hill and Leighton2000; Boyer, Guazzelli & Pouliquen Reference Boyer, Guazzelli and Pouliquen2011; Dbouk, Lobry & Lemaire Reference Dbouk, Lobry and Lemaire2013; Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014; Dagois-Bohy et al. Reference Dagois-Bohy, Hormozi, Guazzelli and Pouliquen2015).
Suspensions flowing through confined geometries such as pipes can see their particle flow obstructed by the formation of a ‘clog’, which can broadly be defined as causing a local or global restriction in particle flux (Dressaire & Sauret Reference Dressaire and Sauret2017; Marin & Souzy Reference Marin and Souzy2024). This could be a local and temporary blockage of particle movement because of particles forming intermittent ‘bridges’ between walls, or a complete clog in which particle flux is zero everywhere. Even in free-flowing regimes, particle–wall interactions are an important factor in suspension flow. It is common to assume that particles and fluid are unable to move past walls because of frictional effects, leading to a non-uniform shear which is highest at the walls and zero around the centreline; in straight pipes, this leads to radially varying material properties and an unyielded central plug (Lyon & Leal Reference Lyon and Leal1998; Ahnert, Münch & Wagner Reference Ahnert, Münch and Wagner2019). In practice, particles are typically able to move past walls to some extent because of a fluid lubrication layer (Barnes Reference Barnes1995), allowing some particle flux even in an unyielded mixture (Szafraniec et al. Reference Szafraniec, Valdez, Iffrig, Lam, Higgins, Pearce and Wood2022, Reference Szafraniec, Bull, Higgins, Stone, Krüger, Pearce and Wood2025). The unique features of suspensions – wall slip, yield stress and clogging – lead to a complex relationship between particle volume fraction and particle flux, and typically the particle flux through channels is maximal at an intermediate volume fraction (Farutin et al. Reference Farutin2018).
The first aim of our work is to understand the implications of these features at a continuum scale as particle suspensions flow through gradually spatially varying geometries. Local geometrical constrictions have widely been found to promote clogging (Dressaire & Sauret Reference Dressaire and Sauret2017; Dincau, Dressaire & Sauret Reference Dincau, Dressaire and Sauret2023; Marin & Souzy Reference Marin and Souzy2024). More generally, changes in geometry often cause non-intuitive variation in particle volume fraction in dense suspensions, as in the phenomenon of ‘self-filtration’ (Haw Reference Haw2004; Kulkarni, Metzger & Morris Reference Kulkarni, Metzger and Morris2010), where driving a high-particle-fraction material through a constriction leads to a lower solid fraction downstream. Theoretical models of suspension flow through constrictions have typically taken an approach that resolves all particles and accounts for their stochastic motion (Parry & Millet Reference Parry and Millet2010; Mondal, Wu & Sharma Reference Mondal, Wu and Sharma2016; Bächer et al. Reference Bächer, Schrack and Gekle2017; Marin & Souzy Reference Marin and Souzy2024). However, in general it is not understood how particle dynamics connects to emergent (spatially varying) material properties, which drive continuum-scale suspension flow profiles and particle transport.
The second aim of our work is to explore whether analogies can be drawn between the effects on suspension flow of spatial variations in channel geometry and of spatial variations in particle properties. Biological suspensions containing certain constituents can exhibit such changes at the particle level, with substantial effects on overall suspension flow and function. An example is human blood, which is essentially a suspension containing deformable, fluid-filled capsules: red blood cells (RBCs). The properties of RBCs change as they age, and are also affected during certain drug treatments (Stathoulopoulos et al. Reference Stathoulopoulos, König, Ramachandran and Balabani2025) and by diseases such as malaria and diabetes. Here, we are interested in spatial variations in RBC properties, which can be drastic in sickle cell disease, a genetic disease that affects the haemoglobin inside RBCs, causing them to stiffen under deoxygenated conditions. Experiments on blood from patients with sickle cell disease have shown that deoxygenation can promote vessel occlusion (Wood et al. Reference Wood, Soriano, Mahadevan, Higgins and Bhatia2012; Szafraniec et al. Reference Szafraniec, Valdez, Iffrig, Lam, Higgins, Pearce and Wood2022), which is linked to poor clinical outcomes. More recently, large spatio-temporal variations in RBC volume fraction (or haematocrit) have been observed to emerge in channels with a deoxygenated region containing stiffened cells (Szafraniec et al. Reference Szafraniec, Bull, Higgins, Stone, Krüger, Pearce and Wood2025) – it is not clear whether this process is analogous to self-filtration in constricted channels. Phenomenological continuum models have been able to capture vessel occlusion in sickle cell disease by imposing the differential velocity between particles (RBCs) and the suspending fluid (Cohen & Mahadevan Reference Cohen and Mahadevan2013), but it is not understood how this differential flow emerges mechanistically, i.e. how changes to differential flow follow from changes to particle stiffness.
To achieve these two aims, here we build a minimal, mechanistic continuum model for the flow of a dense particle suspension down a pipe with a constriction in the axial direction. We also allow for varying particle properties via a spatially varying jamming fraction, which has been found to depend on deformability (Tapia et al. Reference Tapia, Hong, Aussillous and Guazzelli2024). For simplicity, we assume no slip at the walls. We allow for the emergence of variations in volume fraction, mediated by differential flow of the suspending fluid past the particles. We find that, for a given constriction or variation in particle properties, flow with a sufficiently low particle fraction can be sustained with a steady, spatially varying profile down the pipe. However, if the intended particle flux crosses a critical threshold, the model predicts the emergence of a high-particle-fraction, high-resistance, ‘clogged’ state, which develops from a free-flowing, unclogged initial state and propagates upstream.
2. Model
2.1. Mass and momentum conservation
The overall suspension consists of solid particles, filling a volume fraction
$\phi (\boldsymbol{x},t)$
, and interstitial fluid, occupying a fraction
$1-\phi$
. We assume the suspension is dense, in the sense that both hydrodynamic and contact forces contribute to the dynamics (Stickel & Powell Reference Stickel and Powell2005; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). The overall velocity
$\boldsymbol{U}$
in such a suspension can be broken into the respective local solid and fluid velocities
$\boldsymbol{u}_s$
and
$\boldsymbol{u}_f$
as
$ \boldsymbol{U} = \phi \boldsymbol{u}_s + (1-\phi )\boldsymbol{u}_f$
(see schematic in figure 1
a). Assuming each phase is individually incompressible, conservation of mass in each phase requires that
\begin{equation} \frac {\partial \phi }{\partial t} + \boldsymbol{abla}\,\boldsymbol{\cdot}\,(\phi \boldsymbol{u}_s )= 0, \qquad \frac {\partial \phi }{\partial t} - \boldsymbol{abla \boldsymbol{\cdot }} [(1-\phi ) \boldsymbol{u}_f]= 0, \end{equation}
which combine to yield
$\boldsymbol{abla \boldsymbol{\cdot }U} = 0$
. Rather than specifying individual, coupled momentum equations for the separate solid and fluid phases (Baumgarten & Kamrin Reference Baumgarten and Kamrin2019), which can be difficult to compare with rheometric measurements, we instead choose to decompose
$\boldsymbol{U}$
into the motion of the bulk ‘wet solid’ phase, which tracks the solid particles and the interstitial fluid moving at the same speed, and thus has velocity
$\phi \boldsymbol{u}_s + (1-\phi ) \boldsymbol{u}_s = \boldsymbol{u}_s$
, and a Darcy seepage velocity,
$\boldsymbol{u}_D \equiv (1-\phi ) (\boldsymbol{u_f}-\boldsymbol{u}_s)$
, which captures the differential transport of fluid through the moving solid particles. Therefore, the part of the fluid motion that moves with the solid particles is tracked within the bulk velocity
$\boldsymbol{u}_s$
, and the part of the fluid motion that moves past the solid is tracked by
$\boldsymbol{u}_D$
. This decomposition, which is similar to the recent ‘mixture-theory’ approach of Lu et al. (Reference Lu, Li, Chen, Wu, Lei and Wang2025), leads to
\begin{equation} \boldsymbol{U} = \boldsymbol{u}_s + \boldsymbol{u}_D. \end{equation}
(a) A schematic of the ‘wet solid’ phase, with velocity
$\boldsymbol{u}_s$
, and the differential ‘Darcy’ phase, with velocity
$\boldsymbol{u}_f$
, in pipe flow with a total flow rate
$Q$
. (b) Equivalent flow resistance schematic with wet solid phase
$\mathcal{R}_s$
and differential flow
$\mathcal{R}_{\!D}$
in parallel.
Momentum conservation for the overall suspension, neglecting any inertial or gravitational effects, yields
\begin{equation} \boldsymbol{abla \boldsymbol{\cdot }\sigma } = 0, \end{equation}
where
$\boldsymbol{\sigma }$
is the total stress tensor for the suspension. This stress can be divided into its isotropic and deviatoric parts: the former consists of the fluid pressure
$p_{\!f}$
and the excess ‘particle pressure’
$p_s$
that arises due to the interactions of the solid particles, while the latter we label
$\boldsymbol \tau$
, leading to
\begin{equation} \boldsymbol{\sigma } = - (p_{\!f} + p_s ) \boldsymbol{I} + \boldsymbol{\tau }, \end{equation}
where
$\boldsymbol{I}$
is the identity tensor, and thus, from (2.3),
\begin{equation} \boldsymbol{abla \boldsymbol{\cdot }\tau } = \boldsymbol{abla } (p_{\!f} + p_s). \end{equation}
Based on the premise that the suspension remains in a dense regime, we assume that the stress state of the interstitial fluid is dominated by viscous (Darcy) drag on the particle suspension. As such, deviatoric contributions to the total stress from the motion of the differential fluid are neglected; the differential flow is instead controlled by Darcy’s law, which relates fluid pressure gradients to the Darcy seepage velocity
\begin{equation} \boldsymbol{u}_D = -\frac {K(\phi )}{\eta _f}\boldsymbol{abla } p_{\!f}, \end{equation}
where
$\eta _f$
is the fluid viscosity and
$K(\phi )$
describes the permeability of the suspension.
2.2. Suspension rheology and constitutive laws
It has been widely observed that both shear stress and particle pressure in a sheared dense suspension scale linearly with the strain rate
$\boldsymbol{\dot \gamma } = \boldsymbol{abla u}_s + \boldsymbol{abla u}_s^T$
(Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). Given this observation, we follow Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011) by assuming a tensorial rheological model for the components of the stress tensor in (2.4) of the form
\begin{equation} p_s = \eta _f\, \eta _n(\phi ) |\boldsymbol{\dot\gamma}|, \qquad \boldsymbol{\tau } = \eta _f\, \eta _s(\phi ) \boldsymbol{\dot \gamma }, \end{equation}
provided
$\phi \lt \phi _m$
, where
$\phi _m$
is the jamming fraction, with
$\boldsymbol{\dot \gamma } = \boldsymbol{0}$
otherwise. Here,
$\eta _n(\phi )$
and
$\eta _s(\phi )$
are the dimensionless normal and shear viscosity functions, respectively; these are constitutive functions which must diverge as
$\phi \to \phi _m$
. The tensorial magnitude
$|\boldsymbol{\dot \gamma } | = \sqrt {\dot \gamma _{ij}\dot \gamma _{ij}/2}$
denotes the second invariant of that tensor.
This way (2.7) of presenting the rheology gives the particle pressure and stress in terms of the solid fraction and strain rate. An entirely equivalent modelling approach is to write the solid fraction and the stress in terms of the dimensionless ‘viscous number’
$J \equiv \eta _f \left |\boldsymbol{\dot \gamma }\right |/p_s$
, which provides a comparison of characteristic time scales for particle rearrangement and shear (cf. Boyer et al. Reference Boyer, Guazzelli and Pouliquen2011), in the form
\begin{equation} \boldsymbol{\tau } = \frac {\mu (J)}{J} \boldsymbol{\dot \gamma }, \qquad \phi = \eta _n^{-1}(1/J), \end{equation}
where
$\mu (J) = J \eta _s[\phi (J)]$
. The two formulations in (2.7) and (2.8) are entirely equivalent.
Here, we take slightly simplified versions of the constitutive expressions presented by Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011)
\begin{equation} \eta _n(\phi ) = \frac {\phi ^2}{(\phi _m - \phi )^2} , \qquad \eta _s(\phi ) = 1 + \mu _1 \eta _n(\phi ) , \end{equation}
with both viscosity functions diverging quadratically as the solid fraction approaches its jamming value,
$\phi \to \phi _m$
, in agreement with numerous experimental measurements (Zarraga et al. Reference Zarraga, Hill and Leighton2000; Boyer et al. Reference Boyer, Guazzelli and Pouliquen2011; Dbouk et al. Reference Dbouk, Lobry and Lemaire2013; Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014; Dagois-Bohy et al. Reference Dagois-Bohy, Hormozi, Guazzelli and Pouliquen2015). The viscosity ratio, on the other hand, converges to a finite limiting ‘yield friction’
$\eta _s/\eta _n \to \mu _1$
in this limit. With these constitutive laws, (2.8) reduces to
\begin{equation} \left . \begin{array}{c} \boldsymbol{\tau } = ( \eta _f + \mu _1 p_s /|\boldsymbol{\dot \gamma }| ) \boldsymbol{\dot \gamma } \\[6pt] {\phi }/{\phi _m} = 1/( 1 + J^{1/2}) \end{array} \right \} \quad \text{for} \quad \left | \boldsymbol{\tau } \right | \geqslant \mu _1 p_s; \qquad J \equiv \frac {\eta _f |\boldsymbol{\dot\gamma}|}{p_s}, \end{equation}
with
$\boldsymbol{\dot \gamma } = \boldsymbol{0}$
and
$\phi = \phi _m$
otherwise. This representation makes explicit the existence of a ‘yield stress’ in the formulation (at
$\left |\boldsymbol{\tau }\right | = \mu _1 p_s$
), which is somewhat obscured in (2.7). Below this stress, the material is held at its maximum jamming fraction and cannot deform; above this stress the material deforms and dilates, and the particle fraction is forced to be lower than the maximum jamming fraction. In fact, the expression for
$\boldsymbol{\tau }$
in (2.10) resembles the classical Bingham viscoplastic law, with the only difference being that here the yield stress depends on the particle pressure
$p_s$
.
Note that taking the full expressions for the constitutive laws in (2.9) from Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011) leads to a system that qualitatively resembles (2.10), sharing the essential features of a pressure-dependent yield stress with viscous flow above yield, but with a slightly more complicated expression for
$\boldsymbol{\tau }$
that proves rather less analytically tractable in the subsequent analysis. Since our aim here is to provide a general framework to describe suspensions, without placing too much emphasis on specific rheological choices, we retain the simplified form (2.10) here.
In a similar vein, for the permeability of the packing we adopt the commonly used Carman–Kozeny law
\begin{equation} K(\phi ) = \hat {k} k(\phi ) =\frac {\hat {k}(1 - \phi )^3}{\phi ^2}, \end{equation}
which comfortingly vanishes when the void space vanishes and diverges as the particle fraction vanishes. The dimensional prefactor
$\hat {k}$
is expected to scale with the square of the particle size, and sets the magnitude of the resistance to Darcy flow. This constitutive law is an established model for spherical particles, and in the absence of any more specific particle shape is used illustratively throughout this work. In fact, as we will see, the qualitative conclusions of this work do not depend sensitively on this choice, and the behaviour we observe would also occur even if a constant value was used for the permeability throughout.
2.3. Pipe flow, pipe averaging and lubrication approximation
For unidirectional, steady flow in a straight pipe of radius
$R$
aligned with the
$z$
direction, the equations simplify dramatically: only the shear stress
$\tau = \tau _{rz}$
and the pressures
$p_{\!f}$
and
$p_s$
play a role, with (2.3) reducing to
\begin{equation} \frac {1}{r} \frac {\partial }{\partial r}( r \tau ) = \frac {\partial }{\partial z} (p_s + p_{\!f}) \equiv -G, \end{equation}
given a total pressure gradient
$G$
along the pipe. Integrating yields a linear stress profile across the pipe,
$\tau = - {G r}/{2}$
. Comparison of this linear profile with the rheology in (2.8) shows that the stress must inevitably fall below the yield value
$\mu _1 p_s$
at some radius, provided the particle pressure is non-zero, leading to the existence of a ‘plugged’ core in the pipe with particle fraction
$\phi _m$
, as has previously been observed (Ahnert et al. Reference Ahnert, Münch and Wagner2019).
Mass conservation (2.1) is more usefully presented in its integrated, or pipe-averaged, form,
\begin{equation} \frac {\partial \overline {\phi }}{\partial t} + \frac {\partial }{\partial z} \left ( \overline {\phi u_s} \right ) = 0, \qquad \frac {Q}{\pi R^2} = \overline {u_s} + \overline {u_D}, \end{equation}
where
$Q$
is the total flux of the overall suspension down the pipe,
$u_s$
and
$u_D$
denote the respective flows in the
$z$
direction, and the overbar denotes a radial average, so
$\overline {f} = (2/R^2) \int _0^R r f \, {\rm d}r$
. The associated total pressure drop along a stretch of pipe of length
$\hat {L}$
is
$\Delta p = \int _0^{\hat {L}} G \, {\rm d}z$
, which motivates the definition of the overall resistance to flow down the pipe,
\begin{equation} \mathcal{R}_o = \frac {\Delta p}{Q}. \end{equation}
In a similar manner, we can define the individual resistances from the wet solid and Darcy phases, respectively, as
$\mathcal{R}_s = \Delta p/ (\pi R^2 \overline {u_s})$
and
$\mathcal{R}_{\!D} = \Delta p/ (\pi R^2 \overline {u_D})$
; these individual resistances effectively act in parallel (see figure 1
b) to make up the total resistance
\begin{equation} \mathcal{R}_o = \frac {\Delta p}{Q} = \left [\frac {1}{\mathcal{R}_s}+\frac {1}{\mathcal{R}_{\!D}}\right ]^{-1}, \end{equation}
which follows from the expression for the total flux in (2.13). Note that these pipe-averaged equations are one-dimensional, which precludes any issues of ill posedness that can arise in granular and suspension models (Barker et al. Reference Barker, Schaeffer, Bohorquez and Gray2015).
It is straightforward to generalise these results beyond straight pipes, provided that variations in the pipe geometry occur over a length scale
$\hat {L}$
that is much larger than the characteristic radial scale
$\hat {R}$
of the pipe (i.e.
$\partial R/\partial z \ll 1$
). This is the realm of classical lubrication theory. After scaling radial lengths with
$\hat {R}$
, axial lengths with
$\hat {L}$
, velocities with an as-yet undetermined scale
$V$
, stresses and particle pressures with
$\eta _f V/\hat {R}$
and fluid pressure with
$\eta _f V \hat {L}/\hat {R}^2$
, and working under the assumption that
$\hat {R} \ll \hat {L}$
, we arrive at the following equations, written now in terms of dimensionless variables:
\begin{equation} G = -\frac {\partial p_{\!f}}{\partial z}, \qquad \tau = - \frac {G r}{2}, \qquad u_D = -Da \, k(\phi ) \frac {\partial p_{\!f}}{\partial z}, \end{equation}
to leading order in
$\hat {R}/\hat {L}$
, while the rheology reduces to
\begin{equation} p_s = \eta _n(\phi ) \bigg| \frac {\partial u_s}{\partial r} \bigg|, \qquad \tau = \eta _s(\phi ) \frac {\partial u_s}{\partial r}, \end{equation}
provided
$\phi \lt \phi _m$
. The Darcy number
\begin{equation} Da = \frac {\hat {k}}{\hat {R}^2}, \end{equation}
determines the role of the permeability, and compares the square of the typical particle size,
$\hat {k}^{1/2}$
, with the square of the pipe radius,
$\hat {R}$
; it should, therefore, be significantly smaller than unity for a continuum description of the suspension to be reasonable. The Darcy number in our model is a measure of the resistance of the Darcy phase; indeed, at the jamming fraction
$\overline {\phi }=\phi _m$
the suspension resistance
$\mathcal{R}_s$
diverges, leaving the overall resistance (2.15) equal to the Darcy resistance
$\mathcal{R}_{\!D}$
, with
$\mathcal{R}_{\!D}\sim Da^{-1}$
, in line with classical models of flow through porous media (see (A7)).
The pipe-averaged equations (2.13), written in terms of dimensionless variables, become
\begin{equation} \frac {\partial \overline {\phi }}{\partial t} + \frac {1}{R^2}\frac {\partial }{\partial z} \big(R^2 \overline {\phi u_s}\, \big) = 0, \qquad \frac {Q}{\pi R^2} = \overline {u_s} + \overline {u_D}, \end{equation}
where
$R(z)$
is the dimensionless pipe radius. The definitions of the resistances remain as in (2.15). Note that one can set the undetermined velocity scale
$V$
to scale out either the total flux
$Q$
or the total pressure drop
$\Delta p$
from the problem, but in either case the resistance in (2.15), which is the ratio of the two, is unchanged. From this point on, we focus on the case of an imposed total flux, and thus choose the velocity scale
$V$
so as to scale the dimensionless flux
$Q$
to unity (i.e.
$V = Q/\pi \hat {R}^2$
), although we briefly revisit this question of whether the flux or pressure drop is imposed in Appendix C. The solid flux that appears in (2.19) will prove an important quantity, and so we label it as
\begin{equation} \mathcal{F} = R^2 \overline {\phi u_s}. \end{equation}
Note that, given the choice of scaling the dimensionless flux to unity,
$\mathcal{F}$
here represents the fraction of the total volume flux that is made up of solid.
As a side note, a key simplification arising from these lubrication scalings is the first equation in (2.16): the suspension rheology (2.7) enforces that the particle pressure
$p_s$
must scale with the shear stress
$\tau$
, but the usual lubrication scalings enforce that the fluid pressure
$p_{\!f}$
is asymptotically larger (by a factor
$\hat {L}/\hat {R}$
). As such, particle-pressure gradients are negligible compared with fluid pressure gradients, and so can be ignored in the overall pressure gradient
$G$
: the suspension is pushed along the channel by fluid (pore) pressure gradients to leading order. As a result of this difference in scaling for the two pressures, (2.19) becomes a hyperbolic equation for
$\overline {\phi }$
. In principle, one could include the asymptotically small particle-pressure gradients, which would add a weak non-local and non-hyperbolic character to the model that should slightly smooth over sharp fronts and shocks when they arise.
In summary, the model comprises the evolution equation for
$\overline {\phi }(z,t)$
in (2.19), which, in turn, is determined by knowledge of the pipe-averaged quantities
$\overline {\phi }$
,
$\overline {u_s}$
,
$\overline {u_D}$
and
$\overline {\phi u_s}$
. These can be deduced from the local force balances across the pipe and the suspension rheology (2.16)–(2.17), as outlined at the start of the next section. The parameters of the model are the Darcy number
$Da$
, which controls the permeability of the suspension, and the maximum solid fraction
$\phi _m$
of the suspension, together with the other rheological constants and any imposed variation in the pipe radius
$R(z)$
(in § 3.5 we briefly consider other forms of variation along the pipe). Throughout this work, we take the limiting friction coefficient to be
$\mu _1 = 0.3$
, following Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011).
We take a fixed total flux, which in dimensionless form is scaled to give
$Q = 1$
, and consider a fixed inlet solid fraction
$\overline {\phi }_{in} = \overline {\phi }(z=0,t=0)$
at the start of the pipe. In a steady state, the model simply reduces to two algebraic statements of flux conservation: the total flux
$Q \equiv 1$
and the relative solid flux
$\mathcal{F}$
(2.20) are both fixed. The time-dependent model (2.19), on the other hand, is a hyperbolic first-order partial differential equation, details of which are discussed in § 3.4 below and Appendix B.
3. Results
3.1. Flow and volume fraction profiles in a straight pipe
From (2.16)–(2.17), we see that the ratio of the shear stress and the particle pressure yields a function of
$\phi$
alone,
$\left |{\tau }/{p_s}\right | = {\eta _s(\phi )}/{\eta _n(\phi )} = {G r}/{ 2 p_s}$
, provided
$\phi \lt \phi _m$
. Given the specific rheological functions in (2.9), this equation reduces to
\begin{gather} \phi (r,t) = \frac {\phi _m\sqrt {R}}{\sqrt {R} + \sqrt {\mu _w \left (r - R \mu _1/\mu _w\right )}} \quad \text{for} \quad r\gt r_c, \end{gather}
where
\begin{equation} \mu _w(t) \equiv \frac{\textit{GR}}{2 p_s}, \qquad \text{and} \qquad r_c(t) \equiv \frac {R \mu _1}{\mu _w} = \frac {2 p_s \mu _1}{G}, \end{equation}
with
$\phi (r\lt r_c) = \phi _m$
. Here, the grouping
$\mu _w$
combines the local pipe radius, pressure gradient and particle pressure, and can be thought of as an effective friction coefficient for the suspension at the pipe walls. The critical ‘yield’ radius
$r_c$
separates a central core of material with
$\phi = \phi _m$
from a flowing region with
$\phi \lt \phi _m$
given by (3.1). Example plots of the profiles of
$\phi$
predicted by (3.1) are shown in figure 2(a).
Flow properties in a straight pipe, with
$R=1$
unless otherwise specified. (a) Varying particle fraction and solid velocity profiles for increasing shear to normal stress ratios,
$\mu _w = \textit{GR}/2 p_s = 0.5,0.75,1.5$
, and (b) scaled average solid fraction
$\overline {\phi }/\phi _m$
as a function of
$\mu _w=\textit{GR}/p_s$
. (c) Average wet solid and Darcy phase velocities varying with
$\overline {\phi }/\phi _m$
for different
$Da$
. (d) Flow resistances (2.14)–(2.15) for
$Da = 10^{-2}$
varying with
$\overline {\phi }/\phi _m$
; the inset shows how the overall resistance for
$\overline {\phi }$
near to
$\phi _m$
increases as
$Da$
is reduced. (e) Solid flux showing a maximum at a value of
$\overline {\phi } \lt \phi _m$
, for different
$Da$
; the inset shows the scaled solid flux
$\mathcal{F}/\mathcal{R}_o$
, which is independent of
$Da$
, if the overall pressure gradient is fixed, rather than the overall flux. ( f) The maximum achievable solid flux, which decreases as the pipe radius decreases and as the Darcy number (permeability) increases.
Given
$\phi$
, the associated velocity profile
$u_s$
follows from the expression for
$\tau = -\textit{Gr}/2$
in (2.16), and thus the expression for
$\partial u_s /\partial r$
in (2.17). We assume here that the suspension experiences a no-slip condition (
$u_s = 0$
) at the walls of the pipe
$r=R$
, and as such can integrate
$\partial u_s /\partial r$
to give the yield-stresslike flow profile
\begin{gather} u_s(r,t) = \begin{cases} \dfrac {G}{4}(R - r_c)^2 \def\lumina\hspace {0.3cm} &\text{for $r \lt r_c$}\\[8pt] \dfrac {G}{4}\left [(R - r_c)^2 - (r - r_c)^2\right ] \def\lumina\hspace {0.3cm} &\text{for $r \geqslant r_c $}, \end{cases} \end{gather}
(the same expression can equivalently be reached directly by integrating the expression for the stress in (2.10) using (2.17)). Equation (3.3) describes a central plug of unyielded material (corresponding to where
$\phi =\phi _m$
), bordered by sheared regions, as illustrated in figure 2(a). This qualitative flow structure for dense suspensions in a channel or pipe is well known (Phillips et al. Reference Phillips, Armstrong, Brown, Graham and Abbott1992; Koh, Hookham & Leal Reference Koh, Hookham and Leal1994; Lyon & Leal Reference Lyon and Leal1998).
The cross-pipe profiles for
$\phi$
and
$u_s$
provide all the necessary information to determine the pipe-averaged quantities in (2.19). These integrals can be computed analytically, but are algebraically unpleasant: details of this integration and the form of the resultant average quantities are given in Appendix A. Note that the structure of the flow across the pipe depends only on the grouping
$\mu _w = \textit{GR}/2p_s$
, which sets the radius of the central plug region (3.2), rather than on the pressure gradient or particle pressure alone. It is, therefore, this grouping
$\mu _w$
that controls whether the suspension is able to flow or not: the suspension formally clogs (i.e. the solid flux goes to zero) only if the central plug fills the pipe,
$r_c \to R$
, in which limit
$\overline {\phi } \to \phi _m$
.
Figure 2(b) shows the dependence of the average solid fraction
$\overline {\phi }$
on
$\mu _w$
, which follows from the radial integral of (3.1). This plot illustrates directly the above point:
$\overline {\phi } \to \phi _m$
when
$\textit{GR} \leqslant 2\mu _1 p_s$
, as the central plug fills the entire pipe. The associated solid flux will also vanish in this limit, owing to the no-slip condition on the walls.
3.2. Role of the Darcy phase: flow resistance and maximum particle flux
The results so far are independent of the differential (Darcy) flow of fluid through the suspension. Its role is to modulate the total flux
$Q$
, through (2.19), which provides another constraint linking
$G$
and
$p_s$
, and thus determines the behaviour of the system explicitly (see Appendix A for the explicit expressions). Figure 2(c) shows the individual mean ‘wet solid’ and Darcy velocities as a function of the mean solid fraction, given that the total flux is fixed at
$Q=1$
. It is clear that the solid velocity dominates the flux until the solid fraction approaches its maximum packing. In this limit, the flow resistance from the wet solid phase is so great that fluid is forced through the packing instead: both the pressure gradient
$G$
and the particle pressure
$p_s$
grow significantly over a narrow range of
$\overline {\phi }$
in order to maintain the desired solid fraction whilst providing the necessary total flux. As such, the Darcy velocity grows near
$\overline {\phi } = \phi _m$
as the solid flux drops towards zero in this limit. The Darcy number
$Da$
controls the permeability of the packing, and thus the resistance to Darcy flow; hence the region in which the Darcy flow becomes dominant is increasingly localised to be close to
$\phi _m$
as
$Da$
is reduced (figure 2
c). Note that the pressure gradient
$G$
drives both the bulk wet solid motion and the Darcy seepage flow along the pipe, and thus for any
$G$
there is always some Darcy flow, but its contribution is generally negligible unless the solid fraction is near its maximum value (as can be seen in figure 2
c).
The competition between ‘wet solid’ flow and Darcy flow is perhaps better illustrated in figure 2(d), which shows the individual resistances for each phase and the overall flow resistance, as defined in (2.15). The ‘wet solid’ resistance is generally much lower than the Darcy resistance, but it increases dramatically as
$\overline {\phi } \to \phi _m$
. Because the two resistances effectively act in parallel (figure 1
b and (2.15)), the overall resistance (dashed line) becomes limited by the Darcy resistance in that limit instead.
As a consequence of this behaviour, the solid flux down the pipe (as a fraction of the overall flux)
$\mathcal{F} = R^2 \overline {\phi u_s}$
does not simply increase monotonically as
$\overline {\phi }$
is increased (figure 2
e): it has a maximum
$\mathcal{F} = \mathcal{F}_{max}$
at a solid fraction lower than
$\phi _m$
, before reducing again for high solid fraction and vanishing in the limit
$\overline {\phi } \to \phi _m$
, when the particles become stationary (because there is no slip of particles on the walls in our model). This means that any given solid flux (below the maximum) can be attained at two different values of the solid fraction: a low
$\overline {\phi }$
, with a correspondingly high solid velocity and relatively low resistance; or a high
$\overline {\phi }$
, with a correspondingly low solid velocity and high resistance. Again, lowering the Darcy number (permeability) increases the resistance to porous flow and tends to squash the ‘high-
$\overline {\phi }$
’ branch of the flux closer to
$\overline {\phi }= \phi _m$
(figure 2
e). Note that the same qualitative behaviour of a non-monotonic solid flux is also seen if the pressure gradient
$G$
is held fixed, rather than the total flux (inset to figure 2
e).
More generally, these solid-flux profiles
$\mathcal{F}(\overline {\phi }; Da)$
also vary with the pipe radius (figure 2
f), with smaller radii offering a lower maximum solid flux
$\mathcal{F}_{max}$
. The reason for this reduction in solid flux is that the relative contribution of the Darcy flux to the total flux increases in thinner pipes: given a pipe with cross-sectional area
$A\sim R^2$
, the solid flux scales with
$\textit{AGR}^2 \sim \textit{GR}^4$
, because
$u_s \sim \textit{GR}^2$
(from (3.3)), whereas the Darcy velocity is constant across the pipe and the Darcy flux scales like
$\sim AG \sim GR^2$
(see also Appendix A). This reduction in the attainable solid flux with radius will prove a crucial detail for interpreting the behaviour in constricted pipes where
$R$
, and thus all the quantities calculated here, can vary along the pipe.
Illustration of ‘non-clogging’ and a ‘clogging’ constrictions. (a) Solid flux
$\mathcal{F}(\overline {\phi })$
for
$\phi _m = 0.65$
,
$Da = 10^{-2}$
and different values of the radius
$R$
varying between
$R=1$
and
$R = R_{min} = 0.5$
. The red dashed line shows
$\mathcal{F}_{in}$
for
$\overline {\phi }_{in} = 0.35$
. (b) Steady solution
$\overline {\phi }(z)$
for an imposed radial constriction between
$R=1$
and
$R= R_{min} = 0.5$
as shown by the blue dashed line; the green crosses correspond to the crosses in (a). (c) The same flux curves as in (a), but now showing
$\mathcal{F}_{in}$
for a slightly larger inlet solid fraction
$\overline {\phi }_{in} = 0.45$
; there is no steady solution with this value of
$\overline {\phi }_{in}$
because the downstream constriction cannot sustain this flux. (d) Heat map of overall resistance against
$\overline {\phi }_{in}$
and
$Da$
for the same constricted system; the dashed line represents
$Da = 10^{-2}$
and the two stars correspond to the cases illustrated in (a) and (c).
3.3. Flow and particle flux in constricted pipes
We now consider pipes with an imposed radial constriction
$R(z)$
(figure 3), in which the radius decreases smoothly from
$R=1$
to
$R = R_{min}$
as shown in figure 3(b). For all cases shown here, we adopt a profile
\begin{equation} R(z) = \frac { \left (1+R_{min} \right ) - \left (1-R_{min}\right ) \,\text{erf}\left [\left (z-0.5\right )10\right ]}{2}. \end{equation}
Provided that gradients in the true (dimensional) pipe radius are sufficiently small, the evolution of solid fraction is described by (2.19), with the local flux at any location
$R(z)$
given by the flux for a straight pipe of that radius. As such, the steady-state flow is determined simply by a statement of flux conservation: the solid flux through each section of the pipe must be equal, which should allow for determination of the associated average profile
$\overline {\phi }(z)$
.
Figure 3(a,b) show such a steady solution. Given a value of the inlet solid fraction
$\overline {\phi }_{in}$
, the local solid fraction along the pipe
$\overline {\phi }(z)$
follows directly from the expression for the solid flux
$\mathcal{F}(\overline {\phi },R)$
evaluated at each local radius
$R(z)$
(figure 3
a). Since the solid flux is conserved, the variation in the volume fraction along the pipe can simply be read off from the intersection of the different solid-flux curves for each radius with the incident solid flux
$\mathcal{F}_{in}$
(green crosses in figure 3
a), leading to a profile
$\overline {\phi }(z)$
that increases slightly in the constricted part of the channel (figure 3
b) – this effect has been observed qualitatively in particle-scale simulations (Bächer et al. Reference Bächer, Schrack and Gekle2017).
However, this solution construction is not always possible. For sufficiently high inlet volume fractions
$\overline {\phi }_{in}$
, or sufficiently large constrictions, it is no longer possible to determine the local solid fraction using the incident solid flux and the local radius in this manner. This is the case when the solid-flux curve for the most constricted part of the pipe is everywhere lower than the solid flux set by the inlet solid fraction (as illustrated in figure 3
c), and so there is no steady solution to (2.19) that conserves the imposed inlet solid flux in that part of the channel.
A full parameter sweep illustrates that, as
$\overline {\phi }_{in}$
increases, the overall pipe resistance increases until the ‘no solution’ region of parameter space is reached; figure 3(d) shows this parameter space for a fixed radial constriction with
$R_{min} = 0.5$
. Note that, for sufficiently high
$\overline {\phi }_{in}$
, the associated inlet solid flux
$\mathcal{F}_{in}$
reduces (see flux curves in figure 3
a) and a steady solution can again be found. In this region of parameter space,
$\overline {\phi }$
is always on the upper branch of the flux curve: the solid fraction now decreases through the constriction, rather than increasing, and is everywhere fairly close to the jamming fraction, such that the overall flow resistance can be very large (figure 3
d). Physically, this region of parameter space appears because sufficiently near
$\overline {\phi }_{in} = \phi _m$
, the solid flux
$\mathcal{F}_{in}$
in the wider section of the pipe must fall below the maximum flux
$\mathcal{F}_{max}$
that can be sustained in the constricted section of the pipe.
3.4. Emergent spatio-temporal heterogeneity in constricted pipes
To understand the dynamics of the system in the regions of parameter space that do not appear to exhibit a steady solution (black region in figure 3
d), we solve the time-dependent problem (2.19) using the Clawpack finite-volume method for hyperbolic problems (Clawpack Development Team 2024) (see Appendix B for further details). We use an initial condition given by the steady state for
$\overline {\phi }_{in} = 0.2$
, corresponding to a solid flux
$\mathcal{F}_0$
, and consider a step increase in the injected solid fraction at
$z=0$
to a new (higher) value of
$\overline {\phi }_{in}(t=0)$
. We subsequently enforce reflection conditions (
$\partial \overline {\phi }/\partial z = 0$
) at
$z=0$
, rather than holding the inlet particle fraction fixed, in order to allow for variations in
$\overline {\phi }$
to emerge throughout the entire pipe.
Transient evolution of constricted flow. (a) Evolution from a low-solid-fraction steady state, with
$\overline {\phi }_{in} = 0.2$
, when
$\overline {\phi }_{in}(t=0) = 0.35$
(upper) and
$\overline {\phi }_{in}(t=0) = 0.45$
(lower), for
$\phi _m = 0.65$
,
$Da = 10^{-2}$
and
$R_{min} = 0.5$
(with
$0\leqslant z\leqslant 1$
for all panels). These solutions correspond to a ‘not clogged’ (upper) and emergent ‘clogged’ (lower) state, respectively, with symbols corresponding to those marked in (c). (b) Overall resistance
$\mathcal{R}_o$
over time for each case. (c) The solid flux
$\mathcal{F}$
for the widest and narrowest parts of the pipe, for the ‘not-clogged’ (left) and ‘clogged’ (right) examples in panel (a). The pre-existing steady state has flux
$\mathcal{F}_0$
that is initially increased at
$z=0$
to
$\mathcal{F}_{in}$
(green dot). In the right-hand case, the downstream flux (black cross) is too low, and the upstream solid fraction is forced to increase towards the red dot, reducing the eventual solid flux
$\mathcal{F}_{\infty }$
, as outlined in the main text.
The final steady-state (a) upstream (inlet) solid fraction and (b) total resistance as a function of the initial inlet solid fraction, for a constricted pipe with minimum radius as marked,
$Da = 10^{-2}$
and
$\phi _m = 0.65$
. The central line, which matches the parameters in figure 4, corresponds to traversing the yellow dashed line in the phase plot in figure 3(d). The discontinuities in each case represent the development of the upstream-propagating ‘clogged’ state.
We obtain temporal solutions in regions of parameter space both with and without a steady solution. For the former, the system reaches the steady state described above in which the solid fraction increases slightly within the constricted part of the channel (figure 4
a upper panels). In the ‘no solution’ region of parameter space, however, while the solid fraction again increases downstream, it now starts to ‘back-up’ solid, and forms a shock just upstream of the constriction. The shock then travels upstream all the way back to the inlet (figure 4
a lower panels), raising the solid fraction significantly everywhere upstream. Note that the existence of a shock – that is, an axial discontinuity in solid fraction – locally violates our lubrication approximation; in reality we might expect this backwards-propagating region to be smoothed out over a dimensional length comparable to the pipe radius, with locally more complex and non-unidirectional dynamics. The associated overall resistance steadily increases as the shock moves backwards, extending the region where the solid fraction – and thus the flow resistance – is high (figure 4
a,b). Thus, as the shock passes through the inlet the solid fraction
$\overline {\phi }_{in}$
increases to a much higher value, which we designate
$\overline {\phi }_{in}(t \to \infty )$
, and a new steady state is reached.
We refer to this new state – shown in figure 4(a–c) – as ‘clogged’. This is for two reasons: first, because of its very high solid fraction and overall resistance, and second, because of the way it emerges discontinuously to fill the pipe from an apparently free-flowing, moderate-solid-fraction state. In this manner, the emergence of the ‘clogged’ state provides a downstream control on upstream properties (via the backwards-propagating shock), in direct contrast to the ‘non-clogged’ state, in which the imposed upstream solid fraction
$\phi _{in}$
sets the profile down the pipe. For clarity, we note that in our model, the solid flux is not reduced to zero in this ‘clogged’ state, which has
$\overline {\phi }\lt \phi _m$
everywhere. However, the resistance increases so dramatically (figure 4
b) that an experimental apparatus imposing a fixed overall flux (as in our theory) might be expected to fail in this regime. Alternatively, the apparatus may be unable to adjust to a new, higher, inlet particle fraction, as occurs in our model, which would also lead to failure. If, instead, the pressure drop down the pipe were fixed, the overall flux would decrease significantly in the ‘clogged’ state, as briefly discussed in Appendix C.
The transient solutions indicate how to find this ‘clogged’ steady state directly. If the imposed solid flux is too large for the constricted pipe (i.e. in the ‘no-solution’ region of phase space; figure 3
d) then the flux is lowered until a solution can be found. Physically, this occurs by ‘backing up’ the excess solid upstream of the constriction, which pushes the upstream solid fraction up onto the high-
$\overline {\phi }$
branch of its flux curve. This process is illustrated in figure 4(c), which shows how the solid flux evolves in the pipe for each of the two cases considered in panel (a). In each case, the initial steady state carries a flux
$\mathcal{F}_0$
, and at
$t=0$
this is increased at
$z=0$
to
$\mathcal{F}_{in}$
. In the ‘not clogged’ case, this new flux can be sustained all along the pipe, and the suspension simply evolves towards the associated steady state; there is no subsequent change in the solid flux and
$\mathcal{F}_{\infty } = \mathcal{F}(t\to \infty )$
is simply
$\mathcal{F}_{in}$
. In the ‘clogged’ case, however, the flux is limited by the downstream maximum: incoming solid builds up upstream of the constriction, and the only resolution to this mismatch is for the upstream solid fraction to increase to the point where the flux is reduced to match the downstream limiting value. Thus
$\mathcal{F}_{\infty }$
is lower than
$\mathcal{F}_{in}$
, as shown in figure 4(c).
The eventual emergent state in this case can therefore be directly constructed by finding the limiting downstream flux, and following the high-
$\overline {\phi }$
flux branch upstream from this value. Equivalently, this steady state corresponds to the point where a horizontal line in a phase-space plot like figure 3(d) meets the right-hand boundary of the ‘no solution’ region.
We find that this construction holds for any initial condition within the ‘no solution’ region. Figure 5(a) shows how the eventual solid fraction at the start of the pipe,
$\overline {\phi }_{in}(t \to \infty )$
, varies with the initial solid fraction at the start of the pipe. These quantities are equal when the solid fractions are sufficiently small, but above some critical
$\overline {\phi }_{in}$
the fluxes down the pipe can no longer balance, as outlined above, and the system is forced into this new state with a much higher solid fraction throughout the pipe. The associated resistance is also significantly increased, as shown in figure 5(b). If the size of the constriction is increased (i.e.
$R_{min}$
decreased), the transition to this emergent ‘clogged’ state is more dramatic, and occurs at a lower inlet solid fraction (figure 5). Note that both the final inlet solid fraction and the resistance increase again for
$\phi$
very close to
$\phi _m$
; this region of parameter space corresponds to the steady states that lie entirely on the high-
$\overline {\phi }$
branch of the flux curves, as discussed at the end of the previous section and illustrated by the band of steady solutions close to
$\phi _m$
in figure 3(d).
(a) Evolving solutions for a spatially varying jamming fraction
$\phi _m$
decreasing from
$0.8$
to
$0.6$
(dashed line) along a pipe, with
$Da = 10^{-1}$
, uniform radius
$R=1$
and uniform initial conditions
$\overline {\phi }(t=0) = 0.35$
(upper row) and
$\overline {\phi }(t=0) = 0.5$
(lower row). (b) The associated evolution of the flow resistance in each case. (c) The solid flux in each case, comparing the flux for the highest and lowest values of
$\phi _m$
, with fluxes and symbols as in figure 4.
3.5. Emergent heterogeneity in straight pipes with variation in particle properties
Motivated by recent experiments on blood flow from patients with sickle cell disease (Szafraniec et al. Reference Szafraniec, Valdez, Iffrig, Lam, Higgins, Pearce and Wood2022, Reference Szafraniec, Bull, Higgins, Stone, Krüger, Pearce and Wood2025), in which RBCs stiffen under deoxygenated conditions, we also briefly consider suspension flow in straight pipes with spatial variation in particle properties. To represent such variations in the simplest way possible, we allow the maximum particle fraction
$\phi _m(z)$
to decrease smoothly from high to low through the pipe (figure 6
a), qualitatively capturing the reduction in jamming fraction for suspensions containing rigid particles in comparison with more deformable particles (Tapia et al. Reference Tapia, Hong, Aussillous and Guazzelli2024). In so doing, our aim is to explore whether analogies can be drawn between variation in particle properties and the geometrical constrictions studied above. We also aim to provide at least a basic qualitative explanation of experimental observations of large-scale spatial variations in haematocrit (i.e. particle fraction) for blood flow in channels with a downstream deoxygenated region that contains stiffened RBCs (Szafraniec et al. Reference Szafraniec, Bull, Higgins, Stone, Krüger, Pearce and Wood2025).
We find that steady solutions of (2.19) with varying
$\phi _m(z)$
exhibit a bifurcation analogous to the one observed for channel constrictions: if the incident particle fraction
$\overline {\phi }_{in}(t = 0)$
is too large relative to the downstream maximum particle fraction, then the incident solid flux
$\mathcal{F}_{in}$
cannot be maintained, and a shock forms and propagates upstream. This behaviour is illustrated in figure 6(a), which shows results from time-dependent simulations in which an initially uniform solid fraction is exposed at
$t=0$
to a prescribed reduction in the jamming fraction
$\phi _m$
along the pipe. If the initial solid fraction is sufficiently low (upper row), the suspension evolves smoothly to a steady state that carries the same solid flux
$\mathcal{F}_{in}$
as before. However, larger initial solid fractions (lower row) instead result in a shock developing and spreading upstream, leading to a final state with a much higher solid fraction, that carries a lower solid flux
$\mathcal{F}_\infty$
. As in the case of a constriction, the associated overall flow resistance increases markedly in this ‘clogged’ state (figure 6
b). The final states can again be constructed simply by consideration of the solid fluxes, as illustrated in figure 6(c).
These results suggest that the spatial variations observed in previous experiments on blood flow (Szafraniec et al. Reference Szafraniec, Bull, Higgins, Stone, Krüger, Pearce and Wood2025) may be explained by a reduction in maximal RBC flux as cells become stiffened under deoxygenation. More broadly, the results demonstrate that our predictions of emergent heterogeneity in solid fraction apply to a broad range of synthetic and biological particle suspension systems in which either geometry or particle properties vary in space.
4. Discussion
We have built a mechanistic continuum model that allows for spatio-temporal variations in the volume fraction of particle suspensions flowing through a pipe. A key feature of our approach is the separation of suspension flow into a ‘wet solid’ phase, which tracks particles moving with the suspending fluid, and a differential Darcy phase, which tracks flow of suspending fluid past the particles; the constitutive properties of each of these phases can, in principle, be directly measured.
In pipes with a sufficiently small constriction or low solid fraction, we find steady solutions in which a given upstream solid fraction increases as it passes through the constriction, maintaining the same solid and total flux. A key finding here is that significant spatio-temporal volume fraction variations can emerge spontaneously if the constriction is large enough, resulting in an abrupt transition to a high-particle-fraction, high-resistance ‘clogged’ state, which develops at the entrance to the constriction and propagates back upstream. The resultant steady state in this case is controlled by the downstream constriction, rather than the upstream particle fraction, and has a particle fraction that drops, rather than increases, through the constriction (akin to observations of self-filtration; e.g. Kulkarni et al. Reference Kulkarni, Metzger and Morris2010).
The transition to a ‘clogged’ state follows directly from the existence of a maximal solid flux at an intermediate particle fraction that depends on geometry and particle properties, and can therefore vary in space. More specifically, the emergence of ‘clogging’ can occur whenever the flux in one region (e.g. a wider region of the pipe) is higher than the maximal flux in another region (e.g. a constriction; see figure 3). There are numerous parallels and qualitative similarities in these predictions to models of constrictions in traffic flow (e.g. Lighthill & Whitham Reference Lighthill and Whitham1955); these, in their basic form, also comprise hyperbolic advection problems with non-monotonic flux curves, which can result in the upstream build-up of traffic density – that is, traffic jams – via backwards-propagating shocks if the maximum flux in the ‘constricted’ region is not sufficiently large.
In our model, the presence of this maximum solid flux emerges from the constitutive assumptions. In a suspension with fixed overall flux, without a differential Darcy flow, the flow speed of the ‘wet solid’ phase would be fixed and the solid flux would linearly increase with particle fraction (the flux of interstitial fluid that moves at the same speed as the particles would reduce). However, in our model we include differential flow of fluid through the pore space of the bulk suspension, and the proportion of the overall flux taken up by each phase is determined by their relative resistances, which effectively act in parallel (figure 1). As the solid fraction approaches its maximal jamming fraction, the resistance of the wet solid phase increases drastically, and more of the flux is pushed through the differential Darcy phase, leading to a reduction in the particle flux (figure 2). The emergence of the ‘clogged’ state in a constricted channel is a consequence of the fact that the maximal solid flux decreases as the pipe radius decreases. The wet solid resistance has a stronger dependence on radius than the resistance to Darcy flow (see § 3.2); Darcy flux is thus promoted for smaller radii, leading to a reduction in the attainable solid flux there.
In addition to studying constrictions, we have demonstrated that, in principle, variations in particle properties (parameterised here by variation in the jamming fraction
$\phi _m$
) can similarly lead to the emergence of ‘clogging’. This provides a qualitative explanation for recent experimental results on blood flow from patients with sickle cell disease (Szafraniec et al. Reference Szafraniec, Valdez, Iffrig, Lam, Higgins, Pearce and Wood2022, Reference Szafraniec, Bull, Higgins, Stone, Krüger, Pearce and Wood2025), in which spatial variations in haematocrit were observed to emerge because of deoxygenation-induced stiffening of RBCs in downstream regions of flow. Specifically, our results show that the mechanism for these results might be the same as previous results on clogging in constricted channels: the particle (or cell) flux in certain regions of the channel is higher than the maximal flux in other regions. These results will contribute to our understanding of vaso-occlusion and other pathological outcomes related to blood rheology in sickle cell disease.
Our model is continuum and deterministic; it therefore differs in approach from stochastic and discrete theories for the effect of constrictions on particle suspensions (Zuriguel et al. Reference Zuriguel2014; Marin et al. Reference Marin, Lhuissier, Rossi and Kähler2018; Souzy, Zuriguel & Marin Reference Souzy, Zuriguel and Marin2020), but can be interpreted as parameterising particle-scale effects mechanistically over long length and time scales. As an example, previous work on clogging has identified a probability of particles forming a temporary bridge between confining walls (‘bridging’), contributing to a local build up of particles (Dressaire & Sauret Reference Dressaire and Sauret2017; Vani, Escudier & Sauret Reference Vani, Escudier and Sauret2022). To directly capture such effects, including the intermittent nature of particle bridges, our model would need to be extended to include stochasticity. However, our model can be interpreted as reflecting this process at the macroscale by relating the probability of bridging to the resistance of the suspension phase, which allows fluid to flow differentially through the Darcy phase.
Further additional physics could be incorporated into our model for detailed quantitative comparisons with experiments on synthetic and biological suspensions. For example, allowing particle slip on the walls would allow for the prediction of an unyielded but flowing state, as has been observed experimentally (Szafraniec et al. Reference Szafraniec, Valdez, Iffrig, Lam, Higgins, Pearce and Wood2022, Reference Szafraniec, Bull, Higgins, Stone, Krüger, Pearce and Wood2025). Improvements could also be made in capturing the rheology of deformable particles beyond our simple approximation of a change in the jamming fraction: particle deformability would also affect the ‘wet solid’ material properties, as well as the presence of a particle-free layer at the walls. More generally, relaxing the lubrication approximation would allow the model to describe pipes with sharper along-flow variations in their properties. Our model provides a step towards a broadly applicable continuum framework that links particle motion and macroscopic material properties during self-filtration and clogging of suspensions.
Acknowledgements
We would like to acknowledge numerous helpful discussions with V. Ciccone, A. Juel, H. Szafraniec, D. Wood, P. Bazazi, J. Higgins, T. Krüger and H. Stone. P.P. is supported by a UKRI Future Leaders Fellowship [MR/V022385/1].
All data are contained within the manuscript. Code to generate the figures is freely available at https://github.com/anushka07/Emergent_clogging.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Cross-sectional averaging
A.1. Radial integration
To solve (2.19) we need to compute various radial averages, which follow from suitable integration of the expressions for
$\phi$
(3.1) and
$u_s$
(3.3). We can isolate the dependence on the grouping
$\mu _w = GR/2p_s$
in these expression to write the relevant averages in terms of four dimensionless functions
\begin{align} \frac {\overline {\phi }}{\phi _m} &= f_1(\mu _w), & \frac {\overline { u_s\phi }}{G\phi _m R^2} &= f_2(\mu _w),\\ \end{align}
\begin{align} \frac {\overline {u_s}}{\textit{GR}^2} &= f_3(\mu _w), & \frac {\overline {u_D}}{Da \, G} &=f_4(\mu _w, \phi _m), \end{align}
where the functions
$f_i$
can be computed analytically by integrating the expressions in (3.1), (3.3) and the final equation in (2.16). The expressions are rather analytically involved, and are listed below for reference. Note that
$f_4$
is really just the radially averaged permeability function, coming from Darcy’s law (2.16), and thus depends in general on both
$\mu _w$
and
$\phi _m$
.
The average functions
$f_i$
, defined in (A1)–(A2), which depend on the grouping
$\mu _w= GR/2p_s$
, and represent the scaled average solid fraction, solid flux, solid velocity and Darcy velocity, respectively. For
$\mu _w \lt \mu _1$
, the entire pipe is clogged with
$\phi = \phi _m$
everywhere and
$u_s = 0$
. Here, we have taken
$Da = 10^{-1}$
and
$\phi _m = 0.8$
.
Figure 7 shows how these functions vary with
$\mu _w$
. For
$\mu _w\lt \mu _1$
, the entire pipe is clogged:
$\phi = \phi _m$
everywhere and
$u_s = 0$
; hence
$f_1 =1$
and
$f_2 = f_3 = 0$
over that range. As
$\mu _w$
is increased beyond this critical value, the material starts to flow, so
$f_1$
decreases and the other functions increase.
In terms of these functions, the governing equations (2.19) become
\begin{equation} \frac {\partial \phi _m f_1}{\partial t} + \frac {1}{R^2}\frac {\partial }{\partial z}\big(R^4\phi _m G f_2\big) = 0, \qquad \frac {Q}{\pi R^2} = R^2Gf_3 + Da \, G f_4. \end{equation}
The latter expression can be inverted to give an explicit expression for the pressure gradient
$G$
,
\begin{gather} G(\mu _w, R, \phi _m) = \frac {Q}{\pi R^4\big(f_3 + Da f_4/R^2\big)}, \end{gather}
and so (A3) becomes
\begin{equation} \frac {\partial }{\partial t}(\phi _m f_1) + \frac {1}{R^2}\frac {\partial \mathcal{F}}{\partial z} = 0, \qquad \mathcal{F}(\mu _w,\phi _m, R) \equiv \frac {\phi _m f_2 Q}{\pi \big(f_3 + Da f_4/R^2\big)}, \end{equation}
where
$\mathcal{F} = R^2\overline {u_s \phi }$
is the solid flux. The total pressure drop along the channel can be calculated from
\begin{gather} \Delta p = \int _{0}^{1} G \, {\rm d}z = Q/\pi \int _{0}^{1} \big(R^4\big(f_3 + Da f_4/R^2\big)\big)^{-1} {\rm d}z, \end{gather}
which also determines the overall resistance
$\mathcal{R}_o = \Delta p/Q = 1/\pi \int _0^1 (R^2(f_3 + Da f_4/R^2) )^{-1} {\rm d}z$
and the individual resistances of each phase
\begin{align}\mathcal{R}_s = \frac {\Delta p}{\pi R^2 \overline {u_s}} = \left (\frac {\Delta p}{Q} \right ) \frac {f_3 + Da f_4/R^2}{f_3}, \qquad \mathcal{R}_{\!D} = \frac {\Delta p}{\pi R^2 \overline {u_D}} = \left (\frac {\Delta p}{Q}\right ) \frac {f_3 + Da f_4/R^2}{Da f_4/R^2}.\end{align}
In the limit
$\overline {\phi } \to \phi _m$
, the particle velocity, and thus the function
$f_3$
, vanish, such that
$\mathcal{R}_s \to \infty$
and
$\mathcal{R}_{\!D} \to \mathcal{R}_o \sim Da^{-1}$
.
Note that the radius
$R$
enters the expression for the solid flux
$\mathcal{F}$
in (A5) only in the term multiplying the Darcy number. This feature illustrates how the reduction in the flux for smaller pipe radii, which controls the emergent ‘clogging’ behaviour that we observe here, is physically a result of the fact that the Darcy flow is relatively larger when the radius is reduced, as discussed in § 3.2. Relatively speaking, more of the total flux is taken up by differential seepage flow in a thinner pipe.
A.2. Integral functions
Explicit forms of the integral functions defined above are listed below for reference
\begin{gather} \begin{aligned} f_1 &=2(\mu _w\big(\,{-}\,3 \,{+}\, 2(\mu _w \,{-}\, \mu _1)^{1/2}\big) \,{+}\, 3\mu _1 \,{+}\, 2(\mu _w \,{-}\, \mu _1)^{1/2}(3 \,{+}\, 2\mu _1)\\[5pt] &\quad \,{-}\,\, 6(1 \,{+}\, \mu _1){\rm }\big(1 \,{+}\, (\mu _w \,{-}\, \mu _1)^{1/2}\big)/(3\mu _w^2) \,{+}\, (\mu _1/\mu _w)^2\\[5pt] f_2 &= (2\mu _1 \,{-}\, \mu _w)/\big(2\mu _w^2\big) \,{-}\, (\mu _w \,{-}\, \mu _1)^{1/2}/(7\mu _w) \,{-}\, \big(\mu _1^2 \,{-}\, 1\big)/\big(2\mu _w^3\big) \,{-}\, \mu _1/\big(2\mu _w^2\big)\\[5pt]& \quad \,{+}\, 1/(6\mu _w) \,{+}\, (\mu _1 \,{+}\, 1)/\big(4\mu _w^2\big) \,{-}\, \mu _1^3/\big(6\mu _w^4\big) \,{-}\, \big((\mu _w \,{-}\, \mu _1)^{1/2}\big(\,{-}\, 64\mu _1^3 \,{+}\, 140\mu _1^2\mu _w\\[5pt] & \quad \,{-}\, 119\mu _1^2\,{-}\, 70\mu _1\mu _w^2 \,{+}\, 210\mu _1\mu _w \,{+}\, 70\mu _1 \,{-}\, 105\mu _w^2 \,{+}\,105\big)\big)/\big(105\mu _w^4\big)\\[5pt] & \quad \,{+}\, \big((8\mu _1 \,{-}\, 7)(\mu _w \,{-}\, \mu _1)^{1/2}\big)/\big(35\mu _w^2\big)\,{-}\, \big(\mu _1\big(\,{-}\, \mu _1^2 \,{+}\, 2\mu _1\mu _w \,{-}\, \mu _w^2 \,{+}\, 1\big)\big)/\big(2\mu _w^4\big)\\[5pt] & \quad \,{+}\, \big(\mu _1^2(\mu _1 \,{-}\, 1)\big)/\big(4\mu _w^4\big)\,{+}\, \big((\mu _w \,{-}\, \mu _1)^{1/2}\big(32\mu _1^2 \,{-}\, 70\mu _1\mu _w \,{+}\, 7\mu _1 \,{+}\, 35\mu _w^2 \,{-}\, 35\big)\big)/\\[5pt] & \qquad \big(105\mu _w^3\big) \,{+}\, \big(\mu _1^2(\mu _1 \,{-}\, \mu _w)^2\big)/\big(4\mu _w^4\big)\\[5pt] & \quad\,{+}\, \big({\rm log}\big(\big(\mu _w \,{-}\, \mu _1 \,{+}\, 2(\mu _w \,{-}\, \mu _1)^{1/2} \,{+}\, 1\big)\big)(\mu _1 \,{+}\, 1)\\[5pt]&\qquad \,{\times}\,\big({-}\, \mu _1^2 \,{+}\, 2\mu _1\mu _w \,{-}\, \mu _w^2 \,{+}\, 1\big)\big)/\big(2\mu _w^4\big)\end{aligned} \end{gather}
\begin{gather} \begin{aligned} f_3 &= \big((\mu _1/\mu _w )^2(\mu _1/\mu _w - 1/2)\big)/2 - \mu _1/6\mu _w + \big((\mu _1/\mu _w)^2(\mu _1/\mu _w - 1)^2\big)/4\\[4pt] & \quad - \big(5(\mu _1/\mu _w )^4\big)/24+ 1/8\\[4pt] f_4 &=-\big(16\mu _1^2(\mu _w - \mu _1)^{1/2} - 24\mu _w^2(\mu _w - \mu _1)^{1/2} + 60\phi _m^3(\mu _w - \mu _1)^{1/2}\\[4pt] & \quad + 15\mu _1\mu _w^2 + 30\mu _1\phi _m^3 + 45\mu _w^2\phi _m - 30\mu _w\phi _m^3 - 5\mu _1^3 - 15\mu _w^2 - 10\mu _w^3\\[4pt]& \quad - 30\phi _m^3{\rm log}\big(\big(\mu _w - \mu _1 + 2(\mu _w - \mu _1)^{1/2} + 1\big)\big)+ 15\mu _1^2\phi _m^3 - 45\mu _w^2\phi _m^2\\[4pt] & \quad + 8\mu _1\mu _w(\mu _w - \mu _1)^{1/2} - 30\mu _1\phi _m^3{\rm log}\big(\big(\mu _w - \mu _1 + 2(\mu _w - \mu _1)^{1/2} + 1\big)\big)\\[4pt] & \quad - 24\mu _1^2\phi _m(\mu _w - \mu _1)^{1/2}+ 40\mu _1\phi _m^3(\mu _w - \mu _1)^{1/2} + 36\mu _w^2\phi _m(\mu _w - \mu _1)^{1/2}\\[4pt] & \quad + 20\mu _w\phi _m^3(\mu _w - \mu _1)^{1/2}- 12\mu _1 \mu _w\phi _m(\mu _w - \mu _1)^{1/2}\big)/\big(15\mu _w^2\phi _m^2\big).\end{aligned} \end{gather}
Appendix B. Time-dependent numerical method
Equation (2.19) is a quasi-linear hyperbolic partial differential equation which can give rise to shocks for certain initial conditions. To handle these, we use the finite-volume package Clawpack (Clawpack Development Team 2024). Clawpack uses a high-resolution Godunov-type method, which employs Riemann solvers at cell interfaces to capture the propagation of waves accurately. We use second-order Godunov method with a van Leer flux limiter to solve the Riemann problems for the whole domain.
For the boundary conditions, we use zero-order extrapolation by implementing ghost cells that extend the cell values from the upstream and downstream boundaries; these are updated at every time step. This ensures that no spurious waves are generated at the boundaries that might alter the solution. The system of equations that we input in Clawpack are in conservative form to ensure accuracy. Using definitions of integral functions from (A2), we define a new variable
$q$
for computational purposes
\begin{align} q &= R^2\phi _m f_1(\mu _w). \end{align}
We encode a spatially varying radius
$R(z$
) or maximum packing fraction
$\phi _m(z)$
through two extra equations in combination with the explicit form of the conservation equation (2.19)
\begin{equation} \frac {\partial q}{\partial t} + \frac {\partial }{\partial z}\left (\frac {\phi _m f_2}{f_3 + Da f_4/R^2 }\right )=0, \qquad \frac {\partial R}{\partial t} = 0, \qquad \frac {\partial \phi _m}{\partial t} = 0. \end{equation}
The specific form of
$R(z)$
or
$\phi _m(z)$
is then set as an ‘initial’ condition for that variable. The form for radial constrictions was given in (3.4), and for simulations with varying jamming fraction we set
$ \phi _m(z,t=0) = a - (a-b)\left \{1 + \textit{erf}\,[(z-0.5)10]\right \}/2$
, where
$a$
and
$b$
represent the upstream and downstream values.
Steady results for fixed pressure drop, showing how the final steady-state total flux
$Q$
decreases with the initial inlet solid fraction, for (a) varying
$Da$
with
$R_{min}= 0.5$
and (b) varying constriction ratio
$R_{min}$
with
$Da = 10^{-2}$
. Here,
$\phi _m = 0.65$
.
Appendix C. Fixed pressure drop
The time-dependent numerical method would become significantly more numerically expensive if the pressure drop along the pipe, as opposed to the total flux, was fixed (this would require iteration of an integral constraint on the flux at each time step). However, the steady solutions can be straightforwardly inverted to give the case with a fixed pressure drop: results follow directly from the overall resistance in figure 5(b). The resistance provides a measure of the pressure drop required to sustain a given total flux, and so its inverse provides the flux required to sustain a given pressure drop in steady state. This quantity is plotted in figure 8, for various different parameters, which shows how the total flux drops discontinuously upon entering the ‘clogged’ regime.
Open access
