UBC Faculty Research and Publications

Slurry flow, gravitational settling, and a proppant transport model for hydraulic fractures Dontsov, Egor; Peirce, Anthony Aug 13, 2014

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


52383-Dontsov_Peirce_Slurry_Flow_unformatted.pdf [ 1.64MB ]
JSON: 52383-1.0079342.json
JSON-LD: 52383-1.0079342-ld.json
RDF/XML (Pretty): 52383-1.0079342-rdf.xml
RDF/JSON: 52383-1.0079342-rdf.json
Turtle: 52383-1.0079342-turtle.txt
N-Triples: 52383-1.0079342-rdf-ntriples.txt
Original Record: 52383-1.0079342-source.json
Full Text

Full Text

Slurry flow, gravitational settling, and a proppant transport modelfor hydraulic fracturesE.V. Dontsov and A.P. PeirceAugust 8, 2014AbstractThe goal of this study is to analyze the steady flow of a Newtonian fluid mixed with spherical particles ina channel based on a continuum model, where the constitutive behaviour of the slurry is approximatedby an empirical formula. In order to account for the gravitational settling of particles, two-dimensionalflow needs to be considered as the pressure gradient and gravity may not always be collinear. It is shownthat the problem under consideration features a boundary layer, whose size is on the order of the particleradius. The expressions for both the outer (i.e. outside the boundary layer) and inner (i.e. within theboundary layer) solutions are obtained in terms of the particle concentration, particle velocity, and fluidvelocity. Unfortunately, these solutions require numerical solution of an integral equation, depend on theratio between the pressure gradient and the gravity force, and the orientation of the pressure gradientrelative to the gravity. Consequently, the development of a proppant transport model for hydraulicfracturing based on these results is not practical. For this reason, an approximate solution is introduced,where the effect of gravity is accounted for in an approximate fashion, reducing the complexity of theslurry flow solution. To validate the use of this approximation, the error is estimated for different regimesof flow. The approximate solution is then used to calculate the expressions for the slurry flux and theproppant flux, which are the basis for a model that can be used to account for proppant transport withgravitational settling in a fully coupled hydraulic fracturing simulator.1 IntroductionThe flow of a slurry is a problem that is relevant to many natural processes, such as mud flows or landslides,as well as to industrial applications, such as the flow of a cement slurry. This study, however, addressesthe problem of slurry flow in the context of hydraulic fracturing that is relevant to oil and gas reservoirstimulation [1]. Typically, proppant is used to prevent the fracture from closing once the well is depressurized.In this case, modelling the fracture propagation driven just by a viscous fluid is not sufficient, since theproppant, blended with the fracturing fluid, alters the properties of the latter. In this situation, it isnecessary to consider the effects associated with the flow of a slurry as well as the proppant transport inhydraulic fractures. As mentioned in [2], for the purpose of hydraulic fracturing, the slurry is typicallymodelled as a Newtonian fluid, whose viscosity dependence on proppant content is calculated based on anempirical formula. Moreover, the particle distribution across the fracture is assumed to be uniform and onlythe slip velocity due to gravity is considered. While all the aforementioned simplifications could potentiallylead to sufficiently accurate results in some cases, it is nearly impossible to judge the accuracy withouthaving a higher-level, more accurate model. To obtain a more accurate model, this study aims to analyzethe steady flow of a slurry in a channel based on a recently introduced constitutive model [3], and to establisha framework necessary to formulate the problem of hydraulic fracturing by a slurry, which accounts for theproppant transport, gravitational settling, as well as formation and growth of packed regions.To simplify the analysis, it is assumed that the proppant particles are spherical and that they all havethe same size. In order to describe the motion of a viscous fluid with monodisperse spherical particles, it isnecessary to distinguish between different flow regimes. Following [4], there are two primary dimensionless1groups that control the type of motion, namely the Reynolds and Peclet numbers, that are defined asRe = ρfa2γ˙µf , P e =6piµfγ˙a3kBT. (1)Here ρf and µf denote the mass density and the viscosity of the fluid respectively, a is the particle radius, γ˙ isthe magnitude of the shear rate, kB is Boltzmann’s constant, while T is the absolute temperature. Note thatone may replace ρf in (1a) by ρp+12ρf (ρp is the mass density of particles), in which case the Reynolds number,Re, would reflect the ratio between inertia forces (accounting for the added fluid mass) and viscous forces.On the other hand, the Peclet number, Pe, reflects the effects associated with chaotic thermal motion andcontrols the transition from Brownian to non-Brownian motion. By dealing with sufficiently large particlesand a sufficiently viscous fluid, this study is focused on the case wherePe→∞, Re→ 0,which corresponds to non-Brownian motion of particles in highly viscous fluids. This regime has been studiedfrom both theoretical and experimental prospectives [5, 6, 7, 8, 9], where one of the main objectives was thestudy of particle migration in shear flow. In simple words, the particles try to avoid high shearing and tomove towards the regions with smaller shear rates. Note that this phenomenon is different from the Segre-Silberberg effect [10], since the latter is observed in the inertial regime, where the viscosity of the carryingfluid is small. One very important parameter in the analysis of particle migration is the characteristic timescale for reaching steady-state flow. As noted in [6], some earlier experimental studies failed to show theeffect of particle migration due to the fact that the flow was not in the steady-state regime. As also discussedin [6], the characteristic length, required to establish a steady flow for dense suspensions can be estimatedasLw ∼(wa)2, (2)where w is the characteristic width of the channel and a is the particle radius.A great deal of effort has been devoted to developing an accurate constitutive model for a mixture of aNewtonian fluid with spherical particles. It started with [11], who introduced a first-order correction to theviscosity of dilute suspensions. Later, a second-order correction was made by [12]. The most challengingproblem, however, was to establish a constitutive model that captures the behaviour of dense suspensionswith high particle concentrations where the interaction between particles plays a crucial role. For instance, acomparison of various approaches for modelling dense suspensions can be found in [4]. Another constitutivemodel has been recently proposed by [3], where experimental observations were used to establish empiricalrelations between shear and normal stresses versus particle concentration. This, the latest model, appearsto be the most accurate to date, and, for this reason, is chosen for the analysis in this paper.The paper is organized in the following way. First, the governing equations for the slurry flow, based onthe empirical constitutive model are in formulated in Section 2. Then, the resulting equations are used tosolve the problem of slurry flow in a channel, which is described in Section 3. Recognizing the complexityof the solution, Section 3 also introduces an approximate solution and estimates the errors caused by thisapproximation. Finally, Section 4 utilizes the approximate solution of the channel problem to find the massbalance laws for the slurry and the particles, which feature the average flow of the slurry and the proppantthrough the channel.2 Governing equations for motion of fluid with spherical particlesTo formulate the governing equations for the motion of the slurry (i.e. the mixture of viscous fluid andparticles), this study follows an approach described in [8]. The main difference comes from the constitutivebehaviour, which is based on recent experimental results by [3]. The balances of linear momentum and massfor the particles are taken in the following formφρp(∂vp∂t + vp ·∇vp)= −∇·(Qpp) + 2∇·(µp∇s0vp)+ φρpg + fpf,∂φρp∂t +∇·(φρpvp) = 0. (3)2Here φ is the volume fraction of the particles, ρp is the particle mass density, vp is the velocity associatedwith the macroscopic movement of the particles, pp is the pressure due to the particle collisions (definedas the time-averaged momentum transfer per unit area), µp is the effective viscosity associated with theparticles, ∇s0vp = 12 (∇vp+∇Tvp−23∇·vp I) is the deviatoric strain rate tensor, g is the gravity force per unitmass, and fpf is the interaction force between the viscous fluid and the particles. As follows from [8], thesecond-order tensor Q describes the anisotropy of the normal stresses and can be represented asQ =3∑i=1λiei ⊗ ei,where λi are dimensionless constants and ei are unit vectors in the direction of the flow (i= 1), gradient(i= 2), and vorticity (i= 3). Note that one of the λi is always equal to unity since the stress magnitude isreflected in pp. The governing equations for the fluid can be written in a similar fashion as(1−φ)ρf(∂vf∂t + vf ·∇vf)= ∇·σf + (1−φ)ρfg − fpf,∂(1−φ)ρf∂t +∇·((1−φ)ρfvf)= 0, (4)where ρf, pf, and vf denote the mass density, pressure, and macroscopic velocity associated with the fluidrespectively, σf = −pfI +2µf∇s0vf is the fluid stress tensor, while µf is the fluid viscosity. Note that the massdensity of the mixture is ρ = (1−φ)ρf + φρp, the total pressure (i.e. total normal force per unit area actingon a wall or imaginary plane, neglecting viscosity) is p = pf + ppn·Q·n, while the total shear force also hascontributions from both the particles and the fluid. Note that it is very important to understand that theparticle pressure is defined as an average momentum transfer (due to collisions) per unit time per unit area– similar to gases. So, one may regard the slurry as a mixture of two gases with partial pressures pp and pf.This is consistent with the way pp is defined in the experiments [3], since the total force per unit area actingon the top plate from the top (see Fig. 1b in [3]) is pf + pp, and this should match the total force per unitarea coming from the slurry, so the total pressure of the slurry is pf + pp. Similar governing equations werealso used by [13] in the context of sediment dynamics.On considering both the viscous fluid and the particles to be incompressible, the remaining task is tospecify pp, µp, fpf and λi (i= 1, 2, 3). By generalizing the constitutive model, that is proposed in [3], onemay takepp = µfA(φ)−2√2∇s0vp :∇s0vp,µp = µf[52φmA(φ)−1+(µ1+µ2 − µ11+I0A(φ)−2)A(φ)−2], A(φ) = φm/φ− 1, (5)where the maximum volume fraction φm=0.585, and the constants µ1 =0.32, µ2 =0.7 and I0 =5×10−3 areempirical and based on experimental observations. Note that A(φ)n is a short-hand notation for(A(φ))n.Since the flow regime under consideration assumes small Reynolds numbers, the viscous part of the interactionforce can be deduced from Stokes’ law, while the second part of the force comes from the effect of buoyancy,so thatfpf = η(φ)(1−φ)(vf − vp) + φ∇·σf, η(φ) = 9µfφ2a2(1−φ)α¯ , (6)where an additional correction term (1−φ)α¯ in the expression for η(φ) accounts for interaction between theparticles at high values of φ (see e.g. [8]). Note, that these authors defined the slip velocity as vf−〈v〉=vf−φvf−(1−φ)vp =(1−φ)(vf−vp), so that their α is related to α¯ by α¯ = α−1. This parameter was chosenα=4 in [8], while according to [14] and [15], α = 5.1 provides a better fit. For this reason, α¯ = 4.1 is chosenfor further computations. Finally, the values of λi (i=1, 2, 3) are estimated in [16] asλ1 = 1.05, λ2 = 1, λ3 = 0.65.As also noted in [16], λi satisfy λ1 > λ2 > λ3, λ1 ' λ2, and λ3 ' 12 . As will be shown shortly, these valuesof λi (i = 1, 2, 3) are not essential for the problem under consideration.33 Slurry flow in a fracture3.1 Problem formulation and assumptionsThis section aims to analyze the steady motion of a slurry in a fracture, which for our purposes is modelledas a channel of width w. Since the gravity and the pressure gradient may act in different directions and theproblem is nonlinear (i.e. superposition does not apply), it is necessary to consider two-dimensional flowinside the fracture. To this end, we introduce a coordinate system (x, y, z), where x is the horizontal axisalong the fracture, z is the vertical axis, while y is the coordinate across the fracture, as shown in Fig. 1.Assuming a steady flow, the governing equations (3a) and (4a) can be reduced to2axyz gezxzwvpFigure 1: Schematics of the hydraulic fracture (left) and the slurry flow inside it (right).∇·(Qpp) = ∂∂y(µp(φ)∂vp∂y)+ η(φ)∆v − φ(ρp−ρf)gez,∇p˜f = µf ∂2vp∂y2 + µf ∂2∆v∂y2 − η(φ)∆v + 〈φ〉(ρp−ρf)gez, (7)where ∇ = (∂/∂x, ∂/∂z), and vectors vp and ∆v = vf − vp also have only x and z components, 〈φ〉 denotesaverage particle concentration across the channel, whilep˜f = pf + (ρfg + 〈φ〉(ρp−ρf)g)z,is the pressure that accounts for hydrostatic pressure. The particle pressure and effective shear viscosity ofthe particles are given respectively bypp = µfA(φ)−2∣∣∣∂vp∂y∣∣∣,µp = µf(52φmA(φ) + µ1 +µ2 − µ11 + I0A(φ)−2)A(φ)−2, A(φ)=φm/φ− 1. (8)Here p˜f = p˜f(x, z) and pp = pp(x, z) due to the steady-state property of the flow.Assumptions. To account for realistic values of the parameters and to simplify the problem, severalassumptions are made, in particular:1. For a continuum constitutive model to apply, it is implicitly assumed that 2a w.2. The particle pressure gradient is neglected. To justify this assumption, note that the characteristiclength of the problem, L (i.e. the length of the fracture), can be safely assumed to be much bigger thanthe width of the fracture, w. For instance, in typical hydraulic fracturing geometries, w/L = O(10−3)or less. In this case, one can use (7b) and (8a) to obtain estimates for pf and pp as∇p˜f = O(µfvfw−2), ∇pp = O(µfvp(Lw)−1),4which leads to |∇pp/∇pf| = O(w/L)  1, provided that vf = O(vp). This estimate shows, that thepressure gradient associated with the particles has a negligible contribution to the momentum change ofthe mixture. Note that this is true only for φ 6= φm. As soon as particles become packed, i.e. φ = φm,they may sustain notable pressures and the particle pressure may not be neglected. In addition, itshould be noted that ignoring the particle pressure gradient leads to the fact that the values of λ1 andλ3 are no longer important, because the variation of particle pressure becomes relevant only across thechannel (see (8)), i.e. in the direction that corresponds to λ2.3. There is no gravity force in the y direction, i.e. across the channel. The geometry of the hydraulicfracture is determined predominantly by confining stresses. For deep fractures, the vertical stress istypically higher than the horizontal stress, in which case the fracture propagates in a vertical plane,and thus the projection of gravity force on y axis is zero.3.2 Steady solutionTo obtain a steady solution, it is noted that the problem under consideration (7) features a length scalea w, which causes a presence of a boundary layer with width O(a). Let’s first focus on the outer solution,i.e. the one that is valid away from the channel walls. In this situation, the term with the second derivativeof ∆v in (7) can be neglected. Indeed, when the y coordinate scales with w, one can estimate∂∂y(µf ∂∆v∂y)= O(µf∆vw2), η(φ)∆v = O(φµf∆va2),so that the term with the second derivative of ∆v can be neglected as soon as a2/(φw2) 1. In situationswhen φ=O(a2/w2) 1, the viscosity associated with proppant is µp(φ)=O(φ) 1, which makes proppantcontribution to the slurry flow negligible and leads to a well-known parabolic velocity profile. The factthat the term with the derivative of ∆v can be neglected should not be confused with the assumption ofsmallness of |∆v/vp|. For relatively small particle concentrations indeed |∆v/vp| = O(a2/w2)  1, whilehigher concentrations (φ≈ φm) may lead to |∆v/vp| = O(µ1φ2a2/(w2(φm−φ)2)) = O(1), or even vp = 0for φ=φm. As was mentioned before, the term with the derivative of ∆v can be neglected away from theboundary, but at the distances O(a) from the boundary (i.e. in the boundary layer where the y coordinatescales with a), this term may have a significant contribution. As a result of the presence of the boundarylayer, neither of vp nor ∆v (calculated for the outer solution) should satisfy a no slip boundary conditionat channel walls. Instead, ∆v should approach some constant, while vp should approach another constant,whose magnitude is O(∆v). Naturally, these constants should arise from the solution of the boundary layerproblem.By neglecting the particle pressure gradient and the second derivative of ∆v in (7), integrating the sumof the resultant equations, and using (8), one can write∇p˜f(y−y0) + (ρp−ρf)gez∫ yy0(φ−〈φ〉) dy¯ = (µf+µp(φ))∂vp∂y − τ(y0), y0 < y 612w,[∇p˜f + (φm−〈φ〉)(ρp−ρf)gez]y = τ, |τ | 6 µ1pp, 0 6 y 6 y0, (9)where τ denotes shear stress. Here only the solution for y > 0 is considered due to the symmetry consid-erations. As indicated in the solution (9), there is a region with “no failure” (0 6 y 6 y0), i.e. where theparticles form a rigid cluster, and a region with shear motion or failure (y0 < y 6 12w). To find the particlepressure, one needs to evaluate (9b) at y= y0, add the result to (9a) evaluated at y=w/2 and use (8a) toobtainpp = wµf|∇p˜f|2(µf+µp(φw))A(φw)2, (10)where φw = φ|y=w/2 is the particle concentration at the wall. Since |τ | = µ1pp at the plug boundary,equations (10) and (9b) can be used to calculate the plug size asy0 =w2µ1[1 + (φm/〈φ〉−1)2G2ρ + 2(φm/〈φ〉−1)Gρ cosψ]−1/2(1+µp(φw)/µf)A(φw)2, (11)5whereGρ =(ρp−ρf)g〈φ〉|∇p˜f|, cosψ = |∇p˜f|−1 ∂p˜f∂z , (12)are two dimensionless parameters that represent the ratio between the gravitational force and the pressuregradient, and the angle between the pressure gradient and the vertical z axis. respectively. To find a solutionfor the particle concentration φ, it is useful to rewrite (7a) using (8a) and (10) as:[(2yw)2+ 4G2ρ[∫ y/w0(φ/〈φ〉−1) ds]2+ 8yw Gρ cosψ∫ y/w0(φ/〈φ〉−1) ds]1/2= (µf+µp(φ))A(φ)2(µf+µp(φw))A(φw)2, y0 < y 6 12w.(13)Given the values of Gρ, cosψ, and φw, φ can be determined by solving (13) numerically via Newton’s method.In this case the solution can be written in the form φ(s, φw, Gρ, cosψ) or φ(s, 〈φ〉, Gρ, cosψ), where s = y/w.It is interesting to comment on the existence of a solution of (13). Unfortunately, (13) does not alwayshave a solution. A thorough analysis of the existence of the solution is beyond the scope of this study, forthis reason, let’s focus on the simple case when cosψ = −1. In this situation, the left hand side of (13) isproportional to the absolute value of the function h(y) = 1−wGρ/y∫ y/w0 (φ/〈φ〉−1) ds. This function h(y) isequal to h(y0) = 1−Gρ(φm/〈φ〉−1) at the plug boundary, while h(w/2) = 1 at the wall. If h(y0) < 0, thenat some point y0 < y∗ < w/2 this function vanishes (provided that it is continuous). However, the righthand side of (13) cannot vanish, since its minimum value is µ1µf/[(µf +µp(φw)A(φw)2]> 0. In this case,the necessary condition for the existence of a solution is 1−Gρ(φm/〈φ〉−1)> 0. This condition bounds thepressure gradient as−∂p˜f∂z > (ρp−ρf )g(φm−〈φ〉). (14)Note that ∂p˜f/∂x = 0 and −∂p˜f/∂z > 0 since cosψ = −1, in which case the range over which the steadysolution does not exist (according to (14)) is relatively narrow. It is also possible to explain why there is nosteady solution. First, if −∂p˜f/∂z > 0, the slurry should flow in the positive z direction since the averageof the corresponding gravity term in the sum of (7) is zero. Since the shear stress is zero at the centre ofthe channel, the particles tend to form a plug. If the plug is formed, then, if (14) is not satisfied, the plugstarts to sink as the pressure gradient is not sufficient to sustain it. If the plug sinks (i.e. flows in thenegative z direction), and the average velocity is positive, then, since the shear stress should be continuous,there is a point at which the shear stress is zero (as its sign is different at the wall and the boundary ofthe plug). If the shear stress is zero, then the particle pressure is zero as well (or there is another plug,which, if exists, would sink too and cause a similar problem). Since there is some particle pressure at thecentral plug (otherwise, it would not form) and this pressure vanishes somewhere outside of it, then thereis a particle pressure gradient, which moves particles away from the plug. As soon as particles move awayand the concentration in the plug reduces, the particle pressure at the centre becomes zero (since the shearstress is zero due to symmetry), and this generates an opposite particle pressure gradient, which tends toform a plug again. This cyclic process repeats, so no steady solution exists in this case.To find the outer solution (i.e. outside of the boundary layer) for the particle velocity, (9a) can beintegrated asvp = vp|y=w/2 + vp0 (15)= vp|y=w/2 −w2∇p˜fµf∫ 1/2sˆs¯ ds¯1+µp(φ(s¯))/µf −w2(ρp−ρf)gezµf∫ 1/2sˆ11+µp(φ(s¯))/µf∫ s¯0(φ(s)−〈φ〉) ds ds¯,where sˆ = max{y0/w, y/w}, vp0 is the solution that satisfies the no slip condition at the boundary, whilevp|y=w/2 is the integration constant that comes from the solution of the boundary layer problem. The slipvelocity for the outer solution can be found from (7a) and using (9a) as∆v = − ∇p˜fη(φ)∂∂s( µp(φ)sµf + µp(φ))+ (ρp−ρf)gezη(φ)∂∂s( µf∫ s0 φds¯µf + µp(φ) +µp(φ)〈φ〉sµf + µp(φ)), (16)6where s = y/w is the scaled y coordinate.To find the inner or boundary layer solution, it is useful to represent the “full” solution as the sum ofthe outer and the inner, so that φF = φ+ δφ, ∆vF = ∆v + δ∆v, and vp,F = vp0 + δvp. By substituting theabove expressions into (7) and (8), and using a Taylor’s series expansion, one finds that∂µp∂φ∣∣∣∣w∂vp∂y∣∣∣∣w∂δφ∂y + µp(φw)∂2δvp∂y2 + η(φw)δ∆v = 0,µf ∂2δvp∂y2 + µf ∂2δ∆v∂y2 − η(φw)δ∆v = 0, (17)∣∣∣∣∂vp∂y∣∣∣w∣∣∣∣−2(∂δvp∂y ·∂vp∂y∣∣∣w)= 2A(φw)∂A∂φ∣∣∣∣wδφ,where “|w” means evaluation of the outer solution at the wall. Note that since y = O(a) and |∆v/vp| =O(a2/w2) for the inner solution (y = O(w) for the outer solution), it can be concluded that δφ = O(a/w) 1,δvp = O(∆v) vp, which validates the use of a Taylor series expansion. At the same time, δ∆v = O(∆v),but no expansion has been used with ∆v. The solution of the above system of equations (17), that accountsfor both, the no slip boundary conditions and the far-field behaviour isδ∆v = −V1e−γ1y − V2e−γ2y,δvp = C1V1(1− e−γ1y) + C2V2(1− e−γ2y), (18)δφ = C3e−γ1y,whereV1 =∣∣∣∣∂vp∂y∣∣∣w∣∣∣∣−2(∆v|w ·∂vp∂y∣∣∣w)∂vp∂y∣∣∣w, V2 = ∆v|w − V1,γ1 =( η(φw)µf(1−B(φw)))1/2, γ2 =(η(φw)µf(1+ µfµp(φw)))1/2,C1 = −B(φw), C2 = −µfµf + µp(φw), (19)C3 = B(φw)( η(φw)µf(1−B(φw)))1/2∣∣∣∣∂vp∂y∣∣∣w∣∣∣∣−2(∆v|w ·∂vp∂y∣∣∣w)A(φw)φ2w2φm,and the auxiliary function B(φ) isB(φ) = 2A(φ)[2A(φ) + 52φm +2(µ2−µ1)I0A(φ)(I0 +A(φ)2)2]−1. (20)Note that ∆vF = ∆v + δ∆v satisfies the no slip boundary condition since V1 + V2 = ∆v|w. The “full”solution for the particle velocity vp,F = vp0 + δvp satisfies the no slip boundary condition as well. To findvpx|y=w/2 featured in (15), one needs to take a limit in (18b) to findvp|y=w/2 = C1V1 + C2V2. (21)Since C1(φw = 0) = C2(φw = 0) = −1 and C1(φw = φm) = C2(φw = φm) = 0, equation (34) implies thatvf|y=w/2 = 0 for zero proppant concentration and vp|y=w/2 = 0 for the maximum proppant concentration.Indeed, if there are no particles, the fluid velocity (outer solution) should satisfy a no slip boundary condition,at the same time, the maximum particle concentration corresponds to a rigid plug, which does not move,so that the particle velocity (outer solution) should satisfy a no slip condition at the boundary. It isalso important to note that the values of γ1 and γ2 are always O(1/a), and in the limit of small particleconcentrations they reach finite values γ1≈√18/5 a−1, and γ2≈√9/5 a−1.To visualize the solution, Fig. 2 shows the velocity profile vp0 (calculated numerically), normalized byv∗=(12µf)−1w2|∇p˜f| and proppant concentration variation (calculated numerically), normalized by φm, for70 0.5 1 −0.50 0.5−1.5−1−0.50y/wvpx,0/v∗vp z,0/v∗0 0.5 1 −0.50 0.5−1.5−1−0.50y/wvpx,0/v∗vp z,0/v∗0 0.5 1 −0.5 00.5−2−10y/wvpx,0/v∗vp z,0/v∗0 0.5φ/φmy/w0 0.5φ/φmy/w0 0.5φ/φmy/wcos = 0cos =12cos = 1G⇢=0G⇢=1G⇢=4G⇢=10G⇢=100G⇢=0G⇢=1G⇢=4G⇢=10G⇢=100G⇢=0G⇢=1G⇢=4G⇢=10G⇢=100G⇢=0G⇢=1G⇢=4G⇢=10G⇢=100G⇢=0G⇢=1G⇢=4G⇢=10G⇢=100G⇢=0G⇢=1G⇢=4G⇢=10G⇢=100Figure 2: Normalized particle velocity profile and normalized particle volume fraction versus y/w for Gρ ={0, 1, 4, 10, 100}, cosψ = 0 (top row), cosψ = 12 (middle row), cosψ = 1 (bottom row) and 〈φ〉 = 0.2× φm.8different values of the parameter Gρ = {0, 1, 4, 10, 100}, and for different values of cosψ = {0, 12 , 1}, whilethe average proppant concentration is kept fixed φ¯ = 〈φ〉/φm = 0.2. As can be seen from the figure, thepresence of the gravitational force leads to redistribution of the particles across the channel. Moreover, thevalue of cosψ has a small impact on the particle distribution. In the limiting case when Gρ1, the particlesare distributed uniformly everywhere except for the narrow region at the centre (which corresponds to theplug). For cosψ = 1 and Gρ = 100, the velocity profile resembles a wedge. This strange velocity profile canbe explained by the fact that the forcing term for the sum of equations (7) behaves like a delta functionfor Gρ 1, i.e. ∂p˜f/∂z(1+Gρ(φ/〈φ〉−1))→ ∂p˜f/∂z w δ(y), where δ(y) denotes Dirac delta function. Thisdetla-function-like pressure gradient distribution leads to the wedge-like velocity profile.3.3 Approximate steady solutionNumerical solution of (13) together with equations (11), (15), (16) and (34) give a complete solution for theproblem. Unfortunately, this solution relies on the numerical evaluation of the function φ(s, 〈φ〉, Gρ, cosψ),which depends on many parameters and thus is hard to tabulate for the purpose of implementing into ahydraulic fracturing simulator. Moreover, as noted in the previous section, a steady solution does not exist forall ranges of parameters. For these reasons, it is more practical to introduce an approximate solution, whichwould simplify (13) and lead to a solution that can be implemented into a hydraulic fracturing simulator.The “full” solution can then be used to estimate the error in the approximation.Clearly, if there is no gravity (i.e. g= 0 and consequently Gρ = 0), equation (13) becomes an algebraicequation that is easy to solve. Unfortunately, the absence of gravity leads not only to a simplified slurryflow, but the model also loses the ability to capture gravitational settling of particles. However, the latterphenomenon is very important and needs to be accounted for at least approximately. To motivate theapproximation we are about to make, we observe from the plots of φ/φm in Fig. 2 that the proppantdistribution does not change appreciably as Gρ is varied (except for Gρ = 100) or as cosψ varies from 0 to1. Thus to simplify the solution and to keep the gravitational settling at the same time, it is assumed thatthe gravity does not affect the particle distribution, i.e. terms with Gρ can be neglected in (9a). In otherwords, the solution for φ is approximated byφ(s, 〈φ〉, Gρ, cosψ) ≈ φ(s, 〈φ〉, 0, 0), (22)where s = y/w. Note that this assumption is equivalent to replacing 〈φ〉 by φ in the gravity term in (7b),instead of neglecting the whole gravity term when the gravity is assumed to be negligible. As a result,the gravity term that causes the gravitational settling is still accounted for, but in an approximate fashion.In terms of the error, there are two main quantities that are important for hydraulic fracturing, namelythe slurry flux and the proppant flux. For the cases considered on the Fig. 2, i.e. Gρ = {1, 4, 10, 100},the corresponding error (relative to the Gρ = 0 solution) in terms of the absolute value of the slurry fluxand the particle flux is es = {0.02, 0.01, 0.03, 0.14} and ep = {0.04, 0.07, 0.06, 0.01} for cosψ = 0, es ={0.03, 0.11, 0.18, 0.24} and ep = {0.005, 0.04, 0.08, 0.09} for cosψ = 12 , and es = {0.07, 0.20, 0.29, 0.29} andep = {0.06, 0.15, 0.20, 0.15} for cosψ = 1. Here the slurry flux is calculated as an integral of vp0 , while theparticle flux is calculated as an integral of φvp0 . This shows that even in the worst case of Gρ=100, the fluxerror is less than 30%. Also, there is an error in the direction of the flux, but the direction can be correctedby the presence of extra hydrostatic pressure gradient inside the fracture. So the error in the direction of theflow affects the hydrostatic pressure inside the fracture. Note that for other average particle concentrations,the order of magnitude of the error does not change appreciably. For instance, for φ¯ = 〈φ〉/φm = 0.5,the corresponding errors are es = {0.03, 0.11, 0.14, 0.16} and ep = {0.04, 0.13, 0.19, 0.23} for cosψ = 0,es = {0.01, 0.02, 0.03, 0.09} and ep = {0.02, 0.05, 0.07, 0.17} for cosψ = 12 , and es = {0.02, 0.06, 0.08, 0.01}and ep = {0.02, 0.06, 0.06, 0.04} for cosψ = 1.To estimate the error distribution over the fracture, it is useful to provide an estimate for the parameterGρ. One may calculate the pressure gradient as |∇p˜f| = 12µf|V |/w2, where µf is the intrinsic fluid viscosity,w is the width of the channel, and V is the average fluid velocity. By taking µf = 0.1 Pa·s, w = 1 cm,9〈φ〉 = 0.1, ρp−ρf = 1300 kg/m3, |V | = 0.1 m/s, one may estimate Gρ asGρ =〈φ〉(ρp−ρf)gw212µf|V | ≈ 1, (23)which shows that the effect of gravity is relatively small for the considered parameters. However, Gρ mightincrease noticeably for wider and slower fractures. Typically, hydraulic fractures are narrow and propagaterapidly at early stages, while they become wider and slow down later in their evolution. In this case, it isimportant to monitor the criterion (23) for mature fractures. Note that for a better estimate, one needsto insert the effective viscosity in (23), which further decreases the value of Gρ. Also note that both thefracture width w and the average fluid velocity V reach their respective maxima at the wellbore (the velocityis singular if a point source is used). In this case, the parameter Gρ reaches its maximum some distanceaway from the wellbore, and is negligible near the wellbore and near the crack tip. Given the fact, that theparameter Gρ may vary significantly over the fracture and may cause relatively small error only in a localizedregion, the solution (22) can be used throughout the fracture as an adequately accurate approximation.Using the approximation (22) and (10), equation (13) can be inverted to findA(φ) = F(|∇p˜f|y − y0pp), (24)where F is a function that can be calculated numerically (or evaluated using the formula for the roots of a4th degree polynomial). For future reference it is noted that the asymptotic behaviour of the function F isF (s) s→0= F0(s) =2s5φm, F (s) s→∞= F∞(s) =√s. (25)Using the definition of A(φ) from (8b), one may express the particle concentration from (24) asφ(y) = φm1 + F(|∇p˜f|y − y0pp) , y0 < y 6 12w,φ(y) = φm, 0 6 y 6 y0. (26)It is also useful to introduce the averaged particle concentration normalized by φm asφ¯ = 〈φ〉φm= 2w∫ w/20φ(y)φmdy, (27)so that 0 6 φ¯ 6 1. By inverting the latter relation, the particle pressure and the half-width of the rigidcluster zone can be computed aspp = w|∇p˜f|Π(φ¯), y0 = µ1wΠ(φ¯), (28)where Gρ=0 (or φ = 〈φ〉) is used in (11) to calculate y0, and Π(φ¯) is a function that is evaluated numerically.Note that the asymptotic behaviour of the function Π(φ¯) for small and high average volume concentrationsof particles isΠ(φ¯) φ¯→0= Π0(φ¯) = 18 φ¯2, Π(φ¯) φ¯→1= Π1(φ¯) =12µ1−√5φm(1−φ¯)4µ31. (29)Equation (28) allows us to replace the particle pressure pp by the normalized average volume fraction φ¯, inwhich case (26) can be rewritten asφ = φm1 + F(max{ y/wΠ(φ¯)−µ1, 0}) . (30)100 0.2 0.4 0.6 0.8φ/φmy/w0 0.2 0.4 0.6 0.8|vp0|/max{|vp0|}¯=0.4¯=0.8¯=0.95¯=0.2¯=0.4¯=0.8¯=0.95¯=0.2Figure 3: Normalized velocity of the slurry (left panel) and normalized particle volume fraction (right panel)versus y/w for different values of φ¯ = {0.2; 0.4; 0.8; 0.95}. The dashed line on the left panel indicates theparabolic profile associated with φ¯ = 0.The function Π(φ¯) is directly related to the particle concentration at the wall through the relationφw =φm1 + F( 12Π(φ¯) − µ1) . (31)so that the above relation φw=φw(φ¯) can be alternatively used instead of Π(φ¯).To find the velocity distribution, we integrate (8a), use (28) and the fact that vp0 ∝ ∇p˜f to obtainvp = vp|y=w/2 + vp0 = vp|y=w/2 −w2Π(φ¯)µf ∇p˜f∫ 1/2max{y/w,µ1Π(φ¯)}[F( sΠ(φ¯)−µ1)]2ds. (32)As an illustration, Fig. 3 shows the velocity profiles given by (32) and the corresponding particle concentrationprofiles, computed using (30) versus y/w. Fig. 3 clearly indicates the presence of the rigid zone, however,the size of this zone becomes more pronounced only for high values of φ¯, i.e. high concentrations. Also notethat the velocity profile deviates from a parabolic shape for higher φ¯ and has a blunted profile, which wasalso shown in [6, 17].To capture the difference between the fluid and particle velocities, the slip velocity can be deducedfrom (7a), (8a), (10), (13) and (16) as∆v = vf − vp = − 1− B(φ)η(φ) ∇p˜f + (ρp−ρf)gezφη(φ) , (33)where the φ = φ(y) is taken from (30) and B is given in (20). Equation (33) signifies that there are twomechanisms for the slip velocity. The first is due to the fact that the particles collide and introduce shearresistance, which slows them down relative to the fluid. At the same time, the second term in (33) representsgravitational settling. Fig. 4 shows the normalized slip velocity for the case when g= 0, i.e. |∆v∇p˜f |/v∇p˜f∗ ,and for the case when ∇p˜f = 0, i.e. |∆vg|/vg∗ , versus y/w for different values of φ¯. The slip velocities arenormalized respectively byv∇p˜f∗ =2a29µf |∇p˜f|, vg∗ =2a29µf (ρp−ρf)g.As can be seen from the figure, both slip velocities have qualitatively similar profiles, although their magni-tudes are different.110 0.2 0.4 0.6 0.8∆vx/v∆∗y/w0 0.2 0.4 0.6 0.800.∆vgx/vg∗y/w|vrp˜f |/vrp˜f⇤ | |/vg⇤¯=0.2¯=0.4¯=0.8¯=0.95¯=0.4¯=0.8¯=0.95¯=0.2Figure 4: Normalized slip velocity due to the pressure gradient |∆v∇p˜f |/v∇p˜f∗ (left) and normalized settlingvelocity |∆vg|/vg∗ (right) versus scaled y/w coordinate for different values of φ¯ = {0.2; 0.4; 0.8; 0.95}.To find vp|y=w/2, one may simplify (34) using (33) tovp|y=w/2 = −C11− B(φw)η(φw)∇p˜f + C2(ρp−ρf)gezφwη(φw)+ (C1−C2)(ρp−ρf)g|∇p˜f|2∂p˜f∂zφwη(φw)∇p˜f, (34)where the values of C1 and C2 are given in (19).To summarize, the particle, fluid and slurry velocities are given respectively byvp = vp|y=w/2 + vp0 ,vf = vp|y=w/2 + vp0 + ∆v, (35)vs = φvp + (1−φ)vf = vp|y=w/2 + vp0 + (1−φ)∆v,where vp|y=w/2 is given in (34), vp0 stems from (32), while ∆v is calculated according to (33).3.4 Estimation of the characteristic length to establish a steady flowThe governing equations adopted for this study open the possibility to estimate the characteristic lengthrequired to establish a steady flow, which is an alternative to (2). By using the current problem geometry,assuming that initially φ(y)=const., and that the velocity profile is parabolic, i.e. vp = 6w−2(w2/4−y2)〈vp〉,the y component of the particle pressure gradient can be estimated from (8) as∂pp∂y = − µfA(φ)−2 12|〈vp〉|w2 ,where |〈vp〉| is an average velocity of the flow. Using the previous equation together with (3a), the ycomponent of the velocity of the particles can be found asvpy = −83φ(1−φ)α¯(φm−φ)2a2w2 |〈vp〉|.The characteristic time required to establish steady flow is t∗=w/(2|vpy |), while the corresponding length isL= |〈vp〉|t∗, so thatLw =316(φm−φ)2φ(1−φ)α¯w2a2 . (36)12Relation (36) is in agreement with (2), compatible with the adopted model, and provides a more accurateestimate for different values of φ. As an example, for φ = 0.1, a = 0.4 mm and w = 1 cm, the length requiredto establish a steady flow is approximately equal to 4 m. For this set of problem parameters, if the fractureis longer than 4 m (while the width is 1 cm), then the assumption of steady flow is appropriate.4 Mass transport equationsThe goal of this section is to connect the developments of Section 3 and the hydraulic fracturing problem.To achieve this, the fracture is represented by a channel with a variable width w(x, z, t), and it is assumedthat the flow at any point is always in equilibrium or in a steady-state condition. To formulate the governingequations for the mass transport in a one-dimensional setting, it is necessary to integrate (3b) and (4b) overy from −w/2 to w/2 and to derive equations for averaged quantities (note that the coordinate system is thesame as that assumed in Section 3). By using symmetry, relation (27), and noting that for any functionh(t, s)2 ∂∂t∫ w/20h(t, s) ds = 2∫ w/20∂h(t, s)∂t ds+ h(t, w/2)∂w∂t ,equations (3b) and (4b) can be integrated to obtain∂w∂t +∇·qs + g = 0,∂wφ¯∂t +∇·qp = 0, (37)where the following boundary conditions have been usedvpy |y=w/2 =12∂w∂t , vpy |y=0 = 0, vfy|y=w/2 =12∂w∂t +g1−φ|y=w/2, vfy|y=0 = 0.Here the leak-off, g, appearing in (37a), is introduced through the above velocity boundary condition andenables us to account for the fluid losses due to the filtration into surrounding rock.The slurry and particle fluxes are computed on the basis of (32)–(35) asqs = − w312µfQs(φ¯)∇p˜f − a2w12µfD(φ¯)∇p˜f − a2w12µf (ρp−ρf)gezGs(φ¯)−a2w12µf(ρp−ρf)g|∇p˜f|2∂p˜f∂z S(φ¯)∇p˜f,qp = −w(w2−w2cr)12µf Qp(φ¯)∇p˜f − a2w12µf (ρp−ρf)gezGp(φ¯)−a2w12µf(ρp−ρf)g|∇p˜f|2∂p˜f∂z φ¯S(φ¯)∇p˜f, (38)13whereQs(φ¯) = 24 Π(φ¯)∫ 1/20∫ 1/2max{s,µ1Π(φ¯)}[F( s′Π(φ¯)−µ1)]2ds′ ds,Qp(φ¯) = 24 Π(φ¯)∫ 1/20φ(s)φm∫ 1/2max{s,µ1Π(φ¯)}[F( s′Π(φ¯)−µ1)]2ds′ ds,D(φ¯) = 163∫ 1/201−B(φ(s))φ(s) (1−φ(s))α¯+1 ds− 831−B(φw)φw(1−φw)α¯B(φw),Gs(φ¯) = 83µf(1−φw)α¯µf+µp(φw)− 163∫ 1/20(1−φ(s))α¯+1 ds, (39)S(φ¯) = 83(B(φw)−µfµf + µp(φw))(1−φw)α¯,wcr(φ¯) = a√831−B(φw)φwQp(φ¯)φ¯(1−φw)α¯B(φw),Gp(φ¯) = 83µfφ¯(1−φw)α¯µf+µp(φw),φw(φ¯) is given in (31), while B is defined in (20). To understand the nature of the (1−φ|y=w/2)−1 multiplierin the leak-off velocity, one can add a thin virtual layer of pure fluid right next to the boundary. Then they component of the velocity associated with the leak-off in this virtual layer is g. However, since the fluidis incompressible, the jump in the volume fraction φ between the thin layer of fluid and bulk slurry causesa jump in the velocity. To preserve the volume, one should have (1−φ|y=w/2)vˆfy|y=w/2 = g, which impliesvˆfy|y=w/2 = (1−φ|y=w/2)−1g, where vˆfy|y=w/2 is the component of vfy|y=w/2 that accounts for the leak-off.Noting that φw ≈ 12 φ¯φm for small values of φ¯, the asymptotic behaviour of the functions defined in (39)can be computed with the help of (25) and (29) asQs0(φ¯) = 1, Qs1(φ¯) =45φmµ1(1−φ¯)3/2,Qp0(φ¯) =65 φ¯, Qp1(φ¯) =45φmµ1(1−φ¯)3/2,D0(φ¯) = −53((α¯+2)φm−1)φ¯, D1(φ¯) =8(1−φm)α¯3φm,Gs0(φ¯) =23 φ¯φm(2α¯−1), Gs1(φ¯) = −83(1−φm)α¯+1, (40)S0(φ¯) =53φmφ¯, S1(φ¯) =163( 25φm)2√5φmµ1(1−φm)α¯(1−φ¯)1/2,wcr,0(φ¯) =53a, wcr,1(φ¯) =4µ3/41 a31/251/4φ1/2m(1−φm)α¯/2(1−φ¯)−1/2,Gp0(φ¯) =83 φ¯, Gp1(φ¯) =3215(1−φm)α¯ 1−φ¯φm.For the purpose of fast numerical evaluations of the functions (39), their values are first calculated accuratelyfor a small set of φ¯, then divided by the corresponding asymptotic behaviour (40) and finally stored. Inthis case, the stored quantities are all O(1), which allows us to use spline interpolation and the asymptoticformulas (40) to restore the values of all functions in a fast and accurate manner. As an illustration Fig. 5 plotsthe functions Qs, Qp, D, wcr/a, Gs, Gp, S and φw/φm versus the normalized average particle concentrationφ¯.The critical width, appearing in (39), effectively precludes the proppant in areas where the fracture isnarrow. But, as can be seen in the Fig. 5, for most particle concentrations, this width is on the order of 1.5×a,140 0.5 100.51φ¯Qs0 0.5 100.10.2φ¯Qp0 0.5 1−0.6−0.4−0.200.2φ¯D0 0.5 100.511.52w cr/aφ¯0 0.5 1−φ¯Gs0 0.5φ¯Gp0 0.5 100.10.2φ¯S0 0.5 100.51φ w/φmφ¯Figure 5: Variation of functions Qs, Qp, D, wcr/a, Gs, Gp, S and φw/φm versus normalized average particleconcentration φ¯.15while it goes unbounded for a very narrow region near φ¯≈1. The particles, however, can not physically bein a channel whose width is smaller than the particle diameter. Moreover, the bridging of the proppant mayoccur earlier when the fracture width is equal to several particle diameters. For this reason, there is a needto introduce a “blocking” function, which would prevent the proppant in the narrow regions. To this end,let us introduceB(wa)= 12H( w2a−N)H(wB − w2a)[1 + cos(piwB − w2a)]+H(w − wB2a), (41)where N represents “several” particle diameters (i.e. three), H denotes Heaviside step function, while wB=2a(N+1), which provides a continuous vanishing of the function and helps the numerical implementation.To effectively prevent the proppant from moving into the narrow regions, the proppant flux in (38) needs tobe multiplied by (41), in which case the proppant flux that is calculated according toqpB = B(wa)qp (42)accounts for the proppant stalling in the narrow fracture regions.4.1 Comments on the modelIt is instructive to understand the physical meaning of all terms in the expressions for the fluxes in (38).The first term in the slurry flux is a Poiseuille flow type term, where Qs(φ¯) can be understood as theinverse of the effective viscosity of the mixture. As φ¯ approaches 1, the effective viscosity goes to infinityas (1−φ¯)−3/2, as can be seen from (40). The second term in the slurry flow can be related to Darcy’s law.Indeed, according to Darcy’s law, the total flux through the channel is proportional to the pressure gradient,the permeability, and the width of the channel. The permeability is proportional to the average squared poresize, which is, in turn, proportional to a2 and some dimensionless function of φ¯. Since a/w  1, for mostparticle concentrations, the flux is dominated by Poiseuille’s law. However, in the plugged regions (φ¯ ≈ 1)Qs(φ¯) ≈ 0, while D(φ¯) > 0. So that the flux is governed mainly by Darcy’s law. The remaining two termsrepresent the effect of the gravity. In particular, they signify that when the slurry is in a static position, i.e.qs =0, the pressure is not equal to hydrostatic pressure. Instead, there is a pressure gradient in addition tothe hydrostatic pressure, which balances the shear stresses at the wall caused by the gravitational settlingof particles. With regard to the flux of particles, the first term corresponds to the flux of particles “carried”by the fluid, i.e. the advective term. The last two terms of the particle flux capture the effects associatedwith gravitational settling of the particles. One may also observe strange non-monotonic variation of wcr inFig. (5). This is, unfortunately, one drawback of the constitutive model. This constitutive model is basedon the fitting. So while µp(φ¯) may be accurate, the accuracy of its derivative may deteriorate significantly.Function B (20), for instance, depends on the derivative of µp(φ¯), and wcr depends on B, see (39).4.2 Simplified modelDespite the fact that (38) gives a complete answer for a broad range of parameters, the fact that the smallestfracture width where the proppant can be present is bounded by several particle diameters can be furtherused to simplify the result. The function wcr, shown in Fig. 5, is very steep, but monotone in the regionφ¯ ≈ 1, so can be inverted and be written in the form φ¯wcr (wcr/a). Let’s consider the case when N=3 in (41)and try to estimate the quantity w2−w2cr that appears in (38). In the worst case of w= 6a, this quantityvanishes for 1−φ¯wcr (6) = O(10−4). Moreover, 1−φ¯wcr (2) = O(10−3), in which case w2cr/w2 < 1/9 for allpossible values of w (for which proppant can be inside the fracture). This shows that if the term with wcr isneglected, a noticeable error may be introduced only for near-maximum concentrations 1−φ¯ < O(10−3). Toshow whether the variations of φ¯ in the region 1−φ¯ < O(10−3) affect the slurry flow, one needs to comparethe terms with Qs(φ¯) and D(φ¯) in the slurry flux. To find the particle concentration, which characterizesthe transition from Poiseuille’s to Darcy’s flow regime, one needs to express the relation w2Qs(φ¯)=a2D(φ¯)in the form φ¯P−D(w/a). For the worst case w = 6a, it can be shown that 1− φ¯P−D(6) = O(10−2). As aresult, wcr can be neglected as D(φ) already dominates the flow for 1−φ¯ < O(10−3) and its variations in this16range are negligible (so that the error of 0.1% in φ¯ does not lead to noticeable changes in the slurry flow).Moreover, the variations of D are approximately equal to 5% when 1−φ¯ = O(10−2), in which case D(φ¯) canbe approximated by φ¯D(1) (note that D(φ¯) should vanish for φ¯=0). The estimation of smallness of w2cr/w2is valid for the extreme case when w = 6a. In situations when w/a > 6, the asymptotic behaviour of Qsand wcr (40) can be used to find that 1−φ¯P−D(w/a) ∝ (a/w)4/3 for the Poiseuille-to-Darcy transition, while1−φ¯wcr (w/a) ∝ (a/w)2. This shows that the separation between two characteristic values of 1−φ¯ increasesfor bigger fracture widths (while they both decrease), which allows us to safely neglect wcr in (38) and usean approximation for D(φ¯).Since a/w  1, the gravity terms featured in (38) become important only for small values of the verticalpressure gradient, i.e. ∂p˜f/∂z = O(a2/w2(ρp−ρf )g). This situation may occur when the hydraulic fractureis surrounded by stress barriers, since some hydraulic fracturing models even assume constant pressure ineach vertical cross-section, see [18]. Due to the presence of stress barriers, the pressure gradient in thesefractures is predominantly horizontal, which could make the vertical component of the pressure gradientsufficiently small to cause ∂p˜f/∂z = O(a2/w2(ρp−ρf )g). Note, however, that the absolute value of thepressure gradient may be big and so Gρ, defined in (12), could be small. The fact that ∂p˜f/∂z/|∇p˜f|  1allows us to neglect terms in (38) that are proportional to ∂p˜f/∂z/|∇p˜f|  1, which reduces the complexityof the fluxes. Moreover, it is useful to introducepˆf = pf − pfh = pf + ρfgz + (ρp−ρf )gφmφ¯z + (ρp−ρf )g∫ z0a2Gs(φ¯)w2Qs(φ¯) + a2φ¯D(1) dz′, (43)which allows us to simplify (38) (also using the simplifications discussed in the previous paragraph) toqs = − w312µf Qˆs(φ¯, wa)∇pˆf,qpB = B(wa)Qˆp(φ¯, wa)qs −B(wa) a2w12µf (ρp−ρf)gezGˆp(φ¯, wa), (44)whereQˆs(φ¯, wa)= Qs(φ¯) + a2w2 φ¯D(1),Qˆp(φ¯, wa)= w2Qp(φ¯)w2Qs(φ¯) + a2φ¯D(1) , (45)Gˆp(φ¯, wa)= Gp(φ¯)− w2Gs(φ¯)Qp(φ¯)w2Qs(φ¯) + a2φ¯D(1) .Here pfh is the “true” hydrostatic pressure in the sense that there is no motion of the slurry if pf = pfh. Interms of fluid driven fracture problems, hydrostatic pressure is important for buoyancy-driven fractures, seee.g. [19], while is commonly neglected for other industrial hydraulic fracturing problems. In the simplifiedmodel (44), functions (45) determine the motion of the slurry and proppant. The function Qˆs describesthe slurry flow, and, in particular, accounts for the Poiseuille-to-Darcy transition of the flow. Clearly, Qˆpdescribes convective proppant transport, while Gˆp captures gravitational settling. The fact that the modelcaptures Poiseuille-to-Darcy transition of the flow implies that both proppant transport and plugging (i.e.the formation of a zone with nearly maximum proppant concentration) are accounted for autonomously. Asan illustration, Fig. 6 plots the variations of Qˆp and Gˆp versus φ¯ for different values of w/a = {6, 10, 100}.One can see that the parameter w/a is important only in the regions with high particle concentration andleads to a smooth transition of functions Qˆp and Gˆp to zero. It should be noted that there is almost novisual difference between Qs and Qˆs. At the same time, qualitatively, there is a significant difference, sinceQs goes to zeros as φ¯→ 1, while Qˆs approaches a small constant. For the purpose of numerical calculations,the variation of functions Qˆp and Gˆp versus w/a can be approximated by either taking a limit w/a→∞ orw/a=6, depending on the numerical scheme. As Fig. 6 shows, the error induced by such an approximationmay slightly alter the dynamics of the plug (i.e. the region with a nearly maximum concentration), whilethe proppant transport at smaller concentrations will not be affected.170 0.5φ¯Gˆp0 0.5 100.51φ¯Qˆp 0.9 0.95φ¯Gˆpwa= 6wa= 10wa= 100wa= 6wa= 10wa= 100Figure 6: Variation of functions Qˆp and Gˆp versus normalized average particle concentration φ¯ for differentvalues of w/a = {6, 10, 100}.In the current simplified formulation, the functions Qˆs, Qˆp and Gˆp that enter (44) are based on thisspecific model for the slurry flow. In practical applications, however, all assumptions of the model maynot always be met and the model may deviate from the observations. For instance, the particles may beaspherical, or the mixture of particles with different sizes is used. If a similar model cannot be derived due toe.g. lack of a constitutive relation for the mixture, one possible solution is to keep all the terms in (44), butto modify the functions Qˆs, Qˆp and Gˆp to fit the data. For instance, one could measure the effective viscosityand permeability of the mixture and modify the functions accordingly. This formulation represents a broadclass of models, which could satisfy experimental observations, capture the Poiseuille-to-Darcy transition ofthe slurry flow, as well as particle transport and settling.5 SummaryThis paper studies the steady flow of a Newtonian fluid mixed with spherical particles based on a continuumapproach. The constitutive behaviour of the mixture is taken from an empirical model, where both theshear stress and particle pressure are expressed in terms of functions of the particle concentration. Since thepressure gradient and the gravity are not always collinear, and the problem under consideration is nonlinear,two-dimensional flow of the slurry needs to be considered. The solution is obtained in terms of the particleconcentration profile, particle velocity, and fluid velocity. According to the solution, the particles form arigid cluster in the centre, and its size depends on the average concentration of the particles over the width.The velocity profile is shown to deviate from the classical parabolic shape for higher concentrations due toa nonuniform distribution of particles. Also, the solution features a boundary layer, whose size is on theorder of the particle radius. While the effect of the boundary layer is minimal for relatively small particleconcentrations, it becomes crucial when the particle concentration is close to the maximum value. Due tothe nonlinearity of the problem, the solution depends on the ratio between the pressure gradient and gravityforce, as well as on the direction of the pressure gradient relative to the gravity. This, however, makes theresult less valuable for hydraulic fracturing applications, as the implementation of such a complex model intoa hydraulic fracturing simulator becomes extremely cumbersome. For this reason an approximate solutionis introduced, in which the effect of gravity is accounted for approximately. This approximate solutioncoincides with the “full” solution for the case of buoyant particles, while it deviates when the gravity effectsare considerable. To provide an estimate for the error that is introduced by the approximation, the “full”solution is compared to the approximate solution for different regimes of the flow. It is shown that even inthe worst case when the pressure gradient is much weaker than the gravity force, the errors in the slurryflux and proppant flux do not exceed 30%. For more realistic parameters, the error becomes even lower, onthe order of 10%.The simplified approximate solution of the channel flow allows us to calculate the average fluxes for18the slurry and the particles, which is the basis for the analysis of hydraulic fracturing by a slurry and thecorresponding proppant transport problem. The slurry flux has two terms, one Poiseuille-law-type term withan effective viscosity (which goes to infinity as the concentration reaches a critical value), and a Darcy-law-type term, where the average velocity is proportional to the particle size squared and the pressure gradient.The flux of the particles also has two terms, one proportional to the pressure gradient (which can also beexpressed in terms of the slurry flux), a nonlinear function of concentration and width, and another relatedto the gravitational forces. The first term describes the transport of the particles by the fluid, i.e. advectiveterm, while the second term describes gravitational settling. The simplified model (44) is now in a form thatcan be implemented in a hydraulic fracturing model that couples elasticity with slurry flow, and capturesproppant transport, gravitational settling, and plug formation.AcknowledgementsThe authors gratefully acknowledge the support of the British Columbia Oil and Gas Commission and theNSERC discovery grants program. In addition, the authors would like to thank Professor E. Detournay(University of Minnesota) for stimulating the discussions at the inception of this project.References[1] M.J. Economides and K.G. Nolte, editors. Reservoir Stimulation. John Wiley & Sons, Chichester, UK,3rd edition, 2000.[2] J. Adachi, E. Siebrits, A. Peirce, and J. Desroches. Computer simulation of hydraulic fractures. Int. J.Rock Mech. Min. Sci., 44:739–757, 2007.[3] F. Boyer, E. Guazzelli, and O. Pouliquen. Unifying suspension and granular rheology. Phys. Rev. Lett.,107:188301, 2011.[4] J.J. Stickel and R.L. Powell. Fluid mechanics and rheology of dense suspensions. Annu. Rev. FluidMech., 37:129–149, 2005.[5] D. Leighton and A. Acrivos. The shear-induced migration of particles in concentrated suspensions. J.Fluid Mech., 181:415–439, 1987.[6] P.R. Nott and J.F. Brady. Pressure-driven flow of suspensions: simulation and theory. J. Fluid Mech.,275:157–199, 1994.[7] J.F. Brady and J.F. Morris. Microstructure of strongly sheared suspensions and its impact on rheologyand diffusion. J. Fluid Mech., 348:103–139, 1997.[8] J.F. Morris and F. Boulay. Microstructure of strongly sheared suspensions and its impact on rheologyand diffusion. J. Rheol., 43:1213–1237, 1999.[9] D. Lhuillier. Migration of rigid particles in non-brownian viscous suspensions. Phys. Fluids, 21:023302,2009.[10] G. Segre and A. Silberberg. Radial particle displacements in poiseuille ow of suspensions. Nature,189:209–210, 1961.[11] A. Einstein. A new determination of molecular dimensions. Annal. Physik, 4:289–306, 1906.[12] G. Batchelor and J. Green. The determination of the bulk stress in a suspension of spherical particlesto order c2. J. Fluid Mech., 56:401–427, 1972.[13] M. Ouriemi, P. Aussillous, and E. Guazzelli. Sediment dynamics. part 1. bed-load transport by laminarshearing flows. J. Fluid Mech., 636:295–319, 2009.19[14] J. Garside and M. R. Al-Dibouni. Velocity-voidage relationships for fluidization and sedimentation insolid-liquid systems. Ind. Eng. Chem. Process Des. Dev., 16:206–214, 1977.[15] R.H. Davis and A. Acrivos. Sedimentation of noncolloidal particles at low reynolds numbers. Ann. Rev.Fluid Mech., 17:91–118, 1985.[16] F. Boyer. Suspensions concentre´es: expe´riences originales de rhe´ologie. PhD thesis, Aix-MarseilleUniversite´, 2011.[17] D. Eskin and M.J. Miller. A model of non-newtonian slurry flow in a fracture. Powder Technology,182:313–322, 2008.[18] J. I. Adachi, E. Detournay, and A. P. Peirce. An analysis of classical pseudo-3D model for hydraulicfracture with equilibrium height growth across stress barriers. Int. J. of Rock Mech. and Min. Sci.,47:625–639, 2010.[19] J.R. Lister. Buoyancy-driven fluid fracture: the effects of material toughness and of low-viscosityprecursors. J. Fluid. Mech., 210:263–280, 1990.20


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items