Finite-size spherical particles in a square duct flow of an elastoviscoplastic fluid: an experimental study

The present experimental study addresses the flow of a Yield Stress Fluid with some elasticity (Carbopol gel) in a square duct. The behaviour of two fluids with lower and higher yield stress is investigated at multiple Reynolds numbers $Re^* \in$ (1, 200) and Bingham numbers $Bi \in$ (0.01, 0.35). A secondary flow consisting of eight vortices is observed to recirculate the fluid from the corners to the core. Their extent and intensity grows with increasing $Re^*$. The second objective of this study is to explore the change in flow in the presence of particles. Almost neutrally-buoyant finite-size spherical particles with duct height, $2H$, to particle diameter, $d_p$, ratio of 12 are used at two volume fractions $\phi$ = 5 and 10\%. Particle Tracking Velocimetry (PTV) is used to measure the velocity of these refractive-index-matched spheres, and PIV to extract the fluid velocity. Simple shadowgraphy is also used for qualitatively visualising the development of the particle distribution along the streamwise direction. The particle distribution pattern changes from being concentrated at the four corners, at low flow rates, to being focussed along a diffused ring between the center and the corners, at high flow rates. The presence of particles induces streamwise and wall-normal velocity fluctuations in the fluid phase; however, the primary Reynolds shear stress is still very small compared to turbulent flows. The size of the plug in the particle-laden cases appears to be smaller than the corresponding single phase cases. Similar to Newtonian fluids, the friction factor increases due to the presence of particles, almost independently of the suspending fluid matrix. Interestingly, predictions based on an increased effective suspension viscosity agrees quite well with the experimental friction factor for the concentrations used in this study.


YSF in a square duct
Square duct flows of YSF, mostly of the Bingham type, have been a subject of several numerical studies. Following Taylor & Wilson (1997), Van Pham & Mitsoulis (1998) performed simulations of laminar duct flow using the Bingham model and proposed a master curve relating the flow rate to the pressure drop; useful information in the design of extrusion geometries. Saramito & Roquet (2001) used the augmented Lagrangian method for a Bingham model fluid flow in a square duct at varying Bi to compute the critical value of Bi above which the flow stops, known as the stopping criterion (also see Mosolov & Miasnikov (1965)). Huilgol & You (2005) extended the above method to a HB fluid model and found that, for a fixed Bi (or Od in their case), the plug region increases as the power law exponent in the HB model decreases.
A YSF may also display elastic properties both below and above yield stress, thus behaving like a viscoelastic solid before yielding and a viscoelastic liquid after yielding. Recently, Fraggedakis, Dimakopoulos & Tsamopoulos (2016b) compared the prediction from different elastoviscoplastic constitutive models in simple rheometric flows. Letelier, Siginer & González (2017) introduced elasticity in their simulations of a Bingham fluid for an approximate square duct geometry and observed that, for a fixed Bi, the flow rate increases with increasing elasticity (quantified using a Weissenberg number Wi, which is the ratio of elastic stresses in the form of a normal-stress difference to the viscous stress due to shear forces). In a later work by the same authors, by including higher-order terms in Wi, a secondary flow in the form of streamwise vortices is seen (Letelier et al. 2018). Increasing the Bi at the same Wi reduced the intensity of the secondary flow, and displaced the vortices further away from the centre of the duct. Similar observations are also reported in this work, the novelty being the HB nature of the experimental fluid against the Bingham numerical model along with a higher Bi in our case compared to their simulations, as will be described in the following sections. Often a Deborah number De is defined as the ratio of the relaxation time of the fluid λ to the characteristic observation time scale of the flow T f . Under certain circumstances, De = λ/T f is equivalent to the Wi (Poole 2012b), which is more appropriately related to the characteristic rate of deformation (Dealy 2010). Since both Wi and Re * increase with the flow rate for a given elastoviscoplastic fluid, an elasticity number El is defined as Wi/Re * to quantify the effects of elastic forces compared to inertial forces.
Many researchers since Ericksen (1956) have investigated the existence and strength of the aforementioned streamwise vortices in laminar rectangular duct flows. This secondary flow originates from viscoelastic forces in non-circular pipes, specifically the second normal-stress difference (Dodson et al. 1974), and its strength increases with the elasticity (Debbaut & Dooley 1999). Thus, no secondary flow would be observed in a purely viscous or viscoplastic fluid or for a certain special relationship between the second normal-stress difference and the fluid viscosity (Xue, Phan-Thien & Tanner 1995). Despite being very weak (around two orders in magnitude lower than the mean axial velocity), the presence of such a secondary velocity field may have substantial effects on the rate of heat transfer (Kostic 1994;Gao & Hartnett 1996). Its weak nature also presents challenges in measuring it experimentally (Schroeder & Jeffrey 2011). Yue, Dooley & Feng (2008) provided a general criterion to predict the direction of these secondary flows.
1.4. Studies in particle-laden flow Dispersed phases like particles, drops or bubbles may appear as desired or undesired components in the final product made using YSF. Hence, their interaction with YSF is of practical importance and this paper partly addresses this problem.
In the case of particle-laden flows, sedimentation of a single spherical particle in a YSF is the most common study on account of its obvious importance in gaining fundamental understanding and as a precursor for multi-particle problems. A heavy isolated particle can settle inside a YSF only when the buoyancy force exceeds the resistance offered by the yield stress. Accordingly, a dimensionless yield-gravity parameter Y = τ y /((ρ p − ρ f )gd p ) can be defined such that the particle moves only when Y is lower than a threshold critical value ≈0.2 while noting that the exact value is still a topic of investigation (Chhabra 2006). Creeping flow of an ideal yield stress fluid (Bingham or Herschel-Bulkley) around a sphere is reversible and displays fore-aft symmetry (Beaulne & Mitsoulis 1997;Putz & Frigaard 2010).

Need for careful preparation of the solution
In contrast to the above observation, when Putz et al. (2008) performed particle image velocimetry (PIV) measurements of the flow field around a spherical particle sedimenting in a Carbopol solution at low Re < 1, they observed a breaking of the symmetry that is not observed in Newtonian fluids. In particular, the characteristics of the flow regions around the falling sphere were associated with the rheological properties of the fluid, thus calling for numerical models to include both elastic (Fraggedakis, Dimakopoulos & Tsamopoulos 2016a) and hysteresis/thixotropic effects to accurately reproduce experimental observations as well as for experiments to ensure careful preparation of the fluid. On the other hand, Tabuteau, Coussot & de Bruyn (2007) could obtain reproducible results in experiments for the drag on a settling sphere, in agreement with the simulations of Beaulne & Mitsoulis (1997). The authors attributed this to the careful preparation (homogenising for ten days!) of the yield stress fluid.
Indeed, in the recent work by Dinkgreve et al. (2018), it was clearly shown that the discrepancy in the simple yield stress behaviour of Carbopol observed in previous studies was most likely due to lack of optimal mixing. At low shear rates, transient shear banding (Divoux et al. 2010) may occur that can cause apparent hysteresis in the flow curve i.e. the relation between shear stress and shear rate in simple shear, due to a lack of sufficient measurement time. The measurement time required to reach a steady state approaches unrealistically large values in the presence of wall slip, which occurs when a smooth rotating surface is used during rheometry (Poumaere et al. 2014).
1.6. Particle migration in non-Newtonian fluids Since the Carbopol gel used in this study exhibits small but measurable elastic effects both below and above the yield stress, the normal stresses and ensuing secondary flows are expected to affect the particle migration and their equilibrium distribution. Hence, a review of particle migration inside a duct flow of viscoelastic fluid is appropriate.
Cross-stream particle migration in viscoelastic fluids is quite different to their migration in a Newtonian fluid (D'Avino & Maffettone 2015). In the absence of inertia, which is typical of microfluidic applications (e.g. cell focussing and separation), depending on the initial particle position, fluid elasticity drives the particle towards the channel centreline or the closest wall and from there towards the nearest corner in case of a duct geometry (Villone et al. 2013). This can be understood by observing the distribution of the first normal-stress difference, which has local minima at the centre and the corners of the duct cross-section (Yang et al. 2011). The shear-thinning effects reduce the first normal-stress difference (Li, McKinley & Ardekani 2015) thus, augmenting particle migration towards the closest wall (D'Avino et al. 2012). When second normal-stress differences are present (e.g. in a duct), the particle migration velocity, proportional to De(d p /(2H)) 2 , may be higher or lower than the ensuing secondary flow velocity, proportional to De 4 (Yue et al. 2008). Note that the observation time T f in the definition of De above is taken as the inverse of the nominal shear rate: T f = 1/(U Bulk /H) by Villone et al. (2013). Thus, at larger De and smaller particle confinement d p /(2H), the secondary flow may overcome the migration velocity and the particle may also find a stable position near the centre of the streamwise vortices that are characteristic of the secondary flow (Villone et al. 2013).
It is shown by Seo, Kang & Lee (2014) that complex interplay between elasticity, inertia, shear thinning and confinement can lead to multiple stable and unstable equilibrium configurations. Lim et al. (2014a) realised that particles in microfluidic flows can be focussed at a high throughput rate even at very high Re ≈ 10 000 i.e. high inertia, provided that the elastic normal stresses also increase proportionally i.e. the important criterion governing particle focussing is the elasticity number El.
In the Stokes regime i.e. negligible particle inertia, flow of a particle suspension in a YSF inside a tube has been recently simulated by Siqueira & de Souza Mendes (2019) using a regularised viscosity function, with the diffusive flux model (Phillips et al. 1992) being used to describe the shear-induced particle migration. Particles were found to concentrate at the boundary of the plug and with increasing particle concentration, the maxima shifted towards the tube wall. No particles were found inside the plug due to the high but finite viscosity gradient (due to the regularisation) that slowly pushed them radially outwards. A similar method was used to model the particle migration for a tube flow of a Bingham fluid in Lavrenteva & Nir (2016). The maximum particle concentration was found at the interface between the yielded and unyielded regions. On the other hand, in the simulations of Hormozi & Frigaard (2017) to model particleladen flow in a fracture, particles were found to concentrate in the plug. 1.7. Importance of the present study Most of the above experimental studies concerning particle migration are microfluidic in their treatment i.e. they have a very low fluid and particle inertia. Also, they deal with a very low particle concentration (Φ 0.1 %) i.e. negligible inter-particle interactions, that is typical of flow-focusing experiments, and to our knowledge there have not been many experimental studies devoted to multi-particle dynamics with the exception of the impressive observations made by Gauthier, Goldsmith & Mason (1971) and Tehrani (1996) a while ago. A noteworthy contribution of this study lies in extending previous scientific investigations beyond the above regimes using a novel set-up with a large flow cross-section where the flow can be resolved at the particle scale, a feature that is very difficult to capture in microfluidic devices. In addition, the presence of yield stress in the suspending fluid and its interaction with the finite-size dispersed phase is expected to complement studies using the continuum-based approaches e.g. Hormozi & Frigaard (2017).

Outline
In the following sections, we start by describing the experimental set-up and the PIV-PTV (particle tracking velocimetry) measurement techniques. The fluid rheology and particle properties are also mentioned. Later, in the results section, we first describe the single-phase flow of two YSFs, with high (HYS) and low (LYS) yield stresses, at varying Re * and Bi. This is then followed by the results concerning the particle-laden cases. Finally, in the discussion section, we try to interpret the above results by supporting our main observations with simple visual images obtained from shadowgraphy experiments, which are described therein. This is followed by a summary of the main conclusions.

Experimental set-up
The experimental set-up used in this study has been previously used in Zade et al. (2018) to investigate particle-laden turbulent flows of Newtonian and non-Newtonian flows up to φ = 20 %, and is shown in figure 1. Only the most important features are described here and the interested reader is referred to the previous publications for additional details. The flow loop consists of a 5 m long transparent square duct with 50 mm × 50 mm cross-section (refer to figure 1a). The particles are introduced in the conical tank and are circulated through the loop for a long enough time before measurements begin, so that they homogeneously disperse in the carrier fluid. A static mixer (Vortab Company, CA, USA) is installed close to the inlet of the duct to neutralise any swirling motions that may arise from the long 90 • bend at the exit of the tank. The temperature is maintained at nearly 20 • C by means of an immersed-coil heat exchanger in the tank. A very gentle disc pump (Model: 2015-8-2HHD Close coupled, Discflo Corporations, CA, USA) has been chosen, similar to Draad, Kuiken & Nieuwstadt (1998), to minimise the mechanical breakage of the particles and avoid unwanted pulsations in the flow. An electromagnetic flowmeter (Krohne Optiflux 1000 with IFC 300 signal converter, Krohne Messtechnik GmbH, Germany) is used to measure the volume flow rate of the particle-fluid mixture. The pressure loss is measured over a distance of 54H starting at 140H from the inlet using a differential pressure transducer (0-1 kPa, Model: FKC11, Fuji Electric France, S.A.S.). Data acquisition from the camera, flow meter and pressure transducer is performed using a National Instruments NI-6215 DAQ card using Labview TM software. The entry flow problem for YSF has been studied by Ookawara et al. (2000), amongst others, who proposed that the spatial development of the streamwise velocity at the edge of the plug is a suitable indicator for a fully developed flow. The development length was deduced in terms of a modified Reynolds number accounting for the plug radius, and for the maximum Re * encountered in our case, this development length would be less than 30H. Thus, at our measurement location ≈150H, the flow can be considered to be fully developed for the single phase cases. For the particle-laden cases, the evolution of the particle concentration is qualitatively observed using shadowgraphy techniques, described later in the discussion section, and it is observed that an equilibrium concentration is established upstream of the measurement location.

Fluid rheology
An aqueous solution of Carbomer powder supplied as Carbopol NF 980 (Lubrizol Corporation, USA), a commercial thickener, that is neutralised with sodium hydroxide (NaOH), is used as the model yield stress fluid. It is chosen due to its high transparency, small thixotropy and material stability (very slow ageing).
For each of our Carbopol batches, a weighed amount of Carbopol powder (0.25 % by weight) was first dispersed in 24 kg of water at room temperature (∼20 • C). A high shear mixer (Silverson AX5, Silverson Machines, Inc., USA) was used to disperse the powder at a maximum of ∼1200 rpm for up to ∼30 minutes. The dispersion was then left stationary for ∼30 minutes to allow for air bubbles to evacuate. Then, a 18 wt./v % of NaOH solution at 2.3 times the mass of Carbomer powder was used to neutralise the dispersion. The NaOH was gradually added with a pipette whilst stirring the solution gently at low rotational speeds (∼135 rpm) using  . At the end of the neutralisation process, the pH was noted to be ideally within the required range (6.5-7.0) so as to ensure maximum swelling and, hence, a high yield stress. The neutralised solution was further mixed in a large cement mixer (1402 HR, AL-KO, Germany) at a low rotational speed for an ∼5 additional hours to ensure complete homogenisation.
The whole preparation process led to a concentrated solution with a high yield stress. Experiments had to be performed with a thinner (less viscous) fluid in order to facilitate pumping. Thus, the concentrated solution was gradually diluted in the flow loop by mixing with water and recirculating for ∼2 hours before measurements began. Since the shear rates in the flow loop are not extremely high, solution degradation and the resulting change in properties can be considered insignificant. Repeatability of the properties of the final solution was ascertained by (i) measuring the flow curves in a rheometer and (ii) by monitoring the pressure drop at a fixed flow rate. It was found that by systematically following the above mixing protocol, solutions with nearly the same pressure drop exhibiting the same flow curve were obtained repeatedly. In order to minimise the occurrence of bubbles arising due to dissolved gases in the water, tap water was heated so as to remove these dissolved gases and then cooled before being used for preparation of the Carbopol solution. Air bubbles entrained during introduction of the hydrogel particles were removed in the conical tank that is open to the atmosphere.
The rheological tests focused on the determination of the steady-state flow curve. They were performed using a TA AR-2000ex rheometer (TA Instruments, Inc., USA) with a vane and cup geometry. A pre-shear of 300 s −1 was applied for 10 s followed by a rest period of 40 s to establish a reproducible state. Flow sweeps in both ascending and descending controlled shear rate were then carried out in a range of shear rates (0.001 to 80 1 s −1 with 10 points per decade) at room temperature to determine the flow curves and to check for thixotropy. This range of shear rates includes the maximum shear rate in the experiments. For each measurement point, the shear rate was held constant for 30 s. Figure  2017) were visible for a few points at low shear rates during the up sweep i.e. the accumulated strain was perhaps insufficient for the material to flow. Hence, the down sweep curve is used to fit the Herschel-Bulkley model, whose parameters are mentioned in the legend of figure 2. The repeatability of the flow curves was ensured by performing the above tests on three separate solutions prepared independently of each other, both before and after adding the particles, and an average curve is used. The working fluids display a shear-thinning behaviour after yielding i.e. they are yield pseudo-plastic fluids. The dynamic yield stress determined above (static yield stress is determined using creep tests) is low compared to other wall-bounded flow studies of YSF (see Güzel et al. (2009) and Escudier & Presti (1996) amongst others) but for such flows, the non-dimensional Bi is the relevant criterion to assess the importance of plasticity on the flow. This point about the relevance of Bi will be illustrated further by means of velocity profiles that exhibit a solid plug. The viscoelastic behaviour near the yield limit is measured by oscillatory tests using a splined bob and cub attachment in a Kinexus pro+ rheometer (Malvern Panalytical, UK). The linear viscoelastic region that exists at small strain γ amplitudes can be identified by nearly constant elastic G and viscous G moduli over a wide range of strain amplitudes. This is shown in figure 3 where oscillatory measurements are performed at a frequency ω of 1 Hz. The point of cross-over between G and G for each of the two fluids at 1 Hz, Hτ y ∼ 4 Pa and Lτ y ∼ 1 Pa, is of a similar order of magnitude as the yield stress, albeit slightly larger. Figure 4 shows the variation of G and G with the frequency ω at multiple strain amplitudes γ , starting from the linear and extending towards the nonlinear viscoelastic regime associated with yielding. In the linear regime i.e. at lower γ , the elastic modulus G (refer to figure 4a) is approximately flat and around one order of magnitude greater than the viscous modulus G (refer to figure 4b); G and G are nearly independent of γ in this linear regime (γ 0.5 %), as expected. As mentioned in Varges et al. (2019), at these low amplitudes the microgels deform elastically but they do not move significantly relative to each other and hence have a low internal dissipation (also see Piau 2007). As the deformation amplitude increases, the microstructure starts breaking, causing elastohydrodynamic friction forces to rise and dissipation to increase, leading to a higher loss modulus. Finally, at higher γ , outside the The viscoelastic moduli G and G are relevant at small shear rates. In order to quantify the viscoelastic effects at higher shear rates, other steady-state viscometric functions namely, the first and second normal-stress differences: N 1 and N 2 needs to be measured. These measurements are fraught with difficulties in YSF and for N 1 , both positive (Tehrani 1996;Peixinho et al. 2005;Piau 2007) and negative values (Janmey et al. 2007) are reported. In the recent work of De Cagny et al. (2019), the authors investigated common simple i.e. non-thixotropic YSF, including Carbopol, and found that N 1 > 0 and N 2 < 0 with N 1 > N 2 ; they increase quadratically with the shear stress, both below and above the yield stress. The above authors proposed that, in the absence of normal-stress measurements, N 1 can be estimated by using the linear elastic modulus G as N 1 ∼ τ 2 /G . Here, τ is the shear stress and, for our fluids, this implies that N 1 < 6 Pa (for LYS at the maximum shear rate of figure 2). Such small values are quite difficult to measure, especially due to normal force drift and inertia (Poole 2012a). It was hence not possible to accurately measure the normal-stress differences for our weakly elastic fluid. Nevertheless, viscoelasticity of Carbopol after yielding has been observed in many studies (e.g. see Taylor & Bagley 1974;Laurencena & Williams 1974). Peixinho et al. (2005) also observed substantial drag reduction in turbulent flow of Carbopol, a typical viscoelastic effect.
2.3. Particle properties As already noted in our previous work (Zade et al. 2018) the particles are commercially procured super-absorbent (polyacrylamide-based) hydrogel spheres which are delivered in dry form. After grading them into different sizes using sieves, one fairly mono-dispersed fraction is used in these experiments. These dry particles are mixed with tap water that has a very small amount of fluorescent rhodamine dye (a few ppm) dissolved in it, and then left submerged for around one day. Ultimately, they grow to an equilibrium size of 4.2 ± 0.8 mm (3 times standard deviation), thus yielding a duct height to particle diameter ratio 2H/d p of ∼12. The tiny amount of rhodamine is absorbed by each particle, and makes them glow in the laser sheet, so that they can be detected in the PIV images, as will be shown later. The particle size was determined both by a digital imaging system and from the PIV images of particles in the flow. A small Gaussian-like spread in the particle diameter was observed. The fact that a Gaussian-like particle size distribution has a small effect on the flow statistics has already been shown in Fornari, Picano & Brandt (2018).
The density ratio of the particle to tap water ρ p /ρ f is found to be equal to 1.0035 ± 0.0003. Since the final concentration of Carbomer in the working fluid is less than 0.1 %, the change in density of the fluid from tap water is assumed to be insignificant. Since the typical flow velocities are relatively small (U Bulk < 0.3 m s −1 ) in the present flow configuration, the corresponding dynamical forces are low under these conditions and the hydrogel particles did not exhibit any visible deformation (as also confirmed from the time-resolved movies of the particles in flow). Also, as pointed out by Gondret, Lance & Petit (2002), for the low impact Stokes number, collisions between particles or between a particle and the wall are dominated by viscous effects, and particles do not rebound i.e. the effective coefficient of restitution will be very small. Based on the above arguments, the particles can be considered to behave as rigid spheres. For details about the density measurements and calculations pertaining to the restitution coefficient, the reader is referred to our previous work (Zade et al. 2019a).

Velocity measurement technique
The coordinate system is sketched in figure 1(a) with x in the streamwise, y in the wall-normal and z in the spanwise directions. The velocity field is measured using two-dimensional PIV (2D-PIV) in multiple spanwise planes (z/H = 0 to 1) to measure the inhomogeneous velocity field in the square duct. The PIV set-up consists of a continuous wave laser (wavelength = 532 nm, power = 2 W) and a high-speed camera (Phantom Miro 120, Vision Research, NJ, USA), as shown in figure 1(b). The thickness of the laser light sheet is 1 mm. The seeding particles are hollow glass microspheres (Sphericel 110P8, Potters industries Inc., d p = 9-13 µm, ρ p = 1.10 ± 0.05 gm cm 3 ) and polyamide microspheres (d p ≈ 55 µm, ρ p = 1.2 gm cm −3 ). The following paragraphs are briefer reproductions of the PIV-PTV algorithm already described in Zade et al. (2019a).
PIV image pairs are captured at a resolution of approximately 60 mm/1024 pixels; the size of the measurement section is 50 mm × 50 mm. For each flow rate, a time delay was chosen so that the maximum pixel displacement, based on the mean velocity, did not exceed a quarter of the size of the final interrogation window (IW) (Raffel et al. 2013). Thus, it can range from 36 × 10 −3 s for the lowest flow rate = 2 l min −1 (litres per minute), to approximately 1.8 × 10 −3 s for the highest flow rate = 40 l min −1 . The exposure time is kept small enough (<0.1 × 10 −3 s) so that enough light is captured by the camera sensor while pixel displacement within one frame is significantly less than 1 pixel. Images were processed using an in-house, three-step, fast Fourier transform (FFT)-based, cross-correlation algorithm (Kawata & Obi 2014). The degree of overlap is approximately 47 % and can be estimated from the fact that the corresponding final resolution is 1 mm × 1 mm per IW. Between 200 and 400 image pairs have been observed to be sufficient to ensure statistically converged results. Figure 5(a) depicts one image from a typical PIV sequence for particle-laden flow. Raw images captured during the experiment were saved in groups of two different intensity levels. The first group of images, a typical example being figure 5(a), were used for regular PIV processing according to the algorithm mentioned above to find the fluid velocity field. The same first group of images were later contrast enhanced FIGURE 5. Images involved to calculate the fluid and particle velocity fields. The raw image (a) is used for PIV analysis. The enhanced image (b) is used to detect circles and for the subsequent PTV analysis. The combined PIV (for fluid phase denoted by green arrows) and PTV (for particle phase denoted by blue arrows) velocity field is shown in (c). The detected particles are also shown in (b). The above images correspond to LYS fluid, φ = 10 %, z/H = 0.9 and Re * = 156.
in the post-processing step and constituted the second group of images, an example being figure 5(b). These second group of images were used to detect the particles. Again, the details of the image post-processing and the PTV algorithm are explained at length in Zade et al. (2019a). The velocity of both the fluid and particle phases is defined on a spatially fixed Eulerian grid. Consequently, a mask matrix is defined, which assumes the value 1 if the point lies inside the particle and 0 if it lies outside the particle. The particle velocity is determined using PTV at its centre, which is assigned to the grid points inside the particle (mask = 1). The velocity field of the particle phase is now available at the same grid points as that of the fluid and the ensemble averaging, reported later, are phase averaged statistics. Figure 5(c) shows the combined fluid (PIV) and particle (PTV) velocity field. It may be observed that the intensity of the laser sheet is weaker away from the centre of the images and particles are not accurately detected in this region. Hence, PIV and PTV results are extracted from a reduced area neglecting this poorly illuminated region, as shown in figure 5(c).

Results
The results are described first in terms of the pressure drop. The velocity field for single-phase flow is described later, before the velocity and concentration profiles for the particle-laden cases are presented. Figure 6(a) displays the friction factor f as a function of the Reynolds number Re * for all the cases investigated in this study. The single-phase cases seem to agree reasonably well with the 16/Re * curve, which is the analytical solution for Poiseuille flow of a Newtonian fluid. This is noteworthy considering the fact that a complex expression is used to define Re * , see (1.1). There are measurable changes due to the addition of particles and to better appreciate them, the data in the figure are re-plotted as a percentage change compared to the analytical solution 16/Re * in figure 6(b). The deviation for single-phase flow is higher at the lowest Re * for HYS fluid but drops to values very close to the analytical solution for higher Re * . At particle volume  (b) shows the same plot but with a Reynolds number Re e * accounting for the additional suspension viscosity due to particles: the filled symbols are based on the Eilers fit and the unfilled symbols are based on a suspension having a modified yield stress and consistency (Chateau, Ovarlez & Trung 2008).

Pressure drop
fraction φ = 5 %, within the extent of the error bars, the increase in pressure drop is independent of the suspending fluid and Re * . The departure in pressure drop from the single-phase case is more pronounced for φ = 10 %. Note that for φ = 10 %, experiments are conducted only in LYS.
The increase in the viscosity of a suspension due to the addition of hard spheres can be semi-empirically described using the Eilers fit (von Eilers 1941) for a Newtonian suspending fluid under negligible inertia at the particle scale as µ e /µ = (1 + (5/4)φ/(1 − φ/φ Max )) 2 . This effective suspension viscosity µ e is only a function of the nominal particle concentration φ and the maximum packing fraction φ Max , assumed here to be equal to the value for random-close packing ≈65 %. Note that for φ 10 %, changing φ Max even by ±10 % does not produce any substantial change in the estimate of µ e . In the limit of low particle inertia and a uniform particle distribution, it may be used to predict the change in the friction factor for a Newtonian suspending fluid. Accordingly, in the inset of figure 6(b), a new effective Reynolds number Re e * is defined using µ e such that Re e * = Re * µ/µ e and the experimentally measured friction factor f at Re * is compared to the expected friction factor 16/Re e * for an effective fluid (see filled symbols). A good collapse of all the data points within ±5 % of the single-phase results is observed using this approach. A more appropriate rheological description of suspensions in YSF is provided by Chateau et al. (2008) and Mahaut et al. (2008). Their experiments and theoretical arguments suggest that a particle suspension in YSF behaves like a Herschel-Bulkley fluid with the same flow index n as the suspending YSF, but with a different equivalent yield stress τ y

√
(1 − φ)g(φ) and consistency κ (g(φ)) n+1 /(1 − φ) n−1 . The function g(φ) is the Krieger-Dougherty fit for suspension viscosity: g(φ) = (1 − φ/φ Max ) −2.5φ Max (Krieger & Dougherty 1959) and is virtually indistinguishable from the Eilers fit (mentioned above) up to the low φ = 10 % investigated in the present case. Using the above modification in yield stress and consistency, an equivalent Reynolds number, based on (1.1), can be defined and the corresponding laminar friction factor is (also) plotted in the inset of figure 6(b) (see unfilled symbols). Indeed, there are differences between the predictions using a simple effective viscosity (filled symbols) and the more reasonable modification in viscosity for YSF (unfilled symbols) proposed by Chateau et al. (2008). Both the above fits are based on the assumption of a homogeneous particle distribution and, as will be shown later, this is definitely not the case in our present experiments, where particles migrate to specific regions of the flow domain. Thus, there is no strong reason for an effective viscosity formulation to accurately predict the friction factor in laminar flows for a square duct. A more thorough investigation at multiple φ, particle sizes, Re * and flow geometries is therefore required to test the effectiveness of the above approach. In our previous work (Zade et al. 2018), we have observed that the friction factor in a turbulent flow of Newtonian fluid is a function of both the particle volume fraction and size. Also, from the results in Zade, Lundell & Brandt (2019b), we learnt that the predictions from the Eilers fit diverge at higher volume fractions φ 15 % for Newtonian turbulent flows.
A mean secondary flow is also observed in the duct, as will be shown in the coming sections. However, considering the weak amplitude of this flow, it is assumed to have a negligible effect on the friction factor, as also previously observed by Gao & Hartnett (1993). Another point to note is that the Reynolds number, if defined using the viscosity at the wall as Re µ Wall = ρU Bulk (2H)/µ Wall , is very similar in magnitude (within 10 %) to the Re * used in this study. The viscosity at the wall µ Wall , also used in the previous paragraph to estimate the particle scale Reynolds number, can be calculated using the flow curves in figure 2 and the average shear stress at the wall, estimated from the streamwise pressure gradient.

Velocity statistics for single-phase flow
For the single-phase flow, PIV measurements were performed in multiple spanwise planes. Mean profiles have been obtained by averaging in time and along the streamwise direction. An example of the results of such measurements is shown in figures 7 and 8 for HYS and LYS, respectively, for the same flow rate. The streamwise velocity profiles, normalised by the bulk velocity, display a flat region of constant velocity in the centre plane z/H = 0. This region of zero velocity gradient is representative of the solid plug and it shrinks while moving to a plane closer to the side wall z/H = 1. At the same flow rate, the plug is larger for the thicker fluid HYS due to the high yield stress (compare figures 7a and 8a). In the centre plane z/H = 0, the region of zero velocity gradient for the streamwise velocity has also a nearly zero wall-normal velocity i.e. there is no secondary flow inside the plug. The most striking feature comparing the two figures is that the intensity of the wall-normal velocities increases for the thinner fluid LYS. The plug region is smaller for LYS as deduced from the corresponding streamwise velocity profiles, allowing for a larger yielded region and hence, larger and stronger secondary flows. The action of the plug in damping the secondary flow was also observed in the simulations by Letelier et al. (2018).
A rectangular duct flow is symmetric about the orthogonal axes z − z and y − y. In the special case of the square duct, the flow is also symmetric about its diagonals. Thus, measurements of two velocity components in a limited number of spanwise planes (as in the 2D-PIV used here) can be used to calculate all the three velocity components at other locations in the duct. This approach is illustrated in the schematic of figure 9. Thus, 2D-PIV measurements in N planes (in figure 9, N = 5 denoted by darker lines) can yield the three-dimensional velocity at (2N − 1) 2 locations.
The above technique is applied using N = 9 planes for three different flow rates or Re * of the thicker fluid HYS as shown in figure 10. The secondary flow can now be clearly visualised: it carries the high momentum fluid from the yielded regions around the core to the centre of the wall and low momentum fluid from the corners towards the core. The intensity of this secondary flow increases with increasing Re * , coupled with the reduction of the plug size.
The extent of the plug at the duct core is determined as follows: a two-dimensional polynomial of sixth order is fitted to the streamwise velocity data (coefficient of determination R 2 0.99). A no-slip condition is assumed at the boundaries. Higher-order polynomials did not show any significant improvement in the quality of the fit. Using this fit, the streamwise velocity field was smoothly mapped on a plane of higher resolution (400 × 400 points). The norm of the velocity gradient (∂U/∂y) 2 + (∂U/∂z) 2 , using a simple forward difference scheme, is evaluated at each of these points and the plug is designated at those points where the norm is less than 1 % of its maximum value (occurs at the wall centre). This procedure naturally suggests the existence of small unyielded plug-like regions also at the corners where the velocity gradients are small. Despite the likely presence of these stagnant zones, as found in many previous numerical studies (Saramito & Roquet 2001), they are not shown here due to their proximity to the corners where PIV measurements are not performed. Nevertheless, their size is assumed to be very small for the Bi used in this study. The shape of the moving plug in the core region appears to be near circular, more so for higher Re * or, equivalently, lower Bi. Figure 11(a) shows the variation of Bi with Re * . For the same Re * , the thicker fluid HYS has a larger Bi. As mentioned previously in the Introduction, if a Bingham number Bi based on the HB parameters is defined as Bi HB = τ y /(κ(U/L) n ), it ranges from 0.36 to 0.72 for HYS and from 0.14 to 0.43 for LYS i.e. around twice the value of the Bi shown in figure 11(a). From Saramito & Roquet (2001), it can be expected that, with increasing Bi, both the moving plug at the centre and the stagnant plug at the corners grow. The growing plug in the core would become non-circular before finally merging with the plugs at the corner at some critical Bi, whose value for a strictly Bingham fluid was found to be 4/(2 + √ (π)) ≈ 1.06 (Mosolov & Miasnikov 1965). At this stage, the flow stops. Figures 11(b) and 11(c) show the change in the ratio of the maximum to bulk velocity for different Re * and Bi respectively. By virtue of a higher yield stress, at the same Re * , HYS has a larger plug than LYS and, hence, a flatter velocity profile i.e. a smaller U Max /U Bulk (see figure 11b). Figure 11(c) shows that, at the same Bi, HYS has a smaller U Max /U Bulk i.e. a larger plug region. This is in agreement with Huilgol & You (2005), who observed that the plug region is larger for a fluid with a lower power law exponent n, for fixed Bingham number Bi. From figure 2, it is deduced that LYS has a higher n than HYS but a smaller plug, as indicated by the higher U Max /U Bulk .
Before moving to the results for the particle-laden cases, a comparison between the velocity profiles from the simulations of Saramito & Roquet (2001)  scaling for the velocity: U Saramito = H 2 ( P/ x)/(2µ), µ being the Newtonian viscosity after yielding. A practically useful fit between U Saramito and U Bulk was provided, which is used to plot their profiles (see solid lines in 12). In contrast to our experiments with Carbopol, which is a shear-thinning Herschel-Bulkley fluid and has viscoelastic properties both before and after yielding, the simulations of Saramito & Roquet (2001) use a Bingham model i.e. no shear thinning in the yielded zone. Moreover, in the absence of viscoelasticity, they do not observe secondary flow, which is instead seen in our experiments. So, in light of the above differences in the rheology of the fluids, the flow field in our experiments is different than the simulation results of Saramito & Roquet (2001).   (2001), using a Bingham model, are shown (solid lines) for a few representative cases. Note that in order to aid visualisation, similar symbols correspond to nearly the same Bi for both the experimental fluids.
3.3. Particle-laden flow field The greatest advantage of using particles with the same refractive index as the suspending fluid is the ability to assess the spatio-temporal particle concentration even at high volume fractions. The different panels of figure 13 depict the instantaneous distribution of particles for φ = 5 % in a plane close to the side wall of the duct (z/H = 0.9) for increasing Re * . A plane closer to the wall is chosen since most of the particles are seen to migrate towards this region. At low flow rates, most of the particles migrate towards the top and bottom of this plane i.e. to the corners. Additional particles migrating towards the corners surround the existing particles, thus forming layers. This is especially visible at φ = 10 % and will be shown in the next figure. A distinct change in their distribution is observed as the flow rate is increased: particles migrate away from the corners towards the centre of the near-wall plane. Since the total number of particles in this near-wall plane decreases with increasing flow rate, particles also migrate away from the side wall towards the duct core. Particles detected using the algorithm described earlier are outlined with a blue colour. The difference in size of the particles is not due to the variability in their physical diameter but due to the fact that only the portion of the particles intercepting the laser sheet is illuminated. Several instantaneous snapshots as those reported in figure 13 are processed and the average particle concentration profiles shown in figure 14. The first row (figures 14a-14c) corresponds to φ = 5 % with the two symbols corresponding to two spanwise planes z/H = 0 and 0.9, as mentioned in the caption. The high concentration near the corners at the low flow rate (see the plane z/H = 0.9), as seen in the instantaneous snapshots, is quantitatively confirmed in figure 14(a). Particle layering is also observed in the form of multiple local maxima, whose values decrease away from the walls. The wall-normal size of each peak closely matches the diameter of the particle. At the highest flow rate, particles move away from the corners and migrate towards the centre of the side walls, as seen from the nearly flat concentration in this region in figure 14(c).
Gravitational effects are apparent from the asymmetry in the concentration distribution; there are more particles at the bottom than the top. However, this asymmetry is reduced at higher flow rates when the particle settling velocity reduces in comparison with the flow velocity. In the same figures 14(a)-14(c), it is apparent that the concentration in the centre plane (z/H = 0) is comparatively lower than close to the side wall. A small peak in concentration is observed very close to the core, inside the plug, at lower flow rates, which vanishes completely at the highest flow rates, where the plug is very small. At these high flow rates, concentration peaks are seen at the top and bottom of the centre plane. At high flow rates, the overall picture suggests the existence of a ring-like distribution with negligible concentration at the core and corners. On the other hand, at low flow rates, particles are seen to find an equilibrium position predominantly at the corners and in the centre of the relatively larger plug. These observations will be confirmed using shadowgraphy and explained using the secondary flow, the extent of the plug region and inter-particle collisions in the discussion section. The second row of the same figure (panels 14d-14f ) displays data from the flow at volume fraction φ = 10 %; it offers the same picture as above but at a higher concentration. Gravitational settling effects are more pronounced here. In this regard, an interesting point to note is that, when the flow is stopped by switching off the driving pump, particles come to rest in a very short time (a fraction of the time needed to stop when flowing in tap water at the same flow rate) due to the high viscosity of the suspending liquid. If the pump is not switched on, the particles remain in the same positions for a period of days, indicating that the yield stress of the fluid is higher than the buoyancy forces of individual particles. The yield-gravity parameter Y = τ y /((ρ p − ρ f )gd p ) for the fluid-particle system of this study is around 4, which is substantially greater than the critical value for a particle to move (between 0.04 and 0.2) (Chhabra 2006). Thus, the sedimentation observed in figure 14 is due to particles falling down through the suspending material that has yielded under shear during flow, as described in Ovarlez et al. (2012). Also, particles migrating towards the four corners can form packed layers near the wall and make physical contact with each other. This can create a frictional torque at the point of inter-particle contact which, along with the presence of the confining wall, can inhibit particle rotation (due . Change of mean streamwise velocity profiles due to introduction of particles at two spanwise planes: z/H = 0 (centre plane) and z/H = 0.9 (plane close to the side wall). Panels (a-c) corresponds to φ = 5 % and panels (d-f ) corresponds to φ = 10 %. The filled and hollow symbols represent the fluid and particle velocity respectively, and the solid lines represent the fluid velocity when no particles are present i.e. for φ = 0 %. Panels (a-c) share the same legends mentioned in (a,b) and (d-f ) share the same legends mentioned in (d,e).
to mean shear). Such a behaviour is visually observed in our experiments and often, particles in mutual contact intermittently slide at the wall. The flow-dependent non-uniform particle distribution causes a change in the mean streamwise velocity, as shown in figure 15. As in the previous figure, we display the particle and fluid velocities in two spanwise planes with the top row corresponding to φ = 5 % and the bottom row to φ = 10 %. The single-phase reference case at the same Re * is also shown for comparison.
The most striking changes happen for the low flow rates (figures 15a and 15d), when the particles occupy the corners. In these regions of high particle concentration, the fluid velocity is reduced as compared to its single-phase counterpart (see z/H = 0.9 plane in figures 15a and 15d). To make up for the velocity deficit in the near-wall planes, the velocity in the plug region increases (see z/H = 0 plane in figures 15a and 15d) and, as a consequence, the size of the plug reduces. Again, sedimentation also affects the mean velocity profiles, more for the higher φ. The particles move with approximately the velocity of the fluid in the centre plane but lag significantly the fluid phase in the near-wall plane, especially near the wall bisectors. It should be recalled that, at these low flow rates, most of the particles are concentrated at the corners and there are very few particles near the middle of the side walls, as seen from the concentration profiles in figures 14(a) and 14(d); we speculate that these slowly moving particles near the corners occasionally move towards the centre of the side walls, thus contributing to the apparent lag in the particle velocity. With increasing the flow rate, the particles redistribute into a porous ring-like region between the corners and the core, and, as a result, the mean velocity distribution has smaller deviations from the corresponding single-phase cases. Indeed, the plug size reduces with increasing particle concentration. At the highest flow rate considered, the plug in the single phase is already very small and hence, the reduction is not significant. A small but finite slip is also seen near the top and bottom wall of the centre plane due to the finite size of the spherical particles. This slip, due to finite-size effects, is more pronounced at the higher wall shear that occurs in turbulent flows (Zade et al. 2018).
The presence of particles introduces unsteady disturbances in the surrounding fluid field. The streamwise and wall-normal components of the velocity fluctuations associated with these disturbances are displayed in figures 16 and 17, following the same scheme as in the last two figures, as also detailed in the captions. Compared to the single-phase case, the streamwise and wall-normal fluctuations are higher. For the lowest flow rates, sedimentation causes rather high streamwise fluctuations at the bottom of the centre plane (see z/H = 0 in figures 16a and 16d).
At higher flow rates, the fluctuation profiles become increasingly symmetric. Considering the duct centre plane, streamwise fluctuations reach peak values close to the top and bottom walls, which are also the locations of high particle concentration. Interestingly, the fluctuations are higher for the intermediate Re * = 53 than at higher Re * = 156. This is perhaps the result of particles arranging themselves in an intermediate configuration between accumulations at the corners observed at low flow rates and the ring-like structure of the high flow rates. In other words, at such Re * , the spread of the particles in the duct is higher as also illustrated in figure 13 and hence for the same φ, the unsteadiness in the fluid field is larger. At higher flow rates, most of the particles cluster near the wall bisectors, the frequency of particles passing in this region is increased and hence, the unsteadiness in the velocities, or departure from the mean velocity, is not so pronounced leading to lower root mean square values in the near-wall planes. Measurements at additional planes will be needed to prove the above qualitative picture; nevertheless, the presence of increased streamwise disturbances in the presence of particles is clearly demonstrated.
The appearance of streamwise fluctuations is also accompanied by wall-normal velocity fluctuations, as depicted in figure 17. Settling of particles at lower flow rates (refer to figures 17a and 17d) leads to higher wall-normal fluctuations near the bottom of the duct while at the highest flow rates, the profiles are symmetric with nearly the same magnitudes for the two highest Re * (compare figure 17b-17c and 17e-17f ).
Despite the presence of intense streamwise and wall-normal fluid velocity fluctuations, the correlation between these two components, i.e. the Reynolds shear stress u v , does not drastically increase when compared to the corresponding single-phase case, as shown in figure 18. For reference, the peak values of u v /U 2 Bulk in Newtonian turbulence is approximately 3-5 × 10 −3 . Thus, the fluctuations generated in the flow are not a signature of the classical turbulence in a duct distinguished by a larger Reynolds shear stress, but is an unsteady flow perturbed by solid particles.

Discussion
In order to shed more light on the overall particle distribution in the duct, the quantitative measurements in figure 14 are complemented by flow visualisations obtained with shadowgraphy experiments. The set-up is sketched in figure 19, where the duct is illuminated by a point source of light such as an LED. The slight non-uniformity in the refractive index created by the particles inside the suspending medium casts a magnified shadow over the entire volume of a section of the duct on a screen, observable in a dark room. This shadow is captured by a camera, as shown in figure  FIGURE 21. Schematic of the particle distribution in the cross-stream section of the duct. The arrows and streamlines representing the normalised secondary flow and the extent of the plug are from the actual experimental measurements for single-phase flow at the corresponding flow rates. The contours represent the streamwise velocity normalised by the bulk velocity. Note that the colour bar is the same as in figure 10. The particle distribution is schematically reconstructed from the shadowgraphs and PIV measurements.
used to sketch an approximate particle distribution pattern at low and high flow rates, see figures 21(a) and 21(b). The increased particle concentration near the bottom wall for low flow rates, due to gravitational forces, is also reproduced in figure 21(a). For consistency, the same number of particles is depicted in both the figures. Note also that visualisations in additional planes besides those reported in figure 14 using the PIV laser sheet have contributed to draw the sketches in figure 21.
In order to explain the observed concentration distribution, it is essential to acknowledge the presence of competing forces of an inertial and viscoelastic nature. Scaling analysis can be used to demonstrate the significance of the different forces. In a fluid without elasticity, when the Reynolds number at the particle scale Re p,γ = ργ (d p /2) 2 /µ f becomes 1, the shear gradient lift forces F L ∝ ρ fγ 2 d p 4 (Asmolov 1999) will push the particle away from the core, whereas the repulsive wall forces will act to displace the particle away from the walls. As a result, an isolated particle reaches a stable equilibrium position somewhere in between the centre and the wall, a phenomenon famously known as the Segré-Silberberg effect (Segre & Silberberg 1961). In a square duct, due to repulsion from adjacent orthogonal walls, the stable equilibrium position is at the centre of the walls (Stoecklein & Di Carlo 2018). In our experiments, the maximum value of Re p,γ occurs at the wall and is estimated to be ≈0.1 at 2 l min −1 , and ≈6 at 40 l min −1 for LYS fluid flow. Thus, for the lower flow rates the inertial migration towards the centre of the wall is estimated to be very weak. On the other hand, one can estimate the relative magnitudes of viscoelastic forces F VE ∝ (N 1 + 2N 2 )d p 2 (see D'Avino et al. 2012). From De Cagny et al. (2019), the normal-stress differences can be approximated as ∼τ 2 /G . Thus, For a given flow rate, the force ratio is maximum at the wall since the denominator is some form of apparent viscosity that should be lowest at the wall. Thus, if the wall-shear stress τ w = τ y /Bi and a nominal shear rateγ = U Bulk/H are used to evaluate the above force ratio, it ranges from 10 −2 , for the lowest flow rate to 1, for the highest flow rate. This indicates the relative dominance of viscoelastic forces over inertial lift forces at the lowest flow rates. It also explains why particles migrate towards the corners (see illustration in figure 21a) which is a typical behaviour in inertialess (Re p ∼ 0) viscoelastic fluids (Villone et al. 2013), compared to inertial migration towards the centre of the walls in a square geometry. The observed migration is, in particular, due to the non-uniform distribution of the first normal-stress difference, being maximum at the wall centre and minimum at the corners and core (Yang et al. 2011). Migration towards the core is not as pronounced due to the presence of the non-yielded core with its infinite viscosity. Hence, the non-zero concentration inside the plug is most likely due to particles that were trapped in the bulk at the inlet to the duct.
With increasing flow rates, Re p as well as F L /F VE , increases i.e. inertia starts to affect the migration and drives the particles to the centre of the walls. However, our experimental measurements and shadowgraphy observations show that, instead of being concentrated only at the wall centres, the particle distribution is similar to that sketched in figure 21(b). Thus, we have reason to believe that, in addition to inertial migration, the particle distribution is also influenced by the increasing intensity of secondary flow, induced by the second normal-stress differences. The competition between the first normal-stress driven migration and secondary flow driven migration was already explored in Villone et al. (2013), who found that with increasing flow rates, the velocity of the secondary flow increases faster than the migration velocity due to first normal-stress differences. As a consequence, particles may now also be dragged by the helicoidal motion of the secondary flow (also see Lim, Nam & Shin (2014b)). The secondary flow measured for the single-phase cases is most likely changed by the presence of the particles; however, it is not reported in this study since measurements displayed excessive noise and showed consistent behaviour. Nonetheless, in a Newtonian suspending medium, the modifications of the secondary flow due to cross-stream migration of finite-size particles in a square duct have been shown in numerical studies as in Kazerooni et al. (2017). Previously, Ramachandran & Leighton (2008) proved that strong secondary normal-stress differences, and hence secondary flows, can be produced in concentrated suspensions due to an anisotropic particle microstructure. Finally, due to high local concentration and non-negligible Re p at high flow rates, multi-particle interactions may also influence particle dispersion. Non-negligible Re p at higher flow rates may lead to inertial effects like shear thickening (Picano et al. 2013), which can influence the friction factor. We note that in our experiments at the higher flow rates, the plug in the core has a negligible size or is perhaps non-existent, and the flow therefore approaches the case of particles in a shear-thinning viscoelastic fluid.
The simple shadowgraphy visualisations also help in assessing the evolution of the particle distribution along the duct streamwise length. Figures 20(a) and 20(c) show the shadowgraphs of the particle distribution very close to the entrance of the duct (at x/H ≈ 12) at low and high flow rates. At the low flow rate, particles seem to enter the duct while being mostly concentrated in the bottom half. This entry profile may either be due to gravitational forces or because of the secondary Dean motion from the bend at the exit of the tank in our set-up. The installed static mixer meant to cancel these vortices is, perhaps, less effective at such low Re * . Nevertheless, as the flow develops, the initial disturbances are dissipated and particles attain an equilibrium concentration, as shown in figure 20(b), which is captured very close to the location where the PIV measurements are performed. Thus, particles migrate to the top against gravitational forces. Particles travelling at a higher speed in the core are seen to move without any change in their relative positions on account of being locked up inside the solid plug. On many occasions, clusters of two or more particles are seen to move together inside the plug (which was observed before in the case of particle sedimentation by Chaparian, Wachs & Frigaard (2018)). The distribution profile just described is attained much earlier than x/H ≈ 160 and persists along the remaining length of the duct. Projection from only one direction is insufficient to give full information on the particle distribution and, hence, the duct was also illuminated from the bottom to check a second projection. And indeed, a picture very similar to that in figure 20(b) is obtained, thus confirming the preferential accumulation of particles at the corners.
Focussing on the highest flow rates, it appears that the concentration profile is more uniform at the entrance (see figure 20c) and rapidly develops to an equilibrium pattern, as seen in figure 20(d). By visualising fast and slow particles, it appears that very few particles travel near the core, which has the highest velocity, or at the corners, which have the smallest velocity. Most of the particles travel in the intermediate region between the core and the corners. Moreover, the particle motion appears very erratic with instances of rapid change in their velocities. This indicates that strong particle-particle and particle-wall interactions are most likely relevant at these concentrations and flow rates. The presence of finite-size particles also causes an increase in the local strain rates, and the corresponding reduction in the local viscosity for the shear-thinning fluid may intensify the strength of these interactions, which may be an important mechanism for generating the fluid velocity fluctuations shown in figures 16 and 17.
To conclude, we note that from a visualisation perspective, another noteworthy technique can be used, exploiting the property of elastoviscoplastic flows to come to a rather quick stop after switching off the pump, as mentioned earlier. As a result, the particles get trapped in their flowing configuration without much change of the inter-particle distance, thus making it possible to assess their position and hence, the concentration distribution. This can be conveniently done by translating the laser sheet across such a frozen flow and capturing still photos which can be later stacked together to reconstruct the full three-dimensional concentration field. Such experiments are not performed here but can be advantageously used for these kinds of fluids.

Conclusion
We have presented an experimental study of the flow of an elastoviscoplastic fluid, Carbopol gel, at two different concentrations -one with a high yield stress and the other with a low yield stress -and studied the laminar single-phase flow and the flow in the presence of finite-size spherical particles at relatively high concentrations inside a square duct. Pressure drop, mean and fluctuating velocities and concentration fields have been measured using a combination of PIV and PTV; these optical techniques have been possible due to the refractive-index matching of the fluid-particle mixture. To facilitate comparative simulations, the rheological properties of the fluids investigated are described in terms of the steady-state flow curves and viscoelastic moduli. A measure of normal-stress differences is provided from the recent work of De Cagny et al. (2019).
The experimental pressure drop is in agreement, over the range of flow rates considered, with the laminar friction factor based on the semi-empirical Reynolds number Re * (see (1.1)) by Liu & Masliyah (1998). By making use of the symmetries in a square duct, we have presented the full 3-component velocity field (at moderate spatial resolution) for single-phase flow using two-dimensional PIV in multiple spanwise planes (at high spatial resolution), evidencing the existence of a secondary flow due to the viscoelastic properties of the fluid after yielding. Using a low threshold value of the velocity gradient, we have identified the near-circular plug region in the core, whose size reduces with increasing the flow rate or equivalently Re * . At the same Re * , the fluid with a higher yield stress has a higher Bingham number Bi and a larger plug. Also, at the same Bi, the ratio of the maximum streamwise velocity to the bulk velocity is larger for a fluid with a higher power law exponent n. The strength and size of the streamwise vortices representing the secondary flow increases for reducing plug sizes. Shear thinning and secondary flow cause a deviation of the mean streamwise velocity when compared with the simulations of Saramito & Roquet (2001) for a Bingham fluid.
Similar to Newtonian suspending fluids, the friction factor increases by increasing the particle volume fraction φ at fixed Re * . This increase seems independent of the rheology of the suspending fluid. A reasonably good prediction is obtained based on the effective viscosity of suspensions of rigid spheres in Newtonian fluid, like the Eilers fit. The matching with predictions based on the rheological modifications of YSF proposed by Chateau et al. (2008) is not so good. These differences are most likely due to a non-uniform particle distribution in the duct, which violates the assumptions of the above theory. At low flow rates, in the presence of a larger plug and weaker secondary flow, particles migrate inside the yielded zone, defying gravity, towards the four corners under the influence of the first normal-stress difference. Fluid (and particle) velocity reduces in the near-wall planes (where there is a high particle concentration) and increases in the plane of the wall bisector when compared to the corresponding single-phase flow at the same locations. This increase of the maximum velocity reflects in an increased shear in the core region and results in a smaller plug than for the single-phase flow. At high flow rates, in the presence of a negligibly small plug and a stronger secondary flow and inertial lift forces, particles arrange along a diffused ring between the core and the corners. A simple scaling of forces is provided to explain the dominance of viscoelastic forces at low flow rates and inertial forces at high flow rates. Particles are also seen to induce time-dependent fluid velocities, leading to uncorrelated fluctuations in the streamwise and wall-normal directions. Visualisation of shadowgraphs have, even though qualitatively, confirmed the above concentration profiles.
The importance and novelty of this experiment, when compared to other studies with particles in non-Newtonian fluids, is the ability to assess the velocity distribution of both the fluid and particle phases along with particle distribution at high bulk concentrations. This will, hopefully, help in future modelling attempts of such systems and shed light on the complex fluid-particle interactions present therein. These results are especially relevant to microfluidic applications, which share similar values of the relevant non-dimensional parameters, e.g. the Reynolds number Re * and confinement (or blockage) ratios and where resolved measurements are prohibited by the small length scales. As always, further investigation with different particle sizes, volume fractions and fluid rheologies, perhaps coupled with spatio-temporally resolved measurements, will help to further understand these complex systems.