A study of squeezing flow in fracture channels

0148-9062/93 $6.00 + 0.00 Pergamon Press Ltd

Int. J. Rock Mech. Min. Sci. & Geomech. Abstr. Vol.30, No.7, pp. 841-844, 1993

Printed in Great Britain

A Study of Squeezing Flow in Fracture Channels D. LIN*

J.-C. R O E G I E R S *

INTRODUCTION The Fracturing Fluid Characterization Facility (FFCF) is an experimental facility under construction to study the behavior and effects of fracturing fluids and slurries during and after the hydraulic fracturing process. It is envisioned that a high-pressure simulator (Hele-Shaw analog) can provide the means for determining critical parameters under nearly fullscale operating conditions. The parallel plate model is the basis of a mathematical analysis described in this paper. Fracturing fluids contain significant amounts of high-molecular-weight polymers which cause them to be viscoelastic and exhibit many different forms of time-dependent non-Newtonian fluid behavior. The two parallel plates move normal to each other with a slow, time-dependent speed. The velocity field developed in the fracture itself is approximately that of Poiseuille flow. This paper is mainly concerned with the stress (forces) generated in a fracturing fluid, which separates two solid rock surfaces in relative motion. The constitutive equation used is the one proposed by Harnoy (1978) . Unlike the inelastic power-law model, the Harnoy fluid exhibits stress relaxation and normal-stress difference in simple shear flows. Solutions to two basic problems are presented: 1) parallel surfaces with constant approach velocity; and, 2) parallel surfaces with constant load. CONSTITUTIVE

T+A DT Dt = 2 p E ( u )

(2)

where t is time, and w is the rotation tensor of the fluid particle. All theses derivatives are frameindifferent. The Jaumann derivative has a = 0. An upper convected Maxwell model arises from (2) when a - 4-1, while a -- - 1 corresponds to the lower convected Maxwell model and a=0 to the co-rotational model. In elastic flow, the Deborah number De = A/•-, which is the characteristic relaxation time for the fluid A divided by the characteristic time r of the flow. The smaller the De number, the less prominent the elastic effect of the flow. For low De numbers the constitutive equations of viscoelastic fluid may be simplified. For this purpose the deviated stress tensor T from equation (1) can be written in the following operatoric form as a function of E, the rate of strain tensor: T = 2U[1 + A~--t]-lE(u)

(3)

Expanding the operator on the right side of equation (3) to an infinite series in increasing orders of time derivatives of E, and considering only the first two terms, owing to the fact that the fluid time parameter A is very much less than unity, T

EQUATION

It is well-known that the linear viscoelastic equations hold only in the case of small strain rates and stresses. The generalized Maxwell model is (see Tanner 1988):

=

2p[E(u) - AD-~E(u)]

(4)

where p is the dynamic viscosity. Cartesian tensors are adequate for channel flow with low curvature, and are used here. In a rectangular coordinate system, one obtains that:

(a)

where T is the deviatoric stress tensor, A is a relaxation time, # is the viscosity of the fluid, E(u) stands for the strain rate tensor, and D / D t is an invariant derivative given by: *School of P e t r o l e u m & Geological Engineering, The University of Oklahoma, N o r m a n , O k l a h o m a 73019-0628. ~ws ~:7-a

DT OT + (u. V ) T + Tw - w T - a ( E T + T E ) Dt Ot

841

coij

-

D

=

2u(Eij

=

D - P 6 i j + 21a(Eij - A'-~~Eij)

(5)

(6)

where Sij is the total stress tensor, 6ij is the Kronecker delta, and P is the pressure.

842

ROCK MECHANICS IN THE 1990s

In order to separate the stress relaxation from the normal-stre6s effects, the Harnoy material time derivative is used. The time derivative is defined in a rigid rectangular co-ordinate system. The origin is attached to the fluid particle and moves with it, but the directions always coincide with the three principal axes of the rate-of strain tensor. Thus DE OE = + (u. V ) E + Ef~ - f~E Dt Ot

(7)

. Oui

Oui. •

OP - ~

OTq

+ Oxi

(8)

where p is the fluid density, Tq is the stress component, ui is the velocity, and P is the pressure. Narrowing the problem to two-dimensional flow in a narrow gap between two solid surfaces, let L and h denote the reference channel length and gap height, respectively; h0 being the initial gap height at time t = 0. The gapwise velocity is V = dh/dt with the initial gapwise velocity Vo = dho/dt, The narrow gap geometry shows that L > > h. The in-plane velocities are much larger than the gapwise velocity and the primary deformation is shearing across the thickness of the gap. The two-dimensional momentum equations are obtained by substituting the constitutive equations (5,6) into equation (8); hence, Ou

OP

Ou

o(~ + u~ + .~)

v* = v Vono .,~/2a e = P-5;--'

(13) The non-dimensional parameter Re is a Reynolds number and De is a non-dimensional elasticity parameter called the Deborah number. Dropping the star notation for convenience, the following dimensionless governing equations are obtained 02

Ou

?)u

?)P

,)P .L?y .....

-

02u

Ox + 8y--7 0

(14)

together with the dimensionless continuity equation Ou Ov 0-~ + ~ = 0

(15)

In order to study the motion which results from squeezing the fracture flow between two parallel plane surfaces in relative motion, particular attention was first given to a special case where one of the surfaces is fixed and the other is moving with a constant approach velocity or with constant load. The boundary conditions are 1) Entrance pressure condition: P(0, y) = 0

(16)

2) Exit pressure condition: OP(1, y) = 0 Ox

(17)

3) No-slip velocity conditions at the lower surface: u(x,O)

=

v(x, 0) = 0

(18)

4) Velocity conditions at the upper moving surface: (19)

-o---; + ~'-g~y~-

Rearranging the above equation (9) and considering the continuity condition, the equation governing the viscoelastic flow along a long channel becomes Ou

Ou

0P

Oy Ou Ov ~x't'~y

OP

=

0

02u

(10)

0

It should be noticed that for a Newtonian fluid. i.e, De=O, the problem reduces to the conventional Navier-Stokes solution. PERTURBATION SOLUTIONS The first step in a perturbation analysis is to identify the perturbation quantity. This is done by expressing the basic equation (15) in a dimensionless form and, assessing the order of magnitude of different terms and identifying the terms that are small. In the case of a low De and Re situation, asymptotic results can be found, whence one can write u

=

uo + DeoUd +.Reour + O(De~, DeoReo, Re2o)

(20)

v

=

vo + Deovd + Reovr + O(De~, Deoneo, ae2o)

(21)

Po + DeoPd + ReoP~ + O(De~, DeoReo, Re~)

(22)

(11)

Using the following dimensionless quantities: x = -L'

Ou

OH h u(x, H) = 0. v = --~-, If = ho

02u

,, OSu 03u t~atu~ + v-- azoy Oya Ou 02u Ou O:u . (o) Ox ---~ Oy + Oy O-~y )

z"

vbA De =

(Des-j,~ + R e ) ( ~ + u ~ + v-~;,i~) -

where f~ is the rotation tensor of the principal axes of E. The time derivative is objective, namely indifferent under arbitrary rotation of the frame of reference coordinates. Considering steady state flow of an incompressible viscoelastic fluid through a channel; in absence of body forces, the equation of motion is:

P ( W + uj ~ - ~ =

hou u* = VoL'

y ph3 Y* = -~o' P" = -pVoL '---f'

t*

Vot = -~-0' (12)

P

=

ROCK MECHANICS IN THE 1990s Next, this series solution is substituted into the governing equations for the problem. By equating the coefficients of each power of De and Re to zero, one can generate a sequence of subproblems. The problem is solved in succession to obtain the unknown coefficients of the series solution. The zeroth-order term is the classical lubrication approximation:

0 (

H

30Po 120H N

0 (H 3 OPal ) =

0

1/3)

parallel surfaces are separating, H = 1 + t. The pressure represents a correction to the lubrication theory which can be expressed a perturbation as follows

P = P,+DeoPd+neoP~+O(De2o, DeoReo, Re~o) (32) where 1) the lubrication part is represented by

(23)

In plane squeezing flow, non-Newtonian effects alter the pressure field

843

P~=T

[6(1 - (1 - x)2)] H3

(33)

2) the Deborah number effect is represented by p~ =

[72(1 - (1 - x)2)] 5H 3

(24)

(34)

3) the Reynolds number effect is represented by where 302P° 1H4OH(cqP°)2 ~=H 0--~-~+4 Ox Ox

p~ =

[51(1 - (1 - x)2)] 35H 3

(25)

(35)

The pressure can be found as

s = H

OPo O~Po 2

(26) P

Inertia effects on the pressure field can be estimated as follows: a ( 3OP,

1 a [H2(a+lfl)]

=

2 2 ]-~[H-~(H- )] } Pd

=

-3[1-

P~

= 3[1-

d

( 1 - x)2][

Ret

=

Det

=

Re0(1 ± t ) 1 DeO ±~---t l

(37) (38)

The dimensionless load, W, can be scaled as W( 3

W =

pv0

(39)

(28) (H-Z)]

(29)

where ~ is the ration, Lh--~o , and w is the real load. Integrating P(x, t) in eq.(33) the dimensionless load, W, is obtained:

=

//

=

4(1~

=

f01 Pl(x,t) = =F4(1 ±1 t) s

(1- x)2)]H2{d~¢H-2)-

dt 2

W

Hence, by use of the perturbation method, the solution for viscoelastic fluid problems are decomposed into viscous and elastic components. C o n s t a n t a p p r o a c h v e l o c i t y case For the special case of a fixed lower surface and an upper surface which moves impulsively from rest to a constant velocity, the length of this transition period has been discussed by Jones and Wilson (1975) in the case of Newtonian flow. If the parallel surfaces approach or recede with a constant squeezing rate: H = 1± t

(36)

where

3 1 1 - ( 1 - x ) 2 ) ] { d ( H -2) 3

T12De~ + 17 Re,) 4U

(27)

The resulting pressure expressions, subject to the boundary conditions, can be integrated to obtain: P0

6 [(1 - (1 - x)2)] ( + -

(31)

When the upper surface is descending, H = 1 - t , the two surfaces are receding from one another. When two

Wl

P(x,t)

(40)

17 ( 41 ) (q:l - ~ Det + -~--~Ret) (42)

This last load solution, W, is a correction to the classic lubrication results, Wz. At zero Deborah number the pressure or load solutions reduce to the Newtonian case. Figure 1 shows that the inertia always increases the load whereas elasticity always decreases the load. Furthermore, one notes that the coefficient of the Deborah number term in the expressions for the dimensionless loads, is much larger than the corresponding coefficients of the Reynolds terms. Thus one expects that elasticity soon starts to dominate the flow and that the perturbation solution can only be valid for very small Deborah numbers.

844

ROCK MECHANICS IN THE 19905

1.o

. . . . .

-

.....

iii

0.8

•

-

-

,

"-,IV . . . .

. . . . . . . . .

.

~'11.

...................................................:= . . . . . .

II ~ 1

L .................~ ' : t ................ ,~. . . . . . . .

........

,v

'"

,,,

I--

,,I

Z nl

"'~ 0.6

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

I: R e O = 0 . ! , DeO=O.O iiiReO=b~,70eO=d:0 " .III: R e 0 = 0 , 1 , D e 0 = 0 . 0 2 IV: 11o0=0,1., D e O = 0 . 0 4

0.4

\~,

•

...................... :

,, . ~'~

C-

3

i

,

2

"

r~ 0.2

.

~:

I'*

'~

I: ReO=,O.O i D e O = 0 . 4 0 r { II: R e 0 = ~ . 0 . D e 0 = 0 . 2 0 ,~ Ill:Re0=0.2 D e O = 0 . 2 0 .~ / IV: Re0.=-0.f De0=0.O0/i~/~

~ ~t'

i

, i

J~' '

................................................................... . . . . . . . . . . . . . . . . .

o

0

0.25

0.50 LOAD

0.75

I

k.£,"

'.

1.00

(

~. 1.25

RATIO K

0 0.95

:

.

,.

1~o5

I.oo CHANNEL

THICKNESS

1.!o

115

R A T I O K1

Figure 1: Load ratio k = W/Wt in constant approach velocity case

Figure 2: Channel gaswise ratio, kl = H/Ht, in constant load, W, case

Constant load case

of load reduction. Moreover, the analytical results show that descent time is an increasing function of the Deborah number.

If a constant normal force is exerted on the upper surface during the motion we use a new variable y related to the dimensionless height H, y = 1/H ~. In the case of a low De situation, asymptotic results can be found, whence one can write

Y = Yl + Deoyd + Reoy~ + O(De~, DeoReo, Re2o) (43) where 1) the lubrication part is given by: 1 Yt = ~ W t + 1

(44)

2) the Deborah number effect is given by: Yd = -- 3 1 n ( y t ) W

(45)

3) the Reynolds number effect is given by: 1W( 1 Yr = ~-~ 1--'~t)

(46)

In Figure 2, a channel gapwise ratio (H/H~) versus dimensionless time W t / 2 is plotted. From case I (P~o = O.O,Deo = 0.40) and case IV (Reo = 0.1, Deo = 0.00), we found that the inertia effect decreases the descent time while the elasticity serves to increase the descent time. Furthermore, one notes that the effect of the Deborah number term in the expressions for the dimensionless channel gapwise will increase with time and the corresponding effect of the Reynolds terms decreases with time. Further research is required to compare theoretical predictions with experimental results more precisely. Nevertheless, it is poasible to establish that this simple model predicts correctly the trends of the experiments

CONCLUSION This paper proposes an analysis to determine the interaction between boundary exerting forces and viscoelastic fracture flow based on the continuum approach of viscoelastic squeeze flow. The Harnoy model is applied to unsteady (squeeze) fracture flows, including fluid inertia forces. A perturbation solution for a low l~eynolds number and a low Deborah number is obtained. The goal is to study stress relaxation effects in squeeze flow where the fracture thickness varies with time; and to show the difference between Newtonian and non-Newtonian flows. It has been shown that the stress-relaxation effect, which is described by the relaxation time, can reduce normal stress (or load) while the two fracture surfaces approach each other with constant velocity and increase the descent time when a constant load is applied. REFERENCES 1. Haznoy A., An analysis of stress relaxation in e!asticoviscous fluid lubrication of journal bearings. Trans. ASME, vol.lO0, ~87-~95 (1978)

2. Tanner R. L, ENGINEERING RHEOLOGY, O z . ford Engineering Science Seriesl~, Revised Edition (1988). 3. Coleman B.D. and Noll W,, Foundations of linear: viscoelasticity. Rev.Mod.Phys. 331 ~39-~9 (1961).