INTRODUCTION
According to Sarma and Aziz (2006), Naturally Fractured
Reservoir (NFR) is a complex system with irregular fractures network, vugs and
matrix blocks. They added that NFR can be defined as a reservoir having a connected
fractures network which has significant higher permeability than matrix. This
implies that production of hydrocarbons is highly dependent on the matrixfractures
interaction. This paper assessed the influence of different shape factor in
modeling the matrixfracture transfer rate.
NFR can be modeled by using double porosity concept. The double porosity concept
was introduced by Barenblatt et al. (1960), while
Warren and Root (1963) were the first to use double
porosity concept in reservoir simulation. Double porosity concept is having
two separate partial differential equations to define matrix and fractures flow.
Matrix usually have low permeability and high storativity, while fractures have
high permeability but low storativity. This suggest that matrix function as
a main source of hydrocarbons while, fractures become the flow path of hydrocarbon
production. For this reason, interaction between matrix and fractures should
be considered. This interaction can be described by using a transfer function.
The matrixfracture transfer function was given by Warren
and Root (1963) as:
Equation 1 showed that the matrixfracture transfer function
which requires shape factor to governs the flow.
The initial double porosity model of Warren and Root (1963)
used assumption of pseudosteady flow and the NFR system is simplified into blocks
of matrix and fractures set which looks like sugar cubes (Fig.
1).

Fig. 1(ab): 
Actual reservoir block and sugar cubes model 
Each cube is known as matrix that contained in within a systematic array of
identical and rectangular parallelepipeds. Matrix is assumed to be homogenous
and isotropic. All the fractures are continuous and may have different spacing
and width to simulate certain degree of anisotropy.
SHAPE FACTOR
Most NFR modeling require shape factor in the matrixfracture transfer function,
although, some approaches do not require shape factor (Firoozabadi
and Thomas, 1990). Shape factor is commonly used and it is a crucial parameter
in matrixfracture transfer function. Warren and Root (1963)
have defined rectangular shape factor as:
where, n refers to the number of fractures sets and:
In subsequent years, Kazemi et al. (1976) developed
new shape factor for their simulator using finite difference method. Their shape
factor for rectangular geometry is:
Ueda et al. (1989) research concluded that Kazemi’s
shape factor needed to be multiplied by a factor of 2 or 3 in order to get more
realistic pressure distributions. Their works were later supported by Lim
and Aziz (1995), which showed that Kazemi’s shape factor needed to
be adjusted with factor of ~2.5.
Coats (1989) derived shape factor that are doubled
of Kazemi’s shape factor. Coat’s shape factor for rectangular geometry
is:
The method used by Coats (1989) is Fourier finite sine
transform and integration. Fourier transformation was also used by Chang
(1993) and Lim and Aziz (1995) to arrive at another
shape factor different from Coats (1989). Coat’s
work became the main references for both of them. They continued Coat’s
work but with different boundary conditions. By using pressure boundary conditions,
they arrived at similar shape factor for rectangular geometry Eq.
8:
Lim and Aziz (1995) added that the total amount of
mass entered a system at time t, M_{t} with corresponding mass after
infinite time, M_{8} can be expressed as in Eq. 9.
In additions, the matrixfracture transfer rate can be expressed as in Eq.
10:
Equation 1113 are the analytical solutions
given by Lim and Aziz (1995) for single phase flow in
fracture. The solutions can be differentiated with respect to time and related
with Eq. 10 to obtain respective shape factor:
Meanwhile, Chang (1993) has derived another shape factor
using constant flow rate boundary conditions which are:
On the other hand, Quintard and Whitaker (1995) used
assumption of infinite permeability in the fracture to set the boundary value
problem for double porosity flow. By solving using Fourier series for rectangular
geometry, they reached conclusion of shape factors which are:
Table 1 summarized all the shape factor for rectangular geometries.
Table 1: 
Shape factors for rectangular geometry 

Table 2: 
General data used in comparison study 

COMPARISON OF SHAPE FACTORS
The reservoir problem selected for this comparison study was a multiphase depletion
run with 5x3x2 blocks and only one production well at (1,1,1). The 5x3x2 blocks
were selected for this purpose to allow 3ways flow between the blocks during
simulation. There is no injection well and only one production well. This is
to avoid complication of the problems which later would complicate the results
analysis. The production runs with no flow constraints. Table
2 describes the details of the reservoir problem. The additional reservoir
properties used for this comparison are from the Sixth SPE Comparative Solution
Project (Firoozabadi and Thomas, 1990).
COMPARISON APPROACH
Lim and Aziz (1995) has provided analytical derivation
of shape factors for single phase flow.
In the effort of deriving shape factors, they have shown that the total amount
of mass entered a system at time t, M_{t } and the corresponding mass
after an infinite time, M_{8 }can be expressed as in Eq.
9. The Equation 9 is known as a dimensionless pressure,
P_{D}. This dimensionless pressure is a function of dimensionless time
as shown in Eq. 18.
The equivalent fractures length, L_{e} is given in the Eq.
34:
and the analytical expressions for single phase flow is shown in the Eq.
1113:
It is desired to know the effect of different shape factors in multiphase flow
NFR simulation. The comparison is done by representing the simulation results
in P_{D} and t_{D}. A basic double porosity simulator is used
to solve the reservoir problem. The simulator solves the pressure and saturation
by using Implicit Pressures Explicit Saturation (IMPES) method. When P_{D}
is plotted against t_{D}, the gradient (P_{D}/t_{D})
gives indication of the matrixfracture transfer rate. From Eq.
19, it is shown that the matrixfracture transfer rate is proportional to
the gradient (P_{D}/t_{D}). Eq. 19 can be
found by extending the analytical solution given by Lim
and Aziz (1995). The derivation new correlation is shown in the appendix:
This relation is for analyzing the results. In additions, the results are compared
against the Lim and Aziz (1995) analytical solutions.
Single set of fractures: Shape factor for single set of fractures is
used when fractures are assumed to be present only in one direction (Fig.
2). The flow between matrix and fracture is assumed to be in a single direction.
The direction of fractures is not necessary has to be in xaxis, but it can
be at any axis such as straight yaxis or slanted xy axis. This also applies
to double and triple set of fractures.
The comparison result for single set of fractures is presented in Fig.
5. By applying Eq. 19, it is deduced that steeper gradient
shows higher matrixfracture transfer rate. At t_{D}<0.5, Kazemi’s
shape factor has showed much lower transfer rate as compared with other researcher’s
shape factor. It also require twice the t_{D} value to reach the same
P_{D} value. Note that all the curve gradient is very high at early
t_{D} and becoming lower as t_{D} increases. This indicate that
flow is high at early time and decreasing with time. This curve behaviour is
associated with rate of change of flow Eq. 19:

Fig. 2(ab): 
(a) Matrix block with single set fractures and (b) Flow boundary
between matrix and fractures 

Fig. 3(ab): 
(a) Matrix block with double set fractures and (b) Flow boundary
between matrix and fractures 
Conversely, Warren and Root (1963) shape factor yield
the highest Δgradient. The Warrenand Root (1963)
shape factor flow is the highest at early t_{D} and then, it dropped
quickly. In general, all the results converges into P_{D }and the curve
gradient becoming linear. This indicate that the flow is becoming steady. When
the gradient is 0, it literally means that no flow between matrixfracture and
it occurs at P_{D} = 1. At P_{D} = 1, the p_{m} must
be equal to p_{f }since the initial pressure, p_{i} always constant.
This supported by Eq. 1 whereby there must be a pressure difference
to initiate the flow.
As stated earlier, analytical solution from Lim and Aziz
(1995) is based on single phase flow and direct comparison cannot be made.
However, it can be observed that the analytical solutions has much higher Δg
as compared with others. This denotes that the analytical solutions flow is
very high at initial t_{d} and then decreases rapidly. The reservoir
problem is a multiphase flow problem whereby the transfer of fluids are much
more complex. The components that present in a multiphase flow is water, oil
and gas. Total multiphase matrixfracture transfer rate is a summation of all
the components, while single phase matrixfracture transfer rate is having only
the oil components. Analytical solution cannot be compared directly with results
of other shape factors but it can serves as a reference line. This can be used
to detect shape factor results that yield faster transfer rate by comparing
the gradient and P_{D}.
Two and three set of fractures: Shape factor for two sets of fractures
is used when fractures assumed a matchstick model (Fig. 3)
and three sets of fractures is used when fractures assumed sugarcube model (Fig.
4).

Fig. 4(ab): 
(a) Matrix block with triple set fractures (b) Flow boundary
between matrix and fractures 

Fig. 5: 
Single set of fractures 

Fig. 6: 
Two sets of fractures 
The comparison for two and three set of fractures are presented in Fig.
6 and 7, respectively. Both the results can be analyzed
using the same approach discussed for the parallel plate model.

Fig. 7: 
Three sets of fractures 
It is interesting to note that at early t_{D}, Warren
and Root (1963) shape factor yields the highest flow as compared with the
analytical solution.
In general, Warren and Root (1963) has the highest
rate of change of flow while Kazemi et al. (1976)
shape factor produces the lowest rate of change of flow. These two models can
be taken as the two extreme bounds where (Quintard and Whitaker,
1995; Chang, 1993; Lim and Aziz,
1995 and Coats, 1989) shape factors produce intermediate
rate of change of flow.
It is apparent that higher shape factor will results in high rate of change
of flow. The flow behavior is such that there is an initial high flow between
matrixfracture and then followed by quick drop of flow rate. Figure
8 is an example of 15x1x1 line injection model investigated using different
shape factors.
It showed that the oil production using shape factor from Warren
and Root (1963) has the highest productivity before the fifth year, followed
by a quick drop in the productivity index. Kazemi et
al. (1976) shape factor resulted in the lowest productivity among the
models in the initial time but the production picked up after the fifthyear.
Not surprisingly, oil production using shape factors from Coats
(1989) and Lim and Aziz (1995) are bounded in between
the results by Warren and Root (1963) and Kazemi
et al. (1976).
CONCLUSION
Comparison of shape factors is presented in dimensionless parameters, P_{D}
and t_{D}. Gradient P_{D}/t_{D} are proportional to
the matrixfracture transfer rate. Higher gradient P_{D}/t_{D}
indicate higher transfer rate. The dimensionless comparison displays the different
flow behavior of different shape factors.
NOMENCLATURE
t 
= 
time, T 
L 
= 
length, L 
V 
= 
volume, L^{3} 
q 
= 
matrixfracture transfer rate, L^{3}/T 
Q 
= 
matrixfracture transfer rate, M/(L^{3}T) 
p 
= 
pressure, M/(LT^{2}) 
k 
= 
absolute permeability, L^{2} 
σ 
= 
shape factor, 1/L^{2} 
μ 
= 
viscosity, M/(LT) 
n 
= 
set of fractures 
φ 
= 
porosity, fraction 
c 
= 
compressibility, fraction 
ρ 
= 
density, M/L^{3} 
SUBSCRIPTS
i 
= 
Initial 
t 
= 
Total 
m 
= 
Matrix 
∞ 
= 
Infinity 
f 
= 
Fracture 
d 
= 
Dimensionless 
APPENDIX
Derivation of dimensionless correlation: Equation 13
can be differentiated with time to obtain Eq. A1:
By using finite approximations for first degree derivative, Eq.
A1 can be rewritten as Eq. A2:
For Sufficient small points and if p_{f} = Constant, then:
Substitute Eq. A3 into the given Eq. 10,
we get:
If:
is constant/equivalent, while t varies, then Eq. A2 can be
rewrite into Eq. A5:
Substitute Eq. A6 into Eq. A5, we get
Eq. A6: