Characteristics of turbulent core–annular flow with water-lubricated high viscosity oil in a horizontal pipe

Abstract Direct numerical simulations are performed to investigate the characteristics of a turbulent core–annular flow with water-lubricated high viscosity oil in a horizontal pipe. Six different superficial velocity ratios ($\,j_w/j_o = 0.057\unicode{x2013}0.41$) are examined by changing the water superficial velocity $j_w$ for a fixed oil superficial velocity $j_o$. The pressure drops in the pipe and the shapes of the phase interface agree well with those from previous experiments. The oil core flow is almost a plug flow, and the gaps between the phase interface and pipe wall are narrow and wide near the upper and lower surfaces of the pipe, respectively, due to the buoyancy. Within a narrow gap, water is confined mostly in a valley region of the wavy phase interface and hardly goes through its crest. On the other hand, water near the phase interface at a wide gap convects downstream almost at the core speed and the flow near the wall is similar to that of single-phase wall-bounded turbulent flow. The annular flow is characterized by three different regimes depending on the clearance Reynolds number ($Re_c$) based on the core velocity and local gap size: laminar Couette flow driven by the core for $Re_c \le 600$, transitional flow for $600 < Re_c < 2500$ and turbulent flow for $Re_c \ge 2500$. The minimum pressure drop occurs at $j_w/j_o = 0.11$ in the early transition regime. For all $j_w/j_o$ considered, the negative lift force acting on the core comes from the pressure force which balances the buoyancy.


Introduction
Heavy oil reserves occupy approximately 70 % of whole oil reserves (Dusseault 2001), but one of the major challenges with the use of heavy oil is its transport within a pipe † Email address for correspondence: choi@snu.ac.kr because of its high viscosity.Among various drag reduction technologies investigated, water-lubricated transport has been known as an effective tool (Ghosh et al. 2009).This arrangement is called a core-annular flow because high viscosity oil is encapsulated in the pipe core region by less viscous water in the annulus.In experiments, core-annular flows in horizontal and vertical pipes have been considered, and drag reductions by water in the annulus have been observed in wide ranges of oil properties and flow conditions (Charles, Govier & Hodgson 1961;Bai, Chen & Joseph 1992;Prada & Bannwart 2001;Sotgia, Tartarini & Stalio 2008) and also validated in real-scale experiments with several-hundred-metre pipes (Arney et al. 1996;Rodriguez, Bannwart & de Carvalho 2009).The differences in the flow characteristics between the vertical and horizontal pipes are attributed to different gravitational directions.
For a horizontal pipe, the directions of the main flow and gravity are perpendicular to each other and heavy oil rises by the buoyancy because the density of heavy oil is usually lower than that of water (Joseph et al. 1997).The buoyancy force is balanced by the hydrodynamic (negative) lift force, which prevents the core from touching the wall and forms an eccentric core in the pipe.This lift force is caused by the pressure force from water lubrication at low Reynolds number (Ooms et al. 1983) but is dominated by the inertial effect with increasing Reynolds number (Oliemans et al. 1987;Feng, Huang & Joseph 1995;Ooms, Pourquie & Poesio 2012).
Many experimental studies have measured the variations of the pressure drop and oil holdup ( o = V o /V, where V o and V are the oil and whole pipe volumes, respectively) with the superficial velocities of oil and water for core-annular flows in horizontal pipes, and shown that the core-annular flow is effective in reducing the drag on the pipe wall (Charles et al. 1961;Oliemans et al. 1987;Arney et al. 1993;Sotgia et al. 2008;Vuong 2009;Shi, Gourma & Yeung 2017;Tripathi et al. 2017).Vuong (2009; μ o = 0.23-1.07Pa s) and Shi (2015; μ o = 3.3-7.1 Pa s) showed that, when the oil viscosity (μ o ) is high enough, the pressure drop and flow pattern are not significantly changed by μ o at high Reynolds numbers.Arney et al. (1993) observed that a core-annular flow with a non-Newtonian Bingham plastic oil (waxy crude oil/water emulsion: μ o = 200-900 Pa s) produces drag reduction but larger fluctuations in the pressure drop than that of the Newtonian oil (No. 6 fuel oil: μ o = 2.7 Pa s).Based on their own and other experimental data, Arney et al. (1993) suggested a pressure drop model based on the friction factor of single-phase pipe flow and an empirical holdup model, which was later modified by considering the eccentricity and oil fouling effects (Shi et al. 2017).Other pressure drop and oil holdup models have been suggested from several studies (Oliemans, Pots & Trompe 1986;Brauner 1991;Bai et al. 1992;Bannwart 2001;Kim & Choi 2018).
Only a few studies have investigated the characteristics of turbulent core-annular flow in a horizontal pipe by experiments owing to the difficulty in flow measurements.Oliemans et al. (1987) showed that turbulence effects are restricted near the lower surface of the pipe, and Tripathi et al. (2017) observed a broad-banded spectrum of the phase interface wave.They demonstrated that large-wavelength interfacial modes dominate at low Reynolds numbers, while small-wavelength interfacial modes dominate at high Reynolds numbers.However, the dynamics of the phase interface, and the lift and drag forces on the core have been rarely investigated for the turbulent core-annular flow in a horizontal pipe.Thus, the detailed characteristics of turbulent flow and phase interface have been studied mostly by numerical simulations using interface tracking methods.For laminar flow, a few numerical simulations have been conducted to find the pressure distribution across the phase interface between oil and water.For example, Bai, Kelkar & Joseph (1996) observed for axisymmetric laminar core-annular flows that the pressure difference across the phase interface between heavy oil and water increases as the gap size between the interface and pipe wall decreases, and thus suggested that an eccentric core-annular flow can be stably maintained because the gap size below the upper pipe wall is smaller (so, high pressure there) than that above the lower pipe wall (low pressure).This was confirmed by Ooms et al. (2012) for a laminar core-annular flow in a horizontal pipe using numerical simulation.They observed that, for a given wavy phase interface shape, the difference in the water pressures on the upper and lower interfaces increases with increasing eccentricity of the core and balances the buoyancy except at a very low Reynolds number where the direction of the lift force is the opposite to the gravitational direction.
For turbulent flow, studies have been limited to numerical simulations with turbulence models.For example, Shi et al. (2017) showed that the flow statistics change with turbulence models adopted, and a shear stress transport k − ω turbulence model with turbulence damping at the phase interface works better than others.Their numerical results show a 30 % difference in the pressure drop and a 14 % difference in the holdup from those by experiments.Housz et al. (2017) showed using Launder-Sharma k − ω turbulence model that the difference in the pressure drops from experiment and numerical simulation is within 15 %, and their instantaneous amplitude and wavelength of the phase interface wave reasonably agree with those observed in experiments.So far, to the best of our knowledge, the only numerical simulation without using a turbulence model (i.e.direct numerical simulation) for turbulent core-annular flow is that by Kim & Choi (2018) in which the flow in a vertical pipe was considered.Li et al. (2023) conducted Reynolds-Averaged Navier-Stokes simulation using the Launder-Sharma low-Reynolds number k-model at the same condition considered by Kim & Choi (2018).Their simulation observed travelling interfacial waves like those of Kim & Choi (2018), but provided an 18 % lower friction factor and slightly higher holdup ratio.
As summarized above, the understanding of the flow characteristics in turbulent core-annular flow with water-lubricated high viscosity oil in a horizontal pipe is still very limited.Therefore, in the present study, we perform direct numerical simulation of this flow to investigate the asymmetric flow features inside the pipe and the pressure variation across the interface between oil and water.The spatiotemporal deformation of the phase interface is tracked with a level-set method.The details of numerical methods and flow conditions are given in § 2. In § 3, the axial pressure drop is compared with that of the experiment (Sotgia et al. 2008), and the characteristics of water flow in the annulus and high viscosity oil flow in the core are discussed, respectively.The spectral characteristics and dynamics of the phase interface are examined in § § 4.1 and 4.2, respectively, followed by the investigation of near-wall flow dynamics in § 4.3.Finally, a summary is given in § 5.

Computational details
The numerical method to track the phase interface between two immiscible fluids, high viscosity oil and water, is essentially the same as that in our previous study on turbulent core-annular flow in a vertical pipe (Kim & Choi 2018).We use a level-set method to track the interface (Herrmann 2008;Kim & Moin 2011), where x j are the cylindrical coordinates, u j are the corresponding velocity components and φ is the level-set function which is a signed-distance function from the phase interface having positive values in water, negative ones in oil and zero at the phase interface.Equation (2.1) is solved at grids near the phase interface using a third-order total variation diminishing Runge-Kutta method in time (Gottlieb & Shu 1998) and a fifth-order weighted essentially non-oscillatory scheme in space in conjunction with a local Lax-Friedrichs entropy correction (Jiang & Peng 2000).
A partial-differential-equation-based reinitialization method (Sussman, Smereka & Osher 1994;Peng et al. 1999) and a global mass conservation method (Son 2001;Zhang, Zou & Greaves 2010) are used to preserve the signed-distance property, |∇φ| = 1, and compensate the volume loss, respectively.The periodic boundary condition is applied in the axial and azimuthal directions, and the Neumann boundary condition, ∂φ/∂r = 0, is used at the pipe wall.
The governing equations for unsteady incompressible two-phase flow are where −dP/dx is the mean pressure gradient to drive a constant mass flow rate in a pipe, p is the pressure fluctuation, g is the gravitational acceleration, σ is the surface tension coefficient, κ is the curvature, n i is the surface-normal vector on the phase interface, δ ij is the Kronecker delta, i = 1 and 3 for the streamwise (or axial) and vertical (the opposite direction to that of the gravity) directions, respectively, and δ is the delta function (non-zero value at the phase interface and zero otherwise).The density and viscosity are constants for each fluid, but change continuously near the phase interface like where the subscripts w and o denote water and oil, respectively.Here, ψ is the water volume fraction in a control volume calculated from the linearization of the level-set function (van der Pijl et al. 2005;Herrmann 2008).Thus, ψ = 0 or 1 for the cells containing oil or water only, respectively, and 0 < ψ < 1 for the cells containing the phase interface.The density and viscosity of oil used for present simulations are ρ o = 889 kg m −3 and μ o = 0.919 Pa s, and those of water are ρ w = 998 kg m −3 and μ w = 0.001 Pa s, respectively, and the surface tension coefficient is σ = 0.02 N m −1 (Sotgia et al. 2008).We also confirmed that oil used in the experiment was Newtonian (G.Sotgia, private communication).Equations (2.2) and (2.3) are solved using a four-step fractional step method (Choi & Moin 1994), where ρ and μ are the provisional density and viscosity obtained from (2.4), (2.5) and the continuity equation, ∂ρ/∂t + ∂ρu j /∂x j = 0 (see Kim & Choi (2018) for the detail), t is the computational time step size, and the superscript n is the time step index.The second-order central difference scheme is used for all the spatial derivative terms except the convection term near the phase interface where an upwind-type mass-flux scheme is applied (Kim & Moin 2011;Raessi & Pitsch 2012).The Crank-Nicolson method is applied to the convection and diffusion terms to avoid severe time step restriction.To solve the resulting system matrix, a Newton iterative method with an approximate factorization (Choi, Moin & Kim 1993) and an Aitken-type accelerator (Irons & Tuck 1969) are adopted.The surface tension term is explicitly treated with a continuum surface force approach (Brackbill, Kothe & Zemach 1992).An iterative constant-coefficient Poisson solver (Kim & Moin 2011) using a fast Fourier transform is applied to solve the variable-coefficient pressure Poisson equation (2.8).The details of numerical methods are available in Kim & Choi (2018).
Figure 1 shows the schematic diagram of a core-annular flow in a horizontal pipe.The gravity direction is perpendicular to the mean flow direction, and thus the oil core naturally rises by the buoyancy due to its lower density than that of water.The cylindrical coordinates (x, r, θ) and the corresponding velocity components (u, v, w) are used for simulation, and the y and z coordinates are used for convenience to represent the lateral and vertical directions at the pipe cross-section, respectively.We fix the superficial velocity of oil, j o = q o /πR 2 = 0.88 m s −1 , and vary the superficial velocity of water as j w = q w /πR 2 = 0.05, 0.075, 0.1, 0.15, 0.22 and 0.36 m s −1 , where R(= 0.013 m) is the pipe radius, and q o and q w are the volume flow rates of oil and water, respectively.This is one of the experimental conditions in Sotgia et al. (2008), where a core-annular flow is maintained in the pipe.Note that Sotgia et al. (2008) varied the oil and water superficial velocities together with the pipe radius and showed various flow regimes including core-annular and dispersed flows.The numerical method of maintaining constant superficial velocities of both fluids is given in Appendix A.
Table 1 shows the bulk Reynolds number, Re b = u b (2R)/ν w = ( j w + j o )(2R)/ν w , number of grid points, computed friction Reynolds number, Re τ = u τ R/ν w , and oil holdup, o = V o /V, for six superficial velocity ratios, j w /j o .Here, V is the total computational volume, V o is the volume occupied by oil, u τ = √ τw /ρ w is the friction velocity, τw is the mean wall shear stress and ν w is the kinematic viscosity of water.For comparison, we also conduct a numerical simulation of water flow at a similar Reynolds number (Re b = 24 580), whose results agree well with those of Wu & Moin (2008).The streamwise domain size (L x ) is 2πR except that of single-phase flow simulation (10R), and periodic boundary conditions are applied in the streamwise (x) and azimuthal (θ) directions where uniform grids are used.As shown in § 3, large-scale structures in the core-annular flow are confined by the wavelength of the phase interface, and thus the streamwise domain size of L x = 2πR is large enough to contain those structures for the cases considered.The no-slip boundary condition is used at the wall, and in the radial (r) direction dense grids are allocated near the wall and interface, but coarse grids are used in the core region where flow is laminar due to the high viscosity of oil.Grid resolutions in wall units are x + = xu τ /ν w ≈ 6-7, r + min ≈ 0.6-0.7 and (R θ) + ≈ 6-7.The number of grid points used for the level-set equation is twice that for the Navier-Stokes equation in each direction.We simulated flow with one and half times the grid points for the streamwise and azimuthal directions for the case of j w /j o = 0.17, and the result showed the changes in the pressure drop and oil holdup less than 4 % and 1 %, respectively.All the computations are carried out for the non-dimensional time of approximately 250R/j o , and averaging is conducted over 150R/j o to obtain mean values.

Characteristics of the core and annular flows
Figure 2 shows the variation of the mean pressure gradient with the superficial velocity ratio for a fixed superficial oil velocity ( j o = q o /πR 2 = 0.88 m s −1 ), together with the experimental result by Sotgia et al. (2008).An excellent agreement between the present and experimental data is observed for j w /j o ≥ 0.17, but the agreement becomes poorer at lower j w /j o .This relative disagreement at low j w /j o may be attributed to the fact that oil frequently sticks to the wall due to a narrow water region (Bai et al. 1992;Arney et al. 1993), which may change the relative location of the phase interface and thus the mean pressure drop (note that the present numerical method does not prevent the formation of oil drops or oil sticking to the wall, but such phenomena do not occur during the present simulations).A similar discrepancy at low j w /j o was also observed for turbulent core-annular flow in a vertical pipe (Kim & Choi 2018).With increasing j w /j o , the mean pressure gradient first decreases, reaches a minimum and then increases, providing an optimal superficial velocity of water ( j w | opt ) for the lowest mean pressure gradient to transport a given oil flow rate (q o or j o ) in a pipe.This has also been observed in other experiments of core-annular flows in vertical and horizontal pipes (Bai et al. 1992;Prada & Bannwart 2001;Sotgia et al. 2008;Vuong 2009;Housz et al. 2017).The optimal superficial velocity ratios obtained from present numerical simulation and experiment by Sotgia et al. (2008) are slightly different ( j w /j o | opt = 0.11 and 0.085, respectively), because flow conditions are not exactly same due to non-hydrodynamic effects such as chemical adhesion (Arney et al. 1996) as described before.When j w /j o is smaller than the optimal one, the mean pressure gradient increases because of high wall shear stress from the narrow gap between the phase interface and upper wall (for j w /j o = 0.057, the minimum gap size is 8 wall units, where 9 and 18 grid points for the Navier-Stokes and level-set equations are located).On the other hand, when j w /j o is larger than the optimal one, the pressure gradient increases again because of the increase in the bulk Reynolds number, Re b = ( j w + j o )(2R)/ν w .For large superficial velocity ratios, the corresponding water flows are similar to single-phase flows (see later in figure 6c), and thus the mean pressure gradient follows the Blasius friction factor formula, , where u b = j w + j o .Since Re −1/4 b little changes in the range of the Reynolds numbers considered (see table 1), the mean pressure gradient normalized by j 2 o increases with increasing j w /j o for j w /j o ≥ 0.17 (figure 2), i.e.
The mean pressure gradient of the core-annular flow normalized by ρ w and u b (= Figure 3 shows the shapes of the instantaneous phase interface for different superficial velocity ratios, together with snapshots of the phase interface from experiment by Sotgia et al. (2008).Also shown in figure 3(af ) are the contours of the instantaneous relative velocity to the core velocity in water flow, where interface qualitatively agree well with the experimental ones considering the differences in the experimental and simulation set-ups as mentioned above.Due to the buoyancy, the flow characteristics vary in the azimuthal direction.The lower gap between the phase interface and pipe wall rapidly increases with increasing superficial velocity ratio, whereas the upper gap slowly increases (see below).The phase interfaces consist of different streamwise and azimuthal wavenumber components, and the wavelength and wave amplitude increase with increasing superficial velocity ratio (or increasing gap size).When the gap is narrow, the crest of the phase interface almost touches the wall, and no small-scale vortices are observed there (see later in figure 5a).With increasing j w /j o , the flow in the gap changes from laminar to turbulent, especially near the lower pipe wall.
where h is the instantaneous radial location of the phase interface, the superscript .denotes the averaging over the streamwise direction and in time, and T is the integration time.The averaged value of h(θ ) over the azimuthal direction is , and they are h avg /R = 0.96, 0.94, 0.93, 0.90, 0.87 and 0.81 for j w /j o = 0.057, 0.085, 0.11, 0.17, 0.25 and 0.41, respectively.For j w /j o = 0.057, h(θ ) is nearly constant in the azimuthal direction due to the small amount of water flow, but, for higher j w /j o , it has a plateau near the top wall (θ ≈ 0 • ) but rapidly decreases near the bottom wall (θ ≈ 180 • ) because the core rises upward by the buoyancy.The gap size between the phase interface and wall, R − h(θ ), increases with j w /j o for all azimuthal locations (figure 4b).At large superficial velocity ratios, the gap size increases almost linearly with the superficial velocity ratio.This linear growth in the gap size starts to occur at lower superficial velocity ratio for larger θ .
Figure 5 shows the contours of the instantaneous modified pressure fluctuations, p * (= p + ρ w gz − p * wall ), and the instantaneous velocity vectors relative to the core velocity for j w /j o = 0.057 and 0.41, where p * wall = ( p + ρ w gz) dA wall /A wall and A wall is the area of the pipe wall.Here, we add ρ w gz in the modified pressure to remove the hydrostatic effect in the pressure.For the narrow gap of j w /j o = 0.057, there is a recirculating flow in the wave valley, and high-and low-pressure fluctuations are observed ahead of and behind the crest, which is similar to those observed from core-annular flows in vertical pipes (Bai et al. 1996;Li & Renardy 1999;Kim & Choi 2018).However, for j w /j o = 0.41, the gap between the interface and lower pipe wall is sufficiently wide and turbulent flow is observed within an oil-free region near the lower wall.Near the crest, p * is high and low ahead of and behind the crest, respectively.This indicates that the oil flow in the core drags the water flow.Also, p * is overall high and low near the upper and lower surfaces of the pipe, respectively, and this difference is in balance with the buoyancy (Feng et al. 1995;Ooms et al. 2012), which is discussed later in § 4.2. Figure 6 shows the mean streamwise velocity profiles of water (ũ w ) and oil (ũ o ) normalized by the core velocity u core and local friction velocity ũτ (θ ), respectively, for j w /j o = 0.057, 0.17 and 0.41, where ũw , ũo and ũτ are obtained as Since the core flow is almost a plug flow, ũo /u core ≈ 1 but slightly decreases with increasing r.The annular flow may be considered as a Poiseuille-Couette-type flow driven by both the mean pressure gradient and the core.For j w /j o = 0.057 (figure 6a), the gap between the interface and wall is narrow and little varies along the azimuthal direction (see and core (u core ) and annular (u annular ) velocities with j w /j o .The holdup ratio is the ratio of the bulk velocity of oil (core) to that of water (annular) and is greater than 1, because the annular water flow is heavily influenced by the viscous effect.With increasing j w /j o , the holdup ratio decreases from 1.50 to 1.29, indicating that the amount of increase in V w /V o is smaller than that in j w /j o (= q w /q o ).With increasing j w /j o , the oil holdup o decreases, and the core and annular velocities normalized by j o (u core /j o = 1/ o ) increase.
Figure 7(b) shows the contours of the mean water volume fraction ψ(r, θ) and the mean location of the phase interface h(θ ) for j w /j o = 0.41, where ψ(r, θ) = T 0 L x 0 ψ(x, r, θ, t) dx dt/(L x T).The region of 0 < ψ < 1 is wider near the lower pipe wall than that near the upper one, which indicates that the amplitude of the phase interface wave is large and small near the lower and upper pipe walls, respectively.As shown with h(θ ), the core is deformed due to the buoyancy and is eccentric to the pipe centre.Figure 7(c) shows the roundness χ and eccentricity e of the phase interface.The roundness χ is defined as (the ratio of the volume of the core to that of a circular cylinder having the

A19-12
Turbulent core-annular flow in a horizontal pipe same surface area of the core) and the eccentricity is e = z cm /(R − h avg ), where z cm (= is the vertical distance from the pipe centre to the centre of mass of the core, and h avg = 2π 0 h(θ ) dθ/2π is the overall mean radius of the phase interface.To calculate the perimeter h dθ , the trapezoidal rule is used.The roundness of the circle is 1, but that of the core is less than 1 and decreases with increasing j w /j o .For j w /j o = 0.41, the interface is most deformed among the cases considered and the roundness is approximately 0.94, which is similar to that of the ellipse whose major and minor axes are 1.5 and 1, respectively.Instantaneous interface shapes are much more non-circular than the mean shape shown in figure 7(b).The eccentricity is very small for j w /j o = 0.057, indicating that the phase interface is nearly circular because oil occupies almost the whole circular pipe.The eccentricity rapidly increases from j w /j o = 0.057 to 0.085, is nearly constant for 0.085 < j w /j o < 0.25 and then slowly increases with increasing j w /j o , because the core rises upward.

Wave characteristics of the phase interface
The wave characteristics of the phase interface are investigated by computing the streamwise wavenumber (k x ) and frequency (ω) power spectra of the phase interface amplitude, ζ(θ)(= h − h(θ )).The range of the streamwise wavenumber is 0 ≤ k x ≤ 383/R ( k x = 1/R), and the sampling interval is t s = 0.02R/j o during T = 81.92R/jo which is divided into 15 overlapping segments with 50 % overlap (0 ≤ ω ≤ 157j o /R with ω = 0.6136j o /R; for more details, see Choi & Moin (1990)).The one-dimensional wavenumber spectrum ϕ(k x , θ) and frequency spectrum ϕ(ω, θ) of the phase interface amplitude ζ satisfy the following condition: where ζ rms (θ ) is the root-mean-square fluctuations of the phase interface amplitude.
Figure 8 shows the streamwise wavenumber and frequency spectra of the phase interface amplitude at θ = 0 • and 180 • for j w /j o = 0.057, 0.17 and 0.41.The spectra are broad-banded, indicating that the motion of the phase interface amplitude has a turbulent nature.For j w /j o = 0.057, the streamwise wavenumber spectrum is more or less homogeneous in the azimuthal direction, and high powers appear in low wavenumbers, k x R ≤ 12 (λ ≥ πR/6), where λ is the corresponding wavelength.With increasing j w /j o , the power peak moves to lower wavenumbers on the bottom interface (θ = 180 • ), whereas the power on the top interface is relatively insensitive to j w /j o .This is consistent with the variation of large-scale wavy structures on the bottom interface with j w /j o in figure 3: the power peak occurs at k x R = 12, 5 and 3 for j w /j o = 0.057, 0.17 and 0.41, respectively, corresponding to the dominant wavelengths of λ/R = 0.52, 1.26 and 2.09.Note also that the power at low wavenumbers and frequencies is much larger at θ = 180 • than at θ = 0 • for j w /j o = 0.17 and 0.41.The frequency spectra are very similar to the streamwise wavenumber spectra, showing the convective nature of the phase interface.
Figure 9(a) shows the contours of the streamwise wavenumber-frequency power spectrum of the phase interface amplitude, Φ(k x , ω) = 2π 0 Φ(k x , ω, θ) dθ/2π, for j w /j o = 0.17, where Φ(k x , ω, θ) is the two-dimensional power spectrum as a function of the azimuthal angle.As shown, the streamwise wavenumber-frequency power spectrum shows a strong convective nature.The convection velocity of the phase interface can be obtained as (Wills 1971;Choi & Moin 1990) where ω(k x , θ) is obtained from A quadratic polynomial is used to find ω(k x , θ) in (4.3) (Kang, Moin & Iaccarino 2008).Figure 9(b) shows the convection velocity normalized by the core velocity as a function of the streamwise wavenumbers at three different azimuthal angles.The scatters of the data shown in this figure are due to the limited statistical sample available for locating the maxima of the spectrum within a set of discrete frequencies and wavenumbers, and the convection velocities at low wavenumbers are also contaminated by a specific window function (Hanning window) or the sampling period used (see also Choi & Moin (1990) and Kim & Choi (2018)).The convection velocity is smaller than the core velocity, indicating that the core drags the phase interface.At low wavenumbers where the spectrum has a peak (k x R = 1-10), the convection velocities are lower than those at higher wavenumbers ( Ũconv (k x , θ)/u core ≈ 0.94-0.98),suggesting that the energy-containing large-scale structures of the phase interface more strongly interact with the water flow and resist the core flow more than the small-scale structures do.This characteristic is similar to that of the core-annular flow in a vertical pipe (Kim & Choi 2018).Note also that the convection velocity at θ = 0 • is lower than those at θ = 90 • and 180 • due to the viscous effect at the narrow gap, although the difference is not so large.The overall convection velocity U conv of the most energetic wave is obtained from where The overall convection velocity increases with increasing j w /j o ; i.e.U conv /u core ≈ 0.94, 0.95 and 0.96 for j w /j o = 0.057, 0.17 and 0.41, respectively.This indicates that most energetic motions move downstream at a speed slightly lower than the core velocity.

Dynamics of oil core
The stress on the phase interface is obtained as where the subscripts i = 1, 2 and 3 indicate the axial (x), lateral (y) and vertical (z) coordinates, respectively (figure 1).The pressure and velocity gradient on the interface are obtained by a linear interpolation from grid points at the nearest oil region to those at the phase interface.From the stress boundary condition at the phase interface ( jump condition), the stress at the phase interface from the core region is the sum of the stress at the phase interface from the annular region and the surface tension (Brackbill et al. 1992).
Figure 10 shows the contours of the instantaneous stresses on the interface, σ zj n j /(ρ w − ρ o )gR and σ xj n j /ρ w u 2 τ , and their components for j w /j o = 0.41.Here, the viscous normal stress is included in the pressure.Note that the stresses σ zj and σ xj are normalized by the hydrostatic pressure difference and wall shear stress, respectively, to explain σ zj and σ xj in terms of the lift and drag forces, respectively.The stagnation pressure ahead of the crest generates the wall-repulsive (−r direction) and drag (−x direction) forces (in this figure, δ xj n j are positive and negative on the forward and leeward sides, respectively).The pressure on the upper interface is higher than that of the lower one (δ zj n j are positive and negative on the upper and lower interfaces, respectively), providing a negative lift (−z  direction) force on the core which balances its buoyancy force, whereas the contribution of the viscous shear stress to the lift force is very small (figure 10a).For the drag force (figure 10b), both the pressure and viscous shear stress are highly positive and negative on the forward side of the crest, respectively, and thus push the crest in the upstream direction, whereas low pressure and positive viscous shear stress are generated on the leeward side of the crest, indicating that the flow separates across the crest.
Figure 11 and A core is the surface area of the phase interface at the same oil holdup.The analytical mean drag coefficient of C D = √ o χ is obtained by combining the momentum balances on the whole domain and oil core, respectively, and introducing χ = 4πA o /P 2 , o = V o /V = A o /πR 2 and A core = PL x , where A o is the cross-sectional area of the oil core, and P is the wetted perimeter of the phase interface (see Appendix B for the derivation).As shown in this figure, the computations of C L and C D on the phase interface using (4.7) are not very accurate but contain maximum 10 % errors because the surface tension force is smeared out from the phase interface by using the continuum surface force approach (Brackbill et al. 1992).As discussed in figure 10(a), the contribution of the pressure to the lift is dominant and that of the viscous shear is very small (Oliemans et al. 1987;Feng et al. 1995).For the drag coefficient, both the pressure and viscous shear stress on the interface are important.For low j w /j o , the contribution of the viscous shear stress is bigger than that of the pressure owing to the narrow gap between the interface and pipe wall.With increasing j w /j o , the gap size increases and thus the contribution of the viscous shear stress on the interface to the drag decreases, while that of the pressure is nearly unchanged, resulting in the decrease in the drag coefficient.

Flow transition and near-wall dynamics
Figure 12 shows the azimuthal variation of the clearance Reynolds number Re c (θ ) and the variation of the mean wall shear stress τw (θ ) with Re c (θ ) for various j w /j o , where Re c (θ ) = u core (R − havg (θ ))/ν w .The clearance Reynolds numbers for j w /j o = 0.085, 0.11 and 0.17 are mostly in the transitional region.The normalized mean wall shear stresses, τw (θ )/ρ w u 2 core , follow the relations of the laminar flow at low Re c (narrow gap) and of turbulent flow at high Re c (wide gap), respectively, with a transient behaviour in between them.At low Re c , the velocity profile of water flow is almost linear from the wall to the core (see figure 6a).Assuming a linear velocity profile having u core at r = havg (θ ), the wall shear stress at low Re c can be approximated as Although the Blasius friction factor formula applies to a turbulent Poiseuille flow, Orlandi, Bernardini & Pirozzoli (2015) showed that the friction factor of a turbulent Couette flow in a channel is also proportional to Re .The present result suggests that the friction factor (normalized by u core ) also decreases with Re −1/4 c .For 600 < Re c < 2500, τw (θ ) is in the transitional region.It is worth reporting that this Reynolds number dependency of the wall shear stress justifies the use of the pressure drop model by Arney et al. (1993), where their model was based on laminar and turbulent friction factor formulae for single-phase pipe flow with a different definition of the Reynolds number.
Figure 12(a) indicates that flow in water is laminar ( j w /j o = 0.057), laminar at 0 • ≤ θ < 100 • but transitional at 100 • θ ≤ 180 • ( j w /j o = 0.085), early transitional ( j w /j o = 0.11), fully transitional ( j w /j o = 0.17) and transitional and turbulent ( j w /j o = 0.25 and 0.41).The drag on the pipe wall comes from the integration from θ = 0 • to 180 • .The case of j w /j o = 0.11 contains only early transitional flow where the skin friction does not reach that of fully transitional flow, resulting in lowest drag (figure 2).Note that some laminar flow region contains higher skin friction than that of early transitional flow (figure 12b), and thus the drag at j w /j o = 0.057 is larger than those at j w /j o = 0.085 and 0.11.
Figure 13 shows the two-dimensional power spectra of the wall shear stress, Φ(k x , ω, θ), at three different azimuthal angles for j w /j o = 0.17, together with that of single-phase turbulent pipe flow.At the top of the pipe wall (Re c ≈ 750 at θ = 0 • ; figure 13a), the power spectrum, which is similar to that of the phase interface (figure 9a), indicates that the wall shear stress is convection dominated with a nearly constant convection velocity for all streamwise wavenumbers.Since the laminar Couette-type flow is dominant at this azimuthal location, the convection velocity obtained from (4.2) and (4.3) is slightly smaller than u core (Bai 1995;Rodriguez & Bannwart 2006).On the other hand, at the bottom of the pipe (Re c ≈ 2500 at θ = 180 • ; figure 13c), the power spectrum is broad-banded like that of single-phase flow (figure 13d) but has a steeper slope than that of single-phase flow, indicating that the flow in the bottom gap is turbulent but the convection velocity of the wall shear stress is higher than that of single-phase flow.At the side of the pipe (Re c ≈ 1260 at θ = 90 • ; figure 13b), the power spectrum contains both characteristics shown in figure 13(a,c).The two-dimensional power spectra of the wall shear stress for j w /j o = 0.057 and 0.41 are shown in figure 14.For j w /j o = 0.057, the clearance Reynolds numbers are quite small (Re c ≤ 560), and thus the spectra are very similar to that in figure 13(a).On the other hand, for j w /j o = 0.41, Re c ≥ 1400 and thus the spectra are more like those in figure 13(b,c).Note that the power spectrum at θ = 0 • is no longer like that of laminar Couette flow but like the transitional flow due to the wide gap there.

Summary
We numerically investigated the characteristics of the core-annular flow of high viscosity oil and water in a horizontal pipe for six different superficial velocity ratios of water-to-oil.The mean pressure gradient in a pipe and instantaneous shape of the phase interface agreed well with those of the experiment at the same flow conditions.The flow and wave characteristics of the core-annular flow in a horizontal pipe were different from those in a vertical pipe (Kim & Choi 2018) because of the buoyancy.Thus, the core rose up due to the buoyancy and was eccentric to the pipe centre, and the flow characteristics changed along the azimuthal direction.
The flow inside the core was nearly uniform and laminar, whereas the flow in the gap changed significantly with the gap size.Especially, the flow in the gap showed two different characteristics in the annular flow depending on the clearance Reynolds number Re c based on the core velocity and local gap size.For low Re c , small-scale vortices were rarely observed and the annular flow was similar to the laminar Couette flow driven by the core.On the other hand, for high Re c , there existed an oil-free region near the wall and the flow was very similar to that of single-phase turbulent pipe flow.With increasing superficial velocity ratio j w /j o , the gap size increased, and its rate of increase was faster at the bottom wall than that at the top wall.Thus, the transition to turbulence started from the bottom wall to the top one.The wall shear stress was proportional to Re  (600 < Re c < 2500), transition to turbulence occurred and a local minimum of the wall shear stress existed.Since Re c for j w /j o = 0.11 were in this early transition region, the minimum mean pressure gradient occurred at this j w /j o for a given j o (figure 2).The lift and drag forces acting on the core were analysed.The stagnation pressure in the forward side of the crest generated the wall-repulsive force which was higher near the top wall than near the bottom one, which balanced the buoyancy.The drag coefficient from the viscous shear stress was high at low j w /j o due to high wall shear stress in the narrow gap and decreased with increasing j w /j o , but that from the pressure was more or less similar for j w /j o considered.
where ṁn+1 o (= ρ o u n+1 b,o V n+1 o /L x ) and ṁn+1 w (= ρ w u n+1 b,w V n+1 w /L x ) are the mass flow rates of oil and water at the(n + 1)th time step, and u b,o and u b,w are the oil and water bulk velocities, respectively.At t = 0, we provide the analytic velocity profile of a laminar core-annular flow having a cylindrically shaped phase interface (Bai et al. 1992), and adjust the maximum streamwise velocity and radial location of the phase interface to match the given mass flow rates of oil and water.
With the mean pressure gradient obtained from (A3), the total mass flow rate is maintained during simulation; i.e. ṁn+1 To relocate the phase interface to satisfy the volume ratio in (A4), we use an algorithm for the global mass conservation for the level-set method (Son 2001;Zhang et al. 2010), where τ v is the pseudotime for the volume correction iteration.Since the stress and the pressure change smoothly near the phase interface in the present study, a small change in the location of the phase interface has a negligible effect on the real dynamics, and the flow can reach a fully developed state by maintaining the mass flow rates of both fluids.
Figure 15 shows the time histories of the superficial velocities of oil and water and the oil volume for j w /j o = 0.11 using the present numerical method.The oil and water superficial velocities and oil volume change in time, but their magnitudes of the fluctuations are negligible.

Appendix B. Analytic mean drag coefficient on the oil core
Consider a two-dimensional core annular flow where oil is encapsulated in the pipe core by water in the annulus (figure 16).The momentum balances on oil and water and oil only in the axial direction, respectively, provide   drag coefficient on the oil core as where A core and P are the surface area and wetted perimeter of the phase interface, respectively.Note that this analytical solution is valid when the phase interface does not vary along the axial direction.

Figure 1 .
Figure 1.Schematic diagram of the core-annular flow in a horizontal pipe.
Re 2 b (from momentum analysis) ≈ 0.0067 at j w /j o = 0.085 (Re b = 24 780), which is similar to 0.0064 obtained for the single-phase water flow at Re b = 24 580, and is much lower than (−dP/dx)(R/ρ w u 2 b ) = (16/Re b )(ρ o /ρ w ) ≈ 0.64 for single-phase high-viscosity laminar oil flow (u b = j o = 0.88 m s −1 ; Re b = 22.13; normalized by ρ w for comparison).This clearly indicates that the core-annular flow significantly reduces the drag for the heavy-oil delivery.

Figure 4 .
Figure 4. Mean radial locations of the phase interface and mean gap sizes along the azimuthal direction for different superficial velocity ratios: (a) mean radial location, h(θ); (b) mean gap size, R − h(θ).Here, θ = 0 • and 180 • denote the top and bottom locations of the pipe, respectively.

Figure 4
Figure 4(a) shows the mean radial locations of the phase interface in the azimuthal direction for different superficial velocity ratios,

Figure 5 .
Figure 5. Contours of the instantaneous modified pressure fluctuations and relative velocity vectors: (a) j w /j o = 0.057; (b) 0.41.The velocity vectors relative to the core velocity are drawn at every eighth grid points in the streamwise and radial directions, except for the zoomed view in (a) where they are drawn at every other grid point.

Figure 7 .
Figure 7. Mean characteristics of the core flow: (a) holdup ratio (s), oil holdup ( o ), core velocity (u core /j o ) and annular velocity (u annular /j o ); (b) contours of the mean water volume fraction ψ (black solid line) and the mean location of the phase interface h(θ) (red solid line) for j w /j o = 0.41; (c) eccentricity (e) and roundness (χ).Contour levels of ψ are from 0.05 to 0.95 by increments of 0.1.

Figure 7
Figure 7(a) shows the variations of the mean core-flow characteristics including the holdup ratio(s = (q o /q w )/(V o /V w )), oil holdup ( o = V o /(V w + V o ))and core (u core ) and annular (u annular ) velocities with j w /j o .The holdup ratio is the ratio of the bulk velocity of oil (core) to that of water (annular) and is greater than 1, because the annular water flow is heavily influenced by the viscous effect.With increasing j w /j o , the holdup ratio decreases from 1.50 to 1.29, indicating that the amount of increase in V w /V o is smaller than that in j w /j o (= q w /q o ).With increasing j w /j o , the oil holdup o decreases, and the core and annular velocities normalized by j o (u core /j o = 1/ o ) increase.Figure7(b) shows the contours of the mean water volume fraction ψ(r, θ) and the mean location of the phase interface h(θ ) for j w /j o = 0.41, where ψ(r, θ) =

Figure 11 .
Figure 11.Variations of the mean lift and drag force coefficients on the core with j w /j o : (a) lift coefficient; (b) drag coefficient.Here, , viscous stress; +, pressure; •, their sum.In (b), solid circles (•) indicate the analytic mean drag coefficient on the phase interface.

Figure 12 .
Figure 12.Clearance Reynolds number and mean wall shear stress: (a) azimuthal variation of the clearance Reynolds number; (b) mean shear stress versus clearance Reynolds number.The region from laminar to turbulent transition is grey coloured in (a).The dashed lines in (b) are the relations of the wall shear stress with the clearance Reynolds number for laminar and turbulent flows, respectively.
.8) As shown in figure 12(b), (4.8) holds for Re c ≤ 600.For j w /j o = 0.057, all Re c are within this regime.On the other hand, for Re c ≥ 2500, τw (θ ) follows the following relation (similar to the Blasius friction factor formula):

for
Re c ≤ 600 and Re c ≥ 2500, respectively.In between these two regions 986
ΣF o+w = −D + ( p o + p o )A o + ( p w + p w )A w − p o A o − p w A w = 0, (B1) ΣF o = −D i + ( p o + p o )A o − p o A o = 0, (B2)where D and D i are the forces on the pipe wall and phase interface in the axial direction, p o and p w are the oil and water mean pressures, A o and A w are the oil and water cross-sectional areas, respectively, and L x is the pipe length.Here, p o = p w = p, and A o + A w = πR 2 .Then, from (B1) and (B2),D i = DA o /πR 2 .Introducing D = τw 2πRL x , τw = ρ w u 2 τ , o = V o /V = A o /(πR 2), A core = PL x and χ = 4πA o /P 2 provide the mean

Figure 15 .
Figure 15.Time histories of the superficial velocities and oil volume ( j w /j o = 0.11): (a) superficial velocities of oil (black solid line) and water (black dashed line); (b) oil volume.Here, Vo is the time average of oil volume at the fully developed state.

Figure 16 .
Figure 16.Schematic diagram of momentum balances on oil and water and oil only.

Table 1 .
Superficial velocity ratio, bulk Reynolds number, number of grid points, computed friction Reynolds number and oil holdup.