Two-Phase Flow in Pipes: Numerical Improvements and Qualitative Analysis for a Refining Process

Two-phase flow in pipes occurs frequently in refineries, oil and gas production facilities and petrochemical units.The accurate design of such processing plants requires that numerical algorithms be combined with suitable models for predicting expected pressure drops. In performing such calculations, pressure gradientsmaybeobtained fromempirical correlations suchasBeggs andBrill, and theymust be integrated over the total length of the pipe segment, simultaneously with the enthalpy-gradient equation when the temperature profile is unknown. This paper proposes that the set of differential and algebraic equations involved should be solved as aDifferential AlgebraicEquations (DAE)System,which poses a more CPU-efficient alternative to the “marching algorithm” employed by most related work. Demonstrating the use of specific regularization functions in preventing convergence failure in calculations due to discontinuities inherent to such empirical correlations is also a key feature of this study.The developed numerical techniques are then employed to examine the sensitivity to heat-transfer parameters of the results obtained for a typical refinery two-phase flow design problem. Résumé — Écoulements diphasiques dans les conduites : améliorations numériques et analyse qualitative pour un procédé de raffinage — Les écoulements diphasiques apparaissent fréquemment dans les conduites des raffineries, des installations de production de pétrole et de gaz et des unités pétrochimiques. La conception précise de telles unités requière l’intégration d’algorithmes numériques dans des modèles adaptés, afin de prédire les baisses de pression. Dans le cadre de tels calculs, les gradients de pression peuvent être obtenus à partir de corrélations empiriques, telles que celles de Beggs et Brill, et doivent être intégrés sur la longueur totale du segment de la conduite, en utilisant simultanément l’équation du gradient d’enthalpie lorsque le profil de température n’est pas connu. Cet article propose que les systèmes d’équations différentielles et algébriques en question soient résolus sous forme d’un système d’équations différentielles algébriques (DAE, Differential Algebraic Equations) offrant une alternative économe en temps calcul (CPU) aux algorithmes de type « marching algorithm » utilisé dans la plupart des travaux connexes. Nous démontrons aussi que l’utilisation de fonctions de régularisation spécifiques permet d’éviter les erreurs de convergence dans les calculs, imputables aux discontinuités inhérentes à de telles corrélations empiriques. Les techniques numériques développées sont alors utilisées pour examiner la sensibilité des résultats obtenus aux paramètres de transfert de chaleur pour un problème typique de conception d’écoulement diphasique dans une raffinerie. Oil & Gas Science and Technology – Rev. IFP Energies nouvelles, Vol. 70 (2015), No. 3, pp. 497-510 R.G.D. Teixeira et al., published by IFP Energies nouvelles, 2014 DOI: 10.2516/ogst/2013191 This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Re´sume´-E ´coulements diphasiques dans les conduites : ame´liorations nume´riques et analyse qualitative pour un proce´de´de raffinage -Les e´coulements diphasiques apparaissent fre´quemment dans les conduites des raffineries, des installations de production de pe´trole et de gaz et des unite´s pe´trochimiques.La conception pre´cise de telles unite´s requie`re l'inte´gration d'algorithmes nume´riques dans des mode`les adapte´s, afin de pre´dire les baisses de pression.Dans le cadre de tels calculs, les gradients de pression peuvent eˆtre obtenus a`partir de corre´lations empiriques, telles que celles de Beggs et Brill, et doivent eˆtre inte´gre´s sur la longueur totale du segment de la conduite, en utilisant simultane´ment l'e´quation du gradient d'enthalpie lorsque le profil de tempe´rature n'est pas connu.Cet article propose que les syste`mes d'e´quations diffe´rentielles et alge´briques en question soient re´solus sous forme d'un syste`me d'e´quations diffe´rentielles alge´briques (DAE, Differential Algebraic Equations) offrant une alternative e´conome en temps calcul (CPU) aux algorithmes de type « marching algorithm » utilise´dans la plupart des travaux connexes.Nous de´montrons aussi que l'utilisation de fonctions de re´gularisation spe´cifiques permet d'e´viter les erreurs de convergence dans les calculs, imputables aux discontinuite´s inhe´rentes a`de telles corre´lations empiriques.Les techniques nume´riques de´veloppe´es sont alors utilise´es pour examiner la sensibilite´des re´sultats obtenus aux parame`tres de transfert de chaleur pour un proble`me typique de conception d'e´coulement diphasique dans une raffinerie.No-slip liquid holdup r L Surface tension q Density q L Liquid phase's density q V Vapor phase's density q n Two-phase mixture's density (k L as weighting factor) q s Two-phase mixture's density (HL as weighting factor) l L Liquid phase's absolute viscosity l V Vapor phase's absolute viscosity l n Two-phase mixture's absolute viscosity (k L as weighting factor) h

LIST OF SYMBOLS
Pipe's inclination angle from horizontal

INTRODUCTION
Simultaneous flow of liquid and vapor phases through pipes takes place in countless situations in refineries, oil and gas production facilities and petrochemical units.Therefore, design of such petroleum processing plants requires the ability to estimate liquid content (holdup), pressure drops and temperature changes in these pipelines with acceptable engineering accuracy.These calculations must be performed through application of principles of conservation of mass, momentum and energy, which can be most rigorously accomplished through Computational Fluid Dynamics (CFD) simulations.While this may be the most logical choice in theory and for research purposes, it poses practical problems when design of industrial facilities is concerned, due to the large numbers of pipes involved.It is not uncommon for refining units to contain thousands of pipes.Even when two-phase flow is only expected in a small fraction of these, it still is impractical to run CFD simulations for each one of them.
In a typical refinery unit design, sizing from a pressure loss standpoint is the sole requirement for two-phase flow pipes.Only in specific situations (such as severe operating conditions) does the need for rigorous local temperature and pressure results arise.Thus, it has become common practice in industry to employ CFD tools only when particularly critical conditions or requirements must be met.When this is not the case, two-phase flow pressure-gradient correlations are preferred, since these can be more readily solved for given diameters, vapor fractions and physical properties.
Pressure (and possibly, temperature) varies with distance along the pipe length.As a result, vapor fraction and physical properties also change along the pipe, and so does the pressure gradient.Hence, pressure drop calculations in two-phase flow require that the pressure gradient be numerically integrated from one end of the pipe to the other [1].It has been found that obsolete, CPU-inefficient (i.e., computer costly) numerical integration techniques are still widely used for this purpose, in spite of the advances that have been made in this area in the past decades.Given this context, this study aims to demonstrate the application of a numerical integration method which is shown to be better suited for two-phase flow problems.
It is well known that some two-phase flow correlations exhibit large discontinuities (usually when flow conditions change just enough that a different flow pattern gets predicted), and that these discontinuities may lead to convergence problems when integrating the pressure gradient over the pipe length [1][2][3].The use of specific regularization functions in preventing convergence failure in these cases is also addressed here.
Oil and gas production problems have motivated most of the work done on the subject of two-phase flow in pipes so far, including development of pressure gradient correlations.Bibliographical references which prioritize this type of problem are abundant, and can be easily found [1,4].On the other hand, process engineers involved in simulation and design of refineries and petrochemical process units are faced with a lack of studies on two-phase flow in pipes when it comes to petroleum fractions and petroleum downstream industry.As a result, correlations originally developed for oil and gas twophase flow have also been widely adopted in downstream design calculations.
Given the above scenario, all calculations performed in this work consider a typical refinery vapor-liquid mixture.Obtained results are analyzed in a first attempt to determine which heat-transfer considerations are critical (and which are not) in this petroleum downstream industry design problem.In spite of this, it must be emphasized that the numerical techniques presented here are relevant for oil and gas production calculations as well.

PRESSURE-GRADIENT CORRELATIONS
The first correlations ever developed for prediction of pressure gradients in two-phase flow in pipes are empirical in nature, often based on dimensionless groups and data from laboratory test facilities (with a few using field data) [5][6][7].
It has long been recognized that the accuracy of pressure-gradient predictions cannot be improved without the introduction of basic physical mechanisms, which is the main limiting factor of the empirical methods [1].For this reason, over the past decades, efforts have been directed to the development of phenomenological or mechanistic correlations, i.e., correlations in which fundamental laws are used to model the most important flow phenomena and less empiricism is required [3,8].
Empirical correlations were readily adopted in the petroleum industry, as they were the first to be conceived.Although greater accuracy is expected when using a theoretical correlation instead of an empirical one, selecting a mechanistic method for use in design calculations is not a trivial task, specially when it comes to refining or petrochemical processes, for some of these methods only apply to certain pipe inclinations [8] while others [3] have not been validated enough to persuade engineers and companies into updating long-experienced design practices.It follows that empirical methods still enjoy widespread use in industry and also for comparison purposes in the development of new predictive techniques [9][10][11].In particular, the correlation of Beggs and Brill [5] stands out for being able to predict pressure drops for angles other than vertical upward flow.Interestingly enough, the pressure drops calculated by this method in recent studies have been found to show acceptable agreement with corresponding CFD results [12,13].
The Beggs and Brill correlation was adopted for liquid holdup and pressure gradient calculations in this work for the reasons outlined above, and also because it shows the discontinuities that can lead to convergence failure, which this study intends to demonstrate how to avoid.A brief summary of this method is presented in this section, as some of its equations and coefficients need to be referenced in subsequent discussions.

The Correlation of Beggs and Brill
The velocities at which both phases would flow individually through the pipe cross section (superficial velocities) can be calculated from the internal diameter of the pipe (d) and their volumetric flow rates (q L and q V ): In which A t ¼ pd 2 =4 is the pipe cross-sectional area.
A mixture velocity (v m ) and corresponding Froude number (N Fr ) are given by: If the liquid and vapor phases both travelled at the same velocity (no-slippage), the liquid holdup along the entire pipeline would equal the input liquid content.Therefore, no-slip liquid holdup (k L ) can be calculated from: The flow pattern which would exist if the pipe were horizontal can be determined by the location of the pair (k L , N Fr ) in the correlation flow pattern map.The chart area is divided in four regions -each corresponding to one of four possible flow regimes: Segregated, Transition, Intermittent and Distributed -by the following transition boundaries: Thus, each flow pattern is associated with a set of inequalities, as shown in Table 1.
In two-phase flow in pipes, liquid and vapor phases travelling at different velocities (slippage) result in in-situ liquid volume fractions which differ from the original liquid input content.Liquid holdup is defined as the fraction of the volume of a segment of the pipeline which is occupied by liquid.
The holdup which would exist at the given conditions if the pipe were horizontal can be determined from: with the restriction that: The coefficients a, b and c are given in Table 2 for all flow patterns, except for Transition, which is discussed afterwards.
The horizontal-position liquid holdup must be corrected for the effect of pipe inclination by a multiplicative factor w: This correction factor is calculated from the liquid velocity number (N Lv ), the pipe inclination angle from horizontal (hÞ; k L and N Fr : In which q L and r L are the density and the surface tension of the liquid phase.The coefficients e, f , g and h are given in Table 3 for each horizontal-flow regime (except Transition).
When the horizontal-flow pattern is found to be Transition, the liquid holdup should be interpolated between the values calculated for the Segregated (HL Seg ) and Intermittent (HL Int ) regimes: A normalizing friction factor f n can be obtained from the smooth pipe curve on a Moody diagram or from:  The Reynolds number in Equation ( 10) is given by: The two-phase friction factor, f tp , is determined from the relation: in which: and: When 1 < y < 1:2, S is given by: The three components which make up the total pressure gradient are determined next.The pressure change resulting from friction losses at the pipe wall is given by the following equation: The pressure gradient due to hydrostatic head effects is given by: in which: The pressure gradient due to kinetic energy effects is often negligible but should be included for increased accuracy.This can be done by use of a dimensionless acceleration term E K : The total pressure gradient is given by: Payne et al. [14] proposed two well accepted modifications to improve liquid holdup and pressure-drop predictions by this correlation.They found that the Beggs and Brill method overpredicted liquid holdup in uphill and downhill flow and recommended correction factors of 0.924 and 0.685, respectively.They also recommended that the normalizing friction factor f n be obtained from a Moody diagram, as the use of the smooth-pipe equation results in underpredicted values of two-phase friction factors.All calculations performed in this study include these modifications.

CONSERVATION OF ENERGY AND HEAT TRANSFER
It is clear from Section 1.1 that the integration of the pressure gradient over the pipe length requires values of the physical properties of the liquid and vapor phases as a function of the axial coordinate.These physical properties, in turn, are functions of temperature and pressure.Hence, knowledge of temperature as a function of the axial coordinate is also necessary.When this temperature distribution along the pipeline is unknown, which is often the case in design calculations in the petroleum industry, it can be determined through a rigorous enthalpy balance.Application of the energy conservation principle to a steady-state flow through a pipe segment leads to Equation (23) for one-dimensional enthalpy gradient [1]: in which Q is the heat flux between the flowing fluids and the pipe surroundings, w is the total mass flow rate, v is the velocity and h is the specific enthalpy of the fluid.
The heat flux Q is the product of the overall heattransfer coefficient U and the temperature difference between the fluids (at temperature T Þ and the surroundings (T e ): If the pipe is exposed to ambient air, which happens quite often in refining units, U can be calculated from [15]: in which D is the outer diameter of the pipe segment, k is the thermal conductivity and h i and h o are coefficients of convective heat transfer between the flowing fluids and the pipe internal wall and between the pipe outer surface and the surroundings, respectively.

CALCULATION OF PRESSURE TRAVERSES
From the above discussion, it is clear that a typical pressure drop calculation (when the temperature profile along the pipe is not known in advance) requires the numerical integration of the two-phase pressure gradient simultaneously with Equation (23).
A number of publications on the subject of two-phase flow in pipes which have become popular in the petroleum industry [1,2,4] present an iterative technique for carrying out these integrations, which [1] refers to as a "marching algorithm".First, the pipe is divided into m increments of length Áx, small enough that the pressure and enthalpy gradients can be considered constant within each of them.
Then, starting at one end of the pipe, where both pressure and temperature are known, the pressure and temperature at the end of each of the m segments are sequentially solved for, in a trial-and-error fashion.Complete pressure and temperature profiles will have been calculated after this has been done for all m increments, and the total pressure drop can finally be determined.
Figure 1 shows the calculation steps that need to be performed for each segment of the pipe, assuming that the upstream conditions of the increment (denoted by subscript i) have already been determined and that downstream pressure and temperature (subscript i þ 1) are to be calculated.The PVT-calculation (Pressure/ Volume/Temperature) box involves phase-equilibria calculations as well as those of physical and thermodynamic properties (densities, viscosities, surface tension and enthalpies).In short, this strategy, which is very similar to a forward Euler integration, implies solving for the unknown variables using two separate, nested loops with the algebraic thermodynamic and heat-transfer constraints being solved in each iteration of the inner loop, rather than numerically solving all the equations in a true simultaneous fashion.
It is one goal of this study to emphasize and demonstrate that the "marching algorithm" is not the most efficient technique available for coupling the relevant differential equations (pressure gradient and enthalpy gradient).These equations and those required for PVT and heat-transfer calculations all form a Differential Algebraic Equations (DAE) system, and they must be represented as such in order to be simultaneously solved by more efficient, suitable numerical methods.
Pressure, temperature and enthalpy can be considered the main variables of this problem, since all relevant equations can be expressed as functions of them.Thus, a minimum state vector would be: Three equations must be solved along the pipe for these variables: where f 1 is a two-phase pressure-gradient correlation (in this study, Eq. 22, i.e., the Beggs and Brill correlation), f 2 is the right-hand side of Equation ( 23) and g 1 is a  thermodynamic expression for h: Equation (29) can be rewritten as: Equations (27, 28) and (30) can be combined by use of a mass matrix M [16]: in which: and: MATLAB's ODE15S built-in stiff solver was used in this work for solving the DAE System (Eq.31).This approach is expected to prove more efficient in that all differential equations and algebraic constraints will be simultaneously processed, instead of sequentially being iterated for in a nested structure.CPU-cost gains are typical benefits of not separately solving algebraic equations in every step, as verified by [17] in a flash-drum DAE problem.

METHODS
Unless otherwise indicated, all reported calculations were performed considering 115 000 kg/h of naphtha entering a horizontal, 30-meter long carbon steel pipe segment (k = 60.5 W/m/K, d = 202, 7 mm and D = 219,1 mm) at 168°C and 5.0 kgf/cm 2 .A value of 27°C was assumed for ambient air temperature.These latter values represent typical refinery design conditions.Table 4 gives the assumed naphtha composition, which is very similar to that of the gas chromatography example presented in [18].The Peng-Robinson equation of state was used in calculations of Vapor-Liquid Equilibrium (VLE), vapor-phase density and residual enthalpies.Enthalpy calculations by the latter method also require determination of ideal gas state heat capacities at constant pressure, which was done using equations and constants given in [19].Liquid-phase density, viscosities of liquid and vapor phases and surface tensions were all calculated using methods presented in [18] and [20] for petroleum fractions of defined composition.
In two-phase flow, radial convective heat transfer between the flowing fluids and the internal wall of the pipe depends on flow pattern, which makes it more complex than for single-phase flow.But according to [4], in most two-phase flow problems, flow is turbulent enough that convective heat transfer becomes significant and temperature differences between the flowing mixture and the wall are negligible.This last assumption was adopted in all reported calculations, which corresponds to assigning a zero value to the internal convection resistance in the bracketed term in Equation (25).
The convective heat transfer coefficient between the outer surface of the pipe and the surroundings must be calculated from a Nusselt number, which in turn should be determined from equations for natural or forced convection, whichever be the case [15].The correlations of Churchill and Chu [21] and Churchill and Bernstein [22] for natural and forced convective heat transfer were used for estimating h o under stable and unstable (36 km/h wind velocity) atmospheric conditions, respectively.Occurrence of wind was considered by default.

RESULTS AND DISCUSSION
The two-phase flow problem must be solved both by the "marching algorithm" and the DAE approach if these methods are to be compared in terms of results and CPU requirements.Seeing as the "marching algorithm" was expectedly found to be somewhat intensive with respect to computational effort -its execution took 206 seconds when the pipe was subdivided into 100 segments -it was initially sought to find a reasonable value of m (number of pipe increments) representing a compromise between calculated results and CPU cost, in order to ensure that comparisons of both methods were not biased by excessive pipe subdivision here.
Changing m from 100 to 10 drastically reduced execution time (to 25 s) at the cost of very slight changes in the pressure at the outlet of the pipe (À0.0033%) and other calculated results, as shown in Figure 2. Further reduction of m to 5 led to a smaller relative change of execution time (to 14 s) while the difference of downstream pressure to that of m ¼ 100, although still small, increased approximately five times (to 0.0153%).

Discontinuities and Lack of Convergence
Attempting to solve the same problem by the DAE method failed as convergence was not achieved at a distance close to 15 m from the pipe inlet.
The results calculated from some two-phase flow correlations show rather severe discontinuities which often lead to lack of convergence in simultaneous pressure and temperature calculations, as detailed in [2].In the Beggs and Brill method, this occurs when flow conditions change just enough that a different flow regime gets predicted, which causes different values of coefficients a, b, c, e, f , g and h to be used when calculating liquid holdup (Tab.2, 3).
Investigation of the results previously calculated by the "marching algorithm" shows that from x ¼ 15 m to x ¼ 18 m, k L changes from 0:0732 to 0:0723, causing L 1 to vary from 143:47 to 142:94.N Fr changes from 142:87 to 145:89, thus becoming greater than L 1 .Consequently, according to Table 1, the horizontal-flow regime changes from Intermittent to Distributed.Figure 3 shows the resulting abrupt change of liquid holdup.

Usage of Regularization Functions
Large discontinuities in empirical models are known to cause difficulties to integrator solvers indeed, and the use of specific interpolation functions for linking two adjacent discontinuous domains in order to come around this problem has been investigated in recent studies [23,24].The regularization function technique used in [25] when approximating discontinuous jumps in consistent initializations of DAE systems has also been shown by [24] to successfully replace model discontinuities while even reducing computational effort.In [24], discontinuities such as: are replaced by a single algebraic equation: where g arg; d ð Þ is a continuous regularization function evaluated at scalar argument arg and parameter d.In order for Equation (36) to closely reproduce 35, g must be such that: where 0 < n << 1.When Àn arg n, g assumes an intermediate value in the 0; 1 ½ interval.n is a function of user-supplied parameter d, which must be adjusted to each particular application so that the interval Àn; n ½ (inside which Eq. 36 does not exactly reproduce Eq. 35) is small enough that deviations between both equations can be neglected.
This scheme ensures a continuous and smooth transition from z 1 X ð Þ to z 2 X ð Þ while closely resembling the original model.One possible regularization function, which was adopted in this study, is given by: Figure 4 illustrates how the magnitude of n can be adjusted through tuning of d when using regularization function (38).

Regularization Functions and Discontinuities at Flow-Regime Boundaries
The horizontal-flow regime is a function of k L and N Fr , as can be seen from Table 1.Any flow-regime inequalities can be easily combined with Equation (38) into a continuous function of these variables which yields a value of 1 when the given regime is to be predicted, and 0 when it is not.Such a function for Segregated flow-regime would be: Table 5 gives similar functions for all other horizontalflow regimes.
A continuous equation for calculating liquid holdup for any flow regime is: where A is now determined from: In Equation ( 40), HL 1 and HL 2 are given by Equation ( 7), and subscript 2 always denotes Intermittent flow-regime holdup calculation.Subscript 1, in turn, denotes Segregated holdup calculation when Transition flow-regime gets predicted, and whichever flow regime is detected when it is not Transition.Tuning of regularization function d .The empirical coefficient a that is used in calculation of HL 1 should also be determined from its corresponding continuous function:

Segregated
and for calculation of HL 2 : Similar equations for all empirical coefficients of the Beggs and Brill liquid holdup correlation are given in Table 6, where checking for downhill flow is also done through the regularization function approach, so that discontinuities of e, f , g and h with respect to h are dealt with as well.This is accomplished by defining a function which yields 1 for downhill flow, and 0 for horizontal and uphill flow: with h given in degrees.The 0.01 offset in Equation ( 44) ensures that Dow 0 ð Þ is always zero, and not 0:5 (because A value of d ¼ 10 À7 was used throughout this study, as its associated value of n was shown to be much smaller than 10 À3 (Fig. 4), which is negligible enough next to angles in degrees and calculated values of L 1 , L 2 , L 3 and L 4 for k L close to 0:07.
Figure 5 shows empirical coefficients a 1 and b 1 calculated by this method for k L ¼ 0:073 (L 1 ¼ 143:35, also shown in this figure as a continuous function) and N Fr between 142 and 145, where the horizontal-flow regime changes from Intermittent to Distributed, as previously discussed.It can be clearly seen that even though the TABLE 6 Empirical coefficients expressed as continuous functions of k L , N Fr and h developed equations (such as Eq.42) are continuous functions, they closely resemble the original correlation discontinuities.

Results with Regularization Functions
Figure 6 shows a perfect match between "marching algorithm" results obtained with and without the use of regularization functions.No significant changes in previously reported execution times were observed for the "marching algorithm".Most importantly, using regularization functions, convergence was attained when solving the problem through the DAE approach, and the CPU time was 9:65 seconds, corresponding to 38:6% of the "marching algorithm" measured time for m ¼ 10 (and 68:9%, when m ¼ 5).These results (and their perfect agreement with those of the "marching algorithm") are illustrated in Figure 6 as well.Figure 7 confirms that the usage of regularization functions enabled the DAE method to obtain convergence while correctly predicting the discontinuous results of liquid holdup.
These calculations were also performed considering an uphill flow through a vertical pipe.In this problem, the flow pattern also changes from Intermittent to Distributed (this time, close to x ¼ 6 m), and the DAE solver again fails to converge when regularization functions are not used.When they are enabled, identical results to those of the "marching algorithm" are reached after 10 seconds (38:5% of the 26 seconds required by the latter method when m ¼ 10).
All these results and observations expectedly attest to the superior efficiency of the DAE approach in solving the two-phase flow problem, as this takes approximately one third of the time required by the "marching algorithm".The usage of regularization functions has also been proven an effective, easy-to-implement solution for preventing convergence failure of numerical integrators while closely reproducing a discontinuous pressuregradient correlation.Accurate temperature prediction is known to be very important in oil and gas production two-phase flow problems [4], which often involve transportation of two-phase mixtures over long distances and high temperature gradients.Approximate analytical solutions for calculation of temperatures in wells have also been published [1], which provides an alternative to performing rigorous heat transfer calculations.In short, there is a reasonable number of studies on this matter which engineers engaged in well design can call upon.
On the other hand, little work on this subject has been published when it comes to refining and petrochemical two-phase flow.Consequently, process engineers involved in downstream design are often dubious on how strictly they should consider heat transfer calculations, or even if they need to be considered at all for calculating pressure drops with acceptable engineering accuracy.In order to develop a brief first insight into this issue, the numerical techniques discussed above were used to solve the two-phase flow problem with three additional heat-transfer extreme approaches: -natural convection, external air temperature (T e ) of 40°C, -forced convection, external air temperature of 15°C, -constant mixture temperature (isothermal flow).
The first and the second scenarios are typical summer and winter conditions in a Brazilian refinery.
Obtained results are illustrated in Figure 8, where they are also compared with previously calculated values.It can be observed that nearly identical results were predicted in all three calculations where heat transfer was considered (all scenarios except isothermal).Hence, carefully considering exact atmospheric conditions is not important in this problem when calculating pressure drops for design purposes, as long as the energy balance is considered.
Disconsidering heat transfer and assuming isothermal flow greatly affected vapor fractions along the pipe.This is expected because in this case, the pressure drops while the temperature remains constant rather than dropping as well, which leads to more intense vaporization.The corresponding pressure curve deviates from the other three in a similar trend to that of the vapor fraction, but in the end, the total pressure drops do not differ too significantly (4:44%).
Similar results were also predicted for vertical upward flow.It follows that, as long as the energy balance is considered, specifying exact convection conditions does not seem necessary for calculating pressure drops in uninsulated, exposed, downstream-industry pipelines.Consideration of isothermal flow, on the other hand, led to dramatic changes in predicted results such as vapor fractions, even though predicted pressure drops were acceptably close to heat-transfer associated results.Therefore, further investigation should be conducted before isothermal-flow assumption can be said to adequately dismiss rigorous enthalpy balances.

SUMMARY AND CONCLUSIONS
The reasons for choosing practical pressure-gradient correlations over more advanced techniques (such as CFD) in specific situations in industrial design for solving two-phase flow problems nowadays have been discussed.These two-phase flow methods require values of physical properties of both flowing phases (liquid and vapor), which also depend on temperature in addition to pressure.Therefore, the pressure gradients obtained from these correlations must be integrated over the pipe length simultaneously with an enthalpy-gradient equation if pressure drops are to be calculated.It has been proposed that the "marching algorithm", which some publications present for coupling the pressure-gradient and enthalpy-gradient equations, is not the most efficient numerical method for solving this two-phase flow problem, as it requires that pressures and temperatures be separately iterated for in nested loops, with the algebraic thermodynamic and heat-transfer constraints being solved at each step of the inner loop.Representing the relevant equations as a Differential Algebraic Equations (DAE) System, thus enabling them to be simultaneously solved by more advanced, suitable techniques, has been demonstrated.The proposed DAE approach expectedly performed faster than the "marching algorithm" in all reported comparison calculations.
Empirical two-phase flow pressure-gradient correlations such as Beggs and Brill's (which was adopted in this study) have long been known to bear mathematical discontinuities which can lead to convergence failure of numerical integrators, and this did indeed occur when using the DAE solver as first attempted.
The use of regularization functions in closely reproducing discontinuous mathematical models was presented, and this was proposed as an effective technique for coming around this lack of convergence problem.The discontinous equations of the Beggs and Brill correlation were rewritten in terms of the regularization functions, and this approach was successful in preventing convergence failure while not noticeably altering the correlation's original results.
Seeing as little work has been published when it comes to two-phase flow in downstream petroleum industry, the numerical techniques discussed above were used to analyse the influence of different heat-transfer approaches on calculated results.The goal here was to shed some light on the actual need for rigorous heat-transfer calculations in downstream industry design, as engineers are frequently faced with this issue.
It was found that for the given problem, as long as heat transfer was considered, it was not necessary to specify precise ambient air temperature or atmospheric conditions (natural or forced convection), as oppositely extreme considerations for these parameters yielded virtually indistinguishable results.Replacing the energy balance with an isothermal-flow consideration resulted in small differences in predicted pressure drops, but dramatically changed other results, such as vapor fractions along the pipe.Hence, this assumption requires further testing before it can be reported as good design practice.
Finally, it must be stressed that even though a typical refinery problem was considered in all calculations performed in this study, the numerical techniques presented are just as relevant for oil and gas production design, where improvements in efficiency and convergence properties are also expected.

Figure 2 "Figure 3
Figure 2 "Marching algorithm" results along the pipe (m ¼ 10 and m ¼ 100) for pressure and temperature.

Figure 5
Figure 5 Discontinuities in empirical coefficients a 1 and b 1 captured by the regularization function approach.

Figure 6 Results
Figure 6Results along the pipe for a) pressure and temperature, and b) specific enthalpy and vapor mass fraction, calculated by the original correlation ("marching algorithm" with m ¼ 10) and with regularization functions ("marching algorithm" and DAE approaches).

Figure 7 Identical
Figure 7Identical results calculated both by the DAE solver and by the "marching algorithm" for discontinuous liquid holdup.

Figure 8
Figure 8 Effect of several heat-transfer-calculation approaches on results along the pipe for a) pressure, b) temperature, and c) vapor mass fraction.

TABLE 1 Flow
regime limits in Beggs and Brill horizontal-flow-pattern map Segregated k L < 0.01 and N Fr <L 1 or k L ! 0.01 and N Fr < L 2 Transition k L ! 0.01 and L 2 N Fr L 3 Intermittent 0.01 k L < 0.4 and L 3 < N Fr L 1 or k L ! 0.4 and L 3 < N Fr L 4 Distributed k L < 0.4 and N Fr !L 1 or k L ! 0.4 and N Fr > L 4

TABLE 4
R.G.D. Teixeira et al. / Two-Phase Flow in Pipes: Numerical Improvements and Qualitative Analysis for a Refining Process Teixeira et al. / Two-Phase Flow in Pipes: Numerical Improvements and Qualitative Analysis for a Refining Process