Optimum tubing size prediction model for vertical multiphase flow during flow production period of oil wells

Optimum tubing size (OTS) selection was traditionally done by using nodal analysis to perform sensitivity analysis on the different tubing sizes. This approach was found to be both cumbersome and time-consuming. This study developed a user-friendly and time-efficient OTS prediction computer model that could allow Petroleum Production Engineers to select the best tubing size for any vertical oil well. The tubing size selection was based on the present operating flow rate, economic considerations and future operating flow rate as defined by the OTS selection criteria of this study. The robustness of the model was tested using tubing sizes ranging from 0.824 to 6.0 inch in a vertical well producing from both saturated and undersaturated oil reservoirs. The 2.750-inch tubing was found the OTS for both scenarios. In the validation, the results obtained from the novel OTS prediction model and Guo et al. (Petroleum production engineering: a computer-assisted approach, Gulf Professional Publishing, Cambridge) spreadsheet program using the Poetmann–Carpenter method were in excellent agreement for operating flow rate but not for operating pressure. Furthermore, the novel OTS prediction model was in excellent agreement with the same spreadsheet program based on modified Hagedorn–Brown correlation for both operating flow rate and pressure. The results showed that the model developed in this study is reliable and can be used in the field for vertical oil wells. The new model could as well inform the Production Engineer when the well would need artificial lift for economic production of the well. It was recommended that Newton–Raphson and modified Hagedorn–Brown methods be used in future study.


p wft
Flowing tubing bottom-hole pressure (psia) Reservoir pressure at future time (psia) p wh , p bh Wellhead and bottom-hole pressure (psia) p t (k), p b (k) Pressure at top and bottom of the kth tubing segment (psia) Δp Incremental pressure (psia) (k) Average density of kth tubing segment (lb/ cuft) t (k), b (k) Mixture density at top and bottom of kth tubing segment (lb/cuft) F 2F Fanning's two-phase friction factor f w Water cut D pv Numerator of Reynolds number D(i) Diameter of the ith tubing (inches) M Total mass associated with 1 stb of oil (lb) GOR Producing gas-oil ratio, scf/stb WOR Producing water-oil ratio, bbl/stb GLR Producing gas-liquid ratio, scf/stb 1 3 V mt (k), V mb (k) Mixture volume associated with 1 stb of oil at the top and bottom of kth tubing segment (cuft) B ot (k), B ob (k) Oil formation volume factor at the top and bottom of the kth tubing segment (scf/stb) R st (k), R sb (k) Solution gas-oil ratio at the top and bottom of the kth tubing segment (scf/stb) z t (k), z b (k) Z-factor at the top and bottom of kth tubing segment T wh , T bh Wellhead and bottom-hole temperatures (°F) Temperature at top and bottom of the kth tubing segment (°F) g , o , w Specific gravities of gas, oil and water q Input gas fraction u m Mixture velocity (ft/s) OTS Optimum tubing size h tb Tubing shoe depth (feet) Δh Length of tubing segment (feet)

Introduction
The system analysis that combines the inflow performance relationship (IPR) and tubing performance relationship (TPR) curves so as to obtain the operating flow rate and pressure at a specified node is called nodal analysis (Fan and Sarica 2019). Bottom-hole or wellhead is usually used as the solution node in nodal analysis (Guo et al. 2007). In practice, Jansen (2017) considered IPR to be a near-linear relationship between flow rate, q and bottom-hole flowing pressure, p wf at least when p wf is above the bubble point pressure, p b . The productivity index J is the proportionality constant and is independent of the flow rate above the bubble point pressure (Guo et al. 2007). IPR for single-phase reservoir in this case was defined as the ratio of flow rate, q, to the pressure drawdown, ΔP, (Guo et al. 2007(Guo et al. , 2017Jansen 2017) mathematically expressed as Several empirical equations have been developed in the past for modeling the IPR of two-phase reservoirs (Guo 2019;Guo et al. 2007). IPR curves were numerically computed for wells producing in fictitious solution gas drive reservoirs with a wide range of oil PVT properties and relative permeability characteristics including hydraulically fractured wells (Vogel 1968). Vogel (1968) used the computed IPR curves to develop a reference IPR curve or better known as Vogel's equation [Eq. (2)] applicable to two-phase reservoirs.
(1) J = q P r − P wf Standing (1971) improved upon Vogel's equation and developed an equation that can be used to predict the future IPR of the well given the productivity index measured at a given average pressure as given by Eq. (3).
A series of back pressure test were performed by Fetkovich (1973) for reservoirs above the critical gas saturation and found out that the oil well back pressure curves follow the same general form as that used for expressing the rate-pressure relationship for gas wells as given by Eq. (4).
IMEX and ECLIPSE reservoir simulators were used by Bendakhlia and Aziz (1989) to develop IPR curves for horizontal wells producing from solution gas drive reservoirs. They suggested a new equation in dimensionless form that combines Vogel and Fetkovich equations as follows.
BEST-VHS reservoir simulator was employed by Chang (1992) to numerically develop IPR's of horizontal, vertical and slanted wells for solution gas drive reservoirs. The author also analyzed the newly developed IPR curves and found out that: i. Curvature of the IPR of horizontal wells increases with depletion up to 8% and then decreases slightly at later stages of oil recovery. ii. Curvature of the IPR of vertical wells decreases with depletion.
Daoud et al. (2017) employed a combination of 550 rock and fluid properties to generate 550 IPR curves by performing series of simulations on a single well solution gas drive 3D reservoir. The authors then employed a nonparametric regression technique to generate a new IPR correlation. Their correlation showed a high range of application than traditional correlations because the new IPR correlation is a function of rock and fluid properties that vary from one reservoir to the next.
Several empirical models have also been developed for modeling tubing performance relationship (TPR) for multiphase flow in vertical oil wells. These TPR models can be classified as homogeneous flow and separated flow models (Guo et al. 2007(Guo et al. , 2017Shoham 2006). In their pioneering work for homogeneous flow models, Poettmann and Carpentar (1952) assumed no slip of liquid phases (between oil and water) in developing a gas-oil-water three-phase model that uses an estimate of the average mixture density ( ̄ ) and two-phase friction factor ( f 2F ) in computing the pressure losses in vertical wellbores. The acceleration term in their work was neglected, and their equation for calculating the pressure traverse in vertical tubing in field units was given as: Poettmann and Carpentar (1952) correlation was improved upon by Cicchitti et al. (1960) and Dukler et al. (1964) by including the effect of viscosity in the homogeneous model. Furthermore, Guo and Ghalambor (2005) improved on the work of Poettmann and Carpentar (1952) by developing a gas-oil-water-sand four-phase model. Iterations were required to solve the Poettmann-Carpenter model as well as the Guo-Ghalambor model; hence, computer spreadsheet programs were developed by Guo et al. (2007) to this effect.
Several separated flow models have been developed over the years. Griffith and Willis (1961) developed a separated flow model for slug flow regime. Duns and Ros (1963) developed a correlation with its associated calculating procedure for mist flow. This correlation was applicable in gas/ condensate wells as well as in gas/condensate wells with water cut as long as no emulsion is formed. Furthermore, Hagedorn and Brown (1965) developed correlations and equations that permit accurate prediction of flowing pressure gradients for a wide range of flow conditions and liquid properties. The modified Hagedorn-Brown (mH-B) method was developed based on the original work of Hagedorn and Brown (1965). The mH-B method uses either the original Hagedorn and Brown (1965) correlation or the Griffith and Willis (1961) correlation in calculating the liquid holdup depending on whether the predicted flow regime is bubble flow regime or not. Bubble flow regime exists when the input gas fraction (also referred to as no-slip holdup of gas) q is less than a quantity, L B as given by Eq. (7).
where u m , mixture velocity and D, tubing internal diameter.
The quantity L B must be greater than or equal to 0.13 and as such L B is set at 0.13 if the calculated value of L B is less than 0.13 (Economides et al. 1994;Guo et al. 2007Guo et al. , 2017. From the above explanation, it was suggested that if Δh 144 bubble flow regime is predicted by g < L B , then the Griffith and Willis (1961) correlation is used to determine the liquid holdup, otherwise the original Hagedorn and Brown (1965) correlation is used to determine the liquid holdup. Once a well has been drilled and completed, the reservoir fluid can be produced through the casing, tubing or both. Mostly, wells are produced through the tubing in order to isolate the casing from corrosion and for use of artificial lift system (Guo et al. 2007;Guo 2019). It became essential that for any selected production tubing size, the well should flow naturally. Choosing an undersized tubing will result in excessive flow velocity and hence increased friction resistance in the well which limits the well production rate. The undersized tubing may as well restrict the type and size of artificial lift equipment (Clegg 2007). Using oversized tubing on the other hand would result in low flow velocity and hence excessive liquid loss due to gas slippage effect. Large tubing size would also complicate workover operations due to loading of the well resulting from heading and unstable flow (Clegg 2007). Renpu (2011) emphasized the need to select the optimum tubing size that ensures an optimum state for the friction resistance and liquid-phase loss due to gas slippage which in turn would ensure the longest flowing time and lowest lifting energy consumption.
Optimum tubing selection (OTS) is commonly done by using nodal analysis to perform sensitivity analysis on the various tubing sizes during the flowing production period of the well (Renpu 2011). The author further gave two traditional methods that can be used in analyzing the optimum tubing size in which the OTS is taken as that tubing size that will: i. Maximize the production rate for a given surface condition such as wellhead pressure, separator pressure or flowline size. ii. Minimize the producing gas-oil ratio, maximize gas expansion energy utilization efficiency and ensure the longest production time for a given production rate.
For a given wellhead pressure, the sensitivity analysis requires the determination of the IPR curve. Then, any of the TPR models can be used to determine the TPR for all the various tubing sizes. The Poettmann and Carpenter TPR model was selected for use in this work. The point of intersection of the IPR and TPR curves gives the operating point of the various tubing sizes. The tubing size with the maximum flow rate becomes the optimum tubing size.
C# programming language was employed by Li and Xiong (2011) to develop a software for optimizing tubing and production casing sizes for flowing oil and gas wells. Their software was shown to produce the accurate results but does not consider the future IPR and hence would not be able to predict when the selected optimum tubing size will stop flowing.
Furthermore, several total pressure gradient (TPG) models were generated as a function of liquid density, gas density, liquid flow rate, gas flow rate and tubing size after fitting 100 data points from the Niger Delta (Kinate et al. 2018). The authors then selected the best TPG model and used it with the aid of Microsoft Visual Basic to develop an optimum tubing size (OTS) prediction software. Their optimum tubing size prediction software uses the total pressure gradient in selecting the optimum size. The model showed better performance than traditional TPG models. A major drawback of the model is inability to predict when the selected OTS will stop flowing. Hence, time to consider any of the artificial lift methods cannot be predicted.
The modified Hagedorn and Brown (1965) correlation was recently used by Santhosh (2017) to determine TPR curves for different tubing sizes. The author's procedure involved cumbersome and time-consuming sensitivity analysis on the different tubing sizes. This was exacerbated by the fact the author also had to use charts to obtain the liquid holdup and friction factor. The optimum tubing size was selected by just observing the shape of the different TPR curves. Their procedure again was not designed to predict when the well will stop flowing and hence would not be able to predict when to consider any of the artificial lift methods. Agboola et al. (2016) used MS Excel and Visual Basic to develop a simple and user-friendly pressure gradient model that can be used in selecting the optimum tubing size. Their model is designed to predict when the well will stop flowing and hence the time to consider any of the artificial lift methods.
Sensitivity analysis of the tubing size has been observed to be time-consuming and cumbersome. The objective of this work is to use MATLAB in developing a user-friendly and time-efficient computer model for optimum tubing size (OTS) prediction. The model would generate and plot the present and future IPR curves for the well and TPR curves for different tubing sizes on same axis. The optimum tubing size selection criteria would be the basis for selecting the OTS. The model developed in this work would further inform the Production Engineer if the predicted OTS would permit natural flow when reservoir pressure drops to 75% and 50% of its original value, respectively. The novel OTS prediction has similar capabilities as the model developed by Agboola et al. (2016) since both models are capable of predicting the OTS, time the well will stop flowing and time to consider any of the artificial lift methods. Spreadsheet programs have been developed by Guo et al. (2007) for determining the operating point of flowing oil wells. The spreadsheet programs entitled "BottomholeNodalOil-PC. xls" and "BottomholeNodalOil-HB.xls" that are based on the Poettmann-Carpenter and modified Hagedorn-Brown models, respectively, have been used in this study to validate our novel OTS prediction model.

Saturated oil reservoir
In saturated oil reservoir, the reservoir pressure p r is usually less than or equal to the bubble point pressure p b (that is, p r ≤ p b ). Under this scenario, the production flow rate can be determined using Eq. (8); where q max = Jp r ∕1.8 is the maximum flow rate for saturated reservoir scenario as initially derived by Standing (1981).
The values of the flow rate were generated with Eq. (10) and used in Eq. (9) to determine corresponding values of the flowing bottom-hole pressure and the result used to generate the IPR curve for a saturated oil reservoir (Guo 2019). The flow rate data generated would run between zero and maximum flow rate as given in Eq. (10).
where 1 ≤ j ≤ m + 1 i and m are positive integers.

Undersaturated oil reservoir
It has been understood that reservoir pressure continues to decline as oil production progresses. An oil reservoir above bubble point pressure is usually undersaturated, but this does not remain constant in the life of an oil reservoir because as production continues, the reservoir pressure continues to decline and could get to the bubble point pressure or even below it. Guo et al. (2007) expressed a situation of partial twophase oil reservoir, where the flowing bottom-hole pressure below the bubble point pressure still has the reservoir pressure above the bubble point pressure. They proposed a generalized IPR model that combines the straight-line IPR model for single-phase flow with Vogel's IPR model for twophase flow as given in Eq. (11); The maximum flow rate occurs when the flowing bottomhole pressure equals zero, and substituting in Eq. (11) gives a maximum flow rate given by Eq. (12) as the sum of q b and q v expressed as; Equations (8) and (11) can both be rewritten in dimensionless form as given by Eqs. (13) and (14), respectively. Performing an analogy between Eqs. (13) and (14) showed that both equations are similar if q(j) q max and p r in Eq. (13) are, respectively, replaced with (14).
It follows then that replacing q(j) q max and p r in Eq. (9) with and p b , respectively, results in Eq. (15), the Generalized Vogel's equation version of Eq. (9). Equation (15) was used in this work in generating the IPR curve for an undersaturated reservoir with flow rate data generated between zero and the maximum flow rate as given in Eq. (10).

Tubing performance relationship (TPR) generation
This study adopted the Poettmann and Carpentar (1952) model for tubing performance relationship (TPR) since the model considers gas-oil-water three-phase homogeneous reservoir fluid flow. Poettmann and Carpentar (1952) model appears to be the earliest and widely accepted empirical correlation for TPR. This model considered real field stream, though less accurate, but its mechanistic nature gives the model popularity with field calibrations. Being mechanistic allows models like Poettmann and Carpenter to cover all the parameter ranges in field operations utilizing the basic physics of multiphase flow and widen the ranges of application of the model. Guo et al. (2017) preferred Poettmann andCarpentar (1952) model for ease of coding in a computer program. The pressure traverse with fixed length calculation procedure as proposed in Economides et al. (1994) was employed in this study where the tubing was divided into segments of equal lengths, Δh with the wellhead as the top of the first tubing segment and the bottom-hole as the bottom of the last tubing segment. Equations (16) and (17) were applied to each of the segments starting with the first tubing segment with the wellhead pressure as the input pressure.
where 1 ≤ k ≤ ns ; p b (k) , pressure at bottom of the kth tubing segment of length Δh ; p t (k) , pressure at top of the kth tubing segment of length Δh ; ̄(k) , average density of the kth tubing segment of length Δh ; and ns , number of tubing segments of length Δh.
Equations (16) and (17) clearly showed that the Poettmann and Carpentar (1952) model uses the average mixture density ̄ and two-phase friction factor f 2F to determine the pressure p b (k) at the bottom of the tubing segment. The average density, ̄ , and two-phase friction factor, f 2F , would be the input variables into Eqs. (16) and (17).

Two-phase friction factor determination
The two-phase friction factor, f 2F , can be determined from the charts presented by Poettmann and Carpentar (1952). However, for ease of coding in the computer, the correlation developed by Guo and Ghalambor (2002) was employed since it approximates the values of the Poettmann and Carpentar (1952) two-phase friction factor chart given by Eqs. (18) and (19). The two-phase friction factor is a function of flow rate and tubing size; hence, it is same for all tubing segments for a particular flow rate and tubing size.

Average density determination
According to the Poettmann and Carpentar (1952) model, gas-oil-water three-phase mixture in the tubing was treated as homogeneous mixture. In order to simulate this homogeneous mixture, it became necessary that the mixture density at the top and bottom of each tubing segment be determined and then used to determine the average mixture density of each tubing segment. This was achieved by computing the density at the top and bottom of the kth tubing segment based on the mass flow rate, M , and volume flow rate at the top, V mt (k) , and bottom, V mb (k) , of the kth tubing (Guo 2019;Nind 1964). Total mass associated with one stock tank barrel of oil, M (Guo 2019) would be computed first given by Eq. (20).
The gas-oil ratio (GOR), water-oil ratio (WOR) and oilspecific gravity ( o ) are as given in Eqs. (21), (22) and (23), respectively. These quantities would be used in Eq. (20) for the determination of the total mass associated with one stock tank barrel of oil, M.

Methodology
Production facilities used in oil and gas fields to a large extent is an important factor in both economic and technical evaluations of oil production from a well. Imposed pressure at the wellhead through choke size determines the bottomhole pressure in a producing well. In effect, the achievable production rate from a well according to Guo et al. (2007) is determined through wellhead pressure and flow performance of production string (i.e., tubing, casing or both). Proper evaluation of flow performance of a well must establish a mathematical relationship that involves tubing geometry (size), wellhead pressure imposed by the choke size, bottom-hole pressure, fluid (oil, gas and water) properties and production rate. Good knowledge and understanding of this process in the industry could lead to design of production facilities for optimum well production.

Model development
The model in this study would be developed such that it can handle data for both undersaturated and saturated oil reservoirs. This justifies the reason to use the developed model in this study at all stages of the production life of the well no matter the pressure value of the reservoir.
A computer model for optimum tubing size prediction was developed with MATLAB. The model would aid the Production Engineer in selecting the optimum tubing size that would permit natural flow of reservoir fluid to the surface at present and future reservoir pressures. This study adopted Vogel (1968) equation since it has been acclaimed as the most widely used in oil industry and the Poettmann and Carpentar (1952) gas-oil-water three-phase TPR model.
The following assumptions were made in the development of the computer model;

Initialization of tubing segments
For the average density determination, each of the tubing segments must be initialized. The first tubing segment (k = 1) is initialized by setting the pressure, p t (1) and temperature, T t (1) , at the top of the first tubing segment equal to the wellhead pressure, p wh , and temperature, T wh , respectively; The second to the last tubing segments (k = 2 to k = ns) were each initialized by setting the pressure, p t (k) , and temperature, T t (k) , at the top of the current tubing segment equal to the pressure, p b (k − 1) , and temperature, T b (k − 1) , at the bottom of the previous tubing segment, respectively. Consider a case in which the second tubing segment is the current tubing segment and the first tubing segment is the previous tubing segment. It follows from our augment that the pressure, p t (2) , at the top of the second tubing segment is equal to the pressure, p b (1) at the bottom of the first tubing segment. Similarly, the temperature, T t (2) at the top of the second tubing segment is equal to the temperature, T b (1) , at the bottom of the first tubing segment; p t (1) = p wh T t (1) = T wh A general equation relating the pressure, p t (k) , and temperature, T t (k) , at the top of the current tubing segment to the pressure, p b (k − 1) , and temperature, T b (k − 1) , at the bottom of the previous tubing segment for tubing segments k = 2 to k = ns would be given as in Eqs. (24) and (25), respectively.
where p t (k) , pressure at the top of the current tubing segment, psia; p b (k − 1) , pressure at the bottom of the previous tubing segment, psia; T t (k) , temperature at the top of the current tubing segment, °F; T b (k − 1) , temperature at the bottom of the previous tubing segment, °F.

Solution procedure for all tubing segments
For tubing segments k = 1 to k = ns, the pressure, p t (k) , and temperature, T t (k) , at the top of the kth tubing segment were used in computing the solution gas-oil ratio at the top of the kth tubing segment, R st (k) by employing the Standing (1981) solution gas-oil ratio correlation [Eq. (26)] as reported in Guo et al. (2017). Furthermore, R st (k) and T t (k) were applied in the Standing (1981) formation volume factor correlation [Eq. (27)] as reported in Guo et al. (2017) in approximating the oil formation volume factor at the top of the kth tubing segment, B ot (k).
In addition, the z-factor, z t (k) , at the top of the kth tubing segment was computed by employing the Brill and Beggs (1974) z-factor correlation with p t (k) and T t (k) as input data. This correlation was considered one of the best explicit z-factor correlations and requires a non-iterative procedure for convergence (Kareem et al. 2015) unlike Hall and Yarborough (1973) z-factor correlation and has been confirmed to be accurate enough devoid of the Newton-Raphson method by Guo et al. (2007) for many engineering calculations.
The mixture volume, V mt (k) , associated with one stock tank barrel of oil at the top of the kth tubing segment was determined by employing Eq. (28) as initially reported by Poettmann and Carpentar (1952) and later reported in a simplified form by Guo (2019). The mixture density at the top of the kth tubing segment, (29)], was calculated based on the total mass, M and mixture volume associated with one stock tank barrel of oil at the top of the kth tubing segment, V mt (k) , as reported by Poettmann and Carpentar (1952) and Guo (2019).
The computed mixture density, t (k) , and pressure, p t (k) , at the top of the kth tubing segment were used in Eq. (17) to obtain an approximate value of the pressure at the bottom of the kth tubing segment as given in Eq. (30). Economides et al. (1994) reported that the temperature profile between the wellhead temperature, T wh , and bottomhole temperature, T bh , is approximately linear. Therefore, the use of equation of a straight line to obtain the value of the temperature at the bottom of the kth tubing segment, T b (k) as given in Eq. (31) is valid.
The Standing (1981) solution gas-oil ratio correlation [Eq. (32)] as reported in Guo et al. (2007) was also used for solution gas-oil ratio approximation at the bottom of the kth tubing segment, R sb (k) with p b (k) and T b (k) as input data. The difference between R st (k) and R sb (k) gives R st (k) as a function of pressure, p t (k) , and temperature, T t (k) , at the top of the kth tubing segment while R sb (k) as a function of the pressure, p b (k) , and temperature, T b (k) , at the bottom of the kth tubing segment. The oil API and gas-specific gravity g remain unchanged in Eqs. (26) and (32). Similarly, the oil formation volume factor at the bottom of the kth tubing segment, B ob (k) was approximated using Standing (1981) formation volume factor correlation [Eq. (33)] as reported in Guo et al. (2017). The difference between B ot (k) and B ob (k) is that B ot (k) is a function of the solution gas-oil ratio, R st (k) , and temperature, T t (k) , at the top of the kth tubing segment while B ob (k) is a function of the solution gas-oil ratio, R sb (k) , and temperature, T b (k) , at the bottom of the kth tubing segment. The specific gravities
The Brill and Beggs (1974) z-factor correlation was equally used to determine the z-factor, z b (k) , at the bottom of the kth tubing segment with p b (k) and T b (k) as input data.
The mixture volume, V mb (k) , associated with one stock tank barrel of oil at the bottom of the kth tubing segment was computed using Eq. (34) as initially reported by Poettmann and Carpentar (1952) and later reported in a simplified form by Guo (2019). The difference between V mt (k) and V mb (k) is that V mt (k) is a function of solution gas-oil ratio, oil formation volume factor, z-factor, temperature and pressure at the top of the kth tubing segment, while V mb (k) is a function of solution gas-oil ratio, oil formation volume factor, z-factor, temperature and pressure at the bottom of the number kth tubing segment. The water-oil ratio, gas-oil ratio and water formation volume factor remain the same in both Eqs. (28) and (34).
The mixture density at the bottom of the kth tubing segment was computed as the total mass associated with one stock tank barrel of oil, M divided by mixture volume associated with one stock tank barrel of oil at the bottom of the kth tubing segment, V mb (k) as given in Eq. (35) (Guo et al. 2007(Guo et al. , 2017Poettmann and Carpentar 1952).
Gas is compressible and that is why the volume of a given mass of gas and its density changes significantly with pressure and temperature. The z-factors, z t (k) and z b (k) were, respectively, incorporated into V mt (k) and V mb (k) to account for these changes in gas volume and gas density with pressure and temperature in our model. The average mixture density, ̌(k) , of the kth tubing segment was computed as the average of the mixture density at the top, t (k) and bottom, b (k) of the kth tubing segment thereby simulating the assumption of homogeneous mixture (Guo et al. 2007(Guo et al. , 2017Poettmann and Carpentar 1952).
The computed average mixture density, ̌(k) , of the kth tubing segment was used in Eq. (17) to compute the correct value of the pressure, p b (k) , at the bottom of the kth tubing segment as given by Eq. (37).
The pressure at the bottom of the last tubing segment (k = ns) was recorded as the flowing tubing bottom-hole pressure, p wft (j) at flow rate, q(j) . That is; The solution procedure is repeated for the next flow rate data from j = 2 to j = m + 1 , and this results in the TPR curve for a particular tubing size, D(i), i = 1 . The entire process is repeated for the next tubing diameter data, D(i) , from i = 2 to i = n.

Optimum tubing size selection criteria
The selection criteria for optimum tubing size (OTS) were based on the fact that the operating flow rate will increase with increase in tubing size. However, a certain tubing size will be reached where this increase of operating flow rate with increase in tubing size becomes negligible. The tubing size just before this negligible increase in operating flow rate should be considered as the optimum tubing size provided it satisfies all three OTS selection criteria. If it does not satisfy all three selection criteria, the immediate larger tubing size should be considered if it satisfies three OTS selection criteria. This should be repeated until a tubing size is reached that satisfies all three OTS selection criteria. The three OTS selection criteria are as given below: 1. The difference between the operating flow rate of the considered OTS and the immediate larger tubing size must be minimal. 2. The considered OTS must be cheaper than the larger tubing sizes. 3. The considered OTS must be able to produce when the reservoir pressure drops to 75% of its original value.
The flowchart for optimum tubing size prediction model during the flow production period of an oil well is as shown in Fig. 1.

Model test
Data obtained from Guo et al. (2007Guo et al. ( , 2017 as summarized in Tables 1 and 2 were used to test the robustness 1 3 of the novel optimum tubing size prediction model. Internal diameter of production tubing used ranged from 0.824 inch to 6.0 inches. The number of tubing segments, ns in Table 2, was selected such that the incremental length, Δh , of the tubing segments was approximately equal to 200 feet which is a typical length increment value for flow in tubing as reported by Economides et al. (1994).

Results
See Figs. 2, 3, 4, 5, 6 and 7.  Figure 2 shows the present and two future IPR curves of a saturated oil reservoir and TPR curves for different tubing sizes plotted on the same axis. Figure 3, the zoomed copy of Fig. 2, shows that the operating rate increases with increase in tubing size up to the 2.750 inch after which the increase in operating rate became negligible. This became evident from the clustering of the TPR curves from 2.750 inch TPR curve and above. Based on the OTS selection criteria, the 2.750inch tubing was selected as the optimum tubing size because it satisfied all three OTS selection criteria depicted below:

Saturated oil reservoir example
i. The 2.750-inch tubing would permit natural flow to the surface at a rate of 1250 stb/day as compared to the largest tubing size of 6 inch with an operating rate of 1280 stb/day as shown in Fig. 4. The small difference in rate of 30 stb/day does not justify using a 6-inch tubing. ii. From economic point of view, the 2.750-inch tubing was cheaper than any of the larger tubing sizes. iii. The 2.750-inch tubing still permits natural flow at a rate of 540 stb/day when the reservoir pressure drops to 75% of its original value (Fig. 4). This shows that the OTS will permit flow for a long time.
The model developed in this study was designed to check if the selected optimum tubing size will permit flow naturally to the surface when the reservoir pressure drops to 75% and 50% of its original value, respectively. As shown in Fig. 4, when the reservoir pressure drops to 75% of its original value the OTS of 2.75 inch was still flowing naturally at a flow rate of 540 stb/day. However, the OTS of 2.750 inch stopped flowing when the reservoir pressure dropped to 50% of its original value. This would serve as a valuable source of information to the Production Engineer to consider any of the artificial lift methods when the reservoir pressure starts dropping close to 50% of its original value. Figure 5 shows the present and two future IPR curves for an undersaturated oil reservoir as well as TPR curves for different tubing sizes plotted on the same axis. Figure 6, the zoomed version of Fig. 5, shows that the operating rate increases with increase in tubing size, but the increase in operating rate became negligible after the tubing size of 1.867 inch was reached. This is evident in the clustering of the TPR curves starting from the TPR curve of 1.867-inch tubing. Based on the proposed OTS selection criteria, the 1.867-2.441-inch tubing sizes were not considered because they satisfied only two out of the three OTS selection criteria. The 2.750-inch tubing was selected as the optimum tubing size because it satisfied all three OTS selection criteria as depicted below:

Undersaturated oil reservoir example
i. It produced at a rate of 2280 stb/day as compared to the other larger tubing sizes with the 6.0 inch being the largest with operating rate of 2320 stb/day as shown in Fig. 7. ii. It is cheaper than the other larger tubing sizes. iii. It will flow naturally at a rate of 250 stb/day when reservoir pressure drops to 75% its original value (Fig. 7).
As shown in Fig. 7, the selected OTS of 2.750 inch would no longer support natural flow as soon as the reservoir pressure drops below 75% its original value. This served as a

Model validation
The results provided by the OTS prediction model were compared with those of the excel spreadsheets program presented by Guo et al. (2017). Two examples as originally presented in Guo et al. (2017) and summarized in Table 1 were used in the validation of the model.  The BottomholeNodalOil-PC.xls spreadsheet program that is based on the Poetmann-Carpentar method was used by Guo et al. (2017) in solving the saturated oil reservoir example (Table 1) for the operating point of a 1.660-inch tubing. An operating flow rate and pressure of 1127 stb/day and 1873 psia, respectively, were obtained for the 1.660-inch tubing. The new model was equally used to solve the same saturated oil reservoir example, and the result is as shown in Fig. 8. As shown in Fig. 8, operating flow rate and pressure of 1160 stb/day and 1500 psia, respectively, were obtained for the 1.660-inch tubing.
The BottomholeNodalOil-HB.xls spreadsheet program that is based on the modified Hagedorn-Brown correlation was also used by Guo et al. (2017) to solve the undersaturated oil reservoir example (Table 1) for the operating point of a 1.995-inch tubing. The spreadsheet program gave operating rate and pressure of 2200 stb/day and 3500 psia, respectively. The new model was also used to     where q er , operating flow rate error percent; p er , operating pressure error percent; q lit , operating flow rate obtained from spreadsheet program; p lit , operating pressure obtained from spreadsheet program; q mod , operating flow rate obtained from novel OTS prediction model; and p mod , operating pressure obtained from novel OTS prediction model. As shown in Table 3, the results provided by the novel OTS prediction model and the BottomholeNodalOil-PC. xls spreadsheet program based on the Poetmann-Carpenter method were in excellent agreement for the operating flow rate but not for the operating pressure. Furthermore, the novel OTS prediction model was in excellent agreement with the BottomholeNodalOil-HB.xls spreadsheet program that was based on modified Hagedorn-Brown correlation for both operating flow rate and pressure.

3
i. The increase in operating rate between the OTS and the next larger tubing size was minimal. ii. The selected OTS is cheaper than the other larger tubing sizes. iii. The selected OTS would still produce naturally when reservoir pressure drops to 75% its original value.
The results showed that the computer model developed in this study for optimum tubing size prediction was reliable and could be applied to both saturated and undersaturated oil reservoirs. The incorporation of the future IPR in the model developed would serve as a valuable source of information to the Production Engineer on when to consider any of the artificial lift methods.
The results provided by the novel OTS prediction model and the BottomholeNodalOil-PC.xls spreadsheet program based on the Poetmann-Carpenter method were in excellent agreement for the operating flow rate but not for the operating pressure. Furthermore, the novel OTS prediction model was in excellent agreement with the Bottomhole-NodalOil-HB.xls spreadsheet program that was based on modified Hagedorn-Brown correlation for both operating flow rate and pressure.

Recommendations
It was recommended that further research is carried out using the Newton-Raphson method to iteratively determine the operating points and hence the optimum tubing size. Furthermore, a separated flow model such as the modified Hagedorn-Brown method should also be used in place of the Poettmann-Carpenter homogeneous flow model. The results of the OTS model to be developed in the future study should be compared with the current novel OTS model and the BottomholeNodalOil-PC.xls spreadsheet program based on the Poetmann-Carpenter method. This would further investigate why the novel OTS prediction model and the BottomholeNodalOil-PC.xls spreadsheet program based on the Poetmann-Carpenter method were not in excellent agreement for operating pressure.