1 Introduction
One of the interesting features of a yield-stress fluid is that it can hold particles/bubbles of small buoyancy statically suspended. This is practically important for many industries: from producing stable fresh concrete in the construction industry (Roussel Reference Roussel2006) to preventing sedimentation of drilling/rock cuttings in the oil and gas industry (Elgaddafi et al. Reference Elgaddafi, Ahmed, George and Growcock2012). Moreover, transporting suspensions of yield-stress fluids (e.g. pumping concrete (Choi, Kim & Kwon Reference Choi, Kim and Kwon2013) and fracturing flows (Hormozi & Frigaard Reference Hormozi and Frigaard2017)) is another crucial problem to be understood for these applications because of multiple intrinsic complexities from carrying fluid rheology to particle dynamics. Valiant efforts of Segré & Silberberg (Reference Segré and Silberberg1961, Reference Segré and Silberberg1962a,Reference Segré and Silberbergb) and other scholars uncovered the complex behaviour of particles suspended in a Newtonian fluid Poiseuille flow. Nevertheless, for yield-stress fluids, the Poiseuille flow itself is split into two main regions: the sheared yielded regions close to the walls and the core unyielded region. The complexity of particle dynamics in the yielded regions comes from various sources such as the nonlinearity of the fluid behaviour and particle–particle and particle–wall interactions, whereas in the core unyielded region particles might be trapped (which was evidenced experimentally by Merkak et al. (Reference Merkak, Jossic and Magnin2008, Reference Merkak, Jossic and Magnin2009) and Zade et al. (Reference Zade, Shamu, Lundell and Brandt2020)) because of the small stresses compared to the yield stress.
In the present study, we will shed light on the particle behaviour inside the core unyielded region which is important for understanding the homogeneity of the suspension flow. Indeed, we systematically address the stability of particles suspended inside yield-stress fluid Poiseuille flows. The yield number is the most important non-dimensional number for this type of problem. For instance, in studies of static stability of non-neutrally buoyant particles, the yield number (
$Y$) takes the form 
$\hat{\unicode[STIX]{x1D70F}}_{Y}/\unicode[STIX]{x0394}\hat{\unicode[STIX]{x1D70C}}{\hat{g}}\hat{\ell }$, where it represents the ratio of the yield stress (
$\hat{\unicode[STIX]{x1D70F}}_{Y}$) of the fluid to the buoyancy stress (
$\unicode[STIX]{x0394}\hat{\unicode[STIX]{x1D70C}}{\hat{g}}\hat{\ell }$), with 
$\unicode[STIX]{x0394}\hat{\unicode[STIX]{x1D70C}}$ designating the difference between the particle density (
$\hat{\unicode[STIX]{x1D70C}}_{p}$) and the fluid density (
$\hat{\unicode[STIX]{x1D70C}}_{f}$), 
${\hat{g}}$ the gravitational acceleration and 
$\hat{\ell }$ the characteristic length scale of the particle. We review the existing literature concerning the particle stability in yield-stress fluids within different settings in what follows.
1.1 Stability of particles inside a quiescent bath of a yield-stress fluid
The study of particle motion in yield-stress fluids dates back to the 1950s when Volarovich & Gutkin (Reference Volarovich and Gutkin1953), in the very first attempts in this context, found that if a particle moves within a viscoplastic medium then it should do so within a bounded region, since the stress falls below the yield stress sufficiently far from the particle. Almost a decade later (in the early 1960s), Whitmore and colleagues started to look at the same problem in a series of experiments and tried to measure the limiting force (
$\hat{F}_{D,c}$) that should be exceeded to move a particle in a yield-stress fluid (i.e. overcome the yield-stress resistance). To show that 
$\hat{F}_{D,c}$ is not simply equal to the yield stress times the surface area of the particle (a ‘folk-tale’ at that time), Boardman & Whitmore (Reference Boardman and Whitmore1960) measured 
$\hat{F}_{D,c}$ for different orientations of particles and demonstrated that the values are different. Five years later, Valentik & Whitmore (Reference Valentik and Whitmore1965) experimentally observed that there is ‘an envelope of suspension [unyielded material] attaching itself to a moving sphere’. Many interesting attempts by other scholars partly uncovered different features of the problem (e.g. Ansley & Smith Reference Ansley and Smith1967; Yoshioka, Adachi & Ishimura Reference Yoshioka, Adachi and Ishimura1971; Adachi & Yoshioka Reference Adachi and Yoshioka1973), until the remarkable numerical study of Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985) fully revealed the flow pattern about a moving sphere in a yield-stress fluid. Two caps of unyielded fluid were found to be present at the leading and trailing parts of the sphere (close to the stagnation points where the stress goes below the yield stress). They also calculated 
$\hat{F}_{D,c}$ for a sphere. Following Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985), many scholars have worked on single-particle motion in yield-stress fluids, especially on finding the limiting force (
$\hat{F}_{D,c}$) which was an open question since the time of Boardman & Whitmore (Reference Boardman and Whitmore1960). Many examples can be found in the literature for different particle shapes using a resistance formulation (i.e. imposing velocity on the particle and calculating the drag force; see e.g. Tokpavi, Magnin & Jay (Reference Tokpavi, Magnin and Jay2008) and Nirmalkar, Chhabra & Poole (Reference Nirmalkar, Chhabra and Poole2012)). Measuring the limiting force also has been the objective of many experimental studies. Examples can be found in the works of Magnin and co-workers (Tokpavi et al. Reference Tokpavi, Jay, Magnin and Jossic2009).
Putz & Frigaard (Reference Putz and Frigaard2010) formed a mathematical framework for studying the stability of particles by considering the mobility problem (i.e. considering the buoyancy of the particle and seeing if the particle settles or is static):
$$\begin{eqnarray}\displaystyle \left\{\begin{array}{@{}l@{}}\text{Particle settles}\quad \text{iff }Y<Y_{c}^{\ast }\text{(large buoyancy compared to the yield stress)},\\ \text{Particle is static}\quad \text{iff }Y>Y_{c}^{\ast }\text{(large yield stress compared to the buoyancy)},\end{array}\right.\quad & & \displaystyle\end{eqnarray}$$ where 
$Y_{c}^{\ast }$ is the critical value of the yield number which decides the stability of the particle (see § 2.3). Chaparian & Frigaard (Reference Chaparian and Frigaard2017b) followed the same framework and also revealed the relevance of perfect-plasticity theories in the study of yield limit or particle stability. Exploiting the cloaking phenomenon (Chaparian & Frigaard Reference Chaparian and Frigaard2017a), they systematically presented a model for calculating 
$Y_{c}^{\ast }$ for a single symmetric particle by finding the unyielded envelope enclosing the particle and postulating an admissible stress/velocity field about the particle. Chaparian (Reference Chaparian2018) also explored the effect of particle orientation on the yielding behaviour. An assembly of particles could exhibit different yielding characteristics depending on their proximity which was investigated recently by Chaparian, Wachs & Frigaard (Reference Chaparian, Wachs and Frigaard2018) and Koblitz, Lovett & Nikiforakis (Reference Koblitz, Lovett and Nikiforakis2018): a cluster of particles will be formed by connector bridges of unyielded material.
1.2 Particle migration and stability in a Poiseuille flow
In contrast to the particle stability inside a quiescent bath of yield-stress fluids, the problem of entrapment of particles inside an unyielded plug and particle migration within a Poiseuille flow is less understood. This problem was investigated experimentally by Merkak et al. (Reference Merkak, Jossic and Magnin2008, Reference Merkak, Jossic and Magnin2009) in inertialess pipe flows (maximum particle Reynolds number of 
$O(10^{-1})$). The suspensions used were highly stabilized and at rest the spheres could not settle under gravity alone: 
$Y\gg Y_{c}^{\ast }$. The main findings can be highlighted as follows:
(i) Particles situated in the plug zone will move with a translational velocity which is independent of their position in the plug zone. No rotation or migration was reported for this case.
(ii) Particles situated completely or partly in the sheared zones (yielded zones) will translate and also rotate (i.e. migrate) with velocities that are functions of the yield stress and the pressure drop (i.e. the Bingham number).
These observations have been reported for neutrally buoyant particles. For the case in which particles are denser than the fluid (i.e. non-neutrally buoyant particles), particles sediment in the sheared zone by bypassing the plug zone and reach their final position at the bottom of the pipe. This is intuitive as a sheared yield-stress fluid cannot hold non-neutrally buoyant particles static (Ovarlez et al. Reference Ovarlez, Bertrand, Coussot and Chateau2012). In the plug zone, however, different scenarios can happen:
(a) Non-neutrally buoyant particles with large
$Y$ are stable in the core unyielded zone and just translate with a same velocity as the plug zone. This is the same as the phenomenon observed for neutrally buoyant particles.(b) Particle sedimentation can happen for smaller values, yet larger than
$Y_{c}^{\ast }$, of the yield number.
In the present study, we analyse this problem numerically to avoid the inevitable constraints in experiments. Hence, the main objective of the current study is to systematically investigate the conditions under which particles are stable in the core unyielded plug: calculating 
$Y_{c}$ as a function of different effective parameters such as the confinement ratio and the pressure drop. If 
$Y>Y_{c}$ then a particle is stable inside the plug region and otherwise (i.e. 
$Y<Y_{c}$) will break the plug and start to sediment/migrate. Moreover, we address the relevance of 
$Y_{c}$ to 
$Y_{c}^{\ast }$.
Contradictory to the stable particles in the core unyielded region which have been observed experimentally by Merkak et al. (Reference Merkak, Jossic and Magnin2008, Reference Merkak, Jossic and Magnin2009), very recently Siqueira & de Souza Mendes (Reference Siqueira and de Souza Mendes2019) studied the particle migration in an ‘apparent’ yield-stress fluid and reported that eventually (given infinite time) the plug zone is evacuated from the particles, which does not seem intuitive. It may be a consequence of using the regularized effective viscosity of a viscoplastic material with the Krieger–Dougherty expression for the suspension viscosity and a diffusive flux model for the shear-induced particle migration. This contradiction markedly signals the need for direct numerical simulations using unregularized rheology, which is the approach in the present study.
An outline of the paper is as follows. In § 2, we set out the studied problem and briefly review the simple analytical solutions of the basic flow and some features of the implemented numerical method. Then, § 3 presents the main results. It includes firstly a demonstration of the negligible effect of confined geometry on its own. Then follows a study of non-neutrally buoyant particles hosted by a viscoplastic Poiseuille flow, addressing the yield limit, and also possible scalings in this limit. The behaviour of neutrally buoyant particles is considered next and some mathematical hints for finding the furthest stable particles are provided. The final section (§ 4) concludes the present study with a brief summary and some comments on the applicability of the results.
2 Problem statement
An inertialess viscoplastic Poiseuille flow is considered between two infinite parallel plates. We consider a rigid non-neutrally buoyant particle (two-dimensional circle) of diameter 
$\hat{D}$ (and radius 
$\hat{R}$) inside an infinitely long channel of width 
${\hat{H}}$ (see figure 1). The centre of the particle is positioned at 
${\hat{y}}^{p}$ where the 
$x$ axis is aligned and put at the axial centreline of the channel.
2.1 Scalings and remarks
Three important non-dimensional parameters for the considered problem are the Bingham number, the yield number and the confinement ratio, which take the forms
$$\begin{eqnarray}\displaystyle B=\frac{\hat{\unicode[STIX]{x1D70F}}_{Y}\hat{D}}{\hat{\unicode[STIX]{x1D707}}\hat{U} },\quad Y=\frac{\hat{\unicode[STIX]{x1D70F}}_{Y}}{(\hat{\unicode[STIX]{x1D70C}}_{p}-\hat{\unicode[STIX]{x1D70C}}_{f}){\hat{g}}\hat{D}}\quad \text{and}\quad \unicode[STIX]{x1D709}=\frac{{\hat{H}}}{\hat{D}}, & & \displaystyle\end{eqnarray}$$ respectively. The confinement ratio represents the ratio of the channel width (
${\hat{H}}$) to the particle diameter. The Bingham number and yield number show the ratio of the yield stress to the characteristic viscous stress (
$\hat{\unicode[STIX]{x1D707}}\hat{U} /\hat{D}$) and buoyancy stress, respectively. Note that here the particle diameter is used as the length scale to facilitate any direct comparison with the stability of particles in a quiescent background fluid where no other characteristic length exists (e.g. 
${\hat{H}}$ in the current study).
On scaling the velocity vector (
$\hat{\boldsymbol{u}}=(\hat{u} ,\hat{v})$) with the mean velocity at the inlet (
$\hat{U}$), pressure (
$\hat{p}$) and deviatoric stress tensor (
$\hat{\unicode[STIX]{x1D749}}$) with the characteristic viscous stress (
$\hat{\unicode[STIX]{x1D707}}\hat{U} /\hat{D}$), we find the dimensionless Stokes equation and Bingham constitutive equation as
$$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}+G\boldsymbol{e}_{x}-\frac{B}{Y}\frac{\unicode[STIX]{x1D70C}}{1-\unicode[STIX]{x1D70C}}\boldsymbol{e}_{y}=0 & & \displaystyle\end{eqnarray}$$and
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}\unicode[STIX]{x1D749}=\left(1+\displaystyle \frac{B}{\Vert \dot{\unicode[STIX]{x1D738}}\Vert }\right)\dot{\unicode[STIX]{x1D738}}\quad \text{iff }\Vert \unicode[STIX]{x1D749}\Vert >B,\\ \dot{\unicode[STIX]{x1D738}}=0\quad \text{iff }\Vert \unicode[STIX]{x1D749}\Vert \leqslant B,\end{array}\right\} & & \displaystyle\end{eqnarray}$$ respectively, where 
$G$ is the magnitude of the pressure gradient in the axial direction which is required to provide the flow rate of 
$\hat{U} {\hat{H}}$. It can be also absorbed in the pressure and form 
$p^{\prime }$. The gravitational acceleration is aligned with 
$-\boldsymbol{e}_{y}$. The density ratio is represented by 
$\unicode[STIX]{x1D70C}=\hat{\unicode[STIX]{x1D70C}}_{f}/\hat{\unicode[STIX]{x1D70C}}_{p}$ and the rate of strain tensor by 
$\dot{\unicode[STIX]{x1D738}}$. Note that 
$B/Y=(\hat{\unicode[STIX]{x1D70C}}_{p}-\hat{\unicode[STIX]{x1D70C}}_{f}){\hat{g}}\hat{D}^{2}/\hat{\unicode[STIX]{x1D707}}\hat{U}$ is termed the buoyancy number in Newtonian fluid mechanics, and represents the ratio of the particle Stokes settling velocity to the mean channel velocity.
Schematic of the problem.
2.2 Undisturbed (particle-free) background flow
Here we review the solutions to the velocity and stress fields of the undisturbed background flow (i.e. in the absence of any particle) which are useful for the rest of the present study. In this short subsection, we use 
${\hat{H}}$ as the only length scale in the absence of any particle and designate the new variables with a tilde 
$(\tilde{\cdot })$. Hence, two new non-dimensional quantities will be the slot Bingham number, 
$\tilde{B}=\hat{\unicode[STIX]{x1D70F}}_{Y}{\hat{H}}/\hat{\unicode[STIX]{x1D707}}\hat{U}$, and 
$\tilde{G}$ which is the absolute value of the pressure gradient scaled with 
$\hat{\unicode[STIX]{x1D707}}\hat{U} /{\hat{H}}^{2}$. Then, the new variables can be easily converted to the ones defined earlier based on the particle diameter as 
${\tilde{y}}=y/\unicode[STIX]{x1D709}$, 
$\tilde{\unicode[STIX]{x1D749}}=\unicode[STIX]{x1D749}\unicode[STIX]{x1D709}$, 
$\tilde{B}=B\unicode[STIX]{x1D709}$ and 
$\tilde{G}=G\unicode[STIX]{x1D709}^{2}$.
The mean velocity at the inlet is 
$\hat{U} =(\int _{-{\hat{H}}/2}^{{\hat{H}}/2}\hat{u} \,\text{d}{\hat{y}})/{\hat{H}}$. Hence, the continuity equation in non-dimensional form is written as
$$\begin{eqnarray}\displaystyle 1=2\int _{0}^{1/2}\tilde{u} \,\text{d}{\tilde{y}}, & & \displaystyle\end{eqnarray}$$ and 
$\tilde{G}$ will satisfy
$$\begin{eqnarray}\displaystyle 1=\frac{\tilde{G}}{12}-\tilde{B}\left(\frac{1}{4}-\frac{\tilde{B}^{2}}{3\tilde{G}^{2}}\right), & & \displaystyle\end{eqnarray}$$ and 
$\tilde{\unicode[STIX]{x1D70F}}_{xy}=\tilde{G}{\tilde{y}}$. Note that, as a consequence of scaling the velocity vector with the mean velocity at the inlet, we always have a non-zero background flow even in the limit 
$\tilde{B}\rightarrow \infty$ (or 
$B\rightarrow \infty$ for a fixed 
$\unicode[STIX]{x1D709}$). Indeed, in the limit 
$\tilde{B}\rightarrow \infty$, almost the whole gap width of the channel moves as a plug with approximately unity velocity except for the two highly sheared boundary layers close to the walls. Nevertheless, in this limit, there is an intrinsic high shear stress (i.e. 
$\tilde{G}\rightarrow \infty$) in the flow which makes it totally different from the quiescent fluid which is unyielded.
The velocity profile (
$\tilde{\boldsymbol{u}}=(\tilde{u} ,0)$) of the undisturbed flow can be derived as
$$\begin{eqnarray}\displaystyle \tilde{u} ({\tilde{y}})=\left\{\begin{array}{@{}l@{}}\tilde{B}\left(|{\tilde{y}}|-\displaystyle \frac{1}{2}\right)+\displaystyle \frac{\tilde{G}}{2}\left(\frac{1}{4}-{\tilde{y}}^{2}\right)\quad \text{iff }{\tilde{y}}^{\ast }<|{\tilde{y}}|\leqslant \displaystyle \frac{1}{2},\\ \displaystyle \frac{1}{2}\frac{\tilde{B}^{2}}{\tilde{G}}+\frac{\tilde{G}}{8}-\frac{\tilde{B}}{2}\quad \text{iff }0\leqslant |{\tilde{y}}|\leqslant {\tilde{y}}^{\ast }.\end{array}\right. & & \displaystyle\end{eqnarray}$$ The position of the undisturbed yield surface is denoted as 
$\pm {\tilde{y}}^{\ast }$ which is equal to 
$\pm \tilde{B}/\tilde{G}$.
2.3 Stability of particles in the sedimentation problem
In the context of the stability of particles inside a quiescent bath of a yield-stress fluid, the limit 
$Y\rightarrow Y_{c}^{\ast }$ is often studied through a steady Stokes flow and has the meaning of a yield number limit above which the steady velocity solution is zero (no motion). This is sometimes referred to as static stability. The mathematical definition of 
$Y_{c}^{\ast }$ could be extracted from the energy equation (Putz & Frigaard Reference Putz and Frigaard2010; Chaparian & Frigaard Reference Chaparian and Frigaard2017b) as
$$\begin{eqnarray}\displaystyle Y_{c}^{\ast }=\sup _{\boldsymbol{v}\in {\mathcal{V}},\boldsymbol{v}\not =0}\frac{L(\boldsymbol{v})}{j(\boldsymbol{v})}=\sup _{\boldsymbol{v}\in {\mathcal{V}},\boldsymbol{v}\not =0}\frac{-\displaystyle \int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{e}_{y}\,\text{d}A}{\displaystyle \int _{\unicode[STIX]{x1D6FA}\setminus \bar{X}}\Vert \dot{\unicode[STIX]{x1D738}}\Vert \,\text{d}A}, & & \displaystyle\end{eqnarray}$$ where 
$\boldsymbol{v}$ is a velocity test function from the set of admissible velocity fields 
${\mathcal{V}}$, 
$L(\boldsymbol{v})$ is the work done by the body force (in this case by gravity) and 
$j(\boldsymbol{v})$ is the plastic dissipation over 
$B$. For a single particle inside a quiescent bath of a yield-stress fluid, 
$Y>Y_{c}^{\ast }$ also corresponds to a dynamic stability limit above which (even inertial) particle motion decays to zero in finite time as shown both theoretically and computationally by Wachs & Frigaard (Reference Wachs and Frigaard2016).
For a single two-dimensional disc, 
$Y_{c}^{\ast }\approx 0.0658$ has been reported in the literature (Tokpavi et al. Reference Tokpavi, Magnin and Jay2008; Putz & Frigaard Reference Putz and Frigaard2010; Chaparian & Frigaard Reference Chaparian and Frigaard2017b). Note that in the present study we have chosen the diameter of the particle as the length scale. If the radius is chosen (the same as the mentioned references) then the value of 
$Y_{c}^{\ast }$ will be doubled.
2.4 Numerical details
We use a finite element/augmented Lagrangian approach (Roquet & Saramito Reference Roquet and Saramito2003) to attack equations (2.2) and (2.3). The background Poiseuille flow is solved using the same strategy as that of Roustaei & Frigaard (Reference Roustaei and Frigaard2015). Indeed, the desired pressure gradient which makes the steady base flow (Poiseuille flow) is found in a semi-iterative way. The coupling of the particle with the base flow is ensured using the same method as utilized by Putz & Frigaard (Reference Putz and Frigaard2010) and Chaparian & Frigaard (Reference Chaparian and Frigaard2017b). We avoid repeating the details here. For validation studies and details of the mesh adaptation, the reader is refereed to the appendix of Chaparian & Frigaard (Reference Chaparian and Frigaard2017b) and Chaparian et al. (Reference Chaparian, Izbassarov, De Vita, Brandt and Tammisola2019). Regarding the size of the computational domain in the 
$x$ direction, we add that it has been chosen large enough to avoid altering the flow field about and in the vicinity of the particle (e.g. in some cases as large as 60 times the particle diameter). However, in the following figures only a small box close to the particle is shown in post-processing. The implementation has been conducted in an open-source C++ finite element environment, FreeFEM
$++$ (Hecht Reference Hecht2012).
3 Results
In this section, we present the new findings on particle stability in two-dimensional channel flow. It consists of computing 
$Y_{c}$ in the present problem: if 
$Y>Y_{c}$ then the particle is stable inside the plug region and otherwise (i.e. 
$Y<Y_{c}$) it will break the plug and start to sediment/migrate. Different effects are analysed individually. First, we outline the influence of confinement (walls) alone (i.e. in the absence of any background flow) on the stability of the non-neutrally buoyant particles in § 3.1. Then we move further to consider the particles hosted in a Poiseuille flow (§ 3.2). Finally, we decouple the buoyancy effects from the problem in § 3.3. In all the cases, we focus on the situation where the particle losses its stability and study the limiting yield number as well as different behaviours that occur in this limit.
3.1 Confinement
In this subsection, the effect of geometry confinement on the static stability of a particle is investigated. In other words, we try to determine if the confinement, independently of the flow field, has any effect on the particle static stability in a quiescent yield-stress fluid. We use the same scalings as in the rest of the paper, but here there is no background pressure-driven Poiseuille flow. Indeed, 
$G=0$ and the Bingham number does not have any special interpretation in the absence of characteristic viscous stress associated with the background flow. Therefore, in this case, it is only the yield number that decides the particle stability.
Figure 2(a–c) shows the flow fields near the yield limit (
$Y\approx Y_{c}$) for different confinement ratios. Figure 2(d), however, shows the sliplines of the perfectly plastic problem (see Chaparian & Frigaard (Reference Chaparian and Frigaard2017b) for details of the slipline solution) which is relevant to the yield limit in a viscoplastic problem; i.e. when the viscous dissipation is at least one order of magnitude smaller than the plastic dissipation (Putz & Frigaard Reference Putz and Frigaard2010). As can be observed, no significant change is distinguishable in the flow fields, although for the case 
$\unicode[STIX]{x1D709}=2$, the unyielded kidney close to the particle is smaller which manifests small deviation from the velocity field associated with infinite domain (figure 2c).
Velocity magnitude contours at 
$Y=0.0655$: (a) 
$\unicode[STIX]{x1D709}=2$, (b) 
$\unicode[STIX]{x1D709}=5$ and (c) infinite domain (
$\unicode[STIX]{x1D709}\rightarrow \infty$). Note that one quarter of the domain is presented in each panel. (d) The sliplines associated with the perfectly plastic problem when 
$\unicode[STIX]{x1D709}\rightarrow \infty$. In (a–c), the white lines represent the yield surfaces.
Velocity magnitude contours near the yield limit: (a) 
$\unicode[STIX]{x1D709}=1.05$ and 
$Y=0.04$, (b) 
$\unicode[STIX]{x1D709}=1.2$ and 
$Y=0.055$ and (c) 
$\unicode[STIX]{x1D709}=1.5$ and 
$Y=0.065$. Note that one quarter of the domain is presented in each panel. (d) The same as figure 2(d). In (a–c), the white lines represent the yield surfaces. (e) Critical yield number 
$Y_{c}$ versus 
$\unicode[STIX]{x1D709}$ (the circular symbols). The brown solid line marks 
$Y_{c}^{\ast }$ for reference.
Causing negligible differences in the velocity fields, confinement by itself does not change 
$Y_{c}$ of the particle in the range of confinement ratios larger than or equal to 2: for 
$2\leqslant \unicode[STIX]{x1D709}$, 
$Y_{c}=Y_{c}^{\ast }=0.0658$. However, on decreasing 
$\unicode[STIX]{x1D709}$ to smaller values, the channel walls interact considerably with the flow field about the particle (see figure 3). Figure 3(e) clearly demonstrates that 
$Y_{c}$ decreases on decreasing 
$\unicode[STIX]{x1D709}$. This could be explained by considering that by decreasing 
$\unicode[STIX]{x1D709}$, the walls prevent particle sedimentation by resisting against the flow generated by the settling particles. Hence, particles should have larger buoyancy to overcome this hindering effect of the nearby walls.
In what follows, we only consider the case 
$2\leqslant \unicode[STIX]{x1D709}$, in which the effect of the confined geometry by itself is negligible on the particle stability.
3.2 Non-neutrally buoyant particles in Poiseuille flow
We will consider the stability of particles inside a Poiseuille flow in this subsection. For the Stokes flow, the problem is symmetric in the vertical direction (despite buoyancy); therefore, only particles in the top half of the channel will be considered (
$y_{p}\geqslant 0$). In what follows, for a particle at a given height, 
$\bar{B}$ designates the value of the Bingham number at which the undisturbed plug thickness is just sufficient to contain the whole particle inside the plug. In other words, when 
$B=\bar{B}$, then 
$y^{\ast }=y_{p}+R=y_{p}+\frac{1}{2}$. Hence, 
$\bar{B}$ is a value corresponding to the particle position, not a flow characteristic.
3.2.1 The case $y_{p}=0$
We consider the particles that are positioned at the centreline of the channel in this section. For instance, two sample simulations are illustrated in figure 4, for the case 
$\unicode[STIX]{x1D709}=5$ and 
$B=2$ and two yield numbers, one larger than 
$Y_{c}$ and the other one less than that.
Velocity magnitude contours for 
$\unicode[STIX]{x1D709}=5,B=2$: (a) 
$Y=0.75>Y_{c}$ and (b) 
$Y=0.3<Y_{c}$. The solid white lines show the yield surfaces, while the dashed white lines show the position of the yield surfaces in the absence of any particle (i.e. 
$y^{\ast }$ in single-phase flow).
From figure 4, it is clear that for this case 
$Y_{c}\gg Y_{c}^{\ast }$. Indeed, the shear stress which is generated by the background flow helps the buoyancy stress to overcome the yield stress resistance and breaks the plug. Figure 5 presents 
$Y_{c}$ versus the Bingham number for a wide range of confinement ratios. It clearly shows that by increasing the Bingham number, given that 
$\unicode[STIX]{x1D709}$ is constant, 
$Y_{c}$ decreases. This behaviour is not observed for 
$\unicode[STIX]{x1D709}=2$ where almost the same 
$Y_{c}$ value is computed for the range of Bingham numbers that is considered (
$B>\bar{B}$).
Critical yield number versus Bingham number for various confinement ratios. Dashed vertical lines show the values of 
$\bar{B}$ for different confinement ratios with the same colour interpretation. The brown solid line shows the value of 
$Y_{c}^{\ast }$. The values of 
$Y_{c}$ at large Bingham numbers are marked with 
$Y_{c}^{\infty }$.
In what follows in this subsection, we aim to describe the stabilizing effect of increasing Bingham number (figure 5) and also to find a simple model for the limiting value of 
$Y_{c}(B\rightarrow \infty )=Y_{c}^{\infty }$ based on the lucid solutions to the undisturbed flow (§ 2.2).
There are two contradictory effects regarding the particle stability when the Bingham number increases: (i) the plug thickness grows and (ii) the background shear stress will be larger. It is intuitive that the larger the plug thickness compared to the particle diameter, the more stable is the particle – resulting in a smaller 
$Y_{c}$. However, by increasing the Bingham number, the shear stress of the background flow will increase as well (see figure 6a) which has an opposite effect: the background shear stress helps the particle buoyancy stress to overcome the yield stress resistance (resulting in increasing 
$Y_{c}$). Figure 5 clearly demonstrates that the first effect is dominant except for the case 
$\unicode[STIX]{x1D709}=2$. This can be roughly understood by considering different stresses which contribute at the yield limit: one can approximately assume that
$$\begin{eqnarray}\displaystyle \text{Shear stress contribution `}+\text{' buoyancy stress }\sim \text{yield stress} & & \displaystyle\end{eqnarray}$$ at the yield limit. We point out that here the plus sign is not a ‘mathematical’ operator since the stress cannot be linearly superimposed. Hence, at the yield limit (remembering that 
${\hat{G}}$ is the pressure gradient and the shear stress is proportional to that), 
${\hat{G}}\hat{D}\text{ `}+\text{' }\unicode[STIX]{x0394}\hat{\unicode[STIX]{x1D70C}}{\hat{g}}\hat{D}\sim \hat{\unicode[STIX]{x1D70F}}_{Y}$, or in non-dimensional form
$$\begin{eqnarray}\displaystyle G\text{ `}+\text{' }\frac{B}{Y_{c}}\sim B. & & \displaystyle\end{eqnarray}$$ We emphasize that as a ‘naive’ model, in the above equation we use the pressure gradient of the undisturbed flow, expression (2.5). The above equation shows that if the right-hand side increases and the increase in 
$G$ as a function of Bingham number is not large enough to balance the right-hand side, then 
$Y_{c}$ decreases. Figure 6(a) illustrates 
$G$ as a function of the Bingham number. It reveals that the increase in 
$G$ is almost linear with 
$B$ and only for the case 
$\unicode[STIX]{x1D709}=2$ is the slope (i.e. 
$a$) larger than unity. Hence, for this case no decrease in 
$Y_{c}$ is observable with increasing Bingham number since the background shear stress exceedingly contributes in breaking the yield stress resistance.
Exploiting the approximate linearity of 
$G$ over this range of the Bingham number, 
$G\sim aB$ (see figure 6a), and knowing that as 
$B\rightarrow \infty$ then 
$Y_{c}$ reaches its limiting value (i.e. 
$Y_{c}\rightarrow Y_{c}^{\infty }$), from equation (3.2) we can extract that 
$Y_{c}^{\infty }$ is a function of 
$a$. From figure 6(b), we can clearly deduce that 
$Y_{c}^{\infty }\approx 1.1a$, and therefore
$$\begin{eqnarray}\displaystyle G\sim aB\approx \frac{Y_{c}^{\infty }}{1.1}B. & & \displaystyle\end{eqnarray}$$ Using the fact that 
$y^{\ast }=B/G$, we rewrite the above expression as
$$\begin{eqnarray}\displaystyle Y_{c}^{\infty }\approx \frac{1.1}{y^{\ast ,\infty }} & & \displaystyle\end{eqnarray}$$ as a ‘naive’ model for calculating 
$Y_{c}^{\infty }$, where 
$y^{\ast ,\infty }$ shows the position of the yield surface at the limit 
$B\rightarrow \infty$ (i.e. approximately half of the thickness of the plug in the limit 
$B\rightarrow \infty$). Figure 6(c) compares this model with the numerical data where 
$\unicode[STIX]{x1D709}/2$ is used for 
$y^{\ast ,\infty }$: indeed when 
$B\rightarrow \infty$, then 
$y^{\ast }\rightarrow \unicode[STIX]{x1D709}/2$.
(a) Pressure gradient, 
$G$, of the undisturbed flow as a function of the Bingham number for different confinement ratios; equation (2.5). (b) Critical yield number 
$Y_{c}^{\infty }$ versus 
$a$. (c) Critical yield number 
$Y_{c}^{\infty }$ versus 
$\unicode[STIX]{x1D709}$; the blue solid line is the fitted curve on the numerical data (the symbols) and green dashed line shows the expression (3.4).
3.2.2 The case $y_{p}\neq 0$
Knowing the stability of the particles which are positioned at the centreline of the channel, we now move to investigate the particles with 
$y_{p}>0$. Again as an example, the simulations for the case of 
$\unicode[STIX]{x1D709}=5,B=2$ and 
$y_{p}=0.5$ are shown in figure 7 for different yield numbers. Due to the asymmetry that the particles with 
$y_{p}\neq 0$ introduce to the flow, the breakage of the plug may be different from that of the cases investigated in the previous section (see figure 7b).
Velocity magnitude contours for 
$\unicode[STIX]{x1D709}=5,B=2$ and 
$y_{p}=0.5$: (a) 
$Y=0.75>Y_{c}$, (b) 
$Y=0.4<Y_{c}$ and (c) 
$Y=0.1<Y_{c}$. The interpretation of the white lines is the same as in figure 4.
The maps corresponding to the stability of the particles are sketched in figure 8(a–c). One can clearly conclude that the most important parameter here is 
$\bar{B}$. Indeed, when the particle is wholly inside the undisturbed plug (i.e. 
$B>\bar{B}$), then the critical yield number for the offset particles is the same as the computed 
$Y_{c}$ for the case 
$y_{p}=0$ (see figure 8a). However, the critical yield number for the particles that are partially in the undisturbed plug (i.e. the plug in the absence of any particle) increases (see figure 8b). This makes sense since even a small buoyancy stress can break the plug. However, for ultra-small buoyancy (large yield number), we observe that the plug grows locally to enclose the particle and stabilize it (see figure 8d). Increase in 
$Y_{c}$ for the particles which are positioned partially in the undisturbed plug is increasingly dramatic when 
$B\ll \bar{B}$ (i.e. when the particle is closer to the wall). For instance, the particle at 
$y_{p}=1.7$ considered in figure 8(c) is unstable even for 
$B=10$ and 
$Y=100$ (not shown in the figure).
(a–c) Different scenarios that can happen based on different yield numbers, Bingham numbers and positions of the particles (
$y_{p}$) for the case 
$\unicode[STIX]{x1D709}=5$: filled blue circles mark stable cases; squares represent cases similar to figure 7(b); and open blue circles represent the totally detached particles from the plug (e.g. the case in figure 7c). The red circles and lines are borrowed from figure 5 (showing 
$Y_{c}$ for 
$\unicode[STIX]{x1D709}=5$ and 
$y_{p}=0$) and the dashed red lines mark 
$\bar{B}$ associated with 
$y_{p}$ in each panel. (d) The velocity magnitude contours within the top half of the channel for the case 
$\unicode[STIX]{x1D709}=5,Y=1,B=2$ and 
$y_{p}=1$.
It is worth mentioning that although we just presented the results for the case 
$\unicode[STIX]{x1D709}=5$ for the sake of brevity, qualitatively quite similar behaviour was also observed for the other confinement ratios studied (
$2\leqslant \unicode[STIX]{x1D709}\leqslant 20$).
3.3 Neutrally buoyant particles inside Poiseuille flow
In the previous section, we have shown that 
$Y_{c}$ much increases for particles for which 
$\bar{B}\gg B$. In other words, we should expect that non-neutrally buoyant particles will lose their stability and break the plug if a portion of their surface is located outside the unperturbed yield surface, unless that portion is small. This being so, the question that naturally arises is how far from the position of the undisturbed plug can we expect any stable particle (i.e. how far the plug can grow locally to enclose the particles inside it). To answer this question we consider neutrally buoyant particles to address the most extreme limit: the stability map is shown in figure 9 for the cases 
$\unicode[STIX]{x1D709}=5$ and 
$\unicode[STIX]{x1D709}=8$.
Stability map for neutrally buoyant particles: (a,b) 
$\unicode[STIX]{x1D709}=5$ and (c,d) 
$\unicode[STIX]{x1D709}=8$. Circles show the position of the centre of the particles (
$y_{p}$): filled blue circles mark stable cases, while open blue circles mark the particles totally detached from the plug. The solid red lines show 
$y^{\ast }$ and the dashed red lines are 
$y^{\ast }$ shifted down by the radius of the particle (
$R=\frac{1}{2}$): below the dashed red lines particles are totally inside the undisturbed plug. The dashed blue line is a reference line in (a) which is 
$y^{\ast }$ shifted by 
$-0.22$. The green symbols and lines represent the values of 
$B_{c}$. In (b,d), 
$B_{c}$ as a function of 
$\bar{B}$ is plotted and the solid black lines are the fitted lines to the data: (b) 
$B_{c}\sim 0.274\bar{B}$ and (d) 
$B_{c}\sim 0.371\bar{B}$.
The boundary between the filled and open circles (
$y_{c}$ shown by green asterisks and lines) is the maximum distance from the centreline that we should expect neutrally buoyant particles to be stable inside the plug as a function of the Bingham number. Equivalently, the value of the Bingham number at which a particle with a considered position is just stable is designated by 
$B_{c}$. The interesting conclusions are as follows:
(i) Neutrally buoyant particles that feel
$B>\bar{B}$ are definitely stable.(ii) The plug cannot contain neutrally buoyant particles whose centre is out of the undisturbed plug. Not only those, but also some other particles that are closer to the centreline are unstable. For instance, considering figure 9(a), if
$\unicode[STIX]{x1D709}=5$ and $B=6.25$, then $y^{\ast }$ is approximately 1.71 but only particles that are below $y\approx 1.5$ are stable.(iii) As the Bingham number increases (
$y^{\ast }$ increases as well), we should expect particles with larger $y_{p}$ (far from the centreline) to be stabilized as well. On the other hand, the growth in $y_{c}$ (the critical position below which all the particles are stable) will be smaller and smaller as the Bingham number increases since the portion of the particles which are outside the undisturbed plug will experience larger forces as the Bingham number increases due to the highly sheared unyielded regions close to the walls. To explain by an example, considering the figure 9(a) again, increasing the Bingham number from 6.25 to 13.75 (by approximately 7.5) will increase $y_{c}$ from 1.5 to 1.7, while a further increase in $B$ by approximately 20.62 ($B\approx 34.37$) will increase $y_{c}$ by just 0.1 ($y_{c}\approx 1.8$). This is also clearly distinguishable by looking at the dashed blue line ($y^{\ast }-0.22$) in figure 9(a) which is closely approximate to $y_{c}$ up to $B\approx 15$ but will lose its validity beyond that.(iv) The critical Bingham number,
$B_{c}$, increases linearly with $\bar{B}$ as shown in figure 9(b,d).
The last conclusion can be explained using the same argument as expression (3.2) at the yield limit, except that here there is no buoyancy stress. Hence at the yield limit, roughly speaking, the shear stress contribution acting on one side of the particle surface is the same as the yield stress acting on the other side. Hence, in non-dimensional form
$$\begin{eqnarray}\displaystyle b\unicode[STIX]{x1D70F}_{xy}\sim bGy_{p}\sim B_{c}, & & \displaystyle\end{eqnarray}$$ where 
$b$ is the possible factor of proportionality. Considering that
$$\begin{eqnarray}\displaystyle \text{when }B=\bar{B}:y^{\ast }=y_{p}+R=\frac{\bar{B}}{G}, & & \displaystyle\end{eqnarray}$$it yields
$$\begin{eqnarray}\displaystyle B_{c}\sim b\bar{B}. & & \displaystyle\end{eqnarray}$$ The value of 
$b$ is 0.274 and 0.371 for the cases 
$\unicode[STIX]{x1D709}=5$ and 
$\unicode[STIX]{x1D709}=8$, respectively, as shown in figure 9(b,d).
4 Summary and discussion
The stability of neutrally and non-neutrally buoyant particles immersed in yield-stress Poiseuille flows is investigated numerically in the present study: particles can be transported in the core unyielded region without any relative motion compared to the carrying adjacent fluid if the yield stress is large enough compared to the buoyancy stress and the other stresses exerted on the particles. This exceptional feature of yield-stress fluids was previously addressed in the case of particle sedimentation in a quiescent bath of yield-stress fluids (see Putz & Frigaard Reference Putz and Frigaard2010; Chaparian & Frigaard Reference Chaparian and Frigaard2017a,Reference Chaparian and Frigaardb), which resulted in reporting 
$Y_{c}^{\ast }$ – the critical ratio of the yield stress to the buoyancy stress beyond which the particle does not move.
Numerical experiments in the present study revealed that particles in channel flows may lose their stability with substantially smaller buoyancy compared to the sedimentation case inside a bath of quiescent yield-stress fluid because of the inherent shear stress in the Poiseuille flow (i.e. 
$Y_{c}>Y_{c}^{\ast }$). The increase in 
$Y_{c}$ is more marked as the confinement ratio becomes smaller (see figure 5). However, we demonstrated that confinement by itself (in the absence of any pressure gradient/background flow) has an almost negligible effect for confinement ratios beyond 2, which is the main focus of the present study.
We demonstrated that for the particles positioned at the centreline, given that the confinement ratio is fixed, 
$Y_{c}$ decreases with increasing Bingham number. Indeed, the Bingham number has a stabilizing effect. The decrease in 
$Y_{c}$ on increasing 
$B$ is more apparent for large confinement ratios, while for the case 
$\unicode[STIX]{x1D709}=2$, almost no change in 
$Y_{c}$ was observed since the shear stress (which has a destabilizing effect) grows markedly with increasing Bingham number (i.e. 
$a(\unicode[STIX]{x1D709}=2)>1$). It has been also clarified that when 
$B\rightarrow \infty$, then 
$Y_{c}$ asymptotes to its limiting value, 
$Y_{c}^{\infty }$, which is a function of the confinement ratio. An analytical model was extracted from the numerical data for calculating 
$Y_{c}^{\infty }$.
For particles that are not at the centreline of the channel, it was shown that the position of the particle relative to the position of the undisturbed yield surface (i.e. comparison between 
$y_{p}$ and 
$y^{\ast }$ or between 
$\bar{B}$ and 
$B$) is a key parameter for the stability of the particle. As far as the whole particle is located inside the undisturbed core plug region (
$y_{p}\leqslant y^{\ast }-1/2$), the above statements are valid. However, if a particle portion is situated in the undisturbed yielded regions (i.e. 
$y_{p}>y^{\ast }-1/2$), then 
$Y_{c}$ can dramatically increase. The more volume a particle shares with the undisturbed sheared regions, the larger 
$Y_{c}$ will be. Yet, if the yield stress is high enough (i.e. large 
$Y$), the core unyielded plug can grow locally to some extent to contain the adjacent particles. This picture holds for neutrally buoyant particles as well. Indeed, in the absence of any density difference between the yield-stress fluid and the particle, particles that are totally in the undisturbed core plug region are strictly stable. The core unyielded plug can also form bulges to contain the neighbouring particles and prevent them from being released into the sheared regions.
We comment on the current findings, in terms of applicability, as follows:
(i) As figures 5 and 6(c) clearly illustrate, when
$\unicode[STIX]{x1D709}\rightarrow \infty$ then $Y_{c}^{\infty }\rightarrow Y_{c}^{\ast }$. Simultaneously, $Y_{c}$ reaches the limiting value (i.e. $Y_{c}^{\infty }$) at smaller Bingham numbers. Hence, we can conclude that for suspension transport problems for which $\unicode[STIX]{x1D709}$ is large (say $\unicode[STIX]{x1D709}\gtrsim 25$), at adequate Bingham numbers (e.g. $B>O(1)$), one can use $Y_{c}^{\ast }$ as $Y_{c}$. However, if one is dealing with a problem of $\unicode[STIX]{x1D709}\approx O(1)$, then calculating correct $Y_{c}$ (e.g. using figure 5) is of great importance. This is shown schematically in figure 10.(ii) Unstable particles may migrate in the yielded sheared regions as previously evidenced experimentally by Merkak et al. (Reference Merkak, Jossic and Magnin2009). We did not cover this behaviour in the present study, as it requires a detailed investigation of its own. We briefly mention that those non-neutrally particles that are positioned in the bottom half of the channel and are unstable move towards the bottom wall as always
$Y_{c}$ is larger close to the yield surfaces compared to $Y_{c}$ at the centreline. However, those particles starting at the top half of the channel and are unstable may be trapped again in the core unyielded region (e.g. centreline) if $Y$ is large enough (e.g. compared to $Y_{c}$ at the centreline).(iii) We remind the reader that the present study is a two-dimensional one. This restriction was chosen for two main reasons. First of all, less information about the static stability of three-dimensional particles can be found in the literature within the context of sedimentation problems in a bath of a quiescent yield-stress fluid which is the baseline reference point. Moreover, direct numerical simulation of any three-dimensional multiphase flow problem would be extremely expensive using the augmented Lagrangian approach which is vital for capturing the yield limit accurately. As has been experienced by Siqueira & de Souza Mendes (Reference Siqueira and de Souza Mendes2019), the regularized approaches cannot fully address the stability of particles.
The present study provides a baseline/fundamental result from which to look at more complex cases representative of applications in future studies.
Schematic of the suspension transport problem with yield-stress carrying fluid. Panel (a) represents 
$\unicode[STIX]{x1D709}\rightarrow \infty$ and 
$B\geqslant O(1)$ while (b) corresponds to 
$\unicode[STIX]{x1D709}\approx O(1)$.
Acknowledgements
E.C. gratefully acknowledges the support of a Linné FLOW PostDoc grant during the course of this study. O.T. acknowledges the support of the Swedish Research Council through grant VR 2017-0489.
Declaration of interests
The authors report no conflict of interest.
Open access
