Simulation of Spontaneous Imbibition Using Rayleigh-Ritz ... - CiteSeerX

PETROLEUM SOCIETY

PAPER 2004-228

CANADIAN INSTITUTE OF MINING, METALLURGY & PETROLEUM

Simulation of Spontaneous Imbibition Using Rayleigh-Ritz Finite Element Method—A Discrete Fracture Approach S.P. KAUL, E. PUTRA, D.S. SCHECHTER Texas A&M University This paper is to be presented at the Petroleum Society’s 5th Canadian International Petroleum Conference (55th Annual Technical Meeting), Calgary, Alberta, Canada, June 8 – 10, 2004. Discussion of this paper is invited and may be presented at the meeting if filed in writing with the technical program chairman prior to the conclusion of the meeting. This paper and any discussion filed will be considered for publication in Petroleum Society journals. Publication rights are reserved. This is a pre-print and subject to correction.

size example, using discrete fracture approach, was developed and its results compared with a commercial simulator.

Abstract Spontaneous imbibition plays a very important role in the displacement mechanism of non-wetting fluid in naturally fractured reservoirs. We developed a new 2D two-phase finite element numerical model, as available commercial simulators cannot be used to model small-scale experiments with different and complex boundary conditions. For the non-linear diffusion saturation equation we cannot apply Rayleigh-Ritz Finite Element Method (FEM). Traditionally, the way around it is to use Galerkin FEM or Mixed FEM formulation, iterative nature of those, makes them unsuitable for solving large-scale field problems. But if we truncate the non-linear terms, decouple and solve analytically the dependent variables from saturation – the primary variable, this non-linear FEM problem reduces to a simple weighted integral weak form, which can be solved with Rayleigh-Ritz method. The advantage of this method is that it is non-iterative, which reduces computation time. We compared our numerical models with the analytical solution of this diffusion equation. We validated a Finite Difference Method (FDM) numerical model using X-Ray Tomography (CT) experimental data, and then went ahead and compared the results of FEM model to that of FDM model. A two-phase field

Introduction The quest to produce more oil has led various researchers to evaluate more complex reservoirs such as naturally fractured ones. The complexity of fluid flow in these types of reservoirs arises from the fact that there are two media, which can allow fluids to flow through them. This leads to the logical conclusion that there are two principal media of fluid flow. The first is through the porous medium called the matrix flow and the other is the flow through the fracture network. At the heart of this phenomenon lies the problem of interaction of these two flow media with each other. In naturally fractured reservoirs the matrix acts as a source where hydrocarbons are present whereas the fractures facilitate in fast recovery of these hydrocarbons. Hence it is important to study what makes the matrix produce more oil. Water is used as a means to efficiently displace oil. But fluid flow in porous media, which is determined primarily by capillary force, is relatively difficult phenomenon to quantify, as per Craig1, and there has been much research effort directed in this direction like that of Handy2, Garg et al.3,

1

Babadagli and Ershaghi4, Li and Horne5, Akin and Kovscek6, Reis and Cil7, Zhou et al.8, to name a few. Also, depending on the geometry of the fracture, there may or may not be any capillary force in the fracture network. This force is responsible for imparting the spontaneity to fluid flow within naturally fractured reservoirs. Given the complexity of quantifying the spread of fractures, it is even difficult to ascertain the limits where the fracture flow acts as an independent flow entity instead of being a part of porous matrix. Although very wide in its scope, fluid flow studies in

The major difference between finite difference and finite element (Galerkin Method) arises from the fact that in the case of former, all the coefficients of the differential equation are linear and direct solvers can be readily applied whereas in the latter, we have to use iterative solvers to ascertain the values of these coefficients before actual inversion can take place. Thus we reach a conclusion that as long we use Galerkin type of

fractured media, as we have narrowed it down, deals

formulation, finite element would be time consuming and

exclusively with the study of this spontaneous phenomenon that

computer memory taxing.

helps displace oil out of matrix. This paper deals with

In comparison to the Galerkin FEM, the Rayleigh-Ritz FEM

quantifying, with the help of a numerical model, the

derives its simplicity from the fact that its formulation, known

experimental study of spontaneous imbibition process in both

as the weighted-integral variational weak form, incorporates the

unfractured and fractured cores having negligible capillary

boundary condition along with a form of differential equation.

pressure in the fracture.

This makes it easy to solve that form of differential equation,

The challenge associated with undertaking such

without any iteration. Also, in general, higher-order

experimental studies, on a laboratory scale, is that commercially

interpolation functions are to be used in Galerkin FEM since the

available simulation softwares have some limitations. They are

weighted-residual form, Galerkin FEM, does not include any

not designed for such a task. Out of the two popular softwares

specified boundary conditions. Rayleigh-Ritz FEM results in

available to us, one has to have aquifer as a constant pressure

mathematically simpler formulations. An attempt has been

boundary condition, whereas the other would not run without

made to use this Rayleigh-Ritz FEM, in which iterative nature

wells in place, thus undermining the spontaneity of the process

of the solution methodology is completely overcome, in

of fluid flow. This made it amply clear to us that, if we had to

comparison to Galerkin FEM. The resulting matrix, which is

understand the spontaneous imbibition process on the

solved, is also of lesser band-width as only two equation are

laboratory scale then, we would have to come up with our own

solved instead of four. The direct advantage of this is that we

numerical model which would have to be flexible enough so

can incorporate all the geological and geophysical details into

that the affect of various parameters on spontaneous imbibition

an efficient mathematical model which is less mathematically

process can be studied in-depth. As a tool to verify this

intense than existing finite element methods and can simulate

numerical model, we used X-Ray Tomography, as it is the best

fluid flow in porous media. This is achieved without having to

tool available to study fluid flow in porous media.

resort to large-scale refinements, to capture reservoir

Solution to spontaneous imbibition differential equation is

heterogeneities and details, as done by finite difference method.

usually derived with the help of FDM scheme. In case of naturally fractured reservoirs, ever since the model proposed by

Basic Concepts

Warren and Root9, various researchers such Mattax and Kyte10,

There are two basic common concepts that are used in this

deSwaan11, Kazemi et al12, Civan13 to name a few, have

method.

underlined the importance of this process and help understand

The first concept is that multiplication of numbers is

recovery from such reservoirs. Attempts have also been made to

independent of the order of the numbers. This means that if two

solve fluid flow equations with the help of FEM by Javandel

numbers are multiplied, the answer is independent whether we

and Witherspoon14, Settari et al15, Young16, McMichael and

multiply first number with the second or second with the first.

Thomas17, etc. and specifically for naturally fractured reservoirs

The second concept is that if saturation moves from point A

by Asfari and Witherspoon18, Bhatia et al19, Naji and Kazemi20,

to point B (spatial domain) during a given time T, then

Sutopo et al21, Karimi-Fard and Firoozabadi22, etc. Invariably

saturation at each point, between A and B, will vary from 0 to 1

these researchers have approached the simulation problem from

(saturation domain). Any variation of saturation from 0 to 1 is

Galerkin finite element method point of view.

independent of variations of porosity and permeability. The 2

only way this is affected is, if there is another fluid in the pore

and late time solution where we have to consider the effect

spaces. In that case saturation of water will vary from connate water saturation,

S wi

of gravity:

S or .

to connate oil saturation,

⎡ ⎛ ρgR 2 ⎞ ……………………...….(3) h ⎤ ⎟⎟t h + H eq ln ⎢1 − ⎥ = −⎜⎜ ⎝ 8µ ⎠ ⎣⎢ H eq ⎦⎥

Concentrating on the former easier case, mathematically both these concepts can be represented by: T

B 1 = ⎡ ∫ (SpatialDomain )• ∫ (SaturationDomain )⎤ ⎥⎦ 0 ⎢⎣ A 0

The non-dimensional form of above equation is obtained by dividing throughout by Heq:

What is clear from above expression is that the solution of

⎡ ⎛ 1 h h ⎤ + ln ⎢1 − ⎥ = −⎜⎜ H eq ⎣⎢ H eq ⎦⎥ ⎝ H eq

saturation domain is known since we know the limits of the integration for all times, but solution of spatial domain is

⎞⎛ ρgR 2 ⎟⎜ ⎟⎜ 8µ ⎠⎝

⎞ …………………...(4) ⎟⎟t ⎠

unknown and has to be assigned. Thus we have: If we represent the characteristic time, τeq as: T

B 1 = ⎡ ∫ (Answerunknown )• ∫ (Answerknown )⎤ ⎢⎣ A ⎥⎦ 0 0

τ eq =

8 H eq µ

ρgR 2

………………………..………………..……(5)

leads to the solution:

Washburn Law and Spontaneous Imbibition Theory – Validation of Basic Concepts

⎛ t h ⎞⎟ ⎤ ⎡ ⎜− − ⎜ τ eq H eq ⎟ ⎠⎥ ⎝ ⎢ h = H eq 1 − e …………………...…………(6) ⎢ ⎥ ⎣⎢ ⎦⎥

a) Capillary Tube:

These two solutions clearly prove that the early time

In order to apply these basic concepts we have to understand

(diffusive part) solution varies as a square of time whereas the

fluid flow under the influence of capillary forces. This forms the

late time solution varies as exponential of time. This is the

basis of spontaneous imbibition in porous media. Most of what

Washburn’s Law24 and the characteristic time differentiates

follows here is taken from Dubey et al23. Taking up the simplest

between the two times.

case, we study capillary suction in glass tubes as shown in figure 1. If a capillary tube is immersed in water, the resulting

b) Porous Media – Single Phase Flow:

viscous flow pressure drop is the result of the upward capillary

The Washburn’s Law24 can be extended to single-phase flow

pressure and downward pressure due to the gravity (weight) of

in porous media with the assumption that it is a collection of

the fluid column. Rearranging it, for small change in volume,

bundle of tubes. Thus assuming permeability:

we have:

R2 k= 8

Water − Air cos θ ⎛ 8µ ⎞ dh 2σ = − ρgh …………………….(1) ⎜ 2 ⎟h dt R R ⎝ ⎠

…………………………………………………….(7)

we obtain similar solutions. The only difference is in the

This is the basic differential equation. There are two possible

expression of characteristic time, τeq, for porous media:

solutions to this equation, which represents early time (diffusive part) when the influence of gravity can be neglected as:

⎛ σ Water − Air R cos θ h 2 = ⎜⎜ 2µ ⎝

⎛ H eqφ ⎞⎛ µ ⎞ ………………………………………...(8) ⎟⎟⎜ ⎟ τ eq = ⎜⎜ ⎝ ρg ⎠⎝ k ⎠

⎞ ………………………………..(2) ⎟⎟t ⎠

3

Here we take interfacial tension instead of surface tension in

The validity of this law as it is applied to porous media flow 3

can be demonstrated by the data given by Garg et al’s X-Ray

calculating the capillary rise. The validity of this equation was

CT experiment. If we plot this data of the Wasburn’s equation

established with the help of data similar to Garg et al’s3 data

on a semilog paper, the diffusive part will fall on a straight line

and generated with the help of C M G

as shown in figure 2. When we plot the one-dimensional

®

. Figure 4 shows the

deviation of the diffusive data from the exponential data. This

solution of the same equation against the experimental CT data,

can be easily assessed by plotting the data on a semilog paper,

we find that they are in pretty good agreement as shown in

calculating the characteristic time, and then comparing it with

figure 3.

the value given by the formula. Both the values are in pretty good agreement. Our concept is applied here to derive the value

c) Porous Media – Multiphase Flow:

of

The known solution in case of one-dimensional multiphase flow, obtained with the help of Laplace transform, is given by:

⎛ h h = H eq erfc⎜⎜ ⎝2 t

µ Multiphase . Wherever possible this value is calculated

analytically or with the help of Mathematica

®

or Maple

®

,

which provides us with the solution. Thus we have proved what

⎞ ………………………………………...(9) ⎟⎟ ⎠

we started with, and that indeed answer to the problem is possible if we are able to ‘decouple’ the answer into its component parts.

24

There is no known Washburn type

solution and this

equation does not lead to one either. Applying the basic

Dual Porosity versus Discreet Fractures

concepts as discussed previously, the governing differential

Much effort has been directed towards having a rigorous

equation is given by:

mathematical model of fractured oil reservoirs. Principally, these are the two approaches, which presents a solution to fluid

⎞⎛ g ⎞ kk rw k ro dh ⎛ ⎟⎜ ⎟ ρ (H eq − h )− ∆ρh .….(10) h = ⎜⎜ dt ⎝ (µ w k ro + µ o k rw )⎟⎠⎜⎝ φ ⎟⎠

[

]

flow problem in fractured reservoirs. Researchers have, for past many years, defended one over the other. Both approaches have their own merits and demerits and we intend to enumerate a few

Assuming,

here, in order to bring out the complexity of the issue. In discreet fracture numerical model, as per Russell25, the length

1

µ Multiphase

1⎛ ⎞ ……………………....(11) k rw k ro ⎟⎟ = ∫ ⎜⎜ 0 (µ k ⎝ w ro + µ o k rw )⎠

scale of the fracture is very small as compared to the computational grid size and hence a “ continuum hypothesis” is plausible. Simply said, discreet fracture approach is a single porosity approach wherein realistic orientation of the fractures

the solution for multiphase flow is given by:

is taken into account. This is where FEM steps in. These fracture can be continuous or disconnected, which is the

⎡ (ρ + ∆ρ )2 t (ρ + ∆ρ ) h ⎤ ⎤ ⎡ ⎧ − ⎥⎫ ⎢− ρ ρ τ eq ρ H eq ⎦⎥ ⎪⎥ ⎪ ⎣⎢ h = H eq ⎢ 1 e − ⎬⎥ (12) ⎢ (ρ + ∆ρ )⎨ ⎪⎭⎥ ⎪⎩ ⎦ ⎣⎢

principal advantage of FEM. Another advantage is that we don’t have to calculate the transfer functions, as this is a single porosity approach. These fracture networks may or may not be coupled with the matrix. Discreet fracture approach is best

The characteristic time in this case is represented by:

suited where permeability of the matrix is very small and is known to give better results where the transfer mechanism

⎛ H eqφ ⎞⎛ µ Multiphase ⎟⎟⎜⎜ τ eq = ⎜⎜ k ⎝ ρg ⎠⎝

⎞ ……………………………….(13) ⎟⎟ ⎠

between the fracture and the matrix takes lot of time. Although one-dimensional, it can be seen from equations (12) and (13), the physical dimensions of the porous media also play an important part in the imbibition process. The realistic fracture orientations are not possible with FDM approach. Dual 4

porosity Warren and Root9 sugar-cube model approaches these

ρ o ∂S o = − ρ w ∂S w ……………………………….…..….... (15)

orientations as a series of approximations, since grid elements can only be square or rectangular. Kazemi’s model13,15 went

If we define capillary potential, Φc, as:

further and based the transfer mechanism of imbibition process on the amount of surface area of the matrix block available for

Φ c = Pc − ∆ρg sin α …………………………………….(16)

imbibition. The square and rectangular grids helps in assessing, with relative simplicity, the surface area/volume of these sugar

The end result of all mathematical manipulation leads to the

cubes. Thus spontaneous imbibition is a surface dependent

general displacement equation (one-dimesional), which is non-

phenomena. In this paper, we have attempted to use the best of

linear in nature and is given by:

both these approaches. We use rectangular grids and couple the matrix and the fracture. We follow Kazemi’s15 type transfer

⎛ ρ ⎞ ∂S ⎜⎜1 + w ⎟⎟ w ρ o ⎠ ∂t ⎝

functions between the matrix and fracture, for rectangular or square grid blocks, and leave it at that point. Transfer functions are not available for complex geometries. We believe that this idea can then be extended to triangular or any other polygonal

=

shaped grids and simulation performed just by adjusting the value of the transfer fuction.

⎤ ⎡ ⎛ ∂Φ c kk rw k ro ∂ ⎡ ⎢ ⎥ ⎢− ⎜ ∂S w ⎣φ (µ w k ro + µ o k rw )⎦ ⎢⎣ ⎜⎝ ∂S w

⎞⎛ ∂S w ⎞ ⎟⎟⎜ ⎟ ⎠⎝ ∂x ⎠

⎡ ⎤ ∂ ⎡ ⎛ ∂Φ c kk rw k ro +⎢ ⎥ ⎢− ⎜⎜ ⎣φ (µ w k ro + µ o k rw )⎦ ∂x ⎣ ⎝ ∂S w

Mathematical Model In order to arrive at a mathematical model, we have to start

2

⎤ ⎥ ⎥⎦

⎞ ∂S w ⎤ ⎟⎟ ⎥ ⎠ ∂x ⎦

from the basic equations. We take a control volume approach in

⎞ qws ……………………..…….…... (17) ⎛ k ro µ w ⎟⎟ + ⎜⎜ µ k µ k + o rw ⎠ φA ⎝ w ro

this case. Let us assume a control volume of dimensions ∆x, ∆ y, and ∆ z. Fluid flow occurring across this control volume obeys two basic laws of physics, namely:

Here we make the following assumption. For saturation, which varies from 0.0 to 1.0, the value of the first derivative is 1) Law of Conservation of Mass – also known as Continuity

small compared to unity and thus

Equation. 2) Law of Conservation of Momentum – represented by

∂S w