The local dynamics of spatially developing liquid-liquid dispersed flows at low superficial velocities, ranging from 0.2 to 0.8 m s^{−1}, are investigated. The dispersions are generated with an in-line static mixer. Detailed measurements with laser-based diagnostic tools are conducted at two axial pipe locations downstream of the mixer, namely, at 15 and 135 equivalent pipe diameters. Different flow patterns are recorded, and their development along the streamwise direction is shown to depend on the initial size and concentration of the drops as well as the mixture velocity. The drop size is accurately predicted by an empirical formula. The variations in drop concentration over the pipe cross-section along the pipe result in local changes of the physical properties of the mixture and consequently in asymmetrical velocity profiles, with the maxima of the velocity located in the drop-free region. Computational fluid dynamics simulations based on a mixture approach predict the experimental results close to the experimental uncertainties for the majority of the cases. The simulation results reveal that gravity and lift forces, as well as shear-induced diffusion are the most important mechanisms affecting the drop migration. It is found that the drops behave as suspensions of rigid spheres for the conditions investigated, despite the deformation effects, which are found experimentally to be stronger at the densely packed region.

## I. INTRODUCTION

The flow and separation characteristics of unstable liquid-liquid dispersions in pipes involve complex phenomena of fundamental theoretical interest.^{1} Apart from nuclear cooling systems, separators, and two-phase reactors, liquid-liquid dispersions are particularly important for the oil and gas industry where water is often present in oil transportation pipelines, resulting either from condensation of saturated gas or directly produced from the reservoirs as they mature. Depending on the conditions, different flow patterns can form which in turn affect the flow characteristics of the mixture.^{2} Dispersions often form at pieces of equipment such as choke valves or bends, and they tend to separate further downstream. Recently, there have been efforts to develop in-line separators in the pipe and avoid transporting water over long distances, especially in subsea pipelines. Understanding the interactions between the spatial distribution of the two phases and the local flow characteristics is crucial for the design and operation of these facilities.

Toward this direction, a need for high-resolution and detailed measurements of the spatial evolution characteristics of segregating liquid-liquid dispersions has emerged, which can shed light into the underlying physical mechanisms behind the motion of drops in dispersed pipe flows and into the macroscopic behavior of the mixture. Such measurements can also set the foundations for the development and validation of advanced numerical models. While nuclear magnetic resonance (NMR) imaging has mainly been used to simultaneously record the velocity and concentration profiles in suspension flows,^{3} in liquid-liquid systems, optical measurements with refractive index matched fluids are increasingly implemented.^{4–7} Compared to NMR, the optical measurements have robustness and versatility at reduced costs.^{3}

Studies of dispersed Poiseuille flows, where drop deformability and polydispersity are neglected, have shown that the available rheological models developed for liquid-solid suspensions are capable of predicting integral parameters, such as the overall pressure drop in liquid-liquid pipe flows.^{8–11} These relatively recent findings are promising because one might consider investigating the local hydrodynamic behavior of axially developing liquid-liquid dispersions using the well-established models developed for suspensions.

In contrast to liquid-liquid dispersed flows, a relatively large number of works have explored the dynamics of settling in pressure-driven solid-liquid suspension flows.^{12,13} Asymmetrical velocity profiles were observed when the two phases separated.^{13,14} These observations made for suspensions were mainly attributed to local mixture viscosity changes, caused by high concentration gradients of the suspended phase. To predict the dynamics of these flows, simplified phenomenological modeling approaches have been developed. One example is the work of Zhang and Acrivos^{12} who extended the model originally developed by Phillips *et al.*^{15} for unidirectional Couette flows to pipe flows, where settling of the mixture takes place and a bed of the suspended phase forms. By assuming that gravity and particle interactions due to shear and concentration gradients lead to diffusion, this model captured the cross-stream particles motion reasonably well. This modeling approach has not been explored for liquid-liquid systems.

In this article, a combined experimental and computational approach is used to (i) investigate the links between the separation properties and the flow characteristics and (ii) to analyze the mechanisms affecting the motion of the drops in concentrated liquid-liquid dispersions that are generated by a static mixer in a horizontal pipe. The numerical framework is a first attempt to apply a mixture modelling approach developed for suspensions to liquid-liquid dispersions with the aid of Computational Fluid Dynamics (CFD) simulations. As such, another aim of this paper is to (iii) investigate the applicability as well as limitations in using this mixture modeling approach for dispersions in laminar flow. Previous modeling strategies for liquid-liquid flows have implemented both Eulerian-Eulerian^{16,17} and mixture modeling frameworks.^{18} There are, however, two key differences between those models and the one presented herein. On the one hand, the contribution of shear-induced diffusion was neglected in the previous works as the flows modeled were turbulent, and on the other hand, the viscosity closures employed in those works are different to the one employed here.

Experimentally, a combined high-speed planar laser-induced fluorescence (PLIF) and particle image velocimetry (PIV) method is used to simultaneously obtain spatiotemporally resolved information on the drop size distribution, drop concentration, and velocity of the continuous phase. The refractive indices of the fluids are matched to enable the studies at high dispersed phase fractions. A wide range of input dispersed phase volume fractions (10%–50%) and mixture velocities (0.2–0.8 m s^{−1}) are studied, which enable comparisons of drop size characteristics in both oil and water continuous phases. The flow properties of separating oil-in-water dispersions are further considered, where the drops segregate near the pipe wall and are no longer homogeneously distributed along the vertical direction of the pipe.

Using the model, the velocity profiles and development length of the dispersions can be predicted. The results of the simulations are used to understand how the presence of large drop concentration gradients, caused by shear-induced diffusion, lift and gravity, can influence the local rheology and consequently the fluid dynamics of liquid-liquid flows. Importantly, using the model, an approximation for the development length of the dispersion can be predicted, which is important for designing pipelines. Finally, the limitations of the approach in modeling liquid-liquid dispersions are explored by analyzing the deformation effects of the drops along the pipe vertical cross-section.

## II. EXPERIMENTAL METHODOLOGY

### A. Facility and fluids

The experiments were carried out in a large-scale flow facility consisting of two horizontal acrylic 4 m long pipe sections of inner diameter *D* = 26 mm, connected through a U-bend. All the experiments were conducted in the front section of the loop and considerably away from the U-bend to avoid any end-effects. The liquids were introduced in the pipe section cocurrently via a Y-junction with the water coming from the bottom pipe and the oil from the top one, and they then passed through a 6-element helical static mixer. Details of the test section can be found in the work of Voulgaropoulos and Angeli^{7} and Voulgaropoulos *et al.*^{19} To be able to carry out optical measurements at high dispersed phase fractions, refractive index-matching liquids should be used.^{6,7,10,20} In the experiments, a low viscosity (*μ*_{o} = 0.0046 Pa s) and high density (*ρ*_{o} = 913 kg m^{−3}) silicone oil, along with a 52% w/w water/glycerol mixture (*ρ*_{$w$} = 1146 kg m^{−3}, *μ*_{$w$} = 0.0084 Pa s), were chosen as test fluids with a refractive index equal to 1.396. Henceforth, the two phases are abbreviated to “oil” and “water,” respectively. The interfacial tension was *σ* = 0.0316 N m^{−1}. While the temperature is not actively controlled in the system, it was periodically recorded and was found to be about 20 °C; the properties of the fluids were measured at this temperature. The maximum temperature change recorded was below Δ*T* < 5 °C, which translates to a modification of the refractive index ratio by 0.2%. This change was negligible as there was no noticeable optical distortion.

Pressure gradient measurements were also performed with a differential pressure transducer (ABB 266MST), which had the pressure taps located at *x* = 3 and *x* = 3.5 m downstream of the inlet. These two axial locations were selected to avoid both entrance effects from the inlet and back-flow effects from the U-bend. Measurements were taken over a 60 s period at a sampling frequency of 1 kHz for each flow condition.

For the optical measurements, a diode-pumped solid-state continuous laser (Laserglow Technologies), emitting at 532 nm, was placed below the pipe, while a thin laser sheet of approximately 1 mm thickness was generated and focused in the middle of the pipe with appropriate optical lenses. A high-speed camera (Phantom v1212 by Vision Research) was placed perpendicular to the laser sheet in front of the test section. Simultaneous time-resolved PLIF and PIV measurements were conducted. A small amount (∼0.02 ppm) of Rhodamine 6G fluorescent dye (Sigma-Aldrich) was added to the water to distinguish between the two phases and improve the tracking of the interfaces. Quasineutrally buoyant Rhodamine B coated PMMA tracer particles (Sigma-Aldrich) with a size of 1–20 *μ*m were added to the aqueous phase.

Acquisition frequencies ranging from 1 to 4 kHz were used to acquire time-averaged results over 4 s. A high-pass filter cutting below 580 nm was mounted on the lens of the camera to avoid spurious light and reflections of the laser and isolate the emitted light from the dye and tracers. The measurements were conducted at two locations downstream of the inlet mixer, at *x*^{+} = 15 and 135, where *x*^{+} is the normalized distance in the streamwise (axial) direction *x* with the pipe diameter *D*. An acrylic visualization box filled with pure glycerol (which matches relatively closely the refractive index of the pipe material) was placed at these two locations to avoid optical distortions caused by the pipe curvature.

### B. Analysis and uncertainties

The image analysis and the steps followed to reduce measurement errors are explained in detail in the work of Voulgaropoulos and Angeli^{7} and Voulgaropoulos,^{21} but a brief description is given here for completeness. The image analysis methodology followed is summarized in Fig. 1. The raw images obtained [Fig. 1(a)] were adjusted with an adaptive histogram equalization filter to increase their clarity, while a median filter and background removal were implemented to isolate the drops [Fig. 1(b)] and the PIV tracer particles [Fig. 1(c)], respectively. From the local median values of the images, adaptive image thresholding was used to binarize the images for both the drops and the tracers, as shown in Figs. 1(d) and 1(e), respectively. From the binarized images, the drop size and local phase fraction [Fig. 1(f)] as well as the velocity profiles in the aqueous phase [Fig. 1(g)] were calculated.

To calculate the velocity in the aqueous phase with PIV, an iterative discrete Fourier transformation (DFT) algorithm was used,^{21} which calculates the correlation matrix in the frequency domain. The presence of drops in the correlation windows can affect the cross-correlation as the number of tracer particles decreased when the local concentration of the drops increased. For this reason, a three-pass DFT was used with an initial interrogation window of 64 × 64 pixels^{2}. The area of this interrogation window is larger than the typical drop areas observed in the experiments so that there are enough particles in each window to maintain signal to noise ratios above 2, and primary peak ratios, namely, the ratio between the primary correlation peak and the second tallest peak of the correlation map, above 1.3. For most flow conditions, the drops have diameters smaller than 0.1*D*, while the interrogation window during the first iteration was approximately equal to ∼0.15*D*. The relative uncertainties of the streamwise velocity *u* computed from the PIV are of the order of 1%, but close to the pipe walls or at regions of dense-packed drops, the uncertainties were estimated at 30%.

The relative uncertainty of the drop size algorithm was found equal to 10%,^{7} while the smallest drop that could be captured was equal to 4 pixels or 0.24 mm. Spatially resolved *in situ* phase fraction measurements were conducted and their relative uncertainty reached up to 20%, stemming primarily from the adjustment of the sensitivity of the adaptive binarization method. Similar issues have been tackled by Morgan *et al.*^{5} The acquisition time of 4 s in conjunction with the high sampling frequency ensured that a large number of samples were obtained for each set of conditions and that statistical convergence was reached. Over 4000 velocity and concentration profiles were obtained and time-averaged, while at least 600 individual drops (enough to capture very large drops of very low probability^{22,23}) were measured for each flow condition. Any changes to the profiles and distributions for larger samples were found to be negligible, while the sampling errors based on the standard deviations were calculated below 1%.

### C. Conditions

The mixture velocity is equal to $um,in=Qc+Qd/A$, where *A* is the pipe cross-sectional area, while *Q*_{c} and *Q*_{d} are the input volumetric flow rates of the continuous and the dispersed phases, respectively. The input dispersed phase volume fraction is equal to $\epsilon d,in=Qd/Qc+Qd$. Several flow patterns were generated from the static mixer, but since the focus was on dispersed flows, mixture velocities from 0.1 to 1.0 m s^{−1} were used, while the input dispersed phase volume fractions ranged from 0.1 to 0.5 for both oil and water continuous dispersions.

A Reynolds number based on the mixture velocity, the pipe diameter, and the mixture properties is defined as

The mixture density is computed as

where *ε*_{c} ≡ 1 − *ε*_{d} is the volume fraction of the continuous phase, *ε*_{d} is the volume fraction of the dispersed phase, and *ρ*_{c} and *ρ*_{d} denote the density of the continuous and dispersed phases, respectively. The mixture viscosity, *μ*_{m}, can be calculated by a Krieger-Dougherty type of equation,^{24} which was developed to describe the viscosity of a suspension of rigid spheres as

where *μ*_{c} is the viscosity of the continuous phase and *ε*_{max} = 0.74 is the maximum packing. This maximum packing value is higher than the 0.64 normally considered for randomly packed monodisperse spheres and was found by Pouplin *et al.*^{10} to describe well the polydispersity effects expected at concentrated liquid-liquid dispersed systems. For neutrally buoyant particulate pipe flows, the critical Re_{m} for transitioning to turbulence was found by Matas *et al.*^{25} to depend, apart from the mixture Reynolds number [Eq. (1)], on the drop to pipe diameter ratio and *ε*_{d}. Based on their findings and for the drop to pipe diameter ratios and drop concentrations investigated in this study, laminar flows are expected, with Re_{m} varying from 100 to 2000 for 36 out of the 40 flow conditions.

## III. RESULTS

### A. Flow patterns and pressure drop

The two immiscible liquids enter the static mixer in a stratified (ST) manner, with the lighter oil at the top and the heavier water/glycerol at the bottom of the pipe. In all cases studied, dispersions were formed at the exit of the mixer, which then segregated further downstream. Figure 2 illustrates typical high-speed PLIF images obtained for representative flow conditions at the two axial locations studied, namely, at *x*^{+} = 15 (top images) and 135 (bottom images). Figures 2(a)–2(c) show an oil in water dispersion, *o*/$w$, while Figs. 2(d)–2(f) show a water in oil dispersion $w$/*o*. The water contains the fluorescent dye and the PIV tracer particles and appears bright in the images, while the oil appears dark.

For the water continuous cases [Figs. 2(a)–2(c)], the dispersions at the outlet of the mixer at *x*^{+} = 15 appear uniform. As the dispersions evolve, the oil drops start to separate, and a clear aqueous phase layer forms at the bottom of the pipe. This is less pronounced in the case of *u*_{m,in} = 0.58 m s^{−1} and *ε*_{o,in} = 0.12 [Fig. 2(a)], where the initial dispersion contains some very small drops which are still present close to the bottom of the pipe at *x*^{+} = 135. Another possible reason for the location of these drops can be re-entrainment, which has been reported for similar conditions.^{4} Apart from the stratification of the oil fraction, a stratification of drop size also occurs along the pipe, with larger drops located near the upper part of the pipe.

Differences are observed for the oil continuous $w$/*o* cases of Figs. 2(d)–2(f), obtained at the same mixture velocities and dispersed phase volume fractions. Larger drops are generated by the mixer in the oil continuous dispersions compared to the water continuous ones. For the lower mixture velocity case [Fig. 2(f) top], the drops have already started to separate at *x*^{+} = 15, and a clear oil layer is visible at the top of the pipe. In all cases, at the downstream location *x*^{+} = 135, a clear water layer has developed at the bottom of the pipe. This suggests that the settled drops coalesce to form the continuous layer.

From the PLIF observations, flow pattern maps are developed in nondimensional form and presented in Figs. 3(a) and 3(b) for the two axial locations, 15*D* and 135*D*, respectively. The Froude number is taken equal to $Fr=um,in/gD0.5$. In the maps, the transition boundaries (dashed lines) for the flows without the static mixer are also shown for reference. For both axial locations, all the patterns observed when the mixer is present are dispersed, while at the same conditions without the mixer, the flow is stratified with no drops present. The patterns generated by the mixer are not stable and should separate downstream provided the test section is long enough.

As can be seen from Fig. 3 at the initial axial location *x*^{+} = 15, homogeneously dispersed flows (either oil or water continuous) can be achieved for Fr ∼ 1. At the same conditions, at *x*^{+} = 135 downstream of the mixer [Fig. 3(b)], the drops have segregated and a clear drop-free layer of the continuous phase forms. In the oil continuous dispersions, however, a layer of the initial dispersed water phase has already formed, indicating a flow pattern transition to dual continuous (*DC*) flow. The analysis provided later in the manuscript (Sec. III D) shows that the dispersed flows of Fig. 3(b) are developed at *x*^{+} = 135. However, the same statement cannot be argued for the *DC* flows of Fig. 3(b) since the *DC* are more susceptible to separation and coalescence effects as is indicated by the drop size measurements, which are presented in Sec. III B.

Figure 4 illustrates the pressure gradients, Δ*P*/Δ*x*, recorded with the differential pressure transducer at the end of the test section, i.e., ∼125*D*, for the same 40 flow conditions as in Fig. 3. Pressure drop differences were recorded between the different flow patterns. The data exhibit lower Δ*P*/Δ*x* values for the *DC* than the *o*/$w$*&*$w$ flows. The higher viscosity of the water phase combined with the mixture effects of the dense-packed layer at the top of the pipe could be the possible contributing factors. The phase inversion point was found for this work at *ε*_{o} ≈ 0.5 ± 0.05^{21} and is highlighted with the shaded area in Fig. 3(b). As can be seen, pressure drops for these conditions are high, which can be attributed to the increased phase fraction of the drops, and according to Eq. (3), to an exponential increase of the mixture viscosity.

### B. Drop size variations

To characterize the drop sizes of the dispersions in the pipe the Sauter mean diameter is calculated as follows:

where *N*_{d} is the number of drops, while *d*_{i} is the diameter (chord length) of drop *i*. It must be noted that the PLIF drop size measurements of Fig. 1 can only measure the drop chord lengths. The drop chord lengths are used in Eq. (4) and below, while they are expected to be slightly smaller than the actual drop diameters. Also, in the present work a proportionality is found between *d*_{32} and *d*_{max} with R^{2} ≃ 0.97 in agreement with the literature^{26} despite the fact that the PLIF measurements provide chord lengths instead of diameters.

The average drop size of the dispersions measured at the exit of the mixer *x*^{+} = 15 can be compared against the correlation $d32\u221d\rho c/\sigma \u22123/5\u0113\u22122/5$ suggested for static mixers.^{27} The mean rate of energy dissipation *ē* is related to the frictional pressure drop in the mixer^{27} as $\u0113=2um,in3f/D\rho m/\rho c1\u2212\epsilon d,in$, where *f* is the friction factor, which can be calculated from the Blasius equation. For dilute liquid-liquid flows, the equation simplifies to

where $We=\rho cum,in2D/\sigma $ is the Weber number and $d32+=d32/D$.

For dense dispersions with *ε*_{d,in} > 0.10, the following equation applies:^{28}

A comparison between the normalized Sauter mean diameters, $d32+$, measured experimentally and those predicted from Eq. (6), for a fitted proportionality constant of the order of 1, is shown in Fig. 5 for the initial axial measuring location downstream of the static mixer, *x*^{+} = 15. Over 30 conditions were tested, and good agreement with Eq. (6) is found for the cases where $d32+<0.1$ for both the oil in water (*o*/$w$) and the water in oil ($w$/*o*) dispersions. The results suggest that the drop generation in the mixer follows the theory of turbulent breakup. As the flow develops, turbulence dampens and coalescence becomes the dominating mechanism,^{7,29} which results in deviations between the experimental and predicted drop sizes (not shown here).

The time and cross-sectional-averaged drop size distributions of both the *o*/$w$ and $w$/*o* dispersions and their change with axial location are plotted as probability histograms P(*d*) in Fig. 6, for the same cases shown in Fig. 2. A constant bin size equal to 0.25 mm is used, which is found to describe accurately the trends for most of the conditions studied and illustrates in enough detail the characteristics of each distribution. We do not expect drops in the smallest bin as it overlaps with the resolution of the drop size measurement technique (i.e., 0.24 mm). The corresponding log-normal probability density functions PDF(*d*) for both types of dispersions are also plotted. They are computed as

where *m*_{d} is the mean and *s*_{d} is the standard deviation of the distributions. The log-normal density distributions are found to represent well the majority of cases.

The distributions for the *o*/$w$ dispersions, illustrated in Figs. 6(a)–6(c), are narrow with most of the drops being smaller than 1 mm. As the input concentration of the drops in the pipe increases, so does their size along with the probability to find large drops [cf. Figs. 6(a) and 6(b)]. When the mixture velocity decreases, the distribution becomes wider [cf. Figs. 6(b) and 6(c)], with a significant decrease in the number of small drops. The $w$/*o* dispersions are shown in Figs. 6(d)–6(f). The trends are similar, but the distributions are, in general, wider than in the *o*/$w$ dispersions. The difference in the distributions is a result of the lower energy dissipation in the static mixer when the oil is continuous. The oil has a lower viscosity than the water mixture, and the oil continuous dispersions have a lower friction factor and hence larger drop sizes [Eq. (5)] compared to the water continuous ones. The $w$/*o* distributions start deviating from the log-normal PDFs, especially at *x*^{+} = 135, where coalescence effects are more pronounced. This deviation is caused by the presence of very large drops *d* > 2 mm.

The variations of the average Sauter mean drop diameters along the vertical pipe directions for the same conditions are presented for both axial locations in Figs. 7(a)–7(c) and 7(d)–7(f) for the *o*/$w$ and $w$/*o* dispersions, respectively. To calculate the averages, the vertical distance is split in 10 equal horizontal segments of Δ*y*^{+} = 2.6 mm, and *d*_{32} for each segment is computed. The number of drops averaged in time in each segment is *N*_{d} > 200. For the *o*/$w$ dispersions in Figs. 7(a)–7(c), the drop size is almost uniform along the vertical direction at *x*^{+} = 15. There is a change of the profile at *x*^{+} = 135, especially for the high volume fraction and low mixture velocity case [Fig. 7(c)], which reflects the settling of the drops seen in the PLIF images of Fig. 2. Interestingly, a slight reduction in the drop size is recorded near the top in all cases of *o*/$w$, which could be attributed to lift forces close to the pipe walls.

Figures 7(d)–7(f) show the $w$/*o* cases. The dispersions are relatively homogeneous at the initial axial location for the high mixture velocity cases. In the case of Fig. 7(f), the segregation of the drops has already started and the mean drop sizes are larger in the lower part of the pipe, while a continuous layer of the aqueous phase forms further downstream. The drop size is found to increase close to the formed interface, and high gradients in the average drop size profiles occur, as also shown in the PLIF images of Fig. 2(f).

### C. Flow characteristics of the mixtures

The *o*/$w$ dispersions will be considered in more detail here. For these mixtures, the oil remains the dispersed phase along the length of the test section and only segregation of the dispersed drops occurs. The results can then be compared with the modeling approach detailed in Appendix, which is valid for suspension flows.

To predict all of the possible flow patterns, a universal modeling approach, which considers the point-wise fluid dynamics of both phases, would be needed. Only fundamental modeling approaches such as Direct Numerical Simulations (DNS) are able to do this. These approaches, however, are extremely demanding for industrial applications. For the latter, therefore, modelers find it more convenient to use simplified modeling approaches, such as the Eulerian-Eulerian (multifluid) approach. This approach uses a formal mathematical procedure of averaging and is able to predict the evolution of the observable parameters (such as volume fractions) directly. By doing so, the dynamic description of the system, compared to those obtained using fundamental approaches, such as DNS, is considerably reduced, but the information that usually is of interest in engineering applications can be captured with sufficient detail.

For a dispersed liquid-liquid isothermal mixture of monodisperse drops, the multifluid model solves four balance equations, namely, one mass and one linear momentum balance equation for each phase. A simplified, but very practical version of the multifluid approach is called the mixture modeling approach. In it, the mixture is considered as an effective fluid, and consequently only one continuity and one linear momentum balance equation are solved for it. Hence, a velocity field and a pressure field are obtained for the mixture as a whole. To calculate the volume fraction of each phase, the mass balance equation for one of the phases should be solved; this is usually done for the dispersed phase.

In the simulations, it is assumed that drops are monodisperse and their size does not change with axial location, i.e., no breakup or coalescence takes place. The experiments showed that the Sauter mean diameters for the *o*/$w$ dispersions only varied by a maximum of 200 *μ*m between the two measurement locations. With these sizes, the drop Stokes number is well below unity for all conditions investigated, allowing the use of the mixture model in the numerical simulations.

The flow conditions of *o*/$w$ dispersions investigated both experimentally and numerically are in the range of *u*_{m,in} = 0.40–0.60 m s^{−1}, and the input oil volume fractions varied from *ε*_{o,in} = 0.10–0.40. For these flow conditions, homogeneous *o*/$w$ dispersions were formed at the inlet, while at *x*^{+} = 135, segregation of the drops occurred, resulting in a drop-free water layer at the bottom of the pipe. PLIF images at the two axial locations for three typical cases were already illustrated in Figs. 2(a)–2(c). The Reynolds number ranges from Re_{m} = 500–1462, and thus turbulence can be ignored in the simulations.

The effects of gravity and shear-induced diffusion on the motion of the drops and the development of the concentration profiles along the pipe are discussed here. The vertical profiles of the time-averaged *in situ* oil volume fractions [⟨*ε*_{o}(*y*)⟩] (from the PLIF images) divided by the respective input oil volume fractions *ε*_{o,in} at *x*^{+} = 15 are plotted in Fig. 8 for four typical conditions. At the initial measuring location, the ratios of ⟨*ε*_{o}(*y*)⟩/*ε*_{o,in} are very close to unity, and the dispersion is almost homogeneous, particularly in the middle of the pipe. This justifies the assumption in the simulations that the dispersions are initially homogeneously distributed.

For the same four conditions, the experimental normalized vertical profiles of the local oil volume fractions are compared against the CFD simulations at *x*^{+} = 135 in Fig. 9. In the simulations, it is assumed that the drops are uniformly distributed at the pipe inlet (as shown in Fig. 7) and have a size equal to the arithmetic mean of the experimental overall *d*_{32} values measured at 15*D* and 135*D*. The simulated profiles agree reasonably well with the experimental results, which also illustrate asymmetric profiles, where drops have accumulated at the upper part of the pipe leading to a clear water layer at the bottom (*o*/$w$ *and* $w$ pattern). These findings are in agreement with the visualizations of Fig. 2. Similar profiles were observed for relatively higher mixture velocities and Reynolds numbers in the experiments of Conan *et al.*^{4}

The height of the clear water layer decreases with increasing input oil volume fraction and is predicted accurately, with absolute deviations from the experiments of Δ*y*^{+} < 0.05. Some differences are seen for *u*_{m,in} = 0.58 m s^{−1} and *ε*_{o,in} = 0.12 [Fig. 2(a)], with the simulations predicting a clear water layer, while the experiments showed that a few drops are still present in the bottom of the pipe. This difference may be attributed to the polydisperity of the mixture, which is not considered in the CFD model. In the upper region of the pipe, some deviations can be seen for most cases. While the shape of the profiles is qualitatively similar for both simulations and experiments, the actual values at each measuring location *y*^{+} can differ by more than the 20% experimental uncertainty, as can be seen in Fig. 9(c).

Finally, a decrease in the oil drop concentration very close to the wall is observed in all cases within that experimental range of conditions. A similar decrease in the profiles was reported by Ekambara *et al.*^{30} for gas-liquid and Ngan^{9} for liquid-liquid dispersed horizontal flows. In both those studies, a lift force was used to capture this trend in numerical models based on an Eulerian-Eulerian approach. When only gravity and shear-induced diffusion were considered in the current simulations, i.e., *F*_{l} = 0 [in Eq. (A8)], the concentration decrease close to the wall was not predicted. The differences shown in Fig. 9 between the model and the experiments close to the top wall can be attributed to the fact that the simulations consider only an average drop size and not the actual distributions found in the experiments [Figs. 6 and 7]. The drops in the distribution with sizes smaller than the average size used in the simulations experience a smaller lift force and increase the concentration close to the wall. The dependence of the lift force on size is reflected in the simulation results since the cases with the largest drop size exhibit the highest concentration gradients close to the wall [Figs. 9(b) and 9(d)].

For homogeneously dispersed liquid-liquid flows, parabolic profiles of the streamwise velocity component have been reported for input dispersed phase volume fractions below *ε*_{d,in} < 0.5,^{10} while for higher concentrations, the velocity profiles resembled power-law shear-thinning ones.^{11,20} In the present work, the input oil volume fraction is below 0.5, but with the settling of the dispersion, it increases locally and can reach values even close to *ε*_{max}, as was revealed by the profiles of Fig. 9. For *o*/$w$ *and* $w$/*o* dispersions, Conan *et al.*^{4} found that the velocities at the dense-packed layer decrease drastically, while the maximum of the velocity profiles is present in the drop-free water layer.

At the beginning of the pipe, where the dispersions are homogeneous, symmetric parabolic velocity profiles are found experimentally and also predicted from the simulations. This behavior, however, changes as the dispersions separate along the pipe. The time-averaged vertical profiles of the streamwise (axial) velocity component of the continuous water phase at *x*^{+} = 135 are plotted together with the simulation results in Fig. 10 for four typical flow conditions. The velocities are normalized with the average mixture velocity *u*_{m,in}. Agreement, close to the mean relative 20% uncertainty considered for the measurements, is found. The profiles are asymmetrical with low velocities at the upper part of the pipe, where the dispersion has become dense and high velocities at the clear water layer at the bottom. The agreement between the experiments, and the simulations illustrates that the use of a mixture viscosity in the model [Eq. (3)] can represent the phenomena taking place in the pipe.

For most cases simulated, the maximum of the velocity profile is underpredicted by the model—even exceeding the deviations of 20% for the more concentrated cases [see Figs. 10(b) and 10(d)]. The vertical location, however, of the maximum velocity is captured very well, with deviations below Δ*y*^{+} ± 0.06. Finally, the case of *u*_{m,in} = 0.58 m s^{−1} and *ε*_{o,in} = 0.12 [Fig. 10(c), corresponding to the PLIF images of Fig. 2(a)], is presented here to illustrate the limitations of the mixture model and the importance of polydispersity. For this case, smaller drops remain at the bottom of the pipe in the experiments. These drops are found to have a dominant effect on the resulting velocity structure, as it is evidenced by Fig. 10(c).

### D. Prediction of flow pattern development

The length of the test section required for the flow pattern to develop is discussed here. While the separation length scale is dictated by buoyancy and is expected to be much shorter than the diffusion length scale, shear-induced migration counteracts the separation due to buoyancy, moving the drops from regions of high to regions of low shear and concentration. It can be seen from the high-speed PLIF images of Fig. 2 that segregation has already taken place at *x*^{+} = 135 for most cases examined.

The development length can be predicted by using the results of the numerical simulations of this work. First, the time-averaged *in situ* oil volume fractions across the pipe diameter *D* are computed as

while an evolution parameter *E*_{p} for the oil volume fraction is introduced to quantify the development length and is based on suspension flows.^{31,32} The evolution parameter is defined as a scalar measure of the profile evolution and is given by

where the overbars and angular brackets characterize space (over *y*) and time (over *t*) average quantities, respectively, whereas ⟨*ε*_{o,0}(*y*)⟩ and $\u27e8\epsilon \xafo,0\u27e9$ are the time-averaged initial (*x* = 0) local and spatially averaged (over *y*) oil volume fractions, respectively. It follows that $\epsilon o,in\u2261\u27e8\epsilon \xafo,0\u27e9$ is constant for each condition and *E*_{p}(*x*) = 0 for a dispersion that remains homogeneous across *y* and along *x*. The simulations are initialized with a homogeneous drop concentration along *y*, so $\u27e8\epsilon o,0(y)\u27e9/\u27e8\epsilon \xafo,0\u27e9=1$.

The evolution parameter computed for the simulations from Eq. (9) for eight axial locations along the pipe is plotted in Fig. 11 against the normalized axial distance for four typical conditions. The evolution parameter increases asymptotically and reaches a plateau. The evolution parameter reaches higher values for the lower *ε*_{o,in} cases, and the plateau is reached at smaller *x*^{+} values for the lower *u*_{m,in} cases. Hampton *et al.*^{31} fitted the growth of *E*_{p} with *x*^{+} to an exponential function (continuous and dashed lines of Fig. 11), which can be written as

where *c*_{1}, *c*_{2}, and *c*_{3} are fitting parameters computed from the least-squares method with *c*_{2} < 0. The coefficient *c*_{3} should be equal to zero to yield the homogeneous profile in the beginning of the pipe but is fitted with small values instead to better describe the results.

From the fitted curves of Fig. 11, an entrance length $L$ can then be defined as the axial location where the evolution parameter reaches 95% of its asymptotic value, written as $Ep\u221e=c1+c3$, and can be further calculated as $0.95Ep\u221ex+=c11\u2212expc2\u2009L0.8+c3$.^{32} The entrance length (normalized with the pipe diameter) then is

with *c*_{4} ≃ 0.8. The entrance lengths $L$ computed from the oil concentration profiles are $L<135$ for the cases considered in the simulations, translating to fully separated flow conditions before the second axial measuring location is reached at *x*^{+} = 135. The effect of the mixture velocity is clear on the development length, with shorter $L$ calculated for lower mixture velocities. However, it is more difficult to quantify the effect of the dispersed phase volume fraction on $L$ as it appears that the drop size distribution plays an intrinsic role to the development.

### E. Measurements of deformation effects

The experimental data can be used to characterize the collisions between the drops and also evaluate the applicability of the mixture model as it assumes that the collisions are between nondeformable particles. The Capillary number can be used to investigate the deformation of the drops. According to Pouplin *et al.*^{10} for homogeneously distributed drops in dispersed pipe flow, deformation can be expected when the Capillary number is higher than 0.1. The maximum Capillary number is written as

For laminar pipe flow, the maximum shear rate is calculated from the velocity gradient at the wall, which can be approximated to $\gamma \u0307max=8um,in/D$. For the conditions that were investigated both experimentally and numerically, the maximum Capillary number of Eq. (12) is close to unity. However, this formulation for the maximum Capillary number [Eq. (12)] is not directly relevant to inhomogeneous dispersions as the drop size and shear rates vary along the vertical direction and may lead to lower local Capillary numbers.

In inhomogeneous dispersions, a local Capillary number can be defined instead, considering the “particle/drop phase pressure” Π, as defined by Jeffrey *et al.*^{33} for suspensions and followed by Abbas *et al.*^{11} for emulsions,

To use Eq. (13), the vertical profiles of the drop size and drop concentration computed from PLIF and the vertical profiles of the velocity gradients computed from PIV are implemented. The resulting local Capillary numbers are presented in Fig. 12 for three typical cases. It is found that the local Capillary numbers are in the majority of cases below the deformation limit (Ca = 0.1). However, for the dense cases considered, Ca(*y*) reaches values close to 1 at the top of the pipe, where the drop concentrations are close to the maximum packing limit, as also evidenced by the profiles of Fig. 9. Deformation of the drops is expected at this small portion of the pipe for the concentrated conditions. For these regions, the shear-induced diffusion theory considered in the mixture modeling approach might not hold true as the drop collision mechanism changes. From Fig. 12, it is also shown that the local Capillary numbers are more sensitive to the input dispersed phase volume fraction *ε*_{d,in} than the mixture velocity.

From the above, it can be concluded that the mixture modeling approach, which assumes nondeformable drops, can be used with deviations expected for local concentrations close to the maximum packing limit. The deviations in the comparisons between the simulations and the experiments of Figs. 9 and 10 can mainly be attributed to drop size effects. The mixture model works with the mean drop size and does not consider the distributions shown in Fig. 6 nor the variations of the mean drop size along the vertical direction, as shown in Fig. 7. Considering polydispersity effects could significantly influence the model predictions, as the size directly affects the gravity and lift forces and the shear-induced diffusion, and consequently the motion of the drops.

## IV. CONCLUSIONS

In this work, experiments were conducted in a horizontal flow test section with two immiscible liquids with matching refractive index. Dispersions at low velocities were generated in the beginning of the test section, and their spatial evolution was investigated by conducting laser-based optical measurements at two axial locations, namely, *x*^{+} = 15 and 135 from the inlet. The PLIF measurements revealed the local flow structures, concentration profiles, and drop sizes, while the PIV measurements recorded the velocity and its gradients in the aqueous phase. Narrower drop size distributions and smaller drops were measured for the same flow conditions for the *o*/$w$ than for the $w$/*o* dispersions. The Sauter mean drop diameter profiles along a vertical diameter were found to be asymmetric with larger drops close to the densely packed layer. In the *o*/$w$ dispersions large drops were observed close to the top part of the pipe but away from it.

In all cases, the dispersions started to segregate and a clear layer of the continuous phase formed. In the $w$/*o* dispersions, a continuous layer of the water phase also formed by the *x*^{+} = 135 axial location. The characteristics of the *o*/$w$ dispersions were simulated with a numerical model based on the mixture approach. In the model, the mixture viscosity was calculated from the model by Krieger and Dougherty.^{24} The model uses a constant drop size. Agreement was found between the experiments and the simulations, with both illustrating asymmetric velocity profiles with the velocity maxima located in the drop-free layer. The model demonstrated that the motion of the drops was influenced by gravity and lift forces together with migration due to shear and concentration gradients. The results also showed that the dispersed drops separate earlier at low mixture velocities compared to high ones. From the concentration, drop size data, and the velocity profiles obtained with the PLIF and the PIV techniques, the assumption of the drop deformability used in the model was considered by calculating the local Capillary numbers.

This work provides detailed and simultaneous measurements of *in situ* phase fractions, velocities, and drop sizes for various flow conditions and flow patterns, which (i) provide an insight into separating liquid-liquid dispersed flows and (ii) are used to validate a CFD model. The model used here can help predict the development and flow characteristics of these liquid-liquid dispersions and can be used as a computationally inexpensive tool for industrial applications. While these initial findings are encouraging, a wider range of flow conditions, covering various flow structures, would be needed to further explore the importance of the assumptions considered here and of the individual mechanisms affecting the motion of the drops (gravity, diffusion, and lift). Specifically, the diffusion constants of Eq. (A9), which are taken from experiments with solid spheres, will need to be validated. In addition, including the drop size distributions and possibly the changes in the drop sizes in the simulations would further help to model developing liquid-liquid dispersed flows more accurately.

## SUPPLEMENTARY MATERIAL

Considering the approximation level of the constitutive equations employed in the model, the majority of which are based on empirical relations and heuristic arguments developed for 2-D systems, the improvement in accuracy by simulating the actual three-dimensional (3-D) geometry is not critical. However, for completeness, the supplementary material section discusses the extension of the model to a 3-D geometry, as well as the results obtained. It is found that both the 2-D and the 3-D simulation results lie close to the experimental uncertainties and any differences between them are not considerable, as shown in Figs. S2 and S3 of the supplementary material.

## ACKNOWLEDGMENTS

This project was funded under the UK Engineering and Physical Sciences Research Council (EPSRC) Programme Grants MEMPHIS (Grant No. EP/K003976) and CORAL (Grant No. EP/N024915/1). VV is also grateful to Chevron Corporation and to University College London for their support. Data supporting this publication can be obtained on request to p.angeli@ucl.ac.uk.

### APPENDIX: NUMERICAL METHODOLOGY

In this paper, the interest in the modeling is placed in the cases in which the secondary phase remains dispersed in the continuous phase and partial separation occurs. It is assumed that the effect of drop coalescence is negligible and that the drops behave like rigid spheres. The validity of this assumption is tested against the experimental findings in Sec. III E.

It is often convenient to use simplified modeling approaches that deal with averaged variables and are able to predict the evolution of the observable parameters (such as volume fractions) directly. One such approach is the mixture model, in which the mixture of the two phases is considered as one fluid, and one continuity and one linear momentum balance equation are solved for it. Hence, one velocity and one pressure field are obtained for the mixture as a whole. To calculate the volume fraction of each phase, the mass balance equation for one of the phases should be solved; this is usually done for the dispersed phase.^{34}

The basic assumption behind the mixture modeling approach is the local dynamical equilibrium between the phases. The dispersed drops should relax in a time interval much shorter than the characteristic time scale of the flow field, i.e., the Stokes number (the ratio of these two time scales) needs to be much smaller than one.^{35,36} Also, to calculate the velocity field of each phase, an additional dynamic equation is required, which is derived from a force balance equation. The result, which calculates the slip velocity between the phases, is usually an algebraic equation. Theoretically, the mixture approach as described before could be used for (nearly) neutrally buoyant drops or particles for which the terminal settling (or rising) velocity is negligible compared to the velocity scale of the mean flow. To justify the use of the mixture approach for the current system, the different time scales that are important in this study have to be investigated.

The first important time scale to consider is the relaxation time, *τ*_{R}, which indicates the time a drop takes to reach local dynamic equilibrium. This time characterizes the drag force, and in the absence of the virtual mass force can be written, according to Jackson,^{35} as

where *β* is the drag coefficient for a dense system of a particulate phase dispersed in a continuous phase. While there are several expressions available for the drag coefficient in the literature, without losing generality of the results, the one reported by Zhang and Acrivos^{12} is selected,

in which *a* is the radius of the dispersed drops. The function *f*_{h}, which is known as the hindrance function, accounts for the reduction in the settling (or rising) velocity due to the presence of other drops around a single drop in a concentrated system. The above hindrance function is selected because it approaches values close to zero at the maximum packing [according to Eq. (3), the mixture viscosity diverges at this limit]. This prevents *ε*_{d} from exceeding the maximum packing limit and generating numerical problems.^{37} In the limit of very dilute systems, this function approaches unity and the drag coefficient coincides with Stokes’ drag law.

The second important time scale is the gravity time, *τ*_{G}, which is the time effective gravity takes to accelerate a particle (in the vertical direction) making it reach a velocity of magnitude equal to that of the local fluid velocity. This time characterizes the effective gravitational field (effective because the buoyancy force reduces the gravitational acceleration to which a particle is subjected). The gravity time is scaled as *τ*_{G} ∼ *U*_{e}/*g*_{e}, where *U*_{e} is the scale of the total local fluid velocity and *g*_{e} = |1 − *ρ*_{c}/*ρ*_{d}|*g* is the effective gravitational acceleration. The first condition required for the mixture modeling approach to hold is *τ*_{R} ≪ *τ*_{G}. This condition ensures that, as the drop relaxes to local dynamic equilibrium, the effective gravitational field does not have enough time to act significantly, which results in a drop terminal velocity of much smaller magnitude than that of the local velocity of the fluid. In other words, during the relaxation time *τ*_{R}, effective gravity results in a terminal velocity *U*_{T} ∼ *τ*_{R}*g*_{e}, which is much smaller than *U*_{e}. In the current work, it is found that *τ*_{R} ∼ 0.001 s and *τ*_{G} ∼ 0.1 s, thus ensuring that *U*_{T} ≪ *U*_{e}.

The third important time scale is the convection time, *τ*_{C}, which indicates the time convection takes to change significantly the local fluid velocity. This time can be scaled as *τ*_{C} ∼ *L*_{C}/*U*_{e}, where *L*_{C} is the length scale required for the local fluid velocity to change significantly along the trajectory of the drops. In the present case, since *U*_{T} ≪ *U*_{e}, it can be assumed that this trajectory is not affected significantly by gravity. Therefore, *L*_{C} can be estimated by the length required for the fluid velocity to change significantly in the horizontal direction; this is approximately the developing length of the velocity profile in the pipe, which is calculated in Sec. III D. The second condition required for the mixture modeling approach to hold is *τ*_{R} ≪ *τ*_{C}. This condition ensures that the drop does relax to local dynamic equilibrium because convection would otherwise change the local fluid velocity as the particle tries to relax and therefore relaxation would be unattainable. In the current case, it is *τ*_{C} ∼ 4 s. It can, therefore, be concluded that both conditions needed for the applicability of the mixture approach are satisfied. It has to be emphasized that this analysis does not mean that the effective gravity has no effect but that to have an appreciable effect one needs to wait a time that is sufficiently long. Specifically, its effect over a time of order *τ*_{R} is negligible, while over a much longer time, its final effect is significant. The mixture residence time in the pipe is about 10 s that is long enough for gravity effects to take place.

The two incompressible phases, namely, the oil and the water/gycerol mixture, are considered in an isothermal laminar flow. Hence, neither energy equation nor equations for turbulence modeling are reported. The continuity equation for the mixture is written as

where

In the aforementioned expressions, ** v** is the averaged velocity and the subscripts

*c*,

*d*, and

*m*stand for the continuous phase, dispersed phase, and mixture, respectively.

The linear momentum balance equation for the mixture is written as

where *p* is the mixture pressure, ** g** is the gravitational vector, and

*ω*

_{c}and

*ω*

_{d}are the mass fractions of the continuous and dispersed phases, respectively. The last term on the right-hand side of this equation is the divergence of the so-called

*diffusion stress tensor*and appears because the velocities of the two phases are not equal.

The mass balance equation for the dispersed phase is required to track the volume fractions of the phases. This equation can be written as

If the convective flux in this equation is written in terms of the mixture velocity, then

Since the dispersed phase does not move at the same velocity as the mixture, an extra term appears in the mass balance equation of this phase, which can be considered as a diffusive flux. Similar to the diffusive flux in the linear momentum balance equation, this term involves the slip velocity between the two phases. The slip velocity is unclosed, but if the two phases reach rapidly local equilibrium, an approximate equation can be obtained for it. This has been done by Jackson^{35} who obtained the following expression for the slip velocity:

in which **Σ**_{d} is the volume-averaged effective stress tensor of the dispersed phase, $Dt$ is the material derivative $Dt$ ≡ *∂*/*∂t* + *v*_{m} · ∇, and *F*_{l} is the lift force exerted by the continuous phase on the dispersed phase per unit volume of the mixture. This equation states that the gravity, inertial acceleration, effective stress tensor of the dispersed phase, and the lift force are responsible for the diffusion flux appearing in Eqs. (A5) and (A7). There are a few methods that are suggested in the literature to consider the migration induced by the effective stress tensor of the dispersed phase. The “diffusive flux model” is used in the present work, which was initially developed for neutrally buoyant particles.^{15} In this category of models, different sources for the migration of the dispersed phase, such as Brownian diffusion as well as shear and concentration gradients, are considered. By neglecting the Brownian diffusion because of the relatively large size of the drops considered here, the divergence of the effective stress tensor of the dispersed phase in this model can be written as

in which *K*_{c} = 0.43 and *K*_{μ} = 0.65 are constants obtained from experimental data for suspensions in shear flows and $\gamma \u0307$ is the local mixture shear rate. This approach has been used to model particle migration in suspensions of solid particles in liquids and can successfully predict the concentration profiles observed in experiments; see, for instance, Refs. 12, 15, and 37–45.

In Eq. (A8), the lift force appears as one of the sources of the dispersed phase diffusion. This force on spherical particles has been investigated extensively in the literature. It is shown that it pushes the particles away from the wall boundaries in their vicinity where the shear rate is high and this lift-induced migration is more pronounced for large particles.^{46} Analytical solutions for the lift-induced velocity of an isolated spherical particle close to a wall have been derived by several authors^{47–51} for small particles and channel (or pipe) Reynolds numbers, defined based on continuous phase properties. Since these models are developed for isolated particles, unlike the drag force, the hindrance effect caused by other particles is not included in them. However, according to Eq. (A8), this effect is accounted for through the drag coefficient *β*. Therefore, an expression which is derived for an isolated particle in a laminar channel flow is used, and then Eq. (A8) is employed to account for the hindrance effect.

Here, the result obtained by Asmolov^{52} is used, who derived an analytical solution for the lift force on an isolated spherical particle close to a wall as a function of the distance from the wall. This force, which holds for a large channel Reynolds number $O{102}\u2212{103}$ in its most general form, reads

where *U*_{e,max} is the maximum velocity of the continuous phase in the channel and *f*_{l} is the force on one particle, which is related to the volume-averaged lift force as *f*_{l} = *F*_{l}*V*_{p}/*ε*_{d}. Here, *V*_{p} = 4*πa*^{3}/3 is the volume of a single particle. The function *k*(*y*/*H*) is the lift coefficient (*y* is the distance of the particle center from the wall, and *H* is the channel width), which is obtained numerically by Asmolov;^{52} at the center of the channel, it is zero and reaches a maximum adjacent to the wall, as it is shown in Fig. 13 for Re_{m} = 300.

The set of Eqs. (A3), (A5), (A7), and (A8) should be solved numerically to predict the hydrodynamics of the current system and provide velocity and concentration profiles that could be compared with the experimental results.

As stated in the Introduction, the purpose of the modeling is not to capture all the dynamics of the system. Instead, the goal is to gain insight into the physical mechanisms leading to the separation of the drops and the thin layer close the pipe walls where the concentration of the drops decreases. For this reason and to reduce the computational costs, a two-dimensional (2-D) geometry with a height (diameter) *D* and length 150*D* was used as the computational domain for the simulations. The COMSOL Multiphysics software was used for the simulation of the two-phase flow in the domain, and a grid independence study was performed using three different numbers of mesh cells. At the end, a structured grid with 1 482 237 degrees of freedom was selected and all the simulations were carried out using the same grid. At the inlet, the mixture velocity and average concentration (defined in Sec. II C) based on experimental values are set equal to the inlet conditions. At the exit, the magnitude of the gauge pressure is set to zero. The walls of the geometry are considered as no-slip.