Numerical simulations of a sphere settling in simple shear flows of yield stress fluids

Mohammad Sarabian; Marco E. Rosti; Luca Brandt; Sarah Hormozi

doi:10.1017/jfm.2020.316

Numerical simulations of a sphere settling in simple shear flows of yield stress fluids

Published online by Cambridge University Press: 01 June 2020

and

Mohammad Sarabian: Affiliation:
Department of Mechanical Engineering, Ohio University, 251 Stocker Center, Athens, OH 45701, USA
Marco E. Rosti: Affiliation:
Linné Flow Centre and SeRC (Swedish e-Science Research Centre), KTH Mechanics, SE 100 44 Stockholm, SwedenComplex Fluids and Flows Unit, Okinawa Institute of Science and Technology Graduate University, 1919-1 Tancha, Onnason, Okinawa904-0495, Japan
Luca Brandt: Affiliation:
Linné Flow Centre and SeRC (Swedish e-Science Research Centre), KTH Mechanics, SE 100 44 Stockholm, Sweden
Sarah Hormozi*: Affiliation:
Department of Mechanical Engineering, Ohio University, 251 Stocker Center, Athens, OH 45701, USA
*: † Email address for correspondence: hormozi@ohio.edu

Article contents

Abstract
Introduction
Problem definition
Results and discussion
Conclusion
References

Rights & Permissions

Abstract

We perform three-dimensional numerical simulations to investigate the sedimentation of a single sphere in the absence and presence of a simple cross-shear flow in a yield stress fluid with weak inertia. In our simulations, the settling flow is considered to be the primary flow, whereas the linear cross-shear flow is a secondary flow with amplitude 10 % of the primary flow. To study the effects of elasticity and plasticity of the carrying fluid on the sphere drag as well as the flow dynamics, the fluid is modelled using the elastoviscoplastic constitutive laws proposed by Saramito (J. Non-Newtonian Fluid Mech., vol. 158 (1–3), 2009, pp. 154–161). The extra non-Newtonian stress tensor is fully coupled with the flow equation and the solid particle is represented by an immersed boundary method. Our results show that the fore–aft asymmetry in the velocity is less pronounced and the negative wake disappears when a linear cross-shear flow is applied. We find that the drag on a sphere settling in a sheared yield stress fluid is reduced significantly compared to an otherwise quiescent fluid. More importantly, the sphere drag in the presence of a secondary cross-shear flow cannot be derived from the pure sedimentation drag law owing to the nonlinear coupling between the simple shear flow and the uniform flow. Finally, we show that the drag on the sphere settling in a sheared yield stress fluid is reduced at higher material elasticity mainly due to the form and viscous drag reduction.

Information

Type: JFM Papers
Information: Journal of Fluid Mechanics , Volume 896 , 10 August 2020, A17
DOI: https://doi.org/10.1017/jfm.2020.316 [Opens in a new window]
Copyright: © The Author(s), 2020. Published by Cambridge University Press

1 Introduction

Suspensions of dense particles in a yield stress fluid are ubiquitous in many engineering processes (crude oil, foodstuff transport, cosmetics, microfluidics, mineral slurries, cement pastes, fermentation processes, 3-D printing, drilling muds, etc.), natural phenomena (debris flows, lava flows, natural muds, etc.) and biological systems (physiology, biolocomotion, tissue engineering, etc.). Yield stress fluids exhibit a yield stress, beyond which they deform as a non-Newtonian viscous liquid while they behave as solids for lower stress levels. In mentioned practical applications, one is usually dealing with the transport of the suspended particles in yield stress fluids. Hence, in these systems, particle sedimentation occurs in the presence of a shear flow. In Stokes flow of Newtonian fluids, the governing equations of motion are linear thus the settling rate of a spherical particle is not affected by a shear flow superimposed on the settling flow. However, in yield stress fluids, the sphere drag is affected by a cross-shear flow as the linearity of the Stokes equations breaks due to the nonlinear rheology of the yield stress fluids. The objective of this paper is to investigate the flow dynamics and drag laws of a sphere settling in yield stress fluids when a cross-shear flow exists.

Sole plastic properties of yield stress fluids can be addressed via adopting ideal yield stress constitutive laws such as Bingham and Herschel–Bulkley models (also known as viscoplastic models). However, experimental works show that elasticity plays an important role in problems involving inclusions in practical yield stress fluids (Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Firouznia et al. Reference Firouznia, Metzger, Ovarlez and Hormozi2018). Therefore, the goal of the present study is twofold: we do not only study the nonlinear coupling of a settling and shear flow around a sphere but also take into account the plastic as well as the elastic properties characterising many practical yield stress fluids. Hence, in this introduction, we first give a review of several works that have been undertaken on flows around a three-dimensional (3-D) or a two-dimensional (2-D) particle immersed in a viscoplastic fluid (see also Maleki et al. Reference Maleki, Hormozi, Roustaei and Frigaard2015). Second, we review the works addressing the role of elasticity in a yield stress fluid and a pure sedimentation problem. Third, we review a limited recent number of works that investigate shear-induced sedimentation of particles in viscoelastic fluids.

The problem of a single spherical particle settling in an ideal yield stress fluid has been extensively studied theoretically and numerically (see e.g. Andres Reference Andres1960; Yoshioka, Adachi & Ishimura Reference Yoshioka, Adachi and Ishimura1971; Beris et al. Reference Beris, Tsamopoulos, Armstrong and Brown1985; Beaulne & Mitsoulis Reference Beaulne and Mitsoulis1997; Blackery & Mitsoulis Reference Blackery and Mitsoulis1997; Liu, Muller & Denn Reference Liu, Muller and Denn2002; Yu & Wachs Reference Yu and Wachs2007). Two main numerical difficulties arise in solving this problem: handling a freely moving particle, and the numerical treatment of the yield stress constitutive equations as the effective viscosity becomes infinite at the yield surface and within a solid region.

To eliminate the first difficulty, this problem has primarily been studied by fixing the particle and imposing an external flow while the confining walls are translating with the same velocity of the medium. This is called the resistance problem (i.e. the flow past a fixed obstacle). Most of these simulations were performed in two dimensions, either by assuming the obstacle to be an infinitely long circular cylinder (see e.g. Zisis & Mitsoulis Reference Zisis and Mitsoulis2002; De Besses, Magnin & Jay Reference De Besses, Magnin and Jay2003; Roquet & Saramito Reference Roquet and Saramito2003; Mitsoulis Reference Mitsoulis2004; Tokpavi, Magnin & Jay Reference Tokpavi, Magnin and Jay2008; Putz & Frigaard Reference Putz and Frigaard2010; Wachs & Frigaard Reference Wachs and Frigaard2016; Chaparian & Frigaard Reference Chaparian and Frigaard2017a,Reference Chaparian and Frigaardb; Ouattara et al. Reference Ouattara, Jay, Blésès and Magnin2018) or a spherical particle by imposing 2-D axisymmetric boundary conditions (e.g. Beris et al. Reference Beris, Tsamopoulos, Armstrong and Brown1985; Beaulne & Mitsoulis Reference Beaulne and Mitsoulis1997; Blackery & Mitsoulis Reference Blackery and Mitsoulis1997; Liu et al. Reference Liu, Muller and Denn2002). To our knowledge, only two studies exist that solve the problem of a particle settling in viscoplastic fluids (Yu & Wachs Reference Yu and Wachs2007) and creeping flow of Bingham plastics around translating objects (Sverdrup, Almgren & Nikiforakis Reference Sverdrup, Almgren and Nikiforakis2019) using fully 3-D numerical simulations. Recently, the flow past a rotating sphere in a Bingham plastic fluid has been investigated by Pantokratoras (Reference Pantokratoras2018) for a wide range of material properties and Reynolds number.

Two different approaches are implemented to remove the second difficulty: the regularization method and the augmented-Lagrangian method. The former removes the inherent discontinuity of an ideal yield stress constitutive equation (see e.g. Bingham Reference Bingham1922; Herschel & Bulkley Reference Herschel and Bulkley1926) by approximating the solid region as a fluid with an extremely high viscosity. Frigaard & Nouar (Reference Frigaard and Nouar2005) provides a review on different regularization schemes. The augmented-Lagrangian scheme, on the other hand, maintains the true constitutive equation but requires the solution of an often expensive minimization problem. For more details on this method the reader is referred to Glowinski & Wachs (Reference Glowinski and Wachs2011).

When a sphere is falling in a yield stress fluid, the stresses decay as we move away from the surface of the particle and they may become smaller than the yield stress. Therefore, there exists a rigid envelope (a plug or an unyielded region) surrounding the liquid (yielded) zone around the particle (Volarovich & Gutkin Reference Volarovich and Gutkin1953; Ansley & Smith Reference Ansley and Smith1967). Moreover, two unyielded triangular polar caps form at the front and rear stagnation points of the sphere (Beris et al. Reference Beris, Tsamopoulos, Armstrong and Brown1985; Beaulne & Mitsoulis Reference Beaulne and Mitsoulis1997; Blackery & Mitsoulis Reference Blackery and Mitsoulis1997). It is noteworthy to mention that in the literature different terms have been used to refer to the location of the polar caps, such as the leading and trailing parts of the spheres, or north and south poles of the sphere. For the 2-D cases (i.e. cylinders in an infinite domain) it has been shown that, in addition to the rigid envelope and the stagnant regions attached to the front and rear of the cylinder, two counter-rotating solid islands may form at both sides of the cylinder’s equator in the fluid zone (see e.g. De Besses et al. Reference De Besses, Magnin and Jay2003; Chaparian & Frigaard Reference Chaparian and Frigaard2017b). Ansley & Smith (Reference Ansley and Smith1967) postulated the shape of the plug regions around a falling sphere and pointed to plasticity theory as a tool for tackling problems involving yield stress fluids. However, it is only more recently that researchers have made systematic use of plasticity theory and especially the slipline method for analysing the yield limit in 2-D problems involving yield stress fluids (Liu et al. Reference Liu, Balmforth, Hormozi and Hewitt2016; Balmforth et al. Reference Balmforth, Craster, Hewitt, Hormozi and Maleki2017; Chaparian & Frigaard Reference Chaparian and Frigaard2017a,Reference Chaparian and Frigaardb).

Wall effects on the evolution of the yielded/unyielded zones has been investigated numerically and experimentally. These studies include either a particle settling in a tube filled with a viscoplastic fluid (Beaulne & Mitsoulis Reference Beaulne and Mitsoulis1997; Blackery & Mitsoulis Reference Blackery and Mitsoulis1997) or the flow of a viscoplastic fluid past a 2-D cylinder (a long circular cylinder) located between two parallel plates (see e.g. Zisis & Mitsoulis Reference Zisis and Mitsoulis2002; Mitsoulis Reference Mitsoulis2004; Ouattara et al. Reference Ouattara, Jay, Blésès and Magnin2018). Generally, it is found that at a fixed confinement ratio, the yielded zone surrounding the particle shrinks with increasing yield stress leaving thin viscous layers around the particle. A well-resolved series of simulations for the case of a settling circular disk shows that the thin viscous layers resemble a cross-eyed owl (Wachs & Frigaard Reference Wachs and Frigaard2016). At fixed yield stress, increasing the confinement ratio results in the extension of the yielded region around the particle and, eventually, its interaction with the walls.

The drag exerted on a particle settling in an ideal yield stress fluid is found to be a function of both the Bingham number and the confinement ratio. The sphere drag is an increasing function of Bingham number ($Bi$) at fixed confinement ratio (see e.g. Atapattu, Chhabra & Uhlherr Reference Atapattu, Chhabra and Uhlherr1995; Blackery & Mitsoulis Reference Blackery and Mitsoulis1997; Tabuteau, Coussot & de Bruyn Reference Tabuteau, Coussot and de Bruyn2007; Ahonguio, Jossic & Magnin Reference Ahonguio, Jossic and Magnin2014). Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997) reported the Stokes drag coefficient (i.e. the ratio of the total drag force to the Stokes drag on a sphere settling in a Newtonian fluid) as a function of the Bingham number (i.e. the ratio of the yield stress to viscous stresses) and showed an enhancement of the Stokes drag coefficient with an increasing confinement ratio. For large Bingham numbers $(Bi\geqslant 10)$, however, the Stokes drag coefficient is independent of the confinement ratio. At large enough $Bi$ numbers, the envelope of the yielded zone around a particle is encapsulated by an outer plug region attached to the boundaries. Thus, at this point, a further reduction of the confinement does not affect the flow dynamics and drag forces. These results are qualitatively valid for cases of viscoplastic fluids flowing around a 2-D cylinder in a duct (Zisis & Mitsoulis Reference Zisis and Mitsoulis2002; Mitsoulis Reference Mitsoulis2004). Wall slip also alters the flow dynamics and consequently the drag forces in flows of yield stress fluids over solid boundaries. In experimenting with practical viscoplastic materials, the slip at the wall is unavoidable (see e.g. Meeker, Bonnecaze & Cloitre Reference Meeker, Bonnecaze and Cloitre2004; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012). De Besses et al. (Reference De Besses, Magnin and Jay2003) showed that the wall slip reduces the drag force exerted by a viscoplastic fluid on a 2-D cylinder in an infinite medium.

It is of practical interest to know the force required to move a body immersed in a yield stress fluid. We refer the reader to the early and historical work of Boardman & Whitmore (Reference Boardman and Whitmore1960). In sedimentation flows, a particle does not fall when the yield stress resistance is larger than the buoyancy stress. A dimensionless number called a gravity number $Y_{G}$ is defined as the ratio between material yield stress and the buoyancy stress. There exists a critical value of $Y_{G}$ beyond which the particle ceases to move and it is entrapped within the yield stress fluid. The critical value of $Y_{G}$ depends on the geometry of the immersed object (Ovarlez & Hormozi Reference Ovarlez and Hormozi2019) and it has been found numerically for a 2-D circle (Tokpavi et al. Reference Tokpavi, Magnin and Jay2008), a sphere (Beris et al. Reference Beris, Tsamopoulos, Armstrong and Brown1985), a 2-D square (Nirmalkar, Chhabra & Poole Reference Nirmalkar, Chhabra and Poole2012) and a 2-D ellipse (Putz et al. Reference Putz, Burghelea, Frigaard and Martinez2008). Experimental attempts also have been made to determine the critical value of $Y_{G}$ for objects with different shapes (see e.g. Jossic & Magnin Reference Jossic and Magnin2001; Tabuteau et al. Reference Tabuteau, Coussot and de Bruyn2007; Tokpavi et al. Reference Tokpavi, Jay, Magnin and Jossic2009; Ahonguio, Jossic & Magnin Reference Ahonguio, Jossic and Magnin2015). There exists a quantitative disagreement between the numerical results and experimental results due to the unavoidable wall-slip effects in experiments as well as owing to the discrepancy between the rheological behaviour of the practical yield stress fluids and the ideal yield stress laws used in the theoretical and computational work. To this end, one of the objectives of the present work is to conduct a numerical study of the flows of yield stress fluids around a sphere using a more realistic constitutive law, i.e. including some fluid elasticity.

In many soft materials (e.g. emulsions, foams, gels, colloidal pastes, etc.), elasticity is found to play an important role in the dynamics of the flow around a single sphere or 2-D cylinder (see e.g. Atapattu et al. Reference Atapattu, Chhabra and Uhlherr1995; Gueslin et al. Reference Gueslin, Talini, Herzhaft, Peysson and Allain2006; Dollet & Graner Reference Dollet and Graner2007; Tabuteau et al. Reference Tabuteau, Coussot and de Bruyn2007; Putz et al. Reference Putz, Burghelea, Frigaard and Martinez2008; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Ahonguio et al. Reference Ahonguio, Jossic and Magnin2014; Ouattara et al. Reference Ouattara, Jay, Blésès and Magnin2018), particles of various shapes (Jossic & Magnin Reference Jossic and Magnin2001), dilute suspensions (Hariharaputhiran et al. Reference Hariharaputhiran, Subramanian, Campbell and Chhabra1998), concentrated suspensions (Coussot et al. Reference Coussot, Raynaud, Bertrand, Moucheront, Guilbaud, Huynh, Jarny and Lesueur2002) and bubbles rising in yield stress fluids (see e.g. Sikorski, Tabuteau & de Bruyn Reference Sikorski, Tabuteau and de Bruyn2009; Fraggedakis et al. Reference Fraggedakis, Pavlidis, Dimakopoulos and Tsamopoulos2016b). Thus, to accurately predict the behaviour of practical viscoplastic fluids, one must also take the elasticity effects into account. These materials are called elastoviscoplastic (EVP). There are a number of different constitutive equations proposed in the literature to model EVP fluids (e.g. Saramito Reference Saramito2007; de Souza Mendes Reference de Souza Mendes2007; Bénito et al. Reference Bénito, Bruneau, Colin, Gay and Molino2008; Saramito Reference Saramito2009; Dimitriou, Ewoldt & McKinley Reference Dimitriou, Ewoldt and McKinley2013; Dimitriou & McKinley Reference Dimitriou and McKinley2014; Geri et al. Reference Geri, Venkatesan, Sambath and McKinley2017). For more information on the details of these models and their implementation, the reader is referred to Izbassarov et al. (Reference Izbassarov, Rosti, Ardekani, Sarabian, Hormozi, Brandt and Tammisola2018). Several modelling and numerical works have implemented the constitutive laws of EVP materials to quantitatively capture the flow characteristics of practical yield stress fluids (e.g. Cheddadi et al. Reference Cheddadi, Saramito, Dollet, Raufaste and Graner2011; Cheddadi & Saramito Reference Cheddadi and Saramito2013; Fraggedakis, Dimakopoulos & Tsamopoulos Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a).

The problem of a particle settling in an EVP fluid confined in a cylindrical container is studied through axisymmetric finite element computations by Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a). The numerical results show that the loss of fore–aft symmetry in the velocity field around the sphere and the appearance of the negative wake downstream of the sphere settling in a yield stress fluid are due to the elasticity. These phenomena were also observed experimentally (see e.g. Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012). Furthermore, it has been demonstrated that the extent and shape of the yielded/unyielded regions, the particle stoppage criterion and the sphere drag are influenced by the presence of elasticity in laboratory yield stress fluids (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a). The flow field around a single sphere and the trajectories of two spheres in simple shear flows of practical yield stress fluids (Carbopol gel) were also studied experimentally showing the importance of elastic effects (Firouznia et al. Reference Firouznia, Metzger, Ovarlez and Hormozi2018).

There exists a handful experimental studies focusing on the shear-induced sedimentation of suspensions of particles in practical yield stress fluids (see e.g. Merkak, Jossic & Magnin Reference Merkak, Jossic and Magnin2009; Ovarlez, Barral & Coussot Reference Ovarlez, Barral and Coussot2010; Ovarlez et al. Reference Ovarlez, Bertrand, Coussot and Chateau2012). Merkak et al. (Reference Merkak, Jossic and Magnin2009) studied the shear-induced sedimentation of suspensions of particles in a pipe flow configuration and proposed the criterion for particle stability in the sheared yield stress fluid. Ovarlez et al. (Reference Ovarlez, Bertrand, Coussot and Chateau2012) examined the settling rate of suspensions of particles in different yield stress fluids (concentrated emulsions and Carbopol gels). The density mismatched suspensions are sheared in a wide-gap Couette device; the velocity profile and solid volume fraction are measured by employing magnetic resonance imaging techniques. It is found that, for all the systems, the particles that were stable at rest start sedimenting in the yielded zone as the shear is introduced. Moreover, the experimental results of Ovarlez et al. (Reference Ovarlez, Bertrand, Coussot and Chateau2012) show that the particle settling velocity scaled by the Stokes velocity is an increasing function of the inverse Bingham number. Despite all the experimental results falling into a master curve in the limits of high Bingham number (plastic regime), there exists a discrepancy when $Bi\leqslant 1$, which might be due to another governing dimensionless number associated with the presence of elasticity in a practical yield stress fluid.

A larger number of studies focus on the effect of shear flow on the particle settling rate in viscoelastic fluids through experimental measurements (see e.g. van den Brule & Gheissary Reference van den Brule and Gheissary1993; Murch et al. Reference Murch, Krishnan, Shaqfeh and Iaccarino2017), direct numerical simulations (see e.g. Padhy et al. Reference Padhy, Shaqfeh, Iaccarino, Morris and Tonmukayakul2013b,Reference Padhy, Rodriguez, Shaqfeh, Iaccarino, Morris and Tonmukayakula; Murch et al. Reference Murch, Krishnan, Shaqfeh and Iaccarino2017) and theoretical analysis (see e.g. Brunn Reference Brunn1977b,Reference Brunna; Housiadas & Tanner Reference Housiadas and Tanner2012; Vishnampet & Saintillan Reference Vishnampet and Saintillan2012; Housiadas Reference Housiadas2014; Housiadas & Tanner Reference Housiadas and Tanner2014; Einarsson & Mehlig Reference Einarsson and Mehlig2017).

The drag enhancement on the sphere when it is subjected to an externally imposed cross-shear flow is verified through direct numerical simulations both in a weakly (Padhy et al. Reference Padhy, Shaqfeh, Iaccarino, Morris and Tonmukayakul2013b) and in a highly (Murch et al. Reference Murch, Krishnan, Shaqfeh and Iaccarino2017) constant-viscosity viscoelastic fluid. In the case of weakly viscoelastic fluids, it has been shown that adding polymers indirectly increases the drag on the sphere through breaking the symmetry in the polymer and viscous stresses on the sphere surface (Padhy et al. Reference Padhy, Shaqfeh, Iaccarino, Morris and Tonmukayakul2013b). Contrary to the weakly viscoelastic fluids, the form drag is a primary cause of the drag enhancement on the particle in the limit of highly viscoelastic fluids when the cross-shear flow is coupled with the uniform flow (Murch et al. Reference Murch, Krishnan, Shaqfeh and Iaccarino2017). In addition to the elasticity effects, the sphere drag coefficient is also affected by the rheological properties of the viscoelastic fluid such as the degree of shear thinning and the ratio of polymer to total viscosity (Padhy et al. Reference Padhy, Rodriguez, Shaqfeh, Iaccarino, Morris and Tonmukayakul2013a).

The main novel contributions of our study can be summarised as follows. First, we validate a new tool for the 3-D numerical simulations of a sphere settling in a quiescent and sheared yield stress fluid at small particle Reynolds number. We formulate the problem using both an ideal yield stress fluid law (Bingham model) and EVP fluid law. The governing equations are solved using our recently developed 3-D numerical solver to handle the problems including rigid and deformable particles suspended in complex fluids (Izbassarov et al. Reference Izbassarov, Rosti, Ardekani, Sarabian, Hormozi, Brandt and Tammisola2018). The particle motion in the flow is simulated by means of an efficient immersed boundary (IB) method coupled with a flow solver for the equations of motion associated with the suspending complex fluids. The mechanical and mathematical modelling, the computational matrix, the boundary conditions and the numerical scheme of the simulations presented here are reported in § 2. Second, we provide for the first time a 3-D and extensive analysis of the velocity, stress fields and yield surfaces around a sphere settling in yield stress fluids. We consider both quiescent and sheared yield stress fluids and compute the drag force and its individual components, i.e. form drag, viscous drag, polymer drag and inertia drag and discuss their dependency on the elastic and plastic properties of the suspending fluids as well as the flow dynamics. Third, we show that the flow dynamics and drag forces significantly change when a secondary linear cross-shear flow is superimposed. In particular, the nonlinear coupling of the settling and simple shear flow plays a significant role in determining the flow dynamics and consequently the drag force on the sphere, which makes the prediction of sphere drag in a sheared yield stress fluid difficult (§ 3). The main conclusions of this work are summarized in § 4.

Figure 1.

Computational domain and boundary conditions.

2 Problem definition

We consider the sedimentation of a single spherical particle subjected to a linear cross-shear flow in an EVP fluid with slight inertia. The reference frame is attached to the particle with its origin at the particle centre. Thus, the particle remains stationary and the fluid moves with uniform velocity $U_{\infty }$ in the opposite direction to gravity (figure 1). The top and bottom walls are translating with constant and equal velocity, but in the opposite direction to create the linear shear flow in the $x$–$z$ plane. Moreover, the walls move with the uniform velocity $U_{\infty }$ in the streamwise direction $y$. Consequently, the upper and lower plates are moving obliquely. While remaining fixed, the particle rotates freely due to the applied cross-shear flow.

2.1 Mechanical model of EVP material

In this paper, we have adopted the model proposed by Saramito (Reference Saramito2009) to simulate the EVP material. Figure 2 represents a mechanical description of this model, based on the friction $\unicode[STIX]{x1D70F}_{0}$, spring stiffness $G$ and solvent $\unicode[STIX]{x1D702}_{s}$ and polymer $\unicode[STIX]{x1D702}_{p}$ viscosities. When the elastic strain energy of the system exceeds the threshold value (yield value), the friction element breaks allowing viscous deformation. After yielding, the deformation is unbounded in time, thus the system behaves as a nonlinear viscoelastic fluid. Before yielding, however, the material behaves as a nonlinear viscoelastic solid.

Figure 2.

Mechanical model schematic of an EVP material proposed by Saramito (Reference Saramito2009).

2.2 Mathematical model

This section presents the governing equations in dimensionless form. Here, the characteristic velocity scale is the uniform inflow velocity $U_{\infty }$, the characteristic length scale is the particle diameter $D$, the characteristic time scale is $D/U_{\infty }$ (since the uniform flow is the main flow) and $\unicode[STIX]{x1D702}_{0}U_{\infty }/D$ is the characteristic stress. Here, $\unicode[STIX]{x1D702}_{0}$ is the total viscosity at zero shear-rate which is sum of the solvent $\unicode[STIX]{x1D702}_{s}$ and polymer $\unicode[STIX]{x1D702}_{p}$ viscosities $\unicode[STIX]{x1D702}_{0}=\unicode[STIX]{x1D702}_{s}+\unicode[STIX]{x1D702}_{p}$. The continuity and momentum equations are as follows:

(2.1)

$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.2)

$$\begin{eqnarray}\displaystyle & \displaystyle Re_{p}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\unicode[STIX]{x1D735}\boldsymbol{u}\right)=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(-pI+2(1-\unicode[STIX]{x1D6FD})D(\boldsymbol{u})+\unicode[STIX]{x1D749})+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}$ is the fluid velocity vector and $D(\boldsymbol{u})$ is the rate of deformation tensor defined as $D(\boldsymbol{u})=1/2(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})$. It is noteworthy to mention that, we subtract the static pressure of the fluid and scale the remainder with the viscous stress scale. Therefore, the dimensionless pressure $p$ is the dimensionless dynamic pressure. Here, $Re_{p}$ is the particle Reynolds number that is defined as the ratio of the inertial forces to viscous forces: $Re_{p}=(\unicode[STIX]{x1D70C}_{f}U_{\infty }D)/\unicode[STIX]{x1D702}_{0}$, where $\unicode[STIX]{x1D70C}_{f}$ is the fluid density. The retardation parameter, $\unicode[STIX]{x1D6FD}$, is the ratio between the polymer and the total viscosities $\unicode[STIX]{x1D702}_{p}/\unicode[STIX]{x1D702}_{0}$, $\boldsymbol{f}$ is an external body force (IB force) used to model the presence of the particle. Here, $\unicode[STIX]{x1D749}$ is the extra stress tensor due to the plasticity and elasticity of the EVP fluid.

The following constitutive equation is proposed to model the stress-deformation of EVP materials (Saramito Reference Saramito2009):

(2.3)

$$\begin{eqnarray}Wi\overset{\unicode[STIX]{x1D735}}{\unicode[STIX]{x1D749}}+\unicode[STIX]{x1D705}_{n}(|\unicode[STIX]{x1D70F}_{d}|)\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D749}-2\unicode[STIX]{x1D6FD}D(\boldsymbol{u})=0,\end{eqnarray}$$

where $Wi$ is the shear Weissenberg number, the ratio of the material relaxation time scale $\unicode[STIX]{x1D706}$ to the cross-shear flow time scale $1/\dot{\unicode[STIX]{x1D6FE}_{0}}$. Therefore, $Wi=\unicode[STIX]{x1D706}\dot{\unicode[STIX]{x1D6FE}_{0}}$. Note that the material relaxation time $\unicode[STIX]{x1D706}$ is the ratio of polymer viscosity $\unicode[STIX]{x1D702}_{p}$ to the elasticity parameter $G$: $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D702}_{p}/G$. Here, $\unicode[STIX]{x1D6FC}$ is another dimensionless parameter whose magnitude determines the primary and secondary flow; in particular, the ratio of the applied shear rate $\dot{\unicode[STIX]{x1D6FE}_{0}}$ to the shear rate induced by particle settling in the fluid $\dot{\unicode[STIX]{x1D6FE}}_{sett}$ is denoted by $\unicode[STIX]{x1D6FC}=\dot{\unicode[STIX]{x1D6FE}_{0}}/\dot{\unicode[STIX]{x1D6FE}}_{sett}$. The settling shear rate is the ratio of the uniform inflow velocity to the particle diameter, i.e. $\dot{\unicode[STIX]{x1D6FE}}_{sett}=U_{\infty }/D$ and thus $\unicode[STIX]{x1D6FC}=\dot{\unicode[STIX]{x1D6FE}_{0}}D/U_{\infty }$. In the context of EVP fluids, $\unicode[STIX]{x1D6FC}$ can be viewed as the ratio between the shear Weissenberg number $Wi$ and the settling Weissenberg number, which is called ‘$Wi_{\infty }$’ in the present work and is equal to $(\unicode[STIX]{x1D706}U_{\infty })/D$.

The upper convected derivative of the stress field is $\overset{\unicode[STIX]{x1D735}}{\unicode[STIX]{x1D749}}$; it is computed based on the following relation:

(2.4)

$$\begin{eqnarray}\overset{\unicode[STIX]{x1D735}}{\unicode[STIX]{x1D749}}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D749}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D749}-(\unicode[STIX]{x1D735}\boldsymbol{u})^{\text{T}}\boldsymbol{\cdot }\unicode[STIX]{x1D749}-\unicode[STIX]{x1D749}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}.\end{eqnarray}$$

In (2.3) $\unicode[STIX]{x1D705}_{n}(|\unicode[STIX]{x1D70F}_{d}|)$ is the plasticity criterion, defined by the following relation:

(2.5)

$$\begin{eqnarray}\unicode[STIX]{x1D705}_{n}(|\unicode[STIX]{x1D70F}_{d}|)=\max \left(0,\frac{|\unicode[STIX]{x1D70F}_{d}|-Bi}{(2\unicode[STIX]{x1D6FD})^{1-n}|\unicode[STIX]{x1D70F}_{d}|^{n}}\right)^{1/n},\end{eqnarray}$$

where $|\unicode[STIX]{x1D70F}_{d}|$ is the second invariant of the deviatoric part of the extra stress tensor,

(2.6)

$$\begin{eqnarray}\unicode[STIX]{x1D749}_{\boldsymbol{d}}=\unicode[STIX]{x1D749}-{\textstyle \frac{1}{3}}\text{tr}(\unicode[STIX]{x1D749})\unicode[STIX]{x1D644},\end{eqnarray}$$

with $\unicode[STIX]{x1D644}$ the unit tensor and $n$ the power-law index. The Bingham number is $Bi$, which is the ratio of yield stress to viscous stress, $Bi=(\unicode[STIX]{x1D70F}_{0}D)/(\unicode[STIX]{x1D702}_{0}U_{\infty })$. Here, $\unicode[STIX]{x1D70F}_{0}$ is the material yield stress. Note that throughout the paper $Bi$ is the Bingham number defined based on the particle settling shear rate. Polymer viscosity is computed via $\unicode[STIX]{x1D702}_{p}=\unicode[STIX]{x1D705}(U_{\infty }/D)^{n-1}$, where $\unicode[STIX]{x1D705}$ is the consistency parameter.

The rotational particle velocity is computed by solving the Euler equation in the body-fixed reference frame:

(2.7)

$$\begin{eqnarray}I_{s}\frac{\text{d}\unicode[STIX]{x1D714}_{c}}{\text{d}t}=\oint _{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}\boldsymbol{r}\times (\unicode[STIX]{x1D70F}\boldsymbol{\cdot }\boldsymbol{n})\,\text{d}\,A.\end{eqnarray}$$

In (2.7), $I_{s}$ is the moment of inertia of the particle and is equal to $2/5\unicode[STIX]{x1D70C}_{s}V_{s}R^{2}$, where $\unicode[STIX]{x1D70C}_{s}$ is the particle density, $V_{s}$ is the particle volume and $R$ is the particle radius. Here, $\boldsymbol{n}$ is the unit normal vector at the particle surface $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ and $\unicode[STIX]{x1D714}_{c}$ is the particle rotational velocity. Since our simulations are performed in a body-fixed reference frame, the particle angular velocity is tracked in an inertial reference frame by adopting the rotation matrix from Chen et al. (Reference Chen, Phan-Thien, Khoo and Fan2006). For more details on the transformation the reader is referred to Vázquez-Quesada, Bian & Ellero (Reference Vázquez-Quesada, Bian and Ellero2016).

2.3 Boundary conditions

The computational domain is a rectangular box with length $16D$ in the streamwise $y$ direction, $6D$ in the spanwise $x$ direction and $12D$ in the wall-normal $z$ direction. The particle is located in the middle of the box with its centre as the origin of the coordinate system. The inlet velocity boundary condition is a combination of uniform velocity and the simple shear flow. Hence, at the inlet, the fluid velocity has a component in both the streamwise and spanwise directions. The convective outflow boundary condition is applied at the outlet where the fluid is leaving the domain in the $y$ direction. The velocity is extrapolated by solving the convective equation at the location of the outlet (Uhlmann Reference Uhlmann2003):

(2.8)

$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+U_{\infty }\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}\boldsymbol{n}}=0.\end{eqnarray}$$

In order to maintain the overall mass balance and satisfy the compatibility condition, the convective velocity should be equal to the uniform inflow velocity $U_{\infty }$. A no-slip boundary condition is applied at the top and bottom plates with the velocities of the plates equal to $\boldsymbol{u}=-\dot{\unicode[STIX]{x1D6FE}}_{0}L_{z}\hat{x}+U_{\infty }{\hat{y}}$ and $\boldsymbol{u}=\dot{\unicode[STIX]{x1D6FE}}_{0}L_{z}\hat{x}+U_{\infty }{\hat{y}}$. A homogeneous Neumann condition is used for pressure at the two walls as well as the inlet and outlet ($\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}\boldsymbol{n}=0$). A no-flux condition is specified normal to the confinement walls for the extra stress tensor. A periodic boundary condition is applied for the velocity, pressure and extra stress tensor in the spanwise $x$ direction. The no-slip/no-penetration boundary condition is satisfied at the sphere surface implicitly by using the multidirect forcing immersed boundary scheme (Breugem Reference Breugem2012).

We analytically solve the EVP constitutive equations of Saramito (Reference Saramito2009) for the combination of Couette and uniform flow at steady state in the absence of the spherical particle and then apply this solution for the EVP stress tensor at the inlet. At steady state the flow of EVP fluid is entirely yielded. Thus, the plasticity criteria function is positive and non-zero ($\unicode[STIX]{x1D705}_{n}(|\unicode[STIX]{x1D70F}_{d}|)>0$). In this flow configuration, all the components of the extra stress tensor are zero except the normal and shear stress in the spanwise direction of the shear plane $x$–$z$. The non-vanishing value of the normal stress difference, which is observed in practical yield stress fluids (e.g. Janiaud & Graner Reference Janiaud and Graner2005), is the result of incorporating the elastic stresses in the EVP constitutive equation as the classical viscoplastic models (such as Bingham Reference Bingham1922 and Herschel & Bulkley Reference Herschel and Bulkley1926) do not predict the first normal stress difference in shear flows. The analytical solution is obtained by solving the following coupled nonlinear equations for the normal and shear EVP stresses:

(2.9)

$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{|\unicode[STIX]{x1D70F}_{d}|-Bi}{|\unicode[STIX]{x1D70F}_{d}|}\right)\unicode[STIX]{x1D70F}_{xz}=2\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}, & \displaystyle\end{eqnarray}$$

(2.10)

$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{xz}=\sqrt{\frac{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}{2Wi}\unicode[STIX]{x1D70F}_{xx}}, & \displaystyle\end{eqnarray}$$

where $|\unicode[STIX]{x1D70F}_{d}|=\sqrt{1/3\unicode[STIX]{x1D70F}_{xx}^{2}+\unicode[STIX]{x1D70F}_{xz}^{2}}$ as $\unicode[STIX]{x1D70F}_{xx}$ and $\unicode[STIX]{x1D70F}_{xz}$ are the only non-zero elements of the stress tensor. The nonlinear system of (2.9) and (2.10) is easily solved using the ‘fsolve’ routine of MATLAB. For the pure sedimentation simulations, the stress tensor at the inlet and outlet is set equal to zero.

2.4 Computational matrix

In this study, all the simulations are conducted at $\unicode[STIX]{x1D6FC}=0.1$, i.e. the particle settling shear rate is ten times larger than the externally imposed shear rate. We will present a series of simulations for the problem of single sphere settling through EVP fluid in the absence and presence of simple cross-shear flow. The dimensionless parameters used for the simulations presented here are reported in table 1.

Table 1.

Computational matrix.

2.5 Numerical method

The comprehensive details of the numerical algorithm are explained in Izbassarov et al. (Reference Izbassarov, Rosti, Ardekani, Sarabian, Hormozi, Brandt and Tammisola2018). In brief, the equations are solved on a Cartesian, staggered, uniform grid with velocities located on the cell faces and all the other variables (pressure, stress and material component properties) at the cell centres. All the spatial derivatives are approximated with second-order centred finite differences except for the advection terms in (2.4) where the fifth-order weighted essentially non-oscillatory (known as WENO) scheme is adopted (Shu Reference Shu2009). The time integration is performed with a fractional-step method (Kim & Moin Reference Kim and Moin1985), where all the terms in the evolution equations are advanced in time with a third-order explicit Runge–Kutta scheme except for the EVP stress terms which are advanced with the Crank–Nicolson method. Moreover, a fast Poisson solver is used to enforce the condition of zero divergence for the velocity field. The coupling of the fluid and particles is performed with the immersed boundary method proposed by Breugem (Reference Breugem2012).

We have adopted a grid resolution of $32$ Eulerian grid points for particle diameter. Each simulation is performed on $128$ cores working in parallel and the steady-state solution is obtained after approximately $16$ weeks, corresponding to approximately $2600$ central processing unit (known as CPU) hours.

2.6 Code validation

The present three-dimensional numerical solver has been used and extensively validated in the past for particulate flows (Lashgari et al. Reference Lashgari, Picano, Breugem and Brandt2014), non-Newtonian flows (Rosti & Brandt Reference Rosti and Brandt2017; Alghalibi et al. Reference Alghalibi, Lashgari, Brandt and Hormozi2018; Rosti et al. Reference Rosti, Izbassarov, Tammisola, Hormozi and Brandt2018; Shahmardi et al. Reference Shahmardi, Zade, Ardekani, Poole, Lundell, Rosti and Brandt2019) and multiphase problems in non-Newtonian fluids (De Vita et al. Reference De Vita, Rosti, Izbassarov, Duffo, Tammisola, Hormozi and Brandt2018). In particular, the code has been recently validated for suspensions of rigid and soft particles and droplets in EVP and viscoelastic fluids (Izbassarov et al. Reference Izbassarov, Rosti, Ardekani, Sarabian, Hormozi, Brandt and Tammisola2018). Nonetheless, we report here a further validation case for a simple shear flow of EVP fluid at different power-law indexes $n$.

2.6.1 EVP Couette flow

We consider the two-dimensional shear flow (Couette flow) of an EVP fluid. Initially, the material is at rest and starts flowing due to an externally applied constant shear rate $\dot{\unicode[STIX]{x1D6FE}}_{0}$. The objective here is to track the time evolution of the shear stress ($\unicode[STIX]{x1D70F}_{12}$, with $1$ and $2$ being the streamwise and wall-normal directions, respectively) for different values of the power-law index $n$ and compare them with the analytical solution provided by Saramito (Reference Saramito2009).

The shear Weissenberg number is fixed to $Wi=1$, the Bingham number $Bi=1$ and the retardation parameter $\unicode[STIX]{x1D6FD}=1$. The simulations are performed for the power-law index range $0.05\leqslant n\leqslant 1$. The time evolution of the shear stress $\unicode[STIX]{x1D70F}_{12}$ and the comparison with Saramito (Reference Saramito2009) are depicted in figure 3. We observe that, initially, the shear stress grows linearly, but as soon as the stress level achieves a threshold value, which is called ‘yield stress’, the growth stops and a plateau is reached; this is expected in the yielded state. Indeed as soon as the second invariant of the deviatoric stress tensor exceeds the material yield stress (or $Bi$ number in dimensionless context), the plasticity criteria function $\unicode[STIX]{x1D705}_{n}(|\unicode[STIX]{x1D70F}_{d}|)$ becomes positive and non-zero in the constitutive equation precluding the stress growth. Consequently, the solution remains bounded for any positive value of the power-law index. As shown in the figure, we find a very good agreement with the results by Saramito (Reference Saramito2009).

Figure 3.

Time evolution of shear stress in a stationary simple Couette flow. The solid lines and symbols represent the analytical solution given by Saramito (Reference Saramito2009) and the present numerical simulations, respectively.

3 Results and discussion

In this section, we first present the results of the pure sedimentation simulations, i.e. in the absence of any external shear flow, and then those of the shear-induced sedimentation, i.e. from the simulations involving cross-shear flow.

3.1 Pure sedimentation in EVP fluids

For the case of a sphere settling in an EVP fluid in the absence of cross-shear flow, the dimensionless numbers have been chosen in such a way that the gravity number (related to Bingham number, $Y_{g}$, see (3.1)) remains smaller than its critical value $Y_{g}^{c}$ reported by Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a). Above the critical condition the yield stress resistance overcomes the buoyancy force and, consequently, the particle is trapped inside the yield stress fluid. The critical gravity number is found to be $Y_{g}^{c}=0.143$ for Bingham fluids (i.e. in the absence of elasticity) by Beris et al. (Reference Beris, Tsamopoulos, Armstrong and Brown1985). However, numerical simulations have demonstrated that this value increases with the material elasticity.

Here, we show that the gravity numbers for the current simulations are smaller than the critical gravity number presented as the stoppage criteria by Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a). To this end, we must first rescale our parameters following their definition. Thus, for each simulation we must find the corresponding Deborah number, $De$, and gravity number, $Y_{g}$, defined as

(3.1a,b)

$$\begin{eqnarray}De=\frac{\unicode[STIX]{x1D706}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}gR}{\unicode[STIX]{x1D702}_{p}},\quad Y_{g}=1.5\frac{\unicode[STIX]{x1D70F}_{0}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}gR},\end{eqnarray}$$

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ is the density difference between the bead and the fluid. The critical gravity number is obtained using nonlinear regression on the simulation data in Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a):

(3.2)

$$\begin{eqnarray}\frac{1}{Y_{g}^{c}}=1.2+\frac{1}{0.176+0.135De}:\end{eqnarray}$$

The parameters of the simulations, including the Deborah and gravity numbers, and the comparison with the critical gravity number predicted by (3.2) are shown in table 2. According to this table, the gravity number in the present work is always smaller than its critical value.

Figure 4.

(a) Normalized velocity magnitude around the particle settling in the $y$–$z$ centreplane ($x=3D$) through Bingham fluid and (b) EVP material at $Wi_{\infty }=0.1$. (c) Streamwise velocity around a particle settling in a Bingham fluid and an EVP material at the channel centre line ($x=3D$, $z=6D$). In all figures, the Bingham number $Bi=0.5$.

Table 2.

Dimensionless numbers of the pure sedimentation simulations with associated gravity number $Y_{g}$, smaller than its critical value $Y_{g}^{c}$.

3.1.1 Velocity and stress fields

In this section the velocity and stress fields around a sphere are demonstrated to show the effects of elasticity and plasticity on the flow dynamics. To better highlight the elasticity effects on the flow features and the sphere drag, we have performed the additional simulations of particle settling in a Bingham fluid, whose details are presented in appendix A.

Figure 4 shows the contour plot of the velocity magnitude in the mid-$y$–$z$ plane around the surface of the particle settling in an otherwise quiescent EVP fluid. An analogous map for the Bingham fluid is included for the sake of comparison (see the methods and details in appendix A). The $Bi$ number is held fixed and equal to $0.5$ for both cases and the settling Weissenberg number $Wi_{\infty }=0.1$ for the EVP material and $0$ for the Bingham fluid. The magnitude of the velocity is normalized by the particle settling velocity, $u_{sett}$. Note that we depict the flow only in a small window around the particle to highlight the shape of the velocity contours.

We observe that the symmetry of the velocity field around the particle stagnation points (north and south poles of the sphere) is broken and an overshoot in the downstream velocity appears as soon as a small elasticity ($Wi_{\infty }=0.1$) is added to the carrying fluid. This behaviour is quantitatively depicted in figure 4(c), where the normalized streamwise velocity is plotted in the flow direction at the centre line of the channel, i.e. as function of $y$ for $x=3D$, $z=6D$.

The loss of the fore–aft symmetry in a yield stress fluid and the negative wake formation have been previously observed in experiments (Gueslin et al. Reference Gueslin, Talini, Herzhaft, Peysson and Allain2006; Putz et al. Reference Putz, Burghelea, Frigaard and Martinez2008; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Ahonguio et al. Reference Ahonguio, Jossic and Magnin2014). Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a) revealed that the elasticity of realistic yield stress fluids is the primary cause of such phenomena. Indeed, these authors demonstrated that the thixotropy (aging of yield stress materials) is not responsible for such behaviour as conjectured before. The absence of this symmetry was also observed in the case of a neutrally buoyant sphere in Carbopol gels subjected to simple shear flow (Firouznia et al. Reference Firouznia, Metzger, Ovarlez and Hormozi2018). Therefore, to correctly predict the behaviour of practical yield stress fluids, the effect of elasticity should be taken into account. It is also worthwhile to mention that the loss of the fore–aft symmetry is not due to the weak inertia of our simulations ($Re_{p}=1$), see § 3.1.3 for an explanation.

The velocity vector field around the particle settling in an EVP fluid is illustrated in figure 5 for the same case as in figure 4.

Figure 5.

Velocity vectors around the sphere settling in an EVP material in the central $y$–$z$ plane. The dashed red and solid green boxes magnify the recirculation zones and the flow stagnation points. The dimensionless numbers are the same as figure 4.

For clarity, the recirculation zones (inside of the dashed red box) and flow stagnation points (from the solid green box) are magnified. The negative wake is observed where the fluid velocity is opposite the particle velocity, upstream of the rear stagnation point. The overshoot in the stream velocity downstream of the particle shown in figure 4(c) occurs in the same area. For a particle with no-slip boundary condition, the recirculation zones arise close to the particle and in its equatorial plane as previously observed experimentally (Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012), which could also be captured computationally by adopting an extensive mesh refinement near the particle surface (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a).

The interaction between shear and normal stresses upstream of the particle is shown to be the main cause for the negative wake formation (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a). In fact, the negative wake is due to the normal stress relaxing faster than the shear stress away from the particle. This is shown in the colourmap of normal and shear stresses around the sphere in figure 6(a,b). The stress relaxation away from the sphere, in particular the shear and normal stresses, is also depicted along the ray at an angle $\unicode[STIX]{x1D701}=30^{\circ }$, measured from the rear stagnation point (sphere north pole) in figure 6(c) (we define the north pole in the positive $y$ direction).

Figure 6.

Stress fields around the particle settling in an EVP fluid in the $y$–$z$ plane at $x=3D$; (a) the $\unicode[STIX]{x1D70F}_{yz}$ shear stress and (b) the $\unicode[STIX]{x1D70F}_{yy}$ normal stress. The angle $\unicode[STIX]{x1D701}$ is measured from the north pole of the sphere. (c) Stress relaxation versus the distance from the particle, $r$, at an angle $\unicode[STIX]{x1D701}=30^{\circ }$ measured from the north pole of the sphere. Here, $r=0$ corresponds to the particle surface. The absolute value of the stresses are shown for the same case as in figure 4.

Clearly, at $r\approx 0.2D$, the absolute value of the shear stress is twice as large as the normal stress. Thus, the elastic shear stress, which is responsible for the creation of a rotational force on the fluid element, contributes to the negative wake and to the secondary flow upstream of the particle. For more details the reader is referred to Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a).

3.1.2 Yielded/unyielded regions

In this section, we show the shape and extent of the yielded/unyielded zones around a sphere at different $Bi$ numbers in an EVP material. To our knowledge, the unyielded regions have not been previously reported in three dimensions. Figure 7 displays the surfaces delimiting the yielded regions at different $Bi$ numbers and constant elasticity ($Wi_{\infty }=0.1$), where unyielded surfaces are plotted with lower transparency for higher $Bi$ numbers ($Bi=0.5,1$) for the sake of clarity. In the figure, red represents the unyielded regions, while grey depicts the yielded surface boundary and the particle. For completeness, figure 8 shows the projection of unyielded surfaces onto the mid-$y$–$z$ and $x$–$y$ planes. In this figure, blue and grey colours denote the yielded and unyielded zones, respectively.

Figure 7.

Surfaces of unyielded regions around the particle in an EVP material at $Wi_{\infty }=0.1$ in the absence of cross-shear flow for various $Bi$: (a) $Bi=0.05$, (b) $Bi=0.1$, (c) $Bi=0.13$, (d) $Bi=0.5$, (e) $Bi=1$. Red represents the unyielded zone, while grey shows the yield surface boundary.

Figure 8.

Evolution of the unyielded zones for flow of an EVP material around a particle for various $Bi$ numbers at $Wi_{\infty }=0.1$; (a,f) $Bi=0.05$, (b,g) $Bi=0.1$, (c,h) $Bi=0.13$, (d,i) $Bi=0.5$, (e,j) $Bi=1$. The first (a–e) and second rows (f–j) represent the unyielded zones in the central $y$–$z$ and $x$–$y$ planes, respectively. Blue and grey represent the yielded and unyielded regions.

Figure 9.

Colourmaps of the second invariant of the shear rate for the flow of an EVP material around a particle for various $Bi$ numbers at $Wi_{\infty }=0.1$; (a,f) $Bi=0.05$,(b,g) $Bi=0.1$, (c,h) $Bi=0.13$, (d,i) $Bi=0.5$, (e,j) $Bi=1$. The first (a–e) and second rows (f–j) represent the shear rate in the central $y$–$z$ and $x$–$y$ planes, respectively.

The simulations reveal the existence of two unyielded zones: the first is the unyielded envelope that surrounds the fluid zone and the second is the unyielded ring located in the yielded zone surrounding the sphere. The projection of this ring in the central $y$–$z$ and $x$–$y$ planes is shown in figure 8. These are the two unyielded islands in the yielded region located at the equator and on either side of the sphere. It is noteworthy to mention that these solid islands are the regions in which the second invariant of the shear rate is almost zero as shown in figure 9. Therefore, unlike the case of 2-D cylinders, these solid rings are not rotating solid islands; note also that, at steady state, the sphere-rotating velocity is zero in the pure sedimentation cases.

We note that the outer unyielded envelope grows progressively as the $Bi$ number increases and the yield surface boundary approaches the surface of the particle from the equator plane causing the particle to stop settling. A similar arrest mechanism has been captured previously in axisymmetric particle-settling simulations in an EVP material (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a). Moreover, figure 8 shows that there exists yielded regions in the vicinity of the channel walls for $Bi=0.1$ and $0.13$. These are associated with the wall effects and had also previously been observed by Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997) in the Bingham fluid flow past a sphere contained in a tube with a diameter that is 10 times larger than that of the sphere.

3.1.3 Drag coefficients

In our simulations we fix the particle and compute the drag exerted on the particle by the surrounding fluid. Since the settling rate is inversely proportional to the drag, drag enhancement is equivalent to settling rate reduction (Padhy et al. Reference Padhy, Shaqfeh, Iaccarino, Morris and Tonmukayakul2013b). The drag coefficient $C_{d}$ is defined as follows:

(3.3)

$$\begin{eqnarray}C_{d}=\frac{2F_{d}}{\unicode[STIX]{x1D702}_{0}U_{\infty }D},\end{eqnarray}$$

where $F_{d}$ is the total drag force exerted on the particle in the streamwise $(y)$ direction.

To investigate the dependency of the drag change on the flow and fluid parameters, we decompose the total drag force on the particle into its individual components, which are associated with the different stress contributions. As we assume a finite particle Reynolds number ($Re_{p}=1$), the total drag consists of four components: ‘form drag’ ($F_{d}^{f}$), ‘viscous drag’ ($F_{d}^{v}$), ‘polymer drag’ ($F_{d}^{p}$) and ‘inertia drag’ ($F_{d}^{I}$). The different components are computed from the following definitions:

(3.4)

$$\begin{eqnarray}\displaystyle & \displaystyle F_{d}^{f}=-\iint _{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}pn_{y}\,\text{d}S, & \displaystyle\end{eqnarray}$$

(3.5)

$$\begin{eqnarray}\displaystyle & \displaystyle F_{d}^{v}=(1-\unicode[STIX]{x1D6FD})\iint _{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}\left(\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}y}\right)n_{j}\,\text{d}S, & \displaystyle\end{eqnarray}$$

(3.6)

$$\begin{eqnarray}\displaystyle & \displaystyle F_{d}^{p}=\iint _{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D70F}_{yj}n_{j}\,\text{d}S, & \displaystyle\end{eqnarray}$$

(3.7)

$$\begin{eqnarray}\displaystyle & \displaystyle F_{d}^{I}=Re_{p}\iint _{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}u_{j}\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}x_{j}}n_{j}\,\text{d}S, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ is non-dimensionalized by $D^{2}$. The details of the numerical integration procedure are given in appendix B.

The form drag (3.4) is the drag component resulting from the distribution of the dynamic pressure on the sphere. It is noteworthy to mention that in the Saramito’s model (Saramito Reference Saramito2009) the extra elastic stress tensor is not traceless, i.e. $\text{tr}(\unicode[STIX]{x1D70F})$ is not zero. Thus, the pressure field $p$ obtained in the numerical solution is basically the field of $\bar{p}=p+(1/3\text{tr}(\unicode[STIX]{x1D70F}))$. However, in all of our simulations, the absolute value of the $1/3\text{tr}\unicode[STIX]{x1D70F}$ at the surface of the sphere is negligible when compared to the absolute value of $\bar{p}=p+(1/3\text{tr}\unicode[STIX]{x1D70F})$. In addition, unlike $\bar{p}$, the sign of $1/3\text{tr}\unicode[STIX]{x1D70F}$ does not change at the stagnation points. Therefore, we conclude that the contribution of $1/3\text{tr}\unicode[STIX]{x1D70F}$ to the drag is negligible as compared to the contribution of the dynamic pressure field $p$ for the simulations conducted in this work.

To significantly reduce the time of each computation, we have assumed small but finite inertia, $Re_{p}=1$, and still expect inertial effects to be negligible. To a posteriori check this, we have compared the sum of the form, viscous and polymer drag force to the IB force $\boldsymbol{f}_{\boldsymbol{j}}$, i.e. the total drag, and found a difference of approximately $1.5\,\%$, which indeed indicates negligible inertial drag, $F_{d}^{I}$.

Figure 10.

Total drag coefficient and its individual components for a sphere settling in an EVP fluid as a function of the $Bi$ number at $Wi_{\infty }=0.1$.

Figure 10 illustrates the variation of the total drag coefficient and its individual contributions (form, viscous and polymer drag) experienced by the sphere settling in an EVP fluid as a function of the $Bi$ number. The total drag increases as the material yield stress is increased, as observed in the past for the particle settling in a purely viscoplastic fluid (Beris et al. Reference Beris, Tsamopoulos, Armstrong and Brown1985; Atapattu et al. Reference Atapattu, Chhabra and Uhlherr1995; Blackery & Mitsoulis Reference Blackery and Mitsoulis1997; Tabuteau et al. Reference Tabuteau, Coussot and de Bruyn2007; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Ahonguio et al. Reference Ahonguio, Jossic and Magnin2014; Wachs & Frigaard Reference Wachs and Frigaard2016) as well as in an EVP fluid (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a). Here, the key finding of our 3-D simulations is that the drag contribution from the dynamic pressure distribution on the particle (form drag) is dominant. Additionally, the viscous drag increases with $Bi$, whereas the drag contribution due to the polymer stresses remains almost constant with $Bi$.

The dynamic pressure and the three components of the viscous stress distribution ($\unicode[STIX]{x1D70F}_{xy},\unicode[STIX]{x1D70F}_{yy},\unicode[STIX]{x1D70F}_{yz}$) are computed on the surface of the sphere and are shown in figure 11 (the computational procedure is explained in detail in appendix B). Note that in this study we are dealing with a complex flow, i.e. the type of flow is in general different at each point, e.g. predominantly shear, extensional, nearly rigid-body motion or a combination of these. However, our goal is to show the components of the stress that contribute to the drag force, i.e. $F_{i}=\int _{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D70F}_{ij}\cdot n_{j}\,\text{d}A$. For a drag force in the $y$ direction, the relevant components of the stress tensor are $\unicode[STIX]{x1D70F}_{yx},\unicode[STIX]{x1D70F}_{yy}$ and $\unicode[STIX]{x1D70F}_{yz}$, so we just examine these components of the viscous stress tensor, since figure 10 shows the polymer stress contribution to the drag force is negligible.

Table 3.

Comparison the total drag coefficient and its components on the surface of the sphere settling in a viscoplastic and in an EVP fluid at $Bi=0.05$. The settling Weissenberg number $Wi_{\infty }=0.1$ for EVP material.

Figure 11.

Contour plots of the dynamic pressure and the three components of the viscous stress distribution on the surface of the particle settling in an EVP fluid at different $Bi$ numbers and $Wi_{\infty }=0.1$. First, second, third, fourth and fifth rows show the distributions at $Bi=0.05,0.1,0.13,0.5$ and $1$, respectively. Axes information in all panels can be extrapolated from panel (a).

The visualizations in figure 11 confirm that the magnitude of the pressure on the particle surface increases with $Bi$. Moreover, the fore–aft symmetry in the pressure distribution around the particle stagnation points breaks leading to further drag enhancement. Increasing the Bingham number, on the other hand, merely magnifies the value of the shear and normal viscous stresses around the poles and on either side of the sphere and does not break the fore–aft symmetry. Hence, in the presence of elasticity, the sphere drag increases through modification of the dynamic pressure distribution by both breaking its fore–aft symmetry with respect to its north and south poles and by magnifying its magnitude.

Table 3 compares the total drag coefficient and its individual components for a sphere settling in an ideal viscoplastic (Bingham) fluid and an EVP material. Three features are evident here. First, adding a small amount of elasticity to an ideal-yield stress fluid causes a total drag reduction of approximately $20\,\%$, in agreement with the observations by Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a). Second, adding elasticity to the purely viscoplastic fluid (Bingham fluid) modifies the stress fields on the sphere surface. Table 3 indeed shows that the dominant drag component for the Bingham fluid case is the viscous drag, while it is the form drag in an EVP fluid. Third, although the polymer drag contribution is small compared to the other components, the total drag is indirectly affected by the presence of the polymers through modification of the form and viscous stresses.

Given the total drag, we compute the elastic contribution (elastic drag) for a sphere settling in a yield stress fluid using the following relation:

(3.8)

$$\begin{eqnarray}(\unicode[STIX]{x1D702}_{0}U_{\infty }D)C_{d}^{EVP}=(\unicode[STIX]{x1D702}_{0}U_{\infty }D)C_{d}^{VP}-(\unicode[STIX]{x1D702}_{0}U_{\infty }D)C_{d}^{e},\end{eqnarray}$$

where $\unicode[STIX]{x1D702}_{0}U_{\infty }D$ is the viscous force scale, and $C_{d}^{EVP}$ and $C_{d}^{VP}$ the total drag coefficient on the sphere settling in an EVP material and a viscoplastic fluid, respectively. Thus, $(\unicode[STIX]{x1D702}_{0}U_{\infty }D)C_{d}^{e}$ quantifies the indirect effect of adding polymers to a yield stress fluid by modifying the dynamic pressure and viscous stresses.

Note that throughout this manuscript we use the formulation ‘adding polymers to a Bingham fluid’ to indicate adding a finite elasticity to a Bingham fluid, although the elastic effects can be present in practical yield stress fluids also in the absence of polymers, e.g. in emulsions. Moreover, in this manuscript we use the term ‘polymer drag’ to denote the drag contribution due to the extra elastic stress tensor and the term ‘elastic drag’ as introduced by (3.8) to indicate the total change in the drag force (through modifying the pressure and viscous stresses as well as the addition of the polymer stresses) due to the addition of a finite elasticity to the ideal yield stress fluids.

The variation of the total drag for a sphere settling in a Bingham fluid as a function of $Bi$ is then found using the correlation given by Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997), which is described in § A.2, assuming the coefficients in (A 4) $a=2.343$ and $b=0.879$ corresponding to a confinement ratio of $12$ (Blackery & Mitsoulis Reference Blackery and Mitsoulis1997). The total drag coefficient for settling in Bingham and EVP fluids along with the elastic drag $C_{d}^{e}$ defined above are displayed in figure 12 versus the Bingham number. Clearly, the elastic drag is an increasing function of $Bi$ for a particle settling in EVP material at constant $Wi_{\infty }$. Moreover, the drag reduction through introduction of elasticity, which was shown in table 3 for a single $Bi$ number ($Bi=0.05$) applies to a range of $Bi$ numbers ($0.05\leqslant Bi\leqslant 1$).

We can relate the sphere drag reduction in the EVP material to the volume of the yielded region around the sphere; indeed, increasing the Deborah number at a fixed Bingham number results in an increase of the yielded region around the sphere, as previously shown by Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a). Consequently, elasticity helps the sphere to translate faster in the EVP material by reducing the drag at high elasticity.

Figure 12.

Total drag coefficient of a particle settling in a Bingham fluid from Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997) and that of a particle settling in an EVP material (with elasticity $Wi_{\infty }=0.1$) as well as the elastic drag resulting from (3.8).

Figure 13.

Contour plots of the dynamic pressure and of three components of viscous stress distribution on the surface of the particle settling in a Bingham fluid (a–d) at$Bi=0.05$ and in the EVP fluid (e–h) at $Bi=0.05$, $Wi_{\infty }=0.1$. Axes information in all panels can be extrapolated from panel (a).

To gain further insight, we display in figure 13 the contour plots of the dynamic pressure and of the $\unicode[STIX]{x1D70F}_{xy}$, $\unicode[STIX]{x1D70F}_{yy}$ and $\unicode[STIX]{x1D70F}_{yz}$ components of the viscous stress distributions on the sphere surface for the case of a Bingham fluid (a–d) and those of an EVP fluid (e–h). The pressure distribution and its magnitude are not affected by the presence of the polymers as seen in table 3 where the form drags are found to be close to each other. Nevertheless, the magnitudes of the viscous shear and normal stresses drop by adding polymers to the pure Bingham fluid.

Figure 14.

Pressure and viscous stress distributions on the surface of the sphere settling in Bingham and EVP fluids for the same cases depicted in figure 13. (a) Spherical coordinate system, (b) pressure, (c) $\unicode[STIX]{x1D70F}_{xy}$, (d) $\unicode[STIX]{x1D70F}_{yy}$, (e) $\unicode[STIX]{x1D70F}_{yz}$. The pressure, $\unicode[STIX]{x1D70F}_{xy}$ and $\unicode[STIX]{x1D70F}_{yy}$, distributions are plotted for $\unicode[STIX]{x1D713}=0$, while the $\unicode[STIX]{x1D70F}_{yz}$ is plotted at $\unicode[STIX]{x1D719}=0$.

This behaviour can be better understood by looking at figure 14, where the pressure and the $\unicode[STIX]{x1D70F}_{xy}$ and $\unicode[STIX]{x1D70F}_{yy}$ components of the viscous stresses are plotted around the sphere stagnation points (at zero polar angle $\unicode[STIX]{x1D713}=0$ and as a function of azimuthal angle $\unicode[STIX]{x1D719}$, see coordinate system in panel a). Figure 14 also shows the $\unicode[STIX]{x1D70F}_{yz}$ component of the viscous stresses on either side of the sphere at $\unicode[STIX]{x1D719}=0$ and as a function of polar angle $\unicode[STIX]{x1D713}$. Polymers affect the viscous stresses on the sphere surface, which results in total drag reduction compared to the case of a Bingham fluid. Furthermore, the slope of the $\unicode[STIX]{x1D70F}_{yz}$ contribution of the viscous stress changes in the case of an EVP fluid.

3.2 Shear-induced sedimentation

In this section, we present the results of the simulations of shear-induced sedimentation of a sphere in an EVP material. In our set-up, the only non-zero component of the particle rotational velocity is $\unicode[STIX]{x1D714}_{y}$, owing to symmetry. The simulations are performed at a constant ratio of shear to sedimentation Weissenberg number $\unicode[STIX]{x1D6FC}=Wi/Wi_{\infty }=0.1$. Equivalently, the ratio between the externally imposed cross-shear rate to the shear rate induced by settling is kept fixed and equal to $0.1$ ($\dot{\unicode[STIX]{x1D6FE}}_{0}/\dot{\unicode[STIX]{x1D6FE}}_{sett}=0.1$). This means that the cross-shear flow is always a secondary flow and the uniform flow is the primary flow.

3.2.1 Velocity field

The velocity field around a sphere settling in a sheared and quiescent EVP fluid is presented in this section. Figure 15 illustrates the velocity magnitude normalized by the particle settling velocity in the mid-plane between the two walls; results are presented for settling both in the presence and absence of an externally imposed cross-shear flow.

Figure 15.

(a) Normalized velocity colourmaps in the $y$–$z$ centreplane ($x=3D$) for the case of pure sedimentation of a single sphere in an EVP material. (b) The same as panel (a), but for the case of shear-induced sedimentation. (c) Streamwise fluid velocity in the flow direction at the channel centreline ($x=3D$, $z=6D$) for the cases of pure sedimentation (orange markers) and shear-induced sedimentation (blue markers) in an EVP material. In all figures the dimensionless numbers are $Bi=0.5$, $Wi_{\infty }=0.1$.

We observe that the fore–aft asymmetry of the velocity field around the sphere is less pronounced when the shear flow is superimposed on the uniform flow. This feature is further clarified in figure 15(c), which illustrates the differences due to the presence of the cross-shear flow in terms of streamwise velocity. Clearly, the velocity overshoot downstream of the sphere disappears once the secondary orthogonal cross-shear flow is superimposed to the primary uniform flow.

Furthermore, the negative wake is eliminated by imposing the cross-shear flow, as demonstrated in figure 16 by the velocity vector field around the settling sphere.

Figure 16.

Velocity vectors around the sphere settling through an EVP material ($Bi=0.5$ and $Wi_{\infty }=0.1$) in the centre $y$–$z$ plane for the case of (a) pure sedimentation, (b) shear-induced sedimentation. The blue boxes magnify the flow upstream of the sphere.

3.2.2 Yielded/unyielded regions

It should be noted that all the regions that remain unyielded in the case of pure sedimentation in figures 7 and 8 yield as soon as the fluid is sheared. Further, in all of the shear-induced sedimentation simulations, the material (either pure viscoplastic or EVP) is yielded when it enters the computational domain. In other words, the components of the polymer stress tensor satisfy the von Mises yielding criterion at the inlet. These components are obtained by solving (2.1)–(2.3) at steady state for the uniform and Couette flow of Bingham and EVP fluid analytically (see § 2.3 and appendix A for the analytical solution of the stress tensor components in the case of EVP and Bingham fluid, respectively).

3.2.3 Drag coefficients

The objective here is to study how the sphere drag coefficient is influenced by the nonlinear coupling of the primary uniform flow and the secondary linear cross-shear flow. The drag force is computed as the volumetric sum of the IB forces from the numerical implementation and, from this, the drag coefficient using (3.3).

Figure 17.

Drag coefficient for a settling sphere in an EVP fluid ($Wi_{\infty }=0.1$, $Wi=0.01$) for the case of pure sedimentation (blue circular markers) and shear-induced sedimentation (red triangular markers) as a function of the Bingham number. The green square markers represent the drag coefficient deduced from the pure sedimentation simulations data by taking into account the second invariant of the shear rate tensor.

Figure 17 shows the drag coefficients for the cases of pure sedimentation (circular markers) and shear-induced sedimentation (triangular markers) of a single sphere in an EVP material. The total drag is an increasing function of the Bingham number, which is similar to the general trend shown in figure 10 for the case of pure sedimentation. Clearly, superimposing a cross-shear flow on the uniform flow results in the drag reduction on the sphere settling in an EVP fluid (see red triangular markers in figure 17). This is a result of the fact that shearing the yield stress fluid in the orthogonal direction with respect to the uniform flow caused the whole medium to become fully yielded. It is intuitive that the overall viscosity of a yield stress fluid reduces as we impose a cross-shear flow and the settling particle experiences a lower viscous resistance. In addition, a sphere settling in a yield stress fluid without any cross-shear flow experiences a larger confinement compared to settling in the presence of external shear flow because the unyielded zones in the EVP fluid act as elastic walls and it is well known that increasing the confinement ratio results in drag enhancement in Newtonian (see Faxén Reference Faxén1922) and viscoelastic fluids (see e.g. Lunsmann et al. Reference Lunsmann, Genieser, Armstrong and Brown1993; Harlen Reference Harlen2002).

The total and individual drag contributions for the case of shear-induced sedimentation, computed using (3.4)–(3.6), are compared with the results pertaining pure sedimentation in figure 18.

Figure 18.

Total and individual drag components on the sphere settling in an EVP fluid at $Wi_{\infty }=0.1$ in the absence of cross-shear flow (dark markers) and in the presence of cross-shear flow (light markers) as a function of $Bi$ number.

According to the data in the figure, the form drag is also the dominant component in the case of shear-induced sedimentation, as documented for a particle settling in the absence of an externally imposed shear flow in figure 10. The viscous drag is the second largest contribution whereas the polymer drag provides the smallest contribution to the total drag in both cases. The form and viscous drag components are both an increasing function of $Bi$ at constant sedimentation Weissenberg number, $Wi_{\infty }$, in which case the polymer drag remains nearly constant as $Bi$ is increased.

In summary, in pure sedimentation flows, adding small amounts of elasticity to the ideal yield stress fluid indirectly modifies the dynamic pressure and viscous stress distributions on the surface of the sphere such that form and viscous drag contributions are responsible for the observed total drag enhancement. Moreover, by superimposing a relatively weak cross-shear flow to the primary settling flow, the form and viscous drag components drop, while the polymer drag remains almost the same.

To further clarify the source of drag reduction resulting from the external cross-shear flow, we examine in detail the case $Bi=1$. The distribution of dynamic pressure as well as the viscous stress components that cause drag are therefore displayed in figure 19 for pure sedimentation and shear-induced sedimentation.

Figure 19.

Contour plots of the pressure and viscous stress components on the surface of the sphere settling in an EVP fluid at $Bi=1$, $Wi_{\infty }=0.1$. The first row corresponds to pure sedimentation ($Wi=0$) and the second row to shear-induced sedimentation ($Wi=0.01$). Axes information in all panels can be extrapolated from panel (a).

This figure shows that the magnitude of the pressure and viscous shear stresses decrease as soon as the cross-shear flow is superimposed on the uniform flow, whereas, the viscous normal stress, $\unicode[STIX]{x1D70F}_{yy}$, remains almost unaffected.

3.3 Shear-induced sedimentation in a Bingham fluid

In this section, we study the effect of superimposing the simple cross-shear flow on the orthogonal uniform flow of a Bingham fluid over a single sphere. As the material is sheared, the whole medium becomes yielded and there is no rigid zone left in the medium.

Table 4 shows the total drag coefficient and its components (form and viscous drag) on the sphere settling in a Bingham fluid in the absence and the presence of an externally imposed cross-shear flow.

Table 4.

Comparison of the drag coefficient and its individual components on a sphere settling in the absence and the presence of cross-shear flow in a Bingham fluid at $Bi=0.05$.

In both cases, the viscous drag component is the dominant one. Furthermore, the presence of the cross-shear flow results in the reduction of both viscous and form drag. This behaviour can be inferred from the pressure and viscous stress distributions in figure 20.

Figure 20.

Contour plots of the dynamic pressure and of the three components of the viscous stress on the surface of the particle in the case of pure sedimentation (a–d) and shear-induced sedimentation (e–h) in a Bingham fluid at $Bi=0.05$. Axes information in all panels can be extrapolated from panel (a).

We next compare the drag on a sphere settling in a Bingham fluid to that of a sphere settling in an EVP fluid where both fluids are subjected to cross-shear flow at the same $Bi$ number. This would be useful to study the effect of elasticity on the sphere drag settling in sheared yield stress fluid as the elasticity is absent for the Bingham fluid. Table 5 provides this comparison.

Table 5.

Comparison of the drag coefficient and its individual components on a sphere settling in sheared viscoplastic (Bingham) fluid at $Bi=0.05$ and EVP material at $Wi_{\infty }=0.1$, $Bi=0.05$.

By introducing the elasticity, the total drag coefficient is reduced. The drag reduction is mainly due to a decrease in viscous drag, while the form drag is hardly affected by the elasticity. Adding polymers to a Bingham fluid indirectly decreases the sphere drag by reducing the viscous stresses on the sphere surface.

3.4 Nonlinear coupling

The Bingham number can be defined in two ways depending on how one selects the scale of the shear rate that can be either the one induced by the settling of the sphere $(\dot{\unicode[STIX]{x1D6FE}}_{sett})$ or by the scalar shear rate that includes both the shear and the settling flow. For the latter choice, the following second invariant of the deformation rate tensor $(\dot{\unicode[STIX]{x1D6FE}}_{inv})$ can be defined:

(3.9)

$$\begin{eqnarray}\dot{\unicode[STIX]{x1D6FE}}_{inv}=\sqrt{\dot{\unicode[STIX]{x1D6FE}}_{sett}^{2}+\dot{\unicode[STIX]{x1D6FE}}_{0}^{2}}=\frac{U_{\infty }}{D}\sqrt{1+\unicode[STIX]{x1D6FC}^{2}}.\end{eqnarray}$$

This is due to the fact that the cross-shear flow results in an additional off-diagonal element of the strain rate tensor. Therefore, $Bi$ and $Bi_{inv}$ can be defined by the following relations:

(3.10a,b)

$$\begin{eqnarray}Bi=\frac{\unicode[STIX]{x1D70F}_{0}}{\unicode[STIX]{x1D702}_{0}\dot{\unicode[STIX]{x1D6FE}}_{sett}},\quad Bi_{inv}=\frac{\unicode[STIX]{x1D70F}_{0}}{\unicode[STIX]{x1D702}_{0}\dot{\unicode[STIX]{x1D6FE}}_{inv}}.\end{eqnarray}$$

In the shear-induced sedimentation simulations, the shear flow is secondary to the sedimentation flow, and one might try to deduce the drag coefficient from the pure sedimentation simulations data. This means that, having the relationship between $C_{d}$ and $Bi$ for settling flows, one can estimate the drag coefficient for shear-induced sedimentation flows by replacing $Bi$ with $Bi_{inv}$, with the green square markers in figure 17 representing the result of such an estimation. Because in our simulations the shear flow is secondary, and consequently $Bi\sim Bi_{inv}$, the mentioned estimate results in a drag coefficient that is almost the same for the settling and shear induced sedimentation flows (as in figure 17 the blue circles are very close to the green squares). However, our simulations show that the drag coefficients for the shear-induced sedimentation (red triangles in figure 17) are much smaller than the values of $C_{d}$ in the pure sedimentation simulations (blue circles in figure 17). This implies that the coupling of the two orthogonal flows play a significant role in determining the sphere drag and the mentioned estimate is not valid even when one of the flows is one order of magnitude smaller than the other.

To expand on the above discussion, our simulations show that the changes in pressure and viscous stresses are not linear with respect to changes in the second invariant of shear rate. For instance at $Bi=0.5$ for a settling flow, a relative increase of $0.5\,\%$ in $\dot{\unicode[STIX]{x1D6FE}}_{inv}$ by superimposing a secondary shear flow contributes to approximately a $28\,\%$ relative decrease in the dynamic pressure, and $18\,\%$ in all of the viscous stress components. This is due to the nonlinear coupling of the secondary simple cross-shear flow to the primary uniform flow of the EVP fluid.

3.5 Shear-induced sedimentation at higher elasticity

In this section, we show the results of the simulations of a single spherical particle settling in a sheared EVP fluid at higher elasticity, i.e. at $Wi_{\infty }=1$. These simulations were performed for the same range of $Bi$ number, i.e. $0.05\leqslant Bi\leqslant 1$. Note again that the ratio between shear and settling Weissenberg number ($\unicode[STIX]{x1D6FC}=Wi/Wi_{\infty }$) is kept fixed to $\unicode[STIX]{x1D6FC}=0.1$. Therefore, the shear Weissenberg number $Wi=0.1$.

Figure 21.

Drag coefficient versus the Bingham number for a sphere settling in sheared EVP fluid at different settling Weissenberg numbers, $Wi_{\infty }$.

The drag coefficient for a sphere settling in an EVP fluid subject to an externally imposed shear flow at these two different settling Weissenberg numbers is reported in figure 21 as a function of $Bi$. For both values of $Wi_{\infty }$ considered, the drag is an increasing function of $Bi$, while it is a decreasing function of $Wi_{\infty }$ for all $Bi$ examined. In other words, the drag reduces as the material elasticity increases at constant plasticity.

Drag reduction for a particle settling in an EVP material at higher elasticity was previously observed in the computations of Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a) of pure sedimentation. In this case, the particle translates in a yielded envelope; as the material elasticity is increased, the volume of this envelope increases while the rigid polar caps shrink because the von Mises yielding criterion is satisfied more easily due to the larger elastic stresses in the medium. Consequently, the elastic walls move further away from the particle surface at larger values of $Wi_{\infty }$ and the particle experiences lower confinement at higher elasticity, which results in drag reduction (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a).

As far as the mechanism of the drag reduction due to the elasticity in a shear-induced sedimentation case is concerned, we display the total drag as well as its components for $Wi_{\infty }=0.1,1$ at $Bi=1$ in table 6. As the material elasticity is increased, all of the drag components decrease. Nevertheless, our simulations show that the form and viscous drag reduction is more pronounced and they are the dominant drag components.

Figure 22.

Contour plots of pressure and viscous stress components projected on the surface of the sphere settling in sheared EVP fluid at $Bi=1$. The first row corresponds to $Wi_{\infty }=0.1$ ($Wi=0.01$) and the second row to $Wi_{\infty }=1$ ($Wi=0.1$). Axes information in all panels can be extrapolated from panel (a).

Table 6.

Comparison of the drag coefficient and its individual components on a sphere settling in a sheared EVP fluid at $Bi=1$ for two values of $Wi_{\infty }$.

To gain further insight, the pressure and viscous stress distributions on the sphere surface are displayed in figure 22 for the same two cases above, i.e. for $Bi=1$ and $Wi_{\infty }=0.1,1$. It can be seen that the asymmetry in the positive and negative pressure distributions around the sphere stagnation points at zero polar angle is enhanced at higher elasticity. Moreover, the absolute value of the pressure on the sphere surface is reduced causing form drag reduction at higher elasticity. The normal viscous stress is relaxed at higher shear Weissenberg number. The loss of symmetry of the positive and negative viscous shear stresses, $\unicode[STIX]{x1D70F}_{xy}$ and $\unicode[STIX]{x1D70F}_{yz}$, is more noticeable around the north and south poles of the sphere, at zero polar angle, and in the equatorial plane on either side of the sphere at zero azimuthal angle. Note that the breaking of the fore–aft symmetry of the velocity magnitude around the sphere by elasticity in a yield stress fluid was observed previously experimentally (see e.g. Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012; Firouznia et al. Reference Firouznia, Metzger, Ovarlez and Hormozi2018) and computationally (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a).

4 Conclusion

Direct numerical simulations are performed to study the sedimentation of a single sphere confined in a quiescent and a sheared yield stress fluid at small particle Reynolds number ($Re_{p}=1$). The 3-D simulations are performed for both ideal yield stress fluids (using Bingham model) and EVP materials as the carrying fluid. The former fluids exhibit viscous and plastic behaviours while EVP materials also display elastic effects. Consequently, we investigate the impact of elasticity on the flow behaviour and dynamics of particle settling. In all of the simulations the ratio of the constant externally imposed shear rate ($\dot{\unicode[STIX]{x1D6FE}_{0}}$) to the shear rate induced by the particle settling ($\dot{\unicode[STIX]{x1D6FE}}_{sett}$), defined as the ratio of particle settling velocity to its diameter $(U_{\infty }/D)$ is fixed and equal to $0.1$. For the simulations of EVP materials, the ratio of the shear to sedimentation Weissenberg number $(Wi/Wi_{\infty })$ is held constant and equal to $0.1$.

The constitutive equations from Saramito (Reference Saramito2009) are implemented to model the EVP material. The non-Newtonian stress tensor is fully coupled with the flow equations that are solved with a fast and highly scalable finite-volume method with fast Fourier transform-based pressure solver. The immersed boundary method, with a computationally efficient multidirect forcing scheme, is adopted to represent the rigid spherical particle. The no-slip/no-penetration boundary condition on the surface of the particle are therefore implicitly imposed by adding a virtual body force to the right-hand side of the momentum equations.

As concerns the flow and particle dynamics, the fluid velocity distribution around the particle is symmetric with respect to the sphere equatorial plane for pure sedimentation in a Bingham fluid. However, the fore–aft symmetry breaks and the formation of a negative wake observed as the sphere settles in an otherwise quiescent EVP material. Since the thixotropy of the yield stress fluid is not considered in the Saramito (Reference Saramito2009) constitutive equation for EVP fluids, elasticity is the primary cause of both the fore–aft asymmetry in the velocity field and the negative wake. This is in line with the recent computations of Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a) and several previous experiments (see e.g. Gueslin et al. Reference Gueslin, Talini, Herzhaft, Peysson and Allain2006; Putz et al. Reference Putz, Burghelea, Frigaard and Martinez2008; Holenberg et al. Reference Holenberg, Lavrenteva, Shavit and Nir2012). Here, we show that, superimposing the secondary cross-shear flow, the fore–aft asymmetry in the velocity field becomes less pronounced. Furthermore, the negative wake generated downstream of the sphere during sedimentation in an EVP material disappears in shear-induced sedimentation. The present 3-D numerical solver enables us to extract the yielded isosurfaces around the sphere. The yielded surface approaches the particle surface from the equatorial plane as $Bi$ increases, which eventually causes the sphere to stop settling.

We have also examined the total drag on the particle, along with its individual contributions. The total drag is calculated from the numerical data as the volumetric sum of the IB forces and the drag components computed by performing the numerical integration of the corresponding stress fields on the surface of the sphere. As previously reported by Fraggedakis et al. (Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016a), the drag coefficient for the settling in an EVP fluid increases when increasing plasticity and decreases when increasing elasticity. This trend is found to hold also in a sheared EVP fluid. In addition, the drag decreases considerably once the cross-shear flow is superimposed on the uniform flow for both Bingham and EVP material at constant plasticity and elasticity. The key finding is that the drag coefficient for a sphere settling in a sheared EVP fluid cannot be obtained from the drag coefficient pertaining a sphere settling in an otherwise quiescent fluid. This implies that the coupling of the cross-shear flow and the uniform flow is nonlinear and plays a major role in determining the sphere drag in yield stress fluids.

In this study, as the particle settles in a sheared yield stress fluid, the second invariant of the deformation rate tensor is increased by only $0.5\,\%$ compared to settling in the absence of cross-shear flow. This $0.5\,\%$ increase in the shear rate induces an approximately $5\,\%$ decrease in dynamic pressure, $4\,\%$ reduction of viscous normal stress, $\unicode[STIX]{x1D70F}_{yy}$, and $3\,\%$ and $16\,\%$ decrease in the $\unicode[STIX]{x1D70F}_{xy}$ and $\unicode[STIX]{x1D70F}_{yz}$ components of the viscous shear stress. Interestingly, however, the change in the dynamic pressure and the viscous stresses on the surface of the sphere settling through an EVP material is more than $18\,\%$ for the whole range of $Bi$ investigated (except for the change in viscous normal stress, $\unicode[STIX]{x1D70F}_{yy}$, at $Bi=1$ which is approximately $2\,\%$). Therefore, the coupling between the cross-shear flow and the uniform flow past a sphere affects the stress fields and the drag significantly.

By decomposing the total drag coefficient for a sphere settling in an otherwise quiescent EVP fluid into its components, it becomes evident that the form drag (resulting from the dynamic pressure on the particle surface) is the dominant component and the primary cause of drag enhancement with increasing Bingham number. Nevertheless, the viscous drag is the largest component in a Bingham fluid at very small $Bi$ number ($Bi=0.05$). The dynamic pressure and stress fields on the surface of the sphere are comprehensively analysed. We find that, for the EVP fluid, as $Bi$ increases, the magnitude of the dynamic pressure also increases while its symmetry around the sphere stagnation points breaks, resulting in further drag enhancement with the material plasticity. Adding a small degree of elasticity to the purely Bingham fluid modifies the viscous stresses on the surface of the sphere settling in both quiescent and sheared EVP fluids. However, the pressure distribution remains almost unaffected. The viscous stresses are modified so as to give viscous drag reduction. Consequently, the total drag is less in the case of the EVP fluid than for the Bingham fluid, regardless of the existence of an externally imposed cross-shear flow. Thus, adding polymers to an ideal yield stress fluid causes drag reduction through the modification of the viscous stresses on the surface of the particle settling in both quiescent and sheared EVP fluids.

We have also performed simulations of a sphere settling in a sheared EVP material at higher elasticity ($Wi_{\infty }=1$) and found that the drag coefficient is lower than that at lower elasticity ($Wi_{\infty }=0.1$). This drag reduction is mainly due to a decrease of both the form and viscous drag components, while the polymer drag remains almost unaffected. Moreover, the longer relaxation time of the macromolecular chains affects the normal viscous stresses, reducing their magnitude on the sphere surface. Conversely, the magnitude of the viscous shear stresses, $\unicode[STIX]{x1D70F}_{xy}$ and $\unicode[STIX]{x1D70F}_{yz}$, slightly increases at the particle surface; this is, however, not sufficient to overcome the decrease in viscous normal stress, $\unicode[STIX]{x1D70F}_{yy}$, so that the viscous drag is lower at higher $Wi$. Note also that the asymmetry in the dynamic pressure and viscous stress distributions is more pronounced at higher elasticity ($Wi_{\infty }=1$).

This study opens an avenue in answering many fundamental questions involving particles in practical yield stress fluids, i.e. the effects of shear thinning and solvent to polymer viscosity ratio on the drag, the drag laws when multibody interactions between particles are present, confinement effects, lubrication forces between particles, lift forces and particle migration when inertia becomes more relevant, etc. In addition, performing experiments and comparing the results with the simulations, simulations would help us to refine the constitutive laws of practical yield stress fluids as in reality thixotropy, elasticity and plasticity coexist in such materials. We hope this work will provide new insights to help tackling these challenging problems.

Acknowledgements

This work was supported by the National Science Foundation (grant CBET-1554044-CAREER), National Science Foundation’s Engineering Research Centers (grant CBET-1641152 Supplementary CAREER), American Chemical Society Petroleum Research Fund (grant 55661-DNI9). The authors acknowledge the computer time provided by the Ohio Supercomputer Center (OSC). L.B. and M.E.R. acknowledge the support from the European Research Council under grant ERC-2013-CoG-616186, TRITOS.

Declaration of interests

The authors report no conflict of interest.

Appendix A. Bingham fluid modelling

To model the stress-deformation behaviour of a viscoplastic fluid, the Bingham constitutive equation was first proposed (Bingham Reference Bingham1922; Bird, Dai & Yarusso Reference Bird, Dai and Yarusso1983). These relations are recovered by setting $Wi=0$, $\unicode[STIX]{x1D6FD}=1$ and $n=1$ in (2.1)–(2.3).

In our simulations, all boundary conditions remain the same as described in § 2.3 for the case of an EVP material, except the inlet condition where we use the components of the viscoplastic stress tensor. As mentioned in the text, these are obtained by analytically solving the Bingham fluid for the combined Couette and uniform flow at steady state in the absence of the particle. In this case, the only non-vanishing component of the stress tensor is the shear stress in the spanwise direction of the shear plane ($\unicode[STIX]{x1D70F}_{xz}$) which is obtained via

(A 1)

$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{xz}=2\unicode[STIX]{x1D6FC}+Bi.\end{eqnarray}$$

Note that contrary to an EVP material, the first normal stress difference is zero for the case of a Bingham fluid. Here, we perform the simulations at two $Bi$ numbers ($Bi=0.05,0.5$) for pure sedimentation and at one $Bi$ number ($Bi=0.05$) for shear-induced sedimentation of a single sphere, with $\unicode[STIX]{x1D6FC}=0$ and $0.1$, respectively.

A.1 Numerical Method

The Bingham fluid constitutive equation has an inherent discontinuity as the state of stress is undetermined before the material yielding point. To overcome the difficulty associated with the numerical treatment of the yield stress constitutive equation, we use the regularization method with the modification proposed by Papanastasiou (Reference Papanastasiou1987) where the exponential growth of the extra stress tensor is controlled by the material parameter, $m$; consequently, we can apply the same equation to the rigid and deformed regions of the medium. The viscous stress and the deformation rate tensors are related to each other via the apparent viscosity

(A 2)

$$\begin{eqnarray}\unicode[STIX]{x1D702}_{app}=\unicode[STIX]{x1D702}_{p}+\frac{\unicode[STIX]{x1D70F}_{0}}{|\dot{\unicode[STIX]{x1D6FE}}|}(1-\exp (-m|\dot{\unicode[STIX]{x1D6FE}}|)),\end{eqnarray}$$

where $|\dot{\unicode[STIX]{x1D6FE}}|$ is the second invariant of the shear rate tensor and $m$ is the stress growth exponent. The stress tensor is then $\unicode[STIX]{x1D749}=\unicode[STIX]{x1D702}_{app}\boldsymbol{\cdot }\unicode[STIX]{x1D738}$.

In order for this equation to mimic the ideal Bingham fluid, $m$ should be chosen to be sufficiently large (see e.g. Blackery & Mitsoulis Reference Blackery and Mitsoulis1997; Zisis & Mitsoulis Reference Zisis and Mitsoulis2002; Mitsoulis Reference Mitsoulis2004; Tokpavi et al. Reference Tokpavi, Magnin and Jay2008). It is worth mentioning that the yield surface location as predicted by the regularized viscosity function is a strong function of the regularization parameter $m$ and the behaviour of the true yield stress fluid is recovered when $(m\rightarrow \infty )$. Here, we have adopted a value of $1000$ for the regularization parameter $m$ except for the validation case when we chose a value of $200$ for consistency with the simulations performed by Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997). Moreover, we have checked that results for values of $m=500$ and $m=1000$ are in close agreement for the cases considered. The stress tensor components are advanced in time with the same Crank–Nicolson scheme used for the EVP fluid and the spatial derivatives computed using second-order central differences.

A.2 Validation case

We validate our implementation against the results by Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997) for the case of the uniform flow of Bingham fluid past a single sphere held stationary in a rectangular channel. The computational domain has a confinement ratio similar to Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997) with domain size $L_{x}=6D$, $L_{y}=16D$ and $L_{z}=4D$, discretized with $192\times 512\times 128$ points. Note that Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997) performed the simulations in a tube by adopting the axisymmetric boundary condition rather than a rectangular channel. The boundary conditions are described in § 2.3. In our simulation, the particle Reynolds number is chosen to be $1$ and the Bingham number $Bi=0.108$.

Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997) proposed a correlation for the drag of a sphere settling in a Bingham fluid as a function of Bingham number and confinement ratio. Defining the Stokes drag as

(A 3)

$$\begin{eqnarray}C_{s}=\frac{2F_{d}}{6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}_{p}U_{\infty }D},\end{eqnarray}$$

the following correlation is given by Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997):

(A 4)

$$\begin{eqnarray}C_{s}=C_{s,N}+aBi^{b},\end{eqnarray}$$

where $C_{s,N}$ is the Stokes drag coefficient for a sphere settling in a Newtonian fluid at the same confinement ratio, which can be obtained using Bohlin’s approximation (see e.g. Miyamura, Iwasaki & Ishii Reference Miyamura, Iwasaki and Ishii1981; Zheng, Phan-Thien & Tanner Reference Zheng, Phan-Thien and Tanner1991). The coefficients $a$ and $b$ depend on the tube to sphere radius ratio. For a confinement ratio of $4$, these constants take the value of $1.92$ and $0.92$, respectively. Table 7 compares the Stokes drag coefficient resulting from the present work with the predictions of (A 4), showing that the present results predict the Stokes drag coefficient with a relative difference of less than $0.5\,\%$.

Figure 23.

Comparison of the yielded/unyielded zones around a particle settling in a Bingham fluid at the confinement ratio of $4$ and at $Bi=0.108$; (a) computations of Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997) at $Re_{p}=0$, (b) present work at $Re_{p}=1$. White and black colours represent the yielded and unyielded zones, respectively.

Table 7.

Comparison of Stokes drag coefficient of a sphere settling in a Bingham fluid resulting from the present work (second column) and the computations of Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997).

In figure 23, we compare the size and shape of the yielded/unyielded zones around the sphere with the data from Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997).

The agreement between our results at $Re_{p}=1$ and the computations of Blackery & Mitsoulis (Reference Blackery and Mitsoulis1997) at $Re_{p}=0$ is quite satisfactory. The slight zigzagging of the yielded surface boundary results from adopting a large value for the regularization parameter ($m=200$) in the Bingham fluid constitutive equation, which generates non-smooth lines in the solution. For more details, the reader is referred to Zisis & Mitsoulis (Reference Zisis and Mitsoulis2002). As the Stokes drag coefficient and the size and shape of the yielded/unyielded zones at $Re_{p}=1$ are in very good agreement with the results for creeping flow conditions, we expect the effect of inertia to be negligible at $Re_{p}=1$. See § 3.1.3 for further discussion.

Appendix B. Numerical integration procedure

In this section, the surface integration procedure implemented to evaluate the integrals in (3.4)–(3.6) along with the interpolation scheme used to compute the pressure and stress fields on the surface of the particle are described in detail.

In the IB method, the first prediction velocity is interpolated from a Eulerian grid cell (used for the fluid phase) to the Lagrangian grid (used for the particle phase). Then, the IB force, computed on the Lagrangian grid from the difference between the actual particle velocity of each Lagrangian force point and the interpolated first prediction velocity, is interpolated back to the Eulerian grid. These are called ‘interpolation’ and ‘spreading’ operations (Peskin Reference Peskin2002; Uhlmann Reference Uhlmann2005; Breugem Reference Breugem2012) and are performed using the regularized Dirac delta function $(\unicode[STIX]{x1D739}_{\boldsymbol{d}})$ from Roma, Peskin & Berger (Reference Roma, Peskin and Berger1999). This function, which extends over three grid cells in each coordinate direction, is approximated in three-dimensional space by the following product:

(B 1)

$$\begin{eqnarray}\unicode[STIX]{x1D739}_{\boldsymbol{d}}^{3}(\boldsymbol{x}_{ijk}-\boldsymbol{X}_{l})=\unicode[STIX]{x1D6FF}_{d}^{1}(x_{ijk}-X_{l})\unicode[STIX]{x1D6FF}_{d}^{1}(\,y_{ijk}-Y_{l})\unicode[STIX]{x1D6FF}_{d}^{1}(z_{ijk}-Z_{l}),\end{eqnarray}$$

where $\boldsymbol{x}_{ijk}$ and $\boldsymbol{X}_{l}$ denote the Eulerian grid point with index $(i,j,k)$ and the Lagrangian marker point position with index $l$. Here, $\unicode[STIX]{x1D6FF}_{d}^{1}(x_{ijk}-X_{l})$ is an one-dimensional approximation of the delta function

(B 2)

$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{d}^{1}(x_{ijk}-X_{l})=\frac{1}{\unicode[STIX]{x0394}x}\unicode[STIX]{x1D719}\left(\frac{x_{ijk}-X_{l}}{\unicode[STIX]{x0394}x}\right).\end{eqnarray}$$

In (B 2), $\unicode[STIX]{x0394}x$ is the Eulerian grid size and $\unicode[STIX]{x1D719}$ is a continuous function that is chosen to satisfy the discrete version of the Dirac delta function properties and is obtained from the following equation (Roma et al. Reference Roma, Peskin and Berger1999):

(B 3)

$$\begin{eqnarray}\unicode[STIX]{x1D719}(r)=\left\{\begin{array}{@{}ll@{}}\frac{1}{6}(5-3|r|-\sqrt{-3(1-|r|)^{2}+1}), & 0.5\leqslant |r|\leqslant 1.5,\\[4.0pt] \frac{1}{3}(1+\sqrt{-3r^{2}+1}), & |r|\leqslant 0.5,\\[4.0pt] 0, & \text{otherwise},\end{array}\right.\end{eqnarray}$$

with $r=(x_{ijk}-X_{l})/\unicode[STIX]{x0394}x$. The one-dimensional approximations for the delta function in the $y$ and $z$ directions (i.e. $\unicode[STIX]{x1D6FF}_{d}^{1}(\,y_{ijk}-Y_{l})$, $\unicode[STIX]{x1D6FF}_{d}^{1}(z_{ijk}-Z_{l})$) are defined similarly to (B 2) by replacing the grid spacing, and Eulerian and Lagrangian grid point positions in the corresponding coordinate direction.

If the dimension of the Eulerian grid cell in each direction is spatially uniform then the regularized Dirac delta function proposed by Roma et al. (Reference Roma, Peskin and Berger1999) ensures that the total hydrodynamic force and torque that the fluid and particles exert onto each other are preserved in the interpolation and spreading operations (Uhlmann Reference Uhlmann2003; Breugem Reference Breugem2012). Thus, in the present work the Eulerian grid is considered to be a Cartesian grid with uniform size in each coordinate direction, i.e. $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=\unicode[STIX]{x0394}z$. For improved accuracy, the Lagrangian points should be uniformly distributed all over the surface of the sphere, with a spacing of the order of the Eulerian grid. In the present work, $3219$ Lagrangian points are therefore used to match the Eulerian grid resolution.

Since the Lagrangian cells are evenly distributed over the surface of the sphere, the surface of the sphere is partitioned into regions of equal area in such a way that the unity between domains is null, while the union of them constitutes the entire sphere surface area. We use the quadrature rule to evaluate the surface integral (see e.g. Atkinson Reference Atkinson1982; An & Chen Reference An and Chen2016; Reeger & Fornberg Reference Reeger and Fornberg2016). The centre position of each Lagrangian partition is chosen as the quadrature node. Hence, the weight of the corresponding quadrature rule is positive and equal, and the integral of any quantities over the surface of the sphere can be estimated by the following equation:

(B 4)

$$\begin{eqnarray}\iint _{S}fn_{j}\,\text{d}S\approx \mathop{\sum }_{l=1}^{N_{l}}f_{l}n_{j,l}S_{l},\end{eqnarray}$$

where $f_{l}$ can be any quantity (e.g. pressure, viscous and polymer stresses) of the $l\text{th}$ force point, $l$ is the Lagrangian point index, $N_{l}$ is the total number of Lagrangian grid points (total number of equally partitioned subareas on the surface), $n_{j,l}$ is the unit normal vector directed outward on the $l\text{th}$ Lagrangian point of the sphere. Here, $S_{l}$ is the surface area of the $l\text{th}$ partition which is equal to $(4\unicode[STIX]{x03C0}R^{2})/N_{l}$ for a sphere with radius $R$.

Therefore, we only need the pressure, viscous/polymer stresses and the unit normal vector on each Lagrangian force point to compute the surface integrals in (3.4)–(3.6) by using (B 4) to find the drag components.

Computing the form, viscous and polymer drag components is performed in several steps. In the first step, the solid volume fraction ($\unicode[STIX]{x1D6FC}_{ijk}$) in the grid cell with index $(i,j,k)$ around the collocation points of the velocities ($u,v,w$) and the pressure ($p$) are computed. This should be done at each Runge–Kutta step.

The solid volume fraction $\unicode[STIX]{x1D6FC}_{ijk}$ in the Eulerian cubic grid cell with index $(i,j,k)$ is determined from the level-set function $\unicode[STIX]{x1D713}$ obtained by calculating the signed distance of each Eulerian grid point to the particle surface $S$ (Kempe, Schwarz & Fröhlich Reference Kempe, Schwarz and Fröhlich2009),

(B 5)

$$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{ijk}=\frac{\mathop{\sum }_{k=1}^{8}-\unicode[STIX]{x1D713}_{k}H(-\unicode[STIX]{x1D713}_{k})}{\mathop{\sum }_{k=1}^{8}|\unicode[STIX]{x1D713}_{k}|}.\end{eqnarray}$$

In (B 5), the sum is over all eight corner nodes of the cubical Eulerian grid cell and $H$ is the Heaviside step function. The signed distance function $\unicode[STIX]{x1D713}$ represents the position of the Eulerian grid cell with respect to the sphere with $\unicode[STIX]{x1D713}>0$ outside and $\unicode[STIX]{x1D713}<0$ inside the particle. Since we use a fully staggered Cartesian grid, the velocity components are computed at the cell faces while the pressure and the stress fields are calculated at the cell vertex. Hence, the collocation points of the $\unicode[STIX]{x1D701}$-component of the velocity collocation points in the $\unicode[STIX]{x1D6FE}$ direction and the collocation points for the pressure are given by the following equations:

(B 6)

$$\begin{eqnarray}\displaystyle & \displaystyle x_{\unicode[STIX]{x1D6FE}}^{\unicode[STIX]{x1D701}}(k)=(k-1-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D701}})\unicode[STIX]{x0394}x_{\unicode[STIX]{x1D6FE}};\quad k=1,2,\ldots ,n_{\unicode[STIX]{x1D6FE}}+\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D701}}, & \displaystyle\end{eqnarray}$$

(B 7)

$$\begin{eqnarray}\displaystyle & \displaystyle x_{\unicode[STIX]{x1D6FE}}^{p}(k)=(k-1)\unicode[STIX]{x0394}x_{\unicode[STIX]{x1D6FE}};\quad k=1,2,\ldots ,n_{\unicode[STIX]{x1D6FE}}. & \displaystyle\end{eqnarray}$$

In (B 6) and (B 7), $\unicode[STIX]{x0394}x_{\unicode[STIX]{x1D6FE}}$ is the size of the grid cell in the $\unicode[STIX]{x1D6FE}$ direction and $n_{\unicode[STIX]{x1D6FE}}$ is the number of pressure grid points in the $\unicode[STIX]{x1D6FE}$ direction. The use of the staggered grid results in different solid volume fractions depending on the variable considered. Figure 24 depicts three different solid volume fractions (highlighted areas) for the grid cells around $p(i,j)$, $u(i+9/2,j+2)$ and $v(i+2,j+1/2)$. The sphere boundary is shown by the red dashed line.

Figure 24.

Solid volume fractions (highlighted area) for grid cells around $p(i,j)$,$u(i+9/2,j+2)$ and $v(i+2,j+1/2)$. The sphere boundary is shown by the red dashed line.

Finally, the fluid velocity in each coordinate direction ($u,v,w$) and the pressure at each Runge–Kutta step are obtained by

(B 8)

$$\begin{eqnarray}\displaystyle & \displaystyle u_{ijk}^{\,fq}=(1-\unicode[STIX]{x1D6FC}_{ijk}^{u})u_{ijk}^{q}, & \displaystyle\end{eqnarray}$$

(B 9)

$$\begin{eqnarray}\displaystyle & \displaystyle v_{ijk}^{\,fq}=(1-\unicode[STIX]{x1D6FC}_{ijk}^{v})v_{ijk}^{q}, & \displaystyle\end{eqnarray}$$

(B 10)

$$\begin{eqnarray}\displaystyle & \displaystyle w_{ijk}^{\,fq}=(1-\unicode[STIX]{x1D6FC}_{ijk}^{w})w_{ijk}^{q}, & \displaystyle\end{eqnarray}$$

(B 11)

$$\begin{eqnarray}\displaystyle & \displaystyle p_{ijk}^{\,fq}=(1-\unicode[STIX]{x1D6FC}_{ijk}^{p})p_{ijk}^{q}, & \displaystyle\end{eqnarray}$$

where $u_{ijk}^{\,fq}$, $v_{ijk}^{\,fq}$, $w_{ijk}^{\,fq}$ and $p_{ijk}^{\,fq}$ are the velocity components in the $x$, $y$ and $z$ directions and the pressure after direct accounting for the inertia of the fluid contained within the sphere at the Runge–Kutta step $q$. Here, $\unicode[STIX]{x1D6FC}_{ijk}^{u}$, $\unicode[STIX]{x1D6FC}_{ijk}^{v}$, $\unicode[STIX]{x1D6FC}_{ijk}^{w}$ are the solid volume fractions around the velocity collocation points in the $x$, $y$ and $z$ directions (Eulerian cell face) and $\unicode[STIX]{x1D6FC}_{ijk}^{p}$ is the solid volume fraction around the pressure collocation point that is at the Eulerian grid cell vertex (see figure 24). Here, $u_{ijk}^{q}$, $v_{ijk}^{q}$, $w_{ijk}^{q}$ and $p_{ijk}^{q}$ are the velocity components in the $x$, $y$ and $z$ directions and the pressure before direct accounting of the inertia of the fluid in the volume of the particle at the $q\text{th}$ Runge–Kutta step. It should be noted that, as the polymer stresses are computed at the cell centre, the same solid volume fraction that is obtained around the pressure collocation point is used for the polymer stress components.

In the second step, the viscous stresses are computed on each grid cell

(B 12)

$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{ij}^{v}=(1-\unicode[STIX]{x1D6FD})\left(\frac{\unicode[STIX]{x2202}u_{i}^{\,fq}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}^{\,fq}}{\unicode[STIX]{x2202}x_{i}}\right);\quad i,j=1,2,3,\end{eqnarray}$$

where the spatial derivatives are estimated with the central differencing scheme.

So far the pressure, viscous and polymer stresses are computed at each Eulerian grid cell with index $(i,j,k)$ considering the solid volume fraction around the pressure and velocity component collocation points. In the third step, the pressure, viscous and polymer stresses should be projected to the corresponding Lagrangian grid cell. This interpolation is done using the same Dirac delta function defined in (B 1)–(B 3):

(B 13)

$$\begin{eqnarray}\displaystyle & \displaystyle P_{l}=\mathop{\sum }_{ijk}^{}p_{ijk}^{\,fq}\unicode[STIX]{x1D739}_{\boldsymbol{d}}^{3}(\boldsymbol{x}_{ijk}-\boldsymbol{X}_{l})\unicode[STIX]{x0394}x\unicode[STIX]{x0394}y\unicode[STIX]{x0394}z;\quad l=1,2,\ldots ,N_{l}, & \displaystyle\end{eqnarray}$$

(B 14)

$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{pq}^{v,l}=\mathop{\sum }_{ijk}^{}\unicode[STIX]{x1D70F}_{pq}^{v,ijk}\unicode[STIX]{x1D739}_{\boldsymbol{d}}^{3}(\boldsymbol{x}_{ijk}-\boldsymbol{X}_{l})\unicode[STIX]{x0394}x\unicode[STIX]{x0394}y\unicode[STIX]{x0394}z;\quad l=1,2,\ldots ,N_{l};\quad p,q=1,2,3, & \displaystyle\end{eqnarray}$$

(B 15)

$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{pq}^{Po,l}=\mathop{\sum }_{ijk}^{}\unicode[STIX]{x1D70F}_{pq}^{Po,ijk}\unicode[STIX]{x1D739}_{\boldsymbol{d}}^{3}(\boldsymbol{x}_{ijk}-\boldsymbol{X}_{l})\unicode[STIX]{x0394}x\unicode[STIX]{x0394}y\unicode[STIX]{x0394}z;\quad l=1,2,\ldots ,N_{l};\quad p,q=1,2,3, & \displaystyle\end{eqnarray}$$

where $P_{l}$, $\unicode[STIX]{x1D70F}_{pq}^{v,l}$ and $\unicode[STIX]{x1D70F}_{pq}^{Po,l}$ are the interpolated pressure, viscous and polymer stress components at the Lagrangian point $l$.

In the last step, we compute the unit normal vector on the $l\text{th}$ Lagrangian point. This is easily done once the position of the centre of the particle is known

(B 16)

$$\begin{eqnarray}n_{j,l}=\frac{\unicode[STIX]{x1D735}G}{|\unicode[STIX]{x1D735}G|},\end{eqnarray}$$

where

(B 17)

$$\begin{eqnarray}G=(X_{l}-X_{c})^{2}+(Y_{l}-Y_{c})^{2}+(Z_{l}-Z_{c})^{2}-R^{2};\quad l=1,2,\ldots ,N_{l},\end{eqnarray}$$

and $X_{c}$, $Y_{c}$ and $Z_{c}$ are the coordinates of the sphere centre. Finally, we are able to estimate the pressure, viscous and polymer stresses at the Lagrangian points. Given the unit normal vector, the form, viscous and polymer drag components (3.4)–(3.6) are found from the following relations:

(B 18)

$$\begin{eqnarray}\displaystyle & \displaystyle \iint _{S}pn_{y}\,\text{d}S\approx \frac{4\unicode[STIX]{x03C0}R^{2}}{N_{l}}\mathop{\sum }_{l=1}^{N_{l}}P_{l}n_{y,l}, & \displaystyle\end{eqnarray}$$

(B 19)

$$\begin{eqnarray}\displaystyle & \displaystyle (1-\unicode[STIX]{x1D6FD})\iint _{S}\left(\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}x_{q}}+\frac{\unicode[STIX]{x2202}u_{q}}{\unicode[STIX]{x2202}y}\right)n_{q}\,\text{d}S\approx (1-\unicode[STIX]{x1D6FD})\frac{4\unicode[STIX]{x03C0}R^{2}}{N_{l}}\mathop{\sum }_{l=1}^{N_{l}}\unicode[STIX]{x1D70F}_{yq}^{v,l}n_{q,l}, & \displaystyle\end{eqnarray}$$

(B 20)

$$\begin{eqnarray}\displaystyle & \displaystyle \iint _{S}\unicode[STIX]{x1D70F}_{yq}n_{q}\,\text{d}S\approx \frac{4\unicode[STIX]{x03C0}R^{2}}{N_{l}}\mathop{\sum }_{l=1}^{N_{l}}\unicode[STIX]{x1D70F}_{yq}^{Po,l}n_{q,l}. & \displaystyle\end{eqnarray}$$

B.1 Validation case

In order to validate the numerical integration scheme that we have implemented in this study, we compute the drag and its components in the case of the Newtonian fluid past a sphere in a finite domain. For this test, the domain size is $L_{x}=6D$, $L_{y}=8D$ and $L_{z}=5D$, with inflow–outflow boundary conditions applied in the streamwise $y$ direction, no-slip enforced in the wall-normal $z$ direction and the periodic boundary conditions in the spanwise $x$ direction. The simulation is performed at $Re_{p}=1$.

In this simulation the total Stokes drag along with the form and viscous drag components are computed and compared with the theoretical solution (Leal Reference Leal2007). Confinement effects are taken into account with the wall-correction factor given by Miyamura et al. (Reference Miyamura, Iwasaki and Ishii1981). Note that at the confinement ratio of $0.2$, the wall-correction given by Miyamura et al. (Reference Miyamura, Iwasaki and Ishii1981) converges to the Faxen law (Faxén Reference Faxén1922). A comparison between the present results and the Faxén (Reference Faxén1922) approximation is presented in table 8, where the Stokes drag coefficient

(B 21)

$$\begin{eqnarray}C_{s}=\frac{F_{d}}{6\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}U_{\infty }R},\end{eqnarray}$$

with $\unicode[STIX]{x1D702}$ the viscosity of the Newtonian fluid.

Table 8.

Comparison of drag coefficient and drag components from the present work, PW, with Bohlin’s approximation and theoretical solutions for Newtonian fluids.

The agreement in the total Stokes drag is satisfactory with a relative difference of approximately 1.5 %. The relative difference in the form drag and the viscous drag is around 2 % and 3 %, respectively. The small discrepancy between the numerical integration of the dynamic pressure and the stress fields on the surface of the sphere and the theoretical predictions is primarily due to the Lagrangian grid resolution. Increasing the number of Lagrangian force points, the integration scheme becomes more accurate and converges towards the theoretical solution.

References

Ahonguio, F., Jossic, L. & Magnin, A. 2014 Influence of surface properties on the flow of a yield stress fluid around spheres. J. Non-Newtonian Fluid Mech. 206, 57–70.CrossRef Google Scholar

Ahonguio, F., Jossic, L. & Magnin, A. 2015 Motion and stability of cones in a yield stress fluid. AIChE J. 61 (2), 709–717.CrossRef Google Scholar

Alghalibi, D., Lashgari, I., Brandt, L. & Hormozi, S. 2018 Interface-resolved simulations of particle suspensions in Newtonian, shear thinning and shear thickening carrier fluids. J. Fluid Mech. 852, 329–357.CrossRef Google Scholar

An, C. & Chen, S.2016 Numerical integration over the unit sphere by using spherical

$t$-design. arXiv:1611.02785.Google Scholar

Andres, U. T. 1960 Equilibrium and motion of spheres in a viscoplastic liquid. In Dokl. Akad. Nauk, vol. 133, pp. 777–780. Russian Academy of Sciences.Google Scholar

Ansley, R. W. & Smith, T. N. 1967 Motion of spherical particles in a Bingham plastic. AIChE J. 13 (6), 1193–1196.CrossRef Google Scholar

Atapattu, D. D., Chhabra, R. P. & Uhlherr, P. H. T. 1995 Creeping sphere motion in Herschel–Bulkley fluids: flow field and drag. J. Non-Newtonian Fluid Mech. 59 (2–3), 245–265.CrossRef Google Scholar

Atkinson, K. 1982 Numerical integration on the sphere. ANZIAM J. 23 (3), 332–347.Google Scholar

Balmforth, N. J., Craster, R. V., Hewitt, D. R., Hormozi, S. & Maleki, A. 2017 Viscoplastic boundary layers. J. Fluid Mech. 813, 929–954.CrossRef Google Scholar

Beaulne, M. & Mitsoulis, E. 1997 Creeping motion of a sphere in tubes filled with Herschel–Bulkley fluids. J. Non-Newtonian Fluid Mech. 72 (1), 55–71.CrossRef Google Scholar

Bénito, S., Bruneau, C.-H., Colin, T., Gay, C. & Molino, F. 2008 An elasto-visco-plastic model for immortal foams or emulsions. Eur. Phys. J. E 25 (3), 225–251.Google Scholar PubMed

Beris, A. N., Tsamopoulos, J. A., Armstrong, R. C. & Brown, R. A. 1985 Creeping motion of a sphere through a Bingham plastic. J. Fluid Mech. 158, 219–244.CrossRef Google Scholar

Bingham, E. C. 1922 Fluidity and Plasticity, vol. 2. McGraw-Hill.Google Scholar

Bird, R. B., Dai, G. C. & Yarusso, B. J. 1983 The rheology and flow of viscoplastic materials. Rev. Chem. Engng 1 (1), 1–70.CrossRef Google Scholar

Blackery, J. & Mitsoulis, E. 1997 Creeping motion of a sphere in tubes filled with a Bingham plastic material. J. Non-Newtonian Fluid Mech. 70 (1–2), 59–77.CrossRef Google Scholar

Boardman, G. & Whitmore, R. L. 1960 Yield stress exerted on a body immersed in a Bingham fluid. Nature 187 (4731), 50–51.CrossRef Google Scholar

Breugem, W.-P. 2012 A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows. J. Comput. Phys. 231 (13), 4469–4498.CrossRef Google Scholar

van den Brule, B. H. A. A. & Gheissary, G. 1993 Effects of fluid elasticity on the static and dynamic settling of a spherical particle. J. Non-Newtonian Fluid Mech. 49 (1), 123–132.CrossRef Google Scholar

Brunn, P. 1977a Errata to the slow motion of a sphere in a second-order fluid. Rheol. Acta 16 (3), 324–325.CrossRef Google Scholar

Brunn, P. 1977b The slow motion of a rigid particle in a second-order fluid. J. Fluid Mech. 82 (3), 529–547.CrossRef Google Scholar

Chaparian, E. & Frigaard, I. A. 2017a Cloaking: particles in a yield-stress fluid. J. Non-Newtonian Fluid Mech. 243, 47–55.CrossRef Google Scholar

Chaparian, E. & Frigaard, I. A. 2017b Yield limit analysis of particle motion in a yield-stress fluid. J. Fluid Mech. 819, 311–351.CrossRef Google Scholar

Cheddadi, I. & Saramito, P. 2013 A new operator splitting algorithm for elastoviscoplastic flow problems. J. Non-Newtonian Fluid Mech. 202, 13–21.CrossRef Google Scholar

Cheddadi, I., Saramito, P., Dollet, B., Raufaste, C. & Graner, F. 2011 Understanding and predicting viscous, elastic, plastic flows. Eur. Phys. J. E 34 (1), 1.Google Scholar PubMed

Chen, S., Phan-Thien, N., Khoo, B. C. & Fan, X. J. 2006 Flow around spheres by dissipative particle dynamics. Phys. Fluids 18 (10), 103605.CrossRef Google Scholar

Coussot, P., Raynaud, J. S., Bertrand, F., Moucheront, P., Guilbaud, J. P., Huynh, H. T., Jarny, S. & Lesueur, D. 2002 Coexistence of liquid and solid phases in flowing soft-glassy materials. Phys. Rev. Lett. 88 (21), 218301.CrossRef Google Scholar PubMed

De Besses, B. D., Magnin, A. & Jay, P. 2003 Viscoplastic flow around a cylinder in an infinite medium. J. Non-Newtonian Fluid Mech. 115 (1), 27–49.CrossRef Google Scholar

De Vita, F., Rosti, M. E., Izbassarov, D., Duffo, L., Tammisola, O., Hormozi, S. & Brandt, L. 2018 Elastoviscoplastic flows in porous media. J. Non-Newtonian Fluid Mech. 258, 10–21.CrossRef Google Scholar

Dimitriou, C. J., Ewoldt, R. H. & McKinley, G. H. 2013 Describing and prescribing the constitutive response of yield stress fluids using large amplitude oscillatory shear stress (LAOStress). J. Rheol. 57 (1), 27–70.CrossRef Google Scholar

Dimitriou, C. J. & McKinley, G. H. 2014 A comprehensive constitutive law for waxy crude oil: a thixotropic yield stress fluid. Soft Matt. 10 (35), 6619–6644.CrossRef Google Scholar PubMed

Dollet, B. & Graner, F. 2007 Two-dimensional flow of foam around a circular obstacle: local measurements of elasticity, plasticity and flow. J. Fluid Mech. 585, 181–211.CrossRef Google Scholar

Einarsson, J. & Mehlig, B. 2017 Spherical particle sedimenting in weakly viscoelastic shear flow. Phys. Rev. Fluids 2 (6), 063301.CrossRef Google Scholar

Faxén, H. 1922 Der Widerstand gegen die Bewegung einer starren Kugel in einer zähen Flüssigkeit, die zwischen zwei parallelen ebenen Wänden eingeschlossen ist. Ann. Phys. 373 (10), 89–119.CrossRef Google Scholar

Firouznia, M., Metzger, B., Ovarlez, G. & Hormozi, S. 2018 The interaction of two spherical particles in simple-shear flows of yield stress fluids. J. Non-Newtonian Fluid Mech. 255, 19–38.CrossRef Google Scholar

Fraggedakis, D., Dimakopoulos, Y. & Tsamopoulos, J. 2016a Yielding the yield-stress analysis: a study focused on the effects of elasticity on the settling of a single spherical particle in simple yield-stress fluids. Soft Matt. 12 (24), 5378–5401.CrossRef Google Scholar

Fraggedakis, D., Pavlidis, M., Dimakopoulos, Y. & Tsamopoulos, J. 2016b On the velocity discontinuity at a critical volume of a bubble rising in a viscoelastic fluid. J. Fluid Mech. 789, 310–346.CrossRef Google Scholar

Frigaard, I. A. & Nouar, C. 2005 On the usage of viscosity regularisation methods for visco-plastic fluid flow computation. J. Non-Newtonian Fluid Mech. 127 (1), 1–26.CrossRef Google Scholar

Geri, M., Venkatesan, R., Sambath, K. & McKinley, G. H. 2017 Thermokinematic memory and the thixotropic elasto-viscoplasticity of waxy crude oils. J. Rheol. 61 (3), 427–454.CrossRef Google Scholar

Glowinski, R. & Wachs, A. 2011 On the numerical simulation of viscoplastic fluid flow. In Handbook of Numerical Analysis, vol. 16, pp. 483–717. Elsevier.Google Scholar

Gueslin, B., Talini, L., Herzhaft, B., Peysson, Y. & Allain, C. 2006 Flow induced by a sphere settling in an aging yield-stress fluid. Phys. Fluids 18 (10), 103101.CrossRef Google Scholar

Hariharaputhiran, M., Subramanian, R. S., Campbell, G. A. & Chhabra, R. P. 1998 The settling of spheres in a viscoplastic fluid. J. Non-Newtonian Fluid Mech. 79 (1), 87–97.CrossRef Google Scholar

Harlen, O. G. 2002 The negative wake behind a sphere sedimenting through a viscoelastic fluid. J. Non-Newtonian Fluid Mech. 108 (1-3), 411–430.CrossRef Google Scholar

Herschel, W. H. & Bulkley, R. 1926 Measurement of consistency as applied to rubber-benzene solutions. In Am. Soc. Test Proc, vol. 26, pp. 621–633. American Society for Testing Materials Proceeding (ASTM Proceeding).Google Scholar

Holenberg, Y., Lavrenteva, O. M., Shavit, U. & Nir, A. 2012 Particle tracking velocimetry and particle image velocimetry study of the slow motion of rough and smooth solid spheres in a yield-stress fluid. Phys. Rev. E 86 (6), 066301.Google Scholar

Housiadas, K. D. 2014 Stress diffusion and high order viscoelastic effects in the 3D flow past a sedimenting sphere subject to orthogonal shear. Rheol. Acta 53 (7), 537–548.CrossRef Google Scholar

Housiadas, K. D. & Tanner, R. I. 2012 The drag of a freely sendimentating sphere in a sheared weakly viscoelastic fluid. J. Non-Newtonian Fluid Mech. 183, 52–56.CrossRef Google Scholar

Housiadas, K. D. & Tanner, R. I. 2014 Rheological effects in the 3D creeping flow past a sedimenting sphere subject to orthogonal shear. Phys. Fluids 26 (1), 013102.CrossRef Google Scholar

Izbassarov, D., Rosti, M. E., Ardekani, M. N., Sarabian, M., Hormozi, S., Brandt, L. & Tammisola, O. 2018 Computational modeling of multiphase viscoelastic and elastoviscoplastic flows. Intl J. Numer. Meth. Fluids 88 (12), 521–543.CrossRef Google Scholar

Janiaud, E. & Graner, F. 2005 Foam in a two-dimensional Couette shear: a local measurement of bubble deformation. J. Fluid Mech. 532, 243–267.CrossRef Google Scholar

Jossic, L. & Magnin, A. 2001 Drag and stability of objects in a yield stress fluid. AIChE J. 47 (12), 2666–2672.CrossRef Google Scholar

Kempe, T., Schwarz, S. & Fröhlich, J. 2009 Modelling of spheroidal particles in viscous flows. In Proceedings of the Academy Colloquium Immersed Boundary Methods: Current Status and Future Research Directions (KNAW, Amsterdam, The Netherlands, 15–17 June 2009), vol. 845.Google Scholar

Kim, J. & Moin, P. 1985 Application of a fractional-step method to incompressible Navier–Stokes equations. J. Comput. Phys. 59 (2), 308–323.CrossRef Google Scholar

Lashgari, I., Picano, F., Breugem, W.-P. & Brandt, L. 2014 Laminar, turbulent, and inertial shear-thickening regimes in channel flow of neutrally buoyant particle suspensions. Phys. Rev. Lett. 113 (25), 254502.CrossRef Google Scholar PubMed

Leal, L. G. 2007 Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes, vol. 7. Cambridge University Press.CrossRef Google Scholar

Liu, B. T., Muller, S. J. & Denn, M. M. 2002 Convergence of a regularization method for creeping flow of a Bingham material about a rigid sphere. J. Non-Newtonian Fluid Mech. 102 (2), 179–191.CrossRef Google Scholar

Liu, Y., Balmforth, N. J., Hormozi, S. & Hewitt, D. R. 2016 Two-dimensional viscoplastic dambreaks. J. Non-Newtonian Fluid Mech. 238, 65–79.CrossRef Google Scholar

Lunsmann, W. J., Genieser, L., Armstrong, R. C. & Brown, R. A. 1993 Finite element analysis of steady viscoelastic flow around a sphere in a tube: calculations with constant viscosity models. J. Non-Newtonian Fluid Mech. 48 (1-2), 63–99.CrossRef Google Scholar

Maleki, A., Hormozi, S., Roustaei, A. & Frigaard, I. A. 2015 Macro-size drop encapsulation. J. Fluid Mech. 769, 482–521.CrossRef Google Scholar

Meeker, S. P., Bonnecaze, R. T. & Cloitre, M. 2004 Slip and flow in soft particle pastes. Phys. Rev. Lett. 92 (19), 198302.CrossRef Google Scholar PubMed

Merkak, O., Jossic, L. & Magnin, A. 2009 Migration and sedimentation of spherical particles in a yield stress fluid flowing in a horizontal cylindrical pipe. AIChE J. 55 (10), 2515–2525.CrossRef Google Scholar

Mitsoulis, E. 2004 On creeping drag flow of a viscoplastic fluid past a circular cylinder: wall effects. Chem. Engng Sci. 59 (4), 789–800.CrossRef Google Scholar

Miyamura, A., Iwasaki, S. & Ishii, T. 1981 Experimental wall correction factors of single solid spheres in triangular and square cylinders, and parallel plates. Intl J. Multiphase Flow 7 (1), 41–46.CrossRef Google Scholar

Murch, W. L., Krishnan, S., Shaqfeh, E. SG. & Iaccarino, G. 2017 Growth of viscoelastic wings and the reduction of particle mobility in a viscoelastic shear flow. Phys. Rev. Fluids 2 (10), 103302.CrossRef Google Scholar

Nirmalkar, N., Chhabra, R. P. & Poole, R. J. 2012 On creeping flow of a Bingham plastic fluid past a square cylinder. J. Non-Newtonian Fluid Mech. 171, 17–30.CrossRef Google Scholar

Ouattara, Z., Jay, P., Blésès, D. & Magnin, A. 2018 Drag of a cylinder moving near a wall in a yield stress fluid. AIChE J. 64 (11), 4118–4130.CrossRef Google Scholar

Ovarlez, G., Barral, Q. & Coussot, P. 2010 Three-dimensional jamming and flows of soft glassy materials. Nat. Mater. 9 (2), 115.CrossRef Google Scholar PubMed

Ovarlez, G., Bertrand, F., Coussot, P. & Chateau, X. 2012 Shear-induced sedimentation in yield stress fluids. J. Non-Newtonian Fluid Mech. 177, 19–28.CrossRef Google Scholar

Ovarlez, G. & Hormozi, S. 2019 Lectures on Visco-Plastic Fluid Mechanics. Springer.CrossRef Google Scholar

Padhy, S., Rodriguez, M., Shaqfeh, E. S. G., Iaccarino, G., Morris, J. F. & Tonmukayakul, N. 2013a The effect of shear thinning and walls on the sedimentation of a sphere in an elastic fluid under orthogonal shear. J. Non-Newtonian Fluid Mech. 201, 120–129.CrossRef Google Scholar

Padhy, S., Shaqfeh, E. S. G., Iaccarino, G., Morris, J. F. & Tonmukayakul, N. 2013b Simulations of a sphere sedimenting in a viscoelastic fluid with cross shear flow. J. Non-Newtonian Fluid Mech. 197, 48–60.CrossRef Google Scholar

Pantokratoras, A. 2018 Flow past a rotating sphere in a Bingham plastic fluid, up to a Reynolds number of 10 000. Rheol. Acta 57 (8-9), 611–617.CrossRef Google Scholar

Papanastasiou, T. C. 1987 Flows of materials with yield. J. Rheol. 31 (5), 385–404.CrossRef Google Scholar

Peskin, C. S. 2002 The immersed boundary method. Acta Numerica 11, 479–517.CrossRef Google Scholar

Putz, A. M. V., Burghelea, T. I., Frigaard, I. A. & Martinez, D. M. 2008 Settling of an isolated spherical particle in a yield stress shear thinning fluid. Phys. Fluids 20 (3), 033102.CrossRef Google Scholar

Putz, A. & Frigaard, I. A. 2010 Creeping flow around particles in a Bingham fluid. J. Non-Newtonian Fluid Mech. 165 (5–6), 263–280.CrossRef Google Scholar

Reeger, J. A. & Fornberg, B. 2016 Numerical quadrature over the surface of a sphere. Stud. Appl. Maths 137 (2), 174–188.CrossRef Google Scholar

Roma, A. M., Peskin, C. S. & Berger, M. J. 1999 An adaptive version of the immersed boundary method. J. Comput. Phys. 153 (2), 509–534.CrossRef Google Scholar

Roquet, N. & Saramito, P. 2003 An adaptive finite element method for Bingham fluid flows around a cylinder. Comput. Meth. Appl. Mech. Engng 192 (31-32), 3317–3341.CrossRef Google Scholar

Rosti, M. E. & Brandt, L. 2017 Numerical simulation of turbulent channel flow over a viscous hyper-elastic wall. J. Fluid Mech. 830, 708–735.CrossRef Google Scholar

Rosti, M. E., Izbassarov, D., Tammisola, O., Hormozi, S. & Brandt, L. 2018 Turbulent channel flow of an elastoviscoplastic fluid. J. Fluid Mech. 853, 488–514.CrossRef Google Scholar

Saramito, P. 2007 A new constitutive equation for elastoviscoplastic fluid flows. J. Non-Newtonian Fluid Mech. 145 (1), 1–14.CrossRef Google Scholar

Saramito, P. 2009 A new elastoviscoplastic model based on the Herschel–Bulkley viscoplastic model. J. Non-Newtonian Fluid Mech. 158 (1-3), 154–161.CrossRef Google Scholar

Shahmardi, A., Zade, S., Ardekani, M. N., Poole, R. J., Lundell, F., Rosti, M. E. & Brandt, L. 2019 Turbulent duct flow with polymers. J. Fluid Mech. 859, 1057–1083.CrossRef Google Scholar

Shu, C.-W. 2009 High order weighted essentially nonoscillatory schemes for convection dominated problems. SIAM Rev. 51 (1), 82–126.Google Scholar

Sikorski, D., Tabuteau, H. & de Bruyn, J. R. 2009 Motion and shape of bubbles rising through a yield-stress fluid. J. Non-Newtonian Fluid Mech. 159 (1-3), 10–16.CrossRef Google Scholar

de Souza Mendes, P. R. 2007 Dimensionless non-Newtonian fluid mechanics. J. Non-Newtonian Fluid Mech. 147 (1), 109–116.CrossRef Google Scholar

Sverdrup, K., Almgren, A. & Nikiforakis, N. 2019 An embedded boundary approach for efficient simulations of viscoplastic fluids in three dimensions. Phys. Fluids 31 (9), 093102.CrossRef Google Scholar

Tabuteau, H., Coussot, P. & de Bruyn, J. R. 2007 Drag force on a sphere in steady motion through a yield-stress fluid. J. Rheol. 51 (1), 125–137.CrossRef Google Scholar

Tokpavi, D. L., Jay, P., Magnin, A. & Jossic, L. 2009 Experimental study of the very slow flow of a yield stress fluid around a circular cylinder. J. Non-Newtonian Fluid Mech. 164 (1-3), 35–44.CrossRef Google Scholar

Tokpavi, D. L., Magnin, A. & Jay, P. 2008 Very slow flow of Bingham viscoplastic fluid around a circular cylinder. J. Non-Newtonian Fluid Mech. 154 (1), 65–76.CrossRef Google Scholar

Uhlmann, M.2003 First experiments with the simulation of particulate flows. Tech. Rep. Centro de Investigaciones Energeticas.Google Scholar

Uhlmann, M. 2005 An immersed boundary method with direct forcing for the simulation of particulate flows. J. Comput. Phys. 209 (2), 448–476.CrossRef Google Scholar

Vázquez-Quesada, A., Bian, X. & Ellero, M. 2016 Three-dimensional simulations of dilute and concentrated suspensions using smoothed particle hydrodynamics. Comput. Part. Mech. 3 (2), 167–178.CrossRef Google Scholar

Vishnampet, R. & Saintillan, D. 2012 Concentration instability of sedimenting spheres in a second-order fluid. Phys. Fluids 24 (7), 073302.CrossRef Google Scholar

Volarovich, M. P. & Gutkin, A. M. 1953 Theory of flow of a viscoplastic medium. Colloid J. 15, 153.Google Scholar

Wachs, A. & Frigaard, I. A. 2016 Particle settling in yield stress fluids: limiting time, distance and applications. J. Non-Newtonian Fluid Mech. 238, 189–204.CrossRef Google Scholar

Yoshioka, N., Adachi, K. & Ishimura, H. 1971 On creeping flow of a viscoplastic fluid past a sphere. Kagaku Kogaku 10 (1144), 6–31.Google Scholar

Yu, Z. & Wachs, A. 2007 A fictitious domain method for dynamic simulation of particle sedimentation in Bingham fluids. J. Non-Newtonian Fluid Mech. 145 (2-3), 78–91.CrossRef Google Scholar

Zheng, R., Phan-Thien, N. & Tanner, R. I. 1991 The flow past a sphere in a cylindrical tube: effects of intertia, shear-thinning and elasticity. Rheol. Acta 30 (6), 499–510.CrossRef Google Scholar

Zisis, T. & Mitsoulis, E. 2002 Viscoplastic flow around a cylinder kept between parallel plates. J. Non-Newtonian Fluid Mech. 105 (1), 1–20.CrossRef Google Scholar

Figure 1.Computational domain and boundary conditions.

Figure 2.Mechanical model schematic of an EVP material proposed by Saramito (2009).

Table 1.Computational matrix.

Figure 3.Time evolution of shear stress in a stationary simple Couette flow. The solid lines and symbols represent the analytical solution given by Saramito (2009) and the present numerical simulations, respectively.

Figure 4.(a) Normalized velocity magnitude around the particle settling in the $y$–$z$ centreplane ($x=3D$) through Bingham fluid and (b) EVP material at $Wi_{\infty }=0.1$. (c) Streamwise velocity around a particle settling in a Bingham fluid and an EVP material at the channel centre line ($x=3D$, $z=6D$). In all figures, the Bingham number $Bi=0.5$.

Table 2.Dimensionless numbers of the pure sedimentation simulations with associated gravity number $Y_{g}$, smaller than its critical value $Y_{g}^{c}$.

Figure 5.Velocity vectors around the sphere settling in an EVP material in the central $y$–$z$ plane. The dashed red and solid green boxes magnify the recirculation zones and the flow stagnation points. The dimensionless numbers are the same as figure 4.

Figure 6.Stress fields around the particle settling in an EVP fluid in the $y$–$z$ plane at $x=3D$; (a) the $\unicode[STIX]{x1D70F}_{yz}$ shear stress and (b) the $\unicode[STIX]{x1D70F}_{yy}$ normal stress. The angle $\unicode[STIX]{x1D701}$ is measured from the north pole of the sphere. (c) Stress relaxation versus the distance from the particle, $r$, at an angle $\unicode[STIX]{x1D701}=30^{\circ }$ measured from the north pole of the sphere. Here, $r=0$ corresponds to the particle surface. The absolute value of the stresses are shown for the same case as in figure 4.

Figure 7.Surfaces of unyielded regions around the particle in an EVP material at $Wi_{\infty }=0.1$ in the absence of cross-shear flow for various $Bi$: (a) $Bi=0.05$, (b) $Bi=0.1$, (c) $Bi=0.13$, (d) $Bi=0.5$, (e) $Bi=1$. Red represents the unyielded zone, while grey shows the yield surface boundary.

Figure 8.Evolution of the unyielded zones for flow of an EVP material around a particle for various $Bi$ numbers at $Wi_{\infty }=0.1$; (a,f) $Bi=0.05$, (b,g) $Bi=0.1$, (c,h) $Bi=0.13$, (d,i) $Bi=0.5$, (e,j) $Bi=1$. The first (a–e) and second rows (f–j) represent the unyielded zones in the central $y$–$z$ and $x$–$y$ planes, respectively. Blue and grey represent the yielded and unyielded regions.

Figure 9.Colourmaps of the second invariant of the shear rate for the flow of an EVP material around a particle for various $Bi$ numbers at $Wi_{\infty }=0.1$; (a,f) $Bi=0.05$,(b,g) $Bi=0.1$, (c,h) $Bi=0.13$, (d,i) $Bi=0.5$, (e,j) $Bi=1$. The first (a–e) and second rows (f–j) represent the shear rate in the central $y$–$z$ and $x$–$y$ planes, respectively.

Figure 10.Total drag coefficient and its individual components for a sphere settling in an EVP fluid as a function of the $Bi$ number at $Wi_{\infty }=0.1$.

Table 3.Comparison the total drag coefficient and its components on the surface of the sphere settling in a viscoplastic and in an EVP fluid at $Bi=0.05$. The settling Weissenberg number $Wi_{\infty }=0.1$ for EVP material.

Figure 11.Contour plots of the dynamic pressure and the three components of the viscous stress distribution on the surface of the particle settling in an EVP fluid at different $Bi$ numbers and $Wi_{\infty }=0.1$. First, second, third, fourth and fifth rows show the distributions at $Bi=0.05,0.1,0.13,0.5$ and $1$, respectively. Axes information in all panels can be extrapolated from panel (a).

Figure 12.Total drag coefficient of a particle settling in a Bingham fluid from Blackery & Mitsoulis (1997) and that of a particle settling in an EVP material (with elasticity $Wi_{\infty }=0.1$) as well as the elastic drag resulting from (3.8).

Figure 13.Contour plots of the dynamic pressure and of three components of viscous stress distribution on the surface of the particle settling in a Bingham fluid (a–d) at$Bi=0.05$ and in the EVP fluid (e–h) at $Bi=0.05$, $Wi_{\infty }=0.1$. Axes information in all panels can be extrapolated from panel (a).

Figure 14.Pressure and viscous stress distributions on the surface of the sphere settling in Bingham and EVP fluids for the same cases depicted in figure 13. (a) Spherical coordinate system, (b) pressure, (c) $\unicode[STIX]{x1D70F}_{xy}$, (d) $\unicode[STIX]{x1D70F}_{yy}$, (e) $\unicode[STIX]{x1D70F}_{yz}$. The pressure, $\unicode[STIX]{x1D70F}_{xy}$ and $\unicode[STIX]{x1D70F}_{yy}$, distributions are plotted for $\unicode[STIX]{x1D713}=0$, while the $\unicode[STIX]{x1D70F}_{yz}$ is plotted at $\unicode[STIX]{x1D719}=0$.

Figure 15.(a) Normalized velocity colourmaps in the $y$–$z$ centreplane ($x=3D$) for the case of pure sedimentation of a single sphere in an EVP material. (b) The same as panel (a), but for the case of shear-induced sedimentation. (c) Streamwise fluid velocity in the flow direction at the channel centreline ($x=3D$, $z=6D$) for the cases of pure sedimentation (orange markers) and shear-induced sedimentation (blue markers) in an EVP material. In all figures the dimensionless numbers are $Bi=0.5$, $Wi_{\infty }=0.1$.

Figure 16.Velocity vectors around the sphere settling through an EVP material ($Bi=0.5$ and $Wi_{\infty }=0.1$) in the centre $y$–$z$ plane for the case of (a) pure sedimentation, (b) shear-induced sedimentation. The blue boxes magnify the flow upstream of the sphere.

Figure 17.Drag coefficient for a settling sphere in an EVP fluid ($Wi_{\infty }=0.1$, $Wi=0.01$) for the case of pure sedimentation (blue circular markers) and shear-induced sedimentation (red triangular markers) as a function of the Bingham number. The green square markers represent the drag coefficient deduced from the pure sedimentation simulations data by taking into account the second invariant of the shear rate tensor.

Figure 18.Total and individual drag components on the sphere settling in an EVP fluid at $Wi_{\infty }=0.1$ in the absence of cross-shear flow (dark markers) and in the presence of cross-shear flow (light markers) as a function of $Bi$ number.

Figure 19.Contour plots of the pressure and viscous stress components on the surface of the sphere settling in an EVP fluid at $Bi=1$, $Wi_{\infty }=0.1$. The first row corresponds to pure sedimentation ($Wi=0$) and the second row to shear-induced sedimentation ($Wi=0.01$). Axes information in all panels can be extrapolated from panel (a).

Table 4.Comparison of the drag coefficient and its individual components on a sphere settling in the absence and the presence of cross-shear flow in a Bingham fluid at $Bi=0.05$.

Figure 20.Contour plots of the dynamic pressure and of the three components of the viscous stress on the surface of the particle in the case of pure sedimentation (a–d) and shear-induced sedimentation (e–h) in a Bingham fluid at $Bi=0.05$. Axes information in all panels can be extrapolated from panel (a).

Table 5.Comparison of the drag coefficient and its individual components on a sphere settling in sheared viscoplastic (Bingham) fluid at $Bi=0.05$ and EVP material at $Wi_{\infty }=0.1$, $Bi=0.05$.

Figure 21.Drag coefficient versus the Bingham number for a sphere settling in sheared EVP fluid at different settling Weissenberg numbers, $Wi_{\infty }$.

Figure 22.Contour plots of pressure and viscous stress components projected on the surface of the sphere settling in sheared EVP fluid at $Bi=1$. The first row corresponds to $Wi_{\infty }=0.1$ ($Wi=0.01$) and the second row to $Wi_{\infty }=1$ ($Wi=0.1$). Axes information in all panels can be extrapolated from panel (a).

Table 6.Comparison of the drag coefficient and its individual components on a sphere settling in a sheared EVP fluid at $Bi=1$ for two values of $Wi_{\infty }$.

Figure 23.Comparison of the yielded/unyielded zones around a particle settling in a Bingham fluid at the confinement ratio of $4$ and at $Bi=0.108$; (a) computations of Blackery & Mitsoulis (1997) at $Re_{p}=0$, (b) present work at $Re_{p}=1$. White and black colours represent the yielded and unyielded zones, respectively.

Table 7.Comparison of Stokes drag coefficient of a sphere settling in a Bingham fluid resulting from the present work (second column) and the computations of Blackery & Mitsoulis (1997).

Figure 24.Solid volume fractions (highlighted area) for grid cells around $p(i,j)$,$u(i+9/2,j+2)$ and $v(i+2,j+1/2)$. The sphere boundary is shown by the red dashed line.

Table 8.Comparison of drag coefficient and drag components from the present work, PW, with Bohlin’s approximation and theoretical solutions for Newtonian fluids.

Article contents

Numerical simulations of a sphere settling in simple shear flows of yield stress fluids

Abstract

Information

1 Introduction

2 Problem definition

2.1 Mechanical model of EVP material

2.2 Mathematical model

2.3 Boundary conditions

2.4 Computational matrix

2.5 Numerical method

2.6 Code validation

2.6.1 EVP Couette flow

3 Results and discussion

3.1 Pure sedimentation in EVP fluids

3.1.1 Velocity and stress fields

3.1.2 Yielded/unyielded regions

3.1.3 Drag coefficients

3.2 Shear-induced sedimentation

3.2.1 Velocity field

3.2.2 Yielded/unyielded regions

3.2.3 Drag coefficients

3.3 Shear-induced sedimentation in a Bingham fluid

3.4 Nonlinear coupling

3.5 Shear-induced sedimentation at higher elasticity

4 Conclusion

Acknowledgements

Declaration of interests

Appendix A. Bingham fluid modelling

A.1 Numerical Method

A.2 Validation case

Appendix B. Numerical integration procedure

B.1 Validation case

References

Save article to Kindle

Save article to Dropbox

Save article to Google Drive

Reply to: Submit a response

Your details

You have entered the maximum number of contributors

Conflicting interests