Particle transport in laboratory models of bifurcating fractures

Accepted Manuscript Particle transport in laboratory models of bifurcating fractures Sojwal Manoorkar, Omer Sedes, Jeffrey F. Morris PII:

S1875-5100(16)30231-1

DOI:

10.1016/j.jngse.2016.04.008

Reference:

JNGSE 1422

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 27 January 2016 Revised Date:

1 April 2016

Accepted Date: 3 April 2016

Please cite this article as: Manoorkar, S., Sedes, O., Morris, J.F., Particle transport in laboratory models of bifurcating fractures, Journal of Natural Gas Science & Engineering (2016), doi: 10.1016/ j.jngse.2016.04.008. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

RI PT

Particle transport in laboratory models of bifurcating fractures Sojwal Manoorkar1 , Omer Sedes1 , and Jeffrey F. Morris1

Benjamin Levich Institute and Department of Chemical Engineering, The City College of New York, New York, NY10031, USA May 3, 2016

M AN U

Abstract

SC

1

Introduction

EP

1

TE

D

This article provides experimental analysis as well as numerical simulation of flow of liquids and suspensions through laboratory models of basic bifurcating channel geometries. Numerical results of single phase Newtonian and non-Newtonian (power law) fluid flow in 3-dimensional square channels show good agreement with experimental data. The particleladen experiments are carried out with neutrally buoyant suspension of spherical particles for different inlet solid volume fraction and different Dk /D⊥ where Dk is the width of the straight branch and D⊥ is the width of side branch of the channel for a range of Reynolds numbers, 0 < Re < 800; Dk = D0 , with D0 the width of the inlet branch. The focus of the multiphase flow is on the role of the size ratio of the side channel to the primary or straight channel, and the relative size of the particle to the channel widths. The Reynolds number is defined as Re = ρDη0 U , where U is the average velocity at inlet, while ρ and η are density and viscosity of the fluid, respectively.

AC C

Suspension flows through networks of channels or conduits appear in magma flow, in engineering contexts, and in study of hemodynamics. This study considers the pure and particle-laden flow in a laboratory model of a bifurcation in a fractured formation, with motivation from the petroleum exploration industry, where hydraulic fracturing enhances the production rate of a hydrocarbon-bearing zone. The basic process in hydraulic fracturing involves pumping liquid and solid particles at high pressure to induce fractures. The solid particles remain after the pressure is released to create a high-permeability path to the well bore (Garagash et al., 2015; Garagash & Lecampion, 2014; Warpinski et al., 2009). The induced fracture may encounter natural fractures and the issue of the flow at junctions in fracture networks motivates the present study. It is interesting that network flows of 1

ACCEPTED MANUSCRIPT

AC C

EP

TE

D

M AN U

SC

RI PT

suspensions are studied extensively in the context of physiological hemodynamics (Pedley, 2008; Pries et al., 1990), as blood is a suspension of red blood cells. It is difficult to consider the realistic porous medium environment of hydraulic fracturing. Much of the work in this field involves flow in porous media, considering mass transfer and permeability of porous rock (Blunt, 2001), or the mechanism of fracture propagation (Bohloli & Pater, 2006). While often simplified, the fracture structures are threedimensional in nature (Warburton, 1989; Piggot, 1997; Berkowitz & Adler, 1998). Single phase flow and solute transport in networks have been studied using percolation theory (Mo et al., 1998). Some studies find preferential paths for fluid flow in fractures (Kwicklis & Hearly, 1993). The present study has a different objective than these, as we are concerned with flow at bifurcations between distinct fractures in fractured media. We have studied in a laboratory model of a fracture bifurcation how various factors affect the manner in which flow of liquid and particles split at a bifurcation. These include the relative widths of downstream branches, solid particle size and the effect of particle migration upstream of the bifurcation. We present results on single-phase and particulate suspension through a single bifurcation at a right angle (90◦ ) from the upstream channel, i.e. T-bifurcation. Literature addressing study of flows in such a T geometry is predominately related to piping networks and hemodynamics. Early studies of flow in such geometries focused on energy and pressure losses. Iwanami & Suu (1969), Iwanami et al. (1969) studied the characteristics water and slurry flow in a circular cross section pipe fitting for right angles, both numerically and experimentally. The flow in a T-geometry at elevated Reynolds number separates at the sharp corner, forming a first separation zone in the side branch. A second separation zone is formed in the straight channel at the wall opposite the side branch (Pedley, 2008). These zones were reported by Iwanami & Suu (1969) and Sierra-Espinosa et al. (2000). Neary & Sortiropoulos (1996) studied the separation zones numerically to predict longitudinal vorticity contours and streamlines. They compared their numerical results with experiments by Liepsch & Moravec (1982), who carried out experiments using laser-Doppler anemometry to determine the effect of Reynolds number on velocity profile and separation zones. Miranda et al. (2008) studied steady and unsteady flow of not only Newtonian but also non-Newtonian (Carreau model) fluid through a 2D T-bifurcation. The noted separations and closed streamlines raise a number of unresolved issues in particleladen flows, as described in the flow past a bluff body by Haddadi et al. (2014). There is very limited work on suspension flows in bifurcating geometries. Bugliarello & Hsiao (1964) described the solid concentrations in downstream branches of a T geometry (and also limited other angles) for varying imposed flow ratio in the two branches as a function of total flow rate and solid concentration. Roberts & Olbricht (2006) studied bifurcation flows as a particle separation technique, considering very small Reynolds number. Xi & Shapley (2008) measured the velocity and concentration profile in asymmetric bifurcations using nuclear magnetic resonance imaging; similar geometries have been modeled by continuum approaches (Ahmed & Singh, 2011). Pries et al. (1990) developed a model to simulate blood flow through large microcirculatory networks using viscosity model for 2

ACCEPTED MANUSCRIPT

2

AC C

EP

TE

D

M AN U

SC

RI PT

blood. One issue of interest in blood flow arterial network is the increase in red blood cell concentration in the downstream branch with higher flow rate, known as Zweifach-Fung effect. This was examined by Doyeux et al. (2011). Blood flow involves deformable objects, and this leads to a near wall cell-free layer which impacts upon the relative flow rates in the downstream branches. This role of migration is found to be significant in solid particle suspensions, but is due to other mechanisms. In the present study the spherical particles are non-deformable at the flow-induced stresses applied. Under these conditions, there are two primary migration phenomena. One of these is shear induced migration (Koh et al., 1994; Hampton et al., 1997) which is related to suspension normal stresses (Shen et al., 1992; Nott & Brady, 1994; Morris & Boulay, 1999; Boyer et al., 2011) and is thus a factor primarily in relatively concentrated suspensions, which we do not consider in the present work. The second mechanism for migration is inertial focusing, which drives individual particles to preferred equilibrium positions in pressure driven flow in tubes (Segre’ & Silberberg, 1961; Han et al., 1999; Matas et al., 2004) and channels (Bhagat et al., 2008; Chun & Ladd, 2006; Matas et al., 2009). Miura et al. (2014) have considered inertial migration in a straight square channel, and their results are thus directly relevant to the present study. We find that inertial migration in the square channel plays a significant role on the distribution in the two channels as a function of Re. We compute, using the latticeBoltzmann method for suspensions, the particle migration in a straight channel and couple this with computed pure-fluid streamlines in the T-geometry flows in the present study. Such a coupling between the upstream concentration profile and the dividing streamline was previously considered by Yan et al. (1991) to analyze the fluid skimming and particle entrainment effect into a small side pore at low Reynolds number. Enden & Popel (1994) calculated the shape of the dividing streamline at the T-shaped bifurcation of two cylindrical tubes. In both of these studies, the authors considered a simple particle concentration profile, where they assumed a particle-free layer near the wall followed by a uniform concentration of particles, and superimposed this profile with their dividing streamline calculations to predict the concentration of particles downstream. Here we have considered the computed, more complex, particle profile which accounts for inertial migration upstream. The experimental and numerical methods are described in §2 and §3, respectively. Results of the experiments and comparison with numerical calculations are presented in §4 with certain microscopic details of pressure field and flow structure provided. This is followed by concluding remarks.

2.1

Experiment Materials

For the single phase experimental studies, a mixture of 50% by volume glycerol in deionized water with density ρ = 1.143 g cm−3 and viscosity η = 0.0086 Pa.s at 20 ◦ C is used as 3

ACCEPTED MANUSCRIPT

2.2

SC

RI PT

the Newtonian fluid. A non-Newtonian fluid relevant to fracturing flows is prepared by dissolving guar polymer (Halliburton WG-36) in a mixture of 17% by volume glycerol in deionized water. The polymer solution is stirred at 100 rpm for 4-5 hours. The pH of polymer solution was measured to be 7.5 at 25◦ C. We have studied guar polymer in the noted base fluid at weight fractions polymer of 0.0024 and 0.0042; this is 20 and 35 pptg (pounds per thousand gallons) in petroleum production terms. Suspensions consist of spherical particles of two different types which include polystyrene particles of density ρp = 1.05 g cm−3 and diameter d = 150 µm - 360 µm (Maxiblast) and PMMA particles of density ρp = 1.18 g cm−3 and d = 50 µm - 75 µm (Lucite International). These particles are mixed into density matching fluids, 17% glycerol by volume in deionized water for polystyrene or 65% glycerol by volume in deionized water for PMMA particles, so that rising and settling are minimal.

Flow experiments

AC C

EP

TE

D

M AN U

The schematic of the experimental setup is shown in figure 1. The experimental device involves different widths of the side branch of the T-channel. These flow channels are made by two techniques depending upon Dk /D⊥ due to mechanical issues. For Dk /D⊥ = 1 and 3, channels are fabricated using two acrylic sheets: the T-shaped channel is etched out of one of the sheets by a drill bit and the other sheet is glued on the top of the etched sheet to form a channel of width D0 = Dk = 2.4 mm and D⊥ = 2.4 mm or 0.8 mm. For Dk /D⊥ = 5, the channel is fabricated using three acrylic sheets in a sandwich fashion, with the T-shaped channel cut into one of the sheets using a laser cutter. The other two acrylic sheets are glued on the top and bottom of this sheet to form the bounded channel. Inlet and outlets of channel are at 14 cm from the midpoint of the T-junction in both cases so that Lk = L⊥ = L0 = 14 cm, yielding the channel aspect ratio L0 /D0 = 58. A high speed camera (Photron FASTCAM 1024 PCI) is vertically mounted at the T-junction to capture the images. Fluid is driven through the channel by a syringe pump (Harvard Apparatus PhD 2000) at varying flow rates, with 30 ml used for each experiment. The outlet flows from both branches are open to atmospheric pressure, and the flow is collected for the duration of the experiment. For every inlet flow rate (discussed below in terms of Reynolds number), five replicates are performed with the same fluid mixture to achieve statistically significant results. The ratio of the flux to the straight downstream channel relative to the inlet channel β = Qk /Q0 is used to characterize the flow condition, where Q is a volumetric flow rate; this quantity is subscripted to indicate whether bulk suspension, particle phase, or fluid phase is meant. We determine βsuspension as the ratio of suspension weight captured at the outlet of the straight branch to inlet branch, and we confirm the mass balance by summing the two outlets. The suspensions are then filtered to separate the particles, and the particles obtained from each sample are dried at 80 ◦ C for several hours in an oven and weighed. The ratio of the mass of the particles in the straight branch to the total gives 4

ACCEPTED MANUSCRIPT

Camera

L

L0 P0

D0= D

Q0

L

P

Q

SC

D

RI PT

Syringe pump

M AN U

P

Q

Figure 1: Schematic of experimental set-up; Pk = P⊥ = atmospheric pressure

2.3

Rheology

TE

D

βparticle . By difference from the total suspension, we calculate βfluid . We also captured images of the T-junction at 2000 fps with the resolution of 1024 pixels x 512 pixels over a linear dimension of 1.5 cm x 0.75 cm.

AC C

EP

The viscosity has been measured by rheometer. The rheology data is obtained for the nonNewtonian fluid using a strain-controlled rheometer (TA Instruments ARES-G2). The rheometer is equipped with parallel plates of 50 mm diameter with 1 mm gap between the two plates. A sample volume of approximately 2 ml was loaded on the geometry for a shear rate ramp between 1 s−1 and 700 s−1 for total time of 15 minutes. The rheology experiments are carried out at the same temperature as the flow experiments. Viscosity of the measured fluid sample is plotted against shear rate in figure 2 for the two polymer concentrations. A power law fit is made, with the parameters as shown in the plots. For non-Newtonian fluids, the definition of Reynolds number is generalized based on the concept of Metzner & Reed (1955). The ‘generalized Reynolds number’ has different forms depending upon the viscosity behavior of fluid with shear rate. For power law fluids, the following definition is used with ξ = 12 for infinite parallel plates which corresponds to a 5

ACCEPTED MANUSCRIPT

2

Experimental data for 20 pptg guar n-1 Power law fit as µ = Κ γ

Experimental data for 35 pptg guar n-1 Power law fit as µ = Κ γ

1

9 8 7

0.1 9

6

8

5

7

Coefficient values n = 0.68 K = 0.149

3

Coefficient values

3

n

= 0.53

K = 1.63

2

0.1

9 8 7

2

6 5 4

0.01

3 2

1

3

4

5 6 7 8

2

3

4

10

5 6 7 8

2

3

4

5 6 7 8

100

2

1000

3

4

5 6 7 8

1

2

10

Shear rate (1/s)

RI PT

5 4

4

)s.aP( ytisocsiV

)s.aP( ytisocsiV

6

3

4

5 6 7 8

2

3

4

5 6 7 8

100

1000

SC

Shear rate (1/s)

(a)

(b)

M AN U

Figure 2: Viscosity data for guar polymer solution (a) 20 pptg (b) 35pptg 2D channel and ξ = 7.1 for a square cross section conduit (Delplace & Leuliet, 1995): Regen =

ρDn U 2−n

n



K (24ξ + 1)/(24 + ξ)n

(1)

ξ n−1

Numerical simulations

EP

3

TE

D

We have also used oscillatory shear at small amplitude to determine the linear viscoelastic response for the guar solutions. Experiments have been performed at constant strain of 50% for the guar solution. In both cases, at very low frequency G00 > G0 , indicating liquid-like behavior and the cross-over frequency for lower concentration guar is larger which means a smaller relaxation time corresponding to 1/ωc , where ωc is the crossover frequency for G0 and G00 .

AC C

Pressure-driven laminar steady flow of a Newtonian fluid as well as power law non-Newtonian fluid in a 2D channel and 3D-square cross sectional channel with T-geometry was computed R numerically using a commercial solver COMSOL Multiphysics 5.2 (COMSOL, 2015). The computational domain was defined to match the geometric parameters of the experimental device, defining L0 = Lk = L⊥ = 1 , L0 /D0 = 58 and Dk /D⊥ = 1, 3 or 5. Additionally, the numerical calculations were also done for Dk /D⊥ = 0.5, 1.5, 2 and 10. The incompressible Navier-Stokes equations were solved in this computational domain, with fully developed laminar flow at the inlet and zero pressure at each outlet. Specifically, we used the singleR phase laminar incompressible flow solver of the COMSOL CFD module. The velocity field was discretized using quadratic elements while the pressure field was discretized using

6

ACCEPTED MANUSCRIPT

10

7

9 8 7 6

6

loss modulus ( G'') storage modulus (G')

5

loss modulus (G'') storage modulus (G')

5 4

4

3 3

1

9 8 7 6 5 4

1

3

9 8

2 7 6

0.1

5 6

7

8

9

2

3

2

10

0.1

3

4

5

6

RI PT

)aP( ''G ro 'G

)aP( ''G ro 'G

2

2

7

8 9

2

1

3

4

5

6

7

8 9

10

Frequency (Hz)

SC

Frequency (Hz)

(a)

(b)

M AN U

Figure 3: Linear viscoelastic response of the guar polymer solution at 50% strain (a) polymer weight fraction 0.0024 (20 pptg) and (b) polymer weight fraction 0.0042 (35 pptg).

AC C

EP

TE

D

linear elements. The streamline diffusion and crosswind diffusion options were enabled for stability. The discretized equations were solved using COMSOL’s iterative generalized minimum residual algorithm (GMRES). The relative tolerance was set to 0.001.Various computational meshes were used to establish mesh independence, with the finest mesh consisting of 587088 hexahedral elements for 3D and 248766 rectangular elements for 2D. The computations were performed at inlet Reynolds numbers of 0 < Re < 600.The simulation results were used to investigate the detailed flow structure and pressure variation at the bifurcation region. In order to obtain information on the particle distribution in the upstream of the T-junction, the pressure driven flow of a suspension of particles in a straight square cross section conduit were computed using the lattice-Boltzmann method for particle-laden flows (Ladd, 1994a,b). This method has been demonstrated for shear flow (Kulkarni & Morris, 2008; Haddadi & Morris, 2014) and pressure-driven flow (Haddadi et al., 2014) of suspensions up to φ = 0.35. Details of the method can be found in the references provided. We simulated the pressure-driven flow of a suspension of φ0 = 0.12 in a straight channel with a square cross section, with d/D = 0.1 at the range of Re used in experiments. Since the maximum Mach number (M amax = umax /c) observed in the simulations was 0.0487, the flow can be considered incompressible. The cross sectional (y-z) positions and the axial (x) velocities of the particles were recorded as the flow progressed, and we have used here the values when the mean flow had axial (x) travel of the same value of L0 /D as in the experimental inlet section. To analyze the influence of particle migration on βparticle , the flow split, these particle profiles were then superimposed with the cross-sectional pure fluid streamline data in the T-junction calculated at the same Re. Assuming particles at a position in the cross-section just before the bifurcation follow the fluid streamline which would 7

ACCEPTED MANUSCRIPT

otherwise pass through their center, we obtain an estimate of the particle split. This is simply the ratio of the flux of particles lying on streamlines following the straight channel to the particle flux on streamlines which are diverted to the side branch:

(jk · n)dA k βparticle = Z Z

= (j · n)dA

i=1 N X j=1

Results and discussion

4.1

. Upj

(2)

SC

4

Upi

RI PT

N

k X

ZZ

Single phase fluid

TE

D

M AN U

Figure 4a shows the variation of the fractional flow rate to the straight branch with the channel Reynolds number, Re (or generalized Regen as given by (1) for non-Newtonian fluid). At low Re, β is close to 0.5 indicating equal division of the inlet to the two outlets. This is expected for flow approaching the Stokes regime. As Re increases, more fluid follows the straight branch. The flow split of non-Newtonian shear thinning fluid with weak elasticity, which roughly follows power law viscosity, can be described very well with Regen as shown in figure 4: the flow split data for different concentrations of polymer solution collapse on the same curve as the Newtonian fluid data. It is observed that there is a slight dip in the value of β before it starts increasing with Regen . This behavior is seen only in the non-Newtonian fluids studied, as β was monotonically increasing with Re for the Newtonian single-phase liquids. A curve is fitted for the experimental variation of β with Re and it is observed that the dependence of β with increase in Re is describable by

EP

β=

1 , 1 + Ae−qRe

(3)

AC C

with parameter values given in Figure 4a for Dk /D⊥ = 1. These parameters depend on the channel geometry. Figure 4b shows the experimental results for Newtonian fluid for different ratios of the width of side channel to straight channel. For Dk /D⊥ = 3 or 5, β shows very weak dependence on Re as compared to Dk /D⊥ = 1. The fitting parameters for these experimental geometries are listed in table 1. Numerical simulations were done for a wider range of values of Dk /D⊥ and results are shown in figure 5a which again indicate the weaker dependence of β with Re with increase in width of the straight channel. The specific condition where Dk /D⊥ = 0.5 (side branch is wider than straight branch) shows that a higher volume of fluid flows through side branch for all studied values of Re compared to straight branch giving β < 0.5. Though it should 8

ACCEPTED MANUSCRIPT

1.2

RI PT

0.8 /

D|| /D⊥= 1 D|| /D⊥= 3 D|| /D⊥= 5

D|| D⊥ = 1

0.7 1.0

Coefficient values q = 0.00067 A = 1

0.5

β

β

0.6

0.8

0.6 0.3

Newtonian fluid 20pptg polymer solution 35pptg polymer solution Exponential fit

0.2

0.4 200

400

600

800

1000

0

200

M AN U

0

SC

0.4

Re or Regen

(a)

Re

400

600

800

(b)

TE

D

Figure 4: Experimental data for single phase fluid (a) Newtonian and non-Newtonian polymer solution for Dk /D⊥ = 1 (b) Newtonian fluid for different Dk /D⊥ = 1.

A 1 0.069 0.0157

AC C

EP

Dk /D⊥ 1 3 5

q 0.00067 0.00014 0.00008

Table 1: Fit parameters for β = (1 + Ae−qRe )−1 in the different geometries studied.

9

1.0

0.8

0.8

0.6

0.6

β

β

1.0

0.4

0.4 D||/D⊥ = 0.5 D||/D⊥ = 1 D||/D⊥ = 1.5 D||/D⊥ = 2 D||/D⊥ = 3 D||/D⊥ = 5 D||/D⊥ = 10

0.2

0.2

0.0

0.0 0

100

200

300

400

500

0

600

2

4

RI PT

ACCEPTED MANUSCRIPT

D||

6 / D⊥

8

10

12

SC

Re

(a)

(b)

M AN U

Figure 5: Numerical results for 3D channel (a) For different Dk /D⊥ (b) β for Re = 200

AC C

EP

TE

D

be noted that as Re → ∞, inertial effects become larger giving eventually more fluid in straight branch and β > 0.5. In figure 5b, we have plotted the variation of β with Dk /D⊥ at Re = 200. As Dk /D⊥ increases, β increases, approaching saturation value of β = 1 indicating negligible fluid flow in side branch when width of side branch is very small compared to straight branch. This figure is significant for determination of comparative flow rate in fractures with different size. All the above numerical simulations were performed for a 3D square cross sectional channel. We have carried out 2D calculation for the same geometry. Figure 6 shows the comparison for experimental data with 3D and 2D calculations. Experimental data agrees well with 3D numerical results but varies significantly from 2D results indicating the difference in basic structure of the dividing streamline. The numerical value for β is higher in both the cases (Dk /D⊥ = 1 and 3) for 2D as compared to 3D results which is explained by studying detailed structures of streamlines in the channel later in the section. To examine the detailed features of single-phase flow in the T-bifurcation, we obtained the streamlines from COMSOL determined flow solution. Figure 7 shows the flow streamlines in a square cross-section channel. These streamlines are at all depths in the direction normal to the page. At Re = 10, there is no flow separation and streamlines are essentially . equally divided, giving β = 0.5. At Re = 400, a separation zone is seen at the entrance to the side branch. A second separation zone forms at the opposite wall from the side channel. Streamlines at higher Re show that the fluid in the side branch is withdrawn from the viscous layer near the wall region of the inlet branch, whereas the core of the flow continues in the straight branch. The gauge pressure traces along the mid-line of the straight (line AB) and side (line CD) branches are shown in figure 8a. The coordinate X = 0 is the inlet of T-channel, (X = 280) corresponds to outlets of side and straight branches where pressure is equal 10

ACCEPTED MANUSCRIPT

0.70

1.00

/

D|| D⊥ = 3

/

D|| D⊥ = 1

0.65

RI PT

0.98

0.60

β

β

0.96 0.55

0.94 0.50

2D numerical results 3D numerical results Experimental data

0.40

0.90 0

200

Re

400

600

800

0

SC

0.92 0.45

100

200

300

400

2D numerical results 3D numerical results Experimental data

500

600

Re

(b)

M AN U

(a)

AC C

EP

TE

D

Figure 6: Comparison of experiments to numerical results (a) Dk /D⊥ = 1 (b) Dk /D⊥ = 3

(a)

(b)

Figure 7: Fluid streamlines in Newtonian fluid (a) Re = 10 (b) Re = 400

11

ACCEPTED MANUSCRIPT

250

10000 Line AB Line CD

200

Line AB Line CD 8000

85 84

0

P

80

A

C

4000

78

B

136

138

140

142

144

X

50

2000

D 0

0 0

50

100

150

X

200

250

0

50

RI PT

( ⊥

(

81

79

100

6000

⊥

(

82

P ro ||P - P

P ro ||P -

P ro ||P - P

⊥

83

)aP

)aP

)aP

150

100

150

200

250

SC

X

(a)

(b)

M AN U

Figure 8: Pressure from computed flow along the midlines of channels for Newtonian fluid in T-bifurcation of Dk /D⊥ = 1 (a) Re = 10 and (b) Re = 400. The numerical values correspond to the flow of the more viscous suspending fluid (higher glycerol content) in the actual channels used in the experiments.

4.2

Suspension flow

D

(P = 0). At X = 140 which is at the midpoint of T-junction, the straight channel shows an increase in pressure, whereas pressure decreases in the side channel. For higher Re, this change in pressure is much higher, in agreement with the trends expected of a simple Bernoulli equation analysis.

AC C

EP

TE

We carried out experiments for particulate flow at different concentrations of solid, φ0 . Figure 9 shows the results for φ0 = 0.12 for Dk /D⊥ = 1; these results include suspension in both Newtonian fluid and the shear-thinning guar polymer solution. We find that βsuspension differs significantly from the single phase fluid split, which implies the underlying streamlines of the mean mixture flow are altered. Also, βparticle 6= βsuspension , a result due at least in part to the effects of upstream inertial migration to distribute the particles nonuniformly across the channel; we show this through numerical simulation later. In figure 9b it can be seen that βparticle is higher for the suspension in the guar solution (polymer weight fraction 0.0042) at the same Re. This can be explained by the more pronounced migration of particles towards the central core of the square cross section in the guar solution as shown in experimental images in figure 10 (a,b). The image of the guar solution is less clear because the viscosity of the solution is higher, requiring higher velocity to achieve the same Reynolds number, while the frame rates for imaging are the same. Figure 10 (c,d) shows the computed flows for the single phase Newtonian fluid and for the power-law description of the guar solution. These do not differ greatly, and thus the

12

ACCEPTED MANUSCRIPT

0.8

0.8 d/D0 = 0.1

0.7

0.6

0.6

β

β

elcitrap

noisnepsus

0.5

0.5

0.4

0.4 Newtonian suspension Non-Newtonian suspension single phase Newtonian fluid

Newtonian suspension Non-Newtonian suspension

0.3

0.3 0

100

200

300

400

500

RI PT

d/D0 = 0.1

0.7

600

0

100

200

300

400

500

600

Re

SC

Re

(a)

(b)

M AN U

Figure 9: Experimental results for φ0 = 0.12 and Dk /D⊥ = 1 (a) βsuspension and (b) βparticle .

AC C

EP

TE

D

experimental observation of substantial difference in the separated flow area at the corner is believed to arise from migration. The possibility of elastic effects (not captured in these numerical results) in the guar solution at the corner cannot be ruled out. Figure 11 shows βparticle is significantly different than βsuspension for Dk /D⊥ = 3, which again points to the role of inertial migration. For φ0 = 0.05, βfluid agrees to within one per cent with the single phase fluid results, while the split of the fluid and suspension differ more for φ0 = 0.12. To address whether βparticle is affected by the migration of particles away from the wall, suspension flow with a narrower side channel was considered. Images from the experiments are shown in figure 12 for Dk /D⊥ = 5. We observed that no particles entered the side channel when d/D0 = 0.1 as shown in figure 14a; note that D⊥ /d ≈ 2 for the mean particle size and there are other particles that are smaller, so size exclusion does not fully explain the observation, even if excluded fluid volume along the wall due to the size of particle is considered. In case of smaller particle size (d/D0 = 0.027) in the same channel, we do see particles entering the side channel as shown in figure 14b. The potential value of using smaller proppant is highlighted here, as is the need to consider lift force and migration effects in addressing the expected entry of proppant into secondary fractures. Table 2 shows the experimental results for φ0 = 0.09 for smaller particles (d/D0 = 0.027). The difference between βsuspension and βparticle , is minor (0.098 vs 0.099, respectively).

4.3

Numerical analysis of upstream migration

To gain insight on the influence of upstream particle migration on the split of particles at the bifurcation, i.e. βparticle , we superimpose the streamlines obtained from the pure

13

RI PT

ACCEPTED MANUSCRIPT

(a)

(d)

TE

(c)

D

M AN U

SC

(b)

AC C

EP

Figure 10: Images of the suspension in the equal branch size (Dk /D⊥ = 1) T-bifurcation at φ0 = 0.12, d/D0 = 0.1: in (a) Newtonian fluid at Re = 50 and (b) 0.0042 weight fraction guar solution at Regen = 50. Computed streamlines for single phase fluid in T-bifurcation at (c) Re = 50 for a Newtonian fluid and (d) Regen = 50 for the power-law description of the 0.0042 weight fraction guar solution rheology.

β (single phase) βsuspension βparticle βfluid

0.982 0.981 0.99 0.98

Table 2: Average value of split for Dk /D⊥ = 5

14

ACCEPTED MANUSCRIPT

1.00 d/D0 = 0.1

Single phase Newtonain fluid

0.98

0.96

0.96

βsuspension βpartcile βfluid

β

β

0.98

0.94

0.94

0.92

Single phase Newtonian fluid

βsuspension βparticle βfluid

0.90

0.90 0

200

400

Re

600

800

1000

0

SC

0.92

d/D0 = 0.1

RI PT

1.00

200

400

600

800

1000

Re

M AN U

(a)

(b)

EP

TE

D

Figure 11: Experimental data on the various flow splits for suspensions in Newtonian liquid for Dk /D⊥ = 3 (a) φ0 = 0.05 and (b) φ0 = 0.12.

(b)

AC C

(a)

Figure 12: Images of suspension flow at φ0 = 0.05 and Re = 200 in channel of Dk /D⊥ = 5 (a) d/D0 = 0.1 and (b) d/D0 = 0.027. In (a), no particles enter the side channel for the duration of the experiment.

15

(b) Particles

(c) Combined

SC

(a) Streamlines

RI PT

ACCEPTED MANUSCRIPT

Figure 13: Superimposition of particle profile with streamlines

AC C

EP

TE

D

M AN U

fluid COMSOL calculations for the bifurcation geometry with the particle positions in the square cross-section obtained from lattice-Boltzmann method calculations of a pressuredriven suspension flow. In this analysis, we assume simply that the particles whose centers lie on the streamtube entering to a particular branch enters the corresponding branch. Based on this procedure, we obtain an estimate of the βparticle for each Re. This procedure thus estimates how cross-sectional spatial distribution affects the split to the daughter branches. The technique is illustrated in figure 13. For Dk /D⊥ = 1, figure 14 (a,b,c) shows the variation of the streamline regions entering each branch with increasing Re, where the shaded region corresponds to the straight branch. At smaller Re than shown here, the dividing streamline defines almost two equal regions (see figure 8). As Re increases, the region that enters the side branch increasingly covers the side walls of the square conduit, illustrating the three-dimensional nature of the flow. Ultimately, at high enough Re, the straight-branch region occupies the core of the cross-section, while the streamlines entering the side branch come from nearly the entire periphery of the conduit in the viscous layer. For a smaller side channel, Dk /D⊥ = 3, the picture differs as shown in figure 14 (d,e,f). The variation of the streamline regions with increasing Re is less pronounced than in the case of equal-sized daughter branches. The areas occupied by each of the regions remain similar, with the side-branch region mainly confined to the bottom of the channel. It is true that at the side walls the area occupied by the side-branch flow increases, but this is a relatively thin region and is not so extreme in changing the form of the regions as was observed for Dk /D⊥ = 1. Therefore the side branch only draws fluid from the region near the bottom wall of the inlet branch. Lattice-Boltzmann computed pressure-driven suspension flow which has traveled on average L0 /D0 = 58 is used to deduce the predicted concentration profile. This is the same flow distance as the experimental flow attains at the bifurcation. The cross-sectional positions of the particles at the junction entrance for φ0 = 0.12 at three different Re are 16

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

(b) Re = 447

(c) Re = 934

(d) Re = 130

TE

D

(a) Re = 130

(e) Re = 447

(f) Re = 934

AC C

EP

Figure 14: Streamlines of the inlet branch which go into each branch of the T-bifurcation for Dk /D⊥ = 1 and 3. The shaded region of the cross-section proceeds along the straight branch, while the white region goes to the side channel; the side channel extends down the page from the bottom of the square channel in this view.

17

(b)

(c)

SC

(a)

RI PT

ACCEPTED MANUSCRIPT

M AN U

Figure 15: Particle profile in straight channel near T-junction length (a)Re = 130 (b) Re = 447 (c)Re = 934

5

AC C

EP

TE

D

given in figure 15. The main observation is the transformation from a ring structure to an almost uniform distribution of particles to a structure with more particles migrating to the core of the tube at high enough Re. For the Dk /D⊥ = 1 channel, since the streamlines show that side-branch mainly samples the near-wall region of the inlet while the straight-region draws from the core, this implies that one expects to find more particles at the side-branch at low Re (βparticle < 1/2) while the opposite is expected at high enough Re. The numerical prediction of βparticle values based on the described superposition procedure are displayed in figure 16a, and agree with this trend and it shows qualitatively good agreement, with quantitative accuracy at Re > 400, when compared with the experimental results for Dk /D⊥ = 1. The streamline structure for Dk /D⊥ = 3 channels indicate that only the near-wall migration should have an effect on the particle phase split. The predicted βparticle values shown in figure 16b displays an increasing trend, i.e. the fraction of particles entering the side branch decreases with increasing Re. This is consistent with the expectation, as the cross-sectional particle distribution evolves from the ring structure towards center of the inlet tube, leaving the near-wall region increasingly depleted of particles as Re increases. On the other hand, for both Dk /D⊥ = 3 and Dk /D⊥ = 1, the numerical analysis underestimates the experimental βparticle values.

Summary and conclusions

We have studied the pure liquid and particle-laden flow in asymmetric T-bifurcations experimentally for Newtonian and non-Newtonian fluid along with numerical simulations using FEM to solve the Navier-Stokes equations for pure liquid complemented by discrete-particle simulations by the lattice-Boltzmann (LB) method for suspensions. For single-phase New-

18

ACCEPTED MANUSCRIPT

0.8

1.00

0.7 0.98

0.6

elcitrap

β

β

elcitrap

0.94 0.4

0.92 0.3

Experimental data Numerical prediction

Experimental data Numerical prediction

0.2

0.90 200

400

Re

600

800

1000

0

200

400

600

800

1000

Re

SC

0

RI PT

0.96 0.5

(a)

(b)

M AN U

Figure 16: Comparison of predictions by the superposition procedure and experiment for βparticle at φ0 = 0.12: (a) Dk /D⊥ = 1 and (b) Dk /D⊥ = 3.

AC C

EP

TE

D

tonian fluid, as Re increases, β(= Qk /Q0 ) increases while exhibiting weaker dependence on Re as Dk /D⊥ increases and it follows the form β = 1+e1−qRe where parameters are dependent on value of Dk /D⊥ . The flow behavior of guar polymer solution showing power law viscosity can be very well quantified using a generalized Reynolds number, Regen . The 3D complex structure of streamlines will not allow us to predict the quantitative behavior of flow based on 2D geometries. From the practical perspective, the desire to understand ‘tall’ or ‘2D’ fracture flows makes aspects of such 3D channel flows suspect, and care must be taken in carrying concepts here towards interpretation of behavior in such fractures. For suspension flows, the significant difference in particle flux ratio from the total material flux through each downstream branch is due to the inertial particle migration which was assessed by the LB simulations. The experiments of guar solution suspension shows stronger migration than a suspension in Newtonian fluid, resulting in higher βparticle at equivalent Re and Regen . For higher solid fraction (φ0 = 0.12 here), fluid streamlines are modified. This effect, when combined with particle migration makes an analysis of the flow based on some effective fluid with bulk material properties not only suspect but quite misleading: such approaches must account for the phase migration. One approach to this has been demonstrated here using numerical prediction for βparticle , based on the principle of superposition of particle profiles from LB calculation with flow streamlines obtained from COMSOL calculation. The effective excluded volume of particles from the near-wall region is much higher because of particle migration and this impacts on entry to small side channels, even when the channel is sufficiently wide for particles to enter.

19

ACCEPTED MANUSCRIPT

Acknowledgments

Nomenclature

RI PT

This work was supported in part by Halliburton.

β

Fraction of fluid in straight branch for single phase fluid

βfluid

Fraction of fluid in straight branch for particle laden flow

SC

βparticle Fraction of particles in straight branch for particle-laden flow

βsuspension Fraction of total suspension in straight branch for particle late flow Viscosity of fluid (P a.s)

φ0

Solid volume fraction at inlet

ρ

Density of fluid (kg/m3 )

ρp

Density of particles (kg/m3 )

d

Diameter of particle (m)

D0

Width of channel in inlet branch (m)

D⊥

Width of channel in side branch (m)

Dk

Width of channel in straight branch (m)

j

Total particle flux in inlet channel

jk

Particle flux in straight branch

L0

Length of channel from inlet to midpoint of T-junction (m)

L⊥

Length of channel from midpoint of T-junction to side branch outlet(m)

N

D

TE

EP

AC C

Lk

M AN U

η

Length of channel from midpoint of T-junction to straight branch outlet(m) Total number of particles in inlet channel

Nk

Number of particle in straight branch

Q⊥

Volumetric flow rate in side branch (m3 /s)

All units are given in SI

20

ACCEPTED MANUSCRIPT

Volumetric flow rate in straight branch (m3 /s)

Re

Reynolds number based on inlet velocity and inlet channel width

U

Average velocity of fluid at inlet (m/s)

Up

Particle velocity in the flow direction (m/s)

References

RI PT

Qk

SC

Ahmed, G.M.Y. & Singh, A. 2011 Numerical simulation of particle migration in asymmetric bifurcation channel. J. Nonewton Fluid Mech. 166, 42–51. Berkowitz, B. & Adler, PM. 1998 Stereological analysis of fracture network structure in geological formations. J. Geophys. Res. 103 (B7), 15339–60.

M AN U

Bhagat, A. A. S., Kuntaegowdanahalli, S. S. & Papautsky, I. 2008 Enhanced particle filtration in straight microchannels using shear-modulated inertial migration. Phys. Fluids 20 (101702), 1–5. Blunt, M.J. 2001 Flow in porous media pore-network models and multiphase flow. Curr. Opin. Colloid Interface Sci. 6, 197–207. Bohloli, B. & Pater, C.J. de 2006 Experimental study on hydraulic fracturing of soft rocks: Influence of fluid rheology and confining stress. J. Pet. Sci. Eng. 53, 1–12.

TE

D

Boyer, F., Pouliquen, Q. & Guazzelli, E. 2011 Dense suspension in rotating-rod flows: normal stresss and particle migration. J. Fluid Mech. 686, 5–25. Bugliarello, G. & Hsiao, G. C.C 1964 Phase separation in suspensions flowing through bifurcations: A simplified hemodynamic model. Sci. 143 (3605), 469–471.

EP

Chun, B. & Ladd, A.C.J. 2006 Inertial migration of neutrally buoyant particles in a square duct: An investigation of multiple equilibrium positions. Phys. Fluids 18 (031704), 1–5.

AC C

COMSOL 2015 Comsol multiphysics 5.2. https://www.comsol.com Stockholm Sweden. Delplace, G. & Leuliet, J.C. 1995 Generalized reynolds number for the flow of power law fluids in cylindrical ducts of arbitrary cross-section. Chem. Eng. J. 56, 33–37. Doyeux, V., Podgorski, T., Peponas, S., Ismail, M. & G., Coupier 2011 Spheres in the vicinity of a bifurcation: elucidating the zweifach-fung effect. J. Fluid Mech. 674, 359–388.

21

ACCEPTED MANUSCRIPT

Enden, G. & Popel, AS. 1994 A numerical study of plasma skimming in small vascular bifurcations. J. Biomech. Eng. 116 (1), 79–88.

RI PT

Garagash, D., Lecampion, B. & Desroches, J. 2015 Pressure-driven suspension flow near jamming. Phys. Rev. Lett. 114 (088301), 1–6. Garagash, D. I. & Lecampion, B. 2014 Confined flow of suspension modeled by a frictional rheology. J. Fluid Mech. 759, 197–235. Haddadi, H. & Morris, J. F. 2014 Microstructure and rheology of finite inertia neutrally buoyant suspensions. J. Fluid Mech. 749, 431–459.

SC

Haddadi, H., Zadeh, S.S., Connington, K. & Morris, J. F. 2014 Suspension flow past a cylinder: particle interactions with recirculating wakes. J. Fluid Mech. 760, 1–12.

M AN U

Hampton, R.E., Mammoli, A.A., Graham, A.L., Tetlow, W. & Altobelli, S.A. 1997 Migration of particles undergoing pressure-driven flow in a circular conduit. J. Rheol. p. 621. Han, M., Kim, C., M., Kim & Lee, S. 1999 Particle migration in tube flow of suspesnsion. J. Rheol. 43 (197), 1157–1174. Iwanami, S. & Suu, T. 1969 Study on flow characteritics in right-angled pipe fittings. Bull. JSME (2nd Report, On Case of Slurries in Hold up Flow) 12 (53), 1051–1061.

D

Iwanami, S., Suu, T. & Kato, H. 1969 Study on flow characteritics in right-angled pipe fittings. Bull. JSME (1st report, On Case of Water Flow) 12 (53), 1041–1050.

TE

Koh, C. J., Hookham, P. & Leal, L.G 1994 An experimental investigation of concentrated suspension flows in a rectangular channel. J. Fluid Mech. 266, 1–32.

EP

Kulkarni, P.M. & Morris, J.F. 2008 Pair-sphere trajectories in finite-reynolds-number shear flow. J. Fluid Mech. 596, 413–435. Kwicklis, ED. & Hearly, R.W. 1993 Numerical investigation of steady liquid water flow in a variably saturated fractured network. Adv. Water Resour. 29, 4091–102.

AC C

Ladd, A.J. 1994a Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 2. numerical results. J. Fluid Mech. 271 (311-339). Ladd, A. J.C. 1994b Numerical simulations of particulate suspensions via a discretized boltzmann equation. part 1. theoretical foundation. J. Fluid Mech. 21, 285–309. Liepsch, D. & Moravec, S. 1982 Measurement and calculations of laminar flow in a ninety degree bifurcation. J. Biomech. 15 (7), 473–85.

22

ACCEPTED MANUSCRIPT

Matas, J.P., Morris, J. F. & Guazzelli, E. 2004 Inertial migration of rigid spherical particles in poiseuille flow. J. Fluid Mech. 515, 171–195.

RI PT

Matas, J.P., Morris, J. F. & Guazzelli, E. 2009 Lateral fore on a rigid sphere in large-inertia laminar pipe flow. J. Fluid Mech. 621, 59–67. Metzner, A.B. & Reed, J.C. 1955 Flow of non-newtonian fluids—correlation of the laminar, transition, and turbulent-flow regions. AICHE J. 1 (4), 434–440.

SC

Miranda, A. I. P., Oliveira, P. J. & Pinho, F. T. 2008 Steady and unsteady laminar flows of newtonian and generalized newtonian fluids in a planar t-junction. Int. J. Numer. Methods Fluids 57, 295–328.

M AN U

Miura, K., Itano, T. & Sugihara-Seki, M. 2014 Inertial migration of neutrally buoyant spheres in a pressure-driven flow through square channels. J. Fluid Mech. 749, 320– 330. Mo, H., Bai, M., Lin, D. & Roegiers, J.C. 1998 Study of flow and transport in fracture network using percolation theory. Appl. Math. Modell. 22, 277–291. Morris, J.F. & Boulay, F. 1999 Curvilinear flows of noncolloidal suspensions: The role of normal stresses. Soc. Rheol. 48 (5), 1213–1237. Neary, V.S. & Sortiropoulos, F. 1996 numerical investigation of laminar flows through 90-degree diversions of rectangular cross-section. Comput. Fluid 25 (2), 95–118.

D

Nott, P. R. & Brady, J. F. 1994 Pressure-driven flow of suspensions: simulation and theory. J. Fluid Mech. 275, 157–199.

TE

Pedley, T.J. 2008 The fluid mechanics of large blood vessels. Cambridge University Press.

EP

Piggot, AR. 1997 Fractal relations for the diameter and trace length of disc-shaped fractures. J. Geophys. Res. 102 (B8), 18121–5. Pries, A.R., Secomb, P., Gaehtgens, P. & Gross, J.F 1990 Blood flow in microvascular networks. Circ. Res. 67, 826–834.

AC C

Roberts, B. W. & Olbricht, W.L. 2006 Flow-induced particulate separations. AICHE J. 49 (11), 2842–2849. Segre’, G. & Silberberg, A. 1961 Radial particle displacements in poiseuille flow of suspensions. Nature 189, 109–210. Shen, HH., Satake, M., Mehrabadi, M., Chang, C.S. & Campbell, C.S., ed. 1992 Advances in micromechanics of granular materials. Elsevier.

23

ACCEPTED MANUSCRIPT

Sierra-Espinosa, F.Z., Bates, C.J. & T., O’Doherty 2000 Turbulent flow in a 90 ◦ pipe junction: Part 2: reverse flow at the branch exit. Comput. Fluid 29, 215–233.

RI PT

Warburton, PM. 1989 Stereological interpretation of joint trace data: influence of joint shape and implication for geological surveys. J. Rock Mech. Min. Sci. Geomech. Abstr. 17, 305–316. Warpinski, N.R., Mayerhofer, M.J., Vincet, M.C., Cipolla, C.L. & Lolon, E.P. 2009 Stimulating unconventional reservoirs: maximizing network growth while optimizing fracture conductivity. J. Can. Pet. Technol. 48 (10), 39–51.

SC

Xi, C. & Shapley, N. C. 2008 Flow of a concentrated suspension through an axisymmetric bifurcation. J. Rheol. 52 (2), 625–647.

AC C

EP

TE

D

M AN U

Yan, Z-Y., Acrivos, A. & Weinbaum, S. 1991 A three-dimensional analysis of plasma skimming at microvascular bifurcations. Microvasc. Res. 42 (1), 17–38.

24

ACCEPTED MANUSCRIPT

Highlights:

•

AC C

EP

TE D

M AN U

SC

•

The suspension behavior is significantly different from single-phase fluid flow in fractures. The relative size of particle to width of a channel plays an important role in particle distribution in fractures. The ratio of width of side fracture to the main fracture is critical in quantification of particle and fluid flow rates in fracturing networks.

RI PT

•