- Email: [email protected]
PARFES—PARallel Finite Element Solvers for flow-induced fracture
Computing Systems in EngineeringVol. 3, Nos 1-4, pp. 379-392, 1992 Printed in Great Britain.
0956-0521/92$5.00+ 0.00 © 1992PergamonPress Ltd
PARFES--PARALLEL FINITE ELEMENT SOLVERS FOR FLOW-INDUCED FRACTURE AMAURY F. JUNIOR~and M. F. AHMAD~ tDepartment of Mechanical Engineering, MIT, Cambridge, MA 02139, U.S.A. :~National Center for Supercomputing Applications--NCSA, University of Illinois, Urbana-Champaign, Urbana, IL 61801, U.S.A. AMtract--The analysis of flow-induced fracture (also known as hydraulic fracture) propagation is important in various applications, such as petroleum recovery [M. P. Cleary, SPE Distinguished Author Series--Journal of Petroleum Technology 13-21 (January 1988)], nuclear waste disposal, and any process that involves the separation of solid surfaces by the motion of fluids between them. The present work addresses a particular problem in hydraulic fracture, i.e. characterization of gravity-driven motion of multiple stages of immisciblefluids within a narrow fracture cavity. Among capabilities developed to study the present problem [M. P. Cleary and Amaury F. Junior, SPE Journal 24825 (October, 1992)], a suite of numerical algorithms, named PARFES (acronym for PARallel Finite Element Solvers)--designed specially to take advantage of highly parallel computer environments--is introduced here. PARFES is composed of three modules: PARFESl--tracks the interface motion and mesh nodal distribution of a given fluid stage; PARFESAX--models the axisymmetrical multiple stage flow problem; PARFES2--a nonlinear Newton-Raphson algorithm to determine nonsymmetrical motions. Current designs for (expensive)commercial oil and gas fracture treatments tend to greatly overestimate the effective area filled by proppant particles (carried by fluid stages) within the cavity, consequently obtaining less than desirable post-fracture-treatment oil/gas production rates. The PARFES algorithms help to better track the movement of proppant-laden fluid stages within the fracture cavity [Amaury F. Junior, Ph.D. Thesis, MIT (1992)], as well as fluid losses to the adjacent strata.
HYDRAULIC FRACTUREMODEL
(displacement discontinuities) to the change in the three-dimensional stress field. The change in the stress Detailed mathematical modeling of hydraulic fracfield involves distinct contributions: from the fluid(s) ture requires a collection of submodels: partial differpressure field (within the fracture cavity); from the ential equations (PDEs), e.g. describing fluid confining stresses; and from the backstresses (due to transport; integral equations, e.g. elastic displacement fluid losses to the porous matrix strata). and fracture mechanics; history-dependent equations, (4) Stress Intensity Factor--K~c, the critical stress e.g. fluid losses to porous strata; a variety of boundintensity factor in mode I, serves as an integral ary conditions (BCs) for flow and fracture perimeters; constraint to balance the local distribution on the fluid rheologies; and equations of state. A model of wetted (e.g. excess stress field) and nonwetted (e.g. a hydraulic fracture (sketched in Fig. 1) must include confining stresses) regions of the fracture cavity. The the following:1'2'3 effect of incorporating the Klc criterion into the (1) Local and Global Volume Conservation--the model is to determine the ratio RF(x, y, t)/RT(x, y, t), volume of injected fluid stages equals the volume of where RF(X,y,t ) represents the perimeter of the the channel, delimited, in the transverse direction, by outermost fluid stage. Both integral equations, dethe crack opening, + ~(x, y, t)/2, and, in the fracture scribed in (3) and here, have singular integral kernels, plane, by the effective wellbore radius Rw, and the and special numerical integration techniques are recrack tip perimeter, RT(x,y, t). A, co, L, and H are, quired to solve them. respectively, the characteristic scales for crack open(5) Fluid Losses--the effective permeability (coning, nonwetted crack-tip zone, length and height, nected porosity) of the stratum(a) involving the fracwith typical ratios: A/to ~ 1 and A/L ~ A H ~ 10 -5. ture, and a variety of conditions (e.g. viscosities of The global mass conservation provides a constraint fluids inside and outside the fracture, BCs at the for all possible configurations of the fluids interface fracture walls, etc.) determine the amount of fluid and crack front at any given time. losses through the fracture walls (transverse direc(2) Momentum Conservation--coupling of the lo- tion--z plane) to the adjacent stratum(a). cal volume conservation with the momentum conser(6) Boundary Conditions--BCs are, in general, eivation provides a set of PDEs that characterize the ther the weakest or the strongest point of a fracture fluid stages motion within a planar three-dimensional simulator. A typical example is the uncertainty about fracture cavity. an adequate boundary condition for the near-the(3) Elastic Displacement--a set of integral crack-tip region (which has been under scrutiny for equations relating the crack opening displacement many years). For instance, in earlier models, 4 the 379
380
AMAURYF. JUNIOR and M. F. AHMAD
t • '
".' ',
~
""
"~4
• ',,"
, .',~,,x,~. "~,, ~
• "
• ,';
....
e.~ ~,,,~ ~ , ~ o ~ o ~
7~TT~T7777~Y777
x
x
z I
I Non-Wetted Crack Tip Zone I
............t~ y
Fig. 1. Sketch of hydraulic fracture growth driven by multiple fluid stages. excess stress field was assumed to be zero at the fluid front, and to be the (negative) value of the confining stresses within the nonwetted region. Recently, that assumption has been shown to lead to unacceptable results.5 (7) Fluid Rheology (and Equations of State)--The most common fluid rheology equations are the Newtonian linear approximation and the power-law, for weakly non-Newtonian fluids. Although the power-law equation is often used for any fluid that presents a non-Newtonian behavior, it is doubtful, for instance, that it could adequately simulate the flow of some more complex gels (and associated breaker agent effects). Nevertheless, such a simplified fluid model is adopted here for first-order characterization and, especially, for comparison with labora-
tory experiments 2 (where such simple fluids are employed). NUMERICAL MODEL AND PARALLELIZATION ISSUES
Although numerical and experimental simulations of hydraulic fracture have been an ongoing research topic for several years, little progress has been made in modeling the phenomena under realistic conditions. 1 Also, computationally, the problem is not homogeneous, i.e. the fluid flow phenomena can be discretized into a local data structure, but the elasticity and fracture propagation discretization possess an inherently global data structure. By local data structure it is meant that the immediate effect of perturbations in a certain degree of freedom at a
Parallel finite element solvers for flow-induced fracture
381
,....................................................................... .~.y
t
T " \,,, ~'\
F
g
F
B
vz
F
v2
F
s
T !
h
......\ \ ,, ,,,,
x
I X Fig. 2. A blob of heavier and more viscous fluid (stage 2) immersed in a lighter and less viscous fluid (stage 1), within a channel of gap width fi(x). particular node is localized to the neighboring nodes; whereas in a global data structure the immediate effect of perturbations is nonlocal, i.e. it can directly affect all degrees of freedom at all nodes in the model. 6,7 Equations in items (1) and (2) above render good numerical discretization with local data structures algorithms (e.g. Finite Element Method with low order interpolation base), but the coupling of the stress field with the crack opening, see items (3) and (4), requires a global data structure algorithm (there are ways to overcome this disparity--one possibility is presented later in the description of PARFES2). 30.0 27.5 25.0 22.5 20.0 17.5 15.0 12.5
NOON.
DISPLACEMENTS
10.0
7.5 5.0 2.5
S.Tlm
0.0
-2.5 -5.0 -7.5
30.0 27.5
25.O 22.5 2O.O 17.5
M~ i lf4.72 aA4 1LII 10J8 8J~0 7,O4O 4A~ 0~4Q
NODAL DISPLACEMENTS mm i
14.'nz 1=t,44
Fig. 3. Balance of forces acting on fluid stage 2.
Difficulties associated with the numerical solution of hydraulic fracture models are related to the unique physical behavior of the system. For instance: (a) change of fluid pressure, within the crack, at any point affects the displacement field (crack opening) everywhere (global data structure); (b) the integral constraint, K]c, is satisfied continuously in time, so pressure must change throughout the crack when a local perturbation is introduced; (c) one free external moving boundary, the perimeter of the fracture, and one partially constrained internal moving boundary, the perimeter of the fluid front, are simultaneously controlled by mass conservation, fluid(s) rheology, the critical stress 30,0 27,5 25,0 22,5 20,0 17,5 15,0 12.5 10.0 7.5 5.0
-5.0 -7.5
30.0 27.5 25.0 22.5 20.0
1Lie
17.5
15.0
12.5
10,,11o Moo
10.0
UlO
15.0 7.5
5.0 2.5
0.0 -2.5
7,,o4o s.z,m 4,4oo Mm 1,,am
-,5.O -7.5
Fw
NOOAL DISPLACEMENTS
II,ta i 1lM,b744 '~ ILlql lM OIJlOe0 IOAI ~ 7,1,O700 4O 0,1JI1m40O
-7.5 -5.0 -2.5 NODAL
0.0
2.5
DISPLACEM
5.0
7.5 10.0
ENTS M~
i 14,7'1 111.111 10,,1
12.5 10.0
7.5 5.0 2.5
0.0 -2.5 -5.0
7J~O
O.11110
-7.5
l°:q~.O ;.S ~.0:2'.S o'.0 ~,~~'.0*¢.5 i~.o
l°&o.o-r'.S ~.o ~d.S i;.O ~,~ 5'.0 ¢.5 *.io.o
Fig. 4. Sequence of prolated blob distortions modeled by PARFES1.
382 30.0
27.5
AMAURYF. JUNIOR and M. F. AHMAD
NODALDISPLACEMENTS Ja~
25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
i ~,1=
SOJN ~Jm
30.0
NODALDISPLACEMENTS
25.0 22.5 20.0
15.0
18-72
10.0 7.5 5.0 2.5
UAG i
17.5 15.0 12.5
33.12 97.se 24Aa =lJO
10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
17.5 12.5
1.44o
30.0 27.5
25.0
12.N lOae 7.~o
1,44o
0.0
2.5
5.0
7.5
ls.ee scum 7.1oo
5.0 2.5 0.0 -2.5 -5.0 -7.5
10.0
NODAL DISPLACEMENTS M~
i le.12
22.5 20.0 17.5 15.0 12.5 I 0.0 7.5
--10L(~0.0"--7'.5 ~ ' . 0 "2'.5"" OiO'" 2~5*" 510"" i15 "11).0
18JI4
-7.5 -5.0 -2.5
reJ4
1A40
:11,4o 11L72
0.0 -2.5 -5.0 -7.5
se.r~
4,.rod
MAG i 11.12 U IrtJle 14A41
20.0
zsJo
le,,84 12.11e lo,oe 7JaOO 4Jm
NODAL DISPLACEMENTS
22.5
e,t.4e
-10:~0.0"-7.5"'-5.0''-2'.5"'hi()''215 510'"7'.5""i1).0 27.5
30.0 27.5 25.0
l/Jle
!1.110 llL~ 111,84 12Jle lOae 1.t00 4.110 1.4,,0 SI0* rlg"id.0
1°:~.ii-7'.S ~.0 :L5 010 ~
Fig. 5. Sequence of prolated blob distortions modeled by PARFES1 (larger deformations). intensity factor and elastic behavior of the porous matrix; (d) chemical/kinetic characterization of the fluids, gels and proppants injected into the fracture is quite complex. In addition, fluid losses, to the involving stratum(a), can alter the relative composition of an injected mixture; (e) the convective motion of fluid stages within the fracture can be affected by two instability mechanisms: fingering (in-plane--x and y directions) and encapsulation (transverse--z direction). Although fingering can modify the morphology of the fluid stages' interfaces, encapsulation instabilities are most important in the more common practical circumstances of increasing fluid viscosity in successive stages. Encapsulation instabilities can substantially affect the average velocity of the heavier fluid(s) carrying the proppant. For instance, if most of the proppant is carried directly to the bottom of the fracture, post-fracturing oil/gas production can be significantly reduced.
DESCRIPTIONOFPARFES The algorithms forming PARFES are all based on a Euler-Lagrangian approach, e.g. nodes delineating the injection region are constrained while other nodes are free to move, according to the coupling of flow field variables with the elastic stress field.
The tracking of the interface motion and distortion between the multiple fluid stages is performed by PARFES1. A finite element algorithm with a local data structure is employed, and an elliptical conjugate-gradient-type iterative solver utilized. Matrices are sparse and not assembled. The mesh movement/ stretching assumes a simple elastostatic formulation: T,jj = f ,
in f2
u~=a~ inOQu
Tiinj= ~
in
(1)
c~ t
(2)
Tij = 21 ,~oul.i + 22(uia + uj,i) 30.O 27.5 25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
NODAL DISPrEMENTS
i 48,11.1 4802 40/..1 ~4JS IIJ 178.7 ImsJI lm.II IlIOJD M.:I~ 21,44
I°%o ~;S :d., ;.S of(; ~IS ;'.O ilS''I'd.0 Fig. 6. Necking of a prolated blob modeled by PARFESI.
Parallel finite element solvers for flow-induced fracture
30.0 27.51 25,0: 22.5 20,0 17,5 15,0 12.5 10,0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
NODAL DISPLACEMENTS MAG i
0.t47 0.134 0.1~1 0.1M 0.0N 0,0~ O,0"/O 0087 O.044 0,0111 0.00e
..~ . . . .
~. . . .
i ....
i ....
i ....
i.,
i
30.0 27.5 25,0 22.5 20,0 17,5 15.0 12,5 10.0 7.5 5.o
383
NODAL DISPLACEMENTS MAO
.
,q
2.5 0.0 -2.5
0.110 0.100 O.OOe O.OU O,077 O,O07 0067 O.O44 0,0~ 0.010 OJD06
-5.0
•
-?.5
i
"~%0.0"-7~ -5.0 -2.5 0.0 2.5 5.0 7.5 to.o -
30.0 27.5 25.0 22.5 2O.O 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
NODAL
-
DISPLACEMENTS
UA~ i o.me o.ora 0.071 oaea o.m4 o~
•
0.~1 o.m~ o.oo,e.
A
30.0 27.5 25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.O -7.5
NODAL DISPLACEMENTS
iM~
A
0.071 0,Q84 0.O87 0.081 OO44 O~ O~ 0,,~ 0017 0,010 O~
-1°:%.o :;.S-6'.0 :~.5 o'.o 2'.5 5'.o ¢.5 lo.o Fig. 7. Sequence of oblated b/ob distortions modeled by P A R F E S I .
30.0 27.5 25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
N O D A L D I S P L A C E M ENTS M~ i
o.oa o.m.t o.ms o,.ow 0.04t ¢olo o.m4 O.Olo o,o14 ¢.me
o.oos
30.0 27.5 25.0 22.5 20.0 I 7.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0
NOOALDISPLACEMENTS
-2.5 -5.0 -7.5
-1°:9~.0 :;.S :d.o :2'.5 0'.o ¢.5 5'.0 ;.5 ic;.o 30.0 27.5 25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.O 2.5 0.0 -2.5 -5.0 -7.5
N O D A L D I S P L A C E M ENTS
iMAQ
o.olo o,ot2
Fig. 8. Sequence of oblated blob distortions modeled by P A R F E S I (cont'd),
M~ 0.041 0041
O,,NI 0.01O OJ~l OAOI O,,W~
384
AMAURY F. JUNIOR and M. F. AHMAD
:!:! ::i:i:i:i:i:i:!: :.:.:.:.:.:.:.:
!i???i???~??~?i :i:!:i:i:i:i:i: !:!:!:i:i:i:~:i :i:i:i:i:i:i:i: :.:->:,:.:.:.1
:;:: ... ::::
:::::::::::::: Z,:.:4.;.:.
!!!i
!!!!ii!iiii!i!
::::::::::::::
:?:?:?:?:?:?:):
:::~:~
i:i:i:!:~:~:~:~ .:.:.:.:.:.:.:.
iiii!i!iii~
:i:i:i:i:!:~:i: ,.v,.,,,..,,. :2:::::::::::::
[ Onset of instability I
lack
l
..................
n=g ] ,,,..,,,
1.1.1"2 ".,,,.,
>>>i >>:.>
. . '. . ' .." . .
Fig. 9. Onset of encapsulation instability--lighter fluid migrating to the walls.
= Jab ~ wh + f 2 w ~ ) d ~
(3)
c~fZu n dot = 0
+ r (t]w h + t2w~) d(0D) d~f i b
ulj)]wi,~ h
(4)
where the displacements (u) at the boundaries are + [,t2(u ~,2 + u h2,, )]w h1,2+ [;t, (u h~,~+ u ~,2 ) +
coupled to the flow field and the elastic stress field-notice that the elastostatic (local) formulation above is not related to the more complex elastic stress field
22(u~,2 + u2,2)]w2,2 h h + [22(uh, l + uL2)w2,, h h } df~
Frosted--Wedge V e l o c i t y vs. D e n s i t y 10
•
'
'
'
'
I
.
.
.
.
I
.
.
.
.
I
.
.
.
Ratio
.
I
'
'
,"
'
Velocityattached019 Velocityunanachcd 019 Velocityattached020 Velocityunattached020 Velocityattached 021 Velocityunattached021
O 0
.B
E 0,1
0.01
0.001
, 1.4
i
i
i
I 1.45
i
i
T
i
I 1.5
i
i
i
i
[ 1.55
,
i
i
i
I 1,6
i
~
i
i
I 1.65
i
i
i
i 1.7
Rho2/Rhol Fig. I0. Experimental results for the downwards velocity of a blob of higher density fluid (2) displacing a lighter fluid (I). Vertical Hele-Shaw type cell configuration, one float glass wall, the other a "frosted" float glass wall. Three different fluids used as the heavier fluid.
Parallel finite element solvers for flow-induced fracture
385 err_P'W
y
EIToT 19g+1 ,
~
en'_VF
IE+O
err_P,J"D~
IE-I 1rE..2
err_RT
195.3'
tn'_VOL
1"£-4
IZ-y'
\1
1~.6"
"'--.....z
~f
.........
1~-7'
IE.& IE.9-
~x,
IE.I0" 1E.11"
~.~ ,,.~,~ ~ ~.
\~
'V'
V~l e, ! ,4 v
"
IE-12. 1E.13. IE-14' 1E-15.
Fig. 11. Downwards movement of three immiscible fluid stages contained within a channel of variable height 6(x,y, t).
(global) formulation, except for the coupling at the boundaries. The finite element basis functions w h and their first derivatives are assumed to be square integrable in fl, with support zero at the prescribed displacement boundaries (~3~.). The numerical implementation of PARFESI allows the mesh to wrap around itself (the topology may change from convex to nonconvex), in order to track the development of more complex shapes (for instance, tracking the development of fingering patterns). A modified nonlinear Newton-Raphson algorithm, PARFES2, incorporates the physics governing the multiple fluids flow motion. The nodal movement and the fluid flow velocity fields are coupled to the pressure and crack opening distribution fields.
¢~_q,IW Error
1ff*.16' O~g+O
1E+2
2E+2
3"£÷2
4'F.+2
5E+2
6E÷2
Tune (seconds)
Fig. 13. PARFESAX results for 1000 mesh points and 101 time steps. Numerical error (comparison against analytical results) for: pressure at the wellbore; velocities at the wellbore; velocities at the fluid front; mesh velocities at the fluid front; fluid front position; and total volume conservation. Although the algorithm's data structure follows the same ideas implemented in PARFES1, its results depend on the 3-D elastic field and fracture criterion coupling performed by R 3 D H - - a general fracture simulator developed in a separate project) '8,9 R3DH is used to synchronize the global variations in crack opening distribution with the predicted pressure flow field given by PARFES2. The coupling of the two algorithms would represent a serious computational bottleneck (due to the physics of the problem), if calculations were performed for both algorithms in refined meshes. However, a better use of the computational resources may be obtained by immersing the very refined meshes needed in PARFES2 with somewhat coarse discretizations done in R3DH (the global data structure surface integral scheme employed in
IE+I IE+O
,~2
i
k\
led
....... e~._._P,F err VOL
II II
~\
Error IE+O i
IE-4
_,I
lq£-~
IE.3 \ l IE-4
1~-sl I \ IE.IO ~
!!1
i'
1~-; t
IE.11
IE-8
1~-91
'~-'~,~., v
,
V
I
~
1E.101
v-
.~rrRP/XTt J
1E-11 1
1E-I~;
1E-12 1
1E.1C 1000
2000
3000
4000
5000
T~e
6000
~econds)
Fig. 12. PARFESAX results for 100 mesh points and 101 time steps. Numerical error (comparison against analytical results) for: pressure at the wellbore; velocities at the wellbore; velocities at the fluid front; mesh velocities at the fluid front; fluid front position; and total volume conservation.
' '
1E.13 1
IP,.14 1 I~.15 1 IE.I6 0
100
2flO
300
400
~
.zrr_R~DOT_tl01
I
5-00
600
Tune (seconds) Fig. 14. PARFESAX results for 13 and 101 time steps (1000 mesh points). Numerical error (comparison against analytical results) for mesh velocities at the fluid front.
386
AMAURY F. JUNIOR and M. F. AHMAD - -
e r r ' P W tO06
Wellbore Pressure Error
,
,trot - -
1E+O
@-- Pt~ssure_2-F1uids_T-200
err_'1~V_t600
1.0E+0
!
err_'pW_tI5000
Pressure_5-Fluids_T-200
--~
err_q'gV_t6000
........
Pres~ure_l-Fluid_T-200
I --I~
err_PW_t6
•
\ 1.0E-3
. . . . .
1.0E4
. . . .
\
1.0E-5 1.0E~5
.
.
.
.
1.OE-7
\
1.0E-8 IE-5 1E.4
1E.3
rE-2
IE.1
IE+O
1E+1
1E+2
IE+3
IE*4 ffime
F i g . 15. P A R F E S A X results, 101 t i m e s t e p s , for d i f f e r e n t ( 0 . 0 6 , 6, 6 0 0 , 6 0 0 0 a n d 1 5 , 0 0 0 son against analytical results)
1.0E-9
IE+5 (seconds)
u s i n g 1000 m e s h p o i n t s a n d values of total running time s). N u m e r i c a l e r r o r ( c o m p a r i for p r e s s u r e at t h e w e l l b o r e .
R 3 D H , although computationally intensive, does not require refined meshes, except near the crack tip perimeter). P A R F E S A X contains all the intrinsic physical modeling capabilities o f P A R F E S 2 combined with the mesh motion o f P A R F E S 1 , but it is limited to track the movement o f multiple fluid stages for axisymmetrical cases. In spite o f this simplifying symmetry, the accurate solution o f cases o f interest is also computationally intensive. By comparing the results o f this code with derived analytical solutions o f one fluid stage problem, a great amount o f physical insight has been gained (see next section).
1.0E-10
. . . . . . . .
1,0E+O
i
1.0E+I
. . . . . . . .
i
1.0E+2
. . . . . . . .
I
. . . . . . . .
1.0E+3
I
1.0E+4
Total Number of Mesh Points Fig. 17. PARFESAX error analysis for the final predicted wellbore pressure versus the total number of mesh points. The same boundary conditions and fluid characteristics of the one-fluid stage are imposed on both the two-stage and five-stage problems. Result for 200 time steps (notice that the error is weakly dependent on the number of time steps---compare with Fig. 16). The mass and m o m e n t u m conservation formulation (assuming an aspect ratio o f order A / L ~ A / H ~ 10-s), for both P A R F E S 2 and P A R F E S A X , can be expressed as
I
J ~(t)
[(ip'v a),, + ('piv ,~%),k + 2pv ' ~ ' V,o~] '~j d f~ = 0 V '~kj e L2['n(t)]
(5)
Wellbore Pressure Error
Wellbore Pressure Error
l-
Pressure_l -Fluid_R-50
Pressure_l-Fluid T-50
Pressure_2-Fluids R-50
Pressure 2-Fluids_T-50 t .OE+O
~lk
Pressure_5-FluidsT-50
1.0E+O
1.0E-1
1.0E-I
1.0E-2
1.0E-2
1,0E-3
1.0E-3
l .OE4
1.0E4
Prossure_5-Fluids R-50
1.0E-5
1.0E-5
1,0E~6
1.0E-6
1.0E-7
I.OE-7 1.0E-8
l.OE-8
1,0E-9
1.0E-9 1.0E-10
1.0E-10 1.0E+0
1.0E+I
1.0E+2
1.0E+3
1.0E+4
Total Number of Mesh Points Fig. 16. PARFESAX error analysis for the final predicted wellbore pressure versus the total number of mesh points. The same boundary conditions and fluid characteristics of the one-fluid stage are imposed on both the two-stage and five-stage problems. Result for 50 time steps.
1.0E+O
1.0E+l
1.0E+2
1.0E+3
1.0E+4
Total Number of Time Steps Fig. 18. PARFESAX error analysis for the final predicted wellbore pressure versus the total number of time steps. The same boundary conditions and fluid characteristics of the one-fluid stage are imposed on both the two-stage and five-stage problems. Result for 50 mesh points.
Parallel finite element solvers for flow-induced fracture
£n,,~[% - ~ 62 [- (',,p),,+'p'vg,]]%d'fl= o V 'gpjE H~['n(t)]
(6)
where for each fluid stage (prescript i): ~v is the volume fraction; *p the density; ~/~the viscosity; % the average (in-plane) velocity component in the k direction; ~V~o~the (transverse) fluid loss velocity (through both fracture walls); g~ the gravity acceleration component in the k direction. The crack opening field 8 is coupled to the fluid(s) pressure field p by means of a surface integral scheme, s encompassing the elastic stress field around the fracture walls, and the stress intensity factor in mode I. The finite element basis functions q~ and ~ are assumed to be square integrable in ~, for any instant of time t. The first derivative of q~is also assumed to be square integrable with support zero at the boundaries. Assuming the following notation for integration in [1 and its boundaries Of~ (atlb~) - [o
f~
(up to fifth order) Backward Differentiation schemes are used. For instance, expanding Eq. (8) and applying Newton-Raphson for time steps (superscript n) n ~<2, with 6 = d + A 6 and cr = p + Aa, gives the following incremental equation: i-
i-
A
~
- 6 At 4 [ p d2e
L'~
n
n
aalnk]
axk
J
+ 2 At ~s \7~ a ~-~-XkI ~Xk/ i-
~./
\ I~ oxk
n
ICXkl
i -2
In
--2Atli~PXkA6lnkl ~
a~bmd[)
[allnm] =--[" atnmd(df~),
387
(7)
-6At~s~(~dEgkA6 -
07~\ ~ dxk/
the previous Eqs (5) and (6) can be combined as
~(ff
I'/')-~s~=O
L'*'
~--n~
&*
+ ¢tact~ 63Tg, lnk
+ 2 At('~)?k A 31 dxk/
J =2Ato~I~d3"~nk
- [~t~ aq'£,ln,]
1
i - 2 d3"gklnk In - 2 At "B(~2['~fi
2/i~ 3 da
o g~ ~Xk/) + 2('tS'V,o,,[~) = 0
(8)
where
+ 2 At %~2s
- d3gk Ox,/
- 2 AtQ~df(k ~xg)"+ 2 At['~dtltf(~lnk]" AR
i-
~XS = ~ RR
~n =
pRgLR - -
c9 a~\"
-- 2 At ~s2([~- d3 ~
~xk/
(9)
O"R
are the nondimensionalization variables; a the excess pressure field (a =p-ac); cr~ the confining stress field; -£k the nodal velocity component in the k direction; and ~ the finite element weight function. The terms in "{'karise from passing the time differentiation out of the integral (first term on the L.H.S.) in Eq. (5). For the implementation of PARFES2, two implicit A-stable m linear multistep schemes are used to time integrate Eq. (8). An Adams-Moulton of second order is used for the first time step and a Backward Differentiation, also of second order, for the subsequent time steps. For the axisymmetrical problem (the gravity term is dropped in PARFESAX), after the second time step, implicit adaptive variable order
- 4 At('H'V~ossl~) n +4(,~61~u).-i
3(~Hdlq')"
(,~61q,). 2.
(10)
More details of the above derivations and equations can be found elsewhere. 2 RESULTS AND CONCLUSIONS There are two versions of PARFESI: one for coarse grained and the other for fine grained parallel computing environments. The coarse grained version was tuned to work at close-to-peak performance at the NCSA's CRAY Y-MP. For instance, substantial reduction of the number of instructions per second and a ten-fold decrease on wall-clock time was obtained with respect to the original version (various
388
AMAURYF. JUNIORand M. F. AHMAD
performance issues, such as efficient array indexation and smaller loop unrolling were already considered on the original version). The fine grained version was written in CM-Fortran and tuned for use on the NCSA's CM-2 (Connection Machine 2). The code has been almost fully rewritten to take advantage of the Massively Parallel environment. Experimental analysis of the convection motion (due to gravity-driven forces) of a fluid stage blob moving downwards within a narrow cell (vertical Hele-Shaw type cell, filled with a different fluid), has demonstrated 2'3 the need to characterize in-plane deformations of the blob. An immiscible heavier blob displacing a bath of lighter fluid, sketched in Fig. 2, will sink against the following counteracting forces (see Fig. 3): (i) shear stresses at the walls (no encapsulation) or at the interfaces with the lighter fluid in the transverse direction (if encapsulated); (ii) buoyancy; (iii) shear stresses to displace the lighter fluid; (iv) surface/interfacial tension and surface wetting characteristics of both the displaced and the displacing fluids. PARFES1 is used to characterize the distortion of the blob under different in-plane boundary conditions. For instance, Figs 4 ~ present a sequence of increasing rates of prolation for a distorted blob--the vertical axis representing the direction of motion is reversed (the more elongated part of the blob is moving downwards). Figures 7 and 8 show one sequence of increasing rates of oblation - - t h e motion is towards the horizontal positive direction. The tendency of higher viscosity fluids to migrate to regions of lower shear rates may dramatically alter the wetting conditions at the fracture walls Wellbore Pressure Error
Pressure_l -Fluid_R- 1000 0-
Pressure 2-Fluids_R-1000
~
Pressure_5-Fluids_R-10CO
Wellbore Pressure (x 2x 3x 4x 5x) [ [
2.0E+6
t
1.8E+6 1.6E+6
--m~
Pressure (Variable Flux) S-Fluids
--@~
Pressure (Variable Density) 5-Fluids Pressure (Variable Viscosity) 5-Fluids
1.4E+6 1.2E+6 1.0E+6 8.0E+5 6.0E+5 4.0E+5 2.0E÷5
!
O.0E+O
........ 1.0E+O
........ i 1.0E+I
• •
) ........ 1,0E+2
i
i ........ 1.0E+3
i 1.0E+4
Total Number of Mesh Points
Fig. 20. Convergence analysis for the final predicted wellbore pressure versus the total number of mesh points tracking five fluid stages. The injected flux/density/viscosity are increased as x, 2x, 3x, 4x, and 5x for each stage (x is the flux/density/viscosity for the first fluid stage).
(see sketch in Fig. 9). This (transverse) migration mechanism, previously termed encapsulation, may substantially affect the average velocity of the heavier proppant-laden fluid(s) and, consequently, the placement of proppant particles within the fracture cavity. Although fingering (in-plane instability-modeled by PARFESI), can modify the morphology of the fluid stages' interfaces, encapsulation instabilities are more important in hydraulic fracturing. For Injected Flux (x 2x 3x 4x 5x) Flux (Variabel Flux) 5-Fluids
1.0E+O
2.0E+O --0-
1.OE-1
.................................................................
Flux (Variabel Density) 5-Fluids
1 .SE+O Flux (Variabe[ Viscosity) 5-Fluids
1.0E-2
10E-3
*-i
*
1.0E-4
1.6E+0
................
1.4E+O
:
"
...............................
~ ...................................
'
] .2E+O i
1.0E-s
i
1.0E+O
:i 1.OE-6
!
i
'
8,0E-1
1.0E-7
6.0E-I
1.0E4~
4.0E-1
1.0E-9
::
l O~.~o
........ 1.0E+O
i 1.0g+l
'
........
! I,OE+2
........
! 1.0E+3
2.0E-1 O.OE+O
. . . . . . . . .
1.0E+4
Total Number of Time Steps
Fig. 19. PARFESAX error analysis for the final predicted wellbore pressure versus the total number of time steps. The same boundary conditions and fluid characteristics of the one-fluid stage are imposed on both the two-stage and five-stage problems. Result for 1000 mesh points.
........ LOE+O
i ........ 1.0E+1
t ........ 1.0E+2
I ........ 1.0E+3
i 1.0E+4
Total Number of Mesh Points Fig. 21. Convergence analysis for the final injected flux versus the total number of mesh points tracking five fluid stages. The injected flux/density/viscosity are increased as x, 2x, 3x, 4x and 5x for each stage (x is the flux/density/viscosity for the first fluid stage).
Parallel finite element solvers for flow-induced fracture
389
fluid when displaced upwards is also constrained to move within the fracture boundaries. Another set of semi-analytical algorithms (also described elsewhere 3) were developed to characterize the influence of the fracture boundaries on the blob motion. The tracking of one interface motion was verified against derived analytical solutions using PARFESAX. An error analysis for various physical variables of interest, is shown in Figs 12-15. Multiple interface tracking was verified in two ways: (i) imposing exactly the same boundary conditions and fluid characteristics of the one-fluid stage analytical solution to the various fluid stages--
instance, experiments conducted for different fluids 3 (see Fig. 10) have demonstrated a spread of about two orders of magnitude on the downwards blob velocity. This velocity spread led us to develop an analytical algorithm, described elsewhere) to predict the size of the encapsulated layers, as well as an experimental criterion to predict the onset of encapsulation. The downwards motion of multiple, immiscible, fluid stages (see Fig. 11) is far more complex than the two-stage problem, especially when modeling inplane distortion. Excluding fluid losses through the fracture walls, the various stages compete against each other when filling the fracture, i.e. the lighter
6 4 Speed Up 2 0
Paral lel Regions Iml = Serial t,lork
[]
: Par.allel t,lork
sS
ss SS sS s
8
s
s~
s~ 4
7
f
s sS
-t~ t ¢s
6
/
is
--
J /J
Speed 5
s~
--
Up
s t
4
--
s
i
ev e s
J
/ f sS
3
-
,/ •
2
1
-
2.5
f
o o
..."
.
¢-~(-
2.0
~.5. ~
2,3 . , . . , - .r
2.4 . .
.
,,. =, ~ .
•
,IF ~ ' p ~ "
2.3
.
....
/I.0 I
, I
1
2
.
I
I
I
I
I
I
3
4
5
6
7
8
# CPU$ Fig. 22. PARFESI optimization resultsfor a 100 x I00 mesh, usingthe NCSA's CRAY Y-MP. Compare parallel scale-up when the number of mesh points is increased (see Fig. 23).
390
AMAURYF. JUNIOR and M. F. AHMAD
Figs 16--19 show an error analysis for 1, 2 and 5 fluid stages; (ii) doing a convergence analysis in the general case (no analytical solution available), see Figs 20 and 21. With respect to the suitability of the algorithm(s) to parallelization, optimization results for PARFES 1, using the NCSA's CRAY Y-MP, show a good parallel scale-up when the number of mesh points is increased (see Figs 22-24). The fine grained implementations, using the Connection Machine 2 (at various locations), confirmed our expectations, e.g. we were able to run very large problems which almost 98% utilization of the processes. Two points--one related to the physics of the multiple stage flow problem, the other related to the
computational data structure--are worth mentioning: (i) in order to capture the development of relevant physical phenomena, such as fingering and encapsulation instabilities,3,tH3 of the multiple fluids motion, each full iteration loop of PARFES2 requires very refined meshes (provided by PARFES1), and it has to be done often and computationally cheaply (also, the algorithm should allow the mesh to wrap around itself); and (ii) the localized data structure of the PARFES1 algorithm (matrices are not assembled), the use of iterative solvers, and the unique layout of local matrices (despite of an overhead in memory usage) have proved to be very successful, and those ideas will be extended to some modules of PARFES2.
m
Speed Up
Parallel Region: m : Serial Work [ ] : Parallel Work
O0f
~ummar~of residl:resid1~49
tS jJ S~ O
8
S°
SS S
SS
~S jS
7 --
Ip
s~
sI sS / sS / sS
Speed
5
-
,/
Up
4o6
,,/"
4.3 ~ 3.9 ..~,~1V"~4.2 '
/// 4 --
///
,,/" 3
-
/
3,5 ~ , - ~ ¢ ~
4.6
3o9
Z.O~ "¢'- 3'5 3.O J
,/,~,1,. 1.8
2 --
1.07"
I --/0.o I"O0.o 1
1
0.1
0.0
0.0
0.I
0.1
0.0
. . . . _'A . . . . . . . . . . . ~ . . . . . . . . . . . ~ ......... "1 ........... L... ....... ~.......... ~i
2
:3
4
o CF~Js
5
6
7
8
Fig. 23. PARFESI optimization results for a 150 x 150 mesh, using the NCSA's CRAY Y-MP.
Parallel finite element solvers for flow-induced fracture
391
]
SpeedUp Parallel Re9ior~ m
= Serial Work
restdl: _
[]
[] = Parallel IJork
r e s l d l Q 4 9 - 2000
8.0
[]
keyshell
7.9
Ke~
--
Speed
7.0
andahls
---
overhead dedicated maximu~ variable ~Ini~u~ variable linear speed up pro9ram speed up
.....
6.9 6.0
-
-
5.9
...... ......
5.0 4.9
-
Up
4.~4.0 3.% 2.7 2oOy 1,4" 2.0
2"
1.0 0.3 0.1
1
2
0.1
3
4 5 G 7 8 w CPU$ Fig. 24. PARFESI optimization results for a 150 x 150 mesh, using the NCSA's CRAY Y-MP. Result for the most computationally intensive subroutine. The program presents a good parallel scale-up (even in its coarse grained implementation) when the number of mesh points is increased.
Acknowledgements--The authors would like to thank: the Gas Research Institute for supporting this research project at MIT; the National Center for Supercomputing Applications for extensive use of their computers; the Advanced Computing Research Facility at the Argonne National Laboratory for use of their computers; and the Thinking Machines Corporation for use of the CMNS. The following MIT undergraduate students, under our supervision, have contributed to the work described here: Yitwah Cheung, Adeeb Shana'a and Lisa Zhang. Thanks for their many valuable suggestions and help.
REFERENCES 1. M. P. Cleary, "The engineering of hydraulic fractures-State of the art and technology of the future," SPE Distinguished Author Series--Journal of Petroleum Technology 13-21 (January 1988).
2. M. P. Cleary and Amaury F. Junior, "Proppant convection and encapsulation in hydraulic fracturing: practical implications of computee and laboratory simulations,'* SPE Journal 24825 (October, 1992). 3. Amaury F. Junior, "Modelling of flow during fracture growth driven by multiple fluid stages using highperformance parallel computing algorithms," Ph.D. Thesis, MIT, 1992. 4. K. Y. Lam, "The development of a fully threedimensional simulator for analysis and design of hydraulic fracturing," Ph.D. Thesis, MIT, 1985. 5. D. Barr, "Leading-edge analysis for correct simulation of interface separation and hydraulic fracturing," Ph.D. Thesis, MIT, 1991. 6. Amaury F. Junior, "An adapted finite element method for massively parallel processors," The Fourth Conference on Hypercubes, Concurrent Computers, and Applications, Monterey, California, pp. 1329-1333, 1989.
392
AMAURY F. JUNIOR and M. F. AHMAD
7. K. K Mathur, "The finite element method on a data parallel architecture," Fifth International Symposium on Numerical Methods in Engineering (edited by R. Gruber, J. Periaux and R. P. Shaw), pp. 599~11, Springer-Verlag, 1989. 8. M. P. Cleary, M. Kavvadas and K. Y. Lam, "Development of a fully three-dimensional simulator for analysis and design of hydraulic fracturing," SPE/DOE Journal 11631 271-282 (1983). 9. M. P. Cleary, D. T. Barr and R. M. Willis, "Enhancement of real-time hydraulic fracturing models with full 3-D simulation," SPE Journal 17713 79-93 (1988). 10. J. Sand and O. Osterby, "Regions of absolute stability,"
Computer Science Department, Aarhus University, Denmark, 1979. 11. N. R. Anturkar, T. C. Papanastasiou and J. O. Wilkes, "Linear stability analysis of multilayer plane Poiseuille flow," Physics of Fluids A 2(4), 530-541 (1990). 12. D. D. Joseph, K. Nguyen and G. S. Beavers, "Nonuniqueness and stability of the configuration of flow of immiscible fluids with different viscosities," Journal of Fluid Mechanics 141, 319-345 (1984). 13. S. J. Weinstein, V. E. B. Dussan and L. H. Ungar, "A theoretical study of two-phase flow through a narrow gap with a moving contact line: viscous fingering in a Hele-Shaw cell," Journal of Fluid Mechanics 221, 53-76 (1990).
0956-0521/92$5.00+ 0.00 © 1992PergamonPress Ltd
PARFES--PARALLEL FINITE ELEMENT SOLVERS FOR FLOW-INDUCED FRACTURE AMAURY F. JUNIOR~and M. F. AHMAD~ tDepartment of Mechanical Engineering, MIT, Cambridge, MA 02139, U.S.A. :~National Center for Supercomputing Applications--NCSA, University of Illinois, Urbana-Champaign, Urbana, IL 61801, U.S.A. AMtract--The analysis of flow-induced fracture (also known as hydraulic fracture) propagation is important in various applications, such as petroleum recovery [M. P. Cleary, SPE Distinguished Author Series--Journal of Petroleum Technology 13-21 (January 1988)], nuclear waste disposal, and any process that involves the separation of solid surfaces by the motion of fluids between them. The present work addresses a particular problem in hydraulic fracture, i.e. characterization of gravity-driven motion of multiple stages of immisciblefluids within a narrow fracture cavity. Among capabilities developed to study the present problem [M. P. Cleary and Amaury F. Junior, SPE Journal 24825 (October, 1992)], a suite of numerical algorithms, named PARFES (acronym for PARallel Finite Element Solvers)--designed specially to take advantage of highly parallel computer environments--is introduced here. PARFES is composed of three modules: PARFESl--tracks the interface motion and mesh nodal distribution of a given fluid stage; PARFESAX--models the axisymmetrical multiple stage flow problem; PARFES2--a nonlinear Newton-Raphson algorithm to determine nonsymmetrical motions. Current designs for (expensive)commercial oil and gas fracture treatments tend to greatly overestimate the effective area filled by proppant particles (carried by fluid stages) within the cavity, consequently obtaining less than desirable post-fracture-treatment oil/gas production rates. The PARFES algorithms help to better track the movement of proppant-laden fluid stages within the fracture cavity [Amaury F. Junior, Ph.D. Thesis, MIT (1992)], as well as fluid losses to the adjacent strata.
HYDRAULIC FRACTUREMODEL
(displacement discontinuities) to the change in the three-dimensional stress field. The change in the stress Detailed mathematical modeling of hydraulic fracfield involves distinct contributions: from the fluid(s) ture requires a collection of submodels: partial differpressure field (within the fracture cavity); from the ential equations (PDEs), e.g. describing fluid confining stresses; and from the backstresses (due to transport; integral equations, e.g. elastic displacement fluid losses to the porous matrix strata). and fracture mechanics; history-dependent equations, (4) Stress Intensity Factor--K~c, the critical stress e.g. fluid losses to porous strata; a variety of boundintensity factor in mode I, serves as an integral ary conditions (BCs) for flow and fracture perimeters; constraint to balance the local distribution on the fluid rheologies; and equations of state. A model of wetted (e.g. excess stress field) and nonwetted (e.g. a hydraulic fracture (sketched in Fig. 1) must include confining stresses) regions of the fracture cavity. The the following:1'2'3 effect of incorporating the Klc criterion into the (1) Local and Global Volume Conservation--the model is to determine the ratio RF(x, y, t)/RT(x, y, t), volume of injected fluid stages equals the volume of where RF(X,y,t ) represents the perimeter of the the channel, delimited, in the transverse direction, by outermost fluid stage. Both integral equations, dethe crack opening, + ~(x, y, t)/2, and, in the fracture scribed in (3) and here, have singular integral kernels, plane, by the effective wellbore radius Rw, and the and special numerical integration techniques are recrack tip perimeter, RT(x,y, t). A, co, L, and H are, quired to solve them. respectively, the characteristic scales for crack open(5) Fluid Losses--the effective permeability (coning, nonwetted crack-tip zone, length and height, nected porosity) of the stratum(a) involving the fracwith typical ratios: A/to ~ 1 and A/L ~ A H ~ 10 -5. ture, and a variety of conditions (e.g. viscosities of The global mass conservation provides a constraint fluids inside and outside the fracture, BCs at the for all possible configurations of the fluids interface fracture walls, etc.) determine the amount of fluid and crack front at any given time. losses through the fracture walls (transverse direc(2) Momentum Conservation--coupling of the lo- tion--z plane) to the adjacent stratum(a). cal volume conservation with the momentum conser(6) Boundary Conditions--BCs are, in general, eivation provides a set of PDEs that characterize the ther the weakest or the strongest point of a fracture fluid stages motion within a planar three-dimensional simulator. A typical example is the uncertainty about fracture cavity. an adequate boundary condition for the near-the(3) Elastic Displacement--a set of integral crack-tip region (which has been under scrutiny for equations relating the crack opening displacement many years). For instance, in earlier models, 4 the 379
380
AMAURYF. JUNIOR and M. F. AHMAD
t • '
".' ',
~
""
"~4
• ',,"
, .',~,,x,~. "~,, ~
• "
• ,';
....
e.~ ~,,,~ ~ , ~ o ~ o ~
7~TT~T7777~Y777
x
x
z I
I Non-Wetted Crack Tip Zone I
............t~ y
Fig. 1. Sketch of hydraulic fracture growth driven by multiple fluid stages. excess stress field was assumed to be zero at the fluid front, and to be the (negative) value of the confining stresses within the nonwetted region. Recently, that assumption has been shown to lead to unacceptable results.5 (7) Fluid Rheology (and Equations of State)--The most common fluid rheology equations are the Newtonian linear approximation and the power-law, for weakly non-Newtonian fluids. Although the power-law equation is often used for any fluid that presents a non-Newtonian behavior, it is doubtful, for instance, that it could adequately simulate the flow of some more complex gels (and associated breaker agent effects). Nevertheless, such a simplified fluid model is adopted here for first-order characterization and, especially, for comparison with labora-
tory experiments 2 (where such simple fluids are employed). NUMERICAL MODEL AND PARALLELIZATION ISSUES
Although numerical and experimental simulations of hydraulic fracture have been an ongoing research topic for several years, little progress has been made in modeling the phenomena under realistic conditions. 1 Also, computationally, the problem is not homogeneous, i.e. the fluid flow phenomena can be discretized into a local data structure, but the elasticity and fracture propagation discretization possess an inherently global data structure. By local data structure it is meant that the immediate effect of perturbations in a certain degree of freedom at a
Parallel finite element solvers for flow-induced fracture
381
,....................................................................... .~.y
t
T " \,,, ~'\
F
g
F
B
vz
F
v2
F
s
T !
h
......\ \ ,, ,,,,
x
I X Fig. 2. A blob of heavier and more viscous fluid (stage 2) immersed in a lighter and less viscous fluid (stage 1), within a channel of gap width fi(x). particular node is localized to the neighboring nodes; whereas in a global data structure the immediate effect of perturbations is nonlocal, i.e. it can directly affect all degrees of freedom at all nodes in the model. 6,7 Equations in items (1) and (2) above render good numerical discretization with local data structures algorithms (e.g. Finite Element Method with low order interpolation base), but the coupling of the stress field with the crack opening, see items (3) and (4), requires a global data structure algorithm (there are ways to overcome this disparity--one possibility is presented later in the description of PARFES2). 30.0 27.5 25.0 22.5 20.0 17.5 15.0 12.5
NOON.
DISPLACEMENTS
10.0
7.5 5.0 2.5
S.Tlm
0.0
-2.5 -5.0 -7.5
30.0 27.5
25.O 22.5 2O.O 17.5
M~ i lf4.72 aA4 1LII 10J8 8J~0 7,O4O 4A~ 0~4Q
NODAL DISPLACEMENTS mm i
14.'nz 1=t,44
Fig. 3. Balance of forces acting on fluid stage 2.
Difficulties associated with the numerical solution of hydraulic fracture models are related to the unique physical behavior of the system. For instance: (a) change of fluid pressure, within the crack, at any point affects the displacement field (crack opening) everywhere (global data structure); (b) the integral constraint, K]c, is satisfied continuously in time, so pressure must change throughout the crack when a local perturbation is introduced; (c) one free external moving boundary, the perimeter of the fracture, and one partially constrained internal moving boundary, the perimeter of the fluid front, are simultaneously controlled by mass conservation, fluid(s) rheology, the critical stress 30,0 27,5 25,0 22,5 20,0 17,5 15,0 12.5 10.0 7.5 5.0
-5.0 -7.5
30.0 27.5 25.0 22.5 20.0
1Lie
17.5
15.0
12.5
10,,11o Moo
10.0
UlO
15.0 7.5
5.0 2.5
0.0 -2.5
7,,o4o s.z,m 4,4oo Mm 1,,am
-,5.O -7.5
Fw
NOOAL DISPLACEMENTS
II,ta i 1lM,b744 '~ ILlql lM OIJlOe0 IOAI ~ 7,1,O700 4O 0,1JI1m40O
-7.5 -5.0 -2.5 NODAL
0.0
2.5
DISPLACEM
5.0
7.5 10.0
ENTS M~
i 14,7'1 111.111 10,,1
12.5 10.0
7.5 5.0 2.5
0.0 -2.5 -5.0
7J~O
O.11110
-7.5
l°:q~.O ;.S ~.0:2'.S o'.0 ~,~~'.0*¢.5 i~.o
l°&o.o-r'.S ~.o ~d.S i;.O ~,~ 5'.0 ¢.5 *.io.o
Fig. 4. Sequence of prolated blob distortions modeled by PARFES1.
382 30.0
27.5
AMAURYF. JUNIOR and M. F. AHMAD
NODALDISPLACEMENTS Ja~
25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
i ~,1=
SOJN ~Jm
30.0
NODALDISPLACEMENTS
25.0 22.5 20.0
15.0
18-72
10.0 7.5 5.0 2.5
UAG i
17.5 15.0 12.5
33.12 97.se 24Aa =lJO
10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
17.5 12.5
1.44o
30.0 27.5
25.0
12.N lOae 7.~o
1,44o
0.0
2.5
5.0
7.5
ls.ee scum 7.1oo
5.0 2.5 0.0 -2.5 -5.0 -7.5
10.0
NODAL DISPLACEMENTS M~
i le.12
22.5 20.0 17.5 15.0 12.5 I 0.0 7.5
--10L(~0.0"--7'.5 ~ ' . 0 "2'.5"" OiO'" 2~5*" 510"" i15 "11).0
18JI4
-7.5 -5.0 -2.5
reJ4
1A40
:11,4o 11L72
0.0 -2.5 -5.0 -7.5
se.r~
4,.rod
MAG i 11.12 U IrtJle 14A41
20.0
zsJo
le,,84 12.11e lo,oe 7JaOO 4Jm
NODAL DISPLACEMENTS
22.5
e,t.4e
-10:~0.0"-7.5"'-5.0''-2'.5"'hi()''215 510'"7'.5""i1).0 27.5
30.0 27.5 25.0
l/Jle
!1.110 llL~ 111,84 12Jle lOae 1.t00 4.110 1.4,,0 SI0* rlg"id.0
1°:~.ii-7'.S ~.0 :L5 010 ~
Fig. 5. Sequence of prolated blob distortions modeled by PARFES1 (larger deformations). intensity factor and elastic behavior of the porous matrix; (d) chemical/kinetic characterization of the fluids, gels and proppants injected into the fracture is quite complex. In addition, fluid losses, to the involving stratum(a), can alter the relative composition of an injected mixture; (e) the convective motion of fluid stages within the fracture can be affected by two instability mechanisms: fingering (in-plane--x and y directions) and encapsulation (transverse--z direction). Although fingering can modify the morphology of the fluid stages' interfaces, encapsulation instabilities are most important in the more common practical circumstances of increasing fluid viscosity in successive stages. Encapsulation instabilities can substantially affect the average velocity of the heavier fluid(s) carrying the proppant. For instance, if most of the proppant is carried directly to the bottom of the fracture, post-fracturing oil/gas production can be significantly reduced.
DESCRIPTIONOFPARFES The algorithms forming PARFES are all based on a Euler-Lagrangian approach, e.g. nodes delineating the injection region are constrained while other nodes are free to move, according to the coupling of flow field variables with the elastic stress field.
The tracking of the interface motion and distortion between the multiple fluid stages is performed by PARFES1. A finite element algorithm with a local data structure is employed, and an elliptical conjugate-gradient-type iterative solver utilized. Matrices are sparse and not assembled. The mesh movement/ stretching assumes a simple elastostatic formulation: T,jj = f ,
in f2
u~=a~ inOQu
Tiinj= ~
in
(1)
c~ t
(2)
Tij = 21 ,~oul.i + 22(uia + uj,i) 30.O 27.5 25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
NODAL DISPrEMENTS
i 48,11.1 4802 40/..1 ~4JS IIJ 178.7 ImsJI lm.II IlIOJD M.:I~ 21,44
I°%o ~;S :d., ;.S of(; ~IS ;'.O ilS''I'd.0 Fig. 6. Necking of a prolated blob modeled by PARFESI.
Parallel finite element solvers for flow-induced fracture
30.0 27.51 25,0: 22.5 20,0 17,5 15,0 12.5 10,0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
NODAL DISPLACEMENTS MAG i
0.t47 0.134 0.1~1 0.1M 0.0N 0,0~ O,0"/O 0087 O.044 0,0111 0.00e
..~ . . . .
~. . . .
i ....
i ....
i ....
i.,
i
30.0 27.5 25,0 22.5 20,0 17,5 15.0 12,5 10.0 7.5 5.o
383
NODAL DISPLACEMENTS MAO
.
,q
2.5 0.0 -2.5
0.110 0.100 O.OOe O.OU O,077 O,O07 0067 O.O44 0,0~ 0.010 OJD06
-5.0
•
-?.5
i
"~%0.0"-7~ -5.0 -2.5 0.0 2.5 5.0 7.5 to.o -
30.0 27.5 25.0 22.5 2O.O 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
NODAL
-
DISPLACEMENTS
UA~ i o.me o.ora 0.071 oaea o.m4 o~
•
0.~1 o.m~ o.oo,e.
A
30.0 27.5 25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.O -7.5
NODAL DISPLACEMENTS
iM~
A
0.071 0,Q84 0.O87 0.081 OO44 O~ O~ 0,,~ 0017 0,010 O~
-1°:%.o :;.S-6'.0 :~.5 o'.o 2'.5 5'.o ¢.5 lo.o Fig. 7. Sequence of oblated b/ob distortions modeled by P A R F E S I .
30.0 27.5 25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0 -2.5 -5.0 -7.5
N O D A L D I S P L A C E M ENTS M~ i
o.oa o.m.t o.ms o,.ow 0.04t ¢olo o.m4 O.Olo o,o14 ¢.me
o.oos
30.0 27.5 25.0 22.5 20.0 I 7.5 15.0 12.5 10.0 7.5 5.0 2.5 0.0
NOOALDISPLACEMENTS
-2.5 -5.0 -7.5
-1°:9~.0 :;.S :d.o :2'.5 0'.o ¢.5 5'.0 ;.5 ic;.o 30.0 27.5 25.0 22.5 20.0 17.5 15.0 12.5 10.0 7.5 5.O 2.5 0.0 -2.5 -5.0 -7.5
N O D A L D I S P L A C E M ENTS
iMAQ
o.olo o,ot2
Fig. 8. Sequence of oblated blob distortions modeled by P A R F E S I (cont'd),
M~ 0.041 0041
O,,NI 0.01O OJ~l OAOI O,,W~
384
AMAURY F. JUNIOR and M. F. AHMAD
:!:! ::i:i:i:i:i:i:!: :.:.:.:.:.:.:.:
!i???i???~??~?i :i:!:i:i:i:i:i: !:!:!:i:i:i:~:i :i:i:i:i:i:i:i: :.:->:,:.:.:.1
:;:: ... ::::
:::::::::::::: Z,:.:4.;.:.
!!!i
!!!!ii!iiii!i!
::::::::::::::
:?:?:?:?:?:?:):
:::~:~
i:i:i:!:~:~:~:~ .:.:.:.:.:.:.:.
iiii!i!iii~
:i:i:i:i:!:~:i: ,.v,.,,,..,,. :2:::::::::::::
[ Onset of instability I
lack
l
..................
n=g ] ,,,..,,,
1.1.1"2 ".,,,.,
>>>i >>:.>
. . '. . ' .." . .
Fig. 9. Onset of encapsulation instability--lighter fluid migrating to the walls.
= Jab ~ wh + f 2 w ~ ) d ~
(3)
c~fZu n dot = 0
+ r (t]w h + t2w~) d(0D) d~f i b
ulj)]wi,~ h
(4)
where the displacements (u) at the boundaries are + [,t2(u ~,2 + u h2,, )]w h1,2+ [;t, (u h~,~+ u ~,2 ) +
coupled to the flow field and the elastic stress field-notice that the elastostatic (local) formulation above is not related to the more complex elastic stress field
22(u~,2 + u2,2)]w2,2 h h + [22(uh, l + uL2)w2,, h h } df~
Frosted--Wedge V e l o c i t y vs. D e n s i t y 10
•
'
'
'
'
I
.
.
.
.
I
.
.
.
.
I
.
.
.
Ratio
.
I
'
'
,"
'
Velocityattached019 Velocityunanachcd 019 Velocityattached020 Velocityunattached020 Velocityattached 021 Velocityunattached021
O 0
.B
E 0,1
0.01
0.001
, 1.4
i
i
i
I 1.45
i
i
T
i
I 1.5
i
i
i
i
[ 1.55
,
i
i
i
I 1,6
i
~
i
i
I 1.65
i
i
i
i 1.7
Rho2/Rhol Fig. I0. Experimental results for the downwards velocity of a blob of higher density fluid (2) displacing a lighter fluid (I). Vertical Hele-Shaw type cell configuration, one float glass wall, the other a "frosted" float glass wall. Three different fluids used as the heavier fluid.
Parallel finite element solvers for flow-induced fracture
385 err_P'W
y
EIToT 19g+1 ,
~
en'_VF
IE+O
err_P,J"D~
IE-I 1rE..2
err_RT
195.3'
tn'_VOL
1"£-4
IZ-y'
\1
1~.6"
"'--.....z
~f
.........
1~-7'
IE.& IE.9-
~x,
IE.I0" 1E.11"
~.~ ,,.~,~ ~ ~.
\~
'V'
V~l e, ! ,4 v
"
IE-12. 1E.13. IE-14' 1E-15.
Fig. 11. Downwards movement of three immiscible fluid stages contained within a channel of variable height 6(x,y, t).
(global) formulation, except for the coupling at the boundaries. The finite element basis functions w h and their first derivatives are assumed to be square integrable in fl, with support zero at the prescribed displacement boundaries (~3~.). The numerical implementation of PARFESI allows the mesh to wrap around itself (the topology may change from convex to nonconvex), in order to track the development of more complex shapes (for instance, tracking the development of fingering patterns). A modified nonlinear Newton-Raphson algorithm, PARFES2, incorporates the physics governing the multiple fluids flow motion. The nodal movement and the fluid flow velocity fields are coupled to the pressure and crack opening distribution fields.
¢~_q,IW Error
1ff*.16' O~g+O
1E+2
2E+2
3"£÷2
4'F.+2
5E+2
6E÷2
Tune (seconds)
Fig. 13. PARFESAX results for 1000 mesh points and 101 time steps. Numerical error (comparison against analytical results) for: pressure at the wellbore; velocities at the wellbore; velocities at the fluid front; mesh velocities at the fluid front; fluid front position; and total volume conservation. Although the algorithm's data structure follows the same ideas implemented in PARFES1, its results depend on the 3-D elastic field and fracture criterion coupling performed by R 3 D H - - a general fracture simulator developed in a separate project) '8,9 R3DH is used to synchronize the global variations in crack opening distribution with the predicted pressure flow field given by PARFES2. The coupling of the two algorithms would represent a serious computational bottleneck (due to the physics of the problem), if calculations were performed for both algorithms in refined meshes. However, a better use of the computational resources may be obtained by immersing the very refined meshes needed in PARFES2 with somewhat coarse discretizations done in R3DH (the global data structure surface integral scheme employed in
IE+I IE+O
,~2
i
k\
led
....... e~._._P,F err VOL
II II
~\
Error IE+O i
IE-4
_,I
lq£-~
IE.3 \ l IE-4
1~-sl I \ IE.IO ~
!!1
i'
1~-; t
IE.11
IE-8
1~-91
'~-'~,~., v
,
V
I
~
1E.101
v-
.~rrRP/XTt J
1E-11 1
1E-I~;
1E-12 1
1E.1C 1000
2000
3000
4000
5000
T~e
6000
~econds)
Fig. 12. PARFESAX results for 100 mesh points and 101 time steps. Numerical error (comparison against analytical results) for: pressure at the wellbore; velocities at the wellbore; velocities at the fluid front; mesh velocities at the fluid front; fluid front position; and total volume conservation.
' '
1E.13 1
IP,.14 1 I~.15 1 IE.I6 0
100
2flO
300
400
~
.zrr_R~DOT_tl01
I
5-00
600
Tune (seconds) Fig. 14. PARFESAX results for 13 and 101 time steps (1000 mesh points). Numerical error (comparison against analytical results) for mesh velocities at the fluid front.
386
AMAURY F. JUNIOR and M. F. AHMAD - -
e r r ' P W tO06
Wellbore Pressure Error
,
,trot - -
1E+O
@-- Pt~ssure_2-F1uids_T-200
err_'1~V_t600
1.0E+0
!
err_'pW_tI5000
Pressure_5-Fluids_T-200
--~
err_q'gV_t6000
........
Pres~ure_l-Fluid_T-200
I --I~
err_PW_t6
•
\ 1.0E-3
. . . . .
1.0E4
. . . .
\
1.0E-5 1.0E~5
.
.
.
.
1.OE-7
\
1.0E-8 IE-5 1E.4
1E.3
rE-2
IE.1
IE+O
1E+1
1E+2
IE+3
IE*4 ffime
F i g . 15. P A R F E S A X results, 101 t i m e s t e p s , for d i f f e r e n t ( 0 . 0 6 , 6, 6 0 0 , 6 0 0 0 a n d 1 5 , 0 0 0 son against analytical results)
1.0E-9
IE+5 (seconds)
u s i n g 1000 m e s h p o i n t s a n d values of total running time s). N u m e r i c a l e r r o r ( c o m p a r i for p r e s s u r e at t h e w e l l b o r e .
R 3 D H , although computationally intensive, does not require refined meshes, except near the crack tip perimeter). P A R F E S A X contains all the intrinsic physical modeling capabilities o f P A R F E S 2 combined with the mesh motion o f P A R F E S 1 , but it is limited to track the movement o f multiple fluid stages for axisymmetrical cases. In spite o f this simplifying symmetry, the accurate solution o f cases o f interest is also computationally intensive. By comparing the results o f this code with derived analytical solutions o f one fluid stage problem, a great amount o f physical insight has been gained (see next section).
1.0E-10
. . . . . . . .
1,0E+O
i
1.0E+I
. . . . . . . .
i
1.0E+2
. . . . . . . .
I
. . . . . . . .
1.0E+3
I
1.0E+4
Total Number of Mesh Points Fig. 17. PARFESAX error analysis for the final predicted wellbore pressure versus the total number of mesh points. The same boundary conditions and fluid characteristics of the one-fluid stage are imposed on both the two-stage and five-stage problems. Result for 200 time steps (notice that the error is weakly dependent on the number of time steps---compare with Fig. 16). The mass and m o m e n t u m conservation formulation (assuming an aspect ratio o f order A / L ~ A / H ~ 10-s), for both P A R F E S 2 and P A R F E S A X , can be expressed as
I
J ~(t)
[(ip'v a),, + ('piv ,~%),k + 2pv ' ~ ' V,o~] '~j d f~ = 0 V '~kj e L2['n(t)]
(5)
Wellbore Pressure Error
Wellbore Pressure Error
l-
Pressure_l -Fluid_R-50
Pressure_l-Fluid T-50
Pressure_2-Fluids R-50
Pressure 2-Fluids_T-50 t .OE+O
~lk
Pressure_5-FluidsT-50
1.0E+O
1.0E-1
1.0E-I
1.0E-2
1.0E-2
1,0E-3
1.0E-3
l .OE4
1.0E4
Prossure_5-Fluids R-50
1.0E-5
1.0E-5
1,0E~6
1.0E-6
1.0E-7
I.OE-7 1.0E-8
l.OE-8
1,0E-9
1.0E-9 1.0E-10
1.0E-10 1.0E+0
1.0E+I
1.0E+2
1.0E+3
1.0E+4
Total Number of Mesh Points Fig. 16. PARFESAX error analysis for the final predicted wellbore pressure versus the total number of mesh points. The same boundary conditions and fluid characteristics of the one-fluid stage are imposed on both the two-stage and five-stage problems. Result for 50 time steps.
1.0E+O
1.0E+l
1.0E+2
1.0E+3
1.0E+4
Total Number of Time Steps Fig. 18. PARFESAX error analysis for the final predicted wellbore pressure versus the total number of time steps. The same boundary conditions and fluid characteristics of the one-fluid stage are imposed on both the two-stage and five-stage problems. Result for 50 mesh points.
Parallel finite element solvers for flow-induced fracture
£n,,~[% - ~ 62 [- (',,p),,+'p'vg,]]%d'fl= o V 'gpjE H~['n(t)]
(6)
where for each fluid stage (prescript i): ~v is the volume fraction; *p the density; ~/~the viscosity; % the average (in-plane) velocity component in the k direction; ~V~o~the (transverse) fluid loss velocity (through both fracture walls); g~ the gravity acceleration component in the k direction. The crack opening field 8 is coupled to the fluid(s) pressure field p by means of a surface integral scheme, s encompassing the elastic stress field around the fracture walls, and the stress intensity factor in mode I. The finite element basis functions q~ and ~ are assumed to be square integrable in ~, for any instant of time t. The first derivative of q~is also assumed to be square integrable with support zero at the boundaries. Assuming the following notation for integration in [1 and its boundaries Of~ (atlb~) - [o
f~
(up to fifth order) Backward Differentiation schemes are used. For instance, expanding Eq. (8) and applying Newton-Raphson for time steps (superscript n) n ~<2, with 6 = d + A 6 and cr = p + Aa, gives the following incremental equation: i-
i-
A
~
- 6 At 4 [ p d2e
L'~
n
n
aalnk]
axk
J
+ 2 At ~s \7~ a ~-~-XkI ~Xk/ i-
~./
\ I~ oxk
n
ICXkl
i -2
In
--2Atli~PXkA6lnkl ~
a~bmd[)
[allnm] =--[" atnmd(df~),
387
(7)
-6At~s~(~dEgkA6 -
07~\ ~ dxk/
the previous Eqs (5) and (6) can be combined as
~(ff
I'/')-~s~=O
L'*'
~--n~
&*
+ ¢tact~ 63Tg, lnk
+ 2 At('~)?k A 31 dxk/
J =2Ato~I~d3"~nk
- [~t~ aq'£,ln,]
1
i - 2 d3"gklnk In - 2 At "B(~2['~fi
2/i~ 3 da
o g~ ~Xk/) + 2('tS'V,o,,[~) = 0
(8)
where
+ 2 At %~2s
- d3gk Ox,/
- 2 AtQ~df(k ~xg)"+ 2 At['~dtltf(~lnk]" AR
i-
~XS = ~ RR
~n =
pRgLR - -
c9 a~\"
-- 2 At ~s2([~- d3 ~
~xk/
(9)
O"R
are the nondimensionalization variables; a the excess pressure field (a =p-ac); cr~ the confining stress field; -£k the nodal velocity component in the k direction; and ~ the finite element weight function. The terms in "{'karise from passing the time differentiation out of the integral (first term on the L.H.S.) in Eq. (5). For the implementation of PARFES2, two implicit A-stable m linear multistep schemes are used to time integrate Eq. (8). An Adams-Moulton of second order is used for the first time step and a Backward Differentiation, also of second order, for the subsequent time steps. For the axisymmetrical problem (the gravity term is dropped in PARFESAX), after the second time step, implicit adaptive variable order
- 4 At('H'V~ossl~) n +4(,~61~u).-i
3(~Hdlq')"
(,~61q,). 2.
(10)
More details of the above derivations and equations can be found elsewhere. 2 RESULTS AND CONCLUSIONS There are two versions of PARFESI: one for coarse grained and the other for fine grained parallel computing environments. The coarse grained version was tuned to work at close-to-peak performance at the NCSA's CRAY Y-MP. For instance, substantial reduction of the number of instructions per second and a ten-fold decrease on wall-clock time was obtained with respect to the original version (various
388
AMAURYF. JUNIORand M. F. AHMAD
performance issues, such as efficient array indexation and smaller loop unrolling were already considered on the original version). The fine grained version was written in CM-Fortran and tuned for use on the NCSA's CM-2 (Connection Machine 2). The code has been almost fully rewritten to take advantage of the Massively Parallel environment. Experimental analysis of the convection motion (due to gravity-driven forces) of a fluid stage blob moving downwards within a narrow cell (vertical Hele-Shaw type cell, filled with a different fluid), has demonstrated 2'3 the need to characterize in-plane deformations of the blob. An immiscible heavier blob displacing a bath of lighter fluid, sketched in Fig. 2, will sink against the following counteracting forces (see Fig. 3): (i) shear stresses at the walls (no encapsulation) or at the interfaces with the lighter fluid in the transverse direction (if encapsulated); (ii) buoyancy; (iii) shear stresses to displace the lighter fluid; (iv) surface/interfacial tension and surface wetting characteristics of both the displaced and the displacing fluids. PARFES1 is used to characterize the distortion of the blob under different in-plane boundary conditions. For instance, Figs 4 ~ present a sequence of increasing rates of prolation for a distorted blob--the vertical axis representing the direction of motion is reversed (the more elongated part of the blob is moving downwards). Figures 7 and 8 show one sequence of increasing rates of oblation - - t h e motion is towards the horizontal positive direction. The tendency of higher viscosity fluids to migrate to regions of lower shear rates may dramatically alter the wetting conditions at the fracture walls Wellbore Pressure Error
Pressure_l -Fluid_R- 1000 0-
Pressure 2-Fluids_R-1000
~
Pressure_5-Fluids_R-10CO
Wellbore Pressure (x 2x 3x 4x 5x) [ [
2.0E+6
t
1.8E+6 1.6E+6
--m~
Pressure (Variable Flux) S-Fluids
--@~
Pressure (Variable Density) 5-Fluids Pressure (Variable Viscosity) 5-Fluids
1.4E+6 1.2E+6 1.0E+6 8.0E+5 6.0E+5 4.0E+5 2.0E÷5
!
O.0E+O
........ 1.0E+O
........ i 1.0E+I
• •
) ........ 1,0E+2
i
i ........ 1.0E+3
i 1.0E+4
Total Number of Mesh Points
Fig. 20. Convergence analysis for the final predicted wellbore pressure versus the total number of mesh points tracking five fluid stages. The injected flux/density/viscosity are increased as x, 2x, 3x, 4x, and 5x for each stage (x is the flux/density/viscosity for the first fluid stage).
(see sketch in Fig. 9). This (transverse) migration mechanism, previously termed encapsulation, may substantially affect the average velocity of the heavier proppant-laden fluid(s) and, consequently, the placement of proppant particles within the fracture cavity. Although fingering (in-plane instability-modeled by PARFESI), can modify the morphology of the fluid stages' interfaces, encapsulation instabilities are more important in hydraulic fracturing. For Injected Flux (x 2x 3x 4x 5x) Flux (Variabel Flux) 5-Fluids
1.0E+O
2.0E+O --0-
1.OE-1
.................................................................
Flux (Variabel Density) 5-Fluids
1 .SE+O Flux (Variabe[ Viscosity) 5-Fluids
1.0E-2
10E-3
*-i
*
1.0E-4
1.6E+0
................
1.4E+O
:
"
...............................
~ ...................................
'
] .2E+O i
1.0E-s
i
1.0E+O
:i 1.OE-6
!
i
'
8,0E-1
1.0E-7
6.0E-I
1.0E4~
4.0E-1
1.0E-9
::
l O~.~o
........ 1.0E+O
i 1.0g+l
'
........
! I,OE+2
........
! 1.0E+3
2.0E-1 O.OE+O
. . . . . . . . .
1.0E+4
Total Number of Time Steps
Fig. 19. PARFESAX error analysis for the final predicted wellbore pressure versus the total number of time steps. The same boundary conditions and fluid characteristics of the one-fluid stage are imposed on both the two-stage and five-stage problems. Result for 1000 mesh points.
........ LOE+O
i ........ 1.0E+1
t ........ 1.0E+2
I ........ 1.0E+3
i 1.0E+4
Total Number of Mesh Points Fig. 21. Convergence analysis for the final injected flux versus the total number of mesh points tracking five fluid stages. The injected flux/density/viscosity are increased as x, 2x, 3x, 4x and 5x for each stage (x is the flux/density/viscosity for the first fluid stage).
Parallel finite element solvers for flow-induced fracture
389
fluid when displaced upwards is also constrained to move within the fracture boundaries. Another set of semi-analytical algorithms (also described elsewhere 3) were developed to characterize the influence of the fracture boundaries on the blob motion. The tracking of one interface motion was verified against derived analytical solutions using PARFESAX. An error analysis for various physical variables of interest, is shown in Figs 12-15. Multiple interface tracking was verified in two ways: (i) imposing exactly the same boundary conditions and fluid characteristics of the one-fluid stage analytical solution to the various fluid stages--
instance, experiments conducted for different fluids 3 (see Fig. 10) have demonstrated a spread of about two orders of magnitude on the downwards blob velocity. This velocity spread led us to develop an analytical algorithm, described elsewhere) to predict the size of the encapsulated layers, as well as an experimental criterion to predict the onset of encapsulation. The downwards motion of multiple, immiscible, fluid stages (see Fig. 11) is far more complex than the two-stage problem, especially when modeling inplane distortion. Excluding fluid losses through the fracture walls, the various stages compete against each other when filling the fracture, i.e. the lighter
6 4 Speed Up 2 0
Paral lel Regions Iml = Serial t,lork
[]
: Par.allel t,lork
sS
ss SS sS s
8
s
s~
s~ 4
7
f
s sS
-t~ t ¢s
6
/
is
--
J /J
Speed 5
s~
--
Up
s t
4
--
s
i
ev e s
J
/ f sS
3
-
,/ •
2
1
-
2.5
f
o o
..."
.
¢-~(-
2.0
~.5. ~
2,3 . , . . , - .r
2.4 . .
.
,,. =, ~ .
•
,IF ~ ' p ~ "
2.3
.
....
/I.0 I
, I
1
2
.
I
I
I
I
I
I
3
4
5
6
7
8
# CPU$ Fig. 22. PARFESI optimization resultsfor a 100 x I00 mesh, usingthe NCSA's CRAY Y-MP. Compare parallel scale-up when the number of mesh points is increased (see Fig. 23).
390
AMAURYF. JUNIOR and M. F. AHMAD
Figs 16--19 show an error analysis for 1, 2 and 5 fluid stages; (ii) doing a convergence analysis in the general case (no analytical solution available), see Figs 20 and 21. With respect to the suitability of the algorithm(s) to parallelization, optimization results for PARFES 1, using the NCSA's CRAY Y-MP, show a good parallel scale-up when the number of mesh points is increased (see Figs 22-24). The fine grained implementations, using the Connection Machine 2 (at various locations), confirmed our expectations, e.g. we were able to run very large problems which almost 98% utilization of the processes. Two points--one related to the physics of the multiple stage flow problem, the other related to the
computational data structure--are worth mentioning: (i) in order to capture the development of relevant physical phenomena, such as fingering and encapsulation instabilities,3,tH3 of the multiple fluids motion, each full iteration loop of PARFES2 requires very refined meshes (provided by PARFES1), and it has to be done often and computationally cheaply (also, the algorithm should allow the mesh to wrap around itself); and (ii) the localized data structure of the PARFES1 algorithm (matrices are not assembled), the use of iterative solvers, and the unique layout of local matrices (despite of an overhead in memory usage) have proved to be very successful, and those ideas will be extended to some modules of PARFES2.
m
Speed Up
Parallel Region: m : Serial Work [ ] : Parallel Work
O0f
~ummar~of residl:resid1~49
tS jJ S~ O
8
S°
SS S
SS
~S jS
7 --
Ip
s~
sI sS / sS / sS
Speed
5
-
,/
Up
4o6
,,/"
4.3 ~ 3.9 ..~,~1V"~4.2 '
/// 4 --
///
,,/" 3
-
/
3,5 ~ , - ~ ¢ ~
4.6
3o9
Z.O~ "¢'- 3'5 3.O J
,/,~,1,. 1.8
2 --
1.07"
I --/0.o I"O0.o 1
1
0.1
0.0
0.0
0.I
0.1
0.0
. . . . _'A . . . . . . . . . . . ~ . . . . . . . . . . . ~ ......... "1 ........... L... ....... ~.......... ~i
2
:3
4
o CF~Js
5
6
7
8
Fig. 23. PARFESI optimization results for a 150 x 150 mesh, using the NCSA's CRAY Y-MP.
Parallel finite element solvers for flow-induced fracture
391
]
SpeedUp Parallel Re9ior~ m
= Serial Work
restdl: _
[]
[] = Parallel IJork
r e s l d l Q 4 9 - 2000
8.0
[]
keyshell
7.9
Ke~
--
Speed
7.0
andahls
---
overhead dedicated maximu~ variable ~Ini~u~ variable linear speed up pro9ram speed up
.....
6.9 6.0
-
-
5.9
...... ......
5.0 4.9
-
Up
4.~4.0 3.% 2.7 2oOy 1,4" 2.0
2"
1.0 0.3 0.1
1
2
0.1
3
4 5 G 7 8 w CPU$ Fig. 24. PARFESI optimization results for a 150 x 150 mesh, using the NCSA's CRAY Y-MP. Result for the most computationally intensive subroutine. The program presents a good parallel scale-up (even in its coarse grained implementation) when the number of mesh points is increased.
Acknowledgements--The authors would like to thank: the Gas Research Institute for supporting this research project at MIT; the National Center for Supercomputing Applications for extensive use of their computers; the Advanced Computing Research Facility at the Argonne National Laboratory for use of their computers; and the Thinking Machines Corporation for use of the CMNS. The following MIT undergraduate students, under our supervision, have contributed to the work described here: Yitwah Cheung, Adeeb Shana'a and Lisa Zhang. Thanks for their many valuable suggestions and help.
REFERENCES 1. M. P. Cleary, "The engineering of hydraulic fractures-State of the art and technology of the future," SPE Distinguished Author Series--Journal of Petroleum Technology 13-21 (January 1988).
2. M. P. Cleary and Amaury F. Junior, "Proppant convection and encapsulation in hydraulic fracturing: practical implications of computee and laboratory simulations,'* SPE Journal 24825 (October, 1992). 3. Amaury F. Junior, "Modelling of flow during fracture growth driven by multiple fluid stages using highperformance parallel computing algorithms," Ph.D. Thesis, MIT, 1992. 4. K. Y. Lam, "The development of a fully threedimensional simulator for analysis and design of hydraulic fracturing," Ph.D. Thesis, MIT, 1985. 5. D. Barr, "Leading-edge analysis for correct simulation of interface separation and hydraulic fracturing," Ph.D. Thesis, MIT, 1991. 6. Amaury F. Junior, "An adapted finite element method for massively parallel processors," The Fourth Conference on Hypercubes, Concurrent Computers, and Applications, Monterey, California, pp. 1329-1333, 1989.
392
AMAURY F. JUNIOR and M. F. AHMAD
7. K. K Mathur, "The finite element method on a data parallel architecture," Fifth International Symposium on Numerical Methods in Engineering (edited by R. Gruber, J. Periaux and R. P. Shaw), pp. 599~11, Springer-Verlag, 1989. 8. M. P. Cleary, M. Kavvadas and K. Y. Lam, "Development of a fully three-dimensional simulator for analysis and design of hydraulic fracturing," SPE/DOE Journal 11631 271-282 (1983). 9. M. P. Cleary, D. T. Barr and R. M. Willis, "Enhancement of real-time hydraulic fracturing models with full 3-D simulation," SPE Journal 17713 79-93 (1988). 10. J. Sand and O. Osterby, "Regions of absolute stability,"
Computer Science Department, Aarhus University, Denmark, 1979. 11. N. R. Anturkar, T. C. Papanastasiou and J. O. Wilkes, "Linear stability analysis of multilayer plane Poiseuille flow," Physics of Fluids A 2(4), 530-541 (1990). 12. D. D. Joseph, K. Nguyen and G. S. Beavers, "Nonuniqueness and stability of the configuration of flow of immiscible fluids with different viscosities," Journal of Fluid Mechanics 141, 319-345 (1984). 13. S. J. Weinstein, V. E. B. Dussan and L. H. Ungar, "A theoretical study of two-phase flow through a narrow gap with a moving contact line: viscous fingering in a Hele-Shaw cell," Journal of Fluid Mechanics 221, 53-76 (1990).