Reservoir Geomechanics with ABAQUS - Simulia

15 Pages·2004·0.4 MB·English

Reservoir Geomechanics with ABAQUS
B. Bostrøm and E. Skomedal
Statoil ASA, Norway
Abstract: The coupled hydro-mechanical behavior of two North Sea high-pressure/high-
temperature gas-condensate fields during depletion is studied using the soils consolidation
procedure available in ABAQUS. Time dependent reservoir pressure fields obtained from three-
dimensional black oil reservoir models were transferred into finite element reservoir
geomechanical models in ABAQUS as three-phase flow in deforming reservoirs cannot be
performed within ABAQUS.
Keywords: Petroleum Engineering, Reservoir, Coupled Analysis, Geomechanics, Failure.
1. Introduction
With reservoir pressures of about 80 – 90 MPa and temperatures of 150 - 170 degrees C the two
gas condensate fields Kristin and Kvitebjørn are characterized as high-pressure/high-temperature
fields (HP/HT). The fields are located offshore mid Norway. Production start is planned in 2004 –
2005.
The Kristin field comprises of two main reservoirs, Garn and Ile, separated by intra reservoir
shale. Kvitebjørn on the other hand is highly faulted. The possible compartmentalization of the
reservoir associated with the faults is a major uncertainty in developing such a field. Production
from the reservoirs will be obtained by huge depletion as no gas or water injection is planned.
Tremendous effective stress changes are therefore expected, which have an influence on the
various aspects of field development and planning. This necessitates a reservoir geomechanical
study.
The focus of the models presented for the two HP/HT gas–condensate fields are different. A full-
field model is presented for the Kristin field, while the results of a detailed analysis of an internal
fault is given for the Kvitebjørn field. The full-field or global model addresses several issues
important for production, like reservoir compaction, subsidence, well casing integrity and infill
drilling options. This case also focuses on the coupling between the standard reservoir simulator
ECLIPSETM (trademark of Geoquest Schlumberger) and ABAQUS that is necessary due to the
above-mentioned shortcoming in the coupled pore fluid flow and stress analysis procedure in
ABAQUS.
2004 ABAQUS Users’ Conference 117
The local analysis on the other hand is carried out to investigate possible fracturing (shear or
tensile) of the formation through and around faults, crushing, reactivation of faults and changes in
the sealing capacity of the interbedded shale layers during depletion, possible top seal leakage
through faults and fractures.
Finally this paper addresses the necessary enhancements of ABAQUS, so that it will be of even
greater value for the reservoir modelers.
2. Full-field model of the Kristin field
2.1 Geometry and model
A 2D plane strain model is sufficient due to the large extension of the reservoir normal to the
plane of interest. The width of the reservoir at studied section is 4166 m, while the reservoir height
is typically 160 m adding the thickness of the Garn and Ile reservoirs. Our geomechanical model
includes in addition the over-, side- and under-burden. The total width of the model is 28 500 m
and the total height is 10 000m.
The geometry of the reservoir layers is taken from a vertical cross section of the 3D ECLIPSE
model. ABAQUS/CAE allow you to import parts from CAD system using industry-standard
formats. Unfortunately, ECLIPSE does not support these formats, implying that ECLIPSE
geometry printed in ASCII format must be transferred into a CAD format using a translator.
Several minor faults have been mapped in the two main gas-condensate reservoirs. A limited
number of drained triaxial tests on faulted reservoir material indicate that the strength and stiffness
of these faults are higher than the surrounding material. Internal faults (barriers) are therefore
modeled by a rough, small sliding contact formulation. In this way there are no relative motion
between the surfaces, while pressure may still be discontinuous across the fault. Faults outside the
reservoir are not included in the present model for the same reasons.
The entire finite element mesh is shown in Figure 1.a. The bottom boundary, the left and right
vertical boundaries are fixed. Displacement nodes at the top (seabed) boundary are free. The
model consists of a total of 41735 nodes and 20639 elements. Reservoir mesh details are shown in
Figure 1.b and c.
Different analysis procedure and matrix strength criterion are used for the different parts of the
model. The reservoir layers use the pore fluid diffusion/stress analysis procedure, while undrained
(no pore fluid flow) behavior is assumed for the surrounding shale. The 6-node triangular
displacement and pore pressure element CPE6MP is used for the reservoir layers, while the 6-node
triangular displacement elements CPE6M is used for the surrounding shale.
2.2 Material
The model is divided into 8 material layers: The hydrocarbon bearing sandstone formations Garn
and Ile, the underlying water bearing sandstone Tofte, the intra reservoir shales Not and Ror and
the shales in the overburden, sideburden and the underburden.
118 2004 ABAQUS Users’ Conference
Reservoir sandstone stiffness and strength deduced from standard triaxial laboratory testing are
made porosity dependent by a field dependent parameter, where the initial field parameter
distribution is taken equal to the distribution of the initial porosity in the reservoir.
Geomechanical input data for the shale material are taken from sonic logs and existing
laboratory tests on shale samples from the exploration wells at the field. Further recall
that a displacement finite element formulation (element CPE6M) is chosen for the shale
material. This implies that undrained stiffness must be calculated. It may be shown that
the undrained Poisson’s ratio, ν , and the undrained Young’s modulus, E , are related to
u u
the drained modulus through the following expressions
( )( ) Eν
1−2ν 1+ν +
K
ν = a
u ( )( ) E
21−2ν 1+ν +
K
a
( ) E
31−2ν +
E K
u = a
E ( )( ) E
21−2ν 1+ν +
K
a
Where of E is the drained Young’s modulus and ν is the drained Poisson’s ratio. The
apparent compressibility of the pore water, K, is given as
a
K
K = w
a n
where K is the bulk modulus of the pore water and n is the porosity. The spatial varying
w
stiffness is controlled by an artificial temperature distribution, i.e. temperature dependent
stiffness.
2.3 Prescribed conditions
The initial effective stresses are generated in the ABAQUS subroutine, SIGINI.
The initial geostatic stress field must be in equilibrium with the applied loads and boundary
conditions. Ideally, the loads and initial stresses should exactly equilibrate and produce zero
deformations. This state is obtained performing an initial ABAQUS analysis fixing all
displacements degree of freedoms. Calculated reaction forces written to the ABAQUS output file
are then used to create nodal point forces, which are applied in the first step of the actual
ABAQUS analysis.
2004 ABAQUS Users’ Conference 119
The pore pressure depletion history within the reservoir is transferred from the ECLIPSE reservoir
model to ABAQUS utilizing the user subroutine DISP. A file containing the pore pressure in each
ECLIPSE block is read for each time step analyzed by ABAQUS.
The initial porosity distribution is transferred from the ECLIPSE reservoir model to ABAQUS
utilizing the user subroutine VOIDRI for reading initial porosity in ABAQUS. The initial void
ratio e , defined as the ratio between the pore volume and the solid volume, is related to the
0
porosity n through: e =n/(1-n).
0
2.4 Results and discussion
Contour plots of applied pore pressure reduction, displacements, stress path parameter and xy-
plots of deformation along well paths are presented.
2.4.1 Pore pressure
Contour plot of the applied pore pressure reduction for year 2025 is shown in Figure 2. The
maximum reduction is slightly more than 80 MPa in the western part of Garn at the end of
production. The minimum reduction at the end of production is about 30 MPa in the eastern part
of Ile. The maximum differential pressure across barriers is slightly smaller than 20 MPa. The
maximum differential pressure occurs in Ile between year 2009 and 2014. The maximum
differential pressure in Garn is about 16 MPa.
2.4.2 Displacement
Horizontal and vertical displacement profiles at years 2007, 2009 and 2025 are plotted in Figure 3.
Three profiles are chosen and the placement of these is indicated in Figure 1.b.
The predicted seabed subsidence is in general very modest. It is equal to 0.7 cm above the central
part of the reservoir in the year 2025 (Figure 3.b), and it increases slightly towards the reservoir
flanks. Reservoir compaction on the other hand is largest in the western parts and decreases
somewhat towards east. It is equal to 45 cm in the west, 38 cm in the central parts and 27 cm in
east in the year 2025. The overburden expansion is a Poisson’s ratio effect, which is amplified due
to the undrained behavior of the shale (ν is larger than 0.47 in the overburden). The horizontal
u
contraction of the reservoirs may be deduced from Figure 3.a and c. The resulting horizontal
contraction is equal to 107 cm in the year 2025.
2.4.3 Stress path
Stress changes during depletion is a key issue for well planning after start up of production as the
decline of the fracture pressure due to reservoir depletion is deduced form this. A common
parameter to monitor the stress path during production is the stress path coefficient, defined as the
ratio of the change in horizontal stress ∆σ to the change in pore pressure ∆p,
h
∆σ
γ= h
∆p
120 2004 ABAQUS Users’ Conference
The calculated spatial variation of the stress path coefficient in the year 2025 is relative
homogeneous. It is less or equal 0.65 in the two main reservoirs as shown in Figure 4.
Similar numbers found using an established analytical model is 0.82 (uniaxial strain model), i.e.
less favorable with regards to infill drilling. The advanced numerical model includes the effect of
bounding shale (arching) and rock nonlinearities.
2.4.4 Well casing damage
Planned well paths are projected into the studied cross section. Displacements in the horizontal
and vertical direction along this well path are shown in Figure 5.a. The axial strain resulting from
these displacements are shown in Figure 5.b. It is seen that the maximum axial strain along this
path is about 4 mStrain in the year 2025. This strain is larger than the strain (of about 2 mStrain)
that normally causes yielding in the steel. The maximum strain takes place in the lower part of the
Ile formation with high porosity. It is however possible that the full-field model is too coarse to
capture localized deformations within the Ile formation. This means that the axial strain may
locally be even slightly larger.
It is also seen that the strains in the well outside the reservoir layers are small. This means that the
uncertainties in the calculated strains of significance first of all are affected by uncertainties in the
drained bulk modulus of the sandstone. Finally, the change in curvature is too small to affect the
deformations of the well.
3. Detailed model of the Kvitebjørn field
3.1 Geometry and model
A detailed 2D plane strain ABAQUS model of two segments with 40 m relative throw and 40
degrees fault dip is shown in Figure 6. The model consists of the 23 layers defined in the full-
field-reservoir-model established in ECLIPSE. The thickness of the reservoir layers is 175 m. The
model includes a 250 m thick impermeable shale layer above and below the reservoir.
Roller boundaries are assumed along the vertical boundaries, the bottom boundary is fixed
vertically and horizontally and the top boundary is loaded with a constant vertical stress equal to
the in-situ vertical stress of the model. The in-situ stress field is assumed constant throughout the
model. The effective overburden stress is equal to σ′ is equal to 11.5 MPa, while the effective
v
horizontal stresses σ′ = σ′ are equal to 7.5 MPa. The loads and initial stresses equilibrate exactly
h H
and produce zero deformations as expected.
Again the 6-node triangular displacement and pore pressure element respresents the reservoir
sandstone in addition to the fault zone. A linear elastic ideal plastic Mohr Coulomb model gives
strength. The shale material in burdens is represented by 6-node triangular displacement elements.
A linear elastic ideal plastic Mises model gives the strength.
The initial void ratio is an important parameter as the stiffness and strength of the reservoir
sandstone depends on this value.
2004 ABAQUS Users’ Conference 121
The fault is represented by continuum elements. Pressure within an internal fault varies linearly
between the pore fluid pressures on each side of the fault. The top and bottom of the fault zone are
assumed to have the strength of the shale, while the central part where the two reservoir segments
are in contact is here assumed to have identical stiffness and strength as the reservoir sandstone
layer with the lowest porosity.
3.2 Results and discussion
The ultimate pressure reduction for a reservoir segment is approximately 60 MPa, starting up at
77.5 MPa and ending up with 18 MPa. Differential pressure across the fault is studied in Case A
and B. The left segment is depleted in Case A, while unchanged pore fluid pressure is assumed
for the right segment. Opposite assumptions for the two reservoir segments are assumed in Case
B. Both the left and the right segment are depleted in Case C.
3.2.1 Case A – depletion of left segment
Evolution of tensile failure and compaction failure in the form of stress paths for selected points
within the reservoir sandstone layers is shown in Figure 7. The placement of the three points
chosen are indicated by the letters A, B and C in Figure 6. Note that the final pressure drawdown
is equal to 100 MPa in this plot.
Possible compaction failure is predicted for a typical stress path for the depleted left reservoir
segment (curve A)). The stress ratio ∆K defined as ∆σ′ /∆σ′ is equal to 0.49 and the stress path
3 1
coefficient γ defined as ∆σ /∆p is equal to 0.58 for this point. A limited increase in the degree of
3
shear mobilization is obtained for the fault zone (curve B), while tensile failure is experienced for
the right reservoir segment with unchanged pressure (curve C).
The degree of shear mobilization, defined as the ratio of the mobilized friction tanρ and the
friction at failure, tanφ, is plotted in Figure 8 for the left reservoir segment. All values are less than
0.6 after 60 MPa depletion. This implies that the state of stress at any point in the depleted
segment is situated far below the Mohr Coulomb failure envelope. Focus in the rest of chapter will
therefore be on the behavior of the reservoir surrounding shale and the fault zone.
Figure 9 gives the evolution of plastic yielding expressed by the effective plastic strain in the
burdens and fault. The effective plastic strain is given by
2
εp = εplεpl
3 ij ij
Reservoir segments are colored grey in this plot, i.e. they are not displayed. Yielding is initiated in
the shale below the acute corner of the depleted segment at 26 MPa depletion. A continuous
failure zone between the left at right reservoir segment is established after 28 MPa depletion, see
Figure 9.a. Shale failure above the top corner of the depleted reservoir is initiated at 54 MPa
depletion. Failure in the intermediate fault zone is not initiated. An effective stress increase in the
fault zone is actually specified due to the assumption of linear pressure variation through this zone.
122 2004 ABAQUS Users’ Conference
3.2.2 Case B - depletion of right segment
Contour plots of the effective plastic strain in burdens and fault zone are shown in Figure 10.a and
b. Yielding is initiated above the acute corner of the depleted reservoir segment at 28 MPa
depletion according to Figure 10.a.
3.2.3 Case C - depletion of both segments
Effective plastic strain contours at 60 MPa depletion are shown in Figure 11. Plastic flow is
concentrated to the shale part of the fault zone. This is believed to be a consequence of the free
boundary at the top of the model. Uniaxial vertical strain is observed for large parts of the model
and this is reflected in the limited plastification. A fixed boundary at the top will result in high
shear stresses in the shale above and below the depleted reservoir segments.
3.2.4 Concluding remarks for internal faults
The local finite element analyses of internal faults with ABAQUS predict
• Most likely increased communication between two reservoir segments isolated by an
internal fault due to large pressure difference across the fault
- no tensile failure or shear failure in the fault zone
- however, there may be communication through the shale at the top and bottom corner of
the depleted segment
• Top-seal leakage through faults and fractures is not expected
A starting point for incorporating these results in ECLIPSE is to search for internal faults with
differential pressure above 28 MPa, and introduce communication between these at the acute angle
of the segment that is most depleted. This search may either be done manually or by a script.
However increased communication between two segments are not a function of the differential
pressure alone. The absolute pressure certainly also play a role as shown by Case C. It is also
believed that the results are sensitive to fault geometry expressed by the dip. The ABAQUS results
so far (Case A and B) indicate that the differential pressure necessary to lead to increased
communication through the bounding shale decreases with decreasing dip, as plastic yielding is
initiated in the bounding shale close to the acute angle of the depleted segment.
4. Necessary enhancement of ABAQUS
Undrained (no pore fluid flow) total stress analyses are carried out for the shale material. Strength
is given by the Mises criterion. However, undrained effective stress analysis is the preferred
alternative for undrained rock behavior. Zero permeability pore pressure elements may be used for
this purpose. This will be quite cumbersome as only continuous pore pressure elements are
available in ABAQUS. The mesh must be “opened up” at the interface between the reservoir
sandstone and the surrounding shale as pressure degree of freedoms should not be common at this
interface.
2004 ABAQUS Users’ Conference 123
Model generation will be simpler if discontinuous pore pressure elements were available in
ABAQUS as these could be used to model undrained shale behavior. Zienkiewicz et al (1983)
introduced the nine-noded displacement and linear pressure element for 2D application. High
accuracy was obtained for incompressible fluid mechanics problems with this element type. In the
compressible case pore pressure degrees of freedom may be eliminated at element level, i.e.
resulting in a displacement element. Naylor (1974) exploded this formulation, showing how a
standard displacement based finite element program can be modified so that stresses are separated
into pore fluid at rock matrix components. Reduced integration was advocated in the nearly
incompressible case. The computational costs related to this formulation are considerably less than
the discontinuous pressure elements. The K -formulation by Naylor (1974) is widely used in
w
specialist geotechnical finite element programs.
The implication of tensile or shear failure may be studied in more detail introducing interface
elements (discrete crack analysis) at critical places. However, all interface contact elements was
removed in ABAQUS version 6.1. Instead one should use the *CONTACT PAIR option with
SMALL SLIDING parameter to model infinitesimal sliding and pore fluid interaction between
deformable bodies. Unfortunately this new option do not account for tangential gap flow, which is
essential for examining possible hydraulic communication between two isolated reservoir
segments. From our standpoint it seems that reintroducing the interface elements again in
ABAQUS is a good idea.
The user subroutine option in ABAQUS makes the program to a very versatile tool. It is even
possible to define your own finite element in ABAQUS. However developing and maintaining
such computer code is not a part of our strategy. So for these topics we hope that the ABAQUS
organization see them as useful suggestions. Material routines on the other hand may be
something it is easier to get accept for internally to spend time on developing.
5. Conclusions
Results from the studies of two North Sea HP/HT gas-condensate fields have been presented here.
The full-field simulations of the Kristin field focus on reservoir compaction, subsidence, and
possible well casing damage and stress changes during depletion. The evaluation of the well
casing integrity is based on axial strain along the well path. Stress changes during depletion is a
key issue for well planning after start up of production as the decline of the fracture pressure due
to reservoir depletion is deduced form this.
The local Kvitebjørn model on the other hand focuses on the mechanical behavior of faults and
barriers during depletion. A realistic geometry, with interbedded shale layers and inclined fault
zone was established. Evolution of shear failure as a consequence of increasing differential
pressure across the fault zone was analyzed and the limiting value for start of cross between two
separated reservoir segments was found. This contributes to better insight of the flow field during
the production history.
The full-field analysis of Kristin utilizes an external coupling between the reservoir simulator
ECLIPSE and ABAQUS that has been established. Initial pressure, initial porosity variation and
simulated pressure history are transferred from ECLIPSE to ABAQUS utilizing the user
124 2004 ABAQUS Users’ Conference
subroutine facilities in ABAQUS. ECLIPSE reservoir geometry is imported into ABAQUS/CAE
reading AutoCAD files. Later experience indicates that a better approach is to incorporate
geometry information directly in the Python script file that generates your ABAQUS model.
Both commercial and research simulators that take the coupled nature of three-phase-flow and
deformation into account exists today. It is then natural to ask if this will make the partly coupled
approach described here superfluous? Our experience is however that a partly coupled approach
between a conventional reservoir simulator and a stress simulator is the best approach for the near
future when advanced geomechanical issues must be taken into account. The partly coupled
approach benefits from the latest developments in physics and numerical techniques of both
simulators. Large-scale 3D simulations, which a natural next step for us is another issue that
favors the partly coupled approach, as it is believed that computational speed will be crucial.
Summarizing our experience so far with ABAQUS for reservoir geomechanical simulations we
may state that ABAQUS is very well suited for reservoir geomechanical simulations. However
with some “minor” enhancements as pointed out previously ABAQUS will be of even greater
value for the reservoir modelers.
6. Acknowledgements
The management of Statoil’s Kristin and Kvitebjørn projects are thanked for permission to publish
this paper. Hans Petter Jostad at Norwegian Geotechnical Institute is acknowledged for valuable
support on the Kristin analysis.
7. References
1. Naylor, D.J., “Stresses in Nearly Incompressible Materials by Finite Elements With
Application to the Calculation of Excess Pore Pressure,” I.J.N.M.E., Vol. 8, pp. 443-460,
1974.
2. Zienkiewicz, O.C, Taylor, R.L. and Baynham, J.A.W., “Mixed and Irreducible Formulations
in Finite Element Analysis,” Chapter 21 of Hybrid and Mixed Finite Element Methods (eds
S.N. Atluri, R.H. Gallager, and O.C. Zienkiewicz), pp. 405-31, Wiley, 1983.
2004 ABAQUS Users’ Conference 125
WWWeeesssttt MMMiiiddddddllleee
EEEaaasssttt
44116666 mm

Garn
Ile

Figure 1. FE Mesh. (a) Entire model. (b) Close up of the two main gas-condensate
reservoirs and (c) Reservoir details in the west.

Figure 2. Pore pressure depletion in year 2025, i.e. after 20 years of production.
126 2004 ABAQUS Users’ Conference