ISSN 0015-4628, Fluid Dynamics, 2014, Vol. 49, No. 2, pp. 208–221. © Pleiades Publishing, Ltd., 2014. Original Russian Text © S.A. Boronin, A.A. Osiptsov, 2014, published in Izvestiya Rossiiskoi Akademii Nauk, Mekhanika Zhidkosti i Gaza, 2014, Vol. 49, No. 2, pp. 80–94.
Effects of Particle Migration on Suspension Flow in a Hydraulic Fracture S. A. Boronina and A. A. Osiptsovb a Moscow
Research Center of the Schlumberger company, Leningradskoe sh. 16A str. 3, Moscow, 125171 Russia b Moscow State University, Faculty of Mechanics and Mathematics, Leninskiye Gory 1, Main Building, GSP-1, Moscow, 119991 Russia e-mail:
[email protected] Received March 19, 2013
Abstract—An asymptotic model of a hydraulic-fracture flow of a sedimenting concentrated suspension is formulated on the basis of the two-fluid approach with account of transverse particle migration. In the thin-layer approximation, a two-dimensional system of equations averaged across the fracture is constructed with account for a nonuniform distribution of the particle concentration. As compared to the similar model without particle migration, the averaged two-dimensional equations contain modified coefficients which explicitly depend on the width of the flow core occupied by the particles. Using the model constructed, a numerical simulation is performed, which shows that the particle migration towards the fracture center results in the increase in the depth of particle penetration into the fracture and the suppression of gravitational convection in the vicinity of the leading front. The calculations are compared with available experimental data and an analytical formula for the height of the dense packed sediment. A good agreement between the analytical theory, the experiments, and the two-dimensional calculations is attained. Keywords: viscous fluid, hydraulic fracturing, two-fluid approach, migration, sedimentation, suspension, fracture, particle. DOI: 10.1134/S0015462814020094
The aims of this work are to construct an asymptotic two-fluid model of suspension flow in a vertical hydraulic fracture with account of a finite particle volume fraction and velocity slip between the phases, and to study the effect of transverse particle migration on the suspension flow. The study of this phenomenon is of interest due to the need in constructing mathematical models which can be implemented in commercial computer codes used for numerical simulation of suspension flows in a fracture, as applied to the hydraulic fracturing technology in oil and gas reservoirs [1]. The hydraulic fracturing technology is based on pumping a liquid-solid mixture in the well at high pressure to create vertical hydraulic fractures in the reservoir. During pumping the mixture into the fracture, a suspension flow develops, in which the particles are carried with the liquid in the horizontal direction, settle down under the gravity, and may migrate in the transverse direction. When the pumping stops, a closed fracture filled with the granular material of densely packed particles creates a highly permeable porous conduit to filter hydrocarbons from the reservoir to the well. Recently, the simulation of suspension flow in fractures has become particularly topical due to a wide spread of the multistage fracturing technology for the shale gas recovery. In the framework of a two-fluid approach, gravitational convection of suspensions with account of the effects of rapid sedimentation in inclined vessels (the Boycott effect [2]) was first considered in [3, 4]. A study of gravitational convection of suspensions in inclined fractures was performed in [5]. The role of the added mass, Archimedes, and Basset–Boussinesq forces in the problems of gravitational convection of suspensions was investigated in [6]. In [7], in the framework of a two-fluid model of a dilute suspension a numerical simulation of the Boycott effect in a tilted square container was performed in the two-dimensional formulation. 208
EFFECTS OF PARTICLE MIGRATION ON SUSPENSION FLOW
209
A general three-dimensional two-fluid model of gravitational convection with account for the non-stationary components of the interphase force and a finite particle volume fraction was constructed in [8]. We will consider the models of suspension flow in hydraulic fractures, known in the literature. In [9], a general formulation of the model of particle transport in suspension flow through a vertical fracture is given in the thin-layer approximation. The equations are obtained using the ‘effective fluid’ model in terms of the mass fluxes of the particles, fluid, and mixture, the expressions for which should be found from the solution of local problems. In [10], explicit formulas for the mass fluxes, which serve as the closure relations for the general equations of the model [9], are proposed. The particle transport model used in the commercial computer code for modeling the growth of a hydraulic fracture is described in [11]. In these models, the fracture is treated as a vertical flat channel with smooth walls. The channel width is given, it can be variable in space but does not depend on time. In more complex models, used in modern commercial computer codes for modeling the hydraulic fracturing [11], the coupled problem of the fracture growth and the suspension flow in the fracture is considered. In solving the hydrodynamic problem, the variation of the fracture geometry is regarded as quasi-stationary. When modeling a suspension flow in a fracture, three substantially different length scales can be distinguished: the particle size, the fracture width, and the fracture length (height). In the papers mentioned above, 2D models are obtained on the largest length scale by averaging across the fracture thickness under assumption that the particle concentration profile is known. Typically, the particle migration across the fracture is neglected, and the particle concentration profile is assumed to be uniform. The lateral migration of particles in suspension channel flow can be caused by different mechanisms, depending on the order of magnitude of the particle volume fraction. For a dilute suspension, the transverse motion of the particles is induced by an inertial lateral force, which arises due to the fluid inertia, the phase velocity slip (particle settling), and the locally shear nature of the flow around the particle [12]. In [12], an expression is obtained for the inertial lateral force exerted on a particle settling in the horizontal viscous flow in a channel with plane vertical walls. The problem of inertial migration of particles in a dilute suspension flow through the initial sections of a plane channel (circular tube) with neglect of particle settling is considered in [13]. It is shown that the particles migrate towards two planes (rings) located at a certain distance from the channel (tube) walls, with regions devoid of particles being formed near the walls. The particle migration in suspension plane-channel flow with account of the force [12] is investigated in [14]. It is found that under the action of this force the particles migrate from the walls towards the central plane of the channel, leaving pure-fluid regions near the walls. In the case of a concentrated suspension, the migration is diffusive in nature and is attributable to the flow shear and transverse gradient of particle concentration [15], which result in the formation of a nonuniform concentration profile with a maximum on the channel axis. For oil recovery applications, both migration mechanism are of interest, since the range of variation of the particle volume fraction comprises both dilute (fracturing of low-permeability rocks in the recovery of a shale gas) and concentrated slurries (standard hydraulic fracturing, for example, in sandstone). A characteristic feature of particle migration in suspension channel flows, regardless of the order of magnitude of the particle volume fraction, is the formation of near-wall low-concentration layers (or pure-fluid layers when the initial particle concentration is small) and a core flow with a high particle concentration. In general, the crosswise profile of the particle volume fraction can be approximated by a piecewise constant dependence (Fig. 1b), and this fact is essentially used in the present work. In previous publications, no systematic study of the particle migration effect on the suspension flow through the fracture has been performed. Only some attempts to take into account the macroscopic effects of migration were undertaken. In particular, in [10] the closure relations for the model [9] were obtained in the limiting case, when all the particles are accumulated in the neighborhood of the central plane of the fracture, with pure-fluid layers being formed near the walls. In general, the existing models of particle transport in the fracture [9, 11] are postulated within the ‘effective fluid’ approximation, in which the suspension is treated as a viscous fluid with density and viscosity depending on the particle volume fraction. In these models, due to the incompressibility of phase materials FLUID DYNAMICS
Vol. 49
No. 2
2014
210
BORONIN, OSIPTSOV
Fig. 1. Scheme of suspension flow in a fracture with lateral particle migration (a). Uniform (C1 ) and nonuniform (C2 ) transverse distributions of the particle concentration (b).
(fluid and particle material) it is assumed that the suspension can also be described by a model of effective incompressible fluid. A two-fluid model of suspension flow in a hydraulic fracture was first proposed in a short communication [16]. In earlier studies, it was assumed that the particles settle down relative to the mixture (rather than relative to the carrier phase) with a given velocity which is a function of the particle volume fraction. For a suspension flow through a vertical plane slot (fracture), the suspension mean-volume velocity profile was taken equal to Poiseuille one. In [16], it was demonstrated that these assumptions are valid only when the phase velocity slip is zero or the particle volume fraction is negligible. An analysis of existing models in comparison with the proposed two-fluid model was performed and a range of similarity parameters was determined, for which the existing effective-fluid model is applicable. The conditions were found, for which the two-velocity effects are important. It is shown that the differences between the model constructed previously on the basis of a heuristic one-velocity approach and the two-fluid model [16] are significant for flows of a dilute suspension with a low-viscosity carrier fluid (hydraulic fracturing in a shale), while in the case of finite particle volume fractions and high-viscosity fluids (standard fracturing, for instance, in a sandstone) the solutions for the two models practically coincide. At the same time, the coupled problem of fracture growth with suspension flow inside, described using the two-fluid model, has not been considered so far. Such a study will make it possible to estimate the quantitative impact of the two-fluid effects on the law of fracture propagation. Bearing in mind the need in developing models of particle transport in a fracture to meet the challenges of new technologies (such as hydraulic fracturing in shales), in this paper we propose a generalization of the averaged asymptotic two-dimensional equations of suspension flow through a fracture [16] to the case of a flow with a substantial lateral migration of the particles. The two-dimensional system of equations will be obtained by averaging across the fracture with account for a given nonuniform particle concentration profile, formed as a result of transverse migration of the particles. Since the problem of finding the actual FLUID DYNAMICS
Vol. 49
No. 2
2014
EFFECTS OF PARTICLE MIGRATION ON SUSPENSION FLOW
211
crosswise particle concentration profile with account of particle migration is beyond the scope of this work, below we will study the effect of particle migration on the macroscale transport, gravitational convection, and suspension sedimentation in the fracture using a prescribed nonuniform crosswise profile of the particle concentration. 1. FORMULATION OF THE PROBLEM We consider a three-dimensional unsteady isothermal laminar flow of a suspension through a vertical hydraulic fracture in the gravity field g. The carrier phase is a viscous incompressible fluid with density ρ 0f . The dispersed phase consists of identical non-colloidal spherical particles of radius σ , mass m, and material density ρ p0 . The fracture is assumed to be a vertical slot with a variable width w and smooth walls. A coordinate system Oxyz is introduced, with the axes x and y, directed along the horizontal and vertical, and z-axis directed perpendicular to the central plane of the fracture (Fig. 1a). The unit vectors of the coordinate system are signified as e1 , e2 , and e3 . The suspension is regarded as two interpenetrating and interacting continua: the particulate medium and the carrier phase. The particulate medium is characterized by the volume fraction C, the number concentration n p , the density ρ p = Cρ p0 , and the mean-mass velocity v p . The carrier phase is described by the averaged density ρ f = (1 − C)ρ 0f and the mean-mass velocity v f . The velocity components of the phases in the projection on the coordinate axes will be denoted as v = (u, v, w). Within the two-fluid approach [17], the differential equations of the mass and momentum balance laws for each phase in dimensional form read:
∂ρf + ∇(ρ f v f ) = 0, ∂t
∂ ρp + ∇(ρ p v p ) = 0, ∂t
(1.1)
ρf
df vf = −∇p f + ∇ j τ if j ei + ρ f g − n p F p , dt
(1.2)
ρp
dpvp = −∇p p + ∇ j τ pi j ei + ρ p g + n p F p , dt
(1.3)
∂vf df vf = + (v f ∇)v f , dt ∂t
dpvp ∂ vp = + (v p ∇)v p . dt ∂t
The boundary conditions in dimensional form are: x = 0 : v f = v f 0 (t, y, z), z=±
w : u f = v f = 0, 2
C = C0 (t, y, z), wf = ±
1 ∂w ± wl , 2 ∂t
(1.4) w p = 0.
At the inlet of the fracture, the fluid velocity v f 0 and the particle volume concentration C0 are given. Since the fracture is surrounded by a porous medium, on the walls of the fracture the normal fluid velocity wl is specified, which takes into account the outflow of the fluid through the permeable walls of the fracture. Also, the condition of the absence of outflow of particles is specified, since, by assumption, the particle size is greater than the pore size. In applications, the carrier phase is an aqueous solution of a polymer, which has the property of shear thinning. Accordingly, as a rheological model for the suspension as a whole, as a rule, a power model with a yield stress (the Herschel–Bulkley model [18]) is used, or for simplification the model of a linear viscous fluid is adopted. In this paper, we assume that the carrier phase is a Newtonian fluid with viscosity μ0 , however the resulting model obtained below can be easily generalized to the case of a more complex rheology by analogy to known models, following, for example, [11]. FLUID DYNAMICS
Vol. 49
No. 2
2014
212
BORONIN, OSIPTSOV
We introduce the mean-volume and mean-mass velocities of suspension vv and vm :
vv = (1 − C)v f + Cv p ,
vm =
(1 − C)ρ 0f v f + Cρ p0 v p (1 − C)ρ 0f + Cρ p0
.
(1.5)
In the models known in the literature, the suspension is described in terms of the mean-volume velocity of the mixture, which is described by the Poiseuille law. This is valid only when the particle and fluid velocities are identical. Then, the mean-mass and mean-volume velocities of the suspension coincide and the phase velocities are equal. In this case, the flow of the mixture can be regarded as a motion of an effective medium and can be described by the Navier–Stokes equations for a viscous incompressible fluid, whose density and viscosity depend on the local particle volume fraction. As a result, for the velocity of suspension flow in a plane channel one can obtain a solution in the form of the Poiseuille law. In a general case, due to a finite volume fraction of the particles and a velocity slip of the phases (particle sedimentation) the mean-mass and mean-volume velocities of the slurry do not coincide and differ from the velocity of each phase. Then, the Poiseuille law for the mean-volume velocity of the suspension does not follow from the momentum equation and is thus an unjustified assumption. It seems possible to assert that the for a finite phase velocity slip the mixture flow cannot be described as a motion of an effective medium, and it is necessary to consider the motion of each phase separately within the two-continuum approximation, writing the balance equations of mass and momentum for each phase in terms of the mean-mass velocities of the phases with account for the interphase momentum exchange (see (1.1)). Moreover, when the phase velocity slip is nonzero, the densities of the particulate medium, the carrier phase, and the suspension are variable, and therefore the mean-mass velocities of the phases are divergent. Thus, the media of the particles and the carrier phase are compressible. Accordingly, we will assume that the stress tensor in the carrier phase can be taken as in a viscous compressible fluid with a viscosity depending on the particle volume fraction [19]:
τ if j
( ) 1 ij ij = 2μ0 μ (C) e − δ ∇v f , 3
( μ (C) = 1 −
C Cmax
)− β
,
μ (0) = 1,
Cmax = 0.65,
(1.6)
β = 1.89.
The dependence of the viscosity on the particle volume fraction is given by the Scott formula [21]. It is assumed that in the gravity-induced sedimentation the chaotic particle velocity is small, and it is possible to neglect the stress tensor in the particulate medium: p p = 0, τ pi j = 0 [19]. In the expression for the interphase force, we take into account the Stokes force with a correction for a finite particle volume fraction f (C) and the Archimedes force F p = 6πσ μ0 (v f − v p )
4 1 − πσ 3 ρ 0f g. f (C) 3
(1.7)
As shown in [8], the added mass, Basset–Boussinesq, and the non-stationary part of the Archimedes force can be neglected in the approximation of inertialess sedimentation, when the phase velocity relaxation length is much smaller than the length scale of the problem. This condition is justified below in the derivation of equations in dimensionless variables. FLUID DYNAMICS
Vol. 49
No. 2
2014
EFFECTS OF PARTICLE MIGRATION ON SUSPENSION FLOW
213
2. DERIVATION OF THE ASYMPTOTIC EQUATIONS We introduce the dimensionless variables (dimensional variables are denoted by a prime, where it is necessary to distinguish them from the dimensionless variables) x′ = Lx,
ρ p′
=
ρ 0f ρ p ,
v′f = U v f ,
τ
i j′
v′p = U v p ,
μ0U i j = τ , L
L t, U
t′ =
μ0U p, p = L ′
ρ ′f = ρ 0f ρ f , (2.1)
′
w = Lw,
w′l
= U wl .
Here, L is the characteristic length scale and U is the characteristic velocity scale. In dimensionless variables, from (1.1)–(1.3) we obtain the following equations: ] [ ∇ (1 − C)v f + Cv p = 0,
∂C + ∇(Cv p ) = 0, ∂t ] [ ] [ df vf dpvp −1 + ηC = −∇p f + ∇ j τ if j ei − Bu0 1 + C(η − 1) e2 , ε Re (1 − C) dt dt ( ) dpvp St η − 1 1 = (v f − v p ) − 2 ε St e3 , dt f (C) Fr η d ε= , L
St =
mU , 6πσ μ0 d
U Fr = √ , gd
Re =
ρ 0f U d μ0
,
η=
ρ p0 , ρ 0f
Bu0 =
Re
ε 2 Fr2
(2.2)
.
Here, d is the characteristic scale of the fracture width. In the first of Eqs. (2.2), the expression in brackets is the mean-volume velocity of the suspension vv (1.5). Based on the fact that vv is divergence-free, in the previous publications a conclusion (generally incorrect) was made that the suspension flow can be described using the Navier–Stokes equations for an incompressible fluid, written in terms of the mean-volume velocity vv . As indicated above, this is true only in the case when the phase velocity slip is absent or the particle volume fraction is negligible. It is assumed that the fracture width w varies only slightly in space. The thin layer approximation is used:
ε ≪ 1,
∇w ≪ 1,
St ∼ 1,
Fr ∼ 1,
Re ∼ 1,
η ∼ 1.
(2.3)
In particular, it follows that the characteristic phase velocity relaxation scale is much smaller than the characteristic scale of the problem (ε St ≪ 1), which, following [8], makes it possible to neglect the Basset–Boussinesq force, the added mass force, and the non-stationary part of the Archimedes force in the expression for the interphase momentum exchange (1.7). We introduce the new stretched variables, which in the thin layer approximation are of the order of unity: z = ε z0 ,
w = ε w0 ,
wl = ε w0l ,
w f = ε w0f ,
p = ε −2 p0 .
(2.4)
For simplicity, the superscript “0” is omitted below. The asymptotic order of the pressure is found from the condition that the resulting asymptotic equations will be least degenerate and self-consistent [20]. In the case under study, this means that the pressure gradient has the same order of magnitude as the viscous terms. Substituting (2.4) into Eq. (2.2), in the asymptotic limit (2.3) we obtain that the continuity equations FLUID DYNAMICS
Vol. 49
No. 2
2014
214
BORONIN, OSIPTSOV
do not change, while the momentum equations take the form [ ( )] ∂ pf ∂uf ∂ pf ∂ + , = 0, μ 0=− ∂x ∂z ∂z ∂z [ ( )] ∂ pf vf ∂ + − Bu(1 + C(η − 1)), μ 0=− ∂y ∂z ∂z ( ) St η −1 vs = − 2 e2 , v p = v f + vs , η Fr f (C) Bu =
(2.5)
ρ 0f g d 2 Re . = μ0U Fr2
We note that the particle lateral migration process on its own is not considered in this study. In the derivation of the averaged equations, we use a nonuniform cross-sectional profile of the particle volume fraction, which has been formed as a result of migration. The solution for the cross-sectional concentration profile obtained in [14] for the flow of a dilute settling suspension in a vertical slot can be approximated by a piecewise constant dependence (Fig. 1b). Thus, as a result of particle migration, at the channel center a high-concentration core is formed, leaving layers devoid of particles near the walls. The solution of Eqs. (2.5), constructed with account of the cross-sectional concentration profile shown in Fig. 1b, takes the following form. In the range z ∈ [−z0 , z0 ], we have the formulas: [ ( ) ) ( ] 1 w2 z 2 2 2 ∂p z − , + 1 − z0 uf = − 8 μ (C) 0 w/2 ∂x w2 vf = − 8
[(
(2.6)
( ) ) ( ) 1 z 2 2 2 ∂p z − + 1 − z0 μ (C) 0 w/2 ∂y
( ) ) ( )] ( z 2 1 + C(η − 1) 2 2 z0 − + 1 − z0 . + Bu μ (C) w/2 In the range z ∈ [−1, −z0 ] ∪ [z0 , 1], we obtain ( ( ) ) z 2 w2 ∂ p 1− , uf = − 8 ∂x w/2 ( )( ( ) ) z 2 w2 ∂ p + Bu 1 − . vf = − 8 ∂y w/2 Note that in the derivation of these formulas, we have used the condition μ (0) = 1. We introduce new functions, averaged across the fracture using the general formula 1 F(t, x, y) = w
w/2 ∫
f (t, x, y, z) dz. −w/2
From (2.2), using (2.6), the following two-dimensional averaged equations are obtained for the case of a FLUID DYNAMICS
Vol. 49
No. 2
2014
EFFECTS OF PARTICLE MIGRATION ON SUSPENSION FLOW
215
uniform particle concentration profile:
∂ wC + ∇(wCV p ) = 0, ∂t [ ] [ ] ) ( ) w3 ( ∂w =∇ ∇P + Bu 1 + C(η − 1) e2 − wCVs − 2vl 1 − C p , ∂t 12μ (C) Vf = − Vs = −
[ ] ) w2 ( ∇P + Bu 1 + C(η − 1) e2 , 12μ (C)
( ) St η − 1 f (C)e2 , Fr2 η
( f (C) = 1 −
V p = V f + Vs ,
C Cmax
)5 ,
Bu =
(2.7)
(2.8)
(2.9)
Re . Fr2
The expression for the particle settling velocity is taken in the semi-empirical form [11], which takes into account the hindered settling effect in a concentrated suspension. The settling velocity tends to zero as the particle volume fraction tends to that for close packing, which makes it possible to simulate the formation of a deposit of densely packed particles using a shock-capturing method. For a nonuniform particle concentration profile, we obtain:
∂ wz0C + div(wz0CV p ) = 0, ∂t
(2.10)
∂w + div(wV f + hz0CVs ) = −2vl (1 − C p ), ∂t ( [ )] ) ( w2 3 1 + C(η − 1)μ (C) K∇p + Bu 1 + z0 − 1 e2 , Vf = − 12 ) ( 1 3 −1 , K = 1 + z0 μ (C)
(2.11)
(2.12)
(2.13)
] ( [ 2z30 w2 2 ∇p z0 (1 − z0 ) + Vp = − 8z0 3μ (C) ] ) [ ( ) 1 + C(η − 1) 2z30 2 + z0 1 − z0 e2 + Vs + Bu μ (C) 3
(2.14)
To solve Eq. (2.7), we specify the initial and boundary conditions for the particle concentration (1.4): t = 0 : C = 0,
(x, y) ∈ [0, L/H] × [0, 1];
x = 0 : C = C0
The boundary conditions for Eq. (2.8) follow from (1.4): x=0:
12μ (C) ∂p =− , ∂x w2
∂p = 0, ∂x x=
Vol. 49
y ∈ [0, y1 ], [y2 , 1],
∂p L : = −Bu; H ∂y
y = 0, 1 : FLUID DYNAMICS
y ∈ [y1 , y2 ];
No. 2
( ) ∂p = −Bu 1 + C∣y=0, 1 (η − 1) . ∂y 2014
y ∈ [y1 , y2 ].
(2.15)
216
BORONIN, OSIPTSOV
The particle settling velocity is determined from the empirical formula which takes into account the reduction in the velocity with increase in the particle volume fraction. The existing models of effective fluids contain the assumption that the mean-volume velocity of the suspension is described by Poiseuille law, whereas in the present study, using the two-fluid approach based on the balance laws, we show that the Poiseuille law should refer to the mean-mass velocity of the carrier phase (2.9). In early models, the expression for the particle velocity (2.9) contained the mean-volume velocity of the suspension instead of the mean-mass velocity of the liquid phase. As a result, in contrast to the existing models, our two-fluid model contains an additional term −∇(wCVs ) on the right side of the equation for the pressure (2.8), which takes into account the two-velocity effects. We will first consider a one-dimensional suspension flow with effects of particle migration in the fracture but without particle settling, when the dependence of the parameters on the vertical coordinate y can be neglected. We will study the dependence of the parameters on time t and the x-coordinate directed along the fracture. The velocities of the particles and the fluid have only horizontal non-zero components. As above, it is assumed that, due to the transverse migration, the particles form a core of a relative width z0 , with particle-free regions near the walls. The particle and fluid velocities are related by the formulas V p = C0 (z0 , C)V f , [ )][ )]−1 ( ( 3 1 2 2 1 3 −1 1 + z0 −1 . C0 (z0 , C) = 1 + z0 2 3 μ (C) μ (C) We will find the value of particle concentration in the flow core C2 using the inlet value of the concentration averaged over the fracture width C1 . The mass conservation law gives: 1/2 ∫
−1/2
ρ p (C1 )V p (C1 )dz =
z0∫/2
ρ p (C2 (z0 ))V p (C2 (z0 )) dz.
(2.16)
−z0 /2
At the inlet, where the particle concentration is uniform, the average particle velocity is equal to the mean velocity of the carrier phase. At the outlet, where the particle concentration profile is nonuniform, the mean velocity of the particles is related with the mean velocity of the fluid by formula (2.14). We note that in the case considered the particle concentration profile varies along the channel, and hence the coefficient in the dependence of the velocity on the pressure gradient will also vary. From the constancy of the flow rate along the channel for the incompressible steady-state flow, the local longitudinal pressure gradient will also vary. We can find the relation between the pressure gradient at the inlet and outlet sections, with the subscripts “1” and “2” referring to these sections, respectively: [ )]−1 ( 1 − C1 ∂ p2 ∂ p1 1 3 = 1 + z0 . ∂x ∂ x (1 − C2 )μ (C1 ) μ (C2 ) − 1 Using this relation and mass balance law (2.16), we obtain the equation for C2 as a function of z0 , with a dependence on the initial concentration C1 (uniform across the fracture) as a parameter: C2 C1 = C0 (C2 , z0 ) . 1 − C1 1 − C2 The differentiation of this algebraic equation over z0 gives a Cauchy problem for the ordinary differential equation for C(z0 ): FLUID DYNAMICS
Vol. 49
No. 2
2014
EFFECTS OF PARTICLE MIGRATION ON SUSPENSION FLOW
[ ][ ] 3A Ab1 − Ba1 dC z0C(Ab2 Ba2 ) −1 = − , + F + dz0 2B B2 (1 − C)B2
217
(2.17)
2 2μ ′ (C) , a2 = z20 2 , 3μ (C) − 1 3μ (C) ) ( μ ′ (C) 1 2 − 1 , b2 = z30 2 , b1 = 3z0 μ (C) μ (C) ) ) ( ( 2 1 − 1 , B = 1 + z30 −1 , A = 1 + z20 3μ (C) μ (C) ( ) 1 C 3A + , F = z0 2B 1 − C (1 − C)2 ) a1 = 2z0
z0 = 1 : C = C1 .
3. THE RESULTS OF NUMERICAL CALCULATIONS Cauchy problem (2.17) was solved using the 4th order Runge–Kutta method. The results of calculations for different values of the initial concentration, uniformly distributed across the fracture, are presented in Fig. 2a. Each curve in the figure shows the concentration evolution in the flow core due to the transverse migration of the particles in the suspension channel flow (Fig. 1a). Recall that these results are obtained under the assumption that in each cross section the particle concentration profile is described by a piecewise constant function shown in Fig. 1b. Here, the distance along the fracture is a parameter varying along the curve in the graph, and the intensity of the lateral migration is additionally characterized by the dependence z0 (x). Figure 2b shows the calculated dependence of the longitudinal pressure gradient ∂ p/∂ x on z0 (the width of the flow core filled with the particles) and on the initial concentration C1 as a parameter. Depending on C1 , the function ∂ p/∂ x(z0 ) can be either monotonically increasing for small values of C1 , and nonmonotonic with increase in the parameter C1 . This effect is explained by the influence of the following two competing mechanisms. With the appearance of particle-free layers near the walls, the pressure gradient required to maintain a constant flow rate in the fracture decreases. At the same time, the flow viscosity in the core filled with the particles increases with a decrease in the core width z0 and increase in the particle concentration. Depending on the relative magnitude of these two effects, the longitudinal pressure gradient can both decrease and increase with a decrease in the core with z0 . To estimate the effects of particle migration and phase velocity slip on the transport and sedimentation of particles in the hydraulic fracture, we performed calculations of two-dimensional system of Eqs. (2.10)–(2.14) on a uniform rectangular grid. The solution of elliptic equation (2.11) was found by a second-order finitedifference method. The resulting system of linear equations with a five-diagonal matrix is solved by the conjugate gradient method with the ILU(2) preconditioner [22], which is a variant of incomplete decomposition on the upper- and lower-triangular matrices. For solving transport equation (2.10), we used an explicit TVD scheme with a flux limiter [24], which is of the second order with respect to spatial coordinates and the first order with respect to time [23]. This method is based on a combination of the second-order Lax–Wendroff schemes in regions without discontinuities and the first-order Euler scheme in the neighborhood of discontinuities. This makes it possible to avoid oscillations near the discontinuity, typical of the Lax–Wendroff scheme, and obtain a solution with the second order accuracy everywhere outside the discontinuities. As the flux limiter, we used a variant of the Superbee limiter [24], which is optimal for the flows with sharp concentration fronts and gives the smallest error in the integral mass conservation law (about 1% on the 300 × 600 grid after 103 time steps). FLUID DYNAMICS
Vol. 49
No. 2
2014
218
BORONIN, OSIPTSOV
∂p ∂x
z0
z0
Fig. 2. Dependence of particle concentration C in the flow core on the core width z0 for different initial values of the concentration uniformly distributed across the fracture thickness (a); dependence of the longitudinal pressure gradient on the core width z0 for different values of the initial concentration (b): (1–4) C1 = 0.1, 0.2, 0.3, 0.4.
As an illustration of the scheme accuracy, we calculated the transport and sedimentation of particles with a uniform transverse concentration profile z0 = 1 in a channel with flat walls, measuring 1.5 m × 16 cm × 5 mm. At the inlet, the constant velocity 9.63 cm/s and the 2% particle volume fraction are specified. The particles of radius 0.075 mm with material density of 2600 kg/m3 were used. The dimensionless parameters for this case are: Bu = 848, η = 2.6. These parameters correspond to the experiment on particle sedimentation conducted in a flat laboratory channel in the Novosibirsk Schlumberger technology center. Figure 3 shows the results of numerical calculations of two-dimensional system of equations (2.10 )–(2.14), compared with the experimental results and with an analytical formula for the sediment layer height at a time instant of 186 s. Figure 3 shows a good agreement between the numerical calculations, analytical theory, and experiment. A slight discrepancy with the experiment near the outlet section is attributable to the difference in the boundary conditions, since in the experiment there was a small obstacle on the channel bottom at the outlet section, which resulted in the accumulation of particles and increase in the height of the sediment layer. The formula for the growth of the sediment layer is obtained using a quasi-one-dimensional approximation under the assumptions that the height is independent of the longitudinal coordinate and the particle volume fraction above the sediment layer is exactly equal to the concentration at the inlet section ( )α − 1 2a2 (ρ p − ρ f )g C C , α = 5, Cmax = 0.65. 1− H =t 9μ0 Cmax Cmax Figure 4a shows the calculations of the one-dimensional problem of particle transport in the fracture with account for the particle migration effect:
∂C ∂ C0 (z0 , C)C + = 0. ∂t ∂x This equation was solved numerically on a uniform grid using an improved explicit scheme with a flux limiter [24]. The figure shows that the particle migration results in the increase in the depth of particle penetration into the fracture on a fixed time interval. Near the leading edge, a zone of a smooth concentration variation from zero to the initial value is formed, which contradicts to the case of absence of migration z0 = 1 when a sharp concentration front is formed. The formation of the transition zone is attributable to nonlinear nature of the transport equation. It is similar to the Buckley–Leverett problem, where the filteredfluid saturation varies smoothly in the vicinity of the fronts due to the nonlinear dependence of the relative permeabilities on the saturation. FLUID DYNAMICS
Vol. 49
No. 2
2014
EFFECTS OF PARTICLE MIGRATION ON SUSPENSION FLOW
219
Fig. 3. Calculated particle concentration distribution in a laboratory slot, measuring 0.16 × 1.5 m: (1) analytical solution for the deposit thickness, (2) experiment. The particle concentration scale in the numerical solution is shown on the right.
Fig. 4. Positions of the particle concentration front in the fracture at the same instant of time, obtained in one-dimensional calculations for different values of the core width z0 : (1–5) z0 = 1, 0.8, 0.6, 0.4, 0.2 (a). The dependence of the coefficients in the particle transport equations on the core width z0 for a fixed value of the particle concentration in the core C = 0.4 (b): (1) K, (2) ρm , (3) Kp , (4) ρmp , (5) ρm /K, (6) ρmp /Kp .
Figure 4b shows the dependences of the coefficients in the expressions for the velocities of particles and the carrier phase (2.12)–(2.14) on the parameter z0 ) ) ( ( 1 3 3 (1 + C(η − 1)) − 1 , ρm = 1 + z0 −1 , K = 1 + z0 μ (C) μ (C) ( )) ( )) ( ( 3 3 2 2(1 + C(η − 1)) 1 + z20 − 1 , ρmp = 1 + z20 −1 . Kp = 2 3μ (C) 2 3μ (C) The graph shows that the relative magnitude of the term describing the effect of gravitational convection (curves 5 and 6) decreases, as the width of the flow core z0 reduces. Figure 5 shows the results of 2D calculations with account of the effect of particle migration, which confirm this conclusion: the gravitational convection, manifested in the sloughing of the leading edge of the slurry, attenuates with decrease in z0 . These calculations correspond to a fixed particle concentration at the inlet section. We can consider an alternative formulation of the problem of suspension flow in the fracture with account of the evolution of FLUID DYNAMICS
Vol. 49
No. 2
2014
220
BORONIN, OSIPTSOV
Fig. 5. Particle concentration distribution in the fracture for z0 = 0.8 (a) and the positions of the particle concentration front as the function of the core width z0 (b). Fracture dimensions: 70 × 70 m, C = 0.4, t = 60 s.
the particle concentration profile due to the particle migration towards the center. In this flow, the initially uniformly distributed particle concentration increases in the flow core, as the core width reduces (as shown above in Fig. 2a). We performed the calculations taking this effect into account and obtained the results similar to those shown in Fig. 5: the particle migration towards the channel center suppresses the sloughing of the leading edge. Summary. A two-fluid model of a suspension flow in a vertical hydraulic fracture is constructed, which takes into account a nonuniform transverse profile of the particle concentration, formed as a result of particle migration to the central plane of the fracture. As compared to the existing models, which do not take into account the effects of particle migration, the two-dimensional equations of the model obtained by the averaging across the fracture thickness contain modified coefficients, which explicitly depend on the width of the flow core occupied by the particles. To estimate the effect of the particle migration on the particle transport and sedimentation, a series of numerical calculations was performed in the framework of one-dimensional and two-dimensional formulations of the problem. It is shown that the particle migration results in particle accumulation near the central plane. With increase in lateral-migration intensity, the depth of penetration of the particles into the fracture increases, as compared to the case of uniform particle distribution, whereas the effect of gravitational convection (sloughing) in the vicinity of the leading edge reduces. Two-dimensional numerical calculations for the case of a uniform particle concentration profile showed a good agreement with the available experimental results and the one-dimensional analytical theory for the growth of a sediment layer on the fracture bottom. The authors are grateful to M.N. Bulova and S.M. Makarychev-Mikhailov for providing the experimental results. The authors express their gratitude to the management of Schlumberger for the opportunity to publish this study. REFERENCES 1. M.J. Economides and K.G. Nolte, Reservoir Stimulation (John Wiley and Sons, New York, 2000). 2. A.E. Boycott, “Sedimentation of Blood Corpuscles,” Nature 104, 532 (1920). 3. W.D. Hill, R.R. Rothfus, and Kun Li, “Boundary-Enhanced Sedimentation due to Settling Convection,” Int. J. Multiphase Flow 3 (6), 567–583 (1977). 4. A. Acrivos and E. Herbolzheimer, “Enhanced Sedimentation in Settling Tanks with Inclined Walls,” J. Fluid Mech. 92, 435–457 (1979). 5. S.J. McCaffery, L. Elliott, and D.B. Ingham, Enhanced Sedimentation in Inclined Fracture Channels. Topics in Engineering 32 (CD-ROM ISBN 1-85312-546-6) (1997). FLUID DYNAMICS
Vol. 49
No. 2
2014
EFFECTS OF PARTICLE MIGRATION ON SUSPENSION FLOW
221
6. Y.A. Nevskii and A.N. Osiptsov, “On the Role of Non-Stationary and ‘Memory’ Forces in Problems of Gravitational Convection of Suspensions,” Bulletin of Moscow State University. Ser. 1. Mathematics and Mechanics (4), 37–40 (2008). 7. Y.A. Nevskii and A.N. Osiptsov, “Modelling Gravitational Convection in Suspensions,” Technical Physics Letters 35 (7), 98–105 (2009). 8. Y.A. Nevskii and A.N. Osiptsov, “Slow Gravitational Convection of Disperse Systems in Domains with Inclined Boundaries,” Fluid Dynamics 46 (2), 225–239 (2011). 9. J.P.A. Pearson, “On Suspension Transport in a Fracture: Framework for a Global Model,” J. Non-Newtonian Fluid Mech. 1994. V. 54. P. 503–513. 10. P.S. Hammond, “Settling and Slumping in a Newtonian Slurry, and Implications for Proppant Placement During Hydraulic Fracturing of Gas Wells,” Chem. Eng. Sci. 50 (20), 3247–3260 (1995). 11. H. Gu and E. Siebrits, “On the Numerical Solution of Hyperbolic Proppant Transport Problems,” in: Proc. 10th Intern. Conf. Hyperbolic Problems: Theory, Numerics, and Applications (Osaka, Japan, September 13–17, 2004). 12. E.S. Asmolov and A.A. Osiptsov, “The Inertial Lift on a Spherical Particle Settling in a Horizontal Viscous Flow through a Vertical Slot”, Physics of Fluids, 21, 063301 (2009). 13. A.A. Osiptsov and E.S. Asmolov, “Asymptotic Model of the Inertial Migration of Particles in a Dilute Suspension Flow through the Entry Region of a Channel,” Physics of Fluids, 20, 123301 (2008). 14. E.S. Asmolov, N.A. Lebedeva, and A.A. Osiptsov, “Inertial Migration of Sedimenting Particles in a Suspension Flow through a Hele–Shaw Cell,” Fluid Dynamics 44 (3), 405–418 (2009). 15. D. Leighton and A. Acrivos, “Shear-Induced Migration of Particles in Concentrated Suspensions,” J. Fluid Mech. 181 (1), 415–439 (1987). 16. S.A. Boronin and A.A. Osiptsov, “A Two-Fluid Model of Slurry Flow in a Hydraulic Fracture,” Dokl. Ross. Acad. Nauk 431 (6), 1–4 (2010). 17. R.I. Nigmatulin, “Dynamics of Multiphase Media V. 1,2,” Hemisphere, New York, 1990. ¨ 18. W.H. Herschel and R. Bulkley, “Uber die V˙iskosita¨t von Solen,” Amer. Soc. Testing Mater. 26, 621–633 (1926). 19. S.A. Boronin, “Investigation of the Stability of a Plane-Channel Suspension Flow with Account for a Finite Particle Volume Fraction,” Fluid Dynamics 43 (6), 873–884 (2008). 20. M. Van Dyke, Perturbation Methods in Fluid Mechanics (Academic Press, New York, 1974). 21. D.G. Thomas, “Transport Characteristics of Suspensions: VIII. A Note on Viscosity of Newtonian Suspensions of Uniform Spherical Particles,” J. Colloid. Sci. 20, 267–277 (1965). 22. Y. Saad, Iterative Methods for Sparse Linear Systems (2nd Edition, SIAM. 2003). 23. C. Hirsch, Numerical Computation of External and Internal Flows (2nd Ed.). The Fundamentals of Computational Fluid Dynamics, (Elsevier, 2007, ISBN: 978-0-7506-6594-0). 24. R. LeVeque, “High-Resolution Conservative Algorithms for Advection in Incompressible Flow,” SIAM J. Numer. Anal. 33 (2), 627–665 (1996).
FLUID DYNAMICS
Vol. 49
No. 2
2014
