Numerical Modelling of Dispersed Water in Oil Flows Using Eulerian-Eulerian Approach and Population Balance Model

: Formation of a dispersed oil—water ﬂow pattern is a common occurrence in ﬂow lines and pipelines. The capability of predicting the size of droplets, as well as the distribution of dispersed phase volume fraction is of utmost importance for proper design of such systems. The present study aims at modelling dispersed water in oil ﬂows in a horizontal pipe by employing a multi-ﬂuid Eulerian approach along with the population balance model. To this end, momentum and continuity equations are solved for oil and water phases, and the coupling between the phases is achieved by considering the drag, lift, turbulent dispersion, and virtual mass forces. Turbulent effects are modelled by employing the standard k- ε model. Furthermore, a population balance model, based on the method of class, along with the breakup and coalescence kernels is adopted for modelling the droplet size distribution. The obtained numerical results are compared to the experimental data in literature for either the in situ Sauter mean diameter or water volume fraction. A comparison among the obtained numerical results and the published experimental data shows a reasonable agreement.


Introduction
The concurrent flows of two immiscible liquids are observed in many industrial sectors such as petroleum, food, and biotechnology. In the petroleum industry, crude oil is often produced and transmitted along with water.
The deformable interface of two liquids in such flows may exhibit different configurations known as flow patterns or regimes. The type of flow pattern in a pipe strongly depends on several factors, such as the fluid and flow characteristics and the geometry of the pipes. Two-phase oil-water flow patterns in a horizontal pipe are classified into two groups based on the interfacial length scale. The flows of two liquids separated with a distinct large-scale interface are generally referred to as segregated flow patterns in which both phases retain their continuities. However, there exist situations in which one of the phases loses its continuity and the large visible interface is replaced by several small interfaces. Such a flow pattern is referred to as dispersed flow. Segregated flows are often observed at low flow rates in which the superficial velocities of the oil and water are comparatively similar and the gravity force dominates over turbulent dispersion force. In such conditions, the stratified layer of oil phase would flow on top of the stratified water layer. Segregated flow patterns may also happen for viscous oils where a continuous oil phase flows at the core of the pipe section surrounded by a continuous annular water layer in the vicinity of the wall known as core-annular flows. At high flow rates, however, the turbulent effects are dominant, which, in turn, tend to disperse one of the phases in the continuity of the other one. For a certain flow rate, fluid properties and the pipe configuration, the water cut of the system specifies the type of dispersion. For water cuts below the phase inversion revealed the significant influence of the lift force on the distribution of the dispersed droplets. Furthermore, Voulgaropoulos et al. [45] and Santos et al. [46] used the simplified 2-D homogeneous mixture model to simulate the dispersed oil-water flow patterns. Apart from these, Darihaki et al. [60] employed a one-dimensional finite difference approach to predict the water concentration profiles at different pipe cross sections. A comparison of their numerical results with those from the experiments and 3D-CFD simulations revealed a reasonable accuracy of their approach.
Numerical modellings of dispersed flows by using Eulerian-Eulerian, Eulerian-Lagrangian and mixture approaches require a prior knowledge of the size of the dispersed phase, which is usually an unknown. In all aforementioned works, the size of the droplets was assumed to be uniform throughout the pipe. Furthermore, in most of the previous works the considered sizes of the droplets were based on either the experimental measurement or try and error, which adversely limits the application of the models for other working conditions. Pouraria et al. [57,58] showed that the correlation proposed by Brauner [61] and Angeli and Hewitt [62] can be used to estimate the average Sauter mean diameter of the dispersed phase. The combination of the predicted Sauter mean diameter along with adopted closure laws for interfacial forces in an Eulrian-Eulerian framework showed a reasonable accuracy for predicting the flow patterns. Burlutskii [43] also showed that the predicted size of the dispersed droplets by using the Hinze [63] model along with Eulerian-Lagrangian CFD approach results in reasonable predictions.
It is worth mentioning that the experimental studies also indicate a reasonable accuracy of the correlation proposed by Brauner [61] and Hinze [63] models in estimating the maximum droplet size when compared to the other correlations [49,64]. Nevertheless, as stated by Brauner [61], the application of these models for significantly different conditions, such as fluid properties and pipe geometry (diameter and inclination angle) may result in unrealistic predictions. Apart from that, the assumption of uniform droplet size throughout the pipe is not realistic and may adversely reduce the accuracy of the numerical results. Recently, a population balance model was employed by Bourdilon et al. [65] to study the droplet size evolution in dispersed oil-water flows. A comparison of the obtained results of the Sauter mean diameter with the corresponding experimental data showed a reasonable accuracy.
The objective of the present study is to employ the Eulerian-Eulerian approach along with population balance method for modelling the dispersed oil-water flow patterns at different operating conditions. The breakup and coalescence of droplets were modeled to determine the Sauter mean diameter of the dispersed phase. The accuracy of the numerical model was judged based on the comparison among the predicted results of either Sauter mean diameter or the in-situ distribution of the dispersed phase volume fraction with the corresponding experimental data.

Numerical Model
Numerical modelling of the dispersed water in oil flow was carried out by using Eulerian-Eulerian approach. In this approach, the conservation equations are solved for both the continuous and the dispersed phase that share the same pressure field. Furthermore, appropriate closure laws corresponding for different interfacial forces are adopted to take into account the coupling between two phases.

Governing Equations for the Model
The continuity and momentum equations in this study are as follows: where ρ k , α k and µ k are density, volume fraction and effective viscosity, respectively. Although the first term at the right-hand side of the Equation (2) shows the pressure gradient, the second and third terms show the gravity and stress. The last term indicates the interfacial forces applied on phase k. The same equations are solved for the secondary phase p.

Interfacial Forces
In the present study, four different interfacial forces, i.e., the drag, lift, turbulent dispersion, and virtual mass forces were considered to take into account the interaction between two phases. Hence, the last term in the momentum equation reads as follows: The drag force exerted on the dispersed droplets is quantified as follows: where, C d and d p indicate the drag coefficient and the diameter of the dispersed droplets.
In the present study the model proposed by Tomiyama [66] is used to quantify the drag coefficient, C d .
where the relative Reynolds number, Re, and Eotvos number, Eo, are defines as follows: The shear flow surrounding the particles exerts a lateral lift force, which is perpendicular to the direction of the droplet's motion. The lateral lift force is calculated as follows: Several models can be used to quantify the lift coefficient C L . For spherical particles, a positive value is generally adopted for the lift coefficient. Whereas, a negative value has also been adopted to reproduce the migration of large bubbles toward the pipe center. Tomiyama et al. [67] proposed a lift model that considers the influence of droplet or bubble deformation on the magnitude and the direction of the lift force. However, considering the observed discrepancies using such an approach, the present study considers a constant lift coefficient of 0.1 which has been suggested by several authors [68][69][70].
In the turbulent dispersed flows, the dispersed phase generally have a strong interaction with the turbulent eddies. These interactions promote the dispersion of the droplets and bubbles. In order to account such interactions in the Eulerian-Eulerian framework, a turbulent dispersion force is introduced to the momentum equation which is proportional to the gradient of volume fraction of the dispersed phase [71,72] Processes 2021, 9, 1345 5 of 22 In the present study, a model proposed by Lopez de Bertodano [62] is used to quantify the turbulence dispersion force, which reads as follows: → F td,k = C TD ρ k k k gradα p (9) where, C TD is turbulent dispersion coefficient, ρ k and k k denote the density and turbulent kinetic energy of the continuous phase (k) and α p indicates the volume fraction of the dispersed phase. The default value was 1 for the turbulent dispersion coefficient, C TD . The virtual mass force is computed as follows.
The term d k dt indicates the phase material time derivative and C vm is a virtual mass force coefficient which was specified as 0.5 in the present study.

Population Balance Model
A prior estimate of the droplet size is needed to proceed the numerical modelling based on the presented method in the previous section. In contrast to the previous studies where the size of droplets was assumed to be uniform throughout the computational domain, the present study aims at considering the droplet size distribution. Hence, a population balance model based on the method of class is employed to directly model the droplet size distribution in each computational cell.
The general form of population balance equation is written as follows [73]: where, B B and D B indicate the birth and death due to the breakup events and B C and D C denote the birth and death due to the coalescence of droplets. Whereas, n di indicates the number of droplets of the size n i . The source terms are given as follows: indicates the breakup rate of particle volume i into particle volume k. Whereas, Ω C v i , v j is the rate of coalescence between dispersed droplets of size i and size j. Details of the numerical model is described in Hagesaether et al. [73].
Several coalescence and breakup models have been proposed by different researchers to quantify the breakup and coalescence rates. A comprehensive review of the proposed breakup and coalescence models can be found in [74,75]. Furthermore, the accuracy of these models for different flow conditions have been assessed by several researchers. In the present study, a breakup model proposed by Lehr et al. [76] is employed since the model has proven its good performance over other models [77,78]. According to this model, the breakup occurs when the inertial energy of the hitting eddies is greater than that of the interfacial energy of the smallest daughter particle.
The breakup model proposed by Lehr et al. [66] is written as follows: In contrast to the breakup models, the coalescence models are not generally very different since most of them have been derived from the same approach. In the present study, a commonly used coalescence model proposed by Luo and Svendsen [79] is used for modelling the coalescence of the dispersed droplets.
The coalescence kernel proposed by Luo and Svendsen [79] is defined as follows: where ω C v i , v j is the frequency of collision of particles and P C v i , v j indicates the probability that the collisions result in coalescence. For two particles with diameter d i and d j and a characteristic collision velocity of u ij , the frequency and the probability of coalescence are defined as follows: where c 1 is a constant which is usually set to 1 and x ij = d i /d j , ρ 1 and ρ 2 are the density of the primary and the secondary phases. Such discretisation requires a prior estimate on the droplet size distribution (DSD) range. Hence, the maximum diameter and the Sauter mean diameter of the dispersed phase were first estimated using the models proposed by Brauner [61] and Angeli and Hewitt [62], respectively. The maximum diameter of the dispersed phase is determined using the correlation proposed by Brauner [61] which reads as follows: (20) where D and α p indicate the internal diameter of pipe and the volume fraction of dispersed phase, respectively. Furthermore, p, k, and m denote the dispersed phase, the continuous phase and the mixture, respectively. The Sauter mean diameter of the dispersed droplets is estimated as 0.48 of the maximum droplet size according to Angeli and Hewitt [52]. In the present study, the minimum and maximum droplet sizes were assigned sufficiently smaller and larger than the estimated Sauter mean diameter by using Brauner [61] and Angeli and Hewitt [62]. The droplet size was discretised by considering 10-11 classes (bins). As will be shown later, such an approach can be used to span the DSD range. Furthermore, the primary and secondary phases are assigned by considering the water cut so that for input water cuts beyond the phase inversion point the water phase is considered as a primary phase and for lower water cuts, the oil phase is considered as a primary phase. In the present study, the oil and water were assigned as the primary and secondary phases, respectively.

Turbulence Model
Turbulence significantly affects the dispersed droplet size distribution and the interfacial forces. In the present study, the turbulent effects were taken into account by employing the per-phase formulation of the standard k-ε model. In this approach, the transport equations for the turbulent kinetic energy, k and the dissipation rate ε are solved for each phase [80]. The turbulent viscosity is then calculated as follows: where C µ is 0.09. The effect of the dispersed phase on the turbulence of the continuous phase, µ t,k,i , was taken into account by employing the model proposed by Sato et al. [81]. Using this model, the turbulence induced by the dispersed droplets is estimated as follows: The constant was set to C µ,p = 0.6. The effective viscosity of the continuous phase is thus determined as the summation of the turbulent viscosity, the induced viscosity by the dispersed droplets and the molecular viscosity.
The phase coupled SIMPLE (PC-SIMPLE) algorithm implemented in ANSYS-fluent [80] was adopted for the pressure-velocity coupling. A second-order upwind scheme was used to discretize the momentum and turbulence transport equations, while the second order implicit method was employed for the transient formulation. Volume fraction and population balance equations were discretized by first-order upwind scheme. The under relaxation factors were reduced to ensure the stability and convergence of the numerical simulations. The employed under relaxation factors for the pressure, density, body force, momentum, volume fraction were 0.1, 0.3, 0.3, 0.2, and 0.15, respectively. In addition, the under relaxation factors for the turbulent kinetic energy, turbulent dissipation, and the population balance model were also set to 0.25, 0.25, and 0.1, respectively. The convergence criterial for the residuals of different equations were set to 1 × 10 −6 . Moreover, W cycle multigrid was chosen for the momentum and pressure while BCGSTAB stabilization model was activated. Comparatively small time-steps ranging from 5 × 10 −5 s to 5 × 10 −4 s were used for the transient simulations. Typical iteration numbers to reach the convergence for each time step ranged from 70 to 100. The transport equations were solved by using ANSYS FLUENT 16.1 package [80].

Geometry and Boundary Conditions
In the present study, the accuracy of the CFD model was assessed by using two sets of experimental data. The detailed information regarding the geometry of the pipes, as well as the fluid and flow properties are shown in Table 1 [62,82]. Numerical simulations were performed for two horizontal pipes with the internal diameter of 0.0563 m and 0.0243 m, respectively. In the experiment carried out by Elseth [82], oil and water were introduced into the pipe from top and bottom, respectively. Although, in Angeli and Hewitt [62] the oil and water were mixed using a T-junction at the inlet. In order to satisfy the experimental conditions of Elseth [82], oil and water were introduced separately by assigning two separate velocity inlet boundary conditions. Similar to the experiment, the water cuts and the mixture velocities in the pipe were adjusted by regulating the velocity of each phase at the inlets. One velocity inlet, consisting a mixture of oil and water, was considered as an inlet for the experiment of Angeli and Hewitt [62]. The CFD model was initialized by assuming a single phase oil flow while the dispersed phase droplets entered to the physical domain from the corresponding inlets according to the experiment. Furthermore, the hydraulic diameter and turbulent intensity at the inlets were specified based on the pipe diameter and the fluid or flow properties. The turbulent intensity at the inlet was estimated as follows: A pressure outlet boundary condition was considered for the outlet boundary and the gauge pressure was considered zero, which is compatible with the separator pressure in the experiments.
A no-slip boundary condition was applied at the wall of the pipe. A symmetric boundary condition was adopted in the present study to reduce the computational. A grid sensitivity study was performed and a grid consist of 364,000 hexahedral cells was found to be suitable for model the experiments conducted by Elseth [82]. The same grid was scaled for modelling the experiment of Angeli and Hewitt [52]. GAMBIT 2.3.16 [83] was used to generate the geometry and the grid. Figure 1 shows the employed grid, as well as the boundary conditions in the present study. velocity inlet boundary conditions. Similar to the experiment, the water cuts and the mixture velocities in the pipe were adjusted by regulating the velocity of each phase at the inlets. One velocity inlet, consisting a mixture of oil and water, was considered as an inlet for the experiment of Angeli and Hewitt [62]. The CFD model was initialized by assuming a single phase oil flow while the dispersed phase droplets entered to the physical domain from the corresponding inlets according to the experiment. Furthermore, the hydraulic diameter and turbulent intensity at the inlets were specified based on the pipe diameter and the fluid or flow properties. The turbulent intensity at the inlet was estimated as follows: A pressure outlet boundary condition was considered for the outlet boundary and the gauge pressure was considered zero, which is compatible with the separator pressure in the experiments.
A no-slip boundary condition was applied at the wall of the pipe. A symmetric boundary condition was adopted in the present study to reduce the computational. A grid sensitivity study was performed and a grid consist of 364,000 hexahedral cells was found to be suitable for model the experiments conducted by Elseth [82]. The same grid was scaled for modelling the experiment of Angeli and Hewitt [52]. GAMBIT 2.3.16 [83] was used to generate the geometry and the grid. Figure 1 shows the employed grid, as well as the boundary conditions in the present study.

Results and Discussions
Numerical simulations were performed based on the described models in the previous section. Four experimental tests were used to evaluate the accuracy of the CFD model in estimating the Sauter mean diameter and water distribution. The capability of the CFD model in predicting the in-situ droplet size was judged by making a comparison with the experimental data of Angeli and Hewitt [62] corresponding to a dilute dispersion of water in oil flow. Furthermore, the experimental data of Elseth [82] was used to evaluate the accuracy of the CFD model in predicting the distribution of the volume fraction of the dispersed water droplets across the pipe section. The physical properties of the fluids and pipes specifications for all the tests can be found in Table 1.
As mentioned, test 1 and test 2 were chosen as a benchmark for assessing the accuracy of the droplet size predictions. The input water cut of both test 1 and test 2 were 0.05,

Results and Discussions
Numerical simulations were performed based on the described models in the previous section. Four experimental tests were used to evaluate the accuracy of the CFD model in estimating the Sauter mean diameter and water distribution. The capability of the CFD model in predicting the in-situ droplet size was judged by making a comparison with the experimental data of Angeli and Hewitt [62] corresponding to a dilute dispersion of water in oil flow. Furthermore, the experimental data of Elseth [82] was used to evaluate the accuracy of the CFD model in predicting the distribution of the volume fraction of the dispersed water droplets across the pipe section. The physical properties of the fluids and pipes specifications for all the tests can be found in Table 1.
As mentioned, test 1 and test 2 were chosen as a benchmark for assessing the accuracy of the droplet size predictions. The input water cut of both test 1 and test 2 were 0.05, while the flow velocities were 1.25 m/s and 1.7 m/s, respectively. The minimum and maximum droplet sizes were set to 5 × 10 −5 m and 3.2 × 10 −3 m, respectively. Furthermore, the specified range of droplet sizes in the population balance model was divided into 10 intervals or classes (bins). It should be mentioned that in the experiment the endoscope window was placed at the core of the pipe and the maximum depth of view was limited to 4.5 mm upward. Hence, the area weighted average of the predicted Sauter mean diameter by CFD model in this area was determined to compare with the experimental measurement data. Table 2 shows the Sauter mean diameter for test 1 and test 2 as measured in the experiments and predicted by the CFD model. A comparison of the predicted and measured Sauter mean diameters using the present CFD model and the experiment indicates a reasonable accuracy of the CFD model. However, the CFD model slightly over predicted the Sauter mean diameter for both tests. Furthermore, the measurement of Sauter mean diameter at two different velocities suggests that an increase in the mixture velocity results in a decrease in the Sauter mean diameter. This trend is also observed in the predictions of the CFD model. Figure 2 shows the vertical distribution of Sauter mean diameter in whole pipe cross section for test 1 and test 2 as predicted by the CFD model. According to this figure, regardless of the magnitude of the droplet size, the distribution of the droplet size across the pipe section is very similar for both tests. As shown in Figure 2, the Sauter mean diameter at the top of the pipe is very small. As we move downward, the Sauter mean diameter significantly increases until reaching its maximum value at slightly below the center of the pipe. However, by further moving downward a gradual decrease in the droplet size is observed, except in the region close to the bottom wall where the droplet size suddenly increases. This trend generally agrees well with the experimental observations since due to the higher turbulence dissipation rate the breakup rate near the wall region is the highest. The observed sudden increase in the droplet size in flow region adjacent to the bottom wall might be attributed to the low accuracy of the adopted turbulence model in this study, which is assuming the isotropic turbulence while the turbulence near the wall is clearly anisotropic. Furthermore, the high volume fraction of the dispersed phase at the bottom of the pipe may increase the probability of coalescence, thereby, increasing the droplet size. Figure 3 shows the contours of Sauter mean diameter at the pipe cross section as predicted by using the CFD model.
The application of the used breakup and coalescence models in this study is limited to the inertial subrange of turbulence. Furthermore, the model presented by Brauner [61] is an extension of Hinze [63] model, which is also based on the particle breakage in the inertial subrange. Hence, the size of the dispersed droplets should be larger than the size of the smallest eddies that are recognized as the Kolmogorov length scale. Figure 4 shows the Kolmogorov length scale across the pipe cross section, which indicates the size of the smallest eddies. A comparison of the size of the predicted Sauter mean diameter in Figure 2 and the smallest eddies in Figure 4 clearly reveals that the dispersed droplets are larger than the size of the smallest eddies, thereby, more prone to the turbulent inertial subrange. Hence, the influence of viscos breakage can be neglected.    The application of the used breakup and coalescence models in this to the inertial subrange of turbulence. Furthermore, the model presented is an extension of Hinze [63] model, which is also based on the particle inertial subrange. Hence, the size of the dispersed droplets should be larg of the smallest eddies that are recognized as the Kolmogorov length scale. the Kolmogorov length scale across the pipe cross section, which indicate smallest eddies. A comparison of the size of the predicted Sauter mean dia 2 and the smallest eddies in Figure 4 clearly reveals that the dispersed dro than the size of the smallest eddies, thereby, more prone to the turbulent in Hence, the influence of viscos breakage can be neglected.  The application of the used breakup and coalescence models in this study is limited to the inertial subrange of turbulence. Furthermore, the model presented by Brauner [61] is an extension of Hinze [63] model, which is also based on the particle breakage in the inertial subrange. Hence, the size of the dispersed droplets should be larger than the size of the smallest eddies that are recognized as the Kolmogorov length scale. Figure 4 shows the Kolmogorov length scale across the pipe cross section, which indicates the size of the smallest eddies. A comparison of the size of the predicted Sauter mean diameter in Figure  2 and the smallest eddies in Figure 4 clearly reveals that the dispersed droplets are larger than the size of the smallest eddies, thereby, more prone to the turbulent inertial subrange. Hence, the influence of viscos breakage can be neglected.    Figure 5 depicts the vertical distribution of the dispersed water volume fraction from the developed CFD model in this work for test 1 and test 2 with the mixture velocities of 1.25 m/s and 1.7 m/s, respectively. As seen, a relatively homogeneous dispersion is predicted by the CFD model in which the water volume fraction increases as moving to the bottom of the pipe. This figure also reveals that for a constant input volume fraction of the dispersed phase an increase in the mixture velocity results in a more homogeneous dispersion. As mentioned earlier, the mesh sensitivity study was conducted to make sure that the obtained CFD results are independent of the mesh resolution. Figure 6 shows the distribution of water volume fraction across the pipe section for test 2 with the mixture velocity of 1.25 m/s and 5% water cut. As seen in this figure, by increasing the cell number beyond 364,000 the obtained numerical results would not change significantly. Figure 7 shows the contours of water volume fraction across the pipe section as obtained by the CFD model for test 1 and test 2.
dispersed phase an increase in the mixture velocity results in a more homogeneo persion. As mentioned earlier, the mesh sensitivity study was conducted to mak that the obtained CFD results are independent of the mesh resolution. Figure 6 sho distribution of water volume fraction across the pipe section for test 2 with the m velocity of 1.25 m/s and 5% water cut. As seen in this figure, by increasing the cell n beyond 364,000 the obtained numerical results would not change significantly. F shows the contours of water volume fraction across the pipe section as obtained CFD model for test 1 and test 2.   that the obtained CFD results are independent of the mesh resolution. Figure 6 sho distribution of water volume fraction across the pipe section for test 2 with the m velocity of 1.25 m/s and 5% water cut. As seen in this figure, by increasing the cell n beyond 364,000 the obtained numerical results would not change significantly. F shows the contours of water volume fraction across the pipe section as obtained CFD model for test 1 and test 2.   Numerical simulations of dispersed water in oil flow were also conducted for two different tests (test 3 and test 4) with available in-situ dispersed phase distributions. The experimental data were taken from Elseth [82]. The operating conditions of test 3 and test 4 have been reported in Table 1. Figure 8 shows the vertical distribution of Sauter mean diameter for test 3, in which the mixture velocity and water cut are 2 m/s and 0.1, respectively. A similar distribution of Sauter mean diameter as observed in Figure 2 is perceived in this figure. The minimum and maximum droplet sizes of 1 × 10 −4 m and 5.5 × 10 −3 m were chosen to specify the range of the droplet size distribution in population balance model, which was divided into 10 intervals or classes (bins). As seen, the Sauter mean diameter at the top of the pipe is very small and as we move downward it significantly increases, following a decrease after reaching its maximum value. However, similar to the observed phenomenon in Figure 2, a sudden increase in the droplet size is observed closed to the wall. Due to the limited experimental data for the in-situ distribution of the droplet size, the accuracy of the predicted Sauter mean diameter cannot be assessed. Figure 9 shows the contours of Sauter mean diameter across the pipe section by using the CFD model. Numerical simulations of dispersed water in oil flow were also conducted for two different tests (test 3 and test 4) with available in-situ dispersed phase distributions. The experimental data were taken from Elseth [82]. The operating conditions of test 3 and test 4 have been reported in Table 1. Figure 8 shows the vertical distribution of Sauter mean diameter for test 3, in which the mixture velocity and water cut are 2 m/s and 0.1, respectively. A similar distribution of Sauter mean diameter as observed in Figure 2 is perceived in this figure. The minimum and maximum droplet sizes of 1 × 10 −4 m and 5.5 × 10 −3 m were chosen to specify the range of the droplet size distribution in population balance model, which was divided into 10 intervals or classes (bins). As seen, the Sauter mean diameter at the top of the pipe is very small and as we move downward it significantly increases, following a decrease after reaching its maximum value. However, similar to the observed phenomenon in Figure 2, a sudden increase in the droplet size is observed closed to the wall. Due to the limited experimental data for the in-situ distribution of the droplet size, the accuracy of the predicted Sauter mean diameter cannot be assessed. Figure 9 shows the contours of Sauter mean diameter across the pipe section by using the CFD model.   Figure 10 illustrates a comparison of the vertical distribution of water volume fraction in pipe cross section as obtained by the combined CFD-PBM model and the experimental measurement using a gamma densitometer. According to this figure, the CFD model could reasonably predict the water distribution across the pipe section. A small deviation observed at the bottom of the pipe could be attributed to the deficiency of the adopted closure laws for turbulent dispersion or the lift forces. The present study employed the turbulent dispersion force suggested by Lopez de Bertodano [72] in which the turbulent dispersion coefficient was a constant. Using a constant dispersion coefficient throughout the computational domain is a common practice in most of the studies regarding the CFD modelling of the dispersed flows [84,85]. However, from the physical point of view the turbulent dispersion coefficient is a function of the Stokes number which is   Figure 10 illustrates a comparison of the vertical distribution of water volume fraction in pipe cross section as obtained by the combined CFD-PBM model and the experimental measurement using a gamma densitometer. According to this figure, the CFD model could reasonably predict the water distribution across the pipe section. A small deviation observed at the bottom of the pipe could be attributed to the deficiency of the adopted closure laws for turbulent dispersion or the lift forces. The present study employed the turbulent dispersion force suggested by Lopez de Bertodano [72] in which the turbulent dispersion coefficient was a constant. Using a constant dispersion coefficient throughout the computational domain is a common practice in most of the studies regarding the CFD modelling of the dispersed flows [84,85]. However, from the physical point of view the turbulent dispersion coefficient is a function of the Stokes number which is  Figure 10 illustrates a comparison of the vertical distribution of water volume fraction in pipe cross section as obtained by the combined CFD-PBM model and the experimental measurement using a gamma densitometer. According to this figure, the CFD model could reasonably predict the water distribution across the pipe section. A small deviation observed at the bottom of the pipe could be attributed to the deficiency of the adopted closure laws for turbulent dispersion or the lift forces. The present study employed the turbulent dispersion force suggested by Lopez de Bertodano [72] in which the turbulent dispersion coefficient was a constant. Using a constant dispersion coefficient throughout the computational domain is a common practice in most of the studies regarding the CFD modelling of the dispersed flows [84,85]. However, from the physical point of view the turbulent dispersion coefficient is a function of the Stokes number which is defined as the ratio of a response time of the dispersed particles to the eddies of the continuous phase [71,72]. Hence, using a variable turbulent dispersion coefficient could potentially enhance the accuracy of the CFD solution. Apart from that, using other turbulent dispersion force models such as the model proposed by Burns et al. [86] may result in better predictions. defined as the ratio of a response time of the dispersed particles to the eddies of the continuous phase [71,72]. Hence, using a variable turbulent dispersion coefficient could potentially enhance the accuracy of the CFD solution. Apart from that, using other turbulent dispersion force models such as the model proposed by Burns et al. [86] may result in better predictions.  Figure 11 shows the contours of water volume fraction distribution across the pipe section as obtained by the CFD model for test 3. Figure 11. Distribution of water volume fraction across pipe section for test 3 (mixture velocity = 2 m/s, water cut = 0.1). Figure 12 shows the contours of the in-situ water volume fraction distribution for test 4 as obtained by the CFD-PBM model. The mixture velocity and water cut for this test are 1 m/s and 0.1, respectively. In the CFD model, the minimum and maximum droplet sizes  Figure 11 shows the contours of water volume fraction distribution across the pipe section as obtained by the CFD model for test 3. defined as the ratio of a response time of the dispersed particles to the eddies of the continuous phase [71,72]. Hence, using a variable turbulent dispersion coefficient could potentially enhance the accuracy of the CFD solution. Apart from that, using other turbulent dispersion force models such as the model proposed by Burns et al. [86] may result in better predictions.  Figure 11 shows the contours of water volume fraction distribution across the pipe section as obtained by the CFD model for test 3. Figure 11. Distribution of water volume fraction across pipe section for test 3 (mixture velocity = 2 m/s, water cut = 0.1). Figure 12 shows the contours of the in-situ water volume fraction distribution for test 4 as obtained by the CFD-PBM model. The mixture velocity and water cut for this test are 1 m/s and 0.1, respectively. In the CFD model, the minimum and maximum droplet sizes Figure 11. Distribution of water volume fraction across pipe section for test 3 (mixture velocity = 2 m/s, water cut = 0.1). Figure 12 shows the contours of the in-situ water volume fraction distribution for test 4 as obtained by the CFD-PBM model. The mixture velocity and water cut for this test are 1 m/s and 0.1, respectively. In the CFD model, the minimum and maximum droplet sizes of 1 × 10 −4 m and 6.4 × 10 −3 m were chosen to specify the range of the droplet size distribution. Furthermore, the droplet size distribution range was divided into 11 classes.
Processes 2021, 9,1345 16 of 22 of 1 × 10 −4 m and 6.4 × 10 −3 m were chosen to specify the range of the droplet size distribution. Furthermore, the droplet size distribution range was divided into 11 classes.  Figure 13 depicts a comparison of the predicted vertical distribution of water volume fraction by using CFD model and the experimental measurement using gamma densitometer. According to this figure the CFD model could reasonably predict the distribution of dispersed water phase across the pipe section. As seen in this figure, the experimental measurement, as well as the CFD model revealed the occurrence of a dense water dispersion in oil at the bottom of the pipe. The models developed for the lift force can be applicable for relatively dilute systems where a dense packing of the dispersed phase does not occur [80,87]. Hence, using such models for dispersed phase volume fractions higher than 0.5 is questionable and likely to give physically misleading results. Furthermore, when the volume fraction of the dispersed phase increases, the lift coefficient tends to approach zero [87,88]. Apart from these, the influence of the lift force for regions close to the center of the pipe is negligible. Thus, the influence of lift force has been neglected for this case as the droplet volume fraction at the bottom of the pipe exceeds 0.5. Figure 14 shows the vertical distribution of Sauter mean diameter for test 4. It is observed that the distribution of Sauter mean diameter is quite similar to the predicted distribution for test 3. However, due to the lower mixture velocity the probability of breakup decreases, while the coalescence probability increases. Hence, the average Sauter mean diameter for test 4 is larger than that for test 3. Figure 15 shows the contours of the predicted Sauter mean diameter for test 4.  Figure 13 depicts a comparison of the predicted vertical distribution of water volume fraction by using CFD model and the experimental measurement using gamma densitometer. According to this figure the CFD model could reasonably predict the distribution of dispersed water phase across the pipe section. As seen in this figure, the experimental measurement, as well as the CFD model revealed the occurrence of a dense water dispersion in oil at the bottom of the pipe. The models developed for the lift force can be applicable for relatively dilute systems where a dense packing of the dispersed phase does not occur [80,87]. Hence, using such models for dispersed phase volume fractions higher than 0.5 is questionable and likely to give physically misleading results. Furthermore, when the volume fraction of the dispersed phase increases, the lift coefficient tends to approach zero [87,88]. Apart from these, the influence of the lift force for regions close to the center of the pipe is negligible. Thus, the influence of lift force has been neglected for this case as the droplet volume fraction at the bottom of the pipe exceeds 0.5.   Figure 14 shows the vertical distribution of Sauter mean diameter for test 4. It is observed that the distribution of Sauter mean diameter is quite similar to the predicted distribution for test 3. However, due to the lower mixture velocity the probability of breakup decreases, while the coalescence probability increases. Hence, the average Sauter mean diameter for test 4 is larger than that for test 3. Figure 15 shows the contours of the predicted Sauter mean diameter for test 4.

Conclusions
A three-dimensional numerical modelling of turbulent dispersed water in oil flow pattern through a horizontal pipe was carried out in the present study. Numerical modelling was conducted using a two-fluid Eulerian approach. Four interfacial forces namely, the drag, virtual mass, lift and turbulent dispersion were considered to close the two-fluid model. The standard K-epsilon model was employed to model the turbulent effects. Furthermore, the droplet size distribution was modelled by using a population balance model based on the method of class in which the breakup and coalescence events were modelled using Lehr and Luo models, respectively.

Conclusions
A three-dimensional numerical modelling of turbulent dispersed water in oil flow pattern through a horizontal pipe was carried out in the present study. Numerical modelling was conducted using a two-fluid Eulerian approach. Four interfacial forces namely, the drag, virtual mass, lift and turbulent dispersion were considered to close the two-fluid model. The standard K-epsilon model was employed to model the turbulent effects. Furthermore, the droplet size distribution was modelled by using a population balance model based on the method of class in which the breakup and coalescence events were modelled using Lehr and Luo models, respectively.
A comparison of the predicted Sauter mean diameter by using the CFD model with the corresponding experimental data indicates that a combination of the Luo coalescence model and the Lehr breakup model can reasonably predict the Sauther mean diameter.
Numerical results of the in-situ water fraction distributions were also compared with the available experimental data. According to the obtained CFD results, a combination of the adopted interfacial force models and the predicted Sauter mean diameter using the population balance model resulted in a reasonable prediction of the dispersed water distribution.
The obtained results of the present study revealed the capability of the employed CFD model for predicting the main features of complex dispersed oil-water flows in pipeline. However, the adopted closure laws for the interfacial forces require further study and calibration. Moreover, a comparative study of different breakup and coalescence models is essential to identify the most appropriate kernels for modelling the developing dispersed oil-water flows.