# The Framework Dedicated to Three Phase Flows Wellbore Modelling

Authors: Bartlomiej Bielecki; Andrzej Krajka
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

To predict physical properties in a wellbore during oil and gas production, scientists use empirical correlations or mechanistic approach algorithms. The typical research in this field is concentrated on a single property analysis as heat transfer, pressure, temperature, and so forth. Here the most proper correlations, regarding the subject, are presented. And the authors studied how to join all correlations into the full framework which returns all production parameters at every depth in a wellbore. Additionally, the presented simulation results are studied here. Based on presented algorithms, the proper tool has been applied and the results shown in this paper are taken from this application.


---

## Body

## 1. Introduction

All of well physical properties in wellbores are deeply studied with many approaches found in literature. Oil and gas industry is a cutting edge domain with a lot of investment involved. This is the reason why it is hard to find any description of a framework which join correlations altogether into one solution. The solution can describe the wellbore physical properties, at every depth, at every stage of production, even for the single whole oilfield. It means that all calculations have been performed every foot across the wellbore, using standard oilfield unit system. Once a few wellbores are considered, the results object keeps all calculations, still performed for every foot, for all wells independently. The final results set is dedicated for oilfield engineers, so pressure, temperature, and production results are most desirable. The whole model considered as full wellbore simulation was presented as a Ph.D. thesis at The University of Texas [1]. The global simulation model is presented there, but without any description of the framework. Both models are based on the same principles but with different approaches. Our aim is to show the possibility of joining correlations altogether using fast and efficient algorithm. Additionally, our model is extended for gas lift which is a part of this framework and can be used in wireless sensor network well monitoring [2]. Optimisation of gas lift is a well known problem [3, 4], but we are concentrating on the gas injection influence on the final results of simulation. Unfortunately, the authors do not have any possibility to verify the simulation results by the comparison to the real data. This kind of data is protected by the oil companies.Correlations presented in Section2 are typical wellbore simulation models. Three types of correlations described in this section are identified: empirical correlations, mechanistic approach, and homogeneous models. Typical approaches for the vertical wells are presented. Three phase flows are considered. The pressure model was developed based on [5] study with the consideration of flow regimes and their dependence on the features of empirical correlations. The temperature model is based on Sagar et al. [6] model.Section3 presents the framework with the pseudocode and some chart to understand the complexity of wellbore modelling. This is the model with the diversion on reservoir and wellbore calculation. Reservoir calculations use layer pressure and geothermal temperature. Wellbore correlations use the pressure and temperature models, as the iterations where the pressure values from any previous calculation are passed to the calculation object at current depth. This model was developed based on well known principles [5]. Its main advantage is a speed and effectiveness.In Section4, typical case study can be found. The tool which has been developed to present all data at every stage of simulation has been implemented. We can present variables as production rate, pressure, temperature, density, solubility, compressibility, formation volume factor, viscosity, heat capacity, surface tension, gas/oil ratio, gas/water ratio, mass rate, bubble point pressure, flow direction, flow regimes, velocity, hydrostatic and kinetic pressure drop, Reynolds number, and others as function of mentioned. It gives the possibility to present different kind of data, which are comparable to other researches.To summarize, the main contributions of this paper are to show the framework idea, which groups wellbore correlations altogether. This framework gives the reliable simulation results which is proved in this paper.


## 2. Correlations

The authors would like to present the single correlation for every property. There is no clear algorithm that suggests which correlation should be taken under the consideration as the most proper one. It depends on many factors as fluid density, trajectory, and so forth. The best way to gain the accurate modelling results is to compare obtained data with the real one. Good example can be found from Tulsa oilfields analysis [1] where it has been proven by the comparison of relative performance factor that the mechanistic approach gives better results than empirical correlations. Ansari model presented in Figure 1 is the best in this example. It does not necessarily means that the Ansari model is reliable for every wellbore simulation. It is only regarding this case.Figure 1
The different models comparison [1]. The lower value means the better results.Here the most important part of calculation is presented. Three major parameters as pressure, temperature, and inflow exist as a function of the following variables: bubble point pressure, gas solubility, gas compressibility, formation volume factors, viscosities, heat capacities, surface tension, completion data, and reservoir properties.


### 2.1. Pressure

Our method is based on solving the equation using a finite difference approximation method with the continuity, momentum, and energy equation as follows:(1)ϱil+1-ϱilΔt+ϱvi+1/2l+1-ϱvi-1/2l+1Δx=0,ϱvi+1/2l+1-ϱi+1/2lΔt+ϱvi+1/22-ϱvi-1/22lΔx+pi+1l+1-pil+1Δx+ϱi+1/2lg+fϱvi+1/222Dl=0,ϱiEi+1l+1-Eil+1Δt+ϱviEil-Ei-1lΔx+Pil+1vi+1/2-vi-1/2lΔx+2HTr-Twlrw=0.Across the years, the multiphase flow theory has been improved and makes the calculation difficult to predict due to several reasons. In the flow problem, the number of variables affecting pressure drop is enormous. The interfacial tension between phases is taken into account, but all kinds of available models try to reduce the total number of variables with the introduction of nondimensional parameter group. Frictional pressure losses are more difficult considering multiphase flow, so, in addition, new kind of energy losses has been used: slippage loss. All these features make the model so complicated that we cannot describe it in this paper. It is important to emphasize that usually the correlations which are used are based on empirical method.Pressure loss across the wellbore during the production is divided into three components:(2)dpdltotal=dpdlhydrostatic+dpdlkinetic+dpdlfrictional.The cohesive forces (CF) and frictional forces (FF) have the crucial influence on fluid drive, especially regarding phases separation or bind. The friccohesity model has been considered as a CF and FF product. Using the Mansingh equation [7](3)σ=σ0tt0±Btnn0±0.00121-ϱ,the friccohesity as a function of surface tension and cohesive and frictional forces is calculated. In this equation, t0 and t are the reference and sample flow times, respectively, r0 is reference friccohesity, and n0 and n are the pendant drop number of reference and sample, respectively. Reference friccohesity is represented by the following formula:(4)σ0=μ0ST0.Mansingh survisemeter is a device to measuring friccohesity [8]. The friccohesity analysis is dependent on hydrocarbon composition. In this model, we are concentrating on the compressibility factor calculation with pseudo-reduced pressure and temperature solutions, but the composition role is replaced by the density value. It does not mean that the model is not predicted for hydrocarbon composition, but, from the optimization point of view, our solution is faster upon this simplification. The deeper friccohesity analysis is consider as a future work for this study. That is the reason why the Papay, as the best compressibility correlation in the comparison to the others, has been chosen [9]. There are solutions which do not consider kinetic pressure losses [10]. Inside the correlations which determine the flow regime, liquid holdup for more than single phase, the pressure drop equation is known under the different form presented as follows:(5)dpdltotal=dp/dlhydrostatic+dp/dlfrictional1-Ek.Depending on the correlation, this kinetic term has the different form, but usually this is a variation of the equation:(6)Ek=0.000216fvmvsgϱnsp.Velocity of each phase is calculated from(7)V=576ΠD2∑i=1n‍BiQi.Oil formation volume factor algorithms are calculated using three different methods: correlation based [5], artificial neural networks based prediction [11], and hybrid soft computing based techniques [12]. The crucial issue of gas formation volume factor correlation is compressibility. Last study estimated this value from the experimental data taken as a function of pressure and temperature [13]. Water formation volume factor is already described as function of two polynomials: pressure and temperature dependent [14]. Consider(8)Bw=1+ΔVwP1+ΔVwT,ΔVwP=-1.0001×10-2+1.33391×10-4T+5.50654×10-7T2,ΔVwT=-1.95×10-9pT-1.728×10-13p2T-3.589×10-7p-2.253×10-10p2.In (6), friction factor is changing due to flow regimes as a function of Reynolds number, diameter, and pressure as follows [15]:(9)f=fL1-ψb1/31-ψb13+fTψb1/3,where ψ is intermittency factor:(10)ψb=log⁡Re/ReLlog⁡ReT/ReL.Reynolds number is calculated from the following formula: (11)Re=124*ϱmVmDμm.Hence,(12)log⁡ReL300=1.7pD-1,log⁡ReT105=0.7pD-1.In the whole calculation, density and viscosity are the most important flow parameters. Ideally, they are determined experimentally in the laboratory on actual fluid samples taken from the field under study. In many cases, correlations are region dependent, so every simulation should provide options to choose a place of production. Dead, saturated, and undersaturated are considered as oil types and identified using fluid properties. Every type has its own correlations set [16].The last parameter in (6) is no slip density which is described in the following formula [5]: (13)ϱns=VlϱlVm+1-VlVmϱg.Hydrostatic pressure is calculated from(14)dpdlhydrostatic=1144ϱm(T,α,Fr,ϕ,St)dxdl.As it is shown, basically ρm is a function of temperature, inclination, and flow regime. Every correlation has the deep study of mixture density. Pressure drop caused by friction has been studied widely. Usually the variations of the following equation are considered [5]: (15)dpdlf=1.294×10-4fϱnsϱlVl2D.A flow in a pipe is turbulent if the Reynolds number is greater than 4000 and laminar below 2100. For laminar flow, the friction factor is calculated by assumption as 64/Re, but for the transition and turbulent flow Chen correlation is an example of friction factor calculation [17]. Consider(16)f=4log⁡R3.7065D-5.0452log⁡(δ)Re2-1,where (17)δ=12.8257RD1.1098+7.149Re0.8981.The equations which have been presented in this section are examples of pressure loss correlation. Flow regimes and pipe trajectory are the crucial factors in every calculation. A lot of exceptions have been studied and analysed and during the simulation process.In the pressure calculation, mixture density depends on surface tension which can be described as a function of pressure and temperature  [18]. As an example, the presented correlation has been created for crude/oil gas systems [19] (18)ST=γ1-t-74γ1-γ2206,where (19)γ1=75-1.108p0.349,γ2=53-0.1048p0.637.Pressure calculation is very complicated in this model and should be presented in another paper. It is worth to say that parameters as superficial velocities, Froude numbers, volume fractions in different flow patterns, and pipe properties having the influence on frictional forces are considered. We would like to focus on this part in next publication.Pressure correlations are good example of diversity between different studies and approaches. In our tool, four different pressure correlations return different results in the bottom hole, respectively [5], Beggs-Brill 1620 psi, Orkiszewski 1580 psi, Aziz-Govier-Fogarasi 1616 psi, and Duns-Ross 1507 psi.


### 2.2. Temperature

From the simulation point of view, temperature is a parameter of pressure, but from the engineering point it is one of the most important factors regarding the production description and understanding. In the calculation, the efficient correlation is the iteration based on the previous values [6] (20)Tfz=GTz-AsinαJCm+AF+Aggsinα+exp⁡E,where (21)A=ΠDUκMMCκ+D2Clog⁡48TDtD-0.29-1,E=-z-z0AAsinαJCmiii·Tfz0-GTz0+AsinαJCm-AF-Aggsinα.F is a correction factor combined with hydrocarbon expansion from high pressure to low pressure during the temperature change, known as a Joule-Thomson effect. Heat transfer occurs between the wellbore fluid and the formation, overcoming resistances offered by the tubing wall, tubing-casing annulus, casing wall, and cement [19].The mixture heat capacity calculated depends on how many phases are considered. Consider(22)Cm=∑i=1n‍MiCi∑i=1n‍Mi.Heat capacity of each phase is a function of temperature and density usually estimated from the experimental data. Gambill correlations are valid, respectively, for oil and gas as it is shown [20]: (23)C=ϱ-1/20.338+0.00045T-460.


### 2.3. Inflow

The derivation of Darcy’s law is used in inflow calculation to determine the flow through the permeable geological layer. This equation is valid for liquid phase inflow at the perforation points, but bubble point pressure has not been taken under the consideration [5]: (24)Q=0.00708khpr-pwμBlog⁡re/rw-0.75+S.Drainage radius may be multiplied by some coefficients, depends on reservoir shape. Skin factor is estimated based on well test analysis. A positive value indicates there is a pressure decline in the near vicinity of well that is more than expected based on the radial flow equation but does not necessarily indicate formation damage [21]. Considering gas inflow, Darcy equation is adjusted in terms of gas properties: (25)Q=0.00708khPr2-Pw2μβTlog⁡re/rw-0.75+S.Inflow equations do not account for the phase change in solution gas reservoirs. In a solution gas drive, there is expansion of the hydrocarbons below bubble point, which is beneficial because it adds energy to the system. Gas liberation in a solution gas drive is also detrimental to oil production because it lowers the effective permeability of oil [22].


## 2.1. Pressure

Our method is based on solving the equation using a finite difference approximation method with the continuity, momentum, and energy equation as follows:(1)ϱil+1-ϱilΔt+ϱvi+1/2l+1-ϱvi-1/2l+1Δx=0,ϱvi+1/2l+1-ϱi+1/2lΔt+ϱvi+1/22-ϱvi-1/22lΔx+pi+1l+1-pil+1Δx+ϱi+1/2lg+fϱvi+1/222Dl=0,ϱiEi+1l+1-Eil+1Δt+ϱviEil-Ei-1lΔx+Pil+1vi+1/2-vi-1/2lΔx+2HTr-Twlrw=0.Across the years, the multiphase flow theory has been improved and makes the calculation difficult to predict due to several reasons. In the flow problem, the number of variables affecting pressure drop is enormous. The interfacial tension between phases is taken into account, but all kinds of available models try to reduce the total number of variables with the introduction of nondimensional parameter group. Frictional pressure losses are more difficult considering multiphase flow, so, in addition, new kind of energy losses has been used: slippage loss. All these features make the model so complicated that we cannot describe it in this paper. It is important to emphasize that usually the correlations which are used are based on empirical method.Pressure loss across the wellbore during the production is divided into three components:(2)dpdltotal=dpdlhydrostatic+dpdlkinetic+dpdlfrictional.The cohesive forces (CF) and frictional forces (FF) have the crucial influence on fluid drive, especially regarding phases separation or bind. The friccohesity model has been considered as a CF and FF product. Using the Mansingh equation [7](3)σ=σ0tt0±Btnn0±0.00121-ϱ,the friccohesity as a function of surface tension and cohesive and frictional forces is calculated. In this equation, t0 and t are the reference and sample flow times, respectively, r0 is reference friccohesity, and n0 and n are the pendant drop number of reference and sample, respectively. Reference friccohesity is represented by the following formula:(4)σ0=μ0ST0.Mansingh survisemeter is a device to measuring friccohesity [8]. The friccohesity analysis is dependent on hydrocarbon composition. In this model, we are concentrating on the compressibility factor calculation with pseudo-reduced pressure and temperature solutions, but the composition role is replaced by the density value. It does not mean that the model is not predicted for hydrocarbon composition, but, from the optimization point of view, our solution is faster upon this simplification. The deeper friccohesity analysis is consider as a future work for this study. That is the reason why the Papay, as the best compressibility correlation in the comparison to the others, has been chosen [9]. There are solutions which do not consider kinetic pressure losses [10]. Inside the correlations which determine the flow regime, liquid holdup for more than single phase, the pressure drop equation is known under the different form presented as follows:(5)dpdltotal=dp/dlhydrostatic+dp/dlfrictional1-Ek.Depending on the correlation, this kinetic term has the different form, but usually this is a variation of the equation:(6)Ek=0.000216fvmvsgϱnsp.Velocity of each phase is calculated from(7)V=576ΠD2∑i=1n‍BiQi.Oil formation volume factor algorithms are calculated using three different methods: correlation based [5], artificial neural networks based prediction [11], and hybrid soft computing based techniques [12]. The crucial issue of gas formation volume factor correlation is compressibility. Last study estimated this value from the experimental data taken as a function of pressure and temperature [13]. Water formation volume factor is already described as function of two polynomials: pressure and temperature dependent [14]. Consider(8)Bw=1+ΔVwP1+ΔVwT,ΔVwP=-1.0001×10-2+1.33391×10-4T+5.50654×10-7T2,ΔVwT=-1.95×10-9pT-1.728×10-13p2T-3.589×10-7p-2.253×10-10p2.In (6), friction factor is changing due to flow regimes as a function of Reynolds number, diameter, and pressure as follows [15]:(9)f=fL1-ψb1/31-ψb13+fTψb1/3,where ψ is intermittency factor:(10)ψb=log⁡Re/ReLlog⁡ReT/ReL.Reynolds number is calculated from the following formula: (11)Re=124*ϱmVmDμm.Hence,(12)log⁡ReL300=1.7pD-1,log⁡ReT105=0.7pD-1.In the whole calculation, density and viscosity are the most important flow parameters. Ideally, they are determined experimentally in the laboratory on actual fluid samples taken from the field under study. In many cases, correlations are region dependent, so every simulation should provide options to choose a place of production. Dead, saturated, and undersaturated are considered as oil types and identified using fluid properties. Every type has its own correlations set [16].The last parameter in (6) is no slip density which is described in the following formula [5]: (13)ϱns=VlϱlVm+1-VlVmϱg.Hydrostatic pressure is calculated from(14)dpdlhydrostatic=1144ϱm(T,α,Fr,ϕ,St)dxdl.As it is shown, basically ρm is a function of temperature, inclination, and flow regime. Every correlation has the deep study of mixture density. Pressure drop caused by friction has been studied widely. Usually the variations of the following equation are considered [5]: (15)dpdlf=1.294×10-4fϱnsϱlVl2D.A flow in a pipe is turbulent if the Reynolds number is greater than 4000 and laminar below 2100. For laminar flow, the friction factor is calculated by assumption as 64/Re, but for the transition and turbulent flow Chen correlation is an example of friction factor calculation [17]. Consider(16)f=4log⁡R3.7065D-5.0452log⁡(δ)Re2-1,where (17)δ=12.8257RD1.1098+7.149Re0.8981.The equations which have been presented in this section are examples of pressure loss correlation. Flow regimes and pipe trajectory are the crucial factors in every calculation. A lot of exceptions have been studied and analysed and during the simulation process.In the pressure calculation, mixture density depends on surface tension which can be described as a function of pressure and temperature  [18]. As an example, the presented correlation has been created for crude/oil gas systems [19] (18)ST=γ1-t-74γ1-γ2206,where (19)γ1=75-1.108p0.349,γ2=53-0.1048p0.637.Pressure calculation is very complicated in this model and should be presented in another paper. It is worth to say that parameters as superficial velocities, Froude numbers, volume fractions in different flow patterns, and pipe properties having the influence on frictional forces are considered. We would like to focus on this part in next publication.Pressure correlations are good example of diversity between different studies and approaches. In our tool, four different pressure correlations return different results in the bottom hole, respectively [5], Beggs-Brill 1620 psi, Orkiszewski 1580 psi, Aziz-Govier-Fogarasi 1616 psi, and Duns-Ross 1507 psi.


## 2.2. Temperature

From the simulation point of view, temperature is a parameter of pressure, but from the engineering point it is one of the most important factors regarding the production description and understanding. In the calculation, the efficient correlation is the iteration based on the previous values [6] (20)Tfz=GTz-AsinαJCm+AF+Aggsinα+exp⁡E,where (21)A=ΠDUκMMCκ+D2Clog⁡48TDtD-0.29-1,E=-z-z0AAsinαJCmiii·Tfz0-GTz0+AsinαJCm-AF-Aggsinα.F is a correction factor combined with hydrocarbon expansion from high pressure to low pressure during the temperature change, known as a Joule-Thomson effect. Heat transfer occurs between the wellbore fluid and the formation, overcoming resistances offered by the tubing wall, tubing-casing annulus, casing wall, and cement [19].The mixture heat capacity calculated depends on how many phases are considered. Consider(22)Cm=∑i=1n‍MiCi∑i=1n‍Mi.Heat capacity of each phase is a function of temperature and density usually estimated from the experimental data. Gambill correlations are valid, respectively, for oil and gas as it is shown [20]: (23)C=ϱ-1/20.338+0.00045T-460.


## 2.3. Inflow

The derivation of Darcy’s law is used in inflow calculation to determine the flow through the permeable geological layer. This equation is valid for liquid phase inflow at the perforation points, but bubble point pressure has not been taken under the consideration [5]: (24)Q=0.00708khpr-pwμBlog⁡re/rw-0.75+S.Drainage radius may be multiplied by some coefficients, depends on reservoir shape. Skin factor is estimated based on well test analysis. A positive value indicates there is a pressure decline in the near vicinity of well that is more than expected based on the radial flow equation but does not necessarily indicate formation damage [21]. Considering gas inflow, Darcy equation is adjusted in terms of gas properties: (25)Q=0.00708khPr2-Pw2μβTlog⁡re/rw-0.75+S.Inflow equations do not account for the phase change in solution gas reservoirs. In a solution gas drive, there is expansion of the hydrocarbons below bubble point, which is beneficial because it adds energy to the system. Gas liberation in a solution gas drive is also detrimental to oil production because it lowers the effective permeability of oil [22].


## 3. The Framework

The framework represents an idea of how to organize the modelling part to get the complete information about the whole physics during the production process. Solving every correlation and calculation needs to be organized properly. Every correlation is a function of data taken from producer or calculated by other correlations. There are few approaches involved, which usually meet in commercial solution and they are protected by copyrights against publicity. The authors would like to present optimised approach for production simulations, divided in stages.


### 3.1. Initialization

Considering the well length and number of physical data during the production, it is very important to establish the simulation points (meshing) in a project. All correlations should run in these points, so the number of data for one loop of simulation is significant. Transient analysis is highly recommended as well, especially during the gas lift or water injection procedures. Correlations use different unit systems and they create difficulties in the framework as well. There is no standardization to keep the modelling results. Flow area in the tubing affects the algorithms, so three options are available: tubing flow, annular flow, and tubing and annular flow. In some cases more than one tubing in a single well provides the production. There are wells where the tubing is split into two independent items from any depth. Sometimes only one casing string without any tubing is involved. It is highly desirable to have the temperature or pressure data from gauges across the well especially at the wellhead and bottom hole. Then, every meshing point is calculated upon this data. Finally, the complexity of modelling properties is relatively high in this kind of the framework.The producer supports the data which has been called as initial information (II). They are divided between sections as follows. Highly unlikely, all of them are given very precisely (for every depth), but from the simulation point of view the more the data is given, the better the results are obtained. Here, they are presented as follows:(i)
trajectory: this set usually has four values: true vertical depth (TVD), measured depth (MD), inclination (Incl), and azimuth (Az); this is the well geometry description;(ii)
geothermal temperature (GT): this is a trajectory function; some frameworks may be based on geothermal temperature especially in a very first part of calculation; finally, knowing this value is highly recommended in terms of the comparison with the real data; in many cases, the full understanding of geothermal temperature explains the well behaviour;(iii)
completion items: in the reality, any additional item changes the physical data at this particular depth; in particular, it is valid for packers and running: electric submersible pumps, gas lift valves, inflow control devices, sand screens, and so forth;(iv)
completion data: there are conductivity (Conn), roughness (R), inner diameter (ID), outer diameter (OD), number of casing strings, and cement properties; outer diameter, conductivity, and cement properties are important regarding the heat transfer and once the flow is observed in the annulus;(v)
mixture properties: it may contain gas/oil ratio (GOR) in surface condition information which affects the amount of gas out of the solution; alternatively, some calculations use gas/water rate (GWR); it also includes API gravity in surface condition as a measure of how heavy or light a petroleum liquid is compared to water and gas specified gravity (SG) in surface condition, which is the ratio of the density of a substance to the density of air and water density (WD);(vi)
layer data as depth (LD), permeability (k), pressure (LP), and reservoir extend (RE) are values usually given as constant for every layer; layer permeability and pressure have the crucial impact on the production results;(vii)
real data: the subjects involved in industry have the real data which can be used as a reference in the modelling process; it may contain wellhead temperature (WHT), bottom hole temperature (BHT), wellhead pressure (WHP), wellhead temperature (WHT), separator temperature (ST), separator pressure (SP), production time (t), and dorm pressure in the annulus (DormP);(viii)
other information: it is very flexible and cannot be specified to the regular mandatory set; it may contain many factors which has the important impact on the simulation results, that is, gamma ray data and surrounding well information; it is almost impossible to predict this information in the simulation model and usually this data is considered during the modelling interpretation.


### 3.2. Framework Steps

Five steps of framework are identified.Step 1. Get the producer data mentioned as II previously.Step 2. Fill the mandatory data as geothermal temperature, flow diameter, heat transfer, geothermal gradient, and TVD gradient.Step 3. Calculate correlations based on geothermal temperature and layer pressure.Step 4. Calculate rates, densities, and viscosities for every phase.Step 5. Find the final correlations based on temperature and pressure from the wellbore. This level is different from Step 3, because physical properties are calculated using the iteration algorithm instead of layers definition and geothermal temperature data.We present the framework in pseudocode based on object oriented programming.


### 3.3. Correlations

Once the literature review has been done, the most proper correlations regarding this framework have been chosen. We are not interested in showing the correlations algorithm here, but just in putting the emphasis on the dependence. Steps 3 and 5 use correlations with the different entry parameters. For Step 3, temperature (T) is the geothermal temperature (GT); from Step 2, pressure (p) is the layer pressure (LP); from Step 1, densities and other physical properties are taken from Steps 1 and 2 as well. For Step 5, temperature and pressure are iterated and density and other physical properties are taken from Steps 4 and 5. Here, the correlations dependence is shown as follows:(i)
BubblePointPressure(T, SepT, API, SG, SepP, GOR)(ii)
GasSolubility(T, API, p, SG)(iii)
DeadOilViscosity(T, API)(iv)
SaturatedOilViscosity(GasSolubility, DeadOilViscosity)(v)
UnderSaturatedOilViscosity(p, BubblePointPressure, SaturatedOilViscosity)(vi)
OilViscosity(DeadOilViscosity, SaturatedOilViscosity, UnderSaturatedOilViscosity, p, BubblePointPressure)(vii)
WaterViscosity(WD, T)(viii)
GasViscosity(T, SG)(ix)
GasCompressibility(T, SG, p)(x)
OilFormationVolumeFactor(T, API, GOR, SG)(xi)
WaterFormationVolumeFactor(T, p)(xii)
GasFormationVolumeFactor(T, GasCompressibility, p)(xiii)
OilHeatCapacity(API, T)(xiv)
WaterHeatCapacity(T, WD)(xv)
GasHeatCapacity(T, SG)(xvi)
OilInflow(RE, k, LP, p, OilFormationVolumeFactor, S)(xvii)
WaterInflow(RE, k, LP, p, WaterFormationVolumeFactor, S)(xviii)
GasInflow(SG, RE, k, LP, p, GasCompressibility, GasViscosity, GasFormationVolumeFactor, S)(xix)
SurfaceTension(T, p)(xx)
EarthThermalDiffusivity(T)(xxi)
FrictionFactor(ReNS, R, FlowDiameter)The temperature and pressure correlations are considered as a part of Step 5.


### 3.4. Initial Information: Step 1

In Step 1 there is nothing to do with the calculation. Only the producer, trajectory, bore, lift, geology, and density data are loaded and stored in one objectStep 1 (see Figure 2).Figure 2
Step 1: data flow.


### 3.5. Basic Properties: Step 2

This is the very first step when calculation is performed. Geothermal temperature is interpolated based on two crucial values: WHT and BHT with respect to MD and TVD data. As the results, geothermal temperature is given in Listing1.Listing 1:Geothermal temperature algorithm.
temperatureStep=(BHT−WHT)/MD[wellBottom]
TVDTemperature[ 0 ]=WHT
for i=0 To WellBottom
TVDTemperature[i+1]=TVDTemperature[i]+temperatureStep
for i=0 To WellBottom
GT[i]=TVDTemperature[TVD[i]]Then geothermal gradient (GG) is calculated as shown in Listing2.Listing 2:Geothermal gradient algorithm.
for i=0 To WellBottom
GG[i]=(GT[i+1]−GT[i])/(MD[i+1]−MD[i])Heat transfer [HT] calculation is shown in Listing3.Listing 3:Heat transfer algorithm.
for i=0 To WellBottom
HT[i]=(log(OD[i]/ID[i])/Conn[i])+(log(boreDiameter[i]/OD[i])/Conn[i])For the tubing flow, the flow diameter is equal to tubing ID.


### 3.6. Reservoir Properties: Step 3

Here the correlations are used based on reservoir properties. Passing the parameters to proper correlations, the algorithm returns results arrays. Hence, we have Listing4.Listing 4:Step 3 algorithm.
for i=0 To WellBottom [
(1)  BubblePointPressure[i]=correlation.BubblePointPressure(GT[i], SepT, API, SG, SepP, GOR)
(2)  GasSolubility[i]=correlation.GasSolubility(GT[i], API, LP[i], GOR)
(3)  GasViscosity[i]=correlation.GasViscosity(GT[i], SG)
(4)  DeadOilViscosity[i]=correlation.DeadOilViscosity(GT[i], API)
(5)  SaturatedOilViscosity[i]=correlation.SaturatedOilViscosity(GasSolubility[i], DeadOilViscosity[i])
(6)  UnderSaturatedOilViscosity[i]=correlation.UnderSaturatedOilViscosity(LP[i], BubblePointPressure[i],
SaturatedOilViscosity[i]
(7) OilViscosity[i]=correlation.OilViscosity(DeadOilViscosity[i], SaturatedOilViscosity[i],
UnderSaturatedOilViscosity[i], LP[i],  BubblePointPressure[i])
(8)  WaterViscosity[i]=correlation.WaterViscosity(WD, GT[i])
(9)  GasCompressibility[i]=correlation.GasCompressibility(GT[i], SG, LP[i])
(10)  GasFormationVolumeFactor[i]=correlation.GasFormationVolumeFactor(GT[i], GasCompressibility[i], LP[i])
(11)  OilFormationVolumeFactor[i]=correlation.OilFormationVolumeFactor(GT[i], API, GOR, SG)
(12)  WaterFormationVolumeFactor[i]=correlation.WaterFormationVolumeFactor(GT[i], API, GOR, SG)
(13)  OilHeatCapacity[i]=correlation.OilHeatCapacity(API, GT[i])
(14)  WaterHeatCapacity[i]=correlation.WaterHeatCapacity(GT[i], WD)
(15)  GasHeatCapacity[i]=correlation.GasHeatCapacity(GT[i], SG)
]
(16)  return new objectStep 3(BubblePoint, GasSolubility, GasViscosity, OilViscosity, DeadOilViscosity, SaturatedOilViscosity,
UnderSaturatedOilViscosity, WaterViscosity, GasCompressibility, GasFormationVolumeFactor, OilFormationVolumeFactor,
WaterFormationVolumeFactor, GasHeatCapacity, OilHeatCapacity, WaterHeatCapacity);As the results, objectStep3 contains correlations which have been used in inflow/production correlations in Step 4.


### 3.7. Flow Properties: Step 4

Here another portion of reservoir calculation is presented. It contains flow rates, mass flow rates, heat capacities, and densities, for all phases and their mixture with the separation for liquid and gas phase. The algorithm has to take into consideration the different unit systems to avoid any problems with the passing parameters to the next step (see Listing5).Listing 5:Step 4 algorithm.
for i=0 To WellBottom [
(1)TemperatureInflow[i]=GT[i]
(2)OilInflow[i]=correlation.OilInflow(RE[i], k[i], LP[i], OilFormationVolumeFactor[i], OilViscosity[i], S[i])
(3)WaterInflow[i]=correlation.WaterInflow(RE[i], k[i], LP[i],WaterFormationVolumeFactor[i], WaterViscosity[i], S[i])
(4)GasInflow[i]=correlation.GasInflow(SG, RE[i], k[i], GasCompressibility[i], GasViscosity[i], GT[i], S[i], LP[i], P[i])
(5)OilFlowRate[i]=OilInflow[i]
(6)WaterFlowRate[i]=WaterInflow[i]
(7)GasFlowRate[i]=GasInflow[i]*GOR
(8)if (GLVExist == true) GasFlowRate[i]+=GasInflowFromGLV
(9)LiquidFlowRate[i]=WaterFlowRate[i]+OilFlowRate[i]
(10)GasMassFlowRate[i]=GasFlowRate[i]*SG
(11)OilMassFlowRate[i]=OilFlowRate[i]*API
(12)WaterMassFlowRate[i]=WaterFlowRate[i]*WD
(13)LiquidMassFlowRate[i]=OilMassFlowRate[i]+WaterMassFlowRate[i]
(14)MixtureMassFlowRate[i]=LiquidMassFlowRate[i]+GasMassFlowRate[i]
(15)LiquidViscosity[i]=(OilMassFlowRate[i]*Step  3.OilViscosity[i]+WaterMassFlowRate[i]*Step  3.WaterViscosity[i])/
(OilMassFlowRate[i]+WaterMassFlowRate[i])
(16)MixtureViscosity[i]=(OilMassFlowRate[i]*Step  3.OilViscosity[i]+WaterMassFlowRate[i]*Step  3.WaterViscosity[i]+
GasMassFlowRate[i]
*
Step  3.GasViscosity[i])/(OilMassFlowRate[i]+WaterMassFlowRate[i]+GasMassFlowRate[i])
(17)LiquidHeatCapacity[i]=(OilMassFlowRate[i]*Step  3.OilHeatCapacity[i]+WaterMassFlowRate[i]*
Step  3.WaterHeatCapacity[i])/(OilMassFlowRate[i]+WaterMassFlowRate[i])
(18)MixtureHeatCapacity[i]=(OilMassFlowRate[i]*Step  3.OilHeatCapacity[i]+WaterMassFlowRate[i]*
Step  3.WaterHeatCapacity[i]+GasMassFlowRate[i]
*
Step  3.GasHeatCapacity[i])/(OilMassFlowRate[i]+
WaterMassFlowRate[i]+GasMassFlowRate[i])
(19)LiquidDensity[i]=(OilMassFlowRate[i]*API+WaterMassFlowRate[i]*WD)/(OilMassFlowRate[i]+
WaterMassFlowRate[i])
(20)MixtureDensity[i]=(OilMassFlowRate[i]*API+WaterMassFlowRate[i]*WD+GasMassFlowRate[i]*SG)/
(OilMassFlowRate[i]+WaterMassFlowRate[i]+GasMassFlowRate[i])
]
return new objectStep  4(…)Lines 5, 6, and 7 are dependent on producer type. Because the gas lift in current case was involved, only GasFlowRate is the function of GOR. In line 8, the condition is created to adjust the amount of lifted gas by gas lift valve into the wellbore.


### 3.8. Wellbore Properties: Step 5

Based on the simply thermodynamic principles, it is obvious that gaseous and liquid states not only are merged into each other in a continuous manner, but also are in fact similar in nature. Volumes of molecules and the intermolecular forces are necessary in establishing the relationship between pressure, volume, and temperature of gases and liquids. So in the foundation of this model Van der Waals equation is used:(26)p+av2v-b=kT.Here the final step is presented. This algorithm is the iteration, so the initial temperature and pressure values are taken from the bottom hole gauge. In this step, temperature and pressure correlations are dependent as follows: (i)
T(T[i+1], GT[i+1], Step4.MixtureHeatCapacity, Step4.MixtureMassFlowRate, t, EarthThermalDiffusivity, GT[i], BoreRadius, FlowDiameter, HT, Incl, GG)(ii)
P(P[i+1], Step4.OilFlowRate, Step4.WaterFlowRate, Step4.GasFlowRate, Step3.OilFormationVolumeFactor, Step3.WaterFormationVolumeFactor, Step3.GasFormationVolumeFactor, Step5.OilViscosity, Step5.WaterViscosity, Step5.SurfaceTension, depth, FlowDirection, Step5.WaterDensity, Step5.OilDensity, Step5.GasDensity, FlowDiameter, Incl, gradTVD, R)Pressure and temperature are calculated from the bottom to the top of the well, as the flow occurs, so index i+1 means the previous correlations results. Additionally the pressure correlation returns object which contains the following: return new objectPressure(P, FlowDirection, FlowType, LiquidSuperficialVelocity, GasSuperficialVelocity, MixtureVelocity, LiquidVolumeFraction, KineticPressureDrop, ReNS, HydrostaticPressureLoss, FrictionFactor)As a FlowType we consider segregated, intermittent, distributed, and transient. Every type has its own flow regimes which is the part of the pressure correlation. The important part of this step is to find the amount of gas which is out of the solution. The relation between solubility of liquid components and the heat of solution is given as (27)δln⁡xiδ1/TP=-ΔslnhiR,where xi is the mole fraction of ith alkane in water and Δslnhi is the difference between the partial enthalpy of the ith hydrocarbon at infinite dilution and the molar enthalpy of the pure hydrocarbon. The heat of solution includes two effects: positive heat of cavity formation and negative heat of hydrophobic interaction between the hydrocarbon and water. These two effects cancel each other at Tm. So the description of solubility of hydrocarbons in water may be presented as (28)ln⁡xiT=ln⁡xiTm+ΔslnCP,iRln⁡TTm+TmT-1.In this step, phase intermixing of viscosity, heat capacity, and rates is also considered.Finally Step 5 is presenetd as shown in Listing6.Listing 6:Step 5 algorithm.
(1) P[wellBottom]=BHP
(2) T[wellBottom]=BHT
(3) for WellBottom-1 To i=0 [
(4) GasDensity[i]=Step 4.GassMassFlowRate[i]*SG
(5) OilDensity[i]=Step 4.OilMassFlowRate[i]*API
(6) WaterDensity[i]=Step 4.WaterMassFlowRate[i]*WD
(7) SumGasFlowRate[i]=GasFlowRate[i]+SumGasFlowRate[i+1]
(8) SumOilFlowRate[i]=OilFlowRate[i]+SumOilFlowRate[i+1]
(9) SumWaterFlowRate[i]=WaterFlowRate[i]+SumWaterFlowRate[i+1]
(10) SumLiquidFlowRate[i]=SumOilFlowRate[i]+SumWaterFlowRate[i]
(11) GOR[i]=GasFlowRate[i]/OilFlowRate[i]
(12) GasSolubility[i]=correlation.GasSolubility(T[i+1], OilDensity[i], P[i+1], GasDensity[i+1])
(13) GasCompressibility[i]=correlation.gasCompressibility(T[i+1], GasDensity[i], P[i+1])
(14) GasFormationVolumeFactor[i]=correlation.GasFormationVolumeFactor(T[i+1], GasCompressibility[i], P[i+1])
(15) OilFormationVolumeFactor[i]=correlation.OilFormationVolumeFactor(T[i+1], OilDensity[i], GOR[i], GasDensity[i])
(16) WaterFormationVolumeFactor[i]=correlation.WaterFormationVolumeFactor(T[i+1], P[i+1])
(17) EarthThermalDiffusivity[i]=correlation.EarthThermalDiffusivity(T[i+1])
(18) SumGasFlowRate[i]*=GasFormationVolumeFactor[i]
(19) SumOilFlowRate[i]*=OilFormationVolumeFactor[i]
(20) SumWaterFlowRate[i]*=WaterFormationVolumeFactor[i]
(21) SumLiquidFlowRate[i]=SumOilFlowRate[i]+SumWaterFlowRate[i]
(22) if (GasSolubility[i] < GOR[i])
FreeGasFlowRate[i]=GasFlowRate[i]−(GasSolubility[i]*OilFlowRate[i])
else FreeGasFlowRate[i]=0
(23) GasInSolutionFlowRates[i]=Step 4.GasFlowRate[i]−FreeGasFlowRate[i]
(24) BubblePointPressure[i]=correlation.BubblePointPressure(T[i], SepT, OilDensity[i], GasDensity[i], SepP, GOR[i])
(25) GasViscosity[i]=correlation.GasViscosity(T[i+1], GasDensity[i])
(26) DeadOilViscosity[i]=correlation.DeadOilViscosity(T[i+1], OilDesnity[i])
(27) SaturatedOilViscosity[i]=correlation.SaturatedOilViscosity(GasSolubility[i], DeadOilViscosity[i])
(28) UnderSaturatedOilViscosity[i]=correlation.UnderSaturatedOilViscosity(P[i+1], BubblePointPressure[i],
SaturatedOilViscosity[i]
(29) OilViscosity[i]=correlation.OilViscosity(DeadOilViscosity[i], SaturatedOilViscosity[i], UnderSaturatedOilViscosity[i],
P[i+1], BubblePointPressure[i])
(30) WaterViscosity[i]=correlation.WaterViscosity(WaterDensity[i], T[i+1])
(31) OilHeatCapacity[i]=correlation.OilHeatCapacity(OilDensity[i], T[i+1])
(32) WaterHeatCapacity[i]=correlation.WaterHeatCapacity(T[i+1], WaterDensity[i])
(33) GasHeatCapacity[i]=correlation.GasHeatCapacity(T[i+1], GasDensity[i])
(34) LiquidViscosity[i]=(Step 4.OilMassFlowRate[i]*Step 5.OilViscosity[i]+
Step 4.WaterMassFlowRate[i]
*
Step 5.WaterViscosity[i])/(Step 4.OilMassFlowRate[i]+Step 4.WaterMassFlowRate[i])
(35) MixtureViscosity[i]=(Step 4.OilMassFlowRate[i]*Step 5.OilViscosity[i]+
Step 4.WaterMassFlowRate[i]
*
Step 5.WaterViscosity[i]+Step 4.GasMassFlowRate[i]
*
Step 5.GasViscosity[i])/(Step 4.OilMassFlowRate[i]+Step 4.WaterMassFlowRate[i]+Step 4.GasMassFlowRate[i])
(36) LiquidHeatCapacity[i]=(Step 4.OilMassFlowRate[i]*Step 5.OilHeatCapacity[i]+Step 4.WaterMassFlowRate[i]*
Step 5.WaterHeatCapacity[i])/(Step 4.OilMassFlowRate[i]+Step 4.WaterMassFlowRate[i])
(37) MixtureHeatCapacity[i]=(Step 4.OilMassFlowRate[i]*Step 5.OilHeatCapacity[i]+
Step 4.WaterMassFlowRate[i]
*
Step 5.WaterHeatCapacity[i]
+Step 4.GasMassFlowRate[i]
*
Step 5.GasHeatCapacity[i])/(Step 4.OilMassFlowRate[i]
+Step 4.WaterMassFlowRate[i]+Step 4.GasMassFlowRate[i])
(38) LiquidDensity[i]=(Step 4.OilMassFlowRate[i]*OilDensity[i]+Step 4.WaterMassFlowRate[i]*
WaterDensity[i])/(Step 4.OilMassFlowRate[i]+Step 4.WaterMassFlowRate[i])
(39) MixtureDensity[i]=(Step 4.OilMassFlowRate[i]*OilDensity[i]+Step 4.WaterMassFlowRate[i]*
WaterDensity+Step 4.GasMassFlowRate[i]
*
GasDensity[i])/(Step 4.OilMassFlowRate[i]
+Step 4.WaterMassFlowRate[i]+Step 4.GasMassFlowRate[i])
(40) SumGasMassFlowRate[i]+=Step 4.GasMassFlowRate[i]+SumGasMassFlowRate[i+1]
(41) SumOilMassFlowRate[i]+=Step 4.OilMassFlowRate[i]+SumOilMassFlowRate[i+1]
(42) SumWaterMassFlowRate[i]+=Step 4.WaterMassFlowRate[i]+SumWaterMassFlowRate[i+1]
(43) MixtureMassFlowRate[i]=SumGasMassFlowRate[i]+SumOilMassFlowRate[i]+SumWaterMassFlowRate[i]
(44) T[i]=correlation.Temperature(…)
(45) SurfaceTension[i]=correlation.SurfaceTension(T[i], P[i+1])
(46) P[i]=correlation.Pressure(…)
]
return new objectStep 5(accumulatedGasMassFlowRate, accumulatedOilMassFlowRate, accumulatedWaterMassFlowRate, gasD,
oilD, watD, gasSolubility, gor, gasCompressibility, GFVF, OFVF, WFVF, bubblePoint, gasViscosity, oilViscosity, liquidViscosity,
mixtureViscosity, waterViscosity, oilHeatCapacity,  waterHeatCapacity, gasHeatCapacity,  T, SurfaceTension, P);


### 3.9. Final Join

All these steps have to be joined in the final calculation which solves all modelling stages. Considering the whole oilfield, this simulation does not take into consideration limited amount of wells. Every well can be treated as a single thread calculation until the gas lift is considered as an optimisation problem for the oilfield. Then, the proper gas distribution is dependent on every well simulation. This study is planned as a future work.Hence, we have Listing7.Listing 7:Final join algorithm.
for (i=0 To numberOfWells) [
(1)  Step 2.Add(Step 2.Run(Step 1)
(2)  Step 3.Add(Step 2.Run(oilFiledData, Step 2))
(3)  InitP=InitP(Step 1)
(4)  Step 4.Add(Step 4.Run(Step 3, Step 2, InitP, GasLift)
(5)  Step 5.Add(Step 5.Run(Step 1, Step 2, Step 4, SG, API, WD))
]FunctionInitP is the initial pressure calculation based on WHP i BHP data with the consideration of trajectory curvature.


## 3.1. Initialization

Considering the well length and number of physical data during the production, it is very important to establish the simulation points (meshing) in a project. All correlations should run in these points, so the number of data for one loop of simulation is significant. Transient analysis is highly recommended as well, especially during the gas lift or water injection procedures. Correlations use different unit systems and they create difficulties in the framework as well. There is no standardization to keep the modelling results. Flow area in the tubing affects the algorithms, so three options are available: tubing flow, annular flow, and tubing and annular flow. In some cases more than one tubing in a single well provides the production. There are wells where the tubing is split into two independent items from any depth. Sometimes only one casing string without any tubing is involved. It is highly desirable to have the temperature or pressure data from gauges across the well especially at the wellhead and bottom hole. Then, every meshing point is calculated upon this data. Finally, the complexity of modelling properties is relatively high in this kind of the framework.The producer supports the data which has been called as initial information (II). They are divided between sections as follows. Highly unlikely, all of them are given very precisely (for every depth), but from the simulation point of view the more the data is given, the better the results are obtained. Here, they are presented as follows:(i)
trajectory: this set usually has four values: true vertical depth (TVD), measured depth (MD), inclination (Incl), and azimuth (Az); this is the well geometry description;(ii)
geothermal temperature (GT): this is a trajectory function; some frameworks may be based on geothermal temperature especially in a very first part of calculation; finally, knowing this value is highly recommended in terms of the comparison with the real data; in many cases, the full understanding of geothermal temperature explains the well behaviour;(iii)
completion items: in the reality, any additional item changes the physical data at this particular depth; in particular, it is valid for packers and running: electric submersible pumps, gas lift valves, inflow control devices, sand screens, and so forth;(iv)
completion data: there are conductivity (Conn), roughness (R), inner diameter (ID), outer diameter (OD), number of casing strings, and cement properties; outer diameter, conductivity, and cement properties are important regarding the heat transfer and once the flow is observed in the annulus;(v)
mixture properties: it may contain gas/oil ratio (GOR) in surface condition information which affects the amount of gas out of the solution; alternatively, some calculations use gas/water rate (GWR); it also includes API gravity in surface condition as a measure of how heavy or light a petroleum liquid is compared to water and gas specified gravity (SG) in surface condition, which is the ratio of the density of a substance to the density of air and water density (WD);(vi)
layer data as depth (LD), permeability (k), pressure (LP), and reservoir extend (RE) are values usually given as constant for every layer; layer permeability and pressure have the crucial impact on the production results;(vii)
real data: the subjects involved in industry have the real data which can be used as a reference in the modelling process; it may contain wellhead temperature (WHT), bottom hole temperature (BHT), wellhead pressure (WHP), wellhead temperature (WHT), separator temperature (ST), separator pressure (SP), production time (t), and dorm pressure in the annulus (DormP);(viii)
other information: it is very flexible and cannot be specified to the regular mandatory set; it may contain many factors which has the important impact on the simulation results, that is, gamma ray data and surrounding well information; it is almost impossible to predict this information in the simulation model and usually this data is considered during the modelling interpretation.


## 3.2. Framework Steps

Five steps of framework are identified.Step 1. Get the producer data mentioned as II previously.Step 2. Fill the mandatory data as geothermal temperature, flow diameter, heat transfer, geothermal gradient, and TVD gradient.Step 3. Calculate correlations based on geothermal temperature and layer pressure.Step 4. Calculate rates, densities, and viscosities for every phase.Step 5. Find the final correlations based on temperature and pressure from the wellbore. This level is different from Step 3, because physical properties are calculated using the iteration algorithm instead of layers definition and geothermal temperature data.We present the framework in pseudocode based on object oriented programming.


## 3.3. Correlations

Once the literature review has been done, the most proper correlations regarding this framework have been chosen. We are not interested in showing the correlations algorithm here, but just in putting the emphasis on the dependence. Steps 3 and 5 use correlations with the different entry parameters. For Step 3, temperature (T) is the geothermal temperature (GT); from Step 2, pressure (p) is the layer pressure (LP); from Step 1, densities and other physical properties are taken from Steps 1 and 2 as well. For Step 5, temperature and pressure are iterated and density and other physical properties are taken from Steps 4 and 5. Here, the correlations dependence is shown as follows:(i)
BubblePointPressure(T, SepT, API, SG, SepP, GOR)(ii)
GasSolubility(T, API, p, SG)(iii)
DeadOilViscosity(T, API)(iv)
SaturatedOilViscosity(GasSolubility, DeadOilViscosity)(v)
UnderSaturatedOilViscosity(p, BubblePointPressure, SaturatedOilViscosity)(vi)
OilViscosity(DeadOilViscosity, SaturatedOilViscosity, UnderSaturatedOilViscosity, p, BubblePointPressure)(vii)
WaterViscosity(WD, T)(viii)
GasViscosity(T, SG)(ix)
GasCompressibility(T, SG, p)(x)
OilFormationVolumeFactor(T, API, GOR, SG)(xi)
WaterFormationVolumeFactor(T, p)(xii)
GasFormationVolumeFactor(T, GasCompressibility, p)(xiii)
OilHeatCapacity(API, T)(xiv)
WaterHeatCapacity(T, WD)(xv)
GasHeatCapacity(T, SG)(xvi)
OilInflow(RE, k, LP, p, OilFormationVolumeFactor, S)(xvii)
WaterInflow(RE, k, LP, p, WaterFormationVolumeFactor, S)(xviii)
GasInflow(SG, RE, k, LP, p, GasCompressibility, GasViscosity, GasFormationVolumeFactor, S)(xix)
SurfaceTension(T, p)(xx)
EarthThermalDiffusivity(T)(xxi)
FrictionFactor(ReNS, R, FlowDiameter)The temperature and pressure correlations are considered as a part of Step 5.


## 3.4. Initial Information: Step 1

In Step 1 there is nothing to do with the calculation. Only the producer, trajectory, bore, lift, geology, and density data are loaded and stored in one objectStep 1 (see Figure 2).Figure 2
Step 1: data flow.


## 3.5. Basic Properties: Step 2

This is the very first step when calculation is performed. Geothermal temperature is interpolated based on two crucial values: WHT and BHT with respect to MD and TVD data. As the results, geothermal temperature is given in Listing1.Listing 1:Geothermal temperature algorithm.
temperatureStep=(BHT−WHT)/MD[wellBottom]
TVDTemperature[ 0 ]=WHT
for i=0 To WellBottom
TVDTemperature[i+1]=TVDTemperature[i]+temperatureStep
for i=0 To WellBottom
GT[i]=TVDTemperature[TVD[i]]Then geothermal gradient (GG) is calculated as shown in Listing2.Listing 2:Geothermal gradient algorithm.
for i=0 To WellBottom
GG[i]=(GT[i+1]−GT[i])/(MD[i+1]−MD[i])Heat transfer [HT] calculation is shown in Listing3.Listing 3:Heat transfer algorithm.
for i=0 To WellBottom
HT[i]=(log(OD[i]/ID[i])/Conn[i])+(log(boreDiameter[i]/OD[i])/Conn[i])For the tubing flow, the flow diameter is equal to tubing ID.


## 3.6. Reservoir Properties: Step 3

Here the correlations are used based on reservoir properties. Passing the parameters to proper correlations, the algorithm returns results arrays. Hence, we have Listing4.Listing 4:Step 3 algorithm.
for i=0 To WellBottom [
(1)  BubblePointPressure[i]=correlation.BubblePointPressure(GT[i], SepT, API, SG, SepP, GOR)
(2)  GasSolubility[i]=correlation.GasSolubility(GT[i], API, LP[i], GOR)
(3)  GasViscosity[i]=correlation.GasViscosity(GT[i], SG)
(4)  DeadOilViscosity[i]=correlation.DeadOilViscosity(GT[i], API)
(5)  SaturatedOilViscosity[i]=correlation.SaturatedOilViscosity(GasSolubility[i], DeadOilViscosity[i])
(6)  UnderSaturatedOilViscosity[i]=correlation.UnderSaturatedOilViscosity(LP[i], BubblePointPressure[i],
SaturatedOilViscosity[i]
(7) OilViscosity[i]=correlation.OilViscosity(DeadOilViscosity[i], SaturatedOilViscosity[i],
UnderSaturatedOilViscosity[i], LP[i],  BubblePointPressure[i])
(8)  WaterViscosity[i]=correlation.WaterViscosity(WD, GT[i])
(9)  GasCompressibility[i]=correlation.GasCompressibility(GT[i], SG, LP[i])
(10)  GasFormationVolumeFactor[i]=correlation.GasFormationVolumeFactor(GT[i], GasCompressibility[i], LP[i])
(11)  OilFormationVolumeFactor[i]=correlation.OilFormationVolumeFactor(GT[i], API, GOR, SG)
(12)  WaterFormationVolumeFactor[i]=correlation.WaterFormationVolumeFactor(GT[i], API, GOR, SG)
(13)  OilHeatCapacity[i]=correlation.OilHeatCapacity(API, GT[i])
(14)  WaterHeatCapacity[i]=correlation.WaterHeatCapacity(GT[i], WD)
(15)  GasHeatCapacity[i]=correlation.GasHeatCapacity(GT[i], SG)
]
(16)  return new objectStep 3(BubblePoint, GasSolubility, GasViscosity, OilViscosity, DeadOilViscosity, SaturatedOilViscosity,
UnderSaturatedOilViscosity, WaterViscosity, GasCompressibility, GasFormationVolumeFactor, OilFormationVolumeFactor,
WaterFormationVolumeFactor, GasHeatCapacity, OilHeatCapacity, WaterHeatCapacity);As the results, objectStep3 contains correlations which have been used in inflow/production correlations in Step 4.


## 3.7. Flow Properties: Step 4

Here another portion of reservoir calculation is presented. It contains flow rates, mass flow rates, heat capacities, and densities, for all phases and their mixture with the separation for liquid and gas phase. The algorithm has to take into consideration the different unit systems to avoid any problems with the passing parameters to the next step (see Listing5).Listing 5:Step 4 algorithm.
for i=0 To WellBottom [
(1)TemperatureInflow[i]=GT[i]
(2)OilInflow[i]=correlation.OilInflow(RE[i], k[i], LP[i], OilFormationVolumeFactor[i], OilViscosity[i], S[i])
(3)WaterInflow[i]=correlation.WaterInflow(RE[i], k[i], LP[i],WaterFormationVolumeFactor[i], WaterViscosity[i], S[i])
(4)GasInflow[i]=correlation.GasInflow(SG, RE[i], k[i], GasCompressibility[i], GasViscosity[i], GT[i], S[i], LP[i], P[i])
(5)OilFlowRate[i]=OilInflow[i]
(6)WaterFlowRate[i]=WaterInflow[i]
(7)GasFlowRate[i]=GasInflow[i]*GOR
(8)if (GLVExist == true) GasFlowRate[i]+=GasInflowFromGLV
(9)LiquidFlowRate[i]=WaterFlowRate[i]+OilFlowRate[i]
(10)GasMassFlowRate[i]=GasFlowRate[i]*SG
(11)OilMassFlowRate[i]=OilFlowRate[i]*API
(12)WaterMassFlowRate[i]=WaterFlowRate[i]*WD
(13)LiquidMassFlowRate[i]=OilMassFlowRate[i]+WaterMassFlowRate[i]
(14)MixtureMassFlowRate[i]=LiquidMassFlowRate[i]+GasMassFlowRate[i]
(15)LiquidViscosity[i]=(OilMassFlowRate[i]*Step  3.OilViscosity[i]+WaterMassFlowRate[i]*Step  3.WaterViscosity[i])/
(OilMassFlowRate[i]+WaterMassFlowRate[i])
(16)MixtureViscosity[i]=(OilMassFlowRate[i]*Step  3.OilViscosity[i]+WaterMassFlowRate[i]*Step  3.WaterViscosity[i]+
GasMassFlowRate[i]
*
Step  3.GasViscosity[i])/(OilMassFlowRate[i]+WaterMassFlowRate[i]+GasMassFlowRate[i])
(17)LiquidHeatCapacity[i]=(OilMassFlowRate[i]*Step  3.OilHeatCapacity[i]+WaterMassFlowRate[i]*
Step  3.WaterHeatCapacity[i])/(OilMassFlowRate[i]+WaterMassFlowRate[i])
(18)MixtureHeatCapacity[i]=(OilMassFlowRate[i]*Step  3.OilHeatCapacity[i]+WaterMassFlowRate[i]*
Step  3.WaterHeatCapacity[i]+GasMassFlowRate[i]
*
Step  3.GasHeatCapacity[i])/(OilMassFlowRate[i]+
WaterMassFlowRate[i]+GasMassFlowRate[i])
(19)LiquidDensity[i]=(OilMassFlowRate[i]*API+WaterMassFlowRate[i]*WD)/(OilMassFlowRate[i]+
WaterMassFlowRate[i])
(20)MixtureDensity[i]=(OilMassFlowRate[i]*API+WaterMassFlowRate[i]*WD+GasMassFlowRate[i]*SG)/
(OilMassFlowRate[i]+WaterMassFlowRate[i]+GasMassFlowRate[i])
]
return new objectStep  4(…)Lines 5, 6, and 7 are dependent on producer type. Because the gas lift in current case was involved, only GasFlowRate is the function of GOR. In line 8, the condition is created to adjust the amount of lifted gas by gas lift valve into the wellbore.


## 3.8. Wellbore Properties: Step 5

Based on the simply thermodynamic principles, it is obvious that gaseous and liquid states not only are merged into each other in a continuous manner, but also are in fact similar in nature. Volumes of molecules and the intermolecular forces are necessary in establishing the relationship between pressure, volume, and temperature of gases and liquids. So in the foundation of this model Van der Waals equation is used:(26)p+av2v-b=kT.Here the final step is presented. This algorithm is the iteration, so the initial temperature and pressure values are taken from the bottom hole gauge. In this step, temperature and pressure correlations are dependent as follows: (i)
T(T[i+1], GT[i+1], Step4.MixtureHeatCapacity, Step4.MixtureMassFlowRate, t, EarthThermalDiffusivity, GT[i], BoreRadius, FlowDiameter, HT, Incl, GG)(ii)
P(P[i+1], Step4.OilFlowRate, Step4.WaterFlowRate, Step4.GasFlowRate, Step3.OilFormationVolumeFactor, Step3.WaterFormationVolumeFactor, Step3.GasFormationVolumeFactor, Step5.OilViscosity, Step5.WaterViscosity, Step5.SurfaceTension, depth, FlowDirection, Step5.WaterDensity, Step5.OilDensity, Step5.GasDensity, FlowDiameter, Incl, gradTVD, R)Pressure and temperature are calculated from the bottom to the top of the well, as the flow occurs, so index i+1 means the previous correlations results. Additionally the pressure correlation returns object which contains the following: return new objectPressure(P, FlowDirection, FlowType, LiquidSuperficialVelocity, GasSuperficialVelocity, MixtureVelocity, LiquidVolumeFraction, KineticPressureDrop, ReNS, HydrostaticPressureLoss, FrictionFactor)As a FlowType we consider segregated, intermittent, distributed, and transient. Every type has its own flow regimes which is the part of the pressure correlation. The important part of this step is to find the amount of gas which is out of the solution. The relation between solubility of liquid components and the heat of solution is given as (27)δln⁡xiδ1/TP=-ΔslnhiR,where xi is the mole fraction of ith alkane in water and Δslnhi is the difference between the partial enthalpy of the ith hydrocarbon at infinite dilution and the molar enthalpy of the pure hydrocarbon. The heat of solution includes two effects: positive heat of cavity formation and negative heat of hydrophobic interaction between the hydrocarbon and water. These two effects cancel each other at Tm. So the description of solubility of hydrocarbons in water may be presented as (28)ln⁡xiT=ln⁡xiTm+ΔslnCP,iRln⁡TTm+TmT-1.In this step, phase intermixing of viscosity, heat capacity, and rates is also considered.Finally Step 5 is presenetd as shown in Listing6.Listing 6:Step 5 algorithm.
(1) P[wellBottom]=BHP
(2) T[wellBottom]=BHT
(3) for WellBottom-1 To i=0 [
(4) GasDensity[i]=Step 4.GassMassFlowRate[i]*SG
(5) OilDensity[i]=Step 4.OilMassFlowRate[i]*API
(6) WaterDensity[i]=Step 4.WaterMassFlowRate[i]*WD
(7) SumGasFlowRate[i]=GasFlowRate[i]+SumGasFlowRate[i+1]
(8) SumOilFlowRate[i]=OilFlowRate[i]+SumOilFlowRate[i+1]
(9) SumWaterFlowRate[i]=WaterFlowRate[i]+SumWaterFlowRate[i+1]
(10) SumLiquidFlowRate[i]=SumOilFlowRate[i]+SumWaterFlowRate[i]
(11) GOR[i]=GasFlowRate[i]/OilFlowRate[i]
(12) GasSolubility[i]=correlation.GasSolubility(T[i+1], OilDensity[i], P[i+1], GasDensity[i+1])
(13) GasCompressibility[i]=correlation.gasCompressibility(T[i+1], GasDensity[i], P[i+1])
(14) GasFormationVolumeFactor[i]=correlation.GasFormationVolumeFactor(T[i+1], GasCompressibility[i], P[i+1])
(15) OilFormationVolumeFactor[i]=correlation.OilFormationVolumeFactor(T[i+1], OilDensity[i], GOR[i], GasDensity[i])
(16) WaterFormationVolumeFactor[i]=correlation.WaterFormationVolumeFactor(T[i+1], P[i+1])
(17) EarthThermalDiffusivity[i]=correlation.EarthThermalDiffusivity(T[i+1])
(18) SumGasFlowRate[i]*=GasFormationVolumeFactor[i]
(19) SumOilFlowRate[i]*=OilFormationVolumeFactor[i]
(20) SumWaterFlowRate[i]*=WaterFormationVolumeFactor[i]
(21) SumLiquidFlowRate[i]=SumOilFlowRate[i]+SumWaterFlowRate[i]
(22) if (GasSolubility[i] < GOR[i])
FreeGasFlowRate[i]=GasFlowRate[i]−(GasSolubility[i]*OilFlowRate[i])
else FreeGasFlowRate[i]=0
(23) GasInSolutionFlowRates[i]=Step 4.GasFlowRate[i]−FreeGasFlowRate[i]
(24) BubblePointPressure[i]=correlation.BubblePointPressure(T[i], SepT, OilDensity[i], GasDensity[i], SepP, GOR[i])
(25) GasViscosity[i]=correlation.GasViscosity(T[i+1], GasDensity[i])
(26) DeadOilViscosity[i]=correlation.DeadOilViscosity(T[i+1], OilDesnity[i])
(27) SaturatedOilViscosity[i]=correlation.SaturatedOilViscosity(GasSolubility[i], DeadOilViscosity[i])
(28) UnderSaturatedOilViscosity[i]=correlation.UnderSaturatedOilViscosity(P[i+1], BubblePointPressure[i],
SaturatedOilViscosity[i]
(29) OilViscosity[i]=correlation.OilViscosity(DeadOilViscosity[i], SaturatedOilViscosity[i], UnderSaturatedOilViscosity[i],
P[i+1], BubblePointPressure[i])
(30) WaterViscosity[i]=correlation.WaterViscosity(WaterDensity[i], T[i+1])
(31) OilHeatCapacity[i]=correlation.OilHeatCapacity(OilDensity[i], T[i+1])
(32) WaterHeatCapacity[i]=correlation.WaterHeatCapacity(T[i+1], WaterDensity[i])
(33) GasHeatCapacity[i]=correlation.GasHeatCapacity(T[i+1], GasDensity[i])
(34) LiquidViscosity[i]=(Step 4.OilMassFlowRate[i]*Step 5.OilViscosity[i]+
Step 4.WaterMassFlowRate[i]
*
Step 5.WaterViscosity[i])/(Step 4.OilMassFlowRate[i]+Step 4.WaterMassFlowRate[i])
(35) MixtureViscosity[i]=(Step 4.OilMassFlowRate[i]*Step 5.OilViscosity[i]+
Step 4.WaterMassFlowRate[i]
*
Step 5.WaterViscosity[i]+Step 4.GasMassFlowRate[i]
*
Step 5.GasViscosity[i])/(Step 4.OilMassFlowRate[i]+Step 4.WaterMassFlowRate[i]+Step 4.GasMassFlowRate[i])
(36) LiquidHeatCapacity[i]=(Step 4.OilMassFlowRate[i]*Step 5.OilHeatCapacity[i]+Step 4.WaterMassFlowRate[i]*
Step 5.WaterHeatCapacity[i])/(Step 4.OilMassFlowRate[i]+Step 4.WaterMassFlowRate[i])
(37) MixtureHeatCapacity[i]=(Step 4.OilMassFlowRate[i]*Step 5.OilHeatCapacity[i]+
Step 4.WaterMassFlowRate[i]
*
Step 5.WaterHeatCapacity[i]
+Step 4.GasMassFlowRate[i]
*
Step 5.GasHeatCapacity[i])/(Step 4.OilMassFlowRate[i]
+Step 4.WaterMassFlowRate[i]+Step 4.GasMassFlowRate[i])
(38) LiquidDensity[i]=(Step 4.OilMassFlowRate[i]*OilDensity[i]+Step 4.WaterMassFlowRate[i]*
WaterDensity[i])/(Step 4.OilMassFlowRate[i]+Step 4.WaterMassFlowRate[i])
(39) MixtureDensity[i]=(Step 4.OilMassFlowRate[i]*OilDensity[i]+Step 4.WaterMassFlowRate[i]*
WaterDensity+Step 4.GasMassFlowRate[i]
*
GasDensity[i])/(Step 4.OilMassFlowRate[i]
+Step 4.WaterMassFlowRate[i]+Step 4.GasMassFlowRate[i])
(40) SumGasMassFlowRate[i]+=Step 4.GasMassFlowRate[i]+SumGasMassFlowRate[i+1]
(41) SumOilMassFlowRate[i]+=Step 4.OilMassFlowRate[i]+SumOilMassFlowRate[i+1]
(42) SumWaterMassFlowRate[i]+=Step 4.WaterMassFlowRate[i]+SumWaterMassFlowRate[i+1]
(43) MixtureMassFlowRate[i]=SumGasMassFlowRate[i]+SumOilMassFlowRate[i]+SumWaterMassFlowRate[i]
(44) T[i]=correlation.Temperature(…)
(45) SurfaceTension[i]=correlation.SurfaceTension(T[i], P[i+1])
(46) P[i]=correlation.Pressure(…)
]
return new objectStep 5(accumulatedGasMassFlowRate, accumulatedOilMassFlowRate, accumulatedWaterMassFlowRate, gasD,
oilD, watD, gasSolubility, gor, gasCompressibility, GFVF, OFVF, WFVF, bubblePoint, gasViscosity, oilViscosity, liquidViscosity,
mixtureViscosity, waterViscosity, oilHeatCapacity,  waterHeatCapacity, gasHeatCapacity,  T, SurfaceTension, P);


## 3.9. Final Join

All these steps have to be joined in the final calculation which solves all modelling stages. Considering the whole oilfield, this simulation does not take into consideration limited amount of wells. Every well can be treated as a single thread calculation until the gas lift is considered as an optimisation problem for the oilfield. Then, the proper gas distribution is dependent on every well simulation. This study is planned as a future work.Hence, we have Listing7.Listing 7:Final join algorithm.
for (i=0 To numberOfWells) [
(1)  Step 2.Add(Step 2.Run(Step 1)
(2)  Step 3.Add(Step 2.Run(oilFiledData, Step 2))
(3)  InitP=InitP(Step 1)
(4)  Step 4.Add(Step 4.Run(Step 3, Step 2, InitP, GasLift)
(5)  Step 5.Add(Step 5.Run(Step 1, Step 2, Step 4, SG, API, WD))
]FunctionInitP is the initial pressure calculation based on WHP i BHP data with the consideration of trajectory curvature.


## 4. Case Study: Simulation

Based on the authors experience in oil and gas industry the small oilfield model was created. The application tool has also been developed. It gives the possibility to look through all the results. We emphasize on the single well, and the full interpretation of simulated properties is not a part of this study. A few examples of the simulated results are presented below.The well trajectory (Figure3) is very characteristic for relatively shallow wells, commonly found in the middle east. Geothermal temperature profile (Figure 4) meets the trajectory. This is the only model, so the surface temperature of about 350 K is relatively high here. This wellbore is a single casing and tubing string with the ID 12 inches and 9 inches, respectively. Total water and oil inflow chart (Figure 5), oil inflow (Figure 6), and gas inflow (Figure 7) charts confirm that two productive layers are given at this oilfield. Total oil production from this simulation is estimated on 2588 BPD which is very accurate value in comparison to real production data for this kind of wells. Gas inflow chart does not consider gas lift from the annulus and it is part of gas accumulated in this reservoir. The most important physical properties in the wellbore are pressure (Figure 9) and temperature (Figure 8). These values meet the criteria very closely but the difference between the simulated temperature and the geothermal temperature is low (Figure 10). These are the results where this simulation is performed under the early stage of production. Below the reservoir the coefficient between temperature and pressure is equal, so it means that this is to prove that the simulation is correct. In Table 1, some results from the simulation are presented.Table 1
Simulation values.
Property
Wellhead value
Bottom hole value
Water heat capacity [Btu/lb*R]
0,982555750110051
0,960586678396243
Water formation volume factor
1,00467452445733
1,02484739136526
Water density [lb/ft2]
8,324
8,35165343126831
Water viscosity [cP]
0,382879865231827
0,25322317446799
Water rate [BPD]
4052,62648
0,24196459483204
MD [ft]
0
3188
TVD [ft]
0
1160
Wellbore temperature [K]
350,150423209288
369,09601487285
Surface tension [dyn/mm]
51,8881443443093
41,5637927992043
Reynolds number
179,714829320058
0,00218243682336947
Wellbore pressure [psi]
1508,01890934205
2000,11967522
Oil viscosity [cp]
1,39027568317487
1,83235372068421
Oil rate [BPD]
2588,381657
0,0145917902247552
Oil density [API]
30
30,0205298689915
Oil heat capacity [Btu/lb*R]
0,442969196896299
0,486093081722976
Oil formation volume factor
1,25236622654504
1,0949263698297
Liquid volume fraction
0,999655206045957
0,999999596771036
Heat transfer [Btu/D/ft2/F]
355147,698582206
207748,085888011
Gas heat capacity [Btu/lb*R]
0,49558307261851
0,543829017247149
Gas formation volume factor
0,010147402881213
0,00924193736098171
Geothermal temperature [K]
350,015683814304
368,130489335027
Gas viscosity [cP]
0,016290648370466
0,0177988938761987
Gas solubility
236,687450637186
257,906526460697
Gas rate [MMScf/D]
1,294190828
0,0000072958951123
Gas density [SG]
0,7
0,700000001749332
Gas compressibility
0,859040738142422
0,908102165721585
Flow type
Segregated
Distributed
Bubble point pressure [psi]
2398,64226965175
11,1431088951275Figure 3
Well trajectory.Figure 4
Geothermal temperature profile.Figure 5
Inflow.Figure 6
Gas inflow per feet.Figure 7
Oil inflow per feet.Figure 8
Wellbore temperature.Figure 9
Wellbore pressure.Figure 10
The difference between wellbore temperature and geothermal temperature.Presented values are the confirmation that the framework and its implementation work properly. These values behave as a real one. Flow rates at the bottom are almost 0, and the surface values are highly reliable: 2588 BPD of oil, 4052 BPD of water, and 1,3 MMScf/D of gas. Temperature and pressure values are highly accurate even if the wellhead temperature is very high. Formation volume factors, heat capacities, and viscosities are accurate as well. The small disadvantage is observed for the Reynolds number at the bottom. This value is almost 0 and it is little bit unexpected according to Moody diagram and its properties. At the surface the Reynolds number gets the proper value in terms of pipe roughness which is equal to 0,0006 inches.


## 5. Conclusions

The literature study, regarding the single property simulations, gives the possibility to create the framework which allows obtaining the simulation for physical properties during the oil production. The authors choose the correlations which meet the criteria in terms of accuracy and efficiency. The obtained results can be classified as proper ones in the comparison to the real data. Unfortunately, upon the legal procedures the comparison results could not be published. But based on the authors experienced in oil and gas industry the presented correlation results are very reliable. The pressure and temperature profiles, as the results of all physical properties simulations, meet the expectations. In this paper, the authors does not analyse physiochemical properties or chemical coordinates, but this analysis is considered as a subject of future paper. This paper has been created to show the possibility of wellbore simulation, which has been proved as an accurate comparison to the real data at oilfields. It is important fact that the whole simulation is running fast because in some cases multithreading technology is used here. The flexibility of the framework idea gives the opportunity to adjust every correlation according to the latest studies. This idea may be extended on optimisation problems for gas lift managing and distribution for the oilfield.



---

*Source: 183982-2015-03-05.xml*