Phase-field simulation of core-annular pipe flow

TL;DR

This paper evaluates the accuracy of phase-field methods in simulating core-annular pipe flows, demonstrating their effectiveness in handling complex geometries, high viscosity contrasts, and unsteady conditions.

Contribution

It provides a quantitative analysis of phase-field simulation convergence to sharp-interface limits in realistic pipe flow scenarios, validating their applicability.

Findings

Phase-field methods accurately simulate non-rectangular geometries.

They handle strong advection and unsteady fluctuations effectively.

The methods work well with large viscosity contrasts.

Abstract

Phase-field methods have long been used to model the flow of immiscible fluids. Their ability to naturally capture interface topological changes is widely recognized, but their accuracy in simulating flows of real fluids in practical geometries is not established. We here quantitatively investigate the convergence of the phase-field method to the sharp-interface limit with simulations of two-phase pipe flow. We focus on core-annular flows, in which a highly viscous fluid is lubricated by a less viscous fluid, and validate our simulations with an analytic laminar solution, a formal linear stability analysis and also in the fully nonlinear regime. We demonstrate the ability of the phase-field method to accurately deal with non-rectangular geometry, strong advection, unsteady fluctuations and large viscosity contrast. We argue that phase-field methods are very promising for quantitatively studying moderately turbulent flows, especially at high concentrations of the disperse phase.

Tables

Table 1. Table 1 : Parameters of the two cases considered for the linear stability analysis, following Hu & Joseph  [ 11 ] . In both cases, the most unstable disturbance is axisymmetric ( m = 0 𝑚 0 m=0 ), evolving into an interfacial-friction (CAF1) or inertial (CAF2) instability.

Test case	$Re$	$Ca$	$\hat{\mu }$	$\hat{\rho }$	$\eta$	$\hat{V}$
CAF1	37.78	$\infty$	1/2	1	0.7	0.9608
CAF2	500	0.5	1/20	1	0.9	4.2632

Table 2. Table 2 : Dimensionless parameters of CAF in the nonlinear regime (bamboo wave) studied experimentally by Bai et al .  [ 3 ] and numerically by Li and Renardy  [ 17 ] . The pipe radius was R = 0.476 𝑅 0.476 R=0.476 cm, whereas ρ oil subscript 𝜌 oil \rho_{\text{oil}} =0.905 g/cm 3 , ρ water subscript 𝜌 water \rho_{\text{water}} =0.995 g/cm 3 , μ oil = 6.01 subscript 𝜇 oil 6.01 \mu_{\text{oil}}=6.01 dyn ⋅ ⋅ \cdot s/cm 2 , μ water = 0.01 subscript 𝜇 water 0.01 \mu_{\text{water}}=0.01 dyn ⋅ ⋅ \cdot s/cm 2 and σ = 8.54 𝜎 8.54 \sigma=8.54 dyn/cm. From these experimental parameters, Li & Renardy [ 17 ] inferred the position of the interface of the core oil flow ( η = 0.78125 𝜂 0.78125 \eta=0.78125 ) and the driving pressure gradient ( f = 1.7455 𝑓 1.7455 f=1.7455 ) for their simulations.

Test case	$Re$	$Ca$	$\hat{\mu }$	$\hat{\rho }$	$\hat{V}$	$Fr$	$f$	$L$
CAF3	1.19	11.7	1/601	1.1	1.566	0.59	1.7455	2.45

Equations

$$\rho (C)=\frac{{\rho }_1+{\rho }_2}{2}+({\rho }_2-{\rho }_1)C,\mu (C)=\frac{{\mu }_1+{\mu }_2}{2}+({\mu }_2-{\mu }_1)C,$$

$$\mathcal{F}={\int }_{\text{V}}\frac{\alpha }{2}{\left(\nabla C\right)}^2+\beta \Psi (C)d\text{V},$$

$$\Phi =\frac{\delta \mathcal{F}}{\delta C}=\beta {\Psi }^{\prime }(C)-\alpha {\nabla }^{2}C,$$

$$\frac{\partial C}{\partial t}+\mathbf{u}\cdot \nabla C=\kappa {\nabla }^2\Phi ,$$

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)$$

$$\nabla \cdot \mathbf{u}$$

$$U=(f-{\rho }_{o}g)\frac{R_I^2}{4{\mu }_o}+(f-{\rho }_{w}g)\frac{R^2-R_I^2}{4{\mu }_w}-({\rho }_o-{\rho }_w)g\frac{R_I^2}{2{\mu }_w}\log \frac{R}{R_I},$$

$${\mathbf{u}}^{\ast }=\frac{\mathbf{u}}{U},{\mathbf{x}}^{\ast }=\frac{\mathbf{x}}{R},$$

$${\rho }^{\ast }=\frac{\rho }{{\rho }_o},{\mu }^{\ast }=$$

$$\rho Re\left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)$$

$$\nabla \cdot \mathbf{u}$$

$$\frac{\partial C}{\partial t}+\mathbf{u}\cdot \nabla C$$

$$\hat{\rho }=\frac{{\rho }_w}{{\rho }_o},\hat{\mu }=\frac{{\mu }_w}{{\mu }_o},$$

$$Re=\frac{{\rho }_{o}UR}{{\mu }_o},Ca=\frac{{\mu }_{o}U}{\sigma },Fr=\frac{U^2}{Rg},$$

$$Cn=\frac{\sqrt{\alpha /\beta }}{R},Pe=\frac{UR}{\kappa \beta },$$

$$Pe\propto C{n}^{-1}.$$

$$\hat{V}=\frac{V_o}{V_w},$$

$$u_z(r)=\{\begin{array}{ll}\left(\frac{1}{{\eta }^2}-\frac{r^2}{{\eta }^2}-2(K_f-1)\text{log}(r)\right)\frac{1}{A}& \text{if }0<r⩽\eta ,\\ 1-\hat{\mu }\frac{r^2}{{\eta }^2}\frac{K_f}{A}& \text{if }\eta <r⩽1,\end{array}$$

$$h=\frac{\hat{\dot{V}}}{\hat{V}},$$

$$A(r,\theta ,z,t)=\sum _{k=-K}^K\sum _{m=-M}^M{\hat{A}}_{k,m}(r,t)e^{(i{k}_{0}kz+im\theta )},$$

$$f=\frac{1}{Re}\frac{4\hat{\mu }}{1-{\eta }^2+{\eta }^2\hat{\mu }}.$$

$$\dot{V}=\frac{\pi }{2}\frac{1-{\eta }^4(1-\hat{\mu })}{1-{\eta }^2(1-\hat{\mu })}.$$

$$\epsilon =\mathcal{O}(Cn\log {\hat{\mu }}^{-1}).$$

$$\mathbf{u}(r,z,\theta ,t)=\mathbf{u}(r,z-ct,\theta ,0),$$

Full text

††journal: International Journal of Multiphase Flow

1 Introduction

Many numerical methods have been developed to deal with the motion of interfaces between immiscible fluids [21]. These methods can be divided into two broad classes: interface-tracking and interface-capturing methods [5]. In interface-tracking (or sharp-interface) approaches, the governing (Navier–Stokes) equations are solved separately for different phases and the coupling occurs through stress and velocity boundary conditions at interfaces. This requires explicit tracking and meshing of the interfaces [28]. Because of the computational cost incurred in re-meshing the interface, volume-of-fluid (VoF) [10] and level-set (LS) methods [27] were developed to ‘capture’ the evolution of the interface. In these methods, a single set of the Navier–Stokes equation for the whole domain is solved together with an advection equation of a scalar field (the volume fraction in VoF and a distance function in LS). The surface-tension force appears then in the Navier–Stokes equations as a volume-force applied in a narrow region near the interface, based on the Young-Laplace formula. These methods yield good approximations of the sharp discontinuities across interfaces, but require the calculation of the surface-normal vector and can lead to difficulties in resolving large interface curvatures. Furthermore, the cost of computing the curvature of the interface rises linearly with the interfacial area. Thus, in problems with large interfacial areas and deformations, such as in non-dilute turbulent dispersions, the additional computational cost of these methods is substantial [19].

1.1 Cahn–Hilliard–Navier–Stokes equations

An alternative interface-capturing approach to computing multiphase flows is the diffusive-interface (or phase-field) method [2]. Here the sharp changes occurring across the interface are smoothed over a finite-thickness mixing layer of width $\epsilon$ separating the phases. Across this mixing layer, the phase-field variable $C$, which describes the composition of the mixture, transitions from the constant value in one phase to that in the other. In our formulation $-0.5\leqslant C\leqslant 0.5$ is chosen, so that $C=-0.5$ and $0.5$ correspond to fluid 1 and 2, respectively. The density and viscosity of the fluids are typically defined as

[TABLE]

for simplicity. The interfacial dynamics is modeled accounting for the free energy of the two-fluid system

[TABLE]

where the material parameters $\alpha$ and $\beta$ determine the contribution of the gradient energy and bulk energy, respectively. The double-well potential $\Psi(C)=(C+0.5)^{2}(C-0.5)^{2}$ models the immiscibility of the two fluids. Equilibrium interface profiles are those that minimize the free energy (2). The interface thickness $\epsilon$ and surface tension $\sigma$ of the system emerge from the competition between the gradient term, which tends to broaden the interface, and the bulk term, which favors phase-separation. In the isothermal case, $\sigma=\sqrt{\alpha\beta/18}$ and the interface thickness $\epsilon$ is typically defined as the region for which $-0.45\leqslant C\leqslant 0.45$. With this choice $\epsilon=4.164\sqrt{\alpha/\beta}$ and the so-defined interface contains about 98.5% of the surface tension stress [12]. The chemical potential of the system

[TABLE]

is the rate of change of the free energy $\mathcal{F}$ with respect to $C$, and so equilibrium interface profiles satisfy $\Phi$=0 [12]. For a planar interface, $C$ has a hyperbolic tangent profile [12]. The temporal evolution of $C$ is governed by the Cahn–Hilliard equation [4]

[TABLE]

where $\kappa$ is the mobility of the chemical potential. The Cahn–Hilliard equation is solved together with the Navier–Stokes equations (CHNS),

[TABLE]

where $\bm{f}$ represents external body forces and $-C\nabla\Phi$ is the motion-causing component of the surface-tension. Its potential component has been absorbed into the generalized pressure $\tilde{p}=p-C\Phi+\beta\Psi-\alpha/2|\nabla C|^{2}$, where $p$ is the true fluid pressure [12]. Note that $C$ also enters the Navier–Stokes equations through the variable density and viscosity of the mixture (1).

The CHNS have been long used to simulate binary fluid flows [2] and more recently to investigate turbulent multiphase flows [25, 1]. Their ability to deal with topological changes and their thermodynamic consistency make them an appealing alternative to other methods. It has been shown that the classical stress balance at fluid-fluid interfaces (employed in sharp-interface methods) is recovered from the CHNS in the limit of vanishingly small interface width $\epsilon\rightarrow 0$, see [2]. Note that this is not the case for other interface-capturing methods, such as volume-of-fluid methods, which may explain their difficulties in producing grid-converged results in some specific problems [16].

1.2 Interface thickness and mobility

For immiscible fluids, $\epsilon$ is a few nanometers wide (e.g. as measured for water and n-alkyl [20]) and the mobility of the chemical potential, though not directly measurable in experiments, was estimated to be $\kappa\sim 10^{-17}\text{m}^{3}\text{kg}^{-1}\text{s}$ [12]. Such values are far beyond the ability of state-of-the-art numerical simulations for problems outside the domain of nanofluidics. In practice, the interface thickness has to be resolved with a few grid points, and this raises the question of how to choose the value of the mobility $\kappa$ consistently. Starting with Jacqmin [12], who argued that the mobility should scale as $\kappa\propto\epsilon^{\delta}$, with $1\leqslant\delta<2$, several scalings have been proposed and tried, see e.g. [15, 32]. Hence, despite the solid theoretical footing and promise of diffuse-interface methods, their predictive power has remained poor because $\kappa$ has been essentially used as a free numerical parameter in previous studies. An improper choice for the mobility parameter (i.e. not consistent with the sharp-interface limit) can lead to a numerical inaccuracy, which hinders the application of the phase-field method to two-phase flows in engineering.

Magaletti et al. [18] resolved this controversy recently. They carried out an asymptotic analysis of the CHNS system for $\epsilon\to 0$ and showed that only the scaling $\kappa\propto\epsilon^{2}$ can correctly match the inner (interface) and outer (bulk) dynamics. This is necessary to recover the correct surface tension force. Physically, the mobility determines the time scale of the diffusion of the chemical potential and cannot be neither too large nor too small (for a given $\epsilon$). If the mobility is too small, then the inner dynamics responds too slowly to advection by the fluid velocity field and the surface-tension force is incorrect. In contrast, if the mobility is too large, then the interface is hardly deformed by the velocity field. In other words, the Cahn–Hilliard equation has to generate the correct $C$ (and $\Phi$) profile at the right rate, i.e. following the change in the outer flow, and this requires a consistent value of the mobility. Magaletti et al. [18] further verified their theoretical derivation with simulations of capillary waves and drop coalescence in two-dimensional Cartesian geometry for matched density and viscosity. In addition, they pointed out that the inner dynamics of the interface, which is much faster than the outer time scale of the fluid velocity field, must be resolved to correctly recover the interfacial physics. Although this can potentially pose a more stringent restriction on the time-step size than numerical stability, it has been scarcely addressed in the literature [18].

1.3 Core-annular flows in pipes

The purpose of this paper is to investigate the numerical convergence of the phase-field model to the sharp interface limit as $\epsilon\to 0$, thereby assessing the robustness of the scaling $\kappa\propto\epsilon^{2}$, in a situation of engineering interest, i.e. with real fluids, practical geometries and flow rates. We focus on the study of core-annular flows (CAF), which are a particular set of two-phase pipe flow regimes where an inner fluid, the “core”, is surrounded by an outer fluid, the “annulus”. These flows are of interest because they allow the lubricated transport of viscous fluids. In the core-annular configuration, the viscous fluid (normally a heavy oil) tends towards the center of the pipe while the less viscous one (usually water) migrates to the high-shear region close to the walls. This structure results in a large reduction in the friction losses when compared with the flow of the single-phase viscous fluid. The interested reader is referred to Joseph et al. [13] for a comprehensive review of CAF, and to Govindarajan and Sahu [8] for a more recent review focusing on miscible and immiscible viscosity-stratified flows. By assuming an infinite or axially periodic pipe, the Navier–Stokes equations with stress and velocity boundary conditions at the interface between the two fluids admit an analytic steady solution termed laminar core-annular flow, for which the interface between the core and annular fluids is parallel to the pipe wall. Figure 1 shows the velocity profiles of laminar CAF for several viscosity ratios and matched densities. When both fluids have the same viscosity, the classic parabolic Poiseuille profile is recovered. As the viscosity ratio increases, the velocity profile becomes steeper in the annular region, whereas in the core the velocity tends to become uniform, i.e. the core exhibits solid-like behavior.

In this paper, we present a comprehensive investigation of the convergence of the CHNS to the sharp-interface limit for CAF with axially periodic boundary conditions. First, the dependence of the convergence rate on the viscosity contrast is investigated for laminar CAF. Second, we examine the convergence of the linear instability of the laminar CAF as a function of the mobility and time-step size. Finally, we compute flow patterns in the fully nonlinear regime and compare them quantitatively to simulations with the VoF method. Overall, our results highlight the ability of the phase-field method to produce accurate results in practical flow situations.

2 Problem specification and methods

We consider upward CAF where the core fluid is oil (subindex o) and the annular fluid is water (subindex w). Here, the flow is driven upward against gravity by imposing a negative average pressure gradient in the axial direction, $\mathrm{d}\overline{p}/\mathrm{d}z<0$. In this case, the centerline velocity of laminar CAF reads

[TABLE]

where $f=-\mathrm{d}\overline{p}/\mathrm{d}z>0$, $R$ is the radius of the pipe, $R_{I}$ the interface position and $\mu_{o}$ ($\mu_{w}$) and $\rho_{o}$ ($\rho_{w}$) are the dynamic viscosity and the density of oil (water). In this paper, all variables and parameters are rendered dimensionless by using the following transformations

[TABLE]

Hereafter, only dimensionless quantities are used and the superscript ∗ is dropped to simplify the notation.

2.1 Governing equations and parameters

The dimensionless CHNS equations for upward CAF read

[TABLE]

where $\bm{T}=\mu(\nabla\bm{u}+\nabla\bm{u}^{T})$ is the viscous stress tensor. The viscosity ratio and the density ratio

[TABLE]

enter the equations through the dimensionless viscosity $\mu$ and density $\rho$ and are solely set by the choice of fluids, whereas the Reynolds, capillary and Froude numbers

[TABLE]

depend also on the fluid velocity. The surface tension effects in the Cahn-Hilliard approach are covered by three dimensionless numbers, namely the capillary number, the Cahn number and the Peclet number. The latter two, which are defined as

[TABLE]

are the dimensionless interface width and the dimensionless inverse of the mobility, respectively, and determine how the sharp-interface limit is approached. In dimensionless form, the relationship $\kappa\propto\epsilon^{2}$ derived by Magaletti et al. [18] turns into

[TABLE]

2.2 Boundary conditions and hold-up ratio

Bai et al. [3] investigated experimentally CAF of very viscous oils lubricated with water. In their experiments, the input volume flow rates of oil ($\dot{V}_{o}$) and water ($\dot{V}_{w}$) were set independently. This setup could be modeled numerically by using inlet boundary conditions in a pipe, but Li and Renardy [17] pointed out that the slow development of the ensuing flow patterns in the streamwise direction makes such simulations infeasible. Instead, they performed simulations with axially periodic boundary conditions, as done in this paper. Note, however, that prescribing both input volume flow rates is not possible in axially periodic pipes, where the two adjustable parameters are the ratio of oil-to-water volumes,

[TABLE]

and the pressure gradient used to drive the flow. The former is fixed by the initial condition, whereas the latter is imposed with a volume force in the right-hand-side of the Navier–Stokes equation ($f$ in eq. 9).

For laminar CAF, it is easy to relate the parameters in numerical simulations of axially periodic pipes to those used in experiments because of the existence of the following analytical solution

[TABLE]

where $\eta$ is the dimensionless position of the interface, related to the volume ratio as $\hat{V}=\eta^{2}/(1-\eta^{2})$, $K_{f}=(Fr\,f-1)/(Fr\,f-\hat{\rho})$ is the driving force ratio and $A=\hat{\mu}K_{f}+1/\eta^{2}-1+2(K_{f}-1)\log(1/\eta)$. By contrast, an analytical solution is not available for the nonlinear regimes, which makes the comparison between numerical simulations and experiments more challenging. The hold-up ratio

[TABLE]

where $\hat{\dot{V}}=\dot{V}_{\text{o}}/\dot{V}_{\text{w}}$, is frequently used to characterize two-phase pipe flows [3]. It is a measure of the relative velocities of both fluids, and in case of a perfectly mixed flow it should have a value of one. In CAF, it is expected that the liquid in contact with the pipe wall is held back and therefore the holdup ratio should be larger than one. To enable comparison with experiment, Li and Renardy [17] extracted the volume ratio from the experimentally measured hold-up ratio, see eq. (18). Then, by using eq. (17), thereby assuming laminar CAF, they computed the pressure gradient yielding the experimentally imposed oil volume flow rate $\dot{V}_{o}$. Note, however, that with this procedure, the water volume flow rate $\dot{V}_{w}$ is fixed indirectly and is only correct in the laminar CAF regime (i.e. in the absence of interface waves).

At the pipe wall, no-slip boundary conditions were applied for the velocity field and no-flux boundary conditions for the phase-field variable and chemical potential. The latter corresponds to 90∘ contact angle at the pipe wall [12]. Note however that this choice is irrelevant because in our simulations the core fluid never touched the wall.

2.3 Numerical method

We used the strategy of [6] to solve the Cahn-Hilliard equation and to treat the variable coefficient of the Laplacian terms in the Navier-Stokes equations. For the time-integration scheme (Crank–Nicolson) and the enforcement of the incompressibility condition (6) with the influence matrix method, we followed the approach of the open-source pipe flow solver openpipeflow [30]. For further details on the time-integration scheme, we refer the reader to [9]. We solved the CHNS in cylindrical coordinates $(r,\theta,z)$ and used 7-point finite-difference stencils to discretize the radial direction, whereas the Fourier–Galerkin spectral method was used in the periodic axial and azimuthal directions. Variables ($\bm{u}$, $p$ and $C$) were written as

[TABLE]

where $k_{0}k$ and $m$ are the wave numbers of the modes in the axial and the azimuthal directions, respectively. The pipe length is $L_{z}=2\pi/k_{0}$, and $\hat{A}_{k,m}$ is the complex Fourier coefficient of mode $(k,m)$.

In the radial direction, the interfacial region ($-0.45<C<0.45$) was discretized with $8$ points, corresponding to $\Delta r\approx 0.5Cn$. This criterion serves to estimate the required resolution for simulations of the CHNS and poses a stark restriction. However, in problems where the interface position is fixed, it can be alleviated by using a non-uniform grid. For the tests presented here we used various grids, all of them presenting a clustering of the radial points close to the wall and near the interface (since its position varied little through the simulation). This rendered a manageable number of grid points in the radial direction, even for very small values of $Cn$, and ensured that spatial discretization errors were negligible. This approach enabled extensive tests of the convergence of the method as a function of the viscosity ratio, $Cn$, $Pe$ and time step, which was one of the main purposes of this work.

The evaluation of the nonlinear terms was performed using the pseudo-spectral technique (with the $3/2$-rule for de-aliasing), which utilizes the fast Fourier transform (FFT) to convert data between physical and spectral spaces. Solving the CH-equation and the full treatment of the viscous term, due to the variable density and viscosity, implies that this method requires roughly 2.5 times more operations for each evaluation of the non-linear terms than a similar single-phase code using the same discretization and time-stepping scheme. We used the MPI-OpenMP hybrid parallelization strategy of [26], which can efficiently utilize $\mathcal{O}(10^{4})$ processors for high-resolution simulations [22].

3 Laminar core-annular flow in the phase-field formulation

We first examined the convergence of the laminar CAF solution to the CHNS equations to the analytic sharp-interface solution (17). For simplicity, we focused on matched density fluids, for which the volume force required to drive the flow in the right-hand-side of the Navier–Stokes equation (9) reads

[TABLE]

In the phase-field formulation there is, to our knowledge, no analytical expression for laminar core-annular flow, so we computed it numerically. Under base flow conditions (radial dependence only, steady) the left-hand-side of the Cahn–Hilliard equation is zero and the basic flow is independent of $Pe$. Moreover, in this case the surface tension force does not cause fluid motion and is thus zero in the potential formulation used here. Hence the Cahn–Hilliard equation decouples from the Navier–Stokes equations and the convergence of the base flow depends only on the interface thickness and the viscosity ratio. The system was solved in one dimension (radial coordinate), in order to obtain the laminar CAF solution without interfacial waves.

After initial transients, our simulations converged toward the velocity profile shown as a red line in figure 2. This agrees excellently with the analytical solution (black line) on the annular side, but is significantly slower in the core region. In particular, the centerline velocity does not reach the value $u_{z}(0)=1$ expected for the analytic solution. A better numerical approximation to the analytical solution in the core is achieved when the fluid motion is driven by imposing a constant (total) volume flux along the pipe (green line). In this case, the driving volume force $f$ is adjusted every time-step so that the total volume flux is identical to that of the analytical solution, namely

[TABLE]

The relative error of the computed laminar solution with respect to the analytical one is shown in figure 3a–b as a function of the radial coordinate for pressure-driven and volume-flux driven flows, respectively. It is interesting to note that despite the difference in the profiles apparent in figure 2, the maximum error is of similar magnitude and occurs near the interface between the two fluids in both cases. Figure 3a confirms that if the pressure gradient of the analytical solution is imposed, the error is concentrated in the core side and decreases toward the wall. In fact, because of axial momentum conservation, with pressure driving the shear stress at the wall should be equal to the analytical one. However, in the Cahn-Hilliard method, there is a slight shift in $C$ in the bulk phases with respect to the expected values, $C=\pm 0.5$. This shift is proportional to $Cn$ [31], and causes an error in the calculation of the viscosity in eq. (1). Hence, even in the pressure-driven case, there is a discrepancy with the analytical profile in the vicinity of the wall. For the volume-flux driven case, the error near the wall is much larger. In particular, imposing the same volume flux as for the analytical solution results in larger shear stress at the wall (see figure  2b) and thus larger pressure gradient to drive the flow.

The maximum error with respect to the analytical solution is shown in figures 3c–d as a function of $Cn$ for several $\hat{\mu}$. Regardless of the driving and the viscosity ratio used, the computed base flow converges linearly to the analytical solution as $Cn$ decreases. This is consistent with previous studies in which the phase-field method was postulated to be first-order accurate with respect to $Cn$ [12, 18]. The effect of the viscosity ratio $\hat{\mu}$ on the error was investigated quantitatively by determining the proportionality constants of the lines in figure 3c–d. These are shown as squares and circles, respectively, in figure 4. Clearly, as $\hat{\mu}$ decreases toward zero, corresponding to very viscous oils in the core, the error increases monotonously. We also computed another case with $\eta=0.7$ (volume ratio $\hat{V}=0.96$) for constant volume-flux driving. The corresponding error is shown as triangles and is substantially smaller than for $\eta=0.9$ ($\hat{V}=4.26$). This is as expected because the shift in $C$ in the bulk phases becomes worse, the more that $\hat{V}$ departs from unity [31].

In order to qualitatively understand the effect of viscosity ratio on the error shown in figure 4, we developed a simple model. We assumed that the equilibrium profile of $C$ approaches a hyperbolic tangent as $Cn\rightarrow 0$ [12]. Then by using this profile in eq. (1) and integrating the one-dimensional steady Navier–Stokes equation, an approximation of the velocity profile was obtained. This model profile also converges to the analytic laminar CAF as $Cn\rightarrow 0$. In particular, it can be shown that at small $\hat{\mu}$ and small $Cn$, the error of the model velocity profile is

[TABLE]

The dashed lines in figure 4 shows that this approximation is indeed very good for low viscosity ratio and constant volume-flux driving. For pressure-driven flow there is more scatter in our data and the agreement is less good (not shown here). Note also that because the assumed $\tanh$-profile for $C$ in the model does not exhibit a shift in the bulk phases, eq. (22) does not capture the dependence of the error on the volume ratio $\hat{V}$.

4 Linear stability analysis of core-annular flow

Hu and Joseph [11] studied the linear stability of density-matched laminar CAF and showed that this is stable only in certain parameter regimes. Furthermore, they did an energy analysis for $Ca=0.5$, $\hat{\mu}=0.1$ and $\eta=\{0.7,0.8\}$ and elucidated the contribution of all the terms in the energy equation, thereby classifying the underlying instabilities according to the dominant physical mechanisms. At low $Re$, capillary instability dominates, whereas as $Re\gtrsim 50$ the interfacial-friction due to the difference in viscosity becomes predominant. For $Re\gtrsim 200$ the production of energy in the bulk of the less viscous (annular) fluid exceeds its dissipation and the instability is predominantly inertial.

We computed the linear-stability of laminar CAF for two specific cases discussed by Hu and Joseph [11] and specified in table 1. The flow was driven by a constant volume flux and simulations were run with a small number of Fourier modes in the axial and azimuthal directions. The initial conditions consisted of the numerically computed laminar CAF, which was disturbed by exciting all Fourier modes with small random noise. After initial transients, the amplitude of the imposed disturbances either increased or decreased exponentially in time and growth rates could be extracted by using exponential fits to the time series of mode amplitudes. This procedure is equivalent to using the power method to compute the leading eigenvalue. The time step size $\Delta t$ was held small (between $2\times 10^{-6}$ and $2\times 10^{-5}$ depending on $Pe$ and $Cn$) to ensure that temporal discretization errors were negligible. The radial resolution was kept large enough (with $7$-$8$ points in the interface region) to ensure that spatial discretization errors were also negligible. This enabled a study of the convergence to the sharp-interface limit as a function of $Cn$ and $Pe$ only.

4.1 CAF1: Instability driven by interfacial-friction

In this test case, there is no surface tension and $C$ enters the Navier–Stokes equations only via the variable viscosity. Following Hu & Joseph [11], two cases were considered to validate our code: axisymmetric ($m=0$) and non-axisymmetric ($m=5$) perturbations to the base flow. Figure 5a shows the fluid velocity fields (arrows) and phase fields $C$ (colormap) of the most unstable perturbation ($m=0$, according to [11]) with axial wavenumber $k=10$. As there is no capillary force and the Reynolds number is small, $Re=37.7$, the instability is driven here by interfacial friction. A similar instability was also investigated in recent CHNS simulations of layered two-phase channel flow without surface tension [24]. Figure 6 shows the growth rate of the most dangerous disturbance, computed with $Cn=0.001$, $Pe=10^{3}$ and $N=288$ radial points, as a function of its axial wavenumber. Our results are in very good agreement with the sharp-interface computations of Hu & Joseph [11] both for axisymmetric ($m=0$) and non-axisymmetric ($m=5$) disturbances. As the axial wavenumber $k$ increases, the difference between our simulation and [11] grows. Larger values of $k$ imply finer structures and therefore more radial finite difference points would be necessary to keep the same level of accuracy.

According to the theory of Magaletti et al. [18], in the absence of surface tension there is no optimal Peclet number and the correct interfacial dynamics should be recovered for sufficiently large $Pe$. Figure 7a shows the relative error between the growth rate obtained in our simulations and that reported by Hu and Joseph [11] for $k=10$. In all cases, the relative error decreases initially as $Pe^{-1}$ when $Pe$ is increased from $O(1)$ to $O(10^{3})$ and then reaches a plateau. If $Pe$ is further increased, the relative error increases monotonically, which differs from the theoretical arguments of [18]. The exception to this is the smallest $Cn=0.001$ considered, where the relative error seems to converge towards a constant value as $Pe$ is increased. This suggests that in this case the prediction of [18] might be valid only very close to the sharp-interface limit. Overall, our results suggest that in practice $Pe=10^{3}$ is sufficient for reproducing the dynamics without surface tension over the range of $Cn$ considered. In fact, figure 8(a) shows that by fixing $Pe=10^{3}$ a faster than linear convergence rate is obtained, as $Cn$ is decreased.

4.2 CAF2: Instability driven by inertia

In this test case, $Re=500$, $\hat{\mu}=1/20$ and the instability is predominantly inertial. Figure 5(b) shows that the velocity field of the disturbance is most intense in the annular region, in agreement with Hu and Joseph [11], who showed that production is negligible with respect to dissipation in the core fluid because of the large viscosity. Figure 7b shows that in the presence of surface tension there is an intermediate range of $Pe$ for which the relative error is minimized. This minimum, which manifests itself clearly only as $Cn$ becomes sufficiently small ($Cn\lesssim 0.005$), is rather broad. It extends over one order of magnitude around the value of $Pe_{\text{opt}}=1/(3Cn)$ proposed by Magaletti et al. [18] as optimal. They proposed the pre-factor $1/3$ from their capillary-wave simulations with matched density and viscosity in a rectangular domain. Our results indicate that $1/3$ also works well in systems with cylindrical geometry, large viscosity ratio, strong inertia and shear. When $Pe$ is decreased linearly, with pre-factors chosen in the near-optimal range, the convergence to the sharp-interface limit is for CAF2 quadratic in $Cn$ (see figure 8(b)). The nearly quadratic convergences agrees with Jacqmin [12], who pointed out that in practice the phase-field method may be faster than linear as predicted theoretically.

5 Water-lubricated pipe flow

Bai et al. [3] investigated experimentally CAF of very viscous oils lubricated with water ($\hat{\mu}=1/601$, $\hat{\rho}=1.1$) in vertical pipes and revealed a wide variety of flow patterns depending on the flow rates of the two phases. In this section, we focus on a case in which the flow is driven upward against gravity with a constant pressure gradient and closely follow the procedure employed by Li and Renardy [17] in their VoF-simulations of this case (see table 2 for the flow parameters). We first obtained the laminar CAF at the corresponding parameters, then determined its linear instability and finally computed the resulting nonlinear flow pattern.

5.1 Linear stability of laminar core-annular flow

Li and Renardy [17] estimated from the experiments of Bai et al. [3] that the dimensionless pressure gradient was $f=1.7455$. In this paper, the value of the pressure gradient was chosen so that the same volume flux as for the laminar CAF in the sharp-interface limit was obtained. This yielded $f=1.7479$ for $Cn=0.001$, obtained with $N=576$ radial points. Figure 9 shows that the computed laminar CAF renders a reasonable approximation of the sharp-interface solution, with 7.5% maximum error with respect to eq. (17), despite the very high viscosity ratio of the fluids used in the experiments.

Following Li and Renardy [17], we destabilized the laminar CAF by disturbing only the first axisymmetric mode in a pipe of length $L=\pi/1.28\approx 2.45$. The computed growth rate as a function of the time-step size is shown in Figure 10a for $Cn=0.001$. The behavior of the temporal error was examined by performing simulations with two different Peclet numbers, $Pe=Pe_{\text{opt}}=333$ and $Pe=0.1Pe_{\text{opt}}=33.3$. Both cases converged towards the same value as $\Delta t\rightarrow 0$ (within a discrepancy of 0.7%), which was calculated by fitting a power law to the finite $\Delta t$ results. However, for $0.1Pe_{\text{opt}}$ and $\Delta t=10^{-4}$, smaller error than for $Pe=Pe_{\text{opt}}$ and $\Delta t=10^{-5}$ was achieved (see figure 10b). Clearly, $0.1Pe_{\text{opt}}$ allows a much larger time-step size for the same accuracy, so we fixed $Pe=0.1Pe_{\text{opt}}$ in the following nonlinear analysis. We note that our result differs by 10% with respect to the sharp-interface limit [17]. Hence it can be concluded that the error in the linear stability analysis stems mainly from the diffuse-interface approximation of the underlying laminar CAF.

The leading eigenvalue of the instability has non-zero imaginary part, and corresponds to a Hopf bifurcation of the laminar CAF. Because of the mean positive advection in the vertical direction, the bifurcating flow pattern is generically expected to be an upward-traveling wave in the axial direction, satisfying

[TABLE]

with $c>0$. The computed linear wave speed $c=0.83$ was found in good agreement with the linear stability analysis of the sharp-interface limit [17], with $c=0.85$.

5.2 Nonlinear bamboo waves

Figure 11a illustrates the temporal evolution of the flow instability starting from a slightly disturbed laminar CAF. At around $t=20$ nonlinear effects begin to kick in and the characteristic bamboo-wave pattern shown in figure 12 fully emerges by $t=30$. Clearly, the interface between the two fluids ceases to be parallel to the pipe axis, and a large number of axial modes is required to accurately resolve the profile of $C$ across the interface. Here, $2K=3456$ axial Fourier modes were used (the number of radial points was not changed, $N=576$). Because of the high computational cost of these two-dimensional nonlinear simulations, we set $\Delta t=2\times 10^{-5}$ and $Pe=33.33$, which gave a good compromise between computing cost and temporal error (see figure 10).

It is worth noting that the agreement between the hold-up ratio in our simulations and those of Li and Renardy [17] is excellent during the transition and the nonlinear saturation phase, whereas the discrepancy in the laminar phase is much larger (see figure 11a). In the laminar flow, the interface and the fluid velocity are parallel to each other. In the Cahn–Hilliard approximation, the interface is smooth and because of the high viscosity ratio, the magnitude of the velocities in the oil and the water phases is lower than in the analytic sharp-interface solution (see figure 9). This results in a significantly lower hold-up ratio. Although eq. (22) shows that this could be ameliorated by further decreasing $Cn$, the computational cost would be enormous. By contrast, the VoF simulation of Li & Renardy [17] can deal with a sharp interface and therefore does not suffer from this problem, rendering better agreement with the analytic solution. In the nonlinear regime, however, the interface is no longer parallel to the fluid velocity and the hold-up ratio of the bamboo pattern becomes less sensitive to the interface thickness. This was confirmed by doubling the Cahn number in our simulations (see the green line in figure 11a), which worsened the agreement in the laminar flow phase, but did not alter the transition process significantly.

The close-up of figure 12d shows that the oil (core) flows upward (from left to right), and part of the water (annular fluid) is trapped in the trough of the wave and is carried upward as well, thereby decreasing the hold-up ratio. The water near the wave crests exhibits swirls and is carried downward. The large gradients in the very thin gap between interface and pipe wall at the wave crest enhance the flow resistance. As a result the transport of oil is strongly reduced, as illustrated by a drop of the centerline velocity by nearly a factor of $3$ with respect to the laminar CAF (see figure 11b).

Finally, we note that although the flow pattern changed litte after $t=30$, the flow continued to evolve until the end of the simulation. Clearly, neither the hold-up ratio nor the centerline velocity saturated at a constant value. Traveling waves satisfying eq. (23) are relative equilibria, for which integral quantities, such as the hold-up ratio, must remain constant in time. This implies that the computed waves are not exact traveling waves, but have suffered further secondary bifurcations already. This is in qualitative agreement with Bai et al. [3], who noted that the bamboo waves observed experimentally were not completely axisymmetric and consisted of the superposition of waves of different wavelengths traveling at different speeds. Studying these phenomena in a realistic setup would require non-axisymmetric (three-dimensional) simulations with much longer domains than used here and is out of the scope of this work.

6 Conclusions and outlook

Phase-field methods have long been used to simulate multiphase flows, mostly to demonstrate their capability of dealing with topological changes, such as droplet break up and coalescence. By contrast, quantitative investigations of their accuracy have remained scarce to date, especially for realistic geometries and fluid properties. This was undertaken here with pseudo-spectral simulations of the Cahn-Hilliard-Navier-Stokes equations for two-phase pipe flow. We validated our code against solutions of the governing equations in the sharp-interface limit (i.e. the Navier–Stokes equations with stress and velocity boundary conditions at the interface between the two fluids). Our tests included an analytic laminar flow solution, linear stability analysis with axisymmetric and non-axisymmetric modes, and a fully nonlinear flow pattern. In all cases, the relationship between interface thickness and mobility derived by Magaletti et al. [18], $Pe\propto Cn^{-1}$, was shown to converge to the sharp-interface limit. Our results suggest that a smaller pre-factor (e.g. $1/30$ instead of their proposed $1/3$) could allow the same accuracy with an order of magnitude larger time step. In general, we expect the optimal pre-factor to be problem dependent.

The setup chosen here exhibits moderate interface areas and deformations, and no topological changes. In this case, interface tracking methods are clearly superior. For example, in laminar CAF the velocity is parallel to the interface and the model error of the CHNS increases as the viscosity ratio $\hat{\mu}$ decreases. In the limit of very low viscosity ratio, corresponding here to very viscous oils in the core, the error is $\varepsilon=\mathcal{O}(Cn\log\hat{\mu}^{-1})$, suggesting that phase-field methods are less suited than other methods, such as VoF, to simulate two-phase flows with very large difference in the viscosity of the two phases. However, our choice of setup was made in order to stringently test phase-field methods against well-known solutions and for flows that can be realized experimentally. Phase-field methods are specially well suited for the treatment of turbulent multiphase flows, where interfaces stretch, deform and even break, largely increasing its surface area [25]. For other interface-capturing methods, such as VoF or level set, which rely on the calculation of the local interface curvature, the additional computing cost cannot be neglected. By contrast, the cost of the simulation with CHNS is independent of the shape and area of the interfaces, which makes CHNS especially attractive for large volume concentrations of the dispersed phase and for turbulent flows  [23]. For the latter, one obvious requirement to faithfully capture the flow physics is that the interface thickness must be smaller than the Kolmogorov scale. Requirements on the accuracy with respect to time-step size remain however to be determined.

Finally, we note that CHNS enable the application of dynamical-system approaches to multiphase flows. In particular, CHNS allow the direct computation of nonlinear traveling waves satisfying eq. (23), and more complex solutions such as relative periodic orbits [29], by solving (9)–(11) directly with the Newton method [14]. This seems to be more difficult with other methods (such as VoF), because the underlying equations must be formulated as a system of smooth partial differential equations, such as (9)–(11). Dynamical-system approaches have helped in elucidating the transition to turbulence in wall-bounded (single-phase) flows [7], and may prove useful to tackle the rich nonlinear dynamics exhibited by two-phase pipe flows [13].

Acknowledgement

We thank Dr. Ashley P. Willis for his help with our understanding of the influence matrix method and the singularity removal in polar coordinates. We acknowledge the Deutsche Forschungsgemeinschaft (Project No. FOR 1182) for financial support and acknowledge the Norddeutsche Verbund für Hoch- und Höchstleistungsrechnen (HLRN) for computing resources.