Sliding flows of yield-stress fluids

Abstract A theoretical and numerical study of complex sliding flows of yield-stress fluids is presented. Yield-stress fluids are known to slide over solid surfaces if the tangential stress exceeds the sliding yield stress . The sliding may occur due to various microscopic phenomena such as the formation of an infinitesimal lubrication layer of the solvent and/or elastic deformation of the suspended soft particles in the vicinity of the solid surfaces. This leads to a ‘stick–slip’ law which complicates the modelling and analysis of the hydrodynamic characteristics of the yield-stress fluid flow. In the present study, we formulate the problem of sliding flow beyond one-dimensional rheometric flows. Then, a numerical scheme based on the augmented Lagrangian method is presented to attack these kind of problems. Theoretical tools are developed for analysing the flow/no-flow limit. The whole framework is benchmarked in planar Poiseuille flow and validated against analytical solutions. Then two more complex physical problems are investigated: slippery particle sedimentation and pressure-driven sliding flow in porous media. The yield limit is addressed in detail for both flow cases. In the particle sedimentation problem, method of characteristics – slipline method – in the presence of slip is revisited from the perfectly plastic mechanics and used as a helpful tool in addressing the yield limit. Finally, flows through model and randomized porous media are studied. The randomized configuration is chosen to capture more sophisticated aspects of the yield-stress fluid flows in porous media at the yield limit – channelization.


Introduction
Slip behaviour in complex fluids has been attributed to various microscopic phenomena such as chain detachment/desorption at the polymer-wall interface in polymer melts (Hatzikiriakos 2012), migration/depletion of disperse phase away from the vicinity of the solid boundaries in suspensions (Vand 1948), and elastic deformation of the suspended soft particles lying on the smooth solid boundaries (Meeker, Bonnecaze & Cloitre 2004b) in pastes -yield-stress fluids. In the present study, however, we rather investigate the macroscopic consequences of slip on hydrodynamic features of yield-stress fluid flows, in order to form a bridge over these two distinct scales.
In experiments, roughened (e.g. using sandpapers) or chemically treated surfaces are usually used to suppress slip in measurements of yield-stress fluid flows. For instance, Christel et al. (2012) proposed a polymethyl methacrylate (known as PMMA) treatment to reduce slip by exciting positive surface charges providing electrostatic interactions and squashing the lubricating layer. This is important since in a great number of laboratory experiments, polymethyl methacrylate is used due to its transparency properties which enable flow visualizations. In industrial scales/processes, however, chemical treatments are not always feasible, and therefore it is crucial to analyse the hydrodynamic consequences of sliding flows. Moreover, slip alters rheometric measurements, causing rheological properties of the samples to be inaccurately evaluated, sometimes even by one order of magnitude (Poumaere et al. 2014). This issue has been discussed extensively in the literature, for instance recently by Medina-Bañuelos et al. (2017, 2019. Hence, it would be constructive to have a generic analysis framework of sliding flows in yield-stress fluids, and this is the aim of the present paper. In the present study, we formulate the well known 'stick-slip' law for yield-stress fluids beyond one-dimensional (1-D) rheometric flows and then generalize the previously proposed numerical algorithm (Roquet & Saramito 2008;Muravleva 2018) based on the augmented Lagrangian method for attacking this kind of problem. Theoretical tools are developed in the presence of slip. The whole framework is first benchmarked by a simple channel Poiseuille flow. Then we proceed to investigate flows in complex geometries: the creeping flows about a circular cylinder and the pressure-driven flows in model and randomized porous media.

General slip law for yield-stress fluids
Characterizing the slip behaviour of yield-stress fluids has been the main objective of a large number of studies from theoretical attempts to experimental measurements. Meeker et al. (2004b) and Piau (2007) tried to theoretically describe the slip origin in soft particle pastes/microgels from an elastohydrodynamic perspective. Rheometric tests with slippery plates in the cone-plate/plate-plate geometries have been conducted to quantify the slip behaviours of yield-stress fluids ranging from Carbopol gels to emulsions (Meeker, Bonnecaze & Cloitre 2004a;Meeker et al. 2004b;Poumaere et al. 2014;Zhang et al. 2018). Moreover, very recently, in a series of experiments using optical coherence tomography (known as OCT) and particle tracking velocimetry (known as PTV), Daneshi et al. (2019) characterized the slip behaviour of Carbopol gel (a model 'simple' yield-stress fluid) inside a capillary tube. The sliding characteristics of yield-stress fluids observed in all the mentioned 1-D experimental studies can be summarized in a general slip law aŝ usually termed as the 'stick-slip' law, whereû s is the slip velocity on the solid surface andτ w is the shear stress at the solid boundaries. The sliding threshold is called the sliding yield stress -τ s . The slip coefficient and the power index are designated byβ s Figure 1. Schematic of the sliding channel flow of a yield-stress fluid; the x-axis is aligned and placed at the axial centreline of the channel which is formed by the two infinite parallel plates separated by the gapĤ and the y-axis is in the wall-normal direction: (a) whenτ w τ s there is no flow, (b) whenτ s <τ w τ y then the fluid in the whole gap slides as an unyielded plug, (c) whenτ y <τ w then the fluid in the vicinity of the walls (ŷ p <ŷ whereτ y <τ xy ) yields and at the same time slides over the walls. There is a core unyielded region as well (ŷ ŷ p ). and k, respectively. Hence,β s has the dimension of m 2k+1 /N k .s or m/Pa k .s. Its physical interpretation (for the case k = 1) is the slip length over the local effective viscosity of the fluid,ˆ s /μ eff . In the entire paper, quantities with a 'hat' symbol (·) are dimensional and others are dimensionless.
For a deeper understanding of the physical meaning of the slip law (1. whereK(Pa · s n ) is the consistency of the fluid and n the power-index (Herschel-Bulkley model),û the streamwise velocity,τ y the material's yield stress,τ xy the shear stress and sgn(·) and abs(·) are the sign and absolute value functions. Figure 1 represents the slip law schematically. When the applied pressure gradient is small enough, the wall shear stress (τ w ) lies below the sliding yield stress (figure 1a) and the slip velocity is zero. Sinceτ s is always less than the material's yield stress, then there is no flow in the channel. When the pressure gradient is increased beyond (Δp/L) c where subscript ' c ' indicates the critical value, the wall shear stress grows beyondτ s and the fluid slides over the walls. If the wall shear stress is still less thanτ y , then the material moves as a sliding unyielded plug with a constant velocity in the entire gap (see figure 1b). When the pressure gradient is increased, at some point, the wall shear stress exceeds the fluid yield stress and the sheared/yielded regions appear, which at the same time slide over the walls -see figure 1(c). Yet, in the vicinity of the centreline, the shear stress drops below the yield stress and a core unyielded region is formed. The flows corresponding to figures 1(b) and 1(c) are usually termed the fully plugged regime and the deformation regime, respectively. In the present study, the objective is to systematically address the effect of slip on yield-stress fluid flows and its consequences on the yield limit; not only by conducting numerical simulations, but also by forming a theoretical framework for analysing this type of problem. Although some attempts have been made to solve simple rheometric flows in previous years (Philippou, Kountouriotis & Georgiou 2016;Panaseti & Georgiou 2017;Damianou, Panaseti & Georgiou 2019), still the lack of generic numerical implementation frameworks (especially non-regularized methods) and the absence of adequate theoretical frameworks to attack sliding yield-stress fluid flows are discernible. Hence, in what follows, we aim to partially fill this gap in the literature. Beyond that, we study complex sliding flows such as particle sedimentation and pressure-driven flows in porous media within these frameworks.
The outline of the paper is as follows. In § 2, we set out the slip law beyond 1-D for the yield-stress fluid flows and generalize the numerical algorithm previously proposed for simple 1-D problems. We then benchmark the numerical algorithm for simple channel Poiseuille flow in § 3 and derive some theoretical tools for investigating the sliding flows. This will be followed by addressing the sliding flow about a circular particle and the aspects of the particle yield limit in § 4. The next two sections are devoted to the sliding flows in porous media: in § 5 some general hints are made by studying the flow through model porous media and in § 6 some deeper investigations are carried out by performing numerical simulations in randomized porous media. Finally, some conclusions are drawn in the closing section § 7.

Slip law and the numerical algorithm beyond 1-D
In this section, we set out and formulate a sample problem of sliding flow of a yield-stress fluid and generalize the slip law and numerical algorithm for solving such problems beyond 1-D. In what follows, we mainly assume that n = 1 (Bingham model). Later in appendix C, we will address how the present framework can be expanded to encompass the other 'simple' yield-stress rheological models, e.g. Herschel-Bulkley. However, the majority of the physical discussions (e.g. the yield limit) in the next sections are independent of the value of n. We will come back to this point in § 7.
Problem P. We consider the Stokes flow of a Bingham fluid in Ω \X (see figure 2) as follows: wherep is the pressure,τ the deviatoric stress tensor,ρf the body force,ρ the fluid density andγ is the rate of deformation tensor. Without loss of generality, we consider no-slip and no-penetration boundary conditions on ∂Ω asû =û 0 and the slip boundary condition, on ∂X, whereû s is the slip velocity,û ns the velocity of the solid boundary ∂X and δû is the restriction ofû on ∂X:û → δû asx → ∂X. The value of the tangential traction vector (i.e. tangential force per unit area on the solid surface) is shown byΛ where the normal and tangential unit vectors to the solid surface ∂X are represented by n and t, respectively. The no-penetration condition on ∂X reads δû · n =û ns · n. Without loss of generality, we fixed k at unity in the slip law (2.3) and will do so in the rest of the present study. This is supported by experimental observations as well. Seth et al. (2012)  et al. (2019) reported k ≈ 1 for the pressure-driven flows of Carbopol gels with different concentrations and stirring rates during the sample preparation. However, experiments have also revealed that k / = 1 in some cases (Piau 2007;Medina-Bañuelos et al. 2017, 2019. We will discuss these cases in appendix C, as the main objectives/conclusions of the present study are independent of the value of k.

Numerical algorithm
The minimization principle for the sliding problem P can be derived as For the augmented Lagrangian approach we require to relax the rate of strain tensor to the auxiliary tensorq and the restriction of the velocity vector on ∂X to the auxiliary vectorξ . Then the saddle-point problem associated with P is where a and b are arbitrary constants (augmentation parameters);T andλ are the Lagrange multipliers. Hence, the Uzawa algorithm in the presence of slip takes the form of algorithm 1. Upon convergence of algorithm 1 with the free augmentation parameters a and b, the Lagrange multiplierT converges to the true stress field,q to the true rate of strain tensor, Lagrange multiplierλ to the traction vector on ∂X, and auxiliary variableξ to the velocity on ∂X.
A mesh adaptation procedure, the same as the one proposed by Roquet & Saramito (2003), could be coupled with the above algorithm to obtain a fine resolution of the yield surfaces. Based on this procedure, the adapted mesh is stretched anisotropically in the Algorithm 1 1: procedure (solving yield-stress fluid flow with slip yield stress boundary condition) 2: m ← 0 3:q 0 ,T 0 ,λ 0 ,ξ 0 ← 0 (or any other initial guess) 4: loop (Uzawa algorithm):
iff |ˆ | τ s , û ns · n n + û ns · t t +β direction of the eigenvectors of the Hessian of μγ :γ +τ y γ , which is the square root of the local energy dissipation. In this study, we implement the entire numerical algorithm in an open-source C++ finite element environment -FreeFEM++ (Hecht 2012). We have previously validated our numerical implementation and mesh adaptation widely within various studies (Chaparian & Frigaard 2017b;Chaparian & Tammisola 2019;Chaparian et al. 2020). Here we will not go into technical details such as the convergence rate, since such details are well-documented in the previous studies -please see, for example, Roquet & Saramito (2008).
In what follows, we quickly validate the algorithm for the sliding flows and most importantly derive variational tools (Mosolov & Miasnikov 1965) for a sample problem. These tools will be used to analyse more complex flows in the next sections.

Benchmark problem: sliding channel Poiseuille flow
In this section we consider the sliding flow of a Bingham fluid in a channel (Poiseuille flow) -the same as the one shown in figure 1. The walls are denoted by Γ and the full domain by Ω. We consider two classic formulations: (i) the mobility problem ([M]) where the pressure gradient is applied and the flow rate can be computed and (ii) the resistance problem ([R]) in which flow rate is set and as a result pressure gradient can be computed. We also validate the presented numerical algorithm and derive/revisit some variational tools (Huilgol 1998) useful for the rest of this study.

The [M] problem
In this subsection, we useĜĤ 2 /μ as the velocity scale,ĜĤ as the characteristic viscous stress to scale the deviatoric stress tensor andĤ/μ to scale the slip coefficient wherê G is the absolute value of the applied pressure-gradient to drive the flow from left to right (in the positive direction of the x-axis). Hence, the non-dimensional governing and constitutive equations take the forms e x + ∇ · τ = 0 (3.1) and respectively, and the slip law where Od =τ y /ĤĜ and Od s =τ s /ĤĜ. In what follows, the ratioτ s /τ y is denoted by α s . The analytical solutions to the velocity and stress fields are derived in appendix A.1.
The energy balance equation in the presence of slip is where a(u, u) = Ωγ (u) :γ (u) dA and Od j(u) = Od Ω γ (u) dA are the viscous and plastic dissipations, respectively. Please note that the slip dissipation ∂X |σ nt u s | dS can be split into two terms: the 'viscous' slip dissipation (1/β s ) Γ u 2 s dS which is designated by a s (u s , u s ) in this study and the 'plastic' slip dissipation Od s Γ |u s | dS with Od s j s (u s ). Moreover, j(u) is the normalized plastic dissipation and j s (u s ) the normalized 'plastic' slip dissipation.
We can rearrange the energy balance equation and form a set of inequalities as follows: (3.5) to find the definition of the critical Oldroyd number in the presence of slip as where v is a velocity test function from the set of all admissible velocity fields V . Then it enforces that for Od c Od, a(u, u) = 0 or u = 0 (i.e. no flow condition). We may use (3.6) for calculating Od c in simple flows such as Poiseuille flow here: the flow can be postulated by two boundary layers with thickness δ attached to the walls, through which the slip velocity U 1 is connected to the plug velocity U 2 , then, where is the length of the chosen control volume. This suggests, or, . (3.9) The maximal value of the argument is 1/2α s which occurs at U 1 /U 2 = 1. Its physical interpretation is that the supremum occurs in the fully sliding regime, which is intuitive. Hence, Od c = 1/2α s , or in other words, there is no flow if 1/2 Od s , which is in agreement with the solution of the full equations derived in appendix A.1.
We perform a few simulations for the [M] problem to benchmark the presented algorithm. Figure 3 illustrates different regimes: figure 3(a) compares the computed velocity profile (blue line) in the gap with the analytical solution (red circles) in the fully sliding plug regime (Od = 1, α s = 0.2 & β s = 0.1) and figure 3(b) within the deforming regime. Utilizing the mesh adaptation seems more indispensable in the deforming regime compared with the fully sliding plug regime, where velocity gradient is zero. It is clear that the velocity distribution after adaptation (red line) is much closer to the analytical solution (red circles) than the solution of the initial mesh (blue line), and also yields a fine resolution of the yield surfaces -see figure 3(d). We shall mention that the focus of the present study is not on quantifying these improvements by the mesh adaptation, as it has been previously studied in detail (Roquet & Saramito 2003Chaparian & Tammisola 2019).

The [R] problem
Here, we reformulate the problem in the resistance scalings and the quantities are designated with asterisk sign ( * ) to prevent any potential confusion with the [M] problem. In the resistance formulation, the velocity is scaled with the mean bulk velocity,Û =Q/Ĥ, and as a result, there is always a non-zero flow with the flow rate equal to unity (Q * = 1), whether the flow is in the deforming regime or the fully sliding plug regime. The deviatoric stress and pressure are scaled with the characteristic viscous stress,μÛ/Ĥ, hence, with the boundary condition, on the walls govern the [R] problem. Please note that G * is the absolute value of the pressure gradient which enforces the unity flow rate for a given Bingham number (τ yĤ /μÛ), α s and β s . The energy balance equation in this case is (3.13) Frigaard (2019) has shown that at the yield limit (B → ∞), viscous dissipation is at least one order of magnitude less than the plastic dissipation (when B → ∞: a(u * , u * ) Bj(u * )). With a little extra effort, we can show that this applies to a s (u * s , u * s ) and B s j s (u * s ) as well, knowing that u * s 1 (from the continuity equation Q * = 1) and B s → ∞ in this limit. Exploiting that generally the mapping between [M] and [R] problems is attainable via Q = Od/B or G * = B/Od, one can extract the critical Oldroyd number from the [R] formulation through . (3.14) To numerical data to extract the yield limit. This is illustrated in figure 4 by the dotted lines, where at B → ∞, G * asymptotes to 2B and 2α s B = 0.4B for no-slip and the sample sliding flow considered, respectively. In this limit, j(u * ) (dashed lines in figure 4) asymptotes to 2 and 0 for the no-slip and sliding flows: this is equivalent to Od c = 1/2 and 1/2α s , respectively, knowing that j s (u s ) asymptotes to 2 for the sliding flows at sufficiently large Bingham numbers (follow the red dashed-dotted line in figure 4). These results are the same as the ones derived in appendix A.1 considering the [M] formulation. Please note that due to the log-log scaling used in figure 4, j(u * ) corresponds to the sliding flow, which is exactly equal to zero forB B is not visible; however, its fast decay is clear: From a different perspective, the inset in figure 4 shows how mesh adaptation helps us to achieve more accurate results.

Slippery particle motion
Particle motion in yield-stress fluids has been considered in numerous studies ranging from numerical simulations (Beris et al. 1985 Unfortunately, fewer studies are devoted to investigating the effect of slip on the flow field around the particle surface. Although Fraggedakis et al. (2016) used the Navier-slip law in their simulations, the emphasis of their study was on the effect of elasticity of the yield-stress fluid in the particle sedimentation problem. On the experimental side, Jossic & Magnin (2001) reported the drag coefficient of the particle moving in a Carbopol gel and showed that for particles with smooth surfaces, the drag coefficient is smaller compared with the rough particles that prevent slip which is intuitive. However, the mentioned authors did not quantify the slip velocity/law on the particle surface and rather qualitatively addressed the effect of slip. In this section, we will shed light on the effect of slippery motion of the yield-stress fluid about a circular two-dimensional (2-D) particle and in detail will address the yield limit under this condition. Flow configuration: the horizontal symmetry axis of the particle is aligned with the x-axis (positive direction to the right) and the vertical one by y-axis (positive direction upward) where the origin of the coordinate system is fixed at the particle centre. The fluid flow about the particle (X) with its radiusR as the length scale can be described again by two formulations.
(i) The [R] formulation, where the particle is dragged through the fluid with a prescribed velocityV * p as the boundary condition. Hence, the drag forceF * D can be calculated. In this formulation,V * p is used as the velocity scale andμV * p /R as the characteristic viscous stress to scale the deviatoric stress tensor and the pressure, hence the governing equation is with (3.11) and (3.12) as the constitutive and the slip law, respectively, where B = τ yR /μV * p . The drag force on the particle can be computed as a result from where the buoyancy of the particle is known and it falls freely under the action of the gravity −ĝe y . Hence, the settling velocity of the particleV p can be calculated as a result. In this formulation, velocity is scaled with ΔρĝR 2 /μ: the velocity scale obtained via balancing the buoyancy stress with the characteristic viscous stress. Please note that hereĝ is the gravitational acceleration and Δρ is the density difference between the particle and the suspending fluid. Hence, the governing equation is where ρ is the ratio of the fluid density to the particle one (i.e.ρ f /ρ p ). The constitutive and slip laws are the same as (3.2) and (3.3), respectively, with a slight difference: the yield number Y is substituted for the Oldroyd number in the particle sedimentation problem; indeed, Y =τ y /ΔρĝR.
The yield limit can be well understood using either the [R] or [M] formulations knowing the mapping between these two formulations: In the [R] formulation, the particle never stops because the particle velocity is set as the boundary condition, yet, in the yield limit B → ∞, the drag force also goes to infinity. However, in the [M] problem, for a large enough yield number (Y > Y c ), the yield stress can overcome the buoyancy of the particle and makes it stationary suspended. The connection between these two situations can be well discussed by the 'plastic drag coefficient', for problem [R], where ⊥ is the width of the particle perpendicular to the flow direction ( ⊥ = 2) andÂ p is the area of the 2-D particle (i.e. πR 2 here). Then, the critical plastic drag coefficient (i.e. plastic drag coefficient at the yield limit) can be converted easily to the critical yield number, (4.5) For more details please see Chaparian & Frigaard (2017b). For instance, for a 2-D circle, as has been extensively validated (Tokpavi et al. 2008;Putz & Frigaard 2010;Chaparian & Frigaard 2017a,b), the critical plastic drag coefficient under the no-slip condition is 11.94, which yields to Y c = 0.1316. Regarding the computational details, we conduct the same strategy as in our previous studies (Chaparian & Frigaard 2017a,b;Iglesias et al. 2020). The computational box (Ω) is chosen large enough to ensure independency from far-field conditions (Chaparian, Wachs & Frigaard 2018). The boundary condition u * = 0 is enforced on ∂Ω.
The yield limit of particle motion in a quiescent yield-stress fluid is directly relevant to the lateral resistance of an object (e.g. a pile) in a perfectly plastic medium (e.g. soil) which can be investigated by method of characteristics (i.e. slipline solution) due to the hyperbolic nature of the equations. We briefly revisit this problem in the next subsection.

Revisiting the slipline solution and lower-and upper-bound calculations
In a deformed perfectly plastic medium, the second invariant of the stress is equal to B everywhere. The equilibrium equations for the unrestricted 2-D plastic flow then form a closed set of hyperbolic equations for which there are two families of orthogonal characteristic lines: the αand β-families (Hill 1950;Chakrabarty 2012). Physically, these lines (termed as sliplines) show the direction of the maximum shear stress which is equal to B. Properties of these lines facilitate studying the 'yield limit' in various viscoplastic problems (Dubash et al. 2009;Chaparian & Frigaard 2017a,b;Chaparian & Nasouri 2018;Hewitt & Balmforth 2018). However, finding the sliplines itself is not trivial in some complex problems. Indeed, the sliplines assist us in postulating admissible stress and velocity fields which yield to finding the lower and upper bounds, respectively, of the load limit.
Initially investigated by Randolph & Houlsby (1984), a slipline solution was devised for calculating the lateral resistance of a circular pile in soil (considered as a perfectly plastic material), taking into account the effect of soil adhesion at the pile-soil interface, α s . In other words, the tangential shear stress at the pile-soil interface is assumed to be equal toα s B. Although widely believed that the Randolph & Houlsby (1984) solution is exact (lower and upper bounds of the load/drag are the same for the whole range ofα s ), an issue Only a quarter of the whole domain is shown due to symmetry.
was firstly detected by Murff, Wagner & Randolph (1989): forα s < 1, there is a region in the vicinity of the pile surfaces in which the rate of strain is negative and its absolute value should be taken into account in calculating the upper bound to avoid negative plastic dissipation. If one does so, then there is a discrepancy between the lower and upper bounds. Later, Martin & Randolph (2006) proposed two new postulated velocity fields which improved the upper-bound predictions to some extent; one for small (α s → 0) and the other one for large (α s → 1) soil adhesion factors. By combining these two mechanisms, Martin & Randolph (2006) proposed an alternative mechanism which markedly shrinks the uncertainty between the lower and upper bounds predictions for the entire range ofα s . We will quickly review the lower- (Randolph & Houlsby 1984) and upper-bound (Martin & Randolph 2006) solutions here; however, for more details, readers are referred to the original references. The first family of sliplines in this problem are straight lines that make a π/4 − ψ/2 angle with the pile surface (say α-lines) where ψ = sin −1α s . Figure 5 represents these lines in cyan colour. The initial α-line is AB which makes a π/4 angle with the vertical symmetry line sinceτ xy is zero on OB. Hence, the region enclosed between AB and the pile surface is the plug region shown in light grey. The α-lines from AB to CD can be found easily as discussed above; the corresponding β-lines are curved lines which are the involutes unwrapped from an imaginary concentric circle -the evolute -with radius η = cos( 1 2 cos −1α s ) shown with a dashed white line. Please note that indeed α-lines can be introduced as tangents of this imaginary circle as well. From CD to CE, α-lines are spokes of a fan centred at C; hence, the β-lines in this part are arcs of concentric circles at point C.
Hencky's equations (Hill 1950), p + 2Bθ = const. along an α-line andp − 2Bθ = const. along a β-line, (4.6) provide the tool for finding the admissible stress field associated with these sliplines which can be used for calculating the lower bound of the plastic drag coefficient. Here, θ is the anticlockwise orientation of the α-lines made with the x-axis (aligned with OE).
The individual stress contributions can be written as in the new curvilinear coordinates α and β. Then the lower bound can be calculated as which basically consists of four individual contributions in the shown quadrant: shear force acting on AB, B sin(ψ/2); shear force acting on AC, B sin ψ cos(ψ/2); normal force acting on AB, B(c + 3π/2) sin(ψ/2); and normal force acting on AC, where cB is the reference mean stress, (σ x +σ y )/2, on CE (position of the centre of the Mohr circle). Please note that c will be eliminated if all four quadrants are taken into account. Hence, the total lower bound of the plastic drag coefficient will be (4.10) The upper-bound calculation can be performed by postulating admissible velocity fields. Geiringer equations (Hill 1950) make it feasible via slipline solution, dũ α −ũ β dθ = 0 along an α-line & dũ β +ũ α dθ = 0 along a β-line, (4.11) whereũ α andũ β are the velocities in the directions α and β in the new curvilinear coordinate system. Based on these equations, Randolph & Houlsby (1984) proposed a mechanism (see figure 6a) in which velocity along α-lines is zero and along each β-line is a constant, which can be found by considering a no-penetration condition on the pile or fore-and-aft plugs surface. As mentioned above, this mechanism leads to a relatively large uncertainty with the lower-bound solution forα s / = 1. Martin & Randolph (2006) improved this mechanism by substituting a rotating rigid block (blue region in figure 6b) in the middle of the domain with angular velocity ω compatible with the no-penetration boundary condition on the pile surface; the centre of that is determined by the angle γ and the imaginary evolute circle of radius η. By optimizing for γ , one can markedly reduce the uncertainty between the lower and upper bounds of C P D,c . It is worth mentioning that althoughα s and α s do not have exactly the same physical origins, in what follows, we show that both display the same effect in the current problem and can be used interchangeably. Nevertheless, we should note that, α s B in the viscoplastic fluid mechanics context not only marks the boundary between the presence of sliding or no-slip flow on the solid surface, but actually also controls the yield limit (see figure 1): for instance in the case of α s = 1, the limit of flow/no-flow reduces to having stress as large as B on the solid surfaces. Without a doubt, when there is a flow, it slides over solid boundaries depending on β s in this case. Yet, the case of α s = 1 shares the same yield limit characteristics with the no-slip condition. In the perfectly plastic mechanics context, α s B roughly has the same interpretation since if the tangential shear stress on the solid boundary is less than that, then there are no deformed regions. Indeed, this is the strict condition at the yield limit in a deformed perfectly plastic solid.  (Martin & Randolph 2006) with the blue region rotating as an unyielded block. Black arrows show velocity vectors. For a better representation, in this schematic the pile velocity is assumed to be e y . Nevertheless, it should be noted that the flow has fore-and-aft symmetry. Figure 7(a) represents the plastic drag coefficient versus the Bingham number under different conditions. As shown by Putz & Frigaard (2010) and Chaparian & Frigaard (2017b), with the no-slip condition on a circular particle surface, lim B→∞ C p D = C p D,c ≈ 11.94 which is the value that can be obtained from the expression (4.10) as well wheñ α s = 1. Three main points can be drawn from figure 7(a) as follows.

Results
(i) Plastic drag coefficients associated with the sliding flows over the particle surface are smaller compared with the same Bingham number flows with the no-slip condition. This seems intuitive and as mentioned has been validated experimentally by Jossic & Magnin (2001). (ii) The plastic drag coefficient decreases by increasing Bingham number and reaches a plateau when B → ∞. More interestingly, the critical plastic drag coefficient associated with the flows of the same α s asymptote to a same limiting value (see asterisks and circles). Indeed, C p D,c is only a function of α s . However, in the Newtonian limit, B → 0, this is β s which controls C p D (see stars and circles). (iii) The lower-bound estimations from the slipline solutions accurately predicts C p D,c . Figure 7(b) compares the lower and upper bounds of the critical plastic drag coefficient with the data obtained by the numerical computations. It clearly shows that even within the sliding flow context, perfectly plastic theories could be useful in investigating the yield limit. It is worth mentioning that, as investigated by Chaparian & Frigaard (2017b) in detail, although the computations and slipline analysis display same behaviours to some extent for the no-slip case, the envelope of the characteristics does not agree fully with the yield surfaces in the viscoplastic computations. Both the characteristics solution and the viscoplastic computation predict rigid caps fore-and-aft of the particle (although they are not exactly the same size), whereas the viscoplastic solution also has kidney plug regions along the side of the cylinder. The velocity fields are not exactly equal eithersee figure 4 of the mentioned reference. The perfectly plastic velocity has discontinuities which are diffused in the viscoplastic field by the narrow viscous boundary layers. Given the discrepancies in both stress and velocity fields for the no-slip case, it is perhaps interesting that slipline solutions still predict the yield limit very well. Figure 8 reveals that these differences exist in the sliding flows as well. Figure 8 (Martin & Randolph (2006) combined mechanism with continuous red line and Randolph & Houlsby (1984) mechanism with discontinuous dashed red line) over the horizontal symmetry line. Please note that particle surface velocity is v * = −1.
compares velocity distributions along the horizontal symmetry line: the two red lines represent the ones associated with the upper-bound mechanisms by Martin & Randolph (2006) (continuous) and Randolph & Houlsby (1984) (discontinuous), while the computed velocities for different Bingham numbers are shown in blue colours; the darker it is the higher the Bingham number. As B increases, the velocity distributions get closer to the ones predicted by the upper-bound mechanisms, especially the combined mechanism where even the slip velocity on the particle surface at the equator is predicted correctly. Nevertheless, again the discontinuities in the perfectly plastic solution are replaced by the viscous boundary layers in the viscoplastic solutions.

Sliding flow in model porous media
Non-Newtonian fluid flows through porous media, especially yield-stress fluids, are of great practical importance for numerous industries such as filtration. Due to the complexities of these types of geometries and their randomness, it is expensive to conduct numerical simulations/experiments at the full scale. Hence, a common way is to model the medium as arrays of cylinders or other model geometries. For instance, Chevalier et al. (2013)  the channel length), therefore, is an inadequate approach to study yield-stress fluid flows along wavy channels.
It should be mentioned that all the above studies are limited to the no-slip condition. Here, in the continuation of the previous sections, we investigate the effect of slip on the pressure-driven yield-stress fluid flows in porous media.
In this section, we consider flows through two model geometries mimicking porous media -the same ones as studied by Chaparian et al. (2020). The schematic is represented in figure 10: both geometries are designed to have the same porosity φ = 0.38 (i.e.ˆ /R = 0.25 and 1.25 in figures 10a and 10b, respectively). The obstacle's radius,R, is used as the length scale. We use the [R] formulation in this section; hence, for all the cases of the same geometry, the flow rate is constant and the pressure drop is a function of the Bingham number (τ yR /μÛ), α s and β s . We consider an inertialess fully developed flow with the sliding boundary condition (3.12) on the solid topologies (∂X). Fixed flow rate condition is enforced using a Lagrange multiplier and periodic boundary conditions are applied at the inlet and outlet. On the top and bottom sides of the computational cell (horizontal faces), a symmetry boundary condition is imposed. For more details about the computational method and its validations please see Chaparian et al. (2020). Figure 11 compares different flows within a wide range of Bingham numbers. Dead zones (shaded in light grey) are smaller in the sliding flows since the fluid slides over the solid surfaces and 'erodes' the fouling parts. However, the unyielded regions in the middle of the passage (yield surfaces in green) are mostly larger compared with the no-slip condition, which is intuitive. Interestingly, because of the slip, unyielded zones are prone to attach to the solid surfaces (e.g. see figure 11p). Figure 12 reveals, quantitatively, various features of the flow in the model porous media. Figure 12(a) sketches G * with respect to the Bingham number; zoomed-in insets are provided in figure 12(b-e). In the range of small Bingham numbers -Newtonian limitthe pressure gradient is a function of B and β s alone; see figures 12(b) and 12(d) in which the cyan and yellow curves (β s = 0.05) converge together and also red and purple curves (β s = 0.1) together. On the other hand, in the limit of large Bingham numbers (i.e. no flow or yield limit; figures 12c and 12e), the critical pressure gradient is a function of B and B s or indeed B and α s : yellow and red curves with α s = 0.2 are almost indistinguishable; the same for the purple and cyan curves with α s = 0.8. Moreover, as shown by Chaparian et al. (2020), in the yield limit (B → ∞), the pressure gradient scales with ∼B; this is true for both the no-slip and sliding flows, yet with different prefactors. Indeed, as shown in 911 A17-19 E. Chaparian and O. Tammisola Figure 11. The velocity magnitude contour |u * |. Panels (a-d) and (i-l) show the no-slip flow in the regular and staggered geometries, respectively. Panels (e-h) and (m-p) illustrate the sliding flow (α s = 0.2, β s = 0.1) in regular and staggered geometries, respectively. Panels (a,e,i,m), (b, f ,j,n), (c,g,k,o) and (d,h,l,p) are associated with B = 1, 10, 100 and 1000, respectively. The dead zones are shown in light grey and green lines represent the yield surfaces in the middle of the flow passage. § 3.2 for the simple channel flow, which is the case in the complex flows through model porous media as well. Figure 12( f ) confirms that the viscous dissipation, even in sliding flows, is at least one order of magnitude less than the plastic and slip dissipations. Interestingly, for sliding flows, a(u * , u * ) converges to its limiting value when B → ∞ and the limiting value is a function of α s only. Whereas, in the no-slip flow, the viscous dissipation increases  by increasing Bingham number but with a much smaller rate compared with the plastic dissipation. Figure 12( f ) also shows that the leading-order dissipation terms (i.e. Bj(u * ) and |σ * nt u * s | dS) scales with ∼B at the yield limit for sliding flows which is the main reason of yielding to the same scaling for G * at this limit -see (3.13). Figure 12(g) also establishes that the leading-order term in the slip dissipation is the 'plastic' slip dissipation B s j s (u * s ). Please note that although j(u * ) of the no-slip flow is not shown in figure 12(h), the trend is the same as the sliding flows, but of course with larger values as it is clear from figure 12 Under the no-slip condition, the asymptotic values of j(u * ) at B → ∞ are equal to 3.42 and 12.57 for regular and staggered geometries which yield to Od c =τ y /RĜ c ≈ 0.1645 and ≈0.2237, respectively, using (3.14) in the absence of any slip dissipation. The critical  Oldroyd number in the presence of slip is calculated using the same expression in table 1. Please note that as an approximation, one may alternatively use in which j s (u * s ) is approximated by j(u * ) as is suggested by figure 12(h), especially for small values of α s (i.e. more slippery flows). These values are reported in table 1 as well. A negligible difference between Od c calculated using (3.14) and (5.2) for α s = 0.2 shows the validity of this approximation when the 'sliding yield stress' is much smaller than the material yield stress.
We shall mention that, as it is clear from figure 12(c,e,h), Od c only depends on α s as reported in table 1 (at least for 0 β s 0.1).

Towards more complex geometries
Although previous studies on the yield-stress fluid flows in the model porous media uncovered some generic interesting aspects of the problem, yet they were not capable of capturing the whole physics behind these type of complex flows. For instance, Talon & Bauer (2013) utilized the lattice Boltzmann method (known as LBM) to study the Bingham fluid flow in a stochastically reconstructed porous media. These authors showed that due to the yield stress, channelization can happen: the material will flow only in self-selected paths depending on the imposed driving pressure gradient; at the yield limit it is only one path and by increasing the pressure gradient gradually, other paths will open up. This has been supported by the work of Liu et al. (2019) in model pore networks connected by straight tubes.
To capture these kinds of fascinating and more realistic behaviours, in this section we will complement the findings in § 5 by studying the flow in random-designed porous media with a wide range of porosities. The computational domain is a box of dimensions 50 × 50 where the fluid flows around 2-D fixed rigid circles of radius unity with symmetric boundary conditions on the pair of the horizontal and vertical box faces. Again we use the [R] formulation in this section and the same numerical method as in the previous section.
To have a better comparison, in addition to the random porous media, we simulate the flow again in the two introduced model porous media (figure 10); yet with the same porosity as the random cases: the size of the model cell in figure 10 is L = √ π/1 − φ. Figure 13 illustrates the no-slip flow in the case φ = 0.7 for a wide range of Bingham numbers from unity to 10 4 . Panels (a,b), (c) show the velocity magnitude contours and When the Bingham number increases, the flow becomes more and more localized, and at the yield limit (B → ∞), the flow only passes through a single channel. Please note that the black level contours represent the quiescent parts. It is worth mentioning that due to the [R] scaling, the flow rate in all the simulations is the same (for a fixed geometry regardless of the Bingham number, α s and β s ), hence, in figure 13(c) the whole fluid passes through the single open channel which results in extremely high velocities (compare the level of contours in the panels (a,b), (c)). By looking at the panels from left to right sequentially, the channelization characteristic is clearly observable: by increasing the Bingham number (i.e. moving to the yield limit), the flow occurs in lesser paths, until only one path remains at the yield limit. Figure 14 shows the pressure drop against the Bingham number. The lines are computed data in the present study and the symbols are computational and experimental data by Bauer et al. (2019) for three-dimensional (3-D) flows. The blue lines are acquired from the randomized porous media and the red lines from the two model porous media for a wide range of porosities. To have a better comparison, we combined G * from the regular and staggered models in the red lines by averaging. This figure demonstrates that the model porous media closely predicts the pressure drop in the randomized porous media, however, there are small discrepancies, more pronounced in the two limits of the problem -Newtonian and the yield limit. Yet, in the intermediate regime of the Bingham numbers, the predictions are more satisfactory. Table 2 compares the critical Oldroyd number under the no-slip condition computed in different configurations. As can be seen in figure 14, the model porous media slightly overpredicts the pressure drop which results in smaller Od c compared with the randomized geometries.

Baseline no-slip flow
In figure 14, we also compare 2-D simulations in the present study and 3-D simulations by Bauer et al. (2019). They are partially comparable; specifically in the large Bingham number limit where a large portion of domain is filled with static/fouling regions.  However, a closer comparison in the Newtonian limit is suggestive of different trends, which could be the consequence of different flow structures in 2-D and 3-D flows in this limit and also the difference between the viscous functionality in the Bingham model which is used here and the Herschel-Bulkley model used by Bauer et al. (2019). More sophisticated analyses/comparisons of the model and random porous media under the no-slip condition are presented in appendix B, as the main focus of the present study is on the sliding flows.
6.2. Effect of slip In this subsection, we take a close look at the effect of slip in the randomized porous media. Figures 15 and 16 compare the velocity and the rate of strain fields, respectively, between the flows under the no-slip condition and the sliding flows (α s = 0.2 & β s = 0.1) at B = 10 4 . These two figures show that when the flow can slide over the pores' surfaces, the channelization is not strong compared with the no-slip condition. Indeed, the flow slides over the solid surfaces which cover a large portion of the domain in the small porosities limit and consequently there is a flow in most of the passages. However, in the limit of high porosities (i.e. less solid surfaces), the effect of slip is negligible: compare panels (c) and  16f ), but the flow mainly passes through the open channel whose shape and location are identical to the no-slip case. Figure 17 looks at this phenomenon from the pressure drop perspective. The lines (no-slip simulations) are borrowed from figure 14 for reference. The symbols refer to sliding flows (α s = 0.2 & β s = 0.1); the blue ones collected from the simulations in the randomized porous media and the red ones from the model porous media. Please note that again we plot the average value of G * corresponding to the two model porous media (normal and staggered cases). This figure demonstrates that as the porosity increases, the effect of slip on the pressure drop decreases, in agreement with figures 15 and 16 for the model porous media. Indeed, the symbols are much closer to their corresponding lines with the same porosity as φ → 1. Moreover, in this high porosity limit, the model porous media predictions are appropriate approximations of the pressure drop in the randomized porous media: the asterisks and pluses with different colours are not distinguishable. However, in the other limit (low porosities), the effect of slip on the pressure drop is more pronounced: the pressure drop in the sliding flows are smaller compared with no-slip flows, which is intuitive. Moreover, the intensity of the slip effect on the pressure drop depends on the geometric structures since there is a significant discrepancy between the predictions of the model porous media and the randomized porous media: the blue and red full circles are relatively far from each other specifically in the low Bingham number regime (please note the logarithmic scale).

Summary and discussion
The hydrodynamic features of sliding flows of yield-stress fluids are investigated in the present study. It has been well-documented in the literature that yield-stress fluids slide over solid surfaces due to microscopic effects: formation of a lubrication layer of solvent and elastic deformation of soft particles in the vicinity of the solid surface (Meeker et al. 2004b). Firstly, we formulated a general vectorial form of the well known 'stick-slip' law and presented a numerical algorithm based on the augmented Lagrangian scheme combined with anisotropic mesh adaptation to attack these kinds of problems. We derived some theoretical tools as well to formulate the yield limit in the presence of slip. The whole framework was benchmarked in a simple channel Poiseuille flow. Then we addressed more complicated problems including moving solid surfaces and the sliding flow over complex topologies. Firstly, we simulated the slippery particle sedimentation problem and addressed the yield limit in detail. Slipline solutions from the perfectly plastic mechanics were revisited and utilized to find the yield limit in the presence of slip. The slipline method has proved capable of finding the yield limit, even though the velocity/stress fields obtained from slipline solutions can differ from the viscoplastic solutions at large, yet finite Bingham numbers -see Chaparian & Frigaard (2017b). The conclusions from the present study are similar, showing again that the lower and upper bounds of the plastic drag coefficient are excellent estimations of the yield limit of slippery particle motion.
Secondly, we addressed the complex sliding flows in the model and randomized porous media. Different hydrodynamic features of the flow were investigated as well as the critical pressure gradient to ensure continuous non-zero flow was reported. Moreover, some general conclusions were drawn about modelling the flow in porous media by comparing the results computed from the model configurations and more realistic random-designed porous media, both under the no-slip condition and sliding flows. Furthermore, it was shown that the effect of slip is more severe in the low porosity limit where a larger portion of the domain is filled by occupying solids (i.e. a large portion of the flow is in contact with the solid surfaces). In this limit, the pressure drop of a slippery flow is decreased compared with the flow under the no-slip condition. For instance, the pressure drop under the no-slip condition at φ = 0.5 closely follows the pressure drop of the sliding flow for the case φ = 0.38 when α s = 0.2 & β s = 0.1 (see figure 17). Indeed, an effective porosity can be defined for the sliding flows in porous media to predict the pressure drop: G * (φ eff ; no-slip) = G * (φ; α s , β s ) where φ eff > φ. In the high porosity limit, however, φ eff ≈ φ since the solid surfaces are occupied by a small portion of the whole domain.
By investigating different physical problems in the presence of slip, from particle motion to pressure-driven flows in complex geometries, the following general conclusions can be drawn about sliding flows of viscoplastic fluids.
(i) There are three sources of dissipation in sliding flows: viscous, plastic and slip dissipations. Slip dissipation itself can be split into two main contributions: 'viscous' slip dissipation which comes from the slip coefficient β s and 'plastic' slip dissipation which is the contribution of the sliding yield stress. (ii) In the yield limit, the leading-order contribution to the total slip dissipation is the 'plastic' slip dissipation. Exploiting that viscous dissipation is at least one order of magnitude less than the plastic dissipation, we can conclude that the yield limit of sliding flows is indeed controlled by α s or sliding yield stress. This was evidenced both for particle sedimentation and pressure-driven flows in porous media. (iii) In the other limit (Newtonian limit), however, the flow characteristics are determined by the slip coefficient β s . For instance, the plastic drag coefficient in the particle sedimentation problem and the pressure drop in the porous media are the same for the flows with the same β s as B → 0.
Finally, we should mention that, although in the entire study we considered k = 1 and the Bingham model as the rheological representation, the presented framework and most of the conclusions that are drawn are independent of these assumptions. In appendix C, we will show how we can adopt the presented numerical algorithm and the analytical tools to study the cases in which k / = 1 and n / = 1. the flow rate curves and how it differs from the simulations/experiments of Bauer et al. (2019). Moreover, we investigate if the model and randomized porous media share the same characteristics in converging to the yield limit at large B. Figure 18(a) illustrates the convergence to the yield limit (1 − Od/Od c ) versus the Bingham number. All the simulations display more or less the same scalings, where ν ≈ −1. However, two more conclusions can be drawn from figure 18(a).
(i) At small porosities, we need to go to the higher Bingham numbers to get close enough to the yield limit. This can be attributed to the highly localized and heterogeneous/channelized flow at small porosities. Data borrowed from Bauer et al. (2019) also show the same characteristics, yet a slower convergence rate (|ν| < 0.9) with the Bingham number compared with our results, which could be the consequence of 3-D flows. (ii) At high porosities, the model and randomized porous media display closely matched convergence to the yield limit by increasing the Bingham number. For instance, see the dotted blue line (randomized) and the two dotted red lines decorated with symbols (regular and staggered models with circles and squares). However, at the low porosities, the model porous media displays faster convergence compared with the randomized porous media (compare the red continuous lines and symbols (φ = 0.38) with the dashed line (φ = 0.5)).
Another interesting feature to compare, as discussed above, is the flow rate curve -see figure 18(b) which shows the flow rate versus the excessive pressure gradient. Contradictory to what has been reported by Bauer et al. (2019), this figure shows that the model porous media (at least in two dimensions), predicts same behaviour: in the sense that at small and large excessive pressure gradients (yield and Newtonian limits, respectively), the flow rate scales linearly with the excessive pressure gradient; however, at the intermediate regime, the flow rate growth as a function of 1/Od − 1/Od c is faster. This behaviour is manifested by the randomized porous media as well. Though, as is clear from figure 18(b), ζ is different in our data from the one reported by Bauer et al. (2019). Indeed, both Bauer et al. (2019) and Chevalier et al. (2013) demonstrated that ζ is close to the inverse of the power-law index of the Herschel-Bulkley fluid, ζ ≈ 1/n. Considering that in the present problem we have used the Bingham model (n = 1), it makes sense why we found a different exponent (ζ ≈ 1) which is consistent with the ζ ≈ 1/n discussion.
Appendix C. Extensions to the non-ideal cases: k / = 1 and n / = 1 In the main analyses throughout the paper, we have assumed that k = 1. As discussed in § § 1 and 2, this assumption is founded on many experimental observations. However, in this appendix, we intend to expand our understanding for the cases in which k / = 1, and also non-ideal viscoplastic rheology. In the resolution of the presented numerical method, 'simple' viscoplastic models can be easily adopted such as Herschel-Bulkley and Casson models: step 7 needs updating in the presented algorithm, which is clearly explained by Huilgol & You (2005). For instance, for the Herschel-Bulkley model, it takes the form whereΣ =T m +âγ (û m+1 ). In the case of more complex models, e.g. elastoviscoplastic models, one needs more effort -please see Chaparian & Tammisola (2019). In the following, as an example, we consider the 2-D circular Couette flow between two infinitely long cylinders of radiiR 1 andR 2 -see figure 19. We consider this specific geometry to compare our numerical results with the experimental measurements reported by Medina-Bañuelos et al. (2017): they reported the velocity profiles in the gap using rheo-particle image velocimetry measurements and also slip laws on the cylinder surfaces for a Carbopol microgel 0.12 wt. % with the rheological propertiesτ y = 27 Pa, K = 5.5 Pa.s n and n = 0.43. Figure 19 in the case of a sliding plug. As can be seen, there is a close agreement between our simulations and experimental measurements. In the case ofΩ 1 = 18 rad s −1 , the fluid is yielded close to the inner cylinder (up to ≈0.75r/d) and slides over the outer cylinder as a rotating plug, whereas in the caseΩ 1 = 0.2 rad s −1 the whole gap is unyielded and undergoes a solid-body rotation. The only change that we need to implement in algorithm 1 is step 8, ξ m+1 ← û ns · n n + û ns · t t +β s 1 + bΦ t, whereΦ = −(λ m · t)|λ m · t| k−1 + b/β s (δû m+1 · t −û ns · t).
As discussed in § 7, although the hydrodynamics features of the flows are different for different n and k, the yield limit of the problems are not dependent on these power indices. We will discuss this in the next few lines by considering a sample example in the circular Couette flow configuration and then will generalize this conclusion.
If the inner cylinder is slippery and the outer cylinder is jagged/treated to eliminate slip, then at the yield limit, the material fully slides over the inner surface and the whole gap is unyielded. Here, T denotes the applied torque on the inner cylinder per unit length. Hence, in the [R] formulation where the scale of the velocity isR 1Ω1 andd is the length scale and Considering the map between the Bingham and the yield number we can conclude that the critical yield number is Y c = R 2 1 /2πα s . Hence, the yield limit is independent of n and k. This has been demonstrated in the problems studied in the