Hydrodynamic Modeling of Oil–Water Stratiﬁed Smooth Two-Phase Turbulent Flow in Horizontal Circular Pipes

: In the petrochemical industry, multiphase ﬂow, including oil–water two-phase stratiﬁed laminar ﬂow, is more common and can be easily obtained through mathematical analysis. However, there is no mathematical, analytical model for the simulation of oil–water ﬂow under turbulent ﬂow. This paper introduces a two-dimensional (2D) numerical simulation method to investigate the pressure gradient, ﬂow ﬁeld, and oil–water interface height of a pipeline cross-section of horizontal tube in an oil–water stratiﬁed smooth ﬂow, which has ﬁeld information of a pipeline cross-section compared with a one-dimensional (1D) simulation and avoids the signiﬁcant calculation required to conduct a three-dimensional (3D) simulation. Three Reynolds average N–S equation models ( k − ε , k − ω , SST k − ω ) are used to simulate oil–water stratiﬁed smooth ﬂow according to the ﬁnite volume method. The pressure gradient and oil–water interface height can be computed according to the given volume ﬂow rate using the iteration method. The predicted data of oil–water interface height and velocity proﬁle by the model ﬁt well with some available experiment data, except that there is a large error in pressure gradient. The SST k − ω turbulence model has higher accuracy and is more suitable for simulating oil–water two-phase stratiﬁed ﬂow in a horizontal pipe.


Introduction
In petroleum transportation, the oil-water two-phase flow is very economic and common technique [1]. The pressure gradient, oil-water interface height, and velocity field are important factors to design the pipe [2] and separator [3]. Pressure gradient is an important basis for designing the wall thickness of pipelines and pressure vessels. Oil-water interface height is used to design the carrying capacity of downstream facilities. Accurate prediction of the velocity field of the pipeline cross-section is helpful for controlling flow assurance issues, such as wax deposition [4], hydrate formation [5], pipeline corrosion [6], etc. Therefore, in order to optimize piping design and equipment operating parameters, two-phase flow has been widely studied. There are multiple flow patterns in oil-water pipe flow, such as stratified flow, dispersed flow, and annular flow, for different fluid properties and flow conditions. Among the above flow patterns, stratified flow is the most common and flow model of gas-liquid to oil-water liquid-liquid flow to research the non-Newtonian characteristics of oil-water stratified flow on horizontal pipes. However, compared with the gas-liquid two-phase flow, the oil-water two-phase flow has significant differences, mainly due to the small oil-water density difference, the small viscosity difference, and the small velocity difference between the two phases. In fact, the shear stress between the oil-water two-phase interface is stronger than the gas-liquid two-phase for stratified flow. If a method similar to that of the gas-liquid phase being divided into two calculation domains is used to discretely solve and update the interface stress, not only is the solution speed slow, but the velocity field distribution of pipe cross section will also have a large error. Consequently, these established models lack sufficient ability to accurately describe oil-water two-phase stratified pipe flow.
To fill the research gap above, the primary intention of this work is to investigate oil-water two-phase pipe flow through establishing a steady-state 2D oil-water two-phase stratified flow model in a bipolar coordinate system, where three turbulence models (k − ε, k − ω, SST k − ω) are utilized to simulate the turbulence flow. Then, the finite volume method is used to discrete the model and an adaptive mesh method is utilized to track the oil and water interface. After the convergence of velocity and turbulence field, the total mass flow conservation is used to adjust the pressure gradient and oil-water interface height. In addition, pressure gradient, oil-water interface height, and velocity field are compared with experimental data to validate the presented hydraulic model under three different turbulence models.

Problem Description and Model Assumption
In this section, a 2D isothermal horizontal pipe oil-water stratified flow model is demonstrated. The governing formulas are composed of mass conservation equations, momentum conservation equations, and turbulence models. Oil transportation is usually carried out by mixing oil and water as shown in Figure 1. In the later stage of oil exploitation, oil transportation is mostly in a stable two-phase stratified flow of oil and water in the pipeline. Crude oil contains colloidal asphaltenes and other impurities, which will affect the flow assurance issue, and pressure gradient, oil-water interface height, and velocity field of the pipe cross section are important parameters for oil-water flow assurance. Therefore, it is important to establish a steady hydraulic model of oil-water two-phase pipe flow. Assumptions are made as follows: (1) The oil and water are in a smooth and steady stratified flow, and there is no emulsified layer at the interface; (2) the oil phase is homogeneous and does not have non-Newtonian characteristics; (3) the influences of heat transfer are negligible; (4) the flow of oil and water is fully developed at steady-state.

Mass Conservation Equation
For simplified calculation, the model's premise is that the emulsification between oil and water, non-Newtonian characteristics, and thermal transmission are not considered. As the oil and water are incompressible, the total flow rate of the oil and water phases remains unchanged on each section of the pipe axis. Therefore, the mass conservation equation can be obtained by where w oil and w water are the axial velocities of the oil and water phases, m/s. A oil and A water are the circulation areas of the oil and water phases, m. Q oil and Q water are volume flow rates at the inlet of the oil and water phases, m 3 /s.

Momentum Conservation Equation
The Navier-Stokes formula of fully developed oil-water stratified pipe flow along the axial direction of the pipeline is depicted in Equation (2), which assumes that the oil-water two-phase flow is incompressible.
where µ e is the effective viscosity coefficient of the oil or water phase, µ e = µ m + γµ t , and N·/m 2 . dp dz is the flow pressure gradient, Pa/m. A single-pressure model is used to calculate the momentum conservation equation.

Turbulence Model
Oil-water pipeline flow is often turbulent, and the characteristics of turbulence are more complicated than laminar flow. At present, the commonly used turbulence models are large eddy simulation and Reynolds average. Compared with the Reynolds time average, the large eddy simulation requires extensive calculation and significant time. Hence, for quickly obtaining the oil-water flow field information on the cross section of the pipeline, three different Reynolds average models (k − ε, k − ω, SST k − ω) are adopted in this study. The governing equations for the k − ε model [33] are as follows: where the closure coefficients are given by C u = 0.09, C ε1 = 1.44, C ε2 = 1.92, σ k = 1.0, and σ ε = 1.3. The governing equations for the k − ω model [34] are as follows: where µ t is the dynamic eddy viscosity, N·/m 2 ; ρ is the density, m 3 /s; and the closure coefficients are given as α = 5/9, β = 0.075, β * = 0.09, and σ k = σ ω = 2. The governing equations for the SST k − ω model are as follows: where S = 2S ij S ij is the modulus of the mean rate-of-strain tensor and P k is a limited production term given by the following equation: F 1 and F 2 are blending functions defined by and The closure coefficient of the SST k − ω model [35] is obtained by blending those of the k − ω model, denoted by φ 1 , and the standard k − ε model, denoted by φ 2 , via the . The coefficients are given by a 1 = 5/9, a 2 = 0.44, β * = 9/100, β 1 = 3/40, β 2 = 0.0828, σ k1 = 0.85, σ k2 = 1, σ ω1 = 0.5, and σ ω2 = 0.856.

Boundary Conditions
On the pipe wall, the nonslip boundary condition is for velocity, turbulent kinetic energy, and eddy viscosity.
When the Reynolds number is low, the turbulence model can capture the damping effect near the wall. If the Reynolds number exceeds a specific value, the thickness of the viscous sublayer of the boundary layer will become thinner, resulting in a small number of grid points, making it difficult to obtain the results. The wall function [36,37] is employed to avoid the problem. The energy dissipation rate of the point near the wall in the k − ε equation is described as follows: The specific energy dissipation rate in the k − ω equation at the boundary is given by The SST k − ω equation combines two specific energy dissipation rates to obtain the near-wall value.

Governing Equations Discretization
The governing equation is discretized by finite volume in the bipolar coordinate system, and a nonuniformly distributed grid is used in the calculation domain [a, b] and [c, d]. The diffusion term is discretized using the central difference scheme.
The equation form after sorting the discretization is where The coefficient a p is given by The Cartesian coordinate system and the bipolar coordinate system have the following coordinate transformation relationship: where ξ and η are the coordinates of point P in bipolar coordinates. The scalar factor of the bipolar coordinate is given by A bipolar coordinate system is introduced to simulate the two-phase noncircular and irregular regions in the stratified flow more accurately ( Figure 2). Through bipolar coordinate transformation, the oil-water stratified pipe flow zone in Cartesian coordinates is transformed into a regular rectangular calculation zone in bipolar coordinates. The calculation zone of the oil phase is θ 1 < ξ < θ 2 , −∞ < η < ∞. The calculation zone of the water phase is θ 2 < ξ < π + θ 1 , −∞ < η < ∞. In numerical computation, the value of η must be an exact value; Newton and Behnia [24] suggested η max = 6.

Solution Methodology
It is necessary to calculate pressure gradient and oil-water interface height to make the velocity field of pipe cross section match the given volume flow. The pressure gradient affects the solution of the momentum equation, and the height of the oil-water interface controls the generation of the adaptive mesh. As shown in Figure 3, the hydrodynamics model is nonlinear. Hence, an iterative solution is required to obtain the pressure gradient, the oil-water interface height, and the velocity field of the pipe cross section that satisfies the mass conservation equation.
The solution of the model is depicted in the flow chart of Figure 3. First, we set the boundary conditions; then, we estimate the oil-water interface height and pressure gradient, and divide the calculation grid. After the grid is divided, the momentum equation and turbulence model can be solved. Thus, the velocity, kinetic energy, and specific energy dissipation region of each phase can be obtained. If the results meet the mass conservation conditions, the calculation ends; otherwise, the oil-water interface height and pressure drop gradient are adjusted, and the calculation continues. Next, we obtain the pressure gradient dp/dz and dimensionless oil-water interface height h l , which is the ratio of the height of the oil-water interface to the pipe diameter. The first is obtained by external iteration calculation, which adopts the Newton-Raphson method to change dp/dz and h l to satisfy the mass conservation equation. Therefore, we can obtain the oil and water two-phase continuity equations as follows: Subsequently, m = dp/dz and n = h l are introduced to simplify notation. The Newton-Raphson iteration method is applied to calculate the solution of the nonlinear equations represented by Equation (25). The next approximate value is calculated by Based on the analysis of Equation (26), the calculation process is iterative. Consequently, it is critical to set the error range reasonably to solve the equation, including the convergence speed, accuracy, calculation results, and so on. When the continuity equation satisfies a specific error range, a reasonable calculation result is obtained.

Grid-Independence Solution
The governing equations are solved on the structure grid, and there is local refinement at the oil-water interface and pipe wall. In order to eliminate the influence of the grid on numerical results, it is necessary to test the grid independence. Through experimentation, Yusuf [38] observed that the oil-water two-phase flow was stratified flow in a 25.4 mm acrylic pipe, when the water superficial velocity is 0.318 m/s and the oil superficial velocity is 0.14 m/s. As shown in Figure 4 and Table 1, the velocity distribution at vertical centerline on the cross section, pressure gradient, and dimensionless oil-water interface height are calculated and compared by using the proposed algorithm in Yusuf experiments under four sets of grid.
It can be seen from Figure 4 that the velocity distribution under different grid sets has a consistent distribution. However, as shown in Table 1, there are obvious differences in pressure gradient and dimensionless oil-water interface height from 20 × 20 to 50 × 50 mesh. Pressure gradient and dimensionless oil-water interface height only changed by about 1%, when the grid changes from 50 × 50 to 80 × 80, and from 80 × 80 to 100 × 100. Therefore, in this paper, the 80 × 80 mesh is selected to simulate all the cases.

Results and Discussion
The model must only input the pipe diameter d, physical properties, and volume flow rates of oil and water phases. Then, the axial pressure gradient, dimensionless oilwater interface height, and velocity and turbulence field of the pipe cross section can be calculated. This section is divided into three subsections. The first subsection validates the turbulence model under a single-phase flow scenario. The second subsection presents the calculation results of the pressure gradient of the two-phase flow and compares them with the experimental values. The last subsection validates the flow field and interface height of the pipe cross section.

Turbulence Model Validation
Here, we simulate oil-water two-phase turbulence, and verify the reliability of the algorithm and turbulence model. First of all, a different Reynolds single-phase flow hydraulic calculation is utilized to validate the turbulence models. The calculation result of the Darcy formula was adopted to validate turbulence models. Then, the test includes distributing the identical fluid properties and volume flow rates for both phases. Different flow models were used at different flow regimes, such as laminar and turbulence.
As shown in Figure 5, the pressure gradient of single-phase pipe flow is given under different models. When the turbulence model is coupled, the pressure gradient obtained by the numerical calculation algorithm are completely consistent with the theoretical formula under laminar flow conditions. Then, for the prediction of turbulence in singlephase flow through coupled turbulence models, the calculation outcomes match well with expected values from the Darcy formula, but the pressure gradient was underestimated at a high Reynolds number. Thus, it can be proved that the turbulence models have sufficient accuracy and the algorithm can obtain the correct pressure gradient based on the single-phase flow test.

Calculation of Two-Phase Flow Pressure Gradient
Before the flow field of the pipe cross section is investigated, the two-phase pressure gradients under different oil-water flow rates with different turbulence models need to be verified. Angeli [39] researched the oil-water two-phase stratified flow in a pipe with an inner diameter of 24.3 mm and length of 9.7 m. The experimental results illustrate that when the mixing velocity is less than or equal to 0.6 m/s, the two-phase flow pattern in the acrylic pipe is stratified flow. Therefore, the superficial velocities of the two phases are chosen to be between 0.11 and 0.55 m/s to satisfy the stratified flow pattern of oil-water two-phase flows. Fifteen experimental data points are applied to confirm the accuracy of the models. Figure 6 and Table 2 illustrate a comparison between prediction and experimental values of a pipe pressure gradient at oil-water two-phase superficial velocities under different turbulence models. Average percentage deviation (APD) and standard deviation (STD) are employed to evaluate the error [40]. (28) where N is the number of flow conditions. The pressure drop predictions at different turbulence models present differences. All three models accurately forecast the pressure gradient given specific oil-water physical properties and flow rates, although the prediction errors of some points exceed 15%. The prediction result is not sufficiently accurate at low superficial velocities. For example, when the superficial velocity of oil and water both are 0.11 m/s, the relative error is approximately 30%. At lower velocities, the turbulence models predict the fluid transition with low accuracy. The flow pattern is given without a picture. Furthermore, it is difficult to observe the specific state of the separated flow pattern. When the superficial velocities of the oil and water phases are relatively close, the prediction model has a high prediction accuracy. From theoretical analysis, when the superficial velocity distinction between the oil and water phases is greater, the oil and water interface is likely to produce oil or water droplets. The droplets produced by entrapment will increase the local viscosity, increasing the pressure gradient. The low prediction accuracy of No. 7, 12, and 15 may be caused by the generation of droplets. It can be found that the SST k − ω average percentage deviation (APD) and standard deviation (STD) of SST k − ω turbulence model are 8.9% and 7.6% smaller than other models. The SST k − ω turbulence model has the highest accuracy of the three models.

Calculation of Interface Height and Flow Field
In the previous two subsections, the pressure gradients of the single-phase flow and two-phase flow were verified. The results illustrate that the accuracy of the prediction model is within 15%, which can satisfy industrial applications. In this subsection, the oil-water interface height and flow field are studied. Kumara [41,42] utilized particle image velocimetry (PIV) and laser Doppler anemometry (LDA) to measure an oil-water two-phase flow field in a horizontal pipe with a diameter of 56 mm and length of 15 m. Three stratified flow conditions were selected for verification. Figures 7-9 illustrate a comparison between prediction and experimental values of a pipe velocity field at different two-phase mixed velocities and water cut under different turbulence models. Figure 7. Velocity distribution at pipe centerline under different turbulence models given mixed velocity 0.5 m/s and water cut 0.5 [42].
As illustrated in Figure 9, the velocity distribution near the centerline of the pipeline agrees with the measured value of the experiment. However, when the mixed velocity is 0.68 m/s, there are more droplets at the interface, which affects the local viscosity distribution. Thus, a deviation exists between the numerical results and the experimental measurement values.  Although there are still droplets on the interface under a lower mixing velocity, the vortex generated by the flow is not strong. This can be verified due to the agreement between the numerical prediction outcomes and the experimental measurement data. Hence, the influence of local droplets on viscosity can be negligible. Therefore, this model can be used to accurately forecast the pressure gradient, oil-water interface height, and flow field of the oil-water stratified flow.
Furthermore, Figures 10-12 illustrate the velocity field and turbulence field under the SST k − ω turbulence model when mixed velocity is 0.43 m/s and water cut is 0.25. There is a larger turbulent energy and energy dissipation rate near the pipe wall. These two parameters directly affect the turbulent eddy viscosity. There is a strong momentum exchange among the boundary layer and the turbulent core area, and the turbulent intensity is high. The turbulent energy dissipation rate (ε = C u kω) determines the dissipation rate per unit of turbulent kinetic energy. The vortex generated by the turbulence is transferred from the large vortex to the small vortex and gradually collapses. This is also consistent with the theory. Furthermore, the oil and water are assumed to be in a smooth stratified flow, and there is no momentum exchange between the two phases. The vortex structure in the turbulent core area is relatively weak.

Conclusions
The pressure gradient, oil-water interface, and field variables of pipe cross section in designing an transportation system for petroleum exploitation is of great significance. In this paper, a new, numerical algorithm is considered for stratified oil-water stratified flow in a horizontal pipe and compared with the calculation outcomes from experimental data. The results indicate that the model can well predict the oil-water interface height, velocity field, and turbulence variable fields at the pipe cross section given specific oil and water flow rates and physical parameters. The predicted outcomes agree closely with the experimental data. Although the oil-water interface height and flow field of pipe cross section have good accuracy, the prediction of pressure gradient has relatively larger error. The pressure gradient average percentage deviation of k − ε, k − ω, and SST k − ω are 10.8%, 9.6%, and 8.9%, respectively; the standard deviations are 9.8%, 8.1%, and 7.6%. At the same time, comparing the oil-water interface height and velocity distribution at vertical centerline, SST k − ω is better than the other two turbulence models.
The results imply that the SST k − ω model is more accurate than the k − ω and k − ε turbulence models. Further, the program does not need to run on clusters [25] and a result can be obtained within several minutes. The fast calculation speed can be effectively applied to the petroleum industry. Our code is available online, which will help more scholars expand on and verify the model (https://github.com/winsway/StratifiedFlow-KQ, accessed on 16 July 2021). In future work, the phase field model is demanded to be involved about the influence of the droplet emulsification.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
The following nomenclature are used in this manuscript: