Gas lift simulation and experiments in conjunction with the Lyapkov P. D. methodic

(1)

Faculty of Science and Technology

MASTER’S THESIS

Study program/Specialization:

Petroleum Engineering/ Natural Gas Engineering

Spring semester, 2016 Open

Writer:

Emil Gazizullin ………

(Writer’s signature) Faculty supervisor:

Rune W. Time Co-supervisor:

Hermonja A. Rabenjafimanantsoa Thesis title:

Gas lift simulation and experiments in conjunction with the Lyapkov P. D. methodic

Credits (ECTS): 30 Key words:

Nodal analysis Multiphase flow Gas lift optimization PVT properties

Pressure distribution curve

Pages: 130

+ enclosure: including Stavanger, 15.06.2016

Date/year

(2)

ii

Acknowledgements

This Master’s Project was a great final of my study at the University of Stavanger. A lot of new skills and knowledge were gained here.

I consider it an honor to work with the Professor Rune Wiggo Time, who have helped me despite his busyness. I am pleased to acknowledge his helpful comments and suggestions without which my work would not go further.

I owe my deepest gratitude to the Senior Engineer Hermonja A. Rabenjafimanantsoa also known as Benja for encouragement in various ways, valuable guidance and time spent for me.

I also thank the UiS engineer Svein Myhren, who spent his valuable time helping me to deal with the hardware and software for the experiments.

I express many thanks to the lab assistant Nikita Potokin for his assistance, without whom my work would significantly slow down.

I would like to thank Rinat Khabibullin, Gazprom research center worker who gave me a solid foundation and understanding of programming skills that formed a basis for my current work.

Finally, all my gratitude to my mother for her encouragement, without her support I would have never achieved my goals.

(3)

Abstract

The present Master’s thesis reports the study of continuous gas lift for a production oil well.

The objective of the present work is a better understanding of processes occurring in a gas lift well, as well as optimization of the gas lift system. Numerical simulation for a gas lift well was carried out by means of Matlab. The program and a step-by-step guide were developed on the basis on nodal analysis, with implementation of the Lyapkov’s methodic; the program calculates pressure distribution curves along with local mixture properties along the wellbore in tubing and annulus space for steady state flow of multiphase mixture. The gas lift simulation has been performed for different gas injection rates to see how the liquid production will increase.

For experimental purposes, a 5 m experimental loop was used. All experiments were carried out with tap water as reservoir fluid and compressed treated air from the atmosphere as injection gas. A range (from 0,01 to 5 standard liters per minute) of different gas injection rates were chosen to observe a shape of the liquid production curve.

Finally, experimental results are to be confirmed by calculations, based on Bernoulli’s equation for the experimental setup.

During the experiments it was a great luck to install and use a new gas flowmeter provided by Alicat Company and a Sensirion SQT-QL500 liquid flowmeter.

(4)

iv

List of tables

Table 1- Well input data ... 25

Table 2- Input data ... 29

Table 3- Output data ... 30

Table 13- Called functions ... 37

Table 22- Continuation ... 42

Table 34- Choice of mixture type and flow regime ... 49

Table 39- Choice of true gas fraction ... 51

Table 40- Choice of true fraction share of water and oil in the liquid ... 52

(7)

vii

Table 50- Pressure and temperature at the outlet of i stage ... 57

Table 58- Liquid flow summary measurements for experiment 1 ... 128

Table 59- Liquid flow summary measurements for experiment 2 ... 129

(8)

viii

List of figures

Figure 1- Flow regimes and flow regime map in vertical two-phase flow. ... 16

Figure 2- The production system of a continuous flow gas lift well (reproduction) ... 18

Figure 3- Open tubing flow installation with gas injection at the tubing shoe ... 19

Figure 4- Summary flowchart of the program ... 24

Figure 5- Structure of the program ... 27

Figure 6- Gas saturation versus pressure at standard pressure and reservoir temperature for Karakuduk oilfield (reservoir J-1/2)... 30

Figure 7– Oil volume factor versus pressure at reservoir temperature for Karakuduk oilfield (reservoir J-1/2) ... 31

Figure 8–Oil density versus pressure at reservoir temperature for Karakuduk oilfield (reservoir J-1/2) ... 32

Figure 9– Oil viscosity versus pressure at reservoir temperature for Karakuduk oilfield (reservoir J-1/2) ... 33

Figure 10- Liquid viscosity change due to temperature according to Lewis- Squires ... 34

Figure 11– Choice of the surface tension between the liquid phase as the outer phase of the stream and the gas ... 40

Figure 12– Surface tension vs watercut for Karakuduk oilfield (reservoir J-1/2) ... 41

Figure 13- Choice of liquid viscosity ... 43

Figure 14- Water viscosity versus watercut for Karakuduk oilfield (reservoir J-1/2) ... 43

Figure 15- Flowchart for PDC along annulus calculation ... 60

Figure 16- Flowchart for PDC along tubing calculation ... 63

Figure 17- Pressure distribution curves along tubing and annulus ... 67

Figure 18- Pressure and temperature gradients along the well ... 68

Figure 19- PVT properties of oil in the well along tubing and annulus ... 69

Figure 20- PVT properties of gas and gas-oil-water mixture in the well along tubing and annulus ... 70

Figure 21- Fluid viscosity versus watercut ... 71

Figure 22- Different fluid properties in the well along tubing and annulus ... 72

Figure 23- Different fluid properties in the well along tubing and annulus (continuation) ... 73

Figure 24 – PDC along tubing and annulus when no gas lift in the well is installed ... 74

Figure 25 - Liquid production and gas lift operating pressure vs. gas injection rate... 75

Figure 26- Illustration of the gas-lift model ... 77

Figure 27- Flowchart of the data logging setup ... 78

Figure 28- LabVIEW interface ... 81

Figure 29- Sensirion program interface ... 82

Figure 30- Flow pattern at 0,2; 0,5; 1,0; 1,5; 2,0; 2,5; 3,0 SPLM of gas (air) respectively ... 83

Figure 31- liquid rate as a function of gas injected into the system ... 86

Figure 32- differential pressure as a function of gas injected into the system ... 87

Figure 33- Bottomhole pressure as a function of gas injected into the system ... 88

Figure 34- Alicat absolute gas pressure as a function of gas injected into the system ... 89

Figure 35- Relative errors for different measurements ... 90

(9)

ix

Figure 36-Section of the pipe, to which Bernoulli’s equation is applied ... 92

Figure 37-Comparison of calculated values of liquid flowrate and the factual value, obtained during the experiments ... 93

Figure 38- Comparison of calculated values of liquid flowrate and the factual value, obtained during the experiments as a function of differential pressure ... 94

Figure 39- Alicat gas flowmeter, mounted on top of the loop ... 124

Figure 40- Atmospheric pressure gauge ... 124

Figure 41- Valve 2 directing gas towards model ... 124

Figure 42- Switch directing gas towards open atmosphere; gas regulator, and inlet pressure Crystal Digital test; Gauge XP2i manometer ... 125

Figure 43- Bottomhole Crystal Digital test Gauge XP2i manometer and the point of pressure measuring and Sensirion SLQ-QT500 liquid flow meter ... 125

Figure 44- Rosemount 3051C pressure gauge differential pressure measurer ... 125

Figure 45- Open position of Valve 1 ... 126

Figure 46- A Norgren gas drier ... 126

Figure 47-Operator’s workplace ... 126

Figure 48- Ball valve on the lower horizontal pipe ... 127

Figure 49- LabVIEW working scheme ... 127

Figure 50- Flow rate values from Sensirion SQL-QT500 liquid flowmeter for experiment 1 ... 130

(10)

Introduction Gazizullin E. 10

1 Introduction

The author (myself) has created a program for an electrical submersible pump (ESP) selection for a vertical inclined onshore well. This program is based on the Lyapkov’s method, which was developed in 1987 and was widely used in the USSR. The method has been reconsidered and implemented in terms of the nodal analysis. It means that the system “well-near wellbore area- reservoir” is considered as a sum of pressure gradients that are easy determined, when the following parameters are known: well inflow, depression and PVT properties of the fluid. The program was implemented on the basis of VBA in MS Office Excel.

Project work

The objective of the following master thesis is to create a program, which is able to calculate a pressure profile along the well in case of a gaslift well as well as its optimization. The program was realized in Matlab workspace. In optimization section it was necessary to simulate an oil well with certain input parameters. It was also assumed that gas lift was installed into the well and different gas injection rates were simulated. The optimization consisted in defining the optimal gas injection rate based on the production curve. Later an experiment has to be carried out in order to prove the gas lifted well’s behavior for a certain range of gas injection rates, namely for linear section of the production curves at low gas rates.

Motivation for oil production

Gas lift takes a large share in oil production. Many oilfield wells are equipped with gas lift injecting systems. The optimal gas lift design guarantees lower operational exchanges along with high oil production rates.

Due to time limits it is not always possible to test gas lifted wells in order to determine their optimal regimes. Another question is an oil well is a sophisticated system with many variables that

(11)

Introduction Gazizullin E. 11

affect pressure drops along the well. Taking into account the fact, that in most cases reservoir fluid is presented by multiphase flow, the necessity of a program that allows to simulate a well with a continuous gas lift is a task of high importance. The program must simulate multiphase flow occurring in the well.

Many attempts have been made by researchers (Vazquez&Hernandez, 2005), (Chia&Hussain, 1999), (Bahadori, Ayatollahi & Shirazz, 2001) to develop such simulators. The main similarity is that they were developed based on the same principle, which is pressure drops calculations based on nodal analysis. The main difference is various correlations, which vary depending on fluid properties for different oil fields.

The program must be relatively simple with a user friendly interface and could be used by a regular engineer.

Similar programs that were previously developed give only pressure curves. The author [myself] has not met any program that can give a characteristic of the multiphase flow at any point along the well from the bottomhole to the surface. For studying purposes it would be interesting to implement such option, that will be represented in plots, showing well fluids’ PVT data or dynamic parameters, such as fraction velocity or flow patterns.

(12)

Theory Gazizullin E. 12

2 Theory

2.1 Gas Lift Operational Principle

When the reservoir energy is not sufficient for the well’s fluids to flow, or an engineer has desired the production rate to be greater than the reservoir energy can deliver, it becomes vital to install artificial lift into the well to provide the energy to bring the flow to the well surface («Schlumberger Well Completions and Productivity», 1999 ) .

In gas lifting natural gas that is compressed at the surface is injected in the well stream at some downhole point. Gas is injected at a certain depth from the tubing string into the flow string or casing-tubing annulus. As a sequence, the gas significantly reduces the density of the well stream and the flowing pressure losses at the point above the injection point since the major part of the pressure drop is due to potential energy change. Hence, total pressure gradient along the tubing will decrease as well, allowing the bottomhole pressure to overcome the weight of the liquid column and lift the fluid to the surface. For dead wells gas lift installation makes it possible to renew exploitation and wells with low production rates are able to produce greater liquid volumes than before. So, continuous gas lift may be implemented as the continuation of the production period (Takacs, 2005).

2.2 Multiphase flow

Multiphase flow can be observed throughout the entire production system from oil and gas reservoirs to processing facilities at the surface. The production system term in this context refers to the reservoir; the well completion; the annulus and tubing strings that connect the reservoir to the surface ; all surface facilities on land, seabed, or offshore platform; and any pipelines that deliver produced fluids to other processing facilities. The multiphase flow can be any combination of a water phase, a hydrocarbon liquid phase, and a natural gas phase (Brill,1987).

(13)

Therefore, for petroleum engineers it is important to have understanding about multiphase phenomena. The basis of any optimization problem is then pressure distribution curves with which the main parameters of fluid flow can be calculated. The precise calculation of multiphase flow parameters not only improves the engineering work, but plays a significant economic role (Takacs, 2005).

Scientists, researchers and engineers in the petroleum industry are faced with the requirement to predict the relationships for the fluids produced from a reservoir between pressure drops, flow rates, and geometry of pipes (diameter, length, angle, etc.) over the entire life of the field (Brill,1987).

2.3 Basic principles

According to Szilas (Szilas, 1985):

 For any multiphase flow problems pressure drop along the well can be affected by a great amount of variables as well as by thermodynamic parameters for each phase, parameters, describing the interaction between the phases (interfacial tension, etc.) and other factors. In correlations scientists try to reduce the number of variables by introducing non dimensional parameter groups.

 Density of the multiphase flow depends on local thermodynamic conditions. Since pressure and temperature vary in a wide range, the respective effect must be accounted for.

 Frictional pressure losses for multiphase flow are more difficult to describe compared to one phase flow, as more than one phase is in contact with well walls. Also velocity distributing in a cross-sectional area differs from a single phase laminar or turbulent flow, that adds difficulty in pressure drop calculations.

 Pressure drop significantly varies with the spatial arrangement in the pipe for different flow patterns or flow regimes. So, different correlations and procedures must be considered when calculating the pressure drop.

Before flow patterns will be discussed in details, some assumptions must be listed:

• Only steady-state flow is investigated. Most horizontal and vertical well problems involve steady-state flow,- only special cases (slugging in long offshore flowlines, etc.) necessitate the more complex treatment required for transient two-phase flow.

 The temperature distribution along the flow path is known. Although simultaneous calculation of the pressure and temperature distribution is possible, it requires additional data on

(14)

the thermal properties of the flowing fluids and the pipe's environment, and these are seldom available.

 Reservoir fluid has behavior of the black oil type, i.e. the composition of the liquid phase remains constant. In case compositional changes occur in the liquid phase along the well, black oil model cannot be applied; however, in most engineering calculations this model can be used (Szilas, 1985).

2.4 Flow regimes

Pressure gradient and the holdup of the dense phase are the important properties of two-phase flow as depend strongly on the distribution of the phases in the pipe. Thus, a great number of scientists have tried to identify the various flow regimes or flow patterns that take place in a well.

One extreme case is observed when dense liquid with low content of dispersed gas flowing in the pipe. The other extreme case is occurs with a small quantity of the liquid phase dispersed in a continuous light phase. Starting from the first extreme with introduction of an increasing amount of the light phase (gas), all possible flow patterns will be considered. Here only vertical flow patterns will be considered (Zavareh, Hill & Podia, 1988).

Depending on phase velocities and local fluid parameters, the main flow parameters had been investigated by Takacs (Takacs, 2005). These are bubble, slug, transition and mist flows. They are given in the sequence of increasing velocity. It reminds conditions in a producing oil well as the fluid goes up, higher velocities are reached. Near the bottomhole single phase fluid is present (in case bottomhole pressure is greater that saturation pressure). Going upward, more and more gas evolves from oil resulting in increased gas velocity. It increases even more due to reduction of the pressure. The continuous gas velocity increase provides all the possible flow patterns in a well.

Here is a classification proposed by Barnea (Barnea, 1987), given for a vertical well from the bottom upwards to the surface. Later, in the experimental part, frames of the flow regimes for different gas rates will be shown.

Bubble Flow

When the gas velocity is low, gas phase forms bubbles distributed discretely in the continuous liquid phase. Gas bubbles due to their lower density tend to pass the liquid, which is called gas slippage. When calculating gas density, one must include gas slippage effects. Friction loss calculations are in this case relatively simple due to contact of liquid with the walls.

The previously discussed case is valid for low liquid and gas velocities. It is usually observed in pipes with large diameter. When gas velocity is increased without changing liquid velocity, the

(15)

individual small bubbles coalesce into so-called Taylor bubbles. The flow regime increasingly changes into the slug flow. The transition may be observed at gas fractions εg > 0.25. But if liquid velocity is also increased, the process of turbulence will break bigger bubbles and will stop the coalescence of smaller ones. So even at gas fractions, higher than εg > 0.25, slug flow cannot be fully developed, however transition from bubble flow into dispersed bubble will occur.

Dispersed Bubble Flow

When liquid velocity is even higher at low gas velocities, the gas bubbles have a very small size and they are uniformly distributed in the continuous liquid phase. Gas bubbles are evenly taken by the continuous liquid phase, so the two phases rise at the same velocity without slippage.

The two-phase mixture act as a homogeneous phase. Thus, for the dispersed bubble flow regime mixture density can be easily found from the no-slip liquid holdup, since εg = 0.52.

This case is close to real single phase liquid flow as before (there is a strong contact between the wall and the liquid phase). For higher gas velocities small gas bubbles are packed so tightly that they coalesce even at high liquid velocities. The flow regime changes into the churn flow at gas volume fractions of εg > 0,52.

Slug Flow

When the continuous liquid phase is present in bubbles and dispersed bubble flow starts to decrease, slugs of gas and liquid flow arise consequently. The large Taylor gas bubbles having a shape of a bullet contain most of the gas phase and occupy the whole cross-sectional area. They are surrounded by a liquid film, falling down at the walls. Liquid slug also occupy the whole cross- sectional area and separate Taylor bubbles one from another. They contain small dispersed bubbles inside. A different approach is required to calculate the mixture density. For pressure drop calculations a unit of a Taylor bubble and a water slug must be considered.

When gas rate is increased, Taylor bubbles become larger in size and liquid contains more dispersed gas inside. If the critical value ( which is about εg = 0.52) is reached, the liquid slug stars to break and transition or churn flow is observed.

Transition (Churn) Row

At higher gas rates neither liquid slugs nor Taylor bubbles exist any longer, neither of the phases is continuous. Liquid now is lifted by the smaller Taylor bubbles of distorted shape. Hence, the liquid phase makes an up and down motion in vertically alternating directions.

When the gas rates increase even more, churn flow is transformed into mist or annular pattern.

Liquid film stability and bridging are two mechanisms that govern the transition.

(16)

Annular (Mist) Row

At extremely high gas velocities, gas phase becomes continuous along the pipe. Liquid is presented as a wavy film on the pipe walls and as small droplets in the gas flow. The last ones are transported at the same velocity as the gas phase velocity. No slippage occurs between the phases.

Density, therefore, can be calculated from the non-slip holdup. Determination of film thickness on the pipe walls is of high importance because of friction losses occurring on the gas-liquid interface (Barnea, 1987).

Figure 1- Flow regimes and flow regime map in vertical two-phase flow. (Time, 2009)

(17)

2.5 Nodal analysis principle

“Nodal analysis is an approach for applying systems analysis to the complete well system from the outer boundary of the reservoir to the sand face, across the perforations and completion section, up the tubing string, the flow line and separator”. A producing well can be considered as a series of hydraulic connected systems that are bordered by certain points, called nodes (in connection with Figure 2). The node can be classified as a functional node in case there is a pressure differential across it and the pressure or flow rate can be represented by a physical or mathematical function (Stoisits, 1992).

System’s performance evaluation can be done by meeting the following requirements:

• Mass flow rate throughout the system is constant, although with changes in pressure and temperature phase conditions change too.

• Pressure decreases in the direction of flow because of the energy losses occurring in the various system components.

• At node points, input pressure to the next component must equal the output pressure of the previous component.

• System parameters being constant for considerable periods of time are:

• the endpoint pressures at the separator and in the reservoir

• the wellbore and surface geometry data (pipe diameters, lengths, etc.)

• the composition of the fluid entering the well bottom

Flow rate can be determined, taking into account above mentioned features of the production system. It can be described by a procedure in which the system is divided into two subsystems at a node that is called the solution node. Starting from the points with known constant pressures (well boundaries-bottomhole and separator pressures), local pressures in the different hydraulic elements of the production system are calculated.

The system analysis principles can be applied to the gas lift well analysis. This is made possible by the fact that multiphase flow is present in the tubing string. A gas lift well scheme is given in Figure 2 (Takacs, 2005).

(18)

Figure 2- The production system of a continuous flow gas lift well (reproduction) (Takacs, 2005)

2.6 Gas lift installation type

In the considered case in the following master thesis an assumption has been made that gas injected in the gas lift well comes to the lower part of the tubing string. It means that gas lift installation type is open.

In the open gas lift installation the production oil system consists of the tubing, hanging inside the casing string without packer installed (Figure 3).

(19)

Figure 3- Open tubing flow installation with gas injection at the tubing shoe (Takacs, 2005) This type was used in the begging of the gas lift era, when no gas lift valves were installed and gas was injected at the tubing shoe. It leads to enormous gas consumption and hence low effectivity. The proper setting depth selection is a difficult task. It is a function of the surface injection pressure. Any changes in the surface pressure after the system is run may lead to wide range of fluctuations in the injected gas rate. Liquid rate fluctuates accordingly. The requirement of continuous flow of gas at a constant gas injection volume cannot be achieved (Takacs, 2005).

2.7 Pressure gradients calculation

In multi-phase flow pressure drop is probably the quantity that researchers deal with most often. Large number of dimensionless variables is a big problem. For example, the friction factor is identified as a function of the Reynolds number, Froude number, the Weber number, the flowrate ratio, density ratio, and the viscosity ratio (Griffith, 1984).

One may consider the total pressure gradient in the pipe as a function of 3 different terms:

(20)

Frictional pressure gradient

f

dp dx

 

 

  , hydrostatic pressure gradient

h

dp dx

 

 

  and acceleration pressure gradient

a

dp dx

 

 

  . Thus

f h a

dp dp dp dp

dx dx dx dx

       

       

        (2.1) Each of these terms contribute in a different way in single phase and two phase flow.

For two-phase flow calculations one starts with the assumption that the fluids mixture properties are given by the mixing rules. This approach is called the homogeneous two-phase pressure drop model. In order to obtain appropriate fluid fractions calculation may be done by assuming either no-slip (S = 1) or by specifying a certain slip ratio. The Beggs and Brill pressure drop model applies a homogeneous model for friction and hydrostatic pressure, some parts calculated with the no-slip assumption, and some taking into account the real slip.

Frictional pressure drop

As in single phase, the following expression is used:

 

²

4 1

Re 2

n

m m mix

f

dp C U

dx D ^ 

    

 

  (2.2)

Where index “m” means “mixture”. The homogeneous model is realistic essentially in turbulent well mixed flow, which means the selection C=0.046, n=0.2 should be used.

The mixture Reynolds number is calculated as:

Re_m ^m ^mix

m

 U D

  (2.3)

Hydrostatic pressure gradient

The hydrostatic pressure gradient can be calculated as follows:

m cos

h

dp g

dx  

  

 

  (2.4) (Angle  relative to vertical direction)

Where the mixture density is different from the single phase flow.

Acceleration pressure drop

The two most important contributions come from:

• Change in gas density

(21)

• Change in velocity change and pipe cross sectional area The acceleration pressure gradient is similar to the single phase flow:

mix m mix a

dU

dp U

dx  dx

    

 

  (2.5) modifying only the density and using the mixture velocity (Time, 2009).

(22)

Simulation part Gazizullin E. 22

3 Simulation part

3.1 Previous works review

A number of papers and articles on experimental studies were researched in order to obtain a theoretical basis for planned simulation work.

A great work has been done by the authors of (Vazquez & Hernandez, 2005), where a continuous gas lift model has been developed based on nodal analysis. Three cases of study were considered to check if the developed model resembled experimental data. Problems solved included gas lift optimization, namely determination of appropriate amount of gas required to lift reservoir fluid, provided that outlet pressure at surface manifold matched actual pressure.

Researchers used computational approach and concluded, that proposed model had produces appropriate results and could be applied in simulation of optimization of gas lift systems.

Paper (Chia & Hussain, 1999) by Y. C. Chia and Sies Hussain gives an insight of challenges that are encountered during gas lift optimization in Malaysian oil fields. It is mentioned that to meet the challenges, the GOAL model (Gas Lift Optimization Allocation Model) is applied, which is a developed PC based production system model. It is also based on nodal analysis in order to generate well performance curves for a system of wells. Optimization is also possible by defining an objective (maximizing liquid production, oil etc.) Model considers emulsion and sand production and is able to simulate dual completion gas lift design.

Paper (Wang & Litvak, 2004) describes a method for gas lift optimization problem for multiple objects and applies certain algorithms that reduce calculation time because of less CPU is required, which is vital in real oilfields, where the amount of wells is significant. The method is named “GLINC”, and generates results for long-term simulation studies.

Similarly to the present work, a work has been done by the authors of paper (Bahadori, Ayatollahi & Shirazz, 2001) for the Ayatollahi oilfield, where it is mentioned that the choice of appropriate PVT and fluid correlations play an important role. After the author have selected the correlations, performance curves was calculated. An optimum production rate was then found from

(23)

the well performance curve and finally valve installation depth was designed. All calculations were performed in a numerical simulator, developed by the authors.

Paper (Pablano, Camacho & Fairuzov, 2005) develops a concept of stability maps for continuous gas lift flow, when transition from stable to unstable flow occurs. It allows to compare different gas lift stability criteria (injection rate, the injector port size, tubing diameter) and quantify the effects of these criteria. The authors mention that well-stabilizing methods may be developed on the basis of the concepts, that provides minimum CAPEX and OPEX.

3.2 Introduction to the program

The software package developed by the author [myself] in the following thesis consists of separate calculation blocks- calculation of PVT properties section, the calculation of pressure gradients section and gaslift parameters calculation section. In its turn, each section consists of separate functions. Some of the functions refer to other ones, which are more elementary. Thus, every function that is used in the program will be described. The description of every function is comprised of the following elements:

1. “Name”-name of the function;

2. “Variables”- variables of the function;

3. “Subfunctions”- called functions that are used in the body of the function;

4. “Function’s body”- a set of formulae, that is used to obtain the results of the function;

5. “Calculated parameters”- function’s output.

Here is a summary flowchart of the program:

(24)

Figure 4- Summary flowchart of the program

Methodology of a gaslift well optimization consists of execution of the following basic steps:

- Preparation of input data describing the necessary parameters of the drained reservoir, production casing of the well, steady-state flow, as well as the properties and fractions of produced oil, water and gas at different thermodynamic conditions;

- Calculation and making a plot of pressure distribution curves (PDC) along the length of the production casing and tubing in the interval from the well’s bottomhole to the surface for a given production oil rate at standard conditions.

(25)

- Calculation of oil and liquid production depending on amount of gas injected in a gaslift well

- Optimization of gaslift, calculation of optimal gas injection rate.

Input data that is necessary and sufficient if shown below:

Table 1- Well input data

Parameter Value Units

Reservoir pressure 17,4 MPa

Productivity 22 m3/(d*MPa)

Reservoir temperature 103 C Geotharmal gradient 0,02 deg/m Wellbore flowline

pressure 0,84 MPa

Annulus pressure 1 MPa

Wellbore length 1607 m

Inclination 0,2 deg

Annulus diameter 157 mm

Tubing duameter 68 mm

Oil density at SC 819,9 kg/m3 Water density at SC 1120 kg/m3 Associated gas density 1,333 kg/m3 Gas injection rate 0,7 m3/s

Watercut 60 %

Saturation pressure 8,5 MPa Gas injection rates range 0-2 m3/s

3.3 Calculation steps

At first, a pressure distribution curve (PDC) along the casing has to be built. It is performed stepwise. The well length is divided into n=100 segments upwards from bottom to top. PDC is an array of pressures and positions in the well, so the goal is to determine local pressures along the well, namely at the ends of each segment. However, it is more convenient to determine pressure gradients instead.

In general, pressure gradient is defined from Bernoulli’s equation:

2 2 2

(

.

) cos

2

g g

o o w w

o o w w sc g g

o w g

P g

L D

     

       

  

 

                 

(3.1)

(26)

Where:

, ,

o w g

   is oil, water and gas fractions respectively;

, ,

o w g

   is oil, water and gas densities respectively;

, ,

o w g

   is oil, water and gas flowrates respectively;

 is average inclination;

D is annulus (tubing) diameter;

- hydraulic friction coefficient.

Since liquid and gas conditions are determined not only at a certain pressure but at a temperature as well, it is important to obtain the temperature gradient:

.

20 2,67

(0, 0034 0, 79 ) 10

Ql sc

D

T G

L

  



(3.2)

Where;

.

Ql sc - well inflow at SC;

G- Geothermal gradient.

The starting point is bottomhole of the well. However, bottomhole pressure has not been found yet. Here an iteration procedure has to be applied. The idea of the procedure is that pressure at the bottomhole is assumed to be known. Knowing the reservoir pressure, well flow can now be determined. Then further calculations are carried out (they are described below in more details), the final result of which is calculated (not real) wellhead pressure. Comparing this pressure and the factual one, assumed bottomhole pressure is “calibrated” until the real wellhead pressure is equal to the calculated wellhead pressure, or more precisely the difference between them is less than a certain tolerance value. When the equality has been reached, such parameters as the well flow and the bottomhole pressure now become determined.

When bottomhole pressure is determined, the three fundamental values which are bottomhole pressure, temperature and position are established. Using empirical correlations and theoretical calculations presented in the Matlab program it is easy to obtain fluid properties at bottomhole conditions that in its turn makes it simple to obtain local pressure and temperature gradients, according to the formulas above.

Pressure at the second end of the first interval can be calculated as follows:

out in

P P P L

L

   



(3.3)

(27)

Where:

P

in-pressure at the inlet of I stage;

P

out- pressure on the next stage.

So, the pressure and the temperature at the both ends of the first interval are know now.

Pressure at the end of the first section is equal to the one at the beginning of the second section.

Hence the same operation is proceeded for the whole range of length intervals, giving us the pressure-temperature-local position array.

The same procedure is performed to obtain an array of pressures for tubing space. Calculations are done upwards from the lower position of the tubing to the wellhead.

It is worth mentioning that during the calculation one obtains arrays for the oil mixture properties, as local densities, volumetric fractions, oil viscosities and so one. By knowing critical values that are not to be exceeded, one may analyze and optimize the oil flow.

After all, pressure distribution curves for annulus and tubing for a given well are obtained.

To get a better understanding of the program the following figure is given:

Figure 5- Structure of the program

In this figure the structure of the program is shown. The two main block calculate pressure distribution curves (PDC) for casing and tubing space for given well parameters. However, at that step the factual well flow rate at bottomhole conditions is unknown. The iteration process then is required, which calculates a series of PDC for the same well with changing well inflow, unless the

(28)

wellhead pressure (which is the resulting parameter of PDC) coincides with the real wellhead pressure given as the initial datum.

Then, when gas lift is simulated, a range of gas injection rates is tested, which requires an additional loop, containing all previous loops. Then based on results of the previous loop, the production curve is built (it represents liquid production and gas lift operating pressure vs. gas injection rate), from where the optimum injection range of rates at which maximum liquid production can be reached is found.

More precisely the program description is given in the next subchapter. For the program code, see Appendix A. There are comments, describing steps that may be challenging, so it is highly recommended to look through the code.

(29)

3.4. Program description

3.4.1 Calculation of PVT

For calculation of physical and chemical properties of the reservoir fluids a technique, suggested by (Lyapkov, 1987), (Dunushkin & Mishenko, 1982), (Mishenko, 2003) is used.

3.4.1.1 Calculation of gas saturation

Function- gas_saturation

Here and after the correlations, that are used in the PVT properties calculations, dynamic parameters calculations and other calculations are based on the following handbooks:

(Gimatudinov, 1983), (Gimatudinov, 1974), (Buhalenko et al., 1983), (Trebin et al., 1980), (Shtof, 1974).

Table 2- Input data

Parameter Symbol Unit

Pressure P ^MPa

Saturation pressure

P

sat ^MPa

Gas to oil ratio at saturation pressure

GOR

sat m3/m3 GOR coeff.-s

m

GOR ^-

n

GOR ^-

, ,

nGOR

GOR sat

sat sat

m P P P

GOR GOR P P

  

   (3.4)

(30)

Figure 6- Gas saturation versus pressure at standard pressure and reservoir temperature for Karakuduk oilfield (reservoir J-1/2)

Table 3- Output data

Gas to oil ratio

GOR

m3/m3

3.4.1.2 Calculation of oil volume factor

Function- volume_factor

Pressure P ^MPa

Saturation pressure

P

_sat MPa

Volume factor coeff.-s

m

b -

n

b - ,

,

b

n

b sat

o n

b sat sat

m P P P b

m P P P

  

 

 

 (3.5) As can be seen from Figure 7:

0 20 40 60 80 100

0 5 10 15

GOR, m3/m3

Pressure P, MPa

Gas saturation versus pressure at standard pressure and reservoir temperature

(31)

Figure 7– Oil volume factor versus pressure at reservoir temperature for Karakuduk oilfield (reservoir J-1/2)

Oil volume factor bo m3/m3

3.4.1.3 Calculation of oil density

Function- oil_density

Pressure P ^MPa

Saturation pressure

P

_sat MPa

Oil density coeff.-s

m_ - n_ -

/ ,

n

sat

o n

sat sat

m P P P



 _{ }^^ ^

  (3.6) As can be seen from Figure 8:

1,1 1,15 1,2 1,25 1,3 1,35

0 5 10 15

Oil volume factor, m3/m3

Pressure P, MPa Oil volume factor versus pressure at reservoir

temperature for

(32)

Figure 8–Oil density versus pressure at reservoir temperature for Karakuduk oilfield (reservoir J-1/2)

Oil density o kg/m3

3.4.1.4 Calculation of oil viscosity

Function- oil_viscosity

Pressure P ^MPa

Temperature Т ^K

Saturation pressure

P

sat ^MPa

Oil viscosity coeff.-s m_ ^-

n_ ^-

Oil viscosity is related to pressure by the following equation:

/ ,

n

sat

o n

sat sat

m P P P m P P P



_{ }^^ ^

  (3.7)

700 710 720 730 740 750 760 770

0 5 10 15

Oil density, kg/m3

Pressure P, MPa

Oil density versus pressure at reservoir temperature

(33)

As can be seen from Figure 9:

Figure 9– Oil viscosity versus pressure at reservoir temperature for Karakuduk oilfield (reservoir J-1/2)

However it is known that oil viscosity is strongly dependent on the temperature. Lewis- Squires correction, which considers it, needs to be introduced:

0 0,002 0,004 0,006 0,0080,01 0,012 0,014 0,016 0,018

0 5 10 15

Oil viscosity, Pa*sec

Pressure P, MPa Oil viscosity versus pressure at reservoir

temperature

(34)

Figure 10- Liquid viscosity change due to temperature according to Lewis- Squires

Example:



_o₁

 2, 4

MPa*s;

T

₁

 313

K; T₂ 293 K. Determine __o₂.

Solution : Lay off the point _o₁2, 4MPa*s on the ordinate. From this point we go right up to the curve; shift along the curve by  T T₁T₂31329320 K. We see what new value of the viscosity corresponds to this shift. _o₂ 6, 3MPa*s.

Thus, when calculating the viscosity of oil as opposed to gas content, density and volume faction calculations, temperature changes have to be taken into account.

0,1 1 10 100 1000 10000 100000

V is c o s it y, mP a *s

Temperature change, K

Liquid viscosity change due to temperature according to Lewis- Squires

100 K 100 K 100 K 100 K 100 K 100 K

(35)

Oil viscosity o MPa*s

3.4.1.5 Calculation of coefficients for oil average density, viscosity, volume factor and gas saturation calculations

If coefficients for calculations of oil density, oil viscosity, GOR and volume factor are unknown, then corresponding coefficients are determined by solution of the following se of equations:

    ^{ }

line y y line

sat y y sat

lg y lg m n lg P lg y lg m n lg P

  

   (3.8)

Where

y

_line and

y

_sat -values of a considered function (



_o,



_o,b_o,GOR), taken from the corresponding curve at line pressure

P

_line and saturation pressure P_sat.

As the result, output data is m_,n_,m_,n_,m n m_b, _b, _GOR,n_GOR coefficients.

3.4.1.6 Compressibility factor calculation

Function- compressibility_factor

Pressure P ^MPa

Temperature Т ^K

Saturation pressure

P

sat ^MPa

Assosiated gas density



gSC ^kg/m3

Nitrogen volume fraction in the associated gas at SC

N2

y ^m3/m3

We calculate gas compressibility factor having previously calculated:

Relative gas density:

(36)

. 1, 205

gSC g rel

   (3.9)

Relative density of all HC and non-HC components except nitrogen:

2

. .

( 0,97 )

1

g rel N

HC rel

N

y y







^

 (3.10)

Pseudoreduced pressure and temperature:

.

2

10 46,9 2,06

HC rel

pr

P P

 



(3.11)

97 172 . pr

HC rel

T T

 

 (3.12) If the value of T_pr is less than 1,05, assume _T_pr =1,05

Compressibility factor excluding nitrogen is obtained from the following formula:

2

3,45 6,1

1 0, 23 (1,88 1, 67 ) , 0 1, 45 &1, 05 1,17 0,13 (6, 05 6, 25) ,1, 45 4 &1, 05 1,17

1 ( 0,18 0,135) 0, 0161 , 0 4 &1,17 2

0, 73

HC pr pr pr pr pr

pr

HC pr pr pr pr

pr

HC pr pr pr

pr pr

z P T P P T

z P T T P T

P

z P P P T

T T

 

        

 

       

 

        

  

(3.13) z value is obtained from:

(3.14) Table 11- Output data

Compressibility factor _z -

2 2

(1 )

HC N HC N

z  z  y  z y

(37)

3.4.1.7 Calculation of gas density

Function– gas_density

Pressure P ^MPa

Temperature Т ^K

Assosiated gas density



gSC ^kg/m3

Compressibility factor z -

Table 13- Called functions Name of the function Denotation

Compressibility factor calculation compressibility_factor

Gas density at a certain cross-section area we may get from:

.

SC

g g rel

SC

PT

  zP T (3.15)

where P_SC- Pressure of air at SC, 101,3 kPa

Tsc- Temperature of air at SC, 293,2 K

Gas density



g ^kg/m3

(38)

3.4.1.8 Calculation of fluid’s volumetric water fraction

Function– volume_share

Pressure P ^MPa

Water volume fraction in produced fluid at SC

.



w sc ^m3/m3

Oil volume factor bo m3/m3

Table 16- Called functions

Name of the function Denotation

Calculation of oil volume factor volume_factor

Volumetric water content in liquid is calculated from:



.



1

1 1

wl

o w sc

 b

 

  (3.16)

Volumetric oil content in liquid is calculated from:

o

1

wl

   

(3.17)

Volumetric water content in produced fluid at given temperature and pressure conditions



wl ^m3/m3

(39)

3.4.1.9 Calculation of the surface tension between the liquid phase that is the outer phase of the stream and gas



lg

Function– surface_tension

Zone (annular space or tubing) - -

Pressure P ^MPa

Temperature Т ^K

Volumetric water content in produced fluid at given temperature

and pressure conditions



wl _m3/m3

Superficial velocity of the mixture at given temperature and

pressure conditions



mix _m/s

First critical velocity of the mixture at given temperature and

pressure conditions



cr1 _m/s

Calculation of fluid’s volumetric water fraction volume_share

Calculation of oil, water and gas flow rates and first and second critical rates flow_rate

Zone (annular space or tubing) parameter is introduced. It can have 2 values: “Casing”- when calculating mixture parameters in the annular space and “Tubing”- when calculating in tubing area.

We calculate the surface tension value at water-gas contact:

1,19 0,01

10

^P



wg



^ ^ (3.18) Then we calculate the surface tension value at oil-gas contact:

(40)

 

1,58 0,05 6

10

^P

72 10 305

og

T

 

^ ^

 

^



(3.19)

If the value



og turned out to be below zero, then:

og

0  

(3.20)

We calculate the surface tension value at oil-water contact:

ow wg og

    

(3.21)

Eventually, the value of the surface tension at liquid-gas contact



_lg we obtain from the following scheme:

Figure 11– Choice of the surface tension between the liquid phase as the outer phase of the stream and the gas

Surface tension at oil-water contact



_ow N/m

Surface tension at liquid-gas contact



_lg N/m

(41)

Figure 12– Surface tension vs watercut for Karakuduk oilfield (reservoir J-1/2)

3.4.1.10 Calculation of fluid viscosity

Function–fluid_viscosity

Zone (annular space or tubing) - -

Pressure P ^MPa

Temperature Т ^K

Water density at SC

.



w sc ^kg/m3

Annulus (tubing) diameter D ^m

Water volume fraction in liquid



wl ^m3/m3

Superficial velocity of the mixture at given temperature and pressure conditions



mix ^m/s

0 0,01 0,02 0,03 0,04 0,05 0,06

0 0,2 0,4 0,6 0,8 1

Su rface ten si o n at li q u id -g as co n tact , N /m

Watercut

Surface tension vs watercut

(42)

Table 22- Continuation First critical velocity of the mixture at given temperature and

rpressure conditions



^cr¹ m/s

Second critical velocity of the mixture at given temperature and

pressure conditions



^cr² m/s

Oil viscosity at given temperature and pressure conditions



_o MPa*s

Calculation of fluid’s volumetric water fraction volume_share

Calculation of oil, water and gas flow rates and first and second critical rates flow_rate

Calculation of oil density oil_density

We calculate water viscosity:

7 . 0,0065( 273)

0, 0014 38 10 ( 1000) 10

w sc

w T



 ^ ^ ^ _



^ (3.22)

We calculate the value of A parameter:

2 0,48

(1 20 )

8

^wl

wl mix

A

D





 

 

 

 

(3.23)

Eventually, the value of liquid viscosity



_l we get from the scheme:

(43)

Figure 13- Choice of liquid viscosity

Water viscosity



w ^MPa*s

Liquid viscosity



l ^MPa*s

Figure 14- Water viscosity versus watercut for Karakuduk oilfield (reservoir J-1/2) (Gazizullin, 2014)

Gas lift simulation and experiments in conjunction with the Lyapkov P. D. methodic

MASTER’S THESIS

Acknowledgements

Abstract

Contents

List of tables

List of figures

1 Introduction

Project work

Motivation for oil production

2 Theory

2.1 Gas Lift Operational Principle

2.2 Multiphase flow

2.3 Basic principles

2.4 Flow regimes

2.5 Nodal analysis principle

2.6 Gas lift installation type

2.7 Pressure gradients calculation

 

3 Simulation part

3.1 Previous works review

3.2 Introduction to the program

3.3 Calculation steps

(

) cos

2

P g

L D

     

       

  

 

                 

(0, 0034 0, 79 ) 10

T G

L

  



P P P L

L

   



P

P

3.4. Program description

3.4.1 Calculation of PVT

3.4.1.1 Calculation of gas saturation

P

GOR

m

n

GOR

3.4.1.2 Calculation of oil volume factor

P

m

n

3.4.1.3 Calculation of oil density

P

3.4.1.4 Calculation of oil viscosity

P





 2, 4

T

 313

0,1 1 10 100 1000 10000 100000

V is c o s it y, mP a *s

Temperature change, K

Liquid viscosity change due to temperature according to Lewis- Squires

3.4.1.5 Calculation of coefficients for oil average density, viscosity, volume factor and gas saturation calculations

     

     

y

y





P

3.4.1.6 Compressibility factor calculation

P



    ^{ }

    ^{ }